Marking exterior boundary Mesh

mesh
1

Hello all!!

I am using a MATLAB toolbox to generate my mesh for my research in fluid flow through porous rocks. From a stack of binary images I generate a volumetric mesh of the rock.

My problem comes from the exterior surface facets that enclose my volumetric Mesh. The toolbox that I use has the limitation of not being able to generate perfectly flat boundaries for my cuboid mesh. Because of that, I can not mark the boundaries with the usual workflow of:

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

Face1().mark(boundaries, 1)

Can I explicitly identify the Facets which compose each of the six faces of my volumetric mesh of the rock with some other workflow?
I know that I can compute the exterior boundary mesh with the function BoundaryMesh(mesh,'exterior')
but until now I was not able to mark all six faces.

Thanks in advance,

Rafael Ferrão

2

Hi, consider

from dolfin import *

mesh = UnitCubeMesh(5, 5, 5)

f = MeshFunction('size_t', mesh, mesh.topology().dim()-1, 0)
DomainBoundary().mark(f, 1)  # inside method is (x, on_boundary) -> on_boundary

# We get it right
for e in SubsetIterator(f, 1):
    x, y, z = e.midpoint().array()
    assert (near(x*(1-x), 0) or near(y*(1-y), 0) or near(z*(1-z), 0))

area = assemble(Constant(1)*ds(domain=mesh, subdomain_data=f, subdomain_id=1))
assert abs(area - 6) < 1E-13, area