OGS
Mesh.cpp
Go to the documentation of this file.
1
15#include "Mesh.h"
16
17#include <memory>
18#include <range/v3/numeric.hpp>
19#include <range/v3/range/conversion.hpp>
20#include <range/v3/view/enumerate.hpp>
21#include <range/v3/view/indirect.hpp>
22#include <range/v3/view/map.hpp>
23#include <unordered_map>
24#include <utility>
25
26#include "BaseLib/RunTime.h"
27#include "Elements/Element.h"
28#include "Elements/Hex.h"
29#include "Elements/Prism.h"
30#include "Elements/Pyramid.h"
31#include "Elements/Quad.h"
32#include "Elements/Tet.h"
33#include "Elements/Tri.h"
34
35#ifdef USE_PETSC
37#endif
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 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{
77 this->resetNodeIDs();
78 this->resetElementIDs();
79 this->setDimension();
80
82
83 this->setElementNeighbors();
84}
85
86Mesh::Mesh(const Mesh& mesh)
87 : _id(global_mesh_counter++),
88 _mesh_dimension(mesh.getDimension()),
89 _node_distance(mesh._node_distance.first, mesh._node_distance.second),
90 _name(mesh.getName()),
91 _nodes(mesh.getNumberOfNodes()),
92 _elements(mesh.getNumberOfElements()),
93 _properties(mesh._properties)
94{
95 const std::vector<Node*>& nodes(mesh.getNodes());
96 const std::size_t nNodes(nodes.size());
97 for (unsigned i = 0; i < nNodes; ++i)
98 {
99 _nodes[i] = new Node(*nodes[i]);
100 }
101
102 const std::vector<Element*>& elements(mesh.getElements());
103 const std::size_t nElements(elements.size());
104 for (unsigned i = 0; i < nElements; ++i)
105 {
106 _elements[i] = elements[i]->clone();
107 for (auto const& [j, node_id] :
108 elements[i]->nodes() | views::ids | ranges::views::enumerate)
109 {
110 _elements[i]->setNode(static_cast<unsigned>(j), _nodes[node_id]);
111 }
112 }
113
114 if (_mesh_dimension == 0)
115 {
116 this->setDimension();
117 }
119 this->setElementNeighbors();
120}
121
122Mesh::Mesh(Mesh&& mesh) = default;
123
125{
126 _elements.clear();
127 _nodes.clear();
128}
129
131{
132 const std::size_t nElements(_elements.size());
133 for (std::size_t i = 0; i < nElements; ++i)
134 {
135 delete _elements[i];
136 }
137
138 const std::size_t nNodes(_nodes.size());
139 for (std::size_t i = 0; i < nNodes; ++i)
140 {
141 delete _nodes[i];
142 }
143}
144
146{
147 _elements.push_back(elem);
148}
149
151{
152 const std::size_t nNodes(_nodes.size());
153 for (std::size_t i = 0; i < nNodes; ++i)
154 {
155 _nodes[i]->setID(i);
156 }
157}
158
160{
161 const std::size_t nElements(this->_elements.size());
162 for (unsigned i = 0; i < nElements; ++i)
163 {
164 _elements[i]->setID(i);
165 }
166}
167
169{
170 const std::size_t nElements(_elements.size());
171 for (unsigned i = 0; i < nElements; ++i)
172 {
174 {
175 _mesh_dimension = _elements[i]->getDimension();
176 }
177 }
178}
179
180std::pair<double, double> minMaxEdgeLength(
181 std::vector<Element*> const& elements)
182{
183 auto min_max = [](auto const a, auto const b) -> std::pair<double, double> {
184 return {std::min(a.first, b.first), std::max(a.second, b.second)};
185 };
186
187 using limits = std::numeric_limits<double>;
188 auto const bounds = ranges::accumulate(
189 elements, std::pair{limits::infinity(), -limits::infinity()}, min_max,
190 [](Element* const e) { return computeSqrEdgeLengthRange(*e); });
191
192 return {std::sqrt(bounds.first), std::sqrt(bounds.second)};
193}
194
196{
197 std::vector<Element const*> neighbors;
198 for (auto element : _elements)
199 {
200 // create vector with all elements connected to current element
201 // (includes lots of doubles!)
202 const std::size_t nNodes(element->getNumberOfBaseNodes());
203 for (unsigned n(0); n < nNodes; ++n)
204 {
205 auto const& conn_elems(
206 _elements_connected_to_nodes[element->getNode(n)->getID()]);
207 neighbors.insert(neighbors.end(), conn_elems.begin(),
208 conn_elems.end());
209 }
210 std::sort(neighbors.begin(), neighbors.end());
211 auto const neighbors_new_end =
212 std::unique(neighbors.begin(), neighbors.end());
213
214 for (auto neighbor = neighbors.begin(); neighbor != neighbors_new_end;
215 ++neighbor)
216 {
217 std::optional<unsigned> const opposite_face_id =
218 element->addNeighbor(const_cast<Element*>(*neighbor));
219 if (opposite_face_id)
220 {
221 const_cast<Element*>(*neighbor)->setNeighbor(element,
222 *opposite_face_id);
223 }
224 }
225 neighbors.clear();
226 }
227}
228
230{
231 return std::count_if(begin(_nodes), end(_nodes),
232 [this](auto const* const node) {
233 return isBaseNode(
234 *node,
235 _elements_connected_to_nodes[node->getID()]);
236 });
237}
238
240{
241 return std::any_of(
242 std::begin(_elements), std::end(_elements),
243 [](Element const* const e)
244 { return e->getNumberOfNodes() != e->getNumberOfBaseNodes(); });
245}
246
247std::vector<MeshLib::Element const*> const& Mesh::getElementsConnectedToNode(
248 std::size_t const node_id) const
249{
250 return _elements_connected_to_nodes[node_id];
251}
252
253std::vector<MeshLib::Element const*> const& Mesh::getElementsConnectedToNode(
254 Node const& node) const
255{
256 return _elements_connected_to_nodes[node.getID()];
257}
258
260 std::string const& property_name,
261 double factor)
262{
263 if (!mesh.getProperties().existsPropertyVector<double>(property_name))
264 {
265 WARN("Did not find PropertyVector '{:s}' for scaling.", property_name);
266 return;
267 }
268 auto& pv = *mesh.getProperties().getPropertyVector<double>(property_name);
269 std::transform(pv.begin(), pv.end(), pv.begin(),
270 [factor](auto const& v) { return v * factor; });
271}
272
274{
275 auto const& properties = mesh.getProperties();
276 if (properties.existsPropertyVector<int>("MaterialIDs",
278 {
279 return properties.getPropertyVector<int>(
280 "MaterialIDs", MeshLib::MeshItemType::Cell, 1);
281 }
282 if (properties.hasPropertyVector("MaterialIDs"))
283 {
284 WARN(
285 "The 'MaterialIDs' mesh property exists but is either of wrong "
286 "type (must be int), or it is not defined on element / cell data.");
287 }
288 return nullptr;
289}
290
291std::unique_ptr<MeshLib::Mesh> createMeshFromElementSelection(
292 std::string mesh_name, std::vector<MeshLib::Element*> const& elements)
293{
294 auto ids_vector = views::ids | to<std::vector>();
295
296 DBUG("Found {:d} elements in the mesh", elements.size());
297
298 // Store bulk element ids for each of the new elements.
299 auto bulk_element_ids = elements | ids_vector;
300
301 // original node ids to newly created nodes.
302 std::unordered_map<std::size_t, MeshLib::Node*> id_node_hash_map;
303 id_node_hash_map.reserve(
304 elements.size()); // There will be at least one node per element.
305
306 for (auto& e : elements)
307 {
308 // For each node find a cloned node in map or create if there is none.
309 unsigned const n_nodes = e->getNumberOfNodes();
310 for (unsigned i = 0; i < n_nodes; ++i)
311 {
312 const MeshLib::Node* n = e->getNode(i);
313 auto const it = id_node_hash_map.find(n->getID());
314 if (it == id_node_hash_map.end())
315 {
316 auto new_node_in_map = id_node_hash_map[n->getID()] =
317 new MeshLib::Node(*n);
318 e->setNode(i, new_node_in_map);
319 }
320 else
321 {
322 e->setNode(i, it->second);
323 }
324 }
325 }
326
327 std::map<std::size_t, MeshLib::Node*> nodes_map;
328 for (const auto& n : id_node_hash_map)
329 {
330 nodes_map[n.first] = n.second;
331 }
332
333 // Copy the unique nodes pointers.
334 auto element_nodes = nodes_map | ranges::views::values | to<std::vector>;
335
336 // Store bulk node ids for each of the new nodes.
337 auto bulk_node_ids = nodes_map | ranges::views::keys | to<std::vector>;
338
339 auto mesh = std::make_unique<MeshLib::Mesh>(
340 std::move(mesh_name), std::move(element_nodes), std::move(elements));
341 assert(mesh != nullptr);
342
344 MeshLib::MeshItemType::Cell, 1, bulk_element_ids);
346 MeshLib::MeshItemType::Node, 1, bulk_node_ids);
347
348#ifdef USE_PETSC
349 return std::make_unique<MeshLib::NodePartitionedMesh>(*mesh);
350#else
351 return mesh;
352#endif
353}
354
355std::vector<std::vector<Node*>> calculateNodesConnectedByElements(
356 Mesh const& mesh)
357{
358 auto const elements_connected_to_nodes = findElementsConnectedToNodes(mesh);
359
360 std::vector<std::vector<Node*>> nodes_connected_by_elements;
361 auto const& nodes = mesh.getNodes();
362 nodes_connected_by_elements.resize(nodes.size());
363 for (std::size_t i = 0; i < nodes.size(); ++i)
364 {
365 auto& adjacent_nodes = nodes_connected_by_elements[i];
366 auto const node_id = nodes[i]->getID();
367
368 // Get all elements, to which this node is connected.
369 auto const& connected_elements = elements_connected_to_nodes[node_id];
370
371 // And collect all elements' nodes.
372 for (Element const* const element : connected_elements)
373 {
374 Node* const* const single_elem_nodes = element->getNodes();
375 std::size_t const nnodes = element->getNumberOfNodes();
376 for (std::size_t n = 0; n < nnodes; n++)
377 {
378 adjacent_nodes.push_back(single_elem_nodes[n]);
379 }
380 }
381
382 // Make nodes unique and sorted by their ids.
383 // This relies on the node's id being equivalent to it's address.
384 std::sort(adjacent_nodes.begin(), adjacent_nodes.end(),
385 [](Node* a, Node* b) { return a->getID() < b->getID(); });
386 auto const last =
387 std::unique(adjacent_nodes.begin(), adjacent_nodes.end());
388 adjacent_nodes.erase(last, adjacent_nodes.end());
389 }
390 return nodes_connected_by_elements;
391}
392
393bool isBaseNode(Node const& node,
394 std::vector<Element const*> const& elements_connected_to_node)
395{
396 // Check if node is connected.
397 if (elements_connected_to_node.empty())
398 {
399 return true;
400 }
401
402 // In a mesh a node always belongs to at least one element.
403 auto const e = elements_connected_to_node[0];
404
405 auto const n_base_nodes = e->getNumberOfBaseNodes();
406 auto const local_index = getNodeIDinElement(*e, &node);
407 return local_index < n_base_nodes;
408}
409} // namespace MeshLib
Definition of the Element class.
Definition of the Hex class.
void DBUG(fmt::format_string< Args... > fmt, Args &&... args)
Definition: Logging.h:30
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 mesh class for partitioned mesh (by node) for parallel computing within the framework o...
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
Definition: Point3dWithID.h:62
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:101
std::vector< std::vector< Element const * > > _elements_connected_to_nodes
Definition: Mesh.h:157
Mesh(std::string name, std::vector< Node * > nodes, std::vector< Element * > elements, Properties const &properties=Properties())
Definition: Mesh.cpp:63
void addElement(Element *elem)
Add an element to the mesh.
Definition: Mesh.cpp:145
std::vector< Element * > const & getElements() const
Get the element-vector for the mesh.
Definition: Mesh.h:104
unsigned _mesh_dimension
Definition: Mesh.h:149
std::size_t computeNumberOfBaseNodes() const
Get the number of base nodes.
Definition: Mesh.cpp:229
unsigned getDimension() const
Returns the dimension of the mesh (determined by the maximum dimension over all elements).
Definition: Mesh.h:83
void resetNodeIDs()
Resets the IDs of all mesh-nodes to their position in the node vector.
Definition: Mesh.cpp:150
Properties & getProperties()
Definition: Mesh.h:129
std::vector< Element * > _elements
Definition: Mesh.h:154
void setDimension()
Sets the dimension of the mesh.
Definition: Mesh.cpp:168
void resetElementIDs()
Resets the IDs of all mesh-elements to their position in the element vector.
Definition: Mesh.cpp:159
void setElementNeighbors()
Definition: Mesh.cpp:195
virtual ~Mesh()
Destructor.
Definition: Mesh.cpp:130
std::vector< Element const * > const & getElementsConnectedToNode(std::size_t node_id) const
Definition: Mesh.cpp:247
void shallowClean()
Definition: Mesh.cpp:124
bool hasNonlinearElement() const
Check if the mesh contains any nonlinear element.
Definition: Mesh.cpp:239
std::vector< Node * > _nodes
Definition: Mesh.h:153
Property manager on mesh items. Class Properties manages scalar, vector or matrix properties....
Definition: Properties.h:36
bool existsPropertyVector(std::string_view name) const
PropertyVector< T > const * getPropertyVector(std::string_view name) const
constexpr ranges::views::view_closure ids
For an element of a range view return its id.
Definition: Mesh.h:314
std::vector< std::vector< Element const * > > findElementsConnectedToNodes(Mesh const &mesh)
Definition: Mesh.cpp:46
std::unique_ptr< MeshLib::Mesh > createMeshFromElementSelection(std::string mesh_name, std::vector< MeshLib::Element * > const &elements)
Definition: Mesh.cpp:291
std::vector< std::vector< Node * > > calculateNodesConnectedByElements(Mesh const &mesh)
Definition: Mesh.cpp:355
PropertyVector< int > const * materialIDs(Mesh const &mesh)
Definition: Mesh.cpp:273
constexpr std::string_view getBulkIDString(MeshItemType mesh_item_type)
Definition: Properties.h:180
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:180
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:393
void scaleMeshPropertyVector(MeshLib::Mesh &mesh, std::string const &property_name, double factor)
Definition: Mesh.cpp:259
void addPropertyToMesh(Mesh &mesh, std::string_view name, MeshItemType item_type, std::size_t number_of_components, std::vector< T > const &values)
Definition: Mesh.h:195