### Connect DOF of subdomain with DOF of global mesh (speed up)

196

views

2

Hello Everyone,

So the code I wrote below works but is extremely (unusable) slow. I create a 2-D boundary mesh (mesh B) from a sub domain from a 3-D (meshA) mesh to do calculations then apply these values at boundary DOF in mesh A. The code below is so slow it defeats the purpose of my code. I was hoping there was a more efficient way to connect the two domains DOF numbering possibly with a Fenics Function.

```
#define domain
x_min = y_min = z_min = 0.0
x_max = y_max = 1.0
z_max = 1.98/5.0
#define number of elements
nx = ny =10
nz = 5
#create mesh A
mesh = BoxMesh(Point(x_min, y_min, z_min), Point(x_max, y_max, z_max), nx, ny, nz)
coordinates = mesh.coordinates()
# define function space
V = FunctionSpace(mesh, 'P', 1)
# tabulate mesh A's DOF coordinates
dof_mesh = V.tabulate_dof_coordinates().reshape(V.dim(),mesh.geometry().dim())
# Create boundary markers
mesh1 = BoundaryMesh(mesh,'exterior')
class Boundary_top(SubDomain):
def inside(self, x, on_boundary):
tol = 1E-14
return on_boundary or near(x[2],z_max,tol)
# create Mesh B
subdom =Boundary_top()
subdomains = CellFunction('size_t', mesh1)
subdomains.set_all(0)
subdom.mark(subdomains, 1)
bmesh = SubMesh(mesh1, subdomains, 1)
h = []
# calculate needed info on 2-D boundary mesh cells (Mesh B)
for bci, cell in enumerate(cells(bmesh)):
#get vertex coordinates
h.append(cell.get_vertex_coordinates())
h = np.array(h)
n1 = h[:,0:3]
n2 = h[:,3:6]
n3 = h[:,6:9]
node1 = np.zeros(shape=(len(n1),))
node2 = np.zeros(shape=(len(n1),))
node3 = np.zeros(shape=(len(n1),))
# connect Mesh B's DOF numbering with Mesh A's DOF numbering
# this is the extremely slow part
for i in range(len(n1)):
for j in range(len(dof_mesh)):
if np.array_equal(n1[i],dof_mesh[j]):
node1[i] = j
elif np.array_equal(n2[i],dof_mesh[j]):
node2[i] = j
elif np.array_equal(n3[i],dof_mesh[j]):
node3[i] = j
# I now have connection of MESH A's DOFs corresponding to each cell in MEsh B
node1 = np.array(node1)
node2 = np.array(node2)
node3 = np.array(node3)
```

Any help would be greatly appreciated!

Cheers,

Terrence

Community: FEniCS Project

### 1 Answer

5

If it is a matter of mapping the dofs from one domain to another, one solution is to interpolate the dofs from the original domain into the subdomain.

For that, you have to define a new function space restricted to the subdomain.

For example:

For that, you have to define a new function space restricted to the subdomain.

For example:

```
#create mesh A
mesh = BoxMesh(Point(x_min, y_min, z_min), Point(x_max, y_max, z_max), nx, ny, nz)
# create Mesh B
subdom =Boundary_top()
subdomains = CellFunction('size_t', mesh1)
subdomains.set_all(0)
subdom.mark(subdomains, 1)
bmesh = SubMesh(mesh1, subdomains, 1)
# define function space
V = FunctionSpace(mesh, 'P', 1)
# define a function space restricted to subdomain B
Vb = FunctionSpace(bmesh, 'P', 1)
# Degrees of freedom
dofs = V.dofmap().dofs()
d = Function(V)
d.vector()[:] = dofs
#Interpolate d on Vb
mapping = interpolate(d,Vb)
# Map from the dofs of the subdomain to the global domain
print(mapping.vector()[:])
```

Hope it works for you.

This is great thank you!

written
4 months ago by
Terrence

Please login to add an answer/comment or follow this question.

https://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.spatial.KDTree.html

However, I haven't ever tried using this myself and am not sure of the caveats involved. There is some bounding box tree functionality in FEniCS, but another user reported issues with it for a similar task and eventually just used the scipy implementation instead:

https://www.allanswered.com/post/avepk/how-can-i-build-a-bounding-box-tree-using-a-point-cloud-in-fenics-2017-2/