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

[Shubhangi Gupta]

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 ;

} ;
_________________________________________________________________________________________________________________ 

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://lists.dune-project.org/pipermail/dune-pdelab/attachments/20161024/1c0650b9/attachment.htm>