periodic boundary conditions + numpy array

Question

I am trying to solve a problem in 1D with periodic boundary conditions. I know the exact solution of the problem. And Fenics solves it correctly.

But when I try to extract the solution values to a numpy array, it seems to give weird values which are oscillating. Plotting in viper gives perfect results.

from dolfin import *
import pylab as plt
import numpy as np

def exact(x, k):
    return np.sin(2*np.pi*x)/(k**2+4*np.pi**2)

class Source(Expression):
    def eval(self, values, x):
        values[0] = sin(2.0*DOLFIN_PI*(x[0]))

class PeriodicBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return bool(x[0] < DOLFIN_EPS and
                x[0] > -DOLFIN_EPS and on_boundary)
    def map(self, x, y):
        y[0] = x[0] - 1.0

n = 32
k = 2 * DOLFIN_PI
ksqrd=k**2

pbc = PeriodicBoundary()
mesh = UnitIntervalMesh(n)
V = FunctionSpace(mesh, "CG", 1, constrained_domain=pbc)

u = TrialFunction(V)
v = TestFunction(V)
f = Source()
a = inner(nabla_grad(u), nabla_grad(v))*dx + ksqrd * u*v*dx
L = f*v*dx

u = Function(V)
solve(a == L, u)

plot(u, interactive=True)

#--------------
x_array = mesh.coordinates()[:,0]
x_array = x_array[0:len(x_array)-1]

u_array = u.vector().array()

print u_array

x=np.linspace(0,1,4*n)
u_exact = exact(x, k)

plt.figure(figsize=(10,6))
ax=plt.subplot(111)
ax.plot(x_array, u_array, 'ro-',)
ax.plot(x, u_exact, 'g',)
plt.xlabel(r'$x$')
plt.ylabel(r'$u$')
plt.legend((r'$u_c(x)$', r'$u_E(x)$'), loc=0)
plt.show()

Plot from viper which gives the correct answer

and plot in matplotlib extracted as above

And the data values of u_array

[ -6.05208825e-18  -2.46687997e-03   2.46687997e-03  -4.83895913e-03
   4.83895913e-03  -7.02507980e-03   7.02507980e-03  -8.94123060e-03
   8.94123060e-03  -1.05137749e-02   1.05137749e-02  -1.16822808e-02
   1.16822808e-02  -1.24018431e-02   1.24018431e-02  -1.26448096e-02
   1.26448096e-02  -1.24018431e-02   1.24018431e-02  -1.16822808e-02
   1.16822808e-02  -1.05137749e-02   1.05137749e-02  -8.94123060e-03
   8.94123060e-03  -7.02507980e-03   7.02507980e-03  -4.83895913e-03
   4.83895913e-03  -2.46687997e-03   2.46687997e-03   2.80989811e-18]

You can see that every alternate element is negative of the previous one.

What is happening? I donot understand. Have I done some mistake in extracting the data values?

Any help is much appreciated..

Answer 1

Since FEniCS version 1.1 ordering of function dofs may not be related to ordering of respective mesh entities (here vertices). You can instruct DOLFIN not to reorder dofs

parameters['reorder_dofs_serial'] = False

so that dofs will be ordered according to mesh vertices in case of CG1 elements. But this is will have no effect in parallel. Generaly you can use helper member functions of V.dofmap() to access correct dofs. In your case using

u_array[V.dofmap().dof_to_vertex_map(mesh)]

should do the trick. See also http://fenicsproject.org/qa/121/python-how-to-set-expansion-coefficients-for-a-function?show=170#a170.

Comments

Thanks. This works.

Is there something else to be done when running in parallel?

---

Is there something else to be done when running in parallel?

This is not well-formed question. For example your plotting code is not of much utility in parallel, but you should be able to access process-local dofs

u_array[V.dofmap().dof_to_vertex_map(mesh)]

ordered by local vertex numbering.

---

Thanks again. I did not mean about plotting when running in parallel. I would need to get the data into numpy arrays and process them.

---

I think I have answered this - you can access local dofs this way.

---

Yes you have... Thanks