[dune-pdelab] Problem with UGGrid-GMSH physical-groups in parallel version
Shubhangi Gupta
sgupta at geomar.de
Mon Oct 24 14:23:24 CEST 2016
Hi Christian,
Thanks a lot for your answer.
Sorry , but I didn't understand exactly what I should do...
So, something like this?
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() ;
if(helper.rank()==0){
Dune::GmshReader<GridType>::read(factory,grid_file,true,true);
factory.createGrid();
grid.loadBalance();
}
// get view
typedef GridType::LeafGridView GV;
const GV& gv=grid.leafGridView();
using ES = Dune::PDELab::NonOverlappingEntitySet<GV>;
And how can I communicate data along the load balance call?
Is there an example in pdelab howto that I can refer to?
Thanks once again,
warm regards, Shubhangi
On 24.10.2016 13:47, Christian Engwer wrote:
> 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
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://lists.dune-project.org/pipermail/dune-pdelab/attachments/20161024/a9d15803/attachment.htm>
More information about the dune-pdelab
mailing list