DOFVectorTerm.hpp 10.2 KB
Newer Older
1
2
3
4
5
6
7
#pragma once

#include <algorithm>
#include <vector>
#include <type_traits>

#include <dune/istl/bvector.hh>
8
#include <dune/functions/gridfunctions/discreteglobalbasisfunction.hh>
9

10
#include <dune/amdis/terms/ExpressionBase.hpp>
11
#include <dune/amdis/utility/GetDegree.hpp>
12

13
namespace AMDiS
14
{
15
  /**
16
17
18
    * \addtogroup OperatorTerms
    * @{
    **/
19

20
21
22
23
  /// \brief An expression that evalues a DOFVector at a given element, either on the
  /// dofs or on the quadrature points. \see valueOf.
  template <class DiscreteGlobalBasisFct>
  class DOFVectorTerm
24
      : ExpressionBase
25
26
27
28
  {
  public:
    using value_type = typename DiscreteGlobalBasisFct::Range;

29
    explicit DOFVectorTerm(std::shared_ptr<DiscreteGlobalBasisFct> const& dofvector, double factor = 1.0)
30
      : dofvector_(dofvector)
31
      , localFunction_(localFunction(*dofvector_))
32
      , factor_(factor)
33
    {}
34

35
    explicit DOFVectorTerm(DiscreteGlobalBasisFct& dofvector, double factor = 1.0)
36
      : localFunction_(localFunction(dofvector))
37
      , factor_(factor)
38
39
40
41
    {}

    template <class Element>
    void bind(Element const& element)
42
    {
43
44
45
      localFunction_.bind(element);

      // TODO: get localView from localFunction
46
47
48
      auto localView = dofvector_->basis().localView();
      localView.bind(element);
      degree_ = getPolynomialDegree(localView.tree());
49
    }
50

51
52
53
54
55
56
    void unbind()
    {
      localFunction_.unbind();
    }

    template <class Context, class PointList>
57
    void init(Context const& context, PointList const& points)
58
59
60
    {
      values_.resize(points.size());
      for (std::size_t iq = 0; iq < points.size(); ++iq)
61
62
63
64
65
66
67
68
69
70
71
72
73
74
        values_[iq] = factor_ * local(context, points[iq].position());
    }

    template <class Context, class LocalCoordinate>
    value_type local(Context const& /*context*/, LocalCoordinate const& lambda) const
    {
      return localFunction_(lambda);
    }

    template <class GlobalCoordinate>
    value_type global(GlobalCoordinate const& x) const
    {
      error_exit("Not implemented.");
      return value_type(0);
75
76
77
78
79
80
81
    }

    value_type const& operator[](std::size_t iq) const
    {
      return values_[iq];
    }

82
83
    // [expects: dofvector_->basis().localView() is bound to an element]
    int getDegree() const
84
85
86
87
88
    {
      return degree_;
    }

  private:
89
    std::shared_ptr<DiscreteGlobalBasisFct> dofvector_ = nullptr;
90
    typename DiscreteGlobalBasisFct::LocalFunction localFunction_;
91
92
    double factor_;

93
    int degree_ = 0;
94
95
    std::vector<value_type> values_;
  };
96

97
  /// Generator function for \ref DOFVectorTerm expressions.
98
  /** Generates a \ref DOFVectorTerm, that takes a \ref DOFVector \p vector,
99
100
101
   *  an optional factor \p factor the value of the DOFVector at the quadrature
   *  points are scaled with.
   **/
102
103
104
  template <class... Args>
  auto valueOf(std::shared_ptr<Dune::Functions::DiscreteGlobalBasisFunction<Args...>> const& dofvector,
               double factor = 1.0)
105
  {
106
107
    using DiscreteGlobalBasisFct = Dune::Functions::DiscreteGlobalBasisFunction<Args...>;
    return DOFVectorTerm<DiscreteGlobalBasisFct>{dofvector, factor};
108
  }
109
110


111
// -----------------------------------------------------------------------------
112

113
#if 0
114
115

  /// \brief An expression that evalues the gradient of a DOFVector at a given
116
  /// element, either on the DOFs or on the quadrature points. \see gradientOf.
117
118
  template <class DOFVectorType>
  class GradientTerm
119
120
121
122
123
  {
    using LocalFE = typename DOFVectorType::FeSpace::LocalView::Tree::FiniteElement;
    using LocalBasis = typename LocalFE::Traits::LocalBasisType;
    using Jacobian = typename LocalBasis::Traits::JacobianType;

124
125
    static constexpr int dim = DOFVectorType::FeSpace::GridView::dimension;
    static constexpr int dow = DOFVectorType::FeSpace::GridView::dimensionworld;
126
127
128
129
130
131
132

    using R = typename LocalBasis::Traits::RangeFieldType;
    using GradientLocalType = Dune::FieldVector<R, dim>;
    using GradientGlobalType = Dune::FieldVector<R, dow>;

  public:
    using value_type = GradientGlobalType;
133
    using field_type = Value_t<DOFVectorType>;
134

135
    explicit GradientTerm(DOFVectorType const& dofvector, double factor = 1.0)
136
137
      : vector(dofvector.getVector())
      , factor(factor)
138
139
      , localView_(dofvector.getFeSpace().localView())
      , localIndexSet_(dofvector.getFeSpace().localIndexSet())
140
    {}
141

142
143
144
145
146
147
148
149
150
151
152
    template <class Element>
    void bind(Element const& element)
    {
      localView_.bind(element);
      localIndexSet_.bind(localView_);
      int deg = getPolynomialDegree(localView_.tree());

      auto t = element.type();
      degree_ = t.isSimplex() ? Math::max(0, deg-1) : deg;
    }

153
    template <class Element, class PointList>
154
    void init(Element const& element,
155
156
157
              PointList const& points)
    {
      AMDIS_FUNCNAME("DOFVectorTerm::init()");
158
159

      auto geometry = element.geometry();
160
      const auto& localFiniteElem = localView_.tree().finiteElement();
161
      const std::size_t nBasisFct = localFiniteElem.size();
162

163
164
      std::vector<Jacobian> shapeGradients(nBasisFct);
      std::vector<field_type> localVec(nBasisFct);
165

166
      for (std::size_t j = 0; j < nBasisFct; ++j) {
167
        const auto global_idx = localIndexSet_.index(j);
168
169
        localVec[j] = factor * vector[global_idx];
      }
170

171
      values.resize(points.size());
172
      for (std::size_t iq = 0; iq < points.size(); ++iq) {
173
174
        const auto jacobian = geometry.jacobianInverseTransposed(points[iq].position());

175
        localFiniteElem.localBasis().evaluateJacobian(points[iq].position(), shapeGradients);
176
        GradientLocalType grdLocal(R(0));
177
        for (std::size_t j = 0; j < nBasisFct; ++j)
178
179
180
          grdLocal.axpy(localVec[j], shapeGradients[j][0]);

        jacobian.mv(grdLocal, values[iq]);
181
182
      }
    }
183

184
    value_type const& operator[](std::size_t iq) const
185
186
187
    {
      return values[iq];
    }
188

189
    int getDegree() const
190
    {
191
      return degree_;
192
    }
193

194
195
196
  private:
    typename DOFVectorType::BaseVector const& vector;
    double factor;
197

198
    using FeSpace = typename DOFVectorType::FeSpace;
199
200
    typename FeSpace::LocalView localView_;
    typename FeSpace::LocalIndexSet localIndexSet_;
201

202
    int degree_ = 0;
203

204
205
206
    std::vector<value_type> values;
  };

207
208
209

  /// \brief Generator function for an expression of gradients of \ref DOFVectors
  /** Generates a \ref GradientTerm, that takes a \ref DOFVector \p vector,
210
211
212
   *  and an optional factor \p factor the value of the gradient of the
   *  DOFVector at the quadrature points are scaled with.
   **/
213
  template <class DOFVectorType>
214
  auto gradientOf(DOFVectorType&& vector, double factor = 1.0)
215
  {
216
    return GradientTerm<std::decay_t<DOFVectorType>>{std::forward<DOFVectorType>(vector), factor};
217
  }
218
219
220


  /// \brief An expression that evalues a partial derivative of a  DOFVector
221
  /// at a given element, either on the DOFs or on the quadrature points. \see derivativeOf.
222
223
  template <class DOFVectorType>
  class DerivativeTerm
224
225
226
227
  {
    using LocalBasis = typename DOFVectorType::FeSpace::LocalView::Tree::FiniteElement::Traits::LocalBasisType;
    using Jacobian = typename LocalBasis::Traits::JacobianType;

228
229
    static constexpr int dim = DOFVectorType::FeSpace::GridView::dimension;
    static constexpr int dow = DOFVectorType::FeSpace::GridView::dimensionworld;
230
231
232
233
234

    using R = typename LocalBasis::Traits::RangeFieldType;
    using GradientLocalType = Dune::FieldVector<R, dim>;
    using Derivative = R;

235
236
  public:
    using value_type = Derivative;
237
    using field_type = Value_t<DOFVectorType>;
238

239
240
241
242
    DerivativeTerm(DOFVectorType const& dofvector, int component, double factor = 1.0)
      : vector(dofvector.getVector())
      , component(component)
      , factor(factor)
243
244
      , localView_(dofvector.getFeSpace().localView())
      , localIndexSet_(dofvector.getFeSpace().localIndexSet())
245
246
    {}

247
248
249
250
251
252
253
254
255
256
257
    template <class Element>
    void bind(Element const& element)
    {
      localView_.bind(element);
      localIndexSet_.bind(localView_);
      int deg = getPolynomialDegree(localView_.tree());

      auto t = element.type();
      degree_ = t.isSimplex() ? Math::max(0, deg-1) : deg;
    }

258
    template <class Element, class PointList>
259
    void init(Element const& element,
260
261
              PointList const& points)
    {
262
      auto geometry = element.geometry();
263
      const auto& localFiniteElem = localView_.tree().finiteElement();
264
      const std::size_t nBasisFct = localFiniteElem.size();
265
266

      std::vector<Jacobian> shapeGradients(nBasisFct);
267
      std::vector<field_type> localVec(nBasisFct);
268

269
      for (std::size_t j = 0; j < nBasisFct; ++j) {
270
        const auto global_idx = localIndexSet_.index(j);
271
272
        localVec[j] = factor * vector[global_idx];
      }
273

274
      values.resize(points.size());
275
      for (std::size_t iq = 0; iq < points.size(); ++iq) {
276
277
278
279
        const auto jacobian = geometry.jacobianInverseTransposed(points[iq].position());
        localFiniteElem.localBasis().evaluateJacobian(points[iq].position(), shapeGradients);

        GradientLocalType grdLocal(R(0));
280
        for (std::size_t j = 0; j < nBasisFct; ++j)
281
282
283
          grdLocal.axpy(localVec[j], shapeGradients[j][0]);

        values[iq] = 0;
284
        for (std::size_t j = 0; j < dim; ++j)
285
          values[iq] += jacobian[component][j] * grdLocal[j];
286
287
      }
    }
288

289
    value_type const& operator[](std::size_t iq) const
290
291
292
    {
      return values[iq];
    }
293

294
    int getDegree() const
295
    {
296
      return degree_;
297
    }
298

299
300
301
302
  private:
    typename DOFVectorType::BaseVector const& vector;
    int component;
    double factor;
303

304
    using FeSpace = typename DOFVectorType::FeSpace;
305
306
    typename FeSpace::LocalView localView_;
    typename FeSpace::LocalIndexSet localIndexSet_;
307

308
    int degree_ = 0;
309

310
311
312
    std::vector<value_type> values;
  };

313
314
315

  /// \brief Generator function for an expression of partial derivatives of \ref DOFVectors
  /** Generates a \ref DerivativeTerm, that takes a \ref DOFVector \p vector,
316
317
318
319
   *  a index representing the component of the gradient that should be returned,
   *  and an optional factor \p factor the value of the partial derivative of the
   *  DOFVector at the quadrature points are scaled with.
   **/
320
321
322
323
324
325
  template <class DOFVectorType>
  DerivativeTerm<std::decay_t<DOFVectorType>>
  derivativeOf(DOFVectorType&& vector, int component, double factor = 1.0)
  {
    return {std::forward<DOFVectorType>(vector), component, factor};
  }
326

327
328
  #endif

329
  /** @} **/
330

331
} // end namespace AMDiS