As requested I have separated this issue into it’s own post. Initially, this problem was similar to the one described in this post, but my problem has a bit more going on.
Original reply:
I am also a new user working on integrating gmsh, meshio, and fenics. I have been able to successfully run 2D models in the fenics environment using simple 2D geometry created by built-in fenics functions. However, I’m also having difficulty importing .msh files. I’ve been following a variety of tutorials, found here, here, here, and here. I started by generating a simply geometry in FreeCAD similar to the pygimli demo:
(all files I discuss here will be here)
I made sure to use the boolean fragments options in FreeCAD as described in the pygimli demo, which I believe achieves the same effect as the need for boolean fragments in .geo files as described elsewhere on this discourse site.
Next, I converted the .brep file to a .msh file as described in the pygimli demo. I added three dimension 3 physical groups, one for each object, and two dimension 2 physical groups, one for the left wall of the cube, one for the right wall of the cube. For a learning exercise, I intend to make the cylinder volume, the sphere volume, the left wall, and the right wall all different constants in a Lagrangian function space (similar to what I’ve accomplished in 2D), and leaves the bulk of the cube unrestrained. As best as I can tell, this mesh should work fine moving forward.
I then attempted to convert the .msh to the necessary .xdmf files as described in this post and various others. Because I have a 3D geometry, and I would like my FEM simulation to have both surfaces and volumes operate as sources and sinks, I think I need to export both a triangle mesh (for surfaces) and a tetrahedral mesh (for volumes).
import meshio
import numpy
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)
out_mesh = meshio.Mesh(points=mesh.points, cells={cell_type: cells}, cell_data={"name_to_read":[cell_data]})
if prune_z:
out_mesh.prune_z_0()
return out_mesh
msh = meshio.read("d4.msh")
triangleMesh = create_mesh(msh, "triangle", False)
tetraMesh = create_mesh(msh, "tetra", False)
meshio.write("d4_triag.xdmf", triangleMesh)
meshio.write("d4_tetra.xdmf", tetraMesh)
from dolfin import *
mesh = Mesh()
print(mesh.topology().dim())
mvc = MeshValueCollection("size_t", mesh, 3)
with XDMFFile("d4_tetra.xdmf") as infile:
infile.read(mesh)
infile.read(mvc, "name_to_read")
cf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
mvc = MeshValueCollection("size_t", mesh, 2)
with XDMFFile("d4_triag.xdmf") as infile:
infile.read(mvc, "name_to_read")
mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
You’ll notice in the lines mvc = MeshValueCollection("size_t", mesh, 3) and mvc = MeshValueCollection("size_t", mesh, 2) I’ve hardcoded the dimension (at least I think it’s supposed to be the dimension). This is because when I printed print(mesh.topology().dim()) I got a result of 18446744073709551615 which seems incorrect. Even if this choice was incorrect, I encountered all of the subsequent errors when using the function exactly as written above (except I took out the last line ds_custom = Measure("ds", domain=mesh, subdomain_data=mf, subdomain_id=12) since I don’t have that many subdomains, and I’m not exactly sure what it does.
My code above actually runs without errors, however, when I try to import the .xdmf files into paraview to view them, the application crashes, specifically after I click on the little eye icon next to the file. Crashes occur with Xdmf3ReaderS and Xdmf3ReaderT, while XDMF Reader exits with the following error:
ERROR: In C:\glr\builds\paraview\paraview-ci\build\superbuild\paraview\src\VTK\IO\Xdmf2\vtkXdmfReader.cxx, line 478
vtkXdmfReader (000001FD57841620): Failed to read data.
ERROR: In C:\glr\builds\paraview\paraview-ci\build\superbuild\paraview\src\VTK\Common\ExecutionModel\vtkExecutive.cxx, line 753
vtkPVCompositeDataPipeline (000001FD62CE79C0): Algorithm vtkFileSeriesReader(000001FD62DC92A0) returned failure for request: vtkInformation (000001FD09294C20)
Debug: Off
Modified Time: 382076
Reference Count: 1
Registered Events: (none)
Request: REQUEST_DATA
FORWARD_DIRECTION: 0
ALGORITHM_AFTER_FORWARD: 1
FROM_OUTPUT_PORT: 0
I’m using Windows 10 19042, Ubuntu 20 with WSL 2, FreeCAD 0.18, gmsh 4.7, Paraview 5.9.0, fenics 2019.2.0.dev0
Questions:
-
Obviously I’m interested in fixing the paraview crash and being able to see the meshes that I’ve created to see if they are correct. Does anyone have any insights?
-
I’ve read a lot on the forums about how fenics can’t or shouldn’t try to deal with mixed dimensional meshes i.e. meshes with both volumes (tetrahedrons) and surfaces (triangles). If that’s the case, why in this example are we creating separate meshes of both kinds? Shouldn’t we just be trying to create a mesh of the highest dimension of the model?
-
Much of what I’ve read about the above problem relates to the idea that a mesh output (gmsh, meshio, .step files, .brep, anything) will output both a surface mesh and a volume mesh and these mesh points don’t overlap and that this lack of overlap is what ultimately causes problems when we get to Fenics, since we will try to apply a boundary condition to a surface, but because that surface is not also considered a part of a volume, the propagation of this boundary condition will not occur in the FEM. Therefore, what we want is a mesh which is only composed of the highest dimensional construct (i.e. tetrahedrons for 3D) but which also contains enough information such that the surface can be extracted (i.e.
gmsh.model.occ.getBoundary). Is this correct? -
I’ve read much about the need for physical groups in order to correctly apply boundary conditions to a subdomain. In this example, I’ve made 3-dimension physical groups for all three volumes in my object, however, I only made two 2-dimension physical groups for the left and right walls that I wanted to apply boundary conditions to. Is this the correct procedure? If I don’t make any 2-dimension physical groups then the “triangle” parts of the
create_meshfunction fail, as expected. Do I need to add every surface to a physical group? I’ve read that gmsh will not save things that are not a part of a physical groups, however I believe all of these points should still be a part of a physical group since they are all contained within the 3 primary volumes of the project. Is it bad for points to be contained within two physical groups at once (points on the left wall of the cube are both in the left-wall-surface-group and the cube-volume-group)? Where is the most appropriate place to try to identify a subdomain? I seen many examples of users identifying subdomains at the FreeCAD or gmsh steps, and struggle to import in fenics. Is it simpler to identify subdomains in fenics, or is this impossible without also accomplishing it at previous steps?
If this doesn’t work I believe my next step will be to recreate the geometry and mesh in gmsh from the start, and not deal with the FreeCAD process, since that appears to be the workflow of so many on here. However I wonder, are there users on here who have all the bugs sorted out, and create CAD files before starting with a gmsh file, or is it usually simpler to build the geometry with pygmsh, meshio, etc.? If that’s the case, is it not simpler to build the mesh and geometry directly in fenics?
