Rheolef  7.1
an efficient C++ finite element environment
mkgeo_ball_gmsh_fix.cc
Go to the documentation of this file.
1 // 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>
32 using namespace rheolef;
33 point 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++) {
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()) {
62  case reference_element::t: {
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  }
76  case reference_element::q: {
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  }
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()) {
167  case reference_element::t: {
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
192  case reference_element::q: {
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()) {
231  case reference_element::t: {
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  }
252  case reference_element::T: {
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  }
291  case reference_element::q: {
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 }
342 void 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 }
349 int 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 }
rheolef::reference_element_t::ilat2loc_inod
static size_type ilat2loc_inod(size_type order, const point_basic< size_type > &ilat)
Definition: reference_element.cc:413
rheolef::geo_basic
generic mesh with rerefence counting
Definition: geo.h:1089
main
int main(int argc, char **argv)
Definition: mkgeo_ball_gmsh_fix.cc:349
usage
void usage()
Definition: mkgeo_ball_gmsh_fix.cc:342
rheolef::point_basic
Definition: point.h:87
check_macro
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
rheolef::reference_element_e::ilat2loc_inod
static size_type ilat2loc_inod(size_type order, const point_basic< size_type > &ilat)
Definition: reference_element.cc:350
rheolef::reference_element_T::ilat2loc_inod
static size_type ilat2loc_inod(size_type order, const point_basic< size_type > &ilat)
Definition: reference_element.cc:582
fix
geo_basic< Float, sequential > fix(const geo_basic< Float, sequential > &omega)
Definition: mkgeo_ball_gmsh_fix.cc:335
rheolef::curved_ball_T::set_boundary_edge
void set_boundary_edge(size_t loc_iedg_curved)
Definition: geo_element_curved_ball.h:205
project_on_unit_sphere
point project_on_unit_sphere(const point &x)
Definition: mkgeo_ball_gmsh_fix.cc:33
rheolef::curved_ball_q
Definition: geo_element_curved_ball.h:55
rheolef::reference_element::T
static const variant_type T
Definition: reference_element.h:79
mkgeo_ball.order
order
Definition: mkgeo_ball.sh:343
rheolef::geo_element::n_subgeo
size_type n_subgeo(size_type subgeo_dim) const
Definition: geo_element.h:212
rheolef::geo_element::size
size_type size() const
Definition: geo_element.h:168
mkgeo_ball.d
int d
Definition: mkgeo_ball.sh:154
rheolef::geo_element
see the geo_element page for the full documentation
Definition: geo_element.h:102
rheolef::curved_ball_T
Definition: geo_element_curved_ball.h:174
rheolef::size_type
size_t size_type
Definition: basis_get.cc:76
rk::gamma
Float gamma[][pmax+1]
Definition: runge_kutta_semiimplicit.icc:70
rheolef::geo_element::variant
variant_type variant() const
Definition: geo_element.h:161
rheolef::norm
T norm(const vec< T, M > &x)
norm(x): see the expression page for the full documentation
Definition: vec.h:387
rheolef::geo_element::dis_ie
size_type dis_ie() const
Definition: geo_element.h:163
rheolef.h
rheolef - reference manual
rheolef::curved_ball_t
Definition: geo_element_curved_ball.h:111
rheolef::reference_element_q::ilat2loc_inod
static size_type ilat2loc_inod(size_type order, const point_basic< size_type > &ilat)
Definition: reference_element.cc:481
rheolef::din
idiststream din
see the diststream page for the full documentation
Definition: diststream.h:427
rheolef::environment
see the environment page for the full documentation
Definition: environment.h:115
rheolef
This file is part of Rheolef.
Definition: compiler_eigen.h:37
error_macro
#define error_macro(message)
Definition: dis_macros.h:49
rheolef::curved_ball_T::set_boundary_face
void set_boundary_face(size_t loc_ifac_curved)
Definition: geo_element_curved_ball.h:190
rheolef::space_numbering::dis_inod
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)
Definition: space_numbering.cc:177
Float
see the Float page for the full documentation
fix_filled
geo_basic< Float, sequential > fix_filled(const geo_basic< Float, sequential > &omega)
Definition: mkgeo_ball_gmsh_fix.cc:101
fix_s
geo_basic< Float, sequential > fix_s(const geo_basic< Float, sequential > &gamma)
Definition: mkgeo_ball_gmsh_fix.cc:34
point
see the point page for the full documentation
rheolef::disarray
see the disarray page for the full documentation
Definition: disarray.h:459
rheolef::reference_element::q
static const variant_type q
Definition: reference_element.h:78
rheolef::reference_element::t
static const variant_type t
Definition: reference_element.h:77
size_type
field::size_type size_type
Definition: branch.cc:425
rheolef::geo_element::name
char name() const
Definition: geo_element.h:169
rheolef::Float
double Float
see the Float page for the full documentation
Definition: Float.h:143
rheolef::geo_element::face
size_type face(size_type i) const
Definition: geo_element.h:210
rheolef::dout
odiststream dout(cout)
see the diststream page for the full documentation
Definition: diststream.h:430
rheolef::distributed
distributed
Definition: asr.cc:228
rheolef::geo_element::edge
size_type edge(size_type i) const
Definition: geo_element.h:209