1 #include <mystdlib.h> 2 #include "meshing.hpp" 3 4 #define noDEBUG 5 6 7 namespace netgen 8 { 9 class MarkedTet; 10 class MarkedPrism; 11 class MarkedIdentification; 12 class MarkedTri; 13 class MarkedQuad; 14 15 typedef NgArray<MarkedTet> T_MTETS; 16 typedef NgArray<MarkedPrism> T_MPRISMS; 17 typedef NgArray<MarkedIdentification> T_MIDS; 18 typedef NgArray<MarkedTri> T_MTRIS; 19 typedef NgArray<MarkedQuad> T_MQUADS; 20 21 22 23 class MarkedTet 24 { 25 public: 26 /// pnums of tet 27 PointIndex pnums[4]; 28 /// material number 29 int matindex; 30 /// element marked for refinement 31 /// marked = 1: marked by element marker, marked = 2 due to closure 32 unsigned int marked:2; 33 /// flag of Arnold-Mukherjee algorithm 34 unsigned int flagged:1; 35 /// tetedge (local coordinates 0..3) 36 unsigned int tetedge1:3; 37 unsigned int tetedge2:3; 38 // marked edge of faces 39 // face_j : face without node j, 40 // mark_k : edge without node k 41 42 char faceedges[4]; 43 // unsigned char faceedges[4]; 44 bool incorder; 45 unsigned int order:6; 46 47 MarkedTet() = default; 48 /* 49 { 50 for (int i = 0; i < 4; i++) { faceedges[i] = 127; } 51 } 52 */ 53 MarkedTet (const MarkedTet&) = default; 54 MarkedTet (MarkedTet &&) = default; 55 MarkedTet & operator= (const MarkedTet&) = default; 56 MarkedTet & operator= (MarkedTet&&) = default; 57 58 }; 59 operator <<(ostream & ost,const MarkedTet & mt)60 ostream & operator<< (ostream & ost, const MarkedTet & mt) 61 { 62 for(int i=0; i<4; i++) 63 ost << mt.pnums[i] << " "; 64 65 ost << mt.matindex << " " << int(mt.marked) << " " << int(mt.flagged) << " " << int(mt.tetedge1) << " " << int(mt.tetedge2) << " "; 66 67 ost << "faceedges = "; 68 for(int i=0; i<4; i++) 69 ost << int(mt.faceedges[i]) << " "; 70 71 ost << " order = "; 72 ost << mt.incorder << " " << int(mt.order) << "\n"; 73 return ost; 74 } operator >>(istream & ost,MarkedTet & mt)75 istream & operator>> (istream & ost, MarkedTet & mt) 76 { 77 for(int i=0; i<4; i++) 78 ost >> mt.pnums[i]; 79 80 ost >> mt.matindex; 81 82 int auxint; 83 ost >> auxint; 84 mt.marked = auxint; 85 ost >> auxint; 86 mt.flagged = auxint; 87 ost >> auxint; 88 mt.tetedge1 = auxint; 89 ost >> auxint; 90 mt.tetedge2 = auxint; 91 92 char auxchar; 93 94 for(int i=0; i<4; i++) 95 { 96 ost >> auxchar; 97 mt.faceedges[i] = auxchar; 98 } 99 100 ost >> mt.incorder; 101 ost >> auxint; 102 mt.order = auxint; 103 return ost; 104 } 105 106 class MarkedPrism 107 { 108 public: 109 /// 6 point numbers 110 PointIndex pnums[6]; 111 /// material number 112 int matindex; 113 /// marked for refinement 114 int marked; 115 /// edge without node k (0,1,2) 116 int markededge; 117 118 bool incorder; 119 unsigned int order:6; 120 }; 121 122 operator <<(ostream & ost,const MarkedPrism & mp)123 ostream & operator<< (ostream & ost, const MarkedPrism & mp) 124 { 125 for(int i=0; i<6; i++) 126 ost << mp.pnums[i] << " "; 127 128 ost << mp.matindex << " " << mp.marked << " " << mp.markededge << " " << mp.incorder << " " << int(mp.order) << "\n"; 129 return ost; 130 } operator >>(istream & ist,MarkedPrism & mp)131 istream & operator>> (istream & ist, MarkedPrism & mp) 132 { 133 for(int i=0; i<6; i++) 134 ist >> mp.pnums[i]; 135 136 ist >> mp.matindex >> mp.marked >> mp.markededge >> mp.incorder; 137 int auxint; 138 ist >> auxint; 139 mp.order = auxint; 140 return ist; 141 } 142 143 144 class MarkedIdentification 145 { 146 public: 147 // number of points of one face (3 or 4) - or edge (in 2d) 148 int np; 149 /// 6 or 8 point numbers - or 4 in 2d 150 PointIndex pnums[8]; 151 /// marked for refinement 152 int marked; 153 /// edge starting with node k (0,1,2, or 3) 154 int markededge; 155 156 bool incorder; 157 unsigned int order:6; 158 }; 159 160 operator <<(ostream & ost,const MarkedIdentification & mi)161 ostream & operator<< (ostream & ost, const MarkedIdentification & mi) 162 { 163 ost << mi.np << " "; 164 for(int i=0; i<2*mi.np; i++) 165 ost << mi.pnums[i] << " "; 166 ost << mi.marked << " " << mi.markededge << " " << mi.incorder << " " << int(mi.order) << "\n"; 167 return ost; 168 } operator >>(istream & ist,MarkedIdentification & mi)169 istream & operator>> (istream & ist, MarkedIdentification & mi) 170 { 171 ist >> mi.np; 172 for(int i=0; i<2*mi.np; i++) 173 ist >> mi.pnums[i]; 174 ist >> mi.marked >> mi.markededge >> mi.incorder; 175 int auxint; 176 ist >> auxint; 177 mi.order = auxint; 178 return ist; 179 } 180 181 182 183 184 185 class MarkedTri 186 { 187 public: 188 MarkedTri () = default; 189 MarkedTri (const MarkedTri&) = default; 190 MarkedTri (MarkedTri &&) = default; 191 MarkedTri & operator= (const MarkedTri&) = default; 192 MarkedTri & operator= (MarkedTri&&) = default; 193 194 /// three point numbers 195 PointIndex pnums[3]; 196 /// three geominfos 197 PointGeomInfo pgeominfo[3]; 198 /// marked for refinement 199 int marked; 200 /// edge without node k 201 int markededge; 202 /// surface id 203 int surfid; 204 205 bool incorder; 206 unsigned int order:6; 207 }; 208 operator <<(ostream & ost,const MarkedTri & mt)209 ostream & operator<< (ostream & ost, const MarkedTri & mt) 210 { 211 for(int i=0; i<3; i++) 212 ost << mt.pnums[i] << " "; 213 for(int i=0; i<3; i++) 214 ost << mt.pgeominfo[i] << " "; 215 ost << mt.marked << " " << mt.markededge << " " << mt.surfid << " " << mt.incorder << " " << int(mt.order) << "\n"; 216 return ost; 217 } operator >>(istream & ist,MarkedTri & mt)218 istream & operator>> (istream & ist, MarkedTri & mt) 219 { 220 for(int i=0; i<3; i++) 221 ist >> mt.pnums[i]; 222 for(int i=0; i<3; i++) 223 ist >> mt.pgeominfo[i]; 224 ist >> mt.marked >> mt.markededge >> mt.surfid >> mt.incorder; 225 int auxint; 226 ist >> auxint; 227 mt.order = auxint; 228 return ist; 229 } 230 231 232 233 class MarkedQuad 234 { 235 public: 236 /// point numbers 237 PointIndex pnums[4]; 238 /// 239 PointGeomInfo pgeominfo[4]; 240 /// marked for refinement 241 int marked; 242 /// marked edge: 0/2 = vertical, 1/3 = horizontal 243 int markededge; 244 /// surface id 245 int surfid; 246 247 bool incorder; 248 unsigned int order:6; 249 }; 250 operator <<(ostream & ost,const MarkedQuad & mt)251 ostream & operator<< (ostream & ost, const MarkedQuad & mt) 252 { 253 for(int i=0; i<4; i++) 254 ost << mt.pnums[i] << " "; 255 for(int i=0; i<4; i++) 256 ost << mt.pgeominfo[i] << " "; 257 ost << mt.marked << " " << mt.markededge << " " << mt.surfid << " " << mt.incorder << " " << int(mt.order) << "\n"; 258 return ost; 259 } operator >>(istream & ist,MarkedQuad & mt)260 istream & operator>> (istream & ist, MarkedQuad & mt) 261 { 262 for(int i=0; i<4; i++) 263 ist >> mt.pnums[i]; 264 for(int i=0; i<4; i++) 265 ist >> mt.pgeominfo[i]; 266 ist >> mt.marked >> mt.markededge >> mt.surfid >> mt.incorder; 267 int auxint; 268 ist >> auxint; 269 mt.order = auxint; 270 return ist; 271 } 272 273 274 275 PrettyPrint(ostream & ost,const MarkedTet & mt)276 void PrettyPrint(ostream & ost, const MarkedTet & mt) 277 { 278 int te1 = mt.tetedge1; 279 int te2 = mt.tetedge2; 280 int order = mt.order; 281 282 ost << "MT: " << mt.pnums[0] << " - " << mt.pnums[1] << " - " 283 << mt.pnums[2] << " - " << mt.pnums[3] << endl 284 << "marked edge: " << te1 << " - " << te2 285 << ", order = " << order << endl; 286 //for (int k = 0; k < 4; k++) 287 // ost << int(mt.faceedges[k]) << " "; 288 for (int k = 0; k < 4; k++) 289 { 290 ost << "face"; 291 for (int j=0; j<4; j++) 292 if(j != k) 293 ost << " " << mt.pnums[j]; 294 for(int i=0; i<3; i++) 295 for(int j=i+1; j<4; j++) 296 if(i != k && j != k && int(mt.faceedges[k]) == 6-k-i-j) 297 ost << " marked edge " << mt.pnums[i] << " " << mt.pnums[j] << endl; 298 } 299 ost << endl; 300 } 301 302 303 304 BTSortEdges(const Mesh & mesh,const NgArray<NgArray<int,PointIndex::BASE> * > & idmaps,INDEX_2_CLOSED_HASHTABLE<int> & edgenumber)305 int BTSortEdges (const Mesh & mesh, 306 const NgArray< NgArray<int,PointIndex::BASE>* > & idmaps, 307 INDEX_2_CLOSED_HASHTABLE<int> & edgenumber) 308 { 309 PrintMessage(4,"sorting ... "); 310 311 // if (mesh.PureTetMesh()) 312 if (1) 313 { 314 // new, fast version 315 316 NgArray<INDEX_2> edges; 317 NgArray<int> eclasses; 318 319 int i, j, k; 320 int cntedges = 0; 321 int go_on; 322 int ned(0); 323 324 // enumerate edges: 325 for (i = 1; i <= mesh.GetNE(); i++) 326 { 327 const Element & el = mesh.VolumeElement (i); 328 static int tetedges[6][2] = 329 { { 1, 2 }, 330 { 1, 3 }, 331 { 1, 4 }, 332 { 2, 3 }, 333 { 2, 4 }, 334 { 3, 4 } } ; 335 static int prismedges[9][2] = 336 { { 1, 2 }, 337 { 1, 3 }, 338 { 2, 3 }, 339 { 4, 5 }, 340 { 4, 6 }, 341 { 5, 6 }, 342 { 1, 4 }, 343 { 2, 5 }, 344 { 3, 6 } }; 345 int pyramidedges[6][2] = 346 { { 1, 2 }, 347 { 3, 4 }, 348 { 1, 5 }, 349 { 2, 5 }, 350 { 3, 5 }, 351 { 4, 5 } }; 352 353 int (*tip)[2] = NULL; 354 355 switch (el.GetType()) 356 { 357 case TET: 358 case TET10: 359 { 360 tip = tetedges; 361 ned = 6; 362 break; 363 } 364 case PRISM: 365 case PRISM12: 366 { 367 tip = prismedges; 368 ned = 6; 369 break; 370 } 371 case PYRAMID: 372 { 373 tip = pyramidedges; 374 ned = 6; 375 break; 376 } 377 default: 378 throw NgException("Bisect, element type not handled in switch"); 379 } 380 381 for (j = 0; j < ned; j++) 382 { 383 INDEX_2 i2(el.PNum(tip[j][0]), el.PNum(tip[j][1])); 384 i2.Sort(); 385 //(*testout) << "edge " << i2 << endl; 386 if (!edgenumber.Used(i2)) 387 { 388 cntedges++; 389 edges.Append (i2); 390 edgenumber.Set(i2, cntedges); 391 } 392 } 393 } 394 395 // additional surface edges: 396 for (i = 1; i <= mesh.GetNSE(); i++) 397 { 398 const Element2d & el = mesh.SurfaceElement (i); 399 static int trigedges[3][2] = 400 { { 1, 2 }, 401 { 2, 3 }, 402 { 3, 1 } }; 403 404 static int quadedges[4][2] = 405 { { 1, 2 }, 406 { 2, 3 }, 407 { 3, 4 }, 408 { 4, 1 } }; 409 410 411 int (*tip)[2] = NULL; 412 413 switch (el.GetType()) 414 { 415 case TRIG: 416 case TRIG6: 417 { 418 tip = trigedges; 419 ned = 3; 420 break; 421 } 422 case QUAD: 423 case QUAD6: 424 { 425 tip = quadedges; 426 ned = 4; 427 break; 428 } 429 default: 430 { 431 cerr << "Error: Sort for Bisect, SE has " << el.GetNP() << " points" << endl; 432 ned = 0; 433 } 434 } 435 436 for (j = 0; j < ned; j++) 437 { 438 INDEX_2 i2(el.PNum(tip[j][0]), el.PNum(tip[j][1])); 439 i2.Sort(); 440 if (!edgenumber.Used(i2)) 441 { 442 cntedges++; 443 edges.Append (i2); 444 edgenumber.Set(i2, cntedges); 445 } 446 } 447 } 448 449 450 451 452 453 eclasses.SetSize (cntedges); 454 for (i = 1; i <= cntedges; i++) 455 eclasses.Elem(i) = i; 456 457 // identify edges in element stack 458 do 459 { 460 go_on = 0; 461 for (i = 1; i <= mesh.GetNE(); i++) 462 { 463 const Element & el = mesh.VolumeElement (i); 464 465 if (el.GetType() != PRISM && 466 el.GetType() != PRISM12 && 467 el.GetType() != PYRAMID) 468 continue; 469 470 int prismpairs[3][4] = 471 { { 1, 2, 4, 5 }, 472 { 2, 3, 5, 6 }, 473 { 1, 3, 4, 6 } }; 474 475 int pyramidpairs[3][4] = 476 { { 1, 2, 4, 3 }, 477 { 1, 5, 4, 5 }, 478 { 2, 5, 3, 5 } }; 479 480 int (*pairs)[4] = NULL; 481 switch (el.GetType()) 482 { 483 case PRISM: 484 case PRISM12: 485 { 486 pairs = prismpairs; 487 break; 488 } 489 case PYRAMID: 490 { 491 pairs = pyramidpairs; 492 break; 493 } 494 default: 495 throw NgException("Bisect, element type not handled in switch, 2"); 496 } 497 498 for (j = 0; j < 3; j++) 499 { 500 INDEX_2 e1 (el.PNum(pairs[j][0]), 501 el.PNum(pairs[j][1])); 502 INDEX_2 e2 (el.PNum(pairs[j][2]), 503 el.PNum(pairs[j][3])); 504 e1.Sort(); 505 e2.Sort(); 506 507 int eclass1 = edgenumber.Get (e1); 508 int eclass2 = edgenumber.Get (e2); 509 510 // (*testout) << "identify edges " << eclass1 << "-" << eclass2 << endl; 511 512 if (eclasses.Get(eclass1) > 513 eclasses.Get(eclass2)) 514 { 515 eclasses.Elem(eclass1) = 516 eclasses.Get(eclass2); 517 go_on = 1; 518 } 519 else if (eclasses.Get(eclass2) > 520 eclasses.Get(eclass1)) 521 { 522 eclasses.Elem(eclass2) = 523 eclasses.Get(eclass1); 524 go_on = 1; 525 } 526 } 527 } 528 529 for(SurfaceElementIndex sei = 0; sei < mesh.GetNSE(); sei++) 530 { 531 const Element2d & el2d = mesh[sei]; 532 533 for(i = 0; i < el2d.GetNP(); i++) 534 { 535 INDEX_2 e1(el2d[i], el2d[(i+1) % el2d.GetNP()]); 536 e1.Sort(); 537 INDEX_2 e2; 538 539 for(k = 0; k < idmaps.Size(); k++) 540 { 541 e2.I1() = (*idmaps[k])[e1.I1()]; 542 e2.I2() = (*idmaps[k])[e1.I2()]; 543 544 if(e2.I1() == 0 || e2.I2() == 0 || 545 e1.I1() == e2.I1() || e1.I2() == e2.I2()) 546 continue; 547 548 e2.Sort(); 549 if(!edgenumber.Used(e2)) 550 continue; 551 552 553 int eclass1 = edgenumber.Get (e1); 554 int eclass2 = edgenumber.Get (e2); 555 556 if (eclasses.Get(eclass1) > 557 eclasses.Get(eclass2)) 558 { 559 eclasses.Elem(eclass1) = 560 eclasses.Get(eclass2); 561 562 563 go_on = 1; 564 } 565 else if (eclasses.Get(eclass2) > 566 eclasses.Get(eclass1)) 567 { 568 eclasses.Elem(eclass2) = 569 eclasses.Get(eclass1); 570 go_on = 1; 571 } 572 } 573 } 574 575 } 576 577 } 578 while (go_on); 579 580 // for (i = 1; i <= cntedges; i++) 581 // { 582 // (*testout) << "edge " << i << ": " 583 // << edges.Get(i).I1() << "-" << edges.Get(i).I2() 584 // << ", class = " << eclasses.Get(i) << endl; 585 // } 586 587 // compute classlength: 588 NgArray<double> edgelength(cntedges); 589 590 /* 591 for (i = 1; i <= cntedges; i++) 592 edgelength.Elem(i) = 1e20; 593 */ 594 595 for (i = 1; i <= cntedges; i++) 596 { 597 INDEX_2 edge = edges.Get(i); 598 double elen = Dist (mesh.Point(edge.I1()), 599 mesh.Point(edge.I2())); 600 edgelength.Elem (i) = elen; 601 } 602 603 /* 604 for (i = 1; i <= mesh.GetNE(); i++) 605 { 606 const Element & el = mesh.VolumeElement (i); 607 608 if (el.GetType() == TET) 609 { 610 for (j = 1; j <= 3; j++) 611 for (k = j+1; k <= 4; k++) 612 { 613 INDEX_2 i2(el.PNum(j), el.PNum(k)); 614 i2.Sort(); 615 616 int enr = edgenumber.Get(i2); 617 double elen = Dist (mesh.Point (i2.I1()), mesh.Point (i2.I2())); 618 if (elen < edgelength.Get(enr)) 619 edgelength.Set (enr, elen); 620 } 621 } 622 else if (el.GetType() == PRISM) 623 { 624 for (j = 1; j <= 3; j++) 625 { 626 k = (j % 3) + 1; 627 628 INDEX_2 i2(el.PNum(j), el.PNum(k)); 629 i2.Sort(); 630 631 int enr = edgenumber.Get(i2); 632 double elen = Dist (mesh.Point (i2.I1()), mesh.Point (i2.I2())); 633 if (elen < edgelength.Get(enr)) 634 edgelength.Set (enr, elen); 635 636 i2 = INDEX_2(el.PNum(j+3), el.PNum(k+3)); 637 i2.Sort(); 638 639 enr = edgenumber.Get(i2); 640 elen = Dist (mesh.Point (i2.I1()), mesh.Point (i2.I2())); 641 if (elen < edgelength.Get(enr)) 642 edgelength.Set (enr, elen); 643 644 if (!edgenumber.Used(i2)) 645 { 646 cntedges++; 647 edgenumber.Set(i2, cntedges); 648 } 649 i2 = INDEX_2(el.PNum(j), el.PNum(j+3)); 650 i2.Sort(); 651 652 enr = edgenumber.Get(i2); 653 elen = Dist (mesh.Point (i2.I1()), mesh.Point (i2.I2())); 654 if (elen < edgelength.Get(enr)) 655 edgelength.Set (enr, elen); 656 } 657 } 658 } 659 */ 660 661 662 for (i = 1; i <= cntedges; i++) 663 { 664 if (eclasses.Get(i) != i) 665 { 666 if (edgelength.Get(i) < edgelength.Get(eclasses.Get(i))) 667 edgelength.Elem(eclasses.Get(i)) = edgelength.Get(i); 668 edgelength.Elem(i) = 1e20; 669 } 670 } 671 672 673 TABLE<int> eclasstab(cntedges); 674 for (i = 1; i <= cntedges; i++) 675 eclasstab.Add1 (eclasses.Get(i), i); 676 677 678 // sort edges: 679 NgArray<int> sorted(cntedges); 680 681 QuickSort (edgelength, sorted); 682 683 int cnt = 0; 684 for (i = 1; i <= cntedges; i++) 685 { 686 int ii = sorted.Get(i); 687 for (j = 1; j <= eclasstab.EntrySize(ii); j++) 688 { 689 cnt++; 690 edgenumber.Set (edges.Get(eclasstab.Get(ii, j)), cnt); 691 } 692 } 693 return cnt; 694 } 695 696 else 697 698 { 699 // old version 700 701 int i, j; 702 int cnt = 0; 703 int found; 704 double len2, maxlen2; 705 INDEX_2 ep; 706 707 // sort edges by length, parallel edges (on prisms) 708 // are added in blocks 709 710 do 711 { 712 found = 0; 713 maxlen2 = 1e30; 714 715 for (i = 1; i <= mesh.GetNE(); i++) 716 { 717 const Element & el = mesh.VolumeElement (i); 718 int ned; 719 int tetedges[6][2] = 720 { { 1, 2 }, 721 { 1, 3 }, 722 { 1, 4 }, 723 { 2, 3 }, 724 { 2, 4 }, 725 { 3, 4 } }; 726 int prismedges[6][2] = 727 { { 1, 2 }, 728 { 1, 3 }, 729 { 2, 4 }, 730 { 4, 5 }, 731 { 4, 6 }, 732 { 5, 6 } }; 733 int pyramidedges[6][2] = 734 { { 1, 2 }, 735 { 3, 4 }, 736 { 1, 5 }, 737 { 2, 5 }, 738 { 3, 5 }, 739 { 4, 5 } }; 740 741 int (*tip)[2]; 742 743 switch (el.GetType()) 744 { 745 case TET: 746 { 747 tip = tetedges; 748 ned = 6; 749 break; 750 } 751 case PRISM: 752 { 753 tip = prismedges; 754 ned = 6; 755 break; 756 } 757 case PYRAMID: 758 { 759 tip = pyramidedges; 760 ned = 6; 761 break; 762 } 763 default: 764 throw NgException("Bisect, element type not handled in switch, 3"); 765 } 766 767 for (j = 0; j < ned; j++) 768 { 769 INDEX_2 i2(el.PNum(tip[j][0]), el.PNum(tip[j][1])); 770 i2.Sort(); 771 if (!edgenumber.Used(i2)) 772 { 773 len2 = Dist (mesh.Point (i2.I1()), 774 mesh.Point (i2.I2())); 775 if (len2 < maxlen2) 776 { 777 maxlen2 = len2; 778 ep = i2; 779 found = 1; 780 } 781 } 782 } 783 } 784 if (found) 785 { 786 cnt++; 787 edgenumber.Set (ep, cnt); 788 789 790 // find connected edges: 791 int go_on = 0; 792 do 793 { 794 go_on = 0; 795 for (i = 1; i <= mesh.GetNE(); i++) 796 { 797 const Element & el = mesh.VolumeElement (i); 798 if (el.GetNP() != 6) continue; 799 800 int prismpairs[3][4] = 801 { { 1, 2, 4, 5 }, 802 { 2, 3, 5, 6 }, 803 { 1, 3, 4, 6 } }; 804 805 int pyramidpairs[3][4] = 806 { { 1, 2, 4, 3 }, 807 { 1, 5, 4, 5 }, 808 { 2, 5, 3, 5 } }; 809 810 int (*pairs)[4]; 811 switch (el.GetType()) 812 { 813 case PRISM: 814 { 815 pairs = prismpairs; 816 break; 817 } 818 case PYRAMID: 819 { 820 pairs = pyramidpairs; 821 break; 822 } 823 default: 824 throw NgException("Bisect, element type not handled in switch, 3a"); 825 } 826 827 for (j = 0; j < 3; j++) 828 { 829 INDEX_2 e1 (el.PNum(pairs[j][0]), 830 el.PNum(pairs[j][1])); 831 INDEX_2 e2 (el.PNum(pairs[j][2]), 832 el.PNum(pairs[j][3])); 833 e1.Sort(); 834 e2.Sort(); 835 836 int used1 = edgenumber.Used (e1); 837 int used2 = edgenumber.Used (e2); 838 839 if (used1 && !used2) 840 { 841 cnt++; 842 edgenumber.Set (e2, cnt); 843 go_on = 1; 844 } 845 if (used2 && !used1) 846 { 847 cnt++; 848 edgenumber.Set (e1, cnt); 849 go_on = 1; 850 } 851 } 852 } 853 } 854 while (go_on); 855 } 856 } 857 while (found); 858 859 return cnt; 860 } 861 } 862 863 864 865 BTDefineMarkedTet(const Element & el,INDEX_2_CLOSED_HASHTABLE<int> & edgenumber,MarkedTet & mt)866 void BTDefineMarkedTet (const Element & el, 867 INDEX_2_CLOSED_HASHTABLE<int> & edgenumber, 868 MarkedTet & mt) 869 { 870 for (int i = 0; i < 4; i++) 871 mt.pnums[i] = el[i]; 872 873 mt.marked = 0; 874 mt.flagged = 0; 875 876 mt.incorder = 0; 877 mt.order = 1; 878 879 int val = 0; 880 // find marked edge of tet: 881 for (int i = 0; i < 3; i++) 882 for (int j = i+1; j < 4; j++) 883 { 884 INDEX_2 i2(mt.pnums[i], mt.pnums[j]); 885 i2.Sort(); 886 int hval = edgenumber.Get(i2); 887 if (hval > val) 888 { 889 val = hval; 890 mt.tetedge1 = i; 891 mt.tetedge2 = j; 892 } 893 } 894 895 896 // find marked edges of faces: 897 for (int k = 0; k < 4; k++) 898 { 899 val = 0; 900 for (int i = 0; i < 3; i++) 901 for (int j = i+1; j < 4; j++) 902 if (i != k && j != k) 903 { 904 INDEX_2 i2(mt.pnums[i], mt.pnums[j]); 905 i2.Sort(); 906 int hval = edgenumber.Get(i2); 907 if (hval > val) 908 { 909 val = hval; 910 int hi = 6 - k - i - j; 911 mt.faceedges[k] = char(hi); 912 } 913 } 914 } 915 } 916 917 918 919 BTDefineMarkedPrism(const Element & el,INDEX_2_CLOSED_HASHTABLE<int> & edgenumber,MarkedPrism & mp)920 void BTDefineMarkedPrism (const Element & el, 921 INDEX_2_CLOSED_HASHTABLE<int> & edgenumber, 922 MarkedPrism & mp) 923 { 924 if (el.GetType() == PRISM || 925 el.GetType() == PRISM12) 926 { 927 for (int i = 0; i < 6; i++) 928 mp.pnums[i] = el[i]; 929 } 930 else if (el.GetType() == PYRAMID) 931 { 932 static int map[6] = 933 { 1, 2, 5, 4, 3, 5 }; 934 for (int i = 0; i < 6; i++) 935 mp.pnums[i] = el.PNum(map[i]); 936 } 937 else if (el.GetType() == TET || 938 el.GetType() == TET10) 939 { 940 static int map[6] = 941 { 1, 4, 3, 2, 4, 3 }; 942 for (int i = 0; i < 6; i++) 943 mp.pnums[i] = el.PNum(map[i]); 944 945 } 946 else 947 { 948 PrintSysError ("Define marked prism called for non-prism and non-pyramid"); 949 } 950 951 952 953 mp.marked = 0; 954 955 mp.incorder = 0; 956 mp.order = 1; 957 958 int val = 0; 959 for (int i = 0; i < 2; i++) 960 for (int j = i+1; j < 3; j++) 961 { 962 INDEX_2 i2(mp.pnums[i], mp.pnums[j]); 963 i2.Sort(); 964 int hval = edgenumber.Get(i2); 965 if (hval > val) 966 { 967 val = hval; 968 mp.markededge = 3 - i - j; 969 } 970 } 971 } 972 973 974 BTDefineMarkedId(const Element2d & el,INDEX_2_CLOSED_HASHTABLE<int> & edgenumber,const NgArray<int,PointIndex::BASE> & idmap,MarkedIdentification & mi)975 bool BTDefineMarkedId(const Element2d & el, 976 INDEX_2_CLOSED_HASHTABLE<int> & edgenumber, 977 const NgArray<int,PointIndex::BASE> & idmap, 978 MarkedIdentification & mi) 979 { 980 981 bool identified = true; 982 mi.np = el.GetNP(); 983 int min1(0),min2(0); 984 for(int j = 0; identified && j < mi.np; j++) 985 { 986 mi.pnums[j] = el[j]; 987 mi.pnums[j+mi.np] = idmap[el[j]]; 988 989 if(j == 0 || el[j] < min1) 990 min1 = el[j]; 991 if(j == 0 || mi.pnums[j+mi.np] < min2) 992 min2 = mi.pnums[j+mi.np]; 993 994 identified = (mi.pnums[j+mi.np] != 0 && mi.pnums[j+mi.np] != mi.pnums[j]); 995 } 996 997 identified = identified && (min1 < min2); 998 999 if(identified) 1000 { 1001 mi.marked = 0; 1002 1003 mi.incorder = 0; 1004 mi.order = 1; 1005 1006 int val = 0; 1007 for (int i = 0; i < mi.np; i++) 1008 { 1009 INDEX_2 i2(mi.pnums[i], mi.pnums[(i+1)%mi.np]); 1010 i2.Sort(); 1011 int hval = edgenumber.Get(i2); 1012 if (hval > val) 1013 { 1014 val = hval; 1015 mi.markededge = i; 1016 } 1017 } 1018 } 1019 1020 return identified; 1021 } 1022 1023 BTDefineMarkedTri(const Element2d & el,INDEX_2_CLOSED_HASHTABLE<int> & edgenumber,MarkedTri & mt)1024 void BTDefineMarkedTri (const Element2d & el, 1025 INDEX_2_CLOSED_HASHTABLE<int> & edgenumber, 1026 MarkedTri & mt) 1027 { 1028 for (int i = 0; i < 3; i++) 1029 { 1030 mt.pnums[i] = el[i]; 1031 mt.pgeominfo[i] = el.GeomInfoPi (i+1); 1032 } 1033 1034 mt.marked = 0; 1035 mt.surfid = el.GetIndex(); 1036 1037 mt.incorder = 0; 1038 mt.order = 1; 1039 1040 int val = 0; 1041 for (int i = 0; i < 2; i++) 1042 for (int j = i+1; j < 3; j++) 1043 { 1044 INDEX_2 i2(mt.pnums[i], mt.pnums[j]); 1045 i2.Sort(); 1046 int hval = edgenumber.Get(i2); 1047 if (hval > val) 1048 { 1049 val = hval; 1050 mt.markededge = 3 - i - j; 1051 } 1052 } 1053 } 1054 1055 1056 PrettyPrint(ostream & ost,const MarkedTri & mt)1057 void PrettyPrint(ostream & ost, const MarkedTri & mt) 1058 { 1059 ost << "MarkedTrig: " << endl; 1060 ost << " pnums = "; for (int i=0; i<3; i++) ost << mt.pnums[i] << " "; ost << endl; 1061 ost << " marked = " << mt.marked << ", markededge=" << mt.markededge << endl; 1062 for(int i=0; i<2; i++) 1063 for(int j=i+1; j<3; j++) 1064 if(mt.markededge == 3-i-j) 1065 ost << " marked edge pnums = " << mt.pnums[i] << " " << mt.pnums[j] << endl; 1066 } 1067 1068 PrettyPrint(ostream & ost,const MarkedQuad & mq)1069 void PrettyPrint(ostream & ost, const MarkedQuad & mq) 1070 { 1071 ost << "MarkedQuad: " << endl; 1072 ost << " pnums = "; for (int i=0; i<4; i++) ost << mq.pnums[i] << " "; ost << endl; 1073 ost << " marked = " << mq.marked << ", markededge=" << mq.markededge << endl; 1074 } 1075 1076 1077 1078 1079 BTDefineMarkedQuad(const Element2d & el,INDEX_2_CLOSED_HASHTABLE<int> & edgenumber,MarkedQuad & mq)1080 void BTDefineMarkedQuad (const Element2d & el, 1081 INDEX_2_CLOSED_HASHTABLE<int> & edgenumber, 1082 MarkedQuad & mq) 1083 { 1084 for (int i = 0; i < 4; i++) 1085 { 1086 mq.pnums[i] = el[i]; 1087 mq.pgeominfo[i] = el.GeomInfoPi (i+1); 1088 } 1089 Swap (mq.pnums[2], mq.pnums[3]); 1090 Swap (mq.pgeominfo[2], mq.pgeominfo[3]); 1091 1092 mq.marked = 0; 1093 mq.markededge = 0; 1094 mq.surfid = el.GetIndex(); 1095 } 1096 1097 1098 1099 1100 // mark elements due to local h BTMarkTets(T_MTETS & mtets,T_MPRISMS & mprisms,const Mesh & mesh)1101 int BTMarkTets (T_MTETS & mtets, 1102 T_MPRISMS & mprisms, 1103 const Mesh & mesh) 1104 { 1105 int marked = 0; 1106 1107 int np = mesh.GetNP(); 1108 Vector hv(np); 1109 for (int i = 0; i < np; i++) 1110 hv(i) = mesh.GetH (mesh.Point(i+1)); 1111 1112 double hfac = 1; 1113 1114 for (int step = 1; step <= 2; step++) 1115 { 1116 for (int i = 1; i <= mtets.Size(); i++) 1117 { 1118 double h = 0; 1119 1120 for (int j = 0; j < 3; j++) 1121 for (int k = j+1; k < 4; k++) 1122 { 1123 const Point<3> & p1 = mesh.Point (mtets.Get(i).pnums[j]); 1124 const Point<3> & p2 = mesh.Point (mtets.Get(i).pnums[k]); 1125 double hh = Dist2 (p1, p2); 1126 if (hh > h) h = hh; 1127 } 1128 h = sqrt (h); 1129 1130 double hshould = 1e10; 1131 for (int j = 0; j < 4; j++) 1132 { 1133 double hi = hv (mtets.Get(i).pnums[j]-1); 1134 if (hi < hshould) 1135 hshould = hi; 1136 } 1137 1138 1139 if (step == 1) 1140 { 1141 if (h / hshould > hfac) 1142 hfac = h / hshould; 1143 } 1144 else 1145 { 1146 if (h > hshould * hfac) 1147 { 1148 mtets.Elem(i).marked = 1; 1149 marked = 1; 1150 } 1151 else 1152 mtets.Elem(i).marked = 0; 1153 } 1154 1155 } 1156 for (int i = 1; i <= mprisms.Size(); i++) 1157 { 1158 double h = 0; 1159 1160 for (int j = 0; j < 2; j++) 1161 for (int k = j+1; k < 3; k++) 1162 { 1163 const Point<3> & p1 = mesh.Point (mprisms.Get(i).pnums[j]); 1164 const Point<3> & p2 = mesh.Point (mprisms.Get(i).pnums[k]); 1165 double hh = Dist2 (p1, p2); 1166 if (hh > h) h = hh; 1167 } 1168 h = sqrt (h); 1169 1170 double hshould = 1e10; 1171 for (int j = 0; j < 6; j++) 1172 { 1173 double hi = hv (mprisms.Get(i).pnums[j]-1); 1174 if (hi < hshould) 1175 hshould = hi; 1176 } 1177 1178 1179 if (step == 1) 1180 { 1181 if (h / hshould > hfac) 1182 hfac = h / hshould; 1183 } 1184 else 1185 { 1186 if (h > hshould * hfac) 1187 { 1188 mprisms.Elem(i).marked = 1; 1189 marked = 1; 1190 } 1191 else 1192 mprisms.Elem(i).marked = 0; 1193 } 1194 1195 } 1196 1197 1198 1199 if (step == 1) 1200 { 1201 if (hfac > 2) 1202 hfac /= 2; 1203 else 1204 hfac = 1; 1205 } 1206 1207 } 1208 return marked; 1209 } 1210 1211 1212 1213 1214 1215 1216 1217 1218 1219 1220 1221 1222 1223 BTBisectTet(const MarkedTet & oldtet,int newp,MarkedTet & newtet1,MarkedTet & newtet2)1224 void BTBisectTet (const MarkedTet & oldtet, int newp, 1225 MarkedTet & newtet1, MarkedTet & newtet2) 1226 { 1227 #ifdef DEBUG 1228 *testout << "bisect tet " << oldtet << endl; 1229 #endif 1230 1231 1232 // points vis a vis from tet-edge 1233 int vis1, vis2; 1234 vis1 = 0; 1235 while (vis1 == oldtet.tetedge1 || vis1 == oldtet.tetedge2) 1236 vis1++; 1237 vis2 = 6 - vis1 - oldtet.tetedge1 - oldtet.tetedge2; 1238 1239 1240 1241 1242 1243 // is tet of type P ? 1244 int istypep = 0; 1245 for (int i = 0; i < 4; i++) 1246 { 1247 int cnt = 0; 1248 for (int j = 0; j < 4; j++) 1249 if (oldtet.faceedges[j] == i) 1250 cnt++; 1251 if (cnt == 3) 1252 istypep = 1; 1253 } 1254 1255 1256 1257 for (int i = 0; i < 4; i++) 1258 { 1259 newtet1.pnums[i] = oldtet.pnums[i]; 1260 newtet2.pnums[i] = oldtet.pnums[i]; 1261 } 1262 newtet1.flagged = istypep && !oldtet.flagged; 1263 newtet2.flagged = istypep && !oldtet.flagged; 1264 1265 int nm = oldtet.marked - 1; 1266 if (nm < 0) nm = 0; 1267 newtet1.marked = nm; 1268 newtet2.marked = nm; 1269 1270 #ifdef DEBUG 1271 *testout << "newtet1,before = " << newtet1 << endl; 1272 *testout << "newtet2,before = " << newtet2 << endl; 1273 #endif 1274 1275 for (int i = 0; i < 4; i++) 1276 { 1277 if (i == oldtet.tetedge1) 1278 { 1279 newtet2.pnums[i] = newp; 1280 newtet2.faceedges[i] = oldtet.faceedges[i]; // inherited face 1281 newtet2.faceedges[vis1] = i; // cut faces 1282 newtet2.faceedges[vis2] = i; 1283 1284 int j = 0; 1285 while (j == i || j == oldtet.faceedges[i]) 1286 j++; 1287 int k = 6 - i - oldtet.faceedges[i] - j; 1288 newtet2.tetedge1 = j; // tet-edge 1289 newtet2.tetedge2 = k; 1290 1291 // new face: 1292 if (istypep && oldtet.flagged) 1293 { 1294 int hi = 6 - oldtet.tetedge1 - j - k; 1295 newtet2.faceedges[oldtet.tetedge2] = char(hi); 1296 } 1297 else 1298 newtet2.faceedges[oldtet.tetedge2] = oldtet.tetedge1; 1299 1300 #ifdef DEBUG 1301 *testout << "i = " << i << ", j = " << j << " k = " << k 1302 << " oldtet.tetedge1 = " << oldtet.tetedge1 1303 << " oldtet.tetedge2 = " << oldtet.tetedge2 1304 << " 6-oldtet.tetedge1-j-k = " << 6 - oldtet.tetedge1 - j - k 1305 << " 6-oldtet.tetedge1-j-k = " << short(6 - oldtet.tetedge1 - j - k) 1306 << endl; 1307 *testout << "vis1 = " << vis1 << ", vis2 = " << vis2 << endl; 1308 for (int j = 0; j < 4; j++) 1309 if (newtet2.faceedges[j] > 3) 1310 { 1311 *testout << "ERROR1" << endl; 1312 } 1313 #endif 1314 } 1315 1316 if (i == oldtet.tetedge2) 1317 { 1318 newtet1.pnums[i] = newp; 1319 newtet1.faceedges[i] = oldtet.faceedges[i]; // inherited face 1320 newtet1.faceedges[vis1] = i; 1321 newtet1.faceedges[vis2] = i; 1322 int j = 0; 1323 while (j == i || j == oldtet.faceedges[i]) 1324 j++; 1325 int k = 6 - i - oldtet.faceedges[i] - j; 1326 newtet1.tetedge1 = j; 1327 newtet1.tetedge2 = k; 1328 1329 // new face: 1330 if (istypep && oldtet.flagged) 1331 { 1332 int hi = 6 - oldtet.tetedge2 - j - k; 1333 newtet1.faceedges[oldtet.tetedge1] = char(hi); 1334 } 1335 else 1336 newtet1.faceedges[oldtet.tetedge1] = oldtet.tetedge2; 1337 1338 #ifdef DEBUG 1339 for (int j = 0; j < 4; j++) 1340 if (newtet2.faceedges[j] > 3) 1341 { 1342 *testout << "ERROR2" << endl; 1343 } 1344 #endif 1345 } 1346 } 1347 1348 newtet1.matindex = oldtet.matindex; 1349 newtet2.matindex = oldtet.matindex; 1350 newtet1.incorder = 0; 1351 newtet1.order = oldtet.order; 1352 newtet2.incorder = 0; 1353 newtet2.order = oldtet.order; 1354 1355 // *testout << "newtet1 = " << newtet1 << endl; 1356 // *testout << "newtet2 = " << newtet2 << endl; 1357 } 1358 1359 1360 1361 BTBisectPrism(const MarkedPrism & oldprism,int newp1,int newp2,MarkedPrism & newprism1,MarkedPrism & newprism2)1362 void BTBisectPrism (const MarkedPrism & oldprism, int newp1, int newp2, 1363 MarkedPrism & newprism1, MarkedPrism & newprism2) 1364 { 1365 for (int i = 0; i < 6; i++) 1366 { 1367 newprism1.pnums[i] = oldprism.pnums[i]; 1368 newprism2.pnums[i] = oldprism.pnums[i]; 1369 } 1370 1371 int pe1 = 0; 1372 if (pe1 == oldprism.markededge) 1373 pe1++; 1374 int pe2 = 3 - oldprism.markededge - pe1; 1375 1376 newprism1.pnums[pe2] = newp1; 1377 newprism1.pnums[pe2+3] = newp2; 1378 newprism1.markededge = pe2; 1379 newprism2.pnums[pe1] = newp1; 1380 newprism2.pnums[pe1+3] = newp2; 1381 newprism2.markededge = pe1; 1382 1383 newprism1.matindex = oldprism.matindex; 1384 newprism2.matindex = oldprism.matindex; 1385 1386 int nm = oldprism.marked - 1; 1387 if (nm < 0) nm = 0; 1388 newprism1.marked = nm; 1389 newprism2.marked = nm; 1390 1391 newprism1.incorder = 0; 1392 newprism1.order = oldprism.order; 1393 newprism2.incorder = 0; 1394 newprism2.order = oldprism.order; 1395 } 1396 1397 BTBisectIdentification(const MarkedIdentification & oldid,NgArray<PointIndex> & newp,MarkedIdentification & newid1,MarkedIdentification & newid2)1398 void BTBisectIdentification (const MarkedIdentification & oldid, 1399 NgArray<PointIndex> & newp, 1400 MarkedIdentification & newid1, 1401 MarkedIdentification & newid2) 1402 { 1403 for(int i=0; i<2*oldid.np; i++) 1404 { 1405 newid1.pnums[i] = oldid.pnums[i]; 1406 newid2.pnums[i] = oldid.pnums[i]; 1407 } 1408 newid1.np = newid2.np = oldid.np; 1409 1410 if(oldid.np == 2) 1411 { 1412 newid1.pnums[1] = newp[0]; 1413 newid2.pnums[0] = newp[0]; 1414 newid1.pnums[3] = newp[1]; 1415 newid2.pnums[2] = newp[1]; 1416 newid1.markededge = 0; 1417 newid2.markededge = 0; 1418 } 1419 1420 if(oldid.np == 3) 1421 { 1422 newid1.pnums[(oldid.markededge+1)%3] = newp[0]; 1423 newid1.pnums[(oldid.markededge+1)%3+3] = newp[1]; 1424 newid1.markededge = (oldid.markededge+2)%3; 1425 1426 newid2.pnums[oldid.markededge] = newp[0]; 1427 newid2.pnums[oldid.markededge+3] = newp[1]; 1428 newid2.markededge = (oldid.markededge+1)%3; 1429 } 1430 else if(oldid.np == 4) 1431 { 1432 newid1.pnums[(oldid.markededge+1)%4] = newp[0]; 1433 newid1.pnums[(oldid.markededge+2)%4] = newp[2]; 1434 newid1.pnums[(oldid.markededge+1)%4+4] = newp[1]; 1435 newid1.pnums[(oldid.markededge+2)%4+4] = newp[3]; 1436 newid1.markededge = (oldid.markededge+3)%4; 1437 1438 newid2.pnums[oldid.markededge] = newp[0]; 1439 newid2.pnums[(oldid.markededge+3)%4] = newp[2]; 1440 newid2.pnums[oldid.markededge+4] = newp[1]; 1441 newid2.pnums[(oldid.markededge+3)%4+4] = newp[3]; 1442 newid2.markededge = (oldid.markededge+1)%4; 1443 } 1444 1445 1446 int nm = oldid.marked - 1; 1447 if (nm < 0) nm = 0; 1448 newid1.marked = newid2.marked = nm; 1449 1450 newid1.incorder = newid2.incorder = 0; 1451 newid1.order = newid2.order = oldid.order; 1452 } 1453 1454 1455 BTBisectTri(const MarkedTri & oldtri,int newp,const PointGeomInfo & newpgi,MarkedTri & newtri1,MarkedTri & newtri2)1456 void BTBisectTri (const MarkedTri & oldtri, int newp, const PointGeomInfo & newpgi, 1457 MarkedTri & newtri1, MarkedTri & newtri2) 1458 { 1459 for (int i = 0; i < 3; i++) 1460 { 1461 newtri1.pnums[i] = oldtri.pnums[i]; 1462 newtri1.pgeominfo[i] = oldtri.pgeominfo[i]; 1463 newtri2.pnums[i] = oldtri.pnums[i]; 1464 newtri2.pgeominfo[i] = oldtri.pgeominfo[i]; 1465 } 1466 1467 int pe1 = 0; 1468 if (pe1 == oldtri.markededge) 1469 pe1++; 1470 int pe2 = 3 - oldtri.markededge - pe1; 1471 1472 newtri1.pnums[pe2] = newp; 1473 newtri1.pgeominfo[pe2] = newpgi; 1474 newtri1.markededge = pe2; 1475 1476 newtri2.pnums[pe1] = newp; 1477 newtri2.pgeominfo[pe1] = newpgi; 1478 newtri2.markededge = pe1; 1479 1480 1481 newtri1.surfid = oldtri.surfid; 1482 newtri2.surfid = oldtri.surfid; 1483 1484 int nm = oldtri.marked - 1; 1485 if (nm < 0) nm = 0; 1486 newtri1.marked = nm; 1487 newtri2.marked = nm; 1488 1489 newtri1.incorder = 0; 1490 newtri1.order = oldtri.order; 1491 newtri2.incorder = 0; 1492 newtri2.order = oldtri.order; 1493 } 1494 1495 BTBisectQuad(const MarkedQuad & oldquad,int newp1,const PointGeomInfo & npgi1,int newp2,const PointGeomInfo & npgi2,MarkedQuad & newquad1,MarkedQuad & newquad2)1496 void BTBisectQuad (const MarkedQuad & oldquad, 1497 int newp1, const PointGeomInfo & npgi1, 1498 int newp2, const PointGeomInfo & npgi2, 1499 MarkedQuad & newquad1, MarkedQuad & newquad2) 1500 { 1501 for (int i = 0; i < 4; i++) 1502 { 1503 newquad1.pnums[i] = oldquad.pnums[i]; 1504 newquad1.pgeominfo[i] = oldquad.pgeominfo[i]; 1505 newquad2.pnums[i] = oldquad.pnums[i]; 1506 newquad2.pgeominfo[i] = oldquad.pgeominfo[i]; 1507 } 1508 1509 /* if (oldquad.marked==1) // he/sz: 2d quads or 3d prism 1510 { 1511 newquad1.pnums[1] = newp1; 1512 newquad1.pgeominfo[1] = npgi1; 1513 newquad1.pnums[3] = newp2; 1514 newquad1.pgeominfo[3] = npgi2; 1515 1516 newquad2.pnums[0] = newp1; 1517 newquad2.pgeominfo[0] = npgi1; 1518 newquad2.pnums[2] = newp2; 1519 newquad2.pgeominfo[2] = npgi2; 1520 } 1521 1522 else if (oldquad.marked==2) // he/sz: 2d quads only 1523 { 1524 newquad1.pnums[0] = newp1; 1525 newquad1.pnums[1] = newp2; 1526 newquad1.pnums[3] = oldquad.pnums[2]; 1527 newquad1.pnums[2] = oldquad.pnums[0]; 1528 newquad1.pgeominfo[0] = npgi1; 1529 newquad1.pgeominfo[1] = npgi2; 1530 newquad1.pgeominfo[3] = oldquad.pgeominfo[2]; 1531 newquad1.pgeominfo[2] = oldquad.pgeominfo[0]; 1532 1533 newquad2.pnums[0] = newp2; 1534 newquad2.pnums[1] = newp1; 1535 newquad2.pnums[3] = oldquad.pnums[1]; 1536 newquad2.pnums[2] = oldquad.pnums[3]; 1537 newquad2.pgeominfo[0] = npgi2; 1538 newquad2.pgeominfo[1] = npgi1; 1539 newquad2.pgeominfo[3] = oldquad.pgeominfo[1]; 1540 newquad2.pgeominfo[2] = oldquad.pgeominfo[3]; 1541 } 1542 1543 */ 1544 1545 if (oldquad.markededge==0 || oldquad.markededge==2) 1546 { 1547 newquad1.pnums[1] = newp1; 1548 newquad1.pgeominfo[1] = npgi1; 1549 newquad1.pnums[3] = newp2; 1550 newquad1.pgeominfo[3] = npgi2; 1551 1552 newquad2.pnums[0] = newp1; 1553 newquad2.pgeominfo[0] = npgi1; 1554 newquad2.pnums[2] = newp2; 1555 newquad2.pgeominfo[2] = npgi2; 1556 } 1557 else // 1 || 3 1558 { 1559 newquad1.pnums[2] = newp1; 1560 newquad1.pgeominfo[2] = npgi1; 1561 newquad1.pnums[3] = newp2; 1562 newquad1.pgeominfo[3] = npgi2; 1563 1564 newquad2.pnums[0] = newp1; 1565 newquad2.pgeominfo[0] = npgi1; 1566 newquad2.pnums[1] = newp2; 1567 newquad2.pgeominfo[1] = npgi2; 1568 } 1569 newquad1.surfid = oldquad.surfid; 1570 newquad2.surfid = oldquad.surfid; 1571 1572 int nm = oldquad.marked - 1; 1573 if (nm < 0) nm = 0; 1574 1575 newquad1.marked = nm; 1576 newquad2.marked = nm; 1577 1578 if (nm==1) 1579 { 1580 newquad1.markededge=1; 1581 newquad2.markededge=1; 1582 } 1583 else 1584 { 1585 newquad1.markededge=0; 1586 newquad2.markededge=0; 1587 } 1588 1589 } 1590 1591 MarkHangingIdentifications(T_MIDS & mids,const INDEX_2_CLOSED_HASHTABLE<PointIndex> & cutedges)1592 int MarkHangingIdentifications(T_MIDS & mids, 1593 const INDEX_2_CLOSED_HASHTABLE<PointIndex> & cutedges) 1594 { 1595 int hanging = 0; 1596 for (int i = 1; i <= mids.Size(); i++) 1597 { 1598 if (mids.Elem(i).marked) 1599 { 1600 hanging = 1; 1601 continue; 1602 } 1603 1604 const int np = mids.Get(i).np; 1605 for(int j = 0; j < np; j++) 1606 { 1607 INDEX_2 edge1(mids.Get(i).pnums[j], 1608 mids.Get(i).pnums[(j+1) % np]); 1609 INDEX_2 edge2(mids.Get(i).pnums[j+np], 1610 mids.Get(i).pnums[((j+1) % np) + np]); 1611 1612 edge1.Sort(); 1613 edge2.Sort(); 1614 if (cutedges.Used (edge1) || 1615 cutedges.Used (edge2)) 1616 { 1617 mids.Elem(i).marked = 1; 1618 hanging = 1; 1619 } 1620 } 1621 } 1622 1623 return hanging; 1624 } 1625 1626 1627 /* 1628 void IdentifyCutEdges(Mesh & mesh, 1629 INDEX_2_CLOSED_HASHTABLE<int> & cutedges) 1630 { 1631 int i,j,k; 1632 1633 NgArray< NgArray<int,PointIndex::BASE>* > idmaps; 1634 for(i=1; i<=mesh.GetIdentifications().GetMaxNr(); i++) 1635 { 1636 idmaps.Append(new NgArray<int,PointIndex::BASE>); 1637 mesh.GetIdentifications().GetMap(i,*idmaps.Last()); 1638 } 1639 1640 1641 1642 for(SurfaceElementIndex sei = 0; sei < mesh.GetNSE(); sei++) 1643 { 1644 const Element2d & el2d = mesh[sei]; 1645 1646 for(i = 0; i < el2d.GetNP(); i++) 1647 { 1648 INDEX_2 e1(el2d[i], el2d[(i+1) % el2d.GetNP()]); 1649 e1.Sort(); 1650 1651 if(!cutedges.Used(e1)) 1652 continue; 1653 1654 1655 for(k = 0; k < idmaps.Size(); k++) 1656 { 1657 INDEX_2 e2((*idmaps[k])[e1.I1()], 1658 (*idmaps[k])[e1.I2()]); 1659 1660 if(e2.I1() == 0 || e2.I2() == 0 || 1661 e1.I1() == e2.I1() || e1.I2() == e2.I2()) 1662 continue; 1663 1664 e2.Sort(); 1665 1666 if(cutedges.Used(e2)) 1667 continue; 1668 1669 Point3d np = Center(mesh.Point(e2.I1()), 1670 mesh.Point(e2.I2())); 1671 int newp = mesh.AddPoint(np); 1672 cutedges.Set(e2,newp); 1673 (*testout) << "DAAA" << endl; 1674 } 1675 } 1676 } 1677 1678 1679 for(i=0; i<idmaps.Size(); i++) 1680 delete idmaps[i]; 1681 idmaps.DeleteAll(); 1682 } 1683 */ 1684 1685 MarkHangingTets(T_MTETS & mtets,const INDEX_2_CLOSED_HASHTABLE<PointIndex> & cutedges,NgTaskManager tm)1686 int MarkHangingTets (T_MTETS & mtets, 1687 const INDEX_2_CLOSED_HASHTABLE<PointIndex> & cutedges, 1688 NgTaskManager tm) 1689 { 1690 static int timer = NgProfiler::CreateTimer ("MarkHangingTets"); 1691 NgProfiler::RegionTimer reg (timer); 1692 1693 int hanging = 0; 1694 // for (int i = 1; i <= mtets.Size(); i++) 1695 ParallelForRange 1696 (tm, mtets.Size(), [&] (size_t begin, size_t end) 1697 { 1698 bool my_hanging = false; 1699 for (size_t i = begin; i < end; i++) 1700 { 1701 MarkedTet & teti = mtets[i]; 1702 1703 if (teti.marked) 1704 { 1705 my_hanging = true; 1706 continue; 1707 } 1708 1709 for (int j = 0; j < 3; j++) 1710 for (int k = j+1; k < 4; k++) 1711 { 1712 INDEX_2 edge(teti.pnums[j], 1713 teti.pnums[k]); 1714 edge.Sort(); 1715 if (cutedges.Used (edge)) 1716 { 1717 teti.marked = 1; 1718 my_hanging = true; 1719 } 1720 } 1721 } 1722 if (my_hanging) hanging = true; 1723 }); 1724 1725 return hanging; 1726 } 1727 1728 1729 MarkHangingPrisms(T_MPRISMS & mprisms,const INDEX_2_CLOSED_HASHTABLE<PointIndex> & cutedges)1730 int MarkHangingPrisms (T_MPRISMS & mprisms, 1731 const INDEX_2_CLOSED_HASHTABLE<PointIndex> & cutedges) 1732 { 1733 int hanging = 0; 1734 for (int i = 1; i <= mprisms.Size(); i++) 1735 { 1736 if (mprisms.Elem(i).marked) 1737 { 1738 hanging = 1; 1739 continue; 1740 } 1741 1742 for (int j = 0; j < 2; j++) 1743 for (int k = j+1; k < 3; k++) 1744 { 1745 INDEX_2 edge1(mprisms.Get(i).pnums[j], 1746 mprisms.Get(i).pnums[k]); 1747 INDEX_2 edge2(mprisms.Get(i).pnums[j+3], 1748 mprisms.Get(i).pnums[k+3]); 1749 edge1.Sort(); 1750 edge2.Sort(); 1751 if (cutedges.Used (edge1) || 1752 cutedges.Used (edge2)) 1753 { 1754 mprisms.Elem(i).marked = 1; 1755 hanging = 1; 1756 } 1757 } 1758 } 1759 return hanging; 1760 } 1761 1762 1763 MarkHangingTris(T_MTRIS & mtris,const INDEX_2_CLOSED_HASHTABLE<PointIndex> & cutedges,NgTaskManager tm)1764 bool MarkHangingTris (T_MTRIS & mtris, 1765 const INDEX_2_CLOSED_HASHTABLE<PointIndex> & cutedges, 1766 NgTaskManager tm) 1767 { 1768 bool hanging = false; 1769 // for (int i = 1; i <= mtris.Size(); i++) 1770 // for (auto & tri : mtris) 1771 ParallelForRange 1772 (tm, mtris.Size(), [&] (size_t begin, size_t end) 1773 { 1774 bool my_hanging = false; 1775 for (size_t i = begin; i < end; i++) 1776 { 1777 auto & tri = mtris[i]; 1778 if (tri.marked) 1779 { 1780 my_hanging = true; 1781 continue; 1782 } 1783 for (int j = 0; j < 2; j++) 1784 for (int k = j+1; k < 3; k++) 1785 { 1786 INDEX_2 edge(tri.pnums[j], 1787 tri.pnums[k]); 1788 edge.Sort(); 1789 if (cutedges.Used (edge)) 1790 { 1791 tri.marked = 1; 1792 my_hanging = true; 1793 } 1794 } 1795 } 1796 if (my_hanging) hanging = true; 1797 }); 1798 return hanging; 1799 } 1800 1801 1802 MarkHangingQuads(T_MQUADS & mquads,const INDEX_2_CLOSED_HASHTABLE<PointIndex> & cutedges)1803 int MarkHangingQuads (T_MQUADS & mquads, 1804 const INDEX_2_CLOSED_HASHTABLE<PointIndex> & cutedges) 1805 { 1806 int hanging = 0; 1807 for (int i = 1; i <= mquads.Size(); i++) 1808 { 1809 if (mquads.Elem(i).marked) 1810 { 1811 hanging = 1; 1812 continue; 1813 } 1814 1815 INDEX_2 edge1(mquads.Get(i).pnums[0], 1816 mquads.Get(i).pnums[1]); 1817 INDEX_2 edge2(mquads.Get(i).pnums[2], 1818 mquads.Get(i).pnums[3]); 1819 edge1.Sort(); 1820 edge2.Sort(); 1821 if (cutedges.Used (edge1) || 1822 cutedges.Used (edge2)) 1823 { 1824 mquads.Elem(i).marked = 1; 1825 mquads.Elem(i).markededge = 0; 1826 hanging = 1; 1827 continue; 1828 } 1829 1830 // he/sz: second case: split horizontally 1831 INDEX_2 edge3(mquads.Get(i).pnums[1], 1832 mquads.Get(i).pnums[3]); 1833 INDEX_2 edge4(mquads.Get(i).pnums[2], 1834 mquads.Get(i).pnums[0]); 1835 1836 edge3.Sort(); 1837 edge4.Sort(); 1838 if (cutedges.Used (edge3) || 1839 cutedges.Used (edge4)) 1840 { 1841 mquads.Elem(i).marked = 1; 1842 mquads.Elem(i).markededge = 1; 1843 hanging = 1; 1844 continue; 1845 } 1846 1847 } 1848 return hanging; 1849 } 1850 1851 1852 ConnectToNodeRec(int node,int tonode,const TABLE<int> & conto,NgArray<int> & connecttonode)1853 void ConnectToNodeRec (int node, int tonode, 1854 const TABLE<int> & conto, NgArray<int> & connecttonode) 1855 { 1856 // (*testout) << "connect " << node << " to " << tonode << endl; 1857 for (int i = 1; i <= conto.EntrySize(node); i++) 1858 { 1859 int n2 = conto.Get(node, i); 1860 if (!connecttonode.Get(n2)) 1861 { 1862 connecttonode.Elem(n2) = tonode; 1863 ConnectToNodeRec (n2, tonode, conto, connecttonode); 1864 } 1865 } 1866 } 1867 1868 1869 1870 1871 T_MTETS mtets; 1872 T_MPRISMS mprisms; 1873 T_MIDS mids; 1874 T_MTRIS mtris; 1875 T_MQUADS mquads; 1876 1877 WriteMarkedElements(ostream & ost)1878 void WriteMarkedElements(ostream & ost) 1879 { 1880 ost << "Marked Elements\n"; 1881 1882 ost << mtets.Size() << "\n"; 1883 for(int i=0; i<mtets.Size(); i++) 1884 ost << mtets[i]; 1885 1886 ost << mprisms.Size() << "\n"; 1887 for(int i=0; i<mprisms.Size(); i++) 1888 ost << mprisms[i]; 1889 1890 ost << mids.Size() << "\n"; 1891 for(int i=0; i<mids.Size(); i++) 1892 ost << mids[i]; 1893 1894 ost << mtris.Size() << "\n"; 1895 for(int i=0; i<mtris.Size(); i++) 1896 ost << mtris[i]; 1897 1898 ost << mquads.Size() << "\n"; 1899 for(int i=0; i<mquads.Size(); i++) 1900 ost << mquads[i]; 1901 ost << endl; 1902 } 1903 ReadMarkedElements(istream & ist,const Mesh & mesh)1904 bool ReadMarkedElements(istream & ist, const Mesh & mesh) 1905 { 1906 string auxstring(""); 1907 if(ist) 1908 ist >> auxstring; 1909 1910 if(auxstring != "Marked") 1911 return false; 1912 1913 if(ist) 1914 ist >> auxstring; 1915 1916 if(auxstring != "Elements") 1917 return false; 1918 1919 int size; 1920 1921 ist >> size; 1922 mtets.SetSize(size); 1923 for(int i=0; i<size; i++) 1924 { 1925 ist >> mtets[i]; 1926 if(mtets[i].pnums[0] > mesh.GetNV() || 1927 mtets[i].pnums[1] > mesh.GetNV() || 1928 mtets[i].pnums[2] > mesh.GetNV() || 1929 mtets[i].pnums[3] > mesh.GetNV()) 1930 return false; 1931 } 1932 1933 ist >> size; 1934 mprisms.SetSize(size); 1935 for(int i=0; i<size; i++) 1936 ist >> mprisms[i]; 1937 1938 ist >> size; 1939 mids.SetSize(size); 1940 for(int i=0; i<size; i++) 1941 ist >> mids[i]; 1942 1943 ist >> size; 1944 mtris.SetSize(size); 1945 for(int i=0; i<size; i++) 1946 ist >> mtris[i]; 1947 1948 ist >> size; 1949 mquads.SetSize(size); 1950 for(int i=0; i<size; i++) 1951 ist >> mquads[i]; 1952 1953 return true; 1954 } 1955 1956 1957 1958 1959 BisectTetsCopyMesh(Mesh & mesh,const NetgenGeometry *,BisectionOptions & opt,const NgArray<NgArray<int,PointIndex::BASE> * > & idmaps,const string & refinfofile)1960 void BisectTetsCopyMesh (Mesh & mesh, const NetgenGeometry *, 1961 BisectionOptions & opt, 1962 const NgArray< NgArray<int,PointIndex::BASE>* > & idmaps, 1963 const string & refinfofile) 1964 { 1965 if (mesh.GetDimension() < 2) 1966 throw Exception ("Mesh bisection is available in 2D and 3D"); 1967 // mtets.SetName ("bisection, tets"); 1968 // mprisms.SetName ("bisection, prisms"); 1969 // mtris.SetName ("bisection, trigs"); 1970 // nmquads.SetName ("bisection, quads"); 1971 // mids.SetName ("bisection, identifications"); 1972 1973 //int np = mesh.GetNP(); 1974 int ne = mesh.GetNE(); 1975 int nse = mesh.GetNSE(); 1976 1977 /* 1978 if (mtets.Size() + mprisms.Size() == mesh.GetNE()) 1979 return; 1980 */ 1981 1982 bool readok = false; 1983 1984 if(refinfofile != "") 1985 { 1986 PrintMessage(3,"Reading marked-element information from \"",refinfofile,"\""); 1987 ifstream ist(refinfofile.c_str()); 1988 1989 readok = ReadMarkedElements(ist,mesh); 1990 1991 ist.close(); 1992 } 1993 1994 if(!readok) 1995 { 1996 PrintMessage(3,"resetting marked-element information"); 1997 mtets.SetSize(0); 1998 mprisms.SetSize(0); 1999 mids.SetSize(0); 2000 mtris.SetSize(0); 2001 mquads.SetSize(0); 2002 2003 2004 INDEX_2_HASHTABLE<int> shortedges(100); 2005 for (int i = 1; i <= ne; i++) 2006 { 2007 const Element & el = mesh.VolumeElement(i); 2008 if (el.GetType() == PRISM || 2009 el.GetType() == PRISM12) 2010 { 2011 for (int j = 1; j <= 3; j++) 2012 { 2013 INDEX_2 se(el.PNum(j), el.PNum(j+3)); 2014 se.Sort(); 2015 shortedges.Set (se, 1); 2016 } 2017 } 2018 } 2019 2020 2021 2022 // INDEX_2_HASHTABLE<int> edgenumber(np); 2023 INDEX_2_CLOSED_HASHTABLE<int> edgenumber(9*ne+4*nse); 2024 2025 BTSortEdges (mesh, idmaps, edgenumber); 2026 2027 2028 for (int i = 1; i <= ne; i++) 2029 { 2030 const Element & el = mesh.VolumeElement(i); 2031 2032 switch (el.GetType()) 2033 { 2034 case TET: 2035 case TET10: 2036 { 2037 // if tet has short edge, it is handled as degenerated prism 2038 2039 int foundse = 0; 2040 for (int j = 1; j <= 3; j++) 2041 for (int k = j+1; k <= 4; k++) 2042 { 2043 INDEX_2 se(el.PNum(j), el.PNum(k)); 2044 se.Sort(); 2045 if (shortedges.Used (se)) 2046 { 2047 // cout << "tet converted to prism" << endl; 2048 2049 foundse = 1; 2050 int p3 = 1; 2051 while (p3 == j || p3 == k) 2052 p3++; 2053 int p4 = 10 - j - k - p3; 2054 2055 // even permutation ? 2056 int pi[4]; 2057 pi[0] = j; 2058 pi[1] = k; 2059 pi[2] = p3; 2060 pi[3] = p4; 2061 int cnt = 0; 2062 for (int l = 1; l <= 4; l++) 2063 for (int m = 0; m < 3; m++) 2064 if (pi[m] > pi[m+1]) 2065 { 2066 Swap (pi[m], pi[m+1]); 2067 cnt++; 2068 } 2069 if (cnt % 2) 2070 Swap (p3, p4); 2071 2072 Element hel = el; 2073 hel.PNum(1) = el.PNum(j); 2074 hel.PNum(2) = el.PNum(k); 2075 hel.PNum(3) = el.PNum(p3); 2076 hel.PNum(4) = el.PNum(p4); 2077 2078 MarkedPrism mp; 2079 BTDefineMarkedPrism (hel, edgenumber, mp); 2080 mp.matindex = el.GetIndex(); 2081 mprisms.Append (mp); 2082 } 2083 } 2084 if (!foundse) 2085 { 2086 MarkedTet mt; 2087 BTDefineMarkedTet (el, edgenumber, mt); 2088 mt.matindex = el.GetIndex(); 2089 mtets.Append (mt); 2090 } 2091 break; 2092 } 2093 case PYRAMID: 2094 { 2095 // eventually rotate 2096 MarkedPrism mp; 2097 2098 INDEX_2 se(el.PNum(1), el.PNum(2)); 2099 se.Sort(); 2100 if (shortedges.Used (se)) 2101 { 2102 Element hel = el; 2103 hel.PNum(1) = el.PNum(2); 2104 hel.PNum(2) = el.PNum(3); 2105 hel.PNum(3) = el.PNum(4); 2106 hel.PNum(4) = el.PNum(1); 2107 BTDefineMarkedPrism (hel, edgenumber, mp); 2108 } 2109 else 2110 { 2111 BTDefineMarkedPrism (el, edgenumber, mp); 2112 } 2113 2114 mp.matindex = el.GetIndex(); 2115 mprisms.Append (mp); 2116 break; 2117 } 2118 case PRISM: 2119 case PRISM12: 2120 { 2121 MarkedPrism mp; 2122 BTDefineMarkedPrism (el, edgenumber, mp); 2123 mp.matindex = el.GetIndex(); 2124 mprisms.Append (mp); 2125 break; 2126 } 2127 default: 2128 throw NgException("Bisect, element type not handled in switch, 4"); 2129 } 2130 } 2131 2132 for (int i = 1; i <= nse; i++) 2133 { 2134 const Element2d & el = mesh.SurfaceElement(i); 2135 if (el.GetType() == TRIG || 2136 el.GetType() == TRIG6) 2137 { 2138 MarkedTri mt; 2139 BTDefineMarkedTri (el, edgenumber, mt); 2140 mtris.Append (mt); 2141 } 2142 else 2143 { 2144 MarkedQuad mq; 2145 BTDefineMarkedQuad (el, edgenumber, mq); 2146 mquads.Append (mq); 2147 } 2148 2149 if(mesh.GetDimension() == 3) 2150 { 2151 MarkedIdentification mi; 2152 for(int j=0; j<idmaps.Size(); j++) 2153 if(BTDefineMarkedId(el, edgenumber, *idmaps[j], mi)) 2154 mids.Append(mi); 2155 } 2156 } 2157 if(mesh.GetDimension() == 2) 2158 { 2159 for (SegmentIndex j=0; j<mesh.GetNSeg(); j++) 2160 { 2161 auto seg = mesh[j]; 2162 for (auto map : idmaps) 2163 { 2164 if (seg[0] > 0 && seg[1] > 0 && (*map)[seg[0]] && (*map)[seg[1]]) 2165 { 2166 MarkedIdentification mi; 2167 mi.np = 2; 2168 mi.pnums[0] = seg[0]; 2169 mi.pnums[1] = seg[1]; 2170 mi.pnums[2] = (*map)[seg[0]]; 2171 mi.pnums[3] = (*map)[seg[1]]; 2172 auto min1 = mi.pnums[0] < mi.pnums[1] ? mi.pnums[0] : mi.pnums[1]; 2173 auto min2 = mi.pnums[2] < mi.pnums[3] ? mi.pnums[2] : mi.pnums[3]; 2174 if (min1 > min2) 2175 continue; 2176 mi.marked = 0; 2177 mi.markededge = 0; 2178 mi.incorder = 0; 2179 mids.Append(mi); 2180 } 2181 } 2182 } 2183 } 2184 } 2185 2186 2187 2188 2189 mesh.mlparentelement.SetSize(ne); 2190 for (int i = 1; i <= ne; i++) 2191 mesh.mlparentelement.Elem(i) = 0; 2192 mesh.mlparentsurfaceelement.SetSize(nse); 2193 for (int i = 1; i <= nse; i++) 2194 mesh.mlparentsurfaceelement.Elem(i) = 0; 2195 2196 if (printmessage_importance>0) 2197 { 2198 ostringstream str1,str2; 2199 str1 << "copied " << mtets.Size() << " tets, " << mprisms.Size() << " prisms"; 2200 str2 << " " << mtris.Size() << " trigs, " << mquads.Size() << " quads"; 2201 2202 PrintMessage(4,str1.str()); 2203 PrintMessage(4,str2.str()); 2204 } 2205 } 2206 2207 2208 /* 2209 void UpdateEdgeMarks2(Mesh & mesh, 2210 const NgArray< NgArray<int,PointIndex::BASE>* > & idmaps) 2211 { 2212 NgArray< NgArray<MarkedTet>*,PointIndex::BASE > mtets_old(mesh.GetNP()); 2213 NgArray< NgArray<MarkedPrism>*,PointIndex::BASE > mprisms_old(mesh.GetNP()); 2214 NgArray< NgArray<MarkedIdentification>*,PointIndex::BASE > mids_old(mesh.GetNP()); 2215 NgArray< NgArray<MarkedTri>*,PointIndex::BASE > mtris_old(mesh.GetNP()); 2216 NgArray< NgArray<MarkedQuad>*,PointIndex::BASE > mquads_old(mesh.GetNP()); 2217 2218 for(int i=PointIndex::BASE; i<mesh.GetNP()+PointIndex::BASE; i++) 2219 mtets_old[i] = new NgArray<MarkedTet>; 2220 for(int i=PointIndex::BASE; i<mesh.GetNP()+PointIndex::BASE; i++) 2221 mprisms_old[i] = new NgArray<MarkedPrism>; 2222 for(int i=PointIndex::BASE; i<mesh.GetNP()+PointIndex::BASE; i++) 2223 mids_old[i] = new NgArray<MarkedIdentification>; 2224 for(int i=PointIndex::BASE; i<mesh.GetNP()+PointIndex::BASE; i++) 2225 mtris_old[i] = new NgArray<MarkedTri>; 2226 for(int i=PointIndex::BASE; i<mesh.GetNP()+PointIndex::BASE; i++) 2227 mquads_old[i] = new NgArray<MarkedQuad>; 2228 2229 for(int i=0; i<mtets.Size(); i++) 2230 mtets_old[mtets[i].pnums[0]]->Append(mtets[i]); 2231 for(int i=0; i<mprisms.Size(); i++) 2232 mprisms_old[mprisms[i].pnums[0]]->Append(mprisms[i]); 2233 for(int i=0; i<mids.Size(); i++) 2234 mids_old[mids[i].pnums[0]]->Append(mids[i]); 2235 for(int i=0; i<mtris.Size(); i++) 2236 { 2237 (*testout) << "i " << i << endl; 2238 (*testout) << "mtris[i] " << mtris[i].pnums[0] << " " << mtris[i].pnums[1] << " " << mtris[i].pnums[2] << endl; 2239 mtris_old[mtris[i].pnums[0]]->Append(mtris[i]); 2240 } 2241 for(int i=0; i<mquads.Size(); i++) 2242 mquads_old[mquads[i].pnums[0]]->Append(mquads[i]); 2243 2244 2245 2246 int np = mesh.GetNP(); 2247 int ne = mesh.GetNE(); 2248 int nse = mesh.GetNSE(); 2249 int i, j, k, l, m; 2250 2251 2252 // if (mtets.Size() + mprisms.Size() == mesh.GetNE()) 2253 // return; 2254 2255 2256 2257 mtets.SetSize(0); 2258 mprisms.SetSize(0); 2259 mids.SetSize(0); 2260 mtris.SetSize(0); 2261 mquads.SetSize(0); 2262 2263 2264 INDEX_2_HASHTABLE<int> shortedges(100); 2265 for (i = 1; i <= ne; i++) 2266 { 2267 const Element & el = mesh.VolumeElement(i); 2268 if (el.GetType() == PRISM || 2269 el.GetType() == PRISM12) 2270 { 2271 for (j = 1; j <= 3; j++) 2272 { 2273 INDEX_2 se(el.PNum(j), el.PNum(j+3)); 2274 se.Sort(); 2275 shortedges.Set (se, 1); 2276 } 2277 } 2278 } 2279 2280 2281 2282 // INDEX_2_HASHTABLE<int> edgenumber(np); 2283 INDEX_2_CLOSED_HASHTABLE<int> edgenumber(9*ne+4*nse); 2284 2285 BTSortEdges (mesh, idmaps, edgenumber); 2286 2287 2288 for (i = 1; i <= ne; i++) 2289 { 2290 const Element & el = mesh.VolumeElement(i); 2291 2292 switch (el.GetType()) 2293 { 2294 case TET: 2295 case TET10: 2296 { 2297 // if tet has short edge, it is handled as degenerated prism 2298 2299 int foundse = 0; 2300 for (j = 1; j <= 3; j++) 2301 for (k = j+1; k <= 4; k++) 2302 { 2303 INDEX_2 se(el.PNum(j), el.PNum(k)); 2304 se.Sort(); 2305 if (shortedges.Used (se)) 2306 { 2307 // cout << "tet converted to prism" << endl; 2308 2309 foundse = 1; 2310 int p3 = 1; 2311 while (p3 == j || p3 == k) 2312 p3++; 2313 int p4 = 10 - j - k - p3; 2314 2315 // even permutation ? 2316 int pi[4]; 2317 pi[0] = j; 2318 pi[1] = k; 2319 pi[2] = p3; 2320 pi[3] = p4; 2321 int cnt = 0; 2322 for (l = 1; l <= 4; l++) 2323 for (m = 0; m < 3; m++) 2324 if (pi[m] > pi[m+1]) 2325 { 2326 Swap (pi[m], pi[m+1]); 2327 cnt++; 2328 } 2329 if (cnt % 2) 2330 Swap (p3, p4); 2331 2332 Element hel = el; 2333 hel.PNum(1) = el.PNum(j); 2334 hel.PNum(2) = el.PNum(k); 2335 hel.PNum(3) = el.PNum(p3); 2336 hel.PNum(4) = el.PNum(p4); 2337 2338 MarkedPrism mp; 2339 2340 BTDefineMarkedPrism (hel, edgenumber, mp); 2341 mp.matindex = el.GetIndex(); 2342 mprisms.Append (mp); 2343 } 2344 } 2345 if (!foundse) 2346 { 2347 MarkedTet mt; 2348 2349 int oldind = -1; 2350 for(l = 0; oldind < 0 && l<mtets_old[el[0]]->Size(); l++) 2351 if(el[1] == (*mtets_old[el[0]])[l].pnums[1] && 2352 el[2] == (*mtets_old[el[0]])[l].pnums[2] && 2353 el[3] == (*mtets_old[el[0]])[l].pnums[3]) 2354 oldind = l; 2355 2356 if(oldind >= 0) 2357 mtets.Append((*mtets_old[el[0]])[oldind]); 2358 else 2359 { 2360 BTDefineMarkedTet (el, edgenumber, mt); 2361 mt.matindex = el.GetIndex(); 2362 mtets.Append (mt); 2363 } 2364 } 2365 break; 2366 } 2367 case PYRAMID: 2368 { 2369 // eventually rotate 2370 MarkedPrism mp; 2371 2372 INDEX_2 se(el.PNum(1), el.PNum(2)); 2373 se.Sort(); 2374 if (shortedges.Used (se)) 2375 { 2376 Element hel = el; 2377 hel.PNum(1) = el.PNum(2); 2378 hel.PNum(2) = el.PNum(3); 2379 hel.PNum(3) = el.PNum(4); 2380 hel.PNum(4) = el.PNum(1); 2381 BTDefineMarkedPrism (hel, edgenumber, mp); 2382 } 2383 else 2384 { 2385 BTDefineMarkedPrism (el, edgenumber, mp); 2386 } 2387 2388 mp.matindex = el.GetIndex(); 2389 mprisms.Append (mp); 2390 break; 2391 } 2392 case PRISM: 2393 case PRISM12: 2394 { 2395 MarkedPrism mp; 2396 BTDefineMarkedPrism (el, edgenumber, mp); 2397 mp.matindex = el.GetIndex(); 2398 mprisms.Append (mp); 2399 break; 2400 } 2401 } 2402 } 2403 2404 for (i = 1; i <= nse; i++) 2405 { 2406 const Element2d & el = mesh.SurfaceElement(i); 2407 if (el.GetType() == TRIG || 2408 el.GetType() == TRIG6) 2409 { 2410 MarkedTri mt; 2411 BTDefineMarkedTri (el, edgenumber, mt); 2412 mtris.Append (mt); 2413 } 2414 else 2415 { 2416 MarkedQuad mq; 2417 BTDefineMarkedQuad (el, edgenumber, mq); 2418 mquads.Append (mq); 2419 } 2420 2421 MarkedIdentification mi; 2422 2423 2424 2425 for(j=0; j<idmaps.Size(); j++) 2426 if(BTDefineMarkedId(el, edgenumber, *idmaps[j], mi)) 2427 { 2428 mids.Append(mi); 2429 2430 int oldind = -1; 2431 for(l = 0; oldind < 0 && l<mids_old[mi.pnums[0]]->Size(); l++) 2432 { 2433 bool equal = true; 2434 for(int m = 1; equal && m < mi.np; m++) 2435 equal = (mi.pnums[m] == (*mids_old[el[0]])[l].pnums[m]); 2436 if(equal) 2437 oldind = l; 2438 } 2439 2440 if(oldind >= 0) 2441 mids.Last() = (*mids_old[mi.pnums[0]])[oldind]; 2442 } 2443 2444 } 2445 2446 2447 2448 for(int i=PointIndex::BASE; i<mesh.GetNP()+PointIndex::BASE; i++) 2449 delete mtets_old[i]; 2450 for(int i=PointIndex::BASE; i<mesh.GetNP()+PointIndex::BASE; i++) 2451 delete mprisms_old[i]; 2452 for(int i=PointIndex::BASE; i<mesh.GetNP()+PointIndex::BASE; i++) 2453 delete mids_old[i]; 2454 for(int i=PointIndex::BASE; i<mesh.GetNP()+PointIndex::BASE; i++) 2455 delete mtris_old[i]; 2456 for(int i=PointIndex::BASE; i<mesh.GetNP()+PointIndex::BASE; i++) 2457 delete mquads_old[i]; 2458 } 2459 */ 2460 2461 UpdateEdgeMarks(Mesh & mesh,const NgArray<NgArray<int,PointIndex::BASE> * > & idmaps)2462 void UpdateEdgeMarks (Mesh & mesh, 2463 const NgArray< NgArray<int,PointIndex::BASE>* > & idmaps) 2464 //const NgArray < NgArray<Element>* > & elements_before, 2465 //const NgArray < NgArray<int>* > & markedelts_num, 2466 // const NgArray < NgArray<Element2d>* > & surfelements_before, 2467 // const NgArray < NgArray<int>* > & markedsurfelts_num) 2468 { 2469 /* 2470 T_MTETS mtets_old; mtets_old.Copy(mtets); 2471 T_MPRISMS mprisms_old; mprisms_old.Copy(mprisms); 2472 T_MIDS mids_old; mids_old.Copy(mids); 2473 T_MTRIS mtris_old; mtris_old.Copy(mtris); 2474 T_MQUADS mquads_old; mquads_old.Copy(mquads); 2475 */ 2476 T_MTETS mtets_old (mtets); 2477 T_MPRISMS mprisms_old (mprisms); 2478 T_MIDS mids_old (mids); 2479 T_MTRIS mtris_old (mtris); 2480 T_MQUADS mquads_old (mquads); 2481 2482 2483 2484 mtets.SetSize(0); 2485 mprisms.SetSize(0); 2486 mids.SetSize(0); 2487 mtris.SetSize(0); 2488 mquads.SetSize(0); 2489 2490 //int nv = mesh.GetNV(); 2491 2492 2493 INDEX_2_CLOSED_HASHTABLE<int> edgenumber(9*mesh.GetNE()+4*mesh.GetNSE()); 2494 2495 int maxnum = BTSortEdges (mesh, idmaps, edgenumber); 2496 2497 for(int m = 0; m < mtets_old.Size(); m++) 2498 { 2499 MarkedTet & mt = mtets_old[m]; 2500 2501 //(*testout) << "old mt " << mt; 2502 2503 INDEX_2 edge (mt.pnums[mt.tetedge1],mt.pnums[mt.tetedge2]); 2504 edge.Sort(); 2505 if(edgenumber.Used(edge)) 2506 { 2507 int val = edgenumber.Get(edge); 2508 //(*testout) << "set voledge " << edge << " from " << val; 2509 if(val <= maxnum) 2510 { 2511 val += 2*maxnum; 2512 edgenumber.Set(edge,val); 2513 } 2514 else if(val <= 2*maxnum) 2515 { 2516 val += maxnum; 2517 edgenumber.Set(edge,val); 2518 } 2519 //(*testout) << " to " << val << endl; 2520 } 2521 2522 for(int k=0; k<4; k++) 2523 for(int i=0; i<3; i++) 2524 for(int j=i+1; i != k && j<4; j++) 2525 if(j != k && int(mt.faceedges[k]) == 6-k-i-j) 2526 { 2527 edge[0] = mt.pnums[i]; 2528 edge[1] = mt.pnums[j]; 2529 edge.Sort(); 2530 if(edgenumber.Used(edge)) 2531 { 2532 int val = edgenumber.Get(edge); 2533 //(*testout) << "set faceedge " << edge << " from " << val; 2534 if(val <= maxnum) 2535 { 2536 val += maxnum; 2537 edgenumber.Set(edge,val); 2538 } 2539 //(*testout) << " to " << val << endl; 2540 } 2541 } 2542 } 2543 2544 2545 2546 2547 for(ElementIndex ei = 0; ei < mesh.GetNE(); ei++) 2548 { 2549 const Element & el = mesh[ei]; 2550 2551 //int pos = elements_before[el[0]]->Pos(el); 2552 //int elnum = (pos >= 0) ? (*markedelts_num[el[0]])[pos] : -1; 2553 2554 switch (el.GetType()) 2555 { 2556 case TET: 2557 case TET10: 2558 { 2559 //if(elnum >= 0) 2560 // { 2561 // mtets.Append(mtets_old[elnum]); 2562 // } 2563 //else 2564 // { 2565 MarkedTet mt; 2566 BTDefineMarkedTet (el, edgenumber, mt); 2567 mt.matindex = el.GetIndex(); 2568 2569 mtets.Append (mt); 2570 2571 //(*testout) << "mtet " << mtets.Last() << endl; 2572 break; 2573 } 2574 case PYRAMID: 2575 { 2576 cerr << "Refinement :: UpdateEdgeMarks not yet implemented for pyramids" 2577 << endl; 2578 break; 2579 } 2580 2581 case PRISM: 2582 case PRISM12: 2583 { 2584 cerr << "Refinement :: UpdateEdgeMarks not yet implemented for prisms" 2585 << endl; 2586 break; 2587 } 2588 default: 2589 throw NgException("Bisect, element type not handled in switch, 6"); 2590 } 2591 2592 } 2593 2594 2595 2596 for(SurfaceElementIndex sei = 0; sei < mesh.GetNSE(); sei++) 2597 { 2598 const Element2d & el = mesh[sei]; 2599 2600 /* 2601 for(int k=0; k<3; k++) 2602 auxind3[k] = el[k]; 2603 2604 auxind3.Sort(); 2605 2606 int pos = oldfaces[auxind3[0]]->Pos(auxind3); 2607 if(pos < 0) 2608 cout << "UIUIUI" << endl; 2609 */ 2610 2611 switch (el.GetType()) 2612 { 2613 case TRIG: 2614 case TRIG6: 2615 { 2616 MarkedTri mt; 2617 BTDefineMarkedTri (el, edgenumber, mt); 2618 mtris.Append (mt); 2619 break; 2620 } 2621 2622 case QUAD: 2623 case QUAD6: 2624 { 2625 MarkedQuad mt; 2626 BTDefineMarkedQuad (el, edgenumber, mt); 2627 mquads.Append (mt); 2628 break; 2629 } 2630 default: 2631 throw NgException("Bisect, element type not handled in switch, 5"); 2632 } 2633 2634 2635 MarkedIdentification mi; 2636 for(int j=0; j<idmaps.Size(); j++) 2637 if(BTDefineMarkedId(el, edgenumber, *idmaps[j], mi)) 2638 mids.Append(mi); 2639 2640 2641 /* 2642 int pos = surfelements_before[el[0]]->Pos(el); 2643 int elnum = (pos >= 0) ? (*markedsurfelts_num[el[0]])[pos] : -1; 2644 2645 2646 switch (el.GetType()) 2647 { 2648 case TRIG: 2649 case TRIG6: 2650 { 2651 if(elnum >= 0) 2652 mtris.Append(mtris_old[elnum]); 2653 else 2654 { 2655 MarkedTri mt; 2656 BTDefineMarkedTri (el, edgenumber, mt); 2657 mtris.Append (mt); 2658 (*testout) << "(new) "; 2659 } 2660 (*testout) << "mtri " << mtris.Last(); 2661 break; 2662 } 2663 2664 case QUAD: 2665 case QUAD6: 2666 { 2667 if(elnum >= 0) 2668 mquads.Append(mquads_old[elnum]); 2669 else 2670 { 2671 MarkedQuad mt; 2672 BTDefineMarkedQuad (el, edgenumber, mt); 2673 mquads.Append (mt); 2674 } 2675 break; 2676 } 2677 } 2678 */ 2679 } 2680 2681 /* 2682 for(int i=0; i<oldfaces.Size(); i++) 2683 { 2684 delete oldfaces[i]; 2685 delete oldmarkededges[i]; 2686 } 2687 */ 2688 2689 } 2690 2691 2692 2693 Bisect(Mesh & mesh,BisectionOptions & opt,NgArray<double> * quality_loss) const2694 void Refinement :: Bisect (Mesh & mesh, 2695 BisectionOptions & opt, 2696 NgArray<double> * quality_loss) const 2697 { 2698 PrintMessage(1,"Mesh bisection"); 2699 PushStatus("Mesh bisection"); 2700 2701 static int timer = NgProfiler::CreateTimer ("Bisect"); 2702 static int timer1 = NgProfiler::CreateTimer ("Bisect 1"); 2703 static int timer1a = NgProfiler::CreateTimer ("Bisect 1a"); 2704 static int timer1b = NgProfiler::CreateTimer ("Bisect 1b"); 2705 static int timer2 = NgProfiler::CreateTimer ("Bisect 2"); 2706 static int timer2a = NgProfiler::CreateTimer ("Bisect 2a"); 2707 static int timer2b = NgProfiler::CreateTimer ("Bisect 2b"); 2708 // static int timer2c = NgProfiler::CreateTimer ("Bisect 2c"); 2709 static int timer3 = NgProfiler::CreateTimer ("Bisect 3"); 2710 static int timer3a = NgProfiler::CreateTimer ("Bisect 3a"); 2711 static int timer3b = NgProfiler::CreateTimer ("Bisect 3b"); 2712 static int timer_bisecttet = NgProfiler::CreateTimer ("Bisect tets"); 2713 static int timer_bisecttrig = NgProfiler::CreateTimer ("Bisect trigs"); 2714 static int timer_bisectsegms = NgProfiler::CreateTimer ("Bisect segms"); 2715 NgProfiler::RegionTimer reg1 (timer); 2716 2717 (*opt.tracer)("Bisect", false); 2718 2719 NgProfiler::StartTimer (timer1); 2720 NgProfiler::StartTimer (timer1a); 2721 static int localizetimer = NgProfiler::CreateTimer("localize edgepoints"); 2722 NgProfiler::RegionTimer * loct = new NgProfiler::RegionTimer(localizetimer); 2723 LocalizeEdgePoints(mesh); 2724 delete loct; 2725 2726 NgArray< NgArray<int,PointIndex::BASE>* > idmaps; 2727 for(int i=1; i<=mesh.GetIdentifications().GetMaxNr(); i++) 2728 { 2729 if(mesh.GetIdentifications().GetType(i) == Identifications::PERIODIC) 2730 { 2731 idmaps.Append(new NgArray<int,PointIndex::BASE>); 2732 mesh.GetIdentifications().GetMap(i,*idmaps.Last(),true); 2733 } 2734 } 2735 2736 2737 string refelementinfofileread = ""; 2738 string refelementinfofilewrite = ""; 2739 2740 if(opt.refinementfilename) 2741 { 2742 ifstream inf(opt.refinementfilename); 2743 string st; 2744 inf >> st; 2745 if(st == "refinementinfo") 2746 { 2747 while(inf) 2748 { 2749 while(inf && st != "markedelementsfile") 2750 inf >> st; 2751 2752 if(inf) 2753 inf >> st; 2754 2755 if(st == "read" && inf) 2756 ReadEnclString(inf,refelementinfofileread,'\"'); 2757 else if(st == "write" && inf) 2758 ReadEnclString(inf,refelementinfofilewrite,'\"'); 2759 } 2760 } 2761 inf.close(); 2762 } 2763 2764 mesh.ComputeNVertices(); 2765 2766 // if (mesh.mglevels == 1 || idmaps.Size() > 0) 2767 if (mesh.level_nv.Size() == 0) // || idmaps.Size() ???? 2768 { 2769 BisectTetsCopyMesh(mesh, NULL, opt, idmaps, refelementinfofileread); 2770 mesh.level_nv.Append (mesh.GetNV()); 2771 } 2772 2773 int np = mesh.GetNV(); 2774 mesh.SetNP(np); 2775 2776 // int ne = mesh.GetNE(); 2777 // int nse = mesh.GetNSE(); 2778 // int i, j, l; 2779 2780 // int initnp = np; 2781 // int maxsteps = 3; 2782 2783 // mesh.mglevels++; 2784 2785 /* 2786 if (opt.refinementfilename || opt.usemarkedelements) 2787 maxsteps = 3; 2788 */ 2789 2790 2791 if (opt.refine_p) 2792 { 2793 int ne = mesh.GetNE(); 2794 int nse = mesh.GetNSE(); 2795 int ox,oy,oz; 2796 for (ElementIndex ei = 0; ei < ne; ei++) 2797 if (mesh[ei].TestRefinementFlag()) 2798 { 2799 mesh[ei].GetOrder(ox,oy,oz); 2800 mesh[ei].SetOrder (ox+1,oy+1,oz+1); 2801 if (mesh[ei].TestStrongRefinementFlag()) 2802 mesh[ei].SetOrder (ox+2,oy+2,oz+2); 2803 } 2804 for (SurfaceElementIndex sei = 0; sei < nse; sei++) 2805 if (mesh[sei].TestRefinementFlag()) 2806 { 2807 mesh[sei].GetOrder(ox,oy); 2808 mesh[sei].SetOrder(ox+1,oy+1); 2809 if (mesh[sei].TestStrongRefinementFlag()) 2810 mesh[sei].SetOrder(ox+2,oy+2); 2811 } 2812 2813 #ifndef SABINE //Nachbarelemente mit ordx,ordy,ordz 2814 2815 NgArray<int,PointIndex::BASE> v_order (mesh.GetNP()); 2816 v_order = 0; 2817 2818 for (ElementIndex ei = 0; ei < ne; ei++) 2819 for (int j = 0; j < mesh[ei].GetNP(); j++) 2820 if (mesh[ei].GetOrder() > v_order[mesh[ei][j]]) 2821 v_order[mesh[ei][j]] = mesh[ei].GetOrder(); 2822 2823 for (SurfaceElementIndex sei = 0; sei < nse; sei++) 2824 for (int j = 0; j < mesh[sei].GetNP(); j++) 2825 if (mesh[sei].GetOrder() > v_order[mesh[sei][j]]) 2826 v_order[mesh[sei][j]] = mesh[sei].GetOrder(); 2827 2828 for (ElementIndex ei = 0; ei < ne; ei++) 2829 for (int j = 0; j < mesh[ei].GetNP(); j++) 2830 if (mesh[ei].GetOrder() < v_order[mesh[ei][j]]-1) 2831 mesh[ei].SetOrder(v_order[mesh[ei][j]]-1); 2832 2833 for (SurfaceElementIndex sei = 0; sei < nse; sei++) 2834 for (int j = 0; j < mesh[sei].GetNP(); j++) 2835 if (mesh[sei].GetOrder() < v_order[mesh[sei][j]]-1) 2836 mesh[sei].SetOrder(v_order[mesh[sei][j]]-1); 2837 2838 #endif 2839 2840 PopStatus(); 2841 return; 2842 } 2843 2844 2845 2846 // INDEX_2_HASHTABLE<int> cutedges(10 + 5 * (mtets.Size()+mprisms.Size()+mtris.Size()+mquads.Size())); 2847 INDEX_2_CLOSED_HASHTABLE<PointIndex> cutedges(10 + 9 * (mtets.Size()+mprisms.Size()+mtris.Size()+mquads.Size())); 2848 2849 bool noprojection = false; 2850 NgProfiler::StopTimer (timer1a); 2851 2852 for (int l = 1; l <= 1; l++) 2853 { 2854 int marked = 0; 2855 if (opt.refinementfilename) 2856 { 2857 ifstream inf(opt.refinementfilename); 2858 PrintMessage(3,"load refinementinfo from file ",opt.refinementfilename); 2859 2860 string st; 2861 inf >> st; 2862 if(st == "refinementinfo") 2863 // new version 2864 { 2865 for(int i=1; i<=mtets.Size(); i++) 2866 mtets.Elem(i).marked = 0; 2867 for(int i=1; i<=mprisms.Size(); i++) 2868 mprisms.Elem(i).marked = 0; 2869 for(int i=1; i<=mtris.Size(); i++) 2870 mtris.Elem(i).marked = 0; 2871 for(int i=1; i<=mquads.Size(); i++) 2872 mquads.Elem(i).marked = 0; 2873 for(int i=1; i<=mprisms.Size(); i++) 2874 mids.Elem(i).marked = 0; 2875 2876 inf >> st; 2877 while(inf) 2878 { 2879 if(st[0] == '#') 2880 { 2881 inf.ignore(10000,'\n'); 2882 inf >> st; 2883 } 2884 else if(st == "markedelementsfile") 2885 { 2886 inf >> st; 2887 ReadEnclString(inf,st,'\"'); 2888 inf >> st; 2889 } 2890 else if(st == "noprojection") 2891 { 2892 noprojection = true; 2893 inf >> st; 2894 } 2895 else if(st == "refine") 2896 { 2897 inf >> st; 2898 if(st == "elements") 2899 { 2900 inf >> st; 2901 bool isint = true; 2902 for(string::size_type ii=0; isint && ii<st.size(); ii++) 2903 isint = (isdigit(st[ii]) != 0); 2904 2905 while(inf && isint) 2906 { 2907 mtets.Elem(atoi(st.c_str())).marked = 3; 2908 marked = 1; 2909 2910 inf >> st; 2911 isint = true; 2912 for(string::size_type ii=0; isint && ii<st.size(); ii++) 2913 isint = (isdigit(st[ii]) != 0); 2914 } 2915 } 2916 else if(st == "orthobrick") 2917 { 2918 double bounds[6]; 2919 for(int i=0; i<6; i++) 2920 inf >> bounds[i]; 2921 2922 int cnt = 0; 2923 2924 for(ElementIndex ei = 0; ei < mesh.GetNE(); ei++) 2925 { 2926 const Element & el = mesh[ei]; 2927 2928 // 2929 Point<3> center(0,0,0); 2930 for(int i=0; i<el.GetNP(); i++) 2931 { 2932 const MeshPoint & point = mesh[el[i]]; 2933 center(0) += point(0); 2934 center(1) += point(1); 2935 center(2) += point(2); 2936 } 2937 for(int i=0; i<3; i++) 2938 center(i) *= 1./double(el.GetNP()); 2939 if(bounds[0] <= center(0) && center(0) <= bounds[3] && 2940 bounds[1] <= center(1) && center(1) <= bounds[4] && 2941 bounds[2] <= center(2) && center(2) <= bounds[5]) 2942 { 2943 mtets[ei].marked = 3; 2944 cnt++; 2945 } 2946 2947 2948 // bool contained = false; 2949 // for(int i=0; !contained && i<el.GetNP(); i++) 2950 // { 2951 // const MeshPoint & point = mesh[el[i]]; 2952 // contained = (bounds[0] <= point.X() && point.X() <= bounds[3] && 2953 // bounds[1] <= point.Y() && point.Y() <= bounds[4] && 2954 // bounds[2] <= point.Z() && point.Z() <= bounds[5]); 2955 // } 2956 // if(contained) 2957 // { 2958 // mtets[ei].marked = 3; 2959 // cnt++; 2960 // } 2961 } 2962 2963 2964 ostringstream strstr; 2965 strstr.precision(2); 2966 strstr << "marked " << float(cnt)/float(mesh.GetNE())*100. 2967 #ifdef WIN32 2968 << "%%" 2969 #else 2970 << "%" 2971 #endif 2972 <<" of the elements"; 2973 PrintMessage(4,strstr.str()); 2974 2975 if(cnt > 0) 2976 marked = 1; 2977 2978 2979 inf >> st; 2980 } 2981 else 2982 { 2983 throw NgException("something wrong with refinementinfo file"); 2984 } 2985 } 2986 } 2987 } 2988 else 2989 { 2990 inf.close(); 2991 inf.open(opt.refinementfilename); 2992 2993 char ch; 2994 for (int i = 1; i <= mtets.Size(); i++) 2995 { 2996 inf >> ch; 2997 if(!inf) 2998 throw NgException("something wrong with refinementinfo file (old format)"); 2999 mtets.Elem(i).marked = (ch == '1'); 3000 } 3001 marked = 1; 3002 } 3003 inf.close(); 3004 } 3005 3006 else if (opt.usemarkedelements) 3007 { 3008 int cntm = 0; 3009 3010 // all in one ! 3011 if (mprisms.Size()) 3012 { 3013 int cnttet = 0; 3014 int cntprism = 0; 3015 for (int i = 1; i <= mesh.GetNE(); i++) 3016 { 3017 if (mesh.VolumeElement(i).GetType() == TET || 3018 mesh.VolumeElement(i).GetType() == TET10) 3019 { 3020 cnttet++; 3021 mtets.Elem(cnttet).marked = 3022 (opt.onlyonce ? 3 : 1) * mesh.VolumeElement(i).TestRefinementFlag(); 3023 if (mtets.Elem(cnttet).marked) 3024 cntm++; 3025 } 3026 else 3027 { 3028 cntprism++; 3029 mprisms.Elem(cntprism).marked = 3030 2 * mesh.VolumeElement(i).TestRefinementFlag(); 3031 if (mprisms.Elem(cntprism).marked) 3032 cntm++; 3033 } 3034 3035 } 3036 } 3037 else 3038 for (int i = 1; i <= mtets.Size(); i++) 3039 { 3040 mtets.Elem(i).marked = 3041 (opt.onlyonce ? 1 : 3) * mesh.VolumeElement(i).TestRefinementFlag(); 3042 if (mtets.Elem(i).marked) 3043 cntm++; 3044 } 3045 3046 // (*testout) << "mtets = " << mtets << endl; 3047 3048 /* 3049 for (i = 1; i <= mtris.Size(); i++) 3050 mtris.Elem(i).marked = 0; 3051 for (i = 1; i <= mquads.Size(); i++) 3052 mquads.Elem(i).marked = 0; 3053 */ 3054 3055 if (printmessage_importance>0) 3056 { 3057 ostringstream str; 3058 str << "marked elements: " << cntm; 3059 PrintMessage(4,str.str()); 3060 } 3061 3062 int cnttrig = 0; 3063 int cntquad = 0; 3064 for (int i = 1; i <= mesh.GetNSE(); i++) 3065 { 3066 if (mesh.SurfaceElement(i).GetType() == TRIG || 3067 mesh.SurfaceElement(i).GetType() == TRIG6) 3068 { 3069 cnttrig++; 3070 mtris.Elem(cnttrig).marked = 3071 mesh.SurfaceElement(i).TestRefinementFlag() ? (opt.onlyonce ? 1 : 2) : 0; 3072 // mtris.Elem(cnttrig).marked = 0; 3073 if (mtris.Elem(cnttrig).marked) 3074 cntm++; 3075 } 3076 else 3077 { 3078 cntquad++; 3079 // 2d: marked=2, 3d prisms: marked=1 3080 mquads.Elem(cntquad).marked = 3081 mesh.SurfaceElement(i).TestRefinementFlag() ? 4-mesh.GetDimension() : 0 ; 3082 // mquads.Elem(cntquad).marked = 0; 3083 if (mquads.Elem(cntquad).marked) 3084 cntm++; 3085 } 3086 } 3087 3088 if (printmessage_importance>0) 3089 { 3090 ostringstream str; 3091 str << "with surface-elements: " << cntm; 3092 PrintMessage(4,str.str()); 3093 } 3094 3095 // he/sz: das wird oben schon richtig gemacht. 3096 // hier sind die quads vergessen! 3097 /* 3098 if (mesh.GetDimension() == 2) 3099 { 3100 cntm = 0; 3101 for (i = 1; i <= mtris.Size(); i++) 3102 { 3103 mtris.Elem(i).marked = 3104 2 * mesh.SurfaceElement(i).TestRefinementFlag(); 3105 // mtris.Elem(i).marked = 2; 3106 if (mtris.Elem(i).marked) 3107 cntm++; 3108 } 3109 3110 if (!cntm) 3111 { 3112 for (i = 1; i <= mtris.Size(); i++) 3113 { 3114 mtris.Elem(i).marked = 2; 3115 cntm++; 3116 } 3117 } 3118 cout << "trigs: " << mtris.Size() << " "; 3119 cout << "marked: " << cntm << endl; 3120 } 3121 */ 3122 3123 marked = (cntm > 0); 3124 } 3125 else 3126 { 3127 marked = BTMarkTets (mtets, mprisms, mesh); 3128 } 3129 3130 if (!marked) break; 3131 3132 3133 //(*testout) << "mtets " << mtets << endl; 3134 3135 if (opt.refine_p) 3136 { 3137 PrintMessage(3,"refine p"); 3138 3139 for (int i = 1; i <= mtets.Size(); i++) 3140 mtets.Elem(i).incorder = mtets.Elem(i).marked ? 1 : 0; 3141 3142 for (int i = 1; i <= mtets.Size(); i++) 3143 if (mtets.Elem(i).incorder) 3144 mtets.Elem(i).marked = 0; 3145 3146 3147 for (int i = 1; i <= mprisms.Size(); i++) 3148 mprisms.Elem(i).incorder = mprisms.Elem(i).marked ? 1 : 0; 3149 3150 for (int i = 1; i <= mprisms.Size(); i++) 3151 if (mprisms.Elem(i).incorder) 3152 mprisms.Elem(i).marked = 0; 3153 3154 3155 for (int i = 1; i <= mtris.Size(); i++) 3156 mtris.Elem(i).incorder = mtris.Elem(i).marked ? 1 : 0; 3157 3158 for (int i = 1; i <= mtris.Size(); i++) 3159 { 3160 if (mtris.Elem(i).incorder) 3161 mtris.Elem(i).marked = 0; 3162 } 3163 } 3164 3165 if (opt.refine_hp) 3166 { 3167 PrintMessage(3,"refine hp"); 3168 NgBitArray singv(np); 3169 singv.Clear(); 3170 3171 if (mesh.GetDimension() == 3) 3172 { 3173 for (int i = 1; i <= mesh.GetNSeg(); i++) 3174 { 3175 const Segment & seg = mesh.LineSegment(i); 3176 singv.Set (seg[0]); 3177 singv.Set (seg[1]); 3178 } 3179 /* 3180 for ( i=1; i<= mesh.GetNSE(); i++) 3181 { 3182 const Element2d & sel = mesh.SurfaceElement(i); 3183 for(int j=1; j<=sel.GetNP(); j++) 3184 singv.Set(sel.PNum(j)); 3185 } 3186 */ 3187 } 3188 else 3189 { 3190 // vertices with 2 different bnds 3191 NgArray<int> bndind(np); 3192 bndind = 0; 3193 for (int i = 1; i <= mesh.GetNSeg(); i++) 3194 { 3195 const Segment & seg = mesh.LineSegment(i); 3196 for (int j = 0; j < 2; j++) 3197 { 3198 int pi = (j == 0) ? seg[0] : seg[1]; 3199 if (bndind.Elem(pi) == 0) 3200 bndind.Elem(pi) = seg.edgenr; 3201 else if (bndind.Elem(pi) != seg.edgenr) 3202 singv.Set (pi); 3203 } 3204 } 3205 } 3206 3207 3208 3209 for (int i = 1; i <= mtets.Size(); i++) 3210 mtets.Elem(i).incorder = 1; 3211 for (int i = 1; i <= mtets.Size(); i++) 3212 { 3213 if (!mtets.Elem(i).marked) 3214 mtets.Elem(i).incorder = 0; 3215 for (int j = 0; j < 4; j++) 3216 if (singv.Test (mtets.Elem(i).pnums[j])) 3217 mtets.Elem(i).incorder = 0; 3218 } 3219 for (int i = 1; i <= mtets.Size(); i++) 3220 if (mtets.Elem(i).incorder) 3221 mtets.Elem(i).marked = 0; 3222 3223 3224 for (int i = 1; i <= mprisms.Size(); i++) 3225 mprisms.Elem(i).incorder = 1; 3226 for (int i = 1; i <= mprisms.Size(); i++) 3227 { 3228 if (!mprisms.Elem(i).marked) 3229 mprisms.Elem(i).incorder = 0; 3230 for (int j = 0; j < 6; j++) 3231 if (singv.Test (mprisms.Elem(i).pnums[j])) 3232 mprisms.Elem(i).incorder = 0; 3233 } 3234 for (int i = 1; i <= mprisms.Size(); i++) 3235 if (mprisms.Elem(i).incorder) 3236 mprisms.Elem(i).marked = 0; 3237 3238 3239 for (int i = 1; i <= mtris.Size(); i++) 3240 mtris.Elem(i).incorder = 1; 3241 for (int i = 1; i <= mtris.Size(); i++) 3242 { 3243 if (!mtris.Elem(i).marked) 3244 mtris.Elem(i).incorder = 0; 3245 for (int j = 0; j < 3; j++) 3246 if (singv.Test (mtris.Elem(i).pnums[j])) 3247 mtris.Elem(i).incorder = 0; 3248 } 3249 for (int i = 1; i <= mtris.Size(); i++) 3250 { 3251 if (mtris.Elem(i).incorder) 3252 mtris.Elem(i).marked = 0; 3253 } 3254 } 3255 3256 3257 3258 int hangingvol, hangingsurf, hangingedge; 3259 3260 //cout << "write?" << endl; 3261 //string yn; 3262 //cin >> yn; 3263 3264 3265 (*testout) << "refine volume elements" << endl; 3266 do 3267 { 3268 // refine volume elements 3269 NgProfiler::StartTimer (timer_bisecttet); 3270 (*opt.tracer)("bisecttet", false); 3271 int nel = mtets.Size(); 3272 for (int i = 1; i <= nel; i++) 3273 if (mtets.Elem(i).marked) 3274 { 3275 MarkedTet oldtet = mtets.Get(i); 3276 3277 INDEX_2 edge(oldtet.pnums[oldtet.tetedge1], 3278 oldtet.pnums[oldtet.tetedge2]); 3279 edge.Sort(); 3280 3281 PointIndex newp; 3282 if (cutedges.Used (edge)) 3283 { 3284 newp = cutedges.Get(edge); 3285 } 3286 else 3287 { 3288 Point<3> npt = Center (mesh.Point (edge.I1()), 3289 mesh.Point (edge.I2())); 3290 newp = mesh.AddPoint (npt); 3291 cutedges.Set (edge, newp); 3292 } 3293 3294 MarkedTet newtet1, newtet2; 3295 BTBisectTet (oldtet, newp, newtet1, newtet2); 3296 3297 mtets.Elem(i) = newtet1; 3298 mtets.Append (newtet2); 3299 3300 mesh.mlparentelement.Append (i); 3301 } 3302 NgProfiler::StopTimer (timer_bisecttet); 3303 (*opt.tracer)("bisecttet", true); 3304 int npr = mprisms.Size(); 3305 for (int i = 1; i <= npr; i++) 3306 if (mprisms.Elem(i).marked) 3307 { 3308 MarkedPrism oldprism; 3309 MarkedPrism newprism1, newprism2; 3310 PointIndex newp1, newp2; 3311 3312 oldprism = mprisms.Get(i); 3313 int pi1 = 0; 3314 if (pi1 == oldprism.markededge) 3315 pi1++; 3316 int pi2 = 3-pi1-oldprism.markededge; 3317 3318 INDEX_2 edge1(oldprism.pnums[pi1], 3319 oldprism.pnums[pi2]); 3320 INDEX_2 edge2(oldprism.pnums[pi1+3], 3321 oldprism.pnums[pi2+3]); 3322 edge1.Sort(); 3323 edge2.Sort(); 3324 3325 if (cutedges.Used (edge1)) 3326 newp1 = cutedges.Get(edge1); 3327 else 3328 { 3329 Point<3> npt = Center (mesh.Point (edge1.I1()), 3330 mesh.Point (edge1.I2())); 3331 newp1 = mesh.AddPoint (npt); 3332 cutedges.Set (edge1, newp1); 3333 } 3334 if (cutedges.Used (edge2)) 3335 newp2 = cutedges.Get(edge2); 3336 else 3337 { 3338 Point<3> npt = Center (mesh.Point (edge2.I1()), 3339 mesh.Point (edge2.I2())); 3340 newp2 = mesh.AddPoint (npt); 3341 cutedges.Set (edge2, newp2); 3342 } 3343 3344 3345 BTBisectPrism (oldprism, newp1, newp2, newprism1, newprism2); 3346 //if(yn == "y") 3347 // (*testout) << "bisected prism " << oldprism << "and got " << newprism1 << "and " << newprism2 << endl; 3348 mprisms.Elem(i) = newprism1; 3349 mprisms.Append (newprism2); 3350 } 3351 3352 int nid = mids.Size(); 3353 for (int i = 1; i <= nid; i++) 3354 if (mids.Elem(i).marked) 3355 { 3356 MarkedIdentification oldid,newid1,newid2; 3357 NgArray<PointIndex> newp; 3358 3359 oldid = mids.Get(i); 3360 3361 NgArray<INDEX_2> edges; 3362 edges.Append(INDEX_2(oldid.pnums[oldid.markededge], 3363 oldid.pnums[(oldid.markededge+1)%oldid.np])); 3364 edges.Append(INDEX_2(oldid.pnums[oldid.markededge + oldid.np], 3365 oldid.pnums[(oldid.markededge+1)%oldid.np + oldid.np])); 3366 3367 if(oldid.np == 4) 3368 { 3369 edges.Append(INDEX_2(oldid.pnums[(oldid.markededge+2)%oldid.np], 3370 oldid.pnums[(oldid.markededge+3)%oldid.np])); 3371 edges.Append(INDEX_2(oldid.pnums[(oldid.markededge+2)%oldid.np + oldid.np], 3372 oldid.pnums[(oldid.markededge+3)%oldid.np + oldid.np])); 3373 } 3374 for (int j = 0; j < edges.Size(); j++) 3375 { 3376 edges[j].Sort(); 3377 3378 if(cutedges.Used(edges[j])) 3379 newp.Append(cutedges.Get(edges[j])); 3380 else 3381 { 3382 Point<3> npt = Center (mesh.Point (edges[j].I1()), 3383 mesh.Point (edges[j].I2())); 3384 newp.Append(mesh.AddPoint(npt)); 3385 cutedges.Set(edges[j],newp[j]); 3386 } 3387 } 3388 3389 BTBisectIdentification(oldid,newp,newid1,newid2); 3390 mids.Elem(i) = newid1; 3391 mids.Append(newid2); 3392 } 3393 3394 3395 //IdentifyCutEdges(mesh, cutedges); 3396 3397 (*opt.tracer)("mark elements", false); 3398 3399 hangingvol = 3400 MarkHangingTets (mtets, cutedges, opt.task_manager) + 3401 MarkHangingPrisms (mprisms, cutedges) + 3402 MarkHangingIdentifications (mids, cutedges); 3403 3404 (*opt.tracer)("mark elements", true); 3405 3406 size_t nsel = mtris.Size(); 3407 NgProfiler::StartTimer (timer_bisecttrig); 3408 (*opt.tracer)("Bisect trigs", false); 3409 for (size_t i = 0; i < nsel; i++) 3410 if (mtris[i].marked) 3411 { 3412 MarkedTri newtri1, newtri2; 3413 PointIndex newp; 3414 3415 MarkedTri oldtri = mtris[i]; 3416 PointIndex oldpi1 = oldtri.pnums[(oldtri.markededge+1)%3]; 3417 PointIndex oldpi2 = oldtri.pnums[(oldtri.markededge+2)%3]; 3418 INDEX_2 edge(oldpi1, oldpi2); 3419 edge.Sort(); 3420 3421 if (cutedges.Used (edge)) 3422 { 3423 newp = cutedges.Get(edge); 3424 } 3425 else 3426 { 3427 Point<3> npt = Center (mesh.Point (edge.I1()), 3428 mesh.Point (edge.I2())); 3429 newp = mesh.AddPoint (npt); 3430 cutedges.Set (edge, newp); 3431 } 3432 3433 int si = mesh.GetFaceDescriptor (oldtri.surfid).SurfNr(); 3434 PointGeomInfo npgi; 3435 3436 if (mesh[newp].Type() != EDGEPOINT) 3437 geo.PointBetween (mesh.Point (oldpi1), mesh.Point (oldpi2), 3438 0.5, si, 3439 oldtri.pgeominfo[(oldtri.markededge+1)%3], 3440 oldtri.pgeominfo[(oldtri.markededge+2)%3], 3441 mesh.Point (newp), npgi); 3442 3443 BTBisectTri (oldtri, newp, npgi, newtri1, newtri2); 3444 3445 mtris[i] = newtri1; 3446 mtris.Append (newtri2); 3447 mesh.mlparentsurfaceelement.Append (i+1); 3448 } 3449 3450 NgProfiler::StopTimer (timer_bisecttrig); 3451 (*opt.tracer)("Bisect trigs", true); 3452 3453 int nquad = mquads.Size(); 3454 for (int i = 1; i <= nquad; i++) 3455 if (mquads.Elem(i).marked) 3456 { 3457 MarkedQuad oldquad; 3458 MarkedQuad newquad1, newquad2; 3459 PointIndex newp1, newp2; 3460 3461 oldquad = mquads.Get(i); 3462 /* 3463 INDEX_2 edge1(oldquad.pnums[0], 3464 oldquad.pnums[1]); 3465 INDEX_2 edge2(oldquad.pnums[2], 3466 oldquad.pnums[3]); 3467 */ 3468 INDEX_2 edge1, edge2; 3469 PointGeomInfo pgi11, pgi12, pgi21, pgi22; 3470 if (oldquad.markededge==0 || oldquad.markededge==2) 3471 { 3472 edge1.I1()=oldquad.pnums[0]; pgi11=oldquad.pgeominfo[0]; 3473 edge1.I2()=oldquad.pnums[1]; pgi12=oldquad.pgeominfo[1]; 3474 edge2.I1()=oldquad.pnums[2]; pgi21=oldquad.pgeominfo[2]; 3475 edge2.I2()=oldquad.pnums[3]; pgi22=oldquad.pgeominfo[3]; 3476 } 3477 else // 3 || 1 3478 { 3479 edge1.I1()=oldquad.pnums[0]; pgi11=oldquad.pgeominfo[0]; 3480 edge1.I2()=oldquad.pnums[2]; pgi12=oldquad.pgeominfo[2]; 3481 edge2.I1()=oldquad.pnums[1]; pgi21=oldquad.pgeominfo[1]; 3482 edge2.I2()=oldquad.pnums[3]; pgi22=oldquad.pgeominfo[3]; 3483 } 3484 3485 edge1.Sort(); 3486 edge2.Sort(); 3487 3488 if (cutedges.Used (edge1)) 3489 { 3490 newp1 = cutedges.Get(edge1); 3491 } 3492 else 3493 { 3494 Point<3> np1 = Center (mesh.Point (edge1.I1()), 3495 mesh.Point (edge1.I2())); 3496 newp1 = mesh.AddPoint (np1); 3497 cutedges.Set (edge1, newp1); 3498 } 3499 3500 if (cutedges.Used (edge2)) 3501 { 3502 newp2 = cutedges.Get(edge2); 3503 } 3504 else 3505 { 3506 Point<3> np2 = Center (mesh.Point (edge2.I1()), 3507 mesh.Point (edge2.I2())); 3508 newp2 = mesh.AddPoint (np2); 3509 cutedges.Set (edge2, newp2); 3510 } 3511 3512 PointGeomInfo npgi1, npgi2; 3513 3514 int si = mesh.GetFaceDescriptor (oldquad.surfid).SurfNr(); 3515 geo.PointBetween(mesh.Point (edge1.I1()), mesh.Point (edge1.I2()), 3516 0.5, si, 3517 pgi11, 3518 pgi12, 3519 mesh.Point (newp1), npgi1); 3520 geo.PointBetween (mesh.Point (edge2.I1()), mesh.Point (edge2.I2()), 3521 0.5, si, 3522 pgi21, 3523 pgi22, 3524 mesh.Point (newp2), npgi2); 3525 3526 BTBisectQuad (oldquad, newp1, npgi1, newp2, npgi2, 3527 newquad1, newquad2); 3528 3529 mquads.Elem(i) = newquad1; 3530 mquads.Append (newquad2); 3531 } 3532 3533 NgProfiler::StartTimer (timer1b); 3534 hangingsurf = 3535 MarkHangingTris (mtris, cutedges, opt.task_manager) + 3536 MarkHangingQuads (mquads, cutedges); 3537 3538 hangingedge = mesh.GetDimension() == 3 ? 0 : MarkHangingIdentifications(mids, cutedges); 3539 NgProfiler::StopTimer (timer1b); 3540 NgProfiler::StartTimer (timer_bisectsegms); 3541 int nseg = mesh.GetNSeg (); 3542 for (int i = 1; i <= nseg; i++) 3543 { 3544 Segment & seg = mesh.LineSegment (i); 3545 INDEX_2 edge(seg[0], seg[1]); 3546 edge.Sort(); 3547 if (cutedges.Used (edge)) 3548 { 3549 hangingedge = 1; 3550 Segment nseg1 = seg; 3551 Segment nseg2 = seg; 3552 3553 int newpi = cutedges.Get(edge); 3554 3555 nseg1[1] = newpi; 3556 nseg2[0] = newpi; 3557 3558 EdgePointGeomInfo newepgi; 3559 3560 geo.PointBetweenEdge(mesh.Point (seg[0]), mesh.Point (seg[1]), 3561 0.5, seg.surfnr1, seg.surfnr2, 3562 seg.epgeominfo[0], seg.epgeominfo[1], 3563 mesh.Point (newpi), newepgi); 3564 nseg1.epgeominfo[1] = newepgi; 3565 nseg2.epgeominfo[0] = newepgi; 3566 3567 mesh.LineSegment (i) = nseg1; 3568 mesh.AddSegment (nseg2); 3569 } 3570 } 3571 3572 NgProfiler::StopTimer (timer_bisectsegms); 3573 } 3574 while (hangingvol || hangingsurf || hangingedge); 3575 3576 /* 3577 if (printmessage_importance>0) 3578 { 3579 ostringstream strstr; 3580 strstr << mtets.Size() << " tets" << endl 3581 << mtris.Size() << " trigs" << endl; 3582 if (mprisms.Size()) 3583 { 3584 strstr << mprisms.Size() << " prisms" << endl 3585 << mquads.Size() << " quads" << endl; 3586 } 3587 strstr << mesh.GetNP() << " points"; 3588 PrintMessage(4,strstr.str()); 3589 } 3590 */ 3591 PrintMessage (4, mtets.Size(), " tets"); 3592 PrintMessage (4, mtris.Size(), " trigs"); 3593 if (mprisms.Size()) 3594 { 3595 PrintMessage (4, mprisms.Size(), " prisms"); 3596 PrintMessage (4, mquads.Size(), " quads"); 3597 } 3598 PrintMessage (4, mesh.GetNP(), " points"); 3599 } 3600 3601 NgProfiler::StopTimer (timer1); 3602 3603 3604 // (*testout) << "mtets = " << mtets << endl; 3605 3606 if (opt.refine_hp) 3607 { 3608 // 3609 NgArray<int> v_order (mesh.GetNP()); 3610 v_order = 0; 3611 if (mesh.GetDimension() == 3) 3612 { 3613 for (int i = 1; i <= mtets.Size(); i++) 3614 if (mtets.Elem(i).incorder) 3615 mtets.Elem(i).order++; 3616 3617 for (int i = 0; i < mtets.Size(); i++) 3618 for (int j = 0; j < 4; j++) 3619 if (int(mtets[i].order) > v_order.Elem(mtets[i].pnums[j])) 3620 v_order.Elem(mtets[i].pnums[j]) = mtets[i].order; 3621 for (int i = 0; i < mtets.Size(); i++) 3622 for (int j = 0; j < 4; j++) 3623 if (int(mtets[i].order) < v_order.Elem(mtets[i].pnums[j])-1) 3624 mtets[i].order = v_order.Elem(mtets[i].pnums[j])-1; 3625 } 3626 else 3627 { 3628 for (int i = 1; i <= mtris.Size(); i++) 3629 if (mtris.Elem(i).incorder) 3630 { 3631 mtris.Elem(i).order++; 3632 } 3633 3634 for (int i = 0; i < mtris.Size(); i++) 3635 for (int j = 0; j < 3; j++) 3636 if (int(mtris[i].order) > v_order.Elem(mtris[i].pnums[j])) 3637 v_order.Elem(mtris[i].pnums[j]) = mtris[i].order; 3638 for (int i = 0; i < mtris.Size(); i++) 3639 { 3640 for (int j = 0; j < 3; j++) 3641 if (int(mtris[i].order) < v_order.Elem(mtris[i].pnums[j])-1) 3642 mtris[i].order = v_order.Elem(mtris[i].pnums[j])-1; 3643 } 3644 } 3645 } 3646 3647 NgProfiler::StartTimer (timer2); 3648 3649 NgProfiler::StartTimer (timer2a); 3650 3651 mtets.SetAllocSize (mtets.Size()); 3652 mprisms.SetAllocSize (mprisms.Size()); 3653 mids.SetAllocSize (mids.Size()); 3654 mtris.SetAllocSize (mtris.Size()); 3655 mquads.SetAllocSize (mquads.Size()); 3656 3657 (*opt.tracer)("copy tets", false); 3658 mesh.ClearVolumeElements(); 3659 mesh.VolumeElements().SetAllocSize (mtets.Size()+mprisms.Size()); 3660 mesh.VolumeElements().SetSize(mtets.Size()); 3661 /* 3662 for (int i = 1; i <= mtets.Size(); i++) 3663 { 3664 Element el(TET); 3665 el.SetIndex (mtets.Get(i).matindex); 3666 for (int j = 1; j <= 4; j++) 3667 el.PNum(j) = mtets.Get(i).pnums[j-1]; 3668 el.SetOrder (mtets.Get(i).order); 3669 mesh.AddVolumeElement (el); 3670 } 3671 */ 3672 ParallelForRange 3673 (opt.task_manager, mtets.Size(), [&] (size_t begin, size_t end) 3674 { 3675 for (size_t i = begin; i < end; i++) 3676 { 3677 Element el(TET); 3678 auto & tet = mtets[i]; 3679 el.SetIndex (tet.matindex); 3680 el.SetOrder (tet.order); 3681 for (int j = 0; j < 4; j++) 3682 el[j] = tet.pnums[j]; 3683 mesh.SetVolumeElement (ElementIndex(i), el); 3684 } 3685 }); 3686 3687 (*opt.tracer)("copy tets", true); 3688 3689 for (int i = 1; i <= mprisms.Size(); i++) 3690 { 3691 Element el(PRISM); 3692 el.SetIndex (mprisms.Get(i).matindex); 3693 for (int j = 1; j <= 6; j++) 3694 el.PNum(j) = mprisms.Get(i).pnums[j-1]; 3695 el.SetOrder (mprisms.Get(i).order); 3696 3697 // degenerated prism ? 3698 static const int map1[] = { 3, 2, 5, 6, 1 }; 3699 static const int map2[] = { 1, 3, 6, 4, 2 }; 3700 static const int map3[] = { 2, 1, 4, 5, 3 }; 3701 3702 3703 const int * map = NULL; 3704 int deg1 = 0, deg2 = 0, deg3 = 0; 3705 // int deg = 0; 3706 if (el.PNum(1) == el.PNum(4)) { map = map1; deg1 = 1; } 3707 if (el.PNum(2) == el.PNum(5)) { map = map2; deg2 = 1; } 3708 if (el.PNum(3) == el.PNum(6)) { map = map3; deg3 = 1; } 3709 3710 switch (deg1+deg2+deg3) 3711 { 3712 case 1: 3713 { 3714 for (int j = 1; j <= 5; j++) 3715 el.PNum(j) = mprisms.Get(i).pnums[map[j-1]-1]; 3716 3717 el.SetType (PYRAMID); 3718 break; 3719 } 3720 case 2: 3721 { 3722 static const int tetmap1[] = { 1, 2, 3, 4 }; 3723 static const int tetmap2[] = { 2, 3, 1, 5 }; 3724 static const int tetmap3[] = { 3, 1, 2, 6 }; 3725 if (!deg1) map = tetmap1; 3726 if (!deg2) map = tetmap2; 3727 if (!deg3) map = tetmap3; 3728 for (int j = 1; j <= 4; j++) 3729 el.PNum(j) = mprisms.Get(i).pnums[map[j-1]-1]; 3730 /* 3731 if (!deg1) el.PNum(4) = el.PNum(4); 3732 if (!deg2) el.PNum(4) = el.PNum(5); 3733 if (!deg3) el.PNum(4) = el.PNum(6); 3734 */ 3735 el.SetType(TET); 3736 break; 3737 } 3738 default: 3739 ; 3740 } 3741 mesh.AddVolumeElement (el); 3742 } 3743 3744 mesh.ClearSurfaceElements(); 3745 mesh.SurfaceElements().SetAllocSize (mtris.Size()+mquads.Size()); 3746 NgProfiler::StopTimer (timer2a); 3747 NgProfiler::StartTimer (timer2b); 3748 /* 3749 for (auto & trig : mtris) 3750 { 3751 Element2d el(TRIG); 3752 el.SetIndex (trig.surfid); 3753 el.SetOrder (trig.order); 3754 for (int j = 0; j < 3; j++) 3755 { 3756 el[j] = trig.pnums[j]; 3757 el.GeomInfoPi(j+1) = trig.pgeominfo[j]; 3758 } 3759 mesh.AddSurfaceElement (el); 3760 } 3761 */ 3762 3763 mesh.SurfaceElements().SetSize(mtris.Size()); 3764 // for (size_t i = 0; i < mtris.Size(); i++) 3765 ParallelForRange 3766 (opt.task_manager, mtris.Size(), [&] (size_t begin, size_t end) 3767 { 3768 for (size_t i = begin; i < end; i++) 3769 { 3770 Element2d el(TRIG); 3771 auto & trig = mtris[i]; 3772 el.SetIndex (trig.surfid); 3773 el.SetOrder (trig.order); 3774 for (int j = 0; j < 3; j++) 3775 { 3776 el[j] = trig.pnums[j]; 3777 el.GeomInfoPi(j+1) = trig.pgeominfo[j]; 3778 } 3779 mesh.SetSurfaceElement (SurfaceElementIndex(i), el); 3780 } 3781 }); 3782 mesh.RebuildSurfaceElementLists(); 3783 3784 for (int i = 1; i <= mquads.Size(); i++) 3785 { 3786 Element2d el(QUAD); 3787 el.SetIndex (mquads.Get(i).surfid); 3788 for (int j = 1; j <= 4; j++) 3789 el.PNum(j) = mquads.Get(i).pnums[j-1]; 3790 Swap (el.PNum(3), el.PNum(4)); 3791 mesh.AddSurfaceElement (el); 3792 } 3793 NgProfiler::StopTimer (timer2b); 3794 3795 3796 // write multilevel hierarchy to mesh: 3797 np = mesh.GetNP(); 3798 mesh.mlbetweennodes.SetSize(np); 3799 // if (mesh.mglevels <= 2) 3800 if (mesh.level_nv.Size() <= 1) 3801 { 3802 PrintMessage(4,"RESETTING mlbetweennodes"); 3803 for (int i = 1; i <= np; i++) 3804 { 3805 mesh.mlbetweennodes.Elem(i).I1() = 0; 3806 mesh.mlbetweennodes.Elem(i).I2() = 0; 3807 } 3808 } 3809 3810 mesh.level_nv.Append (np); 3811 3812 3813 /* 3814 for (i = 1; i <= cutedges.GetNBags(); i++) 3815 for (j = 1; j <= cutedges.GetBagSize(i); j++) 3816 { 3817 INDEX_2 edge; 3818 int newpi; 3819 cutedges.GetData (i, j, edge, newpi); 3820 mesh.mlbetweennodes.Elem(newpi) = edge; 3821 } 3822 */ 3823 3824 NgBitArray isnewpoint(np); 3825 isnewpoint.Clear(); 3826 3827 for (int i = 0; i < cutedges.Size(); i++) 3828 if (cutedges.UsedPos0(i)) 3829 { 3830 INDEX_2 edge; 3831 PointIndex newpi; 3832 cutedges.GetData0 (i, edge, newpi); 3833 isnewpoint.Set(newpi); 3834 mesh.mlbetweennodes.Elem(newpi) = edge; 3835 } 3836 3837 3838 /* 3839 mesh.PrintMemInfo (cout); 3840 cout << "tets "; 3841 mtets.PrintMemInfo (cout); 3842 cout << "prims "; 3843 mprisms.PrintMemInfo (cout); 3844 cout << "tris "; 3845 mtris.PrintMemInfo (cout); 3846 cout << "quads "; 3847 mquads.PrintMemInfo (cout); 3848 cout << "cutedges "; 3849 cutedges.PrintMemInfo (cout); 3850 */ 3851 3852 3853 /* 3854 3855 // find connected nodes (close nodes) 3856 TABLE<int> conto(np); 3857 for (i = 1; i <= mprisms.Size(); i++) 3858 for (j = 1; j <= 6; j++) 3859 { 3860 int n1 = mprisms.Get(i).pnums[j-1]; 3861 int n2 = mprisms.Get(i).pnums[(j+2)%6]; 3862 // if (n1 != n2) 3863 { 3864 int found = 0; 3865 for (k = 1; k <= conto.EntrySize(n1); k++) 3866 if (conto.Get(n1, k) == n2) 3867 { 3868 found = 1; 3869 break; 3870 } 3871 if (!found) 3872 conto.Add (n1, n2); 3873 } 3874 } 3875 mesh.connectedtonode.SetSize(np); 3876 for (i = 1; i <= np; i++) 3877 mesh.connectedtonode.Elem(i) = 0; 3878 3879 3880 // (*testout) << "connection table: " << endl; 3881 // for (i = 1; i <= np; i++) 3882 // { 3883 // (*testout) << "node " << i << ": "; 3884 // for (j = 1; j <= conto.EntrySize(i); j++) 3885 // (*testout) << conto.Get(i, j) << " "; 3886 // (*testout) << endl; 3887 // } 3888 3889 3890 for (i = 1; i <= np; i++) 3891 if (mesh.connectedtonode.Elem(i) == 0) 3892 { 3893 mesh.connectedtonode.Elem(i) = i; 3894 ConnectToNodeRec (i, i, conto, mesh.connectedtonode); 3895 } 3896 */ 3897 3898 // mesh.BuildConnectedNodes(); 3899 3900 3901 3902 3903 mesh.ComputeNVertices(); 3904 3905 (*opt.tracer)("call RebuildSurfElList", false); 3906 mesh.RebuildSurfaceElementLists(); 3907 (*opt.tracer)("call RebuildSurfElList", true); 3908 3909 3910 // update identification tables 3911 for (int i = 1; i <= mesh.GetIdentifications().GetMaxNr(); i++) 3912 { 3913 NgArray<int,PointIndex::BASE> identmap; 3914 3915 mesh.GetIdentifications().GetMap (i, identmap); 3916 3917 3918 /* 3919 for (j = 1; j <= cutedges.GetNBags(); j++) 3920 for (k = 1; k <= cutedges.GetBagSize(j); k++) 3921 { 3922 INDEX_2 i2; 3923 int newpi; 3924 cutedges.GetData (j, k, i2, newpi); 3925 INDEX_2 oi2(identmap.Get(i2.I1()), 3926 identmap.Get(i2.I2())); 3927 oi2.Sort(); 3928 if (cutedges.Used (oi2)) 3929 { 3930 int onewpi = cutedges.Get(oi2); 3931 mesh.GetIdentifications().Add (newpi, onewpi, i); 3932 } 3933 } 3934 */ 3935 3936 for (int j = 0; j < cutedges.Size(); j++) 3937 if (cutedges.UsedPos0(j)) 3938 { 3939 INDEX_2 i2; 3940 PointIndex newpi; 3941 cutedges.GetData0 (j, i2, newpi); 3942 INDEX_2 oi2(identmap.Get(i2.I1()), 3943 identmap.Get(i2.I2())); 3944 oi2.Sort(); 3945 if (cutedges.Used (oi2)) 3946 { 3947 PointIndex onewpi = cutedges.Get(oi2); 3948 mesh.GetIdentifications().Add (newpi, onewpi, i); 3949 } 3950 } 3951 } 3952 3953 (*opt.tracer)("Bisect", true); 3954 3955 // Repair works only for tets! 3956 bool do_repair = mesh.PureTetMesh (); 3957 3958 do_repair = false; // JS, March 2009: multigrid crashes 3959 3960 //if(mesh.mglevels == 3) 3961 // noprojection = true; 3962 3963 //noprojection = true; 3964 3965 if(noprojection) 3966 { 3967 do_repair = false; 3968 for(int ii=1; ii<=mesh.GetNP(); ii++) 3969 { 3970 if(isnewpoint.Test(ii) && mesh.mlbetweennodes[ii][0] > 0) 3971 { 3972 mesh.Point(ii) = Center(mesh.Point(mesh.mlbetweennodes[ii][0]), 3973 mesh.Point(mesh.mlbetweennodes[ii][1])); 3974 } 3975 } 3976 } 3977 3978 3979 // Check/Repair 3980 3981 static bool repaired_once; 3982 // if(mesh.mglevels == 1) 3983 if(mesh.level_nv.Size() == 1) 3984 repaired_once = false; 3985 3986 //mesh.Save("before.vol"); 3987 3988 static int reptimer = NgProfiler::CreateTimer("check/repair"); 3989 NgProfiler::RegionTimer * regt(NULL); 3990 regt = new NgProfiler::RegionTimer(reptimer); 3991 3992 NgArray<ElementIndex> bad_elts; 3993 NgArray<double> pure_badness; 3994 3995 if(do_repair || quality_loss != NULL) 3996 { 3997 pure_badness.SetSize(mesh.GetNP()+2); 3998 GetPureBadness(mesh,pure_badness,isnewpoint); 3999 } 4000 4001 4002 if(do_repair) // by Markus W 4003 { 4004 const double max_worsening = 1; 4005 4006 const bool uselocalworsening = false; 4007 4008 // bool repaired = false; 4009 4010 Validate(mesh,bad_elts,pure_badness,max_worsening,uselocalworsening); 4011 4012 if (printmessage_importance>0) 4013 { 4014 ostringstream strstr; 4015 for(int ii=0; ii<bad_elts.Size(); ii++) 4016 strstr << "bad element " << bad_elts[ii] << "\n"; 4017 PrintMessage(1,strstr.str()); 4018 } 4019 if(repaired_once || bad_elts.Size() > 0) 4020 { 4021 clock_t t1(clock()); 4022 4023 4024 // update id-maps 4025 int j=0; 4026 for(int i=1; i<=mesh.GetIdentifications().GetMaxNr(); i++) 4027 { 4028 if(mesh.GetIdentifications().GetType(i) == Identifications::PERIODIC) 4029 { 4030 mesh.GetIdentifications().GetMap(i,*idmaps[j],true); 4031 j++; 4032 } 4033 } 4034 4035 4036 // do the repair 4037 try 4038 { 4039 RepairBisection(mesh,bad_elts,isnewpoint,*this, 4040 pure_badness, 4041 max_worsening,uselocalworsening, 4042 idmaps); 4043 // repaired = true; 4044 repaired_once = true; 4045 } 4046 catch(NgException & ex) 4047 { 4048 PrintMessage(1,string("Problem: ") + ex.What()); 4049 } 4050 4051 4052 if (printmessage_importance>0) 4053 { 4054 ostringstream strstr; 4055 strstr << "Time for Repair: " << double(clock() - t1)/double(CLOCKS_PER_SEC) << endl 4056 << "bad elements after repair: " << bad_elts << endl; 4057 PrintMessage(1,strstr.str()); 4058 } 4059 4060 if(quality_loss != NULL) 4061 Validate(mesh,bad_elts,pure_badness,1e100,uselocalworsening,quality_loss); 4062 4063 if(idmaps.Size() == 0) 4064 UpdateEdgeMarks(mesh,idmaps); 4065 4066 /* 4067 if(1==1) 4068 UpdateEdgeMarks(mesh,idmaps); 4069 else 4070 mesh.mglevels = 1; 4071 */ 4072 4073 //mesh.ImproveMesh(); 4074 4075 } 4076 } 4077 delete regt; 4078 4079 4080 4081 for(int i=0; i<idmaps.Size(); i++) 4082 delete idmaps[i]; 4083 idmaps.DeleteAll(); 4084 4085 NgProfiler::StopTimer (timer2); 4086 NgProfiler::StartTimer (timer3); 4087 4088 NgProfiler::StartTimer (timer3a); 4089 (*opt.tracer)("topology from bisect", false); 4090 mesh.UpdateTopology(opt.task_manager, opt.tracer); 4091 (*opt.tracer)("topology from bisect", true); 4092 NgProfiler::StopTimer (timer3a); 4093 4094 4095 4096 if(refelementinfofilewrite != "") 4097 { 4098 PrintMessage(3,"writing marked-elements information to \"",refelementinfofilewrite,"\""); 4099 ofstream ofst(refelementinfofilewrite.c_str()); 4100 4101 WriteMarkedElements(ofst); 4102 4103 ofst.close(); 4104 } 4105 4106 NgProfiler::StartTimer (timer3b); 4107 mesh.CalcSurfacesOfNode(); 4108 NgProfiler::StopTimer (timer3b); 4109 4110 PrintMessage (1, "Bisection done"); 4111 4112 PopStatus(); 4113 NgProfiler::StopTimer (timer3); 4114 } 4115 4116 4117 4118 BisectionOptions()4119 BisectionOptions :: BisectionOptions () 4120 { 4121 outfilename = NULL; 4122 mlfilename = NULL; 4123 refinementfilename = NULL; 4124 femcode = NULL; 4125 maxlevel = 50; 4126 usemarkedelements = 0; 4127 refine_hp = 0; 4128 refine_p = 0; 4129 } 4130 } 4131