Indices of the dofs of `CG=2` space associated with the vertices of the mesh

Question

Actually, I want to 'manually' add data to the paraview output files produced by dolfin.

As I understand, the paraview files contain the mesh topology (the vertices and how they are connected to triangles or tetraeders) and the data at the vertices.

If I have an FEM space of, say, quadratic order, there are not only dofs at the vertices but also at the edges.

How can I extract the dofs associated with the mesh vertices? (Such that I can append them to a vtu file that was written by dolfin)

In other words, I need the the array filter_vertices, in the pseudo-code below

import dolfin

N = 5

# setup the function space and write an example output

sqmesh = dolfin.UnitSquareMesh(N, N+1)
V = dolfin.FunctionSpace(sqmesh, 'CG', 2)
xyfunc = dolfin.Expression('x[0]+x[1]')
v = dolfin.interpolate(xyfunc, V)

vfile = dolfin.File('testit.pvd')
vfile << v   

# so far so good
# now I want to append the data on my own...

vvec = v.vector().array()

... # do some operations with vvec

vvec_at_vertices = vvec[filter_vertices]  # extract the dofs at the mesh vertices

... # append vvec_at_vertices to the "testit.pvd" file manually

PS.: I want to do this not getting back to dolfin, e.g. via v.set_local(vvec). So I really need this filter_vertices vector explicitly.

Answer 1

Hi, consider

from dolfin import * 
import numpy as np

mesh = UnitSquareMesh(20, 20)
V = FunctionSpace(mesh, 'CG', 2)
f = Expression('sin(x[0]*x[0])+cos(x[1]*x[1])', degree=2)
v = interpolate(f, V)
array = v.vector().array()

dofmap = V.dofmap()
nvertices = mesh.ufl_cell().num_vertices()
# For each cell these are indices of the cell dofs which correspond to vertices
# = entities of dim 0
indices = [dofmap.tabulate_entity_dofs(0, i)[0] for i in range(nvertices)]

vertex_2_dof = dict()
[vertex_2_dof.update(dict(vd for vd in zip(cell.entities(0),
                                           dofmap.cell_dofs(cell.index())[indices])))
 for cell in cells(mesh)]

# For check: the dof values at vertices match those extracted from the array
vertex_indices, vertex_dofs = map(list, zip(*vertex_2_dof.iteritems()))
X = mesh.coordinates()
X = X[vertex_indices]
x, y = X.transpose()
# Dof values
values0 = np.sin(x**2) + np.cos(y**2)
# Extracted
values = array[vertex_dofs]
# Crunch time
print np.linalg.norm(values0-values, np.inf)

Comments

It works! Thanks a lot.

---

could you tell me, how to manually add data. I want to add boundary point value with Dirichlet boundary condition.
Sincerely.

---

Paraview files have a section where the data corresponding to the vertices of the mesh is listed. I manually replace this section with my numbers...