Hi, folks!
I’ve been struggling with incorporating gmsh into my fenics workflow. I encounter multiple errors while reading a 3D mesh with physical groups I created using gmsh with python. I’ve been using the thread from Transitioning from mesh.xml to mesh.xdmf, from dolfin-convert to meshio - #166 by dokken to do my code and it doesn’t work.
If I use the lines cell_markers = … and facet_markers = …, it gives an error whose reason is “Illegal topological dimension (3)”; whereas if I comment these lines and use a straightforward CompiledSubDomain command, it says that the facets could not be found to apply boundary conditions. When I run the code using debugger, I get the following error in the reading command: ValueError: cannot create std::vector larger than max_size()
I don’t have the slightliest clue on how to proceed. There follows a simple toy example which presents exactly the same issues I’m reporting:
import gmsh
import meshio
import numpy as np
from dolfin import *
# Initialize GMSH
gmsh.initialize()
gmsh.model.add("cube_with_physical_groups")
# Define the cube geometry
lc = 1.0 # Mesh element size
gmsh.model.occ.addBox(0, 0, 0, 15, 15, 25) # Cube from (0,0,0) to (1,1,1)
# Synchronize the CAD representation
gmsh.model.occ.synchronize()
# Define physical groups
volume = gmsh.model.getEntities(dim=3) # Get the volume entity (dim=3)
surfaces = gmsh.model.getEntities(dim=2) # Get the surface entities (dim=2)
# Add a physical group for the volume
gmsh.model.addPhysicalGroup(3, [volume[0][1]], tag=1) # Tag = 1 for volume
gmsh.model.setPhysicalName(3, 1, "Volume")
# Add physical groups for the surfaces
for i, surface in enumerate(surfaces):
gmsh.model.addPhysicalGroup(2, [surface[1]], tag=101 + i) # Tags: 101, 102, ...
gmsh.model.setPhysicalName(2, 101 + i, f"Surface_{i+1}")
# Generate the mesh
gmsh.model.mesh.generate(3)
# Save the mesh to an MSH file
gmsh.write("cube_with_physical_groups.msh")
gmsh.fltk.run()
# Finalize GMSH API
gmsh.finalize()
# Convert the GMSH mesh to a meshio.Mesh with domains and boundaries
mesh = meshio.read("cube_with_physical_groups.msh")
# Extract tetrahedral (volume) and triangular (surface) cells
tetra_cells = mesh.get_cells_type("tetra")
triangle_cells = mesh.get_cells_type("triangle")
# Extract physical group data
physical_volumes = mesh.get_cell_data("gmsh:physical", "tetra")
physical_surfaces = mesh.get_cell_data("gmsh:physical", "triangle")
# Create a mesh of tetrahedrons
tetra_mesh = meshio.Mesh(points=mesh.points, cells={"tetra": tetra_cells},
cell_data={"domain": [physical_volumes]})
# Create a mesh of triangles
triangle_mesh = meshio.Mesh(points=mesh.points, cells={"triangle": triangle_cells},
cell_data={"boundary": [physical_surfaces]})
# Save the new mesh to a file
meshio.write("cube_with_labels.xdmf", tetra_mesh)
meshio.write("boundary_mesh.xdmf", triangle_mesh)
mesh = Mesh()
domain_data = MeshValueCollection("size_t", mesh, 3)
boundary_data = MeshValueCollection("size_t", mesh, 2)
with XDMFFile("cube_with_labels.xdmf") as infile:
infile.read(mesh)
infile.read(domain_data, "domain")
with XDMFFile("boundary_mesh.xdmf") as infile:
infile.read(mesh)
infile.read(boundary_data, "boundary")
# Sets the mesh functions
cell_markers = cpp.mesh.MeshFunctionSizet(mesh, domain_data)
facet_markers = cpp.mesh.MeshFunctionSizet(mesh, boundary_data)
#dx = Measure("dx", domain=mesh, subdomain_data=cell_markers)
#ds = Measure("ds", domain=mesh, subdomain_data=facet_markers)
# Defines a vectorial function space
V = VectorFunctionSpace(mesh, "Lagrange", degree=2)
# Defines a trial function
u = TrialFunction(V)
# Defines a test function
du = TestFunction(V)
# Defines Dirichlet boundary conditions
#bc_dirichlet = [DirichletBC(V.sub(2), expression_below, facet_markers, 5),
#DirichletBC(V.sub(0), expression_side1, facet_markers, 8),
#DirichletBC(V.sub(1), expression_side2, facet_markers, 9)]
#bc_dirichlet = [DirichletBC(V, Constant("0.0", "0.0", "0.0"),
#facet_markers, 8)]
below = CompiledSubDomain('near(x[2],z_boundary) && on_boundary', z_boundary=0.0)
side1 = CompiledSubDomain('near(x[0],x_boundary) && on_boundary', x_boundary=0.0)
side2 = CompiledSubDomain('near(x[1],y_boundary) && on_boundary', y_boundary=0.0)
expression_below = Constant("0.0")#Expression(("0.0", "0.0", "0.0"), degree=1)
expression_side1 = Constant("0.0")
expression_side2 = Constant("0.0")
bc_dirichlet = [DirichletBC(V.sub(2), expression_below, below),
DirichletBC(V.sub(0), expression_side1, side1),
DirichletBC(V.sub(1), expression_side2, side2)]
# Defines Neumann boundary conditions
above = CompiledSubDomain('near(x[2],z_boundary) && on_boundary', z_boundary=25.0)
boundaries = MeshFunction("size_t", mesh, mesh.topology().dim()-1,0)
above.mark(boundaries,1)
ds = Measure('ds', domain=mesh, subdomain_data=boundaries)
# Defines Neumann boundary conditions
# Defines the traction vector
T = Constant(("0.0", "0.0", "100000.0"))
# Defines the body forces vector
B = Constant(("0.0","0.0","0.0"))
###################################################
# Tensorial #
###################################################
epsilon = sym(grad(u))
# Defines materials properties
E = 200E6
nu = 0.3
mu = E/(2*(1+nu))
lmbda = (E*nu)/((1+nu)*(1-2*nu))
# Defines the stress tensor
I = Identity(3)
cauchy_stress = (lmbda*tr(epsilon)*I)+(2*mu*epsilon)
# Defines the variation of the infinitesimal strain
# tensor
d_epsilon = sym(grad(du))
###################################################
# Variational #
# #################################################
# Defines the bilinear form
a = inner(cauchy_stress,d_epsilon)*dx
# Defines the linear form
L = (dot(B,du)*dx)+(dot(T,du)*ds(1))
# Solves the system
u = Function(V)
solve(a==L, u, bc_dirichlet)
# Creates the displacement file
file = File("Results//displacement.pvd")
file << u
# Evaluates stress and all its components
W = TensorFunctionSpace(mesh, "Lagrange", 1)
cauchy_function = 2*mu*sym(grad(u)) + lmbda*tr(grad(u))*Identity(3)
sigma_w = project(cauchy_function, W)
file_sigma = File("Results//stress.pvd")
file_sigma << sigma_w