config_files/patches/netgen-4.9.13-for-SALOME.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       2013-05-31 11:28:04.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       2013-05-31 11:28:04.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     2013-05-31 11:28:04.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/Partition_Inter3d.cxx netgen-4.9.13_new/libsrc/occ/Partition_Inter3d.cxx
 166 --- netgen-4.9.13_orig/libsrc/occ/Partition_Inter3d.cxx 2009-08-24 06:12:24.000000000 +0400
 167 +++ netgen-4.9.13_new/libsrc/occ/Partition_Inter3d.cxx  2013-05-31 11:28:04.000000000 +0400
 168 @@ -86,6 +86,9 @@
 169  #include <TopOpeBRepTool_OutCurveType.hxx>
 170  #include <TopOpeBRep_DSFiller.hxx>
 171  #include <TopTools_DataMapIteratorOfDataMapOfShapeShape.hxx>
 172 +
 173 +#include <Standard_Version.hxx>
 174 +
 175  #include <stdio.h>
 176
 177  //=======================================================================
 178 @@ -243,7 +246,12 @@
 179        Standard_Integer i, nbExt = anExtPS.NbExt();
 180        Extrema_POnSurf aPOnSurf;
 181        for (i = 1; i <= nbExt; ++i )
 182 +#if (OCC_VERSION_MAJOR << 16 | OCC_VERSION_MINOR << 8 | OCC_VERSION_MAINTENANCE) > 0x060400
 183 +// porting to OCCT6.5.1
 184 +        if (anExtPS.SquareDistance( i ) <= TolE * TolE) {
 185 +#else
 186          if (anExtPS.Value( i ) <= TolE) {
 187 +#endif
 188            aPOnSurf = anExtPS.Point( i );
 189            break;
 190          }
 191 diff -Naur netgen-4.9.13_orig/libsrc/occ/Partition_Loop.cxx netgen-4.9.13_new/libsrc/occ/Partition_Loop.cxx
 192 --- netgen-4.9.13_orig/libsrc/occ/Partition_Loop.cxx    2009-08-24 06:12:24.000000000 +0400
 193 +++ netgen-4.9.13_new/libsrc/occ/Partition_Loop.cxx     2013-05-31 11:28:04.000000000 +0400
 194 @@ -178,7 +178,7 @@
 195        }
 196      }
 197
 198 -    Standard_Real anglemax = - PI;
 199 +    Standard_Real anglemax = - M_PI;
 200      TopoDS_Edge   SelectedEdge;
 201      for ( itl.Initialize(LE); itl.More(); itl.Next()) {
 202        const TopoDS_Edge& E = TopoDS::Edge(itl.Value());
 203 diff -Naur netgen-4.9.13_orig/libsrc/occ/Partition_Loop2d.cxx netgen-4.9.13_new/libsrc/occ/Partition_Loop2d.cxx
 204 --- netgen-4.9.13_orig/libsrc/occ/Partition_Loop2d.cxx  2009-08-24 06:12:24.000000000 +0400
 205 +++ netgen-4.9.13_new/libsrc/occ/Partition_Loop2d.cxx   2013-05-31 11:28:04.000000000 +0400
 206 @@ -22,7 +22,6 @@
 207  #include <BRepAdaptor_Surface.hxx>
 208  #include <BRepAlgo_AsDes.hxx>
 209  #include <BRepAlgo_FaceRestrictor.hxx>
 210 -#include <BRepOffset_DataMapOfShapeReal.hxx>
 211  #include <BRepTopAdaptor_FClass2d.hxx>
 212  #include <BRep_Builder.hxx>
 213  #include <BRep_Tool.hxx>
 214 @@ -51,6 +50,15 @@
 215  #include <gp_Pnt.hxx>
 216  #include <gp_Pnt2d.hxx>
 217
 218 +#include <Standard_Version.hxx>
 219 +
 220 +#if (OCC_VERSION_MAJOR << 16 | OCC_VERSION_MINOR << 8 | OCC_VERSION_MAINTENANCE) > 0x060400
 221 +// porting to OCCT6.5.1
 222 +#include <TopTools_DataMapOfShapeReal.hxx>
 223 +#else
 224 +#include <BRepOffset_DataMapOfShapeReal.hxx>
 225 +#endif
 226 +
 227  //=======================================================================
 228  //function : Partition_Loop2d
 229  //purpose  :
 230 @@ -209,7 +217,7 @@
 231      Cc->D1(uc, PC, CTg1);
 232      if (!isForward) CTg1.Reverse();
 233
 234 -    Standard_Real anglemin = 3 * PI, tolAng = 1.e-8;
 235 +    Standard_Real anglemin = 3 * M_PI, tolAng = 1.e-8;
 236
 237      // select an edge whose first derivative is most left of CTg1
 238      // ie an angle between Tg1 and CTg1 is least
 239 @@ -233,7 +241,7 @@
 240        // -PI < angle < PI
 241        Standard_Real angle = Tg1.Angle(CTg1);
 242
 243 -      if (PI - Abs(angle) <= tolAng)
 244 +      if (M_PI - Abs(angle) <= tolAng)
 245        {
 246          // an angle is too close to PI; assure that an angle sign really
 247          // reflects an edge position: +PI - an edge is worst,
 248 @@ -519,7 +527,12 @@
 249      DC.Initialize( DegEdge, F );
 250
 251    // avoid intersecting twice the same edge
 252 +#if (OCC_VERSION_MAJOR << 16 | OCC_VERSION_MINOR << 8 | OCC_VERSION_MAINTENANCE) > 0x060400
 253 +// porting to OCCT6.5.1
 254 +  TopTools_DataMapOfShapeReal EUMap ( EdgesList.Extent() );
 255 +#else
 256    BRepOffset_DataMapOfShapeReal EUMap ( EdgesList.Extent() );
 257 +#endif
 258
 259    Standard_Real U, f, l;
 260    BRep_Tool::Range (DegEdge, f, l);
 261 diff -Naur netgen-4.9.13_orig/libsrc/occ/Partition_Spliter.cxx netgen-4.9.13_new/libsrc/occ/Partition_Spliter.cxx
 262 --- netgen-4.9.13_orig/libsrc/occ/Partition_Spliter.cxx 2009-08-24 06:12:24.000000000 +0400
 263 +++ netgen-4.9.13_new/libsrc/occ/Partition_Spliter.cxx  2013-05-31 11:28:04.000000000 +0400
 264 @@ -79,6 +79,8 @@
 265  #include <GeomAdaptor_Curve.hxx>
 266  #include <TopOpeBRepTool_CurveTool.hxx>
 267
 268 +#include <Standard_Version.hxx>
 269 +
 270  #ifdef DEB
 271  //# define PART_PERF
 272  #endif
 273 @@ -1169,7 +1171,12 @@
 274            for (; j<=nbj && ok; ++j) {
 275              if (Extrema.IsMin(j)) {
 276               hasMin = Standard_True;
 277 +#if (OCC_VERSION_MAJOR << 16 | OCC_VERSION_MINOR << 8 | OCC_VERSION_MAINTENANCE) > 0x060400
 278 +// porting to OCCT6.5.1
 279 +              ok = Extrema.SquareDistance(j) <= tol * tol;
 280 +#else
 281                ok = Extrema.Value(j) <= tol;
 282 +#endif
 283             }
 284            }
 285          }
 286 diff -Naur netgen-4.9.13_orig/libsrc/occ/occconstruction.cpp netgen-4.9.13_new/libsrc/occ/occconstruction.cpp
 287 --- netgen-4.9.13_orig/libsrc/occ/occconstruction.cpp   2009-08-24 06:32:47.000000000 +0400
 288 +++ netgen-4.9.13_new/libsrc/occ/occconstruction.cpp    2013-05-31 11:28:04.000000000 +0400
 289 @@ -28,7 +28,7 @@
 290  #include <BRepAlgoAPI_Common.hxx>
 291  #include <BRepAlgoAPI_Fuse.hxx>
 292  #include <BRepAlgoAPI_Section.hxx>
 293 -#include <BRepOffsetAPI_Sewing.hxx>
 294 +//#include <BRepOffsetAPI_Sewing.hxx>
 295  //#include <BRepAlgo_Sewing.hxx>
 296  #include <BRepOffsetAPI_MakeOffsetShape.hxx>
 297  #include <ShapeFix_Shape.hxx>
 298 diff -Naur netgen-4.9.13_orig/libsrc/occ/occgenmesh.cpp netgen-4.9.13_new/libsrc/occ/occgenmesh.cpp
 299 --- netgen-4.9.13_orig/libsrc/occ/occgenmesh.cpp        2010-03-16 09:30:07.000000000 +0300
 300 +++ netgen-4.9.13_new/libsrc/occ/occgenmesh.cpp 2013-05-31 11:28:04.000000000 +0400
 301 @@ -15,6 +15,8 @@
 302
 303  #define DIVIDEEDGESECTIONS 1000
 304  #define IGNORECURVELENGTH 1e-4
 305 +// a small value used to avoid FPE
 306 +#define VSMALL 1e-10
 307
 308
 309     bool merge_solids = 1;
 310 @@ -26,7 +28,8 @@
 311        double nq = n*q;
 312
 313        Point<3> p = p0 + 0.5*n;
 314 -      double lambda = (p-l.p0)*n / nq;
 315 +      // double lambda = (p-l.p0)*n / nq;  -- avoid FPE
 316 +      double lambda = (fabs(nq) > 1e-10) ? (p-l.p0)*n / nq : -1;
 317
 318        if (lambda >= 0 && lambda <= 1)
 319        {
 320 @@ -54,6 +57,8 @@
 321
 322
 323
 324 +
 325 +   static // useless out of this file
 326     double ComputeH (double kappa)
 327     {
 328        double hret;
 329 @@ -62,7 +67,8 @@
 330        if (mparam.maxh * kappa < 1)
 331           hret = mparam.maxh;
 332        else
 333 -         hret = 1 / kappa;
 334 +        // hret = 1 / kappa; -- avoid FPE
 335 +        hret = 1 / (kappa + VSMALL);
 336
 337        if (mparam.maxh < hret)
 338           hret = mparam.maxh;
 339 @@ -71,8 +77,7 @@
 340     }
 341
 342
 343 -
 344 -
 345 +   static // useless out of this file
 346     void RestrictHTriangle (gp_Pnt2d & par0, gp_Pnt2d & par1, gp_Pnt2d & par2,
 347                             BRepLProp_SLProps * prop, Mesh & mesh, int depth, double h = 0)
 348     {
 349 @@ -168,8 +173,8 @@
 350           if(h < 1e-4*maxside)
 351              return;
 352
 353 -
 354 -         if (h > 30) return;
 355 +         // commented to restrict H on a large sphere for example
 356 +         //if (h > 30) return;
 357        }
 358
 359        if (h < maxside && depth < 10)
 360 @@ -228,6 +233,7 @@
 361
 362
 363
 364 +   static // useless out of this file
 365     void DivideEdge (TopoDS_Edge & edge, Array<MeshPoint> & ps,
 366                      Array<double> & params, Mesh & mesh)
 367     {
 368 @@ -247,8 +253,8 @@
 369        hvalue[0] = 0;
 370        pnt = c->Value(s0);
 371
 372 -      double olddist = 0;
 373 -      double dist = 0;
 374 +      //double olddist = 0; -- useless variables
 375 +      //double dist = 0;
 376
 377        int tmpVal = (int)(DIVIDEEDGESECTIONS);
 378
 379 @@ -256,15 +262,19 @@
 380        {
 381           oldpnt = pnt;
 382           pnt = c->Value(s0+(i/double(DIVIDEEDGESECTIONS))*(s1-s0));
 383 +         // -- no more than 1 segment per <edge length>/DIVIDEEDGESECTIONS
 384           hvalue[i] = hvalue[i-1] +
 385 -            1.0/mesh.GetH(Point3d(pnt.X(), pnt.Y(), pnt.Z()))*
 386 -            pnt.Distance(oldpnt);
 387 +         //   1.0/mesh.GetH(Point3d(pnt.X(), pnt.Y(), pnt.Z()))*
 388 +         //   pnt.Distance(oldpnt);
 389 +           min( 1.0,
 390 +                1.0/mesh.GetH(Point3d(pnt.X(), pnt.Y(), pnt.Z()))*
 391 +                pnt.Distance(oldpnt));
 392
 393           //(*testout) << "mesh.GetH(Point3d(pnt.X(), pnt.Y(), pnt.Z())) " << mesh.GetH(Point3d(pnt.X(), pnt.Y(), pnt.Z()))
 394           //       <<  " pnt.Distance(oldpnt) " << pnt.Distance(oldpnt) << endl;
 395
 396 -         olddist = dist;
 397 -         dist = pnt.Distance(oldpnt);
 398 +         //olddist = dist; -- useless variables
 399 +         //dist = pnt.Distance(oldpnt);
 400        }
 401
 402        //  nsubedges = int(ceil(hvalue[DIVIDEEDGESECTIONS]));
 403 @@ -279,7 +289,10 @@
 404        {
 405           if (hvalue[i1]/hvalue[DIVIDEEDGESECTIONS]*nsubedges >= i)
 406           {
 407 -            params[i] = s0+(i1/double(DIVIDEEDGESECTIONS))*(s1-s0);
 408 +            // -- for nsubedges comparable to DIVIDEEDGESECTIONS
 409 +            //params[i] = s0+(i1/double(DIVIDEEDGESECTIONS))*(s1-s0);
 410 +            double d1 = i1 - (hvalue[i1] - i*hvalue[DIVIDEEDGESECTIONS]/nsubedges)/(hvalue[i1]-hvalue[i1-1]);
 411 +            params[i] = s0+(d1/double(DIVIDEEDGESECTIONS))*(s1-s0);
 412              pnt = c->Value(params[i]);
 413              ps[i-1] = MeshPoint (Point3d(pnt.X(), pnt.Y(), pnt.Z()));
 414              i++;
 415 @@ -323,6 +336,7 @@
 416        (*testout) << "nedges = " << nedges << endl;
 417
 418        double eps = 1e-6 * geom.GetBoundingBox().Diam();
 419 +      const double eps2 = eps * eps; // -- small optimization
 420
 421        for (int i = 1; i <= nvertices; i++)
 422        {
 423 @@ -332,7 +346,8 @@
 424           bool exists = 0;
 425           if (merge_solids)
 426              for (PointIndex pi = 1; pi <= mesh.GetNP(); pi++)
 427 -               if ( Dist2 (mesh[pi], Point<3>(mp)) < eps*eps)
 428 +               //if ( Dist2 (mesh[pi], Point<3>(mp)) < eps*eps)
 429 +               if ( Dist2 (mesh[pi], Point<3>(mp)) < eps2 ) // -- small optimization
 430                 {
 431                    exists = 1;
 432                    break;
 433 @@ -362,6 +377,7 @@
 434           {
 435              TopoDS_Face face = TopoDS::Face(exp1.Current());
 436              int facenr = geom.fmap.FindIndex(face);
 437 +            if ( facenr < 1 ) continue; // -- to support SALOME sub-meshes
 438
 439              if (face2solid[0][facenr-1] == 0)
 440                 face2solid[0][facenr-1] = solidnr;
 441 @@ -381,6 +397,7 @@
 442        int facenr = 0;
 443        int edgenr = 0;
 444
 445 +      edgenr = mesh.GetNSeg(); // to support SALOME sub-meshes
 446
 447        (*testout) << "faces = " << geom.fmap.Extent() << endl;
 448        int curr = 0;
 449 @@ -442,6 +459,7 @@
 450                    //(*testout) << "ignoring degenerated edge" << endl;
 451                    continue;
 452                 }
 453 +               if ( geom.emap.FindIndex(edge) < 1 ) continue; // to support SALOME sub-meshes
 454
 455                 if (geom.vmap.FindIndex(TopExp::FirstVertex (edge)) ==
 456                    geom.vmap.FindIndex(TopExp::LastVertex (edge)))
 457 @@ -479,15 +497,64 @@
 458                 }
 459                 else
 460                 {
 461 -                  Point<3> fp = occ2ng (BRep_Tool::Pnt (TopExp::FirstVertex (edge)));
 462 -                  Point<3> lp = occ2ng (BRep_Tool::Pnt (TopExp::LastVertex (edge)));
 463 +                 TopoDS_Iterator vIt( edge, false );
 464 +                 TopoDS_Vertex v1 = TopoDS::Vertex( vIt.Value() ); vIt.Next();
 465 +                 TopoDS_Vertex v2 = TopoDS::Vertex( vIt.Value() );
 466 +                 if ( v1.Orientation() == TopAbs_REVERSED )
 467 +                   std::swap( v1, v2 );
 468 +                 const bool isClosedEdge = v1.IsSame( v2 );
 469 +
 470 +                  Point<3> fp = occ2ng (BRep_Tool::Pnt (v1));
 471 +                  Point<3> lp = occ2ng (BRep_Tool::Pnt (v2));
 472 +                  double tol2 = std::min( eps*eps, 1e-6 * Dist2( fp, lp ));
 473 +                  if ( isClosedEdge )
 474 +                    tol2 = BRep_Tool::Tolerance( v1 ) * BRep_Tool::Tolerance( v1 );
 475
 476                    pnums[0] = -1;
 477                    pnums.Last() = -1;
 478                    for (PointIndex pi = 1; pi < first_ep; pi++)
 479                    {
 480 -                     if (Dist2 (mesh[pi], fp) < eps*eps) pnums[0] = pi;
 481 -                     if (Dist2 (mesh[pi], lp) < eps*eps) pnums.Last() = pi;
 482 +                    if (Dist2 (mesh[pi], fp) < tol2) pnums[0] = pi;
 483 +                    if (Dist2 (mesh[pi], lp) < tol2) pnums.Last() = pi;
 484 +                  }
 485 +                  if (( isClosedEdge && pnums[0] != pnums.Last() ) ||
 486 +                      ( !isClosedEdge && pnums[0] == pnums.Last() ))
 487 +                    pnums[0] = pnums.Last() = -1;
 488 +                  if ( pnums[0] == -1 || pnums.Last() == -1 )
 489 +                  {
 490 +                    // take into account a possible large gap between a vertex and an edge curve
 491 +                    // end and a large vertex tolerance covering the whole edge
 492 +                    if ( pnums[0] == -1 )
 493 +                    {
 494 +                      double tol = BRep_Tool::Tolerance( v1 );
 495 +                      for (PointIndex pi = 1; pi < first_ep; pi++)
 496 +                        if (pi != pnums.Last() && Dist2 (mesh[pi], fp) < 2*tol*tol)
 497 +                          pnums[0] = pi;
 498 +
 499 +                      if ( pnums[0] == -1 )
 500 +                        pnums[0] = first_ep-1- nvertices + geom.vmap.FindIndex ( v1 );
 501 +                    }
 502 +                    if ( isClosedEdge )
 503 +                    {
 504 +                      pnums.Last() = pnums[0];
 505 +                    }
 506 +                    else
 507 +                    {
 508 +                      if ( pnums.Last() == -1 )
 509 +                      {
 510 +                        double tol = BRep_Tool::Tolerance( v2 );
 511 +                        for (PointIndex pi = 1; pi < first_ep; pi++)
 512 +                          if (pi != pnums[0] && Dist2 (mesh[pi], lp) < 2*tol*tol)
 513 +                            pnums.Last() = pi;
 514 +
 515 +                        if ( pnums.Last() == -1 )
 516 +                          pnums.Last() = first_ep-1-nvertices + geom.vmap.FindIndex ( v2 );
 517 +                      }
 518 +
 519 +                      if ( Dist2( fp, mesh[PointIndex(pnums[0])]) >
 520 +                           Dist2( lp, mesh[PointIndex(pnums.Last())]))
 521 +                      std::swap( pnums[0], pnums.Last() );
 522 +                    }
 523                    }
 524                 }
 525
 526 @@ -633,7 +700,8 @@
 527           }
 528
 529           (*testout) << "mesh face " << k << endl;
 530 -         multithread.percent = 100 * k / (mesh.GetNFD()+1e-10);
 531 +         // multithread.percent = 100 * k / (mesh.GetNFD()+1e-10); -- avoid FPE
 532 +         multithread.percent = 100 * k / (mesh.GetNFD() + VSMALL);
 533           geom.facemeshstatus[k-1] = -1;
 534
 535
 536 @@ -901,7 +969,8 @@
 537           //      if (k != 36) continue;
 538
 539           //      (*testout) << "optimize face " << k << endl;
 540 -         multithread.percent = 100 * k / (mesh.GetNFD()+1e-10);
 541 +         //multithread.percent = 100 * k / (mesh.GetNFD()+1e-10); -- avoid FPE
 542 +         multithread.percent = 100 * k / (mesh.GetNFD() + VSMALL);
 543
 544           FaceDescriptor & fd = mesh.GetFaceDescriptor(k);
 545
 546 @@ -1456,3 +1525,4 @@
 547  }
 548
 549  #endif
 550 +
 551 diff -Naur netgen-4.9.13_orig/libsrc/occ/occgeom.cpp netgen-4.9.13_new/libsrc/occ/occgeom.cpp
 552 --- netgen-4.9.13_orig/libsrc/occ/occgeom.cpp   2010-03-05 16:16:21.000000000 +0300
 553 +++ netgen-4.9.13_new/libsrc/occ/occgeom.cpp    2013-05-31 11:28:04.000000000 +0400
 554 @@ -8,6 +8,8 @@
 555  #include "ShapeAnalysis_CheckSmallFace.hxx"
 556  #include "ShapeAnalysis_DataMapOfShapeListOfReal.hxx"
 557  #include "ShapeAnalysis_Surface.hxx"
 558 +#include <BRepTopAdaptor_FClass2d.hxx> // -- to optimize Project() and FastProject()
 559 +#include <TopAbs_State.hxx>
 560  #include "BRepAlgoAPI_Fuse.hxx"
 561  #include "BRepCheck_Analyzer.hxx"
 562  #include "BRepLib.hxx"
 563 @@ -16,10 +18,17 @@
 564  #include "ShapeFix_FixSmallFace.hxx"
 565  #include "Partition_Spliter.hxx"
 566
 567 -
 568  namespace netgen
 569  {
 570 -   void OCCGeometry :: PrintNrShapes ()
 571 +  // free data used to optimize Project() and FastProject()
 572 +  OCCGeometry::~OCCGeometry()
 573 +  {
 574 +    NCollection_DataMap<int,BRepTopAdaptor_FClass2d*>::Iterator it(fclsmap);
 575 +    for (; it.More(); it.Next())
 576 +      delete it.Value();
 577 +  }
 578 +
 579 +  void OCCGeometry :: PrintNrShapes ()
 580     {
 581        TopExp_Explorer e;
 582        int count = 0;
 583 @@ -951,25 +960,58 @@
 584     }
 585
 586
 587 +   // returns a projector and a classifier for the given surface
 588 +   void OCCGeometry::GetFaceTools(int surfi, Handle(ShapeAnalysis_Surface)& proj,
 589 +                                  BRepTopAdaptor_FClass2d*& cls) const
 590 +   {
 591 +     //MSV: organize caching projector in the map
 592 +     if (fprjmap.IsBound(surfi))
 593 +     {
 594 +       proj = fprjmap.Find(surfi);
 595 +       cls = fclsmap.Find(surfi);
 596 +     }
 597 +     else
 598 +     {
 599 +       const TopoDS_Face& aFace = TopoDS::Face(fmap(surfi));
 600 +       Handle(Geom_Surface) aSurf = BRep_Tool::Surface(aFace);
 601 +       proj = new ShapeAnalysis_Surface(aSurf);
 602 +       fprjmap.Bind(surfi, proj);
 603 +       cls = new BRepTopAdaptor_FClass2d(aFace,Precision::Confusion());
 604 +       fclsmap.Bind(surfi, cls);
 605 +     }
 606 +   }
 607
 608 -
 609 -   void OCCGeometry :: Project (int surfi, Point<3> & p) const
 610 +   // void OCCGeometry :: Project (int surfi, Point<3> & p) const
 611 +   bool OCCGeometry :: Project (int surfi, Point<3> & p, double& u, double& v) const
 612     {
 613        static int cnt = 0;
 614        if (++cnt % 1000 == 0) cout << "Project cnt = " << cnt << endl;
 615
 616        gp_Pnt pnt(p(0), p(1), p(2));
 617
 618 -      double u,v;
 619 -      Handle( Geom_Surface ) thesurf = BRep_Tool::Surface(TopoDS::Face(fmap(surfi)));
 620 -      Handle( ShapeAnalysis_Surface ) su = new ShapeAnalysis_Surface( thesurf );
 621 -      gp_Pnt2d suval = su->ValueOfUV ( pnt, BRep_Tool::Tolerance( TopoDS::Face(fmap(surfi)) ) );
 622 -      suval.Coord( u, v);
 623 -      pnt = thesurf->Value( u, v );
 624 -
 625 -
 626 +      // -- Optimization: use cached projector and classifier
 627 +      // double u,v;
 628 +      // Handle( Geom_Surface ) thesurf = BRep_Tool::Surface(TopoDS::Face(fmap(surfi)));
 629 +      // Handle( ShapeAnalysis_Surface ) su = new ShapeAnalysis_Surface( thesurf );
 630 +      // gp_Pnt2d suval = su->ValueOfUV ( pnt, BRep_Tool::Tolerance( TopoDS::Face(fmap(surfi)) ) );
 631 +      // suval.Coord( u, v);
 632 +      // pnt = thesurf->Value( u, v );
 633 +
 634 +      Handle(ShapeAnalysis_Surface) proj;
 635 +      BRepTopAdaptor_FClass2d *cls;
 636 +      GetFaceTools(surfi, proj, cls);
 637 +
 638 +      gp_Pnt2d p2d = proj->ValueOfUV(pnt, Precision::Confusion());
 639 +      if (cls->Perform(p2d) == TopAbs_OUT)
 640 +      {
 641 +        return false;
 642 +      }
 643 +      pnt = proj->Value(p2d);
 644 +      p2d.Coord(u, v);
 645 +
 646        p = Point<3> (pnt.X(), pnt.Y(), pnt.Z());
 647
 648 +      return true;
 649     }
 650
 651
 652 @@ -979,54 +1021,69 @@
 653     {
 654        gp_Pnt p(ap(0), ap(1), ap(2));
 655
 656 -      Handle(Geom_Surface) surface = BRep_Tool::Surface(TopoDS::Face(fmap(surfi)));
 657 -
 658 -      gp_Pnt x = surface->Value (u,v);
 659 -
 660 -      if (p.SquareDistance(x) <= sqr(PROJECTION_TOLERANCE)) return true;
 661 -
 662 -      gp_Vec du, dv;
 663 -
 664 -      surface->D1(u,v,x,du,dv);
 665 -
 666 -      int count = 0;
 667 -
 668 -      gp_Pnt xold;
 669 -      gp_Vec n;
 670 -      double det, lambda, mu;
 671 -
 672 -      do {
 673 -         count++;
 674 -
 675 -         n = du^dv;
 676 -
 677 -         det = Det3 (n.X(), du.X(), dv.X(),
 678 -            n.Y(), du.Y(), dv.Y(),
 679 -            n.Z(), du.Z(), dv.Z());
 680 -
 681 -         if (det < 1e-15) return false;
 682 -
 683 -         lambda = Det3 (n.X(), p.X()-x.X(), dv.X(),
 684 -            n.Y(), p.Y()-x.Y(), dv.Y(),
 685 -            n.Z(), p.Z()-x.Z(), dv.Z())/det;
 686 -
 687 -         mu     = Det3 (n.X(), du.X(), p.X()-x.X(),
 688 -            n.Y(), du.Y(), p.Y()-x.Y(),
 689 -            n.Z(), du.Z(), p.Z()-x.Z())/det;
 690 -
 691 -         u += lambda;
 692 -         v += mu;
 693 -
 694 -         xold = x;
 695 -         surface->D1(u,v,x,du,dv);
 696 -
 697 -      } while (xold.SquareDistance(x) > sqr(PROJECTION_TOLERANCE) && count < 50);
 698 -
 699 -      //    (*testout) << "FastProject count: " << count << endl;
 700 -
 701 -      if (count == 50) return false;
 702 -
 703 -      ap = Point<3> (x.X(), x.Y(), x.Z());
 704 +      // -- Optimization: use cached projector and classifier
 705 +      // Handle(Geom_Surface) surface = BRep_Tool::Surface(TopoDS::Face(fmap(surfi)));
 706 +      //
 707 +      // gp_Pnt x = surface->Value (u,v);
 708 +      //
 709 +      // if (p.SquareDistance(x) <= sqr(PROJECTION_TOLERANCE)) return true;
 710 +      //
 711 +      // gp_Vec du, dv;
 712 +      //
 713 +      // surface->D1(u,v,x,du,dv);
 714 +      //
 715 +      // int count = 0;
 716 +      //
 717 +      // gp_Pnt xold;
 718 +      // gp_Vec n;
 719 +      // double det, lambda, mu;
 720 +      //
 721 +      // do {
 722 +      //    count++;
 723 +      //
 724 +      //    n = du^dv;
 725 +      //
 726 +      //    det = Det3 (n.X(), du.X(), dv.X(),
 727 +      //       n.Y(), du.Y(), dv.Y(),
 728 +      //       n.Z(), du.Z(), dv.Z());
 729 +      //
 730 +      //    if (det < 1e-15) return false;
 731 +      //
 732 +      //    lambda = Det3 (n.X(), p.X()-x.X(), dv.X(),
 733 +      //       n.Y(), p.Y()-x.Y(), dv.Y(),
 734 +      //       n.Z(), p.Z()-x.Z(), dv.Z())/det;
 735 +      //
 736 +      //    mu     = Det3 (n.X(), du.X(), p.X()-x.X(),
 737 +      //       n.Y(), du.Y(), p.Y()-x.Y(),
 738 +      //       n.Z(), du.Z(), p.Z()-x.Z())/det;
 739 +      //
 740 +      //    u += lambda;
 741 +      //    v += mu;
 742 +      //
 743 +      //    xold = x;
 744 +      //    surface->D1(u,v,x,du,dv);
 745 +      //
 746 +      // } while (xold.SquareDistance(x) > sqr(PROJECTION_TOLERANCE) && count < 50);
 747 +      //
 748 +      // //    (*testout) << "FastProject count: " << count << endl;
 749 +      //
 750 +      // if (count == 50) return false;
 751 +      //
 752 +      // ap = Point<3> (x.X(), x.Y(), x.Z());
 753 +      Handle(ShapeAnalysis_Surface) proj;
 754 +      BRepTopAdaptor_FClass2d *cls;
 755 +      GetFaceTools(surfi, proj, cls);
 756 +
 757 +      gp_Pnt2d p2d = proj->NextValueOfUV(gp_Pnt2d(u,v), p, Precision::Confusion());
 758 +      if (cls->Perform(p2d) == TopAbs_OUT)
 759 +      {
 760 +        //cout << "Projection fails" << endl;
 761 +        return false;
 762 +      }
 763 +
 764 +      p = proj->Value(p2d);
 765 +      p2d.Coord(u, v);
 766 +      ap = Point<3> (p.X(), p.Y(), p.Z());
 767
 768        return true;
 769     }
 770 diff -Naur netgen-4.9.13_orig/libsrc/occ/occgeom.hpp netgen-4.9.13_new/libsrc/occ/occgeom.hpp
 771 --- netgen-4.9.13_orig/libsrc/occ/occgeom.hpp   2010-01-14 19:56:19.000000000 +0300
 772 +++ netgen-4.9.13_new/libsrc/occ/occgeom.hpp    2013-05-31 11:28:04.000000000 +0400
 773 @@ -15,8 +15,8 @@
 774  #include "Geom_Curve.hxx"
 775  #include "Geom2d_Curve.hxx"
 776  #include "Geom_Surface.hxx"
 777 -#include "GeomAPI_ProjectPointOnSurf.hxx"
 778 -#include "GeomAPI_ProjectPointOnCurve.hxx"
 779 +// #include "GeomAPI_ProjectPointOnSurf.hxx"
 780 +// #include "GeomAPI_ProjectPointOnCurve.hxx"
 781  #include "BRepTools.hxx"
 782  #include "TopExp.hxx"
 783  #include "BRepBuilderAPI_MakeVertex.hxx"
 784 @@ -42,8 +42,8 @@
 785  #include "Geom_Curve.hxx"
 786  #include "Geom2d_Curve.hxx"
 787  #include "Geom_Surface.hxx"
 788 -#include "GeomAPI_ProjectPointOnSurf.hxx"
 789 -#include "GeomAPI_ProjectPointOnCurve.hxx"
 790 +// #include "GeomAPI_ProjectPointOnSurf.hxx"
 791 +// #include "GeomAPI_ProjectPointOnCurve.hxx"
 792  #include "TopoDS_Wire.hxx"
 793  #include "BRepTools_WireExplorer.hxx"
 794  #include "BRepTools.hxx"
 795 @@ -68,7 +68,7 @@
 796  #include "IGESToBRep_Reader.hxx"
 797  #include "Interface_Static.hxx"
 798  #include "GeomAPI_ExtremaCurveCurve.hxx"
 799 -#include "Standard_ErrorHandler.hxx"
 800 +//#include "Standard_ErrorHandler.hxx"
 801  #include "Standard_Failure.hxx"
 802  #include "ShapeUpgrade_ShellSewing.hxx"
 803  #include "ShapeFix_Shape.hxx"
 804 @@ -80,6 +80,10 @@
 805  #include "ShapeAnalysis.hxx"
 806  #include "ShapeBuild_ReShape.hxx"
 807
 808 +// -- Optimization: to use cached projector and classifier
 809 +#include <NCollection_DataMap.hxx>
 810 +class Handle_ShapeAnalysis_Surface;
 811 +class BRepTopAdaptor_FClass2d;
 812
 813  // Philippose - 29/01/2009
 814  // OpenCascade XDE Support
 815 @@ -190,6 +194,9 @@
 816     class OCCGeometry : public NetgenGeometry
 817     {
 818        Point<3> center;
 819 +      // -- Optimization: to use cached projector and classifier
 820 +      mutable NCollection_DataMap<int,Handle_ShapeAnalysis_Surface> fprjmap;
 821 +      mutable NCollection_DataMap<int,BRepTopAdaptor_FClass2d*> fclsmap;
 822
 823     public:
 824        TopoDS_Shape shape;
 825 @@ -241,6 +248,8 @@
 826           vmap.Clear();
 827        }
 828
 829 +      ~OCCGeometry();      // -- to free cached projector and classifier
 830 +
 831        void BuildFMap();
 832
 833        Box<3> GetBoundingBox()
 834 @@ -260,9 +269,14 @@
 835        Point<3> Center()
 836        {  return center;}
 837
 838 -      void Project (int surfi, Point<3> & p) const;
 839 +      // void Project (int surfi, Point<3> & p) const; -- optimization
 840 +      bool Project (int surfi, Point<3> & p, double& u, double& v) const;
 841        bool FastProject (int surfi, Point<3> & ap, double& u, double& v) const;
 842
 843 +      // -- Optimization: to use cached projector and classifier
 844 +      void GetFaceTools(int surfi, Handle(ShapeAnalysis_Surface)& proj,
 845 +                        BRepTopAdaptor_FClass2d*& cls) const;
 846 +
 847        OCCSurface GetSurface (int surfi)
 848        {
 849           cout << "OCCGeometry::GetSurface using PLANESPACE" << endl;
 850 diff -Naur netgen-4.9.13_orig/libsrc/occ/occmeshsurf.cpp netgen-4.9.13_new/libsrc/occ/occmeshsurf.cpp
 851 --- netgen-4.9.13_orig/libsrc/occ/occmeshsurf.cpp       2009-08-24 06:32:47.000000000 +0400
 852 +++ netgen-4.9.13_new/libsrc/occ/occmeshsurf.cpp        2013-05-31 11:28:04.000000000 +0400
 853 @@ -6,6 +6,7 @@
 854  #include <meshing.hpp>
 855  #include <GeomLProp_SLProps.hxx>
 856  #include <ShapeAnalysis_Surface.hxx>
 857 +#include <GeomAPI_ProjectPointOnCurve.hxx> // -- moved here from occgeom.hpp
 858
 859
 860  namespace netgen
 861 @@ -434,23 +435,33 @@
 862
 863    void MeshOptimize2dOCCSurfaces :: ProjectPoint (INDEX surfind, Point<3> & p) const
 864    {
 865 -    geometry.Project (surfind, p);
 866 +    // geometry.Project (surfind, p); -- signature of Project() changed for optimization
 867 +    double u, v;
 868 +    geometry.Project (surfind, p, u, v);
 869    }
 870
 871
 872    int MeshOptimize2dOCCSurfaces :: ProjectPointGI (INDEX surfind, Point<3> & p, PointGeomInfo & gi) const
 873    {
 874 -    double u = gi.u;
 875 -    double v = gi.v;
 876 +    //double u = gi.u;
 877 +    //double v = gi.v;
 878
 879      Point<3> hp = p;
 880 -    if (geometry.FastProject (surfind, hp, u, v))
 881 -      {
 882 -       p = hp;
 883 -       return 1;
 884 -      }
 885 -    ProjectPoint (surfind, p);
 886 -    return CalcPointGeomInfo (surfind, gi, p);
 887 +    // -- u and v are computed by FastProject() and Project(), no need to call CalcPointGeomInfo()
 888 +    // if (geometry.FastProject (surfind, hp, u, v))
 889 +    //   {
 890 +    //    p = hp;
 891 +    //    return 1;
 892 +    //   }
 893 +    // ProjectPoint (surfind, p);
 894 +    // return CalcPointGeomInfo (surfind, gi, p);
 895 +    bool ok;
 896 +    if (gi.trignum > 0)
 897 +      ok = geometry.FastProject (surfind, hp, gi.u, gi.v);
 898 +    else
 899 +      ok = geometry.Project (surfind, hp, gi.u, gi.v);
 900 +    p = hp;
 901 +    return ok;
 902    }
 903
 904
 905 @@ -680,7 +691,8 @@
 906         if (!geometry.FastProject (surfi, hnewp, u, v))
 907           {
 908           //  cout << "Fast projection to surface fails! Using OCC projection" << endl;
 909 -           geometry.Project (surfi, hnewp);
 910 +           // geometry.Project (surfi, hnewp); -- Project() changed for optimization
 911 +           geometry.Project (surfi, hnewp, u, v);
 912           }
 913
 914         newgi.trignum = 1;
 915 @@ -689,7 +701,7 @@
 916        }
 917
 918      newp = hnewp;
 919 -  }
 920 +  }//; -- to compile with -Wall -pedantic
 921
 922
 923    void OCCRefinementSurfaces ::
 924 @@ -708,14 +720,18 @@
 925      hnewp = Point<3> (pnt.X(), pnt.Y(), pnt.Z());
 926      newp = hnewp;
 927      newgi = ap1;
 928 -  };
 929 +  }
 930
 931
 932    void OCCRefinementSurfaces :: ProjectToSurface (Point<3> & p, int surfi)
 933    {
 934      if (surfi > 0)
 935 -      geometry.Project (surfi, p);
 936 -  };
 937 +      // geometry.Project (surfi, p); -- Project() changed for optimization
 938 +    {
 939 +      double u, v;
 940 +      geometry.Project (surfi, p, u, v);
 941 +    }
 942 +  }//; -- to compile with -Wall -pedantic
 943
 944    void OCCRefinementSurfaces :: ProjectToSurface (Point<3> & p, int surfi, PointGeomInfo & gi)
 945    {
 946 @@ -723,9 +739,10 @@
 947        if (!geometry.FastProject (surfi, p, gi.u, gi.v))
 948         {
 949           cout << "Fast projection to surface fails! Using OCC projection" << endl;
 950 -         geometry.Project (surfi, p);
 951 +          double u, v;
 952 +         geometry.Project (surfi, p, u, v);
 953         }
 954 -  };
 955 +  }
 956
 957
 958
 959 diff -Naur netgen-4.9.13_orig/libsrc/occ/utilities.h netgen-4.9.13_new/libsrc/occ/utilities.h
 960 --- netgen-4.9.13_orig/libsrc/occ/utilities.h   2009-08-24 06:12:24.000000000 +0400
 961 +++ netgen-4.9.13_new/libsrc/occ/utilities.h    2013-05-31 11:28:04.000000000 +0400
 962 @@ -33,6 +33,7 @@
 963
 964  #include <string>
 965  #include <iostream>
 966 +#include <iomanip>
 967  #include <cstdlib>
 968  // #include "SALOME_Log.hxx"
 969
 970 diff -Naur netgen-4.9.13_orig/libsrc/stlgeom/stlgeommesh.cpp netgen-4.9.13_new/libsrc/stlgeom/stlgeommesh.cpp
 971 --- netgen-4.9.13_orig/libsrc/stlgeom/stlgeommesh.cpp   2009-08-10 15:40:51.000000000 +0400
 972 +++ netgen-4.9.13_new/libsrc/stlgeom/stlgeommesh.cpp    2013-05-31 11:28:04.000000000 +0400
 973 @@ -1435,7 +1435,8 @@
 974
 975           if (!optstring || strlen(optstring) == 0)
 976             {
 977 -             mparam.optimize2d = "smcm";
 978 +             //mparam.optimize2d = (char*)"smcm";
 979 +              mparam.optimize2d = (char*)"smcm";
 980             }
 981           else
 982             {
 983 @@ -1451,7 +1452,8 @@
 984                                  mparam.grading);
 985               mesh -> LoadLocalMeshSize (mparam.meshsizefilename);
 986               mesh -> CalcLocalHFromSurfaceCurvature (stlparam.resthsurfmeshcurvfac);
 987 -             mparam.optimize2d = "cmsmSm";
 988 +             //mparam.optimize2d = (char*)"cmsmSm";
 989 +              mparam.optimize2d = (char*)"cmsmSm";
 990               STLSurfaceOptimization (*stlgeometry, *mesh, mparam);
 991  #ifdef STAT_STREAM
 992               (*statout) << GetTime() << " & ";
 993 @@ -1557,7 +1559,8 @@
 994
 995           if (!optstring || strlen(optstring) == 0)
 996             {
 997 -             mparam.optimize3d = "cmdmstm";
 998 +              //mparam.optimize3d = "cmdmstm";
 999 +             mparam.optimize3d = (char*)"cmdmstm";
1000             }
1001           else
1002             {
1003 diff -Naur netgen-4.9.13_orig/nglib/nglib.cpp netgen-4.9.13_new/nglib/nglib.cpp
1004 --- netgen-4.9.13_orig/nglib/nglib.cpp  2010-05-18 15:19:15.000000000 +0400
1005 +++ netgen-4.9.13_new/nglib/nglib.cpp   2013-05-31 11:33:44.000000000 +0400
1006 @@ -66,9 +66,14 @@
1007  namespace nglib
1008  {
1009     // initialize, deconstruct Netgen library:
1010 -   DLL_HEADER void Ng_Init ()
1011 +   DLL_HEADER void Ng_Init ( const char* filename )
1012     {
1013 -      mycout = &cout;
1014 +      if ( filename == "" ) {
1015 +        mycout = &cout;
1016 +      } else {
1017 +        mycout = new ofstream ( filename );
1018 +        cout << "NOTE: netgen output was redirected to file " << filename << endl;
1019 +      }
1020        myerr = &cerr;
1021        // netgen::testout->SetOutStream (new ofstream ("test.out"));
1022        testout = new ofstream ("test.out");
1023 diff -Naur netgen-4.9.13_orig/nglib/nglib.h netgen-4.9.13_new/nglib/nglib.h
1024 --- netgen-4.9.13_orig/nglib/nglib.h    2010-05-18 15:20:25.000000000 +0400
1025 +++ netgen-4.9.13_new/nglib/nglib.h     2013-05-31 11:29:26.000000000 +0400
1026 @@ -24,7 +24,7 @@
1027  // Philippose - 14.02.2009
1028  // Modifications for creating a DLL in Windows
1029  #ifdef WIN32
1030 -   #ifdef NGLIB_EXPORTS || nglib_EXPORTS
1031 +   #if defined NGLIB_EXPORTS || defined nglib_EXPORTS
1032        #define DLL_HEADER   __declspec(dllexport)
1033     #else
1034        #define DLL_HEADER   __declspec(dllimport)
1035 @@ -166,7 +166,7 @@
1036      program before beginning to use the other Netgen
1037      specific functions.
1038  */
1039 -DLL_HEADER void Ng_Init ();
1040 +DLL_HEADER void Ng_Init ( const char* filename = "" );
1041
1042
1043  /*! \brief Exit the Netgen meshing kernel in a clean manner