Discontinous Galerkin with periodic boundary conditions

[Nils Schween]

Dear deal.ii community,

first of all, a big thank you for all the effort going into this
software.

I am currently developing an application and I stumbled upon
implementing periodic boundary conditions for discontinuous Lagrange
elements.

I found the following thread
https://groups.google.com/g/dealii/c/WlOiww5UVxc/m/mtQJDUwiBQAJ

In this thread Dr. Arndt explains that when discontinuous elements were
used, it made more sense to enforce boundary conditions weakly instead
of using strong boundary conditions. In deal.ii this seems to be
tantamount to telling the MeshWorker to treat the boundary faces as
internal faces without using an additional constraint matrix enforcing
the equality of DoFs on the periodic parts of the domain.

If I understand correctly the MeshWorker can be informed to treat the
boundary faces as internal faces by modifying the triangulation with a
call of the function triangulation.add_periodocity(periodicity_vector).
Where the periodicity vector is filled by calling the function
collect_periodic_faces.

I tried this and it works as expected in 1D, but in 2D I get the
following error message.

> An error occurred in line <1001> of file </home/schween/.tmp/dealii/src/dealii-9.4.1/include/deal.II/grid/tria_iterator.h> in function
> const Accessor& dealii::TriaRawIterator<Accessor>::operator*() const [with Accessor = dealii::CellAccessor<2, 2>]
> The violated condition was:
> Accessor::structure_dimension != Accessor::dimension || state() == IteratorState::valid
> Additional information:
> You tried to dereference a cell iterator for which this is not
> possible. More information on this iterator: level=-1, index=-1,
> state=past_the_end

I adapted the step 12 tutorial such that it is easy for you to reproduce
the error. I attach the code to this email. Just change the template
argument for the parameter dim from 1 to 2 to see for yourself.

I should mention the adapted tutorial is not working with a distributed
triangulation. In my application I am. So at first I thought, that I was
using a wrong filter for the iterator range which I am handing over to
mesh_loop function. But since the error persist even when using only one
processor, it must be something else.

What am I doing wrong? Any help is appreciated.
Thank you very much.

Regards,
Nils

[Nils Schween]

Hmmm ... I really I do not would like to be impatient.

Can no one help me with setting up periodic boundary conditions for
discontinuous elements? Is it uncommon to do this?

Probably a small example code would be enough to get me started.

Thank you very much.
With best regards,
Nils

> --
> 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 dealii+un...@googlegroups.com.
> To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/87pmadlmn0.fsf%40mpi-hd.mpg.de.
>
> [2. adapted_step12 --- application/octet-stream; periodic_boundary_step12.cpp]...
>
>
> --
> Nils Schween
> PhD Student
>
> Phone: +49 6221 516 557
> Mail: nils.s...@mpi-hd.mpg.de
> PGP-Key: 4DD3DCC0532EE96DB0C1F8B5368DBFA14CB81849
>
> Max Planck Institute for Nuclear Physics
> Astrophysical Plasma Theory (APT)
> Saupfercheckweg 1, D-69117 Heidelberg
> https://www.mpi-hd.mpg.de/mpi/en/research/scientific-divisions-and-groups/independent-research-groups/apt


--
Nils Schween
PhD Student

Phone: +49 6221 516 557
Mail: nils.s...@mpi-hd.mpg.de
PGP-Key: 4DD3DCC0532EE96DB0C1F8B5368DBFA14CB81849

Max Planck Institute for Nuclear Physics
Astrophysical Plasma Theory (APT)
Saupfercheckweg 1, D-69117 Heidelberg
https://www.mpi-hd.mpg.de/mpi/en/research/scientific-divisions-and-groups/independent-research-groups/apt

[Daniel Arndt]

Nils,

Using DG with periodic boundary conditions is not very common to my knowledge. It seems you are hitting a bug in deal.II in Debug mode and it should work with

diff --git a/include/deal.II/meshworker/mesh_loop.h b/include/deal.II/meshworker/mesh_loop.h
index f46c737c2a..699264551e 100644
--- a/include/deal.II/meshworker/mesh_loop.h
+++ b/include/deal.II/meshworker/mesh_loop.h
@@ -570,7 +570,10 @@ namespace MeshWorker
 
                     // Now neighbor is on the same refinement level.
                     // Double check.
-                    Assert(!cell->neighbor_is_coarser(face_no),
+                    Assert((!periodic_neighbor &&
+                            !cell->neighbor_is_coarser(face_no)) ||
+                             (periodic_neighbor &&
+                              !cell->periodic_neighbor_is_coarser(face_no)),
                            ExcInternalError());
 
                     // If we own both cells only do faces from one side (unless

Best,

Daniel

To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/871qmofpp9.fsf%40mpi-hd.mpg.de.

[Nils Schween]

Dear Daniel,

thank you very much! Now it works perfectly.

Best,
Nils