How to assgin values to nodes?

1

Hi, I have a very simple question. How to assign a certain value to certain nodes? Following is my minimum working code trying to assign values of mesh nodes in u to that in h:

from dolfin import *

mesh=UnitSquareMesh(5,5)
V=FunctionSpace(mesh,'CG',1)
u=Expression('pow(x[0],2)+x[1]',degree=1)
u=interpolate(u,V)
h=interpolate(Constant(0),V)

print(u(1,1))
print(h(1,1))

h(1,1)=u(1,1)
print(h(1,1))

However, the code leads to an error of

h(1,1)=u(1,1)
    ^
SyntaxError: cannot assign to function call

Then what’s the correct way to assign values to nodes? Thanks!

1 Like
2

Consider the following:

from dolfin import *
import numpy as np
mesh = UnitSquareMesh(10,10)
V = FunctionSpace(mesh, "CG", 1)
x = V.tabulate_dof_coordinates()
point = [0.2,0.9]
# Locate which row of 
index = np.where((x == point).all(axis=1))[0][0]
u = Function(V)
u.vector().vec().setValueLocal(index, 2)
print(u.vector().get_local())
3 Likes
3

@dokken Hi, dokken. I’m using the method you provided above to solve my PDE problems. This method works perfectly when the FunctionSpace is in first order. But when it comes to second order, the method seems can’t work properly. Following is a simple example showing this problem:

    from dolfin import *
    import numpy as np

    mesh=UnitSquareMesh(2,2)
    V=FunctionSpace(mesh,'CG',2)
    u=TrialFunction(V)
    v=TestFunction(V)

    #define u0 by dirictly interpolate and then solve
    u0=interpolate(Constant(2),V)
    a=u*v*dx+inner(grad(u),grad(v))*dx
    L=u0*v*dx+1*v*ds
    ux=Function(V)
    solve(a==L,ux)
    print(ux.vector().get_local())

    #define u0 by assign values to each nodes and then solve
    u02=Function(V)
    x = V.tabulate_dof_coordinates()
    for point in mesh.coordinates():
        index = np.where((x == point).all(axis=1))[0][0]
        u02.vector().vec().setValueLocal(index, 2)
    L=u02*v*dx+1*v*ds
    solve(a==L,ux)
    print(ux.vector().get_local())

Here I used two ways to define the u0 for the simple Poisson equation: 1.by directly interpolate a constant value; 2. by assign the value to each nodes of a Function(V). Then the u0 was used to solve the Poisson problem and the solutions for the two definition methods were compared. When the FunctionSpace’s degree is 1, the two solution were exactly the same:

[6.33280088 6.00653708 6.00653708 6.17090373 5.76205639 6.17090373
 6.00653708 6.00653708 6.33280088]
[6.33280088 6.00653708 6.00653708 6.17090373 5.76205639 6.17090373
 6.00653708 6.00653708 6.33280088]

However, when the FunctionSpace’s degree comes to 2, different solutions arose:

[6.32825367 6.08298703 6.08298703 5.95836006 6.14296311 6.14296311
 6.32701613 5.83834168 5.9586526  5.89823277 6.14355869 5.89823277
 6.32701613 5.9586526  6.14355869 6.08298703 5.89823277 6.14355869
 6.08298703 5.95836006 5.89823277 6.14355869 6.32825367 6.14296311
 6.14296311]
[4.3303064  4.08990554 4.08990554 3.95619877 4.14206933 4.14206933
 4.33877858 3.84526773 3.95805171 3.89892358 4.14583303 3.89892358
 4.33877858 3.95805171 4.14583303 4.08990554 3.89892358 4.14583303
 4.08990554 3.95619877 3.89892358 4.14583303 4.3303064  4.14206933
 4.14206933]

It seems the above assigning method cannot be directly used in second order space. What can I do to improve this method to make it also workable in higher order space? Many thanks.

4

How are you expecting this to work with a higher order function space?
You can easily see why the results differ by saving and visualizing u0 and u02:

XDMFFile("u0.xdmf").write_checkpoint(u0, "u", 0)
XDMFFile("u02.xdmf").write_checkpoint(u02, "u", 0)

When you interpolate a constant over the whole space, every degree of freedom is set to 2.
When you assign to nodes, you assign to the vertices of your mesh (this means that not every dof is assigned the value 2, as a second order space has degrees of freedom on the facets).

2 Likes
5

Hi @dokken, I have a similar problem, in which I need to assign the value 0 to all DOFs of a vector at a single point (say, the origin of my domain). I have tried to implement it similarly to your answer, but it doesn’t work when running in parallel. Basically, I’ve tried

for dof in origin_dofs: # origin_dofs is a np array that contains the DOF id at the origin
            z.vector().vec().setValueLocal(dof, 0.)
            z_old.vector().vec().setValueLocal(dof, 0.)
            z.vector().apply('insert')
            z_old.vector().apply('insert')

The program stops running when it gets to the apply() method.

Any help with this issue?

6

You need to supply a minimal code example that reproduces the error, as the code you sent is not enough for me to reproduce/debug the issue.

7

Here is a MWE. The program arrives at a deadlock and doesn’t finish.

from ufl import *    
from dolfin import * 
import numpy as np
from mpi4py import MPI as nMPI  # Parallel calculations

Lbox = 20.0     
Hbox = 15.0     
Wbox = 100000.0  

N_el_x, N_el_y, N_el_z = 100, 100, 1
mesh = BoxMesh(Point(0.0, 0.0, 0.0), Point(Lbox, Hbox, Wbox), N_el_x, N_el_y, N_el_z)

comm = nMPI.COMM_WORLD
rank = comm.Get_rank()

V_z     = VectorFunctionSpace(mesh, "CG", 1)
z       = Function(V_z, name="z")

# Get the DOF coordinates 
coords = V_z.tabulate_dof_coordinates()

dofmap = V_z.dofmap()
origin_dofs = np.where((coords == [0., 0., 0.]).all(axis=1))[0]

print("Here is the error!")
for dof in origin_dofs: 
        z.vector().vec().setValueLocal(dof, 0.)
        z.vector().apply('insert')

I ran it with mpirun -np 2 python3 MWE.py

8

simply change this to:

print("Here is the error!")
for dof in origin_dofs: 
        z.vector().vec().setValueLocal(dof, 0.)
z.vector().apply('insert')
1 Like