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