Rheolef  7.1
an efficient C++ finite element environment
adapt.cc
Go to the documentation of this file.
1 #include "rheolef/adapt.h"
22 #include "rheolef/form.h"
23 #include "rheolef/field_component.h"
24 #include "rheolef/field_expr.h"
25 
26 namespace rheolef {
27 
28 // externs:
29 template<class T, class M>
30 geo_basic<T,M> adapt_gmsh (const field_basic<T,M>& uh, const adapt_option& opts);
31 template<class T, class M>
32 geo_basic<T,M> adapt_bamg (const field_basic<T,M>& uh, const adapt_option& opts);
33 
34 // -----------------------------------------
35 // proj
36 // -----------------------------------------
37 template<class T, class M>
38 field_basic<T,M>
39 proj (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 // -----------------------------------------
54 template<class T, class M>
55 field_basic<T,M>
56 smooth (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 // -----------------------------------------
71 template<class T, class M>
72 field_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 //
102 template<class T, class M>
103 field_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  field_component_const<T,M> hh_comp [3][3];
116  field_component<T,M> mh_comp [3][3];
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  mh_comp[i_comp][j_comp].proxy_assign (mh(i_comp,j_comp));
121  }}
122  tensor_basic<T> h_value, m_value, Q, Qt;
124  for (size_type idof = 0, ndof = uh.ndof(); idof < ndof; idof++) {
125  for (size_type i_comp = 0; i_comp < d; i_comp++) {
126  for (size_type j_comp = 0; j_comp < d; j_comp++) {
127  h_value(i_comp,j_comp) = hh_comp[i_comp][j_comp].dof (idof);
128  }}
129  const bool use_svd_when_2d = true;
130  // H = Q*diag(lambda)*Q^T
131  if (use_svd_when_2d && d == 2) {
132  // TODO: eig when d=2 is bad...
133  lambda = h_value.svd (Q, Qt, d);
134  } else {
135  lambda = h_value.eig (Q, d);
136  }
137  T h_min = std::numeric_limits<T>::max();
138  for (size_type i_comp = 0; i_comp < d; i_comp++) {
139  // err = c*h^2
140  h_local[i_comp] = 1/sqrt(fabs(lambda[i_comp]));
141  // yield value:
142  h_local[i_comp] = std::max (opts.hmin, h_local[i_comp]);
143  h_min = std::min (h_min, h_local[i_comp]);
144  }
145  // isotropic: takes h_local = h_min in all directions
146  sh.dof (idof) = h_min;
147 
148  // anisotropic
149  m_value = Q*diag(h_local)*trans(Q);
150  for (size_type i_comp = 0; i_comp < d; i_comp++) {
151  for (size_type j_comp = 0; j_comp < d; j_comp++) {
152  mh_comp[i_comp][j_comp].dof (idof) = m_value(i_comp,j_comp);
153  }}
154  }
155  T cut_off = 1e-5;
156  T uh_scale = std::max(cut_off, fabs(uh.max() - uh.min()));
157  T factor = opts.hcoef*sqrt(uh_scale)*pow(opts.err,1./(k+1));
158  mh *= factor;
159  sh *= factor;
160  if (opts.isotropic) {
161  return smooth (sh, opts.n_smooth_metric);
162  } else {
163  return smooth (mh, opts.n_smooth_metric);
164  }
165 }
166 // -----------------------------------------
167 // adapt
168 // -----------------------------------------
170 template<class T, class M>
171 geo_basic<T,M>
173  const field_basic<T,M>& uh,
174  const adapt_option& opts)
175 {
176  size_t d = uh.get_geo().dimension();
177  if (d == 2 && (opts.generator == "bamg" || opts.generator == "")) {
178  // default generator is bamg when d=2:
179  return adapt_bamg (uh, opts);
180  } else {
181  // use always gmsh when d != 2:
182  return adapt_gmsh (uh, opts);
183  }
184 }
185 // ----------------------------------------------------------------------------
186 // instanciation in library
187 // ----------------------------------------------------------------------------
188 #define _RHEOLEF_instanciation(T,M) \
189 template field_basic<T,M> proj (const field_basic<T,M>&, const std::string&); \
190 template field_basic<T,M> smooth (const field_basic<T,M>&, size_t); \
191 template field_basic<T,M> hessian (const field_basic<T,M>&); \
192 template field_basic<T,M> hessian_criterion (const field_basic<T,M>&, const adapt_option&); \
193 template geo_basic<T,M> adapt (const field_basic<T,M>&, const adapt_option&);
194 
195 _RHEOLEF_instanciation(Float,sequential)
196 #ifdef _RHEOLEF_HAVE_MPI
198 #endif // _RHEOLEF_HAVE_MPI
199 
200 } // namespace rheolef
rheolef::space_numbering::ndof
size_type ndof(const basis_basic< T > &b, const geo_size &gs, size_type map_dim)
Definition: space_numbering.cc:28
rheolef::geo_basic< T, M >
rheolef::field_basic::get_space
const space_type & get_space() const
Definition: field.h:300
rheolef::adapt_bamg
geo_basic< T, M > adapt_bamg(const field_basic< T, M > &uh, const adapt_option &opts)
Definition: adapt_bamg.cc:63
rheolef::adapt
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:172
gh
field gh(Float epsilon, Float t, const field &uh, const test &v)
Definition: burgers_diffusion_operators.icc:37
rheolef::point_basic
Definition: point.h:87
check_macro
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
rheolef::_RHEOLEF_instanciation
_RHEOLEF_instanciation(Float, sequential, std::allocator< Float >) _RHEOLEF_instanciation(Float
rheolef::tensor_basic::eig
point_basic< T > eig(tensor_basic< T > &q, size_t dim=3) const
Definition: tensor.cc:426
rheolef::diag
csr< T, M > diag(const vec< T, M > &d)
Definition: csr.cc:56
rheolef::adapt_option
adapt_option: see the adapt page for the full documentation
Definition: adapt.h:147
rheolef::h_local
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
Definition: field_expr_terminal.h:527
rheolef::field_component_const
Definition: field.h:207
rheolef::solver_basic
Definition: solver.h:188
rheolef::field_basic::ndof
size_type ndof() const
Definition: field.h:341
rheolef::space_basic
the finite element space
Definition: space.h:352
rheolef::tensor_basic
Definition: tensor.h:90
rheolef::field_component::proxy_assign
field_component< T, M > & proxy_assign(field_component< T, M > &&uh_comp)
Definition: field_component.h:290
rheolef::form_basic::uu
const csr< T, M > & uu() const
Definition: form.h:192
rheolef::field_component::dof
T & dof(size_type idof)
Definition: field_component.h:115
rheolef::pow
space_mult_list< T, M > pow(const space_basic< T, M > &X, size_t n)
Definition: space_mult.h:120
rheolef::size_type
size_t size_type
Definition: basis_get.cc:76
rheolef::adapt_option::hcoef
Float hcoef
Definition: adapt.h:153
rheolef::adapt_option::generator
std::string generator
Definition: adapt.h:149
p
Definition: sphere.icc:25
rheolef::field_basic::get_geo
const geo_type & get_geo() const
Definition: field.h:301
mkgeo_contraction.md
md
Definition: mkgeo_contraction.sh:192
rheolef::field_basic::set_u
vec< T, M > & set_u()
Definition: field.h:312
rheolef::adapt_option::hmin
Float hmin
Definition: adapt.h:154
rheolef::field_basic::min
T min() const
Definition: field.h:683
rheolef::adapt_option::isotropic
bool isotropic
Definition: adapt.h:150
rheolef::field_basic::max
T max() const
Definition: field.h:699
rheolef::field_basic
Definition: field_expr_utilities.h:38
rheolef::field_basic::dof
T & dof(size_type idof)
Definition: field.h:658
rheolef::tensor_basic::svd
point_basic< T > svd(tensor_basic< T > &u, tensor_basic< T > &v, size_t dim=3) const
Definition: tensor.cc:470
rheolef::field_component
Definition: field.h:206
mkgeo_sector.m
m
Definition: mkgeo_sector.sh:118
rheolef::adapt_option::n_smooth_metric
size_type n_smooth_metric
Definition: adapt.h:159
rheolef
This file is part of Rheolef.
Definition: compiler_eigen.h:37
rheolef::field_component_const::proxy_assign
field_component_const< T, M > & proxy_assign(const field_component_const< T, M > &uh_comp)
Definition: field_component.h:280
rheolef::adapt_option::err
Float err
Definition: adapt.h:151
rheolef::field_component_const::dof
const T & dof(size_type idof) const
Definition: field_component.h:169
Float
see the Float page for the full documentation
mkgeo_ball.d
d
Definition: mkgeo_ball.sh:154
rheolef::proj
field_basic< T, M > proj(const field_basic< T, M > &uh, const std::string &approx="P1")
Definition: adapt.cc:39
u
Float u(const point &x)
Definition: transmission_error.cc:26
rheolef::form_basic
Definition: form.h:144
rheolef::field_basic::valued
const std::string & valued() const
Definition: field.h:305
mkgeo_ball.n
n
Definition: mkgeo_ball.sh:150
size_type
field::size_type size_type
Definition: branch.cc:425
rheolef::hessian
field_basic< T, M > hessian(const field_basic< T, M > &uh)
Definition: adapt.cc:73
rheolef::smooth
field_basic< T, M > smooth(const field_basic< T, M > &uh, size_t n=1)
Definition: adapt.cc:56
rheolef::trans
csr< T, sequential > trans(const csr< T, sequential > &a)
trans(a): see the form page for the full documentation
Definition: csr.h:455
rheolef::distributed
distributed
Definition: asr.cc:228
rheolef::hessian_criterion
field_basic< T, M > hessian_criterion(const field_basic< T, M > &uh0, const adapt_option &opts)
Definition: adapt.cc:104
rheolef::adapt_gmsh
geo_basic< T, M > adapt_gmsh(const field_basic< T, M > &uh, const adapt_option &opts)
Definition: adapt_gmsh.cc:41
rheolef::solver_basic::solve
vec< T, M > solve(const vec< T, M > &b) const
Definition: solver.h:345
T
Expr1::float_type T
Definition: field_expr.h:218
lambda
Definition: yield_slip_circle.h:34