Results of 1st and 2nd node of IntervalMesh are swapped

When doing a modal analysis on Timoshenko beam (IntervalMesh & MixedElement) and extracting modes from SLEPc, the result of 1st and 2nd node gets swapped. It looks like this: Figure_1 When you swap back the results, then the graph is okay: Figure_2

The whole code is given as follows:

import dolfinx
from dolfinx.io import XDMFFile
import ufl
from mpi4py import MPI
import matplotlib.pyplot as plt
import numpy as np
from petsc4py import PETSc
from slepc4py import SLEPc

# beam properties
n_ele = 12
start_x = 0.0
end_x = 1.0
E = 1e6
G = 1e4
I = 6.66667E-9
kappa = 0.53066
A = 0.0002

# 1D mesh and its function space
mesh = dolfinx.IntervalMesh(MPI.COMM_WORLD, n_ele, [start_x, end_x])
P2 = ufl.FiniteElement("CG", "interval", 2)
P1 = ufl.FiniteElement("CG", "interval", 1)
ME = ufl.MixedElement([P2, P2])
V = dolfinx.FunctionSpace(mesh, ME)

# Generalized strain measures
def chi(Phi):
    return Phi.dx(0)
def gamma(w, Phi):
    return -(w.dx(0) - Phi)

# Constitutive laws for bending and shear
def M(Phi):
    return E * I * chi(Phi)
def Q(w, Phi):
    return kappa * G * A * gamma(w,Phi)

u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
(w, phi) = ufl.split(u)
(w_ , phi_) = ufl.split(v)

# Timoshenko
stiffness_form = (ufl.inner(M(phi), chi(phi_)) + ufl.inner(Q(w, phi), gamma(w_, phi_))) * ufl.dx
KM = dolfinx.fem.assemble_matrix(stiffness_form)
KM.assemble()

mass_form = A * ufl.dot(u, v) * ufl.dx
MM = dolfinx.fem.assemble_matrix(mass_form)
MM.assemble()

pc = PETSc.PC().create(MPI.COMM_WORLD)
pc.setType(pc.Type.LU)

ksp = PETSc.KSP().create(MPI.COMM_WORLD)
ksp.setType(ksp.Type.PREONLY)
ksp.setPC(pc)

st = SLEPc.ST().create(MPI.COMM_WORLD)
st.setType(st.Type.SINVERT) # SINVERT / SHIFT / CAYLEY / PRECOND
st.setKSP(ksp)
st.setShift(0)

eigensolver = SLEPc.EPS().create(MPI.COMM_WORLD)
#eigensolver.setST(st)
eigensolver.setOperators(KM, MM)
eigensolver.setType(eigensolver.Type.KRYLOVSCHUR)
eigensolver.setDimensions(3)
eigensolver.setTolerances(1e-3, 1000)
eigensolver.setWhichEigenpairs(eigensolver.Which.SMALLEST_REAL)
eigensolver.setProblemType(SLEPc.EPS.ProblemType.GHEP)
eigensolver.setFromOptions()
eigensolver.solve()

its = eigensolver.getIterationNumber()
sol_type = eigensolver.getType()
nev, ncv, mpd = eigensolver.getDimensions()
tol, maxit = eigensolver.getTolerances()
nconv = eigensolver.getConverged()
print("")
print("Number of iterations of the method: %i" % its)
print("Solution method: %s" % sol_type)
print("Number of requested eigenvalues: %i" % nev)
print("Stopping condition: tol=%.4g, maxit=%d" % (tol, maxit))
print("Number of converged eigenpairs: %d" % nconv)

vr, vi = KM.getVecs()
if nconv > 0:
    print("")
    print("        k          ||Ax-kx||/||kx|| ")
    print("----------------- ------------------")
    for i in range(nconv):
        k = eigensolver.getEigenpair(i, vr, vi)
        error = eigensolver.computeError(i)
        print(" %12f       %12g" % (np.sqrt(k.real) * 0.15915494327, error)) # rad/s to Hz
        if i < 3:
            er = PETSc.Vec.getArray(vr)
            eigenmode = dolfinx.Function(V)
            eigenmode.vector.array[:] = (er / np.max(np.abs(er)))
            eigenmode.vector.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
            result = eigenmode.sub(0).compute_point_values()
            #tmp = np.copy(result[0])
            #result[0] = result[1]
            #result[1] = tmp
            plt.plot(np.linspace(start_x, end_x, n_ele + 1), result)

    print("")
plt.show()

FEniCS / dolfinx

Results of 1st and 2nd node of IntervalMesh are swapped #1673