Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
mkgeo_ball_gmsh_fix.cc
Go to the documentation of this file.
1
21// fix the blending function for curved elements in gmsh:
22// face boundary interior points are incorrects
23// TODO:
24// - fix_full : filled circle and sphere meshes
25// - fix_s : when meshes by quadrangles
26// - fix_s : when input geometry is an ellipse
27// - support mpirun : node are distributed
28//
29#include "rheolef.h"
30#include "rheolef/geo_element_curved_ball.h"
31#include <cstring>
32using namespace rheolef;
33point project_on_unit_sphere (const point& x) { return x/norm(x); }
35{
37 typedef point_basic<size_type> ilat;
38 size_t order = gamma.order();
39 disarray<point,sequential> node = gamma.get_nodes();
40 // 1) put all nodes exactly on the unit sphere
41 for (size_type iv = 0, nv = gamma.n_vertex(); iv < nv; iv++) {
42 node[iv] = project_on_unit_sphere (node[iv]);
43 }
44 // 2) fix nodes along edges:
45 std::vector<size_type> dis_inod;
46 for (size_type iedg = 0, nedg = gamma.size(1); iedg < nedg; iedg++) {
47 const geo_element& E = gamma.get_geo_element (1, iedg);
48 gamma.dis_inod (E, dis_inod);
49 const point& x0 = node[dis_inod[0]];
50 const point& x1 = node[dis_inod[1]];
51 for (size_type i = 1; i < order; i++) {
52 size_type loc_inod = reference_element_e::ilat2loc_inod (order, ilat(i));
53 point xi = x0 + ((1.0*i)/order)*(x1 - x0);
54 node[dis_inod[loc_inod]] = project_on_unit_sphere (xi);
55 }
56 }
57 // 3) fix nodes inside elements
58 for (size_type ie = 0, ne = gamma.size(2); ie < ne; ie++) {
59 const geo_element& K = gamma.get_geo_element (2, ie);
60 gamma.dis_inod (K, dis_inod);
61 switch (K.variant()) {
63 const point& x0 = node[dis_inod[0]];
64 const point& x1 = node[dis_inod[1]];
65 const point& x2 = node[dis_inod[2]];
66 for (size_type i = 1; i < order; i++) {
67 for (size_type j = 1; i+j < order; j++) {
68 size_type loc_inod = reference_element_t::ilat2loc_inod (order, ilat(i, j));
69 point hat_x ((1.0*i)/order, (1.0*j)/order);
70 point xi = (1 - hat_x[0] - hat_x[1])*x0 + hat_x[0]*x1 + hat_x[1]*x2;
71 node[dis_inod[loc_inod]] = project_on_unit_sphere (xi);
72 }
73 }
74 break;
75 }
77 const point& x0 = node[dis_inod[0]];
78 const point& x1 = node[dis_inod[1]];
79 const point& x2 = node[dis_inod[2]];
80 const point& x3 = node[dis_inod[3]];
81 for (size_type i = 1; i < order; i++) {
82 for (size_type j = 1; j < order; j++) {
83 size_type loc_inod = reference_element_q::ilat2loc_inod (order, ilat(i, j));
84 point hat_x (-1+2.*i/order, -1+2.*j/order);
85 point xi = 0.25*((1 - hat_x[0])*(1 - hat_x[1])*x0
86 + (1 + hat_x[0])*(1 - hat_x[1])*x1
87 + (1 + hat_x[0])*(1 + hat_x[1])*x2
88 + (1 - hat_x[0])*(1 + hat_x[1])*x3);
89 node[dis_inod[loc_inod]] = project_on_unit_sphere (xi);
90 }
91 }
92 break;
93 }
94 default: error_macro ("element variant `"<<K.name()<<"' not yet supported");
95 }
96 }
97 geo_basic<Float,sequential> gamma_out = gamma;
98 gamma_out.set_nodes (node);
99 return gamma_out;
100}
103 typedef point_basic<size_type> ilat;
104 static const Float pi = acos(Float(-1.0));
105 size_t order = omega.order();
106 // 1) loop on the boundary and mark boundary faces, edges & vertices
107 const geo_basic<Float,sequential>& bdry = omega["boundary"];
108 const size_type d = omega.map_dimension();
109 const size_type sid_dim = bdry.map_dimension();
110 check_macro (d == sid_dim+1, "unexpected dimensions");
111 std::vector<bool> bdry_ver_mark (omega.n_node(), false);
112 std::vector<bool> bdry_edg_mark (omega.size(1), false);
113 std::vector<bool> bdry_sid_mark (omega.size(sid_dim), false);
114 for (size_type ioisid = 0, noisid = bdry.size(); ioisid < noisid; ioisid++) {
115 const geo_element& S = bdry.get_geo_element(sid_dim, ioisid);
116 bdry_sid_mark [S.dis_ie()] = true;
117 for (size_type loc_iv = 0, loc_nv = S.size(); loc_iv < loc_nv; loc_iv++) {
118 bdry_ver_mark [S[loc_iv]] = true;
119 }
120 if (d != 3) continue;
121 for (size_type loc_iedg = 0, loc_nedg = S.n_subgeo(1); loc_iedg < loc_nedg; loc_iedg++) {
122 bdry_edg_mark [S.edge(loc_iedg)] = true;
123 }
124 }
125 // 2) fix boundary vertices
126 disarray<point,sequential> node = omega.get_nodes();
127 for (size_type iv = 0, nv = omega.n_vertex(); iv < nv; iv++) {
128 if (! bdry_ver_mark [iv]) continue;
129 node [iv] = project_on_unit_sphere (node[iv]);
130 }
131 // 3) fix boundary edges
132 std::vector<size_type> dis_inod;
133 for (size_type iedg = 0, nedg = omega.size(1); iedg < nedg; iedg++) {
134 const geo_element& E = omega.get_geo_element(1, iedg);
135 omega.dis_inod (E, dis_inod);
136 const point& x0 = node[dis_inod[0]];
137 const point& x1 = node[dis_inod[1]];
138 for (size_type i = 1; i < order; i++) {
139 size_type loc_inod = reference_element_e::ilat2loc_inod (order, ilat(i));
140 point xi = x0 + ((1.0*i)/order)*(x1 - x0);
141 if ((d == 2 && bdry_sid_mark [iedg]) || (d == 3 && bdry_edg_mark [iedg])) {
142 xi = project_on_unit_sphere (xi);
143 }
144 node[dis_inod[loc_inod]] = xi;
145 }
146 }
147 // 4) fix boundary sides and internal sides with one boundary edge
148 if (d == 3) {
149 for (size_type ifac = 0, nfac = omega.size(2); ifac < nfac; ifac++) {
150 const geo_element& S = omega.get_geo_element(2, ifac);
151 omega.dis_inod (S, dis_inod);
152 bool side_has_boundary_edge = false;
153 size_type loc_bdry_iedg = 0;
154 if (! bdry_sid_mark [ifac]) {
155 size_type n_bdry_edg = 0;
156 for (size_type loc_iedg = 0, loc_nedg = S.n_subgeo (1); loc_iedg < loc_nedg; loc_iedg++) {
157 size_type iedg = S.edge(loc_iedg);
158 if (bdry_edg_mark [iedg]) {
159 n_bdry_edg++;
160 loc_bdry_iedg = loc_iedg;
161 side_has_boundary_edge = true;
162 }
163 }
164 check_macro (n_bdry_edg <= 1, "unsupported side with "<<n_bdry_edg<<" boundary edges");
165 }
166 switch (S.variant()) {
168 // side is internal and have a curved edge
169 const point& x0 = node[dis_inod[0]];
170 const point& x1 = node[dis_inod[1]];
171 const point& x2 = node[dis_inod[2]];
172 curved_ball_t<Float> F (x0, x1, x2, loc_bdry_iedg);
173 for (size_type i = 1; i < order; i++) {
174 for (size_type j = 1; i+j < order; j++) {
175 size_type loc_inod = reference_element_t::ilat2loc_inod (order, ilat(i, j));
176 point hat_x ((1.0*i)/order, (1.0*j)/order);
177 point xi;
178 if (side_has_boundary_edge) {
179 xi = F (hat_x);
180 } else {
181 xi = (1 - hat_x[0] - hat_x[1])*x0 + hat_x[0]*x1 + hat_x[1]*x2;
182 if (bdry_sid_mark [ifac]) {
183 xi = project_on_unit_sphere (xi);
184 }
185 }
186 node[dis_inod[loc_inod]] = xi;
187 }
188 }
189 break;
190 }
191#ifdef TODO
193 break;
194 }
195#endif // TODO
196 default: error_macro ("element variant `"<<S.name()<<"' not yet supported");
197 }
198 }
199 }
200 // 5) fix interior nodes of elements that have a curved boundary side
201 for (size_type ie = 0, ne = omega.size(); ie < ne; ie++) {
202 const geo_element& K = omega.get_geo_element(d, ie);
203 omega.dis_inod (K, dis_inod);
204 // 5.1) has K one side exactly on the boundary ?
205 size_type n_bdry_sid = 0;
206 size_type loc_bdry_isid = 0;
207 for (size_type loc_isid = 0, loc_nsid = K.n_subgeo (sid_dim); loc_isid < loc_nsid; loc_isid++) {
208 size_type isid = (sid_dim == 1) ? K.edge(loc_isid) : K.face(loc_isid);
209 if (bdry_sid_mark [isid]) {
210 n_bdry_sid++;
211 loc_bdry_isid = loc_isid;
212 }
213 }
214 check_macro (n_bdry_sid <= 1, "unsupported element with "<<n_bdry_sid<<" boundary sides");
215 bool is_on_bdry = (n_bdry_sid == 1); // n_bdry_sid == 1 i.e. K has exactly one boundary curved side
216 // 5.2) has K one edge exactly on the boundary ?
217 size_type n_bdry_edg = 0;
218 size_type loc_bdry_iedg = 0;
219 if (d == 3 && !is_on_bdry) {
220 for (size_type loc_iedg = 0, loc_nedg = K.n_subgeo (1); loc_iedg < loc_nedg; loc_iedg++) {
221 size_type iedg = K.edge(loc_iedg);
222 if (bdry_edg_mark [iedg]) {
223 n_bdry_edg++;
224 loc_bdry_iedg = loc_iedg;
225 }
226 }
227 check_macro (n_bdry_edg <= 1, "unsupported element with "<<n_bdry_edg<<" boundary edges");
228 }
229 bool has_bdry_edge = (n_bdry_edg == 1);
230 switch (K.variant()) {
232 // has a curved boundary edge:
233 const point& x0 = node[dis_inod[0]];
234 const point& x1 = node[dis_inod[1]];
235 const point& x2 = node[dis_inod[2]];
236 curved_ball_t<Float> F (x0, x1, x2, loc_bdry_isid);
237 for (size_type i = 1; i < order; i++) {
238 for (size_type j = 1; i+j < order; j++) {
239 size_type loc_inod = reference_element_t::ilat2loc_inod (order, ilat(i, j));
240 point hat_x ((1.0*i)/order, (1.0*j)/order);
241 point x;
242 if (is_on_bdry) {
243 x = F (hat_x);
244 } else {
245 x = (1 - hat_x[0] - hat_x[1])*x0 + hat_x[0]*x1 + hat_x[1]*x2;
246 }
247 node[dis_inod[loc_inod]] = x;
248 }
249 }
250 break;
251 }
253 if (is_on_bdry || has_bdry_edge) { // have a triangular face or an edge on the boundary
254 const point& x0 = node[dis_inod[0]];
255 const point& x1 = node[dis_inod[1]];
256 const point& x2 = node[dis_inod[2]];
257 const point& x3 = node[dis_inod[3]];
258 curved_ball_T<Float> F (x0, x1, x2, x3);
259 if (is_on_bdry) {
260 F.set_boundary_face (loc_bdry_isid);
261 } else if (has_bdry_edge) {
262 F.set_boundary_edge (loc_bdry_iedg);
263 }
264 for (size_type i = 1; i < order; i++) {
265 for (size_type j = 1; i+j < order; j++) {
266 for (size_type k = 1; i+j+k < order; k++) {
267 size_type loc_inod = reference_element_T::ilat2loc_inod (order, ilat(i, j, k));
268 point hat_x ((1.0*i)/order, (1.0*j)/order, (1.0*k)/order);
269 node[dis_inod[loc_inod]] = F (hat_x);
270 }
271 }
272 }
273 } else { // fully interior tetra:
274 const point& x0 = node[dis_inod[0]];
275 const point& x1 = node[dis_inod[1]];
276 const point& x2 = node[dis_inod[2]];
277 const point& x3 = node[dis_inod[3]];
278 for (size_type i = 1; i < order; i++) {
279 for (size_type j = 1; i+j < order; j++) {
280 for (size_type k = 1; i+j+k < order; k++) {
281 size_type loc_inod = reference_element_T::ilat2loc_inod (order, ilat(i, j, k));
282 point hat_x ((1.0*i)/order, (1.0*j)/order, (1.0*k)/order);
283 point x = (1 - hat_x[0] - hat_x[1] - hat_x[2])*x0 + hat_x[0]*x1 + hat_x[1]*x2 + hat_x[2]*x3;
284 node[dis_inod[loc_inod]] = x;
285 }
286 }
287 }
288 }
289 break;
290 }
292 size_type coord[4];
293 size_type shift[4];
294 shift[0] = loc_bdry_isid;
295 shift[1] = (loc_bdry_isid+1)%4;
296 shift[2] = (loc_bdry_isid+2)%4;
297 shift[3] = (loc_bdry_isid+3)%4;
298 // [x0,x1] is the curved boundary edge:
299 const point& x0 = node[dis_inod[shift[0]]];
300 const point& x1 = node[dis_inod[shift[1]]];
301 const point& x2 = node[dis_inod[shift[2]]];
302 const point& x3 = node[dis_inod[shift[3]]];
303 curved_ball_q<Float> F (x0, x1, x2, x3);
304 for (size_type i = 1; i < order; i++) {
305 for (size_type j = 1; j < order; j++) {
306 size_type loc_inod = reference_element_q::ilat2loc_inod (order, ilat(i, j));
307 coord [0] = i;
308 coord [1] = j;
309 coord [2] = order - i;
310 coord [3] = order - j;
311 size_type i1 = coord [shift[0]%4];
312 size_type j1 = coord [shift[1]%4];
313 point hat_x (-1+2.0*i1/order, -1+2.0*j1/order);
314 point x;
315 if (is_on_bdry) {
316 x = F (hat_x);
317 } else {
318 x = 0.25*( (1 - hat_x[0])*(1 - hat_x[1])*x0
319 + (1 + hat_x[0])*(1 - hat_x[1])*x1
320 + (1 + hat_x[0])*(1 + hat_x[1])*x2
321 + (1 - hat_x[0])*(1 + hat_x[1])*x3);
322 }
323 node[dis_inod[loc_inod]] = x;
324 }
325 }
326 break;
327 }
328 default: error_macro ("element variant `"<<K.name()<<"' not yet supported");
329 }
330 }
331 geo_basic<Float,sequential> omega_out = omega;
332 omega_out.set_nodes (node);
333 return omega_out;
334}
336 if (omega.map_dimension() < omega.dimension()) {
337 return fix_s (omega);
338 } else {
339 return fix_filled (omega);
340 }
341}
342void usage() {
343 std::cerr << "mkgeo_ball_gmsh_fix: usage:" << std::endl
344 << "mkgeo_ball_gmsh_fix "
345 << "[-order int] < input.geo > output.geo"
346 ;
347 exit (1);
348}
349int main(int argc, char**argv) {
350 environment distributed (argc, argv);
351 check_macro (communicator().size() == 1,
352 argv[0] << ": command may be used as mono-process only");
353 // ----------------------------
354 // scan the command line
355 // ----------------------------
356 size_t order = std::numeric_limits<size_t>::max();
357 for (int i = 1; i < argc; i++) {
358
359 if (strcmp (argv[i], "-order") == 0) {
360 if (i == argc-1) { std::cerr << argv[0] << " -order: option argument missing" << std::endl; usage(); }
361 order = atoi(argv[++i]);
362 }
363 else { usage(); }
364 }
365 // ----------------------------
366 // treatment
367 // ----------------------------
368 geo_basic<Float,sequential> omega; din >> omega;
369 if (order != std::numeric_limits<size_t>::max()) {
370 omega.reset_order (order);
371 }
372 omega = fix (omega);
373 dout << omega;
374}
field::size_type size_type
Definition branch.cc:430
see the Float page for the full documentation
see the communicator page for the full documentation
see the point page for the full documentation
void set_boundary_face(size_t loc_ifac_curved)
void set_boundary_edge(size_t loc_iedg_curved)
see the disarray page for the full documentation
Definition disarray.h:497
see the environment page for the full documentation
generic mesh with rerefence counting
Definition geo.h:1089
see the geo_element page for the full documentation
size_type edge(size_type i) const
size_type face(size_type i) const
size_type size() const
variant_type variant() const
size_type n_subgeo(size_type subgeo_dim) const
size_type dis_ie() const
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 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 q
static const variant_type T
static const variant_type t
#define error_macro(message)
Definition dis_macros.h:49
int main()
Definition field2bb.cc:58
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
geo_basic< Float, sequential > fix_filled(const geo_basic< Float, sequential > &omega)
void usage()
geo_basic< Float, sequential > fix_s(const geo_basic< Float, sequential > &gamma)
geo_basic< Float, sequential > fix(const geo_basic< Float, sequential > &omega)
point project_on_unit_sphere(const point &x)
This file is part of Rheolef.
T norm(const vec< T, M > &x)
norm(x): see the expression page for the full documentation
Definition vec.h:387
rheolef - reference manual