Rheolef  7.2
an efficient C++ finite element environment
form.cc
Go to the documentation of this file.
1 //
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"
29 
30 namespace rheolef {
31 using namespace std;
32 
33 template<class T, class M>
35  const space_type& X,
36  const space_type& Y,
37  const std::string& name,
38  const quadrature_option& qopt)
39  : _X(X),
40  _Y(Y),
41  _uu(),
42  _ub(),
43  _bu(),
44  _bb()
45 {
46  form_init (name, false, field_basic<T,M>(), qopt);
47 }
48 template<class T, class M>
50  const space_type& X,
51  const space_type& Y,
52  const std::string& name,
53  const field_basic<T,M>& wh,
54  const quadrature_option& qopt)
55  : _X(X),
56  _Y(Y),
57  _uu(),
58  _ub(),
59  _bu(),
60  _bb()
61 {
63 }
64 template<class T, class M>
66  const space_type& X,
67  const space_type& Y,
68  const std::string& name,
69  const geo_basic<T,M>& gamma,
70  const quadrature_option& qopt)
71  : _X(X),
72  _Y(Y),
73  _uu(),
74  _ub(),
75  _bu(),
76  _bb()
77 {
78  // example:
79  // form m (V,V,"mass",gamma); e.g. \int_\Gamma trace(u) trace(v) ds
80  // with:
81  // geo omega ("square");
82  // geo gamma = omega["boundary"];
83  // V = space(omega,"P1");
85 }
86 template<class T, class M>
88  const space_type& X,
89  const space_type& Y,
90  const std::string& name,
91  const geo_basic<T,M>& gamma,
92  const field_basic<T,M>& wh,
93  const quadrature_option& qopt)
94  : _X(X),
95  _Y(Y),
96  _uu(),
97  _ub(),
98  _bu(),
99  _bb()
100 {
101  // example:
102  // form m (V,V,"mass",gamma, weight); e.g. \int_\Gamma trace(u) trace(v) weight(x) ds
103  // with:
104  // geo omega ("square");
105  // geo gamma = omega["boundary"];
106  // V = space(omega,"P1");
107  check_macro (
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(),
110  "form on domain \""
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() << "\"");
117 }
118 // ----------------------------------------------------------------------------
119 // blas2
120 // ----------------------------------------------------------------------------
121 template<class T, class M>
124 {
125  // TODO: verif des tailles des espaces ET de tous les vecteurs
126  // si pas les memes cl, on pourrait iterer sur la form... + complique
127  field_basic<T,M> yh (_Y, T(0));
128  yh.set_u() = _uu*xh.u() + _ub*xh.b();
129  yh.set_b() = _bu*xh.u() + _bb*xh.b();
130  return yh;
131 }
132 template<class T, class M>
135 {
136  field_basic<T,M> y(get_first_space(), Float(0));
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());
139  return y;
140 }
141 template<class T, class M>
144 {
145  return dual (operator*(uh), vh);
146 }
147 // ----------------------------------------------------------------------------
148 // blas3
149 // ----------------------------------------------------------------------------
150 template<class T, class M>
153 {
154  form_basic<T,M> b(a.get_second_space(), a.get_first_space());
155  b.set_uu() = trans(a.uu());
156  b.set_ub() = trans(a.bu()); // remark: may swap bu & ub
157  b.set_bu() = trans(a.ub());
158  b.set_bb() = trans(a.bb());
159  return b;
160 }
161 // ----------------------------------------------------------------------------
162 // output: print all four csr as a large sparse matrix in matrix-market format
163 // ----------------------------------------------------------------------------
164 
165 struct id {
166  size_t operator() (size_t i) { return i; }
167 };
168 template<class T, class Permutation1, class Permutation2>
169 static
170 void
171 merge (
172  asr<T,sequential>& a,
173  const csr<T,sequential>& m,
174  Permutation1 dis_im2dis_idof,
175  Permutation2 dis_jm2dis_jdof)
176 {
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++) {
182  size_type dis_im = im + i0;
183  size_type dis_idof = dis_im2dis_idof (dis_im);
184  for (typename csr<T,sequential>::const_data_iterator p = ia[im]; p != ia[im+1]; p++) {
185  const size_type& jm = (*p).first;
186  const T& val = (*p).second;
187  size_type dis_jm = jm + j0;
188  size_type dis_jdof = dis_jm2dis_jdof (dis_jm);
189  a.dis_entry (dis_idof, dis_jdof) += val;
190  }
191  }
192 }
193 #ifdef _RHEOLEF_HAVE_MPI
194 template<class T, class Permutation1, class Permutation2>
195 static
196 void
197 merge (
200  Permutation1 dis_im2dis_idof,
201  Permutation2 dis_jm2dis_jdof)
202 {
204  size_type i0 = m.row_ownership().first_index();
205  size_type j0 = m.col_ownership().first_index();
206  typename csr<T,distributed>::const_iterator ia = m.begin();
207  for (size_type im = 0, nrow = m.nrow(); im < nrow; im++) {
208  size_type dis_im = im + i0;
209  size_type dis_idof = dis_im2dis_idof (dis_im);
210  for (typename csr<T,distributed>::const_data_iterator p = ia[im]; p != ia[im+1]; p++) {
211  const size_type& jm = (*p).first;
212  const T& val = (*p).second;
213  size_type dis_jm = jm + j0;
214  size_type dis_jdof = dis_jm2dis_jdof (dis_jm);
215  a.dis_entry (dis_idof, dis_jdof) += val;
216  }
217  }
218  typename csr<T,distributed>::const_iterator ext_ia = m.ext_begin();
219  for (size_type im = 0, nrow = m.nrow(); im < nrow; im++) {
220  size_type dis_im = im + i0;
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]);
223  for (typename csr<T,distributed>::const_data_iterator p = ext_ia[im]; p != ext_ia[im+1]; p++) {
224  const size_type& jext = (*p).first;
225  const T& val = (*p).second;
226  size_type dis_jm = m.jext2dis_j (jext);
227  size_type dis_jdof = dis_jm2dis_jdof (dis_jm);
228  a.dis_entry (dis_idof, dis_jdof) += val;
229  }
230  }
231 }
232 #endif // _RHEOLEF_HAVE_MPI
233 template<class T, class M>
234 odiststream&
235 form_basic<T,M>::put (odiststream& ops, bool show_partition) const
236 {
237  // put all on io_proc
238  size_type dis_nrow = get_second_space().dis_ndof();
239  size_type dis_ncol = get_first_space().dis_ndof();
240  size_type io_proc = odiststream::io_proc();
241  size_type my_proc = comm().rank();
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);
245 
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());
251  } else {
252  error_macro ("not yet");
253  }
254  a.dis_entry_assembly();
255  ops << "%%MatrixMarket matrix coordinate real general" << std::endl
256  << dis_nrow << " " << dis_ncol << " " << a.dis_nnz() << std::endl
257  << a;
258  return ops;
259 }
260 template <class T, class M>
261 void
262 form_basic<T,M>::dump (std::string name) const
263 {
264  _uu.dump (name + "-uu");
265  _ub.dump (name + "-ub");
266  _bu.dump (name + "-bu");
267  _bb.dump (name + "-bb");
268 }
269 // ----------------------------------------------------------------------------
270 // diagonal part
271 // ----------------------------------------------------------------------------
272 template<class T, class M>
275 {
276  form_basic<T,M> a (dh.get_space(), dh.get_space());
277  a.set_uu() = diag(dh.u());
278  a.set_bb() = diag(dh.b());
279  return a;
280 }
281 template<class T, class M>
282 field_basic<T,M>
284 {
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());
288  field_basic<T,M> dh (a.get_first_space());
289  dh.set_u() = vec<T,M>(diag(a.uu()));
290  dh.set_b() = vec<T,M>(diag(a.bb()));
291  return dh;
292 }
293 // ----------------------------------------------------------------------------
294 // instanciation in library
295 // ----------------------------------------------------------------------------
296 #define _RHEOLEF_instanciate(T,M) \
297 template class form_basic<T,M>; \
298 template class form_basic<T,M> trans (const form_basic<T,M>&); \
299 template class field_basic<T,M> diag (const form_basic<T,M>&); \
300 template class form_basic<T,M> diag (const field_basic<T,M>&);
301 
303 #ifdef _RHEOLEF_HAVE_MPI
305 #endif // _RHEOLEF_HAVE_MPI
306 
307 }// namespace rheolef
field::size_type size_type
Definition: branch.cc:430
see the Float page for the full documentation
rep::const_data_iterator const_data_iterator
Definition: csr.h:336
rep::const_iterator const_iterator
Definition: csr.h:334
see the csr page for the full documentation
Definition: csr.h:317
see the distributor page for the full documentation
Definition: distributor.h:69
const space_type & get_space() const
Definition: field.h:270
vec< T, M > & set_b()
Definition: field.h:285
const geo_type & get_geo() const
Definition: field.h:271
vec< T, M > & set_u()
Definition: field.h:284
const vec< T, M > & u() const
Definition: field.h:282
const vec< T, M > & b() const
Definition: field.h:283
scalar_traits< T >::type float_type
Definition: form.h:152
form_basic< T, M > operator*(const form_basic< T, M > &b) const
Definition: form.h:397
odiststream & put(odiststream &ops, bool show_partition=true) const
Definition: form.cc:235
float_type operator()(const field_basic< T, M > &uh, const field_basic< T, M > &vh) const
Definition: form.cc:143
csr< T, M >::size_type size_type
Definition: form.h:150
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)
field_basic< T, M > trans_mult(const field_basic< T, M > &yh) const
Definition: form.cc:134
void dump(std::string name) const
Definition: form.cc:262
void form_init(const std::string &name, bool has_weight, WeightFunction weight, const quadrature_option &qopt)
odiststream: see the diststream page for the full documentation
Definition: diststream.h:137
static size_type io_proc()
Definition: diststream.cc:79
integrate_option quadrature_option
size_t size_type
Definition: basis_get.cc:76
double Float
see the Float page for the full documentation
Definition: Float.h:143
#define error_macro(message)
Definition: dis_macros.h:49
Expr1::float_type T
Definition: field_expr.h:230
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.
t operator()(const t &a, const t &b)
Definition: space.cc:386
_RHEOLEF_instanciate(Float, sequential) _RHEOLEF_instanciate(Float
csr< T, sequential > trans(const csr< T, sequential > &a)
trans(a): see the form page for the full documentation
Definition: csr.h:455
csr< T, M > diag(const vec< T, M > &d)
Definition: csr.cc:56
rheolef::std enable_if ::type dual const Expr1 expr1, const Expr2 expr2 dual(const Expr1 &expr1, const Expr2 &expr2)
Definition: field_expr.h:229
Float gamma[][pmax+1]
Definition: sphere.icc:25