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 On Sunday, April 4, 2021 at 11:55:57 PM UTC+2 [email protected] wrote: > 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 > <https://github.com/dealii/dealii/wiki/Contributing>. We would be happy > to assist you in preparing that patch :-) > > Best, > Marc > > On Sunday, April 4, 2021 at 2:47:58 PM UTC-6 Charlie Za wrote: > >> 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. >> >> [image: visit00_2d.png][image: visit00_3d.png] >> >> Regards, >> Charlie >> >> On Sunday, April 4, 2021 at 2:37:20 AM UTC+2 Wolfgang Bangerth wrote: >> >>> >>> 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 >>> >>> >>> On 4/3/21 4:06 PM, Charlie Za wrote: >>> > 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. >>> > 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 [email protected] >>> > <mailto:[email protected]>. >>> > 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: [email protected] >>> www: http://www.math.colostate.edu/~bangerth/ >>> >>> -- The deal.II project is located at http://www.dealii.org/ For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en --- 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 [email protected]. To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/bb969595-e030-4df3-b379-437978c710ean%40googlegroups.com.
