Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
form_weighted.h
Go to the documentation of this file.
1# ifndef _RHEOLEF_FORM_WEIGHTED_H
2# define _RHEOLEF_FORM_WEIGHTED_H
23// 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"
30namespace rheolef {
31
32// backward compat with named forms:
33// for each weight, it compiles all cases
34// => not very efficient !
35
36namespace details {
37template<class Expr>
38static
39inline
40form_expr_quadrature_on_element<Expr>
41expr (const Expr& e) {
42 return form_expr_quadrature_on_element<Expr>(e);
43}
44
45template<class T, class M, class WeightFunction>
46bool
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.do_integrate (dom, expr(u*v), qopt);
65 else a.do_integrate (dom, expr(w*(u*v)), qopt);
66 return true;
68 if (!has_weight) a.do_integrate (dom, expr(dot(u,v)), qopt);
69 else fatal_macro ("unsupported vectorial mass with weight (HINT: use integrate())");
70 // a.do_integrate (dom, expr(dot(w*u,v)), qopt); compil pbs ?
71 return true;
72 default:
74 if (!has_weight) a.do_integrate (dom, expr(ddot(u,v)), qopt);
75 else fatal_macro ("unsupported tensorial mass with weight (HINT: use integrate())");
76 // a.do_integrate (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.do_integrate (dom, expr(u*v), fopt); return true;
87 case space_constant::vector: a.do_integrate (dom, expr(dot(u,v)),fopt); return true;
88 default:
89 case space_constant::tensor: a.do_integrate (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.do_integrate (dom, expr(u*v), fopt); return true;
99 case space_constant::vector: a.do_integrate (dom, expr(dot(u,v)), fopt); return true;
100 default:
101 case space_constant::tensor: a.do_integrate (dom, expr(ddot(u,v)), fopt); return true;
102 }
103 // --------------------------------
104 } else if (name == "grad") {
105 // --------------------------------
106 if (!has_weight) a.do_integrate (dom, expr(dot(grad(u),v)), qopt);
107 else a.do_integrate (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.do_integrate (dom, expr(div(u)*v), qopt);
114 return true;
115 // --------------------------------
116 } else if (name == "2D") {
117 // --------------------------------
118 if (!has_weight) a.do_integrate (dom, expr(2.*ddot(D(u),v)), qopt);
119 else a.do_integrate (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.do_integrate (dom, expr(dot(grad(u),grad(v))), qopt);
127 else a.do_integrate (dom, expr(dot(w*grad(u),grad(v))), qopt);
128 return true;
129 default:
131 if (!has_weight) a.do_integrate (dom, expr(ddot(grad(u),grad(v))), qopt);
132 else a.do_integrate (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.do_integrate (dom, expr(2.*ddot(D(u),D(v))), qopt);
139 else a.do_integrate (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.do_integrate (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.do_integrate (dom, expr(curl(u)*v), qopt);
154 else a.do_integrate (dom, expr(dot(curl(u),v)), qopt);
155 return true;
156 }
157 return false;
158}
159
160} // namespace details
161
162template<class T, class M>
163template<class WeightFunction>
164void
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:
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}
198template<class T, class M>
199template<class Function>
200inline
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// ================================================================
219template<class T, class M>
220template<class WeightFunction>
221void
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}
236template<class T, class M>
237template<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);
262
263}// namespace rheolef
264# endif /* _RHEOLEF_FORM_WEIGHTED_H */
scalar_traits< T >::type float_type
Definition form.h:152
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)
void form_init(const std::string &name, bool has_weight, WeightFunction weight, const quadrature_option &qopt)
abstract base interface class
Definition geo.h:248
generic mesh with rerefence counting
Definition geo.h:1089
see the integrate_option page for the full documentation
#define fatal_macro(message)
Definition dis_macros.h:33
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
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)
rheolef::details::is_vec dot
This file is part of Rheolef.
std::enable_if< details::has_field_rdof_interface< Expr >::value, details::field_expr_v2_nonlinear_terminal_field< typenameExpr::scalar_type, typenameExpr::memory_type, details::differentiate_option::gradient > >::type grad(const Expr &expr)
grad(uh): see the expression page for the full documentation
std::enable_if< details::has_field_rdof_interface< Expr >::value, details::field_expr_v2_nonlinear_terminal_field< typenameExpr::scalar_type, typenameExpr::memory_type, details::differentiate_option::gradient > >::type D(const Expr &expr)
D(uh): see the expression page for the full documentation.
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
std::enable_if< details::has_field_rdof_interface< Expr >::value, details::field_expr_v2_nonlinear_terminal_field< typenameExpr::scalar_type, typenameExpr::memory_type, details::differentiate_option::divergence > >::type div(const Expr &expr)
div(uh): see the expression page for the full documentation
geo_basic< T, M > compact(const geo_basic< T, M > &gamma)
Definition geo.cc:219
std::enable_if< details::has_field_rdof_interface< Expr >::value, details::field_expr_v2_nonlinear_terminal_field< typenameExpr::scalar_type, typenameExpr::memory_type, details::differentiate_option::curl > >::type curl(const Expr &expr)
curl(uh): see the expression page for the full documentation
Definition leveque.h:25