← Back to branch summary

~ubuntu-branches/ubuntu/trusty/rheolef/trusty

« back to all changes in this revision

Viewing changes to nfem/sbin/mkgeo_ball_gmsh_fix.cc

  • Committer: Package Import Robot
  • Author(s): Pierre Saramito
  • Date: 2012-04-06 09:12:21 UTC
  • mfrom: (1.1.5)
  • Revision ID: package-import@ubuntu.com-20120406091221-m58me99p1nxqui49
Tags: 6.0-1
http://bugs.debian.org/667356
http://bugs.debian.org/661689
* New upstream release 6.0 (major changes):
  - massively distributed and parallel support
  - full FEM characteristic method (Lagrange-Gakerkin method) support
  - enhanced users documentation 
  - source code supports g++-4.7 (closes: #667356)
* debian/control: dependencies for MPI distributed solvers added
* debian/rules: build commands simplified
* debian/librheolef-dev.install: man1/* to man9/* added
* debian/changelog: package description rewritted (closes: #661689)

Show diffs side-by-side

added added

removed removed

Lines of Context:
nfem/sbin/mkgeo_ball_gmsh_fix.cc
 
1
///
 
2
/// This file is part of Rheolef.
 
3
///
 
4
/// Copyright (C) 2000-2009 Pierre Saramito <Pierre.Saramito@imag.fr>
 
5
///
 
6
/// Rheolef is free software; you can redistribute it and/or modify
 
7
/// it under the terms of the GNU General Public License as published by
 
8
/// the Free Software Foundation; either version 2 of the License, or
 
9
/// (at your option) any later version.
 
10
///
 
11
/// Rheolef is distributed in the hope that it will be useful,
 
12
/// but WITHOUT ANY WARRANTY; without even the implied warranty of
 
13
/// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 
14
/// GNU General Public License for more details.
 
15
///
 
16
/// You should have received a copy of the GNU General Public License
 
17
/// along with Rheolef; if not, write to the Free Software
 
18
/// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 
19
///
 
20
/// =========================================================================
 
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
using namespace rheolef;
 
32
using namespace std;
 
33
point project_on_unit_sphere (const point& x) { return x/norm(x); }
 
34
geo_basic<Float,sequential> fix_s (const geo_basic<Float,sequential>& gamma)
 
35
{
 
36
  typedef geo_basic<Float,sequential>::size_type size_type;
 
37
  typedef point_basic<size_type>                 ilat;
 
38
  size_t order = gamma.order();
 
39
  array<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()) {
 
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
  }
 
97
  geo_basic<Float,sequential> gamma_out = gamma;
 
98
  gamma_out.set_nodes (node);
 
99
  return gamma_out;
 
100
}
 
101
geo_basic<Float,sequential> fix_filled (const geo_basic<Float,sequential>& omega) {
 
102
  typedef geo_basic<Float,sequential>::size_type size_type;
 
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
  array<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 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 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 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 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
}
 
335
geo_basic<Float,sequential> fix (const geo_basic<Float,sequential>& omega) {
 
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
  cerr << "mkgeo_ball_gmsh_fix: usage:" << 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) { cerr << argv[0] << " -order: option argument missing" << 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
}