Here’s the complete code:

```
import numpy as np
import matplotlib.pyplot as plt
from dolfinx import mesh, fem, io, plot, default_scalar_type
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
from mpi4py import MPI
from petsc4py import PETSc
from petsc4py.PETSc import ScalarType
import ufl
import gmsh
import os
from tqdm import tqdm
import csv
import sys
### (1a) define the 3D domain
coords_lower_left = [0.0, 0.0, 0.0]
Lx, Ly, Lz = 1.0, 1.0, 1.0
coords_upper_right = [Lx, Ly, Lz]
num_elements = [1, 1, 1]
domain = mesh.create_box(MPI.COMM_WORLD, [coords_lower_left, coords_upper_right], num_elements, mesh.CellType.hexahedron)
### (1b) set element types
element_u = ufl.VectorElement("CG", domain.ufl_cell(), degree=2) ## vector: displacement
element_p= ufl.FiniteElement("CG", domain.ufl_cell(), degree=1) ## scalar: pressure
element_mix = ufl.MixedElement(element_u, element_p)
V = fem.FunctionSpace(domain, element_mix)
V_disp = V.sub(0)
V_press = V.sub(1)
### displacement component in x direction
V_disp_x = V_disp.sub(0)
def front(x):
return np.isclose(x[0], 1, atol=1.0e-12)
fdim = domain.topology.dim - 1
value_zero= fem.Constant(domain, 0.1)
facet_front = mesh.locate_entities_boundary(domain, fdim, front)
dofs_front_x = fem.locate_dofs_topological(V_disp_x, fdim, facet_front)
bc_front_x = fem.dirichletbc(value_zero, dofs_front_x, V_disp_x)
bcs = [bc_front_x]
```

I tried “fem.locate_dofs_topological” in the above to get dofs.

For plot the figure in Paraview, the code is

```
##### plot boundary conditions to Paraview for mixed formulation
Mixed_Space = V
u_boundary = fem.Function(Mixed_Space, name='boundary')
fem.petsc.set_bc(u_boundary.vector, bcs)
u_boundary.x.scatter_forward()
with io.XDMFFile(domain.comm, "boundary/bc_geometrical.xdmf", "w") as xdmf:
xdmf.write_mesh(domain)
xdmf.write_function(u_boundary.sub(0))
```