ukaea / paramak

Create parametric 3D fusion reactor CAD models
https://paramak.readthedocs.io/en/main/
37 stars 12 forks source link

SALOME: a new open-source route for neutronics meshing ? #793

Closed RemDelaporteMathurin closed 3 years ago

RemDelaporteMathurin commented 3 years ago

Salome is an open-source meshing software. It's widely used by engineers and we've used it for our monoblock and breeding blankets meshes.

It also is scriptable. The following script produces a cylinder mesh (.med format) with marked surfaces and volume.

#!/usr/bin/env python

###
### This file is generated automatically by SALOME v9.6.0 with dump python functionality
###

import sys
import salome

salome.salome_init()
import salome_notebook
notebook = salome_notebook.NoteBook()
sys.path.insert(0, r'C:/Users/PC-REM/Desktop')

###
### GEOM component
###

import GEOM
from salome.geom import geomBuilder
import math
import SALOMEDS

geompy = geomBuilder.New()

O = geompy.MakeVertex(0, 0, 0)
OX = geompy.MakeVectorDXDYDZ(1, 0, 0)
OY = geompy.MakeVectorDXDYDZ(0, 1, 0)
OZ = geompy.MakeVectorDXDYDZ(0, 0, 1)
Cylinder_1 = geompy.MakeCylinderRH(100, 300)
inside = geompy.CreateGroup(Cylinder_1, geompy.ShapeType["SOLID"])
geompy.UnionIDs(inside, [1])
top_and_bottom = geompy.CreateGroup(Cylinder_1, geompy.ShapeType["FACE"])
geompy.UnionIDs(top_and_bottom, [10, 12])
side = geompy.CreateGroup(Cylinder_1, geompy.ShapeType["FACE"])
geompy.UnionIDs(side, [3])
[inside, top_and_bottom, side] = geompy.GetExistingSubObjects(Cylinder_1, False)
geompy.addToStudy( O, 'O' )
geompy.addToStudy( OX, 'OX' )
geompy.addToStudy( OY, 'OY' )
geompy.addToStudy( OZ, 'OZ' )
geompy.addToStudy( Cylinder_1, 'Cylinder_1' )
geompy.addToStudyInFather( Cylinder_1, inside, 'inside' )
geompy.addToStudyInFather( Cylinder_1, top_and_bottom, 'top_and_bottom' )
geompy.addToStudyInFather( Cylinder_1, side, 'side' )

###
### SMESH component
###

import  SMESH, SALOMEDS
from salome.smesh import smeshBuilder

smesh = smeshBuilder.New()
#smesh.SetEnablePublish( False ) # Set to False to avoid publish in study if not needed or in some particular situations:
                                 # multiples meshes built in parallel, complex and numerous mesh edition (performance)

Mesh_1 = smesh.Mesh(Cylinder_1)
NETGEN_1D_2D_3D = Mesh_1.Tetrahedron(algo=smeshBuilder.NETGEN_1D2D3D)
inside_1 = Mesh_1.GroupOnGeom(inside,'inside',SMESH.VOLUME)
top_and_bottom_1 = Mesh_1.GroupOnGeom(top_and_bottom,'top_and_bottom',SMESH.FACE)
side_1 = Mesh_1.GroupOnGeom(side,'side',SMESH.FACE)
isDone = Mesh_1.Compute()
[ inside_1, top_and_bottom_1, side_1 ] = Mesh_1.GetGroups()
smesh.SetName(Mesh_1, 'Mesh_1')
try:
  Mesh_1.ExportMED(r'C:/Users/PC-REM/Desktop/out.med',auto_groups=0,version=41,overwrite=1,meshPart=None,autoDimension=1)
  pass
except:
  print('ExportMED() failed. Invalid file name?')

## Set names of Mesh objects
smesh.SetName(NETGEN_1D_2D_3D.GetAlgorithm(), 'NETGEN 1D-2D-3D')
smesh.SetName(top_and_bottom_1, 'top_and_bottom')
smesh.SetName(side_1, 'side')
smesh.SetName(Mesh_1.GetMesh(), 'Mesh_1')
smesh.SetName(inside_1, 'inside')

if salome.sg.hasDesktop():
  salome.sg.updateObjBrowser()

With meshio, this MED file can be converted to h5m with the following script:

import meshio

mesh = meshio.read('mesh.med')
# mesh.points, mesh.cells, mesh.cells_dict, ...

print("This is the correspondance dict")
print(mesh.cell_tags)
# mesh.cell_tags = {-6: ['inside'], -7: ['top_and_bottom'], -8: ['side']}

# print('edges', mesh.cell_data["cell_tags"][0])
# print('tets', mesh.cell_data["cell_tags"][1])
# print('triangles', mesh.cell_data["cell_tags"][2])

# Export mesh that contains only tets
# along with tags
meshio.write_points_cells(
    "mesh_cylinder_with_tags.h5m",
    mesh.points,
    [mesh.cells[1]],
    cell_data={
        "Tet4": {
            "data": [-1*mesh.cell_data["cell_tags"][1]]
            },
        },
)

I'm afraid I haven't been able to visualise the h5m file ... @Shimwell told me this was somehow doable with moab ?

These are the produced files.

shimwell commented 3 years ago

Thanks Remi, this looks very useful and I can have a go at running the example script and visualize the mesh. Now that I see the information you have here I wonder if we should move the issue to the Dagmc or OpenMC repo where I'm sure open source generation of meshes is a topic they are keen to solve and there are more expertise there.

RemDelaporteMathurin commented 3 years ago

Shall I open an issue there then ?

RemDelaporteMathurin commented 3 years ago

Opened in https://github.com/openmc-dev/openmc/issues/1810 closing this.