rheolef  7.0
Go to the documentation of this file.
2 #include "rheolef/form.h"
3 #include "rheolef/field_component.h"
4 #include "rheolef/field_expr_v2_linear.h"
5
6 namespace rheolef {
7
8 template<class T, class M>
10 template<class T, class M>
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>
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>
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>
66  const field_basic<T,M>& uh0,
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].proxy_assign (hh(i_comp,j_comp));
81  mh_comp[i_comp][j_comp].proxy_assign (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>
124  const field_basic<T,M>& uh,
126 {
127  size_t d = uh.get_geo().dimension();
128  if (d == 2 && (opts.generator == "bamg" || opts.generator == "")) {
130  } else {
131  // use always gmsh when d != 2:
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&); \
141
142 _RHEOLEF_instanciation(Float,sequential)
143 #ifdef _RHEOLEF_HAVE_MPI
144 _RHEOLEF_instanciation(Float,distributed)
145 #endif // _RHEOLEF_HAVE_MPI
146
147 } // namespace rheolef
size_type ndof() const
Definition: field.h:318
field - piecewise polynomial finite element field
point_basic< T > svd(tensor_basic< T > &u, tensor_basic< T > &v, size_t dim=3) const
Definition: tensor.cc:348
T & dof(size_type idof)
field_basic< T, M > hessian_criterion(const field_basic< T, M > &uh0, const adapt_option &opts)
space_mult_list< T, M > pow(const space_basic< T, M > &X, size_t n)
Definition: space_mult.h:87
point - vertex of a mesh
Definition: point.h:22
vec< T, M > solve(const vec< T, M > &b) const
Definition: solver.h:275
T max() const
Definition: field.h:592
irheostream, orheostream - large data streams
Definition: compiler.h:7
details::field_expr_v2_nonlinear_terminal_function< details::h_local_pseudo_function< Float > > h_local()
T & dof(size_type idof)
Definition: field.h:553
const std::string & valued() const
Definition: field.h:288
const geo_type & get_geo() const
Definition: field.h:284
T min() const
Definition: field.h:576
csr< T, M > diag(const vec< T, M > &d)
Definition: csr.cc:33
space – piecewise polynomial finite element space
Definition: space.h:229
#define _RHEOLEF_instanciation(T, M)
#define check_macro(ok_condition, message)
Definition: compiler.h:104
double Float
Definition: compiler.h:160
std::string generator
size_t d
solver - direct or interative solver interface
Definition: solver.h:8
point_basic< T > eig(tensor_basic< T > &q, size_t dim=3) const
Definition: tensor.cc:328
field_component< T, M > & proxy_assign(field_component< T, M > &&uh_comp)
geo_basic< T, M > adapt_bamg(const field_basic< T, M > &uh, const adapt_option &opts)
const csr< T, M > & uu() const
Definition: form.h:140
field_basic< T, M > smooth(const field_basic< T, M > &uh, size_t n=1)
const space_type & get_space() const
Definition: field.h:283
geo_basic< T, M > adapt_gmsh(const field_basic< T, M > &uh, const adapt_option &opts)
field_basic< T, M > proj(const field_basic< T, M > &uh, const std::string &approx="P1")
field_basic< T, M > hessian(const field_basic< T, M > &uh)
size_type n_smooth_metric
form - representation of a finite element bilinear form
Definition: form.h:98
csr< T, sequential > trans(const csr< T, sequential > &a)
Definition: csr.h:381
Float u(const point &x)
geo_basic< T, M > adapt(const field_basic< T, M > &uh, const adapt_option &opts)