 Sum area of facets on surface that have midpoints within specific bounds

Hi everyone,

Currently, I have a UnitCubeMesh and I have marked the face where z = 1 using the following:

class marker(SubDomain):
def inside(self, x, on_boundary):
return near(x, 1) 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]

What I am aiming to do now is to select the facets on the z=1 surface with midpoint coordinates that satisfy x^2 + y^2 < r^2. Then I wish to sum up the area of these facets to get an approximate value of the area of a circle with radius r.

Does anyone know how this can be done? Here is what I have so far:

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

r = 0.3

class marker(SubDomain):
def inside(self, x, on_boundary):
return near(x, 1) 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)

mid_coordinates = [] #list of midpoint coordinates
for facet in surface_facets:
cell = all_cells[f_to_c(facet.index())]
local_facets = c_to_f(cell.index())
local_index = np.flatnonzero(local_facets == facet.index())
mid_coordinates.append(facet.midpoint().array())

area = []
for facet in surfacefacets:
cell = all_cells[f_to_c(facet.index())]
local_facets = c_to_f(cell.index())
local_index = np.flatnonzero(local_facets == facet.index())
area.append(cell.facet_area(local_index))

Consider the following:

from IPython import embed
from dolfin import *
import numpy as np

class marker(SubDomain):
def inside(self, x, on_boundary):
return near(x, 1) and on_boundary

def circle_area(r, N):
mesh = UnitCubeMesh(N, N, N)
tdim = mesh.topology().dim()
fdim = tdim - 1
mf = MeshFunction("size_t", mesh, fdim, 0)
value = 1
marker().mark(mf, value)
facet_indices = np.flatnonzero(mf.array() == value)
surface_facets = np.array(list(facets(mesh)))[facet_indices]

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

circ_cells = []
circ_index = []
circ_facets = []
for facet in surface_facets:
midpoint = facet.midpoint().array()
if midpoint**2 + midpoint**2 < r**2:
cell = f_to_c(facet.index())
local_facets = c_to_f(cell)
local_index = np.flatnonzero(local_facets == facet.index())
circ_cells.append(cell)
circ_index.append(local_index)
circ_facets.append(facet.index())
mf = MeshFunction("size_t", mesh, fdim, 0)
mf.array()[np.asarray(circ_facets)] = 1
File("mf.pvd") << mf
area = []
for cell, index in zip(circ_cells, circ_index):
c = Cell(mesh, cell)
area.append(c.facet_area(index))

print(N, sum(area), 1/4 * np.pi*r ** 2)

if __name__ == "__main__":
r = 0.3
for N in [10, 15, 20, 40]:
circle_area(r, N)

returning:

10 0.08000000000000002 0.07068583470577035
15 0.07111111111111111 0.07068583470577035
20 0.07249999999999998 0.07068583470577035
40 0.07062499999999997 0.07068583470577035
1 Like