On mesh refinement & solution transfer with Raviart-Thomas

[Charlie Za]

Hi there,

I've been playing with mesh refinement and solution transfer with the Raviart-Thomas field, but found a very strange thing: in areas where the mesh is coarsened, the solution becomes wrong.

In a very simple 2D test, I first created a uniform RT field, i.e. (1, 0), with the help of `VectorTools::interpolate`, and then tried to refine the left half of a square domain, while coarsening the rest. Of course, I also tried to transfer the solution onto the new mesh. This function `refine_mesh` is as follows.

```

  template <int dim>

  void

  TestProblem<dim>::refine_mesh()

  {

    unsigned int cell_no = 0;

    for (const auto &cell : dof_handler.active_cell_iterators())

    {

      Point<dim> center = cell->center();

      if (center[0]<0.5)

      {

        cell->set_refine_flag();

      }

      else

      {

        cell->set_coarsen_flag();

      }

    }

    SolutionTransfer<dim> solution_trans(dof_handler);

    Vector<double> previous_solution;

    previous_solution = solution;

    triangulation.prepare_coarsening_and_refinement();

    solution_trans.prepare_for_coarsening_and_refinement(previous_solution);

    triangulation.execute_coarsening_and_refinement();

    setup_dofs();

    solution_trans.interpolate(previous_solution, solution);

    constraints.distribute(solution);

  }

```

As you can see in the attached plot, however, the transferred solution in the coarsened area is wrong but OK in the refined part. I'm not sure what goes wrong here in my code and would very much appreciate it if you could give me a hint. A minimal but complete code to reproduce this plot is also attached.

Thanks.

Charlie

[Wolfgang Bangerth]

Charlie,
I don't have time in the next few days to look into the exact cause, but I
have a suspicion and if you are up for some digging, here are a few pointers.

When you transfer the solution from a cell to its children, then
SolutionTransfer what are called the "prolongation matrices" that each element
stores. If you want, you can output these for the element you have by calling
fe.get_prolongation_matrix(...): each of these matrices is used to compute the
dofs_per_cell unknowns on the child cell from the dofs_per_cell unknowns on
the parent cell. This seems to work correctly.

On the other hand, if you want to restrict from cells to their parent cell,
SolutionTransfer uses the "restriction matrices" that you can get via
fe.get_restriction_matrix(...). From your example, it seems to me that these
matrices are wrong by a factor of two. Why that is so I don't know for sure,
but you could take a look at the place where they are computed: In
source/fe/fe_raviart_thomas.cc in the function
FE_RaviartThomas<dim>::initialize_restriction(). There is one place where we
multiply by a factor
Utilities::fixed_power<dim - 1>(.5)
and it's conceivable that that is wrong. You could try to remove the factor
and see what is happening -- you already have an excellent test case to check
this!

It would be fantastic if you could check whether this actually works!

Best
Wolfgang

> visit0000.png
> Thanks.
>
> Charlie
>
> --
> The deal.II project is located at http://www.dealii.org/
> <https://nam10.safelinks.protection.outlook.com/?url=http%3A%2F%2Fwww.dealii.org%2F&data=04%7C01%7CWolfgang.Bangerth%40colostate.edu%7Ca29600b7daca469a9ace08d8f6eccf4a%7Cafb58802ff7a4bb1ab21367ff2ecfc8b%7C0%7C1%7C637530844248762883%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C2000&sdata=Mx8xDq9h1EZP9qrIZzaRkR%2BXCub4BOOyTUo9qODdNis%3D&reserved=0>
> For mailing list/forum options, see
> https://groups.google.com/d/forum/dealii?hl=en
> <https://nam10.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgroups.google.com%2Fd%2Fforum%2Fdealii%3Fhl%3Den&data=04%7C01%7CWolfgang.Bangerth%40colostate.edu%7Ca29600b7daca469a9ace08d8f6eccf4a%7Cafb58802ff7a4bb1ab21367ff2ecfc8b%7C0%7C1%7C637530844248772873%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C2000&sdata=n04HDAGJBUEqB7w6LgZwpQ3wO3N%2B6oZqifsp5BeEtxk%3D&reserved=0>
> ---
> You received this message because you are subscribed to the Google Groups
> "deal.II User Group" group.
> To unsubscribe from this group and stop receiving emails from it, send an
> email to dealii+un...@googlegroups.com
> <mailto:dealii+un...@googlegroups.com>.
> To view this discussion on the web visit
> https://groups.google.com/d/msgid/dealii/be6e661d-a3a9-4d74-b1a2-015d168ff9f1n%40googlegroups.com
> <https://nam10.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgroups.google.com%2Fd%2Fmsgid%2Fdealii%2Fbe6e661d-a3a9-4d74-b1a2-015d168ff9f1n%2540googlegroups.com%3Futm_medium%3Demail%26utm_source%3Dfooter&data=04%7C01%7CWolfgang.Bangerth%40colostate.edu%7Ca29600b7daca469a9ace08d8f6eccf4a%7Cafb58802ff7a4bb1ab21367ff2ecfc8b%7C0%7C1%7C637530844248772873%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C2000&sdata=HysIkpKppj8Y21tQDkv6YRWMbLMW%2FskOL9mccHPnBCc%3D&reserved=0>.


