1 #ifdef OCCGEOMETRY 2 3 #include <mystdlib.h> 4 5 #include <occgeom.hpp> 6 #include <meshing.hpp> 7 #include <GeomLProp_SLProps.hxx> 8 #include <ShapeAnalysis_Surface.hxx> 9 10 11 namespace netgen 12 { 13 #include "occmeshsurf.hpp" 14 15 16 bool glob_testout(false); 17 GetNormalVector(const Point<3> & p,const PointGeomInfo & geominfo,Vec<3> & n) const18 void OCCSurface :: GetNormalVector (const Point<3> & p, 19 const PointGeomInfo & geominfo, 20 Vec<3> & n) const 21 { 22 gp_Pnt pnt; 23 gp_Vec du, dv; 24 25 /* 26 double gu = geominfo.u; 27 double gv = geominfo.v; 28 29 if (fabs (gu) < 1e-3) gu = 0; 30 if (fabs (gv) < 1e-3) gv = 0; 31 32 occface->D1(gu,gv,pnt,du,dv); 33 */ 34 35 /* 36 occface->D1(geominfo.u,geominfo.v,pnt,du,dv); 37 38 n = Cross (Vec<3>(du.X(), du.Y(), du.Z()), 39 Vec<3>(dv.X(), dv.Y(), dv.Z())); 40 n.Normalize(); 41 */ 42 43 44 45 GeomLProp_SLProps lprop(occface,geominfo.u,geominfo.v,1,1e-5); 46 double setu=geominfo.u,setv=geominfo.v; 47 48 if(lprop.D1U().Magnitude() < 1e-5 || lprop.D1V().Magnitude() < 1e-5) 49 { 50 double ustep = 0.01*(umax-umin); 51 double vstep = 0.01*(vmax-vmin); 52 53 n=0; 54 55 while(setu < umax && (lprop.D1U().Magnitude() < 1e-5 || lprop.D1V().Magnitude() < 1e-5)) 56 setu += ustep; 57 if(setu < umax) 58 { 59 lprop.SetParameters(setu,setv); 60 n(0)+=lprop.Normal().X(); 61 n(1)+=lprop.Normal().Y(); 62 n(2)+=lprop.Normal().Z(); 63 } 64 setu = geominfo.u; 65 while(setu > umin && (lprop.D1U().Magnitude() < 1e-5 || lprop.D1V().Magnitude() < 1e-5)) 66 setu -= ustep; 67 if(setu > umin) 68 { 69 lprop.SetParameters(setu,setv); 70 n(0)+=lprop.Normal().X(); 71 n(1)+=lprop.Normal().Y(); 72 n(2)+=lprop.Normal().Z(); 73 } 74 setu = geominfo.u; 75 76 while(setv < vmax && (lprop.D1U().Magnitude() < 1e-5 || lprop.D1V().Magnitude() < 1e-5)) 77 setv += ustep; 78 if(setv < vmax) 79 { 80 lprop.SetParameters(setu,setv); 81 n(0)+=lprop.Normal().X(); 82 n(1)+=lprop.Normal().Y(); 83 n(2)+=lprop.Normal().Z(); 84 } 85 setv = geominfo.v; 86 while(setv > vmin && (lprop.D1U().Magnitude() < 1e-5 || lprop.D1V().Magnitude() < 1e-5)) 87 setv -= ustep; 88 if(setv > vmin) 89 { 90 lprop.SetParameters(setu,setv); 91 n(0)+=lprop.Normal().X(); 92 n(1)+=lprop.Normal().Y(); 93 n(2)+=lprop.Normal().Z(); 94 } 95 setv = geominfo.v; 96 97 n.Normalize(); 98 } 99 else 100 { 101 n(0)=lprop.Normal().X(); 102 n(1)=lprop.Normal().Y(); 103 n(2)=lprop.Normal().Z(); 104 } 105 106 if(glob_testout) 107 { 108 (*testout) << "u " << geominfo.u << " v " << geominfo.v 109 << " du " << lprop.D1U().X() << " "<< lprop.D1U().Y() << " "<< lprop.D1U().Z() 110 << " dv " << lprop.D1V().X() << " "<< lprop.D1V().Y() << " "<< lprop.D1V().Z() << endl; 111 } 112 113 114 115 if (orient == TopAbs_REVERSED) n = -1*n; 116 // (*testout) << "GetNormalVector" << endl; 117 } 118 119 DefineTangentialPlane(const Point<3> & ap1,const PointGeomInfo & geominfo1,const Point<3> & ap2,const PointGeomInfo & geominfo2)120 void OCCSurface :: DefineTangentialPlane (const Point<3> & ap1, 121 const PointGeomInfo & geominfo1, 122 const Point<3> & ap2, 123 const PointGeomInfo & geominfo2) 124 { 125 if (projecttype == PLANESPACE) 126 { 127 p1 = ap1; p2 = ap2; 128 129 //cout << "p1 = " << p1 << endl; 130 //cout << "p2 = " << p2 << endl; 131 132 GetNormalVector (p1, geominfo1, ez); 133 134 ex = p2 - p1; 135 ex -= (ex * ez) * ez; 136 ex.Normalize(); 137 ey = Cross (ez, ex); 138 139 GetNormalVector (p2, geominfo2, n2); 140 141 nmid = 0.5*(n2+ez); 142 143 ez = nmid; 144 ez.Normalize(); 145 146 ex = (p2 - p1).Normalize(); 147 ez -= (ez * ex) * ex; 148 ez.Normalize(); 149 ey = Cross (ez, ex); 150 nmid = ez; 151 //cout << "ex " << ex << " ey " << ey << " ez " << ez << endl; 152 } 153 else 154 { 155 if ( (geominfo1.u < umin) || 156 (geominfo1.u > umax) || 157 (geominfo2.u < umin) || 158 (geominfo2.u > umax) || 159 (geominfo1.v < vmin) || 160 (geominfo1.v > vmax) || 161 (geominfo2.v < vmin) || 162 (geominfo2.v > vmax) ) throw UVBoundsException(); 163 164 165 p1 = ap1; p2 = ap2; 166 psp1 = Point<2>(geominfo1.u, geominfo1.v); 167 psp2 = Point<2>(geominfo2.u, geominfo2.v); 168 169 Vec<3> n; 170 GetNormalVector (p1, geominfo1, n); 171 172 gp_Pnt pnt; 173 gp_Vec du, dv; 174 occface->D1 (geominfo1.u, geominfo1.v, pnt, du, dv); 175 176 DenseMatrix D1(3,2), D1T(2,3), DDTinv(2,2); 177 D1(0,0) = du.X(); D1(1,0) = du.Y(); D1(2,0) = du.Z(); 178 D1(0,1) = dv.X(); D1(1,1) = dv.Y(); D1(2,1) = dv.Z(); 179 180 /* 181 (*testout) << "DefineTangentialPlane" << endl 182 << "---------------------" << endl; 183 (*testout) << "D1 = " << endl << D1 << endl; 184 */ 185 186 Transpose (D1, D1T); 187 DenseMatrix D1TD1(3,3); 188 189 D1TD1 = D1T*D1; 190 if (D1TD1.Det() == 0) throw SingularMatrixException(); 191 192 CalcInverse (D1TD1, DDTinv); 193 DenseMatrix Y(3,2); 194 Vec<3> y1 = (ap2-ap1).Normalize(); 195 Vec<3> y2 = Cross(n, y1).Normalize(); 196 for (int i = 0; i < 3; i++) 197 { 198 Y(i,0) = y1(i); 199 Y(i,1) = y2(i); 200 } 201 202 DenseMatrix A(2,2); 203 A = DDTinv * D1T * Y; 204 DenseMatrix Ainv(2,2); 205 206 if (A.Det() == 0) throw SingularMatrixException(); 207 208 CalcInverse (A, Ainv); 209 210 for (int i = 0; i < 2; i++) 211 for (int j = 0; j < 2; j++) 212 { 213 Amat(i,j) = A(i,j); 214 Amatinv(i,j) = Ainv(i,j); 215 } 216 217 Vec<2> temp = Amatinv * (psp2-psp1); 218 219 220 double r = temp.Length(); 221 // double alpha = -acos (temp(0)/r); 222 double alpha = -atan2 (temp(1),temp(0)); 223 DenseMatrix R(2,2); 224 R(0,0) = cos (alpha); 225 R(1,0) = -sin (alpha); 226 R(0,1) = sin (alpha); 227 R(1,1) = cos (alpha); 228 229 230 A = A*R; 231 232 if (A.Det() == 0) throw SingularMatrixException(); 233 234 CalcInverse (A, Ainv); 235 236 237 for (int i = 0; i < 2; i++) 238 for (int j = 0; j < 2; j++) 239 { 240 Amat(i,j) = A(i,j); 241 Amatinv(i,j) = Ainv(i,j); 242 } 243 244 temp = Amatinv * (psp2-psp1); 245 246 }; 247 248 } 249 250 ToPlane(const Point<3> & p3d,const PointGeomInfo & geominfo,Point<2> & pplane,double h,int & zone) const251 void OCCSurface :: ToPlane (const Point<3> & p3d, 252 const PointGeomInfo & geominfo, 253 Point<2> & pplane, 254 double h, int & zone) const 255 { 256 if (projecttype == PLANESPACE) 257 { 258 Vec<3> p1p, n; 259 GetNormalVector (p3d, geominfo, n); 260 261 p1p = p3d - p1; 262 pplane(0) = (p1p * ex) / h; 263 pplane(1) = (p1p * ey) / h; 264 265 if (n * nmid < 0) 266 zone = -1; 267 else 268 zone = 0; 269 270 /* 271 if(zone == -1) 272 { 273 (*testout) << "zone = -1 for " << p3d << " 2D: " << pplane << " n " << n << " nmid " << nmid << endl; 274 glob_testout = true; 275 GetNormalVector (p3d, geominfo, n); 276 glob_testout = false; 277 } 278 */ 279 } 280 else 281 { 282 pplane = Point<2>(geominfo.u, geominfo.v); 283 // (*testout) << "(u,v) = " << geominfo.u << ", " << geominfo.v << endl; 284 pplane = Point<2> (1/h * (Amatinv * (pplane-psp1))); 285 // pplane = Point<2> (h * (Amatinv * (pplane-psp1))); 286 // pplane = Point<2> (1/h * ((pplane-psp1))); 287 288 zone = 0; 289 }; 290 } 291 292 FromPlane(const Point<2> & pplane,Point<3> & p3d,PointGeomInfo & gi,double h)293 void OCCSurface :: FromPlane (const Point<2> & pplane, 294 Point<3> & p3d, 295 PointGeomInfo & gi, 296 double h) 297 { 298 if (projecttype == PLANESPACE) 299 { 300 // cout << "2d : " << pplane << endl; 301 p3d = p1 + (h * pplane(0)) * ex + (h * pplane(1)) * ey; 302 // cout << "3d : " << p3d << endl; 303 Project (p3d, gi); 304 // cout << "proj : " << p3d << endl; 305 } 306 else 307 { 308 // Point<2> pspnew = Point<2>(1/h * (Amat * Vec<2>(pplane)) + Vec<2>(psp1)); 309 Point<2> pspnew = Point<2>(h * (Amat * Vec<2>(pplane)) + Vec<2>(psp1)); 310 // Point<2> pspnew = Point<2>(h * (Vec<2>(pplane)) + Vec<2>(psp1)); 311 gi.u = pspnew(0); 312 gi.v = pspnew(1); 313 gi.trignum = 1; 314 gp_Pnt val = occface->Value (gi.u, gi.v); 315 p3d = Point<3> (val.X(), val.Y(), val.Z()); 316 }; 317 } 318 319 320 Project(Point<3> & p,PointGeomInfo & gi)321 void OCCSurface :: Project (Point<3> & p, PointGeomInfo & gi) 322 { 323 // static int cnt = 0; 324 // if (cnt++ % 1000 == 0) cout << "********************************************** OCCSurfce :: Project, cnt = " << cnt << endl; 325 326 gp_Pnt pnt(p(0), p(1), p(2)); 327 328 //(*testout) << "pnt = " << pnt.X() << ", " << pnt.Y() << ", " << pnt.Z() << endl; 329 330 331 /* 332 GeomAPI_ProjectPointOnSurf proj(pnt, occface, umin, umax, vmin, vmax); 333 334 if (!proj.NbPoints()) 335 { 336 cout << "Project Point on Surface FAIL" << endl; 337 throw UVBoundsException(); 338 } 339 */ 340 341 342 343 344 345 /* 346 cout << "NP = " << proj.NbPoints() << endl; 347 348 for (int i = 1; i <= proj.NbPoints(); i++) 349 { 350 gp_Pnt pnt2 = proj.Point(i); 351 Point<3> p2 = Point<3> (pnt2.X(), pnt2.Y(), pnt2.Z()); 352 cout << i << ". p = " << p2 << ", dist = " << (p2-p).Length() << endl; 353 } 354 */ 355 356 /* 357 pnt = proj.NearestPoint(); 358 proj.LowerDistanceParameters (gi.u, gi.v); 359 */ 360 361 double u,v; 362 Handle( ShapeAnalysis_Surface ) su = new ShapeAnalysis_Surface( occface ); 363 gp_Pnt2d suval = su->ValueOfUV ( pnt, BRep_Tool::Tolerance( topods_face ) ); 364 suval.Coord( u, v); 365 pnt = occface->Value( u, v ); 366 367 //(*testout) << "pnt(proj) = " << pnt.X() << ", " << pnt.Y() << ", " << pnt.Z() << endl; 368 gi.u = u; 369 gi.v = v; 370 371 372 gi.trignum = 1; 373 374 p = Point<3> (pnt.X(), pnt.Y(), pnt.Z()); 375 } 376 377 Meshing2OCCSurfaces(const TopoDS_Shape & asurf,const Box<3> & abb,int aprojecttype)378 Meshing2OCCSurfaces :: Meshing2OCCSurfaces (const TopoDS_Shape & asurf, 379 const Box<3> & abb, int aprojecttype) 380 : Meshing2(mparam, Box<3>(abb.PMin(), abb.PMax())), surface(TopoDS::Face(asurf), aprojecttype) 381 { 382 ; 383 } 384 385 DefineTransformation(const Point3d & p1,const Point3d & p2,const PointGeomInfo * geominfo1,const PointGeomInfo * geominfo2)386 void Meshing2OCCSurfaces :: DefineTransformation (const Point3d & p1, const Point3d & p2, 387 const PointGeomInfo * geominfo1, 388 const PointGeomInfo * geominfo2) 389 { 390 ((OCCSurface&)surface).DefineTangentialPlane (p1, *geominfo1, p2, *geominfo2); 391 } 392 TransformToPlain(const Point3d & locpoint,const MultiPointGeomInfo & geominfo,Point2d & planepoint,double h,int & zone)393 void Meshing2OCCSurfaces :: TransformToPlain (const Point3d & locpoint, 394 const MultiPointGeomInfo & geominfo, 395 Point2d & planepoint, 396 double h, int & zone) 397 { 398 Point<2> hp; 399 surface.ToPlane (locpoint, geominfo.GetPGI(1), hp, h, zone); 400 planepoint.X() = hp(0); 401 planepoint.Y() = hp(1); 402 } 403 TransformFromPlain(Point2d & planepoint,Point3d & locpoint,PointGeomInfo & gi,double h)404 int Meshing2OCCSurfaces :: TransformFromPlain (Point2d & planepoint, 405 Point3d & locpoint, 406 PointGeomInfo & gi, 407 double h) 408 { 409 Point<3> hp; 410 Point<2> hp2 (planepoint.X(), planepoint.Y()); 411 surface.FromPlane (hp2, hp, gi, h); 412 locpoint = hp; 413 return 0; 414 } 415 416 417 CalcLocalH(const Point3d & p,double gh) const418 double Meshing2OCCSurfaces :: CalcLocalH (const Point3d & p, double gh) const 419 { 420 return gh; 421 } 422 423 424 425 426 427 MeshOptimize2dOCCSurfaces(const OCCGeometry & ageometry)428 MeshOptimize2dOCCSurfaces :: MeshOptimize2dOCCSurfaces (const OCCGeometry & ageometry) 429 : MeshOptimize2d(), geometry(ageometry) 430 { 431 ; 432 } 433 434 ProjectPoint(INDEX surfind,Point<3> & p) const435 void MeshOptimize2dOCCSurfaces :: ProjectPoint (INDEX surfind, Point<3> & p) const 436 { 437 geometry.Project (surfind, p); 438 } 439 440 ProjectPointGI(INDEX surfind,Point<3> & p,PointGeomInfo & gi) const441 int MeshOptimize2dOCCSurfaces :: ProjectPointGI (INDEX surfind, Point<3> & p, PointGeomInfo & gi) const 442 { 443 double u = gi.u; 444 double v = gi.v; 445 446 Point<3> hp = p; 447 if (geometry.FastProject (surfind, hp, u, v)) 448 { 449 p = hp; 450 return 1; 451 } 452 ProjectPoint (surfind, p); 453 return CalcPointGeomInfo (surfind, gi, p); 454 } 455 456 ProjectPoint2(INDEX surfind,INDEX surfind2,Point<3> & p) const457 void MeshOptimize2dOCCSurfaces :: ProjectPoint2 (INDEX surfind, INDEX surfind2, 458 Point<3> & p) const 459 { 460 TopExp_Explorer exp0, exp1; 461 bool done = false; 462 Handle(Geom_Curve) c; 463 464 for (exp0.Init(geometry.fmap(surfind), TopAbs_EDGE); !done && exp0.More(); exp0.Next()) 465 for (exp1.Init(geometry.fmap(surfind2), TopAbs_EDGE); !done && exp1.More(); exp1.Next()) 466 { 467 if (TopoDS::Edge(exp0.Current()).IsSame(TopoDS::Edge(exp1.Current()))) 468 { 469 done = true; 470 double s0, s1; 471 c = BRep_Tool::Curve(TopoDS::Edge(exp0.Current()), s0, s1); 472 } 473 } 474 475 gp_Pnt pnt(p(0), p(1), p(2)); 476 GeomAPI_ProjectPointOnCurve proj(pnt, c); 477 pnt = proj.NearestPoint(); 478 p(0) = pnt.X(); 479 p(1) = pnt.Y(); 480 p(2) = pnt.Z(); 481 482 } 483 484 void MeshOptimize2dOCCSurfaces :: GetNormalVector(INDEX surfind,const Point<3> & p,PointGeomInfo & geominfo,Vec<3> & n) const485 GetNormalVector(INDEX surfind, const Point<3> & p, PointGeomInfo & geominfo, Vec<3> & n) const 486 { 487 gp_Pnt pnt; 488 gp_Vec du, dv; 489 490 Handle(Geom_Surface) occface; 491 occface = BRep_Tool::Surface(TopoDS::Face(geometry.fmap(surfind))); 492 493 occface->D1(geominfo.u,geominfo.v,pnt,du,dv); 494 495 n = Cross (Vec<3>(du.X(), du.Y(), du.Z()), 496 Vec<3>(dv.X(), dv.Y(), dv.Z())); 497 n.Normalize(); 498 499 if (geometry.fmap(surfind).Orientation() == TopAbs_REVERSED) n = -1*n; 500 501 // GetNormalVector (surfind, p, n); 502 } 503 504 505 void MeshOptimize2dOCCSurfaces :: GetNormalVector(INDEX surfind,const Point<3> & p,Vec<3> & n) const506 GetNormalVector(INDEX surfind, const Point<3> & p, Vec<3> & n) const 507 { 508 // static int cnt = 0; 509 // if (cnt++ % 1000 == 0) cout << "GetNV cnt = " << cnt << endl; 510 Standard_Real u,v; 511 512 gp_Pnt pnt(p(0), p(1), p(2)); 513 514 Handle(Geom_Surface) occface; 515 occface = BRep_Tool::Surface(TopoDS::Face(geometry.fmap(surfind))); 516 517 /* 518 GeomAPI_ProjectPointOnSurf proj(pnt, occface); 519 520 if (proj.NbPoints() < 1) 521 { 522 cout << "ERROR: OCCSurface :: GetNormalVector: GeomAPI_ProjectPointOnSurf failed!" 523 << endl; 524 cout << p << endl; 525 return; 526 } 527 528 proj.LowerDistanceParameters (u, v); 529 */ 530 531 Handle( ShapeAnalysis_Surface ) su = new ShapeAnalysis_Surface( occface ); 532 gp_Pnt2d suval = su->ValueOfUV ( pnt, BRep_Tool::Tolerance( TopoDS::Face(geometry.fmap(surfind)) ) ); 533 suval.Coord( u, v); 534 pnt = occface->Value( u, v ); 535 536 537 538 gp_Vec du, dv; 539 occface->D1(u,v,pnt,du,dv); 540 541 /* 542 if (!occface->IsCNu (1) || !occface->IsCNv (1)) 543 (*testout) << "SurfOpt: Differentiation FAIL" << endl; 544 */ 545 546 n = Cross (Vec3d(du.X(), du.Y(), du.Z()), 547 Vec3d(dv.X(), dv.Y(), dv.Z())); 548 n.Normalize(); 549 550 if (geometry.fmap(surfind).Orientation() == TopAbs_REVERSED) n = -1*n; 551 } 552 553 554 int MeshOptimize2dOCCSurfaces :: CalcPointGeomInfo(int surfind,PointGeomInfo & gi,const Point<3> & p) const555 CalcPointGeomInfo(int surfind, PointGeomInfo& gi, const Point<3> & p) const 556 { 557 Standard_Real u,v; 558 559 gp_Pnt pnt(p(0), p(1), p(2)); 560 561 Handle(Geom_Surface) occface; 562 occface = BRep_Tool::Surface(TopoDS::Face(geometry.fmap(surfind))); 563 564 /* 565 GeomAPI_ProjectPointOnSurf proj(pnt, occface); 566 567 if (proj.NbPoints() < 1) 568 { 569 cout << "ERROR: OCCSurface :: GetNormalVector: GeomAPI_ProjectPointOnSurf failed!" 570 << endl; 571 cout << p << endl; 572 return 0; 573 } 574 575 proj.LowerDistanceParameters (u, v); 576 */ 577 578 Handle( ShapeAnalysis_Surface ) su = new ShapeAnalysis_Surface( occface ); 579 gp_Pnt2d suval = su->ValueOfUV ( pnt, BRep_Tool::Tolerance( TopoDS::Face(geometry.fmap(surfind)) ) ); 580 suval.Coord( u, v); 581 //pnt = occface->Value( u, v ); 582 583 584 gi.u = u; 585 gi.v = v; 586 return 1; 587 } 588 589 590 591 592 593 OCCRefinementSurfaces(const OCCGeometry & ageometry)594 OCCRefinementSurfaces :: OCCRefinementSurfaces (const OCCGeometry & ageometry) 595 : Refinement(), geometry(ageometry) 596 { 597 ; 598 } 599 ~OCCRefinementSurfaces()600 OCCRefinementSurfaces :: ~OCCRefinementSurfaces () 601 { 602 ; 603 } 604 605 /* 606 inline double Det3 (double a00, double a01, double a02, 607 double a10, double a11, double a12, 608 double a20, double a21, double a22) 609 { 610 return a00*a11*a22 + a01*a12*a20 + a10*a21*a02 - a20*a11*a02 - a10*a01*a22 - a21*a12*a00; 611 } 612 613 bool ProjectToSurface (gp_Pnt & p, Handle(Geom_Surface) surface, double& u, double& v) 614 { 615 gp_Pnt x = surface->Value (u,v); 616 617 if (p.SquareDistance(x) <= sqr(PROJECTION_TOLERANCE)) return true; 618 619 gp_Vec du, dv; 620 621 surface->D1(u,v,x,du,dv); 622 623 int count = 0; 624 625 gp_Pnt xold; 626 gp_Vec n; 627 double det, lambda, mu; 628 629 do { 630 count++; 631 632 n = du^dv; 633 634 det = Det3 (n.X(), du.X(), dv.X(), 635 n.Y(), du.Y(), dv.Y(), 636 n.Z(), du.Z(), dv.Z()); 637 638 if (det < 1e-15) return false; 639 640 lambda = Det3 (n.X(), p.X()-x.X(), dv.X(), 641 n.Y(), p.Y()-x.Y(), dv.Y(), 642 n.Z(), p.Z()-x.Z(), dv.Z())/det; 643 644 mu = Det3 (n.X(), du.X(), p.X()-x.X(), 645 n.Y(), du.Y(), p.Y()-x.Y(), 646 n.Z(), du.Z(), p.Z()-x.Z())/det; 647 648 u += lambda; 649 v += mu; 650 651 xold = x; 652 surface->D1(u,v,x,du,dv); 653 654 } while (xold.SquareDistance(x) > sqr(PROJECTION_TOLERANCE) || count > 50); 655 656 if (count > 50) return false; 657 658 p = x; 659 660 return true; 661 } 662 */ 663 664 void OCCRefinementSurfaces :: PointBetween(const Point<3> & p1,const Point<3> & p2,double secpoint,int surfi,const PointGeomInfo & gi1,const PointGeomInfo & gi2,Point<3> & newp,PointGeomInfo & newgi) const665 PointBetween (const Point<3> & p1, const Point<3> & p2, double secpoint, 666 int surfi, 667 const PointGeomInfo & gi1, 668 const PointGeomInfo & gi2, 669 Point<3> & newp, PointGeomInfo & newgi) const 670 { 671 Point<3> hnewp; 672 hnewp = p1+secpoint*(p2-p1); 673 674 if (surfi > 0) 675 { 676 677 double u = gi1.u+secpoint*(gi2.u-gi1.u); 678 double v = gi1.v+secpoint*(gi2.v-gi1.v); 679 680 if (!geometry.FastProject (surfi, hnewp, u, v)) 681 { 682 // cout << "Fast projection to surface fails! Using OCC projection" << endl; 683 geometry.Project (surfi, hnewp); 684 } 685 686 newgi.trignum = 1; 687 newgi.u = u; 688 newgi.v = v; 689 } 690 691 newp = hnewp; 692 } 693 694 695 void OCCRefinementSurfaces :: PointBetween(const Point<3> & p1,const Point<3> & p2,double secpoint,int surfi1,int surfi2,const EdgePointGeomInfo & ap1,const EdgePointGeomInfo & ap2,Point<3> & newp,EdgePointGeomInfo & newgi) const696 PointBetween (const Point<3> & p1, const Point<3> & p2, double secpoint, 697 int surfi1, int surfi2, 698 const EdgePointGeomInfo & ap1, 699 const EdgePointGeomInfo & ap2, 700 Point<3> & newp, EdgePointGeomInfo & newgi) const 701 { 702 double s0, s1; 703 704 Point<3> hnewp = p1+secpoint*(p2-p1); 705 gp_Pnt pnt(hnewp(0), hnewp(1), hnewp(2)); 706 GeomAPI_ProjectPointOnCurve proj(pnt, BRep_Tool::Curve(TopoDS::Edge(geometry.emap(ap1.edgenr)), s0, s1)); 707 pnt = proj.NearestPoint(); 708 hnewp = Point<3> (pnt.X(), pnt.Y(), pnt.Z()); 709 newp = hnewp; 710 newgi = ap1; 711 }; 712 713 ProjectToSurface(Point<3> & p,int surfi) const714 void OCCRefinementSurfaces :: ProjectToSurface (Point<3> & p, int surfi) const 715 { 716 if (surfi > 0) 717 geometry.Project (surfi, p); 718 }; 719 ProjectToSurface(Point<3> & p,int surfi,PointGeomInfo & gi) const720 void OCCRefinementSurfaces :: ProjectToSurface (Point<3> & p, int surfi, PointGeomInfo & gi) const 721 { 722 if (surfi > 0) 723 if (!geometry.FastProject (surfi, p, gi.u, gi.v)) 724 { 725 cout << "Fast projection to surface fails! Using OCC projection" << endl; 726 geometry.Project (surfi, p); 727 } 728 }; 729 730 731 732 } 733 734 735 #endif 736