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       2014-01-15 13:59:16.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       2014-01-15 13:59:16.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     2014-01-15 13:59:16.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  2014-01-15 13:59:16.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     2014-01-15 13:59:16.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   2014-01-15 13:59:16.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  2014-01-15 13:59:16.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    2014-01-15 13:59:16.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 2014-01-15 13:59:16.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 @@ -497,17 +564,20 @@
 527                    bool exists = 0;
 528                    int j;
 529                    for (j = first_ep; j <= mesh.GetNP(); j++)
 530 +                  {
 531 +                     if (!merge_solids && mesh.Point(j).GetLayer() != geomedgenr ) continue; // to support SALOME fuse edges
 532                       if ((mesh.Point(j)-Point<3>(mp[i-1])).Length() < eps)
 533                       {
 534                          exists = 1;
 535                          break;
 536                       }
 537 +                  }
 538
 539                       if (exists)
 540                          pnums[i] = j;
 541                       else
 542                       {
 543 -                        mesh.AddPoint (mp[i-1]);
 544 +                        mesh.AddPoint (mp[i-1], geomedgenr); // to support SALOME fuse edges
 545                          (*testout) << "add meshpoint " << mp[i-1] << endl;
 546                          pnums[i] = mesh.GetNP();
 547                       }
 548 @@ -591,6 +661,8 @@
 549        //               (*testout) << "edge " << mesh.LineSegment(i).edgenr << " face " << mesh.LineSegment(i).si
 550        //                               << " p1 " << mesh.LineSegment(i)[0] << " p2 " << mesh.LineSegment(i)[1] << endl;
 551        //       exit(10);
 552 +      for (int j = 1; j <= mesh.GetNP(); j++) // to support SALOME fuse edges: set level to zero\r
 553 +        mesh.Point(j) = MeshPoint( (Point<3>&) mesh.Point(j) );
 554
 555        mesh.CalcSurfacesOfNode();
 556        multithread.task = savetask;
 557 @@ -633,7 +705,8 @@
 558           }
 559
 560           (*testout) << "mesh face " << k << endl;
 561 -         multithread.percent = 100 * k / (mesh.GetNFD()+1e-10);
 562 +         // multithread.percent = 100 * k / (mesh.GetNFD()+1e-10); -- avoid FPE
 563 +         multithread.percent = 100 * k / (mesh.GetNFD() + VSMALL);
 564           geom.facemeshstatus[k-1] = -1;
 565
 566
 567 @@ -901,7 +974,8 @@
 568           //      if (k != 36) continue;
 569
 570           //      (*testout) << "optimize face " << k << endl;
 571 -         multithread.percent = 100 * k / (mesh.GetNFD()+1e-10);
 572 +         //multithread.percent = 100 * k / (mesh.GetNFD()+1e-10); -- avoid FPE
 573 +         multithread.percent = 100 * k / (mesh.GetNFD() + VSMALL);
 574
 575           FaceDescriptor & fd = mesh.GetFaceDescriptor(k);
 576
 577 @@ -1456,3 +1530,4 @@
 578  }
 579
 580  #endif
 581 +
 582 diff -Naur netgen-4.9.13_orig/libsrc/occ/occgeom.cpp netgen-4.9.13_new/libsrc/occ/occgeom.cpp
 583 --- netgen-4.9.13_orig/libsrc/occ/occgeom.cpp   2010-03-05 16:16:21.000000000 +0300
 584 +++ netgen-4.9.13_new/libsrc/occ/occgeom.cpp    2014-01-15 13:59:16.000000000 +0400
 585 @@ -8,6 +8,8 @@
 586  #include "ShapeAnalysis_CheckSmallFace.hxx"
 587  #include "ShapeAnalysis_DataMapOfShapeListOfReal.hxx"
 588  #include "ShapeAnalysis_Surface.hxx"
 589 +#include <BRepTopAdaptor_FClass2d.hxx> // -- to optimize Project() and FastProject()
 590 +#include <TopAbs_State.hxx>
 591  #include "BRepAlgoAPI_Fuse.hxx"
 592  #include "BRepCheck_Analyzer.hxx"
 593  #include "BRepLib.hxx"
 594 @@ -16,10 +18,17 @@
 595  #include "ShapeFix_FixSmallFace.hxx"
 596  #include "Partition_Spliter.hxx"
 597
 598 -
 599  namespace netgen
 600  {
 601 -   void OCCGeometry :: PrintNrShapes ()
 602 +  // free data used to optimize Project() and FastProject()
 603 +  OCCGeometry::~OCCGeometry()
 604 +  {
 605 +    NCollection_DataMap<int,BRepTopAdaptor_FClass2d*>::Iterator it(fclsmap);
 606 +    for (; it.More(); it.Next())
 607 +      delete it.Value();
 608 +  }
 609 +
 610 +  void OCCGeometry :: PrintNrShapes ()
 611     {
 612        TopExp_Explorer e;
 613        int count = 0;
 614 @@ -951,25 +960,58 @@
 615     }
 616
 617
 618 +   // returns a projector and a classifier for the given surface
 619 +   void OCCGeometry::GetFaceTools(int surfi, Handle(ShapeAnalysis_Surface)& proj,
 620 +                                  BRepTopAdaptor_FClass2d*& cls) const
 621 +   {
 622 +     //MSV: organize caching projector in the map
 623 +     if (fprjmap.IsBound(surfi))
 624 +     {
 625 +       proj = fprjmap.Find(surfi);
 626 +       cls = fclsmap.Find(surfi);
 627 +     }
 628 +     else
 629 +     {
 630 +       const TopoDS_Face& aFace = TopoDS::Face(fmap(surfi));
 631 +       Handle(Geom_Surface) aSurf = BRep_Tool::Surface(aFace);
 632 +       proj = new ShapeAnalysis_Surface(aSurf);
 633 +       fprjmap.Bind(surfi, proj);
 634 +       cls = new BRepTopAdaptor_FClass2d(aFace,Precision::Confusion());
 635 +       fclsmap.Bind(surfi, cls);
 636 +     }
 637 +   }
 638
 639 -
 640 -   void OCCGeometry :: Project (int surfi, Point<3> & p) const
 641 +   // void OCCGeometry :: Project (int surfi, Point<3> & p) const
 642 +   bool OCCGeometry :: Project (int surfi, Point<3> & p, double& u, double& v) const
 643     {
 644        static int cnt = 0;
 645        if (++cnt % 1000 == 0) cout << "Project cnt = " << cnt << endl;
 646
 647        gp_Pnt pnt(p(0), p(1), p(2));
 648
 649 -      double u,v;
 650 -      Handle( Geom_Surface ) thesurf = BRep_Tool::Surface(TopoDS::Face(fmap(surfi)));
 651 -      Handle( ShapeAnalysis_Surface ) su = new ShapeAnalysis_Surface( thesurf );
 652 -      gp_Pnt2d suval = su->ValueOfUV ( pnt, BRep_Tool::Tolerance( TopoDS::Face(fmap(surfi)) ) );
 653 -      suval.Coord( u, v);
 654 -      pnt = thesurf->Value( u, v );
 655 -
 656 -
 657 +      // -- Optimization: use cached projector and classifier
 658 +      // double u,v;
 659 +      // Handle( Geom_Surface ) thesurf = BRep_Tool::Surface(TopoDS::Face(fmap(surfi)));
 660 +      // Handle( ShapeAnalysis_Surface ) su = new ShapeAnalysis_Surface( thesurf );
 661 +      // gp_Pnt2d suval = su->ValueOfUV ( pnt, BRep_Tool::Tolerance( TopoDS::Face(fmap(surfi)) ) );
 662 +      // suval.Coord( u, v);
 663 +      // pnt = thesurf->Value( u, v );
 664 +
 665 +      Handle(ShapeAnalysis_Surface) proj;
 666 +      BRepTopAdaptor_FClass2d *cls;
 667 +      GetFaceTools(surfi, proj, cls);
 668 +
 669 +      gp_Pnt2d p2d = proj->ValueOfUV(pnt, Precision::Confusion());
 670 +      if (cls->Perform(p2d) == TopAbs_OUT)
 671 +      {
 672 +        return false;
 673 +      }
 674 +      pnt = proj->Value(p2d);
 675 +      p2d.Coord(u, v);
 676 +
 677        p = Point<3> (pnt.X(), pnt.Y(), pnt.Z());
 678
 679 +      return true;
 680     }
 681
 682
 683 @@ -979,54 +1021,69 @@
 684     {
 685        gp_Pnt p(ap(0), ap(1), ap(2));
 686
 687 -      Handle(Geom_Surface) surface = BRep_Tool::Surface(TopoDS::Face(fmap(surfi)));
 688 -
 689 -      gp_Pnt x = surface->Value (u,v);
 690 -
 691 -      if (p.SquareDistance(x) <= sqr(PROJECTION_TOLERANCE)) return true;
 692 -
 693 -      gp_Vec du, dv;
 694 -
 695 -      surface->D1(u,v,x,du,dv);
 696 -
 697 -      int count = 0;
 698 -
 699 -      gp_Pnt xold;
 700 -      gp_Vec n;
 701 -      double det, lambda, mu;
 702 -
 703 -      do {
 704 -         count++;
 705 -
 706 -         n = du^dv;
 707 -
 708 -         det = Det3 (n.X(), du.X(), dv.X(),
 709 -            n.Y(), du.Y(), dv.Y(),
 710 -            n.Z(), du.Z(), dv.Z());
 711 -
 712 -         if (det < 1e-15) return false;
 713 -
 714 -         lambda = Det3 (n.X(), p.X()-x.X(), dv.X(),
 715 -            n.Y(), p.Y()-x.Y(), dv.Y(),
 716 -            n.Z(), p.Z()-x.Z(), dv.Z())/det;
 717 -
 718 -         mu     = Det3 (n.X(), du.X(), p.X()-x.X(),
 719 -            n.Y(), du.Y(), p.Y()-x.Y(),
 720 -            n.Z(), du.Z(), p.Z()-x.Z())/det;
 721 -
 722 -         u += lambda;
 723 -         v += mu;
 724 -
 725 -         xold = x;
 726 -         surface->D1(u,v,x,du,dv);
 727 -
 728 -      } while (xold.SquareDistance(x) > sqr(PROJECTION_TOLERANCE) && count < 50);
 729 -
 730 -      //    (*testout) << "FastProject count: " << count << endl;
 731 -
 732 -      if (count == 50) return false;
 733 -
 734 -      ap = Point<3> (x.X(), x.Y(), x.Z());
 735 +      // -- Optimization: use cached projector and classifier
 736 +      // Handle(Geom_Surface) surface = BRep_Tool::Surface(TopoDS::Face(fmap(surfi)));
 737 +      //
 738 +      // gp_Pnt x = surface->Value (u,v);
 739 +      //
 740 +      // if (p.SquareDistance(x) <= sqr(PROJECTION_TOLERANCE)) return true;
 741 +      //
 742 +      // gp_Vec du, dv;
 743 +      //
 744 +      // surface->D1(u,v,x,du,dv);
 745 +      //
 746 +      // int count = 0;
 747 +      //
 748 +      // gp_Pnt xold;
 749 +      // gp_Vec n;
 750 +      // double det, lambda, mu;
 751 +      //
 752 +      // do {
 753 +      //    count++;
 754 +      //
 755 +      //    n = du^dv;
 756 +      //
 757 +      //    det = Det3 (n.X(), du.X(), dv.X(),
 758 +      //       n.Y(), du.Y(), dv.Y(),
 759 +      //       n.Z(), du.Z(), dv.Z());
 760 +      //
 761 +      //    if (det < 1e-15) return false;
 762 +      //
 763 +      //    lambda = Det3 (n.X(), p.X()-x.X(), dv.X(),
 764 +      //       n.Y(), p.Y()-x.Y(), dv.Y(),
 765 +      //       n.Z(), p.Z()-x.Z(), dv.Z())/det;
 766 +      //
 767 +      //    mu     = Det3 (n.X(), du.X(), p.X()-x.X(),
 768 +      //       n.Y(), du.Y(), p.Y()-x.Y(),
 769 +      //       n.Z(), du.Z(), p.Z()-x.Z())/det;
 770 +      //
 771 +      //    u += lambda;
 772 +      //    v += mu;
 773 +      //
 774 +      //    xold = x;
 775 +      //    surface->D1(u,v,x,du,dv);
 776 +      //
 777 +      // } while (xold.SquareDistance(x) > sqr(PROJECTION_TOLERANCE) && count < 50);
 778 +      //
 779 +      // //    (*testout) << "FastProject count: " << count << endl;
 780 +      //
 781 +      // if (count == 50) return false;
 782 +      //
 783 +      // ap = Point<3> (x.X(), x.Y(), x.Z());
 784 +      Handle(ShapeAnalysis_Surface) proj;
 785 +      BRepTopAdaptor_FClass2d *cls;
 786 +      GetFaceTools(surfi, proj, cls);
 787 +
 788 +      gp_Pnt2d p2d = proj->NextValueOfUV(gp_Pnt2d(u,v), p, Precision::Confusion());
 789 +      if (cls->Perform(p2d) == TopAbs_OUT)
 790 +      {
 791 +        //cout << "Projection fails" << endl;
 792 +        return false;
 793 +      }
 794 +
 795 +      p = proj->Value(p2d);
 796 +      p2d.Coord(u, v);
 797 +      ap = Point<3> (p.X(), p.Y(), p.Z());
 798
 799        return true;
 800     }
 801 diff -Naur netgen-4.9.13_orig/libsrc/occ/occgeom.hpp netgen-4.9.13_new/libsrc/occ/occgeom.hpp
 802 --- netgen-4.9.13_orig/libsrc/occ/occgeom.hpp   2010-01-14 19:56:19.000000000 +0300
 803 +++ netgen-4.9.13_new/libsrc/occ/occgeom.hpp    2014-01-15 13:59:16.000000000 +0400
 804 @@ -15,8 +15,8 @@
 805  #include "Geom_Curve.hxx"
 806  #include "Geom2d_Curve.hxx"
 807  #include "Geom_Surface.hxx"
 808 -#include "GeomAPI_ProjectPointOnSurf.hxx"
 809 -#include "GeomAPI_ProjectPointOnCurve.hxx"
 810 +// #include "GeomAPI_ProjectPointOnSurf.hxx"
 811 +// #include "GeomAPI_ProjectPointOnCurve.hxx"
 812  #include "BRepTools.hxx"
 813  #include "TopExp.hxx"
 814  #include "BRepBuilderAPI_MakeVertex.hxx"
 815 @@ -42,8 +42,8 @@
 816  #include "Geom_Curve.hxx"
 817  #include "Geom2d_Curve.hxx"
 818  #include "Geom_Surface.hxx"
 819 -#include "GeomAPI_ProjectPointOnSurf.hxx"
 820 -#include "GeomAPI_ProjectPointOnCurve.hxx"
 821 +// #include "GeomAPI_ProjectPointOnSurf.hxx"
 822 +// #include "GeomAPI_ProjectPointOnCurve.hxx"
 823  #include "TopoDS_Wire.hxx"
 824  #include "BRepTools_WireExplorer.hxx"
 825  #include "BRepTools.hxx"
 826 @@ -68,7 +68,7 @@
 827  #include "IGESToBRep_Reader.hxx"
 828  #include "Interface_Static.hxx"
 829  #include "GeomAPI_ExtremaCurveCurve.hxx"
 830 -#include "Standard_ErrorHandler.hxx"
 831 +//#include "Standard_ErrorHandler.hxx"
 832  #include "Standard_Failure.hxx"
 833  #include "ShapeUpgrade_ShellSewing.hxx"
 834  #include "ShapeFix_Shape.hxx"
 835 @@ -80,6 +80,10 @@
 836  #include "ShapeAnalysis.hxx"
 837  #include "ShapeBuild_ReShape.hxx"
 838
 839 +// -- Optimization: to use cached projector and classifier
 840 +#include <NCollection_DataMap.hxx>
 841 +class Handle_ShapeAnalysis_Surface;
 842 +class BRepTopAdaptor_FClass2d;
 843
 844  // Philippose - 29/01/2009
 845  // OpenCascade XDE Support
 846 @@ -190,6 +194,9 @@
 847     class OCCGeometry : public NetgenGeometry
 848     {
 849        Point<3> center;
 850 +      // -- Optimization: to use cached projector and classifier
 851 +      mutable NCollection_DataMap<int,Handle_ShapeAnalysis_Surface> fprjmap;
 852 +      mutable NCollection_DataMap<int,BRepTopAdaptor_FClass2d*> fclsmap;
 853
 854     public:
 855        TopoDS_Shape shape;
 856 @@ -241,6 +248,8 @@
 857           vmap.Clear();
 858        }
 859
 860 +      ~OCCGeometry();      // -- to free cached projector and classifier
 861 +
 862        void BuildFMap();
 863
 864        Box<3> GetBoundingBox()
 865 @@ -260,9 +269,14 @@
 866        Point<3> Center()
 867        {  return center;}
 868
 869 -      void Project (int surfi, Point<3> & p) const;
 870 +      // void Project (int surfi, Point<3> & p) const; -- optimization
 871 +      bool Project (int surfi, Point<3> & p, double& u, double& v) const;
 872        bool FastProject (int surfi, Point<3> & ap, double& u, double& v) const;
 873
 874 +      // -- Optimization: to use cached projector and classifier
 875 +      void GetFaceTools(int surfi, Handle(ShapeAnalysis_Surface)& proj,
 876 +                        BRepTopAdaptor_FClass2d*& cls) const;
 877 +
 878        OCCSurface GetSurface (int surfi)
 879        {
 880           cout << "OCCGeometry::GetSurface using PLANESPACE" << endl;
 881 diff -Naur netgen-4.9.13_orig/libsrc/occ/occmeshsurf.cpp netgen-4.9.13_new/libsrc/occ/occmeshsurf.cpp
 882 --- netgen-4.9.13_orig/libsrc/occ/occmeshsurf.cpp       2009-08-24 06:32:47.000000000 +0400
 883 +++ netgen-4.9.13_new/libsrc/occ/occmeshsurf.cpp        2014-01-15 13:59:16.000000000 +0400
 884 @@ -6,6 +6,7 @@
 885  #include <meshing.hpp>
 886  #include <GeomLProp_SLProps.hxx>
 887  #include <ShapeAnalysis_Surface.hxx>
 888 +#include <GeomAPI_ProjectPointOnCurve.hxx> // -- moved here from occgeom.hpp
 889
 890
 891  namespace netgen
 892 @@ -434,23 +435,33 @@
 893
 894    void MeshOptimize2dOCCSurfaces :: ProjectPoint (INDEX surfind, Point<3> & p) const
 895    {
 896 -    geometry.Project (surfind, p);
 897 +    // geometry.Project (surfind, p); -- signature of Project() changed for optimization
 898 +    double u, v;
 899 +    geometry.Project (surfind, p, u, v);
 900    }
 901
 902
 903    int MeshOptimize2dOCCSurfaces :: ProjectPointGI (INDEX surfind, Point<3> & p, PointGeomInfo & gi) const
 904    {
 905 -    double u = gi.u;
 906 -    double v = gi.v;
 907 +    //double u = gi.u;
 908 +    //double v = gi.v;
 909
 910      Point<3> hp = p;
 911 -    if (geometry.FastProject (surfind, hp, u, v))
 912 -      {
 913 -       p = hp;
 914 -       return 1;
 915 -      }
 916 -    ProjectPoint (surfind, p);
 917 -    return CalcPointGeomInfo (surfind, gi, p);
 918 +    // -- u and v are computed by FastProject() and Project(), no need to call CalcPointGeomInfo()
 919 +    // if (geometry.FastProject (surfind, hp, u, v))
 920 +    //   {
 921 +    //    p = hp;
 922 +    //    return 1;
 923 +    //   }
 924 +    // ProjectPoint (surfind, p);
 925 +    // return CalcPointGeomInfo (surfind, gi, p);
 926 +    bool ok;
 927 +    if (gi.trignum > 0)
 928 +      ok = geometry.FastProject (surfind, hp, gi.u, gi.v);
 929 +    else
 930 +      ok = geometry.Project (surfind, hp, gi.u, gi.v);
 931 +    p = hp;
 932 +    return ok;
 933    }
 934
 935
 936 @@ -680,7 +691,8 @@
 937         if (!geometry.FastProject (surfi, hnewp, u, v))
 938           {
 939           //  cout << "Fast projection to surface fails! Using OCC projection" << endl;
 940 -           geometry.Project (surfi, hnewp);
 941 +           // geometry.Project (surfi, hnewp); -- Project() changed for optimization
 942 +           geometry.Project (surfi, hnewp, u, v);
 943           }
 944
 945         newgi.trignum = 1;
 946 @@ -689,7 +701,7 @@
 947        }
 948
 949      newp = hnewp;
 950 -  }
 951 +  }//; -- to compile with -Wall -pedantic
 952
 953
 954    void OCCRefinementSurfaces ::
 955 @@ -708,14 +720,18 @@
 956      hnewp = Point<3> (pnt.X(), pnt.Y(), pnt.Z());
 957      newp = hnewp;
 958      newgi = ap1;
 959 -  };
 960 +  }
 961
 962
 963    void OCCRefinementSurfaces :: ProjectToSurface (Point<3> & p, int surfi)
 964    {
 965      if (surfi > 0)
 966 -      geometry.Project (surfi, p);
 967 -  };
 968 +      // geometry.Project (surfi, p); -- Project() changed for optimization
 969 +    {
 970 +      double u, v;
 971 +      geometry.Project (surfi, p, u, v);
 972 +    }
 973 +  }//; -- to compile with -Wall -pedantic
 974
 975    void OCCRefinementSurfaces :: ProjectToSurface (Point<3> & p, int surfi, PointGeomInfo & gi)
 976    {
 977 @@ -723,9 +739,10 @@
 978        if (!geometry.FastProject (surfi, p, gi.u, gi.v))
 979         {
 980           cout << "Fast projection to surface fails! Using OCC projection" << endl;
 981 -         geometry.Project (surfi, p);
 982 +          double u, v;
 983 +         geometry.Project (surfi, p, u, v);
 984         }
 985 -  };
 986 +  }
 987
 988
 989
 990 diff -Naur netgen-4.9.13_orig/libsrc/occ/utilities.h netgen-4.9.13_new/libsrc/occ/utilities.h
 991 --- netgen-4.9.13_orig/libsrc/occ/utilities.h   2009-08-24 06:12:24.000000000 +0400
 992 +++ netgen-4.9.13_new/libsrc/occ/utilities.h    2014-01-15 13:59:16.000000000 +0400
 993 @@ -33,6 +33,7 @@
 994
 995  #include <string>
 996  #include <iostream>
 997 +#include <iomanip>
 998  #include <cstdlib>
 999  // #include "SALOME_Log.hxx"
