PeriodicBC.inc.hpp 7.77 KB
Newer Older
1
2
#pragma once

3
4
5
#include <algorithm>
#include <array>
#include <cmath>
6
7
#include <limits>

8
9
#include <dune/functions/functionspacebases/subentitydofs.hh>

10
#include <amdis/LinearAlgebra.hpp>
11
#include <amdis/Output.hpp>
12
#include <amdis/functions/FunctionFromCallable.hpp>
13
#include <amdis/typetree/Traversal.hpp>
14
15
16

namespace AMDiS {

Praetorius, Simon's avatar
Praetorius, Simon committed
17
18
template <class Basis, class TP>
void PeriodicBC<Basis, TP>::
19
init()
20
21
22
{
  if (!bool(hasConnectedIntersections_)) {
    hasConnectedIntersections_ = true;
23
    const auto& gridView = basis_.gridView();
24
25
    for (auto const& element : elements(gridView)) {
      for (const auto& intersection : intersections(gridView, element)) {
26
        if (!boundarySubset_(intersection))
27
28
29
30
31
32
33
34
35
36
37
          continue;

        if (!intersection.neighbor()) {
          hasConnectedIntersections_ = false;
          break;
        }
      }
    }
  }

  if (*hasConnectedIntersections_)
38
    initAssociations();
39
  else
40
    initAssociations2();
41
42
43
}


Praetorius, Simon's avatar
Praetorius, Simon committed
44
45
template <class Basis, class TP>
void PeriodicBC<Basis, TP>::
46
initAssociations()
47
{
48
  using std::sqrt;
49
50
  static const auto tol = sqrt(std::numeric_limits<typename Domain::field_type>::epsilon());
  periodicNodes_.assign(basis_.dimension(), false);
51

52
53
54
  auto localView = basis_.localView();
  auto seDOFs = Dune::Functions::subEntityDOFs(basis_);
  const auto& gridView = basis_.gridView();
55
56
57
58
59
60
  for (auto const& element : elements(gridView)) {
    if (!element.hasBoundaryIntersections())
      continue;

    localView.bind(element);
    for (const auto& intersection : intersections(gridView, element)) {
61
      if (!boundarySubset_(intersection))
62
63
64
65
66
67
68
69
70
        continue;

      test_exit(intersection.neighbor(),
        "Neighbors of periodic intersection not assigned. Use a periodic Grid instead!");

      // collect DOFs of inside element on intersection
      localView.bind(intersection.inside());
      seDOFs.bind(localView, intersection.indexInInside(), 1);
      std::vector<std::size_t> insideDOFs(seDOFs.begin(), seDOFs.end());
71
      std::vector<MultiIndex> insideGlobalDOFs(insideDOFs.size());
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
      std::transform(insideDOFs.begin(), insideDOFs.end(), insideGlobalDOFs.begin(),
        [&](std::size_t localIndex) { return localView.index(localIndex); });
      auto insideCoords = coords(localView.tree(), insideDOFs);

      // collect DOFs of ouside element on intersection
      localView.bind(intersection.outside());
      seDOFs.bind(localView, intersection.indexInOutside(), 1);
      std::vector<std::size_t> outsideDOFs(seDOFs.begin(), seDOFs.end());
      auto outsideCoords = coords(localView.tree(), outsideDOFs);

      // compare mapped coords of DOFs
      assert(insideDOFs.size() == outsideDOFs.size());
      for (std::size_t i = 0; i < insideCoords.size(); ++i) {
        auto x = faceTrafo_.evaluate(insideCoords[i]);
        for (std::size_t j = 0; j < outsideCoords.size(); ++j) {
          auto const& y = outsideCoords[j];
88
          if (distance(x,y) < tol) {
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
            periodicNodes_[insideGlobalDOFs[i]] = true;
            associations_[insideGlobalDOFs[i]] = localView.index(outsideDOFs[j]);
          }
        }
      }
    }
  }
}

namespace Impl
{
  template <class D>
  class CompareTol
  {
    using ctype = typename D::field_type;

  public:
    CompareTol(ctype tol)
      : tol_(tol)
    {}

    bool operator()(D const& lhs, D const& rhs) const
    {
112
      using std::abs;
113
      for (int i = 0; i < D::dimension; i++) {
114
        if (abs(lhs[i] - rhs[i]) < tol_)
115
116
117
118
119
120
121
122
123
124
125
126
          continue;
        return lhs[i] < rhs[i];
      }
      return false;
    }

