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