From 505c696455849205a984ec7fda590d01710c08e2 Mon Sep 17 00:00:00 2001
From: Klaus <klaus.boehnlein@tu-dresden.de>
Date: Sat, 27 Aug 2022 12:55:50 +0200
Subject: [PATCH] Add/Update prestrainedMaterial.hh
---
dune/microstructure/prestrainedMaterial.hh | 84 +++++++++++++------
geometries/material.py | 19 ++++-
inputs/cellsolver.parset | 2 +-
src/Cell-Problem-New.cc | 96 +++++++++++++++++++++-
4 files changed, 172 insertions(+), 29 deletions(-)
diff --git a/dune/microstructure/prestrainedMaterial.hh b/dune/microstructure/prestrainedMaterial.hh
index 74c8dc61..df3a516d 100644
--- a/dune/microstructure/prestrainedMaterial.hh
+++ b/dune/microstructure/prestrainedMaterial.hh
@@ -8,6 +8,7 @@
#include <dune/fufem/dunepython.hh>
+
using namespace Dune;
using namespace MatrixOperations;
using std::pow;
@@ -16,43 +17,52 @@ using std::sqrt;
using std::sin;
using std::cos;
+using std::shared_ptr;
+using std::make_shared;
-prestrainedMaterial
+template <class GridView> // needed for GridViewFunctions
+class prestrainedMaterial
{
public:
static const int dimworld = 3; //GridView::dimensionworld;
- static const int dim = Basis::GridView::dimension; //const int dim = Domain::dimension;
-
+ static const int dim = 3; //const int dim = Domain::dimension;
- using CellGridType = YaspGrid< dim, EquidistantOffsetCoordinates< double, dim>>;
- using Domain = typename CellGridType::LeafGridView::template Codim<0>::Geometry::GlobalCoordinate;
+ // using CellGridType = YaspGrid< dim, EquidistantOffsetCoordinates< double, dim>>;
+ // using Domain = typename CellGridType::LeafGridView::template Codim<0>::Geometry::GlobalCoordinate;
+ using Domain = typename GridView::template Codim<0>::Geometry::GlobalCoordinate;
using ScalarRT = FieldVector< double, 1>;
using VectorRT = FieldVector< double, dimworld>;
using MatrixRT = FieldMatrix< double, dimworld, dimworld>;
- using Func2Tensor = std::function< MatrixRT(const Domain&) >;
using FuncScalar = std::function< double(const Domain&) >;
+ using Func2Tensor = std::function< MatrixRT(const Domain&) >;
+ using Func2TensorParam = std::function< MatrixRT(const MatrixRT& ,const Domain&) >;
+
-protected
+protected:
+ const GridView& gridView_;
const ParameterTree& parameterSet_;
// const FieldVector<double , ...number of mu-Values/Phases> .. schwierig zur compile-time
- const FuncScalar mu_;
- const FuncScalar lambda_;
- double gamma_;
+ // const FuncScalar mu_;
+ // const FuncScalar lambda_;
+ // double gamma_;
- std::string materialFunctionName_
+ std::string materialFunctionName_;
// --- Number of material phases?
- const int phases_;
+ // const int phases_;
+
+ // Func2Tensor materialFunction_; //actually not needed??
- Func2Tensor materialFunction_;
+ // Func2Tensor elasticityTensor_;
+ Func2TensorParam elasticityTensor_;
// VectorCT x_1_, x_2_, x_3_; // (all) Corrector coefficient vectors
// VectorCT phi_1_, phi_2_, phi_3_; // Corrector phi_i coefficient vectors
@@ -87,22 +97,33 @@ public:
///////////////////////////////
// constructor
///////////////////////////////
- prestrainedMaterial( const ParameterTree& parameterSet) // string: "name of material"? // mu_(mu), muValues? müsste Anzahl Phasen bereits kennen..
+ prestrainedMaterial(const GridView gridView,
+ const ParameterTree& parameterSet) // string: "name of material"? // mu_(mu), muValues? müsste Anzahl Phasen bereits kennen..
+ : gridView_(gridView),
+ parameterSet_(parameterSet)
{
-
- std::string materialFunctionName_ = parameterSet.get<std::string>("materialFunction", "material");
+ std::string materialFunctionName_ = parameterSet.get<std::string>("materialFunction", "material");
Python::Module module = Python::import(materialFunctionName_);
- Func2Tensor materialFunction_ = Python::make_function<double>(module.get("f"));
+
+ elasticityTensor_ = Python::make_function<MatrixRT>(module.get("H"));
+
+ // Func2TensorParam elasticityTensor_ = Python::make_function<double>(module.get("L"));
+ // Func2Tensor materialFunction_ = Python::make_function<double>(module.get("f"));
+ // bool isotropic_ = true; // read from module File TODO
+ // Func2Tensor elasticityTensor_ = Python::make_function<double>(module.get("L"));
+ // Func2Tensor elasticityTensor_ = Python::make_function<MatrixRT>(module.get("L"));
+ }
- // if (parameterSet.get<bool>("stiffnessmatrix_cellproblem_to_csv"))
- // csvSystemMatrix();
- // if (parameterSet.get<bool>("rhs_cellproblem_to_csv"))
- // csvRHSs();
- // if (parameterSet.get<bool>("rhs_cellproblem_to_vtk"))
- // vtkLoads();
- }
+
+ MatrixRT applyElasticityTensor(const MatrixRT& G, const Domain& x)
+ {
+ //--- apply elasticityTensor_ to input Matrix G at position x
+ return elasticityTensor_(G,x);
+
+ }
+
// -----------------------------------------------------------------
@@ -122,6 +143,19 @@ public:
///////////////////////////////
ParameterTree getParameterSet() const {return parameterSet_;}
+
+ // Func2Tensor getElasticityTensor() const {return elasticityTensor_;}
+ Func2TensorParam getElasticityTensor() const {return elasticityTensor_;}
+
+
+ // shared_ptr<Func2TensorParam> getElasticityTensor(){return make_shared<Func2TensorParam>(elasticityTensor_);}
+
+
+ //TODO getLameParameters() .. Throw Exception if isotropic_ = false!
+
+
+
+
// shared_ptr<MatrixCT> getStiffnessMatrix(){return make_shared<MatrixCT>(stiffnessMatrix_);}
// shared_ptr<VectorCT> getLoad_alpha1(){return make_shared<VectorCT>(load_alpha1_);}
// shared_ptr<VectorCT> getLoad_alpha2(){return make_shared<VectorCT>(load_alpha2_);}
@@ -157,7 +191,7 @@ public:
-}
+}; // end class
diff --git a/geometries/material.py b/geometries/material.py
index 99fb2188..4aba62b8 100644
--- a/geometries/material.py
+++ b/geometries/material.py
@@ -98,7 +98,22 @@ def b3(x):
# --- elasticity tensor
+# def L(G,x):
+
+
def L(G):
- # input: symmetric matrix G
+ # input: symmetric matrix G, position x
# output: symmetric matrix LG
- return [[0, 0, 0], [0,0,0], [0,0,0]]
+ return [[1, 1, 1], [1, 1, 1],[1, 1, 1]]
+
+
+
+
+
+def H(G,x):
+ # input: symmetric matrix G, position x
+ # output: symmetric matrix LG
+ if (abs(x[0]) > 0.25):
+ return [[1, 1, 1], [1, 1, 1],[1, 1, 1]]
+ else:
+ return [[0, 0, 0], [0,0,0], [0,0,0]]
diff --git a/inputs/cellsolver.parset b/inputs/cellsolver.parset
index 4fac1d02..9f0192d8 100644
--- a/inputs/cellsolver.parset
+++ b/inputs/cellsolver.parset
@@ -25,7 +25,7 @@ print_debug = true #(default=false)
## {start,finish} computes on all grid from 2^(start) to 2^finish refinement
#----------------------------------------------------
-numLevels= 2 3
+numLevels= 2 2
#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
diff --git a/src/Cell-Problem-New.cc b/src/Cell-Problem-New.cc
index 9107f676..6e4b2cd6 100644
--- a/src/Cell-Problem-New.cc
+++ b/src/Cell-Problem-New.cc
@@ -49,6 +49,7 @@
#include <dune/microstructure/vtk_filler.hh> //TEST
#include <dune/microstructure/CorrectorComputer.hh>
#include <dune/microstructure/EffectiveQuantitiesComputer.hh>
+#include <dune/microstructure/prestrainedMaterial.hh>
#include <dune/solvers/solvers/umfpacksolver.hh> //TEST
#include <dune/istl/eigenvalue/test/matrixinfo.hh> // TEST: compute condition Number
@@ -108,6 +109,11 @@ auto equivalent = [](const FieldVector<double,3>& x, const FieldVector<double,3>
);
};
+
+
+// a function:
+int half(int x, int y) {return x/2+y/2;}
+
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
int main(int argc, char *argv[])
@@ -299,7 +305,6 @@ int main(int argc, char *argv[])
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"));
@@ -324,10 +329,92 @@ int main(int argc, char *argv[])
printvector(std::cout, lambda_ , "lambda_", "--");
+ //////////////////////////////////////////////////////////////////////////////////////////////////////////
+ // TESTING PRESTRAINEDMATERIAL.HH:
+ using Func2TensorParam = std::function< MatrixRT(const MatrixRT& ,const Domain&) >;
+
+ auto material_ = prestrainedMaterial(gridView_CE,parameterSet);
+ // Func2Tensor elasticityTensor = material_.getElasticityTensor();
+ // auto elasticityTensor = material_.getElasticityTensor();
+ // Func2TensorParam elasticityTensor_ = *material_.getElasticityTensor();
+ auto elasticityTensor_ = material_.getElasticityTensor();
+
+ 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';
+
+ // using MatrixFunc = std::function< MatrixRT(const MatrixRT&) >;
+ // using MatrixDomainFunc = std::function< MatrixRT(const MatrixRT&,const Domain&)>;
+ // // MatrixFunc elasticityTensor = Python::make_function<MatrixRT>(module.get("L"));
+ // MatrixDomainFunc TestTensor = Python::make_function<MatrixRT>(module.get("H"));
+
+ // auto elasticityTensorGVF = Dune::Functions::makeGridViewFunction(elasticityTensor , Basis_CE.gridView());
+ // auto localElasticityTensor = localFunction(elasticityTensorGVF);
+
+ // Func2Tensor forceTerm = [] (const Domain& x) {
+ // return MatrixRT{{1.0*x[2], 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}}; //TODO könnte hier sign übergeben?
+ // };
+
+ // auto loadGVF = Dune::Functions::makeGridViewFunction(forceTerm, Basis_CE.gridView());
+ // auto loadFunctional = localFunction(loadGVF);
+
+ MatrixRT G1_ {{1.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0, 0.0}};
+
+ // auto TestTensorGVF = Dune::Functions::makeGridViewFunction(TestTensor , Basis_CE.gridView());
+ // auto localTestTensor = localFunction(TestTensorGVF );
+
+
+ // printmatrix(std::cout, elasticityTensor(G1_), "elasticityTensor(G1_)", "--");
+ // auto temp = elasticityTensor(G1_);
+
+
+
+ for (const auto& element : elements(Basis_CE.gridView()))
+ {
+
+ int orderQR = 2;
+ const auto& quad = QuadratureRules<double,dim>::rule(element.type(), orderQR);
+ for (const auto& quadPoint : quad)
+ {
+ const auto& quadPos = quadPoint.position();
+ // std::cout << "quadPos : " << quadPos << std::endl;
+ 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, temp2, "temp2", "--");
+ }
+ }
+ // for (auto&& vertex : vertices(gridView_CE))
+ // {
+ // std::cout << "vertex.geometry().corner(0):" << vertex.geometry().corner(0)<< std::endl;
+ // auto tmp = vertex.geometry().corner(0);
+ // auto temp = elasticityTensor(tmp);
+ // // std::cout << "materialFunction_(vertex.geometry().corner(0))", materialFunction_(vertex.geometry().corner(0)) << std::endl;
+ // // printmatrix(std::cout, localElasticityTensor(G1_,tmp), "localElasticityTensor(vertex.geometry().corner(0))", "--");
+ // }
+
+
+
+ std::function<int(int,int)> fn1 = half;
+ std::cout << "fn1(60,20): " << fn1(60,20) << '\n';
+
+
+ // std::cout << typeid(elasticityTensorGVF).name() << '\n';
+ // std::cout << typeid(localElasticityTensor).name() << '\n';
+
+
// ParameterTree parameterSet_2;
// ParameterTreeParser::readINITree(geometryFunctionPath + "/"+ materialFunctionName + ".py", parameterSet_2);
@@ -347,9 +434,16 @@ int main(int argc, char *argv[])
// }
// std::cout << "materialFunction_({0.0,0.0,0.0})", materialFunction_({0.0,0.0,0.0}) << std::endl;
+
+
+ // --------------------------------------------------------------
+
//TODO// Phasen anhand von Mu bestimmen?
//TODO: DUNE_THROW(Exception, "Inconsistent choice of interpolation method"); if number of Phases != mu/lambda parameters
+ //FÜR L GARNICHT NÖTIG DENN RÜCKGABETYPE IS IMMER MATRIXRT!?!:
+ // BEi materialfunction (isotopic) reicht auch FieldVector<double,2> für lambda/mu
+
switch (Phases)
{
case 1: //homogeneous material
--
GitLab