Hello community,
I have the following setting: I intend to compute the volume by means of a surface integral over an enclosed boundary on a 3D body, which works.
Now, I have a body where this integration surface is open. So - I use a 2D Delauney triangulation to close the surface. Therefore, I have additional boundary elements (triangles) in my XDMF file, but no additional volume elements (tetrahedras) or nodes. Fine.
Now, if I perform the integration on this part of the surface (the “only triangle” one, meaning triangles that do not belong to a volume element), the value is always zero.
The following MWE reflects this (a 3D profile closed by 2D triangles at each end), where I just compute the area of a face.
#!/usr/bin/env python3
import numpy as np
from mpi4py import MPI
from dolfinx import Function, VectorFunctionSpace
from dolfinx.fem import assemble_scalar
from dolfinx.io import XDMFFile
from ufl import FacetNormal, ds, dot
comm = MPI.COMM_WORLD
# read in xdmf mesh - domain
with XDMFFile(comm,'domain.xdmf','r',encoding=XDMFFile.Encoding.ASCII) as infile:
mesh = infile.read_mesh(name="Grid")
# read in xdmf mesh - boundary
mesh.topology.create_connectivity(2, mesh.topology.dim)
with XDMFFile(comm,'boundary.xdmf','r',encoding=XDMFFile.Encoding.ASCII) as infile:
mt_s = infile.read_meshtags(mesh, name="Grid_surf")
# reference normal on mesh facets
n0 = FacetNormal(mesh)
V = VectorFunctionSpace(mesh, ("CG", 1))
u = Function(V)
one = Function(V)
one.vector.set(1.0)
# 8, 9 are the triangulated closures, should yield area 1
# 4, 5, 6, 7 are the outer faces, should yield 2.25 each
ds_ = ds(subdomain_data=mt_s, subdomain_id=9)
# surface form
surfaceform = -dot(one,n0)*ds_
surface = assemble_scalar(surfaceform)
print(surface)
The domain and boundary xdmf files are:
domain.xdmf
<Xdmf Version="3.0">
<Domain>
<Grid Name="Grid">
<Geometry GeometryType="XYZ">
<DataItem DataType="Float" Dimensions="24 3" Format="XML" Precision="8">
0.0 0.0 1.5
0.0 0.0 0.0
0.0 1.5 1.5
0.0 1.5 0.0
1.5 0.0 0.0
1.5 0.0 1.5
1.5 1.5 1.5
1.25 0.25 1.5
0.25 0.25 1.5
1.25 1.25 1.5
0.25 1.25 1.5
1.5 1.5 0.0
1.25 0.25 0.0
0.25 0.25 0.0
1.25 1.25 0.0
0.25 1.25 0.0
0.0 0.75 0.75
0.75 0.0 0.75
0.75 1.5 0.75
1.5 0.75 0.75
0.75 0.25 0.75
1.25 0.75 0.75
0.75 1.25 0.75
0.25 0.75 0.75
</DataItem>
</Geometry>
<Topology NodesPerElement="4" NumberOfElements="48" TopologyType="tetrahedron">
<DataItem DataType="Int" Dimensions="48 4" Format="XML" Precision="4">
8 17 0 5
15 1 3 16
10 18 6 2
8 16 2 0
11 3 18 14
5 6 9 19
17 13 1 4
12 19 11 4
7 21 12 5
2 15 23 10
20 12 5 7
11 14 9 21
12 21 11 19
12 11 21 14
2 23 8 10
2 8 23 16
5 20 8 17
20 5 8 7
5 21 9 7
21 5 9 19
13 0 17 1
13 17 12 4
13 12 17 20
5 12 17 4
5 17 12 20
15 16 2 23
15 2 16 3
15 16 13 1
15 13 16 23
0 13 16 1
9 19 11 21
9 11 19 6
11 9 18 6
10 9 18 22
10 18 9 6
2 15 18 3
14 15 18 22
14 18 15 3
12 19 5 21
12 5 19 4
9 18 14 11
18 9 14 22
16 13 8 23
16 8 13 0
17 13 8 0
17 8 13 20
10 15 18 2
10 18 15 22
</DataItem>
</Topology>
<Attribute AttributeType="Scalar" Center="Cell" Name="materials">
<DataItem DataType="Int" Dimensions="48" Format="XML" Precision="4">
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
</DataItem>
</Attribute>
</Grid>
</Domain>
</Xdmf>
boundary.xdmf
<Xdmf Version="3.0">
<Domain>
<Grid Name="Grid_surf">
<Geometry GeometryType="XYZ">
<DataItem DataType="Float" Dimensions="24 3" Format="XML" Precision="8">
0.0 0.0 1.5
0.0 0.0 0.0
0.0 1.5 1.5
0.0 1.5 0.0
1.5 0.0 0.0
1.5 0.0 1.5
1.5 1.5 1.5
1.25 0.25 1.5
0.25 0.25 1.5
1.25 1.25 1.5
0.25 1.25 1.5
1.5 1.5 0.0
1.25 0.25 0.0
0.25 0.25 0.0
1.25 1.25 0.0
0.25 1.25 0.0
0.0 0.75 0.75
0.75 0.0 0.75
0.75 1.5 0.75
1.5 0.75 0.75
0.75 0.25 0.75
1.25 0.75 0.75
0.75 1.25 0.75
0.25 0.75 0.75
</DataItem>
</Geometry>
<Topology NodesPerElement="3" NumberOfElements="52" TopologyType="triangle">
<DataItem DataType="Int" Dimensions="52 3" Format="XML" Precision="4">
2 0 8
0 5 8
2 10 6
2 8 10
9 5 6
8 5 7
7 5 9
10 9 6
3 15 1
1 13 4
1 15 13
3 11 14
3 14 15
12 11 4
13 12 4
12 14 11
8 23 10
13 23 8
10 23 15
15 23 13
7 20 8
12 20 7
8 20 13
13 20 12
7 9 21
12 7 21
9 14 21
14 12 21
9 10 22
14 9 22
10 15 22
15 14 22
1 0 16
0 2 16
3 1 16
2 3 16
0 1 17
5 0 17
1 4 17
4 5 17
2 18 3
6 18 2
3 18 11
11 18 6
4 19 5
11 19 4
5 19 6
6 19 11
10 8 9
7 8 9
14 13 12
15 13 14
</DataItem>
</Topology>
<Attribute AttributeType="Scalar" Center="Cell" Name="boundaries_surf">
<DataItem DataType="Int" Dimensions="52" Format="XML" Precision="4">
1
1
1
1
1
1
1
1
2
2
2
2
2
2
2
2
3
3
3
3
3
3
3
3
3
3
3
3
3
3
3
3
4
4
4
4
5
5
5
5
6
6
6
6
7
7
7
7
8
8
9
9
</DataItem>
</Attribute>
</Grid>
</Domain>
</Xdmf>
Setting subdomain_id in ds_ to 4, 5, 6, or 7 (outer faces), we each get the desired area of 2.25. Setting subdomain_id to 8 or 9 yields 0, even though these surfaces should give an area of 1.0.
Does someone have an idea why this is not working? Is there some check internally whether a facet belongs to a volume element, and if not, it is ignored?
I’d appreciate any advice!
Best,
Marc