Rheolef  7.1
an efficient C++ finite element environment
geo_neighbour.cc
Go to the documentation of this file.
1 // Banded level set routines
22 //
23 // Author: Pierre Saramito
24 //
25 #include "rheolef/geo.h"
26 
27 namespace 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
36 template<class T, class M>
37 static
38 void
39 add_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)
45 template<class T>
46 static
47 void
48 add_ball_externals (
49  const geo_base_rep<T,sequential>& omega,
50  const disarray<index_set,sequential>& ball)
51 {
52 }
53 #ifdef _RHEOLEF_HAVE_MPI
54 template<class T>
55 static
56 void
57 add_ball_externals (
58  const geo_base_rep<T,distributed>& omega,
59  const disarray<index_set,distributed>& ball)
60 {
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 // ----------------------------------------------------------------------------
81 template<class T, class M>
82 void
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
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
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
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
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)
232 template<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) \
269 template void geo_base_rep<T,M>::init_neighbour() const; \
270 template geo_base_rep<T,M>::size_type \
271  geo_base_rep<T,M>::neighbour (size_type, size_type) const;
272 
273 _RHEOLEF_instanciation(Float,sequential)
274 #ifdef _RHEOLEF_HAVE_MPI
276 #endif // _RHEOLEF_HAVE_MPI
277 
278 } // namespace
rheolef::reference_element::last_variant_by_dimension
static variant_type last_variant_by_dimension(size_type dim)
Definition: reference_element.h:150
rheolef::geo_base_rep< T, distributed >::size_type
base::size_type size_type
Definition: geo.h:533
rheolef::geo_base_rep::const_iterator
base::const_iterator const_iterator
Definition: geo.h:537
check_macro
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
rheolef::geo_element::master
size_type master(bool i) const
Definition: geo_element.h:165
rheolef::index_set::inplace_intersection
void inplace_intersection(const index_set &b)
Definition: index_set_body.icc:126
rheolef::_RHEOLEF_instanciation
_RHEOLEF_instanciation(Float, sequential, std::allocator< Float >) _RHEOLEF_instanciation(Float
rheolef::distributor
see the distributor page for the full documentation
Definition: distributor.h:62
rheolef::geo_element::n_subgeo
size_type n_subgeo(size_type subgeo_dim) const
Definition: geo_element.h:212
rheolef::geo_element::size
size_type size() const
Definition: geo_element.h:168
rheolef::geo_element
see the geo_element page for the full documentation
Definition: geo_element.h:102
rheolef::index_set
Definition: index_set_header.icc:58
rheolef::size_type
size_t size_type
Definition: basis_get.cc:76
rheolef::geo_element::get_side_orientation
orientation_type get_side_orientation(const geo_element &S) const
Definition: geo_element.cc:240
rheolef::reference_element::first_variant_by_dimension
static variant_type first_variant_by_dimension(size_type dim)
Definition: reference_element.h:148
rheolef::geo_element::dis_ie
size_type dis_ie() const
Definition: geo_element.h:163
mkgeo_ball.variant
int variant
Definition: mkgeo_ball.sh:149
vertex
const point vertex[n_vertex]
Definition: edge.icc:67
rheolef::geo_element::n_node
size_type n_node() const
Definition: geo_element.h:170
rheolef::geo_base_rep::init_neighbour
void init_neighbour() const
Definition: geo_neighbour.cc:83
rheolef::geo_element::orientation_type
geo_element_indirect::orientation_type orientation_type
Definition: geo_element.h:134
rheolef
This file is part of Rheolef.
Definition: compiler_eigen.h:37
rheolef::geo_element::set_master
void set_master(bool i, size_type dis_ie) const
Definition: geo_element.h:174
error_macro
#define error_macro(message)
Definition: dis_macros.h:49
rheolef::space_numbering::dis_inod
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)
Definition: space_numbering.cc:177
Float
see the Float page for the full documentation
rheolef::disarray
see the disarray page for the full documentation
Definition: disarray.h:459
size_type
field::size_type size_type
Definition: branch.cc:425
rheolef::distributor::is_owned
bool is_owned(size_type dis_i, size_type iproc) const
true when dis_i in [first_index(iproc):last_index(iproc)[
Definition: distributor.h:213
rheolef::distributor::first_index
size_type first_index(size_type iproc) const
global index range and local size owned by ip-th process
Definition: distributor.h:151
mkgeo_ball.map_dim
map_dim
Definition: mkgeo_ball.sh:337
rheolef::geo_element::face
size_type face(size_type i) const
Definition: geo_element.h:210
rheolef::distributed
distributed
Definition: asr.cc:228
rheolef::geo_element::edge
size_type edge(size_type i) const
Definition: geo_element.h:209
rheolef::geo_base_rep::neighbour
size_type neighbour(size_type ie, size_type loc_isid) const
Definition: geo_neighbour.cc:234