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

196
views
2
4 months ago by

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
You might be able to use the scipy KD tree implementation here (e.g., with the query() method and a small upper bound distance, to find a matching point):

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:

written 4 months ago by David Kamensky

5
4 months ago by
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:
#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