1 #ifdef OCCGEOMETRY 2 3 #include <mystdlib.h> 4 #include <occgeom.hpp> 5 #include <meshing.hpp> 6 7 8 namespace netgen 9 { 10 11 #include "occmeshsurf.hpp" 12 13 #define TCL_OK 0 14 #define TCL_ERROR 1 15 16 #define DIVIDEEDGESECTIONS 10000 // better solution to come soon 17 #define IGNORECURVELENGTH 1e-4 18 #define VSMALL 1e-10 19 20 21 DLL_HEADER bool merge_solids = false; 22 23 24 // can you please explain what you intend to compute here (JS) !!! Dist(Line l)25 double Line :: Dist (Line l) 26 { 27 Vec<3> n = p1-p0; 28 Vec<3> q = l.p1-l.p0; 29 double nq = n*q; 30 31 Point<3> p = p0 + 0.5*n; 32 double lambda = (p-l.p0)*n / (nq + VSMALL); 33 34 if (lambda >= 0 && lambda <= 1) 35 { 36 double d = (p-l.p0-lambda*q).Length(); 37 // if (d < 1e-3) d = 1e99; 38 return d; 39 } 40 else 41 return 1e99; 42 } 43 44 ComputeH(double kappa,const MeshingParameters & mparam)45 double ComputeH (double kappa, const MeshingParameters & mparam) 46 { 47 kappa *= mparam.curvaturesafety; 48 /* 49 double hret; 50 51 if (mparam.maxh * kappa < 1) 52 hret = mparam.maxh; 53 else 54 hret = 1 / (kappa + VSMALL); 55 56 if (mparam.maxh < hret) 57 hret = mparam.maxh; 58 59 return hret; 60 */ 61 // return min(mparam.maxh, 1/kappa); 62 return (mparam.maxh*kappa < 1) ? mparam.maxh : 1/kappa; 63 } 64 65 66 RestrictHTriangle(gp_Pnt2d & par0,gp_Pnt2d & par1,gp_Pnt2d & par2,BRepLProp_SLProps * prop,BRepLProp_SLProps * prop2,Mesh & mesh,int depth,double h,const MeshingParameters & mparam)67 void RestrictHTriangle (gp_Pnt2d & par0, gp_Pnt2d & par1, gp_Pnt2d & par2, 68 BRepLProp_SLProps * prop, BRepLProp_SLProps * prop2, Mesh & mesh, int depth, double h, 69 const MeshingParameters & mparam) 70 { 71 int ls = -1; 72 73 gp_Pnt pnt0,pnt1,pnt2; 74 75 prop->SetParameters (par0.X(), par0.Y()); 76 pnt0 = prop->Value(); 77 78 prop->SetParameters (par1.X(), par1.Y()); 79 pnt1 = prop->Value(); 80 81 prop->SetParameters (par2.X(), par2.Y()); 82 pnt2 = prop->Value(); 83 84 double aux; 85 double maxside = pnt0.Distance(pnt1); 86 ls = 2; 87 aux = pnt1.Distance(pnt2); 88 if(aux > maxside) 89 { 90 maxside = aux; 91 ls = 0; 92 } 93 aux = pnt2.Distance(pnt0); 94 if(aux > maxside) 95 { 96 maxside = aux; 97 ls = 1; 98 } 99 100 101 102 gp_Pnt2d parmid; 103 104 parmid.SetX( (par0.X()+par1.X()+par2.X()) / 3 ); 105 parmid.SetY( (par0.Y()+par1.Y()+par2.Y()) / 3 ); 106 107 if (depth%3 == 0) 108 { 109 double curvature = 0; 110 111 prop2->SetParameters (parmid.X(), parmid.Y()); 112 if (!prop2->IsCurvatureDefined()) 113 { 114 (*testout) << "curvature not defined!" << endl; 115 return; 116 } 117 curvature = max(fabs(prop2->MinCurvature()), 118 fabs(prop2->MaxCurvature())); 119 120 prop2->SetParameters (par0.X(), par0.Y()); 121 if (!prop2->IsCurvatureDefined()) 122 { 123 (*testout) << "curvature not defined!" << endl; 124 return; 125 } 126 curvature = max(curvature,max(fabs(prop2->MinCurvature()), 127 fabs(prop2->MaxCurvature()))); 128 129 prop2->SetParameters (par1.X(), par1.Y()); 130 if (!prop2->IsCurvatureDefined()) 131 { 132 (*testout) << "curvature not defined!" << endl; 133 return; 134 } 135 curvature = max(curvature,max(fabs(prop2->MinCurvature()), 136 fabs(prop2->MaxCurvature()))); 137 138 prop2->SetParameters (par2.X(), par2.Y()); 139 if (!prop2->IsCurvatureDefined()) 140 { 141 (*testout) << "curvature not defined!" << endl; 142 return; 143 } 144 curvature = max(curvature,max(fabs(prop2->MinCurvature()), 145 fabs(prop2->MaxCurvature()))); 146 147 //(*testout) << "curvature " << curvature << endl; 148 149 if (curvature < 1e-3) 150 { 151 //(*testout) << "curvature too small (" << curvature << ")!" << endl; 152 return; 153 // return war bis 10.2.05 auskommentiert 154 } 155 156 157 158 h = ComputeH (curvature+1e-10, mparam); 159 160 if(h < 1e-4*maxside) 161 return; 162 163 164 // if (h > 30) return; 165 } 166 167 if (h < maxside && depth < 10) 168 { 169 //cout << "\r h " << h << flush; 170 gp_Pnt2d pm; 171 172 //cout << "h " << h << " maxside " << maxside << " depth " << depth << endl; 173 //cout << "par0 " << par0.X() << " " << par0.Y() 174 //<< " par1 " << par1.X() << " " << par1.Y() 175 // << " par2 " << par2.X() << " " << par2.Y()<< endl; 176 177 if(ls == 0) 178 { 179 pm.SetX(0.5*(par1.X()+par2.X())); pm.SetY(0.5*(par1.Y()+par2.Y())); 180 RestrictHTriangle(pm, par2, par0, prop, prop2, mesh, depth+1, h, mparam); 181 RestrictHTriangle(pm, par0, par1, prop, prop2, mesh, depth+1, h, mparam); 182 } 183 else if(ls == 1) 184 { 185 pm.SetX(0.5*(par0.X()+par2.X())); pm.SetY(0.5*(par0.Y()+par2.Y())); 186 RestrictHTriangle(pm, par1, par2, prop, prop2, mesh, depth+1, h, mparam); 187 RestrictHTriangle(pm, par0, par1, prop, prop2, mesh, depth+1, h, mparam); 188 } 189 else if(ls == 2) 190 { 191 pm.SetX(0.5*(par0.X()+par1.X())); pm.SetY(0.5*(par0.Y()+par1.Y())); 192 RestrictHTriangle(pm, par1, par2, prop, prop2, mesh, depth+1, h, mparam); 193 RestrictHTriangle(pm, par2, par0, prop, prop2, mesh, depth+1, h, mparam); 194 } 195 196 } 197 else 198 { 199 gp_Pnt pnt; 200 Point3d p3d; 201 202 prop->SetParameters (parmid.X(), parmid.Y()); 203 pnt = prop->Value(); 204 p3d = Point3d(pnt.X(), pnt.Y(), pnt.Z()); 205 mesh.RestrictLocalH (p3d, h); 206 207 p3d = Point3d(pnt0.X(), pnt0.Y(), pnt0.Z()); 208 mesh.RestrictLocalH (p3d, h); 209 210 p3d = Point3d(pnt1.X(), pnt1.Y(), pnt1.Z()); 211 mesh.RestrictLocalH (p3d, h); 212 213 p3d = Point3d(pnt2.X(), pnt2.Y(), pnt2.Z()); 214 mesh.RestrictLocalH (p3d, h); 215 216 //(*testout) << "p = " << p3d << ", h = " << h << ", maxside = " << maxside << endl; 217 218 } 219 } 220 221 222 DivideEdge(TopoDS_Edge & edge,NgArray<MeshPoint> & ps,Array<double> & params,Mesh & mesh,const MeshingParameters & mparam)223 void DivideEdge (TopoDS_Edge & edge, NgArray<MeshPoint> & ps, 224 Array<double> & params, Mesh & mesh, 225 const MeshingParameters & mparam) 226 { 227 double s0, s1; 228 int nsubedges = 1; 229 gp_Pnt pnt, oldpnt; 230 double svalue[DIVIDEEDGESECTIONS]; 231 232 GProp_GProps system; 233 BRepGProp::LinearProperties(edge, system); 234 double L = system.Mass(); 235 236 Handle(Geom_Curve) c = BRep_Tool::Curve(edge, s0, s1); 237 238 double hvalue[DIVIDEEDGESECTIONS+1]; 239 hvalue[0] = 0; 240 pnt = c->Value(s0); 241 242 int tmpVal = (int)(DIVIDEEDGESECTIONS); 243 244 for (int i = 1; i <= tmpVal; i++) 245 { 246 oldpnt = pnt; 247 pnt = c->Value(s0+(i/double(DIVIDEEDGESECTIONS))*(s1-s0)); 248 hvalue[i] = hvalue[i-1] + 249 1.0/mesh.GetH(Point3d(pnt.X(), pnt.Y(), pnt.Z()))* 250 pnt.Distance(oldpnt); 251 252 //(*testout) << "mesh.GetH(Point3d(pnt.X(), pnt.Y(), pnt.Z())) " << mesh.GetH(Point3d(pnt.X(), pnt.Y(), pnt.Z())) 253 // << " pnt.Distance(oldpnt) " << pnt.Distance(oldpnt) << endl; 254 } 255 256 // nsubedges = int(ceil(hvalue[DIVIDEEDGESECTIONS])); 257 nsubedges = max (1, int(floor(hvalue[DIVIDEEDGESECTIONS]+0.5))); 258 259 ps.SetSize(nsubedges-1); 260 params.SetSize(nsubedges+1); 261 262 int i = 1; 263 int i1 = 0; 264 do 265 { 266 if (hvalue[i1]/hvalue[DIVIDEEDGESECTIONS]*nsubedges >= i) 267 { 268 params[i] = s0+(i1/double(DIVIDEEDGESECTIONS))*(s1-s0); 269 pnt = c->Value(params[i]); 270 ps[i-1] = MeshPoint (Point3d(pnt.X(), pnt.Y(), pnt.Z())); 271 i++; 272 } 273 i1++; 274 if (i1 > DIVIDEEDGESECTIONS) 275 { 276 nsubedges = i; 277 ps.SetSize(nsubedges-1); 278 params.SetSize(nsubedges+1); 279 cout << "divide edge: local h too small" << endl; 280 } 281 } while (i < nsubedges); 282 283 params[0] = s0; 284 params[nsubedges] = s1; 285 286 if (params[nsubedges] <= params[nsubedges-1]) 287 { 288 cout << "CORRECTED" << endl; 289 ps.SetSize (nsubedges-2); 290 params.SetSize (nsubedges); 291 params[nsubedges] = s1; 292 } 293 } 294 295 296 297 OCCFindEdges(const OCCGeometry & geom,Mesh & mesh,const MeshingParameters & mparam)298 void OCCFindEdges (const OCCGeometry & geom, Mesh & mesh, const MeshingParameters & mparam) 299 { 300 static Timer t("OCCFindEdges"); RegionTimer r(t); 301 static Timer tsearch("OCCFindEdges - search point"); 302 const char * savetask = multithread.task; 303 multithread.task = "Edge meshing"; 304 305 (*testout) << "edge meshing" << endl; 306 307 int nvertices = geom.vmap.Extent(); 308 int nedges = geom.emap.Extent(); 309 310 Array<Array<PointIndex>> alledgepnums(nedges); 311 Array<Array<double>> alledgeparams(nedges); 312 313 (*testout) << "nvertices = " << nvertices << endl; 314 (*testout) << "nedges = " << nedges << endl; 315 316 double eps = 1e-6 * geom.GetBoundingBox().Diam(); 317 318 tsearch.Start(); 319 for (int i = 1; i <= nvertices; i++) 320 { 321 gp_Pnt pnt = BRep_Tool::Pnt (TopoDS::Vertex(geom.vmap(i))); 322 MeshPoint mp(occ2ng(pnt)); 323 324 bool exists = false; 325 if (merge_solids) 326 for (PointIndex pi : mesh.Points().Range()) 327 if (Dist2 (mesh[pi], Point<3>(mp)) < eps*eps) 328 { 329 exists = true; 330 break; 331 } 332 333 if (!exists) 334 mesh.AddPoint (mp); 335 336 double maxh = OCCGeometry::global_shape_properties[TopoDS::Vertex(geom.vmap(i)).TShape()].maxh; 337 mesh.RestrictLocalH (occ2ng(pnt), maxh); 338 } 339 tsearch.Stop(); 340 341 (*testout) << "different vertices = " << mesh.GetNP() << endl; 342 343 // int first_ep = mesh.GetNP()+1; 344 // PointIndex first_ep = mesh.Points().End(); 345 PointIndex first_ep = *mesh.Points().Range().end(); 346 auto vertexrange = mesh.Points().Range(); 347 348 NgArray<int> face2solid[2]; 349 for (int i = 0; i < 2; i++) 350 { 351 face2solid[i].SetSize (geom.fmap.Extent()); 352 face2solid[i] = 0; 353 } 354 355 int solidnr = 0; 356 for (TopExp_Explorer exp0(geom.shape, TopAbs_SOLID); exp0.More(); exp0.Next()) 357 { 358 solidnr++; 359 for (TopExp_Explorer exp1(exp0.Current(), TopAbs_FACE); exp1.More(); exp1.Next()) 360 { 361 TopoDS_Face face = TopoDS::Face(exp1.Current()); 362 int facenr = geom.fmap.FindIndex(face); 363 if(facenr < 1) continue; 364 365 if (face2solid[0][facenr-1] == 0) 366 face2solid[0][facenr-1] = solidnr; 367 else 368 face2solid[1][facenr-1] = solidnr; 369 } 370 } 371 372 373 int total = 0; 374 for (int i3 = 1; i3 <= geom.fmap.Extent(); i3++) 375 for (TopExp_Explorer exp2(geom.fmap(i3), TopAbs_WIRE); exp2.More(); exp2.Next()) 376 for (TopExp_Explorer exp3(exp2.Current(), TopAbs_EDGE); exp3.More(); exp3.Next()) 377 total++; 378 379 380 int facenr = 0; 381 // int edgenr = mesh.GetNSeg(); 382 383 (*testout) << "faces = " << geom.fmap.Extent() << endl; 384 int curr = 0; 385 386 for (int i3 = 1; i3 <= geom.fmap.Extent(); i3++) 387 { 388 TopoDS_Face face = TopoDS::Face(geom.fmap(i3)); 389 facenr = geom.fmap.FindIndex (face); // sollte doch immer == i3 sein ??? JS 390 391 int solidnr0 = face2solid[0][i3-1]; 392 int solidnr1 = face2solid[1][i3-1]; 393 394 /* auskommentiert am 3.3.05 von robert 395 for (exp2.Init (geom.somap(solidnr0), TopAbs_FACE); exp2.More(); exp2.Next()) 396 { 397 TopoDS_Face face2 = TopoDS::Face(exp2.Current()); 398 if (geom.fmap.FindIndex(face2) == facenr) 399 { 400 // if (face.Orientation() != face2.Orientation()) swap (solidnr0, solidnr1); 401 } 402 } 403 */ 404 405 mesh.AddFaceDescriptor (FaceDescriptor(facenr, solidnr0, solidnr1, 0)); 406 407 // Philippose - 06/07/2009 408 // Add the face colour to the mesh data 409 Quantity_Color face_colour; 410 411 if(!(geom.face_colours.IsNull()) 412 && (geom.face_colours->GetColor(face,XCAFDoc_ColorSurf,face_colour))) 413 { 414 mesh.GetFaceDescriptor(facenr).SetSurfColour({face_colour.Red(),face_colour.Green(),face_colour.Blue()}); 415 } 416 else 417 { 418 auto it = OCCGeometry::global_shape_properties.find(face.TShape()); 419 if (it != OCCGeometry::global_shape_properties.end() && it->second.col) 420 { 421 Vec<3> col = it->second.col.value_or(Vec<3>(0,1,0)); 422 mesh.GetFaceDescriptor(facenr).SetSurfColour(col); 423 } 424 else 425 mesh.GetFaceDescriptor(facenr).SetSurfColour({0.0,1.0,0.0}); 426 } 427 428 if(geom.fnames.Size()>=facenr) 429 mesh.GetFaceDescriptor(facenr).SetBCName(&geom.fnames[facenr-1]); 430 mesh.GetFaceDescriptor(facenr).SetBCProperty(facenr); 431 // ACHTUNG! STIMMT NICHT ALLGEMEIN (RG) 432 433 434 Handle(Geom_Surface) occface = BRep_Tool::Surface(face); 435 436 for (TopExp_Explorer exp2 (face, TopAbs_WIRE); exp2.More(); exp2.Next()) 437 { 438 TopoDS_Shape wire = exp2.Current(); 439 440 for (TopExp_Explorer exp3 (wire, TopAbs_EDGE); exp3.More(); exp3.Next()) 441 { 442 curr++; 443 (*testout) << "edge nr " << curr << endl; 444 multithread.percent = 100 * curr / double (total); 445 if (multithread.terminate) return; 446 447 TopoDS_Edge edge = TopoDS::Edge (exp3.Current()); 448 if (BRep_Tool::Degenerated(edge)) 449 { 450 //(*testout) << "ignoring degenerated edge" << endl; 451 continue; 452 } 453 454 if(geom.emap.FindIndex(edge) < 1) continue; 455 456 if (geom.vmap.FindIndex(TopExp::FirstVertex (edge)) == 457 geom.vmap.FindIndex(TopExp::LastVertex (edge))) 458 { 459 GProp_GProps system; 460 BRepGProp::LinearProperties(edge, system); 461 462 if (system.Mass() < eps) 463 { 464 cout << "ignoring edge " << geom.emap.FindIndex (edge) 465 << ". closed edge with length < " << eps << endl; 466 continue; 467 } 468 } 469 470 471 Handle(Geom2d_Curve) cof; 472 double s0, s1; 473 cof = BRep_Tool::CurveOnSurface (edge, face, s0, s1); 474 475 int geomedgenr = geom.emap.FindIndex(edge); 476 Array<PointIndex> pnums; 477 Array<double> params; 478 479 480 // check for identifications 481 bool copy_identified = false; 482 if (auto it = geom.identifications.find(edge.TShape()); it != geom.identifications.end()) 483 for (auto & ids : it->second) 484 { 485 cout << "edge has identification with trafo " << ids.name << ", inv = " << ids.inverse << endl; 486 int otherind = geom.emap.FindIndex(ids.other); 487 Array<Segment> othersegs; 488 for (auto & seg : mesh.LineSegments()) 489 if (seg.edgenr == otherind) 490 othersegs.Append (seg); 491 492 if (othersegs.Size()) 493 { 494 cout << "other has already segs" << endl; 495 copy_identified = true; 496 497 Array<PointIndex> pnums_other; 498 pnums_other.Append (othersegs[0][0]); 499 for (auto & seg : othersegs) 500 pnums_other.Append (seg[1]); 501 502 auto inv = ids.trafo.CalcInverse(); 503 // for (auto & pi : pnums) 504 for (auto oi : Range(pnums_other)) 505 { 506 PointIndex piother = pnums_other[pnums_other.Size()-oi-1]; 507 Point<3> pother = mesh[piother]; 508 Point<3> p = inv(pother); 509 510 bool found = false; 511 PointIndex pi; 512 for (PointIndex piv : vertexrange) 513 if (Dist2 (mesh[piv], p) < eps*eps) 514 { 515 pi = piv; 516 found = true; 517 } 518 519 if (!found) 520 pi = mesh.AddPoint (p); 521 522 // params.Add ( find parameter p ); 523 double s0, s1; 524 Handle(Geom_Curve) curve = BRep_Tool::Curve(edge, s0, s1); 525 526 GeomAPI_ProjectPointOnCurve proj(ng2occ(p), curve); 527 params.Append (proj.LowerDistanceParameter()); 528 pnums.Append (pi); 529 mesh.GetIdentifications().Add (pi, piother, geomedgenr); 530 531 } 532 mesh.GetIdentifications().SetType(geomedgenr,Identifications::PERIODIC); 533 534 copy_identified = true; 535 break; 536 } 537 } 538 539 540 if (!copy_identified) 541 { 542 543 544 if (alledgepnums[geomedgenr-1].Size()) 545 { 546 pnums = alledgepnums[geomedgenr-1]; 547 params = alledgeparams[geomedgenr-1]; 548 } 549 else 550 { 551 NgArray <MeshPoint> mp; 552 DivideEdge (edge, mp, params, mesh, mparam); 553 554 pnums.SetSize(mp.Size()+2); 555 if (!merge_solids) 556 { 557 pnums[0] = geom.vmap.FindIndex (TopExp::FirstVertex (edge)) + PointIndex::BASE-1; 558 pnums.Last() = geom.vmap.FindIndex (TopExp::LastVertex (edge)) + PointIndex::BASE-1; 559 } 560 else 561 { 562 Point<3> fp = occ2ng (BRep_Tool::Pnt (TopExp::FirstVertex (edge))); 563 Point<3> lp = occ2ng (BRep_Tool::Pnt (TopExp::LastVertex (edge))); 564 565 pnums[0] = PointIndex::INVALID; 566 pnums.Last() = PointIndex::INVALID; 567 for (PointIndex pi : vertexrange) 568 { 569 if (Dist2 (mesh[pi], fp) < eps*eps) pnums[0] = pi; 570 if (Dist2 (mesh[pi], lp) < eps*eps) pnums.Last() = pi; 571 } 572 } 573 574 for (size_t i = 1; i <= mp.Size(); i++) 575 { 576 bool exists = false; 577 tsearch.Start(); 578 579 /* 580 // for (PointIndex j = first_ep; j < mesh.Points().End(); j++) 581 for (PointIndex j = first_ep; j < *mesh.Points().Range().end(); j++) 582 if ((mesh.Point(j)-Point<3>(mp[i-1])).Length() < eps) 583 { 584 exists = true; 585 pnums[i] = j; 586 break; 587 } 588 */ 589 tsearch.Stop(); 590 591 if (!exists) 592 pnums[i] = mesh.AddPoint (mp[i-1]); 593 } 594 } 595 596 597 alledgepnums[geomedgenr-1] = pnums; 598 alledgeparams[geomedgenr-1] = params; 599 } 600 601 auto name = geom.enames.Size() ? geom.enames[geom.emap.FindIndex(edge)-1] : ""; 602 if(name != "") 603 mesh.SetCD2Name(geomedgenr, name); 604 605 (*testout) << "NP = " << mesh.GetNP() << endl; 606 //(*testout) << pnums[pnums.Size()-1] << endl; 607 608 // for (size_t i = 1; i <= mp.Size()+1; i++) 609 for (size_t i = 1; i < pnums.Size(); i++) 610 { 611 // edgenr++; 612 Segment seg; 613 614 seg[0] = pnums[i-1]; 615 seg[1] = pnums[i]; 616 // seg.edgenr = edgenr; 617 seg.edgenr = geomedgenr; 618 seg.si = facenr; 619 seg.epgeominfo[0].dist = params[i-1]; 620 seg.epgeominfo[1].dist = params[i]; 621 seg.epgeominfo[0].edgenr = geomedgenr; 622 seg.epgeominfo[1].edgenr = geomedgenr; 623 624 double s0 = params[i-1]; 625 double s1 = params[i]; 626 double delta = s1-s0; 627 s0 += 1e-10*delta; // fixes normal-vector roundoff problem when endpoint is cone-tip 628 s1 -= 1e-10*delta; 629 gp_Pnt2d p2d1 = cof->Value(s0); 630 gp_Pnt2d p2d2 = cof->Value(s1); 631 seg.epgeominfo[0].u = p2d1.X(); 632 seg.epgeominfo[0].v = p2d1.Y(); 633 seg.epgeominfo[1].u = p2d2.X(); 634 seg.epgeominfo[1].v = p2d2.Y(); 635 636 /* 637 if (occface->IsUPeriodic()) 638 { 639 cout << "U Periodic" << endl; 640 if (fabs(seg.epgeominfo[1].u-seg.epgeominfo[0].u) > 641 fabs(seg.epgeominfo[1].u- 642 (seg.epgeominfo[0].u-occface->UPeriod()))) 643 seg.epgeominfo[0].u = p2d.X()+occface->UPeriod(); 644 645 if (fabs(seg.epgeominfo[1].u-seg.epgeominfo[0].u) > 646 fabs(seg.epgeominfo[1].u- 647 (seg.epgeominfo[0].u+occface->UPeriod()))) 648 seg.epgeominfo[0].u = p2d.X()-occface->UPeriod(); 649 } 650 651 if (occface->IsVPeriodic()) 652 { 653 cout << "V Periodic" << endl; 654 if (fabs(seg.epgeominfo[1].v-seg.epgeominfo[0].v) > 655 fabs(seg.epgeominfo[1].v- 656 (seg.epgeominfo[0].v-occface->VPeriod()))) 657 seg.epgeominfo[0].v = p2d.Y()+occface->VPeriod(); 658 659 if (fabs(seg.epgeominfo[1].v-seg.epgeominfo[0].v) > 660 fabs(seg.epgeominfo[1].v- 661 (seg.epgeominfo[0].v+occface->VPeriod()))) 662 seg.epgeominfo[0].v = p2d.Y()-occface->VPeriod(); 663 } 664 */ 665 666 if (edge.Orientation() == TopAbs_REVERSED) 667 { 668 swap (seg[0], seg[1]); 669 swap (seg.epgeominfo[0].dist, seg.epgeominfo[1].dist); 670 swap (seg.epgeominfo[0].u, seg.epgeominfo[1].u); 671 swap (seg.epgeominfo[0].v, seg.epgeominfo[1].v); 672 } 673 674 mesh.AddSegment (seg); 675 676 //edgesegments[geomedgenr-1]->Append(mesh.GetNSeg()); 677 678 } 679 } 680 } 681 } 682 683 // for(i=1; i<=mesh.GetNSeg(); i++) 684 // (*testout) << "edge " << mesh.LineSegment(i).edgenr << " face " << mesh.LineSegment(i).si 685 // << " p1 " << mesh.LineSegment(i)[0] << " p2 " << mesh.LineSegment(i)[1] << endl; 686 // exit(10); 687 688 mesh.CalcSurfacesOfNode(); 689 multithread.task = savetask; 690 } 691 692 693 694 OCCMeshSurface(const OCCGeometry & geom,Mesh & mesh,const MeshingParameters & mparam)695 void OCCMeshSurface (const OCCGeometry & geom, Mesh & mesh, 696 const MeshingParameters & mparam) 697 { 698 static Timer t("OCCMeshSurface"); RegionTimer r(t); 699 700 // int i, j, k; 701 // int changed; 702 703 const char * savetask = multithread.task; 704 multithread.task = "Surface meshing"; 705 706 geom.facemeshstatus = 0; 707 708 int noldp = mesh.GetNP(); 709 710 double starttime = GetTime(); 711 712 NgArray<int, PointIndex::BASE> glob2loc(noldp); 713 714 //int projecttype = PARAMETERSPACE; 715 716 int projecttype = PARAMETERSPACE; 717 718 int notrys = 1; 719 720 int surfmesherror = 0; 721 722 for (int k = 1; k <= mesh.GetNFD(); k++) 723 { 724 if(1==0 && !geom.fvispar[k-1].IsDrawable()) 725 { 726 (*testout) << "ignoring face " << k << endl; 727 cout << "ignoring face " << k << endl; 728 continue; 729 } 730 731 (*testout) << "mesh face " << k << endl; 732 multithread.percent = 100 * k / (mesh.GetNFD() + VSMALL); 733 geom.facemeshstatus[k-1] = -1; 734 735 FaceDescriptor & fd = mesh.GetFaceDescriptor(k); 736 737 int oldnf = mesh.GetNSE(); 738 739 Box<3> bb = geom.GetBoundingBox(); 740 741 // int projecttype = PLANESPACE; 742 // int projecttype = PARAMETERSPACE; 743 744 static Timer tinit("init"); 745 tinit.Start(); 746 Meshing2OCCSurfaces meshing(geom, TopoDS::Face(geom.fmap(k)), bb, projecttype, mparam); 747 tinit.Stop(); 748 749 750 static Timer tprint("print"); 751 tprint.Start(); 752 if (meshing.GetProjectionType() == PLANESPACE) 753 PrintMessage (2, "Face ", k, " / ", mesh.GetNFD(), " (plane space projection)"); 754 else 755 PrintMessage (2, "Face ", k, " / ", mesh.GetNFD(), " (parameter space projection)"); 756 tprint.Stop(); 757 if (surfmesherror) 758 cout << "Surface meshing error occurred before (in " << surfmesherror << " faces)" << endl; 759 760 // Meshing2OCCSurfaces meshing(f2, bb); 761 meshing.SetStartTime (starttime); 762 //(*testout) << "Face " << k << endl << endl; 763 764 765 if (meshing.GetProjectionType() == PLANESPACE) 766 { 767 static Timer t("MeshSurface: Find edges and points - Physical"); RegionTimer r(t); 768 int cntp = 0; 769 glob2loc = 0; 770 771 for (Segment & seg : mesh.LineSegments()) 772 if (seg.si == k) 773 for (int j = 0; j < 2; j++) 774 { 775 PointIndex pi = seg[j]; 776 if (glob2loc[pi] == 0) 777 { 778 meshing.AddPoint (mesh.Point(pi), pi); 779 cntp++; 780 glob2loc[pi] = cntp; 781 } 782 } 783 784 /* 785 for (int i = 1; i <= mesh.GetNSeg(); i++) 786 { 787 Segment & seg = mesh.LineSegment(i); 788 */ 789 for (Segment & seg : mesh.LineSegments()) 790 if (seg.si == k) 791 { 792 PointGeomInfo gi0, gi1; 793 gi0.trignum = gi1.trignum = k; 794 gi0.u = seg.epgeominfo[0].u; 795 gi0.v = seg.epgeominfo[0].v; 796 gi1.u = seg.epgeominfo[1].u; 797 gi1.v = seg.epgeominfo[1].v; 798 799 meshing.AddBoundaryElement (glob2loc[seg[0]], glob2loc[seg[1]], gi0, gi1); 800 } 801 } 802 else 803 { 804 static Timer t("MeshSurface: Find edges and points - Parameter"); RegionTimer r(t); 805 806 int cntp = 0; 807 808 /* 809 for (int i = 1; i <= mesh.GetNSeg(); i++) 810 if (mesh.LineSegment(i).si == k) 811 cntp+=2; 812 */ 813 for (Segment & seg : mesh.LineSegments()) 814 if (seg.si == k) 815 cntp += 2; 816 817 NgArray<PointGeomInfo> gis; 818 819 gis.SetAllocSize (cntp); 820 gis.SetSize (0); 821 822 for (int i = 1; i <= mesh.GetNSeg(); i++) 823 { 824 Segment & seg = mesh.LineSegment(i); 825 if (seg.si == k) 826 { 827 PointGeomInfo gi0, gi1; 828 gi0.trignum = gi1.trignum = k; 829 gi0.u = seg.epgeominfo[0].u; 830 gi0.v = seg.epgeominfo[0].v; 831 gi1.u = seg.epgeominfo[1].u; 832 gi1.v = seg.epgeominfo[1].v; 833 834 int locpnum[2] = {0, 0}; 835 836 for (int j = 0; j < 2; j++) 837 { 838 PointGeomInfo gi = (j == 0) ? gi0 : gi1; 839 840 /* 841 int l; 842 for (l = 0; l < gis.Size() && locpnum[j] == 0; l++) 843 { 844 double dist = sqr (gis[l].u-gi.u)+sqr(gis[l].v-gi.v); 845 846 if (dist < 1e-10) 847 locpnum[j] = l+1; 848 } 849 850 if (locpnum[j] == 0) 851 { 852 PointIndex pi = seg[j]; 853 meshing.AddPoint (mesh.Point(pi), pi); 854 855 gis.SetSize (gis.Size()+1); 856 gis[l] = gi; 857 locpnum[j] = l+1; 858 } 859 */ 860 for (int l = 0; l < gis.Size(); l++) 861 { 862 double dist = sqr (gis[l].u-gi.u)+sqr(gis[l].v-gi.v); 863 if (dist < 1e-10) 864 { 865 locpnum[j] = l+1; 866 break; 867 } 868 } 869 870 if (locpnum[j] == 0) 871 { 872 PointIndex pi = seg[j]; 873 meshing.AddPoint (mesh.Point(pi), pi); 874 875 gis.Append (gi); 876 locpnum[j] = gis.Size(); 877 } 878 } 879 880 meshing.AddBoundaryElement (locpnum[0], locpnum[1], gi0, gi1); 881 } 882 } 883 } 884 885 886 887 888 // Philippose - 15/01/2009 889 double maxh = min2(geom.face_maxh[k-1], OCCGeometry::global_shape_properties[TopoDS::Face(geom.fmap(k)).TShape()].maxh); 890 //double maxh = mparam.maxh; 891 // int noldpoints = mesh->GetNP(); 892 int noldsurfel = mesh.GetNSE(); 893 894 static Timer tsurfprop("surfprop"); 895 tsurfprop.Start(); 896 GProp_GProps sprops; 897 BRepGProp::SurfaceProperties(TopoDS::Face(geom.fmap(k)),sprops); 898 tsurfprop.Stop(); 899 meshing.SetMaxArea(2.*sprops.Mass()); 900 901 MESHING2_RESULT res; 902 903 // TODO: check overlap not correctly working here 904 MeshingParameters mparam_without_overlap = mparam; 905 mparam_without_overlap.checkoverlap = false; 906 907 try { 908 static Timer t("GenerateMesh"); RegionTimer reg(t); 909 res = meshing.GenerateMesh (mesh, mparam_without_overlap, maxh, k); 910 } 911 912 catch (SingularMatrixException) 913 { 914 // (*myerr) << "Singular Matrix" << endl; 915 res = MESHING2_GIVEUP; 916 } 917 918 catch (UVBoundsException) 919 { 920 // (*myerr) << "UV bounds exceeded" << endl; 921 res = MESHING2_GIVEUP; 922 } 923 924 projecttype = PARAMETERSPACE; 925 static Timer t1("rest of loop"); RegionTimer reg1(t1); 926 927 if (res != MESHING2_OK) 928 { 929 if (notrys == 1) 930 { 931 for (SurfaceElementIndex sei = noldsurfel; sei < mesh.GetNSE(); sei++) 932 mesh.Delete(sei); 933 934 mesh.Compress(); 935 936 cout << "retry Surface " << k << endl; 937 938 k--; 939 // projecttype*=-1; 940 projecttype = PLANESPACE; 941 notrys++; 942 continue; 943 } 944 else 945 { 946 geom.facemeshstatus[k-1] = -1; 947 PrintError ("Problem in Surface mesh generation"); 948 surfmesherror++; 949 // throw NgException ("Problem in Surface mesh generation"); 950 } 951 } 952 else 953 { 954 geom.facemeshstatus[k-1] = 1; 955 } 956 957 notrys = 1; 958 959 for (SurfaceElementIndex sei = oldnf; sei < mesh.GetNSE(); sei++) 960 mesh[sei].SetIndex (k); 961 962 auto n_illegal_trigs = mesh.FindIllegalTrigs(); 963 PrintMessage (3, n_illegal_trigs, " illegal triangles"); 964 } 965 966 // ofstream problemfile("occmesh.rep"); 967 968 // problemfile << "SURFACEMESHING" << endl << endl; 969 970 if (surfmesherror) 971 { 972 cout << "WARNING! NOT ALL FACES HAVE BEEN MESHED" << endl; 973 cout << "SURFACE MESHING ERROR OCCURRED IN " << surfmesherror << " FACES:" << endl; 974 for (int i = 1; i <= geom.fmap.Extent(); i++) 975 if (geom.facemeshstatus[i-1] == -1) 976 { 977 cout << "Face " << i << endl; 978 // problemfile << "problem with face " << i << endl; 979 // problemfile << "vertices: " << endl; 980 TopExp_Explorer exp0,exp1,exp2; 981 for ( exp0.Init(TopoDS::Face (geom.fmap(i)), TopAbs_WIRE); exp0.More(); exp0.Next() ) 982 { 983 TopoDS_Wire wire = TopoDS::Wire(exp0.Current()); 984 for ( exp1.Init(wire,TopAbs_EDGE); exp1.More(); exp1.Next() ) 985 { 986 TopoDS_Edge edge = TopoDS::Edge(exp1.Current()); 987 for ( exp2.Init(edge,TopAbs_VERTEX); exp2.More(); exp2.Next() ) 988 { 989 TopoDS_Vertex vertex = TopoDS::Vertex(exp2.Current()); 990 gp_Pnt point = BRep_Tool::Pnt(vertex); 991 // problemfile << point.X() << " " << point.Y() << " " << point.Z() << endl; 992 } 993 } 994 } 995 // problemfile << endl; 996 997 } 998 cout << endl << endl; 999 cout << "for more information open IGES/STEP Topology Explorer" << endl; 1000 // problemfile.close(); 1001 throw NgException ("Problem in Surface mesh generation"); 1002 } 1003 else 1004 { 1005 // problemfile << "OK" << endl << endl; 1006 // problemfile.close(); 1007 } 1008 1009 for (int i = 0; i < mesh.GetNFD(); i++) 1010 mesh.SetBCName (i, mesh.GetFaceDescriptor(i+1).GetBCName()); 1011 multithread.task = savetask; 1012 } 1013 OCCOptimizeSurface(OCCGeometry & geom,Mesh & mesh,const MeshingParameters & mparam)1014 void OCCOptimizeSurface(OCCGeometry & geom, Mesh & mesh, 1015 const MeshingParameters & mparam) 1016 { 1017 const char * savetask = multithread.task; 1018 multithread.task = "Optimizing surface"; 1019 1020 static Timer timer_opt2d("Optimization 2D"); 1021 timer_opt2d.Start(); 1022 1023 for (int k = 1; k <= mesh.GetNFD(); k++) 1024 { 1025 // if (k != 42) continue; 1026 // if (k != 36) continue; 1027 1028 // (*testout) << "optimize face " << k << endl; 1029 multithread.percent = 100 * k / (mesh.GetNFD() + VSMALL); 1030 1031 FaceDescriptor & fd = mesh.GetFaceDescriptor(k); 1032 1033 PrintMessage (1, "Optimize Surface ", k); 1034 for (int i = 1; i <= mparam.optsteps2d; i++) 1035 { 1036 // (*testout) << "optstep " << i << endl; 1037 if (multithread.terminate) return; 1038 1039 { 1040 MeshOptimize2d meshopt(mesh); 1041 meshopt.SetFaceIndex (k); 1042 meshopt.SetImproveEdges (0); 1043 meshopt.SetMetricWeight (mparam.elsizeweight); 1044 meshopt.SetWriteStatus (0); 1045 meshopt.EdgeSwapping (i > mparam.optsteps2d/2); 1046 } 1047 1048 if (multithread.terminate) return; 1049 { 1050 MeshOptimize2d meshopt(mesh); 1051 meshopt.SetFaceIndex (k); 1052 meshopt.SetImproveEdges (0); 1053 meshopt.SetMetricWeight (mparam.elsizeweight); 1054 meshopt.SetWriteStatus (0); 1055 meshopt.ImproveMesh (mparam); 1056 } 1057 1058 { 1059 MeshOptimize2d meshopt(mesh); 1060 meshopt.SetFaceIndex (k); 1061 meshopt.SetImproveEdges (0); 1062 meshopt.SetMetricWeight (mparam.elsizeweight); 1063 meshopt.SetWriteStatus (0); 1064 meshopt.CombineImprove (); 1065 } 1066 1067 if (multithread.terminate) return; 1068 { 1069 MeshOptimize2d meshopt(mesh); 1070 meshopt.SetFaceIndex (k); 1071 meshopt.SetImproveEdges (0); 1072 meshopt.SetMetricWeight (mparam.elsizeweight); 1073 meshopt.SetWriteStatus (0); 1074 meshopt.ImproveMesh (mparam); 1075 } 1076 } 1077 1078 } 1079 1080 1081 mesh.CalcSurfacesOfNode(); 1082 mesh.Compress(); 1083 timer_opt2d.Stop(); 1084 1085 multithread.task = savetask; 1086 } 1087 1088 1089 OCCSetLocalMeshSize(const OCCGeometry & geom,Mesh & mesh,const MeshingParameters & mparam,const OCCParameters & occparam)1090 void OCCSetLocalMeshSize(const OCCGeometry & geom, Mesh & mesh, 1091 const MeshingParameters & mparam, const OCCParameters& occparam) 1092 { 1093 static Timer t1("OCCSetLocalMeshSize"); 1094 RegionTimer regt(t1); 1095 mesh.SetGlobalH (mparam.maxh); 1096 mesh.SetMinimalH (mparam.minh); 1097 1098 NgArray<double> maxhdom; 1099 maxhdom.SetSize (geom.NrSolids()); 1100 maxhdom = mparam.maxh; 1101 1102 int dom = 0; 1103 for (TopExp_Explorer e(geom.GetShape(), TopAbs_SOLID); e.More(); e.Next(), dom++) 1104 maxhdom[dom] = min2(maxhdom[dom], OCCGeometry::global_shape_properties[e.Current().TShape()].maxh); 1105 1106 mesh.SetMaxHDomain (maxhdom); 1107 1108 Box<3> bb = geom.GetBoundingBox(); 1109 bb.Increase (bb.Diam()/10); 1110 1111 mesh.SetLocalH (bb.PMin(), bb.PMax(), 0.5); 1112 1113 if (mparam.uselocalh) 1114 { 1115 const char * savetask = multithread.task; 1116 multithread.percent = 0; 1117 1118 mesh.SetLocalH (bb.PMin(), bb.PMax(), mparam.grading); 1119 1120 int nedges = geom.emap.Extent(); 1121 1122 double mincurvelength = IGNORECURVELENGTH; 1123 double maxedgelen = 0; 1124 double minedgelen = 1e99; 1125 1126 if(occparam.resthminedgelenenable) 1127 { 1128 mincurvelength = occparam.resthminedgelen; 1129 if(mincurvelength < IGNORECURVELENGTH) mincurvelength = IGNORECURVELENGTH; 1130 } 1131 1132 multithread.task = "Setting local mesh size (elements per edge)"; 1133 1134 // Philippose - 23/01/2009 1135 // Find all the parent faces of a given edge 1136 // and limit the mesh size of the edge based on the 1137 // mesh size limit of the face 1138 TopTools_IndexedDataMapOfShapeListOfShape edge_face_map; 1139 edge_face_map.Clear(); 1140 TopExp::MapShapesAndAncestors(geom.shape, TopAbs_EDGE, TopAbs_FACE, edge_face_map); 1141 1142 // setting elements per edge 1143 for (int i = 1; i <= nedges && !multithread.terminate; i++) 1144 { 1145 TopoDS_Edge e = TopoDS::Edge (geom.emap(i)); 1146 multithread.percent = 100 * (i-1)/double(nedges); 1147 if (BRep_Tool::Degenerated(e)) continue; 1148 1149 GProp_GProps system; 1150 BRepGProp::LinearProperties(e, system); 1151 double len = system.Mass(); 1152 1153 if (len < mincurvelength) 1154 { 1155 (*testout) << "ignored" << endl; 1156 continue; 1157 } 1158 1159 double localh = len/mparam.segmentsperedge; 1160 double s0, s1; 1161 1162 const TopTools_ListOfShape& parent_faces = edge_face_map.FindFromKey(e); 1163 1164 TopTools_ListIteratorOfListOfShape parent_face_list; 1165 1166 for(parent_face_list.Initialize(parent_faces); parent_face_list.More(); parent_face_list.Next()) 1167 { 1168 TopoDS_Face parent_face = TopoDS::Face(parent_face_list.Value()); 1169 1170 int face_index = geom.fmap.FindIndex(parent_face); 1171 1172 if(face_index >= 1) localh = min(localh,geom.face_maxh[face_index - 1]); 1173 localh = min2(localh, OCCGeometry::global_shape_properties[parent_face.TShape()].maxh); 1174 } 1175 1176 Handle(Geom_Curve) c = BRep_Tool::Curve(e, s0, s1); 1177 1178 localh = min2(localh, OCCGeometry::global_shape_properties[e.TShape()].maxh); 1179 maxedgelen = max (maxedgelen, len); 1180 minedgelen = min (minedgelen, len); 1181 int maxj = max((int) ceil(len/localh), 2); 1182 1183 for (int j = 0; j <= maxj; j++) 1184 { 1185 gp_Pnt pnt = c->Value (s0+double(j)/maxj*(s1-s0)); 1186 mesh.RestrictLocalH (Point3d(pnt.X(), pnt.Y(), pnt.Z()), localh); 1187 } 1188 } 1189 1190 multithread.task = "Setting local mesh size (edge curvature)"; 1191 1192 // setting edge curvature 1193 1194 int nsections = 20; 1195 1196 for (int i = 1; i <= nedges && !multithread.terminate; i++) 1197 { 1198 double maxcur = 0; 1199 multithread.percent = 100 * (i-1)/double(nedges); 1200 TopoDS_Edge edge = TopoDS::Edge (geom.emap(i)); 1201 if (BRep_Tool::Degenerated(edge)) continue; 1202 double s0, s1; 1203 Handle(Geom_Curve) c = BRep_Tool::Curve(edge, s0, s1); 1204 BRepAdaptor_Curve brepc(edge); 1205 BRepLProp_CLProps prop(brepc, 2, 1e-5); 1206 1207 for (int j = 1; j <= nsections; j++) 1208 { 1209 double s = s0 + j/(double) nsections * (s1-s0); 1210 prop.SetParameter (s); 1211 double curvature = 0; 1212 if(prop.IsTangentDefined()) 1213 curvature = prop.Curvature(); 1214 if(curvature> maxcur) maxcur = curvature; 1215 1216 if (curvature >= 1e99) 1217 continue; 1218 1219 gp_Pnt pnt = c->Value (s); 1220 1221 mesh.RestrictLocalH (Point3d(pnt.X(), pnt.Y(), pnt.Z()), ComputeH (fabs(curvature), mparam)); 1222 } 1223 } 1224 1225 multithread.task = "Setting local mesh size (face curvature)"; 1226 1227 // setting face curvature 1228 1229 int nfaces = geom.fmap.Extent(); 1230 1231 for (int i = 1; i <= nfaces && !multithread.terminate; i++) 1232 { 1233 multithread.percent = 100 * (i-1)/double(nfaces); 1234 TopoDS_Face face = TopoDS::Face(geom.fmap(i)); 1235 TopLoc_Location loc; 1236 Handle(Geom_Surface) surf = BRep_Tool::Surface (face); 1237 Handle(Poly_Triangulation) triangulation = BRep_Tool::Triangulation (face, loc); 1238 1239 if (triangulation.IsNull()) 1240 { 1241 BRepTools::Clean (geom.shape); 1242 BRepMesh_IncrementalMesh (geom.shape, 0.01, true); 1243 triangulation = BRep_Tool::Triangulation (face, loc); 1244 } 1245 1246 BRepAdaptor_Surface sf(face, Standard_True); 1247 // one prop for evaluating and one for derivatives 1248 BRepLProp_SLProps prop(sf, 0, 1e-5); 1249 BRepLProp_SLProps prop2(sf, 2, 1e-5); 1250 1251 int ntriangles = triangulation -> NbTriangles(); 1252 for (int j = 1; j <= ntriangles; j++) 1253 { 1254 gp_Pnt p[3]; 1255 gp_Pnt2d par[3]; 1256 1257 for (int k = 1; k <=3; k++) 1258 { 1259 // int n = triangulation->Triangles()(j)(k); 1260 // p[k-1] = triangulation->Nodes()(n).Transformed(loc); 1261 // par[k-1] = triangulation->UVNodes()(n); 1262 // fix for OCC7.6.0-dev 1263 int n = triangulation->Triangle(j)(k); 1264 p[k-1] = triangulation->Node(n).Transformed(loc); 1265 par[k-1] = triangulation->UVNode(n); 1266 } 1267 1268 //double maxside = 0; 1269 //maxside = max (maxside, p[0].Distance(p[1])); 1270 //maxside = max (maxside, p[0].Distance(p[2])); 1271 //maxside = max (maxside, p[1].Distance(p[2])); 1272 //cout << "\rFace " << i << " pos11 ntriangles " << ntriangles << " maxside " << maxside << flush; 1273 1274 RestrictHTriangle (par[0], par[1], par[2], &prop, &prop2, mesh, 0, 0, mparam); 1275 //cout << "\rFace " << i << " pos12 ntriangles " << ntriangles << flush; 1276 } 1277 } 1278 1279 // setting close edges 1280 1281 if (mparam.closeedgefac.has_value()) 1282 { 1283 multithread.task = "Setting local mesh size (close edges)"; 1284 1285 int sections = 100; 1286 1287 NgArray<Line> lines(sections*nedges); 1288 1289 /* 1290 BoxTree<3> * searchtree = 1291 new BoxTree<3> (bb.PMin(), bb.PMax()); 1292 */ 1293 BoxTree<3> searchtree(bb.PMin(), bb.PMax()); 1294 1295 int nlines = 0; 1296 for (int i = 1; i <= nedges && !multithread.terminate; i++) 1297 { 1298 TopoDS_Edge edge = TopoDS::Edge (geom.emap(i)); 1299 if (BRep_Tool::Degenerated(edge)) continue; 1300 1301 double s0, s1; 1302 Handle(Geom_Curve) c = BRep_Tool::Curve(edge, s0, s1); 1303 BRepAdaptor_Curve brepc(edge); 1304 BRepLProp_CLProps prop(brepc, 1, 1e-5); 1305 prop.SetParameter (s0); 1306 1307 gp_Vec d0 = prop.D1().Normalized(); 1308 double s_start = s0; 1309 int count = 0; 1310 for (int j = 1; j <= sections; j++) 1311 { 1312 double s = s0 + (s1-s0)*(double)j/(double)sections; 1313 prop.SetParameter (s); 1314 gp_Vec d1 = prop.D1().Normalized(); 1315 double cosalpha = fabs(d0*d1); 1316 if ((j == sections) || (cosalpha < cos(10.0/180.0*M_PI))) 1317 { 1318 count++; 1319 gp_Pnt p0 = c->Value (s_start); 1320 gp_Pnt p1 = c->Value (s); 1321 lines[nlines].p0 = Point<3> (p0.X(), p0.Y(), p0.Z()); 1322 lines[nlines].p1 = Point<3> (p1.X(), p1.Y(), p1.Z()); 1323 1324 Box3d box; 1325 box.SetPoint (Point3d(lines[nlines].p0)); 1326 box.AddPoint (Point3d(lines[nlines].p1)); 1327 1328 searchtree.Insert (box.PMin(), box.PMax(), nlines+1); 1329 nlines++; 1330 1331 s_start = s; 1332 d0 = d1; 1333 } 1334 } 1335 } 1336 1337 NgArray<int> linenums; 1338 1339 for (int i = 0; i < nlines; i++) 1340 { 1341 multithread.percent = (100*i)/double(nlines); 1342 Line & line = lines[i]; 1343 1344 Box3d box; 1345 box.SetPoint (Point3d(line.p0)); 1346 box.AddPoint (Point3d(line.p1)); 1347 double maxhline = max (mesh.GetH(box.PMin()), 1348 mesh.GetH(box.PMax())); 1349 box.Increase(maxhline); 1350 1351 double mindist = 1e99; 1352 linenums.SetSize(0); 1353 searchtree.GetIntersecting(box.PMin(),box.PMax(),linenums); 1354 1355 for (int j = 0; j < linenums.Size(); j++) 1356 { 1357 int num = linenums[j]-1; 1358 if (i == num) continue; 1359 if ((line.p0-lines[num].p0).Length2() < 1e-15) continue; 1360 if ((line.p0-lines[num].p1).Length2() < 1e-15) continue; 1361 if ((line.p1-lines[num].p0).Length2() < 1e-15) continue; 1362 if ((line.p1-lines[num].p1).Length2() < 1e-15) continue; 1363 mindist = min (mindist, line.Dist(lines[num])); 1364 } 1365 1366 mindist /= (*mparam.closeedgefac + VSMALL); 1367 1368 if (mindist < 1e-3 * bb.Diam()) 1369 { 1370 (*testout) << "extremely small local h: " << mindist 1371 << " --> setting to " << 1e-3 * bb.Diam() << endl; 1372 (*testout) << "somewhere near " << line.p0 << " - " << line.p1 << endl; 1373 mindist = 1e-3 * bb.Diam(); 1374 } 1375 1376 mesh.RestrictLocalHLine(line.p0, line.p1, mindist); 1377 } 1378 } 1379 1380 for (auto mspnt : mparam.meshsize_points) 1381 mesh.RestrictLocalH(mspnt.pnt, mspnt.h); 1382 1383 multithread.task = savetask; 1384 1385 } 1386 1387 mesh.LoadLocalMeshSize (mparam.meshsizefilename); 1388 } 1389 } 1390 1391 #endif 1392