Ngx_Mesh C++ interface

libsrc/interface/nginterface.cpp

@@ -1785,7 +1785,7 @@ void Ng_GetPeriodicVertices (int idnr, int * pairs)
 
 int Ng_GetNPeriodicEdges (int idnr)
 {
-  Array<INDEX,PointIndex::BASE> map;
+  Array<int,PointIndex::BASE> map;
   //const MeshTopology & top = mesh->GetTopology();
   int nse = mesh->GetNSeg();
 
@@ -1797,8 +1797,8 @@ int Ng_GetNPeriodicEdges (int idnr)
 
     for (SegmentIndex si = 0; si < nse; si++)
       {
-	PointIndex other1 = map[(*mesh)[si][0]];
+	PointIndex other1 = PointIndex (map[(*mesh)[si][0]]);
-	PointIndex other2 = map[(*mesh)[si][1]];
+	PointIndex other2 = PointIndex (map[(*mesh)[si][1]]);
 	//  (*testout) << "seg = " << (*mesh)[si] << "; other = " 
 	//     << other1 << "-" << other2 << endl;
 	if (other1 && other2 && mesh->IsSegment (other1, other2))
@@ -1812,7 +1812,7 @@ int Ng_GetNPeriodicEdges (int idnr)
 
 void Ng_GetPeriodicEdges (int idnr, int * pairs)
 {
-  Array<INDEX,PointIndex::BASE> map;
+  Array<int,PointIndex::BASE> map;
   const MeshTopology & top = mesh->GetTopology();
   int nse = mesh->GetNSeg();
 
@@ -1825,8 +1825,8 @@ void Ng_GetPeriodicEdges (int idnr, int * pairs)
 
     for (SegmentIndex si = 0; si < nse; si++)
       {
-	PointIndex other1 = map[(*mesh)[si][0]];
+	PointIndex other1 = PointIndex (map[(*mesh)[si][0]]);
-	PointIndex other2 = map[(*mesh)[si][1]];
+	PointIndex other2 = PointIndex (map[(*mesh)[si][1]]);
 	if (other1 && other2 && mesh->IsSegment (other1, other2))
 	  {
 	    SegmentIndex otherseg = mesh->SegmentNr (other1, other2);

libsrc/interface/nginterface_v2.cpp

@@ -1,18 +1,9 @@
-#include <mystdlib.h>
 #include <meshing.hpp>
 
-
-
 #ifdef SOCKETS
 #include "../sockets/sockets.hpp"
 #endif
 
-/*
-#ifndef NOTCL
-#include <visual.hpp>
-#endif
-*/
-
 #include "nginterface.h"
 #include "nginterface_v2.hpp"
 
@@ -21,15 +12,214 @@
 namespace netgen
 {
 #include "writeuser.hpp"
-
   extern AutoPtr<Mesh> mesh;
 }
 
 
-
 namespace netgen
 {
 
+  Ngx_Mesh * LoadMesh (const string & filename)
+  {
+    netgen::mesh.Ptr() = NULL;
+    Ng_LoadMesh (filename.c_str());
+    return new Ngx_Mesh (netgen::mesh.Ptr());
+  }
+
+  Ngx_Mesh :: Ngx_Mesh (Mesh * amesh)
+    : mesh(amesh)
+  { ; }
+
+  Ngx_Mesh :: ~Ngx_Mesh ()
+  {
+    if (netgen::mesh.Ptr() == mesh) 
+      netgen::mesh.Ptr() = NULL;
+    delete mesh;
+  }
+
+  int Ngx_Mesh :: GetDimension() const
+  {
+    return mesh -> GetDimension();
+  }
+
+  int Ngx_Mesh :: GetNLevels() const
+  {
+    return mesh -> mglevels;
+  }
+  
+  int Ngx_Mesh :: GetNElements (int dim) const
+  {
+    switch (dim)
+      {
+      case 0: return mesh -> GetNV();
+      case 1: return mesh -> GetNSeg();
+      case 2: return mesh -> GetNSE();
+      case 3: return mesh -> GetNE();
+      }
+    return -1;
+  }
+  
+  int Ngx_Mesh :: GetNNodes (int nt) const
+  {
+    switch (nt)
+      {
+      case 0: return mesh -> GetNV();
+      case 1: return mesh->GetTopology().GetNEdges();
+      case 2: return mesh->GetTopology().GetNFaces();
+      case 3: return mesh -> GetNE();
+      }
+    return -1;
+  }
+
+
+  Ng_Point Ngx_Mesh :: GetPoint (int nr) const
+  {
+    Ng_Point ret;
+    ret.pt = &mesh->Point(nr + PointIndex::BASE)(0);
+    return ret;
+  }
+
+  template <> DLL_HEADER Ng_Element Ngx_Mesh :: GetElement<0> (int nr) const
+  {
+    cout << "Netgen does not support 0-D elements" << endl;
+	Ng_Element ret;
+	return ret;
+  }
+
+  template <> DLL_HEADER Ng_Element Ngx_Mesh :: GetElement<1> (int nr) const
+  {
+    const Segment & el = mesh->LineSegment (SegmentIndex(nr));
+
+    Ng_Element ret;
+    ret.type = NG_ELEMENT_TYPE(el.GetType());
+
+    ret.points.num = el.GetNP();
+    ret.points.ptr = (int*)&(el[0]);
+
+    ret.vertices.num = 2;
+    ret.vertices.ptr = (int*)&(el[0]);
+
+    ret.edges.num = 1;
+    ret.edges.ptr = mesh->GetTopology().GetSegmentElementEdgesPtr (nr);
+
+    ret.faces.num = 0;
+    ret.faces.ptr = NULL;
+
+    return ret;
+  }
+
+  template <> DLL_HEADER Ng_Element Ngx_Mesh :: GetElement<2> (int nr) const
+  {
+    const Element2d & el = mesh->SurfaceElement (SurfaceElementIndex (nr));
+  
+    Ng_Element ret;
+    ret.type = NG_ELEMENT_TYPE(el.GetType());
+    ret.points.num = el.GetNP();
+    ret.points.ptr  = (int*)&el[0];
+
+    ret.vertices.num = el.GetNV();
+    ret.vertices.ptr = (int*)&(el[0]);
+
+    ret.edges.num = MeshTopology::GetNEdges (el.GetType());
+    ret.edges.ptr = mesh->GetTopology().GetSurfaceElementEdgesPtr (nr);
+
+    ret.faces.num = MeshTopology::GetNFaces (el.GetType());
+    ret.faces.ptr = mesh->GetTopology().GetSurfaceElementFacesPtr (nr);
+
+    return ret;
+  }
+
+  template <> DLL_HEADER Ng_Element Ngx_Mesh :: GetElement<3> (int nr) const
+  {
+    const Element & el = mesh->VolumeElement (ElementIndex (nr));
+  
+    Ng_Element ret;
+    ret.type = NG_ELEMENT_TYPE(el.GetType());
+    ret.points.num = el.GetNP();
+    ret.points.ptr = (int*)&el[0];
+
+    ret.vertices.num = el.GetNV();
+    ret.vertices.ptr = (int*)&(el[0]);
+
+    ret.edges.num = MeshTopology::GetNEdges (el.GetType());
+    ret.edges.ptr = mesh->GetTopology().GetElementEdgesPtr (nr);
+
+    ret.faces.num = MeshTopology::GetNFaces (el.GetType());
+    ret.faces.ptr = mesh->GetTopology().GetElementFacesPtr (nr);
+
+    return ret;
+  }
+
+
+  template <>
+  DLL_HEADER int Ngx_Mesh :: GetElementIndex<0> (int nr) const
+  {
+    return 0;
+  }
+
+  template <>
+  DLL_HEADER int Ngx_Mesh :: GetElementIndex<1> (int nr) const
+  {
+    return (*mesh)[SegmentIndex(nr)].si;
+  }
+  
+  template <>
+  DLL_HEADER int Ngx_Mesh :: GetElementIndex<2> (int nr) const
+  {
+    int ind = (*mesh)[SurfaceElementIndex(nr)].GetIndex(); 
+    return mesh->GetFaceDescriptor(ind).BCProperty();
+  }
+  
+  template <>
+  DLL_HEADER int Ngx_Mesh :: GetElementIndex<3> (int nr) const
+  {
+    return (*mesh)[ElementIndex(nr)].GetIndex();
+  }
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+  /*
+  DLL_HEADER Ng_Point Ng_GetPoint (int nr)
+  {
+    Ng_Point ret;
+    ret.pt = &mesh->Point(nr + PointIndex::BASE)(0);
+    return ret;
+  }
+
+
+  template <>
+  DLL_HEADER int Ng_GetElementIndex<1> (int nr)
+  {
+    return (*mesh)[SegmentIndex(nr)].si;
+  }
+  
+  template <>
+  DLL_HEADER int Ng_GetElementIndex<2> (int nr)
+  {
+    int ind = (*mesh)[SurfaceElementIndex(nr)].GetIndex(); 
+    return mesh->GetFaceDescriptor(ind).BCProperty();
+  }
+  
+  template <>
+  DLL_HEADER int Ng_GetElementIndex<3> (int nr)
+  {
+    return (*mesh)[ElementIndex(nr)].GetIndex();
+  }
 
   template <> int DLL_HEADER Ng_GetNElements<0> ()
   {
@@ -51,8 +241,6 @@ namespace netgen
     return mesh->GetNE();
   }
 
-
-
   template <> DLL_HEADER Ng_Element Ng_GetElement<0> (int nr)
   {
     cout << "Netgen does not support 0-D elements" << endl;
@@ -123,110 +311,217 @@ namespace netgen
 
     return ret;
   }
+  */
 
 
-  DLL_HEADER Ng_Point Ng_GetPoint (int nr)
+
+
+
+
+
+  
+  
+  template <> DLL_HEADER void Ngx_Mesh :: 
+  ElementTransformation<3,3> (int elnr, 
+                              const double * xi,
+                              double * x, 
+                              double * dxdxi) const
   {
-    Ng_Point ret;
+    Point<3> xl(xi[0], xi[1], xi[2]);
-    ret.pt = &mesh->Point(nr + PointIndex::BASE)(0);
+    Point<3> xg;
-    return ret;
+    Mat<3,3> dx;
+    mesh->GetCurvedElements().CalcElementTransformation (xl, elnr, xg, dx);
+
+    if (x)
+      for (int i = 0; i < 3; i++) x[i] = xg(i);
+    
+    if (dxdxi)
+      for (int i=0; i<3; i++)
+        {
+          dxdxi[3*i] = dx(i,0);
+          dxdxi[3*i+1] = dx(i,1);
+          dxdxi[3*i+2] = dx(i,2);
+        }
+  }
+
+  template <> DLL_HEADER void Ngx_Mesh ::
+  ElementTransformation<2,3> (int elnr,
+                              const double * xi, 
+                              double * x, 
+                              double * dxdxi) const
+  {
+    Point<2> xl(xi[0], xi[1]);
+    Point<3> xg;
+    Mat<3,2> dx;
+
+    mesh->GetCurvedElements().CalcSurfaceTransformation (xl, elnr, xg, dx);
+    
+    if (x)
+      for (int i = 0; i < 3; i++) x[i] = xg(i);
+    
+    if (dxdxi)
+      for (int i=0; i<3; i++)
+        {
+          dxdxi[2*i] = dx(i,0);
+          dxdxi[2*i+1] = dx(i,1);
+        }
   }
 
 
-  template <>
+  template <> DLL_HEADER void Ngx_Mesh ::
-  DLL_HEADER int Ng_GetElementIndex<1> (int nr)
+  ElementTransformation<2,2> (int elnr,
+                                 const double * xi, 
+                                 double * x, 
+                                 double * dxdxi) const
   {
-    return (*mesh)[SegmentIndex(nr)].si;
+    Point<2> xl(xi[0], xi[1]);
-  }
+    Point<3> xg;
+    Mat<3,2> dx;
 
-  template <>
+    mesh->GetCurvedElements().CalcSurfaceTransformation (xl, elnr, xg, dx);
-  DLL_HEADER int Ng_GetElementIndex<2> (int nr)
-  {
-    int ind = (*mesh)[SurfaceElementIndex(nr)].GetIndex(); 
-    return mesh->GetFaceDescriptor(ind).BCProperty();
-  }
     
-  template <>
+    if (x)
-  DLL_HEADER int Ng_GetElementIndex<3> (int nr)
+      for (int i = 0; i < 2; i++) x[i] = xg(i);
+    
+    if (dxdxi)
+      for (int i=0; i<2; i++)
         {
-    return (*mesh)[ElementIndex(nr)].GetIndex();
+          dxdxi[2*i] = dx(i,0);
+          dxdxi[2*i+1] = dx(i,1);
+        }
   }
 
 
-  template <>
+
-  DLL_HEADER void Ng_MultiElementTransformation<3,3> (int elnr, int npts,
+
+  template <> DLL_HEADER void Ngx_Mesh :: 
+  ElementTransformation<1,2> (int elnr,
+                              const double * xi,
+                              double * x,
+                              double * dxdxi) const
+  {
+    Point<3> xg;
+    Vec<3> dx;
+
+    mesh->GetCurvedElements().CalcSegmentTransformation (xi[0], elnr, xg, dx);
+    
+    if (x)
+      for (int i = 0; i < 2; i++) x[i] = xg(i);
+    
+    if (dxdxi)
+      for (int i=0; i < 2; i++)
+        dxdxi[i] = dx(i);
+  }
+
+
+
+  template <> DLL_HEADER void Ngx_Mesh :: 
+  ElementTransformation<1,1> (int elnr,
+                              const double * xi,
+                              double * x,
+                              double * dxdxi) const
+  {
+    cout << "1D not supported" << endl;
+  }
+
+
+  template <> DLL_HEADER void Ngx_Mesh :: 
+  ElementTransformation<0,1> (int elnr, 
+                              const double * xi,
+                              double * x,
+                              double * dxdxi) const
+  {
+    cout << "1D not supported" << endl;
+  }
+
+
+
+
+
+
+
+
+
+  
+  
+  template <> DLL_HEADER void Ngx_Mesh :: 
+  MultiElementTransformation<3,3> (int elnr, int npts,
                                    const double * xi, size_t sxi,
                                    double * x, size_t sx,
-                                                      double * dxdxi, size_t sdxdxi)
+                                   double * dxdxi, size_t sdxdxi) const
   {
     mesh->GetCurvedElements().CalcMultiPointElementTransformation (elnr, npts, xi, sxi, x, sx, dxdxi, sdxdxi);
   }
   
-  template <>
+  template <> DLL_HEADER void Ngx_Mesh ::
-  DLL_HEADER void Ng_MultiElementTransformation<2,2> (int elnr, int npts,
+  MultiElementTransformation<2,2> (int elnr, int npts,
                                    const double * xi, size_t sxi,
                                    double * x, size_t sx,
-                                                      double * dxdxi, size_t sdxdxi)
+                                   double * dxdxi, size_t sdxdxi) const
   {
     mesh->GetCurvedElements().CalcMultiPointSurfaceTransformation<2> (elnr, npts, xi, sxi, x, sx, dxdxi, sdxdxi);
   }
 
-  template <>
+  template <> DLL_HEADER void Ngx_Mesh :: 
-  DLL_HEADER void Ng_MultiElementTransformation<2,3> (int elnr, int npts,
+  MultiElementTransformation<2,3> (int elnr, int npts,
                                    const double * xi, size_t sxi,
                                    double * x, size_t sx,
-                                                      double * dxdxi, size_t sdxdxi)
+                                   double * dxdxi, size_t sdxdxi) const
   {
     mesh->GetCurvedElements().CalcMultiPointSurfaceTransformation<3> (elnr, npts, xi, sxi, x, sx, dxdxi, sdxdxi);
   }
 
-  template <>
+  template <> DLL_HEADER void Ngx_Mesh :: 
-  DLL_HEADER void Ng_MultiElementTransformation<1,2> (int elnr, int npts,
+  MultiElementTransformation<1,2> (int elnr, int npts,
                                    const double * xi, size_t sxi,
                                    double * x, size_t sx,
-                                                      double * dxdxi, size_t sdxdxi)
+                                   double * dxdxi, size_t sdxdxi) const
   {
     mesh->GetCurvedElements().CalcMultiPointSegmentTransformation<2> (elnr, npts, xi, sxi, x, sx, dxdxi, sdxdxi);
   }
 
-  template <>
+  template <> DLL_HEADER void Ngx_Mesh :: 
-  DLL_HEADER void Ng_MultiElementTransformation<1,1> (int elnr, int npts,
+  MultiElementTransformation<1,1> (int elnr, int npts,
                                    const double * xi, size_t sxi,
                                    double * x, size_t sx,
-                                                      double * dxdxi, size_t sdxdxi)
+                                   double * dxdxi, size_t sdxdxi) const
   {
     cout << "1D not supported" << endl;
   }
 
-  template <>
+
-  DLL_HEADER void Ng_MultiElementTransformation<0,1> (int elnr, int npts,
+  template <> DLL_HEADER void Ngx_Mesh :: 
+  MultiElementTransformation<0,1> (int elnr, int npts,
                                    const double * xi, size_t sxi,
                                    double * x, size_t sx,
-                                                      double * dxdxi, size_t sdxdxi)
+                                   double * dxdxi, size_t sdxdxi) const
   {
     cout << "1D not supported" << endl;
   }
 
 
 
-  template <> DLL_HEADER int Ng_GetNNodes<1> ()
+
+
+
+  template <> DLL_HEADER int Ngx_Mesh :: GetNNodes<1> ()
   {
     return mesh->GetTopology().GetNEdges();
   }
 
-  template <> DLL_HEADER int Ng_GetNNodes<2> ()
+  template <> DLL_HEADER int Ngx_Mesh :: GetNNodes<2> ()
   {
     return mesh->GetTopology().GetNFaces();
   }
 
-  template <> DLL_HEADER Ng_Node<1> Ng_GetNode<1> (int nr)
+  template <> DLL_HEADER Ng_Node<1> Ngx_Mesh :: GetNode<1> (int nr)
   {
     Ng_Node<1> node;
     node.vertices.ptr = mesh->GetTopology().GetEdgeVerticesPtr(nr);
     return node;
   }
 
-  template <> DLL_HEADER Ng_Node<2> Ng_GetNode<2> (int nr)
+  template <> DLL_HEADER Ng_Node<2> Ngx_Mesh :: GetNode<2> (int nr)
   {
     Ng_Node<2> node;
     node.vertices.ptr = mesh->GetTopology().GetFaceVerticesPtr(nr);
@@ -234,6 +529,57 @@ namespace netgen
     return node;
   }
 
+
+  template <>
+  DLL_HEADER int Ngx_Mesh :: FindElementOfPoint <2>
+  (double * p, double * lami,
+   bool build_searchtree, 
+   int * const indices, int numind)
+
+  {
+    Array<int> dummy(numind);
+    for (int i = 0; i < numind; i++) dummy[i] = indices[i]+1;
+    
+    double lam3[3];
+    Point<3> p2d(p[0], p[1], 0);
+    int ind = 
+      mesh->GetElementOfPoint(p2d, lam3, &dummy, build_searchtree);
+    
+    if (ind > 0)
+      {
+        if(mesh->SurfaceElement(ind).GetType()==QUAD)
+          {
+            lami[0] = lam3[0];
+            lami[1] = lam3[1];
+          }
+        else 
+          {
+            lami[0] = 1-lam3[0]-lam3[1];
+            lami[1] = lam3[0];
+          }
+      }
+    return ind-1;
+  }
+
+
+  template <>
+  DLL_HEADER int Ngx_Mesh :: FindElementOfPoint <3>
+  (double * p, double * lami,
+   bool build_searchtree, 
+   int * const indices, int numind)
+
+  {
+    Array<int> dummy(numind);
+    for (int i = 0; i < numind; i++) dummy[i] = indices[i]+1;
+    
+    Point<3> p3d(p[0], p[1], p[2]);
+    int ind = 
+      mesh->GetElementOfPoint(p3d, lami, &dummy, build_searchtree);
+    return ind-1;
+  }
+
+
+
 }
 
 

libsrc/interface/writejcm.cpp

@@ -104,7 +104,7 @@ void WriteJCMFormat (const Mesh & mesh,
     if (el.GetNP() == 4)
     {
       for (j = 1; j <= 4; j++)
-        pointsOnTetras.Set(el.PNum(j).GetInt(),1);     
+        pointsOnTetras.Set(int (el.PNum(j)),1);     
     }
   }
 

libsrc/interface/writetet.cpp

@@ -494,7 +494,7 @@ namespace netgen
       }
 
 
-    for(PointIndex i=PointIndex::BASE; i<mesh.GetNP()+PointIndex::BASE; i++)
+    for(PointIndex i = mesh.Points().Begin(); i < mesh.Points().End(); i++)
       {
 	if(nodenum[i] == -1)
 	  continue;
@@ -1063,7 +1063,7 @@ namespace netgen
     
     for(int i=0; i<groups.Size(); i++)
       groups[i]->SetSize(0);
-    for(PointIndex i=PointIndex::BASE; i<mesh.GetNP()+PointIndex::BASE; i++)
+    for(PointIndex i = mesh.Points().Begin(); i < mesh.Points().End(); i++)
       {
 	if(i-PointIndex::BASE < point_ids.Size())
 	  {