preCICE v3.1.2
Loading...
Searching...
No Matches
CreateClustering.hpp
Go to the documentation of this file.
1#pragma once
2
3#include <Eigen/Core>
4
5#include <algorithm>
6
7#include "math/math.hpp"
8#include "mesh/Mesh.hpp"
10#include "query/Index.hpp"
11// required for the squared distance computation
13
14namespace precice {
15namespace mapping {
16namespace impl {
17
19
25// Internal helper functions
26namespace {
27
39{
40 return (ids[0] * nCells[1] + ids[1]) * nCells[2] + ids[2];
41}
42
43// Overload to handle the offsets (via ints) properly
45{
46 return (ids[0] * nCells[1] + ids[1]) * nCells[2] + ids[2];
47}
48
56template <int dim>
57std::array<int, math::pow_int<dim>(3) - 1> zCurveNeighborOffsets(std::array<unsigned int, 3> nCells);
58
60template <>
61std::array<int, 8> zCurveNeighborOffsets<2>(std::array<unsigned int, 3> nCells)
62{
63 // All neighbor directions
64 const std::array<std::array<int, 3>, 8> neighbors2DIndex = {{{-1, -1, 0}, {-1, 0, 0}, {-1, 1, 0}, {0, -1, 0}, {0, 1, 0}, {1, -1, 0}, {1, 0, 0}, {1, 1, 0}}};
65 std::array<int, 8> neighbors2D{{}};
66 // Uses the int overload of the zCurve
67 // @todo: mark the function as constexpr, when moving to C++20 (transform is non-constexpr)
68 std::transform(neighbors2DIndex.begin(), neighbors2DIndex.end(), neighbors2D.begin(), [&nCells](auto &v) { return zCurve(v, nCells); });
69
70 return neighbors2D;
71}
72
74template <>
75std::array<int, 26> zCurveNeighborOffsets<3>(std::array<unsigned int, 3> nCells)
76{
77 // All neighbor directions
78 const std::array<std::array<int, 3>, 26> neighbors3DIndex = {{{-1, -1, -1}, {-1, -1, 0}, {-1, -1, 1}, {-1, 0, 0}, {-1, 0, 1}, {-1, 0, -1}, {-1, 1, 0}, {-1, 1, 1}, {-1, 1, -1}, {0, -1, -1}, {0, -1, 0}, {0, -1, 1}, {0, 0, 1}, {0, 0, -1}, {0, 1, 0}, {0, 1, 1}, {0, 1, -1}, {1, -1, -1}, {1, -1, 0}, {1, -1, 1}, {1, 0, 0}, {1, 0, 1}, {1, 0, -1}, {1, 1, 0}, {1, 1, 1}, {1, 1, -1}}};
79 std::array<int, 26> neighbors3D{{}};
80 // Uses the int overload of the zCurve
81 std::transform(neighbors3DIndex.begin(), neighbors3DIndex.end(), neighbors3D.begin(), [&nCells](auto &v) { return zCurve(v, nCells); });
82
83 return neighbors3D;
84}
98template <typename ArrayType>
99std::vector<VertexID> getNeighborsWithinSphere(const Vertices &vertices, VertexID centerID, double radius, const ArrayType neighborOffsets)
100{
101 // number of neighbors = (3^dim)-1
102 static_assert((neighborOffsets.size() == 26) || (neighborOffsets.size() == 8));
104 for (auto off : neighborOffsets) {
105 // corner case where we have 'dead axis' in the bounding box
106 if (off == 0) {
107 continue;
108 }
109 // compute neighbor ID
110 // @todo target type should be an unsigned int
111 PRECICE_ASSERT(off != 0);
112 // Assumption, vertex position in the container == vertex ID
113 VertexID neighborID = centerID + off;
114 // We might run out of the array bounds along the edges, so we only consider vertices inside
115 if (neighborID >= 0 && neighborID < static_cast<int>(vertices.size())) {
116 auto &neighborVertex = vertices[neighborID];
117 auto &center = vertices[centerID];
118 if (!neighborVertex.isTagged() && (computeSquaredDifference(center.rawCoords(), neighborVertex.rawCoords()) < math::pow_int<2>(radius))) {
119 result.emplace_back(neighborID);
120 }
121 }
122 }
123 return result;
124}
125
134void tagEmptyClusters(Vertices &clusterCenters, double clusterRadius, mesh::PtrMesh mesh)
135{
136 // Alternative implementation: mesh->index().getVerticesInsideBox() == 0
137 std::for_each(clusterCenters.begin(), clusterCenters.end(), [&](auto &v) {
138 if (!v.isTagged() && !mesh->index().isAnyVertexInsideBox(v, clusterRadius)) {
139 v.tag();
140 }
141 });
142}
143
152void projectClusterCentersToinputMesh(Vertices &clusterCenters, mesh::PtrMesh mesh)
153{
154 std::transform(clusterCenters.begin(), clusterCenters.end(), clusterCenters.begin(), [&](auto &v) {
155 if (!v.isTagged()) {
156 auto closestCenter = mesh->index().getClosestVertex(v.getCoords()).index;
157 return mesh::Vertex{mesh->vertex(closestCenter).getCoords(), v.getID()};
158 } else {
159 return v;
160 }
161 });
162}
163
175template <int dim>
176void tagDuplicateCenters(Vertices &centers, std::array<unsigned int, 3> nCells, double threshold)
177{
178 PRECICE_ASSERT(threshold >= 0);
179 if (centers.empty())
180 return;
181
182 // Get the index offsets for the z curve
183 auto neighborOffsets = zCurveNeighborOffsets<dim>(nCells);
184 static_assert((neighborOffsets.size() == 8 && dim == 2) || (neighborOffsets.size() == 26 && dim == 3));
185
186 //For the following to work, the vertexID has to correspond to the position in the centers
187 PRECICE_ASSERT(std::all_of(centers.begin(), centers.end(), [idx = 0](auto &v) mutable { return v.getID() == idx++; }));
188 // we check all neighbors
189 for (auto &v : centers) {
190 if (!v.isTagged()) {
191 auto ids = getNeighborsWithinSphere(centers, v.getID(), threshold, neighborOffsets);
192 // @todo check more selective which one to remove
193 if (ids.size() > 0)
194 v.tag();
195 }
196 }
197}
198
207void removeTaggedVertices(Vertices &container)
208{
209 container.erase(std::remove_if(container.begin(), container.end(), [](auto &v) { return v.isTagged(); }), container.end());
210}
211} // namespace
212
225inline double estimateClusterRadius(unsigned int verticesPerCluster, mesh::PtrMesh inMesh, const mesh::BoundingBox &bb)
226{
227 // Step 1: Generate random samples from the input mesh
228 std::vector<VertexID> randomSamples;
229
230 PRECICE_ASSERT(!bb.empty(), "Center computation won't work right now.");
231 // All samples are 'moved' to the closest vertex from the input mesh in order
232 // to avoid empty clusters when we have shell-like meshes (as in surface coupling)
233 // The first sample: the vertex closest to the bounding box center
234 const auto bbCenter = bb.center();
235 randomSamples.emplace_back(inMesh->index().getClosestVertex(bbCenter).index);
236 // Now we place samples for each space dimension above and below the bounding box center
237 // and look for the closest vertex in the input mesh
238 // In 2D the sampling would like this
239 // +---------------------+
240 // | x |
241 // | |
242 // | x c x |
243 // | |
244 // | x |
245 // +---------------------+
246
247 for (int d = 0; d < inMesh->getDimensions(); ++d) {
248 auto sample = bbCenter;
249 // Lower as the center
250 sample[d] -= bb.getEdgeLength(d) * 0.25;
251 randomSamples.emplace_back(inMesh->index().getClosestVertex(sample).index);
252 // Higher than the center
253 sample[d] += bb.getEdgeLength(d) * 0.5;
254 randomSamples.emplace_back(inMesh->index().getClosestVertex(sample).index);
255 }
256
257 // Step 2: Compute the radius of the randomSamples ('centers'), which would have verticesPerCluster vertices
258 std::vector<double> sampledClusterRadii;
259 for (auto s : randomSamples) {
260 // ask the index tree for the k-nearest neighbors in order to estimate the point density
261 auto kNearestVertexIDs = inMesh->index().getClosestVertices(inMesh->vertex(s).getCoords(), verticesPerCluster);
262 // compute the distance of each point to the center
263 std::vector<double> squaredRadius(kNearestVertexIDs.size());
264 std::transform(kNearestVertexIDs.begin(), kNearestVertexIDs.end(), squaredRadius.begin(), [&inMesh, s](auto i) {
265 return computeSquaredDifference(inMesh->vertex(i).rawCoords(), inMesh->vertex(s).rawCoords());
266 });
267 // Store the maximum distance
268 auto maxRadius = std::max_element(squaredRadius.begin(), squaredRadius.end());
269 sampledClusterRadii.emplace_back(std::sqrt(*maxRadius));
270 }
271
272 // Step 3: Sort (using nth_element) the sampled radii and select the median as cluster radius
273 PRECICE_ASSERT(sampledClusterRadii.size() % 2 != 0, "Median calculation is only valid for odd number of elements.");
274 PRECICE_ASSERT(sampledClusterRadii.size() > 0);
275 unsigned int middle = sampledClusterRadii.size() / 2;
276 std::nth_element(sampledClusterRadii.begin(), sampledClusterRadii.begin() + middle, sampledClusterRadii.end());
277 double clusterRadius = sampledClusterRadii[middle];
278
279 return clusterRadius;
280}
281
305 double relativeOverlap, unsigned int verticesPerCluster,
306 bool projectClustersToInput)
307{
308 precice::logging::Logger _log{"impl::createClustering"};
310 PRECICE_ASSERT(relativeOverlap < 1);
311 PRECICE_ASSERT(verticesPerCluster > 0);
312 PRECICE_ASSERT(inMesh->getDimensions() == outMesh->getDimensions());
313
314 // If we have either no input or no output vertices, we return immediately
315 if (inMesh->nVertices() == 0 || outMesh->nVertices() == 0)
316 return {double{}, Vertices{}};
317
318 PRECICE_ASSERT(!outMesh->empty() && !inMesh->empty());
319
320 // startGridAtEdge boolean switch in order to decide either to start the clustering at the edge of the bounding box in each direction
321 // (true) or start the clustering inside the bounding box (edge + 0.5 radius). The latter approach leads to fewer clusters,
322 // but might result in a worse clustering
323 const bool startGridAtEdge = projectClustersToInput;
324
325 PRECICE_DEBUG("Relative overlap: {}", relativeOverlap);
326 PRECICE_DEBUG("Vertices per cluster: {}", verticesPerCluster);
327
328 // Step 1: Compute the local bounding box of the input mesh manually
329 // Note that we don't use the corresponding bounding box functions from
330 // precice::mesh (e.g. ::getBoundingBox), as the stored bounding box might
331 // have the wrong size (e.g. direct access)
332 // @todo: Which mesh should be used in order to determine the cluster centers:
333 // pro outMesh: we want perfectly fitting clusters around our output vertices
334 // however, this makes the cluster distribution/mapping dependent on the output
335 precice::mesh::BoundingBox localBB = inMesh->index().getRtreeBounds();
336
337#ifndef NDEBUG
338 // Safety check
339 precice::mesh::BoundingBox bb_check(inMesh->getDimensions());
340 for (const mesh::Vertex &vertex : inMesh->vertices()) {
341 bb_check.expandBy(vertex);
342 }
343 PRECICE_ASSERT(bb_check == localBB);
344#endif
345
346 // If we have very few vertices in the domain, (in this case twice our cluster
347 // size as we decompose most probably at least in 4 clusters) we just use a
348 // single cluster. The clustering result of the algorithm further down is in
349 // this case not optimal and might lead to too many clusters.
350 // The single cluster has in principle a radius of inf. We use here twice the
351 // length of the longest bounding box edge length and the center of the bounding
352 // box for the center point.
353 if (inMesh->nVertices() < verticesPerCluster * 2)
354 return {localBB.longestEdgeLength() * 2, Vertices{mesh::Vertex({localBB.center(), 0})}};
355
356 // We define a convenience alias for the localBB. In case we need to synchronize the clustering across ranks later on, we need
357 // to work with the global bounding box of the whole domain.
358 auto globalBB = localBB;
359
360 // Step 2: Now we pick random samples from the input mesh and ask the index tree for the k-nearest neighbors
361 // in order to estimate the point density and determine a proper cluster radius (see also the function documentation)
362 double clusterRadius = estimateClusterRadius(verticesPerCluster, inMesh, localBB);
363 PRECICE_DEBUG("Vertex cluster radius: {}", clusterRadius);
364
365 // maximum distance between cluster centers lying diagonal to each other. The maximum distance takes the overlap condition into
366 // account: example for 2D: if the distance between the centers is sqrt( 4 / 2 ) * radius, we violate the overlap condition between
367 // diagonal clusters
368 const int inDim = inMesh->getDimensions();
369 const double maximumCenterDistance = std::sqrt(4. / inDim) * clusterRadius * (1 - relativeOverlap);
370
371 // Step 3: using the maximum distance and the bounding box, compute the number of clusters in each direction
372 // we ceil the number of clusters in order to guarantee the desired overlap
373 std::array<unsigned int, 3> nClustersGlobal{1, 1, 1};
374 for (int d = 0; d < globalBB.getDimension(); ++d)
375 nClustersGlobal[d] = std::ceil(std::max(1., globalBB.getEdgeLength(d) / maximumCenterDistance));
376
377 // Step 4: Determine the centers of the clusters
378 Vertices centers;
379 // Vector used to temporarily store each center coordinates
380 std::vector<double> centerCoords(inDim);
381 // Vector storing the distances between the cluster in each direction
382 std::vector<double> distances(inDim);
383 // Vector storing the starting coordinates in each direction
384 std::vector<double> start(inDim);
385 // Fill the constant vectos, i.e., the distances and the start
386 for (unsigned int d = 0; d < distances.size(); ++d) {
387 // this distance calculation guarantees that we have at least the minimal specified overlap, as we ceil the division for the number of clusters
388 // as an alternative, once could use distance = const = maximumCenterDistance, which might lead to better results when projecting and removing duplicates (?)
389 distances[d] = globalBB.getEdgeLength(d) / (nClustersGlobal[d]);
390 start[d] = globalBB.minCorner()(d);
391 }
392 PRECICE_DEBUG("Distances: {}", distances);
393
394 // Step 5: Take care of the starting layer: if we start the grid at the globalBB edge we have an additional layer in each direction
395 if (startGridAtEdge) {
396 for (int d = 0; d < inDim; ++d) {
397 if (globalBB.getEdgeLength(d) > math::NUMERICAL_ZERO_DIFFERENCE) {
398 nClustersGlobal[d] += 1;
399 }
400 }
401 } // If we don't start at the edge, we shift the starting layer into the middle by 0.5 * distance
402 else {
403 std::transform(start.begin(), start.end(), distances.begin(), start.begin(), [](auto s, auto d) { return s + 0.5 * d; });
404 }
405
406 // Print some information
407 PRECICE_DEBUG("Global cluster distribution: {}", nClustersGlobal);
408 unsigned int nTotalClustersGlobal = std::accumulate(nClustersGlobal.begin(), nClustersGlobal.end(), 1U, std::multiplies<unsigned int>());
409 PRECICE_DEBUG("Global number of total clusters (tentative): {}", nTotalClustersGlobal);
410
411 // Step 6: transform the global metrics (number of clusters and starting layer) into local ones
412 // Since the local and global bounding boxes are the same in the current implementation, the
413 // global metrics and the local metrics will be the same.
414 // First, we need to update the starting points of our clustering scheme
415 std::array<unsigned int, 3> nClustersLocal{1, 1, 1};
416 for (unsigned int d = 0; d < start.size(); ++d) {
417 // Exclude the case where we have no distances (nClustersGlobal[d] == 1 )
418 if (distances[d] > 0) {
419 start[d] += std::ceil(((localBB.minCorner()(d) - start[d]) / distances[d]) - math::NUMERICAL_ZERO_DIFFERENCE) * distances[d];
420 // One cluster for the starting point and a further one for each distance we can fit into the BB
421 nClustersLocal[d] = 1 + std::floor(((localBB.maxCorner()(d) - start[d]) / distances[d]) + math::NUMERICAL_ZERO_DIFFERENCE);
422 }
423 }
424
425 unsigned int nTotalClustersLocal = std::accumulate(nClustersLocal.begin(), nClustersLocal.end(), 1U, std::multiplies<unsigned int>());
426 centers.resize(nTotalClustersLocal, mesh::Vertex({centerCoords, -1}));
427
428 // Print some information
429 PRECICE_DEBUG("Local cluster distribution: {}", nClustersLocal);
430 PRECICE_DEBUG("Local number of total clusters (tentative): {}", nTotalClustersLocal);
431
432 // Step 7: fill the center container using the zCurve for indexing. For the further processing, it is also important to ensure
433 // that the ID within the center Vertex class coincides with the position of the center vertex in the container.
434 // We start with the (bottom left) corner and iterate over all dimension
435 PRECICE_ASSERT(inDim >= 2);
436 if (inDim == 2) {
437 centerCoords[0] = start[0];
438 for (unsigned int x = 0; x < nClustersLocal[0]; ++x, centerCoords[0] += distances[0]) {
439 centerCoords[1] = start[1];
440 for (unsigned int y = 0; y < nClustersLocal[1]; ++y, centerCoords[1] += distances[1]) {
441 auto id = zCurve(std::array<unsigned int, 3>{{x, y, 0}}, nClustersLocal);
442 PRECICE_ASSERT(id < static_cast<int>(centers.size()) && id >= 0);
443 centers[id] = mesh::Vertex({centerCoords, id});
444 }
445 }
446 } else {
447 centerCoords[0] = start[0];
448 for (unsigned int x = 0; x < nClustersLocal[0]; ++x, centerCoords[0] += distances[0]) {
449 centerCoords[1] = start[1];
450 for (unsigned int y = 0; y < nClustersLocal[1]; ++y, centerCoords[1] += distances[1]) {
451 centerCoords[2] = start[2];
452 for (unsigned int z = 0; z < nClustersLocal[2]; ++z, centerCoords[2] += distances[2]) {
453 auto id = zCurve(std::array<unsigned int, 3>{{x, y, z}}, nClustersLocal);
454 PRECICE_ASSERT(id < static_cast<int>(centers.size()) && id >= 0);
455 centers[id] = mesh::Vertex({centerCoords, id});
456 }
457 }
458 }
459 }
460
461 // Step 8: tag vertices, which should be removed later on
462
463 // Step 8a: tag empty clusters
464 tagEmptyClusters(centers, clusterRadius, inMesh);
465 tagEmptyClusters(centers, clusterRadius, outMesh);
466 PRECICE_DEBUG("Number of non-tagged centers after empty clusters were filtered: {}", std::count_if(centers.begin(), centers.end(), [](auto &v) { return !v.isTagged(); }));
467
468 // Step 9: Move the cluster centers if desired
469 if (projectClustersToInput) {
470 projectClusterCentersToinputMesh(centers, inMesh);
471 // @todo: The duplication should probably be defined in terms of the target distance, not the actual distance. Find good default values
472 const auto duplicateThreshold = 0.4 * (*std::min_element(distances.begin(), distances.end()));
473 PRECICE_DEBUG("Tagging duplicates using the threshold value {} for the regular distances {}", duplicateThreshold, distances);
474 // usually inMesh->getDimensions() == 2
475 // to also cover the case where we have only a single layer in a 3D computation we use the number of clusters in z direction
476 // Step 8b or 9a: Moving cluster centers might lead to duplicate centers or centers being too close to each other. Therefore,
477 // we tag duplicate centers (see also the documentation of \ref tagDuplicateCenters ).
478 if (nClustersLocal[2] == 1) {
479 tagDuplicateCenters<2>(centers, nClustersLocal, duplicateThreshold);
480 } else {
481 tagDuplicateCenters<3>(centers, nClustersLocal, duplicateThreshold);
482 }
483 // the tagging here won't filter out a lot of clusters, but finding them anyway allows us to allocate the right amount of
484 // memory for the cluster vector in the mapping, which would otherwise be expensive
485 // we cannot hit any empty vertices in the inputmesh
486 tagEmptyClusters(centers, clusterRadius, outMesh);
487 PRECICE_DEBUG("Number of non-tagged centers after duplicate tagging : {}", std::count_if(centers.begin(), centers.end(), [](auto &v) { return !v.isTagged(); }));
488 }
489 PRECICE_ASSERT(std::all_of(centers.begin(), centers.end(), [idx = 0](auto &v) mutable { return v.getID() == idx++; }));
490
491 // Step 10: finally, remove the vertices from the center container
492 removeTaggedVertices(centers);
493 PRECICE_ASSERT(std::none_of(centers.begin(), centers.end(), [](auto &v) { return v.isTagged(); }));
494 PRECICE_CHECK(centers.size() > 0, "Too many vertices have been filtered out.");
495
496 return {clusterRadius, centers};
497}
498} // namespace impl
499} // namespace mapping
500} // namespace precice
#define PRECICE_DEBUG(...)
Definition LogMacros.hpp:64
#define PRECICE_TRACE(...)
Definition LogMacros.hpp:95
#define PRECICE_CHECK(check,...)
Definition LogMacros.hpp:35
T accumulate(T... args)
T all_of(T... args)
#define PRECICE_ASSERT(...)
Definition assertion.hpp:87
T begin(T... args)
T ceil(T... args)
This class provides a lightweight logger.
Definition Logger.hpp:16
An axis-aligned bounding box around a (partition of a) mesh.
Eigen::VectorXd maxCorner() const
the max corner of the bounding box
Eigen::VectorXd minCorner() const
the min corner of the bounding box
bool empty() const
Check if every dimension's length is equal to zero.
double longestEdgeLength() const
returns the maximum length of the bounding box in any dimension
void expandBy(const BoundingBox &otherBB)
Expand bounding box using another bounding box.
Eigen::VectorXd center() const
Returns the Center Of Gravity of the mesh.
double getEdgeLength(int axis) const
returns the edge length of a specific axis
Vertex of a mesh.
Definition Vertex.hpp:16
T count_if(T... args)
T emplace_back(T... args)
T end(T... args)
T floor(T... args)
T for_each(T... args)
T max_element(T... args)
T max(T... args)
T min_element(T... args)
std::tuple< double, Vertices > createClustering(mesh::PtrMesh inMesh, mesh::PtrMesh outMesh, double relativeOverlap, unsigned int verticesPerCluster, bool projectClustersToInput)
Creates a clustering as a collection of Vertices (representing the cluster centers) and a cluster rad...
double estimateClusterRadius(unsigned int verticesPerCluster, mesh::PtrMesh inMesh, const mesh::BoundingBox &bb)
Computes an estimate for the cluster radius, which results in approximately verticesPerCluster vertic...
std::vector< mesh::Vertex > Vertices
double computeSquaredDifference(const std::array< double, 3 > &u, std::array< double, 3 > v, const std::array< bool, 3 > &activeAxis={{true, true, true}})
Deletes all dead directions from fullVector and returns a vector of reduced dimensionality.
constexpr double NUMERICAL_ZERO_DIFFERENCE
Main namespace of the precice library.
int VertexID
Definition Types.hpp:13
T nth_element(T... args)
static precice::logging::Logger _log("precicec")
static std::unique_ptr< precice::Participant > impl
Definition preciceC.cpp:21
T remove_if(T... args)
T resize(T... args)
T size(T... args)
T sqrt(T... args)
T transform(T... args)