From 5bea3bb6122a8526dec8e89bbb6799c8b2ce80ea Mon Sep 17 00:00:00 2001 From: Matthias Hochsteger Date: Mon, 8 Jun 2020 10:39:31 +0200 Subject: [PATCH 1/4] Implement and export SplineGeometry2d::SetDomainTensorMeshing --- libsrc/geom2d/geometry2d.hpp | 15 +++++++++++++-- libsrc/geom2d/python_geom2d.cpp | 1 + 2 files changed, 14 insertions(+), 2 deletions(-) diff --git a/libsrc/geom2d/geometry2d.hpp b/libsrc/geom2d/geometry2d.hpp index ff4c459c..f5841b19 100644 --- a/libsrc/geom2d/geometry2d.hpp +++ b/libsrc/geom2d/geometry2d.hpp @@ -133,7 +133,7 @@ namespace netgen NgArray materials; NgArray maxh; NgArray quadmeshing; - NgArray tensormeshing; + Array tensormeshing; NgArray layer; NgArray bcnames; double elto0 = 1.0; @@ -216,9 +216,20 @@ namespace netgen } bool GetDomainTensorMeshing ( int domnr ) { - if ( tensormeshing.Size() ) return tensormeshing[domnr-1]; + if ( tensormeshing.Size()>=domnr ) return tensormeshing[domnr-1]; else return false; } + void SetDomainTensorMeshing ( int domnr, bool tm ) + { + if ( tensormeshing.Size()(), meshingparameter_description.c_str()) + .def("_SetDomainTensorMeshing", &SplineGeometry2d::SetDomainTensorMeshing) ; } From d1c7a16d63ba82071957aa21cae74d999c8bb6c5 Mon Sep 17 00:00:00 2001 From: Matthias Hochsteger Date: Mon, 8 Jun 2020 10:41:22 +0200 Subject: [PATCH 2/4] Do linear interpolation of corresponding edge points in SplineGeometry tensor mesh generation better results for curved domains --- libsrc/geom2d/genmesh2d.cpp | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/libsrc/geom2d/genmesh2d.cpp b/libsrc/geom2d/genmesh2d.cpp index 09ef0260..89eaa098 100644 --- a/libsrc/geom2d/genmesh2d.cpp +++ b/libsrc/geom2d/genmesh2d.cpp @@ -527,11 +527,16 @@ namespace netgen for (PointIndex pix = nextpi[c1], ix = 0; pix != c2; pix = nextpi[pix], ix++) + { + Point<3> px = (*mesh)[pix]; for (PointIndex piy = nextpi[c2], iy = 0; piy != c3; piy = nextpi[piy], iy++) { - Point<3> p = (*mesh)[pix] + ( (*mesh)[piy] - (*mesh)[c2] ); - pts[(nex+1)*(iy+1) + ix+1] = mesh -> AddPoint (p , 1, FIXEDPOINT); + double lam = Dist((*mesh)[piy],(*mesh)[c2]) / Dist((*mesh)[c3],(*mesh)[c2]); + auto pix1 = pts[(nex+1)*ney+ix+1]; + auto pnew = px + lam*((*mesh)[pix1]-px); + pts[(nex+1)*(iy+1) + ix+1] = mesh -> AddPoint (pnew, 1, FIXEDPOINT); } + } for (int i = 0; i < ney; i++) for (int j = 0; j < nex; j++) From d08e2daa060da23965f68bfc4b083ae8dd271539 Mon Sep 17 00:00:00 2001 From: Matthias Hochsteger Date: Mon, 8 Jun 2020 10:42:26 +0200 Subject: [PATCH 3/4] Do edge swapping for faces individually with tensor product meshes If the mesh contains quads, the edge swapping algorithm switches to generic improve, which introduces quads everywhere. This is not intended if one domain contains a tensor product mesh. Thus, call the optimizer for each face if mesh contains quads but mp.quad is not set. --- libsrc/meshing/meshfunc2d.cpp | 45 ++++++++++++++++++++++++++++++++--- 1 file changed, 42 insertions(+), 3 deletions(-) diff --git a/libsrc/meshing/meshfunc2d.cpp b/libsrc/meshing/meshfunc2d.cpp index 88955a36..84f5d399 100644 --- a/libsrc/meshing/meshfunc2d.cpp +++ b/libsrc/meshing/meshfunc2d.cpp @@ -20,6 +20,19 @@ namespace netgen } mesh.Compress(); + bool optimize_swap_separate_faces = false; + if(!mp.quad) + { + bool mixed = false; + ParallelFor( Range(mesh.GetNSE()), [&] (auto i) NETGEN_LAMBDA_INLINE + { + if (mesh[SurfaceElementIndex(i)].GetNP() != 3) + mixed = true; + }); + if(mixed) + optimize_swap_separate_faces = true; + } + const char * optstr = mp.optimize2d.c_str(); int optsteps = mp.optsteps2d; @@ -31,16 +44,42 @@ namespace netgen { case 's': { // topological swap - MeshOptimize2d meshopt(mesh); + MeshOptimize2d meshopt(mesh); meshopt.SetMetricWeight (mp.elsizeweight); - meshopt.EdgeSwapping (0); + + if(optimize_swap_separate_faces) + { + for(auto i : Range(1, mesh.GetNFD()+1)) + { + meshopt.SetFaceIndex(i); + meshopt.EdgeSwapping (0); + } + } + else + { + meshopt.SetFaceIndex(0); + meshopt.EdgeSwapping (0); + } break; } case 'S': { // metric swap MeshOptimize2d meshopt(mesh); meshopt.SetMetricWeight (mp.elsizeweight); - meshopt.EdgeSwapping (1); + + if(optimize_swap_separate_faces) + { + for(auto i : Range(1, mesh.GetNFD()+1)) + { + meshopt.SetFaceIndex(i); + meshopt.EdgeSwapping (1); + } + } + else + { + meshopt.SetFaceIndex(0); + meshopt.EdgeSwapping (1); + } break; } case 'm': From c0f50820cb72d6df5d6411490884dee3a800c834 Mon Sep 17 00:00:00 2001 From: Matthias Hochsteger Date: Mon, 8 Jun 2020 10:54:02 +0200 Subject: [PATCH 4/4] Add test for SetDomainTensorMeshing --- .../test_splinegeo_tensordomainmeshing.py | 22 +++++++++++++++++++ 1 file changed, 22 insertions(+) create mode 100644 tests/pytest/test_splinegeo_tensordomainmeshing.py diff --git a/tests/pytest/test_splinegeo_tensordomainmeshing.py b/tests/pytest/test_splinegeo_tensordomainmeshing.py new file mode 100644 index 00000000..0c00ae7d --- /dev/null +++ b/tests/pytest/test_splinegeo_tensordomainmeshing.py @@ -0,0 +1,22 @@ +from netgen.geom2d import * + +def test_tensordomainmeshing(): + geo = SplineGeometry() + w = 10 + h = 0.01 + + p = [ (0, 0), (w, 0), (w, h), (0, h) ] + p = [geo.AppendPoint(*px) for px in p] + + l0 = geo.Append ( ["line", p[0], p[1]], leftdomain=1, rightdomain=0 ) + l1 = geo.Append ( ["line", p[1], p[2]], leftdomain=1, rightdomain=0) + geo.Append ( ["line", p[3], p[2]], leftdomain=0, rightdomain=1, copy=l0 ) + geo.Append ( ["line", p[0], p[3]], leftdomain=0, rightdomain=1, copy=l1 ) + + geo._SetDomainTensorMeshing(1, True) + + mesh = geo.GenerateMesh(maxh=1) + + for el in mesh.Elements2D(): + print(el.vertices) + assert len(el.vertices) == 4