FirstOrderAssembler.hpp 13.5 KB
Newer Older
1
2
3
4
#pragma once

#include <vector>

5
6
#include <dune/typetree/nodetags.hh>

7
#include <dune/amdis/LocalAssembler.hpp>
8
9
10
11
12
13
#include <dune/amdis/common/Mpl.hpp>
#include <dune/amdis/common/TypeDefs.hpp>
#include <dune/amdis/utility/GetEntity.hpp>

namespace AMDiS
{
14
  template <class GridView, class LocalContext, FirstOrderType type>
15
  class FirstOrderAssembler
16
      : public LocalAssembler<GridView, LocalContext>
17
  {
18
    using Super = LocalAssembler<GridView, LocalContext>;
19
20
21
22
23
24
25

    static constexpr int dim = GridView::dimension;
    static constexpr int dow = GridView::dimensionworld;

    using ElementMatrix = Impl::ElementMatrix;
    using ElementVector = Impl::ElementVector;

26
27
    using OrderTag = std::conditional_t<(type==GRD_PHI), tag::fot_grd_phi, tag::fot_grd_psi>;

28
  public:
29
30
31
32
33
34
    FirstOrderAssembler()
      : Super(1, type)
    {}

    template <class Operator, class... Args>
    void bind(double factor, Operator& op, LocalContext const& localContext, Args&&... args)
35
    {
36
37
38
39
      Super::bind(op, localContext, std::forward<Args>(args)...);

      factor_ = factor;
      op.init(OrderTag{}, localContext, Super::getQuadrature());
40
41
42
    }


43
    // tag dispatching for FirstOrderType...
44
    template <class Operator, class RowNode, class ColNode>
45
    void calculateElementMatrix(Operator& op,
46
47
                                ElementMatrix& elementMatrix,
                                RowNode const& rowNode, ColNode const& colNode)
48
    {
49
      calculateElementMatrix(op, elementMatrix, rowNode, colNode,
50
        typename RowNode::NodeTag{}, typename ColNode::NodeTag{}, FirstOrderType_<type>);
51
52
    }

53
    template <class Operator, class Node>
54
    void calculateElementVector(Operator& op,
55
56
                                ElementVector& elementVector,
                                Node const& node)
57
    {
58
59
      calculateElementVector(op, elementVector, node,
        typename Node::NodeTag{}, FirstOrderType_<type>);
60
61
62
    }


63
    template <class Operator, class RowNode, class ColNode>
64
65
    void calculateElementMatrix(Operator& op,
                                ElementMatrix& elementMatrix,
66
                                RowNode const& rowNode, ColNode const& colNode,
67
                                Dune::TypeTree::LeafNodeTag, Dune::TypeTree::LeafNodeTag,
68
                                FirstOrderType_t<GRD_PHI>)
69
    {
70
71
      auto const& rowLocalFE = rowNode.finiteElement();
      auto const& colLocalFE = colNode.finiteElement();
72
      auto const& quad = Super::getQuadrature().getRule();
73
74
75

      for (std::size_t iq = 0; iq < quad.size(); ++iq) {
        // Position of the current quadrature point in the reference element
76
        decltype(auto) local = Super::getLocal(quad[iq].position());
77
78

        // The transposed inverse Jacobian of the map from the reference element to the element
79
        const auto jacobian = Super::getGeometry().jacobianInverseTransposed(local);
80
81

        // The multiplicative factor in the integral transformation formula
82
        const double factor = Super::getLocalGeometry().integrationElement(quad[iq].position()) * quad[iq].weight() * factor_;
83
84

        std::vector<Dune::FieldVector<double,1> > rowShapeValues;
85
        rowLocalFE.localBasis().evaluateFunction(local, rowShapeValues);
86
87
88

        // The gradients of the shape functions on the reference element
        std::vector<Dune::FieldMatrix<double,1,dim> > colReferenceGradients;
89
        colLocalFE.localBasis().evaluateJacobian(local, colReferenceGradients);
90
91
92
93
94
95
96

        // Compute the shape function gradients on the real element
        std::vector<Dune::FieldVector<double,dow> > colGradients(colReferenceGradients.size());

        for (std::size_t i = 0; i < colGradients.size(); ++i)
          jacobian.mv(colReferenceGradients[i][0], colGradients[i]);

97
98
        for (std::size_t i = 0; i < rowLocalFE.size(); ++i) {
          for (std::size_t j = 0; j < colLocalFE.size(); ++j) {
99
100
            const int local_i = rowNode.localIndex(i);
            const int local_j = colNode.localIndex(j);
101
            op.eval(tag::fot_grd_phi{}, elementMatrix[local_i][local_j], iq, factor, rowShapeValues[i], colGradients[j]);
102
103
104
105
106
          }
        }
      }
    }

107

108
    template <class Operator, class RowNode, class ColNode>
109
110
    void calculateElementMatrix(Operator& op,
                                ElementMatrix& elementMatrix,
111
                                RowNode const& rowNode, ColNode const& colNode,
112
                                Dune::TypeTree::LeafNodeTag, Dune::TypeTree::LeafNodeTag,
113
                                FirstOrderType_t<GRD_PSI>)
114
    {
115
116
      auto const& rowLocalFE = rowNode.finiteElement();
      auto const& colLocalFE = colNode.finiteElement();
117
      auto const& quad = Super::getQuadrature().getRule();
118
119
120

      for (std::size_t iq = 0; iq < quad.size(); ++iq) {
        // Position of the current quadrature point in the reference element
121
        decltype(auto) local = Super::getLocal(quad[iq].position());
122
123

        // The transposed inverse Jacobian of the map from the reference element to the element
124
        const auto jacobian = Super::getGeometry().jacobianInverseTransposed(local);
125
126

        // The multiplicative factor in the integral transformation formula
127
        const double factor = Super::getLocalGeometry().integrationElement(quad[iq].position()) * quad[iq].weight() * factor_;
128
129
130

        // The gradients of the shape functions on the reference element
        std::vector<Dune::FieldMatrix<double,1,dim> > rowReferenceGradients;
131
        rowLocalFE.localBasis().evaluateJacobian(local, rowReferenceGradients);
132
133
134
135
136
137
138
139

        // Compute the shape function gradients on the real element
        std::vector<Dune::FieldVector<double,dow> > rowGradients(rowReferenceGradients.size());

        for (std::size_t i = 0; i < rowGradients.size(); ++i)
          jacobian.mv(rowReferenceGradients[i][0], rowGradients[i]);

        std::vector<Dune::FieldVector<double,1> > colShapeValues;
140
        colLocalFE.localBasis().evaluateFunction(local, colShapeValues);
141

142
143
        for (std::size_t i = 0; i < rowLocalFE.size(); ++i) {
          for (std::size_t j = 0; j < colLocalFE.size(); ++j) {
144
145
            const int local_i = rowNode.localIndex(i);
            const int local_j = colNode.localIndex(j);
146
            op.eval(tag::fot_grd_psi{}, elementMatrix[local_i][local_j], iq, factor, rowGradients[i], colShapeValues[j]);
147
148
149
150
151
          }
        }
      }
    }

152

153
    template <class Operator, class RowNode, class ColNode>
154
155
    void calculateElementMatrix(Operator& op,
                                ElementMatrix& elementMatrix,
156
                                RowNode const& rowNode, ColNode const& colNode,
157
158
159
                                Dune::TypeTree::LeafNodeTag, Dune::TypeTree::PowerNodeTag,
                                FirstOrderType_t<GRD_PHI>)
    {
160
161
      auto const& rowLocalFE = rowNode.finiteElement();
      auto const& colLocalFE = colNode.child(0).finiteElement();
162
163
      auto const& quad = Super::getQuadrature().getRule();

164
      test_exit( dow == degree(colNode), "Number of childs of Power-Node must be DOW");
165
166
167

      for (std::size_t iq = 0; iq < quad.size(); ++iq) {
        // Position of the current quadrature point in the reference element
168
        decltype(auto) local = Super::getLocal(quad[iq].position());
169
170

        // The transposed inverse Jacobian of the map from the reference element to the element
171
        const auto jacobian = Super::getGeometry().jacobianInverseTransposed(local);
172
173

        // The multiplicative factor in the integral transformation formula
174
        const double factor = Super::getLocalGeometry().integrationElement(quad[iq].position()) * quad[iq].weight() * factor_;
175
176

        std::vector<Dune::FieldVector<double,1> > rowShapeValues;
177
        rowLocalFE.localBasis().evaluateFunction(local, rowShapeValues);
178
179
180

        // The gradients of the shape functions on the reference element
        std::vector<Dune::FieldMatrix<double,1,dim> > colReferenceGradients;
181
        colLocalFE.localBasis().evaluateJacobian(local, colReferenceGradients);
182
183
184
185
186
187
188
189
190

        // Compute the shape function gradients on the real element
        std::vector<Dune::FieldVector<double,dow> > colGradients(colReferenceGradients.size());

        for (std::size_t i = 0; i < colGradients.size(); ++i)
          jacobian.mv(colReferenceGradients[i][0], colGradients[i]);

        // <div(u), q>
        for (std::size_t i = 0; i < rowLocalFE.size(); ++i) {
191
          const int local_i = rowNode.localIndex(i);
192
          for (std::size_t j = 0; j < colLocalFE.size(); ++j) {
193
194
            for (std::size_t k = 0; k < degree(colNode); ++k) {
              const int local_j = colNode.child(k).localIndex(j);
195
              op.eval(tag::fot_grd_phi{}, elementMatrix[local_i][local_j], iq, factor, rowShapeValues[i], colGradients[j][k]);
196
197
198
199
200
201
            }
          }
        }
      }
    }

202

203
    template <class Operator, class RowNode, class ColNode>
204
205
    void calculateElementMatrix(Operator& op,
                                ElementMatrix& elementMatrix,
206
                                RowNode const& rowNode, ColNode const& colNode,
207
208
209
                                Dune::TypeTree::PowerNodeTag, Dune::TypeTree::LeafNodeTag,
                                FirstOrderType_t<GRD_PSI>)
    {
210
211
      auto const& rowLocalFE = rowNode.child(0).finiteElement();
      auto const& colLocalFE = colNode.finiteElement();
212
213
      auto const& quad = Super::getQuadrature().getRule();

214
      test_exit( dow == degree(rowNode), "Number of childs of Power-Node must be DOW");
215
216
217

      for (std::size_t iq = 0; iq < quad.size(); ++iq) {
        // Position of the current quadrature point in the reference element
218
        decltype(auto) local = Super::getLocal(quad[iq].position());
219
220

        // The transposed inverse Jacobian of the map from the reference element to the element
221
        const auto jacobian = Super::getGeometry().jacobianInverseTransposed(local);
222
223

        // The multiplicative factor in the integral transformation formula
224
        const double factor = Super::getLocalGeometry().integrationElement(quad[iq].position()) * quad[iq].weight() * factor_;
225
226
227

        // The gradients of the shape functions on the reference element
        std::vector<Dune::FieldMatrix<double,1,dim> > rowReferenceGradients;
228
        rowLocalFE.localBasis().evaluateJacobian(local, rowReferenceGradients);
229
230
231
232
233
234
235
236

        // Compute the shape function gradients on the real element
        std::vector<Dune::FieldVector<double,dow> > rowGradients(rowReferenceGradients.size());

        for (std::size_t i = 0; i < rowGradients.size(); ++i)
          jacobian.mv(rowReferenceGradients[i][0], rowGradients[i]);

        std::vector<Dune::FieldVector<double,1> > colShapeValues;
237
        colLocalFE.localBasis().evaluateFunction(local, colShapeValues);
238
239
240

        // <p, div(v)>
        for (std::size_t j = 0; j < colLocalFE.size(); ++j) {
241
          const int local_j = colNode.localIndex(j);
242
          for (std::size_t i = 0; i < rowLocalFE.size(); ++i) {
243
244
            for (std::size_t k = 0; k < degree(rowNode); ++k) {
              const int local_i = rowNode.child(k).localIndex(i);
245
              op.eval(tag::fot_grd_psi{}, elementMatrix[local_i][local_j], iq, factor, rowGradients[i][k], colShapeValues[j]);
246
247
248
249
250
251
252
            }
          }
        }
      }
    }


253
254
255
    template <class Operator, class RowNode, class ColNode, class RowNodeTag, class ColNodeTag, class FOT>
    void calculateElementMatrix(Operator&, ElementMatrix&,
                                RowNode const&, ColNode const&,
256
                                RowNodeTag, ColNodeTag, FOT)
257
258
259
260
    {
      warning("FirstOrderAssembler::calculateElementMatrix not implemented for RowNode x ColNode");
    }

261

262
    template <class Operator, class Node>
263
264
    void calculateElementVector(Operator& op,
                                ElementVector& elementVector,
265
                                Node const& node,
266
                                Dune::TypeTree::LeafNodeTag,
267
                                FirstOrderType_t<GRD_PSI>)
268
    {
269
      auto const& localFE = node.finiteElement();
270
      auto const& quad = Super::getQuadrature().getRule();
271
272
273

      for (std::size_t iq = 0; iq < quad.size(); ++iq) {
        // Position of the current quadrature point in the reference element
274
        decltype(auto) local = Super::getLocal(quad[iq].position());
275
276

        // The transposed inverse Jacobian of the map from the reference element to the element
277
        const auto jacobian = Super::getGeometry().jacobianInverseTransposed(local);
278
279

        // The multiplicative factor in the integral transformation formula
280
        const double factor = Super::getLocalGeometry().integrationElement(quad[iq].position()) * quad[iq].weight() * factor_;
281
282
283

        // The gradients of the shape functions on the reference element
        std::vector<Dune::FieldMatrix<double,1,dim> > rowReferenceGradients;
284
        localFE.localBasis().evaluateJacobian(local, rowReferenceGradients);
285
286
287
288
289
290
291

        // Compute the shape function gradients on the real element
        std::vector<Dune::FieldVector<double,dow> > rowGradients(rowReferenceGradients.size());

        for (std::size_t i = 0; i < rowGradients.size(); ++i)
          jacobian.mv(rowReferenceGradients[i][0], rowGradients[i]);

292
        for (std::size_t i = 0; i < localFE.size(); ++i) {
293
          const int local_i = node.localIndex(i);
294
          op.eval(tag::fot_grd_psi{}, elementVector[local_i], iq, factor, rowGradients[i]);
295
296
297
        }
      }
    }
298

299

300
301
302
    template <class Operator, class Node, class NodeTag, class FOT>
    void calculateElementVector(Operator&, ElementVector&,
                                Node const&, NodeTag, FOT)
303
    {
304
      warning("FirstOrderAssembler::calculateElementVector not implemented for Node");
305
    }
306
307
308

  private:
    double factor_;
309
310
311
  };

} // end namespace AMDiS