Marking materials in 3D CSGeometry

Dear NGSolve community,

As a new NGSolve user, I have recently been looking into the possibility of marking regions of a domain to indicate different materials. I have successfully been able to mark regions of a 2D domain, but I am experiencing some (to me) surprising behavior in 3D. As a MWE, I consider two stacked cubes, where I want to mark the bottom cube as 'steel' and the top cube as 'foam'. For this I use the following code:

[code]
import ngsolve as ng
import netgen.csg as csg
from ngsolve.webgui import Draw

brick1 = csg.OrthoBrick(csg.Pnt(0,0,0), csg.Pnt(1,1,1)).mat("steel")
brick2 = csg.OrthoBrick(csg.Pnt(0,0,1), csg.Pnt(1,1,2)).mat("foam")

geo = csg.CSGeometry()
geo.Add(brick1)
geo.Add(brick2)
mesh = ng.Mesh(geo.GenerateMesh(maxh=0.2))

cf = mesh.RegionCF(ng.VOL, {"steel":1, "foam":2})
Draw(cf, mesh)
[\code]

Given the definition of my coefficient function, I expected the output visualization to show that one cube with value 1 and one cube with value 2. However, in the actual output (see attachment or link below), one cube face has value 1, another has value 2, and the rest of the domain has value 0.

I do not understand why my coefficient function behaves this way. Is it because of a bug, or have a just set up things incorrectly? If the latter, any tips that could point me in the right direction for fixing my code would be greatly appreciated.

Link to output (since I am not sure my attachment uploaded correctly): [url] github.com/blakseth/images/blob/main/output.PNG [/url]

Best regards,
Sindre

Hi Sindre,

to get what you want is the following:

from ngsolve.comp import BoundaryFromVolumeCF
Draw (BoundaryFromVolumeCF(cf), mesh)


We have material indices, and boundary condition indices.
Before the cf was indexed by the bc-index. Now the evaluation at the boundary takes the values from the attached volume, and so it takes the (expected) material index.

best, Joachim

Excellent, thank you Joachim.

I have now moved on to a slightly more complex geometry with a steel cylinder passing through a foam box. Your solution works well on the outside of the box and on the ends of the cylinder, but there is something not quite right on the curved interface between the cylinder and the foam. As you can see from the output (link below), most of the curved surface is shown as steel, but some elements are shown as foam. Ideally, this should not happen, since the surface is part of the cylinder marked as steel, and the foam is defined as being everything inside the box that is not part of the cylinder. If this is just a visualization error, that is okay, but if some vertices on the cylinder are actually marked "foam", that is something I would want to fix. Is there any way I can confirm that all the vertices on the cylinder's curved surface are actually marked as steel, e.g. printing the material of all vertices on the cylinder surface?

Best regards,
Sindre

Link to output visualization: [url] github.com/blakseth/images/blob/main/cyl...r_with_artifacts.PNG [/url]

My code:
import ngsolve as ng
import netgen.csg as csg
from ngsolve.webgui import Draw
from ngsolve.comp import BoundaryFromVolumeCF

w,d,h = 2,2,1
r = 0.5
assert r < 0.5*w and r < 0.5*d

supp_maxh = 0.05
foam_maxh = 0.2

left  = csg.Plane(csg.Pnt(0,0,0), csg.Vec(-1,0,0)).bc("sides")
right = csg.Plane(csg.Pnt(w,d,h), csg.Vec( 1,0,0)).bc("sides")
front = csg.Plane(csg.Pnt(0,0,0), csg.Vec(0,-1,0)).bc("sides")
back  = csg.Plane(csg.Pnt(w,d,h), csg.Vec(0, 1,0)).bc("sides")
bot   = csg.Plane(csg.Pnt(0,0,0), csg.Vec(0,0,-1)).bc("bot")
top   = csg.Plane(csg.Pnt(w,d,h), csg.Vec(0,0, 1)).bc("top")

box = left * right * front * back * bot * top

inf_cyl = csg.Cylinder(csg.Pnt(w/2,d/2,0), csg.Pnt(w/2,d/2,1), r)  # Defines a cylinder of infinite height with radius r and axis intersecting the given points.

supp = box*inf_cyl # The cylindrical support is the intersection of the box and the infinite cylinder.
foam = box - supp  # The foam fills the part of the prism not occupied by the support.

supp.mat("steel")
foam.mat("foam")

geo = csg.CSGeometry()
geo.Add(supp)
geo.Add(foam)

mesh = ng.Mesh(geo.GenerateMesh(maxh=0.2))
mesh.Curve(3)

vol_cf = mesh.RegionCF(ng.VOL, {"steel":1, "foam":2})
bnd_cf = mesh.RegionCF(ng.BND, {"top":1, "bot":2, "sides"})

Draw(bnd_cf, mesh)
Draw(BoundaryFromVolumeCF(vol_cf), mesh)
 

Hi Sindre,

we can certainly fix it with the original Netgen CSG modeler, but I recommend to switch to occ-geometry, which is much mightier 

Joachim
 

Thanks Joachim, then I will look into the occ-geometry. Thanks for your guidance!
Best,
Sindre