preCICE v3.1.2
Loading...
Searching...
No Matches
PartitionOfUnityMapping.hpp
Go to the documentation of this file.
1#pragma once
2
3#include <Eigen/Core>
4#include <numeric>
5
7#include "io/ExportVTU.hpp"
10#include "mesh/Filter.hpp"
12#include "profiling/Event.hpp"
13#include "query/Index.hpp"
14#include "utils/IntraComm.hpp"
15
16namespace precice {
17extern bool syncMode;
18
19namespace mapping {
20
28template <typename RADIAL_BASIS_FUNCTION_T>
30public:
47 Mapping::Constraint constraint,
48 int dimension,
49 RADIAL_BASIS_FUNCTION_T function,
50 Polynomial polynomial,
51 unsigned int verticesPerCluster,
52 double relativeOverlap,
53 bool projectToInput);
54
64 void computeMapping() final override;
65
67 void clear() final override;
68
70 void tagMeshFirstRound() final override;
71
73 void tagMeshSecondRound() final override;
74
76 std::string getName() const final override;
77
78private:
80 precice::logging::Logger _log{"mapping::PartitionOfUnityMapping"};
81
84
86 RADIAL_BASIS_FUNCTION_T _basisFunction;
87
89
91 const unsigned int _verticesPerCluster;
92
94 const double _relativeOverlap;
95
97 const bool _projectToInput;
98
100 double _clusterRadius = 0;
101
104
106 virtual void mapConservative(const time::Sample &inData, Eigen::VectorXd &outData) override;
107
109 virtual void mapConsistent(const time::Sample &inData, Eigen::VectorXd &outData) override;
110
114};
115
116template <typename RADIAL_BASIS_FUNCTION_T>
118 Mapping::Constraint constraint,
119 int dimension,
120 RADIAL_BASIS_FUNCTION_T function,
121 Polynomial polynomial,
122 unsigned int verticesPerCluster,
123 double relativeOverlap,
124 bool projectToInput)
125 : Mapping(constraint, dimension, false, Mapping::InitialGuessRequirement::None),
126 _basisFunction(function), _verticesPerCluster(verticesPerCluster), _relativeOverlap(relativeOverlap), _projectToInput(projectToInput), _polynomial(polynomial)
127{
128 PRECICE_ASSERT(this->getDimensions() <= 3);
129 PRECICE_ASSERT(_polynomial != Polynomial::ON, "Integrated polynomial is not supported for partition of unity data mappings.");
130 PRECICE_ASSERT(_relativeOverlap < 1, "The relative overlap has to be smaller than one.");
131 PRECICE_ASSERT(_verticesPerCluster > 0, "The number of vertices per cluster has to be greater zero.");
132
133 if (isScaledConsistent()) {
136 } else {
139 }
140}
141
142template <typename RADIAL_BASIS_FUNCTION_T>
144{
146
147 precice::profiling::Event e("map.pou.computeMapping.From" + this->input()->getName() + "To" + this->output()->getName(), profiling::Synchronize);
148
149 // Recompute the whole clustering
150 PRECICE_ASSERT(!this->_hasComputedMapping, "Please clear the mapping before recomputing.");
151
152 mesh::PtrMesh inMesh;
153 mesh::PtrMesh outMesh;
154 if (this->hasConstraint(Mapping::CONSERVATIVE)) {
155 inMesh = this->output();
156 outMesh = this->input();
157 } else { // Consistent or scaled consistent
158 inMesh = this->input();
159 outMesh = this->output();
160 }
161
162 precice::profiling::Event eClusters("map.pou.computeMapping.createClustering.From" + this->input()->getName() + "To" + this->output()->getName());
163 // Step 1: get a tentative clustering consisting of centers and a radius from one of the available algorithms
164 auto [clusterRadius, centerCandidates] = impl::createClustering(inMesh, outMesh, _relativeOverlap, _verticesPerCluster, _projectToInput);
165 eClusters.stop();
166
167 _clusterRadius = clusterRadius;
168 PRECICE_ASSERT(_clusterRadius > 0 || inMesh->nVertices() == 0 || outMesh->nVertices() == 0);
169
170 // Step 2: check, which of the resulting clusters are non-empty and register the cluster centers in a mesh
171 // Here, the VertexCluster computes the matrix decompositions directly in case the cluster is non-empty
172 mesh::Mesh centerMesh("pou-centers-" + inMesh->getName(), this->getDimensions(), mesh::Mesh::MESH_ID_UNDEFINED);
173 auto & meshVertices = centerMesh.vertices();
174
175 meshVertices.clear();
176 _clusters.clear();
177 _clusters.reserve(centerCandidates.size());
178 for (const auto &c : centerCandidates) {
179 // We cannot simply copy the vertex from the container in order to fill the vertices of the centerMesh, as the vertexID of each center needs to match the index
180 // of the cluster within the _clusters vector. That's required for the indexing further down and asserted below
181 const VertexID vertexID = meshVertices.size();
182 mesh::Vertex center(c.getCoords(), vertexID);
183 SphericalVertexCluster<RADIAL_BASIS_FUNCTION_T> cluster(center, _clusterRadius, _basisFunction, _polynomial, inMesh, outMesh);
184
185 // Consider only non-empty clusters (more of a safeguard here)
186 if (!cluster.empty()) {
187 PRECICE_ASSERT(center.getID() == static_cast<int>(_clusters.size()), center.getID(), _clusters.size());
188 meshVertices.emplace_back(std::move(center));
189 _clusters.emplace_back(std::move(cluster));
190 }
191 }
192
193 e.addData("n clusters", _clusters.size());
194 // Log the average number of resulting clusters
195 PRECICE_DEBUG("Partition of unity data mapping between mesh \"{}\" and mesh \"{}\": mesh \"{}\" on rank {} was decomposed into {} clusters.", this->input()->getName(), this->output()->getName(), inMesh->getName(), utils::IntraComm::getRank(), _clusters.size());
196
197 if (_clusters.size() > 0) {
198 PRECICE_DEBUG("Average number of vertices per cluster {}", std::accumulate(_clusters.begin(), _clusters.end(), static_cast<unsigned int>(0), [](auto &acc, auto &val) { return acc += val.getNumberOfInputVertices(); }) / _clusters.size());
199 PRECICE_DEBUG("Maximum number of vertices per cluster {}", std::max_element(_clusters.begin(), _clusters.end(), [](auto &v1, auto &v2) { return v1.getNumberOfInputVertices() < v2.getNumberOfInputVertices(); })->getNumberOfInputVertices());
200 PRECICE_DEBUG("Minimum number of vertices per cluster {}", std::min_element(_clusters.begin(), _clusters.end(), [](auto &v1, auto &v2) { return v1.getNumberOfInputVertices() < v2.getNumberOfInputVertices(); })->getNumberOfInputVertices());
201 }
202
203 precice::profiling::Event eWeights("map.pou.computeMapping.computeWeights");
204 // Log a bounding box of the center mesh
205 centerMesh.computeBoundingBox();
206 PRECICE_DEBUG("Bounding Box of the cluster centers {}", centerMesh.getBoundingBox());
207
208 // Step 3: index the clusters / the center mesh in order to define the output vertex -> cluster ownership
209 // the ownership is required to compute the normalized partition of unity weights (Step 4)
210 query::Index clusterIndex(centerMesh);
211 // Step 4: find all clusters the output vertex lies in, i.e., find all cluster centers which have the distance of a cluster radius from the given output vertex
212 // Here, we do this using the RTree on the centerMesh: VertexID (queried from the centersMesh) == clusterID, by construction above. The loop uses
213 // the vertices to compute the weights required for the partition of unity data mapping.
214 // Note: this could also be done on-the-fly in the map data phase for dynamic queries, which would require to make the mesh as well as the indexTree member variables.
215 PRECICE_DEBUG("Computing cluster-vertex association");
216 for (const auto &vertex : outMesh->vertices()) {
217 // Step 4a: get the relevant clusters for the output vertex
218 auto clusterIDs = clusterIndex.getVerticesInsideBox(vertex, _clusterRadius);
219 const auto localNumberOfClusters = clusterIDs.size();
220
221 // Consider the case where we didn't find any cluster (meshes don't match very well)
222 //
223 // In principle, we could assign the vertex to the closest cluster using clusterIDs.emplace_back(clusterIndex.getClosestVertex(vertex.getCoords()).index);
224 // However, this leads to a conflict with weights already set in the corresponding cluster, since we insert the ID and, later on, map the ID to a local weight index
225 // Of course, we could rearrange the weights, but we want to avoid the case here anyway, i.e., prefer to abort.
226 PRECICE_CHECK(localNumberOfClusters > 0,
227 "Output vertex {} of mesh \"{}\" could not be assigned to any cluster in the rbf-pum mapping. This probably means that the meshes do not match well geometry-wise: Visualize the exported preCICE meshes to confirm."
228 " If the meshes are fine geometry-wise, you can try to increase the number of \"vertices-per-cluster\" (default is 50), the \"relative-overlap\" (default is 0.15),"
229 " or disable the option \"project-to-input\"."
230 "These options are only valid for the <mapping:rbf-pum-direct/> tag.",
231 vertex.getCoords(), outMesh->getName());
232
233 // Next we compute the normalized weights of each output vertex for each partition
234 PRECICE_ASSERT(localNumberOfClusters > 0, "No cluster found for vertex {}", vertex.getCoords());
235
236 // Step 4b: compute the weight in each partition individually and store them in 'weights'
237 std::vector<double> weights(localNumberOfClusters);
238 std::transform(clusterIDs.cbegin(), clusterIDs.cend(), weights.begin(), [&](const auto &ids) { return _clusters[ids].computeWeight(vertex); });
239 double weightSum = std::accumulate(weights.begin(), weights.end(), static_cast<double>(0.));
240 // TODO: This covers the edge case of vertices being at the edge of (several) clusters
241 // In case the sum is equal to zero, we assign equal weights for all clusters
242 if (weightSum <= 0) {
243 PRECICE_ASSERT(weights.size() > 0);
244 std::for_each(weights.begin(), weights.end(), [&weights](auto &w) { w = 1. / weights.size(); });
245 weightSum = 1;
246 }
247 PRECICE_ASSERT(weightSum > 0);
248
249 // Step 4c: scale the weight using the weight sum and store the normalized weight in all associated clusters
250 for (unsigned int i = 0; i < localNumberOfClusters; ++i) {
251 PRECICE_ASSERT(clusterIDs[i] < static_cast<int>(_clusters.size()));
252 _clusters[clusterIDs[i]].setNormalizedWeight(weights[i] / weightSum, vertex.getID());
253 }
254 }
255 eWeights.stop();
256
257 // Uncomment to add a VTK export of the cluster center distribution for visualization purposes
258 // exportClusterCentersAsVTU(centerMesh);
259
260 this->_hasComputedMapping = true;
261}
262
263template <typename RADIAL_BASIS_FUNCTION_T>
265{
267
268 precice::profiling::Event e("map.pou.mapData.From" + input()->getName() + "To" + output()->getName(), profiling::Synchronize);
269
270 // Execute the actual mapping evaluation in all clusters
271 // 1. Assert that all output data values were reset, as we accumulate data in all clusters independently
272 PRECICE_ASSERT(outData.isZero());
273
274 // 2. Iterate over all clusters and accumulate the result in the output data
275 std::for_each(_clusters.begin(), _clusters.end(), [&](auto &cluster) { cluster.mapConservative(inData, outData); });
276}
277
278template <typename RADIAL_BASIS_FUNCTION_T>
280{
282
283 precice::profiling::Event e("map.pou.mapData.From" + input()->getName() + "To" + output()->getName(), profiling::Synchronize);
284
285 // Execute the actual mapping evaluation in all clusters
286 // 1. Assert that all output data values were reset, as we accumulate data in all clusters independently
287 PRECICE_ASSERT(outData.isZero());
288
289 // 2. Execute the actual mapping evaluation in all vertex clusters and accumulate the data
290 std::for_each(_clusters.begin(), _clusters.end(), [&](auto &clusters) { clusters.mapConsistent(inData, outData); });
291}
292
293template <typename RADIAL_BASIS_FUNCTION_T>
295{
297 mesh::PtrMesh filterMesh, outMesh;
298 if (this->hasConstraint(Mapping::CONSERVATIVE)) {
299 filterMesh = this->output(); // remote
300 outMesh = this->input(); // local
301 } else {
302 filterMesh = this->input(); // remote
303 outMesh = this->output(); // local
304 }
305
306 if (outMesh->empty())
307 return; // Ranks not at the interface should never hold interface vertices
308
309 // Note that we don't use the corresponding bounding box functions from
310 // precice::mesh (e.g. ::getBoundingBox), as the stored bounding box might
311 // have the wrong size (e.g. direct access)
312 precice::mesh::BoundingBox bb = filterMesh->index().getRtreeBounds();
313
314#ifndef NDEBUG
315 // Safety check
316 precice::mesh::BoundingBox bb_check(filterMesh->getDimensions());
317 for (const mesh::Vertex &vertex : filterMesh->vertices()) {
318 bb_check.expandBy(vertex);
319 }
320 PRECICE_ASSERT(bb_check == bb);
321#endif
322 // @TODO: This is assert and the function might run into problems if we have only a single.
323 // vertex in the received mesh
324 // However, with the current BB implementation, the expandBy function will just do nothing.
325 PRECICE_ASSERT(!bb.empty());
326 // This function behaves differently when disabling the geometric filtering
327 // without filter, we look at the complete mesh, whereas otherwise, we look at
328 // a fraction of the mesh and we might end up with too few vertices per rank
329 // We cannot prevent too few vertices from being tagged, if we have filtered too much
330 // vertices, but the user could increase the safety-factor or disable the filtering.
331 // When no geometric filter is applid, vertices().size() is here the same as
332 // getGlobalNumberOfVertices
333 if (filterMesh->nVertices() < _verticesPerCluster &&
334 filterMesh->nVertices() < static_cast<std::size_t>(filterMesh->getGlobalNumberOfVertices())) {
335 PRECICE_WARN("The repartitioning of the received mesh \"{}\" resulted in {} vertices on this "
336 "rank, which is less than the desired number of vertices per cluster configured "
337 "in the partition of unity mapping ({}). Consider increasing the safety-factor "
338 "or switching off the geometric filter (<receive-mesh: ... geometric-filter=\"no-filter\" .../>)",
339 filterMesh->getName(), filterMesh->nVertices(), _verticesPerCluster);
340 }
341
342 if (_clusterRadius == 0)
343 _clusterRadius = impl::estimateClusterRadius(_verticesPerCluster, filterMesh, bb);
344
345 PRECICE_DEBUG("Cluster radius estimate: {}", _clusterRadius);
346 PRECICE_ASSERT(_clusterRadius > 0);
347
348 // Get the local bounding boxes
349 auto localBB = outMesh->getBoundingBox();
350 // Now we extend the bounding box by the radius
351 localBB.expandBy(2 * _clusterRadius);
352
353 // ... and tag all affected vertices
354 auto verticesNew = filterMesh->index().getVerticesInsideBox(localBB);
355
356 std::for_each(verticesNew.begin(), verticesNew.end(), [&filterMesh](VertexID v) { filterMesh->vertex(v).tag(); });
357}
358
359template <typename RADIAL_BASIS_FUNCTION_T>
361{
362 // Nothing to be done here. There is no global ownership for matrix entries required and we tag all potentially locally relevant vertices already in the first round.
363}
364
365template <typename RADIAL_BASIS_FUNCTION_T>
367{
369
370 auto dataRadius = centerMesh.createData("radius", 1, -1);
371 auto dataCardinality = centerMesh.createData("number-of-vertices", 1, -1);
372 centerMesh.allocateDataValues();
373 dataRadius->values().fill(_clusterRadius);
374 for (unsigned int i = 0; i < _clusters.size(); ++i) {
375 dataCardinality->values()[i] = static_cast<double>(_clusters[i].getNumberOfInputVertices());
376 }
377
378 // We have to create the global offsets in order to export things in parallel
380 // send number of vertices
381 PRECICE_DEBUG("Send number of vertices: {}", centerMesh.nVertices());
382 int numberOfVertices = centerMesh.nVertices();
383 utils::IntraComm::getCommunication()->send(numberOfVertices, 0);
384
385 // receive vertex offsets
386 mesh::Mesh::VertexOffsets vertexOffsets;
387 utils::IntraComm::getCommunication()->broadcast(vertexOffsets, 0);
388 PRECICE_DEBUG("Vertex offsets: {}", vertexOffsets);
389 PRECICE_ASSERT(centerMesh.getVertexOffsets().empty());
390 centerMesh.setVertexOffsets(std::move(vertexOffsets));
391 } else if (utils::IntraComm::isPrimary()) {
392
394 vertexOffsets[0] = centerMesh.nVertices();
395
396 // receive number of secondary vertices and fill vertex offsets
397 for (int secondaryRank : utils::IntraComm::allSecondaryRanks()) {
398 int numberOfSecondaryRankVertices = -1;
399 utils::IntraComm::getCommunication()->receive(numberOfSecondaryRankVertices, secondaryRank);
400 PRECICE_ASSERT(numberOfSecondaryRankVertices >= 0);
401 vertexOffsets[secondaryRank] = numberOfSecondaryRankVertices + vertexOffsets[secondaryRank - 1];
402 }
403
404 // broadcast vertex offsets
405 PRECICE_DEBUG("Vertex offsets: {}", centerMesh.getVertexOffsets());
406 utils::IntraComm::getCommunication()->broadcast(vertexOffsets);
407 centerMesh.setVertexOffsets(std::move(vertexOffsets));
408 }
409
410 io::ExportVTU exporter;
411 exporter.doExport(centerMesh.getName(), "exports", centerMesh);
412}
413
414template <typename RADIAL_BASIS_FUNCTION_T>
416{
418 _clusters.clear();
419 // TODO: Don't reset this here
420 _clusterRadius = 0;
421 this->_hasComputedMapping = false;
422}
423
424template <typename RADIAL_BASIS_FUNCTION_T>
426{
427 return "partition-of-unity RBF";
428}
429} // namespace mapping
430} // namespace precice
#define PRECICE_WARN(...)
Definition LogMacros.hpp:11
#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)
#define PRECICE_ASSERT(...)
Definition assertion.hpp:87
T begin(T... args)
void doExport(const std::string &name, const std::string &location, const mesh::Mesh &mesh) override
Does export. Has to be implemented in subclass.
Definition ExportXML.cpp:23
This class provides a lightweight logger.
Definition Logger.hpp:16
Abstract base class for mapping of data from one mesh to another.
Definition Mapping.hpp:15
Constraint
Specifies additional constraints for a mapping.
Definition Mapping.hpp:29
bool isScaledConsistent() const
Returns true if mapping is a form of scaled consistent mapping.
Definition Mapping.cpp:257
void setInputRequirement(MeshRequirement requirement)
Sets the mesh requirement for the input mesh.
Definition Mapping.cpp:96
void setOutputRequirement(MeshRequirement requirement)
Sets the mesh requirement for the output mesh.
Definition Mapping.cpp:102
InitialGuessRequirement
Specifies whether the mapping requires an initial guess.
Definition Mapping.hpp:63
void tagMeshSecondRound() final override
nothing to do here
void tagMeshFirstRound() final override
tag the vertices required for the mapping
Polynomial _polynomial
polynomial treatment of the RBF system
precice::logging::Logger _log
logger, as usual
const double _relativeOverlap
overlap of vertex clusters
std::vector< SphericalVertexCluster< RADIAL_BASIS_FUNCTION_T > > _clusters
main data container storing all the clusters, which need to be solved individually
const bool _projectToInput
toggles whether we project the cluster centers to the input mesh
std::string getName() const final override
name of the pum mapping
virtual void mapConsistent(const time::Sample &inData, Eigen::VectorXd &outData) override
Maps data using a consistent constraint.
double _clusterRadius
derived parameter based on the input above: the radius of each cluster
PartitionOfUnityMapping(Mapping::Constraint constraint, int dimension, RADIAL_BASIS_FUNCTION_T function, Polynomial polynomial, unsigned int verticesPerCluster, double relativeOverlap, bool projectToInput)
const unsigned int _verticesPerCluster
Input parameters provided by the user for the clustering algorithm:
RADIAL_BASIS_FUNCTION_T _basisFunction
Radial basis function type used in interpolation.
virtual void mapConservative(const time::Sample &inData, Eigen::VectorXd &outData) override
Maps data using a conservative constraint.
void clear() final override
Clears a computed mapping by deleting the content of the _clusters vector.
An axis-aligned bounding box around a (partition of a) mesh.
bool empty() const
Check if every dimension's length is equal to zero.
void expandBy(const BoundingBox &otherBB)
Expand bounding box using another bounding box.
Container and creator for meshes.
Definition Mesh.hpp:39
VertexContainer & vertices()
Returns modifieable container holding all vertices.
Definition Mesh.cpp:53
const std::string & getName() const
Returns the name of the mesh, as set in the config file.
Definition Mesh.cpp:218
static constexpr MeshID MESH_ID_UNDEFINED
Use if the id of the mesh is not necessary.
Definition Mesh.hpp:58
std::size_t nVertices() const
Returns the number of vertices.
Definition Mesh.cpp:63
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
const VertexOffsets & getVertexOffsets() const
Definition Mesh.hpp:248
void computeBoundingBox()
Computes the boundingBox for the vertices.
Definition Mesh.cpp:242
void setVertexOffsets(VertexOffsets vertexOffsets)
Only used for tests.
Definition Mesh.hpp:254
const BoundingBox & getBoundingBox() const
Returns the bounding box of the mesh.
Definition Mesh.cpp:365
void allocateDataValues()
Allocates memory for the vertex data values and corresponding gradient values.
Definition Mesh.cpp:233
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
void stop()
Stops a running event.
Definition Event.cpp:37
void addData(std::string_view key, int value)
Adds named integer data, associated to an event.
Definition Event.cpp:48
Class to query the index trees of the mesh.
Definition Index.hpp:65
std::vector< VertexID > getVerticesInsideBox(const mesh::Vertex &centerVertex, double radius)
Return all the vertices inside the box formed by vertex and radius (boundary exclusive)
Definition Index.cpp:223
static int getSize()
Number of ranks. This includes ranks from both participants, e.g. minimal size is 2.
Definition IntraComm.cpp:47
static Rank getRank()
Current rank.
Definition IntraComm.cpp:42
static bool isPrimary()
True if this process is running the primary rank.
Definition IntraComm.cpp:52
static auto allSecondaryRanks()
Returns an iterable range over salve ranks [1, _size)
Definition IntraComm.hpp:37
static bool isSecondary()
True if this process is running a secondary rank.
Definition IntraComm.cpp:57
static com::PtrCommunication & getCommunication()
Intra-participant communication.
Definition IntraComm.hpp:31
T clear(T... args)
T empty(T... args)
T end(T... args)
T for_each(T... args)
T max_element(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...
Polynomial
How to handle the polynomial?
static constexpr SynchronizeTag Synchronize
Convenience instance of the SynchronizeTag.
Definition Event.hpp:21
Main namespace of the precice library.
int VertexID
Definition Types.hpp:13
bool syncMode
Enabled further inter- and intra-solver synchronisation.
T size(T... args)
T transform(T... args)