we tested our program using the old and latest version of FEniCS (dolfin-version: 2019.1.0), and seems we found a bug in the vector assignment. For the demo purpose, we reproduced the error using the following codes:
from dolfin import *
parameters["reorder_dofs_serial"] = False
mesh = UnitCubeMesh(2,2,2)
V = FunctionSpace(mesh,'CG',1)
test_expression = Expression('x[0]*x[0]+x[1]*x[1]+x[2]*x[2]',degree=2)
test = interpolate(test_expression,V)
vertex_domains = MeshFunction('size_t',mesh, 0)
vertex_domains.array()[3:6]=1
test.vector()[vertex_domains.array()>0] = 100000
This code segment is supposed to change the 4th, the 5th, and the 6th values of the test vector to 100000. However, this program only produced a result with the 1st and 2nd values being changed. Can anyone explain this or this indeed is a bug? Thanks.
The function space DoF indices and mesh vertex indices are not guaranteed to be the same, even in older versions of DOLFIN. It seems you just got lucky… Consider accessing the information you need from V.dofmap(). Cf. here.
we have used the original order by setting parameters[“reorder_dofs_serial”] to False. Even if the DOF orders are not exactly the same, we are supposed to have three values changed. It is no sense only two were changed.
While I would recommend following Nate’s advice, I did figure out what is going on here. Basically, you are attempting to assign values using the NumPy convention of selecting indices with an array of Booleans. However, when doing this with a DOLFIN vector, the Booleans are being cast to integers, so, because there are some Trues and some Falses, which are cast to 1s and 0s respectively, the value of 100000 is being assigned to the 0-th and 1-st entries of the vector. The following seems to work as intended:
#test.vector()[vertex_domains.array()>0] = 100000
for i in range(0,len(vertex_domains.array())):
if(vertex_domains.array()[i] > 0):
test.vector()[i] = 100000