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

35

views

0

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.

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):
self.center = center
SubDomain.__init__(self)
def inside(self, x, on_boundary):
tol = DOLFIN_EPS
res = 1.0
x_flag = np.logical_and(x[0] >= self.center[0]-res/2-tol,
x[0] <= self.center[0]+res/2+tol)
y_flag = np.logical_and(x[1] >= self.center[1]-res/2-tol,
x[1] <= self.center[1]+res/2+tol)
z_flag = np.logical_and(x[2] >= self.center[2]-res/2-tol,
x[2] <= self.center[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.