Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
geo_subdivide.cc
Go to the documentation of this file.
1
21//
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
31namespace rheolef {
32
33// ----------------------------------------------------------------------------
34// subdivide one element
35// ----------------------------------------------------------------------------
36template <class M>
37static
38void
39build_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}
55template <class M>
56static
57void
58build_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}
101template <class M>
102static
103void
104build_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// ----------------------------------------------------------------------------
127template <class T, class M>
128void
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" + std::to_string(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"+std::to_string(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) {
166 variant < reference_element:: last_variant_by_dimension(d); variant++) {
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;
181 variant < reference_element:: last_variant_by_dimension(d); variant++) {
182 geo_element::parameter_type param (variant, 1);
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];
253 domain_indirect_basic<sequential> dom;
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) \
269template <class T> \
270void \
271geo_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) \
284template void geo_rep<Float,M>::build_by_subdividing ( \
285 const geo_basic<Float,M>& omega, size_t nsub);
287#ifdef _RHEOLEF_HAVE_MPI
289#endif // _RHEOLEF_HAVE_MPI
290#undef _RHEOLEF_instanciate
291
292} // namespace rheolef
field::size_type size_type
Definition branch.cc:430
see the Float page for the full documentation
see the distributor page for the full documentation
Definition distributor.h:69
static const size_type decide
Definition distributor.h:83
generic mesh with rerefence counting
Definition geo.h:1089
see the geo_element page for the full documentation
void set_ios_dis_ie(size_type ios_dis_ie)
variant_type variant() const
void set_dis_ie(size_type dis_ie)
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
the finite element space
Definition space.h:382
rheolef::space_base_rep< T, M > t
#define _RHEOLEF_instanciate(T)
Definition csr_seq.cc:491
#define error_macro(message)
Definition dis_macros.h:49
#define _RHEOLEF_geo_build_by_subdividing(M)
This file is part of Rheolef.
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)
space_mult_list< T, M > pow(const space_basic< T, M > &X, size_t n)
Definition space_mult.h:120