Skip to content
Snippets Groups Projects
Commit 12689d09 authored by Gauthier Becker's avatar Gauthier Becker
Browse files

NonLinearSolver --> projects. Update All CMakeLists

remove main_dgshell.cpp which is now obsolete
parent e9e4398f
Branches
Tags
No related merge requests found
Showing
with 0 additions and 2067 deletions
//
// contact Domain
//
// Description: Domain to solve contact problem
//
//
// Author: <Gauthier BECKER>, (C) 2010
//
// Copyright: See COPYING file that comes with this distribution
//
//
#include "contactDomain.h"
contactDomain::contactDomain(const contactDomain &source){
_tag = source._tag;
_phys = source._phys;
_physSlave = source._physSlave;
_penalty = source._penalty;
gSlave = source.gSlave;
gMaster = source.gMaster;
_dom = source._dom;
_bterm = source._bterm;
_massterm = source._massterm;
_lterm = source._lterm;
_contype = source._contype;
_rigid = source._rigid;
_space = source._space;
_integ = source._integ;
}
void contactDomain::setContactType(const int ct){
switch(ct){
case 0:
_contype = rigidCylinder;
_rigid = true;
break;
case 1:
_contype = rigidSphere;
_rigid = true;
break;
default:
Msg::Error("No contact know for int %d",ct);
}
}
MElement* contactDomain::getFirstElement() const {
groupOfElements::elementContainer::const_iterator it = gSlave->begin();
return (*it);
}
rigidCylinderContactDomain::rigidCylinderContactDomain(const int tag, const int physMaster, const int physSlave, const int physPt1,
const int physPt2, const double penalty,
const double h,const double rho) : contactDomain(tag,physMaster,physSlave,
penalty,rigidCylinder,true),
_thick(h), _density(rho){
// void gauss integration
_integ = new QuadratureVoid();
// creation of group of elements
gMaster = new groupOfElements(2,physMaster);
gSlave = new groupOfElements(2,physSlave);
// use for build --> no save
groupOfElements gpt1 = groupOfElements(0,physPt1);
groupOfElements gpt2 = groupOfElements(0,physPt2);
// find GC and dimension of cylinder
// get vertex
groupOfElements::vertexContainer::iterator itpt1 = gpt1.vbegin();
groupOfElements::vertexContainer::iterator itpt2 = gpt2.vbegin();
MVertex *ver1 = *itpt1;
MVertex *ver2 = *itpt2;
// Vertices coordinates
double x1 = ver1->x(); double y1 = ver1->y(); double z1 = ver1->z();
double x2 = ver2->x(); double y2 = ver2->y(); double z2 = ver2->z();
//creation of GC vertex
double xgc = 0.5*(x1+x2); double ygc = 0.5*(y1+y2); double zgc = 0.5*(z1+z2);
_vergc = new MVertex(xgc,ygc,zgc);
// dimension of cylinder
_length = ver1->distance(ver2);
// Radius search for an extreme pnt (ie the point with the most distance of gc)
double dist=-1.; // initialization
MVertex *vermax;
for(groupOfElements::vertexContainer::iterator it=gMaster->vbegin(); it!=gMaster->vend();++it){
MVertex *ver = *it;
double d = ver->distance(_vergc);
if(d > dist){
vermax = ver;
dist = d;
}
}
// radius = smallest distance of extreme point and a center
double r1 = vermax->distance(ver1);
double r2 = vermax->distance(ver2);
(r1>r2) ? _radius=r2 : _radius=r1;
// vector director of cylinder's axis
_axisDirector = new SVector3(ver1->point(),ver2->point());
_axisDirector->normalize();
}
void rigidCylinderContactDomain::initializeTerms(const unknownField *ufield){
rigidContactSpaceBase *sp = dynamic_cast<rigidContactSpaceBase*>(_space);
_massterm = new massRigidCylinder(this,sp);
_lterm = new forceRigidCylinderContact(this,sp,_thickContact,ufield);
rigidContactLinearTermBase<double> *rlterm = dynamic_cast<rigidContactLinearTermBase<double>*>(_lterm);
_bterm = new stiffnessRigidCylinderContact(sp,rlterm,ufield);
}
//
// contact Domain
//
// Description: Domain to solve contact problem
//
//
// Author: <Gauthier BECKER>, (C) 2010
//
// Copyright: See COPYING file that comes with this distribution
//
//
#ifndef CONTACTDOMAIN_H_
#define CONTACTDOMAIN_H_
#ifndef SWIG
#include "groupOfElements.h"
#include "partDomain.h"
#include "contactTerms.h"
#include "rigidCylinderContactTerms.h"
#include "MVertex.h"
#include "MElement.h"
template<class T2> class contactBilinearTermBase;
template<class T2> class contactLinearTermBase;
#endif
class contactDomain{
public:
enum contact{rigidCylinder, rigidSphere};
protected:
int _phys;
int _physSlave;
int _tag;
double _penalty;
BilinearTermBase* _bterm;
BilinearTermBase* _massterm;
LinearTermBase<double>* _lterm;
partDomain *_dom;
contact _contype;
bool _rigid;
FunctionSpaceBase *_space;
QuadratureBase *_integ;
public:
groupOfElements *gMaster;
groupOfElements *gSlave;
#ifndef SWIG
contactDomain() : gMaster(0), gSlave(0), _phys(0), _physSlave(0), _tag(0), _contype(rigidCylinder), _dom(0), _penalty(0.), _rigid(true){}
contactDomain(const int tag, const int phys, const int physSlave, double pe,
contact conty,const bool rigid=false) : _tag(tag), _phys(phys), _physSlave(physSlave),
_penalty(pe), _contype(conty),
_rigid(rigid), _space(NULL){}
contactDomain(const contactDomain &source);
~contactDomain(){}
virtual void setTag(const int t){ _tag =t;}
virtual void setPhys(const int p){_phys =p;}
virtual void setPenalty(const double pe){_penalty=pe;}
virtual void setContactType(const int ct);
virtual void setContactType(const contact ct){_contype =ct;}
virtual int getTag() const{return _tag;}
virtual int getPhys() const{return _phys;}
virtual int getPhysSlave() const{return _physSlave;}
virtual double getPenalty() const{return _penalty;}
virtual contact getContactType() const{return _contype;}
virtual partDomain* getDomain() const {return _dom;}
virtual MElement * getFirstElement() const;
virtual bool isRigid() const {return _rigid;}
virtual BilinearTermBase* getMassTerm(){return _massterm;}
virtual BilinearTermBase* getStiffnessTerm(){return _bterm;}
virtual LinearTermBase<double>* getForceTerm(){return _lterm;}
virtual FunctionSpaceBase* getSpace() {return _space;}
virtual const FunctionSpaceBase* getSpace() const {return _space;}
virtual QuadratureBase* getGaussIntegration() const{return _integ;}
virtual void setDomainAndFunctionSpace(partDomain *dom)=0;
virtual void initializeTerms(const unknownField *ufield)=0;
// static void registerBindings(binding *b);
#endif
};
class rigidCylinderContactDomain : public contactDomain{
protected:
double _length; // length of cylinder;
double _radius; // outer radius of cylinder;
double _thick; // thickness of cylinder;
double _thickContact; // use for shell (contact with neutral axis of external fiber is !=0)
double _density; // density of cylinder Not a material law for now
MVertex *_vergc; // vertex of gravity center
SVector3 *_axisDirector; // normalized vector director of cylinder's axis
#ifndef SWIG
public:
rigidCylinderContactDomain() : contactDomain(), _length(0.), _radius(0.), _thick(0.){
_contype = rigidCylinder;
_rigid = true;
_integ = new QuadratureVoid();
}
rigidCylinderContactDomain(const int tag, const int physMaster, const int physSlave, const int physpt1,
const int physpt2,const double penalty,const double h,const double rho);
~rigidCylinderContactDomain(){delete _axisDirector;}
virtual void initializeTerms(const unknownField *ufield);
SVector3* getAxisDirector() const{return _axisDirector;}
virtual void setDomainAndFunctionSpace(partDomain *dom)=0;
virtual MVertex* getGC() const{return _vergc;}
double getDensity() const{return _density;}
double getLength() const{return _length;}
double getRadius() const{return _radius;}
double getThickness() const{return _thick;}
#endif
};
#endif // CONTACTDOMAIN_H_
//
// functional space for contact
//
// Description:
//
//
// Author: <Gauthier BECKER>, (C) 2011
//
// Copyright: See COPYING file that comes with this distribution
//
#ifndef CONTACTSPACE_H_
#define CONTACTSPACE_H_
#include "functionSpace.h"
template <class T1> class ContactSpaceBase : public FunctionSpaceBase{
protected:
FunctionSpace<T1> *_space; // functional space of domain where contact is computed
const int _id;
std::vector<int> _comp;
public:
ContactSpaceBase(const int id, FunctionSpace<double> *sp) : FunctionSpaceBase(), _id(id), _space(sp){
_comp.push_back(0);
_comp.push_back(1);
_comp.push_back(2);
}
ContactSpaceBase(const int id, FunctionSpace<double> *sp, const int comp1,
const int comp2 = -1, const int comp3 =-1) : FunctionSpaceBase(), _id(id), _space(sp)
{
_comp.push_back(comp1);
if(comp2 !=-1)
_comp.push_back(comp2);
if(comp3 != -1)
_comp.push_back(comp3);
}
virtual void getKeys(MElement *ele,std::vector<Dof> &keys){
_space->getKeys(ele,keys);
}
virtual int getNumKeys(MElement *ele){
return _space->getNumKeys(ele);
}
};
class rigidContactSpaceBase : public ContactSpaceBase<double>{
protected:
const MVertex *_gc; // Node of gravity center
public:
rigidContactSpaceBase(const int id, FunctionSpace<double> *sp, MVertex *ver) : ContactSpaceBase<double>(id,sp), _gc(ver){}
rigidContactSpaceBase(const int id, FunctionSpace<double> *sp, MVertex *ver, const int comp1,
const int comp2 = -1, const int comp3 =-1) : ContactSpaceBase<double>(id,sp,comp1,comp2,comp3){}
virtual void getKeys(MElement *ele,std::vector<Dof> &keys)=0;
virtual int getNumKeys(MElement *ele)=0;
virtual void getKeysOfGravityCenter(std::vector<Dof> &keys)=0;
virtual int getNumKeysOfGravityCenter()=0;
virtual void getKeys(MElement *ele, const int ind, std::vector<Dof> &keys)=0;
virtual int getNumKeys(MElement *ele, int ind)=0; // number key in one component + number of key of GC
};
#endif // CONTACTSPACE_H_
//
// contact base term
//
// Description: contact term
//
//
// Author: <Gauthier BECKER>, (C) 2011
//
// Copyright: See COPYING file that comes with this distribution
//
#include "contactTerms.h"
#include "unknownField.h"
template<> void rigidContactLinearTermBase<double>::get(MElement *ele, int npts, IntPt *GP,fullVector<double> &m) const
{
nbdofs = _spc->getNumKeys(ele);
nbdofsGC = _spc->getNumKeysOfGravityCenter();
m.resize(nbdofs);
m.setAll(0.);
nbvertex = ele->getNumVertices();
nbcomp = (nbdofs-nbdofsGC)/nbvertex;
nbcomptimenbvertex = nbvertex*nbcomp;
// get displacement of vertices (To know the current node positions)
disp.clear(); R.clear();
_spc->getKeys(ele,R);
_ufield->get(R,disp);
verdisp[3] = disp[nbcomptimenbvertex]; verdisp[4] = disp[nbcomptimenbvertex+1]; verdisp[5] = disp[nbcomptimenbvertex+2];
for(int i=0;i<nbvertex;i++){
for(int j=0;j<6;j++) mvalue[j] =0.;
ver = ele->getVertex(i);
verdisp[0] = disp[i]; verdisp[1] = disp[i+nbvertex]; verdisp[2] = disp[i+2*nbvertex];
this->get(ver,verdisp,mvalue);
// look if there is contact force
int jj=0;
while(jj<6){
if(mvalue[jj] != 0.){
// Msg::Info("Contact detected on elem %d vertex %d",ele->getNum(),ele->getVertex(i)->getNum());
_allcontactNode->insert(ele,i,_lastNormalContact);
for(int j=0;j<nbcomp;j++){
m(i+j*nbvertex) += mvalue[j];
m(nbcomptimenbvertex+j) += mvalue[j+3]; // count two times
}
break;
}
jj++;
}
}
// if(m.norm() != 0.)
// m.print("contact");
}
//
// contact base term
//
// Description: define contact terms
//
//
// Author: <Gauthier BECKER>, (C) 2011
//
// Copyright: See COPYING file that comes with this distribution
//
//
#ifndef CONTACTTERMS_H_
#define CONTACTTERMS_H_
#include "terms.h"
#include "nodeStiffnessContact.h"
#include "contactFunctionSpace.h"
class unknownField;
template<class T1=double> class contactLinearTermBase : public LinearTermBase<T1>{
protected:
contactNodeStiffness *_allcontactNode; // need to compute BilinearTerm with more efficiency (compute only for nodes in contact)
public:
contactLinearTermBase(){
_allcontactNode = new contactNodeStiffness();
}
contactLinearTermBase(const contactLinearTermBase &source){
_allcontactNode = source._allcontactNode;
}
virtual ~contactLinearTermBase(){delete _allcontactNode;}
contactNodeStiffness* getContactNode() const {return _allcontactNode;}
void clearContactNodes(){_allcontactNode->clear();}
// virtual void get(MElement *ele, int npts, IntPt *GP, fullVector<T1> &v)=0;
};
template <class T2=double>class contactBilinearTermBase : public BilinearTermBase
{
protected:
contactLinearTermBase<T2> *_lterm; // because bilinear terms is compute only for nodes
// in contact. These nodes are in the linear terms associated to bilinearTerm
public:
contactBilinearTermBase(contactLinearTermBase<T2> *lterm) : _lterm(lterm){}
virtual ~contactBilinearTermBase(){}
// virtual void get(MElement *ele, int npts, IntPt *GP, fullMatrix<T2> &m) const=0;
// same for all contact
contactLinearTermBase<T2>* getLterm() const{return _lterm;}
contactNodeStiffness * getContactNodes() const{return _lterm->getContactNode();}
void clearContactNodes(){_lterm->clearContactNodes();}
};
template<class T2=double> class rigidContactLinearTermBase : public contactLinearTermBase<T2>{
protected:
const unknownField *_ufield;
rigidContactSpaceBase *_spc;
// data for get function (Allocated once)
mutable SVector3 _lastNormalContact; // allow to use the normal in the 2 get functions not very elegant
mutable int nbdofs,nbdofsGC, nbvertex,nbcomp,nbcomptimenbvertex;
mutable std::vector<double> disp;
mutable std::vector<Dof> R;
mutable double verdisp[6];
mutable double mvalue[6];
mutable MVertex *ver;
public:
rigidContactLinearTermBase(rigidContactSpaceBase *sp,
const unknownField *ufield) : contactLinearTermBase<T2>(),
_ufield(ufield), _spc(sp), _lastNormalContact(0.,0.,0.),
nbdofs(0), nbdofsGC(0), nbvertex(0), nbcomp(0), nbcomptimenbvertex(0)
{
}
rigidContactLinearTermBase(const rigidContactLinearTermBase<T2> &source) : contactLinearTermBase<T2>(source){
_ufield = source._ufield;
_spc = source._spc;
}
~rigidContactLinearTermBase(){}
// virtual void get(MElement *ele, fullVector<double> &m);
virtual void get(const MVertex *ver, const double disp[6],double mvalue[6]) const=0;
// same for ALL RIGID CONTACT
virtual void get(MElement *ele, int npts, IntPt *GP, fullVector<T2> &v) const;
};
#endif // CONTACTTERMS_H_
//
// contact stiffness
//
// Description: Class to manage the Dof for which the stiffness must be computed
//
//
// Author: <Gauthier BECKER>, (C) 2010
//
// Copyright: See COPYING file that comes with this distribution
//
//
#ifndef CONTACTNODESTIFFNESS_H_
#define CONTACTNODESTIFFNESS_H_
#include "MElement.h"
#include "SVector3.h"
struct dataContactStiffnessNode{
public:
MElement *_ele;
int _verIndex;
SVector3 _normal;
dataContactStiffnessNode() : _ele(NULL), _verIndex(0), _normal(0.,0.,0.){};
dataContactStiffnessNode(MElement *ele, const int ind, const SVector3 &nor) : _ele(ele), _verIndex(ind), _normal(nor){}
dataContactStiffnessNode(const dataContactStiffnessNode &source){
_ele = source._ele;
_verIndex = source._verIndex;
_normal = source._normal;
}
};
class contactNodeStiffness{
public:
typedef std::vector<dataContactStiffnessNode> nodeContainer;
protected:
nodeContainer _nodeInContact;
groupOfElements *_g; // group with elements in contact
public:
contactNodeStiffness(){_g = new groupOfElements();}
contactNodeStiffness(const contactNodeStiffness &source){
_nodeInContact = source._nodeInContact;
_g = source._g;
}
~contactNodeStiffness()
{
if(_g !=NULL) delete _g;
};
void insert(MElement *ele,const int i, SVector3 &normal){
_nodeInContact.push_back(dataContactStiffnessNode(ele,i,normal));
_g->insert(ele);
}
void clear(){
_nodeInContact.clear();
if(_g != NULL) delete _g;
_g = new groupOfElements();
}
nodeContainer::const_iterator nBegin() const{return _nodeInContact.begin();}
nodeContainer::const_iterator nEnd() const{return _nodeInContact.end();}
groupOfElements::elementContainer::const_iterator elemBegin() const{return _g->begin();}
groupOfElements::elementContainer::const_iterator elemEnd() const{return _g->end();}
};
#endif // contactNodeStiffness
//
// contact base term
//
// Description: contact with a rigid cylinder
//
//
// Author: <Gauthier BECKER>, (C) 2011
//
// Copyright: See COPYING file that comes with this distribution
//
#include "rigidCylinderContactTerms.h"
#include "contactDomain.h"
#include "unknownField.h"
massRigidCylinder::massRigidCylinder(rigidCylinderContactDomain *cdom,
rigidContactSpaceBase *spc) : _length(cdom->getLength()),_radius(cdom->getRadius()),
_thick(cdom->getThickness()),_rho(cdom->getDensity()),
_spc(spc),nbdofs(0.), masse(0.),Rint(0.){}
void massRigidCylinder::get(MElement *ele,int npts,IntPt *GP,fullMatrix<double> &m) const
{
// get the number of cylinder component (same mass in all direction)
int nbdofs = _spc->getNumKeysOfGravityCenter();
m.resize(nbdofs,nbdofs,true); // true --> setAll(0.)
if(_thick > _radius){ // full cylinder
masse = _rho*M_PI*_radius*_radius*_length;
}
else{ // hollow cylinder
Rint = _radius-_thick;
masse = _rho*M_PI*_thick*(_radius*_radius-Rint*Rint);
}
for(int i=0;i<nbdofs;i++)
m(i,i) = masse;
}
forceRigidCylinderContact::forceRigidCylinderContact(const rigidCylinderContactDomain *cdom,
rigidContactSpaceBase *sp,const double thick,
const unknownField *ufield) : rigidContactLinearTermBase(sp,ufield), _cdom(cdom),
_vergc(cdom->getGC()),
_axisDirector(cdom->getAxisDirector()),
_radius(cdom->getRadius()),
_rcontact(cdom->getRadius() + 0.*thick),
_halflength(0.5*cdom->getLength()),
_GC(cdom->getGC()->point()), d(0.),
penetration(0.), penpen(0.)
{
if(_cdom->getDomain()->IsInterfaceTerms()){
partDomain *dom = dynamic_cast<partDomain*>(_cdom->getDomain());
if(dom->getFormulation()){
_fdgfactor = 1.;
_facGC =0.5;
}
else{
_fdgfactor = 0.5;
_facGC =1.;
}
}
else{
_fdgfactor = 0.5;
_facGC = 1.;
}
_penalty = _fdgfactor*_cdom->getPenalty();
}
void forceRigidCylinderContact::get(const MVertex *ver, const double disp[6],double mvalue[6]) const
{
// B = axis' center
// A = vertex for which we look for contact
B = _GC;
Bdisp.setPosition(disp[3],disp[4],disp[5]);
B+=Bdisp;
// Distance between vertex and cylinder axis
// line BA
A = ver->point();
Adisp.setPosition(disp[0],disp[1],disp[2]);
A+=Adisp;
SVector3 BA(A,B);
dirAC = crossprod(BA,*_axisDirector);
double d = dirAC.norm();
// test to know if contact
if(d < _rcontact){
// check if contact occur before end of cylinder
dirAC.normalize();
_lastNormalContact = crossprod(dirAC,*_axisDirector);
_lastNormalContact.normalize();
C.setPosition(A,_lastNormalContact.point(),-d);
// check normal direction if C is further of B than A the direction is wrong
if(C.distance(B) > A.distance(B)){
_lastNormalContact.negate();
_lastNormalContact.normalize();
C.setPosition(A,_lastNormalContact.point(),-d);
}
if(C.distance(B) < _halflength){
penetration = _rcontact-d;
penpen = penetration*_penalty;
for(int j=0;j<3;j++){
mvalue[j]+= penpen*_lastNormalContact(j);
mvalue[j+3] -=_facGC*penpen*_lastNormalContact(j); // count two times
}
}
}
}
void stiffnessRigidCylinderContact::get(MElement *ele, int npts, IntPt *GP, fullMatrix<double> &m) const
{
int nbdofs = _spc->getNumKeys(ele);
m.resize(nbdofs,nbdofs,true); // true --> setAll(0.)
// search the vertices for which the matrix has to be computed
contactNodeStiffness * contactNodes = _lterm->getContactNode();
for(contactNodeStiffness::nodeContainer::const_iterator itn = contactNodes->nBegin(); itn!= contactNodes->nEnd(); ++itn){
if(itn->_ele == ele){ // contact
this->get(ele,itn->_verIndex,m);
}
}
}
// Protected can only be called by stiffnessRigidCylinderContact::get(MElement *ele, int npts, IntPt *GP, fullMatrix<double> &m)
// NO RESIZE
void stiffnessRigidCylinderContact::get(MElement *ele, const int verIndex, fullMatrix<double> &m) const
{
// perturbation numeric
rigidContactLinearTermBase<double>* rlterm = dynamic_cast<rigidContactLinearTermBase<double>*>(_lterm);
int nbdofselem = _spc->getNumKeys(ele) - _spc->getNumKeysOfGravityCenter();
int nbFF = ele->getNumVertices();
ver = ele->getVertex(verIndex);
disp.clear(); R.clear();
_spc->getKeys(ele,verIndex,R);
_ufield->get(R,disp);
for(int j=0;j<6;j++)
pertdisp[j] = disp[j];
// perturbation on each Dofs 3 for vertex and 3 for GC of cylinder
for(int i=0;i<6;i++){
// set force component to zero
for(int j=0;j<6;j++)
fp[j] = fm[j] = 0.;
// Force -
pertdisp[i] = disp[i] - _perturbation;
rlterm->get(ver,pertdisp,fm);
// Force +
pertdisp[i] = disp[i] + _perturbation;
rlterm->get(ver,pertdisp,fp);
// Assembly in matrix
if(i<3){ // perturbation on vertex
for(int j=0; j<3; j++){
m(verIndex+nbFF*j,verIndex+nbFF*i) -= (fp[j]-fm[j])*_doublepertexpm1;
m(nbdofselem+j,verIndex+nbFF*i) -= (fp[j+3]-fm[j+3])*_doublepertexpm1;
}
}
else{ // perturbation on gravity center
for(int j=0; j<3; j++){
m(verIndex+nbFF*j,nbdofselem+i-3) -= (fp[j]-fm[j])*_doublepertexpm1;
m(nbdofselem+j,nbdofselem+i-3) -= (fp[j+3]-fm[j+3])*_doublepertexpm1;
}
}
// virgin pertdisp
pertdisp[i] = disp[i];
}
// ANALYTIC WAY MUST BE FIXED
/* int nbdofs = _spc->getNumKeys(ele);
int nbvertex = ele->getNumVertices();
int nbdofGC = _spc->getNumKeysOfGravityCenter();
int nbcomp = (nbdofs-nbdofGC)/nbvertex;
int nbcompTimesnbVertex = nbcomp*nbvertex;
m.resize(nbdofs,nbdofs);
m.setAll(0.);
double penalty = _fdgfactor*_cdom->getPenalty();
for(int i=0;i<nbvertex;i++){
SVector3 nor = _allcontactNode->getNormal(ele,i);
if(nor.norm() != 0.){ // otherwise no contact
double dianor[3][3];
diaprod(nor,nor,dianor);
for(int j=0;j<3;j++)
for(int k=0;k<3;k++)
dianor[j][k]*=penalty;
for(int j=0;j<nbcomp;j++)
for(int k=0;k<nbcomp;k++){
m(i+j*nbvertex,i+k*nbvertex) -= dianor[j][k];
m(nbcompTimesnbVertex+j,nbcompTimesnbVertex+k) -=0.5*dianor[j][k];
}
}
}
m.print("Contact Matrix");*/
}
//
// contact base term
//
// Description: contact with a rigid cylinder
//
//
// Author: <Gauthier BECKER>, (C) 2011
//
// Copyright: See COPYING file that comes with this distribution
//
#ifndef RIGIDCONTACTCYLINDER_H_
#define RIGIDCONTACTCYLINDER_H_
#include "contactTerms.h"
#include "contactFunctionSpace.h"
class rigidCylinderContactDomain;
class massRigidCylinder : public BilinearTermBase{
protected:
double _length, _radius, _thick, _rho;
rigidContactSpaceBase *_spc;
// Data for get function (Allocated once)
mutable int nbdofs;
mutable double masse;
mutable double Rint;
public:
massRigidCylinder(rigidCylinderContactDomain *cdom,rigidContactSpaceBase *spc);
massRigidCylinder(rigidContactSpaceBase *spc,double leng, double radius,
double thick,double rho) : _spc(spc), _length(leng), _radius(radius),_thick(thick), _rho(rho){}
// Arguments are useless but it's allow to use classical assemble function
void get(MElement *ele, int npts, IntPt *GP, fullMatrix<double> &m) const;
virtual void get(MElement *ele, int npts, IntPt *GP, std::vector<fullMatrix<double> > &mv) const
{
Msg::Error("Define me get by integration point massRigidCylinder");
}
virtual BilinearTermBase* clone () const
{
return new massRigidCylinder(_spc,_length,_radius,_thick,_rho);
}
};
class forceRigidCylinderContact : public rigidContactLinearTermBase<double>{
protected:
const rigidCylinderContactDomain *_cdom;
double _fdgfactor, _facGC;
const MVertex *_vergc;
const SVector3* _axisDirector;
const double _radius;
const double _rcontact;
const double _halflength;
double _penalty;
const SPoint3 _GC;
// Data of get function (Allocated once)
mutable SPoint3 B,Bdisp,A,Adisp,C;
mutable double d,penetration,penpen;
mutable SVector3 dirAC;
public:
forceRigidCylinderContact(const rigidCylinderContactDomain *cdom,
rigidContactSpaceBase *sp,const double thick,
const unknownField *ufield);
~forceRigidCylinderContact(){}
// void get(MElement *ele,int npts, IntPt *GP,fullVector<double> &m);
// specific function to this contact type ( put in rigidContactSpaceBase ??)
// This one cannot be called by Assemble function
void get(const MVertex *ver, const double disp[6],double mvalue[6]) const;
virtual void get(MElement *ele, int npts, IntPt *GP, std::vector<fullVector<double> > &vv) const
{
Msg::Error("define me get by gauss point forceRigidCylinderContact");
}
virtual LinearTermBase<double>* clone () const
{
return new forceRigidCylinderContact(_cdom,this->_spc,_rcontact-_radius,_ufield);
}
};
class stiffnessRigidCylinderContact : public contactBilinearTermBase<double>{
private:
// can be applied only be applied by classical get function perturbation for 1 vertex
void get(MElement *ele, const int verIndex, fullMatrix<double> &m) const;
protected:
const unknownField *_ufield;
rigidContactSpaceBase *_spc;
const double _perturbation; // numerical perturbation
const double _doublepertexpm1;
// Data for get function (Allocated Once)
mutable double fp[6];
mutable double fm[6];
mutable int nbdofs;
mutable MVertex *ver;
mutable std::vector<Dof> R;
mutable std::vector<double> disp;
mutable double pertdisp[6];
public:
stiffnessRigidCylinderContact(rigidContactSpaceBase *space, contactLinearTermBase<double> *lterm, const unknownField *ufield) : contactBilinearTermBase(lterm),
_ufield(ufield),
_perturbation(1.e-10),
_doublepertexpm1(1./(2.e-10)),
_spc(space), nbdofs(0){}
~stiffnessRigidCylinderContact(){}
void get(MElement *ele, int npts, IntPt *GP, fullMatrix<double> &m) const;
virtual void get(MElement *ele, int npts, IntPt *GP, std::vector<fullMatrix<double> > &mv) const
{
Msg::Error("Define me get by integration point stiffnessRigidCylinderContact");
}
virtual BilinearTermBase* clone () const
{
return new stiffnessRigidCylinderContact(_spc,_lterm,_ufield);
}
};
#endif // RIGIDCONTACTCYLINDER_H_
# Gmsh - Copyright (C) 1997-2010 C. Geuzaine, J.-F. Remacle
#
# See the LICENSE.txt file for license information. Please report all
# bugs and problems to <gmsh@geuz.org>.
set(SRC
elementField.cpp
energyField.cpp
)
file(GLOB HDR RELATIVE ${CMAKE_CURRENT_SOURCE_DIR} *.h)
append_gmsh_src(NonLinearSolver/field "${SRC};${HDR}")
//
// C++ Interface: terms
//
// Description: Base class to manage a field on element (All field are managed by element in full Dg so displacement is also
// considered like a field on element
//
// The base class contains archiving data (which is not dependent of the field)
// Author: <Gauthier BECKER>, (C) 2010
//
// Copyright: See COPYING file that comes with this distribution
//
//
#include "elementField.h"
#include <sstream>
#include "partDomain.h"
void elementField::updateFileName(){
// add a numero to file name
std::ostringstream oss;
oss << ++numfile;
std::string s = oss.str();
std::string snum;
// cut filename and its extension
size_t ext_pos;
if(numfile==1){
ext_pos = fileName.find_last_of('.');
snum = "0";
}
else{
oss.str("");
oss << numfile-1;
snum = oss.str();
ext_pos = fileName.find_last_of(snum);
}
std::string newname(fileName,0,ext_pos-(snum.size()-1));
ext_pos = fileName.find_last_of('.');
std::string ext(fileName,ext_pos+1,fileName.size());
fileName = newname+s+"."+ext;
// create initialization of file
FILE *fp = fopen(fileName.c_str(), "w");
fprintf(fp, "$MeshFormat\n2.1 0 8\n$EndMeshFormat\n");
fclose(fp);
}
void elementField::createFile(){
if(view){
// creation of file to store displacement (delete the previous one if exist)
size_t ext_pos = fileName.find_last_of('.');
std::string newname = fileName;
newname.resize(ext_pos);
std::string rfn = "rm "+newname+"*";
system(rfn.c_str());
FILE *fp = fopen(fileName.c_str(), "w");
fprintf(fp, "$MeshFormat\n2.1 0 8\n$EndMeshFormat\n");
fclose(fp);
}
}
elementField::elementField(const std::string &fnn, const uint32_t fms, const int ncomp, const dataType dt,
const bool view_=true) : numfile(0), fileName(fnn), fmaxsize(fms), view(view_),
totelem(0), numcomp(ncomp), type(dt){
this->createFile();
switch(type){
case ElementData :
dataname = "ElementData";
break;
case ElementNodeData :
dataname = "ElementNodeData";
break;
}
}
void elementField::buildView(std::vector<partDomain*> &vdom,const double time,
const int nstep, const std::string &valuename, const int cc=-1, const bool binary=false){
if(view){
// test file size (and create a new one if needed)
uint32_t fileSize;
FILE *fp = fopen(fileName.c_str(), "r");
if(!fp){
Msg::Error("Unable to open file '%s'", fileName.c_str());
return;
}
fseek (fp, 0, SEEK_END);
fileSize = (uint32_t) (ftell(fp));
if(fileSize > fmaxsize) this->updateFileName();
fclose(fp);
fp = fopen(fileName.c_str(), "a");
// compute the number of element
if(binary) Msg::Warning("Binary write not implemented yet");
fprintf(fp, "$%s\n",dataname.c_str());
fprintf(fp, "1\n\"%s\"\n",valuename.c_str());
fprintf(fp, "1\n%.16g\n", time);
fprintf(fp, "3\n%d\n%d\n%d\n", nstep, numcomp, totelem);
std::vector<double> fieldData;
if(type == ElementNodeData){
// for (unsigned int i = 0; i < elasticFields.size(); ++i)
for(std::vector<partDomain*>::iterator itdom=vdom.begin(); itdom!=vdom.end(); ++itdom){
partDomain *dom = *itdom;
for (groupOfElements::elementContainer::const_iterator it = dom->g->begin(); it != dom->g->end(); ++it){
MElement *ele = *it;
int numv = ele->getNumVertices();
this->get(dom,ele,fieldData,cc);
fprintf(fp, "%d %d",ele->getNum(),numv);
for(int i=0;i<numv;i++)
for(int j=0;j<numcomp;j++)
fprintf(fp, " %.16g",fieldData[i+j*numv]);
fprintf(fp,"\n");
fieldData.clear();
}
}
}
else if(type == ElementData){
for (unsigned int i = 0; i < vdom.size(); ++i)
for (groupOfElements::elementContainer::const_iterator it = vdom[i]->g->begin(); it != vdom[i]->g->end(); ++it){
MElement *ele = *it;
fieldData.resize(numcomp);
this->get(vdom[i],ele,fieldData,cc);
fprintf(fp, "%d",ele->getNum());
for(int j=0;j<numcomp;j++)
fprintf(fp, " %.16g",fieldData[j]);
fprintf(fp,"\n");
}
}
fprintf(fp, "$End%s\n",dataname.c_str());
fclose(fp);
}
else Msg::Warning("No displacement view created because the variable view is set to false for this field\n");
}
//
// C++ Interface: terms
//
// Description: Base class to manage a field on element (All field are managed by element in full Dg so displacement is also
// considered like a field on element
//
// The base class contains archiving data (which is not dependent of the field)
// Author: <Gauthier BECKER>, (C) 2010
//
// Copyright: See COPYING file that comes with this distribution
//
//
#ifndef ELEMENTFIELD_H_
#define ELEMENTFIELD_H_
#include <string>
#include <stdint.h>
#include "MElement.h"
#include <vector>
#include <iostream>
#include <fstream>
class partDomain;
#include "groupOfElements.h"
class elementField{
protected :
// data needed to archive displacement
// enum for chosen is ElementData or ElementNodeData more ??
enum dataType{ElementData, ElementNodeData};
dataType type;
bool view;
std::string fileName; // name of file where the displacement are store
long int totelem; // total number of element present in all elasticFields
uint32_t fmaxsize; // Size max of file in bytes (if size of file is greater a new one is created) TODO Argument for this parameter ?
int numfile; // numero of file
int numcomp;
std::string dataname;
// function to update file name
void updateFileName();
void createFile();
public :
elementField(const std::string &fnn, const uint32_t fms, const int ncomp, const dataType dt,
const bool view_);
~elementField(){};
void setTotElem(const int ne){totelem=ne;}
virtual void get(partDomain* dom, MElement *ele,std::vector<double> &fieldData,const int comp=-1)=0; // comp allow to use an enum
// in derivate class to choose which component to save
virtual void buildView(std::vector<partDomain*> &vdom,const double time,
const int nstep, const std::string &valuename, const int cc,const bool binary);
};
#endif //
//
// C++ Interface: energetic field
//
// Description: Class derived from element field to manage energetic balance
//
//
// Author: <Gauthier BECKER>, (C) 2010
//
// Copyright: See COPYING file that comes with this distribution
//
//
#include "energyField.h"
#include "explicitHulbertChung.h"
#include "explicitDofManager.h"
#include "ipField.h"
#include "SVector3.h"
#include "unknownField.h"
#include "unknownDynamicField.h"
double dot(const std::vector<double> &a, const fullVector<typename dofManager<double>::dataMat> &b){
double c = 0.;
for(int i =0; i<a.size(); i++){
c += a[i]*b(i);
}
return c;
}
double dot(const std::vector<double> &a, const std::vector<double> &b){
double c=0.;
for(int i =0; i<a.size(); i++){
c += a[i]*b[i];
}
return c;
}
double energeticField::kineticEnergy() const {
explicitHCDofManager<double> *ddm = dynamic_cast<explicitHCDofManager<double>*>(_pAssembler);
if(ddm == NULL) return 0.; // system is not dynamic --> kinetic energy = 0
else return ddm->getKineticEnergy();
}
double energeticField::kineticEnergy(MElement *ele, const partDomain *dom) const{
explicitHulbertChung<double> * esys = dynamic_cast<explicitHulbertChung<double>*>(_lsys);
explicitHCDofManager<double> *ddm = dynamic_cast<explicitHCDofManager<double>*>(_pAssembler);
double ener=0.;
if(ddm != NULL){
// get Dof
std::vector<Dof> R;
std::vector<double> velocities;
std::vector<double> vmass;
FunctionSpaceBase *sp = dom->getFunctionSpace();
sp->getKeys(ele,R);
int nkeys = sp->getNumKeys(ele);
ddm->getDofValue(R,velocities,initialCondition::velocity);
ddm->getVertexMass(R,vmass);
for(int i=0; i!=nkeys; i++)
ener+=vmass[i]*velocities[i]*velocities[i];
R.clear(); velocities.clear(); vmass.clear();
}
return ener;
}
double energeticField::deformationEnergy(MElement *ele, const partDomain *dom) const{
return _ipf->computeDeformationEnergy(ele,dom);
}
double energeticField::deformationEnergy() const{
return _ipf->computeDeformationEnergy();
}
double energeticField::externalWork(){
double Deltawext=0.;
for(std::vector<partDomain*>::iterator itdom=_domvec.begin(); itdom!=_domvec.end();++itdom){
partDomain *dom = *itdom;
for(partDomain::neuContainer::iterator itn=dom->neuBegin(); itn!=dom->neuEnd();++itn){
nonLinearNeumannBC *neu = *itn;
FunctionSpaceBase *neuspace = neu->_space;
fullVector<typename dofManager<double>::dataMat> localVector;
std::vector<Dof> R;
std::vector<double> disp;
std::vector<double> deltau;
for(groupOfElements::elementContainer::const_iterator it = neu->g->begin(); it!=neu->g->end(); ++it){
MElement *e = *it;
R.clear(); disp.clear(); deltau.clear();
IntPt *GP;
//int npts=dom->getBulkGaussIntegrationRule()->getIntPoints(e,&GP); ?????????,
int npts = neu->integ->getIntPoints(e,&GP);
neu->_term->get(e,npts,GP,localVector);
fullVector<typename dofManager<double>::dataMat> fextnp1(localVector);
neuspace->getKeys(e,R);
_ufield->get(R,disp);
int ndofs = neuspace->getNumKeys(e);
// DeltaWext = 0.5*(Fn+1+Fn)*(un+1-un)
for(int i=0;i<ndofs;i++){
_fextnp1[R[i]] = localVector(i);
localVector(i) += _fextn[R[i]];
deltau.push_back(disp[i] - _dispn[R[i]]);
_dispnp1[R[i]] = disp[i];
}
Deltawext += dot(deltau,localVector);
}
}
}
// Now compute the work of external forces due to prescribed displacement (-Wint - Winertia)
// ATTENTION _allaef containts -fint so no minus
// Do the scalar product on vector (more efficient) TO DO THE SCALAR PRODUCT WITH BLAS HOW ??
std::vector<double> disp;
std::vector<double> acc;
std::vector<double> _mass;
std::vector<double> forceval;
// std::vector<double> finertval;
std::vector<Dof> R;
for(int i=0; i<_allaef.size(); i++){
R.push_back(_allaef[i].D);
forceval.push_back(_allaef[i].val);
}
_ufield->get(R,disp);
// Compute inertial forces if needed
explicitHCDofManager<double> *ddm = dynamic_cast<explicitHCDofManager<double>*>(_pAssembler);
if(ddm != NULL){ // otherwise static cases and no inertial forces
unknownDynamicField *duf = dynamic_cast<unknownDynamicField*>(_ufield);
duf->get(R,acc,initialCondition::acceleration);
ddm->getVertexMass(R,_mass);
for(int i=0; i<forceval.size(); i++)
forceval[i] -= _mass[i]*acc[i];
}
for(int i=0;i<R.size(); i++){
_fextnp1[R[i]] = forceval[i];
forceval[i] += _fextn[R[i]];
_dispnp1[R[i]] = disp[i];
disp[i] -= _dispn[R[i]];
}
Deltawext += dot(forceval,disp);
// swap value (next step)
_fextn.swap(_fextnp1);
_dispn.swap(_dispnp1);
// New value of Wext W = 0.5*F*u
_wext += 0.5*Deltawext;
return _wext;
}
void energeticField::get(partDomain *dom,MElement *ele,std::vector<double> &ener, const int cc){
switch(cc){
//case 0:
// ener[0] = this->kineticEnergy(ele,dom) + this->deformationEnergy(ele,dom);
// break;
case 1:
ener[0] = this->kineticEnergy(ele,dom);
break;
case 2:
ener[0] = this->deformationEnergy(ele,dom);
break;
case 3:
ener[0] = this->externalWork();
default:
ener[0] = this->kineticEnergy(ele,dom) + this->deformationEnergy(ele,dom);
}
}
void energeticField::archive(const double time){
FILE *fp = fopen(_fname.c_str(),"a");
double ekin,edefo,etot,wext;
ekin = this->kineticEnergy();
edefo= this->deformationEnergy();
wext = this->externalWork();
etot = ekin + edefo;
fprintf(fp,"%e;%e;%e;%e;%e;\n",time,ekin,edefo,wext,etot);
fclose(fp);
}
//
// C++ Interface: energetic field
//
// Description: Class derived from element field to manage energetic balance
//
//
// Author: <Gauthier BECKER>, (C) 2010
//
// Copyright: See COPYING file that comes with this distribution
//
//
#ifndef ENERGETICFIELD_H_
#define ENERGETICFIELD_H_
#include "elementField.h"
#include "nonLinearMechSolver.h"
class energeticField : public elementField{
protected:
linearSystem<double> *_lsys;
IPField *_ipf;
unknownField *_ufield;
dofManager<double> *_pAssembler;
const std::string _fname; // file to store total energy in function of time
std::vector<partDomain*> &_domvec;
std::vector<fextPrescribedDisp> &_allaef;
std::map<Dof,double> _fextn;
std::map<Dof,double> _dispn;
std::map<Dof,double> _fextnp1;
std::map<Dof,double> _dispnp1;
double _wext;
public:
energeticField(linearSystem<double> *lsys, IPField *ipf,
unknownField *ufield, dofManager<double> *pAssembler,
std::vector<partDomain*> &domvec,
std::vector<fextPrescribedDisp> &allaef) : elementField("energy.msh",100000000,1,elementField::ElementData,true),
_lsys(lsys), _ipf(ipf), _ufield(ufield), _pAssembler(pAssembler),
_fname("energy.csv"), _domvec(domvec), _allaef(allaef), _wext(0.){
// initialize file to store energy
FILE *fp = fopen(_fname.c_str(),"w");
fprintf(fp,"Time;Kinetic;Deformation;Wext;Total\n");
fprintf(fp,"0.;0.;0.;0;0.\n"); // false if initial deformation FIX IT HOW ??
fclose(fp);
long int nelem=0;
for(std::vector<partDomain*>::iterator itdom=_domvec.begin(); itdom!=_domvec.end();++itdom){
partDomain *dom = *itdom;
nelem+=dom->g->size();
for(partDomain::neuContainer::iterator itn=dom->neuBegin(); itn!=dom->neuEnd();++itn){
nonLinearNeumannBC *neu = *itn;
FunctionSpaceBase *neuspace = neu->_space;
std::vector<Dof> R;
for(groupOfElements::elementContainer::const_iterator it = neu->g->begin(); it!=neu->g->end(); ++it){
MElement *e = *it;
R.clear();
neuspace->getKeys(e,R);
int ndofs = neuspace->getNumKeys(e);
for(int i=0;i<ndofs;i++){
_fextn.insert(std::pair<Dof,double>(R[i],0.));
_dispn.insert(std::pair<Dof,double>(R[i],0.));
_fextnp1.insert(std::pair<Dof,double>(R[i],0.));
_dispnp1.insert(std::pair<Dof,double>(R[i],0.));
}
}
}
}
this->setTotElem(nelem);
// Wext for prescribed displacement
for(int i=0; i<_allaef.size(); i++){
_fextn.insert(std::pair<Dof,double>(_allaef[i].D,0.));
_dispn.insert(std::pair<Dof,double>(_allaef[i].D,0.));
_fextnp1.insert(std::pair<Dof,double>(_allaef[i].D,0.));
_dispnp1.insert(std::pair<Dof,double>(_allaef[i].D,0.));
}
}
~energeticField(){}
// functions to compute the different parts of energy
double kineticEnergy() const; // More efficient than a loop on element thanks to vector operation via PETSC
double kineticEnergy(MElement *ele, const partDomain *dom) const;
double deformationEnergy(MElement *ele, const partDomain *dom) const;
double deformationEnergy() const;
double externalWork();
void get(partDomain *dom,MElement *ele,std::vector<double> &ener, const int cc = -1);
void archive(const double time);
};
#endif // ENERGETICFIELD_H_
# Gmsh - Copyright (C) 1997-2011 C. Geuzaine, J.-F. Remacle
#
# See the LICENSE.txt file for license information. Please report all
# bugs and problems to <gmsh@geuz.org>.
set(SRC
ipField.cpp
ipstate.cpp
# inline
ipvariable.h
)
file(GLOB HDR RELATIVE ${CMAKE_CURRENT_SOURCE_DIR} *.h)
append_gmsh_src(NonLinearSolver/internalPoints "${SRC};${HDR}")
//
// C++ Interface: terms
//
// Description: Class to compute Internal point
//
//
// Author: <Gauthier BECKER>, (C) 2010
//
// Copyright: See COPYING file that comes with this distribution
//
//
#include "ipField.h"
#include "mlaw.h"
#include "ipstate.h"
#include "nonLinearMechSolver.h"
void IPField::compute1state(IPStateBase::whichState ws){
for(std::vector<partDomain*>::iterator itdom=_efield->begin(); itdom!=_efield->end(); ++itdom){
partDomain *dom = *itdom;
dom->computeIPVariable(_AIPS,_ufield,ws);
}
}
void IPField::initialBroken(MElement *iele, materialLaw* mlaw ){
Msg::Info("Interface element %d is broken at initialization",iele->getNum());
materialLaw2LawsInitializer * mflaw = dynamic_cast<materialLaw2LawsInitializer*>(mlaw);
// get ipstate
AllIPState::ipstateElementContainer *vips;
IPStateBase *ips;
vips = _AIPS->getIPstate(iele->getNum());
for(int i=0;i<vips->size();i++){
ips = (*vips)[i];
mflaw->initialBroken(ips);
}
}
void IPField::initialBroken(GModel* pModel, std::vector<int> &vnumphys){
std::vector<MVertex*> vv;
for(int i=0;i<vnumphys.size();i++){
// get the vertex associated to the physical entities
pModel->getMeshVerticesForPhysicalGroup(1,vnumphys[i],vv);
// find the InterfaceElement associated to these vertex (identify interior node as degree 2 min)
for(std::vector<partDomain*>::iterator itfield = _efield->begin(); itfield != _efield->end(); ++itfield){
dgPartDomain *dgdom = dynamic_cast<dgPartDomain*>(*itfield);
for(groupOfElements::elementContainer::const_iterator it = dgdom->gi->begin(); it!=dgdom->gi->end(); ++it)
for(int k=0;k<vv.size();k++)
if(vv[k] == (*it)->getVertex(2) ) this->initialBroken(*it, dgdom->getMaterialLaw());
}
vv.clear();
}
}
void IPField::archive(const double time){
for(std::vector<ip2archive>::iterator it=viparch.begin(); it!=viparch.end(); ++it){
ip2archive &ipa = *it;
// get ip value
double val;
switch(ipa.comp){
case -1 :
val = this->getStress(ipa.domain,ipa.ele,IPStateBase::current,ipa.numgauss,-1);
break;
case 0 :
val = this->getStress(ipa.domain,ipa.ele,IPStateBase::current,ipa.numgauss,component::xx);
break;
case 1 :
val = this->getStress(ipa.domain,ipa.ele,IPStateBase::current,ipa.numgauss,component::yy);
break;
case 2 :
val = this->getStress(ipa.domain,ipa.ele,IPStateBase::current,ipa.numgauss,component::zz);
break;
case 3 :
val = this->getStress(ipa.domain,ipa.ele,IPStateBase::current,ipa.numgauss,component::xy);
break;
case 4 :
val = this->getStress(ipa.domain,ipa.ele,IPStateBase::current,ipa.numgauss,-1);
}
FILE *fp = fopen(ipa.fname.c_str(),"a");
fprintf(fp,"%e;%e\n",time,val);
fclose(fp);
}
}
double IPField::computeDeformationEnergy(MElement *ele, const partDomain *dom) const
{
IntPt *GP;
double jac[3][3];
int npts = dom->getBulkGaussIntegrationRule()->getIntPoints(ele,&GP);
AllIPState::ipstateElementContainer *vips = _AIPS->getIPstate(ele->getNum());
double ener =0.;
for(int i=0;i<npts;i++){
IPStateBase* ips = (*vips)[i];
IPVariableMechanics *ipv = dynamic_cast<IPVariableMechanics*>(ips->getState(IPStateBase::current));
#ifdef DEBUG_
if(ipv == NULL){
Msg::Error("Compute defo energy on an non mechanics gauss' point");
return 0.;
}
#endif
double enerpt = ipv->defoEnergy();
// gauss point weight
double weight = GP[i].weight;
// Shell Only compute detJ OK
// IPVariableShell *ipvs = dynamic_cast<IPVariableShell*>(ipv);
// double detJ2 = ipvs->getLocalBasis()->getJacobian();
const double u = GP[i].pt[0]; const double v = GP[i].pt[1]; const double w = GP[i].pt[2];
const double detJ = ele->getJacobian(u, v, w, jac);
ener += weight*detJ*enerpt;
}
return ener;
}
IPField::IPField(std::vector<partDomain*> *ef,dofManager<double>* pa,
GModel *pmo, unknownField* uf, std::vector<archiveIPVariable> &vaip): _efield(ef), _dm(pa),
_ufield(uf),
elementField("stress.msh",100000000,1,
elementField::ElementData,true){
// Creation of storage for IP data
_AIPS = new AllIPState(pmo, *_efield);
// compute the number of element
long int nelem=0;
for(std::vector<partDomain*>::iterator itdom=_efield->begin(); itdom!=_efield->end(); ++itdom){
partDomain *dom = *itdom;
nelem+= dom->g->size();
}
this->setTotElem(nelem);
this->buildView(*_efield,0.,0,"VonMises",-1,false);
this->buildView(*_efield,0.,0,"sigmaxx",0,false);
this->buildView(*_efield,0.,0,"sigmayy",1,false);
this->buildView(*_efield,0.,0,"tauxy",3,false);
// Initiate file
FILE *fp = fopen("cohesiveLaw.csv","w");
fprintf(fp,"time;inter num;gauss num;sc;deltac;delta;d_n;d_t;n22;m22;n21;m21\n");
fclose(fp);
fp = fopen("crackPath.csv","w");
fprintf(fp,"time;inter num;x0;y0;z0;x1;y1;z1\n");
fclose(fp);
// Build the viparch vector (find the nearest ipvariable to each given vertex)
for(std::vector<archiveIPVariable>::iterator ita=vaip.begin(); ita!=vaip.end(); ++ita){
archiveIPVariable &aip = *ita;
// get the vertex of the node
groupOfElements g(0,aip.numphys);
groupOfElements::vertexContainer::iterator itv = g.vbegin();
MVertex *ver = *itv;
bool flagfind=false;
MElement *elefind = NULL;
partDomain *domfind;
dgPartDomain *dgdom;
bool samedim = false;
for(std::vector<partDomain*>::iterator itdom=_efield->begin(); itdom!=_efield->end(); ++itdom){
partDomain *dom = *itdom;
dgdom = dynamic_cast<dgPartDomain*>(dom);
if(dgdom != NULL){ // otherwise the domain doesn't contain interface element
for(groupOfElements::elementContainer::iterator it2 = dgdom->gi->begin(); it2 != dgdom->gi->end(); ++it2 ){
MElement *iele = *it2; // two elements on interface
int numv = iele->getNumVertices();
for(int j=0;j<numv;j++){
if(ver == iele->getVertex(j)){
elefind = iele;
domfind = dom;
flagfind = true;
break;
}
}
if(flagfind) break;
}
if(!flagfind){ // the interface element are on boundary
for(groupOfElements::elementContainer::iterator iti = dgdom->gib->begin(); iti!=dgdom->gib->end(); ++iti){
MElement *ie = *iti;
int numv = ie->getNumVertices();
for(int j=0;j<numv;j++){
if(ver == ie->getVertex(j)){
elefind = ie;
domfind = dom;
flagfind = true;
break;
}
}
if(flagfind) break;
}
}
}
if(!flagfind){ // no interface element found (Can happend on boundary)
for(groupOfElements::elementContainer::iterator it2 = dom->g->begin(); it2 != dom->g->end(); ++it2 ){
MElement *ele = *it2;
int numv = ele->getNumVertices();
for(int j=0;j<numv;j++){
if(ver == ele->getVertex(j)){
elefind = ele;
samedim = true;
flagfind = true;
domfind = dom;
break;
}
}
if(flagfind) break;
}
}
}
// Now the element where is situated the node is know find the nearest Gauss Point
IntPt *GP;
QuadratureBase *gq;
if(samedim){ //bulk point
gq = domfind->getBulkGaussIntegrationRule();
}
else{
dgPartDomain *dgfind = dynamic_cast<dgPartDomain*>(domfind);
gq = dgfind->getInterfaceGaussIntegrationRule();
}
int npts = gq->getIntPoints(elefind,&GP);
// coordonate in uvw of vertex
double uvw[3];
double xyz[3];
xyz[0] = ver->x(); xyz[1] = ver->y(); xyz[2] = ver->z();
elefind->xyz2uvw(xyz,uvw);
double distmin = 1.e100; // inf value at itnitialization
double dist, u,v,w;
int nummin;
for(int i=0;i<npts;i++){
u = GP[i].pt[0]; v = GP[i].pt[1]; w = GP[i].pt[2];
dist = sqrt((u-uvw[0])*(u-uvw[0]) + (v-uvw[1])*(v-uvw[1]) + (w-uvw[2])*(w-uvw[2]));
if(dist<distmin){
distmin = dist;
nummin = i;
}
}
// Store information
viparch.push_back(ip2archive(elefind,domfind,aip.numphys,nummin,aip.ipval));
// remove file LINUX ONLY ??
std::ostringstream oss;
oss << aip.numphys;
std::string s = oss.str();
oss.str("");
oss << aip.ipval;
std::string s2 = oss.str();
std::string rfname = "rm IP"+s+"val"+s2+".csv";
system(rfname.c_str());
}
}
double IPField::getStress(const partDomain *ef,MElement *ele, const IPStateBase::whichState ws,
const int num, const int cmp) const{
double svm =0.;
if(num<10000){ // value in a particular gauss point
IPStateBase* ips = (*_AIPS->getIPstate(ele->getNum()))[num];
IPVariableMechanics *ipv = dynamic_cast<IPVariableMechanics*>(ips->getState(ws));
#ifdef DEBUG_
if(ipv != NULL)
#endif
return ipv->stressComp(cmp);
#ifdef DEBUG_
else{
Msg::Error("Compute Von Mises stress on a non mechanical element");
return 0.;
}
#endif
}
else{ // loop on all IPVariable of an element and return a particular value (max, min, mean)
double svmp;
IntPt *GP;
int npts = ef->getBulkGaussIntegrationRule()->getIntPoints(ele,&GP);
AllIPState::ipstateElementContainer *vips = _AIPS->getIPstate(ele->getNum());
for(int i=0;i<npts;i++){
IPStateBase* ips = (*vips)[i];
IPVariableMechanics *ipv = dynamic_cast<IPVariableMechanics*>(ips->getState(ws));
#ifdef DEBUG_
if(ipv == NULL){
Msg::Error("Compute Von Mises stress on a non mechanical element");
return 0.;
}
#endif
svmp = ipv->stressComp(cmp);
if(i==0)
svm = svmp;
else{
switch(num){
case IPField::max :
if(svmp>svm) svm=svmp;
break;
case IPField::min :
if(svmp<svm) svm=svmp;
break;
case IPField::mean :
svm+=svmp;
break;
}
}
}
if(num==IPField::mean) svm/=(double)npts;
}
return svm;
}
//
// C++ Interface: terms
//
// Description: Class to compute Internal point
//
//
// Author: <Gauthier BECKER>, (C) 2010
//
// Copyright: See COPYING file that comes with this distribution
//
//
# ifndef _IPFIELD_H_
# define _IPFIELD_H_
#include<vector>
#include"quadratureRules.h"
#include "unknownField.h"
#include "elementField.h"
#include "GModel.h"
#include "ipstate.h"
#include "partDomain.h"
struct archiveIPVariable;
class IPField : public elementField {
protected :
// Struct for archiving
struct ip2archive{
MElement *ele; // element where the point is situated (if node is common to several element the first one met is chosen. Normally OK continuity)
int numgauss; // local gauss num
std::string fname; // file name where archive IP
int comp; // component to archive
const partDomain *domain; // address of domain needed to access IP
// constructor
ip2archive(MElement *e, const partDomain *dom, const int nump, const int num, const int ipval) : ele(e),
domain(dom),numgauss(num), comp(ipval){
// build file name
std::ostringstream oss;
oss << nump;
std::string s = oss.str();
oss.str("");
oss << ipval;
std::string s2 = oss.str();
fname = "IP"+s+"val"+s2+".csv";
}
ip2archive(const ip2archive &source) : ele(source.ele), domain(source.domain){
numgauss = source.numgauss;
fname = source.fname;
comp =source.comp;
}
~ip2archive(){}
};
std::vector<partDomain*>* _efield;
dofManager<double> *_dm;
AllIPState *_AIPS;
unknownField *_ufield;
std::vector<ip2archive> viparch;
// normalVariation _nvaripfield; // To avoid multiple allocation
public :
enum ElemValue{mean=10001, max=10002, min=10003}; // enum to select a particular value on element
IPField(std::vector<partDomain*> *ef,dofManager<double>* pa,
GModel *pmo, unknownField* uf, std::vector<archiveIPVariable> &vaip);
void archive(const double time);
AllIPState* getAips() const {return _AIPS;}
~IPField(){delete _AIPS;}
void compute1state(IPStateBase::whichState ws);
// On element only ?? Higher level to pass the associated element (pass edge num)
// get value with a operation
double getStress(const partDomain *ef,MElement *ele, const IPStateBase::whichState ws, const int num,
const int cmp) const;
// function to archive
virtual void get(partDomain *dom, MElement *ele, std::vector<double> &stress, const int cc=-1){
switch(cc){
case -1 :
stress[0]= this->getStress(dom,ele,IPStateBase::current,mean,-1);
break;
case 0 :
stress[0] = this->getStress(dom,ele,IPStateBase::current,mean,component::xx);
break;
case 1 :
stress[0] = this->getStress(dom,ele,IPStateBase::current,mean,component::yy);
break;
case 2 :
stress[0] = this->getStress(dom,ele,IPStateBase::current,mean,component::zz);
break;
case 3 :
stress[0] = this->getStress(dom,ele,IPStateBase::current,mean,component::xy);
break;
case 4 :
stress[0]= this->getStress(dom,ele,IPStateBase::current,max,-1);
}
}
// Interaction with Aips
void copy(IPStateBase::whichState source, IPStateBase::whichState dest){_AIPS->copy(source,dest);}
void nextStep(const double time=0.){
_AIPS->nextStep();
}
// initial broken
void initialBroken(GModel* pModel, std::vector<int> &vnumphys);
void initialBroken(MElement *iele, materialLaw *mlaw);
template<class T> void getIPv(const MElement *ele,const T** vipv, const IPStateBase::whichState ws= IPStateBase::current) const
{
AllIPState::ipstateElementContainer* vips = _AIPS->getIPstate(ele->getNum());
for(int i=0; i<vips->size(); i++){
IPStateBase* ips = (*vips)[i];
vipv[i] = dynamic_cast<const T*>(ips->getState(ws));
}
}
template<class T> void getIPv(const MElement *ele,const T** vipv, const T** vipvprev) const
{
AllIPState::ipstateElementContainer *vips = _AIPS->getIPstate(ele->getNum());
for(int i=0; i<vips->size(); i++){
IPStateBase* ips = (*vips)[i];
vipv[i] = dynamic_cast<const T*>(ips->getState(IPStateBase::current));
vipvprev[i] = dynamic_cast<const T*>(ips->getState(IPStateBase::previous));
}
}
template<class T1,class T2> void getIPv(const MInterfaceElement *iele, const T1** vipv_m, const T2** vipv_p,
const IPStateBase::whichState ws=IPStateBase::current) const
{
AllIPState::ipstateElementContainer *vips = _AIPS->getIPstate(iele->getNumber());
// SAME NUMBER OF GAUSS POINTS ON BOTH SIDES
int npts = vips->size()/2;
for(int i=0;i<npts;i++){
IPStateBase *ipsm = (*vips)[i];
IPStateBase *ipsp = (*vips)[i+npts];
vipv_m[i] = dynamic_cast<const T1*>(ipsm->getState(ws));
vipv_p[i] = dynamic_cast<const T2*>(ipsp->getState(ws));
}
}
template<class T1,class T2> void getIPv(const MInterfaceElement *iele, const T1** vipv_m, const T1** vipvprev_m,
const T2** vipv_p, const T2** vipvprev_p) const
{
AllIPState::ipstateElementContainer *vips = _AIPS->getIPstate(iele->getNumber());
// SAME NUMBER OF GAUSS POINTS ON BOTH SIDES
int npts = vips->size()/2;
for(int i=0;i<npts;i++){
IPStateBase *ipsm = (*vips)[i];
IPStateBase *ipsp = (*vips)[i+npts];
vipv_m[i] = dynamic_cast<const T1*>(ipsm->getState(IPStateBase::current));
vipv_p[i] = dynamic_cast<const T2*>(ipsp->getState(IPStateBase::current));
vipvprev_m[i] = dynamic_cast<const T1*>(ipsm->getState(IPStateBase::previous));
vipvprev_p[i] = dynamic_cast<const T2*>(ipsp->getState(IPStateBase::previous));
}
}
template<class T1,class T2> void getIPv(const MInterfaceElement *iele, T1** vipv_m, T2** vipv_p,
const IPStateBase::whichState ws=IPStateBase::current) const
{
AllIPState::ipstateElementContainer *vips = _AIPS->getIPstate(iele->getNumber());
// SAME NUMBER OF GAUSS POINTS ON BOTH SIDES
int npts = vips->size()/2;
for(int i=0;i<npts;i++){
IPStateBase *ipsm = (*vips)[i];
IPStateBase *ipsp = (*vips)[i+npts];
vipv_m[i] = dynamic_cast<T1*>(ipsm->getState(ws));
vipv_p[i] = dynamic_cast<T2*>(ipsp->getState(ws));
}
}
template<class T1,class T2> void getIPv(const MInterfaceElement *iele, T1** vipv_m, T2** vipvprev_m,
T2** vipv_p, T2** vipvprev_p) const
{
AllIPState::ipstateElementContainer *vips = _AIPS->getIPstate(iele->getNumber());
// SAME NUMBER OF GAUSS POINTS ON BOTH SIDES
int npts = vips->size()/2;
for(int i=0;i<npts;i++){
IPStateBase *ipsm = (*vips)[i];
IPStateBase *ipsp = (*vips)[i+npts];
vipv_m[i] = dynamic_cast<T1*>(ipsm->getState(IPStateBase::current));
vipv_p[i] = dynamic_cast<T2*>(ipsp->getState(IPStateBase::current));
vipvprev_m[i] = dynamic_cast<T1*>(ipsm->getState(IPStateBase::previous));
vipvprev_p[i] = dynamic_cast<T2*>(ipsp->getState(IPStateBase::previous));
}
}
template<class T1> void getIPv(const MInterfaceElement *iele, const T1** vipv_m,
const IPStateBase::whichState ws=IPStateBase::current) const
{
AllIPState::ipstateElementContainer *vips = _AIPS->getIPstate(iele->getNumber());
// SAME NUMBER OF GAUSS POINTS ON BOTH SIDES
int npts = vips->size();
for(int i=0;i<npts;i++){
IPStateBase *ipsm = (*vips)[i];
vipv_m[i] = dynamic_cast<const T1*>(ipsm->getState(ws));
}
}
template<class T1> void getIPv(const MInterfaceElement *iele, const T1** vipv_m, const T1** vipvprev_m) const
{
AllIPState::ipstateElementContainer *vips = _AIPS->getIPstate(iele->getNumber());
// SAME NUMBER OF GAUSS POINTS ON BOTH SIDES
int npts = vips->size();
for(int i=0;i<npts;i++){
IPStateBase *ipsm = (*vips)[i];
vipv_m[i] = dynamic_cast<const T1*>(ipsm->getState(IPStateBase::current));
vipvprev_m[i] = dynamic_cast<const T1*>(ipsm->getState(IPStateBase::previous));
}
}
template<class T1> void getIPv(const MInterfaceElement *iele, T1** vipv_m,
const IPStateBase::whichState ws=IPStateBase::current) const
{
AllIPState::ipstateElementContainer *vips = _AIPS->getIPstate(iele->getNumber());
int npts = vips->size();
for(int i=0;i<npts;i++){
IPStateBase *ipsm = (*vips)[i];
vipv_m[i] = dynamic_cast<T1*>(ipsm->getState(ws));
}
}
template<class T1> void getIPv(const MInterfaceElement *iele, T1** vipv_m, T1** vipvprev_m) const
{
AllIPState::ipstateElementContainer *vips = _AIPS->getIPstate(iele->getNumber());
int npts = vips->size();
for(int i=0;i<npts;i++){
IPStateBase *ipsm = (*vips)[i];
vipv_m[i] = dynamic_cast<T1*>(ipsm->getState(IPStateBase::current));
vipvprev_m[i] = dynamic_cast<T1*>(ipsm->getState(IPStateBase::previous));
}
}
double computeDeformationEnergy(MElement *ele, const partDomain *dom) const;
double computeDeformationEnergy() const{
// Msg::Error("Defo energy only for shell");
double ener=0.;
for(std::vector<partDomain*>::const_iterator itdom=_efield->begin(); itdom!=_efield->end(); ++itdom){
partDomain *dom = *itdom;
for(groupOfElements::elementContainer::const_iterator it=dom->g->begin(); it!=dom->g->end(); ++it){
MElement *ele = *it;
ener += this->computeDeformationEnergy(ele,dom);
}
}
return ener;
}
};
#endif // IPField
//
// C++ Interface: terms
//
// Description: Class to store internal variables at gauss point
//
//
// Author: <Gauthier BECKER>, (C) 2010
//
// Copyright: See COPYING file that comes with this distribution
//
//
#include "ipstate.h"
#include "partDomain.h"
#include "MInterfaceElement.h"
#include "MElement.h"
IP3State::~IP3State(){
if(_initial) delete _initial;
if(_step2) delete _step2;
if(_step1) delete _step1;
}
IP3State::IP3State(const IP3State &source) : IPStateBase()
{
_initial = source._initial;
_step1 = source._step1;
_step2 = source._step2;
_st = source._st;
}
IP3State & IP3State::operator = (const IPStateBase &source){
IPStateBase::operator=(source);
const IP3State *src = dynamic_cast<const IP3State*>(&source);
*_initial = *src->_initial;
*_step1 = *src->_step1;
*_step2 = *src->_step2;
_st = src->_st;
return *this;
}
IPVariable* IP3State::getState(const whichState wst) const
{
switch(wst){
case initial :
return _initial;
break;
case previous :
if(*_st) return _step1; else return _step2;
break;
case current :
if(*_st) return _step2; else return _step1;
break;
default : Msg::Error("Impossible to select the desired state for internal variable \n");
}
}
AllIPState::AllIPState(GModel *pModel, std::vector<partDomain*> &vdom)
{
state = true; // at creation of object state is true
IntPt *GP; // needed to know the number of gauss point
for(std::vector<partDomain*>::iterator itdom=vdom.begin(); itdom!=vdom.end(); ++itdom){
partDomain *dom = *itdom;
if(dom->IsInterfaceTerms()){
dgPartDomain *dgdom = dynamic_cast<dgPartDomain*>(dom);
materialLaw *mlawMinus = dgdom->getMaterialLawMinus();
materialLaw *mlawPlus = dgdom->getMaterialLawPlus();
// loop
for(groupOfElements::elementContainer::const_iterator it=dgdom->gi->begin(); it!=dgdom->gi->end();++it){
MElement *ele = *it;
MInterfaceElement* iele = dynamic_cast<MInterfaceElement*>(ele);
// 2* because IP is duplicated for fullDg formulation
int npts_inter=2*dgdom->getInterfaceGaussIntegrationRule()->getIntPoints(ele,&GP);
ipstateElementContainer tp(npts_inter);
for(int i=0;i<npts_inter/2;i++){
mlawMinus->createIPState(tp[i],&state,ele,iele->getElem(0)->getNumVertices());
}
for(int i=npts_inter/2;i<npts_inter;i++){
mlawPlus->createIPState(tp[i],&state,ele,iele->getElem(1)->getNumVertices());
}
_mapall.insert(ipstatePairType(iele->getNumber(),tp));
}
// Virtual interface element (no duplication)
materialLaw *mlaw = dgdom->getMaterialLaw();
for(groupOfElements::elementContainer::const_iterator it=dgdom->gib->begin(); it!=dgdom->gib->end();++it){
MElement *ele = *it;
MInterfaceElement *iele = dynamic_cast<MInterfaceElement*>(ele);
int npts_inter=dgdom->getInterfaceGaussIntegrationRule()->getIntPoints(ele,&GP);
ipstateElementContainer tp(npts_inter);
for(int i=0;i<npts_inter;i++){
mlaw->createIPState(tp[i],&state,ele,iele->getElem(0)->getNumVertices());
}
_mapall.insert(ipstatePairType(iele->getNumber(),tp));
}
}
// bulk element
materialLaw *mlaw = dom->getMaterialLaw();
for (groupOfElements::elementContainer::const_iterator it = dom->g->begin(); it != dom->g->end(); ++it){
MElement *ele = *it;
int npts_bulk=dom->getBulkGaussIntegrationRule()->getIntPoints(ele,&GP);
ipstateElementContainer tp(npts_bulk);
for(int i=0;i<npts_bulk;i++){
mlaw->createIPState(tp[i],&state,ele,ele->getNumVertices());
}
_mapall.insert(ipstatePairType(ele->getNum(),tp));
}
}
}
AllIPState::~AllIPState()
{
for(ipstateContainer::iterator it=_mapall.begin();it!=_mapall.end();++it){
ipstateElementContainer vips = (*it).second;
for(int i=0;i<vips.size();i++)
delete vips[i];
}
};
AllIPState::ipstateElementContainer* AllIPState::getIPstate(const long int num)
{
return &(_mapall.find(num)->second);
}
AllIPState::ipstateElementContainer* AllIPState::operator[](const long int num)
{
return &(_mapall.find(num)->second);
}
//
// C++ Interface: terms
//
// Description: Class to store internal variables at gauss point
//
//
// Author: <Gauthier BECKER>, (C) 2010
//
// Copyright: See COPYING file that comes with this distribution
//
//
#ifndef IPSTATE_H_
#define IPSTATE_H_
#include "GModel.h"
#include "ipvariable.h"
// enum to access to component of stress and deformation GLOBAL ??
struct component{
enum enumcomp{xx,yy,zz,xy,xz,yz};
};
class partDomain;
class dgPartDomain;
class materialLaw;
class IPStateBase{
public:
IPStateBase(){}
IPStateBase(const IPStateBase &source){}
virtual IPStateBase& operator=(const IPStateBase &source){return *this;}
enum whichState{initial, previous, current}; // keep previous OVERLOAD enum in daughter class ???
virtual IPVariable* getState(const whichState wst=IPStateBase::current) const=0;
};
class IP3State : public IPStateBase{
protected :
IPVariable *_initial; // initial state t=0
IPVariable *_step1; // previous step if _st = true and current step otherwise
IPVariable *_step2; // current step if _st = true and previous step otherwise
const bool *_st; // pointer on a bool value to choice what vector is current and what vector is previous
public :
IP3State(bool *st) : IPStateBase(),_st(st), _initial(NULL), _step1(NULL), _step2(NULL){}
IP3State(const bool *st,IPVariable *init,
IPVariable *st1, IPVariable *st2): _st(st), _initial(init), _step1(st1), _step2(st2){}
~IP3State();
IP3State(const IP3State &source);
virtual IP3State & operator = (const IPStateBase &source);
IPVariable* getState(const whichState wst=IPStateBase::current) const;
};
// Class to access to the IPState of all gauss point
class AllIPState{
public:
typedef std::map<long int, std::vector<IPStateBase*> > ipstateContainer;
typedef std::pair<long int, std::vector<IPStateBase*> > ipstatePairType;
typedef std::vector<IPStateBase*> ipstateElementContainer;
protected:
ipstateContainer _mapall;
bool state; // flag to switch previous and current (change the pointer in place of copy all variables)
public :
AllIPState(GModel *pModel, std::vector<partDomain*> &vdom);
~AllIPState();
ipstateElementContainer* getIPstate(const long int num);
ipstateElementContainer* operator[](const long int num);
void nextStep() {state ? state=false : state=true;} // change boolvalue
void copy(const IPStateBase::whichState ws1, const IPStateBase::whichState ws2){
for(ipstateContainer::iterator it=_mapall.begin(); it!=_mapall.end();++it){
std::vector<IPStateBase*> *vips = &((*it).second);
for(int i=0;i<vips->size();i++){
IPVariable *ipv1= (*vips)[i]->getState(ws1);
IPVariable *ipv2= (*vips)[i]->getState(ws2);
*ipv2 =*ipv1;
}
}
}
};
#endif // IPSTATE_H_
//
// C++ Interface: ipvariable
//
// Description: Base class for ipvariable
//
//
// Author: <Gauthier BECKER>, (C) 2011
//
// Copyright: See COPYING file that comes with this distribution
//
//
// class with the variables of IP (stress, deformation and localBasis)
#ifndef IPVARIABLE_H_
#define IPVARIABLE_H_
#include <stdlib.h>
#include <stdio.h>
class IPVariable
{
public :
IPVariable(){}
virtual ~IPVariable(){}
// copie constructor
IPVariable(const IPVariable &source){};
virtual IPVariable &operator=(const IPVariable &source){return *this;};
};
class IPVariableMechanics : public IPVariable{
public:
IPVariableMechanics(): IPVariable(){}
IPVariableMechanics(const IPVariableMechanics &source) : IPVariable(source){};
virtual IPVariableMechanics &operator=(const IPVariable &source){
IPVariable::operator=(source);
return *this;
}
virtual double defoEnergy() const{return 0;}
virtual double stressComp(const int i) const{return 0;}
};
template<class Tbulk,class Tfrac> class IPVariable2ForFracture : virtual public IPVariableMechanics{
protected:
Tbulk *_ipvbulk;
Tfrac *_ipvfrac;
bool _broken; // To know which law to use
public:
IPVariable2ForFracture() : IPVariableMechanics(), _ipvbulk(NULL), _ipvfrac(NULL), _broken(false){}
IPVariable2ForFracture(const IPVariable2ForFracture &source) : IPVariableMechanics(source), _broken(source._broken)
{
(*_ipvbulk) = (*(dynamic_cast<const IPVariable*>(source._ipvbulk)));
if(source._ipvfrac != NULL)
if(_ipvfrac == NULL) _ipvfrac = new Tfrac(*(source._ipvfrac));
else (*_ipvfrac) = (*(dynamic_cast<const IPVariable*>(source._ipvfrac)));
}
IPVariable2ForFracture& operator=(const IPVariable &source)
{
IPVariableMechanics::operator=(source);
const IPVariable2ForFracture *ipvf =dynamic_cast<const IPVariable2ForFracture *> (&source);
_broken = ipvf->_broken;
(*_ipvbulk) = (*(dynamic_cast<const IPVariableMechanics*>(ipvf->_ipvbulk)));
if(ipvf->_ipvfrac != NULL){
if(_ipvfrac == NULL) _ipvfrac = new Tfrac(*(ipvf->_ipvfrac));
else (*_ipvfrac) = (*(dynamic_cast<IPVariableMechanics*>(ipvf->_ipvfrac)));
}
else{
if(_ipvfrac != NULL) delete _ipvfrac;
_ipvfrac =NULL;
}
return *this;
}
~IPVariable2ForFracture()
{
if(_ipvbulk != NULL){
delete _ipvbulk;
}
if(_ipvfrac != NULL){
delete _ipvfrac;
}
}
const bool broken() const{return _broken;}
bool isbroken(){return _broken;}
void broken(){_broken=true;}
void nobroken(){_broken=false;}
Tbulk* getIPvBulk(){return _ipvbulk;}
const Tbulk* getIPvBulk() const{return _ipvbulk;}
Tfrac* getIPvFrac(){return _ipvfrac;}
const Tfrac* getIPvFrac() const {return _ipvfrac;}
void setIPvBulk(IPVariable *ipv){ _ipvbulk = dynamic_cast<Tbulk*>(ipv);}
void setIPvFrac(IPVariable *ipv){ _ipvfrac = dynamic_cast<Tfrac*>(ipv);}
};
#endif //IPVARIABLE_H_
# Gmsh - Copyright (C) 1997-2010 C. Geuzaine, J.-F. Remacle
#
# See the LICENSE.txt file for license information. Please report all
# bugs and problems to <gmsh@geuz.org>.
set(SRC
mlaw.cpp
)
file(GLOB HDR RELATIVE ${CMAKE_CURRENT_SOURCE_DIR} *.h)
append_gmsh_src(NonLinearSolver/materialLaw "${SRC};${HDR}")
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment