]> SALOME platform Git repositories - tools/sat_salome.git/commitdiff
Salome HOME
CRLF LF netgen62ForSalome.patch
authorNabil Ghodbane <nabil.ghodbane@cea.fr>
Mon, 13 Sep 2021 07:33:16 +0000 (09:33 +0200)
committerNabil Ghodbane <nabil.ghodbane@cea.fr>
Fri, 1 Oct 2021 10:22:21 +0000 (12:22 +0200)
products/patches/netgen62ForSalome.patch

index ba90e615d493527acb98610bce185926e7d65695..20c187b9755e250cfdfbbe7fc623895830ca74f4 100644 (file)
@@ -3,212 +3,212 @@ index 314d405a..7c912bcc 100644
 --- a/libsrc/occ/occgenmesh.cpp
 +++ b/libsrc/occ/occgenmesh.cpp
 @@ -49,7 +49,7 @@ namespace netgen
-   double ComputeH (double kappa, const MeshingParameters & mparam)
-   {
-     kappa *= mparam.curvaturesafety;
--    /*
-+    /**/
-     double hret;
-     if (mparam.maxh * kappa < 1)
+   double ComputeH (double kappa, const MeshingParameters & mparam)\r
+   {\r
+     kappa *= mparam.curvaturesafety;\r
+-    /*\r
++    /**/\r
+     double hret;\r
\r
+     if (mparam.maxh * kappa < 1)\r
 @@ -61,7 +61,7 @@ namespace netgen
-       hret = mparam.maxh;
-     return hret;
--    */
-+    /**/
-     // return min(mparam.maxh, 1/kappa);
-     return (mparam.maxh*kappa < 1) ? mparam.maxh : 1/kappa;
-   }
+       hret = mparam.maxh;\r
\r
+     return hret;\r
+-    */\r
++    /**/\r
+     // return min(mparam.maxh, 1/kappa);\r
+     return (mparam.maxh*kappa < 1) ? mparam.maxh : 1/kappa;\r
+   }\r
 @@ -253,10 +253,12 @@ namespace netgen
-       {
-         oldpnt = pnt;
-         pnt = c->Value(s0+(i/double(DIVIDEEDGESECTIONS))*(s1-s0));
-+        // !!! no more than 1 segment per <edge length>/DIVIDEEDGESECTIONS
-         hvalue[i] = hvalue[i-1] +
--          1.0/mesh.GetH(Point3d(pnt.X(), pnt.Y(), pnt.Z()))*
--          pnt.Distance(oldpnt);
-+          min( 1.0,
-+               1.0/mesh.GetH(Point3d(pnt.X(), pnt.Y(), pnt.Z()))*
-+               pnt.Distance(oldpnt));
-         //(*testout) << "mesh.GetH(Point3d(pnt.X(), pnt.Y(), pnt.Z())) " << mesh.GetH(Point3d(pnt.X(), pnt.Y(), pnt.Z()))
-         //       <<  " pnt.Distance(oldpnt) " << pnt.Distance(oldpnt) << endl;
+       {\r
+         oldpnt = pnt;\r
+         pnt = c->Value(s0+(i/double(DIVIDEEDGESECTIONS))*(s1-s0));\r
++        // !!! no more than 1 segment per <edge length>/DIVIDEEDGESECTIONS\r
+         hvalue[i] = hvalue[i-1] +\r
+-          1.0/mesh.GetH(Point3d(pnt.X(), pnt.Y(), pnt.Z()))*\r
+-          pnt.Distance(oldpnt);\r
++          min( 1.0,\r
++               1.0/mesh.GetH(Point3d(pnt.X(), pnt.Y(), pnt.Z()))*\r
++               pnt.Distance(oldpnt));\r
\r
+         //(*testout) << "mesh.GetH(Point3d(pnt.X(), pnt.Y(), pnt.Z())) " << mesh.GetH(Point3d(pnt.X(), pnt.Y(), pnt.Z()))\r
+         //       <<  " pnt.Distance(oldpnt) " << pnt.Distance(oldpnt) << endl;\r
\r
 @@ -299,7 +301,7 @@ namespace netgen
-         cout << "CORRECTED" << endl;
-         ps.SetSize (nsubedges-2);
-         params.SetSize (nsubedges);
--        params[nsubedges] = s1;
-+        params[nsubedges-1] = s1;
-       }
-   }
+         cout << "CORRECTED" << endl;\r
+         ps.SetSize (nsubedges-2);\r
+         params.SetSize (nsubedges);\r
+-        params[nsubedges] = s1;\r
++        params[nsubedges-1] = s1;\r
+       }\r
+   }\r
\r
 @@ -323,6 +325,8 @@ namespace netgen
