diff --git a/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_fixture.hh b/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_fixture.hh index 5766243ff..30c63b5f5 100644 --- a/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_fixture.hh +++ b/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_fixture.hh @@ -1,330 +1,333 @@ /** * @file test_cohesive_fixture.hh * * @author Nicolas Richart * * @date creation: Wed Jan 10 2018 * @date last modification: Tue Feb 20 2018 * * @brief Coehsive element test fixture * * @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 "communicator.hh" #include "solid_mechanics_model_cohesive.hh" #include "test_gtest_utils.hh" /* -------------------------------------------------------------------------- */ #include #include /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_TEST_COHESIVE_FIXTURE_HH__ #define __AKANTU_TEST_COHESIVE_FIXTURE_HH__ using namespace akantu; template <::akantu::AnalysisMethod t> using analysis_method_t = std::integral_constant<::akantu::AnalysisMethod, t>; class StrainIncrement : public BC::Functor { public: StrainIncrement(const Matrix & strain, BC::Axis dir) : strain_inc(strain), dir(dir) {} void operator()(UInt /*node*/, Vector & flags, Vector & primal, const Vector & coord) const { if (std::abs(coord(dir)) < 1e-8) { return; } flags.set(true); primal += strain_inc * coord; } static const BC::Functor::Type type = BC::Functor::_dirichlet; private: Matrix strain_inc; BC::Axis dir; }; template class TestSMMCFixture : public ::testing::Test { public: static constexpr ElementType cohesive_type = std::tuple_element_t<0, param_>::value; static constexpr ElementType type_1 = std::tuple_element_t<1, param_>::value; static constexpr ElementType type_2 = std::tuple_element_t<2, param_>::value; static constexpr size_t dim = ElementClass::getSpatialDimension(); void SetUp() override { mesh = std::make_unique(this->dim); if (Communicator::getStaticCommunicator().whoAmI() == 0) { EXPECT_NO_THROW({ mesh->read(this->mesh_name); }); } mesh->distribute(); } void TearDown() override { model.reset(nullptr); mesh.reset(nullptr); } void createModel() { model = std::make_unique(*mesh); model->initFull(_analysis_method = this->analysis_method, _is_extrinsic = this->is_extrinsic); auto time_step = this->model->getStableTimeStep() * 0.01; this->model->setTimeStep(time_step); auto & mat_el = this->model->getMaterial("body"); - auto speed = mat_el.getShearWaveSpeed(Element()); // the slowest speed + auto speed = mat_el.getPushWaveSpeed(Element()); + try { + speed = mat_el.getShearWaveSpeed(Element()); // the slowest speed if exists + } catch(...) {} auto directon = _y; if(dim == 1) directon = _x; auto length = mesh->getUpperBounds()(directon) - mesh->getLowerBounds()(directon); nb_steps = length / 2. / speed / time_step; if (dim == 1) { surface = 1; - group_size = 1; - return; + group_size = 1; + return; } auto facet_type = mesh->getFacetType(this->cohesive_type); auto & fe_engine = model->getFEEngineBoundary(); auto & group = mesh->getElementGroup("insertion").getElements(facet_type); Array ones(fe_engine.getNbIntegrationPoints(facet_type) * group.size()); ones.set(1.); surface = fe_engine.integrate(ones, facet_type, _not_ghost, group); mesh->getCommunicator().allReduce(surface, SynchronizerOperation::_sum); group_size = group.size(); mesh->getCommunicator().allReduce(group_size, SynchronizerOperation::_sum); #define debug_ 0 #if debug_ this->model->addDumpFieldVector("displacement"); this->model->addDumpFieldVector("velocity"); this->model->addDumpField("stress"); this->model->addDumpField("strain"); this->model->assembleInternalForces(); this->model->setBaseNameToDumper("cohesive elements", "cohesive_elements"); this->model->addDumpFieldVectorToDumper("cohesive elements", "displacement"); this->model->addDumpFieldToDumper("cohesive elements", "damage"); this->model->addDumpFieldToDumper("cohesive elements", "tractions"); this->model->addDumpFieldToDumper("cohesive elements", "opening"); this->model->dump(); this->model->dump("cohesive elements"); #endif } void setInitialCondition(const Matrix & strain) { for (auto && data : zip(make_view(this->mesh->getNodes(), this->dim), make_view(this->model->getDisplacement(), this->dim))) { const auto & pos = std::get<0>(data); auto & disp = std::get<1>(data); disp = strain * pos; } } bool checkDamaged() { UInt nb_damaged = 0; auto & damage = model->getMaterial("insertion").getArray("damage", cohesive_type); for (auto d : damage) { if (d >= .99) ++nb_damaged; } return (nb_damaged == group_size); } void steps(const Matrix & strain) { StrainIncrement functor((1. / 300) * strain, this->dim == 1 ? _x : _y); for (auto _[[gnu::unused]] : arange(nb_steps)) { this->model->applyBC(functor, "loading"); this->model->applyBC(functor, "fixed"); if (this->is_extrinsic) this->model->checkCohesiveStress(); this->model->solveStep(); #if debug_ this->model->dump(); this->model->dump("cohesive elements"); #endif } } void checkInsertion() { auto nb_cohesive_element = this->mesh->getNbElement(cohesive_type); mesh->getCommunicator().allReduce(nb_cohesive_element, SynchronizerOperation::_sum); EXPECT_EQ(nb_cohesive_element, group_size); } void checkDissipated(Real expected_density) { Real edis = this->model->getEnergy("dissipated"); EXPECT_NEAR(this->surface * expected_density, edis, 5e-1); } void testModeI() { // EXPECT_NO_THROW(this->createModel()); this->createModel(); SCOPED_TRACE(std::to_string(this->dim) + "D - " + aka::to_string(type_1) + ":" + aka::to_string(type_2)); auto & mat_co = this->model->getMaterial("insertion"); Real sigma_c = mat_co.get("sigma_c"); auto & mat_el = this->model->getMaterial("body"); Real E = mat_el.get("E"); Real nu = mat_el.get("nu"); Matrix strain; if (dim == 1) { strain = {{1.}}; } else if (dim == 2) { strain = {{-nu, 0.}, {0., 1. - nu}}; strain *= (1. + nu); } else if (dim == 3) { strain = {{-nu, 0., 0.}, {0., 1., 0.}, {0., 0., -nu}}; } strain *= sigma_c / E; this->setInitialCondition((1 - 1e-5) * strain); this->steps(1e-2 * strain); } void testModeII() { Vector direction(this->dim, 0.); direction(_x) = 1.; SCOPED_TRACE(std::to_string(this->dim) + "D - " + aka::to_string(type_1) + ":" + aka::to_string(type_2)); EXPECT_NO_THROW(this->createModel()); if (this->dim > 1) this->model->applyBC(BC::Dirichlet::FlagOnly(_y), "sides"); if (this->dim > 2) this->model->applyBC(BC::Dirichlet::FlagOnly(_z), "sides"); auto & mat_co = this->model->getMaterial("insertion"); Real sigma_c = mat_co.get("sigma_c"); Real beta = mat_co.get("beta"); // Real G_c = mat_co.get("G_c"); auto & mat_el = this->model->getMaterial("body"); Real E = mat_el.get("E"); Real nu = mat_el.get("nu"); Matrix strain; if (dim == 1) { strain = {{1.}}; } else if (dim == 2) { strain = {{0., 1.}, {0., 0.}}; strain *= (1. + nu); } else if (dim == 3) { strain = {{0., 1., 0.}, {0., 0., 0.}, {0., 0., 0.}}; strain *= (1. + nu); } strain *= 2 * beta * beta * sigma_c / E; //nb_steps *= 5; this->setInitialCondition((1. - 1e-5) * strain); this->steps(0.005 * strain); } protected: std::unique_ptr mesh; std::unique_ptr model; std::string mesh_name{aka::to_string(cohesive_type) + aka::to_string(type_1) + (type_1 == type_2 ? "" : aka::to_string(type_2)) + ".msh"}; bool is_extrinsic; AnalysisMethod analysis_method; Real surface{0}; UInt nb_steps{1000}; UInt group_size{10000}; }; /* -------------------------------------------------------------------------- */ template constexpr ElementType TestSMMCFixture::cohesive_type; template constexpr ElementType TestSMMCFixture::type_1; template constexpr ElementType TestSMMCFixture::type_2; template constexpr size_t TestSMMCFixture::dim; /* -------------------------------------------------------------------------- */ using IsExtrinsicTypes = std::tuple; using AnalysisMethodTypes = std::tuple>; using coh_types = gtest_list_t, element_type_t<_segment_2>, element_type_t<_segment_2>>, std::tuple, element_type_t<_triangle_3>, element_type_t<_triangle_3>>, std::tuple, element_type_t<_quadrangle_4>, element_type_t<_quadrangle_4>>, std::tuple, element_type_t<_triangle_3>, element_type_t<_quadrangle_4>>, std::tuple, element_type_t<_triangle_6>, element_type_t<_triangle_6>>, std::tuple, element_type_t<_quadrangle_8>, element_type_t<_quadrangle_8>>, std::tuple, element_type_t<_triangle_6>, element_type_t<_quadrangle_8>>, std::tuple, element_type_t<_tetrahedron_4>, element_type_t<_tetrahedron_4>>, std::tuple, element_type_t<_tetrahedron_10>, element_type_t<_tetrahedron_10>> /*, std::tuple, element_type_t<_hexahedron_8>, element_type_t<_hexahedron_8>>, std::tuple, element_type_t<_hexahedron_20>, element_type_t<_hexahedron_20>>*/>>; TYPED_TEST_CASE(TestSMMCFixture, coh_types); #endif /* __AKANTU_TEST_COHESIVE_FIXTURE_HH__ */