Loading...
Searching...
No Matches
interpolate.h
Go to the documentation of this file.
#ifndef _RHEOLEF_INTERPOLATE_H
#define _RHEOLEF_INTERPOLATE_H
//
// This file is part of Rheolef.
//
// Copyright (C) 2000-2009 Pierre Saramito <Pierre.Saramito@imag.fr>
//
// Rheolef is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 2 of the License, or
// (at your option) any later version.
//
// Rheolef is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with Rheolef; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
//
// =========================================================================
// author: Pierre.Saramito@imag.fr
// date: 15 september 2015
namespace rheolef {
} // namespace rheolef
//
// SUMMARY:
// 1. implementation
// 1.1. scalar-valued result field
// 1.2. vector-valued case
// 1.3. tensor-valued case
// 1.4. undeterminated-valued case: determined at run-time
// 1.5. interface of the internal interpolate() function
// 2. the interpolate() function
// 2.1. nonlinear expression and function/functor
// 2.2. re-interpolation of fields and linear field expressions
#include "rheolef/field.h"
#include "rheolef/field_expr.h"
#include "rheolef/field_expr_terminal.h"
namespace rheolef {
namespace details {
// --------------------------------------------------------------------------
// 1. implementation: general nonlinear expr
// --------------------------------------------------------------------------
// notes:
// * function template partial specialization is not allowed --> use class-function
// * interpolation on boundary domain spaces (geo_domain) could generate
// some communications (see interpolate_dom[234]_tst.cc for tests)
// * interpolation of functions and functor are be performed in all cases
// without communications, see below for this implementation.
// here we suppose a general nonlinear expression that is interpolated
// by using a loop on mesh elements
//
template<class T, class M, class Expr, class Result>
field_basic<T,M> interpolate_generic (
const space_basic<T,M>& Xh,
const Expr& expr0)
{
Expr expr = expr0; // could modify a local copy, eg call expr.initialize()
trace_macro ("Expr="<<pretty_typename_macro(Expr));
typedef typename field_basic<T,M>::size_type size_type;
check_macro (Xh.valued_tag() == space_constant::valued_tag_traits<Result>::value,
"interpolate: incompatible "<<Xh.valued()<<"-valued space and "
<< space_constant::valued_name (space_constant::valued_tag_traits<Result>::value) << "-valued expression");
std::vector<size_type> dis_idof;
Eigen::Matrix<Result,Eigen::Dynamic,1> value;
Eigen::Matrix<T,Eigen::Dynamic,1> udof;
field_basic<T,M> uh (Xh, std::numeric_limits<T>::max());
const geo_basic<T,M>& omega = Xh.get_geo();
const basis_basic<T>& b = Xh.get_basis();
const piola_fem<T>& pf = b.get_piola_fem();
trace_macro ("pf.transform_need_piola="<<pf.transform_need_piola());
integrate_option iopt;
piola_on_pointset<T> pops;
pops.initialize (omega.get_piola_basis(), b, iopt);
expr.initialize (Xh, pops, iopt);
expr.template valued_check<Result>();
for (typename geo_basic<T,M>::const_iterator
iter_ie = omega.begin(),
last_ie = omega.end(); iter_ie != last_ie; ++iter_ie) {
const geo_element& K = *iter_ie;
reference_element hat_K = K;
// 1) get u values at nodes of K
expr.evaluate (omega, K, value);
// 2) u values at nodes of K -> udofs on K
if (pf.transform_need_piola()) {
const Eigen::Matrix<piola<T>,Eigen::Dynamic,1>& piola = pops.get_piola (omega, K);
for (size_type loc_inod = 0, loc_nnod = value.size(); loc_inod < loc_nnod; ++loc_inod) {
// be carefull: piola_fem::inv_transform should support inplace call in the "value" arg
#ifndef TO_CLEAN
Result old_value = value[loc_inod];
#endif // TO_CLEAN
pf.inv_transform (piola[loc_inod], value[loc_inod], value[loc_inod]);
#ifndef TO_CLEAN
trace_macro ("inv_transf(K="<<K.name()<<K.dis_ie()<<",loc_inod="<<loc_inod<<"): old_value="
<<old_value<<", value="<<value[loc_inod]);
#endif // TO_CLEAN
}
}
b.compute_dofs (hat_K, value, udof);
// 3) copy all local dofs into the global field
Xh.dis_idof (K, dis_idof);
check_macro (b.ndof(hat_K) == dis_idof.size() && b.ndof(hat_K) == size_type(udof.size()),
"invalid sizes: basis("<<b.name()<<").size("<<hat_K.name()<<") = "<< b.ndof(hat_K)
<<", dis_idof.size="<<dis_idof.size()<<", udof.size="<<udof.size());
for (size_type loc_idof = 0, loc_ndof = udof.size(); loc_idof < loc_ndof; ++loc_idof) {
uh.dis_dof_entry (dis_idof[loc_idof]) = udof [loc_idof];
#ifndef TO_CLEAN
trace_macro ("uh(K="<<K.name()<<K.dis_ie()<<",loc_idof="<<loc_idof<<") = uh(dis_idof="<<dis_idof[loc_idof]
<< ") = " << udof [loc_idof]);
#endif // TO_CLEAN
}
}
uh.dis_dof_update();
trace_macro ("interpolate done");
return uh;
}
template<class T, class M, class Expr, class Result, class Status = typename details::is_equal<Result,typename Expr::value_type>::type>
struct interpolate_internal_check {
field_basic<T,M>
operator() (
const space_basic<T,M>& Xh,
const Expr& expr) const
{
trace_macro ("Expr="<<pretty_typename_macro(Expr));
trace_macro ("Result="<<typename_macro(Result));
trace_macro ("Status="<<typename_macro(Status));
trace_macro ("Expr::value_type="<<typename_macro(typename Expr::value_type));
fatal_macro ("invalid type resolution");
return field_basic<T,M>();
}};
// 1.1. scalar-valued result field
template<class T, class M, class Expr>
struct interpolate_internal_check<T,M,Expr,T,std::true_type> {
field_basic<T,M>
operator() (
const space_basic<T,M>& Xh,
const Expr& expr) const
{
return interpolate_generic<T,M,Expr,T>(Xh,expr);
}};
// 1.2. vector-valued case
template<class T, class M, class Expr>
struct interpolate_internal_check<T,M,Expr,point_basic<T>,std::true_type> {
field_basic<T,M>
operator() (
const space_basic<T,M>& Xh,
const Expr& expr) const
{
return interpolate_generic<T,M,Expr,point_basic<T>>(Xh,expr);
}};
// 1.3. tensor-valued case
template<class T, class M, class Expr>
struct interpolate_internal_check<T,M,Expr,tensor_basic<T>,std::true_type> {
field_basic<T,M>
operator() (
const space_basic<T,M>& Xh,
const Expr& expr) const
{
return interpolate_generic<T,M,Expr,tensor_basic<T>>(Xh,expr);
}};
// 1.4. undeterminated-valued case: determined at run-time
template<class T, class M, class Expr, class Status>
struct interpolate_internal_check<T,M,Expr,undeterminated_basic<T>,Status> {
field_basic<T,M>
operator() (
const space_basic<T,M>& Xh,
const Expr& expr) const
{
switch (expr.valued_tag()) {
case space_constant::scalar: {
interpolate_internal_check<T,M,Expr,T,std::true_type> eval;
return eval (Xh, expr);
}
case space_constant::vector: {
interpolate_internal_check<T,M,Expr,point_basic<T>,std::true_type> eval;
return eval (Xh, expr);
}
case space_constant::tensor:
case space_constant::unsymmetric_tensor: {
interpolate_internal_check<T,M,Expr,tensor_basic<T>,std::true_type> eval;
return eval (Xh, expr);
}
default:
warning_macro ("Expr="<<pretty_typename_macro(Expr));
warning_macro ("Status="<<typename_macro(Status));
error_macro ("unexpected `"
<< space_constant::valued_name (expr.valued_tag())
<< "' valued expression");
return field_basic<T,M>();
}
}};
// 1.5. interface of the internal interpolate() function
template<class T, class M, class Expr, class Result>
field_basic<T,M>
interpolate_internal (
const space_basic<T,M>& Xh,
const Expr& expr)
{
interpolate_internal_check<T,M,Expr,Result> eval;
return eval (Xh,expr);
}
} // namespace details
// --------------------------------------------------------------------------
// 2. the interpolate() function
// --------------------------------------------------------------------------
// 2.1. general nonlinear expression
template<class T, class M, class Expr>
typename std::enable_if<
std::conjunction<
details::is_field_expr_v2_nonlinear_arg<Expr>
,std::negation<
std::disjunction<
details::is_field<Expr>
,details::has_field_rdof_interface<Expr>
,details::is_field_function<Expr>
>
>
>::value
,field_basic<T,M>
>::type
interpolate (const space_basic<T,M>& Xh, const Expr& expr)
{
typedef typename details::field_expr_v2_nonlinear_terminal_wrapper_traits<Expr>::type wrap_t;
typedef typename wrap_t::value_type result_t;
return details::interpolate_internal<T,M,wrap_t,result_t> (Xh, wrap_t(expr));
}
// 2.2. re-interpolation of fields and linear field expressions
// for change of mesh, of approx, ect
template <class T, class M, class Expr>
inline
typename std::enable_if<
details::has_field_rdof_interface<Expr>::value
&& ! details::is_field<Expr>::value
,field_basic<T,M>
>::type
interpolate (const space_basic<T,M>& Xh, const Expr& expr)
{
return interpolate (Xh, field_basic<T,M>(expr));
}
template<class T, class M>
field_basic<T,M>
interpolate (const space_basic<T,M>& X2h, const field_basic<T,M>& u1h);
// 2.3. function & functor
template <class T, class M, class Expr>
inline
typename std::enable_if<
details::is_field_function<Expr>::value
,field_basic<T,M>
>::type
interpolate (const space_basic<T,M>& Xh, const Expr& expr)
{
typedef typename details::field_expr_v2_nonlinear_terminal_function<Expr> wrap_t;
typedef typename wrap_t::value_type result_t;
return details::interpolate_internal<T,M,wrap_t,result_t> (Xh, wrap_t(expr));
}
}// namespace rheolef
#endif // _RHEOLEF_INTERPOLATE_H
rheolef::basis_basic
Definition basis.h:566
rheolef::details::field_expr_v2_nonlinear_terminal_function
Definition field_expr_terminal.h:244
rheolef::field_basic
Definition field.h:219
rheolef::field_basic::dis_dof_update
void dis_dof_update(const SetOp &=SetOp()) const
Definition field.h:763
rheolef::field_basic::dis_dof_entry
dis_reference dis_dof_entry(size_type dis_idof)
Definition field.cc:398
rheolef::field_basic::size_type
std::size_t size_type
Definition field.h:225
rheolef::geo_basic
generic mesh with rerefence counting
Definition geo.h:1089
rheolef::geo_element
see the geo_element page for the full documentation
Definition geo_element.h:102
rheolef::geo_element::dis_ie
size_type dis_ie() const
Definition geo_element.h:163
rheolef::geo_element::name
char name() const
Definition geo_element.h:169
rheolef::integrate_option
see the integrate_option page for the full documentation
Definition integrate_option.h:125
rheolef::piola_on_pointset
Definition piola_on_pointset.h:142
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
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::point_basic
Definition point.h:86
rheolef::reference_element
see the reference_element page for the full documentation
Definition reference_element.h:66
rheolef::reference_element::name
char name() const
Definition reference_element.h:100
rheolef::space_basic
the finite element space
Definition space.h:382
rheolef::tensor_basic
Definition tensor.h:90
trace_macro
#define trace_macro(message)
Definition dis_macros.h:111
error_macro
#define error_macro(message)
Definition dis_macros.h:49
fatal_macro
#define fatal_macro(message)
Definition dis_macros.h:33
warning_macro
#define warning_macro(message)
Definition dis_macros.h:53
T
Expr1::float_type T
Definition field_expr.h:230
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
Definition compiler_eigen.h:39
rheolef::details::interpolate_internal
field_basic< T, M > interpolate_internal(const space_basic< T, M > &Xh, const Expr &expr)
Definition interpolate.h:264
rheolef::details::interpolate_generic
field_basic< T, M > interpolate_generic(const space_basic< T, M > &Xh, const Expr &expr0)
Definition interpolate.h:123
rheolef::space_constant::tensor
@ tensor
Definition space_constant.h:136
rheolef::space_constant::scalar
@ scalar
Definition space_constant.h:134
rheolef::space_constant::unsymmetric_tensor
@ unsymmetric_tensor
Definition space_constant.h:137
rheolef::space_constant::vector
@ vector
Definition space_constant.h:135
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:39
rheolef::valued_check
bool valued_check() const
Definition field_expr_recursive.h:956
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
std
STL namespace.
rheolef::details::field_expr_v2_nonlinear_terminal_wrapper_traits::type
Expr type
Definition field_expr_terminal.h:84
rheolef::details::interpolate_internal_check
Definition interpolate.h:187
rheolef::details::interpolate_internal_check::operator()
field_basic< T, M > operator()(const space_basic< T, M > &Xh, const Expr &expr) const
Definition interpolate.h:189
rheolef::piola
Definition piola.h:67
rheolef::space_constant::valued_tag_traits::value
static const valued_type value
Definition space_constant.h:159
rheolef::undeterminated_basic
helper for generic field value_type: T, point_basic<T> or tensor_basic<T>
Definition undeterminated.h:32
M
Expr1::memory_type M
Definition vec_expr_v2.h:416