Rheolef  7.1
an efficient C++ finite element environment
band.cc
Go to the documentation of this file.
1 // Banded level set routines
22 //
23 // Author: Pierre Saramito
24 //
25 #include "rheolef/band.h"
26 #include "rheolef/field_expr.h"
27 #include "rheolef/rheostream.h"
28 
29 namespace rheolef {
30 
31 // extern:
32 template<class T, class M>
33 geo_basic<T,M>
35  const field_basic<T,M>&,
36  const level_set_option&,
37  std::vector<size_t>&,
38  disarray<size_t,M>&);
39 
40 // =========================================================================
41 // part 0: utility
42 // =========================================================================
44 
45 template <class T>
46 static
47 bool
48 band_is_zero (const T& x) {
49  return fabs(x) <= band_epsilon;
50 }
51 // =========================================================================
52 // part I: subroutine: build vertex cc
53 // =========================================================================
54 template<class T, class M>
55 static
57 build_vertex_connex_component (
58  const geo_basic<T,M>& band,
59  const std::vector<geo_element::size_type>& zero_iv_list,
60  const std::vector<geo_element::size_type>& isolated_ie_list)
61 {
62  band.neighbour_guard(); // will use neighbours
63 
65  communicator comm = band.sizes().node_ownership.comm();
66  size_type map_dim = band.map_dimension();
67  size_type not_marked = std::numeric_limits<size_type>::max();
68  //
69  // 1) isolated elements are marked 0
70  // and icc-th connex component will be marked 1+icc (icc=0,1..)
71  //
72  distributor element_ownership = band.sizes().ownership_by_dimension[map_dim];
73  disarray<size_type,M> element_mark (element_ownership, not_marked);
74  for (size_type iie = 0, ine = isolated_ie_list.size(); iie < ine; iie++) {
75  element_mark [isolated_ie_list [iie]] = 0;
76  }
77  //
78  // 2) propagate a front that stop at isolated of zero marked vertices
79  //
80  // element_mark [ie] = 0; when isolated
81  // element_mark [ie] = 1+icc; when on the i-th connex component
82  //
83  size_type first_dis_ie = element_ownership.first_index();
84  size_type icc = 0;
85  bool cc_need_dis_fusion = false;
86  std::vector<std::vector<size_type> > cc_ie_list_table;
87  size_type side_dim = map_dim-1;
88  index_set ext_ie_set;
89  for (size_type ie0 = 0, ne = band.size(map_dim); ie0 < ne; ie0++) {
90  if (element_mark [ie0] != not_marked) continue;
91  element_mark [ie0] = 1+icc;
92  std::list<size_type> front;
93  front.push_back (ie0);
94  cc_ie_list_table.push_back (std::vector<size_type>());
95  std::vector<size_type>& cc_ie_list = cc_ie_list_table [icc];
96  cc_ie_list.push_back (ie0);
97  while (front.size() != 0) {
98  size_type ie = *(front.begin());
99  front.pop_front();
100  const geo_element& K = band [ie];
101  for (size_type loc_isid = 0, loc_nsid = K.n_subgeo(side_dim); loc_isid < loc_nsid; loc_isid++) {
102  size_type dis_je = band.neighbour (ie, loc_isid);
103  if (dis_je == std::numeric_limits<size_type>::max()) continue; // no neighbour
104  if (! element_ownership.is_owned (dis_je)) {
105  // K' belongs to another partition: cc is split on several parts
106  ext_ie_set += dis_je;
107  cc_need_dis_fusion = true;
108  continue;
109  }
110  // K' belongs to my partition
111  size_type je = dis_je - first_dis_ie;
112  if (element_mark [je] != not_marked) continue; // K' already studied
113  element_mark[je] = 1+icc;
114  front.push_back(je);
115  cc_ie_list.push_back(je);
116  }
117  }
118  icc++;
119  }
120  size_type ncc = cc_ie_list_table.size();
121  distributor cc_ownership (distributor::decide, comm, ncc);
122  size_type dis_ncc = cc_ownership.dis_size();
123  size_type first_dis_icc = cc_ownership.first_index();
124  //
125  // 4) here, cc fusion across partitions are neded...
126  // => fusion of vertex cc that are split by partition boundaries
127  //
128  // non-zero marks are changed from 1+icc to 1+dis_icc numbering
129  for (size_type ie = 0, ne = element_ownership.size(); ie < ne; ie++) {
130  check_macro (element_mark [ie] != not_marked, "unexpected ie="<<ie<<" not marked");
131  if (element_mark [ie] != 0) { element_mark [ie] += first_dis_icc; }
132  }
133  element_mark.set_dis_indexes (ext_ie_set);
134  //
135  // 5) compute the adjacency matrix A of cc across boundaries:
136  //
137  // send all in a pseudo-sequential matrix on process io_proc
138  size_type io_proc = odiststream::io_proc();
139  size_type my_proc = comm.rank();
140  distributor io_cc_ownership (dis_ncc, comm, (my_proc == io_proc ? dis_ncc : 0));
141  asr<Float,M> a_asr (io_cc_ownership, io_cc_ownership);
142  for (size_type ie = 0, ne = element_ownership.size(); ie < ne; ie++) {
143  if (element_mark [ie] == 0) continue;
144  const geo_element& K1 = band [ie];
145  size_type K1_dis_icc = element_mark [ie] - 1;
146  for (size_type loc_isid = 0, loc_nsid = K1.n_subgeo(side_dim); loc_isid < loc_nsid; loc_isid++) {
147  size_type dis_je = band.neighbour (ie, loc_isid);
148  if (dis_je == std::numeric_limits<size_type>::max()) continue; // no neighbour
149  size_type K2_mark = element_mark.dis_at (dis_je);
150  if (K2_mark == 0) continue;
151  size_type K2_dis_icc = K2_mark - 1;
152  a_asr.dis_entry (K1_dis_icc, K2_dis_icc) += 1;
153  a_asr.dis_entry (K2_dis_icc, K1_dis_icc) += 1;
154  }
155  }
156  a_asr.dis_entry_assembly();
157  // convert to csr:
158  csr<Float,M> a_csr (a_asr);
159  size_type nnz_a = a_csr.dis_nnz();
160  //
161  // 6) compute the closure adjacency matrix: C = lim_{n->infty} A^n
162  // A is a boolean matrix: the sequence A^n converge for a finite n
163  //
164  // convert A_csr to A as dense Eigen::Matrix on io_proc
165  //
166  // note: using "bool" instead of "int" in Eigen works with g++ but not with clang++
167  Eigen::Matrix<int,Eigen::Dynamic,Eigen::Dynamic> a = Eigen::Matrix<int,Eigen::Dynamic,Eigen::Dynamic>::Zero (a_csr.nrow(), a_csr.ncol());
168  typename csr<T,M>::const_iterator dia_ia = a_csr.begin();
169  typename csr<T,M>::const_iterator ext_ia = a_csr.ext_begin();
170  size_type first_i = a_csr.row_ownership().first_index();
171  for (size_type i = 0, n = a_csr.nrow(), q = 0; i < n; i++) {
172  size_type dis_i = first_i + i;
173  a (dis_i,dis_i) = 1; // force auto-connection, for locally empty cc (fix a bug)
174  for (typename csr<T,M>::const_data_iterator p = dia_ia[i]; p < dia_ia[i+1]; p++) {
175  size_type dis_j = first_i + (*p).first;
176  a (dis_i,dis_j) = 1;
177  }
178  }
179  Eigen::Matrix<int,Eigen::Dynamic,Eigen::Dynamic> c = a;
180  size_type nnz_c = nnz_a;
181  for (size_type k = 0, n = a_csr.nrow(); k < n; k++) {
182  Eigen::Matrix<int,Eigen::Dynamic,Eigen::Dynamic> c_tmp = c*a;
183  c = c_tmp;
184  size_type new_nnz_c = 0;
185  for (size_type i = 0; i < n; i++)
186  for (size_type j = 0; j < n; j++)
187  if (c(i,j) != 0) new_nnz_c++;
188  if (new_nnz_c == nnz_c) break; // converged
189  }
190  //
191  // 7) each row of C indicates the cc that are split by a partition bdry
192  // => re-number icc on io_proc into merged_icc
193  //
194  size_type unset = std::numeric_limits<size_type>::max();
195  disarray<size_type,M> icc2merged_dis_icc (cc_ownership, unset);
196  size_type merged_dis_ncc = 0;
197  if (my_proc == io_proc) {
198  std::vector<size_type> tmp_icc2merged_dis_icc (cc_ownership.dis_size(), unset);
199  for (size_type dis_icc = 0; dis_icc < dis_ncc; dis_icc++) {
200  if (tmp_icc2merged_dis_icc [dis_icc] != unset) continue;
201  for (size_type dis_jcc = 0; dis_jcc < dis_ncc; dis_jcc++) {
202  if (! c (dis_icc,dis_jcc)) continue;
203  // here dis_icc and dis_jcc are two parts of the same cc
204  // these two parts have been cutted by a partition boundary: now merge their indexes
205  check_macro (tmp_icc2merged_dis_icc [dis_jcc] == unset, "closure adjacency of connex component failed");
206  tmp_icc2merged_dis_icc [dis_jcc] = merged_dis_ncc;
207  }
208  merged_dis_ncc++;
209  }
210  // distribute the tmp_icc2merged_dis_icc disarray:
211  for (size_type dis_icc = 0; dis_icc < dis_ncc; dis_icc++) {
212  icc2merged_dis_icc.dis_entry (dis_icc) = tmp_icc2merged_dis_icc [dis_icc];
213  }
214  }
215  icc2merged_dis_icc.dis_entry_assembly();
216 #ifdef _RHEOLEF_HAVE_MPI
218  mpi::broadcast (comm, merged_dis_ncc, io_proc);
219  }
220 #endif // _RHEOLEF_HAVE_MPI
221  //
222  // 8) on each proc, merge parts of cc components
223  //
224  std::map<size_type,std::vector<size_type> > merged_cc_ie_list; // merged_cc[merged_icc] -> list of elt indexes
225  for (size_type icc = 0; icc < ncc; icc++) {
226  size_type merged_dis_icc = icc2merged_dis_icc [icc];
227  std::vector<size_type>& merged_cc = merged_cc_ie_list [merged_dis_icc];
228  index_set tmp_cc; // used to have a sorted index list
229  for (std::vector<size_type>::const_iterator
230  iter = cc_ie_list_table[icc].begin(),
231  last = cc_ie_list_table[icc].end(); iter != last; ++iter) {
232  size_type ie = *iter;
233  tmp_cc.insert (ie);
234  }
235  for (index_set::const_iterator
236  iter = tmp_cc.begin(),
237  last = tmp_cc.end(); iter != last; ++iter) {
238 
239  merged_cc.push_back(*iter);
240  }
241  }
242  //
243  // 9) merged icc numbering depends upon nproc: renum by sorting
244  // cc based on the smallest ios element index (nproc independant)
245  //
246  std::map<size_type,size_type> min_ios_dis_ie2dis_icc;
247  for (size_type merged_dis_icc = 0; merged_dis_icc < merged_dis_ncc; merged_dis_icc++) {
248  std::map<size_type,std::vector<size_type> >::const_iterator iter
249  = merged_cc_ie_list.find (merged_dis_icc);
250  size_type min_ios_dis_ie = std::numeric_limits<size_type>::max();
251  if (iter != merged_cc_ie_list.end()) {
252  const std::vector<size_type>& ie_list = (*iter).second;
253  for (size_type iie = 0, ine = ie_list.size(); iie < ine; iie++) {
254  size_type ie = ie_list [iie];
255  size_type ios_dis_ie = band.ige2ios_dis_ige (map_dim, ie);
256  min_ios_dis_ie = std::min (min_ios_dis_ie, ios_dis_ie);
257  }
258  }
259 #ifdef _RHEOLEF_HAVE_MPI
261  min_ios_dis_ie = mpi::all_reduce (comm, min_ios_dis_ie, mpi::minimum<size_type>());
262  }
263 #endif // _RHEOLEF_HAVE_MPI
264  min_ios_dis_ie2dis_icc [min_ios_dis_ie] = merged_dis_icc;
265  }
266  check_macro (min_ios_dis_ie2dis_icc.size() == merged_dis_ncc, "something wrong appends");
267  std::vector<size_type> ios_merged_dis_icc2merged_dis_icc (merged_dis_ncc, unset);
268  {
269  size_type ios_merged_dis_icc = 0;
270  for (std::map<size_type,size_type>::const_iterator
271  iter = min_ios_dis_ie2dis_icc.begin(),
272  last = min_ios_dis_ie2dis_icc.end(); iter != last; ++iter, ++ios_merged_dis_icc) {
273  size_type merged_dis_icc = (*iter).second;
274  ios_merged_dis_icc2merged_dis_icc [ios_merged_dis_icc] = merged_dis_icc;
275  }
276  }
277  // 10) build distributed domains from merged parts of each procs
278  //
279  std::vector<size_type> merged_cc_ie_empty_list;
280  for (size_type ios_merged_dis_icc = 0; ios_merged_dis_icc < merged_dis_ncc; ios_merged_dis_icc++) {
281 
282  size_type merged_dis_icc = ios_merged_dis_icc2merged_dis_icc [ios_merged_dis_icc];
283  check_macro (merged_dis_icc < merged_dis_ncc, "merged_dis_icc="<<merged_dis_icc<<" out of range {0:"<<merged_dis_ncc<<"[");
284  std::string merged_cc_name = "cc" + itos (ios_merged_dis_icc);
285  const std::vector<size_type>* ptr_merged_cc_ie_list;
286  std::map<size_type,std::vector<size_type> >::const_iterator iter
287  = merged_cc_ie_list.find (merged_dis_icc);
288  if (iter != merged_cc_ie_list.end()) {
289  ptr_merged_cc_ie_list = & ((*iter).second);
290  } else {
291  ptr_merged_cc_ie_list = & merged_cc_ie_empty_list;
292  }
293  domain_indirect_basic<M> merged_cc_dom (band, merged_cc_name, map_dim, comm, *ptr_merged_cc_ie_list);
294  band.insert_domain_indirect (merged_cc_dom);
295  }
296  return merged_dis_ncc;
297 }
298 // =========================================================================
299 // part II: the main band constructor
300 // =========================================================================
301 template<class T, class M>
303  const field_basic<T,M>& phi_h,
304  const level_set_option& opt)
305  : _gamma(),
306  _band(),
307  _sid_ie2bnd_ie(),
308  _ncc(0)
309 {
310  using namespace details;
311  band_epsilon = opt.epsilon; // set global variable
312  //
313  // 1) build the level set
314  //
315  std::vector<size_type> bnd_dom_ie_list;
316  _gamma = level_set_internal (phi_h, opt, bnd_dom_ie_list, _sid_ie2bnd_ie);
317  //
318  // 2) build the band
319  //
320  // 2.1) build the band domain in lambda
321  const geo_basic<T,M>& lambda = phi_h.get_geo();
322  check_macro(lambda.order() == 1, "Only order=1 band supported");
323  size_type map_dim = lambda.map_dimension();
324  communicator comm = lambda.geo_element_ownership(map_dim).comm();
325  domain_indirect_basic<M> band_indirect (lambda, "band", map_dim, comm, bnd_dom_ie_list);
326  geo_basic<T,M> band_domain = geo_basic<T,M> (band_indirect, lambda);
327  // 2.2) reduce phi_h to the band domain & get its geo_domain band
328  field_basic<T,M> phi_h_band (phi_h [band_domain]);
329  _band = phi_h_band.get_geo();
330  //
331  // 3) create the "zero" containing vertices xi such that phi(xi) = 0
332  // TODO: Pk and k>=2: zero is a domain of idofs: Bh.block(isolated); Bh.unblock (zeros)
333  //
334  check_macro(phi_h_band.get_approx() == "P1", "Only P1 level set function supported");
335  std::vector<size_type> zero_iv_list;
336  for (size_type iv = 0, nv = _band.n_vertex(); iv < nv; iv++) {
337  size_type idof = iv; // assume phi_h is P1
338  if (band_is_zero (phi_h_band.dof(idof))) zero_iv_list.push_back (iv);
339  }
340  domain_indirect_basic<M> zero_dom (_band, "zero", 0, comm, zero_iv_list);
341  _band.insert_domain_indirect (zero_dom);
342  //
343  // 4) create the "isolated" domain containing vertices xi such that:
344  // (i) phi(xi) != 0
345  // (ii) for all element K that contains xi,
346  // we have phi(xj) = 0 for all vertex xj of K and xj != xi
347  //
348  // 4.1) mark isolated vertices
349  phi_h_band.dis_dof_update();
350  distributor element_ownership = _band.sizes().ownership_by_dimension[map_dim];
351  size_type not_marked = std::numeric_limits<size_type>::max();
352  disarray<size_type,M> is_isolated (element_ownership, 1); // 1=true, for bool
353  std::vector<size_type> isolated_ie_list;
354  size_type ie = 0;
355  for (typename geo_basic<T,M>::const_iterator
356  iter = _band.begin(map_dim),
357  last = _band.end(map_dim); iter != last; ++iter, ++ie) {
358  const geo_element& K = *iter;
359  size_type count_non_zero = 0;
360  // TODO Pk: K is isolated when all non-zero are located inside one face exactly
361  // but on this face some can be non-zero (at least 2 ?)
362  // here: only for P1: only one dof is non-zero
363  size_type loc_inod_isolated = std::numeric_limits<size_type>::max();
364  for (size_type loc_inod = 0, loc_nnod = K.size(); loc_inod < loc_nnod; loc_inod++) {
365  size_type dis_inod = K [loc_inod];
366  if (! band_is_zero (phi_h_band.dis_dof (dis_inod))) {
367  count_non_zero++;
368  loc_inod_isolated = loc_inod;
369  }
370  }
371  if (count_non_zero == 1) {
372  is_isolated[ie] = 1;
373  isolated_ie_list.push_back (ie);
374  }
375  }
376  domain_indirect_basic<M> isolated_dom (_band, "isolated", map_dim, comm, isolated_ie_list);
377  _band.insert_domain_indirect (isolated_dom);
378  //
379  // 5) build cc
380  //
381  _ncc = build_vertex_connex_component (_band, zero_iv_list, isolated_ie_list);
382 }
383 // ----------------------------------------------------------------------------
384 // instanciation in library
385 // ----------------------------------------------------------------------------
386 #define _RHEOLEF_instanciation(T,M) \
387 template class band_basic<T,M>;
388 
389 _RHEOLEF_instanciation(Float,sequential)
390 #ifdef _RHEOLEF_HAVE_MPI
391 _RHEOLEF_instanciation(Float,distributed)
392 #endif // _RHEOLEF_HAVE_MPI
393 
394 } // namespace
field::size_type size_type
Definition: branch.cc:425
see the Float page for the full documentation
see the band page for the full documentation
geo_basic< T, M > _gamma
Definition: band.h:116
geo_basic< T, M > _band
Definition: band.h:117
disarray< size_type, M > _sid_ie2bnd_ie
Definition: band.h:118
size_type _ncc
Definition: band.h:119
geo_basic< T, M >::size_type size_type
Definition: band.h:98
see the distributor page for the full documentation
Definition: distributor.h:62
static const size_type decide
Definition: distributor.h:76
const T & dis_dof(size_type dis_idof) const
Definition: field.cc:376
std::string get_approx() const
Definition: field.h:303
const geo_type & get_geo() const
Definition: field.h:301
void dis_dof_update() const
Definition: field.cc:410
T & dof(size_type idof)
Definition: field.h:658
see the geo_element page for the full documentation
Definition: geo_element.h:102
size_type size() const
Definition: geo_element.h:168
reference_element::size_type size_type
Definition: geo_element.h:125
static size_type io_proc()
Definition: diststream.cc:78
size_t size_type
Definition: basis_get.cc:76
static Float band_epsilon
Definition: band.cc:43
Expr1::float_type T
Definition: field_expr.h:261
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
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.
_RHEOLEF_instanciation(Float) _RHEOLEF_instanciation_evaluate(Float
std::string itos(std::string::size_type i)
itos: see the rheostream page for the full documentation
geo_basic< T, M > level_set_internal(const field_basic< T, M > &, const level_set_option &, std::vector< size_t > &, disarray< size_t, M > &)
Definition: level_set.cc:415
Definition: sphere.icc:25
static const bool value
Definition: distributed.h:37
Float epsilon