1 // Gmsh - Copyright (C) 1997-2021 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 #include <map>
7 #include <queue>
8 #include <array>
9 #include "meshGFace.h"
10 #include "meshGEdge.h"
11 #include "GVertex.h"
12 #include "GEdge.h"
13 #include "GFace.h"
14 #include "MVertex.h"
15 #include "MTriangle.h"
16 #include "MQuadrangle.h"
17 #include "Context.h"
18 #include "GmshMessage.h"
19 #include "Numeric.h"
20 
21 /*
22    s4 +-----c3-----+ s3
23       |            |
24       |            |
25      c4            c2
26       |            |
27       |            |
28    s1 +-----c1-----+ s2
29 */
30 
31 // f(u,v) = (1-u) c4(v) + u c2(v) + (1-v) c1(u) + v c3(u)
32 //          - [ (1-u)(1-v) s1 + u(1-v) s2 + uv s3 + (1-u)v s4 ]
33 #define TRAN_QUA(c1, c2, c3, c4, s1, s2, s3, s4, u, v)                         \
34   (1. - u) * c4 + u * c2 + (1. - v) * c1 + v * c3 -                            \
35     ((1. - u) * (1. - v) * s1 + u * (1. - v) * s2 + u * v * s3 +               \
36      (1. - u) * v * s4)
37 
38 // s1=s4=c4
39 // f(u,v) = u c2 (v) + (1-v) c1(u) + v c3(u) - u(1-v) s2 - uv s3
40 #define TRAN_TRI(c1, c2, c3, s1, s2, s3, u, v)                                 \
41   u * c2 + (1. - v) * c1 + v * c3 - (u * (1. - v) * s2 + u * v * s3)
42 
findTransfiniteCorners(GFace * gf,std::vector<MVertex * > & corners)43 void findTransfiniteCorners(GFace *gf, std::vector<MVertex *> &corners)
44 {
45   if(gf->meshAttributes.corners.size()) {
46     // corners have been specified explicitly
47     for(std::size_t i = 0; i < gf->meshAttributes.corners.size(); i++)
48       corners.push_back(gf->meshAttributes.corners[i]->mesh_vertices[0]);
49   }
50   else {
51     // try to find the corners automatically
52     std::vector<GEdge *> fedges = gf->edges();
53     GEdgeLoop el(fedges);
54     for(auto it = el.begin(); it != el.end(); it++)
55       corners.push_back(it->getBeginVertex()->mesh_vertices[0]);
56 
57     // try reaaally hard for 3-sided faces
58     if(corners.size() == 3) {
59       GEdge *first = nullptr, *last = nullptr;
60       for(auto it = fedges.begin(); it != fedges.end(); it++) {
61         if(((*it)->getBeginVertex()->mesh_vertices[0] == corners[0] &&
62             (*it)->getEndVertex()->mesh_vertices[0] == corners[1]) ||
63            ((*it)->getBeginVertex()->mesh_vertices[0] == corners[1] &&
64             (*it)->getEndVertex()->mesh_vertices[0] == corners[0])) {
65           first = *it;
66         }
67         if(((*it)->getBeginVertex()->mesh_vertices[0] == corners[2] &&
68             (*it)->getEndVertex()->mesh_vertices[0] == corners[0]) ||
69            ((*it)->getBeginVertex()->mesh_vertices[0] == corners[0] &&
70             (*it)->getEndVertex()->mesh_vertices[0] == corners[2])) {
71           last = *it;
72         }
73       }
74       if(first && last) {
75         if(first->mesh_vertices.size() != last->mesh_vertices.size()) {
76           std::vector<MVertex *> corners2(3);
77           corners2[0] = corners[1];
78           corners2[1] = corners[2];
79           corners2[2] = corners[0];
80           corners = corners2;
81         }
82       }
83     }
84   }
85 }
86 
computeEdgeLoops(const GFace * gf,std::vector<MVertex * > & all_mvertices,std::vector<int> & indices)87 static void computeEdgeLoops(const GFace *gf,
88                              std::vector<MVertex *> &all_mvertices,
89                              std::vector<int> &indices)
90 {
91   std::vector<GEdge *> const &e = gf->edges();
92   std::vector<int> const &o = gf->orientations();
93 
94   std::vector<GEdge *> edges;
95   std::vector<int> ori;
96   for(std::size_t i = 0; i < e.size(); i++) {
97     if(!e[i]->degenerate(0)) { // skip degenerate curves
98       edges.push_back(e[i]);
99       ori.push_back(i < o.size() ? o[i] : 1);
100     }
101   }
102 
103   auto it = edges.begin();
104   auto ito = ori.begin();
105   indices.push_back(0);
106   GVertex *start =
107     ((*ito) == 1) ? (*it)->getBeginVertex() : (*it)->getEndVertex();
108   GVertex *v_end =
109     ((*ito) != 1) ? (*it)->getBeginVertex() : (*it)->getEndVertex();
110 
111   all_mvertices.push_back(start->mesh_vertices[0]);
112   if(*ito == 1)
113     for(std::size_t i = 0; i < (*it)->mesh_vertices.size(); i++)
114       all_mvertices.push_back((*it)->mesh_vertices[i]);
115   else
116     for(int i = (*it)->mesh_vertices.size() - 1; i >= 0; i--)
117       all_mvertices.push_back((*it)->mesh_vertices[i]);
118 
119   GVertex *v_start = start;
120   while(1) {
121     ++it;
122     ++ito;
123     if(v_end == start) {
124       indices.push_back(all_mvertices.size());
125       if(it == edges.end()) break;
126       start = ((*ito) == 1) ? (*it)->getBeginVertex() : (*it)->getEndVertex();
127       v_end = ((*ito) != 1) ? (*it)->getBeginVertex() : (*it)->getEndVertex();
128       v_start = start;
129     }
130     else {
131       if(it == edges.end()) {
132         Msg::Error("Something wrong in edge loop computation");
133         return;
134       }
135       v_start = ((*ito) == 1) ? (*it)->getBeginVertex() : (*it)->getEndVertex();
136       if(v_start != v_end) {
137         Msg::Error("Something wrong in edge loop computation");
138         return;
139       }
140       v_end = ((*ito) != 1) ? (*it)->getBeginVertex() : (*it)->getEndVertex();
141     }
142     all_mvertices.push_back(v_start->mesh_vertices[0]);
143     if(*ito == 1)
144       for(std::size_t i = 0; i < (*it)->mesh_vertices.size(); i++)
145         all_mvertices.push_back((*it)->mesh_vertices[i]);
146     else
147       for(int i = (*it)->mesh_vertices.size() - 1; i >= 0; i--)
148         all_mvertices.push_back((*it)->mesh_vertices[i]);
149   }
150 }
151 
maxDistParam(const std::vector<double> & U,const std::vector<double> & V)152 double maxDistParam(const std::vector<double> &U, const std::vector<double> &V)
153 {
154   if(U.size() < 2 || (U.size() != V.size())) return 1e22;
155   double d = sqrt(std::pow(U.back() - U.front(), 2) +
156                   std::pow(V.back() - V.front(), 2));
157   for(std::size_t i = 1; i < U.size(); i++)
158     d = std::max(d,  sqrt(std::pow(U[i] - U[i - 1], 2) +
159                           std::pow(V[i] - V[i - 1], 2)));
160   return d;
161 }
162 
MeshTransfiniteSurface(GFace * gf)163 int MeshTransfiniteSurface(GFace *gf)
164 {
165   if(gf->meshAttributes.method != MESH_TRANSFINITE) return 0;
166 
167   Msg::Info("Meshing surface %d (Transfinite)", gf->tag());
168 
169   // make sure that all bounding edges have begin/end points: everything in here
170   // depends on it
171   const std::vector<GEdge *> &edges = gf->edges();
172   for(auto it = edges.begin(); it != edges.end(); it++) {
173     if(!(*it)->getBeginVertex() || !(*it)->getEndVertex()) {
174       Msg::Error("Transfinite algorithm cannot be applied with curve %d which "
175                  "has no begin or end point",
176                  (*it)->tag());
177       return 0;
178     }
179   }
180 
181   std::vector<MVertex *> corners;
182   findTransfiniteCorners(gf, corners);
183   if(corners.size() != 3 && corners.size() != 4) {
184     Msg::Error("Surface %d is transfinite but has %d corner%s", gf->tag(),
185                corners.size(), corners.size() > 1 ? "s" : "");
186     return 0;
187   }
188 
189   std::vector<MVertex *> d_vertices;
190   std::vector<int> indices;
191   computeEdgeLoops(gf, d_vertices, indices);
192 
193   if(indices.size() != 2) {
194     int nh = indices.size() - 2;
195     Msg::Error("Surface %d is transfinite but has %d hole%s", gf->tag(), nh,
196                nh > 1 ? "s" : "");
197     return 0;
198   }
199 
200   // create a list of all boundary vertices, starting at the first
201   // transfinite corner
202   std::vector<MVertex *> m_vertices;
203   std::size_t I;
204   for(I = 0; I < d_vertices.size(); I++)
205     if(d_vertices[I] == corners[0]) break;
206   for(std::size_t j = 0; j < d_vertices.size(); j++)
207     m_vertices.push_back(d_vertices[(I + j) % d_vertices.size()]);
208 
209   // make the ordering of the list consistent with the ordering of the
210   // first two corners (if the second found corner is not the second
211   // corner, just reverse the list)
212   bool reverse = false;
213   for(std::size_t i = 1; i < m_vertices.size(); i++) {
214     MVertex *v = m_vertices[i];
215     if(v == corners[1] || v == corners[2] ||
216        (corners.size() == 4 && v == corners[3])) {
217       if(v != corners[1]) reverse = true;
218       break;
219     }
220   }
221   if(reverse) {
222     std::vector<MVertex *> tmp;
223     tmp.push_back(m_vertices[0]);
224     for(int i = m_vertices.size() - 1; i > 0; i--) tmp.push_back(m_vertices[i]);
225     m_vertices = tmp;
226   }
227 
228   // get the indices of the interpolation corners as well as the u,v
229   // coordinates of all the boundary vertices
230   int iCorner = 0, N[4] = {0, 0, 0, 0};
231   std::vector<double> U, V, U2, V2;
232   for(std::size_t i = 0; i < m_vertices.size(); i++) {
233     MVertex *v = m_vertices[i];
234     if(v == corners[0] || v == corners[1] || v == corners[2] ||
235        (corners.size() == 4 && v == corners[3])) {
236       if(iCorner > 3) {
237         Msg::Error("Surface %d transfinite parameters are incoherent",
238                    gf->tag());
239         return 0;
240       }
241       N[iCorner++] = i;
242     }
243     SPoint2 param;
244     bool onSeam = !reparamMeshVertexOnFace(v, gf, param);
245     U.push_back(param[0]);
246     V.push_back(param[1]);
247     // if v is on a periodic seam, get the other parametric coordinates
248     if(onSeam) reparamMeshVertexOnFace(v, gf, param, true, true, -1);
249     U2.push_back(param[0]);
250     V2.push_back(param[1]);
251   }
252 
253   // quick fix for surfaces with a single seam: choose the trajectory in the
254   // parametric plane that leads to the smallest maximum distance between two
255   // nodes; this should be generalized if we want to handle cases with more than
256   // one seam
257   if(maxDistParam(U2, V2) < maxDistParam(U, V)) {
258     Msg::Debug("Choosing alternate parametrization on surface %d", gf->tag());
259     U = U2;
260     V = V2;
261   }
262 
263   int N1 = N[0], N2 = N[1], N3 = N[2], N4 = N[3];
264   int L = N2 - N1, H = N3 - N2;
265   if(corners.size() == 4) {
266     int Lb = N4 - N3, Hb = m_vertices.size() - N4;
267     if(Lb != L || Hb != H) {
268       Msg::Error("Surface %d cannot be meshed using the transfinite algo "
269                  "(divisions %d != %d or %d != %d)",
270                  gf->tag(), Lb, L, Hb, H);
271       return 0;
272     }
273   }
274   else {
275     int Lb = m_vertices.size() - N3;
276     if(Lb != L) {
277       Msg::Error("Surface %d cannot be meshed using the transfinite algo "
278                  "(divisions %d != %d)",
279                  gf->tag(), L, Lb);
280       return 0;
281     }
282   }
283 
284   /*
285       2L+H +------------+ L+H
286            |            |
287            |            |
288            |            |
289            |            |
290      2L+2H +------------+
291            0            L
292   */
293 
294   // use the average of the node spacing on the two opposing sides, so that we
295   // generate the same u, v coordinates whatever the ordering of the sides
296   bool symmetric = true;
297 
298   std::vector<double> lengths_i;
299   lengths_i.reserve(L);
300   lengths_i.push_back(0.0);
301   double L_i = 0.0;
302   for(int i = 0; i < L; i++) {
303     MVertex *v1 = m_vertices[i];
304     MVertex *v2 = m_vertices[i + 1];
305     double d1 = v1->distance(v2);
306     if(symmetric) {
307       MVertex *v3 = m_vertices[(corners.size() == 3 && !i) ? 0 : (2 * L + H - i)];
308       MVertex *v4 = m_vertices[2 * L + H - i - 1];
309       double d2 = v3->distance(v4);
310       L_i += 0.5 * (d1 + d2);
311     }
312     else {
313       L_i += d1;
314     }
315     lengths_i.push_back(L_i);
316   }
317 
318   std::vector<double> lengths_j;
319   lengths_j.reserve(H);
320   lengths_j.push_back(0.0);
321   double L_j = 0.0;
322   for(int i = 0; i < H; i++) {
323     MVertex *v1 = m_vertices[L + i];
324     MVertex *v2 = m_vertices[L + i + 1];
325     double d1 = v1->distance(v2);
326     if(symmetric && corners.size() == 4) {
327       MVertex *v3 = m_vertices[(!i) ? 0 : (2 * L + 2 * H - i)];
328       MVertex *v4 = m_vertices[2 * L + 2 * H - i - 1];
329       double d2 = v3->distance(v4);
330       L_j += 0.5 * (d1 + d2);
331     }
332     else{
333       L_j += d1;
334     }
335     lengths_j.push_back(L_j);
336   }
337 
338   std::vector<std::vector<MVertex *> > &tab(gf->transfinite_vertices);
339   tab.resize(L + 1);
340   for(int i = 0; i <= L; i++) tab[i].resize(H + 1);
341 
342   if(corners.size() == 4) {
343     tab[0][0] = m_vertices[0];
344     tab[L][0] = m_vertices[L];
345     tab[L][H] = m_vertices[L + H];
346     tab[0][H] = m_vertices[2 * L + H];
347     for(int i = 1; i < L; i++) {
348       tab[i][0] = m_vertices[i];
349       tab[i][H] = m_vertices[2 * L + H - i];
350     }
351     for(int i = 1; i < H; i++) {
352       tab[L][i] = m_vertices[L + i];
353       tab[0][i] = m_vertices[2 * L + 2 * H - i];
354     }
355   }
356   else {
357     tab[0][0] = m_vertices[0];
358     tab[L][0] = m_vertices[L];
359     tab[L][H] = m_vertices[L + H];
360     // degenerated, only necessary for transfinite volume algo
361     tab[0][H] = m_vertices[0];
362     for(int i = 1; i < L; i++) {
363       tab[i][0] = m_vertices[i];
364       tab[i][H] = m_vertices[2 * L + H - i];
365     }
366     for(int i = 1; i < H; i++) {
367       tab[L][i] = m_vertices[L + i];
368       // degenerated, only necessary for transfinite volume algo
369       tab[0][i] = m_vertices[0];
370     }
371   }
372 
373   double UC1 = U[N1], UC2 = U[N2], UC3 = U[N3];
374   double VC1 = V[N1], VC2 = V[N2], VC3 = V[N3];
375 
376   // create points using transfinite interpolation
377   if(corners.size() == 4) {
378     double UC4 = U[N4];
379     double VC4 = V[N4];
380     for(int i = 1; i < L; i++) {
381       double u = lengths_i[i] / L_i;
382       for(int j = 1; j < H; j++) {
383         double v = lengths_j[j] / L_j;
384         int iP1 = N1 + i;
385         int iP2 = N2 + j;
386         int iP3 = N4 - i;
387         int iP4 = (N4 + (N3 - N2) - j) % m_vertices.size();
388         double Up =
389           TRAN_QUA(U[iP1], U[iP2], U[iP3], U[iP4], UC1, UC2, UC3, UC4, u, v);
390         double Vp =
391           TRAN_QUA(V[iP1], V[iP2], V[iP3], V[iP4], VC1, VC2, VC3, VC4, u, v);
392         GPoint gp = gf->point(SPoint2(Up, Vp));
393         MFaceVertex *newv = new MFaceVertex(gp.x(), gp.y(), gp.z(), gf, Up, Vp);
394         gf->mesh_vertices.push_back(newv);
395         tab[i][j] = newv;
396       }
397     }
398   }
399   else {
400     for(int i = 1; i < L; i++) {
401       double u = lengths_i[i] / L_i;
402       for(int j = 1; j < H; j++) {
403         double v = lengths_j[j] / L_j;
404         int iP1 = N1 + i;
405         int iP2 = N2 + j;
406         int iP3 = ((N3 + N2) - i) % m_vertices.size();
407         double Up, Vp;
408         if(gf->geomType() != GEntity::RuledSurface) {
409           Up = TRAN_TRI(U[iP1], U[iP2], U[iP3], UC1, UC2, UC3, u, v);
410           Vp = TRAN_TRI(V[iP1], V[iP2], V[iP3], VC1, VC2, VC3, u, v);
411         }
412         else {
413           // FIXME: to get nice meshes we would need to make the u,v
414           // coords match with the (degenerate) coordinates of the
415           // underlying ruled surface; so instead we just interpolate
416           // in real space
417           double xp = TRAN_TRI(m_vertices[iP1]->x(), m_vertices[iP2]->x(),
418                                m_vertices[iP3]->x(), m_vertices[N1]->x(),
419                                m_vertices[N2]->x(), m_vertices[N3]->x(), u, v);
420           double yp = TRAN_TRI(m_vertices[iP1]->y(), m_vertices[iP2]->y(),
421                                m_vertices[iP3]->y(), m_vertices[N1]->y(),
422                                m_vertices[N2]->y(), m_vertices[N3]->y(), u, v);
423           double zp = TRAN_TRI(m_vertices[iP1]->z(), m_vertices[iP2]->z(),
424                                m_vertices[iP3]->z(), m_vertices[N1]->z(),
425                                m_vertices[N2]->z(), m_vertices[N3]->z(), u, v);
426           // xp,yp,zp can be off the surface so we cannot use parFromPoint
427           gf->XYZtoUV(xp, yp, zp, Up, Vp, 1.0, false);
428         }
429         GPoint gp = gf->point(SPoint2(Up, Vp));
430         MFaceVertex *newv = new MFaceVertex(gp.x(), gp.y(), gp.z(), gf, Up, Vp);
431         gf->mesh_vertices.push_back(newv);
432         tab[i][j] = newv;
433       }
434     }
435   }
436 
437   // should we apply the elliptic smoother?
438   int numSmooth = 0;
439   if(gf->meshAttributes.transfiniteSmoothing < 0 &&
440      CTX::instance()->mesh.nbSmoothing > 1)
441     numSmooth = CTX::instance()->mesh.nbSmoothing;
442   else if(gf->meshAttributes.transfiniteSmoothing > 0)
443     numSmooth = gf->meshAttributes.transfiniteSmoothing;
444 
445   if(corners.size() == 4 && numSmooth) {
446     std::vector<std::vector<double> > u(L + 1), v(L + 1);
447     for(int i = 0; i <= L; i++) {
448       u[i].resize(H + 1);
449       v[i].resize(H + 1);
450     }
451     for(int i = 0; i <= L; i++) {
452       for(int j = 0; j <= H; j++) {
453         int iP1 = N1 + i;
454         int iP2 = N2 + j;
455         int iP3 = N4 - i;
456         int iP4 = (N4 + (N3 - N2) - j) % m_vertices.size();
457         if(j == 0) {
458           u[i][j] = U[iP1];
459           v[i][j] = V[iP1];
460         }
461         else if(i == L) {
462           u[i][j] = U[iP2];
463           v[i][j] = V[iP2];
464         }
465         else if(j == H) {
466           u[i][j] = U[iP3];
467           v[i][j] = V[iP3];
468         }
469         else if(i == 0) {
470           u[i][j] = U[iP4];
471           v[i][j] = V[iP4];
472         }
473         else {
474           tab[i][j]->getParameter(0, u[i][j]);
475           tab[i][j]->getParameter(1, v[i][j]);
476         }
477       }
478     }
479     for(int IT = 0; IT < numSmooth; IT++) {
480       for(int i = 1; i < L; i++) {
481         for(int j = 1; j < H; j++) {
482           double alpha = 0.25 * (std::pow(u[i][j + 1] - u[i][j - 1], 2) +
483                                  std::pow(v[i][j + 1] - v[i][j - 1], 2));
484           double gamma = 0.25 * (std::pow(u[i + 1][j] - u[i - 1][j], 2) +
485                                  std::pow(v[i + 1][j] - v[i - 1][j], 2));
486           double beta =
487             0.0625 *
488             ((u[i + 1][j] - u[i - 1][j]) * (u[i][j + 1] - u[i][j - 1]) +
489              (v[i + 1][j] - v[i - 1][j]) * (v[i][j + 1] - v[i][j - 1]));
490           u[i][j] = 0.5 *
491                     (alpha * (u[i + 1][j] + u[i - 1][j]) +
492                      gamma * (u[i][j + 1] + u[i][j - 1]) -
493                      2. * beta *
494                        (u[i + 1][j + 1] - u[i - 1][j + 1] - u[i + 1][j - 1] +
495                         u[i - 1][j - 1])) /
496                     (alpha + gamma);
497           v[i][j] = 0.5 *
498                     (alpha * (v[i + 1][j] + v[i - 1][j]) +
499                      gamma * (v[i][j + 1] + v[i][j - 1]) -
500                      2. * beta *
501                        (v[i + 1][j + 1] - v[i - 1][j + 1] - v[i + 1][j - 1] +
502                         v[i - 1][j - 1])) /
503                     (alpha + gamma);
504         }
505       }
506     }
507     for(int i = 1; i < L; i++) {
508       for(int j = 1; j < H; j++) {
509         GPoint p = gf->point(SPoint2(u[i][j], v[i][j]));
510         tab[i][j]->x() = p.x();
511         tab[i][j]->y() = p.y();
512         tab[i][j]->z() = p.z();
513         tab[i][j]->setParameter(0, u[i][j]);
514         tab[i][j]->setParameter(1, v[i][j]);
515       }
516     }
517   }
518 
519   // create elements
520   if(corners.size() == 4) {
521     for(int i = 0; i < L; i++) {
522       for(int j = 0; j < H; j++) {
523         MVertex *v1 = tab[i][j];
524         MVertex *v2 = tab[i + 1][j];
525         MVertex *v3 = tab[i + 1][j + 1];
526         MVertex *v4 = tab[i][j + 1];
527         if(CTX::instance()->mesh.recombineAll || gf->meshAttributes.recombine)
528           gf->quadrangles.push_back(new MQuadrangle(v1, v2, v3, v4));
529         else if(gf->meshAttributes.transfiniteArrangement == 1 ||
530                 (gf->meshAttributes.transfiniteArrangement == 2 &&
531                  ((i % 2 == 0 && j % 2 == 1) || (i % 2 == 1 && j % 2 == 0))) ||
532                 (gf->meshAttributes.transfiniteArrangement == -2 &&
533                  ((i % 2 == 0 && j % 2 == 0) || (i % 2 == 1 && j % 2 == 1)))) {
534           //        else if(rand() % 2 == 0){
535           gf->triangles.push_back(new MTriangle(v1, v2, v3));
536           gf->triangles.push_back(new MTriangle(v3, v4, v1));
537         }
538         else {
539           gf->triangles.push_back(new MTriangle(v1, v2, v4));
540           gf->triangles.push_back(new MTriangle(v4, v2, v3));
541         }
542       }
543     }
544   }
545   else {
546     for(int j = 0; j < H; j++) {
547       MVertex *v1 = tab[0][0];
548       MVertex *v2 = tab[1][j];
549       MVertex *v3 = tab[1][j + 1];
550       gf->triangles.push_back(new MTriangle(v1, v2, v3));
551     }
552     for(int i = 1; i < L; i++) {
553       for(int j = 0; j < H; j++) {
554         MVertex *v1 = tab[i][j];
555         MVertex *v2 = tab[i + 1][j];
556         MVertex *v3 = tab[i + 1][j + 1];
557         MVertex *v4 = tab[i][j + 1];
558         if(CTX::instance()->mesh.recombineAll || gf->meshAttributes.recombine)
559           gf->quadrangles.push_back(new MQuadrangle(v1, v2, v3, v4));
560         else if(gf->meshAttributes.transfiniteArrangement == 1 ||
561                 (gf->meshAttributes.transfiniteArrangement == 2 &&
562                  ((i % 2 == 0 && j % 2 == 1) || (i % 2 == 1 && j % 2 == 0))) ||
563                 (gf->meshAttributes.transfiniteArrangement == -2 &&
564                  ((i % 2 == 0 && j % 2 == 0) || (i % 2 == 1 && j % 2 == 1)))) {
565           gf->triangles.push_back(new MTriangle(v1, v2, v3));
566           gf->triangles.push_back(new MTriangle(v3, v4, v1));
567         }
568         else {
569           gf->triangles.push_back(new MTriangle(v1, v2, v4));
570           gf->triangles.push_back(new MTriangle(v4, v2, v3));
571         }
572       }
573     }
574   }
575 
576   gf->meshStatistics.status = GFace::DONE;
577   return 1;
578 }
579 
id_loop_from_closed_pairs(const std::vector<std::array<int,2>> & pairs,std::vector<int> & loop)580 bool id_loop_from_closed_pairs(const std::vector<std::array<int, 2> > &pairs,
581                                std::vector<int> &loop)
582 {
583   /* warning: does not verify if the loop is closed, eg [1,2], [2,3] will
584    * produce [1,2,3] */
585   loop.clear();
586   if(pairs.size() < 2) return false;
587   loop.reserve(pairs.size());
588   std::vector<bool> done(pairs.size(), false);
589   loop.resize(2);
590   loop[0] = pairs[0][0];
591   loop[1] = pairs[0][1];
592   done[0] = true;
593   while(loop.size() != pairs.size()) {
594     int cv = loop.back();
595     bool found = false;
596     for(std::size_t i = 0; i < pairs.size(); ++i)
597       if(!done[i]) {
598         const std::array<int, 2> &vs = pairs[i];
599         if(vs[0] == cv) {
600           loop.push_back(vs[1]);
601           done[i] = true;
602           found = true;
603         }
604         else if(vs[1] == cv) {
605           loop.push_back(vs[0]);
606           done[i] = true;
607           found = true;
608         }
609       }
610     if(!found) {
611       loop.clear();
612       return false;
613     }
614   }
615 
616   return true;
617 }
618 
faceIsValidQuad(GFace * gf,double angle_threshold_rad)619 bool faceIsValidQuad(GFace *gf, double angle_threshold_rad)
620 {
621   if(gf->edges().size() != 4) return false;
622 
623   std::vector<std::array<int, 2> > vpairs;
624   for(GEdge *ge : gf->edges()) {
625     vpairs.push_back(
626       {{ge->getBeginVertex()->tag(), ge->getEndVertex()->tag()}});
627   }
628   std::vector<int> loop;
629   bool ok = id_loop_from_closed_pairs(vpairs, loop);
630   if(!ok || loop.size() != vpairs.size()) return false;
631 
632   if(angle_threshold_rad > 0.) {
633     /* Check angle at each corner */
634     for(std::size_t i = 0; i < loop.size(); ++i) {
635       int v1 = loop[i];
636       int v2 = loop[(i + 1) % loop.size()];
637       int v3 = loop[(i + 2) % loop.size()];
638       SVector3 t21, t23;
639       for(GEdge *ge : gf->edges()) {
640         int e_v1 = ge->getBeginVertex()->tag();
641         int e_v2 = ge->getEndVertex()->tag();
642         if(e_v1 == v1 && e_v2 == v2) {
643           Range<double> range = ge->parBounds(0);
644           SVector3 t = ge->firstDer(range.high());
645           t.normalize();
646           t21 = -t;
647         }
648         else if(e_v1 == v2 && e_v2 == v1) {
649           Range<double> range = ge->parBounds(0);
650           SVector3 t = ge->firstDer(range.low());
651           t.normalize();
652           t21 = t;
653         }
654         else if(e_v1 == v2 && e_v2 == v3) {
655           Range<double> range = ge->parBounds(0);
656           SVector3 t = ge->firstDer(range.low());
657           t.normalize();
658           t23 = t;
659         }
660         else if(e_v1 == v3 && e_v2 == v2) {
661           Range<double> range = ge->parBounds(0);
662           SVector3 t = ge->firstDer(range.high());
663           t.normalize();
664           t23 = -t;
665         }
666       }
667       double agl = angle(t21, t23);
668       if(agl > angle_threshold_rad) {
669         Msg::Debug("quadrangular face %i rejected for automatic transfinite "
670                    "because corner angle = %.3f deg > threshold = %.3f deg\n",
671                    gf->tag(), 180. / M_PI * agl,
672                    180. / M_PI * angle_threshold_rad);
673         return false;
674       }
675     }
676   }
677 
678   return true;
679 }
680 
quad_face_opposite_edge(GFace * face,GEdge * edge)681 GEdge *quad_face_opposite_edge(GFace *face, GEdge *edge)
682 {
683   if(face->edges().size() != 4) return nullptr;
684   GEdge *op = nullptr;
685   int v1 = edge->getBeginVertex()->tag();
686   int v2 = edge->getEndVertex()->tag();
687   bool edgeInside = false;
688   for(GEdge *ge : face->edges()) {
689     if(ge == edge) {
690       edgeInside = true;
691       continue;
692     }
693     int cv1 = ge->getBeginVertex()->tag();
694     int cv2 = ge->getEndVertex()->tag();
695     if(cv1 != v1 && cv1 != v2 && cv2 != v1 && cv2 != v2) {
696       if(op == nullptr) { op = ge; }
697       else { /* already found ? should not happen */
698         return nullptr;
699       }
700     }
701   }
702   if(!edgeInside) return nullptr;
703   return op;
704 }
705 
build_chords(const std::set<GFace * > & faces,std::vector<std::set<GEdge * >> & chords,double maxDiffRel,bool ignoreEmbedded=false)706 void build_chords(const std::set<GFace *> &faces,
707                   std::vector<std::set<GEdge *> > &chords, double maxDiffRel,
708                   bool ignoreEmbedded = false)
709 {
710   /* Connectivity */
711   std::map<GEdge *, std::vector<GFace *> > edge2faces;
712   for(GFace *gf : faces) {
713     if(!ignoreEmbedded &&
714        (gf->embeddedEdges().size() > 0 || gf->embeddedVertices().size() > 0))
715       continue;
716     for(GEdge *ge : gf->edges()) { edge2faces[ge].push_back(gf); }
717   }
718 
719   Msg::Debug("build chords: %li faces, %li edges", faces.size(),
720              edge2faces.size());
721 
722   std::map<GEdge *, bool> done;
723   for(auto &kv : edge2faces) {
724     GEdge *geInit = kv.first;
725     if(done.find(geInit) != done.end()) continue;
726 
727     /* Breath first search starting from a GEdge */
728     std::queue<GEdge *> Q;
729     Q.push(geInit);
730     done[geInit] = true;
731 
732     std::set<GEdge *> chord;
733     while(Q.size() > 0) {
734       GEdge *ge = Q.front();
735       Q.pop();
736       int n = meshGEdgeTargetNumberOfPoints(ge);
737       chord.insert(ge);
738       for(GFace *gf : edge2faces[ge]) {
739         GEdge *ge2 = quad_face_opposite_edge(gf, ge);
740         if(ge2 && done.find(ge2) == done.end()) {
741           int n2 = meshGEdgeTargetNumberOfPoints(ge2);
742           if(double(std::abs(n2 - n)) / double(std::max(n, n2)) > maxDiffRel) {
743             continue;
744           }
745           Q.push(ge2);
746           done[ge2] = true;
747         }
748       }
749     }
750 
751     if(chord.size() >= 2) { chords.push_back(chord); }
752   }
753 }
754 
MeshSetTransfiniteFacesAutomatic(std::set<GFace * > & candidate_faces,double cornerAngle,bool setRecombine,double maxDiffRel,bool ignoreEmbedded)755 bool MeshSetTransfiniteFacesAutomatic(std::set<GFace *> &candidate_faces,
756                                       double cornerAngle, bool setRecombine,
757                                       double maxDiffRel, bool ignoreEmbedded)
758 {
759   /* Filter with topology and geometry */
760   std::set<GFace *> faces;
761   for(GFace *gf : candidate_faces) {
762     if(!ignoreEmbedded &&
763        (gf->embeddedEdges().size() > 0 || gf->embeddedVertices().size() > 0))
764       continue;
765     if(faceIsValidQuad(gf, cornerAngle)) { faces.insert(gf); }
766   }
767 
768   /* Build the topological chords in quad structure */
769   Msg::Debug(
770     "transfinite automatic: building chords from %li quadrangular faces ...",
771     faces.size());
772   std::vector<std::set<GEdge *> > chords;
773   build_chords(faces, chords, maxDiffRel);
774   Msg::Debug("... found %li chords", chords.size());
775 
776   /* Determine the number of points, set the transfinite curves */
777   Msg::Debug("transfinite automatic: assigning number of points ...");
778   std::size_t ne = 0;
779   for(std::set<GEdge *> &chord : chords) {
780     double avgNbPoints = 0;
781     for(GEdge *ge : chord) {
782       int n = meshGEdgeTargetNumberOfPoints(ge);
783       avgNbPoints += double(n);
784     }
785     avgNbPoints /= chord.size();
786 
787     int N = int(avgNbPoints + 0.5);
788     if(N == 0) N = 2;
789     if(N % 2 == 1) N = N + 1;
790 
791     Msg::Debug("- chord with %li edges -> %i points\n", chord.size(), N);
792 
793     for(GEdge *ge : chord) {
794       ge->meshAttributes.method = MESH_TRANSFINITE;
795       ge->meshAttributes.nbPointsTransfinite = N;
796       ge->meshAttributes.typeTransfinite = 1; /* Progression */
797       ge->meshAttributes.coeffTransfinite = 1.;
798       if(CTX::instance()->mesh.algo2d == ALGO_2D_QUAD_QUASI_STRUCT) {
799         ge->meshAttributes.typeTransfinite = 4; /* Use size map */
800       }
801       ne += 1;
802     }
803   }
804   Msg::Debug("transfinite automatic: transfinite set on %li edges", ne);
805 
806   std::size_t nf = 0;
807   for(GFace *gf : faces) {
808     bool transfinite = true;
809     std::vector<int> nPoints(4, 0);
810     std::size_t count = 0;
811     for(GEdge *ge : gf->edges()) {
812       if(ge->meshAttributes.method != MESH_TRANSFINITE) {
813         transfinite = false;
814         break;
815       }
816       nPoints[count] = ge->meshAttributes.nbPointsTransfinite;
817       count += 1;
818     }
819     if(transfinite && nPoints[0] == nPoints[2] && nPoints[1] == nPoints[3]) {
820       nf += 1;
821       gf->meshAttributes.method = MESH_TRANSFINITE;
822       gf->meshAttributes.transfiniteArrangement = 1; /* Right */
823       if(setRecombine) {
824         gf->meshAttributes.recombine = 1;
825         gf->meshAttributes.recombineAngle = 45.;
826       }
827     }
828   }
829   Msg::Debug("transfinite automatic: transfinite set on %li faces", nf);
830   return true;
831 }
832