1 #include <mystdlib.h> 2 #include <meshing.hpp> 3 #include <csg.hpp> 4 5 // #undef DEVELOP 6 // #define DEVELOP 7 8 namespace netgen 9 { 10 11 EdgeCalculation :: EdgeCalculation(const CSGeometry & ageometry,NgArray<SpecialPoint> & aspecpoints,MeshingParameters & amparam)12 EdgeCalculation (const CSGeometry & ageometry, 13 NgArray<SpecialPoint> & aspecpoints, 14 MeshingParameters & amparam) 15 : geometry(ageometry), specpoints(aspecpoints), mparam(amparam) 16 { 17 Box<3> bbox = geometry.BoundingBox(); 18 19 searchtree = new Point3dTree (bbox.PMin(), bbox.PMax()); 20 meshpoint_tree = new Point3dTree (bbox.PMin(), bbox.PMax()); 21 22 for (int i = 0; i < specpoints.Size(); i++) 23 searchtree->Insert (specpoints[i].p, i); 24 25 ideps = 1e-9; 26 } 27 ~EdgeCalculation()28 EdgeCalculation :: ~EdgeCalculation() 29 { 30 delete searchtree; 31 delete meshpoint_tree; 32 } 33 34 Calc(double h,Mesh & mesh)35 void EdgeCalculation :: Calc(double h, Mesh & mesh) 36 { 37 static int timer = NgProfiler::CreateTimer ("CSG: mesh edges"); 38 NgProfiler::RegionTimer reg (timer); 39 40 41 PrintMessage (1, "Find edges"); 42 PushStatus ("Find edges"); 43 44 for (PointIndex pi : mesh.Points().Range()) 45 meshpoint_tree->Insert (mesh[pi], pi); 46 47 48 // add all special points before edge points (important for periodic identification) 49 // JS, Jan 2007 50 const double di=1e-7*geometry.MaxSize(); 51 NgArray<int> locsearch; 52 53 for (int i = 0; i < specpoints.Size(); i++) 54 if (specpoints[i].unconditional) 55 { 56 Point<3> p = specpoints[i].p; 57 meshpoint_tree -> GetIntersecting (p-Vec<3> (di,di,di), 58 p+Vec<3> (di,di,di), locsearch); 59 60 if (locsearch.Size() == 0) 61 { 62 PointIndex pi = mesh.AddPoint (p, specpoints[i].GetLayer(), FIXEDPOINT); 63 meshpoint_tree -> Insert (p, pi); 64 } 65 } 66 67 68 /* 69 // slow version 70 for (int i = 0; i < specpoints.Size(); i++) 71 if (specpoints[i].unconditional) 72 { 73 Point<3> p = specpoints[i].p; 74 bool found = false; 75 for (int j = 1; j <= mesh.GetNP(); j++) 76 if (Dist (p, mesh.Point(j)) < 1e-8) 77 found = true; 78 if (!found) 79 mesh.AddPoint (p, specpoints[i].GetLayer(), FIXEDPOINT); 80 } 81 */ 82 83 84 85 CalcEdges1 (h, mesh); 86 SplitEqualOneSegEdges (mesh); 87 FindClosedSurfaces (h, mesh); 88 PrintMessage (3, cntedge, " edges found"); 89 90 PopStatus (); 91 } 92 93 94 95 96 CalcEdges1(double h,Mesh & mesh)97 void EdgeCalculation :: CalcEdges1 (double h, Mesh & mesh) 98 { 99 NgArray<int> hsp(specpoints.Size()); 100 NgArray<int> glob2hsp(specpoints.Size()); 101 NgArray<int> startpoints, endpoints; 102 103 104 int pos, ep; 105 int layer; 106 107 Point<3> p, np; 108 int pi1, s1, s2, s1_orig, s2_orig; 109 110 NgArray<Point<3> > edgepoints; 111 NgArray<double> curvelength; 112 int copyedge = 0, copyfromedge = -1, copyedgeidentification = -1; 113 114 NgArray<int> locsurfind, locind; 115 116 int checkedcopy = 0; 117 118 // double size = geometry.MaxSize(); 119 // double epspointdist2 = sqr (size) * 1e-12; 120 121 122 // copy special points to work with 123 for (int i = 0; i < specpoints.Size(); i++) 124 { 125 hsp[i] = i; 126 glob2hsp[i] = i; 127 } 128 129 //for(int i=0; i<hsp.Size(); i++) 130 // (*testout) << "hsp["<<i<<"] ... " << specpoints[hsp[i]].p << endl; 131 132 133 cntedge = 0; 134 INDEX_2_HASHTABLE<int> identification_used(100); // identification i already used for startpoint j 135 136 mesh.GetIdentifications().Delete(); 137 138 TABLE<int> specpoint2surface(specpoints.Size()); 139 if (geometry.identifications.Size()) 140 { 141 for (int i = 0; i < specpoints.Size(); i++) 142 for (int j = 0; j < geometry.GetNSurf(); j++) 143 if (geometry.GetSurface(j)->PointOnSurface (specpoints[i].p)) 144 specpoint2surface.Add (i, j); 145 } 146 147 TABLE<int> specpoint2tlo(specpoints.Size()); 148 if (geometry.identifications.Size()) 149 { 150 for (int i = 0; i < specpoints.Size(); i++) 151 for (int j = 0; j < geometry.GetNTopLevelObjects(); j++) 152 { 153 const TopLevelObject * tlo = geometry.GetTopLevelObject (j); 154 if (tlo->GetSolid() && tlo->GetSolid()->VectorIn (specpoints[i].p,specpoints[i].v)) 155 //if (tlo->GetSolid() && tlo->GetSolid()->IsIn (specpoints[i].p)) 156 { 157 #ifdef DEVELOP 158 (*testout) << "point " << specpoints[i].p << " v " <<specpoints[i].v <<" is in " << tlo->GetSolid()->Name() << endl; 159 #endif 160 specpoint2tlo.Add (i, j); 161 } 162 } 163 } 164 165 for (int i = 0; i < specpoints.Size(); i++) 166 specpoints[i].nr = i; 167 168 while (hsp.Size()) 169 { 170 SetThreadPercent(100 - 100 * double (hsp.Size()) / specpoints.Size()); 171 172 #ifdef DEVELOP 173 (*testout) << "hsp.Size() " << hsp.Size() << " specpoints.Size() " << specpoints.Size() << endl; 174 (*testout) << endl << "edge nr " << cntedge+1 << endl; 175 #endif 176 177 edgepoints.SetSize (0); 178 curvelength.SetSize (0); 179 180 181 pi1 = 0; 182 copyedge = 0; 183 // identifyable point available ? 184 185 186 for (int i = 0; i < geometry.identifications.Size() && !pi1; i++) 187 for (int j = checkedcopy; j < startpoints.Size() && !pi1; j++) 188 { 189 #ifdef DEVELOP 190 (*testout) << "checking point " << specpoints[startpoints[j]].p 191 << ", v = " << specpoints[startpoints[j]].v 192 << " for copying (i,j = " << i << ", " << j << ")" << endl; 193 #endif 194 if (geometry.identifications[i]->IdentifyableCandidate (specpoints[startpoints[j]]) && 195 geometry.identifications[i]->IdentifyableCandidate (specpoints[endpoints[j]])) 196 197 198 { 199 int pi1cand = 0; 200 double mindist = 1e10; 201 202 for (int k = 0; k < hsp.Size() && !pi1; k++) 203 { 204 //(*testout) << " ? identifyable with " << specpoints[hsp[k]].p 205 //<< ", v = " << specpoints[hsp[k]].v 206 // << endl; 207 if (identification_used.Used (INDEX_2(i, startpoints[j])) || 208 identification_used.Used (INDEX_2(i, hsp[k]))) 209 { 210 //(*testout) << "failed at pos0" << endl; 211 continue; 212 } 213 214 if (geometry.identifications[i] 215 ->Identifyable(specpoints[startpoints[j]], specpoints[hsp[k]], specpoint2tlo, specpoint2surface) || 216 geometry.identifications[i] 217 ->Identifyable(specpoints[hsp[k]], specpoints[startpoints[j]], specpoint2tlo, specpoint2surface)) 218 { 219 #ifdef DEVELOP 220 (*testout) << "identifyable: " << specpoints[hsp[k]].p << ", v = " << specpoints[hsp[k]].v 221 << " and " << specpoints[startpoints[j]].p << ", v = " << specpoints[startpoints[j]].v 222 << " (identification " << i+1 << ")" << endl; 223 #endif 224 225 if (Dist (specpoints[startpoints[j]].p, specpoints[hsp[k]].p) < mindist) 226 { 227 mindist = Dist (specpoints[startpoints[j]].p, specpoints[hsp[k]].p); 228 pi1cand = k+1; 229 } 230 } 231 } 232 233 234 if (pi1cand) 235 { 236 pi1 = pi1cand; 237 copyedge = 1; 238 copyfromedge = j+1; 239 copyedgeidentification = i+1; 240 241 identification_used.Set (INDEX_2(i, startpoints[j]), 1); 242 identification_used.Set (INDEX_2(i, hsp.Get(pi1)), 1); 243 } 244 } 245 } 246 247 248 // cannot copy from other ege ? 249 if (!pi1) 250 checkedcopy = startpoints.Size(); 251 252 // unconditional special point available ? 253 if (!pi1) 254 for (int i = 1; i <= hsp.Size(); i++) 255 if (specpoints[hsp.Get(i)].unconditional == 1) 256 { 257 pi1 = i; 258 break; 259 } 260 261 262 if (!pi1) 263 { 264 // no unconditional points available, choose first conitional 265 pi1 = 1; 266 } 267 268 layer = specpoints[hsp.Get(pi1)].GetLayer(); 269 270 271 if (!specpoints[hsp.Get(pi1)].unconditional) 272 { 273 specpoints[hsp.Elem(pi1)].unconditional = 1; 274 for (int i = 1; i <= hsp.Size(); i++) 275 if (i != pi1 && 276 Dist (specpoints[hsp.Get(pi1)].p, specpoints[hsp.Get(i)].p) < 1e-8*geometry.MaxSize() && 277 (specpoints[hsp.Get(pi1)].v + specpoints[hsp.Get(i)].v).Length() < 1e-4) 278 { 279 // opposite direction 280 specpoints[hsp.Elem(i)].unconditional = 1; 281 } 282 } 283 284 cntedge++; 285 startpoints.Append (hsp.Get(pi1)); 286 287 #ifdef DEVELOP 288 (*testout) << "start followedge: p1 = " << specpoints[hsp.Get(pi1)].p 289 << ", v = " << specpoints[hsp.Get(pi1)].v << endl; 290 #endif 291 292 FollowEdge (pi1, ep, pos, hsp, h, mesh, 293 edgepoints, curvelength); 294 295 296 if (multithread.terminate) 297 return; 298 299 if (!ep) 300 { 301 // ignore starting point 302 hsp.DeleteElement (pi1); 303 cout << "yes, this happens" << endl; 304 continue; 305 } 306 307 308 309 endpoints.Append (hsp.Get(ep)); 310 311 312 double elen = 0; 313 for (int i = 1; i <= edgepoints.Size()-1; i++) 314 elen += Dist (edgepoints.Get(i), edgepoints.Get(i+1)); 315 316 317 int shortedge = 0; 318 for (int i = 1; i <= geometry.identifications.Size(); i++) 319 if (geometry.identifications.Get(i)->ShortEdge(specpoints[hsp.Get(pi1)], specpoints[hsp.Get(ep)])) 320 shortedge = 1; 321 // (*testout) << "shortedge = " << shortedge << endl; 322 323 324 if (!shortedge) 325 { 326 mesh.RestrictLocalHLine (Point3d (specpoints[hsp.Get(pi1)].p), 327 Point3d (specpoints[hsp.Get(ep)].p), 328 elen / mparam.segmentsperedge); 329 } 330 331 s1 = specpoints[hsp.Get(pi1)].s1; 332 s2 = specpoints[hsp.Get(pi1)].s2; 333 s1_orig = specpoints[hsp.Get(pi1)].s1_orig; 334 s2_orig = specpoints[hsp.Get(pi1)].s2_orig; 335 336 337 // delete initial, terminal and conditional points 338 339 #ifdef DEVELOP 340 (*testout) << "terminal point: p = " << specpoints[hsp.Get(ep)].p 341 << ", v = " << specpoints[hsp.Get(ep)].v << endl; 342 #endif 343 344 searchtree -> DeleteElement (hsp.Get(ep)); 345 searchtree -> DeleteElement (hsp.Get(pi1)); 346 347 if (ep > pi1) 348 { 349 glob2hsp[hsp[ep-1]] = -1; 350 glob2hsp[hsp.Last()] = ep-1; 351 hsp.DeleteElement (ep); 352 353 glob2hsp[hsp[pi1-1]] = -1; 354 glob2hsp[hsp.Last()] = pi1-1; 355 hsp.DeleteElement (pi1); 356 } 357 else 358 { 359 glob2hsp[hsp[pi1-1]] = -1; 360 glob2hsp[hsp.Last()] = pi1-1; 361 hsp.DeleteElement (pi1); 362 363 glob2hsp[hsp[ep-1]] = -1; 364 glob2hsp[hsp.Last()] = ep-1; 365 hsp.DeleteElement (ep); 366 } 367 368 369 for (int j = 1; j <= edgepoints.Size()-1; j++) 370 { 371 p = edgepoints.Get(j); 372 np = Center (p, edgepoints.Get(j+1)); 373 double hd = Dist (p, np); 374 375 376 Box<3> boxp (np - (1.2 * hd) * Vec<3> (1, 1, 1), 377 np + (1.2 * hd) * Vec<3> (1, 1, 1)); 378 searchtree -> GetIntersecting (boxp.PMin(), boxp.PMax(), locind); 379 380 for (int i = 0; i < locind.Size(); i++) 381 { 382 if ( specpoints[locind[i]].HasSurfaces (s1, s2) && 383 specpoints[locind[i]].unconditional == 0) 384 { 385 searchtree -> DeleteElement (locind[i]); 386 387 int li = glob2hsp[locind[i]]; 388 glob2hsp[locind[i]] = -1; 389 glob2hsp[hsp.Last()] = li; 390 hsp.Delete (li); 391 } 392 } 393 394 395 /* 396 for (int i = 1; i <= hsp.Size(); i++) 397 if ( specpoints[hsp.Get(i)].HasSurfaces (s1, s2) && 398 specpoints[hsp.Get(i)].unconditional == 0 && 399 Dist2 (np, specpoints[hsp.Get(i)].p) < 1.2 * hd) 400 { 401 searchtree -> DeleteElement (hsp.Get(i)+1); 402 hsp.DeleteElement (i); 403 i--; 404 } 405 */ 406 } 407 408 409 NgArray<Segment> refedges; 410 NgArray<bool> refedgesinv; 411 412 413 AnalyzeEdge (s1_orig, s2_orig, s1, s2, pos, layer, 414 edgepoints, 415 refedges, refedgesinv); 416 417 418 for (int i = 0; i < refedges.Size(); i++) 419 refedges[i].edgenr = cntedge; 420 421 422 #ifdef DEVELOP 423 (*testout) << "edge " << cntedge << endl 424 << "startp: " << specpoints[startpoints.Last()].p 425 << ", v = " << specpoints[startpoints.Last()].v << endl 426 << "copy = " << copyedge << endl 427 << refedges.Size() << " refedges: "; 428 for (int i = 1; i <= refedges.Size(); i++) 429 (*testout) << " " << refedges.Get(i).si; 430 (*testout) << endl; 431 if (refedgesinv.Size()) 432 (*testout) << "inv[1] = " << refedgesinv.Get(1) << endl; 433 #endif 434 435 if (refedges.Size() == 0) 436 throw NgException ("Problem in edge detection"); 437 438 439 if (!copyedge) 440 { 441 // (*testout) << "store edge" << endl; 442 // int oldnseg = mesh.GetNSeg(); 443 444 if (!shortedge) 445 StoreEdge (refedges, refedgesinv, 446 edgepoints, curvelength, layer, mesh); 447 else 448 StoreShortEdge (refedges, refedgesinv, 449 edgepoints, curvelength, layer, mesh); 450 451 452 for(int i = 0; i < refedges.Size(); i++) 453 { 454 refedges[i].surfnr1 = geometry.GetSurfaceClassRepresentant(refedges[i].surfnr1); 455 refedges[i].surfnr2 = geometry.GetSurfaceClassRepresentant(refedges[i].surfnr2); 456 } 457 458 /* 459 for (int i = oldnseg+1; i <= mesh.GetNSeg(); i++) 460 for (int j = 1; j <= oldnseg; j++) 461 { 462 const Point<3> & l1p1 = mesh.Point (mesh.LineSegment(i).p1); 463 const Point<3> & l1p2 = mesh.Point (mesh.LineSegment(i).p2); 464 const Point<3> & l2p1 = mesh.Point (mesh.LineSegment(j).p1); 465 const Point<3> & l2p2 = mesh.Point (mesh.LineSegment(j).p2); 466 Vec<3> vl1(l1p1, l1p2); 467 for (double lamk = 0; lamk <= 1; lamk += 0.1) 468 { 469 Point<3> l2p = l1p1 + lamk * vl1; 470 double dist = sqrt (MinDistLP2 (l2p1, l2p2, l2p)); 471 if (dist > 1e-12) 472 mesh.RestrictLocalH (l2p, 3*dist); 473 } 474 } 475 */ 476 } 477 else 478 { 479 CopyEdge (refedges, refedgesinv, 480 copyfromedge, 481 specpoints[startpoints.Get(copyfromedge)].p, 482 specpoints[endpoints.Get(copyfromedge)].p, 483 edgepoints.Get(1), edgepoints.Last(), 484 copyedgeidentification, 485 layer, 486 mesh); 487 } 488 489 490 { 491 // named edge ? 492 // cout << "check edge name, size = " << geometry.named_edges.size() << endl; 493 // for (auto pair : geometry.named_edges) 494 // cout << "key = " << get<0> (pair.first) << "-" << get<1> (pair.first) << ", val = " << pair.second << endl; 495 496 Surface * sp1 = const_cast<Surface*> (geometry.GetSurface(s1)); 497 Surface * sp2 = const_cast<Surface*> (geometry.GetSurface(s2)); 498 // cout << "sp1 = " << sp1 << ", sp2 = " << sp2 << endl; 499 500 auto ptr = geometry.named_edges.find(tuple(sp1, sp2)); 501 if (ptr != geometry.named_edges.end()) 502 for (int i = 0; i < refedges.Size(); i++) 503 mesh.SetCD2Name(refedges[i].edgenr, ptr->second); 504 505 ptr = geometry.named_edges.find(tuple(sp2, sp1)); 506 if (ptr != geometry.named_edges.end()) 507 for (int i = 0; i < refedges.Size(); i++) 508 mesh.SetCD2Name(refedges[i].edgenr, ptr->second); 509 } 510 511 for(int i=0; i<refedges.Size(); i++) 512 { 513 auto splinesurface = dynamic_cast<const SplineSurface*>(geometry.GetSurface(refedges[i].surfnr1)); 514 if(splinesurface) 515 { 516 auto name = splinesurface->GetBCNameOf(specpoints[startpoints.Get(refedges[i].edgenr)].p,specpoints[endpoints.Get(refedges[i].edgenr)].p); 517 mesh.SetCD2Name(refedges[i].edgenr,name); 518 } 519 else 520 { 521 auto splinesurface2 = dynamic_cast<const SplineSurface*>(geometry.GetSurface(refedges[i].surfnr2)); 522 if(splinesurface2) 523 { 524 auto name = splinesurface2->GetBCNameOf(specpoints[startpoints.Get(refedges[i].edgenr)].p,specpoints[endpoints.Get(refedges[i].edgenr)].p); 525 mesh.SetCD2Name(refedges[i].edgenr,name); 526 } 527 528 } 529 530 } 531 532 /* 533 // not available ... 534 for (int i = 0; i < refedges.Size(); i++) 535 { 536 EdgeDescriptor ed; 537 ed.SetSurfNr(0, refedges[i].surfnr1); 538 ed.SetSurfNr(1, refedges[i].surfnr2); 539 int hnr = mesh.AddEdgeDescriptor(ed); 540 if (hnr != refedges[i].edgenr) 541 { 542 cerr << "edgedescriptor index wrong: new : " << hnr << " old = " << refedges[i].edgenr << endl; 543 } 544 } 545 */ 546 547 548 549 // for(int i=0; i<hsp.Size(); i++) 550 // { 551 // (*testout) << "pos2 hsp["<<i<<"] ... " << specpoints[hsp[i]].p << endl; 552 // } 553 } 554 } 555 556 557 558 559 560 561 /* 562 If two or more edges share the same initial and end-points, 563 then they need at least two segments 564 */ 565 void EdgeCalculation :: SplitEqualOneSegEdges(Mesh & mesh)566 SplitEqualOneSegEdges (Mesh & mesh) 567 { 568 // int i, j; 569 SegmentIndex si; 570 PointIndex pi; 571 572 NgArray<int> osedges(cntedge); 573 INDEX_2_HASHTABLE<int> osedgesht (cntedge+1); 574 575 osedges = 2; 576 577 // count segments on edges 578 for (si = 0; si < mesh.GetNSeg(); si++) 579 { 580 const Segment & seg = mesh[si]; 581 if (seg.seginfo && seg.edgenr >= 1 && seg.edgenr <= cntedge) 582 osedges.Elem(seg.edgenr)--; 583 } 584 585 // flag one segment edges 586 for (int i = 0; i < cntedge; i++) 587 osedges[i] = (osedges[i] > 0) ? 1 : 0; 588 589 for (si = 0; si < mesh.GetNSeg(); si++) 590 { 591 const Segment & seg = mesh[si]; 592 if (seg.seginfo && seg.edgenr >= 1 && seg.edgenr <= cntedge) 593 { 594 if (osedges.Get(seg.edgenr)) 595 { 596 INDEX_2 i2(seg[0], seg[1]); 597 i2.Sort (); 598 if (osedgesht.Used (i2)) 599 osedgesht.Set (i2, 2); 600 else 601 osedgesht.Set (i2, 1); 602 } 603 } 604 } 605 606 607 // one edge 1 segment, other 2 segments 608 // yes, it happens ! 609 point_on_edge_problem = 0; 610 for (int i = 1; i <= osedgesht.GetNBags(); i++) 611 for (int j = 1; j <= osedgesht.GetBagSize(i); j++) 612 { 613 INDEX_2 i2; 614 int val; 615 osedgesht.GetData (i, j, i2, val); 616 617 const Point<3> & p1 = mesh[PointIndex(i2.I1())]; 618 const Point<3> & p2 = mesh[PointIndex(i2.I2())]; 619 Vec<3> v = p2 - p1; 620 double vlen = v.Length(); 621 v /= vlen; 622 for (pi = PointIndex::BASE; 623 pi < mesh.GetNP()+PointIndex::BASE; pi++) 624 625 if (pi != i2.I1() && pi != i2.I2()) 626 { 627 const Point<3> & p = mesh[pi]; 628 Vec<3> v2 = p - p1; 629 double lam = (v2 * v); 630 if (lam > 0 && lam < vlen) 631 { 632 Point<3> hp = p1 + lam * v; 633 if (Dist (p, hp) < 1e-4 * vlen) 634 { 635 PrintWarning ("Point on edge !!!"); 636 cout << "seg: " << i2 << ", p = " << pi << endl; 637 osedgesht.Set (i2, 2); 638 point_on_edge_problem = 1; 639 640 (*testout) << "Point on edge" << endl 641 << "seg = " << i2 << ", p = " << pi << endl 642 << "pos = " << p << ", projected = " << hp << endl 643 << "seg is = " 644 << mesh.Point(PointIndex(i2.I1())) << " - " 645 << mesh.Point(PointIndex(i2.I2())) << endl; 646 } 647 } 648 } 649 } 650 651 652 // insert new points 653 osedges = -1; 654 655 int nseg = mesh.GetNSeg(); 656 for (si = 0; si < nseg; si++) 657 { 658 const Segment & seg = mesh[si]; 659 if (seg.seginfo && seg.edgenr >= 1 && seg.edgenr <= cntedge) 660 { 661 INDEX_2 i2(seg[0], seg[1]); 662 i2.Sort (); 663 if (osedgesht.Used (i2) && 664 osedgesht.Get (i2) == 2 && 665 osedges.Elem(seg.edgenr) == -1) 666 { 667 Point<3> newp = Center (mesh[PointIndex(seg[0])], 668 mesh[PointIndex(seg[1])]); 669 670 ProjectToEdge (geometry.GetSurface(seg.surfnr1), 671 geometry.GetSurface(seg.surfnr2), 672 newp); 673 674 osedges.Elem(seg.edgenr) = 675 mesh.AddPoint (newp, mesh[PointIndex(seg[0])].GetLayer(), EDGEPOINT); 676 meshpoint_tree -> Insert (newp, osedges.Elem(seg.edgenr)); 677 } 678 } 679 } 680 681 682 // for (int i = 1; i <= nseg; i++) 683 for (Segment & seg : mesh.LineSegments()) 684 { 685 // Segment & seg = mesh.LineSegment (i); 686 if (seg.edgenr >= 1 && seg.edgenr <= cntedge) 687 { 688 if (osedges.Get(seg.edgenr) != -1) 689 { 690 Segment newseg = seg; 691 newseg[0] = osedges.Get(seg.edgenr); 692 seg[1] = osedges.Get(seg.edgenr); 693 mesh.AddSegment (newseg); 694 } 695 } 696 } 697 698 } 699 700 701 702 void EdgeCalculation :: FollowEdge(int pi1,int & ep,int & pos,const NgArray<int> & hsp,double h,const Mesh & mesh,NgArray<Point<3>> & edgepoints,NgArray<double> & curvelength)703 FollowEdge (int pi1, int & ep, int & pos, 704 const NgArray<int> & hsp, 705 double h, const Mesh & mesh, 706 NgArray<Point<3> > & edgepoints, 707 NgArray<double> & curvelength) 708 { 709 int s1, s2, s1_rep, s2_rep; 710 double len, steplen, cursteplen, loch; 711 Point<3> p, np, pnp; 712 Vec<3> a1, a2, t; 713 714 NgArray<int> locind; 715 716 double size = geometry.MaxSize(); 717 double epspointdist2 = size * 1e-6; 718 epspointdist2 = sqr (epspointdist2); 719 int uselocalh = mparam.uselocalh; 720 721 722 s1_rep = specpoints[hsp.Get(pi1)].s1; 723 s2_rep = specpoints[hsp.Get(pi1)].s2; 724 s1 = specpoints[hsp.Get(pi1)].s1_orig; 725 s2 = specpoints[hsp.Get(pi1)].s2_orig; 726 727 p = specpoints[hsp.Get(pi1)].p; 728 //ProjectToEdge (geometry.GetSurface(s1), 729 // geometry.GetSurface(s2), p); 730 geometry.GetSurface(s1) -> CalcGradient (p, a1); 731 geometry.GetSurface(s2) -> CalcGradient (p, a2); 732 733 t = Cross (a1, a2); 734 t.Normalize(); 735 736 pos = (specpoints[hsp.Get(pi1)].v * t) > 0; 737 if (!pos) t *= -1; 738 739 740 edgepoints.Append (p); 741 curvelength.Append (0); 742 len = 0; 743 744 // (*testout) << "geometry.GetSurface(s1) -> LocH (p, 3, 1, h) " << geometry.GetSurface(s1) -> LocH (p, 3, 1, h) 745 // << " geometry.GetSurface(s2) -> LocH (p, 3, 1, h) " << geometry.GetSurface(s2) -> LocH (p, 3, 1, h) << endl; 746 747 loch = min2 (geometry.GetSurface(s1) -> LocH (p, 3, 1, mparam, h), 748 geometry.GetSurface(s2) -> LocH (p, 3, 1, mparam, h)); 749 750 751 752 if (uselocalh) 753 { 754 double lh = mesh.GetH(p); 755 // (*testout) << "lh " << lh << endl; 756 if (lh < loch) 757 loch = lh; 758 } 759 760 steplen = 0.1 * loch; 761 762 do 763 { 764 if (multithread.terminate) 765 return; 766 767 if (fabs (p(0)) + fabs (p(1)) + fabs (p(2)) > 100000*size) 768 { 769 ep = 0; 770 PrintWarning ("Give up line"); 771 break; 772 } 773 774 if (steplen > 0.1 * loch) steplen = 0.1 * loch; 775 776 steplen *= 2; 777 do 778 { 779 steplen *= 0.5; 780 np = p + steplen * t; 781 pnp = np; 782 ProjectToEdge (geometry.GetSurface(s1), 783 geometry.GetSurface(s2), pnp); 784 } 785 while (Dist (np, pnp) > 0.1 * steplen); 786 787 788 cursteplen = steplen; 789 if (Dist (np, pnp) < 0.01 * steplen) steplen *= 2; 790 791 792 np = pnp; 793 ep = 0; 794 795 double hvtmin = 1.5 * cursteplen; 796 797 Box<3> boxp (p - (2 * cursteplen) * Vec<3> (1, 1, 1), 798 p + (2 * cursteplen) * Vec<3> (1, 1, 1)); 799 800 searchtree -> GetIntersecting (boxp.PMin(), boxp.PMax(), locind); 801 802 for (int i = 0; i < locind.Size(); i++) 803 { 804 Vec<3> hv = specpoints[locind[i]].p - p; 805 if (hv.Length2() > 9 * cursteplen * cursteplen) 806 continue; 807 808 double hvt = hv * t; 809 hv -= hvt * t; 810 811 if (hv.Length() < 0.2 * cursteplen && 812 hvt > 0 && 813 // hvt < 1.5 * cursteplen && 814 hvt < hvtmin && 815 specpoints[locind[i]].unconditional == 1 && 816 (specpoints[locind[i]].v + t).Length() < 0.4 ) 817 { 818 Point<3> hep = specpoints[locind[i]].p; 819 ProjectToEdge (geometry.GetSurface(s1), 820 geometry.GetSurface(s2), hep); 821 822 823 if (Dist2 (hep, specpoints[locind[i]].p) < epspointdist2 ) 824 { 825 geometry.GetSurface(s1) -> CalcGradient (hep, a1); 826 geometry.GetSurface(s2) -> CalcGradient (hep, a2); 827 Vec<3> ept = Cross (a1, a2); 828 ept /= ept.Length(); 829 if (!pos) ept *= -1; 830 831 if ( (specpoints[locind[i]].v + ept).Length() < 1e-4 ) 832 { 833 np = specpoints[locind[i]].p; 834 835 for (int jj = 0; jj < hsp.Size(); jj++) 836 if (hsp[jj] == locind[i]) 837 ep = jj+1; 838 839 if (!ep) 840 cerr << "endpoint not found" << endl; 841 // ep = i; 842 hvtmin = hvt; 843 // break; 844 } 845 } 846 } 847 } 848 849 850 851 852 /* 853 for (int i = 1; i <= hsp.Size(); i++) 854 { 855 if (!boxp.IsIn (specpoints[hsp.Get(i)].p)) 856 continue; 857 858 Vec<3> hv = specpoints[hsp.Get(i)].p - p; 859 if (hv.Length2() > 9 * cursteplen * cursteplen) 860 continue; 861 862 double hvt = hv * t; 863 hv -= hvt * t; 864 865 if (hv.Length() < 0.2 * cursteplen && 866 hvt > 0 && 867 // hvt < 1.5 * cursteplen && 868 hvt < hvtmin && 869 specpoints[hsp.Get(i)].unconditional == 1 && 870 (specpoints[hsp.Get(i)].v + t).Length() < 0.4 ) 871 { 872 Point<3> hep = specpoints[hsp.Get(i)].p; 873 ProjectToEdge (geometry.GetSurface(s1), 874 geometry.GetSurface(s2), hep); 875 876 877 if (Dist2 (hep, specpoints[hsp.Get(i)].p) < epspointdist2 ) 878 { 879 geometry.GetSurface(s1) -> CalcGradient (hep, a1); 880 geometry.GetSurface(s2) -> CalcGradient (hep, a2); 881 Vec<3> ept = Cross (a1, a2); 882 ept /= ept.Length(); 883 if (!pos) ept *= -1; 884 885 if ( (specpoints[hsp.Get(i)].v + ept).Length() < 1e-4 ) 886 { 887 np = specpoints[hsp.Get(i)].p; 888 ep = i; 889 hvtmin = hvt; 890 // break; 891 } 892 } 893 } 894 } 895 */ 896 897 loch = min2 (geometry.GetSurface(s1_rep) -> LocH (np, 3, 1, mparam, h), 898 geometry.GetSurface(s2_rep) -> LocH (np, 3, 1, mparam, h)); 899 loch = max2 (loch, mparam.minh); 900 901 if (uselocalh) 902 { 903 double lh = mesh.GetH(np); 904 if (lh < loch) loch = lh; 905 } 906 907 len += Dist (p, np) / loch; 908 edgepoints.Append (np); 909 curvelength.Append (len); 910 911 p = np; 912 913 geometry.GetSurface(s1) -> CalcGradient (p, a1); 914 geometry.GetSurface(s2) -> CalcGradient (p, a2); 915 t = Cross (a1, a2); 916 t.Normalize(); 917 if (!pos) t *= -1; 918 } 919 while (! ep); 920 } 921 922 923 924 925 926 927 928 void EdgeCalculation :: AnalyzeEdge(int s1,int s2,int s1_rep,int s2_rep,int pos,int layer,const NgArray<Point<3>> & edgepoints,NgArray<Segment> & refedges,NgArray<bool> & refedgesinv)929 AnalyzeEdge (int s1, int s2, int s1_rep, int s2_rep, int pos, int layer, 930 const NgArray<Point<3> > & edgepoints, 931 NgArray<Segment> & refedges, 932 NgArray<bool> & refedgesinv) 933 { 934 Segment seg; 935 NgArray<int> locsurfind, locsurfind2; 936 937 NgArray<int> edges_priority; 938 939 double size = geometry.MaxSize(); 940 bool debug = 0; 941 942 #ifdef DEVELOP 943 debug = 1; 944 #endif 945 946 if (debug) 947 { 948 (*testout) << "debug edge !!!" << endl; 949 (*testout) << "edgepoints = " << edgepoints << endl; 950 (*testout) << "s1, s2 = " << s1 << " - " << s2 << endl; 951 (*testout) << "s1_rep, s2_rep = " << s1_rep << " - " << s2_rep << endl; 952 } 953 954 refedges.SetSize(0); 955 refedgesinv.SetSize(0); 956 Point<3> hp = Center (edgepoints[0], edgepoints[1]); 957 ProjectToEdge (geometry.GetSurface(s1), geometry.GetSurface(s2), hp); 958 959 if (debug) 960 *testout << "hp = " << hp << endl; 961 962 Vec<3> t, a1, a2; 963 geometry.GetSurface(s1) -> CalcGradient (hp, a1); 964 geometry.GetSurface(s2) -> CalcGradient (hp, a2); 965 t = Cross (a1, a2); 966 t.Normalize(); 967 if (!pos) t *= -1; 968 969 970 971 for (int i = 0; i < geometry.GetNTopLevelObjects(); i++) 972 { 973 // Solid * locsol; 974 975 if (geometry.GetTopLevelObject(i)->GetLayer() != layer) 976 continue; 977 978 const Solid * sol = geometry.GetTopLevelObject(i)->GetSolid(); 979 const Surface * surf = geometry.GetTopLevelObject(i)->GetSurface(); 980 981 // sol -> TangentialSolid (hp, locsol, locsurfind, size*ideps); 982 auto locsol = sol -> TangentialSolid (hp, locsurfind, size*ideps); 983 984 //*testout << "hp = " << hp << endl; 985 //(*testout) << "locsol: " << endl; 986 //if (locsol) locsol->Print(*testout); 987 //(*testout) << endl; 988 989 990 if (!locsol) continue; 991 992 BoxSphere<3> boxp (hp, hp); 993 boxp.Increase (1e-8*size); 994 boxp.CalcDiamCenter(); 995 996 ReducePrimitiveIterator rpi(boxp); 997 UnReducePrimitiveIterator urpi; 998 999 // ((Solid*)locsol) -> IterateSolid (rpi); 1000 locsol -> IterateSolid (rpi); 1001 1002 locsol -> CalcSurfaceInverse (); 1003 1004 if (!surf) 1005 { 1006 locsol -> GetTangentialSurfaceIndices (hp,locsurfind,ideps*size); 1007 } 1008 else 1009 { 1010 /* 1011 if (fabs (surf->CalcFunctionValue (hp)) < 1e-6) 1012 continue; 1013 */ 1014 locsurfind.SetSize(1); 1015 locsurfind[0] = -1; 1016 for (int j = 0; j < geometry.GetNSurf(); j++) 1017 if (geometry.GetSurface(j) == surf) 1018 { 1019 locsurfind[0] = j; 1020 // geometry.GetSurfaceClassRepresentant(j); 1021 break; 1022 } 1023 } 1024 1025 // ((Solid*)locsol) -> IterateSolid (urpi); 1026 locsol -> IterateSolid (urpi); 1027 1028 1029 if (debug) 1030 (*testout) << "edge of tlo " << i << ", has " << locsurfind.Size() << " faces." << endl; 1031 1032 1033 for (int j = locsurfind.Size()-1; j >= 0; j--) 1034 if (fabs (geometry.GetSurface(locsurfind[j]) 1035 ->CalcFunctionValue (hp) ) > ideps*size) 1036 locsurfind.Delete(j); 1037 1038 if (debug) 1039 (*testout) << locsurfind.Size() << " faces on hp" << endl; 1040 1041 1042 1043 for (int j = 0; j < locsurfind.Size(); j++) 1044 { 1045 int lsi = locsurfind[j]; 1046 int rlsi = geometry.GetSurfaceClassRepresentant(lsi); 1047 1048 1049 // n is outer normal to solid 1050 Vec<3> n = geometry.GetSurface(lsi) -> GetNormalVector (hp); 1051 if (debug) 1052 *testout << "n1 = " << n << endl; 1053 if (geometry.GetSurface (lsi)->Inverse()) 1054 n *= -1; 1055 1056 if (fabs (t * n) > 1e-4) continue; 1057 if (debug) 1058 { 1059 (*testout) << "face " << locsurfind[j] << ", rep = " << rlsi 1060 << " has (t*n) = " << (t*n) << endl; 1061 (*testout) << "n = " << n << endl; 1062 } 1063 1064 // rn is normal to class representant 1065 Vec<3> rn = geometry.GetSurface(rlsi) -> GetNormalVector (hp); 1066 if (debug) 1067 { 1068 (*testout) << "rn = " << rn << endl; 1069 } 1070 1071 //if( n*rn < 0) 1072 // rn *= -1; 1073 1074 bool sameasref = ((n * rn) > 0); 1075 1076 //m = Cross (t, rn); 1077 Vec<3> m = Cross (t, n); 1078 if(!sameasref) m*=-1.; 1079 1080 m.Normalize(); 1081 1082 1083 if (debug) 1084 (*testout) << "m = " << m << endl; 1085 1086 1087 //bool founddirection = false; 1088 //int k; 1089 double eps = 1e-8*size; 1090 1091 ArrayMem<bool,2> pre_ok(2); 1092 bool flip = false; 1093 1094 do 1095 { 1096 eps *= 0.5; 1097 auto in00 = locsol -> VectorIn2 (hp, m, n, eps); 1098 auto in01 = locsol -> VectorIn2 (hp, m, -1. * n, eps); 1099 pre_ok[0] = in00 == IS_OUTSIDE && in01 == IS_INSIDE; 1100 1101 if(in00 == IS_INSIDE && in01 == IS_OUTSIDE) 1102 pre_ok[0] = flip = true; 1103 1104 auto in10 = locsol -> VectorIn2 (hp, -1.*m, n, eps); 1105 auto in11 = locsol -> VectorIn2 (hp, -1.*m, -1. * n, eps); 1106 pre_ok[1] = (in10 == IS_OUTSIDE && in11 == IS_INSIDE); 1107 1108 if(in10 == IS_INSIDE && in11 == IS_OUTSIDE) 1109 pre_ok[1] = flip = true; 1110 1111 if (debug) 1112 { 1113 *testout << "eps = " << eps << endl; 1114 *testout << "in,1 = " << in00 << endl; 1115 *testout << "in,1 = " << in01 << endl; 1116 *testout << "in,1 = " << in10 << endl; 1117 *testout << "in,1 = " << in11 << endl; 1118 } 1119 } 1120 while(pre_ok[0] && pre_ok[1] && eps > 1e-16*size); 1121 1122 if (debug) 1123 { 1124 *testout << "eps = " << eps << ", size = " << size << endl; 1125 *testout << "pre_ok[0,1] = " << pre_ok[0] << "," << pre_ok[1] << endl; 1126 } 1127 1128 eps = 1e-8*size; 1129 1130 1131 for (int k = 1; k <= 2; k ++) 1132 { 1133 bool edgeinv = (k == 2); 1134 1135 if (debug) 1136 { 1137 (*testout) << "onface(" << hp << ", " << m << ")= " << flush; 1138 (*testout) << locsol->OnFace (hp, m, eps) << flush; 1139 (*testout) << " n " << n << flush; 1140 (*testout) << " vec2in = " 1141 << locsol -> VectorIn2 (hp, m, n, eps) << " and " 1142 << locsol -> VectorIn2 (hp, m, -1 * n, eps) << endl; 1143 } 1144 1145 // if (locsol -> OnFace (hp, m)) 1146 1147 1148 // one side must be inside, the other must be outside 1149 bool ok = (pre_ok[k-1] || 1150 (locsol -> VectorIn2 (hp, m, n, eps) == IS_OUTSIDE && 1151 locsol -> VectorIn2 (hp, m, -1 * n, eps) == IS_INSIDE)); 1152 1153 if (debug) 1154 (*testout) << "ok (before) " << ok << endl; 1155 1156 // compute second order approximation 1157 // curve = hp + t m + t*t/2 m2 1158 Vec<3> grad, m2; 1159 Mat<3> hesse; 1160 geometry.GetSurface(lsi) -> CalcGradient (hp, grad); 1161 geometry.GetSurface(lsi) -> CalcHesse (hp, hesse); 1162 double fac = -(m * (hesse * m)) / (grad * grad); 1163 m2 = fac * grad; 1164 // (*testout) << "hp = " << hp << ", m = " << m << ", m2 = " << m2 << endl; 1165 1166 // Solid * locsol2; 1167 auto locsol2 = locsol -> TangentialSolid3 (hp, m, m2, locsurfind2, ideps*size); 1168 if (!locsol2) ok = 0; 1169 // delete locsol2; 1170 1171 1172 if (ok) 1173 { 1174 if (debug) 1175 (*testout) << "is true" << endl; 1176 int hi = 0; 1177 for (int l = 1; !hi && l <= refedges.Size(); l++) 1178 { 1179 if (refedges.Get(l).si == rlsi && // JS sept 2006 1180 // if (refedges.Get(l).si == lsi && 1181 refedgesinv.Get(l) == edgeinv) 1182 { 1183 hi = l; 1184 } 1185 } 1186 1187 if (!hi) 1188 { 1189 seg.si = rlsi; // JS Sept 2006 1190 // seg.si = lsi; 1191 seg.domin = -1; 1192 seg.domout = -1; 1193 seg.tlosurf = -1; 1194 //seg.surfnr1 = s1_rep; 1195 //seg.surfnr2 = s2_rep; 1196 seg.surfnr1 = s1; 1197 seg.surfnr2 = s2; 1198 refedges.Append (seg); 1199 hi = refedges.Size(); 1200 refedgesinv.Append (edgeinv); 1201 edges_priority.Append((pre_ok[k-1]) ? 1 : 0); 1202 } 1203 else 1204 { 1205 if(edges_priority[hi-1] / 10 == -i-1) 1206 edges_priority[hi-1] = 10*(i+1); 1207 else 1208 edges_priority[hi-1] = -10*(i+1); 1209 } 1210 1211 if (!surf) 1212 { 1213 bool inside = sameasref; 1214 if(flip) 1215 inside = !inside; 1216 if (inside) 1217 refedges.Elem(hi).domin = i; 1218 else 1219 refedges.Elem(hi).domout = i; 1220 } 1221 else 1222 { 1223 refedges.Elem(hi).tlosurf = i; 1224 for(int kk = 0; kk < geometry.GetNTopLevelObjects(); kk++) 1225 { 1226 auto othersolid = geometry.GetTopLevelObject(kk)->GetSolid(); 1227 auto othersurf = geometry.GetTopLevelObject(kk)->GetSurface(); 1228 if(!othersurf && dynamic_cast<SplineSurface*>(othersurf)) 1229 { 1230 if(othersolid->IsIn(edgepoints[0]) && 1231 othersolid->IsIn(edgepoints[edgepoints.Size()-1])) 1232 { 1233 refedges.Elem(hi).domin = kk; 1234 refedges.Elem(hi).domout = kk; 1235 } 1236 } 1237 } 1238 } 1239 1240 if(pre_ok[k-1]) 1241 edges_priority[hi-1] = 1; 1242 1243 1244 if (debug) 1245 (*testout) << "add ref seg:" 1246 << "si = " << refedges.Get(hi).si 1247 << ", domin = " << refedges.Get(hi).domin 1248 << ", domout = " << refedges.Get(hi).domout 1249 << ", surfnr1/2 = " << refedges.Get(hi).surfnr1 1250 << ", " << refedges.Get(hi).surfnr2 1251 << ", inv = " << refedgesinv.Get(hi) 1252 << ", refedgenr = " << hi 1253 << ", priority = " << edges_priority[hi-1] 1254 << ", hi = " << hi 1255 << endl; 1256 } 1257 else 1258 { 1259 if (debug) 1260 (*testout) << "is false" << endl; 1261 } 1262 m *= -1; 1263 } 1264 } 1265 // delete locsol; 1266 } 1267 1268 1269 if (debug) 1270 { 1271 *testout << "Refsegments, before delete: " << endl << refedges << endl; 1272 *testout << "inv: " << endl << refedgesinv << endl; 1273 } 1274 1275 if(refedges.Size() == 0) 1276 throw Exception("No edges found, something wrong."); 1277 1278 NgBitArray todelete(refedges.Size()); 1279 todelete.Clear(); 1280 1281 1282 for(int i=0; i<refedges.Size()-1; i++) 1283 { 1284 for(int j=i+1; !todelete.Test(i) && j<refedges.Size(); j++) 1285 { 1286 if(todelete.Test(j)) 1287 continue; 1288 1289 if(refedges[i].si == refedges[j].si && 1290 refedges[i].domin == refedges[j].domin && 1291 refedges[i].domout == refedges[j].domout && 1292 geometry.GetSurfaceClassRepresentant(refedges[i].surfnr1) == geometry.GetSurfaceClassRepresentant(refedges[j].surfnr1) && 1293 geometry.GetSurfaceClassRepresentant(refedges[i].surfnr2) == geometry.GetSurfaceClassRepresentant(refedges[j].surfnr2) 1294 // && refedgesinv[i] == refedgesinv[j] // JS, 20060802 1295 ) 1296 { 1297 if(debug) 1298 (*testout) << "equal segments: " << refedges[i] << " pri " << edges_priority[i] 1299 << " tlosurf " << refedges[i].tlosurf 1300 << "\n and " << refedges[j] << " pri " << edges_priority[j] 1301 << " tlosurf " << refedges[i].tlosurf << endl; 1302 1303 if(edges_priority[i] < 10 && edges_priority[i] < edges_priority[j]) 1304 { 1305 todelete.Set(i); 1306 } 1307 else if (edges_priority[j] < 10 && edges_priority[i] > edges_priority[j]) 1308 { 1309 todelete.Set(j); 1310 } 1311 } 1312 } 1313 1314 } 1315 1316 int num = refedges.Size(); 1317 1318 for(int i=refedges.Size()-1; num>2 && i>=0; i--) 1319 if(todelete.Test(i)) 1320 { 1321 refedges.Delete(i); 1322 refedgesinv.Delete(i); 1323 num--; 1324 } 1325 1326 1327 if (debug) 1328 { 1329 *testout << "Refsegments: " << endl << refedges << endl; 1330 } 1331 } 1332 1333 1334 1335 void EdgeCalculation :: StoreEdge(const NgArray<Segment> & refedges,const NgArray<bool> & refedgesinv,const NgArray<Point<3>> & edgepoints,const NgArray<double> & curvelength,int layer,Mesh & mesh)1336 StoreEdge (const NgArray<Segment> & refedges, 1337 const NgArray<bool> & refedgesinv, 1338 const NgArray<Point<3> > & edgepoints, 1339 const NgArray<double> & curvelength, 1340 int layer, 1341 Mesh & mesh) 1342 { 1343 1344 // Calculate optimal element-length 1345 int i, j, k; 1346 // PointIndex pi; 1347 int ne; 1348 1349 double len, corr, lam; 1350 PointIndex thispi, lastpi; 1351 Point<3> p, np; 1352 Segment seg; 1353 1354 const Surface * surf1 = geometry.GetSurface (refedges.Get(1).surfnr1); 1355 const Surface * surf2 = geometry.GetSurface (refedges.Get(1).surfnr2); 1356 1357 (*testout) << "s1 " << refedges.Get(1).surfnr1 << " s2 " << refedges.Get(1).surfnr2 1358 << " rs1 " << geometry.GetSurfaceClassRepresentant(refedges.Get(1).surfnr1) 1359 << " rs2 " << geometry.GetSurfaceClassRepresentant(refedges.Get(1).surfnr2) << endl; 1360 1361 len = curvelength.Last(); 1362 ne = int (len + 0.5); 1363 if (ne == 0) ne = 1; 1364 if (Dist (edgepoints.Get(1), edgepoints.Last()) < 1e-8*geometry.MaxSize() && 1365 ne <= 6) 1366 ne = 6; 1367 corr = len / ne; 1368 1369 // generate initial point 1370 p = edgepoints.Get(1); 1371 lastpi = PointIndex::INVALID; 1372 1373 /* 1374 for (pi = PointIndex::BASE; 1375 pi < mesh.GetNP()+PointIndex::BASE; pi++) 1376 if (Dist (mesh[pi], p) < 1e-6) 1377 { 1378 lastpi = pi; 1379 break; 1380 } 1381 */ 1382 1383 const double di=1e-7*geometry.MaxSize(); 1384 1385 NgArray<int> locsearch; 1386 meshpoint_tree -> GetIntersecting (p-Vec<3> (di,di,di), 1387 p+Vec<3> (di,di,di), locsearch); 1388 if (locsearch.Size()) 1389 lastpi = locsearch[0]; 1390 1391 1392 1393 if (!lastpi.IsValid()) 1394 { 1395 lastpi = mesh.AddPoint (p, layer, FIXEDPOINT); 1396 meshpoint_tree -> Insert (p, lastpi); 1397 // (*testout) << "test1, store point " << lastpi << ", p = " << p << endl; 1398 } 1399 1400 j = 1; 1401 for (i = 1; i <= ne; i++) 1402 { 1403 while (curvelength.Get(j) < i * corr && j < curvelength.Size()) j++; 1404 1405 1406 lam = (i * corr - curvelength.Get(j-1)) / 1407 (curvelength.Get(j) - curvelength.Get(j-1)); 1408 1409 np(0) = (1-lam) * edgepoints.Get(j-1)(0) + lam * edgepoints.Get(j)(0); 1410 np(1) = (1-lam) * edgepoints.Get(j-1)(1) + lam * edgepoints.Get(j)(1); 1411 np(2) = (1-lam) * edgepoints.Get(j-1)(2) + lam * edgepoints.Get(j)(2); 1412 1413 thispi = PointIndex::INVALID; 1414 if (i == ne) 1415 { 1416 /* 1417 for (pi = PointIndex::BASE; 1418 pi < mesh.GetNP()+PointIndex::BASE; pi++) 1419 if (Dist(mesh[pi], np) < 1e-6) 1420 thispi = pi; 1421 */ 1422 1423 meshpoint_tree -> GetIntersecting (np-Vec<3> (di,di,di), 1424 np+Vec<3> (di,di,di), locsearch); 1425 if (locsearch.Size()) 1426 thispi = locsearch[0]; 1427 } 1428 1429 if (!thispi.IsValid()) 1430 { 1431 ProjectToEdge (surf1, surf2, np); 1432 thispi = mesh.AddPoint (np, layer, (i==ne) ? FIXEDPOINT : EDGEPOINT); 1433 1434 meshpoint_tree -> Insert (np, thispi); 1435 // (*testout) << "test2, store point " << thispi << ", p = " << np << endl; 1436 } 1437 1438 for (k = 1; k <= refedges.Size(); k++) 1439 { 1440 if (refedgesinv.Get(k)) 1441 { 1442 seg[0] = lastpi; 1443 seg[1] = thispi; 1444 } 1445 else 1446 { 1447 seg[0] = thispi; 1448 seg[1] = lastpi; 1449 } 1450 seg.si = refedges.Get(k).si; 1451 seg.domin = refedges.Get(k).domin; 1452 seg.domout = refedges.Get(k).domout; 1453 seg.tlosurf = refedges.Get(k).tlosurf; 1454 seg.edgenr = refedges.Get(k).edgenr; 1455 seg.surfnr1 = refedges.Get(k).surfnr1; 1456 seg.surfnr2 = refedges.Get(k).surfnr2; 1457 seg.seginfo = 0; 1458 if (k == 1) seg.seginfo = (refedgesinv.Get(k)) ? 2 : 1; 1459 mesh.AddSegment (seg); 1460 //(*testout) << "add seg " << mesh[seg.p1] << "-" << mesh[seg.p2] << endl; 1461 //(*testout) << "refedge " << k << " surf1 " << seg.surfnr1 << " surf2 " << seg.surfnr2 << " inv " << refedgesinv.Get(k) << endl; 1462 1463 double maxh = min2 (geometry.GetSurface(seg.surfnr1)->GetMaxH(), 1464 geometry.GetSurface(seg.surfnr2)->GetMaxH()); 1465 1466 if (seg.domin != -1) 1467 { 1468 const Solid * s1 = 1469 geometry.GetTopLevelObject(seg.domin) -> GetSolid(); 1470 maxh = min2 (maxh, s1->GetMaxH()); 1471 maxh = min2 (maxh, geometry.GetTopLevelObject(seg.domin)->GetMaxH()); 1472 mesh.RestrictLocalH (p, maxh); 1473 mesh.RestrictLocalH (np, maxh); 1474 } 1475 if (seg.domout != -1) 1476 { 1477 const Solid * s1 = 1478 geometry.GetTopLevelObject(seg.domout) -> GetSolid(); 1479 maxh = min2 (maxh, s1->GetMaxH()); 1480 maxh = min2 (maxh, geometry.GetTopLevelObject(seg.domout)->GetMaxH()); 1481 mesh.RestrictLocalH (p, maxh); 1482 mesh.RestrictLocalH (np, maxh); 1483 } 1484 if (seg.tlosurf != -1) 1485 { 1486 double hi = geometry.GetTopLevelObject(seg.tlosurf) -> GetMaxH(); 1487 maxh = min2 (maxh, hi); 1488 mesh.RestrictLocalH (p, maxh); 1489 mesh.RestrictLocalH (np, maxh); 1490 } 1491 } 1492 1493 p = np; 1494 lastpi = thispi; 1495 } 1496 1497 #ifdef DEVELOP 1498 (*testout) << " eplast = " << lastpi << " = " << p << endl; 1499 #endif 1500 } 1501 1502 1503 1504 1505 1506 1507 void EdgeCalculation :: StoreShortEdge(const NgArray<Segment> & refedges,const NgArray<bool> & refedgesinv,const NgArray<Point<3>> & edgepoints,const NgArray<double> & curvelength,int layer,Mesh & mesh)1508 StoreShortEdge (const NgArray<Segment> & refedges, 1509 const NgArray<bool> & refedgesinv, 1510 const NgArray<Point<3> > & edgepoints, 1511 const NgArray<double> & curvelength, 1512 int layer, 1513 Mesh & mesh) 1514 { 1515 1516 // Calculate optimal element-length 1517 PointIndex pi; 1518 // int ne; 1519 Segment seg; 1520 1521 /* 1522 double len, corr, lam; 1523 int thispi, lastpi; 1524 Point<3> p, np; 1525 1526 1527 const Surface * surf1 = geometry.GetSurface (refedges.Get(1).surfnr1); 1528 const Surface * surf2 = geometry.GetSurface (refedges.Get(1).surfnr2); 1529 1530 len = curvelength.Last(); 1531 ne = int (len + 0.5); 1532 if (ne == 0) ne = 1; 1533 if (Dist2 (edgepoints[1], edgepoints.Last()) < 1e-8 && 1534 ne <= 6) 1535 ne = 6; 1536 corr = len / ne; 1537 */ 1538 1539 // generate initial point 1540 Point<3> p = edgepoints[0]; 1541 PointIndex pi1 = PointIndex::INVALID; 1542 for (pi = PointIndex::BASE; 1543 pi < mesh.GetNP()+PointIndex::BASE; pi++) 1544 1545 if (Dist (mesh[pi], p) < 1e-6*geometry.MaxSize()) 1546 { 1547 pi1 = pi; 1548 break; 1549 } 1550 1551 if (!pi1.IsValid()) 1552 { 1553 pi1 = mesh.AddPoint (p, layer, FIXEDPOINT); 1554 meshpoint_tree -> Insert (p, pi1); 1555 // (*testout) << "test3, store point " << pi1 << ", p = " << p << endl; 1556 } 1557 1558 p = edgepoints.Last(); 1559 PointIndex pi2 = PointIndex::INVALID; 1560 for (pi = PointIndex::BASE; 1561 pi < mesh.GetNP()+PointIndex::BASE; pi++) 1562 1563 if (Dist (mesh[pi], p) < 1e-6*geometry.MaxSize()) 1564 { 1565 pi2 = pi; 1566 break; 1567 } 1568 if (!pi2.IsValid()) 1569 { 1570 pi2 = mesh.AddPoint (p, layer, FIXEDPOINT); 1571 meshpoint_tree -> Insert (p, pi2); 1572 // (*testout) << "test4, store point " << pi2 << ", p = " << p << endl; 1573 } 1574 1575 /* 1576 1577 j = 1; 1578 for (i = 1; i <= ne; i++) 1579 { 1580 while (curvelength[j] < i * corr && j < curvelength.Size()) j++; 1581 1582 lam = (i * corr - curvelength[j-1]) / 1583 (curvelength[j] - curvelength[j-1]); 1584 1585 np(0) = (1-lam) * edgepoints[j-1](0) + lam * edgepoints[j](0); 1586 np(1) = (1-lam) * edgepoints[j-1](1) + lam * edgepoints[j](1); 1587 np(2) = (1-lam) * edgepoints[j-1](2) + lam * edgepoints[j](2); 1588 1589 1590 thispi = 0; 1591 if (i == ne) 1592 for (j = 1; j <= mesh.GetNP(); j++) 1593 if (Dist(mesh.Point(j), np) < 1e-6) 1594 thispi = j; 1595 1596 if (!thispi) 1597 { 1598 ProjectToEdge (surf1, surf2, np); 1599 thispi = mesh.AddPoint (np); 1600 } 1601 */ 1602 1603 // (*testout) << "short edge " << pi1 << " - " << pi2 << endl; 1604 1605 for (int k = 1; k <= refedges.Size(); k++) 1606 { 1607 if (refedgesinv.Get(k)) 1608 { 1609 seg[0] = pi1; 1610 seg[1] = pi2; 1611 } 1612 else 1613 { 1614 seg[0] = pi2; 1615 seg[1] = pi1; 1616 } 1617 1618 seg.si = refedges.Get(k).si; 1619 seg.domin = refedges.Get(k).domin; 1620 seg.domout = refedges.Get(k).domout; 1621 seg.tlosurf = refedges.Get(k).tlosurf; 1622 seg.edgenr = refedges.Get(k).edgenr; 1623 seg.surfnr1 = refedges.Get(k).surfnr1; 1624 seg.surfnr2 = refedges.Get(k).surfnr2; 1625 seg.seginfo = 0; 1626 if (k == 1) seg.seginfo = (refedgesinv.Get(k)) ? 2 : 1; 1627 mesh.AddSegment (seg); 1628 // (*testout) << "add seg " << seg[0] << "-" << seg[1] << endl; 1629 } 1630 } 1631 1632 1633 1634 1635 1636 1637 1638 void EdgeCalculation :: CopyEdge(const NgArray<Segment> & refedges,const NgArray<bool> & refedgesinv,int copyfromedge,const Point<3> & fromstart,const Point<3> & fromend,const Point<3> & tostart,const Point<3> & toend,int copyedgeidentification,int layer,Mesh & mesh)1639 CopyEdge (const NgArray<Segment> & refedges, 1640 const NgArray<bool> & refedgesinv, 1641 int copyfromedge, 1642 const Point<3> & fromstart, const Point<3> & fromend, 1643 const Point<3> & tostart, const Point<3> & toend, 1644 int copyedgeidentification, 1645 int layer, 1646 Mesh & mesh) 1647 { 1648 int k; 1649 PointIndex pi; 1650 1651 double size = geometry.MaxSize(); 1652 1653 // copy start and end points 1654 for (int i = 1; i <= 2; i++) 1655 { 1656 Point<3> fromp = 1657 (i == 1) ? fromstart : fromend; 1658 Point<3> top = 1659 (i == 1) ? tostart : toend; 1660 1661 PointIndex frompi = PointIndex::INVALID; 1662 PointIndex topi = PointIndex::INVALID; 1663 for (pi = PointIndex::BASE; 1664 pi < mesh.GetNP()+PointIndex::BASE; pi++) 1665 { 1666 if (Dist2 (mesh[pi], fromp) <= 1e-16*size) 1667 frompi = pi; 1668 if (Dist2 (mesh[pi], top) <= 1e-16*size) 1669 topi = pi; 1670 } 1671 1672 1673 if (!topi.IsValid()) 1674 { 1675 topi = mesh.AddPoint (top, layer, FIXEDPOINT); 1676 meshpoint_tree -> Insert (top, topi); 1677 } 1678 1679 const Identification & csi = 1680 (*geometry.identifications.Get(copyedgeidentification)); 1681 1682 1683 if (csi.Identifyable (mesh[frompi], mesh[topi])) 1684 mesh.GetIdentifications().Add(frompi, topi, copyedgeidentification); 1685 else if (csi.Identifyable (mesh[topi], mesh[frompi])) 1686 mesh.GetIdentifications().Add(topi, frompi, copyedgeidentification); 1687 else 1688 { 1689 cerr << "edgeflw.cpp: should identify, but cannot"; 1690 exit(1); 1691 } 1692 #ifdef DEVELOP 1693 (*testout) << "adding identification " << mesh[frompi] << "; " << mesh[topi] 1694 << " (id " << copyedgeidentification <<")" << endl; 1695 #endif 1696 1697 1698 /* 1699 (*testout) << "Add Identification from CopyEdge, p1 = " 1700 << mesh[PointIndex(frompi)] << ", p2 = " 1701 << mesh[PointIndex(topi)] << endl; 1702 1703 mesh.GetIdentifications().Add(frompi, topi, copyedgeidentification); 1704 */ 1705 } 1706 1707 int oldns = mesh.GetNSeg(); 1708 for (int i = 1; i <= oldns; i++) 1709 { 1710 // real copy, since array might be reallocated !! 1711 const Segment oldseg = mesh.LineSegment(i); 1712 if (oldseg.edgenr != copyfromedge) 1713 continue; 1714 if (oldseg.seginfo == 0) 1715 continue; 1716 1717 int pi1 = oldseg[0]; 1718 int pi2 = oldseg[1]; 1719 1720 int npi1 = geometry.identifications.Get(copyedgeidentification) 1721 -> GetIdentifiedPoint (mesh, pi1); 1722 int npi2 = geometry.identifications.Get(copyedgeidentification) 1723 -> GetIdentifiedPoint (mesh, pi2); 1724 1725 //(*testout) << "copy edge, pts = " << npi1 << " - " << npi2 << endl; 1726 1727 Segment seg; 1728 1729 for (k = 1; k <= refedges.Size(); k++) 1730 { 1731 bool inv = refedgesinv.Get(k); 1732 1733 // other edge is inverse 1734 if (oldseg.seginfo == 1) 1735 inv = !inv; 1736 1737 // (*testout) << "inv, now = " << inv << endl; 1738 1739 if (inv) 1740 { 1741 seg[0] = npi1; 1742 seg[1] = npi2; 1743 } 1744 else 1745 { 1746 seg[0] = npi2; 1747 seg[1] = npi1; 1748 } 1749 seg.si = refedges.Get(k).si; 1750 seg.domin = refedges.Get(k).domin; 1751 seg.domout = refedges.Get(k).domout; 1752 seg.tlosurf = refedges.Get(k).tlosurf; 1753 seg.edgenr = refedges.Get(k).edgenr; 1754 seg.surfnr1 = refedges.Get(k).surfnr1; 1755 seg.surfnr2 = refedges.Get(k).surfnr2; 1756 seg.seginfo = 0; 1757 if (k == 1) seg.seginfo = refedgesinv.Get(k) ? 2 : 1; 1758 mesh.AddSegment (seg); 1759 // (*testout) << "copy seg " << seg[0] << "-" << seg[1] << endl; 1760 #ifdef DEVELOP 1761 1762 (*testout) << "copy seg, face = " << seg.si << ": " 1763 << " inv = " << inv << ", refinv = " << refedgesinv.Get(k) 1764 << mesh.Point(seg[0]) << ", " << mesh.Point(seg[1]) << endl; 1765 #endif 1766 1767 } 1768 1769 } 1770 } 1771 1772 1773 1774 1775 1776 1777 1778 void EdgeCalculation :: FindClosedSurfaces(double h,Mesh & mesh)1779 FindClosedSurfaces (double h, Mesh & mesh) 1780 { 1781 // if there is no special point at a sphere, one has to add a segment pair 1782 1783 int nsurf = geometry.GetNSurf(); 1784 int layer = 0; 1785 1786 // Solid * tansol; 1787 NgArray<int> tansurfind; 1788 1789 double size = geometry.MaxSize(); 1790 int nsol = geometry.GetNTopLevelObjects(); 1791 1792 1793 NgBitArray pointatsurface (nsurf); 1794 pointatsurface.Clear(); 1795 1796 for (int i = 1; i <= mesh.GetNSeg(); i++) 1797 { 1798 const Segment & seg = mesh.LineSegment(i); 1799 1800 #ifdef DEVELOP 1801 (*testout) << seg.surfnr1 << ", " << seg.surfnr2 << ", si = " << seg.si << endl; 1802 #endif 1803 int classrep = geometry.GetSurfaceClassRepresentant (seg.si); 1804 pointatsurface.Set (classrep); 1805 } 1806 1807 1808 for (int i = 0; i < nsurf; i++) 1809 { 1810 int classrep = geometry.GetSurfaceClassRepresentant (i); 1811 1812 if (!pointatsurface.Test(classrep)) 1813 { 1814 const Surface * s = geometry.GetSurface(i); 1815 Point<3> p1 = s -> GetSurfacePoint(); 1816 Vec<3> nv = s -> GetNormalVector (p1); 1817 1818 double hloc = 1819 min2 (s->LocH (p1, 3, 1, mparam, h), mesh.GetH(p1)); 1820 1821 1822 1823 Segment seg1; 1824 seg1.si = i; 1825 seg1.domin = -1; 1826 seg1.domout = -1; 1827 1828 Segment seg2; 1829 seg2.si = i; 1830 seg2.domin = -1; 1831 seg2.domout = -1; 1832 1833 seg1.surfnr1 = i; 1834 seg2.surfnr1 = i; 1835 seg1.surfnr2 = i; 1836 seg2.surfnr2 = i; 1837 1838 for (int j = 0; j < nsol; j++) 1839 { 1840 if (geometry.GetTopLevelObject(j)->GetSurface()) 1841 continue; 1842 1843 const Solid * sol = geometry.GetTopLevelObject(j)->GetSolid(); 1844 // sol -> TangentialSolid (p1, tansol, tansurfind, ideps*size); 1845 auto tansol = sol -> TangentialSolid (p1, tansurfind, ideps*size); 1846 layer = geometry.GetTopLevelObject(j)->GetLayer(); 1847 1848 1849 if (tansol) 1850 { 1851 tansol -> GetSurfaceIndices (tansurfind); 1852 1853 if (tansurfind.Size() == 1 && tansurfind.Get(1) == i) 1854 { 1855 hloc = min2 (hloc, geometry.GetTopLevelObject(j)->GetMaxH()); 1856 1857 if (!tansol->VectorIn(p1, nv)) 1858 { 1859 seg1.domin = j; 1860 seg2.domin = j; 1861 seg1.tlosurf = -1; 1862 seg2.tlosurf = -1; 1863 } 1864 else 1865 { 1866 seg1.domout = j; 1867 seg2.domout = j; 1868 seg1.tlosurf = -1; 1869 seg2.tlosurf = -1; 1870 } 1871 // seg.s2 = i; 1872 // seg.invs1 = surfaces[i] -> Inverse(); 1873 // seg.invs2 = ! (surfaces[i] -> Inverse()); 1874 } 1875 // delete tansol; 1876 } 1877 } 1878 1879 1880 Vec<3> tv = nv.GetNormal (); 1881 tv *= (hloc / tv.Length()); 1882 Point<3> p2 = p1 + tv; 1883 s->Project (p2); 1884 1885 1886 if (seg1.domin != -1 || seg1.domout != -1) 1887 { 1888 mesh.AddPoint (p1, layer, EDGEPOINT); 1889 mesh.AddPoint (p2, layer, EDGEPOINT); 1890 seg1[0] = mesh.GetNP()-1; 1891 seg1[1] = mesh.GetNP(); 1892 seg2[1] = mesh.GetNP()-1; 1893 seg2[0] = mesh.GetNP(); 1894 seg1.geominfo[0].trignum = 1; 1895 seg1.geominfo[1].trignum = 1; 1896 seg2.geominfo[0].trignum = 1; 1897 seg2.geominfo[1].trignum = 1; 1898 mesh.AddSegment (seg1); 1899 mesh.AddSegment (seg2); 1900 1901 PrintMessage (5, "Add line segment to smooth surface"); 1902 1903 #ifdef DEVELOP 1904 (*testout) << "Add segment at smooth surface " << i; 1905 if (i != classrep) (*testout) << ", classrep = " << classrep; 1906 (*testout) << ": " 1907 << mesh.Point (mesh.GetNP()-1) << " - " 1908 << mesh.Point (mesh.GetNP()) << endl; 1909 #endif 1910 } 1911 } 1912 } 1913 } 1914 1915 } 1916