I am still completely stuck, even with the much simpler Gmsh (and dolfin) codes.
First, the very simple gmsh code:
SetFactory("OpenCASCADE");
Box(1) = {0, 0, 0, 4, 4, 3};
Box(2) = {0, 0, 3, 4, 4, 1};
v() = BooleanFragments { Volume{1}; Delete; }{ Volume{2}; Delete;};
Physical Volume("top_block") = {#v()};
Physical Volume("bottom_block") = {#v()-1};
Physical Surface("bottom") = {5};
Physical Surface("top") = {11};
Physical Surface("intermediate") = {6};
where I have used to GUI to mark the surfaces of interest as Physical Surface(). So far so good, gmsh file.geo -3 returns no error and the generated mesh looks fine to me.
Here’s the Python code:
from dolfin import *
import meshio
msh = meshio.read("simple_gmsh_mesh.msh")
for key in msh.cell_data_dict["gmsh:physical"].keys():
if key == "triangle":
triangle_data = msh.cell_data_dict["gmsh:physical"][key]
elif key == "tetra":
tetra_data = msh.cell_data_dict["gmsh:physical"][key]
for cell in msh.cells:
if cell.type == "tetra":
tetra_cells = cell.data
elif cell.type == "triangle":
triangle_cells = cell.data
tetra_mesh = meshio.Mesh(points=msh.points, cells={"tetra": tetra_cells},
cell_data={"name_to_read":[tetra_data]})
triangle_mesh =meshio.Mesh(points=msh.points,
cells=[("triangle", triangle_cells)],
cell_data={"name_to_read":[triangle_data]})
meshio.write("file2.xdmf", tetra_mesh)
meshio.write("mf2.xdmf", triangle_mesh)
mesh = Mesh()
with XDMFFile("file2.xdmf") as infile:
infile.read(mesh)
mvc = MeshValueCollection("size_t", mesh, 2)
with XDMFFile("mf2.xdmf") as infile:
infile.read(mvc, "name_to_read")
mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
mvc2 = MeshValueCollection("size_t", mesh, 3)
with XDMFFile("file2.xdmf") as infile:
infile.read(mvc2, "name_to_read")
cf = cpp.mesh.MeshFunctionSizet(mesh, mvc2)
dx = Measure("dx", domain=mesh,subdomain_data=cf)
ds = Measure("ds", domain=mesh, subdomain_data=mf)
subdomains = cf
boundary_parts = mf
print('test volume: ',assemble(Constant(1)*dx(domain=mesh, subdomain_data=cf)))
print('test volume: ',assemble(Constant(1)*dx(domain=mesh, subdomain_data=cf, subdomain_id=0)))
print('test volume: ',assemble(Constant(1)*dx(domain=mesh, subdomain_data=cf, subdomain_id=1)))
print('test volume: ',assemble(Constant(1)*dx(domain=mesh, subdomain_data=cf, subdomain_id=2)))
print('test volume: ',assemble(Constant(1)*dx(domain=mesh, subdomain_data=cf, subdomain_id=3)))
print('test volume: ',assemble(Constant(1)*dx(domain=mesh, subdomain_data=cf, subdomain_id=4)))
print('test volume: ',assemble(Constant(1)*dx(domain=mesh, subdomain_data=cf, subdomain_id=5)))
print('test surface: ',assemble(Constant(1)*ds(domain=mesh, subdomain_data=mf)))
print('test surface: ',assemble(Constant(1)*dS(domain=mesh, subdomain_data=mf)))
print('test surface: ',assemble(Constant(1)*dS(domain=mesh, subdomain_data=mf, subdomain_id=1)))
print('test surface: ',assemble(Constant(1)*dS(domain=mesh, subdomain_data=mf, subdomain_id=2)))
print('test surface: ',assemble(Constant(1)*dS(domain=mesh, subdomain_data=mf, subdomain_id=3)))
print('test surface: ',assemble(Constant(1)*dS(domain=mesh, subdomain_data=mf, subdomain_id=4)))
print('test surface: ',assemble(Constant(1)*dS(domain=mesh, subdomain_data=mf, subdomain_id=5)))
print('test surface: ',assemble(Constant(1)*dS(domain=mesh, subdomain_data=mf, subdomain_id=6)))
print('test surface: ',assemble(Constant(1)*dS(domain=mesh, subdomain_data=mf, subdomain_id=7)))
print('test surface: ',assemble(Constant(1)*ds(domain=mesh, subdomain_data=mf)))
print('test surface: ',assemble(Constant(1)*ds(domain=mesh, subdomain_data=mf, subdomain_id=1)))
print('test surface: ',assemble(Constant(1)*ds(domain=mesh, subdomain_data=mf, subdomain_id=2)))
print('test surface: ',assemble(Constant(1)*ds(domain=mesh, subdomain_data=mf, subdomain_id=3)))
print('test surface: ',assemble(Constant(1)*ds(domain=mesh, subdomain_data=mf, subdomain_id=4)))
print('test surface: ',assemble(Constant(1)*ds(domain=mesh, subdomain_data=mf, subdomain_id=5)))
print('test surface: ',assemble(Constant(1)*ds(domain=mesh, subdomain_data=mf, subdomain_id=6)))
print('test surface: ',assemble(Constant(1)*ds(domain=mesh, subdomain_data=mf, subdomain_id=7)))
print('test surface: ',assemble(Constant(1)*ds(domain=mesh, subdomain_data=mf, subdomain_id=11)))
The output of running this file is extremely strange to me:
test volume: 16.00000000000001
test volume: 0.0
test volume: 0.0
test volume: 16.00000000000001
test volume: 0.0
test volume: 0.0
test volume: 0.0
test surface: 48.00000000000001
test surface: 156.12965142809307
test surface: 0.0
test surface: 0.0
test surface: 0.0
test surface: 0.0
test surface: 0.0
test surface: 0.0
test surface: 0.0
test surface: 48.00000000000001
test surface: 0.0
test surface: 0.0
test surface: 15.999999999999993
test surface: 0.0
test surface: 0.0
test surface: 0.0
test surface: 0.0
test surface: 0.0
Which indicates that dx is indeed all the volume as it should, but dx(0) through dx(5) is no volume element, except for dx(2) which is also all the volume. I would have expected dx(1) to be the top block volume element and dx(2) to be the bottom block volume element, as is shown with Gmsh. What is going on? I have no idea.
Then external surfaces. ds() yields 48 which is 3 x 16, or, I assume, the 3 Physical Surface() I have defined. This means that the interface (internal surface where I want to spatially average u later on) is considered an external surface. I have no idea why. In particular ds(3) = 16 so it is one of the surfaces I have defined as Physical Surface(), but which one, I am not sure. It should be the bottom surface, but then again, if that was true then ds(4) and ds(5) should have returned 16 instead of 0. Thus, I don’t understand what’s going on.
Then dS() is worth 156, I assume it consists of many internal surfaces that I haven’t defined as Physical Surface(). I have no idea how to select the dS() corresponding to the interface I am interested in.
