Add/Update prestrainedMaterial.hh

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
 
 
 

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

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

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