preCICE v3.1.2
Loading...
Searching...
No Matches
Mesh.cpp
Go to the documentation of this file.
1#include <Eigen/Core>
2#include <Eigen/src/Core/Matrix.h>
3#include <algorithm>
4#include <array>
5#include <boost/container/flat_map.hpp>
6#include <functional>
7#include <memory>
8#include <ostream>
9#include <type_traits>
10#include <utility>
11#include <vector>
12
13#include "Edge.hpp"
14#include "Mesh.hpp"
15#include "Tetrahedron.hpp"
16#include "Triangle.hpp"
17#include "logging/LogMacros.hpp"
18#include "math/geometry.hpp"
19#include "mesh/Data.hpp"
20#include "mesh/Vertex.hpp"
22#include "query/Index.hpp"
23#include "utils/assertion.hpp"
24
25namespace precice::mesh {
26
29 int dimensions,
30 MeshID id)
31 : _name(std::move(name)),
32 _dimensions(dimensions),
33 _id(id),
34 _boundingBox(dimensions),
35 _index(*this)
36{
39}
40
42{
44 return _vertices.at(id);
45}
46
48{
50 return _vertices.at(id);
51}
52
57
59{
60 return _vertices;
61}
62
64{
65 return _vertices.size();
66}
67
72
74{
75 return _edges;
76}
77
82
84{
85 return _triangles;
86}
87
89{
90 return _tetrahedra;
91}
92
97
99{
100 return _dimensions;
101}
102
103Vertex &Mesh::createVertex(const Eigen::VectorXd &coords)
104{
105 PRECICE_ASSERT(coords.size() == _dimensions, coords.size(), _dimensions);
106 auto nextID = _vertices.size();
107 _vertices.emplace_back(coords, nextID);
108 return _vertices.back();
109}
110
112 Vertex &vertexOne,
113 Vertex &vertexTwo)
114{
115 _edges.emplace_back(vertexOne, vertexTwo);
116 return _edges.back();
117}
118
120 Edge &edgeOne,
121 Edge &edgeTwo,
122 Edge &edgeThree)
123{
125 edgeOne.connectedTo(edgeTwo) &&
126 edgeTwo.connectedTo(edgeThree) &&
127 edgeThree.connectedTo(edgeOne));
128 _triangles.emplace_back(edgeOne, edgeTwo, edgeThree);
129 return _triangles.back();
130}
131
133 Vertex &vertexOne,
134 Vertex &vertexTwo,
135 Vertex &vertexThree)
136{
137 _triangles.emplace_back(vertexOne, vertexTwo, vertexThree);
138 return _triangles.back();
139}
140
142 Vertex &vertexOne,
143 Vertex &vertexTwo,
144 Vertex &vertexThree,
145 Vertex &vertexFour)
146{
147 _tetrahedra.emplace_back(vertexOne, vertexTwo, vertexThree, vertexFour);
148 return _tetrahedra.back();
149}
150
152 const std::string &name,
153 int dimension,
154 DataID id,
155 int waveformDegree)
156{
157 PRECICE_TRACE(name, dimension);
158 for (const PtrData &data : _data) {
159 PRECICE_CHECK(data->getName() != name,
160 "Data \"{}\" cannot be created twice for mesh \"{}\". "
161 "Please rename or remove one of the use-data tags with name \"{}\".",
162 name, _name, name);
163 }
164 //#rows = dimensions of current mesh #columns = dimensions of corresponding data set
165 PtrData data(new Data(name, id, dimension, _dimensions, waveformDegree));
167 return _data.back();
168}
169
171{
172 return _data;
173}
174
175bool Mesh::hasDataID(DataID dataID) const
176{
177 auto iter = std::find_if(_data.begin(), _data.end(), [dataID](const auto &dptr) {
178 return dptr->getID() == dataID;
179 });
180 return iter != _data.end(); // if id was not found in mesh, iter == _data.end()
181}
182
183const PtrData &Mesh::data(DataID dataID) const
184{
185 auto iter = std::find_if(_data.begin(), _data.end(), [dataID](const auto &dptr) {
186 return dptr->getID() == dataID;
187 });
188 PRECICE_ASSERT(iter != _data.end(), "Data with id not found in mesh.", dataID, _name);
189 return *iter;
190}
191
193{
194 auto iter = std::find_if(_data.begin(), _data.end(), [&dataName](const auto &dptr) {
195 return dptr->getName() == dataName;
196 });
197 return iter != _data.end(); // if name was not found in mesh, iter == _data.end()
198}
199
201{
203 for (const auto &data : _data) {
204 names.push_back(data->getName());
205 }
206 return names;
207}
208
209const PtrData &Mesh::data(std::string_view dataName) const
210{
211 auto iter = std::find_if(_data.begin(), _data.end(), [&dataName](const auto &dptr) {
212 return dptr->getName() == dataName;
213 });
214 PRECICE_ASSERT(iter != _data.end(), "Data not found in mesh", dataName, _name);
215 return *iter;
216}
217
219{
220 return _name;
221}
222
224{
225 return _id;
226}
227
228bool Mesh::isValidVertexID(VertexID vertexID) const
229{
230 return (0 <= vertexID) && (static_cast<size_t>(vertexID) < nVertices());
231}
232
234{
236 const auto expectedCount = _vertices.size();
237 for (PtrData &data : _data) {
238 data->allocateValues(expectedCount);
239 }
240}
241
243{
245
246 // Keep the bounding box if set via the API function.
248
249 for (const Vertex &vertex : _vertices) {
250 bb.expandBy(vertex);
251 }
252 _boundingBox = std::move(bb);
253 PRECICE_DEBUG("Bounding Box, {}", _boundingBox);
254}
255
257{
259 _edges.clear();
262 _index.clear();
263
264 for (mesh::PtrData &data : _data) {
265 data->values().resize(0);
266 }
267}
268
278
279Eigen::VectorXd Mesh::getOwnedVertexData(const Eigen::VectorXd &values)
280{
281 std::vector<double> ownedDataVector;
282 PRECICE_ASSERT(static_cast<std::size_t>(values.size()) >= nVertices());
283 if (empty()) {
284 return {};
285 }
286 int valueDim = values.size() / nVertices();
287 int index = 0;
288
289 for (const auto &vertex : vertices()) {
290 if (vertex.isOwner()) {
291 for (int dim = 0; dim < valueDim; ++dim) {
292 ownedDataVector.push_back(values[index * valueDim + dim]);
293 }
294 }
295 ++index;
296 }
297 Eigen::Map<Eigen::VectorXd> ownedData(ownedDataVector.data(), ownedDataVector.size());
298
299 return ownedData;
300}
301
303{
304 for (auto &vertex : _vertices) {
305 vertex.tag();
306 }
307}
308
310 Mesh &deltaMesh)
311{
314
315 boost::container::flat_map<VertexID, Vertex *> vertexMap;
316 vertexMap.reserve(deltaMesh.nVertices());
317 Eigen::VectorXd coords(_dimensions);
318 for (const Vertex &vertex : deltaMesh.vertices()) {
319 coords = vertex.getCoords();
320 Vertex &v = createVertex(coords);
322 if (vertex.isTagged())
323 v.tag();
326 vertexMap[vertex.getID()] = &v;
327 }
328
329 // you cannot just take the vertices from the edge and add them,
330 // since you need the vertices from the new mesh
331 // (which may differ in IDs)
332 for (const Edge &edge : deltaMesh.edges()) {
333 VertexID vertexIndex1 = edge.vertex(0).getID();
334 VertexID vertexIndex2 = edge.vertex(1).getID();
335 PRECICE_ASSERT((vertexMap.count(vertexIndex1) == 1) &&
336 (vertexMap.count(vertexIndex2) == 1));
337 createEdge(*vertexMap[vertexIndex1], *vertexMap[vertexIndex2]);
338 }
339
340 for (const Triangle &triangle : deltaMesh.triangles()) {
341 VertexID vertexIndex1 = triangle.vertex(0).getID();
342 VertexID vertexIndex2 = triangle.vertex(1).getID();
343 VertexID vertexIndex3 = triangle.vertex(2).getID();
344 PRECICE_ASSERT((vertexMap.count(vertexIndex1) == 1) &&
345 (vertexMap.count(vertexIndex2) == 1) &&
346 (vertexMap.count(vertexIndex3) == 1));
347 createTriangle(*vertexMap[vertexIndex1], *vertexMap[vertexIndex2], *vertexMap[vertexIndex3]);
348 }
349
350 for (const Tetrahedron &tetra : deltaMesh.tetrahedra()) {
351 VertexID vertexIndex1 = tetra.vertex(0).getID();
352 VertexID vertexIndex2 = tetra.vertex(1).getID();
353 VertexID vertexIndex3 = tetra.vertex(2).getID();
354 VertexID vertexIndex4 = tetra.vertex(3).getID();
355
356 PRECICE_ASSERT((vertexMap.count(vertexIndex1) == 1) &&
357 (vertexMap.count(vertexIndex2) == 1) &&
358 (vertexMap.count(vertexIndex3) == 1) &&
359 (vertexMap.count(vertexIndex4) == 1));
360 createTetrahedron(*vertexMap[vertexIndex1], *vertexMap[vertexIndex2], *vertexMap[vertexIndex3], *vertexMap[vertexIndex4]);
361 }
362 _index.clear();
363}
364
366{
367 return _boundingBox;
368}
369
370void Mesh::expandBoundingBox(const BoundingBox &boundingBox)
371{
372 _boundingBox.expandBy(boundingBox);
373}
374
380
382{
383 // Remove duplicate tetrahedra
384 auto tetrahedraCnt = _tetrahedra.size();
386 auto lastTetrahedron = std::unique(_tetrahedra.begin(), _tetrahedra.end());
387 _tetrahedra = TetraContainer{_tetrahedra.begin(), lastTetrahedron};
388
389 // Remove duplicate triangles
390 auto triangleCnt = _triangles.size();
392 auto lastTriangle = std::unique(_triangles.begin(), _triangles.end());
393 _triangles = TriangleContainer{_triangles.begin(), lastTriangle};
394
395 // Remove duplicate edges
396 auto edgeCnt = _edges.size();
398 auto lastEdge = std::unique(_edges.begin(), _edges.end());
399 _edges = EdgeContainer{_edges.begin(), lastEdge};
400
401 PRECICE_DEBUG("Compression removed {} tetrahedra ({} to {}), {} triangles ({} to {}), and {} edges ({} to {})",
402 tetrahedraCnt - _tetrahedra.size(), tetrahedraCnt, _tetrahedra.size(),
403 triangleCnt - _triangles.size(), triangleCnt, _triangles.size(),
404 edgeCnt - _edges.size(), edgeCnt, _edges.size());
405}
406
407namespace {
408
409template <class Primitive, int... Indices>
410auto sortedVertexPtrsForImpl(Primitive &p, std::integer_sequence<int, Indices...>)
411{
412 std::array<Vertex *, Primitive::vertexCount> vs{&p.vertex(Indices)...};
413 std::sort(vs.begin(), vs.end());
414 return std::tuple_cat(vs);
415}
416
425template <class Primitive>
426auto sortedVertexPtrsFor(Primitive &p)
427{
428 return sortedVertexPtrsForImpl(p, std::make_integer_sequence<int, Primitive::vertexCount>{});
429}
430} // namespace
431
433{
434 if (_triangles.empty() && _tetrahedra.empty()) {
435 PRECICE_DEBUG("No implicit primitives required");
436 return; // no implicit primitives
437 }
438
439 // count explicit primitives for debug
440 const auto explTriangles = _triangles.size();
441 const auto explEdges = _triangles.size();
442
443 // First handle all explicit tetrahedra
444
445 // Build a set of all explicit triangles
446 using ExisitingTriangle = std::tuple<Vertex *, Vertex *, Vertex *>;
448 for (auto &t : _triangles) {
449 triangles.insert(sortedVertexPtrsFor(t));
450 }
451
452 // Generate all missing implicit triangles of explicit tetrahedra
453 // Update the triangles set used by the implicit edge generation
454 auto createTriangleIfMissing = [&](Vertex *a, Vertex *b, Vertex *c) {
455 if (triangles.count({a, b, c}) == 0) {
456 triangles.emplace(a, b, c);
457 createTriangle(*a, *b, *c);
458 };
459 };
460 for (auto &t : _tetrahedra) {
461 auto [a, b, c, d] = sortedVertexPtrsFor(t);
462 // Make sure these are in the same order as above
463 createTriangleIfMissing(a, b, c);
464 createTriangleIfMissing(a, b, d);
465 createTriangleIfMissing(a, c, d);
466 createTriangleIfMissing(b, c, d);
467 }
468
469 // Second handle all triangles, both explicit and implicit from the tetrahedron phase
470 // Build an set of all explicit triangles
471 using ExisitingEdge = std::tuple<Vertex *, Vertex *>;
473 for (auto &e : _edges) {
474 edges.emplace(sortedVertexPtrsFor(e));
475 }
476
477 // generate all missing implicit edges of implicit and explicit triangles
478 auto createEdgeIfMissing = [&](Vertex *a, Vertex *b) {
479 if (edges.count({a, b}) == 0) {
480 edges.emplace(a, b);
481 createEdge(*a, *b);
482 };
483 };
484 for (auto &t : _triangles) {
485 auto [a, b, c] = sortedVertexPtrsFor(t);
486 // Make sure these are in the same order as above
487 createEdgeIfMissing(a, b);
488 createEdgeIfMissing(a, c);
489 createEdgeIfMissing(b, c);
490 }
491
492 PRECICE_DEBUG("Generated {} implicit triangles and {} implicit edges",
493 _triangles.size() - explTriangles,
494 _edges.size() - explEdges);
495}
496
497bool Mesh::operator==(const Mesh &other) const
498{
499 bool equal = true;
500 equal &= _vertices.size() == other._vertices.size() &&
501 std::is_permutation(_vertices.begin(), _vertices.end(), other._vertices.begin());
502 equal &= _edges.size() == other._edges.size() &&
503 std::is_permutation(_edges.begin(), _edges.end(), other._edges.begin());
504 equal &= _triangles.size() == other._triangles.size() &&
505 std::is_permutation(_triangles.begin(), _triangles.end(), other._triangles.begin());
506 return equal;
507}
508
509bool Mesh::operator!=(const Mesh &other) const
510{
511 return !(*this == other);
512}
513
514std::ostream &operator<<(std::ostream &os, const Mesh &m)
515{
516 os << "Mesh \"" << m.getName() << "\", dimensionality = " << m.getDimensions() << ":\n";
517 os << "GEOMETRYCOLLECTION(\n";
518 const auto token = ", ";
519 const auto *sep = "";
520 for (auto &vertex : m.vertices()) {
521 os << sep << vertex;
522 sep = token;
523 }
524 sep = ",\n";
525 for (auto &edge : m.edges()) {
526 os << sep << edge;
527 sep = token;
528 }
529 sep = ",\n";
530 for (auto &triangle : m.triangles()) {
531 os << sep << triangle;
532 sep = token;
533 }
534 os << "\n)";
535 return os;
536}
537
538} // namespace precice::mesh
#define PRECICE_DEBUG(...)
Definition LogMacros.hpp:64
#define PRECICE_TRACE(...)
Definition LogMacros.hpp:95
#define PRECICE_CHECK(check,...)
Definition LogMacros.hpp:35
std::string name
#define PRECICE_ASSERT(...)
Definition assertion.hpp:87
T at(T... args)
T back(T... args)
T begin(T... args)
An axis-aligned bounding box around a (partition of a) mesh.
void expandBy(const BoundingBox &otherBB)
Expand bounding box using another bounding box.
Describes a set of data values belonging to the vertices of a mesh.
Definition Data.hpp:29
Linear edge of a mesh, defined by two Vertex objects.
Definition Edge.hpp:16
bool connectedTo(const Edge &other) const
Checks whether both edges share a vertex.
Definition Edge.cpp:39
Container and creator for meshes.
Definition Mesh.hpp:39
void expandBoundingBox(const BoundingBox &bounding_box)
Definition Mesh.cpp:370
Triangle & createTriangle(Edge &edgeOne, Edge &edgeTwo, Edge &edgeThree)
Creates and initializes a Triangle object.
Definition Mesh.cpp:119
MeshID _id
The ID of this mesh.
Definition Mesh.hpp:346
MeshID getID() const
Returns the base ID of the mesh.
Definition Mesh.cpp:223
std::string _name
Name of the mesh.
Definition Mesh.hpp:340
BoundingBox _boundingBox
Definition Mesh.hpp:390
int _globalNumberOfVertices
Number of unique vertices for complete distributed mesh.
Definition Mesh.hpp:376
int getDimensions() const
Definition Mesh.cpp:98
VertexContainer & vertices()
Returns modifieable container holding all vertices.
Definition Mesh.cpp:53
bool hasDataID(DataID dataID) const
Returns whether Mesh has Data with the matchingID.
Definition Mesh.cpp:175
Eigen::VectorXd getOwnedVertexData(const Eigen::VectorXd &values)
Definition Mesh.cpp:279
void clear()
Removes all mesh elements and data values (does not remove data or the bounding boxes).
Definition Mesh.cpp:256
DataContainer _data
Data hold by the vertices of the mesh.
Definition Mesh.hpp:355
void addMesh(Mesh &deltaMesh)
Definition Mesh.cpp:309
std::vector< std::string > availableData() const
Returns the names of all available data.
Definition Mesh.cpp:200
const std::string & getName() const
Returns the name of the mesh, as set in the config file.
Definition Mesh.cpp:218
void removeDuplicates()
Removes all duplicate connectivity.
Definition Mesh.cpp:381
std::size_t nVertices() const
Returns the number of vertices.
Definition Mesh.cpp:63
VertexDistribution _vertexDistribution
Vertex distribution for the primary rank, holding for each secondary rank all vertex IDs it owns.
Definition Mesh.hpp:362
TetraContainer & tetrahedra()
Returns modifiable container holding all tetrahedra.
Definition Mesh.cpp:93
std::vector< Rank > _connectedRanks
each rank stores list of connected remote ranks. In the m2n package, this is used to create the initi...
Definition Mesh.hpp:382
TetraContainer _tetrahedra
Definition Mesh.hpp:352
Vertex & vertex(VertexID id)
Mutable access to a vertex by VertexID.
Definition Mesh.cpp:41
Mesh(std::string name, int dimensions, MeshID id)
Constructor.
Definition Mesh.cpp:27
query::Index _index
Definition Mesh.hpp:392
VertexContainer _vertices
Holds vertices, edges, triangles and tetrahedra.
Definition Mesh.hpp:349
bool hasDataName(std::string_view dataName) const
Returns whether Mesh has Data with the dataName.
Definition Mesh.cpp:192
PtrData & createData(const std::string &name, int dimension, DataID id, int waveformDegree=time::Time::DEFAULT_WAVEFORM_DEGREE)
Create only data for vertex.
Definition Mesh.cpp:151
CommunicationMap _communicationMap
each rank stores list of connected ranks and corresponding vertex IDs here. In the m2n package,...
Definition Mesh.hpp:388
bool empty() const
Does the mesh contain any vertices?
Definition Mesh.hpp:88
void generateImplictPrimitives()
Definition Mesh.cpp:432
void clearPartitioning()
Clears the partitioning information.
Definition Mesh.cpp:270
void computeBoundingBox()
Computes the boundingBox for the vertices.
Definition Mesh.cpp:242
TriangleContainer _triangles
Definition Mesh.hpp:351
Tetrahedron & createTetrahedron(Vertex &vertexOne, Vertex &vertexTwo, Vertex &vertexThree, Vertex &vertexFour)
Creates and initializes a Tetrahedron object.
Definition Mesh.cpp:141
VertexOffsets _vertexOffsets
Holds the index of the last vertex for each rank.
Definition Mesh.hpp:369
bool isValidVertexID(VertexID vertexID) const
Returns true if the given vertexID is valid.
Definition Mesh.cpp:228
EdgeContainer _edges
Definition Mesh.hpp:350
const DataContainer & data() const
Allows access to all data.
Definition Mesh.cpp:170
TriangleContainer & triangles()
Returns modifiable container holding all triangles.
Definition Mesh.cpp:78
const BoundingBox & getBoundingBox() const
Returns the bounding box of the mesh.
Definition Mesh.cpp:365
query::Index & index()
Call preprocess() before index() to ensure correct projection handling.
Definition Mesh.hpp:321
Edge & createEdge(Vertex &vertexOne, Vertex &vertexTwo)
Creates and initializes an Edge object.
Definition Mesh.cpp:111
EdgeContainer & edges()
Returns modifiable container holding all edges.
Definition Mesh.cpp:68
Vertex & createVertex(const Eigen::VectorXd &coords)
Creates and initializes a Vertex object.
Definition Mesh.cpp:103
void allocateDataValues()
Allocates memory for the vertex data values and corresponding gradient values.
Definition Mesh.cpp:233
int _dimensions
Dimension of mesh.
Definition Mesh.hpp:343
Tetrahedron of a mesh, defined by 4 vertices.
Triangle of a mesh, defined by three vertices.
Definition Triangle.hpp:27
Vertex of a mesh.
Definition Vertex.hpp:16
VertexID getID() const
Returns the unique (among vertices of one mesh on one processor) ID of the vertex.
Definition Vertex.hpp:111
bool isTagged() const
Definition Vertex.cpp:32
void setGlobalIndex(int globalIndex)
Definition Vertex.cpp:17
void setOwner(bool owner)
Definition Vertex.cpp:27
bool isOwner() const
Definition Vertex.cpp:22
Eigen::VectorXd getCoords() const
Returns the coordinates of the vertex.
Definition Vertex.hpp:116
int getGlobalIndex() const
Globally unique index.
Definition Vertex.cpp:12
void clear()
Clear the index.
Definition Index.cpp:391
T clear(T... args)
T count(T... args)
T data(T... args)
T emplace_back(T... args)
T emplace(T... args)
T empty(T... args)
T end(T... args)
T find_if(T... args)
T insert(T... args)
T is_permutation(T... args)
provides Mesh, Data and primitives.
int MeshID
Definition Types.hpp:30
int VertexID
Definition Types.hpp:13
int DataID
Definition Types.hpp:25
STL namespace.
T push_back(T... args)
T resize(T... args)
T size(T... args)
T sort(T... args)
T tuple_cat(T... args)
T unique(T... args)