1 #include <mystdlib.h> 2 #include "meshing.hpp" 3 #include "meshing2.hpp" 4 #include "global.hpp" 5 #include "../geom2d/csg2d.hpp" 6 7 namespace netgen 8 { 9 InsertVirtualBoundaryLayer(Mesh & mesh)10 void InsertVirtualBoundaryLayer (Mesh & mesh) 11 { 12 cout << "Insert virt. b.l." << endl; 13 14 int surfid; 15 16 cout << "Boundary Nr:"; 17 cin >> surfid; 18 19 int i; 20 int np = mesh.GetNP(); 21 22 cout << "Old NP: " << mesh.GetNP() << endl; 23 cout << "Trigs: " << mesh.GetNSE() << endl; 24 25 NgBitArray bndnodes(np); 26 NgArray<int> mapto(np); 27 28 bndnodes.Clear(); 29 for (i = 1; i <= mesh.GetNSeg(); i++) 30 { 31 int snr = mesh.LineSegment(i).edgenr; 32 cout << "snr = " << snr << endl; 33 if (snr == surfid) 34 { 35 bndnodes.Set (mesh.LineSegment(i)[0]); 36 bndnodes.Set (mesh.LineSegment(i)[1]); 37 } 38 } 39 for (i = 1; i <= mesh.GetNSeg(); i++) 40 { 41 int snr = mesh.LineSegment(i).edgenr; 42 if (snr != surfid) 43 { 44 bndnodes.Clear (mesh.LineSegment(i)[0]); 45 bndnodes.Clear (mesh.LineSegment(i)[1]); 46 } 47 } 48 49 for (i = 1; i <= np; i++) 50 { 51 if (bndnodes.Test(i)) 52 mapto.Elem(i) = mesh.AddPoint (mesh.Point (i)); 53 else 54 mapto.Elem(i) = 0; 55 } 56 57 for (i = 1; i <= mesh.GetNSE(); i++) 58 { 59 Element2d & el = mesh.SurfaceElement(i); 60 for (int j = 1; j <= el.GetNP(); j++) 61 if (mapto.Get(el.PNum(j))) 62 el.PNum(j) = mapto.Get(el.PNum(j)); 63 } 64 65 66 int nq = 0; 67 for (i = 1; i <= mesh.GetNSeg(); i++) 68 { 69 int snr = mesh.LineSegment(i).edgenr; 70 if (snr == surfid) 71 { 72 int p1 = mesh.LineSegment(i)[0]; 73 int p2 = mesh.LineSegment(i)[1]; 74 int p3 = mapto.Get (p1); 75 if (!p3) p3 = p1; 76 int p4 = mapto.Get (p2); 77 if (!p4) p4 = p2; 78 79 Element2d el(QUAD); 80 el.PNum(1) = p1; 81 el.PNum(2) = p2; 82 el.PNum(3) = p3; 83 el.PNum(4) = p4; 84 el.SetIndex (2); 85 mesh.AddSurfaceElement (el); 86 nq++; 87 } 88 } 89 90 cout << "New NP: " << mesh.GetNP() << endl; 91 cout << "Quads: " << nq << endl; 92 } 93 GenerateBoundaryLayer(Mesh & mesh,const BoundaryLayerParameters & blp)94 void GenerateBoundaryLayer(Mesh& mesh, const BoundaryLayerParameters& blp) 95 { 96 static Timer timer("Create Boundarylayers"); 97 RegionTimer regt(timer); 98 99 int max_edge_nr = -1; 100 for(const auto& seg : mesh.LineSegments()) 101 if(seg.edgenr > max_edge_nr) 102 max_edge_nr = seg.edgenr; 103 104 int new_mat_nr = mesh.GetNDomains() +1; 105 mesh.SetMaterial(new_mat_nr, blp.new_mat); 106 107 auto domains = blp.domains; 108 if(!blp.outside) 109 domains.Invert(); 110 111 mesh.UpdateTopology(); 112 auto& meshtopo = mesh.GetTopology(); 113 114 int np = mesh.GetNP(); 115 int ne = mesh.GetNE(); 116 int nse = mesh.GetNSE(); 117 int nseg = mesh.GetNSeg(); 118 119 Array<Array<PointIndex>, PointIndex> mapto(np); 120 121 Array<Vec<3>, PointIndex> growthvectors(np); 122 growthvectors = 0.; 123 124 Array<double> surfacefacs(mesh.GetNFD()+1); 125 surfacefacs = 0.; 126 127 auto getSurfaceNormal = [&mesh] (const Element2d& el) 128 { 129 auto v0 = mesh[el[0]]; 130 return Cross(mesh[el[1]]-v0, mesh[el[2]]-v0).Normalize(); 131 }; 132 133 // surface index map 134 Array<int> si_map(mesh.GetNFD()+1); 135 si_map = -1; 136 137 int fd_old = mesh.GetNFD(); 138 139 // create new FaceDescriptors 140 for(auto i : Range(1, fd_old+1)) 141 { 142 const auto& fd = mesh.GetFaceDescriptor(i); 143 string name = fd.GetBCName(); 144 if(blp.surfid.Contains(i)) 145 { 146 if(auto isIn = domains.Test(fd.DomainIn()); isIn != domains.Test(fd.DomainOut())) 147 { 148 int new_si = mesh.GetNFD()+1; 149 surfacefacs[i] = isIn ? 1. : -1.; 150 // -1 surf nr is so that curving does not do anything 151 FaceDescriptor new_fd(-1, isIn ? new_mat_nr : fd.DomainIn(), 152 isIn ? fd.DomainOut() : new_mat_nr, -1); 153 new_fd.SetBCProperty(new_si); 154 mesh.AddFaceDescriptor(new_fd); 155 si_map[i] = new_si; 156 mesh.SetBCName(new_si-1, "mapped_" + name); 157 } 158 } 159 } 160 161 // mark points for remapping 162 for(const auto& sel : mesh.SurfaceElements()) 163 { 164 auto n = surfacefacs[sel.GetIndex()] * getSurfaceNormal(sel); 165 if(n.Length2() != 0.) 166 { 167 for(auto pi : sel.PNums()) 168 { 169 auto & np = growthvectors[pi]; 170 if(np.Length() == 0) { np = n; continue; } 171 auto npn = np * n; 172 auto npnp = np * np; 173 auto nn = n * n; 174 if(nn-npn*npn/npnp == 0) { np = n; continue; } 175 np += (nn - npn)/(nn - npn*npn/npnp) * (n - npn/npnp * np); 176 } 177 } 178 } 179 180 // Bit array to keep track of segments already processed 181 BitArray segs_done(nseg); 182 segs_done.Clear(); 183 184 // map for all segments with same points 185 // points to pair of SegmentIndex, int 186 // int is type of other segment, either: 187 // 0 == adjacent surface grows layer 188 // 1 == adjacent surface doesn't grow layer, but layer ends on it 189 // 2 == adjacent surface is interior surface that ends on layer 190 // 3 == adjacent surface is exterior surface that ends on layer (not allowed yet) 191 Array<Array<pair<SegmentIndex, int>>, SegmentIndex> segmap(mesh.GetNSeg()); 192 193 // moved segments 194 Array<SegmentIndex> moved_segs; 195 196 // boundaries to project endings to 197 BitArray project_boundaries(fd_old+1); 198 BitArray move_boundaries(fd_old+1); 199 project_boundaries.Clear(); 200 move_boundaries.Clear(); 201 202 Array<SurfaceElementIndex, SegmentIndex> seg2surfel(mesh.GetNSeg()); 203 for(auto si : Range(mesh.SurfaceElements())) 204 { 205 NgArray<int> surfeledges; 206 meshtopo.GetSurfaceElementEdges(si+1, surfeledges); 207 for(auto edgenr : surfeledges) 208 for(auto sei : Range(mesh.LineSegments())) 209 if(meshtopo.GetEdge(sei)+1 == edgenr && 210 mesh[sei].si == mesh[si].GetIndex()) 211 seg2surfel[sei] = si; 212 } 213 214 for(auto si : Range(mesh.LineSegments())) 215 { 216 if(segs_done[si]) continue; 217 const auto& segi = mesh[si]; 218 if(si_map[segi.si] == -1) continue; 219 segs_done.SetBit(si); 220 segmap[si].Append(make_pair(si, 0)); 221 moved_segs.Append(si); 222 for(auto sj : Range(mesh.LineSegments())) 223 { 224 if(segs_done.Test(sj)) continue; 225 const auto& segj = mesh[sj]; 226 if((segi[0] == segj[0] && segi[1] == segj[1]) || 227 (segi[0] == segj[1] && segi[1] == segj[0])) 228 { 229 segs_done.SetBit(sj); 230 int type; 231 if(si_map[segj.si] != -1) 232 type = 0; 233 else if(const auto& fd = mesh.GetFaceDescriptor(segj.si); domains.Test(fd.DomainIn()) && domains.Test(fd.DomainOut())) 234 { 235 type = 2; 236 if(fd.DomainIn() == 0 || fd.DomainOut() == 0) 237 project_boundaries.SetBit(segj.si); 238 } 239 else if(const auto& fd = mesh.GetFaceDescriptor(segj.si); !domains.Test(fd.DomainIn()) && !domains.Test(fd.DomainOut())) 240 { 241 type = 3; 242 if(fd.DomainIn() == 0 || fd.DomainOut() == 0) 243 project_boundaries.SetBit(segj.si); 244 move_boundaries.SetBit(segj.si); 245 } 246 else 247 { 248 type = 1; 249 // in case 1 we project the growthvector onto the surface 250 project_boundaries.SetBit(segj.si); 251 } 252 segmap[si].Append(make_pair(sj, type)); 253 } 254 } 255 } 256 257 BitArray in_surface_direction(fd_old+1); 258 in_surface_direction.Clear(); 259 260 // project growthvector on surface for inner angles 261 if(blp.grow_edges) 262 { 263 for(const auto& sel : mesh.SurfaceElements()) 264 if(project_boundaries.Test(sel.GetIndex())) 265 { 266 auto n = getSurfaceNormal(sel); 267 for(auto i : Range(sel.PNums())) 268 { 269 auto pi = sel.PNums()[i]; 270 if(growthvectors[pi].Length2() == 0.) 271 continue; 272 auto next = sel.PNums()[(i+1)%sel.GetNV()]; 273 auto prev = sel.PNums()[i == 0 ? sel.GetNV()-1 : i-1]; 274 auto v1 = (mesh[next] - mesh[pi]).Normalize(); 275 auto v2 = (mesh[prev] - mesh[pi]).Normalize(); 276 auto v3 = growthvectors[pi]; 277 v3.Normalize(); 278 if((v1 * v3 > 1e-12) || (v2 * v3 > 1e-12)) 279 in_surface_direction.SetBit(sel.GetIndex()); 280 281 auto& g = growthvectors[pi]; 282 auto ng = n * g; 283 auto gg = g * g; 284 auto nn = n * n; 285 // if(fabs(ng*ng-nn*gg) < 1e-12 || fabs(ng) < 1e-12) continue; 286 auto a = -ng*ng/(ng*ng-nn * gg); 287 auto b = ng*gg/(ng*ng-nn*gg); 288 g += a*g + b*n; 289 } 290 } 291 } 292 else 293 { 294 for(const auto& seg : mesh.LineSegments()) 295 { 296 int count = 0; 297 for(const auto& seg2 : mesh.LineSegments()) 298 if(((seg[0] == seg2[0] && seg[1] == seg2[1]) || (seg[0] == seg2[1] && seg[1] == seg2[0])) && blp.surfid.Contains(seg2.si)) 299 count++; 300 if(count == 1) 301 { 302 growthvectors[seg[0]] = {0., 0., 0.}; 303 growthvectors[seg[1]] = {0., 0., 0.}; 304 } 305 } 306 } 307 308 // insert new points 309 for (PointIndex pi = 1; pi <= np; pi++) 310 if (growthvectors[pi].Length2() != 0) 311 { 312 Point<3> p = mesh[pi]; 313 for(auto i : Range(blp.heights)) 314 { 315 p += blp.heights[i] * growthvectors[pi]; 316 mapto[pi].Append(mesh.AddPoint(p)); 317 } 318 } 319 320 // add 2d quads on required surfaces 321 map<pair<PointIndex, PointIndex>, int> seg2edge; 322 if(blp.grow_edges) 323 { 324 for(auto sei : moved_segs) 325 { 326 // copy here since we will add segments and this would 327 // invalidate a reference! 328 auto segi = mesh[sei]; 329 for(auto [sej, type] : segmap[sei]) 330 { 331 auto segj = mesh[sej]; 332 if(type == 0) 333 { 334 Segment s; 335 s[0] = mapto[segj[0]].Last(); 336 s[1] = mapto[segj[1]].Last(); 337 s[2] = PointIndex::INVALID; 338 auto pair = s[0] < s[1] ? make_pair(s[0], s[1]) : make_pair(s[1], s[0]); 339 if(seg2edge.find(pair) == seg2edge.end()) 340 seg2edge[pair] = ++max_edge_nr; 341 s.edgenr = seg2edge[pair]; 342 s.si = si_map[segj.si]; 343 mesh.AddSegment(s); 344 } 345 // here we need to grow the quad elements 346 else if(type == 1) 347 { 348 PointIndex pp1 = segj[1]; 349 PointIndex pp2 = segj[0]; 350 if(in_surface_direction.Test(segj.si)) 351 { 352 Swap(pp1, pp2); 353 move_boundaries.SetBit(segj.si); 354 } 355 PointIndex p1 = pp1; 356 PointIndex p2 = pp2; 357 PointIndex p3, p4; 358 Segment s0; 359 s0[0] = p1; 360 s0[1] = p2; 361 s0[2] = PointIndex::INVALID; 362 s0.edgenr = segj.edgenr; 363 s0.si = segj.si; 364 mesh.AddSegment(s0); 365 366 for(auto i : Range(blp.heights)) 367 { 368 Element2d sel(QUAD); 369 p3 = mapto[pp2][i]; 370 p4 = mapto[pp1][i]; 371 sel[0] = p1; 372 sel[1] = p2; 373 sel[2] = p3; 374 sel[3] = p4; 375 sel.SetIndex(segj.si); 376 mesh.AddSurfaceElement(sel); 377 378 // TODO: Too many, would be enough to only add outermost ones 379 Segment s1; 380 s1[0] = p2; 381 s1[1] = p3; 382 s1[2] = PointIndex::INVALID; 383 auto pair = make_pair(p2, p3); 384 if(seg2edge.find(pair) == seg2edge.end()) 385 seg2edge[pair] = ++max_edge_nr; 386 s1.edgenr = seg2edge[pair]; 387 s1.si = segj.si; 388 mesh.AddSegment(s1); 389 Segment s2; 390 s2[0] = p4; 391 s2[1] = p1; 392 s2[2] = PointIndex::INVALID; 393 pair = make_pair(p1, p4); 394 if(seg2edge.find(pair) == seg2edge.end()) 395 seg2edge[pair] = ++max_edge_nr; 396 s2.edgenr = seg2edge[pair]; 397 s2.si = segj.si; 398 mesh.AddSegment(s2); 399 p1 = p4; 400 p2 = p3; 401 } 402 Segment s3; 403 s3[0] = p3; 404 s3[1] = p4; 405 s3[2] = PointIndex::INVALID; 406 auto pair = p3 < p4 ? make_pair(p3, p4) : make_pair(p4, p3); 407 if(seg2edge.find(pair) == seg2edge.end()) 408 seg2edge[pair] = ++max_edge_nr; 409 s3.edgenr = seg2edge[pair]; 410 s3.si = segj.si; 411 mesh.AddSegment(s3); 412 } 413 } 414 } 415 } 416 417 BitArray fixed_points(np+1); 418 fixed_points.Clear(); 419 BitArray moveboundarypoint(np+1); 420 moveboundarypoint.Clear(); 421 for(SurfaceElementIndex si = 0; si < nse; si++) 422 { 423 // copy because surfaceels array will be resized! 424 auto sel = mesh[si]; 425 if(si_map[sel.GetIndex()] != -1) 426 { 427 Array<PointIndex> points(sel.PNums()); 428 if(surfacefacs[sel.GetIndex()] > 0) Swap(points[0], points[2]); 429 for(auto j : Range(blp.heights)) 430 { 431 auto eltype = points.Size() == 3 ? PRISM : HEX; 432 Element el(eltype); 433 for(auto i : Range(points)) 434 el[i] = points[i]; 435 for(auto i : Range(points)) 436 points[i] = mapto[sel.PNums()[i]][j]; 437 if(surfacefacs[sel.GetIndex()] > 0) Swap(points[0], points[2]); 438 for(auto i : Range(points)) 439 el[sel.PNums().Size() + i] = points[i]; 440 el.SetIndex(new_mat_nr); 441 mesh.AddVolumeElement(el); 442 } 443 Element2d newel = sel; 444 for(auto& p : newel.PNums()) 445 p = mapto[p].Last(); 446 newel.SetIndex(si_map[sel.GetIndex()]); 447 mesh.AddSurfaceElement(newel); 448 } 449 else 450 { 451 bool has_moved = false; 452 for(auto p : sel.PNums()) 453 if(mapto[p].Size()) 454 has_moved = true; 455 if(has_moved) 456 for(auto p : sel.PNums()) 457 { 458 if(!mapto[p].Size()) 459 { 460 fixed_points.SetBit(p); 461 if(move_boundaries.Test(sel.GetIndex())) 462 moveboundarypoint.SetBit(p); 463 } 464 } 465 } 466 if(move_boundaries.Test(sel.GetIndex())) 467 { 468 for(auto& p : mesh[si].PNums()) 469 if(mapto[p].Size()) 470 p = mapto[p].Last(); 471 } 472 } 473 474 for(SegmentIndex sei = 0; sei < nseg; sei++) 475 { 476 auto& seg = mesh[sei]; 477 if(move_boundaries.Test(seg.si)) 478 for(auto& p : seg.PNums()) 479 if(mapto[p].Size()) 480 p = mapto[p].Last(); 481 } 482 483 for(ElementIndex ei = 0; ei < ne; ei++) 484 { 485 auto el = mesh[ei]; 486 ArrayMem<PointIndex,4> fixed; 487 ArrayMem<PointIndex,4> moved; 488 bool moved_bnd = false; 489 for(const auto& p : el.PNums()) 490 { 491 if(fixed_points.Test(p)) 492 fixed.Append(p); 493 if(mapto[p].Size()) 494 moved.Append(p); 495 if(moveboundarypoint.Test(p)) 496 moved_bnd = true; 497 } 498 499 bool do_move, do_insert; 500 if(domains.Test(el.GetIndex())) 501 { 502 do_move = fixed.Size() && moved_bnd; 503 do_insert = do_move; 504 } 505 else 506 { 507 do_move = !fixed.Size() || moved_bnd; 508 do_insert = !do_move; 509 } 510 511 if(do_move) 512 { 513 for(auto& p : mesh[ei].PNums()) 514 if(mapto[p].Size()) 515 p = mapto[p].Last(); 516 } 517 if(do_insert) 518 { 519 if(el.GetType() == TET) 520 { 521 if(moved.Size() == 3) // inner corner 522 { 523 PointIndex p1 = moved[0]; 524 PointIndex p2 = moved[1]; 525 PointIndex p3 = moved[2]; 526 auto v1 = mesh[p1]; 527 auto n = Cross(mesh[p2]-v1, mesh[p3]-v1); 528 auto d = mesh[mapto[p1][0]] - v1; 529 if(n*d > 0) 530 Swap(p2,p3); 531 PointIndex p4 = p1; 532 PointIndex p5 = p2; 533 PointIndex p6 = p3; 534 for(auto i : Range(blp.heights)) 535 { 536 Element nel(PRISM); 537 nel[0] = p4; nel[1] = p5; nel[2] = p6; 538 p4 = mapto[p1][i]; p5 = mapto[p2][i]; p6 = mapto[p3][i]; 539 nel[3] = p4; nel[4] = p5; nel[5] = p6; 540 nel.SetIndex(el.GetIndex()); 541 mesh.AddVolumeElement(nel); 542 } 543 } 544 if(moved.Size() == 2) 545 { 546 if(fixed.Size() == 1) 547 { 548 PointIndex p1 = moved[0]; 549 PointIndex p2 = moved[1]; 550 for(auto i : Range(blp.heights)) 551 { 552 PointIndex p3 = mapto[moved[1]][i]; 553 PointIndex p4 = mapto[moved[0]][i]; 554 Element nel(PYRAMID); 555 nel[0] = p1; 556 nel[1] = p2; 557 nel[2] = p3; 558 nel[3] = p4; 559 nel[4] = el[0] + el[1] + el[2] + el[3] - fixed[0] - moved[0] - moved[1]; 560 if(Cross(mesh[p2]-mesh[p1], mesh[p4]-mesh[p1]) * (mesh[nel[4]]-mesh[nel[1]]) > 0) 561 Swap(nel[1], nel[3]); 562 nel.SetIndex(el.GetIndex()); 563 mesh.AddVolumeElement(nel); 564 p1 = p4; 565 p2 = p3; 566 } 567 } 568 } 569 if(moved.Size() == 1 && fixed.Size() == 1) 570 { 571 PointIndex p1 = moved[0]; 572 for(auto i : Range(blp.heights)) 573 { 574 Element nel = el; 575 PointIndex p2 = mapto[moved[0]][i]; 576 for(auto& p : nel.PNums()) 577 { 578 if(p == moved[0]) 579 p = p1; 580 else if(p == fixed[0]) 581 p = p2; 582 } 583 p1 = p2; 584 mesh.AddVolumeElement(nel); 585 } 586 } 587 } 588 else if(el.GetType() == PYRAMID) 589 { 590 if(moved.Size() == 2) 591 { 592 if(fixed.Size() != 2) 593 throw Exception("This case is not implemented yet! Fixed size = " + ToString(fixed.Size())); 594 PointIndex p1 = moved[0]; 595 PointIndex p2 = moved[1]; 596 for(auto i : Range(blp.heights)) 597 { 598 PointIndex p3 = mapto[moved[1]][i]; 599 PointIndex p4 = mapto[moved[0]][i]; 600 Element nel(PYRAMID); 601 nel[0] = p1; 602 nel[1] = p2; 603 nel[2] = p3; 604 nel[3] = p4; 605 nel[4] = el[0] + el[1] + el[2] + el[3] + el[4] - fixed[0] - fixed[1] - moved[0] - moved[1]; 606 if(Cross(mesh[p2] - mesh[p1], mesh[p4]-mesh[p1]) * (mesh[nel[4]]-mesh[nel[1]]) > 0) 607 Swap(nel[1], nel[3]); 608 nel.SetIndex(el.GetIndex()); 609 mesh.AddVolumeElement(nel); 610 p1 = p4; 611 p2 = p3; 612 } 613 } 614 else if(moved.Size() == 1) 615 throw Exception("This case is not implemented yet!"); 616 } 617 else 618 throw Exception("Boundarylayer only implemented for tets and pyramids outside yet!"); 619 } 620 } 621 622 for(auto i : Range(1, fd_old+1)) 623 if(si_map[i] != -1) 624 { 625 if(mesh.GetFaceDescriptor(mesh.GetNFD()).DomainIn() == new_mat_nr) 626 mesh.GetFaceDescriptor(i).SetDomainOut(new_mat_nr); 627 else 628 mesh.GetFaceDescriptor(i).SetDomainIn(new_mat_nr); 629 } 630 mesh.GetTopology().ClearEdges(); 631 mesh.UpdateTopology(); 632 } 633 AddDirection(Vec<3> & a,Vec<3> b)634 void AddDirection( Vec<3> & a, Vec<3> b ) 635 { 636 if(a.Length2()==0.) 637 { 638 a = b; 639 return; 640 } 641 642 if(b.Length2()==0.) 643 return; 644 645 auto ab = a * b; 646 if(fabs(ab)>1-1e-8) 647 return; 648 649 Mat<2> m; 650 m(0,0) = a[0]; 651 m(0,1) = a[1]; 652 m(1,0) = b[0]; 653 m(1,1) = b[1]; 654 Vec<2> lam; 655 Vec<2> rhs; 656 rhs[0] = a[0]-b[0]; 657 rhs[1] = a[1]-b[1]; 658 659 const auto Dot = [](Vec<3> a, Vec<3> b) 660 { return a[0]*b[0] + a[1]*b[1] + a[2]*b[2]; }; 661 662 rhs[0] = Dot(a,a); 663 rhs[1] = Dot(b,b); 664 665 m.Solve(rhs, lam); 666 a[0] = lam[0]; 667 a[1] = lam[1]; 668 a[2] = 0.0; 669 return; 670 } 671 Generate2dMesh(Mesh & mesh,int domain)672 static void Generate2dMesh( Mesh & mesh, int domain ) 673 { 674 Box<3> box{Box<3>::EMPTY_BOX}; 675 for(const auto & seg : mesh.LineSegments()) 676 if (seg.si == domain) 677 for (auto pi : {seg[0], seg[1]}) 678 box.Add(mesh[pi]); 679 680 MeshingParameters mp; 681 Meshing2 meshing (*mesh.GetGeometry(), mp, box); 682 683 Array<PointIndex, PointIndex> compress(mesh.GetNP()); 684 compress = PointIndex::INVALID; 685 686 PointIndex cnt = PointIndex::BASE; 687 for(const auto & seg : mesh.LineSegments()) 688 if (seg.si == domain) 689 for (auto pi : {seg[0], seg[1]}) 690 if (compress[pi]==PointIndex{PointIndex::INVALID}) 691 { 692 meshing.AddPoint(mesh[pi], pi); 693 compress[pi] = cnt++; 694 } 695 696 PointGeomInfo gi; 697 gi.trignum = domain; 698 for(const auto & seg : mesh.LineSegments()) 699 if (seg.si == domain) 700 meshing.AddBoundaryElement (compress[seg[0]], compress[seg[1]], gi, gi); 701 702 auto oldnf = mesh.GetNSE(); 703 auto res = meshing.GenerateMesh (mesh, mp, mp.maxh, domain); 704 for (SurfaceElementIndex sei : Range(oldnf, mesh.GetNSE())) 705 mesh[sei].SetIndex (domain); 706 707 int hsteps = mp.optsteps2d; 708 709 const char * optstr = mp.optimize2d.c_str(); 710 MeshOptimize2d meshopt(mesh); 711 meshopt.SetFaceIndex(domain); 712 meshopt.SetMetricWeight (mp.elsizeweight); 713 for (size_t j = 1; j <= strlen(optstr); j++) 714 { 715 switch (optstr[j-1]) 716 { 717 case 's': 718 { // topological swap 719 meshopt.EdgeSwapping (0); 720 break; 721 } 722 case 'S': 723 { // metric swap 724 meshopt.EdgeSwapping (1); 725 break; 726 } 727 case 'm': 728 { 729 meshopt.ImproveMesh(mp); 730 break; 731 } 732 case 'c': 733 { 734 meshopt.CombineImprove(); 735 break; 736 } 737 default: 738 cerr << "Optimization code " << optstr[j-1] << " not defined" << endl; 739 } 740 } 741 742 mesh.Compress(); 743 mesh.OrderElements(); 744 mesh.SetNextMajorTimeStamp(); 745 746 } 747 GenerateBoundaryLayer2(Mesh & mesh,int domain,const Array<double> & thicknesses,bool should_make_new_domain,const Array<int> & boundaries)748 int GenerateBoundaryLayer2 (Mesh & mesh, int domain, const Array<double> & thicknesses, bool should_make_new_domain, const Array<int> & boundaries) 749 { 750 SegmentIndex first_new_seg = mesh.LineSegments().Range().Next(); 751 752 int np = mesh.GetNP(); 753 int nseg = mesh.GetNSeg(); 754 int ne = mesh.GetNSE(); 755 mesh.UpdateTopology(); 756 757 double total_thickness = 0.0; 758 for(auto thickness : thicknesses) 759 total_thickness += thickness; 760 761 Array<Array<PointIndex>, PointIndex> mapto(np); 762 763 // Bit array to keep track of segments already processed 764 BitArray segs_done(nseg); 765 segs_done.Clear(); 766 767 // moved segments 768 Array<SegmentIndex> moved_segs; 769 770 Array<Vec<3>, PointIndex> growthvectors(np); 771 growthvectors = 0.; 772 773 auto & meshtopo = mesh.GetTopology(); 774 775 Array<SurfaceElementIndex, SegmentIndex> seg2surfel(mesh.GetNSeg()); 776 seg2surfel = -1; 777 NgArray<SurfaceElementIndex> temp_els; 778 for(auto si : Range(mesh.LineSegments())) 779 { 780 meshtopo.GetSegmentSurfaceElements ( si+1, temp_els ); 781 // NgArray<int> surfeledges; 782 // meshtopo.GetSurfaceElementEdges(si+1, surfeledges); 783 for(auto seli : temp_els) 784 if(mesh[seli].GetIndex() == mesh[si].si) 785 seg2surfel[si] = seli; 786 } 787 788 Array<SegmentIndex> segments; 789 790 // surface index map 791 Array<int> si_map(mesh.GetNFD()+1); 792 si_map = -1; 793 794 int fd_old = mesh.GetNFD(); 795 796 int max_edge_nr = -1; 797 int max_domain = -1; 798 799 for(const auto& seg : mesh.LineSegments()) 800 { 801 if(seg.epgeominfo[0].edgenr > max_edge_nr) 802 max_edge_nr = seg.epgeominfo[0].edgenr; 803 if(seg.si > max_domain) 804 max_domain = seg.si; 805 } 806 807 int new_domain = max_domain+1; 808 809 BitArray active_boundaries(max_edge_nr+1); 810 BitArray active_segments(nseg); 811 active_boundaries.Clear(); 812 active_segments.Clear(); 813 814 if(boundaries.Size() == 0) 815 active_boundaries.Set(); 816 else 817 for(auto edgenr : boundaries) 818 active_boundaries.SetBit(edgenr); 819 820 for(auto segi : Range(mesh.LineSegments())) 821 { 822 const auto seg = mesh[segi]; 823 if(active_boundaries.Test(seg.epgeominfo[0].edgenr) && seg.si==domain) 824 active_segments.SetBit(segi); 825 } 826 827 for(auto segi : Range(mesh.LineSegments())) 828 { 829 const auto& seg = mesh[segi]; 830 auto si = seg.si; 831 832 if(si_map[si]!=-1) 833 continue; 834 835 if(!active_segments.Test(segi)) 836 continue; 837 838 FaceDescriptor new_fd(0, 0, 0, -1); 839 new_fd.SetBCProperty(new_domain); 840 int new_fd_index = mesh.AddFaceDescriptor(new_fd); 841 si_map[si] = new_domain; 842 if(should_make_new_domain) 843 mesh.SetBCName(new_domain-1, "mapped_" + mesh.GetBCName(si-1)); 844 } 845 846 for(auto si : Range(mesh.LineSegments())) 847 { 848 if(segs_done[si]) continue; 849 segs_done.SetBit(si); 850 const auto& segi = mesh[si]; 851 if(si_map[segi.si] == -1) continue; 852 if(!active_boundaries.Test(segi.epgeominfo[0].edgenr)) 853 continue; 854 moved_segs.Append(si); 855 } 856 857 // calculate growth vectors (average normal vectors of adjacent segments at each point) 858 for (auto si : moved_segs) 859 { 860 auto & seg = mesh[si]; 861 862 meshtopo.GetSegmentSurfaceElements ( si+1, temp_els ); 863 ArrayMem<int, 10> seg_domains; 864 865 temp_els.SetSize(0); 866 if(seg2surfel[si]!=-1) 867 temp_els.Append(seg2surfel[si]); 868 869 int n_temp_els = temp_els.Size(); 870 if(n_temp_els==0) 871 continue; 872 873 int dom0 = mesh[temp_els[0]].GetIndex(); 874 int dom1 = n_temp_els==2 ? mesh[temp_els[1]].GetIndex() : 0; 875 876 bool in_dom0 = dom0 == domain; 877 bool in_dom1 = dom1 == domain; 878 879 if(!in_dom0 && !in_dom1) 880 continue; 881 882 int side = in_dom0 ? 0 : 1; 883 884 auto & sel = mesh[ temp_els[side] ]; 885 886 int domain = sel.GetIndex(); 887 Vec<3> pcenter = 0.0; 888 for(auto i : IntRange(sel.GetNP())) 889 { 890 for(auto d : IntRange(3)) 891 pcenter[d] += mesh[sel[i]][d]; 892 } 893 pcenter = 1.0/sel.GetNP() * pcenter; 894 895 auto n = mesh[seg[1]] - mesh[seg[0]]; 896 n = {-n[1], n[0], 0}; 897 n.Normalize(); 898 899 Vec<3> p0{mesh[seg[0]]}; 900 Vec<3> p1{mesh[seg[0]]}; 901 902 903 auto v = pcenter -0.5*(p0+p1); 904 905 const auto Dot = [](Vec<3> a, Vec<3> b) 906 { return a[0]*b[0] + a[1]*b[1] + a[2]*b[2]; }; 907 if(Dot(n, v)<0) 908 n = -1*n; 909 910 AddDirection(growthvectors[seg[0]], n); 911 AddDirection(growthvectors[seg[1]], n); 912 } 913 914 ////////////////////////////////////////////////////////////////////////// 915 // average growthvectors along straight lines to avoid overlaps in corners 916 BitArray points_done(np+1); 917 points_done.Clear(); 918 919 for(auto si : moved_segs) 920 { 921 auto current_seg = mesh[si]; 922 auto current_si = si; 923 924 auto first = current_seg[0]; 925 auto current = -1; 926 auto next = current_seg[1]; 927 928 if(points_done.Test(first)) 929 continue; 930 931 Array<PointIndex> chain; 932 chain.Append(first); 933 934 // first find closed loops of segments 935 while(next != current && next != first) 936 { 937 current = next; 938 points_done.SetBit(current); 939 chain.Append(current); 940 for(auto sj : meshtopo.GetVertexSegments( current )) 941 { 942 if(!active_segments.Test(sj)) 943 continue; 944 945 if(sj!=current_si) 946 { 947 current_si = sj; 948 current_seg = mesh[sj]; 949 950 next = current_seg[0] + current_seg[1] - current; 951 break; 952 } 953 } 954 } 955 956 auto ifirst = 0; 957 auto n = chain.Size(); 958 959 // angle of adjacent segments at points a[i-1], a[i], a[i+1] 960 auto getAngle = [&mesh, &growthvectors] (FlatArray<PointIndex> a, size_t i) 961 { 962 auto n = a.Size(); 963 auto v0 = growthvectors[a[(i+n-1)%n]]; 964 auto v1 = growthvectors[a[i]]; 965 auto v2 = growthvectors[a[(i+1)%n]]; 966 967 auto p0 = mesh[a[(i+n-1)%n]]; 968 auto p1 = mesh[a[i]]; 969 auto p2 = mesh[a[(i+1)%n]]; 970 971 v0 = p1-p0; 972 v1 = p2-p1; 973 974 auto angle = abs(atan2(v1[0], v1[1]) - atan2(v0[0], v0[1])); 975 if(angle>M_PI) 976 angle = 2*M_PI-angle; 977 978 return angle; 979 }; 980 981 // find first corner point 982 while(getAngle(chain, ifirst) < 1e-5 ) 983 ifirst = (ifirst+1)%n; 984 985 // Copy points of closed loop in correct order, starting with a corner 986 Array<PointIndex> pis(n+1); 987 pis.Range(0, n-ifirst) = chain.Range(ifirst, n); 988 pis.Range(n-ifirst, n) = chain.Range(0, n-ifirst); 989 pis[n] = pis[0]; 990 991 Array<double> lengths(n); 992 993 for(auto i : Range(n)) 994 lengths[i] = (mesh[pis[(i+1)%n]] - mesh[pis[i]]).Length(); 995 996 auto averageGrowthVectors = [&] (size_t first, size_t last) 997 { 998 if(first+1 >= last) 999 return; 1000 1001 double total_len = 0.0; 1002 for(auto l : lengths.Range(first, last)) 1003 total_len += l; 1004 1005 double len = lengths[first]; 1006 auto v0 = growthvectors[pis[first]]; 1007 auto v1 = growthvectors[pis[last]]; 1008 1009 for(auto i : Range(first+1, last)) 1010 { 1011 auto pi = pis[i]; 1012 growthvectors[pi] = (len/total_len)*v1 + (1.0-len/total_len)*v0; 1013 len += lengths[i]; 1014 } 1015 }; 1016 1017 auto icurrent = 0; 1018 1019 while(icurrent<n) 1020 { 1021 auto ilast = icurrent+1; 1022 1023 while(getAngle(pis, ilast) < 1e-5 && ilast < n) 1024 ilast++; 1025 1026 // found straight line -> average growth vectors between end points 1027 if(icurrent!=ilast) 1028 averageGrowthVectors(icurrent, ilast); 1029 1030 icurrent = ilast; 1031 } 1032 } 1033 1034 ////////////////////////////////////////////////////////////////////// 1035 // reduce growthvectors where necessary to avoid overlaps/slim regions 1036 const auto getSegmentBox = [&] (SegmentIndex segi) 1037 { 1038 PointIndex pi0=mesh[segi][0], pi1=mesh[segi][1]; 1039 Box<3> box( mesh[pi0], mesh[pi1] ); 1040 box.Add( mesh[pi0]+growthvectors[pi0] ); 1041 box.Add( mesh[pi1]+growthvectors[pi1] ); 1042 return box; 1043 }; 1044 1045 Array<double, PointIndex> growth(np); 1046 growth = 1.0; 1047 1048 const auto Dot = [](auto a, auto b) 1049 { return a[0]*b[0] + a[1]*b[1] + a[2]*b[2]; }; 1050 1051 const auto restrictGrowthVectors = [&] (SegmentIndex segi0, SegmentIndex segi1) 1052 { 1053 if(!active_segments.Test(segi0)) 1054 return; 1055 1056 const auto & seg0 = mesh[segi0]; 1057 const auto & seg1 = mesh[segi1]; 1058 1059 if(seg0.si != seg1.si) 1060 return; 1061 1062 if(segi0 == segi1) 1063 return; 1064 1065 if(seg0[0]==seg1[0] || seg0[0]==seg1[1] || seg0[1]==seg1[0] || seg0[1] == seg1[1]) 1066 return; 1067 1068 auto n = mesh[seg0[0]] - mesh[seg0[1]]; 1069 n = {-n[1], n[0], 0}; 1070 n.Normalize(); 1071 if(Dot(n, growthvectors[seg0[0]])<0) n = -n; 1072 if(Dot(n, growthvectors[seg0[1]])<0) n = -n; 1073 1074 auto n1 = mesh[seg1[0]] - mesh[seg1[1]]; 1075 n1 = {-n1[1], n1[0], 0}; 1076 n1.Normalize(); 1077 if(Dot(n1, growthvectors[seg1[0]])<0) n1 = -n; 1078 if(Dot(n1, growthvectors[seg1[1]])<0) n1 = -n; 1079 1080 auto p10 = mesh[seg1[0]]; 1081 auto p11 = mesh[seg1[1]]; 1082 1083 for ( auto pi : {seg0[0], seg0[1]} ) 1084 { 1085 if(growthvectors[pi] == 0.0) 1086 continue; 1087 1088 PointIndex pi1 = seg0[0] + seg0[1] - pi; 1089 auto p1 = mesh[pi1]; 1090 auto p = mesh[pi]; 1091 1092 Point<3> points[] = { p10, p11, p10+total_thickness*growthvectors[seg1[0]], p11+total_thickness*growthvectors[seg1[1]], p1+total_thickness*growthvectors[pi1] }; 1093 1094 Vec<3> gn{ growthvectors[pi][1], -growthvectors[pi][0], 0.0 }; 1095 if(Dot(gn, p1-p) < 0) 1096 gn = -gn; 1097 1098 double d0 = Dot(gn, p); 1099 double d1 = Dot(gn, p1); 1100 if(d0>d1) 1101 Swap(d0,d1); 1102 1103 bool all_left=true, all_right=true; 1104 1105 for (auto i: Range(4)) 1106 { 1107 auto p_other = points[i]; 1108 auto dot = Dot(gn,p_other); 1109 if(dot>d0) all_left = false; 1110 if(dot<d1) all_right = false; 1111 } 1112 1113 if(all_left || all_right) 1114 return; 1115 1116 //for ( auto pi : {seg0[0], seg0[1]} ) 1117 { 1118 double safety = 1.3; 1119 double t = safety*total_thickness; 1120 if(growthvectors[pi] == 0.0) 1121 continue; 1122 1123 Point<3> points[] = { p10, p10+t*growthvectors[seg1[0]], p11, p11+t*growthvectors[seg1[1]] }; 1124 auto p0 = mesh[pi]; 1125 auto p1 = p0 + t*growthvectors[pi]; 1126 auto P2 = [](Point<3> p) { return Point<2>{p[0], p[1]}; }; 1127 ArrayMem<pair<double, double>, 4> intersections; 1128 1129 double alpha, beta; 1130 1131 if(X_INTERSECTION == intersect( P2(p0), P2(p1), P2(points[0]), P2(points[2]), alpha, beta )) 1132 intersections.Append({alpha, 0.0}); 1133 1134 if(X_INTERSECTION == intersect( P2(p0), P2(p1), P2(points[1]), P2(points[3]), alpha, beta )) 1135 intersections.Append({alpha, 1.0}); 1136 1137 if(X_INTERSECTION == intersect( P2(p0), P2(p1), P2(points[0]), P2(points[1]), alpha, beta )) 1138 intersections.Append({alpha, beta}); 1139 1140 if(X_INTERSECTION == intersect( P2(p0), P2(p1), P2(points[2]), P2(points[3]), alpha, beta )) 1141 intersections.Append({alpha, beta}); 1142 1143 QuickSort(intersections); 1144 for(auto [alpha,beta] : intersections) 1145 { 1146 if(!active_segments.Test(segi1)) 1147 growth[pi] = min(growth[pi], alpha); 1148 else 1149 { 1150 double mean = 0.5*(alpha+beta); 1151 growth[pi] = min(growth[pi], mean); 1152 growth[seg1[0]] = min(growth[seg1[0]], mean); 1153 growth[seg1[1]] = min(growth[seg1[1]], mean); 1154 } 1155 } 1156 } 1157 } 1158 }; 1159 1160 Box<3> box(Box<3>::EMPTY_BOX); 1161 for (auto segi : Range(mesh.LineSegments())) 1162 { 1163 auto segbox = getSegmentBox( segi ); 1164 box.Add(segbox.PMin()); 1165 box.Add(segbox.PMax()); 1166 } 1167 BoxTree<3> segtree(box); 1168 1169 for (auto segi : Range(mesh.LineSegments())) 1170 { 1171 auto p2 = [](Point<3> p) { return Point<2>{p[0], p[1]}; }; 1172 1173 auto seg = mesh[segi]; 1174 double alpha,beta; 1175 intersect( p2(mesh[seg[0]]), p2(mesh[seg[0]]+total_thickness*growthvectors[seg[0]]), p2(mesh[seg[1]]), p2(mesh[seg[1]]+total_thickness*growthvectors[seg[1]]), alpha, beta ); 1176 1177 if(beta>0 && alpha>0 && alpha<1.1) 1178 growth[seg[0]] = min(growth[seg[0]], 0.8*alpha); 1179 if(alpha>0 && beta>0 && beta<1.1) 1180 growth[seg[1]] = min(growth[seg[1]], 0.8*beta); 1181 1182 for (auto segj : Range(mesh.LineSegments())) 1183 if(segi!=segj) 1184 restrictGrowthVectors(segi, segj); 1185 } 1186 1187 for( auto pi : Range(growthvectors)) 1188 growthvectors[pi] *= growth[pi]; 1189 1190 1191 // insert new points 1192 for(PointIndex pi : Range(mesh.Points())) 1193 if(growthvectors[pi].Length2()!=0) 1194 { 1195 1196 auto & pnew = mapto[pi]; 1197 auto dist = 0.0; 1198 for(auto t : thicknesses) 1199 { 1200 dist+=t; 1201 pnew.Append( mesh.AddPoint( mesh[pi] + dist*growthvectors[pi] ) ); 1202 mesh[pnew.Last()].SetType(FIXEDPOINT); 1203 } 1204 } 1205 1206 map<pair<PointIndex, PointIndex>, int> seg2edge; 1207 1208 // insert new elements ( and move old ones ) 1209 for(auto si : moved_segs) 1210 { 1211 auto seg = mesh[si]; 1212 1213 bool swap = false; 1214 auto & pm0 = mapto[seg[0]]; 1215 auto & pm1 = mapto[seg[1]]; 1216 1217 auto newindex = si_map[seg.si]; 1218 1219 Segment s = seg; 1220 s.geominfo[0] = {}; 1221 s.geominfo[1] = {}; 1222 s[0] = pm0.Last(); 1223 s[1] = pm1.Last(); 1224 s[2] = PointIndex::INVALID; 1225 auto pair = s[0] < s[1] ? make_pair(s[0], s[1]) : make_pair(s[1], s[0]); 1226 if(seg2edge.find(pair) == seg2edge.end()) 1227 seg2edge[pair] = ++max_edge_nr; 1228 s.edgenr = seg2edge[pair]; 1229 s.si = seg.si; 1230 mesh.AddSegment(s); 1231 1232 Swap(s[0], s[1]); 1233 s.si = newindex; 1234 mesh.AddSegment(s); 1235 1236 for ( auto i : Range(thicknesses)) 1237 { 1238 PointIndex pi0, pi1, pi2, pi3; 1239 1240 if(i==0) 1241 { 1242 pi0 = seg[0]; 1243 pi1 = seg[1]; 1244 } 1245 else 1246 { 1247 pi0 = pm0[i-1]; 1248 pi1 = pm1[i-1]; 1249 } 1250 1251 pi2 = pm1[i]; 1252 pi3 = pm0[i]; 1253 1254 if(i==0) 1255 { 1256 auto p0 = mesh[pi0]; 1257 auto p1 = mesh[pi1]; 1258 auto q0 = mesh[pi2]; 1259 auto q1 = mesh[pi3]; 1260 1261 Vec<2> n = {-p1[1]+p0[1], p1[0]-p0[0]}; 1262 Vec<2> v = { q0[0]-p0[0], q0[1]-p0[1]}; 1263 if(n[0]*v[0]+n[1]*v[1]<0) 1264 swap = true; 1265 } 1266 1267 Element2d newel; 1268 newel.SetType(QUAD); 1269 newel[0] = pi0; 1270 newel[1] = pi1; 1271 newel[2] = pi2; 1272 newel[3] = pi3; 1273 newel.SetIndex(si_map[seg.si]); 1274 newel.GeomInfo() = PointGeomInfo{}; 1275 1276 // if(swap) 1277 // { 1278 // Swap(newel[0], newel[1]); 1279 // Swap(newel[2], newel[3]); 1280 // } 1281 1282 mesh.AddSurfaceElement(newel); 1283 1284 } 1285 // segment now adjacent to new 2d-domain! 1286 mesh[si].si = si_map[seg.si]; 1287 } 1288 1289 for(auto pi : Range(mapto)) 1290 { 1291 if(mapto[pi].Size() == 0) 1292 continue; 1293 auto pnew = mapto[pi].Last(); 1294 for(auto old_sei : meshtopo.GetVertexSurfaceElements( pi )) 1295 { 1296 if(mesh[old_sei].GetIndex() == domain) 1297 { 1298 auto & old_el = mesh[old_sei]; 1299 for(auto i : IntRange(old_el.GetNP())) 1300 if(old_el[i]==pi) 1301 old_el[i] = pnew; 1302 } 1303 } 1304 } 1305 1306 for(auto & sel : mesh.SurfaceElements()) 1307 if(sel.GetIndex() == domain) 1308 sel.Delete(); 1309 1310 mesh.Compress(); 1311 mesh.CalcSurfacesOfNode(); 1312 1313 Generate2dMesh(mesh, domain); 1314 1315 // even without new domain, we need temporarily a new domain to mesh the remaining area, without confusing the meshes with quads -> add segments temporarily and reset domain number and segments afterwards 1316 if(!should_make_new_domain) 1317 { 1318 // map new domain back to old one 1319 for(auto & sel : mesh.SurfaceElements()) 1320 if(sel.GetIndex()==new_domain) 1321 sel.SetIndex(domain); 1322 1323 // remove (temporary) inner segments 1324 for(auto segi : Range(first_new_seg, mesh.LineSegments().Range().Next())) 1325 { 1326 mesh[segi][0].Invalidate(); 1327 mesh[segi][1].Invalidate(); 1328 } 1329 1330 for(auto segi : moved_segs) 1331 mesh[segi].si = domain; 1332 1333 mesh.Compress(); 1334 mesh.CalcSurfacesOfNode(); 1335 } 1336 1337 return new_domain; 1338 } 1339 1340 } 1341