1000
1001 diff -Naur netgen-4.9.13_orig/libsrc/stlgeom/stlgeommesh.cpp netgen-4.9.13_new/libsrc/stlgeom/stlgeommesh.cpp
1002 --- netgen-4.9.13_orig/libsrc/stlgeom/stlgeommesh.cpp   2009-08-10 15:40:51.000000000 +0400
1003 +++ netgen-4.9.13_new/libsrc/stlgeom/stlgeommesh.cpp    2014-01-15 13:59:16.000000000 +0400
1004 @@ -1435,7 +1435,8 @@
1005
1006           if (!optstring || strlen(optstring) == 0)
1007             {
1008 -             mparam.optimize2d = "smcm";
1009 +             //mparam.optimize2d = (char*)"smcm";
1010 +              mparam.optimize2d = (char*)"smcm";
1011             }
1012           else
1013             {
1014 @@ -1451,7 +1452,8 @@
1015                                  mparam.grading);
1016               mesh -> LoadLocalMeshSize (mparam.meshsizefilename);
1017               mesh -> CalcLocalHFromSurfaceCurvature (stlparam.resthsurfmeshcurvfac);
1018 -             mparam.optimize2d = "cmsmSm";
1019 +             //mparam.optimize2d = (char*)"cmsmSm";
1020 +              mparam.optimize2d = (char*)"cmsmSm";
1021               STLSurfaceOptimization (*stlgeometry, *mesh, mparam);
1022  #ifdef STAT_STREAM
1023               (*statout) << GetTime() << " & ";
1024 @@ -1557,7 +1559,8 @@
1025
1026           if (!optstring || strlen(optstring) == 0)
1027             {
1028 -             mparam.optimize3d = "cmdmstm";
1029 +              //mparam.optimize3d = "cmdmstm";
1030 +             mparam.optimize3d = (char*)"cmdmstm";
1031             }
1032           else
1033             {
1034 diff -Naur netgen-4.9.13_orig/nglib/nglib.h netgen-4.9.13_new/nglib/nglib.h
1035 --- netgen-4.9.13_orig/nglib/nglib.h    2010-05-18 15:20:25.000000000 +0400
1036 +++ netgen-4.9.13_new/nglib/nglib.h     2014-01-15 13:59:16.000000000 +0400
1037 @@ -24,7 +24,7 @@
1038  // Philippose - 14.02.2009
1039  // Modifications for creating a DLL in Windows
1040  #ifdef WIN32
1041 -   #ifdef NGLIB_EXPORTS || nglib_EXPORTS
1042 +   #if defined NGLIB_EXPORTS || defined nglib_EXPORTS
1043        #define DLL_HEADER   __declspec(dllexport)
1044     #else
1045        #define DLL_HEADER   __declspec(dllimport)