Rheolef  7.1
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
3 #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& expr,
100  const integrate_option& iopt,
101  bool is_on_band)
102 {
103  // ------------------------------------
104  // 0) initialize the loop
105  // ------------------------------------
106  const space_basic<T,M>& Xh = expr.get_vf_space();
107  const geo_basic<T,M>& bgd_omega = Xh.get_constitution().get_background_geo();
108  if (is_on_band) {
109  check_macro (band.get_background_geo() == bgd_omega,
110  "assembly: incompatible integration domain "<<omega.name() << " and test function based domain "
111  << bgd_omega.name());
112  }
113  resize (Xh, 0);
114 
115  if (is_on_band) {
116  expr.initialize (gh, iopt);
117  } else {
118  expr.initialize (omega, iopt);
119  }
120  expr.template valued_check<T>();
121  vec<T,M>& u = set_u();
122  vec<T,M>& b = set_b();
123  std::vector<size_type> dis_idx;
124  Eigen::Matrix<T,Eigen::Dynamic,1> lk;
125  size_type d = omega.dimension();
126  size_type map_d = omega.map_dimension();
127  if (Xh.get_constitution().is_discontinuous()) Xh.get_constitution().neighbour_guard();
128  // ------------------------------------
129  // 1) loop on elements
130  // ------------------------------------
131  for (size_type ie = 0, ne = omega.size(map_d); ie < ne; ie++) {
132  const geo_element& K = omega.get_geo_element (map_d, ie);
133  // ------------------------------------
134  // 1.1) locally compute dofs
135  // ------------------------------------
136  if (!is_on_band) {
137  Xh.get_constitution().assembly_dis_idof (omega, K, dis_idx);
138  } else { // on a band
139  size_type L_ie = gh.sid_ie2bnd_ie (ie);
140  const geo_element& L = band [L_ie];
141  Xh.dis_idof (L, dis_idx);
142  }
143  // ------------------------------------
144  // 1.2) locally compute values
145  // ------------------------------------
146  expr.evaluate (omega, K, lk);
147  // -----------------------------------------
148  // 1.3) assembly local lk in global field lh
149  // -----------------------------------------
150  check_macro (dis_idx.size() == size_type(lk.size()),
151  "incompatible sizes dis_idx("<<dis_idx.size()<<") and lk("<<lk.size()<<")");
152  for (size_type loc_idof = 0, loc_ndof = lk.size(); loc_idof < loc_ndof; loc_idof++) {
153  const T& value = lk [loc_idof];
154  size_type dis_idof = dis_idx [loc_idof];
155  size_type dis_iub = Xh.dis_idof2dis_iub(dis_idof);
156  if (Xh.dis_is_blocked(dis_idof)) b.dis_entry(dis_iub) += value;
157  else u.dis_entry(dis_iub) += value;
158  }
159  }
160  // -----------------------------------------
161  // 2) finalize the assembly process
162  // -----------------------------------------
163  u.dis_entry_assembly(set_add_op<T,T>());
164  b.dis_entry_assembly(set_add_op<T,T>());
165 }
166 template <class T, class M>
167 template <class Expr>
168 inline
169 void
171  const geo_basic<T,M>& omega,
172  const Expr& expr,
173  const integrate_option& iopt)
174 {
175  assembly_internal (omega, omega, band_basic<T,M>(), expr, iopt, false);
176 }
177 template <class T, class M>
178 template <class Expr>
179 inline
180 void
182  const band_basic<T,M>& gh,
183  const Expr& expr,
184  const integrate_option& iopt)
185 {
186  assembly_internal (gh.level_set(), gh.band(), gh, expr, iopt, true);
187 }
188 
189 }// namespace rheolef
190 #endif // _RHEO_FIELD_VF_ASSEMBLY_H
rheolef::geo_basic
generic mesh with rerefence counting
Definition: geo.h:1089
rheolef::field_basic::assembly
void assembly(const geo_basic< T, M > &dom, const Expr &expr, const integrate_option &iopt)
Definition: field_vf_assembly.h:170
gh
field gh(Float epsilon, Float t, const field &uh, const test &v)
Definition: burgers_diffusion_operators.icc:37
mkgeo_ball.expr
expr
Definition: mkgeo_ball.sh:361
check_macro
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
mkgeo_ball.b
int b
Definition: mkgeo_ball.sh:152
rheolef::value
rheolef::std value
rheolef::space_basic
the finite element space
Definition: space.h:352
band
see the band page for the full documentation
mkgeo_ball.d
int d
Definition: mkgeo_ball.sh:154
rheolef::geo_element
see the geo_element page for the full documentation
Definition: geo_element.h:102
rheolef::space_numbering::dis_idof
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)
Definition: space_numbering.cc:187
rheolef::vec
see the vec page for the full documentation
Definition: vec.h:79
rheolef::integrate_option
see the integrate_option page for the full documentation
Definition: integrate_option.h:125
rheolef::band_basic
Definition: band.h:95
rheolef
This file is part of Rheolef.
Definition: compiler_eigen.h:37
u
Definition: leveque.h:25
mkgeo_obstacle.L
string L
Definition: mkgeo_obstacle.sh:133
rheolef::field_basic::assembly_internal
void assembly_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)
Definition: field_vf_assembly.h:95
rheolef::set_add_op
Definition: msg_util.h:61
rheolef::field_basic::size_type
std::size_t size_type
Definition: field.h:239
T
Expr1::float_type T
Definition: field_expr.h:261