Skip to content
Snippets Groups Projects
Commit 553d722e authored by Klaus Böhnlein's avatar Klaus Böhnlein
Browse files

backup

parent d8b9bf3a
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
dune/microstructure/CorrectorComputer.hh
+ 14
− 14
View file @ 553d722e
...@@ -270,8 +270,8 @@ public: ...@@ -270,8 +270,8 @@ public:
auto elasticityTensor = material_.getElasticityTensor(); auto elasticityTensor = material_.getElasticityTensor();
// auto indicatorFunction = material_.getIndicatorFunction(); // auto indicatorFunction = material_.getIndicatorFunction();
// auto localIndicatorFunction = material_.getLocalIndicatorFunction(); auto localIndicatorFunction = material_.getLocalIndicatorFunction();
// localIndicatorFunction.bind(element); localIndicatorFunction.bind(element);
// auto indicatorFunctionGVF = material_.getIndicatorFunctionGVF(); // auto indicatorFunctionGVF = material_.getIndicatorFunctionGVF();
// auto indicatorFunction = localFunction(indicatorFunctionGVF); // auto indicatorFunction = localFunction(indicatorFunctionGVF);
// indicatorFunction.bind(element); // indicatorFunction.bind(element);
...@@ -377,8 +377,8 @@ public: ...@@ -377,8 +377,8 @@ public:
// std::cout << "material_.applyIndicatorFunction(element.geometry().global(quadPos)): " << material_.applyIndicatorFunction(element.geometry().global(quadPos)) << std::endl; // 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; // 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_.ElasticityTensorGlobal(defGradientU,element.geometry().global(quadPos)),sym(defGradientV));
// double energyDensity= scalarProduct(material_.ElasticityTensorPhase(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,localIndicatorFunction(quadPos)),sym(defGradientV));
// double energyDensity= scalarProduct(elasticityTensor(defGradientU,indicatorFunction(element.geometry().global(quadPos))),sym(defGradientV)); // double energyDensity= scalarProduct(elasticityTensor(defGradientU,indicatorFunction(element.geometry().global(quadPos))),sym(defGradientV));
// double energyDensity= scalarProduct(elasticityTensor(defGradientU,indicatorFunction(quadPos)),sym(defGradientV)); // double energyDensity= scalarProduct(elasticityTensor(defGradientU,indicatorFunction(quadPos)),sym(defGradientV));
...@@ -412,8 +412,8 @@ public: ...@@ -412,8 +412,8 @@ public:
// } // }
double energyDensityGphi = scalarProduct(material_.ElasticityTensorGlobal(basisContainer[m],element.geometry().global(quadPos)),sym(defGradientV)); // 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_.ElasticityTensor(basisContainer[m],localIndicatorFunction(quadPos)),sym(defGradientV));
// double energyDensityGphi= scalarProduct(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 << "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 << "linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), basisContainer[m], defGradientV)" << linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), basisContainer[m], defGradientV) << std::endl;
...@@ -456,8 +456,8 @@ public: ...@@ -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; // 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_.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_.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],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],indicatorFunction(element.geometry().global(quadPos))),sym(basisContainer[n]));
// double energyDensityGG= scalarProduct(elasticityTensor(basisContainer[m],element.geometry().global(quadPos)),basisContainer[n]); // double energyDensityGG= scalarProduct(elasticityTensor(basisContainer[m],element.geometry().global(quadPos)),basisContainer[n]);
...@@ -506,8 +506,8 @@ public: ...@@ -506,8 +506,8 @@ public:
auto elasticityTensor = material_.getElasticityTensor(); auto elasticityTensor = material_.getElasticityTensor();
// auto indicatorFunction = material_.getIndicatorFunction(); // auto indicatorFunction = material_.getIndicatorFunction();
// auto localIndicatorFunction = material_.getLocalIndicatorFunction(); auto localIndicatorFunction = material_.getLocalIndicatorFunction();
// localIndicatorFunction.bind(element); localIndicatorFunction.bind(element);
// // auto indicatorFunction = material_.getLocalIndicatorFunction(); // // auto indicatorFunction = material_.getLocalIndicatorFunction();
// auto indicatorFunctionGVF = material_.getIndicatorFunctionGVF(); // auto indicatorFunctionGVF = material_.getIndicatorFunctionGVF();
// auto indicatorFunction = localFunction(indicatorFunctionGVF); // auto indicatorFunction = localFunction(indicatorFunctionGVF);
...@@ -585,8 +585,8 @@ public: ...@@ -585,8 +585,8 @@ public:
defGradientV = crossSectionDirectionScaling((1.0/gamma_),defGradientV); 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_.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_.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),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(element.geometry().global(quadPos))),sym(defGradientV));
// double energyDensity= scalarProduct(elasticityTensor((-1.0)*forceTerm(quadPos),indicatorFunction(quadPos)),sym(defGradientV)); // double energyDensity= scalarProduct(elasticityTensor((-1.0)*forceTerm(quadPos),indicatorFunction(quadPos)),sym(defGradientV));
...@@ -608,8 +608,8 @@ public: ...@@ -608,8 +608,8 @@ public:
// "f*m"-part // "f*m"-part
for (size_t m=0; m<3; m++) 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_.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_.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),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(element.geometry().global(quadPos))),sym(basisContainer[m]));
// double energyDensityfG= scalarProduct(elasticityTensor((-1.0)*forceTerm(quadPos),indicatorFunction(quadPos)),sym(basisContainer[m])); // double energyDensityfG= scalarProduct(elasticityTensor((-1.0)*forceTerm(quadPos),indicatorFunction(quadPos)),sym(basisContainer[m]));
......
dune/microstructure/prestrain_material_geometry.hh
+ 5
− 2
View file @ 553d722e
...@@ -54,7 +54,8 @@ public: ...@@ -54,7 +54,8 @@ public:
} }
else if (imp == "material_neukamm"){ else if (imp == "material_neukamm"){
// std::array<double,2> mu = parameters.get<std::array<double,2>>("mu", {1.0,3.0}); // 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(parameters.get<std::string>("material_neukamm"));
Python::Module module = Python::import("material_neukamm"); Python::Module module = Python::import("material_neukamm");
...@@ -495,7 +496,9 @@ public: ...@@ -495,7 +496,9 @@ public:
} }
if (imp == "material_neukamm"){ if (imp == "material_neukamm"){
// std::array<double,2> lambda = parameters.get<std::array<double,2>>("lambda", {1.0,3.0}); // 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"); Python::Module module = Python::import("material_neukamm");
auto indicatorFunction = Python::make_function<double>(module.get("f")); auto indicatorFunction = Python::make_function<double>(module.get("f"));
......
dune/microstructure/prestrainedMaterial.hh
+ 235
− 133
View file @ 553d722e
...@@ -27,50 +27,124 @@ using std::make_shared; ...@@ -27,50 +27,124 @@ using std::make_shared;
// ----------------------------------------------------------------------------------
template<class Domain> 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; double theta=0.25;
if (x[0] <(-0.5+theta) and x[2]<(-0.5+theta)) if (x[0] <(-0.5+theta) and x[2]<(-0.5+theta))
return 1; //#Phase1 return 1; //#Phase1
else if (x[1]<(-0.5+theta) and x[2]>(0.5-theta)) else if (x[1]<(-0.5+theta) and x[2]>(0.5-theta))
return 2; //#Phase2 return 2; //#Phase2
else else
return 0; //#Phase3 return 3; //#Phase3
} }
MatrixRT material_neukamm(const MatrixRT& G, const int& phase)
MatrixRT material_1(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}; 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) 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();
if (phase == 2) 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();
}
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 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> 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) MatrixRT two_phase_material_2(const MatrixRT& G, const int& phase)
{ {
const FieldVector<double,1> mu = {80.0}; const FieldVector<double,2> mu = {5.0, 15.0};
const FieldVector<double,1> lambda = {80.0}; const FieldVector<double,2> lambda = {6.0, 8.0};
return 2.0 * mu[0] * sym(G) + lambda[0] * trace(sym(G)) * Id(); 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> template<class Domain>
int indicatorFunction_material_2(const Domain& x) int indicatorFunction_material_2(const Domain& x)
{ {
...@@ -80,16 +154,64 @@ using std::make_shared; ...@@ -80,16 +154,64 @@ using std::make_shared;
else else
return 0; //#Phase2 return 0; //#Phase2
} }
MatrixRT material_2(const MatrixRT& G, const int& phase) 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}; const FieldVector<double,2> lambda = {80.0, 160.0};
if (phase == 1) 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 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; ...@@ -107,8 +229,6 @@ using std::make_shared;
// return 2.0 * mu[2] * sym(G) + lambda[2] * trace(sym(G)) * Id(); // return 2.0 * mu[2] * sym(G) + lambda[2] * trace(sym(G)) * Id();
// } // }
template <class GridView> // needed for GridViewFunctions template <class GridView> // needed for GridViewFunctions
class prestrainedMaterial class prestrainedMaterial
{ {
...@@ -129,38 +249,25 @@ public: ...@@ -129,38 +249,25 @@ public:
using Func2TensorParam = std::function< MatrixRT(const MatrixRT& ,const Domain&) >; using Func2TensorParam = std::function< MatrixRT(const MatrixRT& ,const Domain&) >;
using Func2TensorPhase = std::function< MatrixRT(const MatrixRT& ,const int&) >; using Func2TensorPhase = std::function< MatrixRT(const MatrixRT& ,const int&) >;
using MatrixFunc = std::function< MatrixRT(const MatrixRT&) >; using MatrixFunc = std::function< MatrixRT(const MatrixRT&) >;
using MatrixPhase = std::function< MatrixRT(const int&) >;
protected: protected:
const GridView& gridView_; const GridView& gridView_;
const ParameterTree& parameterSet_; const ParameterTree& parameterSet_;
std::string materialFunctionName_;
int Phases_; // MatrixFunc L1_;
// MatrixFunc L2_;
MatrixFunc L1_; // MatrixFunc L3_;
MatrixFunc L2_;
MatrixFunc L3_;
// const FieldVector<double , ...number of mu-Values/Phases> .. schwierig zur compile-time
// const FuncScalar mu_; // const FuncScalar mu_;
// const FuncScalar lambda_; // const FuncScalar lambda_;
// double gamma_;
std::string materialFunctionName_;
// std::string materialFunctionName_ = parameterSet_.get<std::string>("materialFunction", "material");
// --- Number of material phases?
// const int phases_; // const int phases_;
// Func2Tensor materialFunction_; //actually not needed??
// Func2Tensor elasticityTensor_; // Func2Tensor elasticityTensor_;
// Func2TensorParam elasticityTensor_; // Func2TensorParam elasticityTensor_;
Func2TensorPhase elasticityTensor_; Func2TensorPhase elasticityTensor_;
MatrixPhase prestrain_;
// FuncScalar indicatorFunction_;
Func2int indicatorFunction_; Func2int indicatorFunction_;
GridViewFunction<int(const Domain&), GridView> indicatorFunctionGVF_; GridViewFunction<int(const Domain&), GridView> indicatorFunctionGVF_;
...@@ -178,21 +285,17 @@ public: ...@@ -178,21 +285,17 @@ public:
std::string materialFunctionName_ = parameterSet.get<std::string>("materialFunction", "material_1"); std::string materialFunctionName_ = parameterSet.get<std::string>("materialFunction", "material_1");
std::cout << "materialFunctionName_ : " << materialFunctionName_ << std::endl; std::cout << "materialFunctionName_ : " << materialFunctionName_ << std::endl;
// Python::Module module = Python::import(materialFunctionName_); setupMaterial(materialFunctionName_);
// elasticityTensor_ = Python::make_function<MatrixRT>(module.get("L"));
// elasticityTensor_ = setupElasticityTensor(materialFunctionName_);
setupElasticityTensor(materialFunctionName_);
indicatorFunctionGVF_ = Dune::Functions::makeGridViewFunction(indicatorFunction_, gridView_); indicatorFunctionGVF_ = Dune::Functions::makeGridViewFunction(indicatorFunction_, gridView_);
localIndicatorFunction_ = localFunction(indicatorFunctionGVF_); 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")); // auto indicatorFunction = Python::make_function<int>(module.get("indicatorFunction"));
// indicatorFunction_ = Python::make_function<int>(module.get("indicatorFunction")); // indicatorFunction_ = Python::make_function<int>(module.get("indicatorFunction"));
// indicatorFunction_ = Dune::Functions::makeGridViewFunction(indicatorFunction , gridView_); // indicatorFunction_ = Dune::Functions::makeGridViewFunction(indicatorFunction , gridView_);
// module.get("Phases").toC<int>(Phases_); // module.get("Phases").toC<int>(Phases_);
// L1_ = Python::make_function<MatrixRT>(module.get("L1")); // L1_ = Python::make_function<MatrixRT>(module.get("L1"));
// L2_ = Python::make_function<MatrixRT>(module.get("L2")); // L2_ = Python::make_function<MatrixRT>(module.get("L2"));
...@@ -200,8 +303,52 @@ public: ...@@ -200,8 +303,52 @@ public:
// bool isotropic_ = true; // read from module File // bool isotropic_ = true; // read from module File
// auto Test = Dune::Functions::makeGridViewFunction(MAT<Domain> , gridView_); //not working // 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) // static int indicatorFunction_material_1(const Domain& x)
// { // {
...@@ -224,39 +371,6 @@ public: ...@@ -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 // MatrixRT localElasticityTensor(const MatrixRT& G, const Domain& x, Element& element) const //local Version ... muss binding öfters durchführen
// { // {
// localIndicatorFunction_.bind(*element); // localIndicatorFunction_.bind(*element);
...@@ -265,46 +379,36 @@ public: ...@@ -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 // MatrixRT ElasticityTensorGlobal(const MatrixRT& G, const Domain& x) const //global Version (takes global coordinates)
{ // {
// return indicatorFunction_(x); // return MAT(G,x);
return indicatorFunction_material_1(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: ...@@ -330,7 +434,6 @@ public:
// } // }
// ----------------------------------------------------------------- // -----------------------------------------------------------------
// --- write material (grid)functions to VTK // --- write material (grid)functions to VTK
void write_materialFunctions() void write_materialFunctions()
...@@ -340,24 +443,25 @@ public: ...@@ -340,24 +443,25 @@ public:
return; return;
}; };
/////////////////////////////// ///////////////////////////////
// getter // Getter
/////////////////////////////// ///////////////////////////////
ParameterTree getParameterSet() const {return parameterSet_;} ParameterTree getParameterSet() const {return parameterSet_;}
// Func2Tensor getElasticityTensor() const {return elasticityTensor_;} // Func2Tensor getElasticityTensor() const {return elasticityTensor_;}
// Func2TensorParam getElasticityTensor() const {return elasticityTensor_;} // Func2TensorParam getElasticityTensor() const {return elasticityTensor_;}
Func2TensorPhase getElasticityTensor() const {return elasticityTensor_;} Func2TensorPhase getElasticityTensor() const {return elasticityTensor_;}
MatrixPhase getPrestrainFunction() const {return prestrain_;}
auto getLocalIndicatorFunction() const {return localIndicatorFunction_;} //get as localFunction auto getLocalIndicatorFunction() const {return localIndicatorFunction_;} //get as localFunction
auto getIndicatorFunctionGVF() const {return indicatorFunctionGVF_;} //get as GridViewFunction auto getIndicatorFunctionGVF() const {return indicatorFunctionGVF_;} //get as GridViewFunction
auto getIndicatorFunction() const {return indicatorFunction_;} auto getIndicatorFunction() const {return indicatorFunction_;}
// shared_ptr<Func2int> getIndicatorFunction(){return make_shared<Func2int>(indicatorFunction_);} // shared_ptr<Func2int> getIndicatorFunction(){return make_shared<Func2int>(indicatorFunction_);}
// auto getIndicatorFunction(){return make_shared<Func2int>(indicatorFunction_);} // auto getIndicatorFunction(){return make_shared<Func2int>(indicatorFunction_);}
// shared_ptr<Func2TensorParam> getElasticityTensor(){return make_shared<Func2TensorParam>(elasticityTensor_);} // 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: ...@@ -365,6 +469,4 @@ public:
}; // end class }; // end class
#endif #endif
inputs/cellsolver.parset
+ 14
− 62
View file @ 553d722e
# --- Parameter File as Input for 'Cell-Problem' # --- Parameter File as Input for 'Cell-Problem'
#
# NOTE: define variables without whitespaces in between! i.e. : gamma=1.0 instead of gamma = 1.0 # 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! # since otherwise these cant be read from other Files!
# -------------------------------------------------------- # --------------------------------------------------------
...@@ -13,7 +12,7 @@ ...@@ -13,7 +12,7 @@
outputPath=/home/klaus/Desktop/Dune-Testing/dune-microstructure/outputs outputPath=/home/klaus/Desktop/Dune-Testing/dune-microstructure/outputs
# --- DEBUG Option: # --- DEBUG (Output) Option:
print_debug = true #(default=false) print_debug = true #(default=false)
...@@ -25,83 +24,36 @@ 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 ## {start,finish} computes on all grid from 2^(start) to 2^finish refinement
#---------------------------------------------------- #----------------------------------------------------
numLevels= 2 4 numLevels= 2 2 # computes all levels from first to second entry
#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
############################################# #############################################
# Material / Prestrain parameters and ratios # Material / Prestrain parameters and ratios
############################################# #############################################
#-- stiffness-ratio (ratio between material parameters mu1 & mu2 .... beta = 1.0 corresponds to homogeneous case) # --- Choose material definition:
beta=2.0 #materialFunction = "two_phase_material_1"
#materialFunction = "two_phase_material_2"
$--- strength of prestrain materialFunction = "material_neukamm"
rho1=1.0
#--- prestrain-contrast (ratio between prestrain parameters rho1 & rho2)
alpha=8.0
#materialFunction = "homogeneous"
materialFunction = "material_1"
#materialFunction = "material_2" #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 mu=80 80 60
lambda=80 80 25 lambda=80 80 25
rho1=1.0
#mu=1 beta=3.0
#lambda=4
# ---volume fraction (default value = 1.0/4.0)
theta=0.25 theta=0.25
#theta = 0.3 material_prestrain_imp= "material_neukamm"
#theta = 0.75
#theta = 0.125
#theta = 0.5
#--- choose composite-Implementation: #--- choose composite-Implementation:
#geometryFunctionPath = /home/stefan/DUNE/dune-microstructure/geometries #geometryFunctionPath = /home/stefan/DUNE/dune-microstructure/geometries
geometryFunctionPath = /home/klaus/Desktop/Dune-Testing/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)) )
......
src/Cell-Problem-New.cc
+ 5
− 1
View file @ 553d722e
...@@ -201,9 +201,11 @@ int main(int argc, char *argv[]) ...@@ -201,9 +201,11 @@ int main(int argc, char *argv[])
/////////////////////////////////// ///////////////////////////////////
// Get Prestrain/Parameters // Get Prestrain/Parameters
/////////////////////////////////// ///////////////////////////////////
auto prestrainImp = PrestrainImp<dim>(); //NEW auto prestrainImp = PrestrainImp<dim>();
auto B_Term = prestrainImp.getPrestrain(parameterSet); auto B_Term = prestrainImp.getPrestrain(parameterSet);
log << "----- Input Parameters -----: " << std::endl; log << "----- Input Parameters -----: " << std::endl;
log << "alpha: " << alpha << std::endl; log << "alpha: " << alpha << std::endl;
log << "gamma: " << gamma << std::endl; log << "gamma: " << gamma << std::endl;
...@@ -585,6 +587,8 @@ int main(int argc, char *argv[]) ...@@ -585,6 +587,8 @@ int main(int argc, char *argv[])
if(print_debug) if(print_debug)
correctorComputer.checkSymmetry(); correctorComputer.checkSymmetry();
// auto B_Term = material_.getPrestrainFunction();
//--- compute effective quantities //--- compute effective quantities
auto effectiveQuantitiesComputer = EffectiveQuantitiesComputer(correctorComputer,B_Term,material_); auto effectiveQuantitiesComputer = EffectiveQuantitiesComputer(correctorComputer,B_Term,material_);
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment