* Nglib now supports surface and volume meshing of OCC geometry

* OCC code cleanup
* Added to Nglib source code documentation

libsrc/occ/occgenmesh.cpp

@@ -17,10 +17,34 @@ namespace netgen
 #define IGNORECURVELENGTH 1e-4
 
 
-   extern OCCParameters occparam;
-
    bool merge_solids = 1;
 
+   double Line :: Dist (Line l)
+   {
+      Vec<3> n = p1-p0;
+      Vec<3> q = l.p1-l.p0;
+      double nq = n*q;
+
+      Point<3> p = p0 + 0.5*n;
+      double lambda = (p-l.p0)*n / nq;
+
+      if (lambda >= 0 && lambda <= 1)
+      {
+         double d = (p-l.p0-lambda*q).Length();
+         //        if (d < 1e-3) d = 1e99;
+         return d;
+      }
+      else
+         return 1e99;
+   }
+
+
+
+   double Line :: Length ()
+   {
+      return (p1-p0).Length();
+   }
+
 
 
    inline Point<3> occ2ng (const gp_Pnt & p)
@@ -30,6 +54,179 @@ namespace netgen
 
 
 
+   double ComputeH (double kappa)
+   {
+      double hret;
+      kappa *= mparam.curvaturesafety;
+
+      if (mparam.maxh * kappa < 1)
+         hret = mparam.maxh;
+      else
+         hret = 1 / kappa;
+
+      if (mparam.maxh < hret)
+         hret = mparam.maxh;
+
+      return (hret);
+   }
+
+
+
+
+   void RestrictHTriangle (gp_Pnt2d & par0, gp_Pnt2d & par1, gp_Pnt2d & par2,
+                           BRepLProp_SLProps * prop, Mesh & mesh, int depth, double h = 0)
+   {
+      int ls = -1;
+
+      gp_Pnt pnt0,pnt1,pnt2;
+
+      prop->SetParameters (par0.X(), par0.Y());
+      pnt0 = prop->Value();
+
+      prop->SetParameters (par1.X(), par1.Y());
+      pnt1 = prop->Value();
+
+      prop->SetParameters (par2.X(), par2.Y());
+      pnt2 = prop->Value();
+
+      double aux;
+      double maxside = pnt0.Distance(pnt1);
+      ls = 2;
+      aux = pnt1.Distance(pnt2);
+      if(aux > maxside)
+      {
+         maxside = aux;
+         ls = 0;
+      }
+      aux = pnt2.Distance(pnt0);
+      if(aux > maxside)
+      {
+         maxside = aux;
+         ls = 1;
+      }
+
+
+
+      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%3 == 0)
+      {
+         double curvature = 0;
+
+         prop->SetParameters (parmid.X(), parmid.Y());
+         if (!prop->IsCurvatureDefined())
+         {
+            (*testout) << "curvature not defined!" << endl;
+            return;
+         }
+         curvature = max(fabs(prop->MinCurvature()),
+            fabs(prop->MaxCurvature()));
+
+         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())));
+
+         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())));
+
+         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())));
+
+         //(*testout) << "curvature " << curvature << endl;
+
+         if (curvature < 1e-3)
+         {
+            //(*testout) << "curvature too small (" << curvature << ")!" << endl;
+            return;
+            // return war bis 10.2.05 auskommentiert
+         }
+
+
+
+         h = ComputeH (curvature+1e-10);
+
+         if(h < 1e-4*maxside)
+            return;
+
+
+         if (h > 30) return;
+      }
+
+      if (h < maxside && depth < 10)
+      {
+         //cout << "\r h " << h << flush;
+         gp_Pnt2d pm;
+
+         //cout << "h " << h << " maxside " << maxside << " depth " << depth << endl;
+         //cout << "par0 " << par0.X() << " " << par0.Y()
+         //<< " par1 " << par1.X() << " " << par1.Y()
+         //   << " par2 " << par2.X() << " " << par2.Y()<< endl;
+
+         if(ls == 0)
+         {
+            pm.SetX(0.5*(par1.X()+par2.X())); pm.SetY(0.5*(par1.Y()+par2.Y()));
+            RestrictHTriangle(pm, par2, par0, prop, mesh, depth+1, h);
+            RestrictHTriangle(pm, par0, par1, prop, mesh, depth+1, h);
+         }
+         else if(ls == 1)
+         {
+            pm.SetX(0.5*(par0.X()+par2.X())); pm.SetY(0.5*(par0.Y()+par2.Y()));
+            RestrictHTriangle(pm, par1, par2, prop, mesh, depth+1, h);
+            RestrictHTriangle(pm, par0, par1, prop, mesh, depth+1, h);
+         }
+         else if(ls == 2)
+         {
+            pm.SetX(0.5*(par0.X()+par1.X())); pm.SetY(0.5*(par0.Y()+par1.Y()));
+            RestrictHTriangle(pm, par1, par2, prop, mesh, depth+1, h);
+            RestrictHTriangle(pm, par2, par0, prop, mesh, depth+1, h);
+         }
+
+      }
+      else
+      {
+         gp_Pnt pnt;
+         Point3d p3d;
+
+         prop->SetParameters (parmid.X(), parmid.Y());
+         pnt = prop->Value();
+         p3d = Point3d(pnt.X(), pnt.Y(), pnt.Z());
+         mesh.RestrictLocalH (p3d, h);
+
+         p3d = Point3d(pnt0.X(), pnt0.Y(), pnt0.Z());
+         mesh.RestrictLocalH (p3d, h);
+
+         p3d = Point3d(pnt1.X(), pnt1.Y(), pnt1.Z());
+         mesh.RestrictLocalH (p3d, h);
+
+         p3d = Point3d(pnt2.X(), pnt2.Y(), pnt2.Z());
+         mesh.RestrictLocalH (p3d, h);
+
+         //(*testout) << "p = " << p3d << ", h = " << h << ", maxside = " << maxside << endl;
+
+      }
+   }
+
+
 
    void DivideEdge (TopoDS_Edge & edge, Array<MeshPoint> & ps,
                     Array<double> & params, Mesh & mesh)
