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