Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2025 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 :
4 : // This library is free software; you can redistribute it and/or
5 : // modify it under the terms of the GNU Lesser General Public
6 : // License as published by the Free Software Foundation; either
7 : // version 2.1 of the License, or (at your option) any later version.
8 :
9 : // This library is distributed in the hope that it will be useful,
10 : // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 : // Lesser General Public License for more details.
13 :
14 : // You should have received a copy of the GNU Lesser General Public
15 : // License along with this library; if not, write to the Free Software
16 : // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 :
18 :
19 :
20 : // Local includes
21 : #include "libmesh/radial_basis_interpolation.h"
22 :
23 : #include "libmesh/eigen_core_support.h"
24 : #include "libmesh/int_range.h"
25 : #include "libmesh/libmesh_logging.h"
26 : #include "libmesh/mesh_tools.h" // BoundingBox
27 : #include "libmesh/radial_basis_functions.h"
28 :
29 : #ifdef LIBMESH_HAVE_EIGEN
30 : # include "libmesh/ignore_warnings.h"
31 : # include <Eigen/Dense>
32 : # include "libmesh/restore_warnings.h"
33 : #endif
34 :
35 : // C++ includes
36 : #include <iomanip>
37 :
38 :
39 : namespace libMesh
40 : {
41 : //--------------------------------------------------------------------------------
42 : // RadialBasisInterpolation methods
43 : template <unsigned int KDDim, class RBF>
44 0 : void RadialBasisInterpolation<KDDim,RBF>::clear()
45 : {
46 : // Call base class clear method
47 0 : InverseDistanceInterpolation<KDDim>::clear();
48 0 : }
49 :
50 :
51 :
52 : template <unsigned int KDDim, class RBF>
53 142 : void RadialBasisInterpolation<KDDim,RBF>::prepare_for_use()
54 : {
55 : // Call base class methods for prep
56 142 : InverseDistanceInterpolation<KDDim>::prepare_for_use();
57 142 : InverseDistanceInterpolation<KDDim>::construct_kd_tree();
58 :
59 : #ifndef LIBMESH_HAVE_EIGEN
60 :
61 : libmesh_error_msg("ERROR: this functionality presently requires Eigen!");
62 :
63 : #else
64 8 : LOG_SCOPE ("prepare_for_use()", "RadialBasisInterpolation<>");
65 :
66 : // Construct a bounding box for our source points
67 4 : _src_bbox.invalidate();
68 :
69 146 : const std::size_t n_src_pts = this->_src_pts.size();
70 142 : const unsigned int n_vars = this->n_field_variables();
71 4 : libmesh_assert_equal_to (this->_src_vals.size(), n_src_pts*this->n_field_variables());
72 :
73 : {
74 8 : LOG_SCOPE ("prepare_for_use():bbox", "RadialBasisInterpolation<>");
75 :
76 : Point
77 4 : &p_min(_src_bbox.min()),
78 4 : &p_max(_src_bbox.max());
79 :
80 238063 : for (std::size_t p=0; p<n_src_pts; p++)
81 : {
82 13404 : const Point & p_src(_src_pts[p]);
83 :
84 951684 : for (unsigned int d=0; d<LIBMESH_DIM; d++)
85 : {
86 713763 : p_min(d) = std::min(p_min(d), p_src(d));
87 713763 : p_max(d) = std::max(p_max(d), p_src(d));
88 : }
89 : }
90 : }
91 :
92 : // Debugging code
93 : // libMesh::out << "bounding box is \n"
94 : // << _src_bbox.min() << '\n'
95 : // << _src_bbox.max() << std::endl;
96 :
97 :
98 : // Construct the Radial Basis Function, giving it the size of the domain
99 142 : if (_r_override < 0)
100 142 : _r_bbox = (_src_bbox.max() - _src_bbox.min()).norm();
101 : else
102 0 : _r_bbox = _r_override;
103 :
104 142 : RBF rbf(_r_bbox);
105 :
106 : // libMesh::out << "bounding box is \n"
107 : // << _src_bbox.min() << '\n'
108 : // << _src_bbox.max() << '\n'
109 : // << "r_bbox = " << _r_bbox << '\n'
110 : // << "rbf(r_bbox/2) = " << rbf(_r_bbox/2) << std::endl;
111 :
112 :
113 : // Construct the projection Matrix
114 : typedef Eigen::Matrix<Number, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor> DynamicMatrix;
115 : //typedef Eigen::Matrix<Number, Eigen::Dynamic, 1, Eigen::ColMajor> DynamicVector;
116 :
117 146 : DynamicMatrix A(n_src_pts, n_src_pts), x(n_src_pts,n_vars), b(n_src_pts,n_vars);
118 :
119 : {
120 8 : LOG_SCOPE ("prepare_for_use():mat", "RadialBasisInterpolation<>");
121 :
122 238063 : for (std::size_t i=0; i<n_src_pts; i++)
123 : {
124 13404 : const Point & x_i (_src_pts[i]);
125 :
126 : // Diagonal
127 237921 : A(i,i) = rbf(0.);
128 :
129 375673496 : for (std::size_t j=i+1; j<n_src_pts; j++)
130 : {
131 21151300 : const Point & x_j (_src_pts[j]);
132 :
133 364859925 : const Real r_ij = (x_j - x_i).norm();
134 :
135 396586875 : A(i,j) = A(j,i) = rbf(r_ij);
136 : }
137 :
138 : // set source data
139 482942 : for (unsigned int var=0; var<n_vars; var++)
140 258825 : b(i,var) = _src_vals[i*n_vars + var];
141 : }
142 : }
143 :
144 :
145 : {
146 4 : LOG_SCOPE ("prepare_for_use():solve", "RadialBasisInterpolation<>");
147 :
148 : // Solve the linear system
149 142 : x = A.ldlt().solve(b);
150 : //x = A.fullPivLu().solve(b);
151 : }
152 :
153 : // save the weights for each variable
154 146 : _weights.resize (this->_src_vals.size());
155 :
156 238063 : for (std::size_t i=0; i<n_src_pts; i++)
157 482942 : for (unsigned int var=0; var<n_vars; var++)
158 258825 : _weights[i*n_vars + var] = x(i,var);
159 :
160 : #endif
161 :
162 142 : }
163 :
164 :
165 :
166 : template <unsigned int KDDim, class RBF>
167 36922 : void RadialBasisInterpolation<KDDim,RBF>::interpolate_field_data (const std::vector<std::string> & field_names,
168 : const std::vector<Point> & tgt_pts,
169 : std::vector<Number> & tgt_vals) const
170 : {
171 5900 : LOG_SCOPE ("interpolate_field_data()", "RadialBasisInterpolation<>");
172 :
173 : libmesh_experimental();
174 :
175 : const unsigned int
176 2950 : n_vars = this->n_field_variables();
177 :
178 : const std::size_t
179 5900 : n_src_pts = this->_src_pts.size(),
180 5900 : n_tgt_pts = tgt_pts.size();
181 :
182 2950 : libmesh_assert_equal_to (_weights.size(), this->_src_vals.size());
183 2950 : libmesh_assert_equal_to (field_names.size(), this->n_field_variables());
184 :
185 : // If we already have field variables, we assume we are appending.
186 : // that means the names and ordering better be identical!
187 39872 : libmesh_error_msg_if(this->_names.size() != field_names.size(),
188 : "ERROR: when adding field data to an existing list the \nvariable list must be the same!");
189 :
190 73915 : for (auto v : index_range(this->_names))
191 36993 : libmesh_error_msg_if(_names[v] != field_names[v],
192 : "ERROR: when adding field data to an existing list the \nvariable list must be the same!");
193 :
194 :
195 36922 : RBF rbf(_r_bbox);
196 :
197 36922 : tgt_vals.resize (n_tgt_pts*n_vars); /**/ std::fill (tgt_vals.begin(), tgt_vals.end(), Number(0.));
198 :
199 80873 : for (std::size_t tgt=0; tgt<n_tgt_pts; tgt++)
200 : {
201 6296 : const Point & p (tgt_pts[tgt]);
202 :
203 120556552 : for (std::size_t i=0; i<n_src_pts; i++)
204 : {
205 19207896 : const Point & x_i(_src_pts[i]);
206 : const Real
207 110908653 : r_i = (p - x_i).norm(),
208 9603948 : phi_i = rbf(r_i);
209 :
210 241735202 : for (unsigned int var=0; var<n_vars; var++)
211 140470497 : tgt_vals[tgt*n_vars + var] += _weights[i*n_vars + var]*phi_i;
212 : }
213 : }
214 36922 : }
215 :
216 :
217 :
218 : // ------------------------------------------------------------
219 : // Explicit Instantiations
220 : template class LIBMESH_EXPORT RadialBasisInterpolation<3, WendlandRBF<3,0>>;
221 : template class LIBMESH_EXPORT RadialBasisInterpolation<3, WendlandRBF<3,2>>;
222 : template class LIBMESH_EXPORT RadialBasisInterpolation<3, WendlandRBF<3,4>>;
223 : template class LIBMESH_EXPORT RadialBasisInterpolation<3, WendlandRBF<3,8>>;
224 :
225 : } // namespace libMesh
|
