How to apply a boundary condition to a range of DOFs?

General
1

Hi all,

I’m trying to apply a boundary condition to a node set and the range of DOFs 1-6 as specified by my input .inp file like this:

*Boundary
Internal_Selection-1_Fixed-1, 1, 6, 0

I call locate_dofs_topological on the node set and that produces a list of marked DOFs. However, I get a bit stuck here because I don’t know how to apply the boundary condition to just the DOFs 1-6 as specified in the input file. What does that have to do with the marked DOFs I got from locate_dofs_topological? How do these pieces fit together?

If someone would be able to explain how to apply the BC to just that particular range of DOFs, that would really help me. I think I’m still trying to understand conceptually what is going on.

Thank you so much in advance.

Best regards,
Mia

2

Dear @mb246,

You have not put much effort into making this into a minimal reproducible example.

I am assuming that your dofs in range 1-6 (I am assuming that you are using 1-indexing) and that you are referring to degrees of freedom located at vertex 1-6 in the input mesh.

You need to map the input vertices to the dofs on your process.
See for instance:

and subsequent posts.

As clearly illustrated by the dolfinx demos, Dirichlet boundary conditions are applied by giving the indices of the dofs as input to dolfinx.fem.dirichletbc, you just need to map your global input vertex numbering to the local dof index on your process.

3

Thanks for the quick response.

I just wanted to clarify what you said here. Is this what is meant by “Internal_Selection-1_Fixed-1, 1, 6, 0” in the input file? That the “1, 6” is actually the vertex indices, not the degrees of freedom themselves? Thank you. I just want to make sure I’m interpreting the input file properly so I can move forward and apply the BCs correctly.

4

As you have specified your input file (the .inp file, which is not any format we natively support), I assume that you relate indices to indices of your input mesh. Thus you need to map from the input vertex index to the local degree of freedom ,as specified above.

5

Okay, I’ll try that, thanks. Yeah, I’ve been using an .inp file and converting it to xdmf using meshio.

6

Hi @dokken,

As a follow-up, here’s a similar example I found in the ccx_2.22 manual: https://www.dhondt.de/ccx_2.22.pdf#section.7

Specifically, my .inp file has a section that looks just like this:
*NSET, NSET=FIX
97, 96, 95, 94, 93, 20, 19, 18, 17, 16, 15,
14, 13, 12, 11, 10, 9, 4, 3, 2, 1
*BOUNDARY
FIX,1,3

Which as the example explains, means that degrees of freedom 1 through 3 are to be fixed (i.e. set to zero). Does this help provide more context for you regarding my issue? Specifically that I am provided a node set that the boundary condition needs to be applied to?

And just to be clear, even though nodes “3, 2, 1” appear in the node set, those are different from the “1,3” described in the boundary condition? The ccx manual didn’t clarify that for me enough so I wanted to make sure again.

So getting back to my own example, where the range is given as “1-6”, I wanted to make sure I understand your message earlier. So basically what I need to do is:

  1. Treat the nodes in the node set as vertices.
  2. Use the vertex-to-dofs map and the “1-6” range given in the input file to get the dofs corresponding to the vertices in the node set.
  3. What do I do now? The dirichlet BC function takes in a list of DOFs as one of its arguments. So do I pass that new list of DOFs (the ones located at vertices 1-6 in the input mesh) to the dirichlet BC function?

Thank you for your patience in answering my questions. I appreciate the help!

7

I don’t have time to read that user manual.

