src/NETGEN/netgen62ForSalome.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)\r
   7    {\r
   8      kappa *= mparam.curvaturesafety;\r
   9 -    /*\r
  10 +    /**/\r
  11      double hret;\r
  12  \r
  13      if (mparam.maxh * kappa < 1)\r
  14 @@ -61,7 +61,7 @@ namespace netgen
  15        hret = mparam.maxh;\r
  16  \r
  17      return hret;\r
  18 -    */\r
  19 +    /**/\r
  20      // return min(mparam.maxh, 1/kappa);\r
  21      return (mparam.maxh*kappa < 1) ? mparam.maxh : 1/kappa;\r
  22    }\r
  23 @@ -253,10 +253,12 @@ namespace netgen
  24        {\r
  25          oldpnt = pnt;\r
  26          pnt = c->Value(s0+(i/double(DIVIDEEDGESECTIONS))*(s1-s0));\r
  27 +        // !!! no more than 1 segment per <edge length>/DIVIDEEDGESECTIONS\r
  28          hvalue[i] = hvalue[i-1] +\r
  29 -          1.0/mesh.GetH(Point3d(pnt.X(), pnt.Y(), pnt.Z()))*\r
  30 -          pnt.Distance(oldpnt);\r
  31 +          min( 1.0,\r
  32 +               1.0/mesh.GetH(Point3d(pnt.X(), pnt.Y(), pnt.Z()))*\r
  33 +               pnt.Distance(oldpnt));\r
  34  \r
  35          //(*testout) << "mesh.GetH(Point3d(pnt.X(), pnt.Y(), pnt.Z())) " << mesh.GetH(Point3d(pnt.X(), pnt.Y(), pnt.Z()))\r
  36          //        <<  " pnt.Distance(oldpnt) " << pnt.Distance(oldpnt) << endl;\r
  37  \r
  38 @@ -299,7 +301,7 @@ namespace netgen
  39          cout << "CORRECTED" << endl;\r
  40          ps.SetSize (nsubedges-2);\r
  41          params.SetSize (nsubedges);\r
  42 -        params[nsubedges] = s1;\r
  43 +        params[nsubedges-1] = s1;\r
  44        }\r
  45    }\r
  46  \r
  47 @@ -323,6 +325,8 @@ namespace netgen
  48  \r
  49      double eps = 1e-6 * geom.GetBoundingBox().Diam();\r
  50  \r
  51 +    int first_vp = mesh.GetNP()+1;\r
  52 +\r
  53      tsearch.Start();\r
  54      for (int i = 1; i <= nvertices; i++)\r
  55        {\r
  56 @@ -481,22 +485,107 @@ namespace netgen
  57  \r
  58                  if (!merge_solids)\r
  59                    {\r
  60 -                    pnums[0] = geom.vmap.FindIndex (TopExp::FirstVertex (edge)) + PointIndex::BASE-1;\r
  61 -                    pnums.Last() = geom.vmap.FindIndex (TopExp::LastVertex (edge)) + PointIndex::BASE-1;\r
  62 +                    //pnums[0] = geom.vmap.FindIndex (TopExp::FirstVertex (edge)) + PointIndex::BASE-1;\r
  63 +                    //pnums.Last() = geom.vmap.FindIndex (TopExp::LastVertex (edge)) + PointIndex::BASE-1;\r
  64 +                    const int dpi = PointIndex::BASE-1;\r
  65 +                    MeshPoint dfltP ( Point<3> ( 0, 0, 0 ) );\r
  66 +                    PointIndex *ipp[] = { &pnums[0], &pnums.Last() };\r
  67 +                    TopoDS_Iterator vIt( edge, false );\r
  68 +                    TopoDS_Vertex v[2];\r
  69 +                    v[0] = TopoDS::Vertex( vIt.Value() ); vIt.Next();\r
  70 +                    v[1] = TopoDS::Vertex( vIt.Value() );\r
  71 +                    if ( v[0].Orientation() == TopAbs_REVERSED )\r
  72 +                      std::swap( v[0], v[1] );\r
  73 +                    for ( int i = 0; i < 2; ++i)\r
  74 +                    {\r
  75 +                      PointIndex &ip = *ipp[i];\r
  76 +                      ip = geom.vmap.FindIndex ( v[i] ) + dpi;\r
  77 +                      if ( ip == dpi || ip > nvertices + dpi )\r
  78 +                      {\r
  79 +                        PointIndex iv = ip;\r
  80 +                        if ( ip == dpi )\r
  81 +                          ip = iv = const_cast<OCCGeometry&>(geom).vmap.Add( v[i] );\r
  82 +                        gp_Pnt pnt = BRep_Tool::Pnt( v[i] );\r
  83 +                        MeshPoint mp( Point<3>(pnt.X(), pnt.Y(), pnt.Z()) );\r
  84 +                        for (PointIndex pi = 1; pi < first_vp; pi++)\r
  85 +                          if ( Dist2 (mesh.Point(pi), Point<3>(mp)) < 1e-100 )\r
  86 +                          {\r
  87 +                            ip = pi;\r
  88 +                            if ( mesh.Point(ip).GetLayer() != dfltP.GetLayer() && mesh.Point(ip).GetLayer() != iv )\r
  89 +                              continue;\r
  90 +                            if ( mesh.Point(ip).GetLayer() == dfltP.GetLayer())\r
  91 +                              mesh.Point(ip) = MeshPoint( mesh.Point(ip), iv );\r
  92 +                            break;\r
  93 +                          }\r
  94 +                      }\r
  95 +                      else\r
  96 +                      {\r
  97 +                        ip += first_vp - 1;\r
  98 +                      }\r
  99 +                    }\r
 100                    }\r
 101                  else\r
 102                    {\r
 103 -                    Point<3> fp = occ2ng (BRep_Tool::Pnt (TopExp::FirstVertex (edge)));\r
 104 -                    Point<3> lp = occ2ng (BRep_Tool::Pnt (TopExp::LastVertex (edge)));\r
 105 +                    TopoDS_Iterator vIt( edge, false );\r
 106 +                    TopoDS_Vertex v1 = TopoDS::Vertex( vIt.Value() ); vIt.Next();\r
 107 +                    TopoDS_Vertex v2 = TopoDS::Vertex( vIt.Value() );\r
 108 +                    if ( v1.Orientation() == TopAbs_REVERSED )\r
 109 +                      std::swap( v1, v2 );\r
 110 +                    const bool isClosedEdge = v1.IsSame( v2 );\r
 111 +\r
 112 +                    Point<3> fp = occ2ng (BRep_Tool::Pnt (v1));\r
 113 +                    Point<3> lp = occ2ng (BRep_Tool::Pnt (v2));\r
 114 +                    double tol2 = std::min( eps*eps, 1e-6 * Dist2( fp, lp ));\r
 115 +                    if ( isClosedEdge )\r
 116 +                      tol2 = BRep_Tool::Tolerance( v1 ) * BRep_Tool::Tolerance( v1 );\r
 117  \r
 118                      pnums[0] = PointIndex::INVALID;\r
 119                      pnums.Last() = PointIndex::INVALID;\r
 120                      for (PointIndex pi : vertexrange)\r
 121                        {\r
 122 -                        if (Dist2 (mesh[pi], fp) < eps*eps) pnums[0] = pi;\r
 123 -                        if (Dist2 (mesh[pi], lp) < eps*eps) pnums.Last() = pi;\r
 124 +                        if (Dist2 (mesh[pi], fp) < tol2) pnums[0] = pi;\r
 125 +                        if (Dist2 (mesh[pi], lp) < tol2) pnums.Last() = pi;\r
 126                        }\r
 127 -                  }\r
 128 +                    if (( isClosedEdge && pnums[0] != pnums.Last() ) ||\r
 129 +                        ( !isClosedEdge && pnums[0] == pnums.Last() ))\r
 130 +                      pnums[0] = pnums.Last() = PointIndex::INVALID;\r
 131 +                    if ( pnums[0] < PointIndex::BASE || pnums.Last() < PointIndex::BASE )\r
 132 +                    {\r
 133 +                      // take into account a possible large gap between a vertex and an edge curve\r
 134 +                      // end and a large vertex tolerance covering the whole edge\r
 135 +                      if ( pnums[0] < PointIndex::BASE )\r
 136 +                      {\r
 137 +                        double tol = BRep_Tool::Tolerance( v1 );\r
 138 +                        for (PointIndex pi = 1; pi < first_ep; pi++)\r
 139 +                          if (pi != pnums.Last() && Dist2 (mesh[pi], fp) < 2*tol*tol)\r
 140 +                            pnums[0] = pi;\r
 141 +\r
 142 +                        if ( pnums[0] < PointIndex::BASE )\r
 143 +                          pnums[0] = first_ep-1 - nvertices + geom.vmap.FindIndex ( v1 );\r
 144 +                      }\r
 145 +                      if ( isClosedEdge )\r
 146 +                      {\r
 147 +                        pnums.Last() = pnums[0];\r
 148 +                      }\r
 149 +                      else\r
 150 +                      {\r
 151 +                        if ( pnums.Last() < PointIndex::BASE )\r
 152 +                        {\r
 153 +                          double tol = BRep_Tool::Tolerance( v2 );\r
 154 +                          for (PointIndex pi = 1; pi < first_ep; pi++)\r
 155 +                            if (pi != pnums[0] && Dist2 (mesh[pi], lp) < 2*tol*tol)\r
 156 +                              pnums.Last() = pi;\r
 157 +\r
 158 +                          if ( pnums.Last() < PointIndex::BASE )\r
 159 +                            pnums.Last() = first_ep-1-nvertices + geom.vmap.FindIndex ( v2 );\r
 160 +                        }\r
 161 +\r
 162 +                        if ( Dist2( fp, mesh[PointIndex(pnums[0])]) >\r
 163 +                             Dist2( lp, mesh[PointIndex(pnums.Last())]))\r
 164 +                          std::swap( pnums[0], pnums.Last() );\r
 165 +                      }\r
 166 +                    }\r
 167 +                 }\r
 168  \r
 169                  for (size_t i = 1; i <= mp.Size(); i++)\r
 170                    {\r
 171 @@ -505,17 +594,19 @@ namespace netgen
 172  \r
 173                      // for (PointIndex j = first_ep; j < mesh.Points().End(); j++)\r
 174                      for (PointIndex j = first_ep; j < *mesh.Points().Range().end(); j++)\r
 175 +                    {\r
 176 +                      if (!merge_solids && mesh.Point(j).GetLayer() != geomedgenr ) continue;\r
 177                        if ((mesh.Point(j)-Point<3>(mp[i-1])).Length() < eps)\r
 178                          {\r
 179                            exists = true;\r
 180                            pnums[i] = j;\r
 181                            break;\r
 182                          }\r
 183 -                      \r
 184 +                    }\r
 185                      tsearch.Stop();\r
 186                      \r
 187                      if (!exists)\r
 188 -                      pnums[i] = mesh.AddPoint (mp[i-1]);\r
 189 +                      pnums[i] = mesh.AddPoint (mp[i-1], geomedgenr);\r
 190                    }\r
 191                  if(geom.enames.Size() && geom.enames[curr-1] != "")\r
 192                    mesh.SetCD2Name(geomedgenr, geom.enames[curr-1]);\r
 193 @@ -599,6 +690,9 @@ namespace netgen
 194      //                         << " p1 " << mesh.LineSegment(i)[0] << " p2 " << mesh.LineSegment(i)[1] << endl;\r
 195      // exit(10);\r
 196  \r
 197 +    for (int j = 1; j <= mesh.GetNP(); j++) // support SALOME fuse edges: set level to zero\r
 198 +      mesh.Point(j) = MeshPoint( (Point<3>&) mesh.Point(j) );\r
 199 +\r
 200      mesh.CalcSurfacesOfNode();\r
 201      multithread.task = savetask;\r
 202    }\r
 203 @@ -652,7 +746,7 @@ namespace netgen
 204  \r
 205          Box<3> bb = geom.GetBoundingBox();\r
 206  \r
 207 -        int projecttype = PLANESPACE;\r
 208 +        //int projecttype = PLANESPACE;\r
 209          \r
 210          static Timer tinit("init");\r
 211          tinit.Start();\r
 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