Application of point forces, mapping vertex indices to corresponding DOFs

dolfinx
1

I am working in dolfinx with the 3D Reissner-Mindlin shell implementation here: GitHub - RuruX/shell_analysis_fenicsx

I am applying loads as independent point forces at a collection of vertices (I realize there are probably other options here, but it is important for me to stick with this approach presently). So far, it has proved simplest to directly insert the point forces into the force vector of my linear problem, i.e.:

A = assemble_matrix(form(a), bcs)
A.assemble()
b = assemble_vector(form(L))    #just used to initialize with proper sizing
b.setArray(f_node.vector.getArray())   #insert point forces here!

In order to construct f_node I of course must map my forces (which are vertex-based) to the corresponding degree of freedom. My current approach for this–which I would like to improve–is to do a brute force search between the mesh coordinates and the DOF coordinates to build a vtx_to_dof index map, i.e.:

# map vertex values to DOFs 
w0_coll, submap = W.sub(0).collapse()
dof_coords = w0_coll.tabulate_dof_coordinates()
msh_coords = mesh.geometry.x
submap = np.reshape(submap, (-1,3))
vtx_to_dof = np.zeros(np.shape(msh_coords))
for ii in range(np.shape(msh_coords)[0]):
    d = scipy.spatial.distance.cdist([msh_coords[ii,:]], dof_coords)
    assert(np.min(d) < 1.0e-6)
    idx = np.argmin(d)
    vtx_to_dof[ii,:] = submap[idx,:]
vtx_to_dof = np.reshape(vtx_to_dof, (-1,1))

This yields the desired result but, as you can imagine, is not particularly efficient. For the typical problems I work with this procedure is consuming nearly half the total computation time, and for larger problems it scales very poorly.

I suspect there is a cleaner and more efficient way either of applying the forces or, more likely, of constructing/retrieving this vtx_to_dof map. Can anybody point me in the right direction? I have seen various references to the vertex_to_dof_map (e.g. here) but from what I can tell this doesn’t seem to be implemented the same way in dolfinx (happy to be corrected if that isn’t the case though!).

2

Consider the following code, that can map from any mesh node to any Lagrange space:

import ufl
import dolfinx
from mpi4py import MPI
import numpy as np
import scipy


def f(x):
    return x[0] + 3 * x[1]


N = 4
P = 1
Q = 2
mesh = dolfinx.mesh.create_unit_cube(MPI.COMM_WORLD, N, N, N)

el0 = ufl.VectorElement("Lagrange", mesh.ufl_cell(), P)
el1 = ufl.FiniteElement("Lagrange", mesh.ufl_cell(), Q)
el = ufl.MixedElement([el0, el1])
V = dolfinx.fem.FunctionSpace(mesh, el)
u = dolfinx.fem.Function(V)
V0, V0_to_V = V.sub(0).collapse()
dof_layout = V0.dofmap.dof_layout

num_vertices = mesh.topology.index_map(
    0).size_local + mesh.topology.index_map(0).num_ghosts
vertex_to_dof_map = np.empty(num_vertices, dtype=np.int32)
num_cells = mesh.topology.index_map(
    mesh.topology.dim).size_local + mesh.topology.index_map(
    mesh.topology.dim).num_ghosts
c_to_v = mesh.topology.connectivity(mesh.topology.dim, 0)
for cell in range(num_cells):
    vertices = c_to_v.links(cell)
    dofs = V0.dofmap.cell_dofs(cell)
    for i, vertex in enumerate(vertices):
        vertex_to_dof_map[vertex] = dofs[dof_layout.entity_dofs(0, i)]

geometry_indices = dolfinx.cpp.mesh.entities_to_geometry(
    mesh, 0, np.arange(num_vertices, dtype=np.int32), False)
