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  X part_unity(gfs, 1);
62 
63  for (auto it = gfs.gridView().template begin<0>(); it != gfs.gridView().template end<0>(); ++it) {
64 
65  lfs.bind(*it);
66 
67  auto geo = it->geometry();
68  const auto gt = geo.type();
69  const auto& ref_el = Dune::ReferenceElements<double, dim>::general(gt);
70 
71  auto& coeffs = lfs.finiteElement().localCoefficients();
72 
73  for (std::size_t i = 0; i < coeffs.size(); ++i) {
74 
75  auto local_pos = ref_el.position (coeffs.localKey(i).subEntity(), coeffs.localKey(i).codim());
76 
77  auto global_pos = geo.global(local_pos);
78 
79  auto subindex = gfs.entitySet().indexSet().subIndex(*it, coeffs.localKey(i).subEntity(), coeffs.localKey(i).codim());
80 
81  double Hx = 1.0 / (double)partition_x;
82  double Hy = 1.0 / (double)partition_y;
83  double hx = (double)overlap / cells_x;
84  double hy = (double)overlap / cells_y;
85 
86  int row = std::floor(my_rank / partition_x);
87  int col = my_rank - partition_x * row;
88 
89  double dx1 = (col + 1) * Hx + hx - global_pos[0];
90  double dx2 = global_pos[0] - (col * Hx - hx);
91 
92  double dy1 = (row + 1) * Hy + hy - global_pos[1];
93  double dy2 = global_pos[1] - (row * Hy - hy);
94 
95  if (row == 0) dy2 = 2*Hy;
96  if (row == partition_y - 1) dy1 = 2*Hy;
97  if (col == 0) dx2 = 2*Hx;
98  if (col == partition_x - 1) dx1 = 2*Hx;
99 
100  native(part_unity)[subindex] = std::min(std::min(std::min(dx1, dx2), dy1), dy2);
101  }
102  }
103 
104  X sum_dists(part_unity);
105  Dune::PDELab::AddDataHandle<GFS,X> addh_dists(gfs,sum_dists);
106  gfs.gridView().communicate(addh_dists,Dune::All_All_Interface,Dune::ForwardCommunication);
107 
108  auto iter_sum = sum_dists.begin();
109  for (auto iter = part_unity.begin(); iter != part_unity.end(); iter++) {
110  if (*iter > 0)
111  *iter *= 1.0 / *iter_sum;
112  iter_sum++;
113  }
114 
115  Dune::PDELab::set_constrained_dofs(cc,0.0,part_unity); // Zero on Dirichlet domain boundary
116 
117  return part_unity;
118 }
119 
120 
121 #endif //DUNE_PDELAB_BACKEND_ISTL_GENEO_PARTITIONOFUNITY_HH
Dune::PDELab::Backend::native
std::enable_if< std::is_base_of< impl::WrapperBase, T >::value, Native< T > & >::type native(T &t)
Definition: backend/interface.hh:192
Dune::PDELab::set_constrained_dofs
void set_constrained_dofs(const CG &cg, typename XG::ElementType x, XG &xg)
construct constraints from given boundary condition function
Definition: constraints.hh:796
sarkisPartitionOfUnity
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
standardPartitionOfUnity
X standardPartitionOfUnity(const GFS &gfs, const CC &cc)
Compute a simple partition of unity.
Definition: partitionofunity.hh:17
Dune::PDELab::AddDataHandle
Definition: genericdatahandle.hh:665
dim
static const int dim
Definition: adaptivity.hh:84