partitionofunity.hh

#ifndef DUNE_GENEO_PARTITIONOFUNITY_HH
#define DUNE_GENEO_PARTITIONOFUNITY_HH
 
namespace Dune{
    namespace Geneo{
 
        template<const int dim, class X, class GFS,  class LFS, class CC>
            std::shared_ptr<X> standardPartitionOfUnity(const GFS& gfs, LFS& lfs, const CC& cc) {
 
                auto part_unity = std::make_shared<X>(gfs, 1);
 
                Dune::PDELab::set_constrained_dofs(cc,0.0,*part_unity); // Zero on subdomain boundary
 
                Dune::PDELab::AddDataHandle<GFS,X> parth(gfs,*part_unity);
                gfs.gridView().communicate(parth,Dune::All_All_Interface,Dune::ForwardCommunication);
 
                Dune::PDELab::set_constrained_dofs(cc,0.0,*part_unity); // Zero on subdomain boundary (Need that a 2nd time due to add comm before!)
 
                for (auto iter = part_unity->begin(); iter != part_unity->end(); iter++) {
                    if (*iter > 0)
                        *iter = 1.0 / *iter;
                }
                return part_unity;
            }
 
        template<const int dim, class X, class GFS,  class LFS, class CC>
            std::shared_ptr<X> sarkisPartitionOfUnity(const GFS& gfs, LFS& lfs, const CC& cc) {
                using Dune::PDELab::Backend::native;
 
                int my_rank = gfs.gridView().comm().rank();
 
                auto part_unity = std::make_shared<X>(gfs, 1);
 
                int cells = configuration.get<int>("Grid.cells", 10);
                int overlap = configuration.get<int>("Grid.overlap", 1);
                int partition_x = configuration.get<int>("Grid.partition_x", 1);
                int partition_y = configuration.get<int>("Grid.partition_y", 1);
 
                for (auto it = gfs.gridView().template begin<0>(); it != gfs.gridView().template end<0>(); ++it) {
 
                    lfs.bind(*it);
 
                    auto geo  = it->geometry();
                    const auto gt = geo.type();
                    const auto& ref_el = Dune::ReferenceElements<double, dim>::general(gt);
 
                    auto& coeffs = lfs.finiteElement().localCoefficients();
 
                    for (std::size_t i = 0; i < coeffs.size(); ++i) {
 
                        auto local_pos = ref_el.position (coeffs.localKey(i).subEntity(), coeffs.localKey(i).codim());
 
                        auto global_pos = geo.global(local_pos);
 
                        auto subindex = gfs.entitySet().indexSet().subIndex(*it, coeffs.localKey(i).subEntity(), coeffs.localKey(i).codim());
 
                        double Hx = 1.0 / (double)partition_x;
                        double Hy = 1.0 / (double)partition_y;
                        double h = (double)overlap / cells;
 
                        int row = std::floor(my_rank / partition_x);
                        int col = my_rank - partition_x * row;
 
                        double dx1 = (col + 1) * Hx + h - global_pos[0];
                        double dx2 = global_pos[0] - (col * Hx - h);
 
                        double dy1 = (row + 1) * Hy + h - global_pos[1];
                        double dy2 = global_pos[1] - (row * Hy - h);
 
                        if (row == 0) dy2 = 2*Hy;
                        if (row == partition_y - 1) dy1 = 2*Hy;
                        if (col == 0) dx2 = 2*Hx;
                        if (col == partition_x - 1) dx1 = 2*Hx;
 
                        native(*part_unity)[subindex] = std::min(std::min(std::min(dx1, dx2), dy1), dy2);
                    }
                }
 
                X sum_dists(gfs, 0.0);
                sum_dists = *part_unity;
                Dune::PDELab::AddDataHandle<GFS,X> addh_dists(gfs,sum_dists);
                gfs.gridView().communicate(addh_dists,Dune::All_All_Interface,Dune::ForwardCommunication);
 
                auto iter_sum = sum_dists.begin();
                for (auto iter = part_unity->begin(); iter != part_unity->end(); iter++) {
                    if (*iter > 0)
                        *iter *= 1.0 / *iter_sum;
                    iter_sum++;
                }
 
                Dune::PDELab::set_constrained_dofs(cc,0.0,*part_unity); // Zero on Dirichlet domain boundary
 
                return part_unity;
            }
 
    }
}
 
 
#endif //DUNE_GENEO_PARTITIONOFUNITY_HH