definePG.hh

#ifndef DEFINE_PG__HH
#define DEFINE_PG__HH
 
namespace Dune{
    namespace Composites
 
        template<typename GV, typename MODEL>
        std::vector<int> setPG(const GV& gv, const MODEL& model, int Nelem){
            std::vector<int> elemIndx2PG(Nelem);
            for (const auto& eit : elements(gv)){
                int id = gv.indexSet().index(eit);
                auto point = eit.geometry().center();
                double layer_cnt=0, inc=0;
                int layer_switch=0;
                for(int i = 0; i < model.num_layers; i++){
                    layer_switch = (i+1)/2 - i/2;
                    inc = model.lay_res_h*layer_switch+(1-layer_switch)*model.lay_com_h;
                    if(point[2] > layer_cnt + 1e-06 && point[2] < layer_cnt + inc + 1.e-06){
                        elemIndx2PG[id] = 1000+i;
                    }
                    layer_cnt +=  inc;
            }
        }
    return elemIndx2PG;
}
 
}
}
 
#endif