Periodic boundary conditions

[Jan Skala]

Dear Deal.II community,

I'm new in the Deal.II and recently I'm studing boundary conditions from tutorial 45.  The text and example is very clear and helpfull,
but I have one question related to this:
Is it possible to construct periodic boundary conditions by adding element from one border to the element from opposite border as its neighbor?
This approach could remove necessity of using constraint matrix. as it is used in the tutorial 45.
I can imagine that can be problematic for some Deal.II methods (may be mesh refinement) is the elements have neighbors without sharing nods because they are spatially
separated, but still I want to ask if it is possible to have two spatially separated elements as neighbors?
Thank you in advance.
Best wishes,

Jan Skala

[Matthias Maier]

Hi,

On Mon, Aug 31, 2015, at 07:08 CDT, Jan Skala <jsk...@gmail.com> wrote:

> Dear Deal.II community,
>
> I'm new in the Deal.II and recently I'm studing boundary conditions from
> tutorial 45.

This examle step is outdated and will be removed in the future [1]. We
have builtin support for imposing periodic boundary conditions for a
long time now [2]

> The text and example is very clear and helpfull,

Unfortunately, the whole strategy only works for Q1 elements - the
library function is much more flexible.

> but I have one question related to this:
> Is it possible to construct periodic boundary conditions by adding element
> from one border to the element from opposite border as its neighbor?

