diff --git a/src/mesh_utils/cohesive_element_inserter_helper.cc b/src/mesh_utils/cohesive_element_inserter_helper.cc
index b50fea3d3..22f2777fa 100644
--- a/src/mesh_utils/cohesive_element_inserter_helper.cc
+++ b/src/mesh_utils/cohesive_element_inserter_helper.cc
@@ -1,947 +1,948 @@
/**
* @file cohesive_element_inserter_helper.cc
*
* @author Nicolas Richart
*
* @date creation jeu sep 03 2020
*
* @brief A Documented file.
*
* @section LICENSE
*
* Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "cohesive_element_inserter_helper.hh"
/* -------------------------------------------------------------------------- */
#include "data_accessor.hh"
#include "element_synchronizer.hh"
#include "fe_engine.hh"
#include "mesh_accessor.hh"
/* -------------------------------------------------------------------------- */
#include <queue>
/* -------------------------------------------------------------------------- */
namespace akantu {
/* -------------------------------------------------------------------------- */
CohesiveElementInserterHelper::CohesiveElementInserterHelper(
Mesh & mesh, const ElementTypeMapArray<bool> & facet_insertion)
: doubled_nodes(0, 2, "doubled_nodes"), mesh(mesh),
mesh_facets(mesh.getMeshFacets()) {
auto spatial_dimension = mesh_facets.getSpatialDimension();
for (auto gt : ghost_types) {
for (auto type : mesh_facets.elementTypes(_ghost_type = gt)) {
nb_new_facets(type, gt) = mesh_facets.getNbElement(type, gt);
}
}
std::array<Int, 2> nb_facet_to_insert{0, 0};
// creates a vector of the facets to insert
std::vector<Element> potential_facets_to_double;
for (auto gt_facet : ghost_types) {
for (auto type_facet :
mesh_facets.elementTypes(spatial_dimension - 1, gt_facet)) {
const auto & f_insertion = facet_insertion(type_facet, gt_facet);
auto & counter = nb_facet_to_insert[gt_facet == _not_ghost ? 0 : 1];
for (auto && data : enumerate(f_insertion)) {
if (std::get<1>(data)) {
UInt el = std::get<0>(data);
potential_facets_to_double.push_back({type_facet, el, gt_facet});
++counter;
}
}
}
}
// Defines a global order of insertion oof cohesive element to ensure node
// doubling appens in the smae order, this is necessary for the global node
// numbering
if (mesh_facets.isDistributed()) {
const auto & comm = mesh_facets.getCommunicator();
ElementTypeMapArray<Int> global_orderings;
global_orderings.initialize(mesh_facets,
_spatial_dimension = spatial_dimension - 1,
_with_nb_element = true, _default_value = -1);
auto starting_index = nb_facet_to_insert[0];
comm.exclusiveScan(starting_index);
// define the local numbers for facet to insert
for (auto gt_facet : ghost_types) {
for (auto type_facet :
mesh_facets.elementTypes(spatial_dimension - 1, gt_facet)) {
for (auto data : zip(facet_insertion(type_facet, gt_facet),
global_orderings(type_facet, gt_facet))) {
if (std::get<0>(data)) {
std::get<1>(data) = starting_index;
++starting_index;
}
}
}
}
// retreives the oorder number from neighoring processors
auto && synchronizer = mesh_facets.getElementSynchronizer();
SimpleElementDataAccessor<Int> data_accessor(
global_orderings, SynchronizationTag::_ce_insertion_order);
synchronizer.synchronizeOnce(data_accessor,
SynchronizationTag::_ce_insertion_order);
// sort the facets to double based on this global ordering
std::sort(potential_facets_to_double.begin(),
potential_facets_to_double.end(),
[&global_orderings](auto && el1, auto && el2) {
return global_orderings(el1) < global_orderings(el2);
});
}
for (auto d : arange(spatial_dimension)) {
facets_to_double_by_dim[d] = std::make_unique<Array<Element>>(
0, 2, "facets_to_double_" + std::to_string(d));
}
auto & facets_to_double = *facets_to_double_by_dim[spatial_dimension - 1];
MeshAccessor mesh_accessor(mesh_facets);
auto & elements_to_subelements = mesh_accessor.getElementToSubelement();
for (auto && facet_to_double : potential_facets_to_double) {
auto gt_facet = facet_to_double.ghost_type;
auto type_facet = facet_to_double.type;
auto & elements_to_facets = elements_to_subelements(type_facet, gt_facet);
auto & elements_to_facet = elements_to_facets(facet_to_double.element);
if (elements_to_facet[1].type == _not_defined
#if defined(AKANTU_COHESIVE_ELEMENT)
|| elements_to_facet[1].kind() == _ek_cohesive
#endif
) {
AKANTU_DEBUG_WARNING("attempt to double a facet on the boundary");
continue;
}
auto new_facet = nb_new_facets(type_facet, gt_facet)++;
facets_to_double.push_back(Vector<Element>{
facet_to_double, Element{type_facet, new_facet, gt_facet}});
/// update facet_to_element vector
auto & element_to_update = elements_to_facet[1];
auto el = element_to_update.element;
- auto & facets_to_elements = mesh_facets.getSubelementToElement(
+ const auto & facets_to_elements = mesh_facets.getSubelementToElement(
element_to_update.type, element_to_update.ghost_type);
auto facets_to_element = Vector<Element>(
make_view(facets_to_elements, facets_to_elements.getNbComponent())
.begin()[el]);
auto facet_to_update = std::find(facets_to_element.begin(),
facets_to_element.end(), facet_to_double);
AKANTU_DEBUG_ASSERT(facet_to_update != facets_to_element.end(),
"Facet not found");
facet_to_update->element = new_facet;
/// update elements connected to facet
const auto & first_facet_list = elements_to_facet;
elements_to_facets.push_back(first_facet_list);
/// set new and original facets as boundary facets
elements_to_facets(new_facet)[0] = elements_to_facets(new_facet)[1];
elements_to_facets(new_facet)[1] = ElementNull;
elements_to_facets(facet_to_double.element)[1] = ElementNull;
}
}
/* -------------------------------------------------------------------------- */
inline auto
CohesiveElementInserterHelper::hasElement(const Vector<UInt> & nodes_element,
const Vector<UInt> & nodes) -> bool {
// if one of the nodes of "nodes" is not in "nodes_element" it stops
auto it = std::mismatch(nodes.begin(), nodes.end(), nodes_element.begin(),
[&](auto && node, auto && /*node2*/) -> bool {
auto it = std::find(nodes_element.begin(),
nodes_element.end(), node);
return (it != nodes_element.end());
});
// true if all "nodes" where found in "nodes_element"
return (it.first == nodes.end());
}
/* -------------------------------------------------------------------------- */
inline auto CohesiveElementInserterHelper::removeElementsInVector(
const std::vector<Element> & elem_to_remove,
std::vector<Element> & elem_list) -> bool {
if (elem_list.size() <= elem_to_remove.size()) {
return true;
}
auto el_it = elem_to_remove.begin();
auto el_last = elem_to_remove.end();
std::vector<Element>::iterator el_del;
UInt deletions = 0;
for (; el_it != el_last; ++el_it) {
el_del = std::find(elem_list.begin(), elem_list.end(), *el_it);
if (el_del != elem_list.end()) {
elem_list.erase(el_del);
++deletions;
}
}
AKANTU_DEBUG_ASSERT(deletions == 0 || deletions == elem_to_remove.size(),
"Not all elements have been erased");
return deletions == 0;
}
/* -------------------------------------------------------------------------- */
void CohesiveElementInserterHelper::updateElementalConnectivity(
Mesh & mesh, UInt old_node, UInt new_node,
const std::vector<Element> & element_list,
const std::vector<Element> * facet_list) {
AKANTU_DEBUG_IN();
- for (auto & element : element_list) {
+ for (const auto & element : element_list) {
if (element.type == _not_defined) {
continue;
}
Vector<UInt> connectivity = mesh.getConnectivity(element);
if (element.kind() == _ek_cohesive) {
AKANTU_DEBUG_ASSERT(
facet_list != nullptr,
"Provide a facet list in order to update cohesive elements");
const Vector<Element> facets =
mesh_facets.getSubelementToElement(element);
auto facet_nb_nodes = connectivity.size() / 2;
/// loop over cohesive element's facets
for (const auto & facet : enumerate(facets)) {
/// skip facets if not present in the list
if (std::find(facet_list->begin(), facet_list->end(),
std::get<1>(facet)) == facet_list->end()) {
continue;
}
auto n = std::get<0>(facet);
- auto begin = connectivity.begin() + static_cast<UInt>(n * facet_nb_nodes);
+ auto begin =
+ connectivity.begin() + static_cast<UInt>(n * facet_nb_nodes);
auto end = begin + facet_nb_nodes;
auto it = std::find(begin, end, old_node);
AKANTU_DEBUG_ASSERT(it != end, "Node not found in current element");
*it = new_node;
}
} else {
auto it = std::find(connectivity.begin(), connectivity.end(), old_node);
AKANTU_DEBUG_ASSERT(it != connectivity.end(),
"Node not found in current element");
/// update connectivity
*it = new_node;
}
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void CohesiveElementInserterHelper::updateSubelementToElement(UInt dim,
bool facet_mode) {
auto & facets_to_double = *facets_to_double_by_dim[dim];
- auto & facets_to_subfacets =
- elementsOfDimToElementsOfDim(dim + static_cast<decltype(dim)>(facet_mode), dim);
+ auto & facets_to_subfacets = elementsOfDimToElementsOfDim(
+ dim + static_cast<decltype(dim)>(facet_mode), dim);
for (auto && data :
zip(make_view(facets_to_double, 2), facets_to_subfacets)) {
const auto & old_subfacet = std::get<0>(data)[0];
const auto & new_subfacet = std::get<0>(data)[1];
auto & facet_to_subfacets = std::get<1>(data);
MeshAccessor mesh_accessor(mesh_facets);
for (auto & facet : facet_to_subfacets) {
Vector<Element> subfacets = mesh_accessor.getSubelementToElement(facet);
auto && subfacet_to_update_it =
std::find(subfacets.begin(), subfacets.end(), old_subfacet);
AKANTU_DEBUG_ASSERT(subfacet_to_update_it != subfacets.end(),
"Subfacet not found");
*subfacet_to_update_it = new_subfacet;
}
}
}
/* -------------------------------------------------------------------------- */
void CohesiveElementInserterHelper::updateElementToSubelement(UInt dim,
bool facet_mode) {
auto & facets_to_double = *facets_to_double_by_dim[dim];
auto & facets_to_subfacets = elementsOfDimToElementsOfDim(
dim + static_cast<decltype(dim)>(facet_mode), dim);
MeshAccessor mesh_accessor(mesh_facets);
// resize the arrays
mesh_accessor.getElementToSubelement().initialize(
mesh_facets, _spatial_dimension = dim, _with_nb_element = true);
for (auto && data :
zip(make_view(facets_to_double, 2), facets_to_subfacets)) {
const auto & new_facet = std::get<0>(data)[1];
mesh_accessor.getElementToSubelement(new_facet) = std::get<1>(data);
}
}
/* -------------------------------------------------------------------------- */
void CohesiveElementInserterHelper::updateQuadraticSegments(UInt dim) {
AKANTU_DEBUG_IN();
auto spatial_dimension = mesh.getSpatialDimension();
auto & facets_to_double = *facets_to_double_by_dim[dim];
MeshAccessor mesh_accessor(mesh_facets);
auto & connectivities = mesh_accessor.getConnectivities();
/// this ones matter only for segments in 3D
Array<std::vector<Element>> * element_to_subfacet_double = nullptr;
Array<std::vector<Element>> * facet_to_subfacet_double = nullptr;
if (dim == spatial_dimension - 2) {
element_to_subfacet_double = &elementsOfDimToElementsOfDim(dim + 2, dim);
facet_to_subfacet_double = &elementsOfDimToElementsOfDim(dim + 1, dim);
}
- auto & element_to_subelement = mesh_facets.getElementToSubelement();
+ const auto & element_to_subelement = mesh_facets.getElementToSubelement();
std::vector<UInt> middle_nodes;
for (auto && facet_to_double : make_view(facets_to_double, 2)) {
auto & old_facet = facet_to_double[0];
if (old_facet.type != _segment_3) {
continue;
}
auto node = connectivities(old_facet, 2);
if (not mesh.isPureGhostNode(node)) {
middle_nodes.push_back(node);
}
}
auto n = doubled_nodes.size();
doubleNodes(middle_nodes);
for (auto && data : enumerate(make_view(facets_to_double, 2))) {
auto facet = std::get<0>(data);
auto & old_facet = std::get<1>(data)[0];
if (old_facet.type != _segment_3) {
continue;
}
auto old_node = connectivities(old_facet, 2);
if (mesh.isPureGhostNode(old_node)) {
continue;
}
auto new_node = doubled_nodes(n, 1);
auto & new_facet = std::get<1>(data)[1];
connectivities(new_facet, 2) = new_node;
if (dim == spatial_dimension - 2) {
updateElementalConnectivity(mesh_facets, old_node, new_node,
element_to_subelement(new_facet, 0));
updateElementalConnectivity(mesh, old_node, new_node,
(*element_to_subfacet_double)(facet),
&(*facet_to_subfacet_double)(facet));
} else {
updateElementalConnectivity(mesh, old_node, new_node,
element_to_subelement(new_facet, 0));
}
++n;
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
UInt CohesiveElementInserterHelper::insertCohesiveElement() {
auto nb_new_facets = insertFacetsOnly();
if (nb_new_facets == 0) {
return 0;
}
updateCohesiveData();
return nb_new_facets;
}
/* -------------------------------------------------------------------------- */
UInt CohesiveElementInserterHelper::insertFacetsOnly() {
UInt spatial_dimension = mesh.getSpatialDimension();
if (facets_to_double_by_dim[spatial_dimension - 1]->empty()) {
return 0;
}
if (spatial_dimension == 1) {
doublePointFacet();
} else if (spatial_dimension == 2) {
doubleFacets<1>();
findSubfacetToDouble<1>();
doubleSubfacet<2>();
} else if (spatial_dimension == 3) {
doubleFacets<2>();
findSubfacetToDouble<2>();
doubleFacets<1>();
findSubfacetToDouble<1>();
doubleSubfacet<3>();
}
return facets_to_double_by_dim[spatial_dimension - 1]->size();
}
/* -------------------------------------------------------------------------- */
template <UInt dim> void CohesiveElementInserterHelper::doubleFacets() {
AKANTU_DEBUG_IN();
NewElementsEvent new_facets;
auto spatial_dimension = mesh_facets.getSpatialDimension();
auto & facets_to_double = *facets_to_double_by_dim[dim];
MeshAccessor mesh_accessor(mesh_facets);
for (auto && facet_to_double : make_view(facets_to_double, 2)) {
auto && old_facet = facet_to_double[0];
auto && new_facet = facet_to_double[1];
auto & facets_connectivities =
mesh_accessor.getConnectivity(old_facet.type, old_facet.ghost_type);
facets_connectivities.resize(
nb_new_facets(old_facet.type, old_facet.ghost_type));
auto facets_connectivities_begin =
make_view(facets_connectivities, facets_connectivities.getNbComponent())
.begin();
// copy the connectivities
Vector<UInt> new_conn(facets_connectivities_begin[new_facet.element]);
Vector<UInt> old_conn(facets_connectivities_begin[old_facet.element]);
new_conn = old_conn;
// this will fail if multiple facet types exists for a given element type
// \TODO handle multiple sub-facet types
auto nb_subfacet_per_facet = Mesh::getNbFacetsPerElement(old_facet.type);
auto & subfacets_to_facets = mesh_accessor.getSubelementToElementNC(
old_facet.type, old_facet.ghost_type);
subfacets_to_facets.resize(
nb_new_facets(old_facet.type, old_facet.ghost_type), ElementNull);
auto subfacets_to_facets_begin =
make_view(subfacets_to_facets, nb_subfacet_per_facet).begin();
// copy the subfacet to facets information
Vector<Element> old_subfacets_to_facet(
subfacets_to_facets_begin[old_facet.element]);
Vector<Element> new_subfacet_to_facet(
subfacets_to_facets_begin[new_facet.element]);
new_subfacet_to_facet = old_subfacets_to_facet;
for (auto & subfacet : old_subfacets_to_facet) {
if (subfacet == ElementNull) {
continue;
}
/// update facet_to_subfacet array
mesh_accessor.getElementToSubelement(subfacet).push_back(new_facet);
}
new_facets.getList().push_back(new_facet);
}
/// update facet_to_subfacet and _segment_3 facets if any
if (dim != spatial_dimension - 1) {
updateSubelementToElement(dim, true);
updateElementToSubelement(dim, true);
updateQuadraticSegments(dim);
} else {
updateQuadraticSegments(dim);
}
mesh_facets.sendEvent(new_facets);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template <UInt dim> void CohesiveElementInserterHelper::findSubfacetToDouble() {
AKANTU_DEBUG_IN();
UInt spatial_dimension = mesh_facets.getSpatialDimension();
MeshAccessor mesh_accessor(mesh_facets);
auto & facets_to_double = *facets_to_double_by_dim[spatial_dimension - 1];
auto & subfacets_to_facets = mesh_facets.getSubelementToElement();
auto & elements_to_facets = mesh_accessor.getElementToSubelement();
/// loop on every facet
for (auto f : arange(2)) {
for (auto && facet_to_double : make_view(facets_to_double, 2)) {
auto old_facet = facet_to_double[f];
auto new_facet = facet_to_double[1 - f];
const auto & starting_element = elements_to_facets(new_facet, 0)[0];
auto current_facet = old_facet;
Vector<Element> subfacets_to_facet = subfacets_to_facets.get(old_facet);
/// loop on every subfacet
for (auto & subfacet : subfacets_to_facet) {
if (subfacet == ElementNull) {
continue;
}
if (dim == spatial_dimension - 2) {
Vector<Element> subsubfacets_to_subfacet =
subfacets_to_facets.get(subfacet);
/// loop on every subsubfacet
for (auto & subsubfacet : subsubfacets_to_subfacet) {
if (subsubfacet == ElementNull) {
continue;
}
Vector<UInt> subsubfacet_connectivity =
mesh_facets.getConnectivity(subsubfacet);
std::vector<Element> element_list;
std::vector<Element> facet_list;
std::vector<Element> subfacet_list;
/// check if subsubfacet is to be doubled
if (findElementsAroundSubfacet(
starting_element, current_facet, subsubfacet_connectivity,
element_list, facet_list, &subfacet_list) == false and
removeElementsInVector(
subfacet_list, elements_to_facets(subsubfacet)) == false) {
Element new_subsubfacet{
subsubfacet.type,
nb_new_facets(subsubfacet.type, subsubfacet.ghost_type)++,
subsubfacet.ghost_type};
facets_to_double_by_dim[dim - 1]->push_back(
Vector<Element>{subsubfacet, new_subsubfacet});
elementsOfDimToElementsOfDim(dim - 1, dim - 1)
.push_back(subfacet_list);
- elementsOfDimToElementsOfDim(dim, dim - 1)
- .push_back(facet_list);
+ elementsOfDimToElementsOfDim(dim, dim - 1).push_back(facet_list);
elementsOfDimToElementsOfDim(dim + 1, dim - 1)
.push_back(element_list);
}
}
} else {
std::vector<Element> element_list;
std::vector<Element> facet_list;
Vector<UInt> subfacet_connectivity =
mesh_facets.getConnectivity(subfacet);
/// check if subfacet is to be doubled
if (findElementsAroundSubfacet(starting_element, current_facet,
subfacet_connectivity, element_list,
facet_list) == false and
removeElementsInVector(facet_list,
elements_to_facets(subfacet)) == false) {
Element new_subfacet{
subfacet.type,
nb_new_facets(subfacet.type, subfacet.ghost_type)++,
subfacet.ghost_type};
facets_to_double_by_dim[dim - 1]->push_back(
Vector<Element>{subfacet, new_subfacet});
elementsOfDimToElementsOfDim(dim, dim - 1).push_back(facet_list);
elementsOfDimToElementsOfDim(dim + 1, dim - 1)
.push_back(element_list);
}
}
}
}
}
}
/* -------------------------------------------------------------------------- */
void CohesiveElementInserterHelper::doubleNodes(
const std::vector<UInt> & old_nodes) {
auto & position = mesh.getNodes();
auto spatial_dimension = mesh.getSpatialDimension();
auto old_nb_nodes = position.size();
position.reserve(old_nb_nodes + old_nodes.size());
doubled_nodes.reserve(doubled_nodes.size() + old_nodes.size());
auto position_begin = position.begin(spatial_dimension);
for (auto && data : enumerate(old_nodes)) {
auto n = std::get<0>(data);
auto old_node = std::get<1>(data);
decltype(old_node) new_node = old_nb_nodes + n;
/// store doubled nodes
doubled_nodes.push_back(Vector<UInt>{old_node, new_node});
/// update position
Vector<Real> coords = Vector<Real>(position_begin[old_node]);
position.push_back(coords);
}
}
/* -------------------------------------------------------------------------- */
bool CohesiveElementInserterHelper::findElementsAroundSubfacet(
const Element & starting_element, const Element & end_facet,
const Vector<UInt> & subfacet_connectivity,
std::vector<Element> & element_list, std::vector<Element> & facet_list,
std::vector<Element> * subfacet_list) {
bool facet_matched = false;
element_list.push_back(starting_element);
std::queue<Element> elements_to_check;
elements_to_check.push(starting_element);
/// keep going as long as there are elements to check
while (not elements_to_check.empty()) {
/// check current element
Element & current_element = elements_to_check.front();
const Vector<Element> facets_to_element =
mesh_facets.getSubelementToElement(current_element);
// for every facet of the element
for (auto & current_facet : facets_to_element) {
if (current_facet == ElementNull) {
continue;
}
if (current_facet == end_facet) {
facet_matched = true;
}
auto find_facet =
std::find(facet_list.begin(), facet_list.end(), current_facet);
// facet already listed or subfacet_connectivity is not in the
// connectivity of current_facet;
if ((find_facet != facet_list.end()) or
not hasElement(mesh_facets.getConnectivity(current_facet),
subfacet_connectivity)) {
continue;
}
facet_list.push_back(current_facet);
if (subfacet_list != nullptr) {
const Vector<Element> subfacets_of_facet =
mesh_facets.getSubelementToElement(current_facet);
/// check subfacets
for (const auto & current_subfacet : subfacets_of_facet) {
if (current_subfacet == ElementNull) {
continue;
}
auto find_subfacet = std::find(
subfacet_list->begin(), subfacet_list->end(), current_subfacet);
if ((find_subfacet != subfacet_list->end()) or
not hasElement(mesh_facets.getConnectivity(current_subfacet),
subfacet_connectivity)) {
continue;
}
subfacet_list->push_back(current_subfacet);
}
}
/// consider opposing element
const auto & elements_to_facet =
mesh_facets.getElementToSubelement(current_facet);
UInt opposing = 0;
if (elements_to_facet[0] == current_element) {
opposing = 1;
}
auto & opposing_element = elements_to_facet[opposing];
/// skip null elements since they are on a boundary
if (opposing_element == ElementNull) {
continue;
}
/// skip this element if already added
if (std::find(element_list.begin(), element_list.end(),
opposing_element) != element_list.end()) {
continue;
}
/// only regular elements have to be checked
if (opposing_element.kind() == _ek_regular) {
elements_to_check.push(opposing_element);
}
element_list.push_back(opposing_element);
AKANTU_DEBUG_ASSERT(hasElement(mesh.getConnectivity(opposing_element),
subfacet_connectivity),
"Subfacet doesn't belong to this element");
}
/// erased checked element from the list
elements_to_check.pop();
}
return facet_matched;
}
/* -------------------------------------------------------------------------- */
void CohesiveElementInserterHelper::updateCohesiveData() {
UInt spatial_dimension = mesh.getSpatialDimension();
bool third_dimension = spatial_dimension == 3;
MeshAccessor mesh_accessor(mesh);
MeshAccessor mesh_facets_accessor(mesh_facets);
for (auto ghost_type : ghost_types) {
for (auto cohesive_type : mesh.elementTypes(_element_kind = _ek_cohesive)) {
auto nb_cohesive_elements = mesh.getNbElement(cohesive_type, ghost_type);
nb_new_facets(cohesive_type, ghost_type) = nb_cohesive_elements;
mesh_facets_accessor.getSubelementToElement(cohesive_type, ghost_type);
}
}
auto & facets_to_double = *facets_to_double_by_dim[spatial_dimension - 1];
new_elements.reserve(new_elements.size() + facets_to_double.size());
std::array<Element, 2> facets;
auto & element_to_facet = mesh_facets_accessor.getElementToSubelement();
auto & facets_to_cohesive_elements =
mesh_facets_accessor.getSubelementToElement();
for (auto && facet_to_double : make_view(facets_to_double, 2)) {
auto & old_facet = facet_to_double[0];
/// (in 3D cohesive elements connectivity is inverted)
facets[third_dimension ? 1 : 0] = old_facet;
facets[third_dimension ? 0 : 1] = facet_to_double[1];
auto type_cohesive = FEEngine::getCohesiveElementType(old_facet.type);
auto & facet_connectivities =
mesh_facets.getConnectivity(old_facet.type, old_facet.ghost_type);
auto facet_connectivity_it =
facet_connectivities.begin(facet_connectivities.getNbComponent());
auto cohesive_element = Element{
type_cohesive, nb_new_facets(type_cohesive, old_facet.ghost_type)++,
old_facet.ghost_type};
auto & cohesives_connectivities =
mesh_accessor.getConnectivity(type_cohesive, old_facet.ghost_type);
Matrix<UInt> connectivity(facet_connectivities.getNbComponent(), 2);
Vector<Element> facets_to_cohesive_element(2);
for (auto s : arange(2)) {
/// store doubled facets
facets_to_cohesive_element[s] = facets[s];
// update connectivities
connectivity(s) = Vector<UInt>(facet_connectivity_it[facets[s].element]);
/// update element_to_facet vectors
element_to_facet(facets[s], 0)[1] = cohesive_element;
}
cohesives_connectivities.push_back(connectivity);
facets_to_cohesive_elements(type_cohesive, old_facet.ghost_type)
.push_back(facets_to_cohesive_element);
/// add cohesive element to the element event list
new_elements.push_back(cohesive_element);
}
}
/* -------------------------------------------------------------------------- */
void CohesiveElementInserterHelper::doublePointFacet() {
UInt spatial_dimension = mesh.getSpatialDimension();
if (spatial_dimension != 1) {
return;
}
NewElementsEvent new_facets_event;
auto & facets_to_double = *facets_to_double_by_dim[spatial_dimension - 1];
- auto & element_to_facet = mesh_facets.getElementToSubelement();
+ const auto & element_to_facet = mesh_facets.getElementToSubelement();
auto & position = mesh.getNodes();
MeshAccessor mesh_accessor(mesh_facets);
for (auto ghost_type : ghost_types) {
for (auto facet_type : nb_new_facets.elementTypes(
spatial_dimension - 1, ghost_type, _ek_regular)) {
auto nb_new_element = nb_new_facets(facet_type, ghost_type);
auto & connectivities =
mesh_accessor.getConnectivity(facet_type, ghost_type);
- connectivities.resize(connectivities.size() + nb_new_element);
+ connectivities.resize(nb_new_element);
}
}
+ position.reserve(position.size() + facets_to_double.size());
for (auto facet_to_double : make_view(facets_to_double, 2)) {
auto & old_facet = facet_to_double[0];
auto & new_facet = facet_to_double[1];
auto element = element_to_facet(new_facet)[0];
auto & facet_connectivities =
mesh_accessor.getConnectivity(old_facet.type, old_facet.ghost_type);
auto old_node = facet_connectivities(old_facet.element);
auto new_node = position.size();
/// update position
position.push_back(position(old_node));
facet_connectivities(new_facet.element) = new_node;
Vector<UInt> segment_conectivity = mesh.getConnectivity(element);
/// update facet connectivity
auto it = std::find(segment_conectivity.begin(), segment_conectivity.end(),
old_node);
*it = new_node;
doubled_nodes.push_back(Vector<UInt>{old_node, new_node});
new_facets_event.getList().push_back(new_facet);
}
mesh_facets.sendEvent(new_facets_event);
}
/* -------------------------------------------------------------------------- */
template <UInt spatial_dimension>
void CohesiveElementInserterHelper::doubleSubfacet() {
if (spatial_dimension == 1) {
return;
}
NewElementsEvent new_facets_event;
std::vector<UInt> nodes_to_double;
MeshAccessor mesh_accessor(mesh_facets);
auto & connectivities = mesh_accessor.getConnectivities();
auto & facets_to_double = *facets_to_double_by_dim[0];
ElementTypeMap<Int> nb_new_elements;
for (auto && facet_to_double : make_view(facets_to_double, 2)) {
auto && old_element = facet_to_double[0];
nb_new_elements(old_element.type, old_element.ghost_type) = 0;
}
for (auto && facet_to_double : make_view(facets_to_double, 2)) {
auto && old_element = facet_to_double[0];
++nb_new_elements(old_element.type, old_element.ghost_type);
}
for (auto ghost_type : ghost_types) {
for (auto facet_type :
nb_new_elements.elementTypes(0, ghost_type, _ek_regular)) {
auto & connectivities =
mesh_accessor.getConnectivity(facet_type, ghost_type);
connectivities.resize(connectivities.size() +
nb_new_elements(facet_type, ghost_type));
}
}
for (auto && facet_to_double : make_view(facets_to_double, 2)) {
auto & old_facet = facet_to_double(0);
// auto & new_facet = facet_to_double(1);
AKANTU_DEBUG_ASSERT(
old_facet.type == _point_1,
"Only _point_1 subfacet doubling is supported at the moment");
/// list nodes double
nodes_to_double.push_back(connectivities(old_facet));
}
auto old_nb_doubled_nodes = doubled_nodes.size();
doubleNodes(nodes_to_double);
auto double_nodes_view = make_view(doubled_nodes, 2);
for (auto && data :
zip(make_view(facets_to_double, 2),
range(double_nodes_view.begin() + old_nb_doubled_nodes,
double_nodes_view.end()),
arange(facets_to_double.size()))) {
// auto & old_facet = std::get<0>(data)[0];
auto & new_facet = std::get<0>(data)[1];
new_facets_event.getList().push_back(new_facet);
auto & nodes = std::get<1>(data);
auto old_node = nodes(0);
auto new_node = nodes(1);
auto f = std::get<2>(data);
/// add new nodes in connectivity
connectivities(new_facet) = new_node;
updateElementalConnectivity(mesh, old_node, new_node,
elementsOfDimToElementsOfDim(2, 0)(f),
&elementsOfDimToElementsOfDim(1, 0)(f));
updateElementalConnectivity(mesh_facets, old_node, new_node,
elementsOfDimToElementsOfDim(1, 0)(f));
if (spatial_dimension == 3) {
updateElementalConnectivity(mesh_facets, old_node, new_node,
elementsOfDimToElementsOfDim(0, 0)(f));
}
}
updateSubelementToElement(0, spatial_dimension == 2);
updateElementToSubelement(0, spatial_dimension == 2);
mesh_facets.sendEvent(new_facets_event);
}
} // namespace akantu
diff --git a/src/model/model.cc b/src/model/model.cc
index bc948054c..02c111a8e 100644
--- a/src/model/model.cc
+++ b/src/model/model.cc
@@ -1,347 +1,349 @@
/**
* @file model.cc
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
* @author David Simon Kammer <david.kammer@epfl.ch>
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @date creation: Mon Oct 03 2011
* @date last modification: Tue Feb 20 2018
*
* @brief implementation of model common parts
*
*
* 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 <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "model.hh"
#include "communicator.hh"
#include "data_accessor.hh"
#include "element_group.hh"
#include "element_synchronizer.hh"
#include "synchronizer_registry.hh"
/* -------------------------------------------------------------------------- */
namespace akantu {
/* -------------------------------------------------------------------------- */
Model::Model(Mesh & mesh, const ModelType & type,
std::shared_ptr<DOFManager> dof_manager, UInt dim, const ID & id,
const MemoryID & memory_id)
: Memory(id, memory_id),
- ModelSolver(mesh, type, id, memory_id, dof_manager), mesh(mesh),
+ ModelSolver(mesh, type, id, memory_id, std::move(dof_manager)),
+ mesh(mesh),
spatial_dimension(dim == _all_dimensions ? mesh.getSpatialDimension()
: dim),
parser(getStaticParser()) {
- AKANTU_DEBUG_IN();
-
this->mesh.registerEventHandler(*this, _ehp_model);
-
- AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
Model::Model(Mesh & mesh, const ModelType & type, UInt dim, const ID & id,
const MemoryID & memory_id)
- : Model(mesh, type, nullptr, dim, id, memory_id) {}
+ : Memory(id, memory_id), ModelSolver(mesh, type, id, memory_id), mesh(mesh),
+ spatial_dimension(dim == _all_dimensions ? mesh.getSpatialDimension()
+ : dim),
+ parser(getStaticParser()) {
+ this->mesh.registerEventHandler(*this, _ehp_model);
+}
/* -------------------------------------------------------------------------- */
Model::~Model() = default;
/* -------------------------------------------------------------------------- */
void Model::initFullImpl(const ModelOptions & options) {
AKANTU_DEBUG_IN();
method = options.analysis_method;
if (!this->hasDefaultSolver()) {
this->initNewSolver(this->method);
}
initModel();
initFEEngineBoundary();
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void Model::initNewSolver(const AnalysisMethod & method) {
ID solver_name;
TimeStepSolverType tss_type;
std::tie(solver_name, tss_type) = this->getDefaultSolverID(method);
if (not this->hasSolver(solver_name)) {
ModelSolverOptions options = this->getDefaultSolverOptions(tss_type);
this->getNewSolver(solver_name, tss_type, options.non_linear_solver_type);
for (auto && is_type : options.integration_scheme_type) {
if (!this->hasIntegrationScheme(solver_name, is_type.first)) {
this->setIntegrationScheme(solver_name, is_type.first, is_type.second,
options.solution_type[is_type.first]);
}
}
}
this->method = method;
this->setDefaultSolver(solver_name);
}
/* -------------------------------------------------------------------------- */
void Model::initFEEngineBoundary() {
FEEngine & fem_boundary = getFEEngineBoundary();
fem_boundary.initShapeFunctions(_not_ghost);
fem_boundary.initShapeFunctions(_ghost);
fem_boundary.computeNormalsOnIntegrationPoints(_not_ghost);
fem_boundary.computeNormalsOnIntegrationPoints(_ghost);
}
/* -------------------------------------------------------------------------- */
void Model::dumpGroup(const std::string & group_name) {
ElementGroup & group = mesh.getElementGroup(group_name);
group.dump();
}
/* -------------------------------------------------------------------------- */
void Model::dumpGroup(const std::string & group_name,
const std::string & dumper_name) {
ElementGroup & group = mesh.getElementGroup(group_name);
group.dump(dumper_name);
}
/* -------------------------------------------------------------------------- */
void Model::dumpGroup() {
for (auto & group : mesh.iterateElementGroups()) {
group.dump();
}
}
/* -------------------------------------------------------------------------- */
void Model::setGroupDirectory(const std::string & directory) {
for (auto & group : mesh.iterateElementGroups()) {
group.setDirectory(directory);
}
}
/* -------------------------------------------------------------------------- */
void Model::setGroupDirectory(const std::string & directory,
const std::string & group_name) {
ElementGroup & group = mesh.getElementGroup(group_name);
group.setDirectory(directory);
}
/* -------------------------------------------------------------------------- */
void Model::setGroupBaseName(const std::string & basename,
const std::string & group_name) {
ElementGroup & group = mesh.getElementGroup(group_name);
group.setBaseName(basename);
}
/* -------------------------------------------------------------------------- */
DumperIOHelper & Model::getGroupDumper(const std::string & group_name) {
ElementGroup & group = mesh.getElementGroup(group_name);
return group.getDumper();
}
/* -------------------------------------------------------------------------- */
// DUMPER stuff
/* -------------------------------------------------------------------------- */
void Model::addDumpGroupFieldToDumper(const std::string & field_id,
std::shared_ptr<dumpers::Field> field,
DumperIOHelper & dumper) {
#ifdef AKANTU_USE_IOHELPER
dumper.registerField(field_id, std::move(field));
#endif
}
/* -------------------------------------------------------------------------- */
void Model::addDumpField(const std::string & field_id) {
this->addDumpFieldToDumper(mesh.getDefaultDumperName(), field_id);
}
/* -------------------------------------------------------------------------- */
void Model::addDumpFieldVector(const std::string & field_id) {
this->addDumpFieldVectorToDumper(mesh.getDefaultDumperName(), field_id);
}
/* -------------------------------------------------------------------------- */
void Model::addDumpFieldTensor(const std::string & field_id) {
this->addDumpFieldTensorToDumper(mesh.getDefaultDumperName(), field_id);
}
/* -------------------------------------------------------------------------- */
void Model::setBaseName(const std::string & field_id) {
mesh.setBaseName(field_id);
}
/* -------------------------------------------------------------------------- */
void Model::setBaseNameToDumper(const std::string & dumper_name,
const std::string & basename) {
mesh.setBaseNameToDumper(dumper_name, basename);
}
/* -------------------------------------------------------------------------- */
void Model::addDumpFieldToDumper(const std::string & dumper_name,
const std::string & field_id) {
this->addDumpGroupFieldToDumper(dumper_name, field_id, "all", _ek_regular,
false);
}
/* -------------------------------------------------------------------------- */
void Model::addDumpGroupField(const std::string & field_id,
const std::string & group_name) {
ElementGroup & group = mesh.getElementGroup(group_name);
this->addDumpGroupFieldToDumper(group.getDefaultDumperName(), field_id,
group_name, _ek_regular, false);
}
/* -------------------------------------------------------------------------- */
void Model::removeDumpGroupField(const std::string & field_id,
const std::string & group_name) {
ElementGroup & group = mesh.getElementGroup(group_name);
this->removeDumpGroupFieldFromDumper(group.getDefaultDumperName(), field_id,
group_name);
}
/* -------------------------------------------------------------------------- */
void Model::removeDumpGroupFieldFromDumper(const std::string & dumper_name,
const std::string & field_id,
const std::string & group_name) {
ElementGroup & group = mesh.getElementGroup(group_name);
group.removeDumpFieldFromDumper(dumper_name, field_id);
}
/* -------------------------------------------------------------------------- */
void Model::addDumpFieldVectorToDumper(const std::string & dumper_name,
const std::string & field_id) {
this->addDumpGroupFieldToDumper(dumper_name, field_id, "all", _ek_regular,
true);
}
/* -------------------------------------------------------------------------- */
void Model::addDumpGroupFieldVector(const std::string & field_id,
const std::string & group_name) {
ElementGroup & group = mesh.getElementGroup(group_name);
this->addDumpGroupFieldVectorToDumper(group.getDefaultDumperName(), field_id,
group_name);
}
/* -------------------------------------------------------------------------- */
void Model::addDumpGroupFieldVectorToDumper(const std::string & dumper_name,
const std::string & field_id,
const std::string & group_name) {
this->addDumpGroupFieldToDumper(dumper_name, field_id, group_name,
_ek_regular, true);
}
/* -------------------------------------------------------------------------- */
void Model::addDumpFieldTensorToDumper(const std::string & dumper_name,
const std::string & field_id) {
this->addDumpGroupFieldToDumper(dumper_name, field_id, "all", _ek_regular,
true);
}
/* -------------------------------------------------------------------------- */
void Model::addDumpGroupFieldToDumper(const std::string & dumper_name,
const std::string & field_id,
const std::string & group_name,
ElementKind element_kind,
bool padding_flag) {
this->addDumpGroupFieldToDumper(dumper_name, field_id, group_name,
this->spatial_dimension, element_kind,
padding_flag);
}
/* -------------------------------------------------------------------------- */
void Model::addDumpGroupFieldToDumper(const std::string & dumper_name,
const std::string & field_id,
const std::string & group_name,
UInt spatial_dimension,
ElementKind element_kind,
bool padding_flag) {
#ifdef AKANTU_USE_IOHELPER
std::shared_ptr<dumpers::Field> field;
if (!field) {
field = this->createNodalFieldReal(field_id, group_name, padding_flag);
}
if (!field) {
field = this->createNodalFieldUInt(field_id, group_name, padding_flag);
}
if (!field) {
field = this->createNodalFieldBool(field_id, group_name, padding_flag);
}
if (!field) {
field = this->createElementalField(field_id, group_name, padding_flag,
spatial_dimension, element_kind);
}
if (!field) {
field = this->mesh.createFieldFromAttachedData<UInt>(field_id, group_name,
element_kind);
}
if (!field) {
field = this->mesh.createFieldFromAttachedData<Real>(field_id, group_name,
element_kind);
}
#ifndef AKANTU_NDEBUG
if (!field) {
AKANTU_DEBUG_WARNING("No field could be found based on name: " << field_id);
}
#endif
if (field) {
DumperIOHelper & dumper = mesh.getGroupDumper(dumper_name, group_name);
this->addDumpGroupFieldToDumper(field_id, field, dumper);
}
#endif
}
/* -------------------------------------------------------------------------- */
void Model::dump() { mesh.dump(); }
/* -------------------------------------------------------------------------- */
void Model::setDirectory(const std::string & directory) {
mesh.setDirectory(directory);
}
/* -------------------------------------------------------------------------- */
void Model::setDirectoryToDumper(const std::string & dumper_name,
const std::string & directory) {
mesh.setDirectoryToDumper(dumper_name, directory);
}
/* -------------------------------------------------------------------------- */
void Model::setTextModeToDumper() { mesh.setTextModeToDumper(); }
/* -------------------------------------------------------------------------- */
} // namespace akantu
diff --git a/src/model/solid_mechanics/solid_mechanics_model.cc b/src/model/solid_mechanics/solid_mechanics_model.cc
index 73c9c27c6..f35e1755b 100644
--- a/src/model/solid_mechanics/solid_mechanics_model.cc
+++ b/src/model/solid_mechanics/solid_mechanics_model.cc
@@ -1,1242 +1,1242 @@
/**
* @file solid_mechanics_model.cc
*
* @author Ramin Aghababaei <ramin.aghababaei@epfl.ch>
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
* @author Aurelia Isabel Cuba Ramos <aurelia.cubaramos@epfl.ch>
* @author David Simon Kammer <david.kammer@epfl.ch>
* @author Daniel Pino Muñoz <daniel.pinomunoz@epfl.ch>
* @author Nicolas Richart <nicolas.richart@epfl.ch>
* @author Clement Roux <clement.roux@epfl.ch>
* @author Marco Vocialta <marco.vocialta@epfl.ch>
*
* @date creation: Tue Jul 27 2010
* @date last modification: Wed Feb 21 2018
*
* @brief Implementation of the SolidMechanicsModel class
*
*
* 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 <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "solid_mechanics_model.hh"
#include "integrator_gauss.hh"
#include "shape_lagrange.hh"
#include "solid_mechanics_model_tmpl.hh"
#include "communicator.hh"
#include "element_synchronizer.hh"
#include "sparse_matrix.hh"
#include "synchronizer_registry.hh"
#include "dumpable_inline_impl.hh"
#ifdef AKANTU_USE_IOHELPER
#include "dumper_iohelper_paraview.hh"
#endif
#include "material_non_local.hh"
/* -------------------------------------------------------------------------- */
namespace akantu {
/* -------------------------------------------------------------------------- */
/**
* A solid mechanics model need a mesh and a dimension to be created. the model
* by it self can not do a lot, the good init functions should be called in
* order to configure the model depending on what we want to do.
*
* @param mesh mesh representing the model we want to simulate
* @param dim spatial dimension of the problem, if dim = 0 (default value) the
* dimension of the problem is assumed to be the on of the mesh
* @param id an id to identify the model
* @param memory_id the id of the memory
* @param model_type this is an internal parameter for inheritance purposes
*/
SolidMechanicsModel::SolidMechanicsModel(
Mesh & mesh, UInt dim, const ID & id, const MemoryID & memory_id,
std::shared_ptr<DOFManager> dof_manager, const ModelType model_type)
- : Model(mesh, model_type, dof_manager, dim, id, memory_id),
+ : Model(mesh, model_type, std::move(dof_manager), dim, id, memory_id),
BoundaryCondition<SolidMechanicsModel>(),
material_index("material index", id, memory_id),
material_local_numbering("material local numbering", id, memory_id) {
AKANTU_DEBUG_IN();
this->registerFEEngineObject<MyFEEngineType>("SolidMechanicsFEEngine", mesh,
Model::spatial_dimension);
#if defined(AKANTU_USE_IOHELPER)
this->mesh.registerDumper<DumperParaview>("solid_mechanics_model", id, true);
this->mesh.addDumpMesh(mesh, Model::spatial_dimension, _not_ghost,
_ek_regular);
#endif
material_selector = std::make_shared<DefaultMaterialSelector>(material_index);
this->registerDataAccessor(*this);
if (this->mesh.isDistributed()) {
auto & synchronizer = this->mesh.getElementSynchronizer();
this->registerSynchronizer(synchronizer, SynchronizationTag::_material_id);
this->registerSynchronizer(synchronizer, SynchronizationTag::_smm_mass);
this->registerSynchronizer(synchronizer, SynchronizationTag::_smm_stress);
this->registerSynchronizer(synchronizer, SynchronizationTag::_for_dump);
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
SolidMechanicsModel::~SolidMechanicsModel() = default;
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::setTimeStep(Real time_step, const ID & solver_id) {
Model::setTimeStep(time_step, solver_id);
#if defined(AKANTU_USE_IOHELPER)
this->mesh.getDumper().setTimeStep(time_step);
#endif
}
/* -------------------------------------------------------------------------- */
/* Initialization */
/* -------------------------------------------------------------------------- */
/**
* This function groups many of the initialization in on function. For most of
* basics case the function should be enough. The functions initialize the
* model, the internal vectors, set them to 0, and depending on the parameters
* it also initialize the explicit or implicit solver.
*
* @param options
* \parblock
* contains the different options to initialize the model
* \li \c analysis_method specify the type of solver to use
* \endparblock
*/
void SolidMechanicsModel::initFullImpl(const ModelOptions & options) {
material_index.initialize(mesh, _element_kind = _ek_not_defined,
_default_value = UInt(-1), _with_nb_element = true);
material_local_numbering.initialize(mesh, _element_kind = _ek_not_defined,
_with_nb_element = true);
Model::initFullImpl(options);
// initialize the materials
if (not this->parser.getLastParsedFile().empty()) {
this->instantiateMaterials();
this->initMaterials();
}
this->initBC(*this, *displacement, *displacement_increment, *external_force);
}
/* -------------------------------------------------------------------------- */
TimeStepSolverType SolidMechanicsModel::getDefaultSolverType() const {
return TimeStepSolverType::_dynamic_lumped;
}
/* -------------------------------------------------------------------------- */
ModelSolverOptions SolidMechanicsModel::getDefaultSolverOptions(
const TimeStepSolverType & type) const {
ModelSolverOptions options;
switch (type) {
case TimeStepSolverType::_dynamic_lumped: {
options.non_linear_solver_type = NonLinearSolverType::_lumped;
options.integration_scheme_type["displacement"] =
IntegrationSchemeType::_central_difference;
options.solution_type["displacement"] = IntegrationScheme::_acceleration;
break;
}
case TimeStepSolverType::_static: {
options.non_linear_solver_type = NonLinearSolverType::_newton_raphson;
options.integration_scheme_type["displacement"] =
IntegrationSchemeType::_pseudo_time;
options.solution_type["displacement"] = IntegrationScheme::_not_defined;
break;
}
case TimeStepSolverType::_dynamic: {
if (this->method == _explicit_consistent_mass) {
options.non_linear_solver_type = NonLinearSolverType::_newton_raphson;
options.integration_scheme_type["displacement"] =
IntegrationSchemeType::_central_difference;
options.solution_type["displacement"] = IntegrationScheme::_acceleration;
} else {
options.non_linear_solver_type = NonLinearSolverType::_newton_raphson;
options.integration_scheme_type["displacement"] =
IntegrationSchemeType::_trapezoidal_rule_2;
options.solution_type["displacement"] = IntegrationScheme::_displacement;
}
break;
}
default:
AKANTU_EXCEPTION(type << " is not a valid time step solver type");
}
return options;
}
/* -------------------------------------------------------------------------- */
std::tuple<ID, TimeStepSolverType>
SolidMechanicsModel::getDefaultSolverID(const AnalysisMethod & method) {
switch (method) {
case _explicit_lumped_mass: {
return std::make_tuple("explicit_lumped",
TimeStepSolverType::_dynamic_lumped);
}
case _explicit_consistent_mass: {
return std::make_tuple("explicit", TimeStepSolverType::_dynamic);
}
case _static: {
return std::make_tuple("static", TimeStepSolverType::_static);
}
case _implicit_dynamic: {
return std::make_tuple("implicit", TimeStepSolverType::_dynamic);
}
default:
return std::make_tuple("unknown", TimeStepSolverType::_not_defined);
}
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::initSolver(TimeStepSolverType time_step_solver_type,
NonLinearSolverType /*unused*/) {
auto & dof_manager = this->getDOFManager();
/* ------------------------------------------------------------------------ */
// for alloc type of solvers
this->allocNodalField(this->displacement, spatial_dimension, "displacement");
this->allocNodalField(this->previous_displacement, spatial_dimension,
"previous_displacement");
this->allocNodalField(this->displacement_increment, spatial_dimension,
"displacement_increment");
this->allocNodalField(this->internal_force, spatial_dimension,
"internal_force");
this->allocNodalField(this->external_force, spatial_dimension,
"external_force");
this->allocNodalField(this->blocked_dofs, spatial_dimension, "blocked_dofs");
this->allocNodalField(this->current_position, spatial_dimension,
"current_position");
// initialize the current positions
this->current_position->copy(this->mesh.getNodes());
/* ------------------------------------------------------------------------ */
if (!dof_manager.hasDOFs("displacement")) {
dof_manager.registerDOFs("displacement", *this->displacement, _dst_nodal);
dof_manager.registerBlockedDOFs("displacement", *this->blocked_dofs);
dof_manager.registerDOFsIncrement("displacement",
*this->displacement_increment);
dof_manager.registerDOFsPrevious("displacement",
*this->previous_displacement);
}
/* ------------------------------------------------------------------------ */
// for dynamic
if (time_step_solver_type == TimeStepSolverType::_dynamic ||
time_step_solver_type == TimeStepSolverType::_dynamic_lumped) {
this->allocNodalField(this->velocity, spatial_dimension, "velocity");
this->allocNodalField(this->acceleration, spatial_dimension,
"acceleration");
if (!dof_manager.hasDOFsDerivatives("displacement", 1)) {
dof_manager.registerDOFsDerivative("displacement", 1, *this->velocity);
dof_manager.registerDOFsDerivative("displacement", 2,
*this->acceleration);
}
}
}
/* -------------------------------------------------------------------------- */
/**
* Initialize the model,basically it pre-compute the shapes, shapes derivatives
* and jacobian
*/
void SolidMechanicsModel::initModel() {
/// \todo add the current position as a parameter to initShapeFunctions for
/// large deformation
getFEEngine().initShapeFunctions(_not_ghost);
getFEEngine().initShapeFunctions(_ghost);
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::assembleResidual() {
AKANTU_DEBUG_IN();
/* ------------------------------------------------------------------------ */
// computes the internal forces
this->assembleInternalForces();
/* ------------------------------------------------------------------------ */
this->getDOFManager().assembleToResidual("displacement",
*this->external_force, 1);
this->getDOFManager().assembleToResidual("displacement",
*this->internal_force, 1);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::assembleResidual(const ID & residual_part) {
AKANTU_DEBUG_IN();
if ("external" == residual_part) {
this->getDOFManager().assembleToResidual("displacement",
*this->external_force, 1);
AKANTU_DEBUG_OUT();
return;
}
if ("internal" == residual_part) {
this->assembleInternalForces();
this->getDOFManager().assembleToResidual("displacement",
*this->internal_force, 1);
AKANTU_DEBUG_OUT();
return;
}
AKANTU_CUSTOM_EXCEPTION(
debug::SolverCallbackResidualPartUnknown(residual_part));
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
MatrixType SolidMechanicsModel::getMatrixType(const ID & matrix_id) {
// \TODO check the materials to know what is the correct answer
if (matrix_id == "C") {
return _mt_not_defined;
}
if (matrix_id == "K") {
auto matrix_type = _unsymmetric;
for (auto & material : materials) {
matrix_type = std::max(matrix_type, material->getMatrixType(matrix_id));
}
}
return _symmetric;
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::assembleMatrix(const ID & matrix_id) {
if (matrix_id == "K") {
this->assembleStiffnessMatrix();
} else if (matrix_id == "M") {
this->assembleMass();
}
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::assembleLumpedMatrix(const ID & matrix_id) {
if (matrix_id == "M") {
this->assembleMassLumped();
}
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::beforeSolveStep() {
for (auto & material : materials) {
material->beforeSolveStep();
}
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::afterSolveStep(bool converged) {
for (auto & material : materials) {
material->afterSolveStep(converged);
}
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::predictor() { ++displacement_release; }
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::corrector() { ++displacement_release; }
/* -------------------------------------------------------------------------- */
/**
* This function computes the internal forces as \f$F_{int} = \int_{\Omega} N
* \sigma d\Omega@\f$
*/
void SolidMechanicsModel::assembleInternalForces() {
AKANTU_DEBUG_IN();
AKANTU_DEBUG_INFO("Assemble the internal forces");
this->internal_force->zero();
// compute the stresses of local elements
AKANTU_DEBUG_INFO("Compute local stresses");
for (auto & material : materials) {
material->computeAllStresses(_not_ghost);
}
/* ------------------------------------------------------------------------ */
/* Computation of the non local part */
if (this->non_local_manager) {
this->non_local_manager->computeAllNonLocalStresses();
}
// communicate the stresses
AKANTU_DEBUG_INFO("Send data for residual assembly");
this->asynchronousSynchronize(SynchronizationTag::_smm_stress);
// assemble the forces due to local stresses
AKANTU_DEBUG_INFO("Assemble residual for local elements");
for (auto & material : materials) {
material->assembleInternalForces(_not_ghost);
}
// finalize communications
AKANTU_DEBUG_INFO("Wait distant stresses");
this->waitEndSynchronize(SynchronizationTag::_smm_stress);
// assemble the stresses due to ghost elements
AKANTU_DEBUG_INFO("Assemble residual for ghost elements");
for (auto & material : materials) {
material->assembleInternalForces(_ghost);
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::assembleStiffnessMatrix(bool need_to_reassemble) {
AKANTU_DEBUG_IN();
AKANTU_DEBUG_INFO("Assemble the new stiffness matrix.");
// Check if materials need to recompute the matrix
/*bool need_to_reassemble = true;*/
//bool need_to_reassemble = false;
for (auto & material : materials) {
need_to_reassemble |= material->hasMatrixChanged("K");
}
if (need_to_reassemble) {
this->getDOFManager().getMatrix("K").zero();
// call compute stiffness matrix on each local elements
for (auto & material : materials) {
material->assembleStiffnessMatrix(_not_ghost);
}
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::updateCurrentPosition() {
if (this->current_position_release == this->displacement_release) {
return;
}
this->current_position->copy(this->mesh.getNodes());
auto cpos_it = this->current_position->begin(Model::spatial_dimension);
auto cpos_end = this->current_position->end(Model::spatial_dimension);
auto disp_it = this->displacement->begin(Model::spatial_dimension);
for (; cpos_it != cpos_end; ++cpos_it, ++disp_it) {
*cpos_it += *disp_it;
}
this->current_position_release = this->displacement_release;
}
/* -------------------------------------------------------------------------- */
const Array<Real> & SolidMechanicsModel::getCurrentPosition() {
this->updateCurrentPosition();
return *this->current_position;
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::updateDataForNonLocalCriterion(
ElementTypeMapReal & criterion) {
const ID field_name = criterion.getName();
for (auto & material : materials) {
if (!material->isInternal<Real>(field_name, _ek_regular)) {
continue;
}
for (auto ghost_type : ghost_types) {
material->flattenInternal(field_name, criterion, ghost_type, _ek_regular);
}
}
}
/* -------------------------------------------------------------------------- */
/* Information */
/* -------------------------------------------------------------------------- */
Real SolidMechanicsModel::getStableTimeStep() {
AKANTU_DEBUG_IN();
Real min_dt = getStableTimeStep(_not_ghost);
/// reduction min over all processors
mesh.getCommunicator().allReduce(min_dt, SynchronizerOperation::_min);
AKANTU_DEBUG_OUT();
return min_dt;
}
/* -------------------------------------------------------------------------- */
Real SolidMechanicsModel::getStableTimeStep(GhostType ghost_type) {
AKANTU_DEBUG_IN();
Real min_dt = std::numeric_limits<Real>::max();
this->updateCurrentPosition();
Element elem;
elem.ghost_type = ghost_type;
for (auto type :
mesh.elementTypes(Model::spatial_dimension, ghost_type, _ek_regular)) {
elem.type = type;
UInt nb_nodes_per_element = mesh.getNbNodesPerElement(type);
UInt nb_element = mesh.getNbElement(type);
auto mat_indexes = material_index(type, ghost_type).begin();
auto mat_loc_num = material_local_numbering(type, ghost_type).begin();
Array<Real> X(0, nb_nodes_per_element * Model::spatial_dimension);
FEEngine::extractNodalToElementField(mesh, *current_position, X, type,
_not_ghost);
auto X_el = X.begin(Model::spatial_dimension, nb_nodes_per_element);
for (UInt el = 0; el < nb_element;
++el, ++X_el, ++mat_indexes, ++mat_loc_num) {
elem.element = *mat_loc_num;
Real el_h = getFEEngine().getElementInradius(*X_el, type);
Real el_c = this->materials[*mat_indexes]->getCelerity(elem);
Real el_dt = el_h / el_c;
min_dt = std::min(min_dt, el_dt);
}
}
AKANTU_DEBUG_OUT();
return min_dt;
}
/* -------------------------------------------------------------------------- */
Real SolidMechanicsModel::getKineticEnergy() {
AKANTU_DEBUG_IN();
Real ekin = 0.;
UInt nb_nodes = mesh.getNbNodes();
if (this->getDOFManager().hasLumpedMatrix("M")) {
auto m_it = this->mass->begin(Model::spatial_dimension);
auto m_end = this->mass->end(Model::spatial_dimension);
auto v_it = this->velocity->begin(Model::spatial_dimension);
for (UInt n = 0; m_it != m_end; ++n, ++m_it, ++v_it) {
const auto & v = *v_it;
const auto & m = *m_it;
Real mv2 = 0.;
auto is_local_node = mesh.isLocalOrMasterNode(n);
// bool is_not_pbc_slave_node = !isPBCSlaveNode(n);
auto count_node = is_local_node; // && is_not_pbc_slave_node;
if (count_node) {
for (UInt i = 0; i < Model::spatial_dimension; ++i) {
if (m(i) > std::numeric_limits<Real>::epsilon()) {
mv2 += v(i) * v(i) * m(i);
}
}
}
ekin += mv2;
}
} else if (this->getDOFManager().hasMatrix("M")) {
Array<Real> Mv(nb_nodes, Model::spatial_dimension);
this->getDOFManager().assembleMatMulVectToArray("displacement", "M",
*this->velocity, Mv);
for (auto && data : zip(arange(nb_nodes), make_view(Mv, spatial_dimension),
make_view(*this->velocity, spatial_dimension))) {
ekin += std::get<2>(data).dot(std::get<1>(data)) *
static_cast<Real>(mesh.isLocalOrMasterNode(std::get<0>(data)));
}
} else {
AKANTU_ERROR("No function called to assemble the mass matrix.");
}
mesh.getCommunicator().allReduce(ekin, SynchronizerOperation::_sum);
AKANTU_DEBUG_OUT();
return ekin * .5;
}
/* -------------------------------------------------------------------------- */
Real SolidMechanicsModel::getKineticEnergy(ElementType type, UInt index) {
AKANTU_DEBUG_IN();
UInt nb_quadrature_points = getFEEngine().getNbIntegrationPoints(type);
Array<Real> vel_on_quad(nb_quadrature_points, Model::spatial_dimension);
Array<UInt> filter_element(1, 1, index);
getFEEngine().interpolateOnIntegrationPoints(*velocity, vel_on_quad,
Model::spatial_dimension, type,
_not_ghost, filter_element);
auto vit = vel_on_quad.begin(Model::spatial_dimension);
auto vend = vel_on_quad.end(Model::spatial_dimension);
Vector<Real> rho_v2(nb_quadrature_points);
Real rho = materials[material_index(type)(index)]->getRho();
for (UInt q = 0; vit != vend; ++vit, ++q) {
rho_v2(q) = rho * vit->dot(*vit);
}
AKANTU_DEBUG_OUT();
return .5 * getFEEngine().integrate(rho_v2, type, index);
}
/* -------------------------------------------------------------------------- */
Real SolidMechanicsModel::getExternalWork() {
AKANTU_DEBUG_IN();
auto ext_force_it = external_force->begin(Model::spatial_dimension);
auto int_force_it = internal_force->begin(Model::spatial_dimension);
auto boun_it = blocked_dofs->begin(Model::spatial_dimension);
decltype(ext_force_it) incr_or_velo_it;
if (this->method == _static) {
incr_or_velo_it =
this->displacement_increment->begin(Model::spatial_dimension);
} else {
incr_or_velo_it = this->velocity->begin(Model::spatial_dimension);
}
Real work = 0.;
UInt nb_nodes = this->mesh.getNbNodes();
for (UInt n = 0; n < nb_nodes;
++n, ++ext_force_it, ++int_force_it, ++boun_it, ++incr_or_velo_it) {
const auto & int_force = *int_force_it;
const auto & ext_force = *ext_force_it;
const auto & boun = *boun_it;
const auto & incr_or_velo = *incr_or_velo_it;
bool is_local_node = this->mesh.isLocalOrMasterNode(n);
// bool is_not_pbc_slave_node = !this->isPBCSlaveNode(n);
bool count_node = is_local_node; // && is_not_pbc_slave_node;
if (count_node) {
for (UInt i = 0; i < Model::spatial_dimension; ++i) {
if (boun(i)) {
work -= int_force(i) * incr_or_velo(i);
} else {
work += ext_force(i) * incr_or_velo(i);
}
}
}
}
mesh.getCommunicator().allReduce(work, SynchronizerOperation::_sum);
if (this->method != _static) {
work *= this->getTimeStep();
}
AKANTU_DEBUG_OUT();
return work;
}
/* -------------------------------------------------------------------------- */
Real SolidMechanicsModel::getEnergy(const std::string & energy_id) {
AKANTU_DEBUG_IN();
if (energy_id == "kinetic") {
return getKineticEnergy();
}
if (energy_id == "external work") {
return getExternalWork();
}
Real energy = 0.;
for (auto & material : materials) {
energy += material->getEnergy(energy_id);
}
/// reduction sum over all processors
mesh.getCommunicator().allReduce(energy, SynchronizerOperation::_sum);
AKANTU_DEBUG_OUT();
return energy;
}
/* -------------------------------------------------------------------------- */
Real SolidMechanicsModel::getEnergy(const std::string & energy_id,
ElementType type, UInt index) {
AKANTU_DEBUG_IN();
if (energy_id == "kinetic") {
return getKineticEnergy(type, index);
}
UInt mat_index = this->material_index(type, _not_ghost)(index);
UInt mat_loc_num = this->material_local_numbering(type, _not_ghost)(index);
Real energy =
this->materials[mat_index]->getEnergy(energy_id, type, mat_loc_num);
AKANTU_DEBUG_OUT();
return energy;
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::onElementsAdded(const Array<Element> & element_list,
const NewElementsEvent & event) {
AKANTU_DEBUG_IN();
this->material_index.initialize(mesh, _element_kind = _ek_not_defined,
_with_nb_element = true,
_default_value = UInt(-1));
this->material_local_numbering.initialize(
mesh, _element_kind = _ek_not_defined, _with_nb_element = true,
_default_value = UInt(-1));
ElementTypeMapArray<UInt> filter("new_element_filter", this->getID(),
this->getMemoryID());
for (const auto & elem : element_list) {
if (mesh.getSpatialDimension(elem.type) != spatial_dimension) {
continue;
}
if (!filter.exists(elem.type, elem.ghost_type)) {
filter.alloc(0, 1, elem.type, elem.ghost_type);
}
filter(elem.type, elem.ghost_type).push_back(elem.element);
}
// this fails in parallel if the event is sent on facet between constructor
// and initFull \todo: to debug...
this->assignMaterialToElements(&filter);
for (auto & material : materials) {
material->onElementsAdded(element_list, event);
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::onElementsRemoved(
const Array<Element> & element_list,
const ElementTypeMapArray<UInt> & new_numbering,
const RemovedElementsEvent & event) {
for (auto & material : materials) {
material->onElementsRemoved(element_list, new_numbering, event);
}
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::onNodesAdded(const Array<UInt> & nodes_list,
const NewNodesEvent & event) {
AKANTU_DEBUG_IN();
UInt nb_nodes = mesh.getNbNodes();
if (displacement) {
displacement->resize(nb_nodes, 0.);
++displacement_release;
}
if (mass) {
mass->resize(nb_nodes, 0.);
}
if (velocity) {
velocity->resize(nb_nodes, 0.);
}
if (acceleration) {
acceleration->resize(nb_nodes, 0.);
}
if (external_force) {
external_force->resize(nb_nodes, 0.);
}
if (internal_force) {
internal_force->resize(nb_nodes, 0.);
}
if (blocked_dofs) {
blocked_dofs->resize(nb_nodes, false);
}
if (current_position) {
current_position->resize(nb_nodes, 0.);
}
if (previous_displacement) {
previous_displacement->resize(nb_nodes, 0.);
}
if (displacement_increment) {
displacement_increment->resize(nb_nodes, 0.);
}
for (auto & material : materials) {
material->onNodesAdded(nodes_list, event);
}
need_to_reassemble_lumped_mass = true;
need_to_reassemble_mass = true;
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::onNodesRemoved(const Array<UInt> & /*element_list*/,
const Array<UInt> & new_numbering,
const RemovedNodesEvent & /*event*/) {
if (displacement) {
mesh.removeNodesFromArray(*displacement, new_numbering);
++displacement_release;
}
if (mass) {
mesh.removeNodesFromArray(*mass, new_numbering);
}
if (velocity) {
mesh.removeNodesFromArray(*velocity, new_numbering);
}
if (acceleration) {
mesh.removeNodesFromArray(*acceleration, new_numbering);
}
if (internal_force) {
mesh.removeNodesFromArray(*internal_force, new_numbering);
}
if (external_force) {
mesh.removeNodesFromArray(*external_force, new_numbering);
}
if (blocked_dofs) {
mesh.removeNodesFromArray(*blocked_dofs, new_numbering);
}
// if (increment_acceleration)
// mesh.removeNodesFromArray(*increment_acceleration, new_numbering);
if (displacement_increment) {
mesh.removeNodesFromArray(*displacement_increment, new_numbering);
}
if (previous_displacement) {
mesh.removeNodesFromArray(*previous_displacement, new_numbering);
}
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::printself(std::ostream & stream, int indent) const {
std::string space(indent, AKANTU_INDENT);
stream << space << "Solid Mechanics Model [" << std::endl;
stream << space << " + id : " << id << std::endl;
stream << space << " + spatial dimension : " << Model::spatial_dimension
<< std::endl;
stream << space << " + fem [" << std::endl;
getFEEngine().printself(stream, indent + 2);
stream << space << " ]" << std::endl;
stream << space << " + nodals information [" << std::endl;
displacement->printself(stream, indent + 2);
if (velocity) {
velocity->printself(stream, indent + 2);
}
if (acceleration) {
acceleration->printself(stream, indent + 2);
}
if (mass) {
mass->printself(stream, indent + 2);
}
external_force->printself(stream, indent + 2);
internal_force->printself(stream, indent + 2);
blocked_dofs->printself(stream, indent + 2);
stream << space << " ]" << std::endl;
stream << space << " + material information [" << std::endl;
material_index.printself(stream, indent + 2);
stream << space << " ]" << std::endl;
stream << space << " + materials [" << std::endl;
for (const auto & material : materials) {
material->printself(stream, indent + 2);
}
stream << space << " ]" << std::endl;
stream << space << "]" << std::endl;
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::initializeNonLocal() {
this->non_local_manager->synchronize(*this, SynchronizationTag::_material_id);
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::insertIntegrationPointsInNeighborhoods(
GhostType ghost_type) {
for (auto & mat : materials) {
MaterialNonLocalInterface * mat_non_local;
if ((mat_non_local =
dynamic_cast<MaterialNonLocalInterface *>(mat.get())) == nullptr) {
continue;
}
ElementTypeMapArray<Real> quadrature_points_coordinates(
"quadrature_points_coordinates_tmp_nl", this->id, this->memory_id);
quadrature_points_coordinates.initialize(this->getFEEngine(),
_nb_component = spatial_dimension,
_ghost_type = ghost_type);
for (const auto & type : quadrature_points_coordinates.elementTypes(
Model::spatial_dimension, ghost_type)) {
this->getFEEngine().computeIntegrationPointsCoordinates(
quadrature_points_coordinates(type, ghost_type), type, ghost_type);
}
mat_non_local->initMaterialNonLocal();
mat_non_local->insertIntegrationPointsInNeighborhoods(
ghost_type, quadrature_points_coordinates);
}
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::computeNonLocalStresses(GhostType ghost_type) {
for (auto & mat : materials) {
if (not aka::is_of_type<MaterialNonLocalInterface>(*mat)) {
continue;
}
auto & mat_non_local = dynamic_cast<MaterialNonLocalInterface &>(*mat);
mat_non_local.computeNonLocalStresses(ghost_type);
}
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::updateLocalInternal(
ElementTypeMapReal & internal_flat, GhostType ghost_type,
ElementKind kind) {
const ID field_name = internal_flat.getName();
for (auto & material : materials) {
if (material->isInternal<Real>(field_name, kind)) {
material->flattenInternal(field_name, internal_flat, ghost_type, kind);
}
}
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::updateNonLocalInternal(
ElementTypeMapReal & internal_flat, GhostType ghost_type,
ElementKind kind) {
const ID field_name = internal_flat.getName();
for (auto & mat : materials) {
if (not aka::is_of_type<MaterialNonLocalInterface>(*mat)) {
continue;
}
auto & mat_non_local = dynamic_cast<MaterialNonLocalInterface &>(*mat);
mat_non_local.updateNonLocalInternals(internal_flat, field_name, ghost_type,
kind);
}
}
/* -------------------------------------------------------------------------- */
FEEngine & SolidMechanicsModel::getFEEngineBoundary(const ID & name) {
return getFEEngineClassBoundary<MyFEEngineType>(name);
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::splitElementByMaterial(
const Array<Element> & elements,
std::vector<Array<Element>> & elements_per_mat) const {
for (const auto & el : elements) {
Element mat_el = el;
mat_el.element = this->material_local_numbering(el);
elements_per_mat[this->material_index(el)].push_back(mat_el);
}
}
/* -------------------------------------------------------------------------- */
UInt SolidMechanicsModel::getNbData(const Array<Element> & elements,
const SynchronizationTag & tag) const {
AKANTU_DEBUG_IN();
UInt size = 0;
UInt nb_nodes_per_element = 0;
for (const Element & el : elements) {
nb_nodes_per_element += Mesh::getNbNodesPerElement(el.type);
}
switch (tag) {
case SynchronizationTag::_material_id: {
size += elements.size() * sizeof(UInt);
break;
}
case SynchronizationTag::_smm_mass: {
size += nb_nodes_per_element * sizeof(Real) *
Model::spatial_dimension; // mass vector
break;
}
case SynchronizationTag::_smm_for_gradu: {
size += nb_nodes_per_element * Model::spatial_dimension *
sizeof(Real); // displacement
break;
}
case SynchronizationTag::_smm_boundary: {
// force, displacement, boundary
size += nb_nodes_per_element * Model::spatial_dimension *
(2 * sizeof(Real) + sizeof(bool));
break;
}
case SynchronizationTag::_for_dump: {
// displacement, velocity, acceleration, residual, force
size += nb_nodes_per_element * Model::spatial_dimension * sizeof(Real) * 5;
break;
}
default: {
}
}
if (tag != SynchronizationTag::_material_id) {
splitByMaterial(elements, [&](auto && mat, auto && elements) {
size += mat.getNbData(elements, tag);
});
}
AKANTU_DEBUG_OUT();
return size;
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::packData(CommunicationBuffer & buffer,
const Array<Element> & elements,
const SynchronizationTag & tag) const {
AKANTU_DEBUG_IN();
switch (tag) {
case SynchronizationTag::_material_id: {
packElementalDataHelper(
material_index, buffer, elements, false, getFEEngine());
break;
}
case SynchronizationTag::_smm_mass: {
packNodalDataHelper(*mass, buffer, elements, mesh);
break;
}
case SynchronizationTag::_smm_for_gradu: {
packNodalDataHelper(*displacement, buffer, elements, mesh);
break;
}
case SynchronizationTag::_for_dump: {
packNodalDataHelper(*displacement, buffer, elements, mesh);
packNodalDataHelper(*velocity, buffer, elements, mesh);
packNodalDataHelper(*acceleration, buffer, elements, mesh);
packNodalDataHelper(*internal_force, buffer, elements, mesh);
packNodalDataHelper(*external_force, buffer, elements, mesh);
break;
}
case SynchronizationTag::_smm_boundary: {
packNodalDataHelper(*external_force, buffer, elements, mesh);
packNodalDataHelper(*velocity, buffer, elements, mesh);
packNodalDataHelper(*blocked_dofs, buffer, elements, mesh);
break;
}
default: {
}
}
if (tag != SynchronizationTag::_material_id) {
splitByMaterial(elements, [&](auto && mat, auto && elements) {
mat.packData(buffer, elements, tag);
});
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::unpackData(CommunicationBuffer & buffer,
const Array<Element> & elements,
const SynchronizationTag & tag) {
AKANTU_DEBUG_IN();
switch (tag) {
case SynchronizationTag::_material_id: {
for (auto && element : elements) {
UInt recv_mat_index;
buffer >> recv_mat_index;
UInt & mat_index = material_index(element);
if (mat_index != UInt(-1)) {
continue;
}
// add ghosts element to the correct material
mat_index = recv_mat_index;
UInt index = materials[mat_index]->addElement(element);
material_local_numbering(element) = index;
}
break;
}
case SynchronizationTag::_smm_mass: {
unpackNodalDataHelper(*mass, buffer, elements, mesh);
break;
}
case SynchronizationTag::_smm_for_gradu: {
unpackNodalDataHelper(*displacement, buffer, elements, mesh);
break;
}
case SynchronizationTag::_for_dump: {
unpackNodalDataHelper(*displacement, buffer, elements, mesh);
unpackNodalDataHelper(*velocity, buffer, elements, mesh);
unpackNodalDataHelper(*acceleration, buffer, elements, mesh);
unpackNodalDataHelper(*internal_force, buffer, elements, mesh);
unpackNodalDataHelper(*external_force, buffer, elements, mesh);
break;
}
case SynchronizationTag::_smm_boundary: {
unpackNodalDataHelper(*external_force, buffer, elements, mesh);
unpackNodalDataHelper(*velocity, buffer, elements, mesh);
unpackNodalDataHelper(*blocked_dofs, buffer, elements, mesh);
break;
}
default: {
}
}
if (tag != SynchronizationTag::_material_id) {
splitByMaterial(elements, [&](auto && mat, auto && elements) {
mat.unpackData(buffer, elements, tag);
});
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
UInt SolidMechanicsModel::getNbData(const Array<UInt> & dofs,
const SynchronizationTag & tag) const {
AKANTU_DEBUG_IN();
UInt size = 0;
// UInt nb_nodes = mesh.getNbNodes();
switch (tag) {
case SynchronizationTag::_smm_uv: {
size += sizeof(Real) * Model::spatial_dimension * 2;
break;
}
case SynchronizationTag::_smm_res: /* FALLTHRU */
case SynchronizationTag::_smm_mass: {
size += sizeof(Real) * Model::spatial_dimension;
break;
}
case SynchronizationTag::_for_dump: {
size += sizeof(Real) * Model::spatial_dimension * 5;
break;
}
default: {
AKANTU_ERROR("Unknown ghost synchronization tag : " << tag);
}
}
AKANTU_DEBUG_OUT();
return size * dofs.size();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::packData(CommunicationBuffer & buffer,
const Array<UInt> & dofs,
const SynchronizationTag & tag) const {
AKANTU_DEBUG_IN();
switch (tag) {
case SynchronizationTag::_smm_uv: {
packDOFDataHelper(*displacement, buffer, dofs);
packDOFDataHelper(*velocity, buffer, dofs);
break;
}
case SynchronizationTag::_smm_res: {
packDOFDataHelper(*internal_force, buffer, dofs);
break;
}
case SynchronizationTag::_smm_mass: {
packDOFDataHelper(*mass, buffer, dofs);
break;
}
case SynchronizationTag::_for_dump: {
packDOFDataHelper(*displacement, buffer, dofs);
packDOFDataHelper(*velocity, buffer, dofs);
packDOFDataHelper(*acceleration, buffer, dofs);
packDOFDataHelper(*internal_force, buffer, dofs);
packDOFDataHelper(*external_force, buffer, dofs);
break;
}
default: {
AKANTU_ERROR("Unknown ghost synchronization tag : " << tag);
}
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::unpackData(CommunicationBuffer & buffer,
const Array<UInt> & dofs,
const SynchronizationTag & tag) {
AKANTU_DEBUG_IN();
switch (tag) {
case SynchronizationTag::_smm_uv: {
unpackDOFDataHelper(*displacement, buffer, dofs);
unpackDOFDataHelper(*velocity, buffer, dofs);
break;
}
case SynchronizationTag::_smm_res: {
unpackDOFDataHelper(*internal_force, buffer, dofs);
break;
}
case SynchronizationTag::_smm_mass: {
unpackDOFDataHelper(*mass, buffer, dofs);
break;
}
case SynchronizationTag::_for_dump: {
unpackDOFDataHelper(*displacement, buffer, dofs);
unpackDOFDataHelper(*velocity, buffer, dofs);
unpackDOFDataHelper(*acceleration, buffer, dofs);
unpackDOFDataHelper(*internal_force, buffer, dofs);
unpackDOFDataHelper(*external_force, buffer, dofs);
break;
}
default: {
AKANTU_ERROR("Unknown ghost synchronization tag : " << tag);
}
}
AKANTU_DEBUG_OUT();
}
} // namespace akantu