How to print nodal values of the numerical solution to a Neumann boundary-value problem?

Question

I would like to solve the Neumann boundary-value problem for the Laplace operator (using the methodology from http://fenicsproject.org/documentation/dolfin/dev/python/demo/documented/neumann-poisson/python/documentation.html). In addition, I would like to print the nodal values of the numerically obtained solution $u$.

At the end of this post, you can find my implementation of the methodology from http://fenicsproject.org/documentation/dolfin/dev/python/demo/documented/neumann-poisson/python/documentation.html. The script runs without any problems when the block of code between # BEGIN: PROBLEMATIC BLOCK OF CODE and # END: PROBLEMATIC BLOCK OF CODE is commented out but, when the same block of code is present, I get the following error message:

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at

---

*** fenics@fenicsproject.org

---

*** Remember to include the error message listed below and, if possible,
*** include a minimal running example to reproduce the error.

---

*** -------------------------------------------------------------------------
*** Error: Unable to access vector of degrees of freedom.
*** Reason: Cannot access a non-const vector from a subfunction.
*** Where: This error was encountered inside Function.cpp.
*** Process: unknown

---

*** DOLFIN version: 1.5.0
*** Git changeset: unknown
*** -------------------------------------------------------------------------

My implementation of the methodology from http://fenicsproject.org/documentation/dolfin/dev/python/demo/documented/neumann-poisson/python/documentation.html is as follows:

from dolfin import *

nx = 64
ny = 64
mesh = UnitSquareMesh(nx, ny)

V = FunctionSpace(mesh, "CG", 1) 
R = FunctionSpace(mesh, "R", 0)
W = V * R
(u, c) = TrialFunction(W)
(v, d) = TestFunction(W)

f = Expression("10*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02)")
g = Expression("-sin(5*x[0])")
a = ( inner(grad(u),grad(v)) + c*v + d*u )*dx
L = f*v*dx + g*v*ds

w = Function(W)
n_prob_laplace_operator = LinearVariationalProblem( a, L, w )
solver_np_laplace_operator = LinearVariationalSolver(n_prob_laplace_operator)
solver_np_laplace_operator.solve()
(u, c) = w.split()


# BEGIN: PROBLEMATIC BLOCK OF CODE
# Dump solution to the screen
u_nodal_values = u.vector()
u_array = u_nodal_values.array()
coor = mesh.coordinates()
if coor.shape[0] == u_array.shape[0]:  # degree 1 element
    for i in range(len(u_array)):
        print 'u(%8g,%8g) = %g' % (coor[i][0], coor[i][1], u_array[i])
else:
    print """\
 Cannot print vertex coordinates and corresponding values
 because the element is not a first-order Lagrange element.
 """

# END: PROBLEMATIC BLOCK OF CODE

plot(u, interactive=True)

[Remark: The numerical solution of the Neumann boundary-value problem for the Laplacian operatorion is understood in the sense discussed on http://fenicsproject.org/qa/7090/fenics-demo-demo_neumann-poisson-py ]

Comments

Hi, please note that the order of mesh vertices in coor array in general does not correspond to the order of coefficients in u_array. See e.g. here.

---

@MiroK: Thank you for this remark.

Answer 1

Hi kdma-at-dtu,

The Fenics error says

*** Error: Unable to access vector of degrees of freedom.
*** Reason: Cannot access a non-const vector from a subfunction.

This means that since u was defined by splitting the variable w, then you can not call the method vector() on u.
I believe this was a choice of the Fenics developers to prevent users to "shoot their own foot".

Try to substitute line

u_nodal_values = u.vector()

with

u_P1 = project(u, V)
u_nodal_values = u_P1.vector()

Best,

Umbe

Comments

Thanks for the suggestion, Umbe - that worked. I'm accepting your answer.

In addition, I am wondering about the following three things:

1. Since the program computes and displays a numerical solution $u$, what mesh and with respect to what finite-element basis is this $u$ computed?
2. With regard to the proper interpretation of the plot, it seems to me that the code below displays the origin at the lower-left corner with the positive x-axis pointing to the right and the positive y-axis pointing upwards; could you confirm if this interpretation is correct?
3. Are you able to indicate how I can calculate the integral of the numerically obtained solution to the Neumann boundary-value problem along a curve (say a circle contained within the unit square)? Would I then use $u$ or $u_{P1}$?

The code is:

from dolfin import *

nx = 4
ny = 4
mesh = UnitSquareMesh(nx, ny)

V = FunctionSpace(mesh, "CG", 1) 
R = FunctionSpace(mesh, "R", 0)
W = V * R
(u, c) = TrialFunction(W)
(v, d) = TestFunction(W)

f = Expression("10*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02)")
g = Expression("-sin(5*x[0])")
a = ( inner(grad(u),grad(v)) + c*v + d*u )*dx
L = f*v*dx + g*v*ds

w = Function(W)
n_prob_laplace_operator = LinearVariationalProblem( a, L, w )
solver_np_laplace_operator = LinearVariationalSolver(n_prob_laplace_operator)
solver_np_laplace_operator.solve()
(u, c) = w.split()


# Dump solution to the screen
u_P1 = project(u, V)
u_nodal_values = u_P1.vector()
u_array = u_nodal_values.array()
coor = mesh.coordinates()
if coor.shape[0] == u_array.shape[0]:  # degree 1 element
    for i in range(len(u_array)):
        print 'u(%8g,%8g) = %g' % (coor[i][0], coor[i][1], u_array[i])
else:
    print """\
Cannot print vertex coordinates and corresponding values
because the element is not a first-order Lagrange element.
 """

h = plot(u)
h.elevate(65)
h = plot(u, interactive=True, scale = 0.0)

---

Hi,

1) Since the program computes and displays a numerical solution u, what mesh and with respect to what finite-element basis is this u computed?

mesh = UnitSquareMesh(nx, ny)
is a structured mesh with triangular elements. You can do plot(mesh) to see the actual mesh.

V = FunctionSpace(mesh, "CG", 1)
Is the standard H^1-conforming space of piecewise linear finite elements (P1). You can increase the polynomial order if you wish (last parameter).

2) With regard to the proper interpretation of the plot, it seems to me that the code below displays the origin at the lower-left corner with the positive x-axis pointing to the right and the positive y-axis pointing upwards; could you confirm if this interpretation is correct?

