From 05106cf58bc1c7783957fbe3a24198f74e26760a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Joachim=20Sch=C3=B6berl?= Date: Tue, 5 Apr 2016 17:15:31 +0200 Subject: [PATCH] hex-filling of thin domains (first steps ...) --- libsrc/csg/python_csg.cpp | 11 +++++++++-- libsrc/meshing/CMakeLists.txt | 2 +- libsrc/meshing/curvedelems.cpp | 34 ++++++++++++++++++++++++++++++++++ libsrc/meshing/meshfunc.cpp | 10 ++++++++-- libsrc/meshing/meshtype.cpp | 20 +++++++++++++++++++- libsrc/meshing/meshtype.hpp | 1 + libsrc/meshing/parser3.cpp | 15 +++++++++++++++ 7 files changed, 87 insertions(+), 6 deletions(-) diff --git a/libsrc/csg/python_csg.cpp b/libsrc/csg/python_csg.cpp index 7d06a5d7..45b6f6dd 100644 --- a/libsrc/csg/python_csg.cpp +++ b/libsrc/csg/python_csg.cpp @@ -502,8 +502,15 @@ DLL_HEADER void ExportCSG() SetGlobalMesh (dummy); dummy->SetGeometry(geo); ng_geometry = geo; - geo->FindIdenticSurfaces(1e-8 * geo->MaxSize()); - geo->GenerateMesh (dummy, param, 0, 6); + geo->FindIdenticSurfaces(1e-8 * geo->MaxSize()); + try + { + geo->GenerateMesh (dummy, param, 0, 6); + } + catch (NgException ex) + { + cout << "Caught NgException: " << ex.What() << endl; + } return dummy; })) ; diff --git a/libsrc/meshing/CMakeLists.txt b/libsrc/meshing/CMakeLists.txt index 10324cce..e1610332 100644 --- a/libsrc/meshing/CMakeLists.txt +++ b/libsrc/meshing/CMakeLists.txt @@ -20,7 +20,7 @@ add_library(mesh ${NG_LIB_TYPE} smoothing2.cpp smoothing3.cpp specials.cpp tetrarls.cpp topology.cpp triarls.cpp validate.cpp bcfunctions.cpp parallelmesh.cpp paralleltop.cpp paralleltop.hpp basegeom.cpp - python_mesh.cpp + python_mesh.cpp hexarls.cpp ${mesh_object_libs} ) diff --git a/libsrc/meshing/curvedelems.cpp b/libsrc/meshing/curvedelems.cpp index 5e02fbfe..40ebbe8d 100644 --- a/libsrc/meshing/curvedelems.cpp +++ b/libsrc/meshing/curvedelems.cpp @@ -2590,6 +2590,40 @@ namespace netgen shapes[5] = x *(1-y)*(z); shapes[6] = x * y *(z); shapes[7] = (1-x)* y *(z); + + if (info.order == 1) return; + + double mu[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), + }; + + int ii = 8; + const ELEMENT_EDGE * edges = MeshTopology::GetEdges1 (HEX); + + for (int i = 0; i < 8; i++) + { + int eorder = edgeorder[info.edgenrs[i]]; + if (eorder >= 2) + { + int vi1 = edges[i][0]-1, vi2 = edges[i][1]-1; + if (el[vi1] > el[vi2]) swap (vi1, vi2); + + CalcEdgeShape (eorder, mu[vi1]-mu[vi2], &shapes(ii)); + double lame = shapes(vi1)+shapes(vi2); + for (int j = 0; j < order-1; j++) + shapes(ii+j) *= lame; + ii += eorder-1; + } + } + + break; } diff --git a/libsrc/meshing/meshfunc.cpp b/libsrc/meshing/meshfunc.cpp index f419f718..ab6de3b2 100644 --- a/libsrc/meshing/meshfunc.cpp +++ b/libsrc/meshing/meshfunc.cpp @@ -8,6 +8,7 @@ namespace netgen extern const char * prismrules2[]; extern const char * pyramidrules[]; extern const char * pyramidrules2[]; + extern const char * hexrules[]; // extern double teterrpow; @@ -94,8 +95,9 @@ namespace netgen } domain_bbox.Increase (0.01 * domain_bbox.Diam()); - - for (int qstep = 1; qstep <= 3; qstep++) + + for (int qstep = 1; qstep <= 3; qstep++) + // for (int qstep = 0; qstep <= 3; qstep++) // for hex-filling { // cout << "openquads = " << mesh3d.HasOpenQuads() << endl; if (mesh3d.HasOpenQuads()) @@ -105,6 +107,10 @@ namespace netgen const char ** rulep = NULL; switch (qstep) { + case 0: + // rulefile = "/Users/joachim/gitlab/netgen/rules/hexa.rls"; + rulep = hexrules; + break; case 1: rulefile += "/rules/prisms2.rls"; rulep = prismrules2; diff --git a/libsrc/meshing/meshtype.cpp b/libsrc/meshing/meshtype.cpp index 0f59bf8f..b6eb599f 100644 --- a/libsrc/meshing/meshtype.cpp +++ b/libsrc/meshing/meshtype.cpp @@ -1165,7 +1165,7 @@ namespace netgen { 3, 2, 3, 5, 0 }, { 3, 3, 4, 5, 0 }, { 3, 4, 1, 5, 0 } }; - + static const int prismfaces[][5] = { { 3, 1, 3, 2, 0 }, @@ -1174,6 +1174,17 @@ namespace netgen { 4, 2, 3, 6, 5 }, { 4, 3, 1, 4, 6 } }; + + static const int hexfaces[][5] = + { + { 4, 4, 3, 2, 1 }, + { 4, 3, 7, 6, 2 }, + { 4, 7, 8, 5, 6 }, + { 4, 8, 4, 1, 5 }, + { 4, 1, 2, 6, 5 }, + { 4, 3, 4, 8, 7 } + }; + switch (np) { @@ -1209,6 +1220,13 @@ namespace netgen face.PNum(j) = PNum(prismfaces[i-1][j]); break; } + case 8: + { + face.SetType(QUAD); + for (int j = 1; j <= 4; j++) + face.PNum(j) = PNum(hexfaces[i-1][j]); + break; + } } } diff --git a/libsrc/meshing/meshtype.hpp b/libsrc/meshing/meshtype.hpp index 0bf701f0..37b08460 100644 --- a/libsrc/meshing/meshtype.hpp +++ b/libsrc/meshing/meshtype.hpp @@ -679,6 +679,7 @@ namespace netgen case PYRAMID: return 5; case PRISM: case PRISM12: return 5; + case HEX: return 6; default: #ifdef DEBUG PrintSysError ("element3d::GetNFaces not implemented for typ", typ) diff --git a/libsrc/meshing/parser3.cpp b/libsrc/meshing/parser3.cpp index a661fced..7b69212d 100644 --- a/libsrc/meshing/parser3.cpp +++ b/libsrc/meshing/parser3.cpp @@ -530,6 +530,21 @@ void vnetrule :: LoadRule (istream & ist) ist >> elements.Last().PNum(6); ist >> ch; // ',' } + + if (ch == COMMASIGN) + { + // elements.Last().SetNP(6); + elements.Last().SetType(HEX); + ist >> elements.Last().PNum(7); + ist >> ch; // ',' + } + if (ch == COMMASIGN) + { + // elements.Last().SetNP(6); + elements.Last().SetType(HEX); + ist >> elements.Last().PNum(8); + ist >> ch; // ',' + } /* orientations.Append (fourint());