Recently, an issue occurred while loading an STL file found while searching on Google into dolfinx. The reason I posted to ask
“if there is a way to construct a mesh by loading an STL file in the dolfinx environment.”
After reading other postings, I saw that it was possible to load xdmf or msh format files in dolfinx. (I’m not sure… )
If STL format files cannot be loaded directly in the dolfinx environment, is there a way to change STL format files to xdmf or msh format within the code? If there is no method, should I change it to an external program and load it in the dolfinx environment only?
If I can load a stl file in the dolfinx environment, can I also load a shape made up of multiple parts (a product with multiple physical properties combined)?
Please note that an stl mesh is a surface mesh only, and does not mesh the interior of a 3D domain, thus to solve something on a 3D shape, you have to use a mesh generator to mesh the interior of the mesh.
And I tried to extend the FEniCSx tutorial code from 2D specified on (Homepage) by dokken to 3D as shown in the code below.
(Because the stl files I want to use are 3D shapes.)
import meshio
import gmsh
import numpy as np
from mpi4py import MPI
from pathlib import Path
# make directory folder
Path("foo").mkdir(parents=True, exist_ok=True)
# create mesh function (consider "gdim")
def create_mesh(mesh, cell_type, gdim, prune_z=False):
cells = mesh.get_cells_type(cell_type)
cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
points = mesh.points[:,:gdim] if prune_z else mesh.points
out_mesh = meshio.Mesh(points=points, cells={cell_type: cells}, cell_data={"name_to_read":[cell_data]})
return out_mesh
gmsh.initialize()
proc = MPI.COMM_WORLD.rank
top_marker = 2
bottom_marker = 1
if proc == 0:
# We create one rectangle for each subdomain
# I think this line should load the stl file. (Like "gmsh.model.add"...?)
gmsh.model.occ.addBox(0, 0, 0, 1, 0.5, 1, tag=1)
gmsh.model.occ.addBox(0, 0.5, 0, 1, 0.5, 1, tag=2)
# We fuse the two rectangles and keep the interface between them
# How can I achieve sync when loading an stl file?
gmsh.model.occ.fragment([(2,1)],[(2,2)])
gmsh.model.occ.synchronize()
# Mark the top (2) and bottom (1) box
# If we consider it to be a 3D object, how can we set the interior mesh inside it?
top, bottom = None, None
for bx in gmsh.model.getEntities(dim=3):
com = gmsh.model.occ.getCenterOfMass(bx[0], bx[1])
if np.allclose(com, [0.5, 0.25, 0.5]):
bottom = bx[1]
else:
top = bx[1]
gmsh.model.addPhysicalGroup(2, [bottom], bottom_marker)
gmsh.model.addPhysicalGroup(2, [top], top_marker)
gmsh.write("foo/mesh.msh")
# Read in mesh
msh = meshio.read("foo/mesh.msh")
# Create and save one file for the mesh, and one file for the facets
triangle_mesh = create_mesh(msh, "triangle", 3, prune_z=True)
meshio.write("foo/tri.xdmf", triangle_mesh)
gmsh.finalize()
When I ran the code and checked the results, I could see it as shown in the picture below.
I have additional questions (I’ve commented them in the code).
a. In the code above, can’t a mesh be assigned to the interior mesh of 3D object? Is there another way?
(My idea is that if it is not possible to assign a mesh to the inside of the 3D domain loaded as an stl file, we should declare a structural grid larger than the stl object and then mapping the domain loaded as stl file inside this grid. Just an idea…)
b. And when multiple subdomains exist in one object, can each subdomain be recognized differently?
(For example, I wanted to ask if the entire domain can be set as an active domain, but some parts of the entire domain can be set as a solid domain or rigid domain.)
I thought it would be easy to build the read stl file code, but it’s difficult… Thank you again for the answer!
How do you want to embed this inside a 3D object? Should the 3D mesh align with the 2D surface? If not, how do you want to couple the two meshes? You talk about a mapping, but it is not very concrete.
Depending on your simulation, you can pin degrees of freedom by applying Dirichlet conditions on them.
It was overlooked that was difficult for the geometry points in the 3D stl file to be accurately aligned with the 3D structural grid. It’s an idea that wasn’t thought through in depth. I’m sorry if I offended you.
To ask a more specific question, let’s say I’m trying to import the jet engine bracket example-(link)stl file into dolfinx. (If there is a problem with sharing the link, I will delete the post.)
The entire design domain(jetEngineDesignDomainFine.stl) is as shown below, and if two subdomains exist as shown below,
In the bracket shape, there is a solid domain(jetEngineSolidDomainFine.stl) that will configure the Dirichlet condition and a rigid domain(cylinder in middle of figure, jetEngineRigidDomainFine.stl) that will be recognized as a different material.
Then,
Can the stl file that makes up a solid domain be recognized as a Dirichlet condition for the entire domain?
As in the example written in the comment above, I think it is necessary to know the coordinates of all domains in contact to sync between rigid domain and entire domain. Is there a way to automatically sync?
(Also, This is a minor question, but what command should I use to add a model to an stl file using gmsh to make mesh file?)
Is it possible to extract the three files above into one output file?
Thank you for your help whenever I run into algorithmic difficulties.
STL is still just a surface mesh. If you mesh your domain, you should mark the appropriate boundaries of your volumetric mesh (in for instance GMSH), so that They can be transferred into DOLFINx
Thank you for reply, @dokken!
I created a t13.msh file using the code provided in the under link.
Then, I tried to load the file into dolfinx using the code provided in this tutorial link,
from dolfinx import io
from mpi4py import MPI
proc = MPI.COMM_WORLD.rank
mesh, cell_markers, facet_markers = io.gmshio.read_from_msh("t13.msh", MPI.COMM_WORLD, gdim=3)
but this IndexError occurred.
Info : Reading 't13.msh'...
Info : 36 entities
Info : 4299 nodes
Info : 22196 elements
Info : 11 parametrizations
Info : [ 10%] Discrete surface 4 is planar, simplifying parametrization
Info : [ 60%] Discrete surface 9 is planar, simplifying parametrization
Info : [ 70%] Discrete surface 10 is planar, simplifying parametrization
Info : [ 80%] Discrete surface 11 is planar, simplifying parametrization
Info : [ 90%] Discrete surface 12 is planar, simplifying parametrization
Info : Done reading 't13.msh'
IndexError: index -1 is out of bounds for axis 0 with size 0
I wondered if there was a problem with the t13.msh file, so I converted the msh file to an xdmf file using the corresponding command line,
meshio convert t13.msh t13.xdmf
and when I checked the file through ParaView, I was able to load it as shown below.
Can I ask how to resolve this error? Is there any solution so that loading the converted msh file into dolfinx and executing commands such as FunctionSpace?
Thank you for answer. @dokken
I’m late in replying because I was trying to find out related knowledge.
I was able to assign physical entities to the msh file through
the link you provided.
I have additional question.
a) According to the msh file and the stated in this link “Convert msh-files to XDMF using meshio”
import meshio
# create mesh function
def create_mesh(mesh, cell_type, prune_z=False):
cells = mesh.get_cells_type(cell_type)
cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
points = mesh.points[:,:3] if prune_z else mesh.points
out_mesh = meshio.Mesh(points=points, cells={cell_type: cells}, cell_data={"name_to_read":[cell_data]})
return out_mesh
# read in mesh
msh = meshio.read("t13.msh")
# Create and save one file for the mesh, and one file for the facets
triangle_mesh = create_mesh(msh, "triangle", prune_z=True)
# Is it possible...?
# quad_mesh = create_mesh(msh, "quadrangle", prune_z=True)
# hex_mesh = create_mesh(msh, "hexahedral", prune_z=True)
meshio.write("t13.xdmf", triangle_mesh)
As shown in the code above, the msh file was converted to an xdmf file. Then, can’t I declare a hexahedral mesh(hex_mesh) or quad mesh(quad_mesh) other than a triangular mesh?
)
)
