From ec0a2a5ae890cd990bfc99040108ecfe1cea672c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Joachim=20Sch=C3=B6berl?= Date: Sat, 19 Jan 2019 14:26:17 +0100 Subject: [PATCH] quad8 and hex20 geometric elements --- libsrc/meshing/curvedelems.cpp | 87 ++++++++++++++++++++++++++++++++-- libsrc/meshing/meshtype.cpp | 3 ++ libsrc/meshing/meshtype.hpp | 8 ++-- libsrc/meshing/secondorder.cpp | 24 ++++++++++ libsrc/meshing/topology.cpp | 13 ++--- 5 files changed, 121 insertions(+), 14 deletions(-) diff --git a/libsrc/meshing/curvedelems.cpp b/libsrc/meshing/curvedelems.cpp index 56c276d9..5e3fb052 100644 --- a/libsrc/meshing/curvedelems.cpp +++ b/libsrc/meshing/curvedelems.cpp @@ -2885,6 +2885,46 @@ namespace netgen break; + } + + case HEX20: + { + shapes = 0.0; + T x = xi(0); + T y = xi(1); + T z = xi(2); + + shapes[0] = (1-x)*(1-y)*(1-z); + shapes[1] = x *(1-y)*(1-z); + shapes[2] = x * y *(1-z); + shapes[3] = (1-x)* y *(1-z); + shapes[4] = (1-x)*(1-y)*(z); + shapes[5] = x *(1-y)*(z); + shapes[6] = x * y *(z); + shapes[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 = shapes[e[i][0]]+shapes[e[i][1]]; + T xi = sigma[e[i][1]]-sigma[e[i][0]]; + shapes[8+i] = (1-xi*xi)*lame; + } + for (int i = 0; i < 12; i++) + { + shapes[e[i][0]] -= 0.5 * shapes[8+i]; + shapes[e[i][1]] -= 0.5 * shapes[8+i]; + } + break; } default: @@ -3489,12 +3529,49 @@ namespace netgen *testout << "quad, num dshape = " << endl << dshapes << endl; */ break; - - - - break; } - + case HEX20: + { + 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[20]; + + ad[0] = (1-x)*(1-y)*(1-z); + ad[1] = x *(1-y)*(1-z); + ad[2] = x * y *(1-z); + ad[3] = (1-x)* y *(1-z); + ad[4] = (1-x)*(1-y)*(z); + ad[5] = x *(1-y)*(z); + ad[6] = x * y *(z); + ad[7] = (1-x)* y *(z); + + AutoDiff<3,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++) + { + auto lame = ad[e[i][0]]+ad[e[i][1]]; + auto xi = sigma[e[i][1]]-sigma[e[i][0]]; + ad[8+i] = (1-xi*xi)*lame; + } + for (int i = 0; i < 12; i++) + { + ad[e[i][0]] -= 0.5 * ad[8+i]; + ad[e[i][1]] -= 0.5 * ad[8+i]; + } + for (int i = 0; i < 20; i++) + for (int j = 0; j < 3; j++) + dshapes(i,j) = ad[i].DValue(j); + break; + } default: throw NgException("CurvedElements::CalcDShape 3d, element type not handled"); } diff --git a/libsrc/meshing/meshtype.cpp b/libsrc/meshing/meshtype.cpp index baed182e..494a3879 100644 --- a/libsrc/meshing/meshtype.cpp +++ b/libsrc/meshing/meshtype.cpp @@ -1028,6 +1028,7 @@ namespace netgen case 6: typ = PRISM; break; case 8: typ = HEX; break; case 10: typ = TET10; break; + case 20: typ = HEX20; break; default: cerr << "Element::Element: unknown element with " << np << " points" << endl; } orderx = ordery = orderz = 1; @@ -1108,6 +1109,7 @@ namespace netgen case 6: typ = PRISM; break; case 8: typ = HEX; break; case 10: typ = TET10; break; + case 20: typ = HEX20; break; // default: break; cerr << "Element::SetNP unknown element with " << np << " points" << endl; @@ -1127,6 +1129,7 @@ namespace netgen case HEX: np = 8; break; case TET10: np = 10; break; case PRISM12: np = 12; break; + case HEX20: np = 20; break; default: break; cerr << "Element::SetType unknown type " << int(typ) << endl; diff --git a/libsrc/meshing/meshtype.hpp b/libsrc/meshing/meshtype.hpp index a894d090..1b31d597 100644 --- a/libsrc/meshing/meshtype.hpp +++ b/libsrc/meshing/meshtype.hpp @@ -647,7 +647,7 @@ namespace netgen /// ELEMENT_TYPE typ; /// number of points (4..tet, 5..pyramid, 6..prism, 8..hex, 10..quad tet, 12..quad prism) - int np:5; + int np:6; /// class flagstruct { public: @@ -708,7 +708,7 @@ namespace netgen /// uint8_t GetNV() const { - __assume(typ >= TET && typ <= HEX); + __assume(typ >= TET && typ <= HEX20); switch (typ) { case TET: @@ -720,6 +720,7 @@ namespace netgen case PYRAMID: return 5; case HEX: + case HEX20: return 8; default: // not a 3D element #ifdef DEBUG @@ -797,7 +798,8 @@ namespace netgen case PYRAMID: return 5; case PRISM: case PRISM12: return 5; - case HEX: return 6; + case HEX: case HEX20: + return 6; default: #ifdef DEBUG PrintSysError ("element3d::GetNFaces not implemented for typ", typ) diff --git a/libsrc/meshing/secondorder.cpp b/libsrc/meshing/secondorder.cpp index 861f99ed..6241fd4e 100644 --- a/libsrc/meshing/secondorder.cpp +++ b/libsrc/meshing/secondorder.cpp @@ -223,6 +223,22 @@ namespace netgen { 3, 4, 10 }, { 4, 5, 11 }, }; + static int betw_hex[12][3] = + { + { 0, 1, 8 }, + { 2, 3, 9 }, + { 3, 0, 10 }, + { 1, 2, 11 }, + { 4, 5, 12 }, + { 6, 7, 13 }, + { 7, 4, 14 }, + { 5, 6, 15 }, + { 0, 4, 16 }, + { 1, 5, 17 }, + { 2, 6, 18 }, + { 3, 7, 19 }, + }; + int (*betw)[3] = NULL; switch (el.GetType()) @@ -243,6 +259,14 @@ namespace netgen onp = 6; break; } + case HEX: + case HEX20: + { + betw = betw_hex; + newel.SetType (HEX20); + onp = 8; + break; + } default: PrintSysError ("MakeSecondOrder, illegal vol type ", int(el.GetType())); } diff --git a/libsrc/meshing/topology.cpp b/libsrc/meshing/topology.cpp index cddc731f..04962a57 100644 --- a/libsrc/meshing/topology.cpp +++ b/libsrc/meshing/topology.cpp @@ -543,7 +543,7 @@ namespace netgen cnt = 0; ParallelForRange - (tm, mesh->Points().Size(), + (tm, mesh->GetNV(), // Points().Size(), [&] (size_t begin, size_t end) { INDEX_CLOSED_HASHTABLE v2eht(2*max_edge_on_vertex+10); @@ -581,7 +581,7 @@ namespace netgen // accumulate number of edges int ned = edge2vert.Size(); - for (auto v : mesh->Points().Range()) + for (size_t v = 0; v < mesh->GetNV(); v++) { auto hv = cnt[v]; cnt[v] = ned; @@ -595,7 +595,7 @@ namespace netgen // for (PointIndex v = PointIndex::BASE; v < nv+PointIndex::BASE; v++) ParallelForRange - (tm, mesh->Points().Size(), + (tm, mesh->GetNV(), // Points().Size(), [&] (size_t begin, size_t end) { INDEX_CLOSED_HASHTABLE v2eht(2*max_edge_on_vertex+10); @@ -719,7 +719,7 @@ namespace netgen // for (auto v : mesh.Points().Range()) NgProfiler::StartTimer (timer2b1); ParallelForRange - (tm, mesh->Points().Size(), + (tm, mesh->GetNV(), // Points().Size(), [&] (size_t begin, size_t end) { INDEX_3_CLOSED_HASHTABLE vert2face(2*max_face_on_vertex+10); @@ -754,7 +754,8 @@ namespace netgen // accumulate number of faces int nfa = oldnfa; - for (auto v : mesh->Points().Range()) + // for (auto v : Range(mesh->GetNV())) // Points().Range()) + for (size_t v = 0; v < mesh->GetNV(); v++) { auto hv = cnt[v]; cnt[v] = nfa; @@ -765,7 +766,7 @@ namespace netgen // for (auto v : mesh.Points().Range()) ParallelForRange - (tm, mesh->Points().Size(), + (tm, mesh->GetNV(), // Points().Size(), [&] (size_t begin, size_t end) { INDEX_3_CLOSED_HASHTABLE vert2face(2*max_face_on_vertex+10);