diff --git a/dG3D/src/dG3DMaterialLaw.cpp b/dG3D/src/dG3DMaterialLaw.cpp index 677bdbbe0b68ef4a4f330fd45812fa39a81852ec..4225095d4ad640f04b29b75fd1ddade5e2872a68 100644 --- a/dG3D/src/dG3DMaterialLaw.cpp +++ b/dG3D/src/dG3DMaterialLaw.cpp @@ -1910,7 +1910,7 @@ torchANNBasedDG3DMaterialLaw::torchANNBasedDG3DMaterialLaw(const int num, const _EXXmean(EXXmean), _EXXstd(EXXstd), _EXYmean(EXYmean), _EXYstd(EXYstd), _EYYmean(EYYmean), _EYYstd(EYYstd), _EYZmean(EYZmean), _EYZstd(EYZstd), _EZZmean(EZZmean), _EZZstd(EZZstd), _EZXmean(EZXmean), _EZXstd(EZXstd), _SXXmean(SXXmean), _SXXstd(SXXstd), _SXYmean(SXYmean), _SXYstd(SXYstd), _SYYmean(SYYmean), _SYYstd(SYYstd), _SYZmean(SYZmean), _SYZstd(SYZstd), _SZZmean(SZZmean), _SZZstd(SZZstd), - _SZXmean(SZXmean), _SZXstd(SZXstd), _tangentByPerturbation(pert), _pertTol(tol), _kinematicInput(EGL), _NeedExtraNorm(false), _DoubleInput(false), + _SZXmean(SZXmean), _SZXstd(SZXstd), _tangentByPerturbation(pert), _pertTol(tol), _kinematicInput(EGL), _NeedExtraNorm(false), _numberOfExtraInput(0) { #if defined(HAVE_TORCH) @@ -1921,6 +1921,7 @@ torchANNBasedDG3DMaterialLaw::torchANNBasedDG3DMaterialLaw(const int num, const Msg::Error("error loading the model"); } module.eval(); + if(_tangentByPerturbation){torch::NoGradGuard no_grad;} if(_numberOfInput == 3){ // for 2D case VXX = torch::zeros({1, 1, 4}); @@ -1978,7 +1979,7 @@ torchANNBasedDG3DMaterialLaw::torchANNBasedDG3DMaterialLaw(const torchANNBasedDG _SXXmean(src._SXXmean), _SXXstd(src._SXXstd), _SXYmean(src._SXYmean), _SXYstd(src._SXYstd), _SYYmean(src._SYYmean), _SYYstd(src._SYYstd), _SYZmean(src._SYZmean), _SYZstd(src._SYZstd), _SZZmean(src._SZZmean), _SZZstd(src._SZZstd), _SZXmean(src._SZXmean), _SZXstd(src._SZXstd), _tangentByPerturbation(src._tangentByPerturbation), _pertTol(src._pertTol), _kinematicInput(src._kinematicInput), _NeedExtraNorm(src._NeedExtraNorm), - _DoubleInput(src._DoubleInput), _numberOfExtraInput(src._numberOfExtraInput), _initialValue_ExtraInput(src._initialValue_ExtraInput), + _numberOfExtraInput(src._numberOfExtraInput), _initialValue_ExtraInput(src._initialValue_ExtraInput), _normParams_ExtraInput(src._normParams_ExtraInput) { #if defined(HAVE_TORCH) @@ -2098,7 +2099,7 @@ void torchANNBasedDG3DMaterialLaw::initLaws(const std::map<int,materialLaw*> &ma // Added additional container Extra1 for additonal inputs @Mohib TODO: Remove hardcode fullMatrix<double> E1(1,6), S(1,6), Extra1(1, _numberOfExtraInput); torch::Tensor h_init = _initialValue_h*torch::ones({1, 1, _numberOfInternalVariables}); - static torch::Tensor h_tmp = torch::zeros({1, 1, _numberOfInternalVariables}); + auto h_tmp = torch::zeros({1, 1, _numberOfInternalVariables}, torch::requires_grad(false)); fullMatrix<double> DSDE(6,6); //Set initial values for extra inputs @Mohib @@ -2334,7 +2335,7 @@ void torchANNBasedDG3DMaterialLaw::stress(IPVariable* ipv, const IPVariable* ipv { RNNstress_stiff(E0, E1, h0, h1, S, false, DSDE); static fullMatrix<double> E1_plus(1,6),S_plus(1,6); - static torch::Tensor h_tmp = torch::zeros({1, 1, _numberOfInternalVariables}); + auto h_tmp = torch::zeros({1, 1, _numberOfInternalVariables}, torch::requires_grad(false)); for (int i=0; i< 6; i++) { @@ -2352,7 +2353,7 @@ void torchANNBasedDG3DMaterialLaw::stress(IPVariable* ipv, const IPVariable* ipv RNNstressGeo_stiff(E1, Extra1, h0, h1, S, false, DSDE); //S.print("original"); static fullMatrix<double> E1_plus(1,6), S_plus(1,6); - static torch::Tensor h_tmp = torch::zeros({1, 1, _numberOfInternalVariables}); + torch::Tensor h_tmp = torch::zeros({1, 1, _numberOfInternalVariables}); for (int i=0; i< 6; i++) { @@ -2484,20 +2485,21 @@ void torchANNBasedDG3DMaterialLaw::stress(IPVariable* ipv, const IPVariable* ipv void torchANNBasedDG3DMaterialLaw::RNNstress_stiff(const fullMatrix<double>& E0, const fullMatrix<double>& E1, const torch::Tensor& h0, torch::Tensor& h1, fullMatrix<double>& S, const bool stiff, fullMatrix<double>& DSDE) { - static vector<float> E_vec(_numberOfInput), E0_vec(_numberOfInput); - static torch::Tensor E_norm, E0_norm, dE_norm; + vector<float> E_vec(_numberOfInput), E0_vec(_numberOfInput); + torch::Tensor E_norm, E0_norm, dE_norm; + if(!stiff){torch::NoGradGuard no_grad;} // use RNN to predict stress------------- vector<torch::jit::IValue> inputs; Normalize_strain(E1, E_vec); if(_NeedExtraNorm){ Normalize_strain(E0, E0_vec); } - if(_DoubleInput and _NeedExtraNorm){ + if(_NeedExtraNorm){ if(stiff){ E0_norm = torch::from_blob(E0_vec.data(), {1,1, _numberOfInput}, torch::requires_grad()); } else{ - E0_norm = torch::from_blob(E0_vec.data(), {1,1, _numberOfInput}); + E0_norm = torch::from_blob(E0_vec.data(), {1,1, _numberOfInput}, torch::requires_grad(false)); } inputs.push_back(E0_norm); } @@ -2506,29 +2508,27 @@ void torchANNBasedDG3DMaterialLaw::RNNstress_stiff(const fullMatrix<double>& E0, E_norm = torch::from_blob(E_vec.data(), {1,1, _numberOfInput}, torch::requires_grad()); } else{ - E_norm = torch::from_blob(E_vec.data(), {1,1, _numberOfInput}); + E_norm = torch::from_blob(E_vec.data(), {1,1, _numberOfInput}, torch::requires_grad(false)); } inputs.push_back(E_norm); } - static vector<float> norm_E1(1); - static vector<float> DE_vec(_numberOfInput); - static torch::Tensor norm; + vector<float> norm_E1(1); + vector<float> DE_vec(_numberOfInput); + torch::Tensor norm; if(_kinematicInput == torchANNBasedDG3DMaterialLaw::EGLInc or _NeedExtraNorm){ if(_kinematicInput == torchANNBasedDG3DMaterialLaw::EGL and _NeedExtraNorm){ for (int i = 0; i <_numberOfInput; i++) { DE_vec[i] = E_vec[i]-E0_vec[i]; } - if(_DoubleInput){ - if(stiff){ - dE_norm = torch::from_blob(DE_vec.data(), {1,1, _numberOfInput}, torch::requires_grad()); - } - else{ - dE_norm = torch::from_blob(DE_vec.data(), {1,1, _numberOfInput}); - } - inputs.push_back(dE_norm); - } + if(stiff){ + dE_norm = torch::from_blob(DE_vec.data(), {1,1, _numberOfInput}, torch::requires_grad()); + } + else{ + dE_norm = torch::from_blob(DE_vec.data(), {1,1, _numberOfInput}, torch::requires_grad(false)); + } + inputs.push_back(dE_norm); norm_E1[0] = sqrt(std::inner_product(DE_vec.begin(), DE_vec.end(), DE_vec.begin(), 0.0)); } else{ @@ -2539,7 +2539,7 @@ void torchANNBasedDG3DMaterialLaw::RNNstress_stiff(const fullMatrix<double>& E0, norm = torch::from_blob(norm_E1.data(), {1,1}, torch::requires_grad()); } else{ - norm = torch::from_blob(norm_E1.data(), {1,1}); + norm = torch::from_blob(norm_E1.data(), {1,1}, torch::requires_grad(false)); } inputs.push_back(norm); } diff --git a/dG3D/src/dG3DMaterialLaw.h b/dG3D/src/dG3DMaterialLaw.h index 33162faff50f1c30faf3678ede3de41f5372c87b..d7a115343c9e977626761559d7336e260df6f35c 100644 --- a/dG3D/src/dG3DMaterialLaw.h +++ b/dG3D/src/dG3DMaterialLaw.h @@ -476,7 +476,6 @@ class torchANNBasedDG3DMaterialLaw : public dG3DMaterialLaw{ int _numberOfInternalVariables; double _initialValue_h; bool _NeedExtraNorm; - bool _DoubleInput; int _numberOfExtraInput; // number of additional inputs other than 6 (or 3 for 2D) strain components @Mohib std::vector<double> _initialValue_ExtraInput; // initial values for the extra inputs @Mohib std::vector<double> _normParams_ExtraInput; // Container for additional inputs, it is needed for normalization @Mohib @@ -487,7 +486,6 @@ class torchANNBasedDG3DMaterialLaw : public dG3DMaterialLaw{ // TODO: Read from ipvariable. // double _latticeRadius; - // double _latticeSize; // double _pTime; // double _cTime; @@ -549,10 +547,6 @@ class torchANNBasedDG3DMaterialLaw : public dG3DMaterialLaw{ _NeedExtraNorm = NeedExtraNorm; if(_NeedExtraNorm) Msg::Info("Extra Norm is applied"); }; - void setDoubleInput (const bool DoubleInput){ - _DoubleInput = DoubleInput; - if(_DoubleInput) Msg::Info("Double Input is applied"); - }; void setInitialHValue(double initialValue_h) {_initialValue_h=initialValue_h;}; diff --git a/dG3D/src/nonLocalDamageDG3DIPVariable.cpp b/dG3D/src/nonLocalDamageDG3DIPVariable.cpp index be80f5a1ebd7ae4aba7faead820e6a37f65f6539..84e65dc7c5104fa3cdb8085f3a55e0baf21ce1b0 100644 --- a/dG3D/src/nonLocalDamageDG3DIPVariable.cpp +++ b/dG3D/src/nonLocalDamageDG3DIPVariable.cpp @@ -280,18 +280,29 @@ void nonLocalDamageDG3DIPVariable::restart() /* Beginning IPvariable for nonlocalDamageTorchANNDG3DIPVariable*/ -nonLocalDamageTorchANNDG3DIPVariableDG3DIPVariable::nonLocalDamageTorchANNDG3DIPVariableDG3DIPVariable(const int n, const double initial_h, const int numNonlocal, const std::vector<STensor3>& length, const bool createBodyForceHO, const bool oninter): - dG3DIPVariable(createBodyForceHO, oninter,numNonlocal), _n(n), _initial_h(initial_h), characteristicLength(length), _nlD(1,6), _D(1,6), +nonLocalDamageTorchANNDG3DIPVariableDG3DIPVariable::nonLocalDamageTorchANNDG3DIPVariableDG3DIPVariable(const CLengthLaw *cll, const int n, const double initial_h, const int numNonlocal, const bool createBodyForceHO, const bool oninter): + dG3DIPVariable(createBodyForceHO, oninter,numNonlocal), _n(n), _initial_h(initial_h), _nlD(1,6), _D(1,6), restart_internalVars(1,n), _kinematicVariables(1,6) { localVar.resize(numNonlocal,0.); #if defined(HAVE_TORCH) _internalVars = _initial_h*torch::ones({1, 1, _n}); #endif + ipvCL=NULL; + if(cll ==NULL) Msg::Error("nonLocalDamageTorchANNDG3DIPVariableDG3DIPVariable::nonLocalDamageTorchANNDG3DIPVariableDG3DIPVariable has no cll"); + cll->createIPVariable(ipvCL); + } nonLocalDamageTorchANNDG3DIPVariableDG3DIPVariable::nonLocalDamageTorchANNDG3DIPVariableDG3DIPVariable(const nonLocalDamageTorchANNDG3DIPVariableDG3DIPVariable &source):dG3DIPVariable(source),_n(source._n), _initial_h(source._initial_h), _nlD(source._nlD), - _D(source._D), localVar(source.localVar), characteristicLength(source.characteristicLength), _internalVars(source._internalVars), - _kinematicVariables(source._kinematicVariables){}; + _D(source._D), localVar(source.localVar), _internalVars(source._internalVars), + _kinematicVariables(source._kinematicVariables) +{ + ipvCL = NULL; + if(source.ipvCL != NULL) + { + ipvCL = dynamic_cast<IPCLength*>(source.ipvCL->clone()); + } +}; nonLocalDamageTorchANNDG3DIPVariableDG3DIPVariable& nonLocalDamageTorchANNDG3DIPVariableDG3DIPVariable::operator =(const IPVariable& src) { @@ -308,13 +319,30 @@ nonLocalDamageTorchANNDG3DIPVariableDG3DIPVariable& nonLocalDamageTorchANNDG3DIP _nlD = psrc->_nlD; _D = psrc->_D; localVar = psrc->localVar; + if(psrc->ipvCL != NULL) + { + if (ipvCL != NULL){ + ipvCL->operator=(*dynamic_cast<const IPCLength*>(psrc->ipvCL)); + } + else + ipvCL= dynamic_cast<IPCLength*>(psrc->ipvCL->clone()); + } + else{ + if(ipvCL != NULL) delete ipvCL; ipvCL = NULL; + } + } #else Msg::Error("NOT COMPILED WITH TORCH"); #endif - return *this; } -nonLocalDamageTorchANNDG3DIPVariableDG3DIPVariable::~nonLocalDamageTorchANNDG3DIPVariableDG3DIPVariable(){}; +nonLocalDamageTorchANNDG3DIPVariableDG3DIPVariable::~nonLocalDamageTorchANNDG3DIPVariableDG3DIPVariable() +{ + if(ipvCL != NULL) + { + delete ipvCL; + } +}; void nonLocalDamageTorchANNDG3DIPVariableDG3DIPVariable::restart() @@ -329,6 +357,7 @@ void nonLocalDamageTorchANNDG3DIPVariableDG3DIPVariable::restart() restartManager::restart(_nlD.getDataPtr(),_nlD.size1()*_nlD.size2()); restartManager::restart(_D.getDataPtr(),_nlD.size1()*_nlD.size2()); restartManager::restart(localVar); + restartManager::restart(ipvCL); #else Msg::Error("NOT COMPILED WITH TORCH"); diff --git a/dG3D/src/nonLocalDamageDG3DIPVariable.h b/dG3D/src/nonLocalDamageDG3DIPVariable.h index d25a56a78d0db9bf5182d5b2e146d1ff4eeb0729..7ab92293562454b017f95a76472be580b60825f2 100644 --- a/dG3D/src/nonLocalDamageDG3DIPVariable.h +++ b/dG3D/src/nonLocalDamageDG3DIPVariable.h @@ -1496,7 +1496,7 @@ class nonLocalDamageTorchANNDG3DIPVariableDG3DIPVariable : public dG3DIPVariable fullMatrix<double> _nlD; fullMatrix<double> _D; std::vector<double> localVar; - const std::vector<STensor3>& characteristicLength; + IPCLength* ipvCL; #if defined(HAVE_TORCH) torch::Tensor _internalVars; // internal variable @@ -1507,7 +1507,7 @@ class nonLocalDamageTorchANNDG3DIPVariableDG3DIPVariable : public dG3DIPVariable fullMatrix<double> restart_internalVars; // internal variable public: - nonLocalDamageTorchANNDG3DIPVariableDG3DIPVariable(const int n, const double initial_h, const int numNonlocal, const std::vector<STensor3>& length, const bool createBodyForceHO, const bool oninter); + nonLocalDamageTorchANNDG3DIPVariableDG3DIPVariable(const CLengthLaw *cll, const int n, const double initial_h, const int numNonlocal, const bool createBodyForceHO, const bool oninter); nonLocalDamageTorchANNDG3DIPVariableDG3DIPVariable(const nonLocalDamageTorchANNDG3DIPVariableDG3DIPVariable &source); virtual nonLocalDamageTorchANNDG3DIPVariableDG3DIPVariable& operator =(const IPVariable& src); virtual ~nonLocalDamageTorchANNDG3DIPVariableDG3DIPVariable(); @@ -1533,14 +1533,43 @@ class nonLocalDamageTorchANNDG3DIPVariableDG3DIPVariable : public dG3DIPVariable virtual double getConstRefToNonLocalVariable(const int idex) const {return _nonlocalData->nonLocalVariable(idex);}; virtual double& getRefToNonLocalVariable(const int idex) {return _nonlocalData->nonLocalVariable(idex);}; + + virtual double getConstRefToNonLocalVariableInMaterialLaw(const int idex) const {return _nonlocalData->nonLocalVariable(idex);}; + virtual double& getRefToNonLocalVariableInMaterialLaw(const int idex) {return _nonlocalData->nonLocalVariable(idex);}; + + virtual const fullMatrix<double> &getConstRefToNonLocalDamage() const {return _nlD;}; virtual fullMatrix<double> &getRefToNonLocalDamage() {return _nlD;}; virtual const fullMatrix<double> &getConstRefToLocalDamage() const {return _D;}; virtual fullMatrix<double> &getRefToLocalDamage() {return _D;}; - virtual const STensor3 &getConstRefToCharacteristicLengthMatrix(const int idex) const {return characteristicLength[idex];}; - + virtual const IPCLength &getConstRefToIPCLength() const + { + if(ipvCL==NULL) + Msg::Error("nonLocalDamageTorchANNDG3DIPVariableDG3DIPVariable: ipvCL not initialized"); + return *ipvCL; + } + virtual IPCLength &getRefToIPCLength() + { + if(ipvCL==NULL) + Msg::Error("nonLocalDamageTorchANNDG3DIPVariableDG3DIPVariable: ipvCL not initialized"); + return *ipvCL; + } + virtual const STensor3 &getConstRefToCharacteristicLengthMatrix(const int idex) const + { + if(ipvCL==NULL) + Msg::Error("nonLocalDamageTorchANNDG3DIPVariableDG3DIPVariable: ipvCL not initialized"); + + if (idex == 0) + return ipvCL->getConstRefToCL(); + else + { + Msg::Error("the non-local variable %d is not defined",idex); + return ipvCL->getConstRefToCL(); + } + } + virtual IPVariable* clone() const {return new nonLocalDamageTorchANNDG3DIPVariableDG3DIPVariable(*this);}; }; diff --git a/dG3D/src/nonLocalDamageDG3DMaterialLaw.cpp b/dG3D/src/nonLocalDamageDG3DMaterialLaw.cpp index 1dbfea01c67f9afcc03cd90fccd3caff3d9d2a95..066f6fe37e1230699c51e45efea15ab39d9d660d 100644 --- a/dG3D/src/nonLocalDamageDG3DMaterialLaw.cpp +++ b/dG3D/src/nonLocalDamageDG3DMaterialLaw.cpp @@ -2219,7 +2219,7 @@ NonlocalDamageTorchANNBasedDG3DMaterialLaw::NonlocalDamageTorchANNBasedDG3DMater _SZXmean(SZXmean), _SZXstd(SZXstd), _EffSXXmean(EffSXXmean), _EffSXXstd(EffSXXstd), _EffSXYmean(EffSXYmean), _EffSXYstd(EffSXYstd), _EffSYYmean(EffSYYmean), _EffSYYstd(EffSYYstd), _EffSYZmean(EffSYZmean), _EffSYZstd(EffSYZstd), _EffSZZmean(EffSZZmean), _EffSZZstd(EffSZZstd), _EffSZXmean(EffSZXmean), _EffSZXstd(EffSZXstd), _localVarmean(localVarmean), - _localVarstd(localVarstd), _tangentByPerturbation(pert), _pertTol(tol), _kinematicInput(EGL) + _localVarstd(localVarstd), _tangentByPerturbation(pert), _pertTol(tol), _kinematicInput(EGL),_NeedExtraNorm(false),cLLaw(NULL) { #if defined(HAVE_TORCH) try{ @@ -2305,7 +2305,7 @@ NonlocalDamageTorchANNBasedDG3DMaterialLaw::NonlocalDamageTorchANNBasedDG3DMater _SYZstd(EffSYZstd), _SZZmean(EffSZZmean), _SZZstd(EffSZZstd), _SZXmean(EffSZXmean), _SZXstd(EffSZXstd), _EffSXXmean(EffSXXmean), _EffSXXstd(EffSXXstd), _EffSXYmean(EffSXYmean), _EffSXYstd(EffSXYstd), _EffSYYmean(EffSYYmean), _EffSYYstd(EffSYYstd), _EffSYZmean(EffSYZmean), _EffSYZstd(EffSYZstd), _EffSZZmean(EffSZZmean), _EffSZZstd(EffSZZstd), _EffSZXmean(EffSZXmean), _EffSZXstd(EffSZXstd), _localVarmean(localVarmean), - _localVarstd(localVarstd),_tangentByPerturbation(pert), _pertTol(tol), _kinematicInput(EGL) + _localVarstd(localVarstd),_tangentByPerturbation(pert), _pertTol(tol), _kinematicInput(EGL),_NeedExtraNorm(false),cLLaw(NULL) { #if defined(HAVE_TORCH) try{ @@ -2388,7 +2388,7 @@ NonlocalDamageTorchANNBasedDG3DMaterialLaw::NonlocalDamageTorchANNBasedDG3DMater _EffSXXmean(src._EffSXXmean), _EffSXXstd(src._EffSXXstd), _EffSXYmean(src._EffSXYmean), _EffSXYstd(src._EffSXYstd), _EffSYYmean(src._EffSYYmean), _EffSYYstd(src._EffSYYstd), _EffSYZmean(src._EffSYZmean), _EffSYZstd(src._EffSYZstd), _EffSZZmean(src._EffSZZmean), _EffSZZstd(src._EffSZZstd), _EffSZXmean(src._EffSZXmean), _EffSZXstd(src._EffSZXstd), _localVarmean(src._localVarmean), _localVarstd(src._localVarstd), - _tangentByPerturbation(src._tangentByPerturbation),_pertTol(src._pertTol), _kinematicInput(src._kinematicInput) + _tangentByPerturbation(src._tangentByPerturbation),_pertTol(src._pertTol), _kinematicInput(src._kinematicInput),_NeedExtraNorm(src._NeedExtraNorm) { #if defined(HAVE_TORCH) module = src.module; @@ -2421,18 +2421,20 @@ NonlocalDamageTorchANNBasedDG3DMaterialLaw::NonlocalDamageTorchANNBasedDG3DMater #else Msg::Error("NOT COMPILED WITH TORCH"); #endif + cLLaw = NULL; + if(src.cLLaw != NULL) + { + cLLaw=src.cLLaw->clone(); + } } NonlocalDamageTorchANNBasedDG3DMaterialLaw::~NonlocalDamageTorchANNBasedDG3DMaterialLaw() { + if (cLLaw!= NULL) delete cLLaw; } -void NonlocalDamageTorchANNBasedDG3DMaterialLaw::setNonLocalLengthScale(const int index, const CLengthLaw &cLLaw){ - IPCLength* ipCL = NULL; - cLLaw.createIPVariable(ipCL); - cLLaw.computeCL(0.,*ipCL); - characteristicLength[index] = ipCL->getConstRefToCL(); - characteristicLength[index].print("characteristicLength"); - if (ipCL!= NULL) delete ipCL; +void NonlocalDamageTorchANNBasedDG3DMaterialLaw::setNonLocalLengthScale(const int index, const CLengthLaw &_cLLaw){ + if (cLLaw) delete cLLaw; + cLLaw = _cLLaw.clone(); }; @@ -2490,7 +2492,7 @@ void NonlocalDamageTorchANNBasedDG3DMaterialLaw::initLaws(const std::map<int,mat for (int i=0; i< 6; i++) { E1_plus = E1; E1_plus(0,i) += _pertTol; - RNNstress_stiff(E1_plus, h_init, h_tmp, S_plus, EffS_plus, localVar_plus, false, DSDE, DeffSDE, DlocalVarDE); + RNNstress_stiff(E1, E1_plus, h_init, h_tmp, S_plus, EffS_plus, localVar_plus, false, DSDE, DeffSDE, DlocalVarDE); for (int j=0; j<6; j++) { DSDE(j,i) = (S_plus(0,j) - S(0,j))/_pertTol; @@ -2504,7 +2506,7 @@ void NonlocalDamageTorchANNBasedDG3DMaterialLaw::initLaws(const std::map<int,mat } else { - RNNstress_stiff(E1, h_init, h_tmp, S, EffS, _localVar, true, DSDE, DeffSDE, DlocalVarDE); + RNNstress_stiff(E1, E1, h_init, h_tmp, S, EffS, _localVar, true, DSDE, DeffSDE, DlocalVarDE); } convertFullMatrixToSTensor43(DSDE,elasticStiffness); @@ -2522,9 +2524,9 @@ void NonlocalDamageTorchANNBasedDG3DMaterialLaw::createIPState(IPStateBase* &ips bool inter=true; const MInterfaceElement *iele = dynamic_cast<const MInterfaceElement*>(ele); if(iele==NULL) inter=false; - IPVariable* ipvi = new nonLocalDamageTorchANNDG3DIPVariableDG3DIPVariable(_numberOfInternalVariables,_initialValue_h,_numNonlocalVar, characteristicLength, hasBodyForce,inter); - IPVariable* ipv1 = new nonLocalDamageTorchANNDG3DIPVariableDG3DIPVariable(_numberOfInternalVariables,_initialValue_h,_numNonlocalVar, characteristicLength, hasBodyForce,inter); - IPVariable* ipv2 = new nonLocalDamageTorchANNDG3DIPVariableDG3DIPVariable(_numberOfInternalVariables,_initialValue_h,_numNonlocalVar, characteristicLength, hasBodyForce,inter); + IPVariable* ipvi = new nonLocalDamageTorchANNDG3DIPVariableDG3DIPVariable(cLLaw,_numberOfInternalVariables,_initialValue_h,_numNonlocalVar, hasBodyForce,inter); + IPVariable* ipv1 = new nonLocalDamageTorchANNDG3DIPVariableDG3DIPVariable(cLLaw,_numberOfInternalVariables,_initialValue_h,_numNonlocalVar, hasBodyForce,inter); + IPVariable* ipv2 = new nonLocalDamageTorchANNDG3DIPVariableDG3DIPVariable(cLLaw,_numberOfInternalVariables,_initialValue_h,_numNonlocalVar, hasBodyForce,inter); if(ips != NULL) delete ips; ips = new IP3State(state_,ipvi,ipv1,ipv2); } @@ -2537,7 +2539,7 @@ void NonlocalDamageTorchANNBasedDG3DMaterialLaw::createIPVariable(IPVariable* &i if(iele == NULL){ inter=false; } - ipv = new nonLocalDamageTorchANNDG3DIPVariableDG3DIPVariable(_numberOfInternalVariables,_initialValue_h, _numNonlocalVar, characteristicLength,hasBodyForce,inter); + ipv = new nonLocalDamageTorchANNDG3DIPVariableDG3DIPVariable(cLLaw,_numberOfInternalVariables,_initialValue_h, _numNonlocalVar, hasBodyForce,inter); } @@ -2549,6 +2551,7 @@ void NonlocalDamageTorchANNBasedDG3DMaterialLaw::stress(IPVariable* ipv, const I const nonLocalDamageTorchANNDG3DIPVariableDG3DIPVariable* ipvprev = static_cast<const nonLocalDamageTorchANNDG3DIPVariableDG3DIPVariable*>(ipvp); fullMatrix<double>& E1 = ipvcur->getRefToKinematicVariables(); + static fullMatrix<double> E0(1,6); // compute E1 const STensor3& F = ipvcur->getConstRefToDeformationGradient(); static STensor3 FTF, Ustretch, R, FTF0, Ustretch0, R0; @@ -2561,6 +2564,17 @@ void NonlocalDamageTorchANNBasedDG3DMaterialLaw::stress(IPVariable* ipv, const I E1(0,3) = 0.5*FTF(0,1); E1(0,4) = 0.5*FTF(0,2); E1(0,5) = 0.5*FTF(1,2); + + if(_NeedExtraNorm){ + const STensor3& F0 = ipvprev->getConstRefToDeformationGradient(); + STensorOperation::multSTensor3FirstTranspose(F0,F0,FTF0); + E0(0,0) = 0.5*(FTF0(0,0)-1.); + E0(0,1) = 0.5*(FTF0(1,1)-1.); + E0(0,2) = 0.5*(FTF0(2,2)-1.); + E0(0,3) = 0.5*FTF0(0,1); + E0(0,4) = 0.5*FTF0(0,2); + E0(0,5) = 0.5*FTF0(1,2); + } } else if (_kinematicInput == NonlocalDamageTorchANNBasedDG3DMaterialLaw::U) { @@ -2629,7 +2643,7 @@ void NonlocalDamageTorchANNBasedDG3DMaterialLaw::stress(IPVariable* ipv, const I if (stiff && _tangentByPerturbation) { - RNNstress_stiff(E1, h0, h1, S, EffS, _localVar, false, DSDE, DeffSDE, DlocalVarDE); + RNNstress_stiff(E0, E1, h0, h1, S, EffS, _localVar, false, DSDE, DeffSDE, DlocalVarDE); static fullMatrix<double> E1_plus(1,6), S_plus(1,6), EffS_plus(1,6); static std::vector<double> localVar_plus(_numNonlocalVar); static torch::Tensor h_tmp = torch::zeros({1, 1, _numberOfInternalVariables}); @@ -2637,7 +2651,7 @@ void NonlocalDamageTorchANNBasedDG3DMaterialLaw::stress(IPVariable* ipv, const I { E1_plus = E1; E1_plus(0,i) += _pertTol; - RNNstress_stiff(E1_plus, h0, h_tmp, S_plus, EffS_plus, localVar_plus, false, DSDE, DeffSDE, DlocalVarDE); + RNNstress_stiff(E0, E1_plus, h0, h_tmp, S_plus, EffS_plus, localVar_plus, false, DSDE, DeffSDE, DlocalVarDE); for (int j=0; j<6; j++) { DSDE(j,i) = (S_plus(0,j) - S(0,j))/_pertTol; @@ -2651,7 +2665,7 @@ void NonlocalDamageTorchANNBasedDG3DMaterialLaw::stress(IPVariable* ipv, const I } else { - RNNstress_stiff(E1, h0, h1, S, EffS, _localVar, true, DSDE, DeffSDE, DlocalVarDE); + RNNstress_stiff(E0, E1, h0, h1, S, EffS, _localVar, true, DSDE, DeffSDE, DlocalVarDE); } localVar = _localVar[0]; //Damage evaluation --------- @@ -2666,7 +2680,7 @@ void NonlocalDamageTorchANNBasedDG3DMaterialLaw::stress(IPVariable* ipv, const I double nTlratio = (NonlocalVar -NonlocalVar0)/(localVar -localVar0); for(int j=0; j<6; j++){ - D(0,j) = 1.0-S(0,j)/EffS(0,j); + D(0,j) = 1.0-(S(0,j)+1.0e-20)/(EffS(0,j)+1.0e-20); nlD(0,j) = nlD0(0,j) + (D(0,j)-D0(0,j))*nTlratio; } //------------------------ end damage-------- @@ -2857,36 +2871,69 @@ void NonlocalDamageTorchANNBasedDG3DMaterialLaw::stress(IPVariable* ipv, const I #endif }; #if defined(HAVE_TORCH) -void NonlocalDamageTorchANNBasedDG3DMaterialLaw::RNNstress_stiff(const fullMatrix<double>& E1, const torch::Tensor& h0, torch::Tensor& h1, fullMatrix<double>& S, +void NonlocalDamageTorchANNBasedDG3DMaterialLaw::RNNstress_stiff(const fullMatrix<double>& E0, const fullMatrix<double>& E1, const torch::Tensor& h0, torch::Tensor& h1, fullMatrix<double>& S, fullMatrix<double>& EffS, std::vector<double>& localVar, const bool stiff, fullMatrix<double>& DSDE, fullMatrix<double>& DeffSDE, fullMatrix<double>& DlocalVarDE) { - static vector<float> E_vec(_numberOfInput); - static torch::Tensor E_norm; + vector<float> E_vec(_numberOfInput),E0_vec(_numberOfInput); + torch::Tensor E_norm, E0_norm,dE_norm; + if(!stiff){torch::NoGradGuard no_grad;} // use RNN to predict stress------------- vector<torch::jit::IValue> inputs; Normalize_strain(E1, E_vec); + if(_NeedExtraNorm){ + Normalize_strain(E0, E0_vec); + } - if(stiff){ - E_norm = torch::from_blob(E_vec.data(), {1,1, _numberOfInput}, torch::requires_grad()); + if(_NeedExtraNorm){ + if(stiff){ + E0_norm = torch::from_blob(E0_vec.data(), {1,1, _numberOfInput}, torch::requires_grad()); + } + else{ + E0_norm = torch::from_blob(E0_vec.data(), {1,1, _numberOfInput}, torch::requires_grad(false)); + } + inputs.push_back(E0_norm); } - else{ - E_norm = torch::from_blob(E_vec.data(), {1,1, _numberOfInput}); + else{ + if(stiff){ + E_norm = torch::from_blob(E_vec.data(), {1,1, _numberOfInput}, torch::requires_grad()); + } + else{ + E_norm = torch::from_blob(E_vec.data(), {1,1, _numberOfInput}, torch::requires_grad(false)); + } + inputs.push_back(E_norm); } - static vector<float> norm_E1(1); - norm_E1[0] = E1.norm(); - static torch::Tensor norm; - inputs.push_back(E_norm); - if(_kinematicInput == NonlocalDamageTorchANNBasedDG3DMaterialLaw::EGLInc){ - if(stiff){ - norm = torch::from_blob(norm_E1.data(), {1,1}, torch::requires_grad()); - } - else{ - norm = torch::from_blob(norm_E1.data(), {1,1}); - } - inputs.push_back(norm); + vector<float> norm_E1(1); + vector<float> DE_vec(_numberOfInput); + torch::Tensor norm; + + if(_kinematicInput == torchANNBasedDG3DMaterialLaw::EGLInc or _NeedExtraNorm){ + if(_kinematicInput == torchANNBasedDG3DMaterialLaw::EGL and _NeedExtraNorm){ + for (int i = 0; i <_numberOfInput; i++) { + DE_vec[i] = E_vec[i]-E0_vec[i]; + } + if(stiff){ + dE_norm = torch::from_blob(DE_vec.data(), {1,1, _numberOfInput}, torch::requires_grad()); + } + else{ + dE_norm = torch::from_blob(DE_vec.data(), {1,1, _numberOfInput}, torch::requires_grad(false)); + } + inputs.push_back(dE_norm); + norm_E1[0] = sqrt(std::inner_product(DE_vec.begin(), DE_vec.end(), DE_vec.begin(), 0.0)); + } + else{ + norm_E1[0] = E1.norm(); + } + + if(stiff){ + norm = torch::from_blob(norm_E1.data(), {1,1}, torch::requires_grad()); + } + else{ + norm = torch::from_blob(norm_E1.data(), {1,1}, torch::requires_grad(false)); + } + inputs.push_back(norm); } - inputs.push_back(h0); + inputs.push_back(h0); auto outputs= module.forward(inputs).toTuple(); torch::Tensor S_norm = outputs->elements()[0].toTensor(); diff --git a/dG3D/src/nonLocalDamageDG3DMaterialLaw.h b/dG3D/src/nonLocalDamageDG3DMaterialLaw.h index 806300e0015d50578649e4bbcf7c104a89d3da3d..6df7437ad2ede46b2420f9d0e850162dbe93f65f 100644 --- a/dG3D/src/nonLocalDamageDG3DMaterialLaw.h +++ b/dG3D/src/nonLocalDamageDG3DMaterialLaw.h @@ -794,8 +794,9 @@ class NonlocalDamageTorchANNBasedDG3DMaterialLaw : public dG3DMaterialLaw{ int _numberOfInput; int _numberOfInternalVariables; double _initialValue_h; + bool _NeedExtraNorm; int _numNonlocalVar; - std::vector<STensor3> characteristicLength; + CLengthLaw *cLLaw; #if defined(HAVE_TORCH) torch::jit::script::Module module; #endif @@ -887,7 +888,12 @@ class NonlocalDamageTorchANNBasedDG3DMaterialLaw : public dG3DMaterialLaw{ const double localVarstd,bool pert=false, double tol = 1e-5); void setNonLocalLengthScale(const int index, const CLengthLaw &cLLaw); + virtual const CLengthLaw *getCLengthLaw() const {return cLLaw; }; void setKinematicInput(const int i); + void setNeedExtraNorm (const bool NeedExtraNorm){ + _NeedExtraNorm = NeedExtraNorm; + if(_NeedExtraNorm) Msg::Info("Extra Norm is applied"); + }; #ifndef SWIG NonlocalDamageTorchANNBasedDG3DMaterialLaw(const NonlocalDamageTorchANNBasedDG3DMaterialLaw& src); virtual ~NonlocalDamageTorchANNBasedDG3DMaterialLaw(); @@ -920,8 +926,9 @@ class NonlocalDamageTorchANNBasedDG3DMaterialLaw : public dG3DMaterialLaw{ void Normalize_strain(const fullMatrix<double>& E1, vector<float>& E_norm) const; #if defined(HAVE_TORCH) void InverseNormalize_stress(const torch::Tensor& S_norm, fullMatrix<double>& S,fullMatrix<double>& EffS, std::vector<double>& localVar) const; - void RNNstress_stiff(const fullMatrix<double>& E1, const torch::Tensor& h0, torch::Tensor& h1, fullMatrix<double>& S, fullMatrix<double>& EffS, std::vector<double>& localVar, - const bool stiff, fullMatrix<double>& DSDE, fullMatrix<double>& DeffSDE, fullMatrix<double>& DlocalVarDE); + void RNNstress_stiff(const fullMatrix<double>& E0, const fullMatrix<double>& E1, const torch::Tensor& h0, torch::Tensor& h1, fullMatrix<double>& S, fullMatrix<double>& EffS, + std::vector<double>& localVar, const bool stiff, fullMatrix<double>& DSDE, fullMatrix<double>& DeffSDE, fullMatrix<double>& DlocalVarDE); + #endif #endif //SWIG