Wrong assignation of values in dolfinx vectors

variational formulation
1

Dear FEniCSx community,

I am solving a problem using mixed elements. Once I get the solution vector, I need to change the sign of the components associated with the first subspace and set to zero the components associated with the second subspace (the reason is that this new vector will be the initial condition of another problem). I am trying to achieve this by defining:

q_n = fem.Function(Y)
q_n.sub(0).x.array[:] = -qh.sub(0).x.array 
q_n.sub(1).x.array[:] = qh.sub(1).x.array*0

where Y is a mixed function space, qh is the solution vector of the problem and q_n is the new vector that I am trying to define. However, when I do this, q_n is a vector of all zeros.

Please, find below a minimal code of my problem. To make things easier, Iā€™m just solving a the same Poisson problem at each function space:

from dolfinx import fem, mesh, io, plot, cpp, graph
import ufl
import numpy as np
from mpi4py import MPI
rank = MPI.COMM_WORLD.rank

# Define mesh
domain = mesh.create_unit_square(MPI.COMM_WORLD, 8, 8, mesh.CellType.quadrilateral)

# Define function space
el1 = ufl.FiniteElement("Lagrange", domain.ufl_cell(), 1)
el2 = ufl.FiniteElement("Lagrange", domain.ufl_cell(), 1)
mel = ufl.MixedElement([el1, el2])
Y = fem.FunctionSpace(domain, mel)

# Get subspaces
num_subs = Y.num_sub_spaces
spaces = []
maps = []
for i in range(num_subs):
    space_i, map_i = Y.sub(i).collapse()
    spaces.append(space_i)
    maps.append(map_i)
V = spaces[0]
Q = spaces[1]

import numpy
# Create facet to cell connectivity required to determine boundary facets
tdim = domain.topology.dim
fdim = tdim - 1
domain.topology.create_connectivity(fdim, tdim)
boundary_facets = numpy.flatnonzero(mesh.compute_boundary_facets(domain.topology))

# Define solution vector
qh = fem.Function(Y)
uh = qh.sub(0)
ph = qh.sub(1)

# Finite Element problem
u, p = ufl.TrialFunctions(Y)
v, q = ufl.TestFunctions(Y)

from petsc4py.PETSc import ScalarType
f = fem.Constant(domain, ScalarType(-6))
a = ufl.dot(ufl.grad(u), ufl.grad(v)) * ufl.dx + ufl.dot(ufl.grad(p), ufl.grad(q)) * ufl.dx
L = f * v * ufl.dx + f * q * ufl.dx

# Define boundary conditions
uD = fem.Function(Y).sub(0)
uD.interpolate(lambda x: 1 + x[0]**2 + 2 * x[1]**2)
boundary_dofs_1 = fem.locate_dofs_topological(Y.sub(0), fdim, boundary_facets)
boundary_dofs_2 = fem.locate_dofs_topological(Y.sub(1), fdim, boundary_facets)
bc_1 = fem.dirichletbc(uD, boundary_dofs_1)
bc_2 = fem.dirichletbc(uD, boundary_dofs_2)
bc = [bc_1, bc_2]

# Solve variational problem
problem = fem.petsc.LinearProblem(a, L, bcs=bc, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
qh = problem.solve()

# Now, I want to change the solution values
q_n = fem.Function(Y)
q_n.sub(0).x.array[:] = -qh.sub(0).x.array 
q_n.sub(1).x.array[:] = qh.sub(1).x.array*0

print(qh.x.array[:])
print(q_n.x.array[:])

you can see that q_n is a vector of full zeros, when it should have the components associated to the first subspace equal to minus the value of the solution vector qh in the first subspace.

Does anybody know what I am doing wrong?

Thank you very much in advance!

2

Check the size of this array. It will be the same as q_n.x.array
This is because q_n.sub(0) is a view of the original function, and does not have its own array.

You should collapse the function space, and use the map to assign data from a collapsed function to the Mixed function, see: Assigning a scalar expression to mixed function space - #4 by dokken

3 Likes
3

Thank you very much @dokken for your help. It works.
Below is the working code, in case that anybody else needs it.

from dolfinx import fem, mesh, io, plot, cpp, graph
import ufl
import numpy as np
from mpi4py import MPI
rank = MPI.COMM_WORLD.rank

# Define mesh
domain = mesh.create_unit_square(MPI.COMM_WORLD, 8, 8, mesh.CellType.quadrilateral)

# Define function space
el1 = ufl.FiniteElement("Lagrange", domain.ufl_cell(), 1)
el2 = ufl.FiniteElement("Lagrange", domain.ufl_cell(), 1)
mel = ufl.MixedElement([el1, el2])
Y = fem.FunctionSpace(domain, mel)

# Get subspaces
num_subs = Y.num_sub_spaces
spaces = []
maps = []
for i in range(num_subs):
    space_i, map_i = Y.sub(i).collapse()
    spaces.append(space_i)
    maps.append(map_i)
V = spaces[0]
Q = spaces[1]

import numpy
# Create facet to cell connectivity required to determine boundary facets
tdim = domain.topology.dim
fdim = tdim - 1
domain.topology.create_connectivity(fdim, tdim)
boundary_facets = numpy.flatnonzero(mesh.compute_boundary_facets(domain.topology))

# Define solution vector
qh = fem.Function(Y)
uh = qh.sub(0)
ph = qh.sub(1)

# Finite Element problem
u, p = ufl.TrialFunctions(Y)
v, q = ufl.TestFunctions(Y)

from petsc4py.PETSc import ScalarType
f = fem.Constant(domain, ScalarType(-6))
a = ufl.dot(ufl.grad(u), ufl.grad(v)) * ufl.dx + ufl.dot(ufl.grad(p), ufl.grad(q)) * ufl.dx
L = f * v * ufl.dx + f * q * ufl.dx

# Define boundary conditions
uD = fem.Function(Y).sub(0)
uD.interpolate(lambda x: 1 + x[0]**2 + 2 * x[1]**2)
boundary_dofs_1 = fem.locate_dofs_topological(Y.sub(0), fdim, boundary_facets)
boundary_dofs_2 = fem.locate_dofs_topological(Y.sub(1), fdim, boundary_facets)
bc_1 = fem.dirichletbc(uD, boundary_dofs_1)
bc_2 = fem.dirichletbc(uD, boundary_dofs_2)
bc = [bc_1, bc_2]

# Solve variational problem
problem = fem.petsc.LinearProblem(a, L, bcs=bc, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
qh = problem.solve()

# Now, I want to change the solution values
V0, Vmap = Y.sub(0).collapse()
Q0, Qmap = Y.sub(1).collapse()

q_n = fem.Function(Y)
q_n.x.array[Vmap] = -qh.sub(0).collapse().x.array
q_n.x.array[Qmap] = qh.sub(1).collapse().x.array*0

print(qh.x.array[:])
print(q_n.x.array[:])
1 Like