Output facet area from tagged list of surface facets

I/O
1

Hi everyone,

Currently I have a mesh that I have imported externally from a .xml file in order to examine the surface facets of the mesh which i tagged in a group when making the mesh. Currently, I can output the midpoints of each facet and the normal vectors using:

tag = 1  #Tag for 'surface' group in .msh file
surfacefacetindices = np.where(boundaries.array() == tag)[0]
surfacefacets = np.array(list(facets(mesh)))[surfacefacetindices]
coordinates = []
normalvec = []
for f in surfacefacets:
	coordinates.append(f.midpoint().array())
	normalvec.append(f.normal().array())

What I would like to know is, is there anyway to output the area of each surface facet from this tagged list? Thanks in advance

Edit: Given my problem, it is difficult to replicate using an in built mesh, therefore the link to the required files are here: Sphere - Google Drive

The minimal example is shown below:

from dolfin import *
import numpy as np
from mpi4py import MPI
import matplotlib.pyplot as plt
from collections import defaultdict



data = 'sphere'
mesh = Mesh()
hdf = HDF5File(mesh.mpi_comm(), data+'.h5', 'r')
hdf.read(mesh, '/mesh', False)
cells_f = MeshFunction('size_t', mesh, 3)
hdf.read(cells_f, '/cells')
facets_f = MeshFunction('size_t', mesh, 2)
hdf.read(facets_f, '/facets')

boundaries = MeshFunction('size_t', mesh, 'sphere_facet_region.xml') 


Space = FunctionSpace(mesh,'P', 1)
x_Space = Space.tabulate_dof_coordinates().reshape((-1,mesh.geometry().dim()))
V = FunctionSpace(mesh,'P', 1)

tag = 1  #Tag for 'surface' group in .msh file
surfacefacetindices = np.where(boundaries.array() == tag)[0]
surfacefacets = np.array(list(facets(mesh)))[surfacefacetindices]


coordinates = []
normalvec = []
for f in surfacefacets:
	coordinates.append(f.midpoint().array())
	normalvec.append(f.normal().array())
2

Consider the following minimal example:

from dolfin import *
import numpy as np
mesh = UnitSquareMesh(10, 10)

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

mf = MeshFunction("size_t", mesh, mesh.topology().dim() - 1, 0)
value = 1
marker().mark(mf, value)
facet_indices = np.flatnonzero(mf.array() == value)
surface_facets = np.array(list(facets(mesh)))[facet_indices]

all_cells = np.array(list(cells(mesh)))

mesh.init(1,2)
f_to_c = mesh.topology()(1,2)
c_to_f = mesh.topology()(2,1)

for facet in surface_facets:
    cell = all_cells[f_to_c(facet.index())[0]]
    local_facets = c_to_f(cell.index())
    local_index = np.flatnonzero(local_facets == facet.index())
    print(facet.index(), cell.facet_area(local_index))

For further questions, please consider making a similar example to what I have presented above (using a built-in mesh) as it is then easier for others to contribute and reproduce your issue.

3

I have added a minimal working example and a link to the required files as I could not replicate it using an in built mesh. In addition, I am confused what the role of f_to_c and c_to_f are within this code

4

The code I posted above shows exactly what you would like to produce. It prints the area of each facet marked with the value 1. c_to_f maps a cell index to its facets, while f_to_c maps from a facet to the cell(s) that contains it.

5

If I extend your example code to 3D using UnitCubeMesh and then I try to use the same method to print the facet area for facets that lie on two different boundaries (which is something I am attempting to do in my much larger project), I get this error:

Traceback (most recent call last):
  File "simpletest.py", line 31, in <module>
    cell = all_cells[f_to_c(facet.index())[0]]
IndexError: index 6096 is out of bounds for axis 0 with size 6000

Any idea on how I can get this to successfully print the facet area for both boundaries in 3D? The working example is below:

from dolfin import *
import numpy as np
mesh = UnitCubeMesh(10, 10,10)

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

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

mf = MeshFunction("size_t", mesh, mesh.topology().dim() - 1, 0)
value = 1
value2 = 2
marker().mark(mf, value)
marker2().mark(mf, value2)

facet_indices_1 = np.flatnonzero(mf.array() == value)
facet_indices_2 = np.flatnonzero(mf.array() == value2)
facet_indices = np.append(facet_indices_1,facet_indices_2)
surface_facets = np.array(list(facets(mesh)))[facet_indices]

all_cells = np.array(list(cells(mesh)))

mesh.init(1,2)
f_to_c = mesh.topology()(1,2)
c_to_f = mesh.topology()(2,1)

for facet in surface_facets:
    cell = all_cells[f_to_c(facet.index())[0]]
    local_facets = c_to_f(cell.index())
    local_index = np.flatnonzero(local_facets == facet.index())
    print(facet.index(), cell.facet_area(local_index))
6

You should change this to

tdim = mesh.topology().dim()
fdim = tdim - 1
mesh.init(fdim, tdim)
f_to_c = mesh.topology()(fdim, tdim)
c_to_f = mesh.topology()(tdim, fdim)