How to create SubDomain for each "sub-cube" in an N x N x N BoxMesh?

7 weeks ago by
Ryan W  
I have created a regularly spaced BoxMesh, and would like to integrate the solution to an equation within each sub-cube of a BoxMesh.

In general, I have a box, which is composed of a number of sub-cubes. I also have a list of points, which are the centers of those cubes. I want to use the list of points to create a SubDomain for each "sub-cube", and have it select the two (or more, if higher resolution) elements for integration using the Measure() function later.

Here is the code I've been trying. The plot statement should show the elements which are in the subcube, but it always returns an empty plot, even when some of the members of sc_flag.array() are equal to 1.

from fenics import *
from mshr import *
import numpy as np

nx = 3
ny = 3
nz = 3

mesh = BoxMesh(Point(-1,-1,-1), Point(1,1,1), nx, ny, nz)

class SubCube(SubDomain):
    Grab all elements which correspond to the point at a cube's center

    def __init__(self,center): = center

    def inside(self, x, on_boundary):
        tol = DOLFIN_EPS
        res = 1.0
        x_flag = np.logical_and(x[0] >=[0]-res/2-tol, 
                                x[0] <=[0]+res/2+tol)
        y_flag = np.logical_and(x[1] >=[1]-res/2-tol, 
                                x[1] <=[1]+res/2+tol)
        z_flag = np.logical_and(x[2] >=[2]-res/2-tol, 
                                x[2] <=[2]+res/2+tol)
        return np.logical_and(z_flag, np.logical_and(x_flag, y_flag))

spacing = np.linspace(-0.5, 0.5, 3)

grid = np.meshgrid(spacing,spacing,spacing)
point_list = np.vstack(map(np.ravel, grid)).T

subcube = SubCube(point_list[0]) #create instance of subdomain
sc_flag = CellFunction('size_t', mesh) #create the CellFunction
sc_flag.set_all(0) #mark them all 0 intially
subcube.mark(sc_flag, 1) #mark elements within the cube with 1
plot(sc_flag) #should show elements which are in the subcube, but creates empty plot
Community: FEniCS Project
Please login to add an answer/comment or follow this question.

Similar posts:
Search »