Set ElasticityTensor at prestrainedMaterial.hh

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_;}
 
 
 

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