Save the result of mpirun in a numpy array

I/O
4

Dear dokken,

Your solution worked perfectly.
Since I work with a vector function space (Current density vector field), I applied a minimal change to your solution as follows:

from dolfin import *
import numpy as np
mesh = UnitCubeMesh(MPI.comm_world,4,4,4)
V = VectorFunctionSpace(mesh, "DG", 0)
u = Function(V)
u.interpolate(Expression(('3*x[0]', '2*x[1]', 'x[2]*x[0]'), degree=2))

vertex_values0 = u.sub(0, deepcopy=True).compute_vertex_values()
vertex_values1 = u.sub(1, deepcopy=True).compute_vertex_values()
vertex_values2 = u.sub(2, deepcopy=True).compute_vertex_values()
coordinates = mesh.coordinates()


gathered_coordinates = MPI.comm_world.gather(coordinates, root=0)
gathered_values0 = MPI.comm_world.gather(vertex_values0, root=0)
gathered_values1 = MPI.comm_world.gather(vertex_values1, root=0)
gathered_values2 = MPI.comm_world.gather(vertex_values2, root=0)
global_indices = MPI.comm_world.gather(mesh.topology().global_indices(0))
if MPI.comm_world.rank == 0:
    num_vertices = mesh.num_entities_global(0)
    all_coordinates = np.zeros((num_vertices, coordinates.shape[1]))
    all_values = np.zeros((num_vertices,coordinates.shape[1]))
    for coord, value0, value1, value2, indices in zip(gathered_coordinates, gathered_values0, gathered_values1, gathered_values2, global_indices):
        all_values[indices] = np.transpose([value0, value1, value2])  # dedicate each vertex value to corresponding global index
        all_coordinates[indices] = coord

    for (coord, value) in zip(all_coordinates, all_values):
        print(coord, value)

Best wishes,
-Hassan

P.S. I got help from your answer in this topic:
https://fenicsproject.discourse.group/t/how-to-get-the-solution-obatained-by-fenics-as-a-vector/5377/7