Hello,
I would like to generate a 2d mesh by using gmsh.model.occ but I am having some issues. I generate the mesh with
import meshio
import argparse
import gmsh
import warnings
parser = argparse.ArgumentParser()
parser.add_argument( "resolution" )
args = parser.parse_args()
warnings.filterwarnings( "ignore" )
gmsh.initialize()
gmsh.model.add( "my model" )
c_r = [0, 0, 0]
c_R = [0, 0, 0]
r = 1
R = 2
resolution = (float)( args.resolution )
print( f"Mesh resolution = {resolution}" )
disk_r = gmsh.model.occ.addDisk( c_r[0], c_r[1], c_r[2], r, r )
disk_R = gmsh.model.occ.addDisk( c_R[0], c_R[1], c_R[2], R, R )
ring = gmsh.model.occ.cut( [(2, disk_R)], [(2, disk_r)] )
p_1 = gmsh.model.occ.addPoint( (r + R) / 2.0, 0, 0 )
p_2 = gmsh.model.occ.addPoint( (r + R) / 2.0, r / 10.0, 0 )
line_p_1_p_2 = gmsh.model.occ.addLine( p_1, p_2 )
gmsh.model.occ.synchronize()
# add 2-dimensional objects
surfaces = gmsh.model.occ.getEntities( dim=2 )
assert surfaces == ring[0]
disk_subdomain_id = 6
gmsh.model.addPhysicalGroup( surfaces[0][0], [surfaces[0][1]], disk_subdomain_id )
gmsh.model.setPhysicalName( surfaces[0][0], disk_subdomain_id, "disk" )
# add 1-dimensional objects
lines = gmsh.model.occ.getEntities( dim=1 )
circle_r_subdomain_id = 1
circle_R_subdomain_id = 2
line_p_1_p_2_subdomain_id = 3
gmsh.model.addPhysicalGroup( lines[0][0], [lines[0][1]], circle_r_subdomain_id )
gmsh.model.setPhysicalName( lines[0][0], circle_r_subdomain_id, "circle_r" )
gmsh.model.addPhysicalGroup( lines[1][0], [lines[1][1]], circle_R_subdomain_id )
gmsh.model.setPhysicalName( lines[1][0], circle_R_subdomain_id, "circle_R" )
gmsh.model.addPhysicalGroup( lines[2][0], [lines[2][1]], line_p_1_p_2_subdomain_id )
gmsh.model.setPhysicalName( lines[2][0], line_p_1_p_2_subdomain_id, "line_p_1_p_2" )
# add 0-dimensional objects
vertices = gmsh.model.occ.getEntities( dim=0 )
p_1_subdomain_id = 4
p_2_subdomain_id = 5
gmsh.model.addPhysicalGroup( vertices[0][0], [vertices[0][1]], p_1_subdomain_id )
gmsh.model.setPhysicalName( vertices[0][0], p_1_subdomain_id, "p_1" )
gmsh.model.addPhysicalGroup( vertices[1][0], [vertices[1][1]], p_2_subdomain_id )
gmsh.model.setPhysicalName( vertices[1][0], p_2_subdomain_id, "p_2" )
gmsh.model.occ.synchronize()
gmsh.model.mesh.generate( 2 )
gmsh.write( "mesh.msh" )
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]} )
return out_mesh
mesh_from_file = meshio.read( "mesh.msh" )
# create a triangle mesh in which the surfaces will be stored
triangle_mesh = create_mesh( mesh_from_file, "triangle", prune_z=True )
meshio.write( "triangle_mesh.xdmf", triangle_mesh )
# create a line mesh
line_mesh = create_mesh( mesh_from_file, "line", True )
meshio.write( "line_mesh.xdmf", line_mesh )
# create a vertex mesh
vertex_mesh = create_mesh( mesh_from_file, "vertex", True )
meshio.write( "vertex_mesh.xdmf", vertex_mesh )
which creates the mesh and stores its triangles, lines and vertices in the respective xdmf files:
$ python3 generate_2dmesh_ring.py 0.1
Mesh resolution = 0.1
Info : Meshing 1D...
Info : [ 0%] Meshing curve 1 (Ellipse)
Info : [ 40%] Meshing curve 2 (Ellipse)
Info : [ 70%] Meshing curve 3 (Line)
Info : Done meshing 1D (Wall 0.00192213s, CPU 0.001924s)
Info : Meshing 2D...
Info : Meshing surface 2 (Plane, Frontal-Delaunay)
Info : Done meshing 2D (Wall 0.0109537s, CPU 0.010954s)
Info : 61 nodes 123 elements
Info : Writing 'mesh.msh'...
Info : Done writing 'mesh.msh'
However, when I try to read it back and tag dx and the dss with the following code
from fenics import *
from mshr import *
from dolfin import *
from dolfin import *
import numpy as np
#CHANGE PARAMETERS HERE
c_r = [0, 0, 0]
c_R = [0, 0, 0]
r = 1
R = 2
#CHANGE PARAMETERS HERE
# read the mesh
mesh = Mesh()
xdmf = XDMFFile( mesh.mpi_comm(), "triangle_mesh.xdmf" )
xdmf.read( mesh )
print(f"Mesh dimension = {mesh.topology().dim()}")
# read the triangles
mvc = MeshValueCollection( "size_t", mesh, mesh.topology().dim() )
with XDMFFile( "triangle_mesh.xdmf" ) as infile:
infile.read( mvc, "name_to_read" )
vf = cpp.mesh.MeshFunctionSizet( mesh, mvc )
xdmf.close()
#read the lines
mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim()-1)
with XDMFFile("line_mesh.xdmf") as infile:
infile.read(mvc, "name_to_read")
cf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
xdmf.close()
#read the vertices
mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim()-2)
with XDMFFile("vertex_mesh.xdmf") as infile:
infile.read(mvc, "name_to_read")
pf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
xdmf.close()
#analytical expression for a scalar function used to test the ds
class FunctionTestIntegral(UserExpression):
def eval(self, values, x):
c_test = [0.3, 0.76, 0]
r_test = 0.345
values[0] = np.cos(geo.my_norm(np.subtract(x, c_test)) - r_test)**2.0
def value_shape(self):
return (1,)
dx = Measure( "dx", domain=mesh, subdomain_data=vf, subdomain_id=6 )
ds_r = Measure( "ds", domain=mesh, subdomain_data=cf, subdomain_id=1 )
ds_R = Measure( "ds", domain=mesh, subdomain_data=cf, subdomain_id=2 )
ds_custom = Measure("ds", domain=mesh, subdomain_data=pf) # Point measure for points at the edges of the mesh
dS_custom = Measure("dS", domain=mesh, subdomain_data=pf) # Point measure for points in the mesh
Q = FunctionSpace( mesh, 'P', 1 )
f_test_ds = Function( Q )
f_test_ds.interpolate( FunctionTestIntegral( element=Q.ufl_element() ))
print(f"int f dx = {assemble(f_test_ds * dx)}")
print(f"int f ds_r = {assemble(f_test_ds * ds_r)}")
print(f"int f ds_p_1 = {assemble(f_test_ds * ds_custom)}")
To clarify, in the last line with I would like assemble(f_test_ds * ds_custom) return the value of f at the point p_1 previously defined.
When I run the latter code, I get the following error
Mesh dimension = 2
Traceback (most recent call last):
File "read_2dmesh_ring.py", line 31, in <module>
infile.read(mvc, "name_to_read")
RuntimeError:
*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
*** fenics-support@googlegroups.com
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error: Unable to find entity in map.
*** Reason: Error reading MeshValueCollection.
*** Where: This error was encountered inside HDF5File.cpp.
*** Process: 0
***
*** DOLFIN version: 2019.1.0
*** Git changeset: 74d7efe1e84d65e9433fd96c50f1d278fa3e3f3f
*** -------------------------------------------------------------------------
Do you know what I am doing wrong ? I tried searching online and asking ChatGPT before posting, nothing works.
