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.do_integrate (dom, expr(u*v),   qopt);
        else             a.do_integrate (dom, expr(w*(u*v)), qopt);
        return true;
      case space_constant::vector:
        if (!has_weight) a.do_integrate (dom, expr(dot(u,v)),   qopt);
        else             fatal_macro ("unsupported vectorial mass with weight (HINT: use integrate())");
                         // a.do_integrate (dom, expr(dot(w*u,v)), qopt); compil pbs ?
        return true;
      default:
      case space_constant::tensor:
        if (!has_weight) a.do_integrate (dom, expr(ddot(u,v)),   qopt);
        else             fatal_macro ("unsupported tensorial mass with weight (HINT: use integrate())");
                         // a.do_integrate (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.do_integrate (dom, expr(u*v),     fopt); return true;
      case space_constant::vector: a.do_integrate (dom, expr(dot(u,v)),fopt); return true;
      default:
      case space_constant::tensor: a.do_integrate (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.do_integrate (dom, expr(u*v),   fopt); return true;
      case space_constant::vector: a.do_integrate (dom, expr(dot(u,v)),   fopt); return true;
      default:
      case space_constant::tensor: a.do_integrate (dom, expr(ddot(u,v)),   fopt); return true;
    }
  // --------------------------------
  } else if (name == "grad") {
  // --------------------------------
        if (!has_weight) a.do_integrate (dom, expr(dot(grad(u),v)),   qopt);
        else             a.do_integrate (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.do_integrate (dom, expr(div(u)*v),   qopt);
        return true;
  // --------------------------------
  } else if (name == "2D") {
  // --------------------------------
        if (!has_weight) a.do_integrate (dom, expr(2.*ddot(D(u),v)),   qopt);
        else             a.do_integrate (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.do_integrate (dom, expr(dot(grad(u),grad(v))),   qopt);
        else             a.do_integrate (dom, expr(dot(w*grad(u),grad(v))), qopt);
        return true;
      default:
      case space_constant::vector:
        if (!has_weight) a.do_integrate (dom, expr(ddot(grad(u),grad(v))),   qopt);
        else             a.do_integrate (dom, expr(ddot(w*grad(u),grad(v))), qopt);
        return true;
    }
  // --------------------------------
  } else if (name == "2D_D") {
  // --------------------------------
        if (!has_weight) a.do_integrate (dom, expr(2.*ddot(D(u),D(v))),   qopt);
        else             a.do_integrate (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.do_integrate (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.do_integrate (dom, expr(curl(u)*v), qopt);
        else    a.do_integrate (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 */