How to plot a vector solution (3D) on a 2D section of the domain?

Question

I have the following example code, but it works just for a scalar value, I want adapt it to a vector with 3 components. How can I do? Thanks!

from dolfin import *

mesh = BoxMesh (-4,-4,-4,4,4,4,10,10,10)

V = FunctionSpace (mesh, 'CG',1)

dt = 1 rho = 0.012 #giorni^(-1) G = 0.0013 #giorni^(-1) D = 0.0065 #cm^2 \ giorni tf = 60 #giorni

Initial condition

u_t0 = Expression ('exp(-(x[0]x[0])/(2sig0sig0)-(x[1]x[1])/(2sig0sig0)-(x[2]x[2])/(2sig0*sig0))',sig0=3, cell = tetrahedron)

u = TrialFunction (V) v = TestFunction (V) f = ut0 - 0.5*dt*(G-rho)*ut0 #OK a = (uv + dtD0.5inner(nablagrad(u),nablagrad(v)) + dt0.5(G-rho)uv)dx #OK L = fvdx - dt0.5Dinner(nablagrad(ut0),nabla_grad(v))*dx #OK

u = Function(V)

for i in range(1,tf): solve(a == L, u) if (i!=(tf-1)): ut0 = u f = ut0 - 0.5dt(G-rho)u_t0 L = fvdx - dt0.5Dinner(nablagrad(ut0),nabla_grad(v))*dx

Create mesh & corresponding function spaces

mesh2 = RectangleMesh(-4,-4,4,4,10,10) V2 = FunctionSpace (mesh2, 'CG',1) u2 = Function(V2) u = interpolate(u, V)

coords = mesh2.coordinates() dofmap = vertextodofmap(V2)

Pick your values of z

z = 1.0

for i in range(V2.dim()): u2.vector()[dof_map[i]] = u(coords[i][0],coords[i][1],z)

plot(u2) interactive()

Answer 1

Hi, consider

from dolfin import *
import numpy as np

# The requirement on input is that it is a vector valued function in 3d
# For simplicity it is taken as Expression here
class Foo(Expression):
    def eval(self, values, x):
        values[0] = x[0]*x[2]
        values[1] = x[1]*x[2]
        values[2] = 0.
    def value_shape(self):
        return (3, )

f = Foo()
z = 3.

mesh = UnitSquareMesh(10, 10)
V = VectorFunctionSpace(mesh, 'CG', 1)
x = mesh.coordinates().reshape((-1, 2))
v2d = vertex_to_dof_map(V)

values = np.zeros(V.dim(), dtype=float)
for vertex in range(x.shape[0]):
    y = f(x[vertex, 0], x[vertex, 1], z)
    dofs = [v2d[2*vertex+component] for component in range(2)]
    values[dofs] = y

f = Function(V)
f.vector()[:] = values

plot(f)
interactive()

Also have a look into the code if you want to understand why the above works.

Comments

Dear Mirok, thank you so much for your help. It was really helpful.