Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
form_concat.cc
Go to the documentation of this file.
1
21// build form from initializer list (c++ 2011)
22//
23#include "rheolef/form_concat.h"
24#include "rheolef/space_mult.h"
25#include "rheolef/csr_concat.h"
26
27namespace rheolef { namespace details {
28
29// =========================================================================
30// 1rst case : one-line matrix initializer
31// A = {a, b}; // matrix & vector
32// =========================================================================
33
34template <class T, class M>
35void
36form_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
37{
38 // ------------------------------------------------------------
39 // pass 0 : lazy first space computation, compute second space
40 // ------------------------------------------------------------
42 size_t j_comp = 0;
43 bool have_Yh = false;
44 typename std::vector<std::pair<bool,space_basic<T,M> > >::iterator xh_iter = l_Xh.begin();
45 for (typename std::list<value_type>::const_iterator iter = _l.begin(); iter != _l.end(); ++iter, xh_iter++, j_comp++) {
46 const value_type& x = *iter;
47 switch (x.variant) {
49 check_macro (x.s == 0, "unsupported non-nul scalar `"<<x.s<<"' in form concatenation"
50 << " at ("<<i_comp<<","<<j_comp<<")");
51 break;
52 }
54 if (!(*xh_iter).first) {
55 (*xh_iter).first = true;
56 (*xh_iter).second = IR;
57 } else {
58 check_macro (IR == (*xh_iter).second, "form initializer: invalid field:"
59 << " expect a form with its second space `" << (*xh_iter).second.name() << "'"
60 << " at ("<<i_comp<<","<<j_comp<<")");
61 }
62 if (!have_Yh) {
63 have_Yh = true;
64 Yh = x.v.get_space();
65 } else {
66 check_macro (x.v.get_space() == Yh, "form initializer: invalid field space `"
67 << x.v.get_space().name() << "': expect `" << Yh.name() << "'"
68 << " at ("<<i_comp<<","<<j_comp<<")");
69 }
70 break;
71 }
73 if (!(*xh_iter).first) {
74 (*xh_iter).first = true;
75 (*xh_iter).second = x.v.get_space();
76 } else {
77 check_macro (x.v.get_space() == (*xh_iter).second, "form initializer: invalid trans(field) space `"
78 << x.v.get_space().name() << "': expect `" << (*xh_iter).second.name() << "'"
79 << " at ("<<i_comp<<","<<j_comp<<")");
80 }
81 if (!have_Yh) {
82 have_Yh = true;
83 Yh = IR;
84 } else {
85 check_macro (IR == Yh, "form initializer: invalid trans(field): "
86 << " expect a form with its second space `"
87 << Yh.name() << "'" << " at ("<<i_comp<<","<<j_comp<<")");
88 }
89 break;
90 }
92 size_t n = x.vv.size();
93 check_macro (n != 0, "form initializer: invalid vector<field> with zero-size");
94 space_basic<T,M> IRn = pow(IR,n); // or space::real(n) ?
95 if (!(*xh_iter).first) {
96 (*xh_iter).first = true;
97 (*xh_iter).second = x.vv[0].get_space();
98 } else {
99 check_macro (x.vv[0].get_space() == (*xh_iter).second, "form initializer: invalid 0-th vector<field> space `"
100 << x.vv[0].get_space().name() << "': expect `" << (*xh_iter).second.name() << "'"
101 << " at ("<<i_comp<<","<<j_comp<<")");
102 }
103 for (size_t i = 1; i < n; ++i) {
104 check_macro (x.vv[i].get_space() == x.vv[0].get_space(), "form initializer: invalid "<<i<<"-th vector<field> space `"
105 << x.vv[i].get_space().name() << "': expect `" << (*xh_iter).second.name() << "'"
106 << " at ("<<i_comp<<","<<j_comp<<")");
107 }
108 if (!have_Yh) {
109 have_Yh = true;
110 Yh = IRn;
111 } else {
112 check_macro (IRn == Yh, "form initializer: invalid vector<field>: "
113 << " expect a form with its second space `"
114 << Yh.name() << "'" << " at ("<<i_comp<<","<<j_comp<<")");
115 }
116 break;
117 }
119 size_t n = x.vv.size();
120 check_macro (n != 0, "form initializer: invalid trans(vector<field>) with zero-size");
121 space_basic<T,M> IRn = pow(IR,n); // or space::real(n) ? faster...
122 if (!(*xh_iter).first) {
123 (*xh_iter).first = true;
124 (*xh_iter).second = IRn;
125 } else {
126 check_macro (IRn == (*xh_iter).second, "form initializer: invalid trans(vector<field>):"
127 << " expect a form with its second space `" << (*xh_iter).second.name() << "'"
128 << " at ("<<i_comp<<","<<j_comp<<")");
129 }
130 if (!have_Yh) {
131 have_Yh = true;
132 Yh = x.vv[0].get_space();
133 } else {
134 check_macro (x.vv[0].get_space() == Yh, "form initializer: invalid 0-th trans(vector<field>) space `"
135 << x.vv[0].get_space().name() << "': expect `" << Yh.name() << "'"
136 << " at ("<<i_comp<<","<<j_comp<<")");
137 }
138 for (size_t i = 1; i < n; ++i) {
139 check_macro (x.vv[i].get_space() == x.vv[0].get_space(), "form initializer: invalid "<<i<<"-th trans(vector<field>) space `"
140 << x.vv[i].get_space().name() << "': expect `" << (*xh_iter).second.name() << "'"
141 << " at ("<<i_comp<<","<<j_comp<<")");
142 }
143 break;
144 }
146 if (!(*xh_iter).first) {
147 (*xh_iter).first = true;
148 (*xh_iter).second = x.m.get_first_space();
149 } else {
150 check_macro (x.m.get_first_space() == (*xh_iter).second, "form initializer: invalid second space `"
151 << x.m.get_first_space().name() << "': expect `" << (*xh_iter).second.name() << "'"
152 << " at ("<<i_comp<<","<<j_comp<<")");
153 }
154 if (!have_Yh) {
155 have_Yh = true;
156 Yh = x.m.get_second_space();
157 } else {
158 check_macro (x.m.get_second_space() == Yh, "form initializer: invalid second space `"
159 << x.m.get_second_space().name() << "': expect `" << Yh.name() << "'"
160 << " at ("<<i_comp<<","<<j_comp<<")");
161 }
162 break;
163 }
164 default: error_macro ("concatenation not yet supported"
165 << " at ("<<i_comp<<","<<j_comp<<")");
166 }
167 }
168 check_macro (have_Yh, "form concatenation: "<<i_comp<<"th row space remains undefined");
169}
170template <class T, class M>
171void
172form_concat_line<T,M>::build_first_space (const std::vector<std::pair<bool,space_basic<T,M> > >& l_Xh, space_basic<T,M>& Xh)
173{
174 // ------------------------------------------------------------
175 // pass 0b : first space computation
176 // ------------------------------------------------------------
178 size_t j_comp = 0;
179 for (typename std::vector<std::pair<bool,space_basic<T,M> > >::const_iterator
180 xh_iter = l_Xh.begin(),
181 xh_last = l_Xh.end(); xh_iter != xh_last; xh_iter++, j_comp++) {
182 check_macro ((*xh_iter).first, "form concatenation: "<<j_comp<<"th column space remains undefined");
183 sml_X *= (*xh_iter).second;
184 }
185 Xh = space_basic<T,M>(sml_X);
186}
187template <class T, class M>
188void
190{
191 // --------------------------------
192 // pass 1 : both spaces computation
193 // --------------------------------
194 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>()));
195 build_form_pass0 (l_Xh, Yh);
196 build_first_space (l_Xh, Xh);
197}
198template <class T, class M>
201{
202 // -----------------------
203 // pass 2 : compute values
204 // -----------------------
205 form_basic<T,M> a (Xh, Yh);
206 csr_concat_line<T,M> uu, ub, bu, bb;
207 for(typename std::list<value_type>::const_iterator iter = _l.begin(); iter != _l.end(); ++iter) {
208 const value_type& x = *iter;
209 switch (x.variant) {
211 uu.push_back (x.m.uu());
212 ub.push_back (x.m.ub());
213 bu.push_back (x.m.bu());
214 bb.push_back (x.m.bb());
215 break;
216 }
217 default: error_macro ("non-form concatenation not yet supported");
218 }
219 }
220 a.set_uu() = uu.build_csr();
221 a.set_ub() = ub.build_csr();
222 a.set_bu() = bu.build_csr();
223 a.set_bb() = bb.build_csr();
224 return a;
225}
226template <class T, class M>
229{
230 space_basic<T,M> Xh, Yh;
231 build_form_pass1 (Xh, Yh);
232 form_basic<T,M> a = build_form_pass2 (Xh, Yh);
233 return a;
234}
235// =========================================================================
236// 2nd case : multi-line form initializer
237// A = { {a, b },
238// {c, d} };
239// =========================================================================
240template <class T, class M>
243{
244 // ---------------------------
245 // pass 1 : compute spaces
246 // ---------------------------
247 size_t i_comp = 0;
249 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>()));
250 for (typename std::list<line_type>::const_iterator iter = _l.begin(); iter != _l.end(); ++iter, i_comp++) {
251 const line_type& line = *iter;
253 line.build_form_pass0 (l_Xh, Yih, i_comp);
254 sml_Y *= Yih;
255 }
258 space_basic<T,M> Yh (sml_Y);
259 // ------------------------
260 // pass 2 : copy
261 // ------------------------
262 csr_concat<T,M> uu, ub, bu, bb;
263 sizes_type undefs (undef,undef);
264 sizes_type zeros (zero, zero);
265 typedef sizes_pair_type empty; // for clarity
266 for (typename std::list<line_type>::const_iterator iter = _l.begin(); iter != _l.end(); ++iter) {
267 const line_type& line = *iter;
268 csr_concat_line<T,M> uu_i, ub_i, bu_i, bb_i;
269 for (typename std::list<value_type>::const_iterator jter = line.begin(); jter != line.end(); ++jter) {
270 const value_type& x = *jter;
271 switch (x.variant) {
273 check_macro (x.s == 0, "unsupported non-nul scalar `"<<x.s<<"' in form concatenation");
274 // scalar zero: goes to uu block, nothing on others ub bu bb
275 uu_i.push_back (x.s); // zero block with size=(undef,undef)
276 ub_i.push_back (empty(undefs, undefs));
277 bu_i.push_back (empty(undefs, undefs));
278 bb_i.push_back (empty(undefs, undefs));
279 break;
280 }
282 // vertical vec
283 sizes_type u_sizes (x.v.u().size(), x.v.u().dis_size());
284 sizes_type b_sizes (x.v.b().size(), x.v.b().dis_size());
285 uu_i.push_back (x.v.u());
286 ub_i.push_back (empty(u_sizes, zeros));
287 bu_i.push_back (x.v.b());
288 bb_i.push_back (empty(b_sizes, zeros));
289 break;
290 }
292 // horizontal vec
293 sizes_type u_sizes (x.v.u().size(), x.v.u().dis_size());
294 sizes_type b_sizes (x.v.b().size(), x.v.b().dis_size());
295 uu_i.push_back (trans(x.v.u()));
296 ub_i.push_back (trans(x.v.b()));
297 bu_i.push_back (empty(zeros, u_sizes));
298 bb_i.push_back (empty(zeros, b_sizes));
299 break;
300 }
302 // horizontal vector<vec>
303 size_t n = x.vv.size();
304#ifdef TO_CLEAN
305 sizes_type u_sizes (zero, zero);
306 sizes_type b_sizes (zero, zero);
307 std::vector<vec<T,M>> x_vv_u(n), x_vv_b(n);
308 for (size_t i = 0; i < n; ++i) {
309 x_vv_u[i] = x.vv[i].u();
310 x_vv_b[i] = x.vv[i].b();
311 u_sizes.first += x.vv[i].u().size();
312 u_sizes.second += x.vv[i].u().dis_size();
313 b_sizes.first += x.vv[i].b().size();
314 b_sizes.second += x.vv[i].b().dis_size();
315 }
316#endif // TO_CLEAN
317 sizes_type u_sizes, b_sizes;
318 if (n == 0) {
319 u_sizes = sizes_type (undef, undef);
320 b_sizes = sizes_type (undef, undef);
321 } else {
322 u_sizes.first = x.vv[0].u().size();
323 u_sizes.second = x.vv[0].u().dis_size();
324 b_sizes.first = x.vv[0].b().size();
325 b_sizes.second = x.vv[0].b().dis_size();
326 for (size_t i = 0; i < n; ++i) {
327 check_macro (u_sizes.first == x.vv[i].u().size() &&
328 u_sizes.second == x.vv[i].u().dis_size(),
329 "trans(vector<field>): "<<i<<"-th vector.u size [0:"<<x.vv[i].u().size()<<"|"<<x.vv[i].u().dis_size()<<"]"
330 " is incompatible with 0-th one [" <<x.vv[0].u().size()<<"|"<<x.vv[0].u().dis_size()<<"]");
331 check_macro (b_sizes.first == x.vv[i].b().size() &&
332 b_sizes.second == x.vv[i].b().dis_size(),
333 "trans(vector<field>): "<<i<<"-th vector.b size [0:"<<x.vv[i].b().size()<<"|"<<x.vv[i].b().dis_size()<<"]"
334 " is incompatible with 0-th one [" <<x.vv[0].b().size()<<"|"<<x.vv[0].b().dis_size()<<"]");
335 }
336 }
337 std::vector<vec<T,M>> x_vv_u(n), x_vv_b(n);
338 for (size_t i = 0; i < n; ++i) {
339 x_vv_u[i] = x.vv[i].u();
340 x_vv_b[i] = x.vv[i].b();
341 }
342 uu_i.push_back (x_vv_u);
343 ub_i.push_back (x_vv_b);
344 bu_i.push_back (empty(zeros, u_sizes));
345 bb_i.push_back (empty(zeros, b_sizes));
346 break;
347 }
349 // vertical vector<vec>
350 size_t n = x.vv.size();
351 sizes_type u_sizes, b_sizes;
352 if (n == 0) {
353 u_sizes = sizes_type (undef, undef);
354 b_sizes = sizes_type (undef, undef);
355 } else {
356 u_sizes.first = x.vv[0].u().size();
357 u_sizes.second = x.vv[0].u().dis_size();
358 b_sizes.first = x.vv[0].b().size();
359 b_sizes.second = x.vv[0].b().dis_size();
360 for (size_t i = 0; i < n; ++i) {
361 check_macro (u_sizes.first == x.vv[i].u().size() &&
362 u_sizes.second == x.vv[i].u().dis_size(),
363 "trans(vector<field>): "<<i<<"-th vector.u size [0:"<<x.vv[i].u().size()<<"|"<<x.vv[i].u().dis_size()<<"]"
364 " is incompatible with 0-th one [" <<x.vv[0].u().size()<<"|"<<x.vv[0].u().dis_size()<<"]");
365 check_macro (b_sizes.first == x.vv[i].b().size() &&
366 b_sizes.second == x.vv[i].b().dis_size(),
367 "trans(vector<field>): "<<i<<"-th vector.b size [0:"<<x.vv[i].b().size()<<"|"<<x.vv[i].b().dis_size()<<"]"
368 " is incompatible with 0-th one [" <<x.vv[0].b().size()<<"|"<<x.vv[0].b().dis_size()<<"]");
369 }
370 }
371 std::vector<vec<T,M>> x_vv_u(n), x_vv_b(n);
372 for (size_t i = 0; i < n; ++i) {
373 x_vv_u[i] = x.vv[i].u();
374 x_vv_b[i] = x.vv[i].b();
375 }
376 uu_i.push_back (trans(x_vv_u));
377 ub_i.push_back (empty(u_sizes, zeros));
378 bu_i.push_back (trans(x_vv_b));
379 bb_i.push_back (empty(b_sizes, zeros));
380 break;
381 }
383 uu_i.push_back (x.m.uu());
384 ub_i.push_back (x.m.ub());
385 bu_i.push_back (x.m.bu());
386 bb_i.push_back (x.m.bb());
387 break;
388 }
389 default: error_macro ("non-form or scalar concatenation not yet supported");
390 }
391 }
392 uu.push_back (uu_i);
393 ub.push_back (ub_i);
394 bu.push_back (bu_i);
395 bb.push_back (bb_i);
396 }
397 form_basic<T,M> a(Xh, Yh);
398 a.set_uu() = uu.build_csr();
399 a.set_ub() = ub.build_csr();
400 a.set_bu() = bu.build_csr();
401 a.set_bb() = bb.build_csr();
402 return a;
403}
404// ----------------------------------------------------------------------------
405// instanciation in library
406// ----------------------------------------------------------------------------
407#define _RHEOLEF_instanciation(T,M) \
408template class form_concat_line<T,M>; \
409template class form_concat<T,M>;
410
412#ifdef _RHEOLEF_HAVE_MPI
413_RHEOLEF_instanciation(Float,distributed)
414#endif // _RHEOLEF_HAVE_MPI
415
416}} // namespace rheolef::details
#define _RHEOLEF_instanciation(T, M, A)
Definition asr.cc:223
see the Float page for the full documentation
void push_back(const value_type &x)
Definition csr_concat.h:117
csr< T, M > build_csr() const
void push_back(const line_type &line)
Definition csr_concat.h:200
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
form_basic< T, M > build_form_pass2(const space_basic< T, M > &Xh, const space_basic< T, M > &Yh) const
form_basic< T, M > build_form() const
void build_form_pass1(space_basic< T, M > &Xh, space_basic< T, M > &Yh) const
form_concat_value< T, M > value_type
static void build_first_space(const std::vector< std::pair< bool, space_basic< T, M > > > &l_Xh, space_basic< T, M > &Xh)
form_concat_line< T, M > line_type
csr_concat_value< T, M >::sizes_type sizes_type
form_basic< T, M > build_form() const
csr_concat_value< T, M >::sizes_pair_type sizes_pair_type
form_concat_value< T, M > value_type
const space_type & get_space() const
Definition field.h:270
the finite element space
Definition space.h:382
#define error_macro(message)
Definition dis_macros.h:49
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
This file is part of Rheolef.
space_mult_list< T, M > pow(const space_basic< T, M > &X, size_t n)
Definition space_mult.h:120
csr< T, sequential > trans(const csr< T, sequential > &a)
trans(a): see the form page for the full documentation
Definition csr.h:455