From 854d6a0752c95b06d889297df9a2680d6076c3a4 Mon Sep 17 00:00:00 2001
From: Klaus <klaus.boehnlein@tu-dresden.de>
Date: Sun, 28 Aug 2022 13:43:32 +0200
Subject: [PATCH] Set ElasticityTensor at prestrainedMaterial.hh
---
dune/microstructure/prestrainedMaterial.hh | 134 +++++++++++++++++++--
src/Cell-Problem-New.cc | 12 +-
2 files changed, 131 insertions(+), 15 deletions(-)
diff --git a/dune/microstructure/prestrainedMaterial.hh b/dune/microstructure/prestrainedMaterial.hh
index e255a6b2..274cf421 100644
--- a/dune/microstructure/prestrainedMaterial.hh
+++ b/dune/microstructure/prestrainedMaterial.hh
@@ -24,6 +24,26 @@ using std::shared_ptr;
using std::make_shared;
+// constexpr unsigned int str2int(const std::string* str, int h = 0)
+// {
+// return !str[h] ? 5381 : (str2int(str, h+1) * 33) ^ str[h];
+// }
+
+
+ MatrixRT material_1(const MatrixRT& G, const int& phase)
+ {
+ FieldVector<double,3> mu = {1.0, 2.0, 3.0};
+ FieldVector<double,3> lambda = {1.0 ,2.0 , 5.0};
+
+ if (phase == 0)
+ return 2.0 * mu[0] * sym(G) + lambda[0] * trace(sym(G)) * Id();
+ else
+ return 2.0 * mu[1] * sym(G) + lambda[2] * trace(sym(G)) * Id();
+
+ }
+
+
+
template <class GridView> // needed for GridViewFunctions
class prestrainedMaterial
@@ -41,8 +61,10 @@ public:
using VectorRT = FieldVector< double, dimworld>;
using MatrixRT = FieldMatrix< double, dimworld, dimworld>;
using FuncScalar = std::function< double(const Domain&) >;
+ using Func2int = std::function< int(const Domain&) >;
using Func2Tensor = std::function< MatrixRT(const Domain&) >;
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&) >;
@@ -59,7 +81,6 @@ protected:
-
// const FieldVector<double , ...number of mu-Values/Phases> .. schwierig zur compile-time
// const FuncScalar mu_;
@@ -67,6 +88,7 @@ protected:
// double gamma_;
std::string materialFunctionName_;
+ // std::string materialFunctionName_ = parameterSet_.get<std::string>("materialFunction", "material");
// --- Number of material phases?
// const int phases_;
@@ -74,11 +96,13 @@ protected:
// Func2Tensor materialFunction_; //actually not needed??
// Func2Tensor elasticityTensor_;
- Func2TensorParam elasticityTensor_;
+ // Func2TensorParam elasticityTensor_;
+ Func2TensorPhase elasticityTensor_;
// FuncScalar indicatorFunction_;
- GridViewFunction<double(const Domain&), GridView> indicatorFunction_;
+ // GridViewFunction<double(const Domain&), GridView> indicatorFunction_;
+ GridViewFunction<int(const Domain&), GridView> indicatorFunction_;
// static const auto indicatorFunction_;
// VectorCT x_1_, x_2_, x_3_; // (all) Corrector coefficient vectors
@@ -116,17 +140,26 @@ public:
///////////////////////////////
prestrainedMaterial(const GridView gridView,
const ParameterTree& parameterSet) // string: "name of material"? // mu_(mu), muValues? müsste Anzahl Phasen bereits kennen..
- : gridView_(gridView),
- parameterSet_(parameterSet)
+ : gridView_(gridView),
+ parameterSet_(parameterSet)
{
std::string materialFunctionName_ = parameterSet.get<std::string>("materialFunction", "material");
+
+
+ std::cout << "materialFunctionName_ : " << materialFunctionName_ << std::endl;
+
Python::Module module = Python::import(materialFunctionName_);
- elasticityTensor_ = Python::make_function<MatrixRT>(module.get("L"));
+ // elasticityTensor_ = Python::make_function<MatrixRT>(module.get("L"));
+
+ // elasticityTensor_ = setupElasticityTensor(materialFunctionName_);
+ setupElasticityTensor(materialFunctionName_);
+ // elasticityTensor_ = material_1;
+ // elasticityTensor_ = setupElasticityTensor();
// module.get("Phases").toC<int>(Phases_);
- auto indicatorFunction = Python::make_function<double>(module.get("indicatorFunction"));
+ auto indicatorFunction = Python::make_function<int>(module.get("indicatorFunction"));
indicatorFunction_ = Dune::Functions::makeGridViewFunction(indicatorFunction , gridView_);
@@ -148,14 +181,92 @@ public:
-
+
+ // static MatrixRT material_1(const MatrixRT& G, const int& phase)
+ // {
+ // FieldVector<double,3> mu = {1.0, 2.0, 3.0};
+ // FieldVector<double,3> lambda = {1.0 ,2.0 , 5.0};
+
+ // if (phase == 0)
+ // return 2.0 * mu[0] * sym(G) + lambda[0] * trace(sym(G)) * Id();
+ // else
+ // return 2.0 * mu[1] * sym(G) + lambda[2] * trace(sym(G)) * Id();
+
+ // }
+
+
+//---function that determines elasticity Tensor
+ // auto setupElasticityTensor(const std::string name) // cant use materialFunctionName_ here!?
+ // {
+ // if(name == "material")
+ // {
+ // return material_1;
+ // }
+ // else
+ // DUNE_THROW(Exception, "There exists no material in materialDefinitions.hh with this name ");
+ // }
+
+
+
+ //---function that determines elasticity Tensor
+ void setupElasticityTensor(const std::string name) // cant use materialFunctionName_ here!?
+ {
+ if(name == "material")
+ {
+ elasticityTensor_ = material_1;
+ }
+ else
+ DUNE_THROW(Exception, "There exists no material in materialDefinitions.hh with this name ");
+ }
+
+ // }
+ // void setupElasticityTensor()
+ // {
+ // if (materialFunctionName_=="material")
+ // {
+ // elasticityTensor_ = &material_1;
+ // std::cout << "typeid(elasticityTensor_).name() :" << typeid(elasticityTensor_).name() << '\n';
+ // // std::cout << " type_name<decltype(elasticityTensor_)>() " << type_name<decltype(elasticityTensor_)>() << '\n';
+ // }
+
+ // return;
+
+ // }
+
+ //--- apply elasticityTensor_ to input Matrix G at position x
MatrixRT applyElasticityTensor(const MatrixRT& G, const Domain& x) const
{
- //--- apply elasticityTensor_ to input Matrix G at position x
- return elasticityTensor_(G,x);
+
+ // Python::Module module = Python::import("material");
+ // auto indicatorFunctionTest = Python::make_function<double>(module.get("indicatorFunction"));
+ //
+ // return elasticityTensor_(G,indicatorFunctionTest(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));
}
+
+
+
+//////////////////////////////////////////////////////////////////////////////
+
+
+
+
+ // MatrixRT applyElasticityTensor(const MatrixRT& G, const Domain& x) const
+ // {
+ // //--- apply elasticityTensor_ to input Matrix G at position x
+ // return elasticityTensor_(G,x);
+
+ // }
+
MatrixRT applyElasticityTensorLocal(const MatrixRT& G, const Domain& x) const
{
//--- apply elasticityTensor_ to input Matrix G at position x (local coordinates)
@@ -191,7 +302,8 @@ public:
// Func2Tensor getElasticityTensor() const {return elasticityTensor_;}
- Func2TensorParam getElasticityTensor() const {return elasticityTensor_;}
+ // Func2TensorParam getElasticityTensor() const {return elasticityTensor_;}
+ Func2TensorPhase getElasticityTensor() const {return elasticityTensor_;}
diff --git a/src/Cell-Problem-New.cc b/src/Cell-Problem-New.cc
index 72a3744f..a98ecfbc 100644
--- a/src/Cell-Problem-New.cc
+++ b/src/Cell-Problem-New.cc
@@ -411,11 +411,12 @@ int main(int argc, char *argv[])
// std::cout << "quadPos : " << quadPos << std::endl;
- auto temp = TestTensor(G1_, element.geometry().global(quadPos));
- auto temp2 = elasticityTensor_(G1_, element.geometry().global(quadPos));
+ // auto temp = TestTensor(G1_, element.geometry().global(quadPos));
+ // auto temp2 = elasticityTensor_(G1_, element.geometry().global(quadPos));
// std::cout << "material_.applyElasticityTensor:" << std::endl;
auto tmp3 = material_.applyElasticityTensor(G1_, element.geometry().global(quadPos));
- // printmatrix(std::cout, tmp3, "tmp3", "--");
+ // auto tmp3 = material_.applyElasticityTensor(G1_, quadPos);
+ printmatrix(std::cout, tmp3, "tmp3", "--");
}
}
@@ -556,6 +557,9 @@ int main(int argc, char *argv[])
// };
+ }
+/*
+
//------------------------------------------------------------------------------------------------
//--- compute Correctors
// auto correctorComputer = CorrectorComputer(Basis_CE, muTerm, lambdaTerm, gamma, log, parameterSet);
@@ -783,6 +787,6 @@ int main(int argc, char *argv[])
std::cout << "Total time elapsed: " << globalTimer.elapsed() << std::endl;
-
+*/
}
--
GitLab