AMDiS 2.10
The Adaptive Multi-Dimensional Simulation Toolbox
PeriodicBC.inc.hpp
1#pragma once
2
3#include <algorithm>
4#include <array>
5#include <cmath>
6#include <limits>
7
8#include <dune/functions/functionspacebases/subentitydofs.hh>
9
10#include <amdis/LinearAlgebra.hpp>
11#include <amdis/Output.hpp>
12#include <amdis/typetree/Traversal.hpp>
13
14namespace AMDiS {
15
16template <class Basis, class TP>
18init()
19{
20 if (!bool(hasConnectedIntersections_)) {
21 hasConnectedIntersections_ = true;
22 const auto& gridView = basis_.gridView();
23 for (auto const& element : elements(gridView)) {
24 for (const auto& intersection : intersections(gridView, element)) {
25 if (!boundarySubset_(intersection))
26 continue;
27
28 if (!intersection.neighbor()) {
29 hasConnectedIntersections_ = false;
30 break;
31 }
32 }
33 }
34 }
35
36 if (*hasConnectedIntersections_)
37 initAssociations();
38 else
39 initAssociations2();
40}
41
42
43template <class Basis, class TP>
46{
47 using std::sqrt;
48 static const auto tol = sqrt(std::numeric_limits<typename Domain::field_type>::epsilon());
49 periodicNodes_.assign(basis_.dimension(), false);
50
51 auto localView = basis_.localView();
52 auto seDOFs = Dune::Functions::subEntityDOFs(basis_);
53 const auto& gridView = basis_.gridView();
54 for (auto const& element : elements(gridView)) {
55 if (!element.hasBoundaryIntersections())
56 continue;
57
58 localView.bind(element);
59 for (const auto& intersection : intersections(gridView, element)) {
60 if (!boundarySubset_(intersection))
61 continue;
62
63 test_exit(intersection.neighbor(),
64 "Neighbors of periodic intersection not assigned. Use a periodic Grid instead!");
65
66 // collect DOFs of inside element on intersection
67 localView.bind(intersection.inside());
68 seDOFs.bind(localView, intersection.indexInInside(), 1);
69 std::vector<std::size_t> insideDOFs(seDOFs.begin(), seDOFs.end());
70 std::vector<MultiIndex> insideGlobalDOFs(insideDOFs.size());
71 std::transform(insideDOFs.begin(), insideDOFs.end(), insideGlobalDOFs.begin(),
72 [&](std::size_t localIndex) { return localView.index(localIndex); });
73 auto insideCoords = coords(localView.tree(), insideDOFs);
74
75 // collect DOFs of ouside element on intersection
76 localView.bind(intersection.outside());
77 seDOFs.bind(localView, intersection.indexInOutside(), 1);
78 std::vector<std::size_t> outsideDOFs(seDOFs.begin(), seDOFs.end());
79 auto outsideCoords = coords(localView.tree(), outsideDOFs);
80
81 // compare mapped coords of DOFs
82 assert(insideDOFs.size() == outsideDOFs.size());
83 for (std::size_t i = 0; i < insideCoords.size(); ++i) {
84 auto x = faceTrafo_.evaluate(insideCoords[i]);
85 for (std::size_t j = 0; j < outsideCoords.size(); ++j) {
86 auto const& y = outsideCoords[j];
87 if (distance(x,y) < tol) {
88 periodicNodes_[insideGlobalDOFs[i]] = true;
89 associations_[insideGlobalDOFs[i]] = localView.index(outsideDOFs[j]);
90 }
91 }
92 }
93 }
94 }
95}
96
97namespace Impl
98{
99 template <class D>
100 class CompareTol
101 {
102 using ctype = typename D::field_type;
103
104 public:
105 CompareTol(ctype tol)
106 : tol_(tol)
107 {}
108
109 bool operator()(D const& lhs, D const& rhs) const
110 {
111 using std::abs;
112 for (int i = 0; i < D::dimension; i++) {
113 if (abs(lhs[i] - rhs[i]) < tol_)
114 continue;
115 return lhs[i] < rhs[i];
116 }
117 return false;
118 }
119
120 private:
121 ctype tol_;
122 };
123}
124
125
126template <class Basis, class TP>
127void PeriodicBC<Basis, TP>::
128initAssociations2()
129{
130 using std::sqrt;
131 static const auto tol = sqrt(std::numeric_limits<typename Domain::field_type>::epsilon());
132 periodicNodes_.assign(basis_.dimension(), false);
133
134 // Dune::ReservedVector<std::pair<EntitySeed,int>,2>
135 using EntitySeed = typename SubspaceBasis::GridView::Grid::template Codim<0>::EntitySeed;
136 using CoordAssoc = std::map<Domain, std::vector<std::pair<EntitySeed,int>>, Impl::CompareTol<Domain>>;
137 CoordAssoc coordAssoc(Impl::CompareTol<Domain>{tol});
138
139 // generate periodic element pairs
140 const auto& gridView = basis_.gridView();
141 for (auto const& element : elements(gridView)) {
142 for (const auto& intersection : intersections(gridView, element)) {
143 if (!boundarySubset_(intersection))
144 continue;
145
146 auto x = intersection.geometry().center();
147 auto seed = intersection.inside().seed();
148 coordAssoc[x].push_back({seed,intersection.indexInInside()});
149 coordAssoc[faceTrafo_.evaluate(x)].push_back({seed,intersection.indexInInside()});
150 }
151 }
152
153 auto localView = basis_.localView();
154 auto seDOFs = Dune::Functions::subEntityDOFs(basis_);
155 for (auto const& assoc : coordAssoc) {
156 auto const& seeds = assoc.second;
157 if (seeds.size() != 2)
158 continue;
159
160 // collect DOFs of inside element on intersection
161 auto el1 = gridView.grid().entity(seeds[0].first);
162 localView.bind(el1);
163 seDOFs.bind(localView, seeds[0].second, 1);
164 std::vector<std::size_t> insideDOFs(seDOFs.begin(), seDOFs.end());
165 std::vector<MultiIndex> insideGlobalDOFs(insideDOFs.size());
166 std::transform(insideDOFs.begin(), insideDOFs.end(), insideGlobalDOFs.begin(),
167 [&](std::size_t localIndex) { return localView.index(localIndex); });
168 auto insideCoords = coords(localView.tree(), insideDOFs);
169
170 // collect DOFs of outside element on intersection
171 auto el2 = gridView.grid().entity(seeds[1].first);
172 localView.bind(el2);
173 seDOFs.bind(localView, seeds[1].second, 1);
174 std::vector<std::size_t> outsideDOFs(seDOFs.begin(), seDOFs.end());
175 auto outsideCoords = coords(localView.tree(), outsideDOFs);
176
177 // compare mapped coords of DOFs
178 assert(insideDOFs.size() == outsideDOFs.size());
179 for (std::size_t i = 0; i < insideCoords.size(); ++i) {
180 auto x = faceTrafo_.evaluate(insideCoords[i]);
181 for (std::size_t j = 0; j < outsideCoords.size(); ++j) {
182 auto const& y = outsideCoords[j];
183 if (distance(x,y) < tol) {
184 periodicNodes_[insideGlobalDOFs[i]] = true;
185 associations_[insideGlobalDOFs[i]] = localView.index(outsideDOFs[j]);
186 }
187 }
188 }
189 }
190}
191
192
193template <class Basis, class TP>
194 template <class Node>
195auto PeriodicBC<Basis, TP>::
196coords(Node const& tree, std::vector<std::size_t> const& localIndices) const
197{
198 std::vector<Domain> dofCoords(localIndices.size());
199 Traversal::forEachLeafNode(tree, [&](auto const& node, auto&&)
200 {
201 std::size_t size = node.finiteElement().size();
202 auto geometry = node.element().geometry();
203
204 auto const& localInterpol = node.finiteElement().localInterpolation();
205 using FiniteElement = TYPEOF(node.finiteElement());
206 using DomainType = typename FiniteElement::Traits::LocalBasisType::Traits::DomainType;
207 using RangeType = typename FiniteElement::Traits::LocalBasisType::Traits::RangeType;
208
209 std::array<std::vector<typename RangeType::field_type>, Domain::dimension> coeffs;
210 for (int d = 0; d < Domain::dimension; ++d) {
211 auto evalCoord = [&](DomainType const& local) -> RangeType { return geometry.global(local)[d]; };
212 localInterpol.interpolate(evalCoord, coeffs[d]);
213 }
214
215 for (std::size_t j = 0; j < localIndices.size(); ++j) {
216 Domain x;
217 for (std::size_t i = 0; i < size; ++i) {
218 if (node.localIndex(i) == localIndices[j]) {
219 for (int d = 0; d < Domain::dimension; ++d)
220 x[d] = coeffs[d][i];
221 break;
222 }
223 }
224 dofCoords[j] = std::move(x);
225 }
226 });
227
228 return dofCoords;
229}
230
231
232template <class Basis, class TP>
233 template <class Mat, class Sol, class Rhs>
235apply(Mat& matrix, Sol& solution, Rhs& rhs)
236{
237 periodicBC(matrix, solution, rhs, periodicNodes_, associations_);
238}
239
240} // end namespace AMDiS
Implements a periodic boundary condition.
Definition: PeriodicBC.hpp:62
void init()
Compute the pairs of associated boundary DOFs.
Definition: PeriodicBC.inc.hpp:18
void apply(Mat &matrix, Sol &solution, Rhs &rhs)
Move matrix rows (and columns) of dependant DOFs to the corresponding master DOFs.
Definition: PeriodicBC.inc.hpp:235