Rheolef  7.2
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 namespace details {
77 template <class T> T norm_max (Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>& m);
78 template <class T> bool check_is_symmetric (Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>& m, const T& tol_m_max);
79 template <class T> void local_lump (Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>& m);
80 template <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 // ====================================================================
86 template <class T, class M>
87 template <class Expr>
88 void
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 
240 template <class T, class M>
241 template <class Expr>
242 inline
243 void
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 }
251 template <class T, class M>
252 template <class Expr>
253 inline
254 void
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 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
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)
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: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)
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