Changing boundary names via bcmod argument of CSGeometry::Add

Hello,

In my geometry I have a cylinder embedded in a cube. I also have a plane cutting these cylinder in the transverse direction. I would like to give different names for the surfaces in the plane (surface inside the cylinder and surface outside the cylinder)



That is my resulting code. If I change the lines in which I add coreback and cladback by removing the ,bcmod=..., I do get what I expected: cylinder walls are cyl, plane is plane and the outer boundary is outer. However, with the code as is, bcmod apparently changes all the boundaries of the volume being inserted, and not only the ones between the other volume in the bcmod call.

I would appreciate if you could point me in the direction of understanding what am I missing here.

All the best,

Francisco

Dear Mr Orlandini,

unfortunately I can not answer your questions about bcmod. However, I think that your indention can be solve by geo.AddSurface(a, b.bc("name")), where as far as I understand it the intersection of the boundaries of the solids a and b will be added as an additional surface.

---

Code:

from ngsolve import *
import netgen.gui
from netgen.csg import *

d_box = 1
r_cyl = 0.3

cube = OrthoBrick(Pnt(-d_box,-d_box,-d_box),Pnt(d_box,d_box,d_box)).bc('outer')
infcyl = Cylinder(Pnt(0,0,0), Pnt(0,0,1), r_cyl).bc('cyl')
core = ((cube*infcyl)).bc("outer").mat('core')
clad = ((cube-infcyl)).bc("outer").mat('clad')

geo = CSGeometry()
geo.Add(core)
geo.Add(clad)
for solid, bnd_name in zip((core, clad), ("core_2d", "clad_2d")):
   mid = Plane(Pnt(0,0,0), Vec(0,0,1))
   geo.AddSurface(mid, (mid*solid).bc(bnd_name))

mesh = Mesh(geo.GenerateMesh(maxh=0.1))
print(mesh.GetMaterials())
print(set(mesh.GetBoundaries()))
cf = mesh.BoundaryCF({'core_2d':1,'clad_2d':2,'cyl','outer':4},default=-1)

Draw(cf, mesh, "boundary", draw_vol = False)

---

The for loop is necessary since each solid can only hold one name for bc.

I hope that solves your problem.

Best regards,

valentin

 

Dear Mr. Valentin,

Thank you for your workaround! However, the same approach apparently does not me allow to give a name to the boundary of the 2d plane, i.e., edges on the intersection between this plane and the cube.

Would you have any guesses on how to proceed?

from ngsolve import *
import netgen.gui
from netgen.csg import *

d_box = 1
r_cyl = 0.3

cube = OrthoBrick(Pnt(-d_box,-d_box,-d_box),Pnt(d_box,d_box,d_box)).bc('outer')
infcyl = Cylinder(Pnt(0,0,0), Pnt(0,0,1), r_cyl).bc('cyl')
core = ((cube*infcyl)).bc("outer").mat('core')
clad = ((cube-infcyl)).bc("outer").mat('clad')

geo = CSGeometry()
geo.Add(core)
geo.Add(clad)
for solid, bnd_name in zip((core, clad), ("core_2d", "clad_2d")):
mid = Plane(Pnt(0,0,0), Vec(0,0,1))
geo.AddSurface(mid, (mid*solid).bc(bnd_name))
mid = Plane(Pnt(0,0,0), Vec(0,0,1))
geo.NameEdge(cube,mid,"dirichlet_2d")

mesh = Mesh(geo.GenerateMesh(maxh=0.1))
print(mesh.GetMaterials())
print(mesh.GetBoundaries())
print(mesh.GetBBoundaries())
cf = mesh.BoundaryCF({'core_2d':1,'clad_2d':2,'cyl','outer':4},default=-1)

Draw(cf, mesh, "boundary", draw_vol = False)

Dear Mr. Orlandini,

your task is indeed challenging and full of troubles. After using the "netgen-correct" order and the "netgen-correct" halfspace I came up with that solution:

Code:

---

from ngsolve import *
import netgen.gui
from netgen.csg import *


d_box = 1
r_cyl = 0.3


cube = OrthoBrick(Pnt(-d_box,-d_box,-d_box),Pnt(d_box,d_box,d_box)).bc('outer')

infcyl = Cylinder(Pnt(0,0,0), Pnt(0,0,1), r_cyl).bc('cyl')
core = ((cube*infcyl)).bc("outer").mat('core')
clad = ((cube-infcyl)).bc("outer").mat('clad')

geo = CSGeometry()

# solids
geo.Add(core)
geo.Add(clad)

# surfaces
mid1 = Plane(Pnt(0,0,0), Vec(0,0,1))
geo.AddSurface(mid1, (mid1*core).bc("core_2d"))

mid2 = Plane(Pnt(0,0,0), Vec(0,0,1))
geo.AddSurface(mid2, (mid2*clad).bc("clad_2d"))

# edges
geo.NameEdge(mid1, cube, "dirichlet_2d")

mesh = Mesh(geo.GenerateMesh(maxh=0.1))

#cheks
print(mesh.GetMaterials())
print(set(mesh.GetBoundaries()))
print(set(mesh.GetBBoundaries()))
cf = mesh.BoundaryCF({'core_2d':1,'clad_2d':2,'cyl','outer':4},default=-1)

Draw(cf, mesh, "boundary", draw_vol = False)

g = GridFunction(H1(mesh, dirichlet_bbnd="dirichlet_2d"),name='bbndry')
g.Set(CF([1 if bbnd=="dirichlet_2d" else 0 for bbnd in mesh.GetBBoundaries()]), VOL_or_BND=BBND)
Draw(g)

---

 

In case anyone is facing a similar problem, it was suggested that I tried the OCC kernel for generating the mesh, and it indeed made the code considerably simpler.

I am sharing here a code snippet in case it might be helpful to someone in the future

Thank you for the support.




from ngsolve import *
import netgen.gui
from netgen.occ import *
from numpy import random


d_box = 1
l_domain = 1
r_cyl = 0.3
el_clad = 0.2
el_core = 0.1
cube_back = Box(Pnt(-d_box,-d_box,-l_domain/2),Pnt(d_box,d_box,0)).bc('outer')
cube_front = Box(Pnt(-d_box,-d_box,0),Pnt(d_box,d_box,l_domain/2)).bc('outer')
cyl = Cylinder(Pnt(0,0,-l_domain/2),Z, r=r_cyl, h=l_domain)

clad_front = cube_front - cyl
clad_front.mat('clad').maxh = el_clad
clad_back = cube_back - cyl
clad_back.mat('clad').maxh = el_clad
core_front = cube_front * cyl
core_front.mat('core').maxh = el_core
core_back = cube_back * cyl
core_back.mat('core').maxh = el_core

clad_front.faces.Min(Z).name = 'clad_2d'
core_front.faces.Min(Z).name = 'core_2d'

domain_list = [clad_front,clad_back,core_front,core_back]
geo = OCCGeometry(domain_list)
mesh = Mesh(geo.GenerateMesh(maxh=el_core))

#checking surface domains
surflist = {"core_2d":1, "clad_2d":2, "outer"}
surf = mesh.BoundaryCF(surflist)
gu = GridFunction(H1(mesh),name='surfs')
gu.Set(surf, definedon=~mesh.Boundaries(''))
Draw(gu)