--- 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
