form.cc

Go to the documentation of this file.

//
#include "rheolef/field_expr.h"
#include "rheolef/geo_domain.h"
#include "rheolef/form.h"
#include "rheolef/csr.h"
#include "rheolef/field.h"
#include "rheolef/band.h"
#include "rheolef/form_weighted.h"
 
namespace rheolef { 
using namespace std;
 
template<class T, class M>
form_basic<T,M>::form_basic (
    const space_type& X, 
    const space_type& Y,
    const std::string& name,
    const quadrature_option& qopt)
  : _X(X), 
    _Y(Y), 
    _uu(),
    _ub(),
    _bu(),
    _bb()
{
  form_init (name, false, field_basic<T,M>(), qopt);
}
template<class T, class M>
form_basic<T,M>::form_basic (
    const space_type& X, 
    const space_type& Y,
    const std::string& name,
    const field_basic<T,M>& wh,
    const quadrature_option& qopt)
  : _X(X), 
    _Y(Y), 
    _uu(),
    _ub(),
    _bu(),
    _bb()
{
  form_init (name, true, details::field_expr_v2_nonlinear_terminal_field<T,M,details::differentiate_option::none>(wh), qopt);
}
template<class T, class M>
form_basic<T,M>::form_basic (
    const space_type& X, 
    const space_type& Y,
    const std::string& name,
    const geo_basic<T,M>& gamma,
    const quadrature_option& qopt)
  : _X(X), 
    _Y(Y), 
    _uu(),
    _ub(),
    _bu(),
    _bb()
{
    // example:
    //    form m (V,V,"mass",gamma);  e.g. \int_\Gamma trace(u) trace(v) ds     
    // with:
    //    geo omega ("square");
    //    geo gamma = omega["boundary"];
    //    V = space(omega,"P1");
    form_init_on_domain (name, gamma, false, field_basic<T,M>(), gamma, qopt);
}
template<class T, class M>
form_basic<T,M>::form_basic (
    const space_type& X, 
    const space_type& Y,
    const std::string& name,
    const geo_basic<T,M>& gamma,
    const field_basic<T,M>& wh,
    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");
    check_macro (
        wh.get_geo().get_background_geo().name()    == gamma.get_background_geo().name()
     && wh.get_geo().get_background_domain().name() == gamma.get_background_domain().name(),
        "form on domain \""
        << gamma.get_background_domain().name() << "\" of \""
        << gamma.get_background_geo().name()
        << "\" has incompatible weight, defined on \""
        << wh.get_geo().get_background_domain().name() << "\" of \""
        << wh.get_geo().get_background_geo().name() << "\"");
    form_init_on_domain (name, gamma, true, details::field_expr_v2_nonlinear_terminal_field<T,M,details::differentiate_option::none>(wh), wh.get_geo(), qopt);
}
// ----------------------------------------------------------------------------
// blas2
// ----------------------------------------------------------------------------
template<class T, class M>
field_basic<T,M>
form_basic<T,M>::operator* (const field_basic<T,M>& xh) const
{
    // TODO: verif des tailles des espaces ET de tous les vecteurs
    // si pas les memes cl, on pourrait iterer sur la form... + complique
    field_basic<T,M> yh (_Y, T(0));
    yh.set_u() = _uu*xh.u() + _ub*xh.b();
    yh.set_b() = _bu*xh.u() + _bb*xh.b();
    return yh;
}
template<class T, class M>
field_basic<T,M>
form_basic<T,M>::trans_mult (const field_basic<T,M>& x) const
{
    field_basic<T,M> y(get_first_space(), Float(0));
    y.set_u() = _uu.trans_mult(x.u()) + _bu.trans_mult(x.b());
    y.set_b() = _ub.trans_mult(x.u()) + _bb.trans_mult(x.b());
    return y;
}
template<class T, class M>
typename form_basic<T,M>::float_type
form_basic<T,M>::operator() (const field_basic<T,M>& uh, const field_basic<T,M>& vh) const
{
    return dual (operator*(uh), vh);
}
// ----------------------------------------------------------------------------
// blas3
// ----------------------------------------------------------------------------
template<class T, class M>
form_basic<T,M>
trans (const form_basic<T,M>& a)
{
  form_basic<T,M> b(a.get_second_space(), a.get_first_space());
  b.set_uu() = trans(a.uu());
  b.set_ub() = trans(a.bu()); // remark: may swap bu & ub
  b.set_bu() = trans(a.ub());
  b.set_bb() = trans(a.bb());
  return b;
}
// ----------------------------------------------------------------------------
// output: print all four csr as a large sparse matrix in matrix-market format
// ----------------------------------------------------------------------------
 
struct id {
  size_t operator() (size_t i) { return i; }
};
template<class T, class Permutation1, class Permutation2>
static
void
merge (
    asr<T,sequential>&       a, 
    const csr<T,sequential>& m,
    Permutation1             dis_im2dis_idof,
    Permutation2             dis_jm2dis_jdof)
{
    typedef typename form_basic<T,sequential>::size_type size_type;
    size_type i0 = m.row_ownership().first_index();
    size_type j0 = m.col_ownership().first_index();
    typename csr<T,sequential>::const_iterator ia = m.begin(); 
    for (size_type im = 0, nrow = m.nrow(); im < nrow; im++) {
      size_type dis_im = im + i0;
      size_type dis_idof = dis_im2dis_idof (dis_im);
      for (typename csr<T,sequential>::const_data_iterator p = ia[im]; p != ia[im+1]; p++) {
        const size_type& jm  = (*p).first;
        const T&         val = (*p).second;
        size_type dis_jm     = jm + j0;
        size_type dis_jdof   = dis_jm2dis_jdof (dis_jm);
        a.dis_entry (dis_idof, dis_jdof) += val;
      }
    }
}
#ifdef _RHEOLEF_HAVE_MPI
template<class T, class Permutation1, class Permutation2>
static
void
merge (
    asr<T,distributed>&       a, 
    const csr<T,distributed>& m,
    Permutation1              dis_im2dis_idof,
    Permutation2              dis_jm2dis_jdof)
{
    typedef typename form_basic<T,distributed>::size_type size_type;
    size_type i0 = m.row_ownership().first_index();
    size_type j0 = m.col_ownership().first_index();
    typename csr<T,distributed>::const_iterator ia = m.begin(); 
    for (size_type im = 0, nrow = m.nrow(); im < nrow; im++) {
      size_type dis_im = im + i0;
      size_type dis_idof = dis_im2dis_idof (dis_im);
      for (typename csr<T,distributed>::const_data_iterator p = ia[im]; p != ia[im+1]; p++) {
        const size_type& jm  = (*p).first;
        const T&         val = (*p).second;
        size_type dis_jm     = jm + j0;
        size_type dis_jdof   = dis_jm2dis_jdof (dis_jm);
        a.dis_entry (dis_idof, dis_jdof) += val;
      }
    }
    typename csr<T,distributed>::const_iterator ext_ia = m.ext_begin(); 
    for (size_type im = 0, nrow = m.nrow(); im < nrow; im++) {
      size_type dis_im = im + i0;
      size_type dis_idof = dis_im2dis_idof (dis_im);
      long int ext_size_im = std::distance(ext_ia[im],ext_ia[im+1]);
      for (typename csr<T,distributed>::const_data_iterator p = ext_ia[im]; p != ext_ia[im+1]; p++) {
        const size_type& jext = (*p).first;
        const T&         val  = (*p).second;
        size_type dis_jm      = m.jext2dis_j (jext);
        size_type dis_jdof    = dis_jm2dis_jdof (dis_jm);
        a.dis_entry (dis_idof, dis_jdof) += val;
      }
    }
}
#endif // _RHEOLEF_HAVE_MPI
template<class T, class M>
odiststream& 
form_basic<T,M>::put (odiststream& ops, bool show_partition) const
{
    // put all on io_proc 
    size_type dis_nrow = get_second_space().dis_ndof();
    size_type dis_ncol =  get_first_space().dis_ndof();
    size_type io_proc = odiststream::io_proc();
    size_type my_proc = comm().rank();
    distributor io_row_ownership (dis_nrow, comm(), (my_proc == io_proc ? dis_nrow : 0));
    distributor io_col_ownership (dis_ncol, comm(), (my_proc == io_proc ? dis_ncol : 0));
    asr<T,M> a (io_row_ownership, io_col_ownership);
 
    if (show_partition) {
        merge (a, _uu, id(), id());
        merge (a, _ub, id(), id());
        merge (a, _bu, id(), id());
        merge (a, _bb, id(), id());
    } else {
        error_macro ("not yet");
    }
    a.dis_entry_assembly();
    ops << "%%MatrixMarket matrix coordinate real general" << std::endl
        << dis_nrow << " " << dis_ncol << " " << a.dis_nnz() << std::endl
        << a;
    return ops;
}
template <class T, class M>
void
form_basic<T,M>::dump (std::string name) const
{
    _uu.dump (name + "-uu");
    _ub.dump (name + "-ub");
    _bu.dump (name + "-bu");
    _bb.dump (name + "-bb");
}
// ----------------------------------------------------------------------------
// diagonal part
// ----------------------------------------------------------------------------
template<class T, class M>
form_basic<T,M>
diag (const field_basic<T,M>& dh)
{
  form_basic<T,M> a (dh.get_space(), dh.get_space());
  a.set_uu() = diag(dh.u());
  a.set_bb() = diag(dh.b());
  return a;
}
template<class T, class M>
field_basic<T,M>
diag (const form_basic<T,M>& a)
{
  check_macro (a.get_first_space() == a.get_second_space(),
        "diag(form): incompatible first space "<<a.get_first_space().name()
                    << " and second one "<<a.get_second_space().name());
  field_basic<T,M> dh (a.get_first_space());
  dh.set_u() = vec<T,M>(diag(a.uu()));
  dh.set_b() = vec<T,M>(diag(a.bb()));
  return dh;
}
// ----------------------------------------------------------------------------
// instanciation in library
// ----------------------------------------------------------------------------
#define _RHEOLEF_instanciate(T,M)                                       \
template class form_basic<T,M>;                                         \
template class form_basic<T,M> trans (const form_basic<T,M>&);          \
template class field_basic<T,M> diag (const form_basic<T,M>&);          \
template class form_basic<T,M>  diag (const field_basic<T,M>&);
 
_RHEOLEF_instanciate(Float,sequential)
#ifdef _RHEOLEF_HAVE_MPI
_RHEOLEF_instanciate(Float,distributed)
#endif // _RHEOLEF_HAVE_MPI
 
}// namespace rheolef