+@@ -474,20 +494,104 @@
+
+ if (!merge_solids)
+ {
+- pnums[0] = geom.vmap.FindIndex (TopExp::FirstVertex (edge));
+- pnums[pnums.Size()-1] = geom.vmap.FindIndex (TopExp::LastVertex (edge));
++ //pnums[0] = geom.vmap.FindIndex (TopExp::FirstVertex (edge));
++ //pnums[pnums.Size()-1] = geom.vmap.FindIndex (TopExp::LastVertex (edge));
++ MeshPoint dfltP ( Point<3> ( 0, 0, 0 ) );
++ int *ipp[] = { &pnums[0], &pnums[pnums.Size()-1] };
++ 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)
++ {
++ int &ip = *ipp[i];
++ ip = geom.vmap.FindIndex ( v[i] );
++ if ( ip == 0 || ip > nvertices )
++ {
++ int iv = ip;
++ if ( ip == 0 )
++ ip = iv = 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;
++ }
++ }