-     double eps = 1e-6 * geom.GetBoundingBox().Diam();
-+    int first_vp = mesh.GetNP()+1;
-+
-     tsearch.Start();
-     for (int i = 1; i <= nvertices; i++)
-       {
\r
+     double eps = 1e-6 * geom.GetBoundingBox().Diam();\r
\r
++    int first_vp = mesh.GetNP()+1;\r
++\r
+     tsearch.Start();\r
+     for (int i = 1; i <= nvertices; i++)\r
+       {\r
 @@ -481,22 +485,107 @@ namespace netgen
-                 if (!merge_solids)
-                   {
--                    pnums[0] = geom.vmap.FindIndex (TopExp::FirstVertex (edge)) + PointIndex::BASE-1;
--                    pnums.Last() = geom.vmap.FindIndex (TopExp::LastVertex (edge)) + PointIndex::BASE-1;
-+                    //pnums[0] = geom.vmap.FindIndex (TopExp::FirstVertex (edge)) + PointIndex::BASE-1;
-+                    //pnums.Last() = geom.vmap.FindIndex (TopExp::LastVertex (edge)) + PointIndex::BASE-1;
-+                    const int dpi = PointIndex::BASE-1;
-+                    MeshPoint dfltP ( Point<3> ( 0, 0, 0 ) );
-+                    PointIndex *ipp[] = { &pnums[0], &pnums.Last() };
-+                    TopoDS_Iterator vIt( edge, false );
-+                    TopoDS_Vertex v[2];
-+                    v[0] = TopoDS::Vertex( vIt.Value() ); vIt.Next();
-+                    v[1] = TopoDS::Vertex( vIt.Value() );
-+                    if ( v[0].Orientation() == TopAbs_REVERSED )
-+                      std::swap( v[0], v[1] );
-+                    for ( int i = 0; i < 2; ++i)
-+                    {
-+                      PointIndex &ip = *ipp[i];
-+                      ip = geom.vmap.FindIndex ( v[i] ) + dpi;
-+                      if ( ip == dpi || ip > nvertices + dpi )
-+                      {
-+                        PointIndex iv = ip;
-+                        if ( ip == dpi )
-+                          ip = iv = const_cast<OCCGeometry&>(geom).vmap.Add( v[i] );
-+                        gp_Pnt pnt = BRep_Tool::Pnt( v[i] );
-+                        MeshPoint mp( Point<3>(pnt.X(), pnt.Y(), pnt.Z()) );
-+                        for (PointIndex pi = 1; pi < first_vp; pi++)
-+                          if ( Dist2 (mesh.Point(pi), Point<3>(mp)) < 1e-100 )
-+                          {
-+                            ip = pi;
-+                            if ( mesh.Point(ip).GetLayer() != dfltP.GetLayer() && mesh.Point(ip).GetLayer() != iv )
-+                              continue;
-+                            if ( mesh.Point(ip).GetLayer() == dfltP.GetLayer())
-+                              mesh.Point(ip) = MeshPoint( mesh.Point(ip), iv );
-+                            break;
-+                          }
-+                      }
-+                      else
-+                      {
-+                        ip += first_vp - 1;
-+                      }
-+                    }
-                   }
-                 else
-                   {
--                    Point<3> fp = occ2ng (BRep_Tool::Pnt (TopExp::FirstVertex (edge)));
--                    Point<3> lp = occ2ng (BRep_Tool::Pnt (TopExp::LastVertex (edge)));
-+                    TopoDS_Iterator vIt( edge, false );
-+                    TopoDS_Vertex v1 = TopoDS::Vertex( vIt.Value() ); vIt.Next();
-+                    TopoDS_Vertex v2 = TopoDS::Vertex( vIt.Value() );
-+                    if ( v1.Orientation() == TopAbs_REVERSED )
-+                      std::swap( v1, v2 );
-+                    const bool isClosedEdge = v1.IsSame( v2 );
-+
-+                    Point<3> fp = occ2ng (BRep_Tool::Pnt (v1));
-+                    Point<3> lp = occ2ng (BRep_Tool::Pnt (v2));
-+                    double tol2 = std::min( eps*eps, 1e-6 * Dist2( fp, lp ));
-+                    if ( isClosedEdge )
-+                      tol2 = BRep_Tool::Tolerance( v1 ) * BRep_Tool::Tolerance( v1 );
-                     pnums[0] = PointIndex::INVALID;
-                     pnums.Last() = PointIndex::INVALID;
-                     for (PointIndex pi : vertexrange)
-                       {
--                        if (Dist2 (mesh[pi], fp) < eps*eps) pnums[0] = pi;
--                        if (Dist2 (mesh[pi], lp) < eps*eps) pnums.Last() = pi;
-+                        if (Dist2 (mesh[pi], fp) < tol2) pnums[0] = pi;
-+                        if (Dist2 (mesh[pi], lp) < tol2) pnums.Last() = pi;
-                       }
--                  }
-+                    if (( isClosedEdge && pnums[0] != pnums.Last() ) ||
-+                        ( !isClosedEdge && pnums[0] == pnums.Last() ))
-+                      pnums[0] = pnums.Last() = PointIndex::INVALID;
-+                    if ( pnums[0] < PointIndex::BASE || pnums.Last() < PointIndex::BASE )
-+                    {
-+                      // take into account a possible large gap between a vertex and an edge curve
-+                      // end and a large vertex tolerance covering the whole edge
-+                      if ( pnums[0] < PointIndex::BASE )
-+                      {
-+                        double tol = BRep_Tool::Tolerance( v1 );
-+                        for (PointIndex pi = 1; pi < first_ep; pi++)
-+                          if (pi != pnums.Last() && Dist2 (mesh[pi], fp) < 2*tol*tol)
-+                            pnums[0] = pi;
-+
-+                        if ( pnums[0] < PointIndex::BASE )
-+                          pnums[0] = first_ep-1 - nvertices + geom.vmap.FindIndex ( v1 );
-+                      }
-+                      if ( isClosedEdge )
-+                      {
-+                        pnums.Last() = pnums[0];
-+                      }
-+                      else
-+                      {
-+                        if ( pnums.Last() < PointIndex::BASE )
-+                        {
-+                          double tol = BRep_Tool::Tolerance( v2 );
-+                          for (PointIndex pi = 1; pi < first_ep; pi++)
-+                            if (pi != pnums[0] && Dist2 (mesh[pi], lp) < 2*tol*tol)
-+                              pnums.Last() = pi;
-+
-+                          if ( pnums.Last() < PointIndex::BASE )
-+                            pnums.Last() = first_ep-1-nvertices + geom.vmap.FindIndex ( v2 );
-+                        }
-+
-+                        if ( Dist2( fp, mesh[PointIndex(pnums[0])]) >
-+                             Dist2( lp, mesh[PointIndex(pnums.Last())]))
-+                          std::swap( pnums[0], pnums.Last() );
-+                      }
-+                    }
-+                 }
-                 for (size_t i = 1; i <= mp.Size(); i++)
-                   {
\r
+                 if (!merge_solids)\r
+                   {\r
+-                    pnums[0] = geom.vmap.FindIndex (TopExp::FirstVertex (edge)) + PointIndex::BASE-1;\r
+-                    pnums.Last() = geom.vmap.FindIndex (TopExp::LastVertex (edge)) + PointIndex::BASE-1;\r
++                    //pnums[0] = geom.vmap.FindIndex (TopExp::FirstVertex (edge)) + PointIndex::BASE-1;\r
++                    //pnums.Last() = geom.vmap.FindIndex (TopExp::LastVertex (edge)) + PointIndex::BASE-1;\r
++                    const int dpi = PointIndex::BASE-1;\r
++                    MeshPoint dfltP ( Point<3> ( 0, 0, 0 ) );\r
++                    PointIndex *ipp[] = { &pnums[0], &pnums.Last() };\r
++                    TopoDS_Iterator vIt( edge, false );\r
++                    TopoDS_Vertex v[2];\r
++                    v[0] = TopoDS::Vertex( vIt.Value() ); vIt.Next();\r
++                    v[1] = TopoDS::Vertex( vIt.Value() );\r
++                    if ( v[0].Orientation() == TopAbs_REVERSED )\r
++                      std::swap( v[0], v[1] );\r
++                    for ( int i = 0; i < 2; ++i)\r
++                    {\r
++                      PointIndex &ip = *ipp[i];\r
++                      ip = geom.vmap.FindIndex ( v[i] ) + dpi;\r
++                      if ( ip == dpi || ip > nvertices + dpi )\r
++                      {\r
++                        PointIndex iv = ip;\r
++                        if ( ip == dpi )\r
++                          ip = iv = const_cast<OCCGeometry&>(geom).vmap.Add( v[i] );\r
++                        gp_Pnt pnt = BRep_Tool::Pnt( v[i] );\r
++                        MeshPoint mp( Point<3>(pnt.X(), pnt.Y(), pnt.Z()) );\r
++                        for (PointIndex pi = 1; pi < first_vp; pi++)\r
++                          if ( Dist2 (mesh.Point(pi), Point<3>(mp)) < 1e-100 )\r
++                          {\r
++                            ip = pi;\r
++                            if ( mesh.Point(ip).GetLayer() != dfltP.GetLayer() && mesh.Point(ip).GetLayer() != iv )\r
++                              continue;\r
++                            if ( mesh.Point(ip).GetLayer() == dfltP.GetLayer())\r
++                              mesh.Point(ip) = MeshPoint( mesh.Point(ip), iv );\r
++                            break;\r
++                          }\r
++                      }\r
++                      else\r
++                      {\r
++                        ip += first_vp - 1;\r
++                      }\r
++                    }\r
+                   }\r
+                 else\r
+                   {\r
+-                    Point<3> fp = occ2ng (BRep_Tool::Pnt (TopExp::FirstVertex (edge)));\r
+-                    Point<3> lp = occ2ng (BRep_Tool::Pnt (TopExp::LastVertex (edge)));\r
++                    TopoDS_Iterator vIt( edge, false );\r
++                    TopoDS_Vertex v1 = TopoDS::Vertex( vIt.Value() ); vIt.Next();\r
++                    TopoDS_Vertex v2 = TopoDS::Vertex( vIt.Value() );\r
++                    if ( v1.Orientation() == TopAbs_REVERSED )\r
++                      std::swap( v1, v2 );\r
++                    const bool isClosedEdge = v1.IsSame( v2 );\r
++\r
++                    Point<3> fp = occ2ng (BRep_Tool::Pnt (v1));\r
++                    Point<3> lp = occ2ng (BRep_Tool::Pnt (v2));\r
++                    double tol2 = std::min( eps*eps, 1e-6 * Dist2( fp, lp ));\r
++                    if ( isClosedEdge )\r
++                      tol2 = BRep_Tool::Tolerance( v1 ) * BRep_Tool::Tolerance( v1 );\r
\r
+                     pnums[0] = PointIndex::INVALID;\r
+                     pnums.Last() = PointIndex::INVALID;\r
+                     for (PointIndex pi : vertexrange)\r
+                       {\r
+-                        if (Dist2 (mesh[pi], fp) < eps*eps) pnums[0] = pi;\r
+-                        if (Dist2 (mesh[pi], lp) < eps*eps) pnums.Last() = pi;\r
++                        if (Dist2 (mesh[pi], fp) < tol2) pnums[0] = pi;\r
++                        if (Dist2 (mesh[pi], lp) < tol2) pnums.Last() = pi;\r
+                       }\r
+-                  }\r
++                    if (( isClosedEdge && pnums[0] != pnums.Last() ) ||\r
++                        ( !isClosedEdge && pnums[0] == pnums.Last() ))\r
++                      pnums[0] = pnums.Last() = PointIndex::INVALID;\r
++                    if ( pnums[0] < PointIndex::BASE || pnums.Last() < PointIndex::BASE )\r
++                    {\r
++                      // take into account a possible large gap between a vertex and an edge curve\r
++                      // end and a large vertex tolerance covering the whole edge\r
++                      if ( pnums[0] < PointIndex::BASE )\r
++                      {\r
++                        double tol = BRep_Tool::Tolerance( v1 );\r
++                        for (PointIndex pi = 1; pi < first_ep; pi++)\r
++                          if (pi != pnums.Last() && Dist2 (mesh[pi], fp) < 2*tol*tol)\r
++                            pnums[0] = pi;\r
++\r
++                        if ( pnums[0] < PointIndex::BASE )\r
++                          pnums[0] = first_ep-1 - nvertices + geom.vmap.FindIndex ( v1 );\r
++                      }\r
++                      if ( isClosedEdge )\r
++                      {\r
++                        pnums.Last() = pnums[0];\r
++                      }\r
++                      else\r
++                      {\r
++                        if ( pnums.Last() < PointIndex::BASE )\r
++                        {\r
++                          double tol = BRep_Tool::Tolerance( v2 );\r
++                          for (PointIndex pi = 1; pi < first_ep; pi++)\r
++                            if (pi != pnums[0] && Dist2 (mesh[pi], lp) < 2*tol*tol)\r
++                              pnums.Last() = pi;\r
++\r
++                          if ( pnums.Last() < PointIndex::BASE )\r
++                            pnums.Last() = first_ep-1-nvertices + geom.vmap.FindIndex ( v2 );\r
++                        }\r
++\r
++                        if ( Dist2( fp, mesh[PointIndex(pnums[0])]) >\r
++                             Dist2( lp, mesh[PointIndex(pnums.Last())]))\r
++                          std::swap( pnums[0], pnums.Last() );\r
++                      }\r
++                    }\r
++                 }\r
\r
+                 for (size_t i = 1; i <= mp.Size(); i++)\r
+                   {\r
 @@ -505,17 +594,19 @@ namespace netgen
-                     // for (PointIndex j = first_ep; j < mesh.Points().End(); j++)
-                     for (PointIndex j = first_ep; j < *mesh.Points().Range().end(); j++)
-+                    {
-+                      if (!merge_solids && mesh.Point(j).GetLayer() != geomedgenr ) continue;
-                       if ((mesh.Point(j)-Point<3>(mp[i-1])).Length() < eps)
-                         {
-                           exists = true;
-                           pnums[i] = j;
-                           break;
-                         }
--                      
-+                    }
-                     tsearch.Stop();
-                     
-                     if (!exists)
--                      pnums[i] = mesh.AddPoint (mp[i-1]);
-+                      pnums[i] = mesh.AddPoint (mp[i-1], geomedgenr);
-                   }
-                 if(geom.enames.Size() && geom.enames[curr-1] != "")
-                   mesh.SetCD2Name(geomedgenr, geom.enames[curr-1]);
\r
+                     // for (PointIndex j = first_ep; j < mesh.Points().End(); j++)\r
+                     for (PointIndex j = first_ep; j < *mesh.Points().Range().end(); j++)\r
++                    {\r
++                      if (!merge_solids && mesh.Point(j).GetLayer() != geomedgenr ) continue;\r
+                       if ((mesh.Point(j)-Point<3>(mp[i-1])).Length() < eps)\r
+                         {\r
+                           exists = true;\r
+                           pnums[i] = j;\r
+                           break;\r
+                         }\r
+-                      \r
++                    }\r
+                     tsearch.Stop();\r
+                     \r
+                     if (!exists)\r
+-                      pnums[i] = mesh.AddPoint (mp[i-1]);\r
++                      pnums[i] = mesh.AddPoint (mp[i-1], geomedgenr);\r
+                   }\r
+                 if(geom.enames.Size() && geom.enames[curr-1] != "")\r
+                   mesh.SetCD2Name(geomedgenr, geom.enames[curr-1]);\r
 @@ -599,6 +690,9 @@ namespace netgen
