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 <vector>
7 #include <stack>
8 #include <queue>
9 #include "OS.h"
10 #include "GmshConfig.h"
11 #include "gmshCrossFields.h"
12 #include "GModel.h"
13 #include "GFace.h"
14 #include "MEdge.h"
15 #include "MLine.h"
16 #include "MTriangle.h"
17 #include "GmshMessage.h"
18 #include "Context.h"
19 #include "meshGFaceOptimize.h"
20 #include "discreteEdge.h"
21 #include "Numeric.h"
22 #include "GModelParametrize.h"
23
24 #if defined(HAVE_POST)
25 #include "PView.h"
26 #include "PViewDataGModel.h"
27 #endif
28
29 #if defined(HAVE_SOLVER) && defined(HAVE_POST)
30
31 #include "dofManager.h"
32 #include "laplaceTerm.h"
33 #include "linearSystemGmm.h"
34 #include "linearSystemCSR.h"
35 #include "linearSystemFull.h"
36 #include "linearSystemPETSc.h"
37
lifting(double a,double _a)38 static inline double lifting(double a, double _a)
39 {
40 double D = M_PI * .5;
41 if(fabs(_a - a) < fabs(_a - (a + D)) && fabs(_a - a) < fabs(_a - (a - D))) {
42 return a;
43 }
44 else if(fabs(_a - (a + D)) < fabs(_a - a) &&
45 fabs(_a - (a + D)) < fabs(_a - (a - D))) {
46 return a + D;
47 }
48 else {
49 return a - D;
50 }
51 }
52
compat_orientation_extrinsic(const double * o0,const double * n0,const double * o1,const double * n1,double * a1,double * b1)53 static inline double compat_orientation_extrinsic(const double *o0,
54 const double *n0,
55 const double *o1,
56 const double *n1, double *a1,
57 double *b1)
58 {
59 double t0[3] = {n0[1] * o0[2] - n0[2] * o0[1], n0[2] * o0[0] - n0[0] * o0[2],
60 n0[0] * o0[1] - n0[1] * o0[0]};
61 double t1[3] = {n1[1] * o1[2] - n1[2] * o1[1], n1[2] * o1[0] - n1[0] * o1[2],
62 n1[0] * o1[1] - n1[1] * o1[0]};
63
64 const size_t permuts[8][2] = {{0, 0}, {1, 0}, {2, 0}, {3, 0},
65 {0, 1}, {1, 1}, {2, 1}, {3, 1}};
66 const double A[4][3] = {{o0[0], o0[1], o0[2]},
67 {t0[0], t0[1], t0[2]},
68 {-o0[0], -o0[1], -o0[2]},
69 {-t0[0], -t0[1], -t0[2]}};
70 const double B[2][3] = {{o1[0], o1[1], o1[2]}, {t1[0], t1[1], t1[2]}};
71
72 double maxx = -1;
73 int index = 0;
74 for(size_t i = 0; i < 8; i++) {
75 const size_t II = permuts[i][0];
76 const size_t JJ = permuts[i][1];
77 const double xx =
78 A[II][0] * B[JJ][0] + A[II][1] * B[JJ][1] + A[II][2] * B[JJ][2];
79 if(xx > maxx) {
80 index = i;
81 maxx = xx;
82 }
83 }
84 a1[0] = A[permuts[index][0]][0];
85 a1[1] = A[permuts[index][0]][1];
86 a1[2] = A[permuts[index][0]][2];
87 b1[0] = B[permuts[index][1]][0];
88 b1[1] = B[permuts[index][1]][1];
89 b1[2] = B[permuts[index][1]][2];
90 // b1 = B[permuts[index][1]];
91 return maxx;
92 }
93
compat_orientation_extrinsic(const SVector3 & o0,const SVector3 & n0,const SVector3 & o1,const SVector3 & n1,SVector3 & a1,SVector3 & b1)94 static inline double compat_orientation_extrinsic(const SVector3 &o0,
95 const SVector3 &n0,
96 const SVector3 &o1,
97 const SVector3 &n1,
98 SVector3 &a1, SVector3 &b1)
99 {
100 SVector3 t0 = crossprod(n0, o0);
101 SVector3 t1 = crossprod(n1, o1);
102
103 const size_t permuts[8][2] = {{0, 0}, {1, 0}, {2, 0}, {3, 0},
104 {0, 1}, {1, 1}, {2, 1}, {3, 1}};
105 SVector3 A[4]{o0, t0, -o0, -t0};
106 SVector3 B[2]{o1, t1};
107
108 double maxx = -1;
109 int index = 0;
110 for(size_t i = 0; i < 8; i++) {
111 const double xx = dot(A[permuts[i][0]], B[permuts[i][1]]);
112 if(xx > maxx) {
113 index = i;
114 maxx = xx;
115 }
116 }
117 a1 = A[permuts[index][0]];
118 b1 = B[permuts[index][1]];
119 return maxx;
120 }
121
122 class cross2d {
123 public:
124 MEdge _e;
125 bool inCutGraph;
126 bool inBoundary;
127 bool inInternalBoundary;
128 bool rotation;
129 // int cutGraphPart;
130 size_t counter;
131 SVector3 o_i; // orientation vector
132 SVector3 _nrml, _tgt, _tgt2; // local system of coordinates
133 std::vector<MEdge> _neighbors;
134 std::vector<cross2d *> _cneighbors;
135 double da[4];
136 double _cc[4], _ss[4];
137 // euler angles
138 double _a, _b, _c;
139 double _atemp, _btemp, _ctemp;
140 std::vector<MTriangle *> _t;
cross2d(MEdge & e,MTriangle * r,MEdge & e1,MEdge & e2)141 cross2d(MEdge &e, MTriangle *r, MEdge &e1, MEdge &e2)
142 : _e(e), inCutGraph(false), inBoundary(false), inInternalBoundary(false),
143 _a(0), _b(0), _c(0)
144 {
145 _t.push_back(r);
146 _neighbors.push_back(e1);
147 _neighbors.push_back(e2);
148 }
normalize(double & a)149 void normalize(double &a)
150 {
151 double D = M_PI * .5;
152 if(a < 0)
153 while(a < 0) a += D;
154 if(a >= D)
155 while(a >= D) a -= D;
156 }
finish2()157 void finish2()
158 {
159 if(_cneighbors.size() == 4) {
160 SVector3 x(0, 1, 0);
161 _a = _atemp = atan2(dot(_tgt2, x), dot(_tgt, x));
162 if(!inBoundary && !inInternalBoundary) {
163 _b = _btemp = sin(4 * _a);
164 _c = _ctemp = cos(4 * _a);
165 }
166 else {
167 _b = _btemp = 0;
168 _c = _ctemp = 1;
169 }
170 for(size_t i = 0; i < 4; i++) {
171 da[i] = atan2(dot(_tgt2, _cneighbors[i]->_tgt),
172 dot(_tgt, _cneighbors[i]->_tgt));
173 _cc[i] = cos(4 * da[i]);
174 _ss[i] = sin(4 * da[i]);
175 }
176 }
177 }
178
finish(std::map<MEdge,cross2d,MEdgeLessThan> & C)179 void finish(std::map<MEdge, cross2d, MEdgeLessThan> &C)
180 {
181 _tgt = SVector3(1, 0, 0);
182 _tgt2 = SVector3(0, 1, 0);
183 if(_t.size() <= 2) {
184 SVector3 v10(_t[0]->getVertex(1)->x() - _t[0]->getVertex(0)->x(),
185 _t[0]->getVertex(1)->y() - _t[0]->getVertex(0)->y(),
186 _t[0]->getVertex(1)->z() - _t[0]->getVertex(0)->z());
187 SVector3 v20(_t[0]->getVertex(2)->x() - _t[0]->getVertex(0)->x(),
188 _t[0]->getVertex(2)->y() - _t[0]->getVertex(0)->y(),
189 _t[0]->getVertex(2)->z() - _t[0]->getVertex(0)->z());
190 SVector3 xx = crossprod(v20, v10);
191 xx.normalize();
192 SVector3 yy = xx;
193 if(_t.size() == 2) {
194 SVector3 v10b(_t[1]->getVertex(1)->x() - _t[1]->getVertex(0)->x(),
195 _t[1]->getVertex(1)->y() - _t[1]->getVertex(0)->y(),
196 _t[1]->getVertex(1)->z() - _t[1]->getVertex(0)->z());
197 SVector3 v20b(_t[1]->getVertex(2)->x() - _t[1]->getVertex(0)->x(),
198 _t[1]->getVertex(2)->y() - _t[1]->getVertex(0)->y(),
199 _t[1]->getVertex(2)->z() - _t[1]->getVertex(0)->z());
200 yy = crossprod(v20b, v10b);
201 yy.normalize();
202 // if(dot(xx, yy) < .5) inInternalBoundary = 1;
203 }
204 _nrml = xx + yy;
205 _nrml.normalize();
206 SVector3 vt(_e.getVertex(1)->x() - _e.getVertex(0)->x(),
207 _e.getVertex(1)->y() - _e.getVertex(0)->y(),
208 _e.getVertex(1)->z() - _e.getVertex(0)->z());
209 _tgt = vt;
210 _tgt.normalize();
211 _tgt2 = crossprod(_nrml, _tgt);
212 }
213
214 if(_t.size() == 1) { inBoundary = true; }
215 else if(_t.size() >= 2) {
216 if(inBoundary) {
217 inBoundary = false;
218 inInternalBoundary = true;
219 }
220 }
221
222 for(size_t i = 0; i < _neighbors.size(); i++) {
223 auto it = C.find(_neighbors[i]);
224 if(it == C.end())
225 Msg::Error("impossible situation");
226 else
227 _cneighbors.push_back(&(it->second));
228 }
229 if(_cneighbors.size() != 4) {
230 _a = 0;
231 _atemp = _a;
232 }
233 _neighbors.clear();
234 _b = _btemp = sin(4 * _a);
235 _c = _ctemp = cos(4 * _a);
236 }
average_init()237 double average_init()
238 {
239 if(!inBoundary && !inInternalBoundary) {
240 _btemp = 0;
241 _ctemp = 0;
242 for(int i = 0; i < 4; i++) {
243 _ctemp += (_cneighbors[i]->_c * _cc[i] - _cneighbors[i]->_b * _ss[i]);
244 _btemp += (_cneighbors[i]->_c * _ss[i] + _cneighbors[i]->_b * _cc[i]);
245 }
246 _btemp *= 0.25;
247 _ctemp *= 0.25;
248 _b = _btemp;
249 _c = _ctemp;
250 }
251 return 1;
252 }
253
grad()254 double grad()
255 {
256 if(!inBoundary && !inInternalBoundary) {
257 double D = M_PI * .5;
258 double a[4] = {_cneighbors[0]->_a, _cneighbors[1]->_a, _cneighbors[2]->_a,
259 _cneighbors[3]->_a};
260 for(int i = 0; i < 4; i++) {
261 a[i] += da[i];
262 normalize(a[i]);
263 }
264
265 double b[4];
266 for(int i = 0; i < 4; i++) {
267 if(fabs(_a - a[i]) < fabs(_a - (a[i] + D)) &&
268 fabs(_a - a[i]) < fabs(_a - (a[i] - D))) {
269 b[i] = a[i];
270 }
271 else if(fabs(_a - (a[i] + D)) < fabs(_a - (a[i])) &&
272 fabs(_a - (a[i] + D)) < fabs(_a - (a[i] - D))) {
273 b[i] = a[i] + D;
274 }
275 else {
276 b[i] = a[i] - D;
277 }
278 }
279 return fabs(_a - b[0]) + fabs(_a - b[1]) + fabs(_a - b[2]) +
280 fabs(_a - b[3]);
281 }
282 return 0;
283 }
284
lifting(double a)285 double lifting(double a)
286 {
287 double D = M_PI * .5;
288 if(fabs(_a - a) < fabs(_a - (a + D)) && fabs(_a - a) < fabs(_a - (a - D))) {
289 return a;
290 }
291 else if(fabs(_a - (a + D)) < fabs(_a - a) &&
292 fabs(_a - (a + D)) < fabs(_a - (a - D))) {
293 return a + D;
294 }
295 else {
296 return a - D;
297 }
298 }
299
computeVector()300 void computeVector()
301 {
302 _a = _atemp = 0.25 * atan2(_btemp, _ctemp);
303 normalize(_atemp);
304 o_i = _tgt * cos(_atemp) + _tgt2 * sin(_atemp);
305 }
306
computeAngle()307 void computeAngle()
308 {
309 if(_cneighbors.size() != 4) return;
310 double _anew = atan2(dot(_tgt2, o_i), dot(_tgt, o_i));
311 normalize(_anew);
312 _a = _atemp = _anew;
313 }
314
average2()315 double average2()
316 {
317 // TEMPORARY VERSION, slow but correct
318 // Instant field-aligned meshes
319 if(_cneighbors.size() != 4) return 0;
320 double weight = 0.0;
321 SVector3 o_i_old = o_i;
322 SVector3 n_i = _nrml;
323 for(int i = 0; i < 4; i++) {
324 SVector3 o_j = _cneighbors[i]->o_i;
325 SVector3 n_j = _cneighbors[i]->_nrml;
326 SVector3 x, y;
327 compat_orientation_extrinsic(o_i, n_i, o_j, n_j, x, y);
328 o_i = x * weight + y;
329 o_i -= n_i * dot(o_i, n_i);
330 o_i.normalize();
331 weight += 1.0;
332 }
333 return std::min(1. - fabs(dot(o_i, o_i_old)), fabs(dot(o_i, o_i_old)));
334 }
335
average()336 double average()
337 {
338 if(_cneighbors.size() == 4) {
339 double D = M_PI * .5;
340 double a[4] = {_cneighbors[0]->_a, _cneighbors[1]->_a, _cneighbors[2]->_a,
341 _cneighbors[3]->_a};
342 for(int i = 0; i < 4; i++) {
343 a[i] += da[i];
344 normalize(a[i]);
345 }
346
347 double b[4];
348 double avg = 0.0;
349 for(int i = 0; i < 4; i++) {
350 if(fabs(_a - a[i]) < fabs(_a - (a[i] + D)) &&
351 fabs(_a - a[i]) < fabs(_a - (a[i] - D))) {
352 b[i] = a[i];
353 }
354 else if(fabs(_a - (a[i] + D)) < fabs(_a - (a[i])) &&
355 fabs(_a - (a[i] + D)) < fabs(_a - (a[i] - D))) {
356 b[i] = a[i] + D;
357 }
358 else {
359 b[i] = a[i] - D;
360 }
361 }
362 avg = 0.25 * (b[0] + b[1] + b[2] + b[3]);
363
364 normalize(avg);
365
366 double d = fabs(_a - avg);
367 _a = _atemp = avg;
368
369 return d;
370 }
371 return 0;
372 }
373 };
374
375 struct cross2dPtrLessThan {
376 MEdgeLessThan l;
operator ()cross2dPtrLessThan377 bool operator()(const cross2d *v1, const cross2d *v2) const
378 {
379 return l(v1->_e, v2->_e);
380 }
381 };
382
383 // ---------------------------------------------
384 // TODO : MAKE IT PARALLEL AND SUPERFAST
385 // DO IT ON SURFACES
386 // ---------------------------------------------
387
closest(const cross2d & c1,const cross2d & c2,double & a2,double & diff)388 static void closest(const cross2d &c1, const cross2d &c2, double &a2,
389 double &diff)
390 {
391 SVector3 d = c1._tgt * cos(c1._atemp) + c1._tgt2 * sin(c1._atemp);
392
393 double P = M_PI / 2;
394
395 SVector3 d1 = c2._tgt * cos(c2._atemp) + c2._tgt2 * sin(c2._atemp);
396 SVector3 d2 = c2._tgt * cos(c2._atemp + P) + c2._tgt2 * sin(c2._atemp + P);
397 SVector3 d3 =
398 c2._tgt * cos(c2._atemp + 2 * P) + c2._tgt2 * sin(c2._atemp + 2 * P);
399 SVector3 d4 =
400 c2._tgt * cos(c2._atemp + 3 * P) + c2._tgt2 * sin(c2._atemp + 3 * P);
401
402 double D1 = dot(d, d1);
403 double D2 = dot(d, d2);
404 double D3 = dot(d, d3);
405 double D4 = dot(d, d4);
406 diff = acos(D1);
407 if(D1 > D2 && D1 > D3 && D1 > D4) { a2 = c2._atemp; }
408 else if(D2 > D1 && D2 > D3 && D2 > D4) {
409 a2 = c2._atemp + P;
410 }
411 else if(D3 > D1 && D3 > D2 && D3 > D4) {
412 a2 = c2._atemp + 2 * P;
413 }
414 else {
415 a2 = c2._atemp + 3 * P;
416 }
417 }
418
computeLifting(cross2d * first,int branch,std::set<MEdge,MEdgeLessThan> & cutG,std::set<MVertex *,MVertexPtrLessThan> & sing,std::set<MVertex *,MVertexPtrLessThan> & bnd,std::set<cross2d * > & visited)419 static void computeLifting(cross2d *first, int branch,
420 std::set<MEdge, MEdgeLessThan> &cutG,
421 std::set<MVertex *, MVertexPtrLessThan> &sing,
422 std::set<MVertex *, MVertexPtrLessThan> &bnd,
423 std::set<cross2d *> &visited)
424 {
425 // store in _atemp the branch of the neighbor
426 std::set<MVertex *, MVertexPtrLessThan> cg;
427 {
428 auto it = cutG.begin();
429 for(; it != cutG.end(); ++it) {
430 cg.insert(it->getVertex(0));
431 cg.insert(it->getVertex(1));
432 }
433 }
434
435 FILE *_f = fopen("visited.pos", "w");
436 fprintf(_f, "View\"\"{\n");
437
438 std::queue<cross2d *> _s;
439 _s.push(first);
440 first->_atemp = first->_a + branch * M_PI / 2.0;
441 first->_btemp = 10000.;
442 visited.insert(first);
443 while(!_s.empty()) {
444 cross2d *c = _s.front();
445 _s.pop();
446 if(cutG.find(c->_e) == cutG.end()) {
447 for(size_t i = 0; i < c->_cneighbors.size(); i++) {
448 double a2, diff;
449 cross2d *n = c->_cneighbors[i];
450 closest(*c, *n, a2, diff);
451 if(n->_btemp < 1000) {
452 n->_btemp = 10000;
453
454 bool s0 = sing.find(n->_e.getVertex(0)) != sing.end();
455 bool s1 = sing.find(n->_e.getVertex(1)) != sing.end();
456 bool c0 = cg.find(n->_e.getVertex(0)) != cg.end();
457 bool c1 = cg.find(n->_e.getVertex(1)) != cg.end();
458 bool b0 = bnd.find(n->_e.getVertex(0)) != bnd.end();
459 bool b1 = bnd.find(n->_e.getVertex(1)) != bnd.end();
460
461 s0 = s1 = false;
462 // b0 = b1 = false;
463 if(((s0 && c1) || (s0 && b1)) || ((s1 && c0) || (s1 && b0))) {
464 printf("BLOCKED \n");
465 }
466
467 if((!s0 && !s1)) _s.push(n);
468 // printf("%12.5E %12.5E %12.5E\n",n->a,n->_atemp,c->_atemp);
469 n->_atemp = a2;
470 visited.insert(n);
471 SVector3 d = n->_tgt * cos(n->_atemp) + n->_tgt2 * sin(n->_atemp);
472 fprintf(_f, "VL(%g,%g,%g,%g,%g,%g){%g,%g,%g,%g,%g,%g};\n",
473 n->_e.getVertex(0)->x(), n->_e.getVertex(0)->y(),
474 n->_e.getVertex(0)->z(), n->_e.getVertex(1)->x(),
475 n->_e.getVertex(1)->y(), n->_e.getVertex(1)->z(), d.x(),
476 d.y(), d.z(), d.x(), d.y(), d.z());
477 }
478 }
479 }
480 }
481 fprintf(_f, "};\n");
482 fclose(_f);
483 }
484
485 struct groupOfCross2d {
486 int groupId;
487 bool rot;
488 double mat[2][2];
489 std::vector<MVertex *> vertices;
490 std::vector<MVertex *> singularities;
491 std::vector<MVertex *> left;
492 std::vector<MVertex *> right;
493 std::vector<cross2d *> crosses;
494 std::vector<MTriangle *> side;
printgroupOfCross2d495 void print(FILE *f)
496 {
497 for(size_t i = 0; i < crosses.size(); i++) {
498 fprintf(
499 f, "SL(%g,%g,%g,%g,%g,%g){%d,%d};\n", crosses[i]->_e.getVertex(0)->x(),
500 crosses[i]->_e.getVertex(0)->y(), crosses[i]->_e.getVertex(0)->z(),
501 crosses[i]->_e.getVertex(1)->x(), crosses[i]->_e.getVertex(1)->y(),
502 crosses[i]->_e.getVertex(1)->z(), groupId, groupId);
503 }
504 for(size_t i = 0; i < side.size(); i++) {
505 fprintf(f, "ST(%g,%g,%g,%g,%g,%g,%g,%g,%g){%d,%d,%d};\n",
506 side[i]->getVertex(0)->x(), side[i]->getVertex(0)->y(),
507 side[i]->getVertex(0)->z(), side[i]->getVertex(1)->x(),
508 side[i]->getVertex(1)->y(), side[i]->getVertex(1)->z(),
509 side[i]->getVertex(2)->x(), side[i]->getVertex(2)->y(),
510 side[i]->getVertex(2)->z(), groupId, groupId, groupId);
511 }
512 };
groupOfCross2dgroupOfCross2d513 groupOfCross2d(int id) : groupId(id) {}
514 };
515
unDuplicateNodesInCutGraph(std::vector<GFace * > & f,std::map<MVertex *,MVertex *,MVertexPtrLessThan> & new2old)516 static void unDuplicateNodesInCutGraph(
517 std::vector<GFace *> &f,
518 std::map<MVertex *, MVertex *, MVertexPtrLessThan> &new2old)
519 {
520 for(size_t i = 0; i < f.size(); i++) {
521 for(size_t j = 0; j < f[i]->triangles.size(); j++) {
522 MTriangle *t = f[i]->triangles[j];
523 for(size_t k = 0; k < 3; k++) {
524 auto it = new2old.find(t->getVertex(k));
525 if(it != new2old.end()) t->setVertex(k, it->second);
526 }
527 }
528 }
529 }
530
duplicateNodesInCutGraph(std::vector<GFace * > & f,std::map<MEdge,cross2d,MEdgeLessThan> & C,std::map<MVertex *,MVertex *,MVertexPtrLessThan> & new2old,std::multimap<MVertex *,MVertex *,MVertexPtrLessThan> & old2new,std::map<MEdge,MEdge,MEdgeLessThan> & duplicateEdges,std::set<MVertex *,MVertexPtrLessThan> & sing,v2t_cont & adj,std::vector<groupOfCross2d> & G)531 static void duplicateNodesInCutGraph(
532 std::vector<GFace *> &f, std::map<MEdge, cross2d, MEdgeLessThan> &C,
533 std::map<MVertex *, MVertex *, MVertexPtrLessThan> &new2old,
534 std::multimap<MVertex *, MVertex *, MVertexPtrLessThan> &old2new,
535 std::map<MEdge, MEdge, MEdgeLessThan> &duplicateEdges,
536 std::set<MVertex *, MVertexPtrLessThan> &sing, v2t_cont &adj,
537 std::vector<groupOfCross2d> &G)
538 {
539 FILE *_f = fopen("nodes.pos", "w");
540 fprintf(_f, "View \" nodes \"{\n");
541
542 auto it = adj.begin();
543 std::set<MElement *> touched;
544 std::set<MVertex *, MVertexPtrLessThan> vtouched;
545
546 std::vector<std::pair<MElement *, std::pair<int, MVertex *> > > replacements;
547
548 while(it != adj.end()) {
549 std::vector<MElement *> els = it->second;
550 int ITER = 0;
551 while(!els.empty()) {
552 std::vector<MElement *> _side;
553 std::stack<MElement *> _s;
554 _s.push(els[0]);
555 _side.push_back(els[0]);
556 els.erase(els.begin());
557 while(!_s.empty()) {
558 MElement *t = _s.top();
559 _s.pop();
560 for(int i = 0; i < 3; i++) {
561 MEdge e0 = t->getEdge(i);
562 auto it0 = C.find(e0);
563 if(!it0->second.inCutGraph) {
564 for(size_t j = 0; j < it0->second._t.size(); j++) {
565 auto ite = std::find(els.begin(), els.end(), it0->second._t[j]);
566 if(ite != els.end()) {
567 els.erase(ite);
568 _side.push_back(it0->second._t[j]);
569 _s.push(it0->second._t[j]);
570 }
571 }
572 }
573 }
574 }
575
576 if(ITER) {
577 // if (it->first->getNum() != 268){
578 MVertex *v =
579 new MVertex(it->first->x(), it->first->y(), it->first->z(), f[0]);
580 std::pair<MVertex *, MVertex *> p = std::make_pair(it->first, v);
581 old2new.insert(p);
582 f[0]->mesh_vertices.push_back(v);
583 for(size_t i = 0; i < _side.size(); i++) {
584 for(size_t j = 0; j < 3; j++) {
585 if(_side[i]->getVertex(j) == it->first) {
586 std::pair<int, MVertex *> r = std::make_pair(j, v);
587 std::pair<MElement *, std::pair<int, MVertex *> > r2 =
588 std::make_pair(_side[i], r);
589 replacements.push_back(r2);
590 }
591 }
592 }
593 fprintf(_f, "SP(%g,%g,%g){%d};\n", it->first->x(), it->first->y(),
594 it->first->z(), (int)els.size());
595 // printf("found vertex with %d on one side\n",els.size());
596 // }
597 }
598 ++ITER;
599 }
600 // if (ITER > 1)printf("ITER = %d\n",ITER);
601 ++it;
602 }
603 fprintf(_f, "};\n");
604 fclose(_f);
605
606 for(size_t i = 0; i < replacements.size(); i++) {
607 MElement *e = replacements[i].first;
608 int j = replacements[i].second.first;
609 MVertex *v = replacements[i].second.second;
610 MVertex *old = e->getVertex(j);
611 for(int j = 0; j < e->getNumEdges(); j++) {
612 MEdge ed = e->getEdge(j);
613 if(ed.getVertex(0) == old) {
614 duplicateEdges[ed] = MEdge(v, ed.getVertex(1));
615 }
616 else if(ed.getVertex(1) == old) {
617 duplicateEdges[ed] = MEdge(ed.getVertex(0), v);
618 }
619 }
620 new2old[v] = old;
621 e->setVertex(j, v);
622 }
623 }
624
625 static void
computeUniqueVectorPerTriangle(GModel * gm,std::vector<GFace * > & f,std::map<MEdge,cross2d,MEdgeLessThan> & C,int dir,std::map<MTriangle *,SVector3> & d)626 computeUniqueVectorPerTriangle(GModel *gm, std::vector<GFace *> &f,
627 std::map<MEdge, cross2d, MEdgeLessThan> &C,
628 int dir, std::map<MTriangle *, SVector3> &d)
629 {
630 for(size_t i = 0; i < f.size(); i++) {
631 for(size_t j = 0; j < f[i]->triangles.size(); j++) {
632 MTriangle *t = f[i]->triangles[j];
633 SVector3 a0(0, 0, 0);
634 int n = 0;
635 for(int k = 0; k < 3; k++) {
636 MEdge e = t->getEdge(k);
637 auto it = C.find(e);
638 if(it != C.end()) {
639 if(!it->second.inCutGraph) {
640 n++;
641 double C =
642 (dir == 1) ? cos(it->second._atemp) : sin(it->second._atemp);
643 double S =
644 (dir == 1) ? sin(it->second._atemp) : -cos(it->second._atemp);
645 SVector3 aa = it->second._tgt * C + it->second._tgt2 * S;
646 a0 += aa;
647 }
648 }
649 else {
650 printf("ERROR\n");
651 }
652 }
653 if(n) a0.normalize();
654 d[t] = a0;
655 }
656 }
657 }
658
659 int ZERO_ = -8;
660
createLagrangeMultipliers(dofManager<double> & myAssembler,groupOfCross2d & g)661 static void createLagrangeMultipliers(dofManager<double> &myAssembler,
662 groupOfCross2d &g)
663 {
664 // return;
665 // if (g.crosses[0]->inInternalBoundary)return;
666
667 if(g.singularities.size() == 1) {
668 // printf("group id %d singularity %lu (%lu)\n", g.groupId,
669 // g.singularities[0]->getNum(), g.singularities.size());
670 myAssembler.numberVertex(g.singularities[0], 0, 33);
671 myAssembler.numberVertex(g.singularities[0], 0, 34);
672 }
673 else {
674 // printf("group id %d %lu singularities \n", g.groupId,
675 // g.singularities.size());
676 }
677
678 if(g.groupId == ZERO_) {
679 myAssembler.numberVertex(g.left[0], 0, 3 + 10000 * g.groupId);
680 }
681
682 for(size_t K = 1; K < g.left.size(); K++) {
683 myAssembler.numberVertex(g.left[K], 0, 3 + 100 * g.groupId);
684 myAssembler.numberVertex(g.left[K], 0, 4 + 100 * g.groupId);
685 }
686 }
687
LagrangeMultipliers3(dofManager<double> & myAssembler,groupOfCross2d & g,std::map<MTriangle *,SVector3> & d0,bool assemble)688 static void LagrangeMultipliers3(dofManager<double> &myAssembler,
689 groupOfCross2d &g,
690 std::map<MTriangle *, SVector3> &d0,
691 bool assemble)
692 {
693 return;
694 if(g.crosses[0]->inInternalBoundary) {
695 if(!assemble) {
696 // printf("group %d is internal\n",g.groupId);
697 for(size_t K = 1; K < g.left.size(); K++) {
698 myAssembler.numberVertex(g.left[K], 0, 12 + 100 * g.groupId);
699 myAssembler.numberVertex(g.left[K], 0, 13 + 100 * g.groupId);
700 }
701 return;
702 }
703
704 MTriangle *t0 = g.crosses[0]->_t[0];
705 MTriangle *t1 = g.crosses[0]->_t[1];
706 SVector3 dir0 = d0[t0];
707 SVector3 dir1 = d0[t1];
708
709 if(std::find(g.side.begin(), g.side.end(), t1) != g.side.end()) {
710 dir0 = d0[t1];
711 dir1 = d0[t0];
712 }
713
714 printf("%12.5E %12.5E\n", fabs(dot(g.crosses[0]->_tgt, dir0)),
715 fabs(dot(g.crosses[0]->_tgt, dir1)));
716
717 Dof U1R(g.left[0]->getNum(), Dof::createTypeWithTwoInts(0, 1));
718 Dof U2R(g.right[0]->getNum(), Dof::createTypeWithTwoInts(0, 1));
719 Dof V1R(g.left[0]->getNum(), Dof::createTypeWithTwoInts(0, 2));
720 Dof V2R(g.right[0]->getNum(), Dof::createTypeWithTwoInts(0, 2));
721 for(size_t K = 1; K < g.left.size(); K++) {
722 Dof E1(g.left[K]->getNum(),
723 Dof::createTypeWithTwoInts(0, 12 + 100 * g.groupId));
724 Dof E2(g.left[K]->getNum(),
725 Dof::createTypeWithTwoInts(0, 13 + 100 * g.groupId));
726 Dof U1(g.left[K]->getNum(), Dof::createTypeWithTwoInts(0, 1));
727 Dof U2(g.right[K]->getNum(), Dof::createTypeWithTwoInts(0, 1));
728 Dof V1(g.left[K]->getNum(), Dof::createTypeWithTwoInts(0, 2));
729 Dof V2(g.right[K]->getNum(), Dof::createTypeWithTwoInts(0, 2));
730
731 if(fabs(dot(g.crosses[0]->_tgt, dir0)) > .5) {
732 myAssembler.assemble(E1, U1, 1.0);
733 myAssembler.assemble(U1, E1, 1.0);
734 myAssembler.assemble(E1, U1R, -1.0);
735 myAssembler.assemble(U1R, E1, -1.0);
736 }
737 else {
738 myAssembler.assemble(E1, V1, 1.0);
739 myAssembler.assemble(V1, E1, 1.0);
740 myAssembler.assemble(E1, V1R, -1.0);
741 myAssembler.assemble(V1R, E1, -1.0);
742 }
743
744 if(fabs(dot(g.crosses[0]->_tgt, dir1)) > .5) {
745 // printf("HAHAHA %d\n",g.groupId);
746 myAssembler.assemble(E2, U2, 1.0);
747 myAssembler.assemble(U2, E2, 1.0);
748 myAssembler.assemble(E2, U2R, -1.0);
749 myAssembler.assemble(U2R, E2, -1.0);
750 }
751 else {
752 // printf("HAHAHO %d\n",g.groupId);
753 myAssembler.assemble(E2, V2, 1.0);
754 myAssembler.assemble(V2, E2, 1.0);
755 myAssembler.assemble(E2, V2R, -1.0);
756 myAssembler.assemble(V2R, E2, -1.0);
757 }
758 }
759 }
760 }
761
assembleLagrangeMultipliers(dofManager<double> & myAssembler,groupOfCross2d & g)762 static void assembleLagrangeMultipliers(dofManager<double> &myAssembler,
763 groupOfCross2d &g)
764 {
765 Dof U1R(g.left[0]->getNum(), Dof::createTypeWithTwoInts(0, 1));
766 Dof V1R(g.left[0]->getNum(), Dof::createTypeWithTwoInts(0, 2));
767 Dof U2R(g.right[0]->getNum(), Dof::createTypeWithTwoInts(0, 1));
768 Dof V2R(g.right[0]->getNum(), Dof::createTypeWithTwoInts(0, 2));
769
770 // if (g.groupId == 1){
771 // g.mat[0][0]=-1;
772 // g.mat[1][1]=-1;
773 // }
774
775 // printf("GROUP %d\n",g.groupId);
776 // printf("LEFT --- RIGHT\n");
777 // printf("%3lu %3lu\n",g.left[0]->getNum(),g.right[0]->getNum());
778
779 if(g.singularities.size() == 1) {
780 Dof E1(g.singularities[0]->getNum(), Dof::createTypeWithTwoInts(0, 33));
781 Dof E2(g.singularities[0]->getNum(), Dof::createTypeWithTwoInts(0, 34));
782 Dof U(g.singularities[0]->getNum(), Dof::createTypeWithTwoInts(0, 1));
783 Dof V(g.singularities[0]->getNum(), Dof::createTypeWithTwoInts(0, 2));
784
785 myAssembler.assemble(E1, U, 1.0);
786 myAssembler.assemble(E1, U1R, -1.0);
787
788 myAssembler.assemble(E1, U, -g.mat[0][0]);
789 myAssembler.assemble(E1, V, -g.mat[0][1]);
790 myAssembler.assemble(E1, U2R, g.mat[0][0]);
791 myAssembler.assemble(E1, V2R, g.mat[0][1]);
792
793 myAssembler.assemble(E2, V, 1.0);
794 myAssembler.assemble(E2, V1R, -1.0);
795
796 myAssembler.assemble(E2, U, -g.mat[1][0]);
797 myAssembler.assemble(E2, V, -g.mat[1][1]);
798 myAssembler.assemble(E2, U2R, g.mat[1][0]);
799 myAssembler.assemble(E2, V2R, g.mat[1][1]);
800
801 // sym
802
803 myAssembler.assemble(U, E1, 1.0);
804 myAssembler.assemble(U1R, E1, -1.0);
805
806 myAssembler.assemble(U, E1, -g.mat[0][0]);
807 myAssembler.assemble(V, E1, -g.mat[0][1]);
808 myAssembler.assemble(U2R, E1, g.mat[0][0]);
809 myAssembler.assemble(V2R, E1, g.mat[0][1]);
810
811 myAssembler.assemble(V, E2, 1.0);
812 myAssembler.assemble(V1R, E2, -1.0);
813
814 myAssembler.assemble(U, E2, -g.mat[1][0]);
815 myAssembler.assemble(V, E2, -g.mat[1][1]);
816 myAssembler.assemble(U2R, E2, g.mat[1][0]);
817 myAssembler.assemble(V2R, E2, g.mat[1][1]);
818 }
819
820 // TEST NO JUMP ON U for group 3 ...
821 if(g.groupId == ZERO_) {
822 Dof E1(g.left[0]->getNum(),
823 Dof::createTypeWithTwoInts(0, 3 + 10000 * g.groupId));
824 Dof U1(g.left[0]->getNum(), Dof::createTypeWithTwoInts(0, 2));
825 Dof U2(g.right[0]->getNum(), Dof::createTypeWithTwoInts(0, 1));
826 }
827
828 for(size_t K = 1; K < g.left.size(); K++) {
829 // printf("%3lu %3lu\n",g.left[K]->getNum(),g.right[K]->getNum());
830 // EQUATION IDS (Lagrange multipliers)
831 Dof E1(g.left[K]->getNum(),
832 Dof::createTypeWithTwoInts(0, 3 + 100 * g.groupId));
833 Dof E2(g.left[K]->getNum(),
834 Dof::createTypeWithTwoInts(0, 4 + 100 * g.groupId));
835
836 // DOF IDS
837 Dof U1(g.left[K]->getNum(), Dof::createTypeWithTwoInts(0, 1));
838 Dof V1(g.left[K]->getNum(), Dof::createTypeWithTwoInts(0, 2));
839 Dof U2(g.right[K]->getNum(), Dof::createTypeWithTwoInts(0, 1));
840 Dof V2(g.right[K]->getNum(), Dof::createTypeWithTwoInts(0, 2));
841
842 myAssembler.assemble(E1, U1, 1.0);
843 myAssembler.assemble(E1, U1R, -1.0);
844
845 myAssembler.assemble(E1, U2, -g.mat[0][0]);
846 myAssembler.assemble(E1, V2, -g.mat[0][1]);
847 myAssembler.assemble(E1, U2R, g.mat[0][0]);
848 myAssembler.assemble(E1, V2R, g.mat[0][1]);
849
850 myAssembler.assemble(E2, V1, 1.0);
851 myAssembler.assemble(E2, V1R, -1.0);
852
853 myAssembler.assemble(E2, U2, -g.mat[1][0]);
854 myAssembler.assemble(E2, V2, -g.mat[1][1]);
855 myAssembler.assemble(E2, U2R, g.mat[1][0]);
856 myAssembler.assemble(E2, V2R, g.mat[1][1]);
857
858 // sym
859
860 myAssembler.assemble(U1, E1, 1.0);
861 myAssembler.assemble(U1R, E1, -1.0);
862 myAssembler.assemble(U2, E1, -g.mat[0][0]);
863 myAssembler.assemble(V2, E1, -g.mat[0][1]);
864 myAssembler.assemble(U2R, E1, g.mat[0][0]);
865 myAssembler.assemble(V2R, E1, g.mat[0][1]);
866
867 myAssembler.assemble(V1, E2, 1.0);
868 myAssembler.assemble(V1R, E2, -1.0);
869 myAssembler.assemble(U2, E2, -g.mat[1][0]);
870 myAssembler.assemble(V2, E2, -g.mat[1][1]);
871 myAssembler.assemble(U2R, E2, g.mat[1][0]);
872 myAssembler.assemble(V2R, E2, g.mat[1][1]);
873 }
874 }
875
876 static void
LagrangeMultipliers2(dofManager<double> & myAssembler,int NUMDOF,std::map<MEdge,cross2d,MEdgeLessThan> & C,std::vector<std::vector<cross2d * >> & groups,std::map<MEdge,MEdge,MEdgeLessThan> & duplicateEdges,bool assemble,std::map<MTriangle *,SVector3> & lift)877 LagrangeMultipliers2(dofManager<double> &myAssembler, int NUMDOF,
878 std::map<MEdge, cross2d, MEdgeLessThan> &C,
879 std::vector<std::vector<cross2d *> > &groups,
880 std::map<MEdge, MEdge, MEdgeLessThan> &duplicateEdges,
881 bool assemble, std::map<MTriangle *, SVector3> &lift)
882 {
883 for(size_t i = 0; i < groups.size(); i++) {
884 size_t N = groups[i].size();
885 MEdge ed = groups[i][0]->_e;
886 auto ite = duplicateEdges.find(ed);
887 if(ite != duplicateEdges.end()) ed = ite->second;
888 MVertex *v = ed.getVertex(0);
889 auto it = C.find(groups[i][0]->_e);
890 SVector3 aaa = lift[it->second._t[0]];
891
892 double S = fabs(dot(it->second._tgt, aaa));
893
894 // printf("DIR %d S = %12.5E\n",NUMDOF, S);
895
896 if(S < .2 /*sqrt(2)/2.0*/) {
897 for(size_t j = 0; j < N; j++) {
898 ed = groups[i][j]->_e;
899 ite = duplicateEdges.find(ed);
900 if(ite != duplicateEdges.end()) ed = ite->second;
901 for(int k = 0; k < 2; k++) {
902 MVertex *vk = ed.getVertex(k);
903 if(vk != v) {
904 if(!assemble) { myAssembler.numberVertex(vk, 0, 5 + 100 * i); }
905 else {
906 Dof Eref(vk->getNum(),
907 Dof::createTypeWithTwoInts(0, 5 + 100 * i));
908 Dof Uref(vk->getNum(), Dof::createTypeWithTwoInts(0, NUMDOF));
909 Dof U(v->getNum(), Dof::createTypeWithTwoInts(0, NUMDOF));
910 myAssembler.assemble(Eref, Uref, 1.0);
911 myAssembler.assemble(Eref, U, -1.0);
912 myAssembler.assemble(Uref, Eref, 1.0);
913 myAssembler.assemble(U, Eref, -1.0);
914 }
915 }
916 }
917 }
918 }
919 }
920 }
921
922 struct cutGraphPassage {
923 int _id;
924 int _uv;
925 SPoint3 _p;
cutGraphPassagecutGraphPassage926 cutGraphPassage(int id, int uv, const SPoint3 &p) : _id(id), _uv(uv), _p(p) {}
927 };
928
createDofs(dofManager<double> & myAssembler,int NUMDOF,std::set<MVertex *,MVertexPtrLessThan> & vs)929 static void createDofs(dofManager<double> &myAssembler, int NUMDOF,
930 std::set<MVertex *, MVertexPtrLessThan> &vs)
931 {
932 for(auto it = vs.begin(); it != vs.end(); ++it)
933 myAssembler.numberVertex(*it, 0, NUMDOF);
934 }
935
createExtraConnexions(dofManager<double> & myAssembler,std::vector<groupOfCross2d> & G,std::vector<cutGraphPassage> & passages)936 void createExtraConnexions(dofManager<double> &myAssembler,
937 std::vector<groupOfCross2d> &G,
938 std::vector<cutGraphPassage> &passages)
939 {
940 return;
941 // give a number to the equation ...
942 myAssembler.numberVertex(G[0].left[0], 0, 10201020);
943 }
944
assembleExtraConnexions(dofManager<double> & myAssembler,std::vector<groupOfCross2d> & G,std::vector<cutGraphPassage> & passages)945 void assembleExtraConnexions(dofManager<double> &myAssembler,
946 std::vector<groupOfCross2d> &G,
947 std::vector<cutGraphPassage> &passages)
948 {
949 int nConn = 2;
950 int groups[2][2] = {{14, 1}, {13, 2}};
951
952 Dof E(G[0].left[0]->getNum(), Dof::createTypeWithTwoInts(0, 10201020));
953
954 for(int i = 0; i < nConn; i++) {
955 groupOfCross2d &g = G[groups[i][0]];
956 int index = groups[i][1];
957 Dof U1(g.left[0]->getNum(), Dof::createTypeWithTwoInts(0, index));
958 Dof U2(g.right[0]->getNum(), Dof::createTypeWithTwoInts(0, 1));
959 Dof V2(g.right[0]->getNum(), Dof::createTypeWithTwoInts(0, 2));
960 myAssembler.assemble(E, U1, 1.0);
961 myAssembler.assemble(E, U2, -g.mat[index - 1][0]);
962 myAssembler.assemble(E, V2, -g.mat[index - 1][1]);
963 myAssembler.assemble(U1, E, 1.0);
964 myAssembler.assemble(U2, E, -g.mat[index - 1][0]);
965 myAssembler.assemble(V2, E, -g.mat[index - 1][1]);
966 }
967 }
968
computePotential(GModel * gm,std::vector<GFace * > & f,dofManager<double> & dof,std::map<MEdge,cross2d,MEdgeLessThan> & C,std::map<MVertex *,MVertex *,MVertexPtrLessThan> & new2old,std::vector<std::vector<cross2d * >> & groups,std::map<MEdge,MEdge,MEdgeLessThan> & duplicateEdges,std::map<MTriangle *,SVector3> & lift,std::map<MTriangle *,SVector3> & lift2,std::vector<groupOfCross2d> & G,std::map<MVertex *,double> & res,std::map<MVertex *,double> & res2,std::vector<cutGraphPassage> & passages)969 static void computePotential(
970 GModel *gm, std::vector<GFace *> &f, dofManager<double> &dof,
971 std::map<MEdge, cross2d, MEdgeLessThan> &C,
972 std::map<MVertex *, MVertex *, MVertexPtrLessThan> &new2old,
973 std::vector<std::vector<cross2d *> > &groups,
974 std::map<MEdge, MEdge, MEdgeLessThan> &duplicateEdges,
975 std::map<MTriangle *, SVector3> &lift, std::map<MTriangle *, SVector3> &lift2,
976 std::vector<groupOfCross2d> &G, std::map<MVertex *, double> &res,
977 std::map<MVertex *, double> &res2, std::vector<cutGraphPassage> &passages)
978 {
979 double a[3];
980 std::set<MVertex *, MVertexPtrLessThan> vs;
981 for(size_t i = 0; i < f.size(); i++) {
982 for(size_t j = 0; j < f[i]->triangles.size(); j++) {
983 MTriangle *t = f[i]->triangles[j];
984 for(size_t k = 0; k < 3; k++) { vs.insert(t->getVertex(k)); }
985 }
986 }
987
988 #if defined(HAVE_PETSC)
989 linearSystemPETSc<double> *_lsys = new linearSystemPETSc<double>;
990 #elif defined(HAVE_GMM)
991 // linearSystemFull<double> *_lsys = new linearSystemFull<double>;
992 linearSystemGmm<double> *_lsys = new linearSystemGmm<double>;
993 /// _lsys->setPrec(2.e-7);
994 #else
995 linearSystemFull<double> *_lsys = new linearSystemFull<double>;
996 #endif
997
998 dofManager<double> myAssembler(_lsys);
999
1000 // int NUMDOF = dir+1;
1001
1002 createDofs(myAssembler, 1, vs);
1003 createDofs(myAssembler, 2, vs);
1004 LagrangeMultipliers2(myAssembler, 1, C, groups, duplicateEdges, false, lift);
1005 LagrangeMultipliers2(myAssembler, 2, C, groups, duplicateEdges, false, lift2);
1006
1007 createExtraConnexions(myAssembler, G, passages);
1008
1009 for(size_t i = 0; i < G.size(); i++) {
1010 createLagrangeMultipliers(myAssembler, G[i]);
1011 LagrangeMultipliers3(myAssembler, G[i], lift, false);
1012 }
1013
1014 LagrangeMultipliers2(myAssembler, 1, C, groups, duplicateEdges, true, lift);
1015 LagrangeMultipliers2(myAssembler, 2, C, groups, duplicateEdges, true, lift2);
1016 for(size_t i = 0; i < G.size(); i++) {
1017 assembleLagrangeMultipliers(myAssembler, G[i]);
1018 LagrangeMultipliers3(myAssembler, G[i], lift, true);
1019 }
1020
1021 assembleExtraConnexions(myAssembler, G, passages);
1022
1023 simpleFunction<double> ONE(1.0);
1024 laplaceTerm l(nullptr, 1, &ONE);
1025 laplaceTerm l2(nullptr, 2, &ONE);
1026
1027 for(size_t i = 0; i < f.size(); i++) {
1028 for(size_t j = 0; j < f[i]->triangles.size(); j++) {
1029 MTriangle *t = f[i]->triangles[j];
1030 SElement se(t);
1031 l.addToMatrix(myAssembler, &se);
1032 l2.addToMatrix(myAssembler, &se);
1033 SVector3 a0 = lift[t];
1034 SVector3 a1 = lift2[t];
1035 double va, vb, vc;
1036 auto itx = new2old.find(t->getVertex(0));
1037 dof.getDofValue(itx == new2old.end() ? t->getVertex(0) : itx->second, 0,
1038 1, va);
1039 itx = new2old.find(t->getVertex(1));
1040 dof.getDofValue(itx == new2old.end() ? t->getVertex(1) : itx->second, 0,
1041 1, vb);
1042 itx = new2old.find(t->getVertex(2));
1043 dof.getDofValue(itx == new2old.end() ? t->getVertex(2) : itx->second, 0,
1044 1, vc);
1045
1046 double F = (exp(va) + exp(vb) + exp(vc)) / 3.0;
1047
1048 a0 *= F;
1049 a1 *= F;
1050
1051 SPoint3 pp = t->barycenter();
1052 double G1[3] = {a0.x(), a0.y(), a0.z()};
1053 double G2[3] = {a0.x(), a0.y(), a0.z()};
1054 double G3[3] = {a0.x(), a0.y(), a0.z()};
1055 double G11[3] = {a1.x(), a1.y(), a1.z()};
1056 double G21[3] = {a1.x(), a1.y(), a1.z()};
1057 double G31[3] = {a1.x(), a1.y(), a1.z()};
1058 double g1[3];
1059 a[0] = 1;
1060 a[1] = 0;
1061 a[2] = 0;
1062 t->interpolateGrad(a, 0, 0, 0, g1);
1063 double RHS1 = g1[0] * G1[0] + g1[1] * G1[1] + g1[2] * G1[2];
1064 double RHS11 = g1[0] * G11[0] + g1[1] * G11[1] + g1[2] * G11[2];
1065 a[0] = 0;
1066 a[1] = 1;
1067 a[2] = 0;
1068 t->interpolateGrad(a, 0, 0, 0, g1);
1069 double RHS2 = g1[0] * G2[0] + g1[1] * G2[1] + g1[2] * G2[2];
1070 double RHS21 = g1[0] * G21[0] + g1[1] * G21[1] + g1[2] * G21[2];
1071 a[0] = 0;
1072 a[1] = 0;
1073 a[2] = 1;
1074 t->interpolateGrad(a, 0, 0, 0, g1);
1075 double RHS3 = g1[0] * G3[0] + g1[1] * G3[1] + g1[2] * G3[2];
1076 double RHS31 = g1[0] * G31[0] + g1[1] * G31[1] + g1[2] * G31[2];
1077 int num1 = myAssembler.getDofNumber(l.getLocalDofR(&se, 0));
1078 int num2 = myAssembler.getDofNumber(l.getLocalDofR(&se, 1));
1079 int num3 = myAssembler.getDofNumber(l.getLocalDofR(&se, 2));
1080 int num11 = myAssembler.getDofNumber(l2.getLocalDofR(&se, 0));
1081 int num21 = myAssembler.getDofNumber(l2.getLocalDofR(&se, 1));
1082 int num31 = myAssembler.getDofNumber(l2.getLocalDofR(&se, 2));
1083
1084 double V = t->getVolume();
1085 _lsys->addToRightHandSide(num1, RHS1 * V);
1086 _lsys->addToRightHandSide(num2, RHS2 * V);
1087 _lsys->addToRightHandSide(num3, RHS3 * V);
1088 _lsys->addToRightHandSide(num11, RHS11 * V);
1089 _lsys->addToRightHandSide(num21, RHS21 * V);
1090 _lsys->addToRightHandSide(num31, RHS31 * V);
1091 }
1092 }
1093 double A = TimeOfDay();
1094 _lsys->systemSolve();
1095 double B = TimeOfDay();
1096 Msg::Info("Computing potentials (%d unknowns) in %3lf seconds",
1097 myAssembler.sizeOfR(), B - A);
1098
1099 FILE *F = fopen("map.pos", "w");
1100 FILE *F2 = fopen("mapstr.pos", "w");
1101 fprintf(F, "View \"MAP\"{\n");
1102 fprintf(F2, "View \"MAPSTR\"{\n");
1103 for(size_t i = 0; i < f.size(); i++) {
1104 for(size_t j = 0; j < f[i]->triangles.size(); j++) {
1105 MTriangle *t = f[i]->triangles[j];
1106 double a, b, c;
1107 double a1, b1, c1;
1108 myAssembler.getDofValue(t->getVertex(0), 0, 1, a);
1109 myAssembler.getDofValue(t->getVertex(1), 0, 1, b);
1110 myAssembler.getDofValue(t->getVertex(2), 0, 1, c);
1111 myAssembler.getDofValue(t->getVertex(0), 0, 2, a1);
1112 myAssembler.getDofValue(t->getVertex(1), 0, 2, b1);
1113 myAssembler.getDofValue(t->getVertex(2), 0, 2, c1);
1114 fprintf(F, "ST(%g,%g,%g,%g,%g,%g,%g,%g,%g){%g,%g,%g,%g,%g,%g};\n", a, a1,
1115 0., b, b1, 0., c, c1, 0., a, b, c, a1, b1, c1);
1116 fprintf(F2, "ST(%g,%g,%g,%g,%g,%g,%g,%g,%g){%g,%g,%g,%g,%g,%g};\n",
1117 .2 * a + t->getVertex(0)->x(), -.2 * a1 + t->getVertex(0)->y(),
1118 0., .2 * b + t->getVertex(1)->x(),
1119 -.2 * b1 + t->getVertex(1)->y(), 0.,
1120 .2 * c + t->getVertex(2)->x(), -.2 * c1 + t->getVertex(2)->y(),
1121 0., a, b, c, a1, b1, c1);
1122
1123 res[t->getVertex(0)] = a;
1124 res[t->getVertex(1)] = b;
1125 res[t->getVertex(2)] = c;
1126 res2[t->getVertex(0)] = a1;
1127 res2[t->getVertex(1)] = b1;
1128 res2[t->getVertex(2)] = c1;
1129 }
1130 }
1131 fprintf(F, "};\n");
1132 fclose(F);
1133 fprintf(F2, "};\n");
1134 fclose(F2);
1135 }
1136
distance(MTriangle * t,std::set<MVertex *,MVertexPtrLessThan> & boundaries)1137 static double distance(MTriangle *t,
1138 std::set<MVertex *, MVertexPtrLessThan> &boundaries)
1139 {
1140 // return drand48();
1141 SPoint3 p = t->barycenter();
1142 double dmin = 1.e22;
1143 for(auto it = boundaries.begin(); it != boundaries.end(); ++it) {
1144 SPoint3 pp((*it)->x(), (*it)->y(), (*it)->z());
1145 double d = p.distance(pp);
1146 if(d < dmin) { dmin = d; }
1147 }
1148 return -dmin;
1149 }
1150
1151 struct temp_comp {
1152 cross2d *cr;
1153 double a;
temp_comptemp_comp1154 temp_comp(MVertex *v, cross2d *c, cross2d *ref, SVector3 &n) : cr(c)
1155 {
1156 MVertex *tref =
1157 ref->_e.getVertex(0) == v ? ref->_e.getVertex(1) : ref->_e.getVertex(0);
1158 MVertex *tc =
1159 c->_e.getVertex(0) == v ? c->_e.getVertex(1) : c->_e.getVertex(0);
1160
1161 SVector3 t1(tref->x() - v->x(), tref->y() - v->y(), tref->z() - v->z());
1162 SVector3 t2(tc->x() - v->x(), tc->y() - v->y(), tc->z() - v->z());
1163 t1.normalize();
1164 t2.normalize();
1165
1166 double cosTheta = dot(t1, t2);
1167 double sinTheta;
1168 SVector3 cc = crossprod(t1, t2);
1169 if(dot(cc, n) > 0)
1170 sinTheta = norm(crossprod(t1, t2));
1171 else
1172 sinTheta = -norm(crossprod(t1, t2));
1173 a = atan2(sinTheta, cosTheta);
1174 }
operator <temp_comp1175 bool operator<(const temp_comp &other) const { return a < other.a; }
1176 };
1177
isSingular(MVertex * v,std::vector<cross2d * > & adj,double & MAX)1178 static bool isSingular(MVertex *v, std::vector<cross2d *> &adj, double &MAX)
1179 {
1180 const std::size_t TEST = 0;
1181 if(v->getNum() == TEST) printf("VERTEX %lu\n", v->getNum());
1182 SVector3 n(0, 0, 0);
1183 for(size_t i = 0; i < adj.size(); i++) {
1184 if(adj[i]->inBoundary || adj[i]->inInternalBoundary) return false;
1185 n += adj[i]->_nrml;
1186 }
1187 n.normalize();
1188
1189 std::vector<temp_comp> cc;
1190 for(size_t i = 0; i < adj.size(); i++) {
1191 cc.push_back(temp_comp(v, adj[i], adj[0], n));
1192 }
1193 std::sort(cc.begin(), cc.end());
1194 SVector3 ref = cc[0].cr->_tgt * cos(cc[0].cr->_atemp) +
1195 cc[0].cr->_tgt2 * sin(cc[0].cr->_atemp);
1196 SVector3 ref0 = ref;
1197 for(size_t i = 1; i < cc.size() + 1; i++) {
1198 cross2d &c2 = *(cc[i % cc.size()].cr);
1199 double P = M_PI / 2;
1200 SVector3 d1 = c2._tgt * cos(c2._atemp) + c2._tgt2 * sin(c2._atemp);
1201 SVector3 d2 = c2._tgt * cos(c2._atemp + P) + c2._tgt2 * sin(c2._atemp + P);
1202 SVector3 d3 =
1203 c2._tgt * cos(c2._atemp + 2 * P) + c2._tgt2 * sin(c2._atemp + 2 * P);
1204 SVector3 d4 =
1205 c2._tgt * cos(c2._atemp + 3 * P) + c2._tgt2 * sin(c2._atemp + 3 * P);
1206 double D1 = dot(ref, d1);
1207 double D2 = dot(ref, d2);
1208 double D3 = dot(ref, d3);
1209 double D4 = dot(ref, d4);
1210 if(D1 > D2 && D1 > D3 && D1 > D4)
1211 ref = d1;
1212 else if(D2 > D1 && D2 > D3 && D2 > D4)
1213 ref = d2;
1214 else if(D3 > D1 && D3 > D2 && D3 > D4)
1215 ref = d3;
1216 else
1217 ref = d4;
1218 }
1219
1220 if(v->getNum() == TEST)
1221 printf("VERTEX %lu %12.5E %12.5E %12.5E\n", v->getNum(), n.x(), n.y(),
1222 n.z());
1223 SVector3 t0, b0;
1224 std::vector<double> angles;
1225 for(size_t i = 0; i < adj.size(); i++) {
1226 if(i == 0) {
1227 SVector3 t =
1228 (adj[i]->_e.getVertex(0) == v) ? -adj[i]->_tgt : adj[i]->_tgt;
1229 t -= n * (dot(t, n));
1230 t.normalize();
1231 SVector3 b = crossprod(n, t);
1232 b0 = b;
1233 t0 = t;
1234 }
1235
1236 SVector3 repr = adj[i]->o_i - n * dot(adj[i]->o_i, n);
1237 repr.normalize();
1238 // t * dot (,adj[i]->_tgt) +
1239 // b * dot (adj[i]->o_i,adj[i]->_tgt2) ;
1240 double angle = atan2(dot(repr, t0), dot(repr, b0));
1241 adj[i]->normalize(angle);
1242 angles.push_back(angle);
1243 if(v->getNum() == TEST) {
1244 printf("EDGE %lu %lu\n", adj[i]->_e.getVertex(0)->getNum(),
1245 adj[i]->_e.getVertex(1)->getNum());
1246 printf("o %12.5E %12.5E %12.5E\n", adj[i]->o_i.x(), adj[i]->o_i.y(),
1247 adj[i]->o_i.z());
1248 printf("ANGLE = %12.5E %12.5E\n", angle * 180 / M_PI,
1249 lifting(angles[0], angles[i]) * 180 / M_PI);
1250 }
1251 }
1252
1253 MAX = 0;
1254 for(size_t i = 0; i < angles.size(); i++) {
1255 if(v->getNum() == TEST) printf("%12.5E ", angles[i]);
1256 for(size_t j = 0; j < i; j++) {
1257 MAX = std::max(MAX, fabs(angles[i] - lifting(angles[j], angles[i])));
1258 }
1259 }
1260 if(v->getNum() == TEST) printf("\n");
1261 if(v->getNum() == TEST)
1262 printf("vertex %lu %lu edges %12.5E\n", v->getNum(), adj.size(), MAX);
1263 // if (MAX > .5)printf("vertex %lu %lu edges %12.5E -- new method %12.5E\n",
1264 // v->getNum(), adj.size(), MAX,dot(ref,ref0));
1265 return MAX > .5;
1266 }
1267
isMinMax(MVertex * v,std::vector<cross2d * > adj,dofManager<double> * dof,bool & isMin,bool & isMax)1268 void isMinMax(MVertex *v, std::vector<cross2d *> adj, dofManager<double> *dof,
1269 bool &isMin, bool &isMax)
1270 {
1271 double aa;
1272 isMin = isMax = true;
1273 dof->getDofValue(v, 0, 1, aa);
1274 for(size_t i = 0; i < adj.size(); i++) {
1275 double a;
1276 dof->getDofValue(adj[i]->_e.getVertex(0) == v ? adj[i]->_e.getVertex(1) :
1277 adj[i]->_e.getVertex(0),
1278 0, 1, a);
1279 if(a < aa) isMin = false;
1280 if(a > aa) isMax = false;
1281 }
1282 }
1283
1284 static void
computeSingularities(std::map<MEdge,cross2d,MEdgeLessThan> & C,std::set<MVertex *,MVertexPtrLessThan> & singularities,std::map<MVertex *,int> & indices,dofManager<double> * dof)1285 computeSingularities(std::map<MEdge, cross2d, MEdgeLessThan> &C,
1286 std::set<MVertex *, MVertexPtrLessThan> &singularities,
1287 std::map<MVertex *, int> &indices, dofManager<double> *dof)
1288 {
1289 FILE *f_ = fopen("sing.pos", "w");
1290 fprintf(f_, "View \"S\"{\n");
1291 std::multimap<MVertex *, cross2d *, MVertexPtrLessThan> conn;
1292 for(auto it = C.begin(); it != C.end(); ++it) {
1293 std::pair<MVertex *, cross2d *> p =
1294 std::make_pair(it->first.getVertex(0), &it->second);
1295 conn.insert(p);
1296 p = std::make_pair(it->first.getVertex(1), &it->second);
1297 conn.insert(p);
1298 }
1299 MVertex *v = nullptr;
1300 std::vector<cross2d *> adj;
1301 for(auto it = conn.begin(); it != conn.end(); ++it) {
1302 if(it->first == v) { adj.push_back(it->second); }
1303 else {
1304 double MAX;
1305 if(v && isSingular(v, adj, MAX)) {
1306 singularities.insert(v);
1307 bool isMin, isMax;
1308 isMinMax(v, adj, dof, isMin, isMax);
1309 if(isMax)
1310 indices[v] = 1;
1311 else if(isMin)
1312 indices[v] = -1;
1313 else
1314 printf("ERROR -- \n");
1315 fprintf(f_, "SP(%g,%g,%g){%d};\n", v->x(), v->y(), v->z(), indices[v]);
1316 }
1317 adj.clear();
1318 v = it->first;
1319 adj.push_back(it->second);
1320 }
1321 }
1322 fprintf(f_, "};\n");
1323 fclose(f_);
1324 }
1325
cutGraph(std::map<MEdge,cross2d,MEdgeLessThan> & C,std::set<MEdge,MEdgeLessThan> & cutG,std::set<MVertex *,MVertexPtrLessThan> & singularities,std::set<MVertex *,MVertexPtrLessThan> & boundaries)1326 static void cutGraph(std::map<MEdge, cross2d, MEdgeLessThan> &C,
1327 std::set<MEdge, MEdgeLessThan> &cutG,
1328 std::set<MVertex *, MVertexPtrLessThan> &singularities,
1329 std::set<MVertex *, MVertexPtrLessThan> &boundaries)
1330 {
1331 std::set<MTriangle *, MElementPtrLessThan> touched;
1332 std::vector<cross2d *> tree;
1333 std::vector<MEdge> cotree;
1334 std::set<std::pair<double, MTriangle *> > _distances;
1335 {
1336 auto it = C.begin();
1337 for(; it != C.end(); ++it) {
1338 if(it->second._t.size() == 1) {
1339 boundaries.insert(it->first.getVertex(0));
1340 boundaries.insert(it->first.getVertex(1));
1341 }
1342 }
1343 }
1344 std::set<MVertex *, MVertexPtrLessThan> _all = boundaries;
1345 _all.insert(singularities.begin(), singularities.end());
1346
1347 FILE *fff2 = fopen("tree.pos", "w");
1348 fprintf(fff2, "View \"sides\"{\n");
1349
1350 MTriangle *t = (C.begin())->second._t[0];
1351 std::pair<double, MTriangle *> pp = std::make_pair(distance(t, _all), t);
1352 _distances.insert(pp);
1353 touched.insert(t);
1354 while(!_distances.empty()) {
1355 t = (_distances.begin()->second);
1356 _distances.erase(_distances.begin());
1357
1358 for(int i = 0; i < 3; i++) {
1359 MEdge e = t->getEdge(i);
1360 auto it = C.find(e);
1361 for(size_t j = 0; j < it->second._t.size(); j++) {
1362 MTriangle *tt = it->second._t[j];
1363 if(touched.find(tt) == touched.end()) {
1364 std::pair<double, MTriangle *> pp =
1365 std::make_pair(distance(t, _all), tt);
1366 _distances.insert(pp);
1367 touched.insert(tt);
1368 tree.push_back(&it->second);
1369 }
1370 }
1371 }
1372 }
1373
1374 std::sort(tree.begin(), tree.end());
1375 auto it = C.begin();
1376 std::map<MVertex *, std::vector<MEdge>, MVertexPtrLessThan> _graph;
1377 for(; it != C.end(); ++it) {
1378 if(!std::binary_search(tree.begin(), tree.end(), &it->second)) {
1379 for(int i = 0; i < 2; i++) {
1380 auto it0 = _graph.find(it->first.getVertex(i));
1381 if(it0 == _graph.end()) {
1382 std::vector<MEdge> ee;
1383 ee.push_back(it->first);
1384 _graph[it->first.getVertex(i)] = ee;
1385 }
1386 else
1387 it0->second.push_back(it->first);
1388 }
1389 cotree.push_back(it->first);
1390 fprintf(fff2, "SL(%g,%g,%g,%g,%g,%g){%d,%d};\n",
1391 it->first.getVertex(0)->x(), it->first.getVertex(0)->y(),
1392 it->first.getVertex(0)->z(), it->first.getVertex(1)->x(),
1393 it->first.getVertex(1)->y(), it->first.getVertex(1)->z(), 1, 1);
1394 }
1395 }
1396
1397 fprintf(fff2, "};\n");
1398 fclose(fff2);
1399
1400 // {
1401 // for(auto it = boundaries.begin(); it != boundaries.end(); ++it) {
1402 // auto it2 = singularities.find(*it); if(it2 != singularities.end())
1403 // singularities.erase(it2);
1404 // }
1405 // }
1406
1407 while(1) {
1408 bool somethingDone = false;
1409 auto it = _graph.begin();
1410 for(; it != _graph.end(); ++it) {
1411 if(it->second.size() == 1) {
1412 MVertex *v1 = it->second[0].getVertex(0);
1413 MVertex *v2 = it->second[0].getVertex(1);
1414 if(boundaries.find(it->first) == boundaries.end() &&
1415 singularities.find(it->first) == singularities.end()) {
1416 somethingDone = true;
1417 auto it2 = _graph.find(v1 == it->first ? v2 : v1);
1418 auto position =
1419 std::find(it2->second.begin(), it2->second.end(), it->second[0]);
1420 it2->second.erase(position);
1421 it->second.clear();
1422 }
1423 }
1424 }
1425 if(!somethingDone) break;
1426 }
1427
1428 FILE *fff = fopen("cotree.pos", "w");
1429 fprintf(fff, "View \"sides\"{\n");
1430 {
1431 auto it = _graph.begin();
1432 for(; it != _graph.end(); ++it) {
1433 for(size_t i = 0; i < it->second.size(); i++) {
1434 MEdge e = it->second[i];
1435 if(boundaries.find(e.getVertex(0)) == boundaries.end() ||
1436 boundaries.find(e.getVertex(1)) == boundaries.end()) {
1437 cutG.insert(e);
1438 }
1439 }
1440 }
1441 }
1442
1443 // Add internal boundaries to the cut graph
1444 {
1445 auto it = C.begin();
1446 for(; it != C.end(); ++it) {
1447 if(it->second._t.size() > 1 && it->second.inInternalBoundary) {
1448 cutG.insert(it->second._e);
1449 }
1450 }
1451 }
1452
1453 {
1454 auto it = cutG.begin();
1455 for(; it != cutG.end(); ++it) {
1456 MEdge e = *it;
1457 fprintf(fff, "SL(%g,%g,%g,%g,%g,%g){%d,%d};\n", e.getVertex(0)->x(),
1458 e.getVertex(0)->y(), e.getVertex(0)->z(), e.getVertex(1)->x(),
1459 e.getVertex(1)->y(), e.getVertex(1)->z(), 1, 1);
1460 }
1461 }
1462 fprintf(fff, "};\n");
1463 fclose(fff);
1464 }
1465
analyzeCorner(std::multimap<MVertex *,cross2d * > & conn,MVertex * v)1466 int analyzeCorner(std::multimap<MVertex *, cross2d *> &conn, MVertex *v)
1467 {
1468 // compute if this is an external (1) or internal (-1) corner.
1469 return 1;
1470 }
1471
1472 static void
groupBoundaries(GModel * gm,std::map<MEdge,cross2d,MEdgeLessThan> & C,std::vector<std::vector<cross2d * >> & groups,std::set<MVertex *,MVertexPtrLessThan> singularities,std::set<MVertex *,MVertexPtrLessThan> & corners,bool cutGraph=false)1473 groupBoundaries(GModel *gm, std::map<MEdge, cross2d, MEdgeLessThan> &C,
1474 std::vector<std::vector<cross2d *> > &groups,
1475 std::set<MVertex *, MVertexPtrLessThan> singularities,
1476 std::set<MVertex *, MVertexPtrLessThan> &corners,
1477 bool cutGraph = false)
1478 {
1479 std::set<MVertex *, MVertexPtrLessThan> cutgraph;
1480 std::set<MVertex *, MVertexPtrLessThan> boundaries;
1481 for(auto it = C.begin(); it != C.end(); ++it) {
1482 MVertex *v0 = it->first.getVertex(0);
1483 MVertex *v1 = it->first.getVertex(1);
1484 if(it->second.inBoundary) {
1485 boundaries.insert(v0);
1486 boundaries.insert(v1);
1487 }
1488 else if(it->second.inCutGraph) {
1489 cutgraph.insert(v0);
1490 cutgraph.insert(v1);
1491 }
1492 }
1493
1494 std::set<cross2d *> _all;
1495
1496 std::multimap<MVertex *, cross2d *> conn;
1497 for(auto it = C.begin(); it != C.end(); ++it) {
1498 std::pair<MVertex *, cross2d *> p =
1499 std::make_pair(it->first.getVertex(0), &it->second);
1500 conn.insert(p);
1501 p = std::make_pair(it->first.getVertex(1), &it->second);
1502 conn.insert(p);
1503 }
1504
1505 for(auto it = boundaries.begin(); it != boundaries.end(); ++it) {
1506 MVertex *v = *it;
1507 std::vector<cross2d *> bnd;
1508 int countCutGraph = 0;
1509 for(auto it2 = conn.lower_bound(v); it2 != conn.upper_bound(v); ++it2) {
1510 if(it2->second->inBoundary) { bnd.push_back(it2->second); }
1511 if(it2->second->inCutGraph) { countCutGraph++; }
1512 }
1513 if(bnd.size() == 2) {
1514 if(fabs(dot(bnd[0]->o_i, bnd[1]->o_i)) < .25) {
1515 corners.insert(v);
1516 cutgraph.insert(v);
1517 }
1518 if(countCutGraph == 1) { singularities.insert(v); }
1519 }
1520 if(bnd.size() > 2) { cutgraph.insert(v); }
1521 }
1522
1523 std::string ss = gm->getName();
1524 std::string fn = cutGraph ? ss + "_groups_cg.pos" : ss + "_groups_bnd.pos";
1525
1526 FILE *f = fopen(fn.c_str(), "w");
1527
1528 fprintf(f, "View \" \"{\n");
1529
1530 std::set<MVertex *, MVertexPtrLessThan> endPoints = singularities;
1531 {
1532 for(auto it = conn.begin(); it != conn.end(); ++it) {
1533 int count = 0;
1534 for(auto it2 = conn.lower_bound(it->first);
1535 it2 != conn.upper_bound(it->first); ++it2) {
1536 if(it2->second->inCutGraph) { count++; }
1537 }
1538 if(count > 2) endPoints.insert(it->first);
1539 }
1540 }
1541
1542 for(int AA = 0; AA < 4; AA++) {
1543 if(cutGraph) {
1544 for(auto it = endPoints.begin(); it != endPoints.end(); ++it) {
1545 MVertex *v = *it;
1546 std::vector<cross2d *> group;
1547 do {
1548 MVertex *vnew = nullptr;
1549 for(auto it2 = conn.lower_bound(v); it2 != conn.upper_bound(v);
1550 ++it2) {
1551 if((_all.find(it2->second) == _all.end()) &&
1552 (group.empty() || group[group.size() - 1] != it2->second) &&
1553 it2->second->inCutGraph) {
1554 group.push_back(it2->second);
1555 vnew = (it2->second->_e.getVertex(0) == v) ?
1556 it2->second->_e.getVertex(1) :
1557 it2->second->_e.getVertex(0);
1558 fprintf(f, "SL(%g,%g,%g,%g,%g,%g){%lu,%lu};\n",
1559 it2->second->_e.getVertex(0)->x(),
1560 it2->second->_e.getVertex(0)->y(),
1561 it2->second->_e.getVertex(0)->z(),
1562 it2->second->_e.getVertex(1)->x(),
1563 it2->second->_e.getVertex(1)->y(),
1564 it2->second->_e.getVertex(1)->z(), groups.size(),
1565 groups.size());
1566 break;
1567 }
1568 }
1569 if(vnew == nullptr) break;
1570 v = vnew;
1571 } while((boundaries.find(v) == boundaries.end()) &&
1572 (endPoints.find(v) == endPoints.end()));
1573 if(group.size()) {
1574 groups.push_back(group);
1575 _all.insert(group.begin(), group.end());
1576 }
1577 }
1578 }
1579 else {
1580 for(auto it = boundaries.begin(); it != boundaries.end(); ++it) {
1581 MVertex *v = *it;
1582 if(cutgraph.find(v) != cutgraph.end() ||
1583 singularities.find(v) != singularities.end()) {
1584 // printf("START POINT %lu %d %d\n",v->getNum(),cutgraph.find(v)
1585 //!= cutgraph.end() , singularities.find(v) !=
1586 //! singularities.end());
1587 std::vector<cross2d *> group;
1588 do {
1589 MVertex *vnew = nullptr;
1590 for(auto it2 = conn.lower_bound(v); it2 != conn.upper_bound(v);
1591 ++it2) {
1592 if((_all.find(it2->second) == _all.end()) &&
1593 (group.empty() || group[group.size() - 1] != it2->second) &&
1594 (it2->second->inBoundary)) {
1595 group.push_back(it2->second);
1596 vnew = (it2->second->_e.getVertex(0) == v) ?
1597 it2->second->_e.getVertex(1) :
1598 it2->second->_e.getVertex(0);
1599 // printf("EDGE %lu %lu (%d)\n",
1600 // it2->second->_e.getVertex(0)->getNum(),it2->second->_e.getVertex(1)->getNum(),
1601 // singularities.find(v) == singularities.end());
1602 // printf("v %lu EDGE %lu
1603 //%lu\n",v->getNum(),it2->second->_e.getVertex(0)->getNum(),it2->second->_e.getVertex(1)->getNum());
1604 fprintf(f, "SL(%g,%g,%g,%g,%g,%g){%lu,%lu};\n",
1605 it2->second->_e.getVertex(0)->x(),
1606 it2->second->_e.getVertex(0)->y(),
1607 it2->second->_e.getVertex(0)->z(),
1608 it2->second->_e.getVertex(1)->x(),
1609 it2->second->_e.getVertex(1)->y(),
1610 it2->second->_e.getVertex(1)->z(), groups.size(),
1611 groups.size());
1612 break;
1613 }
1614 }
1615 if(vnew == nullptr) break;
1616 v = vnew;
1617 // printf("NEXT POINT %lu %d
1618 //%lu\n",v->getNum(),singularities.find(v) == singularities.end(),
1619 // singularities.size());
1620 } while(cutgraph.find(v) == cutgraph.end() &&
1621 singularities.find(v) == singularities.end());
1622 if(group.size() && _all.find(group[0]) == _all.end()) {
1623 groups.push_back(group);
1624 _all.insert(group.begin(), group.end());
1625 }
1626 }
1627 }
1628 }
1629 }
1630 fprintf(f, "};\n");
1631 fclose(f);
1632 }
1633
1634 static void
fastImplementationExtrinsic(std::map<MEdge,cross2d,MEdgeLessThan> & C,double tol=1.e-10)1635 fastImplementationExtrinsic(std::map<MEdge, cross2d, MEdgeLessThan> &C,
1636 double tol = 1.e-10)
1637 {
1638 double *data = new double[C.size() * 6];
1639 size_t *graph = new size_t[C.size() * 4];
1640 auto it = C.begin();
1641 int counter = 0;
1642
1643 for(; it != C.end(); ++it) {
1644 data[6 * counter + 0] = it->second.o_i.x();
1645 data[6 * counter + 1] = it->second.o_i.y();
1646 data[6 * counter + 2] = it->second.o_i.z();
1647 data[6 * counter + 3] = it->second._nrml.x();
1648 data[6 * counter + 4] = it->second._nrml.y();
1649 data[6 * counter + 5] = it->second._nrml.z();
1650 it->second.counter = counter;
1651 ++counter;
1652 }
1653
1654 it = C.begin();
1655 counter = 0;
1656 for(; it != C.end(); ++it) {
1657 graph[4 * counter + 0] = graph[4 * counter + 1] = graph[4 * counter + 2] =
1658 graph[4 * counter + 3] = it->second.counter;
1659 for(size_t i = 0; i < it->second._cneighbors.size(); i++) {
1660 graph[4 * counter + i] = it->second._cneighbors[i]->counter;
1661 }
1662 if(it->second.inBoundary || it->second.inInternalBoundary) {
1663 graph[4 * counter + 2] = graph[4 * counter + 3] = it->second.counter;
1664 }
1665
1666 counter++;
1667 }
1668
1669 size_t N = C.size();
1670 int MAXITER = 10000;
1671 int ITER = -1;
1672 while(ITER++ < MAXITER) {
1673 double x[3], y[3];
1674 double RES = 0;
1675 for(size_t i = 0; i < N; i++) {
1676 double *r = &data[6 * i + 0];
1677 double *n = &data[6 * i + 3];
1678 SVector3 ro(r[0], r[1], r[2]);
1679 size_t *neigh = &graph[4 * i];
1680 double weight = 0;
1681 if(neigh[2] != neigh[3]) {
1682 for(int j = 0; j < 4; j++) {
1683 size_t k = neigh[j];
1684 const double *r2 = &data[6 * k + 0];
1685 const double *n2 = &data[6 * k + 3];
1686 compat_orientation_extrinsic(r, n, r2, n2, x, y);
1687 r[0] = x[0] * weight + y[0];
1688 r[1] = x[1] * weight + y[1];
1689 r[2] = x[2] * weight + y[2];
1690 const double dd = r[0] * n[0] + r[1] * n[1] + r[2] * n[2];
1691 r[0] -= n[0] * dd;
1692 r[1] -= n[1] * dd;
1693 r[2] -= n[2] * dd;
1694 double NRM = sqrt(r[0] * r[0] + r[1] * r[1] + r[2] * r[2]);
1695 if(NRM != 0.0) {
1696 r[0] /= NRM;
1697 r[1] /= NRM;
1698 r[2] /= NRM;
1699 }
1700 weight += 1;
1701 }
1702 double dp = r[0] * ro[0] + r[1] * ro[1] + r[2] * ro[2];
1703 RES += std::min(1. - fabs(dp), fabs(dp));
1704 }
1705 // data[6*i+0]=r[0];
1706 // data[6*i+1]=r[1];
1707 // data[6*i+2]=r[2];
1708 }
1709 if(ITER % 1000 == 0)
1710 Msg::Info("NL smooth : iter %6d RES = %12.5E", ITER, RES);
1711 if(RES < tol) break;
1712 }
1713
1714 it = C.begin();
1715 for(; it != C.end(); ++it) {
1716 counter = it->second.counter;
1717 it->second.o_i[0] = data[6 * counter + 0];
1718 it->second.o_i[1] = data[6 * counter + 1];
1719 it->second.o_i[2] = data[6 * counter + 2];
1720 }
1721 delete[] data;
1722 delete[] graph;
1723 }
1724
computeH(GModel * gm,std::vector<GFace * > & f,std::set<MVertex *,MVertexPtrLessThan> & vs,std::map<MEdge,cross2d,MEdgeLessThan> & C)1725 static dofManager<double> *computeH(GModel *gm, std::vector<GFace *> &f,
1726 std::set<MVertex *, MVertexPtrLessThan> &vs,
1727 std::map<MEdge, cross2d, MEdgeLessThan> &C)
1728 {
1729 #if defined(HAVE_PETSC)
1730 linearSystemPETSc<double> *_lsys = new linearSystemPETSc<double>;
1731 #elif defined(HAVE_GMM)
1732 // MUMPS !!!
1733 linearSystemGmm<double> *_lsys = new linearSystemGmm<double>;
1734 #else
1735 linearSystemFull<double> *_lsys = new linearSystemFull<double>;
1736 #endif
1737
1738 dofManager<double> *myAssembler = new dofManager<double>(_lsys);
1739
1740 // myAssembler.fixVertex(*vs.begin(), 0, 1, 0);
1741 for(auto it = vs.begin(); it != vs.end(); ++it)
1742 myAssembler->numberVertex(*it, 0, 1);
1743
1744 std::string ss = gm->getName();
1745 std::string fn = ss + "_grad.pos";
1746
1747 FILE *_f = fopen(fn.c_str(), "w");
1748 fprintf(_f, "View \"grad\"{\n");
1749
1750 simpleFunction<double> ONE(1.0);
1751 laplaceTerm l(nullptr, 1, &ONE);
1752
1753 std::map<MTriangle *, SVector3> gradients_of_theta;
1754
1755 for(size_t i = 0; i < f.size(); i++) {
1756 for(size_t j = 0; j < f[i]->triangles.size(); j++) {
1757 MTriangle *t = f[i]->triangles[j];
1758
1759 SVector3 v10(t->getVertex(1)->x() - t->getVertex(0)->x(),
1760 t->getVertex(1)->y() - t->getVertex(0)->y(),
1761 t->getVertex(1)->z() - t->getVertex(0)->z());
1762 SVector3 v20(t->getVertex(2)->x() - t->getVertex(0)->x(),
1763 t->getVertex(2)->y() - t->getVertex(0)->y(),
1764 t->getVertex(2)->z() - t->getVertex(0)->z());
1765 SVector3 xx = crossprod(v20, v10);
1766 xx.normalize();
1767
1768 SElement se(t);
1769 l.addToMatrix(*myAssembler, &se);
1770 MEdge e0 = t->getEdge(0);
1771 MEdge e1 = t->getEdge(1);
1772 MEdge e2 = t->getEdge(2);
1773 auto it0 = C.find(e0);
1774 auto it1 = C.find(e1);
1775 auto it2 = C.find(e2);
1776
1777 // printf("%g %ag %g\n",a0,a1,a2);
1778
1779 double a0 = it0->second._a;
1780 double A1 =
1781 it1->second._a + atan2(dot(it0->second._tgt2, it1->second._tgt),
1782 dot(it0->second._tgt, it1->second._tgt));
1783 it0->second.normalize(A1);
1784 double a1 = it0->second.lifting(A1);
1785 double A2 =
1786 it2->second._a + atan2(dot(it0->second._tgt2, it2->second._tgt),
1787 dot(it0->second._tgt, it2->second._tgt));
1788 it0->second.normalize(A2);
1789 double a2 = it0->second.lifting(A2);
1790
1791 SVector3 x0, x1, x2, x3;
1792 SVector3 t_i = crossprod(xx, it0->second._tgt);
1793 t_i -= xx * dot(xx, t_i);
1794 t_i.normalize();
1795 SVector3 o_i = it0->second.o_i;
1796 o_i -= xx * dot(xx, o_i);
1797 o_i.normalize();
1798 SVector3 o_1 = it1->second.o_i;
1799 o_1 -= xx * dot(xx, o_1);
1800 o_1.normalize();
1801 SVector3 o_2 = it2->second.o_i;
1802 o_2 -= xx * dot(xx, o_2);
1803 o_2.normalize();
1804 compat_orientation_extrinsic(o_i, xx, o_1, xx, x0, x1);
1805 compat_orientation_extrinsic(o_i, xx, o_2, xx, x2, x3);
1806
1807 // compat_orientation_extrinsic
1808 // (it0->second.o_i,it0->second._nrml,it1->second.o_i,it1->second._nrml,x0,x1);
1809 // compat_orientation_extrinsic
1810 // (it0->second.o_i,it0->second._nrml,it2->second.o_i,it2->second._nrml,x2,x3);
1811 // a0 = atan2(dot(it0->second._tgt2,it0->second.o_i),
1812 // dot(it0->second._tgt,it0->second.o_i));
1813 a0 = atan2(dot(t_i, o_i), dot(it0->second._tgt, o_i));
1814
1815 x0 -= xx * dot(xx, x0);
1816 x0.normalize();
1817 x1 -= xx * dot(xx, x1);
1818 x1.normalize();
1819 x2 -= xx * dot(xx, x2);
1820 x2.normalize();
1821 x3 -= xx * dot(xx, x3);
1822 x3.normalize();
1823
1824 it0->second.normalize(a0);
1825 it0->second._a = a0;
1826 // A1 = atan2(dot(it0->second._tgt2,x1),
1827 // dot(it0->second._tgt,x1));
1828 // A2 = atan2(dot(it0->second._tgt2,x3),
1829 // dot(it0->second._tgt,x3));
1830 A1 = atan2(dot(t_i, x1), dot(it0->second._tgt, x1));
1831 A2 = atan2(dot(t_i, x3), dot(it0->second._tgt, x3));
1832 it0->second.normalize(A1);
1833 a1 = it0->second.lifting(A1);
1834 it0->second.normalize(A2);
1835 a2 = it0->second.lifting(A2);
1836 // it0->second.normalize(a1);
1837 // it0->second.normalize(a2);
1838 // it0->second._a = a0 = 0;
1839 // a1 = it0->second.lifting(A1);
1840 // a2 = it0->second.lifting(A2);
1841
1842 double a[3] = {a0 + a2 - a1, a0 + a1 - a2, a1 + a2 - a0};
1843 double g[3] = {0, 0, 0};
1844 t->interpolateGrad(a, 0, 0, 0, g);
1845 gradients_of_theta[t] = SVector3(g[0], g[1], g[2]);
1846 SPoint3 pp = t->barycenter();
1847
1848 // fprintf(_f,"VP(%g,%g,%g){%g,%g,%g};\n",pp.x(),pp.y(),pp.z(),
1849 // x0.x(),x0.y(),x0.z());
1850 // fprintf(_f,"VP(%g,%g,%g){%g,%g,%g};\n",pp.x(),pp.y(),pp.z(),
1851 // x1.x(),x1.y(),x1.z());
1852 // fprintf(_f,"VP(%g,%g,%g){%g,%g,%g};\n",pp.x(),pp.y(),pp.z(),
1853 // x2.x(),x2.y(),x2.z());
1854 // fprintf(_f,"VP(%g,%g,%g){%g,%g,%g};\n",pp.x(),pp.y(),pp.z(),
1855 // x3.x(),x3.y(),x3.z());
1856
1857 fprintf(_f, "VP(%g,%g,%g){%g,%g,%g};\n", pp.x(), pp.y(), pp.z(), g[0],
1858 g[1], g[2]);
1859 // fprintf(_f, "VP(%g,%g,%g){%g,%g,%g};\n", pp.x(), pp.y(), pp.z(),
1860 // -g[1],
1861 // g[0], g[2]);
1862 // printf("A %g vs %g\n",sqrt(
1863 // g[1]*g[1]+g[0]*g[0]), 1./(sqrt((pp.x())*(pp.x())+(pp.y())*(pp.y()))));
1864
1865 SVector3 G(g[0], g[1], g[2]);
1866 SVector3 GT = crossprod(G, xx);
1867
1868 double g1[3];
1869 a[0] = 1;
1870 a[1] = 0;
1871 a[2] = 0;
1872 t->interpolateGrad(a, 0, 0, 0, g1);
1873 double RHS1 = g1[0] * GT.x() + g1[1] * GT.y() + g1[2] * GT.z();
1874 // double RHS1 = g1[0]*g[0]+g1[1]*g[1]+g1[2]*g[2];
1875 a[0] = 0;
1876 a[1] = 1;
1877 a[2] = 0;
1878 t->interpolateGrad(a, 0, 0, 0, g1);
1879 double RHS2 = g1[0] * GT.x() + g1[1] * GT.y() + g1[2] * GT.z();
1880 // double RHS2 = g1[0]*g[0]+g1[1]*g[1]+g1[2]*g[2];
1881 a[0] = 0;
1882 a[1] = 0;
1883 a[2] = 1;
1884 t->interpolateGrad(a, 0, 0, 0, g1);
1885 double RHS3 = g1[0] * GT.x() + g1[1] * GT.y() + g1[2] * GT.z();
1886 // double RHS3 = g1[0]*g[0]+g1[1]*g[1]+g1[2]*g[2];
1887 int num1 = myAssembler->getDofNumber(l.getLocalDofR(&se, 0));
1888 int num2 = myAssembler->getDofNumber(l.getLocalDofR(&se, 1));
1889 int num3 = myAssembler->getDofNumber(l.getLocalDofR(&se, 2));
1890 double V = -t->getVolume();
1891 _lsys->addToRightHandSide(num1, RHS1 * V);
1892 _lsys->addToRightHandSide(num2, RHS2 * V);
1893 _lsys->addToRightHandSide(num3, RHS3 * V);
1894 }
1895 }
1896 fprintf(_f, "};\n");
1897 fclose(_f);
1898 _lsys->systemSolve();
1899
1900 return myAssembler;
1901 }
1902
1903 /*
1904 1) find u_i and v_i of all singularities
1905 2) compute [MAX,MIN] max of u's and v's
1906 3) Divide this interval into
1907
1908 */
1909
coord1d(double a0,double a1,double a)1910 static double coord1d(double a0, double a1, double a)
1911 {
1912 if(a1 == a0) return 0.0;
1913 return (a - a0) / (a1 - a0);
1914 }
1915
1916 struct edgeCuts {
1917 std::vector<SPoint3> ps;
1918 std::vector<MVertex *> vs;
1919 std::vector<int> indexOfCuts;
1920 std::vector<int> idsOfCuts;
addedgeCuts1921 bool add(const SPoint3 &p, int ind, int id)
1922 {
1923 for(size_t i = 0; i < ps.size(); i++) {
1924 SVector3 v(ps[i], p);
1925 if(v.norm() < 1.e-10) { return false; }
1926 }
1927 ps.push_back(p);
1928 indexOfCuts.push_back(ind);
1929 idsOfCuts.push_back(id);
1930 return true;
1931 }
finishedgeCuts1932 void finish(GModel *gm, FILE *f = nullptr)
1933 {
1934 for(size_t i = 0; i < ps.size(); i++) {
1935 GEdge *ge = gm->getEdgeByTag(indexOfCuts[i]);
1936 if(!ge) {
1937 ge = new discreteEdge(gm, indexOfCuts[i]);
1938 gm->add(ge);
1939 }
1940 MEdgeVertex *v =
1941 new MEdgeVertex(ps[i].x(), ps[i].y(), ps[i].z(), ge, 0.0);
1942 if(f)
1943 fprintf(f, "SP(%g,%g,%g){%d};\n", ps[i].x(), ps[i].y(), ps[i].z(),
1944 ge->tag());
1945 vs.push_back(v);
1946 ge->mesh_vertices.push_back(v);
1947 }
1948 }
edgeCutsedgeCuts1949 edgeCuts() {}
1950 };
1951
addCut(const SPoint3 & p,const MEdge & e,int COUNT,int ID,std::map<MEdge,edgeCuts,MEdgeLessThan> & cuts)1952 static bool addCut(const SPoint3 &p, const MEdge &e, int COUNT, int ID,
1953 std::map<MEdge, edgeCuts, MEdgeLessThan> &cuts)
1954 {
1955 auto itc = cuts.find(e);
1956 if(itc != cuts.end()) {
1957 if(!itc->second.add(p, COUNT, ID)) return false;
1958 return true;
1959 }
1960 else {
1961 edgeCuts cc;
1962 if(!cc.add(p, COUNT, ID)) return false;
1963 cuts[e] = cc;
1964 return true;
1965 }
1966 }
1967
computeOneIsoRecur(MVertex * vsing,v2t_cont & adj,double VAL,MVertex * v0,MVertex * v1,SPoint3 & p,std::map<MVertex *,double> & pot,std::map<MEdge,int,MEdgeLessThan> & visited,std::map<MEdge,std::pair<std::map<MVertex *,double> *,double>,MEdgeLessThan> & cutGraphEnds,std::map<MEdge,MEdge,MEdgeLessThan> & d1,std::vector<groupOfCross2d> & G,FILE * f,int COUNT,std::map<MEdge,edgeCuts,MEdgeLessThan> & cuts,int & NB,int & circular)1968 static void computeOneIsoRecur(
1969 MVertex *vsing, v2t_cont &adj, double VAL, MVertex *v0, MVertex *v1,
1970 SPoint3 &p, std::map<MVertex *, double> &pot,
1971 std::map<MEdge, int, MEdgeLessThan> &visited,
1972 std::map<MEdge, std::pair<std::map<MVertex *, double> *, double>,
1973 MEdgeLessThan> &cutGraphEnds,
1974 std::map<MEdge, MEdge, MEdgeLessThan> &d1, std::vector<groupOfCross2d> &G,
1975 FILE *f, int COUNT, std::map<MEdge, edgeCuts, MEdgeLessThan> &cuts, int &NB,
1976 int &circular)
1977 {
1978 MEdge e(v0, v1);
1979
1980 MVertex vvv(p.x(), p.y(), p.z());
1981
1982 if(distance(&vvv, vsing) < 1.e-10) {
1983 circular++;
1984 return;
1985 }
1986
1987 bool added = addCut(p, e, COUNT, NB, cuts);
1988 if(!added) return;
1989
1990 NB++;
1991
1992 // visited[e] = NB;
1993
1994 if(d1.find(e) != d1.end()) {
1995 std::pair<std::map<MVertex *, double> *, double> aa =
1996 std::make_pair(&pot, VAL);
1997 std::pair<MEdge, std::pair<std::map<MVertex *, double> *, double> > p =
1998 std::make_pair(e, aa);
1999 cutGraphEnds.insert(p);
2000 }
2001 std::vector<MElement *> lst = adj[v0];
2002
2003 MVertex *vs[2] = {nullptr, nullptr};
2004 int count = 0;
2005 for(size_t i = 0; i < lst.size(); i++) {
2006 if((lst[i]->getVertex(0) == v0 && lst[i]->getVertex(1) == v1) ||
2007 (lst[i]->getVertex(0) == v1 && lst[i]->getVertex(1) == v0)) {
2008 vs[count++] = lst[i]->getVertex(2);
2009 }
2010 if((lst[i]->getVertex(0) == v0 && lst[i]->getVertex(2) == v1) ||
2011 (lst[i]->getVertex(0) == v1 && lst[i]->getVertex(2) == v0)) {
2012 vs[count++] = lst[i]->getVertex(1);
2013 }
2014 if((lst[i]->getVertex(1) == v0 && lst[i]->getVertex(2) == v1) ||
2015 (lst[i]->getVertex(1) == v1 && lst[i]->getVertex(2) == v0)) {
2016 vs[count++] = lst[i]->getVertex(0);
2017 }
2018 }
2019
2020 double U[2] = {pot[v0], pot[v1]};
2021 SPoint3 p0(v0->x(), v0->y(), v0->z());
2022 SPoint3 p1(v1->x(), v1->y(), v1->z());
2023 for(int i = 0; i < 2; i++) {
2024 if(vs[i]) {
2025 double U2 = pot[vs[i]];
2026 SPoint3 ppp(vs[i]->x(), vs[i]->y(), vs[i]->z());
2027 if((U[0] - VAL) * (U2 - VAL) <= 0) {
2028 double xi = coord1d(U[0], U2, VAL);
2029 SPoint3 pp = p0 * (1. - xi) + ppp * xi;
2030 fprintf(f, "SL(%g,%g,%g,%g,%g,%g){%d,%d};\n", p.x(), p.y(), p.z(),
2031 pp.x(), pp.y(), pp.z(), COUNT, COUNT);
2032 computeOneIsoRecur(vsing, adj, VAL, v0, vs[i], pp, pot, visited,
2033 cutGraphEnds, d1, G, f, COUNT, cuts, NB, circular);
2034 }
2035 else if((U[1] - VAL) * (U2 - VAL) <= 0) {
2036 double xi = coord1d(U[1], U2, VAL);
2037 SPoint3 pp = p1 * (1. - xi) + ppp * xi;
2038 fprintf(f, "SL(%g,%g,%g,%g,%g,%g){%d,%d};\n", p.x(), p.y(), p.z(),
2039 pp.x(), pp.y(), pp.z(), COUNT, COUNT);
2040 computeOneIsoRecur(vsing, adj, VAL, v1, vs[i], pp, pot, visited,
2041 cutGraphEnds, d1, G, f, COUNT, cuts, NB, circular);
2042 }
2043 else {
2044 printf("strange\n");
2045 }
2046 }
2047 }
2048 return;
2049 }
2050
computeOneIso(MVertex * vsing,v2t_cont & adj,double VAL,MVertex * v0,MVertex * v1,SPoint3 & p,std::map<MVertex *,double> * potU,std::map<MVertex *,double> * potV,std::map<MEdge,MEdge,MEdgeLessThan> & d1,std::vector<groupOfCross2d> & G,FILE * f,int COUNT,std::map<MEdge,edgeCuts,MEdgeLessThan> & cuts,std::vector<cutGraphPassage> & passages)2051 static int computeOneIso(MVertex *vsing, v2t_cont &adj, double VAL, MVertex *v0,
2052 MVertex *v1, SPoint3 &p,
2053 std::map<MVertex *, double> *potU,
2054 std::map<MVertex *, double> *potV,
2055 std::map<MEdge, MEdge, MEdgeLessThan> &d1,
2056 std::vector<groupOfCross2d> &G, FILE *f, int COUNT,
2057 std::map<MEdge, edgeCuts, MEdgeLessThan> &cuts,
2058 std::vector<cutGraphPassage> &passages)
2059 {
2060 std::map<MEdge, int, MEdgeLessThan> visited;
2061 std::map<MEdge, std::pair<std::map<MVertex *, double> *, double>,
2062 MEdgeLessThan>
2063 cutGraphEnds;
2064 int NB = 0;
2065 int circular = 0;
2066
2067 computeOneIsoRecur(vsing, adj, VAL, v0, v1, p, *potU, visited, cutGraphEnds,
2068 d1, G, f, COUNT, cuts, NB, circular);
2069
2070 int XX = 1;
2071 passages.clear();
2072 while(!cutGraphEnds.empty()) {
2073 MEdge e = (*cutGraphEnds.begin()).first;
2074 std::map<MVertex *, double> *POT = (*cutGraphEnds.begin()).second.first;
2075 VAL = (*cutGraphEnds.begin()).second.second;
2076 double xi = coord1d((*POT)[e.getVertex(0)], (*POT)[e.getVertex(1)], VAL);
2077 MEdge o = d1[e];
2078 p[0] = (1. - xi) * e.getVertex(0)->x() + xi * e.getVertex(1)->x();
2079 p[1] = (1. - xi) * e.getVertex(0)->y() + xi * e.getVertex(1)->y();
2080 p[2] = (1. - xi) * e.getVertex(0)->z() + xi * e.getVertex(1)->z();
2081 // printf("cutgaphends %lu
2082 // %lu\n",o.getVertex(0)->getNum(),o.getVertex(1)->getNum()); printf("%lu
2083 // ends to the cutgraph\n",cutGraphEnds.size());
2084 cutGraphEnds.erase(cutGraphEnds.begin());
2085 // visited.clear();
2086
2087 int ROT = 0;
2088 int maxCount = 0;
2089 int cutGraphId = -1;
2090 for(size_t i = 0; i < G.size(); i++) {
2091 int count = 0;
2092 count += (std::find(G[i].left.begin(), G[i].left.end(), o.getVertex(0)) !=
2093 G[i].left.end() ?
2094 1 :
2095 0);
2096 count += (std::find(G[i].left.begin(), G[i].left.end(), o.getVertex(1)) !=
2097 G[i].left.end() ?
2098 1 :
2099 0);
2100 count += (std::find(G[i].right.begin(), G[i].right.end(),
2101 e.getVertex(0)) != G[i].right.end() ?
2102 1 :
2103 0);
2104 count += (std::find(G[i].right.begin(), G[i].right.end(),
2105 e.getVertex(1)) != G[i].right.end() ?
2106 1 :
2107 0);
2108 count += (std::find(G[i].left.begin(), G[i].left.end(), e.getVertex(0)) !=
2109 G[i].left.end() ?
2110 1 :
2111 0);
2112 count += (std::find(G[i].left.begin(), G[i].left.end(), e.getVertex(1)) !=
2113 G[i].left.end() ?
2114 1 :
2115 0);
2116 count += (std::find(G[i].right.begin(), G[i].right.end(),
2117 o.getVertex(0)) != G[i].right.end() ?
2118 1 :
2119 0);
2120 count += (std::find(G[i].right.begin(), G[i].right.end(),
2121 o.getVertex(1)) != G[i].right.end() ?
2122 1 :
2123 0);
2124 if(count > maxCount) {
2125 maxCount = count;
2126 ROT = fabs(G[i].mat[0][0]) > .6 ? 0 : 1;
2127 cutGraphId = i;
2128 }
2129 }
2130 if(maxCount == 0) printf("IMPOSSIBLE\n");
2131
2132 int pot = POT == potU ? 0 : 1;
2133 // printf(" --> cutting cut graph %5d %5d\n",cutGraphId,
2134 // pot,passages.size());
2135 int count = 0;
2136 for(std::size_t k = 0; k < passages.size(); k++) {
2137 if(pot == passages[k]._uv && cutGraphId == passages[k]._id) count++;
2138 }
2139
2140 if(count > 20) {
2141 printf("CYCLE DETECTED for SING %lu : ", vsing->getNum());
2142 for(size_t k = 0; k < passages.size(); k++)
2143 printf("(%d,%d) ", passages[k]._id, passages[k]._uv);
2144 printf("\n");
2145 return -1;
2146 }
2147
2148 if(passages.empty() || passages[passages.size() - 1]._uv != pot ||
2149 passages[passages.size() - 1]._id != cutGraphId) {
2150 passages.push_back(cutGraphPassage(cutGraphId, pot, p));
2151 }
2152
2153 if(ROT) { POT = (POT == potU ? potV : potU); }
2154 else {
2155 }
2156 XX += ROT;
2157 // printf("XX = %d ROT = %d %lu\n",XX,ROT,cutGraphEnds.size());
2158 if(distance(e.getVertex(0), o.getVertex(0)) < 1.e-12)
2159 VAL = (1. - xi) * (*POT)[o.getVertex(0)] + xi * (*POT)[o.getVertex(1)];
2160 else
2161 VAL = (1. - xi) * (*POT)[o.getVertex(1)] + xi * (*POT)[o.getVertex(0)];
2162 computeOneIsoRecur(vsing, adj, VAL, o.getVertex(0), o.getVertex(1), p, *POT,
2163 visited, cutGraphEnds, d1, G, f, COUNT, cuts, NB,
2164 circular);
2165 if(XX > 1200) break;
2166 }
2167 return circular;
2168 }
2169
computeIso(MVertex * vsing,v2t_cont & adj,double u,std::map<MVertex *,double> & potU,std::map<MVertex *,double> & potV,FILE * f,std::map<MEdge,MEdge,MEdgeLessThan> & d1,std::vector<groupOfCross2d> & G,int DIR,std::map<MEdge,edgeCuts,MEdgeLessThan> & cuts,std::vector<cutGraphPassage> & passages)2170 static bool computeIso(MVertex *vsing, v2t_cont &adj, double u,
2171 std::map<MVertex *, double> &potU,
2172 std::map<MVertex *, double> &potV, FILE *f,
2173 std::map<MEdge, MEdge, MEdgeLessThan> &d1,
2174 std::vector<groupOfCross2d> &G, int DIR,
2175 std::map<MEdge, edgeCuts, MEdgeLessThan> &cuts,
2176 std::vector<cutGraphPassage> &passages)
2177 {
2178 int COUNT = 100 * vsing->getNum() + 10 * DIR;
2179 std::vector<MElement *> faces = adj[vsing];
2180 for(size_t i = 0; i < faces.size(); i++) {
2181 MVertex *v0 = faces[i]->getVertex(0);
2182 MVertex *v1 = faces[i]->getVertex(1);
2183 MVertex *v2 = faces[i]->getVertex(2);
2184 double U0 = potU[v0];
2185 double U1 = potU[v1];
2186 double U2 = potU[v2];
2187 SPoint3 p0(v0->x(), v0->y(), v0->z());
2188 SPoint3 p1(v1->x(), v1->y(), v1->z());
2189 SPoint3 p2(v2->x(), v2->y(), v2->z());
2190 int circular = 0;
2191 if(v2 == vsing && (U0 - u) * (U1 - u) <= 0) {
2192 double xi = coord1d(U0, U1, u);
2193 SPoint3 pp = p0 * (1 - xi) + p1 * xi;
2194 circular = computeOneIso(vsing, adj, u, v0, v1, pp, &potU, &potV, d1, G,
2195 f, COUNT++, cuts, passages);
2196 }
2197 else if(v1 == vsing && (U0 - u) * (U2 - u) <= 0) {
2198 double xi = coord1d(U0, U2, u);
2199 SPoint3 pp = p0 * (1 - xi) + p2 * xi;
2200 circular = computeOneIso(vsing, adj, u, v0, v2, pp, &potU, &potV, d1, G,
2201 f, COUNT++, cuts, passages);
2202 }
2203 else if(v0 == vsing && (U1 - u) * (U2 - u) <= 0) {
2204 double xi = coord1d(U1, U2, u);
2205 SPoint3 pp = p1 * (1 - xi) + p2 * xi;
2206 circular = computeOneIso(vsing, adj, u, v1, v2, pp, &potU, &potV, d1, G,
2207 f, COUNT++, cuts, passages);
2208 }
2209 if(circular == 2) printf("ISO %d is circular %d\n", COUNT - 1, circular);
2210 if(circular == -1) {
2211 printf("ISO %d is CYCLIC\n", COUNT - 1);
2212 return false;
2213 }
2214 }
2215 return true;
2216 }
2217
computeIsos(GModel * gm,std::vector<GFace * > & faces,std::set<MVertex *,MVertexPtrLessThan> singularities,std::map<MEdge,cross2d,MEdgeLessThan> & C,std::map<MVertex *,MVertex *,MVertexPtrLessThan> & new2old,std::map<MEdge,MEdge,MEdgeLessThan> & duplicateEdges,std::vector<std::vector<cross2d * >> & groups,std::vector<std::vector<cross2d * >> & groups_cg,std::map<MVertex *,double> & potU,std::map<MVertex *,double> & potV,std::set<MEdge,MEdgeLessThan> & cutG,std::vector<groupOfCross2d> & G,std::map<MEdge,edgeCuts,MEdgeLessThan> & cuts,std::vector<cutGraphPassage> & passages)2218 static bool computeIsos(
2219 GModel *gm, std::vector<GFace *> &faces,
2220 std::set<MVertex *, MVertexPtrLessThan> singularities,
2221 std::map<MEdge, cross2d, MEdgeLessThan> &C,
2222 std::map<MVertex *, MVertex *, MVertexPtrLessThan> &new2old,
2223 std::map<MEdge, MEdge, MEdgeLessThan> &duplicateEdges,
2224 std::vector<std::vector<cross2d *> > &groups,
2225 std::vector<std::vector<cross2d *> > &groups_cg,
2226 std::map<MVertex *, double> &potU, std::map<MVertex *, double> &potV,
2227 std::set<MEdge, MEdgeLessThan> &cutG, std::vector<groupOfCross2d> &G,
2228 std::map<MEdge, edgeCuts, MEdgeLessThan> &cuts,
2229 std::vector<cutGraphPassage> &passages)
2230 {
2231 v2t_cont adj;
2232 for(size_t i = 0; i < faces.size(); i++) {
2233 buildVertexToElement(faces[i]->triangles, adj);
2234 }
2235
2236 {
2237 auto it = new2old.begin();
2238 for(; it != new2old.end(); ++it) {
2239 if(singularities.find(it->second) != singularities.end()) {
2240 singularities.insert(it->first);
2241 }
2242 }
2243 }
2244
2245 std::map<MVertex *, MVertex *, MVertexPtrLessThan> duplicates;
2246 {
2247 auto it = new2old.begin();
2248 for(; it != new2old.end(); ++it) {
2249 duplicates[it->first] = it->second;
2250 duplicates[it->second] = it->first;
2251 }
2252 }
2253
2254 std::map<MEdge, MEdge, MEdgeLessThan> d1;
2255 {
2256 for(size_t i = 0; i < G.size(); i++) {
2257 for(size_t j = 0; j < G[i].side.size(); j++) {
2258 for(size_t k = 0; k < 3; k++) {
2259 MVertex *v0 = G[i].side[j]->getVertex(k);
2260 MVertex *v1 = G[i].side[j]->getVertex((k + 1) % 3);
2261 int J = -1, I = -1;
2262 for(size_t l = 0; l < G[i].left.size(); l++) {
2263 if(G[i].left[l] == v0) { I = l; }
2264 if(G[i].left[l] == v1) { J = l; }
2265 }
2266 if(I >= 0 && J >= 0) {
2267 MEdge l(G[i].left[I], G[i].left[J]);
2268 MEdge r(G[i].right[I], G[i].right[J]);
2269 d1[l] = r;
2270 d1[r] = l;
2271 }
2272 }
2273 }
2274 if(G[i].singularities.size() == 1) {
2275 MEdge l(G[i].singularities[0], G[i].left[G[i].left.size() - 1]);
2276 MEdge r(G[i].singularities[0], G[i].right[G[i].right.size() - 1]);
2277 d1[l] = r;
2278 d1[r] = l;
2279 }
2280 }
2281 }
2282 std::string fn = gm->getName() + "_QLayoutResults.pos";
2283 FILE *f = fopen(fn.c_str(), "w");
2284 fprintf(f, "View\"Big Cut\"{\n");
2285
2286 auto it = singularities.begin();
2287 for(; it != singularities.end(); ++it) {
2288 GEntity *ge = (*it)->onWhat();
2289 if(ge->dim() == 2 || ge->edges().size() == 0) {
2290 printf("%lu %d %d %lu %22.15E %22.15E\n", ge->edges().size(), ge->tag(),
2291 ge->dim(), singularities.size(), potU[*it], potV[*it]);
2292 bool success = computeIso(*it, adj, potU[*it], potU, potV, f, d1, G, 0,
2293 cuts, passages);
2294 if(!success) {
2295 printf("CYCLIC STUFF\n");
2296 // return false;
2297 }
2298 success = computeIso(*it, adj, potV[*it], potV, potU, f, d1, G, 1, cuts,
2299 passages);
2300 if(!success) {
2301 printf("CYCLIC STUFF\n");
2302 // return false;
2303 }
2304 }
2305 }
2306 fprintf(f, "};\n");
2307 fclose(f);
2308 return true;
2309 }
2310
getAllConnectedTriangles(cross2d * start,std::vector<cross2d * > & group,std::set<MVertex *,MVertexPtrLessThan> & isolated_singularities,std::set<MVertex *,MVertexPtrLessThan> & all,std::set<MTriangle * > & t,std::set<MTriangle * > & allTrianglesConsidered)2311 void getAllConnectedTriangles(
2312 cross2d *start, std::vector<cross2d *> &group,
2313 std::set<MVertex *, MVertexPtrLessThan> &isolated_singularities,
2314 std::set<MVertex *, MVertexPtrLessThan> &all, std::set<MTriangle *> &t,
2315 std::set<MTriangle *> &allTrianglesConsidered)
2316 {
2317 std::set<cross2d *> touched;
2318
2319 // printf("group %lu isolated singularities\n",
2320 // isolated_singularities.size());
2321
2322 for(size_t i = 0; i < group.size(); i++) {
2323 if(isolated_singularities.find(group[i]->_e.getVertex(0)) ==
2324 isolated_singularities.end())
2325 all.insert(group[i]->_e.getVertex(0));
2326 if(isolated_singularities.find(group[i]->_e.getVertex(1)) ==
2327 isolated_singularities.end())
2328 all.insert(group[i]->_e.getVertex(1));
2329 }
2330
2331 if(allTrianglesConsidered.find(start->_t[0]) !=
2332 allTrianglesConsidered.end()) {
2333 if(!start->_cneighbors[0]->inCutGraph)
2334 start = start->_cneighbors[0];
2335 else if(!start->_cneighbors[1]->inCutGraph)
2336 start = start->_cneighbors[1];
2337 else
2338 printf("error\n");
2339 }
2340 else if(start->_cneighbors.size() == 4 &&
2341 allTrianglesConsidered.find(start->_t[1]) !=
2342 allTrianglesConsidered.end()) {
2343 if(start->_cneighbors.size() == 4 && !start->_cneighbors[2]->inCutGraph)
2344 start = start->_cneighbors[2];
2345 else if(start->_cneighbors.size() == 4 &&
2346 !start->_cneighbors[3]->inCutGraph)
2347 start = start->_cneighbors[3];
2348 else
2349 printf("error\n");
2350 }
2351 else {
2352 if(!start->_cneighbors[0]->inCutGraph)
2353 start = start->_cneighbors[0];
2354 else if(!start->_cneighbors[1]->inCutGraph)
2355 start = start->_cneighbors[1];
2356 else if(start->_cneighbors.size() == 4 &&
2357 !start->_cneighbors[2]->inCutGraph)
2358 start = start->_cneighbors[2];
2359 else if(start->_cneighbors.size() == 4 &&
2360 !start->_cneighbors[3]->inCutGraph)
2361 start = start->_cneighbors[3];
2362 else
2363 printf("error\n");
2364 }
2365
2366 std::stack<cross2d *> _s;
2367 _s.push(start);
2368
2369 while(!_s.empty()) {
2370 start = _s.top();
2371 touched.insert(start);
2372 _s.pop();
2373 for(size_t i = 0; i < start->_t.size(); i++) {
2374 t.insert(start->_t[i]);
2375 allTrianglesConsidered.insert(start->_t[i]);
2376 }
2377
2378 for(size_t i = 0; i < start->_cneighbors.size(); i++) {
2379 cross2d *c = start->_cneighbors[i];
2380 if(!c->inCutGraph && touched.find(c) == touched.end()) {
2381 if(all.find(c->_e.getVertex(0)) != all.end() ||
2382 all.find(c->_e.getVertex(1)) != all.end()) {
2383 _s.push(c);
2384 }
2385 }
2386 }
2387 }
2388 }
2389
computeLeftRight(groupOfCross2d & g,MVertex ** left,MVertex ** right)2390 static bool computeLeftRight(groupOfCross2d &g, MVertex **left, MVertex **right)
2391 {
2392 for(size_t i = 0; i < g.side.size(); i++) {
2393 if(g.side[i]->getVertex(0) == *right || g.side[i]->getVertex(1) == *right ||
2394 g.side[i]->getVertex(2) == *right) {
2395 MVertex *temp = *left;
2396 *left = *right;
2397 *right = temp;
2398 return true;
2399 }
2400 if(g.side[i]->getVertex(0) == *left || g.side[i]->getVertex(1) == *left ||
2401 g.side[i]->getVertex(2) == *left) {
2402 return true;
2403 }
2404 }
2405 return false;
2406 }
2407
createJumpyPairs(groupOfCross2d & g,std::set<MVertex *,MVertexPtrLessThan> & singularities,std::set<MVertex *,MVertexPtrLessThan> & boundaries,std::multimap<MVertex *,MVertex *,MVertexPtrLessThan> & old2new)2408 static void createJumpyPairs(
2409 groupOfCross2d &g, std::set<MVertex *, MVertexPtrLessThan> &singularities,
2410 std::set<MVertex *, MVertexPtrLessThan> &boundaries,
2411 std::multimap<MVertex *, MVertex *, MVertexPtrLessThan> &old2new)
2412 {
2413 std::set<MVertex *, MVertexPtrLessThan> touched;
2414
2415 // printf("GROUP %d \n",g.groupId);
2416 for(size_t i = 0; i < g.crosses.size(); ++i) {
2417 cross2d *c = g.crosses[i];
2418 for(size_t j = 0; j < 2; j++) {
2419 MVertex *vv = c->_e.getVertex(j);
2420 if(touched.find(vv) == touched.end()) {
2421 touched.insert(vv);
2422 MTriangle *t1 = c->_t[0];
2423 MTriangle *t2 = c->_t[1];
2424 MVertex *v0 = nullptr;
2425 MVertex *v1 = nullptr;
2426 if(t1->getVertex(0) == vv || t1->getVertex(1) == vv ||
2427 t1->getVertex(2) == vv) {
2428 if(v0 == nullptr)
2429 v0 = vv;
2430 else if(v1 == nullptr)
2431 v1 = vv;
2432 else
2433 Msg::Error("error in JumpyPairs 1");
2434 }
2435 if(t2->getVertex(0) == vv || t2->getVertex(1) == vv ||
2436 t2->getVertex(2) == vv) {
2437 if(v0 == nullptr)
2438 v0 = vv;
2439 else if(v1 == nullptr)
2440 v1 = vv;
2441 else
2442 Msg::Error("error in JumpyPairs 1");
2443 }
2444 for(auto it = old2new.lower_bound(vv); it != old2new.upper_bound(vv);
2445 ++it) {
2446 MVertex *vvv = it->second;
2447 if(t1->getVertex(0) == vvv || t1->getVertex(1) == vvv ||
2448 t1->getVertex(2) == vvv) {
2449 if(v0 == nullptr)
2450 v0 = vvv;
2451 else if(v1 == nullptr)
2452 v1 = vvv;
2453 else
2454 Msg::Error("error in JumpyPairs 1");
2455 }
2456 if(t2->getVertex(0) == vvv || t2->getVertex(1) == vvv ||
2457 t2->getVertex(2) == vvv) {
2458 if(v0 == nullptr)
2459 v0 = vvv;
2460 else if(v1 == nullptr)
2461 v1 = vvv;
2462 else
2463 Msg::Error("error in JumpyPairs 2");
2464 }
2465 }
2466 if(!v1 || !v0) Msg::Error("error in JumpyPairs 3");
2467 if(computeLeftRight(g, &v0, &v1)) {
2468 if(boundaries.find(vv) != boundaries.end()) {
2469 g.left.insert(g.left.begin(), v0);
2470 g.right.insert(g.right.begin(), v1);
2471 }
2472 else {
2473 g.left.push_back(v0);
2474 g.right.push_back(v1);
2475 }
2476 }
2477 else
2478 Msg::Error("error in jumpy pairs %lu \n", vv->getNum());
2479 }
2480 else if(singularities.find(vv) != singularities.end()) {
2481 g.singularities.push_back(vv);
2482 }
2483 // else if (singularities.find(vv) == singularities.end()){
2484 // printf("ERROR --> no counterpart vertex in the cut graph\n");
2485 // }
2486 }
2487 }
2488 // printf("GRoup %d mat [%g,%g] [%g,%g] %lu nodes on each side \n"
2489 // ,g.groupId,g.mat[0][0],g.mat[0][1],g.mat[1][0],g.mat[1][1],
2490 // g.left.size());
2491 }
2492
2493 static void
analyzeGroup(std::vector<cross2d * > & group,groupOfCross2d & g,std::map<MTriangle *,SVector3> & d,std::map<MTriangle *,SVector3> & d2,v2t_cont & adj,std::set<MVertex *,MVertexPtrLessThan> & isolated_singularities,std::set<MVertex *,MVertexPtrLessThan> & boundaries,std::set<MTriangle * > & allTrianglesConsidered)2494 analyzeGroup(std::vector<cross2d *> &group, groupOfCross2d &g,
2495 std::map<MTriangle *, SVector3> &d,
2496 std::map<MTriangle *, SVector3> &d2, v2t_cont &adj,
2497 std::set<MVertex *, MVertexPtrLessThan> &isolated_singularities,
2498 std::set<MVertex *, MVertexPtrLessThan> &boundaries,
2499 std::set<MTriangle *> &allTrianglesConsidered)
2500 {
2501 g.crosses = group;
2502 double MAX = 0.0;
2503 for(size_t i = 0; i < g.crosses.size(); ++i) {
2504 cross2d *c = g.crosses[i];
2505 if(c->_t.size() == 2) {
2506 SVector3 t1 = d[c->_t[0]];
2507 SVector3 t2 = d[c->_t[1]];
2508 MAX = std::max(dot(t1, t2), MAX);
2509 }
2510 }
2511 if(MAX > .8)
2512 g.rot = false;
2513 else
2514 g.rot = true;
2515 for(size_t i = 0; i < g.crosses.size(); ++i) {
2516 cross2d *c = g.crosses[i];
2517 c->rotation = g.rot;
2518 }
2519
2520 std::set<MTriangle *> t;
2521 std::set<MVertex *, MVertexPtrLessThan> all;
2522 getAllConnectedTriangles(group[0], group, isolated_singularities, all, t,
2523 allTrianglesConsidered);
2524 g.side.insert(g.side.begin(), t.begin(), t.end());
2525 g.vertices.insert(g.vertices.begin(), all.begin(), all.end());
2526
2527 // compute which rotation ...
2528 for(size_t i = 0; i < g.crosses.size(); ++i) {
2529 cross2d *c = g.crosses[i];
2530 if(c->_t.size() == 2) {
2531 if(t.find(c->_t[0]) != t.end()) {
2532 g.mat[0][0] = dot(d[c->_t[0]], d[c->_t[1]]);
2533 g.mat[0][1] = dot(d[c->_t[0]], d2[c->_t[1]]);
2534 g.mat[1][0] = dot(d2[c->_t[0]], d[c->_t[1]]);
2535 g.mat[1][1] = dot(d2[c->_t[0]], d2[c->_t[1]]);
2536 }
2537 else if(t.find(c->_t[1]) != t.end()) {
2538 g.mat[0][0] = dot(d[c->_t[0]], d[c->_t[1]]);
2539 g.mat[1][0] = dot(d[c->_t[0]], d2[c->_t[1]]);
2540 g.mat[0][1] = dot(d2[c->_t[0]], d[c->_t[1]]);
2541 g.mat[1][1] = dot(d2[c->_t[0]], d2[c->_t[1]]);
2542 }
2543 }
2544 }
2545 for(int j = 0; j < 2; j++) {
2546 for(int k = 0; k < 2; k++) {
2547 if(g.mat[j][k] > .7)
2548 g.mat[j][k] = 1;
2549 else if(g.mat[j][k] < -.7)
2550 g.mat[j][k] = -1;
2551 else
2552 g.mat[j][k] = 0;
2553 }
2554 }
2555 }
2556
2557 ///--- class containing the data
2558 class quadLayoutData {
2559 public:
2560 GModel *gm;
2561 std::vector<GFace *> f;
2562 std::map<MEdge, cross2d, MEdgeLessThan> C;
2563 dofManager<double> *myAssembler;
2564 std::set<MVertex *, MVertexPtrLessThan> vs;
2565 std::set<MEdge, MEdgeLessThan> cutG;
2566 std::set<MVertex *, MVertexPtrLessThan> singularities;
2567 std::map<MVertex *, int> indices;
2568 std::map<MVertex *, double> gaussianCurvatures;
2569 std::set<MVertex *, MVertexPtrLessThan> boundaries;
2570 std::set<MVertex *, MVertexPtrLessThan> corners;
2571 std::vector<std::vector<cross2d *> > groups;
2572 std::vector<std::vector<cross2d *> > groups_cg;
2573 std::map<MVertex *, MVertex *, MVertexPtrLessThan> new2old;
2574 std::string modelName;
2575 std::map<MTriangle *, SVector3> d0, d1;
2576 std::vector<groupOfCross2d> G;
2577
printTheta(const char * name)2578 void printTheta(const char *name)
2579 {
2580 std::string fn = modelName + "_" + name + ".pos";
2581 FILE *of = fopen(fn.c_str(), "w");
2582 fprintf(of, "View \"Theta\"{\n");
2583 for(size_t i = 0; i < f.size(); i++) {
2584 for(size_t j = 0; j < f[i]->triangles.size(); j++) {
2585 MTriangle *t = f[i]->triangles[j];
2586 auto it0 = C.find(t->getEdge(0));
2587 auto it1 = C.find(t->getEdge(1));
2588 auto it2 = C.find(t->getEdge(2));
2589
2590 SVector3 d0 = it0->second.o_i;
2591 SVector3 d1 = it1->second.o_i;
2592 SVector3 d2 = it2->second.o_i;
2593 double a = atan2(d0.y(), d0.x());
2594 double b = atan2(d1.y(), d1.x());
2595 double c = atan2(d2.y(), d2.x());
2596 it0->second.normalize(a);
2597 it0->second.normalize(b);
2598 it0->second.normalize(c);
2599 double A = c + a - b;
2600 double B = a + b - c;
2601 double C = b + c - a;
2602 it0->second.normalize(A);
2603 it0->second.normalize(B);
2604 it0->second.normalize(C);
2605 fprintf(of, "ST(%g,%g,%g,%g,%g,%g,%g,%g,%g){%g,%g,%g};\n",
2606 t->getVertex(0)->x(), t->getVertex(0)->y(),
2607 t->getVertex(0)->z(), t->getVertex(1)->x(),
2608 t->getVertex(1)->y(), t->getVertex(1)->z(),
2609 t->getVertex(2)->x(), t->getVertex(2)->y(),
2610 t->getVertex(2)->z(), A, B, C);
2611 }
2612 }
2613 fprintf(of, "};\n");
2614 fclose(of);
2615 }
2616
printCross(const char * name)2617 void printCross(const char *name)
2618 {
2619 std::string fn = modelName + "_" + name + ".pos";
2620 FILE *of = fopen(fn.c_str(), "w");
2621 fprintf(of, "View \"Direction fields\"{\n");
2622 auto it = C.begin();
2623 for(it = C.begin(); it != C.end(); ++it) {
2624 double a0 = it->second._a;
2625 MEdge e0 = it->second._e;
2626 SVector3 d1 = (it->second._tgt * cos(a0) + it->second._tgt2 * sin(a0));
2627 SVector3 d2 = (it->second._tgt * (-sin(a0)) + it->second._tgt2 * cos(a0));
2628 SVector3 d3 = (it->second._tgt * (-cos(a0)) - it->second._tgt2 * sin(a0));
2629 SVector3 d4 = (it->second._tgt * sin(a0) - it->second._tgt2 * cos(a0));
2630
2631 for(size_t I = 0; I < it->second._t.size(); I++) {
2632 fprintf(of, "VP(%g,%g,%g){%g,%g,%g};\n",
2633 0.5 * (e0.getVertex(0)->x() + e0.getVertex(1)->x()),
2634 0.5 * (e0.getVertex(0)->y() + e0.getVertex(1)->y()),
2635 0.5 * (e0.getVertex(0)->z() + e0.getVertex(1)->z()), d1.x(),
2636 d1.y(), d1.z());
2637 fprintf(of, "VP(%g,%g,%g){%g,%g,%g};\n",
2638 0.5 * (e0.getVertex(0)->x() + e0.getVertex(1)->x()),
2639 0.5 * (e0.getVertex(0)->y() + e0.getVertex(1)->y()),
2640 0.5 * (e0.getVertex(0)->z() + e0.getVertex(1)->z()), d2.x(),
2641 d2.y(), d2.z());
2642 fprintf(of, "VP(%g,%g,%g){%g,%g,%g};\n",
2643 0.5 * (e0.getVertex(0)->x() + e0.getVertex(1)->x()),
2644 0.5 * (e0.getVertex(0)->y() + e0.getVertex(1)->y()),
2645 0.5 * (e0.getVertex(0)->z() + e0.getVertex(1)->z()), d3.x(),
2646 d3.y(), d3.z());
2647 fprintf(of, "VP(%g,%g,%g){%g,%g,%g};\n",
2648 0.5 * (e0.getVertex(0)->x() + e0.getVertex(1)->x()),
2649 0.5 * (e0.getVertex(0)->y() + e0.getVertex(1)->y()),
2650 0.5 * (e0.getVertex(0)->z() + e0.getVertex(1)->z()), d4.x(),
2651 d4.y(), d4.z());
2652 }
2653 }
2654 fprintf(of, "};\n");
2655 fclose(of);
2656 }
2657
computeCrossFieldExtrinsic(double tol,size_t nIterLaplace=2000)2658 int computeCrossFieldExtrinsic(double tol, size_t nIterLaplace = 2000)
2659 {
2660 std::vector<cross2d *> pc;
2661 for(auto it = C.begin(); it != C.end(); ++it) pc.push_back(&(it->second));
2662
2663 size_t ITER = 0;
2664 while(ITER++ < nIterLaplace) {
2665 if(ITER % 200 == 0) std::random_shuffle(pc.begin(), pc.end());
2666 for(size_t i = 0; i < pc.size(); i++) pc[i]->average_init();
2667 if(ITER % 1000 == 0) Msg::Info("Linear smooth : iter %6lu", ITER);
2668 }
2669
2670 for(size_t i = 0; i < pc.size(); i++) pc[i]->computeVector();
2671
2672 fastImplementationExtrinsic(C, tol);
2673
2674 for(size_t i = 0; i < pc.size(); i++) pc[i]->computeAngle();
2675
2676 return 0;
2677 }
2678
printScalar(dofManager<double> * dof,char c)2679 void printScalar(dofManager<double> *dof, char c)
2680 {
2681 std::string fn = modelName + "_" + c + ".pos";
2682
2683 FILE *_f = fopen(fn.c_str(), "w");
2684 fprintf(_f, "View \"H\"{\n");
2685
2686 for(size_t i = 0; i < f.size(); i++) {
2687 for(size_t j = 0; j < f[i]->triangles.size(); j++) {
2688 MTriangle *t = f[i]->triangles[j];
2689 double a, b, c;
2690 dof->getDofValue(t->getVertex(0), 0, 1, a);
2691 dof->getDofValue(t->getVertex(1), 0, 1, b);
2692 dof->getDofValue(t->getVertex(2), 0, 1, c);
2693 fprintf(_f, "ST(%g,%g,%g,%g,%g,%g,%g,%g,%g){%g,%g,%g};\n",
2694 t->getVertex(0)->x(), t->getVertex(0)->y(),
2695 t->getVertex(0)->z(), t->getVertex(1)->x(),
2696 t->getVertex(1)->y(), t->getVertex(1)->z(),
2697 t->getVertex(2)->x(), t->getVertex(2)->y(),
2698 t->getVertex(2)->z(), a, b, c);
2699 }
2700 }
2701 fprintf(_f, "};\n");
2702 fclose(_f);
2703 }
2704
computeCrossFieldAndH()2705 int computeCrossFieldAndH()
2706 {
2707 computeCrossFieldExtrinsic(1.e-9);
2708 myAssembler = computeH(gm, f, vs, C);
2709 // printScalar(myAssembler,'H');
2710 computeSingularities(C, singularities, indices, myAssembler);
2711 // print_H_and_Cross(gm, f, C, *myAssembler, singularities);
2712 return 1;
2713 }
2714
computeHFromSingularities(std::map<MVertex *,int> & sing,int nbTurns)2715 dofManager<double> *computeHFromSingularities(std::map<MVertex *, int> &sing,
2716 int nbTurns)
2717 {
2718 #if defined(HAVE_PETSC)
2719 linearSystemPETSc<double> *_lsys = new linearSystemPETSc<double>;
2720 #elif defined(HAVE_GMM)
2721 linearSystemGmm<double> *_lsys = new linearSystemGmm<double>;
2722 #else
2723 linearSystemFull<double> *_lsys = new linearSystemFull<double>;
2724 #endif
2725
2726 dofManager<double> *dof = new dofManager<double>(_lsys);
2727
2728 auto it = C.begin();
2729 std::vector<MEdge> edges;
2730 for(; it != C.end(); ++it) {
2731 if(it->second.inBoundary) { edges.push_back(it->first); }
2732 }
2733 std::vector<std::vector<MVertex *> > vsorted;
2734 SortEdgeConsecutive(edges, vsorted);
2735
2736 // AVERAGE
2737 dof->numberVertex(*vs.begin(), 1, 1);
2738
2739 for(auto it = vs.begin(); it != vs.end(); ++it) {
2740 dof->numberVertex(*it, 0, 1);
2741 }
2742
2743 simpleFunction<double> ONE(1.0);
2744 laplaceTerm l(nullptr, 1, &ONE);
2745
2746 std::set<GEntity *> firsts;
2747 for(size_t i = 0; i < f.size(); i++) {
2748 std::vector<GEdge *> e = f[i]->edges();
2749 if(e.size()) firsts.insert(e[0]);
2750 // printf("--> %lu\n",e[0]->tag());
2751 for(size_t j = 0; j < f[i]->triangles.size(); j++) {
2752 MTriangle *t = f[i]->triangles[j];
2753 SElement se(t);
2754 l.addToMatrix(*dof, &se);
2755 }
2756 }
2757
2758 for(size_t j = 0; j < vsorted.size(); ++j) {
2759 if(vsorted[j][0] == vsorted[j][vsorted[j].size() - 1]) {
2760 vsorted[j].erase(vsorted[j].begin());
2761 }
2762
2763 std::vector<double> CURVATURE;
2764 CURVATURE.resize(vsorted[j].size());
2765 for(size_t i = 0; i < vsorted[j].size(); ++i) { CURVATURE[i] = 0.0; }
2766 double SUM = 0.0;
2767 for(size_t i = 0; i < vsorted[j].size(); ++i) {
2768 MVertex *vi = vsorted[j][i];
2769 MVertex *vip = vsorted[j][(i + 1) % vsorted[j].size()];
2770 MVertex *vim =
2771 vsorted[j][(i + vsorted[j].size() - 1) % vsorted[j].size()];
2772 SVector3 vv(vip->x() - vi->x(), vip->y() - vi->y(), vip->z() - vi->z());
2773 SVector3 ww(vi->x() - vim->x(), vi->y() - vim->y(), vi->z() - vim->z());
2774 vv.normalize();
2775 ww.normalize();
2776 SVector3 xx = crossprod(vv, ww);
2777 double ccos = dot(vv, ww);
2778 double ANGLE = atan2(xx.norm(), ccos);
2779 xx.normalize();
2780
2781 MEdge edze(vi, vim);
2782 auto itip = C.find(edze);
2783 double sign = 1;
2784 if(itip != C.end()) {
2785 MTriangle *tt = itip->second._t[0];
2786 MVertex *vrv;
2787 if(tt->getVertex(0) != vi && tt->getVertex(0) != vim) {
2788 vrv = tt->getVertex(0);
2789 }
2790 else if(tt->getVertex(1) != vi && tt->getVertex(1) != vim) {
2791 vrv = tt->getVertex(1);
2792 }
2793 else
2794 vrv = tt->getVertex(2);
2795 SVector3 aa(vrv->x() - vim->x(), vrv->y() - vim->y(),
2796 vrv->z() - vim->z());
2797 SVector3 zz = crossprod(aa, ww);
2798 zz.normalize();
2799 sign = -dot(zz, xx); // > 0 ? -1 : 1;
2800 // sign = dot(zz, xx) > 0 ? -1 : 1;
2801 }
2802 else
2803 printf("ARGH\n");
2804 // if (vsorted.size() == 1)sign = -1;
2805 CURVATURE[i] += ANGLE * sign;
2806 // printf("%12.5E\n",sign);
2807 }
2808
2809 for(size_t i = 0; i < vsorted[j].size(); ++i) { SUM += CURVATURE[i]; }
2810
2811 // printf("%22.15E %22.15E\n",SUM, CORR );
2812 for(size_t i = 0; i < vsorted[j].size(); ++i) {
2813 Dof E(vsorted[j][i]->getNum(), Dof::createTypeWithTwoInts(0, 1));
2814 _lsys->addToRightHandSide(dof->getDofNumber(E), CURVATURE[i]);
2815 }
2816 }
2817
2818 for(auto it = gaussianCurvatures.begin(); it != gaussianCurvatures.end();
2819 ++it) {
2820 Dof E(it->first->getNum(), Dof::createTypeWithTwoInts(0, 1));
2821 //_lsys->addToRightHandSide(dof->getDofNumber(E),-it->second);
2822 }
2823
2824 double SSUM = 0;
2825 for(auto it = sing.begin(); it != sing.end(); ++it) {
2826 Dof E(it->first->getNum(), Dof::createTypeWithTwoInts(0, 1));
2827 _lsys->addToRightHandSide(dof->getDofNumber(E),
2828 2.0 * M_PI * (double)it->second / nbTurns);
2829 SSUM += 2.0 * M_PI * (double)it->second / nbTurns;
2830 }
2831
2832 // FIX DE LA MORT
2833 // AVERAGE
2834 Dof EAVG((*vs.begin())->getNum(), Dof::createTypeWithTwoInts(1, 1));
2835
2836 for(auto it = vs.begin(); it != vs.end(); ++it) {
2837 Dof E((*it)->getNum(), Dof::createTypeWithTwoInts(0, 1));
2838 dof->assemble(EAVG, E, 1);
2839 dof->assemble(E, EAVG, 1);
2840 }
2841
2842 _lsys->systemSolve();
2843 return dof;
2844 }
2845
computeHFromSingularities(std::map<MVertex *,int> & s)2846 int computeHFromSingularities(std::map<MVertex *, int> &s)
2847 {
2848 myAssembler = computeHFromSingularities(s, 4);
2849 for(auto it = s.begin(); it != s.end(); ++it) {
2850 singularities.insert(it->first);
2851 }
2852 // printScalar(myAssembler, 'H');
2853 return 1;
2854 }
2855
2856 //---------------------------------------------------------------------------
2857
computeThetaUsingHCrouzeixRaviart(std::map<int,std::vector<double>> & dataTHETA)2858 void computeThetaUsingHCrouzeixRaviart(
2859 std::map<int, std::vector<double> > &dataTHETA)
2860 {
2861 #if defined(HAVE_PETSC)
2862 linearSystemPETSc<double> *_lsys = new linearSystemPETSc<double>;
2863 #elif defined(HAVE_GMM)
2864 linearSystemGmm<double> *_lsys = new linearSystemGmm<double>;
2865 #else
2866 linearSystemFull<double> *_lsys = new linearSystemFull<double>;
2867 #endif
2868 dofManager<double> *theta = new dofManager<double>(_lsys);
2869
2870 std::map<MEdge, size_t, MEdgeLessThan> aaa;
2871 size_t count = 0;
2872 for(size_t i = 0; i < f.size(); i++) {
2873 for(size_t j = 0; j < f[i]->triangles.size(); j++) {
2874 for(size_t k = 0; k < 3; k++) {
2875 if(aaa.find(f[i]->triangles[j]->getEdge(k)) == aaa.end()) {
2876 Dof EdgeDof(count, Dof::createTypeWithTwoInts(0, 1));
2877 theta->numberDof(EdgeDof);
2878 aaa[f[i]->triangles[j]->getEdge(k)] = count++;
2879 }
2880 }
2881 }
2882 }
2883
2884 for(size_t i = 0; i < f.size(); i++) {
2885 for(size_t j = 0; j < f[i]->triangles.size(); j++) {
2886 MTriangle *t = f[i]->triangles[j];
2887 double V = t->getVolume();
2888 double g1[3], g2[3], g3[3];
2889 double a[3];
2890 a[0] = 1;
2891 a[1] = 0;
2892 a[2] = 0;
2893 t->interpolateGrad(a, 0, 0, 0, g1);
2894 a[0] = 0;
2895 a[1] = 1;
2896 a[2] = 0;
2897 t->interpolateGrad(a, 0, 0, 0, g2);
2898 a[0] = 0;
2899 a[1] = 0;
2900 a[2] = 1;
2901 t->interpolateGrad(a, 0, 0, 0, g3);
2902 SVector3 G[3];
2903 G[0] = SVector3(g1[0] + g2[0] - g3[0], g1[1] + g2[1] - g3[1],
2904 g1[2] + g2[2] - g3[2]);
2905 G[1] = SVector3(g2[0] + g3[0] - g1[0], g2[1] + g3[1] - g1[1],
2906 g2[2] + g3[2] - g1[2]);
2907 G[2] = SVector3(g1[0] + g3[0] - g2[0], g1[1] + g3[1] - g2[1],
2908 g1[2] + g3[2] - g2[2]);
2909 SVector3 v10(t->getVertex(1)->x() - t->getVertex(0)->x(),
2910 t->getVertex(1)->y() - t->getVertex(0)->y(),
2911 t->getVertex(1)->z() - t->getVertex(0)->z());
2912 SVector3 v20(t->getVertex(2)->x() - t->getVertex(0)->x(),
2913 t->getVertex(2)->y() - t->getVertex(0)->y(),
2914 t->getVertex(2)->z() - t->getVertex(0)->z());
2915 SVector3 xx = crossprod(v20, v10);
2916 xx.normalize();
2917
2918 double H[3];
2919 for(int k = 0; k < 3; k++) {
2920 auto itk = new2old.find(t->getVertex(k));
2921 if(itk == new2old.end())
2922 myAssembler->getDofValue(t->getVertex(k), 0, 1, H[k]);
2923 else
2924 myAssembler->getDofValue(itk->second, 0, 1, H[k]);
2925 }
2926 double gradH[3];
2927 t->interpolateGrad(H, 0, 0, 0, gradH);
2928
2929 SVector3 temp(gradH[0], gradH[1], gradH[2]);
2930 SVector3 gradHOrtho = crossprod(temp, xx);
2931
2932 double RHS[3] = {dot(gradHOrtho, G[0]), dot(gradHOrtho, G[1]),
2933 dot(gradHOrtho, G[2])};
2934
2935 for(size_t k = 0; k < 3; k++) {
2936 Dof Ek(aaa[t->getEdge(k)], Dof::createTypeWithTwoInts(0, 1));
2937 theta->assemble(Ek, RHS[k] * V);
2938 for(size_t l = 0; l < 3; l++) {
2939 Dof El(aaa[t->getEdge(l)], Dof::createTypeWithTwoInts(0, 1));
2940 theta->assemble(Ek, El, -dot(G[k], G[l]) * V);
2941 }
2942 }
2943 }
2944 }
2945
2946 double SUM = 0.0;
2947 for(size_t i = 0; i < aaa.size(); i++) {
2948 double a;
2949 _lsys->getFromRightHandSide(i, a);
2950 SUM += a;
2951 }
2952 printf("SUM = %12.5E\n", SUM);
2953 SUM /= aaa.size();
2954 for(size_t i = 0; i < aaa.size(); i++) {
2955 // _lsys->addToRightHandSide(i, -SUM);
2956 }
2957
2958 _lsys->systemSolve();
2959 // printScalar(theta, 'T');
2960
2961 double sum = 0;
2962 int count_ = 0;
2963
2964 auto it = aaa.begin();
2965 for(; it != aaa.end(); ++it) {
2966 Dof d(it->second, Dof::createTypeWithTwoInts(0, 1));
2967 double t;
2968 theta->getDofValue(d, t);
2969 MVertex *v0, *v1;
2970 auto it0 = new2old.find(it->first.getVertex(0));
2971 if(it0 == new2old.end())
2972 v0 = it->first.getVertex(0);
2973 else
2974 v0 = it0->second;
2975 it0 = new2old.find(it->first.getVertex(1));
2976 if(it0 == new2old.end())
2977 v1 = it->first.getVertex(1);
2978 else
2979 v1 = it0->second;
2980 MEdge e(v0, v1);
2981 auto itc = C.find(e);
2982 // well... at first ...
2983 itc->second.o_i = SVector3(cos(t), sin(t), 0.0);
2984 // end well
2985 double aa = atan2(dot(itc->second._tgt2, itc->second.o_i),
2986 dot(itc->second._tgt, itc->second.o_i));
2987 itc->second.normalize(aa);
2988 if(!itc->second.inBoundary) { itc->second._a = aa; }
2989 else {
2990 // printf("%12.5E %lu %lu\n",aa,
2991 // itc->second._e.getVertex(0)->getNum(),
2992 // itc->second._e.getVertex(1)->getNum());
2993 itc->second._a = 0;
2994 count_++;
2995 sum += aa;
2996 }
2997 }
2998
2999 sum /= count_;
3000 auto itc = C.begin();
3001 for(; itc != C.end(); ++itc) {
3002 if(!itc->second.inBoundary) {
3003 itc->second._a -= sum;
3004 itc->second._atemp = itc->second._a;
3005 itc->second.normalize(itc->second._a);
3006 }
3007 }
3008
3009 // printCross ("crossCR");
3010 {
3011 // std::string fn = modelName + "_"+"thetaCR.pos";
3012 // FILE *of = fopen(fn.c_str(), "w");
3013 // fprintf(of, "View \"Theta - Crouzeix Raviart\"{\n");
3014 for(size_t i = 0; i < f.size(); i++) {
3015 for(size_t j = 0; j < f[i]->triangles.size(); j++) {
3016 MTriangle *t = f[i]->triangles[j];
3017 Dof d0(aaa[f[i]->triangles[j]->getEdge(0)],
3018 Dof::createTypeWithTwoInts(0, 1));
3019 Dof d1(aaa[f[i]->triangles[j]->getEdge(1)],
3020 Dof::createTypeWithTwoInts(0, 1));
3021 Dof d2(aaa[f[i]->triangles[j]->getEdge(2)],
3022 Dof::createTypeWithTwoInts(0, 1));
3023 double a, b, c;
3024 theta->getDofValue(d0, a);
3025 theta->getDofValue(d1, b);
3026 theta->getDofValue(d2, c);
3027 double A = c + a - b;
3028 double B = a + b - c;
3029 double C = b + c - a;
3030 std::vector<double> ts;
3031 ts.push_back(A);
3032 ts.push_back(B);
3033 ts.push_back(C);
3034 dataTHETA[t->getNum()] = ts;
3035 /* fprintf(of, "ST(%g,%g,%g,%g,%g,%g,%g,%g,%g){%g,%g,%g};\n",
3036 t->getVertex(0)->x(), t->getVertex(0)->y(), t->getVertex(0)->z(),
3037 t->getVertex(1)->x(), t->getVertex(1)->y(), t->getVertex(1)->z(),
3038 t->getVertex(2)->x(), t->getVertex(2)->y(), t->getVertex(2)->z(),
3039 A,B,C);*/
3040 }
3041 }
3042 // fprintf(of, "};\n");
3043 // fclose(of);
3044 }
3045 }
3046 //---------------------------------------------------------------------------
3047
quadLayoutData(GModel * _gm,std::vector<GFace * > & _f,const std::string & name,bool includeFeatureEdges=true)3048 quadLayoutData(GModel *_gm, std::vector<GFace *> &_f, const std::string &name,
3049 bool includeFeatureEdges = true)
3050 : gm(_gm), f(_f), myAssembler(nullptr)
3051 {
3052 modelName = name;
3053 for(size_t i = 0; i < f.size(); i++) {
3054 for(size_t j = 0; j < f[i]->triangles.size(); j++) {
3055 MTriangle *t = f[i]->triangles[j];
3056 for(size_t k = 0; k < 3; k++) {
3057 vs.insert(t->getVertex(k));
3058 MEdge e = t->getEdge(k);
3059 MEdge e1 = t->getEdge((k + 1) % 3);
3060 MEdge e2 = t->getEdge((k + 2) % 3);
3061
3062 // Gaussian Curvatures
3063 MVertex *vk = t->getVertex(k);
3064 MVertex *vk1 = t->getVertex((k + 1) % 3);
3065 MVertex *vk2 = t->getVertex((k + 2) % 3);
3066 SVector3 v1(vk1->x() - vk->x(), vk1->y() - vk->y(),
3067 vk1->z() - vk->z());
3068 SVector3 v2(vk2->x() - vk->x(), vk2->y() - vk->y(),
3069 vk2->z() - vk->z());
3070 double CURV = angle(v1, v2);
3071 auto itg = gaussianCurvatures.find(vk);
3072 if(itg == gaussianCurvatures.end())
3073 gaussianCurvatures[vk] = 2 * M_PI - CURV;
3074 else
3075 itg->second -= CURV;
3076 //---------------------------------------------------------------------
3077
3078 cross2d c(e, t, e1, e2);
3079 auto it = C.find(e);
3080 if(it == C.end())
3081 C.insert(std::make_pair(e, c));
3082 else {
3083 it->second._t.push_back(t);
3084 it->second._neighbors.push_back(e1);
3085 it->second._neighbors.push_back(e2);
3086 }
3087 }
3088 }
3089 }
3090 if(includeFeatureEdges) {
3091 for(size_t i = 0; i < f.size(); i++) {
3092 std::vector<GEdge *> e = f[i]->edges();
3093 for(size_t j = 0; j < e.size(); j++) {
3094 for(size_t k = 0; k < e[j]->lines.size(); k++) {
3095 MLine *l = e[j]->lines[k];
3096 MEdge e = l->getEdge(0);
3097 auto it = C.find(e);
3098 if(it != C.end()) { it->second.inBoundary = true; }
3099 }
3100 }
3101 }
3102 }
3103 auto it = C.begin();
3104 for(; it != C.end(); ++it) it->second.finish(C);
3105 it = C.begin();
3106 for(; it != C.end(); ++it) it->second.finish2();
3107 FILE *F = fopen("gc.pos", "w");
3108 fprintf(F, "View\"\"{\n");
3109 double dd = 0;
3110 for(auto it = gaussianCurvatures.begin(); it != gaussianCurvatures.end();
3111 ++it) {
3112 fprintf(F, "SP(%g,%g,%g){%g};\n", it->first->x(), it->first->y(),
3113 it->first->z(), it->second);
3114 dd += it->second;
3115 }
3116 printf("%22.15E %22.15E\n", dd, dd - 4 * M_PI);
3117 fprintf(F, "};\n");
3118 fclose(F);
3119 }
3120
restoreInitialMesh()3121 void restoreInitialMesh()
3122 {
3123 unDuplicateNodesInCutGraph(f, new2old);
3124 G.clear();
3125 groups.clear();
3126 groups_cg.clear();
3127 cutG.clear();
3128 new2old.clear();
3129 // boundaries.clear();
3130 auto it = C.begin();
3131 for(; it != C.end(); ++it) {
3132 it->second.inCutGraph = false;
3133 it->second._btemp = 0;
3134 }
3135 }
3136
computeUniqueVectorsPerTriangle()3137 int computeUniqueVectorsPerTriangle()
3138 {
3139 // LIFTING
3140 std::set<cross2d *> visited;
3141 while(visited.size() != C.size()) {
3142 for(auto it = C.begin(); it != C.end(); ++it) {
3143 if(visited.find(&(it->second)) == visited.end() &&
3144 cutG.find(it->second._e) == cutG.end()) {
3145 computeLifting(&(it->second), 0, cutG, singularities, boundaries,
3146 visited);
3147 break;
3148 }
3149 }
3150 }
3151 computeUniqueVectorPerTriangle(gm, f, C, 0, d0);
3152 computeUniqueVectorPerTriangle(gm, f, C, 1, d1);
3153 return 0;
3154 }
3155
computeCutGraph(std::map<MEdge,MEdge,MEdgeLessThan> & duplicateEdges)3156 int computeCutGraph(std::map<MEdge, MEdge, MEdgeLessThan> &duplicateEdges)
3157 {
3158 // COMPUTING CUT GRAPH
3159 cutGraph(C, cutG, singularities, boundaries);
3160 for(auto it = C.begin(); it != C.end(); ++it) {
3161 MEdge e0 = it->second._e;
3162 if(cutG.find(e0) != cutG.end()) it->second.inCutGraph = true;
3163 }
3164
3165 groupBoundaries(gm, C, groups, singularities, corners, false);
3166 groupBoundaries(gm, C, groups_cg, singularities, corners, true);
3167
3168 v2t_cont adj;
3169 for(size_t i = 0; i < f.size(); i++) {
3170 buildVertexToElement(f[i]->triangles, adj);
3171 }
3172
3173 std::string fn = modelName + "_groups_analyzed.pos";
3174 FILE *_f = fopen(fn.c_str(), "w");
3175 fprintf(_f, "View \"groups\"{\n");
3176
3177 std::set<MVertex *, MVertexPtrLessThan> isolated_singularities;
3178 {
3179 for(auto it = singularities.begin(); it != singularities.end(); ++it) {
3180 int count = 0;
3181 for(size_t i = 0; i < groups_cg.size(); i++) {
3182 for(size_t k = 0; k < groups_cg[i].size(); k++) {
3183 for(size_t j = 0; j < 2; j++) {
3184 MVertex *v = groups_cg[i][k]->_e.getVertex(j);
3185 if(v == *it) count++;
3186 }
3187 }
3188 }
3189 if(count == 1) { isolated_singularities.insert(*it); }
3190 else {
3191 isolated_singularities.insert(*it);
3192 }
3193 }
3194 }
3195
3196 d0.clear();
3197 d1.clear();
3198 computeUniqueVectorsPerTriangle();
3199
3200 // analyzing groups
3201 {
3202 std::set<MTriangle *> allTrianglesConsidered;
3203 for(size_t i = 0; i < groups_cg.size(); i++) {
3204 groupOfCross2d g(i);
3205 analyzeGroup(groups_cg[i], g, d0, d1, adj, isolated_singularities,
3206 boundaries, allTrianglesConsidered);
3207 g.print(_f);
3208 G.push_back(g);
3209 }
3210 }
3211 fprintf(_f, "};\n");
3212 fclose(_f);
3213
3214 std::multimap<MVertex *, MVertex *, MVertexPtrLessThan> old2new;
3215 duplicateNodesInCutGraph(f, C, new2old, old2new, duplicateEdges,
3216 singularities, adj, G);
3217
3218 for(size_t i = 0; i < groups_cg.size(); i++) {
3219 createJumpyPairs(G[i], singularities, boundaries, old2new);
3220 }
3221 return 0;
3222 }
3223
computeCrossFieldAndH(std::map<MVertex *,int> * s,std::map<int,std::vector<double>> & dataTHETA)3224 int computeCrossFieldAndH(std::map<MVertex *, int> *s,
3225 std::map<int, std::vector<double> > &dataTHETA)
3226 {
3227 double a = TimeOfDay();
3228 computeHFromSingularities(*s);
3229 double b = TimeOfDay();
3230 Msg::Info("Real part H (nodal) has been computed in %4g sec", b - a);
3231
3232 std::map<MEdge, MEdge, MEdgeLessThan> duplicateEdges;
3233
3234 double c = TimeOfDay();
3235 computeCutGraph(duplicateEdges);
3236 Msg::Info("Cut Graph has been computed in %4g sec", c - b);
3237
3238 double d = TimeOfDay();
3239 computeThetaUsingHCrouzeixRaviart(dataTHETA);
3240 Msg::Info("Imaginary part H (crouzeix raviart/multi-valued) has been "
3241 "computed in %4g sec",
3242 d - c);
3243 restoreInitialMesh();
3244 return 0;
3245 }
3246
intersectEdgeEdge(MEdge & e,MVertex * v1,MVertex * v2,GFace * gf)3247 MVertex *intersectEdgeEdge(MEdge &e, MVertex *v1, MVertex *v2, GFace *gf)
3248 {
3249 MVertex *e1 = e.getVertex(0);
3250 MVertex *e2 = e.getVertex(1);
3251
3252 if(e1 == v2) return nullptr;
3253 if(e2 == v2) return nullptr;
3254
3255 SVector3 e1e2(e2->x() - e1->x(), e2->y() - e1->y(), e2->z() - e1->z());
3256 SVector3 e1v1(v1->x() - e1->x(), v1->y() - e1->y(), v1->z() - e1->z());
3257 SVector3 e2v1(v1->x() - e2->x(), v1->y() - e2->y(), v1->z() - e2->z());
3258
3259 SVector3 a = crossprod(e1e2, e1v1);
3260 double b = dot(e1v1, e2v1);
3261 if(a.norm() < 1.e-10 && b < 0) return v1;
3262
3263 if(!v2) {
3264 // Msg::Error("Error In CutMesh");
3265 return nullptr;
3266 }
3267
3268 SVector3 e2v2(v2->x() - e2->x(), v2->y() - e2->y(), v2->z() - e2->z());
3269 SVector3 e1v2(v2->x() - e1->x(), v2->y() - e1->y(), v2->z() - e1->z());
3270 a = crossprod(e1e2, e1v2);
3271 b = dot(e1v2, e2v2);
3272 if(a.norm() < 1.e-10 && b < 0) return v2;
3273
3274 double x[2];
3275
3276 bool inters = intersection_segments(
3277 SPoint3(e.getVertex(0)->x(), e.getVertex(0)->y(), e.getVertex(0)->z()),
3278 SPoint3(e.getVertex(1)->x(), e.getVertex(1)->y(), e.getVertex(1)->z()),
3279 SPoint3(v1->x(), v1->y(), v1->z()), SPoint3(v2->x(), v2->y(), v2->z()),
3280 x);
3281 if(!inters) return nullptr;
3282 return new MEdgeVertex(v2->x() * x[1] + v1->x() * (1. - x[1]),
3283 v2->y() * x[1] + v1->y() * (1. - x[1]),
3284 v2->z() * x[1] + v1->z() * (1. - x[1]), nullptr, 0);
3285 }
3286
cut_edge(std::map<MEdge,int,MEdgeLessThan> & ecuts,MVertex * v0,MVertex * v1,MVertex * mid)3287 void cut_edge(std::map<MEdge, int, MEdgeLessThan> &ecuts, MVertex *v0,
3288 MVertex *v1, MVertex *mid)
3289 {
3290 MEdge e(v0, v1);
3291 auto it = ecuts.find(e);
3292 if(it != ecuts.end()) {
3293 int index = it->second;
3294 ecuts.erase(it);
3295 MEdge e1(v0, mid);
3296 MEdge e2(mid, v1);
3297 ecuts[e1] = index;
3298 ecuts[e2] = index;
3299 }
3300 }
3301
cutTriangles(std::vector<MTriangle * > & ts,GFace * gf,MVertex * v1,MVertex * v2,GEdge * ge,std::map<MEdge,int,MEdgeLessThan> & ecuts,int index,FILE * f)3302 void cutTriangles(std::vector<MTriangle *> &ts, GFace *gf, MVertex *v1,
3303 MVertex *v2, GEdge *ge,
3304 std::map<MEdge, int, MEdgeLessThan> &ecuts, int index,
3305 FILE *f)
3306 {
3307 std::map<MEdge, MVertex *, MEdgeLessThan> e_cut;
3308 std::vector<MTriangle *> newt;
3309
3310 for(size_t i = 0; i < ts.size(); ++i) {
3311 MVertex *vs[3] = {nullptr, nullptr, nullptr};
3312 for(size_t j = 0; j < 3; ++j) {
3313 MEdge e = ts[i]->getEdge(j);
3314 if(e_cut.find(e) == e_cut.end()) {
3315 MVertex *v = intersectEdgeEdge(e, v1, v2, gf);
3316 if(v && v != v1 && v != v2) {
3317 gf->mesh_vertices.push_back(v);
3318 if(f)
3319 fprintf(f, "SP(%g,%g,%g){%d};\n", v->x(), v->y(), v->z(),
3320 gf->tag());
3321 }
3322 e_cut[e] = v;
3323 }
3324 vs[j] = e_cut[e];
3325 }
3326
3327 if(!vs[0] && !vs[1] && !vs[2])
3328 newt.push_back(ts[i]);
3329 else if(vs[0] && !vs[1] && !vs[2]) {
3330 // printf("\n");
3331 newt.push_back(
3332 new MTriangle(ts[i]->getVertex(0), vs[0], ts[i]->getVertex(2)));
3333 newt.push_back(
3334 new MTriangle(vs[0], ts[i]->getVertex(1), ts[i]->getVertex(2)));
3335 cut_edge(ecuts, ts[i]->getVertex(0), ts[i]->getVertex(1), vs[0]);
3336 MEdge ed(ts[i]->getVertex(2), vs[0]);
3337 ecuts[ed] = index;
3338 }
3339 else if(!vs[0] && vs[1] && !vs[2]) {
3340 // printf("coucouc\n");
3341 newt.push_back(
3342 new MTriangle(ts[i]->getVertex(0), ts[i]->getVertex(1), vs[1]));
3343 newt.push_back(
3344 new MTriangle(ts[i]->getVertex(0), vs[1], ts[i]->getVertex(2)));
3345 cut_edge(ecuts, ts[i]->getVertex(2), ts[i]->getVertex(1), vs[1]);
3346 MEdge ed(ts[i]->getVertex(0), vs[1]);
3347 ecuts[ed] = index;
3348 }
3349 else if(!vs[0] && !vs[1] && vs[2]) {
3350 // printf("coucouc\n");
3351 newt.push_back(
3352 new MTriangle(ts[i]->getVertex(0), ts[i]->getVertex(1), vs[2]));
3353 newt.push_back(
3354 new MTriangle(vs[2], ts[i]->getVertex(1), ts[i]->getVertex(2)));
3355 cut_edge(ecuts, ts[i]->getVertex(2), ts[i]->getVertex(0), vs[2]);
3356 MEdge ed(ts[i]->getVertex(1), vs[2]);
3357 ecuts[ed] = index;
3358 }
3359 else if(vs[0] && vs[1] && !vs[2]) { // OK
3360 newt.push_back(new MTriangle(ts[i]->getVertex(0), vs[0], vs[1]));
3361 newt.push_back(new MTriangle(ts[i]->getVertex(1), vs[1], vs[0]));
3362 newt.push_back(
3363 new MTriangle(ts[i]->getVertex(0), vs[1], ts[i]->getVertex(2)));
3364 cut_edge(ecuts, ts[i]->getVertex(0), ts[i]->getVertex(1), vs[0]);
3365 cut_edge(ecuts, ts[i]->getVertex(2), ts[i]->getVertex(1), vs[1]);
3366 MEdge ed(vs[0], vs[1]);
3367 ecuts[ed] = index;
3368 }
3369 else if(vs[0] && !vs[1] && vs[2]) { // OK
3370 newt.push_back(new MTriangle(ts[i]->getVertex(0), vs[0], vs[2]));
3371 newt.push_back(new MTriangle(ts[i]->getVertex(2), vs[2], vs[0]));
3372 newt.push_back(
3373 new MTriangle(ts[i]->getVertex(2), vs[0], ts[i]->getVertex(1)));
3374 cut_edge(ecuts, ts[i]->getVertex(0), ts[i]->getVertex(1), vs[0]);
3375 cut_edge(ecuts, ts[i]->getVertex(2), ts[i]->getVertex(0), vs[2]);
3376 MEdge ed(vs[0], vs[2]);
3377 ecuts[ed] = index;
3378 }
3379 else if(!vs[0] && vs[1] && vs[2]) {
3380 newt.push_back(new MTriangle(ts[i]->getVertex(2), vs[2], vs[1]));
3381 newt.push_back(new MTriangle(ts[i]->getVertex(0), vs[1], vs[2]));
3382 newt.push_back(
3383 new MTriangle(ts[i]->getVertex(1), vs[1], ts[i]->getVertex(0)));
3384 cut_edge(ecuts, ts[i]->getVertex(1), ts[i]->getVertex(2), vs[1]);
3385 cut_edge(ecuts, ts[i]->getVertex(2), ts[i]->getVertex(0), vs[2]);
3386 MEdge ed(vs[1], vs[2]);
3387 ecuts[ed] = index;
3388 }
3389 else if(vs[0] && vs[1] && vs[2]) {
3390 newt.push_back(new MTriangle(vs[0], vs[1], vs[2]));
3391 newt.push_back(new MTriangle(ts[i]->getVertex(0), vs[0], vs[2]));
3392 newt.push_back(new MTriangle(ts[i]->getVertex(1), vs[1], vs[0]));
3393 newt.push_back(new MTriangle(ts[i]->getVertex(2), vs[2], vs[1]));
3394 }
3395 else {
3396 newt.push_back(ts[i]);
3397 }
3398 }
3399 // printf("replcing %lu by %lu\n",ts.size(),newt.size());
3400 ts = newt;
3401 }
3402
cutMesh(std::map<MEdge,edgeCuts,MEdgeLessThan> & cuts)3403 void cutMesh(std::map<MEdge, edgeCuts, MEdgeLessThan> &cuts)
3404 {
3405 auto it = cuts.begin();
3406 std::map<MEdge, int, MEdgeLessThan> ecuts;
3407
3408 FILE *F = fopen("addedpoints.pos", "w");
3409 fprintf(F, "View \"\"{\n");
3410 for(; it != cuts.end(); ++it) { it->second.finish(gm, F); }
3411
3412 for(size_t i = 0; i < f.size(); i++) {
3413 std::vector<MTriangle *> newT;
3414 for(size_t j = 0; j < f[i]->triangles.size(); j++) {
3415 std::set<int> indices;
3416 std::multimap<int, std::pair<MVertex *, std::pair<int, int> > > tcuts;
3417
3418 for(size_t k = 0; k < 3; k++) {
3419 MEdge e = f[i]->triangles[j]->getEdge(k);
3420 auto it = cuts.find(e);
3421 if(it != cuts.end()) {
3422 for(size_t l = 0; l < it->second.vs.size(); ++l) {
3423 std::pair<int, std::pair<MVertex *, std::pair<int, int> > > pp =
3424 std::make_pair(
3425 it->second.indexOfCuts[l],
3426 std::make_pair(it->second.vs[l],
3427 std::make_pair(k, it->second.idsOfCuts[l])));
3428 tcuts.insert(pp);
3429 indices.insert(it->second.indexOfCuts[l]);
3430 }
3431 }
3432 }
3433
3434 auto iti = indices.begin();
3435 std::vector<MTriangle *> ttt;
3436 ttt.push_back(f[i]->triangles[j]);
3437 for(; iti != indices.end(); ++iti) {
3438 // if (*iti != 313310)continue;
3439 GEdge *ge = gm->getEdgeByTag(*iti);
3440 if(tcuts.count(*iti) == 1) {
3441 auto itt = tcuts.lower_bound(*iti);
3442 MVertex *v0 = itt->second.first;
3443 int k = itt->second.second.first;
3444 cutTriangles(ttt, f[i], v0,
3445 f[i]->triangles[j]->getVertex((k + 2) % 3), ge, ecuts,
3446 *iti, F);
3447 }
3448 else if(tcuts.count(*iti) == 2) {
3449 auto itt = tcuts.lower_bound(*iti);
3450 MVertex *v0 = itt->second.first;
3451 ++itt;
3452 MVertex *v1 = itt->second.first;
3453 cutTriangles(ttt, f[i], v0, v1, ge, ecuts, *iti, F);
3454 }
3455 else if(tcuts.count(*iti) == 3) {
3456 auto itt = tcuts.lower_bound(*iti);
3457 int k0 = itt->second.second.first;
3458 int id0 = itt->second.second.second;
3459 MVertex *v0 = itt->second.first;
3460 ++itt;
3461 int k1 = itt->second.second.first;
3462 int id1 = itt->second.second.second;
3463 MVertex *v1 = itt->second.first;
3464 ++itt;
3465 int k2 = itt->second.second.first;
3466 int id2 = itt->second.second.second;
3467 MVertex *v2 = itt->second.first;
3468 ;
3469
3470 if(abs(id0 - id1) <= 2) {
3471 cutTriangles(ttt, f[i], v2,
3472 f[i]->triangles[j]->getVertex((k2 + 2) % 3), ge,
3473 ecuts, *iti, F);
3474 cutTriangles(ttt, f[i], v0, v1, ge, ecuts, *iti, F);
3475 }
3476 else if(abs(id0 - id2) <= 2) {
3477 cutTriangles(ttt, f[i], v1,
3478 f[i]->triangles[j]->getVertex((k1 + 2) % 3), ge,
3479 ecuts, *iti, F);
3480 cutTriangles(ttt, f[i], v0, v2, ge, ecuts, *iti, F);
3481 }
3482 else if(abs(id1 - id2) <= 2) {
3483 cutTriangles(ttt, f[i], v0,
3484 f[i]->triangles[j]->getVertex((k0 + 2) % 3), ge,
3485 ecuts, *iti, F);
3486 cutTriangles(ttt, f[i], v1, v2, ge, ecuts, *iti, F);
3487 }
3488 else {
3489 printf("BAD BEHAVIOR 3\n");
3490 printf("%lu %lu %lu \n", v0->getNum(), v1->getNum(),
3491 v2->getNum());
3492 printf("%d %d %d\n", id0, id1, id2);
3493 }
3494 }
3495 else if(tcuts.count(*iti) == 4) {
3496 auto itt = tcuts.lower_bound(*iti);
3497 int id0 = itt->second.second.second;
3498 MVertex *v0 = itt->second.first;
3499 ++itt;
3500 int id1 = itt->second.second.second;
3501 MVertex *v1 = itt->second.first;
3502 ++itt;
3503 int id2 = itt->second.second.second;
3504 MVertex *v2 = itt->second.first;
3505 ++itt;
3506 int id3 = itt->second.second.second;
3507 MVertex *v3 = itt->second.first;
3508 if(abs(id0 - id1) <= 2 || abs(id2 - id3) <= 2) {
3509 cutTriangles(ttt, f[i], v0, v1, ge, ecuts, *iti, F);
3510 cutTriangles(ttt, f[i], v2, v3, ge, ecuts, *iti, F);
3511 }
3512 else if(abs(id0 - id2) <= 2 || abs(id1 - id3) <= 2) {
3513 cutTriangles(ttt, f[i], v0, v2, ge, ecuts, *iti, F);
3514 cutTriangles(ttt, f[i], v1, v3, ge, ecuts, *iti, F);
3515 }
3516 else if(abs(id0 - id3) <= 2 || abs(id1 - id2) <= 2) {
3517 cutTriangles(ttt, f[i], v0, v3, ge, ecuts, *iti, F);
3518 cutTriangles(ttt, f[i], v1, v2, ge, ecuts, *iti, F);
3519 }
3520 else {
3521 printf("%d %d %d %d\n", id0, id1, id2, id3);
3522 }
3523 }
3524 else if(tcuts.count(*iti) == 6) {
3525 auto itt = tcuts.lower_bound(*iti);
3526 std::pair<int, MVertex *> id[10];
3527 for(std::size_t kk = 0; kk < tcuts.count(*iti); kk++) {
3528 id[kk] =
3529 std::make_pair(itt->second.second.second, itt->second.first);
3530 ++itt;
3531 }
3532 std::sort(id, id + 6);
3533 cutTriangles(ttt, f[i], id[0].second, id[1].second, ge, ecuts, *iti,
3534 F);
3535 cutTriangles(ttt, f[i], id[2].second, id[3].second, ge, ecuts, *iti,
3536 F);
3537 cutTriangles(ttt, f[i], id[4].second, id[5].second, ge, ecuts, *iti,
3538 F);
3539 printf("%d %d %d %d %d %d\n", id[0].first, id[1].first, id[2].first,
3540 id[3].first, id[4].first, id[5].first);
3541 }
3542 else {
3543 Msg::Error("TODO %lu in CutMesh !!!!!!!", tcuts.count(*iti));
3544 }
3545 }
3546 newT.insert(newT.begin(), ttt.begin(), ttt.end());
3547 }
3548 f[i]->triangles = newT;
3549 }
3550
3551 fprintf(F, "};\n");
3552 fclose(F);
3553
3554 F = fopen("edges.pos", "w");
3555 fprintf(F, "View \"\"{\n");
3556
3557 for(auto it = ecuts.begin(); it != ecuts.end(); ++it) {
3558 GEdge *ge = gm->getEdgeByTag(it->second);
3559 ge->lines.push_back(
3560 new MLine(it->first.getVertex(0), it->first.getVertex(1)));
3561 fprintf(F, "SL(%g,%g,%g,%g,%g,%g){1,1};\n", it->first.getVertex(0)->x(),
3562 it->first.getVertex(0)->y(), it->first.getVertex(0)->z(),
3563 it->first.getVertex(1)->x(), it->first.getVertex(1)->y(),
3564 it->first.getVertex(1)->z());
3565 for(int i = 0; i < 2; i++) {
3566 if(std::find(ge->mesh_vertices.begin(), ge->mesh_vertices.end(),
3567 it->first.getVertex(i)) == ge->mesh_vertices.end()) {
3568 ge->mesh_vertices.push_back(it->first.getVertex(i));
3569 it->first.getVertex(i)->setEntity(ge);
3570 }
3571 }
3572 }
3573 fprintf(F, "};\n");
3574 fclose(F);
3575 CTX::instance()->mesh.changed = ENT_ALL;
3576 }
3577
correctionOnCutGraph(std::map<MEdge,edgeCuts,MEdgeLessThan> & cuts,std::map<MVertex *,MVertex *,MVertexPtrLessThan> & new2old)3578 int correctionOnCutGraph(
3579 std::map<MEdge, edgeCuts, MEdgeLessThan> &cuts,
3580 std::map<MVertex *, MVertex *, MVertexPtrLessThan> &new2old)
3581 {
3582 std::map<MEdge, MEdge, MEdgeLessThan> duplicateEdges;
3583
3584 for(auto it = cuts.begin(); it != cuts.end(); ++it) {
3585 MVertex *v0 = it->first.getVertex(0);
3586 MVertex *v1 = it->first.getVertex(1);
3587 MVertex *v2 = new2old.find(v0) == new2old.end() ? v0 : new2old[v0];
3588 MVertex *v3 = new2old.find(v1) == new2old.end() ? v1 : new2old[v1];
3589 MEdge e0(v0, v1);
3590 MEdge e1(v2, v3);
3591 duplicateEdges[e0] = e1;
3592 }
3593
3594 for(auto it2 = duplicateEdges.begin(); it2 != duplicateEdges.end(); ++it2) {
3595 auto it3 = cuts.find(it2->first);
3596 auto it4 = cuts.find(it2->second);
3597 if(it3 != cuts.end() && it4 == cuts.end()) {
3598 cuts[it2->second] = it3->second;
3599 }
3600 else if(it4 != cuts.end() && it3 == cuts.end()) {
3601 cuts[it2->first] = it4->second;
3602 }
3603 else if(it4 != cuts.end() && it3 != cuts.end()) {
3604 it4->second.vs.insert(it4->second.vs.begin(), it3->second.vs.begin(),
3605 it3->second.vs.end());
3606 it4->second.indexOfCuts.insert(it4->second.indexOfCuts.begin(),
3607 it3->second.indexOfCuts.begin(),
3608 it3->second.indexOfCuts.end());
3609 it3->second.vs.insert(it3->second.vs.begin(), it4->second.vs.begin(),
3610 it4->second.vs.end());
3611 it3->second.indexOfCuts.insert(it3->second.indexOfCuts.begin(),
3612 it4->second.indexOfCuts.begin(),
3613 it4->second.indexOfCuts.end());
3614 }
3615 }
3616 return 0;
3617 }
3618
computeQuadLayout(std::map<MVertex *,double> & potU,std::map<MVertex *,double> & potV,std::map<MEdge,MEdge,MEdgeLessThan> & duplicateEdges,std::map<MEdge,edgeCuts,MEdgeLessThan> & cuts)3619 int computeQuadLayout(std::map<MVertex *, double> &potU,
3620 std::map<MVertex *, double> &potV,
3621 std::map<MEdge, MEdge, MEdgeLessThan> &duplicateEdges,
3622 std::map<MEdge, edgeCuts, MEdgeLessThan> &cuts)
3623 {
3624 std::vector<cutGraphPassage> passages;
3625 while(1) {
3626 computePotential(gm, f, *myAssembler, C, new2old, groups, duplicateEdges,
3627 d0, d1, G, potU, potV, passages);
3628
3629 if(computeIsos(gm, f, singularities, C, new2old, duplicateEdges, groups,
3630 groups_cg, potU, potV, cutG, G, cuts, passages) == true) {
3631 printf("COMPUTE ISOS DONE\n");
3632 break;
3633 }
3634 break;
3635 }
3636
3637 correctionOnCutGraph(cuts, new2old);
3638
3639 double MAXX = 0.;
3640 // double SUM_LEFT = 0 , SUM_RIGHT = 0;
3641 for(size_t i = 0; i < groups_cg.size(); i++) {
3642 double MAXD1 = -1.e22, MIND1 = 1.e22, MAXD2 = -1.e22, MIND2 = 1.e22;
3643 for(size_t j = 0; j < G[i].left.size(); j++) {
3644 /* if (G[i].groupId == 5 || G[i].groupId == 6){
3645 printf("%lu %lu %12.5E %12.5E %12.5E
3646 %12.5E\n",G[i].left[j]->getNum(),G[i].right[j]->getNum(),
3647 potU[G[i].left[j]],potU[G[i].right[j]],
3648 potV[G[i].left[j]],potV[G[i].right[j]]);
3649 }
3650 */
3651 double Ul = potU[G[i].left[j]];
3652 double Ur = potU[G[i].right[j]];
3653 double Vl = potV[G[i].left[j]];
3654 double Vr = potV[G[i].right[j]];
3655 double D1 = Ul - G[i].mat[0][0] * Ur - G[i].mat[0][1] * Vr;
3656 double D2 = Vl - G[i].mat[1][0] * Ur - G[i].mat[1][1] * Vr;
3657 MAXD1 = std::max(D1, MAXD1);
3658 MAXD2 = std::max(D2, MAXD2);
3659 MIND1 = std::min(D1, MIND1);
3660 MIND2 = std::min(D2, MIND2);
3661 }
3662 // SUM_LEFT += MAXD1;
3663 // SUM_RIGHT += MAXD2;
3664 // Dof E1(g.left[]->getNum(),
3665 // Dof::createTypeWithTwoInts(0, 3 + 100 * G[i].groupId));
3666
3667 Msg::Debug("group %3d ROT (%12.5E %12.5E) (%12.5E %12.5E)", G[i].groupId,
3668 G[i].mat[0][0], G[i].mat[0][1], G[i].mat[1][0],
3669 G[i].mat[1][1]);
3670
3671 Msg::Debug("group %3d DA(%12.5E %12.5E %12.5E) D2(%12.5E %12.5E %12.5E)",
3672 G[i].groupId, MAXD1, MIND1, MAXD1 - MIND1, MAXD2, MIND2,
3673 MAXD2 - MIND2);
3674 MAXX = std::max(MAXD2 - MIND2, MAXX);
3675 }
3676 if(MAXX < 1.e-09)
3677 Msg::Info("Success in computing potentials (all jumps are OK)");
3678 else
3679 Msg::Warning("Quad Layout Failure");
3680 return 0;
3681 }
3682 };
3683
findPhysicalGroupsForSingularities(GModel * gm,std::vector<GFace * > & f,std::map<MVertex *,int> & temp)3684 static void findPhysicalGroupsForSingularities(GModel *gm,
3685 std::vector<GFace *> &f,
3686 std::map<MVertex *, int> &temp)
3687 {
3688 std::map<int, std::vector<GEntity *> > groups[4];
3689 gm->getPhysicalGroups(groups);
3690 for(auto it = groups[0].begin(); it != groups[0].end(); ++it) {
3691 std::string name = gm->getPhysicalName(0, it->first);
3692 if(name == "SINGULARITY_OF_INDEX_THREE") {
3693 for(size_t j = 0; j < it->second.size(); j++) {
3694 if(!it->second[j]->mesh_vertices.empty())
3695 temp[it->second[j]->mesh_vertices[0]] = 1;
3696 }
3697 }
3698 else if(name == "SINGULARITY_OF_INDEX_FIVE") {
3699 for(size_t j = 0; j < it->second.size(); j++) {
3700 if(!it->second[j]->mesh_vertices.empty())
3701 temp[it->second[j]->mesh_vertices[0]] = -1;
3702 }
3703 }
3704 else if(name == "SINGULARITY_OF_INDEX_SIX") {
3705 for(size_t j = 0; j < it->second.size(); j++) {
3706 if(!it->second[j]->mesh_vertices.empty())
3707 temp[it->second[j]->mesh_vertices[0]] = -2;
3708 }
3709 }
3710 else if(name == "SINGULARITY_OF_INDEX_EIGHT") {
3711 for(size_t j = 0; j < it->second.size(); j++) {
3712 if(!it->second[j]->mesh_vertices.empty())
3713 temp[it->second[j]->mesh_vertices[0]] = -4;
3714 }
3715 }
3716 else if(name == "SINGULARITY_OF_INDEX_TWO") {
3717 for(size_t j = 0; j < it->second.size(); j++) {
3718 if(!it->second[j]->mesh_vertices.empty())
3719 temp[it->second[j]->mesh_vertices[0]] = 2;
3720 }
3721 }
3722 }
3723 }
3724
computeCrossFieldAndH(GModel * gm,std::vector<GFace * > & f,std::vector<int> & tags,bool layout=true)3725 static int computeCrossFieldAndH(GModel *gm, std::vector<GFace *> &f,
3726 std::vector<int> &tags, bool layout = true)
3727 {
3728 quadLayoutData qLayout(gm, f, gm->getName());
3729 std::map<MVertex *, int> temp;
3730 std::map<int, std::vector<double> > dataH;
3731 std::map<int, std::vector<double> > dataTHETA;
3732 std::map<int, std::vector<double> > dataDir;
3733 std::map<int, std::vector<double> > dataDirOrtho;
3734 std::map<int, std::vector<double> > dataU;
3735 std::map<int, std::vector<double> > dataV;
3736 std::map<MVertex *, double> potU, potV;
3737 findPhysicalGroupsForSingularities(gm, f, temp);
3738 std::map<MEdge, MEdge, MEdgeLessThan> duplicateEdges;
3739 std::map<MEdge, edgeCuts, MEdgeLessThan> cuts;
3740 if(temp.size()) {
3741 Msg::Info("Computing cross field from %d prescribed singularities",
3742 temp.size());
3743 qLayout.computeCrossFieldAndH(&temp, dataTHETA);
3744 qLayout.computeCutGraph(duplicateEdges);
3745 }
3746 else {
3747 Msg::Info("Computing a cross field");
3748 qLayout.computeCrossFieldAndH();
3749 qLayout.computeCutGraph(duplicateEdges);
3750 qLayout.computeThetaUsingHCrouzeixRaviart(dataTHETA);
3751 }
3752 if(layout) { qLayout.computeQuadLayout(potU, potV, duplicateEdges, cuts); }
3753
3754 PViewDataGModel *d = new PViewDataGModel;
3755 PViewDataGModel *dt = new PViewDataGModel(PViewDataGModel::ElementNodeData);
3756 PViewDataGModel *dd = new PViewDataGModel(PViewDataGModel::ElementData);
3757 std::string name = gm->getName() + "_H";
3758 d->setName(name);
3759 d->setFileName(name + ".msh");
3760 name = gm->getName() + "_Theta";
3761 dt->setName(name);
3762 dt->setFileName(name + ".msh");
3763 name = gm->getName() + "_Directions";
3764 dd->setName(name);
3765 dd->setFileName(name + ".msh");
3766 PViewDataGModel *U = nullptr;
3767 PViewDataGModel *V = nullptr;
3768 if(layout) {
3769 U = new PViewDataGModel(PViewDataGModel::ElementNodeData);
3770 V = new PViewDataGModel(PViewDataGModel::ElementNodeData);
3771 name = gm->getName() + "_U";
3772 U->setName(name);
3773 U->setFileName(name + ".msh");
3774 name = gm->getName() + "_V";
3775 V->setName(name);
3776 V->setFileName(name + ".msh");
3777
3778 for(size_t i = 0; i < f.size(); i++) {
3779 for(size_t j = 0; j < f[i]->triangles.size(); j++) {
3780 MTriangle *t = f[i]->triangles[j];
3781 double a = potU[f[i]->triangles[j]->getVertex(0)];
3782 double b = potU[f[i]->triangles[j]->getVertex(1)];
3783 double c = potU[f[i]->triangles[j]->getVertex(2)];
3784 std::vector<double> ts;
3785 ts.push_back(a);
3786 ts.push_back(b);
3787 ts.push_back(c);
3788 dataU[t->getNum()] = ts;
3789 a = potV[f[i]->triangles[j]->getVertex(0)];
3790 b = potV[f[i]->triangles[j]->getVertex(1)];
3791 c = potV[f[i]->triangles[j]->getVertex(2)];
3792 ts.clear();
3793 ts.push_back(a);
3794 ts.push_back(b);
3795 ts.push_back(c);
3796 dataV[t->getNum()] = ts;
3797 }
3798 }
3799
3800 U->addData(gm, dataU, 0, 0.0, 1, 1);
3801 U->finalize();
3802 V->addData(gm, dataV, 0, 0.0, 1, 1);
3803 V->finalize();
3804 }
3805 for(auto it = qLayout.d0.begin(); it != qLayout.d0.end(); ++it) {
3806 std::vector<double> jj;
3807 jj.push_back(it->second.x());
3808 jj.push_back(it->second.y());
3809 jj.push_back(it->second.z());
3810 dataDir[it->first->getNum()] = jj;
3811 }
3812 for(auto it = qLayout.d1.begin(); it != qLayout.d1.end(); ++it) {
3813 std::vector<double> jj;
3814 jj.push_back(it->second.x());
3815 jj.push_back(it->second.y());
3816 jj.push_back(it->second.z());
3817 dataDirOrtho[it->first->getNum()] = jj;
3818 }
3819 for(auto it = qLayout.vs.begin(); it != qLayout.vs.end(); ++it) {
3820 double h;
3821 qLayout.myAssembler->getDofValue(*it, 0, 1, h);
3822 std::vector<double> jj;
3823 jj.push_back(h);
3824 // printf("adding data for %lu\n",(*it)->getNum());
3825 dataH[(*it)->getNum()] = jj;
3826 }
3827
3828 d->addData(gm, dataH, 0, 0.0, 1, 1);
3829 d->finalize();
3830 dt->addData(gm, dataTHETA, 0, 0.0, 1, 1);
3831 dt->finalize();
3832 dd->addData(gm, dataDir, 0, 0.0, 1, 3);
3833 dd->addData(gm, dataDirOrtho, 1, 0.0, 1, 3);
3834 dd->finalize();
3835
3836 std::string posout = gm->getName() + "_QLayoutResults.pos";
3837
3838 qLayout.restoreInitialMesh();
3839 dt->writePOS(posout, false, true, true);
3840 dd->writePOS(posout, false, true, true);
3841 d->writePOS(posout, false, true, true);
3842 if(layout) {
3843 U->writePOS(posout, false, true, true);
3844 V->writePOS(posout, false, true, true);
3845 }
3846 // return 0;
3847 Msg::Info("Cutting the mesh");
3848
3849 qLayout.cutMesh(cuts);
3850
3851 Msg::Info("Classifying the model");
3852 discreteEdge *de = new discreteEdge(
3853 GModel::current(), GModel::current()->getMaxElementaryNumber(1) + 1,
3854 nullptr, nullptr);
3855 GModel::current()->add(de);
3856 computeNonManifoldEdges(GModel::current(), de->lines, true);
3857 classifyFaces(GModel::current(), M_PI / 4 * .999);
3858 GModel::current()->remove(de);
3859 // delete de;
3860 GModel::current()->pruneMeshVertexAssociations();
3861
3862 std::string mshout = gm->getName() + "_Cut.msh";
3863 gm->writeMSH(mshout, 4.0, false, true);
3864
3865 int countError = 0;
3866 for(auto it = GModel::current()->firstFace();
3867 it != GModel::current()->lastFace(); it++) {
3868 if((*it)->edges().size() != 4) {
3869 Msg::Warning("quad layout failed : face %lu has %lu boundaries",
3870 (*it)->tag(), (*it)->edges().size());
3871 countError++;
3872 }
3873 }
3874 if(!countError) {
3875 Msg::Info(
3876 "Quad layout success : the model is partitioned in %d master quads",
3877 GModel::current()->getNumFaces());
3878 Msg::Info("Partitioned mesh has been saved in %s", mshout.c_str());
3879 Msg::Info("Result of computations have been saved in %s", posout.c_str());
3880 }
3881 delete d;
3882 delete dd;
3883 delete dt;
3884 if(layout) {
3885 delete U;
3886 delete V;
3887 }
3888
3889 return 0;
3890 }
3891
3892 #endif
3893
getFacesOfTheModel(GModel * gm,std::vector<GFace * > & f)3894 static void getFacesOfTheModel(GModel *gm, std::vector<GFace *> &f)
3895 {
3896 for(auto it = gm->firstFace(); it != gm->lastFace(); ++it) {
3897 GFace *gf = *it;
3898 f.push_back(gf);
3899 }
3900 }
3901
computeCrossFieldAndH(GModel * gm,std::vector<int> & tags)3902 int computeCrossFieldAndH(GModel *gm, std::vector<int> &tags)
3903 {
3904 std::vector<GFace *> f;
3905 getFacesOfTheModel(gm, f);
3906 #if defined(HAVE_SOLVER) && defined(HAVE_POST)
3907 return computeCrossFieldAndH(gm, f, tags);
3908 #else
3909 Msg::Error("Cross field computation requires solver and post-pro module");
3910 return -1;
3911 #endif
3912 }
3913
computeQuadLayout(GModel * gm,std::vector<int> & tags)3914 int computeQuadLayout(GModel *gm, std::vector<int> &tags)
3915 {
3916 std::vector<GFace *> f;
3917 getFacesOfTheModel(gm, f);
3918
3919 #if defined(HAVE_SOLVER) && defined(HAVE_POST)
3920 return computeCrossFieldAndH(gm, f, tags, true);
3921 #else
3922 Msg::Error("Cross field computation requires solver and post-pro module");
3923 return -1;
3924 #endif
3925 }
3926
computeCrossField(GModel * gm,std::vector<int> & tags)3927 int computeCrossField(GModel *gm, std::vector<int> &tags)
3928 {
3929 std::vector<GFace *> f;
3930 getFacesOfTheModel(gm, f);
3931
3932 #if defined(HAVE_SOLVER) && defined(HAVE_POST)
3933 return computeCrossFieldAndH(gm, f, tags, true);
3934 // return computeQuadLayout(gm, f);
3935 #else
3936 Msg::Error("Cross field computation requires solver and post-pro module");
3937 return -1;
3938 #endif
3939 }
3940