diff --git a/NonLinearSolver/materialLaw/mlawTFAMaterialLaws.cpp b/NonLinearSolver/materialLaw/mlawTFAMaterialLaws.cpp
index bda7ad45b8094ce584ecdb6185b3ecf4ab93a4cf..4669ea9bde6c02ecfc9b90f54f553cf023742803 100644
--- a/NonLinearSolver/materialLaw/mlawTFAMaterialLaws.cpp
+++ b/NonLinearSolver/materialLaw/mlawTFAMaterialLaws.cpp
@@ -285,8 +285,14 @@ void mlawGenericTFAMaterialLaws::setTFAMethod(int TFAMethodNum, int dim, int cor
else if(TFAMethodNum == 2)
{
polarization = true;
- Msg::Info("IncrementalTangentTFAPolarization is used !!!");
- _reductionModel = new IncrementalTangentTFAPolarization(solverType, dim, correction, withFrame, polarization, tag);
+ Msg::Info("IncrementalTangentHS is used !!!");
+ _reductionModel = new IncrementalTangentHS(solverType, dim, correction, withFrame, polarization, tag);
+ }
+ else if(TFAMethodNum == 3)
+ {
+ polarization = true;
+ Msg::Info("IncrementalTangentPFA is used !!!");
+ _reductionModel = new IncrementalTangentPFA(solverType, dim, correction, withFrame, polarization, tag);
}
else
{
diff --git a/NonLinearSolver/modelReduction/modelReduction.cpp b/NonLinearSolver/modelReduction/modelReduction.cpp
index 65392c38ea8f599cc7d4825f9849ed5ecaba9791..15d5f97aa03e8ec0aafdb1aeeb0227430c3f16ff 100644
--- a/NonLinearSolver/modelReduction/modelReduction.cpp
+++ b/NonLinearSolver/modelReduction/modelReduction.cpp
@@ -612,7 +612,7 @@ void ReductionTFA::preallocateSystem(const ClusteringData *q0, const Clustering
}
else if(_polarization)
{
- getBiLinearTermPolarization(icl,q0,q1, othersId, m);
+ getBiLinearTermHS(icl,q0,q1, othersId, m);
}
else
{
@@ -868,7 +868,7 @@ void ReductionTFA::assembleRes(const ClusteringData *q0, const ClusteringData *
{
for (int icl = 0; icl < nbClusters; icl++)
{
- getLinearTermPolarization(icl,q0,q1,r);
+ getLinearTermHS(icl,q0,q1,r);
std::vector<Dof> R;
getKeys(icl,R);
for (int j=0; j<R.size(); j++)
@@ -944,7 +944,7 @@ void ReductionTFA::assembleDresDu(const ClusteringData *q0, const ClusteringDat
{
for (int icl = 0; icl < nbClusters; icl++)
{
- getBiLinearTermPolarization(icl,q0,q1,othersId,Ke);
+ getBiLinearTermHS(icl,q0,q1,othersId,Ke);
//Ke.print("Ke");
//
std::vector<Dof> R, C;
@@ -1496,7 +1496,125 @@ void ReductionTFA::getTangentTerm(int clusterId, const ClusteringData *q0, cons
};
-void ReductionTFA::getLinearTermPolarization(int clusterId, const ClusteringData *q0, const ClusteringData *q1, fullVector<double>& res) const
+void ReductionTFA::getLinearTermHS(int clusterId, const ClusteringData *q0, const ClusteringData *q1, fullVector<double>& res) const
+{
+ int nbClusters = _clusteringData->getTotalNumberOfClusters();
+
+ const STensor3& eps1i = q1->getConstRefToStrain(clusterId);
+ const STensor3& eps0i = q0->getConstRefToStrain(clusterId);
+
+ const STensor3& macroeps1 = q1->getConstRefToMacroStrain();
+ const STensor3& macroeps0 = q0->getConstRefToMacroStrain();
+
+ static STensor3 resTen;
+ resTen = eps1i;
+ resTen -= eps0i;
+ resTen -= macroeps1;
+ resTen += macroeps0;
+
+ const STensor43& C0iso = _clusteringData->getReferenceStiffnessIso();
+ double G0 = C0iso(0,1,0,1);
+ const double& Gtan = q1->getConstRefToTangentShearModulus();
+
+ for (int jcl =0; jcl < nbClusters; jcl++)
+ {
+ const STensor43 &Dij = getClusterInteractionTensor(clusterId,jcl);
+
+ const STensor3& sig1j = q1->getConstRefToStress(jcl);
+ const STensor3& eps1j = q1->getConstRefToStrain(jcl);
+ const STensor3& sig0j = q0->getConstRefToStress(jcl);
+ const STensor3& eps0j = q0->getConstRefToStrain(jcl);
+
+ STensor3 dsigj(sig1j);
+ dsigj -= sig0j;
+ STensor3 depsj(eps1j);
+ depsj -= eps0j;
+
+ STensor3 fac;
+ fac = dsigj;
+ fac *= G0/Gtan;
+ STensorOperation::multSTensor43STensor3Add(C0iso,depsj,-1.,fac);
+ STensorOperation::multSTensor43STensor3Add(Dij,fac,-1.,resTen);
+ }
+ STensorOperation::fromSTensor3ToFullVector(resTen,res);
+};
+void ReductionTFA::getBiLinearTermHS(int clusterId, const ClusteringData *q0, const ClusteringData *q1, std::vector<int>& othersId, fullMatrix<double>& dres) const
+{
+ int nbClusters = _clusteringData->getTotalNumberOfClusters();
+ othersId.resize(nbClusters);
+ for (int j=0; j< nbClusters; j++)
+ {
+ othersId[j] = j;
+ }
+ dres.resize(6,6*nbClusters,true);
+
+ double vtot;
+ vtot = totalVolume();
+
+ for (int jcl=0; jcl< nbClusters; jcl++)
+ {
+ static STensor43 J_ij;
+ getJacHS_rs(clusterId,jcl,q0,q1,J_ij);
+ static fullMatrix<double> J_ij_matrix;
+ STensorOperation::fromSTensor43ToFullMatrix(J_ij,J_ij_matrix);
+ for (int row=0; row<6; row++)
+ {
+ for (int col=0; col<6; col++)
+ {
+ dres(row, col+6*jcl) += J_ij_matrix(row,col);
+ }
+ }
+ }
+}
+void ReductionTFA::getJacHS_rs(int r, int s, const ClusteringData *q0, const ClusteringData *q1, STensor43& Jac_rs) const
+{
+ if (r==s)
+ {
+ STensorOperation::unity(Jac_rs);
+ }
+ else
+ {
+ STensorOperation::zero(Jac_rs);
+ }
+
+ const STensor3& macroeps1 = q1->getConstRefToMacroStrain();
+
+ const STensor43& C0iso = _clusteringData->getReferenceStiffnessIso();
+ double G0 = C0iso(0,1,0,1);
+ const double& Gtan = q1->getConstRefToTangentShearModulus();
+ const STensor43 &Drs = getClusterInteractionTensor(r,s);
+ const STensor43 &Calgs = q1->getConstRefToTangent(s);
+
+ const STensor3& sig1s = q1->getConstRefToStress(s);
+ const STensor3& sig0s = q0->getConstRefToStress(s);
+ STensor3 dsigs(sig1s);
+ dsigs -= sig0s;
+
+ STensor43 fac;
+ fac = Calgs;
+ fac *= G0/Gtan;
+ fac -= C0iso;
+ STensorOperation::multSTensor43Add(Drs,fac,-1.,Jac_rs);
+
+ /*STensor3 dGtandeps;
+ dtangentShearModulusOveralldeps(q0,q1,s,dGtandeps);
+ STensor3 fac2_part;
+ int nbClusters = _clusteringData->getTotalNumberOfClusters();
+ for (int p =0; p < nbClusters; p++)
+ {
+ const STensor43 &Drp = getClusterInteractionTensor(r,p);
+ const STensor3& sig1p = q1->getConstRefToStress(p);
+ const STensor3& sig0p = q0->getConstRefToStress(p);
+ STensor3 Dsigp(sig1p);
+ Dsigp -= sig0p;
+ STensorOperation::multSTensor43STensor3Add(Drp,Dsigp,1.,fac2_part);
+ }
+ fac2_part *= G0/(Gtan*Gtan);
+ STensorOperation::prodAdd(fac2_part,dGtandeps,1.,Jac_rs);*/
+}
+
+
+void ReductionTFA::getLinearTermPFA(int clusterId, const ClusteringData *q0, const ClusteringData *q1, fullVector<double>& res) const
{
int nbClusters = _clusteringData->getTotalNumberOfClusters();
@@ -1560,17 +1678,15 @@ void ReductionTFA::getLinearTermPolarization(int clusterId, const ClusteringData
//STensorOperation::multSTensor43STensor3(prefac,dsigj,fac);
STensorOperation::multSTensor43STensor3Add(C0iso,depsj,-1.,fac);
-
//STensorOperation::multSTensor43(C0,C_inv,prefac);
//STensorOperation::multSTensor43STensor3(prefac,dsigj,fac);
//STensorOperation::multSTensor43STensor3Add(C0,depsj,-1.,fac);
-
-
+
STensorOperation::multSTensor43STensor3Add(Dij,fac,-1.,resTen);
}
STensorOperation::fromSTensor3ToFullVector(resTen,res);
};
-void ReductionTFA::getBiLinearTermPolarization(int clusterId, const ClusteringData *q0, const ClusteringData *q1, std::vector<int>& othersId, fullMatrix<double>& dres) const
+void ReductionTFA::getBiLinearTermPFA(int clusterId, const ClusteringData *q0, const ClusteringData *q1, std::vector<int>& othersId, fullMatrix<double>& dres) const
{
int nbClusters = _clusteringData->getTotalNumberOfClusters();
othersId.resize(nbClusters);
@@ -1586,7 +1702,7 @@ void ReductionTFA::getBiLinearTermPolarization(int clusterId, const ClusteringDa
for (int jcl=0; jcl< nbClusters; jcl++)
{
static STensor43 J_ij;
- getJacPolarization_rs(clusterId,jcl,q0,q1,J_ij);
+ getJacHS_rs(clusterId,jcl,q0,q1,J_ij);
static fullMatrix<double> J_ij_matrix;
STensorOperation::fromSTensor43ToFullMatrix(J_ij,J_ij_matrix);
for (int row=0; row<6; row++)
@@ -1598,7 +1714,7 @@ void ReductionTFA::getBiLinearTermPolarization(int clusterId, const ClusteringDa
}
}
}
-void ReductionTFA::getJacPolarization_rs(int r, int s, const ClusteringData *q0, const ClusteringData *q1, STensor43& Jac_rs) const
+void ReductionTFA::getJacPFA_rs(int r, int s, const ClusteringData *q0, const ClusteringData *q1, STensor43& Jac_rs) const
{
if (r==s)
{
@@ -1674,6 +1790,7 @@ void ReductionTFA::getJacPolarization_rs(int r, int s, const ClusteringData *q0,
STensorOperation::prodAdd(fac2_part,dGtandeps,1.,Jac_rs);*/
}
+
void ReductionTFA::elasticStiffnessHom(STensor43& CelHom) const
{
int nbClusters = _clusteringData->getTotalNumberOfClusters();
@@ -2227,20 +2344,18 @@ void IncrementalTangentTFAWithFrame::localConstitutive(const STensor3& F0, const
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-IncrementalTangentTFAPolarization::IncrementalTangentTFAPolarization(int solverType, int dim, int correction, bool withFrame, bool polarization, int tag): ReductionTFA(solverType, dim, correction, withFrame, polarization, tag){}
-IncrementalTangentTFAPolarization::~IncrementalTangentTFAPolarization(){}
+IncrementalTangentHS::IncrementalTangentHS(int solverType, int dim, int correction, bool withFrame, bool polarization, int tag): ReductionTFA(solverType, dim, correction, withFrame, polarization, tag){}
+IncrementalTangentHS::~IncrementalTangentHS(){}
-void IncrementalTangentTFAPolarization::localConstitutive(const STensor3& F0, const STensor3& Fn, const ClusteringData *q0, ClusteringData *q1)
+void IncrementalTangentHS::localConstitutive(const STensor3& F0, const STensor3& Fn, const ClusteringData *q0, ClusteringData *q1)
{
int nbClusters = _clusteringData->getTotalNumberOfClusters();
const STensor43& C0 = _clusteringData->getReferenceStiffness();
const STensor43& C0iso = _clusteringData->getReferenceStiffnessIso();
double G0 = C0iso(0,1,0,1);
- double k0 = C0iso(0,0,0,0) - 4./3.*G0;
- //
const STensor3& eps0 = q0->getConstRefToMacroStrain();
STensor3& eps1 = q1->getRefToMacroStrain();
for (int i=0; i<3; i++)
@@ -2260,15 +2375,88 @@ void IncrementalTangentTFAPolarization::localConstitutive(const STensor3& F0, co
double DepsEq = TFANumericTools::EquivalentStrain(Deps);
const STensor3& P0hom = q0->getConstRefToMacroStress();
- const STensor3& epsEigMacro0 = q0->getConstRefToMacroEigenStrain();
STensor3& Phom = q1->getRefToMacroStress();
- STensor3& epsEigMacro = q1->getRefToMacroEigenStrain();
+
+ // local constitutive behavior
+ static STensor3 Isecond(1.);
+ for (int icl =0; icl < nbClusters; icl++)
+ {
+ const materialLaw *law = getLawForCluster(icl);
+ static STensor3 F0i, F1i;
+ F0i = q0->getConstRefToStrain(icl);
+ F1i = q1->getConstRefToStrain(icl);
+ F0i += Isecond;
+ F1i += Isecond;
+ STensor3& P1i = q1->getRefToStress(icl);
+ STensor43& Tangenti = q1->getRefToTangent(icl);
+ STensor43& depsEig1depsi = q1->getRefTodEigenStraindStrain(icl);
+ STensor3& epsEig1i = q1->getRefToEigenStrain(icl);
+ //
+ const IPVariable* ipv0i = q0->getConstRefToIPv(icl);
+ IPVariable* ipv1i = q1->getRefToIPv(icl);
+ if(_clusterIncOri)
+ {
+ const fullVector<double>& euler_i = getClusterMeanFiberOrientation(icl);
+ law->constitutiveTFA_incOri(F0i,F1i,euler_i,P1i,ipv0i,ipv1i,Tangenti,epsEig1i,depsEig1depsi);
+ }
+ else
+ {
+ law->constitutiveTFA(F0i,F1i,P1i,ipv0i,ipv1i,Tangenti,epsEig1i,depsEig1depsi);
+ }
+ }
- STensor3 PtrHom;
- elasticTotalTrialStressHom(q1,PtrHom);
+ STensor3 DPhom;
+ double DsigEqHom;
+ q1->computeHomogenizedStress(Phom);
+ DPhom = Phom;
+ DPhom -= P0hom;
+ DsigEqHom = TFANumericTools::EquivalentStress(DPhom);
- STensor3 PtrHom_true;
- STensorOperation::multSTensor43STensor3(C0,eps1,PtrHom_true);
+ double& Gel = q1->getRefToElasticShearModulus();
+ Gel = G0;
+ double& Gtan = q1->getRefToTangentShearModulus();
+ Gtan = G0;
+ if(DepsEq>1.e-10)
+ {
+ Gtan = DsigEqHom/(3.*DepsEq);
+ }
+ if(Gtan < 5.){Gtan=5.;};
+};
+
+IncrementalTangentPFA::IncrementalTangentPFA(int solverType, int dim, int correction, bool withFrame, bool polarization, int tag): ReductionTFA(solverType, dim, correction, withFrame, polarization, tag){}
+IncrementalTangentPFA::~IncrementalTangentPFA(){}
+
+
+void IncrementalTangentPFA::localConstitutive(const STensor3& F0, const STensor3& Fn, const ClusteringData *q0, ClusteringData *q1)
+{
+ int nbClusters = _clusteringData->getTotalNumberOfClusters();
+
+ const STensor43& C0 = _clusteringData->getReferenceStiffness();
+ const STensor43& C0iso = _clusteringData->getReferenceStiffnessIso();
+ double G0 = C0iso(0,1,0,1);
+ double k0 = C0iso(0,0,0,0) - 4./3.*G0;
+
+ //
+ const STensor3& eps0 = q0->getConstRefToMacroStrain();
+ STensor3& eps1 = q1->getRefToMacroStrain();
+ for (int i=0; i<3; i++)
+ {
+ for (int j=0; j<3; j++)
+ {
+ eps1(i,j) = (Fn(i,j)+Fn(j,i))*0.5;
+ }
+ }
+ eps1(0,0) -= 1.;
+ eps1(1,1) -= 1.;
+ eps1(2,2) -= 1.;
+ double epsEq = TFANumericTools::EquivalentStrain(eps1);
+ STensor3 Deps;
+ Deps = eps1;
+ Deps -= eps0;
+ double DepsEq = TFANumericTools::EquivalentStrain(Deps);
+
+ const STensor3& P0hom = q0->getConstRefToMacroStress();
+ STensor3& Phom = q1->getRefToMacroStress();
// local constitutive behavior
static STensor3 Isecond(1.);
@@ -2286,8 +2474,7 @@ void IncrementalTangentTFAPolarization::localConstitutive(const STensor3& F0, co
STensor3& epsEig1i = q1->getRefToEigenStrain(icl);
//
const IPVariable* ipv0i = q0->getConstRefToIPv(icl);
- IPVariable* ipv1i = q1->getRefToIPv(icl);
- STensor43& Secanti = q1->getRefToSecant(icl);
+ IPVariable* ipv1i = q1->getRefToIPv(icl);
if(_clusterIncOri)
{
const fullVector<double>& euler_i = getClusterMeanFiberOrientation(icl);
@@ -2300,26 +2487,12 @@ void IncrementalTangentTFAPolarization::localConstitutive(const STensor3& F0, co
}
STensor3 DPhom,DepsEigMacro;
- double sigEqHom,DsigEqHom;
- double DepsEigEqHom;
+ double DsigEqHom;
q1->computeHomogenizedStress(Phom);
- //sigEqHom = TFANumericTools::EquivalentStress(Phom);
DPhom = Phom;
DPhom -= P0hom;
DsigEqHom = TFANumericTools::EquivalentStress(DPhom);
-
- //STensor3 sigEigMacro;
- //sigEigMacro = Phom;
- //sigEigMacro -= PtrHom;
-
- //STensor43 SelHom;
- //STensorOperation::inverseSTensor43(C0,SelHom);
- //STensorOperation::multSTensor43STensor3(SelHom,sigEigMacro,epsEigMacro);
- //epsEigMacro *= -1.;
- //DepsEigMacro = epsEigMacro;
- //DepsEigMacro -= epsEigMacro0;
- //DepsEigEqHom = TFANumericTools::EquivalentStrain(DepsEigMacro);
double& Gel = q1->getRefToElasticShearModulus();
Gel = G0;
@@ -2327,8 +2500,6 @@ void IncrementalTangentTFAPolarization::localConstitutive(const STensor3& F0, co
Gtan = G0;
if(DepsEq>1.e-10)
{
- //DepsEigEqHom = DepsEq - DsigEqHom/(3.*G0);
- //Gtan *= (1.-DepsEigEqHom/DepsEq);
Gtan = DsigEqHom/(3.*DepsEq);
}
if(Gtan < 5){Gtan=5;};
@@ -2347,6 +2518,7 @@ void IncrementalTangentTFAPolarization::localConstitutive(const STensor3& F0, co
};
+//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
HierarchicalReductionTFA::HierarchicalReductionTFA(int solverType, int dim, int correction, bool withFrame, bool polarization, int tag):
ReductionTFA(solverType, dim, 0, false, polarization, tag), _sysHierarchical(NULL) {}
@@ -3433,7 +3605,7 @@ void HierarchicalReductionTFA::getJac_rsHierarchical(int r, int s, const Cluster
STensorOperation::multSTensor43Add(Drs,depsEig1dstrains,1.,Jac_rs);
}
-void HierarchicalReductionTFA::getLinearTermPolarizationHierarchical(int clusterId, const ClusteringData *q0, const ClusteringData *q1, fullVector<double>& res) const
+void HierarchicalReductionTFA::getLinearTermHSHierarchical(int clusterId, const ClusteringData *q0, const ClusteringData *q1, fullVector<double>& res) const
{
int nbClusters = q0->clustersNb;
@@ -3482,7 +3654,7 @@ void HierarchicalReductionTFA::getLinearTermPolarizationHierarchical(int cluster
}
STensorOperation::fromSTensor3ToFullVector(resTen,res);
};
-void HierarchicalReductionTFA::getBiLinearTermPolarizationHierarchical(int clusterId, const ClusteringData *q0, const ClusteringData *q1, std::vector<int>& othersId, fullMatrix<double>& dres)
+void HierarchicalReductionTFA::getBiLinearTermHSHierarchical(int clusterId, const ClusteringData *q0, const ClusteringData *q1, std::vector<int>& othersId, fullMatrix<double>& dres)
const
{
int nbClusters = q0->clustersNb;
@@ -3499,7 +3671,7 @@ const
for (int jcl=0; jcl< nbClusters; jcl++)
{
static STensor43 J_ij;
- getJacPolarization_rsHierarchical(clusterId,jcl,q0,q1,J_ij);
+ getJacHS_rsHierarchical(clusterId,jcl,q0,q1,J_ij);
static fullMatrix<double> J_ij_matrix;
STensorOperation::fromSTensor43ToFullMatrix(J_ij,J_ij_matrix);
for (int row=0; row<6; row++)
@@ -3511,7 +3683,7 @@ const
}
}
}
-void HierarchicalReductionTFA::getJacPolarization_rsHierarchical(int r, int s, const ClusteringData *q0, const ClusteringData *q1, STensor43& Jac_rs) const
+void HierarchicalReductionTFA::getJacHS_rsHierarchical(int r, int s, const ClusteringData *q0, const ClusteringData *q1, STensor43& Jac_rs) const
{
std::vector<int> clHostIDs = q0->upperClusterIDs;
diff --git a/NonLinearSolver/modelReduction/modelReduction.h b/NonLinearSolver/modelReduction/modelReduction.h
index d380eef553f2a574f27db19c80ce6db5798760f0..c09f7948ead0f23dd9b52c2d68554c44c18aa3c5 100644
--- a/NonLinearSolver/modelReduction/modelReduction.h
+++ b/NonLinearSolver/modelReduction/modelReduction.h
@@ -23,8 +23,8 @@ class nonLinearSystemPETSc;
class modelReduction
{
public:
- enum ReductionMethod{TFA_IncrementalTangent=0, TFA_IncrementalTangentWithFrame, TFA_IncrementalTangentPolarization,
- TFA_HierarchicalIncrementalTangent, TFA_HierarchicalIncrementalTangentMixed};
+ enum ReductionMethod{TFA_IncrementalTangent=0, TFA_IncrementalTangentWithFrame, HS_IncrementalTangent, PFA_IncrementalTangent,
+ TFA_HierarchicalIncrementalTangent};
#ifndef SWIG
public:
modelReduction(){}
@@ -174,9 +174,13 @@ class ReductionTFA: public modelReduction
//////////////////////////////////////////////////////////////////////////////////////////
virtual void getJac_rs(int r, int s, const ClusteringData *q0, const ClusteringData *q1, STensor43& Jac_rs) const;
- virtual void getLinearTermPolarization(int clusterId, const ClusteringData *q0, const ClusteringData *q1, fullVector<double>& res) const;
- virtual void getBiLinearTermPolarization(int clusterId, const ClusteringData *q0, const ClusteringData *q1, std::vector<int>& othersId, fullMatrix<double>& dres) const;
- virtual void getJacPolarization_rs(int r, int s, const ClusteringData *q0, const ClusteringData *q1, STensor43& Jac_rs) const;
+ virtual void getLinearTermHS(int clusterId, const ClusteringData *q0, const ClusteringData *q1, fullVector<double>& res) const;
+ virtual void getBiLinearTermHS(int clusterId, const ClusteringData *q0, const ClusteringData *q1, std::vector<int>& othersId, fullMatrix<double>& dres) const;
+ virtual void getJacHS_rs(int r, int s, const ClusteringData *q0, const ClusteringData *q1, STensor43& Jac_rs) const;
+
+ virtual void getLinearTermPFA(int clusterId, const ClusteringData *q0, const ClusteringData *q1, fullVector<double>& res) const;
+ virtual void getBiLinearTermPFA(int clusterId, const ClusteringData *q0, const ClusteringData *q1, std::vector<int>& othersId, fullMatrix<double>& dres) const;
+ virtual void getJacPFA_rs(int r, int s, const ClusteringData *q0, const ClusteringData *q1, STensor43& Jac_rs) const;
virtual void elasticStiffnessHom(STensor43& CelHom) const;
virtual void elasticTotalTrialStressHom(const ClusteringData *q, STensor3& PtrHom) const;
@@ -234,13 +238,25 @@ class IncrementalTangentTFAWithFrame : public ReductionTFA
virtual void localConstitutive(const STensor3& F0, const STensor3& Fn, const ClusteringData *q0, ClusteringData *q1);
#endif // SWIG
};
-class IncrementalTangentTFAPolarization : public ReductionTFA
+class IncrementalTangentHS : public ReductionTFA
+{
+ public:
+ IncrementalTangentHS(int solverType, int dim, int correction, bool withFrame, bool polarization, int tag=1000);
+ #ifndef SWIG
+ virtual ~IncrementalTangentHS();
+ virtual ReductionMethod getMethod() const {return modelReduction::HS_IncrementalTangent;};
+
+ protected:
+ virtual void localConstitutive(const STensor3& F0, const STensor3& Fn, const ClusteringData *q0, ClusteringData *q1);
+ #endif // SWIG
+};
+class IncrementalTangentPFA : public ReductionTFA
{
public:
- IncrementalTangentTFAPolarization(int solverType, int dim, int correction, bool withFrame, bool polarization, int tag=1000);
+ IncrementalTangentPFA(int solverType, int dim, int correction, bool withFrame, bool polarization, int tag=1000);
#ifndef SWIG
- virtual ~IncrementalTangentTFAPolarization();
- virtual ReductionMethod getMethod() const {return modelReduction::TFA_IncrementalTangentPolarization;};
+ virtual ~IncrementalTangentPFA();
+ virtual ReductionMethod getMethod() const {return modelReduction::PFA_IncrementalTangent;};
protected:
virtual void localConstitutive(const STensor3& F0, const STensor3& Fn, const ClusteringData *q0, ClusteringData *q1);
@@ -333,9 +349,9 @@ class HierarchicalReductionTFA : public ReductionTFA
virtual void getBiLinearTermHierarchical(int clusterId, const ClusteringData *q0, const ClusteringData *q1, std::vector<int>& othersId, fullMatrix<double>& dres) const;
virtual void getJac_rsHierarchical(int r, int s, const ClusteringData *q0, const ClusteringData *q1, STensor43& Jac_rs) const;
- virtual void getLinearTermPolarizationHierarchical(int clusterId, const ClusteringData *q0, const ClusteringData *q1, fullVector<double>& res) const;
- virtual void getBiLinearTermPolarizationHierarchical(int clusterId, const ClusteringData *q0, const ClusteringData *q1, std::vector<int>& othersId, fullMatrix<double>& dres) const;
- virtual void getJacPolarization_rsHierarchical(int r, int s, const ClusteringData *q0, const ClusteringData *q1, STensor43& Jac_rs) const;
+ virtual void getLinearTermHSHierarchical(int clusterId, const ClusteringData *q0, const ClusteringData *q1, fullVector<double>& res) const;
+ virtual void getBiLinearTermHSHierarchical(int clusterId, const ClusteringData *q0, const ClusteringData *q1, std::vector<int>& othersId, fullMatrix<double>& dres) const;
+ virtual void getJacHS_rsHierarchical(int r, int s, const ClusteringData *q0, const ClusteringData *q1, STensor43& Jac_rs) const;
void initSolveHierarchical(const ClusteringData *q0, const ClusteringData *q1, nonLinearSystemPETSc<double>* &NLsystem);
void preallocateSystemHierarchical(const ClusteringData *q0, const ClusteringData *q1, nonLinearSystemPETSc<double>* &NLsystem);