Extract solution from nodes of certain coordinates

Question

Say I have computed a nodal finite element solution. I want to do some post-processing work with only certain nodes. How do I extract and/or print out both the solution at the nodes and their respective coordinates?

For example, I want to print the solutions of all nodes that have x coordinate == 0.5 (inside a unit square domain). I also want to print out the respective coordinates. Thanks

Answer 1

Hi, consider

from dolfin import *
import numpy as np

mesh = UnitSquareMesh(20, 20)
V = FunctionSpace(mesh, 'CG', 1)
u = interpolate(Expression('x[0]*x[1]'), V)

dofmap = V.dofmap()
dof_x = dofmap.tabulate_all_coordinates(mesh).reshape((V.dim(), -1))
x = dof_x[:, 0]
# Extract dof indices where some condition on the first coordinate of dof
# is met
indices = np.where(np.logical_and(x > 0.2, x < 0.4))[0]
# Get coordinates of dof
xs = dof_x[indices]
# Get value of dof 
vals = u.vector()[indices]

for x, v in zip(xs, vals):
    print x, v

Comments

For those interested in doing the same thing on a set of nodes defined as a facet mesh function, see MiroK's answer (thanks again to him) to this question: https://fenicsproject.org/qa/11876/extract-solution-at-a-set-of-nodes. Martin

---

For versions after 2016.1 see https://fenicsproject.org/qa/10782/tabulate_all_coordinates-removed-in-2016-1-0 and use dofs = V.tabulate_dof_coordinates().reshape(V.dim(),mesh.geometry().dim()) instead