[dune-pdelab] Problem with UGGrid-GMSH physical-groups in parallel version

[Christian Engwer]

There are two things you have to take care of

a) you use UGGrid. UGGrid must be constructed on rank 0 and the you
   call load-balance to distribute the grid.
b) you collect data from the gmsh file. After load-balance indices
   change and thus you have to take case that the correct material
   properties and the boundary conditions are available on the remote
   rank. This gives two options:
   - you can communicate the data along the loadbalance call
   - you can store the data initially w.r.t. the global ID and not the
     local index, e.g. in a map. Every processor can read the file,
     store the association and then we have all required information
     locally available and remap onto a local index.

Christian

On Mon, Oct 24, 2016 at 11:53:46AM +0200, Shubhangi Gupta wrote:
> Hi all,
> I make grid for my problem using UGGrid with a gmsh file.
> I also prescribe the material properties and the boundary conditions using
> 'physical groups'. I have written the functions below (M and BCType).
> These worked fine in the sequential version, but I am now getting
> segmentation fault for parallel version.
> I use the NonOverlappingEntitySet (similar to the
> convection-diffusion/nonoverlappingsinglephaseflow example in
> pdelab-how-to.)
> Here is how I create the grid:
> 
>         typedef Dune::UGGrid<dim> GridType;
>         GridType grid;
>         // read gmsh file
>         Dune::GridFactory<GridType> factory(&grid);
>         std::vector<int> boundary_index_map;
>         std::vector<int> element_index_map;
>         const std::string grid_file = parameters.getGridName() ;
> Dune::GmshReader<GridType>::read(factory,grid_file,true,true);
>         factory.createGrid();
>         grid.loadBalance();
>         std::cout << " after load balance /" << helper.rank() << "/ " <<
> grid.size(0) << std::endl;
>         //grid.globalRefine(1);
>         // std::cout << " after refinement /" << helper.rank() << "/ " <<
> grid.size(0) << std::endl;
>         // get view
>         typedef GridType::LeafGridView GV;
>         const GV& gv=grid.leafGridView();
>         using ES = Dune::PDELab::NonOverlappingEntitySet<GV>;
>         ES es(gv);
>         problem_driver(es);
> 
> I am not sure why I am getting the seg-fault.
> I will really appreciate if anyone can help me with this issue!
> Thanks a lot in advance,
> warm regards,
> Shubhangi
> 
> PS: Here are the functions M and BCType for reference:
> ___________________________________________________________________________________________________________________
> /** \brief A function for material properties */
> template<typename GV, typename RF, typename PGMap>
> class M : public Dune::PDELab::GridFunctionBase<Dune::PDELab::GridFunctionTraits<GV,RF,1,Dune::FieldVector<RF,1>
> >, M<GV,RF,PGMap> >
> {
> private:
>       const GV& gv;
>       const Dune::MultipleCodimMultipleGeomTypeMapper<GV, P0Layout> mapper ;
>       const PGMap& pg ;
> public:
>   typedef Dune::PDELab::GridFunctionTraits<GV,RF,1,Dune::FieldVector<RF,1> >
> Traits;
>   typedef Dune::PDELab::GridFunctionBase < Dune::PDELab::GridFunctionTraits<
> GV, RF, 1 , Dune::FieldVector< RF,1 > >, M< GV, RF, PGMap > > BaseT ;
> 
>   //! construct from grid view
>   M ( const GV& gv_ , const PGMap& pg_ ) : gv( gv_ ) , mapper( gv_ ) , pg(
> pg_ ) {}
> 
>   //! evaluate extended function on element
>   inline void evaluate (const typename Traits::ElementType& e,
>                                     const typename Traits::DomainType&
> xlocal,
>                                     typename Traits::RangeType& y) const {
> 
>       // retrieve element index and corresponding material index on level 0
>           typename GV::template Codim<0>::EntityPointer ep(e) ;
>           while ( ep->level() != 0 ) ep = ep->father() ;
>           const int ei = mapper.map(ep) ;
>           const int physgroup_index = pg[ ei ] ;
> 
>       // evaluate physical group map and set values accordingly
>           switch ( physgroup_index ){
>                 case 1     : y = y=1.     ; break ;
>                 case 2     : y = y=0.     ; break ;
>                 default     : y = 1.         ; break ; // only one material
> here
>           }
>           return ;
>   }
> 
>   //! get a reference to the grid view
>   inline const GV& getGridView () {return gv;}
> };
> ___________________________________________________________________________________________________________________
> /** \brief A function for boundary types */
> template<typename GV, typename PGMap >
> class BCType : public Dune::PDELab::BoundaryGridFunctionBase<Dune::PDELab::BoundaryGridFunctionTraits<GV,int,1,Dune::FieldVector<int,1>
> >,BCType<GV,PGMap> >
> {
> public :
>     typedef Dune::PDELab::BoundaryGridFunctionTraits<
> GV,int,1,Dune::FieldVector<int,1> > Traits ;
> 
>     // ! construct from gridview
>     BCType ( const GV& gv_ , const PGMap& pgmap_ ) : gv ( gv_ ) , pgmap (
> pgmap_ ) {}
>     // ! return bctype at point on intersection
>     template<typename I>
>     inline void evaluate ( I& i , const typename Traits::DomainType& xlocal
> , typename Traits::RangeType& y ) const {
> 
>          // evaluate with maps
>          int physgroup_index = pgmap[
> i.intersection().boundarySegmentIndex() ] ;
>          switch ( physgroup_index )
>          {
>                  case 3    : y = Indices::BCId_neumann;    break ;
>                  case 4    : y = Indices::BCId_dirichlet;      break ;
>              default : y = Indices::BCId_neumann ; break ; // Neumann
>          }
>          return ;
>     }
> 
>     // ! get a reference to the gridview
>     inline const GV& getGridView () { return gv ; }
> 
> private :
>     const GV& gv ;
>     const PGMap& pgmap ;
> 
> } ;
> _________________________________________________________________________________________________________________
> 

> _______________________________________________
> dune-pdelab mailing list
> dune-pdelab at dune-project.org
> http://lists.dune-project.org/mailman/listinfo/dune-pdelab


-- 
Prof. Dr. Christian Engwer 
Institut für Numerische und Angewandte Mathematik
Fachbereich Mathematik und Informatik der Universität Münster
Einsteinstrasse 62
48149 Münster

E-Mail  christian.engwer at uni-muenster.de
Telefon +49 251 83-35067
FAX     +49 251 83-32729