yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #12211
[Branch ~yade-pkg/yade/git-trunk] Rev 3707: add a new version of capillary law. (pushed by Caroline)
------------------------------------------------------------
revno: 3707
committer: Chao Yuan <chaoyuan2012@xxxxxxxxx>
timestamp: Thu 2015-07-23 19:26:21 +0200
message:
add a new version of capillary law. (pushed by Caroline)
added:
pkg/dem/CapillaryPhys1.cpp
pkg/dem/CapillaryPhys1.hpp
pkg/dem/DelaunayInterpolation.hpp
pkg/dem/Law2_ScGeom_CapillaryPhys_Capillarity1.cpp
pkg/dem/Law2_ScGeom_CapillaryPhys_Capillarity1.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
=== added file 'pkg/dem/CapillaryPhys1.cpp'
--- pkg/dem/CapillaryPhys1.cpp 1970-01-01 00:00:00 +0000
+++ pkg/dem/CapillaryPhys1.cpp 2015-07-23 17:26:21 +0000
@@ -0,0 +1,58 @@
+//keep this #ifdef as long as you don't really want to realize a final version publicly, it will save compilation time for everyone else
+//when you want it compiled, you can just uncomment the following line
+// #define CAPILLARYPHYS1
+#ifdef CAPILLARYPHYS1
+
+#include <pkg/dem/CapillaryPhys1.hpp>
+#include<pkg/dem/ScGeom.hpp>
+
+#include<core/Omega.hpp>
+#include<core/Scene.hpp>
+
+CapillaryPhys1::~CapillaryPhys1()
+{
+}
+
+YADE_PLUGIN((CapillaryPhys1));
+
+
+
+YADE_PLUGIN((Ip2_FrictMat_FrictMat_CapillaryPhys1));
+
+void Ip2_FrictMat_FrictMat_CapillaryPhys1::go( const shared_ptr<Material>& b1 //FrictMat
+ , const shared_ptr<Material>& b2 // FrictMat
+ , const shared_ptr<Interaction>& interaction)
+{
+ ScGeom* geom = YADE_CAST<ScGeom*>(interaction->geom.get());
+ if(geom)
+ {
+
+ if(!interaction->phys)
+ {
+ const shared_ptr<FrictMat>& sdec1 = YADE_PTR_CAST<FrictMat>(b1);
+ const shared_ptr<FrictMat>& sdec2 = YADE_PTR_CAST<FrictMat>(b2);
+
+ if (!interaction->phys) interaction->phys = shared_ptr<CapillaryPhys1>(new CapillaryPhys1());
+ const shared_ptr<CapillaryPhys1>& contactPhysics = YADE_PTR_CAST<CapillaryPhys1>(interaction->phys);
+
+ Real Ea = sdec1->young;
+ Real Eb = sdec2->young;
+ Real Va = sdec1->poisson;
+ Real Vb = sdec2->poisson;
+ Real Da = geom->radius1; // FIXME - multiply by factor of sphere interaction distance (so sphere interacts at bigger range that its geometrical size)
+ Real Db = geom->radius2; // FIXME - as above
+ Real fa = sdec1->frictionAngle;
+ Real fb = sdec2->frictionAngle;
+ Real Kn = 2*Ea*Da*Eb*Db/(Ea*Da+Eb*Db);//harmonic average of two stiffnesses
+ Real Ks = 2*Ea*Da*Va*Eb*Db*Vb/(Ea*Da*Va+Eb*Db*Va);//harmonic average of two stiffnesses with ks=V*kn for each sphere
+
+ contactPhysics->tangensOfFrictionAngle = std::tan(std::min(fa,fb));
+ contactPhysics->kn = Kn;
+ contactPhysics->ks = Ks;
+ contactPhysics->computeBridge =computeDefault;
+ }
+ }
+};
+
+#endif //CAPILLARYPHYS1
+
=== added file 'pkg/dem/CapillaryPhys1.hpp'
--- pkg/dem/CapillaryPhys1.hpp 1970-01-01 00:00:00 +0000
+++ pkg/dem/CapillaryPhys1.hpp 2015-07-23 17:26:21 +0000
@@ -0,0 +1,99 @@
+/*************************************************************************
+* Copyright (C) 2006 by luc Scholtes *
+* luc.scholtes@xxxxxxxxxxx *
+* *
+* This program is free software; it is licensed under the terms of the *
+* GNU General Public License v2 or later. See file LICENSE for details. *
+*************************************************************************/
+#pragma once
+#include<pkg/dem/FrictPhys.hpp>
+#include<pkg/dem/DelaunayInterpolation.hpp>
+#include<pkg/common/Dispatching.hpp>
+#include<pkg/common/ElastMat.hpp>
+
+class MeniscusPhysicalData {
+public:
+ double R;
+ double volume;
+ double distance;
+ double surface;
+ double energy;
+ double force;
+ double succion;
+ double delta1;
+ double delta2;
+ double arcLength;
+ bool ending;
+ //default ctor
+ MeniscusPhysicalData() : R(0), volume(0), distance(0), surface(0), energy(0), force(0), succion(0), delta1(0), delta2(0), arcLength(0), ending(1) {}
+ //ctor with input values
+ MeniscusPhysicalData(const double& r, const double& v, const double& d, const double& s, const double& e, const double& f, const double& p, const double& a1, const double& a2, const double& arc, bool end) : R(r), volume(v), distance(d), surface(s), energy(e), force(f), succion(p), delta1(a1), delta2(a2), arcLength(arc), ending(end) {}
+
+ //a minimal list of operators for the interpolation
+ //these operators are requirements for the DataType template parameter of interpolate()
+ //FIXME: algebra operations must include energy, perimeter, and any other new variable for including them in the interpolation
+ MeniscusPhysicalData& operator+= (const MeniscusPhysicalData& m2) {
+ R+=m2.R; volume+=m2.volume; distance+=m2.distance; surface+=m2.surface; energy+=m2.energy; force+=m2.force; succion+=m2.succion; delta1+=m2.delta1; delta2+=m2.delta2; arcLength+=m2.arcLength; ending=(ending && m2.ending);
+ return *this;}
+
+ MeniscusPhysicalData operator* (const double& fact) const {
+ return MeniscusPhysicalData(fact*R, fact*volume, fact*distance, fact*surface, fact*energy, fact*force, fact*succion, fact*delta1, fact*delta2, fact*arcLength, ending);}
+
+ const MeniscusPhysicalData& operator= (const MeniscusPhysicalData& m1) {
+ R=m1.R; volume=m1.volume; distance=m1.distance; surface=m1.surface; energy=m1.energy; force=m1.force; succion=m1.succion; delta1=m1.delta1; delta2=m1.delta2; arcLength=m1.arcLength; ending=m1.ending;
+ return *this;}
+};
+
+//The structure for the meniscus: physical properties + cached values for fast interpolation
+class Meniscus {
+ public:
+ typedef MeniscusPhysicalData Data;
+ Data data; //the variables of Laplace's problem
+ DT::Cell_handle cell; //pointer to the last location in the triangulation, for faster locate()
+ std::vector<K::Vector_3> normals;// 4 normals relative to the current cell
+
+ Meniscus() : data(), cell(DT::Cell_handle()), normals(std::vector<K::Vector_3>()) {}
+};
+
+class CapillaryPhys1 : public FrictPhys
+{
+ public :
+ int currentIndexes [4]; // used for faster interpolation (stores previous positions in tables)
+ Meniscus m;
+
+ virtual ~CapillaryPhys1();
+
+ YADE_CLASS_BASE_DOC_ATTRS_CTOR(CapillaryPhys1,FrictPhys,"Physics (of interaction) for Law2_ScGeom_CapillaryPhys_Capillarity.",
+ ((bool,meniscus,false,,"Presence of a meniscus if true"))
+ ((bool,isBroken,false,,"If true, capillary force is zero and liquid bridge is inactive."))
+ ((bool,computeBridge,true,,"If true, capillary bridge will be computed if not it will be ignored."))
+ ((Real,capillaryPressure,0.,,"Value of the capillary pressure Uc defines as Ugas-Uliquid"))
+ ((Real,vMeniscus,0.,,"Volume of the menicus"))
+ ((Real,Delta1,0.,,"Defines the surface area wetted by the meniscus on the smallest grains of radius R1 (R1<R2)"))
+ ((Real,Delta2,0.,,"Defines the surface area wetted by the meniscus on the biggest grains of radius R2 (R1<R2)"))
+ ((Vector3r,fCap,Vector3r::Zero(),,"Capillary Force produces by the presence of the meniscus"))
+ ((Real,SInterface,0.,,"Fluid-Gaz Interfacial area"))
+ ((Real,arcLength,0.,,"Arc Length of the Fluid-Gaz Interface"))
+ ((short int,fusionNumber,0.,,"Indicates the number of meniscii that overlap with this one"))
+ ,createIndex();currentIndexes[0]=currentIndexes[1]=currentIndexes[2]=currentIndexes[3]=0;
+ );
+ REGISTER_CLASS_INDEX(CapillaryPhys1,FrictPhys);
+};
+REGISTER_SERIALIZABLE(CapillaryPhys1);
+
+class Ip2_FrictMat_FrictMat_CapillaryPhys1 : public IPhysFunctor
+{
+ public :
+ virtual void go( const shared_ptr<Material>& b1,
+ const shared_ptr<Material>& b2,
+ const shared_ptr<Interaction>& interaction);
+
+ FUNCTOR2D(FrictMat,FrictMat);
+ YADE_CLASS_BASE_DOC_ATTRS_CTOR(Ip2_FrictMat_FrictMat_CapillaryPhys1,IPhysFunctor, "RelationShips to use with Law2_ScGeom_CapillaryPhys_Capillarity1\n\n In these RelationShips all the interaction attributes are computed. \n\n.. warning::\n\tas in the others :yref:`Ip2 functors<IPhysFunctor>`, most of the attributes are computed only once, when the interaction is new.",
+ ((bool,computeDefault,true,,"bool to assign the default value of computeBridge.")),;
+
+ );
+
+};
+REGISTER_SERIALIZABLE(Ip2_FrictMat_FrictMat_CapillaryPhys1);
+
=== added file 'pkg/dem/DelaunayInterpolation.hpp'
--- pkg/dem/DelaunayInterpolation.hpp 1970-01-01 00:00:00 +0000
+++ pkg/dem/DelaunayInterpolation.hpp 2015-07-23 17:26:21 +0000
@@ -0,0 +1,150 @@
+/*************************************************************************
+* Copyright (C) 2013 by Bruno Chareyre <bruno.chareyre@xxxxxxxxxxx> *
+* *
+* This program is free software; it is licensed under the terms of the *
+* GNU General Public License v2 or later. See file LICENSE for details. *
+*************************************************************************/
+
+#pragma once
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+#include <CGAL/Delaunay_triangulation_3.h>
+#include <CGAL/Triangulation_vertex_base_with_info_3.h>
+#include <CGAL/Cartesian.h>
+#include <iostream>
+#include <fstream>
+/////////////////////////////////////////////////////////////////////
+/*
+TYPES:
+Triangulation_vertex_base_with_id_3: we redefine a vertex base including an index for each vertex (available in CGAL in 2D but not in 3D),
+MeniscusPhysicalData: the physical variables describing a capillary bridge, with a few algebraic operators for computing weighted averages
+Meniscus: a structure combining MeniscusPhysicalData with some cached data allowing faster operations in multiple queries (pointer to the last cell found and its normals)
+
+FUNCTIONS:
+getIncidentVtxWeights: an interpolation algorithm which is 20x faster than the natural neighbor interpolation of CGAL.
+ returns a list of vertices incident to a query point and their respective weights
+interpolate: uses the results of getIncidentVtxWeights combined with a data array to return a weighted average,
+ may be used with arbitrary data types provided they have the required algebraic operators
+main: example usage
+*/
+/////////////////////////////////////////////////////////////////////
+
+namespace CGAL {
+
+
+
+//Vertex base including an index for each vertex, adapted from CGAL::Triangulation_vertex_base_with_id_2
+template < typename GT, typename Vb = Triangulation_vertex_base_with_info_3<unsigned,GT> >
+class Triangulation_vertex_base_with_id_3 : public Vb
+{
+ int _id;
+public:
+ typedef typename Vb::Cell_handle Cell_handle;
+ typedef typename Vb::Point Point;
+ template < typename TDS3 >
+ struct Rebind_TDS {
+ typedef typename Vb::template Rebind_TDS<TDS3>::Other Vb3;
+ typedef Triangulation_vertex_base_with_id_3<GT, Vb3> Other;
+ };
+
+ Triangulation_vertex_base_with_id_3(): Vb() {}
+ Triangulation_vertex_base_with_id_3(const Point & p): Vb(p) {}
+ Triangulation_vertex_base_with_id_3(const Point & p, Cell_handle c): Vb(p, c) {}
+ Triangulation_vertex_base_with_id_3(Cell_handle c): Vb(c) {}
+ unsigned int id() const { return this->info(); }
+ unsigned int& id() { return this->info(); }
+};
+
+// The function returning vertices and their weights for an arbitrary point in R3 space.
+// The returned triplet contains:
+// - the output iterator pointing to the filled vector of vertices
+// - the sum of weights for normalisation
+// - a bool telling if the query point was located inside the convex hull of the data points
+template <class Dt, class OutputIterator>
+Triple< OutputIterator, // iterator with value type std::pair<Dt::Vertex_handle, Dt::Geom_traits::FT>
+ typename Dt::Geom_traits::FT, // Should provide 0 and 1
+ bool >
+getIncidentVtxWeights(const Dt& dt,
+ const typename Dt::Geom_traits::Point_3& Q,
+ OutputIterator nn_out, typename Dt::Geom_traits::FT & norm_coeff,
+ std::vector<typename Dt::Geom_traits::Vector_3>& normals,
+ typename Dt::Cell_handle& start = CGAL_TYPENAME_DEFAULT_ARG Dt::Cell_handle())
+{
+ //helpful array for permutations
+ const int comb [6] = {1, 2, 3, 0, 1, 2};
+ typedef typename Dt::Geom_traits Gt;
+ typedef typename Gt::Point_3 Point;
+ typedef typename Dt::Cell_handle Cell_handle;
+ typedef typename Dt::Locate_type Locate_type;
+ typedef typename Gt::FT Coord_type;
+ CGAL_triangulation_precondition (dt.dimension()== 3);
+ Locate_type lt; int li, lj;
+ Cell_handle c = dt.locate( Q, lt, li, lj, start);
+ bool updateNormals = (c!=start || normals.size()<4);
+ if (updateNormals) normals.clear();
+ if ( lt == Dt::VERTEX )
+ {
+ *nn_out++= std::make_pair(c->vertex(li),Coord_type(1));
+ return make_triple(nn_out,norm_coeff=Coord_type(1),true);
+ }
+ else if (dt.is_infinite(c))
+ return make_triple(nn_out, Coord_type(1), false);//point outside the convex-hull
+ norm_coeff=0;
+ for ( int k=0;k<4;k++ )
+ {
+ if (updateNormals) {
+ normals.push_back(cross_product(c->vertex(comb[k])->point()-c->vertex(comb[k+1])->point(),c->vertex(comb[k])->point()-c->vertex(comb[k+2])->point()));
+ normals[k] = normals[k]/
+ ((c->vertex(k)->point()-c->vertex(comb[k])->point())*normals[k]);
+ }
+ Coord_type closeness = ((Q-c->vertex(comb[k])->point())*normals[k]);
+ Coord_type w = closeness;
+ *nn_out++= std::make_pair(c->vertex(k),w);
+ norm_coeff += w;
+ }
+ start = c;
+ return make_triple(nn_out,norm_coeff,true);
+}
+
+
+
+} //END NAMESPACE CGAL
+
+typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
+typedef CGAL::Delaunay_triangulation_3<K>::Geom_traits Traits;
+typedef CGAL::Triangulation_vertex_base_with_id_3<Traits> Vb;
+typedef CGAL::Triangulation_cell_base_3<Traits> Cb;
+typedef CGAL::Triangulation_data_structure_3<Vb, Cb> Tds;
+typedef CGAL::Delaunay_triangulation_3<Traits,Tds> DT;
+typedef std::vector< std::pair< DT::Vertex_handle, K::FT> > Vertex_weight_vector;
+
+template <class Dt, class DataOwner>
+typename DataOwner::Data interpolate1 (const Dt& dt, const typename Dt::Geom_traits::Point_3& Q, DataOwner& owner, const std::vector<typename DataOwner::Data>& rawData, bool reset)
+{
+ K::FT norm;
+ Vertex_weight_vector coords;
+ if (reset) owner.cell = dt.cells_begin();
+ CGAL::Triple<std::back_insert_iterator<Vertex_weight_vector>,K::FT, bool> result = CGAL::getIncidentVtxWeights(dt, Q,std::back_inserter(coords), norm, owner.normals , owner.cell);
+
+ typename DataOwner::Data data = typename DataOwner::Data();//initialize null solution
+ if (!result.third) return data;// out of the convex hull, we return the null solution
+ //else, we compute the weighted sum
+ for (unsigned int k=0; k<coords.size(); k++) data += (rawData[coords[k].first->id()]*coords[k].second);
+ if (!data.ending) return data*(1./result.second);
+ else return typename DataOwner::Data();
+}
+template <class Dt, class DataOwner>
+typename DataOwner::Data interpolate2 (const Dt& dt, const typename Dt::Geom_traits::Point_3& Q, DataOwner& owner, const std::vector<typename DataOwner::Data>& rawData, bool reset)
+{
+ K::FT norm;
+ Vertex_weight_vector coords;
+ if (reset) owner.cell = dt.cells_begin();
+ CGAL::Triple<std::back_insert_iterator<Vertex_weight_vector>,K::FT, bool> result = CGAL::getIncidentVtxWeights(dt, Q,std::back_inserter(coords), norm, owner.normals , owner.cell);
+
+ typename DataOwner::Data data = typename DataOwner::Data();//initialize null solution
+ if (!result.third) return data;// out of the convex hull, we return the null solution
+ //else, we compute the weighted sum
+ for (unsigned int k=0; k<coords.size(); k++) data += (rawData[coords[k].first->id()]*coords[k].second);
+ if (!data.ending) return data*(1./result.second);
+ else return typename DataOwner::Data();
+}
+
=== added file 'pkg/dem/Law2_ScGeom_CapillaryPhys_Capillarity1.cpp'
--- pkg/dem/Law2_ScGeom_CapillaryPhys_Capillarity1.cpp 1970-01-01 00:00:00 +0000
+++ pkg/dem/Law2_ScGeom_CapillaryPhys_Capillarity1.cpp 2015-07-23 17:26:21 +0000
@@ -0,0 +1,606 @@
+/*************************************************************************
+* Copyright (C) 2006 by luc Scholtes *
+* luc.scholtes@xxxxxxxxxxx *
+* *
+* This program is free software; it is licensed under the terms of the *
+* GNU General Public License v2 or later. See file LICENSE for details. *
+*************************************************************************/
+
+//Modifs : Parameters renamed as MeniscusParameters
+//id1/id2 as id1 is the smallest grain, FIXME : wetting angle?
+//FIXME : in triaxialStressController, change test about null force in updateStiffnessccc
+//keep this #ifdef as long as you don't really want to realize a final version publicly, it will save compilation time for everyone else
+//when you want it compiled, you can just uncomment the following line
+// #define LAW2_SCGEOM_CAPILLARYPHYS_Capillarity1
+#ifdef LAW2_SCGEOM_CAPILLARYPHYS_Capillarity1
+
+#include <pkg/dem/Law2_ScGeom_CapillaryPhys_Capillarity1.hpp>
+#include <pkg/common/ElastMat.hpp>
+
+#include <pkg/dem/ScGeom.hpp>
+#include <pkg/dem/HertzMindlin.hpp>
+#include <core/Omega.hpp>
+#include <core/Scene.hpp>
+#include <lib/base/Math.hpp>
+#include <iostream>
+#include <fstream>
+
+
+DT Law2_ScGeom_CapillaryPhys_Capillarity1::dtVbased;
+DT Law2_ScGeom_CapillaryPhys_Capillarity1::dtPbased;
+
+ Real Law2_ScGeom_CapillaryPhys_Capillarity1::intEnergy()
+{
+ Real energy=0;
+ FOREACH(const shared_ptr<Interaction>& I, *scene->interactions){
+ if(!I->isReal()) continue;
+ ScGeom* currentGeometry = static_cast<ScGeom*>(I->geom.get());
+ CapillaryPhys1* phys = dynamic_cast<CapillaryPhys1*>(I->phys.get());
+ if(phys) {
+ if (phys->SInterface!=0){
+ energy += liquidTension*(phys->SInterface-4*3.141592653589793238462643383279502884*(pow(currentGeometry->radius1,2))-4*3.141592653589793238462643383279502884*(pow(currentGeometry->radius2,2)));}
+ }
+ }
+ return energy;
+}
+
+
+ Real Law2_ScGeom_CapillaryPhys_Capillarity1::wnInterface()
+{
+ Real wn=0;
+ FOREACH(const shared_ptr<Interaction>& I, *scene->interactions){
+ if(!I->isReal()) continue;
+ ScGeom* currentGeometry = static_cast<ScGeom*>(I->geom.get());
+ CapillaryPhys1* phys = dynamic_cast<CapillaryPhys1*>(I->phys.get());
+ if(phys) {
+ if (phys->SInterface!=0){
+
+ wn += (phys->SInterface-2*3.141592653589793238462643383279502884*(pow(currentGeometry->radius1,2)*(1+cos(phys->Delta1))+pow(currentGeometry->radius2,2)*(1+cos(phys->Delta2))));
+ }
+ }
+ }
+ return wn;
+}
+
+ Real Law2_ScGeom_CapillaryPhys_Capillarity1::swInterface()
+{
+ Real sw=0;
+ FOREACH(const shared_ptr<Interaction>& I, *scene->interactions){
+ if(!I->isReal()) continue;
+ ScGeom* currentGeometry = static_cast<ScGeom*>(I->geom.get());
+ CapillaryPhys1* phys = dynamic_cast<CapillaryPhys1*>(I->phys.get());
+ if(phys) {
+ sw += (2*3.141592653589793238462643383279502884*(pow(currentGeometry->radius1,2)*(1-cos(phys->Delta1))+pow(currentGeometry->radius2,2)*(1-cos(phys->Delta2))));
+
+ }
+ }
+ return sw;
+}
+
+
+ Real Law2_ScGeom_CapillaryPhys_Capillarity1::waterVolume()
+{
+ Real volume=0;
+ FOREACH(const shared_ptr<Interaction>& I, *scene->interactions){
+ if(!I->isReal()) continue;
+ CapillaryPhys1* phys = dynamic_cast<CapillaryPhys1*>(I->phys.get());
+ if(phys) {
+ volume += phys->vMeniscus;}
+ }
+ return volume;
+}
+
+void Law2_ScGeom_CapillaryPhys_Capillarity1::triangulateData() {
+ /// We get data from a file and input them in triangulations
+ if (solutions.size()>0) {LOG_WARN("Law2_ScGeom_CapillaryPhys_Capillarity1 asking triangulation for the second time. Ignored."); return;}
+ ifstream file (inputFilename.c_str());
+ if (!file.is_open()) { LOG_ERROR("No data file found for capillary law. Check path and inputFilename."); return;}
+
+ // convention R,v,d,s,e,f,p,a1,a2,dummy (just for the example, define your own,
+ // dummy is because has too much values per line - with one useless extra colum,)
+ MeniscusPhysicalData dat;
+ double ending;
+ while ( file.good() ) {
+ file >>dat.succion>>dat.force>>dat.distance>>dat.volume>>dat.surface>>dat.arcLength>>dat.delta1>>dat.delta2>>dat.R>>ending;
+ dat.ending=(bool) ending;
+ solutions.push_back(dat);
+ }
+ file.close();
+ // Make lists of points with index, so we can use range insertion, more efficient
+ // see http://doc.cgal.org/latest/Triangulation_3/index.html#Triangulation_3SettingInformationWhileInserting
+ std::vector< std::pair<K::Point_3,unsigned> > pointsP, pointsV;
+ for (unsigned int k=0; k<solutions.size(); k++) {
+ pointsP.push_back(std::make_pair(K::Point_3(solutions[k].R, solutions[k].succion, solutions[k].distance),k));
+ pointsV.push_back(std::make_pair(K::Point_3(solutions[k].R, solutions[k].volume, solutions[k].distance),k));
+ }
+ // and now range insertion
+ dtPbased.insert(pointsP.begin(), pointsP.end());
+ dtVbased.insert(pointsV.begin(), pointsV.end());
+}
+
+YADE_PLUGIN((Law2_ScGeom_CapillaryPhys_Capillarity1));
+
+int firstIteration=1;
+int x=0;
+
+void Law2_ScGeom_CapillaryPhys_Capillarity1::action()
+{
+ bool switched = (switchTriangulation == (imposePressure or totalVolumeConstant));
+ switchTriangulation = (imposePressure or totalVolumeConstant);
+ InteractionContainer::iterator ii = scene->interactions->begin();
+ InteractionContainer::iterator iiEnd = scene->interactions->end();
+ if (imposePressure) {
+ solver(capillaryPressure,switched);
+ }
+ else{
+ if (((totalVolumeConstant || (!totalVolumeConstant && firstIteration==1)) && totalVolumeofWater!=-1) || (totalVolumeConstant && totalVolumeofWater==-1))
+ {
+ if (!totalVolumeConstant) x=1;
+ totalVolumeConstant=1;
+ Real p0=capillaryPressure;
+ Real slope;
+ Real eps=0.0000000001;
+ solver(p0,switched);
+ Real V0=waterVolume();
+ if (totalVolumeConstant && totalVolumeofWater==-1 && firstIteration==1){
+ totalVolumeofWater=V0;
+ firstIteration+=1;
+ }
+ Real p1=capillaryPressure+0.1;
+ solver(p1,switched);
+ Real V1=waterVolume();
+ while (abs((totalVolumeofWater-V1)/totalVolumeofWater)>eps){
+ slope= (p1-p0)/(V1-V0);
+ p0=p1;
+ V0=V1;
+ p1=p1-slope*(V1-totalVolumeofWater);
+ if (p1<0) {
+ cout<< "The requested volume of water is quite big, the simulation will continue at constant suction.:"<< p0 <<endl;
+ capillaryPressure=p0;
+ imposePressure=1;
+ break;
+ }
+ solver(p1,switched);
+ V1=waterVolume();
+ capillaryPressure=p1;
+ }
+ if (x==1){
+ totalVolumeConstant=0;
+ firstIteration+=1;
+ }
+ }
+ else{
+ if ((!totalVolumeConstant && firstIteration==1) && totalVolumeofWater==-1){
+ totalVolumeConstant=1;
+ solver(capillaryPressure,switched);
+ firstIteration+=1;
+ totalVolumeConstant=0;
+ }
+ else
+ {
+ solver(capillaryPressure,switched);
+ }
+
+ }
+ }
+ for (ii= scene->interactions->begin(); ii!=iiEnd ; ++ii) {
+ CapillaryPhys1* phys = dynamic_cast<CapillaryPhys1*>((*ii)->phys.get());
+ if ((*ii)->isReal() && phys-> computeBridge==true) {
+ CapillaryPhys1* cundallContactPhysics=NULL;
+ MindlinCapillaryPhys* mindlinContactPhysics=NULL;
+ if (!hertzOn) cundallContactPhysics = static_cast<CapillaryPhys1*>((*ii)->phys.get());//use CapillaryPhys for linear model
+ else mindlinContactPhysics = static_cast<MindlinCapillaryPhys*>((*ii)->phys.get());//use MindlinCapillaryPhys for hertz model
+
+ if ((hertzOn && mindlinContactPhysics->meniscus) || (!hertzOn && cundallContactPhysics->meniscus)) {
+ if (fusionDetection) {//version with effect of fusion
+//BINARY VERSION : if fusionNumber!=0 then no capillary force
+ short int& fusionNumber = hertzOn?mindlinContactPhysics->fusionNumber:cundallContactPhysics->fusionNumber;
+ if (binaryFusion) {
+ if (fusionNumber!=0) { //cerr << "fusion" << endl;
+ hertzOn?mindlinContactPhysics->fCap:cundallContactPhysics->fCap = Vector3r::Zero();
+ continue;
+ }
+ }
+//LINEAR VERSION : capillary force is divided by (fusionNumber + 1) - NOTE : any decreasing function of fusionNumber can be considered in fact
+ else if (fusionNumber !=0) hertzOn?mindlinContactPhysics->fCap:cundallContactPhysics->fCap /= (fusionNumber+1.);
+ }
+ scene->forces.addForce((*ii)->getId1(), hertzOn?mindlinContactPhysics->fCap:cundallContactPhysics->fCap);
+ scene->forces.addForce((*ii)->getId2(),-(hertzOn?mindlinContactPhysics->fCap:cundallContactPhysics->fCap));
+ }
+ }
+ }
+
+}
+void Law2_ScGeom_CapillaryPhys_Capillarity1::checkFusion()
+{
+//Reset fusion numbers
+ InteractionContainer::iterator ii = scene->interactions->begin();
+ InteractionContainer::iterator iiEnd = scene->interactions->end();
+ for( ; ii!=iiEnd ; ++ii ) {
+ if ((*ii)->isReal()) {
+ if (!hertzOn) static_cast<CapillaryPhys1*>((*ii)->phys.get())->fusionNumber=0;
+ else static_cast<MindlinCapillaryPhys*>((*ii)->phys.get())->fusionNumber=0;
+ }
+ }
+
+ std::list< shared_ptr<Interaction> >::iterator firstMeniscus, lastMeniscus, currentMeniscus;
+ Real angle1 = -1.0;
+ Real angle2 = -1.0;
+
+ for ( int i=0; i< bodiesMenisciiList.size(); ++i ) { // i is the index (or id) of the body being tested
+ CapillaryPhys1* cundallInteractionPhysics1=NULL;
+ MindlinCapillaryPhys* mindlinInteractionPhysics1=NULL;
+ CapillaryPhys1* cundallInteractionPhysics2=NULL;
+ MindlinCapillaryPhys* mindlinInteractionPhysics2=NULL;
+ if ( !bodiesMenisciiList[i].empty() ) {
+ lastMeniscus = bodiesMenisciiList[i].end();
+ for ( firstMeniscus=bodiesMenisciiList[i].begin(); firstMeniscus!=lastMeniscus; ++firstMeniscus ) { //FOR EACH MENISCUS ON THIS BODY...
+ currentMeniscus = firstMeniscus;
+ ++currentMeniscus;
+ if (!hertzOn) {
+ cundallInteractionPhysics1 = YADE_CAST<CapillaryPhys1*>((*firstMeniscus)->phys.get());
+ if (i == (*firstMeniscus)->getId1()) angle1=cundallInteractionPhysics1->Delta1;//get angle of meniscus1 on body i
+ else angle1=cundallInteractionPhysics1->Delta2;
+ }
+ else {
+ mindlinInteractionPhysics1 = YADE_CAST<MindlinCapillaryPhys*>((*firstMeniscus)->phys.get());
+ if (i == (*firstMeniscus)->getId1()) angle1=mindlinInteractionPhysics1->Delta1;//get angle of meniscus1 on body i
+ else angle1=mindlinInteractionPhysics1->Delta2;
+ }
+ for ( ; currentMeniscus!= lastMeniscus; ++currentMeniscus) { //... CHECK FUSION WITH ALL OTHER MENISCII ON THE BODY
+ if (!hertzOn) {
+ cundallInteractionPhysics2 = YADE_CAST<CapillaryPhys1*>((*currentMeniscus)->phys.get());
+ if (i == (*currentMeniscus)->getId1()) angle2=cundallInteractionPhysics2->Delta1;//get angle of meniscus2 on body i
+ else angle2=cundallInteractionPhysics2->Delta2;
+ }
+ else {
+ mindlinInteractionPhysics2 = YADE_CAST<MindlinCapillaryPhys*>((*currentMeniscus)->phys.get());
+ if (i == (*currentMeniscus)->getId1()) angle2=mindlinInteractionPhysics2->Delta1;//get angle of meniscus2 on body i
+ else angle2=mindlinInteractionPhysics2->Delta2;
+ }
+ if (angle1==0 || angle2==0) cerr << "THIS SHOULD NOT HAPPEN!!"<< endl;
+
+
+ Vector3r normalFirstMeniscus = YADE_CAST<ScGeom*>((*firstMeniscus)->geom.get())->normal;
+ Vector3r normalCurrentMeniscus = YADE_CAST<ScGeom*>((*currentMeniscus)->geom.get())->normal;
+
+ Real normalDot = 0;
+ if ((*firstMeniscus)->getId1() == (*currentMeniscus)->getId1() || (*firstMeniscus)->getId2() == (*currentMeniscus)->getId2()) normalDot = normalFirstMeniscus.dot(normalCurrentMeniscus);
+ else normalDot = - (normalFirstMeniscus.dot(normalCurrentMeniscus));
+
+ Real normalAngle = 0;
+ if (normalDot >= 0 ) normalAngle = Mathr::FastInvCos0(normalDot);
+ else normalAngle = ((Mathr::PI) - Mathr::FastInvCos0(-(normalDot)));
+
+ if ((angle1+angle2)*Mathr::DEG_TO_RAD > normalAngle) {
+ if (!hertzOn) {
+ ++(cundallInteractionPhysics1->fusionNumber); //count +1 if 2 meniscii are overlaping
+ ++(cundallInteractionPhysics2->fusionNumber);
+ }
+ else {
+ ++(mindlinInteractionPhysics1->fusionNumber);
+ ++(mindlinInteractionPhysics2->fusionNumber);
+ }
+ };
+ }
+ }
+ }
+ }
+}
+
+
+BodiesMenisciiList1::BodiesMenisciiList1(Scene * scene)
+{
+ initialized=false;
+ prepare(scene);
+}
+
+bool BodiesMenisciiList1::prepare(Scene * scene)
+{
+ interactionsOnBody.clear();
+ shared_ptr<BodyContainer>& bodies = scene->bodies;
+
+ Body::id_t MaxId = -1;
+ BodyContainer::iterator bi = bodies->begin();
+ BodyContainer::iterator biEnd = bodies->end();
+ for( ; bi!=biEnd ; ++bi )
+ {
+ MaxId=max(MaxId, (*bi)->getId());
+ }
+ interactionsOnBody.resize(MaxId+1);
+ for ( unsigned int i=0; i<interactionsOnBody.size(); ++i )
+ {
+ interactionsOnBody[i].clear();
+ }
+
+ InteractionContainer::iterator ii = scene->interactions->begin();
+ InteractionContainer::iterator iiEnd = scene->interactions->end();
+ for( ; ii!=iiEnd ; ++ii ) {
+ if ((*ii)->isReal()) {
+ if (static_cast<CapillaryPhys1*>((*ii)->phys.get())->meniscus) insert(*ii);
+ }
+ }
+
+ return initialized=true;
+}
+
+bool BodiesMenisciiList1::insert(const shared_ptr< Interaction >& interaction)
+{
+ interactionsOnBody[interaction->getId1()].push_back(interaction);
+ interactionsOnBody[interaction->getId2()].push_back(interaction);
+ return true;
+}
+
+
+bool BodiesMenisciiList1::remove(const shared_ptr< Interaction >& interaction)
+{
+ interactionsOnBody[interaction->getId1()].remove(interaction);
+ interactionsOnBody[interaction->getId2()].remove(interaction);
+ return true;
+}
+
+list< shared_ptr<Interaction> >& BodiesMenisciiList1::operator[] (int index)
+{
+ return interactionsOnBody[index];
+}
+
+int BodiesMenisciiList1::size()
+{
+ return interactionsOnBody.size();
+}
+
+void BodiesMenisciiList1::display()
+{
+ list< shared_ptr<Interaction> >::iterator firstMeniscus;
+ list< shared_ptr<Interaction> >::iterator lastMeniscus;
+ for ( unsigned int i=0; i<interactionsOnBody.size(); ++i )
+ {
+ if ( !interactionsOnBody[i].empty() )
+ {
+ lastMeniscus = interactionsOnBody[i].end();
+ for ( firstMeniscus=interactionsOnBody[i].begin(); firstMeniscus!=lastMeniscus; ++firstMeniscus )
+ {
+ if ( *firstMeniscus ) {
+ if ( firstMeniscus->get() )
+ cerr << "(" << ( *firstMeniscus )->getId1() << ", " << ( *firstMeniscus )->getId2() <<") ";
+ else cerr << "(void)";
+ }
+ }
+ cerr << endl;
+ }
+ else cerr << "empty" << endl;
+ }
+}
+
+BodiesMenisciiList1::BodiesMenisciiList1()
+{
+ initialized=false;
+}
+
+void Law2_ScGeom_CapillaryPhys_Capillarity1::solver(Real suction,bool reset)
+{
+
+ if (!scene) cerr << "scene not defined!";
+ shared_ptr<BodyContainer>& bodies = scene->bodies;
+ if (dtPbased.number_of_vertices ()<1 ) triangulateData();
+ if (fusionDetection && !bodiesMenisciiList.initialized) bodiesMenisciiList.prepare(scene);
+ InteractionContainer::iterator ii = scene->interactions->begin();
+ InteractionContainer::iterator iiEnd = scene->interactions->end();
+ bool hertzInitialized = false;
+ Real i=0;
+
+ for (; ii!=iiEnd ; ++ii) {
+
+ i+=1;
+ CapillaryPhys1* phys = dynamic_cast<CapillaryPhys1*>((*ii)->phys.get());////////////////////////////////////////////////////////////////////////////////////////////
+
+
+ if ((*ii)->isReal() && phys-> computeBridge==true) {
+ const shared_ptr<Interaction>& interaction = *ii;
+ if (!hertzInitialized) {//NOTE: We are assuming that only one type is used in one simulation here
+ if (CapillaryPhys1::getClassIndexStatic()==interaction->phys->getClassIndex()) hertzOn=false;
+ else if (MindlinCapillaryPhys::getClassIndexStatic()==interaction->phys->getClassIndex()) hertzOn=true;
+ else LOG_ERROR("The capillary law is not implemented for interactions using"<<interaction->phys->getClassName());
+ }
+ hertzInitialized = true;
+ CapillaryPhys1* cundallContactPhysics=NULL;
+ MindlinCapillaryPhys* mindlinContactPhysics=NULL;
+
+/// contact physics depends on the contact law, that is used (either linear model or hertz model)
+ if (!hertzOn) cundallContactPhysics = static_cast<CapillaryPhys1*>(interaction->phys.get());//use CapillaryPhys for linear model
+ else mindlinContactPhysics = static_cast<MindlinCapillaryPhys*>(interaction->phys.get());//use MindlinCapillaryPhys for hertz model
+
+ unsigned int id1 = interaction->getId1();
+ unsigned int id2 = interaction->getId2();
+
+/// interaction geometry search (this test is to compute capillarity only between spheres (probably a better way to do that)
+ int geometryIndex1 = (*bodies)[id1]->shape->getClassIndex(); // !!!
+ int geometryIndex2 = (*bodies)[id2]->shape->getClassIndex();
+ if (!(geometryIndex1 == geometryIndex2)) continue;
+
+/// definition of interacting objects (not necessarily in contact)
+ ScGeom* currentContactGeometry = static_cast<ScGeom*>(interaction->geom.get());
+
+/// Interacting Grains:
+// If you want to define a ratio between YADE sphere size and real sphere size
+ Real alpha=1;
+ Real R1 = alpha*std::max(currentContactGeometry->radius2,currentContactGeometry->radius1);
+ Real R2 =alpha*std::min(currentContactGeometry->radius2,currentContactGeometry->radius1);
+
+ shared_ptr<BodyContainer>& bodies = scene->bodies;
+
+ Real N=bodies->size();
+ Real epsilon1,epsilon2;
+ Real ran1= (N-id1+0.5)/(N+1);
+ Real ran2= (N-id2+0.5)/(N+1);
+ epsilon1 = epsilonMean*(2*(ran1-0.5)* disp +1);//addednow
+ epsilon2 = epsilonMean*(2*(ran2-0.5)* disp +1);
+// cout << epsilon1 << "separate" <<epsilon2 <<endl;
+ R1 = R1-epsilon1*R1;
+ R2 =R2-epsilon2*R2;
+
+/// intergranular distance
+ Real D = alpha*(-(currentContactGeometry->penetrationDepth))+epsilon1*R1+epsilon2*R2;
+ if ((currentContactGeometry->penetrationDepth>=0)|| D<=0 || createDistantMeniscii) { //||(scene->iter < 1) ) // a simplified way to define meniscii everywhere
+ if (!hertzOn) {
+ if (fusionDetection && !cundallContactPhysics->meniscus) bodiesMenisciiList.insert((*ii));
+ cundallContactPhysics->meniscus=true;
+ } else {
+ if (fusionDetection && !mindlinContactPhysics->meniscus) bodiesMenisciiList.insert((*ii));
+ mindlinContactPhysics->meniscus=true;
+ }
+
+
+ }
+ Real Dinterpol = D/R1;
+
+/// Suction (Capillary pressure):
+ if (imposePressure || (!imposePressure && totalVolumeConstant)){
+ Real Pinterpol = 0;
+ if (!hertzOn) Pinterpol = cundallContactPhysics->isBroken ? 0 : suction*R1/liquidTension;
+ else Pinterpol = mindlinContactPhysics->isBroken ? 0 : suction*R1/liquidTension;
+ if (!hertzOn) cundallContactPhysics->capillaryPressure = suction;
+ else mindlinContactPhysics->capillaryPressure = suction;
+ /// Capillary solution finder:
+ if ((Pinterpol>=0) && (hertzOn? mindlinContactPhysics->meniscus : cundallContactPhysics->meniscus)) {
+ MeniscusPhysicalData solution = interpolate1(dtPbased,K::Point_3(R2/R1, Pinterpol, Dinterpol), cundallContactPhysics->m, solutions, reset);
+/// capillary adhesion force
+ Real Finterpol = solution.force;
+ Vector3r fCap = Finterpol*R1*liquidTension*currentContactGeometry->normal;
+ if (!hertzOn) cundallContactPhysics->fCap = fCap;
+ else mindlinContactPhysics->fCap = fCap;
+/// meniscus volume
+//FIXME: hardcoding numerical constants is bad practice generaly, and it probably reveals a flaw in that case (Bruno)
+ Real Vinterpol = solution.volume*pow(R1,3);
+ Real SInterface = solution.surface*pow(R1,2);
+ if (!hertzOn) {
+ cundallContactPhysics->vMeniscus = Vinterpol;
+ cundallContactPhysics->SInterface = SInterface;
+ if (Vinterpol > 0) cundallContactPhysics->meniscus = true;
+ else cundallContactPhysics->meniscus = false;
+ } else {
+ mindlinContactPhysics->vMeniscus = Vinterpol;
+ if (Vinterpol > 0) mindlinContactPhysics->meniscus = true;
+ else mindlinContactPhysics->meniscus = false;
+ }
+ if (cundallContactPhysics->meniscus== false) cundallContactPhysics->SInterface=4*3.141592653589793238462643383279502884*(pow(R1,2))+4*3.141592653589793238462643383279502884*(pow(R2,2));
+ if (!Vinterpol) {
+ if ((fusionDetection) || (hertzOn ? mindlinContactPhysics->isBroken : cundallContactPhysics->isBroken)) bodiesMenisciiList.remove((*ii));
+ if (D>((interactionDetectionFactor-1)*(currentContactGeometry->radius2+currentContactGeometry->radius1))) scene->interactions->requestErase(interaction);
+ }
+/// wetting angles
+ if (!hertzOn) {
+ cundallContactPhysics->Delta1 = max(solution.delta1,solution.delta2);
+ cundallContactPhysics->Delta2 = min(solution.delta1,solution.delta2);
+ } else {
+ mindlinContactPhysics->Delta1 = max(solution.delta1,solution.delta2);
+ mindlinContactPhysics->Delta2 = min(solution.delta1,solution.delta2);
+ }
+ }
+
+ }
+ else{
+ if (cundallContactPhysics->vMeniscus==0 and cundallContactPhysics->capillaryPressure!=0){
+ Real Pinterpol = 0;
+ if (!hertzOn) Pinterpol = cundallContactPhysics->isBroken ? 0 : cundallContactPhysics->capillaryPressure*R1/liquidTension;
+ else Pinterpol = mindlinContactPhysics->isBroken ? 0 : cundallContactPhysics->capillaryPressure*R1/liquidTension;
+// if (!hertzOn) cundallContactPhysics->capillaryPressure = suction;
+// else mindlinContactPhysics->capillaryPressure = suction;
+ /// Capillary solution finder:
+ if ((Pinterpol>=0) && (hertzOn? mindlinContactPhysics->meniscus : cundallContactPhysics->meniscus)) {
+ MeniscusPhysicalData solution = interpolate1(dtPbased,K::Point_3(R2/R1, Pinterpol, Dinterpol), cundallContactPhysics->m, solutions, reset);
+/// capillary adhesion force
+ Real Finterpol = solution.force;
+ Vector3r fCap = Finterpol*R1*liquidTension*currentContactGeometry->normal;
+ if (!hertzOn) cundallContactPhysics->fCap = fCap;
+ else mindlinContactPhysics->fCap = fCap;
+/// meniscus volume
+//FIXME: hardcoding numerical constants is bad practice generaly, and it probably reveals a flaw in that case (Bruno)
+ Real Vinterpol = solution.volume*pow(R1,3);
+ Real SInterface = solution.surface*pow(R1,2);
+ if (!hertzOn) {
+ cundallContactPhysics->vMeniscus = Vinterpol;
+ cundallContactPhysics->SInterface = SInterface;
+ if (Vinterpol > 0) cundallContactPhysics->meniscus = true;
+ else cundallContactPhysics->meniscus = false;
+ } else {
+ mindlinContactPhysics->vMeniscus = Vinterpol;
+ if (Vinterpol > 0) mindlinContactPhysics->meniscus = true;
+ else mindlinContactPhysics->meniscus = false;
+ }
+ if (cundallContactPhysics->meniscus== false) cundallContactPhysics->SInterface=4*3.141592653589793238462643383279502884*(pow(R1,2))+4*3.141592653589793238462643383279502884*(pow(R2,2));
+ if (!Vinterpol) {
+ if ((fusionDetection) || (hertzOn ? mindlinContactPhysics->isBroken : cundallContactPhysics->isBroken)) bodiesMenisciiList.remove((*ii));
+ if (D>((interactionDetectionFactor-1)*(currentContactGeometry->radius2+currentContactGeometry->radius1))) scene->interactions->requestErase(interaction);
+ }
+/// wetting angles
+ if (!hertzOn) {
+ cundallContactPhysics->Delta1 = max(solution.delta1,solution.delta2);
+ cundallContactPhysics->Delta2 = min(solution.delta1,solution.delta2);
+ } else {
+ mindlinContactPhysics->Delta1 = max(solution.delta1,solution.delta2);
+ mindlinContactPhysics->Delta2 = min(solution.delta1,solution.delta2);
+ }
+ }
+ }
+ else{
+ Real Vinterpol = 0;
+ if (!hertzOn) { Vinterpol = cundallContactPhysics->vMeniscus/pow(R1,3);
+ }
+ else Vinterpol = mindlinContactPhysics->vMeniscus;
+ /// Capillary solution finder:
+ if ((hertzOn? mindlinContactPhysics->meniscus : cundallContactPhysics->meniscus)) {
+ MeniscusPhysicalData solution = interpolate2(dtVbased,K::Point_3(R2/R1, Vinterpol, Dinterpol), cundallContactPhysics->m, solutions,reset);
+/// capillary adhesion force
+ Real Finterpol = solution.force;
+ Vector3r fCap = Finterpol*R1*liquidTension*currentContactGeometry->normal;
+ if (!hertzOn) cundallContactPhysics->fCap = fCap;
+ else mindlinContactPhysics->fCap = fCap;
+/// suction and interfacial area
+
+ Real Pinterpol = solution.succion*liquidTension/R1;
+ Real SInterface = solution.surface*pow(R1,2);
+ if (!hertzOn) {
+ cundallContactPhysics->capillaryPressure = Pinterpol;
+ cundallContactPhysics->SInterface = SInterface;
+ if (Finterpol > 0) cundallContactPhysics->meniscus = true;
+ else{
+ cundallContactPhysics->vMeniscus=0;
+ cundallContactPhysics->meniscus = false;
+ }
+ } else {
+ mindlinContactPhysics->capillaryPressure = Pinterpol;
+ if (Finterpol > 0) mindlinContactPhysics->meniscus = true;
+ else {
+ cundallContactPhysics->vMeniscus=0;
+ mindlinContactPhysics->meniscus = false;
+ }
+ }
+ if (!Vinterpol) {
+ if ((fusionDetection) || (hertzOn ? mindlinContactPhysics->isBroken : cundallContactPhysics->isBroken)) bodiesMenisciiList.remove((*ii));
+ if (D>((interactionDetectionFactor-1)*(currentContactGeometry->radius2+currentContactGeometry->radius1))) scene->interactions->requestErase(interaction);
+ }
+/// wetting angles
+ if (!hertzOn) {
+ cundallContactPhysics->Delta1 = max(solution.delta1,solution.delta2);
+ cundallContactPhysics->Delta2 = min(solution.delta1,solution.delta2);
+ } else {
+ mindlinContactPhysics->Delta1 = max(solution.delta1,solution.delta2);
+ mindlinContactPhysics->Delta2 = min(solution.delta1,solution.delta2);
+ }
+ }
+ }
+ }
+
+
+
+///interaction is not real //If the interaction is not real, it should not be in the list
+ } else if (fusionDetection) bodiesMenisciiList.remove((*ii));
+
+
+ }
+ if (fusionDetection) checkFusion();
+
+}
+
+#endif //LAW2_SCGEOM_CAPILLARYPHYS_Capillarity1
=== added file 'pkg/dem/Law2_ScGeom_CapillaryPhys_Capillarity1.hpp'
--- pkg/dem/Law2_ScGeom_CapillaryPhys_Capillarity1.hpp 1970-01-01 00:00:00 +0000
+++ pkg/dem/Law2_ScGeom_CapillaryPhys_Capillarity1.hpp 2015-07-23 17:26:21 +0000
@@ -0,0 +1,123 @@
+//
+// C++ Interface: Law2_ScGeom_CapillaryPhys_Capillarity
+/*************************************************************************
+* Copyright (C) 2006 by luc Scholtes *
+* luc.scholtes@xxxxxxxxxxx *
+* *
+* This program is free software; it is licensed under the terms of the *
+* GNU General Public License v2 or later. See file LICENSE for details. *
+*************************************************************************/
+#pragma once
+
+#include <core/GlobalEngine.hpp>
+#include <set>
+#include <boost/tuple/tuple.hpp>
+
+#include <vector>
+#include <list>
+#include <utility>
+#include <pkg/dem/CapillaryPhys1.hpp>
+#include<pkg/common/Dispatching.hpp>
+#include <string>
+#include <iostream>
+#include <fstream>
+
+/**
+This law allows one to take into account capillary forces/effects between spheres coming from the presence of interparticular liquid bridges (menisci).
+refs:
+- (french, lot of documentation) L. Scholtes, PhD thesis -> http://tel.archives-ouvertes.fr/tel-00363961/en/
+- (english, less...) L. Scholtes et al. Micromechanics of granular materials with capillary effects. International Journal of Engineering Science 2009,(47)1, 64-75
+
+The law needs ascii files M(r=i) with i=R1/R2 to work (downloaded from https://yade-dem.org/wiki/CapillaryTriaxialTest). They contain a set of results from the resolution of the Laplace-Young equation for different configurations of the interacting geometry and must be placed in the bin directory (where yade exec file is situated) to be taken into account.
+The control parameter is the capillary pressure (or suction) Delta_u, defined as the difference between gas and liquid pressure: Delta_u = u_gas - u_liquid
+Liquid bridges properties (volume V, extent over interacting grains delta1 and delta2) are computed as a result of Delta_u and the interacting geometry (spheres radii and interparticular distance)
+
+Rk: - the formulation is valid only for pendular menisci involving two grains (pendular regime).
+- an algorithm was developed by B. Chareyre to identify menisci overlaps on each spheres (menisci fusion).
+- some assumptions can be made to reduce capillary forces when menisci overlap (binary->F_cap=0 if at least 1 overlap, linear->F_cap=F_cap/numberOfOverlaps)
+*/
+
+/// !!! This version is deprecated. It should be updated to the new formalism -> ToDo !!!
+
+
+
+/// R = ratio(RadiusParticle1 on RadiusParticle2). Here, 10 R values from interpolation files (yade/extra/capillaryFiles), R = 1, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9
+//const int NB_R_VALUES = 10;
+
+// class capillarylaw1; // fait appel a la classe def plus bas
+class Interaction;
+
+
+///This container class is used to check if meniscii overlap. Wet interactions are put in a series of lists, with one list per body.
+class BodiesMenisciiList1
+{
+private:
+vector< list< shared_ptr<Interaction> > > interactionsOnBody;
+
+//shared_ptr<Interaction> empty;
+
+public:
+BodiesMenisciiList1();
+BodiesMenisciiList1(Scene* body);
+bool prepare(Scene* scene);
+bool insert(const shared_ptr<Interaction>& interaction);
+bool remove(const shared_ptr<Interaction>& interaction);
+list< shared_ptr<Interaction> >& operator[] (int index);
+int size();
+void display();
+
+
+bool initialized;
+};
+
+/// This is the constitutive law
+class Law2_ScGeom_CapillaryPhys_Capillarity1 : public GlobalEngine
+{
+public :
+ void checkFusion();
+
+ static DT dtVbased;
+ static DT dtPbased;
+ std::vector<MeniscusPhysicalData> solutions;
+ int switchTriangulation;//to detect switches between P-based and V-based data
+
+
+ BodiesMenisciiList1 bodiesMenisciiList;
+
+ void action();
+ Real intEnergy();
+ Real swInterface();
+ Real wnInterface();
+ Real waterVolume();
+ void solver(Real suction, bool reset);
+ void triangulateData();
+
+ YADE_CLASS_BASE_DOC_ATTRS_CTOR_PY(Law2_ScGeom_CapillaryPhys_Capillarity1,GlobalEngine,"This law allows one to take into account capillary forces/effects between spheres coming from the presence of interparticular liquid bridges (menisci).\n\nThe control parameter is the capillary pressure (or suction) Uc = ugas - Uliquid. Liquid bridges properties (volume V, extent over interacting grains delta1 and delta2) are computed as a result of the defined capillary pressure and of the interacting geometry (spheres radii and interparticular distance).\n\nReferences: in english [Scholtes2009b]_; more detailed, but in french [Scholtes2009d]_.\n\nThe law needs ascii files M(r=i) with i=R1/R2 to work (see https://yade-dem.org/index.php/CapillaryTriaxialTest). These ASCII files contain a set of results from the resolution of the Laplace-Young equation for different configurations of the interacting geometry."
+ "\n\nIn order to allow capillary forces between distant spheres, it is necessary to enlarge the bounding boxes using :yref:`Bo1_Sphere_Aabb::aabbEnlargeFactor` and make the Ig2 define define distant interactions via:yref:`interactionDetectionFactor<Ig2_Sphere_Sphere_ScGeom::interactionDetectionFactor>`. It is also necessary to disable interactions removal by the constitutive law (:yref:`Law2<Law2_ScGeom_FrictPhys_CundallStrack::neverErase>=True`). The only combinations of laws supported are currently capillary law + :yref:`Law2_ScGeom_FrictPhys_CundallStrack` and capillary law + :yref:`Law2_ScGeom_MindlinPhys_Mindlin` (and the other variants of Hertz-Mindlin).\n\nSee CapillaryPhys-example.py for an example script.",
+ ((Real,capillaryPressure,0.,,"Value of the capillary pressure Uc defines as Uc=Ugas-Uliquid"))
+ ((Real,totalVolumeofWater,-1.,,"Value of imposed water volume"))
+ ((Real,liquidTension,0.073,,"Value of the superficial water tension in N/m"))
+ ((Real,epsilonMean,0.,," Mean Value of the roughness"))//old value was epsilon
+ ((Real,disp,0.,," Dispersion from the mean Value of the roughness"))//added now
+ ((Real,interactionDetectionFactor,1.5,,"defines critical distance for deleting interactions. Must be consistent with the Ig2 value."))
+ ((bool,fusionDetection,false,,"If true potential menisci overlaps are checked"))
+ ((bool,initialized,false,," "))
+ ((bool,binaryFusion,true,,"If true, capillary forces are set to zero as soon as, at least, 1 overlap (menisci fusion) is detected"))
+ ((bool,hertzOn,false,,"|yupdate| true if hertz model is used"))
+ ((string,inputFilename,string("capillaryfile.txt"),,"the file with meniscus solutions, used for interpolation."))
+ ((bool,createDistantMeniscii,false,,"Generate meniscii between distant spheres? Else only maintain the existing one. For modeling a wetting path this flag should always be false. For a drying path it should be true for one step (initialization) then false, as in the logic of [Scholtes2009c]_"))
+ ((bool,imposePressure,true,," If True, suction is imposed and is constant if not Volume is imposed-Undrained test"))
+ ((bool,totalVolumeConstant,true,," in undrained test there are 2 options, If True, the total volume of water is imposed,if false the volume of each meniscus is kept constant: in this case capillary pressure can be imposed for initial distribution of meniscus or it is the total volume that can be imposed initially"))
+ ,switchTriangulation=-1,
+ .def("intEnergy",&Law2_ScGeom_CapillaryPhys_Capillarity1::intEnergy,"define the energy of interfaces in unsaturated pendular state")
+// .def("sninterface",&Law2_ScGeom_CapillaryPhys_Capillarity1::sninterface,"define the amount of solid-non-wetting interfaces in unsaturated pendular state")
+ .def("swInterface",&Law2_ScGeom_CapillaryPhys_Capillarity1::swInterface,"define the amount of solid-wetting interfaces in unsaturated pendular state")
+ .def("wnInterface",&Law2_ScGeom_CapillaryPhys_Capillarity1::wnInterface,"define the amount of wetting-non-wetiing interfaces in unsaturated pendular state")
+ .def("waterVolume",&Law2_ScGeom_CapillaryPhys_Capillarity1::waterVolume,"return the total value of water in the sample")
+ );
+};
+
+
+
+REGISTER_SERIALIZABLE(Law2_ScGeom_CapillaryPhys_Capillarity1);
+
Follow ups
