yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #11397
[Branch ~yade-pkg/yade/git-trunk] Rev 3393: -add solidLine in cell info. partly code for force.
------------------------------------------------------------
revno: 3393
committer: Chao Yuan <chaoyuan2012@xxxxxxxxx>
timestamp: Fri 2014-01-24 19:37:37 +0100
message:
-add solidLine in cell info. partly code for force.
modified:
lib/triangulation/Network.hpp
lib/triangulation/Network.ipp
lib/triangulation/def_types.h
pkg/dem/UnsaturatedEngine.cpp
pkg/dem/UnsaturatedEngine.hpp
--
lp:yade
https://code.launchpad.net/~yade-pkg/yade/git-trunk
Your team Yade developers is subscribed to branch lp:yade.
To unsubscribe from this branch go to https://code.launchpad.net/~yade-pkg/yade/git-trunk/+edit-subscription
=== modified file 'lib/triangulation/Network.hpp'
--- lib/triangulation/Network.hpp 2013-08-22 14:32:01 +0000
+++ lib/triangulation/Network.hpp 2014-01-24 18:37:37 +0000
@@ -83,6 +83,9 @@
Vecteur surface_single_fictious_facet(Vertex_handle fSV1, Vertex_handle SV2, Vertex_handle SV3);
double surface_solid_double_fictious_facet(Vertex_handle SV1, Vertex_handle SV2, Vertex_handle SV3);
double surface_solid_facet(Sphere ST1, Sphere ST2, Sphere ST3);
+
+ void Line_Solid_Pore( Cell_handle cell, int j, bool SLIP_ON_LATERALS, bool reuseFacetData=false);
+ double Line_solid_facet(Sphere ST1, Sphere ST2, Sphere ST3);
int facetF1, facetF2, facetRe1, facetRe2, facetRe3;
int F1, F2, Re1, Re2;
=== modified file 'lib/triangulation/Network.ipp'
--- lib/triangulation/Network.ipp 2013-07-26 13:52:11 +0000
+++ lib/triangulation/Network.ipp 2014-01-24 18:37:37 +0000
@@ -569,6 +569,90 @@
//
// return aire_triangle_spherique;
// }
+
+template<class Tesselation>
+void Network<Tesselation>::Line_Solid_Pore(Cell_handle cell, int j, bool SLIP_ON_LATERALS, bool reuseFacetData)
+{
+ if (!reuseFacetData) facetNFictious=detectFacetFictiousVertices(cell,j);
+ double solidLine = 0; //total of solidLine[j][0], solidLine[j][1], solidLine[j][2].
+ Sphere v [3];
+ Vertex_handle W [3];
+
+ for (int kk=0; kk<3; kk++) {
+ W[kk] = cell->vertex(facetVertices[j][kk]);
+ v[kk] = cell->vertex(facetVertices[j][kk])->point();}
+
+ switch (facetNFictious) {
+ case (0) : {
+ Vertex_handle& SV1 = W[0];
+ Vertex_handle& SV2 = W[1];
+ Vertex_handle& SV3 = W[2];
+
+ cell->info().solidLine[j][0]=Line_solid_facet(SV1->point(), SV2->point(), SV3->point());
+ cell->info().solidLine[j][1]=Line_solid_facet(SV2->point(), SV3->point(), SV1->point());
+ cell->info().solidLine[j][2]=Line_solid_facet(SV3->point(), SV1->point(), SV2->point());
+ }; break;
+ case (1) : {
+ Vertex_handle SV1 = cell->vertex(facetVertices[j][facetRe1]);
+ Vertex_handle SV2 = cell->vertex(facetVertices[j][facetRe2]);
+ Vertex_handle SV3 = cell->vertex(facetVertices[j][facetF1]);
+
+ cell->info().solidLine[j][facetRe1]=Line_solid_facet(SV2->point(), SV3->point(), SV1->point());
+ cell->info().solidLine[j][facetRe2]=Line_solid_facet(SV3->point(), SV1->point(), SV2->point());
+
+ Boundary &bi = boundary(SV3->info().id());
+
+ Vecteur AB= SV1->point() - SV2->point();
+ AB[bi.coordinate] = 0;
+
+ if (bi.flowCondition && ! SLIP_ON_LATERALS) {
+ cell->info().solidLine[j][facetF1]=sqrt(AB.squared_length());
+ } else cell->info().solidLine[j][facetF1]=0;
+ }; break;
+ case (2) : {
+ Vertex_handle SV1 = cell->vertex(facetVertices[j][facetF1]);
+ Vertex_handle SV2 = cell->vertex(facetVertices[j][facetF2]);
+ Vertex_handle SV3 = cell->vertex(facetVertices[j][facetRe1]);
+
+ cell->info().solidLine[j][facetRe1]=0.5*M_PI*sqrt(SV3->point().weight());
+
+ Boundary &bi1 = boundary(SV1->info().id());
+ Boundary &bi2 = boundary(SV2->info().id());
+
+ double d13 = (SV1->point())[bi1.coordinate] - (SV3->point())[bi1.coordinate];
+ double d23 = (SV2->point())[bi2.coordinate] - (SV3->point())[bi2.coordinate];
+ if (bi1.flowCondition && ! SLIP_ON_LATERALS) {
+ cell->info().solidLine[j][facetF1]= abs(d23);
+ } else cell->info().solidLine[j][facetF1]=0;
+
+ if (bi2.flowCondition && ! SLIP_ON_LATERALS) {
+ cell->info().solidLine[j][facetF2]= abs(d13);
+ } else cell->info().solidLine[j][facetF2]=0;
+ }; break;
+ }
+
+ double Lsolid = cell->info().solidLine[j][0] + cell->info().solidLine[j][1] + cell->info().solidLine[j][2];
+ if (Lsolid)
+ cell->info().solidLine[j][3]=1.0/Lsolid;
+ else cell->info().solidLine[j][3]=0;
+}
+
+template<class Tesselation>
+double Network<Tesselation>::Line_solid_facet(Sphere ST1, Sphere ST2, Sphere ST3)
+{
+ double Line;
+ double squared_radius = ST1.weight();
+ Vecteur v12 = ST2.point() - ST1.point();
+ Vecteur v13 = ST3.point() - ST1.point();
+
+ double norme12 = v12.squared_length();
+ double norme13 = v13.squared_length();
+
+ double cosA = v12*v13 / (sqrt(norme13 * norme12));
+ Line = acos(cosA) * sqrt(squared_radius);
+ return Line;
+}
+
} //namespace CGT
#endif //FLOW_ENGINE
=== modified file 'lib/triangulation/def_types.h'
--- lib/triangulation/def_types.h 2014-01-13 16:54:44 +0000
+++ lib/triangulation/def_types.h 2014-01-24 18:37:37 +0000
@@ -224,12 +224,13 @@
bool isWaterReservoir;
bool isAirReservoir;
Real capillaryCellVolume;//for calculating saturation
- std::vector<double> poreRadius;
- //pore throat radius for drainage
+ std::vector<double> poreRadius;//pore throat radius for drainage
+ double solidLine [4][4];//the length of intersecting line between sphere and facet. [i][j] is for sphere facet "i" and sphere facetVertices[i][j]. Last component for 1/sum_Lines in the facet.
UnsatCellInfo (void)
{
poreRadius.resize(4, 0);
isWaterReservoir = false; isAirReservoir = false; capillaryCellVolume = 0;
+ for (int k=0; k<4;k++) for (int l=0; l<3;l++) solidLine[k][l]=0;
}
};
=== modified file 'pkg/dem/UnsaturatedEngine.cpp'
--- pkg/dem/UnsaturatedEngine.cpp 2014-01-14 14:04:00 +0000
+++ pkg/dem/UnsaturatedEngine.cpp 2014-01-24 18:37:37 +0000
@@ -350,7 +350,7 @@
}
}
-template<class Cellhandle>//waiting check
+template<class Cellhandle>
double UnsaturatedEngine::bisection(Cellhandle cell, int j, double a, double b)
{
double m = 0.5*(a+b);
@@ -881,7 +881,7 @@
Real UnsaturatedEngine::getSaturation (Solver& flow )
{
updateAirReservoir(flow);
-// updateWaterReservoir(flow);
+ updateWaterReservoir(flow);
RTriangulation& tri = flow->T[flow->currentTes].Triangulation();
Real capillary_volume = 0.0; //total capillary volume
Real air_volume = 0.0; //air volume
@@ -1346,6 +1346,103 @@
file.close();
}*/
//----------end tempt function for Vahid Joekar-Niasar's data (clear later)---------------------
+
+template <class Solver>
+void UnsaturatedEngine::computeFacetForcesWithCache(Solver& flow, bool onlyCache)
+{
+ RTriangulation& Tri = flow->T[currentTes].Triangulation();
+ Finite_cells_iterator cell_end = Tri.finite_cells_end();
+ CGT::Vecteur nullVect(0,0,0);
+ //reset forces
+ if (!onlyCache) for (Finite_vertices_iterator v = Tri.finite_vertices_begin(); v != Tri.finite_vertices_end(); ++v) v->info().forces=nullVect;
+
+ #ifdef parallel_forces
+ if (noCache) {
+ perVertexUnitForce.clear(); perVertexPressure.clear();
+// vector<const Vecteur*> exf; exf.reserve(20);
+// vector<const Real*> exp; exp.reserve(20);
+ perVertexUnitForce.resize(Tri.number_of_vertices());
+ perVertexPressure.resize(Tri.number_of_vertices());}
+ #endif
+
+ Cell_handle neighbour_cell;
+ Vertex_handle mirror_vertex;
+ CGT::Vecteur tempVect;
+ //FIXME : Ema, be carefull with this (noCache), it needs to be turned true after retriangulation
+ if (noCache) {for (VCell_iterator cell_it=flow->T[currentTes].cellHandles.begin(); cell_it!=flow->T[currentTes].cellHandles.end(); cell_it++){
+// if (noCache) for (Finite_cells_iterator cell = Tri.finite_cells_begin(); cell != cell_end; cell++) {
+ Cell_handle& cell = *cell_it;
+ //reset cache
+ for (int k=0;k<4;k++) cell->info().unitForceVectors[k]=nullVect;
+
+ for (int j=0; j<4; j++) if (!Tri.is_infinite(cell->neighbor(j))) {
+ neighbour_cell = cell->neighbor(j);
+ const CGT::Vecteur& Surfk = cell->info().facetSurfaces[j];
+ //FIXME : later compute that fluidSurf only once in hydraulicRadius, for now keep full surface not modified in cell->info for comparison with other forces schemes
+ //The ratio void surface / facet surface
+ Real area = sqrt(Surfk.squared_length()); if (area<=0) cerr <<"AREA <= 0!!"<<endl;
+ CGT::Vecteur facetNormal = Surfk/area;
+ const std::vector<CGT::Vecteur>& crossSections = cell->info().facetSphereCrossSections;
+ CGT::Vecteur fluidSurfk = cell->info().facetSurfaces[j]*cell->info().facetFluidSurfacesRatio[j];
+ /// handle fictious vertex since we can get the projected surface easily here
+ if (cell->vertex(j)->info().isFictious) {
+ Real projSurf=abs(Surfk[flow->boundary(cell->vertex(j)->info().id()).coordinate]);
+ tempVect=-projSurf*flow->boundary(cell->vertex(j)->info().id()).normal;
+ cell->vertex(j)->info().forces = cell->vertex(j)->info().forces+tempVect*cell->info().p();
+ //define the cached value for later use with cache*p
+ cell->info().unitForceVectors[j]=cell->info().unitForceVectors[j]+ tempVect;
+ }
+ /// Apply weighted forces f_k=sqRad_k/sumSqRad*f
+ CGT::Vecteur Facet_Unit_Force = -fluidSurfk*cell->info().solidLine[j][3];
+ CGT::Vecteur Facet_Force = cell->info().p()*cell->info().poreRadius[j]*Facet_Unit_Force;
+
+
+ for (int y=0; y<3;y++) {
+ cell->vertex(facetVertices[j][y])->info().forces = cell->vertex(facetVertices[j][y])->info().forces + Facet_Force*cell->info().solidLine[j][y];
+ //add to cached value
+ cell->info().unitForceVectors[facetVertices[j][y]]=cell->info().unitForceVectors[facetVertices[j][y]]+Facet_Unit_Force*cell->info().solidLine[j][y];
+ //uncomment to get total force / comment to get only pore tension forces
+ if (!cell->vertex(facetVertices[j][y])->info().isFictious) {
+ cell->vertex(facetVertices[j][y])->info().forces = cell->vertex(facetVertices[j][y])->info().forces -facetNormal*cell->info().p()*crossSections[j][y];
+ //add to cached value
+ cell->info().unitForceVectors[facetVertices[j][y]]=cell->info().unitForceVectors[facetVertices[j][y]]-facetNormal*crossSections[j][y];
+ }
+ }
+ #ifdef parallel_forces
+ perVertexUnitForce[cell->vertex(j)->info().id()].push_back(&(cell->info().unitForceVectors[j]));
+ perVertexPressure[cell->vertex(j)->info().id()].push_back(&(cell->info().p()));
+ #endif
+ }
+ }
+ noCache=false;//cache should always be defined after execution of this function
+ if (onlyCache) return;
+ } else {//use cached values
+ #ifndef parallel_forces
+ for (Finite_cells_iterator cell = Tri.finite_cells_begin(); cell != cell_end; cell++) {
+ for (int yy=0;yy<4;yy++) cell->vertex(yy)->info().forces = cell->vertex(yy)->info().forces + cell->info().unitForceVectors[yy]*cell->info().p();}
+
+ #else
+ #pragma omp parallel for num_threads(ompThreads)
+ for (int vn=0; vn<= flow->T[currentTes].max_id; vn++) {
+ Vertex_handle& v = flow->T[currentTes].vertexHandles[vn];
+// for (Finite_vertices_iterator v = Tri.finite_vertices_begin(); v != Tri.finite_vertices_end(); ++v){
+ const int& id = v->info().id();
+ CGT::Vecteur tf (0,0,0);
+ int k=0;
+ for (vector<const Real*>::iterator c = perVertexPressure[id].begin(); c != perVertexPressure[id].end(); c++)
+ tf = tf + (*(perVertexUnitForce[id][k++]))*(**c);
+ v->info().forces = tf;
+ }
+ #endif
+ }
+ if (flow->DEBUG_OUT) {
+ CGT::Vecteur TotalForce = nullVect;
+ for (Finite_vertices_iterator v = Tri.finite_vertices_begin(); v != Tri.finite_vertices_end(); ++v) {
+ if (!v->info().isFictious) TotalForce = TotalForce + v->info().forces;
+ else if (flow->boundary(v->info().id()).flowCondition==1) TotalForce = TotalForce + v->info().forces; }
+ cout << "TotalForce = "<< TotalForce << endl;}
+}
+
YADE_PLUGIN ( ( UnsaturatedEngine ) );
#endif //FLOW_ENGINE
=== modified file 'pkg/dem/UnsaturatedEngine.hpp'
--- pkg/dem/UnsaturatedEngine.hpp 2014-01-13 16:54:44 +0000
+++ pkg/dem/UnsaturatedEngine.hpp 2014-01-24 18:37:37 +0000
@@ -29,8 +29,10 @@
typedef RTriangulation::Vertex_handle Vertex_handle;
typedef RTriangulation::Point CGALSphere;
typedef CGALSphere::Point Point;
-
+ typedef std::vector<Cell_handle> Vector_Cell;
+ typedef typename Vector_Cell::iterator VCell_iterator;
+
protected:
shared_ptr<FlowSolver> solver;
shared_ptr<FlowSolver> backgroundSolver;
@@ -118,6 +120,11 @@
double getPorePressure(Vector3r pos){return solver->getPorePressure(pos[0], pos[1], pos[2]);}
TPL int getCell(double posX, double posY, double posZ, Solver& flow){return flow->getCell(posX, posY, posZ);}
+ bool noCache;//flag for checking if cached values cell->unitForceVectors have been defined
+ TPL void computeFacetForcesWithCache(Solver& flow, bool onlyCache=false);
+ vector< vector<const CGT::Vecteur*> > perVertexUnitForce;
+ vector< vector<const Real*> > perVertexPressure;
+
void emulateAction(){
scene = Omega::instance().getScene().get();
action();}
@@ -141,12 +148,13 @@
void _savePoreBodyInfo(){savePoreBodyInfo(solver);}
void _savePoreThroatInfo(){savePoreThroatInfo(solver);}
void _debugTemp(){debugTemp(solver);}
-
virtual ~UnsaturatedEngine();
virtual void action();
YADE_CLASS_BASE_DOC_ATTRS_DEPREC_INIT_CTOR_PY(UnsaturatedEngine,PartialEngine,"Preliminary version engine of a model for unsaturated soils",
+// ((bool,drainageActivated,true,,"Activate drainage"))//use later
+// ((bool,imbibitionActivated,true,,"Activate imbibition"))//use later
((bool,first,true,,"Controls the initialization/update phases"))
((bool, Debug, false,,"Activate debug messages"))
((double, wall_thickness,0.001,,"Walls thickness"))
