Backup before elasticityTensor modification

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){

dune/microstructure/prestrainedMaterial.hh

+#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

geometries/material.py

+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]]

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

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;
 /*