diff --git a/python/py_mesh.cc b/python/py_mesh.cc
index f30ba8505..e06154338 100644
--- a/python/py_mesh.cc
+++ b/python/py_mesh.cc
@@ -1,200 +1,207 @@
/**
* @file py_mesh.cc
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
* @author Philip Mueller <philip.paul.mueller@bluemail.ch>
* @author Mohit Pundir <mohit.pundir@epfl.ch>
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @date creation: Sun Jun 16 2019
* @date last modification: Mon Mar 15 2021
*
* @brief pybind11 interface to Mesh
*
*
* @section LICENSE
*
* Copyright (©) 2018-2021 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 "aka_config.hh"
/* -------------------------------------------------------------------------- */
#include "py_aka_array.hh"
/* -------------------------------------------------------------------------- */
#include <mesh.hh>
#include <mesh_accessor.hh>
#include <mesh_utils.hh>
/* -------------------------------------------------------------------------- */
#include <pybind11/pybind11.h>
#include <pybind11/stl.h>
/* -------------------------------------------------------------------------- */
namespace py = pybind11;
/* -------------------------------------------------------------------------- */
namespace akantu {
namespace {
/* ------------------------------------------------------------------------ */
template <typename T>
void register_element_type_map_array(py::module & mod,
const std::string & name) {
py::class_<ElementTypeMapArray<T>, std::shared_ptr<ElementTypeMapArray<T>>>(
mod, ("ElementTypeMapArray" + name).c_str())
.def(
"__call__",
[](ElementTypeMapArray<T> & self, ElementType type,
GhostType ghost_type) -> decltype(auto) {
return self(type, ghost_type);
},
py::arg("type"), py::arg("ghost_type") = _not_ghost,
py::return_value_policy::reference)
.def(
"elementTypes",
[](ElementTypeMapArray<T> & self, UInt _dim, GhostType _ghost_type,
ElementKind _kind) -> std::vector<ElementType> {
auto types = self.elementTypes(_dim, _ghost_type, _kind);
std::vector<ElementType> _types;
for (auto && t : types) {
_types.push_back(t);
}
return _types;
},
py::arg("dim") = _all_dimensions,
py::arg("ghost_type") = _not_ghost, py::arg("kind") = _ek_regular);
}
} // namespace
/* -------------------------------------------------------------------------- */
void register_mesh(py::module & mod) {
register_element_type_map_array<Real>(mod, "Real");
register_element_type_map_array<UInt>(mod, "UInt");
// register_element_type_map_array<std::string>(mod, "String");
+ py::class_<Mesh::PeriodicSlaves>(mod, "PeriodicSlaves")
+ .def(
+ "__iter__",
+ [](Mesh::PeriodicSlaves & _this) {
+ py::make_iterator(_this.begin(), _this.end());
+ },
+ py::keep_alive<0, 1>());
+
py::class_<MeshData>(mod, "MeshData")
.def(
"getElementalDataUInt",
[](MeshData & _this, const ID & name) -> decltype(auto) {
return _this.getElementalData<UInt>(name);
},
py::return_value_policy::reference)
.def(
"getElementalDataReal",
[](MeshData & _this, const ID & name) -> decltype(auto) {
return _this.getElementalData<Real>(name);
},
py::return_value_policy::reference);
py::class_<Mesh, GroupManager, Dumpable, MeshData>(mod, "Mesh",
py::multiple_inheritance())
.def(py::init<UInt, const ID &>(), py::arg("spatial_dimension"),
py::arg("id") = "mesh")
.def("read", &Mesh::read, py::arg("filename"),
py::arg("mesh_io_type") = _miot_auto, "read the mesh from a file")
.def(
"getNodes",
[](Mesh & self) -> decltype(auto) { return self.getNodes(); },
py::return_value_policy::reference)
.def("getNbNodes", &Mesh::getNbNodes)
.def(
"getConnectivity",
[](Mesh & self, ElementType type) -> decltype(auto) {
return self.getConnectivity(type);
},
py::return_value_policy::reference)
.def(
"addConnectivityType",
[](Mesh & self, ElementType type, GhostType ghost_type) -> void {
self.addConnectivityType(type, ghost_type);
},
py::arg("type"), py::arg("ghost_type") = _not_ghost)
.def("distribute", [](Mesh & self) { self.distribute(); })
.def("fillNodesToElements", &Mesh::fillNodesToElements,
py::arg("dimension") = _all_dimensions)
.def("getAssociatedElements",
[](Mesh & self, const UInt & node, py::list list) {
Array<Element> elements;
self.getAssociatedElements(node, elements);
for (auto && element : elements) {
list.append(element);
}
})
.def("makePeriodic",
[](Mesh & self, const SpatialDirection & direction) {
self.makePeriodic(direction);
})
.def(
"getNbElement",
[](Mesh & self, const UInt spatial_dimension, GhostType ghost_type,
ElementKind kind) {
return self.getNbElement(spatial_dimension, ghost_type, kind);
},
py::arg("spatial_dimension") = _all_dimensions,
py::arg("ghost_type") = _not_ghost, py::arg("kind") = _ek_not_defined)
.def(
"getNbElement",
[](Mesh & self, ElementType type, GhostType ghost_type) {
return self.getNbElement(type, ghost_type);
},
py::arg("type"), py::arg("ghost_type") = _not_ghost)
.def_static(
"getSpatialDimension",
[](ElementType & type) { return Mesh::getSpatialDimension(type); })
.def(
"getDataReal",
[](Mesh & _this, const ID & name, ElementType type,
GhostType ghost_type) -> decltype(auto) {
return _this.getData<Real>(name, type, ghost_type);
},
py::arg("name"), py::arg("type"), py::arg("ghost_type") = _not_ghost,
py::return_value_policy::reference)
.def(
"hasDataReal",
[](Mesh & _this, const ID & name, ElementType type,
GhostType ghost_type) -> bool {
return _this.hasData<Real>(name, type, ghost_type);
},
py::arg("name"), py::arg("type"), py::arg("ghost_type") = _not_ghost)
.def("isPeriodic", [](const Mesh & _this) { return _this.isPeriodic(); })
.def("getPeriodicMaster", &Mesh::getPeriodicMaster)
.def("getPeriodicSlaves", &Mesh::getPeriodicSlaves)
.def("isPeriodicSlave", &Mesh::isPeriodicSlave)
- .def("isPeriodicMaster", &Mesh::isPeriodicMaster)
- .def("getPeriodicMasterSlaves", &Mesh::getPeriodicMasterSlaves);
+ .def("isPeriodicMaster", &Mesh::isPeriodicMaster);
/* ------------------------------------------------------------------------ */
py::class_<MeshUtils>(mod, "MeshUtils")
.def_static("buildFacets", &MeshUtils::buildFacets);
py::class_<MeshAccessor>(mod, "MeshAccessor")
.def(py::init<Mesh &>(), py::arg("mesh"))
.def(
"resizeConnectivity",
[](MeshAccessor & self, UInt new_size, ElementType type, GhostType gt)
-> void { self.resizeConnectivity(new_size, type, gt); },
py::arg("new_size"), py::arg("type"),
py::arg("ghost_type") = _not_ghost)
.def(
"resizeNodes",
[](MeshAccessor & self, UInt new_size) -> void {
self.resizeNodes(new_size);
},
py::arg("new_size"))
.def("makeReady", &MeshAccessor::makeReady);
}
} // namespace akantu
diff --git a/src/mesh/mesh_inline_impl.hh b/src/mesh/mesh_inline_impl.hh
index 7b5e3aff2..22d27cfd4 100644
--- a/src/mesh/mesh_inline_impl.hh
+++ b/src/mesh/mesh_inline_impl.hh
@@ -1,764 +1,765 @@
/**
* @file mesh_inline_impl.hh
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
* @author Dana Christen <dana.christen@epfl.ch>
* @author Mohit Pundir <mohit.pundir@epfl.ch>
* @author Nicolas Richart <nicolas.richart@epfl.ch>
* @author Marco Vocialta <marco.vocialta@epfl.ch>
*
* @date creation: Thu Jul 15 2010
* @date last modification: Fri Dec 11 2020
*
* @brief Implementation of the inline functions of the mesh class
*
*
* @section LICENSE
*
* Copyright (©) 2010-2021 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 "aka_iterators.hh"
#include "element_class.hh"
#include "mesh.hh"
/* -------------------------------------------------------------------------- */
#ifndef AKANTU_MESH_INLINE_IMPL_HH_
#define AKANTU_MESH_INLINE_IMPL_HH_
namespace akantu {
/* -------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
inline ElementKind Element::kind() const { return Mesh::getKind(type); }
/* -------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
template <typename... pack>
Mesh::ElementTypesIteratorHelper Mesh::elementTypes(pack &&... _pack) const {
return connectivities.elementTypes(_pack...);
}
/* -------------------------------------------------------------------------- */
inline RemovedNodesEvent::RemovedNodesEvent(const Mesh & mesh,
const std::string & origin)
: MeshEvent<UInt>(origin),
new_numbering(mesh.getNbNodes(), 1, "new_numbering") {}
/* -------------------------------------------------------------------------- */
inline RemovedElementsEvent::RemovedElementsEvent(const Mesh & mesh,
const ID & new_numbering_id,
const std::string & origin)
: MeshEvent<Element>(origin),
new_numbering(new_numbering_id, mesh.getID()) {}
/* -------------------------------------------------------------------------- */
template <>
inline void Mesh::sendEvent<NewElementsEvent>(NewElementsEvent & event) {
this->fillNodesToElements();
EventHandlerManager<MeshEventHandler>::sendEvent(event);
}
/* -------------------------------------------------------------------------- */
template <> inline void Mesh::sendEvent<NewNodesEvent>(NewNodesEvent & event) {
this->computeBoundingBox();
this->nodes_flags->resize(this->nodes->size(), NodeFlag::_normal);
EventHandlerManager<MeshEventHandler>::sendEvent(event);
}
/* -------------------------------------------------------------------------- */
template <>
inline void
Mesh::sendEvent<RemovedElementsEvent>(RemovedElementsEvent & event) {
this->connectivities.onElementsRemoved(event.getNewNumbering());
this->fillNodesToElements();
this->computeBoundingBox();
EventHandlerManager<MeshEventHandler>::sendEvent(event);
}
/* -------------------------------------------------------------------------- */
template <>
inline void Mesh::sendEvent<RemovedNodesEvent>(RemovedNodesEvent & event) {
const auto & new_numbering = event.getNewNumbering();
this->removeNodesFromArray(*nodes, new_numbering);
if (nodes_global_ids and not is_mesh_facets) {
this->removeNodesFromArray(*nodes_global_ids, new_numbering);
}
if (not is_mesh_facets) {
this->removeNodesFromArray(*nodes_flags, new_numbering);
}
if (not nodes_to_elements.empty()) {
std::vector<std::unique_ptr<std::set<Element>>> tmp(
nodes_to_elements.size());
auto it = nodes_to_elements.begin();
UInt new_nb_nodes = 0;
for (auto new_i : new_numbering) {
if (new_i != UInt(-1)) {
tmp[new_i] = std::move(*it);
++new_nb_nodes;
}
++it;
}
tmp.resize(new_nb_nodes);
std::move(tmp.begin(), tmp.end(), nodes_to_elements.begin());
}
computeBoundingBox();
EventHandlerManager<MeshEventHandler>::sendEvent(event);
}
/* -------------------------------------------------------------------------- */
template <typename T>
inline void Mesh::removeNodesFromArray(Array<T> & vect,
const Array<UInt> & new_numbering) {
Array<T> tmp(vect.size(), vect.getNbComponent());
UInt nb_component = vect.getNbComponent();
UInt new_nb_nodes = 0;
for (UInt i = 0; i < new_numbering.size(); ++i) {
UInt new_i = new_numbering(i);
if (new_i != UInt(-1)) {
T * to_copy = vect.storage() + i * nb_component;
std::uninitialized_copy(to_copy, to_copy + nb_component,
tmp.storage() + new_i * nb_component);
++new_nb_nodes;
}
}
tmp.resize(new_nb_nodes);
vect.copy(tmp);
}
/* -------------------------------------------------------------------------- */
inline Array<UInt> & Mesh::getNodesGlobalIdsPointer() {
AKANTU_DEBUG_IN();
if (not nodes_global_ids) {
nodes_global_ids = std::make_shared<Array<UInt>>(
nodes->size(), 1, getID() + ":nodes_global_ids");
for (auto && global_ids : enumerate(*nodes_global_ids)) {
std::get<1>(global_ids) = std::get<0>(global_ids);
}
}
AKANTU_DEBUG_OUT();
return *nodes_global_ids;
}
/* -------------------------------------------------------------------------- */
inline Array<UInt> & Mesh::getConnectivityPointer(ElementType type,
GhostType ghost_type) {
if (connectivities.exists(type, ghost_type)) {
return connectivities(type, ghost_type);
}
if (ghost_type != _not_ghost) {
ghosts_counters.alloc(0, 1, type, ghost_type, 1);
}
AKANTU_DEBUG_INFO("The connectivity vector for the type " << type
<< " created");
UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type);
return connectivities.alloc(0, nb_nodes_per_element, type, ghost_type);
}
/* -------------------------------------------------------------------------- */
inline Array<std::vector<Element>> &
Mesh::getElementToSubelementPointer(ElementType type, GhostType ghost_type) {
return getDataPointer<std::vector<Element>>("element_to_subelement", type,
ghost_type, 1, true);
}
/* -------------------------------------------------------------------------- */
inline Array<Element> &
Mesh::getSubelementToElementPointer(ElementType type, GhostType ghost_type) {
auto & array = getDataPointer<Element>(
"subelement_to_element", type, ghost_type, getNbFacetsPerElement(type),
false, is_mesh_facets, ElementNull);
return array;
}
/* -------------------------------------------------------------------------- */
inline const auto & Mesh::getElementToSubelement() const {
return getData<std::vector<Element>>("element_to_subelement");
}
/* -------------------------------------------------------------------------- */
inline auto & Mesh::getElementToSubelementNC() {
return getData<std::vector<Element>>("element_to_subelement");
}
/* -------------------------------------------------------------------------- */
inline const auto & Mesh::getElementToSubelement(ElementType type,
GhostType ghost_type) const {
return getData<std::vector<Element>>("element_to_subelement", type,
ghost_type);
}
/* -------------------------------------------------------------------------- */
inline auto & Mesh::getElementToSubelementNC(ElementType type,
GhostType ghost_type) {
return getData<std::vector<Element>>("element_to_subelement", type,
ghost_type);
}
/* -------------------------------------------------------------------------- */
inline const auto &
Mesh::getElementToSubelement(const Element & element) const {
return getData<std::vector<Element>>("element_to_subelement")(element, 0);
}
/* -------------------------------------------------------------------------- */
inline auto & Mesh::getElementToSubelementNC(const Element & element) {
return getData<std::vector<Element>>("element_to_subelement")(element, 0);
}
/* -------------------------------------------------------------------------- */
inline const auto & Mesh::getSubelementToElement() const {
return getData<Element>("subelement_to_element");
}
/* -------------------------------------------------------------------------- */
inline auto & Mesh::getSubelementToElementNC() {
return getData<Element>("subelement_to_element");
}
/* -------------------------------------------------------------------------- */
inline const auto & Mesh::getSubelementToElement(ElementType type,
GhostType ghost_type) const {
return getData<Element>("subelement_to_element", type, ghost_type);
}
/* -------------------------------------------------------------------------- */
inline auto & Mesh::getSubelementToElementNC(ElementType type,
GhostType ghost_type) {
return getData<Element>("subelement_to_element", type, ghost_type);
}
/* -------------------------------------------------------------------------- */
inline VectorProxy<Element>
Mesh::getSubelementToElement(const Element & element) const {
return this->getSubelementToElement().get(element);
}
/* -------------------------------------------------------------------------- */
inline VectorProxy<Element>
Mesh::getSubelementToElementNC(const Element & element) const {
return this->getSubelementToElement().get(element);
}
/* -------------------------------------------------------------------------- */
template <typename T>
inline Array<T> &
Mesh::getDataPointer(const ID & data_name, ElementType el_type,
GhostType ghost_type, UInt nb_component,
bool size_to_nb_element, bool resize_with_parent) {
Array<T> & tmp = this->getElementalDataArrayAlloc<T>(
data_name, el_type, ghost_type, nb_component);
if (size_to_nb_element) {
if (resize_with_parent) {
tmp.resize(mesh_parent->getNbElement(el_type, ghost_type));
} else {
tmp.resize(this->getNbElement(el_type, ghost_type));
}
}
return tmp;
}
/* -------------------------------------------------------------------------- */
template <typename T>
inline Array<T> &
Mesh::getDataPointer(const ID & data_name, ElementType el_type,
GhostType ghost_type, UInt nb_component,
bool size_to_nb_element, bool resize_with_parent,
const T & defaul_) {
Array<T> & tmp = this->getElementalDataArrayAlloc<T>(
data_name, el_type, ghost_type, nb_component);
if (size_to_nb_element) {
if (resize_with_parent) {
tmp.resize(mesh_parent->getNbElement(el_type, ghost_type), defaul_);
} else {
tmp.resize(this->getNbElement(el_type, ghost_type), defaul_);
}
}
return tmp;
}
/* -------------------------------------------------------------------------- */
template <typename T>
inline const Array<T> & Mesh::getData(const ID & data_name, ElementType el_type,
GhostType ghost_type) const {
return this->getElementalDataArray<T>(data_name, el_type, ghost_type);
}
/* -------------------------------------------------------------------------- */
template <typename T>
inline Array<T> & Mesh::getData(const ID & data_name, ElementType el_type,
GhostType ghost_type) {
return this->getElementalDataArray<T>(data_name, el_type, ghost_type);
}
/* -------------------------------------------------------------------------- */
template <typename T>
inline const ElementTypeMapArray<T> &
Mesh::getData(const ID & data_name) const {
return this->getElementalData<T>(data_name);
}
/* -------------------------------------------------------------------------- */
template <typename T>
inline ElementTypeMapArray<T> & Mesh::getData(const ID & data_name) {
return this->getElementalData<T>(data_name);
}
/* -------------------------------------------------------------------------- */
inline UInt Mesh::getNbElement(ElementType type, GhostType ghost_type) const {
try {
const Array<UInt> & conn = connectivities(type, ghost_type);
return conn.size();
} catch (...) {
return 0;
}
}
/* -------------------------------------------------------------------------- */
inline UInt Mesh::getNbElement(const UInt spatial_dimension,
GhostType ghost_type, ElementKind kind) const {
AKANTU_DEBUG_ASSERT(spatial_dimension <= 3 || spatial_dimension == UInt(-1),
"spatial_dimension is " << spatial_dimension
<< " and is greater than 3 !");
UInt nb_element = 0;
for (auto type : elementTypes(spatial_dimension, ghost_type, kind)) {
nb_element += getNbElement(type, ghost_type);
}
return nb_element;
}
/* -------------------------------------------------------------------------- */
inline void Mesh::getBarycenter(const Element & element,
Vector<Real> & barycenter) const {
Vector<UInt> conn = getConnectivity(element);
Matrix<Real> local_coord(spatial_dimension, conn.size());
auto node_begin = make_view(*nodes, spatial_dimension).begin();
for (auto && node : enumerate(conn)) {
local_coord(std::get<0>(node)) =
Vector<Real>(node_begin[std::get<1>(node)]);
}
Math::barycenter(local_coord.storage(), conn.size(), spatial_dimension,
barycenter.storage());
}
/* -------------------------------------------------------------------------- */
inline UInt Mesh::getNbNodesPerElement(ElementType type) {
UInt nb_nodes_per_element = 0;
#define GET_NB_NODES_PER_ELEMENT(type) \
nb_nodes_per_element = ElementClass<type>::getNbNodesPerElement()
AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_NB_NODES_PER_ELEMENT);
#undef GET_NB_NODES_PER_ELEMENT
return nb_nodes_per_element;
}
/* -------------------------------------------------------------------------- */
inline ElementType Mesh::getP1ElementType(ElementType type) {
ElementType p1_type = _not_defined;
#define GET_P1_TYPE(type) p1_type = ElementClass<type>::getP1ElementType()
AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_P1_TYPE);
#undef GET_P1_TYPE
return p1_type;
}
/* -------------------------------------------------------------------------- */
inline ElementKind Mesh::getKind(ElementType type) {
ElementKind kind = _ek_not_defined;
#define GET_KIND(type) kind = ElementClass<type>::getKind()
AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_KIND);
#undef GET_KIND
return kind;
}
/* -------------------------------------------------------------------------- */
inline UInt Mesh::getSpatialDimension(ElementType type) {
UInt spatial_dimension = 0;
#define GET_SPATIAL_DIMENSION(type) \
spatial_dimension = ElementClass<type>::getSpatialDimension()
AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_SPATIAL_DIMENSION);
#undef GET_SPATIAL_DIMENSION
return spatial_dimension;
}
/* -------------------------------------------------------------------------- */
inline UInt Mesh::getNaturalSpaceDimension(const ElementType & type) {
UInt natural_dimension = 0;
#define GET_NATURAL_DIMENSION(type) \
natural_dimension = ElementClass<type>::getNaturalSpaceDimension()
AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_NATURAL_DIMENSION);
#undef GET_NATURAL_DIMENSION
return natural_dimension;
}
/* -------------------------------------------------------------------------- */
inline UInt Mesh::getNbFacetTypes(ElementType type, UInt /*t*/) {
UInt nb = 0;
#define GET_NB_FACET_TYPE(type) nb = ElementClass<type>::getNbFacetTypes()
AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_NB_FACET_TYPE);
#undef GET_NB_FACET_TYPE
return nb;
}
/* -------------------------------------------------------------------------- */
inline constexpr auto Mesh::getFacetType(ElementType type, UInt t) {
#define GET_FACET_TYPE(type) return ElementClass<type>::getFacetType(t);
AKANTU_BOOST_ALL_ELEMENT_SWITCH_NO_DEFAULT(GET_FACET_TYPE);
#undef GET_FACET_TYPE
return _not_defined;
}
/* -------------------------------------------------------------------------- */
inline constexpr auto Mesh::getAllFacetTypes(ElementType type) {
#define GET_FACET_TYPE(type) return ElementClass<type>::getFacetTypes();
AKANTU_BOOST_ALL_ELEMENT_SWITCH_NO_DEFAULT(GET_FACET_TYPE);
#undef GET_FACET_TYPE
return ElementClass<_not_defined>::getFacetTypes();
}
/* -------------------------------------------------------------------------- */
inline UInt Mesh::getNbFacetsPerElement(ElementType type) {
AKANTU_DEBUG_IN();
UInt n_facet = 0;
#define GET_NB_FACET(type) n_facet = ElementClass<type>::getNbFacetsPerElement()
AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_NB_FACET);
#undef GET_NB_FACET
AKANTU_DEBUG_OUT();
return n_facet;
}
/* -------------------------------------------------------------------------- */
inline UInt Mesh::getNbFacetsPerElement(ElementType type, UInt t) {
AKANTU_DEBUG_IN();
UInt n_facet = 0;
#define GET_NB_FACET(type) \
n_facet = ElementClass<type>::getNbFacetsPerElement(t)
AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_NB_FACET);
#undef GET_NB_FACET
AKANTU_DEBUG_OUT();
return n_facet;
}
/* -------------------------------------------------------------------------- */
inline auto Mesh::getFacetLocalConnectivity(ElementType type, UInt t) {
AKANTU_DEBUG_IN();
#define GET_FACET_CON(type) \
AKANTU_DEBUG_OUT(); \
return ElementClass<type>::getFacetLocalConnectivityPerElement(t)
AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_FACET_CON);
#undef GET_FACET_CON
AKANTU_DEBUG_OUT();
return ElementClass<_not_defined>::getFacetLocalConnectivityPerElement(0);
// This avoid a compilation warning but will certainly
// also cause a segfault if reached
}
/* -------------------------------------------------------------------------- */
inline auto Mesh::getFacetConnectivity(const Element & element, UInt t) const {
AKANTU_DEBUG_IN();
Matrix<const UInt> local_facets(getFacetLocalConnectivity(element.type, t));
Matrix<UInt> facets(local_facets.rows(), local_facets.cols());
const Array<UInt> & conn = connectivities(element.type, element.ghost_type);
for (UInt f = 0; f < facets.rows(); ++f) {
for (UInt n = 0; n < facets.cols(); ++n) {
facets(f, n) = conn(element.element, local_facets(f, n));
}
}
AKANTU_DEBUG_OUT();
return facets;
}
/* -------------------------------------------------------------------------- */
inline VectorProxy<UInt> Mesh::getConnectivity(const Element & element) const {
return connectivities.get(element);
}
/* -------------------------------------------------------------------------- */
inline VectorProxy<UInt> Mesh::getConnectivityNC(const Element & element) {
return connectivities.get(element);
}
/* -------------------------------------------------------------------------- */
template <typename T>
inline void Mesh::extractNodalValuesFromElement(
const Array<T> & nodal_values, T * local_coord, const UInt * connectivity,
UInt n_nodes, UInt nb_degree_of_freedom) const {
for (UInt n = 0; n < n_nodes; ++n) {
memcpy(local_coord + n * nb_degree_of_freedom,
nodal_values.storage() + connectivity[n] * nb_degree_of_freedom,
nb_degree_of_freedom * sizeof(T));
}
}
/* -------------------------------------------------------------------------- */
inline void Mesh::addConnectivityType(ElementType type, GhostType ghost_type) {
getConnectivityPointer(type, ghost_type);
}
/* -------------------------------------------------------------------------- */
inline bool Mesh::isPureGhostNode(UInt n) const {
return ((*nodes_flags)(n)&NodeFlag::_shared_mask) == NodeFlag::_pure_ghost;
}
/* -------------------------------------------------------------------------- */
inline bool Mesh::isLocalOrMasterNode(UInt n) const {
return ((*nodes_flags)(n)&NodeFlag::_local_master_mask) == NodeFlag::_normal;
}
/* -------------------------------------------------------------------------- */
inline bool Mesh::isLocalNode(UInt n) const {
return ((*nodes_flags)(n)&NodeFlag::_shared_mask) == NodeFlag::_normal;
}
/* -------------------------------------------------------------------------- */
inline bool Mesh::isMasterNode(UInt n) const {
return ((*nodes_flags)(n)&NodeFlag::_shared_mask) == NodeFlag::_master;
}
/* -------------------------------------------------------------------------- */
inline bool Mesh::isSlaveNode(UInt n) const {
return ((*nodes_flags)(n)&NodeFlag::_shared_mask) == NodeFlag::_slave;
}
/* -------------------------------------------------------------------------- */
inline bool Mesh::isPeriodicSlave(UInt n) const {
return ((*nodes_flags)(n)&NodeFlag::_periodic_mask) ==
NodeFlag::_periodic_slave;
}
/* -------------------------------------------------------------------------- */
inline bool Mesh::isPeriodicMaster(UInt n) const {
return ((*nodes_flags)(n)&NodeFlag::_periodic_mask) ==
NodeFlag::_periodic_master;
}
/* -------------------------------------------------------------------------- */
inline NodeFlag Mesh::getNodeFlag(UInt local_id) const {
return (*nodes_flags)(local_id);
}
/* -------------------------------------------------------------------------- */
inline Int Mesh::getNodePrank(UInt local_id) const {
auto it = nodes_prank.find(local_id);
return it == nodes_prank.end() ? -1 : it->second;
}
/* -------------------------------------------------------------------------- */
inline UInt Mesh::getNodeGlobalId(UInt local_id) const {
return nodes_global_ids ? (*nodes_global_ids)(local_id) : local_id;
}
/* -------------------------------------------------------------------------- */
inline UInt Mesh::getNodeLocalId(UInt global_id) const {
if (nodes_global_ids == nullptr) {
return global_id;
}
return nodes_global_ids->find(global_id);
}
/* -------------------------------------------------------------------------- */
inline UInt Mesh::getNbGlobalNodes() const {
return nodes_global_ids ? nb_global_nodes : nodes->size();
}
/* -------------------------------------------------------------------------- */
inline UInt Mesh::getNbNodesPerElementList(const Array<Element> & elements) {
UInt nb_nodes_per_element = 0;
UInt nb_nodes = 0;
ElementType current_element_type = _not_defined;
for (const auto & el : elements) {
if (el.type != current_element_type) {
current_element_type = el.type;
nb_nodes_per_element = Mesh::getNbNodesPerElement(current_element_type);
}
nb_nodes += nb_nodes_per_element;
}
return nb_nodes;
}
/* -------------------------------------------------------------------------- */
inline Mesh & Mesh::getMeshFacets() {
if (this->mesh_facets == nullptr) {
AKANTU_SILENT_EXCEPTION(
"No facet mesh is defined yet! check the buildFacets functions");
}
return *this->mesh_facets;
}
/* -------------------------------------------------------------------------- */
inline const Mesh & Mesh::getMeshFacets() const {
if (this->mesh_facets == nullptr) {
AKANTU_SILENT_EXCEPTION(
"No facet mesh is defined yet! check the buildFacets functions");
}
return *this->mesh_facets;
}
/* -------------------------------------------------------------------------- */
inline const Mesh & Mesh::getMeshParent() const {
if (this->mesh_parent == nullptr) {
AKANTU_SILENT_EXCEPTION(
"No parent mesh is defined! This is only valid in a mesh_facets");
}
return *this->mesh_parent;
}
/* -------------------------------------------------------------------------- */
void Mesh::addPeriodicSlave(UInt slave, UInt master) {
if (master == slave) {
return;
}
// if pair already registered
auto master_slaves = periodic_master_slave.equal_range(master);
auto slave_it =
std::find_if(master_slaves.first, master_slaves.second,
[&](auto & pair) { return pair.second == slave; });
if (slave_it == master_slaves.second) {
// no duplicates
periodic_master_slave.insert(std::make_pair(master, slave));
AKANTU_DEBUG_INFO("adding periodic slave, slave gid:"
<< getNodeGlobalId(slave) << " [lid: " << slave << "]"
<< ", master gid:" << getNodeGlobalId(master)
<< " [lid: " << master << "]");
// std::cout << "adding periodic slave, slave gid:" <<
// getNodeGlobalId(slave)
// << " [lid: " << slave << "]"
// << ", master gid:" << getNodeGlobalId(master)
// << " [lid: " << master << "]" << std::endl;
}
periodic_slave_master[slave] = master;
auto set_flag = [&](auto node, auto flag) {
(*nodes_flags)[node] &= ~NodeFlag::_periodic_mask; // clean periodic flags
(*nodes_flags)[node] |= flag;
};
set_flag(slave, NodeFlag::_periodic_slave);
set_flag(master, NodeFlag::_periodic_master);
}
/* -------------------------------------------------------------------------- */
UInt Mesh::getPeriodicMaster(UInt slave) const {
return periodic_slave_master.at(slave);
}
/* -------------------------------------------------------------------------- */
class Mesh::PeriodicSlaves {
using internal_iterator = std::unordered_multimap<UInt, UInt>::const_iterator;
std::pair<internal_iterator, internal_iterator> pair;
public:
PeriodicSlaves(const Mesh & mesh, UInt master)
: pair(mesh.getPeriodicMasterSlaves().equal_range(master)) {}
PeriodicSlaves(const PeriodicSlaves & other) = default;
PeriodicSlaves(PeriodicSlaves && other) = default;
PeriodicSlaves & operator=(const PeriodicSlaves & other) = default;
class const_iterator {
internal_iterator it;
public:
const_iterator(internal_iterator it) : it(it) {}
const_iterator operator++() {
++it;
return *this;
}
bool operator!=(const const_iterator & other) { return other.it != it; }
+ bool operator==(const const_iterator & other) { return other.it == it; }
auto operator*() { return it->second; }
};
auto begin() const { return const_iterator(pair.first); }
auto end() const { return const_iterator(pair.second); }
};
/* -------------------------------------------------------------------------- */
inline decltype(auto) Mesh::getPeriodicSlaves(UInt master) const {
return PeriodicSlaves(*this, master);
}
/* -------------------------------------------------------------------------- */
inline Vector<UInt>
Mesh::getConnectivityWithPeriodicity(const Element & element) const {
Vector<UInt> conn = getConnectivity(element);
if (not isPeriodic()) {
return conn;
}
for (auto && node : conn) {
if (isPeriodicSlave(node)) {
node = getPeriodicMaster(node);
}
}
return conn;
}
} // namespace akantu
#endif /* AKANTU_MESH_INLINE_IMPL_HH_ */