> Since not all elements are supported, I think this function should
> have a separate name for now, perhaps
>
> interpolate_lagrange()
>
> My suggestion would be that the documentation explains that this
> function works in parallel for non-matching meshes (in contrast to the
> regular interpolate) but does only support Lagrange spaces.
Sounds good to me.
Mikael
>
>> Hi,
>>
>> This is in reference to Issue 162
https://bitbucket.org/fenics-project/dolfin/
>> issue/162/interpolation-nonmatching-mesh-in-parallel
>>
>> Currently it is not possible in parallel to interpolate a
>> dolfin Function u0 defined on one mesh to a Function with the
>> same FunctionSpace defined on a different mesh. I would like to be able to do
>> this
>>
>> mesh0 = UnitSquareMesh(100, 100)
>> mesh1 = UnitSquareMesh(20, 20)
>> V0 = FunctionSpace(mesh0, "CG", 1)
>> V1 = FunctionSpace(mesh1, "CG", 1)
>> u0 = interpolate(Expression("x[0]*x[1]"), V0)
>> u0.update()
>> u = interpolate_nonmatching_mesh(u0, V1)
>>
>> or probably overloaded like this
>>
>> u = interpolate(u0, V1)
>>
>> Now there are probably many different ways this can be implemented. The most
>> obvious approach that comes to mind is simply to overload `Function::eval`
with
>> a parallel version. There are many reasons why this doesn’t quite seem to
work.
>> Consider the current `interpolate` implementation
>>
>> FunctionSpace::interpolate(GenericVector& expansion_coefficients,
>> const GenericFunction& v)
>> …
>> // Iterate over mesh and interpolate on each cell
>> UFCCell ufc_cell(*_mesh);
>> for (CellIterator cell(*_mesh); !cell.end(); ++cell)
>> {
>> // Update to current cell
>> ufc_cell.update(*cell);
>>
>> // Restrict function to cell
>> v.restrict(&cell_coefficients[0], *_element, *cell, ufc_cell);
>>
>> // Tabulate dofs
>> const std::vector<dolfin::la_index>& cell_dofs = _dofmap->cell_dofs(cell->
>> index());
>>
>> // Copy dofs to vector
>> expansion_coefficients.set(&cell_coefficients[0],
>> _dofmap->cell_dimension(cell->index()),
>> &cell_dofs[0]);
>> }
>>
>> The cell loop runs on all processors for the part of the mesh of V1 that
lives
>> on that process. `v.restrict` calls `restrict_as_ufc_function` since V0 and
V1
>> have different meshes, `restrict_as_ufc_function` calls
`element.evaluate_dofs`
>> that again calls `u0.eval` for each dof on the cell. Now the major problem
here
>> that makes it difficult to extend this to parallel is that the cell iteration
>> is not collective and `eval` is not called collectively. Since the
coordinates
>> of dofs on a cell in the V1-mesh could live on any other cpu for the V0-mesh,
>> the `eval` needs to be called collectively with some try/catch clause and
with
>> a lot of broadcasting (of coordinates and result) to obtain what is
>> required. Alternatively, the non-collective `eval` that fails must be
captured
>> and handled using MPI++ later.
>>
>> Faced with these problems my solution has instead been to create a
stand-alone
>> version for parallel interpolation on nonmatching meshes. Not trying to
>> overload `eval` or fit it into the current `interpolate` framework. The code
>> that performs the interpolation is attached to issue 162. The code is fairly
>> well documented and should, hopefully, be possible to follow. The algorithm
is
>> basically like this
>>
>> 1) On each process tabulate all coordinates for dofs in the functionspace
>> we’re interpolating to, i.e., V1 in the above code.
>>
>> x = V1.dofmap()->tabulate_all_coordinates(mesh1)
>>
>> 2) Create a map from dof number in V1 to component number in Mixed Space.
>>
>> This is necessary to be able to handle mixed spaces. I do not think I have
>> found the best solution to this problem, so any help here would be much
>> appreciated. I just realized that the current solution does not work
>> for special spaces like BDM or RT.
>>
>> 3) Compute `u0.eval` for all points x computed in 1).
>>
>> Problem here is that u0 and u have different meshes and as such a point in u
>> will not necessarily, or not even likely, be found on the same processor for
>> u0. Hence the dofs coordinates must be passed around and searched on all
ranks
>> until found. When found the result is returned to the process that owns the
>> dof. (Alternatively, but slower, `u0.eval` could simply be called
>> collectively.) A major issue here for performance is that all points not
found
>> on one process are collected in a vector and passed on to the
>> next process through just one single MPI::send_recv. All points not found are
>> sent and points not found by process with one lower rank are received.
>>
>> 4) Finally, set all values in u using the dof to component map.
>>
>>
>> The problem is inherently unscalable by nature since the efficiency depends
on
>> the distribution of two different meshes. If they happen to overlap much,
then
>> the routine may be faster with 2 than one cpu, but in the end that depends on
>> the mesh partitioning that you cannot really control, at least not just using
>> the default mesh-partitioner. Consider running a problem with 2 CPUs. Here it
>> should be quite obvious that each dof must be searched on average 1.5 times
>> (half are found on this CPU and the rest on the other). Since each process'
>> mesh is half the size of the original, the total cost is then approximately
>> 1.5*0.5 = 0.75. So, not counting any cost of MPI communication, you should
>> expect to reduce the computational time to 75% of the original. However, the
>> MPI communications are extensive (half the mesh is sent back and forth) and
you
>> should honestly not expect any speedup at all. The estimate for N cpus will
>> similarly be (N+1) / N / 2 + MPI communications.
>>
>> If you have any comments to the implementation or, even better, suggestions
on
>> how to improve it, then please let me know by following up on this mail. If
you
>> like the current solution and you’d like to see this implemented soon
>> in dolfin, then please vote for issue 162 on bitbucket:-)
>>
>> Meanwhile, the `interpolate_nonmatching_mesh` function is available on
https://
>> github.com/mikaem/fenicstools
>>
>> Best regards
>>
>> Mikael
>>
>>
>
>> _______________________________________________
>> fenics mailing list
>> [email protected]
>> http://fenicsproject.org/mailman/listinfo/fenics
>
> _______________________________________________
> fenics mailing list
> [email protected]
> http://fenicsproject.org/mailman/listinfo/fenics