Yes, this is the default view. Of course you can rotate the image...
If you are familiar with paraview, you can export the solution in paraview format
File("sol.pvd") << u

3) Are you able to indicate how I can calculate the integral of the numerically obtained solution to the Neumann boundary-value problem along a curve (say a circle contained within the unit square)? Would I then use u or uP1?

This is rather tricky. My suggestion is to use an external mesh generator to create a mesh that has the surface you want to integrate as internal boundary. Then you can use FEniCs to integrate the solution u over that internal boundary.
Paraview may alos be able to compute the integral for you.

---

@umberto: Thank you for the feedback.

I would like to ask two follow-up questions, if I may.

About 1): I am wondering what type of object u becomes after the line

(u, c) = w.split()

and how it is different from a finite-element function over V whose degrees of freedom are given by u_nodal_values (u_nodal_values is computed according to your suggestion above).

About 3): Can you suggest some external mesh generators which I may use for this purpose? (In particular, I would like to be able to specify the internal boundaries analytically.)

---

Hi,

About 1): I am not entirely sure what objects u and c are. But you can think at those as another Function objects that views (without owing) part of the data in w.
Since u and c are partial views of another object, it make sense some methods that may change the structure of u (and therefore of w) are disabled.

About 3): Fenics/dolfin has a tool (dolfin-convert) that can convert a mesh generated with commonly-used mesh generator in dolfin format.
Just type dolfin-convert --help for a list of supported mesh generators.

I never built my own meshes, so for more details you maybe be better off asking a new question :)

---

@umberto: Thank you for the suggestions.