- Thank you received: 24
Forum Message
We have moved the forum to https://forum.ngsolve.org . This is an archived version of the topics until 05/05/23. All the topics were moved to the new forum and conversations can be continued there. This forum is just kept as legacy to not invalidate old links. If you want to continue a conversation just look for the topic in the new forum.
Unstructed mesh with structured subdomain
4 years 6 months ago - 4 years 6 months ago #1975
by christopher
Replied by christopher on topic Unstructed mesh with structured subdomain
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
Hope this helps
Best
Christopher
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
Code:
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)
Hope this helps
Best
Christopher
Attachments:
Time to create page: 0.124 seconds