@@ -53,7 +250,9 @@ namespace netgen
       double olddist = 0;
       double dist = 0;
 
-      for (int i = 1; i <= DIVIDEEDGESECTIONS; i++)
+      int tmpVal = (int)(DIVIDEEDGESECTIONS);
+
+      for (int i = 1; i <= tmpVal; i++)
       {
          oldpnt = pnt;
          pnt = c->Value(s0+(i/double(DIVIDEEDGESECTIONS))*(s1-s0));
@@ -64,7 +263,6 @@ namespace netgen
          //(*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);
       }
@@ -111,7 +309,7 @@ namespace netgen
 
 
 
-   static void FindEdges (OCCGeometry & geom, Mesh & mesh)
+   void OCCFindEdges (OCCGeometry & geom, Mesh & mesh)
    {
       const char * savetask = multithread.task;
       multithread.task = "Edge meshing";
@@ -180,8 +378,6 @@ namespace netgen
                total++;
 
 
-
-
       int facenr = 0;
       int edgenr = 0;
 
@@ -273,7 +469,6 @@ namespace netgen
 
                DivideEdge (edge, mp, params, mesh);
  
-
                Array <int> pnums;
                pnums.SetSize (mp.Size()+2);
 
@@ -404,7 +599,7 @@ namespace netgen
 
 
 
-   static void OCCMeshSurface (OCCGeometry & geom, Mesh & mesh, int perfstepsend)
+   void OCCMeshSurface (OCCGeometry & geom, Mesh & mesh, int perfstepsend)
    {
       int i, j, k;
       int changed;
@@ -774,367 +969,30 @@ namespace netgen
       multithread.task = savetask;
    }
 
-   double ComputeH (double kappa)
+
+
+   void OCCSetLocalMeshSize(OCCGeometry & geom, Mesh & mesh)
    {
-      double hret;
+      mesh.SetGlobalH (mparam.maxh);
-      kappa *= mparam.curvaturesafety;
+      mesh.SetMinimalH (mparam.minh);
-
-      if (mparam.maxh * kappa < 1)
-         hret = mparam.maxh;
-      else
-         hret = 1 / kappa;
-
-      if (mparam.maxh < hret)
-         hret = mparam.maxh;
-
-      return (hret);
-   }
-
-
-   class Line
-   {
-   public:
-      Point<3> p0, p1;
-      double Dist (Line l)
-      {
-         Vec<3> n = p1-p0;
-         Vec<3> q = l.p1-l.p0;
-         double nq = n*q;
-
-         Point<3> p = p0 + 0.5*n;
-         double lambda = (p-l.p0)*n / nq;
-
-         if (lambda >= 0 && lambda <= 1)
-         {
-            double d = (p-l.p0-lambda*q).Length();
-            //        if (d < 1e-3) d = 1e99;
-            return d;
-         }
-         else
-            return 1e99;
-      }
-
-      double Length ()
-      {
-         return (p1-p0).Length();
-      };
-   };
-
-   /*
-   void RestrictHTriangle (gp_Pnt2d & par0, gp_Pnt2d & par1, gp_Pnt2d & par2,
-   BRepLProp_SLProps * prop, Mesh & mesh, const double maxside, int depth, double h = 0)
-   {
-
-
-   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)
-   {
-   double curvature = 0;
-
-   prop->SetParameters (parmid.X(), parmid.Y());
-   if (!prop->IsCurvatureDefined())
-   {
-   (*testout) << "curvature not defined!" << endl;
-   return;
-   }
-   curvature = max(fabs(prop->MinCurvature()),
-   fabs(prop->MaxCurvature()));
-
-   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())));
-
-   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())));
-
-   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())));
-
-   //(*testout) << "curvature " << curvature << endl;
-
-   if (curvature < 1e-3)
-   {
-   //(*testout) << "curvature too small (" << curvature << ")!" << endl;
-   return;
-   // return war bis 10.2.05 auskommentiert
-   }
-
-
-
-   h = ComputeH (curvature+1e-10);
-
-   if(h < 1e-4*maxside)
-   return;
-
-
-   if (h > 30) return;
-   }
-
-   if (h < maxside) // && depth < 5)
-   {
-   //cout << "\r h " << h << flush;
-   gp_Pnt2d pm0;
-   gp_Pnt2d pm1;
-   gp_Pnt2d pm2;
-
-   //cout << "h " << h << " maxside " << maxside << " depth " << depth << endl;
-   //cout << "par0 " << par0.X() << " " << par0.Y()
-   //<< " par1 " << par1.X() << " " << par1.Y()
-   //   << " par2 " << par2.X() << " " << par2.Y()<< endl;
-
-
-   pm0.SetX(0.5*(par1.X()+par2.X())); pm0.SetY(0.5*(par1.Y()+par2.Y()));
-   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);
-   }
-   else
-   {
-   gp_Pnt pnt;
-   Point3d p3d;
-
-   prop->SetParameters (parmid.X(), parmid.Y());
-   pnt = prop->Value();
-   p3d = Point3d(pnt.X(), pnt.Y(), pnt.Z());
-   mesh.RestrictLocalH (p3d, h);
-
-
-   prop->SetParameters (par0.X(), par0.Y());
-   pnt = prop->Value();
-   p3d = Point3d(pnt.X(), pnt.Y(), pnt.Z());
-   mesh.RestrictLocalH (p3d, h);
-
-   prop->SetParameters (par1.X(), par1.Y());
-   pnt = prop->Value();
-   p3d = Point3d(pnt.X(), pnt.Y(), pnt.Z());
-   mesh.RestrictLocalH (p3d, h);
-
-   prop->SetParameters (par2.X(), par2.Y());
-   pnt = prop->Value();
-   p3d = Point3d(pnt.X(), pnt.Y(), pnt.Z());
-   mesh.RestrictLocalH (p3d, h);
-
-   }
-   }
-   */
-
-
-   void RestrictHTriangle (gp_Pnt2d & par0, gp_Pnt2d & par1, gp_Pnt2d & par2,
-      BRepLProp_SLProps * prop, Mesh & mesh, int depth, double h = 0)
-   {
-      int ls = -1;
-
-      gp_Pnt pnt0,pnt1,pnt2;
-
-      prop->SetParameters (par0.X(), par0.Y());
-      pnt0 = prop->Value();
-
-      prop->SetParameters (par1.X(), par1.Y());
-      pnt1 = prop->Value();
-
-      prop->SetParameters (par2.X(), par2.Y());
-      pnt2 = prop->Value();
-
-      double aux;
-      double maxside = pnt0.Distance(pnt1);
-      ls = 2;
-      aux = pnt1.Distance(pnt2);
-      if(aux > maxside)
-      {
-         maxside = aux;
-         ls = 0;
-      }
-      aux = pnt2.Distance(pnt0);
-      if(aux > maxside)
-      {
-         maxside = aux;
-         ls = 1;
-      }
-
-
-
-      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%3 == 0)
-      {
-         double curvature = 0;
-
-         prop->SetParameters (parmid.X(), parmid.Y());
-         if (!prop->IsCurvatureDefined())
-         {
-            (*testout) << "curvature not defined!" << endl;
-            return;
-         }
-         curvature = max(fabs(prop->MinCurvature()),
-            fabs(prop->MaxCurvature()));
-
-         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())));
-
-         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())));
-
-         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())));
-
-         //(*testout) << "curvature " << curvature << endl;
-
-         if (curvature < 1e-3)
-         {
-            //(*testout) << "curvature too small (" << curvature << ")!" << endl;
-            return;
-            // return war bis 10.2.05 auskommentiert
-         }
-
-
-
-         h = ComputeH (curvature+1e-10);
-
-         if(h < 1e-4*maxside)
-            return;
-
-
-         if (h > 30) return;
-      }
-
-      if (h < maxside && depth < 10)
-      {
-         //cout << "\r h " << h << flush;
-         gp_Pnt2d pm;
-
-         //cout << "h " << h << " maxside " << maxside << " depth " << depth << endl;
-         //cout << "par0 " << par0.X() << " " << par0.Y()
-         //<< " par1 " << par1.X() << " " << par1.Y()
-         //   << " par2 " << par2.X() << " " << par2.Y()<< endl;
-
-         if(ls == 0)
-         {
-            pm.SetX(0.5*(par1.X()+par2.X())); pm.SetY(0.5*(par1.Y()+par2.Y()));
-            RestrictHTriangle(pm, par2, par0, prop, mesh, depth+1, h);
-            RestrictHTriangle(pm, par0, par1, prop, mesh, depth+1, h);
-         }
-         else if(ls == 1)
-         {
-            pm.SetX(0.5*(par0.X()+par2.X())); pm.SetY(0.5*(par0.Y()+par2.Y()));
-            RestrictHTriangle(pm, par1, par2, prop, mesh, depth+1, h);
-            RestrictHTriangle(pm, par0, par1, prop, mesh, depth+1, h);
-         }
-         else if(ls == 2)
-         {
-            pm.SetX(0.5*(par0.X()+par1.X())); pm.SetY(0.5*(par0.Y()+par1.Y()));
-            RestrictHTriangle(pm, par1, par2, prop, mesh, depth+1, h);
-            RestrictHTriangle(pm, par2, par0, prop, mesh, depth+1, h);
-         }
-
-      }
-      else
-      {
-         gp_Pnt pnt;
-         Point3d p3d;
-
-         prop->SetParameters (parmid.X(), parmid.Y());
-         pnt = prop->Value();
-         p3d = Point3d(pnt.X(), pnt.Y(), pnt.Z());
-         mesh.RestrictLocalH (p3d, h);
-
-         p3d = Point3d(pnt0.X(), pnt0.Y(), pnt0.Z());
-         mesh.RestrictLocalH (p3d, h);
-
-         p3d = Point3d(pnt1.X(), pnt1.Y(), pnt1.Z());
-         mesh.RestrictLocalH (p3d, h);
-
-         p3d = Point3d(pnt2.X(), pnt2.Y(), pnt2.Z());
-         mesh.RestrictLocalH (p3d, h);
-
-         //(*testout) << "p = " << p3d << ", h = " << h << ", maxside = " << maxside << endl;
-
-      }
-   }
-
-
-
-
-   int OCCGenerateMesh (OCCGeometry & geom,
-      Mesh *& mesh,
-      int perfstepsstart, int perfstepsend,
-      char * optstr)
-   {
-      // int i, j;
-
-      multithread.percent = 0;
-
-      if (perfstepsstart <= MESHCONST_ANALYSE)
-      {
-         delete mesh;
-         mesh = new Mesh();
-         mesh->geomtype = Mesh::GEOM_OCC;
-         mesh->SetGlobalH (mparam.maxh);
-         mesh->SetMinimalH (mparam.minh);
 
       Array<double> maxhdom;
       maxhdom.SetSize (geom.NrSolids());
       maxhdom = mparam.maxh;
 
-         mesh->SetMaxHDomain (maxhdom);
+      mesh.SetMaxHDomain (maxhdom);
 
       Box<3> bb = geom.GetBoundingBox();
       bb.Increase (bb.Diam()/10);
 
-         mesh->SetLocalH (bb.PMin(), bb.PMax(), 0.5);
+      mesh.SetLocalH (bb.PMin(), bb.PMax(), 0.5);
 
       if (mparam.uselocalh)
       {
-
          const char * savetask = multithread.task;
          multithread.percent = 0;
 
-            mesh->SetLocalH (bb.PMin(), bb.PMax(), mparam.grading);
+         mesh.SetLocalH (bb.PMin(), bb.PMax(), mparam.grading);
 
          int nedges = geom.emap.Extent();
 
@@ -1202,7 +1060,7 @@ namespace netgen
             for (int j = 0; j <= maxj; j++)
             {
                gp_Pnt pnt = c->Value (s0+double(j)/maxj*(s1-s0));
-                  mesh->RestrictLocalH (Point3d(pnt.X(), pnt.Y(), pnt.Z()), localh);
+               mesh.RestrictLocalH (Point3d(pnt.X(), pnt.Y(), pnt.Z()), localh);
             }
          }
 
@@ -1235,8 +1093,7 @@ namespace netgen
 
                gp_Pnt pnt = c->Value (s);
 
-                  mesh->RestrictLocalH (Point3d(pnt.X(), pnt.Y(), pnt.Z()),
+               mesh.RestrictLocalH (Point3d(pnt.X(), pnt.Y(), pnt.Z()), ComputeH (fabs(curvature)));
-                     ComputeH (fabs(curvature)));
             }
             // (*testout) << "edge " << i << " max. curvature: " << maxcur << endl;
          }
