Dirichlet Boundary Conditions on Subdomains in FEniCSx

variational formulation
1

Hi everyone,

I am wondering if it is possible to define Dirichlet boundary conditions on a subdomain.

For instance, if one primary variable is not of interest in a certain subdomain, I would like to set a DBC with the value 0 in order to save computational time. In the legacy fenics, this was possible.

I appreciate any suggestions!

2

I guess you are referring to A.ident_zeros() in legacy dolfin?
As this simply calls petsc, you can still do this in DOLFINx, as illustrated in the following code:
(The following code is a refined snippet of @CastriMik code at thomas-fermi-demo/demo_tf_mwe.py at mpi-issue · mikics/thomas-fermi-demo · GitHub which avoid global communciation)

problem._A.zeroEntries()

fem.petsc.assemble_matrix(problem._A, problem._a, bcs=problem.bcs)
problem._A.assemble()

# Get diagonal of assembled A matrix
diagonal = problem._A.getDiagonal()
diagonal_values = diagonal.array

# Get zero rows of assembled A matrix.
zero_rows = problem._A.findZeroRows()
zero_rows_values_global = zero_rows.array
local_start = V.dofmap.index_map.local_range[0]*V.dofmap.index_map_bs

# Maps global numbering to local numbering
zero_rows_values_local = zero_rows_values_global - \
    local_start
diagonal.array[zero_rows_values_local] = np.ones(
    len(zero_rows_values_local), dtype=PETSc.ScalarType)
problem._A.setDiagonal(diagonal, PETSc.InsertMode.INSERT_VALUES)
problem._A.assemble()
3

You can also set DirichletBCs on any part of a subdomain in DOLFINx, as the input to the DirichletBC is a set of degrees of freedom.
I.e.

def some_locator_function(x):
    return x[0]<0.5
deactivate_cells = dolfinx.mesh.locate_entities(mesh, mesh.topology.dim, some_locator_function)
deactivate_dofs = dolfinx.fem.locate_dofs_topological(V, mesh.topology.dim, deactivate_cells)
bc = dolfinx.fem.dirichletbc(PETSc.ScalarType(0), deactivate_dofs, V)
1 Like
4

Thank you @dokken ! Both replies are very helpful.

5

Given V as a mixed function space, how can I find the corresponding coordinates of the DOF with the all zero rows?

My attempt is attached. Any suggestions on improvement?

for i in range(V.num_sub_spaces):
    v_sub, sub_to_parent = V.sub(i).collapse()
    
    if len(zero_rows_values_local) > 0:
        coords = v_sub.tabulate_dof_coordinates()
        for dof in zero_rows_values_local:
            if dof in sub_to_parent:
                index = sub_to_parent.index(dof)
                print(f"Found coordinates: {dof = }, {coords[index] = }", flush = True)
6

This can easily be done without a for-loop for the dof in zero_rows_values_local, if you create the inverse map parent_to_sub mapping from the parent space to the sub-space (where some entries will have nonsensical values, as there is nothing mapping from the parent to the sub-space).

for i in range(V.num_sub_spaces):
    v_sub, sub_to_parent = V.sub(i).collapse()
    num_parent_dofs = V.dofmap.index_map.size_local * V.dofmap.index_map_bs
    parent_to_sub = np.full(num_parent_dofs, -1, dtype=np.int32)
    parent_to_sub[sub_to_parent] = np.arange(len(sub_to_parent))
    sub_indices = parent_to_sub[zero_rows_values_local]
    assert np.allclose(sub_indices >= 0, 1)