| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | // This file is part of Bembel, the higher order C++ boundary element library. | ||
| 2 | // | ||
| 3 | // Copyright (C) 2022 see <http://www.bembel.eu> | ||
| 4 | // | ||
| 5 | // It was written as part of a cooperation of J. Doelz, H. Harbrecht, S. Kurz, | ||
| 6 | // M. Multerer, S. Schoeps, and F. Wolf at Technische Universitaet Darmstadt, | ||
| 7 | // Universitaet Basel, and Universita della Svizzera italiana, Lugano. This | ||
| 8 | // source code is subject to the GNU General Public License version 3 and | ||
| 9 | // provided WITHOUT ANY WARRANTY, see <http://www.bembel.eu> for further | ||
| 10 | // information. | ||
| 11 | #ifndef BEMBEL_SRC_SPLINE_LOCALIZE_HPP_ | ||
| 12 | #define BEMBEL_SRC_SPLINE_LOCALIZE_HPP_ | ||
| 13 | |||
| 14 | namespace Bembel { | ||
| 15 | |||
| 16 | /** | ||
| 17 | * \ingroup Spline | ||
| 18 | * \brief The Spl namespace is the backend for basis functions and B-Spline | ||
| 19 | * computations, and should not be used by the user in a boundary element code. | ||
| 20 | */ | ||
| 21 | namespace Spl { | ||
| 22 | // At first the a==0 and b==1 seems insane. However, this is the case quite | ||
| 23 | // often and actually shaves off a couple of % of runtime, since the comparison | ||
| 24 | // is cheaper than a division. | ||
| 25 | 199812216 | constexpr inline double Rescale(double x, double a, double b) noexcept { | |
| 26 |
4/4✓ Branch 0 taken 199812088 times.
✓ Branch 1 taken 128 times.
✓ Branch 2 taken 128 times.
✓ Branch 3 taken 199811960 times.
|
199812216 | return (a == 0 && b == 1) ? x : (x - a) / (b - a); |
| 27 | } | ||
| 28 | |||
| 29 | // We use equidistant points for our interpolation problems. This is not | ||
| 30 | // optimal, but works sufficiently well. | ||
| 31 | 261 | inline std::vector<double> MakeInterpolationMask( | |
| 32 | int polynomial_degree) noexcept { | ||
| 33 | 261 | std::vector<double> out(polynomial_degree); | |
| 34 | 261 | const double h = 1. / (polynomial_degree + 1); | |
| 35 |
2/2✓ Branch 0 taken 872 times.
✓ Branch 1 taken 261 times.
|
1133 | for (int i = 0; i < polynomial_degree; i++) { |
| 36 | 872 | out[i] = ((i + 1) * h); | |
| 37 | } | ||
| 38 | 261 | return out; | |
| 39 | } | ||
| 40 | |||
| 41 | /** | ||
| 42 | * \ingroup Spline | ||
| 43 | * \brief Creates a T-vector of equidistant itnerpolation points for a given | ||
| 44 | * knot vector without any repetitions. | ||
| 45 | **/ | ||
| 46 | template <typename T> | ||
| 47 | 2 | inline std::vector<T> MakeInterpolationPoints( | |
| 48 | const std::vector<T> &uniq, const std::vector<T> &mask) noexcept { | ||
| 49 | 2 | const int size = uniq.size(); | |
| 50 | 2 | std::vector<T> out; | |
| 51 | |||
| 52 | 2 | out.reserve((size - 1) * mask.size()); | |
| 53 | |||
| 54 |
2/2✓ Branch 0 taken 4 times.
✓ Branch 1 taken 2 times.
|
6 | for (int i = 0; i < size - 1; i++) { |
| 55 |
2/2✓ Branch 4 taken 8 times.
✓ Branch 5 taken 4 times.
|
12 | for (auto m : mask) { |
| 56 | 8 | out.push_back(uniq[i] + m * (uniq[i + 1] - uniq[i])); | |
| 57 | } | ||
| 58 | } | ||
| 59 | 2 | return out; | |
| 60 | } | ||
| 61 | |||
| 62 | /** | ||
| 63 | * \ingroup Spline | ||
| 64 | * \brief returns the coefficients to represent a function in the Bernstein | ||
| 65 | * basis on [0,1]. | ||
| 66 | **/ | ||
| 67 | inline Eigen::Matrix<double, -1, -1> GetInterpolationMatrix( | ||
| 68 | int polynomial_degree, const std::vector<double> &mask) { | ||
| 69 | Eigen::Matrix<double, -1, -1> interpolationMatrix(polynomial_degree + 1, | ||
| 70 | polynomial_degree + 1); | ||
| 71 | |||
| 72 | double val[Constants::MaxP + 1]; | ||
| 73 | for (int j = 0; j < polynomial_degree + 1; j++) { | ||
| 74 | Bembel::Basis::ShapeFunctionHandler::evalBasis(polynomial_degree, val, | ||
| 75 | mask[j]); | ||
| 76 | for (int i = 0; i < polynomial_degree + 1; i++) | ||
| 77 | interpolationMatrix(j, i) = val[i]; | ||
| 78 | } | ||
| 79 | |||
| 80 | return interpolationMatrix.inverse(); | ||
| 81 | } | ||
| 82 | |||
| 83 | /** | ||
| 84 | * \ingroup Spline | ||
| 85 | * \brief this solves a generic interpolation problem. | ||
| 86 | **/ | ||
| 87 | template <typename T> | ||
| 88 | inline void GetCoefficients(const int incr, const std::vector<double> &mask, | ||
| 89 | const std::vector<T> &values, T *coefs) { | ||
| 90 | int polynomial_degree = mask.size() - 1; | ||
| 91 | |||
| 92 | Eigen::Matrix<double, -1, -1> im = | ||
| 93 | GetInterpolationMatrix(polynomial_degree, mask); | ||
| 94 | |||
| 95 | Eigen::Matrix<T, -1, 1> rhs(polynomial_degree + 1); | ||
| 96 | |||
| 97 | for (int i = 0; i < incr; i++) { | ||
| 98 | for (int j = 0; j < polynomial_degree + 1; j++) { | ||
| 99 | rhs[j] = values[i * (polynomial_degree + 1) + j]; | ||
| 100 | } | ||
| 101 | Eigen::Matrix<T, -1, 1> tmp(polynomial_degree + 1); | ||
| 102 | tmp = im * rhs; | ||
| 103 | |||
| 104 | for (int j = 0; j < polynomial_degree + 1; j++) { | ||
| 105 | coefs[i * (polynomial_degree + 1) + j] = tmp(j); | ||
| 106 | } | ||
| 107 | } | ||
| 108 | |||
| 109 | return; | ||
| 110 | } | ||
| 111 | |||
| 112 | /** | ||
| 113 | * \ingroup Spline | ||
| 114 | * \brief this solves a generic interpolation problem. | ||
| 115 | **/ | ||
| 116 | template <typename T> | ||
| 117 | inline std::vector<T> GetCoefficients(const int increments, | ||
| 118 | const std::vector<double> &mask, | ||
| 119 | const std::vector<T> &values) { | ||
| 120 | std::vector<double> out(increments * (mask.size())); | ||
| 121 | GetCoefficients<T>(increments, mask, values, out.data()); | ||
| 122 | return out; | ||
| 123 | } | ||
| 124 | } // namespace Spl | ||
| 125 | } // namespace Bembel | ||
| 126 | #endif // BEMBEL_SRC_SPLINE_LOCALIZE_HPP_ | ||
| 127 |