OperatorTerm.hpp 11 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
#pragma once

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

#include <dune/geometry/quadraturerules.hh>
#include <dune/istl/bvector.hh>
#include <dune/functions/common/functionconcepts.hh>
#include <dune/functions/functionspacebases/pqknodalbasis.hh>

12
13
#include "OperatorTermBase.hpp"

14
namespace AMDiS 
15
{
16
17
18
19
20
21
22
    
  template <class MeshView>
  class OperatorTerm
  {
  protected:
    using Codim0  = typename MeshView::template Codim<0>;
    using Element = typename Codim0::Entity;
23
24
25
26
    
    static constexpr int dim = Element::dimension;
    
    using QuadratureRule = Dune::QuadratureRule<double, dim>;
27
    using PointList = std::vector<Dune::QuadraturePoint<double, dim>>;
28
29
      
  public:
30
    virtual void init(Element const& element, 
31
32
                      PointList const& points) = 0;
                      
33
    virtual double evalZot(size_t iq, 
34
35
36
        Dune::FieldVector<double,1> const& test, 
        Dune::FieldVector<double,1> const trial = 1.0) const = 0;
        
37
    virtual double evalFot1(size_t iq, 
38
39
40
        Dune::FieldVector<double,1> const& test, 
        Dune::FieldVector<double,dim> const& grad_trial) const = 0;
        
41
    virtual double evalFot2(size_t iq, 
42
43
44
        Dune::FieldVector<double,dim> const& grad_test, 
        Dune::FieldVector<double,1> const trial = 1.0) const = 0;
        
45
    virtual double evalSot(size_t iq, 
46
47
48
        Dune::FieldVector<double,dim> const& grad_test, 
        Dune::FieldVector<double,dim> const& grad_trial) const = 0;
                          
49
    virtual int getDegree() const = 0;
50
  };
51
52
53
  
  
  
54
55
56
57
58
59
  /// Base class for all operators based on expressions
  template <class MeshView, class Term, class Traits = _none>
  class GenericOperatorTerm 
      : public OperatorTerm<MeshView>
      , public OperatorEvaluation
  {
60
61
62
    using Super   = OperatorTerm<MeshView>;
    using Element = typename Super::Element;
    using PointList = typename Super::PointList;
63
    
64
65
    static constexpr int dim = Element::dimension;

66
  public:
67
    GenericOperatorTerm(Term const& term)
68
      : term(term)
69
    {}
70
71

    virtual void init(Element const& element, 
72
                      PointList const& points) override
73
    {
74
75
76
77
78
79
      term.init(element, points);
      
      // cache term evaluation
      values.resize(points.size());
      for (size_t iq = 0; iq < points.size(); ++iq)
        values[iq] = term[iq];
80
    }
81
                      
82
    virtual double evalZot(size_t iq, 
83
84
        Dune::FieldVector<double,1> const& test, 
        Dune::FieldVector<double,1> const trial = 1.0) const override
85
    {
86
      return this->evalZotImpl(_cat{}, _traits{}, values[iq], test, trial);
87
    }
88
        
89
    virtual double evalFot1(size_t iq, 
90
91
        Dune::FieldVector<double,1> const& test, 
        Dune::FieldVector<double,dim> const& grad_trial) const override
92
    {
93
      return this->evalFotImpl(_cat{}, _traits{}, values[iq], grad_trial, test);
94
    }
95
        
96
    virtual double evalFot2(size_t iq, 
97
98
        Dune::FieldVector<double,dim> const& grad_test, 
        Dune::FieldVector<double,1> const trial = 1.0) const override
99
    {
100
      return this->evalFotImpl(_cat{}, _traits{}, values[iq], grad_test, trial);
101
    }
102
        
103
    virtual double evalSot(size_t iq, 
104
105
        Dune::FieldVector<double,dim> const& grad_test, 
        Dune::FieldVector<double,dim> const& grad_trial) const override
106
    {
107
      return this->evalSotImpl(_cat{}, _traits{}, values[iq], grad_test, grad_trial);
108
    }
109
                            
110
111
    virtual int getDegree() const override
    {
112
      return term.getDegree();
113
    }
114

115
  private:
116
    Term term;
117
118
119
120
121
122
    
    using value_type = std::decay_t< decltype( std::declval<Term>()[std::declval<size_t>()] ) >;
    using _cat    = ValueCategory_t<value_type>;
    using _traits = Traits;
    
    std::vector<value_type> values;
123
  };
124
125
  

126
127
128
// some example terms
// -----------------------------------------------------------------------------

129
130
131
132
133
  /// An expression representing a constant (arithmetic) value
  template <class ValueType>
  class ConstantTerm
  {
  public:
134
135
136
    using value_type = ValueType;
    
    ConstantTerm(value_type value)
137
138
139
      : value(value)
    {}
    
140
141
    template <class Element, class PointList>
    void init(Element const& element, 
142
              PointList const& points) { /* do nothing */ }
143
144
    
    value_type operator[](size_t iq) const
145
    {
146
      return value;
147
148
149
150
    }
    
    int getDegree() const
    {
151
      return 0;
152
    }
153
154
      
  private:
155
    value_type value;
156
157
158
159
160
161
162
  };
    
    
  /// generator function for \ref ConstantTerm expressions
  template <class T>
  std::enable_if_t< std::is_arithmetic<T>::value, ConstantTerm<T> >
  constant(T value) { return {value}; }
163
164
  
  
165
// -----------------------------------------------------------------------------
166
167
  
  
168
169
170
171
172
173
  /// An expression that evaluates to the current coordinate of a dof or 
  /// quadrature point with given index.
  template <class Functor>
  class CoordsTerm
  {
  public:
174
175
    using value_type = typename std::result_of<Functor(Dune::FieldVector<double, 2>)>::type;
    
176
177
178
179
180
181
182
    template <class F,
      class = std::enable_if_t<std::is_same<Functor, std::decay_t<F>>::value> >
    CoordsTerm(F&& f, int degree = 1)
      : fct(std::forward<F>(f))
      , degree(degree)
    {}
    
183
184
    template <class Element, class PointList>
    void init(Element const& element, 
185
              PointList const& points)
186
    {
187
188
189
190
      AMDIS_FUNCNAME("CoordsTerm::init()");
      values.resize(points.size());
      for (size_t iq = 0; iq < points.size(); ++iq)
        values[iq] = fct(element.geometry().global(points[iq].position()));
191
192
193
    }
    
    value_type const& operator[](size_t iq) const
194
    {
195
      return values[iq];
196
197
198
199
    }
    
    int getDegree() const
    {
200
      return degree;
201
    }
202
203
      
  private:
204
205
    Functor fct;
    int degree;
206
207
    
    std::vector<value_type> values;
208
  };
209
210
  
  
211
212
213
  /// generator function for \ref CoordsTerm expressions
  template <class F>
  CoordsTerm< std::decay_t<F> > eval(F&& f) { return {std::forward<F>(f)}; }
214
215
  
  
216
// -----------------------------------------------------------------------------
217
218
  
  
219
220
221
222
  // helper class to extract the polynomial degree of a pqk nodal basis
  template <class FeSpace> struct GetDegree : int_<1> {};
  template <class GV, int k, class ST>
  struct GetDegree<Dune::Functions::PQkNodalBasis<GV, k, ST> > : int_<k> {};
223
    
224
225
226
227
228
229
230
  
  /// An expression that evalues a DOFVector at a given element, either on the
  /// dofs or on the quadrature points
  template <class DOFVectorType>
  class DOFVectorTerm
  {      
  public:
231
232
    using value_type = typename DOFVectorType::value_type;
    
233
    DOFVectorTerm(DOFVectorType const& dofvector, double factor = 1.0)
234
      : vector(dofvector.getVector())
235
236
237
238
239
      , factor(factor)
      , localView(dofvector.getFeSpace().localView())
      , localIndexSet(dofvector.getFeSpace().localIndexSet())
    {}
    
240
241
    template <class Element, class PointList>
    void init(Element const& element, 
242
              PointList const& points)
243
    {
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
      AMDIS_FUNCNAME("DOFVectorTerm::init()");
      localView.bind(element);
      localIndexSet.bind(localView);
      
      const auto& localFiniteElem = localView.tree().finiteElement();
      const size_t nBasisFct = localFiniteElem.size();
      
      std::vector<Dune::FieldVector<double,1> > shapeValues(nBasisFct);
      std::vector<value_type> localVec(nBasisFct);
      
      for (size_t j = 0; j < nBasisFct; ++j) {
        const auto global_idx = localIndexSet.index(j);
        localVec[j] = factor * vector[global_idx];
      }
      
      values.resize(points.size());
      for (size_t iq = 0; iq < points.size(); ++iq) {
        localFiniteElem.localBasis().evaluateFunction(points[iq].position(), shapeValues);
        value_type result = 0;
        for (size_t j = 0; j < shapeValues.size(); ++j)
          result += localVec[j] * shapeValues[j];
        values[iq] = result;
      }
267
268
269
270
    }
    
    value_type const& operator[](size_t iq) const
    {
271
      return values[iq];
272
273
274
275
    }
    
    int getDegree() const
    {
276
      return degree;
277
    }
278
279
      
  private:
280
    typename DOFVectorType::BaseVector const& vector;
281
282
    double factor;
    
283
    using Basis = typename DOFVectorType::FeSpace;
284
285
    typename Basis::LocalView localView;
    typename Basis::LocalIndexSet localIndexSet;
286
287
    
    int degree = GetDegree<Basis>::value;
288
289
    
    std::vector<value_type> values;
290
291
  };

292
  
293
294
295
296
297
298
299
  /// generator function for \ref DOFVector expressions
  template <class DOFVectorType>
  DOFVectorTerm<std::decay_t<DOFVectorType>>
  valueOf(DOFVectorType&& vector, double factor = 1.0)
  {
    return {std::forward<DOFVectorType>(vector), factor};
  }
300
  
301
// -----------------------------------------------------------------------------
302
  
303
304
305
306
307
308
309
  
  /// An expression that evalues a DOFVector at a given element, either on the
  /// dofs or on the quadrature points, and applies a functor to the value
  template <class DOFVectorType, class Func>
  class DOFVectorFuncTerm
  {      
  public:
310
311
    using value_type = typename DOFVectorType::value_type;
    
312
313
    template <class F_>
    DOFVectorFuncTerm(DOFVectorType const& dofvector, F_&& f, int f_deg)
314
      : vector(dofvector.getVector())
315
      , f(std::forward<F_>(f))
316
317
318
319
320
321
322
      , localView(dofvector.getFeSpace().localView())
      , localIndexSet(dofvector.getFeSpace().localIndexSet())
      , f_deg(f_deg)
    {}
    
    template <class Element, class PointList>
    void init(Element const& element, 
323
              PointList const& points)
324
    {
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
      AMDIS_FUNCNAME("DOFVectorFuncTerm::init()");
      localView.bind(element);
      localIndexSet.bind(localView);
      
      const auto& localFiniteElem = localView.tree().finiteElement();
      const size_t nBasisFct = localFiniteElem.size();
      
      std::vector<Dune::FieldVector<double,1> > shapeValues(nBasisFct);
      std::vector<value_type> localVec(nBasisFct);
      
      for (size_t j = 0; j < nBasisFct; ++j) {
        const auto global_idx = localIndexSet.index(j);
        localVec[j] = vector[global_idx];
      }
      
      values.resize(points.size());
      for (size_t iq = 0; iq < points.size(); ++iq) {      
        localFiniteElem.localBasis().evaluateFunction(points[iq].position(), shapeValues);
        value_type result = 0;
        for (size_t j = 0; j < shapeValues.size(); ++j)
          result += localVec[j] * shapeValues[j];
        
        values[iq] = f(result);
      }
349
350
351
352
    }
    
    value_type const& operator[](size_t iq) const
    {
353
      return values[iq];
354
355
356
357
    }
    
    int getDegree() const
    {
358
      return degree * f_deg;
359
    }
360
361
      
  private:
362
    typename DOFVectorType::BaseVector const& vector;
363
    Func f;
364
    
365
    using Basis = typename DOFVectorType::FeSpace;
366
367
368
369
370
371
372
    typename Basis::LocalView localView;
    typename Basis::LocalIndexSet localIndexSet;
    
    int degree = GetDegree<Basis>::value;
    int f_deg;
    
    std::vector<value_type> values;
373
  };
374

375
376
377
378
379
380
    
  /// generator function for \ref DOFVectorFuncTerm expressions
  template <class DOFVectorType, class F>
  DOFVectorFuncTerm<std::decay_t<DOFVectorType>, std::decay_t<F>>
  valueOfFunc(DOFVectorType&& vector, F&& f, int deg = 1)
  {
381
    return {std::forward<DOFVectorType>(vector), std::forward<F>(f), deg};
382
  }
383
384
385
  
  
} // end namespace AMDiS