dune-pdelab 2.7-git
Loading...
Searching...
No Matches
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
16template<class X, class GFS, class CC>
17X 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
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) {
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
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