Implement array as source term in variational form

variational formulation
1

Hello Fenics community

I have a simple issue. I have to implement a term in my variational problem from an array, I don’t know if this could be implemented with interpolate or Project, i would be grateful if someone can help me with this. Below is the short code of my problem, i want to set a source term with the array kappa.

from fenics import *
from mshr import *
import numpy as np

mesh = Mesh('meshp.xml')

xs = SpatialCoordinate(mesh)
V = VectorFunctionSpace(mesh, "CG", 1)
#print (V.dofmap().entity_dofs(mesh, 0))

mesh_points=mesh.coordinates()

mn=len(mesh_points)
## calculate nodal curvature

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

for nod in range(mn):
    i=nod-1
    if i !=0 & i !=mn-1:
        dx[i]=(mesh_points[i+1,0]-mesh_points[i-1,0])/2
        dy[i]=(mesh_points[i+1,1]-mesh_points[i-1,1])/2
    elif i == 0:
        dx[i]=(mesh_points[i+1,0]-mesh_points[-1,0])/2
        dy[i]=(mesh_points[i+1,1]-mesh_points[-1,1])/2
    else:
        dx[i]=(mesh_points[0,0]-mesh_points[i-1,0])/2
        dy[i]=(mesh_points[0,1]-mesh_points[i-1,1])/2

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

for nod in range(mn):
    i=nod-1
    if i !=0 & i !=mn-1:
        ddx[i]=(mesh_points[i+1,0]-2*mesh_points[i,0]+mesh_points[i-1,0])
        ddy[i]=(mesh_points[i+1,1]-2*mesh_points[i,1]+mesh_points[i-1,1])
    elif i == 0:
        ddx[i]=(mesh_points[i+1,0]-2*mesh_points[i,0]+mesh_points[-1,0])
        ddy[i]=(mesh_points[i+1,1]-2*mesh_points[i,1]+mesh_points[-1,1])
    else:
        ddx[i]=(mesh_points[0,0]-2*mesh_points[i,0]+mesh_points[i-1,0])
        ddy[i]=(mesh_points[0,1]-2*mesh_points[i,1]+mesh_points[i-1,1])

kappa=np.zeros(mn)

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

print(kappa)

Thanks in advance

2

As it seems like all your computations are based on the mesh coordinates, you can use the dof_to_vertex_map to insert values into a function from a CG 1 space, see; Mapping from mesh.coordinates() to u.vector() - #2 by dokken

1 Like