Assign a array as function in the correct node

Hi community,

I have an issue with implementation of a array as a function, the problem is about the assign of the nodal value with the array. The value calculated is correctly, but each value is assigned in wrong node. My code is below,

from fenics import *
from mshr import *
import numpy as np
import matplotlib.pyplot as plt

mesh = Mesh('meshp.xml')


xs = SpatialCoordinate(mesh)
V = FunctionSpace(mesh, "CG", 1)

#print (V.dofmap().entity_dofs(mesh, 0))
fkappa=Function(V)

mesh_points=mesh.coordinates()
#print(mesh_points)

fx_dofs = V.dofmap().dofs()

mn=len(mesh_points)
## calculate nodal curvature

dt=1.0
teta = np.zeros(mn)
#Calculate first derive
dx=np.zeros(mn)
dy=np.zeros(mn)

for nod in range(mn):
    i=nod    
    teta[i]=dt*i
    if i >0 and i < mn-1:
        dx[i]=(mesh_points[i+1,0]-mesh_points[i-1,0])/(2*dt)
        dy[i]=(mesh_points[i+1,1]-mesh_points[i-1,1])/(2*dt)
    elif i == 0:
        dx[i]=(mesh_points[i+1,0]-mesh_points[-1,0])/(2*dt)
        dy[i]=(mesh_points[i+1,1]-mesh_points[-1,1])/(2*dt)
    else:
        dx[i]=(mesh_points[0,0]-mesh_points[i-1,0])/(2*dt)
        dy[i]=(mesh_points[0,1]-mesh_points[i-1,1])/(2*dt)

#Calculate second derive
ddx=np.zeros(mn)
ddy=np.zeros(mn)

for nod in range(mn):
    i=nod
    if i >0 and i < mn-1:
        ddx[i]=(mesh_points[i+1,0]-2*mesh_points[i,0]+mesh_points[i-1,0])/(dt**2)
        ddy[i]=(mesh_points[i+1,1]-2*mesh_points[i,1]+mesh_points[i-1,1])/(dt**2)
    elif i == 0:
        ddx[i]=(mesh_points[i+1,0]-2*mesh_points[i,0]+mesh_points[-1,0])/(dt**2)
        ddy[i]=(mesh_points[i+1,1]-2*mesh_points[i,1]+mesh_points[-1,1])/(dt**2)
    else:
        ddx[i]=(mesh_points[0,0]-2*mesh_points[i,0]+mesh_points[i-1,0])/(dt**2)
        ddy[i]=(mesh_points[0,1]-2*mesh_points[i,1]+mesh_points[i-1,1])/(dt**2)

kappa=np.zeros(mn)

for nod in range(mn):
    i=nod
    kappa[i]=(dx[i]*ddy[i]-dy[i]*ddx[i])/(dx[i]**2+dy[i]**2)**(3/2)

fkappa.vector()[fx_dofs] = kappa
kfunc =fkappa.vector().get_local()

vtkfile = File('curvature.pvd')
vtkfile << fkappa

The value calculated is the curvature and in paraview it looks as follows

and this isn’t correct

thanks in advance for any help

See for instance: Mapping from mesh.coordinates() to u.vector() - #2 by dokken