Assuming that what is in NSET describes the original nodes in the mesh you would like to constrain, you need to:

  1. Create a map from the input node index to a vertex on your process, this is simply a matter of inverting mesh.geometry.input_global_indices.
  2. Once you have a map from the global indices (in your case a map from 97, 96, 95, ... to what is local on the process), you remove all the results that are -1.
  3. You then map from geometry to topology indices, and then from topology to dof indices.
    I’ve illustrated the procedure below using a set of random input nodes (based on @Henrik_Nicolay_Finsb vertorized vertex_to_dofmap:
"""
Compute local dof indices of global input node indices
Author: Jørgen S. Dokken and Henrik N.T. Finsberg
SPDX-License-Identifier: MIT
"""

from mpi4py import MPI
import dolfinx
import numpy as np



# Compute map from vertex to dof
def vertex_to_dof_map_vectorized(V):
    mesh = V.mesh
    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




def input_indices_to_dof(V, global_dofs):
    mesh = V.mesh
    mesh.topology.create_connectivity(0, mesh.topology.dim)
    # Create a mapping from global to local node indices
    num_nodes_global = mesh.geometry.index_map().size_global
    local_to_input_global = mesh.geometry.input_global_indices
    global_to_local = np.full(num_nodes_global, -1, dtype=np.int32)
    global_to_local[local_to_input_global] = np.arange(len(local_to_input_global), dtype=np.int32)
    
    # Compute map from geometry node to vertex
    num_vertices_local = mesh.topology.index_map(0).size_local + mesh.topology.index_map(0).num_ghosts
    vertices_to_geometry = dolfinx.mesh.entities_to_geometry(mesh, 0, np.arange(num_vertices_local, dtype=np.int32)).reshape(-1)
    num_nodes_local = mesh.geometry.index_map().size_local + mesh.geometry.index_map().num_ghosts
    geometry_to_vertices = np.full(num_nodes_local, -1, dtype=np.int32)
    geometry_to_vertices[vertices_to_geometry] = np.arange(len(vertices_to_geometry), dtype=np.int32)
    
    # Compute map from vertex to dof
    v_to_d = vertex_to_dof_map_vectorized(V)

    # Remove all nodes that are not on process
    local_nodes = global_to_local[global_dofs]
    local_nodes = local_nodes[local_nodes >= 0]
    local_vertices = geometry_to_vertices[local_nodes]
    local_dofs = v_to_d[local_vertices]
    return local_dofs


if __name__ == "__main__":

    mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 10, 10)

    # Select a random set of input nodes for illustration
    np.random.seed(0)
    num_nodes_global = mesh.geometry.index_map().size_global

    input_nodes = np.random.randint(0, num_nodes_global, 10)
    print(MPI.COMM_WORLD.rank, input_nodes)

    V = dolfinx.fem.functionspace(mesh, ("Lagrange", 1))
    dofs = input_indices_to_dof(V, input_nodes)
    x_coord = V.tabulate_dof_coordinates()

    # Sanity check
    num_nodes_global = mesh.geometry.index_map().size_global
    local_to_input_global = mesh.geometry.input_global_indices
    global_to_local = np.full(num_nodes_global, -1, dtype=np.int32)
    global_to_local[local_to_input_global] = np.arange(len(local_to_input_global), dtype=np.int32)
    local_geom = global_to_local[input_nodes]
    local_geom = local_geom[local_geom >= 0]
    for local_node, dof in zip(local_geom, dofs):
        assert np.allclose(mesh.geometry.x[local_node], x_coord[dof])
        print(mesh.geometry.x[local_node,], x_coord[dof])
8

Hi @dokken, thank you, this helps clarify it more for me. Can I just confirm my steps with you, then?

# 6.1. Invert `mesh.geometry.input_global_indices` map to map input global node index to local vertex index 
# GLOBAL NODES -> LOCAL GEOMETRY
mesh = V.mesh
input_global_to_local = {v: k for k, v in enumerate(mesh.geometry.input_global_indices)}  

# 6.2. Remove the results that are -1
local_nset = [input_global_to_local[n] for n in nset if n in input_global_to_local]

# 6.3. Create a mapping from global to local node indices
mesh.topology.create_connectivity(0, mesh.topology.dim)
num_nodes_global = mesh.geometry.index_map().size_global
local_to_input_global = mesh.geometry.input_global_indices
global_to_local = np.full(num_nodes_global, -1, dtype=np.int32)
global_to_local[local_to_input_global] = np.arange(len(local_to_input_global), dtype=np.int32)

# 6.4. Compute map from geometry node to vertex
# LOCAL GEOMETRY -> VERTICES
num_vertices_local = mesh.topology.index_map(0).size_local + mesh.topology.index_map(0).num_ghosts
vertices_to_geometry = dolfinx.mesh.entities_to_geometry(mesh, 0, np.arange(num_vertices_local, dtype=np.int32)).reshape(-1)
num_nodes_local = mesh.geometry.index_map().size_local + mesh.geometry.index_map().num_ghosts
geometry_to_vertices = np.full(num_nodes_local, -1, dtype=np.int32)
geometry_to_vertices[vertices_to_geometry] = np.arange(len(vertices_to_geometry), dtype=np.int32)

# 6.5.
# Remove all nodes that are not on process
local_nodes = global_to_local[global_dofs]
local_nodes = local_nodes[local_nodes >= 0]
local_vertices = geometry_to_vertices[local_nodes]

