LocalToGlobalAdapter.hpp 8.15 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
#pragma once

#include <cstddef>
#include <vector>

#include <dune/common/fmatrix.hh>
#include <dune/common/fvector.hh>
#include <dune/common/typetraits.hh>
#include <dune/geometry/type.hh>

#include <amdis/common/DerivativeTraits.hpp>
#include <amdis/common/FieldMatVec.hpp>
#include <amdis/utility/LocalBasisCache.hpp>

namespace AMDiS
{
  /// Traits class for local-to-global basis adaptors
  /**
   * \tparam LocalBasisTraits Traits class of the LocalBasis to be adapted.
   * \tparam dimGlobal        Dimension of the global coordinates,
   *                          i.e. Geometry::coorddimension, if the global
   *                          coordinates are determined by a Geometry.
   *
   * \implements BasisInterface::Traits
   */
  template <class LocalBasisTraits, std::size_t dimGlobal>
  struct LocalToGlobalBasisAdapterTraits
  {
    static const std::size_t dimDomainLocal = LocalBasisTraits::dimDomain;
    static const std::size_t dimDomainGlobal = dimGlobal;
    static const std::size_t dimRange = LocalBasisTraits::dimRange;

    using DomainField = typename LocalBasisTraits::DomainFieldType;
    using DomainLocal = typename LocalBasisTraits::DomainType;
    using DomainGlobal = FieldVector<DomainField, dimDomainGlobal>;

    using RangeField = typename LocalBasisTraits::RangeFieldType;
    using Range = typename LocalBasisTraits::RangeType;

    using GradientRange = typename DerivativeTraits<Range(DomainGlobal), tag::gradient>::Range;
    using PartialRange = typename DerivativeTraits<Range(DomainGlobal), tag::partial>::Range;
  };

  /// \brief Convert a simple (scalar) local basis into a global basis
  /**
   * The local basis must be scalar, i.e. LocalBasis::Traits::dimRange must be 1
   * It's values are not transformed.
   *
   * For scalar function \f$f\f$, the gradient is equivalent to the transposed
   * Jacobian \f$\nabla f|_x = J_f^T(x)\f$.  The Jacobian is thus transformed
   * using
   * \f[
   *   \nabla f|_{\mu(\hat x)} =
   *       \hat J_\mu^{-T}(\hat x) \cdot \hat\nabla\hat f|_{\hat x}
   * \f]
   * Here the hat \f$\hat{\phantom x}\f$ denotes local quantities and
   * \f$\mu\f$ denotes the local-to-global map of the geometry.
   *
   * \tparam BasisNode  Type of the typetree node containing a local basis to adopt.
   * \tparam Geometry   Type of the local-to-global transformation.
   *
   * NOTE: The adapter implements a caching of local basis evaluations at coordinates.
   */
  template <class BasisNode, class Geometry>
  class LocalToGlobalBasisAdapter
  {
    using LocalBasis = typename BasisNode::FiniteElement::Traits::LocalBasisType;

    static_assert(std::is_same<typename LocalBasis::Traits::DomainFieldType, typename Geometry::ctype>::value,
      "LocalToGlobalBasisAdapter: LocalBasis must use the same ctype as Geometry");

    static_assert(std::size_t(LocalBasis::Traits::dimDomain) == std::size_t(Geometry::mydimension),
      "LocalToGlobalBasisAdapter: LocalBasis domain dimension must match local dimension of Geometry");

  public:
    using Cache = NodeCache<BasisNode>;
    using Traits = LocalToGlobalBasisAdapterTraits<typename LocalBasis::Traits, Geometry::coorddimension>;

79
    /// \brief Construct a LocalToGlobalBasisAdapter
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
    /**
     * \param node     The basis node in the typetree containing the local basis to adopt.
     * \param geometry The geometry object to use for adaption.
     *
     * \note This class stores the references passed here.  Any use of this
     *       class after these references have become invalid results in
     *       undefined behaviour.  The exception is that the destructor of
     *       this class may still be called.
     */
    LocalToGlobalBasisAdapter(BasisNode const& node, Geometry const& geometry)
      : localBasis_(node.finiteElement().localBasis())
      , localBasisCache_(localBasis_)
      , geometry_(geometry)
      , size_(localBasis_.size())
    {}

    /// Return the number of local basis functions
    std::size_t size() const { return size_; }

