Rheolef  7.2
an efficient C++ finite element environment
field_vf_assembly.h
Go to the documentation of this file.
1 #ifndef _RHEO_FIELD_VF_ASSEMBLY_H
2 #define _RHEO_FIELD_VF_ASSEMBLY_H
23 #include "rheolef/field.h"
24 #include "rheolef/test.h"
25 #include "rheolef/integrate_option.h"
26 #include "rheolef/field_expr_quadrature.h"
27 /*
28  let:
29  l(v) = int_omega expr(v) dx
30 
31  The integrals are evaluated over each element K of the domain
32  by using a quadrature formulae given by iopt
33 
34  expr(v) is a linear expression with respect to the test-function v
35 
36  The test-function v is replaced by each of the basis function of
37  the corresponding finite element space Xh: (phi_i), i=0..dim(Xh)-1
38 
39  The integrals over K are transformed on the reference element with
40  the piola transformation:
41  F : hat_K ---> K
42  hat_x |--> x = F(hat_x)
43 
44  exemples:
45  1) expr(v) = v
46  int_K phi_i(x) dx
47  = int_{hat_K} hat_phi_i(hat_x) det(DF(hat_x)) d hat_x
48  = sum_q hat_phi_i(hat_xq) det(DF(hat_xq)) hat_wq
49 
50  The value(q,i) = (hat_phi_i(hat_xq)) are evaluated on time for all over the
51  reference element hat_K and for the given quadrature formulae by:
52  expr.initialize (omega, quad);
53  This expression is represented by the 'test' class (see test.h)
54 
55  2) expr(v) = f*v
56  int_K f(x)*phi_i(x) dx
57  = int_{hat_K} f(F(hat_x)*hat_phi_i(hat_x) det(DF(hat_x)) d hat_x
58  = sum_q f(F(hat_xq))*hat_phi_i(hat_xq) det(DF(hat_xq)) hat_wq
59 
60  On each K, the computation needs two imbricated loops in (q,i).
61  The expresion is represented by a tree at compile-time.
62  The first subexpr 'f' is represented by field_vf_expr_terminal<f> : it evaluates :
63  value1(q,i) = f(F(hat_xq)) : the values depends only of q and are independent of i.
64  They could be evaluated juste before the 'i' loop.
65  As the field_vf_expr_terminal<> is general and can handle also more complex expressions
66  involving fields with various approximations, all the values on K are evaluated
67  in a specific 'q' loop by
68  subexpr1.element_initialize (K);
69  The second subexpr 'phi_i' is represdented by the 'test' class (see before).
70  value2(q,i) = (hat_phi_i(hat_xq))
71  The '*' performs on the fly the product
72  value(q,i) = value1(q,i)*value2(q,i)
73 
74  3) expr(v) = dot(f,grad(v)) dx
75  The 'f' function is here vector-valued.
76  The expresion is represented by a tree at compile-time :
77  dot -> f
78  -> grad -> v
79  The 'grad' node returns
80  value(q,i) = trans(inv(DF(hat_wq))*grad_phi_i(hat_xq) that is vector-valued
81  The grad_phi values are obtained by a grad_value(q,i) method on the 'test' class.
82  The 'dot' performs on the fly the product
83  value(q,i) = ot (value1(q,i), value2(q,i))
84 
85  This approch generilize for an expression tree.
86 */
87 namespace rheolef {
88 
89 // ====================================================================
90 // common implementation for integration on a band or an usual domain
91 // ====================================================================
92 template <class T, class M>
93 template <class Expr>
94 void
96  const geo_basic<T,M>& omega,
97  const geo_basic<T,M>& band,
98  const band_basic<T,M>& gh,
99  const Expr& expr0,
100  const integrate_option& iopt,
101  bool is_on_band)
102 {
103  Expr expr = expr0; // so could call expr.initialize()
104  // ------------------------------------
105  // 0) initialize the loop
106  // ------------------------------------
107  const space_basic<T,M>& Xh = expr.get_vf_space();
108  const geo_basic<T,M>& bgd_omega = Xh.get_constitution().get_background_geo();
109  if (is_on_band) {
110  check_macro (band.get_background_geo() == bgd_omega,
111  "do_integrate: incompatible integration domain "<<omega.name() << " and test function based domain "
112  << bgd_omega.name());
113  }
114  resize (Xh, 0);
115 
116  if (is_on_band) {
117  expr.initialize (gh, iopt);
118  } else {
119  expr.initialize (omega, iopt);
120  }
121  expr.template valued_check<T>();
122  vec<T,M>& u = set_u();
123  vec<T,M>& b = set_b();
124  std::vector<size_type> dis_idx;
125  Eigen::Matrix<T,Eigen::Dynamic,1> lk;
126  size_type d = omega.dimension();
127  size_type map_d = omega.map_dimension();
128  if (Xh.get_constitution().is_discontinuous()) Xh.get_constitution().neighbour_guard();
129  // ------------------------------------
130  // 1) loop on elements
131  // ------------------------------------
132  for (size_type ie = 0, ne = omega.size(map_d); ie < ne; ie++) {
133  const geo_element& K = omega.get_geo_element (map_d, ie);
134  // ------------------------------------
135  // 1.1) locally compute dofs
136  // ------------------------------------
137  if (!is_on_band) {
138  Xh.get_constitution().assembly_dis_idof (omega, K, dis_idx);
139  } else { // on a band
140  size_type L_ie = gh.sid_ie2bnd_ie (ie);
141  const geo_element& L = band [L_ie];
142  Xh.dis_idof (L, dis_idx);
143  }
144  // ------------------------------------
145  // 1.2) locally compute values
146  // ------------------------------------
147  expr.evaluate (omega, K, lk);
148  // -----------------------------------------
149  // 1.3) assembly local lk in global field lh
150  // -----------------------------------------
151  check_macro (dis_idx.size() == size_type(lk.size()),
152  "incompatible sizes dis_idx("<<dis_idx.size()<<") and lk("<<lk.size()<<")");
153  for (size_type loc_idof = 0, loc_ndof = lk.size(); loc_idof < loc_ndof; loc_idof++) {
154  const T& value = lk [loc_idof];
155  size_type dis_idof = dis_idx [loc_idof];
156  size_type dis_iub = Xh.dis_idof2dis_iub(dis_idof);
157  if (Xh.dis_is_blocked(dis_idof)) b.dis_entry(dis_iub) += value;
158  else u.dis_entry(dis_iub) += value;
159  }
160  }
161  // -----------------------------------------
162  // 2) finalize the assembly process
163  // -----------------------------------------
164  u.dis_entry_assembly(details::generic_set_plus_op());
165  b.dis_entry_assembly(details::generic_set_plus_op());
166 }
167 template <class T, class M>
168 template <class Expr>
169 inline
170 void
172  const geo_basic<T,M>& omega,
173  const Expr& expr,
174  const integrate_option& iopt)
175 {
176  do_integrate_internal (omega, omega, band_basic<T,M>(), expr, iopt, false);
177 }
178 template <class T, class M>
179 template <class Expr>
180 inline
181 void
183  const band_basic<T,M>& gh,
184  const Expr& expr,
185  const integrate_option& iopt)
186 {
187  do_integrate_internal (gh.level_set(), gh.band(), gh, expr, iopt, true);
188 }
189 
190 }// namespace rheolef
191 #endif // _RHEO_FIELD_VF_ASSEMBLY_H
field gh(Float epsilon, Float t, const field &uh, const test &v)
see the band page for the full documentation
void do_integrate_internal(const geo_basic< T, M > &dom, const geo_basic< T, M > &band, const band_basic< T, M > &gh, const Expr &expr, const integrate_option &qopt, bool is_on_band)
std::size_t size_type
Definition: field.h:225
see the geo_element page for the full documentation
Definition: geo_element.h:102
see the integrate_option page for the full documentation
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)")
void dis_idof(const basis_basic< T > &b, const geo_size &gs, const geo_element &K, typename std::vector< size_type >::iterator dis_idof_tab)
This file is part of Rheolef.
Definition: leveque.h:25