OGS
Mesh.cpp
Go to the documentation of this file.
1
15#include "Mesh.h"
16
17#include <memory>
18#include <range/v3/algorithm/contains.hpp>
19#include <range/v3/numeric.hpp>
20#include <range/v3/range/conversion.hpp>
21#include <range/v3/view/enumerate.hpp>
22#include <range/v3/view/indirect.hpp>
23#include <range/v3/view/map.hpp>
24#include <unordered_map>
25#include <utility>
26
27#include "BaseLib/RunTime.h"
28#include "Elements/Element.h"
29#include "Elements/Hex.h"
30#include "Elements/Prism.h"
31#include "Elements/Pyramid.h"
32#include "Elements/Quad.h"
33#include "Elements/Tet.h"
34#include "Elements/Tri.h"
37
39static std::size_t global_mesh_counter = 0;
40
41namespace MeshLib
42{
43using namespace ranges;
44
45std::vector<std::vector<Element const*>> findElementsConnectedToNodes(
46 Mesh const& mesh)
47{
48 std::vector<std::vector<Element const*>> elements_connected_to_nodes;
49 auto const& nodes = mesh.getNodes();
50 elements_connected_to_nodes.resize(nodes.size());
51
52 for (auto const* element : mesh.getElements())
53 {
54 for (auto const node_id : element->nodes() | views::ids)
55 {
56 elements_connected_to_nodes[node_id].push_back(element);
57 }
58 }
59 return elements_connected_to_nodes;
60}
61
62Mesh::Mesh(std::string name,
63 std::vector<Node*>
64 nodes,
65 std::vector<Element*>
66 elements,
67 bool const compute_element_neighbors,
68 Properties const& properties)
69 : _id(global_mesh_counter++),
70 _mesh_dimension(0),
71 _node_distance(std::numeric_limits<double>::max(), 0),
72 _name(std::move(name)),
73 _nodes(std::move(nodes)),
74 _elements(std::move(elements)),
75 _properties(properties),
76 _compute_element_neighbors(compute_element_neighbors)
77{
78 this->resetNodeIDs();
79 this->resetElementIDs();
80 this->setDimension();
81
83
85 {
86 this->setElementNeighbors();
87 }
88}
89
90Mesh::Mesh(const Mesh& mesh)
91 : _id(global_mesh_counter++),
92 _mesh_dimension(mesh.getDimension()),
93 _node_distance(mesh._node_distance.first, mesh._node_distance.second),
94 _name(mesh.getName()),
95 _nodes(mesh.getNumberOfNodes()),
96 _elements(mesh.getNumberOfElements()),
97 _properties(mesh._properties),
98 _compute_element_neighbors(mesh._compute_element_neighbors)
99{
100 const std::vector<Node*>& nodes(mesh.getNodes());
101 const std::size_t nNodes(nodes.size());
102 for (unsigned i = 0; i < nNodes; ++i)
103 {
104 _nodes[i] = new Node(*nodes[i]);
105 }
106
107 const std::vector<Element*>& elements(mesh.getElements());
108 const std::size_t nElements(elements.size());
109 for (unsigned i = 0; i < nElements; ++i)
110 {
111 _elements[i] = elements[i]->clone();
112 for (auto const& [j, node_id] :
113 elements[i]->nodes() | views::ids | ranges::views::enumerate)
114 {
115 _elements[i]->setNode(static_cast<unsigned>(j), _nodes[node_id]);
116 }
117 }
118
119 if (_mesh_dimension == 0)
120 {
121 this->setDimension();
122 }
123
126 {
127 this->setElementNeighbors();
128 }
129}
130
131Mesh::Mesh(Mesh&& mesh) = default;
132
134{
135 _elements.clear();
136 _nodes.clear();
137}
138
140{
141 const std::size_t nElements(_elements.size());
142 for (std::size_t i = 0; i < nElements; ++i)
143 {
144 delete _elements[i];
145 }
146
147 const std::size_t nNodes(_nodes.size());
148 for (std::size_t i = 0; i < nNodes; ++i)
149 {
150 delete _nodes[i];
151 }
152}
153
155{
156 _elements.push_back(elem);
157}
158
160{
161 const std::size_t nNodes(_nodes.size());
162 for (std::size_t i = 0; i < nNodes; ++i)
163 {
164 _nodes[i]->setID(i);
165 }
166}
167
169{
170 const std::size_t nElements(this->_elements.size());
171 for (unsigned i = 0; i < nElements; ++i)
172 {
173 _elements[i]->setID(i);
174 }
175}
176
178{
179 const std::size_t nElements(_elements.size());
180 for (unsigned i = 0; i < nElements; ++i)
181 {
183 {
184 _mesh_dimension = _elements[i]->getDimension();
185 }
186 }
187}
188
189std::pair<double, double> minMaxEdgeLength(
190 std::vector<Element*> const& elements)
191{
192 auto min_max = [](auto const a, auto const b) -> std::pair<double, double> {
193 return {std::min(a.first, b.first), std::max(a.second, b.second)};
194 };
195
196 using limits = std::numeric_limits<double>;
197 auto const bounds = ranges::accumulate(
198 elements, std::pair{limits::infinity(), -limits::infinity()}, min_max,
199 [](Element const* const e) { return computeSqrEdgeLengthRange(*e); });
200
201 return {std::sqrt(bounds.first), std::sqrt(bounds.second)};
202}
203
205{
206 std::vector<Element const*> neighbors;
207 for (auto element : _elements)
208 {
209 // create vector with all elements connected to current element
210 // (includes lots of doubles!)
211 const std::size_t nNodes(element->getNumberOfBaseNodes());
212 for (unsigned n(0); n < nNodes; ++n)
213 {
214 auto const& conn_elems(
215 _elements_connected_to_nodes[element->getNode(n)->getID()]);
216 neighbors.insert(neighbors.end(), conn_elems.begin(),
217 conn_elems.end());
218 }
219 std::sort(neighbors.begin(), neighbors.end());
220 auto const neighbors_new_end =
221 std::unique(neighbors.begin(), neighbors.end());
222
223 for (auto neighbor = neighbors.begin(); neighbor != neighbors_new_end;
224 ++neighbor)
225 {
226 std::optional<unsigned> const opposite_face_id =
227 element->addNeighbor(const_cast<Element*>(*neighbor));
228 if (opposite_face_id)
229 {
230 const_cast<Element*>(*neighbor)->setNeighbor(element,
231 *opposite_face_id);
232 }
233 }
234 neighbors.clear();
235 }
236}
237
239{
240 return std::count_if(begin(_nodes), end(_nodes),
241 [this](auto const* const node) {
242 return isBaseNode(
243 *node,
244 _elements_connected_to_nodes[node->getID()]);
245 });
246}
247
249{
250 return std::any_of(
251 std::begin(_elements), std::end(_elements),
252 [](Element const* const e)
253 { return e->getNumberOfNodes() != e->getNumberOfBaseNodes(); });
254}
255
256std::vector<MeshLib::Element const*> const& Mesh::getElementsConnectedToNode(
257 std::size_t const node_id) const
258{
259 return _elements_connected_to_nodes[node_id];
260}
261
262std::vector<MeshLib::Element const*> const& Mesh::getElementsConnectedToNode(
263 Node const& node) const
264{
265 return _elements_connected_to_nodes[node.getID()];
266}
267
269{
270 auto const& properties = mesh.getProperties();
271 if (properties.existsPropertyVector<int>("MaterialIDs",
273 {
274 return properties.getPropertyVector<int>(
275 "MaterialIDs", MeshLib::MeshItemType::Cell, 1);
276 }
277 if (properties.hasPropertyVector("MaterialIDs"))
278 {
279 WARN(
280 "The 'MaterialIDs' mesh property exists but is either of wrong "
281 "type (must be int), or it is not defined on element / cell data.");
282 }
283 return nullptr;
284}
285
287{
288 return const_cast<PropertyVector<int>*>(
289 MeshLib::materialIDs(std::as_const(mesh)));
290}
291
293{
294 auto const& properties = mesh.getProperties();
295 return properties.getPropertyVector<std::size_t>(
298}
299
301{
302 auto const& properties = mesh.getProperties();
303 return properties.getPropertyVector<std::size_t>(
306}
307
308std::vector<std::vector<Node*>> calculateNodesConnectedByElements(
309 Mesh const& mesh)
310{
311 auto const elements_connected_to_nodes = findElementsConnectedToNodes(mesh);
312
313 std::vector<std::vector<Node*>> nodes_connected_by_elements;
314 auto const& nodes = mesh.getNodes();
315 nodes_connected_by_elements.resize(nodes.size());
316 for (std::size_t i = 0; i < nodes.size(); ++i)
317 {
318 auto& adjacent_nodes = nodes_connected_by_elements[i];
319 auto const node_id = nodes[i]->getID();
320
321 // Get all elements, to which this node is connected.
322 auto const& connected_elements = elements_connected_to_nodes[node_id];
323
324 // And collect all elements' nodes.
325 for (Element const* const element : connected_elements)
326 {
327 Node* const* const single_elem_nodes = element->getNodes();
328 std::size_t const nnodes = element->getNumberOfNodes();
329 for (std::size_t n = 0; n < nnodes; n++)
330 {
331 adjacent_nodes.push_back(single_elem_nodes[n]);
332 }
333 }
334
335 // Make nodes unique and sorted by their ids.
336 // This relies on the node's id being equivalent to it's address.
337 std::sort(adjacent_nodes.begin(), adjacent_nodes.end(),
339 auto const last =
340 std::unique(adjacent_nodes.begin(), adjacent_nodes.end());
341 adjacent_nodes.erase(last, adjacent_nodes.end());
342 }
343 return nodes_connected_by_elements;
344}
345
346bool isBaseNode(Node const& node,
347 std::vector<Element const*> const& elements_connected_to_node)
348{
349 // Check if node is connected.
350 if (elements_connected_to_node.empty())
351 {
352 return true;
353 }
354
355 // In a mesh a node always belongs to at least one element.
356 auto const e = elements_connected_to_node[0];
357
358 auto const n_base_nodes = e->getNumberOfBaseNodes();
359 auto const local_index = getNodeIDinElement(*e, &node);
360 return local_index < n_base_nodes;
361}
362
363Mesh& findMeshByName(std::vector<std::unique_ptr<Mesh>> const& meshes,
364 std::string_view const name)
365{
367 meshes,
368 [&name](auto const& mesh)
369 {
370 assert(mesh != nullptr);
371 return mesh->getName() == name;
372 },
373 [&]() { OGS_FATAL("Required mesh named {:s} not found.", name); });
374}
375
376} // namespace MeshLib
Definition of Duplicate functions.
Definition of the Element class.
#define OGS_FATAL(...)
Definition Error.h:26
Definition of the Hex class.
void WARN(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:40
static std::size_t global_mesh_counter
Mesh counter used to uniquely identify meshes by id.
Definition Mesh.cpp:39
Definition of the Mesh class.
Definition of the Prism class.
Definition of the Pyramid class.
Definition of the Quad class.
Definition of the RunTime class.
std::string getName(std::string const &line)
Returns the name/title from the "Zone"-description.
Definition of the Tet class.
Definition of the Tri class.
std::size_t getID() const
void setNeighbor(Element *neighbor, unsigned const face_id)
Definition Element.cpp:32
virtual unsigned getNumberOfNodes() const =0
virtual unsigned getNumberOfBaseNodes() const =0
std::vector< Node * > const & getNodes() const
Get the nodes-vector for the mesh.
Definition Mesh.h:106
bool const _compute_element_neighbors
Definition Mesh.h:165
std::vector< std::vector< Element const * > > _elements_connected_to_nodes
Definition Mesh.h:162
void addElement(Element *elem)
Add an element to the mesh.
Definition Mesh.cpp:154
std::vector< Element * > const & getElements() const
Get the element-vector for the mesh.
Definition Mesh.h:109
unsigned _mesh_dimension
Definition Mesh.h:154
std::size_t computeNumberOfBaseNodes() const
Get the number of base nodes.
Definition Mesh.cpp:238
unsigned getDimension() const
Returns the dimension of the mesh (determined by the maximum dimension over all elements).
Definition Mesh.h:88
void resetNodeIDs()
Resets the IDs of all mesh-nodes to their position in the node vector.
Definition Mesh.cpp:159
Properties & getProperties()
Definition Mesh.h:134
std::vector< Element * > _elements
Definition Mesh.h:159
void setDimension()
Sets the dimension of the mesh.
Definition Mesh.cpp:177
void resetElementIDs()
Resets the IDs of all mesh-elements to their position in the element vector.
Definition Mesh.cpp:168
void setElementNeighbors()
Definition Mesh.cpp:204
virtual ~Mesh()
Destructor.
Definition Mesh.cpp:139
Mesh(std::string name, std::vector< Node * > nodes, std::vector< Element * > elements, bool const compute_element_neighbors=false, Properties const &properties=Properties())
Definition Mesh.cpp:62
std::vector< Element const * > const & getElementsConnectedToNode(std::size_t node_id) const
Definition Mesh.cpp:256
void shallowClean()
Definition Mesh.cpp:133
bool hasNonlinearElement() const
Check if the mesh contains any nonlinear element.
Definition Mesh.cpp:248
std::vector< Node * > _nodes
Definition Mesh.h:158
Property manager on mesh items. Class Properties manages scalar, vector or matrix properties....
Definition Properties.h:36
PropertyVector< T > const * getPropertyVector(std::string_view name) const
ranges::range_reference_t< Range > findElementOrError(Range &range, std::predicate< ranges::range_reference_t< Range > > auto &&predicate, std::invocable auto error_callback)
Definition Algorithm.h:77
constexpr ranges::views::view_closure ids
For an element of a range view return its id.
Definition Mesh.h:225
std::vector< std::vector< Element const * > > findElementsConnectedToNodes(Mesh const &mesh)
Definition Mesh.cpp:45
std::vector< std::vector< Node * > > calculateNodesConnectedByElements(Mesh const &mesh)
Definition Mesh.cpp:308
Mesh & findMeshByName(std::vector< std::unique_ptr< Mesh > > const &meshes, std::string_view const name)
Definition Mesh.cpp:363
bool idsComparator(T const a, T const b)
Definition Mesh.h:206
PropertyVector< int > const * materialIDs(Mesh const &mesh)
Definition Mesh.cpp:268
constexpr std::string_view getBulkIDString(MeshItemType mesh_item_type)
Definition Properties.h:188
PropertyVector< std::size_t > const * bulkElementIDs(Mesh const &mesh)
Definition Mesh.cpp:300
std::pair< double, double > computeSqrEdgeLengthRange(Element const &element)
Compute the minimum and maximum squared edge length for this element.
Definition Element.cpp:156
std::pair< double, double > minMaxEdgeLength(std::vector< Element * > const &elements)
Returns the minimum and maximum edge length for given elements.
Definition Mesh.cpp:189
unsigned getNodeIDinElement(Element const &element, const MeshLib::Node *node)
Returns the position of the given node in the node array of this element.
Definition Element.cpp:206
bool isBaseNode(Node const &node, std::vector< Element const * > const &elements_connected_to_node)
Definition Mesh.cpp:346
PropertyVector< std::size_t > const * bulkNodeIDs(Mesh const &mesh)
Definition Mesh.cpp:292