diff --git a/src/mesh_utils/cohesive_element_inserter.cc b/src/mesh_utils/cohesive_element_inserter.cc index 793ec51c2..e001ed0d4 100644 --- a/src/mesh_utils/cohesive_element_inserter.cc +++ b/src/mesh_utils/cohesive_element_inserter.cc @@ -1,274 +1,274 @@ /** * @file cohesive_element_inserter.cc * * @author Marco Vocialta * * @date creation: Wed Dec 04 2013 * @date last modification: Mon Feb 19 2018 * * @brief Cohesive element inserter functions * * @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 "cohesive_element_inserter.hh" #include "communicator.hh" #include "element_group.hh" #include "global_ids_updater.hh" #include "mesh_accessor.hh" #include "mesh_iterators.hh" #include "element_synchronizer.hh" /* -------------------------------------------------------------------------- */ #include #include /* -------------------------------------------------------------------------- */ namespace akantu { CohesiveElementInserter::CohesiveElementInserter(Mesh & mesh, const ID & id) : Parsable(ParserType::_cohesive_inserter), id(id), mesh(mesh), mesh_facets(mesh.initMeshFacets()), insertion_facets("insertion_facets", id), insertion_limits(mesh.getSpatialDimension(), 2), check_facets("check_facets", id) { this->registerParam("cohesive_surfaces", physical_groups, _pat_parsable, "List of groups to consider for insertion"); this->registerParam("bounding_box", insertion_limits, _pat_parsable, "Global limit for insertion"); UInt spatial_dimension = mesh.getSpatialDimension(); MeshUtils::resetFacetToDouble(mesh_facets); /// init insertion limits for (UInt dim = 0; dim < spatial_dimension; ++dim) { insertion_limits(dim, 0) = std::numeric_limits::max() * Real(-1.); insertion_limits(dim, 1) = std::numeric_limits::max(); } insertion_facets.initialize(mesh_facets, _spatial_dimension = spatial_dimension - 1, _with_nb_element = true, _default_value = false); } /* -------------------------------------------------------------------------- */ CohesiveElementInserter::~CohesiveElementInserter() = default; /* -------------------------------------------------------------------------- */ void CohesiveElementInserter::parseSection(const ParserSection & section) { Parsable::parseSection(section); if (is_extrinsic) limitCheckFacets(this->check_facets); } /* -------------------------------------------------------------------------- */ void CohesiveElementInserter::limitCheckFacets() { limitCheckFacets(this->check_facets); } /* -------------------------------------------------------------------------- */ -void CohesiveElementInserter::setLimit(SpacialDirection axis, Real first_limit, +void CohesiveElementInserter::setLimit(SpatialDirection axis, Real first_limit, Real second_limit) { AKANTU_DEBUG_ASSERT( axis < mesh.getSpatialDimension(), "You are trying to limit insertion in a direction that doesn't exist"); insertion_limits(axis, 0) = std::min(first_limit, second_limit); insertion_limits(axis, 1) = std::max(first_limit, second_limit); } /* -------------------------------------------------------------------------- */ UInt CohesiveElementInserter::insertIntrinsicElements() { limitCheckFacets(insertion_facets); return insertElements(); } /* -------------------------------------------------------------------------- */ void CohesiveElementInserter::limitCheckFacets( ElementTypeMapArray & check_facets) { AKANTU_DEBUG_IN(); UInt spatial_dimension = mesh.getSpatialDimension(); check_facets.initialize(mesh_facets, _spatial_dimension = spatial_dimension - 1, _with_nb_element = true, _default_value = true); check_facets.set(true); // remove the pure ghost elements for_each_element( mesh_facets, [&](auto && facet) { const auto & element_to_facet = mesh_facets.getElementToSubelement( facet.type, facet.ghost_type)(facet.element); auto & left = element_to_facet[0]; auto & right = element_to_facet[1]; if (right == ElementNull || (left.ghost_type == _ghost && right.ghost_type == _ghost)) { check_facets(facet) = false; return; } if (left.kind() == _ek_cohesive or right.kind() == _ek_cohesive) { check_facets(facet) = false; } }, _spatial_dimension = spatial_dimension - 1); auto tolerance = Math::getTolerance(); Vector bary_facet(spatial_dimension); // set the limits to the bounding box for_each_element( mesh_facets, [&](auto && facet) { auto & need_check = check_facets(facet); if (not need_check) return; mesh_facets.getBarycenter(facet, bary_facet); UInt coord_in_limit = 0; while (coord_in_limit < spatial_dimension and bary_facet(coord_in_limit) > (insertion_limits(coord_in_limit, 0) - tolerance) and bary_facet(coord_in_limit) < (insertion_limits(coord_in_limit, 1) + tolerance)) ++coord_in_limit; if (coord_in_limit != spatial_dimension) need_check = false; }, _spatial_dimension = spatial_dimension - 1); if (physical_groups.size() == 0) { AKANTU_DEBUG_OUT(); return; } if (not mesh_facets.hasData("physical_names")) { AKANTU_DEBUG_ASSERT( physical_groups.size() == 0, "No physical names in the mesh but insertion limited to a group"); AKANTU_DEBUG_OUT(); return; } const auto & physical_ids = mesh_facets.getData("physical_names"); // set the limits to the physical surfaces for_each_element(mesh_facets, [&](auto && facet) { auto & need_check = check_facets(facet); if (not need_check) return; const auto & physical_id = physical_ids(facet); auto it = find(physical_groups.begin(), physical_groups.end(), physical_id); need_check = (it != physical_groups.end()); }, _spatial_dimension = spatial_dimension - 1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ UInt CohesiveElementInserter::insertElements(bool only_double_facets) { CohesiveNewNodesEvent node_event; NewElementsEvent element_event; Array new_pairs(0, 2); if (mesh_facets.isDistributed()) { mesh_facets.getElementSynchronizer().synchronizeOnce( *this, _gst_ce_groups); } UInt nb_new_elements = MeshUtils::insertCohesiveElements( mesh, mesh_facets, insertion_facets, new_pairs, element_event.getList(), only_double_facets); UInt nb_new_nodes = new_pairs.size(); node_event.getList().reserve(nb_new_nodes); node_event.getOldNodesList().reserve(nb_new_nodes); for (UInt n = 0; n < nb_new_nodes; ++n) { node_event.getList().push_back(new_pairs(n, 1)); node_event.getOldNodesList().push_back(new_pairs(n, 0)); } if (nb_new_elements > 0) { updateInsertionFacets(); } MeshAccessor mesh_accessor(mesh); std::tie(nb_new_nodes, nb_new_elements) = mesh_accessor.updateGlobalData(node_event, element_event); return nb_new_elements; } /* -------------------------------------------------------------------------- */ void CohesiveElementInserter::updateInsertionFacets() { AKANTU_DEBUG_IN(); UInt spatial_dimension = mesh.getSpatialDimension(); for (auto && facet_gt : ghost_types) { for (auto && facet_type : mesh_facets.elementTypes(spatial_dimension - 1, facet_gt)) { auto & ins_facets = insertion_facets(facet_type, facet_gt); // this is the intrinsic case if (not is_extrinsic) continue; auto & f_check = check_facets(facet_type, facet_gt); for (auto && pair : zip(ins_facets, f_check)) { bool & ins = std::get<0>(pair); bool & check = std::get<1>(pair); if (ins) ins = check = false; } } } // resize for the newly added facets insertion_facets.initialize(mesh_facets, _spatial_dimension = spatial_dimension - 1, _with_nb_element = true, _default_value = false); // resize for the newly added facets if (is_extrinsic) { check_facets.initialize(mesh_facets, _spatial_dimension = spatial_dimension - 1, _with_nb_element = true, _default_value = false); } else { insertion_facets.set(false); } AKANTU_DEBUG_OUT(); } } // namespace akantu diff --git a/src/mesh_utils/cohesive_element_inserter.hh b/src/mesh_utils/cohesive_element_inserter.hh index 187125069..4bcac768f 100644 --- a/src/mesh_utils/cohesive_element_inserter.hh +++ b/src/mesh_utils/cohesive_element_inserter.hh @@ -1,171 +1,171 @@ /** * @file cohesive_element_inserter.hh * * @author Fabian Barras * @author Nicolas Richart * @author Marco Vocialta * * @date creation: Wed Dec 04 2013 * @date last modification: Sun Feb 04 2018 * * @brief Cohesive element inserter * * @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 "data_accessor.hh" #include "mesh_utils.hh" #include "parsable.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_COHESIVE_ELEMENT_INSERTER_HH__ #define __AKANTU_COHESIVE_ELEMENT_INSERTER_HH__ namespace akantu { class GlobalIdsUpdater; class FacetSynchronizer; } // akantu namespace akantu { class CohesiveElementInserter : public DataAccessor, public Parsable { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: CohesiveElementInserter(Mesh & mesh, const ID & id = "cohesive_element_inserter"); ~CohesiveElementInserter() override; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// set range limitation for intrinsic cohesive element insertion - void setLimit(SpacialDirection axis, Real first_limit, Real second_limit); + void setLimit(SpatialDirection axis, Real first_limit, Real second_limit); /// insert intrinsic cohesive elements in a predefined range UInt insertIntrinsicElements(); /// insert extrinsic cohesive elements (returns the number of new /// cohesive elements) UInt insertElements(bool only_double_facets = false); /// limit check facets to match given insertion limits void limitCheckFacets(); protected: void parseSection(const ParserSection & section) override; protected: /// internal version of limitCheckFacets void limitCheckFacets(ElementTypeMapArray & check_facets); /// update facet insertion arrays after facets doubling void updateInsertionFacets(); /// functions for parallel communications inline UInt getNbData(const Array & elements, const SynchronizationTag & tag) const override; inline void packData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) const override; inline void unpackData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) override; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: AKANTU_GET_MACRO_NOT_CONST(InsertionFacetsByElement, insertion_facets, ElementTypeMapArray &); AKANTU_GET_MACRO(InsertionFacetsByElement, insertion_facets, const ElementTypeMapArray &); AKANTU_GET_MACRO_BY_ELEMENT_TYPE(InsertionFacets, insertion_facets, bool); AKANTU_GET_MACRO(CheckFacets, check_facets, const ElementTypeMapArray &); AKANTU_GET_MACRO_BY_ELEMENT_TYPE(CheckFacets, check_facets, bool); AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(CheckFacets, check_facets, bool); AKANTU_GET_MACRO(MeshFacets, mesh_facets, const Mesh &); AKANTU_GET_MACRO_NOT_CONST(MeshFacets, mesh_facets, Mesh &); AKANTU_SET_MACRO(IsExtrinsic, is_extrinsic, bool); public: friend class SolidMechanicsModelCohesive; /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: /// object id ID id; /// main mesh where to insert cohesive elements Mesh & mesh; /// mesh containing facets Mesh & mesh_facets; /// list of facets where to insert elements ElementTypeMapArray insertion_facets; /// limits for element insertion Matrix insertion_limits; /// list of groups to consider for insertion, ignored if empty std::vector physical_groups; /// vector containing facets in which extrinsic cohesive elements can be /// inserted ElementTypeMapArray check_facets; /// global connectivity ids updater std::unique_ptr global_ids_updater; /// is this inserter used in extrinsic bool is_extrinsic{false}; }; class CohesiveNewNodesEvent : public NewNodesEvent { public: CohesiveNewNodesEvent() = default; ~CohesiveNewNodesEvent() override = default; AKANTU_GET_MACRO_NOT_CONST(OldNodesList, old_nodes, Array &); AKANTU_GET_MACRO(OldNodesList, old_nodes, const Array &); private: Array old_nodes; }; } // akantu #include "cohesive_element_inserter_inline_impl.cc" #endif /* __AKANTU_COHESIVE_ELEMENT_INSERTER_HH__ */ diff --git a/src/model/solid_mechanics/solid_mechanics_model_cohesive/solid_mechanics_model_cohesive.cc b/src/model/solid_mechanics/solid_mechanics_model_cohesive/solid_mechanics_model_cohesive.cc index 9cf093d9d..c434ffa6d 100644 --- a/src/model/solid_mechanics/solid_mechanics_model_cohesive/solid_mechanics_model_cohesive.cc +++ b/src/model/solid_mechanics/solid_mechanics_model_cohesive/solid_mechanics_model_cohesive.cc @@ -1,707 +1,707 @@ /** * @file solid_mechanics_model_cohesive.cc * * @author Fabian Barras * @author Mauro Corrado * @author Nicolas Richart * @author Marco Vocialta * * @date creation: Tue May 08 2012 * @date last modification: Wed Feb 21 2018 * * @brief Solid mechanics model for cohesive elements * * @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 "solid_mechanics_model_cohesive.hh" #include "aka_iterators.hh" #include "cohesive_element_inserter.hh" #include "element_synchronizer.hh" #include "facet_synchronizer.hh" #include "fe_engine_template.hh" #include "global_ids_updater.hh" #include "integrator_gauss.hh" #include "material_cohesive.hh" #include "mesh_accessor.hh" #include "mesh_global_data_updater.hh" #include "parser.hh" #include "shape_cohesive.hh" /* -------------------------------------------------------------------------- */ #include "dumpable_inline_impl.hh" #ifdef AKANTU_USE_IOHELPER #include "dumper_iohelper_paraview.hh" #endif /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ namespace akantu { class CohesiveMeshGlobalDataUpdater : public MeshGlobalDataUpdater { public: CohesiveMeshGlobalDataUpdater(SolidMechanicsModelCohesive & model) : model(model), mesh(model.getMesh()), global_ids_updater(model.getMesh(), *model.cohesive_synchronizer) {} /* ------------------------------------------------------------------------ */ std::tuple updateData(NewNodesEvent & nodes_event, NewElementsEvent & elements_event) override { auto cohesive_nodes_event = dynamic_cast(&nodes_event); if (not cohesive_nodes_event) return std::make_tuple(nodes_event.getList().size(), elements_event.getList().size()); /// update nodes type auto & new_nodes = cohesive_nodes_event->getList(); auto & old_nodes = cohesive_nodes_event->getOldNodesList(); auto local_nb_new_nodes = new_nodes.size(); auto nb_new_nodes = local_nb_new_nodes; if (mesh.isDistributed()) { MeshAccessor mesh_accessor(mesh); - auto & nodes_type = mesh_accessor.getNodesType(); - UInt nb_old_nodes = nodes_type.size(); - nodes_type.resize(nb_old_nodes + local_nb_new_nodes); + auto & nodes_flags = mesh_accessor.getNodesFlags(); + auto nb_old_nodes = nodes_flags.size(); + nodes_flags.resize(nb_old_nodes + local_nb_new_nodes); for (auto && data : zip(old_nodes, new_nodes)) { UInt old_node, new_node; std::tie(old_node, new_node) = data; - nodes_type(new_node) = nodes_type(old_node); + nodes_flags(new_node) = nodes_flags(old_node); } model.updateCohesiveSynchronizers(); nb_new_nodes = global_ids_updater.updateGlobalIDs(new_nodes.size()); } Vector nb_new_stuff = {nb_new_nodes, elements_event.getList().size()}; const auto & comm = mesh.getCommunicator(); comm.allReduce(nb_new_stuff, SynchronizerOperation::_sum); if (nb_new_stuff(1) > 0) { mesh.sendEvent(elements_event); MeshUtils::resetFacetToDouble(mesh.getMeshFacets()); } if (nb_new_nodes > 0) { mesh.sendEvent(nodes_event); // mesh.sendEvent(global_ids_updater.getChangedNodeEvent()); } return std::make_tuple(nb_new_nodes, nb_new_stuff(1)); } private: SolidMechanicsModelCohesive & model; Mesh & mesh; GlobalIdsUpdater global_ids_updater; }; /* -------------------------------------------------------------------------- */ SolidMechanicsModelCohesive::SolidMechanicsModelCohesive( Mesh & mesh, UInt dim, const ID & id, const MemoryID & memory_id) : SolidMechanicsModel(mesh, dim, id, memory_id, ModelType::_solid_mechanics_model_cohesive), tangents("tangents", id), facet_stress("facet_stress", id), facet_material("facet_material", id) { AKANTU_DEBUG_IN(); registerFEEngineObject("CohesiveFEEngine", mesh, Model::spatial_dimension); auto && tmp_material_selector = std::make_shared(*this); tmp_material_selector->setFallback(this->material_selector); this->material_selector = tmp_material_selector; #if defined(AKANTU_USE_IOHELPER) this->mesh.registerDumper("cohesive elements", id); this->mesh.addDumpMeshToDumper("cohesive elements", mesh, Model::spatial_dimension, _not_ghost, _ek_cohesive); #endif if (this->mesh.isDistributed()) { /// create the distributed synchronizer for cohesive elements this->cohesive_synchronizer = std::make_unique( mesh, "cohesive_distributed_synchronizer"); auto & synchronizer = mesh.getElementSynchronizer(); this->cohesive_synchronizer->split(synchronizer, [](auto && el) { return Mesh::getKind(el.type) == _ek_cohesive; }); this->registerSynchronizer(*cohesive_synchronizer, _gst_material_id); this->registerSynchronizer(*cohesive_synchronizer, _gst_smm_stress); this->registerSynchronizer(*cohesive_synchronizer, _gst_smm_boundary); } this->inserter = std::make_unique( this->mesh, id + ":cohesive_element_inserter"); registerFEEngineObject( "FacetsFEEngine", mesh.getMeshFacets(), Model::spatial_dimension - 1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ SolidMechanicsModelCohesive::~SolidMechanicsModelCohesive() = default; /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::setTimeStep(Real time_step, const ID & solver_id) { SolidMechanicsModel::setTimeStep(time_step, solver_id); #if defined(AKANTU_USE_IOHELPER) this->mesh.getDumper("cohesive elements").setTimeStep(time_step); #endif } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::initFullImpl(const ModelOptions & options) { AKANTU_DEBUG_IN(); const auto & smmc_options = dynamic_cast(options); this->is_extrinsic = smmc_options.is_extrinsic; std::cout << "Extrinsic " << is_extrinsic << std::endl; inserter->setIsExtrinsic(is_extrinsic); if (mesh.isDistributed()) { auto & mesh_facets = inserter->getMeshFacets(); auto & synchronizer = dynamic_cast(mesh_facets.getElementSynchronizer()); synchronizeGhostFacetsConnectivity(); /// create the facet synchronizer for extrinsic simulations if (is_extrinsic) { facet_stress_synchronizer = std::make_unique( synchronizer, id + ":facet_stress_synchronizer"); facet_stress_synchronizer->swapSendRecv(); this->registerSynchronizer(*facet_stress_synchronizer, _gst_smmc_facets_stress); } } MeshAccessor mesh_accessor(mesh); mesh_accessor.registerGlobalDataUpdater( std::make_unique(*this)); ParserSection section; bool is_empty; std::tie(section, is_empty) = this->getParserSection(); if (not is_empty) { auto inserter_section = section.getSubSections(ParserType::_cohesive_inserter); if (inserter_section.begin() != inserter_section.end()) { inserter->parseSection(*inserter_section.begin()); } } SolidMechanicsModel::initFullImpl(options); AKANTU_DEBUG_OUT(); } // namespace akantu /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::initMaterials() { AKANTU_DEBUG_IN(); // make sure the material are instantiated if (!are_materials_instantiated) instantiateMaterials(); /// find the first cohesive material UInt cohesive_index = UInt(-1); for (auto && material : enumerate(materials)) { if (dynamic_cast(std::get<1>(material).get())) { cohesive_index = std::get<0>(material); break; } } if (cohesive_index == UInt(-1)) AKANTU_EXCEPTION("No cohesive materials in the material input file"); material_selector->setFallback(cohesive_index); // set the facet information in the material in case of dynamic insertion // to know what material to call for stress checks const Mesh & mesh_facets = inserter->getMeshFacets(); facet_material.initialize( mesh_facets, _spatial_dimension = spatial_dimension - 1, _with_nb_element = true, _default_value = material_selector->getFallbackValue()); for_each_element( mesh_facets, [&](auto && element) { auto mat_index = (*material_selector)(element); auto & mat = dynamic_cast(*materials[mat_index]); facet_material(element) = mat_index; if (is_extrinsic) { mat.addFacet(element); } }, _spatial_dimension = spatial_dimension - 1, _ghost_type = _not_ghost); SolidMechanicsModel::initMaterials(); if (is_extrinsic) { this->initAutomaticInsertion(); } else { this->insertIntrinsicElements(); } AKANTU_DEBUG_OUT(); } // namespace akantu /* -------------------------------------------------------------------------- */ /** * Initialize the model,basically it pre-compute the shapes, shapes derivatives * and jacobian */ void SolidMechanicsModelCohesive::initModel() { AKANTU_DEBUG_IN(); SolidMechanicsModel::initModel(); /// add cohesive type connectivity ElementType type = _not_defined; for (auto && type_ghost : ghost_types) { for (const auto & tmp_type : mesh.elementTypes(spatial_dimension, type_ghost)) { const auto & connectivity = mesh.getConnectivity(tmp_type, type_ghost); if (connectivity.size() == 0) continue; type = tmp_type; auto type_facet = Mesh::getFacetType(type); auto type_cohesive = FEEngine::getCohesiveElementType(type_facet); mesh.addConnectivityType(type_cohesive, type_ghost); } } AKANTU_DEBUG_ASSERT(type != _not_defined, "No elements in the mesh"); getFEEngine("CohesiveFEEngine").initShapeFunctions(_not_ghost); getFEEngine("CohesiveFEEngine").initShapeFunctions(_ghost); if (is_extrinsic) { getFEEngine("FacetsFEEngine").initShapeFunctions(_not_ghost); getFEEngine("FacetsFEEngine").initShapeFunctions(_ghost); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::insertIntrinsicElements() { AKANTU_DEBUG_IN(); inserter->insertIntrinsicElements(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::initAutomaticInsertion() { AKANTU_DEBUG_IN(); this->inserter->limitCheckFacets(); this->updateFacetStressSynchronizer(); this->resizeFacetStress(); /// compute normals on facets this->computeNormals(); this->initStressInterpolation(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::updateAutomaticInsertion() { AKANTU_DEBUG_IN(); this->inserter->limitCheckFacets(); this->updateFacetStressSynchronizer(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::initStressInterpolation() { Mesh & mesh_facets = inserter->getMeshFacets(); /// compute quadrature points coordinates on facets Array & position = mesh.getNodes(); ElementTypeMapArray quad_facets("quad_facets", id); quad_facets.initialize(mesh_facets, _nb_component = Model::spatial_dimension, _spatial_dimension = Model::spatial_dimension - 1); // mesh_facets.initElementTypeMapArray(quad_facets, Model::spatial_dimension, // Model::spatial_dimension - 1); getFEEngine("FacetsFEEngine") .interpolateOnIntegrationPoints(position, quad_facets); /// compute elements quadrature point positions and build /// element-facet quadrature points data structure ElementTypeMapArray elements_quad_facets("elements_quad_facets", id); elements_quad_facets.initialize( mesh, _nb_component = Model::spatial_dimension, _spatial_dimension = Model::spatial_dimension); // mesh.initElementTypeMapArray(elements_quad_facets, // Model::spatial_dimension, // Model::spatial_dimension); for (auto elem_gt : ghost_types) { for (auto & type : mesh.elementTypes(Model::spatial_dimension, elem_gt)) { UInt nb_element = mesh.getNbElement(type, elem_gt); if (nb_element == 0) continue; /// compute elements' quadrature points and list of facet /// quadrature points positions by element const auto & facet_to_element = mesh_facets.getSubelementToElement(type, elem_gt); auto & el_q_facet = elements_quad_facets(type, elem_gt); auto facet_type = Mesh::getFacetType(type); auto nb_quad_per_facet = getFEEngine("FacetsFEEngine").getNbIntegrationPoints(facet_type); auto nb_facet_per_elem = facet_to_element.getNbComponent(); // small hack in the loop to skip boundary elements, they are silently // initialized to NaN to see if this causes problems el_q_facet.resize(nb_element * nb_facet_per_elem * nb_quad_per_facet, std::numeric_limits::quiet_NaN()); for (auto && data : zip(make_view(facet_to_element), make_view(el_q_facet, spatial_dimension, nb_quad_per_facet))) { const auto & global_facet = std::get<0>(data); auto & el_q = std::get<1>(data); if (global_facet == ElementNull) continue; Matrix quad_f = make_view(quad_facets(global_facet.type, global_facet.ghost_type), spatial_dimension, nb_quad_per_facet) .begin()[global_facet.element]; el_q = quad_f; // for (UInt q = 0; q < nb_quad_per_facet; ++q) { // for (UInt s = 0; s < Model::spatial_dimension; ++s) { // el_q_facet(el * nb_facet_per_elem * nb_quad_per_facet + // f * nb_quad_per_facet + q, // s) = quad_f(global_facet * nb_quad_per_facet + q, // s); // } // } //} } } } /// loop over non cohesive materials for (auto && material : materials) { if (dynamic_cast(material.get())) continue; /// initialize the interpolation function material->initElementalFieldInterpolation(elements_quad_facets); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::assembleInternalForces() { AKANTU_DEBUG_IN(); // f_int += f_int_cohe for (auto & material : this->materials) { try { auto & mat = dynamic_cast(*material); mat.computeTraction(_not_ghost); } catch (std::bad_cast & bce) { } } SolidMechanicsModel::assembleInternalForces(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::computeNormals() { AKANTU_DEBUG_IN(); Mesh & mesh_facets = this->inserter->getMeshFacets(); this->getFEEngine("FacetsFEEngine") .computeNormalsOnIntegrationPoints(_not_ghost); /** * @todo store tangents while computing normals instead of * recomputing them as follows: */ /* ------------------------------------------------------------------------ */ UInt tangent_components = Model::spatial_dimension * (Model::spatial_dimension - 1); tangents.initialize(mesh_facets, _nb_component = tangent_components, _spatial_dimension = Model::spatial_dimension - 1); // mesh_facets.initElementTypeMapArray(tangents, tangent_components, // Model::spatial_dimension - 1); for (auto facet_type : mesh_facets.elementTypes(Model::spatial_dimension - 1)) { const Array & normals = this->getFEEngine("FacetsFEEngine") .getNormalsOnIntegrationPoints(facet_type); Array & tangents = this->tangents(facet_type); Math::compute_tangents(normals, tangents); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::interpolateStress() { ElementTypeMapArray by_elem_result("temporary_stress_by_facets", id); for (auto & material : materials) { auto * mat = dynamic_cast(material.get()); if (mat == nullptr) /// interpolate stress on facet quadrature points positions material->interpolateStressOnFacets(facet_stress, by_elem_result); } this->synchronize(_gst_smmc_facets_stress); } /* -------------------------------------------------------------------------- */ UInt SolidMechanicsModelCohesive::checkCohesiveStress() { AKANTU_DEBUG_IN(); if (not is_extrinsic) { AKANTU_EXCEPTION( "This function can only be used for extrinsic cohesive elements"); } interpolateStress(); for (auto & mat : materials) { auto * mat_cohesive = dynamic_cast(mat.get()); if (mat_cohesive) { /// check which not ghost cohesive elements are to be created mat_cohesive->checkInsertion(); } } /// communicate data among processors // this->synchronize(_gst_smmc_facets); /// insert cohesive elements UInt nb_new_elements = inserter->insertElements(); // if (nb_new_elements > 0) { // this->reinitializeSolver(); // } AKANTU_DEBUG_OUT(); return nb_new_elements; } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::onElementsAdded( const Array & element_list, const NewElementsEvent & event) { AKANTU_DEBUG_IN(); SolidMechanicsModel::onElementsAdded(element_list, event); if (is_extrinsic) resizeFacetStress(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::onNodesAdded(const Array & new_nodes, const NewNodesEvent & event) { AKANTU_DEBUG_IN(); SolidMechanicsModel::onNodesAdded(new_nodes, event); const CohesiveNewNodesEvent * cohesive_event; if ((cohesive_event = dynamic_cast(&event)) == nullptr) return; const auto & old_nodes = cohesive_event->getOldNodesList(); auto copy = [this, &new_nodes, &old_nodes](auto & arr) { UInt new_node, old_node; auto view = make_view(arr, spatial_dimension); auto begin = view.begin(); for (auto && pair : zip(new_nodes, old_nodes)) { std::tie(new_node, old_node) = pair; auto old_ = begin + old_node; auto new_ = begin + new_node; *new_ = *old_; } }; copy(*displacement); copy(*blocked_dofs); if (velocity) copy(*velocity); if (acceleration) copy(*acceleration); if (current_position) copy(*current_position); if (previous_displacement) copy(*previous_displacement); // if (external_force) // copy(*external_force); // if (internal_force) // copy(*internal_force); if (displacement_increment) copy(*displacement_increment); copy(getDOFManager().getSolution("displacement")); // this->assembleMassLumped(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::afterSolveStep() { AKANTU_DEBUG_IN(); /* * This is required because the Cauchy stress is the stress measure that * is used to check the insertion of cohesive elements */ for (auto & mat : materials) { if (mat->isFiniteDeformation()) mat->computeAllCauchyStresses(_not_ghost); } SolidMechanicsModel::afterSolveStep(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::printself(std::ostream & stream, int indent) const { std::string space; for (Int i = 0; i < indent; i++, space += AKANTU_INDENT) ; stream << space << "SolidMechanicsModelCohesive [" << std::endl; SolidMechanicsModel::printself(stream, indent + 1); stream << space << "]" << std::endl; } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::resizeFacetStress() { AKANTU_DEBUG_IN(); this->facet_stress.initialize(getFEEngine("FacetsFEEngine"), _nb_component = 2 * spatial_dimension * spatial_dimension, _spatial_dimension = spatial_dimension - 1); // for (auto && ghost_type : ghost_types) { // for (const auto & type : // mesh_facets.elementTypes(spatial_dimension - 1, ghost_type)) { // UInt nb_facet = mesh_facets.getNbElement(type, ghost_type); // UInt nb_quadrature_points = getFEEngine("FacetsFEEngine") // .getNbIntegrationPoints(type, // ghost_type); // UInt new_size = nb_facet * nb_quadrature_points; // facet_stress(type, ghost_type).resize(new_size); // } // } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::addDumpGroupFieldToDumper( const std::string & dumper_name, const std::string & field_id, const std::string & group_name, const ElementKind & element_kind, bool padding_flag) { AKANTU_DEBUG_IN(); UInt spatial_dimension = Model::spatial_dimension; ElementKind _element_kind = element_kind; if (dumper_name == "cohesive elements") { _element_kind = _ek_cohesive; } else if (dumper_name == "facets") { spatial_dimension = Model::spatial_dimension - 1; } SolidMechanicsModel::addDumpGroupFieldToDumper(dumper_name, field_id, group_name, spatial_dimension, _element_kind, padding_flag); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::onDump() { this->flattenAllRegisteredInternals(_ek_cohesive); SolidMechanicsModel::onDump(); } /* -------------------------------------------------------------------------- */ } // namespace akantu 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 d7c1b4688..d6d1703b0 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,335 +1,335 @@ /** * @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); if (dim == 1) { surface = 1; 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(); auto & mat_el = this->model->getMaterial("body"); auto speed = mat_el.getPushWaveSpeed(Element()); auto direction = _y; if(dim == 1) direction = _x; auto length = mesh->getUpperBounds()(direction) - mesh->getLowerBounds()(direction); nb_steps = length / 2. / speed / model->getTimeStep(); 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"); 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() { - EXPECT_NO_THROW(this->createModel()); + ASSERT_NO_THROW(this->createModel()); auto & mat_el = this->model->getMaterial("body"); Real speed; try { speed = mat_el.getShearWaveSpeed(Element()); // the slowest speed if exists } catch(...) { speed = mat_el.getPushWaveSpeed(Element()); } auto direction = _y; if(dim == 1) direction = _x; auto length = mesh->getUpperBounds()(direction) - mesh->getLowerBounds()(direction); nb_steps = length / 2. / speed / model->getTimeStep(); SCOPED_TRACE(std::to_string(this->dim) + "D - " + aka::to_string(type_1) + ":" + aka::to_string(type_2)); 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"); 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__ */