Get components of a vector valued function using .sub()

dolfinx
1

Hello,

I have a quite naive question regarding the sub method for a vector valued function.
Let’s say for instance that i have a vector valued function u defined on a 2D domain,
i would like to get the x component of u on the boundary defined by x=0.
At first i thought that u.sub(0).vector.array would give the numpy array associated with u_x
especially i would expect both u_on_boundary and u_on_boundary_2 to give the same result in the MWE below:

from mpi4py import MPI
from ufl import VectorElement
from numpy import isclose

from dolfinx.fem import (FunctionSpace, locate_dofs_topological, dirichletbc, Constant, Function, set_bc)
from dolfinx.mesh import locate_entities_boundary, create_unit_square
from petsc4py.PETSc import ScalarType

# Mesh
mesh  = create_unit_square(MPI.COMM_WORLD, 1, 1)
ufl_cell = mesh.ufl_cell()
tdim = mesh.topology.dim

boundary_facets = locate_entities_boundary(mesh, dim=1, marker=lambda x: isclose(x[0], 0.0))
c = Constant(mesh, ScalarType(35))

# FunctionSpace
U_e = VectorElement('CG', ufl_cell, degree = 1)
V = FunctionSpace(mesh, U_e)

# Function
u = Function(V)

boundary_x_dofs = locate_dofs_topological(V.sub(0), tdim - 1, boundary_facets)
boundary_dofs = locate_dofs_topological(V, tdim - 1, boundary_facets)
bc_c = dirichletbc(c, boundary_x_dofs, V.sub(0))

#Set bc
set_bc(u.x.array, [bc_c])
print(u.vector.array)
print(u.sub(0).x.array)

#Value of displacement vector on boundary
u_on_boundary = u.x.array[boundary_x_dofs]
u_on_boundary_2 = u.sub(0).x.array[boundary_dofs]
print(u_on_boundary)
print(u_on_boundary_2)

Could someone explain me why it is not the case ?

Thanks for your time, sincerely,

Paul

2

Degrees of freedom of vector spaces are not unrolled when calling locate_dofs_topological.
This means that you have to “unroll” them for block size, i.e. for space 0:

from mpi4py import MPI
from ufl import VectorElement
from numpy import isclose

from dolfinx.fem import (FunctionSpace, locate_dofs_topological,
                         dirichletbc, Constant, Function, set_bc)
from dolfinx.mesh import locate_entities_boundary, create_unit_square
from petsc4py.PETSc import ScalarType

# Mesh
mesh = create_unit_square(MPI.COMM_WORLD, 1, 1)
ufl_cell = mesh.ufl_cell()
tdim = mesh.topology.dim

boundary_facets = locate_entities_boundary(
    mesh, dim=1, marker=lambda x: isclose(x[0], 0.0))


c = Constant(mesh, ScalarType(35))

# FunctionSpace
U_e = VectorElement('Lagrange', ufl_cell, degree=1)
V = FunctionSpace(mesh, U_e)

# Function
u = Function(V)

boundary_x_dofs = locate_dofs_topological(V.sub(0), tdim - 1, boundary_facets)


boundary_dofs = locate_dofs_topological(V, tdim - 1, boundary_facets)
boundary_dofs *= V.dofmap.index_map_bs

bc_c = dirichletbc(c, boundary_x_dofs, V.sub(0))
set_bc(u.x.array, [bc_c])


print(u.vector.array)
print(u.sub(0).x.array)

# Value of displacement vector on boundary
u_on_boundary = u.x.array[boundary_x_dofs]

u_on_boundary_2 = u.x.array[boundary_dofs]
u_on_boundary_3 = u.sub(0).x.array[boundary_dofs]

print(u_on_boundary)
print(u_on_boundary_2)
print(u_on_boundary_3)

Note that u.x.array and u.sub(0).x.array are the same, as u.sub(0) is not a deep copy of the function.
Similarly, for any other space than the 0th one, you can do:

sub_space = 1
boundary_x_dofs = locate_dofs_topological(
    V.sub(sub_space), tdim - 1, boundary_facets)


boundary_dofs = locate_dofs_topological(V, tdim - 1, boundary_facets)
boundary_dofs *= V.dofmap.index_map_bs
boundary_dofs += sub_space

bc_c = dirichletbc(c, boundary_x_dofs, V.sub(0))
set_bc(u.x.array, [bc_c])


print(u.vector.array)
print(u.sub(0).x.array)

# Value of displacement vector on boundary
u_on_boundary = u.x.array[boundary_x_dofs]

u_on_boundary_2 = u.x.array[boundary_dofs]
u_on_boundary_3 = u.sub(0).x.array[boundary_dofs]

print(u_on_boundary)
print(u_on_boundary_2)
print(u_on_boundary_3)
1 Like
3

That’s perfectly clear, thanks a lot !