@@ -1254,6 +1111,7 @@ namespace netgen
             TopLoc_Location loc;
             Handle(Geom_Surface) surf = BRep_Tool::Surface (face);
             Handle(Poly_Triangulation) triangulation = BRep_Tool::Triangulation (face, loc);
+
             if (triangulation.IsNull()) continue;
 
             BRepAdaptor_Surface sf(face, Standard_True);
@@ -1278,7 +1136,7 @@ namespace netgen
                //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, 0);
+               RestrictHTriangle (par[0], par[1], par[2], &prop, mesh, 0);
                //cout << "\rFace " << i << " pos12 ntriangles " << ntriangles << flush;
             }
          }
@@ -1348,8 +1206,8 @@ namespace netgen
                Box3d box;
                box.SetPoint (Point3d(line.p0));
                box.AddPoint (Point3d(line.p1));
-                  double maxhline = max (mesh->GetH(box.PMin()),
+               double maxhline = max (mesh.GetH(box.PMin()),
-                     mesh->GetH(box.PMax()));
+                  mesh.GetH(box.PMax()));
                box.Increase(maxhline);
 
                double mindist = 1e99;
@@ -1377,7 +1235,7 @@ namespace netgen
                   mindist = 1e-3;
                }
 
-                  mesh->RestrictLocalHLine(line.p0, line.p1, mindist);
+               mesh.RestrictLocalHLine(line.p0, line.p1, mindist);
             }
          }
 
