diff --git a/src/model/dof_manager.cc b/src/model/dof_manager.cc index 2af5c8269..795e3182b 100644 --- a/src/model/dof_manager.cc +++ b/src/model/dof_manager.cc @@ -1,607 +1,604 @@ /** * @file dof_manager.cc * * @author Nicolas Richart * * @date creation: Tue Aug 18 2015 * @date last modification: Wed Feb 21 2018 * * @brief Implementation of the common parts of the DOFManagers * * @section LICENSE * * Copyright (©) 2015-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "dof_manager.hh" #include "communicator.hh" #include "element_group.hh" #include "mesh.hh" #include "mesh_utils.hh" #include "node_group.hh" #include "non_linear_solver.hh" #include "sparse_matrix.hh" #include "time_step_solver.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ DOFManager::DOFManager(const ID & id, const MemoryID & memory_id) : Memory(id, memory_id), communicator(Communicator::getStaticCommunicator()) {} /* -------------------------------------------------------------------------- */ DOFManager::DOFManager(Mesh & mesh, const ID & id, const MemoryID & memory_id) : Memory(id, memory_id), mesh(&mesh), local_system_size(0), pure_local_system_size(0), system_size(0), communicator(mesh.getCommunicator()) { this->mesh->registerEventHandler(*this, _ehp_dof_manager); } /* -------------------------------------------------------------------------- */ DOFManager::~DOFManager() = default; /* -------------------------------------------------------------------------- */ // void DOFManager::getEquationsNumbers(const ID &, Array &) { // AKANTU_TO_IMPLEMENT(); // } /* -------------------------------------------------------------------------- */ 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, const ElementType & type, const GhostType & ghost_type, Real scale_factor, const Array & filter_elements) { AKANTU_DEBUG_IN(); UInt nb_element; UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt nb_degree_of_freedom = elementary_vect.getNbComponent() / nb_nodes_per_element; UInt * filter_it = nullptr; if (filter_elements != empty_filter) { nb_element = filter_elements.size(); filter_it = filter_elements.storage(); } 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."); const Array & connectivity = this->mesh->getConnectivity(type, ghost_type); - // Array::const_vector_iterator conn_begin = - // connectivity.begin(nb_nodes_per_element); - // Array::const_vector_iterator conn_it = conn_begin; Array::const_matrix_iterator elem_it = elementary_vect.begin(nb_degree_of_freedom, nb_nodes_per_element); for (UInt el = 0; el < nb_element; ++el, ++elem_it) { UInt element = el; if (filter_it != nullptr) { // conn_it = conn_begin + *filter_it; element = *filter_it; } // const Vector & conn = *conn_it; const Matrix & elemental_val = *elem_it; for (UInt n = 0; n < nb_nodes_per_element; ++n) { UInt offset_node = connectivity(element, n) * nb_degree_of_freedom; Vector assemble(array_assembeled.storage() + offset_node, nb_degree_of_freedom); Vector elem_val = elemental_val(n); assemble.aXplusY(elem_val, scale_factor); } if (filter_it != nullptr) ++filter_it; // else // ++conn_it; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void DOFManager::assembleElementalArrayToResidual( const ID & dof_id, const Array & elementary_vect, const ElementType & type, const GhostType & ghost_type, Real scale_factor, const Array & filter_elements) { AKANTU_DEBUG_IN(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt 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.clear(); this->assembleElementalArrayLocalArray( elementary_vect, array_localy_assembeled, 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, const ElementType & type, const GhostType & ghost_type, Real scale_factor, const Array & filter_elements) { AKANTU_DEBUG_IN(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt 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.clear(); this->assembleElementalArrayLocalArray( elementary_vect, array_localy_assembeled, type, ghost_type, scale_factor, filter_elements); this->assembleToLumpedMatrix(dof_id, array_localy_assembeled, lumped_mtx, 1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ 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); } } /* -------------------------------------------------------------------------- */ DOFManager::DOFData::DOFData(const ID & dof_id) : support_type(_dst_generic), group_support("__mesh__"), dof(nullptr), blocked_dofs(nullptr), increment(nullptr), previous(nullptr), solution(0, 1, dof_id + ":solution"), local_equation_number(0, 1, dof_id + ":local_equation_number") {} /* -------------------------------------------------------------------------- */ DOFManager::DOFData::~DOFData() = default; /* -------------------------------------------------------------------------- */ DOFManager::DOFData & DOFManager::getNewDOFData(const ID & dof_id) { auto it = this->dofs.find(dof_id); if (it != this->dofs.end()) { AKANTU_EXCEPTION("This dof array has already been registered"); } std::unique_ptr dofs_storage = std::make_unique(dof_id); this->dofs[dof_id] = std::move(dofs_storage); return *dofs_storage; } /* -------------------------------------------------------------------------- */ void DOFManager::registerDOFsInternal(const ID & dof_id, Array & dofs_array) { DOFData & dofs_storage = this->getDOFData(dof_id); dofs_storage.dof = &dofs_array; UInt nb_local_dofs = 0; UInt nb_pure_local = 0; const DOFSupportType & support_type = dofs_storage.support_type; switch (support_type) { case _dst_nodal: { UInt nb_nodes = 0; const ID & group = dofs_storage.group_support; NodeGroup * node_group = nullptr; if (group == "__mesh__") { nb_nodes = this->mesh->getNbNodes(); } else { node_group = &this->mesh->getElementGroup(group).getNodeGroup(); nb_nodes = node_group->size(); } nb_local_dofs = nb_nodes; AKANTU_DEBUG_ASSERT( dofs_array.size() == nb_local_dofs, "The array of dof is too shot to be associated to nodes."); for (UInt n = 0; n < nb_nodes; ++n) { UInt node = n; if (node_group) node = node_group->getNodes()(n); nb_pure_local += this->mesh->isLocalOrMasterNode(node) ? 1 : 0; } nb_pure_local *= dofs_array.getNbComponent(); nb_local_dofs *= dofs_array.getNbComponent(); 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."); } } this->pure_local_system_size += nb_pure_local; this->local_system_size += nb_local_dofs; communicator.allReduce(nb_pure_local, SynchronizerOperation::_sum); this->system_size += nb_pure_local; } /* -------------------------------------------------------------------------- */ void DOFManager::registerDOFs(const ID & dof_id, Array & dofs_array, const DOFSupportType & support_type) { DOFData & dofs_storage = this->getNewDOFData(dof_id); dofs_storage.support_type = support_type; this->registerDOFsInternal(dof_id, dofs_array); } /* -------------------------------------------------------------------------- */ void DOFManager::registerDOFs(const ID & dof_id, Array & dofs_array, const ID & support_group) { DOFData & dofs_storage = this->getNewDOFData(dof_id); dofs_storage.support_type = _dst_nodal; dofs_storage.group_support = support_group; this->registerDOFsInternal(dof_id, dofs_array); } /* -------------------------------------------------------------------------- */ 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, UInt order, Array & dofs_derivative) { DOFData & dof = this->getDOFData(dof_id); std::vector *> & derivatives = dof.dof_derivatives; if (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) { DOFData & 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; } /* -------------------------------------------------------------------------- */ void DOFManager::splitSolutionPerDOFs() { auto it = this->dofs.begin(); auto end = this->dofs.end(); for (; it != end; ++it) { DOFData & dof_data = *it->second; dof_data.solution.resize(dof_data.dof->size() * dof_data.dof->getNbComponent()); this->getSolutionPerDOFs(it->first, dof_data.solution); } } /* -------------------------------------------------------------------------- */ SparseMatrix & DOFManager::registerSparseMatrix(const ID & matrix_id, std::unique_ptr & matrix) { SparseMatricesMap::const_iterator it = this->matrices.find(matrix_id); if (it != this->matrices.end()) { AKANTU_EXCEPTION("The matrix " << matrix_id << " already exists in " << this->id); } SparseMatrix & ret = *matrix; this->matrices[matrix_id] = std::move(matrix); return ret; } /* -------------------------------------------------------------------------- */ /// Get an instance of a new SparseMatrix Array & DOFManager::getNewLumpedMatrix(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_EXCEPTION("The lumped matrix " << matrix_id << " already exists in " << this->id); } auto mtx = std::make_unique>(this->local_system_size, 1, matrix_id); this->lumped_matrices[matrix_id] = std::move(mtx); return *this->lumped_matrices[matrix_id]; } /* -------------------------------------------------------------------------- */ 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(); } /* -------------------------------------------------------------------------- */ Array & 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 Array & 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)); } /* -------------------------------------------------------------------------- */ /* Mesh Events */ /* -------------------------------------------------------------------------- */ std::pair DOFManager::updateNodalDOFs(const ID & dof_id, const Array & nodes_list) { auto & dof_data = this->getDOFData(dof_id); UInt nb_new_local_dofs = 0; UInt nb_new_pure_local = 0; nb_new_local_dofs = nodes_list.size(); for (const auto & node : nodes_list) { nb_new_pure_local += this->mesh->isLocalOrMasterNode(node) ? 1 : 0; } const auto & dof_array = *dof_data.dof; nb_new_pure_local *= dof_array.getNbComponent(); nb_new_local_dofs *= dof_array.getNbComponent(); this->pure_local_system_size += nb_new_pure_local; this->local_system_size += nb_new_local_dofs; UInt nb_new_global = nb_new_pure_local; communicator.allReduce(nb_new_global, SynchronizerOperation::_sum); this->system_size += nb_new_global; dof_data.solution.resize(dof_data.solution.size() + nb_new_local_dofs); return std::make_pair(nb_new_local_dofs, nb_new_pure_local); } /* -------------------------------------------------------------------------- */ 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) != UInt(-1)) new_nodes_list.push_back(node); } this->updateNodalDOFs(dof_id, new_nodes_list); } } } /* -------------------------------------------------------------------------- */ void DOFManager::onNodesRemoved(const Array &, const Array &, const RemovedNodesEvent &) {} /* -------------------------------------------------------------------------- */ void DOFManager::onElementsAdded(const Array &, const NewElementsEvent &) {} /* -------------------------------------------------------------------------- */ void DOFManager::onElementsRemoved(const Array &, const ElementTypeMapArray &, const RemovedElementsEvent &) {} /* -------------------------------------------------------------------------- */ void DOFManager::onElementsChanged(const Array &, const Array &, const ElementTypeMapArray &, const ChangedElementsEvent &) {} /* -------------------------------------------------------------------------- */ } // namespace akantu diff --git a/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_exponential.cc b/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_exponential.cc index a2030ce7f..eaa8b5163 100644 --- a/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_exponential.cc +++ b/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_exponential.cc @@ -1,338 +1,338 @@ /** * @file material_cohesive_exponential.cc * * @author Mauro Corrado * @author Seyedeh Mohadeseh Taheri Mousavi * @author Marco Vocialta * * @date creation: Mon Jul 09 2012 * @date last modification: Tue Feb 20 2018 * * @brief Exponential irreversible cohesive law of mixed mode loading * * @section LICENSE * * Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "material_cohesive_exponential.hh" #include "dof_synchronizer.hh" #include "solid_mechanics_model.hh" #include "sparse_matrix.hh" namespace akantu { /* -------------------------------------------------------------------------- */ template MaterialCohesiveExponential::MaterialCohesiveExponential( SolidMechanicsModel & model, const ID & id) : MaterialCohesive(model, id) { AKANTU_DEBUG_IN(); this->registerParam("beta", beta, Real(0.), _pat_parsable, "Beta parameter"); this->registerParam("exponential_penalty", exp_penalty, true, _pat_parsable, "Is contact penalty following the exponential law?"); this->registerParam( "contact_tangent", contact_tangent, Real(1.0), _pat_parsable, "Ratio of contact tangent over the initial exponential tangent"); // this->initInternalArray(delta_max, 1, _ek_cohesive); use_previous_delta_max = true; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveExponential::initMaterial() { AKANTU_DEBUG_IN(); MaterialCohesive::initMaterial(); if ((exp_penalty) && (contact_tangent != 1)) { contact_tangent = 1; AKANTU_DEBUG_WARNING("The parsed paramter is forced to " "1.0 when the contact penalty follows the exponential " "law"); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveExponential::computeTraction( const Array & normal, ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); /// define iterators auto traction_it = tractions(el_type, ghost_type).begin(spatial_dimension); auto opening_it = opening(el_type, ghost_type).begin(spatial_dimension); auto normal_it = normal.begin(spatial_dimension); auto traction_end = tractions(el_type, ghost_type).end(spatial_dimension); auto delta_max_it = delta_max(el_type, ghost_type).begin(); auto delta_max_prev_it = delta_max.previous(el_type, ghost_type).begin(); /// compute scalars Real beta2 = beta * beta; /// loop on each quadrature point for (; traction_it != traction_end; ++traction_it, ++opening_it, ++normal_it, ++delta_max_it, ++delta_max_prev_it) { /// compute normal and tangential opening vectors Real normal_opening_norm = opening_it->dot(*normal_it); Vector normal_opening(spatial_dimension); normal_opening = (*normal_it); normal_opening *= normal_opening_norm; Vector tangential_opening(spatial_dimension); tangential_opening = *opening_it; tangential_opening -= normal_opening; Real tangential_opening_norm = tangential_opening.norm(); /** * compute effective opening displacement * @f$ \delta = \sqrt{ * \beta^2 \Delta_t^2 + \Delta_n^2 } @f$ */ Real delta = tangential_opening_norm; delta *= delta * beta2; delta += normal_opening_norm * normal_opening_norm; delta = sqrt(delta); if ((normal_opening_norm < 0) && (std::abs(normal_opening_norm) > Math::getTolerance())) { Vector op_n(*normal_it); op_n *= normal_opening_norm; Vector delta_s(*opening_it); delta_s -= op_n; delta = tangential_opening_norm * beta; computeCoupledTraction(*traction_it, *normal_it, delta, delta_s, *delta_max_it, *delta_max_prev_it); computeCompressiveTraction(*traction_it, *normal_it, normal_opening_norm, *opening_it); } else computeCoupledTraction(*traction_it, *normal_it, delta, *opening_it, *delta_max_it, *delta_max_prev_it); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveExponential::computeCoupledTraction( Vector & tract, const Vector & normal, Real delta, const Vector & opening, Real & delta_max_new, Real delta_max) { AKANTU_DEBUG_IN(); /// full damage case if (std::abs(delta) < Math::getTolerance()) { /// set traction to zero tract.clear(); } else { /// element not fully damaged /** * Compute traction loading @f$ \mathbf{T} = * e \sigma_c \frac{\delta}{\delta_c} e^{-\delta/ \delta_c}@f$ */ /** * Compute traction unloading @f$ \mathbf{T} = * \frac{t_{max}}{\delta_{max}} \delta @f$ */ Real beta2 = beta * beta; Real normal_open_norm = opening.dot(normal); Vector op_n_n(spatial_dimension); op_n_n = normal; op_n_n *= (1 - beta2); op_n_n *= normal_open_norm; tract = beta2 * opening; tract += op_n_n; delta_max_new = std::max(delta_max, delta); - tract *= exp(1) * sigma_c * exp(-delta_max_new / delta_c) / delta_c; + tract *= std::exp(1.) * sigma_c * std::exp(-delta_max_new / delta_c) / delta_c; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveExponential::computeCompressiveTraction( Vector & tract, const Vector & normal, Real delta_n, __attribute__((unused)) const Vector & opening) { Vector temp_tract(normal); if (exp_penalty) { temp_tract *= - delta_n * exp(1) * sigma_c * exp(-delta_n / delta_c) / delta_c; + delta_n * std::exp(1) * sigma_c * std::exp(-delta_n / delta_c) / delta_c; } else { - Real initial_tg = contact_tangent * exp(1) * sigma_c * delta_n / delta_c; + Real initial_tg = contact_tangent * std::exp(1.) * sigma_c * delta_n / delta_c; temp_tract *= initial_tg; } tract += temp_tract; } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveExponential::computeTangentTraction( const ElementType & el_type, Array & tangent_matrix, const Array & normal, GhostType ghost_type) { AKANTU_DEBUG_IN(); auto tangent_it = tangent_matrix.begin(spatial_dimension, spatial_dimension); auto tangent_end = tangent_matrix.end(spatial_dimension, spatial_dimension); auto normal_it = normal.begin(spatial_dimension); auto opening_it = opening(el_type, ghost_type).begin(spatial_dimension); auto delta_max_it = delta_max.previous(el_type, ghost_type).begin(); Real beta2 = beta * beta; /** * compute tangent matrix @f$ \frac{\partial \mathbf{t}} * {\partial \delta} = \hat{\mathbf{t}} \otimes * \frac{\partial (t/\delta)}{\partial \delta} * \frac{\hat{\mathbf{t}}}{\delta}+ \frac{t}{\delta} [ \beta^2 \mathbf{I} + * (1-\beta^2) (\mathbf{n} \otimes \mathbf{n})] @f$ **/ /** * In which @f$ * \frac{\partial(t/ \delta)}{\partial \delta} = * \left\{\begin{array} {l l} * -e \frac{\sigma_c}{\delta_c^2 }e^{-\delta / \delta_c} & \quad if * \delta \geq \delta_{max} \\ * 0 & \quad if \delta < \delta_{max}, \delta_n > 0 * \end{array}\right. @f$ **/ for (; tangent_it != tangent_end; ++tangent_it, ++normal_it, ++opening_it, ++delta_max_it) { Real normal_opening_norm = opening_it->dot(*normal_it); Vector normal_opening(spatial_dimension); normal_opening = (*normal_it); normal_opening *= normal_opening_norm; Vector tangential_opening(spatial_dimension); tangential_opening = *opening_it; tangential_opening -= normal_opening; Real tangential_opening_norm = tangential_opening.norm(); Real delta = tangential_opening_norm; delta *= delta * beta2; delta += normal_opening_norm * normal_opening_norm; delta = sqrt(delta); if ((normal_opening_norm < 0) && (std::abs(normal_opening_norm) > Math::getTolerance())) { Vector op_n(*normal_it); op_n *= normal_opening_norm; Vector delta_s(*opening_it); delta_s -= op_n; delta = tangential_opening_norm * beta; computeCoupledTangent(*tangent_it, *normal_it, delta, delta_s, *delta_max_it); computeCompressivePenalty(*tangent_it, *normal_it, normal_opening_norm); } else computeCoupledTangent(*tangent_it, *normal_it, delta, *opening_it, *delta_max_it); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveExponential::computeCoupledTangent( Matrix & tangent, const Vector & normal, Real delta, const Vector & opening, Real) { AKANTU_DEBUG_IN(); Real beta2 = beta * beta; Matrix J(spatial_dimension, spatial_dimension); J.eye(beta2); if (std::abs(delta) < Math::getTolerance()) { delta = Math::getTolerance(); } Real opening_normal; opening_normal = opening.dot(normal); Vector delta_e(normal); delta_e *= opening_normal; - delta_e *= (1 - beta2); + delta_e *= (1. - beta2); delta_e += (beta2 * opening); - Real exponent = exp(1 - delta / delta_c) * sigma_c / delta_c; + Real exponent = std::exp(1. - delta / delta_c) * sigma_c / delta_c; Matrix first_term(spatial_dimension, spatial_dimension); first_term.outerProduct(normal, normal); - first_term *= (1 - beta2); + first_term *= (1. - beta2); first_term += J; Matrix second_term(spatial_dimension, spatial_dimension); second_term.outerProduct(delta_e, delta_e); second_term /= delta; second_term /= delta_c; Matrix diff(first_term); diff -= second_term; tangent = diff; tangent *= exponent; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveExponential::computeCompressivePenalty( Matrix & tangent, const Vector & normal, Real delta_n) { if (!exp_penalty) - delta_n = 0; + delta_n = 0.; Matrix n_outer_n(spatial_dimension, spatial_dimension); n_outer_n.outerProduct(normal, normal); - Real normal_tg = contact_tangent * exp(1) * sigma_c * - exp(-delta_n / delta_c) * (1 - delta_n / delta_c) / delta_c; + Real normal_tg = contact_tangent * std::exp(1.) * sigma_c * + std::exp(-delta_n / delta_c) * (1. - delta_n / delta_c) / delta_c; n_outer_n *= normal_tg; tangent += n_outer_n; } INSTANTIATE_MATERIAL(cohesive_exponential, MaterialCohesiveExponential); } // 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 b388c4a8a..097491d1a 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,706 +1,704 @@ /** * @file solid_mechanics_model_cohesive.cc * * @author Fabian Barras * @author Mauro Corrado * @author Nicolas Richart * @author Marco Vocialta * * @date creation: Tue May 08 2012 * @date last modification: Wed Feb 21 2018 * * @brief Solid mechanics model for cohesive elements * * @section LICENSE * * Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "solid_mechanics_model_cohesive.hh" #include "aka_iterators.hh" #include "cohesive_element_inserter.hh" #include "element_synchronizer.hh" #include "facet_synchronizer.hh" #include "fe_engine_template.hh" #include "global_ids_updater.hh" #include "integrator_gauss.hh" #include "material_cohesive.hh" #include "mesh_accessor.hh" #include "mesh_global_data_updater.hh" #include "parser.hh" #include "shape_cohesive.hh" /* -------------------------------------------------------------------------- */ #include "dumpable_inline_impl.hh" #ifdef AKANTU_USE_IOHELPER #include "dumper_iohelper_paraview.hh" #endif /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ namespace akantu { class CohesiveMeshGlobalDataUpdater : public MeshGlobalDataUpdater { public: CohesiveMeshGlobalDataUpdater(SolidMechanicsModelCohesive & model) : model(model), mesh(model.getMesh()), global_ids_updater(model.getMesh(), *model.cohesive_synchronizer) {} /* ------------------------------------------------------------------------ */ std::tuple updateData(NewNodesEvent & nodes_event, NewElementsEvent & elements_event) override { auto cohesive_nodes_event = dynamic_cast(&nodes_event); if (not cohesive_nodes_event) return std::make_tuple(nodes_event.getList().size(), elements_event.getList().size()); /// update nodes type auto & new_nodes = cohesive_nodes_event->getList(); auto & old_nodes = cohesive_nodes_event->getOldNodesList(); auto local_nb_new_nodes = new_nodes.size(); auto nb_new_nodes = local_nb_new_nodes; if (mesh.isDistributed()) { MeshAccessor mesh_accessor(mesh); auto & nodes_type = mesh_accessor.getNodesType(); UInt nb_old_nodes = nodes_type.size(); nodes_type.resize(nb_old_nodes + local_nb_new_nodes); for (auto && data : zip(old_nodes, new_nodes)) { UInt old_node, new_node; std::tie(old_node, new_node) = data; nodes_type(new_node) = nodes_type(old_node); } model.updateCohesiveSynchronizers(); nb_new_nodes = global_ids_updater.updateGlobalIDs(new_nodes.size()); } Vector nb_new_stuff = {nb_new_nodes, elements_event.getList().size()}; const auto & comm = mesh.getCommunicator(); comm.allReduce(nb_new_stuff, SynchronizerOperation::_sum); if (nb_new_stuff(1) > 0) { mesh.sendEvent(elements_event); MeshUtils::resetFacetToDouble(mesh.getMeshFacets()); } if (nb_new_nodes > 0) { mesh.sendEvent(nodes_event); // mesh.sendEvent(global_ids_updater.getChangedNodeEvent()); } return std::make_tuple(nb_new_nodes, nb_new_stuff(1)); } private: SolidMechanicsModelCohesive & model; Mesh & mesh; GlobalIdsUpdater global_ids_updater; }; /* -------------------------------------------------------------------------- */ SolidMechanicsModelCohesive::SolidMechanicsModelCohesive( Mesh & mesh, UInt dim, const ID & id, const MemoryID & memory_id) : SolidMechanicsModel(mesh, dim, id, memory_id, ModelType::_solid_mechanics_model_cohesive), tangents("tangents", id), facet_stress("facet_stress", id), facet_material("facet_material", id) { AKANTU_DEBUG_IN(); registerFEEngineObject("CohesiveFEEngine", mesh, Model::spatial_dimension); auto && tmp_material_selector = std::make_shared(*this); tmp_material_selector->setFallback(this->material_selector); this->material_selector = tmp_material_selector; #if defined(AKANTU_USE_IOHELPER) this->mesh.registerDumper("cohesive elements", id); this->mesh.addDumpMeshToDumper("cohesive elements", mesh, Model::spatial_dimension, _not_ghost, _ek_cohesive); #endif if (this->mesh.isDistributed()) { /// create the distributed synchronizer for cohesive elements this->cohesive_synchronizer = std::make_unique( mesh, "cohesive_distributed_synchronizer"); auto & synchronizer = mesh.getElementSynchronizer(); this->cohesive_synchronizer->split(synchronizer, [](auto && el) { return Mesh::getKind(el.type) == _ek_cohesive; }); this->registerSynchronizer(*cohesive_synchronizer, _gst_material_id); this->registerSynchronizer(*cohesive_synchronizer, _gst_smm_stress); this->registerSynchronizer(*cohesive_synchronizer, _gst_smm_boundary); } this->inserter = std::make_unique( this->mesh, id + ":cohesive_element_inserter"); registerFEEngineObject( "FacetsFEEngine", mesh.getMeshFacets(), Model::spatial_dimension - 1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ SolidMechanicsModelCohesive::~SolidMechanicsModelCohesive() = default; /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::setTimeStep(Real time_step, const ID & solver_id) { SolidMechanicsModel::setTimeStep(time_step, solver_id); #if defined(AKANTU_USE_IOHELPER) this->mesh.getDumper("cohesive elements").setTimeStep(time_step); #endif } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::initFullImpl(const ModelOptions & options) { AKANTU_DEBUG_IN(); const auto & smmc_options = dynamic_cast(options); this->is_extrinsic = smmc_options.extrinsic; inserter->setIsExtrinsic(is_extrinsic); 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( synchronizer, id + ":facet_stress_synchronizer"); facet_stress_synchronizer->swapSendRecv(); this->registerSynchronizer(*facet_stress_synchronizer, _gst_smmc_facets_stress); } } MeshAccessor mesh_accessor(mesh); mesh_accessor.registerGlobalDataUpdater( std::make_unique(*this)); ParserSection section; bool is_empty; std::tie(section, is_empty) = this->getParserSection(); if (not is_empty) { auto inserter_section = section.getSubSections(ParserType::_cohesive_inserter); if (inserter_section.begin() != inserter_section.end()) { inserter->parseSection(*inserter_section.begin()); } } SolidMechanicsModel::initFullImpl(options); AKANTU_DEBUG_OUT(); } // namespace akantu /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::initMaterials() { AKANTU_DEBUG_IN(); // make sure the material are instantiated if (!are_materials_instantiated) instantiateMaterials(); /// find the first cohesive material UInt cohesive_index = UInt(-1); for (auto && material : enumerate(materials)) { if (dynamic_cast(std::get<1>(material).get())) { cohesive_index = std::get<0>(material); break; } } if (cohesive_index == UInt(-1)) AKANTU_EXCEPTION("No cohesive materials in the material input file"); material_selector->setFallback(cohesive_index); // set the facet information in the material in case of dynamic insertion // to know what material to call for stress checks const Mesh & mesh_facets = inserter->getMeshFacets(); facet_material.initialize( mesh_facets, _spatial_dimension = spatial_dimension - 1, _with_nb_element = true, _default_value = material_selector->getFallbackValue()); for_each_element( mesh_facets, [&](auto && element) { auto mat_index = (*material_selector)(element); auto & mat = dynamic_cast(*materials[mat_index]); facet_material(element) = mat_index; if (is_extrinsic) { mat.addFacet(element); } }, _spatial_dimension = spatial_dimension - 1, _ghost_type = _not_ghost); SolidMechanicsModel::initMaterials(); if (is_extrinsic) { this->initAutomaticInsertion(); } else { this->insertIntrinsicElements(); } AKANTU_DEBUG_OUT(); } // namespace akantu /* -------------------------------------------------------------------------- */ /** * Initialize the model,basically it pre-compute the shapes, shapes derivatives * and jacobian */ void SolidMechanicsModelCohesive::initModel() { AKANTU_DEBUG_IN(); SolidMechanicsModel::initModel(); /// add cohesive type connectivity ElementType type = _not_defined; for (auto && type_ghost : ghost_types) { for (const auto & tmp_type : mesh.elementTypes(spatial_dimension, type_ghost)) { const auto & connectivity = mesh.getConnectivity(tmp_type, type_ghost); if (connectivity.size() == 0) continue; type = tmp_type; auto type_facet = Mesh::getFacetType(type); auto type_cohesive = FEEngine::getCohesiveElementType(type_facet); mesh.addConnectivityType(type_cohesive, type_ghost); } } AKANTU_DEBUG_ASSERT(type != _not_defined, "No elements in the mesh"); getFEEngine("CohesiveFEEngine").initShapeFunctions(_not_ghost); getFEEngine("CohesiveFEEngine").initShapeFunctions(_ghost); if (is_extrinsic) { getFEEngine("FacetsFEEngine").initShapeFunctions(_not_ghost); getFEEngine("FacetsFEEngine").initShapeFunctions(_ghost); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::insertIntrinsicElements() { AKANTU_DEBUG_IN(); inserter->insertIntrinsicElements(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::initAutomaticInsertion() { AKANTU_DEBUG_IN(); this->inserter->limitCheckFacets(); this->updateFacetStressSynchronizer(); this->resizeFacetStress(); /// compute normals on facets this->computeNormals(); this->initStressInterpolation(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::updateAutomaticInsertion() { AKANTU_DEBUG_IN(); this->inserter->limitCheckFacets(); this->updateFacetStressSynchronizer(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::initStressInterpolation() { Mesh & mesh_facets = inserter->getMeshFacets(); /// compute quadrature points coordinates on facets Array & position = mesh.getNodes(); ElementTypeMapArray quad_facets("quad_facets", id); quad_facets.initialize(mesh_facets, _nb_component = Model::spatial_dimension, _spatial_dimension = Model::spatial_dimension - 1); // mesh_facets.initElementTypeMapArray(quad_facets, Model::spatial_dimension, // Model::spatial_dimension - 1); getFEEngine("FacetsFEEngine") .interpolateOnIntegrationPoints(position, quad_facets); /// compute elements quadrature point positions and build /// element-facet quadrature points data structure ElementTypeMapArray elements_quad_facets("elements_quad_facets", id); elements_quad_facets.initialize( mesh, _nb_component = Model::spatial_dimension, _spatial_dimension = Model::spatial_dimension); // mesh.initElementTypeMapArray(elements_quad_facets, // Model::spatial_dimension, // Model::spatial_dimension); for (auto elem_gt : ghost_types) { for (auto & type : mesh.elementTypes(Model::spatial_dimension, elem_gt)) { UInt nb_element = mesh.getNbElement(type, elem_gt); if (nb_element == 0) continue; /// compute elements' quadrature points and list of facet /// quadrature points positions by element const auto & facet_to_element = mesh_facets.getSubelementToElement(type, elem_gt); auto & el_q_facet = elements_quad_facets(type, elem_gt); auto facet_type = Mesh::getFacetType(type); auto nb_quad_per_facet = getFEEngine("FacetsFEEngine").getNbIntegrationPoints(facet_type); auto nb_facet_per_elem = facet_to_element.getNbComponent(); // small hack in the loop to skip boundary elements, they are silently // initialized to NaN to see if this causes problems el_q_facet.resize(nb_element * nb_facet_per_elem * nb_quad_per_facet, std::numeric_limits::quiet_NaN()); for (auto && data : zip(make_view(facet_to_element), make_view(el_q_facet, spatial_dimension, nb_quad_per_facet))) { const auto & global_facet = std::get<0>(data); auto & el_q = std::get<1>(data); if (global_facet == ElementNull) continue; Matrix quad_f = make_view(quad_facets(global_facet.type, global_facet.ghost_type), spatial_dimension, nb_quad_per_facet) .begin()[global_facet.element]; el_q = quad_f; // for (UInt q = 0; q < nb_quad_per_facet; ++q) { // for (UInt s = 0; s < Model::spatial_dimension; ++s) { // el_q_facet(el * nb_facet_per_elem * nb_quad_per_facet + // f * nb_quad_per_facet + q, // s) = quad_f(global_facet * nb_quad_per_facet + q, // s); // } // } //} } } } /// loop over non cohesive materials for (auto && material : materials) { if (dynamic_cast(material.get())) continue; /// initialize the interpolation function material->initElementalFieldInterpolation(elements_quad_facets); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::assembleInternalForces() { AKANTU_DEBUG_IN(); // f_int += f_int_cohe for (auto & material : this->materials) { try { auto & mat = dynamic_cast(*material); mat.computeTraction(_not_ghost); } catch (std::bad_cast & bce) { } } SolidMechanicsModel::assembleInternalForces(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::computeNormals() { AKANTU_DEBUG_IN(); Mesh & mesh_facets = this->inserter->getMeshFacets(); this->getFEEngine("FacetsFEEngine") .computeNormalsOnIntegrationPoints(_not_ghost); /** * @todo store tangents while computing normals instead of * recomputing them as follows: */ /* ------------------------------------------------------------------------ */ UInt tangent_components = Model::spatial_dimension * (Model::spatial_dimension - 1); tangents.initialize(mesh_facets, _nb_component = tangent_components, _spatial_dimension = Model::spatial_dimension - 1); // mesh_facets.initElementTypeMapArray(tangents, tangent_components, // Model::spatial_dimension - 1); for (auto facet_type : mesh_facets.elementTypes(Model::spatial_dimension - 1)) { const Array & normals = this->getFEEngine("FacetsFEEngine") .getNormalsOnIntegrationPoints(facet_type); Array & tangents = this->tangents(facet_type); Math::compute_tangents(normals, tangents); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::interpolateStress() { ElementTypeMapArray by_elem_result("temporary_stress_by_facets", id); for (auto & material : materials) { auto * mat = dynamic_cast(material.get()); if (mat == nullptr) /// interpolate stress on facet quadrature points positions material->interpolateStressOnFacets(facet_stress, by_elem_result); } this->synchronize(_gst_smmc_facets_stress); } /* -------------------------------------------------------------------------- */ UInt SolidMechanicsModelCohesive::checkCohesiveStress() { AKANTU_DEBUG_IN(); if (not is_extrinsic) { AKANTU_EXCEPTION( "This function can only be used for extrinsic cohesive elements"); } interpolateStress(); for (auto & mat : materials) { auto * mat_cohesive = dynamic_cast(mat.get()); if (mat_cohesive) { /// check which not ghost cohesive elements are to be created mat_cohesive->checkInsertion(); } } /// communicate data among processors // this->synchronize(_gst_smmc_facets); /// insert cohesive elements UInt nb_new_elements = inserter->insertElements(); // if (nb_new_elements > 0) { // this->reinitializeSolver(); // } AKANTU_DEBUG_OUT(); return nb_new_elements; } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::onElementsAdded( const Array & element_list, const NewElementsEvent & event) { AKANTU_DEBUG_IN(); SolidMechanicsModel::onElementsAdded(element_list, event); if (is_extrinsic) resizeFacetStress(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::onNodesAdded(const Array & new_nodes, const NewNodesEvent & event) { AKANTU_DEBUG_IN(); SolidMechanicsModel::onNodesAdded(new_nodes, event); const CohesiveNewNodesEvent * cohesive_event; if ((cohesive_event = dynamic_cast(&event)) == nullptr) return; const auto & old_nodes = cohesive_event->getOldNodesList(); auto copy = [this, &new_nodes, &old_nodes](auto & arr) { UInt new_node, old_node; auto view = make_view(arr, spatial_dimension); auto begin = view.begin(); for (auto && pair : zip(new_nodes, old_nodes)) { std::tie(new_node, old_node) = pair; auto old_ = begin + old_node; auto new_ = begin + new_node; *new_ = *old_; } }; copy(*displacement); copy(*blocked_dofs); if (velocity) copy(*velocity); if (acceleration) copy(*acceleration); if (current_position) copy(*current_position); if (previous_displacement) copy(*previous_displacement); // if (external_force) // copy(*external_force); // if (internal_force) // copy(*internal_force); if (displacement_increment) copy(*displacement_increment); copy(getDOFManager().getSolution("displacement")); // this->assembleMassLumped(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::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(); this->facet_stress.initialize(getFEEngine("FacetsFEEngine"), _nb_component = 2 * spatial_dimension * spatial_dimension, _spatial_dimension = spatial_dimension - 1); // for (auto && ghost_type : ghost_types) { // for (const auto & type : // mesh_facets.elementTypes(spatial_dimension - 1, ghost_type)) { // UInt nb_facet = mesh_facets.getNbElement(type, ghost_type); // UInt nb_quadrature_points = getFEEngine("FacetsFEEngine") // .getNbIntegrationPoints(type, // ghost_type); // UInt new_size = nb_facet * nb_quadrature_points; // facet_stress(type, ghost_type).resize(new_size); // } // } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::addDumpGroupFieldToDumper( const std::string & dumper_name, const std::string & field_id, const std::string & group_name, const ElementKind & element_kind, bool padding_flag) { AKANTU_DEBUG_IN(); UInt spatial_dimension = Model::spatial_dimension; ElementKind _element_kind = element_kind; if (dumper_name == "cohesive elements") { _element_kind = _ek_cohesive; } else if (dumper_name == "facets") { spatial_dimension = Model::spatial_dimension - 1; } SolidMechanicsModel::addDumpGroupFieldToDumper(dumper_name, field_id, group_name, spatial_dimension, _element_kind, padding_flag); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::onDump() { this->flattenAllRegisteredInternals(_ek_cohesive); SolidMechanicsModel::onDump(); } /* -------------------------------------------------------------------------- */ } // namespace akantu diff --git a/src/synchronizer/dof_synchronizer.cc b/src/synchronizer/dof_synchronizer.cc index ed0a9bcf2..515d44bac 100644 --- a/src/synchronizer/dof_synchronizer.cc +++ b/src/synchronizer/dof_synchronizer.cc @@ -1,279 +1,284 @@ /** * @file dof_synchronizer.cc * * @author Aurelia Isabel Cuba Ramos * @author Nicolas Richart * * @date creation: Fri Jun 17 2011 * @date last modification: Tue Feb 06 2018 * * @brief DOF synchronizing object implementation * * @section LICENSE * * Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "dof_synchronizer.hh" #include "aka_iterators.hh" #include "dof_manager_default.hh" #include "mesh.hh" #include "node_synchronizer.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ /** * A DOFSynchronizer needs a mesh and the number of degrees of freedom * per node to be created. In the constructor computes the local and global dof * number for each dof. The member * proc_informations (std vector) is resized with the number of mpi * processes. Each entry in the vector is a PerProcInformations object * that contains the interactions of the current mpi process (prank) with the * mpi process corresponding to the position of that entry. Every * ProcInformations object contains one array with the dofs that have * to be sent to prank and a second one with dofs that willl be received form * prank. * This information is needed for the asychronous communications. The * constructor sets up this information. */ DOFSynchronizer::DOFSynchronizer(DOFManagerDefault & dof_manager, const ID & id, MemoryID memory_id) : SynchronizerImpl(dof_manager.getCommunicator(), id, memory_id), root(0), dof_manager(dof_manager), root_dofs(0, 1, "dofs-to-receive-from-master"), dof_changed(true) { std::vector dof_ids = dof_manager.getDOFIDs(); // Transfers nodes to global equation numbers in new schemes for (const ID & dof_id : dof_ids) { registerDOFs(dof_id); } this->initScatterGatherCommunicationScheme(); } /* -------------------------------------------------------------------------- */ DOFSynchronizer::~DOFSynchronizer() = default; /* -------------------------------------------------------------------------- */ void DOFSynchronizer::registerDOFs(const ID & dof_id) { if (this->nb_proc == 1) return; if (dof_manager.getSupportType(dof_id) != _dst_nodal) return; using const_scheme_iterator = Communications::const_scheme_iterator; const auto equation_numbers = dof_manager.getLocalEquationNumbers(dof_id); const auto & associated_nodes = dof_manager.getDOFsAssociatedNodes(dof_id); const auto & node_synchronizer = dof_manager.getMesh().getNodeSynchronizer(); const auto & node_communications = node_synchronizer.getCommunications(); auto transcode_node_to_global_dof_scheme = [this, &associated_nodes, &equation_numbers](const_scheme_iterator it, const_scheme_iterator end, const CommunicationSendRecv & sr) -> void { for (; it != end; ++it) { auto & scheme = communications.createScheme(it->first, sr); const auto & node_scheme = it->second; for (auto & node : node_scheme) { auto an_begin = associated_nodes.begin(); auto an_it = an_begin; auto an_end = associated_nodes.end(); std::vector global_dofs_per_node; while ((an_it = std::find(an_it, an_end, node)) != an_end) { UInt pos = an_it - an_begin; UInt local_eq_num = equation_numbers(pos); UInt global_eq_num = dof_manager.localToGlobalEquationNumber(local_eq_num); global_dofs_per_node.push_back(global_eq_num); ++an_it; } std::sort(global_dofs_per_node.begin(), global_dofs_per_node.end()); std::transform(global_dofs_per_node.begin(), global_dofs_per_node.end(), global_dofs_per_node.begin(), [this](UInt g) -> UInt { UInt l = dof_manager.globalToLocalEquationNumber(g); return l; }); for (auto & leqnum : global_dofs_per_node) { scheme.push_back(leqnum); } } } }; for (auto sr_it = send_recv_t::begin(); sr_it != send_recv_t::end(); ++sr_it) { auto ncs_it = node_communications.begin_scheme(*sr_it); auto ncs_end = node_communications.end_scheme(*sr_it); transcode_node_to_global_dof_scheme(ncs_it, ncs_end, *sr_it); } dof_changed = true; } /* -------------------------------------------------------------------------- */ void DOFSynchronizer::initScatterGatherCommunicationScheme() { AKANTU_DEBUG_IN(); if (this->nb_proc == 1) { dof_changed = false; AKANTU_DEBUG_OUT(); return; } UInt nb_dofs = dof_manager.getLocalSystemSize(); this->root_dofs.clear(); this->master_receive_dofs.clear(); Array dofs_to_send; for (UInt n = 0; n < nb_dofs; ++n) { if (dof_manager.isLocalOrMasterDOF(n)) { auto & receive_per_proc = master_receive_dofs[this->root]; UInt global_dof = dof_manager.localToGlobalEquationNumber(n); root_dofs.push_back(n); receive_per_proc.push_back(global_dof); dofs_to_send.push_back(global_dof); } } if (this->rank == UInt(this->root)) { Array nb_dof_per_proc(this->nb_proc); communicator.gather(dofs_to_send.size(), nb_dof_per_proc); std::vector requests; for (UInt p = 0; p < nb_proc; ++p) { if (p == UInt(this->root)) continue; auto & receive_per_proc = master_receive_dofs[p]; receive_per_proc.resize(nb_dof_per_proc(p)); if (nb_dof_per_proc(p) != 0) requests.push_back(communicator.asyncReceive( receive_per_proc, p, Tag::genTag(p, 0, Tag::_GATHER_INITIALIZATION, this->hash_id))); } communicator.waitAll(requests); communicator.freeCommunicationRequest(requests); } else { communicator.gather(dofs_to_send.size(), this->root); AKANTU_DEBUG(dblDebug, "I have " << nb_dofs << " dofs (" << dofs_to_send.size() << " to send to master proc"); if (dofs_to_send.size() != 0) communicator.send(dofs_to_send, this->root, Tag::genTag(this->rank, 0, Tag::_GATHER_INITIALIZATION, this->hash_id)); } dof_changed = false; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ bool DOFSynchronizer::hasChanged() { communicator.allReduce(dof_changed, SynchronizerOperation::_lor); return dof_changed; } /* -------------------------------------------------------------------------- */ -void DOFSynchronizer::onNodesAdded(const Array & nodes_list) { +void DOFSynchronizer::onNodesAdded(const Array & /*nodes_list*/) { auto dof_ids = dof_manager.getDOFIDs(); - const auto & node_synchronizer = dof_manager.getMesh().getNodeSynchronizer(); - const auto & node_communications = node_synchronizer.getCommunications(); - - std::set relevant_nodes; - std::map> nodes_per_proc[2]; - for (auto sr_it = send_recv_t::begin(); sr_it != send_recv_t::end(); - ++sr_it) { - auto sit = node_communications.begin_scheme(*sr_it); - auto send = node_communications.end_scheme(*sr_it); - - for (; sit != send; ++sit) { - auto proc = sit->first; - const auto & scheme = sit->second; - for (auto node : nodes_list) { - if (scheme.find(node) == UInt(-1)) - continue; - relevant_nodes.insert(node); - nodes_per_proc[*sr_it][proc].push_back(node); - } + for (auto sr : iterate_send_recv) { + for (auto && data : communications.iterateSchemes(sr)) { + auto & scheme = data.second; + scheme.resize(0); } } - std::map> dofs_per_proc[2]; - for (auto & dof_id : dof_ids) { - const auto & associated_nodes = dof_manager.getDOFsAssociatedNodes(dof_id); - const auto & local_equation_numbers = - dof_manager.getEquationsNumbers(dof_id); - - for (auto tuple : zip(associated_nodes, local_equation_numbers)) { - UInt assoc_node; - UInt local_eq_num; - std::tie(assoc_node, local_eq_num) = tuple; - - for (auto sr_it = send_recv_t::begin(); sr_it != send_recv_t::end(); - ++sr_it) { - for (auto & pair : nodes_per_proc[*sr_it]) { - if (std::find(pair.second.end(), pair.second.end(), assoc_node) != - pair.second.end()) { - dofs_per_proc[*sr_it][pair.first].push_back(local_eq_num); - } - } - } - } + for(auto & dof_id : dof_ids) { + registerDOFs(dof_id); } - for (auto sr_it = send_recv_t::begin(); sr_it != send_recv_t::end(); - ++sr_it) { - for (auto & pair : dofs_per_proc[*sr_it]) { - std::sort(pair.second.begin(), pair.second.end(), - [this](UInt la, UInt lb) -> bool { - UInt ga = dof_manager.localToGlobalEquationNumber(la); - UInt gb = dof_manager.localToGlobalEquationNumber(lb); - return ga < gb; - }); - - auto & scheme = communications.getScheme(pair.first, *sr_it); - for (auto leq : pair.second) { - scheme.push_back(leq); - } - } - } + // const auto & node_synchronizer = dof_manager.getMesh().getNodeSynchronizer(); + // const auto & node_communications = node_synchronizer.getCommunications(); + + // std::map> nodes_per_proc[2]; + + // for (auto sr : iterate_send_recv) { + // for (auto && data : node_communications.iterateSchemes(sr)) { + // auto proc = data.first; + // const auto & scheme = data.second; + // for (auto node : scheme) { + // nodes_per_proc[sr][proc].push_back(node); + // } + // } + // } + + // std::map> dofs_per_proc[2]; + // for (auto & dof_id : dof_ids) { + // const auto & associated_nodes = dof_manager.getDOFsAssociatedNodes(dof_id); + // const auto & local_equation_numbers = + // dof_manager.getEquationsNumbers(dof_id); + + // for (auto tuple : zip(associated_nodes, local_equation_numbers)) { + // UInt assoc_node; + // UInt local_eq_num; + // std::tie(assoc_node, local_eq_num) = tuple; + + // for (auto sr_it = send_recv_t::begin(); sr_it != send_recv_t::end(); + // ++sr_it) { + // for (auto & pair : nodes_per_proc[*sr_it]) { + // if (std::find(pair.second.end(), pair.second.end(), assoc_node) != + // pair.second.end()) { + // dofs_per_proc[*sr_it][pair.first].push_back(local_eq_num); + // } + // } + // } + // } + // } + + // for (auto sr_it = send_recv_t::begin(); sr_it != send_recv_t::end(); + // ++sr_it) { + // for (auto & pair : dofs_per_proc[*sr_it]) { + // std::sort(pair.second.begin(), pair.second.end(), + // [this](UInt la, UInt lb) -> bool { + // UInt ga = dof_manager.localToGlobalEquationNumber(la); + // UInt gb = dof_manager.localToGlobalEquationNumber(lb); + // return ga < gb; + // }); + + // auto & scheme = communications.getScheme(pair.first, *sr_it); + // scheme.resize(0); + // for (auto leq : pair.second) { + // scheme.push_back(leq); + // } + // } + // } dof_changed = true; } /* -------------------------------------------------------------------------- */ } // namespace akantu diff --git a/src/synchronizer/facet_synchronizer.cc b/src/synchronizer/facet_synchronizer.cc index 123f47dce..ae247444b 100644 --- a/src/synchronizer/facet_synchronizer.cc +++ b/src/synchronizer/facet_synchronizer.cc @@ -1,209 +1,208 @@ /** * @file facet_synchronizer.cc * * @author Nicolas Richart * @author Marco Vocialta * * @date creation: Wed Nov 05 2014 * @date last modification: Fri Jan 26 2018 * * @brief Facet synchronizer for parallel simulations with cohesive elments * * @section LICENSE * * Copyright (©) 2015-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "facet_synchronizer.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ FacetSynchronizer::FacetSynchronizer( Mesh & mesh, const ElementSynchronizer & element_synchronizer, const ID & id, MemoryID memory_id) : ElementSynchronizer(mesh, id, memory_id) { - const auto & comm = mesh.getCommunicator(); auto spatial_dimension = mesh.getSpatialDimension(); element_to_prank.initialize(mesh, _spatial_dimension = spatial_dimension - 1, _ghost_type = _ghost, _with_nb_element = true, _default_value = rank); // Build element to prank for (auto && scheme_pair : element_synchronizer.communications.iterateSchemes(_recv)) { auto proc = std::get<0>(scheme_pair); const auto & scheme = std::get<1>(scheme_pair); for (auto && elem : scheme) { const auto & facet_to_element = mesh.getSubelementToElement(elem.type, elem.ghost_type); Vector facets = facet_to_element.begin( facet_to_element.getNbComponent())[elem.element]; for (UInt f = 0; f < facets.size(); ++f) { const auto & facet = facets(f); if (facet == ElementNull) continue; if (facet.ghost_type == _not_ghost) continue; auto & facet_rank = element_to_prank(facet); if ((proc < UInt(facet_rank)) || (UInt(facet_rank) == rank)) facet_rank = proc; } } } ElementTypeMapArray facet_global_connectivities( "facet_global_connectivities", id, memory_id); facet_global_connectivities.initialize( mesh, _spatial_dimension = spatial_dimension - 1, _with_nb_element = true, _with_nb_nodes_per_element = true); mesh.getGlobalConnectivity(facet_global_connectivities); // \TODO perhaps a global element numbering might be useful here... - for (auto type : facet_global_connectivities.elementTypes( - _spatial_dimension = _all_dimensions, - _element_kind = _ek_not_defined, _ghost_type = _not_ghost)) { + for (auto type : facet_global_connectivities.elementTypes(_spatial_dimension = + _all_dimensions, + _element_kind = _ek_not_defined, _ghost_type = _not_ghost)) { auto & conn = facet_global_connectivities(type, _not_ghost); auto conn_view = make_view(conn, conn.getNbComponent()); std::for_each(conn_view.begin(), conn_view.end(), [&](auto & conn) { std::sort(conn.storage(), conn.storage() + conn.size()); }); } /// init facet check tracking ElementTypeMapArray facet_checked("facet_checked", id, memory_id); std::map> recv_connectivities; /// Generate the recv scheme and connnectivities to send to the other /// processors for (auto && scheme_pair : element_synchronizer.communications.iterateSchemes(_recv)) { facet_checked.initialize(mesh, _spatial_dimension = spatial_dimension - 1, _ghost_type = _ghost, _with_nb_element = true, _default_value = false); auto proc = scheme_pair.first; const auto & elements = scheme_pair.second; auto & facet_scheme = communications.createScheme(proc, _recv); // this creates empty arrays... auto & connectivities_for_proc = recv_connectivities[proc]; connectivities_for_proc.setID( id + ":connectivities_for_proc:" + std::to_string(proc)); connectivities_for_proc.initialize( mesh, _spatial_dimension = spatial_dimension - 1, _with_nb_nodes_per_element = true, _ghost_type = _ghost); // for every element in the element synchronizer communication scheme, // check the facets to see if they should be communicated and create a // connectivity array to match with the one other processors might send for (auto && element : elements) { const auto & facet_to_element = mesh.getSubelementToElement(element.type, element.ghost_type); Vector facets = facet_to_element.begin( facet_to_element.getNbComponent())[element.element]; for (UInt f = 0; f < facets.size(); ++f) { auto & facet = facets(f); // exclude no valid facets if (facet == ElementNull) continue; // exclude _ghost facet from send scheme and _not_ghost from receive if (facet.ghost_type != _ghost) continue; // exclude facet from other processors then the one of current // interest in case of receive scheme if (UInt(element_to_prank(facet)) != proc) continue; auto & checked = facet_checked(facet); // skip already checked facets if (checked) continue; checked = true; facet_scheme.push_back(facet); auto & global_conn = facet_global_connectivities(facet.type, facet.ghost_type); Vector conn = global_conn.begin(global_conn.getNbComponent())[facet.element]; std::sort(conn.storage(), conn.storage() + conn.size()); connectivities_for_proc(facet.type, facet.ghost_type).push_back(conn); } } } std::vector send_requests; /// do every communication by element type for (auto && type : mesh.elementTypes(spatial_dimension - 1)) { for (auto && pair : recv_connectivities) { auto proc = std::get<0>(pair); const auto & connectivities_for_proc = std::get<1>(pair); auto && tag = Tag::genTag(proc, type, 1337); send_requests.push_back( - comm.asyncSend(connectivities_for_proc(type, _ghost), proc, tag, - CommunicationMode::_synchronous)); + communicator.asyncSend(connectivities_for_proc(type, _ghost), proc, + tag, CommunicationMode::_synchronous)); } auto nb_nodes_per_facet = Mesh::getNbNodesPerElement(type); Array buffer; - comm.receiveAnyNumber( + communicator.receiveAnyNumber( send_requests, buffer, [&](auto && proc, auto && message) { auto & local_connectivities = facet_global_connectivities(type, _not_ghost); auto & send_scheme = communications.createScheme(proc, _send); auto conn_view = make_view(local_connectivities, nb_nodes_per_facet); auto conn_begin = conn_view.begin(); auto conn_end = conn_view.end(); for (const auto & c_to_match : make_view(message, nb_nodes_per_facet)) { auto it = std::find(conn_begin, conn_end, c_to_match); if (it != conn_end) { auto facet = Element{type, UInt(it - conn_begin), _not_ghost}; send_scheme.push_back(facet); } else { AKANTU_EXCEPTION("No local facet found to send to proc " << proc << " corresponding to " << c_to_match); } } }, Tag::genTag(rank, type, 1337)); } } } // namespace akantu diff --git a/src/synchronizer/node_synchronizer.cc b/src/synchronizer/node_synchronizer.cc index 15c889ffb..f347846d3 100644 --- a/src/synchronizer/node_synchronizer.cc +++ b/src/synchronizer/node_synchronizer.cc @@ -1,122 +1,161 @@ /** - * @file node_synchronizer.cc - * - * @author Nicolas Richart - * - * @date creation: Fri Jun 18 2010 - * @date last modification: Wed Nov 15 2017 - * - * @brief Implementation of the node synchronizer - * - * @section LICENSE - * - * Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) - * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) - * - * Akantu is free software: you can redistribute it and/or modify it under the - * terms of the GNU Lesser General Public License as published by the Free - * Software Foundation, either version 3 of the License, or (at your option) any - * later version. - * - * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY - * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR - * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more - * details. - * - * You should have received a copy of the GNU Lesser General Public License - * along with Akantu. If not, see . - * - */ +* @file node_synchronizer.cc +* +* @author Nicolas Richart +* +* @date creation: Fri Jun 18 2010 +* @date last modification: Wed Nov 15 2017 +* +* @brief Implementation of the node synchronizer +* +* @section LICENSE +* +* Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) +* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) +* +* Akantu is free software: you can redistribute it and/or modify it under the +* terms of the GNU Lesser General Public License as published by the Free +* Software Foundation, either version 3 of the License, or (at your option) any +* later version. +* +* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY +* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR +* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more +* details. +* +* You should have received a copy of the GNU Lesser General Public License +* along with Akantu. If not, see . +* +*/ /* -------------------------------------------------------------------------- */ #include "node_synchronizer.hh" #include "mesh.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ NodeSynchronizer::NodeSynchronizer(Mesh & mesh, const ID & id, MemoryID memory_id, const bool register_to_event_manager, EventHandlerPriority event_priority) : SynchronizerImpl(mesh.getCommunicator(), id, memory_id), mesh(mesh) { AKANTU_DEBUG_IN(); if (register_to_event_manager) { this->mesh.registerEventHandler(*this, event_priority); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ NodeSynchronizer::~NodeSynchronizer() = default; /* -------------------------------------------------------------------------- */ -void NodeSynchronizer::onNodesAdded(const Array & nodes_list, +void NodeSynchronizer::onNodesAdded(const Array & /*nodes_list*/, const NewNodesEvent &) { std::map> nodes_per_proc; - Array sizes_per_proc(nb_proc, nb_proc, UInt(0)); - Vector local_sizes_per_proc = sizes_per_proc.begin(nb_proc)[rank]; - for (auto & local_id : nodes_list) { + // recreates fully the schemes due to changes of global ids + // \TODO add an event to handle global id changes + for(auto && data : communications.iterateRecvSchemes()) { + auto & scheme = data.second; + scheme.resize(0); + } + + for (auto && local_id : arange(mesh.getNbNodes())) { auto type = mesh.getNodeType(local_id); if (type < 0) continue; // local, master or pure ghost auto global_id = mesh.getNodeGlobalId(local_id); auto proc = UInt(type); nodes_per_proc[proc].push_back(global_id); - ++local_sizes_per_proc[proc]; - auto & scheme = communications.getScheme(proc, _recv); + auto & scheme = communications.createScheme(proc, _recv); scheme.push_back(local_id); } - communicator.allGather(sizes_per_proc); - std::vector send_request_per_proc, - recv_request_per_proc; - std::map> nodes_needed_by_proc; - for (UInt proc = 0; proc < nb_proc; ++proc) { - auto size = sizes_per_proc(proc, rank); - if (size == 0) - continue; - - nodes_needed_by_proc[proc].resize(size); - recv_request_per_proc.push_back(communicator.asyncReceive( - nodes_needed_by_proc[proc], proc, Tag::genTag(rank, 0, 0))); - } - + std::vector send_requests; for (auto & pair : nodes_per_proc) { auto proc = pair.first; auto & nodes = pair.second; - send_request_per_proc.push_back( - communicator.asyncSend(nodes, proc, Tag::genTag(proc, 0, 0))); + send_requests.push_back( + communicator.asyncSend(nodes, proc, Tag::genTag(proc, 0, 0xcafe))); } - UInt req_nb; - while ((req_nb = communicator.waitAny(recv_request_per_proc)) != UInt(-1)) { - auto & request = recv_request_per_proc[req_nb]; - auto proc = request.getSource(); + Array buffer; + + communicator.receiveAnyNumber( + send_requests, buffer, + [&](auto && proc, auto && nodes) { + auto & scheme = communications.createScheme(proc, _send); + scheme.resize(nodes.size()); + for (auto && data : enumerate(nodes)) { + auto global_id = std::get<1>(data); + auto local_id = mesh.getNodeLocalId(global_id); + AKANTU_DEBUG_ASSERT(local_id != UInt(-1), + "The global node " << global_id + << "is not known on rank " + << rank); + scheme[std::get<0>(data)] = local_id; + } + }, + Tag::genTag(rank, 0, 0xcafe)); +} + +/* -------------------------------------------------------------------------- */ +UInt NodeSynchronizer::sanityCheckDataSize(const Array & nodes, + const SynchronizationTag & tag, + bool from_comm_desc) const { + UInt size = + SynchronizerImpl::sanityCheckDataSize(nodes, tag, from_comm_desc); - auto & nodes = nodes_needed_by_proc[proc]; + // positions + size += mesh.getSpatialDimension() * sizeof(Real) * nodes.size(); - auto & scheme = communications.getScheme(proc, _send); - for (auto global_id : nodes) { - auto local_id = mesh.getNodeLocalId(global_id); - AKANTU_DEBUG_ASSERT(local_id != UInt(-1), - "The global node " - << global_id << "is not known on rank " << rank); - scheme.push_back(local_id); - } - recv_request_per_proc.erase(recv_request_per_proc.begin() + req_nb); + return size; +} + +/* -------------------------------------------------------------------------- */ +void NodeSynchronizer::packSanityCheckData( + CommunicationBuffer & buffer, const Array & nodes, + const SynchronizationTag & /*tag*/) const { + auto dim = mesh.getSpatialDimension(); + for (auto && node : nodes) { + buffer << Vector(mesh.getNodes().begin(dim)[node]); } +} - communicator.waitAll(send_request_per_proc); - communicator.freeCommunicationRequest(send_request_per_proc); +/* -------------------------------------------------------------------------- */ +void NodeSynchronizer::unpackSanityCheckData(CommunicationBuffer & buffer, + const Array & nodes, + const SynchronizationTag & tag, + UInt proc, UInt rank) const { + auto dim = mesh.getSpatialDimension(); + // std::set skip_conn_tags{_gst_smmc_facets_conn, + // _gst_giu_global_conn}; + + // bool is_skip_tag_conn = skip_conn_tags.find(tag) != skip_conn_tags.end(); + + for (auto && node : nodes) { + Vector pos_remote(dim); + buffer >> pos_remote; + Vector pos(mesh.getNodes().begin(dim)[node]); + + auto dist = pos_remote.distance(pos); + if (not Math::are_float_equal(dist, 0.)) { + AKANTU_EXCEPTION("Unpacking an unknown value for the node " + << node << "(position " << pos << " != buffer " + << pos_remote << ") [" << dist << "] - tag: " << tag + << " comm from " << proc << " to " << rank); + } + } } +/* -------------------------------------------------------------------------- */ } // namespace akantu diff --git a/src/synchronizer/node_synchronizer.hh b/src/synchronizer/node_synchronizer.hh index 5378ef648..0679f0559 100644 --- a/src/synchronizer/node_synchronizer.hh +++ b/src/synchronizer/node_synchronizer.hh @@ -1,87 +1,98 @@ /** * @file node_synchronizer.hh * * @author Nicolas Richart * * @date creation: Tue Nov 08 2016 * @date last modification: Tue Feb 20 2018 * * @brief Synchronizer for nodal information * * @section LICENSE * * Copyright (©) 2016-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "mesh_events.hh" #include "synchronizer_impl.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_NODE_SYNCHRONIZER_HH__ #define __AKANTU_NODE_SYNCHRONIZER_HH__ namespace akantu { class NodeSynchronizer : public MeshEventHandler, public SynchronizerImpl { public: NodeSynchronizer(Mesh & mesh, const ID & id = "element_synchronizer", MemoryID memory_id = 0, const bool register_to_event_manager = true, EventHandlerPriority event_priority = _ehp_synchronizer); ~NodeSynchronizer() override; friend class NodeInfoPerProc; + UInt sanityCheckDataSize(const Array & nodes, + const SynchronizationTag & tag, + bool from_comm_desc) const override; + void packSanityCheckData(CommunicationBuffer & buffer, + const Array & nodes, + const SynchronizationTag & /*tag*/) const override; + void unpackSanityCheckData(CommunicationBuffer & buffer, + const Array & nodes, + const SynchronizationTag & tag, UInt proc, + UInt rank) const override; + /// function to implement to react on akantu::NewNodesEvent void onNodesAdded(const Array &, const NewNodesEvent &) override; /// function to implement to react on akantu::RemovedNodesEvent void onNodesRemoved(const Array &, const Array &, const RemovedNodesEvent &) override {} /// function to implement to react on akantu::NewElementsEvent void onElementsAdded(const Array &, const NewElementsEvent &) override {} /// function to implement to react on akantu::RemovedElementsEvent void onElementsRemoved(const Array &, const ElementTypeMapArray &, const RemovedElementsEvent &) override {} /// function to implement to react on akantu::ChangedElementsEvent void onElementsChanged(const Array &, const Array &, const ElementTypeMapArray &, const ChangedElementsEvent &) override {} public: AKANTU_GET_MACRO(Mesh, mesh, Mesh &); protected: Int getRank(const UInt & /*node*/) const final { AKANTU_TO_IMPLEMENT(); } protected: Mesh & mesh; // std::unordered_map node_to_prank; }; } // namespace akantu #endif /* __AKANTU_NODE_SYNCHRONIZER_HH__ */ diff --git a/src/synchronizer/synchronizer_impl_tmpl.hh b/src/synchronizer/synchronizer_impl_tmpl.hh index 2e9fde4d4..d31bfedb3 100644 --- a/src/synchronizer/synchronizer_impl_tmpl.hh +++ b/src/synchronizer/synchronizer_impl_tmpl.hh @@ -1,494 +1,495 @@ /** * @file synchronizer_impl_tmpl.hh * * @author Nicolas Richart * * @date creation: Wed Sep 07 2016 * @date last modification: Tue Feb 20 2018 * * @brief Implementation of the SynchronizerImpl * * @section LICENSE * * Copyright (©) 2016-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "synchronizer_impl.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_SYNCHRONIZER_IMPL_TMPL_HH__ #define __AKANTU_SYNCHRONIZER_IMPL_TMPL_HH__ namespace akantu { /* -------------------------------------------------------------------------- */ template SynchronizerImpl::SynchronizerImpl(const Communicator & comm, const ID & id, MemoryID memory_id) : Synchronizer(comm, id, memory_id), communications(comm) {} /* -------------------------------------------------------------------------- */ template SynchronizerImpl::SynchronizerImpl(const SynchronizerImpl & other, const ID & id) : Synchronizer(other), communications(other.communications) { this->id = id; } /* -------------------------------------------------------------------------- */ template void SynchronizerImpl::communicateOnce( const std::tuple & send_recv_schemes, const Tag::CommTags & comm_tag, DataAccessor & data_accessor, const SynchronizationTag & tag) const { // no need to synchronize if (this->nb_proc == 1) return; CommunicationSendRecv send_dir, recv_dir; std::tie(send_dir, recv_dir) = send_recv_schemes; using CommunicationRequests = std::vector; using CommunicationBuffers = std::map; CommunicationRequests send_requests, recv_requests; CommunicationBuffers send_buffers, recv_buffers; auto postComm = [&](const auto & sr, auto & buffers, auto & requests) -> void { for (auto && pair : communications.iterateSchemes(sr)) { auto & proc = pair.first; const auto & scheme = pair.second; if (scheme.size() == 0) continue; auto & buffer = buffers[proc]; auto buffer_size = data_accessor.getNbData(scheme, tag); if (buffer_size == 0) continue; #ifndef AKANTU_NDEBUG buffer_size += this->sanityCheckDataSize(scheme, tag, false); #endif buffer.resize(buffer_size); if (sr == recv_dir) { requests.push_back(communicator.asyncReceive( - buffer, proc, Tag::genTag(this->rank, 0, comm_tag, this->hash_id))); + buffer, proc, + Tag::genTag(this->rank, tag, comm_tag, this->hash_id))); } else { #ifndef AKANTU_NDEBUG this->packSanityCheckData(buffer, scheme, tag); #endif data_accessor.packData(buffer, scheme, tag); AKANTU_DEBUG_ASSERT( buffer.getPackedSize() == buffer.size(), "The data accessor did not pack all the data it " "promised in communication with tag " << tag << " (Promised: " << buffer.size() << "bytes, packed: " << buffer.getPackedSize() << "bytes [avg: " << Real(buffer.size() - buffer.getPackedSize()) / scheme.size() << "bytes per entity missing])"); send_requests.push_back(communicator.asyncSend( - buffer, proc, Tag::genTag(proc, 0, comm_tag, this->hash_id))); + buffer, proc, Tag::genTag(proc, tag, comm_tag, this->hash_id))); } } }; // post the receive requests postComm(recv_dir, recv_buffers, recv_requests); // post the send data requests postComm(send_dir, send_buffers, send_requests); // treat the receive requests UInt request_ready; while ((request_ready = communicator.waitAny(recv_requests)) != UInt(-1)) { auto & req = recv_requests[request_ready]; auto proc = req.getSource(); auto & buffer = recv_buffers[proc]; const auto & scheme = this->communications.getScheme(proc, recv_dir); #ifndef AKANTU_NDEBUG this->unpackSanityCheckData(buffer, scheme, tag, proc, this->rank); #endif data_accessor.unpackData(buffer, scheme, tag); AKANTU_DEBUG_ASSERT( buffer.getLeftToUnpack() == 0, "The data accessor ignored some data in communication with tag " << tag); req.free(); recv_requests.erase(recv_requests.begin() + request_ready); } communicator.waitAll(send_requests); communicator.freeCommunicationRequest(send_requests); } /* -------------------------------------------------------------------------- */ template void SynchronizerImpl::slaveReductionOnceImpl( DataAccessor & data_accessor, const SynchronizationTag & tag) const { communicateOnce(std::make_tuple(_recv, _send), Tag::_REDUCE, data_accessor, tag); } /* -------------------------------------------------------------------------- */ template void SynchronizerImpl::synchronizeOnceImpl( DataAccessor & data_accessor, const SynchronizationTag & tag) const { communicateOnce(std::make_tuple(_send, _recv), Tag::_SYNCHRONIZE, data_accessor, tag); } /* -------------------------------------------------------------------------- */ template void SynchronizerImpl::asynchronousSynchronizeImpl( const DataAccessor & data_accessor, const SynchronizationTag & tag) { AKANTU_DEBUG_IN(); if (not this->communications.hasCommunicationSize(tag)) this->computeBufferSize(data_accessor, tag); this->communications.incrementCounter(tag); // Posting the receive ------------------------------------------------------- if (this->communications.hasPendingRecv(tag)) { AKANTU_CUSTOM_EXCEPTION_INFO( debug::CommunicationException(), "There must still be some pending receive communications." << " Tag is " << tag << " Cannot start new ones"); } for (auto && comm_desc : this->communications.iterateRecv(tag)) { comm_desc.postRecv(this->hash_id); } // Posting the sends ------------------------------------------------------- if (communications.hasPendingSend(tag)) { AKANTU_CUSTOM_EXCEPTION_INFO( debug::CommunicationException(), "There must be some pending sending communications." << " Tag is " << tag); } for (auto && comm_desc : this->communications.iterateSend(tag)) { comm_desc.resetBuffer(); #ifndef AKANTU_NDEBUG this->packSanityCheckData(comm_desc); #endif comm_desc.packData(data_accessor); comm_desc.postSend(this->hash_id); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void SynchronizerImpl::waitEndSynchronizeImpl( DataAccessor & data_accessor, const SynchronizationTag & tag) { AKANTU_DEBUG_IN(); #ifndef AKANTU_NDEBUG if (this->communications.begin(tag, _recv) != this->communications.end(tag, _recv) && !this->communications.hasPendingRecv(tag)) AKANTU_CUSTOM_EXCEPTION_INFO(debug::CommunicationException(), "No pending communication with the tag \"" << tag); #endif auto recv_end = this->communications.end(tag, _recv); decltype(recv_end) recv_it; while ((recv_it = this->communications.waitAnyRecv(tag)) != recv_end) { auto && comm_desc = *recv_it; #ifndef AKANTU_NDEBUG this->unpackSanityCheckData(comm_desc); #endif comm_desc.unpackData(data_accessor); comm_desc.resetBuffer(); comm_desc.freeRequest(); } this->communications.waitAllSend(tag); this->communications.freeSendRequests(tag); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void SynchronizerImpl::computeAllBufferSizes( const DataAccessor & data_accessor) { for (auto && tag : this->communications.iterateTags()) { this->computeBufferSize(data_accessor, tag); } } /* -------------------------------------------------------------------------- */ template void SynchronizerImpl::computeBufferSizeImpl( const DataAccessor & data_accessor, const SynchronizationTag & tag) { AKANTU_DEBUG_IN(); if (not this->communications.hasCommunication(tag)) { this->communications.initializeCommunications(tag); AKANTU_DEBUG_ASSERT(communications.hasCommunication(tag) == true, "Communications where not properly initialized"); } for (auto sr : iterate_send_recv) { for (auto && pair : this->communications.iterateSchemes(sr)) { auto proc = pair.first; const auto & scheme = pair.second; UInt size = 0; #ifndef AKANTU_NDEBUG size += this->sanityCheckDataSize(scheme, tag); #endif size += data_accessor.getNbData(scheme, tag); AKANTU_DEBUG_INFO("I have " << size << "(" << printMemorySize(size) << " - " << scheme.size() << " element(s)) data to " << std::string(sr == _recv ? "receive from" : "send to") << proc << " for tag " << tag); this->communications.setCommunicationSize(tag, proc, size, sr); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template template void SynchronizerImpl::split(SynchronizerImpl & in_synchronizer, Pred && pred) { AKANTU_DEBUG_IN(); auto filter_list = [&](auto & list, auto & new_list) { auto copy = list; list.resize(0); new_list.resize(0); for (auto && entity : copy) { if (std::forward(pred)(entity)) { new_list.push_back(entity); } else { list.push_back(entity); } } }; for (auto sr : iterate_send_recv) { for (auto & scheme_pair : in_synchronizer.communications.iterateSchemes(sr)) { auto proc = scheme_pair.first; auto & scheme = scheme_pair.second; auto & new_scheme = communications.createScheme(proc, sr); filter_list(scheme, new_scheme); } } in_synchronizer.communications.invalidateSizes(); communications.invalidateSizes(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template template void SynchronizerImpl::updateSchemes(Updater && scheme_updater) { for (auto sr : iterate_send_recv) { for (auto & scheme_pair : communications.iterateSchemes(sr)) { auto proc = scheme_pair.first; auto & scheme = scheme_pair.second; std::forward(scheme_updater)(scheme, proc, sr); } } communications.invalidateSizes(); } /* -------------------------------------------------------------------------- */ template template void SynchronizerImpl::filterScheme(Pred && pred) { std::vector requests; std::unordered_map> keep_entities; auto filter_list = [](const auto & keep, auto & list) { Array new_list; for (const auto & keep_entity : keep) { const Entity & entity = list(keep_entity); new_list.push_back(entity); } list.copy(new_list); }; // loop over send_schemes for (auto & scheme_pair : communications.iterateSchemes(_recv)) { auto proc = scheme_pair.first; auto & scheme = scheme_pair.second; auto & keep_entity = keep_entities[proc]; for (auto && entity : enumerate(scheme)) { if (pred(std::get<1>(entity))) { keep_entity.push_back(std::get<0>(entity)); } } auto tag = Tag::genTag(this->rank, 0, Tag::_MODIFY_SCHEME); AKANTU_DEBUG_INFO("I have " << keep_entity.size() << " elements to still receive from processor " << proc << " (communication tag : " << tag << ")"); filter_list(keep_entity, scheme); requests.push_back(communicator.asyncSend(keep_entity, proc, tag)); } // clean the receive scheme for (auto & scheme_pair : communications.iterateSchemes(_send)) { auto proc = scheme_pair.first; auto & scheme = scheme_pair.second; auto tag = Tag::genTag(proc, 0, Tag::_MODIFY_SCHEME); AKANTU_DEBUG_INFO("Waiting list of elements to keep from processor " << proc << " (communication tag : " << tag << ")"); CommunicationStatus status; communicator.probe(proc, tag, status); Array keep_entity(status.size(), 1, "keep_element"); AKANTU_DEBUG_INFO("I have " << keep_entity.size() << " elements to keep in my send list to processor " << proc << " (communication tag : " << tag << ")"); communicator.receive(keep_entity, proc, tag); filter_list(keep_entity, scheme); } communicator.waitAll(requests); communicator.freeCommunicationRequest(requests); communications.invalidateSizes(); } /* -------------------------------------------------------------------------- */ template void SynchronizerImpl::swapSendRecv() { communications.swapSendRecv(); } /* -------------------------------------------------------------------------- */ template UInt SynchronizerImpl::sanityCheckDataSize(const Array &, const SynchronizationTag &, bool is_comm_desc) const { if (not is_comm_desc) { return 0; } UInt size = 0; size += sizeof(SynchronizationTag); // tag size += sizeof(UInt); // comm_desc.getNbData(); size += sizeof(UInt); // comm_desc.getProc(); size += sizeof(this->rank); // mesh.getCommunicator().whoAmI(); return size; } /* -------------------------------------------------------------------------- */ template void SynchronizerImpl::packSanityCheckData( CommunicationDescriptor & comm_desc) const { auto & buffer = comm_desc.getBuffer(); buffer << comm_desc.getTag(); buffer << comm_desc.getNbData(); buffer << comm_desc.getProc(); buffer << this->rank; const auto & tag = comm_desc.getTag(); const auto & send_element = comm_desc.getScheme(); this->packSanityCheckData(buffer, send_element, tag); } /* -------------------------------------------------------------------------- */ template void SynchronizerImpl::unpackSanityCheckData( CommunicationDescriptor & comm_desc) const { auto & buffer = comm_desc.getBuffer(); const auto & tag = comm_desc.getTag(); auto nb_data = comm_desc.getNbData(); auto proc = comm_desc.getProc(); auto rank = this->rank; decltype(nb_data) recv_nb_data; decltype(proc) recv_proc; decltype(rank) recv_rank; SynchronizationTag t; buffer >> t; buffer >> recv_nb_data; buffer >> recv_proc; buffer >> recv_rank; AKANTU_DEBUG_ASSERT( t == tag, "The tag received does not correspond to the tag expected"); AKANTU_DEBUG_ASSERT( nb_data == recv_nb_data, "The nb_data received does not correspond to the nb_data expected"); AKANTU_DEBUG_ASSERT(UInt(recv_rank) == proc, "The rank received does not correspond to the proc"); AKANTU_DEBUG_ASSERT(recv_proc == UInt(rank), "The proc received does not correspond to the rank"); auto & recv_element = comm_desc.getScheme(); this->unpackSanityCheckData(buffer, recv_element, tag, proc, rank); } /* -------------------------------------------------------------------------- */ } // namespace akantu #endif /* __AKANTU_SYNCHRONIZER_IMPL_TMPL_HH__ */