Loading...
Searching...
No Matches
form_concat.cc
Go to the documentation of this file.
// build form from initializer list (c++ 2011)
//
#include "rheolef/form_concat.h"
#include "rheolef/space_mult.h"
#include "rheolef/csr_concat.h"
namespace rheolef { namespace details {
// =========================================================================
// 1rst case : one-line matrix initializer
// A = {a, b}; // matrix & vector
// =========================================================================
template <class T, class M>
void
form_concat_line<T,M>::build_form_pass0 (std::vector<std::pair<bool,space_basic<T,M> > >& l_Xh, space_basic<T,M>& Yh, size_t i_comp) const
{
// ------------------------------------------------------------
// pass 0 : lazy first space computation, compute second space
// ------------------------------------------------------------
space_basic<T,M> IR = space_basic<T,M>::real();
size_t j_comp = 0;
bool have_Yh = false;
typename std::vector<std::pair<bool,space_basic<T,M> > >::iterator xh_iter = l_Xh.begin();
for (typename std::list<value_type>::const_iterator iter = _l.begin(); iter != _l.end(); ++iter, xh_iter++, j_comp++) {
const value_type& x = *iter;
switch (x.variant) {
case form_concat_value<T,M>::scalar: {
check_macro (x.s == 0, "unsupported non-nul scalar `"<<x.s<<"' in form concatenation"
<< " at ("<<i_comp<<","<<j_comp<<")");
break;
}
case form_concat_value<T,M>::field: {
if (!(*xh_iter).first) {
(*xh_iter).first = true;
(*xh_iter).second = IR;
} else {
check_macro (IR == (*xh_iter).second, "form initializer: invalid field:"
<< " expect a form with its second space `" << (*xh_iter).second.name() << "'"
<< " at ("<<i_comp<<","<<j_comp<<")");
}
if (!have_Yh) {
have_Yh = true;
Yh = x.v.get_space();
} else {
check_macro (x.v.get_space() == Yh, "form initializer: invalid field space `"
<< x.v.get_space().name() << "': expect `" << Yh.name() << "'"
<< " at ("<<i_comp<<","<<j_comp<<")");
}
break;
}
case form_concat_value<T,M>::field_transpose: {
if (!(*xh_iter).first) {
(*xh_iter).first = true;
(*xh_iter).second = x.v.get_space();
} else {
check_macro (x.v.get_space() == (*xh_iter).second, "form initializer: invalid trans(field) space `"
<< x.v.get_space().name() << "': expect `" << (*xh_iter).second.name() << "'"
<< " at ("<<i_comp<<","<<j_comp<<")");
}
if (!have_Yh) {
have_Yh = true;
Yh = IR;
} else {
check_macro (IR == Yh, "form initializer: invalid trans(field): "
<< " expect a form with its second space `"
<< Yh.name() << "'" << " at ("<<i_comp<<","<<j_comp<<")");
}
break;
}
case form_concat_value<T,M>::vector_field: {
size_t n = x.vv.size();
check_macro (n != 0, "form initializer: invalid vector<field> with zero-size");
space_basic<T,M> IRn = pow(IR,n); // or space::real(n) ?
if (!(*xh_iter).first) {
(*xh_iter).first = true;
(*xh_iter).second = x.vv[0].get_space();
} else {
check_macro (x.vv[0].get_space() == (*xh_iter).second, "form initializer: invalid 0-th vector<field> space `"
<< x.vv[0].get_space().name() << "': expect `" << (*xh_iter).second.name() << "'"
<< " at ("<<i_comp<<","<<j_comp<<")");
}
for (size_t i = 1; i < n; ++i) {
check_macro (x.vv[i].get_space() == x.vv[0].get_space(), "form initializer: invalid "<<i<<"-th vector<field> space `"
<< x.vv[i].get_space().name() << "': expect `" << (*xh_iter).second.name() << "'"
<< " at ("<<i_comp<<","<<j_comp<<")");
}
if (!have_Yh) {
have_Yh = true;
Yh = IRn;
} else {
check_macro (IRn == Yh, "form initializer: invalid vector<field>: "
<< " expect a form with its second space `"
<< Yh.name() << "'" << " at ("<<i_comp<<","<<j_comp<<")");
}
break;
}
case form_concat_value<T,M>::vector_field_transpose: {
size_t n = x.vv.size();
check_macro (n != 0, "form initializer: invalid trans(vector<field>) with zero-size");
space_basic<T,M> IRn = pow(IR,n); // or space::real(n) ? faster...
if (!(*xh_iter).first) {
(*xh_iter).first = true;
(*xh_iter).second = IRn;
} else {
check_macro (IRn == (*xh_iter).second, "form initializer: invalid trans(vector<field>):"
<< " expect a form with its second space `" << (*xh_iter).second.name() << "'"
<< " at ("<<i_comp<<","<<j_comp<<")");
}
if (!have_Yh) {
have_Yh = true;
Yh = x.vv[0].get_space();
} else {
check_macro (x.vv[0].get_space() == Yh, "form initializer: invalid 0-th trans(vector<field>) space `"
<< x.vv[0].get_space().name() << "': expect `" << Yh.name() << "'"
<< " at ("<<i_comp<<","<<j_comp<<")");
}
for (size_t i = 1; i < n; ++i) {
check_macro (x.vv[i].get_space() == x.vv[0].get_space(), "form initializer: invalid "<<i<<"-th trans(vector<field>) space `"
<< x.vv[i].get_space().name() << "': expect `" << (*xh_iter).second.name() << "'"
<< " at ("<<i_comp<<","<<j_comp<<")");
}
break;
}
case form_concat_value<T,M>::form: {
if (!(*xh_iter).first) {
(*xh_iter).first = true;
(*xh_iter).second = x.m.get_first_space();
} else {
check_macro (x.m.get_first_space() == (*xh_iter).second, "form initializer: invalid second space `"
<< x.m.get_first_space().name() << "': expect `" << (*xh_iter).second.name() << "'"
<< " at ("<<i_comp<<","<<j_comp<<")");
}
if (!have_Yh) {
have_Yh = true;
Yh = x.m.get_second_space();
} else {
check_macro (x.m.get_second_space() == Yh, "form initializer: invalid second space `"
<< x.m.get_second_space().name() << "': expect `" << Yh.name() << "'"
<< " at ("<<i_comp<<","<<j_comp<<")");
}
break;
}
default: error_macro ("concatenation not yet supported"
<< " at ("<<i_comp<<","<<j_comp<<")");
}
}
check_macro (have_Yh, "form concatenation: "<<i_comp<<"th row space remains undefined");
}
template <class T, class M>
void
form_concat_line<T,M>::build_first_space (const std::vector<std::pair<bool,space_basic<T,M> > >& l_Xh, space_basic<T,M>& Xh)
{
// ------------------------------------------------------------
// pass 0b : first space computation
// ------------------------------------------------------------
space_mult_list<T,M> sml_X;
size_t j_comp = 0;
for (typename std::vector<std::pair<bool,space_basic<T,M> > >::const_iterator
xh_iter = l_Xh.begin(),
xh_last = l_Xh.end(); xh_iter != xh_last; xh_iter++, j_comp++) {
check_macro ((*xh_iter).first, "form concatenation: "<<j_comp<<"th column space remains undefined");
sml_X *= (*xh_iter).second;
}
Xh = space_basic<T,M>(sml_X);
}
template <class T, class M>
void
form_concat_line<T,M>::build_form_pass1 (space_basic<T,M>& Xh, space_basic<T,M>& Yh) const
{
// --------------------------------
// pass 1 : both spaces computation
// --------------------------------
std::vector<std::pair<bool,space_basic<T,M> > > l_Xh (_l.size(), std::pair<bool,space_basic<T,M> >(false, space_basic<T,M>()));
build_form_pass0 (l_Xh, Yh);
build_first_space (l_Xh, Xh);
}
template <class T, class M>
form_basic<T,M>
form_concat_line<T,M>::build_form_pass2 (const space_basic<T,M>& Xh, const space_basic<T,M>& Yh) const
{
// -----------------------
// pass 2 : compute values
// -----------------------
form_basic<T,M> a (Xh, Yh);
csr_concat_line<T,M> uu, ub, bu, bb;
for(typename std::list<value_type>::const_iterator iter = _l.begin(); iter != _l.end(); ++iter) {
const value_type& x = *iter;
switch (x.variant) {
case form_concat_value<T,M>::form: {
uu.push_back (x.m.uu());
ub.push_back (x.m.ub());
bu.push_back (x.m.bu());
bb.push_back (x.m.bb());
break;
}
default: error_macro ("non-form concatenation not yet supported");
}
}
a.set_uu() = uu.build_csr();
a.set_ub() = ub.build_csr();
a.set_bu() = bu.build_csr();
a.set_bb() = bb.build_csr();
return a;
}
template <class T, class M>
form_basic<T,M>
form_concat_line<T,M>::build_form() const
{
space_basic<T,M> Xh, Yh;
build_form_pass1 (Xh, Yh);
form_basic<T,M> a = build_form_pass2 (Xh, Yh);
return a;
}
// =========================================================================
// 2nd case : multi-line form initializer
// A = { {a, b },
// {c, d} };
// =========================================================================
template <class T, class M>
form_basic<T,M>
form_concat<T,M>::build_form() const
{
// ---------------------------
// pass 1 : compute spaces
// ---------------------------
size_t i_comp = 0;
space_mult_list<T,M> sml_Y;
std::vector<std::pair<bool,space_basic<T,M> > > l_Xh (_l.size(), std::pair<bool,space_basic<T,M> >(false, space_basic<T,M>()));
for (typename std::list<line_type>::const_iterator iter = _l.begin(); iter != _l.end(); ++iter, i_comp++) {
const line_type& line = *iter;
space_basic<T,M> Yih;
line.build_form_pass0 (l_Xh, Yih, i_comp);
sml_Y *= Yih;
}
space_basic<T,M> Xh;
form_concat_line<T,M>::build_first_space (l_Xh, Xh);
space_basic<T,M> Yh (sml_Y);
// ------------------------
// pass 2 : copy
// ------------------------
csr_concat<T,M> uu, ub, bu, bb;
sizes_type undefs (undef,undef);
sizes_type zeros (zero, zero);
typedef sizes_pair_type empty; // for clarity
for (typename std::list<line_type>::const_iterator iter = _l.begin(); iter != _l.end(); ++iter) {
const line_type& line = *iter;
csr_concat_line<T,M> uu_i, ub_i, bu_i, bb_i;
for (typename std::list<value_type>::const_iterator jter = line.begin(); jter != line.end(); ++jter) {
const value_type& x = *jter;
switch (x.variant) {
case form_concat_value<T,M>::scalar: {
check_macro (x.s == 0, "unsupported non-nul scalar `"<<x.s<<"' in form concatenation");
// scalar zero: goes to uu block, nothing on others ub bu bb
uu_i.push_back (x.s); // zero block with size=(undef,undef)
ub_i.push_back (empty(undefs, undefs));
bu_i.push_back (empty(undefs, undefs));
bb_i.push_back (empty(undefs, undefs));
break;
}
case form_concat_value<T,M>::field: {
// vertical vec
sizes_type u_sizes (x.v.u().size(), x.v.u().dis_size());
sizes_type b_sizes (x.v.b().size(), x.v.b().dis_size());
uu_i.push_back (x.v.u());
ub_i.push_back (empty(u_sizes, zeros));
bu_i.push_back (x.v.b());
bb_i.push_back (empty(b_sizes, zeros));
break;
}
case form_concat_value<T,M>::field_transpose: {
// horizontal vec
sizes_type u_sizes (x.v.u().size(), x.v.u().dis_size());
sizes_type b_sizes (x.v.b().size(), x.v.b().dis_size());
uu_i.push_back (trans(x.v.u()));
ub_i.push_back (trans(x.v.b()));
bu_i.push_back (empty(zeros, u_sizes));
bb_i.push_back (empty(zeros, b_sizes));
break;
}
case form_concat_value<T,M>::vector_field: {
// horizontal vector<vec>
size_t n = x.vv.size();
#ifdef TO_CLEAN
sizes_type u_sizes (zero, zero);
sizes_type b_sizes (zero, zero);
std::vector<vec<T,M>> x_vv_u(n), x_vv_b(n);
for (size_t i = 0; i < n; ++i) {
x_vv_u[i] = x.vv[i].u();
x_vv_b[i] = x.vv[i].b();
u_sizes.first += x.vv[i].u().size();
u_sizes.second += x.vv[i].u().dis_size();
b_sizes.first += x.vv[i].b().size();
b_sizes.second += x.vv[i].b().dis_size();
}
#endif // TO_CLEAN
sizes_type u_sizes, b_sizes;
if (n == 0) {
u_sizes = sizes_type (undef, undef);
b_sizes = sizes_type (undef, undef);
} else {
u_sizes.first = x.vv[0].u().size();
u_sizes.second = x.vv[0].u().dis_size();
b_sizes.first = x.vv[0].b().size();
b_sizes.second = x.vv[0].b().dis_size();
for (size_t i = 0; i < n; ++i) {
check_macro (u_sizes.first == x.vv[i].u().size() &&
u_sizes.second == x.vv[i].u().dis_size(),
"trans(vector<field>): "<<i<<"-th vector.u size [0:"<<x.vv[i].u().size()<<"|"<<x.vv[i].u().dis_size()<<"]"
" is incompatible with 0-th one [" <<x.vv[0].u().size()<<"|"<<x.vv[0].u().dis_size()<<"]");
check_macro (b_sizes.first == x.vv[i].b().size() &&
b_sizes.second == x.vv[i].b().dis_size(),
"trans(vector<field>): "<<i<<"-th vector.b size [0:"<<x.vv[i].b().size()<<"|"<<x.vv[i].b().dis_size()<<"]"
" is incompatible with 0-th one [" <<x.vv[0].b().size()<<"|"<<x.vv[0].b().dis_size()<<"]");
}
}
std::vector<vec<T,M>> x_vv_u(n), x_vv_b(n);
for (size_t i = 0; i < n; ++i) {
x_vv_u[i] = x.vv[i].u();
x_vv_b[i] = x.vv[i].b();
}
uu_i.push_back (x_vv_u);
ub_i.push_back (x_vv_b);
bu_i.push_back (empty(zeros, u_sizes));
bb_i.push_back (empty(zeros, b_sizes));
break;
}
case form_concat_value<T,M>::vector_field_transpose: {
// vertical vector<vec>
size_t n = x.vv.size();
sizes_type u_sizes, b_sizes;
if (n == 0) {
u_sizes = sizes_type (undef, undef);
b_sizes = sizes_type (undef, undef);
} else {
u_sizes.first = x.vv[0].u().size();
u_sizes.second = x.vv[0].u().dis_size();
b_sizes.first = x.vv[0].b().size();
b_sizes.second = x.vv[0].b().dis_size();
for (size_t i = 0; i < n; ++i) {
check_macro (u_sizes.first == x.vv[i].u().size() &&
u_sizes.second == x.vv[i].u().dis_size(),
"trans(vector<field>): "<<i<<"-th vector.u size [0:"<<x.vv[i].u().size()<<"|"<<x.vv[i].u().dis_size()<<"]"
" is incompatible with 0-th one [" <<x.vv[0].u().size()<<"|"<<x.vv[0].u().dis_size()<<"]");
check_macro (b_sizes.first == x.vv[i].b().size() &&
b_sizes.second == x.vv[i].b().dis_size(),
"trans(vector<field>): "<<i<<"-th vector.b size [0:"<<x.vv[i].b().size()<<"|"<<x.vv[i].b().dis_size()<<"]"
" is incompatible with 0-th one [" <<x.vv[0].b().size()<<"|"<<x.vv[0].b().dis_size()<<"]");
}
}
std::vector<vec<T,M>> x_vv_u(n), x_vv_b(n);
for (size_t i = 0; i < n; ++i) {
x_vv_u[i] = x.vv[i].u();
x_vv_b[i] = x.vv[i].b();
}
uu_i.push_back (trans(x_vv_u));
ub_i.push_back (empty(u_sizes, zeros));
bu_i.push_back (trans(x_vv_b));
bb_i.push_back (empty(b_sizes, zeros));
break;
}
case form_concat_value<T,M>::form: {
uu_i.push_back (x.m.uu());
ub_i.push_back (x.m.ub());
bu_i.push_back (x.m.bu());
bb_i.push_back (x.m.bb());
break;
}
default: error_macro ("non-form or scalar concatenation not yet supported");
}
}
uu.push_back (uu_i);
ub.push_back (ub_i);
bu.push_back (bu_i);
bb.push_back (bb_i);
}
form_basic<T,M> a(Xh, Yh);
a.set_uu() = uu.build_csr();
a.set_ub() = ub.build_csr();
a.set_bu() = bu.build_csr();
a.set_bb() = bb.build_csr();
return a;
}
// ----------------------------------------------------------------------------
// instanciation in library
// ----------------------------------------------------------------------------
#define _RHEOLEF_instanciation(T,M) \
template class form_concat_line<T,M>; \
template class form_concat<T,M>;
_RHEOLEF_instanciation(Float,sequential)
#ifdef _RHEOLEF_HAVE_MPI
_RHEOLEF_instanciation(Float,distributed)
#endif // _RHEOLEF_HAVE_MPI
}} // namespace rheolef::details
_RHEOLEF_instanciation
#define _RHEOLEF_instanciation(T, M, A)
Definition asr.cc:223
Float
see the Float page for the full documentation
rheolef::details::csr_concat_line
Definition csr_concat.h:93
rheolef::details::csr_concat_line::build_csr
csr< T, M > build_csr() const
Definition csr_concat.cc:353
rheolef::details::csr_concat_line::push_back
void push_back(const value_type &x)
Definition csr_concat.h:117
rheolef::details::csr_concat
Definition csr_concat.h:175
rheolef::details::csr_concat::build_csr
csr< T, M > build_csr() const
Definition csr_concat.cc:378
rheolef::details::csr_concat::push_back
void push_back(const line_type &line)
Definition csr_concat.h:200
rheolef::details::form_concat_line::build_form_pass0
void build_form_pass0(std::vector< std::pair< bool, space_basic< T, M > > > &l_Xh, space_basic< T, M > &Yh, size_t i_comp=0) const
Definition form_concat.cc:36
rheolef::details::form_concat_line::begin
const_iterator begin() const
Definition form_concat.h:145
rheolef::details::form_concat_line::_l
std::list< value_type > _l
Definition form_concat.h:165
rheolef::details::form_concat_line::build_form_pass2
form_basic< T, M > build_form_pass2(const space_basic< T, M > &Xh, const space_basic< T, M > &Yh) const
Definition form_concat.cc:200
rheolef::details::form_concat_line::build_form
form_basic< T, M > build_form() const
Definition form_concat.cc:228
rheolef::details::form_concat_line::build_form_pass1
void build_form_pass1(space_basic< T, M > &Xh, space_basic< T, M > &Yh) const
Definition form_concat.cc:189
rheolef::details::form_concat_line::end
const_iterator end() const
Definition form_concat.h:146
rheolef::details::form_concat_line::value_type
form_concat_value< T, M > value_type
Definition form_concat.h:129
rheolef::details::form_concat_line::const_iterator
std::list< value_type >::const_iterator const_iterator
Definition form_concat.h:130
rheolef::details::form_concat_line::build_first_space
static void build_first_space(const std::vector< std::pair< bool, space_basic< T, M > > > &l_Xh, space_basic< T, M > &Xh)
Definition form_concat.cc:172
rheolef::details::form_concat_value::v
field_basic< T, M > v
Definition form_concat.h:111
rheolef::details::form_concat_value::vv
std::vector< field_basic< T, M > > vv
Definition form_concat.h:112
rheolef::details::form_concat_value::m
form_basic< T, M > m
Definition form_concat.h:113
rheolef::details::form_concat_value::field_transpose
@ field_transpose
Definition form_concat.h:71
rheolef::details::form_concat_value::field
@ field
Definition form_concat.h:70
rheolef::details::form_concat_value::vector_field_transpose
@ vector_field_transpose
Definition form_concat.h:73
rheolef::details::form_concat_value::scalar
@ scalar
Definition form_concat.h:69
rheolef::details::form_concat_value::vector_field
@ vector_field
Definition form_concat.h:72
rheolef::details::form_concat_value::form
@ form
Definition form_concat.h:74
rheolef::details::form_concat_value::s
T s
Definition form_concat.h:110
rheolef::details::form_concat_value::variant
variant_type variant
Definition form_concat.h:114
rheolef::details::form_concat::zero
static constexpr size_type zero
Definition form_concat.h:199
rheolef::details::form_concat::line_type
form_concat_line< T, M > line_type
Definition form_concat.h:201
rheolef::details::form_concat::sizes_type
csr_concat_value< T, M >::sizes_type sizes_type
Definition form_concat.h:195
rheolef::details::form_concat::build_form
form_basic< T, M > build_form() const
Definition form_concat.cc:242
rheolef::details::form_concat::_l
std::list< line_type > _l
Definition form_concat.h:226
rheolef::details::form_concat::sizes_pair_type
csr_concat_value< T, M >::sizes_pair_type sizes_pair_type
Definition form_concat.h:196
rheolef::details::form_concat::value_type
form_concat_value< T, M > value_type
Definition form_concat.h:202
rheolef::details::form_concat::undef
static constexpr size_type undef
Definition form_concat.h:198
rheolef::field_basic::b
const vec< T, M > & b() const
Definition field.h:283
rheolef::field_basic::get_space
const space_type & get_space() const
Definition field.h:270
rheolef::field_basic::u
const vec< T, M > & u() const
Definition field.h:282
rheolef::form_basic
Definition form.h:146
rheolef::space_basic
the finite element space
Definition space.h:382
rheolef::space_mult_list
Definition space_mult.h:54
error_macro
#define error_macro(message)
Definition dis_macros.h:49
check_macro
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
rheolef::details
Definition compiler_eigen.h:39
rheolef
This file is part of Rheolef.
Definition compiler_eigen.h:39
rheolef::pow
space_mult_list< T, M > pow(const space_basic< T, M > &X, size_t n)
Definition space_mult.h:120
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