diff --git a/NonLinearSolver/materialLaw/mlawGenericCrackPhaseField.cpp b/NonLinearSolver/materialLaw/mlawGenericCrackPhaseField.cpp index 6d4fe30927d300d30b7b8a91d18e47ea3c6a03ab..e68411c25e75eb634259e0215e826999756165a2 100644 --- a/NonLinearSolver/materialLaw/mlawGenericCrackPhaseField.cpp +++ b/NonLinearSolver/materialLaw/mlawGenericCrackPhaseField.cpp @@ -111,7 +111,7 @@ mlawGenericCrackPhaseField& mlawGenericCrackPhaseField::operator=(const material _rho = src->_rho; // _lc = src->_lc; _gc = src->_gc; - _test_factor = src->_test_factor; + // _test_factor = src->_test_factor; } @@ -259,8 +259,8 @@ void mlawGenericCrackPhaseField::constitutive( IPGenericCrackPhaseField *qPF1 = dynamic_cast<IPGenericCrackPhaseField *> (q1); // get length scale from Prescribed CLengthLaw - _cLLaw->computeCL(qPF0.getConstRefToLocalDamageVariable(), qPF1.getRefToIPCLength()); - characteristicLength[0]=qPF1.getRefToIPCLength().getConstRefToCL(); + _cLLaw->computeCL(qPF0->getConstRefToLocalDamageVariable(), qPF1->getRefToIPCLength()); + characteristicLength[0]=qPF1->getRefToIPCLength().getConstRefToCL(); double _lc = characteristicLength[0](0,0); @@ -318,37 +318,33 @@ void mlawGenericCrackPhaseField::constitutive( // Compute total P STensorOperation::scale(Pmeca, degrade, P); - // Compute total Tangent - STensorOperation::scale(TangentMeca, degrade, Tangent); - if(elasticTangent !=NULL){ - STensorOperation::scale(elasticTangentMeca, degrade, *elasticTangent); - } + // Update the local damage variable // qPF1.setlocalDamageVariable(((-2.0 * _lc) / _gc) * (1.0 - nonLocalVariable(0)) * hh); - qPF1.getRefToLocalDamageVariable()=(((2.0 * _lc) / _gc) * (1.0 - nonLocalVariable(0)) * hh); - localVariable(0)=qPF1.getConstRefToLocalDamageVariable(); + qPF1->getRefToLocalDamageVariable()=(((2.0 * _lc) / _gc) * (1.0 - nonLocalVariable(0)) * qPF1->getRefToeffEner()); + localVariable(0)=qPF1->getConstRefToLocalDamageVariable(); if(getNbNonLocalVariables() == 1){ - double factor = ((test_factor * 2.0 * _lc) / _gc) * (1.0 - qPF1->getConstRefToNonLocalDamageVariable()); - STensorOperation::scale(qPF1->getConstRefToPMax(), factor, dLocalPlasticStrainDStrain[0]); + // double factor = ((test_factor * 2.0 * _lc) / _gc) * (1.0 - qPF1->getConstRefToNonLocalDamageVariable()); + // STensorOperation::scale(qPF1->getConstRefToPMax(), factor, dLocalPlasticStrainDStrain[0]); double factor = ((2.0 * _lc) / _gc) * (1.0 - nonLocalVariable(0)); - STensorOperation::scale(PMax, factor, dLocalPlasticStrainDStrain[0]); + STensorOperation::scale(qPF1->getConstRefToPMax(), factor, dLocalPlasticStrainDStrain[0]); double factor1 = 2.0 * (1.0 - nonLocalVariable(0)); STensorOperation::scale(Pmeca, factor1, dStressDNonLocalPlasticStrain[0]); - dLocalPlasticStrainDNonLocalPlasticStrain(0,0) = ((test_factor * -2.0 * _lc) / _gc) * qPF1->getConstRefToeffEner(); + dLocalPlasticStrainDNonLocalPlasticStrain(0,0) = ((-2.0 * _lc) / _gc) * qPF1->getConstRefToeffEner(); } qPF1->getRefToElasticEnergy() = qPF1->getConstRefToIpMeca()->defoEnergy() * degrade; double dPlasticDisp = qPF1->getConstRefToIpMeca()->plasticEnergy() - qPF0->getConstRefToIpMeca()->plasticEnergy(); qPF1->setplasticEnergy(qPF0->plasticEnergy() + dPlasticDisp* degrade); - double degrade0 = (1.0 - qPF0->getConstRefToNonLocalDamageVariable()) * - (1.0 - qPF0->getConstRefToNonLocalDamageVariable()); + double degrade0 = (1.0 - nonLocalVariable0(0)) * + (1.0 - nonLocalVariable0(0)); double delta_damageEnergy = qPF1->getConstRefToIpMeca()->defoEnergy() * (degrade-degrade0) - qPF0->damageEnergy(); if(delta_damageEnergy > 0.0){ @@ -396,6 +392,8 @@ void mlawGenericCrackPhaseField::constitutive( if (stiff) { + // Compute total Tangent + STensorOperation::scale(TangentMeca, degrade, Tangent); STensor3& DIrrevEnergyDF = qPF1->getRefToDIrreversibleEnergyDF(); double& DIrrevEnergyDNonlocalVar = qPF1->getRefToDIrreversibleEnergyDNonLocalVariable(); @@ -416,13 +414,13 @@ void mlawGenericCrackPhaseField::constitutive( // // DIrrevEnergyDNonlocalVar = -effEner*q1->getDDamageDp(); // DIrrevEnergyDNonlocalVar = qMecan.defoEnergy() * -2 * (1 - nonLocalVariable(0)); // } - DIrrevEnergyDNonlocalVar = qMecan.defoEnergy() * -2 * (1 - nonLocalVariable(0)); + DIrrevEnergyDNonlocalVar = qPF1->getConstRefToIpMeca()->defoEnergy() * -2 * (1 - nonLocalVariable(0)); } else if (this->getMacroSolver()->getPathFollowingLocalIncrementType() == pathFollowingManager::DISSIPATION_ENERGY) { // DIrrevEnergyDF *= (D - q0->getDamage()); - DIrrevEnergyDF = (1.0 - degrade) * PMax; + DIrrevEnergyDF = (1.0 - degrade) * qPF1->getConstRefToPMax(); // if (qPF1->getNonLocalToLocal()){ @@ -435,13 +433,13 @@ void mlawGenericCrackPhaseField::constitutive( // DIrrevEnergyDNonlocalVar = hh * 2.0 * (1.0 - nonLocalVariable(0)); // } - DIrrevEnergyDNonlocalVar = hh * 2.0 * (1.0 - nonLocalVariable(0)); + DIrrevEnergyDNonlocalVar = qPF1->getRefToeffEner()* 2.0 * (1.0 - nonLocalVariable(0)); } else if (this->getMacroSolver()->getPathFollowingLocalIncrementType() == pathFollowingManager::DAMAGE_ENERGY) { DIrrevEnergyDF = qPF1->getConstRefToPMax() * (1.0 - degrade); - DIrrevEnergyDNonlocalVar = qPF1->getConstRefToeffEner() * 2.0 * (1.0 - qPF1->getConstRefToNonLocalDamageVariable()); + DIrrevEnergyDNonlocalVar = qPF1->getConstRefToeffEner() * 2.0 * (1.0 - nonLocalVariable(0)); } else { diff --git a/NonLinearSolver/materialLaw/mlawGenericCrackPhaseField.h b/NonLinearSolver/materialLaw/mlawGenericCrackPhaseField.h index 7fb4aad5b6b929f92b2fddc091d74028b66749c2..27b8c182758b00fb2ed2b95dd37fe65403fc38b3 100644 --- a/NonLinearSolver/materialLaw/mlawGenericCrackPhaseField.h +++ b/NonLinearSolver/materialLaw/mlawGenericCrackPhaseField.h @@ -26,11 +26,11 @@ protected: double _rho; // double _lc; double _gc; - double _test_factor; + // double _test_factor; public: - mlawGenericCrackPhaseField(int num, double rho, const CLengthLaw &cLLaw, double gc, double test_factor=1.0); + mlawGenericCrackPhaseField(int num, double rho, const CLengthLaw &cLLaw, double gc); const materialLaw& getConstRefMechanicalMaterialLaw() const; materialLaw &getRefMechanicalMaterialLaw(); virtual void setMechanicalMaterialLaw(const materialLaw *matlaw); diff --git a/dG3D/benchmarks/CrackPhaseField/pf.py b/dG3D/benchmarks/CrackPhaseField/pf.py index d78dcdf0714a321d84365d04cac7f9b24bb2da9b..c2b8541d74c491c61b27792371bb88b7d810d03f 100644 --- a/dG3D/benchmarks/CrackPhaseField/pf.py +++ b/dG3D/benchmarks/CrackPhaseField/pf.py @@ -11,7 +11,7 @@ from dG3DpyDebug import* # =================================================================================== Flag_Fracture = 1 # Activate cracks and fracture laws ? (1=yes, 0=no) Solving_mode = 1 # QStatic (=1) vs. Implicit (=2) -LoadingMode = 2 +LoadingMode = 0 # Material law creation @@ -66,7 +66,7 @@ law1.setUseBarF(False) # FracLaw1 = FractureByCohesive3DLaw(FracLawNum1,BulkLawNum1,InterfaceLawNum1) -mylaw1 = GenericCrackPhaseFieldDG3DMaterialLaw(lawnum3, rho, Cl_Law1, 2700.0, 1.0) +mylaw1 = GenericCrackPhaseFieldDG3DMaterialLaw(lawnum3, rho, Cl_Law1, 2700.0) mylaw1.setMechanicalMaterialLaw(law1) # Solver parameters @@ -75,7 +75,7 @@ sol = 2 # Library for solving: Gmm=0 (default) Taucs=1 PETsc=2 if Solving_mode == 1: soltype = 1 # StaticLinear=0 (default, StaticNonLinear=1, Explicit=2, # Multi=3, Implicit=4, Eigen=5) - nstep = 1000 # Number of step + nstep = 100 # Number of step ftime =1.0 # Final time tol=1.e-5 # Relative tolerance for NR scheme tolAbs = 1.e-6 # Absolute tolerance for NR scheme @@ -207,7 +207,7 @@ if LoadingMode == 0: #mysolver.displacementBC("Edge",104,0,0.) # face y = 0 - method = 0 + method = 1 mysolver.pathFollowing(True,method) diff --git a/dG3D/benchmarks/CrackPhaseField1D_da/damage.py b/dG3D/benchmarks/CrackPhaseField1D_da/damage.py index 42843d51761cb933e4cf10a0b3e543a5ffe56263..c71c536f8e86019401be1de70c464af7bdd26ecf 100644 --- a/dG3D/benchmarks/CrackPhaseField1D_da/damage.py +++ b/dG3D/benchmarks/CrackPhaseField1D_da/damage.py @@ -27,7 +27,7 @@ mech_law.setUseBarF(False) if mode == "d": cl_law = IsotropicCLengthLaw(cl_law_num, cl) - pf_law = GenericCrackPhaseFieldDG3DMaterialLaw(pf_law_num, rho, cl_law, gc, 1.0) + pf_law = GenericCrackPhaseFieldDG3DMaterialLaw(pf_law_num, rho, cl_law, gc) pf_law.setMechanicalMaterialLaw(mech_law) diff --git a/dG3D/src/nonLocalDamageDG3DMaterialLaw.cpp b/dG3D/src/nonLocalDamageDG3DMaterialLaw.cpp index f9a4c1ba03ae487285ff276a2932e8a1d403c3d4..04add7ee1bdb85f3071fb283da30c3e85485a47c 100644 --- a/dG3D/src/nonLocalDamageDG3DMaterialLaw.cpp +++ b/dG3D/src/nonLocalDamageDG3DMaterialLaw.cpp @@ -6583,10 +6583,10 @@ double NonlocalDamageTorchANNBasedDG3DMaterialLaw::soundSpeed() const // // TO DO WITH PF -GenericCrackPhaseFieldDG3DMaterialLaw::GenericCrackPhaseFieldDG3DMaterialLaw(const int num, const double rho, const CLengthLaw &cLLaw, const double gc, const double test_factor) : +GenericCrackPhaseFieldDG3DMaterialLaw::GenericCrackPhaseFieldDG3DMaterialLaw(const int num, const double rho, const CLengthLaw &cLLaw, const double gc) : dG3DMaterialLaw(num,rho,true) { - _nldlaw = new mlawGenericCrackPhaseField(num, rho, cLLaw, gc, test_factor); + _nldlaw = new mlawGenericCrackPhaseField(num, rho, cLLaw, gc); // double nu = _nldlaw->poissonRatio(); // double mu = _nldlaw->shearModulus(); // double E = 2.*mu*(1.+nu); diff --git a/dG3D/src/nonLocalDamageDG3DMaterialLaw.h b/dG3D/src/nonLocalDamageDG3DMaterialLaw.h index a1e3a3c5a7a681c521d044e5af311ea7fd5f565e..4450c2cab2efa406d1d5f860e24f2bff68641429 100644 --- a/dG3D/src/nonLocalDamageDG3DMaterialLaw.h +++ b/dG3D/src/nonLocalDamageDG3DMaterialLaw.h @@ -1196,7 +1196,7 @@ class GenericCrackPhaseFieldDG3DMaterialLaw : public dG3DMaterialLaw{ mlawGenericCrackPhaseField *_nldlaw; // pointer to allow to choose between LLN style or Gmsh Style. The choice is performed by the constructor (2 constructors with != arguments) public: - GenericCrackPhaseFieldDG3DMaterialLaw(const int num, const double rho, const CLengthLaw &cLLaw, const double gc, const double test_factor=1.0); + GenericCrackPhaseFieldDG3DMaterialLaw(const int num, const double rho, const CLengthLaw &cLLaw, const double gc); // void setNumberOfNonLocalVariables(int nb) {_nldlaw->setNumberOfNonLocalVariables(nb);} // void setInitialDmax_Inc(double d) {_nldlaw->setInitialDmax_Inc(d);};