x = mesh.geometry.x
bs = V0.dofmap.bs
for vertex, geom_index in enumerate(geometry_indices):
    dof = vertex_to_dof_map[vertex]
    for b in range(bs):
        parent_dof = V0_to_V[dof*bs+b]
        u.x.array[parent_dof] = x[geom_index, b]

u0 = u.sub(0).collapse()
with dolfinx.io.VTXWriter(mesh.comm, "u.bp", [u0]) as vtx:
    vtx.write(0.0)

The code is verified by using the dof_to_vertex_map to map the mesh vertices into a Pth order finite element space.

Note that these functions can be made faster by rewriting the code with pure numpy arrays and use numba to JIT compile them

1 Like
3

Thanks for the reply, @dokken, and sorry for the slow response here. I finally got the chance to look into this a little bit–this does seem like a promising direction and appears to be behaving much more efficiently than my previous solution! I do still have one snag though that perhaps you can help with. In my case I require a vertex_to_dof_map for a 3D VectorFunctionSpace, i.e. where you have a map with shape:

vertex_to_dof_map = np.empty(num_vertices, dtype=np.int32)

I require:

vertex_to_dof_map = np.empty(np.shape(mesh.geometry.x), dtype=np.int32)

I modified the routine you shared to the below…I’m not sure if I interpreted the arguments to the dof_layout.entity_dofs correctly:

    vertex_to_dof_map = np.zeros(np.shape(mesh.geometry.x), dtype=np.int32)
    num_cells = mesh.topology.index_map(
        mesh.topology.dim).size_local + mesh.topology.index_map(
        mesh.topology.dim).num_ghosts
    c_to_v = mesh.topology.connectivity(mesh.topology.dim, 0)
    for cell in range(num_cells):
        vertices = c_to_v.links(cell)
        dofs = W0.dofmap.cell_dofs(cell)
        for i, vertex in enumerate(vertices):
            vertex_to_dof_map[vertex,0] = dofs[dof_layout.entity_dofs(0, i)]
            vertex_to_dof_map[vertex,1] = dofs[dof_layout.entity_dofs(1, i)]
            vertex_to_dof_map[vertex,2] = dofs[dof_layout.entity_dofs(2, i)]   #<-error line

This gives me the following error on the last line marked above:

IndexError: vector::_M_range_check: __n (which is 1) >= this->size() (which is 1)

This leads me to believe that the first argument to dof_layout.entity_dofs may not, in fact, be the dimension as I have supposed here…

Any insights on what has gone wrong here and how I might achieve the mapping for a 3D VectorFunctionSpace?

4

As you can see in my code above, it uses a 3D vector function space

so without a minimal example I cannot help you any further

5

As you can see in my code above, it uses a 3D vector function space

My apologies for the misunderstanding–I was looking for the VectorFunctionSpace command and did not see it. You surely know this, but for others’ reference if nothing else I did not realize that the VectorFunctionSpace essentially just performs the VectorElement creation and the vanilla FunctionSpace call as you did above. TIL.

As for a MWE, perhaps a minor modification to the MWE you provided will help to communicate what I’m after. Please see below, where I’ve noted the 2 lines that I modified and 2 lines that I added:

import ufl
import dolfinx
from mpi4py import MPI
import numpy as np
import scipy


def f(x):
    return x[0] + 3 * x[1]

N = 4
P = 1
Q = 2
mesh = dolfinx.mesh.create_unit_cube(MPI.COMM_WORLD, N, N, N)

el0 = ufl.VectorElement("Lagrange", mesh.ufl_cell(), P)
el1 = ufl.FiniteElement("Lagrange", mesh.ufl_cell(), Q)
el = ufl.MixedElement([el0, el1])
V = dolfinx.fem.FunctionSpace(mesh, el)
u = dolfinx.fem.Function(V)
V0, V0_to_V = V.sub(0).collapse()
dof_layout = V0.dofmap.dof_layout

num_vertices = mesh.topology.index_map(
    0).size_local + mesh.topology.index_map(0).num_ghosts
