diff --git a/libsrc/include/nginterface_v2.hpp b/libsrc/include/nginterface_v2.hpp index f825080f..e53d6dc8 100644 --- a/libsrc/include/nginterface_v2.hpp +++ b/libsrc/include/nginterface_v2.hpp @@ -15,11 +15,16 @@ namespace netgen { - - - - - + struct T_EDGE2 + { + int orient:1; + int nr:31; // 0-based + }; + struct T_FACE2 + { + int orient:3; + int nr:29; // 0-based + }; class Ng_Element { @@ -49,20 +54,22 @@ namespace netgen { public: int num; - const int * ptr; + const T_EDGE2 * ptr; int Size() const { return num; } - int operator[] (int i) const { return abs (ptr[i])-1; } + // int operator[] (int i) const { return abs (ptr[i])-1; } + int operator[] (int i) const { return ptr[i].nr; } }; class Ng_Faces { public: int num; - const int * ptr; + const T_FACE2 * ptr; int Size() const { return num; } - int operator[] (int i) const { return (ptr[i]-1) / 8; } + // int operator[] (int i) const { return (ptr[i]-1) / 8; } + int operator[] (int i) const { return ptr[i].nr; } }; public: @@ -78,10 +85,12 @@ namespace netgen class Ng_Point { - public: double * pt; + public: + Ng_Point (double * apt) : pt(apt) { ; } double operator[] (int i) { return pt[i]; } + operator const double * () { return pt; } }; @@ -152,8 +161,13 @@ namespace netgen class Mesh * mesh; public: - Ngx_Mesh(class Mesh * amesh); + Ngx_Mesh () { ; } + Ngx_Mesh(class Mesh * amesh) : mesh(amesh) { ; } + void LoadMesh (const string & filename); + virtual ~Ngx_Mesh(); + + bool Valid () { return mesh != NULL; } int GetDimension() const; int GetNLevels() const; @@ -197,7 +211,7 @@ namespace netgen template - Ng_Node GetNode (int nr); + Ng_Node GetNode (int nr) const; template @@ -208,16 +222,27 @@ namespace netgen int FindElementOfPoint (double * p, double * lami, bool build_searchtrees = false, - int * const indices = NULL, int numind = 0); + int * const indices = NULL, int numind = 0) const; }; DLL_HEADER Ngx_Mesh * LoadMesh (const string & filename); - - - } + + +#ifdef HAVE_NETGEN_SOURCES +#include + +namespace netgen +{ +#define NGX_INLINE inline +#include +} + +#endif + + #endif diff --git a/libsrc/include/nginterface_v2_impl.hpp b/libsrc/include/nginterface_v2_impl.hpp new file mode 100644 index 00000000..37b16ef9 --- /dev/null +++ b/libsrc/include/nginterface_v2_impl.hpp @@ -0,0 +1,91 @@ +NGX_INLINE DLL_HEADER Ng_Point Ngx_Mesh :: GetPoint (int nr) const +{ + return Ng_Point (&mesh->Point(nr + PointIndex::BASE)(0)); +} + + +template <> +NGX_INLINE DLL_HEADER int Ngx_Mesh :: GetElementIndex<1> (int nr) const +{ + return (*mesh)[SegmentIndex(nr)].si; +} + +template <> +NGX_INLINE DLL_HEADER int Ngx_Mesh :: GetElementIndex<2> (int nr) const +{ + int ind = (*mesh)[SurfaceElementIndex(nr)].GetIndex(); + return mesh->GetFaceDescriptor(ind).BCProperty(); +} + +template <> +NGX_INLINE DLL_HEADER int Ngx_Mesh :: GetElementIndex<3> (int nr) const +{ + return (*mesh)[ElementIndex(nr)].GetIndex(); +} + +template <> +NGX_INLINE 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 = (T_EDGE2*)mesh->GetTopology().GetSegmentElementEdgesPtr (nr); + + ret.faces.num = 0; + ret.faces.ptr = NULL; + + return ret; +} + +template <> +NGX_INLINE 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 = (T_EDGE2*)mesh->GetTopology().GetSurfaceElementEdgesPtr (nr); + + ret.faces.num = MeshTopology::GetNFaces (el.GetType()); + ret.faces.ptr = (T_FACE2*)mesh->GetTopology().GetSurfaceElementFacesPtr (nr); + + return ret; +} + +template <> +NGX_INLINE 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 = (T_EDGE2*)mesh->GetTopology().GetElementEdgesPtr (nr); + + ret.faces.num = MeshTopology::GetNFaces (el.GetType()); + ret.faces.ptr = (T_FACE2*)mesh->GetTopology().GetElementFacesPtr (nr); + + return ret; +} diff --git a/libsrc/interface/nginterface_v2.cpp b/libsrc/interface/nginterface_v2.cpp index e116c548..2f682982 100644 --- a/libsrc/interface/nginterface_v2.cpp +++ b/libsrc/interface/nginterface_v2.cpp @@ -18,6 +18,8 @@ namespace netgen namespace netgen { +#define NGX_INLINE +#include "nginterface_v2_impl.hpp" Ngx_Mesh * LoadMesh (const string & filename) { @@ -26,9 +28,18 @@ namespace netgen return new Ngx_Mesh (netgen::mesh.Ptr()); } + void Ngx_Mesh :: LoadMesh (const string & filename) + { + netgen::mesh.Ptr() = NULL; + Ng_LoadMesh (filename.c_str()); + mesh = netgen::mesh.Ptr(); + } + + /* Ngx_Mesh :: Ngx_Mesh (Mesh * amesh) : mesh(amesh) { ; } + */ Ngx_Mesh :: ~Ngx_Mesh () { @@ -71,13 +82,12 @@ namespace netgen return -1; } - + /* Ng_Point Ngx_Mesh :: GetPoint (int nr) const { - Ng_Point ret; - ret.pt = &mesh->Point(nr + PointIndex::BASE)(0); - return ret; + return Ng_Point (&mesh->Point(nr + PointIndex::BASE)(0)); } + */ template <> DLL_HEADER Ng_Element Ngx_Mesh :: GetElement<0> (int nr) const { @@ -86,6 +96,7 @@ namespace netgen return ret; } + /* template <> DLL_HEADER Ng_Element Ngx_Mesh :: GetElement<1> (int nr) const { const Segment & el = mesh->LineSegment (SegmentIndex(nr)); @@ -149,7 +160,7 @@ namespace netgen return ret; } - + */ template <> DLL_HEADER int Ngx_Mesh :: GetElementIndex<0> (int nr) const @@ -157,6 +168,7 @@ namespace netgen return 0; } + /* template <> DLL_HEADER int Ngx_Mesh :: GetElementIndex<1> (int nr) const { @@ -175,7 +187,7 @@ namespace netgen { return (*mesh)[ElementIndex(nr)].GetIndex(); } - + */ @@ -521,14 +533,14 @@ namespace netgen return mesh->GetTopology().GetNFaces(); } - template <> DLL_HEADER Ng_Node<1> Ngx_Mesh :: GetNode<1> (int nr) + template <> DLL_HEADER Ng_Node<1> Ngx_Mesh :: GetNode<1> (int nr) const { Ng_Node<1> node; node.vertices.ptr = mesh->GetTopology().GetEdgeVerticesPtr(nr); return node; } - template <> DLL_HEADER Ng_Node<2> Ngx_Mesh :: GetNode<2> (int nr) + template <> DLL_HEADER Ng_Node<2> Ngx_Mesh :: GetNode<2> (int nr) const { Ng_Node<2> node; node.vertices.ptr = mesh->GetTopology().GetFaceVerticesPtr(nr); @@ -538,10 +550,10 @@ namespace netgen template <> - DLL_HEADER int Ngx_Mesh :: FindElementOfPoint <2> + DLL_HEADER int Ngx_Mesh :: FindElementOfPoint <2> (double * p, double * lami, bool build_searchtree, - int * const indices, int numind) + int * const indices, int numind) const { Array dummy(numind); @@ -573,7 +585,7 @@ namespace netgen DLL_HEADER int Ngx_Mesh :: FindElementOfPoint <3> (double * p, double * lami, bool build_searchtree, - int * const indices, int numind) + int * const indices, int numind) const { Array dummy(numind); diff --git a/libsrc/meshing/curvedelems.cpp b/libsrc/meshing/curvedelems.cpp index 3a6bbfb9..ad322e1c 100644 --- a/libsrc/meshing/curvedelems.cpp +++ b/libsrc/meshing/curvedelems.cpp @@ -239,6 +239,24 @@ namespace netgen values[i+1] = p1; } } + + template + void EvaluateScaled (int n, S x, S y, T * values) + { + S p1 = 1.0, p2 = 0.0, p3; + + if (n >= 0) + p2 = values[0] = 1.0; + if (n >= 1) + p1 = values[1] = a[0]*y+b[0]*x; + + for (int i = 1; i < n; i++) + { + p3 = p2; p2=p1; + p1 = (a[i]*y+b[i]*x)*p2-c[i]*y*y*p3; + values[i+1] = p1; + } + } }; class JacobiRecPol : public RecPol @@ -331,8 +349,18 @@ namespace netgen if (n < 3) return; Tx hx[50], hy[50*50]; - ScaledJacobiPolynomial (n-3, x, 1-y, 2, 2, hx); - + // ScaledJacobiPolynomial (n-3, x, 1-y, 2, 2, hx); + /* + cout << "scaled jacobi, old: " << endl; + for (int i = 0; i <= n-3; i++) + cout << i << ": " << hx[i] << endl; + */ + jacpols2[2] -> EvaluateScaled (n-3, x, 1-y, hx); + /* + cout << "scaled jacobi, new: " << endl; + for (int i = 0; i <= n-3; i++) + cout << i << ": " << hx[i] << endl; + */ // for (int ix = 0; ix <= n-3; ix++) // JacobiPolynomial (n-3, 2*y-1, 2*ix+5, 2, hy+50*ix); @@ -351,7 +379,6 @@ namespace netgen } - static void CalcTrigShapeDxDy (int n, double x, double y, double * dshape) { if (n < 3) return; @@ -366,35 +393,9 @@ namespace netgen dshape[2*i] = res[i].DValue(0); dshape[2*i+1] = res[i].DValue(1); } - - /* - if (n < 3) return; - - int ndof = (n-1)*(n-2)/2; - double h1[1000], h2[1000]; - double eps = 1e-4; - - CalcTrigShape (n, x+eps, y, h1); - CalcTrigShape (n, x-eps, y, h2); - - for (int i = 0; i < ndof; i++) - dshape[2*i] = (h1[i]-h2[i])/(2*eps); - - CalcTrigShape (n, x, y+eps, h1); - CalcTrigShape (n, x, y-eps, h2); - - for (int i = 0; i < ndof; i++) - dshape[2*i+1] = (h1[i]-h2[i])/(2*eps); - */ } - - - - - - // compute face bubbles up to order n, 0 < y, y-x < 1, x+y < 1 template static void CalcScaledTrigShape (int n, Tx x, Ty y, Tt t, Tr * shape) @@ -407,7 +408,8 @@ namespace netgen // ScaledLegendrePolynomial (n-3, (2*y-1), t, hy); for (int ix = 0; ix <= n-3; ix++) - ScaledJacobiPolynomial (n-3, 2*y-1, t, 2*ix+5, 2, hy+50*ix); + jacpols2[2*ix+5] -> EvaluateScaled (n-3, 2*y-1, t, hy+50*ix); + // ScaledJacobiPolynomial (n-3, 2*y-1, t, 2*ix+5, 2, hy+50*ix); int ii = 0; @@ -434,44 +436,9 @@ namespace netgen dshape[3*i+1] = res[i].DValue(1); dshape[3*i+2] = res[i].DValue(2); } - - /* - double dshape1[6000]; - if (n < 3) return; - double hvl[1000], hvr[1000]; - int nd = (n-1)*(n-2)/2; - - double eps = 1e-6; - - CalcScaledTrigShape (n, x-eps, y, t, hvl); - CalcScaledTrigShape (n, x+eps, y, t, hvr); - for (int i = 0; i < nd; i++) - dshape[3*i] = (hvr[i]-hvl[i])/(2*eps); - - CalcScaledTrigShape (n, x, y-eps, t, hvl); - CalcScaledTrigShape (n, x, y+eps, t, hvr); - for (int i = 0; i < nd; i++) - dshape[3*i+1] = (hvr[i]-hvl[i])/(2*eps); - - CalcScaledTrigShape (n, x, y, t-eps, hvl); - CalcScaledTrigShape (n, x, y, t+eps, hvr); - for (int i = 0; i < nd; i++) - dshape[3*i+2] = (hvr[i]-hvl[i])/(2*eps); - */ - - /* - for (int i = 0; i < 3*nd; i++) - if (fabs (dshape[i]-dshape1[i]) > 1e-8 * fabs(dshape[i]) + 1e-6) - { - cerr - cerr << "big difference: " << dshape1[i] << " != " << dshape[i] << endl; - } - */ } - - - + CurvedElements :: CurvedElements (const Mesh & amesh) : mesh (amesh) @@ -2168,6 +2135,8 @@ namespace netgen bool CurvedElements :: IsElementCurved (ElementIndex elnr) const { + if (mesh[elnr].GetType() != TET) return true; + if (mesh.coarsemesh) { const HPRefElement & hpref_el = @@ -2178,7 +2147,7 @@ namespace netgen const Element & el = mesh[elnr]; ELEMENT_TYPE type = el.GetType(); - + int nfaces = MeshTopology::GetNFaces (type); if (nfaces > 4) { // not a tet @@ -2225,6 +2194,42 @@ namespace netgen } + bool CurvedElements :: IsElementHighOrder (ElementIndex elnr) const + { + if (mesh.coarsemesh) + { + const HPRefElement & hpref_el = + (*mesh.hpelements) [mesh[elnr].hp_elnr]; + + return mesh.coarsemesh->GetCurvedElements().IsElementHighOrder (hpref_el.coarse_elnr); + } + + const Element & el = mesh[elnr]; + ELEMENT_TYPE type = el.GetType(); + + ElementInfo info; + info.elnr = elnr; + info.order = order; + info.ndof = info.nv = MeshTopology::GetNPoints (type); + if (info.order > 1) + { + const MeshTopology & top = mesh.GetTopology(); + + info.nedges = top.GetElementEdges (elnr+1, info.edgenrs, 0); + for (int i = 0; i < info.nedges; i++) info.edgenrs[i]--; + + info.nfaces = top.GetElementFaces (elnr+1, info.facenrs, 0); + for (int i = 0; i < info.nfaces; i++) info.facenrs[i]--; + + for (int i = 0; i < info.nedges; i++) + if (edgecoeffsindex[info.edgenrs[i]+1] > edgecoeffsindex[info.edgenrs[i]]) return true; + for (int i = 0; i < info.nfaces; i++) + if (facecoeffsindex[info.facenrs[i]+1] > facecoeffsindex[info.facenrs[i]]) return true; + } + return false; + } + + @@ -3589,6 +3594,8 @@ namespace netgen double * x, size_t sx, double * dxdxi, size_t sdxdxi) { + // static int timer = NgProfiler::CreateTimer ("calcmultipointtrafo, calcshape"); + if (mesh.coarsemesh) { const HPRefElement & hpref_el = @@ -3675,6 +3682,7 @@ namespace netgen const Element & el = mesh[elnr]; ELEMENT_TYPE type = el.GetType(); + ElementInfo info; info.elnr = elnr; info.order = order; diff --git a/libsrc/meshing/curvedelems.hpp b/libsrc/meshing/curvedelems.hpp index 945f20c8..11d281ef 100644 --- a/libsrc/meshing/curvedelems.hpp +++ b/libsrc/meshing/curvedelems.hpp @@ -52,6 +52,7 @@ public: bool IsSegmentCurved (SegmentIndex segnr) const; bool IsSurfaceElementCurved (SurfaceElementIndex sei) const; bool IsElementCurved (ElementIndex ei) const; + bool IsElementHighOrder (ElementIndex ei) const; void CalcSegmentTransformation (double xi, SegmentIndex segnr, diff --git a/libsrc/meshing/topology.cpp b/libsrc/meshing/topology.cpp index 263fec94..26ef075e 100644 --- a/libsrc/meshing/topology.cpp +++ b/libsrc/meshing/topology.cpp @@ -164,10 +164,10 @@ namespace netgen for (int i = 0; i < ne; i++) for (int j = 0; j < 12; j++) - edges[i][j] = 0; + edges[i][j].nr = -1; for (int i = 0; i < nse; i++) for (int j = 0; j < 4; j++) - surfedges[i][j] = 0; + surfedges[i][j].nr = -1; // keep existing edges cnt = 0; @@ -303,10 +303,14 @@ namespace netgen if (edgedir) swap (edge.I1(), edge.I2()); if (edge.I1() != i) continue; - + int edgenum = edgenr[edge.I2()]; + /* if (edgedir) edgenum *= -1; edges[elnr][k] = edgenum; + */ + edges[elnr][k].nr = edgenum-1; + edges[elnr][k].orient = edgedir; } } @@ -328,8 +332,11 @@ namespace netgen if (edge.I1() != i) continue; int edgenum = edgenr[edge.I2()]; - if (edgedir) edgenum *= -1; - surfedges.Elem(elnr)[k] = edgenum; + // if (edgedir) edgenum *= -1; + // surfedges.Elem(elnr)[k] = edgenum; + surfedges.Elem(elnr)[k].nr = edgenum-1; + surfedges.Elem(elnr)[k].orient = edgedir; + } } @@ -346,8 +353,12 @@ namespace netgen if (edge.I1() != i) continue; int edgenum = edgenr[edge.I2()]; + /* if (edgedir) edgenum *= -1; segedges.Elem(elnr) = edgenum; + */ + segedges.Elem(elnr).nr = edgenum-1; + segedges.Elem(elnr).orient = edgedir; } } } @@ -358,13 +369,16 @@ namespace netgen if (buildfaces) { static int timer2 = NgProfiler::CreateTimer ("topology::buildfaces"); + static int timer2a = NgProfiler::CreateTimer ("topology::buildfacesa"); + static int timer2b = NgProfiler::CreateTimer ("topology::buildfacesb"); + static int timer2c = NgProfiler::CreateTimer ("topology::buildfacesc"); + static int timer2d = NgProfiler::CreateTimer ("topology::buildfacesd"); NgProfiler::RegionTimer reg2 (timer2); if (id == 0) PrintMessage (5, "Update faces "); - - + NgProfiler::StartTimer (timer2a); faces.SetSize(ne); surffaces.SetSize(nse); @@ -381,7 +395,7 @@ namespace netgen for (int elnr = 0; elnr < ne; elnr++) for (int j = 0; j < 6; j++) - faces[elnr][j] = 0; + faces[elnr][j].fnr = -1; int max_face_on_vertex = 0; @@ -392,245 +406,16 @@ namespace netgen } - /* - for (int pass = 1; pass <= 2; pass++) - { - nfa = oldnfa; - for (int v = PointIndex::BASE; v < nv+PointIndex::BASE; v++) - { - INDEX_3_CLOSED_HASHTABLE vert2face(2*max_face_on_vertex+10); - - for (int j = 0; j < vert2oldface[v].Size(); j++) - { - int fnr = vert2oldface[v][j]; - INDEX_3 face (face2vert[fnr].I1(), - face2vert[fnr].I2(), - face2vert[fnr].I3()); - vert2face.Set (face, fnr+1); - } - - - // cout << "inherited faces: " << endl << vert2face << endl; - - - for (int j = 0; j < (*vert2element)[v].Size(); j++) - { - int elnr = (*vert2element)[v][j]; - const Element & el = mesh.VolumeElement (elnr); - - int nelfaces = GetNFaces (el.GetType()); - const ELEMENT_FACE * elfaces = GetFaces1 (el.GetType()); - - for (int j = 0; j < nelfaces; j++) - if (elfaces[j][3] == 0) - - { // triangle - int facenum, facedir; - INDEX_3 face(el.PNum(elfaces[j][0]), - el.PNum(elfaces[j][1]), - el.PNum(elfaces[j][2])); - - facedir = 0; - if (face.I1() > face.I2()) - { swap (face.I1(), face.I2()); facedir += 1; } - if (face.I2() > face.I3()) - { swap (face.I2(), face.I3()); facedir += 2; } - if (face.I1() > face.I2()) - { swap (face.I1(), face.I2()); facedir += 4; } - - if (face.I1() != v) continue; - - if (vert2face.Used (face)) - facenum = vert2face.Get(face); - else - { - nfa++; - vert2face.Set (face, nfa); - facenum = nfa; - - INDEX_4 hface(face.I1(),face.I2(),face.I3(),0); - if (pass == 2) face2vert.Append (hface); - } - faces.Elem(elnr)[j] = 8*(facenum-1)+facedir+1; - } - - else - - { - // quad - int facenum, facedir; - INDEX_4Q face4(el.PNum(elfaces[j][0]), - el.PNum(elfaces[j][1]), - el.PNum(elfaces[j][2]), - el.PNum(elfaces[j][3])); - - facedir = 0; - if (min2 (face4.I1(), face4.I2()) > - min2 (face4.I4(), face4.I3())) - { // z - flip - facedir += 1; - swap (face4.I1(), face4.I4()); - swap (face4.I2(), face4.I3()); - } - if (min2 (face4.I1(), face4.I4()) > - min2 (face4.I2(), face4.I3())) - { // x - flip - facedir += 2; - swap (face4.I1(), face4.I2()); - swap (face4.I3(), face4.I4()); - } - if (face4.I2() > face4.I4()) - { // diagonal flip - facedir += 4; - swap (face4.I2(), face4.I4()); - } - - - INDEX_3 face(face4.I1(), face4.I2(), face4.I3()); - - if (face.I1() != v) continue; - - if (vert2face.Used (face)) - { - facenum = vert2face.Get(face); - } - else - { - nfa++; - vert2face.Set (face, nfa); - facenum = nfa; - - INDEX_4 hface(face4.I1(),face4.I2(),face4.I3(),face4.I4()); - if (pass == 2) face2vert.Append (hface); - } - - faces.Elem(elnr)[j] = 8*(facenum-1)+facedir+1; - } - } - - for (int j = 0; j < (*vert2surfelement)[v].Size(); j++) - { - int elnr = (*vert2surfelement)[v][j]; - // cout << "surfelnr = " << elnr << endl; - const Element2d & el = mesh.SurfaceElement (elnr); - - const ELEMENT_FACE * elfaces = GetFaces1 (el.GetType()); - - if (elfaces[0][3] == 0) - - { // triangle - - int facenum; - int facedir; - - INDEX_3 face(el.PNum(elfaces[0][0]), - el.PNum(elfaces[0][1]), - el.PNum(elfaces[0][2])); - - // cout << "face = " << face << endl; - - facedir = 0; - if (face.I1() > face.I2()) - { - swap (face.I1(), face.I2()); - facedir += 1; - } - if (face.I2() > face.I3()) - { - swap (face.I2(), face.I3()); - facedir += 2; - } - if (face.I1() > face.I2()) - { - swap (face.I1(), face.I2()); - facedir += 4; - } - - if (face.I1() != v) continue; - - if (vert2face.Used (face)) - facenum = vert2face.Get(face); - else - { - nfa++; - vert2face.Set (face, nfa); - facenum = nfa; - - INDEX_4 hface(face.I1(),face.I2(),face.I3(),0); - if (pass == 2) face2vert.Append (hface); - } - - // cout << "face = " << face << " selnr = " << elnr << endl; - surffaces.Elem(elnr) = 8*(facenum-1)+facedir+1; - // face2surfel.Elem(facenum) = elnr; - } - - else - - { - // quad - int facenum; - int facedir; - - INDEX_4Q face4(el.PNum(elfaces[0][0]), - el.PNum(elfaces[0][1]), - el.PNum(elfaces[0][2]), - el.PNum(elfaces[0][3])); - - facedir = 0; - if (min2 (face4.I1(), face4.I2()) > - min2 (face4.I4(), face4.I3())) - { // z - orientation - facedir += 1; - swap (face4.I1(), face4.I4()); - swap (face4.I2(), face4.I3()); - } - if (min2 (face4.I1(), face4.I4()) > - min2 (face4.I2(), face4.I3())) - { // x - orientation - facedir += 2; - swap (face4.I1(), face4.I2()); - swap (face4.I3(), face4.I4()); - } - if (face4.I2() > face4.I4()) - { - facedir += 4; - swap (face4.I2(), face4.I4()); - } - - INDEX_3 face(face4.I1(), face4.I2(), face4.I3()); - if (face.I1() != v) continue; - - if (vert2face.Used (face)) - facenum = vert2face.Get(face); - else - { - nfa++; - vert2face.Set (face, nfa); - facenum = nfa; - - INDEX_4 hface(face4.I1(),face4.I2(),face4.I3(),face4.I3()); - if (pass == 2) face2vert.Append (hface); - } - - surffaces.Elem(elnr) = 8*(facenum-1)+facedir+1; - } - } - } - - face2vert.SetAllocSize (nfa); - } - */ - - - - - - - - + + + + NgProfiler::StopTimer (timer2a); + NgProfiler::StartTimer (timer2b); + + /* for (int pass = 1; pass <= 2; pass++) { + cout << "pass = " << pass << endl; nfa = oldnfa; for (int v = PointIndex::BASE; v < nv+PointIndex::BASE; v++) { @@ -649,41 +434,44 @@ namespace netgen if (pass == 2) - for (int j = nfa; j < face2vert.Size(); j++) - { - if (face2vert[j][0] == v) - { - INDEX_3 face (face2vert[j].I1(), - face2vert[j].I2(), - face2vert[j].I3()); - vert2face.Set (face, j+1); - nfa++; - } - else - break; - } - + { + // *testout << "pass 2, nfa = " << nfa << "; face2vert.Size() = " << face2vert.Size() << endl; + for (int j = nfa; j < face2vert.Size(); j++) + { + if (face2vert[j][0] == v) + { + INDEX_3 face (face2vert[j].I1(), + face2vert[j].I2(), + face2vert[j].I3()); + vert2face.Set (face, j+1); + nfa++; + } + else + break; + } + } // cout << "inherited faces: " << endl << vert2face << endl; for (int j = 0; j < (*vert2element)[v].Size(); j++) { + // NgProfiler::RegionTimer reg3 (timer2d); + ElementIndex elnr = (*vert2element)[v][j]; const Element & el = mesh[elnr]; int nelfaces = GetNFaces (el.GetType()); - const ELEMENT_FACE * elfaces = GetFaces1 (el.GetType()); + const ELEMENT_FACE * elfaces = GetFaces0 (el.GetType()); for (int j = 0; j < nelfaces; j++) - if (elfaces[j][3] == 0) - + if (elfaces[j][3] < 0) + { // triangle int facenum, facedir; - INDEX_3 face(el.PNum(elfaces[j][0]), - el.PNum(elfaces[j][1]), - el.PNum(elfaces[j][2])); - + INDEX_3 face(el[elfaces[j][0]], el[elfaces[j][1]], + el[elfaces[j][2]]); + facedir = 0; if (face.I1() > face.I2()) { swap (face.I1(), face.I2()); facedir += 1; } @@ -693,11 +481,12 @@ namespace netgen { swap (face.I1(), face.I2()); facedir += 4; } if (face.I1() != v) continue; - + if (vert2face.Used (face)) facenum = vert2face.Get(face); else { + if (pass == 2) cout << "hier in pass 2" << endl; nfa++; vert2face.Set (face, nfa); facenum = nfa; @@ -707,16 +496,14 @@ namespace netgen } faces[elnr][j] = 8*(facenum-1)+facedir+1; } - + else { // quad int facenum, facedir; - INDEX_4Q face4(el.PNum(elfaces[j][0]), - el.PNum(elfaces[j][1]), - el.PNum(elfaces[j][2]), - el.PNum(elfaces[j][3])); + INDEX_4Q face4(el[elfaces[j][0]], el[elfaces[j][1]], + el[elfaces[j][2]], el[elfaces[j][3]]); facedir = 0; if (min2 (face4.I1(), face4.I2()) > @@ -750,6 +537,7 @@ namespace netgen } else { + if (pass == 2) cout << "hier in pass 2" << endl; nfa++; vert2face.Set (face, nfa); facenum = nfa; @@ -873,7 +661,7 @@ namespace netgen // *testout << "faces = " << face2vert << endl; if (pass == 1) { - // *testout << "sort from " << first_fa << " to " << nfa << endl; + // *testout << "pass1, sort from " << first_fa << " to " << nfa << endl; for (int i = first_fa; i < nfa; i++) for (int j = first_fa+1; j < nfa; j++) if (face2vert[j] < face2vert[j-1]) @@ -883,15 +671,344 @@ namespace netgen } - face2vert.SetAllocSize (nfa); - + face2vert.SetAllocSize (nfa); } + */ + + /* + // timing tests + static int timer_touch = NgProfiler::CreateTimer ("topology::buildfaces - touch els"); + static int timer_touch2 = NgProfiler::CreateTimer ("topology::buildfaces - touch els vert"); + + NgProfiler::StartTimer (timer_touch); + + int sum = 0; + for (ElementIndex ei = 0; ei < mesh.GetNE(); ei++) + { + const Element & el = mesh[ei]; + for (int j = 0; j < el.GetNP(); j++) + sum += el[j]; + } + + NgProfiler::StopTimer (timer_touch); + + NgProfiler::StartTimer (timer_touch2); + + for (int v = PointIndex::BASE; v < nv+PointIndex::BASE; v++) + { + INDEX_3_CLOSED_HASHTABLE vert2face(2*max_face_on_vertex+10); + + for (int pass = 1; pass <= 2; pass++) + for (int j = 0; j < (*vert2element)[v].Size(); j++) + { + // NgProfiler::RegionTimer reg3 (timer2d); + + ElementIndex elnr = (*vert2element)[v][j]; + const Element & el = mesh[elnr]; + + + int nelfaces = GetNFaces (el.GetType()); + const ELEMENT_FACE * elfaces = GetFaces0 (TET); + + for (int j = 0; j < 4; j++) + { // triangle + INDEX_3 face(el[elfaces[j][0]], el[elfaces[j][1]], + el[elfaces[j][2]]); + + int facedir = 0; + if (face.I1() > face.I2()) + { swap (face.I1(), face.I2()); facedir += 1; } + if (face.I2() > face.I3()) + { swap (face.I2(), face.I3()); facedir += 2; } + if (face.I1() > face.I2()) + { swap (face.I1(), face.I2()); facedir += 4; } + + sum += face.I1(); + sum += face.I1() < face.I2(); + } + } + } + + NgProfiler::StopTimer (timer_touch2); + *testout << "sum" << sum << endl; + */ + + + + nfa = oldnfa; + for (int v = PointIndex::BASE; v < nv+PointIndex::BASE; v++) + { + int first_fa = nfa; + + INDEX_3_CLOSED_HASHTABLE vert2face(2*max_face_on_vertex+10); + + for (int j = 0; j < vert2oldface[v].Size(); j++) + { + int fnr = vert2oldface[v][j]; + INDEX_3 face (face2vert[fnr].I1(), + face2vert[fnr].I2(), + face2vert[fnr].I3()); + vert2face.Set (face, fnr+1); + } + + + for (int pass = 1; pass <= 2; pass++) + { + + if (pass == 2) + { + for (int j = first_fa; j < face2vert.Size(); j++) + { + if (face2vert[j][0] == v) + { + INDEX_3 face (face2vert[j].I1(), + face2vert[j].I2(), + face2vert[j].I3()); + vert2face.Set (face, j+1); + } + else + break; + } + } + + + + for (int j = 0; j < (*vert2element)[v].Size(); j++) + { + // NgProfiler::RegionTimer reg3 (timer2d); + + ElementIndex elnr = (*vert2element)[v][j]; + const Element & el = mesh[elnr]; + + int nelfaces = GetNFaces (el.GetType()); + const ELEMENT_FACE * elfaces = GetFaces0 (el.GetType()); + + for (int j = 0; j < nelfaces; j++) + if (elfaces[j][3] < 0) + + { // triangle + INDEX_3 face(el[elfaces[j][0]], el[elfaces[j][1]], + el[elfaces[j][2]]); + + int facedir = 0; + if (face.I1() > face.I2()) + { swap (face.I1(), face.I2()); facedir += 1; } + if (face.I2() > face.I3()) + { swap (face.I2(), face.I3()); facedir += 2; } + if (face.I1() > face.I2()) + { swap (face.I1(), face.I2()); facedir += 4; } + + if (face.I1() != v) continue; + + if (pass == 1) + { + if (!vert2face.Used (face)) + { + nfa++; + vert2face.Set (face, nfa); + INDEX_4 hface(face.I1(),face.I2(),face.I3(),0); + face2vert.Append (hface); + } + } + else + { + int facenum = vert2face.Get(face); + // faces[elnr][j] = 8*(facenum-1)+facedir+1; + faces[elnr][j].fnr = facenum-1; + faces[elnr][j].forient = facedir; + } + } + + else + + { + // quad + int facenum; + INDEX_4Q face4(el[elfaces[j][0]], el[elfaces[j][1]], + el[elfaces[j][2]], el[elfaces[j][3]]); + + int facedir = 0; + if (min2 (face4.I1(), face4.I2()) > + min2 (face4.I4(), face4.I3())) + { // z - flip + facedir += 1; + swap (face4.I1(), face4.I4()); + swap (face4.I2(), face4.I3()); + } + if (min2 (face4.I1(), face4.I4()) > + min2 (face4.I2(), face4.I3())) + { // x - flip + facedir += 2; + swap (face4.I1(), face4.I2()); + swap (face4.I3(), face4.I4()); + } + if (face4.I2() > face4.I4()) + { // diagonal flip + facedir += 4; + swap (face4.I2(), face4.I4()); + } + + + INDEX_3 face(face4.I1(), face4.I2(), face4.I3()); + + if (face.I1() != v) continue; + + if (vert2face.Used (face)) + { + facenum = vert2face.Get(face); + } + else + { + if (pass == 2) cout << "hier in pass 2" << endl; + nfa++; + vert2face.Set (face, nfa); + facenum = nfa; + + INDEX_4 hface(face4.I1(),face4.I2(),face4.I3(),face4.I4()); + face2vert.Append (hface); + } + + // faces[elnr][j] = 8*(facenum-1)+facedir+1; + faces[elnr][j].fnr = facenum-1; + faces[elnr][j].forient = facedir; + } + } + + for (int j = 0; j < (*vert2surfelement)[v].Size(); j++) + { + int elnr = (*vert2surfelement)[v][j]; + // cout << "surfelnr = " << elnr << endl; + const Element2d & el = mesh.SurfaceElement (elnr); + + const ELEMENT_FACE * elfaces = GetFaces1 (el.GetType()); + + if (elfaces[0][3] == 0) + + { // triangle + + int facenum; + int facedir; + + INDEX_3 face(el.PNum(elfaces[0][0]), + el.PNum(elfaces[0][1]), + el.PNum(elfaces[0][2])); + + // cout << "face = " << face << endl; + + facedir = 0; + if (face.I1() > face.I2()) + { + swap (face.I1(), face.I2()); + facedir += 1; + } + if (face.I2() > face.I3()) + { + swap (face.I2(), face.I3()); + facedir += 2; + } + if (face.I1() > face.I2()) + { + swap (face.I1(), face.I2()); + facedir += 4; + } + + if (face.I1() != v) continue; + + if (vert2face.Used (face)) + facenum = vert2face.Get(face); + else + { + nfa++; + vert2face.Set (face, nfa); + facenum = nfa; + + INDEX_4 hface(face.I1(),face.I2(),face.I3(),0); + face2vert.Append (hface); + } + + // surffaces.Elem(elnr) = 8*(facenum-1)+facedir+1; + surffaces.Elem(elnr).fnr = facenum-1; + surffaces.Elem(elnr).forient = facedir; + } + + else + + { + // quad + int facenum; + int facedir; + + INDEX_4Q face4(el.PNum(elfaces[0][0]), + el.PNum(elfaces[0][1]), + el.PNum(elfaces[0][2]), + el.PNum(elfaces[0][3])); + + facedir = 0; + if (min2 (face4.I1(), face4.I2()) > + min2 (face4.I4(), face4.I3())) + { // z - orientation + facedir += 1; + swap (face4.I1(), face4.I4()); + swap (face4.I2(), face4.I3()); + } + if (min2 (face4.I1(), face4.I4()) > + min2 (face4.I2(), face4.I3())) + { // x - orientation + facedir += 2; + swap (face4.I1(), face4.I2()); + swap (face4.I3(), face4.I4()); + } + if (face4.I2() > face4.I4()) + { + facedir += 4; + swap (face4.I2(), face4.I4()); + } + + INDEX_3 face(face4.I1(), face4.I2(), face4.I3()); + if (face.I1() != v) continue; + + if (vert2face.Used (face)) + facenum = vert2face.Get(face); + else + { + nfa++; + vert2face.Set (face, nfa); + facenum = nfa; + + INDEX_4 hface(face4.I1(),face4.I2(),face4.I3(),face4.I3()); + face2vert.Append (hface); + } + + // surffaces.Elem(elnr) = 8*(facenum-1)+facedir+1; + surffaces.Elem(elnr).fnr = facenum-1; + surffaces.Elem(elnr).forient = facedir; + } + } + + // sort faces + if (pass == 1) + { + for (int i = 0; i < nfa-first_fa; i++) + for (int j = first_fa+1; j < nfa-i; j++) + if (face2vert[j] < face2vert[j-1]) + Swap (face2vert[j-1], face2vert[j]); + } + } + } + + + face2vert.SetAllocSize (nfa); + + + // *testout << "face2vert = " << endl << face2vert << endl; + NgProfiler::StopTimer (timer2b); + NgProfiler::StartTimer (timer2c); @@ -922,7 +1039,8 @@ namespace netgen for (int i = 1; i <= ne; i++) for (int j = 0; j < 6; j++) { - int fnum = (faces.Get(i)[j]+7) / 8; + // int fnum = (faces.Get(i)[j]+7) / 8; + int fnum = faces.Get(i)[j].fnr+1; if (fnum > 0 && face2surfel.Elem(fnum)) { int sel = face2surfel.Elem(fnum); @@ -1012,6 +1130,7 @@ namespace netgen if (cnt_err && ntasks == 1) cout << cnt_err << " elements are not matching !!!" << endl; } + NgProfiler::StopTimer (timer2c); } @@ -1155,7 +1274,8 @@ namespace netgen int ned = GetNEdges (mesh.VolumeElement(elnr).GetType()); eledges.SetSize (ned); for (int i = 0; i < ned; i++) - eledges[i] = abs (edges.Get(elnr)[i]); + eledges[i] = edges.Get(elnr)[i].nr+1; + // eledges[i] = abs (edges.Get(elnr)[i]); } void MeshTopology :: GetElementFaces (int elnr, Array & elfaces, bool withorientation) const { @@ -1166,18 +1286,23 @@ namespace netgen for (int i = 1; i <= nfa; i++) { - elfaces.Elem(i) = (faces.Get(elnr)[i-1]-1) / 8 + 1; + // elfaces.Elem(i) = (faces.Get(elnr)[i-1]-1) / 8 + 1; + elfaces.Elem(i) = faces.Get(elnr)[i-1].fnr+1; } else - - for (int i = 1; i <= nfa; i++) - { + { + cerr << "GetElementFaces with orientation currently not supported" << endl; + /* + for (int i = 1; i <= nfa; i++) + { elfaces.Elem(i) = (faces.Get(elnr)[i-1]-1) / 8 + 1; int orient = (faces.Get(elnr)[i-1]-1) % 8; if(orient == 1 || orient == 2 || orient == 4 || orient == 7) - elfaces.Elem(i) *= -1; - } + elfaces.Elem(i) *= -1; + } + */ + } } void MeshTopology :: GetElementEdgeOrientations (int elnr, Array & eorient) const @@ -1185,7 +1310,8 @@ namespace netgen int ned = GetNEdges (mesh.VolumeElement(elnr).GetType()); eorient.SetSize (ned); for (int i = 1; i <= ned; i++) - eorient.Elem(i) = (edges.Get(elnr)[i-1] > 0) ? 1 : -1; + // eorient.Elem(i) = (edges.Get(elnr)[i-1] > 0) ? 1 : -1; + eorient.Elem(i) = (edges.Get(elnr)[i-1].orient) ? -1 : 1; } void MeshTopology :: GetElementFaceOrientations (int elnr, Array & forient) const @@ -1193,7 +1319,8 @@ namespace netgen int nfa = GetNFaces (mesh.VolumeElement(elnr).GetType()); forient.SetSize (nfa); for (int i = 1; i <= nfa; i++) - forient.Elem(i) = (faces.Get(elnr)[i-1]-1) % 8; + forient.Elem(i) = faces.Get(elnr)[i-1].forient; + // forient.Elem(i) = (faces.Get(elnr)[i-1]-1) % 8; } @@ -1201,24 +1328,32 @@ namespace netgen int MeshTopology :: GetElementEdges (int elnr, int * eledges, int * orient) const { // int ned = GetNEdges (mesh.VolumeElement(elnr).GetType()); - + if (mesh.GetDimension()==3 || 1) { - if (orient) + if (orient) { for (int i = 0; i < 12; i++) { + /* if (!edges.Get(elnr)[i]) return i; eledges[i] = abs (edges.Get(elnr)[i]); orient[i] = (edges.Get(elnr)[i] > 0 ) ? 1 : -1; + */ + if (edges.Get(elnr)[i].nr == -1) return i; + eledges[i] = edges.Get(elnr)[i].nr+1; + orient[i] = edges.Get(elnr)[i].orient ? -1 : 1; } } else { for (int i = 0; i < 12; i++) { - if (!edges.Get(elnr)[i]) return i; - eledges[i] = abs (edges.Get(elnr)[i]); + // if (!edges.Get(elnr)[i]) return i; + // eledges[i] = abs (edges.Get(elnr)[i]); + if (edges.Get(elnr)[i].nr == -1) return i; + eledges[i] = edges.Get(elnr)[i].nr+1; + } } return 12; @@ -1255,17 +1390,24 @@ namespace netgen { for (int i = 0; i < 6; i++) { + /* if (!faces.Get(elnr)[i]) return i; elfaces[i] = (faces.Get(elnr)[i]-1) / 8 + 1; orient[i] = (faces.Get(elnr)[i]-1) % 8; + */ + if (faces.Get(elnr)[i].fnr == -1) return i; + elfaces[i] = faces.Get(elnr)[i].fnr+1; + orient[i] = faces.Get(elnr)[i].forient; } } else { for (int i = 0; i < 6; i++) { - if (!faces.Get(elnr)[i]) return i; - elfaces[i] = (faces.Get(elnr)[i]-1) / 8 + 1; + // if (!faces.Get(elnr)[i]) return i; + // elfaces[i] = (faces.Get(elnr)[i]-1) / 8 + 1; + if (faces.Get(elnr)[i].fnr == -1) return i; + elfaces[i] = faces.Get(elnr)[i].fnr+1; } } return 6; @@ -1276,7 +1418,8 @@ namespace netgen int ned = GetNEdges (mesh.SurfaceElement(elnr).GetType()); eledges.SetSize (ned); for (int i = 0; i < ned; i++) - eledges[i] = abs (surfedges.Get(elnr)[i]); + // eledges[i] = abs (surfedges.Get(elnr)[i]); + eledges[i] = surfedges.Get(elnr)[i].nr+1; } void MeshTopology :: GetEdges (SurfaceElementIndex elnr, Array & eledges) const @@ -1284,17 +1427,20 @@ namespace netgen int ned = GetNEdges (mesh[elnr].GetType()); eledges.SetSize (ned); for (int i = 0; i < ned; i++) - eledges[i] = abs (surfedges[elnr][i])-1; + // eledges[i] = abs (surfedges[elnr][i])-1; + eledges[i] = surfedges[elnr][i].nr; } int MeshTopology :: GetSurfaceElementFace (int elnr) const { - return (surffaces.Get(elnr)-1) / 8 + 1; + // return (surffaces.Get(elnr)-1) / 8 + 1; + return surffaces.Get(elnr).fnr+1; } int MeshTopology :: GetFace (SurfaceElementIndex elnr) const { - return (surffaces[elnr]-1) / 8; + // return (surffaces[elnr]-1) / 8; + return surffaces[elnr].fnr; } @@ -1304,12 +1450,14 @@ namespace netgen int ned = GetNEdges (mesh.SurfaceElement(elnr).GetType()); eorient.SetSize (ned); for (int i = 0; i < ned; i++) - eorient[i] = (surfedges.Get(elnr)[i] > 0) ? 1 : -1; + // eorient[i] = (surfedges.Get(elnr)[i] > 0) ? 1 : -1; + eorient[i] = (surfedges.Get(elnr)[i].orient) ? -1 : 1; } int MeshTopology :: GetSurfaceElementFaceOrientation (int elnr) const { - return (surffaces.Get(elnr)-1) % 8; + // return (surffaces.Get(elnr)-1) % 8; + return surffaces.Get(elnr).forient; } int MeshTopology :: GetSurfaceElementEdges (int elnr, int * eledges, int * orient) const @@ -1321,26 +1469,41 @@ namespace netgen { for (i = 0; i < 4; i++) { + /* if (!surfedges.Get(elnr)[i]) return i; eledges[i] = abs (surfedges.Get(elnr)[i]); orient[i] = (surfedges.Get(elnr)[i] > 0 ) ? 1 : -1; + */ + if (surfedges.Get(elnr)[i].nr == -1) return i; + eledges[i] = surfedges.Get(elnr)[i].nr+1; + orient[i] = (surfedges.Get(elnr)[i].orient) ? -1 : 1; + } } else { for (i = 0; i < 4; i++) { + /* if (!surfedges.Get(elnr)[i]) return i; eledges[i] = abs (surfedges.Get(elnr)[i]); + */ + if (surfedges.Get(elnr)[i].nr == -1) return i; + eledges[i] = surfedges.Get(elnr)[i].nr+1; } } return 4; } else { + /* eledges[0] = abs (segedges.Get(elnr)); if (orient) orient[0] = segedges.Get(elnr) > 0 ? 1 : -1; + */ + eledges[0] = segedges.Get(elnr).nr+1; + if (orient) + orient[0] = segedges.Get(elnr).orient ? -1 : 1; } return 1; } diff --git a/libsrc/meshing/topology.hpp b/libsrc/meshing/topology.hpp index 9339784d..fe5e3db0 100644 --- a/libsrc/meshing/topology.hpp +++ b/libsrc/meshing/topology.hpp @@ -13,6 +13,19 @@ */ +struct T_EDGE +{ + int orient:1; + int nr:31; // 0-based +}; + +struct T_FACE +{ + int forient:3; + int fnr:29; // 0-based +}; + + class MeshTopology { const Mesh & mesh; @@ -21,11 +34,11 @@ class MeshTopology Array edge2vert; Array face2vert; - Array edges; - Array faces; - Array surfedges; - Array segedges; - Array surffaces; + Array edges; + Array faces; + Array surfedges; + Array segedges; + Array surffaces; Array surf2volelement; Array face2surfel; TABLE *vert2element; @@ -66,14 +79,20 @@ public: inline static const ELEMENT_FACE * GetFaces0 (ELEMENT_TYPE et); - int GetSegmentEdge (int segnr) const { return abs(segedges[segnr-1]); } - int GetSegmentEdgeOrientation (int segnr) const { return sgn(segedges[segnr-1]); } - int GetEdge (SegmentIndex segnr) const { return abs(segedges[segnr])-1; } + // int GetSegmentEdge (int segnr) const { return abs(segedges[segnr-1]); } + // int GetSegmentEdgeOrientation (int segnr) const { return sgn(segedges[segnr-1]); } + // int GetEdge (SegmentIndex segnr) const { return abs(segedges[segnr])-1; } + int GetSegmentEdge (int segnr) const { return segedges[segnr-1].nr+1; } + int GetEdge (SegmentIndex segnr) const { return segedges[segnr].nr; } void GetSegmentEdge (int segnr, int & enr, int & orient) const { + /* enr = abs(segedges.Get(segnr)); orient = segedges.Get(segnr) > 0 ? 1 : -1; + */ + enr = segedges.Get(segnr).nr+1; + orient = segedges.Get(segnr).orient; } void GetElementEdges (int elnr, Array & edges) const; @@ -103,12 +122,12 @@ public: int GetSurfaceElementEdges (int elnr, int * edges, int * orient) const; - const int * GetElementEdgesPtr (int elnr) const { return &edges[elnr][0]; } - const int * GetSurfaceElementEdgesPtr (int selnr) const { return &surfedges[selnr][0]; } - const int * GetSegmentElementEdgesPtr (int selnr) const { return &segedges[selnr]; } + const T_EDGE * GetElementEdgesPtr (int elnr) const { return &edges[elnr][0]; } + const T_EDGE * GetSurfaceElementEdgesPtr (int selnr) const { return &surfedges[selnr][0]; } + const T_EDGE * GetSegmentElementEdgesPtr (int selnr) const { return &segedges[selnr]; } - const int * GetElementFacesPtr (int elnr) const { return &faces[elnr][0]; } - const int * GetSurfaceElementFacesPtr (int selnr) const { return &surffaces[selnr]; } + const T_FACE * GetElementFacesPtr (int elnr) const { return &faces[elnr][0]; } + const T_FACE * GetSurfaceElementFacesPtr (int selnr) const { return &surffaces[selnr]; } void GetSurface2VolumeElement (int selnr, int & elnr1, int & elnr2) const @@ -141,7 +160,7 @@ public: -int MeshTopology :: GetNVertices (ELEMENT_TYPE et) +inline int MeshTopology :: GetNVertices (ELEMENT_TYPE et) { switch (et) { @@ -172,14 +191,14 @@ int MeshTopology :: GetNVertices (ELEMENT_TYPE et) case HEX: return 8; - default: - cerr << "Ng_ME_GetNVertices, illegal element type " << et << endl; + // default: + // cerr << "Ng_ME_GetNVertices, illegal element type " << et << endl; } return 0; } -int MeshTopology :: GetNPoints (ELEMENT_TYPE et) +inline int MeshTopology :: GetNPoints (ELEMENT_TYPE et) { switch (et) { @@ -213,15 +232,15 @@ int MeshTopology :: GetNPoints (ELEMENT_TYPE et) case HEX: return 8; - default: - cerr << "Ng_ME_GetNVertices, illegal element type " << et << endl; + // default: + // cerr << "Ng_ME_GetNVertices, illegal element type " << et << endl; } return 0; } -int MeshTopology :: GetNEdges (ELEMENT_TYPE et) +inline int MeshTopology :: GetNEdges (ELEMENT_TYPE et) { switch (et) { @@ -252,14 +271,14 @@ int MeshTopology :: GetNEdges (ELEMENT_TYPE et) case HEX: return 12; - default: - cerr << "Ng_ME_GetNEdges, illegal element type " << et << endl; + // default: + // cerr << "Ng_ME_GetNEdges, illegal element type " << et << endl; } return 0; } -int MeshTopology :: GetNFaces (ELEMENT_TYPE et) +inline int MeshTopology :: GetNFaces (ELEMENT_TYPE et) { switch (et) { @@ -290,8 +309,8 @@ int MeshTopology :: GetNFaces (ELEMENT_TYPE et) case HEX: return 6; - default: - cerr << "Ng_ME_GetNVertices, illegal element type " << et << endl; + // default: + // cerr << "Ng_ME_GetNVertices, illegal element type " << et << endl; } return 0; } @@ -391,8 +410,8 @@ const ELEMENT_EDGE * MeshTopology :: GetEdges1 (ELEMENT_TYPE et) case HEX: return hex_edges; - default: - cerr << "Ng_ME_GetEdges, illegal element type " << et << endl; + // default: + // cerr << "Ng_ME_GetEdges, illegal element type " << et << endl; } return 0; } @@ -489,8 +508,8 @@ const ELEMENT_EDGE * MeshTopology :: GetEdges0 (ELEMENT_TYPE et) case HEX: return hex_edges; - default: - cerr << "Ng_ME_GetEdges, illegal element type " << et << endl; + // default: + // cerr << "Ng_ME_GetEdges, illegal element type " << et << endl; } return 0; } @@ -504,7 +523,7 @@ const ELEMENT_EDGE * MeshTopology :: GetEdges0 (ELEMENT_TYPE et) -const ELEMENT_FACE * MeshTopology :: GetFaces1 (ELEMENT_TYPE et) +inline const ELEMENT_FACE * MeshTopology :: GetFaces1 (ELEMENT_TYPE et) { static const int trig_faces[1][4] = { { 1, 2, 3, 0 } }; @@ -576,8 +595,8 @@ const ELEMENT_FACE * MeshTopology :: GetFaces1 (ELEMENT_TYPE et) case HEX: return hex_faces; - default: - cerr << "Ng_ME_GetVertices, illegal element type " << et << endl; + // default: + // cerr << "Ng_ME_GetVertices, illegal element type " << et << endl; } return 0; } @@ -586,7 +605,7 @@ const ELEMENT_FACE * MeshTopology :: GetFaces1 (ELEMENT_TYPE et) -const ELEMENT_FACE * MeshTopology :: GetFaces0 (ELEMENT_TYPE et) +inline const ELEMENT_FACE * MeshTopology :: GetFaces0 (ELEMENT_TYPE et) { static const int trig_faces[1][4] = { { 0, 1, 2, -1 } }; @@ -658,8 +677,8 @@ const ELEMENT_FACE * MeshTopology :: GetFaces0 (ELEMENT_TYPE et) case HEX: return hex_faces; - default: - cerr << "Ng_ME_GetVertices, illegal element type " << et << endl; + // default: + // cerr << "Ng_ME_GetVertices, illegal element type " << et << endl; } return 0; }