diff --git a/CMakeLists.txt b/CMakeLists.txt
index cc9322f017781ff734fe82ad5489092c76b4bc06..2ac234249ef7768010ee7c2f3cb3ef3077032e20 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -26,6 +26,7 @@ add_subdirectory(src)
add_subdirectory(dune)
add_subdirectory(doc)
add_subdirectory(cmake/modules)
+add_subdirectory(test)
# finalize the dune project, e.g. generating config.h etc.
diff --git a/dune/microstructure/matrix_operations.hh b/dune/microstructure/matrix_operations.hh
index 3cd76beddfeda08b1749c784afbe77ee6a9f40bb..ab0a59332e73849964923ba388c31dadff9c4a0e 100644
--- a/dune/microstructure/matrix_operations.hh
+++ b/dune/microstructure/matrix_operations.hh
@@ -116,7 +116,7 @@ namespace MatrixOperations {
{ 0, 0, 1.0}};
}
default:
- DUNE_THROW(Dune::Exception, " axis not feasible. rotationMatrix is only implemented for 3x3-matrices.");
+ DUNE_THROW(Dune::Exception, " axis not feasible. rotationMatrix is only implemented for 3x3-matrices. Choose between 0: x-axis, 1: y-axis, 2: z-axis");
}
}
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..f2e918562b1e2ea68456621c709640e767a8ec12
--- /dev/null
+++ b/test/CMakeLists.txt
@@ -0,0 +1,18 @@
+set(TESTS
+ parametrizedlaminatetest
+ orthotropicrotationtest)
+
+foreach(_test ${TESTS})
+dune_add_test(SOURCES ${_test}.cc)
+add_dune_pythonlibs_flags(${_test})
+target_compile_options(${_test} PRIVATE -g)
+endforeach()
+
+
+# Copy the parsets used for testing into the build dir
+# file(COPY parsets/parametrizedlaminate.parset DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/parsets)
+# file(COPY parsets/orthotropicrotation_1.parset DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/parsets)
+# file(COPY parsets/orthotropicrotation_2.parset DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/parsets)
+
+
+
\ No newline at end of file
diff --git a/test/orthotropicrotation_1.parset b/test/orthotropicrotation_1.parset
new file mode 100644
index 0000000000000000000000000000000000000000..f337f16b1aff991e24cd50ba300199eb382dad87
--- /dev/null
+++ b/test/orthotropicrotation_1.parset
@@ -0,0 +1,20 @@
+# --- Choose material definition:
+materialFunction = "orthotropicrotation_1"
+
+# --- Choose scale ratio gamma:
+gamma=1.0
+
+#############################################
+# Solver Type: #1: CG - SOLVER , #2: GMRES - SOLVER, #3: QR - SOLVER (default), #4: UMFPACK - SOLVER
+#############################################
+
+# Note: Running the test with the QR-Solver (3) on gitlab
+# fails due to completely wrong results.
+# WHY is currently unclear.
+
+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
+
+#print_conditionNumber = true
+#set_IntegralZero = true #(default = false)
+#set_oneBasisFunction_Zero = true #(default = false)
\ No newline at end of file
diff --git a/test/orthotropicrotation_1.py b/test/orthotropicrotation_1.py
new file mode 100644
index 0000000000000000000000000000000000000000..7d15aba6bc3ac5d02c1e560b3dac830f29ebe6c5
--- /dev/null
+++ b/test/orthotropicrotation_1.py
@@ -0,0 +1,53 @@
+import math
+import ctypes
+import os
+import sys
+import numpy as np
+#---------------------------------------------------------------
+
+#--- define indicator function
+def indicatorFunction(x):
+ if x[2] > 0.25:
+ return 1 #Phase1
+ elif (abs(x[2]) < 0.25):
+ return 2 #Phase2
+ else:
+ return 3 #Phase3
+
+########### 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=3
+
+#--- Define different material phases:
+#- PHASE 1
+phase1_type="orthotropic"
+materialParameters_phase1 = [11.2e3,630,1190,700,230,960,0.63,0.49,0.37] # Walnut parameters (values for compliance matrix) see [Dimwoodie; Timber its nature and behavior p.109]
+
+#- PHASE 2
+phase2_type="orthotropic"
+materialParameters_phase2 = [16.3e3,620,1100,910,190,1180,0.43,0.49,0.38] # Birch parameters (values for compliance matrix) see [Dimwoodie; Timber its nature and behavior p.109]
+
+#- PHASE 3
+phase3_type="orthotropic"
+materialParameters_phase3 = [10.7e3,430,710,620,23,500,0.51,0.38,0.31] # Norway spruce parameters (values for compliance matrix) see [Dimwoodie; Timber its nature and behavior p.109]
+
+
+#--- 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]]
diff --git a/test/orthotropicrotation_2.parset b/test/orthotropicrotation_2.parset
new file mode 100644
index 0000000000000000000000000000000000000000..4eabf4cdf0d377d75a79ac41eeda1bd30d6dd23a
--- /dev/null
+++ b/test/orthotropicrotation_2.parset
@@ -0,0 +1,20 @@
+# --- Choose material definition:
+materialFunction = "orthotropicrotation_2"
+
+# --- Choose scale ratio gamma:
+gamma=1.0
+
+#############################################
+# Solver Type: #1: CG - SOLVER , #2: GMRES - SOLVER, #3: QR - SOLVER (default), #4: UMFPACK - SOLVER
+#############################################
+
+# Note: Running the test with the QR-Solver (3) on gitlab
+# fails due to completely wrong results.
+# WHY is currently unclear.
+
+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
+
+#print_conditionNumber = true
+#set_IntegralZero = true #(default = false)
+#set_oneBasisFunction_Zero = true #(default = false)
\ No newline at end of file
diff --git a/test/orthotropicrotation_2.py b/test/orthotropicrotation_2.py
new file mode 100644
index 0000000000000000000000000000000000000000..9a325f52cbe46fd628ce3ea7ec793448655e188c
--- /dev/null
+++ b/test/orthotropicrotation_2.py
@@ -0,0 +1,59 @@
+import math
+import ctypes
+import os
+import sys
+import numpy as np
+#---------------------------------------------------------------
+
+#--- define indicator function
+def indicatorFunction(x):
+ if x[2] > 0.25:
+ return 1 #Phase1
+ elif (abs(x[2]) < 0.25):
+ return 2 #Phase2
+ else:
+ return 3 #Phase3
+
+########### 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=3
+
+# *Rotate material frame of each phase by 180 degrees (pi)
+
+#--- Define different material phases:
+#- PHASE 1
+phase1_type="orthotropic"
+materialParameters_phase1 = [11.2e3,630,1190,700,230,960,0.63,0.49,0.37] # Walnut parameters (values for compliance matrix) see [Dimwoodie; Timber its nature and behavior p.109]
+phase1_axis = 0
+phase1_angle = np.pi
+
+#- PHASE 2
+phase2_type="orthotropic"
+materialParameters_phase2 = [16.3e3,620,1100,910,190,1180,0.43,0.49,0.38] # Birch parameters (values for compliance matrix) see [Dimwoodie; Timber its nature and behavior p.109]
+phase2_axis = 1
+phase2_angle = np.pi
+#- PHASE 3
+phase3_type="orthotropic"
+materialParameters_phase3 = [10.7e3,430,710,620,23,500,0.51,0.38,0.31] # Norway spruce parameters (values for compliance matrix) see [Dimwoodie; Timber its nature and behavior p.109]
+phase3_axis = 2
+phase3_angle = np.pi
+
+#--- 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]]
diff --git a/test/orthotropicrotationtest.cc b/test/orthotropicrotationtest.cc
new file mode 100644
index 0000000000000000000000000000000000000000..962652fa4760272289c99f9f7c04cbd598c13946
--- /dev/null
+++ b/test/orthotropicrotationtest.cc
@@ -0,0 +1,283 @@
+#include <config.h>
+#include <array>
+#include <vector>
+#include <fstream>
+
+#include <iostream>
+#include <dune/common/indices.hh>
+#include <dune/common/bitsetvector.hh>
+#include <dune/common/parametertree.hh>
+#include <dune/common/parametertreeparser.hh>
+#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/io/file/vtk/subsamplingvtkwriter.hh>
+
+#include <dune/istl/matrix.hh>
+#include <dune/istl/bcrsmatrix.hh>
+#include <dune/istl/multitypeblockmatrix.hh>
+#include <dune/istl/multitypeblockvector.hh>
+#include <dune/istl/matrixindexset.hh>
+#include <dune/istl/solvers.hh>
+#include <dune/istl/spqr.hh>
+#include <dune/istl/preconditioners.hh>
+#include <dune/istl/io.hh>
+
+#include <dune/functions/functionspacebases/interpolate.hh>
+#include <dune/functions/backends/istlvectorbackend.hh>
+#include <dune/functions/functionspacebases/powerbasis.hh>
+#include <dune/functions/functionspacebases/compositebasis.hh>
+#include <dune/functions/functionspacebases/lagrangebasis.hh>
+#include <dune/functions/functionspacebases/periodicbasis.hh>
+#include <dune/functions/functionspacebases/subspacebasis.hh>
+#include <dune/functions/functionspacebases/boundarydofs.hh>
+#include <dune/functions/gridfunctions/discreteglobalbasisfunction.hh>
+#include <dune/functions/gridfunctions/gridviewfunction.hh>
+
+#include <dune/common/fvector.hh>
+#include <dune/common/fmatrix.hh>
+
+#include <dune/microstructure/matrix_operations.hh>
+#include <dune/microstructure/CorrectorComputer.hh>
+#include <dune/microstructure/EffectiveQuantitiesComputer.hh>
+#include <dune/microstructure/prestrainedMaterial.hh>
+
+#include <dune/solvers/solvers/umfpacksolver.hh>
+#include <dune/istl/eigenvalue/test/matrixinfo.hh>
+
+#include <dune/fufem/dunepython.hh>
+
+#include <any>
+#include <variant>
+#include <string>
+#include <iomanip>
+#include <filesystem>
+
+using namespace Dune;
+using namespace MatrixOperations;
+
+
+/**
+ * @brief Compute effective quantities for a trilayer setup with 3 orthotropic phases:
+ * Once with a canonical orthonormal material frame and once with a material frame
+ * that is rotated pi(180 degree) around different axes.
+ * Aim of the test is to check the correctness of the material frame transformation.
+ */
+
+
+//////////////////////////////////////////////////
+// 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 + "/orthotropicrotation_1.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 "rotated" effective quantitites
+ *
+ */
+// ParameterTreeParser::readINITree("../test/orthotropicrotation_2.parset", parameterSet);
+ ParameterTreeParser::readINITree(dir_path + "/orthotropicrotation_2.parset", parameterSet);
+ auto material_2 = prestrainedMaterial(gridView_CE,parameterSet);
+
+ //--- Compute Correctors
+ auto correctorComputer_2 = CorrectorComputer(Basis_CE, material_2, gamma, log, parameterSet);
+ correctorComputer_2.solve();
+
+ //--- Compute effective quantities
+ auto effectiveQuantitiesComputer_2 = EffectiveQuantitiesComputer(correctorComputer_2,material_2);
+ effectiveQuantitiesComputer_2.computeEffectiveQuantities();
+
+ //--- get effective quantities
+ auto Qeff_2 = effectiveQuantitiesComputer_2.getQeff();
+ auto Beff_2 = effectiveQuantitiesComputer_2.getBeff();
+
+ /**
+ * @brief Compare results.
+ *
+ */
+ if( std::abs(Qeff[0][0] - Qeff_2[0][0]) > 1e-6)
+ {
+ std::cerr << std::setprecision(9);
+ std::cerr << "Computed Qeff[0][0] is " << Qeff[0][0] << " but with rotated material frame we get : " << Qeff_2[0][0] << " !" << std::endl;
+ return 1;
+ }
+ if( std::abs(Qeff[1][0] - Qeff_2[1][0]) > 1e-6)
+ {
+ std::cerr << std::setprecision(9);
+ std::cerr << "Computed Qeff[1][0] is " << Qeff[1][0] << " but with rotated material frame we get : " << Qeff_2[1][0] << " !" << std::endl;
+ return 1;
+ }
+ if( std::abs(Qeff[2][0] - Qeff_2[2][0]) > 1e-6)
+ {
+ std::cerr << std::setprecision(9);
+ std::cerr << "Computed Qeff[2][0] is " << Qeff[2][0] << " but with rotated material frame we get : " << Qeff_2[2][0] << " !" << std::endl;
+ return 1;
+ }
+ if( std::abs(Qeff[0][1] - Qeff_2[0][1]) > 1e-6)
+ {
+ std::cerr << std::setprecision(9);
+ std::cerr << "Computed Qeff[0][1] is " << Qeff[0][1] << " but with rotated material frame we get : " << Qeff_2[0][1] << " !" << std::endl;
+ return 1;
+ }
+ if( std::abs(Qeff[1][1] - Qeff_2[1][1]) > 1e-6)
+ {
+ std::cerr << std::setprecision(9);
+ std::cerr << "Computed Qeff[1][1] is " << Qeff[1][1] << " but with rotated material frame we get : " << Qeff_2[1][1] << " !" << std::endl;
+ return 1;
+ }
+ if( std::abs(Qeff[2][1] - Qeff_2[2][1]) > 1e-6)
+ {
+ std::cerr << std::setprecision(9);
+ std::cerr << "Computed Qeff[2][1] is " << Qeff[2][1] << " but with rotated material frame we get : " << Qeff_2[2][1] << " !" << std::endl;
+ return 1;
+ }
+ if( std::abs(Qeff[0][2] - Qeff_2[0][2]) > 1e-6)
+ {
+ std::cerr << std::setprecision(9);
+ std::cerr << "Computed Qeff[0][2] is " << Qeff[0][2] << " but with rotated material frame we get : " << Qeff_2[0][2] << " !" << std::endl;
+ return 1;
+ }
+ if( std::abs(Qeff[1][2] - Qeff_2[1][2]) > 1e-6)
+ {
+ std::cerr << std::setprecision(9);
+ std::cerr << "Computed Qeff[1][2] is " << Qeff[1][2] << " but with rotated material frame we get : " << Qeff_2[1][2] << " !" << std::endl;
+ return 1;
+ }
+ if( std::abs(Qeff[2][2] - Qeff_2[2][2]) > 1e-6)
+ {
+ std::cerr << std::setprecision(9);
+ std::cerr << "Computed Qeff[2][2] is " << Qeff[2][2] << " but with rotated material frame we get : " << Qeff_2[2][2] << " !" << std::endl;
+ return 1;
+ }
+
+
+ if( std::abs(Beff[0] - Beff_2[0]) > 1e-6)
+ {
+ std::cerr << std::setprecision(9);
+ std::cerr << "Computed Beff[0] is " << Beff[0] << " but with rotated material frame we get : " << Beff_2[0] << " !" << std::endl;
+ return 1;
+ }
+ if( std::abs(Beff[1] - Beff_2[1]) > 1e-6)
+ {
+ std::cerr << std::setprecision(9);
+ std::cerr << "Computed Beff[1] is " << Beff[1] << " but with rotated material frame we get : " << Beff_2[1] << " !" << std::endl;
+ return 1;
+ }
+ if( std::abs(Beff[2] - Beff_2[2]) > 1e-6)
+ {
+ std::cerr << std::setprecision(9);
+ std::cerr << "Computed Beff[2] is " << Beff[2] << " but with rotated material frame we get : " << Beff_2[2] << " !" << std::endl;
+ return 1;
+ }
+
+ std::cout << "Total time elapsed: " << globalTimer.elapsed() << std::endl;
+ std::cout << "Test: orthotropicrotationtest passed." << std::endl;
+}
diff --git a/test/parametrized_laminate.py b/test/parametrized_laminate.py
new file mode 100644
index 0000000000000000000000000000000000000000..9007b436bcc2e51ca4a57929f46cd42014000605
--- /dev/null
+++ b/test/parametrized_laminate.py
@@ -0,0 +1,77 @@
+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]]
+
diff --git a/test/parametrizedlaminate.parset b/test/parametrizedlaminate.parset
new file mode 100644
index 0000000000000000000000000000000000000000..f02b8e746b934135082bba9aa5a3c913269570f9
--- /dev/null
+++ b/test/parametrizedlaminate.parset
@@ -0,0 +1,20 @@
+# --- 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
+#############################################
+
+# Note: Running the test with the QR-Solver (3) on gitlab
+# fails due to completely wrong results.
+# WHY is currently unclear.
+
+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
+
+#print_conditionNumber = true
+#set_IntegralZero = true #(default = false)
+#set_oneBasisFunction_Zero = true #(default = false)
\ No newline at end of file
diff --git a/test/parametrizedlaminatetest.cc b/test/parametrizedlaminatetest.cc
new file mode 100644
index 0000000000000000000000000000000000000000..636df2beee9289422265b92d3a4bf00e12b059b8
--- /dev/null
+++ b/test/parametrizedlaminatetest.cc
@@ -0,0 +1,274 @@
+#include <config.h>
+#include <array>
+#include <vector>
+#include <fstream>
+
+#include <iostream>
+#include <dune/common/indices.hh>
+#include <dune/common/bitsetvector.hh>
+#include <dune/common/parametertree.hh>
+#include <dune/common/parametertreeparser.hh>
+#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/io/file/vtk/subsamplingvtkwriter.hh>
+
+#include <dune/istl/matrix.hh>
+#include <dune/istl/bcrsmatrix.hh>
+#include <dune/istl/multitypeblockmatrix.hh>
+#include <dune/istl/multitypeblockvector.hh>
+#include <dune/istl/matrixindexset.hh>
+#include <dune/istl/solvers.hh>
+#include <dune/istl/spqr.hh>
+#include <dune/istl/preconditioners.hh>
+#include <dune/istl/io.hh>
+
+#include <dune/functions/functionspacebases/interpolate.hh>
+#include <dune/functions/backends/istlvectorbackend.hh>
+#include <dune/functions/functionspacebases/powerbasis.hh>
+#include <dune/functions/functionspacebases/compositebasis.hh>
+#include <dune/functions/functionspacebases/lagrangebasis.hh>
+#include <dune/functions/functionspacebases/periodicbasis.hh>
+#include <dune/functions/functionspacebases/subspacebasis.hh>
+#include <dune/functions/functionspacebases/boundarydofs.hh>
+#include <dune/functions/gridfunctions/discreteglobalbasisfunction.hh>
+#include <dune/functions/gridfunctions/gridviewfunction.hh>
+
+#include <dune/common/fvector.hh>
+#include <dune/common/fmatrix.hh>
+
+#include <dune/microstructure/matrix_operations.hh>
+#include <dune/microstructure/CorrectorComputer.hh>
+#include <dune/microstructure/EffectiveQuantitiesComputer.hh>
+#include <dune/microstructure/prestrainedMaterial.hh>
+
+#include <dune/solvers/solvers/umfpacksolver.hh>
+#include <dune/istl/eigenvalue/test/matrixinfo.hh>
+
+#include <dune/fufem/dunepython.hh>
+
+#include <any>
+#include <variant>
+#include <string>
+#include <iomanip>
+#include <filesystem>
+
+using namespace Dune;
+using namespace MatrixOperations;
+
+
+/**
+ * @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 + "/parametrizedlaminate.parset", parameterSet);
+ else
+ {
+ ParameterTreeParser::readINITree(argv[1], parameterSet);
+ ParameterTreeParser::readOptions(argc, argv, parameterSet);
+ }
+
+ // ParameterTree parameterSet;
+ // if (argc < 2)
+ // ParameterTreeParser::readINITree("parset/parametrizedlaminate.parset", parameterSet);
+ // else
+ // {
+ // ParameterTreeParser::readINITree(argv[1], parameterSet);
+ // ParameterTreeParser::readOptions(argc, argv, parameterSet);
+ // }
+
+ // ParameterTree parameterSet;
+ // if (argc < 2)
+ // ParameterTreeParser::readINITree("/parset/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::cout << "std::abs(q_1 - Qeff[0][0]):" << std::abs(q_1 - Qeff[0][0]) << std::endl;
+ 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;
+ return 0;
+}