Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
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
26namespace 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"
107namespace rheolef {
108
109namespace 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//
122template<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()
128trace_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();
140trace_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>();
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();
183trace_macro ("interpolate done");
184 return uh;
185}
186template<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
201template<class T, class M, class Expr>
202struct 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
211template<class T, class M, class Expr>
212struct 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
221template<class T, class M, class Expr>
222struct 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
231template<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()) {
241 return eval (Xh, expr);
242 }
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
262template<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
278template<class T, class M, class Expr>
279typename 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
292interpolate (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
301template <class T, class M, class Expr>
302inline
303typename std::enable_if<
304 details::has_field_rdof_interface<Expr>::value
305 && ! details::is_field<Expr>::value
306 ,field_basic<T,M>
307>::type
308interpolate (const space_basic<T,M>& Xh, const Expr& expr)
309{
310 return interpolate (Xh, field_basic<T,M>(expr));
311}
313template<class T, class M>
314field_basic<T,M>
315interpolate (const space_basic<T,M>& X2h, const field_basic<T,M>& u1h);
316
317// 2.3. function & functor
319template <class T, class M, class Expr>
320inline
321typename std::enable_if<
322 details::is_field_function<Expr>::value
323 ,field_basic<T,M>
324>::type
325interpolate (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
std::size_t size_type
Definition field.h:225
generic mesh with rerefence counting
Definition geo.h:1089
see the geo_element page for the full documentation
size_type dis_ie() const
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
the finite element space
Definition space.h:382
#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_internal(const space_basic< T, M > &Xh, const Expr &expr)
field_basic< T, M > interpolate_generic(const space_basic< T, M > &Xh, const Expr &expr0)
const std::string & valued_name(valued_type valued_tag)
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
STL namespace.
field_basic< T, M > operator()(const space_basic< T, M > &Xh, const Expr &expr) const
helper for generic field value_type: T, point_basic<T> or tensor_basic<T>
Expr1::memory_type M