Commit 61ec19ab authored by Bruno Chareyre's avatar Bruno Chareyre
Browse files

some renaming and code cleaning in PFV code

parent 0ad481e8
......@@ -29,15 +29,14 @@ class FlowBoundingSphere : public Network<_Tesselation>
DECLARE_TESSELATION_TYPES(Network<Tesselation>)
//painfull, but we need that for templates inheritance...
using _N::T; using _N::x_min; using _N::x_max; using _N::y_min; using _N::y_max; using _N::z_min; using _N::z_max; using _N::Rmoy; using _N::SectionArea; using _N::Height; using _N::Vtotale; using _N::currentTes; using _N::DEBUG_OUT; using _N::nOfSpheres; using _N::x_min_id; using _N::x_max_id; using _N::y_min_id; using _N::y_max_id; using _N::z_min_id; using _N::z_max_id; using _N::boundsIds; using _N::Corner_min; using _N::Corner_max; using _N::Vsolid_tot; using _N::Vtotalissimo; using _N::Vporale; using _N::Ssolid_tot; using _N::V_porale_porosity; using _N::V_totale_porosity; using _N::boundaries; using _N::id_offset; using _N::vtk_infinite_vertices; using _N::vtk_infinite_cells; using _N::num_particles; using _N::boundingCells; using _N::facetVertices; using _N::facetNFictious;
using _N::T; using _N::xMin; using _N::xMax; using _N::yMin; using _N::yMax; using _N::zMin; using _N::zMax; using _N::Rmoy; using _N::sectionArea; using _N::Height; using _N::vTotal; using _N::currentTes; using _N::debugOut; using _N::nOfSpheres; using _N::xMinId; using _N::xMaxId; using _N::yMinId; using _N::yMaxId; using _N::zMinId; using _N::zMaxId; using _N::boundsIds; using _N::cornerMin; using _N::cornerMax; using _N::VSolidTot; using _N::Vtotalissimo; using _N::vPoral; using _N::sSolidTot; using _N::vPoralPorosity; using _N::vTotalePorosity; using _N::boundaries; using _N::idOffset; using _N::vtkInfiniteVertices; using _N::vtkInfiniteCells; using _N::num_particles; using _N::boundingCells; using _N::facetVertices; using _N::facetNFictious;
//same for functions
using _N::Define_fictious_cells; using _N::AddBoundingPlanes; using _N::boundary;
using _N::defineFictiousCells; using _N::addBoundingPlanes; using _N::boundary;
virtual ~FlowBoundingSphere();
FlowBoundingSphere();
bool SLIP_ON_LATERALS;
// bool areaR2Permeability;
bool slipOnLaterals;
double TOLERANCE;
double RELAX;
double ks; //Hydraulic Conductivity
......@@ -48,17 +47,17 @@ class FlowBoundingSphere : public Network<_Tesselation>
bool pressureChanged;//are imposed pressures modified (on python side)? When it happens, we have to reApplyBoundaryConditions
int errorCode;
//Handling imposed pressures/fluxes on elements in the form of {point,value} pairs, IPCells contains the cell handles corresponding to point
//Handling imposed pressures/fluxes on elements in the form of {point,value} pairs, ipCells contains the cell handles corresponding to point
vector<pair<Point,Real> > imposedP;
vector<Cell_handle> IPCells;
vector<CellHandle> ipCells;
vector<pair<Point,Real> > imposedF;
vector<Cell_handle> IFCells;
vector<CellHandle> ifCells;
void initNewTri () {noCache=true; /*isLinearSystemSet=false; areCellsOrdered=false;*/}//set flags after retriangulation
bool permeability_map;
bool permeabilityMap;
bool computeAllCells;//exececute computeHydraulicRadius for all facets and all spheres (double cpu time but needed for now in order to define crossSections correctly)
double K_opt_factor;
double KOptFactor;
double minKdivKmean;
double maxKdivKmean;
int Iterations;
......@@ -71,15 +70,16 @@ class FlowBoundingSphere : public Network<_Tesselation>
vector< vector<const Vecteur*> > perVertexUnitForce;
vector< vector<const Real*> > perVertexPressure;
#endif
vector <Finite_edges_iterator> Edge_list;
vector <double> Edge_Surfaces;
vector <pair<int,int> > Edge_ids;
vector <double> edgeSurfaces;
vector <pair<int,int> > edgeIds;
vector <Real> edgeNormalLubF;
vector <Vector3r> viscousShearForces;
vector <Vector3r> viscousShearTorques;
vector <Vector3r> normLubForce;
vector <Matrix3r> viscousBodyStress;
vector <Matrix3r> lubBodyStress;
vector <Vector3r> shearLubricationForces;
vector <Vector3r> shearLubricationTorques;
vector <Vector3r> pumpLubricationTorques;
vector <Vector3r> twistLubricationTorques;
vector <Vector3r> normalLubricationForce;
vector <Matrix3r> shearLubricationBodyStress;
vector <Matrix3r> normalLubricationBodyStress;
vector <Vector3r> deltaNormVel;
vector <Vector3r> deltaShearVel;
vector <Vector3r> normalV;
......@@ -89,76 +89,78 @@ class FlowBoundingSphere : public Network<_Tesselation>
vector <Matrix3r> normalStressInteraction;
void Localize();
void Compute_Permeability();
virtual void GaussSeidel (Real dt=0);
virtual void ResetNetwork();
void computePermeability();
virtual void gaussSeidel (Real dt=0);
virtual void resetNetwork();
void Fictious_cells ( );
double k_factor; //permeability moltiplicator
double kFactor; //permeability moltiplicator
std::string key; //to give to consolidation files a name with iteration number
std::vector<double> Pressures; //for automatic write maximum pressures during consolidation
bool tess_based_force; //allow the force computation method to be chosen from FlowEngine
// std::vector<double> pressures; //for automatic write maximum pressures during consolidation
bool tessBasedForce; //allow the force computation method to be chosen from FlowEngine
Real minPermLength; //min branch length for Poiseuille
double P_SUP, P_INF, P_INS, VISCOSITY;
double VISCOSITY;
double fluidBulkModulus;
Tesselation& Compute_Action ( );
Tesselation& Compute_Action ( int argc, char *argv[ ], char *envp[ ] );
Tesselation& LoadPositions(int argc, char *argv[ ], char *envp[ ]);
void SpheresFileCreator ();
void DisplayStatistics();
void Initialize_pressures ( double P_zero );
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();
void computeTetrahedralForces();
/// Define forces spliting drag and buoyancy terms
void ComputeFacetForcesWithCache(bool onlyCache=false);
void saveVtk (const char* folder);
void computeFacetForcesWithCache(bool onlyCache=false);
void saveVtk (const char* folder );
#ifdef XVIEW
void Dessine_Triangulation ( Vue3D &Vue, RTriangulation &T );
void Dessine_Short_Tesselation ( Vue3D &Vue, Tesselation &Tes );
void dessineTriangulation ( Vue3D &Vue, RTriangulation &T );
void dessineShortTesselation ( Vue3D &Vue, Tesselation &Tes );
#endif
double Permeameter ( double P_Inf, double P_Sup, double Section, double DeltaY, const char *file );
double Sample_Permeability( double& x_Min,double& x_Max ,double& y_Min,double& y_Max,double& z_Min,double& z_Max);
double Compute_HydraulicRadius (Cell_handle cell, int j );
double permeameter ( double PInf, double PSup, double Section, double DeltaY, const char *file );
double samplePermeability( double& xMin,double& xMax ,double& yMin,double& yMax,double& zMin,double& zMax);
double computeHydraulicRadius (CellHandle cell, int j );
Real checkSphereFacetOverlap(const Sphere& v0, const Sphere& v1, const Sphere& v2);
double dotProduct ( Vecteur x, Vecteur y );
double Compute_EffectiveRadius(Cell_handle cell, int j);
double Compute_EquivalentRadius(Cell_handle cell, int j);
double computeEffectiveRadius(CellHandle cell, int j);
double computeEquivalentRadius(CellHandle cell, int j);
//return the list of constriction values
vector<double> getConstrictions();
vector<Constriction> getConstrictionsFull();
void GenerateVoxelFile ( );
void generateVoxelFile ( );
void computeEdgesSurfaces();
Vector3r computeViscousShearForce(const Vector3r& deltaV, const int& edge_id, const Real& Rh);
Real computeNormalLubricationForce(const Real& deltaNormV, const Real& dist, const int& edge_id, const Real& eps, const Real& stiffness, const Real& dt, const Real& meanRad);
Vector3r computeShearLubricationForce(const Vector3r& deltaShearV, const Real& dist, const int& edge_id, const Real& eps, const Real& centerDist, const Real& meanRad);
Vector3r computePumpTorque(const Vector3r& deltaShearAngV, const Real& dist, const int& edge_id, const Real& eps, const Real& meanRad );
Vector3r computeTwistTorque(const Vector3r& deltaNormAngV, const Real& dist, const int& edge_id, const Real& eps, const Real& meanRad );
RTriangulation& Build_Triangulation ( Real x, Real y, Real z, Real radius, unsigned const id );
RTriangulation& buildTriangulation ( Real x, Real y, Real z, Real radius, unsigned const id );
bool isInsideSphere ( double& x, double& y, double& z );
void SliceField (const char *filename);
void ComsolField();
void sliceField (const char *filename);
void comsolField();
void Interpolate ( Tesselation& Tes, Tesselation& NewTes );
virtual void Average_Relative_Cell_Velocity();
void Average_Fluid_Velocity();
void ApplySinusoidalPressure(RTriangulation& Tri, double Amplitude, double Average_Pressure, double load_intervals);
void interpolate ( Tesselation& Tes, Tesselation& NewTes );
virtual void averageRelativeCellVelocity();
void averageFluidVelocity();
void applySinusoidalPressure(RTriangulation& Tri, double amplitude, double averagePressure, double loadIntervals);
bool isOnSolid (double X, double Y, double Z);
double getPorePressure (double X, double Y, double Z);
void measurePressureProfile(double Wall_up_y, double Wall_down_y);
void measurePressureProfile(double WallUpy, double WallDowny);
double averageSlicePressure(double Y);
double averagePressure();
double getCell (double X,double Y,double Z);
double boundaryFlux(unsigned int boundaryId);
vector<Real> Average_Fluid_Velocity_On_Sphere(unsigned int Id_sph);
vector<Real> averageFluidVelocityOnSphere(unsigned int Id_sph);
//Solver?
int useSolver;//(0 : GaussSeidel, 1 : TAUCS, 2 : PARDISO, 3:CHOLMOD)
};
......
This diff is collapsed.
......@@ -44,9 +44,9 @@ public:
using FlowType::T;
using FlowType::currentTes;
using FlowType::boundary;
using FlowType::y_min_id;
using FlowType::y_max_id;
using FlowType::DEBUG_OUT;
using FlowType::yMinId;
using FlowType::yMaxId;
using FlowType::debugOut;
using FlowType::TOLERANCE;
using FlowType::RELAX;
using FlowType::fluidBulkModulus;
......@@ -55,7 +55,7 @@ public:
using FlowType::computedOnce;
//! TAUCS DECs
vector<Finite_cells_iterator> orderedCells;
vector<FiniteCellsIterator> orderedCells;
bool isLinearSystemSet;
bool isFullLinearSystemGSSet;
bool areCellsOrdered;//true when orderedCells is filled, turn it false after retriangulation
......@@ -97,7 +97,7 @@ public:
vector<int> T_jn;//(size+1);
vector<int> T_ia;//(size*5);
vector<double> T_f;//(size); // right-hand size vector object
vector<Cell_handle> T_cells;//(size)
vector<CellHandle> T_cells;//(size)
int T_index;
vector<double> T_b;
......@@ -143,35 +143,35 @@ public:
FlowBoundingSphereLinSolv();
///Linear system solve
virtual int SetLinearSystem(Real dt);
void VectorizedGaussSeidel(Real dt);
virtual int SetLinearSystemFullGS(Real dt);
virtual int setLinearSystem(Real dt);
void vectorizedGaussSeidel(Real dt);
virtual int setLinearSystemFullGS(Real dt);
int TaucsSolveTest();
int TaucsSolve(Real dt);
int PardisoSolveTest();
int PardisoSolve(Real dt);
int taucsSolveTest();
int taucsSolve(Real dt);
int pardisoSolveTest();
int pardisoSolve(Real dt);
int eigenSolve(Real dt);
void CopyGsToCells();
void CopyCellsToGs(Real dt);
void copyGsToCells();
void copyCellsToGs(Real dt);
void CopyLinToCells();
void CopyCellsToLin(Real dt);
void swap_fwd (double* v, int i);
void swap_fwd (int* v, int i);
void sort_v(int k1, int k2, int* is, double* ds);
void copyLinToCells();
void copyCellsToLin(Real dt);
void swapFwd (double* v, int i);
void swapFwd (int* v, int i);
void sortV(int k1, int k2, int* is, double* ds);
virtual void GaussSeidel (Real dt) {
virtual void gaussSeidel (Real dt) {
switch (useSolver) {
case 0:
VectorizedGaussSeidel(dt);
vectorizedGaussSeidel(dt);
break;
case 1:
TaucsSolve(dt);
taucsSolve(dt);
break;
case 2:
PardisoSolve(dt);
pardisoSolve(dt);
break;
case 3:
eigenSolve(dt);
......@@ -179,7 +179,7 @@ public:
}
computedOnce=true;
}
virtual void ResetNetwork();
virtual void resetNetwork();
};
} //namespace CGT
......
......@@ -46,46 +46,45 @@ class Network
Tesselation T [2];
bool currentTes;
double x_min, x_max, y_min, y_max, z_min, z_max, Rmoy, SectionArea, Height, Vtotale;
bool DEBUG_OUT;
double xMin, xMax, yMin, yMax, zMin, zMax, Rmoy, sectionArea, Height, vTotal;
bool debugOut;
int nOfSpheres;
int x_min_id, x_max_id, y_min_id, y_max_id, z_min_id, z_max_id;
int xMinId, xMaxId, yMinId, yMaxId, zMinId, zMaxId;
int* boundsIds [6];
vector<Cell_handle> boundingCells [6];
Point Corner_min;
Point Corner_max;
Real Vsolid_tot, Vtotalissimo, Vporale, Ssolid_tot, V_porale_porosity, V_totale_porosity;
vector<CellHandle> boundingCells [6];
Point cornerMin;
Point cornerMax;
Real VSolidTot, Vtotalissimo, vPoral, sSolidTot, vPoralPorosity, vTotalePorosity;
Boundary boundaries [6];
Boundary& boundary (int b) {return boundaries[b-id_offset];}
short id_offset;
int vtk_infinite_vertices, vtk_infinite_cells, num_particles;
void AddBoundingPlanes();
void AddBoundingPlane (Vecteur Normal, int id_wall);
void AddBoundingPlane (Real center[3], double thickness, Vecteur Normal, int id_wall );
void Define_fictious_cells( );
int detectFacetFictiousVertices (Cell_handle& cell, int& j);
double volumeSolidPore (const Cell_handle& cell);
double volume_single_fictious_pore(const Vertex_handle& SV1, const Vertex_handle& SV2, const Vertex_handle& SV3, const Point& PV1, const Point& PV2, Vecteur& facetSurface);
double volume_double_fictious_pore(const Vertex_handle& SV1, const Vertex_handle& SV2, const Vertex_handle& SV3, const Point& PV1, const Point& PV2, Vecteur& facetSurface);
double spherical_triangle_volume(const Sphere& ST1, const Point& PT1, const Point& PT2, const Point& PT3);
Boundary& boundary (int b) {return boundaries[b-idOffset];}
short idOffset;
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 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 sphericalTriangleVolume(const Sphere& ST1, const Point& PT1, const Point& PT2, const Point& PT3);
double fast_spherical_triangle_area(const Sphere& STA1, const Point& STA2, const Point& STA3, const Point& PTA1);
Real fast_solid_angle(const Point& STA1, const Point& PTA1, const Point& PTA2, const Point& PTA3);
double volume_double_fictious_pore(Vertex_handle SV1, Vertex_handle SV2, Vertex_handle SV3, Point PV1);
double volume_single_fictious_pore(Vertex_handle SV1, Vertex_handle SV2, Vertex_handle SV3, Point PV1);
double Volume_Pore_VoronoiFraction ( Cell_handle& cell, int& j, bool reuseFacetData=false);
double Surface_Solid_Pore( Cell_handle cell, int j, bool SLIP_ON_LATERALS, bool reuseFacetData=false);
double spherical_triangle_area ( Sphere STA1, Sphere STA2, Sphere STA3, Point PTA1 );
double fastSphericalTriangleArea(const Sphere& STA1, const Point& STA2, const Point& STA3, const Point& PTA1);
Real fastSolidAngle(const Point& STA1, const Point& PTA1, const Point& PTA2, const Point& PTA3);
double volumeDoubleFictiousPore(VertexHandle SV1, VertexHandle SV2, VertexHandle SV3, Point PV1);
double volumeSingleFictiousPore(VertexHandle SV1, VertexHandle SV2, VertexHandle SV3, Point PV1);
double volumePoreVoronoiFraction ( CellHandle& cell, int& j, bool reuseFacetData=false);
double surfaceSolidPore( CellHandle cell, int j, bool slipOnLaterals, bool reuseFacetData=false);
double sphericalTriangleArea ( Sphere STA1, Sphere STA2, Sphere STA3, Point PTA1 );
Vecteur surface_double_fictious_facet(Vertex_handle fSV1, Vertex_handle fSV2, Vertex_handle SV3);
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);
Vecteur surfaceDoubleFictiousFacet(VertexHandle fSV1, VertexHandle fSV2, VertexHandle SV3);
Vecteur surfaceSingleFictiousFacet(VertexHandle fSV1, VertexHandle SV2, VertexHandle SV3);
double surfaceSolidDoubleFictiousFacet(VertexHandle SV1, VertexHandle SV2, VertexHandle SV3);
double surfaceSolidFacet(Sphere ST1, Sphere ST2, Sphere ST3);
int facetF1, facetF2, facetRe1, facetRe2, facetRe3;
int F1, F2, Re1, Re2;
int facetNFictious;
int real_vertex;
double FAR;
......@@ -99,6 +98,4 @@ class Network
#include "Network.ipp"
#endif //FLOW_ENGINE
This diff is collapsed.
This diff is collapsed.
......@@ -20,11 +20,11 @@ namespace CGT{
class PeriodicFlow : public PeriodicFlowBoundingSphere
{
public:
void Interpolate(Tesselation& Tes, Tesselation& NewTes);
void ComputeFacetForcesWithCache(bool onlyCache=false);
void Compute_Permeability();
void GaussSeidel(Real dt=0);
void DisplayStatistics();
void interpolate(Tesselation& Tes, Tesselation& NewTes);
void computeFacetForcesWithCache(bool onlyCache=false);
void computePermeability();
void gaussSeidel(Real dt=0);
void displayStatistics();
void computeEdgesSurfaces();
double boundaryFlux(unsigned int boundaryId);
#ifdef EIGENSPARSE_LIB
......
......@@ -27,8 +27,8 @@ public:
PeriodicFlowLinSolv();
///Linear system solve
virtual int SetLinearSystem(Real dt=0);
virtual int SetLinearSystemFullGS(Real dt=0);
virtual int setLinearSystem(Real dt=0);
virtual int setLinearSystemFullGS(Real dt=0);
};
} //namespace CGTF
......
......@@ -48,7 +48,7 @@ PeriodicFlowLinSolv::~PeriodicFlowLinSolv()
PeriodicFlowLinSolv::PeriodicFlowLinSolv(): LinSolver() {}
int PeriodicFlowLinSolv::SetLinearSystem(Real dt)
int PeriodicFlowLinSolv::setLinearSystem(Real dt)
{
//WARNING : boundary conditions (Pcondition, p values) must have been set for a correct definition of
RTriangulation& Tri = T[currentTes].Triangulation();
......@@ -64,11 +64,10 @@ int PeriodicFlowLinSolv::SetLinearSystem(Real dt)
unsigned int maxindex = 0;
//FIXME: this is way too large since many cells will be ghosts
T_cells.resize(Tri.number_of_finite_cells()+1);
// indices.resize(Tri.number_of_finite_cells()+1);
///Ordered cells
orderedCells.clear();
const Finite_cells_iterator cell_end = Tri.finite_cells_end();
for (Finite_cells_iterator cell = Tri.finite_cells_begin(); cell != cell_end; cell++) {
const FiniteCellsIterator cellEnd = Tri.finite_cells_end();
for (FiniteCellsIterator cell = Tri.finite_cells_begin(); cell != cellEnd; cell++) {
if (cell->info().Pcondition || cell->info().isGhost) continue;
orderedCells.push_back(cell);
// T_cells[++ncols]=cell;
......@@ -79,7 +78,6 @@ int PeriodicFlowLinSolv::SetLinearSystem(Real dt)
}
// spatial_sort(orderedCells.begin(),orderedCells.end(), CellTraits_for_spatial_sort());//FIXME: ordering will not work with the new "indices", which needs reordering to
T_cells.resize(ncols+1);
// indices.resize(maxindex+1);
isLinearSystemSet=false;
areCellsOrdered=true;
}
......@@ -93,48 +91,45 @@ int PeriodicFlowLinSolv::SetLinearSystem(Real dt)
T_bv.resize(ncols);
bodv.resize(ncols);
xodv.resize(ncols);
// gsB.resize(ncols+1);
T_cells.resize(ncols+1);
T_nnz=0;}
for (int kk=0; kk<ncols;kk++) T_b[kk]=0;
///Ordered cells
int nIndex=0; Cell_handle neighbour_cell;
int nIndex=0; CellHandle neighbourCell;
for (int i=0; i<ncols; i++)
{
Finite_cells_iterator& cell = orderedCells[i];
FiniteCellsIterator& cell = orderedCells[i];
///Non-ordered cells
// for (Finite_cells_iterator cell = Tri.finite_cells_begin(); cell != cell_end; cell++) {
// if (!cell->info().Pcondition && !cell->info().isGhost)
{
const int& index=cell->info().index;
const Cell_Info& cInfo = cell->info();
const CellInfo& cInfo = cell->info();
if (!isLinearSystemSet) {
//Add diagonal term
is[T_nnz] = index;
js[T_nnz] = index;
vs[T_nnz] = (cInfo.k_norm())[0]+ (cInfo.k_norm())[1]+ (cInfo.k_norm())[2]+ (cInfo.k_norm())[3];
vs[T_nnz] = (cInfo.kNorm())[0]+ (cInfo.kNorm())[1]+ (cInfo.kNorm())[2]+ (cInfo.kNorm())[3];
if (vs[T_nnz]<0) cerr<<"!!!! WTF !!!"<<endl;
if (fluidBulkModulus>0) vs[T_nnz] += (1.f/(dt*fluidBulkModulus*cInfo.invVoidVolume()));
T_nnz++;
}
for (int j=0; j<4; j++) {
neighbour_cell = cell->neighbor(j);
if (Tri.is_infinite(neighbour_cell)) continue;
Cell_Info& nInfo = neighbour_cell->info();
neighbourCell = cell->neighbor(j);
if (Tri.is_infinite(neighbourCell)) continue;
CellInfo& nInfo = neighbourCell->info();
nIndex=nInfo.index;
if (nIndex==index) {
cerr<<"ERROR: nIndex==index, the cell is neighbour to itself"<< index<<endl;
errorCode=3;}
if (nInfo.Pcondition) T_b[index-1]+=cInfo.k_norm()[j]*nInfo.shiftedP();
if (nInfo.Pcondition) T_b[index-1]+=cInfo.kNorm()[j]*nInfo.shiftedP();
else {
if (!isLinearSystemSet && index>nIndex) {
is[T_nnz] = index;
js[T_nnz] = nIndex;
vs[T_nnz] = - (cInfo.k_norm())[j];
vs[T_nnz] = - (cInfo.kNorm())[j];
if (vs[T_nnz]>0) cerr<<"!!!! WTF2 !!!"<<endl;
T_nnz++;}
if (nInfo.isGhost) T_b[index-1]+=cInfo.k_norm()[j]*nInfo.pShift();
if (nInfo.isGhost) T_b[index-1]+=cInfo.kNorm()[j]*nInfo.pShift();
}
}
}
......@@ -203,16 +198,9 @@ int PeriodicFlowLinSolv::SetLinearSystem(Real dt)
tripletList.clear(); tripletList.resize(T_nnz);
for(int k=0;k<T_nnz;k++) {
tripletList[k]=ETriplet(is[k]-1,js[k]-1,vs[k]);
// file<<is[k]-1<<" "<<js[k]-1<<" "<<vs[k]<<endl;
// if (is[k]==js[k]) file2<<is[k]-1<<" "<<js[k]-1<<" "<<1.0001*vs[k]<<endl;
// else file2<<is[k]-1<<" "<<js[k]-1<<" "<<vs[k]<<endl;
// if (is[k]==js[k]) file3<<is[k]-1<<" "<<js[k]-1<<" "<<1.00000001*vs[k]<<endl;
// else file3<<is[k]-1<<" "<<js[k]-1<<" "<<vs[k]<<endl;
}
A.resize(ncols,ncols);
A.setFromTriplets(tripletList.begin(), tripletList.end());
// file << A;
// file.close();
#else
cerr<<"yade compiled without CHOLMOD, FlowEngine.useSolver="<< useSolver <<" not supported"<<endl;
#endif
......@@ -226,7 +214,7 @@ int PeriodicFlowLinSolv::SetLinearSystem(Real dt)
/// For Gauss Seidel, we need the full matrix
int PeriodicFlowLinSolv::SetLinearSystemFullGS(Real dt)
int PeriodicFlowLinSolv::setLinearSystemFullGS(Real dt)
{
//WARNING : boundary conditions (Pcondition, p values) must have been set for a correct definition
RTriangulation& Tri = T[currentTes].Triangulation();
......@@ -237,10 +225,10 @@ int PeriodicFlowLinSolv::SetLinearSystemFullGS(Real dt)
T_index=0;//FIXME : no need to clear if we don't re-triangulate
T_nnz=0;
ncols=0;
const Finite_cells_iterator cell_end = Tri.finite_cells_end();
const FiniteCellsIterator cellEnd = Tri.finite_cells_end();
orderedCells.clear();
T_cells.resize(n_cells+1);
for (Finite_cells_iterator cell = Tri.finite_cells_begin(); cell != cell_end; cell++) {
for (FiniteCellsIterator cell = Tri.finite_cells_begin(); cell != cellEnd; cell++) {
if (cell->info().Pcondition || cell->info().isGhost) continue;
++ncols;
T_cells[cell->info().index]=cell;
......@@ -266,53 +254,51 @@ int PeriodicFlowLinSolv::SetLinearSystemFullGS(Real dt)
///Ordered cells
for (int i=1; i<=ncols; i++)
{
Cell_handle& cell = T_cells[i];
CellHandle& cell = T_cells[i];
///Non-ordered cells
if (!cell->info().Pcondition && !cell->info().isGhost) {
//Add diagonal term
fullAvalues[cell->info().index][4] = 1.f/((cell->info().k_norm())[0]+ (cell->info().k_norm())[1]+ (cell->info().k_norm())[2]+ (cell->info().k_norm())[3] + (fluidBulkModulus>0? 1.f/(dt*fluidBulkModulus*cell->info().invVoidVolume()):0));
fullAvalues[cell->info().index][4] = 1.f/((cell->info().kNorm())[0]+ (cell->info().kNorm())[1]+ (cell->info().kNorm())[2]+ (cell->info().kNorm())[3] + (fluidBulkModulus>0? 1.f/(dt*fluidBulkModulus*cell->info().invVoidVolume()):0));
//DUMP
// cout<< cell->info().index<<" "<< cell->info().index<<" "<<fullAvalues[cell->info().index][4] <<endl;
T_nnz++;
for (int j=0; j<4; j++) {
Cell_handle neighbour_cell = cell->neighbor(j);
const Cell_Info& nInfo = neighbour_cell->info();
Cell_Info& cInfo = cell->info();
if (Tri.is_infinite(neighbour_cell)) {
CellHandle neighbourCell = cell->neighbor(j);
const CellInfo& nInfo = neighbourCell->info();
CellInfo& cInfo = cell->info();
if (Tri.is_infinite(neighbourCell)) {
fullAvalues[cInfo.index][j] = 0;
//We point to the pressure of the adjacent cell. If the cell is ghost, then it has the index of the real one, and then the pointer is correct
fullAcolumns[cInfo.index][j] = &gsP[0];
continue;}
if (!nInfo.Pcondition) {
++T_nnz;
fullAvalues[cInfo.index][j] = (cInfo.k_norm())[j];
fullAvalues[cInfo.index][j] = (cInfo.kNorm())[j];
fullAcolumns[cInfo.index][j] = &gsP[nInfo.index];
//DUMP
// cout<< cInfo.index<<" "<< nInfo.index<<" "<<fullAvalues[cInfo.index][j] <<endl;
//if the adjacent cell is ghost, we account for the pressure shift in the RHS
if (nInfo.isGhost){
gsB[cInfo.index]+=cInfo.k_norm()[j]*nInfo.pShift();
gsB[cInfo.index]+=cInfo.kNorm()[j]*nInfo.pShift();
}
} else {
fullAvalues[cInfo.index][j] = 0;
fullAcolumns[cInfo.index][j] = &gsP[0];
gsB[cInfo.index]+=cInfo.k_norm()[j]*nInfo.shiftedP();
gsB[cInfo.index]+=cInfo.kNorm()[j]*nInfo.shiftedP();
}
}
}
}
} else for (int i=1; i<=ncols; i++)
{
Cell_handle& cell = T_cells[i];
CellHandle& cell = T_cells[i];
///Non-ordered cells
if (!cell->info().Pcondition && !cell->info().isGhost) {
for (int j=0; j<4; j++) {
Cell_handle neighbour_cell = cell->neighbor(j);
const Cell_Info& nInfo = neighbour_cell->info();
Cell_Info& cInfo = cell->info();
CellHandle neighbourCe