src/NETGEN/netgen49ForSalome.patch

   1 diff -Naur netgen-4.9.13_orig/libsrc/meshing/meshtype.cpp netgen-4.9.13_new/libsrc/meshing/meshtype.cpp
   2 --- netgen-4.9.13_orig/libsrc/meshing/meshtype.cpp      2009-09-13 14:28:38.000000000 +0400
   3 +++ netgen-4.9.13_new/libsrc/meshing/meshtype.cpp       2011-06-03 17:54:39.000000000 +0400
   4 @@ -1,4 +1,5 @@
   5  #include <mystdlib.h>
   6 +#include <float.h> // to get DBL_MIN defined
   7
   8  #include "meshing.hpp"
   9
  10 @@ -650,7 +651,8 @@
  11
  12          double det = trans.Det();
  13
  14 -        if (det <= 0)
  15 +        // if (det <= 0)
  16 +        if (det <= DBL_MIN) // avoid FPE
  17            err += 1e12;
  18          else
  19            err += frob * frob / det;
  20 @@ -706,7 +708,8 @@
  21
  22              double det = trans(0,0)*trans(1,1)-trans(1,0)*trans(0,1);
  23
  24 -            if (det <= 0)
  25 +            // if (det <= 0)
  26 +            if (det <= DBL_MIN)  // avoid FPE
  27                {
  28                  dd = 0;
  29                  return 1e12;
  30 @@ -790,7 +793,8 @@
  31            = dtrans(0,0) * trans(1,1) - trans(0,1) * dtrans(1,0)
  32            + trans(0,0) * dtrans(1,1) - dtrans(0,1) * trans(1,0);
  33
  34 -        if (det <= 0)
  35 +        // if (det <= 0)
  36 +        if (det <= DBL_MIN) // avoid FPE
  37            err += 1e12;
  38          else
  39            {
  40 @@ -840,7 +844,8 @@
  41          frob /= 2;
  42
  43          double det = trans.Det();
  44 -        if (det <= 0)
  45 +        //if (det <= 0)
  46 +        if (det <= DBL_MIN) // avoid FPE
  47            err += 1e12;
  48          else
  49            err += frob * frob / det;
  50 @@ -1857,7 +1862,8 @@
  51        case PYRAMID:
  52          {
  53            double noz = 1-p(2);
  54 -          if (noz == 0.0) noz = 1e-10;
  55 +          //if (noz == 0.0) noz = 1e-10;
  56 +          if (noz <= DBL_MIN) noz = 1e-10; // avoid FPE
  57
  58            double xi  = p(0) / noz;
  59            double eta = p(1) / noz;
  60 @@ -2035,7 +2041,8 @@
  61
  62          double det = -trans.Det();
  63
  64 -        if (det <= 0)
  65 +        //if (det <= 0)
  66 +        if (det <= DBL_MIN) // avoid FPE
  67            err += 1e12;
  68          else
  69            err += frob * frob * frob / det;
  70 @@ -2107,7 +2114,8 @@
  71          ddet *= -1;
  72
  73
  74 -        if (det <= 0)
  75 +        //if (det <= 0)
  76 +        if (det <= DBL_MIN) // avoid FPE
  77            err += 1e12;
  78          else
  79            {
  80 @@ -2189,7 +2197,7 @@
  81
  82          det *= -1;
  83
  84 -        if (det <= 0)
  85 +        if (det <= DBL_MIN)
  86            err += 1e12;
  87          else
  88            {
  89 @@ -2522,10 +2530,10 @@
  90
  91    MeshingParameters :: MeshingParameters ()
  92    {
  93 -    optimize3d = "cmdmustm";
  94 +    optimize3d = (char*)"cmdmustm"; // optimize3d = "cmdmustm";
  95      //optimize3d = "cmdmstm";
  96      optsteps3d = 3;
  97 -    optimize2d = "smsmsmSmSmSm";
  98 +    optimize2d = (char*)"smsmsmSmSmSm"; // optimize2d = "smsmsmSmSmSm";
  99      optsteps2d = 3;
 100      opterrpow = 2;
 101      blockfill = 1;
 102 diff -Naur netgen-4.9.13_orig/libsrc/meshing/meshtype.hpp netgen-4.9.13_new/libsrc/meshing/meshtype.hpp
 103 --- netgen-4.9.13_orig/libsrc/meshing/meshtype.hpp      2009-11-09 13:50:43.000000000 +0300
 104 +++ netgen-4.9.13_new/libsrc/meshing/meshtype.hpp       2011-06-03 17:56:43.000000000 +0400
 105 @@ -12,6 +12,7 @@
 106      Classes for NETGEN
 107  */
 108
 109 +class Mesh; // added due to compilation errors on some platforms
 110
 111  enum ELEMENT_TYPE {
 112    SEGMENT = 1, SEGMENT3 = 2,
 113 diff -Naur netgen-4.9.13_orig/libsrc/meshing/smoothing2.cpp netgen-4.9.13_new/libsrc/meshing/smoothing2.cpp
 114 --- netgen-4.9.13_orig/libsrc/meshing/smoothing2.cpp    2009-11-09 13:47:09.000000000 +0300
 115 +++ netgen-4.9.13_new/libsrc/meshing/smoothing2.cpp     2011-06-03 18:04:35.000000000 +0400
 116 @@ -302,7 +302,8 @@
 117      vgrad = 0;
 118      badness = 0;
 119
 120 -    meshthis -> GetNormalVector (surfi, sp1, gi1, n);
 121 +    //normal already computed: meshthis -> GetNormalVector (surfi, sp1, gi1, n);
 122 +    n = normal;
 123      pp1 = sp1 + x(0) * t1 + x(1) * t2;
 124
 125      //  meshthis -> ProjectPoint (surfi, pp1);
 126 @@ -360,7 +361,8 @@
 127      vgrad = 0;
 128      badness = 0;
 129
 130 -    meshthis -> GetNormalVector (surfi, sp1, gi1, n);
 131 +    //normal already computed: meshthis -> GetNormalVector (surfi, sp1, gi1, n);
 132 +    n = normal;
 133
 134      pp1 = sp1 + x(0) * t1 + x(1) * t2;
 135
 136 @@ -514,7 +516,8 @@
 137      vgrad = 0;
 138      badness = 0;
 139
 140 -    meshthis -> GetNormalVector (surfi, sp1, gi1, n);
 141 +    //normal already computed: meshthis -> GetNormalVector (surfi, sp1, gi1, n);
 142 +    n = normal;
 143
 144      pp1 = sp1 + x(0) * t1 + x(1) * t2;
 145
 146 @@ -586,7 +589,8 @@
 147      vgrad = 0;
 148      badness = 0;
 149
 150 -    meshthis -> GetNormalVector (surfi, sp1, gi1, n);
 151 +    //normal already computed: meshthis -> GetNormalVector (surfi, sp1, gi1, n);
 152 +    n = normal;
 153
 154      // pp1 = sp1;
 155      //    pp1.Add2 (x.Get(1), t1, x.Get(2), t2);
 156 @@ -973,7 +977,7 @@
 157                 {
 158                   mesh[pi] = Point<3> (origp);
 159                 }
 160 -
 161 +             break; // exit as <fact> is not used anymore
 162             }
 163         }
 164
 165 diff -Naur netgen-4.9.13_orig/libsrc/occ/occconstruction.cpp netgen-4.9.13_new/libsrc/occ/occconstruction.cpp
 166 --- netgen-4.9.13_orig/libsrc/occ/occconstruction.cpp   2009-08-24 06:32:47.000000000 +0400
 167 +++ netgen-4.9.13_new/libsrc/occ/occconstruction.cpp    2011-04-20 13:20:56.000000000 +0400
 168 @@ -28,7 +28,7 @@
 169  #include <BRepAlgoAPI_Common.hxx>
 170  #include <BRepAlgoAPI_Fuse.hxx>
 171  #include <BRepAlgoAPI_Section.hxx>
 172 -#include <BRepOffsetAPI_Sewing.hxx>
 173 +//#include <BRepOffsetAPI_Sewing.hxx>
 174  //#include <BRepAlgo_Sewing.hxx>
 175  #include <BRepOffsetAPI_MakeOffsetShape.hxx>
 176  #include <ShapeFix_Shape.hxx>
 177 diff -Naur netgen-4.9.13_orig/libsrc/occ/occgenmesh.cpp netgen-4.9.13_new/libsrc/occ/occgenmesh.cpp
 178 --- netgen-4.9.13_orig/libsrc/occ/occgenmesh.cpp        2010-03-16 09:30:07.000000000 +0300
 179 +++ netgen-4.9.13_new/libsrc/occ/occgenmesh.cpp 2011-06-28 12:04:32.000000000 +0400
 180 @@ -15,6 +15,8 @@
 181
 182  #define DIVIDEEDGESECTIONS 1000
 183  #define IGNORECURVELENGTH 1e-4
 184 +// a small value used to avoid FPE
 185 +#define VSMALL 1e-10
 186
 187
 188     bool merge_solids = 1;
 189 @@ -26,7 +28,8 @@
 190        double nq = n*q;
 191
 192        Point<3> p = p0 + 0.5*n;
 193 -      double lambda = (p-l.p0)*n / nq;
 194 +      // double lambda = (p-l.p0)*n / nq;  -- avoid FPE
 195 +      double lambda = (fabs(nq) > 1e-10) ? (p-l.p0)*n / nq : -1;
 196
 197        if (lambda >= 0 && lambda <= 1)
 198        {
 199 @@ -54,6 +57,8 @@
 200
 201
 202
 203 +
 204 +   static // useless out of this file
 205     double ComputeH (double kappa)
 206     {
 207        double hret;
 208 @@ -62,7 +67,8 @@
 209        if (mparam.maxh * kappa < 1)
 210           hret = mparam.maxh;
 211        else
 212 -         hret = 1 / kappa;
 213 +        // hret = 1 / kappa; -- avoid FPE
 214 +        hret = 1 / (kappa + VSMALL);
 215
 216        if (mparam.maxh < hret)
 217           hret = mparam.maxh;
 218 @@ -71,8 +77,7 @@
 219     }
 220
 221
 222 -
 223 -
 224 +   static // useless out of this file
 225     void RestrictHTriangle (gp_Pnt2d & par0, gp_Pnt2d & par1, gp_Pnt2d & par2,
 226                             BRepLProp_SLProps * prop, Mesh & mesh, int depth, double h = 0)
 227     {
 228 @@ -168,8 +173,8 @@
 229           if(h < 1e-4*maxside)
 230              return;
 231
 232 -
 233 -         if (h > 30) return;
 234 +         // commented to restrict H on a large sphere for example
 235 +         //if (h > 30) return;
 236        }
 237
 238        if (h < maxside && depth < 10)
 239 @@ -228,6 +233,7 @@
 240
 241
 242
 243 +   static // useless out of this file
 244     void DivideEdge (TopoDS_Edge & edge, Array<MeshPoint> & ps,
 245                      Array<double> & params, Mesh & mesh)
 246     {
 247 @@ -247,8 +253,8 @@
 248        hvalue[0] = 0;
 249        pnt = c->Value(s0);
 250
 251 -      double olddist = 0;
 252 -      double dist = 0;
 253 +      //double olddist = 0; -- useless variables
 254 +      //double dist = 0;
 255
 256        int tmpVal = (int)(DIVIDEEDGESECTIONS);
 257
 258 @@ -256,15 +262,19 @@
 259        {
 260           oldpnt = pnt;
 261           pnt = c->Value(s0+(i/double(DIVIDEEDGESECTIONS))*(s1-s0));
 262 +         // -- no more than 1 segment per <edge length>/DIVIDEEDGESECTIONS
 263           hvalue[i] = hvalue[i-1] +
 264 -            1.0/mesh.GetH(Point3d(pnt.X(), pnt.Y(), pnt.Z()))*
 265 -            pnt.Distance(oldpnt);
 266 +         //   1.0/mesh.GetH(Point3d(pnt.X(), pnt.Y(), pnt.Z()))*
 267 +         //   pnt.Distance(oldpnt);
 268 +           min( 1.0,
 269 +                1.0/mesh.GetH(Point3d(pnt.X(), pnt.Y(), pnt.Z()))*
 270 +                pnt.Distance(oldpnt));
 271
 272           //(*testout) << "mesh.GetH(Point3d(pnt.X(), pnt.Y(), pnt.Z())) " << mesh.GetH(Point3d(pnt.X(), pnt.Y(), pnt.Z()))
 273           //       <<  " pnt.Distance(oldpnt) " << pnt.Distance(oldpnt) << endl;
 274
 275 -         olddist = dist;
 276 -         dist = pnt.Distance(oldpnt);
 277 +         //olddist = dist; -- useless variables
 278 +         //dist = pnt.Distance(oldpnt);
 279        }
 280
 281        //  nsubedges = int(ceil(hvalue[DIVIDEEDGESECTIONS]));
 282 @@ -279,7 +289,10 @@
 283        {
 284           if (hvalue[i1]/hvalue[DIVIDEEDGESECTIONS]*nsubedges >= i)
 285           {
 286 -            params[i] = s0+(i1/double(DIVIDEEDGESECTIONS))*(s1-s0);
 287 +            // -- for nsubedges comparable to DIVIDEEDGESECTIONS
 288 +            //params[i] = s0+(i1/double(DIVIDEEDGESECTIONS))*(s1-s0);
 289 +            double d1 = i1 - (hvalue[i1] - i*hvalue[DIVIDEEDGESECTIONS]/nsubedges)/(hvalue[i1]-hvalue[i1-1]);
 290 +            params[i] = s0+(d1/double(DIVIDEEDGESECTIONS))*(s1-s0);
 291              pnt = c->Value(params[i]);
 292              ps[i-1] = MeshPoint (Point3d(pnt.X(), pnt.Y(), pnt.Z()));
 293              i++;
 294 @@ -323,6 +336,7 @@
 295        (*testout) << "nedges = " << nedges << endl;
 296
 297        double eps = 1e-6 * geom.GetBoundingBox().Diam();
 298 +      const double eps2 = eps * eps; // -- small optimization
 299
 300        for (int i = 1; i <= nvertices; i++)
 301        {
 302 @@ -332,7 +346,8 @@
 303           bool exists = 0;
 304           if (merge_solids)
 305              for (PointIndex pi = 1; pi <= mesh.GetNP(); pi++)
 306 -               if ( Dist2 (mesh[pi], Point<3>(mp)) < eps*eps)
 307 +               //if ( Dist2 (mesh[pi], Point<3>(mp)) < eps*eps)
 308 +               if ( Dist2 (mesh[pi], Point<3>(mp)) < eps2 ) // -- small optimization
 309                 {
 310                    exists = 1;
 311                    break;
 312 @@ -362,6 +377,7 @@
 313           {
 314              TopoDS_Face face = TopoDS::Face(exp1.Current());
 315              int facenr = geom.fmap.FindIndex(face);
 316 +            if ( facenr < 1 ) continue; // -- to support SALOME sub-meshes
 317
 318              if (face2solid[0][facenr-1] == 0)
 319                 face2solid[0][facenr-1] = solidnr;
 320 @@ -381,6 +397,7 @@
 321        int facenr = 0;
 322        int edgenr = 0;
 323
 324 +      edgenr = mesh.GetNSeg(); // to support SALOME sub-meshes
 325
 326        (*testout) << "faces = " << geom.fmap.Extent() << endl;
 327        int curr = 0;
 328 @@ -442,6 +459,7 @@
 329                    //(*testout) << "ignoring degenerated edge" << endl;
 330                    continue;
 331                 }
 332 +               if ( geom.emap.FindIndex(edge) < 1 ) continue; // to support SALOME sub-meshes
 333
 334                 if (geom.vmap.FindIndex(TopExp::FirstVertex (edge)) ==
 335                    geom.vmap.FindIndex(TopExp::LastVertex (edge)))
 336 @@ -481,13 +499,45 @@
 337                 {
 338                    Point<3> fp = occ2ng (BRep_Tool::Pnt (TopExp::FirstVertex (edge)));
 339                    Point<3> lp = occ2ng (BRep_Tool::Pnt (TopExp::LastVertex (edge)));
 340 +                  double tol2 = std::min( eps*eps, 1e-6 * Dist2( fp, lp ));
 341
 342                    pnums[0] = -1;
 343                    pnums.Last() = -1;
 344                    for (PointIndex pi = 1; pi < first_ep; pi++)
 345                    {
 346 -                     if (Dist2 (mesh[pi], fp) < eps*eps) pnums[0] = pi;
 347 -                     if (Dist2 (mesh[pi], lp) < eps*eps) pnums.Last() = pi;
 348 +                    if (Dist2 (mesh[pi], fp) < tol2) pnums[0] = pi;
 349 +                    if (Dist2 (mesh[pi], lp) < tol2) pnums.Last() = pi;
 350 +                  }
 351 +                  if ( pnums[0] == pnums.Last() )
 352 +                    pnums[0] = -1;
 353 +                  if ( pnums[0] == -1 || pnums.Last() == -1 )
 354 +                  {
 355 +                    // take into account a possible large gap between a vertex and an edge curve
 356 +                    // and a large vertex tolerance covering the whole edge
 357 +                    if ( pnums[0] == -1 )
 358 +                    {
 359 +                      double tol = BRep_Tool::Tolerance( TopExp::FirstVertex (edge));
 360 +                      for (PointIndex pi = 1; pi < first_ep; pi++)
 361 +                        if (pi != pnums.Last() && Dist2 (mesh[pi], fp) < 2*tol*tol)
 362 +                          pnums[0] = pi;
 363 +
 364 +                      if ( pnums[0] == -1 )
 365 +                        pnums[0] = geom.vmap.FindIndex (TopExp::FirstVertex (edge));
 366 +                    }
 367 +                    if ( pnums.Last() == -1 )
 368 +                    {
 369 +                      double tol = BRep_Tool::Tolerance( TopExp::LastVertex (edge));
 370 +                      for (PointIndex pi = 1; pi < first_ep; pi++)
 371 +                        if (pi != pnums[0] && Dist2 (mesh[pi], lp) < 2*tol*tol)
 372 +                          pnums.Last() = pi;
 373 +
 374 +                      if ( pnums.Last() == -1 )
 375 +                        pnums.Last() = geom.vmap.FindIndex (TopExp::LastVertex (edge));
 376 +                    }
 377 +
 378 +                    if ( Dist2( fp, mesh[PointIndex(pnums[0])]) >
 379 +                         Dist2( lp, mesh[PointIndex(pnums.Last())]))
 380 +                      std::swap( pnums[0], pnums.Last() );
 381                    }
 382                 }
 383
 384 @@ -633,7 +683,8 @@
 385           }
 386
 387           (*testout) << "mesh face " << k << endl;
 388 -         multithread.percent = 100 * k / (mesh.GetNFD()+1e-10);
 389 +         // multithread.percent = 100 * k / (mesh.GetNFD()+1e-10); -- avoid FPE
 390 +         multithread.percent = 100 * k / (mesh.GetNFD() + VSMALL);
 391           geom.facemeshstatus[k-1] = -1;
 392
 393
 394 @@ -901,7 +952,8 @@
 395           //      if (k != 36) continue;
 396
 397           //      (*testout) << "optimize face " << k << endl;
 398 -         multithread.percent = 100 * k / (mesh.GetNFD()+1e-10);
 399 +         //multithread.percent = 100 * k / (mesh.GetNFD()+1e-10); -- avoid FPE
 400 +         multithread.percent = 100 * k / (mesh.GetNFD() + VSMALL);
 401
 402           FaceDescriptor & fd = mesh.GetFaceDescriptor(k);
 403
 404 @@ -1456,3 +1508,4 @@
 405  }
 406
 407  #endif
 408 +
 409 diff -Naur netgen-4.9.13_orig/libsrc/occ/occgeom.cpp netgen-4.9.13_new/libsrc/occ/occgeom.cpp
 410 --- netgen-4.9.13_orig/libsrc/occ/occgeom.cpp   2010-03-05 16:16:21.000000000 +0300
 411 +++ netgen-4.9.13_new/libsrc/occ/occgeom.cpp    2011-06-03 19:32:06.000000000 +0400
 412 @@ -8,6 +8,8 @@
 413  #include "ShapeAnalysis_CheckSmallFace.hxx"
 414  #include "ShapeAnalysis_DataMapOfShapeListOfReal.hxx"
 415  #include "ShapeAnalysis_Surface.hxx"
 416 +#include <BRepTopAdaptor_FClass2d.hxx> // -- to optimize Project() and FastProject()
 417 +#include <TopAbs_State.hxx>
 418  #include "BRepAlgoAPI_Fuse.hxx"
 419  #include "BRepCheck_Analyzer.hxx"
 420  #include "BRepLib.hxx"
 421 @@ -16,10 +18,17 @@
 422  #include "ShapeFix_FixSmallFace.hxx"
 423  #include "Partition_Spliter.hxx"
 424
 425 -
 426  namespace netgen
 427  {
 428 -   void OCCGeometry :: PrintNrShapes ()
 429 +  // free data used to optimize Project() and FastProject()
 430 +  OCCGeometry::~OCCGeometry()
 431 +  {
 432 +    NCollection_DataMap<int,BRepTopAdaptor_FClass2d*>::Iterator it(fclsmap);
 433 +    for (; it.More(); it.Next())
 434 +      delete it.Value();
 435 +  }
 436 +
 437 +  void OCCGeometry :: PrintNrShapes ()
 438     {
 439        TopExp_Explorer e;
 440        int count = 0;
 441 @@ -951,25 +960,58 @@
 442     }
 443
 444
 445 +   // returns a projector and a classifier for the given surface
 446 +   void OCCGeometry::GetFaceTools(int surfi, Handle(ShapeAnalysis_Surface)& proj,
 447 +                                  BRepTopAdaptor_FClass2d*& cls) const
 448 +   {
 449 +     //MSV: organize caching projector in the map
 450 +     if (fprjmap.IsBound(surfi))
 451 +     {
 452 +       proj = fprjmap.Find(surfi);
 453 +       cls = fclsmap.Find(surfi);
 454 +     }
 455 +     else
 456 +     {
 457 +       const TopoDS_Face& aFace = TopoDS::Face(fmap(surfi));
 458 +       Handle(Geom_Surface) aSurf = BRep_Tool::Surface(aFace);
 459 +       proj = new ShapeAnalysis_Surface(aSurf);
 460 +       fprjmap.Bind(surfi, proj);
 461 +       cls = new BRepTopAdaptor_FClass2d(aFace,Precision::Confusion());
 462 +       fclsmap.Bind(surfi, cls);
 463 +     }
 464 +   }
 465
 466 -
 467 -   void OCCGeometry :: Project (int surfi, Point<3> & p) const
 468 +   // void OCCGeometry :: Project (int surfi, Point<3> & p) const
 469 +   bool OCCGeometry :: Project (int surfi, Point<3> & p, double& u, double& v) const
 470     {
 471        static int cnt = 0;
 472        if (++cnt % 1000 == 0) cout << "Project cnt = " << cnt << endl;
 473
 474        gp_Pnt pnt(p(0), p(1), p(2));
 475
 476 -      double u,v;
 477 -      Handle( Geom_Surface ) thesurf = BRep_Tool::Surface(TopoDS::Face(fmap(surfi)));
 478 -      Handle( ShapeAnalysis_Surface ) su = new ShapeAnalysis_Surface( thesurf );
 479 -      gp_Pnt2d suval = su->ValueOfUV ( pnt, BRep_Tool::Tolerance( TopoDS::Face(fmap(surfi)) ) );
 480 -      suval.Coord( u, v);
 481 -      pnt = thesurf->Value( u, v );
 482 -
 483 -
 484 +      // -- Optimization: use cached projector and classifier
 485 +      // double u,v;
 486 +      // Handle( Geom_Surface ) thesurf = BRep_Tool::Surface(TopoDS::Face(fmap(surfi)));
 487 +      // Handle( ShapeAnalysis_Surface ) su = new ShapeAnalysis_Surface( thesurf );
 488 +      // gp_Pnt2d suval = su->ValueOfUV ( pnt, BRep_Tool::Tolerance( TopoDS::Face(fmap(surfi)) ) );
 489 +      // suval.Coord( u, v);
 490 +      // pnt = thesurf->Value( u, v );
 491 +
 492 +      Handle(ShapeAnalysis_Surface) proj;
 493 +      BRepTopAdaptor_FClass2d *cls;
 494 +      GetFaceTools(surfi, proj, cls);
 495 +
 496 +      gp_Pnt2d p2d = proj->ValueOfUV(pnt, Precision::Confusion());
 497 +      if (cls->Perform(p2d) == TopAbs_OUT)
 498 +      {
 499 +        return false;
 500 +      }
 501 +      pnt = proj->Value(p2d);
 502 +      p2d.Coord(u, v);
 503 +
 504        p = Point<3> (pnt.X(), pnt.Y(), pnt.Z());
 505
 506 +      return true;
 507     }
 508
 509
 510 @@ -979,54 +1021,69 @@
 511     {
 512        gp_Pnt p(ap(0), ap(1), ap(2));
 513
 514 -      Handle(Geom_Surface) surface = BRep_Tool::Surface(TopoDS::Face(fmap(surfi)));
 515 -
 516 -      gp_Pnt x = surface->Value (u,v);
 517 -
 518 -      if (p.SquareDistance(x) <= sqr(PROJECTION_TOLERANCE)) return true;
 519 -
 520 -      gp_Vec du, dv;
 521 -
 522 -      surface->D1(u,v,x,du,dv);
 523 -
 524 -      int count = 0;
 525 -
 526 -      gp_Pnt xold;
 527 -      gp_Vec n;
 528 -      double det, lambda, mu;
 529 -
 530 -      do {
 531 -         count++;
 532 -
 533 -         n = du^dv;
 534 -
 535 -         det = Det3 (n.X(), du.X(), dv.X(),
 536 -            n.Y(), du.Y(), dv.Y(),
 537 -            n.Z(), du.Z(), dv.Z());
 538 -
 539 -         if (det < 1e-15) return false;
 540 -
 541 -         lambda = Det3 (n.X(), p.X()-x.X(), dv.X(),
 542 -            n.Y(), p.Y()-x.Y(), dv.Y(),
 543 -            n.Z(), p.Z()-x.Z(), dv.Z())/det;
 544 -
 545 -         mu     = Det3 (n.X(), du.X(), p.X()-x.X(),
 546 -            n.Y(), du.Y(), p.Y()-x.Y(),
 547 -            n.Z(), du.Z(), p.Z()-x.Z())/det;
 548 -
 549 -         u += lambda;
 550 -         v += mu;
 551 -
 552 -         xold = x;
 553 -         surface->D1(u,v,x,du,dv);
 554 -
 555 -      } while (xold.SquareDistance(x) > sqr(PROJECTION_TOLERANCE) && count < 50);
 556 -
 557 -      //    (*testout) << "FastProject count: " << count << endl;
 558 -
 559 -      if (count == 50) return false;
 560 -
 561 -      ap = Point<3> (x.X(), x.Y(), x.Z());
 562 +      // -- Optimization: use cached projector and classifier
 563 +      // Handle(Geom_Surface) surface = BRep_Tool::Surface(TopoDS::Face(fmap(surfi)));
 564 +      //
 565 +      // gp_Pnt x = surface->Value (u,v);
 566 +      //
 567 +      // if (p.SquareDistance(x) <= sqr(PROJECTION_TOLERANCE)) return true;
 568 +      //
 569 +      // gp_Vec du, dv;
 570 +      //
 571 +      // surface->D1(u,v,x,du,dv);
 572 +      //
 573 +      // int count = 0;
 574 +      //
 575 +      // gp_Pnt xold;
 576 +      // gp_Vec n;
 577 +      // double det, lambda, mu;
 578 +      //
 579 +      // do {
 580 +      //    count++;
 581 +      //
 582 +      //    n = du^dv;
 583 +      //
 584 +      //    det = Det3 (n.X(), du.X(), dv.X(),
 585 +      //       n.Y(), du.Y(), dv.Y(),
 586 +      //       n.Z(), du.Z(), dv.Z());
 587 +      //
 588 +      //    if (det < 1e-15) return false;
 589 +      //
 590 +      //    lambda = Det3 (n.X(), p.X()-x.X(), dv.X(),
 591 +      //       n.Y(), p.Y()-x.Y(), dv.Y(),
 592 +      //       n.Z(), p.Z()-x.Z(), dv.Z())/det;
 593 +      //
 594 +      //    mu     = Det3 (n.X(), du.X(), p.X()-x.X(),
 595 +      //       n.Y(), du.Y(), p.Y()-x.Y(),
 596 +      //       n.Z(), du.Z(), p.Z()-x.Z())/det;
 597 +      //
 598 +      //    u += lambda;
 599 +      //    v += mu;
 600 +      //
 601 +      //    xold = x;
 602 +      //    surface->D1(u,v,x,du,dv);
 603 +      //
 604 +      // } while (xold.SquareDistance(x) > sqr(PROJECTION_TOLERANCE) && count < 50);
 605 +      //
 606 +      // //    (*testout) << "FastProject count: " << count << endl;
 607 +      //
 608 +      // if (count == 50) return false;
 609 +      //
 610 +      // ap = Point<3> (x.X(), x.Y(), x.Z());
 611 +      Handle(ShapeAnalysis_Surface) proj;
 612 +      BRepTopAdaptor_FClass2d *cls;
 613 +      GetFaceTools(surfi, proj, cls);
 614 +
 615 +      gp_Pnt2d p2d = proj->NextValueOfUV(gp_Pnt2d(u,v), p, Precision::Confusion());
 616 +      if (cls->Perform(p2d) == TopAbs_OUT)
 617 +      {
 618 +        //cout << "Projection fails" << endl;
 619 +        return false;
 620 +      }
 621 +
 622 +      p = proj->Value(p2d);
 623 +      p2d.Coord(u, v);
 624 +      ap = Point<3> (p.X(), p.Y(), p.Z());
 625
 626        return true;
 627     }
 628 diff -Naur netgen-4.9.13_orig/libsrc/occ/occgeom.hpp netgen-4.9.13_new/libsrc/occ/occgeom.hpp
 629 --- netgen-4.9.13_orig/libsrc/occ/occgeom.hpp   2010-01-14 19:56:19.000000000 +0300
 630 +++ netgen-4.9.13_new/libsrc/occ/occgeom.hpp    2011-06-28 11:07:52.000000000 +0400
 631 @@ -15,8 +15,8 @@
 632  #include "Geom_Curve.hxx"
 633  #include "Geom2d_Curve.hxx"
 634  #include "Geom_Surface.hxx"
 635 -#include "GeomAPI_ProjectPointOnSurf.hxx"
 636 -#include "GeomAPI_ProjectPointOnCurve.hxx"
 637 +// #include "GeomAPI_ProjectPointOnSurf.hxx"
 638 +// #include "GeomAPI_ProjectPointOnCurve.hxx"
 639  #include "BRepTools.hxx"
 640  #include "TopExp.hxx"
 641  #include "BRepBuilderAPI_MakeVertex.hxx"
 642 @@ -42,8 +42,8 @@
 643  #include "Geom_Curve.hxx"
 644  #include "Geom2d_Curve.hxx"
 645  #include "Geom_Surface.hxx"
 646 -#include "GeomAPI_ProjectPointOnSurf.hxx"
 647 -#include "GeomAPI_ProjectPointOnCurve.hxx"
 648 +// #include "GeomAPI_ProjectPointOnSurf.hxx"
 649 +// #include "GeomAPI_ProjectPointOnCurve.hxx"
 650  #include "TopoDS_Wire.hxx"
 651  #include "BRepTools_WireExplorer.hxx"
 652  #include "BRepTools.hxx"
 653 @@ -68,7 +68,7 @@
 654  #include "IGESToBRep_Reader.hxx"
 655  #include "Interface_Static.hxx"
 656  #include "GeomAPI_ExtremaCurveCurve.hxx"
 657 -#include "Standard_ErrorHandler.hxx"
 658 +//#include "Standard_ErrorHandler.hxx"
 659  #include "Standard_Failure.hxx"
 660  #include "ShapeUpgrade_ShellSewing.hxx"
 661  #include "ShapeFix_Shape.hxx"
 662 @@ -80,6 +80,10 @@
 663  #include "ShapeAnalysis.hxx"
 664  #include "ShapeBuild_ReShape.hxx"
 665
 666 +// -- Optimization: to use cached projector and classifier
 667 +#include <NCollection_DataMap.hxx>
 668 +class Handle_ShapeAnalysis_Surface;
 669 +class BRepTopAdaptor_FClass2d;
 670
 671  // Philippose - 29/01/2009
 672  // OpenCascade XDE Support
 673 @@ -190,6 +194,9 @@
 674     class OCCGeometry : public NetgenGeometry
 675     {
 676        Point<3> center;
 677 +      // -- Optimization: to use cached projector and classifier
 678 +      mutable NCollection_DataMap<int,Handle_ShapeAnalysis_Surface> fprjmap;
 679 +      mutable NCollection_DataMap<int,BRepTopAdaptor_FClass2d*> fclsmap;
 680
 681     public:
 682        TopoDS_Shape shape;
 683 @@ -241,6 +248,8 @@
 684           vmap.Clear();
 685        }
 686
 687 +      ~OCCGeometry();      // -- to free cached projector and classifier
 688 +
 689        void BuildFMap();
 690
 691        Box<3> GetBoundingBox()
 692 @@ -260,9 +269,14 @@
 693        Point<3> Center()
 694        {  return center;}
 695
 696 -      void Project (int surfi, Point<3> & p) const;
 697 +      // void Project (int surfi, Point<3> & p) const; -- optimization
 698 +      bool Project (int surfi, Point<3> & p, double& u, double& v) const;
 699        bool FastProject (int surfi, Point<3> & ap, double& u, double& v) const;
 700
 701 +      // -- Optimization: to use cached projector and classifier
 702 +      void GetFaceTools(int surfi, Handle(ShapeAnalysis_Surface)& proj,
 703 +                        BRepTopAdaptor_FClass2d*& cls) const;
 704 +
 705        OCCSurface GetSurface (int surfi)
 706        {
 707           cout << "OCCGeometry::GetSurface using PLANESPACE" << endl;
 708 diff -Naur netgen-4.9.13_orig/libsrc/occ/occmeshsurf.cpp netgen-4.9.13_new/libsrc/occ/occmeshsurf.cpp
 709 --- netgen-4.9.13_orig/libsrc/occ/occmeshsurf.cpp       2009-08-24 06:32:47.000000000 +0400
 710 +++ netgen-4.9.13_new/libsrc/occ/occmeshsurf.cpp        2011-06-03 19:55:14.000000000 +0400
 711 @@ -6,6 +6,7 @@
 712  #include <meshing.hpp>
 713  #include <GeomLProp_SLProps.hxx>
 714  #include <ShapeAnalysis_Surface.hxx>
 715 +#include <GeomAPI_ProjectPointOnCurve.hxx> // -- moved here from occgeom.hpp
 716
 717
 718  namespace netgen
 719 @@ -434,23 +435,33 @@
 720
 721    void MeshOptimize2dOCCSurfaces :: ProjectPoint (INDEX surfind, Point<3> & p) const
 722    {
 723 -    geometry.Project (surfind, p);
 724 +    // geometry.Project (surfind, p); -- signature of Project() changed for optimization
 725 +    double u, v;
 726 +    geometry.Project (surfind, p, u, v);
 727    }
 728
 729
 730    int MeshOptimize2dOCCSurfaces :: ProjectPointGI (INDEX surfind, Point<3> & p, PointGeomInfo & gi) const
 731    {
 732 -    double u = gi.u;
 733 -    double v = gi.v;
 734 +    //double u = gi.u;
 735 +    //double v = gi.v;
 736
 737      Point<3> hp = p;
 738 -    if (geometry.FastProject (surfind, hp, u, v))
 739 -      {
 740 -       p = hp;
 741 -       return 1;
 742 -      }
 743 -    ProjectPoint (surfind, p);
 744 -    return CalcPointGeomInfo (surfind, gi, p);
 745 +    // -- u and v are computed by FastProject() and Project(), no need to call CalcPointGeomInfo()
 746 +    // if (geometry.FastProject (surfind, hp, u, v))
 747 +    //   {
 748 +    //    p = hp;
 749 +    //    return 1;
 750 +    //   }
 751 +    // ProjectPoint (surfind, p);
 752 +    // return CalcPointGeomInfo (surfind, gi, p);
 753 +    bool ok;
 754 +    if (gi.trignum > 0)
 755 +      ok = geometry.FastProject (surfind, hp, gi.u, gi.v);
 756 +    else
 757 +      ok = geometry.Project (surfind, hp, gi.u, gi.v);
 758 +    p = hp;
 759 +    return ok;
 760    }
 761
 762
 763 @@ -680,7 +691,8 @@
 764         if (!geometry.FastProject (surfi, hnewp, u, v))
 765           {
 766           //  cout << "Fast projection to surface fails! Using OCC projection" << endl;
 767 -           geometry.Project (surfi, hnewp);
 768 +           // geometry.Project (surfi, hnewp); -- Project() changed for optimization
 769 +           geometry.Project (surfi, hnewp, u, v);
 770           }
 771
 772         newgi.trignum = 1;
 773 @@ -689,7 +701,7 @@
 774        }
 775
 776      newp = hnewp;
 777 -  }
 778 +  }//; -- to compile with -Wall -pedantic
 779
 780
 781    void OCCRefinementSurfaces ::
 782 @@ -708,14 +720,18 @@
 783      hnewp = Point<3> (pnt.X(), pnt.Y(), pnt.Z());
 784      newp = hnewp;
 785      newgi = ap1;
 786 -  };
 787 +  }
 788
 789
 790    void OCCRefinementSurfaces :: ProjectToSurface (Point<3> & p, int surfi)
 791    {
 792      if (surfi > 0)
 793 -      geometry.Project (surfi, p);
 794 -  };
 795 +      // geometry.Project (surfi, p); -- Project() changed for optimization
 796 +    {
 797 +      double u, v;
 798 +      geometry.Project (surfi, p, u, v);
 799 +    }
 800 +  }//; -- to compile with -Wall -pedantic
 801
 802    void OCCRefinementSurfaces :: ProjectToSurface (Point<3> & p, int surfi, PointGeomInfo & gi)
 803    {
 804 @@ -723,9 +739,10 @@
 805        if (!geometry.FastProject (surfi, p, gi.u, gi.v))
 806         {
 807           cout << "Fast projection to surface fails! Using OCC projection" << endl;
 808 -         geometry.Project (surfi, p);
 809 +          double u, v;
 810 +         geometry.Project (surfi, p, u, v);
 811         }
 812 -  };
 813 +  }
 814
 815
 816
 817 diff -Naur netgen-4.9.13_orig/libsrc/occ/utilities.h netgen-4.9.13_new/libsrc/occ/utilities.h
 818 --- netgen-4.9.13_orig/libsrc/occ/utilities.h   2009-08-24 06:12:24.000000000 +0400
 819 +++ netgen-4.9.13_new/libsrc/occ/utilities.h    2011-04-20 13:20:56.000000000 +0400
 820 @@ -33,6 +33,7 @@
 821
 822  #include <string>
 823  #include <iostream>
 824 +#include <iomanip>
 825  #include <cstdlib>
 826  // #include "SALOME_Log.hxx"
 827
 828 diff -Naur netgen-4.9.13_orig/libsrc/stlgeom/stlgeommesh.cpp netgen-4.9.13_new/libsrc/stlgeom/stlgeommesh.cpp
 829 --- netgen-4.9.13_orig/libsrc/stlgeom/stlgeommesh.cpp   2009-08-10 15:40:51.000000000 +0400
 830 +++ netgen-4.9.13_new/libsrc/stlgeom/stlgeommesh.cpp    2011-06-03 19:57:16.000000000 +0400
 831 @@ -1435,7 +1435,8 @@
 832
 833           if (!optstring || strlen(optstring) == 0)
 834             {
 835 -             mparam.optimize2d = "smcm";
 836 +             //mparam.optimize2d = (char*)"smcm";
 837 +              mparam.optimize2d = (char*)"smcm";
 838             }
 839           else
 840             {
 841 @@ -1451,7 +1452,8 @@
 842                                  mparam.grading);
 843               mesh -> LoadLocalMeshSize (mparam.meshsizefilename);
 844               mesh -> CalcLocalHFromSurfaceCurvature (stlparam.resthsurfmeshcurvfac);
 845 -             mparam.optimize2d = "cmsmSm";
 846 +             //mparam.optimize2d = (char*)"cmsmSm";
 847 +              mparam.optimize2d = (char*)"cmsmSm";
 848               STLSurfaceOptimization (*stlgeometry, *mesh, mparam);
 849  #ifdef STAT_STREAM
 850               (*statout) << GetTime() << " & ";
 851 @@ -1557,7 +1559,8 @@
 852
 853           if (!optstring || strlen(optstring) == 0)
 854             {
 855 -             mparam.optimize3d = "cmdmstm";
 856 +              //mparam.optimize3d = "cmdmstm";
 857 +             mparam.optimize3d = (char*)"cmdmstm";
 858             }
 859           else
 860             {
 861 --- netgen-4.9.13_orig/nglib/nglib.h    2010-05-18 15:20:25.000000000 +0400
 862 +++ netgen-4.9.13_new/nglib/nglib.h     2010-05-31 13:02:19.000000000 +0400
 863 @@ -24,7 +24,7 @@
 864  // Philippose - 14.02.2009
 865  // Modifications for creating a DLL in Windows
 866  #ifdef WIN32
 867 -   #ifdef NGLIB_EXPORTS || nglib_EXPORTS
 868 +   #if defined NGLIB_EXPORTS || defined nglib_EXPORTS
 869        #define DLL_HEADER   __declspec(dllexport)
 870     #else
 871        #define DLL_HEADER   __declspec(dllimport)