rheolef  6.5
geo_domain_seq.cc
Go to the documentation of this file.
1 
2 #include "rheolef/geo_domain.h"
3 
4 namespace rheolef {
5 
6 template <class T>
7 void
9  const domain_indirect_rep<sequential>& indirect,
10  const geo_abstract_rep<T,sequential>& bgd_omega,
11  size_type sid_dim,
12  array<size_type,sequential>& bgd_isid2dom_dis_isid,
13  array<size_type,sequential>& dom_isid2bgd_isid,
14  array<size_type,sequential>& dom_isid2dom_ios_dis_isid,
15  size_type size_by_variant [reference_element::max_variant])
16 {
17  if (sid_dim != 0 && base::_map_dimension <= sid_dim) return;
18  communicator comm = bgd_omega.geo_element_ownership(sid_dim).comm();
19  array<size_type,sequential> bgd_isid_is_on_domain (bgd_omega.geo_element_ownership(sid_dim), 0); // logical, init to "false"
20  for (size_type ioige = 0, noige = indirect.size(); ioige < noige; ioige++) {
21  size_type ige = indirect.oige (ioige).index();
22  const geo_element& bgd_K = bgd_omega.get_geo_element (base::_map_dimension, ige);
23  for (size_type loc_isid = 0, loc_nsid = bgd_K.n_subgeo(sid_dim); loc_isid < loc_nsid; loc_isid++) {
24  size_type bgd_dis_isid = 0; // TODO: = bgd_K.subgeo (sid_dim, loc_isid);
25  switch (sid_dim) {
26  case 0: {
27  size_type bgd_dis_inod = bgd_K[loc_isid];
28  bgd_dis_isid = bgd_omega.dis_inod2dis_iv (bgd_dis_inod);
29  break;
30  }
31  case 1: bgd_dis_isid = bgd_K.edge(loc_isid); break;
32  case 2: bgd_dis_isid = bgd_K.face(loc_isid); break;
33  default: error_macro ("domain: unexpected side dimension " << sid_dim);
34  }
35  bgd_isid_is_on_domain[bgd_dis_isid] += 1;
36  }
37  }
38  size_type dom_nsid = 0;
39  for (size_type bgd_isid = 0, bgd_nsid = bgd_omega.geo_element_ownership(sid_dim).size(); bgd_isid < bgd_nsid; bgd_isid++) {
40  if (bgd_isid_is_on_domain[bgd_isid] != 0) dom_nsid++ ;
41  }
42  // 1.4 numbering dom_isid & permutation: bgd_isid <--> dom_isid
44  variant < reference_element:: last_variant_by_dimension(sid_dim); variant++) {
45  size_by_variant [variant] = 0;
46  }
47  distributor dom_isid_ownership (distributor::decide, comm, dom_nsid);
48  bgd_isid2dom_dis_isid.resize (bgd_omega.geo_element_ownership(sid_dim), std::numeric_limits<size_type>::max());
49  dom_isid2bgd_isid.resize (dom_isid_ownership, std::numeric_limits<size_type>::max());
50  for (size_type dom_isid = 0, bgd_isid = 0, bgd_nsid = bgd_omega.geo_element_ownership(sid_dim).size(); bgd_isid < bgd_nsid; bgd_isid++) {
51  if (bgd_isid_is_on_domain[bgd_isid] == 0) continue;
52  size_type dom_dis_isid = dom_isid; // sequential case !
53  bgd_isid2dom_dis_isid [bgd_isid] = dom_dis_isid;
54  dom_isid2bgd_isid [dom_isid] = bgd_isid;
55  const geo_element& bgd_S = bgd_omega.get_geo_element(sid_dim,bgd_isid);
56  size_by_variant [bgd_S.variant()]++;
57  dom_isid++;
58  }
59  size_type dom_nge = 0;
60  size_type dom_dis_nge = 0;
62  variant < reference_element:: last_variant_by_dimension(sid_dim); variant++) {
63 
64  distributor dom_igev_ownership (distributor::decide, comm, size_by_variant [variant]);
65  geo_element::parameter_type param (variant, 1);
66  base::_geo_element[variant].resize (dom_igev_ownership, param);
67  base::_gs.ownership_by_variant[variant] = dom_igev_ownership;
68  dom_nge += dom_igev_ownership.size();
69  dom_dis_nge += dom_igev_ownership.dis_size();
70  }
71  base::_gs.ownership_by_dimension[sid_dim] = distributor (dom_dis_nge, comm, dom_nge);
72 }
73 template <class T>
74 void
76  const domain_indirect_rep<sequential>& indirect,
77  const geo_abstract_rep<T,sequential>& bgd_omega,
78  array<size_type,sequential>& bgd_iv2dom_dis_iv,
79  size_type sid_dim,
80  array<size_type,sequential>& bgd_isid2dom_dis_isid,
81  array<size_type,sequential>& dom_isid2bgd_isid,
82  array<size_type,sequential>& dom_isid2dom_ios_dis_isid,
83  size_type size_by_variant [reference_element::max_variant])
84 {
85  if (sid_dim != 0 && base::_map_dimension <= sid_dim) return;
86  communicator comm = bgd_omega.geo_element_ownership(sid_dim).comm();
87  distributor dom_isid_ownership = dom_isid2bgd_isid.ownership();
88  for (size_type dom_isid = 0, dom_nsid = dom_isid_ownership.size(); dom_isid < dom_nsid; dom_isid++) {
89  size_type bgd_isid = dom_isid2bgd_isid [dom_isid];
90  size_type dom_dis_isid = dom_isid;
91  size_type dom_ios_dis_isid = dom_dis_isid;
92  const geo_element& bgd_S = bgd_omega.get_geo_element(sid_dim,bgd_isid);
93  geo_element& dom_S = get_geo_element(sid_dim,dom_isid);
94  dom_S = bgd_S;
95  dom_S.set_dis_ie (dom_dis_isid);
96  dom_S.set_ios_dis_ie (dom_ios_dis_isid);
97  }
98  if (sid_dim > 0) {
99  for (size_type dom_isid = 0, dom_nsid = dom_isid_ownership.size(); dom_isid < dom_nsid; dom_isid++) {
100  geo_element& dom_S = get_geo_element(sid_dim,dom_isid);
101  if (dom_S.dimension() == 0) continue;
102  for (size_type iloc = 0, nloc = dom_S.size(); iloc < nloc; iloc++) {
103  size_type bgd_dis_inod = dom_S[iloc];
104  size_type bgd_dis_iv = bgd_omega.dis_inod2dis_iv (bgd_dis_inod);
105  size_type dom_dis_iv = bgd_iv2dom_dis_iv.dis_at (bgd_dis_iv);
106  size_type dom_dis_inod = base::dis_iv2dis_inod (dom_dis_iv);
107  dom_S[iloc] = dom_dis_inod;
108  }
109  }
110  }
111 }
112 template <class T>
113 void
115  const domain_indirect_rep<sequential>& indirect,
116  const geo_abstract_rep<T,sequential>& bgd_omega,
117  std::map<size_type,size_type>& bgd_ie2dom_ie)
118 {
119  base::_name = bgd_omega.name() + "[" + indirect.name() + "]";
120  base::_version = 4;
121  base::_dimension = bgd_omega.dimension();
122  base::_numbering.set_degree (bgd_omega.get_piola_basis().degree());
123  base::_map_dimension = indirect.map_dimension();
124  base::_have_connectivity = 1;
125  size_type map_dim = base::_map_dimension;
126  size_type size_by_variant [reference_element::max_variant];
127  std::fill (size_by_variant, size_by_variant+reference_element::max_variant, 0);
128  array<size_type,sequential> dom_isid2dom_ios_dis_isid [4];
129  boost::array<array<size_type,sequential>,4> bgd_ige2dom_dis_ige;
130  boost::array<array<size_type,sequential>,4> dom_ige2bgd_ige;
131  domain_set_side_part1 (indirect, bgd_omega, 0,
132  bgd_ige2dom_dis_ige[0], dom_ige2bgd_ige[0],
133  dom_isid2dom_ios_dis_isid [0], size_by_variant);
134  domain_set_side_part2 (indirect, bgd_omega, bgd_ige2dom_dis_ige[0], 0,
135  bgd_ige2dom_dis_ige[0], dom_ige2bgd_ige[0],
136  dom_isid2dom_ios_dis_isid [0], size_by_variant);
138  variant < reference_element:: last_variant_by_dimension(map_dim); variant++) {
139  size_by_variant [variant] = 0;
140  }
141  for (size_type ioige = 0, noige = indirect.size(); ioige < noige; ioige++) {
142  size_type ige = indirect.oige (ioige).index();
143  bgd_ie2dom_ie [ige] = ioige;
144  const geo_element& bgd_K = bgd_omega.get_geo_element (map_dim, ige);
145  size_by_variant [bgd_K.variant()]++;
146  }
147  size_type dis_nge = 0;
148  size_type nge = 0;
150  variant < reference_element:: last_variant_by_dimension(map_dim); variant++) {
151  distributor dom_igev_ownership (distributor::decide, base::comm(), size_by_variant [variant]);
152  geo_element::parameter_type param (variant, 1);
153  base::_geo_element[variant].resize (dom_igev_ownership, param);
154  base::_gs.ownership_by_variant [variant] = dom_igev_ownership;
155  dis_nge += dom_igev_ownership.dis_size();
156  nge += dom_igev_ownership.size();
157  }
158  base::_gs.ownership_by_dimension [map_dim] = distributor (dis_nge, base::comm(), nge);
159  // => then can determine the node count (order > 1)
160  for (size_type sid_dim = 1; sid_dim < base::_map_dimension; sid_dim++) {
161  domain_set_side_part1 (indirect, bgd_omega, sid_dim,
162  bgd_ige2dom_dis_ige[sid_dim], dom_ige2bgd_ige[sid_dim],
163  dom_isid2dom_ios_dis_isid [sid_dim], size_by_variant);
164  }
165  // => then dis_iv2dis_inod works (used at step 6)
166  boost::array<size_type,reference_element::max_variant> loc_nnod_by_variant ;
167  reference_element::init_local_nnode_by_variant (base::order(), loc_nnod_by_variant);
168  size_type nnod = 0;
169  for (size_type variant = 0;
170  variant < reference_element::last_variant_by_dimension(base::map_dimension());
171  variant++) {
172  nnod += base::_gs.ownership_by_variant [variant].size() * loc_nnod_by_variant [variant];
173  }
174  distributor dom_node_ownership (distributor::decide, base::comm(), nnod);
175  base::_gs.node_ownership = dom_node_ownership;
176  size_type first_dis_ioige = indirect.ownership().first_index();
177  for (size_type ioige = 0, noige = indirect.size(); ioige < noige; ioige++) {
178  size_type dis_ioige = first_dis_ioige + ioige;
179  size_type ige = indirect.oige (ioige).index();
180  geo_element& dom_K = get_geo_element (map_dim, ioige);
181  const geo_element& bgd_K = bgd_omega.get_geo_element (map_dim, ige);
182  dom_K = bgd_K;
183  dom_K.set_dis_ie (dis_ioige);
184  size_type ini_dis_ioige = ioige;
185  dom_K.set_ios_dis_ie (ini_dis_ioige);
186  for (size_type iloc = 0, nloc = dom_K.size(); iloc < nloc; iloc++) {
187  size_type bgd_dis_inod = bgd_K[iloc];
188  size_type bgd_dis_iv = bgd_omega.dis_inod2dis_iv (bgd_dis_inod);
189  size_type dom_dis_iv = bgd_ige2dom_dis_ige[0].dis_at (bgd_dis_iv);
190  size_type dom_dis_inod = base::dis_iv2dis_inod (dom_dis_iv);
191  dom_K[iloc] = dom_dis_inod;
192  }
193  }
194  // reset also dom_ige2bgd_ige[map_dim] : used by _node[] for order > 1 geometries
195  if (base::_map_dimension > 0) {
196  dom_ige2bgd_ige [base::_map_dimension].resize(indirect.ownership());
197  for (size_type ioige = 0, noige = indirect.size(); ioige < noige; ioige++) {
198  size_type bgd_ige = indirect.oige (ioige).index();
199  dom_ige2bgd_ige [base::_map_dimension] [ioige] = bgd_ige;
200  }
201  }
202  for (size_type sid_dim = 1; sid_dim < base::_map_dimension; sid_dim++) {
203  domain_set_side_part2 (indirect, bgd_omega, bgd_ige2dom_dis_ige[0], sid_dim,
204  bgd_ige2dom_dis_ige[sid_dim], dom_ige2bgd_ige[sid_dim],
205  dom_isid2dom_ios_dis_isid [sid_dim], size_by_variant);
206  }
207  base::_node.resize (dom_node_ownership);
208  array<size_type,sequential> dom_inod2bgd_inod (dom_node_ownership, std::numeric_limits<size_type>::max());
209  size_type dom_inod = 0;
210  size_type first_bgd_inod_v = 0;
211  for (size_type dim = 0; dim <= base::_map_dimension; dim++) {
212  size_type dom_ige = 0;
215  variant++) {
216  if (loc_nnod_by_variant [variant] == 0) continue;
217  size_type first_bgd_v = bgd_omega.sizes().first_by_variant [variant].size();
218  for (size_type dom_igev = 0, dom_ngev = base::_geo_element [variant].size(); dom_igev < dom_ngev; dom_igev++, dom_ige++) {
219  const geo_element& K = base::_geo_element [variant][dom_igev];
220  size_type bgd_ige = dom_ige2bgd_ige [dim][dom_ige];
221  assert_macro (bgd_ige >= first_bgd_v, "invalid index");
222  size_type bgd_igev = bgd_ige - first_bgd_v;
223  for (size_type loc_inod = 0, loc_nnod = loc_nnod_by_variant [variant]; loc_inod < loc_nnod; loc_inod++, dom_inod++) {
224  size_type bgd_inod = first_bgd_inod_v + bgd_igev * loc_nnod_by_variant [variant] + loc_inod;
225  dom_inod2bgd_inod [dom_inod] = bgd_inod;
226  base::_node [dom_inod] = bgd_omega.node (bgd_inod);
227  }
228  }
229  first_bgd_inod_v += bgd_omega.sizes().ownership_by_variant [variant].size() * loc_nnod_by_variant [variant];
230  }
231  }
232  for (size_type dim = 1; dim < base::_map_dimension; dim++) {
233  set_element_side_index (dim);
234  }
235  for (size_type dom_iv = 0, dom_nv = base::_geo_element[reference_element::p].size(); dom_iv < dom_nv; dom_iv++) {
236  geo_element& P = base::_geo_element[reference_element::p][dom_iv];
237  size_type dom_inod = base::dis_iv2dis_inod (dom_iv);
238  P[0] = dom_inod;
239  }
243 }
244 template class geo_rep<Float,sequential>;
245 
246 } // namespace rheolef
247