Hi all,
I would like to run a lot of things in parallel. An example is a very simple code is here:
from fenics import *
mesh = RectangleMesh(Point(0, 0), Point(1, 1), 10, 10)
F = FunctionSpace(mesh, ‘P’, 1)
c = Function(F)
for i in range(F.dim()): c.vector()[i] = 1.
This is the comment I am running to do this: mpiexec -n 4 test.py
But I am getting the following error:
for i in range(F.dim()): c.vector()[i] = 1.
IndexError: Vector index out of range
How can I make it run?
Thanks
Tunc
F.dim() is the global dimension of the function space. Consider the following:
from dolfin import *
import numpy as np
mesh = RectangleMesh(Point(0, 0), Point(1, 1), 2, 2)
F = FunctionSpace(mesh, "CG", 1)
c = Function(F)
dm = F.dofmap()
local_range = dm.ownership_range()
local_dim = local_range[1] - local_range[0]
info("This process has ownership of DoFs: %s" % str(local_range))
c.vector().set_local(np.array([1.0]*local_dim, dtype=np.double))
Hi nate,
Thanks, it works.
By the way, in that way set_local() needs an array to assign values. But I would like to do that based on a criteria like this:
local_range = F.dofmap().ownership_range()
local_dim = local_range[1] - local_range[0]
for i in range(local_dim):
if c.vector().get_local()[i] <= 0.: c.vector().set_local()[i] = 1.
else : c.vector().set_local()[i] = 1.
That is not working since set_local cant accept indexed assignment. Is there a way to do that?
Thanks
Maybe try going through the PETSc backend:
from fenics import *
mesh = RectangleMesh(Point(0, 0), Point(1, 1), 10, 10)
F = FunctionSpace(mesh, 'P', 1)
c = Function(F)
cv = as_backend_type(c.vector()).vec()
Istart, Iend = cv.getOwnershipRange()
for i in range(Istart,Iend):
cv[i] = 1.0
cv.assemble()
Alternatively, the vector() function can pass arrays in a vectorized format. You can also do this if you’re just trying to assign all local DoFs to a Function:
from fenics import *
mesh = RectangleMesh(Point(0, 0), Point(1, 1), 10, 10)
F = FunctionSpace(mesh, ‘P’, 1)
c = Function(F)
c.vector()[:] = 1.0
print(c.vector().get_local(), MPI.rank(mpi_comm_world() ) )
Thanks all, aim is to generate a circle on the mesh, to do this, I first define a sign distance function (u_0) and try to make the nodal values either 1 or 0 based on the sign of the initial function. I found this way looks works (but not efficient in terms of time consumption and I see a strange function at the end):
u_0 = Expression(‘sqrt((x[0]-8.)(x[0]-8.) + (x[1]-8.)(x[1]-8.))-r’, degree=1, r=0.5)
c.assign(interpolate(u_0,V))
local_range = V.dofmap().ownership_range()
local_dim = local_range[1] - local_range[0]
print(local_dim)
arr = c.vector().get_local()
assign_array = c.vector().get_local()
for i in range(local_dim):
if arr[i] <= 0.:
assign_array[i] = 1.
c.vector().set_local(assign_array)
else:
assign_array[i] = 0.
c.vector().set_local(assign_array)
And what I see is a strange function as the end:
I guess that line is the boundary of where function divided to 2 CPUs. I could not understand the reason of that line. What I was expecting to see something like this (1 CPU case):
Here is the full code for testing:
from fenics import *
import numpy as np
cfile = File(“c.pvd”)
nx = ny = 80
mesh = RectangleMesh(Point(0, 0), Point(10, 10), nx, ny)
V = FunctionSpace(mesh, ‘P’, 1)
c = Function(V)
u_n = Function(V)
u_0 = Expression(‘sqrt((x[0]-8.)(x[0]-8.) + (x[1]-8.)(x[1]-8.))-r’, degree=1, r=0.5)
c.assign(interpolate(u_0,V))
local_range = V.dofmap().ownership_range()
local_dim = local_range[1] - local_range[0]
arr = c.vector().get_local()
assign_array = c.vector().get_local()
for i in range(local_dim):
if arr[i] <= 0.:
assign_array[i] = 1.
c.vector().set_local(assign_array)
else:
assign_array[i] = 0.
c.vector().set_local(assign_array)
cfile << c
Do you have any idea what is going on?
Thanks
I actually have a piece of code that does this. You can do this by defining a Expression. Example below is done a UnitSquareMesh:
import dolfin as dfn
class C0(dfn.Expression):
def eval(self, values, x):
x0 = 0.5 # define center x-coord
y0 = 0.5 # define center y-coord
# define as an ellipse with a = b
a = 0.25
b = a
if dfn.sqrt( ((x[0]- x0)/a)**2 +( (x[1]-y0)/b )**2 ) <= 1.0 + dfn.DOLFIN_EPS:
values[0] = 1.0
else:
values[0] = 0.0nx = ny = 80
mesh = dfn.RectangleMesh(dfn.Point(0, 0), dfn.Point(1, 1), nx, ny)
Q = dfn.FunctionSpace(mesh,“CG”,1)test_exp = C0(degree=1)
test = dfn.interpolate(test_exp, Q)
dfn.File(‘test.pvd’) << test
Here is the output from Paraview:
For your code above, I believe you just need to use apply() after you set_local:
e.g.
c.vector().set_local(assign_array)
c.vector().apply(“insert”)
Yes, apply worked, thank you. Thanks for the code as well. My case was a demo for asking about the problem. In general, I read some images and need to make assignment like 1 or 0 based on the image which at the end results in an arbitrary shape rather than a circle.
Apply() seems to be the solution.
Thank you
