yade-users team mailing list archive
-
yade-users team
-
Mailing list archive
-
Message #00360
without contact interactions Works!
Hi,
The problems I told you are now over :-) . They were due to few mistakes
in the definition of InteractionGeometry, which was not updated in the
case of "no contact". So, everything seems ok, but, in order to prevent
existence of contact forces whithout contact, I have introduced a
condition of interaction (penetration depth (Un) > 0) in the
ElasticContactLaw. Non-existence of boxes for spheres was due to the
condition I used which called for something undefined (radius of the box
!?), but it is now ok with the penetrationDepth (for the rest, how, and
where it is calculated for Sphere/Box interaction?
InteractingBox2InteractingSphere4SphereContactGeometry?).
So, to sum up, capillary interactions and contact interactions are of
the same type (volatile). As they are both potentially active when the
boundingBoxes (1.5*radius) are interpenatred, I have added conditions
related to the penetration depth of interacting elements : Un > 0 for
the Contact Forces, and Un <= 0 for the Capillary ones (when Un > 0,
capillary forces are calculated as if Un = 0).
Modified files are attached.
If you have suggestions. I think there will be other pb when a capillary
interaction will be destructed and rebuild by the sollicitations
(meniscus parameters storage (water volume...)), but It is another question.
So long,
Luc
--
SCHOLTES Luc
Ph.D. student
luc.scholtes@xxxxxxxxxxx
Institut National Polytechnique de Grenoble
Laboratoire 3S (Soils, Solids, Structures)
BP 53 - 38041 GRENOBLE cedex 9 - FRANCE
//
// C++ Interface: CapillaryCohesiveLaw
/*************************************************************************
* Copyright (C) 2004 by Olivier Galizzi *
* olivier.galizzi@xxxxxxx *
* *
* 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. *
*************************************************************************/
#ifndef CAPILLARY_COHESIVE_LAW_HPP
#define CAPILLARY_COHESIVE_LAW_HPP
#include <yade/yade-core/InteractionSolver.hpp>
#include <set>
#include <boost/tuple/tuple.hpp>
// ajouts
#include <vector>
#include <utility>
#include <iostream>
#include <fstream>
#include <string>
using namespace std;
class Parameters
{
public :
Real V;
Real F;
Real delta1;
Real delta2;
Parameters();
Parameters(const Parameters &source); // construction par recopie
~Parameters();
} ;
const int NB_R_VALUES = 10;
class PhysicalAction;
class capillarylaw; // ???
class CapillaryCohesiveLaw : public InteractionSolver
{
private :
shared_ptr<PhysicalAction> actionForce;
shared_ptr<PhysicalAction> actionMomentum;
public :
int sdecGroupMask;
Real CapillaryPressure;
shared_ptr<capillarylaw> capillary;
CapillaryCohesiveLaw();
void action(Body* body);
protected :
void registerAttributes();
REGISTER_CLASS_NAME(CapillaryCohesiveLaw);
REGISTER_BASE_CLASS_NAME(InteractionSolver);
};
class TableauD
{
public:
Real D;
std::vector<std::vector<Real> > data;
// std::pair<Real, Real> Interpolate3(Real P);
Parameters Interpolate3(Real P);
TableauD();
TableauD(std::ifstream& file);
~TableauD();
};
// Fonction d'écriture de tableau, utilisée dans le constructeur pour test
class Tableau;
std::ostream& operator<<(std::ostream& os, Tableau& T);
class Tableau
{
public:
Real R;
std::vector<TableauD> full_data;
// std::pair<Real, Real> Interpolate2(Real D, Real P);
Parameters Interpolate2(Real D, Real P);
std::ifstream& operator<< (std::ifstream& file);
Tableau();
Tableau(const char* filename);
~Tableau();
};
class capillarylaw
{
public:
capillarylaw();
std::vector<Tableau> data_complete;
// std::pair<Real,Real> Interpolate(Real R1, Real R2, Real D, Real P);
Parameters Interpolate(Real R1, Real R2, Real D, Real P);
void fill (const char* filename);
};
REGISTER_SERIALIZABLE(CapillaryCohesiveLaw,false);
#endif // CAPILLARY_COHESIVE_LAW_HPP
/*************************************************************************
* Copyright (C) 2004 by Olivier Galizzi *
* olivier.galizzi@xxxxxxx *
* *
* 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. *
*************************************************************************/
#include "CapillaryCohesiveLaw.hpp"
#include "BodyMacroParameters.hpp"
#include "SpheresContactGeometry.hpp"
#include "SDECLinkGeometry.hpp"
#include "ElasticContactParameters.hpp"
#include "SDECLinkPhysics.hpp"
#include <yade/yade-core/Omega.hpp>
#include <yade/yade-core/MetaBody.hpp>
#include <yade/yade-package-common/Force.hpp>
#include <yade/yade-package-common/Momentum.hpp>
#include <yade/yade-core/PhysicalAction.hpp>
#include <yade/yade-lib-wm3-math/Math.hpp>
#include <iostream>
#include <fstream>
using namespace std;
CapillaryCohesiveLaw::CapillaryCohesiveLaw() : InteractionSolver() , actionForce(new Force) , actionMomentum(new Momentum)
{
sdecGroupMask=1;
capillary = shared_ptr<capillarylaw>(new capillarylaw); // ????????
//importFilename = "../data"; // chemin vers M(r=...)
capillary->fill("M(r=1)");
capillary->fill("M(r=1)");
capillary->fill("M(r=1.1)");
capillary->fill("M(r=1.25)");
capillary->fill("M(r=1.5)");
capillary->fill("M(r=1.75)");
capillary->fill("M(r=2)");
capillary->fill("M(r=3)");
capillary->fill("M(r=4)");
capillary->fill("M(r=5)");
capillary->fill("M(r=10)");
CapillaryPressure=1000;
}
void CapillaryCohesiveLaw::registerAttributes()
{
InteractionSolver::registerAttributes();
REGISTER_ATTRIBUTE(sdecGroupMask);
REGISTER_ATTRIBUTE(CapillaryPressure);
}
Parameters::Parameters()
{
V = 0;
F = 0;
delta1 = 0;
delta2 = 0;
};
Parameters::Parameters(const Parameters &source)
{
V = source.V;
F = source.F;
delta1 = source.delta1;
delta2 = source.delta2;
}
Parameters::~Parameters() {};
//FIXME : remove bool first !!!!!
void CapillaryCohesiveLaw::action(Body* body)
{
MetaBody * ncb = dynamic_cast<MetaBody*>(body);
shared_ptr<BodyContainer>& bodies = ncb->bodies;
/// Non Permanents Links ///
//InteractionContainer::iterator ii = ncb->persistentInteractions->begin();
InteractionContainer::iterator ii = ncb->volatileInteractions->begin();
//InteractionContainer::iterator iiEnd = ncb->persistentInteractions->end();
InteractionContainer::iterator iiEnd = ncb->volatileInteractions->end();
// initialisation du volume avant calcul
Real Vtotal = 0;
for( ; ii!=iiEnd ; ++ii )
{
if ((*ii)->isReal)
{
const shared_ptr<Interaction>& interaction = *ii;
unsigned int id1 = interaction->getId1();
unsigned int id2 = interaction->getId2();
if( !( (*bodies)[id1]->getGroupMask() & (*bodies)[id2]->getGroupMask() & sdecGroupMask) )
continue; // skip other groups, BTW: this is example of a good usage of 'continue' keyword
// d�inition des objets en interaction (non en contact)
BodyMacroParameters* de1 = dynamic_cast<BodyMacroParameters*>((*bodies)[id1]->physicalParameters.get());
BodyMacroParameters* de2 = dynamic_cast<BodyMacroParameters*>((*bodies)[id2]->physicalParameters.get());
SpheresContactGeometry* currentContactGeometry = dynamic_cast<SpheresContactGeometry*>(interaction->interactionGeometry.get());
//SDECLinkGeometry* currentContactGeometry = dynamic_cast<SDECLinkGeometry*>(interaction->interactionGeometry.get());
ElasticContactParameters* currentContactPhysics = dynamic_cast<ElasticContactParameters*>(interaction->interactionPhysics.get());
//SDECLinkPhysics* currentContactPhysics = dynamic_cast<SDECLinkPhysics*>(interaction->interactionPhysics.get());
/// Capillary components definition:
Real liquidTension = 0.073; // superficial tension of water N/m (20°C)
//Real teta = 0; // mouillage parfait (eau pure/billes de verre)
// Interacting Grains:
Real R1 = std::min(currentContactGeometry->radius2,currentContactGeometry->radius1);
Real R2 = std::max(currentContactGeometry->radius2,currentContactGeometry->radius1);
// intergranular distance
Real intergranularDistance = currentContactGeometry->penetrationDepth;
Real D =
(de2->se3.position-de1->se3.position).length()-(currentContactGeometry->radius1+currentContactGeometry->radius2);
if (D < 0) D = 0; // definie Fcap qd spheres interpenetrees!!
Real Dinterpol = D/R2;
//if (intergranularDistance <= 0)
//{
// Suction (Capillary pressure):
Real Pinterpol = CapillaryPressure*(R2/liquidTension);
cerr << "r = "<< R2/R1 << endl ;
cerr << "D = " << D << " | penetrationDepth = " << currentContactGeometry->penetrationDepth << endl;
/// Capillary solution finder:
if (CapillaryPressure>0)
{
// pair<Real, Real> solution = capillary->Interpolate( R1, R2, Dinterpol, Pinterpol);
Parameters solution(capillary->Interpolate( R1, R2, Dinterpol, Pinterpol));
/// capillary adhesion force
//Vector3r Fcap = 10000*currentContactGeometry->normal; //test
// Real Finterpol = solution.second;
Real Finterpol = solution.F;
Vector3r Fcap = Finterpol*(2*Mathr::PI*R2*liquidTension)*currentContactGeometry->normal;
cerr << " Fcap = " << Finterpol ;
/// meniscus volume
// Real Vinterpol = solution.first;
Real Vinterpol = solution.V;
Real Vmeniscus = Vinterpol*(R2*R2*R2);
Vtotal = (Vtotal+Vinterpol);
if (Vmeniscus != 0)
{ currentContactPhysics->meniscus = true;}
else
{ currentContactPhysics->meniscus = false;}
cerr << " | meniscus = "<< currentContactPhysics->meniscus << " | V = " << Vinterpol << endl ;
/// wetting angles
Real delta1 = max(solution.delta1,solution.delta2);
Real delta2 = min(solution.delta1,solution.delta2);
cerr << "delta1 = " << delta1 << "| delta2 = " << delta2 << endl;
static_cast<Force*> (ncb->actionParameters->find( id1 , actionForce ->getClassIndex()).get())->force += Fcap;
static_cast<Force*> (ncb->actionParameters->find( id2 , actionForce ->getClassIndex()).get())->force -= Fcap;
currentContactPhysics->prevNormal = currentContactGeometry->normal;
}
//}
}
/// Fs n'est pas modifi� par Fcap mais doit �re recalcul� (vitesse relative)?
// Vector3r x = ((de1->se3.position+de2->se3.position)*0.5)*currentContactGeometry->normal;
// Vector3r c1x = (x - de1->se3.position);
// Vector3r c2x = (x - de2->se3.position);
// Vector3r relativeVelocity = (de2->velocity+de2->angularVelocity.cross(c2x)) - (de1->velocity+de1->angularVelocity.cross(c1x));
// Vector3r shearVelocity = relativeVelocity-currentContactGeometry->normal.dot(relativeVelocity)*currentContactGeometry->normal;
// Vector3r shearDisplacement = shearVelocity*dt;
// currentContactPhysics->shearForce -= currentContactPhysics->ks*shearDisplacement;
/// modif:Fcap modifie crit�e de rupture?
// Real maxFs = (currentContactPhysics->normalForce.squaredLength() + std::pow(Mathr::PI * sigma * sqrt(currentContactGeometry->radius1*currentContactGeometry->radius2) * (c+exp( b )) ,2) * std::pow(currentContactPhysics->tangensOfFrictionAngle,2);
// maxFs= (//Fn//+ Fcap(D=0)) * (tg(frictionangle))
// if( shearForce.squaredLength() > maxFs )
// {
// maxFs = Mathr::sqRoot(maxFs) / shearForce.length();
// shearForce *= maxFs;
// }
//Vector3r f = currentContactPhysics->normalForce + currentContactPhysics->shearForce ;
}
cerr << "Vtotal = " << Vtotal << endl;
}
capillarylaw::capillarylaw() {} // ?? constructeur
void capillarylaw::fill(const char* filename)
{
data_complete.push_back(Tableau(filename));
}
// pair<Real, Real> capillarylaw::Interpolate(Real R1, Real R2, Real D, Real P)
Parameters capillarylaw::Interpolate(Real R1, Real R2, Real D, Real P)
{
if (R1 > R2) { Real R3 = R1; R1 = R2; R2 = R3;}
Real R = R2/R1;
// pair< Real, Real> result_inf;
// pair< Real, Real> result_sup;
// pair< Real, Real> result;
Parameters result_inf;
Parameters result_sup;
Parameters result;
for (int i=0; i < NB_R_VALUES; i++)
{ if (data_complete[i].R > R) // Attention à l'ordre ds lequel vont être rangés les tableau R (décroissant)
{
Tableau& tab_inf=data_complete[i-1];
Tableau& tab_sup=data_complete[i];
Real r=(R-tab_inf.R)/(tab_sup.R-tab_inf.R);
result_inf = tab_inf.Interpolate2(D,P);
result_sup = tab_sup.Interpolate2(D,P);
// result.first = result_inf.first*(1-r) + r*result_sup.first;
// result.second = result_inf.second*(1-r) + r*result_sup.second;
result.V = result_inf.V*(1-r) + r*result_sup.V;
result.F = result_inf.F*(1-r) + r*result_sup.F;
result.delta1 = result_inf.delta1*(1-r) + r*result_sup.delta1;
result.delta2 = result_inf.delta2*(1-r) + r*result_sup.delta2;
i=NB_R_VALUES;
}
else if (data_complete[i].R == R)
{
result = data_complete[i].Interpolate2(D,P);
i=NB_R_VALUES;
}
}
return result;
}
Tableau::Tableau() {} // ?? constructeur
Tableau::Tableau(const char* filename)
{
ifstream file (filename);
file >> R;
//cout << R << endl;
int n_D; //number of D values
file >> n_D;
if (!file.is_open()) cout << "problem opening file for capillary law" << endl;
for (int i=0; i<n_D; i++) full_data.push_back(TableauD(file));
file.close();
//cerr << *this; // exemple d'utilisation de la fonction d'écriture (this est le pointeur vers l'objet courant)
}
Tableau::~Tableau() {} // ?? destructeur
// pair<Real, Real> Tableau::Interpolate2(Real D, Real P)
Parameters Tableau::Interpolate2(Real D, Real P)
{
// pair< Real, Real> result;
// pair< Real, Real> result_inf;
// pair< Real, Real> result_sup;
Parameters result;
Parameters result_inf;
Parameters result_sup;
for ( int i=0; i < full_data.size(); ++i)
{ if (full_data[i].D > D ) // ok si D rangés ds l'ordre croissant
{
Real rD = (D-full_data[i-1].D)/(full_data[i].D-full_data[i-1].D);
result_inf = full_data[i-1].Interpolate3(P);
result_sup = full_data[i].Interpolate3(P);
// result.first = result_inf.first*(1-rD) + rD*result_sup.first;
// result.second = result_inf.second*(1-rD) + rD*result_sup.second;
result.V = result_inf.V*(1-rD) + rD*result_sup.V;
result.F = result_inf.F*(1-rD) + rD*result_sup.F;
result.delta1 = result_inf.delta1*(1-rD) + rD*result_sup.delta1;
result.delta2 = result_inf.delta2*(1-rD) + rD*result_sup.delta2;
i = full_data.size();
}
else if (full_data[i].D == D)
{
result=full_data[i].Interpolate3(P);
i=full_data.size();
}
}
return result;
}
TableauD::TableauD() {} // ?? constructeur
TableauD::TableauD(ifstream& file)
{
int i=0;
Real x;
int n_lines; //pb: n_lines is real!!!
file >> n_lines;
//cout << n_lines << endl;
file.ignore(200, '\n'); // saute les caractères (200 au maximum) jusque au caractère \n (fin de ligne)*_
if (n_lines!=0) for (i; i<n_lines; ++i)
{
data.push_back(vector<Real> ());
for (int j=0; j < 6; ++j) // [D,P,V,F,delta1,delta2]
{
file >> x;
data[i].push_back(x);
}
}
D = data[i-1][0];
}
// pair<Real, Real> TableauD::Interpolate3(Real P)
Parameters TableauD::Interpolate3(Real P)
{
// pair<Real, Real> result;
Parameters result;
for (int k=0; k < data.size(); ++k) // length(data) ??
{
if ( data[k][1] > P) // OK si P rangés ds l'ordre croissant
{
Real Pinf=data[k-1][1];
Real Finf=data[k-1][3];
Real Vinf=data[k-1][2];
Real Delta1inf=data[k-1][4];
Real Delta2inf=data[k-1][5];
Real Psup=data[k][1];
Real Fsup=data[k][3];
Real Vsup=data[k][2];
Real Delta1sup=data[k-1][4];
Real Delta2sup=data[k-1][5];
// result.first = Vinf+((Vsup-Vinf)/(Psup-Pinf))*(P-Pinf);
// result.second = Finf+((Fsup-Finf)/(Psup-Pinf))*(P-Pinf);
result.V = Vinf+((Vsup-Vinf)/(Psup-Pinf))*(P-Pinf);
result.F = Finf+((Fsup-Finf)/(Psup-Pinf))*(P-Pinf);
result.delta1 = Delta1inf+((Delta1sup-Delta1inf)/(Psup-Pinf))*(P-Pinf);
result.delta2 = Delta2inf+((Delta2sup-Delta2inf)/(Psup-Pinf))*(P-Pinf);
k=data.size();
}
else if (data[k][1]==P)
{
// result.first=data[k][2];
// result.second=data[k][3];
result.V = data[k][2];
result.F = data[k][3];
result.delta1 = data[k][4];
result.delta2 = data[k][5];
k=data.size();
}
}
return result;
}
TableauD::~TableauD() {} // ?? destructeur
std::ostream& operator<<(std::ostream& os, Tableau& T)
{
os << "Tableau : R=" << T.R << endl;
for (int i=0; i<T.full_data.size(); i++)
{
os << "TableauD : D=" << T.full_data[i].D << endl;
for (int j=0; j<T.full_data[i].data.size();j++)
{
for (int k=0; k<T.full_data[i].data[j].size(); k++)
os << T.full_data[i].data[j][k] << " ";
os << endl;
}
}
os << endl;
return os;
}
/*************************************************************************
* Copyright (C) 2004 by Olivier Galizzi *
* olivier.galizzi@xxxxxxx *
* *
* 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. *
*************************************************************************/
#include "ElasticContactLaw.hpp"
#include "BodyMacroParameters.hpp"
#include "SpheresContactGeometry.hpp"
#include "SDECLinkGeometry.hpp"
#include "ElasticContactParameters.hpp"
#include "SDECLinkPhysics.hpp"
#include <yade/yade-core/Omega.hpp>
#include <yade/yade-core/MetaBody.hpp>
#include <yade/yade-package-common/Force.hpp>
#include <yade/yade-package-common/Momentum.hpp>
#include <yade/yade-core/PhysicalAction.hpp>
ElasticContactLaw::ElasticContactLaw() : InteractionSolver() , actionForce(new Force) , actionMomentum(new Momentum)
{
sdecGroupMask=1;
momentRotationLaw = true;
}
void ElasticContactLaw::registerAttributes()
{
InteractionSolver::registerAttributes();
REGISTER_ATTRIBUTE(sdecGroupMask);
REGISTER_ATTRIBUTE(momentRotationLaw);
}
//FIXME : remove bool first !!!!!
void ElasticContactLaw::action(Body* body)
{
MetaBody * ncb = dynamic_cast<MetaBody*>(body);
shared_ptr<BodyContainer>& bodies = ncb->bodies;
Real dt = Omega::instance().getTimeStep();
/// Non Permanents Links ///
InteractionContainer::iterator ii = ncb->volatileInteractions->begin();
InteractionContainer::iterator iiEnd = ncb->volatileInteractions->end();
for( ; ii!=iiEnd ; ++ii )
{
if ((*ii)->isReal)
{ //cerr << "*ii is Real" << endl;
const shared_ptr<Interaction>& contact = *ii;
//cerr << "contact = " << contact;
int id1 = contact->getId1();
int id2 = contact->getId2();
if( !( (*bodies)[id1]->getGroupMask() & (*bodies)[id2]->getGroupMask() & sdecGroupMask) )
continue; // skip other groups, BTW: this is example of a good usage of 'continue' keyword
BodyMacroParameters* de1 = dynamic_cast<BodyMacroParameters*>((*bodies)[id1]->physicalParameters.get());
BodyMacroParameters* de2 = dynamic_cast<BodyMacroParameters*>((*bodies)[id2]->physicalParameters.get());
SpheresContactGeometry* currentContactGeometry = dynamic_cast<SpheresContactGeometry*>(contact->interactionGeometry.get());
ElasticContactParameters* currentContactPhysics = dynamic_cast<ElasticContactParameters*> (contact->interactionPhysics.get());
Vector3r& shearForce = currentContactPhysics->shearForce;
/// test
// Real D = (currentContactGeometry->radius1+currentContactGeometry->radius2)-(de2->se3.position-de1->se3.position).length();
// //
// Real r = std::max(currentContactGeometry->radius2,currentContactGeometry->radius1)/std::min(currentContactGeometry->radius2,currentContactGeometry->radius1);
// cerr << " R1 = " << std::min(currentContactGeometry->radius2,currentContactGeometry->radius1);
// cerr << " R2 = " << std::max(currentContactGeometry->radius2,currentContactGeometry->radius1) << endl;
//
// cerr << "r = " << r;
// cerr << " D = " << D << endl;
//
// if (D >= 0) /// c'est cette condition qui pose pb. la rendre active pour les spheres seulement?
// {
Real un = currentContactGeometry->penetrationDepth;
//cerr << " penetrationdepth = " << un;
//cerr << " Fn = " << currentContactPhysics->normalForce << endl;
if (un > 0) //new test
{
cerr << " contact!" << endl;
//ancien prog
if ( contact->isNew)
shearForce = Vector3r(0,0,0);
currentContactPhysics->normalForce = currentContactPhysics->kn*un*currentContactGeometry->normal;
//cerr << " Fn = " << currentContactPhysics->normalForce << endl;
Vector3r axis;
Real angle;
/// Here is the code with approximated rotations ///
axis = currentContactPhysics->prevNormal.cross(currentContactGeometry->normal);
shearForce -= shearForce.cross(axis);
angle = dt*0.5*currentContactGeometry->normal.dot(de1->angularVelocity+de2->angularVelocity);
axis = angle*currentContactGeometry->normal;
shearForce -= shearForce.cross(axis);
/// Here is the code with exact rotations ///
// Quaternionr q;
//
// axis = currentContactPhysics->prevNormal.cross(currentContactGeometry->normal);
// angle = acos(currentContactGeometry->normal.dot(currentContactPhysics->prevNormal));
// q.fromAngleAxis(angle,axis);
//
// currentContactPhysics->shearForce = currentContactPhysics->shearForce*q;
//
// angle = dt*0.5*currentContactGeometry->normal.dot(de1->angularVelocity+de2->angularVelocity);
// axis = currentContactGeometry->normal;
// q.fromAngleAxis(angle,axis);
// currentContactPhysics->shearForce = q*currentContactPhysics->shearForce;
/// ///
Vector3r x = currentContactGeometry->contactPoint;
Vector3r c1x = (x - de1->se3.position);
Vector3r c2x = (x - de2->se3.position);
Vector3r relativeVelocity = (de2->velocity+de2->angularVelocity.cross(c2x)) - (de1->velocity+de1->angularVelocity.cross(c1x));
Vector3r shearVelocity = relativeVelocity-currentContactGeometry->normal.dot(relativeVelocity)*currentContactGeometry->normal;
Vector3r shearDisplacement = shearVelocity*dt;
shearForce -= currentContactPhysics->ks*shearDisplacement;
// PFC3d SlipModel, is using friction angle. CoulombCriterion
Real maxFs = currentContactPhysics->normalForce.squaredLength() * std::pow(currentContactPhysics->tangensOfFrictionAngle,2);
if( shearForce.squaredLength() > maxFs )
{
maxFs = Mathr::sqRoot(maxFs) / shearForce.length();
shearForce *= maxFs;
}
////////// PFC3d SlipModel
Vector3r f = currentContactPhysics->normalForce + shearForce;
// it will be some macro( body->actionParameters, ActionType , bodyId )
static_cast<Force*> ( ncb->actionParameters->find( id1 , actionForce ->getClassIndex() ).get() )->force -= f;
static_cast<Force*> ( ncb->actionParameters->find( id2 , actionForce ->getClassIndex() ).get() )->force += f;
static_cast<Momentum*>( ncb->actionParameters->find( id1 , actionMomentum->getClassIndex() ).get() )->momentum -= c1x.cross(f);
static_cast<Momentum*>( ncb->actionParameters->find( id2 , actionMomentum->getClassIndex() ).get() )->momentum += c2x.cross(f);
currentContactPhysics->prevNormal = currentContactGeometry->normal;
// fin de nouvelle boucle
}
}
}
}
/*************************************************************************
* Copyright (C) 2004 by Olivier Galizzi *
* olivier.galizzi@xxxxxxx *
* Copyright (C) 2004 by Janek Kozicki *
* cosurgi@xxxxxxxxxx *
* *
* 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. *
*************************************************************************/
#include "InteractingSphere2InteractingSphere4SpheresContactGeometryWater.hpp"
#include "SpheresContactGeometry.hpp"
#include "SDECLinkGeometry.hpp"
#include <yade/yade-package-common/InteractingSphere.hpp>
bool InteractingSphere2InteractingSphere4SpheresContactGeometryWater::go( const shared_ptr<InteractingGeometry>& cm1,
const shared_ptr<InteractingGeometry>& cm2,
const Se3r& se31,
const Se3r& se32,
const shared_ptr<Interaction>& c)
{
InteractingSphere* s1 = static_cast<InteractingSphere*>(cm1.get());
InteractingSphere* s2 = static_cast<InteractingSphere*>(cm2.get());
/// c'est ce calcul qui semble poser pb !!!
Vector3r normal = se32.position-se31.position;
Real penetrationDepth = s1->radius+s2->radius-normal.normalize();
shared_ptr<SpheresContactGeometry> scm;
shared_ptr<SDECLinkGeometry> linkGeometry;
if (c->interactionGeometry)
{
scm = dynamic_pointer_cast<SpheresContactGeometry>(c->interactionGeometry);
}
// BEGIN ....... FIXME FIXME - wrong hack, to make cohesion work.
else
scm = shared_ptr<SpheresContactGeometry>(new SpheresContactGeometry());
//linkGeometry = shared_ptr<SDECLinkGeometry>(new SDECLinkGeometry());
if (penetrationDepth>0)
{
scm->contactPoint = se31.position+(s1->radius-0.5*penetrationDepth)*normal;//0.5*(pt1+pt2);
scm->normal = normal;
scm->penetrationDepth = penetrationDepth;
scm->radius1 = s1->radius;
scm->radius2 = s2->radius;
// ?????
if (!c->interactionGeometry)
c->interactionGeometry = scm;
return true;
}
else //définition de géométrie d'interaction pour le cas ou !contact
scm->normal = normal;
scm->penetrationDepth = penetrationDepth;
scm->radius1 = s1->radius;
scm->radius2 = s2->radius;
//linkGeometry->normal.normalize();
//linkGeometry->radius1 = s1->radius;
//linkGeometry->radius2 = s2->radius;
if (!c->interactionGeometry)
c->interactionGeometry = scm;
//if (!c->interactionGeometry)
// c->interactionGeometry = linkGeometry;
return true;
//return false;
}
bool InteractingSphere2InteractingSphere4SpheresContactGeometryWater::goReverse( const shared_ptr<InteractingGeometry>& cm1,
const shared_ptr<InteractingGeometry>& cm2,
const Se3r& se31,
const Se3r& se32,
const shared_ptr<Interaction>& c)
{
return go(cm1,cm2,se31,se32,c);
}
_______________________________________________
Yade-users mailing list
Yade-users@xxxxxxxxxxxxxxxx
https://lists.berlios.de/mailman/listinfo/yade-users
