Rheolef  7.1
an efficient C++ finite element environment
form_weighted.h
Go to the documentation of this file.
1 # ifndef _RHEOLEF_FORM_WEIGHTED_H
2 # define _RHEOLEF_FORM_WEIGHTED_H
3 // form constructor: old interface with name
24 // maintained for backward compatibility purpose
25 //
26 #include "rheolef/form.h"
27 #include "rheolef/field_expr_variational.h"
28 #include "rheolef/form_expr_variational.h"
29 #include "rheolef/form_vf_assembly.h"
30 namespace rheolef {
31 
32 // backward compat with named forms:
33 // for each weight, it compiles all cases
34 // => not very efficient !
35 
36 namespace details {
37 template<class Expr>
38 static
39 inline
40 form_expr_quadrature_on_element<Expr>
41 expr (const Expr& e) {
42  return form_expr_quadrature_on_element<Expr>(e);
43 }
44 
45 template<class T, class M, class WeightFunction>
46 bool
49  const geo_basic<T,M>& dom,
50  const std::string& name,
51  bool has_weight,
52  WeightFunction w,
53  const quadrature_option& qopt)
54 {
55  // backward compatibility code:
56  // branch to variational formulation routines:
57  test_basic<T,M,details::vf_tag_10> u (a.get_first_space()); // trial
58  test_basic<T,M,details::vf_tag_01> v (a.get_second_space()); // test
59  // --------------------------------
60  if (name == "mass") {
61  // --------------------------------
62  switch (a.get_first_space().valued_tag()) {
64  if (!has_weight) a.assembly (dom, expr(u*v), qopt);
65  else a.assembly (dom, expr(w*(u*v)), qopt);
66  return true;
68  if (!has_weight) a.assembly (dom, expr(dot(u,v)), qopt);
69  else fatal_macro ("unsupported vectorial mass with weight (HINT: use integrate())");
70  // a.assembly (dom, expr(dot(w*u,v)), qopt); compil pbs ?
71  return true;
72  default:
74  if (!has_weight) a.assembly (dom, expr(ddot(u,v)), qopt);
75  else fatal_macro ("unsupported tensorial mass with weight (HINT: use integrate())");
76  // a.assembly (dom, expr(ddot(w*u,v)), qopt); compil pbs ?
77  return true;
78  }
79  // --------------------------------
80  } else if (name == "inv_mass") {
81  // --------------------------------
82  check_macro (!has_weight, "unsupported weighted "<<name<<" form (HINT: use integrate())");
83  integrate_option fopt (qopt);
84  fopt.invert = true;
85  switch (a.get_first_space().valued_tag()) {
86  case space_constant::scalar: a.assembly (dom, expr(u*v), fopt); return true;
87  case space_constant::vector: a.assembly (dom, expr(dot(u,v)),fopt); return true;
88  default:
89  case space_constant::tensor: a.assembly (dom, expr(ddot(u,v)), fopt); return true;
90  }
91  // --------------------------------
92  } else if (name == "lumped_mass") {
93  // --------------------------------
94  check_macro (!has_weight, "unsupported weighted "<<name<<" form (HINT: use integrate())");
95  integrate_option fopt (qopt);
96  fopt.lump = true;
97  switch (a.get_first_space().valued_tag()) {
98  case space_constant::scalar: a.assembly (dom, expr(u*v), fopt); return true;
99  case space_constant::vector: a.assembly (dom, expr(dot(u,v)), fopt); return true;
100  default:
101  case space_constant::tensor: a.assembly (dom, expr(ddot(u,v)), fopt); return true;
102  }
103  // --------------------------------
104  } else if (name == "grad") {
105  // --------------------------------
106  if (!has_weight) a.assembly (dom, expr(dot(grad(u),v)), qopt);
107  else a.assembly (dom, expr(dot(w*grad(u),v)), qopt);
108  return true;
109  // --------------------------------
110  } else if (name == "div") {
111  // --------------------------------
112  check_macro (!has_weight, "unsupported weighted "<<name<<" form (HINT: use integrate())");
113  a.assembly (dom, expr(div(u)*v), qopt);
114  return true;
115  // --------------------------------
116  } else if (name == "2D") {
117  // --------------------------------
118  if (!has_weight) a.assembly (dom, expr(2.*ddot(D(u),v)), qopt);
119  else a.assembly (dom, expr(2.*ddot(w*D(u),v)), qopt);
120  return true;
121  // --------------------------------
122  } else if (name == "grad_grad") {
123  // --------------------------------
124  switch (a.get_first_space().valued_tag()) {
126  if (!has_weight) a.assembly (dom, expr(dot(grad(u),grad(v))), qopt);
127  else a.assembly (dom, expr(dot(w*grad(u),grad(v))), qopt);
128  return true;
129  default:
131  if (!has_weight) a.assembly (dom, expr(ddot(grad(u),grad(v))), qopt);
132  else a.assembly (dom, expr(ddot(w*grad(u),grad(v))), qopt);
133  return true;
134  }
135  // --------------------------------
136  } else if (name == "2D_D") {
137  // --------------------------------
138  if (!has_weight) a.assembly (dom, expr(2.*ddot(D(u),D(v))), qopt);
139  else a.assembly (dom, expr(2.*ddot(w*D(u),D(v))), qopt);
140  return true;
141  // --------------------------------
142  } else if (name == "div_div") {
143  // --------------------------------
144  check_macro (!has_weight, "unsupported weighted "<<name<<" form (HINT: use integrate())");
145  a.assembly (dom, expr(div(u)*div(v)), qopt);
146  return true;
147  // --------------------------------
148  } else if (name == "curl") {
149  // --------------------------------
150  check_macro (!has_weight, "unsupported weighted "<<name<<" form (HINT: use integrate())");
151  if (u.get_vf_space().get_geo().dimension() == 2 &&
152  u.get_vf_space().valued_tag() == space_constant::vector)
153  a.assembly (dom, expr(curl(u)*v), qopt);
154  else a.assembly (dom, expr(dot(curl(u),v)), qopt);
155  return true;
156  }
157  return false;
158 }
159 
160 } // namespace details
161 
162 template<class T, class M>
163 template<class WeightFunction>
164 void
166  const std::string& name,
167  bool has_weight,
168  WeightFunction w,
169  const quadrature_option& qopt)
170 {
171  if (name == "" || name == "nul" || name == "null") {
172  // empty name => nul form, but with declared csr matrix sizes
173  _uu.resize (_Y.iu_ownership(), _X.iu_ownership());
174  _ub.resize (_Y.iu_ownership(), _X.ib_ownership());
175  _bu.resize (_Y.ib_ownership(), _X.iu_ownership());
176  _bb.resize (_Y.ib_ownership(), _X.ib_ownership());
177  return;
178  }
179  // determine the domain of integration:
180  typedef typename form_basic<T,M>::float_type float_type;
181  check_macro (_X.get_geo().get_background_geo() == _Y.get_geo().get_background_geo(),
182  "form("<<name<<") between incompatible geo " << _X.get_geo().name() << " and " << _Y.get_geo().name());
183  bool X_is_on_domain = (_X.get_geo().variant() == geo_abstract_base_rep<float_type>::geo_domain);
184  bool Y_is_on_domain = (_Y.get_geo().variant() == geo_abstract_base_rep<float_type>::geo_domain);
185  geo_basic<T,M> dom;
186  if ((!X_is_on_domain && ! Y_is_on_domain) || (X_is_on_domain && Y_is_on_domain)) {
187  dom = _X.get_geo();
188  } else if (X_is_on_domain) {
189  dom = _X.get_geo().get_background_domain();
190  } else {// Y_is_on_domain
191  dom = _Y.get_geo().get_background_domain();
192  }
193  // try the new form initializer interface:
194  if (details::form_named_init (*this, dom, name, has_weight, w, qopt)) return;
195 
196  fatal_macro ("unsupported form name: \""<<name<<"\" (HINT: use integrate())");
197 }
198 template<class T, class M>
199 template<class Function>
200 inline
202  const space_type& X,
203  const space_type& Y,
204  const std::string& name,
205  Function weight,
206  const quadrature_option& qopt)
207  : _X(X),
208  _Y(Y),
209  _uu(),
210  _ub(),
211  _bu(),
212  _bb()
213 {
214  form_init (name, true, weight, qopt);
215 }
216 // ================================================================
217 // also for Robin boundary conditions with a non-constant coef
218 // ================================================================
219 template<class T, class M>
220 template<class WeightFunction>
221 void
223  const std::string& name,
224  const geo_basic<T,M>& gamma,
225  bool has_weight,
226  WeightFunction w,
227  const geo_basic<T,M>& w_omega,
228  const quadrature_option& qopt)
229 {
230  const geo_basic<T,M>& omega = _X.get_geo().get_background_geo();
231 
232  // try the new interface:
233  if (details::form_named_init (*this, gamma, name, has_weight, w, qopt)) return;
234  fatal_macro ("unsupported form name: \""<<name<<"\"");
235 }
236 template<class T, class M>
237 template<class Function>
239  const space_type& X,
240  const space_type& Y,
241  const std::string& name,
242  const geo_basic<T,M>& gamma,
243  Function weight,
244  const quadrature_option& qopt)
245  : _X(X),
246  _Y(Y),
247  _uu(),
248  _ub(),
249  _bu(),
250  _bb()
251 {
252  // example:
253  // form m (V,V,"mass",gamma, weight); e.g. \int_\Gamma trace(u) trace(v) weight(x) ds
254  // with:
255  // geo omega ("square");
256  // geo gamma = omega["boundary"];
257  // V = space(omega,"P1");
258  // Note: transform the domain gamma into a compacted-mesh for the
259  // evaluation of the weight
260  form_init_on_domain (name, gamma, true, weight, compact(gamma), qopt);
261 }
262 
263 }// namespace rheolef
264 # endif /* _RHEOLEF_FORM_WEIGHTED_H */
rheolef::compact
geo_basic< T, M > compact(const geo_basic< T, M > &gamma)
Definition: geo.cc:219
rheolef::geo_basic
generic mesh with rerefence counting
Definition: domain_indirect.h:64
rheolef::div
std::enable_if< details::is_field_convertible< Expr >::value,details::field_expr_v2_nonlinear_terminal_field< typename Expr::scalar_type,typename Expr::memory_type,details::differentiate_option::divergence >>::type div(const Expr &expr)
div(uh): see the expression page for the full documentation
Definition: field_expr_terminal.h:1031
rheolef::form_basic::form_init_on_domain
void form_init_on_domain(const std::string &name, const geo_basic< T, M > &gamma, bool has_weight, WeightFunction weight, const geo_basic< T, M > &w_omega, const quadrature_option &qopt)
Definition: form_weighted.h:222
mkgeo_ball.expr
expr
Definition: mkgeo_ball.sh:361
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::space_basic< float_type, M >
rheolef::form_basic::float_type
scalar_traits< T >::type float_type
Definition: form.h:150
rheolef::float_type
typename float_traits< value_type >::type float_type
Definition: field_expr_recursive.h:501
rheolef::grad
std::enable_if< details::is_field_convertible< Expr >::value,details::field_expr_v2_nonlinear_terminal_field< typename Expr::scalar_type,typename Expr::memory_type,details::differentiate_option::gradient >>::type grad(const Expr &expr)
grad(uh): see the expression page for the full documentation
Definition: field_expr_terminal.h:911
rk::gamma
Float gamma[][pmax+1]
Definition: runge_kutta_semiimplicit.icc:70
rheolef::integrate_option::lump
bool lump
Definition: integrate_option.h:168
rheolef::geo_abstract_base_rep
abstract base interface class
Definition: geo.h:248
rheolef::form_basic::form_basic
form_basic()
Definition: form.h:300
rheolef::integrate_option
see the integrate_option page for the full documentation
Definition: integrate_option.h:125
rheolef::space_constant::tensor
Definition: space_constant.h:138
rheolef::space_constant::scalar
Definition: space_constant.h:136
rheolef::test_basic
Definition: test.h:207
rheolef::ddot
T ddot(const tensor_basic< T > &a, const tensor_basic< T > &b)
ddot(x,y): see the expression page for the full documentation
Definition: tensor.cc:278
a
Definition: diffusion_isotropic.h:25
rheolef::details::dot
rheolef::details::is_vec dot
rheolef::form_basic::form_init
void form_init(const std::string &name, bool has_weight, WeightFunction weight, const quadrature_option &qopt)
Definition: form_weighted.h:165
fatal_macro
#define fatal_macro(message)
Definition: dis_macros.h:33
rheolef
This file is part of Rheolef.
Definition: compiler_eigen.h:37
rheolef::integrate_option::invert
bool invert
Definition: integrate_option.h:168
u
Definition: leveque.h:25
rheolef::Function
rheolef::std Function
rheolef::curl
std::enable_if< details::is_field_convertible< Expr >::value,details::field_expr_v2_nonlinear_terminal_field< typename Expr::scalar_type,typename Expr::memory_type,details::differentiate_option::curl >>::type curl(const Expr &expr)
curl(uh): see the expression page for the full documentation
Definition: field_expr_terminal.h:1089
rheolef::D
std::enable_if< details::is_field_convertible< Expr >::value,details::field_expr_v2_nonlinear_terminal_field< typename Expr::scalar_type,typename Expr::memory_type,details::differentiate_option::gradient >>::type D(const Expr &expr)
D(uh): see the expression page for the full documentation.
Definition: field_expr_terminal.h:969
u
Float u(const point &x)
Definition: transmission_error.cc:26
rheolef::form_basic
Definition: form.h:144
rheolef::space_constant::vector
Definition: space_constant.h:137
rheolef::details::form_named_init
bool form_named_init(form_basic< T, M > &a, const geo_basic< T, M > &dom, const std::string &name, bool has_weight, WeightFunction w, const quadrature_option &qopt)
Definition: form_weighted.h:47
rheolef::quadrature_option
integrate_option quadrature_option
Definition: integrate_option.h:186
mkgeo_contraction.name
name
Definition: mkgeo_contraction.sh:133