rheolef  6.5
adapt.cc
Go to the documentation of this file.
1 #include "rheolef/adapt.h"
2 #include "rheolef/form.h"
3 #include "rheolef/field_component.h"
4 #include "rheolef/field_expr_ops.h"
5 
6 namespace rheolef {
7 
8 template<class T, class M>
9 geo_basic<T,M> adapt_gmsh (const field_basic<T,M>& uh, const adapt_option_type& opts);
10 template<class T, class M>
11 geo_basic<T,M> adapt_bamg (const field_basic<T,M>& uh, const adapt_option_type& opts);
12 
13 template<class T, class M>
14 field_basic<T,M>
15 proj (const field_basic<T,M>& uh, const std::string& approx = "P1")
16 {
17  const space_basic<T,M>& Uh = uh.get_space();
18  if (Uh.get_approx() == approx) return uh;
19  space_basic<T,M> Vh (uh.get_geo(), approx, uh.valued());
20  form_basic<T,M> m (Vh, Vh, "lumped_mass");
21  form_basic<T,M> p (Uh, Vh, "mass");
22  solver_basic<T,M> sm (m.uu());
23  field_basic<T,M> vh (Vh, 0);
24  vh.set_u() = sm.solve((p*uh).u());
25  return vh;
26 }
27 template<class T, class M>
28 field_basic<T,M>
29 smooth (const field_basic<T,M>& uh, size_t n = 1) {
30  const space_basic<T,M>& Xh = uh.get_space();
31  check_macro (Xh.get_approx() == "P1", "smooth: expect P1 approx, but have " << Xh.get_approx());
32  form_basic<T,M> m (Xh, Xh, "mass");
33  form_basic<T,M> md (Xh, Xh, "lumped_mass");
34  solver_basic<T,M> smd (md.uu());
35  field_basic<T,M> vh = uh;
36  for (size_t i = 0; i < n; i++) {
37  vh.set_u() = smd.solve ((m*vh).u());
38  }
39  return vh;
40 }
41 template<class T, class M>
42 field_basic<T,M>
44 {
45  const geo_basic<T,M>& omega = uh.get_geo();
46  const space_basic<T,M>& Xh = uh.get_space();
47  check_macro (Xh.valued() == "scalar", "hessian: unexpected "<<Xh.valued()<<"-valued field");
48  check_macro (Xh.get_approx() == "P1", "hessian: unexpected "<<Xh.get_approx()<<" approximation");
49  space_basic<T,M> Vh (omega, "P1", "vector");
50  form_basic<T,M> bv (Xh, Vh, "grad");
51  form_basic<T,M> mv (Vh, Vh, "lumped_mass");
52  solver_basic<T,M> smv (mv.uu());
53  field_basic<T,M> gh (Vh, 0);
54  gh.set_u() = smv.solve ((bv*uh).u());
55  space_basic<T,M> Th (omega, "P1", "tensor");
56  form_basic<T,M> bt (Vh, Th, "2D");
57  form_basic<T,M> mt (Th, Th, "lumped_mass");
58  solver_basic<T,M> smt (mt.uu());
59  field_basic<T,M> hh (Th, 0);
60  hh.set_u() = smt.solve ((bt*gh).u());
61  return hh;
62 }
63 template<class T, class M>
64 field_basic<T,M>
66  const field_basic<T,M>& uh0,
67  const adapt_option_type& opts)
68 {
69  typedef typename geo_basic<T,M>::size_type size_type;
70  size_type d = uh0.get_geo().dimension();
71  size_type k = uh0.get_space().degree();
72  field_basic<T,M> uh = proj(uh0);
73  field_basic<T,M> hh = hessian(uh);
74  field_basic<T,M> mh (hh.get_space(), 0);
75  field_basic<T,M> sh (uh.get_space(), 0);
76  field_component_const<T,M> hh_comp [3][3];
77  field_component<T,M> mh_comp [3][3];
78  for (size_type i_comp = 0; i_comp < d; i_comp++) {
79  for (size_type j_comp = 0; j_comp < d; j_comp++) {
80  hh_comp[i_comp][j_comp] = hh(i_comp,j_comp);
81  mh_comp[i_comp][j_comp] = mh(i_comp,j_comp);
82  }}
83  tensor_basic<T> h_value, m_value, Q, Qt;
84  point_basic<T> lambda, h_local;
85  for (size_type idof = 0, ndof = uh.ndof(); idof < ndof; idof++) {
86  for (size_type i_comp = 0; i_comp < d; i_comp++) {
87  for (size_type j_comp = 0; j_comp < d; j_comp++) {
88  h_value(i_comp,j_comp) = hh_comp[i_comp][j_comp].dof (idof);
89  }}
90  const bool use_svd_when_2d = true;
91  if (use_svd_when_2d && d == 2) {
92  lambda = h_value.svd (Q, Qt, d);
93  } else {
94  lambda = h_value.eig (Q, d);
95  }
96  T h_min = std::numeric_limits<T>::max();
97  for (size_type i_comp = 0; i_comp < d; i_comp++) {
98  h_local[i_comp] = 1/sqrt(fabs(lambda[i_comp]));
99  h_local[i_comp] = std::max (opts.hmin, h_local[i_comp]);
100  h_min = std::min (h_min, h_local[i_comp]);
101  }
102  sh.dof (idof) = h_min;
103 
104  m_value = Q*diag(h_local)*trans(Q);
105  for (size_type i_comp = 0; i_comp < d; i_comp++) {
106  for (size_type j_comp = 0; j_comp < d; j_comp++) {
107  mh_comp[i_comp][j_comp].dof (idof) = m_value(i_comp,j_comp);
108  }}
109  }
110  T cut_off = 1e-5;
111  T uh_scale = std::max(cut_off, fabs(uh.max() - uh.min()));
112  T factor = opts.hcoef*sqrt(uh_scale)*pow(opts.err,1./(k+1));
113  mh *= factor;
114  sh *= factor;
115  if (opts.isotropic) {
116  return smooth (sh, opts.n_smooth_metric);
117  } else {
118  return smooth (mh, opts.n_smooth_metric);
119  }
120 }
121 template<class T, class M>
122 geo_basic<T,M>
124  const field_basic<T,M>& uh,
125  const adapt_option_type& opts)
126 {
127  size_t d = uh.get_geo().dimension();
128  if (d == 2 && (opts.generator == "bamg" || opts.generator == "")) {
129  return adapt_bamg (uh, opts);
130  } else {
131  // use always gmsh when d != 2:
132  return adapt_gmsh (uh, opts);
133  }
134 }
135 #define _RHEOLEF_instanciation(T,M) \
136 template field_basic<T,M> proj (const field_basic<T,M>&, const std::string&); \
137 template field_basic<T,M> smooth (const field_basic<T,M>&, size_t); \
138 template field_basic<T,M> hessian (const field_basic<T,M>&); \
139 template field_basic<T,M> hessian_criterion (const field_basic<T,M>&, const adapt_option_type&); \
140 template geo_basic<T,M> adapt (const field_basic<T,M>&, const adapt_option_type&);
141 
142 _RHEOLEF_instanciation(Float,sequential)
143 #ifdef _RHEOLEF_HAVE_MPI
145 #endif // _RHEOLEF_HAVE_MPI
146 
147 } // namespace rheolef
148