diff --git a/src/model/solid_mechanics/materials/material_elastic_linear_anisotropic.cc b/src/model/solid_mechanics/materials/material_elastic_linear_anisotropic.cc index af25dff2e..42ee47945 100644 --- a/src/model/solid_mechanics/materials/material_elastic_linear_anisotropic.cc +++ b/src/model/solid_mechanics/materials/material_elastic_linear_anisotropic.cc @@ -1,263 +1,266 @@ /** * @file material_elastic_linear_anisotropic.cc * * @author Aurelia Isabel Cuba Ramos * @author Till Junge * @author Enrico Milanese * @author Nicolas Richart * * @date creation: Wed Sep 25 2013 * @date last modification: Tue Feb 20 2018 * * @brief Anisotropic elastic material * * @section LICENSE * * Copyright (©) 2014-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ #include "material_elastic_linear_anisotropic.hh" #include "solid_mechanics_model.hh" #include #include namespace akantu { /* -------------------------------------------------------------------------- */ template MaterialElasticLinearAnisotropic::MaterialElasticLinearAnisotropic( SolidMechanicsModel & model, const ID & id, bool symmetric) : Material(model, id), rot_mat(dim, dim), Cprime(dim * dim, dim * dim), C(voigt_h::size, voigt_h::size), eigC(voigt_h::size), symmetric(symmetric), alpha(0), was_stiffness_assembled(false) { AKANTU_DEBUG_IN(); this->dir_vecs.push_back(std::make_unique>(dim)); (*this->dir_vecs.back())[0] = 1.; this->registerParam("n1", *(this->dir_vecs.back()), _pat_parsmod, "Direction of main material axis"); if (dim > 1) { this->dir_vecs.push_back(std::make_unique>(dim)); (*this->dir_vecs.back())[1] = 1.; this->registerParam("n2", *(this->dir_vecs.back()), _pat_parsmod, "Direction of secondary material axis"); } if (dim > 2) { this->dir_vecs.push_back(std::make_unique>(dim)); (*this->dir_vecs.back())[2] = 1.; this->registerParam("n3", *(this->dir_vecs.back()), _pat_parsmod, "Direction of tertiary material axis"); } for (UInt i = 0; i < voigt_h::size; ++i) { UInt start = 0; if (this->symmetric) { start = i; } for (UInt j = start; j < voigt_h::size; ++j) { std::stringstream param("C"); param << "C" << i + 1 << j + 1; this->registerParam(param.str(), this->Cprime(i, j), Real(0.), _pat_parsmod, "Coefficient " + param.str()); } } this->registerParam("alpha", this->alpha, _pat_parsmod, "Proportion of viscous stress"); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialElasticLinearAnisotropic::initMaterial() { AKANTU_DEBUG_IN(); Material::initMaterial(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialElasticLinearAnisotropic::updateInternalParameters() { Material::updateInternalParameters(); if (this->symmetric) { for (UInt i = 0; i < voigt_h::size; ++i) { for (UInt j = i + 1; j < voigt_h::size; ++j) { this->Cprime(j, i) = this->Cprime(i, j); } } } this->rotateCprime(); this->C.eig(this->eigC); this->was_stiffness_assembled = false; } /* -------------------------------------------------------------------------- */ template void MaterialElasticLinearAnisotropic::rotateCprime() { // start by filling the empty parts fo Cprime UInt diff = Dim * Dim - voigt_h::size; for (UInt i = voigt_h::size; i < Dim * Dim; ++i) { for (UInt j = 0; j < Dim * Dim; ++j) { this->Cprime(i, j) = this->Cprime(i - diff, j); } } for (UInt i = 0; i < Dim * Dim; ++i) { for (UInt j = voigt_h::size; j < Dim * Dim; ++j) { this->Cprime(i, j) = this->Cprime(i, j - diff); } } // construction of rotator tensor // normalise rotation matrix for (UInt j = 0; j < Dim; ++j) { Vector rot_vec = this->rot_mat(j); rot_vec = *this->dir_vecs[j]; rot_vec.normalize(); } // make sure the vectors form a right-handed base Vector test_axis(3); - Vector v1(3), v2(3), v3(3); + Vector v1(3), v2(3), v3(3, 0.); if (Dim == 2) { for (UInt i = 0; i < Dim; ++i) { v1[i] = this->rot_mat(0, i); v2[i] = this->rot_mat(1, i); - v3[i] = 0.; } - v3[2] = 1.; - v1[2] = 0.; - v2[2] = 0.; + + v3.crossProduct(v1, v2); + if (v3.norm() < 8 * std::numeric_limits::epsilon()) { + AKANTU_ERROR("The axis vectors parallel."); + } + + v3.normalize(); } else if (Dim == 3) { v1 = this->rot_mat(0); v2 = this->rot_mat(1); v3 = this->rot_mat(2); } test_axis.crossProduct(v1, v2); test_axis -= v3; if (test_axis.norm() > 8 * std::numeric_limits::epsilon()) { AKANTU_ERROR("The axis vectors do not form a right-handed coordinate " << "system. I. e., ||n1 x n2 - n3|| should be zero, but " << "it is " << test_axis.norm() << "."); } // create the rotator and the reverse rotator Matrix rotator(Dim * Dim, Dim * Dim); Matrix revrotor(Dim * Dim, Dim * Dim); for (UInt i = 0; i < Dim; ++i) { for (UInt j = 0; j < Dim; ++j) { for (UInt k = 0; k < Dim; ++k) { for (UInt l = 0; l < Dim; ++l) { UInt I = voigt_h::mat[i][j]; UInt J = voigt_h::mat[k][l]; rotator(I, J) = this->rot_mat(k, i) * this->rot_mat(l, j); revrotor(I, J) = this->rot_mat(i, k) * this->rot_mat(j, l); } } } } // create the full rotated matrix Matrix Cfull(Dim * Dim, Dim * Dim); Cfull = rotator * Cprime * revrotor; for (UInt i = 0; i < voigt_h::size; ++i) { for (UInt j = 0; j < voigt_h::size; ++j) { this->C(i, j) = Cfull(i, j); } } } /* -------------------------------------------------------------------------- */ template void MaterialElasticLinearAnisotropic::computeStress( ElementType el_type, GhostType ghost_type) { // Wikipedia convention: // 2*eps_ij (i!=j) = voigt_eps_I // http://en.wikipedia.org/wiki/Voigt_notation AKANTU_DEBUG_IN(); MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, ghost_type); this->computeStressOnQuad(grad_u, sigma); MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END; } /* -------------------------------------------------------------------------- */ template void MaterialElasticLinearAnisotropic::computeTangentModuli( const ElementType & el_type, Array & tangent_matrix, GhostType ghost_type) { AKANTU_DEBUG_IN(); MATERIAL_TANGENT_QUADRATURE_POINT_LOOP_BEGIN(tangent_matrix); this->computeTangentModuliOnQuad(tangent); MATERIAL_TANGENT_QUADRATURE_POINT_LOOP_END; this->was_stiffness_assembled = true; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialElasticLinearAnisotropic::computePotentialEnergy( ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); Material::computePotentialEnergy(el_type, ghost_type); AKANTU_DEBUG_ASSERT(!this->finite_deformation, "finite deformation not possible in material anisotropic " "(TO BE IMPLEMENTED)"); if (ghost_type != _not_ghost) return; Array::scalar_iterator epot = this->potential_energy(el_type, ghost_type).begin(); MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, ghost_type); computePotentialEnergyOnQuad(grad_u, sigma, *epot); ++epot; MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template Real MaterialElasticLinearAnisotropic::getCelerity( __attribute__((unused)) const Element & element) const { return std::sqrt(this->eigC(0) / rho); } /* -------------------------------------------------------------------------- */ INSTANTIATE_MATERIAL(elastic_anisotropic, MaterialElasticLinearAnisotropic); } // akantu diff --git a/src/model/solid_mechanics/materials/material_elastic_orthotropic.cc b/src/model/solid_mechanics/materials/material_elastic_orthotropic.cc index d701bf394..26abc2e62 100644 --- a/src/model/solid_mechanics/materials/material_elastic_orthotropic.cc +++ b/src/model/solid_mechanics/materials/material_elastic_orthotropic.cc @@ -1,168 +1,168 @@ /** * @file material_elastic_orthotropic.cc * * @author Till Junge * @author Enrico Milanese * @author Nicolas Richart * @author Marco Vocialta * * @date creation: Fri Jun 18 2010 * @date last modification: Tue Feb 20 2018 * * @brief Orthotropic elastic material * * @section LICENSE * * Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ #include "material_elastic_orthotropic.hh" #include "solid_mechanics_model.hh" #include namespace akantu { /* -------------------------------------------------------------------------- */ template MaterialElasticOrthotropic::MaterialElasticOrthotropic( SolidMechanicsModel & model, const ID & id) : MaterialElasticLinearAnisotropic(model, id) { AKANTU_DEBUG_IN(); this->registerParam("E1", E1, Real(0.), _pat_parsmod, "Young's modulus (n1)"); this->registerParam("E2", E2, Real(0.), _pat_parsmod, "Young's modulus (n2)"); this->registerParam("nu12", nu12, Real(0.), _pat_parsmod, "Poisson's ratio (12)"); this->registerParam("G12", G12, Real(0.), _pat_parsmod, "Shear modulus (12)"); if (Dim > 2) { this->registerParam("E3", E3, Real(0.), _pat_parsmod, "Young's modulus (n3)"); this->registerParam("nu13", nu13, Real(0.), _pat_parsmod, "Poisson's ratio (13)"); this->registerParam("nu23", nu23, Real(0.), _pat_parsmod, "Poisson's ratio (23)"); this->registerParam("G13", G13, Real(0.), _pat_parsmod, "Shear modulus (13)"); this->registerParam("G23", G23, Real(0.), _pat_parsmod, "Shear modulus (23)"); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialElasticOrthotropic::initMaterial() { AKANTU_DEBUG_IN(); Material::initMaterial(); + AKANTU_DEBUG_ASSERT(not this->finite_deformation, + "finite deformation not possible in material orthotropic " + "(TO BE IMPLEMENTED)"); + updateInternalParameters(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialElasticOrthotropic::updateInternalParameters() { /* 1) construction of temporary material frame stiffness tensor------------ */ // http://solidmechanics.org/Text/Chapter3_2/Chapter3_2.php#Sect3_2_13 Real nu21 = nu12 * E2 / E1; Real nu31 = nu13 * E3 / E1; Real nu32 = nu23 * E3 / E2; // Full (i.e. dim^2 by dim^2) stiffness tensor in material frame if (Dim == 1) { AKANTU_ERROR("Dimensions 1 not implemented: makes no sense to have " "orthotropy for 1D"); } Real Gamma; if (Dim == 3) Gamma = 1 / (1 - nu12 * nu21 - nu23 * nu32 - nu31 * nu13 - 2 * nu21 * nu32 * nu13); if (Dim == 2) Gamma = 1 / (1 - nu12 * nu21); // Lamé's first parameters this->Cprime(0, 0) = E1 * (1 - nu23 * nu32) * Gamma; this->Cprime(1, 1) = E2 * (1 - nu13 * nu31) * Gamma; if (Dim == 3) this->Cprime(2, 2) = E3 * (1 - nu12 * nu21) * Gamma; // normalised poisson's ratio's this->Cprime(1, 0) = this->Cprime(0, 1) = E1 * (nu21 + nu31 * nu23) * Gamma; if (Dim == 3) { this->Cprime(2, 0) = this->Cprime(0, 2) = E1 * (nu31 + nu21 * nu32) * Gamma; this->Cprime(2, 1) = this->Cprime(1, 2) = E2 * (nu32 + nu12 * nu31) * Gamma; } // Lamé's second parameters (shear moduli) if (Dim == 3) { this->Cprime(3, 3) = G23; this->Cprime(4, 4) = G13; this->Cprime(5, 5) = G12; } else this->Cprime(2, 2) = G12; /* 1) rotation of C into the global frame */ this->rotateCprime(); this->C.eig(this->eigC); } /* -------------------------------------------------------------------------- */ template void MaterialElasticOrthotropic:: computePotentialEnergyByElement(ElementType type, UInt index, Vector & epot_on_quad_points) { - AKANTU_DEBUG_ASSERT(!this->finite_deformation, - "finite deformation not possible in material orthotropic " - "(TO BE IMPLEMENTED)"); - Array::matrix_iterator gradu_it = this->gradu(type).begin(spatial_dimension, spatial_dimension); Array::matrix_iterator gradu_end = this->gradu(type).begin(spatial_dimension, spatial_dimension); Array::matrix_iterator stress_it = this->stress(type).begin(spatial_dimension, spatial_dimension); UInt nb_quadrature_points = this->fem.getNbIntegrationPoints(type); gradu_it += index * nb_quadrature_points; gradu_end += (index + 1) * nb_quadrature_points; stress_it += index * nb_quadrature_points; Real * epot_quad = epot_on_quad_points.storage(); Matrix grad_u(spatial_dimension, spatial_dimension); for (; gradu_it != gradu_end; ++gradu_it, ++stress_it, ++epot_quad) { grad_u.copy(*gradu_it); this->computePotentialEnergyOnQuad(grad_u, *stress_it, *epot_quad); } } /* -------------------------------------------------------------------------- */ INSTANTIATE_MATERIAL(elastic_orthotropic, MaterialElasticOrthotropic); -} // akantu +} // namespace akantu diff --git a/test/test_model/patch_tests/CMakeLists.txt b/test/test_model/patch_tests/CMakeLists.txt index bcb5329d1..781b4c8da 100644 --- a/test/test_model/patch_tests/CMakeLists.txt +++ b/test/test_model/patch_tests/CMakeLists.txt @@ -1,113 +1,124 @@ #=============================================================================== # @file CMakeLists.txt # # @author Guillaume Anciaux # @author David Simon Kammer # @author Nicolas Richart # # @date creation: Fri Oct 22 2010 # @date last modification: Thu Feb 08 2018 # # @brief configuration for patch tests # # @section LICENSE # -# Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) +# Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) +# Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) # -# Akantu is free software: you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. +# Akantu is free software: you can redistribute it and/or modify it under the +# terms of the GNU Lesser General Public License as published by the Free +# Software Foundation, either version 3 of the License, or (at your option) any +# later version. # -# Akantu is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more details. +# Akantu is distributed in the hope that it will be useful, but WITHOUT ANY +# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR +# A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more +# details. # -# You should have received a copy of the GNU Lesser General Public License along with Akantu. If not, see . +# You should have received a copy of the GNU Lesser General Public License +# along with Akantu. If not, see . # # @section DESCRIPTION # #=============================================================================== add_subdirectory(data) register_gtest_sources(SOURCES patch_test_linear_elastic_explicit.cc PACKAGE solid_mechanics FILES_TO_COPY data/material_check_stress_plane_strain.dat data/material_check_stress_plane_stress.dat) register_gtest_sources(SOURCES patch_test_linear_elastic_implicit.cc PACKAGE solid_mechanics implicit FILES_TO_COPY data/material_check_stress_plane_strain.dat data/material_check_stress_plane_stress.dat) -register_gtest_sources(SOURCES patch_test_linear_anisotropic_explicit.cc +register_gtest_sources(SOURCES patch_test_linear_anisotropic.cc PACKAGE solid_mechanics lapack - UNSTABLE - FILES_TO_COPY data/material_anisotropic.dat) + FILES_TO_COPY + data/material_anisotropic_1.dat + data/material_anisotropic_2.dat + data/material_anisotropic_3.dat + ) register_gtest_sources(SOURCES test_lumped_mass.cc PACKAGE solid_mechanics FILES_TO_COPY data/material_lumped.dat) register_gtest_sources( SOURCES patch_test_linear_heat_transfer_explicit.cc FILES_TO_COPY data/heat_transfer_input.dat PACKAGE heat_transfer) register_gtest_sources( SOURCES patch_test_linear_heat_transfer_static.cc FILES_TO_COPY data/heat_transfer_input.dat PACKAGE heat_transfer implicit) register_gtest_sources( SOURCES patch_test_linear_heat_transfer_implicit.cc FILES_TO_COPY data/heat_transfer_input.dat PACKAGE heat_transfer implicit ) register_gtest_test(patch_test_linear FILES_TO_COPY ${PATCH_TEST_MESHES}) register_test(test_linear_elastic_explicit_python SCRIPT patch_test_linear_elastic_explicit.py PYTHON PACKAGE python_interface FILES_TO_COPY patch_test_linear_fixture.py FILES_TO_COPY patch_test_linear_solid_mechanics_fixture.py FILES_TO_COPY ${PATCH_TEST_MESHES}) register_test(test_linear_elastic_static_python SCRIPT patch_test_linear_elastic_static.py PYTHON PACKAGE python_interface implicit FILES_TO_COPY patch_test_linear_fixture.py FILES_TO_COPY patch_test_linear_solid_mechanics_fixture.py) register_test(test_linear_anisotropic_explicit_python SCRIPT patch_test_linear_anisotropic_explicit.py PYTHON UNSTABLE PACKAGE python_interface lapack FILES_TO_COPY patch_test_linear_fixture.py FILES_TO_COPY patch_test_linear_solid_mechanics_fixture.py FILES_TO_COPY data/material_anisotropic.dat) register_test(patch_test_linear_heat_transfer_explicit_python SCRIPT patch_test_linear_heat_transfer_explicit.py PYTHON PACKAGE python_interface heat_transfer FILES_TO_COPY patch_test_linear_fixture.py FILES_TO_COPY patch_test_linear_heat_transfer_fixture.py FILES_TO_COPY data/heat_transfer_input.dat) register_test(patch_test_linear_heat_transfer_static_python SCRIPT patch_test_linear_heat_transfer_static.py PYTHON PACKAGE python_interface heat_transfer implicit FILES_TO_COPY patch_test_linear_fixture.py FILES_TO_COPY patch_test_linear_heat_transfer_fixture.py FILES_TO_COPY data/heat_transfer_input.dat) register_test(patch_test_linear_heat_transfer_implicit_python SCRIPT patch_test_linear_heat_transfer_implicit.py PYTHON PACKAGE python_interface heat_transfer implicit FILES_TO_COPY patch_test_linear_fixture.py FILES_TO_COPY patch_test_linear_heat_transfer_fixture.py FILES_TO_COPY data/heat_transfer_input.dat) diff --git a/test/test_model/patch_tests/data/material_anisotropic_1.dat b/test/test_model/patch_tests/data/material_anisotropic_1.dat new file mode 100644 index 000000000..8d4afe6fc --- /dev/null +++ b/test/test_model/patch_tests/data/material_anisotropic_1.dat @@ -0,0 +1,7 @@ +material elastic_anisotropic [ + name = aluminum + rho = 1.6465129043044597 #gramm/mol/Å + C11 = 105.092023 #eV/Å^3 + + n1 = [-1] +] diff --git a/test/test_model/patch_tests/data/material_anisotropic_2.dat b/test/test_model/patch_tests/data/material_anisotropic_2.dat new file mode 100644 index 000000000..50d2460d1 --- /dev/null +++ b/test/test_model/patch_tests/data/material_anisotropic_2.dat @@ -0,0 +1,13 @@ +material elastic_anisotropic [ + name = aluminum + rho = 1.6465129043044597 #gramm/mol/Å + C11 = 105.092023 #eV/Å^3 + C12 = 59.4637759 + + C22 = 105.092023 + + C33 = 30.6596356 + + n1 = [-1, 1] + n2 = [ 1, 1] +] diff --git a/test/test_model/patch_tests/data/material_anisotropic.dat b/test/test_model/patch_tests/data/material_anisotropic_3.dat similarity index 100% rename from test/test_model/patch_tests/data/material_anisotropic.dat rename to test/test_model/patch_tests/data/material_anisotropic_3.dat diff --git a/test/test_model/patch_tests/patch_test_linear_anisotropic.cc b/test/test_model/patch_tests/patch_test_linear_anisotropic.cc new file mode 100644 index 000000000..6fd9511b6 --- /dev/null +++ b/test/test_model/patch_tests/patch_test_linear_anisotropic.cc @@ -0,0 +1,223 @@ +/** + * @file patch_test_linear_anisotropic_explicit.cc + * + * @author Guillaume Anciaux + * @author Till Junge + * @author David Simon Kammer + * @author Nicolas Richart + * @author Cyprien Wolff + * + * @date creation: Tue Dec 05 2017 + * @date last modification: Tue Feb 13 2018 + * + * @brief patch test for elastic material in solid mechanics model + * + * @section LICENSE + * + * Copyright (©) 2016-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) + * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) + * + * Akantu is free software: you can redistribute it and/or modify it under the + * terms of the GNU Lesser General Public License as published by the Free + * Software Foundation, either version 3 of the License, or (at your option) any + * later version. + * + * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY + * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR + * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more + * details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with Akantu. If not, see . + * + */ + +/* -------------------------------------------------------------------------- */ +#include "patch_test_linear_solid_mechanics_fixture.hh" +#include "non_linear_solver.hh" +/* -------------------------------------------------------------------------- */ + +using namespace akantu; + +// Stiffness tensor, rotated by hand + +/* -------------------------------------------------------------------------- */ +TYPED_TEST(TestPatchTestSMMLinear, AnisotropicExplicit) { + Real C[3][3][3][3] = { + {{{112.93753505, 1.85842452538e-10, -4.47654358027e-10}, + {1.85847317471e-10, 54.2334345331, -3.69840984824}, + {-4.4764768395e-10, -3.69840984824, 56.848605217}}, + {{1.85847781609e-10, 25.429294233, -3.69840984816}, + {25.429294233, 3.31613847493e-10, -8.38797920011e-11}, + {-3.69840984816, -8.38804581349e-11, -1.97875715813e-10}}, + {{-4.47654358027e-10, -3.69840984816, 28.044464917}, + {-3.69840984816, 2.09374961813e-10, 9.4857455224e-12}, + {28.044464917, 9.48308098714e-12, -2.1367885239e-10}}}, + {{{1.85847781609e-10, 25.429294233, -3.69840984816}, + {25.429294233, 3.31613847493e-10, -8.38793479119e-11}, + {-3.69840984816, -8.38795699565e-11, -1.97876381947e-10}}, + {{54.2334345331, 3.31617400207e-10, 2.09372075233e-10}, + {3.3161562385e-10, 115.552705733, -3.15093728886e-10}, + {2.09372075233e-10, -3.15090176173e-10, 54.2334345333}}, + {{-3.69840984824, -8.38795699565e-11, 9.48219280872e-12}, + {-8.38795699565e-11, -3.1509195253e-10, 25.4292942335}, + {9.48441325477e-12, 25.4292942335, 3.69840984851}}}, + {{{-4.47653469848e-10, -3.69840984816, 28.044464917}, + {-3.69840984816, 2.09374073634e-10, 9.48752187924e-12}, + {28.044464917, 9.48552347779e-12, -2.1367885239e-10}}, + {{-3.69840984824, -8.3884899027e-11, 9.48219280872e-12}, + {-8.3884899027e-11, -3.150972816e-10, 25.4292942335}, + {9.48041645188e-12, 25.4292942335, 3.69840984851}}, + {{56.848605217, -1.97875493768e-10, -2.13681516925e-10}, + {-1.97877270125e-10, 54.2334345333, 3.69840984851}, + {-2.13683293282e-10, 3.69840984851, 112.93753505}}}}; + + if (this->dim == 2) { + for (UInt i = 0; i < this->dim; ++i) { + for (UInt j = 0; j < this->dim; ++j) { + for (UInt k = 0; k < this->dim; ++k) { + for (UInt l = 0; l < this->dim; ++l) { + C[i][j][k][l] = 0; + } + } + } + } + C[0][0][0][0] = C[1][1][1][1] = 112.93753504999995; + C[0][0][1][1] = C[1][1][0][0] = 51.618263849999984; + C[0][1][0][1] = C[1][0][0][1] = C[0][1][1][0] = C[1][0][1][0] = 22.814123549999987; + } + + if (this->dim == 1) { + C[0][0][0][0] = 105.092023; + } + this->initModel(_explicit_lumped_mass, + "material_anisotropic_" + std::to_string(this->dim) + ".dat"); + + const auto & coordinates = this->mesh->getNodes(); + auto & displacement = this->model->getDisplacement(); + + // set the position of all nodes to the static solution + for (auto && tuple : zip(make_view(coordinates, this->dim), + make_view(displacement, this->dim))) { + this->setLinearDOF(std::get<1>(tuple), std::get<0>(tuple)); + } + + for (UInt s = 0; s < 100; ++s) { + this->model->solveStep(); + } + + auto ekin = this->model->getEnergy("kinetic"); + EXPECT_NEAR(0, ekin, 1e-16); + + auto & mat = this->model->getMaterial(0); + + this->checkDOFs(displacement); + this->checkGradient(mat.getGradU(this->type), displacement); + + this->result_tolerance = 1e-11; + this->checkResults( + [&](const Matrix & pstrain) { + auto strain = (pstrain + pstrain.transpose()) / 2.; + decltype(strain) stress(this->dim, this->dim); + + for (UInt i = 0; i < this->dim; ++i) { + for (UInt j = 0; j < this->dim; ++j) { + stress(i, j) = 0; + for (UInt k = 0; k < this->dim; ++k) { + for (UInt l = 0; l < this->dim; ++l) { + stress(i, j) += C[i][j][k][l] * strain(k, l); + } + } + } + } + return stress; + }, + mat.getStress(this->type), displacement); +} + + +TYPED_TEST(TestPatchTestSMMLinear, AnisotropicStatic) { + Real C[3][3][3][3] = { + {{{112.93753505, 1.85842452538e-10, -4.47654358027e-10}, + {1.85847317471e-10, 54.2334345331, -3.69840984824}, + {-4.4764768395e-10, -3.69840984824, 56.848605217}}, + {{1.85847781609e-10, 25.429294233, -3.69840984816}, + {25.429294233, 3.31613847493e-10, -8.38797920011e-11}, + {-3.69840984816, -8.38804581349e-11, -1.97875715813e-10}}, + {{-4.47654358027e-10, -3.69840984816, 28.044464917}, + {-3.69840984816, 2.09374961813e-10, 9.4857455224e-12}, + {28.044464917, 9.48308098714e-12, -2.1367885239e-10}}}, + {{{1.85847781609e-10, 25.429294233, -3.69840984816}, + {25.429294233, 3.31613847493e-10, -8.38793479119e-11}, + {-3.69840984816, -8.38795699565e-11, -1.97876381947e-10}}, + {{54.2334345331, 3.31617400207e-10, 2.09372075233e-10}, + {3.3161562385e-10, 115.552705733, -3.15093728886e-10}, + {2.09372075233e-10, -3.15090176173e-10, 54.2334345333}}, + {{-3.69840984824, -8.38795699565e-11, 9.48219280872e-12}, + {-8.38795699565e-11, -3.1509195253e-10, 25.4292942335}, + {9.48441325477e-12, 25.4292942335, 3.69840984851}}}, + {{{-4.47653469848e-10, -3.69840984816, 28.044464917}, + {-3.69840984816, 2.09374073634e-10, 9.48752187924e-12}, + {28.044464917, 9.48552347779e-12, -2.1367885239e-10}}, + {{-3.69840984824, -8.3884899027e-11, 9.48219280872e-12}, + {-8.3884899027e-11, -3.150972816e-10, 25.4292942335}, + {9.48041645188e-12, 25.4292942335, 3.69840984851}}, + {{56.848605217, -1.97875493768e-10, -2.13681516925e-10}, + {-1.97877270125e-10, 54.2334345333, 3.69840984851}, + {-2.13683293282e-10, 3.69840984851, 112.93753505}}}}; + + if (this->dim == 2) { + for (UInt i = 0; i < this->dim; ++i) { + for (UInt j = 0; j < this->dim; ++j) { + for (UInt k = 0; k < this->dim; ++k) { + for (UInt l = 0; l < this->dim; ++l) { + C[i][j][k][l] = 0; + } + } + } + } + C[0][0][0][0] = C[1][1][1][1] = 112.93753504999995; + C[0][0][1][1] = C[1][1][0][0] = 51.618263849999984; + C[0][1][0][1] = C[1][0][0][1] = C[0][1][1][0] = C[1][0][1][0] = 22.814123549999987; + } + + if (this->dim == 1) { + C[0][0][0][0] = 105.092023; + } + + this->initModel(_static, + "material_anisotropic_" + std::to_string(this->dim) + ".dat"); + + auto & solver = this->model->getNonLinearSolver(); + solver.set("max_iterations", 2); + solver.set("threshold", 2e-4); + solver.set("convergence_type", _scc_residual); + + this->model->solveStep(); + + auto & mat = this->model->getMaterial(0); + + const auto & displacement = this->model->getDisplacement(); + this->checkDOFs(displacement); + this->checkGradient(mat.getGradU(this->type), displacement); + + this->result_tolerance = 1e-11; + this->checkResults( + [&](const Matrix & pstrain) { + auto strain = (pstrain + pstrain.transpose()) / 2.; + decltype(strain) stress(this->dim, this->dim); + + for (UInt i = 0; i < this->dim; ++i) { + for (UInt j = 0; j < this->dim; ++j) { + stress(i, j) = 0; + for (UInt k = 0; k < this->dim; ++k) { + for (UInt l = 0; l < this->dim; ++l) { + stress(i, j) += C[i][j][k][l] * strain(k, l); + } + } + } + } + return stress; + }, + mat.getStress(this->type), displacement); +} diff --git a/test/test_model/patch_tests/patch_test_linear_anisotropic_explicit.cc b/test/test_model/patch_tests/patch_test_linear_anisotropic_explicit.cc deleted file mode 100644 index d9b29efee..000000000 --- a/test/test_model/patch_tests/patch_test_linear_anisotropic_explicit.cc +++ /dev/null @@ -1,114 +0,0 @@ -/** - * @file patch_test_linear_anisotropic_explicit.cc - * - * @author Guillaume Anciaux - * @author Till Junge - * @author David Simon Kammer - * @author Nicolas Richart - * @author Cyprien Wolff - * - * @date creation: Tue Dec 05 2017 - * @date last modification: Tue Feb 13 2018 - * - * @brief patch test for elastic material in solid mechanics model - * - * @section LICENSE - * - * Copyright (©) 2016-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) - * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) - * - * Akantu is free software: you can redistribute it and/or modify it under the - * terms of the GNU Lesser General Public License as published by the Free - * Software Foundation, either version 3 of the License, or (at your option) any - * later version. - * - * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY - * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR - * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more - * details. - * - * You should have received a copy of the GNU Lesser General Public License - * along with Akantu. If not, see . - * - */ - -/* -------------------------------------------------------------------------- */ -#include "patch_test_linear_solid_mechanics_fixture.hh" -/* -------------------------------------------------------------------------- */ - -using namespace akantu; - -// Stiffness tensor, rotated by hand -Real C[3][3][3][3] = { - {{{112.93753505, 1.85842452538e-10, -4.47654358027e-10}, - {1.85847317471e-10, 54.2334345331, -3.69840984824}, - {-4.4764768395e-10, -3.69840984824, 56.848605217}}, - {{1.85847781609e-10, 25.429294233, -3.69840984816}, - {25.429294233, 3.31613847493e-10, -8.38797920011e-11}, - {-3.69840984816, -8.38804581349e-11, -1.97875715813e-10}}, - {{-4.47654358027e-10, -3.69840984816, 28.044464917}, - {-3.69840984816, 2.09374961813e-10, 9.4857455224e-12}, - {28.044464917, 9.48308098714e-12, -2.1367885239e-10}}}, - {{{1.85847781609e-10, 25.429294233, -3.69840984816}, - {25.429294233, 3.31613847493e-10, -8.38793479119e-11}, - {-3.69840984816, -8.38795699565e-11, -1.97876381947e-10}}, - {{54.2334345331, 3.31617400207e-10, 2.09372075233e-10}, - {3.3161562385e-10, 115.552705733, -3.15093728886e-10}, - {2.09372075233e-10, -3.15090176173e-10, 54.2334345333}}, - {{-3.69840984824, -8.38795699565e-11, 9.48219280872e-12}, - {-8.38795699565e-11, -3.1509195253e-10, 25.4292942335}, - {9.48441325477e-12, 25.4292942335, 3.69840984851}}}, - {{{-4.47653469848e-10, -3.69840984816, 28.044464917}, - {-3.69840984816, 2.09374073634e-10, 9.48752187924e-12}, - {28.044464917, 9.48552347779e-12, -2.1367885239e-10}}, - {{-3.69840984824, -8.3884899027e-11, 9.48219280872e-12}, - {-8.3884899027e-11, -3.150972816e-10, 25.4292942335}, - {9.48041645188e-12, 25.4292942335, 3.69840984851}}, - {{56.848605217, -1.97875493768e-10, -2.13681516925e-10}, - {-1.97877270125e-10, 54.2334345333, 3.69840984851}, - {-2.13683293282e-10, 3.69840984851, 112.93753505}}}}; - -/* -------------------------------------------------------------------------- */ -TYPED_TEST(TestPatchTestSMMLinear, AnisotropicExplicit) { - this->initModel(_explicit_lumped_mass, "material_anisotropic.dat"); - - const auto & coordinates = this->mesh->getNodes(); - auto & displacement = this->model->getDisplacement(); - - // set the position of all nodes to the static solution - for (auto && tuple : zip(make_view(coordinates, this->dim), - make_view(displacement, this->dim))) { - this->setLinearDOF(std::get<1>(tuple), std::get<0>(tuple)); - } - - for (UInt s = 0; s < 100; ++s) { - this->model->solveStep(); - } - - auto ekin = this->model->getEnergy("kinetic"); - EXPECT_NEAR(0, ekin, 1e-16); - - auto & mat = this->model->getMaterial(0); - - this->checkDOFs(displacement); - this->checkGradient(mat.getGradU(this->type), displacement); - - this->checkResults( - [&](const Matrix & pstrain) { - auto strain = (pstrain + pstrain.transpose()) / 2.; - decltype(strain) stress(this->dim, this->dim); - - for (UInt i = 0; i < this->dim; ++i) { - for (UInt j = 0; j < this->dim; ++j) { - stress(i, j) = 0; - for (UInt k = 0; k < this->dim; ++k) { - for (UInt l = 0; l < this->dim; ++l) { - stress(i, j) += C[i][j][k][l] * strain(k, l); - } - } - } - } - return stress; - }, - mat.getStress(this->type), displacement); -} diff --git a/test/test_model/patch_tests/patch_test_linear_fixture.hh b/test/test_model/patch_tests/patch_test_linear_fixture.hh index 096887109..e866bf5e7 100644 --- a/test/test_model/patch_tests/patch_test_linear_fixture.hh +++ b/test/test_model/patch_tests/patch_test_linear_fixture.hh @@ -1,183 +1,183 @@ /** * @file patch_test_linear_fixture.hh * * @author Nicolas Richart * * @date creation: Tue Jan 30 2018 * @date last modification: Wed Jan 31 2018 * * @brief Fixture for linear patch tests * * @section LICENSE * * Copyright (©) 2016-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "element_group.hh" #include "mesh_utils.hh" #include "model.hh" #include "test_gtest_utils.hh" /* -------------------------------------------------------------------------- */ #include #include /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_PATCH_TEST_LINEAR_FIXTURE_HH__ #define __AKANTU_PATCH_TEST_LINEAR_FIXTURE_HH__ //#define DEBUG_TEST using namespace akantu; template class TestPatchTestLinear : public ::testing::Test { public: static constexpr ElementType type = type_::value; static constexpr size_t dim = ElementClass::getSpatialDimension(); virtual void SetUp() { mesh = std::make_unique(dim); mesh->read(aka::to_string(type) + ".msh"); MeshUtils::buildFacets(*mesh); mesh->createBoundaryGroupFromGeometry(); model = std::make_unique(*mesh); } virtual void TearDown() { model.reset(nullptr); mesh.reset(nullptr); } virtual void initModel(const AnalysisMethod & method, const std::string & material_file) { debug::setDebugLevel(dblError); getStaticParser().parse(material_file); this->model->initFull(_analysis_method = method); this->applyBC(); if (method != _static) this->model->setTimeStep(0.8 * this->model->getStableTimeStep()); } virtual void applyBC() { auto & boundary = this->model->getBlockedDOFs(); for (auto & eg : mesh->getElementGroups()) { for (const auto & node : eg.second->getNodeGroup()) { for (UInt s = 0; s < boundary.getNbComponent(); ++s) { boundary(node, s) = true; } } } } virtual void applyBConDOFs(const Array & dofs) { const auto & coordinates = this->mesh->getNodes(); for (auto & eg : this->mesh->getElementGroups()) { for (const auto & node : eg.second->getNodeGroup()) { this->setLinearDOF(dofs.begin(dofs.getNbComponent())[node], coordinates.begin(this->dim)[node]); } } } template Matrix prescribed_gradient(const V & dof) { Matrix gradient(dof.getNbComponent(), dim); for (UInt i = 0; i < gradient.rows(); ++i) { for (UInt j = 0; j < gradient.cols(); ++j) { gradient(i, j) = alpha(i, j + 1); } } return gradient; } template void checkGradient(const Gradient & gradient, const DOFs & dofs) { auto pgrad = prescribed_gradient(dofs); for (auto & grad : make_view(gradient, gradient.getNbComponent() / dim, dim)) { auto diff = grad - pgrad; auto gradient_error = diff.template norm() / grad.template norm(); EXPECT_NEAR(0, gradient_error, gradient_tolerance); } } - template + template void checkResults(presult_func_t && presult_func, const Result & results, const DOFs & dofs) { auto presult = presult_func(prescribed_gradient(dofs)); for (auto & result : make_view(results, results.getNbComponent() / dim, dim)) { auto diff = result - presult; auto result_error = - diff.template norm() / result.template norm(); + diff.template norm() / presult.template norm(); EXPECT_NEAR(0, result_error, result_tolerance); } } template void setLinearDOF(V1 && dof, V2 && coord) { for (UInt i = 0; i < dof.size(); ++i) { dof(i) = this->alpha(i, 0); for (UInt j = 0; j < coord.size(); ++j) { dof(i) += this->alpha(i, j + 1) * coord(j); } } } template void checkDOFs(V && dofs) { const auto & coordinates = mesh->getNodes(); Vector ref_dof(dofs.getNbComponent()); for (auto && tuple : zip(make_view(coordinates, dim), make_view(dofs, dofs.getNbComponent()))) { setLinearDOF(ref_dof, std::get<0>(tuple)); auto diff = std::get<1>(tuple) - ref_dof; auto dofs_error = diff.template norm(); EXPECT_NEAR(0, dofs_error, dofs_tolerance); } } protected: std::unique_ptr mesh; std::unique_ptr model; Matrix alpha{{0.01, 0.02, 0.03, 0.04}, {0.05, 0.06, 0.07, 0.08}, {0.09, 0.10, 0.11, 0.12}}; Real gradient_tolerance{1e-13}; Real result_tolerance{1e-13}; Real dofs_tolerance{1e-15}; }; template constexpr ElementType TestPatchTestLinear::type; template constexpr size_t TestPatchTestLinear::dim; #endif /* __AKANTU_PATCH_TEST_LINEAR_FIXTURE_HH__ */