[dune-pdelab] periodic boundary conditions

[Christian Engwer]

On Wed, Nov 02, 2011 at 01:12:17PM +0100, Oleh Krehel wrote:
> Hi,
> 
> I tried this code:
> 
>   const int d=2;
>   Dune::FieldVector<double,d> L(1);
>   Dune::FieldVector<int,d> N(2);
>   Dune::FieldVector<bool,d> periodic(true);
>   int overlap=0;
        ^^^^^^^^^

>   Dune::YaspGrid<d> grid(L,N,periodic,overlap);
> 
> Then when I traverse on the leafview the instersections of each element,
> the sets of is->neighbor() and is->boundary() don't overlap, they're
> just like the non-periodic ones.

You asked for a grid without overlap

> Can someone please explain, maybe there's something I've missed?
> 
> regards,
> Oleh
> 
> 
> 
> 
> 
> 2011/10/31 Christian Engwer <christian.engwer at uni-muenster.de>:
> > Hi!
> >
> > first you have to settle the question of the grid. I don't have a
> > complete overview of the support of periodic boundaries of the
> > different grid managers.
> >
> > I know:
> > - YaspGrid handles what you want, but without inclusions
> > - UG does have some kind of periodic support, but I don't know what
> >  exactly is supported and you can't use easily fromm DUNE [patches
> >  wellcome ;-)]
> > - Alugrid from 1.23 on support periodic boundaries in 3D
> >
> > Depending on the grid you are using the periodicity is implemented
> > similar to parallel coupling, which uses domain decomposition with
> > overlap or ghosts. In this case you should use the periodic
> > constraints. If the grid implements periodicity as periodicity in the
> > topology of the grid, you don't have to change anything (not sure if
> > any grid supports this).
> >
> > The last option is implement periodictity in PDELab without support of
> > the grid. You should the write your own constraints-assembler
> > local-operator, which adds dirichlet constraints on one side. There is
> > no real tutorial on this. You will have to dive into it...
> >
> > Christian
> >
> >> I need a fully periodic cube: left-right, top-bottom in 2D,
> >> and additionally front-back in 3D.
> >> The cube has an arbitrarily shaped inclusion(like a circle or a sphere) inside,
> >> with Dirichlet b.c.
> >> I was looking at section 8.4 of the tutorial explaining constraints
> >> for Dirichlet b.c.,
> >> but the transformation matrix there is local, when I need the global one for the
> >> periodic b.c.
> >> Can you point me to the source file that does what you described above?
> >>
> >> regards,
> >> Oleh
> >>
> >>
> >>
> >> 2011/10/31 Christian Engwer <christian.engwer at uni-muenster.de>:
> >> > Dear Oleh,
> >> >
> >> >> I need to implement periodic boundary conditions for a transport problem.
> >> >> Google shows that there's no such functionality in DUNE at the moment.
> >> >> Perhaps someone can give me some pointers or suggestions on how
> >> >> to implement periodic b.c. in a straightforward way.
> >> >
> >> > it depends on which BC you need.
> >> >
> >> > Straight forward: You have a mesh like this:
> >> >
> >> > o-x-x-x-o
> >> > | | | | |
> >> > o-x-x-x-o
> >> > | | | | |
> >> > o-x-x-x-o
> >> >
> >> > where the 'o' nodes are associated with each other. This kind of mesh
> >> > is available with yaspgrid. It is implemented as a periodic mesh with
> >> > overlap on the left and right side.
> >> >
> >> > If you need something more general it will get tricky and you should
> >> > explain in detail what you want to achieve.
> >> >
> >> > Cheers
> >> > Christian
> >> >
> >> >
> >>
> >

> #define DUNE_DEVEL_MODE
> //#define NDEBUG
> #include "../config.h"
> #include <iostream>
> #include <dune/common/mpihelper.hh>
> #include <dune/common/exceptions.hh>
> #include <boost/exception/all.hpp>
> #include <dune/grid/yaspgrid.hh>
> 
> template <class GV>
> void traverseGrid(const GV& gv){
>   typedef typename GV::template Codim<0>::Iterator element_iterator;
>   typedef typename GV::IntersectionIterator intersection_iterator;
>   element_iterator it=gv.template begin<0>(),
>                   end=gv.template end<0>();
>   for(;it!=end;++it){
>     Dune::GeometryType gt=it->type();
>     std::cout<<"\nvisiting leaf "<<gt
>              <<" with first vertex at "<<it->geometry().corner(0);
>     intersection_iterator is=gv.ibegin(*it),
>                        isend=gv.iend(*it);
>     for(;is!=isend;++is){
>       if(is->neighbor())
>         std::cout<<"\nneighbor";
>       if(is->boundary())
>         std::cout<<"\nboundary";
>     }
>   }
> }
> 
> void dowork(int level){
>   const int d=2;
>   Dune::FieldVector<double,d> L(1);
>   Dune::FieldVector<int,d> N(2);
>   Dune::FieldVector<bool,d> periodic(true);
>   int overlap=0;
>   Dune::YaspGrid<d> grid(L,N,periodic,overlap);
>   //grid.globalRefine(level);
>   typedef Dune::YaspGrid<d>::LeafGridView GV;
>   const GV& gv=grid.leafView();
>   traverseGrid(gv);
> }
> 
> int main(int argc, char** argv)
>     try
>     {
>       Dune::MPIHelper::instance(argc,argv);
>       int level=2;
>       if(argc==2)
>         sscanf(argv[1],"%d",&level);
>       dowork(level);
>       std::cout<<"\n;\n";
>       return 0;
>     }
>     catch (Dune::Exception& e){
>       std::cerr << "Dune reported error: " << e << std::endl;
>     }
>     catch (boost::exception& e){
>       std::cerr << boost::diagnostic_information(e);
>     }
>     catch (...){
>       std::cerr << "Unknown exception thrown!" << std::endl;
>     }