25 #include "rheolef/band.h"
26 #include "rheolef/field_expr.h"
27 #include "rheolef/rheostream.h"
32 template<
class T,
class M>
35 const field_basic<T,M>&,
36 const level_set_option&,
48 band_is_zero (
const T& x) {
54 template<
class T,
class M>
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)
62 band.neighbour_guard();
65 communicator comm =
band.sizes().node_ownership.comm();
67 size_type not_marked = std::numeric_limits<size_type>::max();
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;
83 size_type first_dis_ie = element_ownership.first_index();
85 bool cc_need_dis_fusion =
false;
86 std::vector<std::vector<size_type> > cc_ie_list_table;
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) {
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++) {
103 if (dis_je == std::numeric_limits<size_type>::max())
continue;
104 if (! element_ownership.is_owned (dis_je)) {
106 ext_ie_set += dis_je;
107 cc_need_dis_fusion =
true;
112 if (element_mark [je] != not_marked)
continue;
113 element_mark[je] = 1+icc;
115 cc_ie_list.push_back(je);
122 size_type dis_ncc = cc_ownership.dis_size();
123 size_type first_dis_icc = cc_ownership.first_index();
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; }
133 element_mark.set_dis_indexes (ext_ie_set);
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++) {
148 if (dis_je == std::numeric_limits<size_type>::max())
continue;
149 size_type K2_mark = element_mark.dis_at (dis_je);
150 if (K2_mark == 0)
continue;
152 a_asr.dis_entry (K1_dis_icc, K2_dis_icc) += 1;
153 a_asr.dis_entry (K2_dis_icc, K1_dis_icc) += 1;
156 a_asr.dis_entry_assembly();
158 csr<Float,M> a_csr (a_asr);
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++) {
174 for (
typename csr<T,M>::const_data_iterator
p = dia_ia[i];
p < dia_ia[i+1];
p++) {
179 Eigen::Matrix<int,Eigen::Dynamic,Eigen::Dynamic>
c =
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;
187 if (
c(i,j) != 0) new_nnz_c++;
188 if (new_nnz_c == nnz_c)
break;
194 size_type unset = std::numeric_limits<size_type>::max();
195 disarray<size_type,M> icc2merged_dis_icc (cc_ownership, unset);
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;
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;
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];
215 icc2merged_dis_icc.dis_entry_assembly();
216 #ifdef _RHEOLEF_HAVE_MPI
218 mpi::broadcast (comm, merged_dis_ncc, io_proc);
224 std::map<size_type,std::vector<size_type> > merged_cc_ie_list;
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];
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) {
235 for (index_set::const_iterator
236 iter = tmp_cc.begin(),
237 last = tmp_cc.end(); iter != last; ++iter) {
239 merged_cc.push_back(*iter);
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++) {
256 min_ios_dis_ie = std::min (min_ios_dis_ie, ios_dis_ie);
259 #ifdef _RHEOLEF_HAVE_MPI
261 min_ios_dis_ie = mpi::all_reduce (comm, min_ios_dis_ie, mpi::minimum<size_type>());
264 min_ios_dis_ie2dis_icc [min_ios_dis_ie] = merged_dis_icc;
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);
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;
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++) {
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);
291 ptr_merged_cc_ie_list = & merged_cc_ie_empty_list;
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);
296 return merged_dis_ncc;
301 template<
class T,
class M>
310 using namespace details;
315 std::vector<size_type> bnd_dom_ie_list;
324 communicator comm =
lambda.geo_element_ownership(
map_dim).comm();
335 std::vector<size_type> zero_iv_list;
338 if (band_is_zero (phi_h_band.
dof(idof))) zero_iv_list.push_back (iv);
341 _band.insert_domain_indirect (zero_dom);
351 size_type not_marked = std::numeric_limits<size_type>::max();
353 std::vector<size_type> isolated_ie_list;
357 last =
_band.end(
map_dim); iter != last; ++iter, ++ie) {
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++) {
368 loc_inod_isolated = loc_inod;
371 if (count_non_zero == 1) {
373 isolated_ie_list.push_back (ie);
377 _band.insert_domain_indirect (isolated_dom);
381 _ncc = build_vertex_connex_component (
_band, zero_iv_list, isolated_ie_list);
386 #define _RHEOLEF_instanciation(T,M) \
387 template class band_basic<T,M>;
390 #ifdef _RHEOLEF_HAVE_MPI
field::size_type size_type
see the Float page for the full documentation
see the band page for the full documentation
disarray< size_type, M > _sid_ie2bnd_ie
geo_basic< T, M >::size_type size_type
see the distributor page for the full documentation
static const size_type decide
const T & dis_dof(size_type dis_idof) const
std::string get_approx() const
const geo_type & get_geo() const
void dis_dof_update() const
see the geo_element page for the full documentation
reference_element::size_type size_type
static size_type io_proc()
static Float band_epsilon
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.
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 > &)
_RHEOLEF_instanciation(Float, sequential, std::allocator< Float >) _RHEOLEF_instanciation(Float