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))