Thank you, it appears to work now! The vertex tagging is also clear, I tagged the two lines first (ids 1 and 2), then the two vertices at the boundary (ids 3 and 4) and finally the vertex in the middle (id 5).
When I read them out, all the integrals check out:
import numpy
import meshio
import gmsh
import pygmsh
import argparse
import numpy as np
parser = argparse.ArgumentParser()
parser.add_argument("resolution")
args = parser.parse_args()
#mesh resolution
resolution = (float)(args.resolution)
#mesh parameters
#CHANGE PARAMETERS HERE
L = 1.0
h = L
r = 1.0
c_r = [0, 0, 0]
#CHANGE PARAMETERS HERE
print("L = ", L)
print("h = ", h)
print("r = ", r)
print("c_r = ", c_r)
print("resolution = ", resolution)
# Initialize empty geometry using the build in kernel in GMSH
geometry = pygmsh.occ.Geometry()
# Fetch model we would like to add data to
model = geometry.__enter__()
#add a 1d object a set of lines
points = [model.add_point( (0, 0, 0), mesh_size=resolution ),
model.add_point((np.pi/8.0, 0, 0), mesh_size=resolution),
model.add_point((L, 0, 0), mesh_size=resolution)
]
my_lines = [model.add_line( points[0], points[1] ), model.add_line( points[1], points[2] )]
#add a 2d object: a plane surface starting from the 4 lines above
model.synchronize()
print("# of lines added = ", len(my_lines))
model.add_physical([my_lines[0]], "line1")
model.add_physical([my_lines[1]], "line2")
model.add_physical([points[0]], "point_l")
model.add_physical([points[2]], "point_r")
model.add_physical([points[1]], "point_in")
geometry.generate_mesh(dim=3)
gmsh.write("mesh.msh")
model.__exit__()
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 tetrahedron mesh
tetrahedron_mesh = create_mesh(msh, "tetra", True)
meshio.write("tetrahedron_mesh.xdmf", tetrahedron_mesh)
#create a triangle mesh
triangle_mesh = create_mesh(mesh_from_file, "triangle", prune_z=False)
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)
then I read the mesh with
from __future__ import print_function
from fenics import *
from mshr import *
from fenics import *
from mshr import *
import numpy as np
import ufl as ufl
import argparse
from dolfin import *
parser = argparse.ArgumentParser()
parser.add_argument("input_directory")
args = parser.parse_args()
#CHANGE PARAMETERS HERE
L = 1
h = 1
r = 1
c_r = [0, 0]
#CHANGE PARAMETERS HERE
#read the mesh
mesh = Mesh()
xdmf = XDMFFile(mesh.mpi_comm(), "line_mesh.xdmf")
xdmf.read(mesh)
#read the lines
mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim())
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()-1)
with XDMFFile("vertex_mesh.xdmf") as infile:
infile.read(mvc, "name_to_read")
sf = 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):
values[0] = (np.cos(3+x[0]))**2
def value_shape(self):
return (1,)
dv_custom = Measure("dx", domain=mesh, subdomain_data=cf) # Line measure
ds_custom = Measure("ds", domain=mesh, subdomain_data=sf) # Point measure for points at the edges of the mesh
dS_custom = Measure("dS", domain=mesh, subdomain_data=sf) # Point measure for points in the mesh
Q = FunctionSpace( mesh, 'P', 1 )
# f_test_ds is a scalar function defined on the mesh, that will be used to test whether the boundary elements ds_circle, ds_inflow, ds_outflow, .. are defined correclty . This will be done by computing an integral of f_test_ds over these boundary terms and comparing with the exact result
f_test_ds = Function( Q )
f_test_ds.interpolate( FunctionTestIntegral( element=Q.ufl_element() ))
#print out the integrals on the surface elements and compare them with the exact values to double check that the elements are tagged correctly
print(f"Volume = {assemble(Constant(1.0)*dv_custom)}, should be 1.0")
print(f"Integral over the whole domain = {assemble( f_test_ds * dv_custom )}", " should be 0.727324")
print(f"Integral over line #1 = {assemble( f_test_ds * dv_custom(1) )}", "should be 0.373126")
print(f"Integral over line #2 = {assemble( f_test_ds * dv_custom(2) )}", "should be 0.354198")
#this computes \sum_{i \in vertices in ds_custom} f_test_ds (i-th vertex in ds_custom)
print(f"Integral over point_l = {assemble( f_test_ds * ds_custom(3) )} should be 0.980085")
print(f"Integral over point_r = {assemble( f_test_ds * ds_custom(4) )} should be 0.42725")
print(f"Integral over point_in = {assemble( f_test_ds * dS_custom(5) )} should be 0.93826")
which gives
Volume = 1.0, should be 1.0
Integral over the whole domain = 0.8171933306687675 should be 0.817193
Integral over line #1 = 0.38654493390746114 should be 0.386545
Integral over line #2 = 0.43064839676130673 should be 0.430648
Integral over point_l = 0.9800851433251829 should be 0.980085
Integral over point_r = 0.4272499830956933 should be 0.42725
Integral over point_in = 0.9382597571646916 should be 0.93826
I hope that this thread will be useful to other Fenics users 