Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
adapt.cc
Go to the documentation of this file.
1
21#include "rheolef/adapt.h"
22#include "rheolef/form.h"
23#include "rheolef/field_wdof_sliced.h"
24#include "rheolef/field_expr.h"
25
26namespace rheolef {
27
28// externs:
29template<class T, class M>
30geo_basic<T,M> adapt_gmsh (const field_basic<T,M>& uh, const adapt_option& opts);
31template<class T, class M>
32geo_basic<T,M> adapt_bamg (const field_basic<T,M>& uh, const adapt_option& opts);
33
34// -----------------------------------------
35// proj
36// -----------------------------------------
37template<class T, class M>
38field_basic<T,M>
39proj (const field_basic<T,M>& uh, const std::string& approx = "P1")
40{
41 const space_basic<T,M>& Uh = uh.get_space();
42 if (Uh.get_approx() == approx) return uh;
43 space_basic<T,M> Vh (uh.get_geo(), approx, uh.valued());
44 form_basic<T,M> m (Vh, Vh, "lumped_mass");
45 form_basic<T,M> p (Uh, Vh, "mass");
46 solver_basic<T,M> sm (m.uu());
47 field_basic<T,M> vh (Vh, 0);
48 vh.set_u() = sm.solve((p*uh).u());
49 return vh;
50}
51// -----------------------------------------
52// smooth
53// -----------------------------------------
54template<class T, class M>
55field_basic<T,M>
56smooth (const field_basic<T,M>& uh, size_t n = 1) {
57 const space_basic<T,M>& Xh = uh.get_space();
58 check_macro (Xh.get_approx() == "P1", "smooth: expect P1 approx, but have " << Xh.get_approx());
59 form_basic<T,M> m (Xh, Xh, "mass");
60 form_basic<T,M> md (Xh, Xh, "lumped_mass");
61 solver_basic<T,M> smd (md.uu());
62 field_basic<T,M> vh = uh;
63 for (size_t i = 0; i < n; i++) {
64 vh.set_u() = smd.solve ((m*vh).u());
65 }
66 return vh;
67}
68// -----------------------------------------
69// hessian
70// -----------------------------------------
71template<class T, class M>
72field_basic<T,M>
74{
75 // assume that uh is P1 and scalar
76 const geo_basic<T,M>& omega = uh.get_geo();
77 const space_basic<T,M>& Xh = uh.get_space();
78 check_macro (Xh.valued() == "scalar", "hessian: unexpected "<<Xh.valued()<<"-valued field");
79 check_macro (Xh.get_approx() == "P1", "hessian: unexpected "<<Xh.get_approx()<<" approximation");
80 space_basic<T,M> Vh (omega, "P1", "vector");
81 form_basic<T,M> bv (Xh, Vh, "grad");
82 form_basic<T,M> mv (Vh, Vh, "lumped_mass");
83 // TODO: inv_mass optimize: lumped and by components
84 solver_basic<T,M> smv (mv.uu());
85 field_basic<T,M> gh (Vh, 0);
86 gh.set_u() = smv.solve ((bv*uh).u());
87 space_basic<T,M> Th (omega, "P1", "tensor");
88 form_basic<T,M> bt (Vh, Th, "2D");
89 form_basic<T,M> mt (Th, Th, "lumped_mass");
90 solver_basic<T,M> smt (mt.uu());
91 field_basic<T,M> hh (Th, 0);
92 hh.set_u() = smt.solve ((bt*gh).u());
93 return hh;
94}
95// |hessian(uh)|
96// M = ----------------------------
97// err*hcoef^2*(sup(uh)-inf(uh))
98// where
99// |H| = same as H but with absolute value of eigenvalues
100// = Q*diag(|lambda_i|)*Qt
101//
102template<class T, class M>
103field_basic<T,M>
105 const field_basic<T,M>& uh0,
106 const adapt_option& opts)
107{
108 typedef typename geo_basic<T,M>::size_type size_type;
109 size_type d = uh0.get_geo().dimension();
110 size_type k = uh0.get_space().degree();
111 field_basic<T,M> uh = proj(uh0);
112 field_basic<T,M> hh = hessian(uh);
113 field_basic<T,M> mh (hh.get_space(), 0);
114 field_basic<T,M> sh (uh.get_space(), 0);
115#ifdef TO_CLEAN
117 for (size_type i_comp = 0; i_comp < d; i_comp++) {
118 for (size_type j_comp = 0; j_comp < d; j_comp++) {
119 hh_comp[i_comp][j_comp].proxy_assign (hh(i_comp,j_comp));
120 }}
121#endif // TO_CLEAN
122 tensor_basic<T> h_value, m_value, Q, Qt;
124 T cut_off = 1e-5;
125 T uh_scale = std::max(cut_off, fabs(uh.max() - uh.min()));
126 T factor = opts.hcoef*sqrt(uh_scale)*pow(opts.err,1./(k+1));
127 T eps = std::numeric_limits<T>::epsilon();
128 for (size_type idof = 0, ndof = uh.ndof(); idof < ndof; idof++) {
129 for (size_type i_comp = 0; i_comp < d; i_comp++) {
130 for (size_type j_comp = 0; j_comp < d; j_comp++) {
131 h_value(i_comp,j_comp) = hh(i_comp,j_comp).dof (idof);
132 }}
133 const bool use_svd_when_2d = true;
134 // H = Q*diag(lambda)*Q^T
135 if (use_svd_when_2d && d == 2) {
136 // TODO: eig when d=2 is bad...
137 lambda = h_value.svd (Q, Qt, d);
138 } else {
139 lambda = h_value.eig (Q, d);
140 }
141 T h_min = opts.hmax;
142 for (size_type i_comp = 0; i_comp < d; i_comp++) {
143 if (fabs(lambda[i_comp]) < eps) continue;
144 // err = c*h^2
145 h_local[i_comp] = factor/sqrt(fabs(lambda[i_comp]));
146 // yield values
147 h_local[i_comp] = std::min (opts.hmax, h_local[i_comp]);
148 h_local[i_comp] = std::max (opts.hmin, h_local[i_comp]);
149 h_min = std::min (h_min, h_local[i_comp]);
150 }
151 // isotropic: takes h_local = h_min in all directions
152 sh.dof (idof) = h_min;
153
154 // anisotropic
155 m_value = Q*diag(h_local)*trans(Q);
156 for (size_type i_comp = 0; i_comp < d; i_comp++) {
157 for (size_type j_comp = 0; j_comp < d; j_comp++) {
158 mh(i_comp,j_comp).dof (idof) = m_value(i_comp,j_comp);
159 }}
160 }
161#ifdef TO_CLEAN
162 trace_macro ("sh.min="<<sh.min()<<", sh.max="<<sh.max());
163 mh *= factor;
164 sh *= factor;
165 trace_macro ("uh_scale="<<uh_scale<<", factor="<<uh_scale<<" => sh.min="<<sh.min()<<", sh.max="<<sh.max());
166#endif // TO_CLEAN
167 if (opts.isotropic) {
168 return smooth (sh, opts.n_smooth_metric);
169 } else {
170 return smooth (mh, opts.n_smooth_metric);
171 }
172}
173// -----------------------------------------
174// adapt
175// -----------------------------------------
177template<class T, class M>
178geo_basic<T,M>
180 const field_basic<T,M>& uh,
181 const adapt_option& opts)
182{
183 size_t d = uh.get_geo().dimension();
184 if (d == 2 && (opts.generator == "bamg" || opts.generator == "")) {
185 // default generator is bamg when d=2:
186 return adapt_bamg (uh, opts);
187 } else {
188 // use always gmsh when d != 2:
189 return adapt_gmsh (uh, opts);
190 }
191}
192// ----------------------------------------------------------------------------
193// instanciation in library
194// ----------------------------------------------------------------------------
195#define _RHEOLEF_instanciation(T,M) \
196template field_basic<T,M> proj (const field_basic<T,M>&, const std::string&); \
197template field_basic<T,M> smooth (const field_basic<T,M>&, size_t); \
198template field_basic<T,M> hessian (const field_basic<T,M>&); \
199template field_basic<T,M> hessian_criterion (const field_basic<T,M>&, const adapt_option&); \
200template geo_basic<T,M> adapt (const field_basic<T,M>&, const adapt_option&);
201
203#ifdef _RHEOLEF_HAVE_MPI
204_RHEOLEF_instanciation(Float,distributed)
205#endif // _RHEOLEF_HAVE_MPI
206
207} // namespace rheolef
#define _RHEOLEF_instanciation(T, M, A)
Definition asr.cc:223
field::size_type size_type
Definition branch.cc:430
field gh(Float epsilon, Float t, const field &uh, const test &v)
see the Float page for the full documentation
T & dof(size_type idof)
Definition field.h:738
vec< T, M > & set_u()
Definition field.h:284
T min() const
Definition field.h:786
T max() const
Definition field.h:802
const geo_type & get_geo() const
Definition field.h:271
const space_type & get_space() const
Definition field.h:270
size_type ndof() const
Definition field.h:298
const std::string & valued() const
Definition field.h:274
const csr< T, M > & uu() const
Definition form.h:204
generic mesh with rerefence counting
Definition geo.h:1089
vec< T, M > solve(const vec< T, M > &b) const
Definition solver.h:345
the finite element space
Definition space.h:382
point_basic< T > eig(tensor_basic< T > &q, size_t dim=3) const
Definition tensor.cc:426
point_basic< T > svd(tensor_basic< T > &u, tensor_basic< T > &v, size_t dim=3) const
Definition tensor.cc:470
#define trace_macro(message)
Definition dis_macros.h:111
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)")
This file is part of Rheolef.
field_basic< T, M > hessian_criterion(const field_basic< T, M > &uh0, const adapt_option &opts)
Definition adapt.cc:104
geo_basic< T, M > adapt(const field_basic< T, M > &uh, const adapt_option &opts)
adapt(uh,opts): see the adapt page for the full documentation
Definition adapt.cc:179
field_basic< T, M > hessian(const field_basic< T, M > &uh)
Definition adapt.cc:73
space_mult_list< T, M > pow(const space_basic< T, M > &X, size_t n)
Definition space_mult.h:120
field_basic< T, M > proj(const field_basic< T, M > &uh, const std::string &approx="P1")
Definition adapt.cc:39
geo_basic< T, M > adapt_gmsh(const field_basic< T, M > &uh, const adapt_option &opts)
Definition adapt_gmsh.cc:41
field_basic< T, M > smooth(const field_basic< T, M > &uh, size_t n=1)
Definition adapt.cc:56
details::field_expr_v2_nonlinear_terminal_function< details::h_local_pseudo_function< Float > > h_local()
h_local: see the expression page for the full documentation
geo_basic< T, M > adapt_bamg(const field_basic< T, M > &uh, const adapt_option &opts)
Definition adapt_bamg.cc:63
csr< T, M > diag(const vec< T, M > &d)
Definition csr.cc:56
csr< T, sequential > trans(const csr< T, sequential > &a)
trans(a): see the form page for the full documentation
Definition csr.h:455
Definition sphere.icc:25
adapt_option: see the adapt page for the full documentation
Definition adapt.h:147
std::string generator
Definition adapt.h:149
size_type n_smooth_metric
Definition adapt.h:159
Definition leveque.h:25