    /// \brief Return maximum polynomial order of the base function
    /**
     * This is to determine the required quadrature order.  For an affine
     * geometry this is the same order as for the local basis.  For other
     * geometries this returns the order of the local basis plus the global
     * dimension minus 1.  The assumption for non-affine geometries is that
     * they are still multi-linear.
     */
    std::size_t order() const
    {
      if (geometry_.affine())
        // affine linear
        return localBasis_.order();
      else
        // assume at most order dim
        return localBasis_.order() + Traits::dimDomainGlobal - 1;
    }

Praetorius, Simon's avatar
Praetorius, Simon committed
117
    /// Evaluate the local basis functions in the local coordinate `x`
118
119
120
121
122
123
    void evaluateFunction(typename Traits::DomainLocal const& x,
                          std::vector<typename Traits::Range>& out) const
    {
      out = localBasisCache_.evaluateFunction(geometry_.type(), x);
    }

Praetorius, Simon's avatar
Praetorius, Simon committed
124
    /// Evaluate the local basis functions in the local coordinate `x` and
125
126
127
128
129
130
131
    /// return the result using a reference to a thread_local (or static) vector.
    auto const& valuesAt(typename Traits::DomainLocal const& x) const
    {
      return localBasisCache_.evaluateFunction(geometry_.type(), x);
    }

    /// Return the full (global) gradient of the local basis functions in
Praetorius, Simon's avatar
Praetorius, Simon committed
132
    /// the local coordinate `x`
133
134
135
136
137
138
139
140
141
142
143
144
145
    void evaluateGradient(typename Traits::DomainLocal const& x,
                          std::vector<typename Traits::GradientRange>& out) const
    {
      auto const& localJacobian = localBasisCache_.evaluateJacobian(geometry_.type(), x);
      auto const& geoJacobian = geometry_.jacobianInverseTransposed(x);

      out.resize(size_);
      for (std::size_t i = 0; i < size_; ++i)
        geoJacobian.mv(Dune::MatVec::as_vector(localJacobian[i]),
                       Dune::MatVec::as_vector(out[i]));
    }

    /// Return the full (global) gradient of the local basis functions in
Praetorius, Simon's avatar
Praetorius, Simon committed
146
    /// the local coordinate `x` and return the result using a reference
147
148
149
150
151
152
153
154
155
    /// to a thread_local (or static) vector.
    auto const& gradientsAt(typename Traits::DomainLocal const& x) const
    {
      thread_local std::vector<typename Traits::GradientRange> grad;
      evaluateGradient(x, grad);
      return grad;
    }

    /// Return the (global) partial derivative in direction `comp` of the
Praetorius, Simon's avatar
Praetorius, Simon committed
156
    /// local basis functions in the local coordinate `x`
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
    void evaluatePartial(typename Traits::DomainLocal const& x,
                         std::size_t comp,
                         std::vector<typename Traits::PartialRange>& out) const
    {
      auto const& localJacobian = localBasisCache_.evaluateJacobian(geometry_.type(), x);
      // NOTE: geoJacobian might be a Dune::DiagonalMatrix with limited interface!
      auto const& geoJacobian = geometry_.jacobianInverseTransposed(x);

      out.resize(size_);
      typename Traits::GradientRange grad;
      auto&& grad_ = Dune::MatVec::as_vector(grad);
      for (std::size_t i = 0; i < size_; ++i) {
        geoJacobian.mv(Dune::MatVec::as_vector(localJacobian[i]), grad_);
        out[i] = grad_[comp];
      }
    }

    /// Return the (global) partial derivative in direction `comp` of the
Praetorius, Simon's avatar
Praetorius, Simon committed
175
    /// local basis functions in the local coordinate `x` and return the
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
    /// result using a reference to a thread_local (or static) vector.
    auto const& partialsAt(typename Traits::DomainLocal const& x, std::size_t comp) const
    {
      thread_local std::vector<typename Traits::PartialRange> d_comp;
      evaluatePartial(x, comp, d_comp);
      return d_comp;
    }

  private:
    /// Reference to the local basis bound to the adapter
    LocalBasis const& localBasis_;

    /// A basis cache for the evaluation at local coordinates
    Cache localBasisCache_;

    /// Reference to the geometry this adapter is bound to
    Geometry const& geometry_;

    /// The number of basis functions
    std::size_t size_;
  };


  /// Generator function for the \ref LocalToGlobalBasisAdapter. \relates LocalToGlobalBasisAdapter.
  template <class BasisNode, class Geometry>
  LocalToGlobalBasisAdapter<BasisNode, Geometry>
  makeLocalToGlobalBasisAdapter(BasisNode const& node, Geometry const& geometry)
  {
    return LocalToGlobalBasisAdapter<BasisNode, Geometry>{node, geometry};
  }

} // namespace AMDiS