Add Test against parametrized_Laminate (analytical results)

test/parametrized_laminate.py

+import math
+import ctypes
+import os
+import sys
+#---------------------------------------------------------------
+
+#Parameters used:
+theta = 0.25;
+mu_1 = 1.0;
+theta_mu = 2.0;
+mu_2 = theta_mu * mu_1;
+rho_1 = 1.0;
+theta_rho = 2.0;
+rho_2 = rho_1*theta_rho;
+
+
+#--- define indicator function
+def indicatorFunction(x):
+    theta=0.25
+    factor=1
+    if (abs(x[0]) < (theta/2) and x[2] < 0 ):
+        return 1    #Phase1
+    elif (abs(x[0]) > (theta/2) and x[2] > 0 ):
+        return 2    #Phase2
+    elif (abs(x[0]) < (theta/2) and x[2] > 0 ):
+        return 3    #Phase3
+    else :
+        return 4    #Phase4
+
+
+########### Options for material phases: #################################
+#     1. "isotropic"     2. "orthotropic"      3. "transversely_isotropic"   4. "general_anisotropic"
+#########################################################################
+## Notation - Parameter input :
+# isotropic (Lame parameters) : [mu , lambda]
+#         orthotropic         : [E1,E2,E3,G12,G23,G31,nu12,nu13,nu23]   # see https://en.wikipedia.org/wiki/Poisson%27s_ratio with x=1,y=2,z=3
+# transversely_isotropic      : [E1,E2,G12,nu12,nu23]
+# general_anisotropic         : full compliance matrix C
+######################################################################
+
+# --- Number of material phases
+Phases=4
+
+#--- Define different material phases:
+
+#- PHASE 1
+phase1_type="isotropic"
+materialParameters_phase1 = [2.0, 0]   
+
+#- PHASE 2
+phase2_type="isotropic"
+materialParameters_phase2 = [1.0, 0]   
+ 
+
+#- PHASE 3
+phase3_type="isotropic"
+materialParameters_phase3 = [2.0, 0]
+
+#- PHASE 4
+phase4_type="isotropic"
+materialParameters_phase4 = [1.0, 0]
+
+
+#--- define prestrain function for each phase
+# (also works with non-constant values)
+def prestrain_phase1(x):
+    return [[2, 0, 0], [0,2,0], [0,0,2]]
+
+def prestrain_phase2(x):
+    return [[1, 0, 0], [0,1,0], [0,0,1]]
+
+def prestrain_phase3(x):
+    return [[0, 0, 0], [0,0,0], [0,0,0]]
+
+def prestrain_phase4(x):
+    return [[0, 0, 0], [0,0,0], [0,0,0]]
+

test/parametrizedlaminate.parset

+# --- Choose material definition:
+materialFunction = "parametrized_laminate"
+
+# --- Choose scale ratio gamma:
+gamma=1.0
+
+#############################################
+#  Solver Type: #1: CG - SOLVER , #2: GMRES - SOLVER, #3: QR - SOLVER (default), #4: UMFPACK - SOLVER
+#############################################
+#Solvertype = 2        # recommended to use iterative solver (e.g GMRES) for finer grid-levels
+#Solver_verbosity = 0  #(default = 2)  degree of information for solver output
\ No newline at end of file

test/parametrizedlaminatetest.cc

@@ -11,12 +11,10 @@
 #include <dune/common/float_cmp.hh>
 #include <dune/common/math.hh>
 
-
 #include <dune/geometry/quadraturerules.hh>
 
 #include <dune/grid/uggrid.hh>
 #include <dune/grid/yaspgrid.hh>
-// #include <dune/grid/utility/structuredgridfactory.hh> //TEST
 #include <dune/grid/io/file/vtk/subsamplingvtkwriter.hh>
 
 #include <dune/istl/matrix.hh>
@@ -41,37 +39,216 @@
 #include <dune/functions/gridfunctions/gridviewfunction.hh>
 
 #include <dune/common/fvector.hh>
-#include <dune/common/fmatrix.hh> 
+#include <dune/common/fmatrix.hh>
 
-// #include <dune/microstructure/prestrain_material_geometry.hh>
 #include <dune/microstructure/matrix_operations.hh>
-// #include <dune/microstructure/vtk_filler.hh>    //TEST
-#include <dune/microstructure/CorrectorComputer.hh>    
-#include <dune/microstructure/EffectiveQuantitiesComputer.hh>  
-#include <dune/microstructure/prestrainedMaterial.hh>  
+#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 
+#include <dune/solvers/solvers/umfpacksolver.hh>
+#include <dune/istl/eigenvalue/test/matrixinfo.hh>
 
-// #include <dune/fufem/discretizationerror.hh>
 #include <dune/fufem/dunepython.hh>
 
