Obtaining vertex values after loading a mesh

Question

After loading a mesh (.XML format), how may I extract the vertex coordinates and values without suffering interpolation errors.

This seems to work, I just want to be sure the values are direct, and not interpolated.

mesh = Mesh(meshFileName);
V = FunctionSpace(mesh);
x = Function(V, xFileName);

std::size_t count = mesh.num_vertices();

std::vector<double> vertex_values(count);
x.compute_vertex_values ( vertex_values ) ;

for(std::size_t i = 0; i < count; i++)
{
    double X = mesh.geometry().point(i).x();
    double Y = mesh.geometry().point(i).y();
    double f = vertex_values[i];
    // I store the point / value in a list (sort of a cache)
}

I then refine the mesh, and interpolate a function onto x again, where I use cached values of x from the load so I don't have to recalculate the function for the points already known (x is an expensive function).

My concern came with reading compute_vertex_values

void Function::compute_vertex_values(std::vector<double>& vertex_values,
                                     const Mesh& mesh) const
{
    // ....
       element.interpolate_vertex_values(&cell_vertex_values[0],
                                  &coefficients[0],
                                  cell_orientation,
                                  ufc_cell);
    // ...
}

Answer 1

You can access DOFs of vertex-based function using vertex indices by V->dofmap()->dof_to_vertex_map(mesh).

For DOFs located elsewhere use DofMap::tabulate_coordinates or DofMap::tabulate_all_coordinates.

Python example for CG1 space

from dolfin import *
mesh = UnitSquareMesh(2, 2)
coords = mesh.coordinates()

V = FunctionSpace(mesh, 'CG', 1)
vertex_to_dof_map = V.dofmap().vertex_to_dof_map(mesh)
dof_to_vertex_map = V.dofmap().dof_to_vertex_map(mesh)

u = Function(V)
x = u.vector()
dofs_at_vertices = x[dof_to_vertex_map]

# let's iterate over vertices
for v in vertices(mesh):
    print 'vertex index', v.index()
    print 'at point', v.point().str()
    print 'at coordinates', coords[v.index()]
    print 'dof', dofs_at_vertices[v.index()]

# let's iterate over dofs
for i in range(V.dim()):
    vertex_index = vertex_to_dof_map[i]
    print 'vertex index', vertex_index
    print 'at point', Vertex(mesh, vertex_index).point().str()
    print 'at coordinates', coords[vertex_index]
    print 'dof', x[i]

Comments

Yeah, tabulate_coordinates or tabulate_all_coordinates is the right function.

---

Thank you Jan!

---

It is not working with my Dolfin 1.4 on ipython notebook.

File "", line 1, in
runfile('/home/parallels/vertex_values_example.py', wdir='/home/parallels')

File "/usr/lib/python2.7/dist-packages/spyderlib/widgets/externalshell/sitecustomize.py", line 540, in runfile
execfile(filename, namespace)

File "/home/parallels/vertex_values_example.py", line 6, in
vertex_to_dof_map = V.dofmap().vertex_to_dof_map(mesh)

AttributeError: 'GenericDofMap' object has no attribute 'vertex_to_dof_map'

---

Use free functions dof_to_vertex_map or vertex_to_dof_map. Naming has changed, see issue 111.

---

Thank you Jan! I am using Dolfin 1.4 on ipython
of SageMath Cloud. Plot in line with dolfin
is not available at this time so I will try matplotlib.