Select Git revision
meshGFace.cpp
Forked from
gmsh / gmsh
Source project has a limited visibility.
dG3DMaterialLaw.cpp 125.94 KiB
//
// C++ Interface: materialLaw
//
// Description: Class with definition of materialLaw for shell
//
//
// Author: <Gauthier BECKER>, (C) 2011
//
// Copyright: See COPYING file that comes with this distribution
//
//
#include "dG3DMaterialLaw.h"
#include "ipstate.h"
#include "MInterfaceElement.h"
#include "ipField.h"
#include "mlawLinearThermoMechanics.h"
#include "mlawSMP.h"
#include "mlawPhenomenologicalSMP.h"
#include "mlaw.h"
#include "mlawLinearElecTherMech.h"
#include "mlawAnIsotropicElecTherMech.h"
#include "mlawElecSMP.h"
#include <iostream>
#include <fstream>
#include <iomanip>
#include <vector>
#include <string>
#include "dG3DCohesiveIPVariable.h"
#include "FractureCohesiveDG3DIPVariable.h"
void dG3DMaterialLaw::fillElasticStiffness(double E, double nu, STensor43 &K_) const
{
double mu = 0.5*E/(1.+nu);
double lambda = (E*nu)/(1.+nu)/(1.-2.*nu);
double twicemu = mu+mu;
K_*=0.;
K_(0,0,0,0) = lambda + twicemu;
K_(1,1,0,0) = lambda;
K_(2,2,0,0) = lambda;
K_(0,0,1,1) = lambda;
K_(1,1,1,1) = lambda + twicemu;
K_(2,2,1,1) = lambda;
K_(0,0,2,2) = lambda;
K_(1,1,2,2) = lambda;
K_(2,2,2,2) = lambda + twicemu;
if(lambda>=1000.*mu)
{
mu = lambda + mu;
}
K_(1,0,1,0) = mu;
K_(2,0,2,0) = mu;
K_(0,1,0,1) = mu;
K_(2,1,2,1) = mu;
K_(0,2,0,2) = mu;
K_(1,2,1,2) = mu;
K_(0,1,1,0) = mu;
K_(0,2,2,0) = mu;
K_(1,0,0,1) = mu;
K_(1,2,2,1) = mu;
K_(2,0,0,2) = mu;
K_(2,1,1,2) = mu;
}
void dG3DMaterialLaw::initialIPVariable(IPVariable* ipv, bool stiff){ IPVariable* ipvTemp = ipv->clone();
this->stress(ipv,ipvTemp,stiff,false);
delete ipvTemp;
};
dG3DLinearElasticMaterialLaw::dG3DLinearElasticMaterialLaw(const int num, const double rho):
dG3DMaterialLaw(num,rho,false),_elasticLaw(NULL){}
dG3DLinearElasticMaterialLaw::dG3DLinearElasticMaterialLaw(const int num, const double rho, const double E, const double nu):
dG3DMaterialLaw(num,rho,false){
_elasticLaw = new smallStrainLinearElasticPotential(num,E,nu);
}
dG3DLinearElasticMaterialLaw::dG3DLinearElasticMaterialLaw(const dG3DLinearElasticMaterialLaw& src):dG3DMaterialLaw(src){
_elasticLaw = NULL;
if (src._elasticLaw != NULL){
_elasticLaw = src._elasticLaw->clone();
}
};
void dG3DLinearElasticMaterialLaw::setElasticPotential(const elasticPotential* p) {
if (_elasticLaw){
delete _elasticLaw;
}
_elasticLaw = p->clone();
}
void dG3DLinearElasticMaterialLaw::initLaws(const std::map<int,materialLaw*> &maplaw){
if (!_initialized){
elasticStiffness = _elasticLaw->getElasticTensor();
_initialized = true;
};
};
void dG3DLinearElasticMaterialLaw::createIPState(IPStateBase* &ips,const bool* state_,const MElement *ele, const int nbFF_, const IntPt *GP, const int gpt) const
{
// check interface element
bool inter=true;
const MInterfaceElement *iele = dynamic_cast<const MInterfaceElement*>(ele);
if(iele==NULL) inter=false;
IPVariable* ipvi = new LinearElasticDG3DIPVariable(inter);
IPVariable* ipv1 = new LinearElasticDG3DIPVariable(inter);
IPVariable* ipv2 = new LinearElasticDG3DIPVariable(inter);
if(ips != NULL) delete ips;
ips = new IP3State(state_,ipvi,ipv1,ipv2);
}
void dG3DLinearElasticMaterialLaw::createIPVariable(IPVariable* &ipv,const MElement *ele, const int nbFF_, const IntPt *GP, const int gpt) const
{
if(ipv !=NULL) delete ipv;
bool inter=true;
const MInterfaceElement *iele = dynamic_cast<const MInterfaceElement*>(ele);
if(iele == NULL){
inter=false;
}
ipv = new LinearElasticDG3DIPVariable(inter);
}
void dG3DLinearElasticMaterialLaw::stress(IPVariable* ipv, const IPVariable* ipvp, const bool stiff, const bool checkfrac)
{
/* get ipvariable */
dG3DIPVariableBase* ipvcur = static_cast<dG3DIPVariableBase*>(ipv);;
const dG3DIPVariableBase* ipvprev = static_cast<const dG3DIPVariableBase*>(ipvp);;
/* strain xyz */
const STensor3& F = ipvcur->getConstRefToDeformationGradient();
STensor3& P = ipvcur->getRefToFirstPiolaKirchhoffStress();
STensor43& H = ipvcur->getRefToTangentModuli();
_elasticLaw->constitutive(F,P,stiff,&H);
ipvcur->setRefToDGElasticTangentModuli(this->elasticStiffness);
}
double dG3DLinearElasticMaterialLaw::scaleFactor() const{ return _elasticLaw->getYoungModulus();
};
double dG3DLinearElasticMaterialLaw::soundSpeed() const
{
double E = _elasticLaw->getYoungModulus();
double nu = _elasticLaw->getPoissonRatio();
double factornu = (1.-nu)/((1.+nu)*(1.-2.*nu));
return sqrt(E*factornu/_rho);
};
J2SmallStrainDG3DMaterialLaw::J2SmallStrainDG3DMaterialLaw(const int num,const double rho,
double E,const double nu,const double sy0,const double h,
const double tol,const bool tangentByPert, const double pert) : dG3DMaterialLaw(num,rho,true),
_j2law(num,E,nu,rho,sy0,h,tol,tangentByPert,pert)
{
fillElasticStiffness(E, nu, elasticStiffness);
}
J2SmallStrainDG3DMaterialLaw::J2SmallStrainDG3DMaterialLaw(const int num,const double rho,
double E, const double nu, const J2IsotropicHardening &_j2IH,
const double tol, const bool tangentByPert, const double pert) : dG3DMaterialLaw(num,rho,true),
_j2law(num,E,nu,rho,_j2IH,tol,tangentByPert,pert)
{
fillElasticStiffness(E, nu, elasticStiffness);
}
J2SmallStrainDG3DMaterialLaw::J2SmallStrainDG3DMaterialLaw(const J2SmallStrainDG3DMaterialLaw &source) : dG3DMaterialLaw(source), _j2law(source._j2law)
{
}
void J2SmallStrainDG3DMaterialLaw::createIPState(IPStateBase* &ips,const bool* state_,const MElement *ele, const int nbFF_, const IntPt *GP, const int gpt) const
{
// check interface element
bool inter=true;
const MInterfaceElement *iele = dynamic_cast<const MInterfaceElement*>(ele);
if(iele==NULL) inter=false;
IPVariable* ipvi = new J2SmallStrainDG3DIPVariable(_j2law,inter);
IPVariable* ipv1 = new J2SmallStrainDG3DIPVariable(_j2law,inter);
IPVariable* ipv2 = new J2SmallStrainDG3DIPVariable(_j2law,inter);
if(ips != NULL) delete ips;
ips = new IP3State(state_,ipvi,ipv1,ipv2);
}
void J2SmallStrainDG3DMaterialLaw::createIPVariable(IPVariable* &ipv,const MElement *ele, const int nbFF_, const IntPt *GP, const int gpt) const
{
if(ipv !=NULL) delete ipv;
bool inter=true;
const MInterfaceElement *iele = dynamic_cast<const MInterfaceElement*>(ele);
if(iele == NULL){
inter=false;
}
ipv = new J2SmallStrainDG3DIPVariable(_j2law,inter);
}
void J2SmallStrainDG3DMaterialLaw::stress(IPVariable* ipv, const IPVariable* ipvp, const bool stiff, const bool checkfrac)
{
/* get ipvariable */
J2SmallStrainDG3DIPVariable* ipvcur;
const J2SmallStrainDG3DIPVariable* ipvprev;
FractureCohesive3DIPVariable* ipvtmp = dynamic_cast<FractureCohesive3DIPVariable*>(ipv);
if(ipvtmp !=NULL)
{
ipvcur = static_cast<J2SmallStrainDG3DIPVariable*>(ipvtmp->getIPvBulk());
const FractureCohesive3DIPVariable *ipvtmp2 = static_cast<const FractureCohesive3DIPVariable*>(ipvp);
ipvprev = static_cast<const J2SmallStrainDG3DIPVariable*>(ipvtmp2->getIPvBulk()); }
else
{
ipvcur = static_cast<J2SmallStrainDG3DIPVariable*>(ipv);
ipvprev = static_cast<const J2SmallStrainDG3DIPVariable*>(ipvp);
}
/* strain xyz */
const STensor3& F1 = ipvcur->getConstRefToDeformationGradient();
const STensor3& F0 = ipvprev->getConstRefToDeformationGradient();
/* data for J2 law */
IPJ2linear* q1 = ipvcur->getIPJ2linear();
const IPJ2linear* q0 = ipvprev->getIPJ2linear();
/* compute stress */
_j2law.constitutive(F0,F1,ipvcur->getRefToFirstPiolaKirchhoffStress(),q0,q1,ipvcur->getRefToTangentModuli(),stiff);
ipvcur->setRefToDGElasticTangentModuli(this->elasticStiffness);
}
J2LinearDG3DMaterialLaw::J2LinearDG3DMaterialLaw(const int num,const double rho,
double E,const double nu,const double sy0,const double h,
const double tol, const bool matrixBypert, const double tolPert) : dG3DMaterialLaw(num,rho,true),
_j2law(num,E,nu,rho,sy0,h,tol,matrixBypert,tolPert)
{
fillElasticStiffness(E, nu, elasticStiffness);
}
void J2LinearDG3DMaterialLaw::setStrainOrder(const int order){
_j2law.setStrainOrder(order);
};
J2LinearDG3DMaterialLaw::J2LinearDG3DMaterialLaw(const int num,const double rho,
double E, const double nu, const J2IsotropicHardening &_j2IH,
const double tol,const bool matrixBypert, const double tolPert) :
dG3DMaterialLaw(num,rho,true), _j2law(num,E,nu,rho,_j2IH,tol,matrixBypert,tolPert)
{
fillElasticStiffness(E, nu, elasticStiffness);
}
J2LinearDG3DMaterialLaw::J2LinearDG3DMaterialLaw(const J2LinearDG3DMaterialLaw &source) : dG3DMaterialLaw(source), _j2law(source._j2law)
{
}
void J2LinearDG3DMaterialLaw::createIPState(IPStateBase* &ips,const bool* state_,const MElement *ele, const int nbFF_, const IntPt *GP, const int gpt) const
{
// check interface element
bool inter=true;
const MInterfaceElement *iele = dynamic_cast<const MInterfaceElement*>(ele);
if(iele==NULL) inter=false;
IPVariable* ipvi = new J2LinearDG3DIPVariable(_j2law,inter);
IPVariable* ipv1 = new J2LinearDG3DIPVariable(_j2law,inter);
IPVariable* ipv2 = new J2LinearDG3DIPVariable(_j2law,inter);
if(ips != NULL) delete ips;
ips = new IP3State(state_,ipvi,ipv1,ipv2);
}
void J2LinearDG3DMaterialLaw::createIPVariable(IPVariable* &ipv,const MElement *ele, const int nbFF_, const IntPt *GP, const int gpt) const
{
if(ipv !=NULL) delete ipv;
bool inter=true;
const MInterfaceElement *iele = dynamic_cast<const MInterfaceElement*>(ele);
if(iele == NULL){
inter=false;
}
ipv = new J2LinearDG3DIPVariable(_j2law,inter);
}
void J2LinearDG3DMaterialLaw::stress(IPVariable* ipv, const IPVariable* ipvp, const bool stiff, const bool checkfrac)
{
/* get ipvariable */
J2LinearDG3DIPVariable* ipvcur;
const J2LinearDG3DIPVariable* ipvprev;
FractureCohesive3DIPVariable* ipvtmp = dynamic_cast<FractureCohesive3DIPVariable*>(ipv);
if(ipvtmp !=NULL)
{
ipvcur = static_cast<J2LinearDG3DIPVariable*>(ipvtmp->getIPvBulk());
const FractureCohesive3DIPVariable *ipvtmp2 = static_cast<const FractureCohesive3DIPVariable*>(ipvp);
ipvprev = static_cast<const J2LinearDG3DIPVariable*>(ipvtmp2->getIPvBulk());
}
else
{
ipvcur = static_cast<J2LinearDG3DIPVariable*>(ipv);
ipvprev = static_cast<const J2LinearDG3DIPVariable*>(ipvp);
}
/* strain xyz */
const STensor3& F1 = ipvcur->getConstRefToDeformationGradient();
const STensor3& F0 = ipvprev->getConstRefToDeformationGradient();
/* data for J2 law */
IPJ2linear* q1 = ipvcur->getIPJ2linear();
const IPJ2linear* q0 = ipvprev->getIPJ2linear();
STensor43& elasticL = ipvcur->getRefToElasticTangentModuli();
/* compute stress */
_j2law.constitutive(F0,F1,ipvcur->getRefToFirstPiolaKirchhoffStress(),q0,q1,ipvcur->getRefToTangentModuli(),stiff,&elasticL);
ipvcur->setRefToDGElasticTangentModuli(this->elasticStiffness);
}
HyperViscoElasticDG3DMaterialLaw::HyperViscoElasticDG3DMaterialLaw(const int num, const double rho, const double E,const double nu,
const bool matrixbyPerturbation, const double pert): dG3DMaterialLaw(num,rho), _viscoLaw(num,E,nu,rho,matrixbyPerturbation,pert){
fillElasticStiffness(E, nu, elasticStiffness);
};
HyperViscoElasticDG3DMaterialLaw::HyperViscoElasticDG3DMaterialLaw(const HyperViscoElasticDG3DMaterialLaw& src): dG3DMaterialLaw(src),
_viscoLaw(src._viscoLaw){};
void HyperViscoElasticDG3DMaterialLaw::setStrainOrder(const int order){
_viscoLaw.setStrainOrder(order);
};
void HyperViscoElasticDG3DMaterialLaw::setViscoelasticMethod(const int method){
_viscoLaw.setViscoelasticMethod(method);
};
void HyperViscoElasticDG3DMaterialLaw::setViscoElasticNumberOfElement(const int N){
_viscoLaw.setViscoElasticNumberOfElement(N);
};
void HyperViscoElasticDG3DMaterialLaw::setViscoElasticData(const int i, const double Ei, const double taui){
_viscoLaw.setViscoElasticData(i,Ei,taui);
};
void HyperViscoElasticDG3DMaterialLaw::setViscoElasticData_Bulk(const int i, const double Ki, const double ki){
_viscoLaw.setViscoElasticData_Bulk(i,Ki,ki);
};
void HyperViscoElasticDG3DMaterialLaw::setViscoElasticData_Shear(const int i, const double Gi, const double gi){
_viscoLaw.setViscoElasticData_Shear(i,Gi,gi);
};
void HyperViscoElasticDG3DMaterialLaw::setViscoElasticData(const std::string filename){
_viscoLaw.setViscoElasticData(filename);
};
void HyperViscoElasticDG3DMaterialLaw::createIPState(IPStateBase* &ips,const bool* state_,const MElement *ele, const int nbFF_, const IntPt *GP, const int gpt) const
{
// check interface element
bool inter=true;
const MInterfaceElement *iele = dynamic_cast<const MInterfaceElement*>(ele); if(iele==NULL) inter=false;
IPVariable* ipvi = new HyperViscoElasticdG3DIPVariable(_viscoLaw,inter);
IPVariable* ipv1 = new HyperViscoElasticdG3DIPVariable(_viscoLaw,inter);
IPVariable* ipv2 = new HyperViscoElasticdG3DIPVariable(_viscoLaw,inter);
if(ips != NULL) delete ips;
ips = new IP3State(state_,ipvi,ipv1,ipv2);
}
void HyperViscoElasticDG3DMaterialLaw::createIPVariable(IPVariable* &ipv,const MElement *ele, const int nbFF_, const IntPt *GP, const int gpt) const
{
if(ipv !=NULL) delete ipv;
bool inter=true;
const MInterfaceElement *iele = dynamic_cast<const MInterfaceElement*>(ele);
if(iele == NULL){
inter=false;
}
ipv = new HyperViscoElasticdG3DIPVariable(_viscoLaw,inter);
}
void HyperViscoElasticDG3DMaterialLaw::stress(IPVariable* ipv, const IPVariable* ipvp, const bool stiff, const bool checkfrac)
{
/* get ipvariable */
HyperViscoElasticdG3DIPVariable* ipvcur;
const HyperViscoElasticdG3DIPVariable* ipvprev;
FractureCohesive3DIPVariable* ipvtmp = dynamic_cast<FractureCohesive3DIPVariable*>(ipv);
if(ipvtmp !=NULL)
{
ipvcur = static_cast<HyperViscoElasticdG3DIPVariable*>(ipvtmp->getIPvBulk());
const FractureCohesive3DIPVariable *ipvtmp2 = static_cast<const FractureCohesive3DIPVariable*>(ipvp);
ipvprev = static_cast<const HyperViscoElasticdG3DIPVariable*>(ipvtmp2->getIPvBulk());
}
else
{
ipvcur = static_cast<HyperViscoElasticdG3DIPVariable*>(ipv);
ipvprev = static_cast<const HyperViscoElasticdG3DIPVariable*>(ipvp);
}
/* strain xyz */
const STensor3& F1 = ipvcur->getConstRefToDeformationGradient();
const STensor3& F0 = ipvprev->getConstRefToDeformationGradient();
/* data for J2 law */
IPHyperViscoElastic* q1 = ipvcur->getIPHyperViscoElastic();
const IPHyperViscoElastic* q0 = ipvprev->getIPHyperViscoElastic();
/* compute stress */
_viscoLaw.constitutive(F0,F1,ipvcur->getRefToFirstPiolaKirchhoffStress(),q0,q1,ipvcur->getRefToTangentModuli(),stiff);
ipvcur->setRefToDGElasticTangentModuli(this->elasticStiffness);
}
HyperViscoElastoPlasticPowerYieldDG3DMaterialLaw::HyperViscoElastoPlasticPowerYieldDG3DMaterialLaw(const int num, const double rho, const double E,const double nu, const double tol ,
const bool matrixbyPerturbation, const double pert):
dG3DMaterialLaw(num,rho), _viscoLaw(num,E,nu,rho,tol,matrixbyPerturbation,pert){};
HyperViscoElastoPlasticPowerYieldDG3DMaterialLaw::HyperViscoElastoPlasticPowerYieldDG3DMaterialLaw(const HyperViscoElastoPlasticPowerYieldDG3DMaterialLaw& src): dG3DMaterialLaw(src),
_viscoLaw(src._viscoLaw){};
void HyperViscoElastoPlasticPowerYieldDG3DMaterialLaw::setStrainOrder(const int order){
_viscoLaw.setStrainOrder(order);
};
void HyperViscoElastoPlasticPowerYieldDG3DMaterialLaw::setViscoelasticMethod(const int method){
_viscoLaw.setViscoelasticMethod(method);
};
void HyperViscoElastoPlasticPowerYieldDG3DMaterialLaw::setViscoElasticNumberOfElement(const int N){
_viscoLaw.setViscoElasticNumberOfElement(N);
};
void HyperViscoElastoPlasticPowerYieldDG3DMaterialLaw::setViscoElasticData(const int i, const double Ei, const double taui){
_viscoLaw.setViscoElasticData(i,Ei,taui);};
void HyperViscoElastoPlasticPowerYieldDG3DMaterialLaw::setViscoElasticData_Bulk(const int i, const double Ki, const double ki){
_viscoLaw.setViscoElasticData_Bulk(i,Ki,ki);
};
void HyperViscoElastoPlasticPowerYieldDG3DMaterialLaw::setViscoElasticData_Shear(const int i, const double Gi, const double gi){
_viscoLaw.setViscoElasticData_Shear(i,Gi,gi);
};
void HyperViscoElastoPlasticPowerYieldDG3DMaterialLaw::setViscoElasticData(const std::string filename){
_viscoLaw.setViscoElasticData(filename);
};
void HyperViscoElastoPlasticPowerYieldDG3DMaterialLaw::setCompressionHardening(const J2IsotropicHardening& comp){
_viscoLaw.setCompressionHardening(comp);
};
void HyperViscoElastoPlasticPowerYieldDG3DMaterialLaw::setTractionHardening(const J2IsotropicHardening& trac){
_viscoLaw.setTractionHardening(trac);
};
void HyperViscoElastoPlasticPowerYieldDG3DMaterialLaw::setKinematicHardening(const kinematicHardening& kin){
_viscoLaw.setKinematicHardening(kin);
};
void HyperViscoElastoPlasticPowerYieldDG3DMaterialLaw::setViscosityEffect(const viscosityLaw& v, const double p){
_viscoLaw.setViscosityEffect(v,p);
};
void HyperViscoElastoPlasticPowerYieldDG3DMaterialLaw::setYieldPowerFactor(const double n){
_viscoLaw.setPowerFactor(n);
};
void HyperViscoElastoPlasticPowerYieldDG3DMaterialLaw::nonAssociatedFlowRuleFactor(const double b){
_viscoLaw.nonAssociatedFlowRuleFactor(b);
};
void HyperViscoElastoPlasticPowerYieldDG3DMaterialLaw::setPlasticPoissonRatio(const double nup){
_viscoLaw.setPlasticPoissonRatio(nup);
};
void HyperViscoElastoPlasticPowerYieldDG3DMaterialLaw::setNonAssociatedFlow(const bool flag){
_viscoLaw.setNonAssociatedFlow(flag);
};
void HyperViscoElastoPlasticPowerYieldDG3DMaterialLaw::createIPState(IPStateBase* &ips,const bool* state_,const MElement *ele, const int nbFF_, const IntPt *GP, const int gpt) const
{
// check interface element
bool inter=true;
const MInterfaceElement *iele = dynamic_cast<const MInterfaceElement*>(ele);
if(iele==NULL) inter=false;
IPVariable* ipvi = new HyperViscoElastoPlasticdG3DIPVariable(_viscoLaw,inter);
IPVariable* ipv1 = new HyperViscoElastoPlasticdG3DIPVariable(_viscoLaw,inter);
IPVariable* ipv2 = new HyperViscoElastoPlasticdG3DIPVariable(_viscoLaw,inter);
if(ips != NULL) delete ips;
ips = new IP3State(state_,ipvi,ipv1,ipv2);
}
void HyperViscoElastoPlasticPowerYieldDG3DMaterialLaw::createIPVariable(IPVariable* &ipv,const MElement *ele, const int nbFF_, const IntPt *GP, const int gpt) const
{
if(ipv !=NULL) delete ipv;
bool inter=true;
const MInterfaceElement *iele = dynamic_cast<const MInterfaceElement*>(ele);
if(iele == NULL){
inter=false;
}
ipv = new HyperViscoElastoPlasticdG3DIPVariable(_viscoLaw,inter);
}
void HyperViscoElastoPlasticPowerYieldDG3DMaterialLaw::stress(IPVariable* ipv, const IPVariable* ipvp, const bool stiff, const bool checkfrac)
{
/* get ipvariable */
HyperViscoElastoPlasticdG3DIPVariable* ipvcur;
const HyperViscoElastoPlasticdG3DIPVariable* ipvprev;
FractureCohesive3DIPVariable* ipvtmp = dynamic_cast<FractureCohesive3DIPVariable*>(ipv);
if(ipvtmp !=NULL)
{
ipvcur = static_cast<HyperViscoElastoPlasticdG3DIPVariable*>(ipvtmp->getIPvBulk());
const FractureCohesive3DIPVariable *ipvtmp2 = static_cast<const FractureCohesive3DIPVariable*>(ipvp);
ipvprev = static_cast<const HyperViscoElastoPlasticdG3DIPVariable*>(ipvtmp2->getIPvBulk());
}
else
{
ipvcur = static_cast<HyperViscoElastoPlasticdG3DIPVariable*>(ipv);
ipvprev = static_cast<const HyperViscoElastoPlasticdG3DIPVariable*>(ipvp);
}
/* strain xyz */
const STensor3& F1 = ipvcur->getConstRefToDeformationGradient();
const STensor3& F0 = ipvprev->getConstRefToDeformationGradient();
/* data for J2 law */
IPHyperViscoElastoPlastic* q1 = ipvcur->getIPHyperViscoElastoPlastic();
const IPHyperViscoElastoPlastic* q0 = ipvprev->getIPHyperViscoElastoPlastic();
/* compute stress */
_viscoLaw.constitutive(F0,F1,ipvcur->getRefToFirstPiolaKirchhoffStress(),q0,q1,ipvcur->getRefToTangentModuli(),stiff);
ipvcur->setRefToDGElasticTangentModuli(this->elasticStiffness);
}
//=========================================================ViscoelastiDG3DMaterialLaw===============================================================//
std::vector<double> ViscoelasticDG3DMaterialLaw::read_file(const char* file_name, int file_size)
{
std::ifstream file(file_name);
std::vector<double> output;
if(file.is_open())
{
for(int i = 0; i < file_size; ++i)
{
std::string curr_string;
file >> curr_string;
double curr_number = std::atof(curr_string.c_str());
output.push_back(curr_number);
}
}
return output;
}
ViscoelasticDG3DMaterialLaw::ViscoelasticDG3DMaterialLaw(const int num,const double rho,
const double K, const std::string &array_eta, int array_size_eta, const std::string &array_mu,
const double tol, const bool pert, const double eps) :
dG3DMaterialLaw(num,rho,true), _Vislaw(num,K,rho, array_eta, array_size_eta, array_mu, tol,pert,eps)
{
std::vector<double> vec_mu;
vec_mu = read_file(array_mu.c_str(),array_size_eta+1);
double sum_mu = 0.;
for(int i=0;i<vec_mu.size();++i)
{
sum_mu += vec_mu.at(i);
}
double E = 9.*K*sum_mu/(3.*K+sum_mu);
double nu = (3.*K - 2.*sum_mu)/2./(3.*K + sum_mu);
fillElasticStiffness(E, nu, elasticStiffness);
}
ViscoelasticDG3DMaterialLaw::ViscoelasticDG3DMaterialLaw(const ViscoelasticDG3DMaterialLaw &source) : dG3DMaterialLaw(source), _Vislaw(source._Vislaw){
}
void ViscoelasticDG3DMaterialLaw::createIPState(IPStateBase* &ips,const bool* state_,const MElement *ele, const int nbFF_, const IntPt *GP, const int gpt) const
{
bool inter=true;
const MInterfaceElement *iele = dynamic_cast<const MInterfaceElement*>(ele);
if(iele==NULL) inter=false;
IPVariable* ipvi = new ViscoelasticDG3DIPVariable(_Vislaw,inter);
IPVariable* ipv1 = new ViscoelasticDG3DIPVariable(_Vislaw,inter);
IPVariable* ipv2 = new ViscoelasticDG3DIPVariable(_Vislaw,inter);
if(ips != NULL) delete ips;
ips = new IP3State(state_,ipvi,ipv1,ipv2);
}
void ViscoelasticDG3DMaterialLaw::createIPVariable(IPVariable* &ipv,const MElement *ele, const int nbFF_, const IntPt *GP, const int gpt) const
{
if(ipv !=NULL) delete ipv;
bool inter=true;
const MInterfaceElement *iele = dynamic_cast<const MInterfaceElement*>(ele);
if(iele == NULL){
inter=false;
}
ipv = new ViscoelasticDG3DIPVariable(_Vislaw,inter);
}
void ViscoelasticDG3DMaterialLaw::stress(IPVariable* ipv, const IPVariable* ipvp, const bool stiff, const bool checkfrac)
{
// get ipvariable
ViscoelasticDG3DIPVariable* ipvcur;
const ViscoelasticDG3DIPVariable* ipvprev;
ipvcur = static_cast<ViscoelasticDG3DIPVariable*>(ipv);
ipvprev = static_cast<const ViscoelasticDG3DIPVariable*>(ipvp);
const STensor3& F1 = ipvcur->getConstRefToDeformationGradient();
const STensor3& F0 = ipvprev->getConstRefToDeformationGradient();
IPViscoelastic* IPVis = ipvcur->getIPViscoelastic();
const IPViscoelastic* IPVisprev = ipvprev->getIPViscoelastic();
// compute stress
_Vislaw.constitutive(F0,F1,ipvcur->getRefToFirstPiolaKirchhoffStress(),IPVis,IPVisprev,stiff);
ipvcur->setRefToDGElasticTangentModuli(this->elasticStiffness);
}
//-------------------------------------------EOSDG3DMaterialLaw (compatible for all the other material laws)-------------------//
EOSDG3DMaterialLaw::EOSDG3DMaterialLaw(const int num,
const double rho,const double K, materialLaw *_mlaw, const double Pressure0) : dG3DMaterialLaw(num,rho,true),_EOSlaw(num,rho,K,Pressure0)
{
// Accoustic EOS for non-artificial viscosity case
maplaw.insert(std::pair<int,materialLaw*>(_mlaw->getNum(),_mlaw));
}
EOSDG3DMaterialLaw::EOSDG3DMaterialLaw(int flag, const int num,
const double rho,const double P1, double P2, double P3, materialLaw *_mlaw, const double Pressure0) : dG3DMaterialLaw(num,rho,true),_EOSlaw(flag,num,rho,P1,P2,P3,Pressure0)
{
// P1 = s, P2 = C0, P3 = Gamma: Hugoniot EOS for non-artificial viscosity case
// P1 = K, P2 = c1, P3 = cl: Accoustic EOS for artificial viscosity case
maplaw.insert(std::pair<int,materialLaw*>(_mlaw->getNum(),_mlaw));
}
EOSDG3DMaterialLaw::EOSDG3DMaterialLaw(int flag, const int num,
const double rho,const double s, const double C0, double c1, double cl, const double Gamma, materialLaw *_mlaw) : dG3DMaterialLaw(num,rho,true),_EOSlaw(flag, num,rho, s, C0,c1,cl,Gamma)
{
// Hugoniot EOS for artificial viscosity case
maplaw.insert(std::pair<int,materialLaw*>(_mlaw->getNum(),_mlaw));
}
EOSDG3DMaterialLaw::EOSDG3DMaterialLaw(int flag, int flagC, const int num,
const double rho,const double K, double c1, double cl, double ESqrd, double ESlin, materialLaw *_mlaw, const double Pressure0) : dG3DMaterialLaw(num,rho,true),_EOSlaw(flag,num,rho,K,c1,cl,Pressure0)
{
maplaw.insert(std::pair<int,materialLaw*>(_mlaw->getNum(),_mlaw));
}
EOSDG3DMaterialLaw::EOSDG3DMaterialLaw(const EOSDG3DMaterialLaw &source) : dG3DMaterialLaw(source),_EOSlaw(source._EOSlaw),maplaw(source.maplaw)
{
}
void EOSDG3DMaterialLaw::createIPState(IPStateBase* &ips,const bool* state_,const MElement *ele, const int nbFF_, const IntPt *GP, const int gpt) const
{
dG3DMaterialLaw *DG3Dit = dynamic_cast<dG3DMaterialLaw*>(maplaw.begin()->second);
materialLaw::matname mlaw_dev = DG3Dit->getType();
if(mlaw_dev == materialLaw::j2linear )//J2linear
{
J2LinearDG3DMaterialLaw *DG3DJ2linearit = dynamic_cast<J2LinearDG3DMaterialLaw*>(DG3Dit);
bool inter=true;
const MInterfaceElement *iele = dynamic_cast<const MInterfaceElement*>(ele);
if(iele==NULL) inter=false;
double size = 2.*(( const_cast<MElement*>(ele) )->getInnerRadius());
J2LinearDG3DIPVariable* ipvi_J2 = new J2LinearDG3DIPVariable(DG3DJ2linearit->_j2law,inter);
J2LinearDG3DIPVariable* ipv1_J2 = new J2LinearDG3DIPVariable(DG3DJ2linearit->_j2law,inter);
J2LinearDG3DIPVariable* ipv2_J2 = new J2LinearDG3DIPVariable(DG3DJ2linearit->_j2law,inter);
IPVariable* ipvi = new EOSDG3DIPVariable(size, ipvi_J2->getIPJ2linear(),inter);
IPVariable* ipv1 = new EOSDG3DIPVariable(size, ipv1_J2->getIPJ2linear(),inter);
IPVariable* ipv2 = new EOSDG3DIPVariable(size, ipv2_J2->getIPJ2linear(),inter);
if(ips != NULL) delete ips;
ips = new IP3State(state_,ipvi,ipv1,ipv2);
}
else if(mlaw_dev == materialLaw::viscoelastic)//viscoelastic
{
ViscoelasticDG3DMaterialLaw *DG3Dvisit = dynamic_cast<ViscoelasticDG3DMaterialLaw*>(DG3Dit);
bool inter=true;
const MInterfaceElement *iele = dynamic_cast<const MInterfaceElement*>(ele);
if(iele==NULL) inter=false;
double size = 2.*(( const_cast<MElement*>(ele) )->getInnerRadius());
ViscoelasticDG3DIPVariable* ipvi_vis = new ViscoelasticDG3DIPVariable(DG3Dvisit->_Vislaw,inter);
ViscoelasticDG3DIPVariable* ipv1_vis = new ViscoelasticDG3DIPVariable(DG3Dvisit->_Vislaw,inter);
ViscoelasticDG3DIPVariable* ipv2_vis = new ViscoelasticDG3DIPVariable(DG3Dvisit->_Vislaw,inter);
IPVariable* ipvi = new EOSDG3DIPVariable(size, ipvi_vis->getIPViscoelastic(),inter);
IPVariable* ipv1 = new EOSDG3DIPVariable(size, ipv1_vis->getIPViscoelastic(),inter);
IPVariable* ipv2 = new EOSDG3DIPVariable(size, ipv2_vis->getIPViscoelastic(),inter);
if(ips != NULL) delete ips;
ips = new IP3State(state_,ipvi,ipv1,ipv2);
}
}
void EOSDG3DMaterialLaw::createIPVariable(IPVariable* &ipv,const MElement *ele, const int nbFF_, const IntPt *GP, const int gpt) const
{
dG3DMaterialLaw *DG3Dit = dynamic_cast<dG3DMaterialLaw*>(maplaw.begin()->second);
matname mlaw_dev = DG3Dit->getType();
if(mlaw_dev == materialLaw::j2linear)//J2linear
{
J2LinearDG3DMaterialLaw *DG3DJ2linearit = dynamic_cast<J2LinearDG3DMaterialLaw*>(DG3Dit);
bool inter=true;
const MInterfaceElement *iele = dynamic_cast<const MInterfaceElement*>(ele);
if(iele==NULL) inter=false;
double size = 2.*(( const_cast<MElement*>(ele) )->getInnerRadius());
J2LinearDG3DIPVariable* ipv_J2 = new J2LinearDG3DIPVariable(DG3DJ2linearit->_j2law,inter);
ipv = new EOSDG3DIPVariable(size, ipv_J2,inter);
}
else if(mlaw_dev == materialLaw::viscoelastic)//viscoelastic
{
ViscoelasticDG3DMaterialLaw *DG3Dvisit = dynamic_cast<ViscoelasticDG3DMaterialLaw*>(DG3Dit);
bool inter=true;
const MInterfaceElement *iele = dynamic_cast<const MInterfaceElement*>(ele);
if(iele==NULL) inter=false;
double size = 2.*(( const_cast<MElement*>(ele) )->getInnerRadius());
ViscoelasticDG3DIPVariable* ipv_vis = new ViscoelasticDG3DIPVariable(DG3Dvisit->_Vislaw,inter);
ipv = new EOSDG3DIPVariable(size, ipv_vis,inter);
}
}
void EOSDG3DMaterialLaw::stress(IPVariable* ipv, const IPVariable* ipvp, const bool stiff, const bool checkfrac)
{
dG3DMaterialLaw *DG3Dit = dynamic_cast<dG3DMaterialLaw*>(maplaw.begin()->second);
materialLaw::matname mlaw_dev = DG3Dit->getType();
EOSDG3DIPVariable *DG3DEOSipv = dynamic_cast<EOSDG3DIPVariable*>(ipv);
const EOSDG3DIPVariable *DG3DEOSipvp = static_cast<const EOSDG3DIPVariable*>(ipvp);
IPEOS* IPeosCUR = DG3DEOSipv->getIPEOS();
const IPEOS* IPeosPREV = DG3DEOSipvp->getIPEOS();
STensor3 Fn(0.);
Fn = DG3DEOSipv->getConstRefToDeformationGradient();
STensor3 F0(0.);
F0 = DG3DEOSipvp->getConstRefToDeformationGradient();
STensor3 Fdot(0.);
// STensor43 ViselasticStiff(0.);
if(mlaw_dev == materialLaw::j2linear) //J2linear
{
J2LinearDG3DMaterialLaw *DG3DJ2linearit = dynamic_cast<J2LinearDG3DMaterialLaw*>(DG3Dit);
IPJ2linear *J2linearipv = dynamic_cast<IPJ2linear*>(IPeosCUR->_IPV);
IPJ2linear *J2linearipvp = dynamic_cast<IPJ2linear*>(IPeosPREV->_IPV);
/* compute stress */
DG3DJ2linearit->_j2law.constitutive(F0,Fn,DG3DEOSipv->getRefToFirstPiolaKirchhoffStress(),J2linearipv, J2linearipvp,DG3DEOSipv->getRefToTangentModuli(),stiff);
TotalelasticStiffness = DG3DJ2linearit->elasticStiffness;
// Careful, make sure the volumetric properties of your deviatoric law are the same as in your EOS law, as the elastic stiffness is only calcualted through the deviatoric law
}
else if(mlaw_dev == materialLaw::viscoelastic)//viscoelastic
{
ViscoelasticDG3DMaterialLaw *DG3Dvisit = dynamic_cast<ViscoelasticDG3DMaterialLaw*>(DG3Dit);
IPViscoelastic *Visipv = dynamic_cast<IPViscoelastic*>(IPeosCUR->_IPV);
IPViscoelastic *Visipvp = dynamic_cast<IPViscoelastic*>(IPeosPREV->_IPV);
DG3Dvisit->_Vislaw.constitutive(F0,Fn,DG3DEOSipv->getRefToFirstPiolaKirchhoffStress(),Visipv,Visipvp,stiff);
TotalelasticStiffness = DG3Dvisit->elasticStiffness;
// Careful, make sure the volumetric properties of your deviatoric law are the same as in your EOS law, as the elastic stiffness is only calcualted through the deviatoric law
}
_EOSlaw.constitutive(F0,Fn,DG3DEOSipv->getRefToFirstPiolaKirchhoffStress(),DG3DEOSipv->getFirstPiolaViscous(), DG3DEOSipv->getRefToFirstPiolaKirchhoffStressb4AV(), IPeosCUR,IPeosPREV,Fdot);
DG3DEOSipv->setRefToDGElasticTangentModuli(TotalelasticStiffness);
}
VUMATinterfaceDG3DMaterialLaw::VUMATinterfaceDG3DMaterialLaw(const int num, const double rho, const char *propName) :
dG3DMaterialLaw(num,rho,true),
_vumatlaw(num,rho,propName)
{
double nu = 0.3; // NOTE: This is to get an order of magnitude; another way of doing it would be to
// change vumat_soundspeed to vumat_mu_nu which would give mu and nu
// then reconstruct mlawVUMATinterface::soundSpeed to use this instead
// and create mlawVUMATinterface::poissonRation/shearModulus using it too
double soundSpeed = _vumatlaw.soundSpeed();
double lambda_plus_two_mu = soundSpeed*soundSpeed*rho;
double E = lambda_plus_two_mu * ( ((1+nu)*(1-2*nu)) / (1-nu) );
fillElasticStiffness(E, nu, elasticStiffness);
}
VUMATinterfaceDG3DMaterialLaw::VUMATinterfaceDG3DMaterialLaw(const VUMATinterfaceDG3DMaterialLaw &source) :
dG3DMaterialLaw(source), _vumatlaw(source._vumatlaw)
{
}
void VUMATinterfaceDG3DMaterialLaw::createIPState(IPStateBase* &ips,const bool* state_,const MElement *ele, const int nbFF_, const IntPt *GP, const int gpt) const
{
// check interface element
bool inter=true;
const MInterfaceElement *iele = dynamic_cast<const MInterfaceElement*>(ele);
if(iele==NULL) inter=false;
int nsdv = _vumatlaw.getNsdv();
double size = 2.*(( const_cast<MElement*>(ele) )->getInnerRadius());
IPVariable* ipvi = new VUMATinterfaceDG3DIPVariable(nsdv,size,inter);
IPVariable* ipv1 = new VUMATinterfaceDG3DIPVariable(nsdv,size,inter);
IPVariable* ipv2 = new VUMATinterfaceDG3DIPVariable(nsdv,size,inter);
if(ips != NULL) delete ips;
ips = new IP3State(state_,ipvi,ipv1,ipv2);
}
void VUMATinterfaceDG3DMaterialLaw::createIPVariable(IPVariable* &ipv,const MElement *ele, const int nbFF_, const IntPt *GP, const int gpt) const
{
if(ipv !=NULL) delete ipv;
bool inter=true;
const MInterfaceElement *iele = dynamic_cast<const MInterfaceElement*>(ele);
if(iele == NULL){
inter=false;
}
ipv = new VUMATinterfaceDG3DIPVariable(_vumatlaw.getNsdv(),2.*(( const_cast<MElement*>(ele) )->getInnerRadius()),inter);
}
void VUMATinterfaceDG3DMaterialLaw::stress(IPVariable* ipv, const IPVariable* ipvp, const bool stiff, const bool checkfrac)
{
/* get ipvariable */
VUMATinterfaceDG3DIPVariable* ipvcur;
const VUMATinterfaceDG3DIPVariable* ipvprev;
FractureCohesive3DIPVariable* ipvtmp = dynamic_cast<FractureCohesive3DIPVariable*>(ipv);
if(ipvtmp !=NULL)
{
ipvcur = static_cast<VUMATinterfaceDG3DIPVariable*>(ipvtmp->getIPvBulk());
const FractureCohesive3DIPVariable *ipvtmp2 = static_cast<const FractureCohesive3DIPVariable*>(ipvp);
ipvprev = static_cast<const VUMATinterfaceDG3DIPVariable*>(ipvtmp2->getIPvBulk());
}
else {
ipvcur = static_cast<VUMATinterfaceDG3DIPVariable*>(ipv);
ipvprev = static_cast<const VUMATinterfaceDG3DIPVariable*>(ipvp);
}
/* strain xyz */
const STensor3& F1 = ipvcur->getConstRefToDeformationGradient();
const STensor3& F0 = ipvprev->getConstRefToDeformationGradient();
/* data for VUMAT law */
IPVUMATinterface* q1 = ipvcur->getIPVUMATinterface();
const IPVUMATinterface* q0 = ipvprev->getIPVUMATinterface();
/* compute stress */
ipvcur->getRefToFirstPiolaKirchhoffStress()=( const_cast<VUMATinterfaceDG3DIPVariable*>(ipvprev) )->getRefToFirstPiolaKirchhoffStress();
_vumatlaw.constitutive(F0,F1,ipvcur->getRefToFirstPiolaKirchhoffStress(),q0,q1,ipvcur->getRefToTangentModuli(),stiff);
ipvcur->setRefToDGElasticTangentModuli(this->elasticStiffness);
}
LocalDamageHyperelasticDG3DMaterialLaw::LocalDamageHyperelasticDG3DMaterialLaw(const int num, const double rho, const elasticPotential& elP, const DamageLaw &damLaw):
dG3DMaterialLaw(num,rho,true){
_nlLaw = new mlawLocalDamageHyperelastic(num,rho,elP,damLaw);
double E = elP.getYoungModulus();
double nu = elP.getPoissonRatio();
fillElasticStiffness(E, nu, elasticStiffness);
};
LocalDamageHyperelasticDG3DMaterialLaw::LocalDamageHyperelasticDG3DMaterialLaw(const LocalDamageHyperelasticDG3DMaterialLaw& src):
dG3DMaterialLaw(src)
{
_nlLaw = NULL;
if (src._nlLaw != NULL){
_nlLaw = dynamic_cast<mlawLocalDamageHyperelastic*>(src._nlLaw->clone());
}
};
void LocalDamageHyperelasticDG3DMaterialLaw::setTime(const double t,const double dtime)
{
dG3DMaterialLaw::setTime(t,dtime);
_nlLaw->setTime(t,dtime);
}
void LocalDamageHyperelasticDG3DMaterialLaw::createIPState(IPStateBase* &ips,const bool* state_,const MElement *ele, const int nbFF_, const IntPt *GP, const int gpt) const
{
bool inter=true;
const MInterfaceElement *iele = dynamic_cast<const MInterfaceElement*>(ele);
if(iele==NULL) inter=false;
IPVariable* ipvi = new localDamageIsotropicElasticityDG3DIPVariable(*_nlLaw,inter);
IPVariable* ipv1 = new localDamageIsotropicElasticityDG3DIPVariable(*_nlLaw,inter);
IPVariable* ipv2 = new localDamageIsotropicElasticityDG3DIPVariable(*_nlLaw,inter);
if(ips != NULL) delete ips;
ips = new IP3State(state_,ipvi,ipv1,ipv2);
}
void LocalDamageHyperelasticDG3DMaterialLaw::createIPVariable(IPVariable* &ipv,const MElement *ele, const int nbFF_, const IntPt *GP, const int gpt) const
{
bool inter=true;
const MInterfaceElement *iele = dynamic_cast<const MInterfaceElement*>(ele);
if(iele == NULL) inter=false;
if(ipv !=NULL) delete ipv;
ipv = new localDamageIsotropicElasticityDG3DIPVariable(*_nlLaw,inter);
}
void LocalDamageHyperelasticDG3DMaterialLaw::stress(IPVariable* ipv, const IPVariable* ipvp, const bool stiff, const bool checkfrac)
{
/* get ipvariable */ localDamageIsotropicElasticityDG3DIPVariable* ipvcur;
const localDamageIsotropicElasticityDG3DIPVariable* ipvprev;
FractureCohesive3DIPVariable* ipvtmp = dynamic_cast<FractureCohesive3DIPVariable*>(ipv);
if(ipvtmp !=NULL)
{
ipvcur = static_cast<localDamageIsotropicElasticityDG3DIPVariable*>(ipvtmp->getIPvBulk());
const FractureCohesive3DIPVariable *ipvtmp2 = static_cast<const FractureCohesive3DIPVariable*>(ipvp);
ipvprev = static_cast<const localDamageIsotropicElasticityDG3DIPVariable*>(ipvtmp2->getIPvBulk());
}
else
{
ipvcur = static_cast<localDamageIsotropicElasticityDG3DIPVariable*>(ipv);
ipvprev = static_cast<const localDamageIsotropicElasticityDG3DIPVariable*>(ipvp);
}
/* Deformation gradient F */
const STensor3& F1 = ipvcur->getConstRefToDeformationGradient();
const STensor3& F0 = ipvprev->getConstRefToDeformationGradient();
/* data for non local isotropic elastic law */
IPLocalDamageIsotropicElasticity* q1 = ipvcur->getIPLocalDamageIsotropicElasticity();
const IPLocalDamageIsotropicElasticity* q0 = ipvprev->getIPLocalDamageIsotropicElasticity();
/* compute stress */
_nlLaw->constitutive(F0,F1,ipvcur->getRefToFirstPiolaKirchhoffStress(),q0,q1,ipvcur->getRefToTangentModuli(),stiff);
ipvcur->setRefToDGElasticTangentModuli(this->elasticStiffness); // Is it correct ????
}
// *********************
LocalAnisotropicDamageHyperelasticDG3DMaterialLaw::LocalAnisotropicDamageHyperelasticDG3DMaterialLaw(const int num, const double rho,
const elasticPotential& elP,const DamageLaw &damLaw): LocalDamageHyperelasticDG3DMaterialLaw(num,rho,elP,damLaw){
if (_nlLaw != NULL) delete _nlLaw;
_nlLaw = new mlawLocalAnisotropicDamageHyperelastic(num,rho,elP,damLaw);
};
//
LocalDamageJ2HyperDG3DMaterialLaw::LocalDamageJ2HyperDG3DMaterialLaw(const int num,const double rho,
double E,const double nu,const J2IsotropicHardening &_j2IH,
const DamageLaw &_damLaw,const double tol) :
dG3DMaterialLaw(num,rho,true)
{
_nldJ2Hyperlaw = new mlawLocalDamageJ2Hyper(num,E,nu,rho,_j2IH,_damLaw,tol);
fillElasticStiffness(E, nu, elasticStiffness);
}
LocalDamageJ2HyperDG3DMaterialLaw::~LocalDamageJ2HyperDG3DMaterialLaw() {
if (_nldJ2Hyperlaw !=NULL){
delete _nldJ2Hyperlaw;
_nldJ2Hyperlaw = NULL;
}
};
LocalDamageJ2HyperDG3DMaterialLaw::LocalDamageJ2HyperDG3DMaterialLaw(const LocalDamageJ2HyperDG3DMaterialLaw &source) :
dG3DMaterialLaw(source)
{
_nldJ2Hyperlaw = NULL;
if (source._nldJ2Hyperlaw !=NULL){
_nldJ2Hyperlaw = new mlawLocalDamageJ2Hyper(*(source._nldJ2Hyperlaw));
}
}
void
LocalDamageJ2HyperDG3DMaterialLaw::setTime(const double t,const double dtime){
dG3DMaterialLaw::setTime(t,dtime);
_nldJ2Hyperlaw->setTime(t,dtime);
}
materialLaw::matname
LocalDamageJ2HyperDG3DMaterialLaw::getType() const {return _nldJ2Hyperlaw->getType();}
void LocalDamageJ2HyperDG3DMaterialLaw::createIPState(IPStateBase* &ips,const bool* state_,const MElement *ele, const int nbFF_, const IntPt *GP, const int gpt) const
{
// check interface element
bool inter=true;
const MInterfaceElement *iele = dynamic_cast<const MInterfaceElement*>(ele);
if(iele==NULL) inter=false;
IPVariable* ipvi = new localDamageJ2HyperDG3DIPVariable(*_nldJ2Hyperlaw,inter);
IPVariable* ipv1 = new localDamageJ2HyperDG3DIPVariable(*_nldJ2Hyperlaw,inter);
IPVariable* ipv2 = new localDamageJ2HyperDG3DIPVariable(*_nldJ2Hyperlaw,inter);
if(ips != NULL) delete ips;
ips = new IP3State(state_,ipvi,ipv1,ipv2);
}
void LocalDamageJ2HyperDG3DMaterialLaw::createIPVariable(IPVariable* &ipv,const MElement *ele, const int nbFF_, const IntPt *GP, const int gpt) const
{
if(ipv !=NULL) delete ipv;
bool inter=true;
const MInterfaceElement *iele = dynamic_cast<const MInterfaceElement*>(ele);
if(iele == NULL){
inter=false;
}
ipv = new localDamageJ2HyperDG3DIPVariable(*_nldJ2Hyperlaw,inter);
}
double
LocalDamageJ2HyperDG3DMaterialLaw::soundSpeed() const{return _nldJ2Hyperlaw->soundSpeed();} // or change value ??
void LocalDamageJ2HyperDG3DMaterialLaw::stress(IPVariable* ipv, const IPVariable* ipvp, const bool stiff, const bool checkfrac)
{
/* get ipvariable */
localDamageJ2HyperDG3DIPVariable* ipvcur;
const localDamageJ2HyperDG3DIPVariable* ipvprev;
FractureCohesive3DIPVariable* ipvtmp = dynamic_cast<FractureCohesive3DIPVariable*>(ipv);
if(ipvtmp !=NULL)
{
ipvcur = static_cast<localDamageJ2HyperDG3DIPVariable*>(ipvtmp->getIPvBulk());
const FractureCohesive3DIPVariable *ipvtmp2 = static_cast<const FractureCohesive3DIPVariable*>(ipvp);
ipvprev = static_cast<const localDamageJ2HyperDG3DIPVariable*>(ipvtmp2->getIPvBulk());
}
else
{
ipvcur = static_cast<localDamageJ2HyperDG3DIPVariable*>(ipv);
ipvprev = static_cast<const localDamageJ2HyperDG3DIPVariable*>(ipvp);
}
/* strain xyz */
const STensor3& F1 = ipvcur->getConstRefToDeformationGradient();
const STensor3& F0 = ipvprev->getConstRefToDeformationGradient();
/* data for J2 law */
IPLocalDamageJ2Hyper* q1 = ipvcur->getIPLocalDamageJ2Hyper();
const IPLocalDamageJ2Hyper* q0 = ipvprev->getIPLocalDamageJ2Hyper();
/* compute stress */
_nldJ2Hyperlaw->constitutive(F0,F1,ipvcur->getRefToFirstPiolaKirchhoffStress(),q0,q1,ipvcur->getRefToTangentModuli(),stiff);
ipvcur->setRefToDGElasticTangentModuli(this->elasticStiffness);
}
TransverseIsotropicDG3DMaterialLaw::TransverseIsotropicDG3DMaterialLaw(
const int num,const double rho,
double E,const double nu,
const double EA, const double GA, const double nu_minor,
const double Ax, const double Ay, const double Az) : dG3DMaterialLaw(num,rho,true),
_tilaw(num,E,nu,rho,EA,GA,nu_minor, Ax,Ay,Az), _E(E), _nu(nu), _nu_minor(nu_minor), _EA(EA), _GA(GA)
{
fillElasticStiffness(E, nu, elasticStiffness);
}
TransverseIsotropicDG3DMaterialLaw::TransverseIsotropicDG3DMaterialLaw(const TransverseIsotropicDG3DMaterialLaw &source) : dG3DMaterialLaw(source), _tilaw(source._tilaw),
_E(source._E), _nu(source._nu), _EA(source._EA),
_GA(source._GA), _nu_minor(source._nu_minor)
{
}
void TransverseIsotropicDG3DMaterialLaw::createIPState(IPStateBase* &ips,const bool* state_,const MElement *ele, const int nbFF_, const IntPt *GP, const int gpt) const
{
// check interface element
bool inter=true;
const MInterfaceElement *iele = dynamic_cast<const MInterfaceElement*>(ele);
if(iele==NULL) inter=false;
IPVariable* ipvi = new TransverseIsotropicDG3DIPVariable(inter);
IPVariable* ipv1 = new TransverseIsotropicDG3DIPVariable(inter);
IPVariable* ipv2 = new TransverseIsotropicDG3DIPVariable(inter);
if(ips != NULL) delete ips;
ips = new IP3State(state_,ipvi,ipv1,ipv2);
}
void TransverseIsotropicDG3DMaterialLaw::createIPVariable(IPVariable* &ipv,const MElement *ele, const int nbFF_, const IntPt *GP, const int gpt) const
{
if(ipv !=NULL) delete ipv;
bool inter=true;
const MInterfaceElement *iele = dynamic_cast<const MInterfaceElement*>(ele);
if(iele == NULL){
inter=false;
}
ipv = new TransverseIsotropicDG3DIPVariable(inter);
}
void TransverseIsotropicDG3DMaterialLaw::stress(IPVariable* ipv, const IPVariable* ipvp, const bool stiff, const bool checkfrac)
{
/* get ipvariable */
TransverseIsotropicDG3DIPVariable* ipvcur;
const TransverseIsotropicDG3DIPVariable* ipvprev;
FractureCohesive3DIPVariable* ipvtmp = dynamic_cast<FractureCohesive3DIPVariable*>(ipv);
if(ipvtmp !=NULL)
{
ipvcur = static_cast<TransverseIsotropicDG3DIPVariable*>(ipvtmp->getIPvBulk());
const FractureCohesive3DIPVariable *ipvtmp2 = static_cast<const FractureCohesive3DIPVariable*>(ipvp);
ipvprev = static_cast<const TransverseIsotropicDG3DIPVariable*>(ipvtmp2->getIPvBulk());
}
else
{
ipvcur = static_cast<TransverseIsotropicDG3DIPVariable*>(ipv);
ipvprev = static_cast<const TransverseIsotropicDG3DIPVariable*>(ipvp);
}
/* strain xyz */
const STensor3& F1 = ipvcur->getConstRefToDeformationGradient();
const STensor3& F0 = ipvprev->getConstRefToDeformationGradient();
/* data for J2 law */
IPTransverseIsotropic* q1 = ipvcur->getIPTransverseIsotropic();
const IPTransverseIsotropic* q0 = ipvprev->getIPTransverseIsotropic();
STensor43& elasticL = ipvcur->getRefToElasticTangentModuli();
/* compute stress */
_tilaw.constitutive(F0,F1,ipvcur->getRefToFirstPiolaKirchhoffStress(),q0,q1,ipvcur->getRefToTangentModuli(),stiff,&elasticL);
ipvcur->setRefToDGElasticTangentModuli(this->elasticStiffness);
}
TransverseIsoCurvatureDG3DMaterialLaw::TransverseIsoCurvatureDG3DMaterialLaw(
const int num,const double rho,
double E,const double nu,
const double EA, const double GA, const double nu_minor,
const double Centroid_x, const double Centroid_y, const double Centroid_z,
const double PointStart1_x, const double PointStart1_y, const double PointStart1_z,
const double PointStart2_x, const double PointStart2_y, const double PointStart2_z,
const double init_Angle) :
dG3DMaterialLaw(num,rho,true),
_tilaw(num,E,nu,rho,EA,GA,nu_minor, Centroid_x, Centroid_y, Centroid_z, PointStart1_x, PointStart1_y, PointStart1_z, PointStart2_x, PointStart2_y, PointStart2_z,init_Angle), _E(E), _nu(nu), _nu_minor(nu_minor), _EA(EA), _GA(GA)
{
fillElasticStiffness(E, nu, elasticStiffness);
}
TransverseIsoCurvatureDG3DMaterialLaw::TransverseIsoCurvatureDG3DMaterialLaw(const TransverseIsoCurvatureDG3DMaterialLaw &source) : dG3DMaterialLaw(source), _tilaw(source._tilaw),
_E(source._E), _nu(source._nu), _EA(source._EA),
_GA(source._GA), _nu_minor(source._nu_minor)
{
}
void TransverseIsoCurvatureDG3DMaterialLaw::createIPState(IPStateBase* &ips,const bool* state_,const MElement *ele, const int nbFF_, const IntPt *GP, const int gpt) const
{
// check interface element
bool inter=true;
const MInterfaceElement *iele = dynamic_cast<const MInterfaceElement*>(ele);
if(iele==NULL) inter=false;
SPoint3 pt_Global;
ele->pnt(GP[gpt].pt[0],GP[gpt].pt[1],GP[gpt].pt[2],pt_Global);
SVector3 GaussP(pt_Global);
IPVariable* ipvi = new TransverseIsoCurvatureDG3DIPVariable(GaussP, _tilaw.getCentroid(), _tilaw.getPointStart1(), _tilaw.getPointStart2(), _tilaw.getInit_Angle(), inter);
IPVariable* ipv1 = new TransverseIsoCurvatureDG3DIPVariable(GaussP, _tilaw.getCentroid(), _tilaw.getPointStart1(), _tilaw.getPointStart2(), _tilaw.getInit_Angle(), inter);
IPVariable* ipv2 = new TransverseIsoCurvatureDG3DIPVariable(GaussP, _tilaw.getCentroid(), _tilaw.getPointStart1(), _tilaw.getPointStart2(), _tilaw.getInit_Angle(), inter);
if(ips != NULL) delete ips;
ips = new IP3State(state_,ipvi,ipv1,ipv2);
}
void TransverseIsoCurvatureDG3DMaterialLaw::createIPVariable(IPVariable* &ipv,const MElement *ele, const int nbFF_, const IntPt *GP, const int gpt) const
{
if(ipv !=NULL) delete ipv;
bool inter=true;
const MInterfaceElement *iele = dynamic_cast<const MInterfaceElement*>(ele);
if(iele == NULL){
inter=false;
}
SPoint3 pt_Global;
ele->pnt(GP[gpt].pt[0],GP[gpt].pt[1],GP[gpt].pt[2],pt_Global);
SVector3 GaussP(pt_Global);
ipv = new TransverseIsoCurvatureDG3DIPVariable(GaussP, _tilaw.getCentroid(), _tilaw.getPointStart1(), _tilaw.getPointStart2(),_tilaw.getInit_Angle(), inter);
}
void TransverseIsoCurvatureDG3DMaterialLaw::stress(IPVariable* ipv, const IPVariable* ipvp, const bool stiff, const bool checkfrac)
{
/* get ipvariable */
TransverseIsoCurvatureDG3DIPVariable* ipvcur; const TransverseIsoCurvatureDG3DIPVariable* ipvprev;
FractureCohesive3DIPVariable* ipvtmp = dynamic_cast<FractureCohesive3DIPVariable*>(ipv);
if(ipvtmp !=NULL)
{
ipvcur = static_cast<TransverseIsoCurvatureDG3DIPVariable*>(ipvtmp->getIPvBulk());
const FractureCohesive3DIPVariable *ipvtmp2 = static_cast<const FractureCohesive3DIPVariable*>(ipvp);
ipvprev = static_cast<const TransverseIsoCurvatureDG3DIPVariable*>(ipvtmp2->getIPvBulk());
}
else
{
ipvcur = static_cast<TransverseIsoCurvatureDG3DIPVariable*>(ipv);
ipvprev = static_cast<const TransverseIsoCurvatureDG3DIPVariable*>(ipvp);
}
/* strain xyz */
const STensor3& F1 = ipvcur->getConstRefToDeformationGradient();
const STensor3& F0 = ipvprev->getConstRefToDeformationGradient();
/* data for J2 law */
IPTransverseIsoCurvature* q1 = ipvcur->getIPTransverseIsoCurvature();
const IPTransverseIsoCurvature* q0 = ipvprev->getIPTransverseIsoCurvature();
/* compute stress */
_tilaw.constitutive(F0,F1,ipvcur->getRefToFirstPiolaKirchhoffStress(),q0,q1,ipvcur->getRefToTangentModuli(),stiff);
ipvcur->setRefToDGElasticTangentModuli(this->elasticStiffness);
}
//Constructor for yarn
TransverseIsoYarnBDG3DMaterialLaw::TransverseIsoYarnBDG3DMaterialLaw(int lawnum,//tag of the material law
double rho, //material density
int physicnum,//physic number of
char* yarnInfo, //Yarn information file
int isunidir, //
double rDSvDT,//To initialize '_CriticalrDSvDT'
double MaxStrainRate, //The maximal value of the strain rate
double Cxx11,
double Cxx21,
double Cxx22,
double Cxx31,
double Cxx32,
double Cyy11,
double Cyy21,
double Cyy22,
double Cyy31,
double Cyy32,
double Cxy11,
double Cxy21,
double Cxy22,
double Cxy31,
double Cxy32,
double Cyz11,
double Cyz21,
double Cyz22,
double Cyz31,
double Cyz32,
double nuxy,
double nuxz,
double nuyz) :
dG3DMaterialLaw(lawnum,rho,true),
_tilaw(lawnum,rho, physicnum,yarnInfo,
isunidir, rDSvDT, MaxStrainRate,
Cxx11, Cxx21, Cxx22, Cxx31, Cxx32,
Cyy11, Cyy21, Cyy22, Cyy31, Cyy32,
Cxy11, Cxy21, Cxy22, Cxy31, Cxy32,
Cyz11, Cyz21, Cyz22, Cyz31, Cyz32,
nuxy, nuxz, nuyz)
{ double EArray[4] = {Cxx11,Cyy11,Cxy11,Cyz11};
int MaxEArrayIndex = maxArray_pos(EArray);
double E = EArray[MaxEArrayIndex];
double nu = 0.4;
fillElasticStiffness(E,nu,elasticStiffness);
cout<<"Constitutive law of physic volume "<< physicnum <<" has been built!" <<endl<<endl;
}
TransverseIsoYarnBDG3DMaterialLaw::TransverseIsoYarnBDG3DMaterialLaw(const TransverseIsoYarnBDG3DMaterialLaw &source) :
dG3DMaterialLaw(source), _tilaw(source._tilaw)
{
}
//Constructor for matrix
TransverseIsoYarnBDG3DMaterialLaw::TransverseIsoYarnBDG3DMaterialLaw(int lawnum,//tag of the material law
double rho, //material density
int physicnum,//physic number of
double rDSvDT,//To initialize '_CriticalrDSvDT'
double MaxStrainRate, //The maximal value of the strain rate
double Cxx11,
double Cxx21,
double Cxx22,
double Cxx31,
double Cxx32,
double Cyy11,
double Cyy21,
double Cyy22,
double Cyy31,
double Cyy32,
double Cxy11,
double Cxy21,
double Cxy22,
double Cxy31,
double Cxy32,
double Cyz11,
double Cyz21,
double Cyz22,
double Cyz31,
double Cyz32,
double nuxy,
double nuxz,
double nuyz) :
dG3DMaterialLaw(lawnum,rho,true),
_tilaw(lawnum,rho, physicnum,
rDSvDT, MaxStrainRate,
Cxx11, Cxx21, Cxx22, Cxx31, Cxx32,
Cyy11, Cyy21, Cyy22, Cyy31, Cyy32,
Cxy11, Cxy21, Cxy22, Cxy31, Cxy32,
Cyz11, Cyz21, Cyz22, Cyz31, Cyz32,
nuxy, nuxz, nuyz)
{
double EArray[4] = {Cxx11,Cyy11,Cxy11,Cyz11};
int MaxEArrayIndex = maxArray_pos(EArray);
double E = EArray[MaxEArrayIndex];
double nu = 0.4;
fillElasticStiffness(E,nu,elasticStiffness);
cout<<"Constitutive law of physic volume "<< physicnum <<" has been built!" <<endl<<endl;
}
void TransverseIsoYarnBDG3DMaterialLaw::createIPState(IPStateBase* &ips,
const bool* state_,
const MElement *ele,
const int nbFF,
const IntPt *GP, const int gpt) const
{
// check interface element
bool inter=true;
const MInterfaceElement *iele = dynamic_cast<const MInterfaceElement*>(ele);
if(iele==NULL) inter=false;
SPoint3 pt_Global;
ele->pnt(GP[gpt].pt[0],GP[gpt].pt[1],GP[gpt].pt[2],pt_Global);
SVector3 GaussP(pt_Global);
IPVariable* ipvi;
IPVariable* ipv1;
IPVariable* ipv2;
int isUniDir = _tilaw.get_IsUniDir();
switch(isUniDir)
{
case 0:
ipvi = new TransverseIsoYarnBDG3DIPVariable(GaussP, _tilaw.get_Yarn(), inter);
ipv1 = new TransverseIsoYarnBDG3DIPVariable(GaussP, _tilaw.get_Yarn(), inter);
ipv2 = new TransverseIsoYarnBDG3DIPVariable(GaussP, _tilaw.get_Yarn(), inter);
break;
case 1:
ipvi = new TransverseIsoYarnBDG3DIPVariable(GaussP, _tilaw.get_UniDir(), inter);
ipv1 = new TransverseIsoYarnBDG3DIPVariable(GaussP, _tilaw.get_UniDir(), inter);
ipv2 = new TransverseIsoYarnBDG3DIPVariable(GaussP, _tilaw.get_UniDir(), inter);
break;
case 2:
ipvi = new TransverseIsoYarnBDG3DIPVariable(GaussP, inter);
ipv1 = new TransverseIsoYarnBDG3DIPVariable(GaussP, inter);
ipv2 = new TransverseIsoYarnBDG3DIPVariable(GaussP, inter);
break;
}
if(ips != NULL) delete ips;
ips = new IP3State(state_,ipvi,ipv1,ipv2);
}
void TransverseIsoYarnBDG3DMaterialLaw::createIPVariable(IPVariable* &ipv,
const MElement *ele,
const int nbFF,
const IntPt *GP,
const int gpt) const
{
if(ipv !=NULL) delete ipv;
bool inter=true;
const MInterfaceElement *iele = dynamic_cast<const MInterfaceElement*>(ele);
if(iele == NULL){
inter=false;
}
SPoint3 pt_Global;
ele->pnt(GP[gpt].pt[0],GP[gpt].pt[1],GP[gpt].pt[2],pt_Global);
SVector3 GaussP(pt_Global);
int isUniDir = _tilaw.get_IsUniDir();
switch(isUniDir)
{
case 0:
ipv = new TransverseIsoYarnBDG3DIPVariable(GaussP, _tilaw.get_Yarn(), inter);
break;
case 1:
ipv = new TransverseIsoYarnBDG3DIPVariable(GaussP, _tilaw.get_UniDir(), inter);
break;
case 2:
ipv = new TransverseIsoYarnBDG3DIPVariable(GaussP, inter); break;
}
}
void TransverseIsoYarnBDG3DMaterialLaw::stress(IPVariable* ipv, const IPVariable* ipvp, const bool stiff, const bool checkfrac)
{
//get ipvariable
TransverseIsoYarnBDG3DIPVariable* ipvcur;
const TransverseIsoYarnBDG3DIPVariable* ipvprev;
FractureCohesive3DIPVariable* ipvtmp = dynamic_cast<FractureCohesive3DIPVariable*>(ipv);
if(ipvtmp !=NULL)
{
ipvcur = static_cast<TransverseIsoYarnBDG3DIPVariable*>(ipvtmp->getIPvBulk());
const FractureCohesive3DIPVariable *ipvtmp2 = static_cast<const FractureCohesive3DIPVariable*>(ipvp);
ipvprev = static_cast<const TransverseIsoYarnBDG3DIPVariable*>(ipvtmp2->getIPvBulk());
}
else
{
ipvcur = static_cast<TransverseIsoYarnBDG3DIPVariable*>(ipv);
ipvprev = static_cast<const TransverseIsoYarnBDG3DIPVariable*>(ipvp);
}
//strain xyz
const STensor3& F1 = ipvcur->getConstRefToDeformationGradient();
//const STensor3& F0 = ipvprev->getConstRefToDeformationGradient();
//data for J2 law
IPTransverseIsoYarnB* q1 = ipvcur->getIPTransverseIsoYarnB();
const IPTransverseIsoYarnB* q0 = ipvprev->getIPTransverseIsoYarnB();
// compute stress
_tilaw.constitutive(F1,ipvcur->getRefToFirstPiolaKirchhoffStress(),q0,q1,
ipvcur->getRefToTangentModuli(), stiff);
ipvcur->setRefToDGElasticTangentModuli(this->elasticStiffness);
}
AnisotropicDG3DMaterialLaw::AnisotropicDG3DMaterialLaw(
const int num,const double rho,
const double Ex,const double Ey,const double Ez,
const double Vxy,const double Vxz,const double Vyz,
const double MUxy,const double MUxz,const double MUyz,
const double euler[3]) :
dG3DMaterialLaw(num,rho,true),
_tilaw(num,rho,Ex,Ey,Ez,Vxy,Vxz,Vyz,
MUxy,MUxz,MUyz,euler[0],euler[1],euler[2]),
_ISTensor3(1.)
{
double nu = _tilaw.poissonRatio();
double mu = _tilaw.shearModulus();
double E = 2.*mu*(1.+nu); //Where does this come from? To Have the E corresponding to the Vmax, MUmax
//In the linear & homogeneous & isotropic case ?
fillElasticStiffness(E, nu, elasticStiffness); // WTF here? not sure to get it ...
}
AnisotropicDG3DMaterialLaw::AnisotropicDG3DMaterialLaw(
const int num,const double rho,
const double Ex,const double Ey,const double Ez,
const double Vxy,const double Vxz,const double Vyz,
const double MUxy,const double MUxz,const double MUyz,
const double alpha, const double beta, const double gamma) :
dG3DMaterialLaw(num,rho,true),
_tilaw(num,rho,Ex,Ey,Ez,Vxy,Vxz,Vyz,
MUxy,MUxz,MUyz,alpha, beta, gamma),
_ISTensor3(1.)
{
double nu = _tilaw.poissonRatio(); double mu = _tilaw.shearModulus();
double E = 2.*mu*(1.+nu); //Where does this come from? To Have the E corresponding to the Vmax, MUmax
//In the linear & homogeneous & isotropic case ?
fillElasticStiffness(E, nu, elasticStiffness); // WTF here? not sure to get it ...
}
AnisotropicDG3DMaterialLaw::AnisotropicDG3DMaterialLaw(const AnisotropicDG3DMaterialLaw &source) :
dG3DMaterialLaw(source),
_tilaw(source._tilaw),
_ISTensor3(source._ISTensor3)
{
}
void AnisotropicDG3DMaterialLaw::createIPState(IPStateBase* &ips,const bool* state_,const MElement *ele, const int nbFF_, const IntPt *GP, const int gpt) const
{
// check interface element
bool inter=true;
const MInterfaceElement *iele = dynamic_cast<const MInterfaceElement*>(ele);
if(iele==NULL) inter=false;
IPVariable* ipvi = new AnisotropicDG3DIPVariable(inter);
IPVariable* ipv1 = new AnisotropicDG3DIPVariable(inter);
IPVariable* ipv2 = new AnisotropicDG3DIPVariable(inter);
if(ips != NULL) delete ips;
ips = new IP3State(state_,ipvi,ipv1,ipv2);
}
void AnisotropicDG3DMaterialLaw::createIPVariable(IPVariable* &ipv,const MElement *ele, const int nbFF_, const IntPt *GP, const int gpt) const
{
if(ipv !=NULL) delete ipv;
bool inter=true;
const MInterfaceElement *iele = dynamic_cast<const MInterfaceElement*>(ele);
if(iele == NULL){
inter=false;
}
ipv = new AnisotropicDG3DIPVariable(inter);
}
void AnisotropicDG3DMaterialLaw::stress(IPVariable* ipv, const IPVariable* ipvp, const bool stiff, const bool checkfrac)
{
/* get ipvariable */
AnisotropicDG3DIPVariable* ipvcur;
const AnisotropicDG3DIPVariable* ipvprev;
FractureCohesive3DIPVariable* ipvtmp = dynamic_cast<FractureCohesive3DIPVariable*>(ipv);
if(ipvtmp !=NULL)
{
ipvcur = static_cast<AnisotropicDG3DIPVariable*>(ipvtmp->getIPvBulk());
const FractureCohesive3DIPVariable *ipvtmp2 = static_cast<const FractureCohesive3DIPVariable*>(ipvp);
ipvprev = static_cast<const AnisotropicDG3DIPVariable*>(ipvtmp2->getIPvBulk());
}
else
{
ipvcur = static_cast<AnisotropicDG3DIPVariable*>(ipv);
ipvprev = static_cast<const AnisotropicDG3DIPVariable*>(ipvp);
}
/* strain xyz */
const STensor3& F1 = ipvcur->getConstRefToDeformationGradient();
const STensor3& F0 = ipvprev->getConstRefToDeformationGradient();
/* data for J2 law */
IPAnisotropic* q1 = ipvcur->getIPAnisotropic();
const IPAnisotropic* q0 = ipvprev->getIPAnisotropic();
/* compute stress */
_tilaw.constitutive(F0,F1,ipvcur->getRefToFirstPiolaKirchhoffStress(),q0,q1,ipvcur->getRefToTangentModuli(),stiff);
ipvcur->setRefToDGElasticTangentModuli(this->elasticStiffness);
}
AnisotropicStochDG3DMaterialLaw::AnisotropicStochDG3DMaterialLaw(const int num, const double rho, const double alpha,
const double beta, const double gamma, const double OrigX, const double OrigY,
const char *propName, const int intpl):
dG3DMaterialLaw(num,rho,true), _tilaw(num, rho, alpha, beta, gamma, OrigX, OrigY, propName, intpl)
{
double nu = _tilaw.poissonRatio();
double mu = _tilaw.shearModulus();
double E = 2.*mu*(1.+nu);
fillElasticStiffness(E, nu, elasticStiffness);
}
AnisotropicStochDG3DMaterialLaw::AnisotropicStochDG3DMaterialLaw(const AnisotropicStochDG3DMaterialLaw &source) :
dG3DMaterialLaw(source), _tilaw(source._tilaw){}
void AnisotropicStochDG3DMaterialLaw::createIPState(IPStateBase* &ips,const bool* state_,const MElement *ele, const int nbFF_, const IntPt *GP, const int gpt) const
{
// check interface element
bool inter=true;
const MInterfaceElement *iele = dynamic_cast<const MInterfaceElement*>(ele);
if(iele==NULL) inter=false;
SPoint3 pt_Global;
ele->pnt(GP[gpt].pt[0],GP[gpt].pt[1],GP[gpt].pt[2],pt_Global);
SVector3 GaussP(pt_Global);
IPVariable* ipvi = new AnisotropicStochDG3DIPVariable(GaussP, _tilaw.getExMat(), _tilaw.getEyMat(), _tilaw.getEzMat(),
_tilaw.getVxyMat(), _tilaw.getVxzMat(), _tilaw.getVyzMat(), _tilaw.getMUxyMat(), _tilaw.getMUxzMat(),
_tilaw.getMUyzMat(), _tilaw.getDeltaX(), _tilaw.getDeltaY(), _tilaw.getOrigX(), _tilaw.getOrigY(), _tilaw.getInterpoly());
IPVariable* ipv1 = new AnisotropicStochDG3DIPVariable(GaussP, _tilaw.getExMat(), _tilaw.getEyMat(), _tilaw.getEzMat(),
_tilaw.getVxyMat(), _tilaw.getVxzMat(), _tilaw.getVyzMat(), _tilaw.getMUxyMat(), _tilaw.getMUxzMat(),
_tilaw.getMUyzMat(), _tilaw.getDeltaX(), _tilaw.getDeltaY(), _tilaw.getOrigX(), _tilaw.getOrigY(), _tilaw.getInterpoly());
IPVariable* ipv2 = new AnisotropicStochDG3DIPVariable(GaussP, _tilaw.getExMat(), _tilaw.getEyMat(), _tilaw.getEzMat(),
_tilaw.getVxyMat(), _tilaw.getVxzMat(), _tilaw.getVyzMat(), _tilaw.getMUxyMat(), _tilaw.getMUxzMat(),
_tilaw.getMUyzMat(), _tilaw.getDeltaX(), _tilaw.getDeltaY(), _tilaw.getOrigX(), _tilaw.getOrigY(), _tilaw.getInterpoly());
if(ips != NULL) delete ips;
ips = new IP3State(state_,ipvi,ipv1,ipv2);
}
void AnisotropicStochDG3DMaterialLaw::createIPVariable(IPVariable* &ipv,const MElement *ele, const int nbFF_, const IntPt *GP, const int gpt) const
{
if(ipv !=NULL) delete ipv;
bool inter=true;
const MInterfaceElement *iele = dynamic_cast<const MInterfaceElement*>(ele);
if(iele == NULL){ inter=false; }
SPoint3 pt_Global;
ele->pnt(GP[gpt].pt[0],GP[gpt].pt[1],GP[gpt].pt[2],pt_Global);
SVector3 GaussP(pt_Global);
ipv = new AnisotropicStochDG3DIPVariable(GaussP, _tilaw.getExMat(), _tilaw.getEyMat(), _tilaw.getEzMat(),
_tilaw.getVxyMat(), _tilaw.getVxzMat(), _tilaw.getVyzMat(), _tilaw.getMUxyMat(), _tilaw.getMUxzMat(),
_tilaw.getMUyzMat(), _tilaw.getDeltaX(), _tilaw.getDeltaY(), _tilaw.getOrigX(), _tilaw.getOrigY(), _tilaw.getInterpoly());
}
void AnisotropicStochDG3DMaterialLaw::stress(IPVariable* ipv, const IPVariable* ipvp, const bool stiff, const bool checkfrac)
{
/* get ipvariable */
AnisotropicStochDG3DIPVariable* ipvcur;
const AnisotropicStochDG3DIPVariable* ipvprev;
FractureCohesive3DIPVariable* ipvtmp = dynamic_cast<FractureCohesive3DIPVariable*>(ipv);
if(ipvtmp !=NULL)
{
ipvcur = static_cast<AnisotropicStochDG3DIPVariable*>(ipvtmp->getIPvBulk());
const FractureCohesive3DIPVariable *ipvtmp2 = static_cast<const FractureCohesive3DIPVariable*>(ipvp);
ipvprev = static_cast<const AnisotropicStochDG3DIPVariable*>(ipvtmp2->getIPvBulk());
}
else
{
ipvcur = static_cast<AnisotropicStochDG3DIPVariable*>(ipv);
ipvprev = static_cast<const AnisotropicStochDG3DIPVariable*>(ipvp);
}
/* strain xyz */
const STensor3& F1 = ipvcur->getConstRefToDeformationGradient();
const STensor3& F0 = ipvprev->getConstRefToDeformationGradient();
/* data for J2 law */
IPAnisotropicStoch* q1 = ipvcur->getIPAnisotropicStoch();
const IPAnisotropicStoch* q0 = ipvprev->getIPAnisotropicStoch();
/* compute stress */
_tilaw.constitutive(F0,F1,ipvcur->getRefToFirstPiolaKirchhoffStress(),q0,q1,ipvcur->getRefToTangentModuli(),stiff);
ipvcur->setRefToDGElasticTangentModuli(this->elasticStiffness);
}
LinearThermoMechanicsDG3DMaterialLaw::LinearThermoMechanicsDG3DMaterialLaw(const int num,const double rho,const double Ex, const double Ey, const double Ez, const double Vxy,
const double Vxz,const double Vyz, const double MUxy, const double MUxz, const double MUyz,const double alpha,
const double beta, const double gamma ,const double t0,const double Kx,const double Ky, const double Kz,
const double alphax,const double alphay,const double alphaz,const double cp) :
dG3DMaterialLaw(num,rho,true)
{
_lawLinearTM = new mlawLinearThermoMechanics(num,rho,Ex,Ey,Ez,Vxy,Vxz,Vyz,MUxy,MUxz,MUyz,alpha,beta,gamma,t0,Kx,Ky,Kz,alphax,alphay,alphaz,cp);
double nu = _lawLinearTM->poissonRatio();
double mu = _lawLinearTM->shearModulus();
double E = 2.*mu*(1.+nu); //Where does this come from? To Have the E corresponding to the Vmax, MUmax
//In the linear & homogeneous & isotropic case ?
linearK = _lawLinearTM->getConductivityTensor();
Stiff_alphaDilatation=_lawLinearTM->getStiff_alphaDilatation();
fillElasticStiffness(E, nu, elasticStiffness);
}
LinearThermoMechanicsDG3DMaterialLaw::LinearThermoMechanicsDG3DMaterialLaw(const LinearThermoMechanicsDG3DMaterialLaw &source) :
dG3DMaterialLaw(source)
{
_lawLinearTM = NULL;
if (source._lawLinearTM != NULL){
_lawLinearTM = new mlawLinearThermoMechanics(*(source._lawLinearTM));
// _linearTM = new mlawLinearThermoMechanics(*(source._lawLinearTM)); or this
}
linearK = source.linearK;
Stiff_alphaDilatation=source.Stiff_alphaDilatation;
}
void LinearThermoMechanicsDG3DMaterialLaw::createIPState(IPStateBase* &ips,const bool* state_,const MElement *ele, const int nbFF_, const IntPt *GP, const int gpt) const
{
// check interface element
bool inter=true;
const MInterfaceElement *iele = dynamic_cast<const MInterfaceElement*>(ele);
if(iele==NULL) inter=false;
IPVariable* ipvi = new LinearThermoMechanicsDG3DIPVariable(_lawLinearTM->getInitialTemperature(),inter);
IPVariable* ipv1 = new LinearThermoMechanicsDG3DIPVariable(_lawLinearTM->getInitialTemperature(),inter);
IPVariable* ipv2 = new LinearThermoMechanicsDG3DIPVariable(_lawLinearTM->getInitialTemperature(),inter); if(ips != NULL) delete ips;
ips = new IP3State(state_,ipvi,ipv1,ipv2);
}
void LinearThermoMechanicsDG3DMaterialLaw::createIPVariable(IPVariable* &ipv,const MElement *ele, const int nbFF_, const IntPt *GP, const int gpt) const
{
if(ipv !=NULL) delete ipv;
bool inter=true;
const MInterfaceElement *iele = dynamic_cast<const MInterfaceElement*>(ele);
if(iele == NULL){
inter=false;
}
ipv = new LinearThermoMechanicsDG3DIPVariable(_lawLinearTM->getInitialTemperature(),inter);
}
double LinearThermoMechanicsDG3DMaterialLaw::getExtraDofStoredEnergyPerUnitField(double T) const
{
return _lawLinearTM->getExtraDofStoredEnergyPerUnitField(T);
}
double LinearThermoMechanicsDG3DMaterialLaw::getInitialExtraDofStoredEnergyPerUnitField() const{
return _lawLinearTM->getInitialExtraDofStoredEnergyPerUnitField();
};
void LinearThermoMechanicsDG3DMaterialLaw::stress(IPVariable* ipv, const IPVariable* ipvp, const bool stiff, const bool checkfrac)
{
/* get ipvariable */
LinearThermoMechanicsDG3DIPVariable* ipvcur;
const LinearThermoMechanicsDG3DIPVariable* ipvprev;
FractureCohesive3DIPVariable* ipvtmp = dynamic_cast<FractureCohesive3DIPVariable*>(ipv);
if(ipvtmp !=NULL)
{
ipvcur = static_cast<LinearThermoMechanicsDG3DIPVariable*>(ipvtmp->getIPvBulk());
const FractureCohesive3DIPVariable *ipvtmp2 = static_cast<const FractureCohesive3DIPVariable*>(ipvp);
ipvprev = static_cast<const LinearThermoMechanicsDG3DIPVariable*>(ipvtmp2->getIPvBulk());
}
else
{
ipvcur = static_cast<LinearThermoMechanicsDG3DIPVariable*>(ipv);
ipvprev = static_cast<const LinearThermoMechanicsDG3DIPVariable*>(ipvp);
}
/* strain xyz */
const STensor3& F1 = ipvcur->getConstRefToDeformationGradient();
const STensor3& F0 = ipvprev->getConstRefToDeformationGradient();
const double T = ipvcur->getConstRefToTemperature();
const double T0 = ipvprev->getConstRefToTemperature();
const SVector3& gradT= ipvcur->getConstRefToGradT();
STensor3& dPdT=ipvcur->getRefTodPdT();
SVector3& fluxT=ipvcur->getRefToThermalFlux();
STensor3& dqdgradT=ipvcur->getRefTodThermalFluxdGradT();
SVector3& dqdT=ipvcur->getRefTodThermalFluxdT();
STensor33& dqdF=ipvcur->getRefTodThermalFluxdF();
//const STensor3 *linearK()=ipvcur-> getConstRefTolinearK();
STensor3& P=ipvcur->getRefToFirstPiolaKirchhoffStress();
STensor43& Tangent=ipvcur->getRefToTangentModuli();
double &w = ipvcur->getRefToThermalSource();
double &dwdt = ipvcur->getRefTodThermalSourcedField();
STensor3 &dwdf = ipvcur->getRefTodThermalSourcedF();
const SVector3&N=ipvcur->getRefToReferenceOutwardNormal();
const SVector3& ujump=ipvcur->getRefToJump()(0);
const double &Tjump = ipvcur->getRefToFieldJump()(0);
double& Cp = ipvcur->getRefToExtraDofFieldCapacityPerUnitField()(0);
/* data for J2 law */
IPLinearThermoMechanics* q1 = ipvcur->getIPLinearThermoMechanics(); const IPLinearThermoMechanics* q0 = ipvprev->getIPLinearThermoMechanics();
_lawLinearTM-> constitutive(F0,F1, P,q0, q1,Tangent,T0,T, gradT,fluxT,dPdT,dqdgradT,dqdT,dqdF,w,dwdt,dwdf,stiff);
ipvcur->setRefToDGElasticTangentModuli(this->elasticStiffness);
ipvcur->setRefToLinearConductivity(this->linearK);
ipvcur->setRefToStiff_alphaDilatation(this->Stiff_alphaDilatation);
/* if (checkfrac==true) // to get the H1 norm in the interface
q1->_fracEnergy=defoEnergy(ujump,Tjump,N);
else
q1->_elasticEnergy=defoEnergy(F1, gradT);*/
Cp = _lawLinearTM->getExtraDofStoredEnergyPerUnitField(T);
}
double LinearThermoMechanicsDG3DMaterialLaw::defoEnergy(const STensor3& Fn,const SVector3& gradT) const // If i have sigma ...
{
static STensor3 defo;
for (int i=0; i<3; i++)
for (int j=0; j<3; j++)
defo(i,j)=(Fn(j,i)+Fn(i,j))*0.5;
defo(0,0)-=1.;
defo(1,1)-=1.;
defo(2,2)-=1.;
double En=0.;
for(int i=0;i<3;i++)
{
for(int j=0;j<3;j++)
{
// En+=defo(i,j)*defo(i,j);
}
En+=gradT(i)*gradT(i);
}
return En;
}
double LinearThermoMechanicsDG3DMaterialLaw::defoEnergy( const SVector3&ujump,const double &Tjump, const SVector3& N) const // If i have sigma ...
{
double e_interface=0.;
for(int i=0;i<3;i++)
{
// for(int j=0;j<3;j++)
// {
//e_interface+=ujump(i)*ujump(i);//+1./2.*Tjump*N(i)*N(i)*Tjump;
e_interface+=Tjump*N(i)*N(i)*Tjump;
}
// }
return fabs(e_interface);
}
J2ThermoMechanicsDG3DMaterialLaw::J2ThermoMechanicsDG3DMaterialLaw(const int num,const double rho,
double E, const double nu, const J2IsotropicHardening &_j2IH,
const double t0,const double Kther, const double alphather, const double cp,
const double tol) : dG3DMaterialLaw(num,rho,true), _j2law(num,E,nu,rho,_j2IH,t0,Kther,alphather,cp,tol)
{
fillElasticStiffness(E, nu, elasticStiffness);
linearK = _j2law.getConductivityTensor();
Stiff_alphaDilatation=_j2law.getStiff_alphaDilatation();
}
J2ThermoMechanicsDG3DMaterialLaw::J2ThermoMechanicsDG3DMaterialLaw(const J2ThermoMechanicsDG3DMaterialLaw &source) : dG3DMaterialLaw(source), _j2law(source._j2law)
{
linearK = source.linearK;
Stiff_alphaDilatation=source.Stiff_alphaDilatation;}
void J2ThermoMechanicsDG3DMaterialLaw::createIPState(IPStateBase* &ips,const bool* state_,const MElement *ele, const int nbFF_, const IntPt *GP, const int gpt) const
{
// check interface element
bool inter=true;
const MInterfaceElement *iele = dynamic_cast<const MInterfaceElement*>(ele);
if(iele==NULL) inter=false;
IPVariable* ipvi = new J2ThermoMechanicsDG3DIPVariable(_j2law.getInitialTemperature(),_j2law,inter);
IPVariable* ipv1 = new J2ThermoMechanicsDG3DIPVariable(_j2law.getInitialTemperature(),_j2law,inter);
IPVariable* ipv2 = new J2ThermoMechanicsDG3DIPVariable(_j2law.getInitialTemperature(),_j2law,inter);
if(ips != NULL) delete ips;
ips = new IP3State(state_,ipvi,ipv1,ipv2);
}
void J2ThermoMechanicsDG3DMaterialLaw::createIPVariable(IPVariable* &ipv,const MElement *ele, const int nbFF_, const IntPt *GP, const int gpt) const
{
if(ipv !=NULL) delete ipv;
bool inter=true;
const MInterfaceElement *iele = dynamic_cast<const MInterfaceElement*>(ele);
if(iele == NULL){
inter=false;
}
ipv = new J2ThermoMechanicsDG3DIPVariable(_j2law.getInitialTemperature(),_j2law,inter);
}
double J2ThermoMechanicsDG3DMaterialLaw::getExtraDofStoredEnergyPerUnitField(double T) const
{
return _j2law.getExtraDofStoredEnergyPerUnitField(T);
}
double J2ThermoMechanicsDG3DMaterialLaw::getInitialExtraDofStoredEnergyPerUnitField() const{
return _j2law.getInitialExtraDofStoredEnergyPerUnitField();
}
void J2ThermoMechanicsDG3DMaterialLaw::stress(IPVariable* ipv, const IPVariable* ipvp, const bool stiff, const bool checkfrac)
{
/* get ipvariable */
J2ThermoMechanicsDG3DIPVariable* ipvcur;
const J2ThermoMechanicsDG3DIPVariable* ipvprev;
FractureCohesive3DIPVariable* ipvtmp = dynamic_cast<FractureCohesive3DIPVariable*>(ipv);
if(ipvtmp !=NULL)
{
ipvcur = static_cast<J2ThermoMechanicsDG3DIPVariable*>(ipvtmp->getIPvBulk());
const FractureCohesive3DIPVariable *ipvtmp2 = static_cast<const FractureCohesive3DIPVariable*>(ipvp);
ipvprev = static_cast<const J2ThermoMechanicsDG3DIPVariable*>(ipvtmp2->getIPvBulk());
}
else
{
ipvcur = static_cast<J2ThermoMechanicsDG3DIPVariable*>(ipv);
ipvprev = static_cast<const J2ThermoMechanicsDG3DIPVariable*>(ipvp);
}
/* strain xyz */
const STensor3& F1 = ipvcur->getConstRefToDeformationGradient();
const STensor3& F0 = ipvprev->getConstRefToDeformationGradient();
const double T = ipvcur->getConstRefToTemperature();
const double T0 = ipvprev->getConstRefToTemperature();
const SVector3& gradT= ipvcur->getConstRefToGradT();
STensor3& dPdT=ipvcur->getRefTodPdT();
SVector3& fluxT=ipvcur->getRefToThermalFlux();
STensor3& dqdgradT=ipvcur->getRefTodThermalFluxdGradT();
SVector3& dqdT=ipvcur->getRefTodThermalFluxdT();
STensor33& dqdF=ipvcur->getRefTodThermalFluxdF();
//const STensor3 *linearK()=ipvcur-> getConstRefTolinearK();
STensor3& P=ipvcur->getRefToFirstPiolaKirchhoffStress();
STensor43& Tangent=ipvcur->getRefToTangentModuli(); double &w = ipvcur->getRefToThermalSource();
double &dwdt = ipvcur->getRefTodThermalSourcedField();
STensor3 &dwdf = ipvcur->getRefTodThermalSourcedF();
double& Cp = ipvcur->getRefToExtraDofFieldCapacityPerUnitField()(0);
/* data for J2 law */
IPJ2ThermoMechanics* q1 = ipvcur->getIPJ2ThermoMechanics();
const IPJ2ThermoMechanics* q0 = ipvprev->getIPJ2ThermoMechanics();
/* compute stress */
_j2law.constitutive(F0,F1, P,q0, q1,Tangent,T0,T, gradT,fluxT,dPdT,dqdgradT,dqdT,dqdF,w,dwdt,dwdf,stiff);
ipvcur->setRefToDGElasticTangentModuli(this->elasticStiffness);
ipvcur->setRefToLinearConductivity(this->linearK);
ipvcur->setRefToStiff_alphaDilatation(this->Stiff_alphaDilatation);
Cp = _j2law.getExtraDofStoredEnergyPerUnitField(T);
}
FullJ2ThermoMechanicsDG3DMaterialLaw::FullJ2ThermoMechanicsDG3DMaterialLaw(const int num, const double E, const double nu,const double rho,
const double sy0,const double h,const double tol,
const bool matrixbyPerturbation, const double pert):
dG3DMaterialLaw(num,rho),_j2FullThermo(num,E,nu,rho,sy0,h,tol,matrixbyPerturbation,pert,true){
fillElasticStiffness(E, nu, elasticStiffness);
};
FullJ2ThermoMechanicsDG3DMaterialLaw::FullJ2ThermoMechanicsDG3DMaterialLaw(const FullJ2ThermoMechanicsDG3DMaterialLaw &src):
dG3DMaterialLaw(src),_j2FullThermo(src._j2FullThermo){};
void FullJ2ThermoMechanicsDG3DMaterialLaw::setStrainOrder(const int i){
_j2FullThermo.setStrainOrder(i);
};
void FullJ2ThermoMechanicsDG3DMaterialLaw::setReferenceTemperature(const double Tr){
_j2FullThermo.setReferenceTemperature(Tr);
};
void FullJ2ThermoMechanicsDG3DMaterialLaw::setIsotropicHardeningLaw(const J2IsotropicHardening& isoHard){
_j2FullThermo.setIsotropicHardeningLaw(isoHard);
};
void FullJ2ThermoMechanicsDG3DMaterialLaw::setReferenceThermalExpansionCoefficient(const double al){
_j2FullThermo.setReferenceThermalExpansionCoefficient(al);
};
void FullJ2ThermoMechanicsDG3DMaterialLaw::setTemperatureFunction_ThermalExpansionCoefficient(const scalarFunction& Tfunc){
_j2FullThermo.setTemperatureFunction_ThermalExpansionCoefficient(Tfunc);
};
void FullJ2ThermoMechanicsDG3DMaterialLaw::setReferenceThermalConductivity(const double KK){
_j2FullThermo.setReferenceThermalConductivity(KK);
};
void FullJ2ThermoMechanicsDG3DMaterialLaw::setTemperatureFunction_ThermalConductivity(const scalarFunction& Tfunc){
_j2FullThermo.setTemperatureFunction_ThermalConductivity(Tfunc);
};
void FullJ2ThermoMechanicsDG3DMaterialLaw::setReferenceCp(const double Cp){
_j2FullThermo.setReferenceCp(Cp);
};
void FullJ2ThermoMechanicsDG3DMaterialLaw::setTemperatureFunction_Cp(const scalarFunction& Tfunc){
_j2FullThermo.setTemperatureFunction_Cp(Tfunc);
};
void FullJ2ThermoMechanicsDG3DMaterialLaw::setTemperatureFunction_Hardening(const scalarFunction& Tfunc){
_j2FullThermo.setTemperatureFunction_Hardening(Tfunc);
};
void FullJ2ThermoMechanicsDG3DMaterialLaw::setTemperatureFunction_InitialYieldStress(const scalarFunction& Tfunc){
_j2FullThermo.setTemperatureFunction_InitialYieldStress(Tfunc);
};
void FullJ2ThermoMechanicsDG3DMaterialLaw::setTemperatureFunction_BulkModulus(const scalarFunction& Tfunc){
_j2FullThermo.setTemperatureFunction_BulkModulus(Tfunc);
};void FullJ2ThermoMechanicsDG3DMaterialLaw::setTemperatureFunction_ShearModulus(const scalarFunction& Tfunc){
_j2FullThermo.setTemperatureFunction_ShearModulus(Tfunc);
};
void FullJ2ThermoMechanicsDG3DMaterialLaw::setTaylorQuineyFactor(const double f){
_j2FullThermo.setTaylorQuineyFactor(f);
};
void FullJ2ThermoMechanicsDG3DMaterialLaw::setThermalEstimationPreviousConfig(const bool flag){
_j2FullThermo.setThermalEstimationPreviousConfig(flag);
}
void FullJ2ThermoMechanicsDG3DMaterialLaw::createIPState(IPStateBase* &ips,const bool* state_,const MElement *ele, const int nbFF_, const IntPt *GP, const int gpt) const
{
bool inter=true;
const MInterfaceElement *iele = dynamic_cast<const MInterfaceElement*>(ele);
if(iele==NULL) inter=false;
IPVariable* ipvi = new J2FullThermoMechanicsDG3DIPVariable(_j2FullThermo.getJ2IsotropicHardening(),inter);
IPVariable* ipv1 = new J2FullThermoMechanicsDG3DIPVariable(_j2FullThermo.getJ2IsotropicHardening(),inter);
IPVariable* ipv2 = new J2FullThermoMechanicsDG3DIPVariable(_j2FullThermo.getJ2IsotropicHardening(),inter);
if(ips != NULL) delete ips;
ips = new IP3State(state_,ipvi,ipv1,ipv2);
}
void FullJ2ThermoMechanicsDG3DMaterialLaw::createIPVariable(IPVariable* &ipv,const MElement *ele, const int nbFF_, const IntPt *GP, const int gpt) const
{
bool inter=true;
const MInterfaceElement *iele = dynamic_cast<const MInterfaceElement*>(ele);
if(iele==NULL) inter=false;
if(ipv !=NULL) delete ipv;
ipv = new J2FullThermoMechanicsDG3DIPVariable(_j2FullThermo.getJ2IsotropicHardening(),inter);
}
void FullJ2ThermoMechanicsDG3DMaterialLaw::stress(IPVariable* ipv, const IPVariable* ipvp, const bool stiff, const bool checkfrac)
{
J2FullThermoMechanicsDG3DIPVariable* ipvcur = static_cast<J2FullThermoMechanicsDG3DIPVariable*>(ipv);
const J2FullThermoMechanicsDG3DIPVariable* ipvprev = static_cast<const J2FullThermoMechanicsDG3DIPVariable*>(ipvp);
/* strain xyz */
const STensor3& F1 = ipvcur->getConstRefToDeformationGradient();
const STensor3& F0 = ipvprev->getConstRefToDeformationGradient();
const double T = ipvcur->getConstRefToTemperature();
const double T0 = ipvprev->getConstRefToTemperature();
const SVector3& gradT0= ipvprev->getConstRefToGradT();
const SVector3& gradT= ipvcur->getConstRefToGradT();
STensor3& dPdT=ipvcur->getRefTodPdT();
SVector3& fluxT=ipvcur->getRefToThermalFlux();
STensor3& dqdgradT=ipvcur->getRefTodThermalFluxdGradT();
SVector3& dqdT=ipvcur->getRefTodThermalFluxdT();
STensor33& dqdF=ipvcur->getRefTodThermalFluxdF();
STensor3& P=ipvcur->getRefToFirstPiolaKirchhoffStress();
STensor43& Tangent=ipvcur->getRefToTangentModuli();
double &w = ipvcur->getRefToThermalSource();
double &dwdt = ipvcur->getRefTodThermalSourcedField();
STensor3 &dwdf = ipvcur->getRefTodThermalSourcedF();
double& mechSource = ipvcur->getRefToMechanicalSource()(0);
double& dmechSourcedt = ipvcur->getRefTodMechanicalSourcedField()(0,0);
STensor3 & dmechSourcedf = ipvcur->getRefTodMechanicalSourcedF()[0];
double& Cp = ipvcur->getRefToExtraDofFieldCapacityPerUnitField()(0);
/* data for J2 law */
IPJ2ThermoMechanics* q1 = ipvcur->getIPJ2ThermoMechanics();
const IPJ2ThermoMechanics* q0 = ipvprev->getIPJ2ThermoMechanics();
/* compute stress */
_j2FullThermo.constitutive(F0,F1,P,q0,q1,Tangent,T0,T,gradT0,gradT,fluxT,dPdT,dqdgradT,dqdT,dqdF,
w,dwdt,dwdf,mechSource,dmechSourcedt,dmechSourcedf,stiff);
ipvcur->setRefToDGElasticTangentModuli(this->elasticStiffness);
ipvcur->setRefToLinearConductivity(_j2FullThermo.getConductivityTensor());
ipvcur->setRefToStiff_alphaDilatation(_j2FullThermo.getStiff_alphaDilatation());
Cp = _j2FullThermo.getExtraDofStoredEnergyPerUnitField(T);
}
double FullJ2ThermoMechanicsDG3DMaterialLaw::getInitialExtraDofStoredEnergyPerUnitField() const{
return _j2FullThermo.getInitialExtraDofStoredEnergyPerUnitField();
};
double FullJ2ThermoMechanicsDG3DMaterialLaw::getExtraDofStoredEnergyPerUnitField(double T) const
{
return _j2FullThermo.getExtraDofStoredEnergyPerUnitField(T);
}
double FullJ2ThermoMechanicsDG3DMaterialLaw::scaleFactor() const {
return _j2FullThermo.shearModulus();
};
mlawAnIsotropicTherMechDG3DMaterialLaw::mlawAnIsotropicTherMechDG3DMaterialLaw(const int num,const double E, const double nu,const double rho,const double EA, const double GA, const double nu_minor,
const double Ax, const double Ay, const double Az, const double alpha, const double beta, const double gamma,const double t0,
const double kx,const double ky,const double kz,const double cp,const double alphath):
dG3DMaterialLaw(num,rho,true)
{
_lawAnTM = new mlawAnIsotropicTherMech(num,E,nu,rho,EA,GA,nu_minor,Ax,Ay,Az,alpha, beta,gamma,t0, kx, ky,kz,cp, alphath);
linearK = _lawAnTM->getConductivityTensor();
Stiff_alphaDilatation=_lawAnTM->getStiff_alphaDilatation();
fillElasticStiffness(E, nu, elasticStiffness);
}
mlawAnIsotropicTherMechDG3DMaterialLaw::mlawAnIsotropicTherMechDG3DMaterialLaw(const mlawAnIsotropicTherMechDG3DMaterialLaw &source) :
dG3DMaterialLaw (source)// , dG3DMaterialLaw (source)
{
_lawAnTM = NULL;
if (source._lawAnTM !=NULL){
//_lawTM = new mlawAnIsotropicTherMech(*(source._lawTM));
_lawAnTM = new mlawAnIsotropicTherMech(*(source._lawAnTM));
// _linearTM = new mlawLinearThermoMechanics(*(source._lawLinearTM)); or this
}
linearK = source.linearK;
Stiff_alphaDilatation=source. Stiff_alphaDilatation;
}
void mlawAnIsotropicTherMechDG3DMaterialLaw::createIPState(IPStateBase* &ips,const bool* state_,const MElement *ele, const int nbFF_, const IntPt *GP, const int gpt) const
{
// check interface element
bool inter=true;
const MInterfaceElement *iele = dynamic_cast<const MInterfaceElement*>(ele);
if(iele==NULL) inter=false;
IPVariable* ipvi = new AnIsotropicTherMechDG3DIPVariable(_lawAnTM->getInitialTemperature(),inter);
IPVariable* ipv1 = new AnIsotropicTherMechDG3DIPVariable(_lawAnTM->getInitialTemperature(),inter);
IPVariable* ipv2 = new AnIsotropicTherMechDG3DIPVariable(_lawAnTM->getInitialTemperature(),inter);
if(ips != NULL) delete ips;
ips = new IP3State(state_,ipvi,ipv1,ipv2);
}
void mlawAnIsotropicTherMechDG3DMaterialLaw::createIPVariable(IPVariable* &ipv,const MElement *ele, const int nbFF_, const IntPt *GP, const int gpt) const
{
if(ipv !=NULL) delete ipv;
bool inter=true;
const MInterfaceElement *iele = dynamic_cast<const MInterfaceElement*>(ele);
if(iele == NULL){
inter=false;
} ipv = new AnIsotropicTherMechDG3DIPVariable(_lawAnTM->getInitialTemperature(),inter);
}
void mlawAnIsotropicTherMechDG3DMaterialLaw::stress(IPVariable* ipv, const IPVariable* ipvp, const bool stiff, const bool checkfrac)
{
AnIsotropicTherMechDG3DIPVariable* ipvcur;
const AnIsotropicTherMechDG3DIPVariable* ipvprev;
FractureCohesive3DIPVariable* ipvtmp = dynamic_cast<FractureCohesive3DIPVariable*>(ipv);
if(ipvtmp !=NULL)
{
ipvcur = static_cast<AnIsotropicTherMechDG3DIPVariable*>(ipvtmp->getIPvBulk());
const FractureCohesive3DIPVariable *ipvtmp2 = static_cast<const FractureCohesive3DIPVariable*>(ipvp);
ipvprev = static_cast<const AnIsotropicTherMechDG3DIPVariable*>(ipvtmp2->getIPvBulk());
}
else
{
ipvcur = static_cast<AnIsotropicTherMechDG3DIPVariable*>(ipv);
ipvprev = static_cast<const AnIsotropicTherMechDG3DIPVariable*>(ipvp);
}
const STensor3& F1 = ipvcur->getConstRefToDeformationGradient();
const STensor3& F0 = ipvprev->getConstRefToDeformationGradient();
const double T = ipvcur->getConstRefToTemperature();
const double T0 = ipvprev->getConstRefToTemperature();
const SVector3& gradT= ipvcur->getConstRefToGradT();
STensor3& dPdT=ipvcur->getRefTodPdT();
SVector3& fluxT=ipvcur->getRefToThermalFlux();
STensor3& dqdgradT=ipvcur->getRefTodThermalFluxdGradT();
SVector3& dqdT=ipvcur->getRefTodThermalFluxdT();
STensor33& dqdF=ipvcur->getRefTodThermalFluxdF();
//const STensor3 *linearK()=ipvcur-> getConstRefTolinearK();
STensor3& P=ipvcur->getRefToFirstPiolaKirchhoffStress();
STensor43& Tangent=ipvcur->getRefToTangentModuli();
double &w = ipvcur->getRefToThermalSource();
double &dwdt = ipvcur->getRefTodThermalSourcedField();
STensor3 &dwdf = ipvcur->getRefTodThermalSourcedF();
IPAnIsotropicTherMech* q1 = ipvcur->getIPAnIsotropicTherMech();
const IPAnIsotropicTherMech* q0 = ipvprev->getIPAnIsotropicTherMech();
_lawAnTM->constitutive(F0,F1, P,q0, q1,Tangent,T0,T, gradT,fluxT,dPdT,dqdgradT,dqdT,dqdF,w,dwdt,dwdf,stiff);
ipvcur->setRefToLinearConductivity(this->linearK);
ipvcur->setRefToDGElasticTangentModuli(this->elasticStiffness);
ipvcur->setRefToStiff_alphaDilatation(this->Stiff_alphaDilatation);
}
SMPDG3DMaterialLaw::SMPDG3DMaterialLaw(const int num,const double rho,const double alpha, const double beta, const double gamma ,const double t0, const double Kx,const double Ky, const double Kz,
const double mu_groundState3,const double Im3, const double pi,const double Tr, const double Nmu_groundState2 ,const double mu_groundState2 ,
const double Im2,const double Sgl2,const double Sr2,const double Delta,const double m2, const double epsilonr ,const double n,const double epsilonp02,
const double alphar1,const double alphagl1,const double Ggl1,const double Gr1 ,const double Mgl1,const double Mr1,const double Mugl1,
const double Mur1,const double epsilon01,const double Qgl1,const double Qr1,const double epsygl1 ,const double d1,
const double m1,const double V1,const double alphap, const double Sa0,const double ha1,const double b1,const double g1,
const double phia01,const double Z1, const double r1,const double s1,const double Sb01,const double Hgl1,const double Lgl1,const double Hr1,
const double Lr1,const double l1,const double Kb,const double be1,const double c0,const double wp,const double c1):
dG3DMaterialLaw (num,rho, true)
{
_lawTMSMP = new mlawSMP( num,rho,alpha,beta,gamma,t0,Kx,Ky,Kz,mu_groundState3,Im3,pi,Tr,Nmu_groundState2, mu_groundState2, Im2,Sgl2, Sr2, Delta, m2,epsilonr,n,epsilonp02, alphar1, alphagl1, Ggl1, Gr1, Mgl1, Mr1, Mugl1, Mur1, epsilon01, Qgl1,Qr1, epsygl1 , d1, m1, V1, alphap,
Sa0, ha1, b1, g1, phia01, Z1,r1, s1, Sb01, Hgl1, Lgl1, Hr1, Lr1, l1, Kb, be1,c0,wp,c1);
double nu = Mr1;
double mu = Gr1;
double E = 2.*mu*(1.+nu); //Where does this come from? To Have the E corresponding to the Vmax, MUmax
//In the linear & homogeneous & isotropic case ?
linearK = _lawTMSMP->getInitialConductivityTensor();
Stiff_alphaDilatation=_lawTMSMP->getStiff_alphaDilatation();
fillElasticStiffness(E, nu, elasticStiffness);
}
SMPDG3DMaterialLaw::SMPDG3DMaterialLaw(const SMPDG3DMaterialLaw &source) :
dG3DMaterialLaw (source)// , dG3DMaterialLaw (source)
{
_lawTMSMP = NULL;
if (source._lawTMSMP != NULL){
_lawTMSMP = new mlawSMP(*(source._lawTMSMP));
// _linearTM = new mlawLinearThermoMechanics(*(source._lawLinearTM)); or this
}
linearK = source.linearK;
Stiff_alphaDilatation = source.Stiff_alphaDilatation;
}
void SMPDG3DMaterialLaw::setStrainOrder(const int order ){
_lawTMSMP->setStrainOrder(order);
};
void SMPDG3DMaterialLaw::setMechanism2(bool mechanism2) {
_lawTMSMP->setMechanism2(mechanism2);
};
double SMPDG3DMaterialLaw::getExtraDofStoredEnergyPerUnitField(double T) const
{
return _lawTMSMP->getExtraDofStoredEnergyPerUnitField(T);
}
void SMPDG3DMaterialLaw::createIPState(IPStateBase* &ips,const bool* state_,const MElement *ele, const int nbFF_, const IntPt *GP, const int gpt) const
{
// check interface element
bool inter=true;
const MInterfaceElement *iele = dynamic_cast<const MInterfaceElement*>(ele);
if(iele==NULL) inter=false;
IPVariable* ipvi = new SMPDG3DIPVariable(_lawTMSMP->getInitialTemperature(),*_lawTMSMP,inter);
IPVariable* ipv1 = new SMPDG3DIPVariable(_lawTMSMP->getInitialTemperature(),*_lawTMSMP,inter);
IPVariable* ipv2 = new SMPDG3DIPVariable(_lawTMSMP->getInitialTemperature(),*_lawTMSMP,inter);
if(ips != NULL) delete ips;
ips = new IP3State(state_,ipvi,ipv1,ipv2);
}
void SMPDG3DMaterialLaw::createIPVariable(IPVariable* &ipv,const MElement *ele, const int nbFF_, const IntPt *GP, const int gpt) const
{
if(ipv !=NULL) delete ipv;
bool inter=true;
const MInterfaceElement *iele = dynamic_cast<const MInterfaceElement*>(ele);
if(iele == NULL){
inter=false;
}
ipv = new SMPDG3DIPVariable(_lawTMSMP->getInitialTemperature(),*_lawTMSMP,inter);
}
void SMPDG3DMaterialLaw::stress(IPVariable* ipv, const IPVariable* ipvp, const bool stiff, const bool checkfrac)
{
/* get ipvariable */
SMPDG3DIPVariable* ipvcur;
const SMPDG3DIPVariable* ipvprev;
FractureCohesive3DIPVariable* ipvtmp = dynamic_cast<FractureCohesive3DIPVariable*>(ipv);
if(ipvtmp !=NULL)
{
ipvcur = static_cast<SMPDG3DIPVariable*>(ipvtmp->getIPvBulk());
const FractureCohesive3DIPVariable *ipvtmp2 = static_cast<const FractureCohesive3DIPVariable*>(ipvp);
ipvprev = static_cast<const SMPDG3DIPVariable*>(ipvtmp2->getIPvBulk());
}
else
{
ipvcur = static_cast<SMPDG3DIPVariable*>(ipv);
ipvprev = static_cast<const SMPDG3DIPVariable*>(ipvp);
}
/* strain xyz */
const STensor3& F1 = ipvcur->getConstRefToDeformationGradient();
const STensor3& F0 = ipvprev->getConstRefToDeformationGradient();
const double temperature = ipvcur->getConstRefToTemperature();
const double T0 = ipvprev->getConstRefToTemperature();
const SVector3& gradT= ipvcur->getConstRefToGradT();
STensor3& dPdT=ipvcur->getRefTodPdT();
SVector3& fluxT=ipvcur->getRefToThermalFlux();
STensor3& dqdgradT=ipvcur->getRefTodThermalFluxdGradT();
SVector3& dqdT=ipvcur->getRefTodThermalFluxdT();
STensor33& dqdF=ipvcur->getRefTodThermalFluxdF();
STensor3& P=ipvcur->getRefToFirstPiolaKirchhoffStress();
STensor43& Tangent=ipvcur->getRefToTangentModuli();
double &w = ipvcur->getRefToThermalSource();
double &dwdt = ipvcur->getRefTodThermalSourcedField();
STensor3& dwdf = ipvcur->getRefTodThermalSourcedF();
/* data for J2 law */
IPSMP* q1 = ipvcur->getIPSMP();
const IPSMP* q0 = ipvprev->getIPSMP();
_lawTMSMP-> constitutive(F0,F1, P,q0, q1,Tangent,T0,temperature, gradT,fluxT,dPdT,dqdgradT,dqdT,dqdF,stiff,w,dwdt,dwdf);
ipvcur->setRefToDGElasticTangentModuli(this->elasticStiffness);
ipvcur->setRefToLinearConductivity(this->linearK);
ipvcur->setRefToStiff_alphaDilatation(this->Stiff_alphaDilatation);
}
PhenomenologicalSMPDG3DMaterialLaw::PhenomenologicalSMPDG3DMaterialLaw(const int num, const double rho, const double alpha, const double beta, const double gamma ,const double t0, const double Kxg,const double Kyg, const double Kzg, const double Kxr, const double Kyr, const double Kzr,const double w, const double TaylorQuiney, const double c, const double cp, const double cg, const double cr,const double Tg,const double DeltaTg, const double Eg, const double Er, const double nug, const double nur,const double alphag, const double alphar, const double zg0, const double tau_y0):
dG3DMaterialLaw (num,rho, true)
{
_lawTMSMP = new mlawPhenomenologicalSMP(num,rho,alpha,beta,gamma,t0,Kxg,Kyg,Kzg,Kxr,Kyr,Kzr,w,TaylorQuiney,c,cp,cg,cr,Tg,DeltaTg,Eg,Er,nug,nur,alphag,alphar,zg0,tau_y0);
double nu = 0.3; //to be updated
double mu = 1.e6; //to be updated
double E = 2.*mu*(1.+nu); //Where does this come from? To Have the E corresponding to the Vmax, MUmax
//In the linear & homogeneous & isotropic case ?
linearK = _lawTMSMP->getInitialConductivityTensor();
Stiff_alphaDilatation=_lawTMSMP->getStiff_alphaDilatation();
fillElasticStiffness(E, nu, elasticStiffness);
}
PhenomenologicalSMPDG3DMaterialLaw::PhenomenologicalSMPDG3DMaterialLaw(const PhenomenologicalSMPDG3DMaterialLaw &source) :
dG3DMaterialLaw (source)// , dG3DMaterialLaw (source)
{
_lawTMSMP = NULL;
if (source._lawTMSMP != NULL){
_lawTMSMP = new mlawPhenomenologicalSMP(*(source._lawTMSMP));
}
linearK = source.linearK;
Stiff_alphaDilatation = source.Stiff_alphaDilatation;}
void PhenomenologicalSMPDG3DMaterialLaw::setStrainOrder(const int order ){
_lawTMSMP->setStrainOrder(order);
};
double PhenomenologicalSMPDG3DMaterialLaw::getExtraDofStoredEnergyPerUnitField(double zg) const
{
return _lawTMSMP->getExtraDofStoredEnergyPerUnitField(zg);
}
void PhenomenologicalSMPDG3DMaterialLaw::createIPState(IPStateBase* &ips,const bool* state_,const MElement *ele, const int nbFF_, const IntPt *GP, const int gpt) const
{
// check interface element
bool inter=true;
const MInterfaceElement *iele = dynamic_cast<const MInterfaceElement*>(ele);
if(iele==NULL) inter=false;
IPVariable* ipvi = new PhenomenologicalSMPDG3DIPVariable(_lawTMSMP->getInitialTemperature(),*_lawTMSMP,inter);
IPVariable* ipv1 = new PhenomenologicalSMPDG3DIPVariable(_lawTMSMP->getInitialTemperature(),*_lawTMSMP,inter);
IPVariable* ipv2 = new PhenomenologicalSMPDG3DIPVariable(_lawTMSMP->getInitialTemperature(),*_lawTMSMP,inter);
if(ips != NULL) delete ips;
ips = new IP3State(state_,ipvi,ipv1,ipv2);
}
void PhenomenologicalSMPDG3DMaterialLaw::createIPVariable(IPVariable* &ipv,const MElement *ele, const int nbFF_, const IntPt *GP, const int gpt) const
{
if(ipv !=NULL) delete ipv;
bool inter=true;
const MInterfaceElement *iele = dynamic_cast<const MInterfaceElement*>(ele);
if(iele == NULL){
inter=false;
}
ipv = new PhenomenologicalSMPDG3DIPVariable(_lawTMSMP->getInitialTemperature(),*_lawTMSMP,inter);
}
void PhenomenologicalSMPDG3DMaterialLaw::stress(IPVariable* ipv, const IPVariable* ipvp, const bool stiff, const bool checkfrac)
{
/* get ipvariable */
PhenomenologicalSMPDG3DIPVariable* ipvcur; //= static_cast < nonLocalDamageDG3DIPVariable *> (ipv);
const PhenomenologicalSMPDG3DIPVariable* ipvprev; //= static_cast <const nonLocalDamageDG3DIPVariable *> (ipvp);
FractureCohesive3DIPVariable* ipvtmp = dynamic_cast<FractureCohesive3DIPVariable*>(ipv);
if(ipvtmp !=NULL)
{
ipvcur = static_cast<PhenomenologicalSMPDG3DIPVariable*>(ipvtmp->getIPvBulk());
const FractureCohesive3DIPVariable *ipvtmp2 = static_cast<const FractureCohesive3DIPVariable*>(ipvp);
ipvprev = static_cast<const PhenomenologicalSMPDG3DIPVariable*>(ipvtmp2->getIPvBulk());
}
else
{
ipvcur = static_cast<PhenomenologicalSMPDG3DIPVariable*>(ipv);
ipvprev = static_cast<const PhenomenologicalSMPDG3DIPVariable*>(ipvp);
}
/* strain xyz */
const STensor3& F1 = ipvcur->getConstRefToDeformationGradient();
const STensor3& F0 = ipvprev->getConstRefToDeformationGradient();
const double temperature = ipvcur->getConstRefToTemperature();
const double T0 = ipvprev->getConstRefToTemperature();
const SVector3& gradT= ipvcur->getConstRefToGradT();
STensor3& dPdT=ipvcur->getRefTodPdT();
SVector3& fluxT=ipvcur->getRefToThermalFlux();
STensor3& dqdgradT=ipvcur->getRefTodThermalFluxdGradT();
SVector3& dqdT=ipvcur->getRefTodThermalFluxdT();
STensor33& dqdF=ipvcur->getRefTodThermalFluxdF();
STensor3& P=ipvcur->getRefToFirstPiolaKirchhoffStress(); STensor43& Tangent=ipvcur->getRefToTangentModuli();
double &w = ipvcur->getRefToThermalSource();
double &dwdt = ipvcur->getRefTodThermalSourcedField();
STensor3& dwdf = ipvcur->getRefTodThermalSourcedF();
/* data for J2 law */
IPPhenomenologicalSMP* q1 = ipvcur->getIPPhenomenologicalSMP();
const IPPhenomenologicalSMP* q0 = ipvprev->getIPPhenomenologicalSMP();
_lawTMSMP->constitutive(F0,F1, P,q0, q1,Tangent,T0,temperature, gradT,fluxT,dPdT,dqdgradT,dqdT,dqdF,stiff,w,dwdt,dwdf);
ipvcur->setRefToDGElasticTangentModuli(this->elasticStiffness);
ipvcur->setRefToLinearConductivity(this->linearK);
ipvcur->setRefToStiff_alphaDilatation(this->Stiff_alphaDilatation);
}
//
LinearElecTherMechDG3DMaterialLaw::LinearElecTherMechDG3DMaterialLaw(const int num,const double rho,const double Ex, const double Ey, const double Ez, const double Vxy,
const double Vxz, const double Vyz,const double MUxy, const double MUxz, const double MUyz, const double alpha, const double beta, const double gamma,
const double t0,const double Kx,const double Ky,const double Kz,const double alphax,const double alphay,const double alphaz,const double lx,const double ly,
const double lz,const double seebeck,double cp,const double v0):
dG3DMaterialLaw(num,rho,true)
{
_lawLinearETM = new mlawLinearElecTherMech(num,rho, Ex, Ey, Ez, Vxy, Vxz, Vyz, MUxy, MUxz, MUyz,alpha, beta,gamma, t0, Kx, Ky,Kz, alphax, alphay, alphaz,lx,ly,lz,seebeck,cp,v0);
double nu = _lawLinearETM->poissonRatio();
double mu = _lawLinearETM->shearModulus();
double E = 2.*mu*(1.+nu); //Where does this come from? To Have the E corresponding to the Vmax, MUmax
//In the linear & homogeneous & isotropic case ?
fillElasticStiffness(E, nu, elasticStiffness);
Stiff_alphaDilatation=_lawLinearETM->getStiff_alphaDilatation();
}
LinearElecTherMechDG3DMaterialLaw::LinearElecTherMechDG3DMaterialLaw(const LinearElecTherMechDG3DMaterialLaw &source) :
dG3DMaterialLaw(source)
{
_lawLinearETM = NULL;
if (source._lawLinearETM != NULL){
_lawLinearETM = new mlawLinearElecTherMech(*(source._lawLinearETM));
// _linearETM = new mlawLinearElecTherMech(*(source._lawLinearETM)); or this
}
Stiff_alphaDilatation=source.Stiff_alphaDilatation;
}
void LinearElecTherMechDG3DMaterialLaw::createIPState(IPStateBase* &ips,const bool* state_,const MElement *ele, const int nbFF_, const IntPt *GP, const int gpt) const
{
// check interface element
bool inter=true;
const MInterfaceElement *iele = dynamic_cast<const MInterfaceElement*>(ele);
if(iele==NULL) inter=false;
IPVariable* ipvi = new LinearElecTherMechDG3DIPVariable(_lawLinearETM->getInitialTemperature(),_lawLinearETM->getInitialVoltage(),inter);
IPVariable* ipv1 = new LinearElecTherMechDG3DIPVariable(_lawLinearETM->getInitialTemperature(),_lawLinearETM->getInitialVoltage(),inter);
IPVariable* ipv2 = new LinearElecTherMechDG3DIPVariable(_lawLinearETM->getInitialTemperature(),_lawLinearETM->getInitialVoltage(),inter);
if(ips != NULL) delete ips;
ips = new IP3State(state_,ipvi,ipv1,ipv2);
}
void LinearElecTherMechDG3DMaterialLaw::createIPVariable(IPVariable* &ipv,const MElement *ele, const int nbFF_, const IntPt *GP, const int gpt) const
{
if(ipv !=NULL) delete ipv;
bool inter=true;
const MInterfaceElement *iele = dynamic_cast<const MInterfaceElement*>(ele);
if(iele == NULL){
inter=false;
}
ipv = new LinearElecTherMechDG3DIPVariable(_lawLinearETM->getInitialTemperature(),_lawLinearETM->getInitialVoltage(),inter);
}
void LinearElecTherMechDG3DMaterialLaw::stress(IPVariable* ipv, const IPVariable* ipvp, const bool stiff, const bool checkfrac)
{
/* get ipvariable */
LinearElecTherMechDG3DIPVariable* ipvcur;
const LinearElecTherMechDG3DIPVariable* ipvprev;
FractureCohesive3DIPVariable* ipvtmp = dynamic_cast<FractureCohesive3DIPVariable*>(ipv);
if(ipvtmp !=NULL)
{
ipvcur = static_cast<LinearElecTherMechDG3DIPVariable*>(ipvtmp->getIPvBulk());
const FractureCohesive3DIPVariable *ipvtmp2 = static_cast<const FractureCohesive3DIPVariable*>(ipvp);
ipvprev = static_cast<const LinearElecTherMechDG3DIPVariable*>(ipvtmp2->getIPvBulk());
}
else
{
ipvcur = static_cast<LinearElecTherMechDG3DIPVariable*>(ipv);
ipvprev = static_cast<const LinearElecTherMechDG3DIPVariable*>(ipvp);
}
/* strain xyz */
const STensor3& Fn = ipvcur->getConstRefToDeformationGradient();
const STensor3& F0 = ipvprev->getConstRefToDeformationGradient();
const double temperature = ipvcur->getConstRefToTemperature();
const double T0 = ipvprev->getConstRefToTemperature();
const double voltage = ipvcur->getConstRefToVoltage();
const SVector3& gradT= ipvcur->getConstRefToGradT();
const SVector3& gradV= ipvcur->getConstRefToGradV();
STensor3& dPdT=ipvcur->getRefTodPdT();
STensor3& dPdV=ipvcur->getRefTodPdV();
SVector3& fluxT=ipvcur->getRefToFluxT();
SVector3& fluxV=ipvcur->getRefToFluxV();
STensor3& dqdgradT=ipvcur->getRefTodFluxTdGradT();
STensor3& dqdgradV=ipvcur->getRefTodFluxTdGradV();
STensor3& dedgradV=ipvcur->getRefTodFluxVdGradV();
STensor3& dedgradT=ipvcur->getRefTodFluxVdGradT();
SVector3& dqdT=ipvcur->getRefTodFluxTdT();
SVector3& dqdV=ipvcur->getRefTodFluxTdV();
SVector3& dedV=ipvcur->getRefTodFluxVdV();
SVector3& dedT=ipvcur->getRefTodFluxVdT();
STensor33& dqdF=ipvcur->getRefTodFluxTdF();
STensor33& dedF=ipvcur->getRefTodFluxVdF();
STensor3& P=ipvcur->getRefToFirstPiolaKirchhoffStress();
STensor43& Tangent=ipvcur->getRefToTangentModuli();
double &w = ipvcur->getRefToThermalSource();
double &dwdt = ipvcur->getRefTodThermalSourcedT();
STensor3& dwdf = ipvcur->getRefTodThermalSourcedF();
STensor3& k10 = ipvcur->getRefToTherConductivity();
STensor3& k20 = ipvcur->getRefToCrossTherConductivity();
STensor3& l10 = ipvcur->getRefToElecConductivity();
STensor3& l20 = ipvcur->getRefToCrossElecConductivity();
STensor3& dk10dT = ipvcur->getRefTodTherConductivitydT();
STensor3& dk20dT = ipvcur->getRefTodCrossTherConductivitydT();
STensor3& dl10dT = ipvcur->getRefTodElecConductivitydT();
STensor3& dl20dT = ipvcur->getRefTodCrossElecConductivitydT();
STensor3& dk10dv = ipvcur->getRefTodTherConductivitydv();
STensor3& dl20dv = ipvcur->getRefTodCrossElecConductivitydv();
STensor43& dl10dF=ipvcur->getRefTodElecConductivitydF();
STensor43& dl20dF=ipvcur->getRefTodCrossElecConductivitydF();
STensor43& dk10dF=ipvcur->getRefTodTherConductivitydF();
STensor43& dk20dF=ipvcur->getRefTodCrossTherConductivitydF();
SVector3 &fluxjy=ipvcur->getRefToFluxEnergy();
SVector3 &djydV=ipvcur->getRefTodFluxEnergydV();
STensor3 &djydgradV=ipvcur->getRefTodFluxEnergydGradV();
SVector3 &djydT=ipvcur->getRefTodFluxEnergydT();
STensor3 &djydgradT=ipvcur->getRefTodFluxEnergydGradT(); STensor33 &djydF=ipvcur->getRefTodFluxEnergydF();
STensor3 &jy1=ipvcur->getRefToEnergyConductivity();
STensor3 &djy1dT= ipvcur->getRefTodEnergyConductivitydT();
STensor3 &djy1dV= ipvcur->getRefTodEnergyConductivitydV();
STensor43 &djy1dF=ipvcur->getRefTodEnergyConductivitydF();
/* data for J2 law */
IPLinearElecTherMech* q1 = ipvcur->getIPLinearElecTherMech();
const IPLinearElecTherMech* q0 = ipvprev->getIPLinearElecTherMech();
_lawLinearETM->constitutive( F0,Fn,P, q0, q1, Tangent,T0, temperature, voltage,gradT, gradV, fluxT,fluxV, dPdT,dPdV, dqdgradT,dqdT, dqdgradV,dqdV,
dedV,dedgradV,dedT,dedgradT,dqdF,dedF,w,dwdt,dwdf,l10,l20,k10,k20,dl10dT,dl20dT,dk10dT,dk20dT,dl20dv,dk10dv,dl10dF,dl20dF,dk10dF,dk20dF
,fluxjy,djydV,djydgradV,djydT,djydgradT,djydF,jy1,djy1dT,djy1dV,djy1dF,stiff);
ipvcur->setRefToDGElasticTangentModuli(this->elasticStiffness);
ipvcur->setRefToStiff_alphaDilatation(this->Stiff_alphaDilatation);
}
mlawAnIsotropicElecTherMechDG3DMaterialLaw::mlawAnIsotropicElecTherMechDG3DMaterialLaw(const int num,const double E, const double nu,const double rho,const double EA, const double GA, const double nu_minor,
const double Ax, const double Ay, const double Az, const double alpha, const double beta, const double gamma,const double t0,const double lx,
const double ly,const double lz,const double seebeck,const double kx,const double ky,
const double kz,const double cp,const double alphath,const double v0):
dG3DMaterialLaw(num,rho,true)
{
_lawETM = new mlawAnIsotropicElecTherMech(num,E,nu,rho,EA,GA,nu_minor,Ax,Ay,Az,alpha, beta,gamma,t0, lx,ly,lz,seebeck,kx, ky,kz,cp, alphath,v0);
Stiff_alphaDilatation=_lawETM->getStiff_alphaDilatation();
fillElasticStiffness(E, nu, elasticStiffness);
}
mlawAnIsotropicElecTherMechDG3DMaterialLaw::mlawAnIsotropicElecTherMechDG3DMaterialLaw(const mlawAnIsotropicElecTherMechDG3DMaterialLaw &source) :
dG3DMaterialLaw (source)// , dG3DMaterialLaw (source)
{
_lawETM = NULL;
if (source._lawETM != NULL){
_lawETM = new mlawAnIsotropicElecTherMech(*(source._lawETM));
// _linearTM = new mlawLinearThermoMechanics(*(source._lawLinearTM)); or this
}
Stiff_alphaDilatation=source. Stiff_alphaDilatation;
}
void mlawAnIsotropicElecTherMechDG3DMaterialLaw::createIPState(IPStateBase* &ips,const bool* state_,const MElement *ele, const int nbFF_, const IntPt *GP, const int gpt) const
{
// check interface element
bool inter=true;
const MInterfaceElement *iele = dynamic_cast<const MInterfaceElement*>(ele);
if(iele==NULL) inter=false;
IPVariable* ipvi = new AnIsotropicElecTherMechDG3DIPVariable(_lawETM->getInitialTemperature(),_lawETM->getInitialVoltage(),*_lawETM,inter);
IPVariable* ipv1 = new AnIsotropicElecTherMechDG3DIPVariable(_lawETM->getInitialTemperature(),_lawETM->getInitialVoltage(),*_lawETM,inter);
IPVariable* ipv2 = new AnIsotropicElecTherMechDG3DIPVariable(_lawETM->getInitialTemperature(),_lawETM->getInitialVoltage(),*_lawETM,inter);
if(ips != NULL) delete ips;
ips = new IP3State(state_,ipvi,ipv1,ipv2);
}
void mlawAnIsotropicElecTherMechDG3DMaterialLaw::createIPVariable(IPVariable* &ipv,const MElement *ele, const int nbFF_, const IntPt *GP, const int gpt) const
{
if(ipv !=NULL) delete ipv;
bool inter=true;
const MInterfaceElement *iele = dynamic_cast<const MInterfaceElement*>(ele);
if(iele == NULL){
inter=false;
} ipv = new AnIsotropicElecTherMechDG3DIPVariable(_lawETM->getInitialTemperature(),_lawETM->getInitialVoltage(),*_lawETM,inter);
}
void mlawAnIsotropicElecTherMechDG3DMaterialLaw::stress(IPVariable* ipv, const IPVariable* ipvp, const bool stiff, const bool checkfrac)
{
AnIsotropicElecTherMechDG3DIPVariable* ipvcur;
const AnIsotropicElecTherMechDG3DIPVariable* ipvprev;
FractureCohesive3DIPVariable* ipvtmp = dynamic_cast<FractureCohesive3DIPVariable*>(ipv);
if(ipvtmp !=NULL)
{
ipvcur = static_cast<AnIsotropicElecTherMechDG3DIPVariable*>(ipvtmp->getIPvBulk());
const FractureCohesive3DIPVariable *ipvtmp2 = static_cast<const FractureCohesive3DIPVariable*>(ipvp);
ipvprev = static_cast<const AnIsotropicElecTherMechDG3DIPVariable*>(ipvtmp2->getIPvBulk());
}
else
{
ipvcur = static_cast<AnIsotropicElecTherMechDG3DIPVariable*>(ipv);
ipvprev = static_cast<const AnIsotropicElecTherMechDG3DIPVariable*>(ipvp);
}
const STensor3& Fn = ipvcur->getConstRefToDeformationGradient();
const STensor3& F0 = ipvprev->getConstRefToDeformationGradient();
const double temperature = ipvcur->getConstRefToTemperature();
const double T0 = ipvprev->getConstRefToTemperature();
const double voltage = ipvcur->getConstRefToVoltage();
const SVector3& gradT= ipvcur->getConstRefToGradT();
const SVector3& gradV= ipvcur->getConstRefToGradV();
STensor3& dPdT =ipvcur->getRefTodPdT();
STensor3& dPdV =ipvcur->getRefTodPdV();
SVector3& fluxT=ipvcur->getRefToFluxT();
SVector3& fluxV=ipvcur->getRefToFluxV();
STensor3& dqdgradT =ipvcur->getRefTodFluxTdGradT();
STensor3& dqdgradV =ipvcur->getRefTodFluxTdGradV();
STensor3& djedgradV=ipvcur->getRefTodFluxVdGradV();
STensor3& djedgradT=ipvcur->getRefTodFluxVdGradT();
SVector3& dqdT =ipvcur->getRefTodFluxTdT();
SVector3& dqdV =ipvcur->getRefTodFluxTdV();
SVector3& djedV =ipvcur->getRefTodFluxVdV();
SVector3& djedT =ipvcur->getRefTodFluxVdT();
STensor33& dqdF =ipvcur->getRefTodFluxTdF();
STensor33& djedF=ipvcur->getRefTodFluxVdF();
STensor3& P=ipvcur->getRefToFirstPiolaKirchhoffStress();
STensor43& Tangent=ipvcur->getRefToTangentModuli();
double &w = ipvcur->getRefToThermalSource();
double &dwdT = ipvcur->getRefTodThermalSourcedT();
STensor3& dwdf = ipvcur->getRefTodThermalSourcedF();
STensor3& k10 = ipvcur->getRefToTherConductivity();
STensor3& k20 = ipvcur->getRefToCrossTherConductivity();
STensor3& l10 = ipvcur->getRefToElecConductivity();
STensor3& l20 = ipvcur->getRefToCrossElecConductivity();
STensor3& dk10dT = ipvcur->getRefTodTherConductivitydT();
STensor3& dk20dT = ipvcur->getRefTodCrossTherConductivitydT();
STensor3& dl10dT = ipvcur->getRefTodElecConductivitydT();
STensor3& dl20dT = ipvcur->getRefTodCrossElecConductivitydT();
STensor3& dk10dv = ipvcur->getRefTodTherConductivitydv();
STensor3& dl20dv= ipvcur->getRefTodCrossElecConductivitydv();
STensor43& dl10dF=ipvcur->getRefTodElecConductivitydF();
STensor43& dl20dF=ipvcur->getRefTodCrossElecConductivitydF();
STensor43& dk10dF=ipvcur->getRefTodTherConductivitydF();
STensor43& dk20dF=ipvcur->getRefTodCrossTherConductivitydF();
SVector3 &fluxjy=ipvcur->getRefToFluxEnergy();
SVector3 &djydV=ipvcur->getRefTodFluxEnergydV();
STensor3 &djydgradV=ipvcur->getRefTodFluxEnergydGradV(); SVector3 &djydT=ipvcur->getRefTodFluxEnergydT();
STensor3 &djydgradT=ipvcur->getRefTodFluxEnergydGradT();
STensor33 &djydF=ipvcur->getRefTodFluxEnergydF();
STensor3 &jy1=ipvcur->getRefToEnergyConductivity();
STensor3 &djy1dT= ipvcur->getRefTodEnergyConductivitydT();
STensor3 &djy1dV= ipvcur->getRefTodEnergyConductivitydV();
STensor43 &djy1dF=ipvcur->getRefTodEnergyConductivitydF();
IPAnIsotropicElecTherMech* q1 = ipvcur->getIPAnIsotropicElecTherMech();
const IPAnIsotropicElecTherMech* q0 = ipvprev->getIPAnIsotropicElecTherMech();
_lawETM->constitutive( F0,Fn,P, q0, q1, Tangent, T0, temperature, voltage,gradT, gradV, fluxT,fluxV, dPdT,dPdV, dqdgradT,dqdT, dqdgradV,dqdV,
djedV,djedgradV,djedT,djedgradT,dqdF,djedF,w,dwdT,dwdf,l10,l20,k10,k20,dl10dT,dl20dT,dk10dT,dk20dT,
dl20dv,dk10dv,dl10dF,dl20dF,dk10dF,dk20dF,fluxjy,djydV,djydgradV,djydT,djydgradT,djydF,jy1,djy1dT,djy1dV,djy1dF,stiff);
ipvcur->setRefToDGElasticTangentModuli(this->elasticStiffness);
ipvcur->setRefToStiff_alphaDilatation(this->Stiff_alphaDilatation);
}
mlawElecSMPDG3DMaterialLaw::mlawElecSMPDG3DMaterialLaw(const int num,const double rho, const double alpha, const double beta, const double gamma,const double t0,const double Kx,const double Ky,
const double Kz, const double mu_groundState3,const double Im3, const double pi ,const double Tr,
const double Nmu_groundState2,const double mu_groundState2, const double Im2,const double Sgl2,const double Sr2,const double Delta,const double m2, const double epsilonr,
const double n, const double epsilonp02,const double alphar1,const double alphagl1,const double Ggl1,const double Gr1,const double Mgl1,const double Mr1,const double Mugl1
,const double Mur1, const double epsilon01,const double QglOnKb1,const double QrOnKb1,const double epsygl1 ,const double d1,const double m1,const double VOnKb1,
const double alphap,const double Sa0,const double ha1,const double b1,const double g1,const double phia01,const double Z1,const double r1,const double s1,const double Sb01
,const double Hgl1,const double Lgl1,const double Hr1,const double Lr1,const double l1,const double Kb,const double be1,const double c0,const double wp,const double c1
,const double lx,const double ly,const double lz,const double seebeck,const double v0):
dG3DMaterialLaw(num,rho,true)
{
_lawETMSMP = new mlawElecSMP( num,rho,alpha,beta,gamma,t0,Kx,Ky,Kz,mu_groundState3,Im3,pi,
Tr,Nmu_groundState2, mu_groundState2, Im2,Sgl2, Sr2, Delta, m2,epsilonr,n,epsilonp02, alphar1, alphagl1, Ggl1, Gr1, Mgl1, Mr1, Mugl1, Mur1, epsilon01, QglOnKb1,
QrOnKb1, epsygl1 , d1, m1, VOnKb1, alphap, Sa0, ha1, b1, g1, phia01, Z1,r1, s1, Sb01, Hgl1, Lgl1, Hr1, Lr1, l1, Kb, be1,c0,wp,c1,lx,ly,lz,seebeck,v0);
double nu = Mur1;
double mu = Gr1 ;
double E = 2.*mu*(1.+nu); //Where does this come from? To Have the E corresponding to the Vmax, MUmax
//In the linear & homogeneous & isotropic case ?
fillElasticStiffness(E, nu, elasticStiffness);
Stiff_alphaDilatation= _lawETMSMP->getStiff_alphaDilatation();
}
mlawElecSMPDG3DMaterialLaw::mlawElecSMPDG3DMaterialLaw(const mlawElecSMPDG3DMaterialLaw &source) :
dG3DMaterialLaw (source)// , dG3DMaterialLaw (source)
{
_lawETMSMP = NULL;
if (source._lawETMSMP != NULL){
_lawETMSMP = new mlawElecSMP(*(source._lawETMSMP));
}
Stiff_alphaDilatation=source. Stiff_alphaDilatation;
// _linearTM = new mlawLinearThermoMechanics(*(source._lawLinearTM)); or this
/* K10 = source.K10;
K20 = source.K20;
L10 = source.L10;
L20 = source.L20;*/
}
void mlawElecSMPDG3DMaterialLaw::createIPState(IPStateBase* &ips,const bool* state_,const MElement *ele, const int nbFF_, const IntPt *GP, const int gpt) const
{ // check interface element
bool inter=true;
const MInterfaceElement *iele = dynamic_cast<const MInterfaceElement*>(ele);
if(iele==NULL) inter=false;
IPVariable* ipvi = new ElecSMPDG3DIPVariable(_lawETMSMP->getInitialTemperature(),_lawETMSMP->getInitialVoltage(),*_lawETMSMP,inter);
IPVariable* ipv1 = new ElecSMPDG3DIPVariable(_lawETMSMP->getInitialTemperature(),_lawETMSMP->getInitialVoltage(),*_lawETMSMP,inter);
IPVariable* ipv2 = new ElecSMPDG3DIPVariable(_lawETMSMP->getInitialTemperature(),_lawETMSMP->getInitialVoltage(),*_lawETMSMP,inter);
if(ips != NULL) delete ips;
ips = new IP3State(state_,ipvi,ipv1,ipv2);
}
void mlawElecSMPDG3DMaterialLaw::createIPVariable(IPVariable* &ipv,const MElement *ele, const int nbFF_, const IntPt *GP, const int gpt) const
{
if(ipv !=NULL) delete ipv;
bool inter=true;
const MInterfaceElement *iele = dynamic_cast<const MInterfaceElement*>(ele);
if(iele == NULL){
inter=false;
}
ipv = new ElecSMPDG3DIPVariable(_lawETMSMP->getInitialTemperature(),_lawETMSMP->getInitialVoltage(),*_lawETMSMP,inter);
}
void mlawElecSMPDG3DMaterialLaw::stress(IPVariable* ipv, const IPVariable* ipvp, const bool stiff, const bool checkfrac)
{
ElecSMPDG3DIPVariable* ipvcur;
const ElecSMPDG3DIPVariable* ipvprev;
FractureCohesive3DIPVariable* ipvtmp = dynamic_cast<FractureCohesive3DIPVariable*>(ipv);
if(ipvtmp !=NULL)
{
ipvcur = static_cast<ElecSMPDG3DIPVariable*>(ipvtmp->getIPvBulk());
const FractureCohesive3DIPVariable *ipvtmp2 = static_cast<const FractureCohesive3DIPVariable*>(ipvp);
ipvprev = static_cast<const ElecSMPDG3DIPVariable*>(ipvtmp2->getIPvBulk());
}
else
{
ipvcur = static_cast<ElecSMPDG3DIPVariable*>(ipv);
ipvprev = static_cast<const ElecSMPDG3DIPVariable*>(ipvp);
}
const STensor3& Fn = ipvcur->getConstRefToDeformationGradient();
const STensor3& F0 = ipvprev->getConstRefToDeformationGradient();
const double temperature = ipvcur->getConstRefToTemperature();
const double T0 = ipvprev->getConstRefToTemperature();
const double voltage = ipvcur->getConstRefToVoltage();
const SVector3& gradT= ipvcur->getConstRefToGradT();
const SVector3& gradV= ipvcur->getConstRefToGradV();
STensor3& dPdT = ipvcur->getRefTodPdT();
STensor3& dPdV = ipvcur->getRefTodPdV();
SVector3& fluxT= ipvcur->getRefToFluxT();
SVector3& fluxV= ipvcur->getRefToFluxV();
STensor3& dqdgradT= ipvcur->getRefTodFluxTdGradT();
STensor3& dqdgradV= ipvcur->getRefTodFluxTdGradV();
STensor3& dedgradV= ipvcur->getRefTodFluxVdGradV();
STensor3& dedgradT= ipvcur->getRefTodFluxVdGradT();
SVector3& dqdT= ipvcur->getRefTodFluxTdT();
SVector3& dqdV= ipvcur->getRefTodFluxTdV();
SVector3& dedV= ipvcur->getRefTodFluxVdV();
SVector3& dedT= ipvcur->getRefTodFluxVdT();
STensor33& dqdF= ipvcur->getRefTodFluxTdF();
STensor33& dedF= ipvcur->getRefTodFluxVdF();
STensor3& P = ipvcur->getRefToFirstPiolaKirchhoffStress();
STensor43& Tangent=ipvcur->getRefToTangentModuli();
double &w = ipvcur->getRefToThermalSource();
double &dwdT = ipvcur->getRefTodThermalSourcedT();
STensor3& dwdf = ipvcur->getRefTodThermalSourcedF();
STensor3& k10 = ipvcur->getRefToTherConductivity();
STensor3& k20 = ipvcur->getRefToCrossTherConductivity();
STensor3& l10 = ipvcur->getRefToElecConductivity();
STensor3& l20 = ipvcur->getRefToCrossElecConductivity();
STensor3& dk10dT = ipvcur->getRefTodTherConductivitydT();
STensor3& dk20dT = ipvcur->getRefTodCrossTherConductivitydT();
STensor3& dl10dT = ipvcur->getRefTodElecConductivitydT();
STensor3& dl20dT = ipvcur->getRefTodCrossElecConductivitydT();
STensor3& dk10dv = ipvcur->getRefTodTherConductivitydv();
STensor3& dl20dv = ipvcur->getRefTodCrossElecConductivitydv();
STensor43& dl10dF= ipvcur->getRefTodElecConductivitydF();
STensor43& dl20dF= ipvcur->getRefTodCrossElecConductivitydF();
STensor43& dk10dF= ipvcur->getRefTodTherConductivitydF();
STensor43& dk20dF= ipvcur->getRefTodCrossTherConductivitydF();
SVector3 &fluxjy = ipvcur->getRefToFluxEnergy();
SVector3 &djydV = ipvcur->getRefTodFluxEnergydV();
STensor3 &djydgradV= ipvcur->getRefTodFluxEnergydGradV();
SVector3 &djydT = ipvcur->getRefTodFluxEnergydT();
STensor3 &djydgradT= ipvcur->getRefTodFluxEnergydGradT();
STensor33 &djydF = ipvcur->getRefTodFluxEnergydF();
STensor3 &jy1 = ipvcur->getRefToEnergyConductivity();
STensor3 &djy1dT = ipvcur->getRefTodEnergyConductivitydT();
STensor3 &djy1dV = ipvcur->getRefTodEnergyConductivitydV();
STensor43 &djy1dF = ipvcur->getRefTodEnergyConductivitydF();
/* IPElecSMP* q1 = ipvcur->getIPElecSMP();
const IPElecSMP* q0 = ipvprev->getIPElecSMP();*/
IPElecSMP* q1 = ipvcur->getIPElecSMP();
const IPElecSMP* q0 = ipvprev->getIPElecSMP();
_lawETMSMP->constitutive( F0,Fn,P, q0, q1, Tangent, T0, temperature, voltage,gradT, gradV, fluxT,fluxV, dPdT,dPdV, dqdgradT,dqdT, dqdgradV,dqdV,
dedV,dedgradV,dedT,dedgradT,dqdF,dedF,w,dwdT,dwdf,l10,l20,k10,k20,dl10dT,dl20dT,dk10dT,dk20dT
,dl20dv,dk10dv,dl10dF,dl20dF,dk10dF,dk20dF,fluxjy,djydV,djydgradV,djydT,djydgradT,djydF,jy1,djy1dT,djy1dV,djy1dF,stiff);
ipvcur->setRefToDGElasticTangentModuli(this->elasticStiffness);
ipvcur->setRefToStiff_alphaDilatation(this->Stiff_alphaDilatation);
}
void mlawElecSMPDG3DMaterialLaw::setStrainOrder(const int order ){
_lawETMSMP->setStrainOrder(order);
};
void mlawElecSMPDG3DMaterialLaw::setMechanism2(bool mechanism2)
{
_lawETMSMP->setMechanism2(mechanism2);
};