diff --git a/dune/microstructure/matrix_operations.hh b/dune/microstructure/matrix_operations.hh
index 1c89d36fbc2cf07abd97995fb319d17e32a561c9..ed4036047e548c434aaebeafabba51f87d2c5e39 100644
--- a/dune/microstructure/matrix_operations.hh
+++ b/dune/microstructure/matrix_operations.hh
@@ -110,7 +110,7 @@ namespace MatrixOperations {
{
auto t1 = 2.0 * mu * sym(E1) + lambda * trace(sym(E1)) * Id();
auto tmp1 = scalarProduct(t1,sym(E2));
-
+ // auto tmp1 = scalarProduct(t1,E2);
// auto t2 = 2.0 * mu * sym(E2) + lambda * trace(sym(E2)) * Id();
// auto tmp2 = scalarProduct(t2,sym(E1));
//
@@ -124,6 +124,14 @@ namespace MatrixOperations {
}
+ // static MatrixRT elasticityTensor()
+ // {
+
+
+
+
+ // }
+
// --- Generalization: Define Quadratic QuadraticForm
static double QuadraticForm(const double mu, const double lambda, const MatrixRT M){
diff --git a/dune/microstructure/prestrainedMaterial.hh b/dune/microstructure/prestrainedMaterial.hh
new file mode 100644
index 0000000000000000000000000000000000000000..74c8dc61d25fe901c42b5a6befc9ddd2b501039b
--- /dev/null
+++ b/dune/microstructure/prestrainedMaterial.hh
@@ -0,0 +1,165 @@
+#ifndef DUNE_MICROSTRUCTURE_PRESTRAINEDMATERIAL_HH
+#define DUNE_MICROSTRUCTURE_PRESTRAINEDMATERIAL_HH
+
+
+#include <dune/grid/uggrid.hh>
+#include <dune/grid/yaspgrid.hh>
+#include <dune/microstructure/matrix_operations.hh>
+
+#include <dune/fufem/dunepython.hh>
+
+using namespace Dune;
+using namespace MatrixOperations;
+using std::pow;
+using std::abs;
+using std::sqrt;
+using std::sin;
+using std::cos;
+
+
+
+
+prestrainedMaterial
+{
+
+public:
+ static const int dimworld = 3; //GridView::dimensionworld;
+ static const int dim = Basis::GridView::dimension; //const int dim = Domain::dimension;
+
+
+ using CellGridType = YaspGrid< dim, EquidistantOffsetCoordinates< double, dim>>;
+ using Domain = typename CellGridType::LeafGridView::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&) >;
+
+protected
+
+ const ParameterTree& parameterSet_;
+
+
+ // const FieldVector<double , ...number of mu-Values/Phases> .. schwierig zur compile-time
+
+ const FuncScalar mu_;
+ const FuncScalar lambda_;
+ double gamma_;
+
+ std::string materialFunctionName_
+
+ // --- Number of material phases?
+ const int phases_;
+
+ Func2Tensor materialFunction_;
+
+// VectorCT x_1_, x_2_, x_3_; // (all) Corrector coefficient vectors
+// VectorCT phi_1_, phi_2_, phi_3_; // Corrector phi_i coefficient vectors
+// FieldVector<double,3> m_1_, m_2_, m_3_; // Corrector m_i coefficient vectors
+
+// MatrixRT M1_, M2_, M3_; // (assembled) corrector matrices M_i
+// const std::array<MatrixRT*, 3 > mContainer = {&M1_ , &M2_, &M3_};
+// const std::array<VectorCT, 3> phiContainer = {phi_1_,phi_2_,phi_3_};
+
+ // ---- Basis for R_sym^{2x2}
+ MatrixRT G1_ {{1.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0, 0.0}};
+ MatrixRT G2_ {{0.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0, 0.0, 0.0}};
+ MatrixRT G3_ {{0.0, 1.0/sqrt(2.0), 0.0}, {1.0/sqrt(2.0), 0.0, 0.0}, {0.0, 0.0, 0.0}};
+ std::array<MatrixRT,3 > MatrixBasisContainer_ = {G1_, G2_, G3_};
+
+ Func2Tensor x3G_1_ = [] (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?
+ };
+
+ Func2Tensor x3G_2_ = [] (const Domain& x) {
+ return MatrixRT{{0.0, 0.0, 0.0}, {0.0, 1.0*x[2], 0.0}, {0.0, 0.0, 0.0}};
+ };
+
+ Func2Tensor x3G_3_ = [] (const Domain& x) {
+ return MatrixRT{{0.0, (1.0/sqrt(2.0))*x[2], 0.0}, {(1.0/sqrt(2.0))*x[2], 0.0, 0.0}, {0.0, 0.0, 0.0}};
+ };
+
+ const std::array<Func2Tensor, 3> x3MatrixBasisContainer_ = {x3G_1_, x3G_2_, x3G_3_};
+
+
+public:
+ ///////////////////////////////
+ // constructor
+ ///////////////////////////////
+ prestrainedMaterial( const ParameterTree& parameterSet) // string: "name of material"? // mu_(mu), muValues? müsste Anzahl Phasen bereits kennen..
+ {
+
+ std::string materialFunctionName_ = parameterSet.get<std::string>("materialFunction", "material");
+ Python::Module module = Python::import(materialFunctionName_);
+ Func2Tensor materialFunction_ = Python::make_function<double>(module.get("f"));
+
+
+
+ // 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();
+ }
+
+
+ // -----------------------------------------------------------------
+ // --- write material (grid)functions to VTK
+ void write_materialFunctions()
+ {
+
+
+ return;
+
+ };
+
+
+
+ ///////////////////////////////
+ // getter
+ ///////////////////////////////
+ ParameterTree getParameterSet() const {return parameterSet_;}
+
+ // 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_);}
+ // shared_ptr<VectorCT> getLoad_alpha3(){return make_shared<VectorCT>(load_alpha3_);}
+
+ // shared_ptr<FuncScalar> getMu(){return make_shared<FuncScalar>(mu_);}
+ // shared_ptr<FuncScalar> getLambda(){return make_shared<FuncScalar>(lambda_);}
+
+
+ // --- Get Correctors
+ // shared_ptr<VectorCT> getMcontainer(){return make_shared<VectorCT>(mContainer);}
+ // auto getMcontainer(){return make_shared<std::array<MatrixRT*, 3 > >(mContainer);}
+ // auto getMcontainer(){return mContainer;}
+ // shared_ptr<std::array<VectorCT, 3>> getPhicontainer(){return make_shared<std::array<VectorCT, 3>>(phiContainer);}
+
+
+ // // shared_ptr<std::array<VectorRT, 3>> getBasiscontainer(){return make_shared<std::array<VectorRT, 3>>(basisContainer_);}
+ // auto getMatrixBasiscontainer(){return make_shared<std::array<MatrixRT,3 >>(MatrixBasisContainer_);}
+ // // auto getx3MatrixBasiscontainer(){return make_shared<std::array<Func2Tensor, 3>>(x3MatrixBasisContainer_);}
+ // auto getx3MatrixBasiscontainer(){return x3MatrixBasisContainer_;}
+
+
+
+
+
+ // shared_ptr<VectorCT> getCorr_a(){return make_shared<VectorCT>(x_a_);}
+ // shared_ptr<VectorCT> getCorr_phi1(){return make_shared<VectorCT>(phi_1_);}
+ // shared_ptr<VectorCT> getCorr_phi2(){return make_shared<VectorCT>(phi_2_);}
+ // shared_ptr<VectorCT> getCorr_phi3(){return make_shared<VectorCT>(phi_3_);}
+
+
+
+
+
+
+}
+
+
+
+
+#endif
diff --git a/geometries/material.py b/geometries/material.py
new file mode 100644
index 0000000000000000000000000000000000000000..99fb21886c2b15606b7dedee5657c77d4bf21e90
--- /dev/null
+++ b/geometries/material.py
@@ -0,0 +1,104 @@
+import math
+
+
+
+Phases = 3
+
+lu = [1,2,3]
+mu_ = [3, 5, 6]
+lambda_ = [9, 7, 8]
+
+
+# lu = mu[0] mu[1] mu[2]
+
+
+
+#Indicator function that determines both phases
+# x[0] : y1-component
+# x[1] : y2-component
+# x[2] : x3-component
+# --- replace with your definition of indicatorFunction:
+###############
+# Wood
+###############
+def f(x):
+ theta=0.25
+ # mu_ = [3, 5, 6]
+ # lambda_ = [9, 7, 8]
+ # mu_ = 3 5 6
+ # lambda_ = 9 7 8
+
+ if ((abs(x[0]) < theta/2) and x[2]<0.25):
+ return [mu_[0], lambda_[0]] #latewood
+ # return 5 #latewood
+ elif ((abs(x[0]) > theta/2) and x[2]<0.25):
+ return [mu_[1], lambda_[1]] #latewood
+ # return 2
+ else :
+ return [mu_[2],lambda_[2]] #latewood #Phase3
+ # return 1
+
+
+
+#Workaround
+def muValue(x):
+ return mu_
+
+def lambdaValue(x):
+ return lambda_
+
+# def b1(x):
+# return [[.5, 0, 0], [0,1,0], [0,0,0]]
+
+# def b2(x):
+# return [[.4, 0, 0], [0,.4,0], [0,0,0]]
+
+# def b3(x):
+# return [[0, 0, 0], [0,0,0], [0,0,0]]
+
+
+
+
+###############
+# Cross
+###############
+# def f(x):
+# theta=0.25
+# factor=1
+# if (x[0] <-1/2+theta and x[2]<-1/2+theta):
+# return 1 #Phase1
+# elif (x[1]< -1/2+theta and x[2]>1/2-theta):
+# return 2 #Phase2
+# else :
+# return 0 #Phase3
+
+
+
+
+
+# def f(x):
+# # --- replace with your definition of indicatorFunction:
+# if (abs(x[0]) > 0.25):
+# return 1 #Phase1
+# else :
+# return 0 #Phase2
+
+def b1(x):
+ return [[1, 0, 0], [0,1,0], [0,0,1]]
+
+def b2(x):
+ return [[1, 0, 0], [0,1,0], [0,0,1]]
+
+def b3(x):
+ return [[0, 0, 0], [0,0,0], [0,0,0]]
+
+# mu=80 80 60
+
+# lambda=80 80 25
+
+
+# --- elasticity tensor
+def L(G):
+ # input: symmetric matrix G
+ # output: symmetric matrix LG
+ return [[0, 0, 0], [0,0,0], [0,0,0]]
diff --git a/inputs/cellsolver.parset b/inputs/cellsolver.parset
index 9f0192d8d07c83221b7cba98b10ee07a1da8d569..4fac1d0232c3f450664a966d2f0f1f04cecef87d 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 2
+numLevels= 2 3
#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 b7ae0127d8f60344f93fbd3223c1d3c7a5078487..9107f67674666b7c41d8480192d402bf2d44ef56 100644
--- a/src/Cell-Problem-New.cc
+++ b/src/Cell-Problem-New.cc
@@ -246,8 +246,9 @@ int main(int argc, char *argv[])
std::cout << "Host grid has " << gridView_CE.size(dim) << " vertices." << std::endl;
// //not needed
- // using MatrixRT = FieldMatrix< double, dimWorld, dimWorld>;
- // using Domain = GridView::Codim<0>::Geometry::GlobalCoordinate;
+ using MatrixRT = FieldMatrix< double, dimWorld, dimWorld>;
+ using Domain = GridView::Codim<0>::Geometry::GlobalCoordinate;
+ using Func2Tensor = std::function< MatrixRT(const Domain&) >;
// using Func2Tensor = std::function< MatrixRT(const Domain&) >;
// using VectorCT = BlockVector<FieldVector<double,1> >;
// using MatrixCT = BCRSMatrix<FieldMatrix<double,1,1> >;
@@ -292,25 +293,91 @@ int main(int argc, char *argv[])
-
+ //TEST
//Read from Parset...
- int Phases = parameterSet.get<int>("Phases", 3);
- std::string materialFunction = parameterSet.get<std::string>("materialFunction", "material");
+ // int Phases = parameterSet.get<int>("Phases", 3);
+
+
+ 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;
+
+
+ // std::cout << typeid(mu_).name() << '\n';
- Python::Module module = Python::import("material");
- auto indicatorFunction = Python::make_function<double>(module.get("f"));
- std::cout << "Phases:" << Phases << std::endl;
+ //---- 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_", "--");
+
+
+
+
+
+
+ // ParameterTree parameterSet_2;
+ // ParameterTreeParser::readINITree(geometryFunctionPath + "/"+ materialFunctionName + ".py", parameterSet_2);
+
+ // auto lu = parameterSet_2.get<FieldVector<double,3>>("lu", {1.0,3.0,2.0});
+ // std::cout <<"lu[1]: " << lu[1]<< std::endl;
+ // std::cout <<"lu: " << parameterSet_2.get<std::array<double,3>>("lu", {1.0,3.0,2.0}) << std::endl;
+
+ // auto mU_ = module.evaluate(parameterSet_2.get<std::string>("lu", "[1,2,3]"));
+ // std::cout << "typeid(mU_).name()" << typeid(mU_.operator()()).name() << '\n';
+
+
+ // for (auto&& vertex : vertices(gridView_CE))
+ // {
+ // std::cout << "vertex.geometry().corner(0):" << vertex.geometry().corner(0)<< std::endl;
+ // // std::cout << "materialFunction_(vertex.geometry().corner(0))", materialFunction_(vertex.geometry().corner(0)) << std::endl;
+ // printvector(std::cout, materialFunction_(vertex.geometry().corner(0)), "materialFunction_(vertex.geometry().corner(0))", "--");
+ // }
+ // 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
+ switch (Phases)
+ {
+ case 1: //homogeneous material
+ {
+ std::cout << "Phase - 1" << std::endl;
+ auto materialFunction_ = Python::make_function<FieldVector<double,2>>(module.get("f"));
+ break;
+ }
+ case 2:
+ {
+ std::cout << "Phase - 1" << std::endl;
+ auto materialFunction_ = Python::make_function<FieldVector<double,2>>(module.get("f"));
+ break;
+ }
+ case 3:
+ {
+ std::cout << "Phase - 3" << std::endl;
+ auto materialFunction_ = Python::make_function<FieldVector<double,2>>(module.get("f"));
+ break;
+ }
+ }
+
// switch (Phases)
// {
// case 1: //homogeneous material
// {
- // std::cout << "Phase - 1" << std::endl;
+ // std::cout << "Phases - 1" << std::endl;
// std::array<double,1> mu_ = parameterSet.get<std::array<double,1>>("mu", {1.0});
// Python::Module module = Python::import(materialFunction);
// auto indicatorFunction = Python::make_function<double>(module.get("f")); // get indicator function
@@ -319,12 +386,12 @@ int main(int argc, char *argv[])
// }
// case 2:
// {
- // std::cout << "Phase - 2" << std::endl;
+ // std::cout << "Phases - 2" << std::endl;
// std::array<double,2> mu_ = parameterSet.get<std::array<double,2>>("mu", {1.0,3.0});
// Python::Module module = Python::import(materialFunction);
// auto indicatorFunction = Python::make_function<double>(module.get("f")); // get indicator function
// auto muTerm = [mu_,indicatorFunction] (const Domain& x)
- // {
+ // {
// if (indicatorFunction(x) == 1)
// return mu_[0];
// else
@@ -334,7 +401,7 @@ int main(int argc, char *argv[])
// }
// case 3:
// {
- // std::cout << "Phase - 3" << std::endl;
+ // std::cout << "Phases - 3" << std::endl;
// std::array<double,3> mu_ = parameterSet.get<std::array<double,3>>("mu", {1.0,3.0, 5.0});
// Python::Module module = Python::import(materialFunction);
// auto indicatorFunction = Python::make_function<double>(module.get("f")); // get indicator function
@@ -352,6 +419,8 @@ int main(int argc, char *argv[])
// }
+
+
//TEST
// std::cout << "Test crossSectionDirectionScaling:" << std::endl;
/*