Hi all,
I am trying to mark the boundaries of my 2D mesh which I create with a gmsh script. I combined a couple of solutions from various topics within the community. Unfortunately, I end up with a strange result. My 2D mesh somehow turn in to a 3D mesh.
Step 1) I create “mesh.msh”, “mesh.xml” using gmsh.
Step 2) Read the “mesh.msh” file using meshio.
Step 3) Create “mesh.xdmf” and “mf.xdmf” files to obtain the marked boundaries that I already defined within the gmsh script.
Step 4) Read the “mesh.xdmf” and “mf.xdmf” files and write them to a .pvd file so that I can check if my boundaries are marked correctly.
Step 5) Create my mixed function space and check if the dimension of the mesh is alright.
Here is a simplified working code that shows my problem:
(you will get the 3d result)
from fenics import *
from mshr import *
import numpy as np
import matplotlib.pyplot as plt
from subprocess import call
%matplotlib inline
# vertices of the domain (in anticlockwise order)
V=[[],[],[],[]]
V[0]=[-1,-1]
V[1]=[ 1,-1]
V[2]=[ 1, 1]
V[3]=[-1, 1]
# Corners of the unit cell
point_script="\n// Mesh definition\n// Outer boundary\n"
for i in np.arange(4):
point_script+="Point({0})={{{1},{2},0,gridsize}};\n".format(i+1,V[i][0],V[i][1])
# Convert to a gmsh script, run gmsh and then use dolfin-convert to get the mesh into xml format
param_script="// Inputs\n\ngridsize=1/{0};\n".format(10)
line_script="Line(1) = {1, 2};\nLine(2) = {2, 3};\nLine(3) = {3, 4};\nLine(4) = {4, 1};\n"
line_loop_script="Line Loop(5) = {4, 1, 2, 3};\n"
periodic_bc_script="Periodic Line {1}={3};\n"
surface_script="Plane Surface(6) = {5};\n"
physical_script="Physical Surface(7) = {6};\nPhysical Line(8) = {1,2};\nPhysical Line(9) = {3,4};"
# Total script
mesh_script=param_script+point_script+line_script+line_loop_script+periodic_bc_script+surface_script+physical_script
# Run gmsh and dolfin-convert to get the xml mesh file
print("Running Gmsh and dolfin-convert...")
geo_file=open("mesh.geo",'w')
geo_file.write(mesh_script)
geo_file.close()
gmsh_err=call(["gmsh", "-2", "mesh.geo", "-o", "mesh.msh"])
if gmsh_err:
print("...something bad happened when the mesh was being generated...")
else:
dolfin_err=call(["dolfin-convert", "mesh.msh", "mesh.xml"])
if dolfin_err:
print("...something bad happened when the mesh was being generated...")
else:
print("...mesh generated successfully!")
import meshio
import numpy as np
msh = meshio.read("mesh.msh")
line_cells = []
for cell in msh.cells:
if cell.type == "triangle":
triangle_cells = cell.data
elif cell.type == "line":
if len(line_cells) == 0:
line_cells = cell.data
else:
line_cells = np.vstack([line_cells, cell.data])
line_data = []
for key in msh.cell_data_dict["gmsh:physical"].keys():
if key == "line":
if len(line_data) == 0:
line_data = msh.cell_data_dict["gmsh:physical"][key]
else:
line_data = np.vstack([line_data, msh.cell_data_dict["gmsh:physical"][key]])
elif key == "triangle":
triangle_data = msh.cell_data_dict["gmsh:physical"][key]
triangle_mesh = meshio.Mesh(points=msh.points, cells={"triangle": triangle_cells})
line_mesh =meshio.Mesh(points=msh.points,
cells=[("line", line_cells)],
cell_data={"name_to_read":[line_data]})
meshio.write("mesh.xdmf", triangle_mesh)
meshio.xdmf.write("mf.xdmf", line_mesh)
from dolfin import *
mesh = Mesh()
with XDMFFile("mesh.xdmf") as infile:
infile.read(mesh)
mvc = MeshValueCollection("size_t", mesh, 2)
with XDMFFile("mf.xdmf") as infile:
infile.read(mvc, "name_to_read")
mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
File("dolfinmesh.pvd").write(mesh)
File("dolfincellfunc.pvd").write(mf)
# Check if the dimensions of the marked boundaries are OK
ds_custom = Measure("ds", domain=mesh, subdomain_data=mf, subdomain_id=8)
print(assemble(1*ds_custom))
ds_custom = Measure("ds", domain=mesh, subdomain_data=mf, subdomain_id=9)
print(assemble(1*ds_custom))
# Create the mixed function-space
V_1 = FiniteElement("Lagrange", mesh.ufl_cell(), 2, 2)
V_2 = FiniteElement("Lagrange", mesh.ufl_cell(), 2, 2)
V_Mixed = FunctionSpace(mesh, V_1*V_2)
# HERE IS THE STRANGE RESULT: Now it turned into a 3D mesh?
print("dim =", V_Mixed.ufl_cell().geometric_dimension())
I am definitely making a major mistake but I am a little bit lost with the codes of meshio right now. I just took whatever made sense to me for my purpose from the discussion within the community. Also, I would really appreciate if there is a much simpler way to do this because I end up with 13 files in my folder just to do this operation. I just want to mark my boundaries in the end.
Thank you!