From 73c230e6245b2b7a58fe4e91816d9d80a00bf4c9 Mon Sep 17 00:00:00 2001 From: jleclerc <julien.leclerc@ulg.ac.be> Date: Fri, 18 Jan 2019 11:15:27 +0100 Subject: [PATCH 01/11] add some message in debug mode only --- NonLinearSolver/materialLaw/STensorOperations.cpp | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/NonLinearSolver/materialLaw/STensorOperations.cpp b/NonLinearSolver/materialLaw/STensorOperations.cpp index 968fc48c1..3af562176 100644 --- a/NonLinearSolver/materialLaw/STensorOperations.cpp +++ b/NonLinearSolver/materialLaw/STensorOperations.cpp @@ -12,7 +12,9 @@ bool STensorOperation::polynomial(const STensor3& a, const int order, const full static const STensor43 I4(1.,1.); if (STensorOperation::isnan(a)){ + #ifdef _DEBUG Msg::Error("STensorOperation::polynomial of a NAN tensor"); + #endif // _DEBUG STensorOperation::zero(F); if (dF){ (*dF) = I4; @@ -127,7 +129,9 @@ void STensorOperation::logSTensor3(const STensor3& a, const int order, STensor3& bool ok = STensorOperation::polynomial(ami,order,coeffs,loga,dloga,ddloga); if (!ok){ + #ifdef _DEBUG Msg::Error("error in computing logSTensor3"); + #endif // _DEBUG } } else if (order == -1) @@ -235,7 +239,9 @@ void STensorOperation::expSTensor3(const STensor3& a, const int order, STensor3& bool ok = STensorOperation::polynomial(a,order,coeffs,expa,dexpa,ddexpa); if (!ok){ + #ifdef _DEBUG Msg::Error("error in computing expSTensor3"); + #endif // _DEBUG } } else if (order == -1) -- GitLab From 57b705b173eb35c4b12fe4e3a8ab26bebdd25ec6 Mon Sep 17 00:00:00 2001 From: jleclerc <julien.leclerc@ulg.ac.be> Date: Fri, 18 Jan 2019 11:16:09 +0100 Subject: [PATCH 02/11] update checkFailed function and penalty term --- dG3D/src/dG3DCohesiveBandMaterialLaw.cpp | 48 +++++++++++++----------- 1 file changed, 26 insertions(+), 22 deletions(-) diff --git a/dG3D/src/dG3DCohesiveBandMaterialLaw.cpp b/dG3D/src/dG3DCohesiveBandMaterialLaw.cpp index 9522593c1..db83cbf7b 100644 --- a/dG3D/src/dG3DCohesiveBandMaterialLaw.cpp +++ b/dG3D/src/dG3DCohesiveBandMaterialLaw.cpp @@ -1062,7 +1062,11 @@ void PorosityCohesiveBand3DLaw::checkInternalState(IPVariable* ipv, const IPVari // interfaceIP should be deleted after satisfying failure criterion if (ipvbulk->getIPNonLocalPorosity()->isFailed()){ fipv->Delete(true); - Msg::Warning("ipv is deleted"); + Msg::Info("ipv is deleted"); + } + else if (ipvbulk->getIPNonLocalPorosity()->getYieldPorosity() < 1.e-5){ + fipv->Delete(true); + Msg::Info("ipv is deleted due to compression"); } else{ fipv->Delete(false); @@ -1092,17 +1096,17 @@ void PorosityCohesiveBand3DLaw::stress(IPVariable* ipv, const IPVariable* ipvpre const SVector3& ujump = fipv->getConstRefToJump(); const SVector3& ujumpDG = fipv->getConstRefToDGJump(); double realUJumpNormal = dot(ujump,curN) - dot(ujumpDG,curN); - if (realUJumpNormal < 0){ - fipv->setTensionFlag(true); - if(!fipvprev->ifTension()) Msg::Info("Penalty term added"); + if (realUJumpNormal < 0.){ + fipv->setTensionFlag(false); + if(fipvprev->ifTension()) Msg::Info("Penalty term added"); for (int i=0; i<3; i++){ - //T(i) += _Kp*realUJumpNormal*curN(i); - T(i) -= _Kp*realUJumpNormal*realUJumpNormal*curN(i); + T(i) += _Kp*realUJumpNormal*curN(i); + //T(i) -= _Kp*realUJumpNormal*realUJumpNormal*curN(i); } } else{ - fipv->setTensionFlag(false); - if(fipvprev->ifTension()) Msg::Info("Penalty term deleted"); + fipv->setTensionFlag(true); + if(!fipvprev->ifTension()) Msg::Info("Penalty term deleted"); } @@ -1112,12 +1116,12 @@ void PorosityCohesiveBand3DLaw::stress(IPVariable* ipv, const IPVariable* ipvpre STensorOperation::zero(fipv->getRefToDInterfaceForceDjump()); STensorOperation::zero(fipv->getRefToDInterfaceForceDjumpGradient()); - if (realUJumpNormal < 0){ + if (realUJumpNormal < 0.){ STensor3& dTdjump = fipv->getRefToDInterfaceForceDjump(); for (int i=0; i<3; i++){ for (int j=0; j<3; j++){ - //dTdjump(i,j) += _Kp*curN(i)*curN(j); - dTdjump(i,j) -= 2.*_Kp*realUJumpNormal*curN(i)*curN(j); + dTdjump(i,j) += _Kp*curN(i)*curN(j); + //dTdjump(i,j) -= 2.*_Kp*realUJumpNormal*curN(i)*curN(j); } } } @@ -1162,17 +1166,17 @@ void PorosityCohesiveBand3DLaw::stress(IPVariable* ipv, const IPVariable* ipvpre const SVector3& ujumpDG = fipv->getConstRefToDGJump(); double realUJumpNormal = dot(ujump,curN) - dot(ujumpDG,curN); - if (realUJumpNormal < 0){ - fipv->setTensionFlag(true); - if(!fipvprev->ifTension()) Msg::Info("Add penalty term"); + if (realUJumpNormal < 0.){ + fipv->setTensionFlag(false); + if(fipvprev->ifTension()) Msg::Info("Add penalty term"); for (int i=0; i<3; i++){ - T(i) -= _Kp*realUJumpNormal*realUJumpNormal*curN(i); - //T(i) += _Kp*realUJumpNormal*curN(i); + T(i) += _Kp*realUJumpNormal*curN(i); + //T(i) -= _Kp*realUJumpNormal*realUJumpNormal*curN(i); } } else{ - fipv->setTensionFlag(false); - if(fipvprev->ifTension()) Msg::Info("Delete penalty term"); + fipv->setTensionFlag(true); + if(!fipvprev->ifTension()) Msg::Info("Delete penalty term"); } @@ -1214,12 +1218,12 @@ void PorosityCohesiveBand3DLaw::stress(IPVariable* ipv, const IPVariable* ipvpre } } - if (realUJumpNormal < 0){ + if (realUJumpNormal < 0.){ for (int i=0; i<3; i++){ for (int j=0; j<3; j++){ - dTdjump(i,j) -= 2.*_Kp*realUJumpNormal*curN(i)*curN(j); - //dTdjump(i,j) += _Kp*curN(i)*curN(j); - }; + dTdjump(i,j) += _Kp*curN(i)*curN(j); + //dTdjump(i,j) -= 2.*_Kp*realUJumpNormal*curN(i)*curN(j); + } } }; -- GitLab From 9a66de9af25f7d5156106e7c021f684623513b22 Mon Sep 17 00:00:00 2001 From: jleclerc <julien.leclerc@ulg.ac.be> Date: Fri, 18 Jan 2019 11:27:07 +0100 Subject: [PATCH 03/11] add a warning message in transferInterfaceDataToBulk in CBMlaw --- dG3D/src/dG3DCohesiveBandMaterialLaw.cpp | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/dG3D/src/dG3DCohesiveBandMaterialLaw.cpp b/dG3D/src/dG3DCohesiveBandMaterialLaw.cpp index db83cbf7b..5940c6996 100644 --- a/dG3D/src/dG3DCohesiveBandMaterialLaw.cpp +++ b/dG3D/src/dG3DCohesiveBandMaterialLaw.cpp @@ -1450,7 +1450,9 @@ void PorosityCohesiveBand3DLaw::transferInterfaceDataToBulk(IPVariable* ipv, con F(i,j) += _gradFactor*(Gradjump(i,j) - Gradjump0(i,j)); } } - + + double detF = STensorOperation::determinantSTensor3(F); + if (detF < 1.e-10){Msg::Error("PorosityCohesiveBand3DLaw:: Negative jacobian= %e", detF);}; // Compute derivatives if (stiff) -- GitLab From 63f779d00d7400795b96e6dc1eedd7995652a51e Mon Sep 17 00:00:00 2001 From: jleclerc <julien.leclerc@ulg.ac.be> Date: Wed, 23 Jan 2019 10:39:57 +0100 Subject: [PATCH 04/11] add something to delete problematic IPVs --- dG3D/src/dG3DCohesiveBandMaterialLaw.cpp | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/dG3D/src/dG3DCohesiveBandMaterialLaw.cpp b/dG3D/src/dG3DCohesiveBandMaterialLaw.cpp index 5940c6996..a2ba69cfa 100644 --- a/dG3D/src/dG3DCohesiveBandMaterialLaw.cpp +++ b/dG3D/src/dG3DCohesiveBandMaterialLaw.cpp @@ -1064,14 +1064,18 @@ void PorosityCohesiveBand3DLaw::checkInternalState(IPVariable* ipv, const IPVari fipv->Delete(true); Msg::Info("ipv is deleted"); } - else if (ipvbulk->getIPNonLocalPorosity()->getYieldPorosity() < 1.e-5){ + else if (ipvbulk->getIPNonLocalPorosity()->getYieldPorosity() < 0.75*_mlawPorous->getInitialPorosity()){ + fipv->Delete(true); + Msg::Error("ipv is deleted due to too much compression: fv= %e", ipvbulk->getIPNonLocalPorosity()->getYieldPorosity()); + } + else if(ipvbulk->getJ() < 1.e-10){ + Msg::Error("ipv is deleted due to negative jacobian: J= %e", ipvbulk->getJ()); fipv->Delete(true); - Msg::Info("ipv is deleted due to compression"); } else{ fipv->Delete(false); - }; - }; + } + } }; -- GitLab From 8b84181770771b71564de918d3081020d829407b Mon Sep 17 00:00:00 2001 From: jleclerc <julien.leclerc@ulg.ac.be> Date: Wed, 23 Jan 2019 14:13:51 +0100 Subject: [PATCH 05/11] modify tol on cl --- NonLinearSolver/materialLaw/CLengthLaw.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/NonLinearSolver/materialLaw/CLengthLaw.cpp b/NonLinearSolver/materialLaw/CLengthLaw.cpp index 00ad7a5d9..1fc3e56c6 100644 --- a/NonLinearSolver/materialLaw/CLengthLaw.cpp +++ b/NonLinearSolver/materialLaw/CLengthLaw.cpp @@ -72,7 +72,7 @@ CLengthLaw * ZeroCLengthLaw::clone() const IsotropicCLengthLaw::IsotropicCLengthLaw(const int num, const double _cl, const bool init) : CLengthLaw(num,init), cl(_cl) { - double tol=1.e-10; + double tol=1.e-12; if(cl < tol) cl=tol; } IsotropicCLengthLaw::IsotropicCLengthLaw(const IsotropicCLengthLaw &source) : -- GitLab From a3c7ca9899161d487c2949ee9b8baa4957bf9d42 Mon Sep 17 00:00:00 2001 From: jleclerc <julien.leclerc@ulg.ac.be> Date: Tue, 29 Jan 2019 11:23:38 +0100 Subject: [PATCH 06/11] add a limited contribution of negative jump --- dG3D/src/dG3DCohesiveBandMaterialLaw.cpp | 87 +++++++++++++++++------- dG3D/src/dG3DCohesiveBandMaterialLaw.h | 23 ++++--- 2 files changed, 79 insertions(+), 31 deletions(-) diff --git a/dG3D/src/dG3DCohesiveBandMaterialLaw.cpp b/dG3D/src/dG3DCohesiveBandMaterialLaw.cpp index a2ba69cfa..e4432258b 100644 --- a/dG3D/src/dG3DCohesiveBandMaterialLaw.cpp +++ b/dG3D/src/dG3DCohesiveBandMaterialLaw.cpp @@ -990,12 +990,12 @@ void NonlocalPorosityCohesiveBand3DLaw::checkInternalState(IPVariable* ipv, cons - - +// -------------- PorosityCohesiveBand3DLaw PorosityCohesiveBand3DLaw::PorosityCohesiveBand3DLaw(const int num, const double h, const double Kp, const bool accountNormalBulk): Cohesive3DLaw(num,true), _h(h),_Kp(Kp), - _accountNormalBulk(accountNormalBulk), _coalescenceLaw(NULL), _mlawPorous(NULL) + _accountNormalBulk(accountNormalBulk), _enabledNegativeJumpN_limit(true),_maxPenetrationjumpN(-1.e-8),_enableDeletionOfBadConditionedIPVs(false), + _coalescenceLaw(NULL), _mlawPorous(NULL) { }; @@ -1003,7 +1003,9 @@ PorosityCohesiveBand3DLaw::PorosityCohesiveBand3DLaw(const int num, const double PorosityCohesiveBand3DLaw::PorosityCohesiveBand3DLaw(const PorosityCohesiveBand3DLaw& src): Cohesive3DLaw(src), - _h(src._h),_Kp(src._Kp), _accountNormalBulk(src._accountNormalBulk) + _h(src._h),_Kp(src._Kp), + _accountNormalBulk(src._accountNormalBulk), _enabledNegativeJumpN_limit(src._enabledNegativeJumpN_limit), + _maxPenetrationjumpN(src._maxPenetrationjumpN),_enableDeletionOfBadConditionedIPVs(src._enableDeletionOfBadConditionedIPVs) { _coalescenceLaw = NULL; if (src._coalescenceLaw!=NULL){ @@ -1046,42 +1048,56 @@ void PorosityCohesiveBand3DLaw::createIPVariable(IPVariable* &ipv,const MElement ipv = new CohesiveBand3DIPVariable(); } -void PorosityCohesiveBand3DLaw::checkInternalState(IPVariable* ipv, const IPVariable* ipvp) const{ +void PorosityCohesiveBand3DLaw::checkInternalState(IPVariable* ipv, const IPVariable* ipvp) const +{ + // Determine if the ipv is deleted or not. + /// if deleted, nothing is computed, except strains and penalty terms + /// otherwise, compute interface force and their derivatives from ipv bulk information + + // get needed ipv FractureCohesive3DIPVariable *fipv = static_cast < FractureCohesive3DIPVariable *> (ipv); const FractureCohesive3DIPVariable *fipvprev = static_cast < const FractureCohesive3DIPVariable *> (ipvp); - if (fipvprev->isDeleted()) - { + + if (fipvprev->isDeleted()){ + // If ipv is already deleted, it always stays deleted fipv->Delete(true); } else{ - // If the element is not deleted, compute interface force and their derivatives for ipv bulk information - // get ipvbulk + // If the ipv is not deleted, check if a criterion is satisfied + + // get needed ipvbulk const nonLocalPorosityDG3DIPVariable* ipvbulk = dynamic_cast<nonLocalPorosityDG3DIPVariable*>(fipv->getIPvBulk()); - // Deletion criterion: + // 1st deletion criterion: // interfaceIP should be deleted after satisfying failure criterion if (ipvbulk->getIPNonLocalPorosity()->isFailed()){ fipv->Delete(true); Msg::Info("ipv is deleted"); } - else if (ipvbulk->getIPNonLocalPorosity()->getYieldPorosity() < 0.75*_mlawPorous->getInitialPorosity()){ - fipv->Delete(true); - Msg::Error("ipv is deleted due to too much compression: fv= %e", ipvbulk->getIPNonLocalPorosity()->getYieldPorosity()); - } - else if(ipvbulk->getJ() < 1.e-10){ - Msg::Error("ipv is deleted due to negative jacobian: J= %e", ipvbulk->getJ()); - fipv->Delete(true); + // 2nd deletion criterion: + // interfaceIP should be deleted if considred as bad conditionned + else if (_enableDeletionOfBadConditionedIPVs){ + if (ipvbulk->getIPNonLocalPorosity()->getYieldPorosity() < 0.75*_mlawPorous->getInitialPorosity()){ + fipv->Delete(true); + Msg::Error("ipv is deleted due to too much compression: fv= %e", ipvbulk->getIPNonLocalPorosity()->getYieldPorosity()); + } + else if(ipvbulk->getJ() < 1.e-10){ + Msg::Error("ipv is deleted due to negative jacobian: J= %e", ipvbulk->getJ()); + fipv->Delete(true); + } } else{ + // if none of criterion is satisfied : the ipv is not deleted fipv->Delete(false); } + } }; void PorosityCohesiveBand3DLaw::stress(IPVariable* ipv, const IPVariable* ipvprev, const bool stiff, const bool checkfrac) { - // Get ipvfrac + // get ipvfrac FractureCohesive3DIPVariable *fipv = static_cast < FractureCohesive3DIPVariable *> (ipv); const FractureCohesive3DIPVariable *fipvprev = static_cast < const FractureCohesive3DIPVariable *> (ipvprev); @@ -1245,14 +1261,14 @@ void PorosityCohesiveBand3DLaw::stress(IPVariable* ipv, const IPVariable* ipvpre dirrEnergDF(k,l) += _h*dirrEnergdFinf(i,j)*dFinterfacedF(i,j,k,l); } } - } } } - + } } } + }; @@ -1391,19 +1407,21 @@ void PorosityCohesiveBand3DLaw::checkCohesiveInsertion(IPStateBase* ipsm, IPStat }; }; + + void PorosityCohesiveBand3DLaw::transferInterfaceDataToBulk(IPVariable* ipv, const IPVariable* ipvprev, const bool stiff) const { // Update kinematic data if the cohesive zone is active (== if broken), otherwise, nothing to do FractureCohesive3DIPVariable *fipv = static_cast < FractureCohesive3DIPVariable *> (ipv); const FractureCohesive3DIPVariable *fipvprev = static_cast < const FractureCohesive3DIPVariable *> (ipvprev); - if (fipv->isbroken()) - { + if (fipv->isbroken()){ CohesiveBand3DIPVariable* cohIpv = dynamic_cast<CohesiveBand3DIPVariable*>(fipv->getIPvFrac()); const CohesiveBand3DIPVariable* cohIpvprev = dynamic_cast<const CohesiveBand3DIPVariable*>(fipvprev->getIPvFrac()); // Update deformation gradient for bulk law STensor3& F = fipv->getRefToDeformationGradient(); + // - start from previous F = fipvprev->getConstRefToDeformationGradient(); @@ -1444,6 +1462,22 @@ void PorosityCohesiveBand3DLaw::transferInterfaceDataToBulk(IPVariable* ipv, con F(i,j) += (ujump(i) - ujump0(i))/_h*refN(j); } } + + // - substract normal part if asked + static SVector3 curN; + curN = fipvprev->getConstRefToCurrentOutwardNormal(); + curN.normalize(); + double jumpN = dot(ujump,curN)- dot(ujump0,curN); + + if(_enabledNegativeJumpN_limit){ + if(jumpN < _maxPenetrationjumpN){ + for (int i=0; i<3; i++){ + for (int j=0; j<3; j++){ + F(i,j) -= (jumpN-_maxPenetrationjumpN)*curN(i)/_h*refN(j); + } + } + } + } // - add jump gradient contribution double _gradFactor = 0.5; @@ -1474,6 +1508,13 @@ void PorosityCohesiveBand3DLaw::transferInterfaceDataToBulk(IPVariable* ipv, con for (int j=0; j<3; j++){ for (int k=0; k<3; k++){ dFinterfacedJump(i,j,k) += I(i,k)*refN(j)/_h; + + if(_enabledNegativeJumpN_limit){ + if(jumpN < _maxPenetrationjumpN){ + dFinterfacedJump(i,j,k) -= curN(k)*curN(i)*refN(j)/_h; + } + } + for (int l=0; l<3; l++){ dFinterfacedF(i,j,k,l) += I(i,k)*I(j,l); dFinterfacedgradJump(i,j,k,l) += _gradFactor*I(i,k)*I(j,l); @@ -1484,8 +1525,8 @@ void PorosityCohesiveBand3DLaw::transferInterfaceDataToBulk(IPVariable* ipv, con } } } - } + } }; diff --git a/dG3D/src/dG3DCohesiveBandMaterialLaw.h b/dG3D/src/dG3DCohesiveBandMaterialLaw.h index 2ac5c2d34..fee88b59d 100644 --- a/dG3D/src/dG3DCohesiveBandMaterialLaw.h +++ b/dG3D/src/dG3DCohesiveBandMaterialLaw.h @@ -81,18 +81,25 @@ class NonlocalPorosityCohesiveBand3DLaw : public CohesiveBand3DLaw{ class PorosityCohesiveBand3DLaw : public Cohesive3DLaw{ #ifndef SWIG - protected: - double _h; // cohesive band thickness - double _Kp; // penalty parameter to avoid penetration - bool _accountNormalBulk; //false if the normal part of bulk defo is not accounted for in interfacical defo - // true otherwise - CoalescenceLaw* _coalescenceLaw; // will be deleted sooon - mlawNonLocalPorosity* _mlawPorous; +protected: + // Cohesive band parameters + double _h; // cohesive band thickness + double _Kp; // penalty parameter to avoid penetration + + // Options and numerical parameters + bool _accountNormalBulk; //false if the normal part of bulk defo is not accounted for in interfacical defo + // true otherwise (==all the components of the bulk F is used...) + bool _enabledNegativeJumpN_limit; // true if the normal jump is limited in the negative values + double _maxPenetrationjumpN; // maximal (negative!) value of admissible penetration jump + bool _enableDeletionOfBadConditionedIPVs; // true if bad-conditionned ipvs are deleted for convergence purpose + // Porous law pointer (for mlaw-dependent case) + mlawNonLocalPorosity* _mlawPorous; + CoalescenceLaw* _coalescenceLaw; // will be deleted sooon #endif public: - PorosityCohesiveBand3DLaw(const int num, const double h, const double Kp, const bool accountNormalBulk); + PorosityCohesiveBand3DLaw(const int num, const double h, const double Kp, const bool accountNormalBulk=true); void setCoalescenceLaw(const CoalescenceLaw& coalesLaw); void setPorousMaterialLaw(const NonLocalPorousDuctileDamageDG3DMaterialLaw& dG3DmlawPorous); #ifndef SWIG -- GitLab From 3f04a7a58df3690f854f8474572104d01d9c836f Mon Sep 17 00:00:00 2001 From: jleclerc <julien.leclerc@ulg.ac.be> Date: Tue, 29 Jan 2019 11:34:04 +0100 Subject: [PATCH 07/11] add options settings --- dG3D/src/dG3DCohesiveBandMaterialLaw.cpp | 25 ++++++++++++++++++++++++ dG3D/src/dG3DCohesiveBandMaterialLaw.h | 2 ++ 2 files changed, 27 insertions(+) diff --git a/dG3D/src/dG3DCohesiveBandMaterialLaw.cpp b/dG3D/src/dG3DCohesiveBandMaterialLaw.cpp index e4432258b..7b9187129 100644 --- a/dG3D/src/dG3DCohesiveBandMaterialLaw.cpp +++ b/dG3D/src/dG3DCohesiveBandMaterialLaw.cpp @@ -1030,6 +1030,31 @@ void PorosityCohesiveBand3DLaw::setPorousMaterialLaw(const NonLocalPorousDuctile _coalescenceLaw = _mlawPorous->getCoalescenceLaw()->clone(); }; +void PorosityCohesiveBand3DLaw::setMaximalAdmissiblePenetrationJump(bool activationFlag, double maxNegJump){ + _enabledNegativeJumpN_limit = activationFlag; + _maxPenetrationjumpN = maxNegJump; + if (maxNegJump > 0.){ + Msg::Fatal("PorosityCohesiveBand3DLaw::the admissible value shloud be negative instead of %e",maxNegJump); + } + + if (_enabledNegativeJumpN_limit){ + Msg::Info("PorosityCohesiveBand3DLaw:: the normal penetration jump is limited to %e",maxNegJump); + } + else{ + Msg::Info("PorosityCohesiveBand3DLaw:: the normal penetration jump is not limited"); + } +} + +void PorosityCohesiveBand3DLaw::setDeletionOfBadConditionnedIPV(bool activationFlag){ + _enableDeletionOfBadConditionedIPVs = activationFlag; + if (_enableDeletionOfBadConditionedIPVs){ + Msg::Info("PorosityCohesiveBand3DLaw:: the bad-conditionned ipvs will be deleted"); + } + else{ + Msg::Info("PorosityCohesiveBand3DLaw:: the bad-conditionned ipvs will NOT be deleted"); + } +} + void PorosityCohesiveBand3DLaw::createIPState(IPStateBase* &ips,const bool* state_, const MElement *ele, const int nbFF_, const IntPt *GP, diff --git a/dG3D/src/dG3DCohesiveBandMaterialLaw.h b/dG3D/src/dG3DCohesiveBandMaterialLaw.h index fee88b59d..964853aa0 100644 --- a/dG3D/src/dG3DCohesiveBandMaterialLaw.h +++ b/dG3D/src/dG3DCohesiveBandMaterialLaw.h @@ -102,6 +102,8 @@ protected: PorosityCohesiveBand3DLaw(const int num, const double h, const double Kp, const bool accountNormalBulk=true); void setCoalescenceLaw(const CoalescenceLaw& coalesLaw); void setPorousMaterialLaw(const NonLocalPorousDuctileDamageDG3DMaterialLaw& dG3DmlawPorous); + void setMaximalAdmissiblePenetrationJump(bool activationFlag, double maxNegJump=-1.); + void setDeletionOfBadConditionnedIPV(bool activationFlag); #ifndef SWIG PorosityCohesiveBand3DLaw(const PorosityCohesiveBand3DLaw& src); virtual ~PorosityCohesiveBand3DLaw(){ -- GitLab From 9589d0f64a8cca6e3a11c5532ab8033069fcfbff Mon Sep 17 00:00:00 2001 From: jleclerc <julien.leclerc@ulg.ac.be> Date: Tue, 29 Jan 2019 11:35:06 +0100 Subject: [PATCH 08/11] adapt printed message --- dG3D/src/dG3DCohesiveBandMaterialLaw.cpp | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/dG3D/src/dG3DCohesiveBandMaterialLaw.cpp b/dG3D/src/dG3DCohesiveBandMaterialLaw.cpp index 7b9187129..e6aaedd79 100644 --- a/dG3D/src/dG3DCohesiveBandMaterialLaw.cpp +++ b/dG3D/src/dG3DCohesiveBandMaterialLaw.cpp @@ -1514,9 +1514,11 @@ void PorosityCohesiveBand3DLaw::transferInterfaceDataToBulk(IPVariable* ipv, con } } - double detF = STensorOperation::determinantSTensor3(F); - if (detF < 1.e-10){Msg::Error("PorosityCohesiveBand3DLaw:: Negative jacobian= %e", detF);}; - + if(!fipvprev->isDeleted()){ + double detF = STensorOperation::determinantSTensor3(F); + if (detF < 1.e-10){Msg::Error("PorosityCohesiveBand3DLaw:: Negative jacobian= %e", detF);}; + } + // Compute derivatives if (stiff) { -- GitLab From 0c637c984dd805734cf32c638972dc2e80f6566a Mon Sep 17 00:00:00 2001 From: jleclerc <julien.leclerc@ulg.ac.be> Date: Fri, 1 Feb 2019 15:54:25 +0100 Subject: [PATCH 09/11] update first version of predictor update --- NonLinearSolver/Domain/partDomain.h | 3 + .../nlsolver/implicitHulbertChungPETSc.h | 26 +++++- .../nlsolver/nonLinearMechSolver.cpp | 93 +++++++++++++++---- .../nlsolver/nonLinearMechSolver.h | 2 + NonLinearSolver/nlsolver/nonLinearSystems.h | 1 + dG3D/src/dG3DDomain.cpp | 44 +++++++++ dG3D/src/dG3DDomain.h | 1 + 7 files changed, 151 insertions(+), 19 deletions(-) diff --git a/NonLinearSolver/Domain/partDomain.h b/NonLinearSolver/Domain/partDomain.h index d347b1ae6..113e49771 100644 --- a/NonLinearSolver/Domain/partDomain.h +++ b/NonLinearSolver/Domain/partDomain.h @@ -198,6 +198,9 @@ public: virtual FunctionSpaceBase* getFunctionSpace() =0; virtual const FunctionSpaceBase* getFunctionSpace() const=0; + + virtual void getDofKeyswithZeroPredictor(std::vector<Dof>& DofVector) const{}; + // some data of BC have to be set by domain virtual FunctionSpaceBase* getSpaceForBC(const nonLinearBoundaryCondition::type bc_type, const nonLinearBoundaryCondition::location bc_location, const int dof_comp,const groupOfElements *groupBC) const=0; // for dirichlet, neumann and initial diff --git a/NonLinearSolver/nlsolver/implicitHulbertChungPETSc.h b/NonLinearSolver/nlsolver/implicitHulbertChungPETSc.h index c9b8fd678..03c2126fe 100644 --- a/NonLinearSolver/nlsolver/implicitHulbertChungPETSc.h +++ b/NonLinearSolver/nlsolver/implicitHulbertChungPETSc.h @@ -85,10 +85,14 @@ class implicitHulbertChungPetsc : public nonLinearSystem<scalar>, bool _whichStep; double _alpham,_alphaf, _beta,_gamma; // parameters of scheme double _timeStep; // value of time step - int _nbRows; // To know the system size. - + int _nbRows; // To know the system size implicitStepPetsc<scalar> *_currentStep, *_previousStep, *_initialStep, *_step1, *_step2; Mat _M; // mass matrix + + // Predictor management + int _nbZeroPredictorDofs; // number of dof for which the predictor is equal to the previous step + std::vector<int> _zeroPredictorDofKeys; // Keys of concerned dof + std::vector<scalar> _zeros; // Vectors of zeros int _numIteration; double OneMinusAlphaf; // 1 - alphaf @@ -110,6 +114,7 @@ class implicitHulbertChungPetsc : public nonLinearSystem<scalar>, implicitHulbertChungPetsc(double alpham=0., double alphaf = 1., double beta=0., double gamma=0.5): _alphaf(alphaf), _alpham(alpham),_beta(beta),_gamma(gamma), linearSystemPETSc<scalar>(),_whichStep(true),_timeStep(0.), _nbRows(0), + _nbZeroPredictorDofs(0), _zeroPredictorDofKeys(), _zeros(), _numIteration(0),_flagb(false),_initialStep(NULL){ OneMinusAlphaf = 1. - _alphaf; AlphafMinusOne = -1.* OneMinusAlphaf; @@ -154,6 +159,16 @@ class implicitHulbertChungPetsc : public nonLinearSystem<scalar>, implicitStepPetsc<scalar>& dst = this->getState(destination); dst = src; } + + virtual void addZeroPredictorDofKeys(const std::vector<int>& dofKeys){ + // Copy the vector dof keys and set _nbZeroPredictorDofs and _zeros accordingly + _zeroPredictorDofKeys.clear(); + _zeroPredictorDofKeys = dofKeys; + _nbZeroPredictorDofs = _zeroPredictorDofKeys.size(); + _zeros.clear(); + _zeros.assign(_nbZeroPredictorDofs,0.); + } + void nextStep(){ if(_whichStep){ @@ -463,6 +478,13 @@ class implicitHulbertChungPetsc : public nonLinearSystem<scalar>, _try(VecAssemblyBegin(_currentStep->_xddot)); _try(VecAssemblyEnd(_currentStep->_xddot)); _try(VecZeroEntries(this->_currentStep->_xddot)); + + if(_nbZeroPredictorDofs > 0){ + // Reset the predictor state to zero for some of the dofs + /* check if this is at the good place */ + _try(VecSetValues(_currentStep->_xdot , _nbZeroPredictorDofs, &_zeroPredictorDofKeys[0], &_zeros[0], INSERT_VALUES)); + _try(VecSetValues(_currentStep->_xddot, _nbZeroPredictorDofs, &_zeroPredictorDofKeys[0], &_zeros[0], INSERT_VALUES)); + } } else{ _try(MatAssemblyBegin(this->_a,MAT_FINAL_ASSEMBLY)); diff --git a/NonLinearSolver/nlsolver/nonLinearMechSolver.cpp b/NonLinearSolver/nlsolver/nonLinearMechSolver.cpp index 98c025c48..b4f3bed24 100644 --- a/NonLinearSolver/nlsolver/nonLinearMechSolver.cpp +++ b/NonLinearSolver/nlsolver/nonLinearMechSolver.cpp @@ -3695,6 +3695,12 @@ void nonLinearMechSolver::velocityBC(std::string onwhat, const int numphys, cons allDirichlet.back()._mycondition = nonLinearBoundaryCondition::velocity; } +void nonLinearMechSolver::accelerationBC(std::string onwhat, const int numphys, const int comp, const double value) +{ + this->displacementBC(onwhat,numphys,comp,value); + allDirichlet.back()._mycondition = nonLinearBoundaryCondition::acceleration; +} + void nonLinearMechSolver::displacementBC(std::string onwhat, const int numphys, const int comp, const double value, elementFilter *filter) { simpleFunction<double> *fct = new simpleFunctionTime<double>(value); @@ -5798,15 +5804,68 @@ void nonLinearMechSolver::initializeStaticScheme() AssembleMass(cdom->getMassTerm(),rcsp,pAssembler); } } - } + } } - std::string name = "A"; linearSystem<double>* lsys =pAssembler->getLinearSystem(name); - nonLinearSystem<double>* nonsys = dynamic_cast<nonLinearSystem<double>*>(lsys); + nonLinearSystem<double>* nlsys = dynamic_cast<nonLinearSystem<double>*>(lsys); pathFollowingSystem<double>* pathSys = dynamic_cast<pathFollowingSystem<double>*>(lsys); + if (whatScheme == Implicit){ //now only for implicit and dgdomain + + // Collect all dofs number associated with components for which the predictor + // (in the time integration scheme) needs to be zero (== equal to previous state) + std::vector<int> DofKeysContainer; // vector of all the considered Dof + std::vector<Dof> DofForOneDomain; // vector for one domain + + /// Loop on all domains + for(std::vector<partDomain*>::iterator it=domainVector.begin(); it!=domainVector.end(); ++it){ + partDomain *dom = *it; + // take back the keys of dofs + dom->getDofKeyswithZeroPredictor(DofForOneDomain); + for(int i = 0; i < DofForOneDomain.size(); ++i){ + // get the position of the dof and store it in the vector + int posDoF = pAssembler->getDofNumber(DofForOneDomain[i]); + if(posDoF >= 0){ + DofKeysContainer.push_back(posDoF); + } + } + } + // transfer dof to non-linear system + nlsys->addZeroPredictorDofKeys(DofKeysContainer); + + + /* + /// Loop on all domains and elements + for(std::vector<partDomain*>::iterator it=domainVector.begin(); it!=domainVector.end(); ++it){ + partDomain *dom = *it; + for(groupOfElements::elementContainer::iterator el_it=dom->g_cbegin(); el_it!= dom->g_cend(); ++el_it){ + MElement* ele = *el_it; + /// get all the dof keys of an element + std::vector<Dof> R; + dom->getFunctionSpace()->getKeys(ele,R); + + /// Loop on the considered components + for(std::vector<int>::iterator comp_it=_noMassMechanicalDofComp.begin(); comp_it!=_noMassMechanicalDofComp.end(); ++comp_it){ + /// create the filter for the desired components + FilterDof* compDofFilter = dom->createFilterDof(*comp_it);/// (maybe do it only once, when the domain is called and create a vector ?) + for(int i = 0; i < R.size(); ++i){ + // if the actual considered dof is associated with the sought component, add its dof number in the container + if(compDofFilter->operator()(R[i])){ + int posDoF = pAssembler->getDofNumber(R[i]); + if(posDoF >= 0){ + DofKeysContainer.push_back(posDoF); + } + } + } + } + } + } + */ + + } + if (_pathFollowing){ this->computeLoadVector(_ufield); } @@ -6595,7 +6654,7 @@ int nonLinearMechSolver::NewtonRaphson(dofManager<double> *pmanager,unknownField break; } } - + //this->fixNodalDofs(); // reset values of converged state in accordance to BC return iter; } @@ -12142,20 +12201,20 @@ void nonLinearMechSolver::readPBCDataFromFile(const std::string fn){ void nonLinearMechSolver::saveStiffnessMatrixToFile(const int iter){ -linearSystem<double>* lsys=NULL; -std::string name="A"; -lsys=pAssembler->getLinearSystem(name); + linearSystem<double>* lsys=NULL; + std::string name="A"; + lsys=pAssembler->getLinearSystem(name); #if defined(HAVE_PETSC) -linearSystemPETSc<double>* lpet=dynamic_cast<linearSystemPETSc<double>*>(lsys); -if(lpet==NULL)return; -else{ - functionPETSc::MatToFile(lpet->getMatrix(),"data/Mat_Stiffness.txt"); -} -implicitHulbertChungPetsc<double>* hcsys=dynamic_cast<implicitHulbertChungPetsc<double>*>(lsys); -if(hcsys==NULL)return; -else{ - functionPETSc::MatToFile(hcsys->getMassMatrix(),"data/Mat_Mass.txt"); -} + linearSystemPETSc<double>* lpet=dynamic_cast<linearSystemPETSc<double>*>(lsys); + if(lpet==NULL)return; + else{ + functionPETSc::MatToFile(lpet->getMatrix(),"data/Mat_Stiffness.txt"); + } + implicitHulbertChungPetsc<double>* hcsys=dynamic_cast<implicitHulbertChungPetsc<double>*>(lsys); + if(hcsys==NULL)return; + else{ + functionPETSc::MatToFile(hcsys->getMassMatrix(),"data/Mat_Mass.txt"); + } #endif }; diff --git a/NonLinearSolver/nlsolver/nonLinearMechSolver.h b/NonLinearSolver/nlsolver/nonLinearMechSolver.h index 1bc13d9d2..f2d8232f3 100644 --- a/NonLinearSolver/nlsolver/nonLinearMechSolver.h +++ b/NonLinearSolver/nlsolver/nonLinearMechSolver.h @@ -358,6 +358,7 @@ class nonLinearMechSolver selectiveUpdateBase* _selectiveObject; + #ifndef SWIG void commonModel(); // common operation for loadModel and createModel void initMapRanks(const std::set<int>& parts); @@ -601,6 +602,7 @@ class nonLinearMechSolver virtual void velocityBC(std::string onwhat, const int numphys, const int comp, const double value,elementFilter *filter); virtual void velocityBC(std::string onwhat, const int numphys, const int comp, simpleFunction<double> *fct); virtual void velocityBC(std::string onwhat, const int numphys, const int comp, elementFilter *filter, simpleFunction<double> *fct); + virtual void accelerationBC(std::string onwhat, const int numphys, const int comp, const double value); virtual void constraintBC(std::string onwhat, const int numphys, const int comp); virtual void periodicBC(std::string onwhat, const int phys1, const int phys2, const int v1, const int v2, const int comp); diff --git a/NonLinearSolver/nlsolver/nonLinearSystems.h b/NonLinearSolver/nlsolver/nonLinearSystems.h index bf2e26675..6da50cfb5 100644 --- a/NonLinearSolver/nlsolver/nonLinearSystems.h +++ b/NonLinearSolver/nlsolver/nonLinearSystems.h @@ -80,6 +80,7 @@ class nonLinearSystem{ virtual void getFromSolution(int row, scalar &val, const nonLinearBoundaryCondition::whichCondition wv) const=0; virtual void getFromMatrix(int row, int col, scalar &val,const typename nonLinearSystem<scalar>::whichMatrix wm) const=0; virtual void printAll() {}; + virtual void addZeroPredictorDofKeys(const std::vector<int>& dofKeys){Msg::Fatal("nonLinearSystem::addZeroPredictorDofKeys should not be called");}; #if defined(HAVE_MPI) virtual double* getArrayIndexPreviousRightHandSidePlusMPI(const int index)=0; // virtual void setPreviousRightHandSidePlusMPI(const int row, const scalar &val)=0; diff --git a/dG3D/src/dG3DDomain.cpp b/dG3D/src/dG3DDomain.cpp index c90dcff51..bd96a00df 100644 --- a/dG3D/src/dG3DDomain.cpp +++ b/dG3D/src/dG3DDomain.cpp @@ -28,6 +28,7 @@ #include "nonLocalDamageDG3DIPVariable.h" + dG3DDomain::dG3DDomain (const int tag,const int phys, const int ws, const int lnum, const int fdg, const int dim, int nonLocalVar, int nonConstitutiveExtraDOFVar) : dgPartDomain(tag,phys,ws,fdg,nonLocalVar, nonConstitutiveExtraDOFVar), _lnum(lnum), @@ -456,6 +457,8 @@ QuadratureBase* dG3DDomain::getQuadratureRulesForNeumannBC(nonLinearNeumannBC &n return integ; } + + void dG3DDomain::createTerms(unknownField *uf,IPField*ip) { FunctionSpace<double>* dgspace = static_cast<FunctionSpace<double>*>(_space); @@ -693,6 +696,47 @@ BilinearTermBase* dG3DDomain::createNeumannMatrixTerm(FunctionSpace<double> *spn return new BiNonLinearTermVoid(); } } + + + +void dG3DDomain::getDofKeyswithZeroPredictor(std::vector<Dof>& DofVector) const +{ + // Send back the dof for which the predictor must be set to zero + DofVector.clear(); // clean it first + + // Create a vector of filters for desired dofs + /// /!\ now, it is only non-local dofs + std::vector<FilterDof*> FilterVec; + for(int i = 0; i < getNumNonLocalVariable(); ++i){ + FilterVec.push_back(this->createFilterDof(i+3)); + } + + // Loop on all elements + std::vector<Dof> R; + const partDomain *dom = static_cast<const partDomain*>(this); + for(groupOfElements::elementContainer::const_iterator it=gi->begin(); it!=gi->end();++it){ + // get dofs of an element + MElement *ele = *it; + R.clear(); + _space->getKeys(ele,R); + + // Loop and all dof and apply all the filters + int nbDof = _space->getNumKeys(ele); + for(int i = 0; i < nbDof; ++i){ + for(std::vector<FilterDof*>::const_iterator it = FilterVec.cbegin(); it != FilterVec.cend(); ++it){ + FilterDof* filter = *it; + if(filter->operator()(R[i])){ + // if the dof belongs to the desired components, add it to the list and break the filter loop + DofVector.push_back(R[i]); + break; + } + } + } + } +} + + + void dG3DDomain::computeAllIPStrain(AllIPState *aips,const unknownField *ufield,const IPStateBase::whichState ws, const bool virt){ IntPt *GP; fullVector<double> dispm; diff --git a/dG3D/src/dG3DDomain.h b/dG3D/src/dG3DDomain.h index d0804496e..e872885d6 100644 --- a/dG3D/src/dG3DDomain.h +++ b/dG3D/src/dG3DDomain.h @@ -178,6 +178,7 @@ class dG3DDomain : public dgPartDomain{ virtual int getMinusLawNum() const{return _lnum;} virtual int getPlusLawNum() const{return _lnum;} virtual void setGaussIntegrationRule(); + virtual void getDofKeyswithZeroPredictor(std::vector<Dof>& DofVector) const; // get back the dof for which the predictor must be set to zero virtual void createTerms(unknownField *uf,IPField*ip); virtual void initializeTerms(unknownField *uf,IPField *ip); virtual LinearTermBase<double>* createNeumannTerm(FunctionSpace<double> *spneu, const groupOfElements* g, -- GitLab From 21fdcd5603403c27e479eb6454c234ee1845e7a8 Mon Sep 17 00:00:00 2001 From: jleclerc <julien.leclerc@ulg.ac.be> Date: Mon, 4 Feb 2019 11:44:14 +0100 Subject: [PATCH 10/11] correct mistake in parenthesis and indentation --- dG3D/src/dG3DTerms.cpp | 34 +++++++++++++++------------------- 1 file changed, 15 insertions(+), 19 deletions(-) diff --git a/dG3D/src/dG3DTerms.cpp b/dG3D/src/dG3DTerms.cpp index a7e7c3b4d..14d7a30e1 100644 --- a/dG3D/src/dG3DTerms.cpp +++ b/dG3D/src/dG3DTerms.cpp @@ -1257,7 +1257,7 @@ void dG3DForceInter::get(MElement *ele, int npts, IntPt *GP, fullVector<double> if( ipvm->getNumberNonLocalVariable()>getNumNonLocalVariable() ) Msg::Fatal("Your minus material law uses more non local variable than your domain dG3DForceInter::get"); if( ipvp->getNumberNonLocalVariable()>getNumNonLocalVariable() ) Msg::Fatal("Your plus material law uses more non local variable than your domain dG3DForceInter::get"); if( ipvm->getNumberNonLocalVariable()!= ipvm->getNumberNonLocalVariable() ) Msg::Fatal("Your minus and plus material laws do not have the same number of non local variables dG3DForceInter::get"); - if( getNonLocalContinuity() && !broken) + if( getNonLocalContinuity()) { // characteristic size const double nbetahs = getNonLocalStabilityParameter()/hs; @@ -1361,7 +1361,7 @@ void dG3DForceInter::get(MElement *ele, int npts, IntPt *GP, fullVector<double> Msg::Fatal("Your plus and minus material laws do not use the same number of constitutive extra dof variables dG3DForceInter::get"); for (int extraDOFField = 0; extraDOFField< ipvp->getNumConstitutiveExtraDofDiffusionVariable(); extraDOFField++){ for (int extraDOFField2 = 0; extraDOFField2< ipvp->getNumConstitutiveExtraDofDiffusionVariable(); extraDOFField2++){ - if( getConstitutiveExtraDofDiffusionContinuity() && !broken) + if( getConstitutiveExtraDofDiffusionContinuity()) { const double nbetahs = getConstitutiveExtraDofDiffusionStabilityParameter()/hs; double fieldJump; @@ -1494,19 +1494,19 @@ void dG3DForceInter::get(MElement *ele, int npts, IntPt *GP, fullVector<double> } //new if(extraDOFField==0){ - double gamma=0.; - for(int j=0;j<nbFFm;j++) - { - m(j+indexField*nbFFm) += gamma*( Stiff_alphadialitationm_JumpN*Valsm[j+indexField*nbFFm])*(wJ/2.)*(-1.); - } - for(int j=0;j<nbFFp;j++) - { - m(j+indexField*nbFFp+nbdofm) += gamma*( Stiff_alphadialitationp_JumpN*Valsp[j+indexField*nbFFp])*(wJ/2.)*(-1.); + double gamma=0.; + for(int j=0;j<nbFFm;j++){ + m(j+indexField*nbFFm) += gamma*( Stiff_alphadialitationm_JumpN*Valsm[j+indexField*nbFFm])*(wJ/2.)*(-1.); + } + for(int j=0;j<nbFFp;j++){ + m(j+indexField*nbFFp+nbdofm) += gamma*( Stiff_alphadialitationp_JumpN*Valsp[j+indexField*nbFFp])*(wJ/2.)*(-1.); + } } - } + } + } + } } - else //broken case - { + else{ // broken case const FractureCohesive3DIPVariable *ipvmf = static_cast<const FractureCohesive3DIPVariable*>(ipvm); const FractureCohesive3DIPVariable *ipvpf = static_cast<const FractureCohesive3DIPVariable*>(ipvp); @@ -1529,14 +1529,10 @@ void dG3DForceInter::get(MElement *ele, int npts, IntPt *GP, fullVector<double> m(j+k*nbFFp+nbdofm) += interfaceForce[k]*(Valsp[j+0*nbFFp]*wJ); } } - } } - } - } - } } - //} - //m.print("Force inter"); + } +//m.print("Force inter"); } void dG3DStiffnessInter::get(MElement *ele,int npts,IntPt *GP,fullMatrix<double> &stiff) const -- GitLab From 046b20159f32d1fbda66eb543d71b3e7a901c0cd Mon Sep 17 00:00:00 2001 From: jleclerc <julien.leclerc@ulg.ac.be> Date: Mon, 4 Feb 2019 13:41:43 +0100 Subject: [PATCH 11/11] clean code for the selection of code --- .../nlsolver/implicitHulbertChungPETSc.h | 30 +++++++++++++++-- .../nlsolver/nonLinearMechSolver.cpp | 33 ++----------------- dG3D/src/dG3DDomain.cpp | 32 ++++++++++++++---- 3 files changed, 54 insertions(+), 41 deletions(-) diff --git a/NonLinearSolver/nlsolver/implicitHulbertChungPETSc.h b/NonLinearSolver/nlsolver/implicitHulbertChungPETSc.h index 03c2126fe..32a91d908 100644 --- a/NonLinearSolver/nlsolver/implicitHulbertChungPETSc.h +++ b/NonLinearSolver/nlsolver/implicitHulbertChungPETSc.h @@ -467,6 +467,29 @@ class implicitHulbertChungPetsc : public nonLinearSystem<scalar>, _try(VecAssemblyEnd(this->_previousStep->_xdot)); } #endif // HAVE_MPI +/* + // Version with predictor = 0 + _try(VecCopy(this->_previousStep->_x,this->_currentStep->_x)); + + _try(VecCopy(this->_previousStep->_xdot,this->_currentStep->_xdot)); + if(_nbZeroPredictorDofs > 0){ + _try(VecSetValues(_currentStep->_xdot, _nbZeroPredictorDofs, &_zeroPredictorDofKeys[0], &_zeros[0], INSERT_VALUES)); + } + _try(VecAXPY(this->_currentStep->_x,_timeStep,this->_currentStep->_xdot)); + + _try(VecAssemblyBegin(_currentStep->_xddot)); + _try(VecAssemblyEnd(_currentStep->_xddot)); + _try(VecCopy(this->_previousStep->_xddot,this->_currentStep->_xddot)); + if(_nbZeroPredictorDofs > 0){ + _try(VecSetValues(_currentStep->_xddot, _nbZeroPredictorDofs, &_zeroPredictorDofKeys[0], &_zeros[0], INSERT_VALUES)); + } + _try(VecAXPY(this->_currentStep->_x,Deltat2halfMinusbeta,this->_currentStep->_xddot)); + _try(VecAXPY(this->_currentStep->_xdot,DeltatOneMinusGamma,this->_currentStep->_xddot)); + + _try(VecZeroEntries(this->_currentStep->_xddot)); + +*/ + _try(VecCopy(this->_previousStep->_x,this->_currentStep->_x)); _try(VecAXPY(this->_currentStep->_x,_timeStep,this->_previousStep->_xdot)); @@ -478,13 +501,14 @@ class implicitHulbertChungPetsc : public nonLinearSystem<scalar>, _try(VecAssemblyBegin(_currentStep->_xddot)); _try(VecAssemblyEnd(_currentStep->_xddot)); _try(VecZeroEntries(this->_currentStep->_xddot)); - + +/* if(_nbZeroPredictorDofs > 0){ - // Reset the predictor state to zero for some of the dofs - /* check if this is at the good place */ + // Reset the predictor state to zero for some of the dofs (xdot and xddot) _try(VecSetValues(_currentStep->_xdot , _nbZeroPredictorDofs, &_zeroPredictorDofKeys[0], &_zeros[0], INSERT_VALUES)); _try(VecSetValues(_currentStep->_xddot, _nbZeroPredictorDofs, &_zeroPredictorDofKeys[0], &_zeros[0], INSERT_VALUES)); } +*/ } else{ _try(MatAssemblyBegin(this->_a,MAT_FINAL_ASSEMBLY)); diff --git a/NonLinearSolver/nlsolver/nonLinearMechSolver.cpp b/NonLinearSolver/nlsolver/nonLinearMechSolver.cpp index b4f3bed24..537e11774 100644 --- a/NonLinearSolver/nlsolver/nonLinearMechSolver.cpp +++ b/NonLinearSolver/nlsolver/nonLinearMechSolver.cpp @@ -5813,7 +5813,7 @@ void nonLinearMechSolver::initializeStaticScheme() pathFollowingSystem<double>* pathSys = dynamic_cast<pathFollowingSystem<double>*>(lsys); if (whatScheme == Implicit){ //now only for implicit and dgdomain - + Msg::Info("start dof checking"); // Collect all dofs number associated with components for which the predictor // (in the time integration scheme) needs to be zero (== equal to previous state) std::vector<int> DofKeysContainer; // vector of all the considered Dof @@ -5834,36 +5834,7 @@ void nonLinearMechSolver::initializeStaticScheme() } // transfer dof to non-linear system nlsys->addZeroPredictorDofKeys(DofKeysContainer); - - - /* - /// Loop on all domains and elements - for(std::vector<partDomain*>::iterator it=domainVector.begin(); it!=domainVector.end(); ++it){ - partDomain *dom = *it; - for(groupOfElements::elementContainer::iterator el_it=dom->g_cbegin(); el_it!= dom->g_cend(); ++el_it){ - MElement* ele = *el_it; - /// get all the dof keys of an element - std::vector<Dof> R; - dom->getFunctionSpace()->getKeys(ele,R); - - /// Loop on the considered components - for(std::vector<int>::iterator comp_it=_noMassMechanicalDofComp.begin(); comp_it!=_noMassMechanicalDofComp.end(); ++comp_it){ - /// create the filter for the desired components - FilterDof* compDofFilter = dom->createFilterDof(*comp_it);/// (maybe do it only once, when the domain is called and create a vector ?) - for(int i = 0; i < R.size(); ++i){ - // if the actual considered dof is associated with the sought component, add its dof number in the container - if(compDofFilter->operator()(R[i])){ - int posDoF = pAssembler->getDofNumber(R[i]); - if(posDoF >= 0){ - DofKeysContainer.push_back(posDoF); - } - } - } - } - } - } - */ - + Msg::Info("Dof checking for predictor done: %i DoF are with a predictor = 0",DofKeysContainer.size()); } if (_pathFollowing){ diff --git a/dG3D/src/dG3DDomain.cpp b/dG3D/src/dG3DDomain.cpp index bd96a00df..7670448da 100644 --- a/dG3D/src/dG3DDomain.cpp +++ b/dG3D/src/dG3DDomain.cpp @@ -704,25 +704,25 @@ void dG3DDomain::getDofKeyswithZeroPredictor(std::vector<Dof>& DofVector) const // Send back the dof for which the predictor must be set to zero DofVector.clear(); // clean it first +/* // Create a vector of filters for desired dofs /// /!\ now, it is only non-local dofs std::vector<FilterDof*> FilterVec; for(int i = 0; i < getNumNonLocalVariable(); ++i){ - FilterVec.push_back(this->createFilterDof(i+3)); + FilterVec.push_back(this->createFilterDof(getDim()+i)); } - // Loop on all elements std::vector<Dof> R; const partDomain *dom = static_cast<const partDomain*>(this); - for(groupOfElements::elementContainer::const_iterator it=gi->begin(); it!=gi->end();++it){ + for(groupOfElements::elementContainer::const_iterator it=g->begin(); it!=g->end();++it){ // get dofs of an element MElement *ele = *it; R.clear(); _space->getKeys(ele,R); - // Loop and all dof and apply all the filters - int nbDof = _space->getNumKeys(ele); - for(int i = 0; i < nbDof; ++i){ + // Loop and all dof and apply all the filters + //int nbDof = _space->getNumKeys(ele); + for(int i = 0; i < R.size(); ++i){ for(std::vector<FilterDof*>::const_iterator it = FilterVec.cbegin(); it != FilterVec.cend(); ++it){ FilterDof* filter = *it; if(filter->operator()(R[i])){ @@ -733,8 +733,26 @@ void dG3DDomain::getDofKeyswithZeroPredictor(std::vector<Dof>& DofVector) const } } } +*/ + // Loop on all elements + std::vector<Dof> R; + const partDomain *dom = static_cast<const partDomain*>(this); + for(groupOfElements::elementContainer::const_iterator it=g->begin(); it!=g->end();++it){ + // get dofs of an element + MElement *ele = *it; + R.clear(); + _space->getKeys(ele,R); + + int nbFF = ele->getNumShapeFunctions(); + + // Loop on all non-local var and shape functions + for(int nl = 0; nl < getNumNonLocalVariable(); ++nl){ + for(int i = 0; i < nbFF; ++i){ + DofVector.push_back(R[i+(3+nl)*nbFF]); + } + } + } } - void dG3DDomain::computeAllIPStrain(AllIPState *aips,const unknownField *ufield,const IPStateBase::whichState ws, const bool virt){ -- GitLab