src/NETGEN/netgen50ForSalome.patch

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