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){
                        //Define physical groups for square
                        //if(model.radius == 0){
                        elemIndx2PG[id] = 1000+i;
                    }
                    //}
                    /*else if(config.get<std::string>("grid_name").compare("wingbox")==0){
                    //first limb
                    double leftpoint = 0;
                    double rightpoint = model.limb_a;
                    if(point[0] >= leftpoint && point[0] < rightpoint) {     elemIndx2PG[id] = 1000+i;  }
                    //first radial section
                    leftpoint = rightpoint;  rightpoint += model.radius;
                    if(point[0] >= leftpoint && point[0] < rightpoint) {     elemIndx2PG[id] = 2000+i;  }
                    //second limb
                    leftpoint = rightpoint;  rightpoint += model.limb_b;
                    if(point[0] >= leftpoint && point[0] < rightpoint) {     elemIndx2PG[id] = 3000+i;  }
                    //second radial section
                    leftpoint = rightpoint;  rightpoint += model.radius;
                    if(point[0] >= leftpoint && point[0] < rightpoint) {     elemIndx2PG[id] = 4000+i;  }
                    //third limb
                    leftpoint = rightpoint;  rightpoint += model.limb_a;
                    if(point[0] >= leftpoint && point[0] < rightpoint) {     elemIndx2PG[id] = 5000+i;  }
                    //third radial section
                    leftpoint = rightpoint;  rightpoint += model.radius;
                    if(point[0] >= leftpoint && point[0] < rightpoint) {     elemIndx2PG[id] = 6000+i;  }
                    //fourth limb
                    leftpoint = rightpoint;  rightpoint += model.limb_b;
                    if(point[0] >= leftpoint && point[0] < rightpoint) {     elemIndx2PG[id] = 7000+i;  }
                    //fourth radial section
                    leftpoint = rightpoint;  rightpoint += model.radius;
                    if(point[0] >= leftpoint && point[0] < rightpoint) {     elemIndx2PG[id] = 8000+i;  }
                    if(config.get<int>("mpc",0) ){
                    for(int c = 0; c < eit.geometry().corners(); c++){
                    auto corner = eit.geometry().corner(c);
                    if(corner[1] > model.Y - 10){
                    elemIndx2PG[id] = 9000+i;
                    }
                    }
                    }
 
                    }
                    //Define physical groups for lshape
                    else{
                    if(model.num_limbs == 2){
                    if(point[0] < model.limb+ 1e-06)                                              {     elemIndx2PG[id] = 1000+i;  }
                    if(point[0] > model.limb+ 1e-06 && point[0] < model.limb+model.radius+ 1e-06) {     elemIndx2PG[id] = 2000+i;  }
                    if(point[0] >= model.limb+model.radius)                                       {     elemIndx2PG[id] = 3000+i;  }
                    if(config.get<int>("mpc",0)){
                    for(int c = 0; c < eit.geometry().corners(); c++){
                    auto corner = eit.geometry().corner(c);
                    if(corner[0] > model.radius+model.num_limbs*model.limb - 1e-6){
                    elemIndx2PG[id] = 4000+i;
                    }
                    }
                    }
                    if(config.get<int>("edge_treatment",0)){
                    if(point[1] < config.get<double>("edge_width",0) || point[1] > model.Y-config.get<double>("edge_width")){
                    elemIndx2PG[id] = 5000; //resin at edges
                    }
                    }
                    }
                    else{
                    if(point[0] <  model.radius+ 1e-06) {     elemIndx2PG[id] = 1000+i;  }
                    if(point[0] >= model.radius+ 1e-06) {     elemIndx2PG[id] = 2000+i;  }
                    if(config.get<int>("mpc",0)){
                    for(int c = 0; c < eit.geometry().corners(); c++){
                    auto corner = eit.geometry().corner(c);
                    if(corner[0] > model.radius+model.limb - 1e-6){
                    elemIndx2PG[id] = 3000+i;
                    }
                    }
                    }
                    }
                    }
                    }*/
                    layer_cnt +=  inc;
            }
        }
    return elemIndx2PG;
}
 
}
}
 
#endif