preCICE v3.1.2
Loading...
Searching...
No Matches
SphericalVertexCluster.hpp
Go to the documentation of this file.
1#pragma once
2
3#include <Eigen/Core>
4
5#include <boost/container/flat_set.hpp>
6
10#include "mesh/Filter.hpp"
12#include "profiling/Event.hpp"
13
14namespace precice {
15
16namespace mapping {
17
34template <typename RADIAL_BASIS_FUNCTION_T>
36public:
56 double radius,
57 RADIAL_BASIS_FUNCTION_T function,
58 Polynomial polynomial,
59 mesh::PtrMesh inputMesh,
60 mesh::PtrMesh outputMesh);
61
63 void mapConservative(const time::Sample &inData, Eigen::VectorXd &outData) const;
64
66 void mapConsistent(const time::Sample &inData, Eigen::VectorXd &outData) const;
67
69 void setNormalizedWeight(double normalizedWeight, VertexID vertexID);
70
72 double computeWeight(const mesh::Vertex &v) const;
73
75 unsigned int getNumberOfInputVertices() const;
76
79
81 void clear();
82
85 bool empty() const;
86
87private:
89 precice::logging::Logger _log{"mapping::SphericalVertexCluster"};
90
93
95 const double _radius;
96
99
100 // Stores the global IDs of the vertices so that we can apply a binary
101 // search in order to query specific objects. Here we have logarithmic
102 // complexity for setting the weights (only performed once), and constant
103 // complexity when traversing through the IDs (performed in each iteration).
104
106 boost::container::flat_set<VertexID> _inputIDs;
107
109 boost::container::flat_set<VertexID> _outputIDs;
110
113 Eigen::VectorXd _normalizedWeights;
114
117
120
123};
124
125// --------------------------------------------------- HEADER IMPLEMENTATIONS
126
127template <typename RADIAL_BASIS_FUNCTION_T>
129 mesh::Vertex center,
130 double radius,
131 RADIAL_BASIS_FUNCTION_T function,
132 Polynomial polynomial,
133 mesh::PtrMesh inputMesh,
134 mesh::PtrMesh outputMesh)
135 : _center(center), _radius(radius), _polynomial(polynomial), _weightingFunction(radius)
136{
138 precice::profiling::Event eq("map.pou.computeMapping.queryVertices");
139 // Disable integrated polynomial, as it might cause locally singular matrices
140 PRECICE_ASSERT(_polynomial != Polynomial::ON, "Integrated polynomial is not supported for partition of unity data mappings.");
141
142 // Get vertices to be mapped
143 // Subtract a safety margin to exclude the vertices at the edge
144 auto outIDs = outputMesh->index().getVerticesInsideBox(center, radius - math::NUMERICAL_ZERO_DIFFERENCE);
145 // Constructing the partition when we don't have evaluation points is pointless
146 auto inIDs = inputMesh->index().getVerticesInsideBox(center, radius);
147
148 // Transform the vector to the appropriate boost data structure
149 // The IDs are sorted in the boost flat_set, hence, the function here has N log(N) complexity
150 _inputIDs.insert(inIDs.begin(), inIDs.end());
151 _outputIDs.insert(outIDs.begin(), outIDs.end());
152 eq.stop();
153 // If the cluster is empty, we return immediately
154 if (empty()) {
155 return;
156 }
157
158 PRECICE_DEBUG("SphericalVertexCluster input size: {}", inIDs.size());
159 PRECICE_DEBUG("SphericalVertexCluster output size: {}", outIDs.size());
160
161 // The polynomial system is underdetermined if inIDs.size() < dimension + 1. However, the dynamic adoption of the axis in the RBF solver
162 // disables axis, if necessary. Hence, we don't disable the complete polynomial here for underdetermined systems. The case should anyway
163 // only occur for an almost unreasonable small vertices-per-cluster configuration.
164
165 // Construct the solver. Here, the constructor of the RadialBasisFctSolver computes already the decompositions etc, such that we can mark the
166 // mapping in this cluster as computed (mostly for debugging purpose)
167 std::vector<bool> deadAxis(inputMesh->getDimensions(), false);
168 precice::profiling::Event e("map.pou.computeMapping.rbfSolver");
169 _rbfSolver = RadialBasisFctSolver<RADIAL_BASIS_FUNCTION_T>{function, *inputMesh.get(), _inputIDs, *outputMesh.get(), _outputIDs, deadAxis, _polynomial};
170 _hasComputedMapping = true;
171}
172
173template <typename RADIAL_BASIS_FUNCTION_T>
175{
176 PRECICE_ASSERT(_outputIDs.size() > 0);
177 PRECICE_ASSERT(_outputIDs.contains(id), id);
178 PRECICE_ASSERT(normalizedWeight > 0);
179
180 if (_normalizedWeights.size() == 0)
181 _normalizedWeights.resize(_outputIDs.size());
182
183 // The find method of boost flat_set comes with O(log(N)) complexity (the more expensive part here)
184 auto localID = _outputIDs.index_of(_outputIDs.find(id));
185
186 PRECICE_ASSERT(static_cast<Eigen::Index>(localID) < _normalizedWeights.size(), localID, _normalizedWeights.size());
187 _normalizedWeights[localID] = normalizedWeight;
188}
189
190template <typename RADIAL_BASIS_FUNCTION_T>
192{
193 // First, a few sanity checks. Empty partitions shouldn't be stored at all
194 PRECICE_ASSERT(!empty());
195 PRECICE_ASSERT(_hasComputedMapping);
196 PRECICE_ASSERT(_normalizedWeights.size() == static_cast<Eigen::Index>(_outputIDs.size()));
197
198 // Define an alias for data dimension in order to avoid ambiguity
199 const unsigned int nComponents = inData.dataDims;
200 const auto & localInData = inData.values;
201
202 // TODO: We can probably reduce the temporary allocations here
203 Eigen::VectorXd in(_rbfSolver.getOutputSize());
204
205 // Now we perform the data mapping component-wise
206 for (unsigned int c = 0; c < nComponents; ++c) {
207 // Step 1: extract the relevant input data from the global input data and store
208 // it in a contiguous array, which is required for the RBF solver
209 for (unsigned int i = 0; i < _outputIDs.size(); ++i) {
210 const auto dataIndex = *(_outputIDs.nth(i));
211 PRECICE_ASSERT(dataIndex * nComponents + c < localInData.size(), dataIndex * nComponents + c, localInData.size());
212 PRECICE_ASSERT(_normalizedWeights[i] > 0, _normalizedWeights[i], i);
213 // here, we also directly apply the weighting, i.e., we split the input data
214 in[i] = localInData[dataIndex * nComponents + c] * _normalizedWeights[i];
215 }
216
217 // Step 2: solve the system using a conservative constraint
218 auto result = _rbfSolver.solveConservative(in, _polynomial);
219 PRECICE_ASSERT(result.size() == static_cast<Eigen::Index>(_inputIDs.size()));
220
221 // Step 3: now accumulate the result into our global output data
222 for (unsigned int i = 0; i < _inputIDs.size(); ++i) {
223 const auto dataIndex = *(_inputIDs.nth(i));
224 PRECICE_ASSERT(dataIndex * nComponents + c < outData.size(), dataIndex * nComponents + c, outData.size());
225 outData[dataIndex * nComponents + c] += result(i);
226 }
227 }
228}
229
230template <typename RADIAL_BASIS_FUNCTION_T>
231void SphericalVertexCluster<RADIAL_BASIS_FUNCTION_T>::mapConsistent(const time::Sample &inData, Eigen::VectorXd &outData) const
232{
233 // First, a few sanity checks. Empty partitions shouldn't be stored at all
234 PRECICE_ASSERT(!empty());
235 PRECICE_ASSERT(_hasComputedMapping);
236 PRECICE_ASSERT(_normalizedWeights.size() == static_cast<Eigen::Index>(_outputIDs.size()));
237
238 // Define an alias for data dimension in order to avoid ambiguity
239 const unsigned int nComponents = inData.dataDims;
240 const auto & localInData = inData.values;
241
242 Eigen::VectorXd in(_rbfSolver.getInputSize());
243
244 // Now we perform the data mapping component-wise
245 for (unsigned int c = 0; c < nComponents; ++c) {
246 // Step 1: extract the relevant input data from the global input data and store
247 // it in a contiguous array, which is required for the RBF solver (last polyparams entries remain zero)
248 for (unsigned int i = 0; i < _inputIDs.size(); i++) {
249 const auto dataIndex = *(_inputIDs.nth(i));
250 PRECICE_ASSERT(dataIndex * nComponents + c < localInData.size(), dataIndex * nComponents + c, localInData.size());
251 in[i] = localInData[dataIndex * nComponents + c];
252 }
253
254 // Step 2: solve the system using a consistent constraint
255 auto result = _rbfSolver.solveConsistent(in, _polynomial);
256 PRECICE_ASSERT(static_cast<Eigen::Index>(_outputIDs.size()) == result.size());
257
258 // Step 3: now accumulate the result into our global output data
259 for (unsigned int i = 0; i < _outputIDs.size(); ++i) {
260 const auto dataIndex = *(_outputIDs.nth(i));
261 PRECICE_ASSERT(dataIndex * nComponents + c < outData.size(), dataIndex * nComponents + c, outData.size());
262 PRECICE_ASSERT(_normalizedWeights[i] > 0);
263 // here, we also directly apply the weighting, i.e., split the result data
264 outData[dataIndex * nComponents + c] += result(i) * _normalizedWeights[i];
265 }
266 }
267}
268
269template <typename RADIAL_BASIS_FUNCTION_T>
271{
272 // Assume that the local interpolant and the weighting function are the same
273 // TODO: We don't need to reduce the dead coordinates here as the values should reduce anyway
274 auto res = computeSquaredDifference(_center.rawCoords(), v.rawCoords(), {{true, true, true}});
275 return _weightingFunction.evaluate(std::sqrt(res));
276}
277
278template <typename RADIAL_BASIS_FUNCTION_T>
280{
281 return _inputIDs.size();
282}
283
284template <typename RADIAL_BASIS_FUNCTION_T>
289
290template <typename RADIAL_BASIS_FUNCTION_T>
292{
293 return _outputIDs.size() == 0 || _inputIDs.size() == 0;
294}
295
296template <typename RADIAL_BASIS_FUNCTION_T>
298{
300 _inputIDs.clear();
301 _outputIDs.clear();
302 _rbfSolver.clear();
303 _hasComputedMapping = false;
304}
305} // namespace mapping
306} // namespace precice
#define PRECICE_DEBUG(...)
Definition LogMacros.hpp:64
#define PRECICE_TRACE(...)
Definition LogMacros.hpp:95
#define PRECICE_ASSERT(...)
Definition assertion.hpp:87
This class provides a lightweight logger.
Definition Logger.hpp:16
Wendland radial basis function with compact support.
boost::container::flat_set< VertexID > _outputIDs
Global VertexIDs associated to the output mesh (see constructor)
unsigned int getNumberOfInputVertices() const
Number of input vertices this partition operates on.
mesh::Vertex _center
center vertex of the cluster
const double _radius
radius of the vertex cluster
std::array< double, 3 > getCenterCoords() const
The center coordinate of this cluster.
double computeWeight(const mesh::Vertex &v) const
Compute the weight for a given vertex.
bool _hasComputedMapping
Boolean switch in order to indicate that a mapping was computed.
CompactPolynomialC2 _weightingFunction
The weighting function.
void setNormalizedWeight(double normalizedWeight, VertexID vertexID)
Set the normalized weight for the given vertexID in the outputMesh.
Polynomial _polynomial
Polynomial treatment in the RBF solver.
RadialBasisFctSolver< RADIAL_BASIS_FUNCTION_T > _rbfSolver
The RBF solver.
void mapConsistent(const time::Sample &inData, Eigen::VectorXd &outData) const
Evaluates a consistent mapping and agglomerates the result in the given output data.
void clear()
Invalidates and erases data structures the cluster holds.
SphericalVertexCluster(mesh::Vertex center, double radius, RADIAL_BASIS_FUNCTION_T function, Polynomial polynomial, mesh::PtrMesh inputMesh, mesh::PtrMesh outputMesh)
boost::container::flat_set< VertexID > _inputIDs
Global VertexIDs associated to the input mesh (see constructor)
void mapConservative(const time::Sample &inData, Eigen::VectorXd &outData) const
Evaluates a conservative mapping and agglomerates the result in the given output data.
precice::logging::Logger _log
logger, as usual
Vertex of a mesh.
Definition Vertex.hpp:16
const RawCoords & rawCoords() const
Direct access to the coordinates.
Definition Vertex.hpp:123
Eigen::VectorXd getCoords() const
Returns the coordinates of the vertex.
Definition Vertex.hpp:116
void stop()
Stops a running event.
Definition Event.cpp:37
T get(T... args)
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.
Polynomial
How to handle the polynomial?
constexpr double NUMERICAL_ZERO_DIFFERENCE
Main namespace of the precice library.
int VertexID
Definition Types.hpp:13
T sqrt(T... args)
int dataDims
The dimensionality of the data.
Definition Sample.hpp:45
Eigen::VectorXd values
Definition Sample.hpp:49