gmsh-TingyuanDoc  0.1
An Open-Source Timing-driven Analytical Mixed-size FPGA Placer
meshGRegionTransfinite.cpp
Go to the documentation of this file.
1 // Gmsh - Copyright (C) 1997-2022 C. Geuzaine, J.-F. Remacle
2 //
3 // See the LICENSE.txt file in the Gmsh root directory for license information.
4 // Please report all issues on https://gitlab.onelab.info/gmsh/gmsh/issues.
5 //
6 // Contributor(s):
7 // Trevor S. Strickler
8 // Michael Ermakov (ermakov@ipmnet.ru)
9 //
10 
11 #include <map>
12 #include "GmshConfig.h"
13 #include "GmshMessage.h"
14 #include "meshGFace.h"
15 #include "GFace.h"
16 #include "GRegion.h"
17 #include "MVertex.h"
18 #include "MTetrahedron.h"
19 #include "MHexahedron.h"
20 #include "MPrism.h"
21 #include "Context.h"
22 
23 #if defined(HAVE_QUADTRI)
24 #include "QuadTriTransfinite3D.h"
25 #endif
26 
27 /*
28  Transfinite volume meshes
29 
30  a0 s0 s1 f0 s0 s1 s5 s4 s6
31  s7 s6 a1 s1 s2 f1 s1 s2 s6 s5 *
32  *-------* a2 s3 s2 f2 s3 s2 s6 s7 /|\
33  |\s4 |\ a3 s0 s3 f3 s0 s3 s7 s4 / | \
34  | *-------* s5 a4 s4 s5 f4 s0 s1 s2 s3 s7/s4/ |s2\
35  | | s2| | a5 s5 s6 f5 s4 s5 s6 s7 *---*---* s5
36  s3 *-|-----* | a6 s7 s6 | / \ |
37  \| \| a7 s4 s7 | / \ |
38  *-------* a8 s0 s4 |/ \|
39  v w s0 s1 a9 s1 s5 *-------*
40  \| a10 s2 s6 v w s3/s0 s1
41  *--u a11 s3 s7 \|
42  *--u
43  Important limitations:
44 
45  - only works with 5- or 6-face volumes
46 
47  - the faces have to be meshed with the 2D transfinite algorithm
48 
49  - the definition of a 5-face volume has to follow the ordering given
50  in the figure above (degenerescence has to be along s0,s4)
51 
52  - meshing a volume with prisms or tetrahedra assumes that the
53  triangular mesh is consistent with the volume mesh: there is no
54  coherence check in the volume algorithm to ensure that edges will
55  match.
56 */
57 
58 #define CREATE_HEX \
59  new MHexahedron(tab[i][j][k], tab[i + 1][j][k], tab[i + 1][j + 1][k], \
60  tab[i][j + 1][k], tab[i][j][k + 1], tab[i + 1][j][k + 1], \
61  tab[i + 1][j + 1][k + 1], tab[i][j + 1][k + 1])
62 
63 #define CREATE_PRISM_1 \
64  new MPrism(tab[i][j][k], tab[i + 1][j][k], tab[i][j + 1][k], \
65  tab[i][j][k + 1], tab[i + 1][j][k + 1], tab[i][j + 1][k + 1])
66 
67 #define CREATE_PRISM_2 \
68  new MPrism(tab[i + 1][j + 1][k], tab[i][j + 1][k], tab[i + 1][j][k], \
69  tab[i + 1][j + 1][k + 1], tab[i][j + 1][k + 1], \
70  tab[i + 1][j][k + 1])
71 
72 #define CREATE_PRISM_3 \
73  new MPrism(tab[i][j][k], tab[i][j + 1][k], tab[i + 1][j + 1][k], \
74  tab[i][j][k + 1], tab[i][j + 1][k + 1], tab[i + 1][j + 1][k + 1])
75 
76 #define CREATE_PRISM_4 \
77  new MPrism(tab[i][j][k], tab[i + 1][j][k], tab[i + 1][j + 1][k], \
78  tab[i][j][k + 1], tab[i + 1][j][k + 1], tab[i + 1][j + 1][k + 1])
79 
80 #define CREATE_SIM_1 \
81  new MTetrahedron(tab[i][j][k], tab[i + 1][j][k], tab[i][j + 1][k], \
82  tab[i][j][k + 1])
83 
84 #define CREATE_SIM_2 \
85  new MTetrahedron(tab[i + 1][j][k], tab[i][j + 1][k], tab[i][j][k + 1], \
86  tab[i + 1][j][k + 1])
87 
88 #define CREATE_SIM_3 \
89  new MTetrahedron(tab[i][j][k + 1], tab[i + 1][j][k + 1], tab[i][j + 1][k], \
90  tab[i][j + 1][k + 1])
91 
92 #define CREATE_SIM_4 \
93  new MTetrahedron(tab[i + 1][j][k], tab[i][j + 1][k], tab[i + 1][j][k + 1], \
94  tab[i + 1][j + 1][k])
95 
96 #define CREATE_SIM_5 \
97  new MTetrahedron(tab[i][j + 1][k], tab[i][j + 1][k + 1], \
98  tab[i + 1][j][k + 1], tab[i + 1][j + 1][k])
99 
100 #define CREATE_SIM_6 \
101  new MTetrahedron(tab[i + 1][j][k + 1], tab[i][j + 1][k + 1], \
102  tab[i + 1][j + 1][k + 1], tab[i + 1][j + 1][k])
103 
104 #define CREATE_SIM_7 \
105  new MTetrahedron(tab[i][j][k], tab[i][j + 1][k], tab[i][j][k + 1], \
106  tab[i + 1][j + 1][k])
107 
108 #define CREATE_SIM_8 \
109  new MTetrahedron(tab[i][j + 1][k], tab[i][j + 1][k + 1], tab[i][j][k + 1], \
110  tab[i + 1][j + 1][k])
111 
112 #define CREATE_SIM_9 \
113  new MTetrahedron(tab[i][j][k + 1], tab[i][j + 1][k + 1], \
114  tab[i + 1][j + 1][k + 1], tab[i + 1][j + 1][k])
115 
116 #define CREATE_SIM_10 \
117  new MTetrahedron(tab[i][j][k], tab[i + 1][j][k], tab[i + 1][j + 1][k], \
118  tab[i][j][k + 1])
119 
120 #define CREATE_SIM_11 \
121  new MTetrahedron(tab[i + 1][j][k], tab[i + 1][j + 1][k], tab[i][j][k + 1], \
122  tab[i + 1][j][k + 1])
123 
124 #define CREATE_SIM_12 \
125  new MTetrahedron(tab[i][j][k + 1], tab[i + 1][j][k + 1], \
126  tab[i + 1][j + 1][k], tab[i + 1][j + 1][k + 1])
127 
128 static double transfiniteHex(double f1, double f2, double f3, double f4,
129  double f5, double f6, double c1, double c2,
130  double c3, double c4, double c5, double c6,
131  double c7, double c8, double c9, double c10,
132  double c11, double c12, double s1, double s2,
133  double s3, double s4, double s5, double s6,
134  double s7, double s8, double u, double v, double w)
135 {
136  return (1 - u) * f4 + u * f2 + (1 - v) * f1 + v * f3 + (1 - w) * f5 + w * f6 -
137  ((1 - u) * (1 - v) * c9 + (1 - u) * v * c12 + u * (1 - v) * c10 +
138  u * v * c11) -
139  ((1 - v) * (1 - w) * c1 + (1 - v) * w * c5 + v * (1 - w) * c3 +
140  v * w * c7) -
141  ((1 - u) * (1 - w) * c4 + (1 - w) * u * c2 + w * (1 - u) * c8 +
142  u * w * c6) +
143  (1 - u) * (1 - v) * (1 - w) * s1 + u * (1 - v) * (1 - w) * s2 +
144  u * v * (1 - w) * s3 + (1 - u) * v * (1 - w) * s4 +
145  (1 - u) * (1 - v) * w * s5 + u * (1 - v) * w * s6 + u * v * w * s7 +
146  (1 - u) * v * w * s8;
147 }
148 
149 static MVertex *
150 transfiniteHex(GRegion *gr, MVertex *f1, MVertex *f2, MVertex *f3, MVertex *f4,
151  MVertex *f5, MVertex *f6, MVertex *c1, MVertex *c2, MVertex *c3,
152  MVertex *c4, MVertex *c5, MVertex *c6, MVertex *c7, MVertex *c8,
153  MVertex *c9, MVertex *c10, MVertex *c11, MVertex *c12,
154  MVertex *s1, MVertex *s2, MVertex *s3, MVertex *s4, MVertex *s5,
155  MVertex *s6, MVertex *s7, MVertex *s8, double u, double v,
156  double w)
157 {
158  double x = transfiniteHex(
159  f1->x(), f2->x(), f3->x(), f4->x(), f5->x(), f6->x(), c1->x(), c2->x(),
160  c3->x(), c4->x(), c5->x(), c6->x(), c7->x(), c8->x(), c9->x(), c10->x(),
161  c11->x(), c12->x(), s1->x(), s2->x(), s3->x(), s4->x(), s5->x(), s6->x(),
162  s7->x(), s8->x(), u, v, w);
163  double y = transfiniteHex(
164  f1->y(), f2->y(), f3->y(), f4->y(), f5->y(), f6->y(), c1->y(), c2->y(),
165  c3->y(), c4->y(), c5->y(), c6->y(), c7->y(), c8->y(), c9->y(), c10->y(),
166  c11->y(), c12->y(), s1->y(), s2->y(), s3->y(), s4->y(), s5->y(), s6->y(),
167  s7->y(), s8->y(), u, v, w);
168  double z = transfiniteHex(
169  f1->z(), f2->z(), f3->z(), f4->z(), f5->z(), f6->z(), c1->z(), c2->z(),
170  c3->z(), c4->z(), c5->z(), c6->z(), c7->z(), c8->z(), c9->z(), c10->z(),
171  c11->z(), c12->z(), s1->z(), s2->z(), s3->z(), s4->z(), s5->z(), s6->z(),
172  s7->z(), s8->z(), u, v, w);
173  return new MVertex(x, y, z, gr);
174 }
175 
176 class GOrientedTransfiniteFace {
177 private:
178  GFace *_gf;
179  int _ll, _hh;
180  int _permutation, _index;
181  std::vector<MVertex *> _list;
182 
183 public:
184  GOrientedTransfiniteFace()
185  : _gf(nullptr), _ll(0), _hh(0), _permutation(-1), _index(-1)
186  {
187  }
188  GOrientedTransfiniteFace(GFace *gf, std::vector<MVertex *> &corners)
189  : _gf(gf), _ll(0), _hh(0), _permutation(-1), _index(-1)
190  {
191  _ll = gf->transfinite_vertices.size() - 1;
192  if(_ll <= 0) return;
193  _hh = gf->transfinite_vertices[0].size() - 1;
194  if(_hh <= 0) return;
195  Msg::Debug("Face %d: L = %d H = %d", gf->tag(), _ll, _hh);
196 
197  // get the corners of the transfinite volume interpolation
198  std::vector<MVertex *> s(8);
199  if(corners.size() == 8) {
200  for(int i = 0; i < 8; i++) s[i] = corners[i];
201  }
202  else if(corners.size() == 6) {
203  s[0] = corners[0];
204  s[1] = corners[1];
205  s[2] = corners[2];
206  s[3] = corners[0];
207  s[4] = corners[3];
208  s[5] = corners[4];
209  s[6] = corners[5];
210  s[7] = corners[3];
211  }
212  else
213  return;
214 
215  // get the corners of the transfinite surface mesh
216  std::vector<MVertex *> c(4);
217  if(_gf->meshAttributes.corners.empty() ||
218  _gf->meshAttributes.corners.size() == 4) {
219  c[0] = _gf->transfinite_vertices[0][0];
220  c[1] = _gf->transfinite_vertices[_ll][0];
221  c[2] = _gf->transfinite_vertices[_ll][_hh];
222  c[3] = _gf->transfinite_vertices[0][_hh];
223  }
224  else if(_gf->meshAttributes.corners.size() == 3) {
225  c[0] = _gf->transfinite_vertices[0][0];
226  c[1] = _gf->transfinite_vertices[_ll][0];
227  c[2] = _gf->transfinite_vertices[_ll][_hh];
228  c[3] = _gf->transfinite_vertices[0][0];
229  }
230  else
231  return;
232 
233  // map the surface mesh onto the canonical transfinite hexahedron
234  int faces[] = {0, 1, 5, 4, 1, 2, 6, 5, 3, 2, 6, 7,
235  0, 3, 7, 4, 0, 1, 2, 3, 4, 5, 6, 7};
236  int permutations[] = {0, 1, 2, 3, 1, 2, 3, 0, 2, 3, 0, 1, 3, 0, 1, 2,
237  3, 2, 1, 0, 2, 1, 0, 3, 1, 0, 3, 2, 0, 3, 2, 1};
238  for(int p = 0; p < 8; p++) {
239  for(int f = 0; f < 6; f++) {
240  if(s[faces[4 * f + 0]] == c[permutations[4 * p + 0]] &&
241  s[faces[4 * f + 1]] == c[permutations[4 * p + 1]] &&
242  s[faces[4 * f + 2]] == c[permutations[4 * p + 2]] &&
243  s[faces[4 * f + 3]] == c[permutations[4 * p + 3]]) {
244  _index = f;
245  _permutation = p;
246  break;
247  }
248  }
249  }
250  Msg::Debug("Found face index %d (permutation = %d)", _index, _permutation);
251  for(int i = 0; i <= _ll; i++)
252  for(int j = 0; j <= _hh; j++)
253  _list.push_back(_gf->transfinite_vertices[i][j]);
254  }
255  // returns the surface
256  GFace *getSurface() { return _gf; }
257  // returns the index of the face in the reference hexahedron
258  int index() const { return _index; }
259  // returns true if the face is recombined
260  int recombined() const
261  {
262  return CTX::instance()->mesh.recombineAll || _gf->meshAttributes.recombine;
263  }
264  // returns the number or points in the transfinite mesh in both
265  // parameter directions
266  int getNumU() { return (_permutation % 2) ? _hh + 1 : _ll + 1; }
267  int getNumV() { return (_permutation % 2) ? _ll + 1 : _hh + 1; }
268  // returns the (i,j) vertex in the face, i and j being defined in
269  // the coordinate system of the reference transfinite hexahedron
270  MVertex *getVertex(int i, int j)
271  {
272  int index = -1, m = i, n = j, M = getNumU(), N = getNumV();
273  switch(_permutation) {
274  case 0: index = (n + N * m); break;
275  case 1: index = (M * N - M * (n + 1) + m); break;
276  case 2: index = (M * N - (n + N * m) - 1); break;
277  case 3: index = (M + n * M - m - 1); break;
278  case 4: index = (N + m * N - n - 1); break;
279  case 5: index = (M * N - (m + M * n) - 1); break;
280  case 6: index = (M * N - N * (m + 1) + n); break;
281  case 7: index = (m + M * n); break;
282  }
283  MVertex *v = nullptr;
284  if(index >= 0 && index < (int)_list.size()) v = _list[index];
285  if(index < 0 || index >= (int)_list.size() || !v) {
286  Msg::Error("Wrong index in transfinite mesh of surface %d: "
287  "m=%d n=%d M=%d N=%d perm=%d",
288  _gf->tag(), m, n, M, N, _permutation);
289  return _list[0];
290  }
291  return v;
292  }
293 };
294 
295 void findTransfiniteCorners(GRegion *gr, std::vector<MVertex *> &corners)
296 {
297  if(gr->meshAttributes.corners.size()) {
298  // corners have been specified explicitly
299  for(std::size_t i = 0; i < gr->meshAttributes.corners.size(); i++)
300  corners.push_back(gr->meshAttributes.corners[i]->mesh_vertices[0]);
301  }
302  else {
303  // try to find the corners automatically
304  std::vector<GFace *> faces = gr->faces();
305  GFace *gf = nullptr;
306  if(faces.size() == 6) {
307  // any face will do as a starting face
308  gf = faces.front();
309  }
310  else if(faces.size() == 5) {
311  // we need to start with a triangular face
312  auto found_it = std::find_if(begin(faces), end(faces), [](GFace *face) {
313  return face->edges().size() == 3 ||
314  face->meshAttributes.corners.size() == 3;
315  });
316  if(found_it != end(faces)) { gf = *found_it; }
317  }
318  if(gf) {
319  std::vector<GEdge *> redges = gr->edges();
320  for(auto *fedge : gf->edges()) {
321  const auto found_it = std::find(begin(redges), end(redges), fedge);
322  if(found_it != end(redges)) { redges.erase(found_it); }
323  }
324  findTransfiniteCorners(gf, corners);
325  std::size_t N = corners.size();
326  for(std::size_t i = 0; i < N; i++) {
327  for(auto it = redges.begin(); it != redges.end(); it++) {
328  if((*it)->getBeginVertex()->mesh_vertices[0] == corners[i]) {
329  corners.push_back((*it)->getEndVertex()->mesh_vertices[0]);
330  break;
331  }
332  else if((*it)->getEndVertex()->mesh_vertices[0] == corners[i]) {
333  corners.push_back((*it)->getBeginVertex()->mesh_vertices[0]);
334  break;
335  }
336  }
337  }
338  }
339  }
340 }
341 
342 int MeshTransfiniteVolume(GRegion *gr)
343 {
344  if(gr->meshAttributes.method != MESH_TRANSFINITE) return 0;
345 
346  Msg::Info("Meshing volume %d (Transfinite)", gr->tag());
347 
348  std::vector<GFace *> faces = gr->faces();
349  if(faces.size() != 5 && faces.size() != 6) {
350  Msg::Error(
351  "Transfinite algorithm only available for 5- and 6-face volumes");
352  return 0;
353  }
354 
355  // make sure that all bounding edges have begin/end points: everything in here
356  // depends on it
357  const std::vector<GEdge *> &edges = gr->edges();
358  for(auto it = edges.begin(); it != edges.end(); it++) {
359  if(!(*it)->getBeginVertex() || !(*it)->getEndVertex()) {
360  Msg::Error("Transfinite algorithm cannot be applied with curve %d which "
361  "has no begin or end point",
362  (*it)->tag());
363  return 0;
364  }
365  }
366 
367  std::vector<MVertex *> corners;
368  findTransfiniteCorners(gr, corners);
369  if(corners.size() != 6 && corners.size() != 8) {
370  Msg::Error("Volume %d is transfinite but has %d corners", gr->tag(),
371  corners.size());
372  return 0;
373  }
374 
375  std::vector<GOrientedTransfiniteFace> orientedFaces(6);
376  for(auto it = faces.begin(); it != faces.end(); ++it) {
377  GOrientedTransfiniteFace f(*it, corners);
378  if(f.index() < 0) {
379  Msg::Error("Incompatible surface %d in transfinite volume %d",
380  (*it)->tag(), gr->tag());
381  return 0;
382  }
383  else if(orientedFaces[f.index()].getSurface()) {
384  Msg::Error("Surfaces %d and %d mapped to the same index in transfinite "
385  "volume %d",
386  orientedFaces[f.index()].getSurface()->tag(), (*it)->tag(),
387  gr->tag());
388  return 0;
389  }
390  orientedFaces[f.index()] = f;
391  }
392 
393  int N_i = orientedFaces[4].getNumU();
394  int N_j = orientedFaces[4].getNumV();
395  int N_k = orientedFaces[1].getNumV();
396 
397  std::vector<double> lengths_i, lengths_j, lengths_k;
398  double L_i = 0., L_j = 0., L_k = 0.;
399  lengths_i.push_back(0.);
400  lengths_j.push_back(0.);
401  lengths_k.push_back(0.);
402  for(int i = 0; i < N_i - 1; i++) {
403  MVertex *v1 = orientedFaces[4].getVertex(i, 0);
404  MVertex *v2 = orientedFaces[4].getVertex(i + 1, 0);
405  L_i += v1->distance(v2);
406  lengths_i.push_back(L_i);
407  }
408  for(int i = 0; i < N_j - 1; i++) {
409  MVertex *v1 = orientedFaces[1].getVertex(i, 0);
410  MVertex *v2 = orientedFaces[1].getVertex(i + 1, 0);
411  L_j += v1->distance(v2);
412  lengths_j.push_back(L_j);
413  }
414  for(int i = 0; i < N_k - 1; i++) {
415  MVertex *v1 = orientedFaces[1].getVertex(0, i);
416  MVertex *v2 = orientedFaces[1].getVertex(0, i + 1);
417  L_k += v1->distance(v2);
418  lengths_k.push_back(L_k);
419  }
420 
421  // create points using transfinite interpolation
422 
423  MVertex *s0 = orientedFaces[4].getVertex(0, 0);
424  MVertex *s1 = orientedFaces[4].getVertex(N_i - 1, 0);
425  MVertex *s2 = orientedFaces[4].getVertex(N_i - 1, N_j - 1);
426  MVertex *s3 = orientedFaces[4].getVertex(0, N_j - 1);
427 
428  MVertex *s4 = orientedFaces[5].getVertex(0, 0);
429  MVertex *s5 = orientedFaces[5].getVertex(N_i - 1, 0);
430  MVertex *s6 = orientedFaces[5].getVertex(N_i - 1, N_j - 1);
431  MVertex *s7 = orientedFaces[5].getVertex(0, N_j - 1);
432 
433  std::vector<std::vector<std::vector<MVertex *> > > &tab(
434  gr->transfinite_vertices);
435  tab.resize(N_i);
436  for(int i = 0; i < N_i; i++) {
437  tab[i].resize(N_j);
438  for(int j = 0; j < N_j; j++) { tab[i][j].resize(N_k); }
439  }
440 
441  for(int i = 0; i < N_i; i++) {
442  double u = lengths_i[i] / L_i;
443 
444  for(int j = 0; j < N_j; j++) {
445  double v = lengths_j[j] / L_j;
446 
447  MVertex *c0 = orientedFaces[4].getVertex(i, 0);
448  MVertex *c1 = orientedFaces[4].getVertex(N_i - 1, j);
449  MVertex *c2 = orientedFaces[4].getVertex(i, N_j - 1);
450  MVertex *c3 = orientedFaces[4].getVertex(0, j);
451 
452  MVertex *c4 = orientedFaces[5].getVertex(i, 0);
453  MVertex *c5 = orientedFaces[5].getVertex(N_i - 1, j);
454  MVertex *c6 = orientedFaces[5].getVertex(i, N_j - 1);
455  MVertex *c7 = orientedFaces[5].getVertex(0, j);
456 
457  MVertex *f4 = orientedFaces[4].getVertex(i, j);
458  MVertex *f5 = orientedFaces[5].getVertex(i, j);
459 
460  for(int k = 0; k < N_k; k++) {
461  double w = lengths_k[k] / L_k;
462 
463  MVertex *c8 = orientedFaces[0].getVertex(0, k);
464  MVertex *c9 = orientedFaces[0].getVertex(N_i - 1, k);
465  MVertex *c10 = orientedFaces[2].getVertex(N_i - 1, k);
466  MVertex *c11 = orientedFaces[2].getVertex(0, k);
467 
468  MVertex *f0 = orientedFaces[0].getVertex(i, k);
469  MVertex *f1 = orientedFaces[1].getVertex(j, k);
470  MVertex *f2 = orientedFaces[2].getVertex(i, k);
471  MVertex *f3;
472  if(corners.size() == 8)
473  f3 = orientedFaces[3].getVertex(j, k);
474  else
475  f3 = c8;
476 
477  if(i && j && k && i != N_i - 1 && j != N_j - 1 && k != N_k - 1) {
478  MVertex *newv = transfiniteHex(
479  gr, f0, f1, f2, f3, f4, f5, c0, c1, c2, c3, c4, c5, c6, c7, c8, c9,
480  c10, c11, s0, s1, s2, s3, s4, s5, s6, s7, u, v, w);
481  gr->mesh_vertices.push_back(newv);
482  tab[i][j][k] = newv;
483  }
484  else if(!i) {
485  tab[i][j][k] = f3;
486  }
487  else if(!j) {
488  tab[i][j][k] = f0;
489  }
490  else if(!k) {
491  tab[i][j][k] = f4;
492  }
493  else if(i == N_i - 1) {
494  tab[i][j][k] = f1;
495  }
496  else if(j == N_j - 1) {
497  tab[i][j][k] = f2;
498  }
499  else if(k == N_k - 1) {
500  tab[i][j][k] = f5;
501  }
502  }
503  }
504  }
505 
506 #if defined(HAVE_QUADTRI)
507  // for QuadTri, get external boundary diagonals for element subdivision
508  // purposes
509  std::set<std::pair<MVertex *, MVertex *> > boundary_diags;
510  if(gr->meshAttributes.QuadTri) {
511  if(!getTransfiniteBoundaryDiags(gr, &boundary_diags)) {
512  Msg::Error(
513  "In MeshTransfiniteVolume(), getTransfiniteBoundaryDiags() failed. "
514  "Aborting mesh of region %d.",
515  gr->tag());
516  return 0;
517  }
518  }
519 #endif
520 
521  // create elements
522 
523  if(faces.size() == 6) {
524  for(int i = 0; i < N_i - 1; i++) {
525  for(int j = 0; j < N_j - 1; j++) {
526  for(int k = 0; k < N_k - 1; k++) {
527 #if defined(HAVE_QUADTRI)
528  if(gr->meshAttributes.QuadTri) {
529  // create vertex array
530  std::vector<MVertex *> verts;
531  verts.resize(8);
532  verts[0] = tab[i][j][k];
533  verts[1] = tab[i + 1][j][k];
534  verts[2] = tab[i + 1][j + 1][k];
535  verts[3] = tab[i][j + 1][k];
536  verts[4] = tab[i][j][k + 1];
537  verts[5] = tab[i + 1][j][k + 1];
538  verts[6] = tab[i + 1][j + 1][k + 1];
539  verts[7] = tab[i][j + 1][k + 1];
540  if((!orientedFaces[3].recombined() && i == 0) ||
541  (!orientedFaces[1].recombined() && i == N_i - 2) ||
542  (!orientedFaces[0].recombined() && j == 0) ||
543  (!orientedFaces[2].recombined() && j == N_j - 2) ||
544  (!orientedFaces[4].recombined() && k == 0) ||
545  (!orientedFaces[5].recombined() && k == N_k - 2)) {
546  // make subdivided element
547  meshTransfElemWithInternalVertex(gr, verts, &boundary_diags);
548  }
549  // if not adjacent to unrecombined edge
550  else
551  gr->hexahedra.push_back(
552  new MHexahedron(verts[0], verts[1], verts[2], verts[3],
553  verts[4], verts[5], verts[6], verts[7]));
554  // continue, skipping the rest which is for non-divided elements
555  continue;
556  }
557 #endif
558  if(orientedFaces[0].recombined() && orientedFaces[1].recombined() &&
559  orientedFaces[2].recombined() && orientedFaces[3].recombined() &&
560  orientedFaces[4].recombined() && orientedFaces[5].recombined()) {
561  gr->hexahedra.push_back(CREATE_HEX);
562  }
563  else if(!orientedFaces[0].recombined() &&
564  orientedFaces[1].recombined() &&
565  !orientedFaces[2].recombined() &&
566  orientedFaces[3].recombined() &&
567  orientedFaces[4].recombined() &&
568  orientedFaces[5].recombined()) {
569  gr->prisms.push_back(new MPrism(
570  tab[i][j][k], tab[i + 1][j][k], tab[i][j][k + 1],
571  tab[i][j + 1][k], tab[i + 1][j + 1][k], tab[i][j + 1][k + 1]));
572  gr->prisms.push_back(
573  new MPrism(tab[i + 1][j][k + 1], tab[i][j][k + 1],
574  tab[i + 1][j][k], tab[i + 1][j + 1][k + 1],
575  tab[i][j + 1][k + 1], tab[i + 1][j + 1][k]));
576  }
577  else if(orientedFaces[0].recombined() &&
578  !orientedFaces[1].recombined() &&
579  orientedFaces[2].recombined() &&
580  !orientedFaces[3].recombined() &&
581  orientedFaces[4].recombined() &&
582  orientedFaces[5].recombined()) {
583  gr->prisms.push_back(new MPrism(
584  tab[i + 1][j][k], tab[i + 1][j + 1][k], tab[i + 1][j][k + 1],
585  tab[i][j][k], tab[i][j + 1][k], tab[i][j][k + 1]));
586  gr->prisms.push_back(
587  new MPrism(tab[i + 1][j + 1][k + 1], tab[i + 1][j][k + 1],
588  tab[i + 1][j + 1][k], tab[i][j + 1][k + 1],
589  tab[i][j][k + 1], tab[i][j + 1][k]));
590  }
591  else if(orientedFaces[0].recombined() &&
592  orientedFaces[1].recombined() &&
593  orientedFaces[2].recombined() &&
594  orientedFaces[3].recombined() &&
595  !orientedFaces[4].recombined() &&
596  !orientedFaces[5].recombined()) {
597  gr->prisms.push_back(CREATE_PRISM_1);
598  gr->prisms.push_back(CREATE_PRISM_2);
599  }
600  else if(!orientedFaces[0].recombined() &&
601  !orientedFaces[1].recombined() &&
602  !orientedFaces[2].recombined() &&
603  !orientedFaces[3].recombined() &&
604  !orientedFaces[4].recombined() &&
605  !orientedFaces[5].recombined()) {
606  gr->tetrahedra.push_back(CREATE_SIM_1);
607  gr->tetrahedra.push_back(CREATE_SIM_2);
608  gr->tetrahedra.push_back(CREATE_SIM_3);
609  gr->tetrahedra.push_back(CREATE_SIM_4);
610  gr->tetrahedra.push_back(CREATE_SIM_5);
611  gr->tetrahedra.push_back(CREATE_SIM_6);
612  }
613  else {
614  Msg::Error("Wrong surface recombination in transfinite volume %d",
615  gr->tag());
616  return 0;
617  }
618  }
619  }
620  }
621  }
622  else if(faces.size() == 5) {
623  bool standard_algo = true;
624  for(int i = 0; i < 5; i++) {
625  if(orientedFaces[i].getSurface() &&
626  orientedFaces[i].getSurface()->meshAttributes.transfinite3) {
627  standard_algo = false;
628  break;
629  }
630  }
631  if(standard_algo) {
632  for(int j = 0; j < N_j - 1; j++) {
633  for(int k = 0; k < N_k - 1; k++) {
634 #if defined(HAVE_QUADTRI)
635  if(gr->meshAttributes.QuadTri) {
636  std::vector<MVertex *> verts;
637  verts.resize(6);
638  verts[0] = tab[0][j][k];
639  verts[1] = tab[1][j][k];
640  verts[2] = tab[1][j + 1][k];
641  verts[3] = tab[0][j][k + 1];
642  verts[4] = tab[1][j][k + 1];
643  verts[5] = tab[1][j + 1][k + 1];
644  if((!orientedFaces[0].recombined() && j == 0) ||
645  (!orientedFaces[2].recombined() && j == N_j - 2) ||
646  (!orientedFaces[4].recombined() && k == 0) ||
647  (!orientedFaces[5].recombined() && k == N_k - 2)) {
648  // make subdivided element
649  meshTransfElemWithInternalVertex(gr, verts, &boundary_diags);
650  }
651  else
652  gr->prisms.push_back(new MPrism(verts[0], verts[1], verts[2],
653  verts[3], verts[4], verts[5]));
654  // continue, skipping the rest which is for non-divided elements
655  continue;
656  }
657 #endif
658  if((orientedFaces[0].recombined() && orientedFaces[1].recombined() &&
659  orientedFaces[2].recombined() && orientedFaces[4].recombined() &&
660  orientedFaces[5].recombined()) ||
661  (orientedFaces[0].recombined() && orientedFaces[1].recombined() &&
662  orientedFaces[2].recombined() && !orientedFaces[4].recombined() &&
663  !orientedFaces[5].recombined())) {
664  gr->prisms.push_back(new MPrism(
665  tab[0][j][k], tab[1][j][k], tab[1][j + 1][k], tab[0][j][k + 1],
666  tab[1][j][k + 1], tab[1][j + 1][k + 1]));
667  }
668  else if(!orientedFaces[0].recombined() &&
669  !orientedFaces[1].recombined() &&
670  !orientedFaces[2].recombined() &&
671  !orientedFaces[4].recombined() &&
672  !orientedFaces[5].recombined()) {
673  gr->tetrahedra.push_back(new MTetrahedron(
674  tab[0][j][k], tab[1][j][k], tab[1][j + 1][k], tab[0][j][k + 1]));
675  gr->tetrahedra.push_back(
676  new MTetrahedron(tab[1][j][k], tab[1][j + 1][k], tab[0][j][k + 1],
677  tab[1][j][k + 1]));
678  gr->tetrahedra.push_back(
679  new MTetrahedron(tab[0][j][k + 1], tab[1][j + 1][k + 1],
680  tab[1][j][k + 1], tab[1][j + 1][k]));
681  }
682  else {
683  Msg::Error("Wrong surface recombination in transfinite volume %d",
684  gr->tag());
685  return 0;
686  }
687  }
688  }
689  for(int i = 1; i < N_i - 1; i++) {
690  for(int j = 0; j < N_j - 1; j++) {
691  for(int k = 0; k < N_k - 1; k++) {
692 #if defined(HAVE_QUADTRI)
693  if(gr->meshAttributes.QuadTri) {
694  // create vertex array
695  std::vector<MVertex *> verts;
696  verts.resize(8);
697  verts[0] = tab[i][j][k];
698  verts[1] = tab[i + 1][j][k];
699  verts[2] = tab[i + 1][j + 1][k];
700  verts[3] = tab[i][j + 1][k];
701  verts[4] = tab[i][j][k + 1];
702  verts[5] = tab[i + 1][j][k + 1];
703  verts[6] = tab[i + 1][j + 1][k + 1];
704  verts[7] = tab[i][j + 1][k + 1];
705  if((!orientedFaces[1].recombined() && i == N_i - 2) ||
706  (!orientedFaces[0].recombined() && j == 0) ||
707  (!orientedFaces[2].recombined() && j == N_j - 2) ||
708  (!orientedFaces[4].recombined() && k == 0) ||
709  (!orientedFaces[5].recombined() && k == N_k - 2)) {
710  // make subdivided element
711  meshTransfElemWithInternalVertex(gr, verts, &boundary_diags);
712  }
713  else
714  gr->hexahedra.push_back(
715  new MHexahedron(verts[0], verts[1], verts[2], verts[3],
716  verts[4], verts[5], verts[6], verts[7]));
717  // continue, skipping the rest which is for non-divided elements
718  continue;
719  }
720 #endif
721  if(orientedFaces[0].recombined() && orientedFaces[1].recombined() &&
722  orientedFaces[2].recombined() && orientedFaces[4].recombined() &&
723  orientedFaces[5].recombined()) {
724  gr->hexahedra.push_back(CREATE_HEX);
725  }
726  else if(orientedFaces[0].recombined() &&
727  orientedFaces[1].recombined() &&
728  orientedFaces[2].recombined() &&
729  !orientedFaces[4].recombined() &&
730  !orientedFaces[5].recombined()) {
731  gr->prisms.push_back(CREATE_PRISM_1);
732  gr->prisms.push_back(CREATE_PRISM_2);
733  }
734  else if(!orientedFaces[0].recombined() &&
735  !orientedFaces[1].recombined() &&
736  !orientedFaces[2].recombined() &&
737  !orientedFaces[4].recombined() &&
738  !orientedFaces[5].recombined()) {
739  gr->tetrahedra.push_back(CREATE_SIM_1);
740  gr->tetrahedra.push_back(CREATE_SIM_2);
741  gr->tetrahedra.push_back(CREATE_SIM_3);
742  gr->tetrahedra.push_back(CREATE_SIM_4);
743  gr->tetrahedra.push_back(CREATE_SIM_5);
744  gr->tetrahedra.push_back(CREATE_SIM_6);
745  }
746  else {
747  Msg::Error("Wrong surface recombination in transfinite volume %d",
748  gr->tag());
749  return 0;
750  }
751  }
752  }
753  }
754  }
755  else { // if surfaces meshed with specific algo for 3-sided surfaces
756  for(int i = 0; i < N_i - 1; i++) {
757  int j = i;
758  for(int k = 0; k < N_k - 1; k++) {
759  if(orientedFaces[0].recombined() && orientedFaces[1].recombined() &&
760  orientedFaces[2].recombined()) {
761  gr->prisms.push_back(CREATE_PRISM_4);
762  }
763  else if(!orientedFaces[0].recombined() &&
764  !orientedFaces[1].recombined() &&
765  !orientedFaces[2].recombined() &&
766  !orientedFaces[4].recombined() &&
767  !orientedFaces[5].recombined()) {
768  gr->tetrahedra.push_back(CREATE_SIM_10);
769  gr->tetrahedra.push_back(CREATE_SIM_11);
770  gr->tetrahedra.push_back(CREATE_SIM_12);
771  }
772  else {
773  Msg::Error("Wrong surface recombination in transfinite volume %d",
774  gr->tag());
775  return 0;
776  }
777  }
778  }
779  for(int i = 1; i < N_i - 1; i++) {
780  for(int j = 0; j < i; j++) {
781  for(int k = 0; k < N_k - 1; k++) {
782  if(orientedFaces[0].recombined() && orientedFaces[1].recombined() &&
783  orientedFaces[2].recombined()) {
784  gr->prisms.push_back(CREATE_PRISM_3);
785  gr->prisms.push_back(CREATE_PRISM_4);
786  }
787  else if(!orientedFaces[0].recombined() &&
788  !orientedFaces[1].recombined() &&
789  !orientedFaces[2].recombined() &&
790  !orientedFaces[4].recombined() &&
791  !orientedFaces[5].recombined()) {
792  gr->tetrahedra.push_back(CREATE_SIM_7);
793  gr->tetrahedra.push_back(CREATE_SIM_8);
794  gr->tetrahedra.push_back(CREATE_SIM_9);
795  gr->tetrahedra.push_back(CREATE_SIM_10);
796  gr->tetrahedra.push_back(CREATE_SIM_11);
797  gr->tetrahedra.push_back(CREATE_SIM_12);
798  }
799  else {
800  Msg::Error("Wrong surface recombination in transfinite volume %d",
801  gr->tag());
802  return 0;
803  }
804  }
805  }
806  }
807  }
808  }
809 
810  return 1;
811 }
GOrientedTransfiniteFace::GOrientedTransfiniteFace
GOrientedTransfiniteFace(GFace *gf, std::vector< MVertex * > &corners)
Definition: meshGRegionTransfinite.cpp:188
contextMeshOptions::recombineAll
int recombineAll
Definition: Context.h:30
CREATE_PRISM_2
#define CREATE_PRISM_2
Definition: meshGRegionTransfinite.cpp:67
GOrientedTransfiniteFace::getSurface
GFace * getSurface()
Definition: meshGRegionTransfinite.cpp:256
GFace::edges
virtual std::vector< GEdge * > const & edges() const
Definition: GFace.h:121
CREATE_SIM_9
#define CREATE_SIM_9
Definition: meshGRegionTransfinite.cpp:112
CREATE_SIM_8
#define CREATE_SIM_8
Definition: meshGRegionTransfinite.cpp:108
GRegion::method
char method
Definition: GRegion.h:146
GFace::recombine
int recombine
Definition: GFace.h:344
GFace.h
MTetrahedron
Definition: MTetrahedron.h:34
GOrientedTransfiniteFace::index
int index() const
Definition: meshGRegionTransfinite.cpp:258
CREATE_PRISM_3
#define CREATE_PRISM_3
Definition: meshGRegionTransfinite.cpp:72
GFace
Definition: GFace.h:33
Msg::Info
static void Info(const char *fmt,...)
Definition: GmshMessage.cpp:587
GOrientedTransfiniteFace::recombined
int recombined() const
Definition: meshGRegionTransfinite.cpp:260
CREATE_SIM_4
#define CREATE_SIM_4
Definition: meshGRegionTransfinite.cpp:92
Msg::Debug
static void Debug(const char *fmt,...)
Definition: GmshMessage.cpp:752
c
static double c(int i, int j, fullMatrix< double > &CA, const std::vector< SPoint3 > &P, const std::vector< SPoint3 > &Q)
Definition: discreteFrechetDistance.cpp:15
MESH_TRANSFINITE
#define MESH_TRANSFINITE
Definition: GmshDefines.h:259
MVertex
Definition: MVertex.h:24
Msg::Error
static void Error(const char *fmt,...)
Definition: GmshMessage.cpp:482
MVertex::z
double z() const
Definition: MVertex.h:62
LegendrePolynomials::f
void f(int n, double u, double *val)
Definition: orthogonalBasis.cpp:77
GFace::meshAttributes
struct GFace::@18 meshAttributes
GRegion::transfinite_vertices
std::vector< std::vector< std::vector< MVertex * > > > transfinite_vertices
Definition: GRegion.h:161
MPrism
Definition: MPrism.h:34
GmshMessage.h
edges
static int edges[6][2]
Definition: meshGRegionLocalMeshMod.cpp:23
GRegion::corners
std::vector< GVertex * > corners
Definition: GRegion.h:150
GOrientedTransfiniteFace::_index
int _index
Definition: meshGRegionTransfinite.cpp:180
CREATE_SIM_1
#define CREATE_SIM_1
Definition: meshGRegionTransfinite.cpp:80
GOrientedTransfiniteFace::getNumV
int getNumV()
Definition: meshGRegionTransfinite.cpp:267
CTX::instance
static CTX * instance()
Definition: Context.cpp:122
GRegion::faces
virtual std::vector< GFace * > faces() const
Definition: GRegion.h:64
GRegion::hexahedra
std::vector< MHexahedron * > hexahedra
Definition: GRegion.h:164
CREATE_SIM_12
#define CREATE_SIM_12
Definition: meshGRegionTransfinite.cpp:124
GRegion.h
MeshTransfiniteVolume
int MeshTransfiniteVolume(GRegion *gr)
Definition: meshGRegionTransfinite.cpp:342
GEntity::mesh_vertices
std::vector< MVertex * > mesh_vertices
Definition: GEntity.h:56
CREATE_PRISM_4
#define CREATE_PRISM_4
Definition: meshGRegionTransfinite.cpp:76
MHexahedron.h
CREATE_SIM_3
#define CREATE_SIM_3
Definition: meshGRegionTransfinite.cpp:88
GOrientedTransfiniteFace::getVertex
MVertex * getVertex(int i, int j)
Definition: meshGRegionTransfinite.cpp:270
MVertex.h
findTransfiniteCorners
void findTransfiniteCorners(GRegion *gr, std::vector< MVertex * > &corners)
Definition: meshGRegionTransfinite.cpp:295
GOrientedTransfiniteFace::_hh
int _hh
Definition: meshGRegionTransfinite.cpp:179
GOrientedTransfiniteFace::_permutation
int _permutation
Definition: meshGRegionTransfinite.cpp:180
MHexahedron
Definition: MHexahedron.h:28
GRegion::edges
virtual std::vector< GEdge * > const & edges() const
Definition: GRegion.cpp:459
CTX::mesh
contextMeshOptions mesh
Definition: Context.h:313
GEntity::tag
int tag() const
Definition: GEntity.h:280
CREATE_SIM_6
#define CREATE_SIM_6
Definition: meshGRegionTransfinite.cpp:100
CREATE_HEX
#define CREATE_HEX
Definition: meshGRegionTransfinite.cpp:58
GRegion
Definition: GRegion.h:28
GRegion::prisms
std::vector< MPrism * > prisms
Definition: GRegion.h:165
GRegion::QuadTri
int QuadTri
Definition: GRegion.h:152
CREATE_SIM_2
#define CREATE_SIM_2
Definition: meshGRegionTransfinite.cpp:84
GOrientedTransfiniteFace::_ll
int _ll
Definition: meshGRegionTransfinite.cpp:179
CREATE_SIM_11
#define CREATE_SIM_11
Definition: meshGRegionTransfinite.cpp:120
meshGFace.h
GOrientedTransfiniteFace::getNumU
int getNumU()
Definition: meshGRegionTransfinite.cpp:266
Context.h
MTetrahedron.h
CREATE_PRISM_1
#define CREATE_PRISM_1
Definition: meshGRegionTransfinite.cpp:63
transfiniteHex
static double transfiniteHex(double f1, double f2, double f3, double f4, double f5, double f6, double c1, double c2, double c3, double c4, double c5, double c6, double c7, double c8, double c9, double c10, double c11, double c12, double s1, double s2, double s3, double s4, double s5, double s6, double s7, double s8, double u, double v, double w)
Definition: meshGRegionTransfinite.cpp:128
z
const double z
Definition: GaussQuadratureQuad.cpp:56
GRegion::meshAttributes
struct GRegion::@20 meshAttributes
CREATE_SIM_5
#define CREATE_SIM_5
Definition: meshGRegionTransfinite.cpp:96
GOrientedTransfiniteFace::GOrientedTransfiniteFace
GOrientedTransfiniteFace()
Definition: meshGRegionTransfinite.cpp:184
MPrism.h
GOrientedTransfiniteFace
Definition: meshGRegionTransfinite.cpp:176
GFace::transfinite_vertices
std::vector< std::vector< MVertex * > > transfinite_vertices
Definition: GFace.h:420
MVertex::y
double y() const
Definition: MVertex.h:61
CREATE_SIM_7
#define CREATE_SIM_7
Definition: meshGRegionTransfinite.cpp:104
GOrientedTransfiniteFace::_gf
GFace * _gf
Definition: meshGRegionTransfinite.cpp:178
MVertex::distance
double distance(MVertex *const v)
Definition: MVertex.h:105
GOrientedTransfiniteFace::_list
std::vector< MVertex * > _list
Definition: meshGRegionTransfinite.cpp:181
c1
const double c1
Definition: GaussQuadratureHex.cpp:16
GRegion::tetrahedra
std::vector< MTetrahedron * > tetrahedra
Definition: GRegion.h:163
MVertex::x
double x() const
Definition: MVertex.h:60
CREATE_SIM_10
#define CREATE_SIM_10
Definition: meshGRegionTransfinite.cpp:116
faces
static int faces[4][3]
Definition: meshGRegionDelaunayInsertion.cpp:165
GFace::corners
std::vector< GVertex * > corners
Definition: GFace.h:350