Size mismatch Extracting vector and mesh points

1

Having computed my solution u to the stokes equations, I plotted it on a courser mesh using

```python
# plotting 

L = Constant(10)

mesh_plot = RectangleMesh(Point(0, 0), Point(1, L), 2, 5)

W2 = FunctionSpace(mesh_plot, V)
u_plot = interpolate(u, W2)

plot(u_plot)



However when I tried to extract the meshpoints and vector values to process in python the sizes of the 
meshpoints array and the vector values don't align



ur=np.array(u_plot.sub(0,deepcopy=True))
uz= np.array(u_plot.sub(1,deepcopy=True))

coor = np.array(mesh_plot.coordinates() )

print(np.size(coor))
print(np.size(uz))



How can I get the vector values evaluated at each of the meshpoints in the plot.
Thanks

Full code to reproduce



from dolfin import *
import numpy as np

L = Constant(10)
# Define mesh and 
mesh = RectangleMesh(Point(0, 0), Point(1, L), 100, 100)


boundaries = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
boundaries.set_all(0)


# Define Taylor--Hood function space W
V = VectorElement("CG", triangle, 2)
Q = FiniteElement("CG", triangle, 1)
W = FunctionSpace(mesh, MixedElement([V, Q]))

# Define Function and TestFunction(s)
w = Function(W)
(u, p) = split(w)
(v, q) = split(TestFunction(W))

# Boundary Conditions
u_out = Expression(('-x[0]*0.5*std::log(30)', '0'), degree=1)

# manually put in L
inl = 'near(x[1], 10)'
FB = 'near(x[0], 1.0)'
cent = 'near(x[0], 0.0)'
out = 'near(x[1], 0.0)'
bcu_inflow = DirichletBC(W.sub(0), (0.0,0.0),  inl)
bcu_cent = DirichletBC(W.sub(0).sub(0), (0.0), cent)  
bcu_outflow = DirichletBC(W.sub(0), u_out, out)
bcs = [bcu_inflow, bcu_cent, bcu_outflow]


# weak form
x = SpatialCoordinate(mesh)


def div_cyl(v):
    return (1/x[0])*(x[0]*v[0]).dx(0) + v[1].dx(1)



a1 = (inner( grad(u), grad(v) ) -  p*div_cyl(v) )  *dx
a2 = div_cyl(u)*q*dx    

F = a1 + a2


# Solve problem
solve(F == 0, w, bcs)


(u, p) = w.split()

# plotting 

L = Constant(10)

mesh_plot = RectangleMesh(Point(0, 0), Point(1, L), 2, 5)

W2 = FunctionSpace(mesh_plot, V)
u_plot = interpolate(u, W2)

plot(u_plot)


# mismatched sizes

ur=np.array(u_plot.sub(0,deepcopy=True))
uz= np.array(u_plot.sub(1,deepcopy=True))

coor = np.array(mesh_plot.coordinates() )

print(np.size(coor))
print(np.size(uz))
2

u_plot.sub(0) is the velocity approximation.
u_plot.sub(0).sub(0) is the radial component of the velocity approximation.
u_plot.sub(0).sub(1) is the z component of the velocity approximation.
u_plot.sub(1) is the pressure approximation.

Also, I think you’re missing the (\nabla\vec{u})_{\theta,\theta} component in your stress tensor.

3

Thanks for the reply I really appreciate it.

When I try

ur = np.array(u_plot.sub(0).sub(0))

I get the following error

RuntimeError                              Traceback (most recent call last)
<ipython-input-9-e86fae803159> in <module>
----> 1 ur = np.array(u_plot.sub(0).sub(0))
      2 uz= np.array(u_plot.sub(0).sub(1))
      3 
      4 coor = np.array(mesh_plot.coordinates() )
      5 

/usr/lib/petsc/lib/python3/dist-packages/dolfin/function/function.py in sub(self, i, deepcopy)
    542             raise RuntimeError("No subfunctions to extract")
    543         if not i < num_sub_spaces:
--> 544             raise RuntimeError("Can only extract subfunctions with i = 0..%d"
    545                                % num_sub_spaces)
    546 

RuntimeError: Can only extract subfunctions with i = 0..0

Does it work for you?

1 Like
4

I read your question too fast:

uz= u_plot.sub(1,deepcopy=True).vector().get_local()

will get you an array of the DoFs associated with the z component of your problem. But keep in mind you’re discretising with P2 Lagrange function space. This has DoFs at the vertices and the midpoints of facets.