Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
space.cc
Go to the documentation of this file.
1
21
22#include "rheolef/space.h"
23#include "rheolef/space_numbering.h"
24#include "rheolef/space_mult.h"
25#include "rheolef/piola_util.h"
26#include "rheolef/piola_on_pointset.h"
27
28namespace rheolef {
29
30// ======================================================================
31// space_base_rep
32// ======================================================================
33template <class T, class M>
35 const geo_basic<T,M>& omega_in,
36 std::string approx,
37 std::string valued)
38: _constit(omega_in,
39 (approx == "" || valued == "scalar") ? approx : valued+"("+approx+")"),
40 _xdof(),
41 _have_freezed(false),
42 _idof2blk_dis_iub(),
43 _has_nt_basis(),
44 _normal(),
45 _iu_ownership(),
46 _ib_ownership(),
47 _parent_subgeo_owner_initialized(false),
48 _parent_subgeo_owner_reattribution_required(false),
49 _parent_subgeo_owner()
50{
51 // omega_in is compressed by constitution_terminal() allocator
52 // when it is a domain: then it becomes a geo_domain, with compressed numbering
53 // so, omega_in can be different from omega:
54 if (approx == "") return; // empty element => default cstor
55 _idof2blk_dis_iub.resize (_constit.ownership());
56 _has_nt_basis.resize (_constit.ownership()); // TODO: _constit.scalar_ownership()
57 _normal.resize (_constit.ownership());
58 init_xdof();
59}
60template <class T, class M>
61space_base_rep<T,M>::space_base_rep (
62 const geo_basic<T,M>& omega_in,
63 const basis_basic<T>& b)
64: _constit(omega_in, b.name()),
65 _xdof(),
66 _have_freezed(false),
67 _idof2blk_dis_iub(),
68 _has_nt_basis(),
69 _normal(),
70 _iu_ownership(),
71 _ib_ownership(),
72 _parent_subgeo_owner_initialized(false),
73 _parent_subgeo_owner_reattribution_required(false),
74 _parent_subgeo_owner()
75{
76 // omega_in is compressed by constitution_terminal() allocator
77 // when it is a domain: then it becomes a geo_domain, with compressed numbering
78 // so, omega_in can be different from omega:
79 if (! b.is_initialized()) return; // empty element => default cstor
80 _idof2blk_dis_iub.resize (_constit.ownership());
81 _has_nt_basis.resize (_constit.ownership());
82 _normal.resize (_constit.ownership());
83 init_xdof();
84}
85template <class T, class M>
86space_base_rep<T,M>::space_base_rep (const space_constitution<T,M>& constit)
87: _constit(constit),
88 _xdof(),
89 _have_freezed(false),
90 _idof2blk_dis_iub(),
91 _has_nt_basis(),
92 _normal(),
93 _iu_ownership(),
94 _ib_ownership(),
95 _parent_subgeo_owner_initialized(false),
96 _parent_subgeo_owner_reattribution_required(false),
97 _parent_subgeo_owner()
98{
99 distributor idof_ownership = _constit.ownership();
100 _idof2blk_dis_iub.resize (idof_ownership);
101 _has_nt_basis.resize (idof_ownership);
102 _normal.resize (idof_ownership);
103 init_xdof();
104}
105template <class T, class M>
106space_base_rep<T,M>::space_base_rep (const space_mult_list<T,M>& expr)
107: _constit(expr),
108 _xdof(),
109 _have_freezed(false),
110 _idof2blk_dis_iub(),
111 _has_nt_basis(),
112 _normal(),
113 _iu_ownership(),
114 _ib_ownership(),
115 _parent_subgeo_owner_initialized(false),
116 _parent_subgeo_owner_reattribution_required(false),
117 _parent_subgeo_owner()
118{
119 _idof2blk_dis_iub.resize (_constit.ownership());
120 _has_nt_basis.resize (_constit.ownership());
121 _normal.resize (_constit.ownership());
122 init_xdof();
123}
124template <class T, class M>
125void
126space_base_rep<T,M>::init_xdof()
127{
128 if (_constit.is_hierarchical()) {
129 trace_macro ("init_xdof: hierarchical space: xdof undefined");
130 return;
131 }
132 const geo_basic<T,M>& omega = get_geo();
133 const basis_basic<T>& b = get_basis();
134 size_type dis_nnod = space_numbering::dis_nnod (b, omega.sizes(), omega.map_dimension());
135 size_type nnod = space_numbering:: nnod (b, omega.sizes(), omega.map_dimension());
136 communicator comm = ownership().comm();
137 distributor nod_ownership (dis_nnod, comm, nnod);
138 _xdof.resize (nod_ownership);
139 piola_on_pointset<T> pops;
140 pops.initialize (omega.get_piola_basis(), b, integrate_option());
141 std::vector<size_type> dis_inod_tab;
142 for (size_type ie = 0, ne = omega.size(); ie < ne; ie++) {
143 const geo_element& K = omega [ie];
144 const Eigen::Matrix<piola<T>,Eigen::Dynamic,1>& piola = pops.get_piola (omega, K);
145 space_numbering::dis_inod (b, omega.sizes(), K, dis_inod_tab);
146 for (size_type loc_inod = 0, loc_nnod = piola.size(); loc_inod < loc_nnod; ++loc_inod) {
147 _xdof.dis_entry (dis_inod_tab[loc_inod]) = piola [loc_inod].F;
148 }
149 }
150 _xdof.dis_entry_assembly();
151 // add external nodes on current element:
152 if (is_distributed<M>::value && comm.size() > 1) {
153 index_set dis_inod_set;
154 for (size_type ie = 0, ne = omega.size(); ie < ne; ++ie) {
155 const geo_element& K = omega [ie];
156 space_numbering::dis_inod (b, omega.sizes(), K, dis_inod_tab);
157 for (size_type loc_inod = 0, loc_nnod = dis_inod_tab.size(); loc_inod < loc_nnod; ++loc_inod) {
158 size_type dis_inod = dis_inod_tab[loc_inod];
159 if (nod_ownership.is_owned(dis_inod)) continue;
160 check_macro (dis_inod < dis_nnod, "invalid dis_inod="<<dis_inod<<" out of range [0:"<<dis_nnod<<"[");
161 dis_inod_set.insert (dis_inod);
162 }
163 }
164 _xdof.set_dis_indexes (dis_inod_set);
165 }
166}
167template <class T, class M>
168void
169space_base_rep<T,M>::base_freeze_body () const
170{
171 _constit.neighbour_guard(); // geo could uses S.master(i), a neighbour connnectivity
172 // -----------------------------------------------------------------------
173 // 1) loop on domains: mark blocked dofs
174 // -----------------------------------------------------------------------
175 disarray<size_type,M> blocked_flag = _constit.build_blocked_flag();
176
177 // copy the blocked_flag into the dis_iub disarray, as the "blocked" bit:
178 for (size_type idof = 0, ndof = blocked_flag.size(); idof < ndof; idof++) {
179 _idof2blk_dis_iub [idof].set_blocked (blocked_flag[idof]);
180 }
181 // -----------------------------------------------------------------------
182 // 2) init numbering
183 // -----------------------------------------------------------------------
184 size_type n_unknown = 0;
185 size_type n_blocked = 0;
186 for (size_type idof = 0, ndof = _idof2blk_dis_iub.size(); idof < ndof; idof++) {
187 bool blk = _idof2blk_dis_iub [idof].is_blocked();
188 if (! blk) {
189 _idof2blk_dis_iub[idof].set_dis_iub (n_unknown);
190 n_unknown++;
191 } else {
192 _idof2blk_dis_iub[idof].set_dis_iub (n_blocked);
193 n_blocked++;
194 }
195 }
196 size_type dis_n_unknown = n_unknown;
197 size_type dis_n_blocked = n_blocked;
198#ifdef _RHEOLEF_HAVE_MPI
199 if (is_distributed<M>::value) {
200 dis_n_unknown = mpi::all_reduce (comm(), dis_n_unknown, std::plus<T>());
201 dis_n_blocked = mpi::all_reduce (comm(), dis_n_blocked, std::plus<T>());
202 }
203#endif // // _RHEOLEF_HAVE_MPI
204 _iu_ownership = distributor (dis_n_unknown, comm(), n_unknown);
205 _ib_ownership = distributor (dis_n_blocked, comm(), n_blocked);
206}
207// ----------------------------------------------------------------------------
208// still used by geo_subdivide.cc and internally in space.cc
209// ----------------------------------------------------------------------------
210template <class T, class M>
211void
212space_base_rep<T,M>::dis_idof (const geo_element& K, std::vector<size_type>& dis_idof) const
213{
214 freeze_guard();
215 _constit.assembly_dis_idof (get_geo(), K, dis_idof);
216}
217// ----------------------------------------------------------------------------
218// name : e.g. "P1(square)", for field_expr<Expr> checks
219// ----------------------------------------------------------------------------
220template <class T, class M>
221std::string
222space_base_rep<T,M>::name() const
223{
224 return _constit.name();
225}
226// ----------------------------------------------------------------------------
227// u["left"]
228// => build here the requiresd temporary indirect disarray
229// ----------------------------------------------------------------------------
239template <class T, class M>
240disarray<typename space_base_rep<T,M>::size_type, M>
241space_base_rep<T,M>::build_dom_dis_idof2bgd_dis_idof (
242 const space_base_rep<T,M>& Wh, // Wh = space on domain
243 const std::string& dom_name // redundant: contained in Wh.get_geo()
244 ) const
245{
246 const geo_basic<T,M>& bgd_gamma = get_geo()[dom_name];
247 return build_dom_dis_idof2bgd_dis_idof (Wh, bgd_gamma);
248}
249template <class T, class M>
250disarray<typename space_base_rep<T,M>::size_type, M>
251space_base_rep<T,M>::build_dom_dis_idof2bgd_dis_idof (
252 const space_base_rep<T,M>& Wh, // Wh = space on domain
253 const geo_basic<T,M>& bgd_gamma2
254 ) const
255{
256 // TODO: move it to the call ?
257 const geo_basic<T,M>& bgd_gamma = bgd_gamma2.get_background_domain();
258
259 // TODO: verifier que le domaine bgd_gamma est compatible:
260 // => il doit y avoir un meme maillage de base a bgd_gamma & Wh.get_geo
261 const space_base_rep<T,M>& Vh = *this; // Vh = space on background mesh
262 Vh.freeze_guard();
263 Wh.freeze_guard();
264 const geo_basic<T,M>& dom_gamma = Wh.get_geo();
265 size_type map_dim = dom_gamma.map_dimension();
266 check_macro (dom_gamma.size() == bgd_gamma.size(), "incompatible domains");
267 distributor dom_ownership = Wh.ownership();
268 distributor bgd_ownership = Vh.ownership();
269 size_type my_proc = dom_ownership.comm().rank();
270trace_macro("Vh="<<Vh.name()<<", Wh="<<Wh.name()<<", dom_dis_idof2bgd_dis_idof.size=Wh.size="<<dom_ownership.size()<<"...");
271 size_type first_dom_dis_idof = dom_ownership.first_index();
272 size_type first_bgd_dis_idof = bgd_ownership.first_index();
273 std::vector<size_type> dom_dis_idofs, bgd_dis_idofs;
274 disarray<size_type, M> dom_dis_idof2bgd_dis_idof (dom_ownership, std::numeric_limits<size_type>::max());
275 std::set<size_type> ext_dom_dis_idof;
276 for (size_type ige = 0, nge = dom_gamma.size(); ige < nge; ige++) {
277 const geo_element& dom_S = dom_gamma[ige];
278 const geo_element& bgd_S = bgd_gamma[ige];
279 Wh.dis_idof (dom_S, dom_dis_idofs);
280 Vh.dis_idof (bgd_S, bgd_dis_idofs);
281trace_macro("ige="<<ige<<": dom_S="<<dom_S.name()<<dom_S.dis_ie()<<", bgd_S="<<bgd_S.name()<<bgd_S.dis_ie()
282 << ": dom_dis_idofs.size="<<dom_dis_idofs.size()
283 << ", bgd_dis_idofs.size="<<bgd_dis_idofs.size());
284 for (size_type loc_idof = 0, loc_ndof = dom_dis_idofs.size(); loc_idof < loc_ndof; loc_idof++) {
285 size_type dom_dis_idof = dom_dis_idofs [loc_idof];
286 size_type bgd_dis_idof = bgd_dis_idofs [loc_idof];
287 dom_dis_idof2bgd_dis_idof.dis_entry (dom_dis_idof) = bgd_dis_idof;
288 trace_macro("ige="<<ige<<": dom_dis_idof2bgd_dis_idof["<<dom_dis_idof<<"] = "<<bgd_dis_idof);
289 if (! dom_ownership.is_owned (dom_dis_idof, my_proc)) {
290 ext_dom_dis_idof.insert (dom_dis_idof);
291 }
292 }
293 }
294 trace_macro ("ext_dom_dis_idof.size="<<ext_dom_dis_idof.size());
295 dom_dis_idof2bgd_dis_idof.dis_entry_assembly();
296 dom_dis_idof2bgd_dis_idof.set_dis_indexes (ext_dom_dis_idof); // doit etre en deuxieme...?
297#ifdef TO_CLEAN
298 // move to local numbering:
299 for (size_type dom_idof = 0, dom_ndof = dom_dis_idof2bgd_dis_idof.size(); dom_idof < dom_ndof; dom_idof++) {
300 size_type bgd_dis_idof = dom_dis_idof2bgd_dis_idof [dom_idof];
301 size_type dom_dis_idof = dom_idof + first_dom_dis_idof;
302 check_macro (bgd_ownership.is_owned (bgd_dis_idof), "bgd_dis_idof="<<bgd_dis_idof<<" is out of range ["<<first_bgd_dis_idof<<":"<< bgd_ownership.last_index()<<"[");
303 size_type bgd_idof = bgd_dis_idof - first_bgd_dis_idof;
304 dom_dis_idof2bgd_dis_idof [dom_idof] = bgd_idof;
305 }
306trace_macro("Vh="<<Vh.name()<<", Wh="<<Wh.name()<<", dom_dis_idof2bgd_dis_idof.size=Wh.size="<<dom_ownership.size()<<" done");
307#endif // TO_CLEAN
308 return dom_dis_idof2bgd_dis_idof;
309}
310#define _RHEOLEF_space_real(M) \
311template <class T> \
312space_basic<T,M> \
313space_basic<T,M>::real() \
314{ \
315 return space_basic<T,M> (geo_basic<T,M>(details::zero_dimension()), "P1"); \
316}
318#ifdef _RHEOLEF_HAVE_MPI
320#endif // _RHEOLEF_HAVE_MPI
321#undef _RHEOLEF_space_real
322// ======================================================================
323// parent_subgeo_owner
324// ======================================================================
325/* IMPLEMENTATION NOTES:
326
327 Consider the following situation:
328
329 P0 S0 top
330 +-----+------+--
331 |\ K0 |
332 S1 | \ |
333 | K1 \|
334 P1 +-----+
335 |
336 left |
337 +
338
339 vertex P0 is shared by edges S0 and S1 in triangle K0
340 edge S0 belongs to domain "top"
341 edge S1 belongs to domain "left"
342 assume that the mesh partition attributes
343 K0 on partition iproc=0
344 K1 on partition iproc=1
345 then, sides and vertex are also attributed to a proc
346 via a recursive descent on subgeos:
347 iproc=0 K0, S0, P0
348 iproc=1 K1, S1, P1
349 when explorating S1, its vertex P0 is already attributed to iproc=0
350
351 Then, consider:
352 space Wh (omega["left"], "P1");
353 The element edge S1 contains an **orphean** dof associated to vertex P0 :
354 it do not belong to S1 owner iproc=1 and is not referenced by any others sides
355 of "left" that belongs to its owner iproc=0
356 So, in order to interpolate in a lazy element-by-element loop, we have to
357 re-attribute this dof to S0 owner, iproc=1.
358
359 notes:
360 - this idof iproc owner reattribution required only when
361 * nproc > 1
362 * continuous elements
363 * omega is a geo_domain, i.e. some parent elements are missing
364 otherwise could return usual dis_idof owner
365
366 - all dofs are associated to a subgeo
367 - when continuous approx, this subgeo corresponds to a mesh subgeo
368 - when discontinuous, this subgeo de-multiplicated from mesh:
369 * verticies 1D, faces 3D or edges 2D : doubled when internals
370 * edges 3D and verticies 2D and 3D : de-multiplied
371 in all cases, dofs conceptually conserve the memory of the parent subgeo
372*/
373template <class T, class M>
376{
377 parent_subgeo_owner_check();
378 if (! _parent_subgeo_owner_reattribution_required) {
379 return ownership().find_owner (dis_idof);
380 }
381 return _parent_subgeo_owner.dis_at (dis_idof);
382}
383namespace details {
384struct pair_dof_parent_owners_min_op {
385 using t = std::pair<size_t,size_t>; // pair: first=parent_owner & second=dof_owner
386 t operator() (const t& a, const t& b) {
387 // if "a" or "b" parent owner match the dof owner: select it
388 // otherwise, choose the parent owner arbitrarily
389 // avoid unset=max(size_t) by tacking the min
390 if (a.first == a.second) return a;
391 if (b.first == b.second) return b;
392 size_t parent_owner = std::min (a.first, b.first); // avoid unset=max
393 return std::make_pair (parent_owner, a.second);
394 }
395};
396} // namespace details
397
398template <class T, class M>
399bool
400space_base_rep<T,M>::parent_subgeo_owner_guard() const
401{
402 if (_parent_subgeo_owner_initialized) {
403 return _parent_subgeo_owner_reattribution_required;
404 }
405 _parent_subgeo_owner_initialized = true;
406 freeze_guard();
407
408 check_macro (! get_constitution().is_hierarchical(), "parent_subgeo_owner_reattribution: invalid product-space");
409 const space_constitution_terminal<T,M>& terminal_constit = get_constitution().get_terminal();
410 size_type nproc = ownership().comm().size();
411 size_type my_proc = ownership().comm().rank();
412 const geo_basic<T,M>& omega = terminal_constit.get_geo();
413 _parent_subgeo_owner_reattribution_required
414 = terminal_constit.get_basis().is_continuous()
415 && nproc > 1
416 && omega.variant() == geo_abstract_base_rep<T>::geo_domain;
417 if (! _parent_subgeo_owner_reattribution_required) {
418 return _parent_subgeo_owner_reattribution_required;
419 }
420 // try to reduce comms: when a parent owner is the dof owner, prefer it
421 // --> see pair_dof_parent_owners_min_op
422 std::vector<size_type> dis_idx;
423 size_type unset = std::numeric_limits<size_type>::max();
424 std::pair<size_type,size_type> unset_pair (unset, unset);
425 disarray<std::pair<size_type,size_type>> dof_parent_choice (ownership(), unset_pair);
426 size_type parent_proc = my_proc;
427 std::set<size_type> ext_dis_idof;
428 for (size_type ie = 0, ne = omega.size(); ie < ne; ++ie) {
429 const geo_element& K = omega [ie];
430 get_constitution().assembly_dis_idof (omega, K, dis_idx);
431 for (size_type loc_idof = 0, loc_ndof = dis_idx.size(); loc_idof < loc_ndof; loc_idof++) {
432 size_type dis_idof = dis_idx[loc_idof];
433 size_type dof_proc = ownership().find_owner (dis_idof);
434 dof_parent_choice.dis_entry (dis_idof) = std::make_pair (parent_proc, dof_proc);
435 if (dof_proc != my_proc) {
436 ext_dis_idof.insert (dis_idof);
437 }
438 }
439 }
440 dof_parent_choice.dis_entry_assembly (details::pair_dof_parent_owners_min_op());
441 dof_parent_choice.set_dis_indexes (ext_dis_idof);
442
443 _parent_subgeo_owner.resize (ownership());
444 std::fill (_parent_subgeo_owner.begin(), _parent_subgeo_owner.end(), std::numeric_limits<size_type>::max());
445 for (size_type ie = 0, ne = omega.size(); ie < ne; ++ie) {
446 const geo_element& K = omega [ie];
447 get_constitution().assembly_dis_idof (omega, K, dis_idx);
448 for (size_type loc_idof = 0, loc_ndof = dis_idx.size(); loc_idof < loc_ndof; loc_idof++) {
449 size_type dis_idof = dis_idx[loc_idof];
450 size_type chosen_parent_proc = dof_parent_choice.dis_at(dis_idof).first;
451 // when chosen_parent_proc == unset : it means that this dof is orphean
452 // ie there are no element K in the geo_domain containing this dof and such that K belongs to the dof_owner
453 size_type reattributed_dof_proc = (chosen_parent_proc != unset) ? chosen_parent_proc : my_proc;
454 _parent_subgeo_owner.dis_entry (dis_idof) = reattributed_dof_proc;
455#ifdef TO_CLEAN
456 size_type dof_proc = ownership().find_owner (dis_idof);
457 if (dof_proc != my_proc) {
458 warning_macro ("K="<<K.name()<<K.dis_ie()<<", dis_idof="<<dis_idof<<": dof_proc="<<dof_proc<<" != my_proc="<<my_proc<<": chosen_parent_proc="<<chosen_parent_proc<<", reattributed_dof_proc="<<reattributed_dof_proc);
459 }
460#endif // TO_CLEAN
461 }
462 }
463 // choose for each dis_idof the proc owner that has the min index:
464 _parent_subgeo_owner.dis_entry_assembly (details::generic_set_min_op());
465 _parent_subgeo_owner.set_dis_indexes (ext_dis_idof);
466 return _parent_subgeo_owner_reattribution_required;
467}
468// ----------------------------------------------------------------------------
469// instanciation in library
470// ----------------------------------------------------------------------------
471#define _RHEOLEF_instanciation(T,M) \
472template class space_base_rep<T,M>; \
473template space_basic<T,M> space_basic<T,M>::real();
474
476#ifdef _RHEOLEF_HAVE_MPI
477_RHEOLEF_instanciation(Float,distributed)
478#endif // _RHEOLEF_HAVE_MPI
479#undef _RHEOLEF_instanciation
480
481} // namespace rheolef
field::size_type size_type
Definition branch.cc:430
see the Float page for the full documentation
see the communicator page for the full documentation
space_pair_type::size_type size_type
Definition space.h:164
#define trace_macro(message)
Definition dis_macros.h:111
#define warning_macro(message)
Definition dis_macros.h:53
void get_geo(istream &in, my_geo &omega)
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)
size_type nnod(const basis_basic< T > &b, const geo_size &gs, size_type map_dim)
size_type dis_nnod(const basis_basic< T > &b, const geo_size &gs, size_type map_dim)
size_type ndof(const basis_basic< T > &b, const geo_size &gs, size_type map_dim)
void dis_inod(const basis_basic< T > &b, const geo_size &gs, const geo_element &K, typename std::vector< size_type >::iterator dis_inod_tab)
This file is part of Rheolef.
#define _RHEOLEF_instanciation(T, M)
Definition space.cc:471
#define _RHEOLEF_space_real(M)
Definition space.cc:310
Expr1::memory_type M