diff --git a/NonLinearSolver/materialLaw/mlawTFAMaterialLaws.cpp b/NonLinearSolver/materialLaw/mlawTFAMaterialLaws.cpp
index 24843d1a15a10537de969f3225dafb6c1c66a914..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
{
@@ -523,12 +529,6 @@ void mlawTFAHierarchicalMaterialLaws::setTFAMethod(int TFAMethodNum, int dim, in
Msg::Info("HierarchicalIncrementalTangentTFA is used !!!");
_reductionModel = new HierarchicalIncrementalTangentTFA(solverType, dim, 0, false, false, tag);
}
- /*else if(TFAMethodNum == 1)
- {
- polarization = true;
- Msg::Info("HierarchicalIncrementalTangentTFAMixed is used !!!");
- _reductionModel = new HierarchicalIncrementalTangentTFAMixed(solverType, dim, correction, false, true, tag);
- } */
}
mlawTFAHierarchicalMaterialLaws::~mlawTFAHierarchicalMaterialLaws()
diff --git a/NonLinearSolver/modelReduction/Clustering.cpp b/NonLinearSolver/modelReduction/Clustering.cpp
index 68f61e4be3653c3e82b51e09053ef4ae6d5c9780..c79eca2fe8d9156e5fafd29bca2664fab0e07d25 100644
--- a/NonLinearSolver/modelReduction/Clustering.cpp
+++ b/NonLinearSolver/modelReduction/Clustering.cpp
@@ -756,7 +756,8 @@ void Clustering::setNumberOfClustersForMaterialHierarchical(int matNum, int nbCl
}
void Clustering::solveStrainConcentration(nonLinearMechSolver& solver,
- int computeHomStiffness, int homStiffnessFromLocal,
+ int computeHomStiffness,
+ int homStiffnessFromLocal,
int getLocalOrientation,
bool saveToFile, const std::string strainConcentrationFileName,
bool saveLocalOrientationDataToFile, const std::string localOrientationFileName)
@@ -1280,7 +1281,7 @@ void Clustering::solveStrainConcentration(nonLinearMechSolver& solver,
fclose(f);
printf("done writing strain concentration tensor to file\n");
}
- if (saveLocalOrientationDataToFile and getLocalOrientation)
+ if (saveLocalOrientationDataToFile)
{
printf("start writing local orientation to file\n");
FILE* f = fopen(localOrientationFileName.c_str(),"w");
@@ -2298,8 +2299,6 @@ bool Clustering::computeClustersElasticity(int minNumberElementEachCluster, int
bool saveToFile, std::string clusterFileName,
bool saveSummaryFile, std::string summaryFileName,
bool saveMostStrainIPFile, const std::string saveMostStrainIPFileName,
- bool saveStrainConcStDFile, const std::string saveStrainConcStDFileName,
- bool saveClusterPositionFile, const std::string saveClusterPositionFileName,
bool saveClusterMeanFiberOrientationFile, const std::string saveClusterMeanFiberOrientationFileName
)
{
@@ -2368,17 +2367,7 @@ bool Clustering::computeClustersElasticity(int minNumberElementEachCluster, int
// most strain IP
mostStrainIP[clusterNum] = clutering_impls::maxConcentrationCluster(cluster);
- Msg::Info("most strain in cluster = %d is IP %d",clusterNum,mostStrainIP[clusterNum]);
-
- /*fullVector<double> cluster_position(3);
- double x_mean(0.), y_mean(0.), z_mean(0.);
- clutering_impls::meanClusterPosition(cluster,_mapGaussPoints,x_mean,y_mean,z_mean);
- cluster_position(0) = x_mean;
- cluster_position(1) = y_mean;
- cluster_position(2) = z_mean;
- _clusterPositionMap[clusterNum] = cluster_position;
- Msg::Info("position of cluster = %d is %.16g,%.16g,%.16g ",clusterNum,x_mean,y_mean,z_mean);*/
-
+ Msg::Info("most strain in cluster = %d is IP %d",clusterNum,mostStrainIP[clusterNum]);
/*fullVector<double> cluster_fiberOrientation(3);
double E1_mean(0.), E2_mean(0.), E3_mean(0.);
@@ -2417,28 +2406,7 @@ bool Clustering::computeClustersElasticity(int minNumberElementEachCluster, int
fclose(f);
printf("done writing data to file: %s \n",fname.c_str());
}
-
- if (saveClusterPositionFile)
- {
- std::vector<std::string> splitFileName = SplitFileName(saveClusterPositionFileName);
- std::string fname = splitFileName[1]+"_nbClusters_"+std::to_string(_totalNumberOfClusters)+splitFileName[2];
- FILE* f = fopen(fname.c_str(),"w");
- printf("start writing cluster position to file: %s\n",fname.c_str());
- std::string str ="CLUSTER_ID X Y Z";
- fprintf(f,"%s\n",str.c_str());
- for (std::map<int, fullVector<double> >::const_iterator it = _clusterPositionMap.begin(); it!= _clusterPositionMap.end(); it++)
- {
- int clusterId = it->first;
- fullVector<double> xyz = it->second;
- double x = xyz(0);
- double y = xyz(1);
- double z = xyz(2);
- fprintf(f, "%d %.16g %.16g %.16g\n",clusterId,x,y,z);
- }
- fclose(f);
- printf("done writing cluster position to file: %s \n",fname.c_str());
- }
-
+
if (saveClusterMeanFiberOrientationFile)
{
std::vector<std::string> splitFileName = SplitFileName(saveClusterMeanFiberOrientationFileName);
@@ -2523,48 +2491,6 @@ bool Clustering::computeClustersElasticity(int minNumberElementEachCluster, int
fclose(f);
printf("done writing cluster summary to file: %s \n",fname.c_str());
}
-
-
- if(saveStrainConcStDFile)
- {
- std::vector<std::string> splitFileName = SplitFileName(saveStrainConcStDFileName);
- std::string fname = splitFileName[1]+"_nbClusters_"+std::to_string(_totalNumberOfClusters)+splitFileName[2];
- FILE* f = fopen(fname.c_str(),"w");
- printf("start writing cluster std to file: %s\n",fname.c_str());
- std::string str ="CLUSTER_ID MAT_NUM VOLUME";
- // write data
- for (int i=0; i<6; i++)
- {
- for (int j=0; j<6; j++)
- {
- str += " STD_A_"+std::to_string(i)+std::to_string(j);
- }
- }
- fprintf(f,"%s\n",str.c_str());
- // write data
- for (std::map<int, STensor43>::const_iterator it = _clusterSTDTensorMap.begin(); it!= _clusterSTDTensorMap.end(); it++)
- {
- int cluster = it->first;
- int matNum = _clusterMaterialMap[cluster];
- double volume = _clusterVolumeMap[cluster];
- fprintf(f,"%d %d %.16g",cluster, matNum, volume);
- const STensor43& Astd = it->second;
- static fullMatrix<double> mat(6,6);
- STensorOperation::fromSTensor43ToFullMatrix(Astd,mat);
- printf("writing data for cluster %d\n",cluster);
- mat.print("std strain concentration tensor");
- for (int i=0; i<6; i++)
- {
- for (int j=0; j<6; j++)
- {
- fprintf(f, " %.16g",mat(i,j));
- }
- }
- fprintf(f, "\n");
- }
- fclose(f);
- printf("done writing cluster std to file: %s \n",fname.c_str());
- }
return true;
};
@@ -2846,9 +2772,6 @@ bool Clustering::computeClustersPlasticity(int minNumberElementEachCluster, int
bool saveSummaryFile, std::string summaryFileName,
bool savePlastStrainConcFile, std::string plastStrainconcFileName,
bool savePlastEqStrainConcFile, std::string plastEqstrainconcFileName,
- bool savePlastStrainConcSTDFile, std::string plastStrainconcstdevFileName,
- bool saveMostStrainIPFile, const std::string saveMostStrainIPFileName,
- bool saveClusterPositionFile, const std::string saveClusterPositionFileName,
bool saveClusterMeanFiberOrientationFile, const std::string saveClusterMeanFiberOrientationFileName
)
{
@@ -3040,45 +2963,6 @@ bool Clustering::computeClustersPlasticity(int minNumberElementEachCluster, int
}
}
//
- if (saveMostStrainIPFile)
- {
- std::vector<std::string> splitFileName = SplitFileName(saveMostStrainIPFileName);
- std::string fname = splitFileName[1]+"_nbClusters_"+std::to_string(_totalNumberOfClusters)+splitFileName[2];
- FILE* f = fopen(fname.c_str(),"w");
- std::string str ="CLUSTER_ID ELE_NUM GP_NUM";
- fprintf(f,"%s\n",str.c_str());
- for (std::map<int, int >::const_iterator it = mostStrainIP.begin(); it!= mostStrainIP.end(); it++)
- {
- int clusterId = it->first;
- int ipType = it->second;
- int ele, gp;
- numericalMaterialBase::getTwoIntsFromType(ipType,ele,gp);
- fprintf(f, "%d %d %d\n",clusterId, ele, gp);
- }
- fclose(f);
- printf("done writing data to file: %s \n",fname.c_str());
- }
-
- if (saveClusterPositionFile)
- {
- std::vector<std::string> splitFileName = SplitFileName(saveClusterPositionFileName);
- std::string fname = splitFileName[1]+"_nbClusters_"+std::to_string(_totalNumberOfClusters)+splitFileName[2];
- FILE* f = fopen(fname.c_str(),"w");
- printf("start writing cluster position to file: %s\n",fname.c_str());
- std::string str ="CLUSTER_ID X Y Z";
- fprintf(f,"%s\n",str.c_str());
- for (std::map<int, fullVector<double> >::const_iterator it = _clusterPositionMap.begin(); it!= _clusterPositionMap.end(); it++)
- {
- int clusterId = it->first;
- fullVector<double> xyz = it->second;
- double x = xyz(0);
- double y = xyz(1);
- double z = xyz(2);
- fprintf(f, "%d %.16g %.16g %.16g\n",clusterId,x,y,z);
- }
- fclose(f);
- printf("done writing cluster position to file: %s \n",fname.c_str());
- }
if (saveClusterMeanFiberOrientationFile)
{
@@ -3288,36 +3172,6 @@ bool Clustering::computeClustersPlasticity(int minNumberElementEachCluster, int
fclose(fplPowMean);
printf("done writing cluster plastic eq. power mean summary to file: %s \n",fnamePlPowMean.c_str());
}
-
- if (savePlastStrainConcSTDFile)
- {
- std::vector<std::string> splitFileNameSTD = SplitFileName(plastStrainconcstdevFileName);
- std::string fnamePlStd = splitFileNameSTD[1]+"_nbClusters_"+std::to_string(_totalNumberOfClusters)+splitFileNameSTD[2];
- FILE* fplstd = fopen(fnamePlStd.c_str(),"w");
- std::string strPlSTD ="CLUSTER_ID MAT_NUM VOLUME";
- // write data
- for (int i=0; i<36; i++)
- {
- strPlSTD += " STD_A_pl_"+std::to_string(i);
- }
- fprintf(fplstd,"%s\n",strPlSTD.c_str());
- // write data
- for (std::map<int, fullVector<double>>::const_iterator it = _clusterPlSTDVectorMap.begin(); it!= _clusterPlSTDVectorMap.end(); it++)
- {
- int cluster = it->first;
- int matNum = _clusterMaterialMap[cluster];
- double volume = _clusterVolumeMap[cluster];
- fprintf(fplstd,"%d %d %.16g",cluster, matNum, volume);
- const fullVector<double>& plasticSTD = it->second;
- for (int i=0; i<36; i++)
- {
- fprintf(fplstd, " %.16g",plasticSTD(i));
- }
- fprintf(fplstd, "\n");
- }
- fclose(fplstd);
- printf("done writing cluster plastic eq. std to file: %s \n",fnamePlStd.c_str());
- }
return true;
};
diff --git a/NonLinearSolver/modelReduction/Clustering.h b/NonLinearSolver/modelReduction/Clustering.h
index 0288a3466037e3153372782f29010be22bbeab4a..7fa9082f840066cc2f7972ef36a364fc9ca61dd4 100644
--- a/NonLinearSolver/modelReduction/Clustering.h
+++ b/NonLinearSolver/modelReduction/Clustering.h
@@ -217,8 +217,9 @@ class Clustering
Clustering();
// solve strain concentrartion tensor at gauss points and save to file
void solveStrainConcentration(nonLinearMechSolver& solver,
- int computeHomStiffness = 0, int homStiffnessFromLocal=0,
- int getLocalOrientation = 0,
+ int computeHomStiffness = 0, // compute homogenized elastic stiffness: 0=no; 1=yes
+ int homStiffnessFromLocal=0, // 0= compute hom el. stiffness from hom. stress, strain; 1=compute hom. el. stiffness as local average over C(x):A(x)
+ int getLocalOrientation = 0, // for materials with microstructure, e.g. woven composite
bool saveToFile=true, const std::string strainConcentrationFileName = "strainConcentrationData.csv",
bool saveLocalOrientationDataToFile=true, const std::string localOrientationFileName = "orientationData.csv");
void solveStrainConcentrationInelastic(nonLinearMechSolver& solverInelastic,bool saveToFile=true,
@@ -249,10 +250,8 @@ class Clustering
bool computeClustersElasticity(int minNumberElementEachCluster, int maxNbIterations=10000,
bool saveToFile=true, const std::string clusterFileName = "ClusterData.csv",
- bool saveSummaryFile=true, const std::string summaryFileName = "ClusterDataSummary.csv",
- bool saveMostStrainIPFile = true, const std::string saveMostStrainIPFileName = "MostStrainingIPSummary.csv",
- bool saveStrainConcStDFile = true, const std::string saveStrainConcStDFileName = "ClusterStrainConcStD.csv",
- bool saveClusterPositionFile = true, const std::string saveClusterPositionFileName = "ClusterPosition.csv",
+ bool saveSummaryFile=true, const std::string summaryFileName = "ClusterDataSummary.csv",
+ bool saveMostStrainIPFile=true, const std::string saveMostStrainIPFileName = "MostStrainingIPSummary.csv",
bool saveClusterMeanFiberOrientationFile=true, const std::string saveClusterMeanFiberOrientationFileName = "ClusterMeanFiberOrientation.csv"
);
@@ -266,9 +265,6 @@ class Clustering
bool saveSummaryFile=true, const std::string summaryFileName = "ClusterDataSummary.csv",
bool savePlastStrainConcFile=true, std::string plastStrainconcFileName = "ClusterPlasticSummary.csv",
bool savePlastEqStrainConcFile=true, std::string plastEqstrainconcFileName = "ClusterPlasticEqSummary.csv",
- bool savePlastStrainConcSTDFile=true, std::string plastStrainconcstdevFileName = "ClusterPlasticStrainTensorSTD.csv",
- bool saveMostStrainIPFile = true, const std::string saveMostStrainIPFileName = "MostStrainingIPSummary.csv",
- bool saveClusterPositionFile = true, const std::string saveClusterPositionFileName = "ClusterPosition.csv",
bool saveClusterMeanFiberOrientationFile=true, const std::string saveClusterMeanFiberOrientationFileName = "ClusterMeanFiberOrientation.csv"
);
diff --git a/NonLinearSolver/modelReduction/modelReduction.cpp b/NonLinearSolver/modelReduction/modelReduction.cpp
index 81a0d59da0e13b8101bb1bad2dc7d1cb2a5c57c5..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();
@@ -2160,21 +2277,6 @@ void IncrementalTangentTFA::localConstitutive(const STensor3& F0, const STensor3
{
law->constitutiveTFA(F0i,F1i,P1i,ipv0i,ipv1i,Tangenti,epsEig1i,depsEig1depsi);
}
-
- /*STensor3 Deps(Fn);
- Deps -= F0;
- const STensor43 &Aelistd = getClusterStrainConcentrationStDTensor(icl);
- const double &aelistd = getClusterStrainConcentrationStDScalar(icl);
- //Msg::Info("ASTD %.16g", aelistd);
- STensor3 epsi_std,F1i_max,F1i_min;
- //STensorOperation::multSTensor43STensor3(Aelistd,Deps,epsi_std);
- //epsi_std *= 2.;
- epsi_std = 2.*aelistd*Deps;
- F1i_max = F1i;
- F1i_max += epsi_std;
- F1i_min = F1i;
- F1i_min -= epsi_std;
- law->constitutiveTFA_corr(F0i,F1i,F1i_min,F1i_max,P1i,ipv0i,ipv1i,Tangenti,epsEig1i,depsEig1depsi);*/
}
};
@@ -2242,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++)
@@ -2275,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.);
@@ -2301,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);
@@ -2315,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;
@@ -2342,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;};
@@ -2362,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) {}
@@ -3448,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;
@@ -3497,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;
@@ -3514,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++)
@@ -3526,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);