Assign a array as function in the correct node

Errors
1

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

imagen
imagen783×680 15.6 KB

and this isn’t correct

thanks in advance for any help

2

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