------------------------------------------------------------ revno: 3418 committer: Chao Yuan <[email protected]> timestamp: Fri 2014-05-09 16:29:16 +0200 message: -big change, inherit from FlowEngine. added: pkg/pfv/UnsaturatedEngine.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
=== added file 'pkg/pfv/UnsaturatedEngine.cpp' --- pkg/pfv/UnsaturatedEngine.cpp 1970-01-01 00:00:00 +0000 +++ pkg/pfv/UnsaturatedEngine.cpp 2014-05-09 14:29:16 +0000 @@ -0,0 +1,1065 @@ +/************************************************************************* +* Copyright (C) 2012 by Chao Yuan <[email protected]> * +* Copyright (C) 2012 by Bruno Chareyre <[email protected]> * +* * +* 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. * +*************************************************************************/ +#ifdef FLOW_ENGINE + +//keep this #ifdef for commited versions unless you really have stable version that should be compiled by default +//it will save compilation time for everyone else +//when you want it compiled, you can pass -DDFNFLOW to cmake, or just uncomment the following line +#define UNSATURATED_FLOW +#ifdef UNSATURATED_FLOW +#define TEMPLATE_FLOW_NAME UnsaturatedEngineT +#include <yade/pkg/pfv/FlowEngine.hpp> + +/// We can add data to the Info types by inheritance +class UnsatCellInfo : public FlowCellInfo { + public: + bool isWaterReservoir; + bool isAirReservoir; + Real capillaryCellVolume;//abs(cell.volume) - abs(cell.solid.volume) + std::vector<double> poreRadius;//pore throat radius for drainage + double solidLine [4][4];//the length of intersecting line between sphere and facet. [i][j] is for sphere facet "i" and sphere facetVertices[i][j]. Last component for 1/sumLines in the facet. + double trapCapP;//for calculating the pressure of trapped phase, cell.pressureTrapped = pressureAir - trapCapP. + int windowsID;//a temp cell info for experiment comparison + UnsatCellInfo (void) + { + poreRadius.resize(4, 0); + isWaterReservoir = false; isAirReservoir = false; capillaryCellVolume = 0; + for (int k=0; k<4;k++) for (int l=0; l<3;l++) solidLine[k][l]=0; + trapCapP = 0; + windowsID = 0; + } +}; + +class UnsatVertexInfo : public FlowVertexInfo { +// add later; +public: +// UnsatVertexInfo (void) +}; + +typedef TemplateFlowEngine<UnsatCellInfo,UnsatVertexInfo> UnsaturatedEngineT; +REGISTER_SERIALIZABLE(UnsaturatedEngineT); +YADE_PLUGIN((UnsaturatedEngineT)); + +class UnsaturatedEngine : public UnsaturatedEngineT +{ + protected: + void testFunction(); + + + public : + void initializeCellIndex(); + void updatePoreRadius(); + void updateVolumeCapillaryCell(); + double computeEffPoreRadius(CellHandle cell, int j); + double computeEffPoreRadiusFine(CellHandle cell, int j); + double bisection(CellHandle cell, int j, double a, double b); + double computeDeltaForce(CellHandle cell,int j, double rCap); + void computeSolidLine(); + void computeFacetPoreForcesWithCache(bool onlyCache=false); + + void invade(); + Real getMinEntryValue(); + Real getSaturation(); + + void invadeSingleCell1(CellHandle cell, double pressure); + void invade1(); + void checkTrap(double pressure);//check trapped phase, define trapCapP. + Real getMinEntryValue1(); + void updatePressure(); + void updatePressureReservoir(); + void initReservoirBound(); + void initWaterReservoirBound(); + void updateWaterReservoir(); + void waterReservoirRecursion(CellHandle cell); + void initAirReservoirBound(); + void updateAirReservoir(); + void airReservoirRecursion(CellHandle cell); + Real getSaturation1(); + + void invadeSingleCell2(CellHandle cell, double pressure); + void invade2(); + void updatePressure2(); + Real getMinEntryValue2(); + Real getSaturation2(); + + //record and test functions + void checkCellsConnection(); + void checkEntryCapillaryPressure(); + void checkLatticeNodeY(double y); + void checkReservoirInfo(int boundN); + void checkBoundingCellsInfo(); + void saveVtk(const char* folder) {bool initT=solver->noCache; solver->noCache=false; solver->saveVtk(folder); solver->noCache=initT;} + //temp functions + void initializeCellWindowsID(); + Real getWindowsSaturation(int windowsID); + Real getWindowsSaturation1(int i); + Real getWindowsSaturation2(int i); + double getRadiusMin(CellHandle cell, int j); + void debugTemp(); + + virtual ~UnsaturatedEngine(); + + virtual void action(); + + YADE_CLASS_BASE_DOC_ATTRS_INIT_CTOR_PY(UnsaturatedEngine,UnsaturatedEngineT,"Preliminary version engine of a drainage model for unsaturated soils. Note:Air reservoir is on the top; water reservoir is on the bottom.", + ((bool, isPhaseTrapped, true,,"Activate invade mode. If True, the wetting phase can be trapped, activate invade mode 1; if false, the wetting phase cann't be trapped, activate invade mode 2.")) + ((double,gasPressure,0,,"Invasion pressure")) + ((double,surfaceTension,0.0728,,"Water Surface Tension in contact with air at 20 Degrees Celsius is: 0.0728(N/m)")) + + ((bool, computeForceActivated, true,,"Activate capillary force computation. WARNING: turning off means capillary force is not computed at all, but the drainage can still work.")) + ((bool, invadeBoundary, false,,"Invade from boundaries.")) + ((int, windowsNo, 10,, "Number of genrated windows/(zoomed samples).")) + ,,, + + .def("saveVtk",&UnsaturatedEngine::saveVtk,(python::arg("folder")="./VTK"),"Save pressure field in vtk format. Specify a folder name for output.") + .def("testFunction",&UnsaturatedEngine::testFunction,"The playground for Chao's experiments.") + .def("getSaturation",&UnsaturatedEngine::getSaturation,"get saturation.") + .def("getWindowsSaturation",&UnsaturatedEngine::getWindowsSaturation,(python::arg("windowsID")), "get saturation of windowsID") + .def("getMinEntryValue",&UnsaturatedEngine::getMinEntryValue,"get the minimum air entry pressure for the next invade step.") + .def("checkCellsConnection",&UnsaturatedEngine::checkCellsConnection,"Check cell connections.") + .def("checkEntryCapillaryPressure",&UnsaturatedEngine::checkEntryCapillaryPressure,"Check entry capillary pressure between neighbor cells.") + .def("checkLatticeNodeY",&UnsaturatedEngine::checkLatticeNodeY,(python::arg("y")),"Check the slice of lattice nodes for yNormal(y). 0: out of sphere; 1: inside of sphere.") + .def("checkReservoirInfo",&UnsaturatedEngine::checkReservoirInfo,(python::arg("boundN")),"Check reservoir cells(N=2,3) statement and export to 'waterReservoirBoundInfo.txt' and 'airReservoirBoundInfo.txt'.") + .def("checkBoundingCellsInfo",&UnsaturatedEngine::checkBoundingCellsInfo,"Check boundary cells (without reservoirs) statement and export to 'boundInfo.txt'.") + .def("debugTemp",&UnsaturatedEngine::debugTemp,"debug temp file.(temporary)") + .def("initializeCellWindowsID",&UnsaturatedEngine::initializeCellWindowsID,"Initialize cell windows index. A temp function for comparison with experiments, will delete soon") + .def("invade",&UnsaturatedEngine::invade,"Run the drainage invasion.") + .def("computeForce",&UnsaturatedEngine::computeFacetPoreForcesWithCache,"Compute capillary force. ") + ) + DECLARE_LOGGER; +}; + +REGISTER_SERIALIZABLE(UnsaturatedEngine); +YADE_PLUGIN((UnsaturatedEngine)); + +UnsaturatedEngine::~UnsaturatedEngine(){} + +void UnsaturatedEngine::testFunction() +{ + cout<<"This is UnsaturatedEngine test program"<<endl; + RTriangulation& tri = solver->T[solver->currentTes].Triangulation(); + if (tri.number_of_vertices()==0) { + cout<< "triangulation is empty: building a new one" << endl; + scene = Omega::instance().getScene().get();//here define the pointer to Yade's scene + setPositionsBuffer(true);//copy sphere positions in a buffer... + buildTriangulation(pZero,*solver);//create a triangulation and initialize pressure in the elements, everything will be contained in "solver" + initializeCellIndex();//initialize cell index + updatePoreRadius();//save all pore radii before invade + updateVolumeCapillaryCell();//save capillary volume of all cells, for fast calculating saturation + computeSolidLine();//save cell->info().solidLine[j][y] + } + solver->noCache = true; +} + +void UnsaturatedEngine::action() +{/* + if ( !isActivated ) return; + RTriangulation& tri = solver->T[solver->currentTes].Triangulation(); + if ( (tri.number_of_vertices()==0) || (updateTriangulation) ) { + cout<< "triangulation is empty: building a new one" << endl; + scene = Omega::instance().getScene().get();//here define the pointer to Yade's scene + setPositionsBuffer(true);//copy sphere positions in a buffer... + buildTriangulation(pZero,solver);//create a triangulation and initialize pressure in the elements, everything will be contained in "solver" + initializeCellIndex();//initialize cell index + updatePoreRadius();//save all pore radii before invade + updateVolumeCapillaryCell();//save capillary volume of all cells, for calculating saturation + computeSolidLine();//save cell->info().solidLine[j][y] + solver->noCache = true; + } + ///compute invade + if (pressureForce) { + invade(); + } + + ///compute force + if(computeForceActivated){ + computeFacetPoreForcesWithCache(solver); + Vector3r force; + FiniteVerticesIterator vertices_end = solver->T[solver->currentTes].Triangulation().finite_vertices_end(); + for ( FiniteVerticesIterator V_it = solver->T[solver->currentTes].Triangulation().finite_vertices_begin(); V_it != vertices_end; V_it++ ) { + force = pressureForce ? Vector3r ( V_it->info().forces[0],V_it->info().forces[1],V_it->info().forces[2] ): Vector3r(0,0,0); + scene->forces.addForce ( V_it->info().id(), force); + } +}*/ +} + +void UnsaturatedEngine::initializeCellIndex() +{ + int k=0; + RTriangulation& tri = solver->T[solver->currentTes].Triangulation(); + FiniteCellsIterator cellEnd = tri.finite_cells_end(); + for ( FiniteCellsIterator cell = tri.finite_cells_begin(); cell != cellEnd; cell++ ) { + cell->info().index=k++; + } +} + +void UnsaturatedEngine::updatePoreRadius() +{ + RTriangulation& tri = solver->T[solver->currentTes].Triangulation(); + FiniteCellsIterator cellEnd = tri.finite_cells_end(); + CellHandle neighbourCell; + for (FiniteCellsIterator cell = tri.finite_cells_begin(); cell != cellEnd; cell++) { + for (int j=0; j<4; j++) { + neighbourCell = cell->neighbor(j); + if (!tri.is_infinite(neighbourCell)) { + cell->info().poreRadius[j]=computeEffPoreRadius(cell, j); + neighbourCell->info().poreRadius[tri.mirror_index(cell, j)]= cell->info().poreRadius[j]; + } + } + } +} + +void UnsaturatedEngine::updateVolumeCapillaryCell() +{ + initializeVolumes(*solver); + RTriangulation& tri = solver->T[solver->currentTes].Triangulation(); + FiniteCellsIterator cellEnd = tri.finite_cells_end(); + CellHandle neighbourCell; + for (FiniteCellsIterator cell = tri.finite_cells_begin(); cell != cellEnd; cell++) { + cell->info().capillaryCellVolume = abs( cell->info().volume() ) - solver->volumeSolidPore(cell); + } +} + +void UnsaturatedEngine::invade() +{ + if (isPhaseTrapped) invade1(); + else invade2(); +} + +Real UnsaturatedEngine::getMinEntryValue() +{ + if (isPhaseTrapped) return getMinEntryValue1(); + else return getMinEntryValue2(); +} + +Real UnsaturatedEngine::getSaturation() +{ + if (isPhaseTrapped) return getSaturation1(); + else return getSaturation2(); +} + +///invade mode 1. update phase reservoir before invasion. Consider no viscous effects, and invade gradually. +void UnsaturatedEngine::invadeSingleCell1(CellHandle cell, double pressure) +{ + if (invadeBoundary==true) { + for (int facet = 0; facet < 4; facet ++) { + if (solver->T[solver->currentTes].Triangulation().is_infinite(cell->neighbor(facet))) continue; + if (cell->neighbor(facet)->info().Pcondition) continue; + if (cell->neighbor(facet)->info().isWaterReservoir == true) { + double nCellP = surfaceTension/cell->info().poreRadius[facet]; + if (pressure > nCellP) { + CellHandle nCell = cell->neighbor(facet); + nCell->info().p() = pressure; + nCell->info().isAirReservoir=true; + nCell->info().isWaterReservoir=false; + invadeSingleCell1(nCell, pressure);}}}} + else { + for (int facet = 0; facet < 4; facet ++) { + if (solver->T[solver->currentTes].Triangulation().is_infinite(cell->neighbor(facet))) continue; + if (cell->neighbor(facet)->info().Pcondition) continue;//FIXME:defensive + if (cell->neighbor(facet)->info().isFictious) continue; + if (cell->neighbor(facet)->info().isWaterReservoir == true) { + double nCellP = surfaceTension/cell->info().poreRadius[facet]; + if (pressure > nCellP) { + CellHandle nCell = cell->neighbor(facet); + nCell->info().p() = pressure; + nCell->info().isAirReservoir=true; + nCell->info().isWaterReservoir=false; + invadeSingleCell1(nCell, pressure);}}}} +} + +void UnsaturatedEngine::invade1() +{ + updatePressureReservoir(); + RTriangulation& tri = solver->T[solver->currentTes].Triangulation(); + FiniteCellsIterator cellEnd = tri.finite_cells_end(); + for ( FiniteCellsIterator cell = tri.finite_cells_begin(); cell != cellEnd; cell++ ) { + if(cell->info().isAirReservoir == true) + invadeSingleCell1(cell,cell->info().p()); + } + checkTrap(bndCondValue[3]); + FiniteCellsIterator ncellEnd = tri.finite_cells_end(); + for ( FiniteCellsIterator ncell = tri.finite_cells_begin(); ncell != ncellEnd; ncell++ ) { + if( (ncell->info().isWaterReservoir) || (ncell->info().isAirReservoir) ) continue; + ncell->info().p() = bndCondValue[3] - ncell->info().trapCapP; + } + initReservoirBound(); +} + +//check trapped phase, define trapCapP. +void UnsaturatedEngine::checkTrap(double pressure) +{ + RTriangulation& tri = solver->T[solver->currentTes].Triangulation(); + FiniteCellsIterator cellEnd = tri.finite_cells_end(); + for ( FiniteCellsIterator cell = tri.finite_cells_begin(); cell != cellEnd; cell++ ) { + if( (cell->info().isWaterReservoir) || (cell->info().isAirReservoir) ) continue; + if(cell->info().trapCapP!=0) continue; + cell->info().trapCapP=pressure; + } +} + +Real UnsaturatedEngine::getMinEntryValue1() +{ + updatePressureReservoir(); + Real nextEntry = 1e50; + RTriangulation& tri = solver->T[solver->currentTes].Triangulation(); + if (invadeBoundary==true) { + FiniteCellsIterator cellEnd = tri.finite_cells_end(); + for ( FiniteCellsIterator cell = tri.finite_cells_begin(); cell != cellEnd; cell++ ) { + if (cell->info().isAirReservoir == true) { + for (int facet=0; facet<4; facet ++) { + if (tri.is_infinite(cell->neighbor(facet))) continue; + if (cell->neighbor(facet)->info().Pcondition) continue; + if ( cell->neighbor(facet)->info().isWaterReservoir == true ) { + double nCellP = surfaceTension/cell->info().poreRadius[facet]; + nextEntry = min(nextEntry,nCellP);}}}}} + else { + FiniteCellsIterator cellEnd = tri.finite_cells_end(); + for ( FiniteCellsIterator cell = tri.finite_cells_begin(); cell != cellEnd; cell++ ) { + if (cell->info().isAirReservoir == true) { + for (int facet=0; facet<4; facet ++) { + if (tri.is_infinite(cell->neighbor(facet))) continue; + if (cell->neighbor(facet)->info().Pcondition) continue;//FIXME:defensive + if (cell->neighbor(facet)->info().isFictious) continue; + if ( cell->neighbor(facet)->info().isWaterReservoir == true ) { + double nCellP = surfaceTension/cell->info().poreRadius[facet]; + nextEntry = min(nextEntry,nCellP);}}}}} + if (nextEntry==1e50) { + cout << "End drainage !" << endl; + return nextEntry=0;} + else return nextEntry; +} + +void UnsaturatedEngine::updatePressure() +{ + boundaryConditions(*solver); + solver->pressureChanged=true; + solver->reApplyBoundaryConditions(); + RTriangulation& tri = solver->T[solver->currentTes].Triangulation(); + FiniteCellsIterator cellEnd = tri.finite_cells_end(); + for ( FiniteCellsIterator cell = tri.finite_cells_begin(); cell != cellEnd; cell++ ) { + if (cell->info().isWaterReservoir==true) cell->info().p()=bndCondValue[2]; + if (cell->info().isAirReservoir==true) cell->info().p()=bndCondValue[3]; + } +} + +//update pressure and reservoir attr +void UnsaturatedEngine::updatePressureReservoir() +{ + updatePressure();//NOTE:updatePressure must be run before update reservoirs. + updateAirReservoir(); + updateWaterReservoir(); +} + +//NOTE:keep boundingCells[2],boundingCells[3] always being reservoirs. +void UnsaturatedEngine::initReservoirBound() +{ + initWaterReservoirBound(); + initAirReservoirBound(); +} + +//boundingCells[2] is water reservoir. +void UnsaturatedEngine::initWaterReservoirBound() +{ + if (solver->boundingCells[2].size()==0) { + cerr<<"ERROR! set bndCondIsPressure[2] true. boundingCells.size=0!"; + } + else { + vector<CellHandle>::iterator it = solver->boundingCells[2].begin(); + for ( it ; it != solver->boundingCells[2].end(); it++) { + if ((*it)->info().index == 0) continue; + (*it)->info().isWaterReservoir = true; + } + } +} +void UnsaturatedEngine::updateWaterReservoir() +{ + RTriangulation& tri = solver->T[solver->currentTes].Triangulation(); + FiniteCellsIterator cellEnd = tri.finite_cells_end(); + for ( FiniteCellsIterator cell = tri.finite_cells_begin(); cell != cellEnd; cell++ ) { + cell->info().isWaterReservoir = false; + } + initWaterReservoirBound(); + vector<CellHandle>::iterator it = solver->boundingCells[2].begin(); + for ( it ; it != solver->boundingCells[2].end(); it++) { + if ((*it)->info().index == 0) continue; + waterReservoirRecursion(*it); + } +} +void UnsaturatedEngine::waterReservoirRecursion(CellHandle cell) +{ + for (int facet = 0; facet < 4; facet ++) { + if (solver->T[solver->currentTes].Triangulation().is_infinite(cell->neighbor(facet))) continue; + if (cell->neighbor(facet)->info().p() != bndCondValue[2]) continue; + if (cell->neighbor(facet)->info().isWaterReservoir==true) continue; + CellHandle nCell = cell->neighbor(facet); + nCell->info().isWaterReservoir = true; + waterReservoirRecursion(nCell); + } +} +//boundingCells[3] is air reservoir +void UnsaturatedEngine::initAirReservoirBound() +{ + if (solver->boundingCells[3].size()==0) { + cerr<<"ERROR! set bndCondIsPressure[3] true. boundingCells.size=0!"; + } + else { + vector<CellHandle>::iterator it = solver->boundingCells[3].begin(); + for ( it ; it != solver->boundingCells[3].end(); it++) { + if((*it)->info().index == 0) continue; + (*it)->info().isAirReservoir = true; + } + } +} + +void UnsaturatedEngine::updateAirReservoir() +{ + RTriangulation& tri = solver->T[solver->currentTes].Triangulation(); + FiniteCellsIterator cellEnd = tri.finite_cells_end(); + for ( FiniteCellsIterator cell = tri.finite_cells_begin(); cell != cellEnd; cell++ ) { + cell->info().isAirReservoir = false; + } + initAirReservoirBound(); + vector<CellHandle>::iterator it = solver->boundingCells[3].begin(); + for ( it ; it != solver->boundingCells[3].end(); it++) { + if ((*it)->info().index == 0) continue; + airReservoirRecursion(*it); + } +} + +void UnsaturatedEngine::airReservoirRecursion(CellHandle cell) +{ + for (int facet = 0; facet < 4; facet ++) { + if (solver->T[solver->currentTes].Triangulation().is_infinite(cell->neighbor(facet))) continue; + if (cell->neighbor(facet)->info().p() != bndCondValue[3]) continue; + if (cell->neighbor(facet)->info().isAirReservoir == true) continue; + CellHandle nCell = cell->neighbor(facet); + nCell->info().isAirReservoir = true; + airReservoirRecursion(nCell); + } +} + +Real UnsaturatedEngine::getSaturation1() +{ + updatePressureReservoir(); + RTriangulation& tri = solver->T[solver->currentTes].Triangulation(); + Real capillaryVolume = 0.0; //total capillary volume + Real airVolume = 0.0; //air volume + FiniteCellsIterator cellEnd = tri.finite_cells_end(); + + if (invadeBoundary==true) { + for ( FiniteCellsIterator cell = tri.finite_cells_begin(); cell != cellEnd; cell++ ) { + if (tri.is_infinite(cell)) continue; + if (cell->info().Pcondition) continue;//NOTE:reservoirs cells should not be included in saturation + capillaryVolume = capillaryVolume + cell->info().capillaryCellVolume; + if (cell->info().isAirReservoir==true) { + airVolume = airVolume + cell->info().capillaryCellVolume;}}} + else { + for ( FiniteCellsIterator cell = tri.finite_cells_begin(); cell != cellEnd; cell++ ) { + if (tri.is_infinite(cell)) continue; + if (cell->info().Pcondition) continue; + if (cell->info().isFictious) continue; + capillaryVolume = capillaryVolume + cell->info().capillaryCellVolume; + if (cell->info().isAirReservoir==true) { + airVolume = airVolume + cell->info().capillaryCellVolume;}}} + Real saturation = 1 - airVolume/capillaryVolume; + return saturation; +} + +///invade mode 2. Consider no trapped phase. +void UnsaturatedEngine::invadeSingleCell2(CellHandle cell, double pressure) +{ + if (invadeBoundary==true) { + for (int facet = 0; facet < 4; facet ++) { + if (solver->T[solver->currentTes].Triangulation().is_infinite(cell->neighbor(facet))) continue; + if (cell->neighbor(facet)->info().Pcondition) continue; + if (cell->neighbor(facet)->info().p()!=0) continue; + double nCellP = surfaceTension/cell->info().poreRadius[facet]; + if (pressure > nCellP) { + CellHandle nCell = cell->neighbor(facet); + nCell->info().p() = pressure; + invadeSingleCell2(nCell, pressure);}}} + else { + for (int facet = 0; facet < 4; facet ++) { + if (solver->T[solver->currentTes].Triangulation().is_infinite(cell->neighbor(facet))) continue; + if (cell->neighbor(facet)->info().Pcondition) continue;//FIXME:defensive + if (cell->neighbor(facet)->info().isFictious) continue; + if (cell->neighbor(facet)->info().p()!=0) continue; + double nCellP = surfaceTension/cell->info().poreRadius[facet]; + if (pressure > nCellP) { + CellHandle nCell = cell->neighbor(facet); + nCell->info().p() = pressure; + invadeSingleCell2(nCell, pressure);}}} +} + +void UnsaturatedEngine::invade2() +{ + updatePressure2(); + RTriangulation& tri = solver->T[solver->currentTes].Triangulation(); + FiniteCellsIterator cellEnd = tri.finite_cells_end(); + for ( FiniteCellsIterator cell = tri.finite_cells_begin(); cell != cellEnd; cell++ ) { + if (cell->info().p()!=0) { + invadeSingleCell2(cell, cell->info().p()); + } + } +} + +void UnsaturatedEngine::updatePressure2() +{ + boundaryConditions(*solver); + solver->pressureChanged=true; + solver->reApplyBoundaryConditions(); + RTriangulation& tri = solver->T[solver->currentTes].Triangulation(); + FiniteCellsIterator cellEnd = tri.finite_cells_end(); + for ( FiniteCellsIterator cell = tri.finite_cells_begin(); cell != cellEnd; cell++ ) { + if (cell->info().p()!=0) cell->info().p()=bndCondValue[3]; + } +} + +Real UnsaturatedEngine::getMinEntryValue2() +{ + Real nextEntry = 1e50; + RTriangulation& tri = solver->T[solver->currentTes].Triangulation(); + FiniteCellsIterator cellEnd = tri.finite_cells_end(); + + if (invadeBoundary==true) { + for ( FiniteCellsIterator cell = tri.finite_cells_begin(); cell != cellEnd; cell++ ) { + if (cell->info().p()!=0) { + for (int facet=0; facet<4; facet ++) { + if (tri.is_infinite(cell->neighbor(facet))) continue; + if (cell->neighbor(facet)->info().Pcondition) continue; + if (cell->neighbor(facet)->info().p()==0) { + double nCellP = surfaceTension/cell->info().poreRadius[facet]; + nextEntry = min(nextEntry,nCellP);}}}}} + else { + for ( FiniteCellsIterator cell = tri.finite_cells_begin(); cell != cellEnd; cell++ ) { + if (cell->info().p()!=0) { + for (int facet=0; facet<4; facet ++) { + if (tri.is_infinite(cell->neighbor(facet))) continue; + if (cell->neighbor(facet)->info().Pcondition) continue; + if (cell->neighbor(facet)->info().isFictious) continue;//FIXME:defensive + if (cell->neighbor(facet)->info().p()==0) { + double nCellP = surfaceTension/cell->info().poreRadius[facet]; + nextEntry = min(nextEntry,nCellP);}}}}} + if (nextEntry==1e50) { + cout << "End drainage !" << endl; + return nextEntry=0; + } + else { + return nextEntry; + } +} + +Real UnsaturatedEngine::getSaturation2() +{ + RTriangulation& tri = solver->T[solver->currentTes].Triangulation(); + Real capillaryVolume = 0.0; + Real waterVolume = 0.0; + FiniteCellsIterator cellEnd = tri.finite_cells_end(); + + if (invadeBoundary==true) { + for ( FiniteCellsIterator cell = tri.finite_cells_begin(); cell != cellEnd; cell++ ) { + if (tri.is_infinite(cell)) continue; + if (cell->info().Pcondition) continue; + capillaryVolume = capillaryVolume + cell->info().capillaryCellVolume; + if (cell->info().p()==0) { + waterVolume = waterVolume + cell->info().capillaryCellVolume;}}} + else { + for ( FiniteCellsIterator cell = tri.finite_cells_begin(); cell != cellEnd; cell++ ) { + if (tri.is_infinite(cell)) continue; + if (cell->info().Pcondition) continue; + if (cell->info().isFictious) continue; + capillaryVolume = capillaryVolume + cell->info().capillaryCellVolume; + if (cell->info().p()==0) { + waterVolume = waterVolume + cell->info().capillaryCellVolume;}}} + Real saturation = waterVolume/capillaryVolume; + return saturation; +} + +double UnsaturatedEngine::computeEffPoreRadius(CellHandle cell, int j) +{ + double rInscribe = abs(solver->computeEffectiveRadius(cell, j)); + CellHandle cellh = CellHandle(cell); + int facetNFictious = solver->detectFacetFictiousVertices (cellh,j); + switch (facetNFictious) { + case (0) : { return computeEffPoreRadiusFine(cell,j); }; break; + case (1) : { return rInscribe; }; break; + case (2) : { return rInscribe; }; break; + } +} +//seperate rmin=getMinPoreRadius(cell,j) later; +double UnsaturatedEngine::computeEffPoreRadiusFine(CellHandle cell, int j) +{ + RTriangulation& tri = solver->T[solver->currentTes].Triangulation(); + if (tri.is_infinite(cell->neighbor(j))) return 0; + + Vector3r posA = makeVector3r(cell->vertex(facetVertices[j][0])->point().point()); + Vector3r posB = makeVector3r(cell->vertex(facetVertices[j][1])->point().point()); + Vector3r posC = makeVector3r(cell->vertex(facetVertices[j][2])->point().point()); + double rA = sqrt(cell->vertex(facetVertices[j][0])->point().weight()); + double rB = sqrt(cell->vertex(facetVertices[j][1])->point().weight()); + double rC = sqrt(cell->vertex(facetVertices[j][2])->point().weight()); + double AB = (posA-posB).norm(); + double AC = (posA-posC).norm(); + double BC = (posB-posC).norm(); + double A = acos(((posB-posA).dot(posC-posA))/((posB-posA).norm()*(posC-posA).norm())); + double B = acos(((posA-posB).dot(posC-posB))/((posA-posB).norm()*(posC-posB).norm())); + double C = acos(((posA-posC).dot(posB-posC))/((posA-posC).norm()*(posB-posC).norm())); + double rAB = 0.5*(AB-rA-rB); if (rAB<0) { rAB=0; } + double rBC = 0.5*(BC-rB-rC); if (rBC<0) { rBC=0; } + double rAC = 0.5*(AC-rA-rC); if (rAC<0) { rAC=0; } + + double rmin = max(rAB,max(rBC,rAC)); if (rmin==0) { rmin= 1.0e-10; } + double rmax = abs(solver->computeEffectiveRadius(cell, j));//rmin>rmax ? + if(rmin>rmax) { cerr<<"WARNING! rmin>rmax. rmin="<<rmin<<" ,rmax="<<rmax<<endl; } + + double deltaForceRMin = computeDeltaForce(cell,j,rmin); + double deltaForceRMax = computeDeltaForce(cell,j,rmax); + if(deltaForceRMax<0) { + double effPoreRadius = rmax; + cerr<<"deltaForceRMax Negative. cellID: "<<cell->info().index<<". posA="<<posA<<"; posB="<<posB<<"; posC="<< posC<<"; rA="<< rA<<"; rB="<< rB<<"; rC="<<rC<<endl; + return effPoreRadius;} + else if(deltaForceRMin<0) { + double effPoreRadius = bisection(cell,j,rmin,rmax);// cerr<<"1";//we suppose most cases should be this. + return effPoreRadius;} + else if( (deltaForceRMin>0) && (deltaForceRMin<deltaForceRMax) ) { + double effPoreRadius = rmin;// cerr<<"2"; + return effPoreRadius;} + else if(deltaForceRMin>deltaForceRMax) { + double effPoreRadius = rmax; + cerr<<"WARNING! deltaForceRMin>deltaForceRMax. cellID: "<<cell->info().index<<"; deltaForceRMin="<<deltaForceRMin<<"; deltaForceRMax="<<deltaForceRMax<<". posA="<<posA<<"; posB="<<posB<<"; posC="<< posC<<"; rA="<< rA<<"; rB="<< rB<<"; rC="<<rC<<endl; + return effPoreRadius;} +} + +double UnsaturatedEngine::bisection(CellHandle cell, int j, double a, double b) +{ + double m = 0.5*(a+b); + if (abs(b-a)>abs((solver->computeEffectiveRadius(cell, j)*1.0e-6))) { + if ( computeDeltaForce(cell,j,m) * computeDeltaForce(cell,j,a) < 0 ) { + b = m; + return bisection(cell,j,a,b);} + else { + a = m; + return bisection(cell,j,a,b);}} + else return m; +} + +double UnsaturatedEngine::computeDeltaForce(CellHandle cell,int j, double rCap) +{ + RTriangulation& tri = solver->T[solver->currentTes].Triangulation(); + if (tri.is_infinite(cell->neighbor(j))) return 0; + + Vector3r posA = makeVector3r(cell->vertex(facetVertices[j][0])->point().point()); + Vector3r posB = makeVector3r(cell->vertex(facetVertices[j][1])->point().point()); + Vector3r posC = makeVector3r(cell->vertex(facetVertices[j][2])->point().point()); + double rA = sqrt(cell->vertex(facetVertices[j][0])->point().weight()); + double rB = sqrt(cell->vertex(facetVertices[j][1])->point().weight()); + double rC = sqrt(cell->vertex(facetVertices[j][2])->point().weight()); + double AB = (posA-posB).norm(); + double AC = (posA-posC).norm(); + double BC = (posB-posC).norm(); + double A = acos(((posB-posA).dot(posC-posA))/((posB-posA).norm()*(posC-posA).norm())); + double B = acos(((posA-posB).dot(posC-posB))/((posA-posB).norm()*(posC-posB).norm())); + double C = acos(((posA-posC).dot(posB-posC))/((posA-posC).norm()*(posB-posC).norm())); + double rAB = 0.5*(AB-rA-rB); if (rAB<0) { rAB=0; } + double rBC = 0.5*(BC-rB-rC); if (rBC<0) { rBC=0; } + double rAC = 0.5*(AC-rA-rC); if (rAC<0) { rAC=0; } + + //In triangulation ArB,rCap is the radius of sphere r; + double _AB = (pow((rA+rCap),2)+pow(AB,2)-pow((rB+rCap),2))/(2*(rA+rCap)*AB); if(_AB>1.0) {_AB=1.0;} if(_AB<-1.0) {_AB=-1.0;} + double alphaAB = acos(_AB); + double _BA = (pow((rB+rCap),2)+pow(AB,2)-pow((rA+rCap),2))/(2*(rB+rCap)*AB); if(_BA>1.0) {_BA=1.0;} if(_BA<-1.0) {_BA=-1.0;} + double alphaBA = acos(_BA); + double _ArB = (pow((rA+rCap),2)+pow((rB+rCap),2)-pow(AB,2))/(2*(rA+rCap)*(rB+rCap)); if(_ArB>1.0) {_ArB=1.0;} if(_ArB<-1.0) {_ArB=-1.0;} + double alphaArB = acos(_ArB); + + double lengthLiquidAB = alphaArB*rCap; + double AreaArB = 0.5*(rA+rCap)*(rB+rCap)*sin(alphaArB); + double areaLiquidAB = AreaArB-0.5*alphaAB*pow(rA,2)-0.5*alphaBA*pow(rB,2)-0.5*alphaArB*pow(rCap,2); + + double _AC = (pow((rA+rCap),2)+pow(AC,2)-pow((rC+rCap),2))/(2*(rA+rCap)*AC); if(_AC>1.0) {_AC=1.0;} if(_AC<-1.0) {_AC=-1.0;} + double alphaAC = acos(_AC); + double _CA = (pow((rC+rCap),2)+pow(AC,2)-pow((rA+rCap),2))/(2*(rC+rCap)*AC); if(_CA>1.0) {_CA=1.0;} if(_CA<-1.0) {_CA=-1.0;} + double alphaCA = acos(_CA); + double _ArC = (pow((rA+rCap),2)+pow((rC+rCap),2)-pow(AC,2))/(2*(rA+rCap)*(rC+rCap)); if(_ArC>1.0) {_ArC=1.0;} if(_ArC<-1.0) {_ArC=-1.0;} + double alphaArC = acos(_ArC); + + double lengthLiquidAC = alphaArC*rCap; + double AreaArC = 0.5*(rA+rCap)*(rC+rCap)*sin(alphaArC); + double areaLiquidAC = AreaArC-0.5*alphaAC*pow(rA,2)-0.5*alphaCA*pow(rC,2)-0.5*alphaArC*pow(rCap,2); + + double _BC = (pow((rB+rCap),2)+pow(BC,2)-pow((rC+rCap),2))/(2*(rB+rCap)*BC); if(_BC>1.0) {_BC=1.0;} if(_BC<-1.0) {_BC=-1.0;} + double alphaBC = acos(_BC); + double _CB = (pow((rC+rCap),2)+pow(BC,2)-pow((rB+rCap),2))/(2*(rC+rCap)*BC); if(_CB>1.0) {_CB=1.0;} if(_CB<-1.0) {_CB=-1.0;} + double alphaCB = acos(_CB); + double _BrC = (pow((rB+rCap),2)+pow((rC+rCap),2)-pow(BC,2))/(2*(rB+rCap)*(rC+rCap)); if(_BrC>1.0) {_BrC=1.0;} if(_BrC<-1.0) {_BrC=-1.0;} + double alphaBrC = acos(_BrC); + + double lengthLiquidBC = alphaBrC*rCap; + double AreaBrC = 0.5*(rB+rCap)*(rC+rCap)*sin(alphaBrC); + double areaLiquidBC = AreaBrC-0.5*alphaBC*pow(rB,2)-0.5*alphaCB*pow(rC,2)-0.5*alphaBrC*pow(rCap,2); + + double areaCap = sqrt(cell->info().facetSurfaces[j].squared_length()) * cell->info().facetFluidSurfacesRatio[j]; + double areaPore = areaCap - areaLiquidAB - areaLiquidAC - areaLiquidBC; + + //FIXME:rethink here,areaPore Negative, Flat facets, do nothing ? + if(areaPore<0) {cerr<<"ERROR! areaPore Negative. cellID: "<<cell->info().index<<". posA="<<posA<<"; posB="<<posB<<"; posC="<< posC<<"; rA="<< rA<<"; rB="<< rB<<"; rC="<<rC<<endl;} + double perimeterPore = lengthLiquidAB + lengthLiquidAC + lengthLiquidBC + (A - alphaAB - alphaAC)*rA + (B - alphaBA - alphaBC)*rB + (C - alphaCA - alphaCB)*rC; + if(perimeterPore<0) {cerr<<"ERROR! perimeterPore Negative. cellID: "<<cell->info().index<<". posA="<<posA<<"; posB="<<posB<<"; posC="<< posC<<"; rA="<< rA<<"; rB="<< rB<<"; rC="<<rC<<endl;} + + double deltaForce = perimeterPore - areaPore/rCap;//deltaForce=surfaceTension*(perimeterPore - areaPore/rCap) + return deltaForce; +} + +void UnsaturatedEngine::checkCellsConnection() +{ + RTriangulation& tri = solver->T[solver->currentTes].Triangulation(); + ofstream file; + file.open("cellsConnection.txt"); + file << "cell" << " " << "neighborCells" << endl; + FiniteCellsIterator cellEnd = tri.finite_cells_end(); + for ( FiniteCellsIterator cell = tri.finite_cells_begin(); cell != cellEnd; cell++ ) { + file << cell->info().index << " " <<cell->neighbor(0)->info().index << " " << cell->neighbor(1)->info().index << " " << cell->neighbor(2)->info().index << " " << cell->neighbor(3)->info().index << endl; + } + file.close(); +} + +void UnsaturatedEngine::checkEntryCapillaryPressure() +{ + RTriangulation& tri = solver->T[solver->currentTes].Triangulation(); + ofstream file; + file.open("entryCapillaryPressure.txt"); + file << "#List of Connections \n"; + file << "cell" << " " << "neighborCell" << " " << "entryCapillaryPressure" << endl; + FiniteCellsIterator cellEnd = tri.finite_cells_end(); + for ( FiniteCellsIterator cell = tri.finite_cells_begin(); cell != cellEnd; cell++ ) { + if (tri.is_infinite(cell)) continue; + file << cell->info().index << " " <<cell->neighbor(0)->info().index << " " << surfaceTension/cell->info().poreRadius[0] << endl; + file << cell->info().index << " " <<cell->neighbor(1)->info().index << " " << surfaceTension/cell->info().poreRadius[1] << endl; + file << cell->info().index << " " <<cell->neighbor(2)->info().index << " " << surfaceTension/cell->info().poreRadius[2] << endl; + file << cell->info().index << " " <<cell->neighbor(3)->info().index << " " << surfaceTension/cell->info().poreRadius[3] << endl; + } + file.close(); +} + +void UnsaturatedEngine::checkLatticeNodeY(double y) +{ + RTriangulation& tri = solver->T[solver->currentTes].Triangulation(); + if((y<solver->yMin)||(y>solver->yMax)) { + cerr<<"y is out of range! "<<"pleas set y between "<<solver->yMin<<" and "<<solver->yMax<<endl;} + else { + int N=100;// the default Node number for each slice is 100X100 + ofstream file; + std::ostringstream fileNameStream(".txt"); + fileNameStream << "LatticeNodeY_"<< y; + std::string fileName = fileNameStream.str(); + file.open(fileName.c_str()); +// file << "#Slice Of LatticeNodes: 0: out of sphere; 1: inside of sphere \n"; + Real deltaX = (solver->xMax-solver->xMin)/N; + Real deltaZ = (solver->zMax-solver->zMin)/N; + for (int j=0; j<N+1; j++) { + for (int k=0; k<N+1; k++) { + double x=solver->xMin+j*deltaX; + double z=solver->zMin+k*deltaZ; + int M=0; + Vector3r LatticeNode = Vector3r(x,y,z); + for (FiniteVerticesIterator V_it = tri.finite_vertices_begin(); V_it != tri.finite_vertices_end(); V_it++) { + if(V_it->info().isFictious) continue; + Vector3r SphereCenter = makeVector3r(V_it->point().point()); + if ((LatticeNode-SphereCenter).squaredNorm() < V_it->point().weight()) { + M=1; + break;}} + file << M;} + file << "\n";} + file.close();} +} + +void UnsaturatedEngine::checkBoundingCellsInfo() +{ + ofstream file; + file.open("boundInfo.txt"); + file << "#Checking the boundingCells statement"; + file << "CellID" << " CellPressure" << " isAirReservoir" << " isWaterReservoir" <<endl; + RTriangulation& tri = solver->T[solver->currentTes].Triangulation(); + FiniteCellsIterator cellEnd = tri.finite_cells_end(); + for ( FiniteCellsIterator cell = tri.finite_cells_begin(); cell != cellEnd; cell++ ) { + if (tri.is_infinite(cell)) continue; + if (cell->info().index==0) continue; + if ((cell->info().isFictious==true)&&(cell->info().Pcondition==false)) { + file << cell->info().index <<" "<<cell->info().p()<<" "<<cell->info().isAirReservoir<<" "<<cell->info().isWaterReservoir<<endl; + } + } +} + +void UnsaturatedEngine::checkReservoirInfo(int boundN) +{ + if(solver->boundingCells[boundN].size()==0) { + cerr << "please set corresponding bndCondIsPressure[bound] to be true ."<< endl; + } + else { + if (boundN==2) { + ofstream file; + file.open("waterReservoirBoundInfo.txt"); + file << "#Checking the water reservoir cells statement"; + file << "CellID" << " CellPressure" << " isAirReservoir" << " isWaterReservoir" <<endl; + vector<CellHandle>::iterator it = solver->boundingCells[boundN].begin(); + for ( it ; it != solver->boundingCells[boundN].end(); it++) { + if ((*it)->info().index == 0) continue; + file << (*it)->info().index <<" "<<(*it)->info().p()<<" "<<(*it)->info().isAirReservoir<<" "<<(*it)->info().isWaterReservoir<<endl; + } + file.close(); + } + else if (boundN==3) { + ofstream file; + file.open("airReservoirBoundInfo.txt"); + file << "#Checking the air reservoir cells statement"; + file << "CellID"<<" CellPressure"<<" isAirReservoir"<<" isWaterReservoir"<<endl; + vector<CellHandle>::iterator it = solver->boundingCells[boundN].begin(); + for ( it ; it != solver->boundingCells[boundN].end(); it++) { + if ((*it)->info().index == 0) continue; + file << (*it)->info().index <<" "<<(*it)->info().p()<<" "<<(*it)->info().isAirReservoir<<" "<<(*it)->info().isWaterReservoir<<endl; + } + file.close(); + } + else { + cerr<<"This is not a reservoir boundary. Please set boundN to be 2(waterReservoirBound) or 3(airReservoirBound)."<<endl; + } + } +} + +//----temp function for Vahid Joekar-Niasar's data---- +double UnsaturatedEngine::getRadiusMin(CellHandle cell, int j) +{ + RTriangulation& tri = solver->T[solver->currentTes].Triangulation(); + if (tri.is_infinite(cell->neighbor(j))) return 0; + + Vector3r posA = makeVector3r(cell->vertex(facetVertices[j][0])->point().point()); + Vector3r posB = makeVector3r(cell->vertex(facetVertices[j][1])->point().point()); + Vector3r posC = makeVector3r(cell->vertex(facetVertices[j][2])->point().point()); + double rA = sqrt(cell->vertex(facetVertices[j][0])->point().weight()); + double rB = sqrt(cell->vertex(facetVertices[j][1])->point().weight()); + double rC = sqrt(cell->vertex(facetVertices[j][2])->point().weight()); + double AB = (posA-posB).norm(); + double AC = (posA-posC).norm(); + double BC = (posB-posC).norm(); + double A = acos(((posB-posA).dot(posC-posA))/((posB-posA).norm()*(posC-posA).norm())); + double B = acos(((posA-posB).dot(posC-posB))/((posA-posB).norm()*(posC-posB).norm())); + double C = acos(((posA-posC).dot(posB-posC))/((posA-posC).norm()*(posB-posC).norm())); + double rAB = 0.5*(AB-rA-rB); if (rAB<0) { rAB=0; } + double rBC = 0.5*(BC-rB-rC); if (rBC<0) { rBC=0; } + double rAC = 0.5*(AC-rA-rC); if (rAC<0) { rAC=0; } + + double rmin = max(rAB,max(rBC,rAC)); if (rmin==0) { rmin= 1.0e-10; } + return rmin; +} + +void UnsaturatedEngine::debugTemp() +{ + RTriangulation& tri = solver->T[solver->currentTes].Triangulation(); + ofstream file; + file.open("bugTemp.txt"); + FiniteCellsIterator cellEnd = tri.finite_cells_end(); + for ( FiniteCellsIterator cell = tri.finite_cells_begin(); cell != cellEnd; cell++ ) { + if (tri.is_infinite(cell)) continue; + file << cell->info().index << " " <<cell->info().poreRadius[0]<<" "<<getRadiusMin(cell,0)<<" "<<abs(solver->computeEffectiveRadius(cell, 0))<<endl; + file << cell->info().index << " " <<cell->info().poreRadius[1]<<" "<<getRadiusMin(cell,1)<<" "<<abs(solver->computeEffectiveRadius(cell, 1))<<endl; + file << cell->info().index << " " <<cell->info().poreRadius[2]<<" "<<getRadiusMin(cell,2)<<" "<<abs(solver->computeEffectiveRadius(cell, 2))<<endl; + file << cell->info().index << " " <<cell->info().poreRadius[3]<<" "<<getRadiusMin(cell,3)<<" "<<abs(solver->computeEffectiveRadius(cell, 3))<<endl; + } + file.close(); +}//----------end temp function for Vahid Joekar-Niasar's data (clear later)--------------------- + +//----------temp functions for comparison with experiment----------------------- +void UnsaturatedEngine::initializeCellWindowsID() +{ + RTriangulation& tri = solver->T[solver->currentTes].Triangulation(); + FiniteCellsIterator cellEnd = tri.finite_cells_end(); + for ( FiniteCellsIterator cell = tri.finite_cells_begin(); cell != cellEnd; cell++ ) { + for (int i=1; i<(windowsNo+1); i++) { + if ( (cell->info()[1]>(solver->yMin+(i-1)*(solver->yMax-solver->yMin)/windowsNo) ) && (cell->info()[1] < (solver->yMin+i*(solver->yMax-solver->yMin)/windowsNo)) ) + {cell->info().windowsID=i; break;} + } + } +} + +Real UnsaturatedEngine::getWindowsSaturation(int windowsID) +{ + RTriangulation& tri = solver->T[solver->currentTes].Triangulation(); + FiniteCellsIterator cellEnd = tri.finite_cells_end(); + for ( FiniteCellsIterator cell = tri.finite_cells_begin(); cell != cellEnd; cell++ ) { + if (cell->info().windowsID==0) {cerr<<"Please initialize windowsID"<<endl;break;} + } + if (isPhaseTrapped) { + return getWindowsSaturation1(windowsID); + } + else { + return getWindowsSaturation2(windowsID); + } +} +Real UnsaturatedEngine::getWindowsSaturation1(int i) +{ + updatePressureReservoir(); + RTriangulation& tri = solver->T[solver->currentTes].Triangulation(); + Real capillaryVolume = 0.0; //total capillary volume + Real airVolume = 0.0; //air volume + FiniteCellsIterator cellEnd = tri.finite_cells_end(); + + if (invadeBoundary==true) { + for ( FiniteCellsIterator cell = tri.finite_cells_begin(); cell != cellEnd; cell++ ) { + if (tri.is_infinite(cell)) continue; + if (cell->info().Pcondition) continue;//NOTE:reservoirs cells should not be included in saturation + if (cell->info().windowsID != i) continue; + capillaryVolume = capillaryVolume + cell->info().capillaryCellVolume; + if (cell->info().isAirReservoir==true) { + airVolume = airVolume + cell->info().capillaryCellVolume;}}} + else { + for ( FiniteCellsIterator cell = tri.finite_cells_begin(); cell != cellEnd; cell++ ) { + if (tri.is_infinite(cell)) continue; + if (cell->info().Pcondition) continue; + if (cell->info().isFictious) continue; + if (cell->info().windowsID != i) continue; + capillaryVolume = capillaryVolume + cell->info().capillaryCellVolume; + if (cell->info().isAirReservoir==true) { + airVolume = airVolume + cell->info().capillaryCellVolume;}}} + Real saturation = 1 - airVolume/capillaryVolume; + return saturation; +} +Real UnsaturatedEngine::getWindowsSaturation2(int i) +{ + RTriangulation& tri = solver->T[solver->currentTes].Triangulation(); + Real capillaryVolume = 0.0; + Real waterVolume = 0.0; + FiniteCellsIterator cellEnd = tri.finite_cells_end(); + if (invadeBoundary==true) { + for ( FiniteCellsIterator cell = tri.finite_cells_begin(); cell != cellEnd; cell++ ) { + if (tri.is_infinite(cell)) continue; + if (cell->info().Pcondition) continue; + if (cell->info().windowsID != i) continue; + capillaryVolume = capillaryVolume + cell->info().capillaryCellVolume; + if (cell->info().p()==0) { + waterVolume = waterVolume + cell->info().capillaryCellVolume;}}} + else { + for ( FiniteCellsIterator cell = tri.finite_cells_begin(); cell != cellEnd; cell++ ) { + if (tri.is_infinite(cell)) continue; + if (cell->info().Pcondition) continue; + if (cell->info().isFictious) continue; + if (cell->info().windowsID != i) continue; + capillaryVolume = capillaryVolume + cell->info().capillaryCellVolume; + if (cell->info().p()==0) { + waterVolume = waterVolume + cell->info().capillaryCellVolume;}}} + Real saturation = waterVolume/capillaryVolume; + return saturation; +}//----------end temp functions for comparison with experiment------------------- + +///compute forces +void UnsaturatedEngine::computeSolidLine() +{ + RTriangulation& Tri = solver->T[solver->currentTes].Triangulation(); + FiniteCellsIterator cellEnd = Tri.finite_cells_end(); + for (FiniteCellsIterator cell = Tri.finite_cells_begin(); cell != cellEnd; cell++) { + for(int j=0; j<4; j++) { + solver->lineSolidPore(cell, j); + } + } +} + +void UnsaturatedEngine::computeFacetPoreForcesWithCache(bool onlyCache) +{ + RTriangulation& Tri = solver->T[solver->currentTes].Triangulation(); + FiniteCellsIterator cellEnd = Tri.finite_cells_end(); + 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; + + #ifdef parallel_forces + if (solver->noCache) { + solver->perVertexUnitForce.clear(); solver->perVertexPressure.clear(); +// vector<const CVector*> exf; exf.reserve(20); +// vector<const Real*> exp; exp.reserve(20); + solver->perVertexUnitForce.resize(Tri.number_of_vertices()); + solver->perVertexPressure.resize(Tri.number_of_vertices());} + #endif + + CellHandle neighbourCell; + VertexHandle mirrorVertex; + CVector tempVect; + //FIXME : Ema, be carefull with this (noCache), it needs to be turned true after retriangulation + if (solver->noCache) {for (FlowSolver::VCellIterator cellIt=solver->T[currentTes].cellHandles.begin(); cellIt!=solver->T[currentTes].cellHandles.end(); cellIt++){ +// if (noCache) for (FiniteCellsIterator cell = Tri.finite_cells_begin(); cell != cellEnd; cell++) { + CellHandle& cell = *cellIt; + //reset cache + for (int k=0;k<4;k++) cell->info().unitForceVectors[k]=nullVect; + for (int j=0; j<4; j++) if (!Tri.is_infinite(cell->neighbor(j))) { + neighbourCell = cell->neighbor(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; + 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[solver->boundary(cell->vertex(j)->info().id()).coordinate]); + tempVect=-projSurf*solver->boundary(cell->vertex(j)->info().id()).normal; + cell->vertex(j)->info().forces = cell->vertex(j)->info().forces+tempVect*cell->info().p(); + //define the cached value for later use with cache*p + cell->info().unitForceVectors[j]=cell->info().unitForceVectors[j]+ tempVect; + } + /// Apply weighted forces f_k=sqRad_k/sumSqRad*f + CVector Facet_Unit_Force = -fluidSurfk*cell->info().solidLine[j][3]; + CVector Facet_Force = cell->info().p()*Facet_Unit_Force; + + for (int y=0; y<3;y++) { + cell->vertex(facetVertices[j][y])->info().forces = cell->vertex(facetVertices[j][y])->info().forces + Facet_Force*cell->info().solidLine[j][y]; + //add to cached value + cell->info().unitForceVectors[facetVertices[j][y]]=cell->info().unitForceVectors[facetVertices[j][y]]+Facet_Unit_Force*cell->info().solidLine[j][y]; + //uncomment to get total force / comment to get only pore tension forces + if (!cell->vertex(facetVertices[j][y])->info().isFictious) { + cell->vertex(facetVertices[j][y])->info().forces = cell->vertex(facetVertices[j][y])->info().forces -facetNormal*cell->info().p()*crossSections[j][y]; + //add to cached value + cell->info().unitForceVectors[facetVertices[j][y]]=cell->info().unitForceVectors[facetVertices[j][y]]-facetNormal*crossSections[j][y]; + } + } + #ifdef parallel_forces + solver->perVertexUnitForce[cell->vertex(j)->info().id()].push_back(&(cell->info().unitForceVectors[j])); + solver->perVertexPressure[cell->vertex(j)->info().id()].push_back(&(cell->info().p())); + #endif + } + } + solver->noCache=false;//cache should always be defined after execution of this function + if (onlyCache) return; + } else {//use cached values when triangulation doesn't change + #ifndef parallel_forces + for (FiniteCellsIterator cell = Tri.finite_cells_begin(); cell != cellEnd; cell++) { + for (int yy=0;yy<4;yy++) cell->vertex(yy)->info().forces = cell->vertex(yy)->info().forces + cell->info().unitForceVectors[yy]*cell->info().p();} + + #else + #pragma omp parallel for num_threads(ompThreads) + for (int vn=0; vn<= solver->T[currentTes].maxId; vn++) { + VertexHandle& v = solver->T[currentTes].vertexHandles[vn]; +// for (FiniteVerticesIterator v = Tri.finite_vertices_begin(); v != Tri.finite_vertices_end(); ++v){ + const int& id = v->info().id(); + CVector tf (0,0,0); + int k=0; + for (vector<const Real*>::iterator c = solver->perVertexPressure[id].begin(); c != solver->perVertexPressure[id].end(); c++) + tf = tf + (*(solver->perVertexUnitForce[id][k++]))*(**c); + v->info().forces = tf; + } + #endif + } + if (solver->debugOut) { + 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 (solver->boundary(v->info().id()).flowCondition==1) TotalForce = TotalForce + v->info().forces; } + cout << "TotalForce = "<< TotalForce << endl;} +} + +#endif //UNSATURATED_FLOW +#endif //FLOW_ENGINE \ No newline at end of file
_______________________________________________ Mailing list: https://launchpad.net/~yade-dev Post to : [email protected] Unsubscribe : https://launchpad.net/~yade-dev More help : https://help.launchpad.net/ListHelp
