1 #include <mystdlib.h> 2 #include <atomic> 3 #include <set> 4 #include "meshing.hpp" 5 6 #ifdef NG_PYTHON 7 // must be included to instantiate Archive::Shallow(NetgenGeometry&) 8 #include <core/python_ngcore.hpp> 9 #endif 10 11 namespace netgen 12 { 13 Find3dElement(const Mesh & mesh,const netgen::Point<3> & p,double * lami,const NgArray<int> * const indices,BoxTree<3> * searchtree,const bool allowindex=true)14 int Find3dElement (const Mesh& mesh, 15 const netgen::Point<3> & p, 16 double * lami, 17 const NgArray<int> * const indices, 18 BoxTree<3> * searchtree, 19 const bool allowindex = true) 20 { 21 int ne = 0; 22 NgArray<int> locels; 23 if (searchtree) 24 { 25 searchtree->GetIntersecting (p, p, locels); 26 ne = locels.Size(); 27 } 28 else 29 ne = mesh.GetNE(); 30 31 for (int i = 1; i <= ne; i++) 32 { 33 int ii; 34 35 if (searchtree) 36 ii = locels.Get(i); 37 else 38 ii = i; 39 40 if(indices != NULL && indices->Size() > 0) 41 { 42 bool contained = indices->Contains(mesh.VolumeElement(ii).GetIndex()); 43 if((allowindex && !contained) || (!allowindex && contained)) continue; 44 } 45 46 if(mesh.PointContainedIn3DElement(p,lami,ii)) 47 return ii; 48 } 49 50 // Not found, try uncurved variant: 51 for (int i = 1; i <= ne; i++) 52 { 53 int ii; 54 55 if (searchtree) 56 ii = locels.Get(i); 57 else 58 ii = i; 59 60 if(indices != NULL && indices->Size() > 0) 61 { 62 bool contained = indices->Contains(mesh.VolumeElement(ii).GetIndex()); 63 if((allowindex && !contained) || (!allowindex && contained)) continue; 64 } 65 66 67 if(mesh.PointContainedIn3DElementOld(p,lami,ii)) 68 { 69 (*testout) << "WARNING: found element of point " << p <<" only for uncurved mesh" << endl; 70 return ii; 71 } 72 } 73 return 0; 74 } 75 Find2dElement(const Mesh & mesh,const netgen::Point<3> & p,double * lami,const NgArray<int> * const indices,BoxTree<3> * searchtree,const bool allowindex=true)76 int Find2dElement (const Mesh& mesh, 77 const netgen::Point<3> & p, 78 double * lami, 79 const NgArray<int> * const indices, 80 BoxTree<3> * searchtree, 81 const bool allowindex = true) 82 { 83 double vlam[3]; 84 int velement = 0; 85 86 if(mesh.GetNE()) 87 velement = Find3dElement(mesh, p,vlam,NULL,searchtree,allowindex); 88 89 //(*testout) << "p " << p << endl; 90 //(*testout) << "velement " << velement << endl; 91 92 // first try to find a volume element containing p and project to face 93 if(velement!=0) 94 { 95 auto & topology = mesh.GetTopology(); 96 NgArray<int> faces; 97 topology.GetElementFaces(velement,faces); 98 99 //(*testout) << "faces " << faces << endl; 100 101 for(int i=0; i<faces.Size(); i++) 102 faces[i] = topology.GetFace2SurfaceElement(faces[i]); 103 104 //(*testout) << "surfel " << faces << endl; 105 106 for(int i=0; i<faces.Size(); i++) 107 { 108 if(faces[i] == 0) 109 continue; 110 auto sel = mesh.SurfaceElement(faces[i]); 111 if(indices && indices->Size() != 0 && !indices->Contains(sel.GetIndex())) 112 continue; 113 114 auto & el = mesh.VolumeElement(velement); 115 116 if (el.GetType() == TET) 117 { 118 double lam4[4] = { vlam[0], vlam[1], vlam[2], 1.0-vlam[0]-vlam[1]-vlam[2] }; 119 double face_lam = lam4[i]; 120 if(face_lam < 1e-5) 121 { 122 // found volume point very close to a face -> use barycentric coordinates directly 123 lami[2] = 0.0; 124 for(auto j : Range(1,3)) 125 for(auto k : Range(4)) 126 if(sel[j] == el[k]) 127 lami[j-1] = lam4[k]/(1.0-face_lam); 128 return faces[i]; 129 } 130 } 131 132 if(mesh.PointContainedIn2DElement(p,lami,faces[i],true)) 133 return faces[i]; 134 } 135 } 136 137 // Did't find any matching face of a volume element, search 2d elements directly 138 int ne; 139 140 NgArray<int> locels; 141 // TODO: build search tree for surface elements 142 if (!mesh.GetNE() && searchtree) 143 { 144 searchtree->GetIntersecting (p, p, locels); 145 ne = locels.Size(); 146 } 147 else 148 ne = mesh.GetNSE(); 149 150 for (int i = 1; i <= ne; i++) 151 { 152 int ii; 153 154 if (locels.Size()) 155 ii = locels.Get(i); 156 else 157 ii = i; 158 159 if(indices != NULL && indices->Size() > 0) 160 { 161 bool contained = indices->Contains(mesh.SurfaceElement(ii).GetIndex()); 162 if((allowindex && !contained) || (!allowindex && contained)) continue; 163 } 164 165 if(mesh.PointContainedIn2DElement(p,lami,ii)) return ii; 166 167 } 168 return 0; 169 } 170 Find1dElement(const Mesh & mesh,const netgen::Point<3> & p,double * lami,const NgArray<int> * const indices,BoxTree<3> * searchtree,const bool allowindex=true)171 int Find1dElement (const Mesh& mesh, 172 const netgen::Point<3> & p, 173 double * lami, 174 const NgArray<int> * const indices, 175 BoxTree<3> * searchtree, 176 const bool allowindex = true) 177 { 178 double vlam[3]; 179 int velement = Find2dElement(mesh, p, vlam, NULL, searchtree, allowindex); 180 if(velement == 0) 181 return 0; 182 183 vlam[2] = 1.-vlam[0] - vlam[1]; 184 NgArray<int> edges; 185 auto & topology = mesh.GetTopology(); 186 topology.GetSurfaceElementEdges(velement, edges); 187 Array<SegmentIndex> segs(edges.Size()); 188 for(auto i : Range(edges)) 189 segs[i] = topology.GetSegmentOfEdge(edges[i]); 190 191 for(auto i : Range(segs)) 192 { 193 if(IsInvalid(segs[i])) 194 continue; 195 auto& el = mesh.SurfaceElement(velement); 196 if(el.GetType() == TRIG) 197 { 198 double seg_lam; 199 double lam; 200 auto seg = mesh.LineSegment(segs[i]); 201 for(auto k : Range(3)) 202 { 203 if(seg[0] == el[k]) 204 lam = vlam[k]; 205 if(seg[1] == el[k]) 206 seg_lam = vlam[k]; 207 } 208 if(1.- seg_lam - lam < 1e-5) 209 { 210 // found point close to segment -> use barycentric coordinates directly 211 lami[0] = lam; 212 return int(segs[i])+1; 213 } 214 } 215 else 216 throw NgException("Quad not implemented yet!"); 217 } 218 219 return 0; 220 } 221 222 static mutex buildsearchtree_mutex; 223 Mesh()224 Mesh :: Mesh () 225 : topology(*this), surfarea(*this) 226 { 227 boundaryedges = nullptr; 228 surfelementht = nullptr; 229 segmentht = nullptr; 230 231 lochfunc = nullptr; 232 // mglevels = 1; 233 elementsearchtree = nullptr; 234 elementsearchtreets = NextTimeStamp(); 235 majortimestamp = timestamp = NextTimeStamp(); 236 hglob = 1e10; 237 hmin = 0; 238 numvertices = -1; 239 dimension = 3; 240 241 curvedelems = make_unique<CurvedElements> (*this); 242 clusters = make_unique<AnisotropicClusters> (*this); 243 ident = make_unique<Identifications> (*this); 244 245 hpelements = NULL; 246 coarsemesh = NULL; 247 248 ps_startelement = 0; 249 250 geomtype = NO_GEOM; 251 252 bcnames.SetSize(0); 253 cd2names.SetSize(0); 254 255 // this->comm = netgen :: ng_comm; 256 #ifdef PARALLEL 257 paralleltop = make_unique<ParallelMeshTopology> (*this); 258 #endif 259 } 260 261 ~Mesh()262 Mesh :: ~Mesh() 263 { 264 // delete lochfunc; 265 // delete boundaryedges; 266 // delete surfelementht; 267 // delete segmentht; 268 // delete curvedelems; 269 // delete clusters; 270 // delete ident; 271 // delete elementsearchtree; 272 // delete coarsemesh; 273 // delete hpelements; 274 275 for (int i = 0; i < materials.Size(); i++) 276 delete materials[i]; 277 for(int i = 0; i < userdata_int.Size(); i++) 278 delete userdata_int[i]; 279 for(int i = 0; i < userdata_double.Size(); i++) 280 delete userdata_double[i]; 281 282 for (int i = 0; i < bcnames.Size(); i++ ) 283 delete bcnames[i]; 284 285 for (int i = 0; i < cd2names.Size(); i++) 286 delete cd2names[i]; 287 288 // #ifdef PARALLEL 289 // delete paralleltop; 290 // #endif 291 } 292 SetCommunicator(NgMPI_Comm acomm)293 void Mesh :: SetCommunicator(NgMPI_Comm acomm) 294 { 295 this->comm = acomm; 296 } 297 operator =(const Mesh & mesh2)298 Mesh & Mesh :: operator= (const Mesh & mesh2) 299 { 300 geometry = mesh2.geometry; 301 dimension = mesh2.dimension; 302 points = mesh2.points; 303 segments = mesh2.segments; 304 surfelements = mesh2.surfelements; 305 volelements = mesh2.volelements; 306 lockedpoints = mesh2.lockedpoints; 307 facedecoding = mesh2.facedecoding; 308 dimension = mesh2.dimension; 309 310 311 materials.SetSize( mesh2.materials.Size() ); 312 for ( int i = 0; i < mesh2.materials.Size(); i++ ) 313 if ( mesh2.materials[i] ) materials[i] = new string ( *mesh2.materials[i] ); 314 else materials[i] = 0; 315 316 std::map<const string*, string*> bcmap; 317 bcnames.SetSize( mesh2.bcnames.Size() ); 318 for ( int i = 0; i < mesh2.bcnames.Size(); i++ ) 319 { 320 if ( mesh2.bcnames[i] ) bcnames[i] = new string ( *mesh2.bcnames[i] ); 321 else bcnames[i] = 0; 322 bcmap[mesh2.bcnames[i]] = bcnames[i]; 323 } 324 325 // Remap string* members in FaceDescriptor to new mesh 326 for (auto & f : facedecoding) 327 f.SetBCName( bcmap[&f.GetBCName()] ); 328 329 330 cd2names.SetSize(mesh2.cd2names.Size()); 331 for (int i=0; i < mesh2.cd2names.Size(); i++) 332 if (mesh2.cd2names[i]) cd2names[i] = new string(*mesh2.cd2names[i]); 333 else cd2names[i] = 0; 334 335 cd3names.SetSize(mesh2.cd3names.Size()); 336 for (int i=0; i < mesh2.cd3names.Size(); i++) 337 if (mesh2.cd3names[i]) cd3names[i] = new string(*mesh2.cd3names[i]); 338 else cd3names[i] = 0; 339 340 numvertices = mesh2.numvertices; 341 342 return *this; 343 } 344 345 DeleteMesh()346 void Mesh :: DeleteMesh() 347 { 348 NgLock lock(mutex); 349 lock.Lock(); 350 points.SetSize(0); 351 segments.SetSize(0); 352 surfelements.SetSize(0); 353 volelements.SetSize(0); 354 lockedpoints.SetSize(0); 355 // surfacesonnode.SetSize(0); 356 357 // delete boundaryedges; 358 boundaryedges = nullptr; 359 segmentht = nullptr; 360 surfelementht = nullptr; 361 362 openelements.SetSize(0); 363 facedecoding.SetSize(0); 364 365 ident = make_unique<Identifications> (*this); 366 topology = MeshTopology (*this); 367 curvedelems = make_unique<CurvedElements> (*this); 368 clusters = make_unique<AnisotropicClusters> (*this); 369 370 for ( int i = 0; i < bcnames.Size(); i++ ) 371 if ( bcnames[i] ) delete bcnames[i]; 372 for (int i= 0; i< cd2names.Size(); i++) 373 if (cd2names[i]) delete cd2names[i]; 374 375 #ifdef PARALLEL 376 paralleltop = make_unique<ParallelMeshTopology> (*this); 377 #endif 378 379 lock.UnLock(); 380 381 timestamp = NextTimeStamp(); 382 } 383 384 ClearSurfaceElements()385 void Mesh :: ClearSurfaceElements() 386 { 387 surfelements.SetSize(0); 388 /* 389 for (int i = 0; i < facedecoding.Size(); i++) 390 facedecoding[i].firstelement = -1; 391 */ 392 for (auto & fd : facedecoding) 393 fd.firstelement = -1; 394 395 timestamp = NextTimeStamp(); 396 } 397 398 399 AddPoint(const Point3d & p,int layer)400 PointIndex Mesh :: AddPoint (const Point3d & p, int layer) 401 { 402 return AddPoint (p, layer, INNERPOINT); 403 } 404 AddPoint(const Point3d & p,int layer,POINTTYPE type)405 PointIndex Mesh :: AddPoint (const Point3d & p, int layer, POINTTYPE type) 406 { 407 408 // PointIndex pi = points.End(); 409 PointIndex pi = *points.Range().end(); 410 if (points.Size() == points.AllocSize()) 411 { 412 NgLock lock(mutex); 413 lock.Lock(); 414 points.Append ( MeshPoint (p, layer, type) ); 415 lock.UnLock(); 416 } 417 else 418 { 419 points.Append ( MeshPoint (p, layer, type) ); 420 } 421 422 timestamp = NextTimeStamp(); 423 424 return pi; 425 } 426 427 AddSegment(const Segment & s)428 SegmentIndex Mesh :: AddSegment (const Segment & s) 429 { 430 NgLock lock(mutex); 431 lock.Lock(); 432 timestamp = NextTimeStamp(); 433 434 int maxn = max2 (s[0], s[1]); 435 maxn += 1-PointIndex::BASE; 436 437 /* 438 if (maxn > ptyps.Size()) 439 { 440 int maxo = ptyps.Size(); 441 ptyps.SetSize (maxn); 442 for (int i = maxo; i < maxn; i++) 443 ptyps[i] = INNERPOINT; 444 } 445 446 if (ptyps[s[0]] > EDGEPOINT) ptyps[s[0]] = EDGEPOINT; 447 if (ptyps[s[1]] > EDGEPOINT) ptyps[s[1]] = EDGEPOINT; 448 */ 449 450 if (maxn <= points.Size()) 451 { 452 if (points[s[0]].Type() > EDGEPOINT) 453 points[s[0]].SetType (EDGEPOINT); 454 if (points[s[1]].Type() > EDGEPOINT) 455 points[s[1]].SetType (EDGEPOINT); 456 } 457 /* 458 else 459 { 460 cerr << "edge points nrs > points.Size" << endl; 461 } 462 */ 463 464 SegmentIndex si = segments.Size(); 465 segments.Append (s); 466 467 lock.UnLock(); 468 return si; 469 } 470 AddSurfaceElement(const Element2d & el)471 SurfaceElementIndex Mesh :: AddSurfaceElement (const Element2d & el) 472 { 473 timestamp = NextTimeStamp(); 474 475 PointIndex maxn = el[0]; 476 for (int i = 1; i < el.GetNP(); i++) 477 if (el[i] > maxn) maxn = el[i]; 478 479 /* 480 maxn += 1-PointIndex::BASE; 481 if (maxn <= points.Size()) 482 { 483 for (int i = 0; i < el.GetNP(); i++) 484 if (points[el[i]].Type() > SURFACEPOINT) 485 points[el[i]].SetType(SURFACEPOINT); 486 } 487 */ 488 // if (maxn < points.End()) 489 if (maxn < *points.Range().end()) 490 for (PointIndex pi : el.PNums()) 491 if (points[pi].Type() > SURFACEPOINT) 492 points[pi].SetType(SURFACEPOINT); 493 494 495 SurfaceElementIndex si = surfelements.Size(); 496 if (surfelements.AllocSize() == surfelements.Size()) 497 { 498 NgLock lock(mutex); 499 lock.Lock(); 500 surfelements.Append (el); 501 lock.UnLock(); 502 } 503 else 504 { 505 surfelements.Append (el); 506 } 507 508 if (el.index<=0 || el.index > facedecoding.Size()) 509 cerr << "has no facedecoding: fd.size = " << facedecoding.Size() << ", ind = " << el.index << endl; 510 511 surfelements.Last().next = facedecoding[el.index-1].firstelement; 512 facedecoding[el.index-1].firstelement = si; 513 514 if (SurfaceArea().Valid()) 515 SurfaceArea().Add (el); 516 517 return si; 518 } 519 SetSurfaceElement(SurfaceElementIndex sei,const Element2d & el)520 void Mesh :: SetSurfaceElement (SurfaceElementIndex sei, const Element2d & el) 521 { 522 int maxn = el[0]; 523 for (int i = 1; i < el.GetNP(); i++) 524 if (el[i] > maxn) maxn = el[i]; 525 526 maxn += 1-PointIndex::BASE; 527 528 if (maxn <= points.Size()) 529 { 530 for (int i = 0; i < el.GetNP(); i++) 531 if (points[el[i]].Type() > SURFACEPOINT) 532 points[el[i]].SetType(SURFACEPOINT); 533 } 534 535 surfelements[sei] = el; 536 if (el.index > facedecoding.Size()) 537 cerr << "has no facedecoding: fd.size = " << facedecoding.Size() << ", ind = " << el.index << endl; 538 539 // add lock-free to list ... slow, call RebuildSurfaceElementLists later 540 /* 541 surfelements[sei].next = facedecoding[el.index-1].firstelement; 542 auto & head = reinterpret_cast<atomic<SurfaceElementIndex>&> (facedecoding[el.index-1].firstelement); 543 while (!head.compare_exchange_weak (surfelements[sei].next, sei)) 544 ; 545 */ 546 547 /* 548 if (SurfaceArea().Valid()) 549 SurfaceArea().Add (el); 550 */ 551 } 552 553 AddVolumeElement(const Element & el)554 ElementIndex Mesh :: AddVolumeElement (const Element & el) 555 { 556 /* 557 int maxn = el[0]; 558 for (int i = 1; i < el.GetNP(); i++) 559 if (el[i] > maxn) maxn = el[i]; 560 561 maxn += 1-PointIndex::BASE; 562 */ 563 564 /* 565 if (maxn > ptyps.Size()) 566 { 567 int maxo = ptyps.Size(); 568 ptyps.SetSize (maxn); 569 for (i = maxo+PointIndex::BASE; 570 i < maxn+PointIndex::BASE; i++) 571 ptyps[i] = INNERPOINT; 572 } 573 */ 574 /* 575 if (maxn > points.Size()) 576 { 577 cerr << "add vol element before point" << endl; 578 } 579 */ 580 581 int ve = volelements.Size(); 582 583 if (volelements.Size() == volelements.AllocSize()) 584 { 585 NgLock lock(mutex); 586 lock.Lock(); 587 volelements.Append (el); 588 lock.UnLock(); 589 } 590 else 591 { 592 volelements.Append (el); 593 } 594 volelements.Last().flags.illegal_valid = 0; 595 volelements.Last().flags.fixed = 0; 596 volelements.Last().flags.deleted = 0; 597 598 // while (volelements.Size() > eltyps.Size()) 599 // eltyps.Append (FREEELEMENT); 600 601 timestamp = NextTimeStamp(); 602 603 return ve; 604 } 605 SetVolumeElement(ElementIndex ei,const Element & el)606 void Mesh :: SetVolumeElement (ElementIndex ei, const Element & el) 607 { 608 /* 609 int maxn = el[0]; 610 for (int i = 1; i < el.GetNP(); i++) 611 if (el[i] > maxn) maxn = el[i]; 612 613 maxn += 1-PointIndex::BASE; 614 */ 615 616 volelements[ei] = el; 617 volelements[ei].flags.illegal_valid = 0; 618 volelements[ei].flags.fixed = 0; 619 volelements[ei].flags.deleted = 0; 620 } 621 622 623 624 625 Save(const string & filename) const626 void Mesh :: Save (const string & filename) const 627 { 628 if (filename.find(".vol.bin") != string::npos) 629 { 630 BinaryOutArchive in(filename); 631 in & const_cast<Mesh&>(*this); 632 return; 633 } 634 635 ostream * outfile; 636 if (filename.find(".vol.gz")!=string::npos) 637 outfile = new ogzstream(filename.c_str()); 638 else if (filename.find(".vol")!=string::npos) 639 outfile = new ofstream(filename.c_str()); 640 else 641 outfile = new ogzstream((filename+".vol.gz").c_str()); 642 643 Save(*outfile); 644 delete outfile; 645 } 646 647 648 Save(ostream & outfile) const649 void Mesh :: Save (ostream & outfile) const 650 { 651 static Timer timer("Mesh::Save"); RegionTimer rt(timer); 652 int i, j; 653 654 double scale = 1; // globflags.GetNumFlag ("scale", 1); 655 int inverttets = 0; // globflags.GetDefineFlag ("inverttets"); 656 int invertsurf = 0; // globflags.GetDefineFlag ("invertsurfacemesh"); 657 658 outfile << "# Generated by NETGEN " << GetLibraryVersion("netgen") << endl << endl; 659 660 661 outfile << "mesh3d" << "\n"; 662 663 outfile << "dimension\n" << GetDimension() << "\n"; 664 665 outfile << "geomtype\n" << int(geomtype) << "\n"; 666 667 668 outfile << "\n"; 669 outfile << "# surfnr bcnr domin domout np p1 p2 p3" 670 << "\n"; 671 672 673 switch (geomtype) 674 { 675 case GEOM_STL: 676 outfile << "surfaceelementsgi" << "\n"; 677 break; 678 case GEOM_OCC: case GEOM_ACIS: 679 outfile << "surfaceelementsuv" << "\n"; 680 break; 681 default: 682 outfile << "surfaceelements" << "\n"; 683 } 684 685 outfile << GetNSE() << "\n"; 686 687 for (SurfaceElementIndex sei = 0; sei < GetNSE(); sei++) 688 { 689 if ((*this)[sei].GetIndex()) 690 { 691 outfile << " " << GetFaceDescriptor((*this)[sei].GetIndex ()).SurfNr()+1; 692 outfile << " " << GetFaceDescriptor((*this)[sei].GetIndex ()).BCProperty(); 693 outfile << " " << GetFaceDescriptor((*this)[sei].GetIndex ()).DomainIn(); 694 outfile << " " << GetFaceDescriptor((*this)[sei].GetIndex ()).DomainOut(); 695 } 696 else 697 outfile << " 0 0 0"; 698 699 Element2d sel = (*this)[sei]; 700 if (invertsurf) 701 sel.Invert(); 702 703 outfile << " " << sel.GetNP(); 704 for (j = 0; j < sel.GetNP(); j++) 705 outfile << " " << sel[j]; 706 707 switch (geomtype) 708 { 709 case GEOM_STL: 710 for (j = 1; j <= sel.GetNP(); j++) 711 outfile << " " << sel.GeomInfoPi(j).trignum; 712 break; 713 case GEOM_OCC: case GEOM_ACIS: 714 for (j = 1; j <= sel.GetNP(); j++) 715 { 716 outfile << " " << sel.GeomInfoPi(j).u; 717 outfile << " " << sel.GeomInfoPi(j).v; 718 } 719 break; 720 default: 721 ; 722 } 723 outfile << "\n"; 724 } 725 726 outfile << "\n" << "\n"; 727 outfile << "# matnr np p1 p2 p3 p4" << "\n"; 728 outfile << "volumeelements" << "\n"; 729 outfile << GetNE() << "\n"; 730 731 for (ElementIndex ei = 0; ei < GetNE(); ei++) 732 { 733 outfile << (*this)[ei].GetIndex(); 734 outfile << " " << (*this)[ei].GetNP(); 735 736 Element el = (*this)[ei]; 737 if (inverttets) el.Invert(); 738 739 for (j = 0; j < el.GetNP(); j++) 740 outfile << " " << el[j]; 741 outfile << "\n"; 742 } 743 744 745 outfile << "\n" << "\n"; 746 // outfile << " surf1 surf2 p1 p2" << "\n"; 747 outfile << "# surfid 0 p1 p2 trignum1 trignum2 domin/surfnr1 domout/surfnr2 ednr1 dist1 ednr2 dist2 \n"; 748 outfile << "edgesegmentsgi2" << "\n"; 749 outfile << GetNSeg() << "\n"; 750 751 for (i = 1; i <= GetNSeg(); i++) 752 { 753 const Segment & seg = LineSegment (i); 754 outfile.width(8); 755 outfile << seg.si; // 2D: bc number, 3D: wievielte Kante 756 outfile.width(8); 757 outfile << 0; 758 outfile.width(8); 759 outfile << seg[0]; 760 outfile.width(8); 761 outfile << seg[1]; 762 outfile << " "; 763 outfile.width(8); 764 outfile << seg.geominfo[0].trignum; // stl dreiecke 765 outfile << " "; 766 outfile.width(8); 767 outfile << seg.geominfo[1].trignum; // << endl; // stl dreieck 768 769 if (dimension == 3) 770 { 771 outfile << " "; 772 outfile.width(8); 773 outfile << seg.surfnr1+1; 774 outfile << " "; 775 outfile.width(8); 776 outfile << seg.surfnr2+1; 777 } 778 else 779 { 780 outfile << " "; 781 outfile.width(8); 782 outfile << seg.domin; 783 outfile << " "; 784 outfile.width(8); 785 outfile << seg.domout; 786 } 787 788 outfile << " "; 789 outfile.width(8); 790 outfile << seg.edgenr; 791 outfile << " "; 792 outfile.width(12); 793 outfile.precision(16); 794 outfile << seg.epgeominfo[0].dist; // splineparameter (2D) 795 outfile << " "; 796 outfile.width(8); 797 outfile.precision(16); 798 outfile << seg.epgeominfo[1].edgenr; // geometry dependent 799 outfile << " "; 800 outfile.width(12); 801 outfile << seg.epgeominfo[1].dist; 802 803 outfile << "\n"; 804 } 805 806 807 outfile << "\n" << "\n"; 808 outfile << "# X Y Z" << "\n"; 809 outfile << "points" << "\n"; 810 outfile << GetNP() << "\n"; 811 outfile.precision(16); 812 outfile.setf (ios::fixed, ios::floatfield); 813 outfile.setf (ios::showpoint); 814 815 PointIndex pi; 816 for (pi = PointIndex::BASE; 817 pi < GetNP()+PointIndex::BASE; pi++) 818 { 819 outfile.width(22); 820 outfile << (*this)[pi](0)/scale << " "; 821 outfile.width(22); 822 outfile << (*this)[pi](1)/scale << " "; 823 outfile.width(22); 824 outfile << (*this)[pi](2)/scale << "\n"; 825 } 826 827 outfile << "\n" << "\n"; 828 outfile << "# pnum index" << "\n"; 829 outfile << "pointelements" << "\n"; 830 outfile << pointelements.Size() << "\n"; 831 832 for (i = 0; i < pointelements.Size(); i++) 833 { 834 outfile.width(8); 835 outfile << pointelements[i].pnum << " "; 836 outfile.width(8); 837 outfile << pointelements[i].index << "\n"; 838 } 839 840 if (ident -> GetMaxNr() > 0) 841 { 842 outfile << "identifications\n"; 843 NgArray<INDEX_2> identpairs; 844 int cnt = 0; 845 for (i = 1; i <= ident -> GetMaxNr(); i++) 846 { 847 ident -> GetPairs (i, identpairs); 848 cnt += identpairs.Size(); 849 } 850 outfile << cnt << "\n"; 851 for (i = 1; i <= ident -> GetMaxNr(); i++) 852 { 853 ident -> GetPairs (i, identpairs); 854 for (j = 1; j <= identpairs.Size(); j++) 855 { 856 outfile.width (8); 857 outfile << identpairs.Get(j).I1(); 858 outfile.width (8); 859 outfile << identpairs.Get(j).I2(); 860 outfile.width (8); 861 outfile << i << "\n"; 862 } 863 } 864 865 outfile << "identificationtypes\n"; 866 outfile << ident -> GetMaxNr() << "\n"; 867 for (i = 1; i <= ident -> GetMaxNr(); i++) 868 { 869 int type = ident -> GetType(i); 870 outfile << " " << type; 871 } 872 outfile << "\n"; 873 } 874 875 int cntmat = 0; 876 for (i = 1; i <= materials.Size(); i++) 877 if (materials.Get(i) && materials.Get(i)->length()) 878 cntmat++; 879 880 if (cntmat) 881 { 882 outfile << "materials" << endl; 883 outfile << cntmat << endl; 884 for (i = 1; i <= materials.Size(); i++) 885 if (materials.Get(i) && materials.Get(i)->length()) 886 outfile << i << " " << *materials.Get(i) << endl; 887 } 888 889 890 int cntbcnames = 0; 891 for ( int ii = 0; ii < bcnames.Size(); ii++ ) 892 if ( bcnames[ii] ) cntbcnames++; 893 894 if ( cntbcnames ) 895 { 896 outfile << "\n\nbcnames" << endl << bcnames.Size() << endl; 897 for ( i = 0; i < bcnames.Size(); i++ ) 898 outfile << i+1 << "\t" << GetBCName(i) << endl; 899 outfile << endl << endl; 900 } 901 int cntcd2names = 0; 902 for (int ii = 0; ii<cd2names.Size(); ii++) 903 if(cd2names[ii]) cntcd2names++; 904 905 if(cntcd2names) 906 { 907 outfile << "\n\ncd2names" << endl << cd2names.Size() << endl; 908 for (i=0; i<cd2names.Size(); i++) 909 outfile << i+1 << "\t" << GetCD2Name(i) << endl; 910 outfile << endl << endl; 911 } 912 913 int cntcd3names = 0; 914 for (int ii = 0; ii<cd3names.Size(); ii++) 915 if(cd3names[ii]) cntcd3names++; 916 917 if(cntcd3names) 918 { 919 outfile << "\n\ncd3names" << endl << cd3names.Size() << endl; 920 for (i=0; i<cd3names.Size(); i++) 921 outfile << i+1 << "\t" << GetCD3Name(i) << endl; 922 outfile << endl << endl; 923 } 924 925 /* 926 if ( GetDimension() == 2 ) 927 { 928 for (i = 1; i <= GetNSeg(); i++) 929 { 930 const Segment & seg = LineSegment (i); 931 if ( ! bcprops.Contains(seg.si) && seg.GetBCName() != "" ) 932 { 933 bcprops.Append(seg.si); 934 cntbcnames++; 935 } 936 } 937 } 938 else 939 { 940 for (sei = 0; sei < GetNSE(); sei++) 941 { 942 if ((*this)[sei].GetIndex()) 943 { 944 int bcp = GetFaceDescriptor((*this)[sei].GetIndex ()).BCProperty(); 945 string name = GetFaceDescriptor((*this)[sei].GetIndex ()).BCName(); 946 if ( !bcprops.Contains(bcp) && 947 name != "" ) 948 { 949 bcprops.Append(bcp); 950 cntbcnames++; 951 } 952 } 953 } 954 } 955 956 bcprops.SetSize(0); 957 if ( cntbcnames ) 958 { 959 outfile << "\nbcnames" << endl << cntbcnames << endl; 960 if ( GetDimension() == 2 ) 961 { 962 for (i = 1; i <= GetNSeg(); i++) 963 { 964 const Segment & seg = LineSegment (i); 965 if ( ! bcprops.Contains(seg.si) && seg.GetBCName() != "" ) 966 { 967 bcprops.Append(seg.si); 968 outfile << seg.si << "\t" << seg.GetBCName() << endl; 969 } 970 } 971 } 972 else 973 { 974 for (sei = 0; sei < GetNSE(); sei++) 975 { 976 if ((*this)[sei].GetIndex()) 977 { 978 int bcp = GetFaceDescriptor((*this)[sei].GetIndex ()).BCProperty(); 979 string name = GetFaceDescriptor((*this)[sei].GetIndex ()).BCName(); 980 if ( !bcprops.Contains(bcp) && 981 name != "" ) 982 { 983 bcprops.Append(bcp); 984 outfile << bcp << "\t" << name << endl; 985 } 986 } 987 } 988 } 989 outfile << endl << endl; 990 } 991 */ 992 993 int cnt_sing = 0; 994 // for (PointIndex pi = points.Begin(); pi < points.End(); pi++) 995 // if ((*this)[pi].Singularity()>=1.) cnt_sing++; 996 for (auto & p : points) 997 if (p.Singularity() >= 1.) cnt_sing++; 998 999 if (cnt_sing) 1000 { 1001 outfile << "singular_points" << endl << cnt_sing << endl; 1002 // for (PointIndex pi = points.Begin(); pi < points.End(); pi++) 1003 for (PointIndex pi : points.Range()) 1004 if ((*this)[pi].Singularity()>=1.) 1005 outfile << int(pi) << "\t" << (*this)[pi].Singularity() << endl; 1006 } 1007 1008 cnt_sing = 0; 1009 for (SegmentIndex si = 0; si < GetNSeg(); si++) 1010 if ( segments[si].singedge_left ) cnt_sing++; 1011 if (cnt_sing) 1012 { 1013 outfile << "singular_edge_left" << endl << cnt_sing << endl; 1014 for (SegmentIndex si = 0; si < GetNSeg(); si++) 1015 if ( segments[si].singedge_left ) 1016 outfile << int(si) << "\t" << segments[si].singedge_left << endl; 1017 } 1018 1019 cnt_sing = 0; 1020 for (SegmentIndex si = 0; si < GetNSeg(); si++) 1021 if ( segments[si].singedge_right ) cnt_sing++; 1022 if (cnt_sing) 1023 { 1024 outfile << "singular_edge_right" << endl << cnt_sing << endl; 1025 for (SegmentIndex si = 0; si < GetNSeg(); si++) 1026 if ( segments[si].singedge_right ) 1027 outfile << int(si) << "\t" << segments[si].singedge_right << endl; 1028 } 1029 1030 1031 cnt_sing = 0; 1032 for (SurfaceElementIndex sei = 0; sei < GetNSE(); sei++) 1033 if ( GetFaceDescriptor ((*this)[sei].GetIndex()).domin_singular) 1034 cnt_sing++; 1035 1036 if (cnt_sing) 1037 { 1038 outfile << "singular_face_inside" << endl << cnt_sing << endl; 1039 for (SurfaceElementIndex sei = 0; sei < GetNSE(); sei++) 1040 if ( GetFaceDescriptor ((*this)[sei].GetIndex()).domin_singular) 1041 outfile << int(sei) << "\t" << 1042 GetFaceDescriptor ((*this)[sei].GetIndex()).domin_singular << endl; 1043 } 1044 1045 cnt_sing = 0; 1046 for (SurfaceElementIndex sei = 0; sei < GetNSE(); sei++) 1047 if ( GetFaceDescriptor ((*this)[sei].GetIndex()).domout_singular) cnt_sing++; 1048 if (cnt_sing) 1049 { 1050 outfile << "singular_face_outside" << endl << cnt_sing << endl; 1051 for (SurfaceElementIndex sei = 0; sei < GetNSE(); sei++) 1052 if ( GetFaceDescriptor ((*this)[sei].GetIndex()).domout_singular) 1053 outfile << int(sei) << "\t" 1054 << GetFaceDescriptor ((*this)[sei].GetIndex()).domout_singular << endl; 1055 } 1056 1057 1058 // Philippose - 09/07/2009 1059 // Add mesh face colours to Netgen Vol file format 1060 // The colours are saved in RGB triplets 1061 int cnt_facedesc = GetNFD(); 1062 if (cnt_facedesc) 1063 { 1064 outfile << endl << endl << "# Surfnr Red Green Blue" << endl; 1065 outfile << "face_colours" << endl << cnt_facedesc << endl; 1066 1067 outfile.precision(8); 1068 outfile.setf(ios::fixed, ios::floatfield); 1069 outfile.setf(ios::showpoint); 1070 1071 for(i = 1; i <= cnt_facedesc; i++) 1072 { 1073 outfile.width(8); 1074 outfile << GetFaceDescriptor(i).SurfNr()+1 << " "; 1075 outfile.width(12); 1076 outfile << GetFaceDescriptor(i).SurfColour()[0] << " "; 1077 outfile.width(12); 1078 outfile << GetFaceDescriptor(i).SurfColour()[1] << " "; 1079 outfile.width(12); 1080 outfile << GetFaceDescriptor(i).SurfColour()[2]; 1081 outfile << endl; 1082 } 1083 } 1084 1085 outfile << endl << endl << "endmesh" << endl << endl; 1086 if (geometry) 1087 geometry -> SaveToMeshFile (outfile); 1088 } 1089 1090 1091 Load(const string & filename)1092 void Mesh :: Load (const string & filename) 1093 { 1094 cout << "filename = " << filename << endl; 1095 1096 if (filename.find(".vol.bin") != string::npos) 1097 { 1098 BinaryInArchive in(filename); 1099 in & (*this); 1100 return; 1101 } 1102 1103 istream * infile = NULL; 1104 1105 if (filename.find(".vol.gz") != string::npos) 1106 infile = new igzstream (filename.c_str()); 1107 else 1108 infile = new ifstream (filename.c_str()); 1109 1110 // ifstream infile(filename.c_str()); 1111 if (! (infile -> good()) ) 1112 throw NgException ("mesh file not found"); 1113 1114 Load(*infile); 1115 delete infile; 1116 } 1117 1118 1119 1120 // Reads mandatory integer and optional string token from input stream 1121 // used for parsing bcnames, cd2names etc. ReadNumberAndName(istream & infile,int & i,string & s)1122 void ReadNumberAndName( istream & infile, int & i, string & s ) 1123 { 1124 string line; 1125 std::istringstream iline; 1126 1127 bool empty_line = true; 1128 1129 while(empty_line && infile) 1130 { 1131 std::getline(infile, line); 1132 iline = std::istringstream{line}; 1133 iline >> i; 1134 1135 if(iline) 1136 empty_line = false; 1137 1138 iline >> s; 1139 } 1140 1141 if(!infile) 1142 throw Exception("Reached end of file while parsing"); 1143 } 1144 Load(istream & infile)1145 void Mesh :: Load (istream & infile) 1146 { 1147 static Timer timer("Mesh::Load"); RegionTimer rt(timer); 1148 if (! (infile.good()) ) 1149 { 1150 cout << "cannot load mesh" << endl; 1151 throw NgException ("mesh file not found"); 1152 } 1153 1154 int rank = GetCommunicator().Rank(); 1155 int ntasks = GetCommunicator().Size(); 1156 1157 char str[100]; 1158 int i, n; 1159 1160 double scale = 1; // globflags.GetNumFlag ("scale", 1); 1161 int inverttets = 0; // globflags.GetDefineFlag ("inverttets"); 1162 int invertsurf = 0; // globflags.GetDefineFlag ("invertsurfacemesh"); 1163 1164 1165 facedecoding.SetSize(0); 1166 1167 bool endmesh = false; 1168 1169 1170 while (infile.good() && !endmesh) 1171 { 1172 infile >> str; 1173 if (strcmp (str, "dimension") == 0) 1174 { 1175 infile >> dimension; 1176 } 1177 1178 if (strcmp (str, "geomtype") == 0) 1179 { 1180 int hi; 1181 infile >> hi; 1182 geomtype = GEOM_TYPE(hi); 1183 } 1184 1185 1186 if (strcmp (str, "surfaceelements") == 0 || strcmp (str, "surfaceelementsgi")==0 || strcmp (str, "surfaceelementsuv") == 0) 1187 { 1188 static Timer t1("read surface elements"); RegionTimer rt1(t1); 1189 infile >> n; 1190 PrintMessage (3, n, " surface elements"); 1191 1192 bool geominfo = strcmp (str, "surfaceelementsgi") == 0; 1193 bool uv = strcmp (str, "surfaceelementsuv") == 0; 1194 1195 1196 for (i = 1; i <= n; i++) 1197 { 1198 int surfnr, bcp, domin, domout, nep, faceind = 0; 1199 1200 infile >> surfnr >> bcp >> domin >> domout; 1201 surfnr--; 1202 1203 bool invert_el = false; 1204 /* 1205 if (domin == 0) 1206 { 1207 invert_el = true; 1208 Swap (domin, domout); 1209 } 1210 */ 1211 1212 for (int j = 1; j <= facedecoding.Size(); j++) 1213 if (GetFaceDescriptor(j).SurfNr() == surfnr && 1214 GetFaceDescriptor(j).BCProperty() == bcp && 1215 GetFaceDescriptor(j).DomainIn() == domin && 1216 GetFaceDescriptor(j).DomainOut() == domout) 1217 faceind = j; 1218 1219 // if (facedecoding.Size()) faceind = 1; // for timing 1220 1221 if (!faceind) 1222 { 1223 faceind = AddFaceDescriptor (FaceDescriptor(surfnr, domin, domout, 0)); 1224 GetFaceDescriptor(faceind).SetBCProperty (bcp); 1225 } 1226 1227 infile >> nep; 1228 if (!nep) nep = 3; 1229 1230 Element2d tri(nep); 1231 tri.SetIndex(faceind); 1232 1233 for (int j = 1; j <= nep; j++) 1234 infile >> tri.PNum(j); 1235 1236 if (geominfo) 1237 for (int j = 1; j <= nep; j++) 1238 infile >> tri.GeomInfoPi(j).trignum; 1239 1240 if (uv) 1241 for (int j = 1; j <= nep; j++) 1242 infile >> tri.GeomInfoPi(j).u >> tri.GeomInfoPi(j).v; 1243 1244 if (invertsurf) tri.Invert(); 1245 if (invert_el) tri.Invert(); 1246 1247 AddSurfaceElement (tri); 1248 } 1249 } 1250 1251 if (strcmp (str, "volumeelements") == 0) 1252 { 1253 static Timer t1("read volume elements"); RegionTimer rt1(t1); 1254 infile >> n; 1255 PrintMessage (3, n, " volume elements"); 1256 for (i = 1; i <= n; i++) 1257 { 1258 Element el(TET); 1259 int hi, nep; 1260 infile >> hi; 1261 if (hi == 0) hi = 1; 1262 el.SetIndex(hi); 1263 infile >> nep; 1264 el.SetNP(nep); 1265 el.SetCurved (nep != 4); 1266 for (int j = 0; j < nep; j++) 1267 infile >> (int&)(el[j]); 1268 1269 if (inverttets) 1270 el.Invert(); 1271 1272 AddVolumeElement (el); 1273 } 1274 } 1275 1276 1277 if (strcmp (str, "edgesegments") == 0) 1278 { 1279 static Timer t1("read edge segments"); RegionTimer rt1(t1); 1280 infile >> n; 1281 for (i = 1; i <= n; i++) 1282 { 1283 Segment seg; 1284 int hi; 1285 infile >> seg.si >> hi >> seg[0] >> seg[1]; 1286 AddSegment (seg); 1287 } 1288 } 1289 1290 1291 1292 if (strcmp (str, "edgesegmentsgi") == 0) 1293 { 1294 static Timer t1("read edge segmentsgi"); RegionTimer rt1(t1); 1295 infile >> n; 1296 for (i = 1; i <= n; i++) 1297 { 1298 Segment seg; 1299 int hi; 1300 infile >> seg.si >> hi >> seg[0] >> seg[1] 1301 >> seg.geominfo[0].trignum 1302 >> seg.geominfo[1].trignum; 1303 AddSegment (seg); 1304 } 1305 } 1306 1307 if (strcmp (str, "edgesegmentsgi2") == 0) 1308 { 1309 static Timer t1("read edge segmentsgi2"); RegionTimer rt1(t1); 1310 int a; 1311 infile >> a; 1312 n=a; 1313 1314 PrintMessage (3, n, " curve elements"); 1315 1316 for (i = 1; i <= n; i++) 1317 { 1318 Segment seg; 1319 int hi; 1320 infile >> seg.si >> hi >> seg[0] >> seg[1] 1321 >> seg.geominfo[0].trignum 1322 >> seg.geominfo[1].trignum 1323 >> seg.surfnr1 >> seg.surfnr2 1324 >> seg.edgenr 1325 >> seg.epgeominfo[0].dist 1326 >> seg.epgeominfo[1].edgenr 1327 >> seg.epgeominfo[1].dist; 1328 1329 seg.epgeominfo[0].edgenr = seg.epgeominfo[1].edgenr; 1330 1331 seg.domin = seg.surfnr1; 1332 seg.domout = seg.surfnr2; 1333 1334 seg.surfnr1--; 1335 seg.surfnr2--; 1336 1337 AddSegment (seg); 1338 } 1339 } 1340 1341 if (strcmp (str, "points") == 0) 1342 { 1343 static Timer t1("read points"); RegionTimer rt1(t1); 1344 infile >> n; 1345 PrintMessage (3, n, " points"); 1346 for (i = 1; i <= n; i++) 1347 { 1348 Point3d p; 1349 infile >> p.X() >> p.Y() >> p.Z(); 1350 p.X() *= scale; 1351 p.Y() *= scale; 1352 p.Z() *= scale; 1353 AddPoint (p); 1354 } 1355 PrintMessage (3, n, " points done"); 1356 } 1357 1358 if (strcmp (str, "pointelements") == 0) 1359 { 1360 static Timer t1("read point elements"); RegionTimer rt1(t1); 1361 infile >> n; 1362 PrintMessage (3, n, " pointelements"); 1363 for (i = 1; i <= n; i++) 1364 { 1365 Element0d el; 1366 infile >> el.pnum >> el.index; 1367 pointelements.Append (el); 1368 } 1369 PrintMessage (3, n, " pointelements done"); 1370 } 1371 1372 if (strcmp (str, "identifications") == 0) 1373 { 1374 infile >> n; 1375 PrintMessage (3, n, " identifications"); 1376 for (i = 1; i <= n; i++) 1377 { 1378 PointIndex pi1, pi2; 1379 int ind; 1380 infile >> pi1 >> pi2 >> ind; 1381 ident -> Add (pi1, pi2, ind); 1382 } 1383 } 1384 1385 if (strcmp (str, "identificationtypes") == 0) 1386 { 1387 infile >> n; 1388 PrintMessage (3, n, " identificationtypes"); 1389 for (i = 1; i <= n; i++) 1390 { 1391 int type; 1392 infile >> type; 1393 ident -> SetType(i,Identifications::ID_TYPE(type)); 1394 } 1395 } 1396 1397 if (strcmp (str, "materials") == 0) 1398 { 1399 infile >> n; 1400 for ( auto i : Range(n) ) 1401 { 1402 int nr; 1403 string mat; 1404 ReadNumberAndName( infile, nr, mat ); 1405 SetMaterial (nr, mat.c_str()); 1406 } 1407 } 1408 1409 if ( strcmp (str, "bcnames" ) == 0 ) 1410 { 1411 infile >> n; 1412 Array<int> bcnrs(n); 1413 SetNBCNames(n); 1414 for ( auto i : Range(n) ) 1415 { 1416 string nextbcname; 1417 ReadNumberAndName( infile, bcnrs[i], nextbcname ); 1418 bcnames[bcnrs[i]-1] = new string(nextbcname); 1419 } 1420 1421 if ( GetDimension() == 3 ) 1422 { 1423 for (SurfaceElementIndex sei = 0; sei < GetNSE(); sei++) 1424 { 1425 if ((*this)[sei].GetIndex()) 1426 { 1427 int bcp = GetFaceDescriptor((*this)[sei].GetIndex ()).BCProperty(); 1428 if ( bcp <= n ) 1429 GetFaceDescriptor((*this)[sei].GetIndex ()).SetBCName(bcnames[bcp-1]); 1430 else 1431 GetFaceDescriptor((*this)[sei].GetIndex ()).SetBCName(0); 1432 1433 } 1434 } 1435 1436 } 1437 } 1438 1439 if ( strcmp (str, "cd2names" ) == 0) 1440 { 1441 infile >> n; 1442 Array<int> cd2nrs(n); 1443 SetNCD2Names(n); 1444 for ( auto i : Range(n) ) 1445 { 1446 string nextcd2name; 1447 ReadNumberAndName( infile, cd2nrs[i], nextcd2name ); 1448 cd2names[cd2nrs[i]-1] = new string(nextcd2name); 1449 } 1450 if (GetDimension() < 2) 1451 { 1452 throw NgException("co dim 2 elements not implemented for dimension < 2"); 1453 } 1454 } 1455 1456 if ( strcmp (str, "cd3names" ) == 0) 1457 { 1458 infile >> n; 1459 Array<int> cd3nrs(n); 1460 SetNCD3Names(n); 1461 for( auto i : Range(n) ) 1462 { 1463 string nextcd3name; 1464 ReadNumberAndName( infile, cd3nrs[i], nextcd3name ); 1465 infile >> cd3nrs[i-1] >> nextcd3name; 1466 cd3names[cd3nrs[i-1]-1] = new string(nextcd3name); 1467 } 1468 if (GetDimension() < 3) 1469 { 1470 throw NgException("co dim 3 elements not implemented for dimension < 3"); 1471 } 1472 } 1473 1474 if (strcmp (str, "singular_points") == 0) 1475 { 1476 infile >> n; 1477 for (i = 1; i <= n; i++) 1478 { 1479 PointIndex pi; 1480 double s; 1481 infile >> pi; 1482 infile >> s; 1483 (*this)[pi].Singularity (s); 1484 } 1485 } 1486 1487 if (strcmp (str, "singular_edge_left") == 0) 1488 { 1489 infile >> n; 1490 for (i = 1; i <= n; i++) 1491 { 1492 SegmentIndex si; 1493 double s; 1494 infile >> si; 1495 infile >> s; 1496 (*this)[si].singedge_left = s; 1497 } 1498 } 1499 if (strcmp (str, "singular_edge_right") == 0) 1500 { 1501 infile >> n; 1502 for (i = 1; i <= n; i++) 1503 { 1504 SegmentIndex si; 1505 double s; 1506 infile >> si; 1507 infile >> s; 1508 (*this)[si].singedge_right = s; 1509 } 1510 } 1511 1512 if (strcmp (str, "singular_face_inside") == 0) 1513 { 1514 infile >> n; 1515 for (i = 1; i <= n; i++) 1516 { 1517 SurfaceElementIndex sei; 1518 double s; 1519 infile >> sei; 1520 infile >> s; 1521 GetFaceDescriptor((*this)[sei].GetIndex()).domin_singular = s; 1522 } 1523 } 1524 1525 if (strcmp (str, "singular_face_outside") == 0) 1526 { 1527 infile >> n; 1528 for (i = 1; i <= n; i++) 1529 { 1530 SurfaceElementIndex sei; 1531 double s; 1532 infile >> sei; 1533 infile >> s; 1534 GetFaceDescriptor((*this)[sei].GetIndex()).domout_singular = s; 1535 } 1536 } 1537 1538 // Philippose - 09/07/2009 1539 // Add mesh face colours to Netgen Vol file format 1540 // The colours are read in as RGB triplets 1541 if (strcmp (str, "face_colours") == 0) 1542 { 1543 int cnt_facedesc = GetNFD(); 1544 infile >> n; 1545 if(n == cnt_facedesc) 1546 { 1547 for(i = 1; i <= n; i++) 1548 { 1549 int surfnr = 0; 1550 Vec3d surfcolour(0.0,1.0,0.0); 1551 1552 infile >> surfnr 1553 >> surfcolour.X() 1554 >> surfcolour.Y() 1555 >> surfcolour.Z(); 1556 1557 surfnr--; 1558 1559 if(surfnr > 0) 1560 { 1561 for(int facedesc = 1; facedesc <= cnt_facedesc; facedesc++) 1562 { 1563 if(surfnr == GetFaceDescriptor(facedesc).SurfNr()) 1564 { 1565 GetFaceDescriptor(facedesc).SetSurfColour(surfcolour); 1566 } 1567 } 1568 } 1569 } 1570 } 1571 } 1572 1573 if (strcmp (str, "endmesh") == 0) 1574 endmesh = true; 1575 1576 1577 1578 strcpy (str, ""); 1579 } 1580 1581 1582 1583 1584 CalcSurfacesOfNode (); 1585 1586 if (ntasks == 1) // sequential run only 1587 { 1588 topology.Update(); 1589 clusters -> Update(); 1590 } 1591 1592 SetNextMajorTimeStamp(); 1593 // PrintMemInfo (cout); 1594 } 1595 1596 DoArchive(Archive & archive)1597 void Mesh :: DoArchive (Archive & archive) 1598 { 1599 static Timer t("Mesh::Archive"); RegionTimer r(t); 1600 1601 #ifdef PARALLEL 1602 auto comm = GetCommunicator(); 1603 if (archive.IsParallel() && comm.Size() > 1) 1604 { // parallel pickling supported only for output archives 1605 if (comm.Rank() == 0) 1606 archive & dimension; 1607 1608 auto rank = comm.Rank(); 1609 1610 auto & partop = GetParallelTopology(); 1611 1612 // global enumration of points: 1613 // not used now, but will be needed for refined meshes 1614 // GridFunciton pickling is not compatible, now 1615 // should go to paralleltopology 1616 1617 1618 1619 // merge points 1620 Array<PointIndex, PointIndex> globnum(points.Size()); 1621 PointIndex maxglob = -1; 1622 for (auto pi : Range(points)) 1623 { 1624 globnum[pi] = partop.GetGlobalPNum(pi); 1625 // globnum[pi] = global_pnums[pi]; 1626 maxglob = max(globnum[pi], maxglob); 1627 } 1628 1629 maxglob = comm.AllReduce (maxglob, MPI_MAX); 1630 int numglob = maxglob+1-PointIndex::BASE; 1631 if (comm.Rank() > 0) 1632 { 1633 comm.Send (globnum, 0, 200); 1634 comm.Send (points, 0, 200); 1635 } 1636 else 1637 { 1638 Array<PointIndex, PointIndex> globnumi; 1639 Array<MeshPoint, PointIndex> pointsi; 1640 Array<MeshPoint, PointIndex> globpoints(numglob); 1641 for (int j = 1; j < comm.Size(); j++) 1642 { 1643 comm.Recv (globnumi, j, 200); 1644 comm.Recv (pointsi, j, 200); 1645 for (auto i : Range(globnumi)) 1646 globpoints[globnumi[i]] = pointsi[i]; 1647 } 1648 archive & globpoints; 1649 } 1650 1651 1652 // sending surface elements 1653 auto copy_el2d (surfelements); 1654 for (auto & el : copy_el2d) 1655 for (auto & pi : el.PNums()) 1656 pi = globnum[pi]; 1657 1658 if (comm.Rank() > 0) 1659 comm.Send(copy_el2d, 0, 200); 1660 else 1661 { 1662 Array<Element2d, SurfaceElementIndex> el2di; 1663 for (int j = 1; j < comm.Size(); j++) 1664 { 1665 comm.Recv(el2di, j, 200); 1666 for (auto & el : el2di) 1667 copy_el2d += el; 1668 } 1669 archive & copy_el2d; 1670 } 1671 1672 1673 // sending volume elements 1674 auto copy_el3d (volelements); 1675 for (auto & el : copy_el3d) 1676 for (auto & pi : el.PNums()) 1677 pi = globnum[pi]; 1678 1679 if (comm.Rank() > 0) 1680 comm.Send(copy_el3d, 0, 200); 1681 else 1682 { 1683 Array<Element, ElementIndex> el3di; 1684 for (int j = 1; j < comm.Size(); j++) 1685 { 1686 comm.Recv(el3di, j, 200); 1687 for (auto & el : el3di) 1688 copy_el3d += el; 1689 } 1690 archive & copy_el3d; 1691 } 1692 1693 1694 // sending 1D elements 1695 auto copy_el1d (segments); 1696 for (auto & el : copy_el1d) 1697 for (auto & pi : el.pnums) 1698 if (pi != PointIndex(PointIndex::INVALID)) 1699 pi = globnum[pi]; 1700 1701 if (comm.Rank() > 0) 1702 comm.Send(copy_el1d, 0, 200); 1703 else 1704 { 1705 Array<Segment, SegmentIndex> el1di; 1706 for (int j = 1; j < comm.Size(); j++) 1707 { 1708 comm.Recv(el1di, j, 200); 1709 for (auto & el : el1di) 1710 copy_el1d += el; 1711 } 1712 archive & copy_el1d; 1713 } 1714 1715 if (comm.Rank() == 0) 1716 { 1717 archive & facedecoding; 1718 archive & materials & bcnames & cd2names & cd3names; 1719 auto mynv = numglob; 1720 archive & mynv; // numvertices; 1721 archive & *ident; 1722 1723 if(archive.GetVersion("netgen") >= "v6.2.2103-1") 1724 { 1725 archive.NeedsVersion("netgen", "v6.2.2103-1"); 1726 archive & vol_partition & surf_partition & seg_partition; 1727 } 1728 1729 archive.Shallow(geometry); 1730 archive & *curvedelems; 1731 } 1732 1733 if (comm.Rank() == 0) 1734 return; 1735 } 1736 #endif 1737 1738 1739 archive & dimension; 1740 archive & points; 1741 archive & surfelements; 1742 archive & volelements; 1743 archive & segments; 1744 archive & facedecoding; 1745 archive & materials & bcnames & cd2names & cd3names; 1746 archive & numvertices; 1747 1748 archive & *ident; 1749 1750 // cout << "archive, ngsversion = " << archive.GetVersion("netgen") << endl; 1751 if(archive.GetVersion("netgen") >= "v6.2.2103-1") 1752 { 1753 // cout << "do the partition" << endl; 1754 archive.NeedsVersion("netgen", "v6.2.2103-1"); 1755 archive & vol_partition & surf_partition & seg_partition; 1756 } 1757 // else 1758 // cout << "no partition" << endl; 1759 1760 archive.Shallow(geometry); 1761 archive & *curvedelems; 1762 1763 if (archive.Input()) 1764 { 1765 int rank = GetCommunicator().Rank(); 1766 int ntasks = GetCommunicator().Size(); 1767 1768 RebuildSurfaceElementLists(); 1769 1770 CalcSurfacesOfNode (); 1771 if (ntasks == 1) // sequential run only 1772 { 1773 topology.Update(); 1774 clusters -> Update(); 1775 } 1776 SetNextMajorTimeStamp(); 1777 } 1778 } 1779 1780 Merge(const string & filename,const int surfindex_offset)1781 void Mesh :: Merge (const string & filename, const int surfindex_offset) 1782 { 1783 ifstream infile(filename.c_str()); 1784 if (!infile.good()) 1785 throw NgException ("mesh file not found"); 1786 1787 Merge(infile,surfindex_offset); 1788 1789 } 1790 1791 1792 Merge(istream & infile,const int surfindex_offset)1793 void Mesh :: Merge (istream & infile, const int surfindex_offset) 1794 { 1795 char str[100]; 1796 int i, n; 1797 1798 1799 int inverttets = 0; // globflags.GetDefineFlag ("inverttets"); 1800 1801 int oldnp = GetNP(); 1802 int oldne = GetNSeg(); 1803 int oldnd = GetNDomains(); 1804 1805 for(SurfaceElementIndex si = 0; si < GetNSE(); si++) 1806 for(int j=1; j<=(*this)[si].GetNP(); j++) (*this)[si].GeomInfoPi(j).trignum = -1; 1807 1808 int max_surfnr = 0; 1809 for (i = 1; i <= GetNFD(); i++) 1810 max_surfnr = max2 (max_surfnr, GetFaceDescriptor(i).SurfNr()); 1811 max_surfnr++; 1812 1813 if(max_surfnr < surfindex_offset) max_surfnr = surfindex_offset; 1814 1815 1816 bool endmesh = false; 1817 1818 while (infile.good() && !endmesh) 1819 { 1820 infile >> str; 1821 1822 if (strcmp (str, "surfaceelementsgi") == 0 || strcmp (str, "surfaceelements") == 0) 1823 { 1824 infile >> n; 1825 PrintMessage (3, n, " surface elements"); 1826 for (i = 1; i <= n; i++) 1827 { 1828 int j; 1829 int surfnr, bcp, domin, domout, nep, faceind = 0; 1830 infile >> surfnr >> bcp >> domin >> domout; 1831 1832 surfnr--; 1833 1834 if(domin > 0) domin += oldnd; 1835 if(domout > 0) domout += oldnd; 1836 surfnr += max_surfnr; 1837 1838 1839 for (j = 1; j <= facedecoding.Size(); j++) 1840 if (GetFaceDescriptor(j).SurfNr() == surfnr && 1841 GetFaceDescriptor(j).BCProperty() == bcp && 1842 GetFaceDescriptor(j).DomainIn() == domin && 1843 GetFaceDescriptor(j).DomainOut() == domout) 1844 faceind = j; 1845 1846 if (!faceind) 1847 { 1848 faceind = AddFaceDescriptor (FaceDescriptor(surfnr, domin, domout, 0)); 1849 if(GetDimension() == 2) bcp++; 1850 GetFaceDescriptor(faceind).SetBCProperty (bcp); 1851 } 1852 1853 infile >> nep; 1854 if (!nep) nep = 3; 1855 1856 Element2d tri(nep); 1857 tri.SetIndex(faceind); 1858 1859 for (j = 1; j <= nep; j++) 1860 { 1861 infile >> tri.PNum(j); 1862 tri.PNum(j) = tri.PNum(j) + oldnp; 1863 } 1864 1865 1866 if (strcmp (str, "surfaceelementsgi") == 0) 1867 for (j = 1; j <= nep; j++) 1868 { 1869 infile >> tri.GeomInfoPi(j).trignum; 1870 tri.GeomInfoPi(j).trignum = -1; 1871 } 1872 1873 AddSurfaceElement (tri); 1874 } 1875 } 1876 1877 1878 if (strcmp (str, "edgesegments") == 0) 1879 { 1880 infile >> n; 1881 for (i = 1; i <= n; i++) 1882 { 1883 Segment seg; 1884 int hi; 1885 infile >> seg.si >> hi >> seg[0] >> seg[1]; 1886 seg[0] = seg[0] + oldnp; 1887 seg[1] = seg[1] + oldnp; 1888 AddSegment (seg); 1889 } 1890 } 1891 1892 1893 1894 if (strcmp (str, "edgesegmentsgi") == 0) 1895 { 1896 infile >> n; 1897 for (i = 1; i <= n; i++) 1898 { 1899 Segment seg; 1900 int hi; 1901 infile >> seg.si >> hi >> seg[0] >> seg[1] 1902 >> seg.geominfo[0].trignum 1903 >> seg.geominfo[1].trignum; 1904 seg[0] = seg[0] + oldnp; 1905 seg[1] = seg[1] + oldnp; 1906 AddSegment (seg); 1907 } 1908 } 1909 if (strcmp (str, "edgesegmentsgi2") == 0) 1910 { 1911 infile >> n; 1912 PrintMessage (3, n, " curve elements"); 1913 1914 for (i = 1; i <= n; i++) 1915 { 1916 Segment seg; 1917 int hi; 1918 infile >> seg.si >> hi >> seg[0] >> seg[1] 1919 >> seg.geominfo[0].trignum 1920 >> seg.geominfo[1].trignum 1921 >> seg.surfnr1 >> seg.surfnr2 1922 >> seg.edgenr 1923 >> seg.epgeominfo[0].dist 1924 >> seg.epgeominfo[1].edgenr 1925 >> seg.epgeominfo[1].dist; 1926 seg.epgeominfo[0].edgenr = seg.epgeominfo[1].edgenr; 1927 1928 seg.surfnr1--; 1929 seg.surfnr2--; 1930 1931 if(seg.surfnr1 >= 0) seg.surfnr1 = seg.surfnr1 + max_surfnr; 1932 if(seg.surfnr2 >= 0) seg.surfnr2 = seg.surfnr2 + max_surfnr; 1933 seg[0] = seg[0] +oldnp; 1934 seg[1] = seg[1] +oldnp; 1935 *testout << "old edgenr: " << seg.edgenr << endl; 1936 seg.edgenr = seg.edgenr + oldne; 1937 *testout << "new edgenr: " << seg.edgenr << endl; 1938 seg.epgeominfo[1].edgenr = seg.epgeominfo[1].edgenr + oldne; 1939 1940 AddSegment (seg); 1941 } 1942 } 1943 1944 if (strcmp (str, "volumeelements") == 0) 1945 { 1946 infile >> n; 1947 PrintMessage (3, n, " volume elements"); 1948 for (i = 1; i <= n; i++) 1949 { 1950 Element el(TET); 1951 int hi, nep; 1952 infile >> hi; 1953 if (hi == 0) hi = 1; 1954 el.SetIndex(hi+oldnd); 1955 infile >> nep; 1956 el.SetNP(nep); 1957 1958 for (int j = 0; j < nep; j++) 1959 { 1960 infile >> (int&)(el[j]); 1961 el[j] = el[j]+oldnp; 1962 } 1963 1964 if (inverttets) 1965 el.Invert(); 1966 1967 AddVolumeElement (el); 1968 } 1969 } 1970 1971 1972 if (strcmp (str, "points") == 0) 1973 { 1974 infile >> n; 1975 PrintMessage (3, n, " points"); 1976 for (i = 1; i <= n; i++) 1977 { 1978 Point3d p; 1979 infile >> p.X() >> p.Y() >> p.Z(); 1980 AddPoint (p); 1981 } 1982 } 1983 1984 1985 if (strcmp (str, "endmesh") == 0) 1986 { 1987 endmesh = true; 1988 } 1989 1990 1991 if (strcmp (str, "materials") == 0) 1992 { 1993 infile >> n; 1994 for (i = 1; i <= n; i++) 1995 { 1996 int nr; 1997 string mat; 1998 infile >> nr >> mat; 1999 SetMaterial (nr+oldnd, mat.c_str()); 2000 } 2001 } 2002 2003 2004 strcpy (str, ""); 2005 } 2006 2007 CalcSurfacesOfNode (); 2008 2009 topology.Update(); 2010 clusters -> Update(); 2011 2012 SetNextMajorTimeStamp(); 2013 } 2014 2015 2016 2017 2018 2019 2020 2021 2022 2023 TestOk() const2024 bool Mesh :: TestOk () const 2025 { 2026 for (ElementIndex ei = 0; ei < volelements.Size(); ei++) 2027 { 2028 for (int j = 0; j < 4; j++) 2029 if ( (*this)[ei][j] <= PointIndex::BASE-1) 2030 { 2031 (*testout) << "El " << ei << " has 0 nodes: "; 2032 for (int k = 0; k < 4; k++) 2033 (*testout) << (*this)[ei][k]; 2034 break; 2035 } 2036 } 2037 CheckMesh3D (*this); 2038 return 1; 2039 } 2040 SetAllocSize(int nnodes,int nsegs,int nsel,int nel)2041 void Mesh :: SetAllocSize(int nnodes, int nsegs, int nsel, int nel) 2042 { 2043 points.SetAllocSize(nnodes); 2044 segments.SetAllocSize(nsegs); 2045 surfelements.SetAllocSize(nsel); 2046 volelements.SetAllocSize(nel); 2047 } 2048 BuildBoundaryEdges(bool rebuild)2049 void Mesh :: BuildBoundaryEdges(bool rebuild) 2050 { 2051 static Timer t("Mesh::BuildBoundaryEdges"); RegionTimer reg(t); 2052 2053 if(!rebuild && boundaryedges) 2054 return; 2055 2056 boundaryedges = make_unique<INDEX_2_CLOSED_HASHTABLE<int>> 2057 (3 * (GetNSE() + GetNOpenElements()) + GetNSeg() + 1); 2058 2059 2060 for (SurfaceElementIndex sei = 0; sei < GetNSE(); sei++) 2061 { 2062 const Element2d & sel = surfelements[sei]; 2063 if (sel.IsDeleted()) continue; 2064 2065 // int si = sel.GetIndex(); 2066 2067 if (sel.GetNP() <= 4) 2068 for (int j = 0; j < sel.GetNP(); j++) 2069 { 2070 INDEX_2 i2; 2071 i2.I1() = sel.PNumMod(j+1); 2072 i2.I2() = sel.PNumMod(j+2); 2073 i2.Sort(); 2074 boundaryedges->Set (i2, 1); 2075 } 2076 else if (sel.GetType()==TRIG6) 2077 { 2078 for (int j = 0; j < 3; j++) 2079 { 2080 INDEX_2 i2; 2081 i2.I1() = sel[j]; 2082 i2.I2() = sel[(j+1)%3]; 2083 i2.Sort(); 2084 boundaryedges->Set (i2, 1); 2085 } 2086 } 2087 else 2088 cerr << "illegal element for buildboundaryedges" << endl; 2089 } 2090 2091 2092 for (int i = 0; i < openelements.Size(); i++) 2093 { 2094 const Element2d & sel = openelements[i]; 2095 for (int j = 0; j < sel.GetNP(); j++) 2096 { 2097 INDEX_2 i2; 2098 i2.I1() = sel.PNumMod(j+1); 2099 i2.I2() = sel.PNumMod(j+2); 2100 i2.Sort(); 2101 boundaryedges->Set (i2, 1); 2102 2103 points[sel[j]].SetType(FIXEDPOINT); 2104 } 2105 } 2106 2107 for (int i = 0; i < GetNSeg(); i++) 2108 { 2109 const Segment & seg = segments[i]; 2110 INDEX_2 i2(seg[0], seg[1]); 2111 i2.Sort(); 2112 2113 boundaryedges -> Set (i2, 2); 2114 //segmentht -> Set (i2, i); 2115 } 2116 2117 2118 } 2119 CalcSurfacesOfNode()2120 void Mesh :: CalcSurfacesOfNode () 2121 { 2122 static Timer t("Mesh::CalcSurfacesOfNode"); RegionTimer reg (t); 2123 static Timer tn2se("Mesh::CalcSurfacesOfNode - surf on node"); 2124 static Timer tht("Mesh::CalcSurfacesOfNode - surfelementht"); 2125 // surfacesonnode.SetSize (GetNP()); 2126 TABLE<int,PointIndex::BASE> surfacesonnode(GetNP()); 2127 2128 // delete boundaryedges; 2129 // boundaryedges = NULL; 2130 boundaryedges = nullptr; 2131 2132 // delete surfelementht; 2133 // surfelementht = nullptr; 2134 surfelementht = nullptr; 2135 // delete segmentht; 2136 2137 /* 2138 surfelementht = new INDEX_3_HASHTABLE<int> (GetNSE()/4 + 1); 2139 segmentht = new INDEX_2_HASHTABLE<int> (GetNSeg() + 1); 2140 */ 2141 2142 if (dimension == 3) 2143 surfelementht = make_unique<INDEX_3_CLOSED_HASHTABLE<int>> (3*GetNSE() + 1); 2144 segmentht = make_unique<INDEX_2_CLOSED_HASHTABLE<int>> (3*GetNSeg() + 1); 2145 2146 tn2se.Start(); 2147 if (dimension == 3) 2148 /* 2149 for (SurfaceElementIndex sei = 0; sei < GetNSE(); sei++) 2150 { 2151 const Element2d & sel = surfelements[sei]; 2152 */ 2153 for (const Element2d & sel : surfelements) 2154 { 2155 if (sel.IsDeleted()) continue; 2156 2157 int si = sel.GetIndex(); 2158 2159 /* 2160 for (int j = 0; j < sel.GetNP(); j++) 2161 { 2162 PointIndex pi = sel[j]; 2163 */ 2164 for (PointIndex pi : sel.PNums()) 2165 { 2166 if (!surfacesonnode[pi].Contains(si)) 2167 surfacesonnode.Add (pi, si); 2168 /* 2169 bool found = 0; 2170 for (int k = 0; k < surfacesonnode[pi].Size(); k++) 2171 if (surfacesonnode[pi][k] == si) 2172 { 2173 found = 1; 2174 break; 2175 } 2176 2177 if (!found) 2178 surfacesonnode.Add (pi, si); 2179 */ 2180 } 2181 } 2182 /* 2183 for (sei = 0; sei < GetNSE(); sei++) 2184 { 2185 const Element2d & sel = surfelements[sei]; 2186 if (sel.IsDeleted()) continue; 2187 2188 INDEX_3 i3; 2189 i3.I1() = sel.PNum(1); 2190 i3.I2() = sel.PNum(2); 2191 i3.I3() = sel.PNum(3); 2192 i3.Sort(); 2193 surfelementht -> PrepareSet (i3); 2194 } 2195 2196 surfelementht -> AllocateElements(); 2197 */ 2198 tn2se.Stop(); 2199 2200 tht.Start(); 2201 if (dimension==3) 2202 for (SurfaceElementIndex sei = 0; sei < GetNSE(); sei++) 2203 { 2204 const Element2d & sel = surfelements[sei]; 2205 if (sel.IsDeleted()) continue; 2206 2207 INDEX_3 i3; 2208 i3.I1() = sel.PNum(1); 2209 i3.I2() = sel.PNum(2); 2210 i3.I3() = sel.PNum(3); 2211 i3.Sort(); 2212 surfelementht -> Set (i3, sei); // war das wichtig ??? sel.GetIndex()); 2213 } 2214 tht.Stop(); 2215 2216 // int np = GetNP(); 2217 2218 if (dimension == 3) 2219 { 2220 static Timer t("Mesh::CalcSurfacesOfNode, pointloop"); RegionTimer reg (t); 2221 /* 2222 for (PointIndex pi = points.Begin(); pi < points.End(); pi++) 2223 points[pi].SetType (INNERPOINT); 2224 */ 2225 for (auto & p : points) 2226 p.SetType (INNERPOINT); 2227 2228 if (GetNFD() == 0) 2229 { 2230 for (SurfaceElementIndex sei = 0; sei < GetNSE(); sei++) 2231 { 2232 const Element2d & sel = surfelements[sei]; 2233 if (sel.IsDeleted()) continue; 2234 for (int j = 0; j < sel.GetNP(); j++) 2235 { 2236 PointIndex pi = SurfaceElement(sei)[j]; 2237 points[pi].SetType(FIXEDPOINT); 2238 } 2239 } 2240 } 2241 else 2242 { 2243 for (SurfaceElementIndex sei = 0; sei < GetNSE(); sei++) 2244 { 2245 const Element2d & sel = surfelements[sei]; 2246 if (sel.IsDeleted()) continue; 2247 for (int j = 0; j < sel.GetNP(); j++) 2248 { 2249 PointIndex pi = sel[j]; 2250 int ns = surfacesonnode[pi].Size(); 2251 if (ns == 1) 2252 points[pi].SetType(SURFACEPOINT); 2253 if (ns == 2) 2254 points[pi].SetType(EDGEPOINT); 2255 if (ns >= 3) 2256 points[pi].SetType(FIXEDPOINT); 2257 } 2258 } 2259 } 2260 } 2261 2262 /* 2263 for (int i = 0; i < segments.Size(); i++) 2264 { 2265 const Segment & seg = segments[i]; 2266 */ 2267 for (const Segment & seg : segments) 2268 { 2269 for (int j = 1; j <= 2; j++) 2270 { 2271 PointIndex hi = (j == 1) ? seg[0] : seg[1]; 2272 if (points[hi].Type() == INNERPOINT || 2273 points[hi].Type() == SURFACEPOINT) 2274 points[hi].SetType(EDGEPOINT); 2275 } 2276 } 2277 2278 for (int i = 0; i < lockedpoints.Size(); i++) 2279 points[lockedpoints[i]].SetType(FIXEDPOINT); 2280 2281 2282 /* 2283 for (i = 0; i < openelements.Size(); i++) 2284 { 2285 const Element2d & sel = openelements[i]; 2286 for (j = 0; j < sel.GetNP(); j++) 2287 { 2288 INDEX_2 i2; 2289 i2.I1() = sel.PNumMod(j+1); 2290 i2.I2() = sel.PNumMod(j+2); 2291 i2.Sort(); 2292 boundaryedges->Set (i2, 1); 2293 2294 points[sel[j]].SetType(FIXEDPOINT); 2295 } 2296 } 2297 */ 2298 2299 // eltyps.SetSize (GetNE()); 2300 // eltyps = FREEELEMENT; 2301 2302 for (int i = 0; i < GetNSeg(); i++) 2303 { 2304 const Segment & seg = segments[i]; 2305 INDEX_2 i2(seg[0], seg[1]); 2306 i2.Sort(); 2307 2308 //boundaryedges -> Set (i2, 2); 2309 segmentht -> Set (i2, i); 2310 } 2311 } 2312 2313 // NgBitArray base is PointIndex::BASE ... FixPoints(const NgBitArray & fixpoints)2314 void Mesh :: FixPoints (const NgBitArray & fixpoints) 2315 { 2316 if (fixpoints.Size() != GetNP()) 2317 { 2318 cerr << "Mesh::FixPoints: sizes don't fit" << endl; 2319 return; 2320 } 2321 int np = GetNP(); 2322 /* 2323 for (int i = 1; i <= np; i++) 2324 if (fixpoints.Test(i)) 2325 { 2326 points.Elem(i).SetType (FIXEDPOINT); 2327 } 2328 */ 2329 for (PointIndex pi : points.Range()) 2330 if (fixpoints.Test(pi)) 2331 points[pi].SetType(FIXEDPOINT); 2332 } 2333 2334 FindOpenElements(int dom)2335 void Mesh :: FindOpenElements (int dom) 2336 { 2337 static Timer t("Mesh::FindOpenElements"); RegionTimer reg (t); 2338 static Timer t_table("Mesh::FindOpenElements - build table"); 2339 static Timer t_pointloop("Mesh::FindOpenElements - pointloop"); 2340 2341 int np = GetNP(); 2342 int ne = GetNE(); 2343 int nse = GetNSE(); 2344 2345 t_table.Start(); 2346 2347 auto elsonpoint = ngcore::CreateSortedTable<ElementIndex, PointIndex>( volelements.Range(), 2348 [&](auto & table, ElementIndex ei) 2349 { 2350 const Element & el = (*this)[ei]; 2351 if (dom == 0 || dom == el.GetIndex()) 2352 { 2353 if (el.GetNP() == 4) 2354 { 2355 INDEX_4 i4(el[0], el[1], el[2], el[3]); 2356 i4.Sort(); 2357 table.Add (PointIndex(i4.I1()), ei); 2358 table.Add (PointIndex(i4.I2()), ei); 2359 } 2360 else 2361 { 2362 for (PointIndex pi : el.PNums()) 2363 table.Add(pi, ei); 2364 } 2365 } 2366 }, GetNP()); 2367 2368 2369 NgArray<int,PointIndex::BASE> numonpoint(np); 2370 /* 2371 numonpoint = 0; 2372 for (ElementIndex ei = 0; ei < ne; ei++) 2373 { 2374 const Element & el = (*this)[ei]; 2375 if (dom == 0 || dom == el.GetIndex()) 2376 { 2377 if (el.GetNP() == 4) 2378 { 2379 INDEX_4 i4(el[0], el[1], el[2], el[3]); 2380 i4.Sort(); 2381 numonpoint[i4.I1()]++; 2382 numonpoint[i4.I2()]++; 2383 } 2384 else 2385 for (int j = 0; j < el.GetNP(); j++) 2386 numonpoint[el[j]]++; 2387 } 2388 } 2389 2390 TABLE<ElementIndex,PointIndex::BASE> elsonpoint(numonpoint); 2391 for (ElementIndex ei = 0; ei < ne; ei++) 2392 { 2393 const Element & el = (*this)[ei]; 2394 if (dom == 0 || dom == el.GetIndex()) 2395 { 2396 if (el.GetNP() == 4) 2397 { 2398 INDEX_4 i4(el[0], el[1], el[2], el[3]); 2399 i4.Sort(); 2400 elsonpoint.Add (i4.I1(), ei); 2401 elsonpoint.Add (i4.I2(), ei); 2402 } 2403 else 2404 for (int j = 0; j < el.GetNP(); j++) 2405 elsonpoint.Add (el[j], ei); 2406 } 2407 } 2408 */ 2409 t_table.Stop(); 2410 2411 2412 NgArray<bool, 1> hasface(GetNFD()); 2413 2414 for (int i = 1; i <= GetNFD(); i++) 2415 { 2416 int domin = GetFaceDescriptor(i).DomainIn(); 2417 int domout = GetFaceDescriptor(i).DomainOut(); 2418 hasface[i] = 2419 ( dom == 0 && (domin != 0 || domout != 0) ) || 2420 ( dom != 0 && (domin == dom || domout == dom) ); 2421 } 2422 2423 numonpoint = 0; 2424 for (SurfaceElementIndex sii = 0; sii < nse; sii++) 2425 { 2426 int ind = surfelements[sii].GetIndex(); 2427 /* 2428 if ( 2429 GetFaceDescriptor(ind).DomainIn() && 2430 (dom == 0 || dom == GetFaceDescriptor(ind).DomainIn()) 2431 || 2432 GetFaceDescriptor(ind).DomainOut() && 2433 (dom == 0 || dom == GetFaceDescriptor(ind).DomainOut()) 2434 ) 2435 */ 2436 if (hasface[ind]) 2437 { 2438 /* 2439 Element2d hel = surfelements[i]; 2440 hel.NormalizeNumbering(); 2441 numonpoint[hel[0]]++; 2442 */ 2443 const Element2d & hel = surfelements[sii]; 2444 int mini = 0; 2445 for (int j = 1; j < hel.GetNP(); j++) 2446 if (hel[j] < hel[mini]) 2447 mini = j; 2448 numonpoint[hel[mini]]++; 2449 } 2450 } 2451 2452 TABLE<SurfaceElementIndex,PointIndex::BASE> selsonpoint(numonpoint); 2453 for (SurfaceElementIndex sii = 0; sii < nse; sii++) 2454 { 2455 int ind = surfelements[sii].GetIndex(); 2456 2457 /* 2458 if ( 2459 GetFaceDescriptor(ind).DomainIn() && 2460 (dom == 0 || dom == GetFaceDescriptor(ind).DomainIn()) 2461 || 2462 GetFaceDescriptor(ind).DomainOut() && 2463 (dom == 0 || dom == GetFaceDescriptor(ind).DomainOut()) 2464 ) 2465 */ 2466 if (hasface[ind]) 2467 { 2468 /* 2469 Element2d hel = surfelements[i]; 2470 hel.NormalizeNumbering(); 2471 selsonpoint.Add (hel[0], i); 2472 */ 2473 const Element2d & hel = surfelements[sii]; 2474 int mini = 0; 2475 for (int j = 1; j < hel.GetNP(); j++) 2476 if (hel[j] < hel[mini]) 2477 mini = j; 2478 selsonpoint.Add (hel[mini], sii); 2479 } 2480 } 2481 2482 2483 // PointIndex pi; 2484 // SurfaceElementIndex sei; 2485 // Element2d hel; 2486 2487 struct tval { int index; PointIndex p4; }; 2488 openelements.SetSize(0); 2489 2490 t_pointloop.Start(); 2491 2492 /* 2493 INDEX_3_CLOSED_HASHTABLE<tval> faceht(100); 2494 2495 for (PointIndex pi : points.Range()) 2496 if (selsonpoint[pi].Size()+elsonpoint[pi].Size()) 2497 { 2498 faceht.SetSize (2 * selsonpoint[pi].Size() + 4 * elsonpoint[pi].Size()); 2499 2500 for (SurfaceElementIndex sei : selsonpoint[pi]) 2501 { 2502 Element2d hel = SurfaceElement(sei); 2503 if (hel.GetType() == TRIG6) hel.SetType(TRIG); 2504 int ind = hel.GetIndex(); 2505 2506 if (GetFaceDescriptor(ind).DomainIn() && 2507 (dom == 0 || dom == GetFaceDescriptor(ind).DomainIn()) ) 2508 { 2509 hel.NormalizeNumbering(); 2510 if (hel.PNum(1) == pi) 2511 { 2512 INDEX_3 i3(hel[0], hel[1], hel[2]); 2513 tval i2; 2514 i2.index = GetFaceDescriptor(ind).DomainIn(); 2515 i2.p4 = (hel.GetNP() == 3) 2516 ? PointIndex (PointIndex::INVALID) 2517 : hel.PNum(4); 2518 faceht.Set (i3, i2); 2519 } 2520 } 2521 if (GetFaceDescriptor(ind).DomainOut() && 2522 (dom == 0 || dom == GetFaceDescriptor(ind).DomainOut()) ) 2523 { 2524 hel.Invert(); 2525 hel.NormalizeNumbering(); 2526 if (hel.PNum(1) == pi) 2527 { 2528 INDEX_3 i3(hel[0], hel[1], hel[2]); 2529 tval i2; 2530 i2.index = GetFaceDescriptor(ind).DomainOut(); 2531 i2.p4 = (hel.GetNP() == 3) 2532 ? PointIndex (PointIndex::INVALID) 2533 : hel.PNum(4); 2534 faceht.Set (i3, i2); 2535 } 2536 } 2537 } 2538 2539 for (ElementIndex ei : elsonpoint[pi]) 2540 { 2541 const Element & el = VolumeElement(ei); 2542 2543 if (dom == 0 || el.GetIndex() == dom) 2544 { 2545 for (int j = 1; j <= el.GetNFaces(); j++) 2546 { 2547 Element2d hel(TRIG); 2548 el.GetFace (j, hel); 2549 hel.Invert(); 2550 hel.NormalizeNumbering(); 2551 2552 if (hel[0] == pi) 2553 { 2554 INDEX_3 i3(hel[0], hel[1], hel[2]); 2555 2556 if (faceht.Used (i3)) 2557 { 2558 tval i2 = faceht.Get(i3); 2559 if (i2.index == el.GetIndex()) 2560 { 2561 i2.index = PointIndex::BASE-1; 2562 faceht.Set (i3, i2); 2563 } 2564 else 2565 { 2566 if (i2.index == 0) 2567 { 2568 PrintSysError ("more elements on face"); 2569 (*testout) << "more elements on face!!!" << endl; 2570 (*testout) << "el = " << el << endl; 2571 (*testout) << "hel = " << hel << endl; 2572 (*testout) << "face = " << i3 << endl; 2573 (*testout) << "points = " << endl; 2574 for (int jj = 1; jj <= 3; jj++) 2575 (*testout) << "p = " << Point(i3.I(jj)) << endl; 2576 } 2577 } 2578 } 2579 else 2580 { 2581 hel.Invert(); 2582 hel.NormalizeNumbering(); 2583 INDEX_3 i3(hel[0], hel[1], hel[2]); 2584 2585 tval i2; 2586 i2.index = el.GetIndex(); 2587 i2.p4 = (hel.GetNP() == 3) 2588 ? PointIndex (PointIndex::INVALID) 2589 : hel[3]; 2590 faceht.Set (i3, i2); 2591 } 2592 } 2593 } 2594 } 2595 } 2596 2597 for (int i = 0; i < faceht.Size(); i++) 2598 if (faceht.UsedPos (i)) 2599 { 2600 INDEX_3 i3; 2601 //INDEX_2 i2; 2602 tval i2; 2603 faceht.GetData (i, i3, i2); 2604 if (i2.index != PointIndex::BASE-1) 2605 { 2606 Element2d tri ( (i2.p4 == PointIndex::BASE-1) ? TRIG : QUAD); 2607 for (int l = 0; l < 3; l++) 2608 tri[l] = i3.I(l+1); 2609 tri.PNum(4) = i2.p4; 2610 tri.SetIndex (i2.index); 2611 openelements.Append (tri); 2612 } 2613 } 2614 } 2615 2616 */ 2617 2618 size_t numtasks = 4*ngcore::TaskManager::GetNumThreads(); 2619 Array<Array<Element2d>> thread_openelements(numtasks); 2620 ParallelJob 2621 ( [&](TaskInfo & ti) 2622 { 2623 auto myrange = points.Range().Split(ti.task_nr, ti.ntasks); 2624 INDEX_3_CLOSED_HASHTABLE<tval> faceht(100); 2625 for (PointIndex pi : myrange) 2626 if (selsonpoint[pi].Size()+elsonpoint[pi].Size()) 2627 { 2628 faceht.SetSize (2 * selsonpoint[pi].Size() + 4 * elsonpoint[pi].Size()); 2629 2630 for (SurfaceElementIndex sei : selsonpoint[pi]) 2631 { 2632 Element2d hel = SurfaceElement(sei); 2633 if (hel.GetType() == TRIG6) hel.SetType(TRIG); 2634 int ind = hel.GetIndex(); 2635 2636 if (GetFaceDescriptor(ind).DomainIn() && 2637 (dom == 0 || dom == GetFaceDescriptor(ind).DomainIn()) ) 2638 { 2639 hel.NormalizeNumbering(); 2640 if (hel.PNum(1) == pi) 2641 { 2642 INDEX_3 i3(hel[0], hel[1], hel[2]); 2643 tval i2; 2644 i2.index = GetFaceDescriptor(ind).DomainIn(); 2645 i2.p4 = (hel.GetNP() == 3) 2646 ? PointIndex (PointIndex::INVALID) 2647 : hel.PNum(4); 2648 faceht.Set (i3, i2); 2649 } 2650 } 2651 if (GetFaceDescriptor(ind).DomainOut() && 2652 (dom == 0 || dom == GetFaceDescriptor(ind).DomainOut()) ) 2653 { 2654 hel.Invert(); 2655 hel.NormalizeNumbering(); 2656 if (hel.PNum(1) == pi) 2657 { 2658 INDEX_3 i3(hel[0], hel[1], hel[2]); 2659 tval i2; 2660 i2.index = GetFaceDescriptor(ind).DomainOut(); 2661 i2.p4 = (hel.GetNP() == 3) 2662 ? PointIndex (PointIndex::INVALID) 2663 : hel.PNum(4); 2664 faceht.Set (i3, i2); 2665 } 2666 } 2667 } 2668 2669 for (ElementIndex ei : elsonpoint[pi]) 2670 { 2671 const Element & el = VolumeElement(ei); 2672 2673 if (dom == 0 || el.GetIndex() == dom) 2674 { 2675 for (int j = 1; j <= el.GetNFaces(); j++) 2676 { 2677 Element2d hel(TRIG); 2678 el.GetFace (j, hel); 2679 hel.Invert(); 2680 hel.NormalizeNumbering(); 2681 2682 if (hel[0] == pi) 2683 { 2684 INDEX_3 i3(hel[0], hel[1], hel[2]); 2685 2686 if (faceht.Used (i3)) 2687 { 2688 tval i2 = faceht.Get(i3); 2689 if (i2.index == el.GetIndex()) 2690 { 2691 i2.index = PointIndex::BASE-1; 2692 faceht.Set (i3, i2); 2693 } 2694 else 2695 { 2696 if (i2.index == 0) 2697 { 2698 PrintSysError ("more elements on face"); 2699 (*testout) << "more elements on face!!!" << endl; 2700 (*testout) << "el = " << el << endl; 2701 (*testout) << "hel = " << hel << endl; 2702 (*testout) << "face = " << i3 << endl; 2703 (*testout) << "points = " << endl; 2704 for (int jj = 1; jj <= 3; jj++) 2705 (*testout) << "p = " << Point(i3.I(jj)) << endl; 2706 } 2707 } 2708 } 2709 else 2710 { 2711 hel.Invert(); 2712 hel.NormalizeNumbering(); 2713 INDEX_3 i3(hel[0], hel[1], hel[2]); 2714 2715 tval i2; 2716 i2.index = el.GetIndex(); 2717 i2.p4 = (hel.GetNP() == 3) 2718 ? PointIndex (PointIndex::INVALID) 2719 : hel[3]; 2720 faceht.Set (i3, i2); 2721 } 2722 } 2723 } 2724 } 2725 } 2726 2727 for (int i = 0; i < faceht.Size(); i++) 2728 if (faceht.UsedPos (i)) 2729 { 2730 INDEX_3 i3; 2731 tval i2; 2732 faceht.GetData (i, i3, i2); 2733 if (i2.index != PointIndex::BASE-1) 2734 { 2735 Element2d tri ( (i2.p4 == PointIndex::BASE-1) ? TRIG : QUAD); 2736 for (int l = 0; l < 3; l++) 2737 tri[l] = i3.I(l+1); 2738 tri.PNum(4) = i2.p4; 2739 tri.SetIndex (i2.index); 2740 thread_openelements[ti.task_nr].Append (tri); 2741 } 2742 } 2743 }}, numtasks); 2744 2745 for (auto & a : thread_openelements) 2746 for (auto & el : a) 2747 openelements.Append (el); 2748 2749 t_pointloop.Stop(); 2750 2751 int cnt3 = 0; 2752 for (int i = 0; i < openelements.Size(); i++) 2753 if (openelements[i].GetNP() == 3) 2754 cnt3++; 2755 2756 int cnt4 = openelements.Size() - cnt3; 2757 2758 2759 MyStr treequad; 2760 if (cnt4) 2761 treequad = MyStr(" (") + MyStr(cnt3) + MyStr (" + ") + 2762 MyStr(cnt4) + MyStr(")"); 2763 2764 PrintMessage (5, openelements.Size(), treequad, " open elements"); 2765 2766 BuildBoundaryEdges(); 2767 2768 2769 for (int i = 1; i <= openelements.Size(); i++) 2770 { 2771 const Element2d & sel = openelements.Get(i); 2772 2773 if (boundaryedges) 2774 for (int j = 1; j <= sel.GetNP(); j++) 2775 { 2776 INDEX_2 i2; 2777 i2.I1() = sel.PNumMod(j); 2778 i2.I2() = sel.PNumMod(j+1); 2779 i2.Sort(); 2780 boundaryedges->Set (i2, 1); 2781 } 2782 2783 for (int j = 1; j <= 3; j++) 2784 { 2785 PointIndex pi = sel.PNum(j); 2786 // if (pi < points.End()) 2787 if (pi < *points.Range().end()) 2788 points[pi].SetType (FIXEDPOINT); 2789 } 2790 } 2791 2792 2793 2794 /* 2795 for (i = 1; i <= GetNSeg(); i++) 2796 { 2797 const Segment & seg = LineSegment(i); 2798 INDEX_2 i2(seg[0], seg[1]); 2799 i2.Sort(); 2800 2801 if (!boundaryedges->Used (i2)) 2802 cerr << "WARNING: no boundedge, but seg edge: " << i2 << endl; 2803 2804 boundaryedges -> Set (i2, 2); 2805 segmentht -> Set (i2, i-1); 2806 } 2807 */ 2808 } 2809 HasOpenQuads() const2810 bool Mesh :: HasOpenQuads () const 2811 { 2812 int no = GetNOpenElements(); 2813 for (int i = 0; i < no; i++) 2814 if (openelements[i].GetNP() == 4) 2815 return true; 2816 return false; 2817 } 2818 2819 2820 2821 2822 FindOpenSegments(int surfnr)2823 void Mesh :: FindOpenSegments (int surfnr) 2824 { 2825 // int i, j, k; 2826 2827 // new version, general elements 2828 // hash index: pnum1-2, surfnr 2829 // hash data : surfel-nr (pos) or segment nr(neg) 2830 INDEX_3_HASHTABLE<int> faceht(4 * GetNSE()+GetNSeg()+1); 2831 2832 PrintMessage (5, "Test Opensegments"); 2833 for (int i = 1; i <= GetNSeg(); i++) 2834 { 2835 const Segment & seg = LineSegment (i); 2836 2837 if (surfnr == 0 || seg.si == surfnr) 2838 { 2839 INDEX_3 key(seg[0], seg[1], seg.si); 2840 int data = -i; 2841 2842 if (faceht.Used (key)) 2843 { 2844 cerr << "ERROR: Segment " << seg << " already used" << endl; 2845 (*testout) << "ERROR: Segment " << seg << " already used" << endl; 2846 } 2847 2848 faceht.Set (key, data); 2849 } 2850 } 2851 2852 2853 /* 2854 // not possible with surfnr as hash-index 2855 for (int i = 1; i <= GetNSeg(); i++) 2856 { 2857 const Segment & seg = LineSegment (i); 2858 2859 if (surfnr == 0 || seg.si == surfnr) 2860 { 2861 INDEX_2 key(seg[1], seg[0]); 2862 if (!faceht.Used(key)) 2863 { 2864 cerr << "ERROR: Segment " << seg << " brother not used" << endl; 2865 (*testout) << "ERROR: Segment " << seg << " brother not used" << endl; 2866 } 2867 } 2868 } 2869 */ 2870 2871 2872 // bool buggy = false; 2873 // ofstream bout("buggy.out"); 2874 2875 for (int i = 1; i <= GetNSE(); i++) 2876 { 2877 const Element2d & el = SurfaceElement(i); 2878 if (el.IsDeleted()) continue; 2879 2880 if (surfnr == 0 || el.GetIndex() == surfnr) 2881 { 2882 for (int j = 1; j <= el.GetNP(); j++) 2883 { 2884 INDEX_3 seg (el.PNumMod(j), el.PNumMod(j+1), el.GetIndex()); 2885 int data; 2886 2887 if (seg.I1() < PointIndex::BASE || seg.I2() < PointIndex::BASE) 2888 cerr << "seg = " << seg << endl; 2889 2890 if (faceht.Used(seg)) 2891 { 2892 faceht.Set (seg, 0); 2893 /* 2894 data = faceht.Get(seg); 2895 2896 if (data.I1() == el.GetIndex()) 2897 { 2898 data.I1() = 0; 2899 faceht.Set (seg, data); 2900 } 2901 else 2902 { 2903 // buggy = true; 2904 PrintWarning ("hash table si not fitting for segment: ", 2905 seg.I1(), "-", seg.I2(), " other = ", 2906 data.I2(), ", surfnr = ", surfnr); 2907 } 2908 */ 2909 } 2910 else 2911 { 2912 Swap (seg.I1(), seg.I2()); 2913 // data.I1() = el.GetIndex(); 2914 // data.I2() = i; 2915 faceht.Set (seg, i); 2916 } 2917 } 2918 } 2919 } 2920 2921 /* 2922 if (buggy) 2923 { 2924 for (int i = 1; i <= GetNSeg(); i++) 2925 bout << "seg" << i << " " << LineSegment(i) << endl; 2926 2927 for (int i = 1; i <= GetNSE(); i++) 2928 bout << "sel" << i << " " << SurfaceElement(i) << " ind = " 2929 << SurfaceElement(i).GetIndex() << endl; 2930 2931 bout << "hashtable: " << endl; 2932 for (int j = 1; j <= faceht.GetNBags(); j++) 2933 { 2934 bout << "bag " << j << ":" << endl; 2935 for (int k = 1; k <= faceht.GetBagSize(j); k++) 2936 { 2937 INDEX_2 i2, data; 2938 faceht.GetData (j, k, i2, data); 2939 bout << "key = " << i2 << ", data = " << data << endl; 2940 } 2941 } 2942 exit(1); 2943 } 2944 */ 2945 2946 (*testout) << "open segments: " << endl; 2947 opensegments.SetSize(0); 2948 for (int i = 1; i <= faceht.GetNBags(); i++) 2949 for (int j = 1; j <= faceht.GetBagSize(i); j++) 2950 { 2951 INDEX_3 i2; 2952 int data; 2953 faceht.GetData (i, j, i2, data); 2954 if (data) // surfnr 2955 { 2956 Segment seg; 2957 seg[0] = i2.I1(); 2958 seg[1] = i2.I2(); 2959 seg.si = i2.I3(); 2960 2961 // find geomdata: 2962 if (data > 0) 2963 { 2964 // segment due to triangle 2965 const Element2d & el = SurfaceElement (data); 2966 for (int k = 1; k <= el.GetNP(); k++) 2967 { 2968 if (seg[0] == el.PNum(k)) 2969 seg.geominfo[0] = el.GeomInfoPi(k); 2970 if (seg[1] == el.PNum(k)) 2971 seg.geominfo[1] = el.GeomInfoPi(k); 2972 } 2973 2974 (*testout) << "trig seg: "; 2975 } 2976 else 2977 { 2978 // segment due to line 2979 const Segment & lseg = LineSegment (-data); 2980 seg.geominfo[0] = lseg.geominfo[0]; 2981 seg.geominfo[1] = lseg.geominfo[1]; 2982 2983 (*testout) << "line seg: "; 2984 } 2985 2986 (*testout) << seg[0] << " - " << seg[1] 2987 << " len = " << Dist (Point(seg[0]), Point(seg[1])) 2988 << endl; 2989 2990 opensegments.Append (seg); 2991 if (seg.geominfo[0].trignum <= 0 || seg.geominfo[1].trignum <= 0) 2992 { 2993 (*testout) << "Problem with open segment: " << seg << endl; 2994 } 2995 2996 } 2997 } 2998 2999 PrintMessage (3, opensegments.Size(), " open segments found"); 3000 (*testout) << opensegments.Size() << " open segments found" << endl; 3001 3002 /* 3003 ptyps.SetSize (GetNP()); 3004 for (i = 1; i <= ptyps.Size(); i++) 3005 ptyps.Elem(i) = SURFACEPOINT; 3006 3007 for (i = 1; i <= GetNSeg(); i++) 3008 { 3009 const Segment & seg = LineSegment (i); 3010 ptyps.Elem(seg[0]) = EDGEPOINT; 3011 ptyps.Elem(seg[1]) = EDGEPOINT; 3012 } 3013 for (i = 1; i <= GetNOpenSegments(); i++) 3014 { 3015 const Segment & seg = GetOpenSegment (i); 3016 ptyps.Elem(seg[0]) = EDGEPOINT; 3017 ptyps.Elem(seg[1]) = EDGEPOINT; 3018 } 3019 */ 3020 /* 3021 for (int i = 1; i <= points.Size(); i++) 3022 points.Elem(i).SetType(SURFACEPOINT); 3023 */ 3024 for (auto & p : points) 3025 p.SetType (SURFACEPOINT); 3026 3027 for (int i = 1; i <= GetNSeg(); i++) 3028 { 3029 const Segment & seg = LineSegment (i); 3030 points[seg[0]].SetType(EDGEPOINT); 3031 points[seg[1]].SetType(EDGEPOINT); 3032 } 3033 for (int i = 1; i <= GetNOpenSegments(); i++) 3034 { 3035 const Segment & seg = GetOpenSegment (i); 3036 points[seg[0]].SetType (EDGEPOINT); 3037 points[seg[1]].SetType (EDGEPOINT); 3038 } 3039 3040 3041 3042 /* 3043 3044 for (i = 1; i <= openelements.Size(); i++) 3045 { 3046 const Element2d & sel = openelements.Get(i); 3047 3048 if (boundaryedges) 3049 for (j = 1; j <= sel.GetNP(); j++) 3050 { 3051 INDEX_2 i2; 3052 i2.I1() = sel.PNumMod(j); 3053 i2.I2() = sel.PNumMod(j+1); 3054 i2.Sort(); 3055 boundaryedges->Set (i2, 1); 3056 } 3057 3058 for (j = 1; j <= 3; j++) 3059 { 3060 int pi = sel.PNum(j); 3061 if (pi <= ptyps.Size()) 3062 ptyps.Elem(pi) = FIXEDPOINT; 3063 } 3064 } 3065 */ 3066 } 3067 3068 RemoveOneLayerSurfaceElements()3069 void Mesh :: RemoveOneLayerSurfaceElements () 3070 { 3071 int np = GetNP(); 3072 3073 FindOpenSegments(); 3074 NgBitArray frontpoints(np+1); // for 0- and 1-based 3075 frontpoints.Clear(); 3076 3077 for (int i = 1; i <= GetNOpenSegments(); i++) 3078 { 3079 const Segment & seg = GetOpenSegment(i); 3080 frontpoints.Set (seg[0]); 3081 frontpoints.Set (seg[1]); 3082 } 3083 3084 for (int i = 1; i <= GetNSE(); i++) 3085 { 3086 Element2d & sel = surfelements[i-1]; 3087 bool remove = false; 3088 for (int j = 1; j <= sel.GetNP(); j++) 3089 if (frontpoints.Test(sel.PNum(j))) 3090 remove = true; 3091 if (remove) 3092 sel.PNum(1).Invalidate(); 3093 } 3094 3095 for (int i = surfelements.Size(); i >= 1; i--) 3096 { 3097 if (!surfelements[i-1].PNum(1).IsValid()) 3098 { 3099 surfelements[i-1] = surfelements.Last(); 3100 surfelements.DeleteLast(); 3101 } 3102 } 3103 3104 RebuildSurfaceElementLists (); 3105 /* 3106 for (int i = 0; i < facedecoding.Size(); i++) 3107 facedecoding[i].firstelement = -1; 3108 for (int i = surfelements.Size()-1; i >= 0; i--) 3109 { 3110 int ind = surfelements[i].GetIndex(); 3111 surfelements[i].next = facedecoding[ind-1].firstelement; 3112 facedecoding[ind-1].firstelement = i; 3113 } 3114 */ 3115 3116 timestamp = NextTimeStamp(); 3117 // Compress(); 3118 } 3119 3120 3121 3122 3123 FreeOpenElementsEnvironment(int layers)3124 void Mesh :: FreeOpenElementsEnvironment (int layers) 3125 { 3126 static Timer timer("FreeOpenElementsEnvironment"); RegionTimer rt(timer); 3127 int i, j, k; 3128 PointIndex pi; 3129 const int large = 9999; 3130 NgArray<int,PointIndex::BASE> dist(GetNP()); 3131 3132 dist = large; 3133 3134 for (int i = 1; i <= GetNOpenElements(); i++) 3135 { 3136 const Element2d & face = OpenElement(i); 3137 for (j = 0; j < face.GetNP(); j++) 3138 dist[face[j]] = 1; 3139 } 3140 3141 for (k = 1; k <= layers; k++) 3142 for (i = 1; i <= GetNE(); i++) 3143 { 3144 const Element & el = VolumeElement(i); 3145 if (el[0] == -1 || el.IsDeleted()) continue; 3146 3147 int elmin = large; 3148 for (j = 0; j < el.GetNP(); j++) 3149 if (dist[el[j]] < elmin) 3150 elmin = dist[el[j]]; 3151 3152 if (elmin < large) 3153 { 3154 for (j = 0; j < el.GetNP(); j++) 3155 if (dist[el[j]] > elmin+1) 3156 dist[el[j]] = elmin+1; 3157 } 3158 } 3159 3160 int cntfree = 0; 3161 for (i = 1; i <= GetNE(); i++) 3162 { 3163 Element & el = VolumeElement(i); 3164 if (el[0] == -1 || el.IsDeleted()) continue; 3165 3166 int elmin = large; 3167 for (j = 0; j < el.GetNP(); j++) 3168 if (dist[el[j]] < elmin) 3169 elmin = dist[el[j]]; 3170 3171 el.flags.fixed = elmin > layers; 3172 // eltyps.Elem(i) = (elmin <= layers) ? 3173 // FREEELEMENT : FIXEDELEMENT; 3174 if (elmin <= layers) 3175 cntfree++; 3176 } 3177 3178 PrintMessage (5, "free: ", cntfree, ", fixed: ", GetNE()-cntfree); 3179 (*testout) << "free: " << cntfree << ", fixed: " << GetNE()-cntfree << endl; 3180 3181 for (pi = PointIndex::BASE; 3182 pi < GetNP()+PointIndex::BASE; pi++) 3183 { 3184 if (dist[pi] > layers+1) 3185 points[pi].SetType(FIXEDPOINT); 3186 } 3187 } 3188 3189 3190 SetLocalH(netgen::Point<3> pmin,netgen::Point<3> pmax,double grading)3191 void Mesh :: SetLocalH (netgen::Point<3> pmin, netgen::Point<3> pmax, double grading) 3192 { 3193 using netgen::Point; 3194 Point<3> c = Center (pmin, pmax); 3195 double d = max3 (pmax(0)-pmin(0), 3196 pmax(1)-pmin(1), 3197 pmax(2)-pmin(2)); 3198 d /= 2; 3199 Point<3> pmin2 = c - Vec<3> (d, d, d); 3200 Point<3> pmax2 = c + Vec<3> (d, d, d); 3201 3202 lochfunc = make_unique<LocalH> (pmin2, pmax2, grading, dimension); 3203 } 3204 RestrictLocalH(const Point3d & p,double hloc)3205 void Mesh :: RestrictLocalH (const Point3d & p, double hloc) 3206 { 3207 if(hloc < hmin) 3208 hloc = hmin; 3209 3210 //cout << "restrict h in " << p << " to " << hloc << endl; 3211 if (!lochfunc) 3212 { 3213 PrintWarning("RestrictLocalH called, creating mesh-size tree"); 3214 3215 Point3d boxmin, boxmax; 3216 GetBox (boxmin, boxmax); 3217 SetLocalH (boxmin, boxmax, 0.8); 3218 } 3219 3220 lochfunc -> SetH (p, hloc); 3221 } 3222 RestrictLocalHLine(const Point3d & p1,const Point3d & p2,double hloc)3223 void Mesh :: RestrictLocalHLine (const Point3d & p1, 3224 const Point3d & p2, 3225 double hloc) 3226 { 3227 if(hloc < hmin) 3228 hloc = hmin; 3229 3230 // cout << "restrict h along " << p1 << " - " << p2 << " to " << hloc << endl; 3231 int i; 3232 int steps = int (Dist (p1, p2) / hloc) + 2; 3233 Vec3d v(p1, p2); 3234 3235 for (i = 0; i <= steps; i++) 3236 { 3237 Point3d p = p1 + (double(i)/double(steps) * v); 3238 RestrictLocalH (p, hloc); 3239 } 3240 } 3241 3242 SetMinimalH(double h)3243 void Mesh :: SetMinimalH (double h) 3244 { 3245 hmin = h; 3246 } 3247 3248 SetGlobalH(double h)3249 void Mesh :: SetGlobalH (double h) 3250 { 3251 hglob = h; 3252 } 3253 MaxHDomain(int dom) const3254 double Mesh :: MaxHDomain (int dom) const 3255 { 3256 if (maxhdomain.Size()) 3257 return maxhdomain.Get(dom); 3258 else 3259 return 1e10; 3260 } 3261 SetMaxHDomain(const NgArray<double> & mhd)3262 void Mesh :: SetMaxHDomain (const NgArray<double> & mhd) 3263 { 3264 maxhdomain.SetSize(mhd.Size()); 3265 for (int i = 1; i <= mhd.Size(); i++) 3266 maxhdomain.Elem(i) = mhd.Get(i); 3267 } 3268 3269 GetH(const Point3d & p) const3270 double Mesh :: GetH (const Point3d & p) const 3271 { 3272 double hmin = hglob; 3273 if (lochfunc) 3274 { 3275 double hl = lochfunc->GetH (p); 3276 if (hl < hglob) 3277 hmin = hl; 3278 } 3279 return hmin; 3280 } 3281 GetMinH(const Point3d & pmin,const Point3d & pmax)3282 double Mesh :: GetMinH (const Point3d & pmin, const Point3d & pmax) 3283 { 3284 double hmin = hglob; 3285 if (lochfunc) 3286 { 3287 double hl = lochfunc->GetMinH (pmin, pmax); 3288 if (hl < hmin) 3289 hmin = hl; 3290 } 3291 return hmin; 3292 } 3293 3294 3295 3296 3297 AverageH(int surfnr) const3298 double Mesh :: AverageH (int surfnr) const 3299 { 3300 int i, j, n; 3301 double hi, hsum; 3302 double maxh = 0, minh = 1e10; 3303 3304 hsum = 0; 3305 n = 0; 3306 for (i = 1; i <= GetNSE(); i++) 3307 { 3308 const Element2d & el = SurfaceElement(i); 3309 if (surfnr == 0 || el.GetIndex() == surfnr) 3310 { 3311 for (j = 1; j <= 3; j++) 3312 { 3313 hi = Dist (Point (el.PNumMod(j)), 3314 Point (el.PNumMod(j+1))); 3315 3316 hsum += hi; 3317 3318 if (hi > maxh) maxh = hi; 3319 if (hi < minh) minh = hi; 3320 n++; 3321 } 3322 } 3323 } 3324 3325 PrintMessage (5, "minh = ", minh, " avh = ", (hsum/n), " maxh = ", maxh); 3326 return (hsum / n); 3327 } 3328 3329 3330 CalcLocalH(double grading)3331 void Mesh :: CalcLocalH (double grading) 3332 { 3333 static Timer t("Mesh::CalcLocalH"); RegionTimer reg(t); 3334 3335 if (!lochfunc) 3336 { 3337 Point3d pmin, pmax; 3338 GetBox (pmin, pmax); 3339 // SetLocalH (pmin, pmax, mparam.grading); 3340 SetLocalH (pmin, pmax, grading); 3341 } 3342 3343 PrintMessage (3, 3344 "CalcLocalH: ", 3345 GetNP(), " Points ", 3346 GetNE(), " Elements ", 3347 GetNSE(), " Surface Elements"); 3348 3349 3350 for (int i = 0; i < GetNSE(); i++) 3351 { 3352 const Element2d & el = surfelements[i]; 3353 int j; 3354 3355 if (el.GetNP() == 3) 3356 { 3357 double hel = -1; 3358 for (j = 1; j <= 3; j++) 3359 { 3360 const Point3d & p1 = points[el.PNumMod(j)]; 3361 const Point3d & p2 = points[el.PNumMod(j+1)]; 3362 3363 /* 3364 INDEX_2 i21(el.PNumMod(j), el.PNumMod(j+1)); 3365 INDEX_2 i22(el.PNumMod(j+1), el.PNumMod(j)); 3366 if (! identifiedpoints->Used (i21) && 3367 ! identifiedpoints->Used (i22) ) 3368 */ 3369 if (!ident -> UsedSymmetric (el.PNumMod(j), 3370 el.PNumMod(j+1))) 3371 { 3372 double hedge = Dist (p1, p2); 3373 if (hedge > hel) 3374 hel = hedge; 3375 // lochfunc->SetH (Center (p1, p2), 2 * Dist (p1, p2)); 3376 // (*testout) << "trigseth, p1,2 = " << el.PNumMod(j) << ", " << el.PNumMod(j+1) 3377 // << " h = " << (2 * Dist(p1, p2)) << endl; 3378 } 3379 } 3380 3381 if (hel > 0) 3382 { 3383 const Point3d & p1 = points[el.PNum(1)]; 3384 const Point3d & p2 = points[el.PNum(2)]; 3385 const Point3d & p3 = points[el.PNum(3)]; 3386 lochfunc->SetH (Center (p1, p2, p3), hel); 3387 } 3388 } 3389 else 3390 { 3391 { 3392 const Point3d & p1 = points[el.PNum(1)]; 3393 const Point3d & p2 = points[el.PNum(2)]; 3394 lochfunc->SetH (Center (p1, p2), 2 * Dist (p1, p2)); 3395 } 3396 { 3397 const Point3d & p1 = points[el.PNum(3)]; 3398 const Point3d & p2 = points[el.PNum(4)]; 3399 lochfunc->SetH (Center (p1, p2), 2 * Dist (p1, p2)); 3400 } 3401 } 3402 } 3403 3404 for (int i = 0; i < GetNSeg(); i++) 3405 { 3406 const Segment & seg = segments[i]; 3407 const Point3d & p1 = points[seg[0]]; 3408 const Point3d & p2 = points[seg[1]]; 3409 /* 3410 INDEX_2 i21(seg[0], seg[1]); 3411 INDEX_2 i22(seg[1], seg[0]); 3412 if (identifiedpoints) 3413 if (!identifiedpoints->Used (i21) && !identifiedpoints->Used (i22)) 3414 */ 3415 if (!ident -> UsedSymmetric (seg[0], seg[1])) 3416 { 3417 lochfunc->SetH (Center (p1, p2), Dist (p1, p2)); 3418 } 3419 } 3420 /* 3421 cerr << "do vol" << endl; 3422 for (i = 1; i <= GetNE(); i++) 3423 { 3424 const Element & el = VolumeElement(i); 3425 if (el.GetType() == TET) 3426 { 3427 int j, k; 3428 for (j = 2; j <= 4; j++) 3429 for (k = 1; k < j; k++) 3430 { 3431 const Point3d & p1 = Point (el.PNum(j)); 3432 const Point3d & p2 = Point (el.PNum(k)); 3433 lochfunc->SetH (Center (p1, p2), 2 * Dist (p1, p2)); 3434 (*testout) << "set vol h to " << (2 * Dist (p1, p2)) << endl; 3435 3436 } 3437 } 3438 } 3439 */ 3440 3441 /* 3442 const char * meshsizefilename = 3443 globflags.GetStringFlag ("meshsize", NULL); 3444 if (meshsizefilename) 3445 { 3446 ifstream msf(meshsizefilename); 3447 if (msf) 3448 { 3449 int nmsp; 3450 msf >> nmsp; 3451 for (i = 1; i <= nmsp; i++) 3452 { 3453 Point3d pi; 3454 double hi; 3455 msf >> pi.X() >> pi.Y() >> pi.Z(); 3456 msf >> hi; 3457 lochfunc->SetH (pi, hi); 3458 } 3459 } 3460 } 3461 */ 3462 // lochfunc -> Convexify(); 3463 // lochfunc -> PrintMemInfo (cout); 3464 } 3465 3466 CalcLocalHFromPointDistances(double grading)3467 void Mesh :: CalcLocalHFromPointDistances(double grading) 3468 { 3469 PrintMessage (3, "Calculating local h from point distances"); 3470 3471 if (!lochfunc) 3472 { 3473 Point3d pmin, pmax; 3474 GetBox (pmin, pmax); 3475 3476 // SetLocalH (pmin, pmax, mparam.grading); 3477 SetLocalH (pmin, pmax, grading); 3478 } 3479 3480 PointIndex i,j; 3481 double hl; 3482 3483 3484 for (i = PointIndex::BASE; 3485 i < GetNP()+PointIndex::BASE; i++) 3486 { 3487 for(j=i+1; j<GetNP()+PointIndex::BASE; j++) 3488 { 3489 const Point3d & p1 = points[i]; 3490 const Point3d & p2 = points[j]; 3491 hl = Dist(p1,p2); 3492 RestrictLocalH(p1,hl); 3493 RestrictLocalH(p2,hl); 3494 //cout << "restricted h at " << p1 << " and " << p2 << " to " << hl << endl; 3495 } 3496 } 3497 3498 3499 } 3500 3501 CalcLocalHFromSurfaceCurvature(double grading,double elperr)3502 void Mesh :: CalcLocalHFromSurfaceCurvature (double grading, double elperr) 3503 { 3504 PrintMessage (3, "Calculating local h from surface curvature"); 3505 3506 if (!lochfunc) 3507 { 3508 Point3d pmin, pmax; 3509 GetBox (pmin, pmax); 3510 3511 // SetLocalH (pmin, pmax, mparam.grading); 3512 SetLocalH (pmin, pmax, grading); 3513 } 3514 3515 3516 INDEX_2_HASHTABLE<int> edges(3 * GetNP() + 2); 3517 INDEX_2_HASHTABLE<int> bedges(GetNSeg() + 2); 3518 int i, j; 3519 3520 for (i = 1; i <= GetNSeg(); i++) 3521 { 3522 const Segment & seg = LineSegment(i); 3523 INDEX_2 i2(seg[0], seg[1]); 3524 i2.Sort(); 3525 bedges.Set (i2, 1); 3526 } 3527 for (i = 1; i <= GetNSE(); i++) 3528 { 3529 const Element2d & sel = SurfaceElement(i); 3530 if (!sel.PNum(1)) 3531 continue; 3532 for (j = 1; j <= 3; j++) 3533 { 3534 INDEX_2 i2(sel.PNumMod(j), sel.PNumMod(j+1)); 3535 i2.Sort(); 3536 if (bedges.Used(i2)) continue; 3537 3538 if (edges.Used(i2)) 3539 { 3540 int other = edges.Get(i2); 3541 3542 const Element2d & elother = SurfaceElement(other); 3543 3544 int pi3 = 1; 3545 while ( (sel.PNum(pi3) == i2.I1()) || 3546 (sel.PNum(pi3) == i2.I2())) 3547 pi3++; 3548 pi3 = sel.PNum(pi3); 3549 3550 int pi4 = 1; 3551 while ( (elother.PNum(pi4) == i2.I1()) || 3552 (elother.PNum(pi4) == i2.I2())) 3553 pi4++; 3554 pi4 = elother.PNum(pi4); 3555 3556 double rad = ComputeCylinderRadius (Point (PointIndex(i2.I1())), 3557 Point (PointIndex(i2.I2())), 3558 Point (PointIndex(pi3)), 3559 Point (PointIndex(pi4))); 3560 3561 RestrictLocalHLine (Point(PointIndex(i2.I1())), Point(PointIndex(i2.I2())), rad/elperr); 3562 3563 3564 /* 3565 (*testout) << "pi1,2, 3, 4 = " << i2.I1() << ", " << i2.I2() << ", " << pi3 << ", " << pi4 3566 << " p1 = " << Point(i2.I1()) 3567 << ", p2 = " << Point(i2.I2()) 3568 // << ", p3 = " << Point(pi3) 3569 // << ", p4 = " << Point(pi4) 3570 << ", rad = " << rad << endl; 3571 */ 3572 } 3573 else 3574 edges.Set (i2, i); 3575 } 3576 } 3577 3578 3579 // Restrict h due to line segments 3580 3581 for (i = 1; i <= GetNSeg(); i++) 3582 { 3583 const Segment & seg = LineSegment(i); 3584 const Point3d & p1 = Point(seg[0]); 3585 const Point3d & p2 = Point(seg[1]); 3586 RestrictLocalH (Center (p1, p2), Dist (p1, p2)); 3587 } 3588 3589 3590 3591 /* 3592 3593 3594 int i, j; 3595 int np = GetNP(); 3596 int nseg = GetNSeg(); 3597 int nse = GetNSE(); 3598 3599 NgArray<Vec3d> normals(np); 3600 NgBitArray linepoint(np); 3601 3602 linepoint.Clear(); 3603 for (i = 1; i <= nseg; i++) 3604 { 3605 linepoint.Set (LineSegment(i)[0]); 3606 linepoint.Set (LineSegment(i)[1]); 3607 } 3608 3609 for (i = 1; i <= np; i++) 3610 normals.Elem(i) = Vec3d(0,0,0); 3611 3612 for (i = 1; i <= nse; i++) 3613 { 3614 Element2d & el = SurfaceElement(i); 3615 Vec3d nf = Cross (Vec3d (Point (el.PNum(1)), Point(el.PNum(2))), 3616 Vec3d (Point (el.PNum(1)), Point(el.PNum(3)))); 3617 for (j = 1; j <= 3; j++) 3618 normals.Elem(el.PNum(j)) += nf; 3619 } 3620 3621 for (i = 1; i <= np; i++) 3622 normals.Elem(i) /= (1e-12 + normals.Elem(i).Length()); 3623 3624 for (i = 1; i <= nse; i++) 3625 { 3626 Element2d & el = SurfaceElement(i); 3627 Vec3d nf = Cross (Vec3d (Point (el.PNum(1)), Point(el.PNum(2))), 3628 Vec3d (Point (el.PNum(1)), Point(el.PNum(3)))); 3629 nf /= nf.Length(); 3630 Point3d c = Center (Point(el.PNum(1)), 3631 Point(el.PNum(2)), 3632 Point(el.PNum(3))); 3633 3634 for (j = 1; j <= 3; j++) 3635 { 3636 if (!linepoint.Test (el.PNum(j))) 3637 { 3638 double dist = Dist (c, Point(el.PNum(j))); 3639 double dn = (nf - normals.Get(el.PNum(j))).Length(); 3640 3641 RestrictLocalH (Point(el.PNum(j)), dist / (dn+1e-12) /elperr); 3642 } 3643 } 3644 } 3645 */ 3646 } 3647 3648 RestrictLocalH(resthtype rht,int nr,double loch)3649 void Mesh :: RestrictLocalH (resthtype rht, int nr, double loch) 3650 { 3651 int i; 3652 switch (rht) 3653 { 3654 case RESTRICTH_FACE: 3655 { 3656 for (i = 1; i <= GetNSE(); i++) 3657 { 3658 const Element2d & sel = SurfaceElement(i); 3659 if (sel.GetIndex() == nr) 3660 RestrictLocalH (RESTRICTH_SURFACEELEMENT, i, loch); 3661 } 3662 break; 3663 } 3664 case RESTRICTH_EDGE: 3665 { 3666 for (i = 1; i <= GetNSeg(); i++) 3667 { 3668 const Segment & seg = LineSegment(i); 3669 if (seg.edgenr == nr) 3670 RestrictLocalH (RESTRICTH_SEGMENT, i, loch); 3671 } 3672 break; 3673 } 3674 case RESTRICTH_POINT: 3675 { 3676 RestrictLocalH (Point (nr), loch); 3677 break; 3678 } 3679 3680 case RESTRICTH_SURFACEELEMENT: 3681 { 3682 const Element2d & sel = SurfaceElement(nr); 3683 Point3d p = Center (Point(sel.PNum(1)), 3684 Point(sel.PNum(2)), 3685 Point(sel.PNum(3))); 3686 RestrictLocalH (p, loch); 3687 break; 3688 } 3689 case RESTRICTH_SEGMENT: 3690 { 3691 const Segment & seg = LineSegment(nr); 3692 RestrictLocalHLine (Point (seg[0]), Point(seg[1]), loch); 3693 break; 3694 } 3695 } 3696 } 3697 3698 LoadLocalMeshSize(const string & meshsizefilename)3699 void Mesh :: LoadLocalMeshSize (const string & meshsizefilename) 3700 { 3701 // Philippose - 10/03/2009 3702 // Improve error checking when loading and reading 3703 // the local mesh size file 3704 3705 if (meshsizefilename.empty()) return; 3706 3707 ifstream msf(meshsizefilename.c_str()); 3708 3709 // Philippose - 09/03/2009 3710 // Adding print message information in case the specified 3711 // does not exist, or does not load successfully due to 3712 // other reasons such as access rights, etc... 3713 if (!msf) 3714 { 3715 PrintMessage(3, "Error loading mesh size file: ", meshsizefilename, "....","Skipping!"); 3716 return; 3717 } 3718 3719 PrintMessage (3, "Load local mesh-size file: ", meshsizefilename); 3720 3721 int nmsp = 0; 3722 int nmsl = 0; 3723 3724 msf >> nmsp; 3725 if(!msf.good()) 3726 throw NgException ("Mesh-size file error: No points found\n"); 3727 3728 if(nmsp > 0) 3729 PrintMessage (4, "Number of mesh-size restriction points: ", nmsp); 3730 3731 for (int i = 0; i < nmsp; i++) 3732 { 3733 Point3d pi; 3734 double hi; 3735 msf >> pi.X() >> pi.Y() >> pi.Z(); 3736 msf >> hi; 3737 if (!msf.good()) 3738 throw NgException ("Mesh-size file error: Number of points don't match specified list size\n"); 3739 RestrictLocalH (pi, hi); 3740 } 3741 3742 msf >> nmsl; 3743 if(!msf.good()) 3744 throw NgException ("Mesh-size file error: No line definitions found\n"); 3745 3746 if(nmsl > 0) 3747 PrintMessage (4, "Number of mesh-size restriction lines: ", nmsl); 3748 3749 for (int i = 0; i < nmsl; i++) 3750 { 3751 Point3d p1, p2; 3752 double hi; 3753 msf >> p1.X() >> p1.Y() >> p1.Z(); 3754 msf >> p2.X() >> p2.Y() >> p2.Z(); 3755 msf >> hi; 3756 if (!msf.good()) 3757 throw NgException ("Mesh-size file error: Number of line definitions don't match specified list size\n"); 3758 RestrictLocalHLine (p1, p2, hi); 3759 } 3760 3761 msf.close(); 3762 } 3763 3764 3765 GetBox(Point3d & pmin,Point3d & pmax,int dom) const3766 void Mesh :: GetBox (Point3d & pmin, Point3d & pmax, int dom) const 3767 { 3768 if (points.Size() == 0) 3769 { 3770 pmin = pmax = Point3d(0,0,0); 3771 return; 3772 } 3773 3774 if (dom <= 0) 3775 { 3776 pmin = Point3d (1e10, 1e10, 1e10); 3777 pmax = Point3d (-1e10, -1e10, -1e10); 3778 3779 // for (PointIndex pi = points.Begin(); pi < points.End(); pi++) 3780 for (PointIndex pi : points.Range()) 3781 { 3782 pmin.SetToMin ( (*this) [pi] ); 3783 pmax.SetToMax ( (*this) [pi] ); 3784 } 3785 } 3786 else 3787 { 3788 int j, nse = GetNSE(); 3789 SurfaceElementIndex sei; 3790 3791 pmin = Point3d (1e10, 1e10, 1e10); 3792 pmax = Point3d (-1e10, -1e10, -1e10); 3793 for (sei = 0; sei < nse; sei++) 3794 { 3795 const Element2d & el = (*this)[sei]; 3796 if (el.IsDeleted() ) continue; 3797 3798 if (dom == -1 || el.GetIndex() == dom) 3799 { 3800 for (j = 0; j < 3; j++) 3801 { 3802 pmin.SetToMin ( (*this) [el[j]] ); 3803 pmax.SetToMax ( (*this) [el[j]] ); 3804 } 3805 } 3806 } 3807 } 3808 3809 if (pmin.X() > 0.5e10) 3810 { 3811 pmin = pmax = Point3d(0,0,0); 3812 } 3813 } 3814 3815 3816 3817 GetBox(Point3d & pmin,Point3d & pmax,POINTTYPE ptyp) const3818 void Mesh :: GetBox (Point3d & pmin, Point3d & pmax, POINTTYPE ptyp) const 3819 { 3820 if (points.Size() == 0) 3821 { 3822 pmin = pmax = Point3d(0,0,0); 3823 return; 3824 } 3825 3826 pmin = Point3d (1e10, 1e10, 1e10); 3827 pmax = Point3d (-1e10, -1e10, -1e10); 3828 3829 // for (PointIndex pi = points.Begin(); pi < points.End(); pi++) 3830 for (PointIndex pi : points.Range()) 3831 if (points[pi].Type() <= ptyp) 3832 { 3833 pmin.SetToMin ( (*this) [pi] ); 3834 pmax.SetToMax ( (*this) [pi] ); 3835 } 3836 } 3837 3838 3839 3840 ElementError(int eli,const MeshingParameters & mp) const3841 double Mesh :: ElementError (int eli, const MeshingParameters & mp) const 3842 { 3843 const Element & el = volelements[eli-1]; 3844 return CalcTetBadness (points[el[0]], points[el[1]], 3845 points[el[2]], points[el[3]], -1, mp); 3846 } 3847 AddLockedPoint(PointIndex pi)3848 void Mesh :: AddLockedPoint (PointIndex pi) 3849 { 3850 lockedpoints.Append (pi); 3851 } 3852 ClearLockedPoints()3853 void Mesh :: ClearLockedPoints () 3854 { 3855 lockedpoints.SetSize (0); 3856 } 3857 3858 3859 Compress()3860 void Mesh :: Compress () 3861 { 3862 static Timer t("Mesh::Compress"); RegionTimer reg(t); 3863 NgLock lock(mutex); 3864 lock.Lock(); 3865 3866 Array<PointIndex,PointIndex> op2np(GetNP()); 3867 Array<bool, PointIndex> pused(GetNP()); 3868 3869 /* 3870 (*testout) << "volels: " << endl; 3871 for (i = 1; i <= volelements.Size(); i++) 3872 { 3873 for (j = 1; j <= volelements.Get(i).GetNP(); j++) 3874 (*testout) << volelements.Get(i).PNum(j) << " "; 3875 (*testout) << endl; 3876 } 3877 (*testout) << "np: " << GetNP() << endl; 3878 */ 3879 3880 for (int i = 0; i < volelements.Size(); i++) 3881 if (volelements[i][0] <= PointIndex::BASE-1 || 3882 volelements[i].IsDeleted()) 3883 { 3884 volelements.DeleteElement(i); 3885 i--; 3886 } 3887 3888 3889 for (int i = 0; i < surfelements.Size(); i++) 3890 if (surfelements[i].IsDeleted()) 3891 { 3892 surfelements.DeleteElement(i); 3893 i--; 3894 } 3895 3896 for (int i = 0; i < segments.Size(); i++) 3897 if (segments[i][0] <= PointIndex::BASE-1) 3898 { 3899 segments.DeleteElement(i); 3900 i--; 3901 } 3902 3903 for(int i=0; i < segments.Size(); i++) 3904 if(segments[i].edgenr < 0) 3905 segments.DeleteElement(i--); 3906 3907 pused = false; 3908 /* 3909 for (int i = 0; i < volelements.Size(); i++) 3910 { 3911 const Element & el = volelements[i]; 3912 for (int j = 0; j < el.GetNP(); j++) 3913 pused[el[j]] = true; 3914 } 3915 */ 3916 /* 3917 for (const Element & el : volelements) 3918 for (PointIndex pi : el.PNums()) 3919 pused[pi] = true; 3920 */ 3921 3922 ParallelForRange 3923 (volelements.Range(), [&] (auto myrange) 3924 { 3925 for (const Element & el : volelements.Range(myrange)) 3926 for (PointIndex pi : el.PNums()) 3927 pused[pi] = true; 3928 }); 3929 3930 /* 3931 for (int i = 0; i < surfelements.Size(); i++) 3932 { 3933 const Element2d & el = surfelements[i]; 3934 for (int j = 0; j < el.GetNP(); j++) 3935 pused[el[j]] = true; 3936 } 3937 */ 3938 ParallelForRange 3939 (surfelements.Range(), [&] (auto myrange) 3940 { 3941 for (const Element2d & el : surfelements.Range(myrange)) 3942 for (PointIndex pi : el.PNums()) 3943 pused[pi] = true; 3944 }); 3945 3946 for (int i = 0; i < segments.Size(); i++) 3947 { 3948 const Segment & seg = segments[i]; 3949 for (int j = 0; j < seg.GetNP(); j++) 3950 pused[seg[j]] = true; 3951 } 3952 3953 for (int i = 0; i < openelements.Size(); i++) 3954 { 3955 const Element2d & el = openelements[i]; 3956 for (int j = 0; j < el.GetNP(); j++) 3957 pused[el[j]] = true; 3958 } 3959 3960 for (int i = 0; i < lockedpoints.Size(); i++) 3961 pused[lockedpoints[i]] = true; 3962 3963 3964 /* 3965 // compress points doesn't work for identified points ! 3966 if (identifiedpoints) 3967 { 3968 for (i = 1; i <= identifiedpoints->GetNBags(); i++) 3969 if (identifiedpoints->GetBagSize(i)) 3970 { 3971 pused.Set (); 3972 break; 3973 } 3974 } 3975 */ 3976 // pused.Set(); 3977 3978 3979 { 3980 Array<MeshPoint> hpoints; 3981 int npi = PointIndex::BASE; 3982 for (PointIndex pi : points.Range()) 3983 if (pused[pi]) 3984 { 3985 op2np[pi] = npi; 3986 npi++; 3987 hpoints.Append (points[pi]); 3988 } 3989 else 3990 { 3991 op2np[pi].Invalidate(); 3992 } 3993 3994 points.SetSize(0); 3995 for (int i = 0; i < hpoints.Size(); i++) 3996 points.Append (hpoints[i]); 3997 } 3998 3999 /* 4000 for (int i = 1; i <= volelements.Size(); i++) 4001 { 4002 Element & el = VolumeElement(i); 4003 for (int j = 0; j < el.GetNP(); j++) 4004 el[j] = op2np[el[j]]; 4005 } 4006 */ 4007 ParallelForRange 4008 (volelements.Range(), [&] (auto myrange) 4009 { 4010 for (Element & el : volelements.Range(myrange)) 4011 for (PointIndex & pi : el.PNums()) 4012 pi = op2np[pi]; 4013 }); 4014 4015 /* 4016 for (int i = 1; i <= surfelements.Size(); i++) 4017 { 4018 Element2d & el = SurfaceElement(i); 4019 for (int j = 0; j < el.GetNP(); j++) 4020 el[j] = op2np[el[j]]; 4021 } 4022 */ 4023 ParallelForRange 4024 (surfelements.Range(), [&] (auto myrange) 4025 { 4026 for (Element2d & el : surfelements.Range(myrange)) 4027 for (PointIndex & pi : el.PNums()) 4028 pi = op2np[pi]; 4029 }); 4030 4031 4032 for (int i = 0; i < segments.Size(); i++) 4033 { 4034 Segment & seg = segments[i]; 4035 for (int j = 0; j < seg.GetNP(); j++) 4036 seg[j] = op2np[seg[j]]; 4037 } 4038 4039 for (int i = 1; i <= openelements.Size(); i++) 4040 { 4041 Element2d & el = openelements.Elem(i); 4042 for (int j = 0; j < el.GetNP(); j++) 4043 el[j] = op2np[el[j]]; 4044 } 4045 4046 4047 for (int i = 0; i < lockedpoints.Size(); i++) 4048 lockedpoints[i] = op2np[lockedpoints[i]]; 4049 /* 4050 for (int i = 0; i < facedecoding.Size(); i++) 4051 facedecoding[i].firstelement = -1; 4052 for (int i = surfelements.Size()-1; i >= 0; i--) 4053 { 4054 int ind = surfelements[i].GetIndex(); 4055 surfelements[i].next = facedecoding[ind-1].firstelement; 4056 facedecoding[ind-1].firstelement = i; 4057 } 4058 */ 4059 RebuildSurfaceElementLists (); 4060 CalcSurfacesOfNode(); 4061 4062 4063 // FindOpenElements(); 4064 timestamp = NextTimeStamp(); 4065 lock.UnLock(); 4066 } 4067 OrderElements()4068 void Mesh :: OrderElements() 4069 { 4070 for (auto & el : surfelements) 4071 { 4072 if (el.GetType() == TRIG) 4073 while (el[0] > el[1] || el[0] > el[2]) 4074 { // rotate element 4075 auto hp = el[0]; 4076 el[0] = el[1]; 4077 el[1] = el[2]; 4078 el[2] = hp; 4079 auto hgi = el.GeomInfoPi(1); 4080 el.GeomInfoPi(1) = el.GeomInfoPi(2); 4081 el.GeomInfoPi(2) = el.GeomInfoPi(3); 4082 el.GeomInfoPi(3) = hgi; 4083 } 4084 } 4085 4086 for (auto & el : volelements) 4087 if (el.GetType() == TET) 4088 { 4089 // lowest index first ... 4090 int mini = 0; 4091 for (int i = 1; i < 4; i++) 4092 if (el[i] < el[mini]) mini = i; 4093 if (mini != 0) 4094 { // swap 0 with mini, and the other two ... 4095 int i3 = -1, i4 = -1; 4096 for (int i = 1; i < 4; i++) 4097 if (i != mini) 4098 { 4099 i4 = i3; 4100 i3 = i; 4101 } 4102 swap (el[0], el[mini]); 4103 swap (el[i3], el[i4]); 4104 } 4105 4106 while (el[1] > el[2] || el[1] > el[3]) 4107 { // rotate element to move second index to second position 4108 auto hp = el[1]; 4109 el[1] = el[2]; 4110 el[2] = el[3]; 4111 el[3] = hp; 4112 } 4113 } 4114 } 4115 CheckConsistentBoundary() const4116 int Mesh :: CheckConsistentBoundary () const 4117 { 4118 int nf = GetNOpenElements(); 4119 INDEX_2_HASHTABLE<int> edges(nf+2); 4120 INDEX_2 i2, i2s, edge; 4121 int err = 0; 4122 4123 for (int i = 1; i <= nf; i++) 4124 { 4125 const Element2d & sel = OpenElement(i); 4126 4127 for (int j = 1; j <= sel.GetNP(); j++) 4128 { 4129 i2.I1() = sel.PNumMod(j); 4130 i2.I2() = sel.PNumMod(j+1); 4131 4132 int sign = (i2.I2() > i2.I1()) ? 1 : -1; 4133 i2.Sort(); 4134 if (!edges.Used (i2)) 4135 edges.Set (i2, 0); 4136 edges.Set (i2, edges.Get(i2) + sign); 4137 } 4138 } 4139 4140 for (int i = 1; i <= edges.GetNBags(); i++) 4141 for (int j = 1; j <= edges.GetBagSize(i); j++) 4142 { 4143 int cnt = 0; 4144 edges.GetData (i, j, i2, cnt); 4145 if (cnt) 4146 { 4147 PrintError ("Edge ", i2.I1() , " - ", i2.I2(), " multiple times in surface mesh"); 4148 4149 (*testout) << "Edge " << i2 << " multiple times in surface mesh" << endl; 4150 i2s = i2; 4151 i2s.Sort(); 4152 for (int k = 1; k <= nf; k++) 4153 { 4154 const Element2d & sel = OpenElement(k); 4155 for (int l = 1; l <= sel.GetNP(); l++) 4156 { 4157 edge.I1() = sel.PNumMod(l); 4158 edge.I2() = sel.PNumMod(l+1); 4159 edge.Sort(); 4160 4161 if (edge == i2s) 4162 (*testout) << "edge of element " << sel << endl; 4163 } 4164 } 4165 4166 4167 err = 2; 4168 } 4169 } 4170 4171 return err; 4172 } 4173 4174 4175 CheckOverlappingBoundary()4176 int Mesh :: CheckOverlappingBoundary () 4177 { 4178 static Timer t("Mesh::CheckOverlappingBoundary"); RegionTimer reg(t); 4179 4180 Point3d pmin, pmax; 4181 GetBox (pmin, pmax); 4182 BoxTree<3, SurfaceElementIndex> setree(pmin, pmax); 4183 // NgArray<SurfaceElementIndex> inters; 4184 4185 bool overlap = 0; 4186 bool incons_layers = 0; 4187 4188 for (Element2d & el : SurfaceElements()) 4189 el.badel = false; 4190 4191 for (SurfaceElementIndex sei : Range(SurfaceElements())) 4192 { 4193 const Element2d & tri = SurfaceElement(sei); 4194 4195 Box<3> box(Box<3>::EMPTY_BOX); 4196 for (PointIndex pi : tri.PNums()) 4197 box.Add (Point(pi)); 4198 4199 box.Increase(1e-3*box.Diam()); 4200 setree.Insert (box, sei); 4201 } 4202 4203 std::mutex m; 4204 // for (SurfaceElementIndex sei : Range(SurfaceElements())) 4205 ParallelForRange 4206 (Range(SurfaceElements()), [&] (auto myrange) 4207 { 4208 for (SurfaceElementIndex sei : myrange) 4209 { 4210 const Element2d & tri = SurfaceElement(sei); 4211 4212 Box<3> box(Box<3>::EMPTY_BOX); 4213 for (PointIndex pi : tri.PNums()) 4214 box.Add (Point(pi)); 4215 4216 setree.GetFirstIntersecting 4217 (box.PMin(), box.PMax(), 4218 [&] (SurfaceElementIndex sej) 4219 { 4220 const Element2d & tri2 = SurfaceElement(sej); 4221 4222 if ( (*this)[tri[0]].GetLayer() != (*this)[tri2[0]].GetLayer()) 4223 return false; 4224 4225 if ( (*this)[tri[0]].GetLayer() != (*this)[tri[1]].GetLayer() || 4226 (*this)[tri[0]].GetLayer() != (*this)[tri[2]].GetLayer()) 4227 { 4228 incons_layers = 1; 4229 // cout << "inconsistent layers in triangle" << endl; 4230 } 4231 4232 const netgen::Point<3> *trip1[3], *trip2[3]; 4233 for (int k = 0; k < 3; k++) 4234 { 4235 trip1[k] = &Point (tri[k]); 4236 trip2[k] = &Point (tri2[k]); 4237 } 4238 4239 if (IntersectTriangleTriangle (&trip1[0], &trip2[0])) 4240 { 4241 overlap = 1; 4242 lock_guard<std::mutex> guard(m); 4243 if(!incons_layers) 4244 { 4245 PrintWarning ("Intersecting elements " 4246 ,int(sei), " and ", int(sej)); 4247 4248 (*testout) << "Intersecting: " << endl; 4249 (*testout) << "openelement " << sei << " with open element " << sej << endl; 4250 4251 cout << "el1 = " << tri << endl; 4252 cout << "el2 = " << tri2 << endl; 4253 cout << "layer1 = " << (*this)[tri[0]].GetLayer() << endl; 4254 cout << "layer2 = " << (*this)[tri2[0]].GetLayer() << endl; 4255 } 4256 4257 for (int k = 1; k <= 3; k++) 4258 (*testout) << tri.PNum(k) << " "; 4259 (*testout) << endl; 4260 for (int k = 1; k <= 3; k++) 4261 (*testout) << tri2.PNum(k) << " "; 4262 (*testout) << endl; 4263 4264 for (int k = 0; k <= 2; k++) 4265 (*testout) << *trip1[k] << " "; 4266 (*testout) << endl; 4267 for (int k = 0; k <= 2; k++) 4268 (*testout) << *trip2[k] << " "; 4269 (*testout) << endl; 4270 4271 (*testout) << "Face1 = " << GetFaceDescriptor(tri.GetIndex()) << endl; 4272 (*testout) << "Face1 = " << GetFaceDescriptor(tri2.GetIndex()) << endl; 4273 4274 SurfaceElement(sei).badel = 1; 4275 SurfaceElement(sej).badel = 1; 4276 } 4277 return false; 4278 }); 4279 } 4280 }); 4281 // bug 'fix' 4282 if (incons_layers) overlap = 0; 4283 4284 return overlap; 4285 } 4286 4287 CheckVolumeMesh() const4288 int Mesh :: CheckVolumeMesh () const 4289 { 4290 PrintMessage (3, "Checking volume mesh"); 4291 4292 int ne = GetNE(); 4293 DenseMatrix dtrans(3,3); 4294 int i, j; 4295 4296 PrintMessage (5, "elements: ", ne); 4297 for (i = 1; i <= ne; i++) 4298 { 4299 Element & el = (Element&) VolumeElement(i); 4300 el.flags.badel = 0; 4301 int nip = el.GetNIP(); 4302 for (j = 1; j <= nip; j++) 4303 { 4304 el.GetTransformation (j, Points(), dtrans); 4305 double det = dtrans.Det(); 4306 if (det > 0) 4307 { 4308 PrintError ("Element ", i , " has wrong orientation"); 4309 el.flags.badel = 1; 4310 } 4311 } 4312 } 4313 4314 return 0; 4315 } 4316 4317 // Search for surface trigs with same vertices ( may happen for instance with close surfaces in stl geometies ) FindIllegalTrigs()4318 int Mesh :: FindIllegalTrigs () 4319 { 4320 // Temporary table to store the vertex numbers of all triangles 4321 INDEX_3_CLOSED_HASHTABLE<int> temp_tab(3*GetNSE() + 1); 4322 size_t cnt = 0; 4323 for (SurfaceElementIndex sei = 0; sei < GetNSE(); sei++) 4324 { 4325 const Element2d & sel = surfelements[sei]; 4326 if (sel.IsDeleted()) continue; 4327 4328 INDEX_3 i3(sel[0], sel[1], sel[2]); 4329 i3.Sort(); 4330 if(temp_tab.Used(i3)) 4331 { 4332 temp_tab.Set (i3, -1); 4333 cnt++; 4334 } 4335 else 4336 { 4337 temp_tab.Set (i3, sei); 4338 } 4339 } 4340 4341 illegal_trigs = make_unique<INDEX_3_CLOSED_HASHTABLE<int>> (2*cnt+1); 4342 for (SurfaceElementIndex sei = 0; sei < GetNSE(); sei++) 4343 { 4344 const Element2d & sel = surfelements[sei]; 4345 if (sel.IsDeleted()) continue; 4346 4347 INDEX_3 i3(sel[0], sel[1], sel[2]); 4348 i3.Sort(); 4349 if(temp_tab.Get(i3)==-1) 4350 illegal_trigs -> Set (i3, 1); 4351 } 4352 return cnt; 4353 } 4354 LegalTrig(const Element2d & el) const4355 bool Mesh :: LegalTrig (const Element2d & el) const 4356 { 4357 if(illegal_trigs) 4358 { 4359 INDEX_3 i3 (el[0], el[1], el[2]); 4360 i3.Sort(); 4361 if(illegal_trigs->Used(i3)) 4362 return false; 4363 } 4364 4365 return 1; 4366 if ( /* hp */ 1) // needed for old, simple hp-refinement 4367 { 4368 // trigs with 2 or more segments are illegal 4369 int i; 4370 int nseg = 0; 4371 4372 if (!segmentht) 4373 { 4374 cerr << "no segmentht allocated" << endl; 4375 return 0; 4376 } 4377 4378 // Point3d cp(0.5, 0.5, 0.5); 4379 for (i = 1; i <= 3; i++) 4380 { 4381 INDEX_2 i2(el.PNumMod (i), el.PNumMod (i+1)); 4382 i2.Sort(); 4383 if (segmentht -> Used (i2)) 4384 nseg++; 4385 } 4386 if (nseg >= 2) 4387 return 0; 4388 } 4389 return 1; 4390 } 4391 CalcTotalBad(const MeshingParameters & mp)4392 double Mesh :: CalcTotalBad (const MeshingParameters & mp ) 4393 { 4394 static Timer t("CalcTotalBad"); RegionTimer reg(t); 4395 static constexpr int n_classes = 20; 4396 4397 double sum = 0; 4398 4399 tets_in_qualclass.SetSize(n_classes); 4400 tets_in_qualclass = 0; 4401 4402 ParallelForRange( IntRange(volelements.Size()), [&] (auto myrange) 4403 { 4404 double local_sum = 0.0; 4405 double teterrpow = mp.opterrpow; 4406 4407 std::array<int,n_classes> classes_local{}; 4408 4409 for (auto i : myrange) 4410 { 4411 double elbad = pow (max2(CalcBad (points, volelements[i], 0, mp),1e-10), 1/teterrpow); 4412 4413 int qualclass = int (n_classes / elbad + 1); 4414 if (qualclass < 1) qualclass = 1; 4415 if (qualclass > n_classes) qualclass = n_classes; 4416 classes_local[qualclass-1]++; 4417 4418 local_sum += elbad; 4419 } 4420 4421 AtomicAdd(sum, local_sum); 4422 4423 for (auto i : Range(n_classes)) 4424 AsAtomic(tets_in_qualclass[i]) += classes_local[i]; 4425 }); 4426 4427 return sum; 4428 } 4429 4430 4431 4432 4433 /// LegalTet2(Element & el) const4434 bool Mesh :: LegalTet2 (Element & el) const 4435 { 4436 // static int timer1 = NgProfiler::CreateTimer ("Legaltet2"); 4437 4438 // Test, whether 4 points have a common surface plus 4439 // at least 4 edges at the boundary 4440 4441 if(!boundaryedges) 4442 const_cast<Mesh *>(this)->BuildBoundaryEdges(); 4443 4444 4445 // non-tets are always legal 4446 if (el.GetType() != TET) 4447 { 4448 el.SetLegal (1); 4449 return 1; 4450 } 4451 4452 POINTTYPE pointtype[4]; 4453 for(int i = 0; i < 4; i++) 4454 pointtype[i] = (*this)[el[i]].Type(); 4455 4456 4457 4458 // element has at least 2 inner points ---> legal 4459 int cnti = 0; 4460 for (int j = 0; j < 4; j++) 4461 if ( pointtype[j] == INNERPOINT) 4462 { 4463 cnti++; 4464 if (cnti >= 2) 4465 { 4466 el.SetLegal (1); 4467 return 1; 4468 } 4469 } 4470 4471 4472 4473 // which faces are boundary faces ? 4474 int bface[4]; 4475 for (int i = 0; i < 4; i++) 4476 { 4477 bface[i] = surfelementht->Used (INDEX_3::Sort(el[gftetfacesa[i][0]], 4478 el[gftetfacesa[i][1]], 4479 el[gftetfacesa[i][2]])); 4480 } 4481 4482 int bedge[4][4]; 4483 int segedge[4][4]; 4484 static const int pi3map[4][4] = { { -1, 2, 1, 1 }, 4485 { 2, -1, 0, 0 }, 4486 { 1, 0, -1, 0 }, 4487 { 1, 0, 0, -1 } }; 4488 4489 static const int pi4map[4][4] = { { -1, 3, 3, 2 }, 4490 { 3, -1, 3, 2 }, 4491 { 3, 3, -1, 1 }, 4492 { 2, 2, 1, -1 } }; 4493 4494 4495 for (int i = 0; i < 4; i++) 4496 for (int j = 0; j < i; j++) 4497 { 4498 bool sege = false, be = false; 4499 4500 int pos = boundaryedges -> Position0(INDEX_2::Sort(el[i], el[j])); 4501 if (pos != -1) 4502 { 4503 be = true; 4504 if (boundaryedges -> GetData0(pos) == 2) 4505 sege = true; 4506 } 4507 4508 segedge[j][i] = segedge[i][j] = sege; 4509 bedge[j][i] = bedge[i][j] = be; 4510 } 4511 4512 // two boundary faces and no edge is illegal 4513 for (int i = 0; i < 3; i++) 4514 for (int j = i+1; j < 4; j++) 4515 { 4516 if (bface[i] && bface[j]) 4517 if (!segedge[pi3map[i][j]][pi4map[i][j]]) 4518 { 4519 // 2 boundary faces withoud edge in between 4520 el.SetLegal (0); 4521 return 0; 4522 } 4523 } 4524 4525 // three boundary edges meeting in a Surface point 4526 for (int i = 0; i < 4; i++) 4527 { 4528 if ( pointtype[i] == SURFACEPOINT) 4529 { 4530 bool alledges = 1; 4531 for (int j = 0; j < 4; j++) 4532 if (j != i && !bedge[i][j]) 4533 { 4534 alledges = 0; 4535 break; 4536 } 4537 if (alledges) 4538 { 4539 // cout << "tet illegal due to unmarked node" << endl; 4540 el.SetLegal (0); 4541 return 0; 4542 } 4543 } 4544 } 4545 4546 4547 4548 for (int fnr = 0; fnr < 4; fnr++) 4549 if (!bface[fnr]) 4550 for (int i = 0; i < 4; i++) 4551 if (i != fnr) 4552 { 4553 int pi1 = pi3map[i][fnr]; 4554 int pi2 = pi4map[i][fnr]; 4555 4556 if ( pointtype[i] == SURFACEPOINT) 4557 { 4558 // two connected edges on surface, but no face 4559 if (bedge[i][pi1] && bedge[i][pi2]) 4560 { 4561 el.SetLegal (0); 4562 return 0; 4563 } 4564 } 4565 4566 if ( pointtype[i] == EDGEPOINT) 4567 { 4568 // connected surface edge and edge edge, but no face 4569 if ( (bedge[i][pi1] && segedge[i][pi2]) || 4570 (bedge[i][pi2] && segedge[i][pi1]) ) 4571 { 4572 el.SetLegal (0); 4573 return 0; 4574 } 4575 } 4576 4577 } 4578 4579 4580 el.SetLegal (1); 4581 return 1; 4582 4583 } 4584 4585 4586 GetNDomains() const4587 int Mesh :: GetNDomains() const 4588 { 4589 int ndom = 0; 4590 4591 for (int k = 0; k < facedecoding.Size(); k++) 4592 { 4593 if (facedecoding[k].DomainIn() > ndom) 4594 ndom = facedecoding[k].DomainIn(); 4595 if (facedecoding[k].DomainOut() > ndom) 4596 ndom = facedecoding[k].DomainOut(); 4597 } 4598 4599 return ndom; 4600 } 4601 SetDimension(int dim)4602 void Mesh :: SetDimension (int dim) 4603 { 4604 if (dimension == 3 && dim == 2) 4605 { 4606 // change mesh-dim from 3 to 2 (currently needed for OCC) 4607 for (auto str : materials) 4608 delete str; 4609 materials.SetSize(0); 4610 for (auto str : bcnames) 4611 materials.Append(str); 4612 bcnames.SetSize(0); 4613 for (auto str : cd2names) 4614 bcnames.Append(str); 4615 cd2names.SetSize(0); 4616 for (auto str : cd3names) 4617 cd2names.Append(str); 4618 cd3names.SetSize(0); 4619 4620 for (auto & seg : LineSegments()) 4621 seg.si = seg.edgenr; 4622 } 4623 dimension = dim; 4624 } 4625 SurfaceMeshOrientation()4626 void Mesh :: SurfaceMeshOrientation () 4627 { 4628 int i, j; 4629 int nse = GetNSE(); 4630 4631 NgBitArray used(nse); 4632 used.Clear(); 4633 INDEX_2_HASHTABLE<int> edges(nse+1); 4634 4635 bool haschanged = 0; 4636 4637 4638 const Element2d & tri = SurfaceElement(1); 4639 for (j = 1; j <= 3; j++) 4640 { 4641 INDEX_2 i2(tri.PNumMod(j), tri.PNumMod(j+1)); 4642 edges.Set (i2, 1); 4643 } 4644 used.Set(1); 4645 4646 bool unused; 4647 do 4648 { 4649 bool changed; 4650 do 4651 { 4652 changed = 0; 4653 for (i = 1; i <= nse; i++) 4654 if (!used.Test(i)) 4655 { 4656 Element2d & el = surfelements[i-1]; 4657 int found = 0, foundrev = 0; 4658 for (j = 1; j <= 3; j++) 4659 { 4660 INDEX_2 i2(el.PNumMod(j), el.PNumMod(j+1)); 4661 if (edges.Used(i2)) 4662 foundrev = 1; 4663 swap (i2.I1(), i2.I2()); 4664 if (edges.Used(i2)) 4665 found = 1; 4666 } 4667 4668 if (found || foundrev) 4669 { 4670 if (foundrev) 4671 swap (el.PNum(2), el.PNum(3)); 4672 4673 changed = 1; 4674 for (j = 1; j <= 3; j++) 4675 { 4676 INDEX_2 i2(el.PNumMod(j), el.PNumMod(j+1)); 4677 edges.Set (i2, 1); 4678 } 4679 used.Set (i); 4680 } 4681 } 4682 if (changed) 4683 haschanged = 1; 4684 } 4685 while (changed); 4686 4687 4688 unused = 0; 4689 for (i = 1; i <= nse; i++) 4690 if (!used.Test(i)) 4691 { 4692 unused = 1; 4693 const Element2d & tri = SurfaceElement(i); 4694 for (j = 1; j <= 3; j++) 4695 { 4696 INDEX_2 i2(tri.PNumMod(j), tri.PNumMod(j+1)); 4697 edges.Set (i2, 1); 4698 } 4699 used.Set(i); 4700 break; 4701 } 4702 } 4703 while (unused); 4704 4705 if (haschanged) 4706 timestamp = NextTimeStamp(); 4707 } 4708 4709 Split2Tets()4710 void Mesh :: Split2Tets() 4711 { 4712 PrintMessage (1, "Split To Tets"); 4713 bool has_prisms = 0; 4714 4715 int oldne = GetNE(); 4716 for (int i = 1; i <= oldne; i++) 4717 { 4718 Element el = VolumeElement(i); 4719 4720 if (el.GetType() == PRISM) 4721 { 4722 // prism, to 3 tets 4723 4724 // make minimal node to node 1 4725 int minpi=0; 4726 PointIndex minpnum; 4727 minpnum = GetNP() + 1; 4728 4729 for (int j = 1; j <= 6; j++) 4730 { 4731 if (el.PNum(j) < minpnum) 4732 { 4733 minpnum = el.PNum(j); 4734 minpi = j; 4735 } 4736 } 4737 4738 if (minpi >= 4) 4739 { 4740 for (int j = 1; j <= 3; j++) 4741 swap (el.PNum(j), el.PNum(j+3)); 4742 minpi -= 3; 4743 } 4744 4745 while (minpi > 1) 4746 { 4747 int hi = 0; 4748 for (int j = 0; j <= 3; j+= 3) 4749 { 4750 hi = el.PNum(1+j); 4751 el.PNum(1+j) = el.PNum(2+j); 4752 el.PNum(2+j) = el.PNum(3+j); 4753 el.PNum(3+j) = hi; 4754 } 4755 minpi--; 4756 } 4757 4758 /* 4759 version 1: edge from pi2 to pi6, 4760 version 2: edge from pi3 to pi5, 4761 */ 4762 4763 static const int ntets[2][12] = 4764 { { 1, 4, 5, 6, 1, 2, 3, 6, 1, 2, 5, 6 }, 4765 { 1, 4, 5, 6, 1, 2, 3, 5, 3, 1, 5, 6 } }; 4766 4767 const int * min2pi; 4768 4769 if (min2 (el.PNum(2), el.PNum(6)) < 4770 min2 (el.PNum(3), el.PNum(5))) 4771 { 4772 min2pi = &ntets[0][0]; 4773 // (*testout) << "version 1 "; 4774 } 4775 else 4776 { 4777 min2pi = &ntets[1][0]; 4778 // (*testout) << "version 2 "; 4779 } 4780 4781 4782 int firsttet = 1; 4783 for (int j = 1; j <= 3; j++) 4784 { 4785 Element nel(TET); 4786 for (int k = 1; k <= 4; k++) 4787 nel.PNum(k) = el.PNum(min2pi[4 * j + k - 5]); 4788 nel.SetIndex (el.GetIndex()); 4789 4790 int legal = 1; 4791 for (int k = 1; k <= 3; k++) 4792 for (int l = k+1; l <= 4; l++) 4793 if (nel.PNum(k) == nel.PNum(l)) 4794 legal = 0; 4795 4796 // (*testout) << nel << " "; 4797 if (legal) 4798 { 4799 if (firsttet) 4800 { 4801 VolumeElement(i) = nel; 4802 firsttet = 0; 4803 } 4804 else 4805 { 4806 AddVolumeElement(nel); 4807 } 4808 } 4809 } 4810 if (firsttet) cout << "no legal"; 4811 (*testout) << endl; 4812 } 4813 4814 4815 4816 else if (el.GetType() == HEX) 4817 { 4818 // hex to A) 2 prisms or B) to 5 tets 4819 4820 // make minimal node to node 1 4821 int minpi=0; 4822 PointIndex minpnum; 4823 minpnum = GetNP() + 1; 4824 4825 for (int j = 1; j <= 8; j++) 4826 { 4827 if (el.PNum(j) < minpnum) 4828 { 4829 minpnum = el.PNum(j); 4830 minpi = j; 4831 } 4832 } 4833 4834 if (minpi >= 5) 4835 { 4836 for (int j = 1; j <= 4; j++) 4837 swap (el.PNum(j), el.PNum(j+4)); 4838 minpi -= 4; 4839 } 4840 4841 while (minpi > 1) 4842 { 4843 int hi = 0; 4844 for (int j = 0; j <= 4; j+= 4) 4845 { 4846 hi = el.PNum(1+j); 4847 el.PNum(1+j) = el.PNum(2+j); 4848 el.PNum(2+j) = el.PNum(3+j); 4849 el.PNum(3+j) = el.PNum(4+j); 4850 el.PNum(4+j) = hi; 4851 } 4852 minpi--; 4853 } 4854 4855 4856 4857 static const int to_prisms[3][12] = 4858 { { 0, 1, 2, 4, 5, 6, 0, 2, 3, 4, 6, 7 }, 4859 { 0, 1, 5, 3, 2, 6, 0, 5, 4, 3, 6, 7 }, 4860 { 0, 7, 4, 1, 6, 5, 0, 3, 7, 1, 2, 6 }, 4861 }; 4862 4863 const int * min2pi = 0; 4864 if (min2 (el[4], el[6]) < min2 (el[5], el[7])) 4865 min2pi = &to_prisms[0][0]; 4866 else if (min2 (el[3], el[6]) < min2 (el[2], el[7])) 4867 min2pi = &to_prisms[1][0]; 4868 else if (min2 (el[1], el[6]) < min2 (el[2], el[5])) 4869 min2pi = &to_prisms[2][0]; 4870 4871 if (min2pi) 4872 { 4873 has_prisms = 1; 4874 for (int j = 0; j < 2; j++) 4875 { 4876 Element nel(PRISM); 4877 for (int k = 0; k < 6; k++) 4878 nel[k] = el[min2pi[6*j + k]]; 4879 nel.SetIndex (el.GetIndex()); 4880 4881 if (j == 0) 4882 VolumeElement(i) = nel; 4883 else 4884 AddVolumeElement(nel); 4885 } 4886 } 4887 else 4888 { 4889 // split to 5 tets 4890 4891 static const int to_tets[20] = 4892 { 4893 1, 2, 0, 5, 4894 3, 0, 2, 7, 4895 4, 5, 7, 0, 4896 6, 7, 5, 2, 4897 0, 2, 7, 5 4898 }; 4899 4900 for (int j = 0; j < 5; j++) 4901 { 4902 Element nel(TET); 4903 for (int k = 0; k < 4; k++) 4904 nel[k] = el[to_tets[4*j + k]]; 4905 nel.SetIndex (el.GetIndex()); 4906 4907 if (j == 0) 4908 VolumeElement(i) = nel; 4909 else 4910 AddVolumeElement(nel); 4911 } 4912 4913 } 4914 } 4915 4916 4917 4918 4919 4920 else if (el.GetType() == PYRAMID) 4921 { 4922 // pyramid, to 2 tets 4923 4924 // cout << "pyramid: " << el << endl; 4925 4926 static const int ntets[2][8] = 4927 { { 1, 2, 3, 5, 1, 3, 4, 5 }, 4928 { 1, 2, 4, 5, 4, 2, 3, 5 }}; 4929 4930 const int * min2pi; 4931 4932 if (min2 (el[0], el[2]) < min2 (el[1], el[3])) 4933 min2pi = &ntets[0][0]; 4934 else 4935 min2pi = &ntets[1][0]; 4936 4937 bool firsttet = 1; 4938 for (int j = 0; j < 2; j++) 4939 { 4940 Element nel(TET); 4941 for (int k = 0; k < 4; k++) 4942 nel[k] = el[min2pi[4*j + k]-1]; 4943 nel.SetIndex (el.GetIndex()); 4944 4945 // cout << "pyramid-tet: " << nel << endl; 4946 4947 bool legal = 1; 4948 for (int k = 0; k < 3; k++) 4949 for (int l = k+1; l < 4; l++) 4950 if (nel[k] == nel[l]) 4951 legal = 0; 4952 4953 if (legal) 4954 { 4955 (*testout) << nel << " "; 4956 if (firsttet) 4957 VolumeElement(i) = nel; 4958 else 4959 AddVolumeElement(nel); 4960 4961 firsttet = 0; 4962 } 4963 } 4964 if (firsttet) cout << "no legal"; 4965 (*testout) << endl; 4966 } 4967 } 4968 4969 4970 int oldnse = GetNSE(); 4971 for (int i = 1; i <= oldnse; i++) 4972 { 4973 Element2d el = SurfaceElement(i); 4974 if (el.GetNP() == 4) 4975 { 4976 (*testout) << "split el: " << el << " to "; 4977 4978 static const int ntris[2][6] = 4979 { { 1, 2, 3, 1, 3, 4 }, 4980 { 1, 2, 4, 4, 2, 3 }}; 4981 4982 const int * min2pi; 4983 4984 if (min2 (el.PNum(1), el.PNum(3)) < 4985 min2 (el.PNum(2), el.PNum(4))) 4986 min2pi = &ntris[0][0]; 4987 else 4988 min2pi = &ntris[1][0]; 4989 4990 for (int j = 0; j <6; j++) 4991 (*testout) << min2pi[j] << " "; 4992 4993 4994 int firsttri = 1; 4995 for (int j = 1; j <= 2; j++) 4996 { 4997 Element2d nel(3); 4998 for (int k = 1; k <= 3; k++) 4999 nel.PNum(k) = el.PNum(min2pi[3 * j + k - 4]); 5000 nel.SetIndex (el.GetIndex()); 5001 5002 int legal = 1; 5003 for (int k = 1; k <= 2; k++) 5004 for (int l = k+1; l <= 3; l++) 5005 if (nel.PNum(k) == nel.PNum(l)) 5006 legal = 0; 5007 5008 if (legal) 5009 { 5010 (*testout) << nel << " "; 5011 if (firsttri) 5012 { 5013 SurfaceElement(i) = nel; 5014 firsttri = 0; 5015 } 5016 else 5017 { 5018 AddSurfaceElement(nel); 5019 } 5020 } 5021 } 5022 (*testout) << endl; 5023 5024 } 5025 } 5026 5027 5028 if (has_prisms) 5029 5030 Split2Tets(); 5031 5032 else 5033 { 5034 for (int i = 1; i <= GetNE(); i++) 5035 { 5036 Element & el = VolumeElement(i); 5037 const Point3d & p1 = Point (el.PNum(1)); 5038 const Point3d & p2 = Point (el.PNum(2)); 5039 const Point3d & p3 = Point (el.PNum(3)); 5040 const Point3d & p4 = Point (el.PNum(4)); 5041 5042 double vol = (Vec3d (p1, p2) * 5043 Cross (Vec3d (p1, p3), Vec3d(p1, p4))); 5044 if (vol > 0) 5045 swap (el.PNum(3), el.PNum(4)); 5046 } 5047 5048 5049 5050 UpdateTopology(); 5051 timestamp = NextTimeStamp(); 5052 } 5053 5054 RebuildSurfaceElementLists(); 5055 } 5056 BuildElementSearchTree()5057 void Mesh :: BuildElementSearchTree () 5058 { 5059 if (elementsearchtreets == GetTimeStamp()) return; 5060 5061 { 5062 std::lock_guard<std::mutex> guard(buildsearchtree_mutex); 5063 if (elementsearchtreets != GetTimeStamp()) 5064 { 5065 NgLock lock(mutex); 5066 lock.Lock(); 5067 5068 PrintMessage (4, "Rebuild element searchtree"); 5069 5070 elementsearchtree = nullptr; 5071 5072 int ne = (dimension == 2) ? GetNSE() : GetNE(); 5073 if (dimension == 3 && !GetNE() && GetNSE()) 5074 ne = GetNSE(); 5075 5076 if (ne) 5077 { 5078 if (dimension == 2 || (dimension == 3 && !GetNE()) ) 5079 { 5080 Box<3> box (Box<3>::EMPTY_BOX); 5081 for (SurfaceElementIndex sei = 0; sei < ne; sei++) 5082 // box.Add (points[surfelements[sei].PNums()]); 5083 for (auto pi : surfelements[sei].PNums()) 5084 box.Add (points[pi]); 5085 5086 box.Increase (1.01 * box.Diam()); 5087 elementsearchtree = make_unique<BoxTree<3>> (box); 5088 5089 for (SurfaceElementIndex sei = 0; sei < ne; sei++) 5090 { 5091 // box.Set (points[surfelements[sei].PNums()]); 5092 5093 Box<3> box (Box<3>::EMPTY_BOX); 5094 for (auto pi : surfelements[sei].PNums()) 5095 box.Add (points[pi]); 5096 5097 auto & el = surfelements[sei]; 5098 if(el.IsCurved() && curvedelems->IsSurfaceElementCurved(sei)) 5099 { 5100 netgen::Point<2> lami [4] = {netgen::Point<2>(0.5,0), netgen::Point<2>(0,0.5), netgen::Point<2>(0.5,0.5), netgen::Point<2>(1./3,1./3)}; 5101 for (auto lam : lami) 5102 { 5103 netgen::Point<3> x; 5104 Mat<3,2> Jac; 5105 5106 curvedelems->CalcSurfaceTransformation(lam,sei,x,Jac); 5107 box.Add (x); 5108 } 5109 box.Scale(1.2); 5110 } 5111 elementsearchtree -> Insert (box, sei+1); 5112 } 5113 } 5114 else 5115 { 5116 Box<3> box (Box<3>::EMPTY_BOX); 5117 for (ElementIndex ei = 0; ei < ne; ei++) 5118 // box.Add (points[volelements[ei].PNums()]); 5119 for (auto pi : volelements[ei].PNums()) 5120 box.Add (points[pi]); 5121 5122 box.Increase (1.01 * box.Diam()); 5123 elementsearchtree = make_unique<BoxTree<3>> (box); 5124 5125 for (ElementIndex ei = 0; ei < ne; ei++) 5126 { 5127 // box.Set (points[volelements[ei].PNums()]); 5128 5129 Box<3> box (Box<3>::EMPTY_BOX); 5130 for (auto pi : volelements[ei].PNums()) 5131 box.Add (points[pi]); 5132 5133 auto & el = volelements[ei]; 5134 if(el.IsCurved() && curvedelems->IsElementCurved(ei)) 5135 box.Scale(1.2); 5136 5137 5138 elementsearchtree -> Insert (box, ei+1); 5139 } 5140 } 5141 5142 elementsearchtreets = GetTimeStamp(); 5143 } 5144 } 5145 } 5146 } 5147 5148 SolveLinearSystemLS(const Vec3d & col1,const Vec3d & col2,const Vec3d & rhs,Vec2d & sol)5149 int SolveLinearSystemLS (const Vec3d & col1, 5150 const Vec3d & col2, 5151 const Vec3d & rhs, 5152 Vec2d & sol) 5153 { 5154 double a11 = col1 * col1; 5155 double a12 = col1 * col2; 5156 double a22 = col2 * col2; 5157 5158 double det = a11 * a22 - a12 * a12; 5159 5160 if (det*det <= 1e-24 * a11 * a22) 5161 { 5162 sol = Vec2d (0, 0); 5163 return 1; 5164 } 5165 5166 Vec2d aTrhs; 5167 aTrhs.X() = col1*rhs; 5168 aTrhs.Y() = col2*rhs; 5169 5170 sol.X() = ( a22 * aTrhs.X() - a12 * aTrhs.Y()) / det; 5171 sol.Y() = (-a12 * aTrhs.X() + a11 * aTrhs.Y()) / det; 5172 return 0; 5173 } 5174 ValidBarCoord(double lami[3],double eps=1e-12)5175 bool ValidBarCoord(double lami[3], double eps=1e-12) 5176 { 5177 return (lami[0]<=1.+eps && lami[0]>=0.-eps && lami[1]<=1.+eps && lami[1]>=0.-eps && lami[2]<=1.+eps && lami[2]>=0.-eps ); 5178 } 5179 PointContainedIn2DElement(const Point3d & p,double lami[3],const int element,bool consider3D) const5180 bool Mesh :: PointContainedIn2DElement(const Point3d & p, 5181 double lami[3], 5182 const int element, 5183 bool consider3D) const 5184 { 5185 Vec3d col1, col2, col3; 5186 Vec3d rhs, sol; 5187 const double eps = 1e-6; 5188 5189 NgArray<Element2d> loctrigs; 5190 5191 5192 //SZ 5193 if(SurfaceElement(element).GetType()==QUAD) 5194 { 5195 const Element2d & el = SurfaceElement(element); 5196 5197 const Point3d & p1 = Point(el.PNum(1)); 5198 const Point3d & p2 = Point(el.PNum(2)); 5199 const Point3d & p3 = Point(el.PNum(3)); 5200 const Point3d & p4 = Point(el.PNum(4)); 5201 5202 // Coefficients of Bilinear Mapping from Ref-Elem to global Elem 5203 // X = a + b x + c y + d x y 5204 Vec3d a = p1; 5205 Vec3d b = p2 - a; 5206 Vec3d c = p4 - a; 5207 Vec3d d = p3 - a - b - c; 5208 5209 /*cout << "p = " << p << endl; 5210 cout << "p1 = " << p1 << endl; 5211 cout << "p2 = " << p2 << endl; 5212 cout << "p3 = " << p3 << endl; 5213 cout << "p4 = " << p4 << endl; 5214 5215 cout << "a = " << a << endl; 5216 cout << "b = " << b << endl; 5217 cout << "c = " << c << endl; 5218 cout << "d = " << d << endl;*/ 5219 5220 5221 Vec3d pa = p-a; 5222 double dxb = d.X()*b.Y()-d.Y()*b.X(); 5223 double dxc = d.X()*c.Y()-d.Y()*c.X(); 5224 double bxc = b.X()*c.Y()-b.Y()*c.X(); 5225 double bxpa = b.X()*pa.Y()-b.Y()*pa.X(); 5226 double cxpa = c.X()*pa.Y()-c.Y()*pa.X(); 5227 double dxpa = d.X()*pa.Y()-d.Y()*pa.X(); 5228 5229 /*cout << "dxb = " << dxb << endl; 5230 cout << "dxc = " << dxc << endl; 5231 cout << "bxc = " << bxc << endl; 5232 cout << "bxpa = " << bxpa << endl; 5233 cout << "cxpa = " << cxpa << endl; 5234 cout << "dxpa = " << dxpa << endl;*/ 5235 5236 /* 5237 P = a + b x + c y + d x y 5238 1) P1 = a1 + b1 x + c1 y + d1 x y 5239 2) P2 = a2 + b2 x + c2 y + d2 x y 5240 5241 -> det(x,d) = det(a,d) + det(b,d) x + det(c,d) y 5242 -> x = 1/det(b,d) *( det(P-a,d)-det(c,d) y ) 5243 -> y = 1/det(c,d) *( det(P-a,d)-det(b,d) x ) 5244 5245 -> x = (P1 - a1 - c1 y)/(b1 + d1 y) 5246 -> det(c,d) y**2 + [det(d,P-a) + det(c,b)] y + det(b,P-a) = 0 5247 ( same if we express x = (P2 - a2 - c2 y)/(b2 + d2 y) ) 5248 5249 -> y = (P1 - a1 - b1 x)/(c1 + d1 x) 5250 -> det(b,d) x**2 + [det(d,P-a) + det(b,c)] x + det(c,P-a) = 0 5251 ( same if we express y = (P2 - a2 - b2 x)/(c2 + d2 x) 5252 */ 5253 5254 lami[2]=0.; 5255 double eps = 1.E-12; 5256 double c1,c2,r; 5257 5258 //First check if point is "exactly" a vertex point 5259 Vec3d d1 = p-p1; 5260 Vec3d d2 = p-p2; 5261 Vec3d d3 = p-p3; 5262 Vec3d d4 = p-p4; 5263 5264 //cout << " d1 = " << d1 << ", d2 = " << d2 << ", d3 = " << d3 << ", d4 = " << d4 << endl; 5265 5266 if (d1.Length2() < sqr(eps)*d2.Length2() && d1.Length2() < sqr(eps)*d3.Length2() && d1.Length2() < sqr(eps)*d4.Length2()) 5267 { 5268 lami[0] = lami[1] = 0.; 5269 return true; 5270 } 5271 else if (d2.Length2() < sqr(eps)*d1.Length2() && d2.Length2() < sqr(eps)*d3.Length2() && d2.Length2() < sqr(eps)*d4.Length2()) 5272 { 5273 lami[0] = 1.; 5274 lami[1] = 0.; 5275 return true; 5276 } 5277 else if (d3.Length2() < sqr(eps)*d1.Length2() && d3.Length2() < sqr(eps)*d2.Length2() && d3.Length2() < sqr(eps)*d4.Length2()) 5278 { 5279 lami[0] = lami[1] = 1.; 5280 return true; 5281 } 5282 else if (d4.Length2() < sqr(eps)*d1.Length2() && d4.Length2() < sqr(eps)*d2.Length2() && d4.Length2() < sqr(eps)*d3.Length2()) 5283 { 5284 lami[0] = 0.; 5285 lami[1] = 1.; 5286 return true; 5287 }//if d is nearly 0: solve resulting linear system 5288 else if (d.Length2() < sqr(eps)*b.Length2() && d.Length2() < sqr(eps)*c.Length2()) 5289 { 5290 Vec2d sol; 5291 SolveLinearSystemLS (b, c, p-a, sol); 5292 lami[0] = sol.X(); 5293 lami[1] = sol.Y(); 5294 return ValidBarCoord(lami, eps); 5295 }// if dxc is nearly 0: solve resulting linear equation for y and compute x 5296 else if (fabs(dxc) < sqr(eps)) 5297 { 5298 lami[1] = -bxpa/(dxpa-bxc); 5299 lami[0] = (dxpa-dxc*lami[1])/dxb; 5300 return ValidBarCoord(lami, eps); 5301 }// if dxb is nearly 0: solve resulting linear equation for x and compute y 5302 else if (fabs(dxb) < sqr(eps)) 5303 { 5304 lami[0] = -cxpa/(dxpa+bxc); 5305 lami[1] = (dxpa-dxb*lami[0])/dxc; 5306 return ValidBarCoord(lami, eps); 5307 }//if dxb >= dxc: solve quadratic equation in y and compute x 5308 else if (fabs(dxb) >= fabs(dxc)) 5309 { 5310 c1 = (bxc-dxpa)/dxc; 5311 c2 = -bxpa/dxc; 5312 r = c1*c1/4.0-c2; 5313 5314 //quadratic equation has only 1 (unstable) solution 5315 if (fabs(r) < eps) //not eps^2! 5316 { 5317 lami[1] = -c1/2; 5318 lami[0] = (dxpa-dxc*lami[1])/dxb; 5319 return ValidBarCoord(lami, eps); 5320 } 5321 if (r < 0) return false; 5322 5323 lami[1] = -c1/2+sqrt(r); 5324 lami[0] = (dxpa-dxc*lami[1])/dxb; 5325 5326 if (ValidBarCoord(lami, eps)) 5327 return true; 5328 else 5329 { 5330 lami[1] = -c1/2-sqrt(r); 5331 lami[0] = (dxpa-dxc*lami[1])/dxb; 5332 return ValidBarCoord(lami, eps); 5333 } 5334 }//if dxc > dxb: solve quadratic equation in x and compute y 5335 else 5336 { 5337 c1 = (-bxc-dxpa)/dxb; 5338 c2 = -cxpa/dxb; 5339 r = c1*c1/4.0-c2; 5340 5341 //quadratic equation has only 1 (unstable) solution 5342 if (fabs(r) < eps) //not eps^2! 5343 { 5344 lami[0] = -c1/2; 5345 lami[1] = (dxpa-dxb*lami[0])/dxc; 5346 return ValidBarCoord(lami, eps); 5347 } 5348 if (r < 0) return false; 5349 5350 lami[0] = -c1/2+sqrt(r); 5351 lami[1] = (dxpa-dxb*lami[0])/dxc; 5352 5353 if (ValidBarCoord(lami, eps)) 5354 return true; 5355 else 5356 { 5357 lami[0] = -c1/2-sqrt(r); 5358 lami[1] = (dxpa-dxb*lami[0])/dxc; 5359 return ValidBarCoord(lami, eps); 5360 } 5361 } 5362 5363 /* 5364 double dxa = d.X()*a.Y()-d.Y()*a.X(); 5365 double dxp = d.X()*p.Y()-d.Y()*p.X(); 5366 5367 5368 double c0,c1,c2; // ,rt; 5369 5370 5371 Vec3d dp13 = p3-p1; 5372 Vec3d dp24 = p4-p2; 5373 double d1 = dp13.Length2(); 5374 double d2 = dp24.Length2(); 5375 5376 // if(fabs(d.X()) <= eps && fabs(d.Y())<= eps) 5377 //if (d.Length2() < sqr(eps)) 5378 if (d.Length2() < sqr(eps)*d1 && d.Length2() < sqr(eps)*d2) 5379 { 5380 //Solve Linear System 5381 Vec2d sol; 5382 SolveLinearSystemLS (b, c, p-a, sol); 5383 lami[0] = sol.X(); 5384 lami[1] = sol.Y(); 5385 5386 if(lami[1]<=1.+eps && lami[1]>=0.-eps && lami[0]<=1.+eps && lami[0]>=0.-eps) 5387 return true; 5388 5389 5390 //lami[0]=(c.Y()*(p.X()-a.X())-c.X()*(p.Y()-a.Y()))/ 5391 //(b.X()*c.Y() -b.Y()*c.X()); 5392 //lami[1]=(-b.Y()*(p.X()-a.X())+b.X()*(p.Y()-a.Y()))/ 5393 // (b.X()*c.Y() -b.Y()*c.X()); 5394 5395 } 5396 else 5397 if(fabs(dxb) <= eps*fabs(dxc)) 5398 { 5399 lami[1] = (dxp-dxa)/dxc; 5400 if(fabs(b.X()+d.X()*lami[1])>=fabs(b.Y()+d.Y()*lami[1])) 5401 lami[0] = (p.X()-a.X() - c.X()*lami[1])/(b.X()+d.X()*lami[1]); 5402 else 5403 lami[0] = (p.Y()-a.Y() - c.Y()*lami[1])/(b.Y()+d.Y()*lami[1]); 5404 5405 if(lami[1]<=1.+eps && lami[1]>=0.-eps && lami[0]<=1.+eps && lami[0]>=0.-eps) 5406 return true; 5407 } 5408 else 5409 if(fabs(dxc) <= eps*fabs(dxb)) 5410 { 5411 lami[0] = (dxp-dxa)/dxb; 5412 if(fabs(c.X()+d.X()*lami[0])>=fabs(c.Y()+d.Y()*lami[0])) 5413 lami[1] = (p.X()-a.X() - b.X()*lami[0])/(c.X()+d.X()*lami[0]); 5414 else 5415 lami[1] = (p.Y()-a.Y() - b.Y()*lami[0])/(c.Y()+d.Y()*lami[0]); 5416 5417 if(lami[1]<=1.+eps && lami[1]>=0.-eps && lami[0]<=1.+eps && lami[0]>=0.-eps) 5418 return true; 5419 } 5420 else //Solve quadratic equation 5421 { 5422 c2 = -d.X()*dxb; 5423 c1 = b.X()*dxc - c.X()*dxb + d.X()*(dxp-dxa); 5424 c0 = c.X()*(dxp-dxa) + (a.X()-p.X())*dxc; 5425 double rt = c1*c1 - 4*c2*c0; 5426 5427 if (rt < 0.) return false; 5428 lami[1] = (-c1 + sqrt(rt))/2/c2; 5429 5430 5431 if(lami[1]<=1.+eps && lami[1]>=0.-eps) 5432 { 5433 lami[0] = (dxp - dxa -dxb*lami[1])/dxc; 5434 5435 if(lami[0]<=1.+eps && lami[0]>=0.-eps) 5436 return true; 5437 } 5438 lami[1] = (-c1 - sqrt(rt))/2/c2; 5439 5440 lami[0] = (dxp - dxa -dxb*lami[1])/dxc; 5441 5442 if(lami[1]<=1.+eps && lami[1]>=0.-eps && lami[0]<=1.+eps && lami[0]>=0.-eps) 5443 return true; 5444 5445 c2 = d.Y()*dxb; 5446 c1 = b.Y()*dxc - c.Y()*dxb + d.Y()*(dxp-dxa); 5447 c0 = c.Y()*(dxp -dxa) + (a.Y()-p.Y())*dxc; 5448 rt = c1*c1 - 4*c2*c0; 5449 5450 if (rt < 0.) return false; 5451 lami[1] = (-c1 + sqrt(rt))/2/c2; 5452 5453 if(lami[1]<=1.+eps && lami[1]>=0.-eps) 5454 { 5455 lami[0] = (dxp - dxa -dxb*lami[1])/dxc; 5456 5457 if(lami[0]<=1.+eps && lami[0]>=0.-eps) 5458 return true; 5459 } 5460 lami[1] = (-c1 - sqrt(rt))/2/c2; 5461 5462 lami[0] = (dxp - dxa -dxb*lami[1])/dxc; 5463 5464 if(lami[1]<=1.+eps && lami[1]>=0.-eps && lami[0]<=1.+eps && lami[0]>=0.-eps) 5465 return true; 5466 5467 c2 = -d.X()*dxc; 5468 c1 = -b.X()*dxc + c.X()*dxb + d.X()*(dxp-dxa); 5469 c0 = b.X()*(dxp -dxa) + (a.X()-p.X())*dxb; 5470 rt = c1*c1 - 4*c2*c0; 5471 5472 if (rt < 0.) return false; 5473 lami[1] = (-c1 + sqrt(rt))/2/c2; 5474 5475 if(lami[1]<=1.+eps && lami[1]>=0.-eps) 5476 { 5477 lami[0] = (dxp - dxa -dxc*lami[1])/dxb; 5478 5479 if(lami[0]<=1.+eps && lami[0]>=0.-eps) 5480 return true; 5481 } 5482 lami[1] = (-c1 - sqrt(rt))/2/c2; 5483 5484 lami[0] = (dxp - dxa -dxc*lami[1])/dxb; 5485 5486 if(lami[1]<=1.+eps && lami[1]>=0.-eps && lami[0]<=1.+eps && lami[0]>=0.-eps) 5487 return true; 5488 }*/ 5489 5490 5491 //cout << "lam0,1 = " << lami[0] << ", " << lami[1] << endl; 5492 5493 /*if( lami[0] <= 1.+eps && lami[0] >= -eps && lami[1]<=1.+eps && lami[1]>=-eps) 5494 { 5495 if(consider3D) 5496 { 5497 Vec3d n = Cross(b,c); 5498 lami[2] = 0; 5499 for(int i=1; i<=3; i++) 5500 lami[2] +=(p.X(i)-a.X(i)-lami[0]*b.X(i)-lami[1]*c.X(i)) * n.X(i); 5501 if(lami[2] >= -eps && lami[2] <= eps) 5502 return true; 5503 } 5504 else 5505 return true; 5506 }*/ 5507 5508 return false; 5509 5510 } 5511 else 5512 { 5513 // SurfaceElement(element).GetTets (loctets); 5514 loctrigs.SetSize(1); 5515 loctrigs.Elem(1) = SurfaceElement(element); 5516 5517 5518 5519 for (int j = 1; j <= loctrigs.Size(); j++) 5520 { 5521 const Element2d & el = loctrigs.Get(j); 5522 5523 5524 const Point3d & p1 = Point(el.PNum(1)); 5525 const Point3d & p2 = Point(el.PNum(2)); 5526 const Point3d & p3 = Point(el.PNum(3)); 5527 /* 5528 Box3d box; 5529 box.SetPoint (p1); 5530 box.AddPoint (p2); 5531 box.AddPoint (p3); 5532 box.AddPoint (p4); 5533 if (!box.IsIn (p)) 5534 continue; 5535 */ 5536 col1 = p2-p1; 5537 col2 = p3-p1; 5538 col3 = Cross(col1,col2); 5539 //col3 = Vec3d(0, 0, 1); 5540 rhs = p - p1; 5541 5542 // int retval = 5543 SolveLinearSystem (col1, col2, col3, rhs, sol); 5544 5545 //(*testout) << "retval " << retval << endl; 5546 5547 //(*testout) << "col1 " << col1 << " col2 " << col2 << " col3 " << col3 << " rhs " << rhs << endl; 5548 //(*testout) << "sol " << sol << endl; 5549 5550 if (SurfaceElement(element).GetType() ==TRIG6 || curvedelems->IsSurfaceElementCurved(element-1)) 5551 { 5552 netgen::Point<2> lam(1./3,1./3); 5553 Vec<3> rhs; 5554 Vec<2> deltalam; 5555 netgen::Point<3> x; 5556 Mat<3,2> Jac,Jact; 5557 5558 double delta=1; 5559 5560 bool retval; 5561 5562 int i = 0; 5563 5564 const int maxits = 30; 5565 while(delta > 1e-16 && i<maxits) 5566 { 5567 curvedelems->CalcSurfaceTransformation(lam,element-1,x,Jac); 5568 rhs = p-x; 5569 Jac.Solve(rhs,deltalam); 5570 5571 lam += deltalam; 5572 5573 delta = deltalam.Length2(); 5574 5575 i++; 5576 //(*testout) << "pcie i " << i << " delta " << delta << " p " << p << " x " << x << " lam " << lam << endl; 5577 //<< "Jac " << Jac << endl; 5578 } 5579 5580 if(i==maxits) 5581 return false; 5582 5583 sol.X() = lam(0); 5584 sol.Y() = lam(1); 5585 5586 if (SurfaceElement(element).GetType() !=TRIG6 ) 5587 { 5588 sol.Z() = sol.X(); 5589 sol.X() = sol.Y(); 5590 sol.Y() = 1.0 - sol.Z() - sol.X(); 5591 } 5592 5593 } 5594 if (sol.X() >= -eps && sol.Y() >= -eps && 5595 sol.X() + sol.Y() <= 1+eps) 5596 { 5597 if(!consider3D || (sol.Z() >= -eps && sol.Z() <= eps)) 5598 { 5599 lami[0] = sol.X(); 5600 lami[1] = sol.Y(); 5601 lami[2] = sol.Z(); 5602 5603 return true; 5604 } 5605 } 5606 } 5607 } 5608 5609 return false; 5610 5611 } 5612 5613 5614 5615 PointContainedIn3DElement(const Point3d & p,double lami[3],const int element) const5616 bool Mesh :: PointContainedIn3DElement(const Point3d & p, 5617 double lami[3], 5618 const int element) const 5619 { 5620 //bool oldresult = PointContainedIn3DElementOld(p,lami,element); 5621 //(*testout) << "old result: " << oldresult 5622 // << " lam " << lami[0] << " " << lami[1] << " " << lami[2] << endl; 5623 5624 //if(!curvedelems->IsElementCurved(element-1)) 5625 // return PointContainedIn3DElementOld(p,lami,element); 5626 5627 5628 const double eps = 1.e-4; 5629 const Element & el = VolumeElement(element); 5630 5631 netgen::Point<3> lam = 0.0; 5632 5633 if (el.GetType() == TET || el.GetType() == TET10) 5634 { 5635 lam = 0.25; 5636 } 5637 else if (el.GetType() == PRISM) 5638 { 5639 lam(0) = 0.33; lam(1) = 0.33; lam(2) = 0.5; 5640 } 5641 else if (el.GetType() == PYRAMID) 5642 { 5643 lam(0) = 0.4; lam(1) = 0.4; lam(2) = 0.2; 5644 } 5645 else if (el.GetType() == HEX) 5646 { 5647 lam = 0.5; 5648 } 5649 5650 5651 Vec<3> deltalam,rhs; 5652 netgen::Point<3> x; 5653 Mat<3,3> Jac,Jact; 5654 5655 double delta=1; 5656 5657 bool retval; 5658 5659 int i = 0; 5660 5661 const int maxits = 30; 5662 while(delta > 1e-16 && i<maxits) 5663 { 5664 curvedelems->CalcElementTransformation(lam,element-1,x,Jac); 5665 rhs = p-x; 5666 Jac.Solve(rhs,deltalam); 5667 5668 lam += deltalam; 5669 5670 delta = deltalam.Length2(); 5671 5672 i++; 5673 //(*testout) << "pcie i " << i << " delta " << delta << " p " << p << " x " << x << " lam " << lam << endl; 5674 //<< "Jac " << Jac << endl; 5675 } 5676 5677 if(i==maxits) 5678 return false; 5679 5680 for(i=0; i<3; i++) 5681 lami[i] = lam(i); 5682 5683 5684 5685 if (el.GetType() == TET || el.GetType() == TET10) 5686 { 5687 retval = (lam(0) > -eps && 5688 lam(1) > -eps && 5689 lam(2) > -eps && 5690 lam(0) + lam(1) + lam(2) < 1+eps); 5691 } 5692 else if (el.GetType() == PRISM || el.GetType() == PRISM15) 5693 { 5694 retval = (lam(0) > -eps && 5695 lam(1) > -eps && 5696 lam(2) > -eps && 5697 lam(2) < 1+eps && 5698 lam(0) + lam(1) < 1+eps); 5699 } 5700 else if (el.GetType() == PYRAMID || el.GetType() == PYRAMID13) 5701 { 5702 retval = (lam(0) > -eps && 5703 lam(1) > -eps && 5704 lam(2) > -eps && 5705 lam(0) + lam(2) < 1+eps && 5706 lam(1) + lam(2) < 1+eps); 5707 } 5708 else if (el.GetType() == HEX || el.GetType() == HEX20) 5709 { 5710 retval = (lam(0) > -eps && lam(0) < 1+eps && 5711 lam(1) > -eps && lam(1) < 1+eps && 5712 lam(2) > -eps && lam(2) < 1+eps); 5713 } 5714 else 5715 throw NgException("Da haun i wos vagessn"); 5716 5717 return retval; 5718 } 5719 5720 5721 PointContainedIn3DElementOld(const Point3d & p,double lami[3],const int element) const5722 bool Mesh :: PointContainedIn3DElementOld(const Point3d & p, 5723 double lami[3], 5724 const int element) const 5725 { 5726 Vec3d col1, col2, col3; 5727 Vec3d rhs, sol; 5728 const double eps = 1.e-4; 5729 5730 NgArray<Element> loctets; 5731 5732 VolumeElement(element).GetTets (loctets); 5733 5734 for (int j = 1; j <= loctets.Size(); j++) 5735 { 5736 const Element & el = loctets.Get(j); 5737 5738 const Point3d & p1 = Point(el.PNum(1)); 5739 const Point3d & p2 = Point(el.PNum(2)); 5740 const Point3d & p3 = Point(el.PNum(3)); 5741 const Point3d & p4 = Point(el.PNum(4)); 5742 5743 Box3d box; 5744 box.SetPoint (p1); 5745 box.AddPoint (p2); 5746 box.AddPoint (p3); 5747 box.AddPoint (p4); 5748 if (!box.IsIn (p)) 5749 continue; 5750 5751 col1 = p2-p1; 5752 col2 = p3-p1; 5753 col3 = p4-p1; 5754 rhs = p - p1; 5755 5756 SolveLinearSystem (col1, col2, col3, rhs, sol); 5757 5758 if (sol.X() >= -eps && sol.Y() >= -eps && sol.Z() >= -eps && 5759 sol.X() + sol.Y() + sol.Z() <= 1+eps) 5760 { 5761 NgArray<Element> loctetsloc; 5762 NgArray<netgen::Point<3> > pointsloc; 5763 5764 VolumeElement(element).GetTetsLocal (loctetsloc); 5765 VolumeElement(element).GetNodesLocalNew (pointsloc); 5766 5767 const Element & le = loctetsloc.Get(j); 5768 5769 5770 Point3d pp = 5771 pointsloc.Get(le.PNum(1)) 5772 + sol.X() * Vec3d (pointsloc.Get(le.PNum(1)), pointsloc.Get(le.PNum(2))) 5773 + sol.Y() * Vec3d (pointsloc.Get(le.PNum(1)), pointsloc.Get(le.PNum(3))) 5774 + sol.Z() * Vec3d (pointsloc.Get(le.PNum(1)), pointsloc.Get(le.PNum(4))) ; 5775 5776 lami[0] = pp.X(); 5777 lami[1] = pp.Y(); 5778 lami[2] = pp.Z(); 5779 return true; 5780 } 5781 } 5782 return false; 5783 } 5784 5785 GetElementOfPoint(const netgen::Point<3> & p,double lami[3],bool build_searchtree,const int index,const bool allowindex) const5786 int Mesh :: GetElementOfPoint (const netgen::Point<3> & p, 5787 double lami[3], 5788 bool build_searchtree, 5789 const int index, 5790 const bool allowindex) const 5791 { 5792 if(index != -1) 5793 { 5794 NgArray<int> dummy(1); 5795 dummy[0] = index; 5796 return GetElementOfPoint(p,lami,&dummy,build_searchtree,allowindex); 5797 } 5798 else 5799 return GetElementOfPoint(p,lami,NULL,build_searchtree,allowindex); 5800 } 5801 5802 5803 5804 GetElementOfPoint(const netgen::Point<3> & p,double lami[3],const NgArray<int> * const indices,bool build_searchtree,const bool allowindex) const5805 int Mesh :: GetElementOfPoint (const netgen::Point<3> & p, 5806 double lami[3], 5807 const NgArray<int> * const indices, 5808 bool build_searchtree, 5809 const bool allowindex) const 5810 { 5811 if ( (dimension == 2 && !GetNSE()) || 5812 (dimension == 3 && !GetNE() && !GetNSE()) ) 5813 return -1; 5814 5815 if (build_searchtree) 5816 const_cast<Mesh&>(*this).BuildElementSearchTree (); 5817 5818 if (dimension == 2 || (dimension==3 && !GetNE() && GetNSE())) 5819 return Find2dElement(*this, p, lami, indices, elementsearchtree.get(), allowindex); 5820 5821 return Find3dElement(*this, p, lami, indices, elementsearchtree.get(), allowindex); 5822 } 5823 5824 5825 GetSurfaceElementOfPoint(const netgen::Point<3> & p,double lami[3],bool build_searchtree,const int index,const bool allowindex) const5826 int Mesh :: GetSurfaceElementOfPoint (const netgen::Point<3> & p, 5827 double lami[3], 5828 bool build_searchtree, 5829 const int index, 5830 const bool allowindex) const 5831 { 5832 if(index != -1) 5833 { 5834 NgArray<int> dummy(1); 5835 dummy[0] = index; 5836 return GetSurfaceElementOfPoint(p,lami,&dummy,build_searchtree,allowindex); 5837 } 5838 else 5839 return GetSurfaceElementOfPoint(p,lami,NULL,build_searchtree,allowindex); 5840 } 5841 5842 5843 5844 GetSurfaceElementOfPoint(const netgen::Point<3> & p,double lami[3],const NgArray<int> * const indices,bool build_searchtree,const bool allowindex) const5845 int Mesh :: GetSurfaceElementOfPoint (const netgen::Point<3> & p, 5846 double lami[3], 5847 const NgArray<int> * const indices, 5848 bool build_searchtree, 5849 const bool allowindex) const 5850 { 5851 if (!GetNE() && build_searchtree) 5852 const_cast<Mesh&>(*this).BuildElementSearchTree (); 5853 5854 if (dimension == 2) 5855 return Find1dElement(*this, p, lami, indices, elementsearchtree.get(), allowindex); 5856 else 5857 return Find2dElement(*this, p, lami, indices, elementsearchtree.get(), allowindex); 5858 return 0; 5859 } 5860 5861 GetIntersectingVolEls(const Point3d & p1,const Point3d & p2,NgArray<int> & locels) const5862 void Mesh::GetIntersectingVolEls(const Point3d& p1, const Point3d& p2, 5863 NgArray<int> & locels) const 5864 { 5865 elementsearchtree->GetIntersecting (p1, p2, locels); 5866 } 5867 SplitIntoParts()5868 void Mesh :: SplitIntoParts() 5869 { 5870 int i, j, dom; 5871 int ne = GetNE(); 5872 int np = GetNP(); 5873 int nse = GetNSE(); 5874 5875 NgBitArray surfused(nse); 5876 NgBitArray pused (np); 5877 5878 surfused.Clear(); 5879 5880 dom = 0; 5881 5882 while (1) 5883 { 5884 int cntd = 1; 5885 5886 dom++; 5887 5888 pused.Clear(); 5889 5890 int found = 0; 5891 for (i = 1; i <= nse; i++) 5892 if (!surfused.Test(i)) 5893 { 5894 SurfaceElement(i).SetIndex (dom); 5895 for (j = 1; j <= 3; j++) 5896 pused.Set (SurfaceElement(i).PNum(j)); 5897 found = 1; 5898 cntd = 1; 5899 surfused.Set(i); 5900 break; 5901 } 5902 5903 if (!found) 5904 break; 5905 5906 int change; 5907 do 5908 { 5909 change = 0; 5910 for (i = 1; i <= nse; i++) 5911 { 5912 int is = 0, isnot = 0; 5913 for (j = 1; j <= 3; j++) 5914 if (pused.Test(SurfaceElement(i).PNum(j))) 5915 is = 1; 5916 else 5917 isnot = 1; 5918 5919 if (is && isnot) 5920 { 5921 change = 1; 5922 for (j = 1; j <= 3; j++) 5923 pused.Set (SurfaceElement(i).PNum(j)); 5924 } 5925 5926 if (is) 5927 { 5928 if (!surfused.Test(i)) 5929 { 5930 surfused.Set(i); 5931 SurfaceElement(i).SetIndex (dom); 5932 cntd++; 5933 } 5934 } 5935 } 5936 5937 5938 for (i = 1; i <= ne; i++) 5939 { 5940 int is = 0, isnot = 0; 5941 for (j = 1; j <= 4; j++) 5942 if (pused.Test(VolumeElement(i).PNum(j))) 5943 is = 1; 5944 else 5945 isnot = 1; 5946 5947 if (is && isnot) 5948 { 5949 change = 1; 5950 for (j = 1; j <= 4; j++) 5951 pused.Set (VolumeElement(i).PNum(j)); 5952 } 5953 5954 if (is) 5955 { 5956 VolumeElement(i).SetIndex (dom); 5957 } 5958 } 5959 } 5960 while (change); 5961 5962 PrintMessage (3, "domain ", dom, " has ", cntd, " surfaceelements"); 5963 } 5964 5965 /* 5966 facedecoding.SetSize (dom); 5967 for (i = 1; i <= dom; i++) 5968 { 5969 facedecoding.Elem(i).surfnr = 0; 5970 facedecoding.Elem(i).domin = i; 5971 facedecoding.Elem(i).domout = 0; 5972 } 5973 */ 5974 ClearFaceDescriptors(); 5975 for (i = 1; i <= dom; i++) 5976 AddFaceDescriptor (FaceDescriptor (0, i, 0, 0)); 5977 CalcSurfacesOfNode(); 5978 timestamp = NextTimeStamp(); 5979 } 5980 SplitSeparatedFaces()5981 void Mesh :: SplitSeparatedFaces () 5982 { 5983 PrintMessage (3, "SplitSeparateFaces"); 5984 int fdi; 5985 int np = GetNP(); 5986 5987 NgBitArray usedp(np); 5988 Array<SurfaceElementIndex> els_of_face; 5989 5990 fdi = 1; 5991 while (fdi <= GetNFD()) 5992 { 5993 GetSurfaceElementsOfFace (fdi, els_of_face); 5994 5995 if (els_of_face.Size() == 0) continue; 5996 5997 SurfaceElementIndex firstel = els_of_face[0]; 5998 5999 usedp.Clear(); 6000 for (int j = 1; j <= SurfaceElement(firstel).GetNP(); j++) 6001 usedp.Set (SurfaceElement(firstel).PNum(j)); 6002 6003 bool changed; 6004 do 6005 { 6006 changed = false; 6007 6008 for (int i = 0; i < els_of_face.Size(); i++) 6009 { 6010 const Element2d & el = SurfaceElement(els_of_face[i]); 6011 6012 bool has = 0; 6013 bool hasno = 0; 6014 for (int j = 0; j < el.GetNP(); j++) 6015 { 6016 if (usedp.Test(el[j])) 6017 has = true; 6018 else 6019 hasno = true; 6020 } 6021 6022 if (has && hasno) 6023 changed = true; 6024 6025 if (has) 6026 for (int j = 0; j < el.GetNP(); j++) 6027 usedp.Set (el[j]); 6028 } 6029 } 6030 while (changed); 6031 6032 int nface = 0; 6033 for (int i = 0; i < els_of_face.Size(); i++) 6034 { 6035 Element2d & el = SurfaceElement(els_of_face[i]); 6036 6037 int hasno = 0; 6038 for (int j = 1; j <= el.GetNP(); j++) 6039 if (!usedp.Test(el.PNum(j))) 6040 hasno = 1; 6041 6042 if (hasno) 6043 { 6044 if (!nface) 6045 { 6046 FaceDescriptor nfd = GetFaceDescriptor(fdi); 6047 nface = AddFaceDescriptor (nfd); 6048 } 6049 6050 el.SetIndex (nface); 6051 } 6052 } 6053 6054 // reconnect list 6055 if (nface) 6056 { 6057 facedecoding[nface-1].firstelement = -1; 6058 facedecoding[fdi-1].firstelement = -1; 6059 6060 for (int i = 0; i < els_of_face.Size(); i++) 6061 { 6062 int ind = SurfaceElement(els_of_face[i]).GetIndex(); 6063 SurfaceElement(els_of_face[i]).next = facedecoding[ind-1].firstelement; 6064 facedecoding[ind-1].firstelement = els_of_face[i]; 6065 } 6066 6067 // map the segments 6068 for(auto& seg : segments) 6069 if(!usedp.Test(seg[0]) || !usedp.Test(seg[1])) 6070 if(seg.si == fdi) 6071 seg.si = nface; 6072 } 6073 6074 fdi++; 6075 } 6076 6077 6078 /* 6079 fdi = 1; 6080 while (fdi <= GetNFD()) 6081 { 6082 int firstel = 0; 6083 for (int i = 1; i <= GetNSE(); i++) 6084 if (SurfaceElement(i).GetIndex() == fdi) 6085 { 6086 firstel = i; 6087 break; 6088 } 6089 if (!firstel) continue; 6090 6091 usedp.Clear(); 6092 for (int j = 1; j <= SurfaceElement(firstel).GetNP(); j++) 6093 usedp.Set (SurfaceElement(firstel).PNum(j)); 6094 6095 int changed; 6096 do 6097 { 6098 changed = 0; 6099 for (int i = 1; i <= GetNSE(); i++) 6100 { 6101 const Element2d & el = SurfaceElement(i); 6102 if (el.GetIndex() != fdi) 6103 continue; 6104 6105 int has = 0; 6106 int hasno = 0; 6107 for (int j = 1; j <= el.GetNP(); j++) 6108 { 6109 if (usedp.Test(el.PNum(j))) 6110 has = 1; 6111 else 6112 hasno = 1; 6113 } 6114 if (has && hasno) 6115 changed = 1; 6116 6117 if (has) 6118 for (int j = 1; j <= el.GetNP(); j++) 6119 usedp.Set (el.PNum(j)); 6120 } 6121 } 6122 while (changed); 6123 6124 int nface = 0; 6125 for (int i = 1; i <= GetNSE(); i++) 6126 { 6127 Element2d & el = SurfaceElement(i); 6128 if (el.GetIndex() != fdi) 6129 continue; 6130 6131 int hasno = 0; 6132 for (int j = 1; j <= el.GetNP(); j++) 6133 { 6134 if (!usedp.Test(el.PNum(j))) 6135 hasno = 1; 6136 } 6137 6138 if (hasno) 6139 { 6140 if (!nface) 6141 { 6142 FaceDescriptor nfd = GetFaceDescriptor(fdi); 6143 nface = AddFaceDescriptor (nfd); 6144 } 6145 6146 el.SetIndex (nface); 6147 } 6148 } 6149 fdi++; 6150 } 6151 */ 6152 } 6153 6154 6155 RebuildSurfaceElementLists()6156 void Mesh :: RebuildSurfaceElementLists () 6157 { 6158 static Timer t("Mesh::LinkSurfaceElements"); RegionTimer reg (t); 6159 6160 for (int i = 0; i < facedecoding.Size(); i++) 6161 facedecoding[i].firstelement = -1; 6162 for (int i = surfelements.Size()-1; i >= 0; i--) 6163 { 6164 int ind = surfelements[i].GetIndex(); 6165 surfelements[i].next = facedecoding[ind-1].firstelement; 6166 facedecoding[ind-1].firstelement = i; 6167 } 6168 } 6169 GetSurfaceElementsOfFace(int facenr,Array<SurfaceElementIndex> & sei) const6170 void Mesh :: GetSurfaceElementsOfFace (int facenr, Array<SurfaceElementIndex> & sei) const 6171 { 6172 static int timer = NgProfiler::CreateTimer ("GetSurfaceElementsOfFace"); 6173 NgProfiler::RegionTimer reg (timer); 6174 6175 /* 6176 sei.SetSize (0); 6177 for (SurfaceElementIndex i = 0; i < GetNSE(); i++) 6178 { 6179 if ( (*this)[i].GetIndex () == facenr && (*this)[i][0] >= PointIndex::BASE && 6180 !(*this)[i].IsDeleted() ) 6181 { 6182 sei.Append (i); 6183 } 6184 } 6185 */ 6186 6187 /* Philippose - 01/10/2009 6188 Commented out the following lines, and activated the originally 6189 commented out lines above because of a bug which causes corruption 6190 of the variable "facedecoding" when a mesh is converted to second order 6191 */ 6192 6193 // int size1 = sei.Size(); 6194 sei.SetSize(0); 6195 6196 SurfaceElementIndex si = facedecoding[facenr-1].firstelement; 6197 while (si != -1) 6198 { 6199 if ( (*this)[si].GetIndex () == facenr && (*this)[si][0] >= PointIndex::BASE && 6200 !(*this)[si].IsDeleted() ) 6201 { 6202 sei.Append (si); 6203 } 6204 6205 si = (*this)[si].next; 6206 } 6207 6208 /* 6209 // *testout << "with list = " << endl << sei << endl; 6210 6211 if (size1 != sei.Size()) 6212 { 6213 cout << "size mismatch" << endl; 6214 exit(1); 6215 } 6216 */ 6217 } 6218 6219 6220 6221 CalcMinMaxAngle(double badellimit,double * retvalues)6222 void Mesh :: CalcMinMaxAngle (double badellimit, double * retvalues) 6223 { 6224 int i, j; 6225 int lpi1, lpi2, lpi3, lpi4; 6226 double phimax = 0, phimin = 10; 6227 double facephimax = 0, facephimin = 10; 6228 int illegaltets = 0, negativetets = 0, badtets = 0; 6229 6230 for (i = 1; i <= GetNE(); i++) 6231 { 6232 int badel = 0; 6233 6234 Element & el = VolumeElement(i); 6235 6236 if (el.GetType() != TET) 6237 { 6238 VolumeElement(i).flags.badel = 0; 6239 continue; 6240 } 6241 6242 if (el.Volume(Points()) < 0) 6243 { 6244 badel = 1; 6245 negativetets++; 6246 } 6247 6248 6249 if (!LegalTet (el)) 6250 { 6251 badel = 1; 6252 illegaltets++; 6253 (*testout) << "illegal tet: " << i << " "; 6254 for (j = 1; j <= el.GetNP(); j++) 6255 (*testout) << el.PNum(j) << " "; 6256 (*testout) << endl; 6257 } 6258 6259 6260 // angles between faces 6261 for (lpi1 = 1; lpi1 <= 3; lpi1++) 6262 for (lpi2 = lpi1+1; lpi2 <= 4; lpi2++) 6263 { 6264 lpi3 = 1; 6265 while (lpi3 == lpi1 || lpi3 == lpi2) 6266 lpi3++; 6267 lpi4 = 10 - lpi1 - lpi2 - lpi3; 6268 6269 const Point3d & p1 = Point (el.PNum(lpi1)); 6270 const Point3d & p2 = Point (el.PNum(lpi2)); 6271 const Point3d & p3 = Point (el.PNum(lpi3)); 6272 const Point3d & p4 = Point (el.PNum(lpi4)); 6273 6274 Vec3d n(p1, p2); 6275 n /= n.Length(); 6276 Vec3d v1(p1, p3); 6277 Vec3d v2(p1, p4); 6278 6279 v1 -= (n * v1) * n; 6280 v2 -= (n * v2) * n; 6281 6282 double cosphi = (v1 * v2) / (v1.Length() * v2.Length()); 6283 double phi = acos (cosphi); 6284 if (phi > phimax) phimax = phi; 6285 if (phi < phimin) phimin = phi; 6286 6287 if ((180/M_PI) * phi > badellimit) 6288 badel = 1; 6289 } 6290 6291 6292 // angles in faces 6293 for (j = 1; j <= 4; j++) 6294 { 6295 Element2d face(TRIG); 6296 el.GetFace (j, face); 6297 for (lpi1 = 1; lpi1 <= 3; lpi1++) 6298 { 6299 lpi2 = lpi1 % 3 + 1; 6300 lpi3 = lpi2 % 3 + 1; 6301 6302 const Point3d & p1 = Point (el.PNum(lpi1)); 6303 const Point3d & p2 = Point (el.PNum(lpi2)); 6304 const Point3d & p3 = Point (el.PNum(lpi3)); 6305 6306 Vec3d v1(p1, p2); 6307 Vec3d v2(p1, p3); 6308 double cosphi = (v1 * v2) / (v1.Length() * v2.Length()); 6309 double phi = acos (cosphi); 6310 if (phi > facephimax) facephimax = phi; 6311 if (phi < facephimin) facephimin = phi; 6312 6313 if ((180/M_PI) * phi > badellimit) 6314 badel = 1; 6315 6316 } 6317 } 6318 6319 6320 VolumeElement(i).flags.badel = badel; 6321 if (badel) badtets++; 6322 } 6323 6324 if (!GetNE()) 6325 { 6326 phimin = phimax = facephimin = facephimax = 0; 6327 } 6328 6329 if (!retvalues) 6330 { 6331 PrintMessage (1, ""); 6332 PrintMessage (1, "between planes: phimin = ", (180/M_PI) * phimin, 6333 " phimax = ", (180/M_PI) *phimax); 6334 PrintMessage (1, "inside planes: phimin = ", (180/M_PI) * facephimin, 6335 " phimax = ", (180/M_PI) * facephimax); 6336 PrintMessage (1, ""); 6337 } 6338 else 6339 { 6340 retvalues[0] = (180/M_PI) * facephimin; 6341 retvalues[1] = (180/M_PI) * facephimax; 6342 retvalues[2] = (180/M_PI) * phimin; 6343 retvalues[3] = (180/M_PI) * phimax; 6344 } 6345 PrintMessage (3, "negative tets: ", negativetets); 6346 PrintMessage (3, "illegal tets: ", illegaltets); 6347 PrintMessage (3, "bad tets: ", badtets); 6348 } 6349 6350 MarkIllegalElements()6351 int Mesh :: MarkIllegalElements () 6352 { 6353 if(!boundaryedges) 6354 BuildBoundaryEdges(); 6355 6356 atomic<int> cnt = 0; 6357 ParallelForRange( Range(volelements), [&] (auto myrange) 6358 { 6359 int cnt_local = 0; 6360 for(auto & el : volelements.Range(myrange)) 6361 if (!LegalTet (el)) 6362 cnt_local++; 6363 cnt += cnt_local; 6364 }); 6365 return cnt; 6366 } 6367 6368 // #ifdef NONE 6369 // void Mesh :: AddIdentification (int pi1, int pi2, int identnr) 6370 // { 6371 // INDEX_2 pair(pi1, pi2); 6372 // // pair.Sort(); 6373 // identifiedpoints->Set (pair, identnr); 6374 // if (identnr > maxidentnr) 6375 // maxidentnr = identnr; 6376 // timestamp = NextTimeStamp(); 6377 // } 6378 6379 // int Mesh :: GetIdentification (int pi1, int pi2) const 6380 // { 6381 // INDEX_2 pair(pi1, pi2); 6382 // if (identifiedpoints->Used (pair)) 6383 // return identifiedpoints->Get(pair); 6384 // else 6385 // return 0; 6386 // } 6387 6388 // int Mesh :: GetIdentificationSym (int pi1, int pi2) const 6389 // { 6390 // INDEX_2 pair(pi1, pi2); 6391 // if (identifiedpoints->Used (pair)) 6392 // return identifiedpoints->Get(pair); 6393 6394 // pair = INDEX_2 (pi2, pi1); 6395 // if (identifiedpoints->Used (pair)) 6396 // return identifiedpoints->Get(pair); 6397 6398 // return 0; 6399 // } 6400 6401 6402 // void Mesh :: GetIdentificationMap (int identnr, NgArray<int> & identmap) const 6403 // { 6404 // int i, j; 6405 6406 // identmap.SetSize (GetNP()); 6407 // for (i = 1; i <= identmap.Size(); i++) 6408 // identmap.Elem(i) = 0; 6409 6410 // for (i = 1; i <= identifiedpoints->GetNBags(); i++) 6411 // for (j = 1; j <= identifiedpoints->GetBagSize(i); j++) 6412 // { 6413 // INDEX_2 i2; 6414 // int nr; 6415 // identifiedpoints->GetData (i, j, i2, nr); 6416 6417 // if (nr == identnr) 6418 // { 6419 // identmap.Elem(i2.I1()) = i2.I2(); 6420 // } 6421 // } 6422 // } 6423 6424 6425 // void Mesh :: GetIdentificationPairs (int identnr, NgArray<INDEX_2> & identpairs) const 6426 // { 6427 // int i, j; 6428 6429 // identpairs.SetSize(0); 6430 6431 // for (i = 1; i <= identifiedpoints->GetNBags(); i++) 6432 // for (j = 1; j <= identifiedpoints->GetBagSize(i); j++) 6433 // { 6434 // INDEX_2 i2; 6435 // int nr; 6436 // identifiedpoints->GetData (i, j, i2, nr); 6437 6438 // if (identnr == 0 || nr == identnr) 6439 // identpairs.Append (i2); 6440 // } 6441 // } 6442 // #endif 6443 IdentifyPeriodicBoundaries(const string & s1,const string & s2,const Transformation<3> & mapping,double pointTolerance)6444 int Mesh::IdentifyPeriodicBoundaries(const string &s1, 6445 const string &s2, 6446 const Transformation<3> &mapping, 6447 double pointTolerance) 6448 { 6449 auto nr = ident->GetMaxNr() + 1; 6450 ident->SetType(nr, Identifications::PERIODIC); 6451 double lami[4]; 6452 set<int> identified_points; 6453 if(pointTolerance < 0.) 6454 { 6455 Point3d pmin, pmax; 6456 GetBox(pmin, pmax); 6457 pointTolerance = 1e-8 * (pmax-pmin).Length(); 6458 } 6459 for(const auto& se : surfelements) 6460 { 6461 if(GetBCName(se.index-1) != s1) 6462 continue; 6463 6464 for(const auto& pi : se.PNums()) 6465 { 6466 if(identified_points.find(pi) != identified_points.end()) 6467 continue; 6468 auto pt = (*this)[pi]; 6469 auto mapped_pt = mapping(pt); 6470 auto other_nr = GetElementOfPoint(mapped_pt, lami, true); 6471 int index = -1; 6472 if(other_nr != 0) 6473 { 6474 auto other_el = VolumeElement(other_nr); 6475 for(auto i : Range(other_el.PNums().Size())) 6476 if((mapped_pt - (*this)[other_el.PNums()[i]]).Length() < pointTolerance) 6477 { 6478 index = i; 6479 break; 6480 } 6481 if(index == -1) 6482 { 6483 cout << "point coordinates = " << pt << endl; 6484 cout << "mapped coordinates = " << mapped_pt << endl; 6485 throw Exception("Did not find mapped point with nr " + ToString(pi) + ", are you sure your mesh is periodic?"); 6486 } 6487 auto other_pi = other_el.PNums()[index]; 6488 identified_points.insert(pi); 6489 ident->Add(pi, other_pi, nr); 6490 } 6491 else 6492 { 6493 cout << "point coordinates = " << pt << endl; 6494 cout << "mapped coordinates = " << mapped_pt << endl; 6495 throw Exception("Mapped point with nr " + ToString(pi) + " is outside of mesh, are you sure your mesh is periodic?"); 6496 } 6497 } 6498 } 6499 return nr; 6500 } 6501 InitPointCurve(double red,double green,double blue) const6502 void Mesh :: InitPointCurve(double red, double green, double blue) const 6503 { 6504 pointcurves_startpoint.Append(pointcurves.Size()); 6505 pointcurves_red.Append(red); 6506 pointcurves_green.Append(green); 6507 pointcurves_blue.Append(blue); 6508 } AddPointCurvePoint(const Point3d & pt) const6509 void Mesh :: AddPointCurvePoint(const Point3d & pt) const 6510 { 6511 pointcurves.Append(pt); 6512 } GetNumPointCurves(void) const6513 int Mesh :: GetNumPointCurves(void) const 6514 { 6515 return pointcurves_startpoint.Size(); 6516 } GetNumPointsOfPointCurve(int curve) const6517 int Mesh :: GetNumPointsOfPointCurve(int curve) const 6518 { 6519 if(curve == pointcurves_startpoint.Size()-1) 6520 return (pointcurves.Size() - pointcurves_startpoint.Last()); 6521 else 6522 return (pointcurves_startpoint[curve+1]-pointcurves_startpoint[curve]); 6523 } 6524 GetPointCurvePoint(int curve,int n) const6525 Point3d & Mesh :: GetPointCurvePoint(int curve, int n) const 6526 { 6527 return pointcurves[pointcurves_startpoint[curve]+n]; 6528 } 6529 GetPointCurveColor(int curve,double & red,double & green,double & blue) const6530 void Mesh :: GetPointCurveColor(int curve, double & red, double & green, double & blue) const 6531 { 6532 red = pointcurves_red[curve]; 6533 green = pointcurves_green[curve]; 6534 blue = pointcurves_blue[curve]; 6535 } 6536 6537 ComputeNVertices()6538 void Mesh :: ComputeNVertices () 6539 { 6540 numvertices = 0; 6541 6542 for (const Element & el : VolumeElements()) 6543 for (PointIndex v : el.Vertices()) 6544 if (v > numvertices) numvertices = v; 6545 6546 for (const Element2d & el : SurfaceElements()) 6547 for (PointIndex v : el.Vertices()) 6548 if (v > numvertices) numvertices = v; 6549 6550 numvertices += 1-PointIndex::BASE; 6551 } 6552 GetNV() const6553 int Mesh :: GetNV () const 6554 { 6555 if (numvertices < 0) 6556 return GetNP(); 6557 else 6558 return numvertices; 6559 } 6560 SetNP(int np)6561 void Mesh :: SetNP (int np) 6562 { 6563 points.SetSize(np); 6564 // ptyps.SetSize(np); 6565 6566 int mlold = mlbetweennodes.Size(); 6567 mlbetweennodes.SetSize(np); 6568 if (np > mlold) 6569 for (int i = mlold+PointIndex::BASE; 6570 i < np+PointIndex::BASE; i++) 6571 { 6572 mlbetweennodes[i].I1() = PointIndex::BASE-1; 6573 mlbetweennodes[i].I2() = PointIndex::BASE-1; 6574 } 6575 6576 GetIdentifications().SetMaxPointNr (np + PointIndex::BASE-1); 6577 } 6578 6579 CreatePoint2ElementTable(std::optional<BitArray> points) const6580 Table<ElementIndex, PointIndex> Mesh :: CreatePoint2ElementTable(std::optional<BitArray> points) const 6581 { 6582 if(points) 6583 { 6584 const auto & free_points = *points; 6585 return ngcore::CreateSortedTable<ElementIndex, PointIndex>( volelements.Range(), 6586 [&](auto & table, ElementIndex ei) 6587 { 6588 const auto & el = (*this)[ei]; 6589 if(el.IsDeleted()) 6590 return; 6591 6592 for (PointIndex pi : el.PNums()) 6593 if(free_points[pi]) 6594 table.Add (pi, ei); 6595 }, GetNP()); 6596 } 6597 else 6598 return ngcore::CreateSortedTable<ElementIndex, PointIndex>( volelements.Range(), 6599 [&](auto & table, ElementIndex ei) 6600 { 6601 const auto & el = (*this)[ei]; 6602 if(el.IsDeleted()) 6603 return; 6604 6605 for (PointIndex pi : el.PNums()) 6606 table.Add (pi, ei); 6607 }, GetNP()); 6608 } 6609 CreatePoint2SurfaceElementTable(int faceindex) const6610 Table<SurfaceElementIndex, PointIndex> Mesh :: CreatePoint2SurfaceElementTable( int faceindex ) const 6611 { 6612 static Timer timer("Mesh::CreatePoint2SurfaceElementTable"); RegionTimer rt(timer); 6613 6614 if(faceindex==0) 6615 { 6616 return ngcore::CreateSortedTable<SurfaceElementIndex, PointIndex>( surfelements.Range(), 6617 [&](auto & table, SurfaceElementIndex ei) 6618 { 6619 for (PointIndex pi : (*this)[ei].PNums()) 6620 table.Add (pi, ei); 6621 }, GetNP()); 6622 } 6623 6624 Array<SurfaceElementIndex> face_els; 6625 GetSurfaceElementsOfFace(faceindex, face_els); 6626 return ngcore::CreateSortedTable<SurfaceElementIndex, PointIndex>( face_els.Range(), 6627 [&](auto & table, size_t i) 6628 { 6629 for (PointIndex pi : (*this)[face_els[i]].PNums()) 6630 table.Add (pi, face_els[i]); 6631 }, GetNP()); 6632 } 6633 6634 6635 6636 /* 6637 void Mesh :: BuildConnectedNodes () 6638 { 6639 if (PureTetMesh()) 6640 { 6641 connectedtonode.SetSize(0); 6642 return; 6643 } 6644 6645 6646 int i, j, k; 6647 int np = GetNP(); 6648 int ne = GetNE(); 6649 TABLE<int> conto(np); 6650 for (i = 1; i <= ne; i++) 6651 { 6652 const Element & el = VolumeElement(i); 6653 6654 if (el.GetType() == PRISM) 6655 { 6656 for (j = 1; j <= 6; j++) 6657 { 6658 int n1 = el.PNum (j); 6659 int n2 = el.PNum ((j+2)%6+1); 6660 // if (n1 != n2) 6661 { 6662 int found = 0; 6663 for (k = 1; k <= conto.EntrySize(n1); k++) 6664 if (conto.Get(n1, k) == n2) 6665 { 6666 found = 1; 6667 break; 6668 } 6669 if (!found) 6670 conto.Add (n1, n2); 6671 } 6672 } 6673 } 6674 else if (el.GetType() == PYRAMID) 6675 { 6676 for (j = 1; j <= 4; j++) 6677 { 6678 int n1, n2; 6679 switch (j) 6680 { 6681 case 1: n1 = 1; n2 = 4; break; 6682 case 2: n1 = 4; n2 = 1; break; 6683 case 3: n1 = 2; n2 = 3; break; 6684 case 4: n1 = 3; n2 = 2; break; 6685 } 6686 6687 int found = 0; 6688 for (k = 1; k <= conto.EntrySize(n1); k++) 6689 if (conto.Get(n1, k) == n2) 6690 { 6691 found = 1; 6692 break; 6693 } 6694 if (!found) 6695 conto.Add (n1, n2); 6696 } 6697 } 6698 } 6699 6700 connectedtonode.SetSize(np); 6701 for (i = 1; i <= np; i++) 6702 connectedtonode.Elem(i) = 0; 6703 6704 for (i = 1; i <= np; i++) 6705 if (connectedtonode.Elem(i) == 0) 6706 { 6707 connectedtonode.Elem(i) = i; 6708 ConnectToNodeRec (i, i, conto); 6709 } 6710 6711 6712 6713 } 6714 6715 void Mesh :: ConnectToNodeRec (int node, int tonode, 6716 const TABLE<int> & conto) 6717 { 6718 int i, n2; 6719 // (*testout) << "connect " << node << " to " << tonode << endl; 6720 for (i = 1; i <= conto.EntrySize(node); i++) 6721 { 6722 n2 = conto.Get(node, i); 6723 if (!connectedtonode.Get(n2)) 6724 { 6725 connectedtonode.Elem(n2) = tonode; 6726 ConnectToNodeRec (n2, tonode, conto); 6727 } 6728 } 6729 } 6730 */ 6731 6732 PureTrigMesh(int faceindex) const6733 bool Mesh :: PureTrigMesh (int faceindex) const 6734 { 6735 // if (!faceindex) return !mparam.quad; 6736 6737 if (!faceindex) 6738 { 6739 for (int i = 1; i <= GetNSE(); i++) 6740 if (SurfaceElement(i).GetNP() != 3) 6741 return false; 6742 return true; 6743 } 6744 6745 for (int i = 1; i <= GetNSE(); i++) 6746 if (SurfaceElement(i).GetIndex() == faceindex && 6747 SurfaceElement(i).GetNP() != 3) 6748 return false; 6749 return true; 6750 } 6751 PureTetMesh() const6752 bool Mesh :: PureTetMesh () const 6753 { 6754 for (ElementIndex ei = 0; ei < GetNE(); ei++) 6755 if (VolumeElement(ei).GetNP() != 4) 6756 return 0; 6757 return 1; 6758 } 6759 UpdateTopology(NgTaskManager tm,NgTracer tracer)6760 void Mesh :: UpdateTopology (NgTaskManager tm, 6761 NgTracer tracer) 6762 { 6763 static Timer t("Update Topology"); RegionTimer reg(t); 6764 topology.Update(tm, tracer); 6765 (*tracer)("call update clusters", false); 6766 clusters->Update(); 6767 (*tracer)("call update clusters", true); 6768 #ifdef PARALLEL 6769 if (paralleltop) 6770 { 6771 paralleltop->Reset(); 6772 paralleltop->UpdateCoarseGrid(); 6773 } 6774 #endif 6775 updateSignal.Emit(); 6776 } 6777 BuildCurvedElements(const Refinement * ref,int aorder,bool arational)6778 void Mesh :: BuildCurvedElements (const Refinement * ref, int aorder, bool arational) 6779 { 6780 GetCurvedElements().BuildCurvedElements (ref, aorder, arational); 6781 6782 for (SegmentIndex seg = 0; seg < GetNSeg(); seg++) 6783 (*this)[seg].SetCurved (GetCurvedElements().IsSegmentCurved (seg)); 6784 for (SurfaceElementIndex sei = 0; sei < GetNSE(); sei++) 6785 (*this)[sei].SetCurved (GetCurvedElements().IsSurfaceElementCurved (sei)); 6786 for (ElementIndex ei = 0; ei < GetNE(); ei++) 6787 (*this)[ei].SetCurved (GetCurvedElements().IsElementCurved (ei)); 6788 6789 SetNextMajorTimeStamp(); 6790 } 6791 BuildCurvedElements(int aorder)6792 void Mesh :: BuildCurvedElements (int aorder) 6793 { 6794 if (!GetGeometry()) 6795 throw NgException ("don't have a geometry for mesh curving"); 6796 6797 GetCurvedElements().BuildCurvedElements (&GetGeometry()->GetRefinement(), aorder, false); 6798 6799 for (SegmentIndex seg = 0; seg < GetNSeg(); seg++) 6800 (*this)[seg].SetCurved (GetCurvedElements().IsSegmentCurved (seg)); 6801 for (SurfaceElementIndex sei = 0; sei < GetNSE(); sei++) 6802 (*this)[sei].SetCurved (GetCurvedElements().IsSurfaceElementCurved (sei)); 6803 for (ElementIndex ei = 0; ei < GetNE(); ei++) 6804 (*this)[ei].SetCurved (GetCurvedElements().IsElementCurved (ei)); 6805 6806 SetNextMajorTimeStamp(); 6807 } 6808 SetMaterial(int domnr,const string & mat)6809 void Mesh :: SetMaterial (int domnr, const string & mat) 6810 { 6811 if (domnr > materials.Size()) 6812 { 6813 int olds = materials.Size(); 6814 materials.SetSize (domnr); 6815 for (int i = olds; i < domnr-1; i++) 6816 materials[i] = new string("default"); 6817 } 6818 /* 6819 materials.Elem(domnr) = new char[strlen(mat)+1]; 6820 strcpy (materials.Elem(domnr), mat); 6821 */ 6822 materials.Elem(domnr) = new string(mat); 6823 } 6824 6825 string Mesh :: defaultmat = "default"; GetMaterial(int domnr) const6826 const string & Mesh :: GetMaterial (int domnr) const 6827 { 6828 if (domnr <= materials.Size()) 6829 return *materials.Get(domnr); 6830 static string emptystring("default"); 6831 return emptystring; 6832 } 6833 SetNBCNames(int nbcn)6834 void Mesh ::SetNBCNames ( int nbcn ) 6835 { 6836 if ( bcnames.Size() ) 6837 for ( int i = 0; i < bcnames.Size(); i++) 6838 if ( bcnames[i] ) delete bcnames[i]; 6839 bcnames.SetSize(nbcn); 6840 bcnames = 0; 6841 } 6842 SetBCName(int bcnr,const string & abcname)6843 void Mesh ::SetBCName ( int bcnr, const string & abcname ) 6844 { 6845 if (bcnr >= bcnames.Size()) 6846 { 6847 int oldsize = bcnames.Size(); 6848 bcnames.SetSize (bcnr+1); // keeps contents 6849 for (int i = oldsize; i <= bcnr; i++) 6850 bcnames[i] = nullptr; 6851 } 6852 6853 if ( bcnames[bcnr] ) delete bcnames[bcnr]; 6854 if ( abcname != "default" ) 6855 bcnames[bcnr] = new string ( abcname ); 6856 else 6857 bcnames[bcnr] = nullptr; 6858 6859 for (auto & fd : facedecoding) 6860 if (fd.BCProperty() <= bcnames.Size()) 6861 fd.SetBCName (bcnames[fd.BCProperty()-1]); 6862 } 6863 GetBCName(int bcnr) const6864 const string & Mesh ::GetBCName ( int bcnr ) const 6865 { 6866 static string defaultstring = "default"; 6867 6868 if ( !bcnames.Size() ) 6869 return defaultstring; 6870 6871 if (bcnr < 0 || bcnr >= bcnames.Size()) 6872 throw RangeException("Illegal bc number ", bcnr, 0, bcnames.Size()); 6873 6874 if ( bcnames[bcnr] ) 6875 return *bcnames[bcnr]; 6876 else 6877 return defaultstring; 6878 } 6879 SetNCD2Names(int ncd2n)6880 void Mesh :: SetNCD2Names( int ncd2n ) 6881 { 6882 if (cd2names.Size()) 6883 for(int i=0; i<cd2names.Size(); i++) 6884 if(cd2names[i]) delete cd2names[i]; 6885 cd2names.SetSize(ncd2n); 6886 cd2names = 0; 6887 } 6888 SetCD2Name(int cd2nr,const string & abcname)6889 void Mesh :: SetCD2Name ( int cd2nr, const string & abcname ) 6890 { 6891 cd2nr--; 6892 (*testout) << "setCD2Name on edge " << cd2nr << " to " << abcname << endl; 6893 if (cd2nr >= cd2names.Size()) 6894 { 6895 int oldsize = cd2names.Size(); 6896 cd2names.SetSize(cd2nr+1); 6897 for(int i= oldsize; i<= cd2nr; i++) 6898 cd2names[i] = nullptr; 6899 } 6900 //if (cd2names[cd2nr]) delete cd2names[cd2nr]; 6901 if (abcname != "default") 6902 cd2names[cd2nr] = new string(abcname); 6903 else 6904 cd2names[cd2nr] = nullptr; 6905 } 6906 6907 string Mesh :: cd2_default_name = "default"; 6908 string Mesh :: default_bc = "default"; GetCD2Name(int cd2nr) const6909 const string & Mesh :: GetCD2Name (int cd2nr) const 6910 { 6911 static string defaultstring = "default"; 6912 if (!cd2names.Size()) 6913 return defaultstring; 6914 6915 if (cd2nr < 0 || cd2nr >= cd2names.Size()) 6916 return defaultstring; 6917 6918 if (cd2names[cd2nr]) 6919 return *cd2names[cd2nr]; 6920 else 6921 return defaultstring; 6922 } 6923 SetNCD3Names(int ncd3n)6924 void Mesh :: SetNCD3Names( int ncd3n ) 6925 { 6926 if (cd3names.Size()) 6927 for(int i=0; i<cd3names.Size(); i++) 6928 if(cd3names[i]) delete cd3names[i]; 6929 cd3names.SetSize(ncd3n); 6930 cd3names = 0; 6931 } 6932 SetCD3Name(int cd3nr,const string & abcname)6933 void Mesh :: SetCD3Name ( int cd3nr, const string & abcname ) 6934 { 6935 cd3nr--; 6936 (*testout) << "setCD3Name on vertex " << cd3nr << " to " << abcname << endl; 6937 if (cd3nr >= cd3names.Size()) 6938 { 6939 int oldsize = cd3names.Size(); 6940 cd3names.SetSize(cd3nr+1); 6941 for(int i= oldsize; i<= cd3nr; i++) 6942 cd3names[i] = nullptr; 6943 } 6944 if (abcname != "default") 6945 cd3names[cd3nr] = new string(abcname); 6946 else 6947 cd3names[cd3nr] = nullptr; 6948 } 6949 AddCD3Name(const string & aname)6950 int Mesh :: AddCD3Name (const string & aname) 6951 { 6952 for (int i = 0; i < cd3names.Size(); i++) 6953 if (*cd3names[i] == aname) 6954 return i; 6955 cd3names.Append (new string(aname)); 6956 return cd3names.Size()-1; 6957 } 6958 6959 string Mesh :: cd3_default_name = "default"; GetCD3Name(int cd3nr) const6960 const string & Mesh :: GetCD3Name (int cd3nr) const 6961 { 6962 static string defaultstring = "default"; 6963 if (!cd3names.Size()) 6964 return defaultstring; 6965 6966 if (cd3nr < 0 || cd3nr >= cd3names.Size()) 6967 return defaultstring; 6968 6969 if (cd3names[cd3nr]) 6970 return *cd3names[cd3nr]; 6971 else 6972 return defaultstring; 6973 } 6974 6975 GetRegionNamesCD(int codim)6976 NgArray<string*> & Mesh :: GetRegionNamesCD (int codim) 6977 { 6978 switch (codim) 6979 { 6980 case 0: return materials; 6981 case 1: return bcnames; 6982 case 2: return cd2names; 6983 case 3: return cd3names; 6984 default: throw Exception("don't have regions of co-dimension "+ToString(codim)); 6985 } 6986 } 6987 6988 6989 SetUserData(const char * id,NgArray<int> & data)6990 void Mesh :: SetUserData(const char * id, NgArray<int> & data) 6991 { 6992 if(userdata_int.Used(id)) 6993 delete userdata_int[id]; 6994 6995 NgArray<int> * newdata = new NgArray<int>(data); 6996 6997 userdata_int.Set(id,newdata); 6998 } GetUserData(const char * id,NgArray<int> & data,int shift) const6999 bool Mesh :: GetUserData(const char * id, NgArray<int> & data, int shift) const 7000 { 7001 if(userdata_int.Used(id)) 7002 { 7003 if(data.Size() < (*userdata_int[id]).Size()+shift) 7004 data.SetSize((*userdata_int[id]).Size()+shift); 7005 for(int i=0; i<(*userdata_int[id]).Size(); i++) 7006 data[i+shift] = (*userdata_int[id])[i]; 7007 return true; 7008 } 7009 else 7010 { 7011 data.SetSize(0); 7012 return false; 7013 } 7014 } SetUserData(const char * id,NgArray<double> & data)7015 void Mesh :: SetUserData(const char * id, NgArray<double> & data) 7016 { 7017 if(userdata_double.Used(id)) 7018 delete userdata_double[id]; 7019 7020 NgArray<double> * newdata = new NgArray<double>(data); 7021 7022 userdata_double.Set(id,newdata); 7023 } GetUserData(const char * id,NgArray<double> & data,int shift) const7024 bool Mesh :: GetUserData(const char * id, NgArray<double> & data, int shift) const 7025 { 7026 if(userdata_double.Used(id)) 7027 { 7028 if(data.Size() < (*userdata_double[id]).Size()+shift) 7029 data.SetSize((*userdata_double[id]).Size()+shift); 7030 for(int i=0; i<(*userdata_double[id]).Size(); i++) 7031 data[i+shift] = (*userdata_double[id])[i]; 7032 return true; 7033 } 7034 else 7035 { 7036 data.SetSize(0); 7037 return false; 7038 } 7039 } 7040 7041 7042 PrintMemInfo(ostream & ost) const7043 void Mesh :: PrintMemInfo (ostream & ost) const 7044 { 7045 ost << "Mesh Mem:" << endl; 7046 7047 ost << GetNP() << " Points, of size " 7048 << sizeof (Point3d) << " + " << sizeof(POINTTYPE) << " = " 7049 << GetNP() * (sizeof (Point3d) + sizeof(POINTTYPE)) << endl; 7050 7051 ost << GetNSE() << " Surface elements, of size " 7052 << sizeof (Element2d) << " = " 7053 << GetNSE() * sizeof(Element2d) << endl; 7054 7055 ost << GetNE() << " Volume elements, of size " 7056 << sizeof (Element) << " = " 7057 << GetNE() * sizeof(Element) << endl; 7058 7059 // ost << "surfs on node:"; 7060 // surfacesonnode.PrintMemInfo (cout); 7061 7062 ost << "boundaryedges: "; 7063 if (boundaryedges) 7064 boundaryedges->PrintMemInfo (cout); 7065 7066 ost << "surfelementht: "; 7067 if (surfelementht) 7068 surfelementht->PrintMemInfo (cout); 7069 } 7070 Mirror(netgen::Point<3> p_plane,Vec<3> n_plane)7071 shared_ptr<Mesh> Mesh :: Mirror ( netgen::Point<3> p_plane, Vec<3> n_plane ) 7072 { 7073 Mesh & m = *this; 7074 auto nm_ = make_shared<Mesh>(); 7075 Mesh & nm = *nm_; 7076 nm = m; 7077 7078 Point3d pmin, pmax; 7079 GetBox(pmin, pmax); 7080 auto v = pmax-pmin; 7081 double eps = v.Length()*1e-8; 7082 7083 auto onPlane = [&] (const MeshPoint & p) -> bool 7084 { 7085 auto v = p_plane-p; 7086 auto l = v.Length(); 7087 if(l<eps) return true; 7088 7089 auto ret = fabs(v*n_plane)/l; 7090 return fabs(v*n_plane) < eps; 7091 }; 7092 7093 auto mirror = [&] (PointIndex pi) -> PointIndex 7094 { 7095 auto & p = m[pi]; 7096 7097 auto v = p_plane-p; 7098 auto l = v.Length(); 7099 if(l<eps) 7100 return pi; 7101 7102 if(fabs(v*n_plane)/l < eps) 7103 return pi; 7104 7105 auto new_point = p + 2*(v*n_plane)*n_plane; 7106 return nm.AddPoint( new_point, p.GetLayer(), p.Type() ); 7107 }; 7108 7109 Array<PointIndex, PointIndex> point_map; 7110 point_map.SetSize(GetNP()); 7111 point_map = -1; 7112 7113 for(auto pi : Range(points)) 7114 point_map[pi] = mirror(pi); 7115 7116 for(auto & el : VolumeElements()) 7117 { 7118 auto nel = el; 7119 for(auto i : Range(el.GetNP())) 7120 nel[i] = point_map[el[i]]; 7121 nm.AddVolumeElement(nel); 7122 } 7123 7124 for (auto ei : Range(SurfaceElements())) 7125 { 7126 auto & el = m[ei]; 7127 auto nel = el; 7128 for(auto i : Range(el.GetNP())) 7129 nel[i] = point_map[el[i]]; 7130 7131 if(!(nel==el)) 7132 nm.AddSurfaceElement(nel); 7133 } 7134 7135 for (auto ei : Range(LineSegments())) 7136 { 7137 auto & el = LineSegments()[ei]; 7138 auto nel = el; 7139 bool is_same = true; 7140 7141 for(auto i : Range(el.GetNP())) 7142 { 7143 auto pi = el[i]; 7144 nel[i] = point_map[pi]; 7145 if(point_map[pi]!=pi) 7146 is_same = false; 7147 } 7148 7149 if(!is_same) 7150 nm.AddSegment(nel); 7151 } 7152 7153 return nm_; 7154 } 7155 7156 } 7157