Rheolef  7.1
an efficient C++ finite element environment
field_element.h
Go to the documentation of this file.
1 #ifndef _RHEO_FIELD_ELEMENT_H
2 #define _RHEO_FIELD_ELEMENT_H
3 //
24 // for checking the spectral convergence exp(-ak)
25 // on the reference element
26 //
27 // author: Pierre.Saramito@imag.fr
28 //
29 // date: 1 october 2017
30 //
31 #include "rheolef/quadrature.h"
32 #include "rheolef/basis.h"
33 #include "rheolef/field_expr_utilities.h"
34 #include "equispaced.icc"
35 
36 namespace rheolef {
37 // ---------------------------------------------------
38 struct space_element {
39 // ---------------------------------------------------
41  space_element (std::string name, reference_element hat_K)
42  : _b(name), _hat_K(hat_K) {}
44  : _b(b), _hat_K(hat_K) {}
45  valued_type valued_tag() const { return _b.valued_tag(); }
46 // data:
49 };
50 // ---------------------------------------------------
51 struct field_element {
52 // ---------------------------------------------------
54  : _V(V), _dof()
55  {
56  _dof.resize (_V._b.ndof (_V._hat_K));
57  }
58  template <class Function>
59  void interpolate (const Function& f)
60  {
61  _V._b.compute_dof (_V._hat_K, f, _dof);
62  }
63  void interpolate (const field_element& uk);
64 
65  template<class Value>
66  Value evaluate (const point& hat_x) const
67  {
68  Eigen::Matrix<Value,Eigen::Dynamic,1> phi;
69  _V._b.evaluate (_V._hat_K, hat_x, phi);
70  check_macro (phi.size() == _dof.size(), "invalid sizes");
71  Value result = Value(0);
72  for (size_t loc_idof = 0, loc_ndof = _dof.size(); loc_idof < loc_ndof; ++loc_idof) {
73  result = result + _dof[loc_idof]*phi[loc_idof];
74  }
75  return result;
76  }
77  template<class Value>
78  Value grad_evaluate (const point& hat_x) const
79  {
80  Eigen::Matrix<Value,Eigen::Dynamic,1> grad_phi;
81  _V._b.grad_evaluate (_V._hat_K, hat_x, grad_phi);
82  check_macro (grad_phi.size() == _dof.size(), "invalid sizes");
83  Value result = Value();
84  for (size_t loc_idof = 0, loc_ndof = _dof.size(); loc_idof < loc_ndof; ++loc_idof) {
85  result = result + _dof[loc_idof]*grad_phi[loc_idof];
86  }
87  return result;
88  }
89  Float operator() (const point& hat_x) const
90  {
91  return (*this).template evaluate<Float> (hat_x);
92  }
93  size_t ndof() const { return _dof.size(); }
94  Float dof(size_t i) const { return _dof(i); }
95 // data:
97  Eigen::Matrix<Float,Eigen::Dynamic,1> _dof;
98 };
99 // re-interpolate a field_element, for commutation checks
100 template<class Value>
102  Value operator() (const point& hat_x) const
103  { return _uk.template evaluate<Value> (hat_x); }
106 };
107 inline
108 void
110 {
111  _V = fk._V;
112  switch (fk._V.valued_tag()) {
113  case space_constant::scalar: {
115  (*this).interpolate (f);
116  break;
117  }
118  case space_constant::vector: {
120  (*this).interpolate (f);
121  break;
122  }
123  default: error_macro ("unsupported valued field_element::interpolate(field_element)");
124  }
125 }
126 // ---------------------------------------------------
127 // error_l2
128 // ---------------------------------------------------
129 Float generic_norm2 (const Float& x) { return sqr(x); }
130 Float generic_norm2 (const point& x) { return norm2(x); }
131 Float generic_norm2 (const tensor& x) { return norm2(x); }
132 
133 template <class Function>
134 Float
135 error_l2 (const field_element& uk, const Function& u, size_t qorder = size_t(0)) {
137  if (qorder == size_t(0)) qorder = max(size_t(5), 2*uk._V._b.degree() + 2);
138  quadrature_option qopt;
139  qopt.set_order (qorder);;
140  quadrature<Float> quad (qopt);
141  Float err2 = 0;
142  for (typename quadrature<Float>::const_iterator iter_q = quad.begin(uk._V._hat_K),
143  last_q = quad.end(uk._V._hat_K); iter_q != last_q; iter_q++) {
144  const point& hat_xq = (*iter_q).x;
145  Float hat_wq = (*iter_q).w;
146  err2 += generic_norm2(u(hat_xq) - uk.template evaluate<result_type>(hat_xq))*hat_wq;
147  }
148  return sqrt(err2/measure(uk._V._hat_K));
149 }
150 // ---------------------------------------------------
151 // error_linf
152 // ---------------------------------------------------
153 template <class Function>
154 Float
155 error_linf (const field_element& uk, const Function& u, size_t nsub = 0) {
157  if (nsub == 0) nsub = max(size_t(5), 2*uk._V._b.degree() + 2);
158  Eigen::Matrix<point_basic<Float>,Eigen::Dynamic,1> hat_xnod;
159  pointset_lagrange_equispaced (uk._V._hat_K, nsub, hat_xnod);
160  Float err2 = 0;
161  for (size_t loc_inod = 0, loc_nnod = hat_xnod.size(); loc_inod < loc_nnod; ++loc_inod) {
162  err2 = max (err2, generic_norm2(uk.template evaluate<result_type>(hat_xnod[loc_inod]) - u(hat_xnod[loc_inod])));
163  }
164  return sqrt(err2);
165 }
166 // ---------------------------------------------------
167 // error_div_l2
168 // ---------------------------------------------------
169 template <class Function>
170 Float
171 error_div_l2 (const field_element& uk, const Function& u, size_t qorder = size_t(0)) {
173  using grad_type = typename space_constant::rank_up<result_type>::type;
174  using div_type = typename space_constant::rank_down<result_type>::type;
175  if (qorder == size_t(0)) qorder = max(size_t(5), 2*uk._V._b.degree() + 2);
176  quadrature_option qopt;
177  qopt.set_order (qorder);
178  quadrature<Float> quad (qopt);
179  Float err2 = 0;
180  for (typename quadrature<Float>::const_iterator iter_q = quad.begin(uk._V._hat_K),
181  last_q = quad.end(uk._V._hat_K); iter_q != last_q; iter_q++) {
182  const point& hat_xq = (*iter_q).x;
183  Float hat_wq = (*iter_q).w;
184  grad_type grad_u_xq = uk.template grad_evaluate<grad_type>(hat_xq);
185  div_type div_u_xq = tr(grad_u_xq);
186  err2 += generic_norm2(u.div(hat_xq) - div_u_xq)*hat_wq;
187  }
188  return sqrt(err2/measure(uk._V._hat_K));
189 }
190 // ---------------------------------------------------
191 // error_div_linf
192 // ---------------------------------------------------
193 template <class Function>
194 Float
195 error_div_linf (const field_element& uk, const Function& u, size_t qorder = size_t(0)) {
197  using grad_type = typename space_constant::rank_up<result_type>::type;
198  using div_type = typename space_constant::rank_down<result_type>::type;
199  if (qorder == size_t(0)) qorder = max(size_t(5), 2*uk._V._b.degree() + 2);
200  quadrature_option qopt;
201  qopt.set_order (qorder);
202  quadrature<Float> quad (qopt);
203  Float err2 = 0;
204  for (typename quadrature<Float>::const_iterator iter_q = quad.begin(uk._V._hat_K),
205  last_q = quad.end(uk._V._hat_K); iter_q != last_q; iter_q++) {
206  const point& hat_xq = (*iter_q).x;
207  Float hat_wq = (*iter_q).w;
208  grad_type grad_u_xq = uk.template grad_evaluate<grad_type>(hat_xq);
209  div_type div_u_xq = tr(grad_u_xq);
210  err2 = max (err2, generic_norm2(u.div(hat_xq) - div_u_xq));
211  }
212  return sqrt(err2);
213 }
214 
215 } // namespace rheolef
216 #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