Saving values in a VectorFunctionSpace function

Errors
1

I am getting the following error:

IndexError: Values to set must be a 1D array

The lines of code causing the error are these ones:

self.Vect = VectorFunctionSpace(mesh, ‘P’, 1)
q_function = Function(self.Vect)
q_function.vector()[:] = q

Isn’t it possible to use .vector()[:] with 2D arrays?

2

please read Read before posting: How do I get my question answered?. The code you provide is not reproducible, and not correctly formatted.

In general, the answer to

Isn’t it possible to use .vector()[:] with 2D arrays?

is no. It would be good to clarify why you would want it to be a 2D array. Is it because you interpret the first column of the 2D array is the first component of q, and the second column is the second component?

3

Yes, I want the first column of the 2D array to be the first component of q, and the second column to be the second component. Here is a MWE:

import numpy as np
from fenics import *
  
width = 1000 # Width of the mesh (um)
height = 1000 # Height of the mesh (um)
c_0 = 1.4e-05
dt = 0.1
dt_environment = 1
t = 0.0

mesh = RectangleMesh(Point(0, 0), Point(width,height), 10, 10)
W_D1 = FunctionSpace(mesh, 'P', 1)
Vect = VectorFunctionSpace(mesh, 'P', 1)

v = TestFunction(W_D1)
u = Function(W_D1)
u_old = Function(W_D1)
        
u_initial= f"{c_0}"
u_aux = Expression(u_initial,degree=1, L = width)
u_0 = interpolate(u_aux, W_D1)

u.assign(u_0)
u_old.assign(u_0)

q = np.zeros([mesh.num_vertices(),2])
for vertex in vertices(mesh):
    q[vertex_to_dof_map(W_D1)[vertex.index()],:] = np.random.uniform()
q_function = Function(Vect)

# q_function.vector()[:] = q
q_function = project(q,Vect) 

F_D1 = (
     inner(u - u_old, v) * dx
     + dt * q_function * grad(v) * dx
 )

# Solver
for time in np.arange(0, dt_environment, dt):
    set_log_level(30)
    solver_parameters = {
        "newton_solver": {
            "linear_solver": "gmres",
            "preconditioner": "ilu",
            "maximum_iterations": 1000,
            "relaxation_parameter": 0.7,
            "relative_tolerance": 1e-7,
            "absolute_tolerance": 1e-8
        }
    }
    
    solve(F_D1 == 0, u, solver_parameters=solver_parameters) 
    u.vector()[:] = np.maximum(u.vector().get_local(), 0)
        
    t += dt
    u_old.assign(u)
4

Consider the following approach:

from fenics import *


mesh = RectangleMesh(Point(0, 0), Point(width, height), 10, 10)
Vect = VectorFunctionSpace(mesh, "P", 1)
u = Function(Vect)

u1 = u.sub(1, deepcopy=True)
print(u1.vector().get_local())
# Modify u1
u1.vector()[:] = 2

# Create assigner from u1 to u
assigner = FunctionAssigner(Vect.sub(1), u1.function_space())
assigner.assign(u.sub(1), u1)

# Assign to u
print(u.vector().get_local())
File("u.pvd") << u