Rheolef  7.1
an efficient C++ finite element environment
geo_subdivide.cc
Go to the documentation of this file.
1 //
22 // build a new mesh by subdividing each mesh edge in "nsub" sub-edges
23 //
24 // TODO:
25 // * add qTPH elts
26 // * add domains
27 // * distributed: initialize _inod2ios_dis_inod all all geo_rep<dis> data
28 #include "rheolef/geo.h"
29 #include "rheolef/space.h"
30 
31 namespace rheolef {
32 
33 // ----------------------------------------------------------------------------
34 // subdivide one element
35 // ----------------------------------------------------------------------------
36 template <class M>
37 static
38 void
39 build_edge (
40  const std::vector<size_t>& new_dis_inod,
41  size_t nsub,
42  hack_array<geo_element_hack,M>& ge,
43  size_t& ie)
44 {
46  typedef point_basic<size_type> ilat;
47  for (size_type i = 0; i < nsub; i++) {
48  size_type loc_inod0 = reference_element_e::ilat2loc_inod (nsub, ilat(i));
49  size_type loc_inod1 = reference_element_e::ilat2loc_inod (nsub, ilat(i+1));
50  ge[ie][0] = new_dis_inod[loc_inod0];
51  ge[ie][1] = new_dis_inod[loc_inod1];
52  ++ie;
53  }
54 }
55 template <class M>
56 static
57 void
58 build_triangle (
59  const std::vector<size_t>& new_dis_inod,
60  size_t nsub,
61  std::array<hack_array<geo_element_hack,M>,reference_element::max_variant>&
62  ge,
63  std::array<size_t, reference_element::max_variant>&
64  count_by_variant)
65 {
66  constexpr size_t e = reference_element::e;
67  constexpr size_t t = reference_element::t;
68  size_t& ie = count_by_variant[t];
69  size_t& e_ie = count_by_variant[e];
71  typedef point_basic<size_type> ilat;
72  for (size_type i = 0; i < nsub; i++) {
73  for (size_type j = 0; i+j < nsub; j++) {
74  // add two newly created triangles
75  size_type loc_inod00 = reference_element_t::ilat2loc_inod (nsub, ilat(i, j));
76  size_type loc_inod10 = reference_element_t::ilat2loc_inod (nsub, ilat(i+1, j));
77  size_type loc_inod01 = reference_element_t::ilat2loc_inod (nsub, ilat(i, j+1));
78  ge[t][ie][0] = new_dis_inod[loc_inod00];
79  ge[t][ie][1] = new_dis_inod[loc_inod10];
80  ge[t][ie][2] = new_dis_inod[loc_inod01];
81  ++ie;
82  if (i+j+1 >= nsub) continue;
83  size_type loc_inod11 = reference_element_t::ilat2loc_inod (nsub, ilat(i+1, j+1));
84  ge[t][ie][0] = new_dis_inod[loc_inod10];
85  ge[t][ie][1] = new_dis_inod[loc_inod11];
86  ge[t][ie][2] = new_dis_inod[loc_inod01];
87  ++ie;
88  // add also newly createe internal edges
89  ge[e][e_ie][0] = new_dis_inod[loc_inod10];
90  ge[e][e_ie][1] = new_dis_inod[loc_inod11];
91  ++e_ie;
92  ge[e][e_ie][0] = new_dis_inod[loc_inod11];
93  ge[e][e_ie][1] = new_dis_inod[loc_inod01];
94  ++e_ie;
95  ge[e][e_ie][0] = new_dis_inod[loc_inod01];
96  ge[e][e_ie][1] = new_dis_inod[loc_inod10];
97  ++e_ie;
98  }
99  }
100 }
101 template <class M>
102 static
103 void
104 build_element (
105  const reference_element& hat_K,
106  const std::vector<size_t>& new_dis_inod,
107  size_t nsub,
108  std::array<hack_array<geo_element_hack,M>,reference_element::max_variant>&
109  ge,
110  std::array<size_t, reference_element::max_variant>&
111  count_by_variant)
112 {
113  switch (hat_K.variant()) {
115  build_edge (new_dis_inod, nsub, ge[reference_element::e], count_by_variant[reference_element::e]);
116  break;
118  build_triangle (new_dis_inod, nsub, ge, count_by_variant);
119  break;
120  default:
121  error_macro ("unexpected element variant");
122  }
123 }
124 // ----------------------------------------------------------------------------
125 // mesh subdivide
126 // ----------------------------------------------------------------------------
127 template <class T, class M>
128 void
130  geo_rep<T,M>& new_omega,
131  const geo_basic<T,M>& old_omega,
132  typename geo_rep<T,M>::size_type nsub)
133 {
134  typedef typename geo_rep<T,M>::size_type size_type;
135  // ------------------------------------------------------------------------
136  // 1) store header infos in geo
137  // ------------------------------------------------------------------------
138  new_omega._version = 4;
139  new_omega._have_connectivity = true;
140  new_omega._name = old_omega.name() + "-P" + itos(nsub);
141  new_omega._dimension = old_omega.dimension();
142  new_omega._gs._map_dimension = old_omega.map_dimension();
143  new_omega._sys_coord = old_omega.coordinate_system();
144  new_omega._serial_number = 0;
145  new_omega._piola_basis = basis_basic<T>("P1");
146  // P1: subdivided mesh has order=1 for simplicity,
147  // even when original old_omega is curved, i.e. Pl, l>=1
148  // ------------------------------------------------------------------------
149  // 2) count & build nodes
150  // ------------------------------------------------------------------------
151  space_basic<T,M> old_Xh_sub (old_omega, "P"+itos(nsub));
152  new_omega._node = old_Xh_sub.get_xdofs();
153  new_omega._gs.ownership_by_variant[0]
154  = new_omega._gs.ownership_by_dimension[0]
155  = new_omega._gs.node_ownership
156  = new_omega._node.ownership();
157  // ------------------------------------------------------------------------
158  // 3) count elements
159  // ------------------------------------------------------------------------
160  size_type map_d = old_omega.map_dimension();
161  communicator comm = old_omega.comm();
162  std::array<size_type, reference_element::max_variant> size_by_variant;
163  std::fill (size_by_variant.begin(), size_by_variant.end(), 0);
164  for (size_type d = 1; d <= map_d; ++d) {
167  size_type loc_nge = pow(nsub,d);
168  size_type old_nge = old_omega.sizes().ownership_by_variant [variant].size();
169  size_by_variant[variant] = loc_nge*old_nge;
170  if (variant == reference_element::t) {
171  // count also newly locally created internal edges
172  size_type loc_nedg = 3*nsub*(nsub-1)/2;
173  size_by_variant[reference_element::e] += loc_nedg*old_nge;
174  }
175  // TODO: add internal also for q, T, P, H
176  }
177  }
178  for (size_type d = 1; d <= map_d; ++d) {
179  size_type ne = 0, dis_ne = 0;
183  new_omega._gs.ownership_by_variant [variant].resize (distributor::decide, comm, size_by_variant[variant]);
184  new_omega._geo_element [variant].resize (new_omega._gs.ownership_by_variant [variant], param);
185  ne += size_by_variant[variant];
186  }
187  new_omega._gs.ownership_by_dimension [d].resize (distributor::decide, comm, ne);
188  }
189  // ------------------------------------------------------------------------
190  // 4) compute elements
191  // ------------------------------------------------------------------------
192  std::array<size_type, reference_element::max_variant> count_by_variant;
193  std::fill (count_by_variant.begin(), count_by_variant.end(), 0);
194  std::vector<size_type> dis_inod;
195  for (size_type d = 1; d <= map_d; ++d) {
196  for (size_type ie = 0, ne = old_omega.size(d); ie < ne; ++ie) {
197  const geo_element& old_K = old_omega.get_geo_element (d, ie);
198  old_Xh_sub.dis_idof (old_K, dis_inod);
199  size_type variant = old_K.variant();
200  build_element (old_K, dis_inod, nsub, new_omega._geo_element, count_by_variant);
201  }
202  }
203  // ------------------------------------------------------------------------
204  // 5) create vertex-element (0d elements)
205  // ------------------------------------------------------------------------
207  new_omega._geo_element [reference_element::p].resize (new_omega._gs.ownership_by_dimension[0], param);
208  size_type first_iv = new_omega._gs.ownership_by_dimension[0].first_index();
209  {
210  size_type iv = 0;
211  for (auto iter = new_omega.begin(0), last = new_omega.end(0); iter != last; iter++, iv++) {
212  geo_element& P = *iter;
213  P[0] = first_iv + iv;
214  P.set_dis_ie (first_iv + iv); // TODO: P[0] & dis_ie redundant for `p'
215  P.set_ios_dis_ie (first_iv + iv);
216  }
217  }
218  // ------------------------------------------------------------------------
219  // 5b) set faces & edge in new_omega
220  // ------------------------------------------------------------------------
221  if (new_omega._gs._map_dimension > 0) {
222 
223  for (size_type side_dim = new_omega._gs._map_dimension - 1; side_dim >= 1; side_dim--) {
224 #ifdef TO_CLEAN
225  size_type nsid = new_omega._gs.ownership_by_dimension[side_dim].dis_size();
226  new_omega._gs.ownership_by_dimension [side_dim] = distributor (nsid, new_omega.comm(), nsid);
227 #endif // TO_CLEAN
228  size_type isid = 0;
229  for (typename geo_rep<T,M>::iterator iter = new_omega.begin(side_dim), last = new_omega.end(side_dim); iter != last; iter++, isid++) {
230  geo_element& S = *iter;
231  S.set_ios_dis_ie (isid);
232  S.set_dis_ie (isid);
233  }
234  }
235  }
236  // ------------------------------------------------------------------------
237  // 6) K.set_dis_ie
238  // ------------------------------------------------------------------------
239  for (size_type d = 0; d <= map_d; ++d) {
240  size_type first_dis_ie = new_omega._gs.ownership_by_dimension[d].first_index();
241  for (size_type ie = 0, ne = new_omega.size(d); ie < ne; ++ie) {
242  geo_element& K = new_omega.get_geo_element(d,ie);
243  K.set_dis_ie (first_dis_ie + ie);
244  K.set_ios_dis_ie (first_dis_ie + ie);
245  }
246  }
247 #ifdef TODO
248  // ------------------------------------------------------------------------
249  // 8) get domain, until end-of-file
250  // ------------------------------------------------------------------------
251  // TODO: build domin in new_omega from those in old_omega
252  vector<index_set> ball [4];
254  while (dom.get (ips, *this, ball)) {
255  base::_domains.push_back (dom);
256  }
257 #endif // TODO
258  // ------------------------------------------------------------------------
259  // 9) set indexes on faces and edges of elements, for P2 approx
260  // ------------------------------------------------------------------------
261  new_omega.set_element_side_index (1);
262  new_omega.set_element_side_index (2);
263  // ------------------------------------------------------------------------
264  // 10) bounding box: _xmin, _xmax
265  // ------------------------------------------------------------------------
266  new_omega.compute_bbox();
267 }
268 #define _RHEOLEF_geo_build_by_subdividing(M) \
269 template <class T> \
270 void \
271 geo_rep<T,M>::build_by_subdividing (const geo_basic<T,M>& omega, size_type nsub) \
272 { \
273  geo_build_by_subdividing (*this, omega, nsub); \
274 }
276 #ifdef _RHEOLEF_HAVE_MPI
278 #endif // _RHEOLEF_HAVE_MPI
279 #undef _RHEOLEF_geo_build_by_subdividing
280 // ----------------------------------------------------------------------------
281 // instanciation in library
282 // ----------------------------------------------------------------------------
283 #define _RHEOLEF_instanciate(T,M) \
284 template void geo_rep<Float,M>::build_by_subdividing ( \
285  const geo_basic<Float,M>& omega, size_t nsub);
287 #ifdef _RHEOLEF_HAVE_MPI
288 _RHEOLEF_instanciate(Float,distributed)
289 #endif // _RHEOLEF_HAVE_MPI
290 #undef _RHEOLEF_instanciate
291 
292 } // namespace rheolef
field::size_type size_type
Definition: branch.cc:425
see the Float page for the full documentation
see the distributor page for the full documentation
Definition: distributor.h:62
static const size_type decide
Definition: distributor.h:76
idiststream & get(idiststream &ips, const geo_rep< T, sequential > &omega, std::vector< index_set > *ball)
see the geo_element page for the full documentation
Definition: geo_element.h:102
void set_ios_dis_ie(size_type ios_dis_ie)
Definition: geo_element.h:173
variant_type variant() const
Definition: geo_element.h:161
void set_dis_ie(size_type dis_ie)
Definition: geo_element.h:172
sequential mesh representation
Definition: geo.h:778
std::vector< T, A >::size_type size_type
Definition: hack_array.h:349
static size_type ilat2loc_inod(size_type order, const point_basic< size_type > &ilat)
static size_type ilat2loc_inod(size_type order, const point_basic< size_type > &ilat)
static const variant_type e
static const variant_type max_variant
static variant_type last_variant_by_dimension(size_type dim)
static variant_type first_variant_by_dimension(size_type dim)
static const variant_type p
static const variant_type t
size_t size_type
Definition: basis_get.cc:76
#define _RHEOLEF_instanciate(T, M)
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_geo_build_by_subdividing(sequential) _RHEOLEF_instanciate(Float
std::string itos(std::string::size_type i)
itos: see the rheostream page for the full documentation
space_mult_list< T, M > pow(const space_basic< T, M > &X, size_t n)
Definition: space_mult.h:120
void geo_build_by_subdividing(geo_rep< T, M > &new_omega, const geo_basic< T, M > &old_omega, typename geo_rep< T, M >::size_type k)