yade-dev team mailing list archive

[noreply@xxxxxxxxxxxxx]

------------------------------------------------------------
revno: 3929
committer: Bruno Chareyre <bruno.chareyre@xxxxxxxxxxxxxxx>
timestamp: Thu 2014-04-24 23:37:44 +0200
message:
  Add user defined boundary condition to FlowEngine
modified:
  lib/triangulation/FlowBoundingSphere.hpp
  lib/triangulation/FlowBoundingSphere.ipp
  pkg/pfv/FlowEngine.hpp
  pkg/pfv/FlowEngine.ipp


--
lp:yade
https://code.launchpad.net/~yade-pkg/yade/git-trunk

Your team Yade developers is subscribed to branch lp:yade.
To unsubscribe from this branch go to https://code.launchpad.net/~yade-pkg/yade/git-trunk/+edit-subscription

=== modified file 'lib/triangulation/FlowBoundingSphere.hpp'
--- lib/triangulation/FlowBoundingSphere.hpp	2014-04-07 09:33:19 +0000
+++ lib/triangulation/FlowBoundingSphere.hpp	2014-04-24 21:37:44 +0000
@@ -52,6 +52,8 @@
 		vector<CellHandle> IPCells;
 		vector<pair<Point,Real> > imposedF;
 		vector<CellHandle> IFCells;
+		//Pointers to vectors used for user defined boundary pressure
+		vector<Real> *pxpos, *ppval;
 		
 		void initNewTri () {noCache=true; /*isLinearSystemSet=false; areCellsOrdered=false;*/}//set flags after retriangulation
 		bool permeabilityMap;
@@ -146,6 +148,7 @@
 		virtual void averageRelativeCellVelocity();
 		void averageFluidVelocity();
 		void applySinusoidalPressure(RTriangulation& Tri, double amplitude, double averagePressure, double loadIntervals);
+		void applyUserDefinedPressure(RTriangulation& Tri, vector<Real>& xpos, vector<Real>& pval);
 		bool isOnSolid  (double X, double Y, double Z);
 		double getPorePressure (double X, double Y, double Z);
 		void measurePressureProfile(double WallUpy, double WallDowny);

=== modified file 'lib/triangulation/FlowBoundingSphere.ipp'
--- lib/triangulation/FlowBoundingSphere.ipp	2014-04-15 16:39:53 +0000
+++ lib/triangulation/FlowBoundingSphere.ipp	2014-04-24 21:37:44 +0000
@@ -86,6 +86,7 @@
 	maxKdivKmean=100.;
 	ompThreads=1;
 	errorCode=0;
+	pxpos=ppval=NULL;
 }
 
 template <class Tesselation> 
@@ -397,6 +398,29 @@
 	  }
 	}
 }
+
+template <class Tesselation> 
+void FlowBoundingSphere<Tesselation>::applyUserDefinedPressure(RTriangulation& Tri, vector<Real>& xpos, vector<Real>& pval)
+{
+	if (!(xpos.size() && xpos.size()==pval.size())) {cerr << "Wrong definition of boundary pressure, check input" <<endl; return;}
+	pxpos=&xpos; ppval=&pval;
+	Real dx = xpos[1] - xpos[0]; Real xinit=xpos[0]; Real xlast=xpos.back();
+	VectorCell tmpCells; tmpCells.resize(10000);
+	VCellIterator cellsEnd = Tri.incident_cells(T[currentTes].vertexHandles[yMaxId],tmpCells.begin());
+	for (VCellIterator it = tmpCells.begin(); it != cellsEnd; it++)
+	{
+		if(Tri.is_infinite(*it)) continue;
+		Point& p1 = (*it)->info();
+		CellHandle& cell = *it;
+		if (p1.x()<xinit || p1.x()>xlast) cerr<<"udef pressure: cell out of range"<<endl;
+		else {
+			Real frac, intg;
+			frac=modf((p1.x()-xinit)/dx,&intg);
+			cell->info().p() = pval[intg]*(1-frac) + pval[intg+1]*frac;
+		}
+	}
+}
+
 template <class Tesselation> 
 void FlowBoundingSphere<Tesselation>::interpolate(Tesselation& Tes, Tesselation& NewTes)
 {
@@ -725,6 +749,8 @@
 			}
                 }
         }
+        if (ppval && pxpos) applyUserDefinedPressure(Tri,*pxpos,*ppval);
+        
         IPCells.clear();
         for (unsigned int n=0; n<imposedP.size();n++) {
 		CellHandle cell=Tri.locate(imposedP[n].first);
@@ -763,6 +789,7 @@
 			}
                 }
         }
+        if (ppval && pxpos) applyUserDefinedPressure(T[currentTes].Triangulation(),*pxpos,*ppval);
         for (unsigned int n=0; n<imposedP.size();n++) {
 		IPCells[n]->info().p()=imposedP[n].second;
 		IPCells[n]->info().Pcondition=true;}

=== modified file 'pkg/pfv/FlowEngine.hpp'
--- pkg/pfv/FlowEngine.hpp	2014-04-17 15:10:23 +0000
+++ pkg/pfv/FlowEngine.hpp	2014-04-24 21:37:44 +0000
@@ -292,6 +292,8 @@
 		((int, numSolveThreads, 1,,"number of openblas threads in the solve phase."))
 		((int, numFactorizeThreads, 1,,"number of openblas threads in the factorization phase"))
 		#endif
+		((vector<Real>, boundaryPressure,vector<Real>(),,"values defining pressure along x-axis for the top surface. See also :yref:`TEMPLATE_FLOW_NAME::boundaryXPos`"))
+		((vector<Real>, boundaryXPos,vector<Real>(),,"values of the x-coordinate for which pressure is defined. See also :yref:`TEMPLATE_FLOW_NAME::boundaryPressure`"))
 		,
 		/*deprec*/
 		((meanK_opt,clampKValues,"the name changed"))

=== modified file 'pkg/pfv/FlowEngine.ipp'
--- pkg/pfv/FlowEngine.ipp	2014-04-17 15:10:23 +0000
+++ pkg/pfv/FlowEngine.ipp	2014-04-24 21:37:44 +0000
@@ -278,6 +278,7 @@
 	
         if ( !first && !multithread && (useSolver==0 || fluidBulkModulus>0 || doInterpolate)) flow.interpolate ( flow.T[!flow.currentTes], flow.T[flow.currentTes] );
         if ( waveAction ) flow.applySinusoidalPressure ( flow.T[flow.currentTes].Triangulation(), sineMagnitude, sineAverage, 30 );
+	else if (boundaryPressure.size()!=0) flow.applyUserDefinedPressure ( flow.T[flow.currentTes].Triangulation(), boundaryXPos , boundaryPressure);
         if (normalLubrication || shearLubrication || viscousShear) flow.computeEdgesSurfaces();
 }
 template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >