From 274ab1abf0c6499177115bf6e6cd7cdae438de48 Mon Sep 17 00:00:00 2001
From: Klaus <klaus.boehnlein@tu-dresden.de>
Date: Mon, 29 Aug 2022 10:27:00 +0200
Subject: [PATCH] Add more material definitions
---
dune/microstructure/CorrectorComputer.hh | 20 ++++----
dune/microstructure/prestrainedMaterial.hh | 59 +++++++++++++++++++++-
inputs/cellsolver.parset | 8 ++-
src/Cell-Problem-New.cc | 31 +++++++-----
4 files changed, 90 insertions(+), 28 deletions(-)
diff --git a/dune/microstructure/CorrectorComputer.hh b/dune/microstructure/CorrectorComputer.hh
index 323ddbc8..f4b52c9c 100644
--- a/dune/microstructure/CorrectorComputer.hh
+++ b/dune/microstructure/CorrectorComputer.hh
@@ -350,14 +350,14 @@ public:
// double energyDensity= scalarProduct(etmp,defGradientV);
// double energyDensity= scalarProduct(material_.applyElasticityTensor(defGradientU,element.geometry().global(quadPos)),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(material_.applyElasticityTensor(defGradientU,element.geometry().global(quadPos)),sym(defGradientV));
// double energyDensity= scalarProduct(material_.applyElasticityTensor(defGradientU,quadPos),sym(defGradientV));
// double energyDensity= scalarProduct(material_.applyElasticityTensorLocal(defGradientU,quadPos),defGradientV);
- double energyDensity = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), defGradientU, defGradientV);
+ // double energyDensity = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), defGradientU, defGradientV);
// double energyDensity = linearizedStVenantKirchhoffDensity(mu(element.geometry().global(quadPos)), lambda(element.geometry().global(quadPos)), defGradientU, defGradientV); //TEST
// double energyDensity = generalizedDensity(mu(quadPos), lambda(quadPos), defGradientU, defGradientV); // also works..
@@ -371,13 +371,13 @@ public:
{
// double energyDensityGphi= scalarProduct(material_.applyElasticityTensor(basisContainer[m],element.geometry().global(quadPos)),defGradientV);
- // double energyDensityGphi= scalarProduct(elasticityTensor(basisContainer[m],indicatorFunction(element.geometry().global(quadPos))),sym(defGradientV));
+ double energyDensityGphi= scalarProduct(elasticityTensor(basisContainer[m],indicatorFunction(element.geometry().global(quadPos))),sym(defGradientV));
// double energyDensityGphi= scalarProduct(elasticityTensor(basisContainer[m],indicatorFunction(quadPos)),sym(defGradientV));
// double energyDensityGphi= scalarProduct(material_.applyElasticityTensor(basisContainer[m],element.geometry().global(quadPos)),sym(defGradientV));
// double energyDensityGphi= scalarProduct(material_.applyElasticityTensor(basisContainer[m],quadPos),sym(defGradientV));
// double energyDensityGphi= scalarProduct(elasticityTensor(basisContainer[m],element.geometry().global(quadPos)),defGradientV);
// double energyDensityGphi= scalarProduct(material_.applyElasticityTensorLocal(basisContainer[m],quadPos),defGradientV);
- double energyDensityGphi = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), basisContainer[m], defGradientV);
+ // double energyDensityGphi = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), basisContainer[m], defGradientV);
// double energyDensityGphi = linearizedStVenantKirchhoffDensity(mu(element.geometry().global(quadPos)), lambda(element.geometry().global(quadPos)), basisContainer[m], defGradientV); //TEST
auto value = energyDensityGphi * quadPoint.weight() * integrationElement;
elementMatrix[row][localPhiOffset+m] += value;
@@ -399,7 +399,7 @@ public:
// std::cout << "mu(quadPos): " << mu(quadPos) << std::endl;
// std::cout << "lambda(quadPos): " << lambda(quadPos) << std::endl;
- // 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(material_.applyElasticityTensor(basisContainer[m],element.geometry().global(quadPos)),basisContainer[n]);
// double energyDensityGG= scalarProduct(elasticityTensor(basisContainer[m],indicatorFunction(quadPos)),sym(basisContainer[n]));
@@ -407,7 +407,7 @@ public:
// double energyDensityGG= scalarProduct(material_.applyElasticityTensor(basisContainer[m],quadPos),sym(basisContainer[n]));
// double energyDensityGG= scalarProduct(material_.applyElasticityTensorLocal(basisContainer[m],quadPos),basisContainer[n]);
- double energyDensityGG = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), basisContainer[m], basisContainer[n]);
+ // double energyDensityGG = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), basisContainer[m], basisContainer[n]);
// double energyDensityGG = linearizedStVenantKirchhoffDensity(mu(element.geometry().global(quadPos)), lambda(element.geometry().global(quadPos)), basisContainer[m], basisContainer[n]); //TEST
elementMatrix[localPhiOffset+m][localPhiOffset+n] += energyDensityGG * quadPoint.weight() * integrationElement; // += !!!!! (Fixed-Bug)
@@ -520,14 +520,14 @@ public:
defGradientV = crossSectionDirectionScaling((1.0/gamma_),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(material_.applyElasticityTensor((-1.0)*forceTerm(quadPos),element.geometry().global(quadPos)),sym(defGradientV));
// double energyDensity= scalarProduct(material_.applyElasticityTensor((-1.0)*forceTerm(quadPos),quadPos),sym(defGradientV));
// double energyDensity= scalarProduct(material_.applyElasticityTensor((-1.0)*forceTerm(quadPos),element.geometry().global(quadPos)),defGradientV);
// double energyDensity= scalarProduct(material_.applyElasticityTensorLocal((-1.0)*forceTerm(quadPos),quadPos),defGradientV);
- double energyDensity = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos),(-1.0)*forceTerm(quadPos), defGradientV );
+ // double energyDensity = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos),(-1.0)*forceTerm(quadPos), defGradientV );
// double energyDensity = linearizedStVenantKirchhoffDensity(mu(element.geometry().global(quadPos)), lambda(element.geometry().global(quadPos)),forceTerm(quadPos), defGradientV ); //TEST
// double energyDensity = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos),(-1.0)*forceTerm(quadPos), defGradientV ); //TEST
// double energyDensity = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos),forceTerm(element.geometry().global(quadPos)), defGradientV ); //TEST
@@ -539,13 +539,13 @@ public:
// "f*m"-part
for (size_t m=0; m<3; 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(material_.applyElasticityTensor((-1.0)*forceTerm(quadPos),element.geometry().global(quadPos)),sym(basisContainer[m]));
// double energyDensityfG = scalarProduct(material_.applyElasticityTensor((-1.0)*forceTerm(quadPos),quadPos),sym(basisContainer[m]));
// double energyDensityfG = scalarProduct(material_.applyElasticityTensor((-1.0)*forceTerm(quadPos),element.geometry().global(quadPos)),basisContainer[m]);
// double energyDensityfG = scalarProduct(material_.applyElasticityTensor((-1.0)*forceTerm(quadPos),quadPos),basisContainer[m]);
- double energyDensityfG = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), (-1.0)*forceTerm(quadPos),basisContainer[m] );
+ // double energyDensityfG = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), (-1.0)*forceTerm(quadPos),basisContainer[m] );
// double energyDensityfG = linearizedStVenantKirchhoffDensity(mu(element.geometry().global(quadPos)), lambda(element.geometry().global(quadPos)), forceTerm(quadPos),basisContainer[m] ); //TEST
// double energyDensityfG = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), (-1.0)*forceTerm(quadPos),basisContainer[m] ); //TEST
// double energyDensityfG = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), forceTerm(element.geometry().global(quadPos)),basisContainer[m] );//TEST
diff --git a/dune/microstructure/prestrainedMaterial.hh b/dune/microstructure/prestrainedMaterial.hh
index ac423af3..210d46ca 100644
--- a/dune/microstructure/prestrainedMaterial.hh
+++ b/dune/microstructure/prestrainedMaterial.hh
@@ -40,6 +40,7 @@ using std::make_shared;
else
return 0; //#Phase3
}
+
MatrixRT material_1(const MatrixRT& G, const int& phase)
{
const FieldVector<double,3> mu = {80.0, 80.0, 60.0};
@@ -53,6 +54,44 @@ using std::make_shared;
}
+ 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();
+ }
+
+
+ template<class Domain>
+ int indicatorFunction_material_2(const Domain& x)
+ {
+ double theta=0.25;
+ if(abs(x[0]) > (theta/2.0))
+ return 1; //#Phase1
+ 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> lambda = {80.0, 160.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();
+ }
+
+
+
+
@@ -236,6 +275,22 @@ public:
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_1<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 ");
}
@@ -252,8 +307,8 @@ public:
// auto tmp = Dune::Functions::makeGridViewFunction(indicatorFunction_material_1, gridView_);
// auto tmpLocal = LocalF
// return elasticityTensor_(G,IFunction(x));
- // return elasticityTensor_(G,indicatorFunction_(x));
- return elasticityTensor_(G,indicatorFunctionGVF_(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));
diff --git a/inputs/cellsolver.parset b/inputs/cellsolver.parset
index 4fac1d02..3dcbfe05 100644
--- a/inputs/cellsolver.parset
+++ b/inputs/cellsolver.parset
@@ -56,6 +56,10 @@ rho1=1.0
alpha=8.0
+#materialFunction = "homogeneous"
+#materialFunction = "material"
+materialFunction = "material_2"
+
#--- Lame-Parameters
mu1=80.0
lambda1=80.0
@@ -86,10 +90,10 @@ theta=0.25
#--- 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= "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= "parametrized_Laminate"
#material_prestrain_imp= "homogeneous"
#material_prestrain_imp= "analytical_Example"
#material_prestrain_imp="isotropic_bilayer"
diff --git a/src/Cell-Problem-New.cc b/src/Cell-Problem-New.cc
index 4b024380..537c9c44 100644
--- a/src/Cell-Problem-New.cc
+++ b/src/Cell-Problem-New.cc
@@ -306,29 +306,32 @@ int main(int argc, char *argv[])
// int Phases = parameterSet.get<int>("Phases", 3);
- std::string materialFunctionName = parameterSet.get<std::string>("materialFunction", "material");
- Python::Module module = Python::import(materialFunctionName);
+ // std::string materialFunctionName = parameterSet.get<std::string>("materialFunction", "material");
+ // Python::Module module = Python::import(materialFunctionName);
+
+
+
// auto indicatorFunction = Python::make_function<double>(module.get("f"));
// Func2Tensor indicatorFunction = Python::make_function<double>(module.get("f"));
// auto materialFunction_ = Python::make_function<double>(module.get("f"));
// auto materialFunction_ = Python::make_function<double>(module.get("f"));
// auto materialFunction_ = Python::make_function<FieldVector<double,2>>(module.get("f"));
- int Phases;
- module.get("Phases").toC<int>(Phases);
- std::cout << "Number of Phases used:" << Phases << std::endl;
+ // int Phases;
+ // module.get("Phases").toC<int>(Phases);
+ // std::cout << "Number of Phases used:" << Phases << std::endl;
// std::cout << typeid(mu_).name() << '\n';
//---- Get mu/lambda values (for isotropic material) from Material-file
- FieldVector<double,3> mu_(0);
- module.get("mu_").toC<FieldVector<double,3>>(mu_);
- printvector(std::cout, mu_ , "mu_", "--");
- FieldVector<double,3> lambda_(0);
- module.get("lambda_").toC<FieldVector<double,3>>(lambda_);
- printvector(std::cout, lambda_ , "lambda_", "--");
+ // FieldVector<double,3> mu_(0);
+ // module.get("mu_").toC<FieldVector<double,3>>(mu_);
+ // printvector(std::cout, mu_ , "mu_", "--");
+ // FieldVector<double,3> lambda_(0);
+ // module.get("lambda_").toC<FieldVector<double,3>>(lambda_);
+ // printvector(std::cout, lambda_ , "lambda_", "--");
//////////////////////////////////////////////////////////////////////////////////////////////////////////
@@ -341,13 +344,13 @@ int main(int argc, char *argv[])
// Func2TensorParam elasticityTensor_ = *material_.getElasticityTensor();
auto elasticityTensor_ = material_.getElasticityTensor();
- Func2TensorParam TestTensor = Python::make_function<MatrixRT>(module.get("H"));
+ // Func2TensorParam TestTensor = Python::make_function<MatrixRT>(module.get("H"));
// std::cout << "decltype(elasticityTensor_) " << decltype(elasticityTensor_) << std::endl;
- std::cout << "typeid(elasticityTensor).name() :" << typeid(elasticityTensor_).name() << '\n';
- std::cout << "typeid(TestTensor).name() :" << typeid(TestTensor).name() << '\n';
+ // std::cout << "typeid(elasticityTensor).name() :" << typeid(elasticityTensor_).name() << '\n';
+ // std::cout << "typeid(TestTensor).name() :" << typeid(TestTensor).name() << '\n';
using MatrixFunc = std::function< MatrixRT(const MatrixRT&) >;
// std::cout << "Import NOW:" << std::endl;
--
GitLab