  private:
    ctype tol_;
  };
}


Praetorius, Simon's avatar
Praetorius, Simon committed
127
128
template <class Basis, class TP>
void PeriodicBC<Basis, TP>::
129
initAssociations2()
130
{
131
  using std::sqrt;
132
133
  static const auto tol = sqrt(std::numeric_limits<typename Domain::field_type>::epsilon());
  periodicNodes_.assign(basis_.dimension(), false);
134
135

  // Dune::ReservedVector<std::pair<EntitySeed,int>,2>
136
137
138
  using EntitySeed = typename SubspaceBasis::GridView::Grid::template Codim<0>::EntitySeed;
  using CoordAssoc = std::map<Domain, std::vector<std::pair<EntitySeed,int>>, Impl::CompareTol<Domain>>;
  CoordAssoc coordAssoc(Impl::CompareTol<Domain>{tol});
139
140

  // generate periodic element pairs
141
  const auto& gridView = basis_.gridView();
142
143
  for (auto const& element : elements(gridView)) {
    for (const auto& intersection : intersections(gridView, element)) {
144
      if (!boundarySubset_(intersection))
145
146
147
148
149
150
151
152
153
        continue;

      auto x = intersection.geometry().center();
      auto seed = intersection.inside().seed();
      coordAssoc[x].push_back({seed,intersection.indexInInside()});
      coordAssoc[faceTrafo_.evaluate(x)].push_back({seed,intersection.indexInInside()});
    }
  }

154
155
  auto localView = basis_.localView();
  auto seDOFs = Dune::Functions::subEntityDOFs(basis_);
156
157
158
159
160
161
162
163
164
165
  for (auto const& assoc : coordAssoc) {
    auto const& seeds = assoc.second;
    if (seeds.size() != 2)
      continue;

    // collect DOFs of inside element on intersection
    auto el1 = gridView.grid().entity(seeds[0].first);
    localView.bind(el1);
    seDOFs.bind(localView, seeds[0].second, 1);
    std::vector<std::size_t> insideDOFs(seDOFs.begin(), seDOFs.end());
166
    std::vector<MultiIndex> insideGlobalDOFs(insideDOFs.size());
167
168
169
170
    std::transform(insideDOFs.begin(), insideDOFs.end(), insideGlobalDOFs.begin(),
      [&](std::size_t localIndex) { return localView.index(localIndex); });
    auto insideCoords = coords(localView.tree(), insideDOFs);

171
    // collect DOFs of outside element on intersection
172
173
174
175
176
177
178
179
180
181
182
183
    auto el2 = gridView.grid().entity(seeds[1].first);
    localView.bind(el2);
    seDOFs.bind(localView, seeds[1].second, 1);
    std::vector<std::size_t> outsideDOFs(seDOFs.begin(), seDOFs.end());
    auto outsideCoords = coords(localView.tree(), outsideDOFs);

    // compare mapped coords of DOFs
    assert(insideDOFs.size() == outsideDOFs.size());
    for (std::size_t i = 0; i < insideCoords.size(); ++i) {
      auto x = faceTrafo_.evaluate(insideCoords[i]);
      for (std::size_t j = 0; j < outsideCoords.size(); ++j) {
        auto const& y = outsideCoords[j];
184
        if (distance(x,y) < tol) {
185
186
187
188
189
190
191
192
193
          periodicNodes_[insideGlobalDOFs[i]] = true;
          associations_[insideGlobalDOFs[i]] = localView.index(outsideDOFs[j]);
        }
      }
    }
  }
}


Praetorius, Simon's avatar
Praetorius, Simon committed
194
template <class Basis, class TP>
195
  template <class Node>
Praetorius, Simon's avatar
Praetorius, Simon committed
196
auto PeriodicBC<Basis, TP>::
197
198
coords(Node const& tree, std::vector<std::size_t> const& localIndices) const
{
199
  std::vector<Domain> dofCoords(localIndices.size());
200
  Traversal::forEachLeafNode(tree, [&](auto const& node, auto&&)
201
202
203
204
205
  {
    std::size_t size = node.finiteElement().size();
    auto geometry = node.element().geometry();

    auto const& localInterpol = node.finiteElement().localInterpolation();
206
    using FiniteElement = TYPEOF(node.finiteElement());
207
208
209
    using DomainType = typename FiniteElement::Traits::LocalBasisType::Traits::DomainType;
    using RangeType = typename FiniteElement::Traits::LocalBasisType::Traits::RangeType;

210
211
    std::array<std::vector<typename RangeType::field_type>, Domain::dimension> coeffs;
    for (int d = 0; d < Domain::dimension; ++d) {
212
213
214
      auto evalCoord = [&](DomainType const& local) -> RangeType { return geometry.global(local)[d]; };
      auto evalCoordFct = functionFromCallable<RangeType(DomainType)>(evalCoord);
      localInterpol.interpolate(evalCoordFct, coeffs[d]);
215
216
217
    }

    for (std::size_t j = 0; j < localIndices.size(); ++j) {
218
      Domain x;
219
220
      for (std::size_t i = 0; i < size; ++i) {
        if (node.localIndex(i) == localIndices[j]) {
221
          for (int d = 0; d < Domain::dimension; ++d)
222
223
224
225
226
227
228
229
230
231
232
233
            x[d] = coeffs[d][i];
          break;
        }
      }
      dofCoords[j] = std::move(x);
    }
  });

  return dofCoords;
}


Praetorius, Simon's avatar
Praetorius, Simon committed
234
235
236
template <class Basis, class TP>
  template <class Mat, class Sol, class Rhs>
void PeriodicBC<Basis, TP>::
237
apply(Mat& matrix, Sol& solution, Rhs& rhs)
238
{
239
  periodicBC(matrix, solution, rhs, periodicNodes_, associations_);
240
241
242
}

} // end namespace AMDiS