products/patches/netgen62ForSalome_windows.patch

   1 diff --git a/libsrc/occ/occgenmesh.cpp b/libsrc/occ/occgenmesh.cpp
   2 index 314d405a..7c912bcc 100644
   3 --- a/libsrc/occ/occgenmesh.cpp
   4 +++ b/libsrc/occ/occgenmesh.cpp
   5 @@ -49,7 +49,7 @@ namespace netgen
   6    double ComputeH (double kappa, const MeshingParameters & mparam)
   7    {
   8      kappa *= mparam.curvaturesafety;
   9 -    /*
  10 +    /**/
  11      double hret;
  12
  13      if (mparam.maxh * kappa < 1)
  14 @@ -61,7 +61,7 @@ namespace netgen
  15        hret = mparam.maxh;
  16
  17      return hret;
  18 -    */
  19 +    /**/
  20      // return min(mparam.maxh, 1/kappa);
  21      return (mparam.maxh*kappa < 1) ? mparam.maxh : 1/kappa;
  22    }
  23 @@ -253,10 +253,12 @@ namespace netgen
  24        {
  25          oldpnt = pnt;
  26          pnt = c->Value(s0+(i/double(DIVIDEEDGESECTIONS))*(s1-s0));
  27 +        // !!! no more than 1 segment per <edge length>/DIVIDEEDGESECTIONS
  28          hvalue[i] = hvalue[i-1] +
  29 -          1.0/mesh.GetH(Point3d(pnt.X(), pnt.Y(), pnt.Z()))*
  30 -          pnt.Distance(oldpnt);
  31 +          min( 1.0,
  32 +               1.0/mesh.GetH(Point3d(pnt.X(), pnt.Y(), pnt.Z()))*
  33 +               pnt.Distance(oldpnt));
  34
  35          //(*testout) << "mesh.GetH(Point3d(pnt.X(), pnt.Y(), pnt.Z())) " << mesh.GetH(Point3d(pnt.X(), pnt.Y(), pnt.Z()))
  36          //        <<  " pnt.Distance(oldpnt) " << pnt.Distance(oldpnt) << endl;
  37
  38 @@ -299,7 +301,7 @@ namespace netgen
  39          cout << "CORRECTED" << endl;
  40          ps.SetSize (nsubedges-2);
  41          params.SetSize (nsubedges);
  42 -        params[nsubedges] = s1;
  43 +        params[nsubedges-1] = s1;
  44        }
  45    }
  46
  47 @@ -323,6 +325,8 @@ namespace netgen
  48
  49      double eps = 1e-6 * geom.GetBoundingBox().Diam();
  50
  51 +    int first_vp = mesh.GetNP()+1;
  52 +
  53      tsearch.Start();
  54      for (int i = 1; i <= nvertices; i++)
  55        {
  56 @@ -481,22 +485,107 @@ namespace netgen
  57
  58                  if (!merge_solids)
  59                    {
  60 -                    pnums[0] = geom.vmap.FindIndex (TopExp::FirstVertex (edge)) + PointIndex::BASE-1;
  61 -                    pnums.Last() = geom.vmap.FindIndex (TopExp::LastVertex (edge)) + PointIndex::BASE-1;
  62 +                    //pnums[0] = geom.vmap.FindIndex (TopExp::FirstVertex (edge)) + PointIndex::BASE-1;
  63 +                    //pnums.Last() = geom.vmap.FindIndex (TopExp::LastVertex (edge)) + PointIndex::BASE-1;
  64 +                    const int dpi = PointIndex::BASE-1;
  65 +                    MeshPoint dfltP ( Point<3> ( 0, 0, 0 ) );
  66 +                    PointIndex *ipp[] = { &pnums[0], &pnums.Last() };
  67 +                    TopoDS_Iterator vIt( edge, false );
  68 +                    TopoDS_Vertex v[2];
  69 +                    v[0] = TopoDS::Vertex( vIt.Value() ); vIt.Next();
  70 +                    v[1] = TopoDS::Vertex( vIt.Value() );
  71 +                    if ( v[0].Orientation() == TopAbs_REVERSED )
  72 +                      std::swap( v[0], v[1] );
  73 +                    for ( int i = 0; i < 2; ++i)
  74 +                    {
  75 +                      PointIndex &ip = *ipp[i];
  76 +                      ip = geom.vmap.FindIndex ( v[i] ) + dpi;
  77 +                      if ( ip == dpi || ip > nvertices + dpi )
  78 +                      {
  79 +                        PointIndex iv = ip;
  80 +                        if ( ip == dpi )
  81 +                          ip = iv = const_cast<OCCGeometry&>(geom).vmap.Add( v[i] );
  82 +                        gp_Pnt pnt = BRep_Tool::Pnt( v[i] );
  83 +                        MeshPoint mp( Point<3>(pnt.X(), pnt.Y(), pnt.Z()) );
  84 +                        for (PointIndex pi = 1; pi < first_vp; pi++)
  85 +                          if ( Dist2 (mesh.Point(pi), Point<3>(mp)) < 1e-100 )
  86 +                          {
  87 +                            ip = pi;
  88 +                            if ( mesh.Point(ip).GetLayer() != dfltP.GetLayer() && mesh.Point(ip).GetLayer() != iv )
  89 +                              continue;
  90 +                            if ( mesh.Point(ip).GetLayer() == dfltP.GetLayer())
  91 +                              mesh.Point(ip) = MeshPoint( mesh.Point(ip), iv );
  92 +                            break;
  93 +                          }
  94 +                      }
  95 +                      else
  96 +                      {
  97 +                        ip += first_vp - 1;
  98 +                      }
  99 +                    }
 100                    }
 101                  else
 102                    {
 103 -                    Point<3> fp = occ2ng (BRep_Tool::Pnt (TopExp::FirstVertex (edge)));
 104 -                    Point<3> lp = occ2ng (BRep_Tool::Pnt (TopExp::LastVertex (edge)));
 105 +                    TopoDS_Iterator vIt( edge, false );
 106 +                    TopoDS_Vertex v1 = TopoDS::Vertex( vIt.Value() ); vIt.Next();
 107 +                    TopoDS_Vertex v2 = TopoDS::Vertex( vIt.Value() );
 108 +                    if ( v1.Orientation() == TopAbs_REVERSED )
 109 +                      std::swap( v1, v2 );
 110 +                    const bool isClosedEdge = v1.IsSame( v2 );
 111 +
 112 +                    Point<3> fp = occ2ng (BRep_Tool::Pnt (v1));
 113 +                    Point<3> lp = occ2ng (BRep_Tool::Pnt (v2));
 114 +                    double tol2 = std::min( eps*eps, 1e-6 * Dist2( fp, lp ));
 115 +                    if ( isClosedEdge )
 116 +                      tol2 = BRep_Tool::Tolerance( v1 ) * BRep_Tool::Tolerance( v1 );
 117
 118                      pnums[0] = PointIndex::INVALID;
 119                      pnums.Last() = PointIndex::INVALID;
 120                      for (PointIndex pi : vertexrange)
 121                        {
 122 -                        if (Dist2 (mesh[pi], fp) < eps*eps) pnums[0] = pi;
 123 -                        if (Dist2 (mesh[pi], lp) < eps*eps) pnums.Last() = pi;
 124 +                        if (Dist2 (mesh[pi], fp) < tol2) pnums[0] = pi;
 125 +                        if (Dist2 (mesh[pi], lp) < tol2) pnums.Last() = pi;
 126                        }
 127 -                  }
 128 +                    if (( isClosedEdge && pnums[0] != pnums.Last() ) ||
 129 +                        ( !isClosedEdge && pnums[0] == pnums.Last() ))
 130 +                      pnums[0] = pnums.Last() = PointIndex::INVALID;
 131 +                    if ( pnums[0] < PointIndex::BASE || pnums.Last() < PointIndex::BASE )
 132 +                    {
 133 +                      // take into account a possible large gap between a vertex and an edge curve
 134 +                      // end and a large vertex tolerance covering the whole edge
 135 +                      if ( pnums[0] < PointIndex::BASE )
 136 +                      {
 137 +                        double tol = BRep_Tool::Tolerance( v1 );
 138 +                        for (PointIndex pi = 1; pi < first_ep; pi++)
 139 +                          if (pi != pnums.Last() && Dist2 (mesh[pi], fp) < 2*tol*tol)
 140 +                            pnums[0] = pi;
 141 +
 142 +                        if ( pnums[0] < PointIndex::BASE )
 143 +                          pnums[0] = first_ep-1 - nvertices + geom.vmap.FindIndex ( v1 );
 144 +                      }
 145 +                      if ( isClosedEdge )
 146 +                      {
 147 +                        pnums.Last() = pnums[0];
 148 +                      }
 149 +                      else
 150 +                      {
 151 +                        if ( pnums.Last() < PointIndex::BASE )
 152 +                        {
 153 +                          double tol = BRep_Tool::Tolerance( v2 );
 154 +                          for (PointIndex pi = 1; pi < first_ep; pi++)
 155 +                            if (pi != pnums[0] && Dist2 (mesh[pi], lp) < 2*tol*tol)
 156 +                              pnums.Last() = pi;
 157 +
 158 +                          if ( pnums.Last() < PointIndex::BASE )
 159 +                            pnums.Last() = first_ep-1-nvertices + geom.vmap.FindIndex ( v2 );
 160 +                        }
 161 +
 162 +                        if ( Dist2( fp, mesh[PointIndex(pnums[0])]) >
 163 +                             Dist2( lp, mesh[PointIndex(pnums.Last())]))
 164 +                          std::swap( pnums[0], pnums.Last() );
 165 +                      }
 166 +                    }
 167 +                 }
 168
 169                  for (size_t i = 1; i <= mp.Size(); i++)
 170                    {
 171 @@ -505,17 +594,19 @@ namespace netgen
 172
 173                      // for (PointIndex j = first_ep; j < mesh.Points().End(); j++)
 174                      for (PointIndex j = first_ep; j < *mesh.Points().Range().end(); j++)
 175 +                    {
 176 +                      if (!merge_solids && mesh.Point(j).GetLayer() != geomedgenr ) continue;
 177                        if ((mesh.Point(j)-Point<3>(mp[i-1])).Length() < eps)
 178                          {
 179                            exists = true;
 180                            pnums[i] = j;
 181                            break;
 182                          }
 183 -
 184 +                    }
 185                      tsearch.Stop();
 186
 187                      if (!exists)
 188 -                      pnums[i] = mesh.AddPoint (mp[i-1]);
 189 +                      pnums[i] = mesh.AddPoint (mp[i-1], geomedgenr);
 190                    }
 191                  if(geom.enames.Size() && geom.enames[curr-1] != "")
 192                    mesh.SetCD2Name(geomedgenr, geom.enames[curr-1]);
 193 @@ -599,6 +690,9 @@ namespace netgen
 194      //                         << " p1 " << mesh.LineSegment(i)[0] << " p2 " << mesh.LineSegment(i)[1] << endl;
 195      // exit(10);
 196
 197 +    for (int j = 1; j <= mesh.GetNP(); j++) // support SALOME fuse edges: set level to zero
 198 +      mesh.Point(j) = MeshPoint( (Point<3>&) mesh.Point(j) );
 199 +
 200      mesh.CalcSurfacesOfNode();
 201      multithread.task = savetask;
 202    }
 203 @@ -652,7 +746,7 @@ namespace netgen
 204
 205          Box<3> bb = geom.GetBoundingBox();
 206
 207 -        int projecttype = PLANESPACE;
 208 +        //int projecttype = PLANESPACE;
 209
 210          static Timer tinit("init");
 211          tinit.Start();
 212 diff --git a/libsrc/occ/occmeshsurf.cpp b/libsrc/occ/occmeshsurf.cpp
 213 index 7fa5d237..9e05fd95 100644
 214 --- a/libsrc/occ/occmeshsurf.cpp
 215 +++ b/libsrc/occ/occmeshsurf.cpp
 216 @@ -173,50 +173,6 @@ namespace netgen
 217         gp_Vec du, dv;
 218         occface->D1 (geominfo1.u, geominfo1.v, pnt, du, dv);
 219
 220 -        // static Timer t("occ-defintangplane calculations");
 221 -        // RegionTimer reg(t);
 222 -
 223 -        Mat<3,2> D1_;
 224 -       D1_(0,0) = du.X(); D1_(1,0) = du.Y(); D1_(2,0) = du.Z();
 225 -       D1_(0,1) = dv.X(); D1_(1,1) = dv.Y(); D1_(2,1) = dv.Z();
 226 -        auto D1T_ = Trans(D1_);
 227 -       auto D1TD1_ = D1T_*D1_;
 228 -       if (Det (D1TD1_) == 0) throw SingularMatrixException();
 229 -        Mat<2,2> DDTinv_;
 230 -        CalcInverse (D1TD1_, DDTinv_);
 231 -
 232 -        Mat<3,2> Y_;
 233 -       Vec<3> y1_ = (ap2-ap1).Normalize();
 234 -       Vec<3> y2_ = Cross(n, y1_).Normalize();
 235 -       for (int i = 0; i < 3; i++)
 236 -         {
 237 -           Y_(i,0) = y1_(i);
 238 -           Y_(i,1) = y2_(i);
 239 -         }
 240 -
 241 -        auto A_ = DDTinv_ * D1T_ * Y_;
 242 -       Mat<2,2> Ainv_;
 243 -       if (Det(A_) == 0) throw SingularMatrixException();
 244 -       CalcInverse (A_, Ainv_);
 245 -
 246 -       Vec<2> temp_ = Ainv_ * (psp2-psp1);
 247 -       double r_ = temp_.Length();
 248 -        Mat<2,2> R_;
 249 -        R_(0,0) = temp_(0)/r_;
 250 -        R_(1,0) = temp_(1)/r_;
 251 -        R_(0,1) = -R_(1,0);
 252 -        R_(1,1) = R_(0,0);
 253 -
 254 -        A_ = A_ * R_;
 255 -        Ainv_ = Trans(R_) * Ainv_;
 256 -
 257 -        Amat = A_;
 258 -        Amatinv = Ainv_;
 259 -
 260 -       // temp = Amatinv * (psp2-psp1);
 261 -
 262 -
 263 -#ifdef OLD
 264         DenseMatrix D1(3,2), D1T(2,3), DDTinv(2,2);
 265         D1(0,0) = du.X(); D1(1,0) = du.Y(); D1(2,0) = du.Z();
 266         D1(0,1) = dv.X(); D1(1,1) = dv.Y(); D1(2,1) = dv.Z();
 267 @@ -289,8 +245,7 @@ namespace netgen
 268             }
 269          // cout << "=?= Ainv = " << endl << Ainv << endl;
 270         temp = Amatinv * (psp2-psp1);
 271 -        cout << " =?= Amatinv = " << Amatinv << endl;
 272 -#endif
 273 +        // cout << " =?= Amatinv = " << Amatinv << endl;
 274        };
 275
 276    }
 277 @@ -380,6 +335,7 @@ namespace netgen
 278
 279      double u = gi.u;
 280      double v = gi.v;
 281 +    if ( 0 ) { // exclude Newton's method
 282      gp_Pnt x = occface->Value (u,v);
 283
 284      if (p.SquareDistance(x) <= sqr(PROJECTION_TOLERANCE)) return;
 285 @@ -429,7 +385,7 @@ namespace netgen
 286            }
 287        }
 288      while (count < 20);
 289 -
 290 +    }
 291
 292      // Newton did not converge, use OCC projection
 293