diff --git a/src/model/common/dof_manager/dof_manager.cc b/src/model/common/dof_manager/dof_manager.cc index 1e8028ac5..fbce811ba 100644 --- a/src/model/common/dof_manager/dof_manager.cc +++ b/src/model/common/dof_manager/dof_manager.cc @@ -1,1062 +1,1058 @@ /** * Copyright (©) 2015-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * This file is part of Akantu * * 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.hh" #include "communicator.hh" #include "mesh.hh" #include "mesh_utils.hh" #include "node_group.hh" #include "node_synchronizer.hh" #include "non_linear_solver.hh" #include "periodic_node_synchronizer.hh" #include "time_step_solver.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ DOFManager::DOFManager(const ID & id) : id(id), dofs_flag(0, 1, std::string(id + ":dofs_type")), global_equation_number(0, 1, "global_equation_number"), communicator(Communicator::getStaticCommunicator()) {} /* -------------------------------------------------------------------------- */ DOFManager::DOFManager(Mesh & mesh, const ID & id) : id(id), mesh(&mesh), dofs_flag(0, 1, std::string(id + ":dofs_type")), global_equation_number(0, 1, "global_equation_number"), communicator(mesh.getCommunicator()) { this->mesh->registerEventHandler(*this, _ehp_dof_manager); } /* -------------------------------------------------------------------------- */ DOFManager::~DOFManager() = default; /* -------------------------------------------------------------------------- */ std::vector DOFManager::getDOFIDs() const { std::vector keys; for (const auto & dof_data : this->dofs) { keys.push_back(dof_data.first); } return keys; } /* -------------------------------------------------------------------------- */ void DOFManager::assembleElementalArrayLocalArray( const Array & elementary_vect, Array & array_assembeled, ElementType type, GhostType ghost_type, Real scale_factor, const Array & filter_elements) { assembleElementalArrayLocalArray( elementary_vect, array_assembeled, this->mesh->getConnectivity(type, ghost_type), type, ghost_type, scale_factor, filter_elements); } /* -------------------------------------------------------------------------- */ void DOFManager::assembleElementalArrayLocalArray( const Array & elementary_vect, Array & array_assembeled, const Array & connectivity, ElementType type, GhostType ghost_type, Real scale_factor, const Array & filter_elements) { AKANTU_DEBUG_IN(); Int nb_element; auto nb_nodes_per_element = Mesh::getNbNodesPerElement(type); auto nb_degree_of_freedom = elementary_vect.getNbComponent() / nb_nodes_per_element; const Idx * filter_it = nullptr; if (filter_elements != empty_filter) { nb_element = filter_elements.size(); filter_it = filter_elements.data(); } else { nb_element = this->mesh->getNbElement(type, ghost_type); } AKANTU_DEBUG_ASSERT(elementary_vect.size() == nb_element, "The vector elementary_vect(" << elementary_vect.getID() << ") has not the good size."); auto elem_it = make_view(elementary_vect, nb_degree_of_freedom, nb_nodes_per_element) .begin(); auto assemble_it = make_view(array_assembeled, nb_degree_of_freedom).begin(); for (Int el = 0; el < nb_element; ++el, ++elem_it) { auto element = el; if (filter_it != nullptr) { element = *filter_it; } // const Vector & conn = *conn_it; const auto & elemental_val = *elem_it; for (Int n = 0; n < nb_nodes_per_element; ++n) { auto node = connectivity(element, n); auto && assemble = assemble_it[node]; assemble += scale_factor * elemental_val(n); } if (filter_it != nullptr) { ++filter_it; } // else // ++conn_it; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void DOFManager::assembleElementalArrayToResidual( const ID & dof_id, const Array & elementary_vect, ElementType type, GhostType ghost_type, Real scale_factor, const Array & filter_elements) { assembleElementalArrayToResidual( dof_id, elementary_vect, this->mesh->getConnectivity(type, ghost_type), type, ghost_type, scale_factor, filter_elements); } /* -------------------------------------------------------------------------- */ void DOFManager::assembleElementalArrayToResidual( const ID & dof_id, const Array & elementary_vect, const Array & connectivity, ElementType type, GhostType ghost_type, Real scale_factor, const Array & filter_elements) { AKANTU_DEBUG_IN(); auto nb_nodes_per_element = Mesh::getNbNodesPerElement(type); auto nb_degree_of_freedom = elementary_vect.getNbComponent() / nb_nodes_per_element; Array array_localy_assembeled(this->mesh->getNbNodes(), nb_degree_of_freedom); array_localy_assembeled.zero(); this->assembleElementalArrayLocalArray( elementary_vect, array_localy_assembeled, connectivity, type, ghost_type, scale_factor, filter_elements); this->assembleToResidual(dof_id, array_localy_assembeled, 1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void DOFManager::assembleElementalArrayToLumpedMatrix( const ID & dof_id, const Array & elementary_vect, const ID & lumped_mtx, ElementType type, GhostType ghost_type, Real scale_factor, const Array & filter_elements) { assembleElementalArrayToLumpedMatrix( dof_id, elementary_vect, lumped_mtx, this->mesh->getConnectivity(type, ghost_type), type, ghost_type, scale_factor, filter_elements); } /* -------------------------------------------------------------------------- */ void DOFManager::assembleElementalArrayToLumpedMatrix( const ID & dof_id, const Array & elementary_vect, const ID & lumped_mtx, const Array & connectivity, ElementType type, GhostType ghost_type, Real scale_factor, const Array & filter_elements) { AKANTU_DEBUG_IN(); auto nb_nodes_per_element = Mesh::getNbNodesPerElement(type); auto nb_degree_of_freedom = elementary_vect.getNbComponent() / nb_nodes_per_element; Array array_localy_assembeled(this->mesh->getNbNodes(), nb_degree_of_freedom); array_localy_assembeled.zero(); this->assembleElementalArrayLocalArray( elementary_vect, array_localy_assembeled, connectivity, type, ghost_type, scale_factor, filter_elements); this->assembleToLumpedMatrix(dof_id, array_localy_assembeled, lumped_mtx, 1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void DOFManager::assembleElementalMatricesToMatrix( const ID & matrix_id, const ID & dof_id, const Array & elementary_mat, ElementType type, GhostType ghost_type, const MatrixType & elemental_matrix_type, const Array & filter_elements) { assembleElementalMatricesToMatrix(matrix_id, dof_id, elementary_mat, mesh->getConnectivity(type, ghost_type), type, ghost_type, elemental_matrix_type, filter_elements); } /* -------------------------------------------------------------------------- */ void DOFManager::assembleMatMulDOFsToResidual(const ID & A_id, Real scale_factor) { for (auto & pair : this->dofs) { const auto & dof_id = pair.first; auto & dof_data = *pair.second; this->assembleMatMulVectToResidual(dof_id, A_id, *dof_data.dof, scale_factor); } } /* -------------------------------------------------------------------------- */ void DOFManager::splitSolutionPerDOFs() { for (auto && data : this->dofs) { auto & dof_data = *data.second; dof_data.solution.resize(dof_data.dof->size() * dof_data.dof->getNbComponent()); this->getSolutionPerDOFs(data.first, dof_data.solution); } } /* -------------------------------------------------------------------------- */ void DOFManager::getSolutionPerDOFs(const ID & dof_id, Array & solution_array) { AKANTU_DEBUG_IN(); this->getArrayPerDOFs(dof_id, this->getSolution(), solution_array); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void DOFManager::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 DOFManager::assembleToResidual(const ID & dof_id, Array & array_to_assemble, Real scale_factor) { AKANTU_DEBUG_IN(); // this->makeConsistentForPeriodicity(dof_id, array_to_assemble); this->assembleToGlobalArray(dof_id, array_to_assemble, this->getResidual(), scale_factor); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void DOFManager::assembleToLumpedMatrix(const ID & dof_id, Array & array_to_assemble, const ID & lumped_mtx, Real scale_factor) { AKANTU_DEBUG_IN(); // this->makeConsistentForPeriodicity(dof_id, array_to_assemble); auto & lumped = this->getLumpedMatrix(lumped_mtx); this->assembleToGlobalArray(dof_id, array_to_assemble, lumped, scale_factor); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ DOFManager::DOFData::DOFData(const ID & dof_id) : support_type(_dst_generic), group_support("__mesh__"), solution(0, 1, dof_id + ":solution"), local_equation_number(0, 1, dof_id + ":local_equation_number"), associated_nodes(0, 1, dof_id + "associated_nodes") {} /* -------------------------------------------------------------------------- */ DOFManager::DOFData::~DOFData() = default; /* -------------------------------------------------------------------------- */ template auto DOFManager::countDOFsForNodes(const DOFData & dof_data, Int nb_nodes, Func && getNode) { auto nb_local_dofs = nb_nodes; decltype(nb_local_dofs) nb_pure_local = 0; for (auto n : arange(nb_nodes)) { UInt node = getNode(n); // http://www.open-std.org/jtc1/sc22/open/n2356/conv.html // bool are by convention casted to 0 and 1 when promoted to int nb_pure_local += this->mesh->isLocalOrMasterNode(node); nb_local_dofs -= this->mesh->isPeriodicSlave(node); } const auto & dofs_array = *dof_data.dof; nb_pure_local *= dofs_array.getNbComponent(); nb_local_dofs *= dofs_array.getNbComponent(); return std::make_pair(nb_local_dofs, nb_pure_local); } /* -------------------------------------------------------------------------- */ auto DOFManager::getNewDOFDataInternal(const ID & dof_id) -> DOFData & { auto it = this->dofs.find(dof_id); if (it != this->dofs.end()) { AKANTU_EXCEPTION("This dof array has already been registered"); } std::unique_ptr dof_data_ptr = this->getNewDOFData(dof_id); DOFData & dof_data = *dof_data_ptr; this->dofs[dof_id] = std::move(dof_data_ptr); return dof_data; } /* -------------------------------------------------------------------------- */ void DOFManager::registerDOFs(const ID & dof_id, Array & dofs_array, DOFSupportType support_type) { auto & dofs_storage = this->getNewDOFDataInternal(dof_id); dofs_storage.support_type = support_type; this->registerDOFsInternal(dof_id, dofs_array); resizeGlobalArrays(); } /* -------------------------------------------------------------------------- */ void DOFManager::registerDOFs(const ID & dof_id, Array & dofs_array, const ID & support_group) { auto & dofs_storage = this->getNewDOFDataInternal(dof_id); dofs_storage.support_type = _dst_nodal; dofs_storage.group_support = support_group; this->registerDOFsInternal(dof_id, dofs_array); resizeGlobalArrays(); } /* -------------------------------------------------------------------------- */ std::tuple DOFManager::registerDOFsInternal(const ID & dof_id, Array & dofs_array) { DOFData & dof_data = this->getDOFData(dof_id); dof_data.dof = &dofs_array; Int nb_local_dofs = 0; Int nb_pure_local = 0; const auto & support_type = dof_data.support_type; switch (support_type) { case _dst_nodal: { const auto & group = dof_data.group_support; std::function getNode; if (group == "__mesh__") { AKANTU_DEBUG_ASSERT( dofs_array.size() == this->mesh->getNbNodes(), "The array of dof is too short to be associated to nodes."); std::tie(nb_local_dofs, nb_pure_local) = countDOFsForNodes( dof_data, this->mesh->getNbNodes(), [](auto && n) { return n; }); } else { const auto & node_group = this->mesh->getElementGroup(group).getNodeGroup().getNodes(); AKANTU_DEBUG_ASSERT( dofs_array.size() == node_group.size(), "The array of dof is too shot to be associated to nodes."); std::tie(nb_local_dofs, nb_pure_local) = countDOFsForNodes(dof_data, node_group.size(), [&node_group](auto && n) { return node_group(n); }); } break; } case _dst_generic: { nb_local_dofs = nb_pure_local = dofs_array.size() * dofs_array.getNbComponent(); break; } default: { AKANTU_EXCEPTION("This type of dofs is not handled yet."); } } dof_data.local_nb_dofs = nb_local_dofs; dof_data.pure_local_nb_dofs = nb_pure_local; dof_data.ghosts_nb_dofs = nb_local_dofs - nb_pure_local; this->pure_local_system_size += nb_pure_local; this->local_system_size += nb_local_dofs; auto nb_total_pure_local = nb_pure_local; communicator.allReduce(nb_total_pure_local, SynchronizerOperation::_sum); this->system_size += nb_total_pure_local; // updating the dofs data after counting is finished switch (support_type) { case _dst_nodal: { const auto & group = dof_data.group_support; if (group != "__mesh__") { auto & support_nodes = this->mesh->getElementGroup(group).getNodeGroup().getNodes(); this->updateDOFsData( dof_data, nb_local_dofs, nb_pure_local, support_nodes.size(), [&support_nodes](Idx node) -> Idx { return support_nodes[node]; }); } else { this->updateDOFsData(dof_data, nb_local_dofs, nb_pure_local, mesh->getNbNodes(), [](Idx node) -> Idx { return node; }); } break; } case _dst_generic: { this->updateDOFsData(dof_data, nb_local_dofs, nb_pure_local); break; } } return std::make_tuple(nb_local_dofs, nb_pure_local, nb_total_pure_local); } /* -------------------------------------------------------------------------- */ void DOFManager::registerDOFsPrevious(const ID & dof_id, Array & array) { DOFData & dof = this->getDOFData(dof_id); if (dof.previous != nullptr) { AKANTU_EXCEPTION("The previous dofs array for " << dof_id << " has already been registered"); } dof.previous = &array; } /* -------------------------------------------------------------------------- */ void DOFManager::registerDOFsIncrement(const ID & dof_id, Array & array) { DOFData & dof = this->getDOFData(dof_id); if (dof.increment != nullptr) { AKANTU_EXCEPTION("The dofs increment array for " << dof_id << " has already been registered"); } dof.increment = &array; } /* -------------------------------------------------------------------------- */ void DOFManager::registerDOFsDerivative(const ID & dof_id, Int order, Array & dofs_derivative) { auto & dof = this->getDOFData(dof_id); auto & derivatives = dof.dof_derivatives; if (Int(derivatives.size()) < order) { derivatives.resize(order, nullptr); } else { if (derivatives[order - 1] != nullptr) { AKANTU_EXCEPTION("The dof derivatives of order " << order << " already been registered for this dof (" << dof_id << ")"); } } derivatives[order - 1] = &dofs_derivative; } /* -------------------------------------------------------------------------- */ void DOFManager::registerBlockedDOFs(const ID & dof_id, Array & blocked_dofs) { auto & dof = this->getDOFData(dof_id); if (dof.blocked_dofs != nullptr) { AKANTU_EXCEPTION("The blocked dofs array for " << dof_id << " has already been registered"); } dof.blocked_dofs = &blocked_dofs; } /* -------------------------------------------------------------------------- */ SparseMatrix & DOFManager::registerSparseMatrix(const ID & matrix_id, std::unique_ptr & matrix) { auto it = this->matrices.find(matrix_id); if (it != this->matrices.end()) { AKANTU_EXCEPTION("The matrix " << matrix_id << " already exists in " << this->id); } auto & ret = *matrix; this->matrices[matrix_id] = std::move(matrix); return ret; } /* -------------------------------------------------------------------------- */ /// Get an instance of a new SparseMatrix SolverVector & DOFManager::registerLumpedMatrix(const ID & matrix_id, std::unique_ptr & matrix) { auto it = this->lumped_matrices.find(matrix_id); if (it != this->lumped_matrices.end()) { AKANTU_EXCEPTION("The lumped matrix " << matrix_id << " already exists in " << this->id); } auto & ret = *matrix; this->lumped_matrices[matrix_id] = std::move(matrix); ret.resize(); return ret; } /* -------------------------------------------------------------------------- */ NonLinearSolver & DOFManager::registerNonLinearSolver( const ID & non_linear_solver_id, std::unique_ptr & non_linear_solver) { NonLinearSolversMap::const_iterator it = this->non_linear_solvers.find(non_linear_solver_id); if (it != this->non_linear_solvers.end()) { AKANTU_EXCEPTION("The non linear solver " << non_linear_solver_id << " already exists in " << this->id); } NonLinearSolver & ret = *non_linear_solver; this->non_linear_solvers[non_linear_solver_id] = std::move(non_linear_solver); return ret; } /* -------------------------------------------------------------------------- */ TimeStepSolver & DOFManager::registerTimeStepSolver( const ID & time_step_solver_id, std::unique_ptr & time_step_solver) { TimeStepSolversMap::const_iterator it = this->time_step_solvers.find(time_step_solver_id); if (it != this->time_step_solvers.end()) { AKANTU_EXCEPTION("The non linear solver " << time_step_solver_id << " already exists in " << this->id); } TimeStepSolver & ret = *time_step_solver; this->time_step_solvers[time_step_solver_id] = std::move(time_step_solver); return ret; } /* -------------------------------------------------------------------------- */ SparseMatrix & DOFManager::getMatrix(const ID & id) { ID matrix_id = this->id + ":mtx:" + id; SparseMatricesMap::const_iterator it = this->matrices.find(matrix_id); if (it == this->matrices.end()) { AKANTU_SILENT_EXCEPTION("The matrix " << matrix_id << " does not exists in " << this->id); } return *(it->second); } /* -------------------------------------------------------------------------- */ bool DOFManager::hasMatrix(const ID & id) const { ID mtx_id = this->id + ":mtx:" + id; auto it = this->matrices.find(mtx_id); return it != this->matrices.end(); } /* -------------------------------------------------------------------------- */ SolverVector & DOFManager::getLumpedMatrix(const ID & id) { ID matrix_id = this->id + ":lumped_mtx:" + id; LumpedMatricesMap::const_iterator it = this->lumped_matrices.find(matrix_id); if (it == this->lumped_matrices.end()) { AKANTU_SILENT_EXCEPTION("The lumped matrix " << matrix_id << " does not exists in " << this->id); } return *(it->second); } /* -------------------------------------------------------------------------- */ const SolverVector & DOFManager::getLumpedMatrix(const ID & id) const { ID matrix_id = this->id + ":lumped_mtx:" + id; auto it = this->lumped_matrices.find(matrix_id); if (it == this->lumped_matrices.end()) { AKANTU_SILENT_EXCEPTION("The lumped matrix " << matrix_id << " does not exists in " << this->id); } return *(it->second); } /* -------------------------------------------------------------------------- */ bool DOFManager::hasLumpedMatrix(const ID & id) const { ID mtx_id = this->id + ":lumped_mtx:" + id; auto it = this->lumped_matrices.find(mtx_id); return it != this->lumped_matrices.end(); } /* -------------------------------------------------------------------------- */ NonLinearSolver & DOFManager::getNonLinearSolver(const ID & id) { ID non_linear_solver_id = this->id + ":nls:" + id; NonLinearSolversMap::const_iterator it = this->non_linear_solvers.find(non_linear_solver_id); if (it == this->non_linear_solvers.end()) { AKANTU_EXCEPTION("The non linear solver " << non_linear_solver_id << " does not exists in " << this->id); } return *(it->second); } /* -------------------------------------------------------------------------- */ bool DOFManager::hasNonLinearSolver(const ID & id) const { ID solver_id = this->id + ":nls:" + id; auto it = this->non_linear_solvers.find(solver_id); return it != this->non_linear_solvers.end(); } /* -------------------------------------------------------------------------- */ TimeStepSolver & DOFManager::getTimeStepSolver(const ID & id) { ID time_step_solver_id = this->id + ":tss:" + id; TimeStepSolversMap::const_iterator it = this->time_step_solvers.find(time_step_solver_id); if (it == this->time_step_solvers.end()) { AKANTU_EXCEPTION("The non linear solver " << time_step_solver_id << " does not exists in " << this->id); } return *(it->second); } /* -------------------------------------------------------------------------- */ bool DOFManager::hasTimeStepSolver(const ID & solver_id) const { ID time_step_solver_id = this->id + ":tss:" + solver_id; auto it = this->time_step_solvers.find(time_step_solver_id); return it != this->time_step_solvers.end(); } /* -------------------------------------------------------------------------- */ void DOFManager::savePreviousDOFs(const ID & dofs_id) { this->getPreviousDOFs(dofs_id).copy(this->getDOFs(dofs_id)); } /* -------------------------------------------------------------------------- */ void DOFManager::zeroResidual() { this->residual->zero(); } /* -------------------------------------------------------------------------- */ void DOFManager::zeroMatrix(const ID & mtx) { this->getMatrix(mtx).zero(); } /* -------------------------------------------------------------------------- */ void DOFManager::zeroLumpedMatrix(const ID & mtx) { this->getLumpedMatrix(mtx).zero(); } /* -------------------------------------------------------------------------- */ /* Mesh Events */ /* -------------------------------------------------------------------------- */ std::pair DOFManager::updateNodalDOFs(const ID & dof_id, const Array & nodes_list) { auto & dof_data = this->getDOFData(dof_id); Int nb_new_local_dofs, nb_new_pure_local; std::tie(nb_new_local_dofs, nb_new_pure_local) = countDOFsForNodes(dof_data, nodes_list.size(), [&nodes_list](auto && n) { return nodes_list(n); }); this->pure_local_system_size += nb_new_pure_local; this->local_system_size += nb_new_local_dofs; auto nb_new_global = nb_new_pure_local; communicator.allReduce(nb_new_global, SynchronizerOperation::_sum); this->system_size += nb_new_global; dof_data.solution.resize(local_system_size); updateDOFsData(dof_data, nb_new_local_dofs, nb_new_pure_local, nodes_list.size(), [&nodes_list](auto pos) -> UInt { return nodes_list[pos]; }); return std::make_pair(nb_new_local_dofs, nb_new_pure_local); } /* -------------------------------------------------------------------------- */ void DOFManager::resizeGlobalArrays() { // resize all relevant arrays this->residual->resize(); this->solution->resize(); this->data_cache->resize(); for (auto & lumped_matrix : lumped_matrices) { lumped_matrix.second->resize(); } for (auto & matrix : matrices) { matrix.second->clearProfile(); } } /* -------------------------------------------------------------------------- */ void DOFManager::onNodesAdded(const Array & nodes_list, const NewNodesEvent &) { for (auto & pair : this->dofs) { const auto & dof_id = pair.first; auto & dof_data = this->getDOFData(dof_id); if (dof_data.support_type != _dst_nodal) { continue; } const auto & group = dof_data.group_support; if (group == "__mesh__") { this->updateNodalDOFs(dof_id, nodes_list); } else { const auto & node_group = this->mesh->getElementGroup(group).getNodeGroup(); Array new_nodes_list; for (const auto & node : nodes_list) { if (node_group.find(node) != Int(-1)) { new_nodes_list.push_back(node); } } this->updateNodalDOFs(dof_id, new_nodes_list); } } this->resizeGlobalArrays(); } /* -------------------------------------------------------------------------- */ void DOFManager::onMeshIsDistributed(const MeshIsDistributedEvent & /*event*/) { AKANTU_DEBUG_ASSERT(this->mesh != nullptr, "The `Mesh` pointer is not set."); // check if the distributed state of the residual and the mesh are the same. if (this->mesh->isDistributed() != this->residual->isDistributed()) { // TODO: Allow to reallocate the internals, in that case one could actually // react on that event. auto is_or_is_not = [](bool q) { return ((q) ? std::string("is") : std::string("is not")); }; AKANTU_EXCEPTION("There is an inconsistency about the distribution state " "of the `DOFManager`." " It seams that the `Mesh` " << is_or_is_not(this->mesh->isDistributed()) << " distributed, but the `DOFManager`'s residual " << is_or_is_not(this->residual->isDistributed()) << ", which is of type " << debug::demangle(typeid(this->residual).name()) << "."); } } /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ class GlobalDOFInfoDataAccessor : public DataAccessor { public: GlobalDOFInfoDataAccessor(DOFManager::DOFData & dof_data, DOFManager & dof_manager) : dof_data(dof_data), dof_manager(dof_manager) { - for (auto && pair : + for (auto && [dof, node] : zip(dof_data.local_equation_number, dof_data.associated_nodes)) { - Idx node; - Idx dof; - std::tie(dof, node) = pair; dofs_per_node[node].push_back(dof); } } Int getNbData(const Array & nodes, const SynchronizationTag & tag) const override { if (tag == SynchronizationTag::_ask_nodes or tag == SynchronizationTag::_giu_global_conn) { return nodes.size() * dof_data.dof->getNbComponent() * sizeof(Idx); } return 0; } void packData(CommunicationBuffer & buffer, const Array & nodes, const SynchronizationTag & tag) const override { if (tag == SynchronizationTag::_ask_nodes or tag == SynchronizationTag::_giu_global_conn) { for (const auto & node : nodes) { const auto & dofs = dofs_per_node.at(node); for (const auto & dof : dofs) { buffer << dof_manager.global_equation_number(dof); } } } } void unpackData(CommunicationBuffer & buffer, const Array & nodes, const SynchronizationTag & tag) override { if (tag == SynchronizationTag::_ask_nodes or tag == SynchronizationTag::_giu_global_conn) { for (const auto & node : nodes) { const auto & dofs = dofs_per_node[node]; for (const auto & dof : dofs) { Idx global_dof; buffer >> global_dof; AKANTU_DEBUG_ASSERT( (dof_manager.global_equation_number(dof) == -1 or dof_manager.global_equation_number(dof) == global_dof), "This dof already had a global_dof_id which is different from " "the received one. " << dof_manager.global_equation_number(dof) << " != " << global_dof); dof_manager.global_equation_number(dof) = global_dof; dof_manager.global_to_local_mapping[global_dof] = dof; } } } } protected: std::unordered_map> dofs_per_node; DOFManager::DOFData & dof_data; DOFManager & dof_manager; }; /* -------------------------------------------------------------------------- */ auto DOFManager::computeFirstDOFIDs(Int nb_new_local_dofs, Int nb_new_pure_local) { // determine the first local/global dof id to use Int offset = 0; this->communicator.exclusiveScan(nb_new_pure_local, offset); auto first_global_dof_id = this->first_global_dof_id + offset; auto first_local_dof_id = this->local_system_size - nb_new_local_dofs; offset = nb_new_pure_local; this->communicator.allReduce(offset); this->first_global_dof_id += offset; return std::make_pair(first_local_dof_id, first_global_dof_id); } /* -------------------------------------------------------------------------- */ void DOFManager::updateDOFsData(DOFData & dof_data, Int nb_new_local_dofs, Int nb_new_pure_local, Int nb_node, const std::function & getNode) { auto nb_local_dofs_added = nb_node * dof_data.dof->getNbComponent(); auto first_dof_pos = dof_data.local_equation_number.size(); dof_data.local_equation_number.reserve(dof_data.local_equation_number.size() + nb_local_dofs_added); dof_data.associated_nodes.reserve(dof_data.associated_nodes.size() + nb_local_dofs_added); this->dofs_flag.resize(this->local_system_size, NodeFlag::_normal); this->global_equation_number.resize(this->local_system_size, -1); std::unordered_map, Idx> masters_dofs; // update per dof info Int local_eq_num, first_global_dof_id; std::tie(local_eq_num, first_global_dof_id) = computeFirstDOFIDs(nb_new_local_dofs, nb_new_pure_local); for (auto d : arange(nb_local_dofs_added)) { auto node = getNode(d / dof_data.dof->getNbComponent()); auto dof_flag = this->mesh->getNodeFlag(node); dof_data.associated_nodes.push_back(node); auto is_local_dof = this->mesh->isLocalOrMasterNode(node); auto is_periodic_slave = this->mesh->isPeriodicSlave(node); auto is_periodic_master = this->mesh->isPeriodicMaster(node); if (is_periodic_slave) { dof_data.local_equation_number.push_back(-1); continue; } // update equation numbers this->dofs_flag(local_eq_num) = dof_flag; dof_data.local_equation_number.push_back(local_eq_num); if (is_local_dof) { this->global_equation_number(local_eq_num) = first_global_dof_id; this->global_to_local_mapping[first_global_dof_id] = local_eq_num; ++first_global_dof_id; } else { this->global_equation_number(local_eq_num) = -1; } if (is_periodic_master) { auto node = getNode(d / dof_data.dof->getNbComponent()); auto dof = d % dof_data.dof->getNbComponent(); masters_dofs.insert( std::make_pair(std::make_pair(node, dof), local_eq_num)); } ++local_eq_num; } // correct periodic slave equation numbers if (this->mesh->isPeriodic()) { auto assoc_begin = dof_data.associated_nodes.begin(); for (auto d : arange(nb_local_dofs_added)) { auto node = dof_data.associated_nodes(first_dof_pos + d); if (not this->mesh->isPeriodicSlave(node)) { continue; } auto master_node = this->mesh->getPeriodicMaster(node); auto dof = d % dof_data.dof->getNbComponent(); dof_data.local_equation_number(first_dof_pos + d) = masters_dofs[std::make_pair(master_node, dof)]; } } // synchronize the global numbering for slaves nodes if (this->mesh->isDistributed()) { GlobalDOFInfoDataAccessor data_accessor(dof_data, *this); if (this->mesh->isPeriodic()) { mesh->getPeriodicNodeSynchronizer().synchronizeOnce( data_accessor, SynchronizationTag::_giu_global_conn); } auto & node_synchronizer = this->mesh->getNodeSynchronizer(); node_synchronizer.synchronizeOnce(data_accessor, SynchronizationTag::_ask_nodes); } } /* -------------------------------------------------------------------------- */ void DOFManager::updateDOFsData(DOFData & dof_data, Int nb_new_local_dofs, Int nb_new_pure_local) { dof_data.local_equation_number.reserve(dof_data.local_equation_number.size() + nb_new_local_dofs); - Int first_local_dof_id, first_global_dof_id; - std::tie(first_local_dof_id, first_global_dof_id) = + auto [first_local_dof_id, first_global_dof_id] = computeFirstDOFIDs(nb_new_local_dofs, nb_new_pure_local); this->dofs_flag.resize(this->local_system_size, NodeFlag::_normal); this->global_equation_number.resize(this->local_system_size, -1); // update per dof info for (auto _ [[gnu::unused]] : arange(nb_new_local_dofs)) { // update equation numbers this->dofs_flag(first_local_dof_id) = NodeFlag::_normal; dof_data.local_equation_number.push_back(first_local_dof_id); this->global_equation_number(first_local_dof_id) = first_global_dof_id; this->global_to_local_mapping[first_global_dof_id] = first_local_dof_id; ++first_global_dof_id; ++first_local_dof_id; } } /* -------------------------------------------------------------------------- */ void DOFManager::onNodesRemoved(const Array &, const Array &, const RemovedNodesEvent &) {} /* -------------------------------------------------------------------------- */ void DOFManager::onElementsAdded(const Array & /*unused*/, const NewElementsEvent & /*unused*/) {} /* -------------------------------------------------------------------------- */ void DOFManager::onElementsRemoved(const Array &, const ElementTypeMapArray &, const RemovedElementsEvent &) {} /* -------------------------------------------------------------------------- */ void DOFManager::onElementsChanged(const Array &, const Array &, const ElementTypeMapArray &, const ChangedElementsEvent &) {} /* -------------------------------------------------------------------------- */ void DOFManager::updateGlobalBlockedDofs() { this->previous_global_blocked_dofs.copy(this->global_blocked_dofs); this->global_blocked_dofs.reserve(this->local_system_size, 0); this->previous_global_blocked_dofs_release = this->global_blocked_dofs_release; for (auto & pair : dofs) { if (not this->hasBlockedDOFs(pair.first)) { continue; } DOFData & dof_data = *pair.second; for (auto && data : zip(dof_data.getLocalEquationsNumbers(), make_view(*dof_data.blocked_dofs))) { const auto & dof = std::get<0>(data); const auto & is_blocked = std::get<1>(data); if (is_blocked) { this->global_blocked_dofs.push_back(dof); } } } std::sort(this->global_blocked_dofs.begin(), this->global_blocked_dofs.end()); auto last = std::unique(this->global_blocked_dofs.begin(), this->global_blocked_dofs.end()); this->global_blocked_dofs.resize(last - this->global_blocked_dofs.begin()); auto are_equal = global_blocked_dofs.size() == previous_global_blocked_dofs.size() and std::equal(global_blocked_dofs.begin(), global_blocked_dofs.end(), previous_global_blocked_dofs.begin()); if (not are_equal) { ++this->global_blocked_dofs_release; } } /* -------------------------------------------------------------------------- */ void DOFManager::applyBoundary(const ID & matrix_id) { auto & J = this->getMatrix(matrix_id); if (this->jacobian_release == J.getRelease()) { if (this->hasBlockedDOFsChanged()) { J.applyBoundary(); } previous_global_blocked_dofs.copy(global_blocked_dofs); } else { J.applyBoundary(); } this->jacobian_release = J.getRelease(); this->previous_global_blocked_dofs_release = this->global_blocked_dofs_release; } /* -------------------------------------------------------------------------- */ void DOFManager::assembleMatMulVectToGlobalArray(const ID & dof_id, const ID & A_id, const Array & x, SolverVector & array, Real scale_factor) { auto & A = this->getMatrix(A_id); data_cache->resize(); data_cache->zero(); this->assembleToGlobalArray(dof_id, x, *data_cache, 1.); A.matVecMul(*data_cache, array, scale_factor, 1.); } /* -------------------------------------------------------------------------- */ void DOFManager::assembleMatMulVectToResidual(const ID & dof_id, const ID & A_id, const Array & x, Real scale_factor) { assembleMatMulVectToGlobalArray(dof_id, A_id, x, *residual, scale_factor); } } // 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 68b5a47ff..a3068710e 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,741 +1,741 @@ /** * Copyright (©) 2012-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * This file is part of Akantu * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . */ /* -------------------------------------------------------------------------- */ #include "solid_mechanics_model_cohesive.hh" #include "aka_iterators.hh" #include "cohesive_element_inserter.hh" #include "element_synchronizer.hh" #include "facet_synchronizer.hh" #include "fe_engine_template.hh" #include "global_ids_updater.hh" #include "integrator_gauss.hh" #include "material_cohesive.hh" #include "mesh_accessor.hh" #include "mesh_global_data_updater.hh" #include "parser.hh" #include "shape_cohesive.hh" /* -------------------------------------------------------------------------- */ #include "dumper_iohelper_paraview.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ namespace akantu { class CohesiveMeshGlobalDataUpdater : public MeshGlobalDataUpdater { public: CohesiveMeshGlobalDataUpdater(SolidMechanicsModelCohesive & model) : model(model), mesh(model.getMesh()), global_ids_updater(model.getMesh(), model.cohesive_synchronizer.get()) { } /* ------------------------------------------------------------------------ */ std::tuple updateData(NewNodesEvent & nodes_event, NewElementsEvent & elements_event) override { auto * cohesive_nodes_event = dynamic_cast(&nodes_event); if (cohesive_nodes_event == nullptr) { return std::make_tuple(nodes_event.getList().size(), elements_event.getList().size()); } /// update nodes type auto & new_nodes = cohesive_nodes_event->getList(); auto & old_nodes = cohesive_nodes_event->getOldNodesList(); auto local_nb_new_nodes = new_nodes.size(); auto nb_new_nodes = local_nb_new_nodes; if (mesh.isDistributed()) { MeshAccessor mesh_accessor(mesh); auto & nodes_flags = mesh_accessor.getNodesFlags(); auto nb_old_nodes = nodes_flags.size(); nodes_flags.resize(nb_old_nodes + local_nb_new_nodes); for (auto && [old_node, new_node] : zip(old_nodes, new_nodes)) { nodes_flags(new_node) = nodes_flags(old_node); } model.updateCohesiveSynchronizers(elements_event); nb_new_nodes = global_ids_updater.updateGlobalIDs(new_nodes.size()); } auto nb_new_elements = elements_event.getList().size(); const auto & comm = mesh.getCommunicator(); comm.allReduce(nb_new_elements, SynchronizerOperation::_sum); if (nb_new_elements > 0) { mesh.sendEvent(elements_event); } if (nb_new_nodes > 0) { mesh.sendEvent(nodes_event); } return std::make_tuple(nb_new_nodes, nb_new_elements); } private: SolidMechanicsModelCohesive & model; Mesh & mesh; GlobalIdsUpdater global_ids_updater; }; /* -------------------------------------------------------------------------- */ SolidMechanicsModelCohesive::SolidMechanicsModelCohesive( Mesh & mesh, Int dim, const ID & id, const std::shared_ptr & dof_manager) : SolidMechanicsModel(mesh, dim, id, dof_manager, ModelType::_solid_mechanics_model_cohesive), tangents("tangents", id), facet_stress("facet_stress", id), facet_material("facet_material", id) { AKANTU_DEBUG_IN(); registerFEEngineObject("CohesiveFEEngine", mesh, Model::spatial_dimension); auto && tmp_material_selector = std::make_shared(*this); tmp_material_selector->setFallback(this->getConstitutiveLawSelector()); this->setConstitutiveLawSelector(tmp_material_selector); this->mesh.registerDumper("cohesive elements", id); this->mesh.addDumpMeshToDumper("cohesive elements", mesh, Model::spatial_dimension, _not_ghost, _ek_cohesive); if (this->mesh.isDistributed()) { /// create the distributed synchronizer for cohesive elements this->cohesive_synchronizer = std::make_unique( mesh, "cohesive_distributed_synchronizer"); auto & synchronizer = mesh.getElementSynchronizer(); this->cohesive_synchronizer->split(synchronizer, [](auto && el) { return Mesh::getKind(el.type) == _ek_cohesive; }); this->registerSynchronizer(*cohesive_synchronizer, SynchronizationTag::_constitutive_law_id); this->registerSynchronizer(*cohesive_synchronizer, SynchronizationTag::_smm_stress); this->registerSynchronizer(*cohesive_synchronizer, SynchronizationTag::_smm_boundary); } this->inserter = std::make_unique( this->mesh, id + ":cohesive_element_inserter"); registerFEEngineObject( "FacetsFEEngine", mesh.getMeshFacets(), Model::spatial_dimension - 1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::setTimeStep(Real time_step, const ID & solver_id) { SolidMechanicsModel::setTimeStep(time_step, solver_id); this->mesh.getDumper("cohesive elements").setTimeStep(time_step); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::initFullImpl(const ModelOptions & options) { AKANTU_DEBUG_IN(); const auto & smmc_options = aka::as_type(options); this->is_extrinsic = smmc_options.is_extrinsic; inserter->setIsExtrinsic(is_extrinsic); if (mesh.isDistributed()) { auto & mesh_facets = inserter->getMeshFacets(); auto & synchronizer = aka::as_type(mesh_facets.getElementSynchronizer()); // synchronizeGhostFacetsConnectivity(); /// create the facet synchronizer for extrinsic simulations if (is_extrinsic) { facet_stress_synchronizer = std::make_unique( synchronizer, id + ":facet_stress_synchronizer"); facet_stress_synchronizer->swapSendRecv(); this->registerSynchronizer(*facet_stress_synchronizer, SynchronizationTag::_smmc_facets_stress); } } MeshAccessor mesh_accessor(mesh); mesh_accessor.registerGlobalDataUpdater( std::make_unique(*this)); auto && [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::initConstitutiveLaws() { AKANTU_DEBUG_IN(); // make sure the material are instantiated instantiateMaterials(); auto & material_selector = this->getConstitutiveLawSelector(); // set the facet information in the material in case of dynamic insertion // to know what material to call for stress checks const Mesh & mesh_facets = inserter->getMeshFacets(); facet_material.initialize( mesh_facets, _spatial_dimension = spatial_dimension - 1, _with_nb_element = true, _default_value = DefaultMaterialCohesiveSelector::getDefaultCohesiveMaterial(*this)); for_each_element( mesh_facets, [&](auto && element) { auto mat_index = material_selector(element); if (not mat_index) { return; } auto & mat = aka::as_type(this->getConstitutiveLaw(mat_index)); facet_material(element) = mat_index; if (is_extrinsic) { mat.addFacet(element); } }, _spatial_dimension = spatial_dimension - 1, _ghost_type = _not_ghost); SolidMechanicsModel::initConstitutiveLaws(); auto & initial_nodes = mesh.getNodalData("initial_nodes_match"); initial_nodes.resize(mesh.getNbNodes()); for (auto && [node, init_node] : enumerate(initial_nodes)) { init_node = node; } lambda = std::make_unique>(0, 1, "cohesive lambda"); if (lambda) { mesh.getElementalData("initial_nodes_connectivities"); mesh.getElementalData("lambda_connectivities"); - auto nodes_to_lambda = mesh.getNodalData("nodes_to_lambda"); + auto & nodes_to_lambda = mesh.getNodalData("nodes_to_lambda"); nodes_to_lambda.resize(mesh.getNbNodes(), -1); } if (is_extrinsic) { this->initAutomaticInsertion(); } else { this->insertIntrinsicElements(); } /// TODO : INITIALIZE LAMBDA HERE AKANTU_DEBUG_OUT(); } // namespace akantu /* -------------------------------------------------------------------------- */ /** * Initialize the model,basically it pre-compute the shapes, shapes derivatives * and jacobian */ void SolidMechanicsModelCohesive::initModel() { AKANTU_DEBUG_IN(); SolidMechanicsModel::initModel(); /// add cohesive type connectivity ElementType type = _not_defined; for (auto && type_ghost : ghost_types) { for (const auto & tmp_type : mesh.elementTypes(spatial_dimension, type_ghost)) { const auto & connectivity = mesh.getConnectivity(tmp_type, type_ghost); if (connectivity.empty()) { continue; } type = tmp_type; auto type_facet = Mesh::getFacetType(type); auto type_cohesive = FEEngine::getCohesiveElementType(type_facet); AKANTU_DEBUG_ASSERT(Mesh::getKind(type_cohesive) == _ek_cohesive, "The element type " << type_cohesive << " is not a cohesive type"); mesh.addConnectivityType(type_cohesive, type_ghost); } } AKANTU_DEBUG_ASSERT(type != _not_defined, "No elements in the mesh"); this->allocNodalField(this->displacement, spatial_dimension, "displacement"); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::insertIntrinsicElements() { AKANTU_DEBUG_IN(); inserter->insertIntrinsicElements(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::initAutomaticInsertion() { AKANTU_DEBUG_IN(); this->inserter->limitCheckFacets(); this->updateFacetStressSynchronizer(); this->resizeFacetStress(); /// compute normals on facets this->computeNormals(); this->initStressInterpolation(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::updateAutomaticInsertion() { AKANTU_DEBUG_IN(); this->inserter->limitCheckFacets(); this->updateFacetStressSynchronizer(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::initStressInterpolation() { Mesh & mesh_facets = inserter->getMeshFacets(); /// compute quadrature points coordinates on facets Array & position = mesh.getNodes(); ElementTypeMapArray quad_facets("quad_facets", id); quad_facets.initialize(mesh_facets, _nb_component = Model::spatial_dimension, _spatial_dimension = Model::spatial_dimension - 1); 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); for (auto elem_gt : ghost_types) { for (const auto & type : mesh.elementTypes(Model::spatial_dimension, elem_gt)) { auto nb_element = mesh.getNbElement(type, elem_gt); if (nb_element == 0) { continue; } /// compute elements' quadrature points and list of facet /// quadrature points positions by element const auto & facet_to_element = mesh_facets.getSubelementToElement(type, elem_gt); auto & el_q_facet = elements_quad_facets(type, elem_gt); auto facet_type = Mesh::getFacetType(type); auto nb_quad_per_facet = getFEEngine("FacetsFEEngine").getNbIntegrationPoints(facet_type); auto nb_facet_per_elem = facet_to_element.getNbComponent(); // small hack in the loop to skip boundary elements, they are silently // initialized to NaN to see if this causes problems el_q_facet.resize(nb_element * nb_facet_per_elem * nb_quad_per_facet, std::numeric_limits::quiet_NaN()); for (auto && data : zip(make_view(facet_to_element), make_view(el_q_facet, spatial_dimension, nb_quad_per_facet))) { const auto & global_facet = std::get<0>(data); auto & el_q = std::get<1>(data); if (global_facet == ElementNull) { continue; } auto && quad_f = make_view(quad_facets(global_facet.type, global_facet.ghost_type), spatial_dimension, nb_quad_per_facet) .begin()[global_facet.element]; el_q = quad_f; } } } /// loop over non cohesive materials this->for_each_constitutive_law([&](auto && material) { if (aka::is_of_type(material)) { return; } /// initialize the interpolation function material.initElementalFieldInterpolation(elements_quad_facets); }); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::assembleInternalForces() { AKANTU_DEBUG_IN(); // f_int += f_int_cohe this->for_each_constitutive_law([&](auto && material) { try { auto & mat = aka::as_type(material); mat.computeTraction(_not_ghost); } catch (std::bad_cast & bce) { } }); SolidMechanicsModel::assembleInternalForces(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::computeNormals() { AKANTU_DEBUG_IN(); Mesh & mesh_facets = this->inserter->getMeshFacets(); this->getFEEngine("FacetsFEEngine") .computeNormalsOnIntegrationPoints(_not_ghost); /** * @todo store tangents while computing normals instead of * recomputing them as follows: */ /* ------------------------------------------------------------------------ */ UInt tangent_components = Model::spatial_dimension * (Model::spatial_dimension - 1); tangents.initialize(mesh_facets, _nb_component = tangent_components, _spatial_dimension = Model::spatial_dimension - 1); 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); this->for_each_constitutive_law([&](auto && material) { if (not aka::is_of_type(material)) { /// interpolate stress on facet quadrature points positions material.interpolateStressOnFacets(facet_stress, by_elem_result); } }); this->synchronize(SynchronizationTag::_smmc_facets_stress); } /* -------------------------------------------------------------------------- */ UInt SolidMechanicsModelCohesive::checkCohesiveStress() { AKANTU_DEBUG_IN(); if (not is_extrinsic) { AKANTU_EXCEPTION( "This function can only be used for extrinsic cohesive elements"); } interpolateStress(); this->for_each_constitutive_law([&](auto && material) { if (aka::is_of_type(material)) { /// check which not ghost cohesive elements are to be created auto & mat_cohesive = aka::as_type(material); mat_cohesive.checkInsertion(); } }); /// insert cohesive elements auto nb_new_elements = inserter->insertElements(); AKANTU_DEBUG_OUT(); return nb_new_elements; } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::onElementsAdded( const Array & element_list, const NewElementsEvent & event) { AKANTU_DEBUG_IN(); SolidMechanicsModel::onElementsAdded(element_list, event); if (lambda) { auto & initial_nodes_connectivities = mesh.getElementalData("initial_nodes_connectivities"); auto & lambda_connectivities = mesh.getElementalData("lambda_connectivities"); for (auto ghost_type : ghost_types) { - for (auto type : mesh.elementTypes(_element_kind = _ek_cohesive)) { + for (auto type : mesh.elementTypes(_element_kind = _ek_cohesive, + _ghost_type = ghost_type)) { auto size = mesh.getConnectivity(type, ghost_type).size(); if (not initial_nodes_connectivities.exists(type, ghost_type)) { auto underlying_type = Mesh::getFacetType(type); initial_nodes_connectivities.alloc( size, Mesh::getNbNodesPerElement(underlying_type), type, ghost_type, -1); lambda_connectivities.alloc( size, Mesh::getNbNodesPerElement(underlying_type), type, ghost_type, -1); } else { initial_nodes_connectivities(type, ghost_type).resize(size, -1); lambda_connectivities(type, ghost_type).resize(size, -1); } } } // Set some connectivities, it will be corrected at on nodes added for (const auto & el : element_list) { if (Mesh::getKind(el.type) != _ek_cohesive) { continue; } auto && conn = mesh.getConnectivity(el); auto && iconn = initial_nodes_connectivities.get(el); iconn = conn.block(0, 0, iconn.size(), 1); } } if (is_extrinsic) { resizeFacetStress(); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::onNodesAdded(const Array & new_nodes, const NewNodesEvent & event) { AKANTU_DEBUG_IN(); SolidMechanicsModel::onNodesAdded(new_nodes, event); auto & initial_nodes = mesh.getNodalData("initial_nodes_match"); auto old_max_nodes = initial_nodes.size(); initial_nodes.resize(mesh.getNbNodes()); const auto * cohesive_event = dynamic_cast(&event); if (cohesive_event == nullptr) { for (auto && node : new_nodes) { initial_nodes(node) = node; } return; } auto old_nodes = cohesive_event->getOldNodesList(); // getting a corrected old_nodes for multiple doubling for (auto & onode : old_nodes) { while (onode >= old_max_nodes) { auto nnode = new_nodes.find(onode); AKANTU_DEBUG_ASSERT(nnode != -1, "If the old node is bigger than old_max_nodes it " "should also be a new node"); onode = nnode; } } auto copy = [&new_nodes, &old_nodes](auto & arr) { auto it = make_view(arr, arr.getNbComponent()).begin(); for (auto [new_node, old_node] : zip(new_nodes, old_nodes)) { it[new_node] = it[old_node]; } }; 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 (displacement_increment) { copy(*displacement_increment); } copy(getDOFManager().getSolution("displacement")); copy(initial_nodes); // correct connectivities if (lambda) { auto & initial_nodes_connectivities = mesh.getElementalData("initial_nodes_connectivities"); auto & lambda_connectivities = mesh.getElementalData("lambda_connectivities"); auto & nodes_to_lambda = mesh.getNodalData("nodes_to_lambda"); auto lambda_id = lambda->size(); for (auto ghost_type : ghost_types) { for (auto type : initial_nodes_connectivities.elementTypes( - _element_kind = _ek_cohesive)) { + _element_kind = _ek_cohesive, _ghost_type = ghost_type)) { auto & initial_nodes_connectivity = initial_nodes_connectivities(type, ghost_type); auto & lambda_connectivity = lambda_connectivities(type, ghost_type); for (auto && [conn, lambda_conn] : zip(make_view(initial_nodes_connectivity), make_view(lambda_connectivity))) { conn = initial_nodes(conn); auto & ntl = nodes_to_lambda(conn); if (ntl == -1) { ntl = lambda_id; ++lambda_id; } lambda_conn = ntl; } } } lambda->resize(lambda_id); - copy(*lambda); // Not sure if necessary } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::afterSolveStep(bool converged) { AKANTU_DEBUG_IN(); /* * This is required because the Cauchy stress is the stress measure that * is used to check the insertion of cohesive elements */ if (converged) { this->for_each_constitutive_law([](auto && mat) { if (mat.isFiniteDeformation()) { mat.computeAllCauchyStresses(_not_ghost); } }); } SolidMechanicsModel::afterSolveStep(converged); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::printself(std::ostream & stream, int indent) const { std::string space(indent, AKANTU_INDENT); stream << space << "SolidMechanicsModelCohesive [" << "\n"; SolidMechanicsModel::printself(stream, indent + 2); stream << space << "]\n"; } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::resizeFacetStress() { AKANTU_DEBUG_IN(); this->facet_stress.initialize(getFEEngine("FacetsFEEngine"), _nb_component = 2 * spatial_dimension * spatial_dimension, _spatial_dimension = spatial_dimension - 1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::addDumpGroupFieldToDumper( const std::string & dumper_name, const std::string & field_id, const std::string & group_name, ElementKind element_kind, bool padding_flag) { AKANTU_DEBUG_IN(); Int 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