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
23 //
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
Float phi(const point &nu, Float a, Float b)
see the Float page for the full documentation
see the basis page for the full documentation
see the point page for the full documentation
rep::const_iterator const_iterator
Definition: quadrature.h:195
const_iterator end(reference_element hat_K) const
Definition: quadrature.h:219
const_iterator begin(reference_element hat_K) const
Definition: quadrature.h:218
see the reference_element page for the full documentation
integrate_option quadrature_option
typename details::generic_binary_traits< BinaryFunction >::template result_hint< typename Expr1::result_type, typename Expr2::result_type >::type result_type
rheolef::std Function
see the tensor page for the full documentation
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
This file is part of Rheolef.
Float error_div_linf(const field_element &uk, const Function &u, size_t qorder=size_t(0))
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
Float error_l2(const field_element &uk, const Function &u, size_t qorder=size_t(0))
Float error_linf(const field_element &uk, const Function &u, size_t nsub=0)
Float measure(reference_element hat_K)
T norm2(const vec< T, M > &x)
norm2(x): see the expression page for the full documentation
Definition: vec.h:379
U tr(const tensor_basic< U > &a, size_t d=3)
Float generic_norm2(const Float &x)
Float error_div_l2(const field_element &uk, const Function &u, size_t qorder=size_t(0))
Definition: cavity_dg.h:29
Definition: phi.h:25
Value operator()(const point &hat_x) const
field_element_function_wrapper(const field_element &uk)
Eigen::Matrix< Float, Eigen::Dynamic, 1 > _dof
Definition: field_element.h:97
Value grad_evaluate(const point &hat_x) const
Definition: field_element.h:78
Float operator()(const point &hat_x) const
Definition: field_element.h:89
Float dof(size_t i) const
Definition: field_element.h:94
size_t ndof() const
Definition: field_element.h:93
void interpolate(const Function &f)
Definition: field_element.h:59
Value evaluate(const point &hat_x) const
Definition: field_element.h:66
field_element(const space_element &V)
Definition: field_element.h:53
space_element(std::string name, reference_element hat_K)
Definition: field_element.h:41
valued_type valued_tag() const
Definition: field_element.h:45
space_constant::valued_type valued_type
Definition: field_element.h:40
space_element(const basis &b, reference_element hat_K)
Definition: field_element.h:43
reference_element _hat_K
Definition: field_element.h:48
helper for generic field value_type: T, point_basic<T> or tensor_basic<T>
Definition: leveque.h:25
Float u(const point &x)