vertex_to_dof_map = np.empty((num_vertices,3), dtype=np.int32)   #modified line
num_cells = mesh.topology.index_map(
    mesh.topology.dim).size_local + mesh.topology.index_map(
    mesh.topology.dim).num_ghosts
c_to_v = mesh.topology.connectivity(mesh.topology.dim, 0)
for cell in range(num_cells):
    vertices = c_to_v.links(cell)
    dofs = V0.dofmap.cell_dofs(cell)
    for i, vertex in enumerate(vertices):
        vertex_to_dof_map[vertex,0] = dofs[dof_layout.entity_dofs(0, i)]   #modified line
        vertex_to_dof_map[vertex,1] = dofs[dof_layout.entity_dofs(1, i)]   #added line
        vertex_to_dof_map[vertex,2] = dofs[dof_layout.entity_dofs(2, i)]   #added line

geometry_indices = dolfinx.cpp.mesh.entities_to_geometry(
    mesh, 0, np.arange(num_vertices, dtype=np.int32), False)
x = mesh.geometry.x
bs = V0.dofmap.bs
for vertex, geom_index in enumerate(geometry_indices):
    dof = vertex_to_dof_map[vertex]
    for b in range(bs):
        parent_dof = V0_to_V[dof*bs+b]
        u.x.array[parent_dof] = x[geom_index, b]

u0 = u.sub(0).collapse()
with dolfinx.io.VTXWriter(mesh.comm, "u.bp", [u0]) as vtx:
    vtx.write(0.0)

This code errors on the vertex_to_dof_map[vertex,1] = dofs[dof_layout.entity_dofs(1, i)], I presume because–as I mentioned above–I don’t think I am properly understanding the arguments to entity_dofs.

I hope this makes it more clear what I am trying to achieve here. In the structural analysis context, for a given vertex, I would like to apply independent forces in each of the [X,Y,Z] directions, each of which should refer to a different DOF, right? The original MWE gives me a single DOF for each vertex whereas I would expect x3 DOFs for each vertex.

I have been digging into this and will continue to do so as I think I might have some fundamental misunderstanding here. Even looking at the line dofs = V0.dofmap.cell_dofs(cell), I am seeing only x4 DOFs per cell but I expect more than this. Not sure what I am missing…

6

Consider:

Here I am assigning a coordinate x=(x,y,z) using the vertex_to_dof_map. Note that I am using the dofmap block size to assign this data

This is because a vector function space uses this to go trhough every vertex. As you can see, the dof in dof_to_vertex_map is expanded by the block size (in this case 3), to make it possible to assign a vector to the vector function space.

7

Indeed, I had not fully grasped that portion of your code and wandered off in the wrong direction. Thank you for providing further explanation on some of these details. The fact that the “real” DOF index can be determined based on a “parent” DOF along with the dofmap block size makes sense in hindsight, but is not straightforward to infer IMO.

Regardless, to fully close the loop on my original code example, I’ve used the code from @dokken to update my code for constructing the vtx_to_dof map having the same structure (i.e. with shape [nVertices, 3]) as what I provided originally:

W0, W0_to_W = W.sub(0).collapse()
dof_layout = W0.dofmap.dof_layout

num_vertices = mesh.topology.index_map(
    0).size_local + mesh.topology.index_map(0).num_ghosts
vertex_to_par_dof_map = np.zeros(num_vertices, dtype=np.int32)
num_cells = mesh.topology.index_map(
    mesh.topology.dim).size_local + mesh.topology.index_map(
    mesh.topology.dim).num_ghosts
c_to_v = mesh.topology.connectivity(mesh.topology.dim, 0)
for cell in range(num_cells):
    vertices = c_to_v.links(cell)
    dofs = W0.dofmap.cell_dofs(cell)
    for i, vertex in enumerate(vertices):
        vertex_to_par_dof_map[vertex] = dofs[dof_layout.entity_dofs(0, i)]

