Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
form.cc
Go to the documentation of this file.
1
21//
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
30namespace rheolef {
31using namespace std;
32
33template<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}
48template<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}
64template<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");
84 form_init_on_domain (name, gamma, false, field_basic<T,M>(), gamma, qopt);
85}
86template<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");
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// ----------------------------------------------------------------------------
121template<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}
132template<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}
141template<class T, class M>
144{
145 return dual (operator*(uh), vh);
146}
147// ----------------------------------------------------------------------------
148// blas3
149// ----------------------------------------------------------------------------
150template<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
165struct id {
166 size_t operator() (size_t i) { return i; }
167};
168template<class T, class Permutation1, class Permutation2>
169static
170void
171merge (
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
194template<class T, class Permutation1, class Permutation2>
195static
196void
197merge (
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
233template<class T, class M>
235form_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();
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}
260template <class T, class M>
261void
262form_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// ----------------------------------------------------------------------------
272template<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}
281template<class T, class M>
282field_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) \
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>&);
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
see the csr page for the full documentation
Definition csr.h:317
see the distributor page for the full documentation
Definition distributor.h:69
vec< T, M > & set_b()
Definition field.h:285
const vec< T, M > & b() const
Definition field.h:283
vec< T, M > & set_u()
Definition field.h:284
const geo_type & get_geo() const
Definition field.h:271
const space_type & get_space() const
Definition field.h:270
const vec< T, M > & u() const
Definition field.h:282
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)
generic mesh with rerefence counting
Definition geo.h:1089
odiststream: see the diststream page for the full documentation
Definition diststream.h:137
static size_type io_proc()
Definition diststream.cc:79
see the vec page for the full documentation
Definition vec.h:79
#define _RHEOLEF_instanciate(T)
Definition csr_seq.cc:491
#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.
rheolef::std enable_if ::type dual const Expr1 expr1, const Expr2 expr2 dual(const Expr1 &expr1, const Expr2 &expr2)
Definition field_expr.h:229
csr< T, M > diag(const vec< T, M > &d)
Definition csr.cc:56
csr< T, sequential > trans(const csr< T, sequential > &a)
trans(a): see the form page for the full documentation
Definition csr.h:455
STL namespace.
Definition sphere.icc:25