From d5f1076e173e9e312c4563f7df11ec3fe58b87ff Mon Sep 17 00:00:00 2001 From: Christopher Lackner Date: Wed, 6 Feb 2019 19:13:51 +0100 Subject: [PATCH] read tet10, pyramid13, prism15 and hex20 from gmsh --- libsrc/general/array.hpp | 7 +- libsrc/include/nginterface.h | 4 +- libsrc/include/nginterface_v2.hpp | 2 +- libsrc/meshing/curvedelems.cpp | 127 ++++++++++++++++++++++++++++ libsrc/meshing/meshclass.cpp | 6 +- libsrc/meshing/meshtype.cpp | 89 ++++++++++++++++++- libsrc/meshing/meshtype.hpp | 15 ++-- libsrc/meshing/python_mesh.cpp | 51 +++-------- libsrc/meshing/secondorder.cpp | 38 +++++++++ libsrc/meshing/topology.hpp | 25 +++++- libsrc/visualization/vsmesh.cpp | 2 +- libsrc/visualization/vssolution.cpp | 11 ++- python/read_gmsh.py | 108 ++++++++++++++++------- 13 files changed, 396 insertions(+), 89 deletions(-) diff --git a/libsrc/general/array.hpp b/libsrc/general/array.hpp index 16bb6e78..df0cd0a5 100644 --- a/libsrc/general/array.hpp +++ b/libsrc/general/array.hpp @@ -243,10 +243,11 @@ namespace netgen } explicit Array(size_t asize) - : FlatArray (asize, new T[asize]) + : FlatArray (asize, asize ? new T[asize] : nullptr) { - allocsize = asize; - ownmem = 1; + allocsize = asize; + if(asize) + ownmem = 1; } /// Generate array in user data diff --git a/libsrc/include/nginterface.h b/libsrc/include/nginterface.h index 38e38ec2..24e90cc1 100644 --- a/libsrc/include/nginterface.h +++ b/libsrc/include/nginterface.h @@ -32,7 +32,7 @@ // max number of nodes per element -#define NG_ELEMENT_MAXPOINTS 12 +#define NG_ELEMENT_MAXPOINTS 20 // max number of nodes per surface element #define NG_SURFACE_ELEMENT_MAXPOINTS 8 @@ -49,7 +49,7 @@ enum NG_ELEMENT_TYPE { NG_SEGM = 1, NG_SEGM3 = 2, NG_TRIG = 10, NG_QUAD=11, NG_TRIG6 = 12, NG_QUAD6 = 13, NG_QUAD8 = 14, NG_TET = 20, NG_TET10 = 21, - NG_PYRAMID = 22, NG_PRISM = 23, NG_PRISM12 = 24, + NG_PYRAMID = 22, NG_PRISM = 23, NG_PRISM12 = 24, NG_PRISM15 = 27, NG_PYRAMID13 = 28, NG_HEX = 25, NG_HEX20 = 26 }; diff --git a/libsrc/include/nginterface_v2.hpp b/libsrc/include/nginterface_v2.hpp index d6d4c5a8..962f061c 100644 --- a/libsrc/include/nginterface_v2.hpp +++ b/libsrc/include/nginterface_v2.hpp @@ -19,7 +19,7 @@ enum NG_ELEMENT_TYPE { NG_SEGM = 1, NG_SEGM3 = 2, NG_TRIG = 10, NG_QUAD=11, NG_TRIG6 = 12, NG_QUAD6 = 13, NG_QUAD8 = 14, NG_TET = 20, NG_TET10 = 21, - NG_PYRAMID = 22, NG_PRISM = 23, NG_PRISM12 = 24, + NG_PYRAMID = 22, NG_PRISM = 23, NG_PRISM12 = 24, NG_PRISM15 = 27, NG_PYRAMID13 = 28, NG_HEX = 25, NG_HEX20 = 26 }; diff --git a/libsrc/meshing/curvedelems.cpp b/libsrc/meshing/curvedelems.cpp index a3e7f209..182945d3 100644 --- a/libsrc/meshing/curvedelems.cpp +++ b/libsrc/meshing/curvedelems.cpp @@ -2790,6 +2790,32 @@ namespace netgen break; } + case PRISM15: + { + shapes = 0.0; + T x = xi(0); + T y = xi(1); + T z = xi(2); + T lam = 1-x-y; + T lamz = 1-z; + shapes[0] = (2*x*x-x) * (2*lamz*lamz-lamz); + shapes[1] = (2*y*y-y) * (2*lamz*lamz-lamz); + shapes[2] = (2*lam*lam-lam) * (2*lamz*lamz-lamz); + shapes[3] = (2*x*x-x) * (2*z*z-z); + shapes[4] = (2*y*y-y) * (2*z*z-z); + shapes[5] = (2*lam*lam-lam) * (2*z*z-z); + shapes[6] = 4 * x * y * (2*lamz*lamz-lamz); + shapes[7] = 4 * x * lam * (2*lamz*lamz-lamz); + shapes[8] = 4 * y * lam * (2*lamz*lamz-lamz); + shapes[9] = x * 4 * z * (1-z); + shapes[10] = y * 4 * z * (1-z); + shapes[11] = lam * 4 * z * (1-z); + shapes[12] = 4 * x * y * (2*z*z-z); + shapes[13] = 4 * x * lam * (2*z*z-z); + shapes[14] = 4 * y * lam * (2*z*z-z); + break; + } + case PYRAMID: { shapes = 0.0; @@ -2839,6 +2865,29 @@ namespace netgen break; } + case PYRAMID13: + { + shapes = 0.0; + T x = xi(0); + T y = xi(1); + T z = xi(2); + z *= 1-1e-12; + shapes[0] = (-z + z*(2*x + z - 1)*(2*y + z - 1)/(-z + 1) + (-2*x - z + 2)*(-2*y - z + 2))*(-0.5*x - 0.5*y - 0.5*z + 0.25); + shapes[1] = (0.5*x - 0.5*y - 0.25)*(-z - z*(2*x + z - 1)*(2*y + z - 1)/(-z + 1) + (2*x + z)*(-2*y - z + 2)); + shapes[2] = (-z + z*(2*x + z - 1)*(2*y + z - 1)/(-z + 1) + (2*x + z)*(2*y + z))*(0.5*x + 0.5*y + 0.5*z - 0.75); + shapes[3] = (-0.5*x + 0.5*y - 0.25)*(-z - z*(2*x + z - 1)*(2*y + z - 1)/(-z + 1) + (2*y + z)*(-2*x - z + 2)); + shapes[4] = z*(2*z - 1); + shapes[5] = 2*x*(-2*x - 2*z + 2)*(-2*y - 2*z + 2)/(-2*z + 2); + shapes[6] = 4*x*y*(-2*x - 2*z + 2)/(-2*z + 2); + shapes[7] = 2*y*(-2*x - 2*z + 2)*(-2*y - 2*z + 2)/(-2*z + 2); + shapes[8] = 4*x*y*(-2*y - 2*z + 2)/(-2*z + 2); + shapes[9] = z*(-2*x - 2*z + 2)*(-2*y - 2*z + 2)/(-z + 1); + shapes[10] = 2*x*z*(-2*y - 2*z + 2)/(-z + 1); + shapes[11] = 4*x*y*z/(-z + 1); + shapes[12] = 2*y*z*(-2*x - 2*z + 2)/(-z + 1); + break; + } + case HEX: { shapes = 0.0; @@ -3300,6 +3349,36 @@ namespace netgen } + case PRISM15: + { + AutoDiff<3,T> x(xi(0), 0); + AutoDiff<3,T> y(xi(1), 1); + AutoDiff<3,T> z(xi(2), 2); + AutoDiff<3,T> ad[15]; + AutoDiff<3,T> lam = 1-x-y; + AutoDiff<3,T> lamz = 1-z; + + ad[0] = (2*x*x-x) * (2*lamz*lamz-lamz); + ad[1] = (2*y*y-y) * (2*lamz*lamz-lamz); + ad[2] = (2*lam*lam-lam) * (2*lamz*lamz-lamz); + ad[3] = (2*x*x-x) * (2*z*z-z); + ad[4] = (2*y*y-y) * (2*z*z-z); + ad[5] = (2*lam*lam-lam) * (2*z*z-z); + ad[6] = 4 * x * y * (2*lamz*lamz-lamz); + ad[7] = 4 * x * lam * (2*lamz*lamz-lamz); + ad[8] = 4 * y * lam * (2*lamz*lamz-lamz); + ad[9] = x * 4 * z * (1-z); + ad[10] = y * 4 * z * (1-z); + ad[11] = lam * 4 * z * (1-z); + ad[12] = 4 * x * y * (2*z*z-z); + ad[13] = 4 * x * lam * (2*z*z-z); + ad[14] = 4 * y * lam * (2*z*z-z); + + for(int i=0; i<15; i++) + for(int j=0; j<3; j++) + dshapes(i,j) = ad[i].DValue(j); + break; + } case PYRAMID: { // if (typeid(T) == typeid(SIMD)) return; @@ -3397,6 +3476,54 @@ namespace netgen break; } + case PYRAMID13: + { + T x = xi(0); + T y = xi(1); + T z = xi(2); + z *= 1-1e-12; + dshapes(0,0) = 0.5*z - 0.5*z*(2*x + z - 1)*(2*y + z - 1)/(-z + 1) - 0.5*(-2*x - z + 2)*(-2*y - z + 2) + (-0.5*x - 0.5*y - 0.5*z + 0.25)*(4*y + 2*z + 2*z*(2*y + z - 1)/(-z + 1) - 4); + dshapes(0,1) = 0.5*z - 0.5*z*(2*x + z - 1)*(2*y + z - 1)/(-z + 1) - 0.5*(-2*x - z + 2)*(-2*y - z + 2) + (-0.5*x - 0.5*y - 0.5*z + 0.25)*(4*x + 2*z + 2*z*(2*x + z - 1)/(-z + 1) - 4); + dshapes(0,2) = 0.5*z - 0.5*z*(2*x + z - 1)*(2*y + z - 1)/(-z + 1) - 0.5*(-2*x - z + 2)*(-2*y - z + 2) + (-0.5*x - 0.5*y - 0.5*z + 0.25)*(2*x + 2*y + 2*z + z*(2*x + z - 1)/(-z + 1) + z*(2*y + z - 1)/(-z + 1) + z*(2*x + z - 1)*(2*y + z - 1)/((-z + 1)*(-z + 1)) - 5 + (2*x + z - 1)*(2*y + z - 1)/(-z + 1)); + dshapes(1,0) = -0.5*z - 0.5*z*(2*x + z - 1)*(2*y + z - 1)/(-z + 1) + 0.5*(2*x + z)*(-2*y - z + 2) + (0.5*x - 0.5*y - 0.25)*(-4*y - 2*z - 2*z*(2*y + z - 1)/(-z + 1) + 4); + dshapes(1,1) = 0.5*z + 0.5*z*(2*x + z - 1)*(2*y + z - 1)/(-z + 1) - 0.5*(2*x + z)*(-2*y - z + 2) + (-4*x - 2*z - 2*z*(2*x + z - 1)/(-z + 1))*(0.5*x - 0.5*y - 0.25); + dshapes(1,2) = (0.5*x - 0.5*y - 0.25)*(-2*x - 2*y - 2*z - z*(2*x + z - 1)/(-z + 1) - z*(2*y + z - 1)/(-z + 1) - z*(2*x + z - 1)*(2*y + z - 1)/((-z + 1)*(-z + 1)) + 1 - (2*x + z - 1)*(2*y + z - 1)/(-z + 1)); + dshapes(2,0) = -0.5*z + 0.5*z*(2*x + z - 1)*(2*y + z - 1)/(-z + 1) + 0.5*(2*x + z)*(2*y + z) + (4*y + 2*z + 2*z*(2*y + z - 1)/(-z + 1))*(0.5*x + 0.5*y + 0.5*z - 0.75); + dshapes(2,1) = -0.5*z + 0.5*z*(2*x + z - 1)*(2*y + z - 1)/(-z + 1) + 0.5*(2*x + z)*(2*y + z) + (4*x + 2*z + 2*z*(2*x + z - 1)/(-z + 1))*(0.5*x + 0.5*y + 0.5*z - 0.75); + dshapes(2,2) = -0.5*z + 0.5*z*(2*x + z - 1)*(2*y + z - 1)/(-z + 1) + 0.5*(2*x + z)*(2*y + z) + (0.5*x + 0.5*y + 0.5*z - 0.75)*(2*x + 2*y + 2*z + z*(2*x + z - 1)/(-z + 1) + z*(2*y + z - 1)/(-z + 1) + z*(2*x + z - 1)*(2*y + z - 1)/((-z + 1)*(-z + 1)) - 1 + (2*x + z - 1)*(2*y + z - 1)/(-z + 1)); + dshapes(3,0) = 0.5*z + 0.5*z*(2*x + z - 1)*(2*y + z - 1)/(-z + 1) - 0.5*(2*y + z)*(-2*x - z + 2) + (-0.5*x + 0.5*y - 0.25)*(-4*y - 2*z - 2*z*(2*y + z - 1)/(-z + 1)); + dshapes(3,1) = -0.5*z - 0.5*z*(2*x + z - 1)*(2*y + z - 1)/(-z + 1) + 0.5*(2*y + z)*(-2*x - z + 2) + (-0.5*x + 0.5*y - 0.25)*(-4*x - 2*z - 2*z*(2*x + z - 1)/(-z + 1) + 4); + dshapes(3,2) = (-0.5*x + 0.5*y - 0.25)*(-2*x - 2*y - 2*z - z*(2*x + z - 1)/(-z + 1) - z*(2*y + z - 1)/(-z + 1) - z*(2*x + z - 1)*(2*y + z - 1)/((-z + 1)*(-z + 1)) + 1 - (2*x + z - 1)*(2*y + z - 1)/(-z + 1)); + dshapes(4,0) = 0; + dshapes(4,1) = 0; + dshapes(4,2) = 4*z - 1; + dshapes(5,0) = -4*x*(-2*y - 2*z + 2)/(-2*z + 2) + 2*(-2*x - 2*z + 2)*(-2*y - 2*z + 2)/(-2*z + 2); + dshapes(5,1) = -4*x*(-2*x - 2*z + 2)/(-2*z + 2); + dshapes(5,2) = -4*x*(-2*x - 2*z + 2)/(-2*z + 2) - 4*x*(-2*y - 2*z + 2)/(-2*z + 2) + 4*x*(-2*x - 2*z + 2)*(-2*y - 2*z + 2)/((-2*z + 2)*(-2*z + 2)); + dshapes(6,0) = -8*x*y/(-2*z + 2) + 4*y*(-2*x - 2*z + 2)/(-2*z + 2); + dshapes(6,1) = 4*x*(-2*x - 2*z + 2)/(-2*z + 2); + dshapes(6,2) = -8*x*y/(-2*z + 2) + 8*x*y*(-2*x - 2*z + 2)/((-2*z + 2)*(-2*z + 2)); + dshapes(7,0) = -4*y*(-2*y - 2*z + 2)/(-2*z + 2); + dshapes(7,1) = -4*y*(-2*x - 2*z + 2)/(-2*z + 2) + 2*(-2*x - 2*z + 2)*(-2*y - 2*z + 2)/(-2*z + 2); + dshapes(7,2) = -4*y*(-2*x - 2*z + 2)/(-2*z + 2) - 4*y*(-2*y - 2*z + 2)/(-2*z + 2) + 4*y*(-2*x - 2*z + 2)*(-2*y - 2*z + 2)/((-2*z + 2)*(-2*z + 2)); + dshapes(8,0) = 4*y*(-2*y - 2*z + 2)/(-2*z + 2); + dshapes(8,1) = -8*x*y/(-2*z + 2) + 4*x*(-2*y - 2*z + 2)/(-2*z + 2); + dshapes(8,2) = -8*x*y/(-2*z + 2) + 8*x*y*(-2*y - 2*z + 2)/((-2*z + 2)*(-2*z + 2)); + dshapes(9,0) = -2*z*(-2*y - 2*z + 2)/(-z + 1); + dshapes(9,1) = -2*z*(-2*x - 2*z + 2)/(-z + 1); + dshapes(9,2) = -2*z*(-2*x - 2*z + 2)/(-z + 1) - 2*z*(-2*y - 2*z + 2)/(-z + 1) + z*(-2*x - 2*z + 2)*(-2*y - 2*z + 2)/((-z + 1)*(-z + 1)) + (-2*x - 2*z + 2)*(-2*y - 2*z + 2)/(-z + 1); + dshapes(10,0) = 2*z*(-2*y - 2*z + 2)/(-z + 1); + dshapes(10,1) = -4*x*z/(-z + 1); + dshapes(10,2) = -4*x*z/(-z + 1) + 2*x*z*(-2*y - 2*z + 2)/((-z + 1)*(-z + 1)) + 2*x*(-2*y - 2*z + 2)/(-z + 1); + dshapes(11,0) = 4*y*z/(-z + 1); + dshapes(11,1) = 4*x*z/(-z + 1); + dshapes(11,2) = 4*x*y*z/((-z + 1)*(-z + 1)) + 4*x*y/(-z + 1); + dshapes(12,0) = -4*y*z/(-z + 1); + dshapes(12,1) = 2*z*(-2*x - 2*z + 2)/(-z + 1); + dshapes(12,2) = -4*y*z/(-z + 1) + 2*y*z*(-2*x - 2*z + 2)/((-z + 1)*(-z + 1)) + 2*y*(-2*x - 2*z + 2)/(-z + 1); + break; + } + case HEX: { // if (typeid(T) == typeid(SIMD)) return; diff --git a/libsrc/meshing/meshclass.cpp b/libsrc/meshing/meshclass.cpp index 0d8f5879..807de27e 100644 --- a/libsrc/meshing/meshclass.cpp +++ b/libsrc/meshing/meshclass.cpp @@ -4712,7 +4712,7 @@ namespace netgen lam(2) > -eps && lam(0) + lam(1) + lam(2) < 1+eps); } - else if (el.GetType() == PRISM) + else if (el.GetType() == PRISM || el.GetType() == PRISM15) { retval = (lam(0) > -eps && lam(1) > -eps && @@ -4720,7 +4720,7 @@ namespace netgen lam(2) < 1+eps && lam(0) + lam(1) < 1+eps); } - else if (el.GetType() == PYRAMID) + else if (el.GetType() == PYRAMID || el.GetType() == PYRAMID13) { retval = (lam(0) > -eps && lam(1) > -eps && @@ -4728,7 +4728,7 @@ namespace netgen lam(0) + lam(2) < 1+eps && lam(1) + lam(2) < 1+eps); } - else if (el.GetType() == HEX) + else if (el.GetType() == HEX || el.GetType() == HEX20) { retval = (lam(0) > -eps && lam(0) < 1+eps && lam(1) > -eps && lam(1) < 1+eps && diff --git a/libsrc/meshing/meshtype.cpp b/libsrc/meshing/meshtype.cpp index 494a3879..b28cc92d 100644 --- a/libsrc/meshing/meshtype.cpp +++ b/libsrc/meshing/meshtype.cpp @@ -1028,6 +1028,8 @@ namespace netgen case 6: typ = PRISM; break; case 8: typ = HEX; break; case 10: typ = TET10; break; + case 13: typ = PYRAMID13; break; + case 15: typ = PRISM15; break; case 20: typ = HEX20; break; default: cerr << "Element::Element: unknown element with " << np << " points" << endl; } @@ -1109,6 +1111,8 @@ namespace netgen case 6: typ = PRISM; break; case 8: typ = HEX; break; case 10: typ = TET10; break; + case 13: typ = PYRAMID13; break; + case 15: typ = PRISM15; break; case 20: typ = HEX20; break; // default: break; @@ -1128,7 +1132,9 @@ namespace netgen case PRISM: np = 6; break; case HEX: np = 8; break; case TET10: np = 10; break; + case PYRAMID13: np = 13; break; case PRISM12: np = 12; break; + case PRISM15: np = 15; break; case HEX20: np = 20; break; default: break; @@ -1958,7 +1964,27 @@ namespace netgen shape(4) = p(2); break; } - + case PYRAMID13: + { + T x = p(0); + T y = p(1); + T z = p(2); + z *= 1-1e-12; + shape[0] = (-z + z*(2*x + z - 1)*(2*y + z - 1)/(-z + 1) + (-2*x - z + 2)*(-2*y - z + 2))*(-0.5*x - 0.5*y - 0.5*z + 0.25); + shape[1] = (0.5*x - 0.5*y - 0.25)*(-z - z*(2*x + z - 1)*(2*y + z - 1)/(-z + 1) + (2*x + z)*(-2*y - z + 2)); + shape[2] = (-z + z*(2*x + z - 1)*(2*y + z - 1)/(-z + 1) + (2*x + z)*(2*y + z))*(0.5*x + 0.5*y + 0.5*z - 0.75); + shape[3] = (-0.5*x + 0.5*y - 0.25)*(-z - z*(2*x + z - 1)*(2*y + z - 1)/(-z + 1) + (2*y + z)*(-2*x - z + 2)); + shape[4] = z*(2*z - 1); + shape[5] = 2*x*(-2*x - 2*z + 2)*(-2*y - 2*z + 2)/(-2*z + 2); + shape[6] = 4*x*y*(-2*x - 2*z + 2)/(-2*z + 2); + shape[7] = 2*y*(-2*x - 2*z + 2)*(-2*y - 2*z + 2)/(-2*z + 2); + shape[8] = 4*x*y*(-2*y - 2*z + 2)/(-2*z + 2); + shape[9] = z*(-2*x - 2*z + 2)*(-2*y - 2*z + 2)/(-z + 1); + shape[10] = 2*x*z*(-2*y - 2*z + 2)/(-z + 1); + shape[11] = 4*x*y*z/(-z + 1); + shape[12] = 2*y*z*(-2*x - 2*z + 2)/(-z + 1); + break; + } case PRISM: { shape(0) = p(0) * (1-p(2)); @@ -1969,6 +1995,30 @@ namespace netgen shape(5) = (1-p(0)-p(1)) * p(2); break; } + case PRISM15: + { + T x = p(0); + T y = p(1); + T z = p(2); + T lam = 1-x-y; + T lamz = 1-z; + shape[0] = (2*x*x-x) * (2*lamz*lamz-lamz); + shape[1] = (2*y*y-y) * (2*lamz*lamz-lamz); + shape[2] = (2*lam*lam-lam) * (2*lamz*lamz-lamz); + shape[3] = (2*x*x-x) * (2*z*z-z); + shape[4] = (2*y*y-y) * (2*z*z-z); + shape[5] = (2*lam*lam-lam) * (2*z*z-z); + shape[6] = 4 * x * y * (2*lamz*lamz-lamz); + shape[7] = 4 * x * lam * (2*lamz*lamz-lamz); + shape[8] = 4 * y * lam * (2*lamz*lamz-lamz); + shape[9] = x * 4 * z * (1-z); + shape[10] = y * 4 * z * (1-z); + shape[11] = lam * 4 * z * (1-z); + shape[12] = 4 * x * y * (2*z*z-z); + shape[13] = 4 * x * lam * (2*z*z-z); + shape[14] = 4 * y * lam * (2*z*z-z); + break; + } case HEX: { shape(0) = (1-p(0))*(1-p(1))*(1-p(2)); @@ -1981,6 +2031,43 @@ namespace netgen shape(7) = (1-p(0))*( p(1))*( p(2)); break; } + case HEX20: + { + T x = p(0); + T y = p(1); + T z = p(2); + shape[0] = (1-x)*(1-y)*(1-z); + shape[1] = x *(1-y)*(1-z); + shape[2] = x * y *(1-z); + shape[3] = (1-x)* y *(1-z); + shape[4] = (1-x)*(1-y)*(z); + shape[5] = x *(1-y)*(z); + shape[6] = x * y *(z); + shape[7] = (1-x)* y *(z); + + T sigma[8]={(1-x)+(1-y)+(1-z),x+(1-y)+(1-z),x+y+(1-z),(1-x)+y+(1-z), + (1-x)+(1-y)+z,x+(1-y)+z,x+y+z,(1-x)+y+z}; + + static const int e[12][2] = + { + { 0, 1 }, { 2, 3 }, { 3, 0 }, { 1, 2 }, + { 4, 5 }, { 6, 7 }, { 7, 4 }, { 5, 6 }, + { 0, 4 }, { 1, 5 }, { 2, 6 }, { 3, 7 }, + }; + + for (int i = 0; i < 12; i++) + { + T lame = shape[e[i][0]]+shape[e[i][1]]; + T xi = sigma[e[i][1]]-sigma[e[i][0]]; + shape[8+i] = (1-xi*xi)*lame; + } + for (int i = 0; i < 12; i++) + { + shape[e[i][0]] -= 0.5 * shape[8+i]; + shape[e[i][1]] -= 0.5 * shape[8+i]; + } + break; + } default: throw NgException("Element :: GetNewShape not implemented for that element"); } diff --git a/libsrc/meshing/meshtype.hpp b/libsrc/meshing/meshtype.hpp index 1b31d597..df997e36 100644 --- a/libsrc/meshing/meshtype.hpp +++ b/libsrc/meshing/meshtype.hpp @@ -21,7 +21,7 @@ namespace netgen SEGMENT = 1, SEGMENT3 = 2, TRIG = 10, QUAD=11, TRIG6 = 12, QUAD6 = 13, QUAD8 = 14, TET = 20, TET10 = 21, - PYRAMID = 22, PRISM = 23, PRISM12 = 24, + PYRAMID = 22, PRISM = 23, PRISM12 = 24, PRISM15 = 27, PYRAMID13 = 28, HEX = 25, HEX20 = 26 }; @@ -708,16 +708,18 @@ namespace netgen /// uint8_t GetNV() const { - __assume(typ >= TET && typ <= HEX20); + __assume(typ >= TET && typ <= PYRAMID13); switch (typ) { case TET: case TET10: return 4; - case PRISM12: - case PRISM: + case PRISM12: + case PRISM15: + case PRISM: return 6; case PYRAMID: + case PYRAMID13: return 5; case HEX: case HEX20: @@ -795,8 +797,9 @@ namespace netgen { case TET: case TET10: return 4; - case PYRAMID: return 5; - case PRISM: + case PYRAMID: case PYRAMID13: return 5; + case PRISM: + case PRISM15: case PRISM12: return 5; case HEX: case HEX20: return 6; diff --git a/libsrc/meshing/python_mesh.cpp b/libsrc/meshing/python_mesh.cpp index dbc498cb..45929e4b 100644 --- a/libsrc/meshing/python_mesh.cpp +++ b/libsrc/meshing/python_mesh.cpp @@ -256,44 +256,19 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m) py::class_(m, "Element3D") .def(py::init([](int index, py::list vertices) { - Element * newel = nullptr; - if (py::len(vertices) == 4) - { - newel = new Element(TET); - for (int i = 0; i < 4; i++) - (*newel)[i] = py::extract(vertices[i])(); - newel->SetIndex(index); - } - else if (py::len(vertices) == 5) - { - newel = new Element(PYRAMID); - for (int i = 0; i < 5; i++) - (*newel)[i] = py::extract(vertices[i])(); - newel->SetIndex(index); - } - else if (py::len(vertices) == 6) - { - newel = new Element(PRISM); - for (int i = 0; i < 6; i++) - (*newel)[i] = py::extract(vertices[i])(); - newel->SetIndex(index); - } - else if (py::len(vertices) == 8) - { - newel = new Element(HEX); - for (int i = 0; i < 8; i++) - (*newel)[i] = py::extract(vertices[i])(); - newel->SetIndex(index); - } - else if (py::len(vertices) == 20) - { - newel = new Element(HEX20); - for (int i = 0; i < 20; i++) - (*newel)[i] = py::extract(vertices[i])(); - newel->SetIndex(index); - } - else - throw NgException ("cannot create element"); + std::map types = {{4, TET}, + {5, PYRAMID}, + {6, PRISM}, + {8, HEX}, + {10, TET10}, + {13, PYRAMID13}, + {15, PRISM15}, + {20, HEX20}}; + int np = py::len(vertices); + auto newel = new Element(types[np]); + for(int i=0; i(vertices[i]); + newel->SetIndex(index); return newel; }), py::arg("index")=1,py::arg("vertices"), diff --git a/libsrc/meshing/secondorder.cpp b/libsrc/meshing/secondorder.cpp index 6241fd4e..49f91662 100644 --- a/libsrc/meshing/secondorder.cpp +++ b/libsrc/meshing/secondorder.cpp @@ -223,6 +223,29 @@ namespace netgen { 3, 4, 10 }, { 4, 5, 11 }, }; + static int betw_prism15[9][3] = + { + { 0, 1, 6 }, + { 0, 2, 7 }, + { 1, 2, 8 }, + { 0, 3, 9 }, + { 1, 4, 10 }, + { 2, 5, 11 }, + { 3, 4, 12 }, + { 3, 5, 13 }, + { 4, 5, 14 } + }; + static int betw_pyramid[8][3] = + { + { 0, 1, 5 }, + { 3, 2, 6 }, + { 3, 0, 7 }, + { 1, 2, 8 }, + { 0, 4, 9 }, + { 1, 4, 10 }, + { 2, 4, 11 }, + { 3, 4, 12 } + }; static int betw_hex[12][3] = { { 0, 1, 8 }, @@ -259,6 +282,21 @@ namespace netgen onp = 6; break; } + case PRISM15: + { + betw = betw_prism15; + newel.SetType(PRISM15); + onp = 6; + break; + } + case PYRAMID: + case PYRAMID13: + { + betw = betw_pyramid; + newel.SetType(PYRAMID13); + onp = 5; + break; + } case HEX: case HEX20: { diff --git a/libsrc/meshing/topology.hpp b/libsrc/meshing/topology.hpp index 625b15d7..b4e646d2 100644 --- a/libsrc/meshing/topology.hpp +++ b/libsrc/meshing/topology.hpp @@ -212,10 +212,12 @@ inline short int MeshTopology :: GetNVertices (ELEMENT_TYPE et) return 4; case PYRAMID: + case PYRAMID13: return 5; case PRISM: case PRISM12: + case PRISM15: return 6; case HEX: @@ -257,10 +259,15 @@ inline short int MeshTopology :: GetNPoints (ELEMENT_TYPE et) case PYRAMID: return 5; + case PYRAMID13: + return 13; case PRISM: - case PRISM12: return 6; + case PRISM12: + return 12; + case PRISM15: + return 15; case HEX: return 8; @@ -277,7 +284,7 @@ inline short int MeshTopology :: GetNPoints (ELEMENT_TYPE et) inline short int MeshTopology :: GetNEdges (ELEMENT_TYPE et) { - __assume(et >= SEGMENT && et <= HEX); + __assume(et >= SEGMENT && et <= PYRAMID13); switch (et) { case SEGMENT: @@ -298,10 +305,12 @@ inline short int MeshTopology :: GetNEdges (ELEMENT_TYPE et) return 6; case PYRAMID: + case PYRAMID13: return 8; case PRISM: case PRISM12: + case PRISM15: return 9; case HEX: @@ -318,7 +327,7 @@ inline short int MeshTopology :: GetNEdges (ELEMENT_TYPE et) inline short int MeshTopology :: GetNFaces (ELEMENT_TYPE et) { - __assume(et >= SEGMENT && et <= HEX); + __assume(et >= SEGMENT && et <= PYRAMID13); switch (et) { case SEGMENT: @@ -339,10 +348,12 @@ inline short int MeshTopology :: GetNFaces (ELEMENT_TYPE et) return 4; case PYRAMID: + case PYRAMID13: return 5; case PRISM: case PRISM12: + case PRISM15: return 5; case HEX: @@ -443,10 +454,12 @@ const ELEMENT_EDGE * MeshTopology :: GetEdges1 (ELEMENT_TYPE et) return tet_edges; case PYRAMID: + case PYRAMID13: return pyramid_edges; case PRISM: case PRISM12: + case PRISM15: return prism_edges; case HEX: @@ -542,10 +555,12 @@ const ELEMENT_EDGE * MeshTopology :: GetEdges0 (ELEMENT_TYPE et) return tet_edges; case PYRAMID: + case PYRAMID13: return pyramid_edges; case PRISM: case PRISM12: + case PRISM15: return prism_edges; case HEX: @@ -627,9 +642,11 @@ inline const ELEMENT_FACE * MeshTopology :: GetFaces1 (ELEMENT_TYPE et) case PRISM: case PRISM12: + case PRISM15: return prism_faces; case PYRAMID: + case PYRAMID13: return pyramid_faces; case SEGMENT: @@ -710,9 +727,11 @@ inline const ELEMENT_FACE * MeshTopology :: GetFaces0 (ELEMENT_TYPE et) case PRISM: case PRISM12: + case PRISM15: return prism_faces; case PYRAMID: + case PYRAMID13: return pyramid_faces; case SEGMENT: diff --git a/libsrc/visualization/vsmesh.cpp b/libsrc/visualization/vsmesh.cpp index d758c626..a6c5ef4c 100644 --- a/libsrc/visualization/vsmesh.cpp +++ b/libsrc/visualization/vsmesh.cpp @@ -2682,7 +2682,7 @@ namespace netgen for (ElementIndex ei = 0; ei < mesh->GetNE(); ei++) { const Element & el = (*mesh)[ei]; - if (el.GetType() == PYRAMID && !el.IsDeleted()) + if ((el.GetType() == PYRAMID || el.GetType() == PYRAMID13) && !el.IsDeleted()) { int i = ei + 1; diff --git a/libsrc/visualization/vssolution.cpp b/libsrc/visualization/vssolution.cpp index 201f1fba..3b0b61cc 100644 --- a/libsrc/visualization/vssolution.cpp +++ b/libsrc/visualization/vssolution.cpp @@ -2857,6 +2857,7 @@ namespace netgen } case PRISM: case PRISM12: + case PRISM15: { lami[0] = (1-lam3) * (1-lam1-lam2); lami[1] = (1-lam3) * lam1; @@ -2907,6 +2908,7 @@ namespace netgen } case PRISM: case PRISM12: + case PRISM15: { lami[0] = (1-lam3) * (1-lam1-lam2); lami[1] = (1-lam3) * lam1; @@ -2918,6 +2920,7 @@ namespace netgen break; } case PYRAMID: + case PYRAMID13: { if (lam3 > 1-1e-5) { @@ -4019,7 +4022,7 @@ namespace netgen if(vispar.donotclipdomain > 0 && vispar.donotclipdomain == (*mesh)[ei].GetIndex()) continue; ELEMENT_TYPE type = (*mesh)[ei].GetType(); - if (type == HEX || type == PRISM || type == TET || type == TET10 || type == PYRAMID) + if (type == HEX || type == PRISM || type == TET || type == TET10 || type == PYRAMID || type == PYRAMID13 || type == PRISM15 || type == HEX20) { const Element & el = (*mesh)[ei]; @@ -4078,6 +4081,8 @@ namespace netgen switch (type) { case PRISM: + case PRISM12: + case PRISM15: if (ix+iy <= n) { ploc = Point<3> (double(ix) / n, double(iy) / n, double(iz) / n); @@ -4088,9 +4093,11 @@ namespace netgen compress[ii] = -1; break; case HEX: + case HEX20: ploc = Point<3> (double(ix) / n, double(iy) / n, double(iz) / n); break; case PYRAMID: + case PYRAMID13: ploc = Point<3> (double(ix) / n * (1-double(iz)/n), double(iy) / n * (1-double(iz)/n), double(iz)/n); @@ -4104,7 +4111,7 @@ namespace netgen locgrid[compress[ii]] = ploc; } - if (type != TET && type != TET10 && type != PRISM) cnt_valid = n3; + if (type != TET && type != TET10 && type != PRISM && type != PRISM12 && type != PRISM15) cnt_valid = n3; locgrid.SetSize(cnt_valid); diff --git a/python/read_gmsh.py b/python/read_gmsh.py index 0617732d..23a7ec52 100644 --- a/python/read_gmsh.py +++ b/python/read_gmsh.py @@ -1,9 +1,10 @@ from netgen.meshing import * - def ReadGmsh(filename): + if not filename.endswith(".msh"): + filename += ".msh" meshdim = 1 - with open(filename + ".msh", 'r') as f: + with open(filename, 'r') as f: while f.readline().split()[0] != "$Elements": pass nelem = int(f.readline()) @@ -16,15 +17,58 @@ def ReadGmsh(filename): meshdim = 3 break - f = open(filename + ".msh", 'r') + f = open(filename, 'r') mesh = Mesh(dim=meshdim) pointmap = {} facedescriptormap = {} - namemap = {} + namemap = { 0 : "default" } materialmap = {} bbcmap = {} + segm = 1 + trig = 2 + quad = 3 + tet = 4 + hex = 5 + prism = 6 + pyramid = 7 + segm3 = 8 # 2nd order line + trig6 = 9 # 2nd order trig + tet10 = 11 # 2nd order tet + point = 15 + quad8 = 16 # 2nd order quad + hex20 = 17 # 2nd order hex + prism15 = 18 # 2nd order prism + pyramid13 = 19 # 2nd order pyramid + segms = [segm, segm3] + trigs = [trig, trig6] + quads = [quad, quad8] + tets = [tet, tet10] + hexes = [hex, hex20] + prisms = [prism, prism15] + pyramids = [pyramid, pyramid13] + elem0d = [point] + elem1d = segms + elem2d = trigs + quads + elem3d = tets + hexes + prisms + pyramids + + num_nodes_map = { segm : 2, + trig : 3, + quad : 4, + tet : 4, + hex : 8, + prism : 6, + pyramid : 5, + segm3 : 3, + trig6 : 6, + tet10 : 10, + point : 1, + quad8 : 8, + hex20 : 20, + prism15 : 18, + pyramid13 : 19 } + while True: line = f.readline() if line == "": @@ -57,29 +101,14 @@ def ReadGmsh(filename): # the first tag is the physical group nr, the second tag is the group nr of the dim tags = [int(line[3 + k]) for k in range(numtags)] - if elmtype == 1: # 2-node line - num_nodes = 2 - elif elmtype == 2: # 3-node trig - num_nodes = 3 - elif elmtype == 3: # 4-node quad - num_nodes = 4 - elif elmtype == 4: # 4-node tet - num_nodes = 4 - elif elmtype == 5: # 8-node hex - num_nodes = 8 - elif elmtype == 6: # 6-node prism - num_nodes = 6 - elif elmtype == 7: # 5-node pyramid - num_nodes = 5 - elif elmtype == 15: # 1-node point - num_nodes = 1 - else: + if elmtype not in num_nodes_map: raise Exception("element type", elmtype, "not implemented") + num_nodes = num_nodes_map[elmtype] nodenums = line[3 + numtags:3 + numtags + num_nodes] nodenums2 = [pointmap[int(nn)] for nn in nodenums] - if elmtype == 1: + if elmtype in elem1d: if meshdim == 3: if tags[1] in bbcmap: index = bbcmap[tags[1]] @@ -117,7 +146,7 @@ def ReadGmsh(filename): mesh.Add(Element1D(index=index, vertices=nodenums2)) - if elmtype in [2, 3]: # 2d elements + if elmtype in elem2d: # 2d elements if meshdim == 3: if tags[1] in facedescriptormap.keys(): index = facedescriptormap[tags[1]] @@ -142,9 +171,17 @@ def ReadGmsh(filename): mesh.SetMaterial(index, "surf" + str(tags[1])) materialmap[tags[1]] = index - mesh.Add(Element2D(index, nodenums2)) + if elmtype in trigs: + ordering = [i for i in range(3)] + if elmtype == trig6: + ordering += [4,5,3] + if elmtype in quads: + ordering = [i for i in range(4)] + if elmtype == quad8: + ordering += [4, 6, 7, 5] + mesh.Add(Element2D(index, [nodenums2[i] for i in ordering])) - if elmtype in [4, 5, 6, 7]: # volume elements + if elmtype in elem3d: # volume elements if tags[1] in materialmap: index = materialmap[tags[1]] else: @@ -157,9 +194,22 @@ def ReadGmsh(filename): nodenums2 = [pointmap[int(nn)] for nn in nodenums] - if elmtype == 4: - mesh.Add(Element3D(index, [nodenums2[0], nodenums2[1], nodenums2[3], nodenums2[2]])) - elif elmtype in [5, 6, 7]: - mesh.Add(Element3D(index, nodenums2)) + if elmtype in tets: + ordering = [0,1,2,3] + if elmtype == tet10: + ordering += [4,6,7,5,9,8] + elif elmtype in hexes: + ordering = [0,1,5,4,3,2,6,7] + if elmtype == hex20: + ordering += [8,16,10,12,13,19,15,14,9,11,18,17] + elif elmtype in prisms: + ordering = [0,2,1,3,5,4] + if elmtype == prism15: + ordering += [7,6,9,8,11,10,13,12,14] + elif elmtype in pyramids: + ordering = [3,2,1,0,4] + if elmtype == pyramid13: + ordering += [10,5,6,8,12,11,9,7] + mesh.Add(Element3D(index, [nodenums2[i] for i in ordering])) return mesh