geometry_indices = dolfinx.cpp.mesh.entities_to_geometry(
    mesh, 0, np.arange(num_vertices, dtype=np.int32), False)
bs = W0.dofmap.bs
vtx_to_dof = np.zeros((num_vertices,bs), dtype=np.int32)
for vertex, geom_index in enumerate(geometry_indices):
    par_dof = vertex_to_par_dof_map[vertex]
    for b in range(bs):
        vtx_to_dof[vertex, b] = W0_to_W[par_dof*bs+b]
vtx_to_dof = np.reshape(vtx_to_dof, (-1,1))

I’ve confirmed that vtx_to_dof here is identical vtx_to_dof in my original post, but, for my typical problem size, is now generated in ~2.5s instead of ~33s using the original method; this is a significant difference in my context. Thanks @dokken!

8

Here is a vectorized version of the code above (credit to @dokken)

def vertex_to_dof_map_original(mesh, V):
    num_vertices = (
        mesh.topology.index_map(0).size_local + mesh.topology.index_map(0).num_ghosts
    )
    vertex_to_dof_map = np.empty(num_vertices, dtype=np.int32)
    num_cells = (
        mesh.topology.index_map(mesh.topology.dim).size_local
        + mesh.topology.index_map(mesh.topology.dim).num_ghosts
    )
    c_to_v = mesh.topology.connectivity(mesh.topology.dim, 0)
    vertex_to_dof_map = np.empty(num_vertices, dtype=np.int32)
    num_cells = (
        mesh.topology.index_map(mesh.topology.dim).size_local
        + mesh.topology.index_map(mesh.topology.dim).num_ghosts
    )
    for cell in range(num_cells):
        vertices = c_to_v.links(cell)
        dofs = V.dofmap.cell_dofs(cell)
        for i, vertex in enumerate(vertices):
            vertex_to_dof_map[vertex] = dofs[V.dofmap.dof_layout.entity_dofs(0, i)]
    return vertex_to_dof_map


def vertex_to_dof_map_vectorized(mesh, V):
    num_vertices_per_cell = dolfinx.cpp.mesh.cell_num_entities(
        mesh.topology.cell_type, 0
    )

    dof_layout2 = np.empty((num_vertices_per_cell,), dtype=np.int32)
    for i in range(num_vertices_per_cell):
        var = V.dofmap.dof_layout.entity_dofs(0, i)
        assert len(var) == 1
        dof_layout2[i] = var[0]

    num_vertices = (
        mesh.topology.index_map(0).size_local + mesh.topology.index_map(0).num_ghosts
    )

    c_to_v = mesh.topology.connectivity(mesh.topology.dim, 0)
    assert (
        c_to_v.offsets[1:] - c_to_v.offsets[:-1] == c_to_v.offsets[1]
    ).all(), "Single cell type supported"

    vertex_to_dof_map = np.empty(num_vertices, dtype=np.int32)
    vertex_to_dof_map[c_to_v.array] = V.dofmap.list[:, dof_layout2].reshape(-1)
    return vertex_to_dof_map

The output of the following program

import time
import dolfinx
  
mesh = dolfinx.mesh.create_unit_cube(MPI.COMM_WORLD, 40, 40, 40)

V = dolfinx.fem.functionspace(mesh, ("P", 2))
u = dolfinx.fem.Function(V)


t0 = time.perf_counter()
v2d_original = vertex_to_dof_map_original(mesh, V)
t1 = time.perf_counter()
print("Original: ", t1 - t0)
v2d_vectorized = vertex_to_dof_map_vectorized(mesh, V)
t3 = time.perf_counter()
print("Vectorized: ", t3 - t1)

np.testing.assert_allclose(v2d_original, v2d_vectorized)

is

Original:  4.220731544000955
Vectorized:  0.006541958000525483
1 Like