-     //                                << " p1 " << mesh.LineSegment(i)[0] << " p2 " << mesh.LineSegment(i)[1] << endl;
-     //        exit(10);
-+    for (int j = 1; j <= mesh.GetNP(); j++) // support SALOME fuse edges: set level to zero
-+      mesh.Point(j) = MeshPoint( (Point<3>&) mesh.Point(j) );
-+
-     mesh.CalcSurfacesOfNode();
-     multithread.task = savetask;
-   }
+     //                                << " p1 " << mesh.LineSegment(i)[0] << " p2 " << mesh.LineSegment(i)[1] << endl;\r
+     //        exit(10);\r
\r
++    for (int j = 1; j <= mesh.GetNP(); j++) // support SALOME fuse edges: set level to zero\r
++      mesh.Point(j) = MeshPoint( (Point<3>&) mesh.Point(j) );\r
++\r
+     mesh.CalcSurfacesOfNode();\r
+     multithread.task = savetask;\r
+   }\r
 @@ -652,7 +746,7 @@ namespace netgen
-         Box<3> bb = geom.GetBoundingBox();
--        int projecttype = PLANESPACE;
-+        //int projecttype = PLANESPACE;
-         
-         static Timer tinit("init");
-         tinit.Start();
\r
+         Box<3> bb = geom.GetBoundingBox();\r
\r
+-        int projecttype = PLANESPACE;\r
++        //int projecttype = PLANESPACE;\r
+         \r
+         static Timer tinit("init");\r
+         tinit.Start();\r
 diff --git a/libsrc/occ/occmeshsurf.cpp b/libsrc/occ/occmeshsurf.cpp
 index 7fa5d237..9e05fd95 100644
 --- a/libsrc/occ/occmeshsurf.cpp