diff --git a/dune/microstructure/CorrectorComputer.hh b/dune/microstructure/CorrectorComputer.hh
index 7aef0e24ad7ceb45a3f0ddd3508627afc1b2e944..f5a7e7efcaddb67fe2d5c3dd2cc031df43538dcc 100644
--- a/dune/microstructure/CorrectorComputer.hh
+++ b/dune/microstructure/CorrectorComputer.hh
@@ -321,7 +321,7 @@ public:
// defGradientV[l][2] = (1.0/gamma)*gradients[j][2]; //X3
defGradientV[l][2] = gradients[j][2]; //X3
- defGradientV = crossSectionDirectionScaling((1.0/gamma_),defGradientV);
+ defGradientV = sym(crossSectionDirectionScaling((1.0/gamma_),defGradientV));
// "phi*phi"-part
for (size_t k=0; k < dimworld; k++)
for (size_t i=0; i < nSf; i++)
@@ -333,7 +333,8 @@ public:
// defGradientU[k][2] = (1.0/gamma)*gradients[i][2]; //X3
defGradientU[k][2] = gradients[i][2]; //X3
// printmatrix(std::cout, defGradientU , "defGradientU", "--");
- defGradientU = crossSectionDirectionScaling((1.0/gamma_),defGradientU);
+ defGradientU = sym(crossSectionDirectionScaling((1.0/gamma_),defGradientU));
+
// auto etmp = material_.applyElasticityTensor(defGradientU,element.geometry().global(quadPos));
@@ -365,7 +366,19 @@ public:
// std::cout << "indicatorFunction_material_1: " << indicatorFunction_material_1(element.geometry().global(quadPos)) << std::endl;
// }
// double energyDensity= scalarProduct(material_.ElasticityTensorGlobal(defGradientU,element.geometry().global(quadPos)),sym(defGradientV));
- double energyDensity= scalarProduct(material_.ElasticityTensor(defGradientU,localIndicatorFunction(quadPos)),sym(defGradientV));
+
+
+ // std::cout << "scalarProduct(elasticityTensor(basisContainer[m],indicatorFunction(element.geometry().global(quadPos))),sym(defGradientV)) " << scalarProduct(elasticityTensor(basisContainer[m],indicatorFunction(element.geometry().global(quadPos))),sym(defGradientV)) << std::endl;
+ // std::cout << "linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), basisContainer[m], defGradientV)" << linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), basisContainer[m], defGradientV) << std::endl;
+ // }
+ std::cout << "--------------- \n";
+ printmatrix(std::cout, defGradientU, "defGradientU", "--");
+ printmatrix(std::cout, material_.applyL(defGradientU,localIndicatorFunction(quadPos)), "material_.applyL(defGradientU,localIndicatorFunction(quadPos))", "--");
+ printmatrix(std::cout, material_.ElasticityTensor(defGradientU,localIndicatorFunction(quadPos)), "material_.ElasticityTensor(defGradientU,localIndicatorFunction(quadPos)", "--");
+
+
+ double energyDensity= scalarProduct(material_.applyL(defGradientU,localIndicatorFunction(quadPos)),sym(defGradientV));
+ // double energyDensity= scalarProduct(material_.ElasticityTensor(defGradientU,localIndicatorFunction(quadPos)),sym(defGradientV));
// double energyDensity= scalarProduct(elasticityTensor(defGradientU,localIndicatorFunction(quadPos)),sym(defGradientV));
// double energyDensity= scalarProduct(elasticityTensor(defGradientU,indicatorFunction(element.geometry().global(quadPos))),sym(defGradientV));
// double energyDensity= scalarProduct(elasticityTensor(defGradientU,indicatorFunction(quadPos)),sym(defGradientV));
@@ -400,7 +413,13 @@ public:
// double energyDensityGphi = scalarProduct(material_.ElasticityTensorGlobal(basisContainer[m],element.geometry().global(quadPos)),sym(defGradientV));
- double energyDensityGphi = scalarProduct(material_.ElasticityTensor(basisContainer[m],localIndicatorFunction(quadPos)),sym(defGradientV));
+
+
+
+
+
+ double energyDensityGphi = scalarProduct(material_.applyL(basisContainer[m],localIndicatorFunction(quadPos)),sym(defGradientV));
+ // double energyDensityGphi = scalarProduct(material_.ElasticityTensor(basisContainer[m],localIndicatorFunction(quadPos)),sym(defGradientV));
// double energyDensityGphi= scalarProduct(elasticityTensor(basisContainer[m],localIndicatorFunction(quadPos)),sym(defGradientV));
// std::cout << "scalarProduct(elasticityTensor(basisContainer[m],indicatorFunction(element.geometry().global(quadPos))),sym(defGradientV))" << scalarProduct(elasticityTensor(basisContainer[m],indicatorFunction(element.geometry().global(quadPos))),sym(defGradientV)) <<std::endl;
// std::cout << "linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), basisContainer[m], defGradientV)" << linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), basisContainer[m], defGradientV) << std::endl;
@@ -444,7 +463,12 @@ public:
// }
// double energyDensityGG = scalarProduct(material_.ElasticityTensorGlobal(basisContainer[m],element.geometry().global(quadPos)),sym(basisContainer[n]));
- double energyDensityGG = scalarProduct(material_.ElasticityTensor(basisContainer[m],localIndicatorFunction(quadPos)),sym(basisContainer[n]));
+
+
+
+
+ double energyDensityGG = scalarProduct(material_.applyL(basisContainer[m],localIndicatorFunction(quadPos)),sym(basisContainer[n]));
+ // double energyDensityGG = scalarProduct(material_.ElasticityTensor(basisContainer[m],localIndicatorFunction(quadPos)),sym(basisContainer[n]));
// double energyDensityGG= scalarProduct(elasticityTensor(basisContainer[m],localIndicatorFunction(quadPos)),sym(basisContainer[n]));
// double energyDensityGG= scalarProduct(elasticityTensor(basisContainer[m],indicatorFunction(element.geometry().global(quadPos))),sym(basisContainer[n]));
// double energyDensityGG= scalarProduct(elasticityTensor(basisContainer[m],element.geometry().global(quadPos)),basisContainer[n]);
@@ -575,10 +599,14 @@ public:
// defGradientV[k][2] = (1.0/gamma_)*gradients[i][2]; //
defGradientV[k][2] = gradients[i][2]; // X3
- defGradientV = crossSectionDirectionScaling((1.0/gamma_),defGradientV);
+ defGradientV = sym(crossSectionDirectionScaling((1.0/gamma_),defGradientV));
// double energyDensity= scalarProduct(material_.ElasticityTensorGlobal((-1.0)*forceTerm(quadPos),element.geometry().global(quadPos)),sym(defGradientV));
- double energyDensity= scalarProduct(material_.ElasticityTensor((-1.0)*forceTerm(quadPos),localIndicatorFunction(quadPos)),sym(defGradientV));
+
+
+
+ double energyDensity= scalarProduct(material_.applyL((-1.0)*forceTerm(quadPos),localIndicatorFunction(quadPos)),sym(defGradientV));
+ // double energyDensity= scalarProduct(material_.ElasticityTensor((-1.0)*forceTerm(quadPos),localIndicatorFunction(quadPos)),sym(defGradientV));
// double energyDensity= scalarProduct(elasticityTensor((-1.0)*forceTerm(quadPos),localIndicatorFunction(quadPos)),sym(defGradientV));
// double energyDensity= scalarProduct(elasticityTensor((-1.0)*forceTerm(quadPos),indicatorFunction(element.geometry().global(quadPos))),sym(defGradientV));
// double energyDensity= scalarProduct(elasticityTensor((-1.0)*forceTerm(quadPos),indicatorFunction(quadPos)),sym(defGradientV));
@@ -601,7 +629,10 @@ public:
for (size_t m=0; m<3; m++)
{
// double energyDensityfG= scalarProduct(material_.ElasticityTensorGlobal((-1.0)*forceTerm(quadPos),element.geometry().global(quadPos)),sym(basisContainer[m]));
- double energyDensityfG= scalarProduct(material_.ElasticityTensor((-1.0)*forceTerm(quadPos),localIndicatorFunction(quadPos)),sym(basisContainer[m]));
+
+
+ double energyDensityfG= scalarProduct(material_.applyL((-1.0)*forceTerm(quadPos),localIndicatorFunction(quadPos)),sym(basisContainer[m]));
+ // double energyDensityfG= scalarProduct(material_.ElasticityTensor((-1.0)*forceTerm(quadPos),localIndicatorFunction(quadPos)),sym(basisContainer[m]));
// double energyDensityfG= scalarProduct(elasticityTensor((-1.0)*forceTerm(quadPos),localIndicatorFunction(quadPos)),sym(basisContainer[m]));
// double energyDensityfG= scalarProduct(elasticityTensor((-1.0)*forceTerm(quadPos),indicatorFunction(element.geometry().global(quadPos))),sym(basisContainer[m]));
// double energyDensityfG= scalarProduct(elasticityTensor((-1.0)*forceTerm(quadPos),indicatorFunction(quadPos)),sym(basisContainer[m]));
@@ -659,11 +690,14 @@ public:
//TEST
//Check Element-Stiffness-Symmetry:
- for (size_t i=0; i<localPhiOffset; i++)
- for (size_t j=0; j<localPhiOffset; j++ )
+ if(parameterSet_.get<bool>("print_debug", false))
{
- if(abs(elementMatrix[i][j] - elementMatrix[j][i]) > 1e-12 )
- std::cout << "ELEMENT-STIFFNESS MATRIX NOT SYMMETRIC!!!" << std::endl;
+ for (size_t i=0; i<localPhiOffset; i++)
+ for (size_t j=0; j<localPhiOffset; j++ )
+ {
+ if(abs(elementMatrix[i][j] - elementMatrix[j][i]) > 1e-12 )
+ std::cout << "ELEMENT-STIFFNESS MATRIX NOT SYMMETRIC!!!" << std::endl;
+ }
}
//////////////////////////////////////////////////////////////////////////////
// GLOBAL STIFFNES ASSEMBLY
diff --git a/dune/microstructure/materialDefinitions.hh b/dune/microstructure/materialDefinitions.hh
index f4e6ee147907235ef1ce01212139d2e0a668a59f..3f13916da25492e484fcee17bea00b28d4256f5e 100644
--- a/dune/microstructure/materialDefinitions.hh
+++ b/dune/microstructure/materialDefinitions.hh
@@ -18,6 +18,11 @@ using MatrixRT = FieldMatrix< double, 3, 3>;
using VectorRT = FieldVector< double, 3>;
+
+
+
+
+
// ----------------------------------------------------------------------------------
template<class Domain>
int indicatorFunction_material_neukamm(const Domain& x)
diff --git a/dune/microstructure/prestrainedMaterial.hh b/dune/microstructure/prestrainedMaterial.hh
index 9ce64e42022b1f616825ee17695e6a3c70a6fd8b..5bfb83eddc4ffac00957366fdcb9743628949167 100644
--- a/dune/microstructure/prestrainedMaterial.hh
+++ b/dune/microstructure/prestrainedMaterial.hh
@@ -23,6 +23,59 @@ using std::shared_ptr;
using std::make_shared;
+
+
+
+
+
+
+
+// // ----------------------------------------------------------------------------------
+// template<class Domain>
+// int indicatorFunction_material(const Domain& x)
+// {
+// double theta=0.25;
+// if (x[0] <(-0.5+theta) and x[2]<(-0.5+theta))
+// return 1; //#Phase1
+// else if (x[1]<(-0.5+theta) and x[2]>(0.5-theta))
+// return 2; //#Phase2
+// else
+// return 3; //#Phase3
+// }
+// MatrixRT getL(const int& phase)
+// {
+// const FieldVector<double,3> mu = {80.0, 80.0, 60.0};
+// const FieldVector<double,3> lambda = {80.0, 80.0, 25.0};
+
+
+
+// if (phase == 1)
+// return 2.0 * mu[0] * sym(G) + lambda[0] * trace(sym(G)) * Id(); //#Phase1 //Isotrop(mu,lambda)
+// else if (phase == 2)
+// return 2.0 * mu[1] * sym(G) + lambda[1] * trace(sym(G)) * Id(); //#Phase2
+// else
+// return 2.0 * mu[2] * sym(G) + lambda[2] * trace(sym(G)) * Id(); //#Phase3
+// }
+// MatrixRT prestrain_material(const int& phase)
+// {
+// if (phase == 1)
+// return {{1.0, 0.0, 0.0}, //#Phase1
+// {0.0, 1.0, 0.0},
+// {0.0, 0.0, 1.0}
+// };
+// else if (phase == 2)
+// return {{1.0, 0.0, 0.0}, //#Phase2
+// {0.0, 1.0, 0.0},
+// {0.0, 0.0, 1.0}
+// };
+// else
+// return {{0.0, 0.0, 0.0}, //#Phase3
+// {0.0, 0.0, 0.0},
+// {0.0, 0.0, 0.0}
+// };
+// }
+
+
//--------------------------------------------------------
template <class GridView>
class prestrainedMaterial
@@ -45,6 +98,8 @@ public:
using MatrixFunc = std::function< MatrixRT(const MatrixRT&) >;
using MatrixPhase = std::function< MatrixRT(const int&) >;
+ // using VoigtMatrix = FieldMatrix< double, 6, 6>;
+
protected:
const GridView& gridView_;
@@ -54,6 +109,11 @@ protected:
// MatrixFunc L1_;
// MatrixFunc L2_;
// MatrixFunc L3_;
+
+ FieldMatrix<double,6,6> L1_;
+ FieldMatrix<double,6,6> L2_;
+ FieldMatrix<double,6,6> L3_;
+
// const FuncScalar mu_;
// const FuncScalar lambda_;
// const int phases_;
@@ -67,6 +127,25 @@ protected:
GridViewFunction<int(const Domain&), GridView> indicatorFunctionGVF_;
LocalFunction<int(const Domain&), typename GridViewEntitySet<GridView, 0>::Element > localIndicatorFunction_;
+
+ //Transformation matrix for Voigt notation
+ FieldMatrix<double,6,6> D = {{1.0, 0.0, 0.0, 0.0 , 0.0, 0.0},
+ {0.0, 1.0, 0.0, 0.0 , 0.0, 0.0},
+ {0.0, 0.0, 1.0, 0.0 , 0.0, 0.0},
+ {0.0, 0.0, 0.0, sqrt(2.0) , 0.0, 0.0},
+ {0.0, 0.0, 0.0, 0.0 , sqrt(2.0), 0.0},
+ {0.0, 0.0, 0.0, 0.0 , 0.0, sqrt(2.0)}
+ };
+ FieldMatrix<double,6,6> D_inv = {{1.0, 0.0, 0.0, 0.0 , 0.0, 0.0},
+ {0.0, 1.0, 0.0, 0.0 , 0.0, 0.0},
+ {0.0, 0.0, 1.0, 0.0 , 0.0, 0.0},
+ {0.0, 0.0, 0.0, 1.0/sqrt(2.0) , 0.0, 0.0},
+ {0.0, 0.0, 0.0, 0.0 , 1.0/sqrt(2.0), 0.0},
+ {0.0, 0.0, 0.0, 0.0 , 0.0, 1.0/sqrt(2.0)}
+ };
+
+
+
public:
///////////////////////////////
// constructor
@@ -84,6 +163,41 @@ public:
indicatorFunctionGVF_ = Dune::Functions::makeGridViewFunction(indicatorFunction_, gridView_);
localIndicatorFunction_ = localFunction(indicatorFunctionGVF_);
+ const FieldVector<double,3> mu = {80.0, 80.0, 60.0};
+ const FieldVector<double,3> lambda = {80.0, 80.0, 25.0};
+
+ L1_ = {{lambda[0]+2.0*mu[0], lambda[0], lambda[0], 0.0 , 0.0, 0.0},
+ {lambda[0], lambda[0]+2*mu[0], lambda[0], 0.0 ,0.0 ,0.0},
+ {lambda[0], lambda[0], lambda[0]+2.0*mu[0], 0.0, 0.0, 0.0},
+ { 0.0 ,0.0, 0.0, 2.0*mu[0], 0.0, 0.0},
+ { 0.0 ,0.0, 0.0, 0.0, 2.0*mu[0], 0.0},
+ { 0.0 ,0.0, 0.0, 0.0, 0.0, 2.0*mu[0]}
+ };
+
+ L2_ = {{lambda[1]+2.0*mu[1], lambda[1], lambda[1], 0.0 , 0.0, 0.0},
+ {lambda[1], lambda[1]+2*mu[1], lambda[1], 0.0 ,0.0 ,0.0},
+ {lambda[1], lambda[1], lambda[1]+2.0*mu[1], 0.0, 0.0, 0.0},
+ { 0.0 ,0.0, 0.0, 2.0*mu[1], 0.0, 0.0},
+ { 0.0 ,0.0, 0.0, 0.0, 2.0*mu[1], 0.0},
+ { 0.0 ,0.0, 0.0, 0.0, 0.0, 2.0*mu[1]}
+ };
+
+ L3_ = {{lambda[2]+2.0*mu[2], lambda[2], lambda[2], 0.0 , 0.0, 0.0},
+ {lambda[2], lambda[2]+2.0*mu[2], lambda[2], 0.0 ,0.0 ,0.0},
+ {lambda[2], lambda[2], lambda[2]+2.0*mu[2], 0.0, 0.0, 0.0},
+ { 0.0 ,0.0, 0.0, 2.0*mu[2], 0.0, 0.0},
+ { 0.0 ,0.0, 0.0, 0.0, 2.0*mu[2], 0.0},
+ { 0.0 ,0.0, 0.0, 0.0, 0.0, 2.0*mu[2]}
+ };
+
+
+ // printmatrix(std::cout, L1_, "L1_", "--");
+ // L1_.leftmultiply(D_inv);
+ // printmatrix(std::cout, L1_, "L1_", "--");
+ // L1_.rightmultiply(D);
+ // printmatrix(std::cout, L1_, "L1_", "--");
+
+
// Python::Module module = Python::import(materialFunctionName_);
// elasticityTensor_ = Python::make_function<MatrixRT>(module.get("L"));
// elasticityTensor_ = setupElasticityTensor(materialFunctionName_);
@@ -130,6 +244,51 @@ public:
DUNE_THROW(Exception, "There exists no material with this name ");
}
+//////////// TESTING //////////////////////////////////////
+
+
+ MatrixRT applyL(const MatrixRT& G, const int& phase) const
+ {
+ // FieldVector<double,6> G_tmp = {G[0][0], G[1][1], G[2][2], sqrt(2.0)*G[1][2], sqrt(2.0)*G[0][2], sqrt(2.0)*G[0][1]};
+ FieldVector<double,6> G_tmp = {G[0][0], G[1][1], G[2][2], G[1][2], G[0][2], G[0][1]};
+ FieldVector<double,6> out(0);
+
+ printvector(std::cout, G_tmp, "G_tmp", "--");
+
+ //
+ // printmatrix(std::cout, L2_, "L2_", "--");
+ // printmatrix(std::cout, L3_, "L3_", "--");
+
+ if (phase == 1)
+ {
+ printmatrix(std::cout, L1_, "L1_", "--");
+ L1_.mv(G_tmp,out);
+ }
+ else if (phase == 2)
+ {
+ printmatrix(std::cout, L2_, "L2_", "--");
+ L2_.mv(G_tmp,out);
+ }
+ else
+ {
+ printmatrix(std::cout, L3_, "L3_", "--");
+ L3_.mv(G_tmp,out);
+ }
+
+
+ printvector(std::cout, out, "out", "--");
+
+ // return image as Matrix again
+ // return {{out[0], (1.0/sqrt(2.0))*out[5], (1.0/sqrt(2.0))*out[4]}, {(1.0/sqrt(2.0))*out[5], out[1], (1.0/sqrt(2.0))*out[3]}, {(1.0/sqrt(2.0))*out[4], (1.0/sqrt(2.0))*out[3], out[2]}};
+ // return {{out[0], (1.0/2.0)*out[5], (1.0/2.0)*out[4]}, {(1.0/2.0)*out[5], out[1], (1.0/2.0)*out[3]}, {(1.0/2.0)*out[4], (1.0/2.0)*out[3], out[2]}};
+ return {{out[0], out[5], out[4]}, {out[5], out[1], out[3]}, {out[4], out[3], out[2]}};
+ }
+
+
+
+
+////////////////////////////////////////////////////////
+
//--- Apply elasticityTensor_ to input Matrix G at material phase
// input: Matrix G, material-phase