libMesh
radial_basis_interpolation.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2017 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 // C++ includes
21 #include <iomanip>
22 
23 // Local includes
24 #include "libmesh/radial_basis_interpolation.h"
25 #include "libmesh/radial_basis_functions.h"
26 #include "libmesh/mesh_tools.h" // BoundingBox
27 #include "libmesh/libmesh_logging.h"
28 #include "libmesh/eigen_core_support.h"
29 
30 #ifdef LIBMESH_HAVE_EIGEN
31 #include <Eigen/Dense>
32 #endif
33 
34 
35 namespace libMesh
36 {
37 //--------------------------------------------------------------------------------
38 // RadialBasisInterpolation methods
39 template <unsigned int KDDim, class RBF>
41 {
42  // Call base class clear method
44 }
45 
46 
47 
48 template <unsigned int KDDim, class RBF>
50 {
51  // Call base class methods for prep
54 
55 #ifndef LIBMESH_HAVE_EIGEN
56 
57  libmesh_error_msg("ERROR: this functionality presently requires Eigen!");
58 
59 #else
60  LOG_SCOPE ("prepare_for_use()", "RadialBasisInterpolation<>");
61 
62  // Construct a bounding box for our source points
63  _src_bbox.invalidate();
64 
65  const std::size_t n_src_pts = this->_src_pts.size();
66  const unsigned int n_vars = this->n_field_variables();
67  libmesh_assert_equal_to (this->_src_vals.size(), n_src_pts*this->n_field_variables());
68 
69  {
70  Point
71  &p_min(_src_bbox.min()),
72  &p_max(_src_bbox.max());
73 
74  for (std::size_t p=0; p<n_src_pts; p++)
75  {
76  const Point & p_src(_src_pts[p]);
77 
78  for (unsigned int d=0; d<LIBMESH_DIM; d++)
79  {
80  p_min(d) = std::min(p_min(d), p_src(d));
81  p_max(d) = std::max(p_max(d), p_src(d));
82  }
83  }
84  }
85 
86  libMesh::out << "bounding box is \n"
87  << _src_bbox.min() << '\n'
88  << _src_bbox.max() << std::endl;
89 
90 
91  // Construct the Radial Basis Function, giving it the size of the domain
92  if (_r_override < 0)
93  _r_bbox = (_src_bbox.max() - _src_bbox.min()).norm();
94  else
95  _r_bbox = _r_override;
96 
97  RBF rbf(_r_bbox);
98 
99  libMesh::out << "bounding box is \n"
100  << _src_bbox.min() << '\n'
101  << _src_bbox.max() << '\n'
102  << "r_bbox = " << _r_bbox << '\n'
103  << "rbf(r_bbox/2) = " << rbf(_r_bbox/2) << std::endl;
104 
105 
106  // Construct the projection Matrix
107  typedef Eigen::Matrix<Number, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor> DynamicMatrix;
108  //typedef Eigen::Matrix<Number, Eigen::Dynamic, 1, Eigen::ColMajor> DynamicVector;
109 
110  DynamicMatrix A(n_src_pts, n_src_pts), x(n_src_pts,n_vars), b(n_src_pts,n_vars);
111 
112  for (std::size_t i=0; i<n_src_pts; i++)
113  {
114  const Point & x_i (_src_pts[i]);
115 
116  // Diagonal
117  A(i,i) = rbf(0.);
118 
119  for (std::size_t j=i+1; j<n_src_pts; j++)
120  {
121  const Point & x_j (_src_pts[j]);
122 
123  const Real r_ij = (x_j - x_i).norm();
124 
125  A(i,j) = A(j,i) = rbf(r_ij);
126  }
127 
128  // set source data
129  for (unsigned int var=0; var<n_vars; var++)
130  b(i,var) = _src_vals[i*n_vars + var];
131  }
132 
133 
134  // Solve the linear system
135  x = A.ldlt().solve(b);
136  //x = A.fullPivLu().solve(b);
137 
138  // save the weights for each variable
139  _weights.resize (this->_src_vals.size());
140 
141  for (std::size_t i=0; i<n_src_pts; i++)
142  for (unsigned int var=0; var<n_vars; var++)
143  _weights[i*n_vars + var] = x(i,var);
144 
145 #endif
146 
147 }
148 
149 
150 
151 template <unsigned int KDDim, class RBF>
152 void RadialBasisInterpolation<KDDim,RBF>::interpolate_field_data (const std::vector<std::string> & field_names,
153  const std::vector<Point> & tgt_pts,
154  std::vector<Number> & tgt_vals) const
155 {
156  LOG_SCOPE ("interpolate_field_data()", "RadialBasisInterpolation<>");
157 
158  libmesh_experimental();
159 
160  const unsigned int
161  n_vars = this->n_field_variables();
162 
163  const std::size_t
164  n_src_pts = this->_src_pts.size(),
165  n_tgt_pts = tgt_pts.size();
166 
167  libmesh_assert_equal_to (_weights.size(), this->_src_vals.size());
168  libmesh_assert_equal_to (field_names.size(), this->n_field_variables());
169 
170  // If we already have field variables, we assume we are appending.
171  // that means the names and ordering better be identical!
172  if (this->_names.size() != field_names.size())
173  libmesh_error_msg("ERROR: when adding field data to an existing list the \nvariable list must be the same!");
174 
175  for (std::size_t v=0; v<this->_names.size(); v++)
176  if (_names[v] != field_names[v])
177  libmesh_error_msg("ERROR: when adding field data to an existing list the \nvariable list must be the same!");
178 
179 
180  RBF rbf(_r_bbox);
181 
182  tgt_vals.resize (n_tgt_pts*n_vars); std::fill (tgt_vals.begin(), tgt_vals.end(), Number(0.));
183 
184  for (std::size_t tgt=0; tgt<n_tgt_pts; tgt++)
185  {
186  const Point & p (tgt_pts[tgt]);
187 
188  for (std::size_t i=0; i<n_src_pts; i++)
189  {
190  const Point & x_i(_src_pts[i]);
191  const Real
192  r_i = (p - x_i).norm(),
193  phi_i = rbf(r_i);
194 
195  for (unsigned int var=0; var<n_vars; var++)
196  tgt_vals[tgt*n_vars + var] += _weights[i*n_vars + var]*phi_i;
197  }
198  }
199 }
200 
201 
202 
203 // ------------------------------------------------------------
204 // Explicit Instantiations
209 
210 } // namespace libMesh
Radial Basis Function interpolation.
virtual void interpolate_field_data(const std::vector< std::string > &field_names, const std::vector< Point > &tgt_pts, std::vector< Number > &tgt_vals) const libmesh_override
Interpolate source data at target points.
virtual void prepare_for_use() libmesh_override
Prepares data structures for use.
const unsigned int n_vars
Definition: tecplot_io.C:68
The libMesh namespace provides an interface to certain functionality in the library.
long double max(long double a, double b)
PetscErrorCode Vec x
virtual void construct_kd_tree()
Build & initialize the KD tree, if needed.
virtual void prepare_for_use()
Prepares data structures for use.
virtual void clear() libmesh_override
Clears all internal data structures and restores to a pristine state.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
static PetscErrorCode Mat * A
OStreamProxy out
virtual void clear() libmesh_override
Clears all internal data structures and restores to a pristine state.
long double min(long double a, double b)
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:38