There has been some discussion on bitbucket about how to best support
evaluation and interpolation in parallel.
At present, user code like this:
mesh = UnitSquareMesh(10, 10)
Q = FunctionSpace(mesh, "CG", 1)
F = Function(Q)
F.interpolate(Expression("x[0]"))
print F(0.2, 0.2)
crashes in parallel. e.g. with mpirun -n 2
*** Error: Unable to evaluate function at point.
*** Reason: The point is not inside the domain. Consider setting
"allow_extrapolation" to allow extrapolation.
*** Where: This error was encountered inside Function.cpp.
Clearly "allow_extrapolation" is not needed, and this is actually
confusing to users.
Also, something like:
mesh = refine(mesh)
Q2 = FunctionSpace(mesh, "CG", 1)
F2 = Function(Q2)
F2.interpolate(F)
will crash in the same way, because the Mesh gets redistributed during
refine.
It should be possible to write the code for Function::eval to send the
results to the processes that need it. Clearly, however
it cannot always be a collective call, as local evaluation is an equally
common task.
One idea is to have two separate "Function::eval" methods, one
collective, and one local, or possibly differentiate using an MPI_Comm.
Alternatively, perhaps "Function::eval" should just return
std::pair<bool, Array<double> > to indicate whether evaluation is valid
on this process,
and delegate communication/extrapolation elsewhere.
It would be good to have get a wide range of ideas on how to fix this
problem in the best possible way, hopefully enhancing user code, and not
breaking it(!)
Thanks,
Chris R.
_______________________________________________
fenics mailing list
[email protected]
http://fenicsproject.org/mailman/listinfo/fenics
