Rheolef  7.2
an efficient C++ finite element environment
form_lazy_convert.h
Go to the documentation of this file.
1 #ifndef _RHEO_FORM_LAZY_CONVERT_H
2 #define _RHEO_FORM_LAZY_CONVERT_H
23 // convert to form from an un-assembled lazy_form expression
24 // AUTHOR: Pierre.Saramito@imag.fr
25 // DATE: 30 march 2020
26 
27 /*
28  a(u,v) = int_omega expr dx
29  = sum_K int_K expr dx
30 
31  The integrals are evaluated over each element K in omega
32  and expr is a bilinear expression, as returned by form_lazy.
33  u & v are trial and test functions, involved in expr.
34 */
35 
36 #include "rheolef/form.h"
37 
38 namespace rheolef {
39 
40 namespace details {
41 // external utilities:
42 template <class T> T norm_max (Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>& m);
43 template <class T> bool check_is_symmetric (Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>& m, const T& tol_m_max);
44 } // namespace details
45 
46 // ====================================================================
47 // common implementation for integration on a band or an usual domain
48 // ====================================================================
49 template <class T, class M>
50 template <class Expr, class Sfinae>
51 void
53 {
54  Expr expr = expr0;
55  // ----------------------------------------
56  // 0) init assembly loop
57  // ----------------------------------------
58  _X = expr.get_trial_space();
59  _Y = expr.get_test_space();
60  geo_basic<T,M> omega = expr.get_geo();
61  expr.initialize (omega);
62  bool is_on_band = expr.is_on_band();
64  if (is_on_band) gh = expr.get_band();
65  asr<T,M> auu (_Y.iu_ownership(), _X.iu_ownership()),
66  aub (_Y.iu_ownership(), _X.ib_ownership()),
67  abu (_Y.ib_ownership(), _X.iu_ownership()),
68  abb (_Y.ib_ownership(), _X.ib_ownership());
69  std::vector<size_type> dis_idy, dis_jdx;
70  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> ak;
71  size_type map_d = omega.map_dimension();
72  if (_X.get_constitution().is_discontinuous()) _X.get_constitution().neighbour_guard();
73  if (_Y.get_constitution().is_discontinuous()) _Y.get_constitution().neighbour_guard();
74  bool is_sym = true;
75  const T eps = 1e3*std::numeric_limits<T>::epsilon();
76  for (size_type ie = 0, ne = omega.size(map_d); ie < ne; ie++) {
77  // ----------------------------------------
78  // 1) compute local form ak
79  // ----------------------------------------
80  const geo_element& K = omega.get_geo_element (map_d, ie);
81  if (! is_on_band) {
82  _X.get_constitution().assembly_dis_idof (omega, K, dis_jdx);
83  _Y.get_constitution().assembly_dis_idof (omega, K, dis_idy);
84  } else {
85  size_type L_ie = gh.sid_ie2bnd_ie (ie);
86  const geo_element& L = gh.band() [L_ie];
87  _X.dis_idof (L, dis_jdx);
88  _Y.dis_idof (L, dis_idy);
89  }
90  expr.evaluate (omega, K, ak);
91  // ----------------------------------------
92  // 2) optional local post-traitement
93  // ----------------------------------------
94  T ak_max = details::norm_max (ak);
95  T eps_ak_max = eps*ak_max;
96  if (is_sym) is_sym = details::check_is_symmetric (ak, eps_ak_max);
97  // ----------------------------------------
98  // 3) assembly local ak in global form a
99  // ----------------------------------------
100  check_macro (size_type(ak.rows()) == dis_idy.size() && size_type(ak.cols()) == dis_jdx.size(),
101  "invalid sizes ak("<<ak.rows()<<","<<ak.cols()
102  <<") with dis_idy("<<dis_idy.size()<<") and dis_jdx("<<dis_jdx.size()<<")");
103  for (size_type loc_idof = 0, ny = ak.rows(); loc_idof < ny; loc_idof++) {
104  for (size_type loc_jdof = 0, nx = ak.cols(); loc_jdof < nx; loc_jdof++) {
105 
106  const T& value = ak (loc_idof, loc_jdof);
107  // filter too small values
108  // * reason to perform:
109  // - efficient : lumped mass, structured meshes => sparsity increases
110  // * reason to avoid:
111  // - conserve the sparsity pattern, even with some zero coefs
112  // useful when dealing with solver::update_values()
113  // - also solver_pastix: assume sparsity pattern symmetry
114  // and failed when a_ij=0 (skipped) and a_ji=1e-15 (conserved) i.e. non-sym pattern
115  // note: this actual pastix wrapper limitation could be suppressed
116  if (fabs(value) < eps_ak_max) continue;
117  size_type dis_idof = dis_idy [loc_idof];
118  size_type dis_jdof = dis_jdx [loc_jdof];
119 
120  size_type dis_iub = _Y.dis_idof2dis_iub (dis_idof);
121  size_type dis_jub = _X.dis_idof2dis_iub (dis_jdof);
122 
123  if (_Y.dis_is_blocked(dis_idof))
124  if (_X.dis_is_blocked(dis_jdof)) abb.dis_entry (dis_iub, dis_jub) += value;
125  else abu.dis_entry (dis_iub, dis_jub) += value;
126  else
127  if (_X.dis_is_blocked(dis_jdof)) aub.dis_entry (dis_iub, dis_jub) += value;
128  else auu.dis_entry (dis_iub, dis_jub) += value;
129  }}
130  }
131  // ----------------------------------------
132  // 4) finalize the assembly process
133  // ----------------------------------------
134  //
135  // since all is local, axx.dis_entry_assembly() compute only axx.dis_nnz
136  //
137  auu.dis_entry_assembly();
138  aub.dis_entry_assembly();
139  abu.dis_entry_assembly();
140  abb.dis_entry_assembly();
141  //
142  // convert dynamic matrix asr to fixed-size one csr
143  //
144  _uu = csr<T,M>(auu);
145  _ub = csr<T,M>(aub);
146  _bu = csr<T,M>(abu);
147  _bb = csr<T,M>(abb);
148  //
149  // set pattern dimension to uu:
150  // => used by solvers, for efficiency: direct(d<3) or iterative(d=3)
151  //
152  _uu.set_pattern_dimension (map_d);
153  _ub.set_pattern_dimension (map_d);
154  _bu.set_pattern_dimension (map_d);
155  _bb.set_pattern_dimension (map_d);
156  //
157  // symmetry is used by solvers, for efficiency: LDL^t or LU, CG or GMRES
158  //
159  // Implementation note: cannot be set at compile time
160  // ex: expression=(eta*u)*v is structurally unsym, but numerical sym
161  // expression=(eta_h*grad(u))*(nu_h*grad(v)) is structurally sym,
162  // but numerical unsym when eta and nu are different tensors
163  // So, test it numerically, at element level:
164 #ifdef _RHEOLEF_HAVE_MPI
165  if (omega.comm().size() > 1 && is_distributed<M>::value) {
166  is_sym = mpi::all_reduce (omega.comm(), size_type(is_sym), mpi::minimum<size_type>());
167  }
168 #endif // _RHEOLEF_HAVE_MPI
169  _uu.set_symmetry (is_sym);
170  _bb.set_symmetry (is_sym);
171  // when sym, the main matrix is set definite and positive by default
172  _uu.set_definite_positive (is_sym);
173  _bb.set_definite_positive (is_sym);
174 }
175 
176 }// namespace rheolef
177 #endif // _RHEO_FORM_LAZY_CONVERT_H
field gh(Float epsilon, Float t, const field &uh, const test &v)
void dis_entry_assembly()
Definition: asr.h:112
dis_reference dis_entry(size_type dis_i, size_type dis_j)
Definition: asr.h:193
void convert_from_form_lazy(const Expr &expr)
csr< T, M >::size_type size_type
Definition: form.h:150
see the geo_element page for the full documentation
Definition: geo_element.h:102
size_t size_type
Definition: basis_get.cc:76
rheolef::std value
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)")
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 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.
Float epsilon