------------------------------------------------------------ revno: 3419 committer: Chao Yuan <[email protected]> timestamp: Fri 2014-05-09 16:40:53 +0200 message: remove old files removed: pkg/dem/UnsaturatedEngine.cpp pkg/dem/UnsaturatedEngine.hpp
-- lp:yade https://code.launchpad.net/~yade-pkg/yade/git-trunk Your team Yade developers is subscribed to branch lp:yade. To unsubscribe from this branch go to https://code.launchpad.net/~yade-pkg/yade/git-trunk/+edit-subscription
=== removed file 'pkg/dem/UnsaturatedEngine.cpp' --- pkg/dem/UnsaturatedEngine.cpp 2014-04-04 17:12:06 +0000 +++ pkg/dem/UnsaturatedEngine.cpp 1970-01-01 00:00:00 +0000 @@ -1,1579 +0,0 @@ -/************************************************************************* -* 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 YADE_CGAL - -#ifdef FLOW_ENGINE -#include<yade/core/Scene.hpp> -#include<yade/lib/base/Math.hpp> -#include<yade/pkg/dem/TesselationWrapper.hpp> -#include<yade/pkg/common/Sphere.hpp> -#include<yade/pkg/common/Wall.hpp> -#include<yade/pkg/common/Box.hpp> -#include <sys/stat.h> -#include <sys/types.h> -#include <boost/thread.hpp> -#include <boost/date_time.hpp> -#include <boost/date_time/posix_time/posix_time.hpp> - -#include "UnsaturatedEngine.hpp" -#include <CGAL/Plane_3.h> -#include <CGAL/Plane_3.h> - -CREATE_LOGGER ( UnsaturatedEngine ); - -CGT::Vecteur makeCgVect2 ( const Vector3r& yv ) {return CGT::Vecteur ( yv[0],yv[1],yv[2] );} -CGT::Point makeCgPoint2 ( const Vector3r& yv ) {return CGT::Point ( yv[0],yv[1],yv[2] );} -Vector3r makeVector3r2 ( const CGT::Point& yv ) {return Vector3r ( yv[0],yv[1],yv[2] );} -Vector3r makeVector3r2 ( const CGT::Vecteur& yv ) {return Vector3r ( yv[0],yv[1],yv[2] );} - - -UnsaturatedEngine::~UnsaturatedEngine() -{ -} - -const int facetVertices [4][3] = {{1,2,3},{0,2,3},{0,1,3},{0,1,2}}; - -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... - Build_Triangulation(P_zero,solver);//create a triangulation and initialize pressure in the elements, everything will be contained in "solver" - initializeCellIndex(solver);//initialize cell index - initializePoreRadius(solver);//save all pore radii before invade - updateVolumeCapillaryCell(solver);//save capillary volume of all cells, for calculating saturation - computeSolidLine(solver);//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... - Build_Triangulation(P_zero,solver);//create a triangulation and initialize pressure in the elements, everything will be contained in "solver" - initializeCellIndex(solver);//initialize cell index - initializePoreRadius(solver);//save all pore radii before invade - updateVolumeCapillaryCell(solver);//save capillary volume of all cells, for calculating saturation - computeSolidLine(solver);//save cell->info().solidLine[j][y] - solver->noCache = true; - } - ///compute invade - if (pressureForce) { - invade(solver); - } - - ///compute force - if(computeForceActivated){ - computeFacetPoreForcesWithCache(solver); - Vector3r force; - Finite_vertices_iterator vertices_end = solver->T[solver->currentTes].Triangulation().finite_vertices_end(); - for ( Finite_vertices_iterator 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); - } -} -} - -template<class Solver> -void UnsaturatedEngine::invade(Solver& flow) -{ - if (isPhaseTrapped) { - invade1(solver); - } - else { - invade2(solver); - } -} -template<class Solver> -Real UnsaturatedEngine::getMinEntryValue(Solver& flow ) -{ - if (isPhaseTrapped) { - return getMinEntryValue1(solver); - } - else { - return getMinEntryValue2(solver); - } -} -template<class Solver> -Real UnsaturatedEngine::getSaturation(Solver& flow ) -{ - if (isPhaseTrapped) { - return getSaturation1(flow); - } - else { - return getSaturation2(flow); - } -} - -///invade mode 1. update phase reservoir before invasion. Consider no viscous effects, and invade gradually. -template<class Solver> -void UnsaturatedEngine::invadeSingleCell1(Cell_handle cell, double pressure, Solver& flow) -{ - if (invadeBoundary==true) { - for (int facet = 0; facet < 4; facet ++) { - if (flow->T[flow->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) { - Cell_handle nCell = cell->neighbor(facet); - nCell->info().p() = pressure; - nCell->info().isAirReservoir=true; - nCell->info().isWaterReservoir=false; - invadeSingleCell1(nCell, pressure, flow);}}}} - else { - for (int facet = 0; facet < 4; facet ++) { - if (flow->T[flow->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) { - Cell_handle nCell = cell->neighbor(facet); - nCell->info().p() = pressure; - nCell->info().isAirReservoir=true; - nCell->info().isWaterReservoir=false; - invadeSingleCell1(nCell, pressure, flow);}}}} -} - -template<class Solver> -void UnsaturatedEngine::invade1(Solver& flow) -{ - updatePressureReservoir(flow); - RTriangulation& tri = flow->T[flow->currentTes].Triangulation(); - Finite_cells_iterator cell_end = tri.finite_cells_end(); - for ( Finite_cells_iterator cell = tri.finite_cells_begin(); cell != cell_end; cell++ ) { - if(cell->info().isAirReservoir == true) - invadeSingleCell1(cell,cell->info().p(),flow); - } - checkTrap(flow,bndCondValue[3]); - Finite_cells_iterator _cell_end = tri.finite_cells_end(); - for ( Finite_cells_iterator _cell = tri.finite_cells_begin(); _cell != _cell_end; _cell++ ) { - if( (_cell->info().isWaterReservoir) || (_cell->info().isAirReservoir) ) continue; - _cell->info().p() = bndCondValue[3] - _cell->info().trapCapP; - } - initReservoirBound(flow); -} - -template<class Solver>//check trapped phase, define trapCapP. -void UnsaturatedEngine::checkTrap(Solver& flow, double pressure) -{ - RTriangulation& tri = flow->T[flow->currentTes].Triangulation(); - Finite_cells_iterator cell_end = tri.finite_cells_end(); - for ( Finite_cells_iterator cell = tri.finite_cells_begin(); cell != cell_end; cell++ ) { - if( (cell->info().isWaterReservoir) || (cell->info().isAirReservoir) ) continue; - if(cell->info().trapCapP!=0) continue; - cell->info().trapCapP=pressure; - } -} - -template<class Solver> -Real UnsaturatedEngine::getMinEntryValue1(Solver& flow ) -{ - updatePressureReservoir(flow); - Real nextEntry = 1e50; - RTriangulation& tri = flow->T[flow->currentTes].Triangulation(); - if (invadeBoundary==true) { - Finite_cells_iterator cell_end = tri.finite_cells_end(); - for ( Finite_cells_iterator cell = tri.finite_cells_begin(); cell != cell_end; 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 n_cell_pe = surfaceTension/cell->info().poreRadius[facet]; - nextEntry = min(nextEntry,n_cell_pe);}}}}} - else { - Finite_cells_iterator cell_end = tri.finite_cells_end(); - for ( Finite_cells_iterator cell = tri.finite_cells_begin(); cell != cell_end; 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 n_cell_pe = surfaceTension/cell->info().poreRadius[facet]; - nextEntry = min(nextEntry,n_cell_pe);}}}}} - if (nextEntry==1e50) { - cout << "End drainage !" << endl; - return nextEntry=0;} - else return nextEntry; -} - -template<class Solver> -void UnsaturatedEngine::updatePressure(Solver& flow) -{ - BoundaryConditions(flow); - flow->pressureChanged=true; - flow->reApplyBoundaryConditions(); - RTriangulation& tri = flow->T[flow->currentTes].Triangulation(); - Finite_cells_iterator cell_end = tri.finite_cells_end(); - for ( Finite_cells_iterator cell = tri.finite_cells_begin(); cell != cell_end; cell++ ) { - if (cell->info().isWaterReservoir==true) cell->info().p()=bndCondValue[2]; - if (cell->info().isAirReservoir==true) cell->info().p()=bndCondValue[3]; - } -} - -template<class Solver>//update pressure and reservoir attr -void UnsaturatedEngine::updatePressureReservoir(Solver& flow) -{ - updatePressure(flow);//NOTE:updatePressure must be run before update reservoirs. - updateAirReservoir(flow); - updateWaterReservoir(flow); -} - -template<class Solver>//NOTE:keep boundingCells[2],boundingCells[3] always being reservoirs. -void UnsaturatedEngine::initReservoirBound(Solver& flow) -{ - initWaterReservoirBound(flow); - initAirReservoirBound(flow); -} - -template<class Solver>//boundingCells[2] is water reservoir. -void UnsaturatedEngine::initWaterReservoirBound(Solver& flow) -{ - if (flow->boundingCells[2].size()==0) { - cerr<<"ERROR! set bndCondIsPressure[2] true. boundingCells.size=0!"; - } - else { - vector<Cell_handle>::iterator it = flow->boundingCells[2].begin(); - for ( it ; it != flow->boundingCells[2].end(); it++) { - if ((*it)->info().index == 0) continue; - (*it)->info().isWaterReservoir = true; - } - } -} -template<class Solver> -void UnsaturatedEngine::updateWaterReservoir(Solver& flow) -{ - RTriangulation& tri = flow->T[flow->currentTes].Triangulation(); - Finite_cells_iterator cell_end = tri.finite_cells_end(); - for ( Finite_cells_iterator cell = tri.finite_cells_begin(); cell != cell_end; cell++ ) { - cell->info().isWaterReservoir = false; - } - initWaterReservoirBound(flow); - vector<Cell_handle>::iterator it = flow->boundingCells[2].begin(); - for ( it ; it != flow->boundingCells[2].end(); it++) { - if ((*it)->info().index == 0) continue; - waterReservoirRecursion((*it),flow); - } -} -template<class Solver> -void UnsaturatedEngine::waterReservoirRecursion(Cell_handle cell, Solver& flow) -{ - for (int facet = 0; facet < 4; facet ++) { - if (flow->T[flow->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; - Cell_handle n_cell = cell->neighbor(facet); - n_cell->info().isWaterReservoir = true; - waterReservoirRecursion(n_cell,flow); - } -} - -template<class Solver>//boundingCells[3] is air reservoir -void UnsaturatedEngine::initAirReservoirBound(Solver& flow) -{ - if (flow->boundingCells[3].size()==0) { - cerr<<"ERROR! set bndCondIsPressure[3] true. boundingCells.size=0!"; - } - else { - vector<Cell_handle>::iterator it = flow->boundingCells[3].begin(); - for ( it ; it != flow->boundingCells[3].end(); it++) { - if((*it)->info().index == 0) continue; - (*it)->info().isAirReservoir = true; - } - } -} -template<class Solver> -void UnsaturatedEngine::updateAirReservoir(Solver& flow) -{ - RTriangulation& tri = flow->T[flow->currentTes].Triangulation(); - Finite_cells_iterator cell_end = tri.finite_cells_end(); - for ( Finite_cells_iterator cell = tri.finite_cells_begin(); cell != cell_end; cell++ ) { - cell->info().isAirReservoir = false; - } - initAirReservoirBound(flow); - vector<Cell_handle>::iterator it = flow->boundingCells[3].begin(); - for ( it ; it != flow->boundingCells[3].end(); it++) { - if ((*it)->info().index == 0) continue; - airReservoirRecursion((*it),flow); - } -} -template<class Solver> -void UnsaturatedEngine::airReservoirRecursion(Cell_handle cell, Solver& flow) -{ - for (int facet = 0; facet < 4; facet ++) { - if (flow->T[flow->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; - Cell_handle n_cell = cell->neighbor(facet); - n_cell->info().isAirReservoir = true; - airReservoirRecursion(n_cell,flow); - } -} - -template<class Solver> -Real UnsaturatedEngine::getSaturation1 (Solver& flow ) -{ - updatePressureReservoir(flow); - RTriangulation& tri = flow->T[flow->currentTes].Triangulation(); - Real capillary_volume = 0.0; //total capillary volume - Real air_volume = 0.0; //air volume - Finite_cells_iterator cell_end = tri.finite_cells_end(); - - if (invadeBoundary==true) { - for ( Finite_cells_iterator cell = tri.finite_cells_begin(); cell != cell_end; cell++ ) { - if (tri.is_infinite(cell)) continue; - if (cell->info().Pcondition) continue;//NOTE:reservoirs cells should not be included in saturation - capillary_volume = capillary_volume + cell->info().capillaryCellVolume; - if (cell->info().isAirReservoir==true) { - air_volume = air_volume + cell->info().capillaryCellVolume;}}} - else { - for ( Finite_cells_iterator cell = tri.finite_cells_begin(); cell != cell_end; cell++ ) { - if (tri.is_infinite(cell)) continue; - if (cell->info().Pcondition) continue; - if (cell->info().isFictious) continue; - capillary_volume = capillary_volume + cell->info().capillaryCellVolume; - if (cell->info().isAirReservoir==true) { - air_volume = air_volume + cell->info().capillaryCellVolume;}}} - Real saturation = 1 - air_volume/capillary_volume; - return saturation; -} - -///invade mode 2. Consider no trapped phase. -template<class Solver> -void UnsaturatedEngine::invadeSingleCell2(Cell_handle cell, double pressure, Solver& flow) -{ - if (invadeBoundary==true) { - for (int facet = 0; facet < 4; facet ++) { - if (flow->T[flow->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) { - Cell_handle nCell = cell->neighbor(facet); - nCell->info().p() = pressure; - invadeSingleCell2(nCell, pressure, flow);}}} - else { - for (int facet = 0; facet < 4; facet ++) { - if (flow->T[flow->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) { - Cell_handle nCell = cell->neighbor(facet); - nCell->info().p() = pressure; - invadeSingleCell2(nCell, pressure, flow);}}} -} - -template<class Solver> -void UnsaturatedEngine::invade2(Solver& flow ) -{ - updatePressure2(flow); - RTriangulation& tri = flow->T[flow->currentTes].Triangulation(); - Finite_cells_iterator cell_end = tri.finite_cells_end(); - for ( Finite_cells_iterator cell = tri.finite_cells_begin(); cell != cell_end; cell++ ) { - if (cell->info().p()!=0) { - invadeSingleCell2(cell,cell->info().p(),flow); - } - } -} - -template<class Solver> -void UnsaturatedEngine::updatePressure2(Solver& flow) -{ - BoundaryConditions(flow); - flow->pressureChanged=true; - flow->reApplyBoundaryConditions(); - RTriangulation& tri = flow->T[flow->currentTes].Triangulation(); - Finite_cells_iterator cell_end = tri.finite_cells_end(); - for ( Finite_cells_iterator cell = tri.finite_cells_begin(); cell != cell_end; cell++ ) { - if (cell->info().p()!=0) cell->info().p()=bndCondValue[3]; - } -} - -template<class Solver> -Real UnsaturatedEngine::getMinEntryValue2(Solver& flow ) -{ - Real nextEntry = 1e50; - RTriangulation& tri = flow->T[flow->currentTes].Triangulation(); - Finite_cells_iterator cell_end = tri.finite_cells_end(); - - if (invadeBoundary==true) { - for ( Finite_cells_iterator cell = tri.finite_cells_begin(); cell != cell_end; 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 n_cell_pe = surfaceTension/cell->info().poreRadius[facet]; - nextEntry = min(nextEntry,n_cell_pe);}}}}} - else { - for ( Finite_cells_iterator cell = tri.finite_cells_begin(); cell != cell_end; 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 n_cell_pe = surfaceTension/cell->info().poreRadius[facet]; - nextEntry = min(nextEntry,n_cell_pe);}}}}} - if (nextEntry==1e50) { - cout << "End drainage !" << endl; - return nextEntry=0; - } - else { - return nextEntry; - } -} - -template<class Solver> -Real UnsaturatedEngine::getSaturation2(Solver& flow ) -{ - RTriangulation& tri = flow->T[flow->currentTes].Triangulation(); - Real capillary_volume = 0.0; - Real water_volume = 0.0; - Finite_cells_iterator cell_end = tri.finite_cells_end(); - - if (invadeBoundary==true) { - for ( Finite_cells_iterator cell = tri.finite_cells_begin(); cell != cell_end; cell++ ) { - if (tri.is_infinite(cell)) continue; - if (cell->info().Pcondition) continue; - capillary_volume = capillary_volume + cell->info().capillaryCellVolume; - if (cell->info().p()==0) { - water_volume = water_volume + cell->info().capillaryCellVolume;}}} - else { - for ( Finite_cells_iterator cell = tri.finite_cells_begin(); cell != cell_end; cell++ ) { - if (tri.is_infinite(cell)) continue; - if (cell->info().Pcondition) continue; - if (cell->info().isFictious) continue; - capillary_volume = capillary_volume + cell->info().capillaryCellVolume; - if (cell->info().p()==0) { - water_volume = water_volume + cell->info().capillaryCellVolume;}}} - Real saturation = water_volume/capillary_volume; - return saturation; -} - -template<class Cellhandle> -double UnsaturatedEngine::computeEffPoreRadius(Cellhandle cell, int j) -{ - double rInscribe = abs(solver->Compute_EffectiveRadius(cell, j)); - Cell_handle cellh = Cell_handle(cell); - int facetNFictious = solver->detectFacetFictiousVertices (cellh,j); - switch (facetNFictious) { - case (0) : { return computeEffPoreRadiusNormal(cell,j); }; break; - case (1) : { return rInscribe; }; break; - case (2) : { return rInscribe; }; break; - } -} -template<class Cellhandle>//seperate rmin=getMinPoreRadius(cell,j) later; -double UnsaturatedEngine::computeEffPoreRadiusNormal(Cellhandle cell, int j) -{ - RTriangulation& tri = solver->T[solver->currentTes].Triangulation(); - if (tri.is_infinite(cell->neighbor(j))) return 0; - - Vector3r posA = makeVector3r2(cell->vertex(facetVertices[j][0])->point().point()); - Vector3r posB = makeVector3r2(cell->vertex(facetVertices[j][1])->point().point()); - Vector3r posC = makeVector3r2(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->Compute_EffectiveRadius(cell, j));//rmin>rmax ? - if(rmin>rmax) { cerr<<"WARNING! rmin>rmax. rmin="<<rmin<<" ,rmax="<<rmax<<endl; } - - double deltaForce_rMin = computeDeltaForce(cell,j,rmin); - double deltaForce_rMax = computeDeltaForce(cell,j,rmax); - if(deltaForce_rMax<0) { - double EffPoreRadius = rmax; - cerr<<"deltaForce_rMax Negative. cellID: "<<cell->info().index<<". posA="<<posA<<"; posB="<<posB<<"; posC="<< posC<<"; rA="<< rA<<"; rB="<< rB<<"; rC="<<rC<<endl; - return EffPoreRadius;} - else if(deltaForce_rMin<0) { - double effPoreRadius = bisection(cell,j,rmin,rmax);// cerr<<"1";//we suppose most cases should be this. - return effPoreRadius;} - else if( (deltaForce_rMin>0) && (deltaForce_rMin<deltaForce_rMax) ) { - double EffPoreRadius = rmin;// cerr<<"2"; - return EffPoreRadius;} - else if(deltaForce_rMin>deltaForce_rMax) { - double EffPoreRadius = rmax; - cerr<<"WARNING! deltaForce_rMin>deltaForce_rMax. cellID: "<<cell->info().index<<"; deltaForce_rMin="<<deltaForce_rMin<<"; deltaForce_rMax="<<deltaForce_rMax<<". posA="<<posA<<"; posB="<<posB<<"; posC="<< posC<<"; rA="<< rA<<"; rB="<< rB<<"; rC="<<rC<<endl; - return EffPoreRadius;} -} - -template<class Cellhandle> -double UnsaturatedEngine::bisection(Cellhandle cell, int j, double a, double b) -{ - double m = 0.5*(a+b); - if (abs(b-a)>abs((solver->Compute_EffectiveRadius(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; -} - -template<class Cellhandle> -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 = makeVector3r2(cell->vertex(facetVertices[j][0])->point().point()); - Vector3r posB = makeVector3r2(cell->vertex(facetVertices[j][1])->point().point()); - Vector3r posC = makeVector3r2(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 length_liquidAB = alphaArB*rcap; - double AreaArB = 0.5*(rA+rcap)*(rB+rcap)*sin(alphaArB); - double Area_liquidAB = 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 length_liquidAC = alphaArC*rcap; - double AreaArC = 0.5*(rA+rcap)*(rC+rcap)*sin(alphaArC); - double Area_liquidAC = 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 length_liquidBC = alphaBrC*rcap; - double AreaBrC = 0.5*(rB+rcap)*(rC+rcap)*sin(alphaBrC); - double Area_liquidBC = 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 - Area_liquidAB - Area_liquidAC - Area_liquidBC; - - //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 = length_liquidAB + length_liquidAC + length_liquidBC + (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; -} - -template<class Solver> -unsigned int UnsaturatedEngine::imposePressure(Vector3r pos, Real p,Solver& flow) -{ - if (!flow) LOG_ERROR("no flow defined yet, run at least one iter"); - flow->imposedP.push_back( pair<CGT::Point,Real>(CGT::Point(pos[0],pos[1],pos[2]),p) ); - return flow->imposedP.size()-1; -} - -template<class Solver> -void UnsaturatedEngine::BoundaryConditions ( Solver& flow ) -{ - - for (int k=0;k<6;k++) { - flow->boundary (wallIds[k]).flowCondition=!bndCondIsPressure[k]; - flow->boundary (wallIds[k]).value=bndCondValue[k]; - flow->boundary (wallIds[k]).velocity = boundaryVelocity[k];//FIXME: needs correct implementation, maybe update the cached pos/vel? - } -} - -template<class Solver> -void UnsaturatedEngine::initSolver ( Solver& flow ) -{ - flow->Vtotalissimo=0; flow->Vsolid_tot=0; flow->Vporale=0; flow->Ssolid_tot=0; - flow->T[flow->currentTes].Clear(); - flow->T[flow->currentTes].max_id=-1; - flow->x_min = 1000.0, flow->x_max = -10000.0, flow->y_min = 1000.0, flow->y_max = -10000.0, flow->z_min = 1000.0, flow->z_max = -10000.0; -} - -template<class Solver> -void UnsaturatedEngine::Build_Triangulation ( Solver& flow ) -{ - Build_Triangulation ( 0.f,flow ); -} - -template<class Solver> -void UnsaturatedEngine::Build_Triangulation ( double P_zero, Solver& flow ) -{ - flow->ResetNetwork(); - if (first) flow->currentTes=0; - else { - flow->currentTes=!flow->currentTes; - if (Debug) cout << "--------RETRIANGULATION-----------" << endl; - } - - initSolver(flow); - - AddBoundary ( flow ); - Triangulate ( flow ); - if ( Debug ) cout << endl << "Tesselating------" << endl << endl; - flow->T[flow->currentTes].Compute(); - - flow->Define_fictious_cells(); - // For faster loops on cells define this vector - flow->T[flow->currentTes].cellHandles.clear(); - flow->T[flow->currentTes].cellHandles.reserve(flow->T[flow->currentTes].Triangulation().number_of_finite_cells()); - Finite_cells_iterator cell_end = flow->T[flow->currentTes].Triangulation().finite_cells_end(); - for ( Finite_cells_iterator cell = flow->T[flow->currentTes].Triangulation().finite_cells_begin(); cell != cell_end; cell++ ) - flow->T[flow->currentTes].cellHandles.push_back(cell); - flow->DisplayStatistics (); - flow->Compute_Permeability ( ); - - porosity = flow->V_porale_porosity/flow->V_totale_porosity; - - BoundaryConditions ( flow ); - flow->Initialize_pressures ( P_zero ); -} - -void UnsaturatedEngine::setPositionsBuffer(bool current) -{ - vector<posData>& buffer = current? positionBufferCurrent : positionBufferParallel; - buffer.clear(); - buffer.resize(scene->bodies->size()); - shared_ptr<Sphere> sph ( new Sphere ); - const int Sph_Index = sph->getClassIndexStatic(); - FOREACH ( const shared_ptr<Body>& b, *scene->bodies ) { - posData& dat = buffer[b->getId()]; - dat.id=b->getId(); - dat.pos=b->state->pos; - dat.isSphere= (b->shape->getClassIndex() == Sph_Index); - if (dat.isSphere) dat.radius = YADE_CAST<Sphere*>(b->shape.get())->radius; - dat.exists=true; - } -} - -template<class Solver> -void UnsaturatedEngine::AddBoundary ( Solver& flow ) -{ - vector<posData>& buffer = positionBufferCurrent; - solver->x_min = Mathr::MAX_REAL, solver->x_max = -Mathr::MAX_REAL, solver->y_min = Mathr::MAX_REAL, solver->y_max = -Mathr::MAX_REAL, solver->z_min = Mathr::MAX_REAL, solver->z_max = -Mathr::MAX_REAL; - FOREACH ( const posData& b, buffer ) { - if ( !b.exists ) continue; - if ( b.isSphere ) { - const Real& rad = b.radius; - const Real& x = b.pos[0]; - const Real& y = b.pos[1]; - const Real& z = b.pos[2]; - flow->x_min = min ( flow->x_min, x-rad ); - flow->x_max = max ( flow->x_max, x+rad ); - flow->y_min = min ( flow->y_min, y-rad ); - flow->y_max = max ( flow->y_max, y+rad ); - flow->z_min = min ( flow->z_min, z-rad ); - flow->z_max = max ( flow->z_max, z+rad ); - } - } - //FIXME id_offset must be set correctly, not the case here (always 0), then we need walls first or it will fail - id_offset = flow->T[flow->currentTes].max_id+1; - flow->id_offset = id_offset; - flow->SectionArea = ( flow->x_max - flow->x_min ) * ( flow->z_max-flow->z_min ); - flow->Vtotale = ( flow->x_max-flow->x_min ) * ( flow->y_max-flow->y_min ) * ( flow->z_max-flow->z_min ); - flow->y_min_id=wallIds[ymin]; - flow->y_max_id=wallIds[ymax]; - flow->x_max_id=wallIds[xmax]; - flow->x_min_id=wallIds[xmin]; - flow->z_min_id=wallIds[zmin]; - flow->z_max_id=wallIds[zmax]; - - //FIXME: Id's order in boundsIds is done according to the enumeration of boundaries from TXStressController.hpp, line 31. DON'T CHANGE IT! - flow->boundsIds[0]= &flow->x_min_id; - flow->boundsIds[1]= &flow->x_max_id; - flow->boundsIds[2]= &flow->y_min_id; - flow->boundsIds[3]= &flow->y_max_id; - flow->boundsIds[4]= &flow->z_min_id; - flow->boundsIds[5]= &flow->z_max_id; - - for (int k=0;k<6;k++) flow->boundary ( *flow->boundsIds[k] ).useMaxMin = boundaryUseMaxMin[k]; - - flow->Corner_min = CGT::Point ( flow->x_min, flow->y_min, flow->z_min ); - flow->Corner_max = CGT::Point ( flow->x_max, flow->y_max, flow->z_max ); - - //assign BCs types and values - BoundaryConditions ( flow ); - - 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() ); - 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]; -// cerr << "id="<<*flow->boundsIds[i] <<" center="<<center[0]<<","<<center[1]<<","<<center[2]<<endl; - flow->AddBoundingPlane ( center, wall_thickness, Normal,*flow->boundsIds[i] ); - } - } -} - -template<class Solver> -void UnsaturatedEngine::Triangulate ( Solver& flow ) -{ -///Using Tesselation wrapper (faster) -// TesselationWrapper TW; -// if (TW.Tes) delete TW.Tes; -// TW.Tes = &(flow->T[flow->currentTes]);//point to the current Tes we have in Flowengine -// TW.insertSceneSpheres();//TW is now really inserting in UnsaturatedEngine, using the faster insert(begin,end) -// TW.Tes = NULL;//otherwise, Tes would be deleted by ~TesselationWrapper() at the end of the function. -///Using one-by-one insertion - vector<posData>& buffer = positionBufferCurrent; - FOREACH ( const posData& b, buffer ) { - if ( !b.exists ) continue; - if ( b.isSphere ) { -// if (b.id==ignoredBody) continue; - flow->T[flow->currentTes].insert ( b.pos[0], b.pos[1], b.pos[2], b.radius, b.id );} - } - flow->T[flow->currentTes].redirected=true;//By inserting one-by-one, we already redirected -} - -template<class Solver> -void UnsaturatedEngine::Initialize_volumes ( Solver& flow ) -{ - typedef typename Solver::element_type Flow; - typedef typename Flow::Finite_vertices_iterator Finite_vertices_iterator; - typedef typename Solver::element_type Flow; - typedef typename Flow::Finite_cells_iterator Finite_cells_iterator; - - Finite_vertices_iterator vertices_end = flow->T[flow->currentTes].Triangulation().finite_vertices_end(); - CGT::Vecteur Zero(0,0,0); - for (Finite_vertices_iterator V_it = flow->T[flow->currentTes].Triangulation().finite_vertices_begin(); V_it!= vertices_end; V_it++) V_it->info().forces=Zero; - - FOREACH(Cell_handle& cell, flow->T[flow->currentTes].cellHandles) - { - switch ( cell->info().fictious() ) - { - case ( 0 ) : cell->info().volume() = Volume_cell ( cell ); break; - case ( 1 ) : cell->info().volume() = Volume_cell_single_fictious ( cell ); break; - case ( 2 ) : cell->info().volume() = Volume_cell_double_fictious ( cell ); break; - case ( 3 ) : cell->info().volume() = Volume_cell_triple_fictious ( cell ); break; - default: break; - } - - if (flow->fluidBulkModulus>0) { cell->info().invVoidVolume() = 1 / ( abs(cell->info().volume()) - flow->volumeSolidPore(cell) ); } - } - if (Debug) cout << "Volumes initialised." << endl; -} - -template<class Solver> -void UnsaturatedEngine::initializeCellIndex(Solver& flow) -{ - int k=0; - RTriangulation& tri = flow->T[flow->currentTes].Triangulation(); - Finite_cells_iterator cell_end = tri.finite_cells_end(); - for ( Finite_cells_iterator cell = tri.finite_cells_begin(); cell != cell_end; cell++ ) { - cell->info().index=k++; - } -} - -template<class Solver> -void UnsaturatedEngine::initializePoreRadius(Solver& flow) -{ - RTriangulation& tri = flow->T[flow->currentTes].Triangulation(); - Finite_cells_iterator cell_end = tri.finite_cells_end(); - Cell_handle neighbour_cell; - for (Finite_cells_iterator cell = tri.finite_cells_begin(); cell != cell_end; cell++) { - for (int j=0; j<4; j++) { - neighbour_cell = cell->neighbor(j); - if (!tri.is_infinite(neighbour_cell)) { - cell->info().poreRadius[j]=computeEffPoreRadius(cell, j); - neighbour_cell->info().poreRadius[tri.mirror_index(cell, j)]= cell->info().poreRadius[j]; - } - } - } -} - -template<class Cellhandle> -Real UnsaturatedEngine::Volume_cell_single_fictious ( Cellhandle cell ) -{ - //With buffer - Vector3r V[3]; - int b=0; - int w=0; - cell->info().volumeSign=1; - Real Wall_coordinate=0; - - for ( int y=0;y<4;y++ ) { - if ( ! ( cell->vertex ( y )->info().isFictious ) ) { -// const shared_ptr<Body>& sph = Body::byId ( cell->vertex ( y )->info().id(), scene ); - V[w]=positionBufferCurrent[cell->vertex ( y )->info().id()].pos; - w++; - } else { - b = cell->vertex ( y )->info().id(); - const shared_ptr<Body>& wll = Body::byId ( b , scene ); - if ( !solver->boundary ( b ).useMaxMin ) Wall_coordinate = wll->state->pos[solver->boundary ( b ).coordinate]+ ( solver->boundary ( b ).normal[solver->boundary ( b ).coordinate] ) *wall_thickness/2; - else Wall_coordinate = solver->boundary ( b ).p[solver->boundary ( b ).coordinate]; - } - } - Real Volume = 0.5* ( ( V[0]-V[1] ).cross ( V[0]-V[2] ) ) [solver->boundary ( b ).coordinate] * ( 0.33333333333* ( V[0][solver->boundary ( b ).coordinate]+ V[1][solver->boundary ( b ).coordinate]+ V[2][solver->boundary ( b ).coordinate] ) - Wall_coordinate ); - return abs ( Volume ); -} -template<class Cellhandle> -Real UnsaturatedEngine::Volume_cell_double_fictious ( Cellhandle cell ) -{ - Vector3r A=Vector3r::Zero(), AS=Vector3r::Zero(),B=Vector3r::Zero(), BS=Vector3r::Zero(); - - cell->info().volumeSign=1; - int b[2]; - int coord[2]; - Real Wall_coordinate[2]; - int j=0; - bool first_sph=true; - - for ( int g=0;g<4;g++ ) { - if ( cell->vertex ( g )->info().isFictious ) { - b[j] = cell->vertex ( g )->info().id(); - coord[j]=solver->boundary ( b[j] ).coordinate; - if ( !solver->boundary ( b[j] ).useMaxMin ) Wall_coordinate[j] = positionBufferCurrent[b[j]].pos[coord[j]] + ( solver->boundary ( b[j] ).normal[coord[j]] ) *wall_thickness/2; - else Wall_coordinate[j] = solver->boundary ( b[j] ).p[coord[j]]; - j++; - } else if ( first_sph ) { - A=AS=/*AT=*/ positionBufferCurrent[cell->vertex(g)->info().id()].pos; - first_sph=false; - } else { - B=BS=/*BT=*/ positionBufferCurrent[cell->vertex(g)->info().id()].pos;; - } - } - AS[coord[0]]=BS[coord[0]] = Wall_coordinate[0]; - - //first pyramid with triangular base (A,B,BS) - Real Vol1=0.5* ( ( A-BS ).cross ( B-BS ) ) [coord[1]]* ( 0.333333333* ( 2*B[coord[1]]+A[coord[1]] )-Wall_coordinate[1] ); - //second pyramid with triangular base (A,AS,BS) - Real Vol2=0.5* ( ( AS-BS ).cross ( A-BS ) ) [coord[1]]* ( 0.333333333* ( B[coord[1]]+2*A[coord[1]] )-Wall_coordinate[1] ); - return abs ( Vol1+Vol2 ); -} -template<class Cellhandle> -Real UnsaturatedEngine::Volume_cell_triple_fictious ( Cellhandle cell ) -{ - Vector3r A; - - int b[3]; - int coord[3]; - Real Wall_coordinate[3]; - int j=0; - cell->info().volumeSign=1; - - for ( int g=0;g<4;g++ ) { - if ( cell->vertex ( g )->info().isFictious ) { - b[j] = cell->vertex ( g )->info().id(); - coord[j]=solver->boundary ( b[j] ).coordinate; - const shared_ptr<Body>& wll = Body::byId ( b[j] , scene ); - if ( !solver->boundary ( b[j] ).useMaxMin ) Wall_coordinate[j] = wll->state->pos[coord[j]] + ( solver->boundary ( b[j] ).normal[coord[j]] ) *wall_thickness/2; - else Wall_coordinate[j] = solver->boundary ( b[j] ).p[coord[j]]; - j++; - } else { - const shared_ptr<Body>& sph = Body::byId ( cell->vertex ( g )->info().id(), scene ); - A= ( sph->state->pos ); - } - } - Real Volume = ( A[coord[0]]-Wall_coordinate[0] ) * ( A[coord[1]]-Wall_coordinate[1] ) * ( A[coord[2]]-Wall_coordinate[2] ); - return abs ( Volume ); -} -template<class Cellhandle> -Real UnsaturatedEngine::Volume_cell ( Cellhandle cell ) -{ - static const Real inv6 = 1/6.; - const Vector3r& p0 = positionBufferCurrent[cell->vertex ( 0 )->info().id()].pos; - const Vector3r& p1 = positionBufferCurrent[cell->vertex ( 1 )->info().id()].pos; - const Vector3r& p2 = positionBufferCurrent[cell->vertex ( 2 )->info().id()].pos; - const Vector3r& p3 = positionBufferCurrent[cell->vertex ( 3 )->info().id()].pos; - Real volume = inv6 * ((p1-p0).cross(p2-p0)).dot(p3-p0); - if ( ! ( cell->info().volumeSign ) ) cell->info().volumeSign= ( volume>0 ) ?1:-1; - return volume; -} - -template<class Solver> -void UnsaturatedEngine::updateVolumeCapillaryCell ( Solver& flow) -{ - Initialize_volumes(flow); - RTriangulation& tri = flow->T[flow->currentTes].Triangulation(); - Finite_cells_iterator cell_end = tri.finite_cells_end(); - Cell_handle neighbour_cell; - for (Finite_cells_iterator cell = tri.finite_cells_begin(); cell != cell_end; cell++) { - cell->info().capillaryCellVolume = abs( cell->info().volume() ) - solver->volumeSolidPore(cell); -// if (cell->info().capillaryCellVolume<0) {cerr<<"volumeCapillaryCell Negative. cell ID: " << cell->info().index << "cell volume: " << cell->info().volume() << " volumeSolidPore: " << solver->volumeSolidPore(cell) << endl; } - } -} - -template<class Solver> -void UnsaturatedEngine::saveListNodes(Solver& flow) -{ - RTriangulation& tri = flow->T[flow->currentTes].Triangulation(); - ofstream file; - file.open("ListOfNodes.txt"); - file << "#List Of Nodes. For one cell,there are four neighbour cells \n"; - file << "Cell_ID" << " " << "NeighborCell_ID" << endl; - Finite_cells_iterator cell_end = tri.finite_cells_end(); - for ( Finite_cells_iterator cell = tri.finite_cells_begin(); cell != cell_end; 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(); -} - -template<class Solver> -void UnsaturatedEngine::saveListConnection(Solver& flow) -{ - RTriangulation& tri = flow->T[flow->currentTes].Triangulation(); - ofstream file; - file.open("ListConnection.txt"); - file << "#List of Connections \n"; - file << "Cell_ID" << " " << "NeighborCell_ID" << " " << "EntryValue" << " " << "poreRadius" << " " << "poreArea" << " " << "porePerimeter" << endl; - double surface_tension = surfaceTension ; //Surface Tension in contact with air at 20 Degrees Celsius is:0.0728(N/m) - Finite_cells_iterator cell_end = tri.finite_cells_end(); - for ( Finite_cells_iterator cell = tri.finite_cells_begin(); cell != cell_end; cell++ ) { - if (tri.is_infinite(cell)) continue; - file << cell->info().index << " " <<cell->neighbor(0)->info().index << " " << surface_tension/cell->info().poreRadius[0] << " " << cell->info().poreRadius[0] << " " << computePoreArea(cell,0) << " " << computePorePerimeter(cell,0) << endl; - file << cell->info().index << " " <<cell->neighbor(1)->info().index << " " << surface_tension/cell->info().poreRadius[1] << " " << cell->info().poreRadius[1] << " " << computePoreArea(cell,1) << " " << computePorePerimeter(cell,1) << endl; - file << cell->info().index << " " <<cell->neighbor(2)->info().index << " " << surface_tension/cell->info().poreRadius[2] << " " << cell->info().poreRadius[2] << " " << computePoreArea(cell,2) << " " << computePorePerimeter(cell,2) << endl; - file << cell->info().index << " " <<cell->neighbor(3)->info().index << " " << surface_tension/cell->info().poreRadius[3] << " " << cell->info().poreRadius[3] << " " << computePoreArea(cell,3) << " " << computePorePerimeter(cell,3) << endl; - } - file.close(); -} - -template<class Solver> -void UnsaturatedEngine::saveLatticeNodeX(Solver& flow, double x) -{ - RTriangulation& tri = flow->T[solver->currentTes].Triangulation(); - if((x<flow->x_min)||(x>flow->x_max)) { - cerr<<"x is out of range! "<<"pleas set x between "<<flow->x_min<<" and "<<flow->x_max<<endl; - } - else { - int N=100;// the default Node number for each slice is 100X100 - ofstream file; - std::ostringstream fileNameStream(".txt"); - fileNameStream << "LatticeNodeX_"<< x; - std::string fileName = fileNameStream.str(); - file.open(fileName.c_str()); -// file << "#Slice Of LatticeNodes: 0: out of sphere; 1: inside of sphere \n"; - Real delta_y = (flow->y_max-flow->y_min)/N; - Real delta_z = (flow->z_max-flow->z_min)/N; - for (int j=0; j<N+1; j++) { - for (int k=0; k<N+1; k++) { - double y=flow->y_min+j*delta_y; - double z=flow->z_min+k*delta_z; - int M=0; - Vector3r LatticeNode = Vector3r(x,y,z); - for (Finite_vertices_iterator V_it = tri.finite_vertices_begin(); V_it != tri.finite_vertices_end(); V_it++) { - if(V_it->info().isFictious) continue; - Vector3r SphereCenter = makeVector3r2(V_it->point().point()); - if ((LatticeNode-SphereCenter).squaredNorm() < V_it->point().weight()) { - M=1; - break; - } - } - file << M; - } - file << "\n"; - } - file.close(); - } -} - -template<class Solver> -void UnsaturatedEngine::saveLatticeNodeY(Solver& flow, double y) -{ - RTriangulation& tri = flow->T[solver->currentTes].Triangulation(); - if((y<flow->y_min)||(y>flow->y_max)) { - cerr<<"y is out of range! "<<"pleas set y between "<<flow->y_min<<" and "<<flow->y_max<<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 delta_x = (flow->x_max-flow->x_min)/N; - Real delta_z = (flow->z_max-flow->z_min)/N; - for (int j=0; j<N+1; j++) { - for (int k=0; k<N+1; k++) { - double x=flow->x_min+j*delta_x; - double z=flow->z_min+k*delta_z; - int M=0; - Vector3r LatticeNode = Vector3r(x,y,z); - for (Finite_vertices_iterator V_it = tri.finite_vertices_begin(); V_it != tri.finite_vertices_end(); V_it++) { - if(V_it->info().isFictious) continue; - Vector3r SphereCenter = makeVector3r2(V_it->point().point()); - if ((LatticeNode-SphereCenter).squaredNorm() < V_it->point().weight()) { - M=1; - break; - } - } - file << M; - } - file << "\n"; - } - file.close(); - } -} - -template<class Solver> -void UnsaturatedEngine::saveLatticeNodeZ(Solver& flow, double z) -{ - RTriangulation& tri = flow->T[solver->currentTes].Triangulation(); - if((z<flow->z_min)||(z>flow->z_max)) { - cerr<<"z is out of range! "<<"pleas set z between "<<flow->z_min<<" and "<<flow->z_max<<endl; - } - else { - int N=100;// the default Node number for each slice is 100X100 - ofstream file; - std::ostringstream fileNameStream(".txt"); - fileNameStream << "LatticeNodeZ_"<< z; - std::string fileName = fileNameStream.str(); - file.open(fileName.c_str()); -// file << "#Slice Of LatticeNodes: 0: out of sphere; 1: inside of sphere \n"; - Real delta_x = (flow->x_max-flow->x_min)/N; - Real delta_y = (flow->y_max-flow->y_min)/N; - for (int j=0; j<N+1; j++) { - for (int k=0; k<N+1; k++) { - double x=flow->x_min+j*delta_x; - double y=flow->z_min+k*delta_y; - int M=0; - Vector3r LatticeNode = Vector3r(x,y,z); - for (Finite_vertices_iterator V_it = tri.finite_vertices_begin(); V_it != tri.finite_vertices_end(); V_it++) { - if(V_it->info().isFictious) continue; - Vector3r SphereCenter = makeVector3r2(V_it->point().point()); - if ((LatticeNode-SphereCenter).squaredNorm() < V_it->point().weight()) { - M=1; - break; - } - } - file << M; - } - file << "\n"; - } - file.close(); - } -} - -template<class Solver> -void UnsaturatedEngine::saveBoundingCellsInfo(Solver& flow) -{ - ofstream file; - file.open("boundInfo.txt"); - file << "#Checking the boundingCells statement"; - file << "Cell_ID"<<" Cell_Pressure"<<" isAirReservoir"<<" isWaterReservoir"<<endl; - RTriangulation& tri = flow->T[solver->currentTes].Triangulation(); - Finite_cells_iterator cell_end = tri.finite_cells_end(); - for ( Finite_cells_iterator cell = tri.finite_cells_begin(); cell != cell_end; 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; - } - } -} -template<class Solver> -void UnsaturatedEngine::saveReservoirInfo(Solver& flow,int boundN) -{ - if(flow->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 << "Cell_ID"<<" Cell_Pressure"<<" isAirReservoir"<<" isWaterReservoir"<<endl; - vector<Cell_handle>::iterator it = flow->boundingCells[boundN].begin(); - for ( it ; it != flow->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 << "Cell_ID"<<" Cell_Pressure"<<" isAirReservoir"<<" isWaterReservoir"<<endl; - vector<Cell_handle>::iterator it = flow->boundingCells[boundN].begin(); - for ( it ; it != flow->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; - } - } -} - -template<class Solver> -void UnsaturatedEngine::setImposedPressure ( unsigned int cond, Real p,Solver& flow ) -{ - if ( cond>=flow->imposedP.size() ) LOG_ERROR ( "Setting p with cond higher than imposedP size." ); - flow->imposedP[cond].second=p; - //force immediate update of boundary conditions - flow->pressureChanged=true; -} - -template<class Solver> -void UnsaturatedEngine::clearImposedPressure ( Solver& flow ) { flow->imposedP.clear(); flow->IPCells.clear();} - -//----temp function for Vahid Joekar-Niasar's data---- -template<class Cellhandle > -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 = makeVector3r2(cell->vertex(facetVertices[j][0])->point().point()); - Vector3r posB = makeVector3r2(cell->vertex(facetVertices[j][1])->point().point()); - Vector3r posC = makeVector3r2(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; -} -template<class Solver> -void UnsaturatedEngine::debugTemp(Solver& flow) -{ - RTriangulation& tri = flow->T[solver->currentTes].Triangulation(); - ofstream file; - file.open("bugTemp.txt"); - Finite_cells_iterator cell_end = tri.finite_cells_end(); - for ( Finite_cells_iterator cell = tri.finite_cells_begin(); cell != cell_end; cell++ ) { - if (tri.is_infinite(cell)) continue; - file << cell->info().index << " " <<cell->info().poreRadius[0]<<" "<<getRadiusMin(cell,0)<<" "<<abs(flow->Compute_EffectiveRadius(cell, 0))<<endl; - file << cell->info().index << " " <<cell->info().poreRadius[1]<<" "<<getRadiusMin(cell,1)<<" "<<abs(flow->Compute_EffectiveRadius(cell, 1))<<endl; - file << cell->info().index << " " <<cell->info().poreRadius[2]<<" "<<getRadiusMin(cell,2)<<" "<<abs(flow->Compute_EffectiveRadius(cell, 2))<<endl; - file << cell->info().index << " " <<cell->info().poreRadius[3]<<" "<<getRadiusMin(cell,3)<<" "<<abs(flow->Compute_EffectiveRadius(cell, 3))<<endl; - } - file.close(); -} - -template<class Cellhandle> -Real UnsaturatedEngine::computePoreArea(Cellhandle cell, int j) -{ - double rInscribe = abs(solver->Compute_EffectiveRadius(cell, j)); - Cell_handle cellh = Cell_handle(cell); - int facetNFictious = solver->detectFacetFictiousVertices (cellh,j); - switch (facetNFictious) { - case (0) : { - RTriangulation& tri = solver->T[solver->currentTes].Triangulation(); - if (tri.is_infinite(cell->neighbor(j))) return 0; - - Vector3r posA = makeVector3r2(cell->vertex(facetVertices[j][0])->point().point()); - Vector3r posB = makeVector3r2(cell->vertex(facetVertices[j][1])->point().point()); - Vector3r posC = makeVector3r2(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 rcap = cell->info().poreRadius[j]; - //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 D1=alphaAB + alphaBA + alphaArB; cerr<<D<<" "; - double length_liquidAB = alphaArB*rcap; - double AreaArB = 0.5*(rA+rcap)*(rB+rcap)*sin(alphaArB); - double Area_liquidAB = 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 D2=alphaAC + alphaCA + alphaArC; cerr<<D<<" "; - double length_liquidAC = alphaArC*rcap; - double AreaArC = 0.5*(rA+rcap)*(rC+rcap)*sin(alphaArC); - double Area_liquidAC = 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 D3=alphaBC + alphaCB + alphaBrC; cerr<<D<<" "; - double length_liquidBC = alphaBrC*rcap; - double AreaBrC = 0.5*(rB+rcap)*(rC+rcap)*sin(alphaBrC); - double Area_liquidBC = 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 - Area_liquidAB - Area_liquidAC - Area_liquidBC; - if(areaPore<0) {cerr<<"areaPore Negative. cellID: "<<cell->info().index<<". posA="<<posA<<"; posB="<<posB<<"; posC="<< posC<<"; rA="<< rA<<"; rB="<< rB<<"; rC="<<rC<<endl; - areaPore=Mathr::PI*pow(rInscribe,2);} - return areaPore; - }; break; - case (1) : { return Mathr::PI*pow(rInscribe,2); }; break; - case (2) : { return Mathr::PI*pow(rInscribe,2); }; break; - } -} - -template<class Cellhandle> -Real UnsaturatedEngine::computePorePerimeter(Cellhandle cell, int j) -{ - double rInscribe = abs(solver->Compute_EffectiveRadius(cell, j)); - Cell_handle cellh = Cell_handle(cell); - int facetNFictious = solver->detectFacetFictiousVertices (cellh,j); - switch (facetNFictious) { - case (0) : { - RTriangulation& tri = solver->T[solver->currentTes].Triangulation(); - if (tri.is_infinite(cell->neighbor(j))) return 0; - - Vector3r posA = makeVector3r2(cell->vertex(facetVertices[j][0])->point().point()); - Vector3r posB = makeVector3r2(cell->vertex(facetVertices[j][1])->point().point()); - Vector3r posC = makeVector3r2(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 rcap = cell->info().poreRadius[j]; - //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 D1=alphaAB + alphaBA + alphaArB; cerr<<D<<" "; - double length_liquidAB = alphaArB*rcap; - double AreaArB = 0.5*(rA+rcap)*(rB+rcap)*sin(alphaArB); - double Area_liquidAB = 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 D2=alphaAC + alphaCA + alphaArC; cerr<<D<<" "; - double length_liquidAC = alphaArC*rcap; - double AreaArC = 0.5*(rA+rcap)*(rC+rcap)*sin(alphaArC); - double Area_liquidAC = 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 D3=alphaBC + alphaCB + alphaBrC; cerr<<D<<" "; - double length_liquidBC = alphaBrC*rcap; - double AreaBrC = 0.5*(rB+rcap)*(rC+rcap)*sin(alphaBrC); - double Area_liquidBC = AreaBrC-0.5*alphaBC*pow(rB,2)-0.5*alphaCB*pow(rC,2)-0.5*alphaBrC*pow(rcap,2); - double perimeterPore = length_liquidAB + length_liquidAC + length_liquidBC + (A - alphaAB - alphaAC)*rA + (B - alphaBA - alphaBC)*rB + (C - alphaCA - alphaCB)*rC; - if(perimeterPore<0) {cerr<<"perimeterPore Negative. perimeterPore="<<perimeterPore<<endl;perimeterPore=2.0*Mathr::PI*rInscribe;} - return perimeterPore; - }; break; - case (1) : { return 2.0*Mathr::PI*rInscribe; }; break; - case (2) : { return 2.0*Mathr::PI*rInscribe; }; break; - } -} - -template<class Solver> -void UnsaturatedEngine::savePoreBodyInfo(Solver& flow) -{ - RTriangulation& tri = flow->T[solver->currentTes].Triangulation(); - ofstream file; - file.open("PoreBodyInfo.txt"); - file << "#pore bodies position (or Voronoi centers) and size (volume) \n"; - file << "Cell_ID " << " x " << " y " << " z " << " volume "<< endl; - Finite_cells_iterator cell_end = tri.finite_cells_end(); - for ( Finite_cells_iterator cell = tri.finite_cells_begin(); cell != cell_end; cell++ ) { - if (tri.is_infinite(cell)) continue; - file << cell->info().index << " " << cell->info() << " " << cell->info().capillaryCellVolume << endl; - } - file.close(); -} - -template<class Solver> -void UnsaturatedEngine::savePoreThroatInfo(Solver& flow) -{ - RTriangulation& tri = flow->T[solver->currentTes].Triangulation(); - ofstream file; - file.open("PoreThroatInfo.txt"); - file << "#pore throat area, inscribeRadius and perimeter. \n"; - file << "Cell_ID " << " NeighborCell_ID "<<" area " << " inscribeRadius " << " perimeter " << endl; - Finite_cells_iterator cell_end = tri.finite_cells_end(); - for ( Finite_cells_iterator cell = tri.finite_cells_begin(); cell != cell_end; cell++ ) { - if (tri.is_infinite(cell)) continue; - file << cell->info().index <<" "<< cell->neighbor(0)->info().index <<" "<< computePoreArea(cell, 0) <<" "<< abs(flow->Compute_EffectiveRadius(cell, 0)) <<" "<< computePorePerimeter(cell,0) <<endl; - file << cell->info().index <<" "<< cell->neighbor(1)->info().index <<" "<< computePoreArea(cell, 1) <<" "<< abs(flow->Compute_EffectiveRadius(cell, 1)) <<" "<< computePorePerimeter(cell,1) <<endl; - file << cell->info().index <<" "<< cell->neighbor(2)->info().index <<" "<< computePoreArea(cell, 2) <<" "<< abs(flow->Compute_EffectiveRadius(cell, 2)) <<" "<< computePorePerimeter(cell,2) <<endl; - file << cell->info().index <<" "<< cell->neighbor(3)->info().index <<" "<< computePoreArea(cell, 3) <<" "<< abs(flow->Compute_EffectiveRadius(cell, 3)) <<" "<< computePorePerimeter(cell,3) <<endl; - } - file.close(); -} - -//----------end temp function for Vahid Joekar-Niasar's data (clear later)--------------------- -//----------temp functions for comparison with experiment----------------------- -template <class Solver> -void UnsaturatedEngine::initializeCellWindowsID(Solver&flow) -{ - RTriangulation& tri = flow->T[flow->currentTes].Triangulation(); - Finite_cells_iterator cell_end = tri.finite_cells_end(); - for ( Finite_cells_iterator cell = tri.finite_cells_begin(); cell != cell_end; cell++ ) { - for (int i=1; i<(windowsNo+1); i++) { - if ( (cell->info()[1]>(flow->y_min+(i-1)*(flow->y_max-flow->y_min)/windowsNo) ) && (cell->info()[1] < (flow->y_min+i*(flow->y_max-flow->y_min)/windowsNo)) ) - {cell->info().windowsID=i; break;} - } - } -} -template<class Solver> -Real UnsaturatedEngine::getWindowsSaturation(Solver&flow, int windowsID) -{ - RTriangulation& tri = flow->T[flow->currentTes].Triangulation(); - Finite_cells_iterator cell_end = tri.finite_cells_end(); - for ( Finite_cells_iterator cell = tri.finite_cells_begin(); cell != cell_end; cell++ ) { - if (cell->info().windowsID==0) {cerr<<"Please initialize windowsID"<<endl;break;} - } - if (isPhaseTrapped) { - return getWindowsSaturation1(flow,windowsID); - } - else { - return getWindowsSaturation2(flow,windowsID); - } -} -template<class Solver> -Real UnsaturatedEngine::getWindowsSaturation1(Solver&flow, int i) -{ - updatePressureReservoir(flow); - RTriangulation& tri = flow->T[flow->currentTes].Triangulation(); - Real capillary_volume = 0.0; //total capillary volume - Real air_volume = 0.0; //air volume - Finite_cells_iterator cell_end = tri.finite_cells_end(); - - if (invadeBoundary==true) { - for ( Finite_cells_iterator cell = tri.finite_cells_begin(); cell != cell_end; 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; - capillary_volume = capillary_volume + cell->info().capillaryCellVolume; - if (cell->info().isAirReservoir==true) { - air_volume = air_volume + cell->info().capillaryCellVolume;}}} - else { - for ( Finite_cells_iterator cell = tri.finite_cells_begin(); cell != cell_end; cell++ ) { - if (tri.is_infinite(cell)) continue; - if (cell->info().Pcondition) continue; - if (cell->info().isFictious) continue; - if (cell->info().windowsID != i) continue; - capillary_volume = capillary_volume + cell->info().capillaryCellVolume; - if (cell->info().isAirReservoir==true) { - air_volume = air_volume + cell->info().capillaryCellVolume;}}} - Real saturation = 1 - air_volume/capillary_volume; - return saturation; -} -template<class Solver> -Real UnsaturatedEngine::getWindowsSaturation2(Solver& flow,int i) -{ - RTriangulation& tri = flow->T[flow->currentTes].Triangulation(); - Real capillary_volume = 0.0; - Real water_volume = 0.0; - Finite_cells_iterator cell_end = tri.finite_cells_end(); - if (invadeBoundary==true) { - for ( Finite_cells_iterator cell = tri.finite_cells_begin(); cell != cell_end; cell++ ) { - if (tri.is_infinite(cell)) continue; - if (cell->info().Pcondition) continue; - if (cell->info().windowsID != i) continue; - capillary_volume = capillary_volume + cell->info().capillaryCellVolume; - if (cell->info().p()==0) { - water_volume = water_volume + cell->info().capillaryCellVolume;}}} - else { - for ( Finite_cells_iterator cell = tri.finite_cells_begin(); cell != cell_end; cell++ ) { - if (tri.is_infinite(cell)) continue; - if (cell->info().Pcondition) continue; - if (cell->info().isFictious) continue; - if (cell->info().windowsID != i) continue; - capillary_volume = capillary_volume + cell->info().capillaryCellVolume; - if (cell->info().p()==0) { - water_volume = water_volume + cell->info().capillaryCellVolume;}}} - Real saturation = water_volume/capillary_volume; - return saturation; -} -//----------end temp functions for comparison with experiment------------------- - -template <class Solver> -void UnsaturatedEngine::computeFacetPoreForcesWithCache(Solver& flow, bool onlyCache) -{ - RTriangulation& Tri = flow->T[solver->currentTes].Triangulation(); - Finite_cells_iterator cell_end = Tri.finite_cells_end(); - CGT::Vecteur nullVect(0,0,0); - //reset forces - if (!onlyCache) for (Finite_vertices_iterator v = Tri.finite_vertices_begin(); v != Tri.finite_vertices_end(); ++v) v->info().forces=nullVect; - - #ifdef parallel_forces - if (solver->noCache) { - solver->perVertexUnitForce.clear(); solver->perVertexPressure.clear(); -// vector<const Vecteur*> 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 - - Cell_handle neighbour_cell; - Vertex_handle mirror_vertex; - CGT::Vecteur tempVect; - //FIXME : Ema, be carefull with this (noCache), it needs to be turned true after retriangulation - if (solver->noCache) {for (VCell_iterator cell_it=flow->T[currentTes].cellHandles.begin(); cell_it!=flow->T[currentTes].cellHandles.end(); cell_it++){ -// if (noCache) for (Finite_cells_iterator cell = Tri.finite_cells_begin(); cell != cell_end; cell++) { - Cell_handle& cell = *cell_it; - //reset cache - for (int k=0;k<4;k++) cell->info().unitForceVectors[k]=nullVect; - for (int j=0; j<4; j++) if (!Tri.is_infinite(cell->neighbor(j))) { - neighbour_cell = cell->neighbor(j); - const CGT::Vecteur& Surfk = cell->info().facetSurfaces[j]; - //FIXME : later compute that fluidSurf only once in hydraulicRadius, for now keep full surface not modified in cell->info for comparison with other forces schemes - //The ratio void surface / facet surface - Real area = sqrt(Surfk.squared_length()); if (area<=0) cerr <<"AREA <= 0!!"<<endl; - CGT::Vecteur facetNormal = Surfk/area; - const std::vector<CGT::Vecteur>& crossSections = cell->info().facetSphereCrossSections; - CGT::Vecteur fluidSurfk = cell->info().facetSurfaces[j]*cell->info().facetFluidSurfacesRatio[j]; - /// handle fictious vertex since we can get the projected surface easily here - if (cell->vertex(j)->info().isFictious) { - Real projSurf=abs(Surfk[flow->boundary(cell->vertex(j)->info().id()).coordinate]); - tempVect=-projSurf*flow->boundary(cell->vertex(j)->info().id()).normal; - cell->vertex(j)->info().forces = cell->vertex(j)->info().forces+tempVect*cell->info().p(); - //define the cached value for later use with cache*p - cell->info().unitForceVectors[j]=cell->info().unitForceVectors[j]+ tempVect; - } - /// Apply weighted forces f_k=sqRad_k/sumSqRad*f - CGT::Vecteur Facet_Unit_Force = -fluidSurfk*cell->info().solidLine[j][3]; - CGT::Vecteur Facet_Force = cell->info().p()*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 (Finite_cells_iterator cell = Tri.finite_cells_begin(); cell != cell_end; cell++) { - for (int yy=0;yy<4;yy++) cell->vertex(yy)->info().forces = cell->vertex(yy)->info().forces + cell->info().unitForceVectors[yy]*cell->info().p();} - - #else - #pragma omp parallel for num_threads(ompThreads) - for (int vn=0; vn<= flow->T[currentTes].max_id; vn++) { - Vertex_handle& v = flow->T[currentTes].vertexHandles[vn]; -// for (Finite_vertices_iterator v = Tri.finite_vertices_begin(); v != Tri.finite_vertices_end(); ++v){ - const int& id = v->info().id(); - CGT::Vecteur tf (0,0,0); - int k=0; - for (vector<const Real*>::iterator c = solver->perVertexPressure[id].begin(); c != solver->perVertexPressure[id].end(); c++) - tf = tf + (*(solver->perVertexUnitForce[id][k++]))*(**c); - v->info().forces = tf; - } - #endif - } - if (flow->DEBUG_OUT) { - CGT::Vecteur TotalForce = nullVect; - for (Finite_vertices_iterator v = Tri.finite_vertices_begin(); v != Tri.finite_vertices_end(); ++v) { - if (!v->info().isFictious) TotalForce = TotalForce + v->info().forces; - else if (flow->boundary(v->info().id()).flowCondition==1) TotalForce = TotalForce + v->info().forces; } - cout << "TotalForce = "<< TotalForce << endl;} -} - -template<class Solver> -void UnsaturatedEngine::computeSolidLine(Solver& flow) -{ - RTriangulation& Tri = flow->T[solver->currentTes].Triangulation(); - Finite_cells_iterator cell_end = Tri.finite_cells_end(); - for (Finite_cells_iterator cell = Tri.finite_cells_begin(); cell != cell_end; cell++) { - for(int j=0; j<4; j++) { - solver->Line_Solid_Pore(cell, j); - } - } -} - -YADE_PLUGIN ( ( UnsaturatedEngine ) ); - -#endif //FLOW_ENGINE - -#endif /* YADE_CGAL */ === removed file 'pkg/dem/UnsaturatedEngine.hpp' --- pkg/dem/UnsaturatedEngine.hpp 2014-04-04 17:12:06 +0000 +++ pkg/dem/UnsaturatedEngine.hpp 1970-01-01 00:00:00 +0000 @@ -1,241 +0,0 @@ -/************************************************************************* -* 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. * -*************************************************************************/ -#pragma once -#include<yade/core/PartialEngine.hpp> -#include<yade/pkg/dem/TriaxialCompressionEngine.hpp> -#include<yade/pkg/dem/TesselationWrapper.hpp> -#include<yade/lib/triangulation/FlowBoundingSphere.hpp> - -class Flow; -class TesselationWrapper; -#define _FlowSolver CGT::FlowBoundingSphere<UnsatTesselation> -#define TPL template<class Solver> - -class UnsaturatedEngine : public PartialEngine -{ - public : - typedef _FlowSolver FlowSolver; - typedef FlowTesselation Tesselation; - typedef FlowSolver::RTriangulation RTriangulation; - typedef FlowSolver::Finite_vertices_iterator Finite_vertices_iterator; - typedef FlowSolver::Finite_cells_iterator Finite_cells_iterator; - typedef FlowSolver::Cell_handle Cell_handle; - typedef RTriangulation::Finite_edges_iterator Finite_edges_iterator; - typedef RTriangulation::Vertex_handle Vertex_handle; - typedef RTriangulation::Point CGALSphere; - typedef CGALSphere::Point Point; - - typedef std::vector<Cell_handle> Vector_Cell; - typedef typename Vector_Cell::iterator VCell_iterator; - - protected: - shared_ptr<FlowSolver> solver; - shared_ptr<FlowSolver> backgroundSolver; - volatile bool backgroundCompleted; - Cell cachedCell; - struct posData {Body::id_t id; Vector3r pos; Real radius; bool isSphere; bool exists; posData(){exists=0;}}; - vector<posData> positionBufferCurrent;//reflect last known positions before we start computations - vector<posData> positionBufferParallel;//keep the positions from a given step for multithread factorization -// //copy positions in a buffer for faster and/or parallel access - void setPositionsBuffer(bool current); - void testFunction(); - - public : - enum {wall_xmin, wall_xmax, wall_ymin, wall_ymax, wall_zmin, wall_zmax}; - Vector3r normal [6]; - bool currentTes; - int id_offset; - TPL void initSolver (Solver& flow); - TPL void Triangulate (Solver& flow); - TPL void AddBoundary (Solver& flow); - TPL void Build_Triangulation (double P_zero, Solver& flow); - TPL void Build_Triangulation (Solver& flow); - TPL void Initialize_volumes (Solver& flow); - TPL void BoundaryConditions(Solver& flow); - TPL void initializeCellIndex(Solver& flow); - TPL void initializePoreRadius(Solver& flow); - - TPL void invade (Solver& flow ); - TPL Real getMinEntryValue (Solver& flow); - TPL Real getSaturation(Solver& flow); - TPL Real getWindowsSaturation(Solver&flow, int windowsID); - TPL Real getWindowsSaturation1(Solver&flow, int i); - TPL Real getWindowsSaturation2(Solver&flow, int i); - - TPL void invadeSingleCell1(Cell_handle cell, double pressure, Solver& flow); - TPL void invade1 (Solver& flow ); - TPL void checkTrap(Solver& flow, double pressure);//check trapped phase, define trapCapP. - TPL Real getMinEntryValue1 (Solver& flow); - TPL void updatePressure(Solver& flow); - TPL void updatePressureReservoir(Solver& flow); - TPL void initReservoirBound(Solver& flow); - TPL void initWaterReservoirBound(Solver& flow); - TPL void updateWaterReservoir(Solver& flow); - TPL void waterReservoirRecursion(Cell_handle cell, Solver& flow); - TPL void initAirReservoirBound(Solver& flow); - TPL void updateAirReservoir(Solver& flow); - TPL void airReservoirRecursion(Cell_handle cell, Solver& flow); - TPL Real getSaturation1(Solver& flow); - - TPL void invadeSingleCell2(Cell_handle cell, double pressure, Solver& flow); - TPL void invade2 (Solver& flow ); - TPL void updatePressure2(Solver& flow); - TPL Real getMinEntryValue2 (Solver& flow); - TPL Real getSaturation2(Solver& flow); - - TPL unsigned int imposePressure(Vector3r pos, Real p,Solver& flow); - TPL void setImposedPressure(unsigned int cond, Real p,Solver& flow); - TPL void clearImposedPressure(Solver& flow); - TPL void saveListNodes(Solver& flow); - TPL void saveListConnection(Solver& flow); - TPL void saveLatticeNodeX(Solver& flow,double x); - TPL void saveLatticeNodeY(Solver& flow,double y); - TPL void saveLatticeNodeZ(Solver& flow,double z); - TPL void saveReservoirInfo(Solver& flow,int boundN); - TPL void saveBoundingCellsInfo(Solver& flow); - TPL void savePoreBodyInfo(Solver& flow); - TPL void savePoreThroatInfo(Solver& flow); - TPL void debugTemp(Solver& flow); - TPL void initializeCellWindowsID(Solver&flow); - TPL void computeSolidLine(Solver& flow); - TPL void computeFacetPoreForcesWithCache(Solver& flow, bool onlyCache=false); - - 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]);} - - template<class Cellhandle > - double getRadiusMin(Cellhandle cell, int j); - template<class Cellhandle> - Real Volume_cell_single_fictious (Cellhandle cell); - template<class Cellhandle> - Real Volume_cell_double_fictious (Cellhandle cell); - template<class Cellhandle> - Real Volume_cell_triple_fictious (Cellhandle cell); - template<class Cellhandle> - Real Volume_cell (Cellhandle cell); - template<class Solver> - void updateVolumeCapillaryCell ( Solver& flow); - template<class Cellhandle> - double computeEffPoreRadius(Cellhandle cell, int j); - template<class Cellhandle> - double computeEffPoreRadiusNormal(Cellhandle cell, int j); - template<class Cellhandle> - double bisection(Cellhandle cell, int j, double a, double b); - template<class Cellhandle> - double computeDeltaForce(Cellhandle cell,int j, double rcap); - template<class Cellhandle> - Real computePoreArea(Cellhandle cell, int j); - template<class Cellhandle> - Real computePorePerimeter(Cellhandle cell, int j); - void saveVtk(const char* folder) {bool initT=solver->noCache; solver->noCache=false; solver->saveVtk(folder); solver->noCache=initT;} - python::list getConstrictions() { - vector<Real> csd=solver->getConstrictions(); python::list pycsd; - for (unsigned int k=0;k<csd.size();k++) pycsd.append(csd[k]); return pycsd;} - double getPorePressure(Vector3r pos){return solver->getPorePressure(pos[0], pos[1], pos[2]);} - TPL int getCell(double posX, double posY, double posZ, Solver& flow){return flow->getCell(posX, posY, posZ);} - - void emulateAction(){ - scene = Omega::instance().getScene().get(); - action();} - unsigned int _imposePressure(Vector3r pos, Real p) {return imposePressure(pos,p,solver);} - void _setImposedPressure(unsigned int cond, Real p) {setImposedPressure(cond,p,solver);} - void _clearImposedPressure() {clearImposedPressure(solver);} - int _getCell(Vector3r pos) {return getCell(pos[0],pos[1],pos[2],solver);} - void _buildTriangulation() {setPositionsBuffer(true); Build_Triangulation(solver);} - void _invade() {invade(solver);} - Real _getMinEntryValue() {return getMinEntryValue(solver);} - Real _getSaturation() {return getSaturation(solver);} - Real _getWindowsSaturation(int windowsID) {return getWindowsSaturation(solver,windowsID);} - void _saveListNodes() {saveListNodes(solver);} - void _saveListConnection() {saveListConnection(solver);} - void _saveLatticeNodeX(double x) {saveLatticeNodeX(solver,x);} - void _saveLatticeNodeY(double y) {saveLatticeNodeY(solver,y);} - void _saveLatticeNodeZ(double z) {saveLatticeNodeZ(solver,z);} - void _saveReservoirInfo(int boundN) {saveReservoirInfo(solver,boundN);} - void _saveBoundingCellsInfo(){saveBoundingCellsInfo(solver);} - void _savePoreBodyInfo(){savePoreBodyInfo(solver);} - void _savePoreThroatInfo(){savePoreThroatInfo(solver);} - void _debugTemp(){debugTemp(solver);} - void _initializeCellWindowsID(){initializeCellWindowsID(solver);} - void _computeFacetPoreForcesWithCache(){computeFacetPoreForcesWithCache(solver);} - Vector3r _fluidForce(unsigned int id_sph) {return fluidForce(id_sph,solver);} - - virtual ~UnsaturatedEngine(); - - virtual void action(); - - YADE_CLASS_BASE_DOC_ATTRS_DEPREC_INIT_CTOR_PY(UnsaturatedEngine,PartialEngine,"Preliminary version engine of a drainage model for unsaturated soils. Note:Air reservoir is on the top; water reservoir is on the bottom.", - ((bool,isActivated,true,,"Activate UnsaturatedEngine.")) - ((bool,first,true,,"Controls the initialization/update phases")) - ((bool, Debug, false,,"Activate debug messages")) - ((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, wall_thickness,0.001,,"Walls thickness")) - ((double,P_zero,0,,"The value used for initializing pore pressure. It is useless for incompressible fluid, but important for compressible model.")) - ((bool, updateTriangulation, 0,,"If true the medium is retriangulated.")) - ((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)")) - ((double, porosity, 0,,"Porosity computed at each retriangulation")) - ((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`.")) - ((int, ymax,3,(Attr::readonly),"See :yref:`FlowEngine::xmin`.")) - ((int, zmin,4,(Attr::readonly),"See :yref:`FlowEngine::xmin`.")) - ((int, zmax,5,(Attr::readonly),"See :yref:`FlowEngine::xmin`.")) - - ((vector<bool>, bndCondIsPressure, vector<bool>(6,false),,"defines the type of boundary condition for each side. True if pressure is imposed, False for no-flux. Indexes can be retrieved with :yref:`FlowEngine::xmin` and friends.")) - ((vector<double>, bndCondValue, vector<double>(6,0),,"Imposed value of a boundary condition. Only applies if the boundary condition is imposed pressure, else the imposed flux is always zero presently (may be generalized to non-zero imposed fluxes in the future).")) - //FIXME: to be implemented: - ((vector<Vector3r>, boundaryVelocity, vector<Vector3r>(6,Vector3r::Zero()),, "velocity on top boundary, only change it using :yref:`FlowEngine::setBoundaryVel`")) - ((vector<int>, wallIds,vector<int>(6),,"body ids of the boundaries (default values are ok only if aabbWalls are appended before spheres, i.e. numbered 0,...,5)")) - ((vector<bool>, boundaryUseMaxMin, vector<bool>(6,true),,"If true (default value) bounding sphere is added as function of max/min sphere coord, if false as function of yade wall position")) - - ((bool, pressureForce, true,,"Compute the pressure field and associated fluid forces. WARNING: turning off means fluid flow is not computed at all.")) - ((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).")) - , - /*deprec*/ - ,, - for (int i=0; i<6; ++i){normal[i]=Vector3r::Zero(); wallIds[i]=i;} - normal[wall_ymin].y()=normal[wall_xmin].x()=normal[wall_zmin].z()=1; - normal[wall_ymax].y()=normal[wall_xmax].x()=normal[wall_zmax].z()=-1; - solver = shared_ptr<FlowSolver> (new FlowSolver); - first=true; - , - .def("imposePressure",&UnsaturatedEngine::_imposePressure,(python::arg("pos"),python::arg("p")),"Impose pressure in cell of location 'pos'. The index of the condition is returned (for multiple imposed pressures at different points).") - .def("setImposedPressure",&UnsaturatedEngine::_setImposedPressure,(python::arg("cond"),python::arg("p")),"Set pressure value at the point indexed 'cond'.") - .def("clearImposedPressure",&UnsaturatedEngine::_clearImposedPressure,"Clear the list of points with pressure imposed.") - .def("getConstrictions",&UnsaturatedEngine::getConstrictions,"Get the list of constrictions (inscribed circle) for all finite facets.") - .def("saveVtk",&UnsaturatedEngine::saveVtk,(python::arg("folder")="./VTK"),"Save pressure field in vtk format. Specify a folder name for output.") - .def("getPorePressure",&UnsaturatedEngine::getPorePressure,(python::arg("pos")),"Measure pore pressure in position pos[0],pos[1],pos[2]") - .def("emulateAction",&UnsaturatedEngine::emulateAction,"get scene and run action (may be used to manipulate engine outside the main loop).") - .def("getCell",&UnsaturatedEngine::_getCell,(python::arg("pos")),"get id of the cell containing (X,Y,Z).") - .def("testFunction",&UnsaturatedEngine::testFunction,"The playground for Chao's experiments.") - .def("buildTriangulation",&UnsaturatedEngine::_buildTriangulation,"Triangulate spheres of the current scene.") - .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("saveListNodes",&UnsaturatedEngine::_saveListNodes,"Save the list of nodes.") - .def("saveListConnection",&UnsaturatedEngine::_saveListConnection,"Save the connections between cells.") - .def("saveLatticeNodeX",&UnsaturatedEngine::_saveLatticeNodeX,(python::arg("x")),"Save the slice of lattice nodes for x_normal(x). 0: out of sphere; 1: inside of sphere.") - .def("saveLatticeNodeY",&UnsaturatedEngine::_saveLatticeNodeY,(python::arg("y")),"Save the slice of lattice nodes for y_normal(y). 0: out of sphere; 1: inside of sphere.") - .def("saveLatticeNodeZ",&UnsaturatedEngine::_saveLatticeNodeZ,(python::arg("z")),"Save the slice of lattice nodes for z_normal(z). 0: out of sphere; 1: inside of sphere.") - .def("saveReservoirInfo",&UnsaturatedEngine::_saveReservoirInfo,(python::arg("boundN")),"Test reservoir cells statement.(temporary)") - .def("saveBoundingCellsInfo",&UnsaturatedEngine::_saveBoundingCellsInfo,"Test boundary cells (without reservoir) statement.(temporary)") - .def("savePoreBodyInfo",&UnsaturatedEngine::_savePoreBodyInfo,"Save pore bodies positions/Voronoi centers and size/volume.") - .def("savePoreThroatInfo",&UnsaturatedEngine::_savePoreThroatInfo,"Save pore throat area, inscribed radius and perimeter.") - .def("debugTemp",&UnsaturatedEngine::_debugTemp,"debug temp file.") - .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. ") - .def("fluidForce",&UnsaturatedEngine::_fluidForce,(python::arg("Id_sph")),"Return the fluid force on sphere Id_sph.") - ) - DECLARE_LOGGER; -}; - -REGISTER_SERIALIZABLE(UnsaturatedEngine);
_______________________________________________ Mailing list: https://launchpad.net/~yade-dev Post to : [email protected] Unsubscribe : https://launchpad.net/~yade-dev More help : https://help.launchpad.net/ListHelp
