1#ifndef DUNE_PDELAB_BACKEND_ISTL_GENEO_PARTITIONOFUNITY_HH
2#define DUNE_PDELAB_BACKEND_ISTL_GENEO_PARTITIONOFUNITY_HH
16template<
class X,
class GFS,
class CC>
24 gfs.gridView().communicate(parth,Dune::All_All_Interface,Dune::ForwardCommunication);
28 for (
auto iter = part_unity.begin(); iter != part_unity.end(); iter++) {
54template<
class X,
class GFS,
class LFS,
class CC>
55X
sarkisPartitionOfUnity(
const GFS& gfs, LFS& lfs,
const CC& cc,
int cells_x,
int cells_y,
int overlap,
int partition_x,
int partition_y) {
58 int my_rank = gfs.gridView().comm().rank();
62 typedef typename GFS::Traits::GridView::Grid::ctype DF;
64 if (gfs.gridView().grid().dimension !=
dim)
65 DUNE_THROW(Dune::NotImplemented,
"Currently, Sarkis partition of unity is only supported for 2 dimensional grids.");
67 if (GFS::Traits::FiniteElement::Traits::LocalBasisType::order() != 1)
68 DUNE_THROW(Dune::NotImplemented,
"Currently, Sarkis partition of unity is only supported for polynomial bases of order 1");
71 if (gfs.gridView().grid().domainSize() != Dune::FieldVector<DF, 2>(1.0))
72 DUNE_THROW(Dune::NotImplemented,
"Currently, Sarkis partition of unity is only supported for the (0,1) x (0,1) unit square.");
76 for (
auto it = gfs.gridView().template begin<0>(); it != gfs.gridView().template end<0>(); ++it) {
80 auto geo = it->geometry();
81 const auto gt = geo.type();
82 const auto& ref_el = Dune::ReferenceElements<double, dim>::general(gt);
84 auto& coeffs = lfs.finiteElement().localCoefficients();
86 for (std::size_t i = 0; i < coeffs.size(); ++i) {
88 auto local_pos = ref_el.position (coeffs.localKey(i).subEntity(), coeffs.localKey(i).codim());
90 auto global_pos = geo.global(local_pos);
92 auto subindex = gfs.entitySet().indexSet().subIndex(*it, coeffs.localKey(i).subEntity(), coeffs.localKey(i).codim());
94 double Hx = 1.0 / (double)partition_x;
95 double Hy = 1.0 / (double)partition_y;
96 double hx = (double)overlap / cells_x;
97 double hy = (double)overlap / cells_y;
99 int row = std::floor(my_rank / partition_x);
100 int col = my_rank - partition_x * row;
102 double dx1 = (col + 1) * Hx + hx - global_pos[0];
103 double dx2 = global_pos[0] - (col * Hx - hx);
105 double dy1 = (row + 1) * Hy + hy - global_pos[1];
106 double dy2 = global_pos[1] - (row * Hy - hy);
108 if (row == 0) dy2 = 2*Hy;
109 if (row == partition_y - 1) dy1 = 2*Hy;
110 if (col == 0) dx2 = 2*Hx;
111 if (col == partition_x - 1) dx1 = 2*Hx;
113 native(part_unity)[subindex] = std::min(std::min(std::min(dx1, dx2), dy1), dy2);
117 X sum_dists(part_unity);
119 gfs.gridView().communicate(addh_dists,Dune::All_All_Interface,Dune::ForwardCommunication);
121 auto iter_sum = sum_dists.begin();
122 for (
auto iter = part_unity.begin(); iter != part_unity.end(); iter++) {
124 *iter *= 1.0 / *iter_sum;
static const int dim
Definition: adaptivity.hh:84
X standardPartitionOfUnity(const GFS &gfs, const CC &cc)
Compute a simple partition of unity.
Definition: partitionofunity.hh:17
X sarkisPartitionOfUnity(const GFS &gfs, LFS &lfs, const CC &cc, int cells_x, int cells_y, int overlap, int partition_x, int partition_y)
Compute a partition of unity according to Sarkis.
Definition: partitionofunity.hh:55
void set_constrained_dofs(const CG &cg, typename XG::ElementType x, XG &xg)
construct constraints from given boundary condition function
Definition: constraints.hh:796
std::enable_if< std::is_base_of< impl::WrapperBase, T >::value, Native< T > & >::type native(T &t)
Definition: backend/interface.hh:192
Definition: genericdatahandle.hh:667