Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
form_vf_assembly.h
Go to the documentation of this file.
1#ifndef _RHEO_FORM_VF_ASSEMBLY_H
2#define _RHEO_FORM_VF_ASSEMBLY_H
23#include "rheolef/form.h"
24#include "rheolef/test.h"
25#include "rheolef/quadrature.h"
26#include "rheolef/field_vf_assembly.h" // for dg_dis_idof()
27#include "rheolef/form_expr_quadrature.h"
28namespace rheolef {
29/*
30 let:
31 a(u,v) = int_domain expr(u,v) dx
32
33 The integrals are evaluated over each element K of the domain
34 by using a quadrature formulae given by iopt
35
36 expr(u,v) is a bilinear expression with respect to the
37 trial and test functions u and v
38
39 The trial function u is replaced by each of the basis function of
40 the corresponding finite element space Xh: (phi_j), j=0..dim(Xh)-1
41
42 The test function v is replaced by each of the basis function of
43 the corresponding finite element space Yh: (psi_i), i=0..dim(Yh)-1
44
45 The integrals over the domain omega is the sum of integrals over K.
46
47 The integrals over K are transformed on the reference element with
48 the piola transformation:
49 F : hat_K ---> K
50 hat_x |--> x = F(hat_x)
51
52 exemples:
53 1) expr(v) = u*v
54 int_K phi_j(x)*psi_i(x) dx
55 = int_{hat_K} hat_phi_j(hat_x)*hat_psi_i(hat_x) det(DF(hat_x)) d hat_x
56 = sum_q hat_phi_j(hat_xq)*hat_psi_i(hat_xq) det(DF(hat_xq)) hat_wq
57
58 The value(q,i,j) = (hat_phi_j(hat_xq)*hat_psi_i(hat_xq))
59 refers to basis values on the reference element.
60 There are evaluated on time for all over the reference element hat_K
61 and for the given quadrature formulae by:
62 expr.initialize (omega, quad);
63 This expression is represented by the 'test' class (see test.h)
64
65 3) expr(v) = dot(grad(u),grad(v)) dx
66 The 'grad' node returns
67 value(q,i) = trans(inv(DF(hat_wq))*grad_phi_i(hat_xq) that is vector-valued
68 The grad_phi values are obtained by a grad_value(q,i) method on the 'test' class.
69 The 'dot' performs on the fly the product
70 value(q,i,j) = dot (value1(q,i), value2(q,j))
71
72 This approch generalizes for an expression tree.
73*/
74
75// external utilities:
76namespace details {
77template <class T> T norm_max (Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>& m);
78template <class T> bool check_is_symmetric (Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>& m, const T& tol_m_max);
79template <class T> void local_lump (Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>& m);
80template <class T> void local_invert (Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>& m, bool is_diag);
81} // namespace details
82
83// ====================================================================
84// common implementation for integration on a band or an usual domain
85// ====================================================================
86template <class T, class M>
87template <class Expr>
88void
90 const geo_basic<T,M>& omega,
91 const geo_basic<T,M>& band,
92 const band_basic<T,M>& gh,
93 const Expr& expr0,
94 const integrate_option& iopt,
95 bool is_on_band)
96{
97 Expr expr = expr0; // so could call expr.initialize()
98 // ----------------------------------------
99 // 0) init assembly loop
100 // ----------------------------------------
101 _X = expr.get_trial_space();
102 _Y = expr.get_test_space();
103 if (is_on_band) {
104 check_macro (band.get_background_geo() == _X.get_geo().get_background_geo(),
105 "do_integratessembly: incompatible integration domain "<<band.name() << " and trial function based domain "
106 << _X.get_geo().name());
107 check_macro (band.get_background_geo() == _Y.get_geo().get_background_geo(),
108 "do_integratessembly: incompatible integration domain "<<band.name() << " and test function based domain "
109 << _Y.get_geo().name());
110 }
111 size_type n_derivative = expr.n_derivative();
112
113 if (iopt.invert) check_macro (
114 _X.get_constitution().have_compact_support_inside_element() &&
115 _Y.get_constitution().have_compact_support_inside_element(),
116 "local inversion requires compact support inside elements (e.g. discontinuous or bubble)");
117 if (iopt.lump) check_macro (n_derivative == 0,
118 "local mass lumping requires no derivative operators");
119
120 iopt._is_on_interface = false;
121 iopt._is_inside_on_local_sides = false;
122 if (!is_on_band) {
123 expr.initialize (omega, iopt);
124 } else {
125 expr.initialize (gh, iopt);
126 }
127 expr.template valued_check<T>();
128 asr<T,M> auu (_Y.iu_ownership(), _X.iu_ownership()),
129 aub (_Y.iu_ownership(), _X.ib_ownership()),
130 abu (_Y.ib_ownership(), _X.iu_ownership()),
131 abb (_Y.ib_ownership(), _X.ib_ownership());
132 std::vector<size_type> dis_idy, dis_jdx;
133 Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> ak;
134 size_type map_d = omega.map_dimension();
135 if (_X.get_constitution().is_discontinuous()) _X.get_constitution().neighbour_guard();
136 if (_Y.get_constitution().is_discontinuous()) _Y.get_constitution().neighbour_guard();
137 bool is_sym = true;
138 const T eps = 1e3*std::numeric_limits<T>::epsilon();
139 for (size_type ie = 0, ne = omega.size(map_d); ie < ne; ie++) {
140 // ----------------------------------------
141 // 1) compute local form ak
142 // ----------------------------------------
143 const geo_element& K = omega.get_geo_element (map_d, ie);
144 if (! is_on_band) {
145 _X.get_constitution().assembly_dis_idof (omega, K, dis_jdx);
146 _Y.get_constitution().assembly_dis_idof (omega, K, dis_idy);
147 } else {
148 size_type L_ie = gh.sid_ie2bnd_ie (ie);
149 const geo_element& L = band [L_ie];
150 _X.dis_idof (L, dis_jdx);
151 _Y.dis_idof (L, dis_idy);
152 }
153 expr.evaluate (omega, K, ak);
154 // ----------------------------------------
155 // 2) optional local post-traitement
156 // ----------------------------------------
157 T ak_max = details::norm_max (ak);
158 T eps_ak_max = eps*ak_max;
159 if (is_sym) is_sym = details::check_is_symmetric (ak, eps_ak_max);
160 if (iopt.lump ) details::local_lump (ak);
161 if (iopt.invert) details::local_invert (ak, iopt.lump);
162 // ----------------------------------------
163 // 3) assembly local ak in global form a
164 // ----------------------------------------
165 check_macro (size_type(ak.rows()) == dis_idy.size() && size_type(ak.cols()) == dis_jdx.size(),
166 "invalid sizes ak("<<ak.rows()<<","<<ak.cols()
167 <<") with dis_idy("<<dis_idy.size()<<") and dis_jdx("<<dis_jdx.size()<<")");
168 for (size_type loc_idof = 0, ny = ak.rows(); loc_idof < ny; loc_idof++) {
169 for (size_type loc_jdof = 0, nx = ak.cols(); loc_jdof < nx; loc_jdof++) {
170
171 const T& value = ak (loc_idof, loc_jdof);
172 // reason to perform:
173 // - efficient : lumped mass, structured meshes => sparsity increases
174 // reason to avoid:
175 // - conserve the sparsity pattern, even with some zero coefs
176 // useful when dealing with solver::update_values()
177 // - also solver_pastix: assume sparsity pattern symmetry
178 // and failed when a_ij=0 (skipped) and a_ji=1e-15 (conserved) i.e. non-sym pattern
179 // note: this actual pastix wrapper limitation could be suppressed
180 if (fabs(value) < eps_ak_max) continue;
181 size_type dis_idof = dis_idy [loc_idof];
182 size_type dis_jdof = dis_jdx [loc_jdof];
183
184 size_type dis_iub = _Y.dis_idof2dis_iub (dis_idof);
185 size_type dis_jub = _X.dis_idof2dis_iub (dis_jdof);
186
187 if (_Y.dis_is_blocked(dis_idof))
188 if (_X.dis_is_blocked(dis_jdof)) abb.dis_entry (dis_iub, dis_jub) += value;
189 else abu.dis_entry (dis_iub, dis_jub) += value;
190 else
191 if (_X.dis_is_blocked(dis_jdof)) aub.dis_entry (dis_iub, dis_jub) += value;
192 else auu.dis_entry (dis_iub, dis_jub) += value;
193 }}
194 }
195 // ----------------------------------------
196 // 4) finalize the assembly process
197 // ----------------------------------------
198 //
199 // since all is local, axx.dis_entry_assembly() compute only axx.dis_nnz
200 //
201 auu.dis_entry_assembly();
202 aub.dis_entry_assembly();
203 abu.dis_entry_assembly();
204 abb.dis_entry_assembly();
205 //
206 // convert dynamic matrix asr to fixed-size one csr
207 //
208 _uu = csr<T,M>(auu);
209 _ub = csr<T,M>(aub);
210 _bu = csr<T,M>(abu);
211 _bb = csr<T,M>(abb);
212 //
213 // set pattern dimension to uu:
214 // => used by solvers, for efficiency: direct(d<3) or iterative(d=3)
215 //
216 _uu.set_pattern_dimension (map_d);
217 _ub.set_pattern_dimension (map_d);
218 _bu.set_pattern_dimension (map_d);
219 _bb.set_pattern_dimension (map_d);
220 //
221 // symmetry is used by solvers, for efficiency: LDL^t or LU, CG or GMRES
222 //
223 // Implementation note: cannot be set at compile time
224 // ex: expression=(eta*u)*v is structurally unsym, but numerical sym
225 // expression=(eta_h*grad(u))*(nu_h*grad(v)) is structurally sym,
226 // but numerical unsym when eta and nu are different tensors
227 // So, test it numerically, at element level:
228#ifdef _RHEOLEF_HAVE_MPI
229 if (omega.comm().size() > 1 && is_distributed<M>::value) {
230 is_sym = mpi::all_reduce (omega.comm(), size_type(is_sym), mpi::minimum<size_type>());
231 }
232#endif // _RHEOLEF_HAVE_MPI
233 _uu.set_symmetry (is_sym);
234 _bb.set_symmetry (is_sym);
235 // when sym, the main matrix is set definite and positive by default
236 _uu.set_definite_positive (is_sym);
237 _bb.set_definite_positive (is_sym);
238}
239
240template <class T, class M>
241template <class Expr>
242inline
243void
245 const geo_basic<T,M>& omega,
246 const Expr& expr,
247 const integrate_option& iopt)
248{
249 do_integrate_internal (omega, omega, band_basic<T,M>(), expr, iopt, false);
250}
251template <class T, class M>
252template <class Expr>
253inline
254void
256 const band_basic<T,M>& gh,
257 const Expr& expr,
258 const integrate_option& iopt)
259{
260 do_integrate_internal (gh.level_set(), gh.band(), gh, expr, iopt, true);
261}
262
263}// namespace rheolef
264#endif // _RHEO_FORM_VF_ASSEMBLY_H
field::size_type size_type
Definition branch.cc:430
field gh(Float epsilon, Float t, const field &uh, const test &v)
see the band page for the full documentation
void dis_entry_assembly()
Definition asr.h:112
dis_reference dis_entry(size_type dis_i, size_type dis_j)
Definition asr.h:193
see the csr page for the full documentation
Definition csr.h:317
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 &fopt, bool is_on_band)
csr< T, M >::size_type size_type
Definition form.h:150
void do_integrate(const geo_basic< T, M > &domain, const Expr &expr, const integrate_option &fopt)
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
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 local_lump(Eigen::Matrix< T, Eigen::Dynamic, Eigen::Dynamic > &m)
bool check_is_symmetric(Eigen::Matrix< T, Eigen::Dynamic, Eigen::Dynamic > &m, const T &tol_m_max)
T norm_max(Eigen::Matrix< T, Eigen::Dynamic, Eigen::Dynamic > &m)
void local_invert(Eigen::Matrix< T, Eigen::Dynamic, Eigen::Dynamic > &m, bool is_diag)
This file is part of Rheolef.