-diff -Nru netgen-4.5_orig/libsrc/csg/meshsurf.cpp netgen-4.5_patch/libsrc/csg/meshsurf.cpp\r
---- netgen-4.5_orig/libsrc/csg/meshsurf.cpp 2006-02-14 10:54:35.000000000 +0300\r
-+++ netgen-4.5_patch/libsrc/csg/meshsurf.cpp 2006-10-25 16:05:59.000000000 +0400\r
-@@ -77,11 +77,12 @@\r
- }\r
- \r
- \r
--void MeshOptimize2dSurfaces :: ProjectPoint (INDEX surfind, Point3d & p) const\r
-+bool MeshOptimize2dSurfaces :: ProjectPoint (INDEX surfind, Point3d & p) const\r
- {\r
- Point<3> hp = p;\r
- geometry.GetSurface(surfind)->Project (hp);\r
- p = hp;\r
-+ return true;\r
- }\r
- \r
- void MeshOptimize2dSurfaces :: ProjectPoint2 (INDEX surfind, INDEX surfind2, \r
-diff -Nru netgen-4.5_orig/libsrc/csg/meshsurf.hpp netgen-4.5_patch/libsrc/csg/meshsurf.hpp\r
---- netgen-4.5_orig/libsrc/csg/meshsurf.hpp 2004-01-20 13:49:44.000000000 +0300\r
-+++ netgen-4.5_patch/libsrc/csg/meshsurf.hpp 2006-10-25 16:08:05.000000000 +0400\r
-@@ -45,7 +45,7 @@\r
- MeshOptimize2dSurfaces (const CSGeometry & ageometry); \r
- \r
- ///\r
-- virtual void ProjectPoint (INDEX surfind, Point3d & p) const;\r
-+ virtual bool ProjectPoint (INDEX surfind, Point3d & p) const;\r
- ///\r
- virtual void ProjectPoint2 (INDEX surfind, INDEX surfind2, Point3d & p) const;\r
- ///\r
-diff -Nru netgen-4.5_orig/libsrc/interface/Makefile netgen-4.5_patch/libsrc/interface/Makefile\r
---- netgen-4.5_orig/libsrc/interface/Makefile 2005-08-09 18:14:59.000000000 +0400\r
-+++ netgen-4.5_patch/libsrc/interface/Makefile 2006-04-27 13:12:54.000000000 +0400\r
-@@ -1,4 +1,5 @@\r
--src = nginterface.cpp writeuser.cpp writediffpack.cpp writeabaqus.cpp writefluent.cpp writepermas.cpp writetochnog.cpp writetecplot.cpp wuchemnitz.cpp writetochnog.cpp writefeap.cpp writeelmer.cpp writegmsh.cpp writejcm.cpp readuser.cpp importsolution.cpp \r
-+#src = nginterface.cpp writeuser.cpp writediffpack.cpp writeabaqus.cpp writefluent.cpp writepermas.cpp writetochnog.cpp writetecplot.cpp wuchemnitz.cpp writetochnog.cpp writefeap.cpp writeelmer.cpp writegmsh.cpp writejcm.cpp readuser.cpp importsolution.cpp\r
-+src = writeuser.cpp writediffpack.cpp writeabaqus.cpp writefluent.cpp writepermas.cpp writetochnog.cpp writetecplot.cpp wuchemnitz.cpp writetochnog.cpp writefeap.cpp writeelmer.cpp writegmsh.cpp writejcm.cpp readuser.cpp nglib.cpp ngnewdelete.cpp\r
- #\r
- lib = nginterface\r
- libpath = libsrc/interface\r
-diff -Nru netgen-4.5_orig/libsrc/interface/nglib.cpp netgen-4.5_patch/libsrc/interface/nglib.cpp\r
---- netgen-4.5_orig/libsrc/interface/nglib.cpp 2005-10-18 17:53:18.000000000 +0400\r
-+++ netgen-4.5_patch/libsrc/interface/nglib.cpp 2006-04-27 13:12:54.000000000 +0400\r
-@@ -56,7 +56,8 @@\r
- \r
- void Ng_Exit ()\r
- {\r
-- ;\r
-+ delete testout;\r
-+ testout = NULL;\r
- }\r
- \r
- \r
-diff -Nru netgen-4.5_orig/libsrc/makefile.inc netgen-4.5_patch/libsrc/makefile.inc\r
---- netgen-4.5_orig/libsrc/makefile.inc 2005-09-02 17:17:51.000000000 +0400\r
-+++ netgen-4.5_patch/libsrc/makefile.inc 2006-10-26 09:33:58.000000000 +0400\r
-@@ -8,17 +8,14 @@\r
- LIBSRC_DIR=$(CPP_DIR)/libsrc\r
- LIB_DIR=$(CPP_DIR)/lib/$(MACHINE)\r
- \r
--#OCC_DIR=../../occ\r
--#OCCINC_DIR=$(OCC_DIR)/inc\r
--#OCCLIB_DIR=$(OCC_DIR)/lib\r
--# OCC_DIR=/opt/OpenCASCADE5.2/ros\r
--# OCC_DIR=/home/joachim/download/occ/Linux\r
--# OCCINC_DIR=$(OCC_DIR)/inc -I$(OCC_DIR)/ros/inc\r
--# OCCLIB_DIR=$(OCC_DIR)/Linux/lib\r
-+OCC_DIR=$(CASROOT)\r
-+OCCINC_DIR=$(OCC_DIR)/inc\r
-+OCCLIB_DIR=$(OCC_DIR)/Linux/lib\r
- #\r
- include $(LIBSRC_DIR)/makefile.mach.$(MACHINE)\r
- #\r
--CPLUSPLUSFLAGS1 = -c -I$(LIBSRC_DIR)/include -I$(OCCINC_DIR) \r
-+CPLUSPLUSFLAGS1 = -c -I$(LIBSRC_DIR)/include -I$(OCCINC_DIR) \\r
-+ -DOCCGEOMETRY -DOCC52 -DHAVE_IOSTREAM -DHAVE_LIMITS_H\r
- #\r
- ARFLAGS = r\r
- #\r
-diff -Nru netgen-4.5_orig/libsrc/makefile.mach.LINUX netgen-4.5_patch/libsrc/makefile.mach.LINUX\r
---- netgen-4.5_orig/libsrc/makefile.mach.LINUX 2004-10-11 23:49:26.000000000 +0400\r
-+++ netgen-4.5_patch/libsrc/makefile.mach.LINUX 2006-04-27 13:12:54.000000000 +0400\r
-@@ -16,7 +16,7 @@\r
- #\r
- CFLAGS2 =\r
- \r
--CPLUSPLUSFLAGS2 = -O2 -I/usr/include/GL3.5 -DLINUX -DOPENGL \\r
-+CPLUSPLUSFLAGS2 = -O2 -I/usr/include/GL3.5 -DLINUX \\r
- -ftemplate-depth-99 -finline-limit=10000 \\r
- -Wdisabled-optimization -funroll-loops -DnoNGSOLVE\r
- \r
-diff -Nru netgen-4.5_orig/libsrc/meshing/improve2.cpp netgen-4.5_patch/libsrc/meshing/improve2.cpp\r
---- netgen-4.5_orig/libsrc/meshing/improve2.cpp 2006-01-11 18:08:19.000000000 +0300\r
-+++ netgen-4.5_patch/libsrc/meshing/improve2.cpp 2006-04-27 13:12:54.000000000 +0400\r
-@@ -4,7 +4,7 @@\r
- #include <opti.hpp>\r
- \r
- #ifndef SMALLLIB\r
--#include <visual.hpp>\r
-+//#include <visual.hpp>\r
- #endif\r
- \r
- namespace netgen\r
-diff -Nru netgen-4.5_orig/libsrc/meshing/improve2.hpp netgen-4.5_patch/libsrc/meshing/improve2.hpp\r
---- netgen-4.5_orig/libsrc/meshing/improve2.hpp 2004-10-12 23:22:55.000000000 +0400\r
-+++ netgen-4.5_patch/libsrc/meshing/improve2.hpp 2006-10-25 16:09:37.000000000 +0400\r
-@@ -32,17 +32,16 @@\r
- ///\r
- virtual void SelectSurfaceOfPoint (const Point3d & p,\r
- const PointGeomInfo & gi);\r
-- ///\r
-- virtual void ProjectPoint (INDEX /* surfind */, Point3d & /* p */) const { };\r
-+\r
-+ /// project point on surface, returns true if success\r
-+ virtual bool ProjectPoint (INDEX /* surfind */, Point3d & /* p */) const { return false; }\r
-+ /// fast project point on surface using point geom info of a neighboring point\r
-+ /// if gi.trignum != 0,\r
-+ /// returns true if success, gi is updated\r
-+ virtual bool ProjectPoint (INDEX surfind, Point3d & p, PointGeomInfo& gi) const\r
-+ { gi.trignum = 1; return ProjectPoint (surfind, p); }\r
- ///\r
- virtual void ProjectPoint2 (INDEX /* surfind */, INDEX /* surfind2 */, Point3d & /* p */) const { };\r
-- /// liefert zu einem 3d-Punkt die geominfo (Dreieck) und liefert 1, wenn erfolgreich, \r
-- /// 0, wenn nicht (Punkt ausserhalb von chart)\r
-- virtual int CalcPointGeomInfo(PointGeomInfo& gi, const Point3d& /*p3*/) const\r
-- { gi.trignum = 1; return 1;};\r
--\r
-- virtual int CalcPointGeomInfo(int /* surfind */, PointGeomInfo& gi, const Point3d& p3) const\r
-- { return CalcPointGeomInfo (gi, p3); }\r
- \r
- ///\r
- virtual void GetNormalVector(INDEX surfind, const Point3d & p, PointGeomInfo & gi, Vec3d & n) const;\r
-diff -Nru netgen-4.5_orig/libsrc/meshing/smoothing2.cpp netgen-4.5_patch/libsrc/meshing/smoothing2.cpp\r
---- netgen-4.5_orig/libsrc/meshing/smoothing2.cpp 2006-01-11 18:08:20.000000000 +0300\r
-+++ netgen-4.5_patch/libsrc/meshing/smoothing2.cpp 2006-10-25 16:10:46.000000000 +0400\r
-@@ -300,7 +300,7 @@\r
- double Opti2SurfaceMinFunction :: \r
- FuncGrad (const Vector & x, Vector & grad) const\r
- {\r
-- Vec3d n, vgrad;\r
-+ Vec3d vgrad;\r
- Point3d pp1;\r
- double g1x, g1y;\r
- double badness, hbadness;\r
-@@ -308,8 +308,6 @@\r
- vgrad = 0;\r
- badness = 0;\r
- \r
-- meshthis -> GetNormalVector (surfi, sp1, gi1, n);\r
--\r
- pp1 = sp1;\r
- pp1.Add2 (x.Get(1), t1, x.Get(2), t2);\r
- \r
-@@ -360,7 +358,7 @@\r
- double Opti2SurfaceMinFunction :: \r
- FuncDeriv (const Vector & x, const Vector & dir, double & deriv) const\r
- {\r
-- Vec3d n, vgrad;\r
-+ Vec3d vgrad;\r
- Point3d pp1;\r
- double g1x, g1y;\r
- double badness, hbadness;\r
-@@ -368,8 +366,6 @@\r
- vgrad = 0;\r
- badness = 0;\r
- \r
-- meshthis -> GetNormalVector (surfi, sp1, gi1, n);\r
--\r
- pp1 = sp1;\r
- pp1.Add2 (x.Get(1), t1, x.Get(2), t2);\r
- \r
-@@ -520,7 +516,7 @@\r
- // from 2d:\r
- \r
- int j, k, lpi, gpi;\r
-- Vec3d n, vgrad;\r
-+ Vec3d vgrad;\r
- Point3d pp1;\r
- Vec2d g1, vdir;\r
- double badness, hbadness, hbad, hderiv;\r
-@@ -528,8 +524,6 @@\r
- vgrad = 0;\r
- badness = 0;\r
- \r
-- meshthis -> GetNormalVector (surfi, sp1, gi1, n);\r
--\r
- pp1 = sp1;\r
- pp1.Add2 (x.Get(1), t1, x.Get(2), t2);\r
- \r
-@@ -593,7 +587,7 @@\r
- // from 2d:\r
- \r
- int j, k, lpi, gpi;\r
-- Vec3d n, vgrad;\r
-+ Vec3d vgrad;\r
- Point3d pp1;\r
- Vec2d g1, vdir;\r
- double badness, hbadness, hbad, hderiv;\r
-@@ -601,8 +595,6 @@\r
- vgrad = 0;\r
- badness = 0;\r
- \r
-- meshthis -> GetNormalVector (surfi, sp1, gi1, n);\r
--\r
- pp1 = sp1;\r
- pp1.Add2 (x.Get(1), t1, x.Get(2), t2);\r
- \r
-@@ -859,19 +851,21 @@\r
- locelements.SetSize(0);\r
- locrots.SetSize (0);\r
- lochs.SetSize (0);\r
-+ ngi.trignum = 0;\r
- \r
- for (j = 0; j < elementsonpoint[pi].Size(); j++)\r
- {\r
- sei = elementsonpoint[pi][j];\r
- const Element2d & bel = mesh[sei];\r
- surfi = mesh.GetFaceDescriptor(bel.GetIndex()).SurfNr();\r
-- \r
-+\r
- locelements.Append (sei);\r
- \r
- for (k = 1; k <= bel.GetNP(); k++)\r
- if (bel.PNum(k) == pi)\r
- {\r
- locrots.Append (k);\r
-+ ngi = bel.GeomInfoPi(k);\r
- break;\r
- }\r
- \r
-@@ -942,7 +936,7 @@\r
- }\r
- \r
- //optimizer loop (if not whole distance is not possible, move only a bit!!!!)\r
-- while (loci <= 5 && !moveisok)\r
-+ while (loci <= 5 && !moveisok)\r
- {\r
- loci ++;\r
- mesh[pi].X() = origp.X() + (x.Get(1) * t1.X() + x.Get(2) * t2.X())*fact;\r
-@@ -951,11 +945,9 @@\r
- fact = fact/2.;\r
- \r
- \r
-- ProjectPoint (surfi, mesh[pi]);\r
-+ moveisok = ProjectPoint (surfi, mesh[pi], ngi);\r
- \r
-- moveisok = CalcPointGeomInfo(surfi, ngi, mesh[pi]); \r
-- // point lies on same chart in stlsurface\r
-- \r
-+ // point lies on same chart in stlsurface\r
- if (moveisok)\r
- {\r
- for (j = 0; j < locelements.Size(); j++)\r
-diff -Nru netgen-4.5_orig/libsrc/occ/occconstruction.cpp netgen-4.5_patch/libsrc/occ/occconstruction.cpp\r
---- netgen-4.5_orig/libsrc/occ/occconstruction.cpp 2005-12-06 17:15:53.000000000 +0300\r
-+++ netgen-4.5_patch/libsrc/occ/occconstruction.cpp 2006-04-27 13:12:54.000000000 +0400\r
-@@ -28,8 +28,8 @@\r
- #include <BRepAlgoAPI_Common.hxx>\r
- #include <BRepAlgoAPI_Fuse.hxx>\r
- #include <BRepAlgoAPI_Section.hxx>\r
--#include <BRepOffsetAPI_Sewing.hxx>\r
--#include <BRepAlgo_Sewing.hxx>\r
-+//#include <BRepOffsetAPI_Sewing.hxx>\r
-+//#include <BRepAlgo_Sewing.hxx>\r
- #include <BRepOffsetAPI_MakeOffsetShape.hxx>\r
- #include <ShapeFix_Shape.hxx>\r
- namespace netgen\r
-diff -Nru netgen-4.5_orig/libsrc/occ/occgenmesh.cpp netgen-4.5_patch/libsrc/occ/occgenmesh.cpp\r
---- netgen-4.5_orig/libsrc/occ/occgenmesh.cpp 2006-02-07 12:12:48.000000000 +0300\r
-+++ netgen-4.5_patch/libsrc/occ/occgenmesh.cpp 2006-10-25 16:14:48.000000000 +0400\r
-@@ -28,7 +28,7 @@\r
- return Point<3> (p.X(), p.Y(), p.Z());\r
- }\r
- \r
-- void DivideEdge (TopoDS_Edge & edge,\r
-+ static void DivideEdge (TopoDS_Edge & edge,\r
- ARRAY<MeshPoint> & ps,\r
- ARRAY<double> & params,\r
- Mesh & mesh)\r
-@@ -49,23 +49,19 @@\r
- hvalue[0] = 0;\r
- pnt = c->Value(s0);\r
- \r
-- double olddist = 0;\r
-- double dist = 0;\r
--\r
-- for (int i = 1; i <= DIVIDEEDGESECTIONS; i++)\r
-+ int i;\r
-+ for (i = 1; i <= DIVIDEEDGESECTIONS; i++)\r
- {\r
- oldpnt = pnt;\r
- pnt = c->Value(s0+(i/double(DIVIDEEDGESECTIONS))*(s1-s0));\r
-+ double dist = pnt.Distance(oldpnt);\r
- hvalue[i] = hvalue[i-1] +\r
- 1.0/mesh.GetH(Point3d(pnt.X(), pnt.Y(), pnt.Z()))*\r
-- pnt.Distance(oldpnt);\r
-+ dist;\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
--\r
-- olddist = dist;\r
-- dist = pnt.Distance(oldpnt);\r
- }\r
- \r
- // nsubedges = int(ceil(hvalue[DIVIDEEDGESECTIONS]));\r
-@@ -74,7 +70,7 @@\r
- ps.SetSize(nsubedges-1);\r
- params.SetSize(nsubedges+1);\r
- \r
-- int i = 1;\r
-+ i = 1;\r
- int i1 = 0;\r
- do\r
- {\r
-@@ -112,7 +108,7 @@\r
- \r
- static void FindEdges (OCCGeometry & geom, Mesh & mesh)\r
- {\r
-- char * savetask = multithread.task;\r
-+ const char * savetask = multithread.task;\r
- multithread.task = "Edge meshing";\r
- \r
- (*testout) << "edge meshing" << endl;\r
-@@ -124,6 +120,7 @@\r
- (*testout) << "nedges = " << nedges << endl;\r
- \r
- double eps = 1e-6 * geom.GetBoundingBox().Diam();\r
-+ double eps2 = eps * eps;\r
- \r
- for (int i = 1; i <= nvertices; i++)\r
- {\r
-@@ -133,7 +130,7 @@\r
- bool exists = 0;\r
- if (merge_solids)\r
- for (PointIndex pi = 1; pi <= mesh.GetNP(); pi++)\r
-- if ( Dist2 (mesh[pi], Point<3>(mp)) < eps*eps)\r
-+ if ( Dist2 (mesh[pi], Point<3>(mp)) < eps2)\r
- {\r
- exists = 1;\r
- break;\r
-@@ -276,8 +273,8 @@\r
- pnums.Last() = -1;\r
- for (PointIndex pi = 1; pi < first_ep; pi++)\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) < eps2) pnums[0] = pi;\r
-+ if (Dist2 (mesh[pi], lp) < eps2) pnums.Last() = pi;\r
- }\r
- }\r
- \r
-@@ -287,7 +284,7 @@\r
- bool exists = 0;\r
- int j;\r
- for (j = first_ep; j <= mesh.GetNP(); j++)\r
-- if ((mesh.Point(j)-Point<3>(mp[i-1])).Length() < eps) \r
-+ if (Dist2(mesh.Point(j), Point<3>(mp[i-1])) < eps2)\r
- {\r
- exists = 1;\r
- break;\r
-@@ -394,7 +391,7 @@\r
- int i, j, k;\r
- int changed;\r
- \r
-- char * savetask = multithread.task;\r
-+ const char * savetask = multithread.task;\r
- multithread.task = "Surface meshing";\r
- \r
- geom.facemeshstatus = 0;\r
-@@ -751,7 +748,7 @@\r
- multithread.task = savetask;\r
- }\r
- \r
-- double ComputeH (double kappa)\r
-+ static double ComputeH (double kappa)\r
- {\r
- double hret;\r
- kappa *= mparam.curvaturesafety;\r
-@@ -779,7 +776,7 @@\r
- double nq = n*q;\r
- \r
- Point<3> p = p0 + 0.5*n;\r
-- double lambda = (p-l.p0)*n / nq;\r
-+ double lambda = (fabs(nq) > 1e-10 ? (p-l.p0)*n / nq : -1);\r
- \r
- if (lambda >= 0 && lambda <= 1)\r
- {\r
-@@ -799,55 +796,55 @@\r
- \r
- \r
- \r
-- void RestrictHTriangle (gp_Pnt2d & par0, gp_Pnt2d & par1, gp_Pnt2d & par2,\r
-- BRepLProp_SLProps * prop, Mesh & mesh, const double maxside, int depth, double h = 0)\r
-+ static void RestrictHTriangle (gp_Pnt2d & par0, gp_Pnt2d & par1, gp_Pnt2d & par2,\r
-+ BRepAdaptor_Surface& surf, Mesh & mesh, const double maxside, int depth, double h = 0)\r
- {\r
--\r
-+ BRepLProp_SLProps prop(surf, 2, 1e-5);\r
- \r
- gp_Pnt2d parmid;\r
- \r
- parmid.SetX(0.3*(par0.X()+par1.X()+par2.X()));\r
- parmid.SetY(0.3*(par0.Y()+par1.Y()+par2.Y()));\r
- \r
-- if (depth == 0)\r
-+ //if (depth == 0)\r
- {\r
- double curvature = 0;\r
- \r
-- prop->SetParameters (parmid.X(), parmid.Y());\r
-- if (!prop->IsCurvatureDefined())\r
-+ prop.SetParameters (parmid.X(), parmid.Y());\r
-+ if (!prop.IsCurvatureDefined())\r
- {\r
- (*testout) << "curvature not defined!" << endl;\r
- return;\r
- }\r
-- curvature = max(fabs(prop->MinCurvature()),\r
-- fabs(prop->MaxCurvature()));\r
-+ curvature = max(fabs(prop.MinCurvature()),\r
-+ fabs(prop.MaxCurvature()));\r
- \r
-- prop->SetParameters (par0.X(), par0.Y());\r
-- if (!prop->IsCurvatureDefined())\r
-+ prop.SetParameters (par0.X(), par0.Y());\r
-+ if (!prop.IsCurvatureDefined())\r
- {\r
- (*testout) << "curvature not defined!" << endl;\r
- return;\r
- }\r
-- curvature = max(curvature,max(fabs(prop->MinCurvature()),\r
-- fabs(prop->MaxCurvature())));\r
-+ curvature = max(curvature,max(fabs(prop.MinCurvature()),\r
-+ fabs(prop.MaxCurvature())));\r
- \r
-- prop->SetParameters (par1.X(), par1.Y());\r
-- if (!prop->IsCurvatureDefined())\r
-+ prop.SetParameters (par1.X(), par1.Y());\r
-+ if (!prop.IsCurvatureDefined())\r
- {\r
- (*testout) << "curvature not defined!" << endl;\r
- return;\r
- }\r
-- curvature = max(curvature,max(fabs(prop->MinCurvature()),\r
-- fabs(prop->MaxCurvature())));\r
-+ curvature = max(curvature,max(fabs(prop.MinCurvature()),\r
-+ fabs(prop.MaxCurvature())));\r
- \r
-- prop->SetParameters (par2.X(), par2.Y());\r
-- if (!prop->IsCurvatureDefined())\r
-+ prop.SetParameters (par2.X(), par2.Y());\r
-+ if (!prop.IsCurvatureDefined())\r
- {\r
- (*testout) << "curvature not defined!" << endl;\r
- return;\r
- }\r
-- curvature = max(curvature,max(fabs(prop->MinCurvature()),\r
-- fabs(prop->MaxCurvature())));\r
-+ curvature = max(curvature,max(fabs(prop.MinCurvature()),\r
-+ fabs(prop.MaxCurvature())));\r
- \r
- //(*testout) << "curvature " << curvature << endl;\r
- \r
-@@ -886,51 +883,47 @@\r
- pm1.SetX(0.5*(par0.X()+par2.X())); pm1.SetY(0.5*(par0.Y()+par2.Y()));\r
- pm2.SetX(0.5*(par1.X()+par0.X())); pm2.SetY(0.5*(par1.Y()+par0.Y()));\r
- \r
-- RestrictHTriangle (pm0, pm1, pm2, prop, mesh, 0.5*maxside, depth+1, h);\r
-- RestrictHTriangle (par0, pm1, pm2, prop, mesh, 0.5*maxside, depth+1, h);\r
-- RestrictHTriangle (par1, pm0, pm2, prop, mesh, 0.5*maxside, depth+1, h);\r
-- RestrictHTriangle (par2, pm1, pm0, prop, mesh, 0.5*maxside, depth+1, h);\r
-+ RestrictHTriangle (pm0, pm1, pm2, surf, mesh, 0.5*maxside, depth+1, h);\r
-+ RestrictHTriangle (par0, pm1, pm2, surf, mesh, 0.5*maxside, depth+1, h);\r
-+ RestrictHTriangle (par1, pm0, pm2, surf, mesh, 0.5*maxside, depth+1, h);\r
-+ RestrictHTriangle (par2, pm1, pm0, surf, mesh, 0.5*maxside, depth+1, h);\r
- }\r
- else\r
- {\r
- gp_Pnt pnt;\r
- Point3d p3d;\r
- \r
-- prop->SetParameters (parmid.X(), parmid.Y());\r
-- pnt = prop->Value();\r
-+ surf.D0(parmid.X(), parmid.Y(), pnt);\r
- p3d = Point3d(pnt.X(), pnt.Y(), pnt.Z());\r
- mesh.RestrictLocalH (p3d, h);\r
- \r
- \r
-- prop->SetParameters (par0.X(), par0.Y());\r
-- pnt = prop->Value();\r
-+ surf.D0(par0.X(), par0.Y(), pnt);\r
- p3d = Point3d(pnt.X(), pnt.Y(), pnt.Z());\r
- mesh.RestrictLocalH (p3d, h);\r
- \r
-- prop->SetParameters (par1.X(), par1.Y());\r
-- pnt = prop->Value();\r
-+ surf.D0(par1.X(), par1.Y(), pnt);\r
- p3d = Point3d(pnt.X(), pnt.Y(), pnt.Z());\r
- mesh.RestrictLocalH (p3d, h);\r
- \r
-- prop->SetParameters (par2.X(), par2.Y());\r
-- pnt = prop->Value();\r
-+ surf.D0(par2.X(), par2.Y(), pnt);\r
- p3d = Point3d(pnt.X(), pnt.Y(), pnt.Z());\r
- mesh.RestrictLocalH (p3d, h);\r
- \r
-- (*testout) << "p = " << p3d << ", h = " << h << ", maxside = " << maxside << endl;\r
-+ //(*testout) << "p = " << p3d << ", h = " << h << ", maxside = " << maxside << endl;\r
- /*\r
- (*testout) << pnt.X() << " " << pnt.Y() << " " << pnt.Z() << endl;\r
- \r
-- prop->SetParameters (par0.X(), par0.Y());\r
-- pnt = prop->Value();\r
-+ prop.SetParameters (par0.X(), par0.Y());\r
-+ pnt = prop.Value();\r
- (*testout) << pnt.X() << " " << pnt.Y() << " " << pnt.Z() << endl;\r
- \r
-- prop->SetParameters (par1.X(), par1.Y());\r
-- pnt = prop->Value();\r
-+ prop.SetParameters (par1.X(), par1.Y());\r
-+ pnt = prop.Value();\r
- (*testout) << pnt.X() << " " << pnt.Y() << " " << pnt.Z() << endl;\r
- \r
-- prop->SetParameters (par2.X(), par2.Y());\r
-- pnt = prop->Value();\r
-+ prop.SetParameters (par2.X(), par2.Y());\r
-+ pnt = prop.Value();\r
- (*testout) << pnt.X() << " " << pnt.Y() << " " << pnt.Z() << endl;\r
- */\r
- }\r
-@@ -970,7 +963,7 @@\r
- if (mparam.uselocalh)\r
- {\r
- \r
-- char * savetask = multithread.task;\r
-+ const char * savetask = multithread.task;\r
- multithread.percent = 0;\r
- \r
- mesh->SetLocalH (bb.PMin(), bb.PMax(), mparam.grading);\r
-@@ -1075,7 +1068,6 @@\r
- if (triangulation.IsNull()) continue;\r
- \r
- BRepAdaptor_Surface sf(face, Standard_True);\r
-- BRepLProp_SLProps prop(sf, 2, 1e-5);\r
- \r
- int ntriangles = triangulation -> NbTriangles();\r
- for (int j = 1; j <= ntriangles; j++)\r
-@@ -1096,7 +1088,7 @@\r
- maxside = max (maxside, p[1].Distance(p[2]));\r
- //cout << "\rFace " << i << " pos11 ntriangles " << ntriangles << " maxside " << maxside << flush;\r
- \r
-- RestrictHTriangle (par[0], par[1], par[2], &prop, *mesh, maxside, 0);\r
-+ RestrictHTriangle (par[0], par[1], par[2], sf, *mesh, maxside, 0);\r
- //cout << "\rFace " << i << " pos12 ntriangles " << ntriangles << flush;\r
- }\r
- }\r
-diff -Nru netgen-4.5_orig/libsrc/occ/occgeom.cpp netgen-4.5_patch/libsrc/occ/occgeom.cpp\r
---- netgen-4.5_orig/libsrc/occ/occgeom.cpp 2006-01-25 15:35:50.000000000 +0300\r
-+++ netgen-4.5_patch/libsrc/occ/occgeom.cpp 2006-10-25 16:15:24.000000000 +0400\r
-@@ -7,6 +7,8 @@\r
- #include "ShapeAnalysis_ShapeContents.hxx"\r
- #include "ShapeAnalysis_CheckSmallFace.hxx"\r
- #include "ShapeAnalysis_DataMapOfShapeListOfReal.hxx"\r
-+#include <ShapeAnalysis_Surface.hxx>\r
-+#include <BRepTopAdaptor_FClass2d.hxx>\r
- #include "BRepAlgoAPI_Fuse.hxx"\r
- #include "BRepCheck_Analyzer.hxx"\r
- #include "BRepLib.hxx"\r
-@@ -16,11 +18,19 @@\r
- #include "Partition_Spliter.hxx"\r
- //#include "VrmlAPI.hxx"\r
- //#include "StlAPI.hxx"\r
-+#include <TopAbs_State.hxx>\r
- \r
- \r
- namespace netgen\r
- {\r
- \r
-+ OCCGeometry::~OCCGeometry()\r
-+ {\r
-+ NCollection_DataMap<int,BRepTopAdaptor_FClass2d*>::Iterator it(fclsmap);\r
-+ for (; it.More(); it.Next())\r
-+ delete it.Value();\r
-+ }\r
-+\r
- void OCCGeometry :: PrintNrShapes ()\r
- {\r
- TopExp_Explorer e;\r
-@@ -947,13 +957,13 @@\r
- \r
- void OCCGeometry :: BuildVisualizationMesh ()\r
- {\r
--\r
-- cout << "Preparing visualization (deflection = " << vispar.occdeflection << ") ... " << flush;\r
-+ double vispar_occdeflection = 0.01;\r
-+ cout << "Preparing visualization (deflection = " << vispar_occdeflection << ") ... " << flush;\r
- \r
- \r
- BRepTools::Clean (shape);\r
- //WriteOCC_STL("test.stl");\r
-- BRepMesh_IncrementalMesh::BRepMesh_IncrementalMesh (shape, vispar.occdeflection, true);\r
-+ BRepMesh_IncrementalMesh::BRepMesh_IncrementalMesh (shape, vispar_occdeflection, true);\r
- cout << "done" << endl;\r
- \r
- \r
-@@ -973,8 +983,27 @@\r
- \r
- }\r
- \r
-+ void OCCGeometry::GetFaceTools(int surfi, Handle(ShapeAnalysis_Surface)& proj,\r
-+ BRepTopAdaptor_FClass2d*& cls) const\r
-+ {\r
-+ //MSV: organize caching projector in the map\r
-+ if (fprjmap.IsBound(surfi))\r
-+ {\r
-+ proj = fprjmap.Find(surfi);\r
-+ cls = fclsmap.Find(surfi);\r
-+ }\r
-+ else\r
-+ {\r
-+ const TopoDS_Face& aFace = TopoDS::Face(fmap(surfi));\r
-+ Handle(Geom_Surface) aSurf = BRep_Tool::Surface(aFace);\r
-+ proj = new ShapeAnalysis_Surface(aSurf);\r
-+ fprjmap.Bind(surfi, proj);\r
-+ cls = new BRepTopAdaptor_FClass2d(aFace,Precision::Confusion());\r
-+ fclsmap.Bind(surfi, cls);\r
-+ }\r
-+ }\r
- \r
-- void OCCGeometry :: Project (int surfi, Point<3> & p) const\r
-+ bool OCCGeometry :: Project (int surfi, Point<3> & p, double& u, double& v) const\r
- {\r
- static int cnt = 0;\r
- if (++cnt % 1000 == 0) cout << "Project cnt = " << cnt << endl;\r
-@@ -983,18 +1012,22 @@\r
- \r
- //(*testout) << "before " << pnt.X() << " "<< pnt.Y() << " "<< pnt.Z() << " " << endl;\r
- \r
-- GeomAPI_ProjectPointOnSurf proj(pnt, BRep_Tool::Surface(TopoDS::Face(fmap(surfi))));\r
-- if (proj.NbPoints() == 0)\r
-- {\r
-- cout << "Projection fails" << endl;\r
-- }\r
-- else\r
-- {\r
-- pnt = proj.NearestPoint();\r
-- //(*testout) << "after " << pnt.X() << " "<< pnt.Y() << " "<< pnt.Z() << " " << endl;\r
-+ Handle(ShapeAnalysis_Surface) proj;\r
-+ BRepTopAdaptor_FClass2d *cls;\r
-+ GetFaceTools(surfi, proj, cls);\r
- \r
-- p = Point<3> (pnt.X(), pnt.Y(), pnt.Z());\r
-- }\r
-+ gp_Pnt2d p2d = proj->ValueOfUV(pnt, Precision::Confusion());\r
-+ if (cls->Perform(p2d) == TopAbs_OUT)\r
-+ {\r
-+ //cout << "Projection fails" << endl;\r
-+ return false;\r
-+ }\r
-+ pnt = proj->Value(p2d);\r
-+ p2d.Coord(u, v);\r
-+ //(*testout) << "after " << pnt.X() << " "<< pnt.Y() << " "<< pnt.Z() << " " << endl;\r
-+\r
-+ p = Point<3> (pnt.X(), pnt.Y(), pnt.Z());\r
-+ return true;\r
- }\r
- \r
- \r
-@@ -1002,54 +1035,20 @@\r
- {\r
- gp_Pnt p(ap(0), ap(1), ap(2));\r
- \r
-- Handle(Geom_Surface) surface = BRep_Tool::Surface(TopoDS::Face(fmap(surfi)));\r
-+ Handle(ShapeAnalysis_Surface) proj;\r
-+ BRepTopAdaptor_FClass2d *cls;\r
-+ GetFaceTools(surfi, proj, cls);\r
- \r
-- gp_Pnt x = surface->Value (u,v);\r
-- \r
-- if (p.SquareDistance(x) <= sqr(PROJECTION_TOLERANCE)) return true;\r
-- \r
-- gp_Vec du, dv;\r
-- \r
-- surface->D1(u,v,x,du,dv);\r
-- \r
-- int count = 0;\r
-- \r
-- gp_Pnt xold;\r
-- gp_Vec n;\r
-- double det, lambda, mu;\r
-- \r
-- do {\r
-- count++;\r
-- \r
-- n = du^dv;\r
-- \r
-- det = Det3 (n.X(), du.X(), dv.X(),\r
-- n.Y(), du.Y(), dv.Y(),\r
-- n.Z(), du.Z(), dv.Z());\r
-- \r
-- if (det < 1e-15) return false; \r
-- \r
-- lambda = Det3 (n.X(), p.X()-x.X(), dv.X(),\r
-- n.Y(), p.Y()-x.Y(), dv.Y(),\r
-- n.Z(), p.Z()-x.Z(), dv.Z())/det;\r
-- \r
-- mu = Det3 (n.X(), du.X(), p.X()-x.X(),\r
-- n.Y(), du.Y(), p.Y()-x.Y(),\r
-- n.Z(), du.Z(), p.Z()-x.Z())/det;\r
-- \r
-- u += lambda;\r
-- v += mu;\r
-- \r
-- xold = x;\r
-- surface->D1(u,v,x,du,dv);\r
-- \r
-- } while (xold.SquareDistance(x) > sqr(PROJECTION_TOLERANCE) && count < 50);\r
--\r
-- // (*testout) << "FastProject count: " << count << endl;\r
-- \r
-- if (count == 50) return false;\r
-+ gp_Pnt2d p2d = proj->NextValueOfUV(gp_Pnt2d(u,v), p, Precision::Confusion());\r
-+ if (cls->Perform(p2d) == TopAbs_OUT)\r
-+ {\r
-+ //cout << "Projection fails" << endl;\r
-+ return false;\r
-+ }\r
- \r
-- ap = Point<3> (x.X(), x.Y(), x.Z());\r
-+ p = proj->Value(p2d);\r
-+ p2d.Coord(u, v);\r
-+ ap = Point<3> (p.X(), p.Y(), p.Z());\r
- \r
- return true;\r
- }\r
-diff -Nru netgen-4.5_orig/libsrc/occ/occgeom.hpp netgen-4.5_patch/libsrc/occ/occgeom.hpp\r
---- netgen-4.5_orig/libsrc/occ/occgeom.hpp 2006-01-25 15:35:50.000000000 +0300\r
-+++ netgen-4.5_patch/libsrc/occ/occgeom.hpp 2006-10-25 16:16:01.000000000 +0400\r
-@@ -15,8 +15,6 @@\r
- #include "Geom_Curve.hxx"\r
- #include "Geom2d_Curve.hxx"\r
- #include "Geom_Surface.hxx"\r
--#include "GeomAPI_ProjectPointOnSurf.hxx"\r
--#include "GeomAPI_ProjectPointOnCurve.hxx"\r
- #include "BRepTools.hxx"\r
- #include "TopExp.hxx"\r
- #include "BRepBuilderAPI_MakeVertex.hxx"\r
-@@ -41,8 +39,6 @@\r
- #include "Geom_Curve.hxx"\r
- #include "Geom2d_Curve.hxx"\r
- #include "Geom_Surface.hxx"\r
--#include "GeomAPI_ProjectPointOnSurf.hxx"\r
--#include "GeomAPI_ProjectPointOnCurve.hxx"\r
- #include "TopoDS_Wire.hxx"\r
- #include "BRepTools_WireExplorer.hxx"\r
- #include "BRepTools.hxx"\r
-@@ -69,7 +65,7 @@\r
- #include "IGESToBRep_Reader.hxx"\r
- #include "Interface_Static.hxx"\r
- #include "GeomAPI_ExtremaCurveCurve.hxx"\r
--#include "Standard_ErrorHandler.hxx"\r
-+//#include "Standard_ErrorHandler.hxx"\r
- #include "Standard_Failure.hxx"\r
- #include "ShapeUpgrade_ShellSewing.hxx"\r
- #include "ShapeFix_Shape.hxx"\r
-@@ -84,11 +80,15 @@\r
- #include "STEPControl_Writer.hxx"\r
- #include "StlAPI_Writer.hxx"\r
- #include "STEPControl_StepModelType.hxx"\r
-+#include <NCollection_DataMap.hxx>\r
-+\r
-+class Handle_ShapeAnalysis_Surface;\r
-+class BRepTopAdaptor_FClass2d;\r
- \r
- namespace netgen\r
- {\r
- \r
--#include "../visualization/vispar.hpp"\r
-+ //#include "../visualization/vispar.hpp"\r
- // class VisualizationParameters;\r
- // extern VisualizationParameters vispar;\r
- \r
-@@ -159,6 +159,8 @@\r
- class OCCGeometry\r
- {\r
- Point<3> center;\r
-+ mutable NCollection_DataMap<int,Handle_ShapeAnalysis_Surface> fprjmap;\r
-+ mutable NCollection_DataMap<int,BRepTopAdaptor_FClass2d*> fclsmap;\r
- \r
- public:\r
- TopoDS_Shape shape;\r
-@@ -189,6 +191,7 @@\r
- vmap.Clear();\r
- }\r
- \r
-+ ~OCCGeometry();\r
- \r
- void BuildFMap();\r
- \r
-@@ -204,10 +207,12 @@\r
- Point<3> Center()\r
- { return center; }\r
- \r
-- void Project (int surfi, Point<3> & p) const;\r
-+ bool Project (int surfi, Point<3> & p, double& u, double& v) const;\r
- bool FastProject (int surfi, Point<3> & ap, double& u, double& v) const;\r
- \r
-- \r
-+ void GetFaceTools(int surfi, Handle(ShapeAnalysis_Surface)& proj,\r
-+ BRepTopAdaptor_FClass2d*& cls) const;\r
-+\r
- OCCSurface GetSurface (int surfi)\r
- {\r
- cout << "OCCGeometry::GetSurface using PLANESPACE" << endl;\r
-diff -Nru netgen-4.5_orig/libsrc/occ/occmeshsurf.cpp netgen-4.5_patch/libsrc/occ/occmeshsurf.cpp\r
---- netgen-4.5_orig/libsrc/occ/occmeshsurf.cpp 2006-01-25 15:36:26.000000000 +0300\r
-+++ netgen-4.5_patch/libsrc/occ/occmeshsurf.cpp 2006-10-25 16:16:26.000000000 +0400\r
-@@ -5,6 +5,8 @@\r
- #include <occgeom.hpp>\r
- #include <meshing.hpp>\r
- #include <GeomLProp_SLProps.hxx>\r
-+#include <GeomAPI_ProjectPointOnSurf.hxx>\r
-+#include <GeomAPI_ProjectPointOnCurve.hxx>\r
- \r
- \r
- namespace netgen\r
-@@ -411,11 +413,16 @@\r
- }\r
- \r
- \r
-- void MeshOptimize2dOCCSurfaces :: ProjectPoint (INDEX surfind, Point3d & p) const\r
-+ bool MeshOptimize2dOCCSurfaces :: ProjectPoint (INDEX surfind, Point3d & p, PointGeomInfo& gi) const\r
- {\r
- Point<3> hp = p;\r
-- geometry.Project (surfind, hp);\r
-+ bool ok;\r
-+ if (gi.trignum > 0)\r
-+ ok = geometry.FastProject (surfind, hp, gi.u, gi.v);\r
-+ else\r
-+ ok = geometry.Project (surfind, hp, gi.u, gi.v);\r
- p = hp;\r
-+ return ok;\r
- }\r
- \r
- void MeshOptimize2dOCCSurfaces :: ProjectPoint2 (INDEX surfind, INDEX surfind2, \r
-@@ -506,38 +513,6 @@\r
- }\r
- \r
- \r
-- int MeshOptimize2dOCCSurfaces :: \r
-- CalcPointGeomInfo(int surfind, PointGeomInfo& gi, const Point3d& p) const\r
-- {\r
-- Standard_Real u,v;\r
--\r
-- gp_Pnt pnt(p.X(), p.Y(), p.Z());\r
--\r
-- Handle(Geom_Surface) occface;\r
-- occface = BRep_Tool::Surface(TopoDS::Face(geometry.fmap(surfind)));\r
--\r
-- GeomAPI_ProjectPointOnSurf proj(pnt, occface);\r
--\r
-- if (proj.NbPoints() < 1)\r
-- {\r
-- cout << "ERROR: OCCSurface :: GetNormalVector: GeomAPI_ProjectPointOnSurf failed!"\r
-- << endl;\r
-- cout << p << endl;\r
-- return 0;\r
-- }\r
-- \r
-- proj.LowerDistanceParameters (u, v); \r
--\r
-- gi.u = u;\r
-- gi.v = v;\r
-- return 1;\r
-- }\r
--\r
--\r
--\r
--\r
--\r
--\r
- OCCRefinementSurfaces :: OCCRefinementSurfaces (const OCCGeometry & ageometry)\r
- : Refinement(), geometry(ageometry)\r
- {\r
-@@ -627,10 +602,11 @@\r
- if (!geometry.FastProject (surfi, hnewp, u, v))\r
- {\r
- cout << "Fast projection to surface fails! Using OCC projection" << endl;\r
-- geometry.Project (surfi, hnewp);\r
-+ double u, v;\r
-+ geometry.Project (surfi, hnewp, u, v);\r
- }\r
- \r
-- newgi.trignum = 1;\r
-+ newgi.trignum = surfi;\r
- }\r
- \r
- newp = hnewp;\r
-@@ -653,14 +629,17 @@\r
- hnewp = Point<3> (pnt.X(), pnt.Y(), pnt.Z());\r
- newp = hnewp;\r
- newgi = ap1;\r
-- };\r
-+ }\r
- \r
- \r
- void OCCRefinementSurfaces :: ProjectToSurface (Point<3> & p, int surfi)\r
- {\r
- if (surfi > 0)\r
-- geometry.Project (surfi, p);\r
-- };\r
-+ {\r
-+ double u, v;\r
-+ geometry.Project (surfi, p, u, v);\r
-+ }\r
-+ }\r
- \r
- void OCCRefinementSurfaces :: ProjectToSurface (Point<3> & p, int surfi, PointGeomInfo & gi)\r
- {\r
-@@ -668,9 +647,10 @@\r
- if (!geometry.FastProject (surfi, p, gi.u, gi.v))\r
- {\r
- cout << "Fast projection to surface fails! Using OCC projection" << endl;\r
-- geometry.Project (surfi, p);\r
-+ double u, v;\r
-+ geometry.Project (surfi, p, u, v);\r
- }\r
-- };\r
-+ }\r
- \r
- \r
- \r
-diff -Nru netgen-4.5_orig/libsrc/occ/occmeshsurf.hpp netgen-4.5_patch/libsrc/occ/occmeshsurf.hpp\r
---- netgen-4.5_orig/libsrc/occ/occmeshsurf.hpp 2005-06-09 18:51:10.000000000 +0400\r
-+++ netgen-4.5_patch/libsrc/occ/occmeshsurf.hpp 2006-10-25 16:17:22.000000000 +0400\r
-@@ -151,7 +151,7 @@\r
- MeshOptimize2dOCCSurfaces (const OCCGeometry & ageometry); \r
- \r
- ///\r
-- virtual void ProjectPoint (INDEX surfind, Point3d & p) const;\r
-+ virtual bool ProjectPoint (INDEX surfind, Point3d & p, PointGeomInfo& gi) const;\r
- ///\r
- virtual void ProjectPoint2 (INDEX surfind, INDEX surfind2, Point3d & p) const;\r
- ///\r
-@@ -159,9 +159,6 @@\r
- ///\r
- virtual void GetNormalVector(INDEX surfind, const Point3d & p, PointGeomInfo & gi, Vec3d & n) const;\r
- \r
-- \r
-- virtual int CalcPointGeomInfo(int surfind, PointGeomInfo& gi, const Point3d& p3) const;\r
--\r
- };\r
- \r
- \r
-diff -Nru netgen-4.5_orig/libsrc/stlgeom/meshstlsurface.cpp netgen-4.5_patch/libsrc/stlgeom/meshstlsurface.cpp\r
---- netgen-4.5_orig/libsrc/stlgeom/meshstlsurface.cpp 2006-01-11 18:08:20.000000000 +0300\r
-+++ netgen-4.5_patch/libsrc/stlgeom/meshstlsurface.cpp 2006-10-25 16:17:47.000000000 +0400\r
-@@ -946,20 +946,23 @@\r
- }\r
- \r
- \r
--void MeshOptimizeSTLSurface :: ProjectPoint (INDEX surfind, Point3d & p) const\r
-+bool MeshOptimizeSTLSurface :: ProjectPoint (INDEX surfind, Point3d & p, PointGeomInfo& gi) const\r
- {\r
- Point<3> hp = p;\r
-- if (!geom.Project (hp))\r
-+ if (gi.trignum > 0)\r
-+ ((STLGeometry&)geom).SelectChartOfTriangle (gi.trignum);\r
-+ if (!(gi.trignum = geom.Project (hp)))\r
- {\r
- PrintMessage(7,"project failed");\r
- \r
-- if (!geom.ProjectOnWholeSurface(hp)) \r
-+ if (!(gi.trignum = geom.ProjectOnWholeSurface(hp))) \r
- {\r
- PrintMessage(7, "project on whole surface failed");\r
- }\r
- }\r
- p = hp;\r
- // geometry.GetSurface(surfind)->Project (p);\r
-+ return gi.trignum > 0;\r
- }\r
- \r
- void MeshOptimizeSTLSurface :: ProjectPoint2 (INDEX surfind, INDEX surfind2, Point3d & p) const\r
-@@ -970,20 +973,6 @@\r
- */\r
- }\r
- \r
--int MeshOptimizeSTLSurface :: CalcPointGeomInfo(PointGeomInfo& gi, const Point3d& p3) const\r
--{\r
-- Point<3> hp = p3;\r
-- gi.trignum = geom.Project (hp);\r
--\r
-- if (gi.trignum)\r
-- {\r
-- return 1;\r
-- }\r
--\r
-- return 0;\r
-- \r
--}\r
--\r
- void MeshOptimizeSTLSurface :: GetNormalVector(INDEX surfind, const Point3d & p, Vec3d & n) const\r
- {\r
- n = geom.GetChartNormalVector();\r
-diff -Nru netgen-4.5_orig/libsrc/stlgeom/meshstlsurface.hpp netgen-4.5_patch/libsrc/stlgeom/meshstlsurface.hpp\r
---- netgen-4.5_orig/libsrc/stlgeom/meshstlsurface.hpp 2004-09-30 17:13:56.000000000 +0400\r
-+++ netgen-4.5_patch/libsrc/stlgeom/meshstlsurface.hpp 2006-10-25 16:17:59.000000000 +0400\r
-@@ -79,12 +79,10 @@\r
- virtual void SelectSurfaceOfPoint (const Point3d & p,\r
- const PointGeomInfo & gi);\r
- ///\r
-- virtual void ProjectPoint (INDEX surfind, Point3d & p) const;\r
-+ virtual bool ProjectPoint (INDEX surfind, Point3d & p, PointGeomInfo& gi) const;\r
- ///\r
- virtual void ProjectPoint2 (INDEX surfind, INDEX surfind2, Point3d & p) const;\r
- ///\r
-- virtual int CalcPointGeomInfo(PointGeomInfo& gi, const Point3d& p3) const;\r
-- ///\r
- virtual void GetNormalVector(INDEX surfind, const Point3d & p, Vec3d & n) const;\r
- };\r
- \r
-diff -Nru netgen-4.5_orig/makeForSalome.sh netgen-4.5_patch/makeForSalome.sh\r
---- netgen-4.5_orig/makeForSalome.sh 1970-01-01 03:00:00.000000000 +0300\r
-+++ netgen-4.5_patch/makeForSalome.sh 2006-04-27 13:12:54.000000000 +0400\r
-@@ -0,0 +1,31 @@\r
-+#! /bin/sh\r
-+cp ngtcltk/ngnewdelete.* libsrc/interface/\r
-+\r
-+MACHINE=LINUX\r
-+export MACHINE\r
-+make -C libsrc/csg\r
-+make -C libsrc/general\r
-+make -C libsrc/geom2d\r
-+make -C libsrc/gprim\r
-+make -C libsrc/interface\r
-+make -C libsrc/linalg\r
-+make -C libsrc/meshing\r
-+make -C libsrc/opti\r
-+make -C libsrc/stlgeom\r
-+make -C libsrc/occ\r
-+\r
-+if [ ! -d install ] ; then\r
-+ mkdir install\r
-+fi\r
-+\r
-+cp -r lib install/\r
-+\r
-+if [ ! -d install/include ] ; then\r
-+ mkdir install/include\r
-+fi\r
-+\r
-+cp libsrc/interface/nglib.h libsrc/general/*.hpp libsrc/csg/*.hpp libsrc/geom2d/*.hpp \\r
-+ libsrc/gprim/*.hpp libsrc/linalg/*.hpp libsrc/meshing/*.hpp \\r
-+ libsrc/occ/*.hpp libsrc/opti/*.hpp libsrc/include/mydefs.hpp \\r
-+ libsrc/stlgeom/*.hpp libsrc/include/mystdlib.h \\r
-+ install/include\r
+diff -Nru netgen-4.5_orig/libsrc/csg/meshsurf.cpp netgen-4.5_patch/libsrc/csg/meshsurf.cpp
+--- netgen-4.5_orig/libsrc/csg/meshsurf.cpp 2006-02-14 10:54:35.000000000 +0300
++++ netgen-4.5_patch/libsrc/csg/meshsurf.cpp 2006-10-25 16:05:59.000000000 +0400
+@@ -77,11 +77,12 @@
+ }
+
+
+-void MeshOptimize2dSurfaces :: ProjectPoint (INDEX surfind, Point3d & p) const
++bool MeshOptimize2dSurfaces :: ProjectPoint (INDEX surfind, Point3d & p) const
+ {
+ Point<3> hp = p;
+ geometry.GetSurface(surfind)->Project (hp);
+ p = hp;
++ return true;
+ }
+
+ void MeshOptimize2dSurfaces :: ProjectPoint2 (INDEX surfind, INDEX surfind2,
+diff -Nru netgen-4.5_orig/libsrc/csg/meshsurf.hpp netgen-4.5_patch/libsrc/csg/meshsurf.hpp
+--- netgen-4.5_orig/libsrc/csg/meshsurf.hpp 2004-01-20 13:49:44.000000000 +0300
++++ netgen-4.5_patch/libsrc/csg/meshsurf.hpp 2006-10-25 16:08:05.000000000 +0400
+@@ -45,7 +45,7 @@
+ MeshOptimize2dSurfaces (const CSGeometry & ageometry);
+
+ ///
+- virtual void ProjectPoint (INDEX surfind, Point3d & p) const;
++ virtual bool ProjectPoint (INDEX surfind, Point3d & p) const;
+ ///
+ virtual void ProjectPoint2 (INDEX surfind, INDEX surfind2, Point3d & p) const;
+ ///
+diff -Nru netgen-4.5_orig/libsrc/interface/Makefile netgen-4.5_patch/libsrc/interface/Makefile
+--- netgen-4.5_orig/libsrc/interface/Makefile 2005-08-09 18:14:59.000000000 +0400
++++ netgen-4.5_patch/libsrc/interface/Makefile 2006-04-27 13:12:54.000000000 +0400
+@@ -1,4 +1,5 @@
+-src = nginterface.cpp writeuser.cpp writediffpack.cpp writeabaqus.cpp writefluent.cpp writepermas.cpp writetochnog.cpp writetecplot.cpp wuchemnitz.cpp writetochnog.cpp writefeap.cpp writeelmer.cpp writegmsh.cpp writejcm.cpp readuser.cpp importsolution.cpp
++#src = nginterface.cpp writeuser.cpp writediffpack.cpp writeabaqus.cpp writefluent.cpp writepermas.cpp writetochnog.cpp writetecplot.cpp wuchemnitz.cpp writetochnog.cpp writefeap.cpp writeelmer.cpp writegmsh.cpp writejcm.cpp readuser.cpp importsolution.cpp
++src = writeuser.cpp writediffpack.cpp writeabaqus.cpp writefluent.cpp writepermas.cpp writetochnog.cpp writetecplot.cpp wuchemnitz.cpp writetochnog.cpp writefeap.cpp writeelmer.cpp writegmsh.cpp writejcm.cpp readuser.cpp nglib.cpp ngnewdelete.cpp
+ #
+ lib = nginterface
+ libpath = libsrc/interface
+diff -Nru netgen-4.5_orig/libsrc/interface/nglib.cpp netgen-4.5_patch/libsrc/interface/nglib.cpp
+--- netgen-4.5_orig/libsrc/interface/nglib.cpp 2005-10-18 17:53:18.000000000 +0400
++++ netgen-4.5_patch/libsrc/interface/nglib.cpp 2006-04-27 13:12:54.000000000 +0400
+@@ -56,7 +56,8 @@
+
+ void Ng_Exit ()
+ {
+- ;
++ delete testout;
++ testout = NULL;
+ }
+
+
+diff -Nru netgen-4.5_orig/libsrc/makefile.inc netgen-4.5_patch/libsrc/makefile.inc
+--- netgen-4.5_orig/libsrc/makefile.inc 2005-09-02 17:17:51.000000000 +0400
++++ netgen-4.5_patch/libsrc/makefile.inc 2006-10-26 09:33:58.000000000 +0400
+@@ -8,17 +8,14 @@
+ LIBSRC_DIR=$(CPP_DIR)/libsrc
+ LIB_DIR=$(CPP_DIR)/lib/$(MACHINE)
+
+-#OCC_DIR=../../occ
+-#OCCINC_DIR=$(OCC_DIR)/inc
+-#OCCLIB_DIR=$(OCC_DIR)/lib
+-# OCC_DIR=/opt/OpenCASCADE5.2/ros
+-# OCC_DIR=/home/joachim/download/occ/Linux
+-# OCCINC_DIR=$(OCC_DIR)/inc -I$(OCC_DIR)/ros/inc
+-# OCCLIB_DIR=$(OCC_DIR)/Linux/lib
++OCC_DIR=$(CASROOT)
++OCCINC_DIR=$(OCC_DIR)/inc
++OCCLIB_DIR=$(OCC_DIR)/Linux/lib
+ #
+ include $(LIBSRC_DIR)/makefile.mach.$(MACHINE)
+ #
+-CPLUSPLUSFLAGS1 = -c -I$(LIBSRC_DIR)/include -I$(OCCINC_DIR)
++CPLUSPLUSFLAGS1 = -c -I$(LIBSRC_DIR)/include -I$(OCCINC_DIR) \
++ -DOCCGEOMETRY -DOCC52 -DHAVE_IOSTREAM -DHAVE_LIMITS_H
+ #
+ ARFLAGS = r
+ #
+diff -Nru netgen-4.5_orig/libsrc/makefile.mach.LINUX netgen-4.5_patch/libsrc/makefile.mach.LINUX
+--- netgen-4.5_orig/libsrc/makefile.mach.LINUX 2004-10-11 23:49:26.000000000 +0400
++++ netgen-4.5_patch/libsrc/makefile.mach.LINUX 2006-04-27 13:12:54.000000000 +0400
+@@ -16,7 +16,7 @@
+ #
+ CFLAGS2 =
+
+-CPLUSPLUSFLAGS2 = -O2 -I/usr/include/GL3.5 -DLINUX -DOPENGL \
++CPLUSPLUSFLAGS2 = -O2 -I/usr/include/GL3.5 -DLINUX \
+ -ftemplate-depth-99 -finline-limit=10000 \
+ -Wdisabled-optimization -funroll-loops -DnoNGSOLVE
+
+diff -Nru netgen-4.5_orig/libsrc/meshing/improve2.cpp netgen-4.5_patch/libsrc/meshing/improve2.cpp
+--- netgen-4.5_orig/libsrc/meshing/improve2.cpp 2006-01-11 18:08:19.000000000 +0300
++++ netgen-4.5_patch/libsrc/meshing/improve2.cpp 2006-04-27 13:12:54.000000000 +0400
+@@ -4,7 +4,7 @@
+ #include <opti.hpp>
+
+ #ifndef SMALLLIB
+-#include <visual.hpp>
++//#include <visual.hpp>
+ #endif
+
+ namespace netgen
+diff -Nru netgen-4.5_orig/libsrc/meshing/improve2.hpp netgen-4.5_patch/libsrc/meshing/improve2.hpp
+--- netgen-4.5_orig/libsrc/meshing/improve2.hpp 2004-10-12 23:22:55.000000000 +0400
++++ netgen-4.5_patch/libsrc/meshing/improve2.hpp 2006-10-25 16:09:37.000000000 +0400
+@@ -32,17 +32,16 @@
+ ///
+ virtual void SelectSurfaceOfPoint (const Point3d & p,
+ const PointGeomInfo & gi);
+- ///
+- virtual void ProjectPoint (INDEX /* surfind */, Point3d & /* p */) const { };
++
++ /// project point on surface, returns true if success
++ virtual bool ProjectPoint (INDEX /* surfind */, Point3d & /* p */) const { return false; }
++ /// fast project point on surface using point geom info of a neighboring point
++ /// if gi.trignum != 0,
++ /// returns true if success, gi is updated
++ virtual bool ProjectPoint (INDEX surfind, Point3d & p, PointGeomInfo& gi) const
++ { gi.trignum = 1; return ProjectPoint (surfind, p); }
+ ///
+ virtual void ProjectPoint2 (INDEX /* surfind */, INDEX /* surfind2 */, Point3d & /* p */) const { };
+- /// liefert zu einem 3d-Punkt die geominfo (Dreieck) und liefert 1, wenn erfolgreich,
+- /// 0, wenn nicht (Punkt ausserhalb von chart)
+- virtual int CalcPointGeomInfo(PointGeomInfo& gi, const Point3d& /*p3*/) const
+- { gi.trignum = 1; return 1;};
+-
+- virtual int CalcPointGeomInfo(int /* surfind */, PointGeomInfo& gi, const Point3d& p3) const
+- { return CalcPointGeomInfo (gi, p3); }
+
+ ///
+ virtual void GetNormalVector(INDEX surfind, const Point3d & p, PointGeomInfo & gi, Vec3d & n) const;
+diff -Nru netgen-4.5_orig/libsrc/meshing/smoothing2.cpp netgen-4.5_patch/libsrc/meshing/smoothing2.cpp
+--- netgen-4.5_orig/libsrc/meshing/smoothing2.cpp 2006-01-11 18:08:20.000000000 +0300
++++ netgen-4.5_patch/libsrc/meshing/smoothing2.cpp 2006-10-25 16:10:46.000000000 +0400
+@@ -300,7 +300,7 @@
+ double Opti2SurfaceMinFunction ::
+ FuncGrad (const Vector & x, Vector & grad) const
+ {
+- Vec3d n, vgrad;
++ Vec3d vgrad;
+ Point3d pp1;
+ double g1x, g1y;
+ double badness, hbadness;
+@@ -308,8 +308,6 @@
+ vgrad = 0;
+ badness = 0;
+
+- meshthis -> GetNormalVector (surfi, sp1, gi1, n);
+-
+ pp1 = sp1;
+ pp1.Add2 (x.Get(1), t1, x.Get(2), t2);
+
+@@ -360,7 +358,7 @@
+ double Opti2SurfaceMinFunction ::
+ FuncDeriv (const Vector & x, const Vector & dir, double & deriv) const
+ {
+- Vec3d n, vgrad;
++ Vec3d vgrad;
+ Point3d pp1;
+ double g1x, g1y;
+ double badness, hbadness;
+@@ -368,8 +366,6 @@
+ vgrad = 0;
+ badness = 0;
+
+- meshthis -> GetNormalVector (surfi, sp1, gi1, n);
+-
+ pp1 = sp1;
+ pp1.Add2 (x.Get(1), t1, x.Get(2), t2);
+
+@@ -520,7 +516,7 @@
+ // from 2d:
+
+ int j, k, lpi, gpi;
+- Vec3d n, vgrad;
++ Vec3d vgrad;
+ Point3d pp1;
+ Vec2d g1, vdir;
+ double badness, hbadness, hbad, hderiv;
+@@ -528,8 +524,6 @@
+ vgrad = 0;
+ badness = 0;
+
+- meshthis -> GetNormalVector (surfi, sp1, gi1, n);
+-
+ pp1 = sp1;
+ pp1.Add2 (x.Get(1), t1, x.Get(2), t2);
+
+@@ -593,7 +587,7 @@
+ // from 2d:
+
+ int j, k, lpi, gpi;
+- Vec3d n, vgrad;
++ Vec3d vgrad;
+ Point3d pp1;
+ Vec2d g1, vdir;
+ double badness, hbadness, hbad, hderiv;
+@@ -601,8 +595,6 @@
+ vgrad = 0;
+ badness = 0;
+
+- meshthis -> GetNormalVector (surfi, sp1, gi1, n);
+-
+ pp1 = sp1;
+ pp1.Add2 (x.Get(1), t1, x.Get(2), t2);
+
+@@ -859,19 +851,21 @@
+ locelements.SetSize(0);
+ locrots.SetSize (0);
+ lochs.SetSize (0);
++ ngi.trignum = 0;
+
+ for (j = 0; j < elementsonpoint[pi].Size(); j++)
+ {
+ sei = elementsonpoint[pi][j];
+ const Element2d & bel = mesh[sei];
+ surfi = mesh.GetFaceDescriptor(bel.GetIndex()).SurfNr();
+-
++
+ locelements.Append (sei);
+
+ for (k = 1; k <= bel.GetNP(); k++)
+ if (bel.PNum(k) == pi)
+ {
+ locrots.Append (k);
++ ngi = bel.GeomInfoPi(k);
+ break;
+ }
+
+@@ -942,7 +936,7 @@
+ }
+
+ //optimizer loop (if not whole distance is not possible, move only a bit!!!!)
+- while (loci <= 5 && !moveisok)
++ while (loci <= 5 && !moveisok)
+ {
+ loci ++;
+ mesh[pi].X() = origp.X() + (x.Get(1) * t1.X() + x.Get(2) * t2.X())*fact;
+@@ -951,11 +945,9 @@
+ fact = fact/2.;
+
+
+- ProjectPoint (surfi, mesh[pi]);
++ moveisok = ProjectPoint (surfi, mesh[pi], ngi);
+
+- moveisok = CalcPointGeomInfo(surfi, ngi, mesh[pi]);
+- // point lies on same chart in stlsurface
+-
++ // point lies on same chart in stlsurface
+ if (moveisok)
+ {
+ for (j = 0; j < locelements.Size(); j++)
+diff -Nru netgen-4.5_orig/libsrc/occ/occconstruction.cpp netgen-4.5_patch/libsrc/occ/occconstruction.cpp
+--- netgen-4.5_orig/libsrc/occ/occconstruction.cpp 2005-12-06 17:15:53.000000000 +0300
++++ netgen-4.5_patch/libsrc/occ/occconstruction.cpp 2006-04-27 13:12:54.000000000 +0400
+@@ -28,8 +28,8 @@
+ #include <BRepAlgoAPI_Common.hxx>
+ #include <BRepAlgoAPI_Fuse.hxx>
+ #include <BRepAlgoAPI_Section.hxx>
+-#include <BRepOffsetAPI_Sewing.hxx>
+-#include <BRepAlgo_Sewing.hxx>
++//#include <BRepOffsetAPI_Sewing.hxx>
++//#include <BRepAlgo_Sewing.hxx>
+ #include <BRepOffsetAPI_MakeOffsetShape.hxx>
+ #include <ShapeFix_Shape.hxx>
+ namespace netgen
+diff -Nru netgen-4.5_orig/libsrc/occ/occgenmesh.cpp netgen-4.5_patch/libsrc/occ/occgenmesh.cpp
+--- netgen-4.5_orig/libsrc/occ/occgenmesh.cpp 2006-02-07 12:12:48.000000000 +0300
++++ netgen-4.5_patch/libsrc/occ/occgenmesh.cpp 2006-10-25 16:14:48.000000000 +0400
+@@ -28,7 +28,7 @@
+ return Point<3> (p.X(), p.Y(), p.Z());
+ }
+
+- void DivideEdge (TopoDS_Edge & edge,
++ static void DivideEdge (TopoDS_Edge & edge,
+ ARRAY<MeshPoint> & ps,
+ ARRAY<double> & params,
+ Mesh & mesh)
+@@ -49,23 +49,19 @@
+ hvalue[0] = 0;
+ pnt = c->Value(s0);
+
+- double olddist = 0;
+- double dist = 0;
+-
+- for (int i = 1; i <= DIVIDEEDGESECTIONS; i++)
++ int i;
++ for (i = 1; i <= DIVIDEEDGESECTIONS; i++)
+ {
+ oldpnt = pnt;
+ pnt = c->Value(s0+(i/double(DIVIDEEDGESECTIONS))*(s1-s0));
++ double dist = pnt.Distance(oldpnt);
+ hvalue[i] = hvalue[i-1] +
+ 1.0/mesh.GetH(Point3d(pnt.X(), pnt.Y(), pnt.Z()))*
+- pnt.Distance(oldpnt);
++ dist;
+
+ //(*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;
+
+-
+- olddist = dist;
+- dist = pnt.Distance(oldpnt);
+ }
+
+ // nsubedges = int(ceil(hvalue[DIVIDEEDGESECTIONS]));
+@@ -74,7 +70,7 @@
+ ps.SetSize(nsubedges-1);
+ params.SetSize(nsubedges+1);
+
+- int i = 1;
++ i = 1;
+ int i1 = 0;
+ do
+ {
+@@ -112,7 +108,7 @@
+
+ static void FindEdges (OCCGeometry & geom, Mesh & mesh)
+ {
+- char * savetask = multithread.task;
++ const char * savetask = multithread.task;
+ multithread.task = "Edge meshing";
+
+ (*testout) << "edge meshing" << endl;
+@@ -124,6 +120,7 @@
+ (*testout) << "nedges = " << nedges << endl;
+
+ double eps = 1e-6 * geom.GetBoundingBox().Diam();
++ double eps2 = eps * eps;
+
+ for (int i = 1; i <= nvertices; i++)
+ {
+@@ -133,7 +130,7 @@
+ bool exists = 0;
+ if (merge_solids)
+ for (PointIndex pi = 1; pi <= mesh.GetNP(); pi++)
+- if ( Dist2 (mesh[pi], Point<3>(mp)) < eps*eps)
++ if ( Dist2 (mesh[pi], Point<3>(mp)) < eps2)
+ {
+ exists = 1;
+ break;
+@@ -276,8 +273,8 @@
+ pnums.Last() = -1;
+ for (PointIndex pi = 1; pi < first_ep; pi++)
+ {
+- 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) < eps2) pnums[0] = pi;
++ if (Dist2 (mesh[pi], lp) < eps2) pnums.Last() = pi;
+ }
+ }
+
+@@ -287,7 +284,7 @@
+ bool exists = 0;
+ int j;
+ for (j = first_ep; j <= mesh.GetNP(); j++)
+- if ((mesh.Point(j)-Point<3>(mp[i-1])).Length() < eps)
++ if (Dist2(mesh.Point(j), Point<3>(mp[i-1])) < eps2)
+ {
+ exists = 1;
+ break;
+@@ -394,7 +391,7 @@
+ int i, j, k;
+ int changed;
+
+- char * savetask = multithread.task;
++ const char * savetask = multithread.task;
+ multithread.task = "Surface meshing";
+
+ geom.facemeshstatus = 0;
+@@ -751,7 +748,7 @@
+ multithread.task = savetask;
+ }
+
+- double ComputeH (double kappa)
++ static double ComputeH (double kappa)
+ {
+ double hret;
+ kappa *= mparam.curvaturesafety;
+@@ -779,7 +776,7 @@
+ double nq = n*q;
+
+ Point<3> p = p0 + 0.5*n;
+- double lambda = (p-l.p0)*n / nq;
++ double lambda = (fabs(nq) > 1e-10 ? (p-l.p0)*n / nq : -1);
+
+ if (lambda >= 0 && lambda <= 1)
+ {
+@@ -799,55 +796,55 @@
+
+
+
+- void RestrictHTriangle (gp_Pnt2d & par0, gp_Pnt2d & par1, gp_Pnt2d & par2,
+- BRepLProp_SLProps * prop, Mesh & mesh, const double maxside, int depth, double h = 0)
++ static void RestrictHTriangle (gp_Pnt2d & par0, gp_Pnt2d & par1, gp_Pnt2d & par2,
++ BRepAdaptor_Surface& surf, Mesh & mesh, const double maxside, int depth, double h = 0)
+ {
+-
++ BRepLProp_SLProps prop(surf, 2, 1e-5);
+
+ gp_Pnt2d parmid;
+
+ parmid.SetX(0.3*(par0.X()+par1.X()+par2.X()));
+ parmid.SetY(0.3*(par0.Y()+par1.Y()+par2.Y()));
+
+- if (depth == 0)
++ //if (depth == 0)
+ {
+ double curvature = 0;
+
+- prop->SetParameters (parmid.X(), parmid.Y());
+- if (!prop->IsCurvatureDefined())
++ prop.SetParameters (parmid.X(), parmid.Y());
++ if (!prop.IsCurvatureDefined())
+ {
+ (*testout) << "curvature not defined!" << endl;
+ return;
+ }
+- curvature = max(fabs(prop->MinCurvature()),
+- fabs(prop->MaxCurvature()));
++ curvature = max(fabs(prop.MinCurvature()),
++ fabs(prop.MaxCurvature()));
+
+- prop->SetParameters (par0.X(), par0.Y());
+- if (!prop->IsCurvatureDefined())
++ prop.SetParameters (par0.X(), par0.Y());
++ if (!prop.IsCurvatureDefined())
+ {
+ (*testout) << "curvature not defined!" << endl;
+ return;
+ }
+- curvature = max(curvature,max(fabs(prop->MinCurvature()),
+- fabs(prop->MaxCurvature())));
++ curvature = max(curvature,max(fabs(prop.MinCurvature()),
++ fabs(prop.MaxCurvature())));
+
+- prop->SetParameters (par1.X(), par1.Y());
+- if (!prop->IsCurvatureDefined())
++ prop.SetParameters (par1.X(), par1.Y());
++ if (!prop.IsCurvatureDefined())
+ {
+ (*testout) << "curvature not defined!" << endl;
+ return;
+ }
+- curvature = max(curvature,max(fabs(prop->MinCurvature()),
+- fabs(prop->MaxCurvature())));
++ curvature = max(curvature,max(fabs(prop.MinCurvature()),
++ fabs(prop.MaxCurvature())));
+
+- prop->SetParameters (par2.X(), par2.Y());
+- if (!prop->IsCurvatureDefined())
++ prop.SetParameters (par2.X(), par2.Y());
++ if (!prop.IsCurvatureDefined())
+ {
+ (*testout) << "curvature not defined!" << endl;
+ return;
+ }
+- curvature = max(curvature,max(fabs(prop->MinCurvature()),
+- fabs(prop->MaxCurvature())));
++ curvature = max(curvature,max(fabs(prop.MinCurvature()),
++ fabs(prop.MaxCurvature())));
+
+ //(*testout) << "curvature " << curvature << endl;
+
+@@ -886,51 +883,47 @@
+ pm1.SetX(0.5*(par0.X()+par2.X())); pm1.SetY(0.5*(par0.Y()+par2.Y()));
+ pm2.SetX(0.5*(par1.X()+par0.X())); pm2.SetY(0.5*(par1.Y()+par0.Y()));
+
+- RestrictHTriangle (pm0, pm1, pm2, prop, mesh, 0.5*maxside, depth+1, h);
+- RestrictHTriangle (par0, pm1, pm2, prop, mesh, 0.5*maxside, depth+1, h);
+- RestrictHTriangle (par1, pm0, pm2, prop, mesh, 0.5*maxside, depth+1, h);
+- RestrictHTriangle (par2, pm1, pm0, prop, mesh, 0.5*maxside, depth+1, h);
++ RestrictHTriangle (pm0, pm1, pm2, surf, mesh, 0.5*maxside, depth+1, h);
++ RestrictHTriangle (par0, pm1, pm2, surf, mesh, 0.5*maxside, depth+1, h);
++ RestrictHTriangle (par1, pm0, pm2, surf, mesh, 0.5*maxside, depth+1, h);
++ RestrictHTriangle (par2, pm1, pm0, surf, mesh, 0.5*maxside, depth+1, h);
+ }
+ else
+ {
+ gp_Pnt pnt;
+ Point3d p3d;
+
+- prop->SetParameters (parmid.X(), parmid.Y());
+- pnt = prop->Value();
++ surf.D0(parmid.X(), parmid.Y(), pnt);
+ p3d = Point3d(pnt.X(), pnt.Y(), pnt.Z());
+ mesh.RestrictLocalH (p3d, h);
+
+
+- prop->SetParameters (par0.X(), par0.Y());
+- pnt = prop->Value();
++ surf.D0(par0.X(), par0.Y(), pnt);
+ p3d = Point3d(pnt.X(), pnt.Y(), pnt.Z());
+ mesh.RestrictLocalH (p3d, h);
+
+- prop->SetParameters (par1.X(), par1.Y());
+- pnt = prop->Value();
++ surf.D0(par1.X(), par1.Y(), pnt);
+ p3d = Point3d(pnt.X(), pnt.Y(), pnt.Z());
+ mesh.RestrictLocalH (p3d, h);
+
+- prop->SetParameters (par2.X(), par2.Y());
+- pnt = prop->Value();
++ surf.D0(par2.X(), par2.Y(), pnt);
+ p3d = Point3d(pnt.X(), pnt.Y(), pnt.Z());
+ mesh.RestrictLocalH (p3d, h);
+
+- (*testout) << "p = " << p3d << ", h = " << h << ", maxside = " << maxside << endl;
++ //(*testout) << "p = " << p3d << ", h = " << h << ", maxside = " << maxside << endl;
+ /*
+ (*testout) << pnt.X() << " " << pnt.Y() << " " << pnt.Z() << endl;
+
+- prop->SetParameters (par0.X(), par0.Y());
+- pnt = prop->Value();
++ prop.SetParameters (par0.X(), par0.Y());
++ pnt = prop.Value();
+ (*testout) << pnt.X() << " " << pnt.Y() << " " << pnt.Z() << endl;
+
+- prop->SetParameters (par1.X(), par1.Y());
+- pnt = prop->Value();
++ prop.SetParameters (par1.X(), par1.Y());
++ pnt = prop.Value();
+ (*testout) << pnt.X() << " " << pnt.Y() << " " << pnt.Z() << endl;
+
+- prop->SetParameters (par2.X(), par2.Y());
+- pnt = prop->Value();
++ prop.SetParameters (par2.X(), par2.Y());
++ pnt = prop.Value();
+ (*testout) << pnt.X() << " " << pnt.Y() << " " << pnt.Z() << endl;
+ */
+ }
+@@ -970,7 +963,7 @@
+ if (mparam.uselocalh)
+ {
+
+- char * savetask = multithread.task;
++ const char * savetask = multithread.task;
+ multithread.percent = 0;
+
+ mesh->SetLocalH (bb.PMin(), bb.PMax(), mparam.grading);
+@@ -1075,7 +1068,6 @@
+ if (triangulation.IsNull()) continue;
+
+ BRepAdaptor_Surface sf(face, Standard_True);
+- BRepLProp_SLProps prop(sf, 2, 1e-5);
+
+ int ntriangles = triangulation -> NbTriangles();
+ for (int j = 1; j <= ntriangles; j++)
+@@ -1096,7 +1088,7 @@
+ maxside = max (maxside, p[1].Distance(p[2]));
+ //cout << "\rFace " << i << " pos11 ntriangles " << ntriangles << " maxside " << maxside << flush;
+
+- RestrictHTriangle (par[0], par[1], par[2], &prop, *mesh, maxside, 0);
++ RestrictHTriangle (par[0], par[1], par[2], sf, *mesh, maxside, 0);
+ //cout << "\rFace " << i << " pos12 ntriangles " << ntriangles << flush;
+ }
+ }
+diff -Nru netgen-4.5_orig/libsrc/occ/occgeom.cpp netgen-4.5_patch/libsrc/occ/occgeom.cpp
+--- netgen-4.5_orig/libsrc/occ/occgeom.cpp 2006-01-25 15:35:50.000000000 +0300
++++ netgen-4.5_patch/libsrc/occ/occgeom.cpp 2006-10-25 16:15:24.000000000 +0400
+@@ -7,6 +7,8 @@
+ #include "ShapeAnalysis_ShapeContents.hxx"
+ #include "ShapeAnalysis_CheckSmallFace.hxx"
+ #include "ShapeAnalysis_DataMapOfShapeListOfReal.hxx"
++#include <ShapeAnalysis_Surface.hxx>
++#include <BRepTopAdaptor_FClass2d.hxx>
+ #include "BRepAlgoAPI_Fuse.hxx"
+ #include "BRepCheck_Analyzer.hxx"
+ #include "BRepLib.hxx"
+@@ -16,11 +18,19 @@
+ #include "Partition_Spliter.hxx"
+ //#include "VrmlAPI.hxx"
+ //#include "StlAPI.hxx"
++#include <TopAbs_State.hxx>
+
+
+ namespace netgen
+ {
+
++ OCCGeometry::~OCCGeometry()
++ {
++ NCollection_DataMap<int,BRepTopAdaptor_FClass2d*>::Iterator it(fclsmap);
++ for (; it.More(); it.Next())
++ delete it.Value();
++ }
++
+ void OCCGeometry :: PrintNrShapes ()
+ {
+ TopExp_Explorer e;
+@@ -947,13 +957,13 @@
+
+ void OCCGeometry :: BuildVisualizationMesh ()
+ {
+-
+- cout << "Preparing visualization (deflection = " << vispar.occdeflection << ") ... " << flush;
++ double vispar_occdeflection = 0.01;
++ cout << "Preparing visualization (deflection = " << vispar_occdeflection << ") ... " << flush;
+
+
+ BRepTools::Clean (shape);
+ //WriteOCC_STL("test.stl");
+- BRepMesh_IncrementalMesh::BRepMesh_IncrementalMesh (shape, vispar.occdeflection, true);
++ BRepMesh_IncrementalMesh::BRepMesh_IncrementalMesh (shape, vispar_occdeflection, true);
+ cout << "done" << endl;
+
+
+@@ -973,8 +983,27 @@
+
+ }
+
++ void OCCGeometry::GetFaceTools(int surfi, Handle(ShapeAnalysis_Surface)& proj,
++ BRepTopAdaptor_FClass2d*& cls) const
++ {
++ //MSV: organize caching projector in the map
++ if (fprjmap.IsBound(surfi))
++ {
++ proj = fprjmap.Find(surfi);
++ cls = fclsmap.Find(surfi);
++ }
++ else
++ {
++ const TopoDS_Face& aFace = TopoDS::Face(fmap(surfi));
++ Handle(Geom_Surface) aSurf = BRep_Tool::Surface(aFace);
++ proj = new ShapeAnalysis_Surface(aSurf);
++ fprjmap.Bind(surfi, proj);
++ cls = new BRepTopAdaptor_FClass2d(aFace,Precision::Confusion());
++ fclsmap.Bind(surfi, cls);
++ }
++ }
+
+- void OCCGeometry :: Project (int surfi, Point<3> & p) const
++ bool OCCGeometry :: Project (int surfi, Point<3> & p, double& u, double& v) const
+ {
+ static int cnt = 0;
+ if (++cnt % 1000 == 0) cout << "Project cnt = " << cnt << endl;
+@@ -983,18 +1012,22 @@
+
+ //(*testout) << "before " << pnt.X() << " "<< pnt.Y() << " "<< pnt.Z() << " " << endl;
+
+- GeomAPI_ProjectPointOnSurf proj(pnt, BRep_Tool::Surface(TopoDS::Face(fmap(surfi))));
+- if (proj.NbPoints() == 0)
+- {
+- cout << "Projection fails" << endl;
+- }
+- else
+- {
+- pnt = proj.NearestPoint();
+- //(*testout) << "after " << pnt.X() << " "<< pnt.Y() << " "<< pnt.Z() << " " << endl;
++ Handle(ShapeAnalysis_Surface) proj;
++ BRepTopAdaptor_FClass2d *cls;
++ GetFaceTools(surfi, proj, cls);
+
+- p = Point<3> (pnt.X(), pnt.Y(), pnt.Z());
+- }
++ gp_Pnt2d p2d = proj->ValueOfUV(pnt, Precision::Confusion());
++ if (cls->Perform(p2d) == TopAbs_OUT)
++ {
++ //cout << "Projection fails" << endl;
++ return false;
++ }
++ pnt = proj->Value(p2d);
++ p2d.Coord(u, v);
++ //(*testout) << "after " << pnt.X() << " "<< pnt.Y() << " "<< pnt.Z() << " " << endl;
++
++ p = Point<3> (pnt.X(), pnt.Y(), pnt.Z());
++ return true;
+ }
+
+
+@@ -1002,54 +1035,20 @@
+ {
+ gp_Pnt p(ap(0), ap(1), ap(2));
+
+- Handle(Geom_Surface) surface = BRep_Tool::Surface(TopoDS::Face(fmap(surfi)));
++ Handle(ShapeAnalysis_Surface) proj;
++ BRepTopAdaptor_FClass2d *cls;
++ GetFaceTools(surfi, proj, cls);
+
+- gp_Pnt x = surface->Value (u,v);
+-
+- if (p.SquareDistance(x) <= sqr(PROJECTION_TOLERANCE)) return true;
+-
+- gp_Vec du, dv;
+-
+- surface->D1(u,v,x,du,dv);
+-
+- int count = 0;
+-
+- gp_Pnt xold;
+- gp_Vec n;
+- double det, lambda, mu;
+-
+- do {
+- count++;
+-
+- n = du^dv;
+-
+- det = Det3 (n.X(), du.X(), dv.X(),
+- n.Y(), du.Y(), dv.Y(),
+- n.Z(), du.Z(), dv.Z());
+-
+- if (det < 1e-15) return false;
+-
+- lambda = Det3 (n.X(), p.X()-x.X(), dv.X(),
+- n.Y(), p.Y()-x.Y(), dv.Y(),
+- n.Z(), p.Z()-x.Z(), dv.Z())/det;
+-
+- mu = Det3 (n.X(), du.X(), p.X()-x.X(),
+- n.Y(), du.Y(), p.Y()-x.Y(),
+- n.Z(), du.Z(), p.Z()-x.Z())/det;
+-
+- u += lambda;
+- v += mu;
+-
+- xold = x;
+- surface->D1(u,v,x,du,dv);
+-
+- } while (xold.SquareDistance(x) > sqr(PROJECTION_TOLERANCE) && count < 50);
+-
+- // (*testout) << "FastProject count: " << count << endl;
+-
+- if (count == 50) return false;
++ gp_Pnt2d p2d = proj->NextValueOfUV(gp_Pnt2d(u,v), p, Precision::Confusion());
++ if (cls->Perform(p2d) == TopAbs_OUT)
++ {
++ //cout << "Projection fails" << endl;
++ return false;
++ }
+
+- ap = Point<3> (x.X(), x.Y(), x.Z());
++ p = proj->Value(p2d);
++ p2d.Coord(u, v);
++ ap = Point<3> (p.X(), p.Y(), p.Z());
+
+ return true;
+ }
+diff -Nru netgen-4.5_orig/libsrc/occ/occgeom.hpp netgen-4.5_patch/libsrc/occ/occgeom.hpp
+--- netgen-4.5_orig/libsrc/occ/occgeom.hpp 2006-01-25 15:35:50.000000000 +0300
++++ netgen-4.5_patch/libsrc/occ/occgeom.hpp 2006-10-25 16:16:01.000000000 +0400
+@@ -15,8 +15,6 @@
+ #include "Geom_Curve.hxx"
+ #include "Geom2d_Curve.hxx"
+ #include "Geom_Surface.hxx"
+-#include "GeomAPI_ProjectPointOnSurf.hxx"
+-#include "GeomAPI_ProjectPointOnCurve.hxx"
+ #include "BRepTools.hxx"
+ #include "TopExp.hxx"
+ #include "BRepBuilderAPI_MakeVertex.hxx"
+@@ -41,8 +39,6 @@
+ #include "Geom_Curve.hxx"
+ #include "Geom2d_Curve.hxx"
+ #include "Geom_Surface.hxx"
+-#include "GeomAPI_ProjectPointOnSurf.hxx"
+-#include "GeomAPI_ProjectPointOnCurve.hxx"
+ #include "TopoDS_Wire.hxx"
+ #include "BRepTools_WireExplorer.hxx"
+ #include "BRepTools.hxx"
+@@ -69,7 +65,7 @@
+ #include "IGESToBRep_Reader.hxx"
+ #include "Interface_Static.hxx"
+ #include "GeomAPI_ExtremaCurveCurve.hxx"
+-#include "Standard_ErrorHandler.hxx"
++//#include "Standard_ErrorHandler.hxx"
+ #include "Standard_Failure.hxx"
+ #include "ShapeUpgrade_ShellSewing.hxx"
+ #include "ShapeFix_Shape.hxx"
+@@ -84,11 +80,15 @@
+ #include "STEPControl_Writer.hxx"
+ #include "StlAPI_Writer.hxx"
+ #include "STEPControl_StepModelType.hxx"
++#include <NCollection_DataMap.hxx>
++
++class Handle_ShapeAnalysis_Surface;
++class BRepTopAdaptor_FClass2d;
+
+ namespace netgen
+ {
+
+-#include "../visualization/vispar.hpp"
++ //#include "../visualization/vispar.hpp"
+ // class VisualizationParameters;
+ // extern VisualizationParameters vispar;
+
+@@ -159,6 +159,8 @@
+ class OCCGeometry
+ {
+ Point<3> center;
++ mutable NCollection_DataMap<int,Handle_ShapeAnalysis_Surface> fprjmap;
++ mutable NCollection_DataMap<int,BRepTopAdaptor_FClass2d*> fclsmap;
+
+ public:
+ TopoDS_Shape shape;
+@@ -189,6 +191,7 @@
+ vmap.Clear();
+ }
+
++ ~OCCGeometry();
+
+ void BuildFMap();
+
+@@ -204,10 +207,12 @@
+ Point<3> Center()
+ { return center; }
+
+- void Project (int surfi, Point<3> & p) const;
++ bool Project (int surfi, Point<3> & p, double& u, double& v) const;
+ bool FastProject (int surfi, Point<3> & ap, double& u, double& v) const;
+
+-
++ void GetFaceTools(int surfi, Handle(ShapeAnalysis_Surface)& proj,
++ BRepTopAdaptor_FClass2d*& cls) const;
++
+ OCCSurface GetSurface (int surfi)
+ {
+ cout << "OCCGeometry::GetSurface using PLANESPACE" << endl;
+diff -Nru netgen-4.5_orig/libsrc/occ/occmeshsurf.cpp netgen-4.5_patch/libsrc/occ/occmeshsurf.cpp
+--- netgen-4.5_orig/libsrc/occ/occmeshsurf.cpp 2006-01-25 15:36:26.000000000 +0300
++++ netgen-4.5_patch/libsrc/occ/occmeshsurf.cpp 2006-10-25 16:16:26.000000000 +0400
+@@ -5,6 +5,8 @@
+ #include <occgeom.hpp>
+ #include <meshing.hpp>
+ #include <GeomLProp_SLProps.hxx>
++#include <GeomAPI_ProjectPointOnSurf.hxx>
++#include <GeomAPI_ProjectPointOnCurve.hxx>
+
+
+ namespace netgen
+@@ -411,11 +413,16 @@
+ }
+
+
+- void MeshOptimize2dOCCSurfaces :: ProjectPoint (INDEX surfind, Point3d & p) const
++ bool MeshOptimize2dOCCSurfaces :: ProjectPoint (INDEX surfind, Point3d & p, PointGeomInfo& gi) const
+ {
+ Point<3> hp = p;
+- geometry.Project (surfind, hp);
++ bool ok;
++ if (gi.trignum > 0)
++ ok = geometry.FastProject (surfind, hp, gi.u, gi.v);
++ else
++ ok = geometry.Project (surfind, hp, gi.u, gi.v);
+ p = hp;
++ return ok;
+ }
+
+ void MeshOptimize2dOCCSurfaces :: ProjectPoint2 (INDEX surfind, INDEX surfind2,
+@@ -506,38 +513,6 @@
+ }
+
+
+- int MeshOptimize2dOCCSurfaces ::
+- CalcPointGeomInfo(int surfind, PointGeomInfo& gi, const Point3d& p) const
+- {
+- Standard_Real u,v;
+-
+- gp_Pnt pnt(p.X(), p.Y(), p.Z());
+-
+- Handle(Geom_Surface) occface;
+- occface = BRep_Tool::Surface(TopoDS::Face(geometry.fmap(surfind)));
+-
+- GeomAPI_ProjectPointOnSurf proj(pnt, occface);
+-
+- if (proj.NbPoints() < 1)
+- {
+- cout << "ERROR: OCCSurface :: GetNormalVector: GeomAPI_ProjectPointOnSurf failed!"
+- << endl;
+- cout << p << endl;
+- return 0;
+- }
+-
+- proj.LowerDistanceParameters (u, v);
+-
+- gi.u = u;
+- gi.v = v;
+- return 1;
+- }
+-
+-
+-
+-
+-
+-
+ OCCRefinementSurfaces :: OCCRefinementSurfaces (const OCCGeometry & ageometry)
+ : Refinement(), geometry(ageometry)
+ {
+@@ -627,10 +602,11 @@
+ if (!geometry.FastProject (surfi, hnewp, u, v))
+ {
+ cout << "Fast projection to surface fails! Using OCC projection" << endl;
+- geometry.Project (surfi, hnewp);
++ double u, v;
++ geometry.Project (surfi, hnewp, u, v);
+ }
+
+- newgi.trignum = 1;
++ newgi.trignum = surfi;
+ }
+
+ newp = hnewp;
+@@ -653,14 +629,17 @@
+ hnewp = Point<3> (pnt.X(), pnt.Y(), pnt.Z());
+ newp = hnewp;
+ newgi = ap1;
+- };
++ }
+
+
+ void OCCRefinementSurfaces :: ProjectToSurface (Point<3> & p, int surfi)
+ {
+ if (surfi > 0)
+- geometry.Project (surfi, p);
+- };
++ {
++ double u, v;
++ geometry.Project (surfi, p, u, v);
++ }
++ }
+
+ void OCCRefinementSurfaces :: ProjectToSurface (Point<3> & p, int surfi, PointGeomInfo & gi)
+ {
+@@ -668,9 +647,10 @@
+ if (!geometry.FastProject (surfi, p, gi.u, gi.v))
+ {
+ cout << "Fast projection to surface fails! Using OCC projection" << endl;
+- geometry.Project (surfi, p);
++ double u, v;
++ geometry.Project (surfi, p, u, v);
+ }
+- };
++ }
+
+
+
+diff -Nru netgen-4.5_orig/libsrc/occ/occmeshsurf.hpp netgen-4.5_patch/libsrc/occ/occmeshsurf.hpp
+--- netgen-4.5_orig/libsrc/occ/occmeshsurf.hpp 2005-06-09 18:51:10.000000000 +0400
++++ netgen-4.5_patch/libsrc/occ/occmeshsurf.hpp 2006-10-25 16:17:22.000000000 +0400
+@@ -151,7 +151,7 @@
+ MeshOptimize2dOCCSurfaces (const OCCGeometry & ageometry);
+
+ ///
+- virtual void ProjectPoint (INDEX surfind, Point3d & p) const;
++ virtual bool ProjectPoint (INDEX surfind, Point3d & p, PointGeomInfo& gi) const;
+ ///
+ virtual void ProjectPoint2 (INDEX surfind, INDEX surfind2, Point3d & p) const;
+ ///
+@@ -159,9 +159,6 @@
+ ///
+ virtual void GetNormalVector(INDEX surfind, const Point3d & p, PointGeomInfo & gi, Vec3d & n) const;
+
+-
+- virtual int CalcPointGeomInfo(int surfind, PointGeomInfo& gi, const Point3d& p3) const;
+-
+ };
+
+
+diff -Nru netgen-4.5_orig/libsrc/stlgeom/meshstlsurface.cpp netgen-4.5_patch/libsrc/stlgeom/meshstlsurface.cpp
+--- netgen-4.5_orig/libsrc/stlgeom/meshstlsurface.cpp 2006-01-11 18:08:20.000000000 +0300
++++ netgen-4.5_patch/libsrc/stlgeom/meshstlsurface.cpp 2006-10-25 16:17:47.000000000 +0400
+@@ -946,20 +946,23 @@
+ }
+
+
+-void MeshOptimizeSTLSurface :: ProjectPoint (INDEX surfind, Point3d & p) const
++bool MeshOptimizeSTLSurface :: ProjectPoint (INDEX surfind, Point3d & p, PointGeomInfo& gi) const
+ {
+ Point<3> hp = p;
+- if (!geom.Project (hp))
++ if (gi.trignum > 0)
++ ((STLGeometry&)geom).SelectChartOfTriangle (gi.trignum);
++ if (!(gi.trignum = geom.Project (hp)))
+ {
+ PrintMessage(7,"project failed");
+
+- if (!geom.ProjectOnWholeSurface(hp))
++ if (!(gi.trignum = geom.ProjectOnWholeSurface(hp)))
+ {
+ PrintMessage(7, "project on whole surface failed");
+ }
+ }
+ p = hp;
+ // geometry.GetSurface(surfind)->Project (p);
++ return gi.trignum > 0;
+ }
+
+ void MeshOptimizeSTLSurface :: ProjectPoint2 (INDEX surfind, INDEX surfind2, Point3d & p) const
+@@ -970,20 +973,6 @@
+ */
+ }
+
+-int MeshOptimizeSTLSurface :: CalcPointGeomInfo(PointGeomInfo& gi, const Point3d& p3) const
+-{
+- Point<3> hp = p3;
+- gi.trignum = geom.Project (hp);
+-
+- if (gi.trignum)
+- {
+- return 1;
+- }
+-
+- return 0;
+-
+-}
+-
+ void MeshOptimizeSTLSurface :: GetNormalVector(INDEX surfind, const Point3d & p, Vec3d & n) const
+ {
+ n = geom.GetChartNormalVector();
+diff -Nru netgen-4.5_orig/libsrc/stlgeom/meshstlsurface.hpp netgen-4.5_patch/libsrc/stlgeom/meshstlsurface.hpp
+--- netgen-4.5_orig/libsrc/stlgeom/meshstlsurface.hpp 2004-09-30 17:13:56.000000000 +0400
++++ netgen-4.5_patch/libsrc/stlgeom/meshstlsurface.hpp 2006-10-25 16:17:59.000000000 +0400
+@@ -79,12 +79,10 @@
+ virtual void SelectSurfaceOfPoint (const Point3d & p,
+ const PointGeomInfo & gi);
+ ///
+- virtual void ProjectPoint (INDEX surfind, Point3d & p) const;
++ virtual bool ProjectPoint (INDEX surfind, Point3d & p, PointGeomInfo& gi) const;
+ ///
+ virtual void ProjectPoint2 (INDEX surfind, INDEX surfind2, Point3d & p) const;
+ ///
+- virtual int CalcPointGeomInfo(PointGeomInfo& gi, const Point3d& p3) const;
+- ///
+ virtual void GetNormalVector(INDEX surfind, const Point3d & p, Vec3d & n) const;
+ };
+
+diff -Nru netgen-4.5_orig/makeForSalome.sh netgen-4.5_patch/makeForSalome.sh
+--- netgen-4.5_orig/makeForSalome.sh 1970-01-01 03:00:00.000000000 +0300
++++ netgen-4.5_patch/makeForSalome.sh 2006-04-27 13:12:54.000000000 +0400
+@@ -0,0 +1,31 @@
++#! /bin/sh
++cp ngtcltk/ngnewdelete.* libsrc/interface/
++
++MACHINE=LINUX
++export MACHINE
++make -C libsrc/csg
++make -C libsrc/general
++make -C libsrc/geom2d
++make -C libsrc/gprim
++make -C libsrc/interface
++make -C libsrc/linalg
++make -C libsrc/meshing
++make -C libsrc/opti
++make -C libsrc/stlgeom
++make -C libsrc/occ
++
++if [ ! -d install ] ; then
++ mkdir install
++fi
++
++cp -r lib install/
++
++if [ ! -d install/include ] ; then
++ mkdir install/include
++fi
++
++cp libsrc/interface/nglib.h libsrc/general/*.hpp libsrc/csg/*.hpp libsrc/geom2d/*.hpp \
++ libsrc/gprim/*.hpp libsrc/linalg/*.hpp libsrc/meshing/*.hpp \
++ libsrc/occ/*.hpp libsrc/opti/*.hpp libsrc/include/mydefs.hpp \
++ libsrc/stlgeom/*.hpp libsrc/include/mystdlib.h \
++ install/include
--- /dev/null
+// Copyright (C) 2005 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
+// CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS, L3S, LJLL, MENSI
+//
+// This library is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 2.1 of the License.
+//
+// This library is distributed in the hope that it will be useful
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+// Lesser General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License along with this library; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+//
+// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+//
+// File : NETGENPlugin_NETGEN_2D_ONLY.cxx
+// Author : Edward AGAPOV (OCC)
+// Project : SALOME
+
+
+#include "NETGENPlugin_NETGEN_2D_ONLY.hxx"
+
+#include "NETGENPlugin_Mesher.hxx"
+
+#include "SMDS_MeshElement.hxx"
+#include "SMDS_MeshNode.hxx"
+#include "SMESHDS_Mesh.hxx"
+#include "SMESH_Comment.hxx"
+#include "SMESH_Gen.hxx"
+#include "SMESH_Mesh.hxx"
+#include "SMESH_MesherHelper.hxx"
+#include "StdMeshers_FaceSide.hxx"
+#include "StdMeshers_MaxElementArea.hxx"
+#include "StdMeshers_LengthFromEdges.hxx"
+#include "StdMeshers_QuadranglePreference.hxx"
+
+#include <Standard_ErrorHandler.hxx>
+#include <Standard_Failure.hxx>
+
+#include "utilities.h"
+
+#include <list>
+#include <vector>
+
+/*
+ Netgen include files
+*/
+namespace nglib {
+#include <nglib.h>
+}
+#define OCCGEOMETRY
+#include <occgeom.hpp>
+#include <meshing.hpp>
+namespace netgen {
+ extern int OCCGenerateMesh (OCCGeometry&, Mesh*&, int, int, char*);
+ extern MeshingParameters mparam;
+}
+
+using namespace std;
+using namespace netgen;
+using namespace nglib;
+
+//=============================================================================
+/*!
+ *
+ */
+//=============================================================================
+
+NETGENPlugin_NETGEN_2D_ONLY::NETGENPlugin_NETGEN_2D_ONLY(int hypId, int studyId,
+ SMESH_Gen* gen)
+ : SMESH_2D_Algo(hypId, studyId, gen)
+{
+ MESSAGE("NETGENPlugin_NETGEN_2D_ONLY::NETGENPlugin_NETGEN_2D_ONLY");
+ _name = "NETGEN_2D_ONLY";
+
+ _shapeType = (1 << TopAbs_FACE);// 1 bit /shape type
+
+ _compatibleHypothesis.push_back("MaxElementArea");
+ _compatibleHypothesis.push_back("LengthFromEdges");
+ _compatibleHypothesis.push_back("QuadranglePreference");
+
+ _hypMaxElementArea = 0;
+ _hypLengthFromEdges = 0;
+ _hypQuadranglePreference = 0;
+}
+
+//=============================================================================
+/*!
+ *
+ */
+//=============================================================================
+
+NETGENPlugin_NETGEN_2D_ONLY::~NETGENPlugin_NETGEN_2D_ONLY()
+{
+ MESSAGE("NETGENPlugin_NETGEN_2D_ONLY::~NETGENPlugin_NETGEN_2D_ONLY");
+}
+
+//=============================================================================
+/*!
+ *
+ */
+//=============================================================================
+
+bool NETGENPlugin_NETGEN_2D_ONLY::CheckHypothesis (SMESH_Mesh& aMesh,
+ const TopoDS_Shape& aShape,
+ Hypothesis_Status& aStatus)
+{
+ _hypMaxElementArea = 0;
+ _hypLengthFromEdges = 0;
+ _hypQuadranglePreference = 0;
+
+ aStatus = HYP_MISSING;
+
+ const list<const SMESHDS_Hypothesis*>& hyps = GetUsedHypothesis(aMesh, aShape, false);
+
+ if (hyps.empty()) return false; // can't work with no hypothesis
+
+ list<const SMESHDS_Hypothesis*>::const_iterator ith;
+ for (ith = hyps.begin(); ith != hyps.end(); ++ith )
+ {
+ const SMESHDS_Hypothesis* hyp = (*ith);
+
+ string hypName = hyp->GetName();
+
+ if ( hypName == "MaxElementArea")
+ _hypMaxElementArea = static_cast<const StdMeshers_MaxElementArea*> (hyp);
+ else if ( hypName == "LengthFromEdges" )
+ _hypLengthFromEdges = static_cast<const StdMeshers_LengthFromEdges*> (hyp);
+ else if ( hypName == "QuadranglePreference" )
+ _hypQuadranglePreference = static_cast<const StdMeshers_QuadranglePreference*>(hyp);
+ else {
+ aStatus = HYP_INCOMPATIBLE;
+ return false;
+ }
+ }
+
+ if ( _hypMaxElementArea && _hypLengthFromEdges ) {
+ aStatus = HYP_CONCURENT;
+ return false;
+ }
+
+ if ( _hypMaxElementArea || _hypLengthFromEdges )
+ aStatus = HYP_OK;
+
+ return aStatus == HYP_OK;
+}
+
+//================================================================================
+/*!
+ * \brief Fill netgen mesh with segments
+ * \retval SMESH_ComputeErrorPtr - error description
+ */
+//================================================================================
+
+static TError AddSegmentsToMesh(netgen::Mesh& ngMesh,
+ OCCGeometry& geom,
+ const TSideVector& wires,
+ SMESH_MesherHelper& helper,
+ vector< const SMDS_MeshNode* > & nodeVec)
+{
+ // ----------------------------
+ // Check wires and count nodes
+ // ----------------------------
+ int nbNodes = 0;
+ for ( int iW = 0; iW < wires.size(); ++iW )
+ {
+ StdMeshers_FaceSidePtr wire = wires[ iW ];
+ if ( wire->MissVertexNode() )
+ return TError
+ (new SMESH_ComputeError(COMPERR_BAD_INPUT_MESH, "Missing nodes on vertices"));
+
+ const vector<UVPtStruct>& uvPtVec = wire->GetUVPtStruct();
+ if ( uvPtVec.size() != wire->NbPoints() )
+ return TError
+ (new SMESH_ComputeError(COMPERR_BAD_INPUT_MESH,
+ SMESH_Comment("Unexpected nb of points on wire ") << iW
+ << ": " << uvPtVec.size()<<" != "<<wire->NbPoints()));
+ nbNodes += wire->NbSegments();
+ }
+ nodeVec.reserve( nbNodes );
+
+ // -----------------
+ // Fill netgen mesh
+ // -----------------
+
+// netgen::Box<3> bb = geom.GetBoundingBox();
+// bb.Increase (bb.Diam()/10);
+// ngMesh.SetLocalH (bb.PMin(), bb.PMax(), 0.5); // set grading
+
+ const int faceID = 1, solidID = 0;
+ ngMesh.AddFaceDescriptor (FaceDescriptor(faceID, solidID, solidID, 0));
+
+ for ( int iW = 0; iW < wires.size(); ++iW )
+ {
+ StdMeshers_FaceSidePtr wire = wires[ iW ];
+ const vector<UVPtStruct>& uvPtVec = wire->GetUVPtStruct();
+
+ int firstPointID = ngMesh.GetNP() + 1;
+ int edgeID = 1, posID = -2;
+ for ( int i = 0; i < wire->NbSegments(); ++i ) // loop on segments
+ {
+ // Add the first point of a segment
+ const SMDS_MeshNode * n = uvPtVec[ i ].node;
+ const int posShapeID = n->GetPosition()->GetShapeId();
+
+ // skip nodes on degenerated edges
+ if ( helper.IsDegenShape( posShapeID ) &&
+ helper.IsDegenShape( uvPtVec[ i+1 ].node->GetPosition()->GetShapeId() ))
+ continue;
+
+ nodeVec.push_back( n );
+
+ MeshPoint mp( Point<3> (n->X(), n->Y(), n->Z()) );
+ ngMesh.AddPoint ( mp, 1, EDGEPOINT );
+
+ // Add the segment
+ Segment seg;
+
+ seg.p1 = ngMesh.GetNP(); // ng node id
+ seg.p2 = seg.p1 + 1; // ng node id
+ seg.edgenr = ngMesh.GetNSeg() + 1;// segment id
+ seg.si = faceID; // = geom.fmap.FindIndex (face);
+
+ for ( int iEnd = 0; iEnd < 2; ++iEnd)
+ {
+ const UVPtStruct& pnt = uvPtVec[ i + iEnd ];
+
+ seg.epgeominfo[ iEnd ].dist = pnt.param; // param on curve
+ seg.epgeominfo[ iEnd ].u = pnt.u;
+ seg.epgeominfo[ iEnd ].v = pnt.v;
+
+ // find out edge id and node parameter on edge
+ bool onVertex = ( pnt.node->GetPosition()->GetTypeOfPosition() == SMDS_TOP_VERTEX );
+ if ( onVertex || posShapeID != posID )
+ {
+ // get edge id
+ double normParam = pnt.normParam;
+ if ( onVertex )
+ normParam = 0.5 * ( uvPtVec[ i ].normParam + uvPtVec[ i+1 ].normParam );
+ const TopoDS_Edge& edge = wire->Edge( wire->EdgeIndex( normParam ));
+ edgeID = geom.emap.FindIndex( edge );
+ posID = posShapeID;
+ if ( onVertex ) // param on curve is different on each of two edges
+ seg.epgeominfo[ iEnd ].dist = helper.GetNodeU( edge, pnt.node );
+ }
+ seg.epgeominfo[ iEnd ].edgenr = edgeID; // = geom.emap.FindIndex(edge);
+ }
+
+ ngMesh.AddSegment (seg);
+
+// cout << "Segment: " << seg.edgenr << endl
+// << "\tp1: " << seg.p1 << endl
+// << "\tp2: " << seg.p2 << endl
+// << "\tp0 param: " << seg.epgeominfo[ 0 ].dist << endl
+// << "\tp0 uv: " << seg.epgeominfo[ 0 ].u <<", "<< seg.epgeominfo[ 0 ].v << endl
+// << "\tp0 edge: " << seg.epgeominfo[ 0 ].edgenr << endl
+// << "\tp1 param: " << seg.epgeominfo[ 1 ].dist << endl
+// << "\tp1 uv: " << seg.epgeominfo[ 1 ].u <<", "<< seg.epgeominfo[ 1 ].v << endl
+// << "\tp1 edge: " << seg.epgeominfo[ 1 ].edgenr << endl;
+ }
+ Segment& seg = ngMesh.LineSegment( ngMesh.GetNSeg() );
+ seg.p2 = firstPointID;
+ }
+
+ ngMesh.CalcSurfacesOfNode();
+
+ return TError();
+}
+
+//=============================================================================
+/*!
+ *Here we are going to use the NETGEN mesher
+ */
+//=============================================================================
+
+bool NETGENPlugin_NETGEN_2D_ONLY::Compute(SMESH_Mesh& aMesh,
+ const TopoDS_Shape& aShape)
+{
+ MESSAGE("NETGENPlugin_NETGEN_2D_ONLY::Compute()");
+
+ SMESHDS_Mesh* meshDS = aMesh.GetMeshDS();
+ int faceID = meshDS->ShapeToIndex( aShape );
+
+ SMESH_MesherHelper helper(aMesh);
+ _quadraticMesh = helper.IsQuadraticSubMesh(aShape);
+ helper.SetElementsOnShape( true );
+ const bool ignoreMediumNodes = _quadraticMesh;
+
+ // ------------------------
+ // get all edges of a face
+ // ------------------------
+ const TopoDS_Face F = TopoDS::Face( aShape.Oriented( TopAbs_FORWARD ));
+ TError problem;
+ TSideVector wires = StdMeshers_FaceSide::GetFaceWires( F, aMesh, ignoreMediumNodes, problem );
+ if ( problem && !problem->IsOK() )
+ return error( problem );
+ int nbWires = wires.size();
+ if ( nbWires == 0 )
+ return error( "Problem in StdMeshers_FaceSide::GetFaceWires()");
+ if ( wires[0]->NbSegments() < 3 ) // ex: a circle with 2 segments
+ return error(COMPERR_BAD_INPUT_MESH,
+ SMESH_Comment("Too few segments: ")<<wires[0]->NbSegments());
+
+ // -------------------------
+ // Make input netgen mesh
+ // -------------------------
+
+ Ng_Init();
+ netgen::Mesh * ngMesh = new netgen::Mesh ();
+
+ netgen::OCCGeometry occgeo;
+ NETGENPlugin_Mesher::PrepareOCCgeometry( occgeo, F );
+
+ vector< const SMDS_MeshNode* > nodeVec;
+ problem = AddSegmentsToMesh( *ngMesh, occgeo, wires, helper, nodeVec );
+ if ( problem && !problem->IsOK() ) {
+ delete ngMesh; Ng_Exit();
+ return error( problem );
+ }
+
+ // --------------------
+ // compute edge length
+ // --------------------
+
+ double edgeLength = 0;
+ if (_hypLengthFromEdges)
+ {
+ int nbSegments = 0;
+ for ( int iW = 0; iW < nbWires; ++iW )
+ {
+ edgeLength += wires[ iW ]->Length();
+ nbSegments += wires[ iW ]->NbSegments();
+ }
+ if ( nbSegments )
+ edgeLength /= nbSegments;
+ }
+ if ( _hypMaxElementArea )
+ {
+ double maxArea = _hypMaxElementArea->GetMaxArea();
+ edgeLength = sqrt(2. * maxArea/sqrt(3.0));
+ }
+ if ( edgeLength < DBL_MIN )
+ edgeLength = occgeo.GetBoundingBox().Diam();
+
+ //cout << " edgeLength = " << edgeLength << endl;
+
+ netgen::mparam.maxh = edgeLength;
+ netgen::mparam.quad = _hypQuadranglePreference ? 1 : 0;
+ //ngMesh->SetGlobalH ( edgeLength );
+
+ // -------------------------
+ // Generate surface mesh
+ // -------------------------
+
+ char *optstr;
+ int startWith = MESHCONST_MESHSURFACE;
+ int endWith = MESHCONST_OPTSURFACE;
+ int err = 1;
+
+ try {
+#if (OCC_VERSION_MAJOR << 16 | OCC_VERSION_MINOR << 8 | OCC_VERSION_MAINTENANCE) > 0x060100
+ OCC_CATCH_SIGNALS;
+#endif
+ err = netgen::OCCGenerateMesh(occgeo, ngMesh, startWith, endWith, optstr);
+ }
+ catch (Standard_Failure& ex) {
+ string comment = ex.DynamicType()->Name();
+ if ( ex.GetMessageString() && strlen( ex.GetMessageString() )) {
+ comment += ": ";
+ comment += ex.GetMessageString();
+ }
+ error(COMPERR_OCC_EXCEPTION, comment);
+ }
+ catch (NgException exc) {
+ error( SMESH_Comment("NgException: ") << exc.What() );
+ }
+ catch (...) {
+ error(COMPERR_EXCEPTION,"Exception in netgen::OCCGenerateMesh()");
+ }
+
+ // ----------------------------------------------------
+ // Fill the SMESHDS with the generated nodes and faces
+ // ----------------------------------------------------
+
+ int nbNodes = ngMesh->GetNP();
+ int nbFaces = ngMesh->GetNSE();
+
+ int nbInputNodes = nodeVec.size();
+ nodeVec.resize( nbNodes, 0 );
+
+ // add nodes
+ for ( int i = nbInputNodes + 1; i <= nbNodes; ++i )
+ {
+ const MeshPoint& ngPoint = ngMesh->Point(i);
+ SMDS_MeshNode * node = meshDS->AddNode(ngPoint.X(), ngPoint.Y(), ngPoint.Z());
+ nodeVec[ i-1 ] = node;
+ }
+
+ // create faces
+ bool reverse = ( aShape.Orientation() == TopAbs_REVERSED );
+ for ( int i = 1; i <= nbFaces ; ++i )
+ {
+ const Element2d& elem = ngMesh->SurfaceElement(i);
+ vector<const SMDS_MeshNode*> nodes( elem.GetNP() );
+ for (int j=1; j <= elem.GetNP(); ++j)
+ {
+ int pind = elem.PNum(j);
+ const SMDS_MeshNode* node = nodeVec.at(pind-1);
+ if ( reverse )
+ nodes[ nodes.size()-j ] = node;
+ else
+ nodes[ j-1 ] = node;
+ if ( node->GetPosition()->GetTypeOfPosition() == SMDS_TOP_3DSPACE )
+ {
+ const PointGeomInfo& pgi = elem.GeomInfoPi(j);
+ meshDS->SetNodeOnFace((SMDS_MeshNode*)node, faceID, pgi.u, pgi.v);
+ }
+ }
+ SMDS_MeshFace* face = 0;
+ if ( elem.GetType() == TRIG )
+ face = helper.AddFace(nodes[0],nodes[1],nodes[2]);
+ else
+ face = helper.AddFace(nodes[0],nodes[1],nodes[2],nodes[3]);
+ }
+
+ Ng_DeleteMesh((nglib::Ng_Mesh*)ngMesh);
+ Ng_Exit();
+
+ NETGENPlugin_Mesher::RemoveTmpFiles();
+
+ return !err;
+}
