diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index f3c48b611d52fd589b30fc0f24ed8665064460d6..f2e918562b1e2ea68456621c709640e767a8ec12 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -1,9 +1,18 @@ set(TESTS - parametrizedlaminatetest) + 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/parametrizedlaminate.parset b/test/parametrizedlaminate.parset index 08e2ee453f5f0124a8e17cb47bd57754bf22e045..f02b8e746b934135082bba9aa5a3c913269570f9 100644 --- a/test/parametrizedlaminate.parset +++ b/test/parametrizedlaminate.parset @@ -7,5 +7,14 @@ 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 + +# 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 index 544e8aaae1177cfd92d24202d001e8e826f8f1c3..636df2beee9289422265b92d3a4bf00e12b059b8 100644 --- a/test/parametrizedlaminatetest.cc +++ b/test/parametrizedlaminatetest.cc @@ -98,13 +98,31 @@ int main(int argc, char *argv[]) ParameterTree parameterSet; if (argc < 2) - ParameterTreeParser::readINITree(dir_path + "/parsets/parametrizedlaminate.parset", parameterSet); + 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__"); @@ -175,8 +193,8 @@ int main(int argc, char *argv[]) //--- get effective quantities auto Qeff = effectiveQuantitiesComputer.getQeff(); auto Beff = effectiveQuantitiesComputer.getBeff(); -// printmatrix(std::cout, Qeff, "Matrix Qeff", "--"); -// printvector(std::cout, Beff, "Beff", "--"); + printmatrix(std::cout, Qeff, "Matrix Qeff", "--"); + printvector(std::cout, Beff, "Beff", "--"); /** * @brief Compute analytical effective quantitites @@ -209,6 +227,7 @@ int main(int argc, char *argv[]) */ 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; @@ -251,4 +270,5 @@ int main(int argc, char *argv[]) std::cout << "Total time elapsed: " << globalTimer.elapsed() << std::endl; std::cout << "Test parametrizedlaminatetest passed." << std::endl; + return 0; }