field_element.h
Go to the documentation of this file.
#ifndef _RHEO_FIELD_ELEMENT_H
#define _RHEO_FIELD_ELEMENT_H
//
// for checking the spectral convergence exp(-ak)
// on the reference element
//
// author: Pierre.Saramito@imag.fr
//
// date: 1 october 2017
//
#include "rheolef/quadrature.h"
#include "rheolef/basis.h"
#include "rheolef/field_expr_utilities.h"
#include "equispaced.icc"
namespace rheolef {
// ---------------------------------------------------
struct space_element {
// ---------------------------------------------------
typedef space_constant::valued_type valued_type;
space_element (std::string name, reference_element hat_K)
: _b(name), _hat_K(hat_K) {}
space_element (const basis& b, reference_element hat_K)
: _b(b), _hat_K(hat_K) {}
valued_type valued_tag() const { return _b.valued_tag(); }
// data:
basis _b;
reference_element _hat_K;
};
// ---------------------------------------------------
struct field_element {
// ---------------------------------------------------
field_element (const space_element& V)
: _V(V), _dof()
{
_dof.resize (_V._b.ndof (_V._hat_K));
}
template <class Function>
void interpolate (const Function& f)
{
_V._b.compute_dof (_V._hat_K, f, _dof);
}
void interpolate (const field_element& uk);
template<class Value>
Value evaluate (const point& hat_x) const
{
Eigen::Matrix<Value,Eigen::Dynamic,1> phi;
_V._b.evaluate (_V._hat_K, hat_x, phi);
check_macro (phi.size() == _dof.size(), "invalid sizes");
Value result = Value(0);
for (size_t loc_idof = 0, loc_ndof = _dof.size(); loc_idof < loc_ndof; ++loc_idof) {
result = result + _dof[loc_idof]*phi[loc_idof];
}
return result;
}
template<class Value>
Value grad_evaluate (const point& hat_x) const
{
Eigen::Matrix<Value,Eigen::Dynamic,1> grad_phi;
_V._b.grad_evaluate (_V._hat_K, hat_x, grad_phi);
check_macro (grad_phi.size() == _dof.size(), "invalid sizes");
Value result = Value();
for (size_t loc_idof = 0, loc_ndof = _dof.size(); loc_idof < loc_ndof; ++loc_idof) {
result = result + _dof[loc_idof]*grad_phi[loc_idof];
}
return result;
}
Float operator() (const point& hat_x) const
{
return (*this).template evaluate<Float> (hat_x);
}
size_t ndof() const { return _dof.size(); }
Float dof(size_t i) const { return _dof(i); }
// data:
space_element _V;
Eigen::Matrix<Float,Eigen::Dynamic,1> _dof;
};
// re-interpolate a field_element, for commutation checks
template<class Value>
struct field_element_function_wrapper {
Value operator() (const point& hat_x) const
{ return _uk.template evaluate<Value> (hat_x); }
field_element_function_wrapper (const field_element& uk) : _uk(uk) {}
field_element _uk;
};
inline
void
field_element::interpolate (const field_element& fk)
{
_V = fk._V;
switch (fk._V.valued_tag()) {
case space_constant::scalar: {
field_element_function_wrapper<Float> f (fk);
(*this).interpolate (f);
break;
}
case space_constant::vector: {
field_element_function_wrapper<point> f (fk);
(*this).interpolate (f);
break;
}
default: error_macro ("unsupported valued field_element::interpolate(field_element)");
}
}
// ---------------------------------------------------
// error_l2
// ---------------------------------------------------
Float generic_norm2 (const Float& x) { return sqr(x); }
Float generic_norm2 (const point& x) { return norm2(x); }
Float generic_norm2 (const tensor& x) { return norm2(x); }
template <class Function>
Float
error_l2 (const field_element& uk, const Function& u, size_t qorder = size_t(0)) {
using result_type = typename details::function_traits<Function>::result_type;
if (qorder == size_t(0)) qorder = max(size_t(5), 2*uk._V._b.degree() + 2);
quadrature_option qopt;
qopt.set_order (qorder);;
quadrature<Float> quad (qopt);
Float err2 = 0;
for (typename quadrature<Float>::const_iterator iter_q = quad.begin(uk._V._hat_K),
last_q = quad.end(uk._V._hat_K); iter_q != last_q; iter_q++) {
const point& hat_xq = (*iter_q).x;
Float hat_wq = (*iter_q).w;
err2 += generic_norm2(u(hat_xq) - uk.template evaluate<result_type>(hat_xq))*hat_wq;
}
return sqrt(err2/measure(uk._V._hat_K));
}
// ---------------------------------------------------
// error_linf
// ---------------------------------------------------
template <class Function>
Float
error_linf (const field_element& uk, const Function& u, size_t nsub = 0) {
using result_type = typename details::function_traits<Function>::result_type;
if (nsub == 0) nsub = max(size_t(5), 2*uk._V._b.degree() + 2);
Eigen::Matrix<point_basic<Float>,Eigen::Dynamic,1> hat_xnod;
pointset_lagrange_equispaced (uk._V._hat_K, nsub, hat_xnod);
Float err2 = 0;
for (size_t loc_inod = 0, loc_nnod = hat_xnod.size(); loc_inod < loc_nnod; ++loc_inod) {
err2 = max (err2, generic_norm2(uk.template evaluate<result_type>(hat_xnod[loc_inod]) - u(hat_xnod[loc_inod])));
}
return sqrt(err2);
}
// ---------------------------------------------------
// error_div_l2
// ---------------------------------------------------
template <class Function>
Float
error_div_l2 (const field_element& uk, const Function& u, size_t qorder = size_t(0)) {
using result_type = typename details::function_traits<Function>::result_type;
using grad_type = typename space_constant::rank_up<result_type>::type;
using div_type = typename space_constant::rank_down<result_type>::type;
if (qorder == size_t(0)) qorder = max(size_t(5), 2*uk._V._b.degree() + 2);
quadrature_option qopt;
qopt.set_order (qorder);
quadrature<Float> quad (qopt);
Float err2 = 0;
for (typename quadrature<Float>::const_iterator iter_q = quad.begin(uk._V._hat_K),
last_q = quad.end(uk._V._hat_K); iter_q != last_q; iter_q++) {
const point& hat_xq = (*iter_q).x;
Float hat_wq = (*iter_q).w;
grad_type grad_u_xq = uk.template grad_evaluate<grad_type>(hat_xq);
div_type div_u_xq = tr(grad_u_xq);
err2 += generic_norm2(u.div(hat_xq) - div_u_xq)*hat_wq;
}
return sqrt(err2/measure(uk._V._hat_K));
}
// ---------------------------------------------------
// error_div_linf
// ---------------------------------------------------
template <class Function>
Float
error_div_linf (const field_element& uk, const Function& u, size_t qorder = size_t(0)) {
using result_type = typename details::function_traits<Function>::result_type;
using grad_type = typename space_constant::rank_up<result_type>::type;
using div_type = typename space_constant::rank_down<result_type>::type;
if (qorder == size_t(0)) qorder = max(size_t(5), 2*uk._V._b.degree() + 2);
quadrature_option qopt;
qopt.set_order (qorder);
quadrature<Float> quad (qopt);
Float err2 = 0;
for (typename quadrature<Float>::const_iterator iter_q = quad.begin(uk._V._hat_K),
last_q = quad.end(uk._V._hat_K); iter_q != last_q; iter_q++) {
const point& hat_xq = (*iter_q).x;
Float hat_wq = (*iter_q).w;
grad_type grad_u_xq = uk.template grad_evaluate<grad_type>(hat_xq);
div_type div_u_xq = tr(grad_u_xq);
err2 = max (err2, generic_norm2(u.div(hat_xq) - div_u_xq));
}
return sqrt(err2);
}
} // namespace rheolef
#endif // _RHEO_FIELD_ELEMENT_H
rheolef::space_element::valued_tag
valued_type valued_tag() const
Definition: field_element.h:45
tensor
see the tensor page for the full documentation
rheolef::space_element::space_element
space_element(std::string name, reference_element hat_K)
Definition: field_element.h:41
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::space_element::_b
basis _b
Definition: field_element.h:47
rheolef::pointset_lagrange_equispaced
void pointset_lagrange_equispaced(reference_element hat_K, size_t order_in, Eigen::Matrix< point_basic< T >, Eigen::Dynamic, 1 > &hat_xnod, size_t internal=0)
Definition: equispaced.icc:44
rheolef::integrate_option::set_order
void set_order(size_t r)
Definition: integrate_option.h:254
rheolef::Function
rheolef::std Function
rheolef::field_element::field_element
field_element(const space_element &V)
Definition: field_element.h:53
phi
Definition: phi.h:25
rheolef::quadrature::const_iterator
rep::const_iterator const_iterator
Definition: quadrature.h:195
mkgeo_ball.f
f
Definition: mkgeo_ball.sh:221
phi
Float phi(const point &nu, Float a, Float b)
Definition: burgers_flux_godunov.icc:25
rheolef::space_constant::valued_type
valued_type
Definition: space_constant.h:135
rheolef::undeterminated_basic
helper for generic field value_type: T, point_basic<T> or tensor_basic<T>
Definition: undeterminated.h:32
rheolef::field_element_function_wrapper::operator()
Value operator()(const point &hat_x) const
Definition: field_element.h:102
rheolef::error_div_linf
Float error_div_linf(const field_element &uk, const Function &u, size_t qorder=size_t(0))
Definition: field_element.h:195
rheolef::quadrature::end
const_iterator end(reference_element hat_K) const
Definition: quadrature.h:219
rheolef::norm2
T norm2(const vec< T, M > &x)
norm2(x): see the expression page for the full documentation
Definition: vec.h:379
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::result_type
typename details::generic_binary_traits< BinaryFunction >::template result_hint< typename Expr1::result_type, typename Expr2::result_type >::type result_type
Definition: field_expr_recursive.h:498
rheolef::space_constant::scalar
@ scalar
Definition: space_constant.h:136
rheolef::space_element::_hat_K
reference_element _hat_K
Definition: field_element.h:48
rheolef::field_element::dof
Float dof(size_t i) const
Definition: field_element.h:94
basis
see the basis page for the full documentation
rheolef::details::function_traits
Definition: field_expr_utilities.h:188
rheolef::generic_norm2
Float generic_norm2(const Float &x)
Definition: field_element.h:129
rheolef::field_element::interpolate
void interpolate(const Function &f)
Definition: field_element.h:59
rheolef::measure
Float measure(reference_element hat_K)
Definition: reference_element.cc:65
rheolef::field_element::grad_evaluate
Value grad_evaluate(const point &hat_x) const
Definition: field_element.h:78
rheolef::error_l2
Float error_l2(const field_element &uk, const Function &u, size_t qorder=size_t(0))
Definition: field_element.h:135
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_element::valued_type
space_constant::valued_type valued_type
Definition: field_element.h:40
u
Definition: leveque.h:25
rheolef::field_element_function_wrapper
Definition: field_element.h:101
rheolef::error_div_l2
Float error_div_l2(const field_element &uk, const Function &u, size_t qorder=size_t(0))
Definition: field_element.h:171
rheolef::space_element
Definition: field_element.h:38
Float
see the Float page for the full documentation
rheolef::field_element_function_wrapper::_uk
field_element _uk
Definition: field_element.h:105
point
see the point page for the full documentation
u
Float u(const point &x)
Definition: transmission_error.cc:26
rheolef::field_element_function_wrapper::field_element_function_wrapper
field_element_function_wrapper(const field_element &uk)
Definition: field_element.h:104
mkgeo_ball.b
b
Definition: mkgeo_ball.sh:152
equispaced.icc
rheolef::field_element::operator()
Float operator()(const point &hat_x) const
Definition: field_element.h:89
rheolef::space_constant::vector
@ vector
Definition: space_constant.h:137
rheolef::field_element::evaluate
Value evaluate(const point &hat_x) const
Definition: field_element.h:66
rheolef::field_element::ndof
size_t ndof() const
Definition: field_element.h:93
rheolef::field_element::_V
space_element _V
Definition: field_element.h:96
f
Definition: cavity_dg.h:29
rheolef::field_element
Definition: field_element.h:51
rheolef::quadrature
Definition: quadrature.h:185
rheolef::error_linf
Float error_linf(const field_element &uk, const Function &u, size_t nsub=0)
Definition: field_element.h:155
mkgeo_contraction.name
name
Definition: mkgeo_contraction.sh:133
rheolef::quadrature::begin
const_iterator begin(reference_element hat_K) const
Definition: quadrature.h:218
rheolef::field_element::_dof
Eigen::Matrix< Float, Eigen::Dynamic, 1 > _dof
Definition: field_element.h:97
rheolef::tr
U tr(const tensor_basic< U > &a, size_t d=3)
rheolef::space_element::space_element
space_element(const basis &b, reference_element hat_K)
Definition: field_element.h:43