# 6.6. Compute map from vertex to DOF
# LOCAL VERTICES -> LOCAL DOFS
v_to_d = vertex_to_dof_map_vectorized(V)

# 6.7. Get the local DOFs
local_dofs = v_to_d[local_vertices]

What I’m wondering about is, if we already have the local_nset from Step 6.2, then do we still need Steps 6.3-6.5, because isn’t it already mapped from global to local at that point? Or is that needed because we need to go from: GLOBAL NODES → LOCAL GEOMETRY → LOCAL VERTICES → LOCAL DOFs?

Also, can you please clarify what global_dofs should contain (in Step 6.5)? Should that be replaced by my local_nset variable from step 6.2? Because that needs to fit in somewhere.

Thank you.

9

Do I really need to re-read your code and compare it to the one I created above?
Why not simply use the code I supplied, i.e. the function:

There is a reason for me adding all the steps above, hich is that you need to go through:

Of course, if you have a higher order mesh, and want to map nodes from say a midpoint of a facet to the dof at a midpoint of a facet, something slightly different should be done.

"""
Compute local dof indices of global input node indices
Author: Jørgen S. Dokken and Henrik N.T. Finsberg
SPDX-License-Identifier: MIT
"""

from mpi4py import MPI
import dolfinx
import numpy as np

def input_indices_to_dof(V, global_dofs):
    mesh = V.mesh
    mesh.topology.create_connectivity(0, mesh.topology.dim)
    # Create a mapping from global to local node indices
    num_nodes_global = mesh.geometry.index_map().size_global
    local_to_input_global = mesh.geometry.input_global_indices
    global_to_local = np.full(num_nodes_global, -1, dtype=np.int32)
    global_to_local[local_to_input_global] = np.arange(len(local_to_input_global), dtype=np.int32)
    

    geom_layout =  mesh.geometry.cmap.create_dof_layout()
    V_layout = V.dofmap.dof_layout
    for i in range(mesh.topology.dim):
        num_entities_per_cell = dolfinx.cpp.mesh.cell_num_entities(
            mesh.topology.cell_type, i
    )
        for j in range(num_entities_per_cell):
            gl = geom_layout.entity_dofs(i, j)
            vl = V_layout.entity_dofs(i, j)
            assert vl == gl, "Mismatch in dof layout"
    geom_to_dof_layout = np.empty(V.dofmap.index_map.size_local+ V.dofmap.index_map.num_ghosts, dtype=np.int32)
    geom_to_dof_layout[mesh.geometry.dofmap.reshape(-1)] = V.dofmap.list.reshape(-1)


    # Remove all nodes that are not on process
    local_nodes = global_to_local[global_dofs]
    local_nodes = local_nodes[local_nodes >= 0]
    return geom_to_dof_layout[local_nodes]


if __name__ == "__main__":

    mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 10, 10)

    # Select a random set of input nodes for illustration
    np.random.seed(0)
    num_nodes_global = mesh.geometry.index_map().size_global

    input_nodes = np.random.randint(0, num_nodes_global, 10)
    print(MPI.COMM_WORLD.rank, input_nodes)

    V = dolfinx.fem.functionspace(mesh, ("Lagrange", 1))
    dofs = input_indices_to_dof(V, input_nodes)
    x_coord = V.tabulate_dof_coordinates()

    # Sanity check
    num_nodes_global = mesh.geometry.index_map().size_global
    local_to_input_global = mesh.geometry.input_global_indices
    global_to_local = np.full(num_nodes_global, -1, dtype=np.int32)
    global_to_local[local_to_input_global] = np.arange(len(local_to_input_global), dtype=np.int32)
    local_geom = global_to_local[input_nodes]
    local_geom = local_geom[local_geom >= 0]
    for local_node, dof in zip(local_geom, dofs):
        assert np.allclose(mesh.geometry.x[local_node], x_coord[dof])
        print(mesh.geometry.x[local_node,], x_coord[dof])

This will only work if you use the same finite element to describe the geometry of your mesh as you would use for the function space you want to store data in.

10

Okay, I ended up replacing global_dofs with local_nset in the workflow, as described in my previous post, so you can disregard this question:

I’ve appreciated your help so far, though your tone has been unnecessarily condescending at times (especially towards a newcomer to the forum) despite my efforts to improve the clarity of my questions and provide code where needed.

I will revisit this task later and will reach out again if anything else comes up.