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
8 kappa *= mparam.curvaturesafety;
\r
13 if (mparam.maxh * kappa < 1)
\r
14 @@ -61,7 +61,7 @@ namespace netgen
20 // return min(mparam.maxh, 1/kappa);
\r
21 return (mparam.maxh*kappa < 1) ? mparam.maxh : 1/kappa;
\r
23 @@ -253,10 +253,12 @@ namespace netgen
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
32 + 1.0/mesh.GetH(Point3d(pnt.X(), pnt.Y(), pnt.Z()))*
\r
33 + pnt.Distance(oldpnt));
\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
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
47 @@ -323,6 +325,8 @@ namespace netgen
49 double eps = 1e-6 * geom.GetBoundingBox().Diam();
\r
51 + int first_vp = mesh.GetNP()+1;
\r
54 for (int i = 1; i <= nvertices; i++)
\r
56 @@ -481,22 +485,107 @@ namespace netgen
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
75 + PointIndex &ip = *ipp[i];
\r
76 + ip = geom.vmap.FindIndex ( v[i] ) + dpi;
\r
77 + if ( ip == dpi || ip > nvertices + dpi )
\r
79 + PointIndex iv = ip;
\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
88 + if ( mesh.Point(ip).GetLayer() != dfltP.GetLayer() && mesh.Point(ip).GetLayer() != iv )
\r
90 + if ( mesh.Point(ip).GetLayer() == dfltP.GetLayer())
\r
91 + mesh.Point(ip) = MeshPoint( mesh.Point(ip), iv );
\r
97 + ip += first_vp - 1;
\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
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
118 pnums[0] = PointIndex::INVALID;
\r
119 pnums.Last() = PointIndex::INVALID;
\r
120 for (PointIndex pi : vertexrange)
\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
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
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
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
142 + if ( pnums[0] < PointIndex::BASE )
\r
143 + pnums[0] = first_ep-1 - nvertices + geom.vmap.FindIndex ( v1 );
\r
145 + if ( isClosedEdge )
\r
147 + pnums.Last() = pnums[0];
\r
151 + if ( pnums.Last() < PointIndex::BASE )
\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
158 + if ( pnums.Last() < PointIndex::BASE )
\r
159 + pnums.Last() = first_ep-1-nvertices + geom.vmap.FindIndex ( v2 );
\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
169 for (size_t i = 1; i <= mp.Size(); i++)
\r
171 @@ -505,17 +594,19 @@ namespace netgen
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
176 + if (!merge_solids && mesh.Point(j).GetLayer() != geomedgenr ) continue;
\r
177 if ((mesh.Point(j)-Point<3>(mp[i-1])).Length() < eps)
\r
188 - pnums[i] = mesh.AddPoint (mp[i-1]);
\r
189 + pnums[i] = mesh.AddPoint (mp[i-1], geomedgenr);
\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
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
200 mesh.CalcSurfacesOfNode();
\r
201 multithread.task = savetask;
\r
203 @@ -652,7 +746,7 @@ namespace netgen
205 Box<3> bb = geom.GetBoundingBox();
\r
207 - int projecttype = PLANESPACE;
\r
208 + //int projecttype = PLANESPACE;
\r
210 static Timer tinit("init");
\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
218 occface->D1 (geominfo1.u, geominfo1.v, pnt, du, dv);
220 - // static Timer t("occ-defintangplane calculations");
221 - // RegionTimer reg(t);
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();
230 - CalcInverse (D1TD1_, DDTinv_);
233 - Vec<3> y1_ = (ap2-ap1).Normalize();
234 - Vec<3> y2_ = Cross(n, y1_).Normalize();
235 - for (int i = 0; i < 3; i++)
241 - auto A_ = DDTinv_ * D1T_ * Y_;
243 - if (Det(A_) == 0) throw SingularMatrixException();
244 - CalcInverse (A_, Ainv_);
246 - Vec<2> temp_ = Ainv_ * (psp2-psp1);
247 - double r_ = temp_.Length();
249 - R_(0,0) = temp_(0)/r_;
250 - R_(1,0) = temp_(1)/r_;
251 - R_(0,1) = -R_(1,0);
255 - Ainv_ = Trans(R_) * Ainv_;
260 - // temp = Amatinv * (psp2-psp1);
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
269 // cout << "=?= Ainv = " << endl << Ainv << endl;
270 temp = Amatinv * (psp2-psp1);
271 - cout << " =?= Amatinv = " << Amatinv << endl;
273 + // cout << " =?= Amatinv = " << Amatinv << endl;
277 @@ -380,6 +335,7 @@ namespace netgen
281 + if ( 0 ) { // exclude Newton's method
282 gp_Pnt x = occface->Value (u,v);
284 if (p.SquareDistance(x) <= sqr(PROJECTION_TOLERANCE)) return;
285 @@ -429,7 +385,7 @@ namespace netgen
292 // Newton did not converge, use OCC projection