How to extract vertex values of *x,y coordinates* of vector-valued function?

Question

I am feeling extremely anxious and frustrated bc I cannot find the answer to what I think is a simple question?

I just want the x-coordinate vertex values of a vector-valued function so I can make a streamplot in matplot lib.

Specifically I want just the x-coordinates of the vertex values of sigma below in this MWE:

from dolfin import *
%matplotlib inline
# Create mesh
mesh = UnitSquareMesh(32, 32)

# Define function spaces and mixed (product) space
BDM = FunctionSpace(mesh, "BDM", 1)
DG = FunctionSpace(mesh, "DG", 0)

element_bdm = FiniteElement("BDM", triangle, 1)
element_dg  = FiniteElement("DG", triangle, 0)
mixed_element = element_bdm * element_dg
W = FunctionSpace(mesh,mixed_element)

# Define trial and test functions
(sigma, u) = TrialFunctions(W)
(tau, v) = TestFunctions(W)

# Define source function
f = Expression("10*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02)",degree=1)

# Define variational form
a = (dot(sigma, tau) + div(tau)*u + div(sigma)*v)*dx
L = - f*v*dx

# Define function G such that G \cdot n = g
class BoundarySource(Expression):
    def __init__(self, mesh, **kwargs):
        self.mesh = mesh
    def eval_cell(self, values, x, ufc_cell):
        cell = Cell(self.mesh, ufc_cell.index)
        n = cell.normal(ufc_cell.local_facet)
        g = sin(5*x[0])
        values[0] = g*n[0]
        values[1] = g*n[1]
    def value_shape(self):
        return (2,)

G = BoundarySource(mesh,degree=1)

# Define essential boundary
def boundary(x):
    return x[1] < DOLFIN_EPS or x[1] > 1.0 - DOLFIN_EPS

bc = DirichletBC(W.sub(0), G, boundary)

# Compute solution
w = Function(W)
solve(a == L, w, bc)
(sigma, u) = w.split()

For scalar values I could just use compute_vertex_values(). But here we see that this concatenates the values into a single vector so I can't tell how to extract just the x-components.

coordinates = mesh.coordinates()
print coordinates.shape
print coordinates[1]
print(sigma(coordinates[1]))
vertex_values = sigma.compute_vertex_values()
print vertex_values.shape

in the Fenics manual supposedly using split should help, this gives me an error

x_sigma,y_sigma = sigma.split()

I also tried the hack on the previous question but you can't do uy.compute_vertex_values() as it throws an error also...please help if you can see at all where I am going wrong! Thank you!

Answer 1

Similar to the Nedlec elements I use, you cannot directly access the x- and y- components of a Function defined on a BDM space (maybe have a look at the definition of BDM elements).

Solution: Interpolate your sigma on a Lagrange VectroFunctionSpace:

.... your code until:
(sigma, u) = w.split()
V_cg = VectorFunctionSpace(mesh, 'CG', 1)
sigma_cg = interpolate(sigma, V_cg)

sigma_x_array = sigma_cg.vector.array()[::2]
sigma_y_array = sigma_cg.vector.array()[1::2]

These values belong to

V_cg.tabulate_dof_coordinates().reshape(-1,2)

Appendix: Interpolation errors are neglegtible. Try the following code snippet

import dolfin as df

mesh=df.UnitSquareMesh(1, 1)
V_cg = df.VectorFunctionSpace(mesh, 'CG', 1)
V_bdm = df.FunctionSpace(mesh, "BDM", 1)

u = df.Expression(('x[0]','x[1]'), element=V_cg.ufl_element())
u_bdm = df.interpolate(u, V)
u_cg_from_u_bdm = df.interpolate(u_bdm, V_cg)
u_cg_from_u_directly = df.interpolate(u, V_cg)

print(u_cg_from_u_bdm.vector().array() -
      u_cg_from_u_directly.vector().array())