Calculate area of marked subdomain

variational formulation
1

Hey,

Is it possible to calculate the area of the cell faces that are marked in a subdomain?

Below is some test code using a basic hex mesh cube domain. I want to extend this to a complex tetrahedral mesh, so if possible I’d like to avoid iterating through the marked faces and calculating their areas manually.

from dolfin import *
from dolfin_adjoint import *


n = 10
Nx, Ny, Nz = n,n,n

mesh = BoxMesh.create(MPI.comm_world, [Point(0,0,0), Point(1,1,1)], [Nx, Ny, Nz],
                          CellType.Type.hexahedron)
DG_FS = FunctionSpace(mesh, "DG", 0)
CG_FS = FunctionSpace(mesh, "CG", 1)
CG_VFS = VectorFunctionSpace(mesh, "CG", 1)

class top(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[1], 1.)

class bottom(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[1], 0)

top = top()
bottom = bottom()

materials = MeshFunction('size_t', mesh, 2)
materials.set_all(0)
top.mark(materials, 1)
bottom.mark(materials, 2)
ds = Measure("ds")(subdomain_data=materials)

# Total area of top
top_area = .... # Exact value is 1 in this case
bottom_area = .... # Exact value is 1 in this case

Thanks,
Tom

2

The following resolves the issue:

ds = Measure("ds",domain=mesh, subdomain_data=materials)

# Total area of top
top_area = assemble(1*ds(1)) # Exact value is 1 in this case
bottom_area =assemble(1*ds(2)) # Exact value is 1 in this case
print(top_area, bottom_area)

Please note that support for hexahedral meshes is very limited in dolfin, and has been improved and extended in dolfinx

2 Likes
3

Thanks for your quick response!

However, I’m getting a UFL error in assemble(1*ds(1)):

UFLException: This integral is missing an integration domain.

Do you know how I can fix this? Apologies if this is a trivial question, I’m quite new to FEniCS.

4

You need to define ds as I did,

Setting domain=mesh removes your error

5

Ah I didn’t notice that - work’s now! Thank you.

6

Is it possible to solve this problem in dolfinx? It seems, this solution doesn’t work for dolfinx

7

Could you elaborate, what exactly does not work for you? Could you provide a minimal example?

8

Sorry. For dolfinx I did the following

   one = Constant(mesh, ScalarType(1.))
   f = form(one*ds(2))
   area = assemble_scalar(f)

It works.