Set Material in custom 2d mesh

Netgen
1

Hello, everyone,

I am trying to construct a custom mesh in 2d on a circle shaped domain with special properties. I am nearly finished, but I noticed that I am not able to add different materials to my mesh and even though the visualization shows 2 different domains mesh.GetNDomains() results in 1. I have already tried to set all domain IDs correctly and also added FaceDescriptors for the boundaries so I don’t know where the problem could be. Here is a shortened version of my code that showcases my problem:

import netgen.meshing as ngm
import netgen.geom2d as geom2d
from ngsolve import *

###desired geometry
radius=0.9
center=(0,0)
geo = geom2d.SplineGeometry()
geo.AddCircle(c=center,r=1,leftdomain=3,bc='o')
geo.AddCircle(c=center, r=radius, leftdomain=4,rightdomain=3)
###geometry inner circle
geo_ref= geom2d.SplineGeometry()
geo_ref.AddCircle(c=center, r=radius)
mesh_ref=geo_ref.GenerateMesh(maxh=0.1)
###Generate custom mesh
mesh=ngm.Mesh()
mesh.SetGeometry(geo)
pmap = { }
pmap1= { }

fd_inter = mesh.Add(ngm.FaceDescriptor(bc=2))
fd_insideb = mesh.Add(ngm.FaceDescriptor(bc=1))
f_in=1
f_out=2
#Build custom mesh
points=mesh_ref.Points()
for e in mesh_ref.Elements1D():
    for v in e.vertices:
        if (v not in pmap):
            pmap[v] = mesh.Add(mesh_ref[v])
            point=points[v]
            alpha=1+0.1/radius
            new_point=(alpha*point[0]+(1-alpha)*center[0],alpha*point[1]+(1-alpha)*center[1])
            pmap1[v]=mesh.Add(ngm.MeshPoint(ngm.Pnt(new_point[0], new_point[1],0)))
    nr=e.edgenr
    verts=[pmap[v] for v in e.vertices]
    mesh.Add(ngm.Element1D(verts,index=fd_insideb,edgenr=nr+4))
    verts1=[pmap1[v] for v in e.vertices]
    mesh.Add(ngm.Element1D(verts1,index=fd_inter,edgenr=nr))
    verts.reverse()
    mesh.Add(ngm.Element2D(f_out, verts+verts1))

for f in mesh_ref.Elements2D():
    for v in f.vertices:
        if (v not in pmap):
            pmap[v] = mesh.Add(mesh_ref[v])
    mesh.Add(ngm.Element2D(f_in,[pmap[v] for v in f.vertices]))

#set materials (seems to be doing nothing)
mesh.SetMaterial(1,"Mat1")
mesh.SetMaterial(2,"Mat2")

ngmesh=Mesh(mesh)
print(ngmesh.GetMaterials())
print(mesh.GetNDomains())

I appreciate your help.

Sincerely,
Florian

2

Hi,
the GetNDomain function in netgen gives the max number that is in domin or domout of facedescriptors and in 2d meaningless.
To have the right number in ngsolve in ngmesh.GetMaterials you have to specify the correct dimension when creating your netgen mesh:

mesh = Mesh(dim=2)

Note, if you want to create a circle with a thin layer of anisotropic quad mesh, there is also the option to do it in occ with close surface identification:

from netgen.occ import *
from ngsolve import *
from netgen.meshing import IdentificationType

c1 = Circle((0,0),1).Face()
c2 = Circle((0,0),1.1).Face()

trf = gp_Trsf.Scale((0,0,0), 1.1)

# 3 is close surface identification
c1.edges.Identify(c2.edges, "cs", 3, trf)

c2 -= c1

shape = Glue([c1, c2])

geo = OCCGeometry(shape)
mesh = Mesh(geo.GenerateMesh(maxh=0.1))
Draw(mesh)

Best
Christopher

3

Hi Christopher,

thank you very much. That solves my problem and now everything works as expected!

Best,
Florian