1 // -*- Mode : c++ -*- 2 // 3 // SUMMARY : 4 // USAGE : 5 // ORG : 6 // AUTHOR : Frederic Hecht 7 // E-MAIL : hecht@ann.jussieu.fr 8 // 9 // ********** DO NOT REMOVE THIS BANNER ********** 10 /* 11 12 This file is part of Freefem++ 13 14 Freefem++ is free software; you can redistribute it and/or modify 15 it under the terms of the GNU Lesser General Public License as published by 16 the Free Software Foundation; either version 2.1 of the License, or 17 (at your option) any later version. 18 19 Freefem++ is distributed in the hope that it will be useful, 20 but WITHOUT ANY WARRANTY; without even the implied warranty of 21 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 22 GNU Lesser General Public License for more details. 23 24 You should have received a copy of the GNU Lesser General Public License 25 along with Freefem++; if not, write to the Free Software 26 Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA 27 */ 28 29 // 30 // SUMMARY: Bamg: Bidimensional Anisotrope Mesh Generator 31 // RELEASE: 0 32 // AUTHOR: F. Hecht, 33 // ORG : UMPC 34 // E-MAIL : Frederic.Hecht@Inria.fr 35 // 36 // ORIG-DATE: frev 98 37 //--------------------------------------------------------- 38 // to make quad 39 // ------------------- 40 41 #include <cstdio> 42 #include <cstring> 43 #include <cmath> 44 #include <ctime> 45 46 #include "Meshio.h" 47 #include "Mesh2.h" 48 #include "QuadTree.h" 49 #include "SetOfE4.h" 50 51 namespace bamg { 52 53 static const Direction NoDirOfSearch = Direction( ); 54 55 #ifdef __MWERKS__ 56 #pragma global_optimizer on 57 #pragma optimization_level 1 58 //#pragma inline_depth 0 59 #endif 60 61 class DoubleAndInt4 { 62 public: 63 double q; 64 Int4 i3j; operator <(DoubleAndInt4 a)65 int operator<(DoubleAndInt4 a) { return q > a.q; } 66 }; // to sort by decroissant 67 68 template< class T > HeapSort(T * c,long n)69 inline void HeapSort(T *c, long n) { 70 long l, j, r, i; 71 T crit; 72 c--; // on decale de 1 pour que le tableau commence a 1 73 if (n <= 1) return; 74 l = n / 2 + 1; 75 r = n; 76 while (1) { // label 2 77 if (l <= 1) { // label 20 78 crit = c[r]; 79 c[r--] = c[1]; 80 if (r == 1) { 81 c[1] = crit; 82 return; 83 } 84 } else 85 crit = c[--l]; 86 j = l; 87 while (1) { // label 4 88 i = j; 89 j = 2 * j; 90 if (j > r) { 91 c[i] = crit; 92 break; 93 } // L8 -> G2 94 if ((j < r) && (c[j] < c[j + 1])) j++; // L5 95 if (crit < c[j]) 96 c[i] = c[j]; // L6+1 G4 97 else { 98 c[i] = crit; 99 break; 100 } // L8 -> G2 101 } 102 } 103 } 104 105 Triangle *swapTest(Triangle *t1, Int2 a); 106 // 107 QuadQuality(const Vertex & a,const Vertex & b,const Vertex & c,const Vertex & d)108 double QuadQuality(const Vertex &a, const Vertex &b, const Vertex &c, const Vertex &d) { 109 // calcul de 4 angles -- 110 R2 A((R2)a), B((R2)b), C((R2)c), D((R2)d); 111 R2 AB(B - A), BC(C - B), CD(D - C), DA(A - D); 112 // Move(A),Line(B),Line(C),Line(D),Line(A); 113 const Metric &Ma = a; 114 const Metric &Mb = b; 115 const Metric &Mc = c; 116 const Metric &Md = d; 117 118 double lAB = Norme2(AB); 119 double lBC = Norme2(BC); 120 double lCD = Norme2(CD); 121 double lDA = Norme2(DA); 122 AB /= lAB; 123 BC /= lBC; 124 CD /= lCD; 125 DA /= lDA; 126 // version aniso 127 double cosDAB = Ma(DA, AB) / (Ma(DA) * Ma(AB)), sinDAB = Det(DA, AB); 128 double cosABC = Mb(AB, BC) / (Mb(AB) * Mb(BC)), sinABC = Det(AB, BC); 129 double cosBCD = Mc(BC, CD) / (Mc(BC) * Mc(CD)), sinBCD = Det(BC, CD); 130 double cosCDA = Md(CD, DA) / (Md(CD) * Md(DA)), sinCDA = Det(CD, DA); 131 double sinmin = Min(Min(sinDAB, sinABC), Min(sinBCD, sinCDA)); 132 // cout << A << B << C << D ; 133 // cout << " sinmin " << sinmin << " " << cosDAB << " " << cosABC << " " << cosBCD << " " << 134 // cosCDA << endl; rattente(1); 135 if (sinmin <= 0) return sinmin; 136 return 1.0 - Max(Max(Abs(cosDAB), Abs(cosABC)), Max(Abs(cosBCD), Abs(cosCDA))); 137 } 138 ProjectOnCurve(Edge & BhAB,Vertex & vA,Vertex & vB,Real8 theta,Vertex & R,VertexOnEdge & BR,VertexOnGeom & GR)139 GeometricalEdge *Triangles::ProjectOnCurve(Edge &BhAB, Vertex &vA, Vertex &vB, Real8 theta, 140 Vertex &R, VertexOnEdge &BR, VertexOnGeom &GR) 141 142 { 143 void *pA = 0, *pB = 0; 144 Real8 tA = 0, tB = 0; 145 R2 A = vA, B = vB; 146 Vertex *pvA = &vA, *pvB = &vB; 147 if (vA.vint == IsVertexOnVertex) { 148 // cout << " Debut vertex = " << BTh.Number(vA.onbv) ; 149 pA = vA.onbv; 150 } else if (vA.vint == IsVertexOnEdge) { 151 pA = vA.onbe->be; 152 tA = vA.onbe->abcisse; 153 // cout << " Debut edge = " << BTh.Number(vA.onbv) << " " << tA ; 154 155 } else { 156 cerr << "ProjectOnCurve On Vertex " << BTh.Number(vA) << " " << endl; 157 cerr << " forget call to SetVertexFieldOnBTh" << endl; 158 MeshError(-1); 159 } 160 161 if (vB.vint == IsVertexOnVertex) { 162 // cout << " Fin vertex = " << BTh.Number(vB.onbv) << endl; 163 pB = vB.onbv; 164 } else if (vB.vint == IsVertexOnEdge) { 165 pB = vB.onbe->be; 166 tB = vB.onbe->abcisse; 167 // cout << " Fin edge = " << BTh.Number(vB.onbe->be) << " " << tB ; 168 169 } else { 170 cerr << "ProjectOnCurve On Vertex " << BTh.Number(vB) << " " << endl; 171 cerr << " forget call to SetVertexFieldOnBTh" << endl; 172 MeshError(-1); 173 } 174 Edge *e = &BhAB; 175 assert(pA && pB && e); 176 // be carefull the back ground edge e is on same geom edge 177 // of the initiale edge def by the 2 vertex A B; 178 assert(e >= BTh.edges && e < BTh.edges + BTh.nbe); // Is a background Mesh; 179 // walk on BTh edge 180 // assert(0 /* not finish ProjectOnCurve with BackGround Mesh*/); 181 // 1 first find a back ground edge contening the vertex A 182 // 2 walk n back gound boundary to find the final vertex B 183 184 #ifdef DEBUG 185 // we suppose if the vertex A is on a background edge and 186 // the abcisse is 0 or 1 when the related point is not 187 // a end of curve <=> !IsRequiredVertex 188 if (vA.vint == IsVertexOnEdge) 189 if (tA <= 0) 190 assert(!(*vA.onbe->be)[0].on->IsRequiredVertex( )); 191 else if (tA >= 1) 192 assert(!(*vA.onbe->be)[1].on->IsRequiredVertex( )); 193 #endif 194 195 if (vA.vint == IsVertexOnEdge) { // find the start edge 196 e = vA.onbe->be; 197 198 } else if (vB.vint == IsVertexOnEdge) { 199 theta = 1 - theta; 200 Exchange(tA, tB); 201 Exchange(pA, pB); 202 Exchange(pvA, pvB); 203 Exchange(A, B); 204 e = vB.onbe->be; 205 206 // cout << " EXCHANGE A et B) " << endl; 207 } else { // do the search by walking 208 assert(0 /* A FAIRE */); 209 } 210 211 // find the direction of walking with sens of edge and pA,PB; 212 R2 AB = B - A; 213 214 Real8 cosE01AB = (((R2)(*e)[1] - (R2)(*e)[0]), AB); 215 int kkk = 0; 216 int sens = (cosE01AB > 0) ? 1 : 0; 217 218 // Real8 l=0; // length of the edge AB 219 Real8 abscisse = -1; 220 221 for (int cas = 0; cas < 2; cas++) { // 2 times algo: 222 // 1 for computing the length l 223 // 2 for find the vertex 224 int iii; 225 Vertex *v0 = pvA, *v1; 226 Edge *neee, *eee; 227 Real8 lg = 0; // length of the curve 228 Real8 te0; 229 // we suppose take the curve's abcisse 230 // cout << kkk << " e = " << BTh.Number(e) << " v0= " 231 // << BTh.Number(v0) << " v1 = " << BTh.Number((*e)[sens]) << endl; 232 for (eee = e, iii = sens, te0 = tA; 233 eee && (((void *)eee) != pB) && ((void *)(v1 = &((*eee)[iii]))) != pB; 234 neee = eee->adj[iii], iii = 1 - neee->Intersection(*eee), eee = neee, v0 = v1, 235 te0 = 1 - iii) { 236 // cout << kkk << " eee = " << BTh.Number(eee) << " v0= " 237 // << BTh.Number(v0) << " v1 = " << BTh.Number(v1) << endl; 238 239 assert(kkk++ < 100); 240 assert(eee); 241 Real8 lg0 = lg; 242 Real8 dp = LengthInterpole(v0->m, v1->m, (R2)*v1 - (R2)*v0); 243 lg += dp; 244 if (cas && abscisse <= lg) { // ok we find the geom edge 245 Real8 sss = (abscisse - lg0) / dp; 246 Real8 thetab = te0 * (1 - sss) + sss * iii; 247 assert(thetab >= 0 && thetab <= 1); 248 BR = VertexOnEdge(&R, eee, thetab); 249 250 // cout << Number(R) << " = " << thetab << " on " << BTh.Number(eee) 251 // << " = " << R << endl; 252 253 return Gh.ProjectOnCurve(*eee, thetab, R, GR); 254 } 255 } 256 // we find the end 257 if (v1 != pvB) { 258 if ((void *)v1 == pB) tB = iii; 259 260 Real8 lg0 = lg; 261 assert(eee); 262 v1 = pvB; 263 Real8 dp = LengthInterpole(v0->m, v1->m, (R2)*v1 - (R2)*v0); 264 lg += dp; 265 abscisse = lg * theta; 266 if (abscisse <= lg && 267 abscisse >= lg0) // small optimisation we know the lenght because end 268 { // ok we find the geom edge 269 Real8 sss = (abscisse - lg0) / dp; 270 Real8 thetab = te0 * (1 - sss) + sss * tB; 271 assert(thetab >= 0 && thetab <= 1); 272 BR = VertexOnEdge(&R, eee, thetab); 273 274 // cout << kkk << " eee = " << BTh.Number(eee) << " v0= " 275 // << BTh.Number(v0) << " " << te0 276 // << " v1 = " << BTh.Number(v1) << " " << tB << endl; 277 278 // out << Number(R) << " Opt = " << thetab << " on " << BTh.Number(eee) 279 // << " = " << R << endl; 280 281 return Gh.ProjectOnCurve(*eee, thetab, R, GR); 282 } 283 } 284 abscisse = lg * theta; 285 } 286 cerr << " Big Bug" << endl; 287 MeshError(678); 288 return 0; // just for the compiler 289 } 290 MakeQuadrangles(double costheta)291 void Triangles::MakeQuadrangles(double costheta) { 292 if (verbosity > 2) cout << " -- MakeQuadrangles costheta = " << costheta << endl; 293 if (verbosity > 5) 294 cout << " (in) Nb of Quadrilaterals = " << NbOfQuad 295 << " Nb Of Triangles = " << nbt - NbOutT - NbOfQuad * 2 296 << " Nb of outside triangles = " << NbOutT << endl; 297 298 if (costheta > 1) { 299 if (verbosity > 5) cout << " do nothing costheta >1 " << endl; 300 return; 301 } 302 303 Int4 nbqq = (nbt * 3) / 2; 304 DoubleAndInt4 *qq = new DoubleAndInt4[nbqq]; 305 306 Int4 i, ij; 307 int j; 308 Int4 k = 0; 309 for (i = 0; i < nbt; i++) 310 for (j = 0; j < 3; j++) 311 if ((qq[k].q = triangles[i].QualityQuad(j)) >= costheta) qq[k++].i3j = i * 3 + j; 312 // sort qq 313 HeapSort(qq, k); 314 315 Int4 kk = 0; 316 for (ij = 0; ij < k; ij++) { 317 // cout << qq[ij].q << " " << endl; 318 i = qq[ij].i3j / 3; 319 j = (int)(qq[ij].i3j % 3); 320 // optisamition no float computation 321 if (triangles[i].QualityQuad(j, 0) >= costheta) triangles[i].SetHidden(j), kk++; 322 } 323 NbOfQuad = kk; 324 if (verbosity > 2) { 325 cout << " (out) Nb of Quadrilaterals = " << NbOfQuad 326 << " Nb Of Triangles = " << nbt - NbOutT - NbOfQuad * 2 327 << " Nb of outside triangles = " << NbOutT << endl; 328 } 329 delete[] qq; 330 #ifdef DRAWING2 331 Draw( ); 332 inquire( ); 333 #endif 334 } 335 /* 336 Triangles::BThBoundary(Edge e,Real 8) const 337 { 338 // pointeur of the background must be on 339 // 340 Edge be = e.on; 341 } 342 */ SplitElement(int choice)343 int Triangles::SplitElement(int choice) { 344 Direction NoDirOfSearch; 345 const int withBackground = &BTh != this; 346 if (verbosity > 2) 347 cout << " -- SplitElement " << (choice ? " Q->4Q and T->4T " : " Q->4Q or T->3Q ") << endl; 348 ; 349 if (verbosity > 5) 350 cout << endl 351 << " (in) Nb of Quadrilaterals = " << NbOfQuad 352 << " Nb Of Triangles = " << nbt - NbOutT - NbOfQuad * 2 353 << " Nb of outside triangles = " << NbOutT << endl; 354 355 ReNumberingTheTriangleBySubDomain( ); 356 #ifdef DRAWING2 357 Draw( ); 358 inquire( ); 359 #endif 360 // int nswap =0; 361 const Int4 nfortria(choice ? 4 : 6); 362 if (withBackground) { 363 BTh.SetVertexFieldOn( ); 364 SetVertexFieldOnBTh( ); 365 } else 366 BTh.SetVertexFieldOn( ); 367 368 Int4 newnbt = 0, newnbv = 0; 369 Int4 *kedge = 0; 370 Int4 newNbOfQuad = NbOfQuad; 371 Int4 *ksplit = 0, *ksplitarray = 0; 372 Int4 kkk = 0; 373 int ret = 0; 374 if (nbvx < nbv + nbe) return 1; // 375 Triangles *OCurrentTh = CurrentTh; 376 CurrentTh = this; 377 // 1) create the new points by spliting the internal edges 378 // set the 379 Int4 nbvold = nbv; 380 Int4 nbtold = nbt; 381 Int4 NbOutTold = NbOutT; 382 Int4 NbEdgeOnGeom = 0; 383 Int4 i; 384 385 nbt = nbt - NbOutT; // remove all the the ouside triangles 386 Int4 nbtsave = nbt; 387 Triangle *lastT = triangles + nbt; 388 for (i = 0; i < nbe; i++) 389 if (edges[i].on) NbEdgeOnGeom++; 390 Int4 newnbe = nbe + nbe; 391 // Int4 newNbVerticesOnGeomVertex=NbVerticesOnGeomVertex; 392 Int4 newNbVerticesOnGeomEdge = NbVerticesOnGeomEdge + NbEdgeOnGeom; 393 // Int4 newNbVertexOnBThVertex=NbVertexOnBThVertex; 394 Int4 newNbVertexOnBThEdge = withBackground ? NbVertexOnBThEdge + NbEdgeOnGeom : 0; 395 396 // do allocation for pointeur to the geometry and background 397 VertexOnGeom *newVerticesOnGeomEdge = new VertexOnGeom[newNbVerticesOnGeomEdge]; 398 VertexOnEdge *newVertexOnBThEdge = 399 newNbVertexOnBThEdge ? new VertexOnEdge[newNbVertexOnBThEdge] : 0; 400 if (NbVerticesOnGeomEdge) 401 memcpy(newVerticesOnGeomEdge, VerticesOnGeomEdge, 402 sizeof(VertexOnGeom) * NbVerticesOnGeomEdge); 403 if (NbVertexOnBThEdge) 404 memcpy(newVertexOnBThEdge, VertexOnBThEdge, sizeof(VertexOnEdge) * NbVertexOnBThEdge); 405 Edge *newedges = new Edge[newnbe]; 406 // memcpy(newedges,edges,sizeof(Edge)*nbe); 407 SetOfEdges4 *edge4 = new SetOfEdges4(nbe, nbv); 408 #ifdef DEBUG 409 for (i = 0; i < nbt; i++) triangles[i].check( ); 410 #endif 411 #ifdef DRAWING2 412 reffecran( ); 413 #endif 414 Int4 k = nbv; 415 Int4 kk = 0; 416 Int4 kvb = NbVertexOnBThEdge; 417 Int4 kvg = NbVerticesOnGeomEdge; 418 Int4 ie = 0; 419 Edge **edgesGtoB = 0; 420 if (withBackground) edgesGtoB = BTh.MakeGeometricalEdgeToEdge( ); 421 Int4 ferr = 0; 422 for (i = 0; i < nbe; i++) newedges[ie].on = 0; 423 424 for (i = 0; i < nbe; i++) { 425 GeometricalEdge *ong = edges[i].on; 426 427 newedges[ie] = edges[i]; 428 newedges[ie].adj[0] = newedges + (edges[i].adj[0] - edges); 429 newedges[ie].adj[1] = newedges + ie + 1; 430 R2 A = edges[i][0], B = edges[i][1]; 431 // cout << " ie = " << ie <<" v0 = " << Number(newedges[ie][0]) << endl; 432 433 kk += (i == edge4->addtrie(Number(edges[i][0]), Number(edges[i][1]))); 434 if (ong) // a geometrical edges 435 { 436 if (withBackground) { 437 // walk on back ground mesh 438 // newVertexOnBThEdge[ibe++] = VertexOnEdge(vertices[k],bedge,absicsseonBedge); 439 // a faire -- difficile 440 // the first PB is to now a background edge between the 2 vertices 441 assert(edgesGtoB); 442 // cout << " ie = " << ie <<" v0 = " << Number(newedges[ie][0]) << endl; 443 ong = ProjectOnCurve(*edgesGtoB[Gh.Number(edges[i].on)], edges[i][0], edges[i][1], 0.5, 444 vertices[k], newVertexOnBThEdge[kvb], newVerticesOnGeomEdge[kvg++]); 445 vertices[k].ReferenceNumber = edges[i].ref; 446 vertices[k].DirOfSearch = NoDirOfSearch; 447 ; 448 // get the Info on background mesh 449 Real8 s = newVertexOnBThEdge[kvb]; 450 Vertex &bv0 = newVertexOnBThEdge[kvb][0]; 451 Vertex &bv1 = newVertexOnBThEdge[kvb][1]; 452 // compute the metrix of the new points 453 vertices[k].m = Metric(1 - s, bv0, s, bv1); 454 kvb++; 455 // cout << " ie = " << ie <<" v0 = " << Number(newedges[ie][0]) << endl; 456 } else { 457 ong = Gh.ProjectOnCurve(edges[i], 0.5, vertices[k], newVerticesOnGeomEdge[kvg++]); 458 // vertices[k].i = toI2( vertices[k].r); 459 vertices[k].ReferenceNumber = edges[i].ref; 460 vertices[k].DirOfSearch = NoDirOfSearch; 461 vertices[k].m = Metric(0.5, edges[i][0], 0.5, edges[i][1]); 462 } 463 } else // straigth line edge --- 464 { 465 vertices[k].r = ((R2)edges[i][0] + (R2)edges[i][1]) * 0.5; 466 vertices[k].m = Metric(0.5, edges[i][0], 0.5, edges[i][1]); 467 vertices[k].on = 0; 468 } 469 // vertices[k].i = toI2( vertices[k].r); 470 R2 AB = vertices[k].r; 471 R2 AA = (A + AB) * 0.5; 472 R2 BB = (AB + B) * 0.5; 473 vertices[k].ReferenceNumber = edges[i].ref; 474 vertices[k].DirOfSearch = NoDirOfSearch; 475 476 newedges[ie].on = Gh.Contening(AA, ong); 477 newedges[ie++].v[1] = vertices + k; 478 479 newedges[ie] = edges[i]; 480 newedges[ie].adj[0] = newedges + ie - 1; 481 newedges[ie].adj[1] = newedges + (edges[i].adj[1] - edges); 482 newedges[ie].on = Gh.Contening(BB, ong); 483 newedges[ie++].v[0] = vertices + k; 484 // cout << " ie = " << ie-2 << " vm " << k << " v0 = " << Number(newedges[ie-2][0]) 485 // << " v1 = " << Number(newedges[ie-1][1]) 486 // << " ong =" << ong-Gh.edges 487 // << " on 0 =" << newedges[ie-2].on -Gh.edges << AA 488 // << " on 1 =" << newedges[ie-1].on -Gh.edges << BB 489 // << endl; 490 k++; 491 } 492 #ifdef DEBUG 493 assert(kvb == newNbVertexOnBThEdge); 494 // verif edge 495 { 496 Vertex *v0 = vertices, *v1 = vertices + k; 497 for (Int4 i = 0; i < ie; i++) { 498 assert(&newedges[i][0] >= v0 && &newedges[i][0] < v1); 499 assert(&newedges[i][1] >= v0 && &newedges[i][1] < v1); 500 } 501 } 502 #endif 503 if (edgesGtoB) delete[] edgesGtoB; 504 edgesGtoB = 0; 505 506 newnbv = k; 507 newNbVerticesOnGeomEdge = kvg; 508 if (newnbv > nbvx) goto Error; // bug 509 510 nbv = k; 511 512 kedge = new Int4[3 * nbt + 1]; 513 ksplitarray = new Int4[nbt + 1]; 514 ksplit = ksplitarray + 1; // because ksplit[-1] == ksplitarray[0] 515 516 for (i = 0; i < 3 * nbt; i++) kedge[i] = -1; 517 518 // 519 520 for (i = 0; i < nbt; i++) { 521 522 Triangle &t = triangles[i]; 523 assert(t.link); 524 for (int j = 0; j < 3; j++) { 525 const TriangleAdjacent ta = t.Adj(j); 526 const Triangle &tt = ta; 527 if (&tt >= lastT) t.SetAdj2(j, 0, 0); // unset adj 528 const Vertex &v0 = t[VerticesOfTriangularEdge[j][0]]; 529 const Vertex &v1 = t[VerticesOfTriangularEdge[j][1]]; 530 Int4 ke = edge4->findtrie(Number(v0), Number(v1)); 531 if (ke > 0) { 532 Int4 ii = Number(tt); 533 int jj = ta; 534 Int4 ks = ke + nbvold; 535 kedge[3 * i + j] = ks; 536 if (ii < nbt) // good triangle 537 kedge[3 * ii + jj] = ks; 538 Vertex &A = vertices[ks]; 539 Real8 aa = 0, bb = 0, cc = 0, dd = 0; 540 if ((dd = Area2(v0.r, v1.r, A.r)) >= 0) { // warning PB roundoff error 541 if (t.link && ((aa = Area2(A.r, t[1].r, t[2].r)) < 0.0 || 542 (bb = Area2(t[0].r, A.r, t[2].r)) < 0.0 || 543 (cc = Area2(t[0].r, t[1].r, A.r)) < 0.0)) 544 ferr++, cerr << " Error : " << ke + nbvold << " not in triangle " << i 545 << " In=" << !!t.link << " " << aa << " " << bb << " " << cc << " " << dd 546 << endl; 547 548 } 549 550 else { 551 if (tt.link && ((aa = Area2(A.r, tt[1].r, tt[2].r)) < 0 || 552 (bb = Area2(tt[0].r, A.r, tt[2].r)) < 0 || 553 (cc = Area2(tt[0].r, tt[1].r, A.r)) < 0)) 554 ferr++, cerr << " Warning : " << ke + nbvold << " not in triangle " << ii 555 << " In=" << !!tt.link << " " << aa << " " << bb << " " << cc << " " 556 << dd << endl; 557 } 558 } 559 } 560 } 561 if (ferr) { 562 cerr << " Number of triangles with P2 interpolation Probleme " << ferr << endl; 563 ; 564 MeshError(9); 565 } 566 567 for (i = 0; i < nbt; i++) { 568 ksplit[i] = 1; // no split by default 569 const Triangle &t = triangles[i]; 570 // cout << " Triangle " << i << " " << t << !!t.link << ":: " ; 571 int nbsplitedge = 0; 572 int nbinvisible = 0; 573 int invisibleedge = 0; 574 int kkk[3]; 575 for (int j = 0; j < 3; j++) { 576 if (t.Hidden(j)) invisibleedge = j, nbinvisible++; 577 578 const TriangleAdjacent ta = t.Adj(j); 579 const Triangle &tt = ta; 580 581 const Vertex &v0 = t[VerticesOfTriangularEdge[j][0]]; 582 const Vertex &v1 = t[VerticesOfTriangularEdge[j][1]]; 583 // cout << " ke = " << kedge[3*i+j] << " " << Number(v0) << " " << Number(v1) << "/ "; 584 if (kedge[3 * i + j] < 0) { 585 Int4 ke = edge4->findtrie(Number(v0), Number(v1)); 586 // cout << ":" << ke << "," << !!t.link << " " << &tt ; 587 if (ke < 0) // new 588 { 589 if (ta.t) // internal triangles all the boundary 590 { // new internal edges 591 Int4 ii = Number(tt); 592 int jj = ta; 593 594 kedge[3 * i + j] = k; // save the vertex number 595 kedge[3 * ii + jj] = k; 596 if (k < nbvx) { 597 vertices[k].r = ((R2)v0 + (R2)v1) / 2; 598 // vertices[k].i = toI2( vertices[k].r); 599 vertices[k].ReferenceNumber = 0; 600 vertices[k].DirOfSearch = NoDirOfSearch; 601 vertices[k].m = Metric(0.5, v0, 0.5, v1); 602 } 603 k++; 604 kkk[nbsplitedge++] = j; 605 } // tt 606 else 607 cerr << endl << " Bug " << i << " " << j << " t=" << t << endl; 608 609 } // ke<0 610 else { // ke >=0 611 kedge[3 * i + j] = nbvold + ke; 612 kkk[nbsplitedge++] = j; // previously splited 613 } 614 } else 615 kkk[nbsplitedge++] = j; // previously splited 616 } 617 assert(nbinvisible < 2); 618 // cout << " " << nbinvisible << " " << nbsplitedge << endl; 619 switch (nbsplitedge) { 620 case 0: 621 ksplit[i] = 10; 622 newnbt++; 623 break; // nosplit 624 case 1: 625 ksplit[i] = 20 + kkk[0]; 626 newnbt += 2; 627 break; // split in 2 628 case 2: 629 ksplit[i] = 30 + 3 - kkk[0] - kkk[1]; 630 newnbt += 3; 631 break; // split in 3 632 case 3: 633 if (nbinvisible) 634 ksplit[i] = 40 + invisibleedge, newnbt += 4; 635 else 636 ksplit[i] = 10 * nfortria, newnbt += nfortria; 637 break; 638 } 639 assert(ksplit[i] >= 40); 640 } 641 // now do the element split 642 newNbOfQuad = 4 * NbOfQuad; 643 nbv = k; 644 #ifdef DRAWING2 645 inquire( ); 646 #endif 647 // cout << " Nbv = " << nbv << endl; 648 kkk = nbt; 649 ksplit[-1] = nbt; 650 // look on old true triangles 651 652 for (i = 0; i < nbtsave; i++) { 653 // cout << "Triangle " << i << " " << ksplit[i] << ":" << triangles[i] 654 // << " ----------------------------------------------- " <<endl; 655 // Triangle * tc=0; 656 int nbmkadj = 0; 657 Int4 mkadj[100]; 658 mkadj[0] = i; 659 Int4 kk = ksplit[i] / 10; 660 int ke = (int)(ksplit[i] % 10); 661 assert(kk < 7 && kk > 0); 662 663 // def the numbering k (edge) i vertex 664 int k0 = ke; 665 int k1 = NextEdge[k0]; 666 int k2 = PreviousEdge[k0]; 667 int i0 = OppositeVertex[k0]; 668 int i1 = OppositeVertex[k1]; 669 int i2 = OppositeVertex[k2]; 670 671 Triangle &t0 = triangles[i]; 672 Vertex *v0 = t0(i0); 673 Vertex *v1 = t0(i1); 674 Vertex *v2 = t0(i2); 675 676 // cout << "nbmkadj " << nbmkadj << " it=" << i <<endl; 677 assert(nbmkadj < 10); 678 // -------------------------- 679 TriangleAdjacent ta0(t0.Adj(i0)), ta1(t0.Adj(i1)), ta2(t0.Adj(i2)); 680 // save the flag Hidden 681 int hid[] = {t0.Hidden(0), t0.Hidden(1), t0.Hidden(2)}; 682 // un set all adj -- save Hidden flag -- 683 t0.SetAdj2(0, 0, hid[0]); 684 t0.SetAdj2(1, 0, hid[1]); 685 t0.SetAdj2(2, 0, hid[2]); 686 // -- remake 687 switch (kk) { 688 case 1: 689 break; // nothing 690 case 2: // 691 { 692 Triangle &t1 = triangles[kkk++]; 693 t1 = t0; 694 assert(kedge[3 * i + i0] >= 0); 695 Vertex *v3 = vertices + kedge[3 * i + k0]; 696 697 t0(i2) = v3; 698 t1(i1) = v3; 699 t0.SetAllFlag(k2, 0); 700 t1.SetAllFlag(k1, 0); 701 } break; 702 case 3: // 703 { 704 Triangle &t1 = triangles[kkk++]; 705 Triangle &t2 = triangles[kkk++]; 706 t2 = t1 = t0; 707 assert(kedge[3 * i + k1] >= 0); 708 assert(kedge[3 * i + k2] >= 0); 709 710 Vertex *v01 = vertices + kedge[3 * i + k2]; 711 Vertex *v02 = vertices + kedge[3 * i + k1]; 712 t0(i1) = v01; 713 t0(i2) = v02; 714 t1(i2) = v02; 715 t1(i0) = v01; 716 t2(i0) = v02; 717 t0.SetAllFlag(k0, 0); 718 t1.SetAllFlag(k1, 0); 719 t1.SetAllFlag(k0, 0); 720 t2.SetAllFlag(k2, 0); 721 } break; 722 case 4: // 723 case 6: // split in 4 724 { 725 Triangle &t1 = triangles[kkk++]; 726 Triangle &t2 = triangles[kkk++]; 727 Triangle &t3 = triangles[kkk++]; 728 t3 = t2 = t1 = t0; 729 assert(kedge[3 * i + k0] >= 0 && kedge[3 * i + k1] >= 0 && kedge[3 * i + k2] >= 0); 730 Vertex *v12 = vertices + kedge[3 * i + k0]; 731 Vertex *v02 = vertices + kedge[3 * i + k1]; 732 Vertex *v01 = vertices + kedge[3 * i + k2]; 733 // cout << Number(t0(i0)) << " " << Number(t0(i1)) 734 // << " " << Number(t0(i2)) 735 // << " " << kedge[3*i+k0] 736 // << " " << kedge[3*i+k1] 737 // << " " << kedge[3*i+k2] << endl; 738 t0(i1) = v01; 739 t0(i2) = v02; 740 t0.SetAllFlag(k0, hid[k0]); 741 742 t1(i0) = v01; 743 t1(i2) = v12; 744 t0.SetAllFlag(k1, hid[k1]); 745 746 t2(i0) = v02; 747 t2(i1) = v12; 748 t2.SetAllFlag(k2, hid[k2]); 749 750 t3(i0) = v12; 751 t3(i1) = v02; 752 t3(i2) = v01; 753 754 t3.SetAllFlag(0, hid[0]); 755 t3.SetAllFlag(1, hid[1]); 756 t3.SetAllFlag(2, hid[2]); 757 758 if (kk == 6) { 759 760 Triangle &t4 = triangles[kkk++]; 761 Triangle &t5 = triangles[kkk++]; 762 763 t4 = t3; 764 t5 = t3; 765 766 t0.SetHidden(k0); 767 t1.SetHidden(k1); 768 t2.SetHidden(k2); 769 t3.SetHidden(0); 770 t4.SetHidden(1); 771 t5.SetHidden(2); 772 773 if (nbv < nbvx) { 774 vertices[nbv].r = ((R2)*v01 + (R2)*v12 + (R2)*v02) / 3.0; 775 vertices[nbv].ReferenceNumber = 0; 776 vertices[nbv].DirOfSearch = NoDirOfSearch; 777 // vertices[nbv].i = toI2(vertices[nbv].r); 778 Real8 a3[] = {1. / 3., 1. / 3., 1. / 3.}; 779 vertices[nbv].m = Metric(a3, v0->m, v1->m, v2->m); 780 Vertex *vc = vertices + nbv++; 781 t3(i0) = vc; 782 t4(i1) = vc; 783 t5(i2) = vc; 784 785 } else 786 goto Error; 787 } 788 789 } break; 790 } 791 792 // cout << " -- " << i << " " << nbmkadj << " " << kkk << " " << tc << endl; 793 // t0.SetDetf(); 794 // save all the new triangles 795 mkadj[nbmkadj++] = i; 796 Int4 jj; 797 if (t0.link) 798 for (jj = nbt; jj < kkk; jj++) { 799 triangles[jj].link = t0.link; 800 t0.link = triangles + jj; 801 mkadj[nbmkadj++] = jj; 802 // triangles[jj].SetDet(); 803 } 804 // cout << " -- " << i << " " << nbmkadj << endl; 805 assert(nbmkadj <= 13); // 13 = 6 + 4 + 3 806 807 if (kk == 6) newNbOfQuad += 3; 808 // triangles[i].Draw(); 809 810 for (jj = ksplit[i - 1]; jj < kkk; jj++) 811 // triangles[jj].SetDet(); 812 // triangles[jj].Draw(); 813 814 nbt = kkk; 815 ksplit[i] = nbt; // save last adresse of the new triangles 816 kkk = nbt; 817 } 818 819 // cout << " nv = " << nbv << " nbt = " << nbt << endl; 820 for (i = 0; i < nbv; i++) vertices[i].m = vertices[i].m * 2.; 821 // 822 if (withBackground) 823 for (i = 0; i < BTh.nbv; i++) BTh.vertices[i].m = BTh.vertices[i].m * 2.; 824 #ifdef DRAWING2 825 Draw( ); 826 inquire( ); 827 #endif 828 829 ret = 2; 830 if (nbt >= nbtx) goto Error; // bug 831 if (nbv >= nbvx) goto Error; // bug 832 // generation of the new triangles 833 834 SetIntCoor("In SplitElement"); 835 836 ReMakeTriangleContainingTheVertex( ); 837 if (withBackground) BTh.ReMakeTriangleContainingTheVertex( ); 838 839 delete[] edges; 840 edges = newedges; 841 nbe = newnbe; 842 NbOfQuad = newNbOfQuad; 843 844 for (i = 0; i < NbSubDomains; i++) { 845 Int4 k = subdomains[i].edge - edges; 846 subdomains[i].edge = edges + 2 * k; // spilt all edge in 2 847 } 848 849 if (ksplitarray) delete[] ksplitarray; 850 if (kedge) delete[] kedge; 851 if (edge4) delete edge4; 852 if (VerticesOnGeomEdge) delete[] VerticesOnGeomEdge; 853 VerticesOnGeomEdge = newVerticesOnGeomEdge; 854 if (VertexOnBThEdge) delete[] VertexOnBThEdge; 855 VertexOnBThEdge = newVertexOnBThEdge; 856 NbVerticesOnGeomEdge = newNbVerticesOnGeomEdge; 857 NbVertexOnBThEdge = newNbVertexOnBThEdge; 858 // ReMakeTriangleContainingTheVertex(); 859 860 FillHoleInMesh( ); 861 862 #ifdef DEBUG 863 for (i = 0; i < nbt; i++) triangles[i].check( ); 864 #endif 865 866 if (verbosity > 2) 867 cout << " (out) Nb of Quadrilaterals = " << NbOfQuad 868 << " Nb Of Triangles = " << nbt - NbOutT - NbOfQuad * 2 869 << " Nb of outside triangles = " << NbOutT << endl; 870 871 CurrentTh = OCurrentTh; 872 return 0; // ok 873 874 Error: 875 nbv = nbvold; 876 nbt = nbtold; 877 NbOutT = NbOutTold; 878 // cleaning memory --- 879 if (newedges) delete[] newedges; 880 if (ksplitarray) delete[] ksplitarray; 881 if (kedge) delete[] kedge; 882 if (newVerticesOnGeomEdge) delete[] newVerticesOnGeomEdge; 883 if (edge4) delete edge4; 884 if (newVertexOnBThEdge) delete[] newVertexOnBThEdge; 885 886 CurrentTh = OCurrentTh; 887 return ret; // ok 888 } 889 890 } // namespace bamg 891