ForkMaster/netgen
1
0
Fork 0
mirror of https://github.com/NGSolve/netgen.git synced 2024-12-25 05:20:34 +05:00
Code Issues Packages Projects Releases Wiki Activity

Merge branch 'gmsh_2nd_order_elements' into 'master'

Browse Source 
read tet10, pyramid13, prism15 and hex20 from gmsh

See merge request jschoeberl/netgen!127
This commit is contained in:
Joachim Schöberl 2019-02-06 23:29:21 +00:00
commit 01fe04fe61
13 changed files with 396 additions and 89 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
libsrc/general/array.hpp
View File

@ -243,9 +243,10 @@ namespace netgen
}
explicit Array(size_t asize)
: FlatArray<T, BASE, TIND> (asize, new T[asize])
: FlatArray<T, BASE, TIND> (asize, asize ? new T[asize] : nullptr)
{
allocsize = asize;
if(asize)
ownmem = 1;
}

4
libsrc/include/nginterface.h
View File

@ -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
};

2
libsrc/include/nginterface_v2.hpp
View File

@ -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
};

127
libsrc/meshing/curvedelems.cpp
View File

@ -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<double>)) 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<double>)) return;

6
libsrc/meshing/meshclass.cpp
View File

@ -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 &&

89
libsrc/meshing/meshtype.cpp
View File

@ -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");
}

9
libsrc/meshing/meshtype.hpp
View File

@ -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 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 PYRAMID: case PYRAMID13: return 5;
case PRISM:
case PRISM15:
case PRISM12: return 5;
case HEX: case HEX20:
return 6;

49
libsrc/meshing/python_mesh.cpp
View File

@ -256,44 +256,19 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m)
py::class_<Element>(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<PointIndex>(vertices[i])();
std::map<int, ELEMENT_TYPE> 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<np; i++)
(*newel)[i] = py::cast<PointIndex>(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<PointIndex>(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<PointIndex>(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<PointIndex>(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<PointIndex>(vertices[i])();
newel->SetIndex(index);
}
else
throw NgException ("cannot create element");
return newel;
}),
py::arg("index")=1,py::arg("vertices"),

38
libsrc/meshing/secondorder.cpp
View File

@ -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:
{

25
libsrc/meshing/topology.hpp
View File

@ -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:

2
libsrc/visualization/vsmesh.cpp
View File

@ -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;

11
libsrc/visualization/vssolution.cpp
View File

@ -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);

108
python/read_gmsh.py
View File

@ -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