Calculate area of marked subdomain

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

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

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.

You need to define ds as I did,

Setting domain=mesh removes your error

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

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

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

Sorry. For dolfinx I did the following

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

It works.