22#include "rheolef/field_expr.h"
23#include "rheolef/geo_domain.h"
24#include "rheolef/form.h"
25#include "rheolef/csr.h"
26#include "rheolef/field.h"
27#include "rheolef/band.h"
28#include "rheolef/form_weighted.h"
33template<
class T,
class M>
37 const std::string& name,
38 const quadrature_option& qopt)
48template<
class T,
class M>
52 const std::string& name,
54 const quadrature_option& qopt)
64template<
class T,
class M>
68 const std::string& name,
70 const quadrature_option& qopt)
86template<
class T,
class M>
90 const std::string& name,
93 const quadrature_option& qopt)
108 wh.
get_geo().get_background_geo().name() == gamma.get_background_geo().name()
109 && wh.
get_geo().get_background_domain().name() == gamma.get_background_domain().name(),
111 << gamma.get_background_domain().name() <<
"\" of \""
112 << gamma.get_background_geo().name()
113 <<
"\" has incompatible weight, defined on \""
114 << wh.
get_geo().get_background_domain().name() <<
"\" of \""
115 << wh.
get_geo().get_background_geo().name() <<
"\"");
121template<
class T,
class M>
128 yh.
set_u() = _uu*xh.
u() + _ub*xh.
b();
129 yh.
set_b() = _bu*xh.
u() + _bb*xh.
b();
132template<
class T,
class M>
137 y.
set_u() = _uu.trans_mult(x.
u()) + _bu.trans_mult(x.
b());
138 y.
set_b() = _ub.trans_mult(x.
u()) + _bb.trans_mult(x.
b());
141template<
class T,
class M>
145 return dual (
operator*(uh), vh);
150template<
class T,
class M>
155 b.set_uu() =
trans(a.uu());
156 b.set_ub() =
trans(a.bu());
157 b.set_bu() =
trans(a.ub());
158 b.set_bb() =
trans(a.bb());
166 size_t operator() (
size_t i) {
return i; }
168template<
class T,
class Permutation1,
class Permutation2>
172 asr<T,sequential>& a,
173 const csr<T,sequential>& m,
174 Permutation1 dis_im2dis_idof,
175 Permutation2 dis_jm2dis_jdof)
178 size_type i0 =
m.row_ownership().first_index();
179 size_type j0 =
m.col_ownership().first_index();
180 typename csr<T,sequential>::const_iterator ia =
m.begin();
181 for (
size_type im = 0, nrow =
m.nrow(); im < nrow; im++) {
184 for (
typename csr<T,sequential>::const_data_iterator
p = ia[im];
p != ia[im+1];
p++) {
186 const T& val = (*p).second;
188 size_type dis_jdof = dis_jm2dis_jdof (dis_jm);
189 a.dis_entry (dis_idof, dis_jdof) += val;
193#ifdef _RHEOLEF_HAVE_MPI
194template<
class T,
class Permutation1,
class Permutation2>
200 Permutation1 dis_im2dis_idof,
201 Permutation2 dis_jm2dis_jdof)
204 size_type i0 = m.row_ownership().first_index();
205 size_type j0 = m.col_ownership().first_index();
207 for (
size_type im = 0, nrow = m.nrow(); im < nrow; im++) {
209 size_type dis_idof = dis_im2dis_idof (dis_im);
212 const T& val = (*p).second;
214 size_type dis_jdof = dis_jm2dis_jdof (dis_jm);
215 a.dis_entry (dis_idof, dis_jdof) += val;
219 for (
size_type im = 0, nrow = m.nrow(); im < nrow; im++) {
221 size_type dis_idof = dis_im2dis_idof (dis_im);
222 long int ext_size_im = std::distance(ext_ia[im],ext_ia[im+1]);
225 const T& val = (*p).second;
227 size_type dis_jdof = dis_jm2dis_jdof (dis_jm);
228 a.dis_entry (dis_idof, dis_jdof) += val;
233template<
class T,
class M>
238 size_type dis_nrow = get_second_space().dis_ndof();
239 size_type dis_ncol = get_first_space().dis_ndof();
242 distributor io_row_ownership (dis_nrow, comm(), (my_proc == io_proc ? dis_nrow : 0));
243 distributor io_col_ownership (dis_ncol, comm(), (my_proc == io_proc ? dis_ncol : 0));
244 asr<T,M> a (io_row_ownership, io_col_ownership);
246 if (show_partition) {
247 merge (a, _uu,
id(),
id());
248 merge (a, _ub,
id(),
id());
249 merge (a, _bu,
id(),
id());
250 merge (a, _bb,
id(),
id());
254 a.dis_entry_assembly();
255 ops <<
"%%MatrixMarket matrix coordinate real general" << std::endl
256 << dis_nrow <<
" " << dis_ncol <<
" " << a.dis_nnz() << std::endl
260template <
class T,
class M>
264 _uu.dump (name +
"-uu");
265 _ub.dump (name +
"-ub");
266 _bu.dump (name +
"-bu");
267 _bb.dump (name +
"-bb");
272template<
class T,
class M>
277 a.set_uu() =
diag(dh.
u());
278 a.set_bb() =
diag(dh.
b());
281template<
class T,
class M>
285 check_macro (a.get_first_space() == a.get_second_space(),
286 "diag(form): incompatible first space "<<a.get_first_space().name()
287 <<
" and second one "<<a.get_second_space().name());
296#define _RHEOLEF_instanciate(T,M) \
297template class form_basic<T,M>; \
298template class form_basic<T,M> trans (const form_basic<T,M>&); \
299template class field_basic<T,M> diag (const form_basic<T,M>&); \
300template class form_basic<T,M> diag (const field_basic<T,M>&);
303#ifdef _RHEOLEF_HAVE_MPI
field::size_type size_type
see the Float page for the full documentation
see the csr page for the full documentation
see the distributor page for the full documentation
const vec< T, M > & b() const
const geo_type & get_geo() const
const space_type & get_space() const
const vec< T, M > & u() const
generic mesh with rerefence counting
odiststream: see the diststream page for the full documentation
static size_type io_proc()
see the vec page for the full documentation
#define _RHEOLEF_instanciate(T)
#define error_macro(message)
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
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)
This file is part of Rheolef.
rheolef::std enable_if ::type dual const Expr1 expr1, const Expr2 expr2 dual(const Expr1 &expr1, const Expr2 &expr2)
csr< T, M > diag(const vec< T, M > &d)
csr< T, sequential > trans(const csr< T, sequential > &a)
trans(a): see the form page for the full documentation