diff --git a/dune/microstructure/CorrectorComputer.hh b/dune/microstructure/CorrectorComputer.hh
index 180552f301d6ec515bdd517ee95a24f3baf0c589..e2864169e49bb6362de64613e8fc49a71a572d15 100644
--- a/dune/microstructure/CorrectorComputer.hh
+++ b/dune/microstructure/CorrectorComputer.hh
@@ -270,8 +270,8 @@ public:
auto elasticityTensor = material_.getElasticityTensor();
// auto indicatorFunction = material_.getIndicatorFunction();
- // auto localIndicatorFunction = material_.getLocalIndicatorFunction();
- // localIndicatorFunction.bind(element);
+ auto localIndicatorFunction = material_.getLocalIndicatorFunction();
+ localIndicatorFunction.bind(element);
// auto indicatorFunctionGVF = material_.getIndicatorFunctionGVF();
// auto indicatorFunction = localFunction(indicatorFunctionGVF);
// indicatorFunction.bind(element);
@@ -377,8 +377,8 @@ public:
// std::cout << "material_.applyIndicatorFunction(element.geometry().global(quadPos)): " << material_.applyIndicatorFunction(element.geometry().global(quadPos)) << std::endl;
// 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_.ElasticityTensorPhase(defGradientU,localIndicatorFunction(quadPos)),sym(defGradientV));
+ // double energyDensity= scalarProduct(material_.ElasticityTensorGlobal(defGradientU,element.geometry().global(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));
@@ -412,8 +412,8 @@ public:
// }
- double energyDensityGphi = scalarProduct(material_.ElasticityTensorGlobal(basisContainer[m],element.geometry().global(quadPos)),sym(defGradientV));
- // double energyDensityGphi = scalarProduct(material_.ElasticityTensorPhase(basisContainer[m],localIndicatorFunction(quadPos)),sym(defGradientV));
+ // 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(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;
@@ -456,8 +456,8 @@ public:
// std::cout << "linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), basisContainer[m], basisContainer[n])" << linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), basisContainer[m], basisContainer[n])<< std::endl;
// }
- double energyDensityGG = scalarProduct(material_.ElasticityTensorGlobal(basisContainer[m],element.geometry().global(quadPos)),sym(basisContainer[n]));
- // double energyDensityGG = scalarProduct(material_.ElasticityTensorPhase(basisContainer[m],localIndicatorFunction(quadPos)),sym(basisContainer[n]));
+ // 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(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]);
@@ -506,8 +506,8 @@ public:
auto elasticityTensor = material_.getElasticityTensor();
// auto indicatorFunction = material_.getIndicatorFunction();
- // auto localIndicatorFunction = material_.getLocalIndicatorFunction();
- // localIndicatorFunction.bind(element);
+ auto localIndicatorFunction = material_.getLocalIndicatorFunction();
+ localIndicatorFunction.bind(element);
// // auto indicatorFunction = material_.getLocalIndicatorFunction();
// auto indicatorFunctionGVF = material_.getIndicatorFunctionGVF();
// auto indicatorFunction = localFunction(indicatorFunctionGVF);
@@ -585,8 +585,8 @@ public:
defGradientV = crossSectionDirectionScaling((1.0/gamma_),defGradientV);
- double energyDensity= scalarProduct(material_.ElasticityTensorGlobal((-1.0)*forceTerm(quadPos),element.geometry().global(quadPos)),sym(defGradientV));
- // double energyDensity= scalarProduct(material_.ElasticityTensorPhase((-1.0)*forceTerm(quadPos),localIndicatorFunction(quadPos)),sym(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(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));
@@ -608,8 +608,8 @@ public:
// "f*m"-part
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_.ElasticityTensorPhase((-1.0)*forceTerm(quadPos),localIndicatorFunction(quadPos)),sym(basisContainer[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(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]));
diff --git a/dune/microstructure/prestrain_material_geometry.hh b/dune/microstructure/prestrain_material_geometry.hh
index 792215f38214dd9c63f371f4f48bbae7f31d4c02..55def935a783ea2ae1fb3db94317ac503ae686b8 100644
--- a/dune/microstructure/prestrain_material_geometry.hh
+++ b/dune/microstructure/prestrain_material_geometry.hh
@@ -54,7 +54,8 @@ public:
}
else if (imp == "material_neukamm"){
// std::array<double,2> mu = parameters.get<std::array<double,2>>("mu", {1.0,3.0});
- std::array<double,3> mu = parameters.get<std::array<double,3>>("mu", {1.0,3.0,2.0});
+ // std::array<double,3> mu = parameters.get<std::array<double,3>>("mu", {1.0,3.0,2.0});
+ std::array<double,3> mu = parameters.get<std::array<double,3>>("mu", {80.0, 80.0, 60.0});
// Python::Module module = Python::import(parameters.get<std::string>("material_neukamm"));
Python::Module module = Python::import("material_neukamm");
@@ -495,7 +496,9 @@ public:
}
if (imp == "material_neukamm"){
// std::array<double,2> lambda = parameters.get<std::array<double,2>>("lambda", {1.0,3.0});
- std::array<double,3> lambda = parameters.get<std::array<double,3>>("lambda", {1.0,3.0,2.0});
+ // std::array<double,3> lambda = parameters.get<std::array<double,3>>("lambda", {1.0,3.0,2.0});
+ std::array<double,3> lambda = parameters.get<std::array<double,3>>("lambda", {80.0, 80.0, 25.0});
+
Python::Module module = Python::import("material_neukamm");
auto indicatorFunction = Python::make_function<double>(module.get("f"));
diff --git a/dune/microstructure/prestrainedMaterial.hh b/dune/microstructure/prestrainedMaterial.hh
index 90d09e54ca8a22915b8b012d644b0c408baf6a42..2bfc9fa37174096b8c546977b0b0f0c7f8cccaab 100644
--- a/dune/microstructure/prestrainedMaterial.hh
+++ b/dune/microstructure/prestrainedMaterial.hh
@@ -27,50 +27,124 @@ using std::make_shared;
-
-
-
-
+ // ----------------------------------------------------------------------------------
template<class Domain>
- int indicatorFunction_material_1(const Domain& x)
+ int indicatorFunction_material_neukamm(const Domain& x)
{
- // std::cout << "x[0] , x[1] , x[2]" << x[0] << " , " << x[1] << ", "<< x[2] << std::endl;
- // std::cout << "-1/2: " << -1/2 << std::endl;
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 0; //#Phase3
+ return 3; //#Phase3
}
-
- MatrixRT material_1(const MatrixRT& G, const int& phase)
+ MatrixRT material_neukamm(const MatrixRT& G, const int& phase)
{
- const FieldVector<double,3> mu = {80.0, 80.0, 60.0};
+ 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
+ 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_neukamm(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 Domain>
+ int indicatorFunction_two_phase_material_1(const Domain& x)
+ {
+ double theta=0.25;
+ if(abs(x[0]) > (theta/2.0))
+ return 1; //#Phase1
+ else
+ return 0; //#Phase2
+ }
+ MatrixRT two_phase_material_1(const MatrixRT& G, const int& phase)
+ {
+ const FieldVector<double,2> mu = {80.0, 160.0};
+ const FieldVector<double,2> lambda = {80.0, 160.0};
if (phase == 1)
return 2.0 * mu[0] * sym(G) + lambda[0] * trace(sym(G)) * Id();
- if (phase == 2)
+ else
return 2.0 * mu[1] * sym(G) + lambda[1] * trace(sym(G)) * Id();
+ }
+ MatrixRT prestrain_two_phase_material_1(const int& phase)
+ {
+ const FieldVector<double,2> rho = {3.0, 5.0};
+ if (phase == 1)
+ return {{rho[0], 0.0, 0.0}, //#Phase1
+ {0.0, rho[0], 0.0},
+ {0.0, 0.0, rho[0]}
+ };
else
- return 2.0 * mu[2] * sym(G) + lambda[2] * trace(sym(G)) * Id();
+ return {{rho[1], 0.0, 0.0}, //#Phase2
+ {0.0, rho[1], 0.0},
+ {0.0, 0.0, rho[1]}
+ };
}
+ // ----------------------------------------------------------------------------------
+ // ----------------------------------------------------------------------------------
template<class Domain>
- int indicatorFunction_material_homogeneous(const Domain& x)
+ int indicatorFunction_two_phase_material_2(const Domain& x)
{
- return 0;
+ double theta=0.25;
+ if(abs(x[0]) > (theta/2.0))
+ return 1; //#Phase1
+ else
+ return 0; //#Phase2
}
- MatrixRT material_homogeneous(const MatrixRT& G, const int& phase)
- {
- const FieldVector<double,1> mu = {80.0};
- const FieldVector<double,1> lambda = {80.0};
- return 2.0 * mu[0] * sym(G) + lambda[0] * trace(sym(G)) * Id();
- }
+ MatrixRT two_phase_material_2(const MatrixRT& G, const int& phase)
+ {
+ const FieldVector<double,2> mu = {5.0, 15.0};
+ const FieldVector<double,2> lambda = {6.0, 8.0};
+ if (phase == 1)
+ return 2.0 * mu[0] * sym(G) + lambda[0] * trace(sym(G)) * Id();
+ else
+ return 2.0 * mu[1] * sym(G) + lambda[1] * trace(sym(G)) * Id();
+ }
+ MatrixRT prestrain_two_phase_material_2(const int& phase)
+ {
+ const FieldVector<double,2> rho = {3.0, 5.0};
+ if (phase == 1)
+ return {{rho[0], 0.0, 0.0}, //#Phase1
+ {0.0, rho[0], 0.0},
+ {0.0, 0.0, rho[0]}
+ };
+ else
+ return {{rho[1], 0.0, 0.0}, //#Phase2
+ {0.0, rho[1], 0.0},
+ {0.0, 0.0, rho[1]}
+ };
+ }
+ // ----------------------------------------------------------------------------------
+
+ // ----------------------------------------------------------------------------------
template<class Domain>
int indicatorFunction_material_2(const Domain& x)
{
@@ -80,16 +154,64 @@ using std::make_shared;
else
return 0; //#Phase2
}
-
MatrixRT material_2(const MatrixRT& G, const int& phase)
{
- const FieldVector<double,2> mu = {80.0, 160.0};
+ const FieldVector<double,2> mu = {80.0, 160.0};
const FieldVector<double,2> lambda = {80.0, 160.0};
if (phase == 1)
- return 2.0 * mu[0] * sym(G) + lambda[0] * trace(sym(G)) * Id();
+ return 2.0 * mu[0] * sym(G) + lambda[0] * trace(sym(G)) * Id(); //#Phase1
else
- return 2.0 * mu[1] * sym(G) + lambda[1] * trace(sym(G)) * Id();
+ return 2.0 * mu[1] * sym(G) + lambda[1] * trace(sym(G)) * Id(); //#Phase2
+ }
+ MatrixRT prestrain_material_2(const int& phase)
+ {
+ const FieldVector<double,2> rho = {3.0, 5.0};
+ if (phase == 1)
+ return {{rho[0], 0.0, 0.0}, //#Phase1
+ {0.0, rho[0], 0.0},
+ {0.0, 0.0, rho[0]}
+ };
+ else
+ return {{rho[1], 0.0, 0.0}, //#Phase2
+ {0.0, rho[1], 0.0},
+ {0.0, 0.0, rho[1]}
+ };
+ }
+ // ----------------------------------------------------------------------------------
+
+
+ // ----------------------------------------------------------------------------------
+ template<class Domain>
+ int indicatorFunction_material_homogeneous(const Domain& x)
+ {
+ return 0;
+ }
+ MatrixRT material_homogeneous(const MatrixRT& G, const int& phase)
+ {
+ const FieldVector<double,1> mu = {80.0};
+ const FieldVector<double,1> lambda = {80.0};
+ return 2.0 * mu[0] * sym(G) + lambda[0] * trace(sym(G)) * Id();
+ }
+ MatrixRT prestrain_homogeneous(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 {{1.0, 0.0, 0.0}, //#Phase3
+ {0.0, 1.0, 0.0},
+ {0.0, 0.0, 1.0}
+ };
}
+ // ----------------------------------------------------------------------------------
+
@@ -107,8 +229,6 @@ using std::make_shared;
// return 2.0 * mu[2] * sym(G) + lambda[2] * trace(sym(G)) * Id();
// }
-
-
template <class GridView> // needed for GridViewFunctions
class prestrainedMaterial
{
@@ -129,38 +249,25 @@ public:
using Func2TensorParam = std::function< MatrixRT(const MatrixRT& ,const Domain&) >;
using Func2TensorPhase = std::function< MatrixRT(const MatrixRT& ,const int&) >;
using MatrixFunc = std::function< MatrixRT(const MatrixRT&) >;
+ using MatrixPhase = std::function< MatrixRT(const int&) >;
protected:
const GridView& gridView_;
const ParameterTree& parameterSet_;
+ std::string materialFunctionName_;
- int Phases_;
-
- MatrixFunc L1_;
- MatrixFunc L2_;
- MatrixFunc L3_;
-
-
- // const FieldVector<double , ...number of mu-Values/Phases> .. schwierig zur compile-time
-
+ // MatrixFunc L1_;
+ // MatrixFunc L2_;
+ // MatrixFunc L3_;
// const FuncScalar mu_;
// const FuncScalar lambda_;
- // double gamma_;
-
- std::string materialFunctionName_;
- // std::string materialFunctionName_ = parameterSet_.get<std::string>("materialFunction", "material");
-
- // --- Number of material phases?
// const int phases_;
-
- // Func2Tensor materialFunction_; //actually not needed??
-
// Func2Tensor elasticityTensor_;
// Func2TensorParam elasticityTensor_;
Func2TensorPhase elasticityTensor_;
+ MatrixPhase prestrain_;
- // FuncScalar indicatorFunction_;
Func2int indicatorFunction_;
GridViewFunction<int(const Domain&), GridView> indicatorFunctionGVF_;
@@ -178,21 +285,17 @@ public:
std::string materialFunctionName_ = parameterSet.get<std::string>("materialFunction", "material_1");
std::cout << "materialFunctionName_ : " << materialFunctionName_ << std::endl;
- // Python::Module module = Python::import(materialFunctionName_);
- // elasticityTensor_ = Python::make_function<MatrixRT>(module.get("L"));
-
- // elasticityTensor_ = setupElasticityTensor(materialFunctionName_);
- setupElasticityTensor(materialFunctionName_);
-
+ setupMaterial(materialFunctionName_);
indicatorFunctionGVF_ = Dune::Functions::makeGridViewFunction(indicatorFunction_, gridView_);
localIndicatorFunction_ = localFunction(indicatorFunctionGVF_);
+ // Python::Module module = Python::import(materialFunctionName_);
+ // elasticityTensor_ = Python::make_function<MatrixRT>(module.get("L"));
+ // elasticityTensor_ = setupElasticityTensor(materialFunctionName_);
// auto indicatorFunction = Python::make_function<int>(module.get("indicatorFunction"));
// indicatorFunction_ = Python::make_function<int>(module.get("indicatorFunction"));
// indicatorFunction_ = Dune::Functions::makeGridViewFunction(indicatorFunction , gridView_);
-
-
// module.get("Phases").toC<int>(Phases_);
// L1_ = Python::make_function<MatrixRT>(module.get("L1"));
// L2_ = Python::make_function<MatrixRT>(module.get("L2"));
@@ -200,8 +303,52 @@ public:
// bool isotropic_ = true; // read from module File
// auto Test = Dune::Functions::makeGridViewFunction(MAT<Domain> , gridView_); //not working
}
-/////////////////////////////////////////////////////////////////////////////////////////////////////////
+ //---function that determines elasticity Tensor
+ void setupMaterial(const std::string name) // cant use materialFunctionName_ here!?
+ {
+ // std::cout << "Using material definition:" << name << std::endl;
+ if(name == "material_neukamm")
+ {
+ indicatorFunction_ = indicatorFunction_material_neukamm<Domain>;
+ elasticityTensor_ = material_neukamm;
+ prestrain_ = prestrain_material_neukamm;
+ }
+ else if(name == "two_phase_material_1")
+ {
+ indicatorFunction_ = indicatorFunction_two_phase_material_1<Domain>;
+ elasticityTensor_ = two_phase_material_1;
+ prestrain_ = prestrain_two_phase_material_1;
+ }
+ else if(name == "two_phase_material_2")
+ {
+ indicatorFunction_ = indicatorFunction_two_phase_material_2<Domain>;
+ elasticityTensor_ = two_phase_material_2;
+ prestrain_ = prestrain_two_phase_material_2;
+ }
+ else if(name == "homogeneous")
+ {
+ indicatorFunction_ = indicatorFunction_material_homogeneous<Domain>;
+ elasticityTensor_ = material_homogeneous;
+ prestrain_ = prestrain_homogeneous;
+ }
+ else
+ DUNE_THROW(Exception, "There exists no material with this name ");
+ }
+
+
+ //--- apply elasticityTensor_ to input Matrix G at position x
+ // input: Matrix G, material-phase
+ MatrixRT ElasticityTensor(const MatrixRT& G, const int& phase) const
+ {
+ return elasticityTensor_(G,phase);
+
+ }
+
+ int IndicatorFunction(const Domain& x) const
+ {
+ return indicatorFunction_(x);
+ }
// static int indicatorFunction_material_1(const Domain& x)
// {
@@ -224,39 +371,6 @@ public:
- //---function that determines elasticity Tensor
- void setupElasticityTensor(const std::string name) // cant use materialFunctionName_ here!?
- {
- std::cout << "Using material definition:" << name << std::endl;
- if(name == "material_1")
- {
- // indicatorFunctionGVF_ = Dune::Functions::makeGridViewFunction(indicatorFunction_material_1, gridView_);
- // indicatorFunction_ = indicatorFunction_material_1;
-
- indicatorFunction_ = indicatorFunction_material_1<Domain>;
- elasticityTensor_ = material_1;
- }
- else if(name == "homogeneous")
- {
- // indicatorFunctionGVF_ = Dune::Functions::makeGridViewFunction(indicatorFunction_material_1, gridView_);
- // indicatorFunction_ = indicatorFunction_material_1;
- // indicatorFunctionGVF_ = Dune::Functions::makeGridViewFunction(indicatorFunction_material_homogeneous<Domain>, gridView_);
- indicatorFunction_ = indicatorFunction_material_homogeneous<Domain>;
- elasticityTensor_ = material_homogeneous;
- }
- else if(name == "material_2")
- {
- // indicatorFunctionGVF_ = Dune::Functions::makeGridViewFunction(indicatorFunction_material_1, gridView_);
- // indicatorFunction_ = indicatorFunction_material_1;
- // indicatorFunctionGVF_ = Dune::Functions::makeGridViewFunction(indicatorFunction_material_2<Domain>, gridView_);
- indicatorFunction_ = indicatorFunction_material_2<Domain>;
- elasticityTensor_ = material_2;
- }
- else
- DUNE_THROW(Exception, "There exists no material in materialDefinitions.hh with this name ");
- }
-
-
// MatrixRT localElasticityTensor(const MatrixRT& G, const Domain& x, Element& element) const //local Version ... muss binding öfters durchführen
// {
// localIndicatorFunction_.bind(*element);
@@ -265,46 +379,36 @@ public:
- //--- apply elasticityTensor_ to input Matrix G at position x
- MatrixRT ElasticityTensorPhase(const MatrixRT& G, const int& phase) const //global Version (takes global coordinates)
- {
- return elasticityTensor_(G,phase);
-
- }
- MatrixRT ElasticityTensorGlobal(const MatrixRT& G, const Domain& x) const //global Version (takes global coordinates)
- {
- return MAT(G,x);
-
- }
- //--- apply elasticityTensor_ to input Matrix G at position x
- MatrixRT applyElasticityTensor(const MatrixRT& G, const Domain& x) const //global Version (takes global coordinates)
- {
- // Python::Module module = Python::import("material");
-
- // auto indicatorFunctionTest = Python::make_function<double>(module.get("indicatorFunction"));
- // return elasticityTensor_(G,indicatorFunctionTest(x));
- // auto IFunction = localFunction(indicatorFunction_);
- // auto IFunction = this->getIndicatorFunction();
- // auto tmp = Dune::Functions::makeGridViewFunction(indicatorFunction_material_1, gridView_);
- // auto tmpLocal = LocalF
- // return elasticityTensor_(G,IFunction(x));
- // return elasticityTensor_(G,localIndicatorFunction_(x));
- return elasticityTensor_(G,indicatorFunction_(x));
- // return elasticityTensor_(G,indicatorFunctionGVF_(x));
- // int phase = indicatorFunction_(x);
- // auto tmp = elasticityTensor_(G,phase);
- // auto tmp = material_1(G,indicatorFunction_(x));
- // printmatrix(std::cout, material_1(G,indicatorFunction_(x)), "material_1(G,indicatorFunction_(x))", "--");
- // printmatrix(std::cout, elasticityTensor_(G,phase), "elasticityTensor_(G,phase)", "--");
- // return tmp;
- // return material_1(G,indicatorFunction_(x));
- }
- int applyIndicatorFunction(const Domain& x) const
- {
- // return indicatorFunction_(x);
- return indicatorFunction_material_1(x);
- }
+ // MatrixRT ElasticityTensorGlobal(const MatrixRT& G, const Domain& x) const //global Version (takes global coordinates)
+ // {
+ // return MAT(G,x);
+ // }
+
+
+ // //--- apply elasticityTensor_ to input Matrix G at position x
+ // MatrixRT applyElasticityTensor(const MatrixRT& G, const Domain& x) const //global Version (takes global coordinates)
+ // {
+ // // Python::Module module = Python::import("material");
+ // // auto indicatorFunctionTest = Python::make_function<double>(module.get("indicatorFunction"));
+ // // return elasticityTensor_(G,indicatorFunctionTest(x));
+ // // auto IFunction = localFunction(indicatorFunction_);
+ // // auto IFunction = this->getIndicatorFunction();
+ // // auto tmp = Dune::Functions::makeGridViewFunction(indicatorFunction_material_1, gridView_);
+ // // auto tmpLocal = LocalF
+ // // return elasticityTensor_(G,IFunction(x));
+ // // return elasticityTensor_(G,localIndicatorFunction_(x));
+ // return elasticityTensor_(G,indicatorFunction_(x));
+ // // return elasticityTensor_(G,indicatorFunctionGVF_(x));
+ // // int phase = indicatorFunction_(x);
+ // // auto tmp = elasticityTensor_(G,phase);
+ // // auto tmp = material_1(G,indicatorFunction_(x));
+ // // printmatrix(std::cout, material_1(G,indicatorFunction_(x)), "material_1(G,indicatorFunction_(x))", "--");
+ // // printmatrix(std::cout, elasticityTensor_(G,phase), "elasticityTensor_(G,phase)", "--");
+ // // return tmp;
+ // // return material_1(G,indicatorFunction_(x));
+ // }
+
//////////////////////////////////////////////////////////////////////////////
@@ -330,7 +434,6 @@ public:
// }
-
// -----------------------------------------------------------------
// --- write material (grid)functions to VTK
void write_materialFunctions()
@@ -340,24 +443,25 @@ public:
return;
};
-
///////////////////////////////
- // getter
+ // Getter
///////////////////////////////
ParameterTree getParameterSet() const {return parameterSet_;}
// Func2Tensor getElasticityTensor() const {return elasticityTensor_;}
// Func2TensorParam getElasticityTensor() const {return elasticityTensor_;}
Func2TensorPhase getElasticityTensor() const {return elasticityTensor_;}
+ MatrixPhase getPrestrainFunction() const {return prestrain_;}
auto getLocalIndicatorFunction() const {return localIndicatorFunction_;} //get as localFunction
auto getIndicatorFunctionGVF() const {return indicatorFunctionGVF_;} //get as GridViewFunction
auto getIndicatorFunction() const {return indicatorFunction_;}
+
// shared_ptr<Func2int> getIndicatorFunction(){return make_shared<Func2int>(indicatorFunction_);}
// auto getIndicatorFunction(){return make_shared<Func2int>(indicatorFunction_);}
// shared_ptr<Func2TensorParam> getElasticityTensor(){return make_shared<Func2TensorParam>(elasticityTensor_);}
- //TODO getLameParameters() .. Throw Exception if isotropic_ = false!
+ //getLameParameters() .. Throw Exception if isotropic_ = false!
@@ -365,6 +469,4 @@ public:
}; // end class
-
-
#endif
diff --git a/inputs/cellsolver.parset b/inputs/cellsolver.parset
index 46f1ffae0eba37b1b1cac55ac4c8f108e3382992..46160e703402550e1d8a8bbb172f5ec601b61a61 100644
--- a/inputs/cellsolver.parset
+++ b/inputs/cellsolver.parset
@@ -1,5 +1,4 @@
# --- Parameter File as Input for 'Cell-Problem'
-#
# NOTE: define variables without whitespaces in between! i.e. : gamma=1.0 instead of gamma = 1.0
# since otherwise these cant be read from other Files!
# --------------------------------------------------------
@@ -13,7 +12,7 @@
outputPath=/home/klaus/Desktop/Dune-Testing/dune-microstructure/outputs
-# --- DEBUG Option:
+# --- DEBUG (Output) Option:
print_debug = true #(default=false)
@@ -25,83 +24,36 @@ print_debug = true #(default=false)
## {start,finish} computes on all grid from 2^(start) to 2^finish refinement
#----------------------------------------------------
-numLevels= 2 4
-#numLevels = 0 0 # computes all levels from first to second entry
-#numLevels = 2 2 # computes all levels from first to second entry
-#numLevels = 1 3 # computes all levels from first to second entry
-#numLevels = 4 4 # computes all levels from first to second entry
-#numLevels = 5 5 # computes all levels from first to second entry
-#numLevels = 6 6 # computes all levels from first to second entry
-#numLevels = 1 6
-
-
-
-########################################################################################
-
-# --- Choose scale ratio gamma:
-gamma=1.0
+numLevels= 2 2 # computes all levels from first to second entry
#############################################
# Material / Prestrain parameters and ratios
#############################################
-#-- stiffness-ratio (ratio between material parameters mu1 & mu2 .... beta = 1.0 corresponds to homogeneous case)
-beta=2.0
-
-$--- strength of prestrain
-rho1=1.0
-
-#--- prestrain-contrast (ratio between prestrain parameters rho1 & rho2)
-alpha=8.0
-
-
-#materialFunction = "homogeneous"
-materialFunction = "material_1"
+# --- Choose material definition:
+#materialFunction = "two_phase_material_1"
+#materialFunction = "two_phase_material_2"
+materialFunction = "material_neukamm"
#materialFunction = "material_2"
+#materialFunction = "homogeneous"
+
+# --- Choose scale ratio gamma:
+gamma=1.0
-#--- Lame-Parameters
-mu1=80.0
-lambda1=80.0
-# better: pass material parameters as a vector
-# Poisson ratio nu = 1/4 => lambda = 0.4*mu
-#mu=1.2 1.2 1
-#lambda=0.48 0.48 0.4
-#rho=1.0 0
mu=80 80 60
lambda=80 80 25
-
-#mu=1
-#lambda=4
-
-
-
-# ---volume fraction (default value = 1.0/4.0)
+rho1=1.0
+beta=3.0
theta=0.25
-#theta = 0.3
-#theta = 0.75
-#theta = 0.125
-#theta = 0.5
-
-
+material_prestrain_imp= "material_neukamm"
#--- choose composite-Implementation:
#geometryFunctionPath = /home/stefan/DUNE/dune-microstructure/geometries
geometryFunctionPath = /home/klaus/Desktop/Dune-Testing/dune-microstructure/geometries
-material_prestrain_imp= "material_neukamm" #(Python-indicator-function with same name determines material phases)
-#material_prestrain_imp= "two_phase_material_2" #(Python-indicator-function with same name determines material phases)
-#material_prestrain_imp= "two_phase_material_3" #(Python-indicator-function with same name determines material phases)
-#material_prestrain_imp= "parametrized_Laminate"
-#material_prestrain_imp= "homogeneous"
-#material_prestrain_imp= "analytical_Example"
-#material_prestrain_imp="isotropic_bilayer"
-#material_prestrain_imp= "circle_fiber" #TEST
-#material_prestrain_imp= "matrix_material_circles"
-#material_prestrain_imp= "matrix_material_squares"
-#nF = 10 # Number of Fibers (default = 3)
-#rF = 0.1 # Fiber-Radius (default = 0.5*(width/(2.0*nF)) //half of the max-fiber-radius mrF = (width/(2.0*nF)) )
+
diff --git a/src/Cell-Problem-New.cc b/src/Cell-Problem-New.cc
index fdbf71640e52166182508a0f886fa11116449dda..22fe2eaf1c788753a83a9618d7f7cf7cc368adbf 100644
--- a/src/Cell-Problem-New.cc
+++ b/src/Cell-Problem-New.cc
@@ -201,9 +201,11 @@ int main(int argc, char *argv[])
///////////////////////////////////
// Get Prestrain/Parameters
///////////////////////////////////
- auto prestrainImp = PrestrainImp<dim>(); //NEW
+ auto prestrainImp = PrestrainImp<dim>();
auto B_Term = prestrainImp.getPrestrain(parameterSet);
+
+
log << "----- Input Parameters -----: " << std::endl;
log << "alpha: " << alpha << std::endl;
log << "gamma: " << gamma << std::endl;
@@ -585,6 +587,8 @@ int main(int argc, char *argv[])
if(print_debug)
correctorComputer.checkSymmetry();
+ // auto B_Term = material_.getPrestrainFunction();
+
//--- compute effective quantities
auto effectiveQuantitiesComputer = EffectiveQuantitiesComputer(correctorComputer,B_Term,material_);