form_weighted.h
Go to the documentation of this file.
# ifndef _RHEOLEF_FORM_WEIGHTED_H
# define _RHEOLEF_FORM_WEIGHTED_H
// form constructor: old interface with name
// maintained for backward compatibility purpose
//
#include "rheolef/form.h"
#include "rheolef/field_expr_variational.h"
#include "rheolef/form_expr_variational.h"
#include "rheolef/form_vf_assembly.h"
namespace rheolef {
// backward compat with named forms:
// for each weight, it compiles all cases
// => not very efficient !
namespace details {
template<class Expr>
static
inline
form_expr_quadrature_on_element<Expr>
expr (const Expr& e) {
return form_expr_quadrature_on_element<Expr>(e);
}
template<class T, class M, class WeightFunction>
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)
{
// backward compatibility code:
// branch to variational formulation routines:
test_basic<T,M,details::vf_tag_10> u (a.get_first_space()); // trial
test_basic<T,M,details::vf_tag_01> v (a.get_second_space()); // test
// --------------------------------
if (name == "mass") {
// --------------------------------
switch (a.get_first_space().valued_tag()) {
case space_constant::scalar:
if (!has_weight) a.assembly (dom, expr(u*v), qopt);
else a.assembly (dom, expr(w*(u*v)), qopt);
return true;
case space_constant::vector:
if (!has_weight) a.assembly (dom, expr(dot(u,v)), qopt);
else fatal_macro ("unsupported vectorial mass with weight (HINT: use integrate())");
// a.assembly (dom, expr(dot(w*u,v)), qopt); compil pbs ?
return true;
default:
case space_constant::tensor:
if (!has_weight) a.assembly (dom, expr(ddot(u,v)), qopt);
else fatal_macro ("unsupported tensorial mass with weight (HINT: use integrate())");
// a.assembly (dom, expr(ddot(w*u,v)), qopt); compil pbs ?
return true;
}
// --------------------------------
} else if (name == "inv_mass") {
// --------------------------------
check_macro (!has_weight, "unsupported weighted "<<name<<" form (HINT: use integrate())");
integrate_option fopt (qopt);
fopt.invert = true;
switch (a.get_first_space().valued_tag()) {
case space_constant::scalar: a.assembly (dom, expr(u*v), fopt); return true;
case space_constant::vector: a.assembly (dom, expr(dot(u,v)),fopt); return true;
default:
case space_constant::tensor: a.assembly (dom, expr(ddot(u,v)), fopt); return true;
}
// --------------------------------
} else if (name == "lumped_mass") {
// --------------------------------
check_macro (!has_weight, "unsupported weighted "<<name<<" form (HINT: use integrate())");
integrate_option fopt (qopt);
fopt.lump = true;
switch (a.get_first_space().valued_tag()) {
case space_constant::scalar: a.assembly (dom, expr(u*v), fopt); return true;
case space_constant::vector: a.assembly (dom, expr(dot(u,v)), fopt); return true;
default:
case space_constant::tensor: a.assembly (dom, expr(ddot(u,v)), fopt); return true;
}
// --------------------------------
} else if (name == "grad") {
// --------------------------------
if (!has_weight) a.assembly (dom, expr(dot(grad(u),v)), qopt);
else a.assembly (dom, expr(dot(w*grad(u),v)), qopt);
return true;
// --------------------------------
} else if (name == "div") {
// --------------------------------
check_macro (!has_weight, "unsupported weighted "<<name<<" form (HINT: use integrate())");
a.assembly (dom, expr(div(u)*v), qopt);
return true;
// --------------------------------
} else if (name == "2D") {
// --------------------------------
if (!has_weight) a.assembly (dom, expr(2.*ddot(D(u),v)), qopt);
else a.assembly (dom, expr(2.*ddot(w*D(u),v)), qopt);
return true;
// --------------------------------
} else if (name == "grad_grad") {
// --------------------------------
switch (a.get_first_space().valued_tag()) {
case space_constant::scalar:
if (!has_weight) a.assembly (dom, expr(dot(grad(u),grad(v))), qopt);
else a.assembly (dom, expr(dot(w*grad(u),grad(v))), qopt);
return true;
default:
case space_constant::vector:
if (!has_weight) a.assembly (dom, expr(ddot(grad(u),grad(v))), qopt);
else a.assembly (dom, expr(ddot(w*grad(u),grad(v))), qopt);
return true;
}
// --------------------------------
} else if (name == "2D_D") {
// --------------------------------
if (!has_weight) a.assembly (dom, expr(2.*ddot(D(u),D(v))), qopt);
else a.assembly (dom, expr(2.*ddot(w*D(u),D(v))), qopt);
return true;
// --------------------------------
} else if (name == "div_div") {
// --------------------------------
check_macro (!has_weight, "unsupported weighted "<<name<<" form (HINT: use integrate())");
a.assembly (dom, expr(div(u)*div(v)), qopt);
return true;
// --------------------------------
} else if (name == "curl") {
// --------------------------------
check_macro (!has_weight, "unsupported weighted "<<name<<" form (HINT: use integrate())");
if (u.get_vf_space().get_geo().dimension() == 2 &&
u.get_vf_space().valued_tag() == space_constant::vector)
a.assembly (dom, expr(curl(u)*v), qopt);
else a.assembly (dom, expr(dot(curl(u),v)), qopt);
return true;
}
return false;
}
} // namespace details
template<class T, class M>
template<class WeightFunction>
void
form_basic<T,M>::form_init (
const std::string& name,
bool has_weight,
WeightFunction w,
const quadrature_option& qopt)
{
if (name == "" || name == "nul" || name == "null") {
// empty name => nul form, but with declared csr matrix sizes
_uu.resize (_Y.iu_ownership(), _X.iu_ownership());
_ub.resize (_Y.iu_ownership(), _X.ib_ownership());
_bu.resize (_Y.ib_ownership(), _X.iu_ownership());
_bb.resize (_Y.ib_ownership(), _X.ib_ownership());
return;
}
// determine the domain of integration:
typedef typename form_basic<T,M>::float_type float_type;
check_macro (_X.get_geo().get_background_geo() == _Y.get_geo().get_background_geo(),
"form("<<name<<") between incompatible geo " << _X.get_geo().name() << " and " << _Y.get_geo().name());
bool X_is_on_domain = (_X.get_geo().variant() == geo_abstract_base_rep<float_type>::geo_domain);
bool Y_is_on_domain = (_Y.get_geo().variant() == geo_abstract_base_rep<float_type>::geo_domain);
geo_basic<T,M> dom;
if ((!X_is_on_domain && ! Y_is_on_domain) || (X_is_on_domain && Y_is_on_domain)) {
dom = _X.get_geo();
} else if (X_is_on_domain) {
dom = _X.get_geo().get_background_domain();
} else {// Y_is_on_domain
dom = _Y.get_geo().get_background_domain();
}
// try the new form initializer interface:
if (details::form_named_init (*this, dom, name, has_weight, w, qopt)) return;
fatal_macro ("unsupported form name: \""<<name<<"\" (HINT: use integrate())");
}
template<class T, class M>
template<class Function>
inline
form_basic<T,M>::form_basic (
const space_type& X,
const space_type& Y,
const std::string& name,
Function weight,
const quadrature_option& qopt)
: _X(X),
_Y(Y),
_uu(),
_ub(),
_bu(),
_bb()
{
form_init (name, true, weight, qopt);
}
// ================================================================
// also for Robin boundary conditions with a non-constant coef
// ================================================================
template<class T, class M>
template<class WeightFunction>
void
form_basic<T,M>::form_init_on_domain (
const std::string& name,
const geo_basic<T,M>& gamma,
bool has_weight,
WeightFunction w,
const geo_basic<T,M>& w_omega,
const quadrature_option& qopt)
{
const geo_basic<T,M>& omega = _X.get_geo().get_background_geo();
// try the new interface:
if (details::form_named_init (*this, gamma, name, has_weight, w, qopt)) return;
fatal_macro ("unsupported form name: \""<<name<<"\"");
}
template<class T, class M>
template<class Function>
form_basic<T,M>::form_basic (
const space_type& X,
const space_type& Y,
const std::string& name,
const geo_basic<T,M>& gamma,
Function weight,
const quadrature_option& qopt)
: _X(X),
_Y(Y),
_uu(),
_ub(),
_bu(),
_bb()
{
// example:
// form m (V,V,"mass",gamma, weight); e.g. \int_\Gamma trace(u) trace(v) weight(x) ds
// with:
// geo omega ("square");
// geo gamma = omega["boundary"];
// V = space(omega,"P1");
// Note: transform the domain gamma into a compacted-mesh for the
// evaluation of the weight
form_init_on_domain (name, gamma, true, weight, compact(gamma), qopt);
}
}// namespace rheolef
# endif /* _RHEOLEF_FORM_WEIGHTED_H */
rheolef::form_basic
Definition: form.h:144
rheolef::form_basic::float_type
scalar_traits< T >::type float_type
Definition: form.h:150
rheolef::form_basic::form_basic
form_basic()
Definition: form.h:300
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
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
rheolef::geo_abstract_base_rep
abstract base interface class
Definition: geo.h:248
rheolef::geo_basic< T, M >
rheolef::integrate_option
see the integrate_option page for the full documentation
Definition: integrate_option.h:125
rheolef::integrate_option::lump
bool lump
Definition: integrate_option.h:168
rheolef::integrate_option::invert
bool invert
Definition: integrate_option.h:168
rheolef::space_basic< float_type, M >
rheolef::test_basic
Definition: test.h:207
rheolef::float_type
typename float_traits< value_type >::type float_type
Definition: field_expr_recursive.h:501
rheolef::quadrature_option
integrate_option quadrature_option
Definition: integrate_option.h:186
rheolef::Function
rheolef::std Function
fatal_macro
#define fatal_macro(message)
Definition: dis_macros.h:33
check_macro
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
mkgeo_ball.expr
expr
Definition: mkgeo_ball.sh:361
mkgeo_contraction.name
string name
Definition: mkgeo_contraction.sh:133
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::details::dot
rheolef::details::is_vec dot
rheolef::space_constant::tensor
@ tensor
Definition: space_constant.h:138
rheolef::space_constant::scalar
@ scalar
Definition: space_constant.h:136
rheolef::space_constant::vector
@ vector
Definition: space_constant.h:137
rheolef
This file is part of Rheolef.
Definition: compiler_eigen.h:37
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::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
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
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
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::compact
geo_basic< T, M > compact(const geo_basic< T, M > &gamma)
Definition: geo.cc:219
rk::gamma
Float gamma[][pmax+1]
Definition: runge_kutta_semiimplicit.icc:70
a
Definition: diffusion_isotropic.h:25
u
Definition: leveque.h:25
u
Float u(const point &x)
Definition: transmission_error.cc:26