Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
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*/
87namespace rheolef {
88
89// ====================================================================
90// common implementation for integration on a band or an usual domain
91// ====================================================================
92template <class T, class M>
93template <class Expr>
94void
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}
167template <class T, class M>
168template <class Expr>
169inline
170void
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}
178template <class T, class M>
179template <class Expr>
180inline
181void
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)
void do_integrate(const geo_basic< T, M > &dom, const Expr &expr, const integrate_option &iopt)
std::size_t size_type
Definition field.h:225
generic mesh with rerefence counting
Definition geo.h:1089
see the geo_element page for the full documentation
see the integrate_option page for the full documentation
the finite element space
Definition space.h:382
see the vec page for the full documentation
Definition vec.h:79
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)")
This file is part of Rheolef.
Definition leveque.h:25