Projecting solutions between distributed meshes

[Daniel Rivlin]

Good morning everyone,

I am attempting to write a program which projects a solution vector from mesh A to mesh B. Both meshes are parallel distributed and have different geometries (the first is adaptively refined, the latter has uniform grid size). In addition, the solution uses Q1 elements on mesh A and is to be piecewise constant on mesh B.

I have attempted to use VectorTools::project() for this purpose.

-        -   I create a child of the Function class which takes in a solution vector and its corresponding mesh (both parallel distributed). This function allows the solution vector to be evaluated at any coordinate within the domain (and should also consider ghost elements).

-         -    I insert this as a parameter into the project() function. This is projected onto a non-ghosted vector, which is then assigned to a ghosted vector and is subsequently compressed.

The resulting vector correctly projects to some elements, however other elements appear to be set to 0 (incorrectly). In addition, the number of elements that are zero and their locations vary with the number of processors used. I suspect that it may be something to do with project() not being able to handle ghost elements but I’m not sure. It is very possible that the Function used to evaluate the input vector at any coordinate is incorrect, i have included it in a text file below. I have tried using both FEFieldFunction() and later find_all_active_cells_around_point() but to no avail (the latter appears to return empty for these zeroed elements, even after I perturb the coordinate by a significant amount when no cell is found or use a high tolerance).

I then attempted the following crude workaround:

-            serialize the vector and triangulation

-            project them to a new serial vector using project()

-            scatter them back onto the parallel mesh

Unfortunately, this also failed. I am having difficulty copying the parallel mesh onto a serial one, which appears far more involved than I originally thought.

On a sidenote I have attempted to use interpolate() as well and was still unsuccessful (I encountered the same issue with elements set to 0 at the same locations for a given number of cores used).

Any help would be very, very much appreciated

 

Best wishes,

Daniel Rivlin

 

[Jan Philipp Thiele]

Hello Daniel,

for interpolate I know that it is between two different finite element spaces on the same mesh,

e.g. from Q2 to Q1 elements, or in your case from Q1 to dG0.

What you will need to do is construct a Q1 DoFHandler on each mesh and use
`interpolate_to_different_mesh` to get a Q1 solution vector on mesh B from the one on mesh A.

In a second step you can then use `interpolate` to get a dG0 solution on mesh B from the
temporary solution vector.

I don't think there is a function to do this in one go as it would be difficult to programmatically
decide which DoFHandler to use as an in-between.

Also note that both meshes have to share the same coarse mesh for this to work,

which should not be a problem if one is a refinement of the other.

And if i recall correctly the difference between interpolate and project is that interpolate is done locally on the element level,
while project works with the global solution and I think requires a linear equation system solve. But other people know more about that specific part than I.

Best,

Philipp

[Wolfgang Bangerth]

Daniel,
in addition to Philip's answer that I think both summarizes how to do
this in serial and the difference between interpolation and projections,
here are a couple other considerations:

> I am attempting to write a program which projects a solution vector from
> mesh A to mesh B. Both meshes are parallel distributed and have
> different geometries (the first is adaptively refined, the latter has
> uniform grid size).

This sentence already captures the key issue you are up against: If you
have two parallel triangulations, even if they cover the same geometry,
they are generally partitioned in different ways among the processes.
What that implies that if a process with mesh B needs information about
a point x on one of its locally owned cells, it needs to ask A for that
information. But x may not be part of the locally owned cells of mesh A,
and so the information is not available on the current process.

In other words, you cannot do interpolation (or projection) from one
parallel mesh to another without (a lot of) communication. Machinery
such as FEFieldFunction or...

> -      - I create a child of the Function class which takes in a
> solution vector and its corresponding mesh (both parallel distributed).
> This function allows the solution vector to be evaluated at any
> coordinate within the domain (and should also consider ghost elements).

...the function you describe is likely not prepared to do this kind of
communication. As a consequence...

> The resulting vector correctly projects to some elements, however other

> *elements appear to be set to 0* (incorrectly).

...I suspect that the zeros you get have something to do with querying
the other mesh at points that are simply not locally owned, and about
which you therefore don't know anything. I don't know, without seeing
the implementation, why you get zeros when I would have assumed that the
right answer should be to get an exception that simply terminates the
program because there is nothing you can do about the situation. But you
probably get what I'm getting at.


> In addition, the number
> of elements that are zero and their *locations vary with the number of
> processors used*. I suspect that it may be something to do with

> project() not being able to handle ghost elements but I’m not sure. It
> is very possible that the Function used to evaluate the input vector at
> any coordinate is incorrect, i have included it in a text file below. I
> have tried using both FEFieldFunction() and later
> find_all_active_cells_around_point() but to no avail (the latter appears
> to return empty for these zeroed elements, even after I perturb the
> coordinate by a significant amount when no cell is found or use a high
> tolerance).

Right. If you're looking for the cells that surround a point, but there
is no locally owned (or ghost) cell that does so, it rightly returns an
empty list.

What you want to do does not have an easy solution. The
FERemovePointEvaluation framework can help you with it, however. You
probably want to read through step-87 and step-89 and see how that can
be helpful to you.

Best
W.