dune-pdelab  2.7-git
partitionofunity.hh
Go to the documentation of this file.
1 #ifndef DUNE_PDELAB_BACKEND_ISTL_GENEO_PARTITIONOFUNITY_HH
2 #define DUNE_PDELAB_BACKEND_ISTL_GENEO_PARTITIONOFUNITY_HH
3 
16 template<class X, class GFS, class CC>
17 X standardPartitionOfUnity(const GFS& gfs, const CC& cc) {
18 
19  X part_unity(gfs, 1);
20 
21  Dune::PDELab::set_constrained_dofs(cc,0.0,part_unity); // Zero on subdomain boundary
22 
23  Dune::PDELab::AddDataHandle<GFS,X> parth(gfs,part_unity);
24  gfs.gridView().communicate(parth,Dune::All_All_Interface,Dune::ForwardCommunication);
25 
26  Dune::PDELab::set_constrained_dofs(cc,0.0,part_unity); // Zero on subdomain boundary (Need that a 2nd time due to add comm before!)
27 
28  for (auto iter = part_unity.begin(); iter != part_unity.end(); iter++) {
29  if (*iter > 0)
30  *iter = 1.0 / *iter;
31  }
32  return part_unity;
33 }
34 
54 template<class X, class GFS, class LFS, class CC>
55 X sarkisPartitionOfUnity(const GFS& gfs, LFS& lfs, const CC& cc, int cells_x, int cells_y, int overlap, int partition_x, int partition_y) {
57 
58  int my_rank = gfs.gridView().comm().rank();
59  const int dim = 2;
60 
61  // throw exception if the setup is currently not supported by the Sarkis partition of unity
62  typedef typename GFS::Traits::GridView::Grid::ctype DF;
63 
64  if (gfs.gridView().grid().dimension != dim)
65  DUNE_THROW(Dune::NotImplemented, "Currently, Sarkis partition of unity is only supported for 2 dimensional grids.");
66 
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");
69 
70  // we assume the grid is instantiated as a (0.0, 1.0) x (0.0, 1.0) grid, so we do an exact flaoting point comparison here
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.");
73 
74  X part_unity(gfs, 1);
75 
76  for (auto it = gfs.gridView().template begin<0>(); it != gfs.gridView().template end<0>(); ++it) {
77 
78  lfs.bind(*it);
79 
80  auto geo = it->geometry();
81  const auto gt = geo.type();
82  const auto& ref_el = Dune::ReferenceElements<double, dim>::general(gt);
83 
84  auto& coeffs = lfs.finiteElement().localCoefficients();
85 
86  for (std::size_t i = 0; i < coeffs.size(); ++i) {
87 
88  auto local_pos = ref_el.position (coeffs.localKey(i).subEntity(), coeffs.localKey(i).codim());
89 
90  auto global_pos = geo.global(local_pos);
91 
92  auto subindex = gfs.entitySet().indexSet().subIndex(*it, coeffs.localKey(i).subEntity(), coeffs.localKey(i).codim());
93 
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;
98 
99  int row = std::floor(my_rank / partition_x);
100  int col = my_rank - partition_x * row;
101 
102  double dx1 = (col + 1) * Hx + hx - global_pos[0];
103  double dx2 = global_pos[0] - (col * Hx - hx);
104 
105  double dy1 = (row + 1) * Hy + hy - global_pos[1];
106  double dy2 = global_pos[1] - (row * Hy - hy);
107 
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;
112 
113  native(part_unity)[subindex] = std::min(std::min(std::min(dx1, dx2), dy1), dy2);
114  }
115  }
116 
117  X sum_dists(part_unity);
118  Dune::PDELab::AddDataHandle<GFS,X> addh_dists(gfs,sum_dists);
119  gfs.gridView().communicate(addh_dists,Dune::All_All_Interface,Dune::ForwardCommunication);
120 
121  auto iter_sum = sum_dists.begin();
122  for (auto iter = part_unity.begin(); iter != part_unity.end(); iter++) {
123  if (*iter > 0)
124  *iter *= 1.0 / *iter_sum;
125  iter_sum++;
126  }
127 
128  Dune::PDELab::set_constrained_dofs(cc,0.0,part_unity); // Zero on Dirichlet domain boundary
129 
130  return part_unity;
131 }
132 
133 
134 #endif //DUNE_PDELAB_BACKEND_ISTL_GENEO_PARTITIONOFUNITY_HH
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