Problem with assigning new values to a function on specific vertices.


140
views
1
3 months ago by
Nima  
Hi,

So I have the following code:
from dolfin import *
import numpy as np

mesh = UnitSquareMesh(40,40)
coords = mesh.coordinates()

V = FunctionSpace(mesh, 'P', 1)

u = Expression('x[0]',degree=1)

U=project(u,V)


for v in vertices(mesh):
    if coords[v.index()][0]<=0.4 :
        U.vector()[v.index()] = 0

plot(U)​
and I want to set the value of the function U equal to zero on the points p=(x,y) for which x<=0.4.
I expect to see a plot with everything to the left of the line x=0.4 in blue (U=0). Instead, i get the following:

What am I doing wrong in here?

Thank you

Community: FEniCS Project

1 Answer


5
3 months ago by
You wrongly assume that the index of vertex is the same as index of degree of freedom of the function space.
Numbering of vertices do not coincide with numbering of degrees of freedom (for example for higher order function spaces there are degrees of freedom not "based" on any vertex).

Just add vertex to dof mapping, e.g.

from dolfin import *
import numpy as np
import matplotlib.pyplot as plt

mesh = UnitSquareMesh(40, 40)
coords = mesh.coordinates()

V = FunctionSpace(mesh, 'P', 1)

u = Expression('x[0]', degree=1)

U = project(u, V)

# Extract vertex to dof mapping
# It takes index of vertex and returns 
# degree of freedom index
vtdmap = vertex_to_dof_map(V)

for v in vertices(mesh):
    if coords[v.index()][0] <= 0.4:
        # Apply vtdmap
        U.vector()[vtdmap[v.index()]] = 0

plot(U)
plt.show()​
Thank you so  much. I was always wondering what those mapping do. Now I know thanks to you.
written 3 months ago by Nima  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »