Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
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
38namespace rheolef {
39
40namespace details {
41// external utilities:
42template <class T> T norm_max (Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>& m);
43template <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// ====================================================================
49template <class T, class M>
50template <class Expr, class Sfinae>
51void
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::size_type size_type
Definition branch.cc:430
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
see the csr page for the full documentation
Definition csr.h:317
void convert_from_form_lazy(const Expr &expr)
csr< T, M >::size_type size_type
Definition form.h:150
generic mesh with rerefence counting
Definition geo.h:1089
see the geo_element 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)")
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)
This file is part of Rheolef.