From 641d40639d6bd76061e9862c68d6919ce22cec19 Mon Sep 17 00:00:00 2001 From: Joachim Schoeberl Date: Tue, 2 Apr 2013 20:27:35 +0000 Subject: [PATCH] Ngx_Mesh C++ interface --- libsrc/interface/nginterface.cpp | 12 +- libsrc/interface/nginterface_v2.cpp | 474 ++++++++++++++++++++++++---- libsrc/interface/writejcm.cpp | 2 +- libsrc/interface/writetet.cpp | 4 +- 4 files changed, 419 insertions(+), 73 deletions(-) diff --git a/libsrc/interface/nginterface.cpp b/libsrc/interface/nginterface.cpp index c6cd4e9b..2f4d3ae1 100644 --- a/libsrc/interface/nginterface.cpp +++ b/libsrc/interface/nginterface.cpp @@ -1785,7 +1785,7 @@ void Ng_GetPeriodicVertices (int idnr, int * pairs) int Ng_GetNPeriodicEdges (int idnr) { - Array map; + Array 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 other2 = map[(*mesh)[si][1]]; + PointIndex other1 = PointIndex (map[(*mesh)[si][0]]); + 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 map; + Array 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 other2 = map[(*mesh)[si][1]]; + PointIndex other1 = PointIndex (map[(*mesh)[si][0]]); + PointIndex other2 = PointIndex (map[(*mesh)[si][1]]); if (other1 && other2 && mesh->IsSegment (other1, other2)) { SegmentIndex otherseg = mesh->SegmentNr (other1, other2); diff --git a/libsrc/interface/nginterface_v2.cpp b/libsrc/interface/nginterface_v2.cpp index a347f94d..6d4caf3f 100644 --- a/libsrc/interface/nginterface_v2.cpp +++ b/libsrc/interface/nginterface_v2.cpp @@ -1,18 +1,9 @@ -#include #include - - #ifdef SOCKETS #include "../sockets/sockets.hpp" #endif -/* -#ifndef NOTCL -#include -#endif -*/ - #include "nginterface.h" #include "nginterface_v2.hpp" @@ -21,15 +12,214 @@ namespace netgen { #include "writeuser.hpp" - extern AutoPtr 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; @@ -103,7 +291,7 @@ namespace netgen return ret; } - template <> DLL_HEADER Ng_Element Ng_GetElement<3> (int nr) + template <> DLL_HEADER Ng_Element Ng_GetElement<3> (int nr) { const Element & el = mesh->VolumeElement (ElementIndex (nr)); @@ -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; - ret.pt = &mesh->Point(nr + PointIndex::BASE)(0); - return ret; + Point<3> xl(xi[0], xi[1], xi[2]); + Point<3> xg; + 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 <> - DLL_HEADER int Ng_GetElementIndex<1> (int nr) + template <> DLL_HEADER void Ngx_Mesh :: + 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; + + mesh->GetCurvedElements().CalcSurfaceTransformation (xl, 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[2*i] = dx(i,0); + dxdxi[2*i+1] = dx(i,1); + } } - - template <> - DLL_HEADER int Ng_GetElementIndex<2> (int nr) + + + + + template <> DLL_HEADER void Ngx_Mesh :: + ElementTransformation<1,2> (int elnr, + const double * xi, + double * x, + double * dxdxi) const { - int ind = (*mesh)[SurfaceElementIndex(nr)].GetIndex(); - return mesh->GetFaceDescriptor(ind).BCProperty(); + 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 int Ng_GetElementIndex<3> (int nr) + + + + template <> DLL_HEADER void Ngx_Mesh :: + ElementTransformation<1,1> (int elnr, + const double * xi, + double * x, + double * dxdxi) const { - return (*mesh)[ElementIndex(nr)].GetIndex(); + 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 Ng_MultiElementTransformation<3,3> (int elnr, int npts, - const double * xi, size_t sxi, - double * x, size_t sx, - double * dxdxi, size_t sdxdxi) + 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) const { mesh->GetCurvedElements().CalcMultiPointElementTransformation (elnr, npts, xi, sxi, x, sx, dxdxi, sdxdxi); } - template <> - DLL_HEADER void Ng_MultiElementTransformation<2,2> (int elnr, int npts, - const double * xi, size_t sxi, - double * x, size_t sx, - double * dxdxi, size_t sdxdxi) + template <> DLL_HEADER void Ngx_Mesh :: + MultiElementTransformation<2,2> (int elnr, int npts, + const double * xi, size_t sxi, + double * x, size_t sx, + double * dxdxi, size_t sdxdxi) const { mesh->GetCurvedElements().CalcMultiPointSurfaceTransformation<2> (elnr, npts, xi, sxi, x, sx, dxdxi, sdxdxi); } - template <> - DLL_HEADER void Ng_MultiElementTransformation<2,3> (int elnr, int npts, - const double * xi, size_t sxi, - double * x, size_t sx, - double * dxdxi, size_t sdxdxi) + template <> DLL_HEADER void Ngx_Mesh :: + MultiElementTransformation<2,3> (int elnr, int npts, + const double * xi, size_t sxi, + double * x, size_t sx, + double * dxdxi, size_t sdxdxi) const { mesh->GetCurvedElements().CalcMultiPointSurfaceTransformation<3> (elnr, npts, xi, sxi, x, sx, dxdxi, sdxdxi); } - template <> - DLL_HEADER void Ng_MultiElementTransformation<1,2> (int elnr, int npts, - const double * xi, size_t sxi, - double * x, size_t sx, - double * dxdxi, size_t sdxdxi) + template <> DLL_HEADER void Ngx_Mesh :: + MultiElementTransformation<1,2> (int elnr, int npts, + const double * xi, size_t sxi, + double * x, size_t sx, + double * dxdxi, size_t sdxdxi) const { mesh->GetCurvedElements().CalcMultiPointSegmentTransformation<2> (elnr, npts, xi, sxi, x, sx, dxdxi, sdxdxi); } - template <> - DLL_HEADER void Ng_MultiElementTransformation<1,1> (int elnr, int npts, - const double * xi, size_t sxi, - double * x, size_t sx, - double * dxdxi, size_t sdxdxi) + template <> DLL_HEADER void Ngx_Mesh :: + MultiElementTransformation<1,1> (int elnr, int npts, + const double * xi, size_t sxi, + double * x, size_t sx, + double * dxdxi, size_t sdxdxi) const { cout << "1D not supported" << endl; } - template <> - DLL_HEADER void Ng_MultiElementTransformation<0,1> (int elnr, int npts, - const double * xi, size_t sxi, - double * x, size_t sx, - double * dxdxi, size_t sdxdxi) + + 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) 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 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 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; + } + + + } diff --git a/libsrc/interface/writejcm.cpp b/libsrc/interface/writejcm.cpp index fc00955d..7517b179 100644 --- a/libsrc/interface/writejcm.cpp +++ b/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); } } diff --git a/libsrc/interface/writetet.cpp b/libsrc/interface/writetet.cpp index 221f0301..bcab41c2 100644 --- a/libsrc/interface/writetet.cpp +++ b/libsrc/interface/writetet.cpp @@ -494,7 +494,7 @@ namespace netgen } - for(PointIndex i=PointIndex::BASE; iSetSize(0); - for(PointIndex i=PointIndex::BASE; i