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
size_type
field::size_type size_type
Definition: branch.cc:425
Float
see the Float page for the full documentation
rheolef::asr
Definition: asr.h:49
rheolef::csr< T, distributed >::const_data_iterator
rep::const_data_iterator const_data_iterator
Definition: csr.h:479
rheolef::csr< T, distributed >::const_iterator
rep::const_iterator const_iterator
Definition: csr.h:477
rheolef::csr< T, sequential >::const_iterator
rep::const_iterator const_iterator
Definition: csr.h:334
rheolef::csr
see the csr page for the full documentation
Definition: csr.h:317
rheolef::details::field_expr_v2_nonlinear_terminal_field< T, M, details::differentiate_option::none >
rheolef::distributor
see the distributor page for the full documentation
Definition: distributor.h:62
rheolef::field_basic< T, M >
rheolef::field_basic::get_space
const space_type & get_space() const
Definition: field.h:300
rheolef::field_basic::set_b
vec< T, M > & set_b()
Definition: field.h:313
rheolef::field_basic::get_geo
const geo_type & get_geo() const
Definition: field.h:301
rheolef::field_basic::set_u
vec< T, M > & set_u()
Definition: field.h:312
rheolef::field_basic::u
const vec< T, M > & u() const
Definition: field.h:310
rheolef::field_basic::b
const vec< T, M > & b() const
Definition: field.h:311
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::operator*
form_basic< T, M > operator*(const form_basic< T, M > &b) const
Definition: form.h:381
rheolef::form_basic::put
odiststream & put(odiststream &ops, bool show_partition=true) const
Definition: form.cc:235
rheolef::form_basic::operator()
float_type operator()(const field_basic< T, M > &uh, const field_basic< T, M > &vh) const
Definition: form.cc:143
rheolef::form_basic::form_basic
form_basic()
Definition: form.h:300
rheolef::form_basic::size_type
csr< T, M >::size_type size_type
Definition: form.h:148
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::trans_mult
field_basic< T, M > trans_mult(const field_basic< T, M > &yh) const
Definition: form.cc:134
rheolef::form_basic::dump
void dump(std::string name) const
Definition: form.cc:262
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_basic< T, M >
rheolef::odiststream
odiststream: see the diststream page for the full documentation
Definition: diststream.h:126
rheolef::odiststream::io_proc
static size_type io_proc()
Definition: diststream.cc:78
rheolef::space_basic< float_type, M >
rheolef::vec< T, M >
rheolef::quadrature_option
integrate_option quadrature_option
Definition: integrate_option.h:186
rheolef::size_type
size_t size_type
Definition: basis_get.cc:76
rheolef::Float
double Float
see the Float page for the full documentation
Definition: Float.h:143
error_macro
#define error_macro(message)
Definition: dis_macros.h:49
T
Expr1::float_type T
Definition: field_expr.h:261
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.b
int b
Definition: mkgeo_ball.sh:152
mkgeo_ball.a
int a
Definition: mkgeo_ball.sh:151
mkgeo_contraction.name
string name
Definition: mkgeo_contraction.sh:133
mkgeo_sector.m
m
Definition: mkgeo_sector.sh:118
rheolef::space_numbering::dis_idof
void dis_idof(const basis_basic< T > &b, const geo_size &gs, const geo_element &K, typename std::vector< size_type >::iterator dis_idof_tab)
Definition: space_numbering.cc:187
rheolef
This file is part of Rheolef.
Definition: compiler_eigen.h:37
rheolef::_RHEOLEF_instanciate
_RHEOLEF_instanciate(Float, sequential) _RHEOLEF_instanciate(Float
rheolef::trans
csr< T, sequential > trans(const csr< T, sequential > &a)
trans(a): see the form page for the full documentation
Definition: csr.h:455
rheolef::diag
csr< T, M > diag(const vec< T, M > &d)
Definition: csr.cc:56
rheolef::dual
rheolef::std enable_if ::type dual const Expr1 expr1, const Expr2 expr2 dual(const Expr1 &expr1, const Expr2 &expr2)
Definition: field_expr.h:260
rk::gamma
Float gamma[][pmax+1]
Definition: runge_kutta_semiimplicit.icc:70
a
Definition: diffusion_isotropic.h:25
p
Definition: sphere.icc:25
rheolef::distributed
Definition: distributed.h:88
rheolef::sequential
Definition: distributed.h:28