Possible to set boundary conditions on a geometry from an STL source?

Main Forum Old-Forum - do not post here
1

Hi,
I have a CSG geometry defined in an STL file which includes a small detailed model of a sailplane inside a large bounding box. I’m able to load this geometry in Netgen.

I have been looking at I-Tutorial 4.2 on CSG geometries in 3D. In that tutorial, I can see how to set boundary and material labels for a CSG object such as a sphere, then add that object to the geometry before meshing. I would like to be able to do something similar with my geometry from the STL file, for example, to set one boundary label on the surface of the airplane and a different label on the surface of the bounding box.

When I look at the contents of the STL file, I see only a list of facet definitions. Is there some way to assign boundary labels to specific facets either in Netgen/Python or directly in the STL file? I may not be thinking about this problem the right way – any suggestions would be most welcome! I’m attaching the STL file in case that is helpful…

Thanks!

Attachment: Ak8.stl

2

Hi,
It seems like your model is placed on the surface of your bounding box and not inside, there may be a mistake in the stl file. But anyhow I think it’s easier to mesh the stl model and the bounding CSG seperately and merge them afterwards. There is a (not fully up to date) documentation on how to do that here:
https://ngsolve.org/docu/latest/netgen_tutorials/working_with_meshes.html
When adding the facedescriptors you can give each of the surfaces a different number and a different name, so that they represent different boundary conditions.
Some new stuff that is not in that docu example:
You can mesh only surfaces (this saves quite some time) by:

import netgen.meshing as meshing
stlgeo.GenerateMesh(maxh=maxh,perfstepsend= meshing.MeshingStep.MESHSURFACE)

Then copy all the surface elements together into one mesh, set the maximum meshsize and name for each domain and generate the volume mesh:

# material numbers are 1 based
mesh.SetMaterial(1,"mat1")
mesh.SetMaterial(2,"mat2")
meshsize = [maxh[dom] for dom in mesh.GetMaterials()]
mesh.SetMaxHDomain(meshsize)
mesh.CalcLocalH(0.5)
mesh.GenerateVolumeMesh()

The 0.5 in CalcLocalH is a grading factor.
You can set names for the bc numbers with

# boundary condition numbers are 0 based! - yes this is a bit annoying... ;)
mesh.SetBCName(number,name)

I’ve attached a file that can be used to merge the surface meshes by calling

add_surfaces(mesh1,mesh2,outerdomain)

If you have any questions or if the code isn’t working (I use it in a bigger module, there may be something missing…) feel free to ask :wink:

PS:
To check if everything is set correctly I use a short snippet:

    from ngsolve import *
    from ngsolve.internal import *
    boundaries = ["bc1",...]
    materials = ["mat1",...]
    fes = H1(mesh)
    u = GridFunction(fes,"u")
    viewoptions.clipping.enable = 1
    Draw(u)
    visoptions.clipsolution = "scal"
    for bc in boundaries:
        u.Set(CoefficientFunction(1),definedon=mesh.Boundaries(bc))
        Redraw()
        input(bc)
    for mat in materials:
        u.Set(CoefficientFunction(1),definedon=mesh.Materials(mat))
        Redraw()
        input(mat)

Best
Christopher

Attachment: merge_meshes.py

3

Hi Christopher,

In my version 6.2.1806 is missing:

from ngsolve_mri.meshtools import max_surfnr, nr_bcs, nr_materials

Where can I find this missing part?

4

Sorry that file was missing, I’ve attached it. Just put it in the same folder and change

from ngsolve_mri.meshtools import max_surfnr, nr_bcs, nr_materials

to

from .meshtools import max_surfnr, nr_bcs, nr_materials

Attachment: meshtools.py

5

Dear Christopher,

I encountered a segmentation fault error during the conversion from Netgen Mesh to NGsolve Mesh. However, saving the mesh and then importing it resolved the issue. Any thoughts on why this occurred?

Best regards,
Vien Che

6

phu, that is quite old code and likely that something changed in structure. I’d guess on face descriptor boundary names or so… can you get a traceback from the debugger?

7

Dear Christopher,

This is what I observed:

Conversion to NGSolve mesh…
set global mesh
Update mesh topology
Update edges
Update faces
double free or corruption (out)
Aborted (core dumped)

The code to reproduce the issue is given below:

from netgen.meshing import MeshingStep

from merge_meshes import add_surfaces

import netgen.occ as occ

from ngsolve import ngsglobals, Mesh

import time

ts = time.time()

ngsglobals.msg_level = 10


box = occ.OCCGeometry(occ.Box(occ.Pnt(0,0,0), occ.Pnt(1.0,1.0,1.0)))
mesh1 = box.GenerateMesh(maxh=0.5, perfstepsend= MeshingStep.MESHSURFACE)
mesh1.Refine()
sphere = occ.OCCGeometry(occ.Sphere(occ.Pnt(0.5,0.5,0.5), 0.2))
mesh2 = sphere.GenerateMesh(maxh=0.5, perfstepsend= MeshingStep.MESHSURFACE)
mesh2.Refine()
add_surfaces(mesh1, mesh2)
print("Meshing volume...")
mesh1.GenerateVolumeMesh()
print("Saving volume mesh...")
mesh1.Save("test_simple_mesh.vol.gz")
print("Time meshing: ", time.time()-ts)
print("Conversion to NGSolve mesh...")
mesh = Mesh(mesh1) #segmenation fault? Why?
mesh = Mesh("test_simple_mesh.vol.gz")

Best regards,
Vien Che

8

Dear Christopher,

I can convert Netgen Mesh to NGSolve mesh if I create an empty mesh and then use

add_surfaces(ngmesh, mesh1)
add_surfaces(ngmesh, mesh2)

The boundaries are correct. However, the materials are not. I got only one material instead of two.

from netgen.meshing import MeshingStep
from merge_meshes import add_surfaces
import netgen.occ as occ
from ngsolve import ngsglobals, Mesh
from netgen.meshing import Mesh as ngMesh
import time
from collections import Counter

ts = time.time()
ngsglobals.msg_level = 2
box = occ.Box(occ.Pnt(0,0,0), occ.Pnt(1.0,1.0,1.0))
box.bc("boxbc")
box.mat("boxmat")
box = occ.OCCGeometry(box)
mesh1 = box.GenerateMesh(maxh=0.5, perfstepsend= MeshingStep.MESHSURFACE)
mesh1.Refine()
sphere = occ.Sphere(occ.Pnt(0.5,0.5,0.5), 0.2)
sphere.bc("spherebc")
sphere.mat("spheremat")
sphere = occ.OCCGeometry(sphere)
mesh2 = sphere.GenerateMesh(maxh=0.5, perfstepsend= MeshingStep.MESHSURFACE)
mesh2.Refine()
ngmesh = ngMesh() # create an empty mesh
add_surfaces(ngmesh, mesh1)
add_surfaces(ngmesh, mesh2)
print("Meshing volume...")
ngmesh.SetMaterial(1,"boxmat")
ngmesh.SetMaterial(2,"spheremat")
ngmesh.GenerateVolumeMesh()
print("Conversion to NGSolve mesh...")
mesh = Mesh(ngmesh)

bc_list = mesh.GetBoundaries()
mat_list = mesh.GetMaterials()
print(Counter(mat_list).keys())
print(Counter(mat_list).values())