Rheolef  7.2
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.h"
105 #include "rheolef/field_expr.h"
106 #include "rheolef/field_expr_terminal.h"
107 namespace rheolef {
108 
109 namespace details {
110 // --------------------------------------------------------------------------
111 // 1. implementation: general nonlinear expr
112 // --------------------------------------------------------------------------
113 // notes:
114 // * function template partial specialization is not allowed --> use class-function
115 // * interpolation on boundary domain spaces (geo_domain) could generate
116 // some communications (see interpolate_dom[234]_tst.cc for tests)
117 // * interpolation of functions and functor are be performed in all cases
118 // without communications, see below for this implementation.
119 // here we suppose a general nonlinear expression that is interpolated
120 // by using a loop on mesh elements
121 //
122 template<class T, class M, class Expr, class Result>
124  const space_basic<T,M>& Xh,
125  const Expr& expr0)
126 {
127  Expr expr = expr0; // could modify a local copy, eg call expr.initialize()
128 trace_macro ("Expr="<<pretty_typename_macro(Expr));
129  typedef typename field_basic<T,M>::size_type size_type;
131  "interpolate: incompatible "<<Xh.valued()<<"-valued space and "
133  std::vector<size_type> dis_idof;
134  Eigen::Matrix<Result,Eigen::Dynamic,1> value;
135  Eigen::Matrix<T,Eigen::Dynamic,1> udof;
136  field_basic<T,M> uh (Xh, std::numeric_limits<T>::max());
137  const geo_basic<T,M>& omega = Xh.get_geo();
138  const basis_basic<T>& b = Xh.get_basis();
139  const piola_fem<T>& pf = b.get_piola_fem();
140 trace_macro ("pf.transform_need_piola="<<pf.transform_need_piola());
141  integrate_option iopt;
143  pops.initialize (omega.get_piola_basis(), b, iopt);
144  expr.initialize (Xh, pops, iopt);
145  expr.template valued_check<Result>();
146  for (typename geo_basic<T,M>::const_iterator
147  iter_ie = omega.begin(),
148  last_ie = omega.end(); iter_ie != last_ie; ++iter_ie) {
149  const geo_element& K = *iter_ie;
150  reference_element hat_K = K;
151  // 1) get u values at nodes of K
152  expr.evaluate (omega, K, value);
153  // 2) u values at nodes of K -> udofs on K
154  if (pf.transform_need_piola()) {
155  const Eigen::Matrix<piola<T>,Eigen::Dynamic,1>& piola = pops.get_piola (omega, K);
156  for (size_type loc_inod = 0, loc_nnod = value.size(); loc_inod < loc_nnod; ++loc_inod) {
157  // be carefull: piola_fem::inv_transform should support inplace call in the "value" arg
158 #ifndef TO_CLEAN
159  Result old_value = value[loc_inod];
160 #endif // TO_CLEAN
161  pf.inv_transform (piola[loc_inod], value[loc_inod], value[loc_inod]);
162 #ifndef TO_CLEAN
163  trace_macro ("inv_transf(K="<<K.name()<<K.dis_ie()<<",loc_inod="<<loc_inod<<"): old_value="
164  <<old_value<<", value="<<value[loc_inod]);
165 #endif // TO_CLEAN
166  }
167  }
168  b.compute_dofs (hat_K, value, udof);
169  // 3) copy all local dofs into the global field
170  Xh.dis_idof (K, dis_idof);
171  check_macro (b.ndof(hat_K) == dis_idof.size() && b.ndof(hat_K) == size_type(udof.size()),
172  "invalid sizes: basis("<<b.name()<<").size("<<hat_K.name()<<") = "<< b.ndof(hat_K)
173  <<", dis_idof.size="<<dis_idof.size()<<", udof.size="<<udof.size());
174  for (size_type loc_idof = 0, loc_ndof = udof.size(); loc_idof < loc_ndof; ++loc_idof) {
175  uh.dis_dof_entry (dis_idof[loc_idof]) = udof [loc_idof];
176 #ifndef TO_CLEAN
177  trace_macro ("uh(K="<<K.name()<<K.dis_ie()<<",loc_idof="<<loc_idof<<") = uh(dis_idof="<<dis_idof[loc_idof]
178  << ") = " << udof [loc_idof]);
179 #endif // TO_CLEAN
180  }
181  }
182  uh.dis_dof_update();
183 trace_macro ("interpolate done");
184  return uh;
185 }
186 template<class T, class M, class Expr, class Result, class Status = typename details::is_equal<Result,typename Expr::value_type>::type>
190  const space_basic<T,M>& Xh,
191  const Expr& expr) const
192 {
193  trace_macro ("Expr="<<pretty_typename_macro(Expr));
194  trace_macro ("Result="<<typename_macro(Result));
195  trace_macro ("Status="<<typename_macro(Status));
196  trace_macro ("Expr::value_type="<<typename_macro(typename Expr::value_type));
197  fatal_macro ("invalid type resolution");
198  return field_basic<T,M>();
199 }};
200 // 1.1. scalar-valued result field
201 template<class T, class M, class Expr>
202 struct interpolate_internal_check<T,M,Expr,T,std::true_type> {
205  const space_basic<T,M>& Xh,
206  const Expr& expr) const
207 {
208  return interpolate_generic<T,M,Expr,T>(Xh,expr);
209 }};
210 // 1.2. vector-valued case
211 template<class T, class M, class Expr>
212 struct interpolate_internal_check<T,M,Expr,point_basic<T>,std::true_type> {
215  const space_basic<T,M>& Xh,
216  const Expr& expr) const
217 {
218  return interpolate_generic<T,M,Expr,point_basic<T>>(Xh,expr);
219 }};
220 // 1.3. tensor-valued case
221 template<class T, class M, class Expr>
222 struct interpolate_internal_check<T,M,Expr,tensor_basic<T>,std::true_type> {
225  const space_basic<T,M>& Xh,
226  const Expr& expr) const
227 {
228  return interpolate_generic<T,M,Expr,tensor_basic<T>>(Xh,expr);
229 }};
230 // 1.4. undeterminated-valued case: determined at run-time
231 template<class T, class M, class Expr, class Status>
235  const space_basic<T,M>& Xh,
236  const Expr& expr) const
237 {
238  switch (expr.valued_tag()) {
239  case space_constant::scalar: {
241  return eval (Xh, expr);
242  }
243  case space_constant::vector: {
245  return eval (Xh, expr);
246  }
250  return eval (Xh, expr);
251  }
252  default:
253  warning_macro ("Expr="<<pretty_typename_macro(Expr));
254  warning_macro ("Status="<<typename_macro(Status));
255  error_macro ("unexpected `"
256  << space_constant::valued_name (expr.valued_tag())
257  << "' valued expression");
258  return field_basic<T,M>();
259  }
260 }};
261 // 1.5. interface of the internal interpolate() function
262 template<class T, class M, class Expr, class Result>
265  const space_basic<T,M>& Xh,
266  const Expr& expr)
267 {
269  return eval (Xh,expr);
270 }
271 
272 } // namespace details
273 // --------------------------------------------------------------------------
274 // 2. the interpolate() function
275 // --------------------------------------------------------------------------
276 // 2.1. general nonlinear expression
278 template<class T, class M, class Expr>
279 typename std::enable_if<
280  std::conjunction<
281  details::is_field_expr_v2_nonlinear_arg<Expr>
282  ,std::negation<
283  std::disjunction<
284  details::is_field<Expr>
285  ,details::has_field_rdof_interface<Expr>
286  ,details::is_field_function<Expr>
287  >
288  >
289  >::value
290  ,field_basic<T,M>
291 >::type
292 interpolate (const space_basic<T,M>& Xh, const Expr& expr)
293 {
295  typedef typename wrap_t::value_type result_t;
296  return details::interpolate_internal<T,M,wrap_t,result_t> (Xh, wrap_t(expr));
297 }
298 // 2.2. re-interpolation of fields and linear field expressions
299 // for change of mesh, of approx, ect
301 template <class T, class M, class Expr>
302 inline
303 typename std::enable_if<
304  details::has_field_rdof_interface<Expr>::value
305  && ! details::is_field<Expr>::value
306  ,field_basic<T,M>
307 >::type
308 interpolate (const space_basic<T,M>& Xh, const Expr& expr)
309 {
310  return interpolate (Xh, field_basic<T,M>(expr));
311 }
313 template<class T, class M>
314 field_basic<T,M>
315 interpolate (const space_basic<T,M>& X2h, const field_basic<T,M>& u1h);
316 
317 // 2.3. function & functor
319 template <class T, class M, class Expr>
320 inline
321 typename std::enable_if<
322  details::is_field_function<Expr>::value
323  ,field_basic<T,M>
324 >::type
325 interpolate (const space_basic<T,M>& Xh, const Expr& expr)
326 {
328  typedef typename wrap_t::value_type result_t;
329  return details::interpolate_internal<T,M,wrap_t,result_t> (Xh, wrap_t(expr));
330 }
331 
332 
333 }// namespace rheolef
334 #endif // _RHEOLEF_INTERPOLATE_H
field::size_type size_type
Definition: branch.cc:430
void dis_dof_update(const SetOp &=SetOp()) const
Definition: field.h:763
dis_reference dis_dof_entry(size_type dis_idof)
Definition: field.cc:398
see the geo_element page for the full documentation
Definition: geo_element.h:102
size_type dis_ie() const
Definition: geo_element.h:163
char name() const
Definition: geo_element.h:169
see the integrate_option page for the full documentation
void initialize(const basis_basic< T > &piola_basis, const quadrature< T > &quad, const integrate_option &iopt)
const Eigen::Matrix< piola< T >, Eigen::Dynamic, 1 > & get_piola(const geo_basic< T, M > &omega, const geo_element &K) const
see the reference_element page for the full documentation
rheolef::std type
size_t size_type
Definition: basis_get.cc:76
rheolef::std value
#define trace_macro(message)
Definition: dis_macros.h:111
#define error_macro(message)
Definition: dis_macros.h:49
#define fatal_macro(message)
Definition: dis_macros.h:33
#define warning_macro(message)
Definition: dis_macros.h:53
Expr1::float_type T
Definition: field_expr.h:230
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
field_basic< T, M > interpolate_generic(const space_basic< T, M > &Xh, const Expr &expr0)
Definition: interpolate.h:123
field_basic< T, M > interpolate_internal(const space_basic< T, M > &Xh, const Expr &expr)
Definition: interpolate.h:264
const std::string & valued_name(valued_type valued_tag)
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)
This file is part of Rheolef.
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
field_basic< T, M > operator()(const space_basic< T, M > &Xh, const Expr &expr) const
Definition: interpolate.h:189
helper for generic field value_type: T, point_basic<T> or tensor_basic<T>
Expr1::memory_type M
Definition: vec_expr_v2.h:416