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