------------------------------------------------------------ revno: 2764 committer: Bruno Chareyre <bruno.chare...@hmg.inpg.fr> branch nick: yade timestamp: Fri 2011-02-25 20:07:36 +0100 message: - per-point imposed pressure mechanism. modified: pkg/dem/FlowEngine.cpp pkg/dem/FlowEngine.hpp
-- lp:yade https://code.launchpad.net/~yade-dev/yade/trunk Your team Yade developers is subscribed to branch lp:yade. To unsubscribe from this branch go to https://code.launchpad.net/~yade-dev/yade/trunk/+edit-subscription
=== modified file 'pkg/dem/FlowEngine.cpp' --- pkg/dem/FlowEngine.cpp 2011-02-17 17:43:37 +0000 +++ pkg/dem/FlowEngine.cpp 2011-02-25 19:07:36 +0000 @@ -19,6 +19,8 @@ #include <sys/stat.h> #include <sys/types.h> +int dd=0; + CREATE_LOGGER (FlowEngine); FlowEngine::~FlowEngine() @@ -54,7 +56,7 @@ currentStrain = triaxialCompressionEngine->strain[1]; timingDeltas->start(); - + if (first) Build_Triangulation(P_zero); timingDeltas->checkpoint("Triangulating"); @@ -110,7 +112,7 @@ std::ofstream settle("settle.txt", std::ios::app); settle << j << " " << time << " " << currentStrain << endl; } - + if (slice_pressures){ char slifile [30]; mkdir("./Slices", S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH); @@ -120,23 +122,20 @@ char *f = "slifile"; flow->SliceField(f); } - - if (save_vtk) {flow->save_vtk_file();} + if (save_vtk &&dd++>50) {flow->save_vtk_file(); dd=0;} } // if ( scene->iter % PermuteInterval == 0 ) // { Update_Triangulation = true; } if (Update_Triangulation && !first) { Build_Triangulation(); - Update_Triangulation = false; - } - + Update_Triangulation = false;} + if (velocity_profile) /*flow->FluidVelocityProfile();*/flow->Average_Fluid_Velocity(); if (first && liquefaction){ wall_up_y = flow->y_max; wall_down_y = flow->y_min;} if (liquefaction) flow->Measure_Pore_Pressure(wall_up_y, wall_down_y); - first=false; // } } @@ -159,6 +158,24 @@ } +void FlowEngine::imposePressure(Vector3r pos, Real p) +{ + flow->imposedP.push_back( pair<CGT::Point,Real>(CGT::Point(pos[0],pos[1],pos[2]),p) ); +} + +void FlowEngine::clearImposedPressure() { flow->imposedP.clear();} + +Real FlowEngine::getFlux(int cond) { + if (cond>=flow->imposedP.size()) LOG_ERROR("Getting flux with cond higher than imposedP size."); + CGT::RTriangulation& Tri = flow->T[flow->currentTes].Triangulation(); + double flux=0; + CGT::Cell_handle cell= Tri.locate(flow->imposedP[cond].first); + for (int ngb=0;ngb<4;ngb++) { + if (!cell->neighbor(ngb)->info().Pcondition) flux+= cell->info().k_norm()[ngb]*(cell->info().p()-cell->neighbor(ngb)->info().p());} + return flux; +} + + void FlowEngine::Build_Triangulation () { Build_Triangulation (0.f); @@ -188,14 +205,14 @@ Triangulate ( ); if (Debug) cout << endl << "Tesselating------" << endl << endl; flow->T[flow->currentTes].Compute(); - + flow->Define_fictious_cells(); flow->DisplayStatistics (); flow->meanK_LIMIT = meanK_correction; flow->meanK_STAT = meanK_opt; flow->Compute_Permeability ( ); - + porosity = flow->V_porale_porosity/flow->V_totale_porosity; if (first) @@ -251,28 +268,28 @@ flow->z_max = max ( flow->z_max, z+rad); } } - + id_offset = flow->T[flow->currentTes].max_id+1; - + flow->id_offset = id_offset; - + flow->SectionArea = ( flow->x_max - flow->x_min ) * ( flow->z_max-flow->z_min ); flow->Vtotale = (flow->x_max-flow->x_min) * (flow->y_max-flow->y_min) * (flow->z_max-flow->z_min); - + flow->y_min_id=triaxialCompressionEngine->wall_bottom_id; flow->y_max_id=triaxialCompressionEngine->wall_top_id; flow->x_max_id=triaxialCompressionEngine->wall_right_id; flow->x_min_id=triaxialCompressionEngine->wall_left_id; flow->z_min_id=triaxialCompressionEngine->wall_back_id; flow->z_max_id=triaxialCompressionEngine->wall_front_id; - + flow->boundary ( flow->y_min_id ).useMaxMin = BOTTOM_Boundary_MaxMin; flow->boundary ( flow->y_max_id ).useMaxMin = TOP_Boundary_MaxMin; flow->boundary ( flow->x_max_id ).useMaxMin = RIGHT_Boundary_MaxMin; flow->boundary ( flow->x_min_id ).useMaxMin = LEFT_Boundary_MaxMin; flow->boundary ( flow->z_max_id ).useMaxMin = FRONT_Boundary_MaxMin; - flow->boundary ( flow->z_min_id ).useMaxMin = BACK_Boundary_MaxMin; - + flow->boundary ( flow->z_min_id ).useMaxMin = BACK_Boundary_MaxMin; + //FIXME: Id's order in boundsIds is done according to the enumerotation of boundaries from TXStressController.hpp, line 31. DON'T CHANGE IT! flow->boundsIds[0]= &flow->y_min_id; flow->boundsIds[1]= &flow->y_max_id; @@ -280,12 +297,12 @@ flow->boundsIds[3]= &flow->x_max_id; flow->boundsIds[4]= &flow->z_max_id; flow->boundsIds[5]= &flow->z_min_id; - + wall_thickness = triaxialCompressionEngine->thickness; - + flow->Corner_min = CGT::Point(flow->x_min, flow->y_min, flow->z_min); flow->Corner_max = CGT::Point(flow->x_max, flow->y_max, flow->z_max); - + if (Debug) { cout << "Section area = " << flow->SectionArea << endl; @@ -298,10 +315,10 @@ cout << "z_min = " << flow->z_min << endl; cout << "z_max = " << flow->z_max << endl; cout << endl << "Adding Boundary------" << endl;} - + double center[3]; - - for (int i=0; i<6; i++) + + for (int i=0; i<6; i++) { CGT::Vecteur Normal (triaxialCompressionEngine->normal[i].x(), triaxialCompressionEngine->normal[i].y(), triaxialCompressionEngine->normal[i].z()); if (flow->boundary(*flow->boundsIds[i]).useMaxMin) flow->AddBoundingPlane (true, Normal, *flow->boundsIds[i]); @@ -309,7 +326,7 @@ const shared_ptr<Body>& wll = Body::byId ( *flow->boundsIds[i] , scene ); for ( int h=0;h<3;h++ ){center[h] = wll->state->pos[h];} flow->AddBoundingPlane (center, wall_thickness, Normal,*flow->boundsIds[i]);}} - + // flow->AddBoundingPlanes(true); } @@ -451,7 +468,7 @@ Real A[3]={0, 0, 0}, AS[3]={0, 0, 0}, AT[3]={0, 0, 0}; Real B[3]={0, 0, 0}, BS[3]={0, 0, 0}, BT[3]={0, 0, 0}; Real C[3]={0, 0, 0}, CS[3]={0, 0, 0}, CT[3]={0, 0, 0}; - + int b[2]; Real Wall_coordinate[2]; int j=0; === modified file 'pkg/dem/FlowEngine.hpp' --- pkg/dem/FlowEngine.hpp 2011-02-17 17:43:37 +0000 +++ pkg/dem/FlowEngine.hpp 2011-02-25 19:07:36 +0000 @@ -48,12 +48,15 @@ Real Volume_cell (CGT::Cell_handle cell); void Oedometer_Boundary_Conditions(); void BoundaryConditions(); + void imposePressure(Vector3r pos, Real p); + void clearImposedPressure(); + Real getFlux(int cond); virtual ~FlowEngine(); virtual void action(); - YADE_CLASS_BASE_DOC_ATTRS_CTOR(FlowEngine,PartialEngine,"An engine to solve flow problem in saturated granular media", + YADE_CLASS_BASE_DOC_ATTRS_INIT_CTOR_PY(FlowEngine,PartialEngine,"An engine to solve flow problem in saturated granular media", ((bool,isActivated,true,,"Activates Flow Engine")) ((bool,first,true,,"Controls the initialization/update phases")) ((bool,save_vtk,false,,"Enable/disable vtk files creation for visualization")) @@ -111,7 +114,13 @@ ((bool, LEFT_Boundary_MaxMin, 1,,"If true bounding sphere is added as function fo max/min sphere coord, if false as function of yade wall position")) ((bool, FRONT_Boundary_MaxMin, 1,,"If true bounding sphere is added as function fo max/min sphere coord, if false as function of yade wall position")) ((bool, BACK_Boundary_MaxMin, 1,,"If true bounding sphere is added as function fo max/min sphere coord, if false as function of yade wall position")) - ,timingDeltas=shared_ptr<TimingDeltas>(new TimingDeltas)); + ,, + timingDeltas=shared_ptr<TimingDeltas>(new TimingDeltas); + , + .def("imposePressure",&FlowEngine::imposePressure,(python::arg("pos"),python::arg("p")),"Impose pressure in cell of location 'pos'.") + .def("clearImposedPressure",&FlowEngine::clearImposedPressure,"Clear the list of points with pressure imposed.") + .def("getFlux",&FlowEngine::getFlux,(python::arg("cond")),"Get influx in cell associated to an imposed P (indexed using 'cond').") + ) DECLARE_LOGGER; };
_______________________________________________ Mailing list: https://launchpad.net/~yade-dev Post to : yade-dev@lists.launchpad.net Unsubscribe : https://launchpad.net/~yade-dev More help : https://help.launchpad.net/ListHelp