Interpolate.hpp 6.13 KB
Newer Older
1
2
3
#pragma once

#include <memory>
4
#include <optional>
5
6
7
8
9
10
11
#include <vector>

#include <dune/common/tuplevector.hh>
#include <dune/common/version.hh>
#include <dune/functions/functionspacebases/interpolate.hh>

#include <amdis/common/Concepts.hpp>
12
#include <amdis/common/FakeContainer.hpp>
13
14
#include <amdis/common/FieldMatVec.hpp>
#include <amdis/common/Logical.hpp>
Praetorius, Simon's avatar
Praetorius, Simon committed
15
#include <amdis/common/Tags.hpp>
16
#include <amdis/functions/FunctionFromCallable.hpp>
17
#include <amdis/functions/HierarchicNodeToRangeMap.hpp>
18
19
#include <amdis/functions/NodeIndices.hpp>
#include <amdis/gridfunctions/GridFunction.hpp>
Praetorius, Simon's avatar
Praetorius, Simon committed
20
#include <amdis/linearalgebra/Traits.hpp>
21
#include <amdis/operations/Assigner.hpp>
22
23
#include <amdis/typetree/Traversal.hpp>

24
25
26
namespace AMDiS
{
  namespace Impl
27
  {
28
29
30
    template <class B, class Vec, class GF, class TP, class C, class BV, class NTRE, class Assign>
    void interpolateTreeSubset(B const& basis, Vec& vector, GF const& gf, TP const& treePath,
                               C& counter, BV const& bitVec, NTRE const& nodeToRangeEntry, Assign assign)
31
    {
32
      auto localView = basis.localView();
33

34
      // set vector to zero at subtree
35
      if (!std::is_same_v<Assign, Assigner::assign>) {
Praetorius, Simon's avatar
Praetorius, Simon committed
36
        for (const auto& e : elements(basis.gridView(), typename BackendTraits<B>::PartitionSet{}))
37
38
39
40
41
42
43
        {
          localView.bind(e);
          auto&& subTree = Dune::TypeTree::child(localView.tree(), treePath);
          vector.forEach(nodeIndices(localView, subTree), [&](auto dof, auto& coeff) {
            if (bitVec[dof]) { coeff = 0; }
          });
        }
44
45
      }

46
47
48
      // Obtain a local view of f
      auto lf = localFunction(gf);

49
50
      vector.init(sizeInfo(basis), false);
      counter.init(sizeInfo(basis), true); // set to zero
Praetorius, Simon's avatar
Praetorius, Simon committed
51
      for (const auto& e : elements(basis.gridView(), typename BackendTraits<B>::PartitionSet{}))
52
      {
53
54
        localView.bind(e);
        lf.bind(e);
55

56
57
58
59
60
        auto&& subTree = Dune::TypeTree::child(localView.tree(),treePath);
        for_each_leaf_node(subTree, [&](auto const& node, auto const& tp)
        {
          using Traits = typename TYPEOF(node)::FiniteElement::Traits::LocalBasisType::Traits;
          using RangeField = typename Traits::RangeFieldType;
61

62
63
          auto&& fe = node.finiteElement();
          std::size_t feSize = fe.localBasis().size();
64

Praetorius, Simon's avatar
Praetorius, Simon committed
65
          auto bitVecRange = mappedRangeView(Dune::range(feSize), [&](std::size_t i) -> bool {
66
67
            return bitVec[localView.index(node.localIndex(i))];
          });
68

69
70
71
          std::vector<RangeField> visit(bitVecRange.begin(), bitVecRange.end());
          if (std::all_of(visit.begin(), visit.end(), [](auto i) { return i == 0; }))
            return;
72

73
          // extract component of local function result corresponding to node in tree
74
          auto localFj = [&](auto const& local)
75
76
77
          {
            const auto& tmp = lf(local);
            return nodeToRangeEntry(node, tp, Dune::MatVec::as_vector(tmp));
78
          };
79

80
          thread_local std::vector<RangeField> interpolationCoeff;
81
82
          auto localFjFct = functionFromCallable<Traits>(localFj);
          fe.localInterpolation().interpolate(localFjFct, interpolationCoeff);
83

84
85
86
          counter.scatter(localView, node, visit, assign);
          vector.scatter(localView, node, interpolationCoeff, visit, assign);
        });
87
      }
88
89
90
      vector.finish();
      counter.finish();
    }
91

92
  } // namespace Impl
93
94


95
96
  template <class T, class Default>
  decltype(auto) value_or(T&& value, Default&&) { return FWD(value); }
97

98
99
  template <class Default>
  decltype(auto) value_or(tag::defaulted, Default&& def) { return FWD(def); }
100
101


102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
  /// \brief Interpolate given function in discrete function space
  /**
   * Interpolation is done wrt the leaf node of the ansatz tree
   * corresponding to the given tree path.
   *
   * Notice that this will only work if the range type of f and
   * the block type of coeff are compatible and supported by
   * flatVectorView.
   *
   * \param basis  Global function space basis of discrete function space
   * \param vec    Coefficient vector to represent the interpolation
   * \param gf     GridFunction to interpolate
   * \param tp_    Tree path specifying the part of the ansatz tree to use [RootTreePath]
   * \param c_     Vector that counts for the number of value assignments [FakeContainer]
   * \param bv_    A vector with flags marking all DOFs that should be interpolated [FakeContainer]
   * \param a_     Assignment mode [Assigner::assign]
   */
  template <class Basis, class Vec, class GF, class TP, class C, class BV, class Assign>
  void interpolate(Basis const& basis, Vec& vec, GF const& gf, TP const& tp_, C&& c_, BV&& bv_, Assign&& a_)
121
  {
122
123
124
125
126
127
128
    auto&& tp = value_or(FWD(tp_), Dune::TypeTree::hybridTreePath());
    auto&& c  = value_or(FWD(c_),  FakeContainer<int,1>());
    auto&& bv = value_or(FWD(bv_), FakeContainer<bool,true>());
    auto&& a  = value_or(FWD(a_),  Assigner::assign());

    auto ntrm = AMDiS::HierarchicNodeToRangeMap();
    AMDiS::Impl::interpolateTreeSubset(basis, vec, gf, tp, c, bv, ntrm, a);
129
130
  }

131
132
  template <class B, class Vec, class GF, class TP, class C, class BV>
  void interpolate(B const& basis, Vec& vec, GF const& gf, TP&& tp, C&& c, BV const& bitVec)
133
  {
134
    static_assert(not std::is_same_v<BV, tag::defaulted>, "");
135
    AMDiS::interpolate(basis, vec, gf, FWD(tp), FWD(c), bitVec, tag::defaulted{});
136
137
  }

138
139
  template <class B, class Vec, class GF, class TP, class C>
  void interpolate(B const& basis, Vec& vec, GF const& gf, TP&& tp, C& counter)
140
  {
141
    static_assert(not std::is_same_v<C, tag::defaulted>, "");
142
    AMDiS::interpolate(basis, vec, gf, FWD(tp), counter, tag::defaulted{}, Assigner::plus_assign{});
143
144
  }

145
146
147
  template <class B, class Vec, class GF, class TreePath>
  void interpolate(B const& basis, Vec& vec, GF const& gf, TreePath const& treePath)
  {
148
    static_assert(not std::is_same_v<TreePath, tag::defaulted>, "");
149
150
    AMDiS::interpolate(basis, vec, gf, treePath, tag::defaulted{}, tag::defaulted{}, Assigner::assign{});
  }
151

152
153
154
155
156
  template <class B, class Vec, class GF>
  void interpolate(B const& basis, Vec& vec, GF const& gf)
  {
    AMDiS::interpolate(basis, vec, gf, tag::defaulted{}, tag::defaulted{}, tag::defaulted{}, Assigner::assign{});
  }
157
158

} // end namespace AMDiS