> This approach could remove necessity of using constraint matrix. [...][

We plan to provide this feature as a side-effect of rewriting the
triangulation backend in deal.II. But this is a major task. So, it won't
be available in the near future.

Best,
Matthias

[1] https://github.com/dealii/dealii/pull/1421
[2] http://dealii.org/developer/doxygen/deal.II/DEALGlossary.html#GlossPeriodicConstraints

[Pavel Kus]

Dear Matthias,

thank you for your reply. The suggested approach is unfortunately not suitable
for us. Our goal is to implement periodic boundary condition in the frame of
discontinuous Galerkin method. We use implementation similar to tutorial
step 12. What we therefore need is to handle faces on the periodic boundary
in integrate_face_term() rather than in integrate_boundary_term(), with a
correct "neighbouring" element provided. Is there any possible workaround
for this problem? Or does it refrain us from the use of MeshWorker framework,
as it is used in step 12?

Thank you.

Best,

Pavel

Dne pondělí 31. srpna 2015 14:23:26 UTC+2 Matthias Maier napsal(a):

[Sungho Yoon]

I see, It's what I have done for periodicity in dg. It's not difficult to do. I have done with adding two steps manually, Deal.II does not provide with this currently

1. Add the non-zero entries on the periodic faces in "make_flux_sparsity_pattern ()"
  - You need to modify the source code inside deal.II, not difficult to modify manually. must be readable to do.

2. And you search for neighbor cell also by manully in intergrating the "assemble_face" part.
  - Yes, you still have to do manully, do not follow the step12. it was functioned in the framework in meshwork.... namespace.
  - would recommend to follow the step30. it is also written based on the DG in old style, but give more practical ways.

However, coupling with AMR for this is a bit challenged but I belive it should be reachable in the context of deal.II

Sungho

2015년 8월 31일 월요일 오전 9시 20분 46초 UTC-7, Pavel Kus 님의 말:

[Andreas Krämer]

Hi there,

I have a similar idea as Sungho proposed implemented in my code.

• Concretely, I use a std::map<active_cell_iterator, std::pair<cell_iterator, face_id> > global_map that stores cells on opposite boundaries (and the corresponding face)
• It is created in the following way:
• The opposite boundaries have two different boundary indicators
• Create two std::map<some_key, active_cell_iterator> (one for each of the two opposite boundaries
• For all cells: if the cell has one of the two boundary indicators, put it into the map
• The key is determined by the coordinates of the face centers (e.g. the y-coordinate for a periodic 2d boundary in x-direction)
• Iterate over the two maps to associate cells with similar keys to each other and store them into global_map
• I modified make_flux_sparsity_pattern to couple the DoFs of cells that are in the global_map
• In the assembly: use the global_map to determine the neighbor cell

It works well on shared memory, but extension to distributed memory is not straightforward.

If you want to have a look at the code, feel free to contact me.

@Sungho Yoon: Do you have an idea how to extend that to a distributed triangulation?

[Daniel Arndt]

Hi Andreas,

what you describe resembles in many parts the already implemented for applying periodic boundary conditions.

Basically, you can call  GridTools::collect_periodic_faces() to get all periodic faces on the coarsest level 

and then proceed to the active cells (and faces) similar to the implementation of DoFTools::make_periodicity_constraints().
To get the required information on a parallel::distributed::Triangulation you can use parallel::distributed::Triangulation::add_periodicity().

As Matthias already said, more information (for cG) can be found in [1] and the new step-45 [2].

If this is not sufficient: What are the problems you are facing with a distributed triangulation?

Bests,

Daniel

[1] http://dealii.org/developer/doxygen/deal.II/DEALGlossary.html#GlossPeriodicConstraints

[2] https://github.com/dealii/dealii/pull/1421 

[Andreas Krämer]

Hi Daniel,

thanks for your very helpful comment.

I am now looking at  DoFTools::make_periodicity_constraints(), trying to transfer it to DG.

Would it suit the deal.II concept to replace the constraint matrix by a map periodic_neighbor between opposite cells

in order to use fluxes rather than constraints?

In building the sparsity pattern and assembling the system, I could just add periodic_neigbor.at(cell) along the periodic boundary, according to the usual cell.get_neighbor() calls.

Best,

Andreas

[Daniel Arndt]

Hi Andreas,

Would it suit the deal.II concept to replace the constraint matrix by a map periodic_neighbor between opposite cells

in order to use fluxes rather than constraints?

In building the sparsity pattern and assembling the system, I could just add periodic_neigbor.at(cell) along the periodic boundary, according to the usual cell.get_neighbor() calls.

Certainly, you can create such a map for the active faces from the information provided by make_flux_sparsity_pattern similar to the recursive call in make_periodicity_constraints.

It should then be not too difficult to create the correct sparsity pattern using this map as new (additional) optional parameter in make_flux_sparsity_pattern. 

We still might want to consider additional constraints, so replacing is no option.

Bests,

Daniel

[Sungho Yoon]

Well.. me neither. I do not yet implement do for distributed mesh for this. In DG, rather than using the constraint matrix in the format in deal.II, connecting the fluxes between the periodic cells would be suitable to DG formulation, since DG is based on the flux. Well. I just did not feel the necessary to do this in distributed mesh, just guessing this one how to implement.

2015년 10월 6일 화요일 오전 7시 35분 13초 UTC-7, Andreas Krämer 님의 말:

[Andreas Krämer]

Hi folks,

based on what you proposed, I have implemented a function
make_periodicity_map_dg that figures out the cells at a periodic boundary recursively.

The files are attached (the other make_sparser_flux_sparsity_pattern can be ignored)

No warranty, but it works fine for my applications.

Best
Andreas

[Daniel Arndt]

Hey Andreas,

I am quite sure that a version of make_flux_sparsity_pattern that accepts additional face pairs is interesting for others.

Why not create a pull request for this? If you need assistance, I happily help.

Bests

Daniel

[Andreas Krämer]

Hi Daniel,

thanks for the kind offer. I have some concers about adding this function to deal.ii, because my version of make_flux_sparsity_pattern depends on the class BoundaryCollection.

Maybe I can provide a low-level interface. If it succeeds, I will let you know.

Best,

Andreas

[Nectarios]

Hi all. I have read Andreas' approach and code with great interest since my application requires DG elements for a problem with periodic boundary conditions. Since I am rather new to deal.ii I would really appreciate it if Andreas or anyone else from the community had a code snippet or some guidelines on how to implement Andreas' approach to a simple example such as a Poisson equation.   

Best Regards,

Nectarios

[Andreas Krämer]

...

[Andreas Krämer]

Hi Nectarios,

the  way I did it works as follows (starting with an unrefined grid):

• Collect periodic faces to create a halo layer of ghost cells across the periodic boundary:   
  // collect periodic faces and extend halo layer of ghost cells across periodic boundaries
      // take care: collect_periodic_faces must not be called for a refined grid
      if (m_triangulation->n_levels() > 1) {
          throw PeriodicBoundaryNotPossible(
                  "The periodic boundaries have to be set before the refinement.");
      }
      std::vector<
              dealii::GridTools::PeriodicFacePair<
                      typename Mesh<dim>::cell_iterator> > matched_pairs;
      dealii::GridTools::collect_periodic_faces(*m_triangulation,
              m_boundaryIndicator1, m_boundaryIndicator2, m_direction,
              matched_pairs);
      // add halo layer
      m_triangulation->add_periodicity(matched_pairs);
  
  
  
  
  
• Refine grid (if you want to) and distribute DoFs
• Create a cell map that identifies faces at the two sides of a periodic boundary
  
  PeriodicCellMap<dim> m_cells;    
  DealIIExtensions::make_periodicity_map_dg(doFHandler, m_boundaryIndicator1, m_boundaryIndicator2, m_direction, m_cells);
  
  
• Make the sparsity pattern, incorporating the PeriodicCellMap (see function make_sparser_flux_sparsity_pattern to get the idea). In general you could just call deal.II's make_flux_sparsity_pattern and iterate over the PeriodicCellMap to add the missing entries.
  
• Assemble the system, incorporating the PeriodicCellMap to identify neighbors across periodic faces. I did something like
  // loop over all cells
  ...
  // loop over all faces
      for (size_t j = 0; j < dealii::GeometryInfo<dim>::faces_per_cell; j++) {
          //Faces at boundary
          if (cell->face(j)->at_boundary()) {
              size_t boundaryIndicator = cell->face(j)->boundary_id();
              if (m_boundaries->isPeriodic(boundaryIndicator)) {
                  // Apply periodic boundaries
                  const boost::shared_ptr<PeriodicBoundary<dim> >& periodicBoundary =
                          m_boundaries->getPeriodicBoundary(boundaryIndicator);
                  assert(periodicBoundary->isFaceInBoundary(cell, j));
                  typename dealii::DoFHandler<dim>::cell_iterator neighborCell;
                  size_t opposite_face =
                          periodicBoundary->getOppositeCellAtPeriodicBoundary(
                                  cell, neighborCell); // this function extracts the neighbor cell from the PeriodicCellMap
                  // now assemble And Distribute Internal Face with the periodic neighbor cell
  
  
  
  The names of the functions should give you an idea what actually happens here.
  

PeriodicCellMap, make_periodicity_map_dg and make_sparser_flux_sparsity pattern are defined in the attached file to give you a rough idea (this version works better than the one posted before).

The code is part of a bigger framework, so you will have to change it quite a bit before using it. However, the most important function, make_periodicity_map_dg, should also work for your code.

If you have further questions, I'm happy to help.

Best,

Andreas

[Nectarios Papanicolaou]

Dear Andreas,

Thank you very much for your help. I will look at it and let you know how things go.

Best Regards,

Nectarios 

--
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 a topic in the Google Groups "deal.II User Group" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/dealii/zdiXh-NRcYk/unsubscribe.
To unsubscribe from this group and all its topics, send an email to dealii+un...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Lukas Korous]

Hello,

I came across this old thread, and successfully used the provided code in my own code (https://github.com/l-korous/mhdeal) - thank you Andreas - and I want to share how I used it:

- here are the (slightly changed) files provided by Andreas:
    https://github.com/l-korous/mhdeal/blob/master/DealiiExtensions.h
https://github.com/l-korous/mhdeal/blob/master/DealiiExtensions.cpp
* I changed the custom class PeriodicBoundary's functionality into simple arrays of boundary_ids and directions - https://github.com/l-korous/mhdeal/blob/master/DealiiExtensions.h#L54

- here is the usage of collect_periodic_faces & add_periodicity: https://github.com/l-korous/mhdeal/blob/master/parameters.cpp#L65
- and here is the usage of the custom functions from Andreas: https://github.com/l-korous/mhdeal/blob/master/problem.cpp#L50

Many thanks