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