Rheolef  7.1
an efficient C++ finite element environment
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"
28 namespace 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:
76 template <class T> T norm_max (Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>& m);
77 template <class T> bool check_is_symmetric (Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>& m, const T& tol_m_max);
78 template <class T> void local_lump (Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>& m);
79 template <class T> void local_invert (Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>& m, bool is_diag);
80 
81 // ====================================================================
82 // common implementation for integration on a band or an usual domain
83 // ====================================================================
84 template <class T, class M>
85 template <class Expr>
86 void
88  const geo_basic<T,M>& omega,
89  const geo_basic<T,M>& band,
90  const band_basic<T,M>& gh,
91  const Expr& expr,
92  const integrate_option& iopt,
93  bool is_on_band)
94 {
95  // ----------------------------------------
96  // 0) init assembly loop
97  // ----------------------------------------
98  _X = expr.get_trial_space();
99  _Y = expr.get_test_space();
100  if (is_on_band) {
101  check_macro (band.get_background_geo() == _X.get_geo().get_background_geo(),
102  "assembly: incompatible integration domain "<<band.name() << " and trial function based domain "
103  << _X.get_geo().name());
104  check_macro (band.get_background_geo() == _Y.get_geo().get_background_geo(),
105  "assembly: incompatible integration domain "<<band.name() << " and test function based domain "
106  << _Y.get_geo().name());
107  }
108  size_type n_derivative = expr.n_derivative();
109 
110  if (iopt.invert) check_macro (
111  _X.get_constitution().have_compact_support_inside_element() &&
112  _Y.get_constitution().have_compact_support_inside_element(),
113  "local inversion requires compact support inside elements (e.g. discontinuous or bubble)");
114  if (iopt.lump) check_macro (n_derivative == 0,
115  "local mass lumping requires no derivative operators");
116 
117  iopt._is_on_interface = false;
118  iopt._is_inside_on_local_sides = false;
119  if (!is_on_band) {
120  expr.initialize (omega, iopt);
121  } else {
122  expr.initialize (gh, iopt);
123  }
124  expr.template valued_check<T>();
125  asr<T,M> auu (_Y.iu_ownership(), _X.iu_ownership()),
126  aub (_Y.iu_ownership(), _X.ib_ownership()),
127  abu (_Y.ib_ownership(), _X.iu_ownership()),
128  abb (_Y.ib_ownership(), _X.ib_ownership());
129  std::vector<size_type> dis_idy, dis_jdx;
130  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> ak;
131  size_type map_d = omega.map_dimension();
132  if (_X.get_constitution().is_discontinuous()) _X.get_constitution().neighbour_guard();
133  if (_Y.get_constitution().is_discontinuous()) _Y.get_constitution().neighbour_guard();
134  bool is_sym = true;
135  const T eps = 1e3*std::numeric_limits<T>::epsilon();
136  for (size_type ie = 0, ne = omega.size(map_d); ie < ne; ie++) {
137  // ----------------------------------------
138  // 1) compute local form ak
139  // ----------------------------------------
140  const geo_element& K = omega.get_geo_element (map_d, ie);
141  if (! is_on_band) {
142  _X.get_constitution().assembly_dis_idof (omega, K, dis_jdx);
143  _Y.get_constitution().assembly_dis_idof (omega, K, dis_idy);
144  } else {
145  size_type L_ie = gh.sid_ie2bnd_ie (ie);
146  const geo_element& L = band [L_ie];
147  _X.dis_idof (L, dis_jdx);
148  _Y.dis_idof (L, dis_idy);
149  }
150  expr.evaluate (omega, K, ak);
151  // ----------------------------------------
152  // 2) optional local post-traitement
153  // ----------------------------------------
154  T ak_max = norm_max (ak);
155  T eps_ak_max = eps*ak_max;
156  if (is_sym) is_sym = check_is_symmetric (ak, eps_ak_max);
157  if (iopt.lump ) local_lump (ak);
158  if (iopt.invert) local_invert (ak, iopt.lump);
159  // ----------------------------------------
160  // 3) assembly local ak in global form a
161  // ----------------------------------------
162  check_macro (size_type(ak.rows()) == dis_idy.size() && size_type(ak.cols()) == dis_jdx.size(),
163  "invalid sizes ak("<<ak.rows()<<","<<ak.cols()
164  <<") with dis_idy("<<dis_idy.size()<<") and dis_jdx("<<dis_jdx.size()<<")");
165  for (size_type loc_idof = 0, ny = ak.rows(); loc_idof < ny; loc_idof++) {
166  for (size_type loc_jdof = 0, nx = ak.cols(); loc_jdof < nx; loc_jdof++) {
167 
168  const T& value = ak (loc_idof, loc_jdof);
169  // reason to perform:
170  // - efficient : lumped mass, structured meshes => sparsity increases
171  // reason to avoid:
172  // - conserve the sparsity pattern, even with some zero coefs
173  // useful when dealing with solver::update_values()
174  // - also solver_pastix: assume sparsity pattern symmetry
175  // and failed when a_ij=0 (skipped) and a_ji=1e-15 (conserved) i.e. non-sym pattern
176  // note: this actual pastix wrapper limitation could be suppressed
177  if (fabs(value) < eps_ak_max) continue;
178  size_type dis_idof = dis_idy [loc_idof];
179  size_type dis_jdof = dis_jdx [loc_jdof];
180 
181  size_type dis_iub = _Y.dis_idof2dis_iub (dis_idof);
182  size_type dis_jub = _X.dis_idof2dis_iub (dis_jdof);
183 
184  if (_Y.dis_is_blocked(dis_idof))
185  if (_X.dis_is_blocked(dis_jdof)) abb.dis_entry (dis_iub, dis_jub) += value;
186  else abu.dis_entry (dis_iub, dis_jub) += value;
187  else
188  if (_X.dis_is_blocked(dis_jdof)) aub.dis_entry (dis_iub, dis_jub) += value;
189  else auu.dis_entry (dis_iub, dis_jub) += value;
190  }}
191  }
192  // ----------------------------------------
193  // 4) finalize the assembly process
194  // ----------------------------------------
195  //
196  // since all is local, axx.dis_entry_assembly() compute only axx.dis_nnz
197  //
198  auu.dis_entry_assembly();
199  aub.dis_entry_assembly();
200  abu.dis_entry_assembly();
201  abb.dis_entry_assembly();
202  //
203  // convert dynamic matrix asr to fixed-size one csr
204  //
205  _uu = csr<T,M>(auu);
206  _ub = csr<T,M>(aub);
207  _bu = csr<T,M>(abu);
208  _bb = csr<T,M>(abb);
209  //
210  // set pattern dimension to uu:
211  // => used by solvers, for efficiency: direct(d<3) or iterative(d=3)
212  //
213  _uu.set_pattern_dimension (map_d);
214  _ub.set_pattern_dimension (map_d);
215  _bu.set_pattern_dimension (map_d);
216  _bb.set_pattern_dimension (map_d);
217  //
218  // symmetry is used by solvers, for efficiency: LDL^t or LU, CG or GMRES
219  //
220  // Implementation note: cannot be set at compile time
221  // ex: expression=(eta*u)*v is structurally unsym, but numerical sym
222  // expression=(eta_h*grad(u))*(nu_h*grad(v)) is structurally sym,
223  // but numerical unsym when eta and nu are different tensors
224  // So, test it numerically, at element level:
225 #ifdef _RHEOLEF_HAVE_MPI
226  if (omega.comm().size() > 1 && is_distributed<M>::value) {
227  is_sym = mpi::all_reduce (omega.comm(), size_type(is_sym), mpi::minimum<size_type>());
228  }
229 #endif // _RHEOLEF_HAVE_MPI
230  _uu.set_symmetry (is_sym);
231  _bb.set_symmetry (is_sym);
232  // when sym, the main matrix is set definite and positive by default
233  _uu.set_definite_positive (is_sym);
234  _bb.set_definite_positive (is_sym);
235 }
236 
237 template <class T, class M>
238 template <class Expr>
239 inline
240 void
242  const geo_basic<T,M>& omega,
243  const Expr& expr,
244  const integrate_option& iopt)
245 {
246  assembly_internal (omega, omega, band_basic<T,M>(), expr, iopt, false);
247 }
248 template <class T, class M>
249 template <class Expr>
250 inline
251 void
253  const band_basic<T,M>& gh,
254  const Expr& expr,
255  const integrate_option& iopt)
256 {
257  assembly_internal (gh.level_set(), gh.band(), gh, expr, iopt, true);
258 }
259 
260 }// namespace rheolef
261 #endif // _RHEO_FORM_VF_ASSEMBLY_H
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
csr< T, M >::size_type size_type
Definition: form.h:148
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 &fopt, bool is_on_band)
void assembly(const geo_basic< T, M > &domain, const Expr &expr, const integrate_option &fopt)
see the geo_element page for the full documentation
Definition: geo_element.h:102
see the integrate_option page for the full documentation
size_t size_type
Definition: basis_get.cc:76
rheolef::std value
Expr1::float_type T
Definition: field_expr.h:261
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.
void local_lump(Eigen::Matrix< T, Eigen::Dynamic, Eigen::Dynamic > &m)
void local_invert(Eigen::Matrix< T, Eigen::Dynamic, Eigen::Dynamic > &m, bool is_diag)
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)
Float epsilon