@@ -1396,7 +1254,24 @@ namespace netgen
       // ** To use the mesh size file as the sole source 
       // for defining the mesh size, uncheck the "uselocalh"
       // option.
-         mesh->LoadLocalMeshSize (mparam.meshsizefilename);
+      mesh.LoadLocalMeshSize (mparam.meshsizefilename);
+   }
+
+
+
+   int OCCGenerateMesh (OCCGeometry & geom, Mesh *& mesh,
+                        int perfstepsstart, int perfstepsend,
+                        char * optstr)
+   {
+      multithread.percent = 0;
+
+      if (perfstepsstart <= MESHCONST_ANALYSE)
+      {
+         delete mesh;
+         mesh = new Mesh();
+         mesh->geomtype = Mesh::GEOM_OCC;
+
+         OCCSetLocalMeshSize(geom,*mesh);
       }
 
       if (multithread.terminate || perfstepsend <= MESHCONST_ANALYSE)
@@ -1404,7 +1279,7 @@ namespace netgen
 
       if (perfstepsstart <= MESHCONST_MESHEDGES)
       {
-         FindEdges (geom, *mesh);
+         OCCFindEdges (geom, *mesh);
 
          /*
          cout << "Removing redundant points" << endl;
@@ -1503,8 +1378,7 @@ namespace netgen
       {
          multithread.task = "Volume meshing";
 
-         MESHING3_RESULT res =
+         MESHING3_RESULT res = MeshVolume (mparam, *mesh);
-            MeshVolume (mparam, *mesh);
 
          ofstream problemfile("occmesh.rep",ios_base::app);
 

libsrc/occ/occgeom.cpp

@@ -1603,8 +1603,6 @@ namespace netgen
     return false;
   }
 
-  extern int OCCGenerateMesh (OCCGeometry & occgeometry, Mesh*& mesh,
-			      int perfstepsstart, int perfstepsend, char* optstring);
 
   int OCCGeometry :: GenerateMesh (Mesh*& mesh,
 				   int perfstepsstart, int perfstepsend, char* optstring)

libsrc/occ/occgeom.hpp

@@ -80,6 +80,7 @@
 #include "ShapeAnalysis.hxx"
 #include "ShapeBuild_ReShape.hxx"
 
+
 // Philippose - 29/01/2009
 // OpenCascade XDE Support
 // Include support for OpenCascade XDE Features
@@ -164,6 +165,17 @@ namespace netgen
 
 
 
+   class Line
+   {
+   public:
+      Point<3> p0, p1;
+
+      double Dist (Line l);
+
+      double Length ();
+   };
+
+
 
    inline double Det3 (double a00, double a01, double a02,
       double a10, double a11, double a12,
@@ -386,6 +398,20 @@ namespace netgen
    OCCGeometry * LoadOCC_BREP (const char * filename);
 
    extern OCCParameters occparam;
+
+
+   // Philippose - 31.09.2009
+   // External access to the mesh generation functions within the OCC
+   // subsystem (Not sure if this is the best way to implement this....!!)
+   extern int OCCGenerateMesh (OCCGeometry & occgeometry, Mesh*& mesh,
+	                            int perfstepsstart, int perfstepsend, 
+                               char* optstring);
+
+   extern void OCCSetLocalMeshSize(OCCGeometry & geom, Mesh & mesh);
+
+   extern void OCCMeshSurface (OCCGeometry & geom, Mesh & mesh, int perfstepsend);
+
+   extern void OCCFindEdges (OCCGeometry & geom, Mesh & mesh);
 }
 
 #endif

nglib/nglib.cpp

@@ -270,8 +270,10 @@ namespace nglib
    {
       Mesh * m = (Mesh*)mesh;
 
-
+      // Philippose - 30/08/2009
-      MeshingParameters mparam;
+      // Do not locally re-define "mparam" here... "mparam" is a global 
+      // object 
+      //MeshingParameters mparam;
       mparam.maxh = mp->maxh;
       mparam.meshsizefilename = mp->meshsize_filename;
 
@@ -369,6 +371,9 @@ namespace nglib
       return (Ng_Geometry_2D *)geom;
    }
 
+
+
+   
    DLL_HEADER Ng_Result Ng_GenerateMesh_2D (Ng_Geometry_2D * geom,
                                             Ng_Mesh ** mesh,
                                             Ng_Meshing_Parameters * mp)
@@ -388,6 +393,9 @@ namespace nglib
       return NG_OK;
    }
 
+
+
+   
    DLL_HEADER void Ng_HP_Refinement (Ng_Geometry_2D * geom,
       Ng_Mesh * mesh,
       int levels)
@@ -396,6 +404,9 @@ namespace nglib
       HPRefinement (*(Mesh*)mesh, &ref, levels);
    }
 
+
+
+
    DLL_HEADER void Ng_HP_Refinement (Ng_Geometry_2D * geom,
       Ng_Mesh * mesh,
       int levels, double parameter)
@@ -461,12 +472,16 @@ namespace nglib
       return geo2;
    }
 
+
+
    // generate new STL Geometry
    DLL_HEADER Ng_STL_Geometry * Ng_STL_NewGeometry ()
    {
       return (Ng_STL_Geometry*)(void*)new STLGeometry;
    } 
 
+
+
    // after adding triangles (and edges) initialize
    DLL_HEADER Ng_Result Ng_STL_InitSTLGeometry (Ng_STL_Geometry * geom)
    {
@@ -489,6 +504,8 @@ namespace nglib
       return NG_SURFACE_INPUT_ERROR;
    }
 
+
+
    // automatically generates edges:
    DLL_HEADER Ng_Result Ng_STL_MakeEdges (Ng_STL_Geometry * geom,
                                           Ng_Mesh* mesh,
@@ -593,7 +610,8 @@ namespace nglib
    // positive orientation
    // normal vector may be null-pointer
    DLL_HEADER void Ng_STL_AddTriangle (Ng_STL_Geometry * geom, 
-      double * p1, double * p2, double * p3, double * nv)
+                                       double * p1, double * p2, double * p3, 
+                                       double * nv)
    {
       Point<3> apts[3];
       apts[0] = Point<3>(p1[0],p1[1],p1[2]);
@@ -621,7 +639,6 @@ namespace nglib
 
 
 
-
 #ifdef OCCGEOMETRY
    // --------------------- OCC Geometry / Meshing Utility Functions -------------------
 
@@ -632,46 +649,166 @@ namespace nglib
    } 
 
 
+   
    // Loads geometry from STEP File
    DLL_HEADER Ng_OCC_Geometry * Ng_OCC_Load_STEP (const char * filename)
    {
       Ng_OCC_Geometry * geo = Ng_OCC_NewGeometry();
 
-      geo = (Ng_OCC_Geometry *)LoadOCC_STEP(filename);
+      OCCGeometry * occgeo = LoadOCC_STEP(filename);
+
+      // Create the initial triangulation for the OCC 
+      BRepMesh_IncrementalMesh::BRepMesh_IncrementalMesh(occgeo->shape,0.1);
+
+      geo = (Ng_OCC_Geometry *)occgeo;
 
       return (geo);
    }
 
 
+   
    // Loads geometry from IGES File
    DLL_HEADER Ng_OCC_Geometry * Ng_OCC_Load_IGES (const char * filename)
    {
       Ng_OCC_Geometry * geo = Ng_OCC_NewGeometry();
 
-      geo = (Ng_OCC_Geometry *)LoadOCC_IGES(filename);
+      OCCGeometry * occgeo = LoadOCC_IGES(filename);
+
+      // Create the initial triangulation for the OCC 
+      BRepMesh_IncrementalMesh::BRepMesh_IncrementalMesh(occgeo->shape,0.1);
+
+      geo = (Ng_OCC_Geometry *)occgeo;
 
       return (geo);
    }
 
 
+   
    // Loads geometry from BREP File
    DLL_HEADER Ng_OCC_Geometry * Ng_OCC_Load_BREP (const char * filename)
    {
       Ng_OCC_Geometry * geo = Ng_OCC_NewGeometry();
 
-      geo = (Ng_OCC_Geometry *)LoadOCC_BREP(filename);
+      OCCGeometry * occgeo = LoadOCC_BREP(filename);
+
+      // Create the initial triangulation for the OCC 
+      BRepMesh_IncrementalMesh::BRepMesh_IncrementalMesh(occgeo->shape,0.1);
+
+      geo = (Ng_OCC_Geometry *)occgeo;
 
       return (geo);
    }
 
 
-   // Extract a 
+
+   //
+   DLL_HEADER Ng_Result Ng_OCC_SetLocalMeshSize (Ng_OCC_Geometry * geom,
+                                                 Ng_Mesh * mesh,
+                                                 Ng_Meshing_Parameters * mp)
+   {
+      OCCGeometry * occgeom = (OCCGeometry*)geom;
+      Mesh * me = (Mesh*)mesh;
+
+      mparam.uselocalh = mp->uselocalh;
+
+      mparam.maxh = mp->maxh;
+      mparam.minh = mp->minh;
+
+      mparam.segmentsperedge = mp->elementsperedge;
+      mparam.curvaturesafety = mp->elementspercurve;
+
+      mparam.grading = mp->grading;
+      mparam.meshsizefilename = mp->meshsize_filename;
+      
+      occparam.resthcloseedgeenable = mp->closeedgeenable;
+      occparam.resthcloseedgefac = mp->closeedgefact;
+
+      delete me;
+      me = new Mesh;
+
+      OCCSetLocalMeshSize(*occgeom, *me);
+
+      return(NG_OK);
+   }
+
+
+   
+   // Mesh the edges and add Face descriptors to prepare for surface meshing
+   DLL_HEADER Ng_Result Ng_OCC_GenerateEdgeMesh (Ng_OCC_Geometry * geom,
+                                                 Ng_Mesh * mesh,
+                                                 Ng_Meshing_Parameters * mp)
+   {
+      OCCGeometry * occgeom = (OCCGeometry*)geom;
+      Mesh * me = (Mesh*)mesh;
+
+      mparam.uselocalh = mp->uselocalh;
+
+      OCCFindEdges(*occgeom, *me);
+
+      if((me->GetNP()) && (me->GetNFD()))
+      {
+         return NG_OK;
+      }
+      else
+      {
+         return NG_ERROR;
+      }
+   }
+
+
+
+   
+   // Mesh the edges and add Face descriptors to prepare for surface meshing
+   DLL_HEADER Ng_Result Ng_OCC_GenerateSurfaceMesh (Ng_OCC_Geometry * geom,
+                                                    Ng_Mesh * mesh,
+                                                    Ng_Meshing_Parameters * mp)
+   {
+      int numpoints = 0;
+
+      OCCGeometry * occgeom = (OCCGeometry*)geom;
+      Mesh * me = (Mesh*)mesh;
+
+      mparam.uselocalh = mp->uselocalh;
+
+      // Only go into surface meshing if the face descriptors have already been added
+      if(!me->GetNFD())
+         return NG_ERROR;
+
+      numpoints = me->GetNP();
+
+      // Initially set up only for surface meshing without any optimisation
+      int perfstepsend = MESHCONST_MESHSURFACE;
+
+      // Check and if required, enable surface mesh optimisation step
+      if(mp->optsurfmeshenable)
+      {
+         perfstepsend = MESHCONST_OPTSURFACE;
+      }
+
+      OCCMeshSurface(*occgeom, *me, perfstepsend);
+
+      me->CalcSurfacesOfNode();
+      
+      if(me->GetNP() <= numpoints)
+         return NG_ERROR;
+
+      if(me->GetNSE() <= 0)
+         return NG_ERROR;
+
+      return NG_OK;
+   }
+
+
+   // Extract the face map from the OCC geometry
+   // The face map basically gives an index to each face in the geometry, 
+   // which can be used to access a specific face
    DLL_HEADER Ng_Result Ng_OCC_GetFMap(Ng_OCC_Geometry * geom, 
                                        Ng_OCC_TopTools_IndexedMapOfShape * FMap)
    {
       OCCGeometry* occgeom = (OCCGeometry*)geom;
       TopTools_IndexedMapOfShape *occfmap = (TopTools_IndexedMapOfShape *)FMap;
 
+      // Copy the face map from the geometry to the given variable
       occfmap->Assign(occgeom->fmap);
 
       if(occfmap->Extent())
@@ -694,13 +831,28 @@ namespace nglib
 
    DLL_HEADER Ng_Meshing_Parameters :: Ng_Meshing_Parameters()
    {
+      uselocalh = 1;
+
       maxh = 1000;
+      minh = 0.0;
+
       fineness = 0.5;
       grading = 0.3;
+
+      elementsperedge = 2.0;
+      elementspercurve = 2.0;
+
+      closeedgeenable = 0;
+      closeedgefact = 2.0;
+
       secondorder = 0;
-      meshsize_filename = 0;
       quad_dominated = 0;
 
+      meshsize_filename = 0;
+
+      optsurfmeshenable = 1;
+      optvolmeshenable = 1;
+
       optsteps_2d = 3;
       optsteps_3d = 3;
    }
@@ -709,13 +861,28 @@ namespace nglib
 
    DLL_HEADER void Ng_Meshing_Parameters :: Reset_Parameters()
    {
+      uselocalh = 1;
+
       maxh = 1000;
+      minh = 0;
+
       fineness = 0.5;
       grading = 0.3;
+
+      elementsperedge = 2.0;
+      elementspercurve = 2.0;
+
+      closeedgeenable = 0;
+      closeedgefact = 2.0;
+
       secondorder = 0;
-      meshsize_filename = 0;
       quad_dominated = 0;
 
+      meshsize_filename = 0;
+
+      optsurfmeshenable = 1;
+      optvolmeshenable = 1;
+
       optsteps_2d = 3;
       optsteps_3d = 3;
    }
@@ -723,6 +890,8 @@ namespace nglib
 } // End of namespace nglib
 
 
+
+
 // compatibility functions:
 
 namespace netgen 
@@ -762,3 +931,4 @@ namespace netgen
    }
 
 } // End of namespace netgen
+

nglib/nglib.h

@@ -92,12 +92,28 @@ enum Ng_Result
 class Ng_Meshing_Parameters 
 {
 public:
-   double maxh;                //!< Maximum global mesh size limit 
+   int uselocalh;              //!< Switch to enable / disable usage of local mesh size modifiers
+
+   double maxh;                //!< Maximum global mesh size allowed
+   double minh;                //!< Minimum global mesh size allowed
+
    double fineness;            //!< Mesh density: 0...1 (0 => coarse; 1 => fine)
    double grading;             //!< Mesh grading: 0...1 (0 => uniform mesh; 1 => aggressive local grading)
+
+   double elementsperedge;     //!< Number of elements to generate per edge of the geometry
+   double elementspercurve;    //!< Elements to generate per curvature radius
+
+   int closeedgeenable;        //!< Enable / Disable mesh refinement at close edges
+   double closeedgefact;       //!< Factor to use for refinement at close edges (STL: larger => finer ; OCC: larger => coarser)
+
    int secondorder;            //!< Generate second-order surface and volume elements
-   char * meshsize_filename;   //!< Optional external mesh size file 
    int quad_dominated;         //!< Creates a Quad-dominated mesh 
+
+   char * meshsize_filename;   //!< Optional external mesh size file 
+
+   int optsurfmeshenable;      //!< Enable / Disable automatic surface mesh optimization
+   int optvolmeshenable;       //!< Enable / Disable automatic volume mesh optimization
+
    int optsteps_3d;            //!< Number of optimize steps to use for 3-D mesh optimization
    int optsteps_2d;            //!< Number of optimize steps to use for 2-D mesh optimization
 
@@ -107,12 +123,19 @@ public:
 
       Note: This constructor initialises the variables in the 
       class with the following default values
+      - #uselocalh: 1
       - #maxh: 1000.0
       - #fineness: 0.5
       - #grading: 0.3
+      - #elementsperedge: 2.0
+      - #elementspercurve: 2.0
+      - #closeedgeenable: 0
+      - #closeedgefact: 2.0
       - #secondorder: 0.0
       - #meshsize_filename: null
       - #quad_dominated: 0
+      - #optsurfmeshenable: 1
+      - #optvolmeshenable: 1
       - #optsteps_2d: 3
       - #optsteps_3d: 3
    */
@@ -587,10 +610,25 @@ DLL_HEADER Ng_OCC_Geometry * Ng_OCC_Load_IGES (const char * filename);
 // Loads geometry from BREP file
 DLL_HEADER Ng_OCC_Geometry * Ng_OCC_Load_BREP (const char * filename);
 
+// Set the local mesh size based on geometry / topology
+DLL_HEADER Ng_Result Ng_OCC_SetLocalMeshSize (Ng_OCC_Geometry * geom,
+                                              Ng_Mesh * mesh,
+                                              Ng_Meshing_Parameters * mp);
+
+// Mesh the edges and add Face descriptors to prepare for surface meshing
+DLL_HEADER Ng_Result Ng_OCC_GenerateEdgeMesh (Ng_OCC_Geometry * geom,
+                                              Ng_Mesh * mesh,
+                                              Ng_Meshing_Parameters * mp);
+
+// Mesh the surfaces of an OCC geometry
+DLL_HEADER Ng_Result Ng_OCC_GenerateSurfaceMesh (Ng_OCC_Geometry * geom,
+                                                 Ng_Mesh * mesh,
+                                                 Ng_Meshing_Parameters * mp); 
+
 // Get the face map of an already loaded OCC geometry
 DLL_HEADER Ng_Result Ng_OCC_GetFMap(Ng_OCC_Geometry * geom, 
                                     Ng_OCC_TopTools_IndexedMapOfShape * FMap);
 
-#endif
+#endif // OCCGEOMETRY
 
-#endif
+#endif // NGLIB