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[2], 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[2], 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())[0]]
    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())[0]]
	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[2], 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[0]**2 + midpoint[1]**2 < r**2:
            cell = f_to_c(facet.index())[0]
            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