diff --git a/NonLinearSolver/nlsolver/homogenizedData.cpp b/NonLinearSolver/nlsolver/homogenizedData.cpp
index 34ed31d52cc3db307917ab816cecae8e844857e8..1bcdb430fac8146939d81070e59f7341c8936550 100644
--- a/NonLinearSolver/nlsolver/homogenizedData.cpp
+++ b/NonLinearSolver/nlsolver/homogenizedData.cpp
@@ -143,7 +143,7 @@ deformationEnergy(src.deformationEnergy),plasticEnergy(src.plasticEnergy),irreve
dIrreversibleEnergydF(src.dIrreversibleEnergydF), dIrreversibleEnergydG(src.dIrreversibleEnergydG),
dIrreversibleEnergydgradT(src.dIrreversibleEnergydgradT),dIrreversibleEnergydT(src.dIrreversibleEnergydT),
dIrreversibleEnergydgradV(src.dIrreversibleEnergydgradV),irreversibleEnergyExtraction(src.irreversibleEnergyExtraction),
-cohesiveLawExtraction(src.cohesiveLawExtraction), dUdF(src.dUdF), dQdG_El(src.dQdG_El),dQdG_correction(src.dQdG_correction)
+cohesiveLawExtraction(src.cohesiveLawExtraction), dUdF(src.dUdF),dQdG_correction(src.dQdG_correction)
{};
homogenizedData& homogenizedData::operator = (const homogenizedData& src){
P = src.P;
@@ -208,7 +208,6 @@ homogenizedData& homogenizedData::operator = (const homogenizedData& src){
irreversibleEnergyExtraction = src.irreversibleEnergyExtraction;
cohesiveLawExtraction = src.cohesiveLawExtraction;
dUdF = src.dUdF;
- dQdG_El = src.dQdG_El;
dQdG_correction = src.dQdG_correction;
return *this;
};
@@ -324,8 +323,6 @@ void allSecondOrderMechanicsFiles::openFiles(const nonLinearMechSolver* solver,c
filename = solver->getFileSavingPrefix()+"second_tangent.csv";
dQdG_File = Tensor63::createFile(filename);
- filename = solver->getFileSavingPrefix()+"second_ElasticTangent.csv";
- dQdG_El_File = Tensor63::createFile(filename);
}
};
@@ -345,9 +342,6 @@ void allSecondOrderMechanicsFiles::dataToFile(const double time, const nonLinear
if (dQdG_File!=NULL){
Tensor63::writeData(dQdG_File,data->getHomogenizedTangentOperator_G_G(),time);
}
- if (dQdG_El_File!=NULL){
- Tensor63::writeData(dQdG_El_File,data->getHomogenizedElasticOperator_G_G(),time);
- }
};
void allConExtraDofFiles::openFiles(const nonLinearMechSolver* solver,const bool strainToFile, const bool stressToFile, const bool failureToFile){
diff --git a/NonLinearSolver/nlsolver/homogenizedData.h b/NonLinearSolver/nlsolver/homogenizedData.h
index 7d3acfe5cfe29746be05e2f0d09fefc3b08e0f5d..c0a7137165afb9e647ae4adb02d3ed46ada08228 100644
--- a/NonLinearSolver/nlsolver/homogenizedData.h
+++ b/NonLinearSolver/nlsolver/homogenizedData.h
@@ -58,7 +58,6 @@ class homogenizedData {
STensor53 dQdF;
STensor63 dQdG;
- STensor63 dQdG_El;
STensor63 dQdG_correction;
std::vector<STensor43> dQdgradT;
std::vector<STensor33> dQdT;
@@ -175,8 +174,6 @@ class homogenizedData {
const STensor53& getHomogenizedTangentOperator_G_F() const {return dQdF;};
STensor63& getHomogenizedTangentOperator_G_G() {return dQdG;};
const STensor63& getHomogenizedTangentOperator_G_G() const {return dQdG;};
- STensor63& getHomogenizedElasticOperator_G_G() {return dQdG_El;};
- const STensor63& getHomogenizedElasticOperator_G_G() const {return dQdG_El;};
STensor63& getHomogenizedCorrectionTerm_G_G() {return dQdG_correction;};
const STensor63& getHomogenizedCorrectionTerm_G_G() const {return dQdG_correction;};
diff --git a/NonLinearSolver/nlsolver/nonLinearMechSolver.cpp b/NonLinearSolver/nlsolver/nonLinearMechSolver.cpp
index 5c805e9798e3e13f16ee12c5d40256418e8f4684..2dc54959cba126d78f47aeee1b22a1d3ea4f830b 100644
--- a/NonLinearSolver/nlsolver/nonLinearMechSolver.cpp
+++ b/NonLinearSolver/nlsolver/nonLinearMechSolver.cpp
@@ -12584,7 +12584,7 @@ void nonLinearMechSolver::extractAverageStressByVolumeIntegral(homogenizedData*
}
BmX.print("BmX");
BmXoX.print("bmXoX- BMoJ");*/
- const STensor43& dPdF = homoData->getHomogenizedTangentOperator_F_F();
+
const STensor63& dQdG_correction = homoData->getHomogenizedCorrectionTerm_G_G();
static STensor33 GM;
static STensor33 modTerm_Q;
@@ -12596,7 +12596,9 @@ void nonLinearMechSolver::extractAverageStressByVolumeIntegral(homogenizedData*
}
STensorOperation::STensor63ContractionSTensor33(dQdG_correction, GM, modTerm_Q);
- for (int i=0; i<3; i++){
+ if(usePresentModuli()){
+ const STensor43& dPdF = homoData->getHomogenizedTangentOperator_F_F();
+ for (int i=0; i<3; i++){
for (int j=0; j<3; j++)
for (int k=0; k<3; k++)
for (int l=0; l<3; l++)
@@ -12604,11 +12606,23 @@ void nonLinearMechSolver::extractAverageStressByVolumeIntegral(homogenizedData*
for (int q=0; q<3; q++){
Q(i,j,q) -= 0.5*(dPdF(i,j,k,l)*(GM(k,l,p))*_rveGeoInertia(p,q)+dPdF(i,q,k,l)*(GM(k,l,p))*_rveGeoInertia(p,j));
}
+ }
}
+ else{
+ const STensor43& dPdF = this->getHomogenizationState(IPStateBase::previous)->getHomogenizedTangentOperator_F_F();
+ for (int i=0; i<3; i++){
+ for (int j=0; j<3; j++)
+ for (int k=0; k<3; k++)
+ for (int l=0; l<3; l++)
+ for (int p=0; p<3; p++)
+ for (int q=0; q<3; q++){
+ Q(i,j,q) -= 0.5*(dPdF(i,j,k,l)*(GM(k,l,p))*_rveGeoInertia(p,q)+dPdF(i,q,k,l)*(GM(k,l,p))*_rveGeoInertia(p,j));
+ }
+ }
+ }
Q += modTerm_Q*this->getRVEVolume();
}
- Q *= (volumeRVEInv);
-
+ Q *= (volumeRVEInv);
}
// extract cohesive law
@@ -13051,24 +13065,7 @@ void nonLinearMechSolver::extractAveragePropertiesByInSystemCondensation(homogen
// extract stress
this->extractAverageStress(homoData);
- if(computedUdF()){
- const STensor43& dPdF = homoData->getHomogenizedTangentOperator_F_F();
- for (int i=0; i<3; i++)
- for (int j=0; j<3; j++)
- for (int k=0; k<3; k++)
- for (int l=0; l<3; l++)
- for (int p=0; p<3; p++)
- for (int q=0; q<3; q++){
- homoData->getHomogenizedElasticOperator_G_G()(i,j,k,l,p,q) = 0.25*(dPdF(i,j,l,p)*_rveGeoInertia(q,k)+dPdF(i,k,l,p)*_rveGeoInertia(q,j))/vRVE
- +0.25*(dPdF(i,j,l,q)*_rveGeoInertia(p,k)+dPdF(i,k,l,q)*_rveGeoInertia(p,j))/vRVE;
- }
- }
-
-
-
-
-
-
+
// extract cohesive law
if (withEnergyDissipation() and getDamageToCohesiveJumpFlag()){
STensor33& DcjumpDF = homoData->getHomogenizedTangentOperator_CohesiveJump_F();
@@ -15256,9 +15253,9 @@ double nonLinearMechSolver::solveMicroSolverSNL(){
nonsys->resetUnknownsToPreviousTimeStep();
_ipf->copy(IPStateBase::previous,IPStateBase::current);
- if (_pathFollowing and computedUdF()){
+ /*if (_pathFollowing and computedUdF()){
pathSys->resetScaleParameter(_timeManager->getTimeStepFactorReduction());
- }
+ }*/
this->getHomogenizationState(IPStateBase::current)->operator=(*this->getHomogenizationState(IPStateBase::previous));
@@ -15668,10 +15665,6 @@ int nonLinearMechSolver::microNewtonRaphsonPathFollowing(const int numstep){
if(computedUdF()){
_ipf->setDeformationGradientGradient(this->getMicroBC()->getSecondOrderKinematicalVariable(), IPStateBase::current);
}
- // (this->getMicroBC()->getSecondOrderKinematicalVariable()).print("KinematicalVariabl G");
-
- // (this->getMicroBC()->getCurrentDeformationGradientGradient()).print("currentVariabl G");
-
_ipf->compute1state(IPStateBase::current,true);
diff --git a/NonLinearSolver/pathFollowing/pbcPathFollowingSystemPETSc.h b/NonLinearSolver/pathFollowing/pbcPathFollowingSystemPETSc.h
index 136d274a754a432f8279da14c2fc16b18a54970c..a3e266cf5ed63842c3eea8610357aa8aaf5ddbb4 100644
--- a/NonLinearSolver/pathFollowing/pbcPathFollowingSystemPETSc.h
+++ b/NonLinearSolver/pathFollowing/pbcPathFollowingSystemPETSc.h
@@ -42,6 +42,7 @@ class pbcPathFollowingSystemPETSC : public pbcNonLinearSystemPETSc<scalar>,
bool _isAllocatedBodyForceVector;
double _loadReduceRatio =1.0;
+
public:
pbcPathFollowingSystemPETSC(MPI_Comm com = PETSC_COMM_SELF) : pbcNonLinearSystemPETSc<scalar>(com),pathFollowingSystem<scalar>(),
_controlParameter(0.),_controlParameterPrev(0.),_controlStep(0), _pseudoTimeIncrement(0.), _pseudoTimeIncrement_prev(0.),
@@ -200,7 +201,8 @@ class pbcPathFollowingSystemPETSC : public pbcNonLinearSystemPETSc<scalar>,
_controlParameterPrev = 0.;
_controlStep = 0.;
_setScale = false;
-
+ _loadReduceRatio = 1.0;
+
_try(VecAssemblyBegin(_stateStep));
_try(VecAssemblyEnd(_stateStep));
_try(VecZeroEntries(_stateStep));
@@ -246,7 +248,7 @@ class pbcPathFollowingSystemPETSC : public pbcNonLinearSystemPETSc<scalar>,
_try(VecDestroy(&Rq));
_try(VecDestroy(&KRq));
- if (_iter == 0){
+ if (_iter == 0) {
// solve linear system
_try(KSPSolve(this->_ksp, this->_qeff, this->_x));
PetscScalar vTv;
@@ -255,32 +257,31 @@ class pbcPathFollowingSystemPETSC : public pbcNonLinearSystemPETSc<scalar>,
_uScale = sqrt(vTv);
uScale2 = vTv;
_setScale = true;
- // Msg::Error("Uscal2 = %f",uScale2);
}
- PetscScalar A = _aU*vTv/uScale2+_aL;
+ PetscScalar A0 = _aU*vTv/uScale2+_aL;
- double a1 = _loadReduceRatio*_pseudoTimeIncrement/sqrt(A);
+ double a1 = _loadReduceRatio*_pseudoTimeIncrement/sqrt(A0);
double a2 = -1.*a1;
PetscScalar uprevTu;
_try(VecDot(_stateStepPrevious,this->_x,&uprevTu));
- //Msg::Error("cre = %f",uprevTu);
+
if (uprevTu>= 0){
- _controlStep = a1;
- _controlParameter += a1/_loadReduceRatio;
+ _controlStep = a1;
+ _controlParameter += a1/_loadReduceRatio;
_try(VecAXPY(_stateStep,a1,this->_x));
_try(VecAXPY(this->_xsol,a1,this->_x));
}
else{
_controlStep = a2;
- _controlParameter += a2/_loadReduceRatio;
+ _controlParameter += a2/_loadReduceRatio;
_try(VecAXPY(_stateStep,a2,this->_x));
_try(VecAXPY(this->_xsol,a2,this->_x));
}
-
}
else{
+
/** compute right hand side**/
if (!this->_rhsflag){
Vec Rrc; // Rrc = R*rc;
@@ -354,8 +355,6 @@ class pbcPathFollowingSystemPETSC : public pbcNonLinearSystemPETSc<scalar>,
if (_beta>1) _beta = 1.;
}
-
-
/** du += dr**/
_try(VecAXPY(_stateStep,_beta,dr));
@@ -413,7 +412,6 @@ class pbcPathFollowingSystemPETSC : public pbcNonLinearSystemPETSc<scalar>,
Msg::Error("relaxation error delta = %e",delta);
return 0;
}
-
_try(VecDestroy(&dr));
}
_iter++;
diff --git a/dG3D/src/dG3DDomain.cpp b/dG3D/src/dG3DDomain.cpp
index bcb0fa64a6d0ded4c8978604c34e25d4b45f8c88..769042e90a7d772660e392dcea23738b15d2a1f4 100644
--- a/dG3D/src/dG3DDomain.cpp
+++ b/dG3D/src/dG3DDomain.cpp
@@ -2131,12 +2131,12 @@ void dG3DDomain::computeIpv(AllIPState *aips,MElement *e, IPStateBase::whichStat
bool usePreviousModuli = false;
if (hasBodyForceForHO()){
- if(usePresentModuliHO()){
- needDTangent=true;
- }
- else {
+ if(usePresentModuliHO()){
+ needDTangent=true;
+ }
+ else{
usePreviousModuli = true;
- }
+ }
}
if(_evalStiff) { this->computeStrain(e,npts_bulk,GP,vips,ws,disp,useBarF);};
for(int j=0;j<npts_bulk;j++)
@@ -3246,13 +3246,13 @@ void dG3DDomain::computeIpv(AllIPState *aips,MInterfaceElement *ie, IntPt *GP,co
bool usePreviousModuli = false;
if (hasBodyForceForHO()){
- if(usePresentModuliHO()){
- needDTangent=true;
- }
- else {
- usePreviousModuli = true;
- }
- }
+ if(usePresentModuliHO()){
+ needDTangent=true;
+ }
+ else{
+ usePreviousModuli = true;
+ }
+ }
AllIPState::ipstateElementContainer *vips = aips->getIPstate(ele->getNum());
//Msg::Info("nb interf gp %d",npts);
for(int j=0;j<npts;j++){
@@ -5044,13 +5044,13 @@ void dG3DDomain::computeIPVariableMPI(AllIPState *aips,const unknownField *ufiel
bool usePreviousModuli = false;
if (hasBodyForceForHO()){
- if(usePresentModuliHO()){
- needDTangent=true;
- }
- else {
+ if(usePresentModuliHO()){
+ needDTangent=true;
+ }
+ else{
usePreviousModuli = true;
- }
- }
+ }
+ }
for (std::set<int>::iterator it = _domainIPBulk.begin(); it != _domainIPBulk.end(); it++){
int num = *it;
int elnum, gnum;
diff --git a/dG3D/src/dG3DIPVariable.cpp b/dG3D/src/dG3DIPVariable.cpp
index 4d351e89e74649145223eb4e88b1d637b7173d17..f73b15f8db64a72a9ad07f9e66a154218d9988c6 100644
--- a/dG3D/src/dG3DIPVariable.cpp
+++ b/dG3D/src/dG3DIPVariable.cpp
@@ -1192,15 +1192,15 @@ void dG3DIPVariable::setValueForBodyForce(const STensor43& K, const bool usePres
STensorOperation::zero(_bodyForceForHO->dBmdFm);
STensorOperation::zero(_bodyForceForHO->dPMGdFm);
if(usePresentModuli){
- STensorOperation::multSTensor43STensor33(_bodyForceForHO->dFmdFM, _bodyForceForHO->GM, tmp1);
- STensorOperation::multSTensor43STensor33RightTranspose(_bodyForceForHO->dFmdFM, _bodyForceForHO->GM, tmp2);
+ STensorOperation::multSTensor43STensor33(_bodyForceForHO->dFmdFM, _bodyForceForHO->GM, tmp1);
+ STensorOperation::multSTensor43STensor33RightTranspose(_bodyForceForHO->dFmdFM, _bodyForceForHO->GM, tmp2);
for (int i=0; i<3; i++){
for (int j=0; j<3; j++){
for (int k=0; k<3; k++){
for (int l=0; l<3; l++){
for (int p=0; p<3; p++){
- for (int q=0; q<3; q++){
- _bodyForceForHO->dBmdFm(i,p,q) -= _bodyForceForHO->dCmdFm(i,j,k,l,p,q)*(tmp1(k,l,j)+tmp2(k,l,j))/2.0;
+ for (int q=0; q<3; q++){
+ _bodyForceForHO->dBmdFm(i,p,q) -= _bodyForceForHO->dCmdFm(i,j,k,l,p,q)*(tmp1(k,l,j)+tmp2(k,l,j))/2.0;
for (int s=0; s<3; s++){
_bodyForceForHO->dPMGdFm(i,j,s,p,q) += _bodyForceForHO->dCmdFm(i,j,k,l,p,q)*(tmp1(k,l,s)+tmp2(k,l,s))/2.0;
}
@@ -1209,8 +1209,8 @@ void dG3DIPVariable::setValueForBodyForce(const STensor43& K, const bool usePres
}
}
}
- }
- }
+ }
+ }
}
else{
Msg::Error("Body Force is not defined");
diff --git a/dG3D/src/hoDGMaterialLaw.cpp b/dG3D/src/hoDGMaterialLaw.cpp
index b250daeeec92d12e119600557f60f93393805a3d..bb30ccd9e7abb4d0bc851c2e66efadd0f5257775 100644
--- a/dG3D/src/hoDGMaterialLaw.cpp
+++ b/dG3D/src/hoDGMaterialLaw.cpp
@@ -386,7 +386,18 @@ void hoDGMultiscaleMaterialLaw::fillElasticStiffness(double & rho, STensor43& LF
JF = solver->getHomogenizationState(IPStateBase::current)->getHomogenizedTangentOperator_G_F();
if (_microSolver->getMicroBC()->getOrder() == 2){
if(solver->computedUdF()){
- JG = solver->getHomogenizationState(IPStateBase::current)->getHomogenizedElasticOperator_G_G();
+ STensorOperation::zero(JG);
+ const STensor3& _rveGeoInertia= solver->getRVEGeometricalInertia();
+ double vRVE = solver->getRVEVolume();
+ for (int i=0; i<3; i++)
+ for (int j=0; j<3; j++)
+ for (int k=0; k<3; k++)
+ for (int l=0; l<3; l++)
+ for (int p=0; p<3; p++)
+ for (int q=0; q<3; q++)
+ JG(i,j,k,l,p,q) = 0.25*(LF(i,j,l,p)*_rveGeoInertia(q,k)+LF(i,k,l,p)*_rveGeoInertia(q,j))/vRVE
+ +0.25*(LF(i,j,l,q)*_rveGeoInertia(p,k)+LF(i,k,l,q)*_rveGeoInertia(p,j))/vRVE;
+
}
else{
JG = solver->getHomogenizationState(IPStateBase::current)->getHomogenizedTangentOperator_G_G();