--
------------------------------------------------------------------------
Wolfgang Bangerth email: bang...@colostate.edu
www: http://www.math.colostate.edu/~bangerth/

[Charlie Za]

Hi Wolfgang,

Thanks very much for your explanation.

I followed your advice by removing the fixed_power of (.5) in `initialize_restriction`, line 329 at [fe_raviart_thomas.cc](https://www.dealii.org/current/doxygen/deal.II/fe__raviart__thomas_8cc_source.html ), and am happy to tell you that it solved the problem: as the attached plots show, solutions are correctly transferred in both 2D and 3D cases.

Regards,

Charlie

[Marc Fehling]

Hi Charlie!

It looks like you found the cause for the issue! The transferred solution now looks like one would expect.

It seems like this is a general problem with the restriction matrices in the Raviart-Thomas elements. Would you mind writing a patch for the deal.II library so that your solution and your testcase can be part of it? You can find a tutorial on how to do it here. We would be happy to assist you in preparing that patch :-)

Best,

Marc

[Charlie Za]

Hi Marc,

Thanks for your reminder. I'll have a good look at the tutorial and come back to you when I'm ready. As far as I can see now, this might involve changes in 3 places:

1) Remove fixed_power of (.5) in `initialize_restriction`, line 329 at [fe_raviart_thomas.cc](https://www.dealii.org/current/doxygen/deal.II/fe__raviart__thomas_8cc_source.html ): this is simple.

2) Update [dealii/tests/fe/rt_5.output](https://github.com/dealii/dealii/blob/981214b60b491f465dd8cf04e24989a2d806229c/tests/fe/rt_5.output ). [dealii/tests/fe/rt_5.cc](https://github.com/dealii/dealii/blob/981214b60b491f465dd8cf04e24989a2d806229c/tests/fe/rt_5.cc ) is designed to verify correctness of the restriction matrix of a Raviart-Thomas element, and thus we need to update rt_5.output.

3) As to "patch...your testcase can be part of it", I wonder where is a proper place to put this testcase? Several tests exist in `dealii/tests/hp` to test solution transfer with hp elements, but I didn't find any equivalent tests for normal/standard elements (Please forgive me if I'm wrong here).

Regards,

Charlie

[Marc Fehling]

Hi Charlie,

to 2) Yes, this update would be necessary. It should be okay to overwrite the .output file with the actual output of the test after you applied 1).

to 3) I now wonder to where to place it as well! I'm surprised that we have lots of tests for SolutionTransfer in the hp context or the parallel distributed context, but not many for single finite elements in serial applications. There is one in `tests/bits`, but that's just to verify the functionality. Since the SolutionTransfer class definition is located in the `numerics` folder, I would suggest to create a test `solution_transfer_01` or `solution_transfer_rt` in `tests/numerics`. But that's just a suggestion, feel free to use a more intuitive location :-)

Best,
Marc