From f15c04cd43409a6b37805bd8c1345431f7369a63 Mon Sep 17 00:00:00 2001 From: FLE_Knight <ujwalkishore.jinaga@mail.polimi.it> Date: Wed, 27 Apr 2022 14:21:30 +0200 Subject: [PATCH] April27 --- .../materialLaw/mlawHyperelastic.cpp | 50 ++++++++----------- .../materialLaw/mlawHyperelastic.h | 11 ++-- .../mlawNonLocalDamageHyperelastic.cpp | 6 +-- dG3D/src/dG3DIPVariable.cpp | 4 +- dG3D/src/dG3DIPVariable.h | 5 +- dG3D/src/dG3DMaterialLaw.cpp | 12 ++--- dG3D/src/dG3DMaterialLaw.h | 4 +- 7 files changed, 42 insertions(+), 50 deletions(-) diff --git a/NonLinearSolver/materialLaw/mlawHyperelastic.cpp b/NonLinearSolver/materialLaw/mlawHyperelastic.cpp index df3c8ed49..a33770601 100644 --- a/NonLinearSolver/materialLaw/mlawHyperelastic.cpp +++ b/NonLinearSolver/materialLaw/mlawHyperelastic.cpp @@ -314,7 +314,7 @@ mlawHyperViscoElastic::mlawHyperViscoElastic(const mlawHyperViscoElastic& src): _I4(src._I4), _I(src._I), _Idev(src._Idev) { - _temFunc_mu = NULL; // shear modulus; +/* _temFunc_mu = NULL; // shear modulus; if (src._temFunc_mu!=NULL){ _temFunc_mu = src._temFunc_mu->clone(); } @@ -329,7 +329,7 @@ mlawHyperViscoElastic::mlawHyperViscoElastic(const mlawHyperViscoElastic& src): for(int i=0;i<_N;i++){ if (src._temFunc_Gi[i]!=NULL) _temFunc_Gi[i] = src._temFunc_Gi[i]->clone(); - } + }*/ }; @@ -839,8 +839,7 @@ void mlawPowerYieldHyper::tangent_full_perturbation( const STensor3 &P, const STensor3 &F, const IPHyperViscoElastoPlastic* q0, - IPHyperViscoElastoPlastic* q1, - const double Tc + IPHyperViscoElastoPlastic* q1 ) const{ static STensor43 tmpSTensor43; @@ -851,7 +850,7 @@ void mlawPowerYieldHyper::tangent_full_perturbation( for (int j=0; j<3; j++){ Fplus = F; Fplus(i,j)+=_perturbationfactor; - this->predictorCorrector(Fplus,q0,&q11,Pplus,false,tmpSTensor43,tmpSTensor43,tmpSTensor43,NULL, Tc); + this->predictorCorrector(Fplus,q0,&q11,Pplus,false,tmpSTensor43,tmpSTensor43,tmpSTensor43,NULL); q1->_DgammaDF(i,j) = (q11._epspbarre - q1->_epspbarre)/_perturbationfactor; q1->_DirreversibleEnergyDF(i,j) = (q11._irreversibleEnergy - q1->_irreversibleEnergy)/_perturbationfactor; for (int k=0; k<3; k++){ @@ -882,23 +881,18 @@ void mlawPowerYieldHyper::constitutive( const bool dTangent, STensor63* dCalgdeps) const{ -//FLE - // double Tc =_Tref; - double Tc = 298.; -//FLE - const IPHyperViscoElastoPlastic *q0=dynamic_cast<const IPHyperViscoElastoPlastic *>(q0i); IPHyperViscoElastoPlastic *q1 = dynamic_cast<IPHyperViscoElastoPlastic *>(q1i); static STensor43 dFedF, dFpdF; if (_tangentByPerturbation){ - this->predictorCorrector(Fn,q0,q1,P,false,Tangent,dFedF,dFpdF,elasticTangent,Tc); + this->predictorCorrector(Fn,q0,q1,P,false,Tangent,dFedF,dFpdF,elasticTangent); if (stiff) - this->tangent_full_perturbation(Tangent,dFedF,dFpdF,P,Fn,q0,q1,Tc); + this->tangent_full_perturbation(Tangent,dFedF,dFpdF,P,Fn,q0,q1); } else{ - this->predictorCorrector(Fn,q0,q1,P,stiff,Tangent,dFedF,dFpdF,elasticTangent,Tc); + this->predictorCorrector(Fn,q0,q1,P,stiff,Tangent,dFedF,dFpdF,elasticTangent); } }; @@ -950,7 +944,7 @@ void mlawPowerYieldHyper::getYieldCoefficientDerivatives(const IPHyperViscoElast }; void mlawPowerYieldHyper::predictorCorrector_nonAssociatedFlow(const STensor3& F, const IPHyperViscoElastoPlastic *q0, IPHyperViscoElastoPlastic *q1, - STensor3&P, const bool stiff, STensor43& Tangent, STensor43& dFedF, STensor43& dFpdF, const double Tc) const{ + STensor3&P, const bool stiff, STensor43& Tangent, STensor43& dFedF, STensor43& dFpdF) const{ /* compute elastic predictor */ STensor3& Fp1 = q1->_Fp; const STensor3& Fp0 = q0->_Fp; @@ -985,7 +979,7 @@ void mlawPowerYieldHyper::predictorCorrector_nonAssociatedFlow(const STensor3& F // update A, B double Ge, Ke; - viscoElasticPredictor(Ee,q0->_Ee,q0,q1,Ke,Ge,Tc); // This should update the moduli values to the current temperature right? + viscoElasticPredictor(Ee,q0->_Ee,q0,q1,Ke,Ge); // This should update the moduli values to the current temperature right? static STensor3 PhiPr; PhiPr = q1->_kirchhoff; @@ -1149,7 +1143,7 @@ void mlawPowerYieldHyper::predictorCorrector_nonAssociatedFlow(const STensor3& F STensorOperation::logSTensor3(Ce,_order,Ee,&DlnDCe,&DDlnDDCe); Ee *= 0.5; // update A, B - updateViscoElasticFlow(q0,q1,Ke,Ge,Tc); // This wont cause a problem because the moduli are reinitiated inside the function. + updateViscoElasticFlow(q0,q1,Ke,Ge); // This wont cause a problem because the moduli are reinitiated inside the function. // backstress static STensor3 DB; // increment @@ -1503,7 +1497,7 @@ void mlawPowerYieldHyper::dPlasticEnergydF(const STensor3 &N, double Gamma, doub } void mlawPowerYieldHyper::predictorCorrector_associatedFlow(const STensor3& F, const IPHyperViscoElastoPlastic *q0, IPHyperViscoElastoPlastic *q1, - STensor3&P, const bool stiff, STensor43& Tangent, STensor43& dFedF, STensor43& dFpdF, const double Tc) const{ + STensor3&P, const bool stiff, STensor43& Tangent, STensor43& dFedF, STensor43& dFpdF) const{ /* compute elastic predictor */ STensor3& Fp1 = q1->_Fp; const STensor3& Fp0 = q0->_Fp; @@ -1537,7 +1531,7 @@ void mlawPowerYieldHyper::predictorCorrector_associatedFlow(const STensor3& F, c // update A, B double Ge, Ke; - viscoElasticPredictor(Ee,q0->_Ee,q0,q1,Ke,Ge,Tc); //Moduli being updated here. + viscoElasticPredictor(Ee,q0->_Ee,q0,q1,Ke,Ge); //Moduli being updated here. static STensor3 devKpr; // dev corotational kirchoff stress predictor static double ppr; // pressure predictor @@ -1710,7 +1704,7 @@ void mlawPowerYieldHyper::predictorCorrector_associatedFlow(const STensor3& F, c STensorOperation::logSTensor3(Ce,_order,Ee,&DlnDCe,&DDlnDDCe); Ee *= 0.5; - updateViscoElasticFlow(q0,q1,Ke,Ge,Tc); // Moduli being updated here. + updateViscoElasticFlow(q0,q1,Ke,Ge); // Moduli being updated here. } else{ q1->getRefToDissipationActive() = false; @@ -2040,27 +2034,27 @@ void mlawPowerYieldHyper::predictorCorrector_associatedFlow(const STensor3& F, c void mlawPowerYieldHyper::predictorCorrector(const STensor3& F, const IPHyperViscoElastoPlastic *q0, IPHyperViscoElastoPlastic *q1, STensor3&P, const bool stiff, STensor43& Tangent,STensor43& dFedF, STensor43& dFpdF, - STensor43* elasticTangent, const double Tc) const{ + STensor43* elasticTangent) const{ if (_tangentByPerturbation){ if (_nonAssociatedFlow){ - this->predictorCorrector_nonAssociatedFlow(F,q0,q1,P,false,Tangent,dFedF,dFpdF,Tc); + this->predictorCorrector_nonAssociatedFlow(F,q0,q1,P,false,Tangent,dFedF,dFpdF); } else{ - this->predictorCorrector_associatedFlow(F,q0,q1,P,false,Tangent,dFedF,dFpdF,Tc); + this->predictorCorrector_associatedFlow(F,q0,q1,P,false,Tangent,dFedF,dFpdF); } if (stiff){ - tangent_full_perturbation(Tangent,dFedF,dFpdF,P,F,q0,q1,Tc); + tangent_full_perturbation(Tangent,dFedF,dFpdF,P,F,q0,q1); } } else{ if (_nonAssociatedFlow){ - this->predictorCorrector_nonAssociatedFlow(F,q0,q1,P,stiff,Tangent,dFedF,dFpdF,Tc); + this->predictorCorrector_nonAssociatedFlow(F,q0,q1,P,stiff,Tangent,dFedF,dFpdF); } else{ - this->predictorCorrector_associatedFlow(F,q0,q1,P,stiff,Tangent,dFedF,dFpdF,Tc); + this->predictorCorrector_associatedFlow(F,q0,q1,P,stiff,Tangent,dFedF,dFpdF); } } // compute mechanical tengent @@ -2070,7 +2064,7 @@ void mlawPowerYieldHyper::predictorCorrector(const STensor3& F, const IPHyperVis static IPHyperViscoElastoPlastic q1tmp(*q1); q1tmp.operator=(*(dynamic_cast<const IPVariable*> (q1))); mlawHyperViscoElastic::predictorCorrector_ViscoElastic(q0->getConstRefToFe(), q1->getConstRefToFe(), Ptmp, q0, &q1tmp, - true, *elasticTangent, Tc); + true, *elasticTangent); } }; @@ -2193,8 +2187,8 @@ void mlawPowerYieldHyperWithFailure::checkInternalState(IPVariable* ipv, const I void mlawPowerYieldHyperWithFailure::predictorCorrector(const STensor3& F, const IPHyperViscoElastoPlastic *q0, IPHyperViscoElastoPlastic *q1, STensor3&P, const bool stiff, STensor43& Tangent, STensor43& dFedF, STensor43& dFpdF, - STensor43* elasticTangent, const double Tc) const{ - mlawPowerYieldHyper::predictorCorrector(F,q0,q1,P,stiff,Tangent,dFedF,dFpdF, elasticTangent,Tc); + STensor43* elasticTangent) const{ + mlawPowerYieldHyper::predictorCorrector(F,q0,q1,P,stiff,Tangent,dFedF,dFpdF, elasticTangent); // compute failure r and DrDF double& gF = q1->getRefToFailurePlasticity(); diff --git a/NonLinearSolver/materialLaw/mlawHyperelastic.h b/NonLinearSolver/materialLaw/mlawHyperelastic.h index d2b482280..3876ee4c2 100644 --- a/NonLinearSolver/materialLaw/mlawHyperelastic.h +++ b/NonLinearSolver/materialLaw/mlawHyperelastic.h @@ -279,16 +279,16 @@ class mlawPowerYieldHyper : public mlawHyperViscoElastic{ const STensor43 &DdevNDCepr, const STensor3 &DtrNDCepr, const STensor3 &KS, const STensor43 &CeprToF, STensor3 &dPlasticEnergydF) const; virtual void predictorCorrector_nonAssociatedFlow(const STensor3& F, const IPHyperViscoElastoPlastic *q0_, IPHyperViscoElastoPlastic *q1, - STensor3&P, const bool stiff, STensor43& Tangent, STensor43& dFedF, STensor43& dFpdF, const double Tc) const; + STensor3&P, const bool stiff, STensor43& Tangent, STensor43& dFedF, STensor43& dFpdF) const; virtual void predictorCorrector_associatedFlow(const STensor3& F, const IPHyperViscoElastoPlastic *q0_, IPHyperViscoElastoPlastic *q1, - STensor3&P, const bool stiff, STensor43& Tangent, STensor43& dFedF, STensor43& dFpdF, const double Tc) const; + STensor3&P, const bool stiff, STensor43& Tangent, STensor43& dFedF, STensor43& dFpdF) const; protected: virtual void predictorCorrector(const STensor3& F, const IPHyperViscoElastoPlastic *q0_, IPHyperViscoElastoPlastic *q1, STensor3&P, const bool stiff, STensor43& Tangent, STensor43& dFedF, STensor43& dFpdF, - STensor43* elasticTangent, const double Tc) const; + STensor43* elasticTangent) const; virtual void tangent_full_perturbation( // compute tangent by perturbation STensor43 &T, // current tangent STensor43& dFedF, @@ -296,8 +296,7 @@ class mlawPowerYieldHyper : public mlawHyperViscoElastic{ const STensor3 &P, // current PK stress const STensor3 &F, // current deformation gradient const IPHyperViscoElastoPlastic* q0, // previous internal variables - IPHyperViscoElastoPlastic* q1, // current internal variables - const double Tc + IPHyperViscoElastoPlastic* q1 ) const; #endif // SWIG @@ -355,7 +354,7 @@ class mlawPowerYieldHyperWithFailure : public mlawPowerYieldHyper{ protected: virtual void predictorCorrector(const STensor3& F, const IPHyperViscoElastoPlastic *q0_, IPHyperViscoElastoPlastic *q1, STensor3&P, const bool stiff, STensor43& Tangent, STensor43& dFedF, STensor43& dFpdF, - STensor43* elasticTangent, const double Tc) const; + STensor43* elasticTangent) const; #endif // SWIG public: mlawPowerYieldHyperWithFailure(const int num,const double E,const double nu, const double rho, diff --git a/NonLinearSolver/materialLaw/mlawNonLocalDamageHyperelastic.cpp b/NonLinearSolver/materialLaw/mlawNonLocalDamageHyperelastic.cpp index 9a94d0a96..0ca5f6a90 100644 --- a/NonLinearSolver/materialLaw/mlawNonLocalDamageHyperelastic.cpp +++ b/NonLinearSolver/materialLaw/mlawNonLocalDamageHyperelastic.cpp @@ -135,7 +135,7 @@ void mlawNonLocalDamagePowerYieldHyper::constitutive( static STensor43 dFedF, dFpdF; static STensor3 Peff; - mlawPowerYieldHyper::predictorCorrector(Fn,ipvprev,ipvcur,Peff,stiff,Tangent,dFedF,dFpdF,elasticTangent,Tc); + mlawPowerYieldHyper::predictorCorrector(Fn,ipvprev,ipvcur,Peff,stiff,Tangent,dFedF,dFpdF,elasticTangent); const STensor3& Fe = ipvcur->_Fe; const STensor3& dgammadF = ipvcur->_DgammaDF; @@ -324,7 +324,7 @@ void mlawLocalDamagePowerYieldHyperWithFailure::constitutive( double Tc = _Tref; static STensor43 dFedF, dFpdF; static STensor3 Peff; - mlawPowerYieldHyperWithFailure::predictorCorrector(Fn,ipvprev,ipvcur,Peff,stiff,Tangent,dFedF,dFpdF,NULL,Tc); + mlawPowerYieldHyperWithFailure::predictorCorrector(Fn,ipvprev,ipvcur,Peff,stiff,Tangent,dFedF,dFpdF,NULL); // get result from effective law const STensor3& Fe = ipvcur->_Fe; @@ -621,7 +621,7 @@ void mlawNonLocalDamagePowerYieldHyperWithFailure::predictorCorrector( // FLE double Tc = _Tref; - mlawPowerYieldHyperWithFailure::predictorCorrector(Fn,ipvprev,ipvcur,Peff,stiff,Tangent,dFedF,dFpdF, elasticTangent,Tc); + mlawPowerYieldHyperWithFailure::predictorCorrector(Fn,ipvprev,ipvcur,Peff,stiff,Tangent,dFedF,dFpdF, elasticTangent); // get result from effective law const STensor3& Fe = ipvcur->_Fe; diff --git a/dG3D/src/dG3DIPVariable.cpp b/dG3D/src/dG3DIPVariable.cpp index 9a65fbb63..a31a273b0 100644 --- a/dG3D/src/dG3DIPVariable.cpp +++ b/dG3D/src/dG3DIPVariable.cpp @@ -4133,8 +4133,8 @@ NonLinearTVMDG3DIPVariable::NonLinearTVMDG3DIPVariable(const mlawNonLinearTVM& v _ipViscoElastic = new IPNonLinearTVM(viscoLaw.getViscoElasticNumberOfElement()); }; */ -NonLinearTVMDG3DIPVariable::NonLinearTVMDG3DIPVariable(const mlawNonLinearTVM& viscoLaw, double Tinitial, const bool oninter) : - ThermoMechanicsDG3DIPVariableBase(oninter){ +NonLinearTVMDG3DIPVariable::NonLinearTVMDG3DIPVariable(const mlawNonLinearTVM& viscoLaw, double Tinitial, const bool createBodyForceHO, const bool oninter) : + ThermoMechanicsDG3DIPVariableBase(createBodyForceHO, oninter){ _ipViscoElastic = new IPNonLinearTVM(viscoLaw.getViscoElasticNumberOfElement()); this->getRefToTemperature()=Tinitial; }; diff --git a/dG3D/src/dG3DIPVariable.h b/dG3D/src/dG3DIPVariable.h index 65f58ade1..58274723b 100644 --- a/dG3D/src/dG3DIPVariable.h +++ b/dG3D/src/dG3DIPVariable.h @@ -4753,13 +4753,12 @@ class NonLinearTVMDG3DIPVariable : public ThermoMechanicsDG3DIPVariableBase{ / protected: IPNonLinearTVM * _ipViscoElastic; - NonLinearTVMDG3DIPVariable(const bool oninter=false) {Msg::Error("NonLinearTVMDG3DIPVariable: Cannot use the default constructor of ThermoMechanicsDG3DIPVariableBase");} + NonLinearTVMDG3DIPVariable(const bool createBodyForceHO, const bool oninter):ThermoMechanicsDG3DIPVariableBase(createBodyForceHO, oninter) {Msg::Error("NonLinearTVMDG3DIPVariable: Cannot use the default constructor of ThermoMechanicsDG3DIPVariableBase");} public: // Constructor Changed // NonLinearTVMDG3DIPVariable(const mlawNonLinearTVM& viscoLaw, const bool oninter=false, int nbExtraDof=0); // ?? - WHAT to add/change? - NonLinearTVMDG3DIPVariable(const mlawNonLinearTVM& viscoLaw, double Tinitial, const bool oninter=false); // Changed - + NonLinearTVMDG3DIPVariable(const mlawNonLinearTVM& viscoLaw, double Tinitial, const bool createBodyForceHO, const bool oninter); // Changed NonLinearTVMDG3DIPVariable(const NonLinearTVMDG3DIPVariable& src); virtual NonLinearTVMDG3DIPVariable& operator = (const IPVariable& src); virtual ~NonLinearTVMDG3DIPVariable(); diff --git a/dG3D/src/dG3DMaterialLaw.cpp b/dG3D/src/dG3DMaterialLaw.cpp index 4b1ae0d27..7cdf57ca0 100644 --- a/dG3D/src/dG3DMaterialLaw.cpp +++ b/dG3D/src/dG3DMaterialLaw.cpp @@ -7423,20 +7423,20 @@ void NonLinearTVMDG3DMaterialLaw::setTemperatureFunction_BranchThermalDilationCo }; // Added extra argument in the NonLinearTVMDG3DIPVariable contructor - Tinitial (added from LinearThermoMech) -void NonLinearTVMDG3DMaterialLaw::createIPState(IPStateBase* &ips,const bool* state_,const MElement *ele, const int nbFF_, const IntPt *GP, const int gpt) const +void NonLinearTVMDG3DMaterialLaw::createIPState(IPStateBase* &ips, bool hasBodyForce, 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 NonLinearTVMDG3DIPVariable(_viscoLaw,_viscoLaw.getInitialTemperature(),inter); // Why 3 times? --- WHAT?? Add getInitialTemperature in Material Law - IPVariable* ipv1 = new NonLinearTVMDG3DIPVariable(_viscoLaw,_viscoLaw.getInitialTemperature(),inter); - IPVariable* ipv2 = new NonLinearTVMDG3DIPVariable(_viscoLaw,_viscoLaw.getInitialTemperature(),inter); + IPVariable* ipvi = new NonLinearTVMDG3DIPVariable(_viscoLaw,_viscoLaw.getInitialTemperature(),hasBodyForce,inter); // Why 3 times? --- WHAT?? Add getInitialTemperature in Material Law + IPVariable* ipv1 = new NonLinearTVMDG3DIPVariable(_viscoLaw,_viscoLaw.getInitialTemperature(),hasBodyForce,inter); + IPVariable* ipv2 = new NonLinearTVMDG3DIPVariable(_viscoLaw,_viscoLaw.getInitialTemperature(),hasBodyForce,inter); if(ips != NULL) delete ips; ips = new IP3State(state_,ipvi,ipv1,ipv2); } -void NonLinearTVMDG3DMaterialLaw::createIPVariable(IPVariable* &ipv,const MElement *ele, const int nbFF_, const IntPt *GP, const int gpt) const +void NonLinearTVMDG3DMaterialLaw::createIPVariable(IPVariable* &ipv, bool hasBodyForce, const MElement *ele, const int nbFF_, const IntPt *GP, const int gpt) const { if(ipv !=NULL) delete ipv; bool inter=true; @@ -7444,7 +7444,7 @@ void NonLinearTVMDG3DMaterialLaw::createIPVariable(IPVariable* &ipv,const MEleme if(iele == NULL){ inter=false; } - ipv = new NonLinearTVMDG3DIPVariable(_viscoLaw,_viscoLaw.getInitialTemperature(),inter); + ipv = new NonLinearTVMDG3DIPVariable(_viscoLaw,_viscoLaw.getInitialTemperature(),hasBodyForce,inter); } // 2 Member Functions added from LinearThermoMechanicsDG3DMaterialLaw -- Add these functions in the material law diff --git a/dG3D/src/dG3DMaterialLaw.h b/dG3D/src/dG3DMaterialLaw.h index 0f2bfd394..f62764b42 100644 --- a/dG3D/src/dG3DMaterialLaw.h +++ b/dG3D/src/dG3DMaterialLaw.h @@ -2580,8 +2580,8 @@ class NonLinearTVMDG3DMaterialLaw : public dG3DMaterialLaw{ // public Material _viscoLaw.setTime(t,dtime); } virtual materialLaw::matname getType() const {return _viscoLaw.getType();} - virtual void createIPState(IPStateBase* &ips,const bool* state_=NULL,const MElement *ele=NULL, const int nbFF_=0, const IntPt *GP=NULL, const int gpt = 0) const; - virtual void createIPVariable(IPVariable* &ipv,const MElement *ele, const int nbFF_, const IntPt *GP, const int gpt) const; + virtual void createIPState(IPStateBase* &ips, bool hasBodyForce, const bool* state_=NULL,const MElement *ele=NULL, const int nbFF_=0, const IntPt *GP=NULL, const int gpt = 0) const; + virtual void createIPVariable(IPVariable* &ipv, bool hasBodyForce, const MElement *ele, const int nbFF_, const IntPt *GP, const int gpt) const; virtual void initLaws(const std::map<int,materialLaw*> &maplaw){} virtual double soundSpeed() const{return _viscoLaw.soundSpeed();} // or change value ?? virtual void stress(IPVariable*ipv, const IPVariable*ipvprev, const bool stiff=true, const bool checkfrac=true, const bool dTangent=false); -- GitLab