yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #10590
[Branch ~yade-pkg/yade/git-trunk] Rev 3844: some renaming and code cleaning in PFV code (part2)
------------------------------------------------------------
revno: 3844
committer: Bruno Chareyre <bruno.chareyre@xxxxxxxxxxxxxxx>
timestamp: Fri 2014-03-21 19:45:24 +0100
message:
some renaming and code cleaning in PFV code (part2)
removed:
lib/triangulation/def_types.h
modified:
lib/triangulation/FlowBoundingSphere.cpp
lib/triangulation/FlowBoundingSphere.hpp
lib/triangulation/FlowBoundingSphere.ipp
lib/triangulation/FlowBoundingSphereLinSolv.hpp
lib/triangulation/FlowBoundingSphereLinSolv.ipp
lib/triangulation/KinematicLocalisationAnalyser.cpp
lib/triangulation/KinematicLocalisationAnalyser.hpp
lib/triangulation/Network.hpp
lib/triangulation/Network.ipp
lib/triangulation/PeriodicFlow.cpp
lib/triangulation/RegularTriangulation.h
lib/triangulation/Tenseur3.cpp
lib/triangulation/Tenseur3.h
lib/triangulation/Tesselation.h
lib/triangulation/TriaxialState.cpp
lib/triangulation/TriaxialState.h
pkg/dem/FlowEngine.cpp
pkg/dem/FlowEngine.hpp
pkg/dem/MicroMacroAnalyser.cpp
--
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/FlowBoundingSphere.cpp'
--- lib/triangulation/FlowBoundingSphere.cpp 2013-06-26 18:15:55 +0000
+++ lib/triangulation/FlowBoundingSphere.cpp 2014-03-21 18:45:24 +0000
@@ -9,9 +9,9 @@
#include "yade/lib/triangulation/FlowBoundingSphere.hpp"
namespace CGT {
-Vecteur PeriodicCellInfo::gradP;
-Vecteur PeriodicCellInfo::hSize[3];
-Vecteur PeriodicCellInfo::deltaP;
+CVector PeriodicCellInfo::gradP;
+CVector PeriodicCellInfo::hSize[3];
+CVector PeriodicCellInfo::deltaP;
}
=== modified file 'lib/triangulation/FlowBoundingSphere.hpp'
--- lib/triangulation/FlowBoundingSphere.hpp 2014-03-21 18:45:24 +0000
+++ lib/triangulation/FlowBoundingSphere.hpp 2014-03-21 18:45:24 +0000
@@ -37,8 +37,8 @@
FlowBoundingSphere();
bool slipOnLaterals;
- double TOLERANCE;
- double RELAX;
+ double tolerance;
+ double relax;
double ks; //Hydraulic Conductivity
bool clampKValues, meanKStat, distance_correction;
bool OUTPUT_BOUDARIES_RADII;
@@ -67,7 +67,7 @@
#define parallel_forces
#ifdef parallel_forces
int ompThreads;
- vector< vector<const Vecteur*> > perVertexUnitForce;
+ vector< vector<const CVector*> > perVertexUnitForce;
vector< vector<const Real*> > perVertexPressure;
#endif
vector <double> edgeSurfaces;
@@ -100,18 +100,12 @@
bool tessBasedForce; //allow the force computation method to be chosen from FlowEngine
Real minPermLength; //min branch length for Poiseuille
- double VISCOSITY;
+ double viscosity;
double fluidBulkModulus;
- Tesselation& computeAction ( );
- Tesselation& computeAction ( int argc, char *argv[ ], char *envp[ ] );
- Tesselation& loadPositions(int argc, char *argv[ ], char *envp[ ]);
void displayStatistics();
void initializePressures ( double pZero );
bool reApplyBoundaryConditions ();
- /// Define forces using the same averaging volumes as for permeability
- void computeTetrahedralForces();
- /// Define forces spliting drag and buoyancy terms
void computeFacetForcesWithCache(bool onlyCache=false);
void saveVtk (const char* folder );
#ifdef XVIEW
@@ -123,7 +117,7 @@
double computeHydraulicRadius (CellHandle cell, int j );
Real checkSphereFacetOverlap(const Sphere& v0, const Sphere& v1, const Sphere& v2);
- double dotProduct ( Vecteur x, Vecteur y );
+ double dotProduct ( CVector x, CVector y );
double computeEffectiveRadius(CellHandle cell, int j);
double computeEquivalentRadius(CellHandle cell, int j);
//return the list of constriction values
@@ -162,7 +156,7 @@
vector<Real> averageFluidVelocityOnSphere(unsigned int Id_sph);
//Solver?
- int useSolver;//(0 : GaussSeidel, 1 : TAUCS, 2 : PARDISO, 3:CHOLMOD)
+ int useSolver;//(0 : GaussSeidel, 1:CHOLMOD)
};
} //namespace CGT
=== modified file 'lib/triangulation/FlowBoundingSphere.ipp'
--- lib/triangulation/FlowBoundingSphere.ipp 2014-03-21 18:45:24 +0000
+++ lib/triangulation/FlowBoundingSphere.ipp 2014-03-21 18:45:24 +0000
@@ -7,11 +7,6 @@
* GNU General Public License v2 or later. See file LICENSE for details. *
*************************************************************************/
#ifdef FLOW_ENGINE
-// #include "def_types.h"
-// #include "def_flow_types.h"
-// #include "CGAL/constructions/constructions_on_weighted_points_cartesian_3.h"
-// #include <CGAL/Width_3.h>
-
// #define XVIEW
#include "FlowBoundingSphere.hpp"//include after #define XVIEW
@@ -23,7 +18,6 @@
#include <assert.h>
#include <sys/stat.h>
#include <sys/types.h>
-// #include "Network.hpp"
#include <omp.h>
@@ -31,10 +25,6 @@
// #include "Vue3D.h" //FIXME implicit dependencies will look for this class (out of tree) even ifndef XVIEW
#endif
-#define FAST
-#define TESS_BASED_FORCES
-#define FACET_BASED_FORCES 1
-
#ifdef YADE_OPENMP
// #define GS_OPEN_MP //It should never be defined if Yade is not using openmp
#endif
@@ -70,14 +60,14 @@
nOfSpheres = 0;
sectionArea = 0, Height=0, vTotal=0;
vtkInfiniteVertices=0, vtkInfiniteCells=0;
- VISCOSITY = 1;
+ viscosity = 1;
fluidBulkModulus = 0;
tessBasedForce = true;
for (int i=0;i<6;i++) boundsIds[i] = 0;
minPermLength=-1;
slipOnLaterals = false;//no-slip/symmetry conditions on lateral boundaries
- TOLERANCE = 1e-07;
- RELAX = 1.9;
+ tolerance = 1e-07;
+ relax = 1.9;
ks=0;
distance_correction = true;
clampKValues = true;
@@ -101,189 +91,6 @@
template <class Tesselation>
void FlowBoundingSphere<Tesselation>::resetNetwork() {T[0].Clear();noCache=true;}
-template <class Tesselation>
-Tesselation& FlowBoundingSphere<Tesselation>::computeAction()
-{
- return computeAction(0,NULL,NULL);
-}
-template <class Tesselation>
-Tesselation& FlowBoundingSphere<Tesselation>::computeAction(int argc, char *argv[ ], char *envp[ ])
-{
- double factor = 1.001;
- VectorR X, Y, Z, R;
- Real_timer clock;
- clock.start();
- clock.top("start");
- Tesselation& Tes = T[0];
- RTriangulation& Tri = Tes.Triangulation();
-
- /** READING SPHERES POSITIONS FROM A TEXT FILE CONTAINING COORDINATES **/
- double x, y, z, r;
- ifstream loadFile(argc==1 ? "cube" : argv[1]); // cree l'objet loadFile de la classe ifstream qui va permettre de lire importFilename
- while (!loadFile.eof()) {
- loadFile >> x >> y >> z >> r;
- X.push_back(x);
- Y.push_back(y);
- Z.push_back(z);
- R.push_back(factor*r);
- nOfSpheres++;
- Rmoy += r;
- xMin = min(xMin,x-r);
- xMax = max(xMax,x+r);
- yMin = min(yMin,y-r);
- yMax = max(yMax,y+r);
- zMin = min(zMin,z-r);
- zMax = max(zMax,z+r);
- }
- Rmoy /= nOfSpheres;
- minPermLength = Rmoy*minLength;
- if (debugOut) cout << "Rmoy = " << Rmoy << endl;
- if (debugOut) cout << "xMin = " << xMin << " xMax = " << xMax << " yMin = " << yMin << " yMax = " << yMax << " zMin = " << zMin << " zMax = " << zMax << endl;
-
- VertexHandle Vh;
- CellHandle neighbourCell, cell, location;
-
- int V = X.size();
- if (debugOut) cout << "V =" << V << "nOfSpheres = " << nOfSpheres << endl;
- if (debugOut) cout << Tes.Max_id() << endl;
- clock.top("loading spheres");
-
-
- vector<Sphere> vs; RTriangulation testT;
- for (int i=0; i<V; i++) {
- vs.push_back(Sphere(Point(X[i],Y[i],Z[i]), R[i]));
- }
- clock.top("make a spheres vector");
- testT.insert(vs.begin(),vs.end());
- clock.top("test speed");
-
- addBoundingPlanes();
- for (int i=0; i<V; i++) {
- int id = Tes.Max_id() +1;
- Vh = Tes.insert(X[i],Y[i],Z[i],R[i],id); /** EMPILEMENT QUELCONQUE **/
-#ifdef XVIEW
- Vue1.SetSpheresColor(0.8,0.6,0.6,1);
- Vue1.Dessine_Sphere(X[i],Y[i],Z[i], R[i], 15);
-#endif
- }
- Height = yMax-yMin;
- sectionArea = (xMax-xMin) * (zMax-zMin);
- vTotal = (xMax-xMin) * (yMax-yMin) * (zMax-zMin);
- clock.top("Triangulation");
-
- Tes.Compute();
- clock.top("tesselation");
-
- boundary(yMinId).flowCondition=0;
- boundary(yMaxId).flowCondition=0;
- boundary(yMinId).value=0;
- boundary(yMaxId).value=1;
- defineFictiousCells();
- clock.top("boundaryConditions");
-
- /** INITIALIZATION OF VOLUMES AND PRESSURES **/
- FiniteCellsIterator cellEnd = Tri.finiteCellsEnd();
- for (FiniteCellsIterator cell = Tri.finite_cells_begin(); cell != cellEnd; cell++) {
- cell->info().volume() = ( std::abs ( ( CGT::Tetrahedron ( cell->vertex(0)->point(),cell->vertex(1)->point(),cell->vertex(2)->point(),cell->vertex(3)->point()).volume() ) ) );
- cell->info().dv() = 0;
- }
-
- clock.top("initializing delta_volumes");
-
- /** PERMEABILITY **/
- /** START PERMEABILITY CALCULATION**/
- kFactor = 1;
- computePermeability();
- clock.top("computePermeability");
- /** END PERMEABILITY CALCULATION**/
-
- if(debugOut) cerr << "TOTAL VOID VOLUME: " << vPoral <<endl;
- if(debugOut) cerr << "Porosity = " << vPoralPorosity / vTotalePorosity << endl;
-
- /** STATISTICS **/
- displayStatistics();
- clock.top("displayStatistics");
- /** START GAUSS SEIDEL */
- // Boundary_Conditions ( Tri );
- double pZero = abs((boundary(yMinId).value-boundary(yMaxId).value)/2);
- initializePressures( pZero );
- clock.top("initializePressures");
- gaussSeidel();
- clock.top("gaussSeidel");
- /** END GAUSS SEIDEL */
- const char* file ="Permeability";
- ks = permeameter(boundary(yMinId).value, boundary(yMaxId).value, sectionArea, Height, file);
- clock.top("permeameter");
-
- computeFacetForcesWithCache();
- clock.top("Compute_Forces");
-
- ///*** VUE 3D ***///
-
-#ifdef XVIEW
- Vue1.SetCouleurSegments(0.1,0,1);
- dessineShortTesselation(Vue1, Tes);
- Vue1.Affiche();
-#endif
- if (slipOnLaterals && debugOut) cout << "SLIP CONDITION IS ACTIVATED" << endl;
- else if (debugOut) cout << "NOSLIP CONDITION IS ACTIVATED" << endl;
-// }
- return Tes;
-}
-
-template <class Tesselation>
-Tesselation& FlowBoundingSphere<Tesselation>::loadPositions(int argc, char *argv[ ], char *envp[ ])
-{
- double factor = 1.001;
- VectorR X, Y, Z, R;
- Tesselation& Tes = T[0];
-// RTriangulation& Tri = Tes.Triangulation();
- /** READING SPHERES POSITIONS FROM A TEXT FILE CONTAINING COORDINATES **/
- double x, y, z, r;
- ifstream loadFile(argc==1 ? "cube" : argv[1]); // cree l'objet loadFile de la classe ifstream qui va permettre de lire importFilename
- while (!loadFile.eof()) {
- loadFile >> x >> y >> z >> r;
- X.push_back(x);
- Y.push_back(y);
- Z.push_back(z);
- R.push_back(factor*r);
- nOfSpheres++;
- Rmoy += r;
- xMin = min(xMin,x-r);
- xMax = max(xMax,x+r);
- yMin = min(yMin,y-r);
- yMax = max(yMax,y+r);
- zMin = min(zMin,z-r);
- zMax = max(zMax,z+r);
- }
- Rmoy /= nOfSpheres;
- minPermLength = Rmoy*minLength;
- VertexHandle Vh;
- CellHandle neighbourCell, cell, location;
-
- int V = X.size();
- if (debugOut) cout << "V =" << V << "nOfSpheres = " << nOfSpheres << endl;
- if (debugOut) cout << Tes.Max_id() << endl;
-
- addBoundingPlanes();
- for (int i=0; i<V; i++) {
- int id = Tes.Max_id() +1;
- Vh = Tes.insert(X[i],Y[i],Z[i],R[i],id); /** EMPILEMENT QUELCONQUE **/
-#ifdef XVIEW
- Vue1.SetSpheresColor(0.8,0.6,0.6,1);
- Vue1.Dessine_Sphere(X[i],Y[i],Z[i], R[i], 15);
-#endif
- }
- Height = yMax-yMin;
- sectionArea = (xMax-xMin) * (zMax-zMin);
- vTotal = (xMax-xMin) * (yMax-yMin) * (zMax-zMin);
-
- Tes.Compute();
- defineFictiousCells();
-
- return Tes;
-}
-
template <class Tesselation>
void FlowBoundingSphere<Tesselation>::averageRelativeCellVelocity()
{
@@ -298,10 +105,10 @@
numCells++;
Real totFlowRate = 0;//used to acount for influxes in elements where pressure is imposed
for ( int i=0; i<4; i++ ) if (!Tri.is_infinite(cell->neighbor(i))){
- Vecteur Surfk = cell->info()-cell->neighbor(i)->info();
+ CVector Surfk = cell->info()-cell->neighbor(i)->info();
Real area = sqrt ( Surfk.squared_length() );
Surfk = Surfk/area;
- Vecteur branch = cell->vertex ( facetVertices[i][0] )->point() - cell->info();
+ CVector branch = cell->vertex ( facetVertices[i][0] )->point() - cell->info();
posAvFacet = (Point) cell->info() + ( branch*Surfk ) *Surfk;
facetFlowRate = (cell->info().kNorm())[i] * (cell->info().shiftedP() - cell->neighbor (i)->info().shiftedP());
totFlowRate += facetFlowRate;
@@ -341,7 +148,7 @@
numVertex++;}
vector<Real> volumes;
- vector<CGT::Vecteur> velocityVolumes;
+ vector<CGT::CVector> velocityVolumes;
velocityVolumes.resize(numVertex);
volumes.resize(numVertex);
@@ -364,7 +171,7 @@
double Rz = (zMax-zMin) /20;
CellHandle cellula;
- Vecteur velocity = CGAL::NULL_VECTOR;
+ CVector velocity = CGAL::NULL_VECTOR;
int i=0;
for(double X=xMin+Rx;X<xMax;X+=Rx){
for (double Y=yMin+Ry;Y<yMax;Y+=Ry){
@@ -380,12 +187,10 @@
vector<Real> FlowBoundingSphere<Tesselation>::averageFluidVelocityOnSphere(unsigned int Id_sph)
{
averageRelativeCellVelocity();
- RTriangulation& Tri = T[noCache?(!currentTes):currentTes].Triangulation();
-
- Real volumes; CGT::Vecteur velocityVolumes;
+ RTriangulation& Tri = T[noCache?(!currentTes):currentTes].Triangulation();
+ Real volumes; CGT::CVector velocityVolumes;
vector<Real> result;
- result.resize(3);
-
+ result.resize(3);
velocityVolumes=CGAL::NULL_VECTOR;
volumes=0.f;
@@ -430,7 +235,6 @@
int captures = 6;
double Rz = (zMax-zMin)/intervals;
double Ry = (WallUpy-WallDowny)/captures;
-
double X=(xMax+xMin)/2;
double Y = WallDowny;
double pressure = 0.f;
@@ -485,7 +289,7 @@
void FlowBoundingSphere<Tesselation>::computeFacetForcesWithCache(bool onlyCache)
{
RTriangulation& Tri = T[currentTes].Triangulation();
- Vecteur nullVect(0,0,0);
+ CVector nullVect(0,0,0);
//reset forces
if (!onlyCache) for (FiniteVerticesIterator v = Tri.finite_vertices_begin(); v != Tri.finite_vertices_end(); ++v) v->info().forces=nullVect;
@@ -497,7 +301,7 @@
#endif
CellHandle neighbourCell;
VertexHandle mirrorVertex;
- Vecteur tempVect;
+ CVector tempVect;
//FIXME : Ema, be carefull with this (noCache), it needs to be turned true after retriangulation
if (noCache) {for (VCellIterator cellIt=T[currentTes].cellHandles.begin(); cellIt!=T[currentTes].cellHandles.end(); cellIt++){
CellHandle& cell = *cellIt;
@@ -506,13 +310,13 @@
for (int j=0; j<4; j++) if (!Tri.is_infinite(cell->neighbor(j))) {
neighbourCell = cell->neighbor(j);
- const Vecteur& Surfk = cell->info().facetSurfaces[j];
+ const CVector& 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;
- Vecteur facetNormal = Surfk/area;
- const std::vector<Vecteur>& crossSections = cell->info().facetSphereCrossSections;
- Vecteur fluidSurfk = cell->info().facetSurfaces[j]*cell->info().facetFluidSurfacesRatio[j];
+ CVector facetNormal = Surfk/area;
+ const std::vector<CVector>& crossSections = cell->info().facetSphereCrossSections;
+ CVector 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[boundary(cell->vertex(j)->info().id()).coordinate]);
@@ -522,10 +326,9 @@
cell->info().unitForceVectors[j]=cell->info().unitForceVectors[j]+ tempVect;
}
/// Apply weighted forces f_k=sqRad_k/sumSqRad*f
- Vecteur facetUnitForce = -fluidSurfk*cell->info().solidSurfaces[j][3];
- Vecteur facetForce = cell->info().p()*facetUnitForce;
-
-
+ CVector facetUnitForce = -fluidSurfk*cell->info().solidSurfaces[j][3];
+ CVector facetForce = cell->info().p()*facetUnitForce;
+
for (int y=0; y<3;y++) {
cell->vertex(facetVertices[j][y])->info().forces = cell->vertex(facetVertices[j][y])->info().forces + facetForce*cell->info().solidSurfaces[j][y];
//add to cached value
@@ -556,7 +359,7 @@
if (T[currentTes].vertexHandles[vn]==NULL) continue;
VertexHandle& v = T[currentTes].vertexHandles[vn];
const int& id = v->info().id();
- Vecteur tf (0,0,0);
+ CVector 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);
@@ -565,51 +368,13 @@
#endif
}
if (debugOut) {
- Vecteur totalForce = nullVect;
+ CVector totalForce = nullVect;
for (FiniteVerticesIterator v = Tri.finite_vertices_begin(); v != Tri.finite_vertices_end(); ++v) {
if (!v->info().isFictious) totalForce = totalForce + v->info().forces;
else if (boundary(v->info().id()).flowCondition==1) totalForce = totalForce + v->info().forces; }
cout << "totalForce = "<< totalForce << endl;}
}
-template <class Tesselation>
-void FlowBoundingSphere<Tesselation>::computeTetrahedralForces()
-{
- RTriangulation& Tri = T[currentTes].Triangulation();
- FiniteCellsIterator cellEnd = Tri.finiteCellsEnd();
- Vecteur nullVect(0,0,0);
- bool ref = Tri.finite_cells_begin()->info().isvisited;
-
- for (FiniteVerticesIterator v = Tri.finite_vertices_begin(); v != Tri.finite_vertices_end(); ++v) {
- v->info().forces=nullVect;
- }
-
- CellHandle neighbourCell;
- VertexHandle mirrorVertex;
- for (FiniteCellsIterator cell = Tri.finite_cells_begin(); cell != cellEnd; cell++) {
- for (int j=0; j<4; j++) if (!Tri.is_infinite(cell->neighbor(j)) && cell->neighbor(j)->info().isvisited==ref) {
- neighbourCell = cell->neighbor(j);
- const Vecteur& Surfk = cell->info().facetSurfaces[j];
- /// handle fictious vertex since we can get the projected surface easily here
- if (cell->vertex(j)->info().isFictious) {
- Real projSurf=abs(Surfk[boundary(cell->vertex(j)->info().id()).coordinate]);
- cell->vertex(j)->info().forces = cell->vertex(j)->info().forces -projSurf*boundary(cell->vertex(j)->info().id()).normal*cell->info().p();
- }
- /// handle the opposite fictious vertex (remember each facet is seen only once)
- mirrorVertex = neighbourCell->vertex(Tri.mirror_index(cell,j));
- VertexInfo& info = neighbourCell->vertex(Tri.mirror_index(cell,j))->info();
- if (info.isFictious) {
- Real projSurf=abs(Surfk[boundary(info.id()).coordinate]);
- info.forces = info.forces - projSurf*boundary(info.id()).normal*neighbourCell->info().p();
- }
- /// Apply weighted forces f_k=sqRad_k/sumSqRad*f
- Vecteur facetForce = (neighbourCell->info().p()-cell->info().p())*Surfk*cell->info().solidSurfaces[j][3];
- for (int y=0; y<3;y++) {
- cell->vertex(facetVertices[j][y])->info().forces = cell->vertex(facetVertices[j][y])->info().forces + facetForce*cell->info().solidSurfaces[j][y];
- }
- }
- cell->info().isvisited=!ref;
- }
-}
+
template <class Tesselation>
void FlowBoundingSphere<Tesselation>::applySinusoidalPressure(RTriangulation& Tri, double amplitude, double averagePressure, double loadIntervals)
{
@@ -640,7 +405,7 @@
for (typename VectorCell::iterator cellIt=NewTes.cellHandles.begin(); cellIt!=NewTes.cellHandles.end(); cellIt++){
CellHandle& newCell = *cellIt;
if (newCell->info().Pcondition || newCell->info().isGhost) continue;
- Vecteur center ( 0,0,0 );
+ CVector center ( 0,0,0 );
if (newCell->info().fictious()==0) for ( int k=0;k<4;k++ ) center= center + 0.25* (Tes.vertex(newCell->vertex(k)->info().id())->point()-CGAL::ORIGIN);
else {
Real boundPos=0; int coord=0;
@@ -651,7 +416,7 @@
boundPos=boundary (newCell->vertex(k)->info().id()).p[coord];
}
}
- center=Vecteur(coord==0?boundPos:center[0],coord==1?boundPos:center[1],coord==2?boundPos:center[2]);
+ center=CVector(coord==0?boundPos:center[0],coord==1?boundPos:center[1],coord==2?boundPos:center[2]);
}
oldCell = Tri.locate(Point(center[0],center[1],center[2]));
newCell->info().p() = oldCell->info().shiftedP();
@@ -687,7 +452,7 @@
CellHandle neighbourCell;
- double k=0, distance = 0, radius = 0, viscosity = VISCOSITY;
+ double k=0, distance = 0, radius = 0;
int surfneg=0;
int NEG=0, POS=0, pass=0;
@@ -695,7 +460,6 @@
Real meanK=0, STDEV=0, meanRadius=0, meanDistance=0;
Real infiniteK=1e10;
-
for (VCellIterator cellIt=T[currentTes].cellHandles.begin(); cellIt!=T[currentTes].cellHandles.end(); cellIt++){
CellHandle& cell = *cellIt;
Point& p1 = cell->info();
@@ -712,13 +476,13 @@
Sphere& v1 = W[1]->point();
Sphere& v2 = W[2]->point();
- cell->info().facetSphereCrossSections[j]=Vecteur(
+ cell->info().facetSphereCrossSections[j]=CVector(
W[0]->info().isFictious ? 0 : 0.5*v0.weight()*acos((v1-v0)*(v2-v0)/sqrt((v1-v0).squared_length()*(v2-v0).squared_length())),
W[1]->info().isFictious ? 0 : 0.5*v1.weight()*acos((v0-v1)*(v2-v1)/sqrt((v1-v0).squared_length()*(v2-v1).squared_length())),
W[2]->info().isFictious ? 0 : 0.5*v2.weight()*acos((v0-v2)*(v1-v2)/sqrt((v1-v2).squared_length()*(v2-v0).squared_length())));
pass+=1;
- Vecteur l = p1 - p2;
+ CVector l = p1 - p2;
distance = sqrt(l.squared_length());
if (!RAVERAGE) radius = 2* computeHydraulicRadius(cell, j);
else radius = (computeEffectiveRadius(cell, j)+computeEquivalentRadius(cell,j))*0.5;
@@ -730,28 +494,23 @@
Real fluidArea=0;
if (distance!=0) {
if (minPermLength>0 && distance_correction) distance=max(minPermLength,distance);
- const Vecteur& Surfk = cell->info().facetSurfaces[j];
+ const CVector& Surfk = cell->info().facetSurfaces[j];
Real area = sqrt(Surfk.squared_length());
- const Vecteur& crossSections = cell->info().facetSphereCrossSections[j];
- Real S0=0;
- S0=checkSphereFacetOverlap(v0,v1,v2);
- if (S0==0) S0=checkSphereFacetOverlap(v1,v2,v0);
- if (S0==0) S0=checkSphereFacetOverlap(v2,v0,v1);
- //take absolute value, since in rare cases the surface can be negative (overlaping spheres)
- fluidArea=abs(area-crossSections[0]-crossSections[1]-crossSections[2]+S0);
- cell->info().facetFluidSurfacesRatio[j]=fluidArea/area;
- k=(fluidArea * pow(radius,2)) / (8*viscosity*distance);
-// } else {
-// cout << "WARNING! if !areaR2Permeability, facetFluidSurfacesRatio will not be defined correctly. Don't use that."<<endl;
-// k = (M_PI * pow(radius,4)) / (8*viscosity*distance);}
-
+ const CVector& crossSections = cell->info().facetSphereCrossSections[j];
+ Real S0=0;
+ S0=checkSphereFacetOverlap(v0,v1,v2);
+ if (S0==0) S0=checkSphereFacetOverlap(v1,v2,v0);
+ if (S0==0) S0=checkSphereFacetOverlap(v2,v0,v1);
+ //take absolute value, since in rare cases the surface can be negative (overlaping spheres)
+ fluidArea=abs(area-crossSections[0]-crossSections[1]-crossSections[2]+S0);
+ cell->info().facetFluidSurfacesRatio[j]=fluidArea/area;
+ k=(fluidArea * pow(radius,2)) / (8*viscosity*distance);
meanDistance += distance;
meanRadius += radius;
meanK += k*kFactor;
if (k<0 && debugOut) {surfneg+=1;
- cout<<"__ k<0 __"<<k<<" "<<" fluidArea "<<fluidArea<<" area "<<area<<" "<<crossSections[0]<<" "<<crossSections[1]<<" "<<crossSections[2] <<" "<<W[0]->info().id()<<" "<<W[1]->info().id()<<" "<<W[2]->info().id()<<" "<<p1<<" "<<p2<<" test "<<endl;}
-
+ cout<<"__ k<0 __"<<k<<" "<<" fluidArea "<<fluidArea<<" area "<<area<<" "<<crossSections[0]<<" "<<crossSections[1]<<" "<<crossSections[2] <<" "<<W[0]->info().id()<<" "<<W[1]->info().id()<<" "<<W[2]->info().id()<<" "<<p1<<" "<<p2<<" test "<<endl;}
} else {cout <<"infinite K1!"<<endl; k = infiniteK;}//Will be corrected in the next loop
(cell->info().kNorm())[j]= k*kFactor;
@@ -789,7 +548,6 @@
}
}
if (debugOut) cout << "PassKcorrect = " << pass << endl;
-
if (debugOut) cout << "POS = " << POS << " NEG = " << NEG << " pass = " << pass << endl;
// A loop to compute the standard deviation of the local K distribution, and use it to include/exclude K values higher then (meanK +/- K_opt_factor*STDEV)
@@ -808,11 +566,9 @@
}cell->info().isvisited = !ref;
}
STDEV = sqrt(STDEV/pass);
- if (debugOut) cout << "PassSTDEV = " << pass << endl;
- cout << "STATISTIC K" << endl;
+ if (debugOut) cout << "PassSTDEV = " << pass << endl << "STATISTIC K" << endl;
double k_min = 0, k_max = meanK + KOptFactor*STDEV;
- cout << "Kmoy = " << meanK << " Standard Deviation = " << STDEV << endl;
- cout << "kmin = " << k_min << " kmax = " << k_max << endl;
+ cout << "Kmoy = " << meanK << " Standard Deviation = " << STDEV << endl<< "kmin = " << k_min << " kmax = " << k_max << endl;
ref = Tri.finite_cells_begin()->info().isvisited;
pass=0;
for (FiniteCellsIterator cell = Tri.finite_cells_begin(); cell != cellEnd; cell++) {
@@ -830,8 +586,6 @@
}
if (debugOut) cout << "PassKopt = " << pass << endl;
}
-
-
if (debugOut) {
FiniteVerticesIterator verticesEnd = Tri.finite_vertices_end();
Real Vgrains = 0;
@@ -873,7 +627,7 @@
|| ((f_it->first->info().index > f_it->first->neighbor(f_it->second)->info().index) && f_it->first->neighbor(f_it->second)->info().isGhost)
|| f_it->first->info().index == 0 || f_it->first->neighbor(f_it->second)->info().index == 0) continue;
vector<double> rn;
- const Vecteur& normal = f_it->first->info().facetSurfaces[f_it->second];
+ const CVector& normal = f_it->first->info().facetSurfaces[f_it->second];
if (!normal[0] && !normal[1] && !normal[2]) continue;
rn.push_back(computeEffectiveRadius(f_it->first, f_it->second));
rn.push_back(normal[0]);
@@ -891,11 +645,11 @@
RTriangulation& Tri = T[currentTes].Triangulation();
if (Tri.is_infinite(cell->neighbor(j))) return 0;
- Vecteur B = cell->vertex(facetVertices[j][1])->point().point()-cell->vertex(facetVertices[j][0])->point().point();
- Vecteur x = B/sqrt(B.squared_length());
- Vecteur C = cell->vertex(facetVertices[j][2])->point().point()-cell->vertex(facetVertices[j][0])->point().point();
- Vecteur z = CGAL::cross_product(x,C);
- Vecteur y = CGAL::cross_product(x,z);
+ CVector B = cell->vertex(facetVertices[j][1])->point().point()-cell->vertex(facetVertices[j][0])->point().point();
+ CVector x = B/sqrt(B.squared_length());
+ CVector C = cell->vertex(facetVertices[j][2])->point().point()-cell->vertex(facetVertices[j][0])->point().point();
+ CVector z = CGAL::cross_product(x,C);
+ CVector y = CGAL::cross_product(x,z);
y = y/sqrt(y.squared_length());
double b1[2]; b1[0] = B*x; b1[1] = B*y;
@@ -1025,8 +779,6 @@
double compFlowFactor=0;
vector<Real> previousP;
previousP.resize(Tri.number_of_finite_cells());
- double tolerance = TOLERANCE;
- double relax = RELAX;
const int num_threads=1;
bool compressible= (fluidBulkModulus>0);
#ifdef GS_OPEN_MP
@@ -1205,7 +957,7 @@
}}
double density = 1;
- double viscosity = VISCOSITY;
+ double viscosity = viscosity;
double gravity = 1;
double Vdarcy = Q1/Section;
double DeltaP = abs(PInf-PSup);
@@ -1443,62 +1195,7 @@
consFile << endl;
consFile.close();
}
-template <class Tesselation>
-void FlowBoundingSphere<Tesselation>::comsolField()
-{
- //Compute av. velocity first, because in the following permeabilities will be overwritten with "junk" (in fact velocities from comsol)
- averageRelativeCellVelocity();
-
- RTriangulation& Tri = T[noCache?(!currentTes):currentTes].Triangulation();
- CellHandle c;
- ifstream loadFile("vx_grid_03_07_ns.txt");
- ifstream loadFileY("vy_grid_03_07_ns.txt");
- ifstream loadFileZ("vz_grid_03_07_ns.txt");
- int Nx=100; int Ny=10; int Nz=100;
- std::ofstream consFile("velComp",std::ios::out);
-
- FiniteCellsIterator cellEnd = Tri.finiteCellsEnd();
- for (FiniteCellsIterator cell = Tri.finite_cells_begin(); cell != cellEnd; cell++){
- cell->info().dv()=0;
- cell->info().modulePermeability[0]=cell->info().modulePermeability[1]=cell->info().modulePermeability[2]=0;
- }
-
- vector<Real> X, Y, Z;
- Real buffer;
- for (int xi=0;xi<Nx;xi++) {loadFile >> buffer; X.push_back(buffer); loadFileY >> buffer; loadFileZ >> buffer;}
- for (int yi=0;yi<Ny;yi++) {loadFile >> buffer; Y.push_back(buffer); loadFileY >> buffer; loadFileZ >> buffer;}
- for (int zi=0;zi<Nz;zi++) {loadFile >> buffer; Z.push_back(buffer); loadFileY >> buffer; loadFileZ >> buffer;}
-
- Real vx, vy, vz;
- Real meanCmsVel=0; int totCmsPoints = 0;
- for (int zi=0;zi<Nz;zi++)
- for (int yi=0;yi<Ny;yi++)
- for (int xi=0;xi<Nx;xi++) {
- loadFile >> vx; loadFileY >> vy; loadFileZ >> vz;
- if (!isInsideSphere(X[xi], Y[yi], Z[zi]) && vx!=0) {
- c = Tri.locate(Point(X[xi], Y[yi], Z[zi]));
- c->info().modulePermeability[0]+=vx;
- c->info().modulePermeability[1]+=vy;
- c->info().modulePermeability[2]+=vz;
- c->info().dv()+=1;
- meanCmsVel+=vy; totCmsPoints++;}
- }
- int kk=0;
- Vecteur diff;
- for (FiniteCellsIterator cell = Tri.finite_cells_begin(); cell != cellEnd && kk<10000; cell++){
- if (cell->info().fictious() || cell->info().dv()<60) continue;
- for (int k=0;k<3;k++) cell->info().modulePermeability[k]/=cell->info().dv();
- cerr << cell->info().modulePermeability[0]<<" "<< cell->info().modulePermeability[1]<<" "<< cell->info().modulePermeability[2]<<" "<<cell->info().dv()<<" "<< cell->info().averageVelocity()<<endl;
- Real m=sqrt(pow(cell->info().modulePermeability[0],2)+pow(cell->info().modulePermeability[1],2)+pow(cell->info().modulePermeability[2],2));
- Vecteur comFlow (cell->info().modulePermeability[0],cell->info().modulePermeability[1],cell->info().modulePermeability[2]);
- Real angle=asin(sqrt(cross_product(comFlow,cell->info().averageVelocity()).squared_length())/(sqrt(comFlow.squared_length())*sqrt(cell->info().averageVelocity().squared_length())));
- cerr<<"norms : "<<m<<" vs. "<<sqrt(cell->info().averageVelocity().squared_length())<<" angle "<<180*angle/3.1415<<endl;
- consFile<<m<<" "<<sqrt(cell->info().averageVelocity().squared_length())<<" "<<180*angle/3.1415<<endl;
- diff = diff + (comFlow - cell->info().averageVelocity());
- kk++;
- }
- cerr << "meanCmsVel "<<meanCmsVel/totCmsPoints<<" mean diff "<<diff/kk<<endl;
-}
+
template <class Tesselation>
void FlowBoundingSphere<Tesselation>::computeEdgesSurfaces()
{
@@ -1542,7 +1239,7 @@
template <class Tesselation>
Vector3r FlowBoundingSphere<Tesselation>::computeViscousShearForce(const Vector3r& deltaV, const int& edge_id, const Real& Rh)
{
- Vector3r tau = deltaV*VISCOSITY/Rh;
+ Vector3r tau = deltaV*viscosity/Rh;
return tau * edgeSurfaces[edge_id];
}
@@ -1550,7 +1247,7 @@
Vector3r FlowBoundingSphere<Tesselation>::computeShearLubricationForce(const Vector3r& deltaShearV, const Real& dist, const int& edge_id, const Real& eps, const Real& centerDist, const Real& meanRad )
{
Real d = max(dist,0.) + 2.*eps*meanRad;
- Vector3r viscLubF = 0.5 * Mathr::PI * VISCOSITY * (-2*meanRad + centerDist*log(centerDist/d)) * deltaShearV;
+ Vector3r viscLubF = 0.5 * Mathr::PI * viscosity * (-2*meanRad + centerDist*log(centerDist/d)) * deltaShearV;
return viscLubF;
}
@@ -1558,7 +1255,7 @@
Vector3r FlowBoundingSphere<Tesselation>::computePumpTorque(const Vector3r& deltaShearAngV, const Real& dist, const int& edge_id, const Real& eps, const Real& meanRad )
{
Real d = max(dist,0.) + 2.*eps*meanRad;
- Vector3r viscPumpC = Mathr::PI * VISCOSITY * pow(meanRad,3) *(3./20. * log(meanRad/d) + 63./500. * (d/meanRad) * log(meanRad/d)) * deltaShearAngV;
+ Vector3r viscPumpC = Mathr::PI * viscosity * pow(meanRad,3) *(3./20. * log(meanRad/d) + 63./500. * (d/meanRad) * log(meanRad/d)) * deltaShearAngV;
return viscPumpC;
}
@@ -1566,7 +1263,7 @@
Vector3r FlowBoundingSphere<Tesselation>::computeTwistTorque(const Vector3r& deltaNormAngV, const Real& dist, const int& edge_id, const Real& eps, const Real& meanRad )
{
Real d = max(dist,0.) + 2.*eps*meanRad;
- Vector3r twistC = Mathr::PI * VISCOSITY * pow(meanRad,2) * d * log(meanRad/d) * deltaNormAngV;
+ Vector3r twistC = Mathr::PI * viscosity * pow(meanRad,2) * d * log(meanRad/d) * deltaNormAngV;
return twistC;
}
@@ -1579,12 +1276,12 @@
if (stiffness>0) {
const Real k = stiffness*meanRad;
Real prevForce = edgeNormalLubF[edge_id];
- Real instantVisc = 1.5*Mathr::PI*pow(meanRad,2)*VISCOSITY/(d-prevForce/k);
+ Real instantVisc = 1.5*Mathr::PI*pow(meanRad,2)*viscosity/(d-prevForce/k);
Real normLubF = instantVisc*(deltaNormV + prevForce/(k*dt))/(1+instantVisc/(k*dt));
edgeNormalLubF[edge_id]=normLubF;
return normLubF;
} else {
- Real normLubF = (1.5*Mathr::PI*pow(meanRad,2)* VISCOSITY* deltaNormV)/d;
+ Real normLubF = (1.5*Mathr::PI*pow(meanRad,2)* viscosity* deltaNormV)/d;
return normLubF;
}
}
=== modified file 'lib/triangulation/FlowBoundingSphereLinSolv.hpp'
--- lib/triangulation/FlowBoundingSphereLinSolv.hpp 2014-03-21 18:45:24 +0000
+++ lib/triangulation/FlowBoundingSphereLinSolv.hpp 2014-03-21 18:45:24 +0000
@@ -47,8 +47,8 @@
using FlowType::yMinId;
using FlowType::yMaxId;
using FlowType::debugOut;
- using FlowType::TOLERANCE;
- using FlowType::RELAX;
+ using FlowType::tolerance;
+ using FlowType::relax;
using FlowType::fluidBulkModulus;
using FlowType::reApplyBoundaryConditions;
using FlowType::pressureChanged;
=== modified file 'lib/triangulation/FlowBoundingSphereLinSolv.ipp'
--- lib/triangulation/FlowBoundingSphereLinSolv.ipp 2014-03-21 18:45:24 +0000
+++ lib/triangulation/FlowBoundingSphereLinSolv.ipp 2014-03-21 18:45:24 +0000
@@ -9,8 +9,6 @@
// #define XVIEW
#include "FlowBoundingSphereLinSolv.hpp"//include after #define XVIEW
-// #include "def_types.h"
-// #include "def_flow_types.h"
#include "CGAL/constructions/constructions_on_weighted_points_cartesian_3.h"
#include <CGAL/Width_3.h>
#include <iostream>
@@ -426,8 +424,6 @@
int j = 0;
double dp_max, p_max, sum_p, p_moy, dp_moy, sum_dp;
- double tolerance = TOLERANCE;
- double relax = RELAX;
#ifdef GS_OPEN_MP
const int num_threads=1;
=== modified file 'lib/triangulation/KinematicLocalisationAnalyser.cpp'
--- lib/triangulation/KinematicLocalisationAnalyser.cpp 2012-03-12 19:57:25 +0000
+++ lib/triangulation/KinematicLocalisationAnalyser.cpp 2014-03-21 18:45:24 +0000
@@ -172,7 +172,7 @@
&& TS1->inside(T.segment(*ed_it).source())
&& TS1->inside(T.segment(*ed_it).target())) {
Segment s = T.segment(*ed_it);
- Vecteur v = s.to_vector();
+ CVector v = s.to_vector();
Real ny = abs(v.y()/sqrt(s.squared_length()));
if (Nymin < ny && ny <= Nymax) filteredList.push_back(ed_it);
@@ -272,7 +272,7 @@
for (Edge_iterator ed_it = T.edges_begin(); ed_it != ed_end; ++ed_it) {
if (!T.is_infinite(*ed_it)) {
Segment s = T.segment(*ed_it);
- Vecteur v = s.to_vector();
+ CVector v = s.to_vector();
Real xx = abs(v.z()/sqrt(s.squared_length()));
if (xx>0.95) edges.push_back(ed_it);
@@ -284,7 +284,7 @@
for (Edge_iterator ed_it = T.edges_begin(); ed_it != ed_end; ++ed_it) {
if (!T.is_infinite(*ed_it)) {
Segment s = T.segment(*ed_it);
- Vecteur v = s.to_vector();
+ CVector v = s.to_vector();
Real xx = abs(v.z()/sqrt(s.squared_length()));
if (xx<0.05) edges.push_back(ed_it);
@@ -296,7 +296,7 @@
for (Edge_iterator ed_it = T.edges_begin(); ed_it != ed_end; ++ed_it) {
if (!T.is_infinite(*ed_it)) {
Segment s = T.segment(*ed_it);
- Vecteur v = s.to_vector();
+ CVector v = s.to_vector();
Real xx = abs(v.z()/sqrt(s.squared_length()));
if (xx>0.65 && xx<0.75) edges.push_back(ed_it);
@@ -369,7 +369,7 @@
RTriangulation& T = state.tesselation().Triangulation();
Edge_iterator ed_end = T.edges_end();
Tenseur_sym3 Tens;
- Vecteur v;
+ CVector v;
Segment s;
for (Edge_iterator ed_it = T.edges_begin(); ed_it != ed_end; ++ed_it) {
if (!T.is_infinite(*ed_it)) {
@@ -397,7 +397,7 @@
state)
{
Tenseur_sym3 Tens;
- Vecteur v;
+ CVector v;
TriaxialState::ContactIterator cend = state.contacts_end();
for (TriaxialState::ContactIterator cit = state.contacts_begin();
@@ -438,7 +438,7 @@
vector<Real> Un_values;
Un_values.resize(edges.size());
Real UNmin(100000), UNmax(-100000);
- Vecteur branch, U;
+ CVector branch, U;
Real Un;
Vertex_handle Vh1, Vh2;
vector<Edge_iterator>::iterator ed_end = edges.end();
@@ -598,7 +598,7 @@
RTriangulation& T = (*TS1).tesselation().Triangulation();
Segment s;
- Vecteur v;
+ CVector v;
for (Edge_iterator ed_it = T.edges_begin(); ed_it != T.edges_end(); ed_it++) {
if (!T.is_infinite(*ed_it)) {
s = T.segment(*ed_it);
@@ -717,22 +717,22 @@
return output_file;
}
-Vecteur KinematicLocalisationAnalyser::Deplacement(Finite_cells_iterator cell, int facet)
+CVector KinematicLocalisationAnalyser::Deplacement(Finite_cells_iterator cell, int facet)
{
- Vecteur v(0.f, 0.f, 0.f);
+ CVector v(0.f, 0.f, 0.f);
int id;// ident. de la particule
- Vecteur fixedPoint = 0.5*((TS0->box.base-CGAL::ORIGIN)+(TS0->box.sommet-CGAL::ORIGIN));
+ CVector fixedPoint = 0.5*((TS0->box.base-CGAL::ORIGIN)+(TS0->box.sommet-CGAL::ORIGIN));
for (int i=0; i<4; i++) {
// char msg [256];
if (i!=facet) {
id = cell->vertex(i)->info().id();
- Vecteur meanFieldDisp =Vecteur(TS0->grain(id).sphere.point().x(), TS0->grain(id).sphere.point().y(), TS0->grain(id).sphere.point().z())-fixedPoint;
+ CVector meanFieldDisp =CVector(TS0->grain(id).sphere.point().x(), TS0->grain(id).sphere.point().y(), TS0->grain(id).sphere.point().z())-fixedPoint;
if (1){//fluctuations
- meanFieldDisp = Vecteur(
+ meanFieldDisp = CVector(
meanFieldDisp[0]*Delta_epsilon(0,0),
meanFieldDisp[1]*Delta_epsilon(1,1),
meanFieldDisp[2]*Delta_epsilon(2,2));
- } else meanFieldDisp=Vecteur(0,0,0);
+ } else meanFieldDisp=CVector(0,0,0);
if (consecutive) v = v + TS1->grain(id).translation-meanFieldDisp;
else v = v + (TS1->grain(id).sphere.point() - TS0->grain(id).sphere.point()-meanFieldDisp);
}
@@ -741,9 +741,9 @@
return v;
}
-void KinematicLocalisationAnalyser::Grad_u(Finite_cells_iterator cell, int facet, Vecteur &V, Tenseur3& T)
+void KinematicLocalisationAnalyser::Grad_u(Finite_cells_iterator cell, int facet, CVector &V, Tenseur3& T)
{
- Vecteur S = cross_product((cell->vertex(l_vertices[facet][1])->point())
+ CVector S = cross_product((cell->vertex(l_vertices[facet][1])->point())
- (cell->vertex(l_vertices[facet][0])->point()),
(cell->vertex(l_vertices[facet][2])->point()) -
(cell->vertex(l_vertices[facet][1])->point())) /2.f;
@@ -754,7 +754,7 @@
Tenseur3& T, bool vol_divide)// Calcule le gradient de d�p.
{
T.reset();
- Vecteur v;
+ CVector v;
for (int facet=0; facet<4; facet++) {
v = Deplacement(cell, facet);
Grad_u(cell, facet, v, T);
=== modified file 'lib/triangulation/KinematicLocalisationAnalyser.hpp'
--- lib/triangulation/KinematicLocalisationAnalyser.hpp 2012-01-20 17:31:56 +0000
+++ lib/triangulation/KinematicLocalisationAnalyser.hpp 2014-03-21 18:45:24 +0000
@@ -78,7 +78,7 @@
void SetDisplacementIncrements (void);
///Add surface*displacement to T
- void Grad_u ( Finite_cells_iterator cell, int facet, Vecteur &V, Tenseur3& T );
+ void Grad_u ( Finite_cells_iterator cell, int facet, CVector &V, Tenseur3& T );
///Compute grad_u in cell (by default, T= average grad_u in cell, if !vol_divide, T=grad_u*volume
void Grad_u ( Finite_cells_iterator cell, Tenseur3& T, bool vol_divide=true );
/// Compute grad_u for all particles, by summing grad_u of all adjaent cells using current states
@@ -88,8 +88,8 @@
///Compute porisity from cumulated spheres volumes and positions of boxes
Real ComputeMacroPorosity (void );
- Vecteur Deplacement ( Cell_handle cell ); //donne le d�placement d'un sommet de voronoi
- Vecteur Deplacement ( Finite_cells_iterator cell, int facet ); //mean displacement on a facet
+ CVector Deplacement ( Cell_handle cell ); //donne le d�placement d'un sommet de voronoi
+ CVector Deplacement ( Finite_cells_iterator cell, int facet ); //mean displacement on a facet
// Calcul du tenseur d'orientation des voisins
//Tenseur_sym3 Orientation_voisins (Tesselation& Tes);
=== modified file 'lib/triangulation/Network.hpp'
--- lib/triangulation/Network.hpp 2014-03-21 18:45:24 +0000
+++ lib/triangulation/Network.hpp 2014-03-21 18:45:24 +0000
@@ -26,7 +26,7 @@
struct Boundary
{
Point p;//position
- Vecteur normal;//orientation
+ CVector normal;//orientation
Vector3r velocity;//motion
int coordinate;//the axis perpendicular to the boundary
bool flowCondition;//flowCondition=0, pressure is imposed // flowCondition=1, flow is imposed
@@ -61,14 +61,14 @@
int vtkInfiniteVertices, vtkInfiniteCells, num_particles;
void addBoundingPlanes();
- void addBoundingPlane (Vecteur Normal, int id_wall);
- void addBoundingPlane (Real center[3], double thickness, Vecteur Normal, int id_wall );
+ void addBoundingPlane (CVector Normal, int id_wall);
+ void addBoundingPlane (Real center[3], double thickness, CVector Normal, int id_wall );
void defineFictiousCells( );
int detectFacetFictiousVertices (CellHandle& cell, int& j);
double volumeSolidPore (const CellHandle& cell);
- double volumeSingleFictiousPore(const VertexHandle& SV1, const VertexHandle& SV2, const VertexHandle& SV3, const Point& PV1, const Point& PV2, Vecteur& facetSurface);
- double volumeDoubleFictiousPore(const VertexHandle& SV1, const VertexHandle& SV2, const VertexHandle& SV3, const Point& PV1, const Point& PV2, Vecteur& facetSurface);
+ double volumeSingleFictiousPore(const VertexHandle& SV1, const VertexHandle& SV2, const VertexHandle& SV3, const Point& PV1, const Point& PV2, CVector& facetSurface);
+ double volumeDoubleFictiousPore(const VertexHandle& SV1, const VertexHandle& SV2, const VertexHandle& SV3, const Point& PV1, const Point& PV2, CVector& facetSurface);
double sphericalTriangleVolume(const Sphere& ST1, const Point& PT1, const Point& PT2, const Point& PT3);
double fastSphericalTriangleArea(const Sphere& STA1, const Point& STA2, const Point& STA3, const Point& PTA1);
@@ -79,8 +79,8 @@
double surfaceSolidPore( CellHandle cell, int j, bool slipOnLaterals, bool reuseFacetData=false);
double sphericalTriangleArea ( Sphere STA1, Sphere STA2, Sphere STA3, Point PTA1 );
- Vecteur surfaceDoubleFictiousFacet(VertexHandle fSV1, VertexHandle fSV2, VertexHandle SV3);
- Vecteur surfaceSingleFictiousFacet(VertexHandle fSV1, VertexHandle SV2, VertexHandle SV3);
+ CVector surfaceDoubleFictiousFacet(VertexHandle fSV1, VertexHandle fSV2, VertexHandle SV3);
+ CVector surfaceSingleFictiousFacet(VertexHandle fSV1, VertexHandle SV2, VertexHandle SV3);
double surfaceSolidDoubleFictiousFacet(VertexHandle SV1, VertexHandle SV2, VertexHandle SV3);
double surfaceSolidFacet(Sphere ST1, Sphere ST2, Sphere ST3);
=== modified file 'lib/triangulation/Network.ipp'
--- lib/triangulation/Network.ipp 2014-03-21 18:45:24 +0000
+++ lib/triangulation/Network.ipp 2014-03-21 18:45:24 +0000
@@ -104,7 +104,7 @@
}
template<class Tesselation>
-double Network<Tesselation>::volumeSingleFictiousPore(const VertexHandle& SV1, const VertexHandle& SV2, const VertexHandle& SV3, const Point& PV1, const Point& PV2, Vecteur& facetSurface)
+double Network<Tesselation>::volumeSingleFictiousPore(const VertexHandle& SV1, const VertexHandle& SV2, const VertexHandle& SV3, const Point& PV1, const Point& PV2, CVector& facetSurface)
{
double A [3], B[3];
@@ -138,7 +138,7 @@
}
template<class Tesselation>
-double Network<Tesselation>::volumeDoubleFictiousPore(const VertexHandle& SV1, const VertexHandle& SV2, const VertexHandle& SV3, const Point& PV1, const Point& PV2, Vecteur& facetSurface)
+double Network<Tesselation>::volumeDoubleFictiousPore(const VertexHandle& SV1, const VertexHandle& SV2, const VertexHandle& SV3, const Point& PV1, const Point& PV2, CVector& facetSurface)
{
double A [3], B[3];
@@ -177,9 +177,6 @@
double Network<Tesselation>::fastSphericalTriangleArea(const Sphere& STA1, const Point& STA2, const Point& STA3, const Point& PTA1)
{
using namespace CGAL;
-#ifndef FAST
- return sphericalTriangleArea(STA1, STA2, STA3, PTA1);
-#endif
double rayon2 = STA1.weight();
if (rayon2 == 0.0) return 0.0;
return rayon2 * fastSolidAngle(STA1,STA2,STA3,PTA1);
@@ -191,9 +188,9 @@
double rayon = STA1.weight();
if ( rayon == 0.0 ) return 0.0;
- Vecteur v12 = STA2.point() - STA1.point();
- Vecteur v13 = STA3.point() - STA1.point();
- Vecteur v14 = PTA1 - STA1.point();
+ CVector v12 = STA2.point() - STA1.point();
+ CVector v13 = STA3.point() - STA1.point();
+ CVector v14 = PTA1 - STA1.point();
double norme12 = ( v12.squared_length() );
double norme13 = ( v13.squared_length() );
@@ -334,7 +331,7 @@
Ssolid1n = fastSphericalTriangleArea(SV3->point(), BB, p1, p2);
cell->info().solidSurfaces[j][facetRe1]=Ssolid1+Ssolid1n;
//area vector of triangle (p1,sphere,p2)
- Vecteur p1p2v1Surface = 0.5*CGAL::cross_product(p1-p2,SV3->point()-p2);
+ CVector p1p2v1Surface = 0.5*CGAL::cross_product(p1-p2,SV3->point()-p2);
if (bi1.flowCondition && ! slipOnLaterals) {
//projection on boundary 1
Ssolid2 = abs(p1p2v1Surface[bi1.coordinate]);
@@ -361,7 +358,7 @@
}
template<class Tesselation>
-Vecteur Network<Tesselation>::surfaceDoubleFictiousFacet(VertexHandle fSV1, VertexHandle fSV2, VertexHandle SV3)
+CVector Network<Tesselation>::surfaceDoubleFictiousFacet(VertexHandle fSV1, VertexHandle fSV2, VertexHandle SV3)
{
//This function is correct only with axis-aligned boundaries
const Boundary &bi1 = boundary(fSV1->info().id());
@@ -370,16 +367,16 @@
double surf [3] = {1,1,1};
surf[bi1.coordinate]=0;
surf[bi2.coordinate]=0;
- return area*Vecteur(surf[0],surf[1],surf[2]);
+ return area*CVector(surf[0],surf[1],surf[2]);
}
template<class Tesselation>
-Vecteur Network<Tesselation>::surfaceSingleFictiousFacet(VertexHandle fSV1, VertexHandle SV2, VertexHandle SV3)
+CVector Network<Tesselation>::surfaceSingleFictiousFacet(VertexHandle fSV1, VertexHandle SV2, VertexHandle SV3)
{
//This function is correct only with axis-aligned boundaries
const Boundary &bi1 = boundary(fSV1->info().id());
// const Boundary &bi2 = boundary ( fSV2->info().id() );
- Vecteur mean_height = (bi1.p[bi1.coordinate]-0.5*(SV3->point()[bi1.coordinate]+SV2->point()[bi1.coordinate]))*bi1.normal;
+ CVector mean_height = (bi1.p[bi1.coordinate]-0.5*(SV3->point()[bi1.coordinate]+SV2->point()[bi1.coordinate]))*bi1.normal;
return CGAL::cross_product(mean_height,SV3->point()-SV2->point());
}
@@ -418,8 +415,8 @@
double squared_radius = ST1.weight();
- Vecteur v12 = ST2.point() - ST1.point();
- Vecteur v13 = ST3.point() - ST1.point();
+ CVector v12 = ST2.point() - ST1.point();
+ CVector v13 = ST3.point() - ST1.point();
double norme12 = v12.squared_length();
double norme13 = v13.squared_length();
@@ -454,18 +451,18 @@
idOffset = Tes.Max_id() +1;//so that boundaries[vertex->id - offset] gives the ordered boundaries (also see function Boundary& boundary(int b))
- addBoundingPlane (Vecteur(0,1,0) , yMinId);
- addBoundingPlane (Vecteur(0,-1,0) , yMaxId);
- addBoundingPlane (Vecteur(-1,0,0) , xMaxId);
- addBoundingPlane (Vecteur(1,0,0) , xMinId);
- addBoundingPlane (Vecteur(0,0,1) , zMinId);
- addBoundingPlane (Vecteur(0,0,-1) , zMaxId);
+ addBoundingPlane (CVector(0,1,0) , yMinId);
+ addBoundingPlane (CVector(0,-1,0) , yMaxId);
+ addBoundingPlane (CVector(-1,0,0) , xMaxId);
+ addBoundingPlane (CVector(1,0,0) , xMinId);
+ addBoundingPlane (CVector(0,0,1) , zMinId);
+ addBoundingPlane (CVector(0,0,-1) , zMaxId);
// addBoundingPlanes(true);
}
template<class Tesselation>
-void Network<Tesselation>::addBoundingPlane (Vecteur Normal, int id_wall)
+void Network<Tesselation>::addBoundingPlane (CVector Normal, int id_wall)
{
// Tesselation& Tes = T[currentTes];
//FIXME: pre-condition: the normal is axis-aligned
@@ -482,7 +479,7 @@
}
template<class Tesselation>
-void Network<Tesselation>::addBoundingPlane (Real center[3], double thickness, Vecteur Normal, int id_wall )
+void Network<Tesselation>::addBoundingPlane (Real center[3], double thickness, CVector Normal, int id_wall )
{
Tesselation& Tes = T[currentTes];
@@ -537,9 +534,9 @@
// double rayon = STA1.weight();
// if ( rayon == 0.0 ) return 0.0;
//
-// Vecteur v12 = STA2.point() - STA1.point();
-// Vecteur v13 = STA3.point() - STA1.point();
-// Vecteur v14 = PTA1 - STA1.point();
+// CVector v12 = STA2.point() - STA1.point();
+// CVector v13 = STA3.point() - STA1.point();
+// CVector v14 = PTA1 - STA1.point();
//
// double norme12 = ( v12.squared_length() );
// double norme13 = ( v13.squared_length() );
=== modified file 'lib/triangulation/PeriodicFlow.cpp'
--- lib/triangulation/PeriodicFlow.cpp 2014-03-21 18:45:24 +0000
+++ lib/triangulation/PeriodicFlow.cpp 2014-03-21 18:45:24 +0000
@@ -2,10 +2,6 @@
#include"PeriodicFlow.hpp"
-#define FAST
-#define TESS_BASED_FORCES
-#define FACET_BASED_FORCES 1
-
#ifdef YADE_OPENMP
// #define GS_OPEN_MP //It should never be defined if Yade is not using openmp
#endif
@@ -20,7 +16,7 @@
for (VectorCell::iterator cellIt=NewTes.cellHandles.begin(); cellIt!=NewTes.cellHandles.end(); cellIt++){
CellHandle& newCell = *cellIt;
if (newCell->info().Pcondition || newCell->info().isGhost) continue;
- Vecteur center ( 0,0,0 );
+ CVector center ( 0,0,0 );
if (newCell->info().fictious()==0) for ( int k=0;k<4;k++ ) center= center + 0.25* (Tes.vertex(newCell->vertex(k)->info().id())->point()-CGAL::ORIGIN);
else {
Real boundPos=0; int coord=0;
@@ -31,7 +27,7 @@
boundPos=boundary (newCell->vertex(k)->info().id()).p[coord];
}
}
- center=Vecteur(coord==0?boundPos:center[0],coord==1?boundPos:center[1],coord==2?boundPos:center[2]);
+ center=CVector(coord==0?boundPos:center[0],coord==1?boundPos:center[1],coord==2?boundPos:center[2]);
}
oldCell = Tri.locate(Point(center[0],center[1],center[2]));
newCell->info().p() = oldCell->info().shiftedP();
@@ -44,8 +40,8 @@
void PeriodicFlow::computeFacetForcesWithCache(bool onlyCache)
{
RTriangulation& Tri = T[currentTes].Triangulation();
- Vecteur nullVect(0,0,0);
- static vector<Vecteur> oldForces;
+ CVector nullVect(0,0,0);
+ static vector<CVector> oldForces;
if (oldForces.size()<=Tri.number_of_vertices()) oldForces.resize(Tri.number_of_vertices()+1);
//reset forces
for (FiniteVerticesIterator v = Tri.finite_vertices_begin(); v != Tri.finite_vertices_end(); ++v) {
@@ -54,7 +50,7 @@
}
CellHandle neighbourCell;
VertexHandle mirrorVertex;
- Vecteur tempVect;
+ CVector tempVect;
//FIXME : Ema, be carefull with this (noCache), it needs to be turned true after retriangulation
if (noCache) {
@@ -66,13 +62,13 @@
// if (!cell->info().Pcondition) ++nPCells;
// #endif
neighbourCell = cell->neighbor(j);
- const Vecteur& Surfk = cell->info().facetSurfaces[j];
+ const CVector& 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());
- Vecteur facetNormal = Surfk/area;
- const std::vector<Vecteur>& crossSections = cell->info().facetSphereCrossSections;
- Vecteur fluidSurfk = cell->info().facetSurfaces[j]*cell->info().facetFluidSurfacesRatio[j];
+ CVector facetNormal = Surfk/area;
+ const std::vector<CVector>& crossSections = cell->info().facetSphereCrossSections;
+ CVector 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[boundary(cell->vertex(j)->info().id()).coordinate]);
@@ -81,7 +77,7 @@
cell->info().unitForceVectors[j]=cell->info().unitForceVectors[j]+ tempVect;
}
/// Apply weighted forces f_k=sqRad_k/sumSqRad*f
- Vecteur facetUnitForce = -fluidSurfk*cell->info().solidSurfaces[j][3];
+ CVector facetUnitForce = -fluidSurfk*cell->info().solidSurfaces[j][3];
for (int y=0; y<3;y++) {
//add to cached value
cell->info().unitForceVectors[facetVertices[j][y]]=cell->info().unitForceVectors[facetVertices[j][y]]+facetUnitForce*cell->info().solidSurfaces[j][y];
@@ -114,7 +110,7 @@
}
}
if (debugOut) {
- Vecteur totalForce = nullVect;
+ CVector totalForce = nullVect;
for (FiniteVerticesIterator v = Tri.finite_vertices_begin(); v != Tri.finite_vertices_end(); v++)
{
if (!v->info().isFictious /*&& !v->info().isGhost*/ ){
@@ -133,13 +129,10 @@
VSolidTot = 0, Vtotalissimo = 0, vPoral = 0, sSolidTot = 0, vTotalePorosity=0, vPoralPorosity=0;
CellHandle neighbourCell;
- double k=0, distance = 0, radius = 0, viscosity = VISCOSITY;
- int surfneg=0;
- int NEG=0, POS=0, pass=0;
-
+ double k=0, distance = 0, radius = 0;
+ int surfneg=0; int NEG=0, POS=0, pass=0;
Real meanK=0, STDEV=0, meanRadius=0, meanDistance=0;
Real infiniteK=1e3;
-
double volume_sub_pore = 0.f;
VectorCell& cellHandles= T[currentTes].cellHandles;
for (VectorCell::iterator cellIt=T[currentTes].cellHandles.begin(); cellIt!=T[currentTes].cellHandles.end(); cellIt++){
@@ -163,13 +156,13 @@
for (int jj=0;jj<3;jj++)
cell->info().facetSphereCrossSections[j][jj]=0.5*W[jj]->point().weight()*Wm3::FastInvCos1((W[permut3[jj][1]]->point()-W[permut3[jj][0]]->point())*(W[permut3[jj][2]]->point()-W[permut3[jj][0]]->point()));
#else
- cell->info().facetSphereCrossSections[j]=Vecteur(
+ cell->info().facetSphereCrossSections[j]=CVector(
W[0]->info().isFictious ? 0 : 0.5*v0.weight()*acos((v1-v0)*(v2-v0)/sqrt((v1-v0).squared_length()*(v2-v0).squared_length())),
W[1]->info().isFictious ? 0 : 0.5*v1.weight()*acos((v0-v1)*(v2-v1)/sqrt((v1-v0).squared_length()*(v2-v1).squared_length())),
W[2]->info().isFictious ? 0 : 0.5*v2.weight()*acos((v0-v2)*(v1-v2)/sqrt((v1-v2).squared_length()*(v2-v0).squared_length())));
#endif
pass+=1;
- Vecteur l = p1 - p2;
+ CVector l = p1 - p2;
distance = sqrt(l.squared_length());
if (!RAVERAGE) radius = 2* computeHydraulicRadius(cell, j);
else radius = (computeEffectiveRadius(cell, j)+computeEquivalentRadius(cell,j))*0.5;
@@ -182,9 +175,9 @@
int test=0;
if (distance!=0) {
if (minPermLength>0 && distance_correction) distance=max(minPermLength,distance);
- const Vecteur& Surfk = cell->info().facetSurfaces[j];
+ const CVector& Surfk = cell->info().facetSurfaces[j];
Real area = sqrt(Surfk.squared_length());
- const Vecteur& crossSections = cell->info().facetSphereCrossSections[j];
+ const CVector& crossSections = cell->info().facetSphereCrossSections[j];
Real S0=0;
S0=checkSphereFacetOverlap(v0,v1,v2);
if (S0==0) S0=checkSphereFacetOverlap(v1,v2,v0);
@@ -200,12 +193,11 @@
if (k<0 && debugOut) {surfneg+=1;
cout<<"__ k<0 __"<<k<<" "<<" fluidArea "<<fluidArea<<" area "<<area<<" "<<crossSections[0]<<" "<<crossSections[1]<<" "<<crossSections[2] <<" "<<W[0]->info().id()<<" "<<W[1]->info().id()<<" "<<W[2]->info().id()<<" "<<p1<<" "<<p2<<" test "<<test<<endl;}
- } else {
- cout <<"infinite K2! surfaces will be missing (FIXME)"<<endl; k = infiniteK;
- }//Will be corrected in the next loop
+ } else cout <<"infinite K2! surfaces will be missing (FIXME)"<<endl; k = infiniteK;
+ //Will be corrected in the next loop
(cell->info().kNorm())[j]= k*kFactor;
if (!neighbourCell->info().isGhost) (neighbourCell->info().kNorm())[Tri.mirror_index(cell, j)]= (cell->info().kNorm())[j];
- //The following block is correct but very usefull, since all values are clamped below with MIN and MAX, skip for now
+ //The following block is correct but not very usefull, since all values are clamped below with MIN and MAX, skip for now
// else {//find the real neighbor connected to our cell through periodicity
// CellHandle true_neighbour_cell = cellHandles[neighbour_cell->info().baseIndex];
// for (int ii=0;ii<4;ii++)
@@ -305,8 +297,6 @@
}
if (debugOut) cout << "PassKopt = " << pass << endl;
}
-
-
if (debugOut) {
FiniteVerticesIterator verticesEnd = Tri.finite_vertices_end();
Real Vgrains = 0;
@@ -338,8 +328,6 @@
double compFlowFactor=0;
vector<Real> previousP;
previousP.resize(Tri.number_of_finite_cells());
- double tolerance = TOLERANCE;
- double relax = RELAX;
const int num_threads=1;
bool compressible= fluidBulkModulus>0;
=== modified file 'lib/triangulation/RegularTriangulation.h'
--- lib/triangulation/RegularTriangulation.h 2014-03-21 18:45:24 +0000
+++ lib/triangulation/RegularTriangulation.h 2014-03-21 18:45:24 +0000
@@ -5,13 +5,217 @@
* This program is free software; it is licensed under the terms of the *
* GNU General Public License v2 or later. See file LICENSE for details. *
*************************************************************************/
+
+//Define basic types from CGAL templates
#pragma once
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+#include <CGAL/Cartesian.h>
+#include <CGAL/Regular_triangulation_3.h>
+#include <CGAL/Regular_triangulation_euclidean_traits_3.h>
+#include <CGAL/Triangulation_vertex_base_with_info_3.h>
+#include <CGAL/Triangulation_cell_base_with_info_3.h>
+#include <CGAL/Delaunay_triangulation_3.h>
+#include <CGAL/circulator.h>
+#include <CGAL/number_utils.h>
+#include <boost/static_assert.hpp>
-#include "def_types.h"
-#include <cassert>
+//This include from yade let us use Eigen types
+#include <yade/lib/base/Math.hpp>
namespace CGT {
-
+//Robust kernel
+typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
+//A bit faster, but gives crash eventualy
+// typedef CGAL::Cartesian<double> K;
+
+typedef CGAL::Regular_triangulation_euclidean_traits_3<K> Traits;
+typedef K::Point_3 Point;
+typedef Traits::Vector_3 CVector;
+typedef Traits::Segment_3 Segment;
+#ifndef NO_REAL_CHECK
+/** compilation inside yade: check that Real in yade is the same as Real we will define; otherwise it might make things go wrong badly (perhaps) **/
+BOOST_STATIC_ASSERT(sizeof(Traits::RT)==sizeof(Real));
+#endif
+typedef Traits::RT Real; //Dans cartesian, RT = FT
+typedef Traits::Weighted_point Sphere;
+typedef Traits::Plane_3 Plane;
+typedef Traits::Triangle_3 Triangle;
+typedef Traits::Tetrahedron_3 Tetrahedron;
+
+class SimpleCellInfo : public Point {
+ public:
+ //"id": unique identifier of each cell, independant of other numberings used in the fluid types.
+ // Care to initialize it if you need it, there is no magic numbering to rely on
+ unsigned int id;
+ Real s;
+ bool isFictious;
+ SimpleCellInfo (void) {isFictious=false; s=0;}
+ SimpleCellInfo& operator= (const Point &p) { Point::operator= (p); return *this; }
+ SimpleCellInfo& operator= (const float &scalar) { s=scalar; return *this; }
+ inline Real x (void) {return Point::x();}
+ inline Real y (void) {return Point::y();}
+ inline Real z (void) {return Point::z();}
+ inline Real& f (void) {return s;}
+ //virtual function that will be defined for all classes, allowing shared function (e.g. for display of periodic and non-periodic with the unique function saveVTK)
+ virtual bool isReal (void) {return !isFictious;}
+};
+
+class SimpleVertexInfo : public CVector {
+protected:
+ Real s;
+ unsigned int i;
+ Real vol;
+public:
+ bool isFictious;
+ SimpleVertexInfo& operator= (const CVector &u) { CVector::operator= (u); return *this; }
+ SimpleVertexInfo& operator= (const float &scalar) { s=scalar; return *this; }
+ SimpleVertexInfo& operator= (const unsigned int &id) { i= id; return *this; }
+ inline Real ux (void) {return CVector::x();}
+ inline Real uy (void) {return CVector::y();}
+ inline Real uz (void) {return CVector::z();}
+ inline Real& f (void) {return s;}
+ inline Real& v (void) {return vol;}
+ inline const unsigned int& id (void) const {return i;}
+ SimpleVertexInfo (void) {isFictious=false; s=0; i=0; vol=-1;}
+ //virtual function that will be defined for all classes, allowing shared function (e.g. for display)
+ virtual bool isReal (void) {return !isFictious;}
+};
+
+class FlowCellInfo : public SimpleCellInfo {
+
+ public:
+ //For vector storage of all cells
+ unsigned int index;
+ int volumeSign;
+ bool Pcondition;
+ Real invVoidV;
+ Real t;
+ int fict;
+ Real volumeVariation;
+ double pression;
+ //average relative (fluid - facet translation) velocity defined for a single cell as 1/Volume * SUM_ON_FACETS(x_average_facet*average_facet_flow_rate)
+ CVector averageCellVelocity;
+ // Surface vectors of facets, pointing from outside toward inside the cell
+ std::vector<CVector> facetSurfaces;
+ //Ratio between fluid surface and facet surface
+ std::vector<Real> facetFluidSurfacesRatio;
+ // Reflects the geometrical property of the cell, so that the force by cell fluid on grain "i" is pressure*unitForceVectors[i]
+ std::vector<CVector> unitForceVectors;
+ // Store the area of triangle-sphere intersections for each facet (used in forces definition)
+ std::vector<CVector> facetSphereCrossSections;
+ std::vector<CVector> cellForce;
+ std::vector<double> rayHydr;
+ std::vector<double> modulePermeability;
+ // Partial surfaces of spheres in the double-tetrahedron linking two voronoi centers. [i][j] is for sphere facet "i" and sphere facetVertices[i][j]. Last component for 1/sum_surfaces in the facet.
+ double solidSurfaces [4][4];
+
+ FlowCellInfo (void)
+ {
+ modulePermeability.resize(4, 0);
+ cellForce.resize(4);
+ facetSurfaces.resize(4);
+ facetFluidSurfacesRatio.resize(4);
+ facetSphereCrossSections.resize(4);
+ unitForceVectors.resize(4);
+ for (int k=0; k<4;k++) for (int l=0; l<3;l++) solidSurfaces[k][l]=0;
+ rayHydr.resize(4, 0);
+// isInside = false;
+ invSumK=0;
+ isFictious=false; Pcondition = false; isGhost = false;
+// isInferior = false; isSuperior = false; isLateral = false; isExternal=false;
+ isvisited = false;
+ index=0;
+ volumeSign=0;
+ s=0;
+ volumeVariation=0;
+ pression=0;
+ invVoidV=0;
+ fict=0;
+ isGhost=false;
+ }
+ bool isGhost;
+ double invSumK;
+ bool isvisited;
+
+ FlowCellInfo& operator= (const std::vector<double> &v) { for (int i=0; i<4;i++) modulePermeability[i]= v[i]; return *this; }
+ FlowCellInfo& operator= (const Point &p) { Point::operator= (p); return *this; }
+ FlowCellInfo& operator= (const float &scalar) { s=scalar; return *this; }
+
+ inline Real& volume (void) {return t;}
+ inline const Real& invVoidVolume (void) const {return invVoidV;}
+ inline Real& invVoidVolume (void) {return invVoidV;}
+ inline Real& dv (void) {return volumeVariation;}
+ inline int& fictious (void) {return fict;}
+ inline double& p (void) {return pression;}
+ //For compatibility with the periodic case
+ inline const double shiftedP (void) const {return pression;}
+ inline const std::vector<double>& kNorm (void) const {return modulePermeability;}
+ inline std::vector<double>& kNorm (void) {return modulePermeability;}
+ inline std::vector< CVector >& facetSurf (void) {return facetSurfaces;}
+ inline std::vector<CVector>& force (void) {return cellForce;}
+ inline std::vector<double>& Rh (void) {return rayHydr;}
+ inline CVector& averageVelocity (void) {return averageCellVelocity;}
+};
+
+class FlowVertexInfo : public SimpleVertexInfo {
+ CVector grainVelocity;
+ Real volumeIncidentCells;
+public:
+ FlowVertexInfo& operator= (const CVector &u) { CVector::operator= (u); return *this; }
+ FlowVertexInfo& operator= (const float &scalar) { s=scalar; return *this; }
+ FlowVertexInfo& operator= (const unsigned int &id) { i= id; return *this; }
+ CVector forces;
+ bool isGhost;
+ FlowVertexInfo (void) {isGhost=false;}
+ inline CVector& force (void) {return forces;}
+ inline CVector& vel (void) {return grainVelocity;}
+ inline Real& volCells (void) {return volumeIncidentCells;}
+ inline const CVector ghostShift (void) {return CGAL::NULL_VECTOR;}
+};
+
+class PeriodicCellInfo : public FlowCellInfo
+{
+ public:
+ static CVector gradP;
+ //for real cell, baseIndex is the rank of the cell in cellHandles. For ghost cells, it is the baseIndex of the corresponding real cell.
+ //Unlike ordinary index, baseIndex is also indexing cells with imposed pressures
+ int baseIndex;
+ int period[3];
+ static CVector hSize[3];
+ static CVector deltaP;
+ int ghost;
+ Real* _pression;
+ PeriodicCellInfo (void){
+ _pression=&pression;
+ period[0]=period[1]=period[2]=0;
+ baseIndex=-1;
+ volumeSign=0;}
+ ~PeriodicCellInfo (void) {}
+ PeriodicCellInfo& operator= (const Point &p) { Point::operator= (p); return *this; }
+ PeriodicCellInfo& operator= (const float &scalar) { s=scalar; return *this; }
+
+ inline const double shiftedP (void) const {return isGhost? (*_pression)+pShift() :(*_pression) ;}
+ inline const double pShift (void) const {return deltaP[0]*period[0] + deltaP[1]*period[1] +deltaP[2]*period[2];}
+// inline const double p (void) {return shiftedP();}
+ inline void setP (const Real& p) {pression=p;}
+ virtual bool isReal (void) {return !(isFictious || isGhost);}
+};
+
+class PeriodicVertexInfo : public FlowVertexInfo {
+ public:
+ PeriodicVertexInfo& operator= (const CVector &u) { CVector::operator= (u); return *this; }
+ PeriodicVertexInfo& operator= (const float &scalar) { s=scalar; return *this; }
+ PeriodicVertexInfo& operator= (const unsigned int &id) { i= id; return *this; }
+ int period[3];
+ //FIXME: the name is misleading, even non-ghost can be out of the period and therefore they need to be shifted as well
+ inline const CVector ghostShift (void) {
+ return period[0]*PeriodicCellInfo::hSize[0]+period[1]*PeriodicCellInfo::hSize[1]+period[2]*PeriodicCellInfo::hSize[2];}
+ PeriodicVertexInfo (void) {isFictious=false; s=0; i=0; period[0]=period[1]=period[2]=0; isGhost=false;}
+ virtual bool isReal (void) {return !(isFictious || isGhost);}
+};
+
+
+///
template<class vertex_info, class cell_info>
class TriangulationTypes {
@@ -35,16 +239,15 @@
typedef typename RTriangulation::Cell_handle CellHandle;
typedef typename RTriangulation::Facet Facet;
-typedef typename RTriangulation::Facet_iterator FacetIterator;
+typedef typename RTriangulation::Facet_iterator FacetIterator;
typedef typename RTriangulation::Facet_circulator FacetCirculator;
-typedef typename RTriangulation::Finite_facets_iterator FiniteFacetsIterator;
+typedef typename RTriangulation::Finite_facets_iterator FiniteFacetsIterator;
typedef typename RTriangulation::Locate_type LocateType;
typedef typename RTriangulation::Edge_iterator EdgeIterator;
typedef typename RTriangulation::Finite_edges_iterator FiniteEdgesIterator;
};
-typedef CGAL::To_double<double> W_TO_DOUBLE; // foncteur Weight to Real
typedef TriangulationTypes<SimpleVertexInfo,SimpleCellInfo> SimpleTriangulationTypes;
typedef TriangulationTypes<FlowVertexInfo,FlowCellInfo> FlowTriangulationTypes;
typedef TriangulationTypes<PeriodicVertexInfo,PeriodicCellInfo> PeriFlowTriangulationTypes;
=== modified file 'lib/triangulation/Tenseur3.cpp'
--- lib/triangulation/Tenseur3.cpp 2014-03-21 18:45:24 +0000
+++ lib/triangulation/Tenseur3.cpp 2014-03-21 18:45:24 +0000
@@ -1,6 +1,6 @@
#include "Tenseur3.h"
-#include "def_types.h"
+#include "RegularTriangulation.h"
using namespace std;
namespace CGT {
@@ -14,16 +14,16 @@
return N;
}
-Vecteur operator* (Tens& tens, Vecteur& vect)
+CVector operator* (Tens& tens, CVector& vect)
{
- Vecteur result;
- result = Vecteur( tens(1,1)*vect.x()+ tens(1,2)*vect.y()+ tens(1,3)*vect.z(),
+ CVector result;
+ result = CVector( tens(1,1)*vect.x()+ tens(1,2)*vect.y()+ tens(1,3)*vect.z(),
tens(2,1)*vect.x()+ tens(2,2)*vect.y()+ tens(2,3)*vect.z(),
tens(3,1)*vect.x()+ tens(3,2)*vect.y()+ tens(3,3)*vect.z() );
return result;
}
-Vecteur& NormalizedVecteur (Vecteur& vect)
+CVector& NormalizedCVector (CVector& vect)
{
vect = vect*(1.0/sqrt(pow(vect.x(),2)+pow(vect.y(),2)+pow(vect.z(),2)));
return vect;
@@ -233,7 +233,7 @@
}
}
-void Tenseur_produit (Vecteur &v1, Vecteur &v2, Tenseur3 &result)
+void Tenseur_produit (CVector &v1, CVector &v2, Tenseur3 &result)
{
result(1,1) = v1.x()*v2.x();
result(1,2) = v1.x()*v2.y();
@@ -246,7 +246,7 @@
result(3,3) = v1.z()*v2.z();
}
-void Somme (Tenseur3 &result, Vecteur &v1, Vecteur &v2)
+void Somme (Tenseur3 &result, CVector &v1, CVector &v2)
{
result(1,1) += v1.x()*v2.x();
result(1,2) += v1.x()*v2.y();
=== modified file 'lib/triangulation/Tenseur3.h'
--- lib/triangulation/Tenseur3.h 2013-08-22 14:32:01 +0000
+++ lib/triangulation/Tenseur3.h 2014-03-21 18:45:24 +0000
@@ -1,8 +1,8 @@
#pragma once
-#include "def_types.h"
#include <iostream>
#include <fstream>
+#include "RegularTriangulation.h"
namespace CGT {
@@ -13,12 +13,12 @@
class Tenseur_sym3;
class Tenseur_anti3;
-Vecteur operator* ( Tens& tens, Vecteur& vect );
-Vecteur& NormalizedVecteur ( Vecteur& vect );
-
-
-void Tenseur_produit ( Vecteur &v1, Vecteur &v2, Tenseur3 &result );
-void Somme ( Tenseur3 &result, Vecteur &v1, Vecteur &v2 );
+CVector operator* ( Tens& tens, CVector& vect );
+CVector& NormalizedCVector ( CVector& vect );
+
+
+void Tenseur_produit ( CVector &v1, CVector &v2, Tenseur3 &result );
+void Somme ( Tenseur3 &result, CVector &v1, CVector &v2 );
std::ostream& operator<< ( std::ostream& os,const Tenseur3& T );
std::ostream& operator<< ( std::ostream& os,const Tenseur_sym3& T );
=== modified file 'lib/triangulation/Tesselation.h'
--- lib/triangulation/Tesselation.h 2014-03-21 18:45:24 +0000
+++ lib/triangulation/Tesselation.h 2014-03-21 18:45:24 +0000
@@ -1,5 +1,5 @@
/*************************************************************************
-* Copyright (C) 2006 by Bruno Chareyre *
+* Copyright (C) 2006 by Bruno Chareyre *
* bruno.chareyre@xxxxxxxxxxx *
* *
* This program is free software; it is licensed under the terms of the *
=== modified file 'lib/triangulation/TriaxialState.cpp'
--- lib/triangulation/TriaxialState.cpp 2012-07-10 21:05:26 +0000
+++ lib/triangulation/TriaxialState.cpp 2014-03-21 18:45:24 +0000
@@ -144,7 +144,7 @@
z <= (box.sommet.z()-filter_distance*mean_radius) );
}
-bool TriaxialState::inside(Vecteur v)
+bool TriaxialState::inside(CVector v)
{
return TriaxialState::inside(v.x(), v.y(), v.z());
}
@@ -182,7 +182,7 @@
//Real tx, ty, tz;
Point pos(CGAL::ORIGIN);
mean_radius=0;
- Vecteur trans, rot;
+ CVector trans, rot;
Real rad; //coordonn�es/rayon
bool isSphere;
@@ -221,7 +221,7 @@
long id1, id2;
int stat;
- Vecteur c_pos, normal, old_fs, fs;
+ CVector c_pos, normal, old_fs, fs;
Real old_fn, fn, frictional_work;
Statefile >> Nc;
contacts.resize(Nc);
=== modified file 'lib/triangulation/TriaxialState.h'
--- lib/triangulation/TriaxialState.h 2013-08-22 14:32:01 +0000
+++ lib/triangulation/TriaxialState.h 2014-03-21 18:45:24 +0000
@@ -41,8 +41,8 @@
int id;
bool isSphere;
Sphere sphere;
- Vecteur translation;
- Vecteur rotation;
+ CVector translation;
+ CVector rotation;
VectorContact contacts;
Grain(void) {id=-1; isSphere=true;}
@@ -52,12 +52,12 @@
Grain* grain1;
Grain* grain2;
- Vecteur position;
- Vecteur normal;
+ CVector position;
+ CVector normal;
Real fn;
- Vecteur fs;
+ CVector fs;
Real old_fn;
- Vecteur old_fs;
+ CVector old_fs;
Real frictional_work;
bool visited;
Status status;
@@ -70,7 +70,7 @@
bool from_file(const char* filename, bool bz2=false);
bool to_file(const char* filename, bool bz2=false);
bool inside(Real x, Real y, Real z);
- bool inside(Vecteur v);
+ bool inside(CVector v);
bool inside(Point p);
static Real find_parameter (const char* parameter_name, const char* filename);
static Real find_parameter (const char* parameter_name, boost::iostreams::filtering_istream& file);
=== removed file 'lib/triangulation/def_types.h'
--- lib/triangulation/def_types.h 2014-03-21 18:45:24 +0000
+++ lib/triangulation/def_types.h 1970-01-01 00:00:00 +0000
@@ -1,220 +0,0 @@
-/*************************************************************************
-* Copyright (C) 2006 by Bruno Chareyre, bruno.chareyre@xxxxxxxxxxx *
-* Copyright (C) 2010 by Emanuele Catalano, emanuele.catalano@xxxxxxxxxxx*
-* This program is free software; it is licensed under the terms of the *
-* GNU General Public License v2 or later. See file LICENSE for details. *
-*************************************************************************/
-
-#include <yade/lib/base/Math.hpp>
-
-//Define basic types from CGAL templates
-#pragma once
-
-#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
-#include <CGAL/Cartesian.h>
-#include <CGAL/Regular_triangulation_3.h>
-#include <CGAL/Regular_triangulation_euclidean_traits_3.h>
-#include <CGAL/Triangulation_vertex_base_with_info_3.h>
-#include <CGAL/Triangulation_cell_base_with_info_3.h>
-#include <CGAL/Delaunay_triangulation_3.h>
-#include <CGAL/circulator.h>
-#include <CGAL/number_utils.h>
-#include <boost/static_assert.hpp>
-
-//This include from yade let us use Eigen types
-#include <yade/lib/base/Math.hpp>
-
-namespace CGT {
-//Robust kernel
-typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
-//A bit faster, but gives crash eventualy
-// typedef CGAL::Cartesian<double> K;
-
-typedef CGAL::Regular_triangulation_euclidean_traits_3<K> Traits;
-typedef K::Point_3 Point;
-typedef Traits::Vector_3 Vecteur;
-typedef Traits::Segment_3 Segment;
-#ifndef NO_REAL_CHECK
-/** compilation inside yade: check that Real in yade is the same as Real we will define; otherwise it might make things go wrong badly (perhaps) **/
-BOOST_STATIC_ASSERT(sizeof(Traits::RT)==sizeof(Real));
-#endif
-typedef Traits::RT Real; //Dans cartesian, RT = FT
-typedef Traits::Weighted_point Sphere;
-typedef Traits::Plane_3 Plane;
-typedef Traits::Triangle_3 Triangle;
-typedef Traits::Tetrahedron_3 Tetrahedron;
-
-class SimpleCellInfo : public Point {
- public:
- //"id": unique identifier of each cell, independant of other numberings used in the fluid types.
- // Care to initialize it if you need it, there is no magic numbering to rely on
- unsigned int id;
- Real s;
- bool isFictious;
- SimpleCellInfo (void) {isFictious=false; s=0;}
- SimpleCellInfo& operator= (const Point &p) { Point::operator= (p); return *this; }
- SimpleCellInfo& operator= (const float &scalar) { s=scalar; return *this; }
- inline Real x (void) {return Point::x();}
- inline Real y (void) {return Point::y();}
- inline Real z (void) {return Point::z();}
- inline Real& f (void) {return s;}
- //virtual function that will be defined for all classes, allowing shared function (e.g. for display of periodic and non-periodic with the unique function saveVTK)
- virtual bool isReal (void) {return !isFictious;}
-};
-
-class SimpleVertexInfo : public Vecteur {
-protected:
- Real s;
- unsigned int i;
- Real vol;
-public:
- bool isFictious;
- SimpleVertexInfo& operator= (const Vecteur &u) { Vecteur::operator= (u); return *this; }
- SimpleVertexInfo& operator= (const float &scalar) { s=scalar; return *this; }
- SimpleVertexInfo& operator= (const unsigned int &id) { i= id; return *this; }
- inline Real ux (void) {return Vecteur::x();}
- inline Real uy (void) {return Vecteur::y();}
- inline Real uz (void) {return Vecteur::z();}
- inline Real& f (void) {return s;}
- inline Real& v (void) {return vol;}
- inline const unsigned int& id (void) const {return i;}
- SimpleVertexInfo (void) {isFictious=false; s=0; i=0; vol=-1;}
- //virtual function that will be defined for all classes, allowing shared function (e.g. for display)
- virtual bool isReal (void) {return !isFictious;}
-};
-
-class FlowCellInfo : public SimpleCellInfo {
-
- public:
- //For vector storage of all cells
- unsigned int index;
- int volumeSign;
- bool Pcondition;
- Real invVoidV;
- Real t;
- int fict;
- Real volumeVariation;
- double pression;
- //average relative (fluid - facet translation) velocity defined for a single cell as 1/Volume * SUM_ON_FACETS(x_average_facet*average_facet_flow_rate)
- Vecteur averageCellVelocity;
- // Surface vectors of facets, pointing from outside toward inside the cell
- std::vector<Vecteur> facetSurfaces;
- //Ratio between fluid surface and facet surface
- std::vector<Real> facetFluidSurfacesRatio;
- // Reflects the geometrical property of the cell, so that the force by cell fluid on grain "i" is pressure*unitForceVectors[i]
- std::vector<Vecteur> unitForceVectors;
- // Store the area of triangle-sphere intersections for each facet (used in forces definition)
- std::vector<Vecteur> facetSphereCrossSections;
- std::vector<Vecteur> cellForce;
- std::vector<double> rayHydr;
- std::vector<double> modulePermeability;
- // Partial surfaces of spheres in the double-tetrahedron linking two voronoi centers. [i][j] is for sphere facet "i" and sphere facetVertices[i][j]. Last component for 1/sum_surfaces in the facet.
- double solidSurfaces [4][4];
-
- FlowCellInfo (void)
- {
- modulePermeability.resize(4, 0);
- cellForce.resize(4);
- facetSurfaces.resize(4);
- facetFluidSurfacesRatio.resize(4);
- facetSphereCrossSections.resize(4);
- unitForceVectors.resize(4);
- for (int k=0; k<4;k++) for (int l=0; l<3;l++) solidSurfaces[k][l]=0;
- rayHydr.resize(4, 0);
-// isInside = false;
- invSumK=0;
- isFictious=false; Pcondition = false; isGhost = false;
-// isInferior = false; isSuperior = false; isLateral = false; isExternal=false;
- isvisited = false;
- index=0;
- volumeSign=0;
- s=0;
- volumeVariation=0;
- pression=0;
- invVoidV=0;
- fict=0;
- isGhost=false;
- }
- bool isGhost;
- double invSumK;
- bool isvisited;
-
- FlowCellInfo& operator= (const std::vector<double> &v) { for (int i=0; i<4;i++) modulePermeability[i]= v[i]; return *this; }
- FlowCellInfo& operator= (const Point &p) { Point::operator= (p); return *this; }
- FlowCellInfo& operator= (const float &scalar) { s=scalar; return *this; }
-
- inline Real& volume (void) {return t;}
- inline const Real& invVoidVolume (void) const {return invVoidV;}
- inline Real& invVoidVolume (void) {return invVoidV;}
- inline Real& dv (void) {return volumeVariation;}
- inline int& fictious (void) {return fict;}
- inline double& p (void) {return pression;}
- //For compatibility with the periodic case
- inline const double shiftedP (void) const {return pression;}
- inline const std::vector<double>& kNorm (void) const {return modulePermeability;}
- inline std::vector<double>& kNorm (void) {return modulePermeability;}
- inline std::vector< Vecteur >& facetSurf (void) {return facetSurfaces;}
- inline std::vector<Vecteur>& force (void) {return cellForce;}
- inline std::vector<double>& Rh (void) {return rayHydr;}
- inline Vecteur& averageVelocity (void) {return averageCellVelocity;}
-};
-
-class FlowVertexInfo : public SimpleVertexInfo {
- Vecteur grainVelocity;
- Real volumeIncidentCells;
-public:
- FlowVertexInfo& operator= (const Vecteur &u) { Vecteur::operator= (u); return *this; }
- FlowVertexInfo& operator= (const float &scalar) { s=scalar; return *this; }
- FlowVertexInfo& operator= (const unsigned int &id) { i= id; return *this; }
- Vecteur forces;
- bool isGhost;
- FlowVertexInfo (void) {isGhost=false;}
- inline Vecteur& force (void) {return forces;}
- inline Vecteur& vel (void) {return grainVelocity;}
- inline Real& volCells (void) {return volumeIncidentCells;}
- inline const Vecteur ghostShift (void) {return CGAL::NULL_VECTOR;}
-};
-
-class PeriodicCellInfo : public FlowCellInfo
-{
- public:
- static Vecteur gradP;
- //for real cell, baseIndex is the rank of the cell in cellHandles. For ghost cells, it is the baseIndex of the corresponding real cell.
- //Unlike ordinary index, baseIndex is also indexing cells with imposed pressures
- int baseIndex;
- int period[3];
- static Vecteur hSize[3];
- static Vecteur deltaP;
- int ghost;
- Real* _pression;
- PeriodicCellInfo (void){
- _pression=&pression;
- period[0]=period[1]=period[2]=0;
- baseIndex=-1;
- volumeSign=0;}
- ~PeriodicCellInfo (void) {}
- PeriodicCellInfo& operator= (const Point &p) { Point::operator= (p); return *this; }
- PeriodicCellInfo& operator= (const float &scalar) { s=scalar; return *this; }
-
- inline const double shiftedP (void) const {return isGhost? (*_pression)+pShift() :(*_pression) ;}
- inline const double pShift (void) const {return deltaP[0]*period[0] + deltaP[1]*period[1] +deltaP[2]*period[2];}
-// inline const double p (void) {return shiftedP();}
- inline void setP (const Real& p) {pression=p;}
- virtual bool isReal (void) {return !(isFictious || isGhost);}
-};
-
-class PeriodicVertexInfo : public FlowVertexInfo {
- public:
- PeriodicVertexInfo& operator= (const Vecteur &u) { Vecteur::operator= (u); return *this; }
- PeriodicVertexInfo& operator= (const float &scalar) { s=scalar; return *this; }
- PeriodicVertexInfo& operator= (const unsigned int &id) { i= id; return *this; }
- int period[3];
- //FIXME: the name is misleading, even non-ghost can be out of the period and therefore they need to be shifted as well
- inline const Vecteur ghostShift (void) {
- return period[0]*PeriodicCellInfo::hSize[0]+period[1]*PeriodicCellInfo::hSize[1]+period[2]*PeriodicCellInfo::hSize[2];}
- PeriodicVertexInfo (void) {isFictious=false; s=0; i=0; period[0]=period[1]=period[2]=0; isGhost=false;}
- virtual bool isReal (void) {return !(isFictious || isGhost);}
-};
-
-
-} // namespace CGT
=== modified file 'pkg/dem/FlowEngine.cpp'
--- pkg/dem/FlowEngine.cpp 2014-03-21 18:45:24 +0000
+++ pkg/dem/FlowEngine.cpp 2014-03-21 18:45:24 +0000
@@ -30,10 +30,9 @@
CREATE_LOGGER ( FlowEngine );
CREATE_LOGGER ( PeriodicFlowEngine );
-CGT::Vecteur makeCgVect ( const Vector3r& yv ) {return CGT::Vecteur ( yv[0],yv[1],yv[2] );}
-CGT::Point makeCgPoint ( const Vector3r& yv ) {return CGT::Point ( yv[0],yv[1],yv[2] );}
+CGT::CVector makeCgVect ( const Vector3r& yv ) {return CGT::CVector ( yv[0],yv[1],yv[2] );}
Vector3r makeVector3r ( const CGT::Point& yv ) {return Vector3r ( yv[0],yv[1],yv[2] );}
-Vector3r makeVector3r ( const CGT::Vecteur& yv ) {return Vector3r ( yv[0],yv[1],yv[2] );}
+Vector3r makeVector3r ( const CGT::CVector& yv ) {return Vector3r ( yv[0],yv[1],yv[2] );}
FlowEngine::~FlowEngine()
@@ -213,9 +212,9 @@
flow->numFactorizeThreads = numFactorizeThreads;
#endif
flow->meanKStat = meanKStat;
- flow->VISCOSITY = viscosity;
- flow->TOLERANCE=tolerance;
- flow->RELAX=relax;
+ flow->viscosity = viscosity;
+ flow->tolerance=tolerance;
+ flow->relax=relax;
flow->clampKValues = clampKValues;
flow->maxKdivKmean = maxKdivKmean;
flow->minKdivKmean = minKdivKmean;
@@ -353,7 +352,7 @@
double center[3];
for ( int i=0; i<6; i++ ) {
if ( *flow->boundsIds[i]<0 ) continue;
- CGT::Vecteur Normal ( normal[i].x(), normal[i].y(), normal[i].z() );
+ CGT::CVector Normal ( normal[i].x(), normal[i].y(), normal[i].z() );
if ( flow->boundary ( *flow->boundsIds[i] ).useMaxMin ) flow->addBoundingPlane(Normal, *flow->boundsIds[i] );
else {
for ( int h=0;h<3;h++ ) center[h] = buffer[*flow->boundsIds[i]].pos[h];
@@ -397,7 +396,7 @@
typedef typename Solver::element_type Flow;
FiniteVerticesIterator verticesEnd = flow->T[flow->currentTes].Triangulation().finite_vertices_end();
- CGT::Vecteur Zero(0,0,0);
+ CGT::CVector Zero(0,0,0);
for (FiniteVerticesIterator vIt = flow->T[flow->currentTes].Triangulation().finite_vertices_begin(); vIt!= verticesEnd; vIt++) vIt->info().forces=Zero;
FOREACH(CellHandle& cell, flow->T[flow->currentTes].cellHandles)
@@ -408,9 +407,9 @@
case ( 1 ) : cell->info().volume() = volumeCellSingleFictious ( cell ); break;
case ( 2 ) : cell->info().volume() = volumeCellDoubleFictious ( cell ); break;
case ( 3 ) : cell->info().volume() = volumeCellTripleFictious ( cell ); break;
- default: break;
+ default: break;
}
- if (flow->fluidBulkModulus>0) { cell->info().invVoidVolume() = 1 / ( abs(cell->info().volume()) - flow->volumeSolidPore(cell) ); }
+ if (flow->fluidBulkModulus>0) { cell->info().invVoidVolume() = 1. / ( abs(cell->info().volume()) - flow->volumeSolidPore(cell) ); }
}
if (debug) cout << "Volumes initialised." << endl;
}
@@ -431,14 +430,14 @@
RTriangulation& Tri = solver->T[solver->currentTes].Triangulation();
CGT::Point posAvFacet;
double volumeFacetTranslation = 0;
- CGT::Vecteur velAv ( Vel[0], Vel[1], Vel[2] );
+ CGT::CVector velAv ( Vel[0], Vel[1], Vel[2] );
for ( int i=0; i<4; i++ ) {
volumeFacetTranslation = 0;
if ( !Tri.is_infinite ( cell->neighbor ( i ) ) ) {
- CGT::Vecteur Surfk = cell->info()-cell->neighbor ( i )->info();
+ CGT::CVector Surfk = cell->info()-cell->neighbor ( i )->info();
Real area = sqrt ( Surfk.squared_length() );
Surfk = Surfk/area;
- CGT::Vecteur branch = cell->vertex ( facetVertices[i][0] )->point() - cell->info();
+ CGT::CVector branch = cell->vertex ( facetVertices[i][0] )->point() - cell->info();
posAvFacet = ( CGT::Point ) cell->info() + ( branch*Surfk ) *Surfk;
volumeFacetTranslation += velAv*cell->info().facetSurfaces[i];
cell->info().averageVelocity() = cell->info().averageVelocity() - volumeFacetTranslation/cell->info().volume() * ( posAvFacet-CGAL::ORIGIN );
@@ -1077,7 +1076,7 @@
void PeriodicFlowEngine::initializeVolumes (shared_ptr<FlowSolver>& flow)
{
FiniteVerticesIterator verticesEnd = flow->T[flow->currentTes].Triangulation().finite_vertices_end();
- CGT::Vecteur Zero ( 0,0,0 );
+ CGT::CVector Zero ( 0,0,0 );
for ( FiniteVerticesIterator vIt = flow->T[flow->currentTes].Triangulation().finite_vertices_begin(); vIt!= verticesEnd; vIt++ ) vIt->info().forces=Zero;
FOREACH(CellHandle& cell, flow->T[flow->currentTes].cellHandles){
@@ -1131,7 +1130,7 @@
if ( debug ) cout << endl << "locateCell------" << endl << endl;
flow->computePermeability ( );
porosity = flow->vPoralPorosity/flow->vTotalePorosity;
- flow->TOLERANCE=tolerance;flow->RELAX=relax;
+ flow->tolerance=tolerance;flow->relax=relax;
flow->displayStatistics ();
//FIXME: check interpolate() for the periodic case, at least use the mean pressure from previous step.
@@ -1150,7 +1149,7 @@
CGT::PeriodicCellInfo::hSize[0] = makeCgVect ( scene->cell->hSize.col ( 0 ) );
CGT::PeriodicCellInfo::hSize[1] = makeCgVect ( scene->cell->hSize.col ( 1 ) );
CGT::PeriodicCellInfo::hSize[2] = makeCgVect ( scene->cell->hSize.col ( 2 ) );
- CGT::PeriodicCellInfo::deltaP=CGT::Vecteur (
+ CGT::PeriodicCellInfo::deltaP=CGT::CVector (
CGT::PeriodicCellInfo::hSize[0]*CGT::PeriodicCellInfo::gradP,
CGT::PeriodicCellInfo::hSize[1]*CGT::PeriodicCellInfo::gradP,
CGT::PeriodicCellInfo::hSize[2]*CGT::PeriodicCellInfo::gradP );
=== modified file 'pkg/dem/FlowEngine.hpp'
--- pkg/dem/FlowEngine.hpp 2014-03-21 18:45:24 +0000
+++ pkg/dem/FlowEngine.hpp 2014-03-21 18:45:24 +0000
@@ -6,6 +6,30 @@
* This program is free software; it is licensed under the terms of the *
* GNU General Public License v2 or later. See file LICENSE for details. *
*************************************************************************/
+
+/*!
+
+FlowEngine is mainly an interface to the fluid solvers defined in lib/triangulation and using the PFV scheme.
+There are also a few non-trivial functions defined here, such has building triangulation and computating of elements volumes.
+
+PeriodicFlowEngine is a variant for periodic boundary conditions.
+
+Which solver will be actually used internally to obtain pore pressure will depend on compile time flags (libcholmod available and #define LINSOLV),
+and on runtime settings (useSolver=0: Gauss-Seidel (iterative), useSolver=1: CHOLMOD if available (direct sparse cholesky)).
+
+The files defining lower level classes are in yade/lib. The code uses CGAL::Triangulation3 for managing the mesh and storing data; eigen3::Sparse, suitesparse::cholmod, and metis for direct resolution of the linear systems. The Gauss-Seidel iterative solver OTOH is implemented directly in Yade.
+
+Most classes in lib/triangulation are templates, and therefore completely defined in header files.
+A pseudo hpp/cpp split is reproduced for clarity, with hpp/ipp extensions (but again, in the end they are all included).
+The files are
+- RegularTriangulation.h : declaration of info types embeded in the mesh (template parameters of CGAL::Triangulation3), and instanciation of RTriangulation classes
+- Tesselation.h/ipp : encapsulate RegularTriangulation and adds functions to manipulate the dual (Voronoi) graph of the triangulation
+- Network.hpp/ipp : specialized for PFV model (the two former are used independently by TesselationWrapper), a set of functions to determine volumes and surfaces of intersections between spheres and various subdomains. Contains two triangulations for smooth transitions while remeshing - e.g. interpolating values in the new mesh using the previous one.
+- FlowBoundingSphere.hpp/ipp and PeriodicFlow.hpp/ipp + LinSolv variants: implement the solver in itself (mesh, boundary conditions, solving, defining fluid-particles interactions)
+- FlowEngine.hpp/cpp (this file)
+
+*/
+
#pragma once
#include<yade/core/PartialEngine.hpp>
#include<yade/pkg/dem/TriaxialCompressionEngine.hpp>
@@ -84,7 +108,7 @@
TPL Real getCellFlux(unsigned int cond, const shared_ptr<Solver>& flow);
TPL Real getBoundaryFlux(unsigned int boundary,Solver& flow) {return flow->boundaryFlux(boundary);}
TPL Vector3r fluidForce(unsigned int id_sph, Solver& flow) {
- const CGT::Vecteur& f=flow->T[flow->currentTes].vertex(id_sph)->info().forces; return Vector3r(f[0],f[1],f[2]);}
+ const CGT::CVector& f=flow->T[flow->currentTes].vertex(id_sph)->info().forces; return Vector3r(f[0],f[1],f[2]);}
TPL Vector3r shearLubForce(unsigned int id_sph, Solver& flow) {
return (flow->shearLubricationForces.size()>id_sph)?flow->shearLubricationForces[id_sph]:Vector3r::Zero();}
TPL Vector3r shearLubTorque(unsigned int id_sph, Solver& flow) {
@@ -157,10 +181,10 @@
double averageSlicePressure(double posY){return solver->averageSlicePressure(posY);}
double averagePressure(){return solver->averagePressure();}
#ifdef LINSOLV
- TPL void exportMatrix(string filename,Solver& flow) {if (useSolver==3) flow->exportMatrix(filename.c_str());
- else cerr<<"available for Cholmod solver (useSolver==3)"<<endl;}
- TPL void exportTriplets(string filename,Solver& flow) {if (useSolver==3) flow->exportTriplets(filename.c_str());
- else cerr<<"available for Cholmod solver (useSolver==3)"<<endl;}
+ TPL void exportMatrix(string filename,Solver& flow) {if (useSolver>0) flow->exportMatrix(filename.c_str());
+ else LOG_ERROR("available for Cholmod solver (useSolver>0)");}
+ TPL void exportTriplets(string filename,Solver& flow) {if (useSolver>0) flow->exportTriplets(filename.c_str());
+ else LOG_ERROR("available for Cholmod solver (useSolver>0)");}
#endif
void emulateAction(){
@@ -239,7 +263,7 @@
((double,permeabilityFactor,1.0,,"permability multiplier"))
((double,viscosity,1.0,,"viscosity of the fluid"))
((double,stiffness, 10000,,"equivalent contact stiffness used in the lubrication model"))
- ((int, useSolver, 0,, "Solver to use 0=G-Seidel, 1=Taucs, 2-Pardiso, 3-CHOLMOD"))
+ ((int, useSolver, 0,, "Solver to use: 0=Gauss-Seidel (iterative), 1=CHOLMOD (direct sparse cholesky)"))
((int, xmin,0,(Attr::readonly),"Index of the boundary $x_{min}$. This index is not equal the the id of the corresponding body in general, it may be used to access the corresponding attributes (e.g. flow.bndCondValue[flow.xmin], flow.wallId[flow.xmin],...)."))
((int, xmax,1,(Attr::readonly),"See :yref:`FlowEngine::xmin`."))
((int, ymin,2,(Attr::readonly),"See :yref:`FlowEngine::xmin`."))
=== modified file 'pkg/dem/MicroMacroAnalyser.cpp'
--- pkg/dem/MicroMacroAnalyser.cpp 2012-07-10 21:05:26 +0000
+++ pkg/dem/MicroMacroAnalyser.cpp 2014-03-21 18:45:24 +0000
@@ -130,7 +130,7 @@
// TS.grains[Idg].translation = trans;
AngleAxisr aa((*bi)->state->ori);
Vector3r rotVec=aa.axis()*aa.angle();
- TS.grains[Idg].rotation = CGT::Vecteur(rotVec[0],rotVec[1],rotVec[2]);
+ TS.grains[Idg].rotation = CGT::CVector(rotVec[0],rotVec[1],rotVec[2]);
TS.box.base = CGT::Point(min(TS.box.base.x(), pos.x()-rad),
min(TS.box.base.y(), pos.y()-rad),
min(TS.box.base.z(), pos.z()-rad));
@@ -173,12 +173,12 @@
c->grain2 = & (TS.grains[id2]);
grains[id1].contacts.push_back(c);
grains[id2].contacts.push_back(c);
- c->normal = CGT::Vecteur((YADE_CAST<ScGeom*> ((*ii)->geom.get()))->normal.x(),
+ c->normal = CGT::CVector((YADE_CAST<ScGeom*> ((*ii)->geom.get()))->normal.x(),
(YADE_CAST<ScGeom*> ((*ii)->geom.get()))->normal.y(),
(YADE_CAST<ScGeom*> ((*ii)->geom.get()))->normal.z());
// c->normal = ( grains[id2].sphere.point()-grains[id1].sphere.point() );
// c->normal = c->normal/sqrt ( pow ( c->normal.x(),2 ) +pow ( c->normal.y(),2 ) +pow ( c->normal.z(),2 ) );
- c->position = CGT::Vecteur((YADE_CAST<ScGeom*> ((*ii)->geom.get()))->contactPoint.x(),
+ c->position = CGT::CVector((YADE_CAST<ScGeom*> ((*ii)->geom.get()))->contactPoint.x(),
(YADE_CAST<ScGeom*> ((*ii)->geom.get()))->contactPoint.y(),
(YADE_CAST<ScGeom*> ((*ii)->geom.get()))->contactPoint.z());
// c->position = 0.5* ( ( grains[id1].sphere.point()-CGAL::ORIGIN ) +
@@ -187,7 +187,7 @@
// ( grains[id2].sphere.weight() *c->normal ) );
c->fn = YADE_CAST<FrictPhys*> (((*ii)->phys.get()))->normalForce.dot((YADE_CAST<ScGeom*> ((*ii)->geom.get()))->normal);
Vector3r fs = YADE_CAST<FrictPhys*> ((*ii)->phys.get())->shearForce;
- c->fs = CGT::Vecteur(fs.x(),fs.y(),fs.z());
+ c->fs = CGT::CVector(fs.x(),fs.y(),fs.z());
c->old_fn = c->fn;
c->old_fs = c->fs;
c->frictional_work = 0;
