8#include <dune/functions/functionspacebases/subentitydofs.hh>
10#include <amdis/LinearAlgebra.hpp>
11#include <amdis/Output.hpp>
12#include <amdis/typetree/Traversal.hpp>
16template <
class Basis,
class TP>
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))
28 if (!intersection.neighbor()) {
29 hasConnectedIntersections_ =
false;
36 if (*hasConnectedIntersections_)
43template <
class Basis,
class TP>
48 static const auto tol = sqrt(std::numeric_limits<typename Domain::field_type>::epsilon());
49 periodicNodes_.assign(basis_.dimension(),
false);
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())
58 localView.bind(element);
59 for (
const auto& intersection : intersections(gridView, element)) {
60 if (!boundarySubset_(intersection))
63 test_exit(intersection.neighbor(),
64 "Neighbors of periodic intersection not assigned. Use a periodic Grid instead!");
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);
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);
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]);
102 using ctype =
typename D::field_type;
105 CompareTol(ctype tol)
109 bool operator()(D
const& lhs, D
const& rhs)
const
112 for (
int i = 0; i < D::dimension; i++) {
113 if (abs(lhs[i] - rhs[i]) < tol_)
115 return lhs[i] < rhs[i];
126template <
class Basis,
class TP>
127void PeriodicBC<Basis, TP>::
131 static const auto tol = sqrt(std::numeric_limits<typename Domain::field_type>::epsilon());
132 periodicNodes_.assign(basis_.dimension(),
false);
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});
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))
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()});
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)
161 auto el1 = gridView.grid().entity(seeds[0].first);
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);
171 auto el2 = gridView.grid().entity(seeds[1].first);
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);
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]);
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
198 std::vector<Domain> dofCoords(localIndices.size());
199 Traversal::forEachLeafNode(tree, [&](
auto const& node,
auto&&)
201 std::size_t size = node.finiteElement().size();
202 auto geometry = node.element().geometry();
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;
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]);
215 for (std::size_t j = 0; j < localIndices.size(); ++j) {
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)
224 dofCoords[j] = std::move(x);
232template <
class Basis,
class TP>
233 template <
class Mat,
class Sol,
class Rhs>
235apply(Mat& matrix, Sol& solution, Rhs& rhs)
237 periodicBC(matrix, solution, rhs, periodicNodes_, associations_);
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