From 4b989f17458c686a5104a80beef7111b8e172d44 Mon Sep 17 00:00:00 2001
From: Klaus <klaus.boehnlein@tu-dresden.de>
Date: Wed, 4 Oct 2023 16:53:43 +0200
Subject: [PATCH] Add Test against parametrized_Laminate (analytical results)
---
test/parametrized_laminate.py | 77 +++++++++++
test/parametrizedlaminate.parset | 11 ++
test/parametrizedlaminatetest.cc | 221 ++++++++++++++++++++++++++++---
3 files changed, 287 insertions(+), 22 deletions(-)
create mode 100644 test/parametrized_laminate.py
create mode 100644 test/parametrizedlaminate.parset
diff --git a/test/parametrized_laminate.py b/test/parametrized_laminate.py
new file mode 100644
index 00000000..9007b436
--- /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 00000000..08e2ee45
--- /dev/null
+++ b/test/parametrizedlaminate.parset
@@ -0,0 +1,11 @@
+# --- 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
diff --git a/test/parametrizedlaminatetest.cc b/test/parametrizedlaminatetest.cc
index 12ef53d8..544e8aaa 100644
--- a/test/parametrizedlaminatetest.cc
+++ b/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;
+}
--
GitLab