Fenicsx: use an array of data for boundary conditions

5

This depends on how you do the interpolation into the appropriate space.
In the code above, i did not specify which cells to use interpolation for, i.e.

as opposed to

u_bc.interpolate(f, np.array([0,1,2], dtype=np.int32))

which would only interpolate over the cells with local index 0, 1 and 2.
In general, for the problem you are considering, where you have only data at a sub-set of the dofs in a cell, you should do something like the following:

from IPython import embed
import dolfinx
from mpi4py import MPI
import numpy as np
Nx = 5
Ny = 6
mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, Nx, Ny)

V = dolfinx.fem.FunctionSpace(mesh, ("Lagrange", 1))
mesh.topology.create_connectivity(mesh.topology.dim-1, mesh.topology.dim)
boundary_facets = dolfinx.mesh.exterior_facet_indices(mesh.topology)

u_bc = dolfinx.fem.Function(V)
mesh.topology.create_connectivity(mesh.topology.dim-1, mesh.topology.dim)
f_to_c = mesh.topology.connectivity(mesh.topology.dim-1, mesh.topology.dim)
mesh.topology.create_connectivity(mesh.topology.dim, mesh.topology.dim-1)
c_to_f = mesh.topology.connectivity(mesh.topology.dim, mesh.topology.dim-1)

dof_layout = V.dofmap.dof_layout
coords = V.tabulate_dof_coordinates()
num_dofs = 0
for facet in boundary_facets:
    # Find local index for each facet
    cells = f_to_c.links(facet)
    assert len(cells) == 1
    facets = c_to_f.links(cells[0])
    local_index = np.flatnonzero(facets == facet)
    # Find local index for dofs in cell
    closure_dofs = dof_layout.entity_closure_dofs(
        mesh.topology.dim-1,  local_index)
    cell_dofs = V.dofmap.cell_dofs(cells[0])
    for dof in closure_dofs:
        local_dof = cell_dofs[dof]
        dof_coordinate = coords[local_dof]
        print(dof, dof_coordinate)
        for b in range(V.dofmap.bs):
            num_dofs += 1
            u_bc.x.array[local_dof * V.dofmap.bs + b] = dof_coordinate[b]
print(num_dofs, 2 * 2 * Nx + 2 * 2 * Ny)
with dolfinx.io.VTXWriter(mesh.comm, "u.bp", u_bc) as vtx:
    vtx.write(0.0)

which you can see only evaluates 44 times in the example I’ve presented, which is equal to going through each cell with a boundary facet once, and only accessing the dofs located on the boundary.

1 Like