-// #include <dune/fufem/functions/virtualgridfunction.hh> //TEST 
-
-// #include <boost/multiprecision/cpp_dec_float.hpp>
 #include <any>
 #include <variant>
 #include <string>
-#include <iomanip>   // needed when working with relative paths e.g. from python-scripts
+#include <iomanip>
+#include <filesystem>
 
 using namespace Dune;
 using namespace MatrixOperations;
 
-#include <iostream>
-
-int main(){
-
-
-    std::cout << "small little program." << std::endl;
 
-}
+/**
\ No newline at end of file
+ * @brief In the special case of a parametrized laminate (Lemma 4.5) in
+ *        [Böhnlein,Neukamm,Padilla-Garza,Sander - A homogenized bending theory for prestrained plates]
+ *        some of the effective quantities can be computed analytically.
+ *        This test compares the analytical quantities with numerical results.
+ */
+
+
+//////////////////////////////////////////////////
+//   Infrastructure for handling periodicity
+//////////////////////////////////////////////////
+// Check whether two points are equal on R/Z x R/Z x R
+auto equivalent = [](const FieldVector<double,3>& x, const FieldVector<double,3>& y)
+                  {
+                    return ( (FloatCmp::eq(x[0],y[0]) or FloatCmp::eq(x[0]+1,y[0]) or FloatCmp::eq(x[0]-1,y[0]))
+                             and (FloatCmp::eq(x[1],y[1]) or FloatCmp::eq(x[1]+1,y[1]) or FloatCmp::eq(x[1]-1,y[1]))
+                             and (FloatCmp::eq(x[2],y[2]))
+                             );
+                  };
+
+
+int main(int argc, char *argv[])
+{
+  MPIHelper::instance(argc, argv);
+
+  Dune::Timer globalTimer;
+  //--- setup Log-File 
+  std::fstream log;
+
+  std::cout << "Current path is " << std::filesystem::current_path() << '\n';
+  std::filesystem::path file_path = (__FILE__);
+  std::cout<< "File path: " << file_path<<std::endl;
+  std::cout << "dir_path: " << file_path.parent_path() << std::endl;
+  std::string dir_path  = file_path.parent_path();
+
+  ParameterTree parameterSet;
+  if (argc < 2)
+    ParameterTreeParser::readINITree(dir_path + "/parsets/parametrizedlaminate.parset", parameterSet);
+  else
+  {
+    ParameterTreeParser::readINITree(argv[1], parameterSet);
+    ParameterTreeParser::readOptions(argc, argv, parameterSet);
+  }
+
+  //--- Start Python interpreter
+  Python::start();
+  Python::Reference main = Python::import("__main__");
+  Python::run("import math");
+  Python::runStream()
+    << std::endl << "import sys"
+    << std::endl << "sys.path.append('" << dir_path  << "')"
+    << std::endl;
+
+
+  constexpr int dim = 3;
+  constexpr int dimWorld = 3;
+  ///////////////////////////////////
+  // Generate the grid
+  ///////////////////////////////////
+  // --- Corrector Problem Domain (-1/2,1/2)^3:
+  FieldVector<double,dim> lower({-1.0/2.0, -1.0/2.0, -1.0/2.0});
+  FieldVector<double,dim> upper({1.0/2.0, 1.0/2.0, 1.0/2.0});
+
+  int gridLevel = 3;
+
+  std::array<int, dim> nElements = {(int)std::pow(2,gridLevel) ,(int)std::pow(2,gridLevel) ,(int)std::pow(2,gridLevel)};
+  std::cout << "Number of Grid-Elements in each direction: " << nElements << std::endl;
+
+  using CellGridType = YaspGrid<dim, EquidistantOffsetCoordinates<double, dim> >;
+  CellGridType grid_CE(lower,upper,nElements);
+  using GridView = CellGridType::LeafGridView;
+  const GridView gridView_CE = grid_CE.leafGridView();
+
+  //--- Choose a finite element space for Cell Problem
+  using namespace Functions::BasisFactory;
+  Functions::BasisFactory::Experimental::PeriodicIndexSet periodicIndices;
+
+  //--- get PeriodicIndices for periodicBasis (Don't do the following in real life: It has quadratic run-time in the number of vertices.)
+  for (const auto& v1 : vertices(gridView_CE))
+    for (const auto& v2 : vertices(gridView_CE))
+      if (equivalent(v1.geometry().corner(0), v2.geometry().corner(0)))
+      {
+        periodicIndices.unifyIndexPair({gridView_CE.indexSet().index(v1)}, {gridView_CE.indexSet().index(v2)});
+      }
+
+  //--- setup first order periodic Lagrange-Basis
+  auto Basis_CE = makeBasis(
+    gridView_CE,
+    power<dim>(
+      Functions::BasisFactory::Experimental::periodic(lagrange<1>(), periodicIndices),
+      flatLexicographic()
+      ));
+
+  ///////////////////////////////////
+  //  Create prestrained material object
+  ///////////////////////////////////
+  auto material_ = prestrainedMaterial(gridView_CE,parameterSet);
+
+  // --- Get scale ratio
+  double gamma = parameterSet.get<double>("gamma",1.0);
+
+  //------------------------------------------------------------------------------------------------
+  //--- Compute Correctors
+  auto correctorComputer = CorrectorComputer(Basis_CE, material_, gamma, log, parameterSet);
+  correctorComputer.solve();
+
+
+  //--- Compute effective quantities
+  auto effectiveQuantitiesComputer = EffectiveQuantitiesComputer(correctorComputer,material_);
+  effectiveQuantitiesComputer.computeEffectiveQuantities();
+
+  //--- get effective quantities
+  auto Qeff = effectiveQuantitiesComputer.getQeff();
+  auto Beff = effectiveQuantitiesComputer.getBeff();
+//   printmatrix(std::cout, Qeff, "Matrix Qeff", "--");
+//   printvector(std::cout, Beff, "Beff", "--");
+
+  /**
+   * @brief Compute analytical effective quantitites
+   *
+   */
+  //Parameters
+  double theta = 0.25;
+  double mu_1 = 1.0;
+  double theta_mu = 2.0;
+  double mu_2 = theta_mu * mu_1;
+  double rho_1 = 1.0;
+  double theta_rho = 2.0;
+  double rho_2 = rho_1*theta_rho;
+
+  //Effective quantities
+  double q_1 = mu_1 * (theta_mu / (6*(theta + (1 - theta)*theta_mu)));
+  double q_2 = (mu_1/6.0) * ((1-theta) + theta*theta_mu);
+  double Beff_1 = (3.0*rho_1/2.0) * (1-(theta*(1+theta_rho)));
+  double Beff_2 = (3.0*rho_1/2.0) * ((1-theta*(1+theta_mu*theta_rho))/(1-theta+theta*theta_mu));
+  double Beff_3 = 0.0;
+  //  std::cout << "q_1: " << q_1 << std::endl;
+  //  std::cout << "q_2: " << q_2 << std::endl;
+  //  std::cout << "Beff_1: " << Beff_1 << std::endl;
+  //  std::cout << "Beff_2: " << Beff_2 << std::endl;
+  //  std::cout << "Beff_3: " << Beff_3 << std::endl;
+
+  /**
+   * @brief  Compare numerical and analytical results.
+   *
+   */
+  if( std::abs(q_1 - Qeff[0][0]) > 1e-2)
+  {
+    std::cerr << std::setprecision(9);
+    std::cerr << "Computed Qeff[0][0] is " << Qeff[0][0] << " but " << q_1 << " was expected!" << std::endl;
+
+    return 1;
+  }
+  if( std::abs(q_2 - Qeff[1][1]) > 1e-2)
+  {
+    std::cerr << std::setprecision(9);
+    std::cerr << "Computed Qeff[1][1] is " << Qeff[1][1] << " but " << q_2 << " was expected!" << std::endl;
+
+    return 1;
+  }
+  if( std::abs(Beff_1 - Beff[0]) > 1e-2)
+  {
+    std::cerr << std::setprecision(9);
+    std::cerr << "Computed Beff[0] is " << Beff[0] << " but " << Beff_1 << " was expected!" << std::endl;
+
+    return 1;
+  }
+  if( std::abs(Beff_2 - Beff[1]) > 1e-2)
+  {
+    std::cerr << std::setprecision(9);
+    std::cerr << "Computed Beff[1] is " << Beff[1] << " but " << Beff_2 << " was expected!" << std::endl;
+
+    return 1;
+  }
+  if( std::abs(Beff_3 - Beff[2]) > 1e-2)
+  {
+    std::cerr << std::setprecision(9);
+    std::cerr << "Computed Beff[2] is " << Beff[2] << " but " << Beff_3 << " was expected!" << std::endl;
+    return 1;
+  }
+  if( (Qeff[2][2] < Qeff[0][0]) || (Qeff[2][2] > Qeff[1][1]) )
+  {
+    std::cerr << std::setprecision(9);
+    std::cerr << "Computed Qeff[2][2] is " << Qeff[2][2] << " does not lie between  Qeff[0][0]:"
+              << Qeff[0][0] << " and Qeff[1][1]:" << Qeff[1][1] << " as was expected!" << std::endl;
+    return 1;
+  }
+
+  std::cout << "Total time elapsed: " << globalTimer.elapsed() << std::endl;
+  std::cout << "Test parametrizedlaminatetest passed." << std::endl;
+}