diff --git a/src/mesh_utils/cohesive_element_inserter.cc b/src/mesh_utils/cohesive_element_inserter.cc index 221683278..d71b9bfe9 100644 --- a/src/mesh_utils/cohesive_element_inserter.cc +++ b/src/mesh_utils/cohesive_element_inserter.cc @@ -1,393 +1,394 @@ /** * @file cohesive_element_inserter.cc * * @author Marco Vocialta * * @date creation: Wed Dec 04 2013 * @date last modification: Sun Oct 04 2015 * * @brief Cohesive element inserter functions * * @section LICENSE * * Copyright (©) 2014, 2015 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_iterators.hh" /* -------------------------------------------------------------------------- */ #include #include /* -------------------------------------------------------------------------- */ namespace akantu { CohesiveElementInserter::CohesiveElementInserter(Mesh & mesh, bool is_extrinsic, 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("groups", physical_groups, _pat_parsable, "List of groups to consider for insertion"); this->registerParam("bounding_box", insertion_limits, _pat_parsable, "Global limit for insertion"); init(is_extrinsic); } /* -------------------------------------------------------------------------- */ CohesiveElementInserter::~CohesiveElementInserter() = default; /* -------------------------------------------------------------------------- */ void CohesiveElementInserter::init(bool is_extrinsic) { AKANTU_DEBUG_IN(); UInt spatial_dimension = mesh.getSpatialDimension(); MeshUtils::resetFacetToDouble(mesh_facets); /// initialize facet insertion array insertion_facets.initialize(mesh_facets, _spatial_dimension = (spatial_dimension - 1), _with_nb_element = true); /// init insertion limits for (UInt dim = 0; dim < spatial_dimension; ++dim) { insertion_limits(dim, 0) = std::numeric_limits::max() * (-1.); insertion_limits(dim, 1) = std::numeric_limits::max(); } if (is_extrinsic) { initFacetsCheck(); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void CohesiveElementInserter::initFacetsCheck() { 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); for_each_elements( 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; } }, _spatial_dimension = spatial_dimension - 1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void CohesiveElementInserter::limitCheckFacets() { AKANTU_DEBUG_IN(); UInt spatial_dimension = mesh.getSpatialDimension(); Vector bary_facet(spatial_dimension); for_each_elements(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 && bary_facet(coord_in_limit) > insertion_limits(coord_in_limit, 0) && bary_facet(coord_in_limit) < insertion_limits(coord_in_limit, 1)) ++coord_in_limit; if (coord_in_limit != spatial_dimension) need_check = false; }, _spatial_dimension = spatial_dimension - 1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void CohesiveElementInserter::setLimit(SpacialDirection 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() { AKANTU_DEBUG_IN(); UInt spatial_dimension = mesh.getSpatialDimension(); Vector bary_facet(spatial_dimension); for_each_elements( mesh_facets, [&](auto && facet) { const auto & element_to_facet = mesh_facets.getElementToSubelement( facet.type, facet.ghost_type)(facet.element); if (element_to_facet[1] == ElementNull) return; mesh_facets.getBarycenter(facet, bary_facet); UInt coord_in_limit = 0; while (coord_in_limit < spatial_dimension && bary_facet(coord_in_limit) > insertion_limits(coord_in_limit, 0) && bary_facet(coord_in_limit) < insertion_limits(coord_in_limit, 1)) ++coord_in_limit; if (coord_in_limit == spatial_dimension) insertion_facets(facet) = true; }, _spatial_dimension = spatial_dimension - 1); AKANTU_DEBUG_OUT(); return insertElements(); } /* -------------------------------------------------------------------------- */ void CohesiveElementInserter::insertIntrinsicElements( const std::string & physname, UInt material_index) { AKANTU_DEBUG_IN(); UInt spatial_dimension = mesh.getSpatialDimension(); ElementTypeMapArray * phys_data; try { phys_data = &(mesh_facets.getData("physical_names")); } catch (...) { phys_data = &(mesh_facets.registerData("physical_names")); phys_data->initialize(mesh_facets, _spatial_dimension = spatial_dimension - 1, _with_nb_element = true); // mesh_facets.initElementTypeMapArray(*phys_data, 1, spatial_dimension - // 1, // false, _ek_regular, true); } Vector bary_facet(spatial_dimension); mesh_facets.createElementGroup(physname); GhostType ghost_type = _not_ghost; for (auto && type_facet : mesh_facets.elementTypes(spatial_dimension - 1, ghost_type)) { auto & f_insertion = insertion_facets(type_facet, ghost_type); auto & element_to_facet = mesh_facets.getElementToSubelement(type_facet, ghost_type); UInt nb_facet = mesh_facets.getNbElement(type_facet, ghost_type); UInt coord_in_limit = 0; auto & group = mesh.getElementGroup(physname); auto & group_facet = mesh_facets.getElementGroup(physname); Vector bary_physgroup(spatial_dimension); Real norm_bary; for (auto && e : group.getElements(type_facet, ghost_type)) { mesh.getBarycenter(Element{type_facet, e, ghost_type}, bary_physgroup); #ifndef AKANTU_NDEBUG bool find_a_partner = false; #endif norm_bary = bary_physgroup.norm(); Array & material_id = (*phys_data)(type_facet, ghost_type); for (UInt f = 0; f < nb_facet; ++f) { if (element_to_facet(f)[1] == ElementNull) continue; mesh_facets.getBarycenter(Element{type_facet, f, ghost_type}, bary_facet); coord_in_limit = 0; while (coord_in_limit < spatial_dimension && (std::abs(bary_facet(coord_in_limit) - bary_physgroup(coord_in_limit)) / norm_bary < Math::getTolerance())) ++coord_in_limit; if (coord_in_limit == spatial_dimension) { f_insertion(f) = true; #ifndef AKANTU_NDEBUG find_a_partner = true; #endif group_facet.add(type_facet, f, ghost_type, false); material_id(f) = material_index; break; } } AKANTU_DEBUG_ASSERT(find_a_partner, "The element nO " << e << " of physical group " << physname << " did not find its associated facet!" << " Try to decrease math tolerance. " << std::endl); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ UInt CohesiveElementInserter::insertElements(bool only_double_facets) { CohesiveNewNodesEvent node_event; NewElementsEvent element_event; Array new_pairs(0, 2); 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 (mesh.isDistributed()) { /// update nodes type updateNodesType(mesh, node_event); updateNodesType(mesh_facets, node_event); /// update global ids nb_new_nodes = this->updateGlobalIDs(node_event); /// compute total number of new elements const auto & comm = mesh.getCommunicator(); comm.allReduce(nb_new_elements, SynchronizerOperation::_sum); } if (nb_new_nodes > 0) mesh.sendEvent(node_event); if (nb_new_elements > 0) { updateInsertionFacets(); // mesh.updateTypesOffsets(_not_ghost); mesh.sendEvent(element_event); MeshUtils::resetFacetToDouble(mesh_facets); } if (mesh.isDistributed()) { /// update global ids this->synchronizeGlobalIDs(node_event); } return nb_new_elements; } /* -------------------------------------------------------------------------- */ void CohesiveElementInserter::updateInsertionFacets() { AKANTU_DEBUG_IN(); UInt spatial_dimension = mesh.getSpatialDimension(); for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) { GhostType facet_gt = *gt; Mesh::type_iterator it = mesh_facets.firstType(spatial_dimension - 1, facet_gt); Mesh::type_iterator last = mesh_facets.lastType(spatial_dimension - 1, facet_gt); for (; it != last; ++it) { ElementType facet_type = *it; Array & ins_facets = insertion_facets(facet_type, facet_gt); // this is the extrinsic case if (check_facets.exists(facet_type, facet_gt)) { Array & f_check = check_facets(facet_type, facet_gt); UInt nb_facets = f_check.size(); for (UInt f = 0; f < ins_facets.size(); ++f) { if (ins_facets(f)) { ++nb_facets; ins_facets(f) = false; f_check(f) = false; } } f_check.resize(nb_facets); } else { // and this the intrinsic one ins_facets.resize(mesh_facets.getNbElement(facet_type, facet_gt)); ins_facets.set(false); } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void CohesiveElementInserter::printself(std::ostream & stream, int indent) const { std::string space; for (Int i = 0; i < indent; i++, space += AKANTU_INDENT) ; stream << space << "CohesiveElementInserter [" << std::endl; stream << space << " + mesh [" << std::endl; mesh.printself(stream, indent + 2); stream << space << AKANTU_INDENT << "]" << std::endl; stream << space << " + mesh_facets [" << std::endl; mesh_facets.printself(stream, indent + 2); stream << space << AKANTU_INDENT << "]" << std::endl; stream << space << "]" << std::endl; } } // namespace akantu diff --git a/src/model/dof_manager.hh b/src/model/dof_manager.hh index 82e8fd0a2..47b4516ef 100644 --- a/src/model/dof_manager.hh +++ b/src/model/dof_manager.hh @@ -1,467 +1,468 @@ /** * @file dof_manager.hh * * @author Nicolas Richart * * @date Wed Jul 22 11:43:43 2015 * * @brief Class handling the different types of dofs * * @section LICENSE * * Copyright (©) 2010-2011 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 "aka_memory.hh" #include "aka_factory.hh" #include "mesh.hh" /* -------------------------------------------------------------------------- */ #include #include /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_DOF_MANAGER_HH__ #define __AKANTU_DOF_MANAGER_HH__ namespace akantu { class TermsToAssemble; class NonLinearSolver; class TimeStepSolver; class SparseMatrix; } // namespace akantu namespace akantu { class DOFManager : protected Memory, protected MeshEventHandler { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: DOFManager(const ID & id = "dof_manager", const MemoryID & memory_id = 0); DOFManager(Mesh & mesh, const ID & id = "dof_manager", const MemoryID & memory_id = 0); ~DOFManager() override; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ private: /// common function to help registering dofs void registerDOFsInternal(const ID & dof_id, Array & dofs_array); public: /// register an array of degree of freedom virtual void registerDOFs(const ID & dof_id, Array & dofs_array, const DOFSupportType & support_type); /// the dof as an implied type of _dst_nodal and is defined only on a subset /// of nodes virtual void registerDOFs(const ID & dof_id, Array & dofs_array, const ID & group_support); /// register an array of previous values of the degree of freedom virtual void registerDOFsPrevious(const ID & dof_id, Array & dofs_array); /// register an array of increment of degree of freedom virtual void registerDOFsIncrement(const ID & dof_id, Array & dofs_array); /// register an array of derivatives for a particular dof array virtual void registerDOFsDerivative(const ID & dof_id, UInt order, Array & dofs_derivative); /// register array representing the blocked degree of freedoms virtual void registerBlockedDOFs(const ID & dof_id, Array & blocked_dofs); /// Assemble an array to the global residual array virtual void assembleToResidual(const ID & dof_id, const Array & array_to_assemble, Real scale_factor = 1.) = 0; /// Assemble an array to the global lumped matrix array virtual void assembleToLumpedMatrix(const ID & dof_id, const Array & array_to_assemble, const ID & lumped_mtx, Real scale_factor = 1.) = 0; /** * Assemble elementary values to a local array of the size nb_nodes * * nb_dof_per_node. The dof number is implicitly considered as * conn(el, n) * nb_nodes_per_element + d. * With 0 < n < nb_nodes_per_element and 0 < d < nb_dof_per_node **/ virtual void assembleElementalArrayLocalArray( const Array & elementary_vect, Array & array_assembeled, const ElementType & type, const GhostType & ghost_type, Real scale_factor = 1., const Array & filter_elements = empty_filter); /** * Assemble elementary values to the global residual array. The dof number is * implicitly considered as conn(el, n) * nb_nodes_per_element + d. * With 0 < n < nb_nodes_per_element and 0 < d < nb_dof_per_node **/ virtual void assembleElementalArrayToResidual( const ID & dof_id, const Array & elementary_vect, const ElementType & type, const GhostType & ghost_type, Real scale_factor = 1., const Array & filter_elements = empty_filter); /** * Assemble elementary values to a global array corresponding to a lumped * matrix */ virtual void assembleElementalArrayToLumpedMatrix( const ID & dof_id, const Array & elementary_vect, const ID & lumped_mtx, const ElementType & type, const GhostType & ghost_type, Real scale_factor = 1., const Array & filter_elements = empty_filter); /** * Assemble elementary values to the global residual array. The dof number is * implicitly considered as conn(el, n) * nb_nodes_per_element + d. With 0 < * n < nb_nodes_per_element and 0 < d < nb_dof_per_node **/ virtual void assembleElementalMatricesToMatrix( const ID & matrix_id, const ID & dof_id, const Array & elementary_mat, const ElementType & type, const GhostType & ghost_type = _not_ghost, const MatrixType & elemental_matrix_type = _symmetric, const Array & filter_elements = empty_filter) = 0; /// multiply a vector by a matrix and assemble the result to the residual virtual void assembleMatMulVectToResidual(const ID & dof_id, const ID & A_id, const Array & x, Real scale_factor = 1) = 0; /// multiply a vector by a lumped matrix and assemble the result to the /// residual virtual void assembleLumpedMatMulVectToResidual(const ID & dof_id, const ID & A_id, const Array & x, Real scale_factor = 1) = 0; /// assemble coupling terms between to dofs virtual void assemblePreassembledMatrix(const ID & dof_id_m, const ID & dof_id_n, const ID & matrix_id, const TermsToAssemble & terms) = 0; /// sets the residual to 0 virtual void clearResidual() = 0; /// sets the matrix to 0 virtual void clearMatrix(const ID & mtx) = 0; /// sets the lumped matrix to 0 virtual void clearLumpedMatrix(const ID & mtx) = 0; /// splits the solution storage from a global view to the per dof storages void splitSolutionPerDOFs(); /// extract a lumped matrix part corresponding to a given dof virtual void getLumpedMatrixPerDOFs(const ID & dof_id, const ID & lumped_mtx, Array & lumped) = 0; protected: /// minimum functionality to implement per derived version of the DOFManager /// to allow the splitSolutionPerDOFs function to work virtual void getSolutionPerDOFs(const ID & dof_id, Array & solution_array) = 0; protected: /* ------------------------------------------------------------------------ */ /// register a matrix SparseMatrix & registerSparseMatrix(const ID & matrix_id, std::unique_ptr & matrix); /// register a non linear solver instantiated by a derived class NonLinearSolver & registerNonLinearSolver(const ID & non_linear_solver_id, std::unique_ptr & non_linear_solver); /// register a time step solver instantiated by a derived class TimeStepSolver & registerTimeStepSolver(const ID & time_step_solver_id, std::unique_ptr & time_step_solver); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// Get the equation numbers corresponding to a dof_id. This might be used to /// access the matrix. inline const Array & getEquationsNumbers(const ID & dof_id) const; /// Global number of dofs AKANTU_GET_MACRO(SystemSize, this->system_size, UInt); /// Local number of dofs AKANTU_GET_MACRO(LocalSystemSize, this->local_system_size, UInt); /// Retrieve all the registered DOFs std::vector getDOFIDs() const; /* ------------------------------------------------------------------------ */ /* DOFs and derivatives accessors */ /* ------------------------------------------------------------------------ */ /// Get a reference to the registered dof array for a given id inline Array & getDOFs(const ID & dofs_id); /// Get the support type of a given dof inline DOFSupportType getSupportType(const ID & dofs_id) const; /// are the dofs registered inline bool hasDOFs(const ID & dofs_id) const; /// Get a reference to the registered dof derivatives array for a given id inline Array & getDOFsDerivatives(const ID & dofs_id, UInt order); /// Does the dof has derivatives inline bool hasDOFsDerivatives(const ID & dofs_id, UInt order) const; /// Get a reference to the blocked dofs array registered for the given id inline const Array & getBlockedDOFs(const ID & dofs_id) const; /// Does the dof has a blocked array inline bool hasBlockedDOFs(const ID & dofs_id) const; /// Get a reference to the registered dof increment array for a given id inline Array & getDOFsIncrement(const ID & dofs_id); /// Does the dof has a increment array inline bool hasDOFsIncrement(const ID & dofs_id) const; /// Does the dof has a previous array inline Array & getPreviousDOFs(const ID & dofs_id); /// Get a reference to the registered dof array for previous step values a /// given id inline bool hasPreviousDOFs(const ID & dofs_id) const; /// saves the values from dofs to previous dofs virtual void savePreviousDOFs(const ID & dofs_id); /// Get a reference to the solution array registered for the given id inline const Array & getSolution(const ID & dofs_id) const; /* ------------------------------------------------------------------------ */ /* Matrices accessors */ /* ------------------------------------------------------------------------ */ /// Get an instance of a new SparseMatrix virtual SparseMatrix & getNewMatrix(const ID & matrix_id, const MatrixType & matrix_type) = 0; /// Get an instance of a new SparseMatrix as a copy of the SparseMatrix /// matrix_to_copy_id virtual SparseMatrix & getNewMatrix(const ID & matrix_id, const ID & matrix_to_copy_id) = 0; /// Get the reference of an existing matrix SparseMatrix & getMatrix(const ID & matrix_id); /// check if the given matrix exists bool hasMatrix(const ID & matrix_id) const; /// Get an instance of a new lumped matrix virtual Array & getNewLumpedMatrix(const ID & matrix_id); /// Get the lumped version of a given matrix const Array & getLumpedMatrix(const ID & matrix_id) const; /// Get the lumped version of a given matrix Array & getLumpedMatrix(const ID & matrix_id); /// check if the given matrix exists bool hasLumpedMatrix(const ID & matrix_id) const; /* ------------------------------------------------------------------------ */ /* Non linear system solver */ /* ------------------------------------------------------------------------ */ /// Get instance of a non linear solver virtual NonLinearSolver & getNewNonLinearSolver( const ID & nls_solver_id, const NonLinearSolverType & _non_linear_solver_type) = 0; /// get instance of a non linear solver virtual NonLinearSolver & getNonLinearSolver(const ID & nls_solver_id); /// check if the given solver exists bool hasNonLinearSolver(const ID & solver_id) const; /* ------------------------------------------------------------------------ */ /* Time-Step Solver */ /* ------------------------------------------------------------------------ */ /// Get instance of a time step solver virtual TimeStepSolver & getNewTimeStepSolver(const ID & time_step_solver_id, const TimeStepSolverType & type, NonLinearSolver & non_linear_solver) = 0; /// get instance of a time step solver virtual TimeStepSolver & getTimeStepSolver(const ID & time_step_solver_id); /// check if the given solver exists bool hasTimeStepSolver(const ID & solver_id) const; /* ------------------------------------------------------------------------ */ const Mesh & getMesh() { if (mesh) { return *mesh; } else { AKANTU_EXCEPTION("No mesh registered in this dof manager"); } } /* ------------------------------------------------------------------------ */ AKANTU_GET_MACRO(Communicator, communicator, const auto &); AKANTU_GET_MACRO_NOT_CONST(Communicator, communicator, auto &); /* ------------------------------------------------------------------------ */ /* MeshEventHandler interface */ /* ------------------------------------------------------------------------ */ protected: /// helper function for the DOFManager::onNodesAdded method virtual std::pair updateNodalDOFs(const ID & dof_id, const Array & nodes_list); public: /// function to implement to react on akantu::NewNodesEvent void onNodesAdded(const Array & nodes_list, const NewNodesEvent & event) override; /// function to implement to react on akantu::RemovedNodesEvent void onNodesRemoved(const Array & nodes_list, const Array & new_numbering, const RemovedNodesEvent & event) override; /// function to implement to react on akantu::NewElementsEvent void onElementsAdded(const Array & elements_list, const NewElementsEvent & event) override; /// function to implement to react on akantu::RemovedElementsEvent void onElementsRemoved(const Array & elements_list, const ElementTypeMapArray & new_numbering, const RemovedElementsEvent & event) override; /// function to implement to react on akantu::ChangedElementsEvent void onElementsChanged(const Array & old_elements_list, const Array & new_elements_list, const ElementTypeMapArray & new_numbering, const ChangedElementsEvent & event) override; protected: struct DOFData; inline DOFData & getDOFData(const ID & dof_id); inline const DOFData & getDOFData(const ID & dof_id) const; template inline _DOFData & getDOFDataTyped(const ID & dof_id); template inline const _DOFData & getDOFDataTyped(const ID & dof_id) const; virtual DOFData & getNewDOFData(const ID & dof_id); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// dof representations in the dof manager struct DOFData { DOFData() = delete; explicit DOFData(const ID & dof_id); virtual ~DOFData(); /// DOF support type (nodal, general) this is needed to determine how the /// dof are shared among processors DOFSupportType support_type; ID group_support; /// Degree of freedom array Array * dof; /// Blocked degree of freedoms array Array * blocked_dofs; /// Degree of freedoms increment Array * increment; /// Degree of freedoms at previous step Array * previous; /// Solution associated to the dof Array solution; /// local numbering equation numbers Array local_equation_number; /* ---------------------------------------------------------------------- */ /* data for dynamic simulations */ /* ---------------------------------------------------------------------- */ /// Degree of freedom derivatives arrays std::vector *> dof_derivatives; }; /// type to store dofs information using DOFStorage = std::map>; /// type to store all the matrices using SparseMatricesMap = std::map>; /// type to store all the lumped matrices using LumpedMatricesMap = std::map>>; /// type to store all the non linear solver using NonLinearSolversMap = std::map>; /// type to store all the time step solver using TimeStepSolversMap = std::map>; /// store a reference to the dof arrays DOFStorage dofs; /// list of sparse matrices that where created SparseMatricesMap matrices; /// list of lumped matrices LumpedMatricesMap lumped_matrices; /// non linear solvers storage NonLinearSolversMap non_linear_solvers; /// time step solvers storage TimeStepSolversMap time_step_solvers; /// reference to the underlying mesh Mesh * mesh{nullptr}; /// Total number of degrees of freedom (size with the ghosts) UInt local_system_size{0}; /// Number of purely local dofs (size without the ghosts) UInt pure_local_system_size{0}; /// Total number of degrees of freedom UInt system_size{0}; /// Communicator used for this manager, should be the same as in the mesh if a /// mesh is registered Communicator & communicator; }; using DefaultDOFManagerFactory = Factory; +using DOFManagerFactory = Factory; } // namespace akantu #include "dof_manager_inline_impl.cc" #endif /* __AKANTU_DOF_MANAGER_HH__ */ diff --git a/src/model/dof_manager_default.cc b/src/model/dof_manager_default.cc index 3d5323a16..93a7fe339 100644 --- a/src/model/dof_manager_default.cc +++ b/src/model/dof_manager_default.cc @@ -1,910 +1,919 @@ /** * @file dof_manager_default.cc * * @author Nicolas Richart * * @date Tue Aug 11 16:21:01 2015 * * @brief Implementation of the default DOFManager * * @section LICENSE * * Copyright (©) 2010-2011 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 "dof_manager_default.hh" +#include "communicator.hh" #include "dof_synchronizer.hh" #include "element_group.hh" #include "node_synchronizer.hh" #include "non_linear_solver_default.hh" #include "sparse_matrix_aij.hh" -#include "communicator.hh" #include "terms_to_assemble.hh" #include "time_step_solver_default.hh" /* -------------------------------------------------------------------------- */ #include #include #include #include /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ inline void DOFManagerDefault::addSymmetricElementalMatrixToSymmetric( SparseMatrixAIJ & matrix, const Matrix & elementary_mat, const Vector & equation_numbers, UInt max_size) { for (UInt i = 0; i < elementary_mat.rows(); ++i) { UInt c_irn = equation_numbers(i); if (c_irn < max_size) { for (UInt j = i; j < elementary_mat.cols(); ++j) { UInt c_jcn = equation_numbers(j); if (c_jcn < max_size) { matrix(c_irn, c_jcn) += elementary_mat(i, j); } } } } } /* -------------------------------------------------------------------------- */ inline void DOFManagerDefault::addUnsymmetricElementalMatrixToSymmetric( SparseMatrixAIJ & matrix, const Matrix & elementary_mat, const Vector & equation_numbers, UInt max_size) { for (UInt i = 0; i < elementary_mat.rows(); ++i) { UInt c_irn = equation_numbers(i); if (c_irn < max_size) { for (UInt j = 0; j < elementary_mat.cols(); ++j) { UInt c_jcn = equation_numbers(j); if (c_jcn < max_size) { if (c_jcn >= c_irn) { matrix(c_irn, c_jcn) += elementary_mat(i, j); } } } } } } /* -------------------------------------------------------------------------- */ inline void DOFManagerDefault::addElementalMatrixToUnsymmetric( SparseMatrixAIJ & matrix, const Matrix & elementary_mat, const Vector & equation_numbers, UInt max_size) { for (UInt i = 0; i < elementary_mat.rows(); ++i) { UInt c_irn = equation_numbers(i); if (c_irn < max_size) { for (UInt j = 0; j < elementary_mat.cols(); ++j) { UInt c_jcn = equation_numbers(j); if (c_jcn < max_size) { matrix(c_irn, c_jcn) += elementary_mat(i, j); } } } } } /* -------------------------------------------------------------------------- */ DOFManagerDefault::DOFManagerDefault(const ID & id, const MemoryID & memory_id) : DOFManager(id, memory_id), residual(0, 1, std::string(id + ":residual")), global_residual(nullptr), global_solution(0, 1, std::string(id + ":global_solution")), global_blocked_dofs(0, 1, std::string(id + ":global_blocked_dofs")), previous_global_blocked_dofs( 0, 1, std::string(id + ":previous_global_blocked_dofs")), dofs_type(0, 1, std::string(id + ":dofs_type")), data_cache(0, 1, std::string(id + ":data_cache_array")), global_equation_number(0, 1, "global_equation_number"), synchronizer(nullptr) {} /* -------------------------------------------------------------------------- */ DOFManagerDefault::DOFManagerDefault(Mesh & mesh, const ID & id, const MemoryID & memory_id) : DOFManager(mesh, id, memory_id), residual(0, 1, std::string(id + ":residual")), global_residual(nullptr), global_solution(0, 1, std::string(id + ":global_solution")), global_blocked_dofs(0, 1, std::string(id + ":global_blocked_dofs")), previous_global_blocked_dofs( 0, 1, std::string(id + ":previous_global_blocked_dofs")), dofs_type(0, 1, std::string(id + ":dofs_type")), data_cache(0, 1, std::string(id + ":data_cache_array")), jacobian_release(0), global_equation_number(0, 1, "global_equation_number"), first_global_dof_id(0), synchronizer(nullptr) { if (this->mesh->isDistributed()) this->synchronizer = std::make_unique( *this, this->id + ":dof_synchronizer", this->memory_id); } /* -------------------------------------------------------------------------- */ DOFManagerDefault::~DOFManagerDefault() = default; /* -------------------------------------------------------------------------- */ template void DOFManagerDefault::assembleToGlobalArray( const ID & dof_id, const Array & array_to_assemble, Array & global_array, T scale_factor) { AKANTU_DEBUG_IN(); const Array & equation_number = this->getLocalEquationNumbers(dof_id); UInt nb_degree_of_freedoms = array_to_assemble.size() * array_to_assemble.getNbComponent(); AKANTU_DEBUG_ASSERT(equation_number.size() == nb_degree_of_freedoms, "The array to assemble does not have a correct size." << " (" << array_to_assemble.getID() << ")"); typename Array::const_scalar_iterator arr_it = array_to_assemble.begin_reinterpret(nb_degree_of_freedoms); Array::const_scalar_iterator equ_it = equation_number.begin(); for (UInt d = 0; d < nb_degree_of_freedoms; ++d, ++arr_it, ++equ_it) { global_array(*equ_it) += scale_factor * (*arr_it); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ DOFManagerDefault::DOFDataDefault::DOFDataDefault(const ID & dof_id) : DOFData(dof_id), associated_nodes(0, 1, dof_id + "associated_nodes") {} /* -------------------------------------------------------------------------- */ DOFManager::DOFData & DOFManagerDefault::getNewDOFData(const ID & dof_id) { this->dofs[dof_id] = std::make_unique(dof_id); return *this->dofs[dof_id]; } /* -------------------------------------------------------------------------- */ class GlobalDOFInfoDataAccessor : public DataAccessor { public: using size_type = typename std::unordered_map>::size_type; GlobalDOFInfoDataAccessor() = default; void addDOFToNode(UInt node, UInt dof) { dofs_per_node[node].push_back(dof); } UInt getNthDOFForNode(UInt nth_dof, UInt node) const { return dofs_per_node.find(node)->second[nth_dof]; } UInt getNbData(const Array & nodes, const SynchronizationTag & tag) const override { if (tag == _gst_size) { return nodes.size() * sizeof(size_type); } if (tag == _gst_update) { UInt total_size = 0; for (auto node : nodes) { auto it = dofs_per_node.find(node); if (it != dofs_per_node.end()) total_size += CommunicationBuffer::sizeInBuffer(it->second); } return total_size; } return 0; } void packData(CommunicationBuffer & buffer, const Array & nodes, const SynchronizationTag & tag) const override { for (auto node : nodes) { auto it = dofs_per_node.find(node); if (tag == _gst_size) { if (it != dofs_per_node.end()) { buffer << it->second.size(); } else { buffer << 0; } } else if (tag == _gst_update) { if (it != dofs_per_node.end()) buffer << it->second; } } } void unpackData(CommunicationBuffer & buffer, const Array & nodes, const SynchronizationTag & tag) override { for (auto node : nodes) { auto it = dofs_per_node.find(node); if (tag == _gst_size) { size_type size; buffer >> size; if (size != 0) dofs_per_node[node].resize(size); } else if (tag == _gst_update) { if (it != dofs_per_node.end()) buffer >> it->second; } } } protected: std::unordered_map> dofs_per_node; }; /* -------------------------------------------------------------------------- */ void DOFManagerDefault::registerDOFsInternal(const ID & dof_id, UInt nb_dofs, UInt nb_pure_local_dofs) { // auto prank = this->communicator.whoAmI(); // auto psize = this->communicator.getNbProc(); // access the relevant data to update auto & dof_data = this->getDOFDataTyped(dof_id); const auto & support_type = dof_data.support_type; const auto & group = dof_data.group_support; if (support_type == _dst_nodal and group != "__mesh__") { auto & support_nodes = this->mesh->getElementGroup(group).getNodeGroup().getNodes(); this->updateDOFsData( dof_data, nb_dofs, nb_pure_local_dofs, [&support_nodes](UInt node) -> UInt { return support_nodes[node]; }); } else { this->updateDOFsData(dof_data, nb_dofs, nb_pure_local_dofs, [](UInt node) -> UInt { return node; }); } // update the synchronizer if needed if (this->synchronizer) this->synchronizer->registerDOFs(dof_id); } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::registerDOFs(const ID & dof_id, Array & dofs_array, const DOFSupportType & support_type) { // stores the current numbers of dofs UInt nb_dofs_old = this->local_system_size; UInt nb_pure_local_dofs_old = this->pure_local_system_size; // update or create the dof_data DOFManager::registerDOFs(dof_id, dofs_array, support_type); UInt nb_dofs = this->local_system_size - nb_dofs_old; UInt nb_pure_local_dofs = this->pure_local_system_size - nb_pure_local_dofs_old; this->registerDOFsInternal(dof_id, nb_dofs, nb_pure_local_dofs); } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::registerDOFs(const ID & dof_id, Array & dofs_array, const ID & group_support) { // stores the current numbers of dofs UInt nb_dofs_old = this->local_system_size; UInt nb_pure_local_dofs_old = this->pure_local_system_size; // update or create the dof_data DOFManager::registerDOFs(dof_id, dofs_array, group_support); UInt nb_dofs = this->local_system_size - nb_dofs_old; UInt nb_pure_local_dofs = this->pure_local_system_size - nb_pure_local_dofs_old; this->registerDOFsInternal(dof_id, nb_dofs, nb_pure_local_dofs); } /* -------------------------------------------------------------------------- */ SparseMatrix & DOFManagerDefault::getNewMatrix(const ID & id, const MatrixType & matrix_type) { ID matrix_id = this->id + ":mtx:" + id; std::unique_ptr sm = std::make_unique(*this, matrix_type, matrix_id); return this->registerSparseMatrix(matrix_id, sm); } /* -------------------------------------------------------------------------- */ SparseMatrix & DOFManagerDefault::getNewMatrix(const ID & id, const ID & matrix_to_copy_id) { ID matrix_id = this->id + ":mtx:" + id; SparseMatrixAIJ & sm_to_copy = this->getMatrix(matrix_to_copy_id); std::unique_ptr sm = std::make_unique(sm_to_copy, matrix_id); return this->registerSparseMatrix(matrix_id, sm); } /* -------------------------------------------------------------------------- */ SparseMatrixAIJ & DOFManagerDefault::getMatrix(const ID & id) { SparseMatrix & matrix = DOFManager::getMatrix(id); return dynamic_cast(matrix); } /* -------------------------------------------------------------------------- */ NonLinearSolver & DOFManagerDefault::getNewNonLinearSolver(const ID & id, const NonLinearSolverType & type) { ID non_linear_solver_id = this->id + ":nls:" + id; std::unique_ptr nls; switch (type) { #if defined(AKANTU_IMPLICIT) case _nls_newton_raphson: case _nls_newton_raphson_modified: { nls = std::make_unique( *this, type, non_linear_solver_id, this->memory_id); break; } case _nls_linear: { nls = std::make_unique( *this, type, non_linear_solver_id, this->memory_id); break; } #endif case _nls_lumped: { nls = std::make_unique( *this, type, non_linear_solver_id, this->memory_id); break; } default: AKANTU_EXCEPTION("The asked type of non linear solver is not supported by " "this dof manager"); } return this->registerNonLinearSolver(non_linear_solver_id, nls); } /* -------------------------------------------------------------------------- */ TimeStepSolver & DOFManagerDefault::getNewTimeStepSolver(const ID & id, const TimeStepSolverType & type, NonLinearSolver & non_linear_solver) { ID time_step_solver_id = this->id + ":tss:" + id; std::unique_ptr tss = std::make_unique( *this, type, non_linear_solver, time_step_solver_id, this->memory_id); return this->registerTimeStepSolver(time_step_solver_id, tss); } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::clearResidual() { this->residual.resize(this->local_system_size); this->residual.clear(); } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::clearMatrix(const ID & mtx) { this->getMatrix(mtx).clear(); } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::clearLumpedMatrix(const ID & mtx) { this->getLumpedMatrix(mtx).clear(); } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::updateGlobalBlockedDofs() { auto it = this->dofs.begin(); auto end = this->dofs.end(); this->previous_global_blocked_dofs.copy(this->global_blocked_dofs); this->global_blocked_dofs.resize(this->local_system_size); this->global_blocked_dofs.clear(); for (; it != end; ++it) { if (!this->hasBlockedDOFs(it->first)) continue; DOFData & dof_data = *it->second; this->assembleToGlobalArray(it->first, *dof_data.blocked_dofs, this->global_blocked_dofs, true); } } /* -------------------------------------------------------------------------- */ template void DOFManagerDefault::getArrayPerDOFs(const ID & dof_id, const Array & global_array, Array & local_array) const { AKANTU_DEBUG_IN(); const Array & equation_number = this->getLocalEquationNumbers(dof_id); UInt nb_degree_of_freedoms = equation_number.size(); local_array.resize(nb_degree_of_freedoms / local_array.getNbComponent()); auto loc_it = local_array.begin_reinterpret(nb_degree_of_freedoms); auto equ_it = equation_number.begin(); for (UInt d = 0; d < nb_degree_of_freedoms; ++d, ++loc_it, ++equ_it) { (*loc_it) = global_array(*equ_it); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ // void DOFManagerDefault::getEquationsNumbers(const ID & dof_id, // Array & equation_numbers) { // AKANTU_DEBUG_IN(); // this->getArrayPerDOFs(dof_id, this->global_equation_number, // equation_numbers); AKANTU_DEBUG_OUT(); // } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::getSolutionPerDOFs(const ID & dof_id, Array & solution_array) { AKANTU_DEBUG_IN(); this->getArrayPerDOFs(dof_id, this->global_solution, solution_array); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::getLumpedMatrixPerDOFs(const ID & dof_id, const ID & lumped_mtx, Array & lumped) { AKANTU_DEBUG_IN(); this->getArrayPerDOFs(dof_id, this->getLumpedMatrix(lumped_mtx), lumped); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::assembleToResidual( const ID & dof_id, const Array & array_to_assemble, Real scale_factor) { AKANTU_DEBUG_IN(); this->assembleToGlobalArray(dof_id, array_to_assemble, this->residual, scale_factor); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::assembleToLumpedMatrix( const ID & dof_id, const Array & array_to_assemble, const ID & lumped_mtx, Real scale_factor) { AKANTU_DEBUG_IN(); Array & lumped = this->getLumpedMatrix(lumped_mtx); this->assembleToGlobalArray(dof_id, array_to_assemble, lumped, scale_factor); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::assembleMatMulVectToResidual(const ID & dof_id, const ID & A_id, const Array & x, Real scale_factor) { SparseMatrixAIJ & A = this->getMatrix(A_id); // Array data_cache(this->local_system_size, 1, 0.); this->data_cache.resize(this->local_system_size); this->data_cache.clear(); this->assembleToGlobalArray(dof_id, x, data_cache, 1.); A.matVecMul(data_cache, this->residual, scale_factor, 1.); } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::assembleLumpedMatMulVectToResidual( const ID & dof_id, const ID & A_id, const Array & x, Real scale_factor) { const Array & A = this->getLumpedMatrix(A_id); // Array data_cache(this->local_system_size, 1, 0.); this->data_cache.resize(this->local_system_size); this->data_cache.clear(); this->assembleToGlobalArray(dof_id, x, data_cache, scale_factor); Array::const_scalar_iterator A_it = A.begin(); Array::const_scalar_iterator A_end = A.end(); Array::const_scalar_iterator x_it = data_cache.begin(); Array::scalar_iterator r_it = this->residual.begin(); for (; A_it != A_end; ++A_it, ++x_it, ++r_it) { *r_it += *A_it * *x_it; } } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::assembleElementalMatricesToMatrix( const ID & matrix_id, const ID & dof_id, const Array & elementary_mat, const ElementType & type, const GhostType & ghost_type, const MatrixType & elemental_matrix_type, const Array & filter_elements) { AKANTU_DEBUG_IN(); DOFData & dof_data = this->getDOFData(dof_id); AKANTU_DEBUG_ASSERT(dof_data.support_type == _dst_nodal, "This function applies only on Nodal dofs"); this->addToProfile(matrix_id, dof_id, type, ghost_type); const Array & equation_number = this->getLocalEquationNumbers(dof_id); SparseMatrixAIJ & A = this->getMatrix(matrix_id); UInt nb_element; UInt * filter_it = nullptr; if (filter_elements != empty_filter) { nb_element = filter_elements.size(); filter_it = filter_elements.storage(); } else { if (dof_data.group_support != "__mesh__") { const Array & group_elements = this->mesh->getElementGroup(dof_data.group_support) .getElements(type, ghost_type); nb_element = group_elements.size(); filter_it = group_elements.storage(); } else { nb_element = this->mesh->getNbElement(type, ghost_type); } } AKANTU_DEBUG_ASSERT(elementary_mat.size() == nb_element, "The vector elementary_mat(" << elementary_mat.getID() << ") has not the good size."); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt nb_degree_of_freedom = dof_data.dof->getNbComponent(); const Array & connectivity = this->mesh->getConnectivity(type, ghost_type); auto conn_begin = connectivity.begin(nb_nodes_per_element); auto conn_it = conn_begin; UInt size_mat = nb_nodes_per_element * nb_degree_of_freedom; Vector element_eq_nb(nb_degree_of_freedom * nb_nodes_per_element); Array::const_matrix_iterator el_mat_it = elementary_mat.begin(size_mat, size_mat); for (UInt e = 0; e < nb_element; ++e, ++el_mat_it) { if (filter_it) conn_it = conn_begin + *filter_it; this->extractElementEquationNumber(equation_number, *conn_it, nb_degree_of_freedom, element_eq_nb); std::transform(element_eq_nb.storage(), element_eq_nb.storage() + element_eq_nb.size(), element_eq_nb.storage(), [&](UInt & local) -> UInt { return this->localToGlobalEquationNumber(local); }); if (filter_it) ++filter_it; else ++conn_it; if (A.getMatrixType() == _symmetric) if (elemental_matrix_type == _symmetric) - this->addSymmetricElementalMatrixToSymmetric( - A, *el_mat_it, element_eq_nb, A.size()); + this->addSymmetricElementalMatrixToSymmetric(A, *el_mat_it, + element_eq_nb, A.size()); else - this->addUnsymmetricElementalMatrixToSymmetric( - A, *el_mat_it, element_eq_nb, A.size()); + this->addUnsymmetricElementalMatrixToSymmetric(A, *el_mat_it, + element_eq_nb, A.size()); else this->addElementalMatrixToUnsymmetric(A, *el_mat_it, element_eq_nb, A.size()); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::assemblePreassembledMatrix( const ID & dof_id_m, const ID & dof_id_n, const ID & matrix_id, const TermsToAssemble & terms) { const Array & equation_number_m = this->getLocalEquationNumbers(dof_id_m); const Array & equation_number_n = this->getLocalEquationNumbers(dof_id_n); SparseMatrixAIJ & A = this->getMatrix(matrix_id); for (const auto & term : terms) { UInt gi = this->localToGlobalEquationNumber(equation_number_m(term.i())); UInt gj = this->localToGlobalEquationNumber(equation_number_n(term.j())); A.add(gi, gj, term); } } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::addToProfile(const ID & matrix_id, const ID & dof_id, const ElementType & type, const GhostType & ghost_type) { AKANTU_DEBUG_IN(); const DOFData & dof_data = this->getDOFData(dof_id); if (dof_data.support_type != _dst_nodal) return; auto mat_dof = std::make_pair(matrix_id, dof_id); auto type_pair = std::make_pair(type, ghost_type); auto prof_it = this->matrix_profiled_dofs.find(mat_dof); if (prof_it != this->matrix_profiled_dofs.end() && std::find(prof_it->second.begin(), prof_it->second.end(), type_pair) != prof_it->second.end()) return; UInt nb_degree_of_freedom_per_node = dof_data.dof->getNbComponent(); const Array & equation_number = this->getLocalEquationNumbers(dof_id); SparseMatrixAIJ & A = this->getMatrix(matrix_id); UInt size = A.size(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); const auto & connectivity = this->mesh->getConnectivity(type, ghost_type); auto cbegin = connectivity.begin(nb_nodes_per_element); auto cit = cbegin; UInt nb_elements = connectivity.size(); UInt * ge_it = nullptr; if (dof_data.group_support != "__mesh__") { const Array & group_elements = this->mesh->getElementGroup(dof_data.group_support) .getElements(type, ghost_type); ge_it = group_elements.storage(); nb_elements = group_elements.size(); } UInt size_mat = nb_nodes_per_element * nb_degree_of_freedom_per_node; Vector element_eq_nb(size_mat); for (UInt e = 0; e < nb_elements; ++e) { if (ge_it) cit = cbegin + *ge_it; this->extractElementEquationNumber( equation_number, *cit, nb_degree_of_freedom_per_node, element_eq_nb); std::transform(element_eq_nb.storage(), element_eq_nb.storage() + element_eq_nb.size(), element_eq_nb.storage(), [&](UInt & local) -> UInt { return this->localToGlobalEquationNumber(local); }); if (ge_it) ++ge_it; else ++cit; for (UInt i = 0; i < size_mat; ++i) { UInt c_irn = element_eq_nb(i); if (c_irn < size) { for (UInt j = 0; j < size_mat; ++j) { UInt c_jcn = element_eq_nb(j); if (c_jcn < size) { A.add(c_irn, c_jcn); } } } } } this->matrix_profiled_dofs[mat_dof].push_back(type_pair); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::applyBoundary(const ID & matrix_id) { SparseMatrixAIJ & J = this->getMatrix(matrix_id); if (this->jacobian_release == J.getValueRelease()) { auto are_equal = std::equal(global_blocked_dofs.begin(), global_blocked_dofs.end(), previous_global_blocked_dofs.begin()); if (are_equal) J.applyBoundary(); previous_global_blocked_dofs.copy(global_blocked_dofs); } else { J.applyBoundary(); } this->jacobian_release = J.getValueRelease(); } /* -------------------------------------------------------------------------- */ const Array & DOFManagerDefault::getGlobalResidual() { if (this->synchronizer) { if (not this->global_residual) { this->global_residual = std::make_unique>(0, 1, "global_residual"); } if (this->synchronizer->getCommunicator().whoAmI() == 0) { this->global_residual->resize(this->system_size); this->synchronizer->gather(this->residual, *this->global_residual); } else { this->synchronizer->gather(this->residual); } return *this->global_residual; } else { return this->residual; } } /* -------------------------------------------------------------------------- */ const Array & DOFManagerDefault::getResidual() const { return this->residual; } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::setGlobalSolution(const Array & solution) { if (this->synchronizer) { if (this->synchronizer->getCommunicator().whoAmI() == 0) { this->synchronizer->scatter(this->global_solution, solution); } else { this->synchronizer->scatter(this->global_solution); } } else { AKANTU_DEBUG_ASSERT(solution.size() == this->global_solution.size(), "Sequential call to this function needs the solution " "to be the same size as the global_solution"); this->global_solution.copy(solution); } } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::onNodesAdded(const Array & nodes_list, const NewNodesEvent & event) { DOFManager::onNodesAdded(nodes_list, event); if (this->synchronizer) this->synchronizer->onNodesAdded(nodes_list); } /* -------------------------------------------------------------------------- */ std::pair DOFManagerDefault::updateNodalDOFs(const ID & dof_id, const Array & nodes_list) { UInt nb_new_local_dofs, nb_new_pure_local; std::tie(nb_new_local_dofs, nb_new_pure_local) = DOFManager::updateNodalDOFs(dof_id, nodes_list); auto & dof_data = this->getDOFDataTyped(dof_id); updateDOFsData(dof_data, nb_new_local_dofs, nb_new_pure_local, [&nodes_list](UInt pos) -> UInt { return nodes_list[pos]; }); return std::make_pair(nb_new_local_dofs, nb_new_pure_local); } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::updateDOFsData( DOFDataDefault & dof_data, UInt nb_new_local_dofs, UInt nb_new_pure_local, const std::function & getNode) { auto prank = this->communicator.whoAmI(); auto psize = this->communicator.getNbProc(); // access the relevant data to update const auto & support_type = dof_data.support_type; GlobalDOFInfoDataAccessor data_accessor; // resize all relevant arrays this->residual.resize(this->local_system_size); this->dofs_type.resize(this->local_system_size); this->global_solution.resize(this->local_system_size); this->global_blocked_dofs.resize(this->local_system_size); this->previous_global_blocked_dofs.resize(this->local_system_size); this->global_equation_number.resize(this->local_system_size); for (auto & lumped_matrix : lumped_matrices) lumped_matrix.second->resize(this->local_system_size); - dof_data.local_equation_number.reserve( - dof_data.local_equation_number.size() + nb_new_local_dofs); + dof_data.local_equation_number.reserve(dof_data.local_equation_number.size() + + nb_new_local_dofs); // determine the first local/global dof id to use Array nb_dofs_per_proc(psize); nb_dofs_per_proc(prank) = nb_new_pure_local; this->communicator.allGather(nb_dofs_per_proc); this->first_global_dof_id += std::accumulate( nb_dofs_per_proc.begin(), nb_dofs_per_proc.begin() + prank, 0); UInt first_dof_id = this->local_system_size - nb_new_local_dofs; if (support_type == _dst_nodal) { dof_data.associated_nodes.reserve(dof_data.associated_nodes.size() + nb_new_local_dofs); } // update per dof info for (UInt d = 0; d < nb_new_local_dofs; ++d) { UInt local_eq_num = first_dof_id + d; dof_data.local_equation_number.push_back(local_eq_num); bool is_local_dof = true; // determine the dof type switch (support_type) { case _dst_nodal: { UInt node = getNode(d / dof_data.dof->getNbComponent()); this->dofs_type(local_eq_num) = this->mesh->getNodeType(node); dof_data.associated_nodes.push_back(node); is_local_dof = this->mesh->isLocalOrMasterNode(node); if (is_local_dof) { data_accessor.addDOFToNode(node, first_global_dof_id); } break; } case _dst_generic: { this->dofs_type(local_eq_num) = _nt_normal; break; } default: { AKANTU_EXCEPTION("This type of dofs is not handled yet."); } } // update global id for local dofs if (is_local_dof) { this->global_equation_number(local_eq_num) = this->first_global_dof_id; this->global_to_local_mapping[this->first_global_dof_id] = local_eq_num; ++this->first_global_dof_id; } else { this->global_equation_number(local_eq_num) = 0; } } // synchronize the global numbering for slaves if (support_type == _dst_nodal && this->synchronizer) { auto nb_dofs_per_node = dof_data.dof->getNbComponent(); auto & node_synchronizer = this->mesh->getNodeSynchronizer(); node_synchronizer.synchronizeOnce(data_accessor, _gst_size); node_synchronizer.synchronizeOnce(data_accessor, _gst_update); std::vector counters(nb_new_local_dofs / nb_dofs_per_node); for (UInt d = 0; d < nb_new_local_dofs; ++d) { UInt local_eq_num = first_dof_id + d; if (not this->isSlaveDOF(local_eq_num)) continue; UInt node = d / nb_dofs_per_node; UInt dof_count = counters[node]++; node = getNode(node); UInt global_dof_id = data_accessor.getNthDOFForNode(dof_count, node); this->global_equation_number(local_eq_num) = global_dof_id; this->global_to_local_mapping[global_dof_id] = local_eq_num; } } } /* -------------------------------------------------------------------------- */ // register in factory -static bool dof_manager_is_registered = DefaultDOFManagerFactory::getInstance().registerAllocator( - "default", [](const ID & id, const MemoryID & mem_id) -> std::unique_ptr { - return std::make_unique(id + ":dof_manager_default:", mem_id); - }); +static bool default_dof_manager_is_registered = + DefaultDOFManagerFactory::getInstance().registerAllocator( + "default", [](const ID & id, + const MemoryID & mem_id) -> std::unique_ptr { + return std::make_unique(id, mem_id); + }); + +static bool dof_manager_is_registered = + DOFManagerFactory::getInstance().registerAllocator( + "default", [](Mesh & mesh, const ID & id, + const MemoryID & mem_id) -> std::unique_ptr { + return std::make_unique(mesh, id, mem_id); + }); } // namespace akantu diff --git a/src/model/model_solver.cc b/src/model/model_solver.cc index 89441b74c..5b39127fe 100644 --- a/src/model/model_solver.cc +++ b/src/model/model_solver.cc @@ -1,364 +1,366 @@ /** * @file model_solver.cc * * @author Nicolas Richart * * @date Wed Aug 12 13:31:56 2015 * * @brief Implementation of ModelSolver * * @section LICENSE * * Copyright (©) 2010-2011 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 "model_solver.hh" #include "dof_manager.hh" #include "dof_manager_default.hh" #include "mesh.hh" #include "non_linear_solver.hh" #include "time_step_solver.hh" #if defined(AKANTU_USE_PETSC) #include "dof_manager_petsc.hh" #endif /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ template static T getOptionToType(const std::string & opt_str) { std::stringstream sstr(opt_str); T opt; sstr >> opt; return opt; } /* -------------------------------------------------------------------------- */ ModelSolver::ModelSolver(Mesh & mesh, const ModelType & type, const ID & id, UInt memory_id) : Parsable(ParserType::_model, id), SolverCallback(), model_type(type), parent_id(id), parent_memory_id(memory_id), mesh(mesh), dof_manager(nullptr), default_solver_id("") {} /* -------------------------------------------------------------------------- */ ModelSolver::~ModelSolver() = default; /* -------------------------------------------------------------------------- */ std::tuple ModelSolver::getParserSection() { - auto - sub_sections = getStaticParser().getSubSections(ParserType::_model); + auto sub_sections = getStaticParser().getSubSections(ParserType::_model); - auto it = std::find_if(sub_sections.begin(), sub_sections.end(), [&](auto && section) { - ModelType type = getOptionToType(section.getName()); - // default id should be the model type if not defined - std::string name = section.getParameter("name", this->parent_id); - return type == model_type and name == this->parent_id; - }); + auto it = std::find_if( + sub_sections.begin(), sub_sections.end(), [&](auto && section) { + ModelType type = getOptionToType(section.getName()); + // default id should be the model type if not defined + std::string name = section.getParameter("name", this->parent_id); + return type == model_type and name == this->parent_id; + }); - if(it == sub_sections.end()) + if (it == sub_sections.end()) return {ParserSection(), true}; return {*it, false}; } - /* -------------------------------------------------------------------------- */ void ModelSolver::initDOFManager() { // default without external solver activated at compilation same as mumps that // is the historical solver but with only the lumped solver ID solver_type = "default"; #if defined(AKANTU_USE_MUMPS) solver_type = "default"; #elif defined(AKANTU_USE_PETSC) solver_type = "petsc"; #endif ParserSection section; bool is_empty; std::tie(section, is_empty) = this->getParserSection(); if (not is_empty) { solver_type = section.getOption(solver_type); this->initDOFManager(section, solver_type); } else { this->initDOFManager(solver_type); } } /* -------------------------------------------------------------------------- */ void ModelSolver::initDOFManager(const ID & solver_type) { try { - this->dof_manager = DefaultDOFManagerFactory::getInstance().allocate(solver_type, - this->parent_id + ":dof_manager" + solver_type, this->parent_memory_id); + this->dof_manager = DOFManagerFactory::getInstance().allocate( + solver_type, mesh, this->parent_id + ":dof_manager" + solver_type, + this->parent_memory_id); } catch (...) { AKANTU_EXCEPTION( "To use the solver " << solver_type << " you will have to code it. This is an unknown solver type."); } this->setDOFManager(*this->dof_manager); } /* -------------------------------------------------------------------------- */ void ModelSolver::initDOFManager(const ParserSection & section, const ID & solver_type) { this->initDOFManager(solver_type); auto sub_sections = section.getSubSections(ParserType::_time_step_solver); // parsing the time step solvers for (auto && section : sub_sections) { ID type = section.getName(); ID solver_id = section.getParameter("name", type); auto tss_type = getOptionToType(type); auto tss_options = this->getDefaultSolverOptions(tss_type); - auto sub_solvers_sect = section.getSubSections(ParserType::_non_linear_solver); + auto sub_solvers_sect = + section.getSubSections(ParserType::_non_linear_solver); auto nb_non_linear_solver_section = section.getNbSubSections(ParserType::_non_linear_solver); auto nls_type = tss_options.non_linear_solver_type; if (nb_non_linear_solver_section == 1) { auto && nls_section = *(sub_solvers_sect.first); nls_type = getOptionToType(nls_section.getName()); } else if (nb_non_linear_solver_section > 0) { AKANTU_EXCEPTION("More than one non linear solver are provided for the " "time step solver " << solver_id); } this->getNewSolver(solver_id, tss_type, nls_type); if (nb_non_linear_solver_section == 1) { const auto & nls_section = *(sub_solvers_sect.first); this->dof_manager->getNonLinearSolver(solver_id).parseSection( nls_section); } auto sub_integrator_sections = section.getSubSections(ParserType::_integration_scheme); for (auto && is_section : sub_integrator_sections) { const auto & dof_type_str = is_section.getName(); ID dof_id; try { ID tmp = is_section.getParameter("name"); dof_id = tmp; } catch (...) { AKANTU_EXCEPTION("No degree of freedom name specified for the " - "integration scheme of type " << dof_type_str); + "integration scheme of type " + << dof_type_str); } - auto it_type = - getOptionToType(dof_type_str); + auto it_type = getOptionToType(dof_type_str); IntegrationScheme::SolutionType s_type = is_section.getParameter( "solution_type", tss_options.solution_type[dof_id]); this->setIntegrationScheme(solver_id, dof_id, it_type, s_type); } for (auto & is_type : tss_options.integration_scheme_type) { if (!this->hasIntegrationScheme(solver_id, is_type.first)) { this->setIntegrationScheme(solver_id, is_type.first, is_type.second, tss_options.solution_type[is_type.first]); } } } if (section.hasParameter("default_solver")) { ID default_solver = section.getParameter("default_solver"); if (this->hasSolver(default_solver)) { this->setDefaultSolver(default_solver); } else AKANTU_EXCEPTION( - "The solver \"" << default_solver - <<"\" was not created, it cannot be set as default solver"); + "The solver \"" + << default_solver + << "\" was not created, it cannot be set as default solver"); } } /* -------------------------------------------------------------------------- */ TimeStepSolver & ModelSolver::getSolver(const ID & solver_id) { ID tmp_solver_id = solver_id; if (tmp_solver_id == "") tmp_solver_id = this->default_solver_id; TimeStepSolver & tss = this->dof_manager->getTimeStepSolver(tmp_solver_id); return tss; } /* -------------------------------------------------------------------------- */ const TimeStepSolver & ModelSolver::getSolver(const ID & solver_id) const { ID tmp_solver_id = solver_id; if (solver_id == "") tmp_solver_id = this->default_solver_id; const TimeStepSolver & tss = this->dof_manager->getTimeStepSolver(tmp_solver_id); return tss; } /* -------------------------------------------------------------------------- */ TimeStepSolver & ModelSolver::getTimeStepSolver(const ID & solver_id) { return this->getSolver(solver_id); } /* -------------------------------------------------------------------------- */ const TimeStepSolver & ModelSolver::getTimeStepSolver(const ID & solver_id) const { return this->getSolver(solver_id); } /* -------------------------------------------------------------------------- */ NonLinearSolver & ModelSolver::getNonLinearSolver(const ID & solver_id) { return this->getSolver(solver_id).getNonLinearSolver(); } /* -------------------------------------------------------------------------- */ const NonLinearSolver & ModelSolver::getNonLinearSolver(const ID & solver_id) const { return this->getSolver(solver_id).getNonLinearSolver(); } /* -------------------------------------------------------------------------- */ bool ModelSolver::hasSolver(const ID & solver_id) const { ID tmp_solver_id = solver_id; if (solver_id == "") tmp_solver_id = this->default_solver_id; if (not this->dof_manager) { AKANTU_EXCEPTION("No DOF manager was initialized"); } return this->dof_manager->hasTimeStepSolver(tmp_solver_id); } /* -------------------------------------------------------------------------- */ void ModelSolver::setDefaultSolver(const ID & solver_id) { AKANTU_DEBUG_ASSERT( this->hasSolver(solver_id), "Cannot set the default solver to a solver that does not exists"); this->default_solver_id = solver_id; } /* -------------------------------------------------------------------------- */ void ModelSolver::solveStep(const ID & solver_id) { AKANTU_DEBUG_IN(); TimeStepSolver & tss = this->getSolver(solver_id); // make one non linear solve tss.solveStep(*this); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ModelSolver::getNewSolver(const ID & solver_id, TimeStepSolverType time_step_solver_type, NonLinearSolverType non_linear_solver_type) { if (this->default_solver_id == "") { this->default_solver_id = solver_id; } if (non_linear_solver_type == _nls_auto) { switch (time_step_solver_type) { case _tsst_dynamic: case _tsst_static: non_linear_solver_type = _nls_newton_raphson; break; case _tsst_dynamic_lumped: non_linear_solver_type = _nls_lumped; break; case _tsst_not_defined: AKANTU_EXCEPTION(time_step_solver_type << " is not a valid time step solver type"); break; } } this->initSolver(time_step_solver_type, non_linear_solver_type); NonLinearSolver & nls = this->dof_manager->getNewNonLinearSolver( solver_id, non_linear_solver_type); this->dof_manager->getNewTimeStepSolver(solver_id, time_step_solver_type, nls); } /* -------------------------------------------------------------------------- */ Real ModelSolver::getTimeStep(const ID & solver_id) const { const TimeStepSolver & tss = this->getSolver(solver_id); return tss.getTimeStep(); } /* -------------------------------------------------------------------------- */ void ModelSolver::setTimeStep(Real time_step, const ID & solver_id) { TimeStepSolver & tss = this->getSolver(solver_id); return tss.setTimeStep(time_step); } /* -------------------------------------------------------------------------- */ void ModelSolver::setIntegrationScheme( const ID & solver_id, const ID & dof_id, const IntegrationSchemeType & integration_scheme_type, IntegrationScheme::SolutionType solution_type) { TimeStepSolver & tss = this->dof_manager->getTimeStepSolver(solver_id); tss.setIntegrationScheme(dof_id, integration_scheme_type, solution_type); } /* -------------------------------------------------------------------------- */ bool ModelSolver::hasDefaultSolver() const { return (this->default_solver_id != ""); } /* -------------------------------------------------------------------------- */ bool ModelSolver::hasIntegrationScheme(const ID & solver_id, const ID & dof_id) const { TimeStepSolver & tss = this->dof_manager->getTimeStepSolver(solver_id); return tss.hasIntegrationScheme(dof_id); } /* -------------------------------------------------------------------------- */ void ModelSolver::predictor() {} /* -------------------------------------------------------------------------- */ void ModelSolver::corrector() {} /* -------------------------------------------------------------------------- */ TimeStepSolverType ModelSolver::getDefaultSolverType() const { return _tsst_dynamic_lumped; } /* -------------------------------------------------------------------------- */ ModelSolverOptions ModelSolver::getDefaultSolverOptions(__attribute__((unused)) const TimeStepSolverType & type) const { ModelSolverOptions options; options.non_linear_solver_type = _nls_auto; return options; } } // namespace akantu diff --git a/src/model/solid_mechanics/solid_mechanics_model.cc b/src/model/solid_mechanics/solid_mechanics_model.cc index 36eeb4a94..2d769cc40 100644 --- a/src/model/solid_mechanics/solid_mechanics_model.cc +++ b/src/model/solid_mechanics/solid_mechanics_model.cc @@ -1,1202 +1,1204 @@ /** * @file solid_mechanics_model.cc * * @author Ramin Aghababaei * @author Guillaume Anciaux * @author Aurelia Isabel Cuba Ramos * @author David Simon Kammer * @author Daniel Pino Muñoz * @author Nicolas Richart * @author Clement Roux * @author Marco Vocialta * * @date creation: Tue Jul 27 2010 * @date last modification: Tue Jan 19 2016 * * @brief Implementation of the SolidMechanicsModel class * * @section LICENSE * * Copyright (©) 2010-2012, 2014, 2015 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.hh" #include "integrator_gauss.hh" #include "shape_lagrange.hh" #include "solid_mechanics_model_tmpl.hh" #include "communicator.hh" #include "element_synchronizer.hh" #include "sparse_matrix.hh" #include "synchronizer_registry.hh" #include "dumpable_inline_impl.hh" #ifdef AKANTU_USE_IOHELPER #include "dumper_paraview.hh" #endif #include "material_non_local.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ /** * A solid mechanics model need a mesh and a dimension to be created. the model * by it self can not do a lot, the good init functions should be called in * order to configure the model depending on what we want to do. * * @param mesh mesh representing the model we want to simulate * @param dim spatial dimension of the problem, if dim = 0 (default value) the * dimension of the problem is assumed to be the on of the mesh * @param id an id to identify the model */ SolidMechanicsModel::SolidMechanicsModel(Mesh & mesh, UInt dim, const ID & id, const MemoryID & memory_id) : Model(mesh, ModelType::_solid_mechanics_model, dim, id, memory_id), BoundaryCondition(), f_m2a(1.0), displacement(nullptr), previous_displacement(nullptr), displacement_increment(nullptr), mass(nullptr), velocity(nullptr), acceleration(nullptr), external_force(nullptr), internal_force(nullptr), blocked_dofs(nullptr), current_position(nullptr), material_index("material index", id, memory_id), material_local_numbering("material local numbering", id, memory_id), increment_flag(false), are_materials_instantiated(false) { //, pbc_synch(nullptr) { AKANTU_DEBUG_IN(); this->registerFEEngineObject("SolidMechanicsFEEngine", mesh, Model::spatial_dimension); #if defined(AKANTU_USE_IOHELPER) this->mesh.registerDumper("paraview_all", id, true); this->mesh.addDumpMesh(mesh, Model::spatial_dimension, _not_ghost, _ek_regular); #endif material_selector = std::make_shared(material_index), this->initDOFManager(); this->registerDataAccessor(*this); if (this->mesh.isDistributed()) { auto & synchronizer = this->mesh.getElementSynchronizer(); this->registerSynchronizer(synchronizer, _gst_material_id); this->registerSynchronizer(synchronizer, _gst_smm_mass); this->registerSynchronizer(synchronizer, _gst_smm_stress); this->registerSynchronizer(synchronizer, _gst_smm_boundary); this->registerSynchronizer(synchronizer, _gst_for_dump); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ SolidMechanicsModel::~SolidMechanicsModel() { AKANTU_DEBUG_IN(); for (auto & internal : this->registered_internals) { delete internal.second; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::setTimeStep(Real time_step, const ID & solver_id) { Model::setTimeStep(time_step, solver_id); #if defined(AKANTU_USE_IOHELPER) this->mesh.getDumper().setTimeStep(time_step); #endif } /* -------------------------------------------------------------------------- */ /* Initialization */ /* -------------------------------------------------------------------------- */ /** * This function groups many of the initialization in on function. For most of * basics case the function should be enough. The functions initialize the * model, the internal vectors, set them to 0, and depending on the parameters * it also initialize the explicit or implicit solver. * * @param material_file the file containing the materials to use * @param method the analysis method wanted. See the akantu::AnalysisMethod for * the different possibilities */ void SolidMechanicsModel::initFullImpl(const ModelOptions & options) { material_index.initialize(mesh, _element_kind = _ek_not_defined, _default_value = UInt(-1), _with_nb_element = true); material_local_numbering.initialize(mesh, _element_kind = _ek_not_defined, _with_nb_element = true); Model::initFullImpl(options); // initialize pbc if (this->pbc_pair.size() != 0) this->initPBC(); // initialize the materials if (this->parser.getLastParsedFile() != "") { this->instantiateMaterials(); } this->initMaterials(); this->initBC(*this, *displacement, *displacement_increment, *external_force); } /* -------------------------------------------------------------------------- */ TimeStepSolverType SolidMechanicsModel::getDefaultSolverType() const { return _tsst_dynamic_lumped; } /* -------------------------------------------------------------------------- */ ModelSolverOptions SolidMechanicsModel::getDefaultSolverOptions( const TimeStepSolverType & type) const { ModelSolverOptions options; switch (type) { case _tsst_dynamic_lumped: { options.non_linear_solver_type = _nls_lumped; options.integration_scheme_type["displacement"] = _ist_central_difference; options.solution_type["displacement"] = IntegrationScheme::_acceleration; break; } case _tsst_static: { options.non_linear_solver_type = _nls_newton_raphson; options.integration_scheme_type["displacement"] = _ist_pseudo_time; options.solution_type["displacement"] = IntegrationScheme::_not_defined; break; } case _tsst_dynamic: { if (this->method == _explicit_consistent_mass) { options.non_linear_solver_type = _nls_newton_raphson; options.integration_scheme_type["displacement"] = _ist_central_difference; options.solution_type["displacement"] = IntegrationScheme::_acceleration; } else { options.non_linear_solver_type = _nls_newton_raphson; options.integration_scheme_type["displacement"] = _ist_trapezoidal_rule_2; options.solution_type["displacement"] = IntegrationScheme::_displacement; } break; } default: AKANTU_EXCEPTION(type << " is not a valid time step solver type"); } return options; } /* -------------------------------------------------------------------------- */ std::tuple SolidMechanicsModel::getDefaultSolverID(const AnalysisMethod & method) { switch (method) { case _explicit_lumped_mass: { return std::make_tuple("explicit_lumped", _tsst_dynamic_lumped); } case _explicit_consistent_mass: { return std::make_tuple("explicit", _tsst_dynamic); } case _static: { return std::make_tuple("static", _tsst_static); } case _implicit_dynamic: { return std::make_tuple("implicit", _tsst_dynamic); } default: return std::make_tuple("unknown", _tsst_not_defined); } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::initSolver(TimeStepSolverType time_step_solver_type, NonLinearSolverType) { auto & dof_manager = this->getDOFManager(); /* ------------------------------------------------------------------------ */ // for alloc type of solvers this->allocNodalField(this->displacement, spatial_dimension, "displacement"); this->allocNodalField(this->previous_displacement, spatial_dimension, "previous_displacement"); this->allocNodalField(this->displacement_increment, spatial_dimension, "displacement_increment"); this->allocNodalField(this->internal_force, spatial_dimension, "internal_force"); this->allocNodalField(this->external_force, spatial_dimension, "external_force"); this->allocNodalField(this->blocked_dofs, spatial_dimension, "blocked_dofs"); this->allocNodalField(this->current_position, spatial_dimension, "current_position"); // initialize the current positions this->current_position->copy(this->mesh.getNodes()); /* ------------------------------------------------------------------------ */ if (!dof_manager.hasDOFs("displacement")) { dof_manager.registerDOFs("displacement", *this->displacement, _dst_nodal); dof_manager.registerBlockedDOFs("displacement", *this->blocked_dofs); dof_manager.registerDOFsIncrement("displacement", *this->displacement_increment); dof_manager.registerDOFsPrevious("displacement", *this->previous_displacement); } /* ------------------------------------------------------------------------ */ // for dynamic if (time_step_solver_type == _tsst_dynamic || time_step_solver_type == _tsst_dynamic_lumped) { this->allocNodalField(this->velocity, spatial_dimension, "velocity"); this->allocNodalField(this->acceleration, spatial_dimension, "acceleration"); if (!dof_manager.hasDOFsDerivatives("displacement", 1)) { dof_manager.registerDOFsDerivative("displacement", 1, *this->velocity); dof_manager.registerDOFsDerivative("displacement", 2, *this->acceleration); } } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::initFEEngineBoundary() { FEEngine & fem_boundary = getFEEngineBoundary(); fem_boundary.initShapeFunctions(_not_ghost); fem_boundary.initShapeFunctions(_ghost); fem_boundary.computeNormalsOnIntegrationPoints(_not_ghost); fem_boundary.computeNormalsOnIntegrationPoints(_ghost); } /* -------------------------------------------------------------------------- */ /** * Initialize the model,basically it pre-compute the shapes, shapes derivatives * and jacobian */ void SolidMechanicsModel::initModel() { /// \todo add the current position as a parameter to initShapeFunctions for /// large deformation getFEEngine().initShapeFunctions(_not_ghost); getFEEngine().initShapeFunctions(_ghost); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::assembleResidual() { AKANTU_DEBUG_IN(); /* ------------------------------------------------------------------------ */ // computes the internal forces this->assembleInternalForces(); /* ------------------------------------------------------------------------ */ this->getDOFManager().assembleToResidual("displacement", *this->external_force, 1); this->getDOFManager().assembleToResidual("displacement", *this->internal_force, 1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ MatrixType SolidMechanicsModel::getMatrixType(const ID & matrix_id) { // \TODO check the materials to know what is the correct answer if (matrix_id == "C") return _mt_not_defined; return _symmetric; } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::assembleMatrix(const ID & matrix_id) { if (matrix_id == "K") { this->assembleStiffnessMatrix(); } else if (matrix_id == "M") { this->assembleMass(); } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::assembleLumpedMatrix(const ID & matrix_id) { if (matrix_id == "M") { this->assembleMassLumped(); } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::beforeSolveStep() { for (auto & material : materials) material->beforeSolveStep(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::afterSolveStep() { for (auto & material : materials) material->afterSolveStep(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::predictor() { ++displacement_release; } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::corrector() { ++displacement_release; } /* -------------------------------------------------------------------------- */ /** * This function computes the internal forces as F_{int} = \int_{\Omega} N * \sigma d\Omega@f$ */ void SolidMechanicsModel::assembleInternalForces() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Assemble the internal forces"); this->internal_force->clear(); // compute the stresses of local elements AKANTU_DEBUG_INFO("Compute local stresses"); for (auto & material : materials) { material->computeAllStresses(_not_ghost); } /* ------------------------------------------------------------------------ */ /* Computation of the non local part */ if (this->non_local_manager) this->non_local_manager->computeAllNonLocalStresses(); // communicate the stresses AKANTU_DEBUG_INFO("Send data for residual assembly"); this->asynchronousSynchronize(_gst_smm_stress); // assemble the forces due to local stresses AKANTU_DEBUG_INFO("Assemble residual for local elements"); for (auto & material : materials) { material->assembleInternalForces(_not_ghost); } // finalize communications AKANTU_DEBUG_INFO("Wait distant stresses"); this->waitEndSynchronize(_gst_smm_stress); // assemble the stresses due to ghost elements AKANTU_DEBUG_INFO("Assemble residual for ghost elements"); for (auto & material : materials) { material->assembleInternalForces(_ghost); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::assembleStiffnessMatrix() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Assemble the new stiffness matrix."); // Check if materials need to recompute the matrix bool need_to_reassemble = false; for (auto & material : materials) { need_to_reassemble |= material->hasStiffnessMatrixChanged(); } if (need_to_reassemble) { this->getDOFManager().getMatrix("K").clear(); // call compute stiffness matrix on each local elements for (auto & material : materials) { material->assembleStiffnessMatrix(_not_ghost); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::updateCurrentPosition() { if (this->current_position_release == this->displacement_release) return; this->current_position->copy(this->mesh.getNodes()); auto cpos_it = this->current_position->begin(Model::spatial_dimension); auto cpos_end = this->current_position->end(Model::spatial_dimension); auto disp_it = this->displacement->begin(Model::spatial_dimension); for (; cpos_it != cpos_end; ++cpos_it, ++disp_it) { *cpos_it += *disp_it; } this->current_position_release = this->displacement_release; } /* -------------------------------------------------------------------------- */ const Array & SolidMechanicsModel::getCurrentPosition() { this->updateCurrentPosition(); return *this->current_position; } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::updateDataForNonLocalCriterion( ElementTypeMapReal & criterion) { const ID field_name = criterion.getName(); for (auto & material : materials) { if (!material->isInternal(field_name, _ek_regular)) continue; for (auto ghost_type : ghost_types) { material->flattenInternal(field_name, criterion, ghost_type, _ek_regular); } } } /* -------------------------------------------------------------------------- */ /* Information */ /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getStableTimeStep() { AKANTU_DEBUG_IN(); Real min_dt = getStableTimeStep(_not_ghost); /// reduction min over all processors mesh.getCommunicator().allReduce(min_dt, SynchronizerOperation::_min); AKANTU_DEBUG_OUT(); return min_dt; } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getStableTimeStep(const GhostType & ghost_type) { AKANTU_DEBUG_IN(); Real min_dt = std::numeric_limits::max(); this->updateCurrentPosition(); Element elem; elem.ghost_type = ghost_type; for (auto type : mesh.elementTypes(Model::spatial_dimension, ghost_type, _ek_regular)) { elem.type = type; UInt nb_nodes_per_element = mesh.getNbNodesPerElement(type); UInt nb_element = mesh.getNbElement(type); auto mat_indexes = material_index(type, ghost_type).begin(); auto mat_loc_num = material_local_numbering(type, ghost_type).begin(); Array X(0, nb_nodes_per_element * Model::spatial_dimension); FEEngine::extractNodalToElementField(mesh, *current_position, X, type, _not_ghost); auto X_el = X.begin(Model::spatial_dimension, nb_nodes_per_element); for (UInt el = 0; el < nb_element; ++el, ++X_el, ++mat_indexes, ++mat_loc_num) { elem.element = *mat_loc_num; Real el_h = getFEEngine().getElementInradius(*X_el, type); Real el_c = this->materials[*mat_indexes]->getCelerity(elem); Real el_dt = el_h / el_c; min_dt = std::min(min_dt, el_dt); } } AKANTU_DEBUG_OUT(); return min_dt; } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getKineticEnergy() { AKANTU_DEBUG_IN(); Real ekin = 0.; UInt nb_nodes = mesh.getNbNodes(); if (this->getDOFManager().hasLumpedMatrix("M")) { auto m_it = this->mass->begin(Model::spatial_dimension); auto m_end = this->mass->end(Model::spatial_dimension); auto v_it = this->velocity->begin(Model::spatial_dimension); for (UInt n = 0; m_it != m_end; ++n, ++m_it, ++v_it) { const Vector & v = *v_it; const Vector & m = *m_it; Real mv2 = 0; bool is_local_node = mesh.isLocalOrMasterNode(n); // bool is_not_pbc_slave_node = !isPBCSlaveNode(n); bool count_node = is_local_node; // && is_not_pbc_slave_node; if (count_node) { for (UInt i = 0; i < Model::spatial_dimension; ++i) { if (m(i) > std::numeric_limits::epsilon()) mv2 += v(i) * v(i) * m(i); } } ekin += mv2; } } else if (this->getDOFManager().hasMatrix("M")) { Array Mv(nb_nodes, Model::spatial_dimension); this->getDOFManager().getMatrix("M").matVecMul(*this->velocity, Mv); auto mv_it = Mv.begin(Model::spatial_dimension); auto mv_end = Mv.end(Model::spatial_dimension); auto v_it = this->velocity->begin(Model::spatial_dimension); for (; mv_it != mv_end; ++mv_it, ++v_it) { ekin += v_it->dot(*mv_it); } } else { AKANTU_DEBUG_ERROR("No function called to assemble the mass matrix."); } mesh.getCommunicator().allReduce(ekin, SynchronizerOperation::_sum); AKANTU_DEBUG_OUT(); return ekin * .5; } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getKineticEnergy(const ElementType & type, UInt index) { AKANTU_DEBUG_IN(); UInt nb_quadrature_points = getFEEngine().getNbIntegrationPoints(type); Array vel_on_quad(nb_quadrature_points, Model::spatial_dimension); Array filter_element(1, 1, index); getFEEngine().interpolateOnIntegrationPoints(*velocity, vel_on_quad, Model::spatial_dimension, type, _not_ghost, filter_element); Array::vector_iterator vit = vel_on_quad.begin(Model::spatial_dimension); Array::vector_iterator vend = vel_on_quad.end(Model::spatial_dimension); Vector rho_v2(nb_quadrature_points); Real rho = materials[material_index(type)(index)]->getRho(); for (UInt q = 0; vit != vend; ++vit, ++q) { rho_v2(q) = rho * vit->dot(*vit); } AKANTU_DEBUG_OUT(); return .5 * getFEEngine().integrate(rho_v2, type, index); } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getExternalWork() { AKANTU_DEBUG_IN(); auto ext_force_it = external_force->begin(Model::spatial_dimension); auto int_force_it = internal_force->begin(Model::spatial_dimension); auto boun_it = blocked_dofs->begin(Model::spatial_dimension); decltype(ext_force_it) incr_or_velo_it; if (this->method == _static) { incr_or_velo_it = this->displacement_increment->begin(Model::spatial_dimension); } else { incr_or_velo_it = this->velocity->begin(Model::spatial_dimension); } Real work = 0.; UInt nb_nodes = this->mesh.getNbNodes(); for (UInt n = 0; n < nb_nodes; ++n, ++ext_force_it, ++int_force_it, ++boun_it, ++incr_or_velo_it) { const auto & int_force = *int_force_it; const auto & ext_force = *ext_force_it; const auto & boun = *boun_it; const auto & incr_or_velo = *incr_or_velo_it; bool is_local_node = this->mesh.isLocalOrMasterNode(n); // bool is_not_pbc_slave_node = !this->isPBCSlaveNode(n); bool count_node = is_local_node; // && is_not_pbc_slave_node; if (count_node) { for (UInt i = 0; i < Model::spatial_dimension; ++i) { if (boun(i)) work -= int_force(i) * incr_or_velo(i); else work += ext_force(i) * incr_or_velo(i); } } } mesh.getCommunicator().allReduce(work, SynchronizerOperation::_sum); if (this->method != _static) work *= this->getTimeStep(); AKANTU_DEBUG_OUT(); return work; } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getEnergy(const std::string & energy_id) { AKANTU_DEBUG_IN(); if (energy_id == "kinetic") { return getKineticEnergy(); } else if (energy_id == "external work") { return getExternalWork(); } Real energy = 0.; for (auto & material : materials) energy += material->getEnergy(energy_id); /// reduction sum over all processors mesh.getCommunicator().allReduce(energy, SynchronizerOperation::_sum); AKANTU_DEBUG_OUT(); return energy; } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getEnergy(const std::string & energy_id, const ElementType & type, UInt index) { AKANTU_DEBUG_IN(); if (energy_id == "kinetic") { return getKineticEnergy(type, index); } UInt mat_index = this->material_index(type, _not_ghost)(index); UInt mat_loc_num = this->material_local_numbering(type, _not_ghost)(index); Real energy = this->materials[mat_index]->getEnergy(energy_id, type, mat_loc_num); AKANTU_DEBUG_OUT(); return energy; } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::onElementsAdded(const Array & element_list, const NewElementsEvent & event) { AKANTU_DEBUG_IN(); this->material_index.initialize(mesh, _element_kind = _ek_not_defined, + _with_nb_element = true, _default_value = UInt(-1)); this->material_local_numbering.initialize( - mesh, _element_kind = _ek_not_defined, _default_value = UInt(-1)); + mesh, _element_kind = _ek_not_defined, _with_nb_element = true, + _default_value = UInt(-1)); ElementTypeMapArray filter("new_element_filter", this->getID(), this->getMemoryID()); for (auto & elem : element_list) { if (!filter.exists(elem.type, elem.ghost_type)) filter.alloc(0, 1, elem.type, elem.ghost_type); filter(elem.type, elem.ghost_type).push_back(elem.element); } this->assignMaterialToElements(&filter); for (auto & material : materials) material->onElementsAdded(element_list, event); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::onElementsRemoved( const Array & element_list, const ElementTypeMapArray & new_numbering, const RemovedElementsEvent & event) { for (auto & material : materials) { material->onElementsRemoved(element_list, new_numbering, event); } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::onNodesAdded(const Array & nodes_list, const NewNodesEvent & event) { AKANTU_DEBUG_IN(); UInt nb_nodes = mesh.getNbNodes(); if (displacement) { displacement->resize(nb_nodes, 0.); ++displacement_release; } if (mass) mass->resize(nb_nodes, 0.); if (velocity) velocity->resize(nb_nodes, 0.); if (acceleration) acceleration->resize(nb_nodes, 0.); if (external_force) external_force->resize(nb_nodes, 0.); if (internal_force) internal_force->resize(nb_nodes, 0.); if (blocked_dofs) blocked_dofs->resize(nb_nodes, 0.); if (current_position) current_position->resize(nb_nodes, 0.); if (previous_displacement) previous_displacement->resize(nb_nodes, 0.); if (displacement_increment) displacement_increment->resize(nb_nodes, 0.); for (auto & material : materials) { material->onNodesAdded(nodes_list, event); } need_to_reassemble_lumped_mass = true; need_to_reassemble_mass = true; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::onNodesRemoved(const Array & /*element_list*/, const Array & new_numbering, const RemovedNodesEvent & /*event*/) { if (displacement) { mesh.removeNodesFromArray(*displacement, new_numbering); ++displacement_release; } if (mass) mesh.removeNodesFromArray(*mass, new_numbering); if (velocity) mesh.removeNodesFromArray(*velocity, new_numbering); if (acceleration) mesh.removeNodesFromArray(*acceleration, new_numbering); if (internal_force) mesh.removeNodesFromArray(*internal_force, new_numbering); if (external_force) mesh.removeNodesFromArray(*external_force, new_numbering); if (blocked_dofs) mesh.removeNodesFromArray(*blocked_dofs, new_numbering); // if (increment_acceleration) // mesh.removeNodesFromArray(*increment_acceleration, new_numbering); if (displacement_increment) mesh.removeNodesFromArray(*displacement_increment, new_numbering); if (previous_displacement) mesh.removeNodesFromArray(*previous_displacement, new_numbering); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::printself(std::ostream & stream, int indent) const { std::string space; for (Int i = 0; i < indent; i++, space += AKANTU_INDENT) ; stream << space << "Solid Mechanics Model [" << std::endl; stream << space << " + id : " << id << std::endl; stream << space << " + spatial dimension : " << Model::spatial_dimension << std::endl; stream << space << " + fem [" << std::endl; getFEEngine().printself(stream, indent + 2); stream << space << AKANTU_INDENT << "]" << std::endl; stream << space << " + nodals information [" << std::endl; displacement->printself(stream, indent + 2); if (velocity) velocity->printself(stream, indent + 2); if (acceleration) acceleration->printself(stream, indent + 2); if (mass) mass->printself(stream, indent + 2); external_force->printself(stream, indent + 2); internal_force->printself(stream, indent + 2); blocked_dofs->printself(stream, indent + 2); stream << space << AKANTU_INDENT << "]" << std::endl; stream << space << " + material information [" << std::endl; material_index.printself(stream, indent + 2); stream << space << AKANTU_INDENT << "]" << std::endl; stream << space << " + materials [" << std::endl; for (auto & material : materials) { material->printself(stream, indent + 1); } stream << space << AKANTU_INDENT << "]" << std::endl; stream << space << "]" << std::endl; } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::initializeNonLocal() { this->non_local_manager->synchronize(*this, _gst_material_id); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::insertIntegrationPointsInNeighborhoods( const GhostType & ghost_type) { for (auto & mat : materials) { MaterialNonLocalInterface * mat_non_local; if ((mat_non_local = dynamic_cast(mat.get())) == nullptr) continue; ElementTypeMapArray quadrature_points_coordinates( "quadrature_points_coordinates_tmp_nl", this->id, this->memory_id); quadrature_points_coordinates.initialize(this->getFEEngine(), _nb_component = spatial_dimension, _ghost_type = ghost_type); for (auto & type : quadrature_points_coordinates.elementTypes( Model::spatial_dimension, ghost_type)) { this->getFEEngine().computeIntegrationPointsCoordinates( quadrature_points_coordinates(type, ghost_type), type, ghost_type); } mat_non_local->initMaterialNonLocal(); mat_non_local->insertIntegrationPointsInNeighborhoods( ghost_type, quadrature_points_coordinates); } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::computeNonLocalStresses( const GhostType & ghost_type) { for (auto & mat : materials) { if (dynamic_cast(mat.get()) == nullptr) continue; auto & mat_non_local = dynamic_cast(*mat); mat_non_local.computeNonLocalStresses(ghost_type); } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::updateLocalInternal( ElementTypeMapReal & internal_flat, const GhostType & ghost_type, const ElementKind & kind) { const ID field_name = internal_flat.getName(); for (auto & material : materials) { if (material->isInternal(field_name, kind)) material->flattenInternal(field_name, internal_flat, ghost_type, kind); } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::updateNonLocalInternal( ElementTypeMapReal & internal_flat, const GhostType & ghost_type, const ElementKind & kind) { const ID field_name = internal_flat.getName(); for (auto & mat : materials) { if (dynamic_cast(mat.get()) == nullptr) continue; auto & mat_non_local = dynamic_cast(*mat); mat_non_local.updateNonLocalInternals(internal_flat, field_name, ghost_type, kind); } } /* -------------------------------------------------------------------------- */ FEEngine & SolidMechanicsModel::getFEEngineBoundary(const ID & name) { return dynamic_cast( getFEEngineClassBoundary(name)); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::splitElementByMaterial( const Array & elements, std::vector> & elements_per_mat) const { ElementType current_element_type = _not_defined; GhostType current_ghost_type = _casper; const Array * mat_indexes = nullptr; const Array * mat_loc_num = nullptr; for (auto && el : elements) { if (el.type != current_element_type || el.ghost_type != current_ghost_type) { current_element_type = el.type; current_ghost_type = el.ghost_type; mat_indexes = &(this->material_index(el.type, el.ghost_type)); mat_loc_num = &(this->material_local_numbering(el.type, el.ghost_type)); } UInt old_id = el.element; el.element = (*mat_loc_num)(old_id); elements_per_mat[(*mat_indexes)(old_id)].push_back(el); } } /* -------------------------------------------------------------------------- */ UInt SolidMechanicsModel::getNbData(const Array & elements, const SynchronizationTag & tag) const { AKANTU_DEBUG_IN(); UInt size = 0; UInt nb_nodes_per_element = 0; for (const Element & el : elements) { nb_nodes_per_element += Mesh::getNbNodesPerElement(el.type); } switch (tag) { case _gst_material_id: { size += elements.size() * sizeof(UInt); break; } case _gst_smm_mass: { size += nb_nodes_per_element * sizeof(Real) * Model::spatial_dimension; // mass vector break; } case _gst_smm_for_gradu: { size += nb_nodes_per_element * Model::spatial_dimension * sizeof(Real); // displacement break; } case _gst_smm_boundary: { // force, displacement, boundary size += nb_nodes_per_element * Model::spatial_dimension * (2 * sizeof(Real) + sizeof(bool)); break; } case _gst_for_dump: { // displacement, velocity, acceleration, residual, force size += nb_nodes_per_element * Model::spatial_dimension * sizeof(Real) * 5; break; } default: {} } if (tag != _gst_material_id) { splitByMaterial(elements, [&](auto && mat, auto && elements) { size += mat.getNbData(elements, tag); }); } AKANTU_DEBUG_OUT(); return size; } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::packData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) const { AKANTU_DEBUG_IN(); switch (tag) { case _gst_material_id: { this->packElementalDataHelper(material_index, buffer, elements, false, getFEEngine()); break; } case _gst_smm_mass: { packNodalDataHelper(*mass, buffer, elements, mesh); break; } case _gst_smm_for_gradu: { packNodalDataHelper(*displacement, buffer, elements, mesh); break; } case _gst_for_dump: { packNodalDataHelper(*displacement, buffer, elements, mesh); packNodalDataHelper(*velocity, buffer, elements, mesh); packNodalDataHelper(*acceleration, buffer, elements, mesh); packNodalDataHelper(*internal_force, buffer, elements, mesh); packNodalDataHelper(*external_force, buffer, elements, mesh); break; } case _gst_smm_boundary: { packNodalDataHelper(*external_force, buffer, elements, mesh); packNodalDataHelper(*velocity, buffer, elements, mesh); packNodalDataHelper(*blocked_dofs, buffer, elements, mesh); break; } default: {} } if (tag != _gst_material_id) { splitByMaterial(elements, [&](auto && mat, auto && elements) { mat.packData(buffer, elements, tag); }); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::unpackData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) { AKANTU_DEBUG_IN(); switch (tag) { case _gst_material_id: { for (auto && element : elements) { UInt recv_mat_index; buffer >> recv_mat_index; UInt & mat_index = material_index(element); if (mat_index != UInt(-1)) continue; // add ghosts element to the correct material mat_index = recv_mat_index; UInt index = materials[mat_index]->addElement(element); material_local_numbering(element) = index; } break; } case _gst_smm_mass: { unpackNodalDataHelper(*mass, buffer, elements, mesh); break; } case _gst_smm_for_gradu: { unpackNodalDataHelper(*displacement, buffer, elements, mesh); break; } case _gst_for_dump: { unpackNodalDataHelper(*displacement, buffer, elements, mesh); unpackNodalDataHelper(*velocity, buffer, elements, mesh); unpackNodalDataHelper(*acceleration, buffer, elements, mesh); unpackNodalDataHelper(*internal_force, buffer, elements, mesh); unpackNodalDataHelper(*external_force, buffer, elements, mesh); break; } case _gst_smm_boundary: { unpackNodalDataHelper(*external_force, buffer, elements, mesh); unpackNodalDataHelper(*velocity, buffer, elements, mesh); unpackNodalDataHelper(*blocked_dofs, buffer, elements, mesh); break; } default: {} } if (tag != _gst_material_id) { splitByMaterial(elements, [&](auto && mat, auto && elements) { mat.unpackData(buffer, elements, tag); }); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ UInt SolidMechanicsModel::getNbData(const Array & dofs, const SynchronizationTag & tag) const { AKANTU_DEBUG_IN(); UInt size = 0; // UInt nb_nodes = mesh.getNbNodes(); switch (tag) { case _gst_smm_uv: { size += sizeof(Real) * Model::spatial_dimension * 2; break; } case _gst_smm_res: { size += sizeof(Real) * Model::spatial_dimension; break; } case _gst_smm_mass: { size += sizeof(Real) * Model::spatial_dimension; break; } case _gst_for_dump: { size += sizeof(Real) * Model::spatial_dimension * 5; break; } default: { AKANTU_DEBUG_ERROR("Unknown ghost synchronization tag : " << tag); } } AKANTU_DEBUG_OUT(); return size * dofs.size(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::packData(CommunicationBuffer & buffer, const Array & dofs, const SynchronizationTag & tag) const { AKANTU_DEBUG_IN(); switch (tag) { case _gst_smm_uv: { packDOFDataHelper(*displacement, buffer, dofs); packDOFDataHelper(*velocity, buffer, dofs); break; } case _gst_smm_res: { packDOFDataHelper(*internal_force, buffer, dofs); break; } case _gst_smm_mass: { packDOFDataHelper(*mass, buffer, dofs); break; } case _gst_for_dump: { packDOFDataHelper(*displacement, buffer, dofs); packDOFDataHelper(*velocity, buffer, dofs); packDOFDataHelper(*acceleration, buffer, dofs); packDOFDataHelper(*internal_force, buffer, dofs); packDOFDataHelper(*external_force, buffer, dofs); break; } default: { AKANTU_DEBUG_ERROR("Unknown ghost synchronization tag : " << tag); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::unpackData(CommunicationBuffer & buffer, const Array & dofs, const SynchronizationTag & tag) { AKANTU_DEBUG_IN(); switch (tag) { case _gst_smm_uv: { unpackDOFDataHelper(*displacement, buffer, dofs); unpackDOFDataHelper(*velocity, buffer, dofs); break; } case _gst_smm_res: { unpackDOFDataHelper(*internal_force, buffer, dofs); break; } case _gst_smm_mass: { unpackDOFDataHelper(*mass, buffer, dofs); break; } case _gst_for_dump: { unpackDOFDataHelper(*displacement, buffer, dofs); unpackDOFDataHelper(*velocity, buffer, dofs); unpackDOFDataHelper(*acceleration, buffer, dofs); unpackDOFDataHelper(*internal_force, buffer, dofs); unpackDOFDataHelper(*external_force, buffer, dofs); break; } default: { AKANTU_DEBUG_ERROR("Unknown ghost synchronization tag : " << tag); } } AKANTU_DEBUG_OUT(); } } // namespace akantu 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 d0f6997bf..bcd0a3548 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,752 +1,745 @@ /** * @file solid_mechanics_model_cohesive.cc * * @author Mauro Corrado * @author Nicolas Richart * @author Marco Vocialta * * @date creation: Tue May 08 2012 * @date last modification: Wed Jan 13 2016 * * @brief Solid mechanics model for cohesive elements * * @section LICENSE * * Copyright (©) 2010-2012, 2014, 2015 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 "integrator_gauss.hh" #include "material_cohesive.hh" #include "parser.hh" #include "shape_cohesive.hh" #include "dumpable_inline_impl.hh" #ifdef AKANTU_USE_IOHELPER #include "dumper_paraview.hh" #endif /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ namespace akantu { const SolidMechanicsModelCohesiveOptions default_solid_mechanics_model_cohesive_options(_explicit_lumped_mass, false); /* -------------------------------------------------------------------------- */ SolidMechanicsModelCohesive::SolidMechanicsModelCohesive( Mesh & mesh, UInt dim, const ID & id, const MemoryID & memory_id) : SolidMechanicsModel(mesh, dim, id, memory_id), tangents("tangents", id), facet_stress("facet_stress", id), facet_material("facet_material", id) { AKANTU_DEBUG_IN(); this->model_type = ModelType::_solid_mechanics_model_cohesive; auto && tmp_material_selector = std::make_shared(*this); tmp_material_selector->setFallback(this->material_selector); 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 cohesive_synchronizer = std::make_unique( mesh, "cohesive_distributed_synchronizer"); auto & synchronizer = mesh.getElementSynchronizer(); 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); } 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.extrinsic; inserter = std::make_unique( mesh, is_extrinsic, id + ":cohesive_element_inserter"); if (mesh.isDistributed()) { auto & mesh_facets = inserter->getMeshFacets(); auto & synchronizer = dynamic_cast(mesh_facets.getElementSynchronizer()); this->registerSynchronizer(synchronizer, _gst_smmc_facets); this->registerSynchronizer(synchronizer, _gst_smmc_facets_conn); synchronizeGhostFacetsConnectivity(); /// create the facet synchronizer for extrinsic simulations if (is_extrinsic) { facet_stress_synchronizer = std::make_unique(mesh_facets); facet_stress_synchronizer->updateSchemes( [&](auto & scheme, auto & proc, auto & direction) { scheme.copy(const_cast(synchronizer) .getCommunications() .getScheme(proc, direction)); }); this->registerSynchronizer(*facet_stress_synchronizer, _gst_smmc_facets_stress); } inserter->initParallel(*cohesive_synchronizer); } 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 = 0; while ((dynamic_cast(materials[cohesive_index].get()) == nullptr) && cohesive_index <= materials.size()) ++cohesive_index; AKANTU_DEBUG_ASSERT(cohesive_index != materials.size(), "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 if (is_extrinsic) { this->initExtrinsicMaterials(); } else { this->initIntrinsicMaterials(); } AKANTU_DEBUG_OUT(); } // namespace akantu /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::initExtrinsicMaterials() { const Mesh & mesh_facets = inserter->getMeshFacets(); - facet_material.initialize(mesh_facets, - _spatial_dimension = Model::spatial_dimension - 1); + facet_material.initialize( + mesh_facets, _spatial_dimension = spatial_dimension - 1, + _with_nb_element = true, + _default_value = material_selector->getFallbackValue()); for_each_elements(mesh_facets, [&](auto && element) { auto mat_index = (*material_selector)(element); auto & mat = dynamic_cast( *materials[mat_index]); facet_material(element) = mat_index; mat.addFacet(element); }, _spatial_dimension = spatial_dimension - 1); SolidMechanicsModel::initMaterials(); this->initAutomaticInsertion(); } /* -------------------------------------------------------------------------- */ #if 0 void SolidMechanicsModelCohesive::initIntrinsicCohesiveMaterials( const std::string & cohesive_surfaces) { AKANTU_DEBUG_IN(); #if defined(AKANTU_PARALLEL_COHESIVE_ELEMENT) if (facet_synchronizer != nullptr) inserter->initParallel(facet_synchronizer, cohesive_element_synchronizer); // inserter->initParallel(facet_synchronizer, synch_parallel); #endif std::istringstream split(cohesive_surfaces); std::string physname; while (std::getline(split, physname, ',')) { AKANTU_DEBUG_INFO( "Pre-inserting cohesive elements along facets from physical surface: " << physname); insertElementsFromMeshData(physname); } synchronizeInsertionData(); SolidMechanicsModel::initMaterials(); auto && tmp = std::make_shared(*this); tmp->setFallback(material_selector->getFallbackValue()); tmp->setFallback(material_selector->getFallbackSelector()); material_selector = tmp; // UInt nb_new_elements = inserter->insertElements(); // if (nb_new_elements > 0) { // this->reinitializeSolver(); // } AKANTU_DEBUG_OUT(); } #endif /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::synchronizeInsertionData() { if (inserter->getMeshFacets().isDistributed()) { inserter->getMeshFacets().getElementSynchronizer().synchronizeOnce( *inserter, _gst_ce_groups); } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::initIntrinsicMaterials() { AKANTU_DEBUG_IN(); SolidMechanicsModel::initMaterials(); inserter->insertElements(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /** * Initialize the model,basically it pre-compute the shapes, shapes derivatives * and jacobian */ void SolidMechanicsModelCohesive::initModel() { AKANTU_DEBUG_IN(); SolidMechanicsModel::initModel(); registerFEEngineObject("CohesiveFEEngine", mesh, Model::spatial_dimension); /// 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 Array & connectivity = - // mesh.getConnectivity(tmp_type, type_ghost); - // if (connectivity.size() == 0) - // continue; - - // type = tmp_type; - // ElementType type_facet = Mesh::getFacetType(type); - // ElementType type_cohesive = - // FEEngine::getCohesiveElementType(type_facet); - // mesh.addConnectivityType(type_cohesive, type_ghost); - // } - // } - // AKANTU_DEBUG_ASSERT(type != _not_defined, "No elements in the mesh"); + 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); registerFEEngineObject( "FacetsFEEngine", mesh.getMeshFacets(), Model::spatial_dimension - 1); 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::updateFacetStressSynchronizer() { -if (facet_stress_synchronizer != nullptr) { + if (facet_stress_synchronizer != nullptr) { const auto & rank_to_element = mesh.getElementSynchronizer().getElementToRank(); const auto & facet_checks = inserter->getCheckFacets(); const auto & element_to_facet = mesh.getElementToSubelement(); UInt rank = mesh.getCommunicator().whoAmI(); facet_stress_synchronizer->updateSchemes( - [&](auto & scheme, auto & proc, auto & direction) { + [&](auto & scheme, auto & proc, auto & /*direction*/) { UInt el = 0; for (auto && element : scheme) { if (not facet_checks(element)) continue; const auto & next_el = element_to_facet(element); UInt rank_left = rank_to_element(next_el[0]); UInt rank_right = rank_to_element(next_el[1]); if ((rank_left == rank and rank_right == proc) or (rank_left == proc and rank_right == rank)) { scheme[el] = element; ++el; } } scheme.resize(el); }); } } - /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::initAutomaticInsertion() { AKANTU_DEBUG_IN(); this->updateFacetStressSynchronizer(); - this->facet_stress.initialize(inserter->getMeshFacets(), - _nb_component = - 2 * spatial_dimension * spatial_dimension, - _spatial_dimension = spatial_dimension - 1); - resizeFacetStress(); + this->resizeFacetStress(); /// compute normals on facets - computeNormals(); + this->computeNormals(); - initStressInterpolation(); + this->initStressInterpolation(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::updateAutomaticInsertion() { AKANTU_DEBUG_IN(); - inserter->limitCheckFacets(); + 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 Array & facet_to_element = mesh_facets.getSubelementToElement(type, elem_gt); UInt nb_facet_per_elem = facet_to_element.getNbComponent(); - Array & el_q_facet = elements_quad_facets(type, elem_gt); - ElementType facet_type = Mesh::getFacetType(type); - UInt nb_quad_per_facet = getFEEngine("FacetsFEEngine").getNbIntegrationPoints(facet_type); el_q_facet.resize(nb_element * nb_facet_per_elem * nb_quad_per_facet); for (UInt el = 0; el < nb_element; ++el) { for (UInt f = 0; f < nb_facet_per_elem; ++f) { Element global_facet_elem = facet_to_element(el, f); UInt global_facet = global_facet_elem.element; GhostType facet_gt = global_facet_elem.ghost_type; const Array & quad_f = quad_facets(facet_type, facet_gt); 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 (UInt m = 0; m < materials.size(); ++m) { - try { - auto & mat __attribute__((unused)) = - dynamic_cast(*materials[m]); - } catch (std::bad_cast &) { - /// initialize the interpolation function - materials[m]->initElementalFieldInterpolation(elements_quad_facets); - } + 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(); if (isDefaultSolverExplicit()) { for (auto & material : materials) { try { auto & mat = dynamic_cast(*material); mat.computeEnergies(); } catch (std::bad_cast & bce) { } } } 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) { try { auto & mat __attribute__((unused)) = dynamic_cast(*material); } catch (std::bad_cast &) { /// interpolate stress on facet quadrature points positions material->interpolateStressOnFacets(facet_stress, by_elem_result); } } #if defined(AKANTU_DEBUG_TOOLS) debug::element_manager.printData( debug::_dm_model_cohesive, "Interpolated stresses before", facet_stress); #endif this->synchronize(_gst_smmc_facets_stress); #if defined(AKANTU_DEBUG_TOOLS) debug::element_manager.printData(debug::_dm_model_cohesive, "Interpolated stresses", facet_stress); #endif } /* -------------------------------------------------------------------------- */ UInt SolidMechanicsModelCohesive::checkCohesiveStress() { interpolateStress(); for (auto & mat : materials) { try { auto & mat_cohesive = dynamic_cast(*mat); /// check which not ghost cohesive elements are to be created mat_cohesive.checkInsertion(); } catch (std::bad_cast &) { } } /// 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(); // } return nb_new_elements; } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::onElementsAdded( const Array & element_list, const NewElementsEvent & event) { AKANTU_DEBUG_IN(); updateCohesiveSynchronizers(); SolidMechanicsModel::onElementsAdded(element_list, event); if (cohesive_synchronizer != nullptr) cohesive_synchronizer->computeAllBufferSizes(*this); if (is_extrinsic) resizeFacetStress(); /// if (method != _explicit_lumped_mass) { /// this->initSolver(); /// } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::onNodesAdded(const Array & new_nodes, const NewNodesEvent & event) { AKANTU_DEBUG_IN(); // Array nodes_list(nb_new_nodes); // for (UInt n = 0; n < nb_new_nodes; ++n) // nodes_list(n) = doubled_nodes(n, 1); SolidMechanicsModel::onNodesAdded(new_nodes, event); UInt new_node, old_node; try { const auto & cohesive_event = dynamic_cast(event); const auto & old_nodes = cohesive_event.getOldNodesList(); auto copy = [this, &new_node, &old_node](auto & arr) { for (UInt s = 0; s < spatial_dimension; ++s) { arr(new_node, s) = arr(old_node, s); } }; for (auto && pair : zip(new_nodes, old_nodes)) { std::tie(new_node, old_node) = pair; 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 (this->getDOFManager().hasMatrix("M")) { // this->assembleMass(old_nodes); // } // if (this->getDOFManager().hasLumpedMatrix("M")) { // this->assembleMassLumped(old_nodes); // } } catch (std::bad_cast &) { } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::onEndSolveStep(const AnalysisMethod &) { 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); } 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(); - Mesh & mesh_facets = inserter->getMeshFacets(); + 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); + // 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 nb_quadrature_points = getFEEngine("FacetsFEEngine") + // .getNbIntegrationPoints(type, + // ghost_type); - UInt new_size = nb_facet * nb_quadrature_points; + // UInt new_size = nb_facet * nb_quadrature_points; - facet_stress(type, ghost_type).resize(new_size); - } - } + // 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/src/model/structural_mechanics/structural_mechanics_model_mass.cc b/src/model/structural_mechanics/structural_mechanics_model_mass.cc index 82046c06e..d3908068b 100644 --- a/src/model/structural_mechanics/structural_mechanics_model_mass.cc +++ b/src/model/structural_mechanics/structural_mechanics_model_mass.cc @@ -1,77 +1,77 @@ /** * @file structural_mechanics_model_mass.cc * * @author Sébastien Hartmann * * @date creation: Mon Jul 07 2014 * @date last modification: Thu Oct 15 2015 * * @brief function handling mass computation * * @section LICENSE * * Copyright (©) 2014, 2015 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 "structural_mechanics_model.hh" #include "material.hh" #include "integrator_gauss.hh" #include "shape_structural.hh" /* -------------------------------------------------------------------------- */ namespace akantu { -class ComputeRhoFunctor { +class ComputeRhoFunctorStruct { public: - explicit ComputeRhoFunctor(const StructuralMechanicsModel & model) + explicit ComputeRhoFunctorStruct(const StructuralMechanicsModel & model) : model(model){}; void operator()(Matrix & rho, const Element & element) const { Real mat_rho = model.getMaterial(element).rho; rho.set(mat_rho); } private: const StructuralMechanicsModel & model; }; /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::assembleMass(){ AKANTU_DEBUG_IN(); assembleMass(_not_ghost); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::assembleMass(GhostType ghost_type) { AKANTU_DEBUG_IN(); MyFEEngineType & fem = getFEEngineClass(); - ComputeRhoFunctor compute_rho(*this); + ComputeRhoFunctorStruct compute_rho(*this); for (auto type : mesh.elementTypes(spatial_dimension, ghost_type, _ek_structural)) { fem.assembleFieldMatrix(compute_rho, "M", "displacement", this->getDOFManager(), type, ghost_type); } AKANTU_DEBUG_OUT(); } } // akantu