Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
geo_neighbour.cc
Go to the documentation of this file.
1
21// Banded level set routines
22//
23// Author: Pierre Saramito
24//
25#include "rheolef/geo.h"
26
27namespace rheolef {
28
29// ----------------------------------------------------------------------------
30// 1) utilities for neighbour initialization
31// ----------------------------------------------------------------------------
32// Implementation note: as member function in geo_base_rep<T,M> could be simpler
33// but function member of template class cannot be specialized
34// and the sequential case do not define scatter_map_type and such
35// so use a specialization of a non-member function
36template<class T, class M>
37static
38void
39add_ball_externals (
40 const geo_base_rep<T,M>& omega,
41 const disarray<index_set,M>& ball)
42{
43}
44// explicit M=sequential specialization (otherwise: linker problems)
45template<class T>
46static
47void
48add_ball_externals (
49 const geo_base_rep<T,sequential>& omega,
50 const disarray<index_set,sequential>& ball)
51{
52}
53#ifdef _RHEOLEF_HAVE_MPI
54template<class T>
55static
56void
57add_ball_externals (
58 const geo_base_rep<T,distributed>& omega,
59 const disarray<index_set,distributed>& ball)
60{
61 typedef typename geo_base_rep<T,distributed>::size_type size_type;
62 index_set ext_iv_set;
63 // faudrait attraper _geo_element[0].get_dis_map_entries() : suit les vertices au lieu nodes
64 const hack_array<geo_element_hack,distributed>& vertex = omega._geo_element[0];
65 typedef typename hack_array<geo_element_hack,distributed>::scatter_map_type map_t;
66 const map_t& map_vertex = vertex.get_dis_map_entries();
67 size_type dis_nv = omega.dis_size (0);
68 for (typename map_t::const_iterator
69 iter = map_vertex.begin(),
70 last = map_vertex.end(); iter != last; ++iter) {
71 size_type dis_iv = (*iter).first;
72 check_macro (dis_iv < dis_nv, "dis_iv="<<dis_iv<<" out of range [0:"<<dis_nv<<"[");
73 ext_iv_set += dis_iv;
74 }
75 ball.set_dis_indexes (ext_iv_set);
76}
77#endif // _RHEOLEF_HAVE_MPI
78// ----------------------------------------------------------------------------
79// 2) (lazy) initialization of neighbour informations
80// ----------------------------------------------------------------------------
81template<class T, class M>
82void
84{
85 if (is_broken()) return; // skip broken domain e.g. "sides" or "internal_sides": no neighbours
86 size_type map_dim = map_dimension();
87 if (map_dim <= 0) return; // a set of points have no neighbours
88 size_type sid_dim = map_dim-1;
89
90 // 0.0) clean-up neighbour information (could come from a copy of elements, from domain, e.g. band)
91 for (const_iterator iter_S = begin(sid_dim), last_S = end(sid_dim); iter_S != last_S; ++iter_S) {
92 const geo_element& S = *iter_S;
93 S.set_master (0, std::numeric_limits<size_type>::max());
94 S.set_master (1, std::numeric_limits<size_type>::max());
95 }
96 // 0.1) manage all external sides S that appears on some local elements K
97 distributor is_ownership = _gs.ownership_by_dimension[map_dim-1];
98 std::array<index_set, reference_element::max_variant> ext_isv_set; // external sides, by variants
99 index_set ext_is_set; // external sides, all variants merged
100 for (const_iterator iter_K = begin(map_dim), last_K = end(map_dim); iter_K != last_K; ++iter_K) {
101 const geo_element& K = *iter_K;
102 for (size_type loc_is = 0, loc_ns = K.n_subgeo(sid_dim); loc_is < loc_ns; ++loc_is) {
103 size_type dis_is = (map_dim == 1) ? K[loc_is] : ((map_dim == 2) ? K.edge(loc_is) : K.face(loc_is));
104 if (is_ownership.is_owned (dis_is)) continue;
105 // have an external side S: convert its dis_is to dis_isv and get its variant
106 size_type variant;
107 size_type dis_isv = _gs.dis_ige2dis_igev_by_dimension (sid_dim, dis_is, variant);
108 ext_isv_set [variant] += dis_isv;
109 ext_is_set += dis_is;
110 }
111 }
112 // for external side S that appears in a local element K may be available
113 // => enlarge the external set of _geo_element for sid_dim
114 // but ext_is_set may be classified by variants and all dis_is may be
115 // converted to dis_isv
117 variant < reference_element:: last_variant_by_dimension(sid_dim); variant++) {
118 _geo_element [variant].append_dis_indexes (ext_isv_set [variant]);
119 }
120 // 1) build ball(xi) := { K; xi is a vertex of K }
121 index_set emptyset;
122 distributor vertex_ownership = sizes().ownership_by_dimension[0];
123 disarray<index_set,M> ball (vertex_ownership, emptyset);
124 for (const_iterator iter = begin(map_dim), last = end(map_dim); iter != last; ++iter) {
125 const geo_element& K = *iter;
126 index_set dis_ie_set;
127 dis_ie_set += K.dis_ie();
128 for (size_type loc_inod = 0, loc_nnod = K.size(); loc_inod < loc_nnod; loc_inod++) {
129 size_type dis_inod = K [loc_inod];
130 size_type dis_iv = dis_inod2dis_iv (dis_inod);
131 ball.dis_entry (dis_iv) += dis_ie_set; // union with {dis_ie}
132 }
133 }
134 ball.dis_entry_assembly();
135 // 2) set ball available at partition boundaries
136 add_ball_externals (*this, ball);
137 // 3) set master elements on all sides
138 distributor ie_ownership = _gs.ownership_by_dimension[map_dim];
139 distributor inod_ownership = _gs.node_ownership;
140 size_type first_dis_ie = ie_ownership.first_index();
141 std::array<index_set, reference_element::max_variant> ext_igev_set; // external elts, by variants
142 index_set ext_ie_set; // external elements, all variants merged
143 index_set ext_inod_set;
144 disarray<size_type, M> side_marked (is_ownership, 0); // boolean: 0,1
145 side_marked.set_dis_indexes (ext_is_set);
146 for (const_iterator iter_K = begin(map_dim), last_K = end(map_dim); iter_K != last_K; ++iter_K) {
147 const geo_element& K = *iter_K;
148 for (size_type loc_is = 0, loc_ns = K.n_subgeo(sid_dim); loc_is < loc_ns; ++loc_is) {
149 size_type dis_is = (sid_dim == 0) ? K[loc_is] : ((sid_dim == 1) ? K.edge(loc_is) : K.face(loc_is));
150 if (side_marked.dis_at (dis_is)) continue;
151 side_marked.dis_entry (dis_is) = 1;
152 const geo_element& S = dis_get_geo_element (sid_dim, dis_is);
153 size_type dis_iv0 = dis_inod2dis_iv (S[0]);
154 index_set ball_intersect = ball.dis_at (dis_iv0);
155 for (size_type loc_iv = 1, loc_nv = S.size(); loc_iv < loc_nv; loc_iv++) {
156 size_type dis_iv = dis_inod2dis_iv (S[loc_iv]);
157 ball_intersect.inplace_intersection (ball.dis_at (dis_iv));
158 }
159 switch (ball_intersect.size()) {
160 case 1: { // one master K: S is on the boundary
161 index_set::const_iterator iter = ball_intersect.begin();
162 size_type dis_ie = *iter;
163 S.set_master (0, dis_ie);
164 S.set_master (1, std::numeric_limits<size_type>::max());
165 break;
166 }
167 case 2: { // two masters K1, K2: S is an internal side
168 index_set::const_iterator iter = ball_intersect.begin();
169 size_type dis_ie1 = *iter++;
170 size_type dis_ie2 = *iter;
171 // one of K1 and K2 is available in distributed environment
172 // the other is not a priori available
173 if (! ie_ownership.is_owned (dis_ie1)) std::swap (dis_ie1, dis_ie2);
174 if (! ie_ownership.is_owned (dis_ie2)) {
175 // convert dis_ie2 to dis_igev2 and get its variant2:
176 size_type variant2;
177 size_type dis_igev2 = _gs.dis_ige2dis_igev_by_dimension (map_dim, dis_ie2, variant2);
178 ext_igev_set [variant2] += dis_igev2;
179 ext_ie_set += dis_ie2;
180 }
181 check_macro (ie_ownership.is_owned (dis_ie1), "unexpected external K1");
182 // let K1 indicated by the normal on S given by S orient
183 // let K2, the other one
184 size_type ie1 = dis_ie1 - first_dis_ie;
185 const geo_element& K1 = get_geo_element (map_dim, ie1);
187 if (orient < 0) std::swap (dis_ie1, dis_ie2);
188 S.set_master (0, dis_ie1);
189 S.set_master (1, dis_ie2);
190 break;
191 }
192 default: {
193 error_macro ("neighbour: side #"<< ball_intersect.size()<<" connectivity problem");
194 }
195 }
196 }
197 }
198 // side values have been changed:
199 // propagate these modifs for sides at the partition boundaries
201 variant < reference_element:: last_variant_by_dimension(sid_dim); variant++) {
202 _geo_element [variant].update_dis_entries();
203 }
204 // for all S in the partition boundary, both K1 & K2 may be available
205 // => enlarge the external set of _geo_element for map_dim
206 // but ext_ie_set may be classified by variants and all dis_ie may be
207 // converted to dis_igev
209 variant < reference_element:: last_variant_by_dimension(map_dim); variant++) {
210 _geo_element [variant].append_dis_indexes (ext_igev_set [variant]);
211 }
212 // manage also external nodes of external elements K:
213 for (index_set::const_iterator iter = ext_ie_set.begin(), last = ext_ie_set.end(); iter != last; ++iter) {
214 size_type dis_ie = *iter;
215 const geo_element& K = dis_get_geo_element (map_dim, dis_ie);
216 for (size_type loc_inod = 0, loc_nnod = K.n_node(); loc_inod < loc_nnod; loc_inod++) {
217 size_type dis_inod = K[loc_inod];
218 if (! inod_ownership.is_owned (dis_inod)) {
219 ext_inod_set += dis_inod;
220 }
221 }
222 }
223 _node.append_dis_indexes (ext_inod_set);
224}
225// ----------------------------------------------------------------------------
226// 3) compute the neighbour K2 of an element K1 across a side S
227// ----------------------------------------------------------------------------
228// returns the dis_ie2 global index associated to K2,
229// the neighbour element of K1, with index ie1,
230// across the side S of local index loc_isid in K1
231// when S is on the boundary, returns size_t(-1)
232template<class T, class M>
235{
236 check_macro (_have_neighbour, "neighbour_guard() may be called globally before local neighbour()");
237 // 1) get the i-th side S in K1
238 size_type map_dim = map_dimension();
239 check_macro (ie1 < size(map_dim), "neighbour: element index ie=" << ie1
240 << " is out of range [0:" << size(map_dim) << "[");
241 const geo_element& K1 = get_geo_element (map_dim, ie1);
242 size_type dis_isid;
243 switch (map_dim) {
244 case 3: dis_isid = K1.face (loc_isid); break;
245 case 2: dis_isid = K1.edge (loc_isid); break;
246 case 1: {
247 size_type dis_inod = K1[loc_isid];
248 size_type dis_iv = dis_inod2dis_iv (dis_inod);
249 dis_isid = dis_iv;
250 break;
251 }
252 case 0:
253 default: return std::numeric_limits<size_type>::max();
254 }
255 size_type sid_dim = map_dim-1;
256 const geo_element& S = dis_get_geo_element (sid_dim, dis_isid);
257
258 // 2) get the neighbour element K2 to K1 across S
259 size_type dis_ie2 = S.master (0);
260 size_type first_dis_ie = _gs.ownership_by_dimension[map_dim].first_index();
261 size_type dis_ie1 = first_dis_ie + ie1;
262 if (dis_ie2 != dis_ie1) return dis_ie2;
263 return S.master(1);
264}
265// ----------------------------------------------------------------------------
266// instanciation in library
267// ----------------------------------------------------------------------------
268#define _RHEOLEF_instanciation(T,M) \
269template void geo_base_rep<T,M>::init_neighbour() const; \
270template geo_base_rep<T,M>::size_type \
271 geo_base_rep<T,M>::neighbour (size_type, size_type) const;
272
274#ifdef _RHEOLEF_HAVE_MPI
275_RHEOLEF_instanciation(Float,distributed)
276#endif // _RHEOLEF_HAVE_MPI
277
278} // namespace
#define _RHEOLEF_instanciation(T, M, A)
Definition asr.cc:223
field::size_type size_type
Definition branch.cc:430
see the Float page for the full documentation
see the disarray page for the full documentation
Definition disarray.h:497
see the distributor page for the full documentation
Definition distributor.h:69
size_type first_index(size_type iproc) const
global index range and local size owned by ip-th process
bool is_owned(size_type dis_i, size_type iproc) const
true when dis_i in [first_index(iproc):last_index(iproc)[
base::size_type size_type
Definition geo.h:533
void init_neighbour() const
size_type neighbour(size_type ie, size_type loc_isid) const
base::const_iterator const_iterator
Definition geo.h:537
see the geo_element page for the full documentation
geo_element_indirect::orientation_type orientation_type
size_type edge(size_type i) const
size_type face(size_type i) const
size_type n_node() const
size_type size() const
orientation_type get_side_orientation(const geo_element &S) const
size_type master(bool i) const
size_type n_subgeo(size_type subgeo_dim) const
size_type dis_ie() const
void set_master(bool i, size_type dis_ie) const
void inplace_intersection(const index_set &b)
static variant_type last_variant_by_dimension(size_type dim)
static variant_type first_variant_by_dimension(size_type dim)
#define error_macro(message)
Definition dis_macros.h:49
const point vertex[n_vertex]
Definition edge.icc:67
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.