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