preCICE v3.1.2
Loading...
Searching...
No Matches
RadialBasisFctSolver.hpp
Go to the documentation of this file.
1#pragma once
2
3#include <Eigen/Cholesky>
4#include <Eigen/QR>
5#include <Eigen/SVD>
6#include <boost/range/adaptor/indexed.hpp>
7#include <boost/range/irange.hpp>
8#include <numeric>
10#include "mesh/Mesh.hpp"
12
13namespace precice {
14namespace mapping {
15
22template <typename RADIAL_BASIS_FUNCTION_T>
24public:
25 using DecompositionType = std::conditional_t<RADIAL_BASIS_FUNCTION_T::isStrictlyPositiveDefinite(), Eigen::LLT<Eigen::MatrixXd>, Eigen::ColPivHouseholderQR<Eigen::MatrixXd>>;
26 using BASIS_FUNCTION_T = RADIAL_BASIS_FUNCTION_T;
29
37 template <typename IndexContainer>
38 RadialBasisFctSolver(RADIAL_BASIS_FUNCTION_T basisFunction, const mesh::Mesh &inputMesh, const IndexContainer &inputIDs,
39 const mesh::Mesh &outputMesh, const IndexContainer &outputIDs, std::vector<bool> deadAxis, Polynomial polynomial);
40
42 Eigen::VectorXd solveConsistent(Eigen::VectorXd &inputData, Polynomial polynomial) const;
43
45 Eigen::VectorXd solveConservative(const Eigen::VectorXd &inputData, Polynomial polynomial) const;
46
47 // Clear all stored matrices
48 void clear();
49
50 // Returns the size of the input data
51 Eigen::Index getInputSize() const;
52
53 // Returns the size of the input data
54 Eigen::Index getOutputSize() const;
55
56private:
57 precice::logging::Logger _log{"mapping::RadialBasisFctSolver"};
58
61
63 Eigen::ColPivHouseholderQR<Eigen::MatrixXd> _qrMatrixQ;
64
66 Eigen::MatrixXd _matrixQ;
67
69 Eigen::MatrixXd _matrixV;
70
72 Eigen::MatrixXd _matrixA;
73};
74
75// ------- Non-Member Functions ---------
76
79 const std::array<double, 3> &u,
81 const std::array<bool, 3> & activeAxis = {{true, true, true}})
82{
83 // Subtract the values and multiply out dead dimensions
84 for (unsigned int d = 0; d < v.size(); ++d) {
85 v[d] = (u[d] - v[d]) * static_cast<int>(activeAxis[d]);
86 }
87 // @todo: this can be replaced by std::hypot when moving to C++17
88 return std::accumulate(v.begin(), v.end(), static_cast<double>(0.), [](auto &res, auto &val) { return res + val * val; });
89}
90
92template <typename IndexContainer>
93constexpr void reduceActiveAxis(const mesh::Mesh &mesh, const IndexContainer &IDs, std::array<bool, 3> &axis)
94{
95 // make a pair of the axis and the difference
97
98 // Compute the difference magnitude per direction
99 for (std::size_t d = 0; d < axis.size(); ++d) {
100 // Ignore dead axis here, i.e., apply the max value such that they are sorted on the last position(s)
101 if (axis[d] == false) {
102 differences[d] = std::make_pair<int, double>(d, std::numeric_limits<double>::max());
103 } else {
104 auto res = std::minmax_element(IDs.begin(), IDs.end(), [&](const auto &a, const auto &b) { return mesh.vertex(a).coord(d) < mesh.vertex(b).coord(d); });
105 // Check if we are above or below the threshold
106 differences[d] = std::make_pair<int, double>(d, std::abs(mesh.vertex(*res.second).coord(d) - mesh.vertex(*res.first).coord(d)));
107 }
108 }
109
110 std::sort(differences.begin(), differences.end(), [](const auto &d1, const auto &d2) { return d1.second < d2.second; });
111 // Disable the axis having the smallest expansion
112 axis[differences[0].first] = false;
113}
114
115// Fill in the polynomial entries
116template <typename IndexContainer>
117inline void fillPolynomialEntries(Eigen::MatrixXd &matrix, const mesh::Mesh &mesh, const IndexContainer &IDs, Eigen::Index startIndex, std::array<bool, 3> activeAxis)
118{
119 // Loop over all vertices in the mesh
120 for (const auto &i : IDs | boost::adaptors::indexed()) {
121
122 // 1. the constant contribution
123 matrix(i.index(), startIndex) = 1.0;
124
125 // 2. the linear contribution
126 const auto & u = mesh.vertex(i.value()).rawCoords();
127 unsigned int k = 0;
128 // Loop over all three space dimension and ignore dead axis
129 for (unsigned int d = 0; d < activeAxis.size(); ++d) {
130 if (activeAxis[d]) {
131 PRECICE_ASSERT(matrix.rows() > i.index(), matrix.rows(), i.index());
132 PRECICE_ASSERT(matrix.cols() > startIndex + 1 + k, matrix.cols(), startIndex + 1 + k);
133 matrix(i.index(), startIndex + 1 + k) = u[d];
134 ++k;
135 }
136 }
137 }
138}
139
140template <typename RADIAL_BASIS_FUNCTION_T, typename IndexContainer>
141Eigen::MatrixXd buildMatrixCLU(RADIAL_BASIS_FUNCTION_T basisFunction, const mesh::Mesh &inputMesh, const IndexContainer &inputIDs,
142 std::array<bool, 3> activeAxis, Polynomial polynomial)
143{
144 // Treat the 2D case as 3D case with dead axis
145 const unsigned int deadDimensions = std::count(activeAxis.begin(), activeAxis.end(), false);
146 const unsigned int dimensions = 3;
147 const unsigned int polyparams = polynomial == Polynomial::ON ? 1 + dimensions - deadDimensions : 0;
148
149 // Add linear polynom degrees if polynomial requires this
150 const auto inputSize = inputIDs.size();
151 const auto n = inputSize + polyparams;
152
153 PRECICE_ASSERT((inputMesh.getDimensions() == 3) || activeAxis[2] == false);
154 PRECICE_ASSERT((inputSize >= 1 + polyparams) || polynomial != Polynomial::ON, inputSize);
155
156 Eigen::MatrixXd matrixCLU(n, n);
157
158 // Required to fill the poly -> poly entries in the matrix, which remain otherwise untouched
159 if (polynomial == Polynomial::ON) {
160 matrixCLU.setZero();
161 }
162
163 // Compute RBF matrix entries
164 auto i_iter = inputIDs.begin();
165 Eigen::Index i_index = 0;
166 for (; i_iter != inputIDs.end(); ++i_iter, ++i_index) {
167 const auto &u = inputMesh.vertex(*i_iter).rawCoords();
168 auto j_iter = i_iter;
169 auto j_index = i_index;
170 for (; j_iter != inputIDs.end(); ++j_iter, ++j_index) {
171 const auto &v = inputMesh.vertex(*j_iter).rawCoords();
172 double squaredDifference = computeSquaredDifference(u, v, activeAxis);
173 matrixCLU(i_index, j_index) = basisFunction.evaluate(std::sqrt(squaredDifference));
174 }
175 }
176
177 // Add potentially the polynomial contribution in the matrix
178 if (polynomial == Polynomial::ON) {
179 fillPolynomialEntries(matrixCLU, inputMesh, inputIDs, inputSize, activeAxis);
180 }
181 matrixCLU.triangularView<Eigen::Lower>() = matrixCLU.transpose();
182 return matrixCLU;
183}
184
185template <typename RADIAL_BASIS_FUNCTION_T, typename IndexContainer>
186Eigen::MatrixXd buildMatrixA(RADIAL_BASIS_FUNCTION_T basisFunction, const mesh::Mesh &inputMesh, const IndexContainer &inputIDs,
187 const mesh::Mesh &outputMesh, const IndexContainer outputIDs, std::array<bool, 3> activeAxis, Polynomial polynomial)
188{
189 // Treat the 2D case as 3D case with dead axis
190 const unsigned int deadDimensions = std::count(activeAxis.begin(), activeAxis.end(), false);
191 const unsigned int dimensions = 3;
192 const unsigned int polyparams = polynomial == Polynomial::ON ? 1 + dimensions - deadDimensions : 0;
193
194 const auto inputSize = inputIDs.size();
195 const auto outputSize = outputIDs.size();
196 const auto n = inputSize + polyparams;
197
198 PRECICE_ASSERT((inputMesh.getDimensions() == 3) || activeAxis[2] == false);
199 PRECICE_ASSERT((inputSize >= 1 + polyparams) || polynomial != Polynomial::ON, inputSize);
200
201 Eigen::MatrixXd matrixA(outputSize, n);
202
203 // Compute RBF values for matrix A
204 for (const auto &i : outputIDs | boost::adaptors::indexed()) {
205 const auto &u = outputMesh.vertex(i.value()).rawCoords();
206 for (const auto &j : inputIDs | boost::adaptors::indexed()) {
207 const auto &v = inputMesh.vertex(j.value()).rawCoords();
208 double squaredDifference = computeSquaredDifference(u, v, activeAxis);
209 matrixA(i.index(), j.index()) = basisFunction.evaluate(std::sqrt(squaredDifference));
210 }
211 }
212
213 // Add potentially the polynomial contribution in the matrix
214 if (polynomial == Polynomial::ON) {
215 fillPolynomialEntries(matrixA, outputMesh, outputIDs, inputSize, activeAxis);
216 }
217 return matrixA;
218}
219
220template <typename RADIAL_BASIS_FUNCTION_T>
221template <typename IndexContainer>
222RadialBasisFctSolver<RADIAL_BASIS_FUNCTION_T>::RadialBasisFctSolver(RADIAL_BASIS_FUNCTION_T basisFunction, const mesh::Mesh &inputMesh, const IndexContainer &inputIDs,
223 const mesh::Mesh &outputMesh, const IndexContainer &outputIDs, std::vector<bool> deadAxis, Polynomial polynomial)
224{
225 PRECICE_ASSERT(!(RADIAL_BASIS_FUNCTION_T::isStrictlyPositiveDefinite() && polynomial == Polynomial::ON), "The integrated polynomial (polynomial=\"on\") is not supported for the selected radial-basis function. Please select another radial-basis function or change the polynomial configuration.");
226 // Convert dead axis vector into an active axis array so that we can handle the reduction more easily
227 std::array<bool, 3> activeAxis({{false, false, false}});
228 std::transform(deadAxis.begin(), deadAxis.end(), activeAxis.begin(), [](const auto ax) { return !ax; });
229
230 // First, assemble the interpolation matrix and check the invertability
231 bool decompositionSuccessful = false;
232 if constexpr (RADIAL_BASIS_FUNCTION_T::isStrictlyPositiveDefinite()) {
233 _decMatrixC = buildMatrixCLU(basisFunction, inputMesh, inputIDs, activeAxis, polynomial).llt();
234 decompositionSuccessful = _decMatrixC.info() == Eigen::ComputationInfo::Success;
235 } else {
236 _decMatrixC = buildMatrixCLU(basisFunction, inputMesh, inputIDs, activeAxis, polynomial).colPivHouseholderQr();
237 decompositionSuccessful = _decMatrixC.isInvertible();
238 }
239
240 PRECICE_CHECK(decompositionSuccessful,
241 "The interpolation matrix of the RBF mapping from mesh \"{}\" to mesh \"{}\" is not invertable. "
242 "This means that the mapping problem is not well-posed. "
243 "Please check if your coupling meshes are correct (e.g. no vertices are duplicated) or reconfigure "
244 "your basis-function (e.g. reduce the support-radius).",
245 inputMesh.getName(), outputMesh.getName());
246
247 // Second, assemble evaluation matrix
248 _matrixA = buildMatrixA(basisFunction, inputMesh, inputIDs, outputMesh, outputIDs, activeAxis, polynomial);
249
250 // In case we deal with separated polynomials, we need dedicated matrices for the polynomial contribution
251 if (polynomial == Polynomial::SEPARATE) {
252
253 // 4 = 1 + dimensions(3) = maximum number of polynomial parameters
254 auto localActiveAxis = activeAxis;
255 unsigned int polyParams = 4 - std::count(localActiveAxis.begin(), localActiveAxis.end(), false);
256
257 do {
258 // First, build matrix Q and check for the condition number
259 _matrixQ.resize(inputIDs.size(), polyParams);
260 fillPolynomialEntries(_matrixQ, inputMesh, inputIDs, 0, localActiveAxis);
261
262 // Compute the condition number
263 Eigen::JacobiSVD<Eigen::MatrixXd> svd(_matrixQ);
264 PRECICE_ASSERT(svd.singularValues().size() > 0);
265 PRECICE_DEBUG("Singular values in polynomial solver: {}", svd.singularValues());
266 const double conditionNumber = svd.singularValues()(0) / std::max(svd.singularValues()(svd.singularValues().size() - 1), math::NUMERICAL_ZERO_DIFFERENCE);
267 PRECICE_DEBUG("Condition number: {}", conditionNumber);
268
269 // Disable one axis
270 if (conditionNumber > 1e5) {
271 reduceActiveAxis(inputMesh, inputIDs, localActiveAxis);
272 polyParams = 4 - std::count(localActiveAxis.begin(), localActiveAxis.end(), false);
273 } else {
274 break;
275 }
276 } while (true);
277
278 // allocate and fill matrix V for the outputMesh
279 _matrixV.resize(outputIDs.size(), polyParams);
280 fillPolynomialEntries(_matrixV, outputMesh, outputIDs, 0, localActiveAxis);
281
282 // 3. compute decomposition
283 _qrMatrixQ = _matrixQ.colPivHouseholderQr();
284 }
285}
286
287template <typename RADIAL_BASIS_FUNCTION_T>
288Eigen::VectorXd RadialBasisFctSolver<RADIAL_BASIS_FUNCTION_T>::solveConservative(const Eigen::VectorXd &inputData, Polynomial polynomial) const
289{
290 PRECICE_ASSERT((_matrixV.size() > 0 && polynomial == Polynomial::SEPARATE) || _matrixV.size() == 0, _matrixV.size());
291 // TODO: Avoid temporary allocations
292 // Au is equal to the eta in our PETSc implementation
293 PRECICE_ASSERT(inputData.size() == _matrixA.rows());
294 Eigen::VectorXd Au = _matrixA.transpose() * inputData;
295 PRECICE_ASSERT(Au.size() == _matrixA.cols());
296
297 // mu in the PETSc implementation
298 Eigen::VectorXd out = _decMatrixC.solve(Au);
299
300 if (polynomial == Polynomial::SEPARATE) {
301 Eigen::VectorXd epsilon = _matrixV.transpose() * inputData;
302 PRECICE_ASSERT(epsilon.size() == _matrixV.cols());
303
304 // epsilon = Q^T * mu - epsilon (tau in the PETSc impl)
305 epsilon -= _matrixQ.transpose() * out;
306 PRECICE_ASSERT(epsilon.size() == _matrixQ.cols());
307
308 // out = out - solveTranspose tau (sigma in the PETSc impl)
309 // Newer version of eigen provide the solve() for transpose() matrix decopmositions
310#if EIGEN_VERSION_AT_LEAST(3, 4, 0)
311 out -= static_cast<Eigen::VectorXd>(_qrMatrixQ.transpose().solve(-epsilon));
312#else
313 // Backwards compatible version
314 Eigen::VectorXd sigma(_matrixQ.rows());
315 const Eigen::Index nonzero_pivots = _qrMatrixQ.nonzeroPivots();
316
317 if (nonzero_pivots == 0) {
318 sigma.setZero();
319 } else {
320 Eigen::VectorXd c(_qrMatrixQ.colsPermutation().transpose() * (-epsilon));
321
322 _qrMatrixQ.matrixQR().topLeftCorner(nonzero_pivots, nonzero_pivots).template triangularView<Eigen::Upper>().transpose().conjugate().solveInPlace(c.topRows(nonzero_pivots));
323
324 sigma.topRows(nonzero_pivots) = c.topRows(nonzero_pivots);
325 sigma.bottomRows(_qrMatrixQ.rows() - nonzero_pivots).setZero();
326
327 sigma.applyOnTheLeft(_qrMatrixQ.householderQ().setLength(nonzero_pivots));
328 out -= sigma;
329 }
330#endif
331 }
332 return out;
333}
334
335template <typename RADIAL_BASIS_FUNCTION_T>
336Eigen::VectorXd RadialBasisFctSolver<RADIAL_BASIS_FUNCTION_T>::solveConsistent(Eigen::VectorXd &inputData, Polynomial polynomial) const
337{
338 PRECICE_ASSERT((_matrixQ.size() > 0 && polynomial == Polynomial::SEPARATE) || _matrixQ.size() == 0);
339 Eigen::VectorXd polynomialContribution;
340 // Solve polynomial QR and subtract it from the input data
341 if (polynomial == Polynomial::SEPARATE) {
342 polynomialContribution = _qrMatrixQ.solve(inputData);
343 inputData -= (_matrixQ * polynomialContribution);
344 }
345
346 // Integrated polynomial (and separated)
347 PRECICE_ASSERT(inputData.size() == _matrixA.cols());
348 Eigen::VectorXd p = _decMatrixC.solve(inputData);
349 PRECICE_ASSERT(p.size() == _matrixA.cols());
350 Eigen::VectorXd out = _matrixA * p;
351
352 // Add the polynomial part again for separated polynomial
353 if (polynomial == Polynomial::SEPARATE) {
354 out += (_matrixV * polynomialContribution);
355 }
356 return out;
357}
358
359template <typename RADIAL_BASIS_FUNCTION_T>
361{
362 _matrixA = Eigen::MatrixXd();
363 _decMatrixC = DecompositionType();
364}
365
366template <typename RADIAL_BASIS_FUNCTION_T>
368{
369 return _matrixA.cols();
370}
371
372template <typename RADIAL_BASIS_FUNCTION_T>
374{
375 return _matrixA.rows();
376}
377} // namespace mapping
378} // namespace precice
std::ostream & out
#define PRECICE_DEBUG(...)
Definition LogMacros.hpp:64
#define PRECICE_CHECK(check,...)
Definition LogMacros.hpp:35
T accumulate(T... args)
#define PRECICE_ASSERT(...)
Definition assertion.hpp:87
T begin(T... args)
This class provides a lightweight logger.
Definition Logger.hpp:16
Eigen::MatrixXd _matrixA
Evaluation matrix (output x input)
Eigen::VectorXd solveConsistent(Eigen::VectorXd &inputData, Polynomial polynomial) const
Maps the given input data.
Eigen::VectorXd solveConservative(const Eigen::VectorXd &inputData, Polynomial polynomial) const
Maps the given input data.
Eigen::ColPivHouseholderQR< Eigen::MatrixXd > _qrMatrixQ
Decomposition of the polynomial (for separate polynomial)
RadialBasisFctSolver()=default
Default constructor.
Eigen::MatrixXd _matrixV
Polynomial matrix of the output mesh (for separate polynomial)
DecompositionType _decMatrixC
Decomposition of the interpolation matrix.
Eigen::MatrixXd _matrixQ
Polynomial matrix of the input mesh (for separate polynomial)
Container and creator for meshes.
Definition Mesh.hpp:39
int getDimensions() const
Definition Mesh.cpp:98
const std::string & getName() const
Returns the name of the mesh, as set in the config file.
Definition Mesh.cpp:218
Vertex & vertex(VertexID id)
Mutable access to a vertex by VertexID.
Definition Mesh.cpp:41
double coord(int index) const
Returns a coordinate of a vertex.
Definition Vertex.hpp:128
const RawCoords & rawCoords() const
Direct access to the coordinates.
Definition Vertex.hpp:123
T count(T... args)
T end(T... args)
T max(T... args)
T minmax_element(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.
void fillPolynomialEntries(Eigen::MatrixXd &matrix, const mesh::Mesh &mesh, const IndexContainer &IDs, Eigen::Index startIndex, std::array< bool, 3 > activeAxis)
constexpr void reduceActiveAxis(const mesh::Mesh &mesh, const IndexContainer &IDs, std::array< bool, 3 > &axis)
given the active axis, computes sets the axis with the lowest spatial expansion to dead
Eigen::MatrixXd buildMatrixCLU(RADIAL_BASIS_FUNCTION_T basisFunction, const mesh::Mesh &inputMesh, const IndexContainer &inputIDs, std::array< bool, 3 > activeAxis, Polynomial polynomial)
Eigen::MatrixXd buildMatrixA(RADIAL_BASIS_FUNCTION_T basisFunction, const mesh::Mesh &inputMesh, const IndexContainer &inputIDs, const mesh::Mesh &outputMesh, const IndexContainer outputIDs, std::array< bool, 3 > activeAxis, Polynomial polynomial)
Polynomial
How to handle the polynomial?
constexpr double NUMERICAL_ZERO_DIFFERENCE
Main namespace of the precice library.
T size(T... args)
T sort(T... args)
T sqrt(T... args)
T transform(T... args)