Rheolef  7.1
an efficient C++ finite element environment
interpolate.h
Go to the documentation of this file.
1 #ifndef _RHEOLEF_INTERPOLATE_H
2 #define _RHEOLEF_INTERPOLATE_H
3 //
4 // This file is part of Rheolef.
5 //
6 // Copyright (C) 2000-2009 Pierre Saramito <Pierre.Saramito@imag.fr>
7 //
8 // Rheolef is free software; you can redistribute it and/or modify
9 // it under the terms of the GNU General Public License as published by
10 // the Free Software Foundation; either version 2 of the License, or
11 // (at your option) any later version.
12 //
13 // Rheolef is distributed in the hope that it will be useful,
14 // but WITHOUT ANY WARRANTY; without even the implied warranty of
15 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 // GNU General Public License for more details.
17 //
18 // You should have received a copy of the GNU General Public License
19 // along with Rheolef; if not, write to the Free Software
20 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
21 //
22 // =========================================================================
23 // author: Pierre.Saramito@imag.fr
24 // date: 15 september 2015
25 
26 namespace rheolef {
90 } // namespace rheolef
91 
92 //
93 // SUMMARY:
94 // 1. implementation
95 // 1.1. scalar-valued result field
96 // 1.2. vector-valued case
97 // 1.3. tensor-valued case
98 // 1.4. undeterminated-valued case: determined at run-time
99 // 1.5. interface of the internal interpolate() function
100 // 2. the interpolate() function
101 // 2.1. nonlinear expression and function/functor
102 // 2.2. re-interpolation of fields and linear field expressions
103 
104 #include "rheolef/field_expr.h"
105 namespace rheolef {
106 
107 namespace details {
108 // --------------------------------------------------------------------------
109 // 1. implementation: general nonlinear expr
110 // --------------------------------------------------------------------------
111 // notes:
112 // * function template partial specialization is not allowed --> use class-function
113 // * interpolation on boundary domain spaces (geo_domain) could generate
114 // some communications (see interpolate_dom[234]_tst.cc for tests)
115 // * interpolation of functions and functor are be performed in all cases
116 // without communications, see below for this implementation.
117 // here we suppose a general nonlinear expression that is interpolated
118 // by using a loop on mesh elements
119 //
120 template<class T, class M, class Expr, class Result>
122  const space_basic<T,M>& Xh,
123  const Expr& expr)
124 {
125 trace_macro ("Expr="<<pretty_typename_macro(Expr));
126  typedef typename field_basic<T,M>::size_type size_type;
128  "interpolate: incompatible "<<Xh.valued()<<"-valued space and "
130  std::vector<size_type> dis_idof;
131  Eigen::Matrix<Result,Eigen::Dynamic,1> value;
132  Eigen::Matrix<T,Eigen::Dynamic,1> udof;
133  field_basic<T,M> uh (Xh, std::numeric_limits<T>::max());
134  const geo_basic<T,M>& omega = Xh.get_geo();
135  const basis_basic<T>& b = Xh.get_basis();
136  const piola_fem<T>& pf = b.get_piola_fem();
137 trace_macro ("pf.transform_need_piola="<<pf.transform_need_piola());
138  integrate_option iopt;
140  pops.initialize (omega.get_piola_basis(), b, iopt);
141  expr.initialize (Xh, pops, iopt);
142  expr.template valued_check<Result>();
143  for (typename geo_basic<T,M>::const_iterator
144  iter_ie = omega.begin(),
145  last_ie = omega.end(); iter_ie != last_ie; ++iter_ie) {
146  const geo_element& K = *iter_ie;
147  reference_element hat_K = K;
148  // 1) get u values at nodes of K
149  expr.evaluate (omega, K, value);
150  // 2) u values at nodes of K -> udofs on K
151  if (pf.transform_need_piola()) {
152  const Eigen::Matrix<piola<T>,Eigen::Dynamic,1>& piola = pops.get_piola (omega, K);
153  for (size_type loc_inod = 0, loc_nnod = value.size(); loc_inod < loc_nnod; ++loc_inod) {
154  // be carefull: piola_fem::inv_transform should support inplace call in the "value" arg
155 #ifndef TO_CLEAN
156  Result old_value = value[loc_inod];
157 #endif // TO_CLEAN
158  pf.inv_transform (piola[loc_inod], value[loc_inod], value[loc_inod]);
159 #ifndef TO_CLEAN
160  trace_macro ("inv_transf(K="<<K.name()<<K.dis_ie()<<",loc_inod="<<loc_inod<<"): old_value="
161  <<old_value<<", value="<<value[loc_inod]);
162 #endif // TO_CLEAN
163  }
164  }
165  b.compute_dofs (hat_K, value, udof);
166  // 3) copy all local dofs into the global field
167  Xh.dis_idof (K, dis_idof);
168  check_macro (b.ndof(hat_K) == dis_idof.size() && b.ndof(hat_K) == size_type(udof.size()),
169  "invalid sizes: basis("<<b.name()<<").size("<<hat_K.name()<<") = "<< b.ndof(hat_K)
170  <<", dis_idof.size="<<dis_idof.size()<<", udof.size="<<udof.size());
171  for (size_type loc_idof = 0, loc_ndof = udof.size(); loc_idof < loc_ndof; ++loc_idof) {
172  uh.dis_dof_entry (dis_idof[loc_idof]) = udof [loc_idof];
173 #ifndef TO_CLEAN
174  trace_macro ("uh(K="<<K.name()<<K.dis_ie()<<",loc_idof="<<loc_idof<<") = uh(dis_idof="<<dis_idof[loc_idof]
175  << ") = " << udof [loc_idof]);
176 #endif // TO_CLEAN
177  }
178  }
179  uh.dis_dof_update();
180 trace_macro ("interpolate done");
181  return uh;
182 }
183 template<class T, class M, class Expr, class Result, class Status = typename details::is_equal<Result,typename Expr::value_type>::type>
187  const space_basic<T,M>& Xh,
188  const Expr& expr) const
189 {
190  trace_macro ("Expr="<<pretty_typename_macro(Expr));
191  trace_macro ("Result="<<typename_macro(Result));
192  trace_macro ("Status="<<typename_macro(Status));
193  trace_macro ("Expr::value_type="<<typename_macro(typename Expr::value_type));
194  fatal_macro ("invalid type resolution");
195  return field_basic<T,M>();
196 }};
197 // 1.1. scalar-valued result field
198 template<class T, class M, class Expr>
199 struct interpolate_internal_check<T,M,Expr,T,std::true_type> {
202  const space_basic<T,M>& Xh,
203  const Expr& expr) const
204 {
205  return interpolate_generic<T,M,Expr,T>(Xh,expr);
206 }};
207 // 1.2. vector-valued case
208 template<class T, class M, class Expr>
209 struct interpolate_internal_check<T,M,Expr,point_basic<T>,std::true_type> {
212  const space_basic<T,M>& Xh,
213  const Expr& expr) const
214 {
215  return interpolate_generic<T,M,Expr,point_basic<T>>(Xh,expr);
216 }};
217 // 1.3. tensor-valued case
218 template<class T, class M, class Expr>
219 struct interpolate_internal_check<T,M,Expr,tensor_basic<T>,std::true_type> {
222  const space_basic<T,M>& Xh,
223  const Expr& expr) const
224 {
225  return interpolate_generic<T,M,Expr,tensor_basic<T>>(Xh,expr);
226 }};
227 // 1.4. undeterminated-valued case: determined at run-time
228 template<class T, class M, class Expr, class Status>
232  const space_basic<T,M>& Xh,
233  const Expr& expr) const
234 {
235  switch (expr.valued_tag()) {
236  case space_constant::scalar: {
238  return eval (Xh, expr);
239  }
240  case space_constant::vector: {
242  return eval (Xh, expr);
243  }
247  return eval (Xh, expr);
248  }
249  default:
250  warning_macro ("Expr="<<pretty_typename_macro(Expr));
251  warning_macro ("Status="<<typename_macro(Status));
252  error_macro ("unexpected `"
253  << space_constant::valued_name (expr.valued_tag())
254  << "' valued expression");
255  return field_basic<T,M>();
256  }
257 }};
258 // 1.5. interface of the internal interpolate() function
259 template<class T, class M, class Expr, class Result>
262  const space_basic<T,M>& Xh,
263  const Expr& expr)
264 {
266  return eval (Xh,expr);
267 }
268 
269 } // namespace details
270 // --------------------------------------------------------------------------
271 // 2. the interpolate() function
272 // --------------------------------------------------------------------------
273 // 2.1. general nonlinear expression
275 template<class T, class M, class Expr>
276 typename std::enable_if<
277  details::is_field_expr_v2_nonlinear_arg<Expr>::value
278  && ! details::is_field_convertible<Expr>::value
279  && ! details::is_field_function<Expr>::value
280  ,field_basic<T,M>
281 >::type
282 interpolate (const space_basic<T,M>& Xh, const Expr& expr)
283 {
285  typedef typename wrap_t::value_type result_t;
286  return details::interpolate_internal<T,M,wrap_t,result_t> (Xh, wrap_t(expr));
287 }
288 // 2.2. function & functor
290 template <class T, class M, class Expr>
291 inline
292 typename std::enable_if<
293  details::is_field_function<Expr>::value
294  ,field_basic<T,M>
295 >::type
296 interpolate (const space_basic<T,M>& Xh, const Expr& expr)
297 {
299  typedef typename wrap_t::value_type result_t;
300  return details::interpolate_internal<T,M,wrap_t,result_t> (Xh, wrap_t(expr));
301 }
302 
303 // 2.3. re-interpolation of fields and linear field expressions
304 // for change of mesh, of approx, ect
306 template<class T, class M>
307 field_basic<T,M>
308 interpolate (const space_basic<T,M>& X2h, const field_basic<T,M>& u1h);
309 
311 template <class T, class M, class Expr>
312 inline
313 typename std::enable_if<
314  details::is_field_convertible<Expr>::value
315  && ! details::is_field<Expr>::value
316  ,field_basic<T,M>
317 >::type
318 interpolate (const space_basic<T,M>& Xh, const Expr& expr)
319 {
320  return interpolate (Xh, field_basic<T,M>(expr));
321 }
322 
323 }// namespace rheolef
324 #endif // _RHEOLEF_INTERPOLATE_H
rheolef::geo_basic
generic mesh with rerefence counting
Definition: geo.h:1089
rheolef::details::interpolate_generic
field_basic< T, M > interpolate_generic(const space_basic< T, M > &Xh, const Expr &expr)
Definition: interpolate.h:121
rheolef::piola_on_pointset::initialize
void initialize(const basis_basic< T > &piola_basis, const quadrature< T > &quad, const integrate_option &iopt)
Definition: piola_on_pointset.h:190
warning_macro
#define warning_macro(message)
Definition: dis_macros.h:53
mkgeo_ball.expr
expr
Definition: mkgeo_ball.sh:361
rheolef::point_basic
Definition: point.h:87
check_macro
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
rheolef::details::interpolate_internal_check::operator()
field_basic< T, M > operator()(const space_basic< T, M > &Xh, const Expr &expr) const
Definition: interpolate.h:186
mkgeo_ball.b
int b
Definition: mkgeo_ball.sh:152
rheolef::details::interpolate_internal
field_basic< T, M > interpolate_internal(const space_basic< T, M > &Xh, const Expr &expr)
Definition: interpolate.h:261
rheolef::details::interpolate_internal_check< T, M, Expr, T, std::true_type >
Definition: interpolate.h:199
rheolef::details::field_expr_v2_nonlinear_terminal_function
Definition: field_expr_terminal.h:258
rheolef::details::field_expr_v2_nonlinear_terminal_wrapper_traits::type
Expr type
Definition: field_expr_terminal.h:98
rheolef::value
rheolef::std value
rheolef::space_basic
the finite element space
Definition: space.h:352
rheolef::tensor_basic
Definition: tensor.h:90
rheolef::geo_element
see the geo_element page for the full documentation
Definition: geo_element.h:102
rheolef::size_type
size_t size_type
Definition: basis_get.cc:76
rheolef::space_numbering::dis_idof
void dis_idof(const basis_basic< T > &b, const geo_size &gs, const geo_element &K, typename std::vector< size_type >::iterator dis_idof_tab)
Definition: space_numbering.cc:187
rheolef::undeterminated_basic
helper for generic field value_type: T, point_basic<T> or tensor_basic<T>
Definition: undeterminated.h:32
rheolef::geo_element::dis_ie
size_type dis_ie() const
Definition: geo_element.h:163
rheolef::field_basic::dis_dof_entry
dis_reference dis_dof_entry(size_type dis_idof)
Definition: field.cc:397
rheolef::basis_basic
Definition: basis.h:532
rheolef::reference_element
see the reference_element page for the full documentation
Definition: reference_element.h:66
rheolef::integrate_option
see the integrate_option page for the full documentation
Definition: integrate_option.h:125
rheolef::space_constant::tensor
@ tensor
Definition: space_constant.h:138
rheolef::space_constant::scalar
@ scalar
Definition: space_constant.h:136
rheolef::field_basic::dis_dof_update
void dis_dof_update() const
Definition: field.cc:410
rheolef::interpolate
field_basic< T, M > interpolate(const space_basic< T, M > &V2h, const field_basic< T, M > &u1h)
see the interpolate page for the full documentation
Definition: interpolate.cc:233
rheolef::reference_element::name
char name() const
Definition: reference_element.h:100
rheolef::type
rheolef::std type
rheolef::field_basic
Definition: field.h:235
rheolef::space_constant::unsymmetric_tensor
@ unsymmetric_tensor
Definition: space_constant.h:139
fatal_macro
#define fatal_macro(message)
Definition: dis_macros.h:33
rheolef::piola_on_pointset
Definition: piola_on_pointset.h:142
rheolef::space_constant::valued_name
const std::string & valued_name(valued_type valued_tag)
Definition: space_constant.cc:43
rheolef
This file is part of Rheolef.
Definition: compiler_eigen.h:37
error_macro
#define error_macro(message)
Definition: dis_macros.h:49
rheolef::space_constant::vector
@ vector
Definition: space_constant.h:137
size_type
field::size_type size_type
Definition: branch.cc:425
rheolef::details::interpolate_internal_check
Definition: interpolate.h:184
rheolef::geo_element::name
char name() const
Definition: geo_element.h:169
rheolef::space_constant::valued_tag_traits
Definition: space_constant.h:161
rheolef::std
Definition: vec_expr_v2.h:402
M
Expr1::memory_type M
Definition: vec_expr_v2.h:416
rheolef::field_basic::size_type
std::size_t size_type
Definition: field.h:239
T
Expr1::float_type T
Definition: field_expr.h:261
trace_macro
#define trace_macro(message)
Definition: dis_macros.h:111
rheolef::piola_on_pointset::get_piola
const Eigen::Matrix< piola< T >, Eigen::Dynamic, 1 > & get_piola(const geo_basic< T, M > &omega, const geo_element &K) const
Definition: piola_on_pointset.h:239
rheolef::piola
Definition: piola.h:67