Unstructed mesh with structured subdomain

Main Forum Old-Forum - do not post here
1

Hi everyone,

is it possible to create a mesh, where an inner subdomain contains a structured mesh (e.g. with meshes.MakeStructured2DMesh() ) which in turn is inside a larger domain with a coarser, unstructured mesh?

So something like this:

geo = SplineGeometry()
geo.AddRectangle([0,0],[1,1])
geo.AddRectangle([0.33,0.33],[0.66,0.66],leftdomain=2, rightdomain=1)
geo.SetDomainMaxH(1,0.1)
geo.SetDomainMaxH(2,0.02)
mesh = Mesh(geo.GenerateMesh())

but with a structured mesh in domain 2.

Best wishes,
Henry

2

It is, but not straight forward:

create a mesh for the outside, where the boundaries of the top and bottom edge of the enclosed subdomain have copied edges. Then you need to fill the interior with a structured mesh. The following code uses some knowledge of how edges and so on are treated internally and moving code lines may result in unexpected behaviour :wink:

from netgen.geom2d import *
from netgen.meshing import *

geo = SplineGeometry()
geo.AddRectangle([0,0],[1,1])

sd = ((0.33,0.33),(0.66,0.66))
pts = [geo.AppendPoint(*p) for p in [sd[0], (sd[1][0], sd[0][1]), sd[1], (sd[0][1], sd[1][0])]]

index = {}
index["bot"] = geo.Append( ["line", pts[0], pts[1]], leftdomain=0, rightdomain=1, maxh=0.01)
index["top"] = geo.Append( ["line", pts[3], pts[2]], copy=index["bot"], maxh=0.01)
index["right"] = geo.Append( ["line", pts[1], pts[2]], leftdomain=0, rightdomain=1)
index["left"] = geo.Append( ["line", pts[0], pts[3]], copy=index["right"])

geo.SetDomainMaxH(1,0.1)

mesh = geo.GenerateMesh()

# inverse map with index shift
inv_map = {}
for key,val in index.items():
    inv_map[val+1] = key

segments = {key : [] for key in index}

for el in mesh.Elements1D():
    if el.index in inv_map:
        segments[inv_map[el.index]].append(el.vertices)

pnts = {}
for key, val in segments.items():
    pnts[key] = [mesh[el[0]] for el in val]
    pnts[key].append(mesh[val[-1][1]])

rows = [[seg[0] for seg in segments["bot"]]]
rows[0].append(segments["bot"][-1][1])
for i in range(1,len(segments["left"])):
    new_row = [segments["left"][i][0]]
    for j in range(1, len(segments["bot"])):
        new_row.append(mesh.Add(MeshPoint(Pnt(pnts["bot"][j][0], pnts["left"][i][1], 0))))
    new_row.append(segments["right"][i][0])
    rows.append(new_row)
rows.append([seg[0] for seg in segments["top"]])
rows[-1].append(segments["top"][-1][1])

mesh.Add(FaceDescriptor(surfnr = 2, bc=2))

for i in range(len(rows)-1):
    for j in range(len(rows[0])-1):
        el = Element2D(index=2, vertices=[rows[i][j], rows[i][j+1], rows[i+1][j+1], rows[i+1][j]])
        mesh.Add(el)

from ngsolve import *
msh = Mesh(mesh)
Draw(msh)

800×799 56.5 KB

Hope this helps :slight_smile:

Best
Christopher

3

Hi Christopher,

thank you very much for your help! I have changed the 2D element loop to get triangular elements, but is working perfectly!

Best wishes,
Henry