Failed to split a vector

dolfinx
1

Hi, I tried to split a vector created by ufl.nabla_grad into its components but I failed to get the correct shape. I’m providing a minimal work example that’s adapted from the tutorial:

from mpi4py import MPI
from dolfinx import mesh
from dolfinx import fem
import numpy as np
import ufl
from petsc4py import PETSc
from dolfinx.fem.petsc import LinearProblem

domain = mesh.create_unit_square(MPI.COMM_WORLD, 8, 8, mesh.CellType.quadrilateral)
V = fem.FunctionSpace(domain, ("Lagrange", 1))
uD = fem.Function(V)
uD.interpolate(lambda x: 1 + x[0]**2 + 2 * x[1]**2)
tdim = domain.topology.dim
fdim = tdim - 1
domain.topology.create_connectivity(fdim, tdim)
boundary_facets = mesh.exterior_facet_indices(domain.topology)
boundary_dofs = fem.locate_dofs_topological(V, fdim, boundary_facets)
bc = fem.dirichletbc(uD, boundary_dofs)
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
default_scalar_type = PETSc.ScalarType
f = fem.Constant(domain, default_scalar_type(-6))
a = ufl.dot(ufl.grad(u), ufl.grad(v)) * ufl.dx
L = f * v * ufl.dx
problem = LinearProblem(a, L, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()
# solve nabla_uh
grad_uh = ufl.nabla_grad(uh) 
V_ = fem.VectorFunctionSpace(domain, ("Lagrange", 1)) 
dx_ = ufl.dx(V_.mesh)
u_ = ufl.TrialFunction(V_)
v_ = ufl.TestFunction(V_)
a_ = ufl.inner(u_,v_) * dx_
L_ =  ufl.inner(grad_uh,v_) * dx_
problem = fem.petsc.LinearProblem(a_, L_, bcs=[], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
grad_uh_value = problem.solve()
grad_uh_x, grad_uh_y = grad_uh_value.split()
print(uh.x.array.shape, grad_uh_value.x.array.shape, grad_uh_y.x.array.shape)

uh.x.array.shape gives (81,), and grad_uh_value.x.array.shape gives (162,) these are expected. However, grad_uh_x.x.array.shape also gives (162,) but I was expecting it to be (81,). Anyone knows how I can split the grad_uh_value correctly? Thanks

2

Read Read before posting: How do I get my question answered?. The code you posted is not reproducible, hence it’s unlikely anyone will be able to help.

3

Hi Francesco, thanks for the reply and sorry for making these confusions. I have provided code adapted from the fenicsx tutorial (and I confirm it can be run on my computer) and edited my question. Can you please take a look again? Thanks a lot

4

Hi.

When working with mixed formulations, like VectorFunctionSpace, you need to collapse the functions (see for example Meaning of .collapse() or here):

Taking

grad_uh_x=grad_uh_x.collapse()
grad_uh_y=grad_uh_y.collapse()
print(uh.x.array.shape, grad_uh_value.x.array.shape, grad_uh_y.x.array.shape)

gives

(81,) (162,) (81,)

as expected.

Cheers.

1 Like