FirstOrderAssembler.hpp 9.37 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
26
27

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

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

  public:
    template <class Operator>
28
    FirstOrderAssembler(Operator& op, LocalContext const& element)
29
30
31
      : Super(op, element, 1, type)
    {
      if (type == GRD_PHI)
32
        op.initFot1(element, Super::getQuadrature());
33
      else
34
        op.initFot2(element, Super::getQuadrature());
35
36
37
    }


38
    // tag dispatching for FirstOrderType...
39
    template <class Operator, class RowView, class ColView>
40
41
42
    void calculateElementMatrix(Operator& op,
                                RowView const& rowView, ColView const& colView,
                                ElementMatrix& elementMatrix, double fac)
43
    {
44
45
46
47
      using RowNode = typename RowView::Tree;
      using ColNode = typename ColView::Tree;
      calculateElementMatrix(op, rowView, colView, elementMatrix, fac,
        typename RowNode::NodeTag{}, typename ColNode::NodeTag{}, FirstOrderType_<type>);
48
49
50
    }

    template <class Operator, class RowView>
51
52
53
    void calculateElementVector(Operator& op,
                                RowView const& rowView,
                                ElementVector& elementVector, double fac)
54
    {
55
56
57
      using RowNode = typename RowView::Tree;
      calculateElementVector(op, rowView, elementVector, fac,
        typename RowNode::NodeTag{}, FirstOrderType_<type>);
58
59
60
61
62
63
64
65
    }


    template <class Operator, class RowView, class ColView>
    void calculateElementMatrix(Operator& op,
                                RowView const& rowView,
                                ColView const& colView,
                                ElementMatrix& elementMatrix,
66
67
68
                                double fac,
                                Dune::TypeTree::LeafNodeTag,
                                Dune::TypeTree::LeafNodeTag,
69
                                FirstOrderType_t<GRD_PHI>)
70
71
72
73
74
75
    {
      auto geometry = rowView.element().geometry();
      auto const& quad_geometry = Super::getGeometry();

      auto const& rowLocalBasis = rowView.tree().finiteElement().localBasis();
      auto const& colLocalBasis = colView.tree().finiteElement().localBasis();
76
      auto const& quad = Super::getQuadrature().getRule();
77
78
79

      for (std::size_t iq = 0; iq < quad.size(); ++iq) {
        // Position of the current quadrature point in the reference element
80
        decltype(auto) quadPos = Super::map(quad[iq].position());
81
82
83
84
85

        // The transposed inverse Jacobian of the map from the reference element to the element
        const auto jacobian = geometry.jacobianInverseTransposed(quadPos);

        // The multiplicative factor in the integral transformation formula
86
        const double factor = quad_geometry.integrationElement(quad[iq].position()) * quad[iq].weight() * fac;
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

        std::vector<Dune::FieldVector<double,1> > rowShapeValues;
        rowLocalBasis.evaluateFunction(quadPos, rowShapeValues);

        // The gradients of the shape functions on the reference element
        std::vector<Dune::FieldMatrix<double,1,dim> > colReferenceGradients;
        colLocalBasis.evaluateJacobian(quadPos, colReferenceGradients);

        // 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]);

        for (std::size_t i = 0; i < num_rows(elementMatrix); ++i) {
          for (std::size_t j = 0; j < num_cols(elementMatrix); ++j) {
            const int local_i = rowView.tree().localIndex(i);
            const int local_j = colView.tree().localIndex(j);
            op.evalFot1(elementMatrix[local_i][local_j], iq, factor, rowShapeValues[i], colGradients[j]);
          }
        }
      }
    }

    template <class Operator, class RowView, class ColView>
    void calculateElementMatrix(Operator& op,
                                RowView const& rowView,
                                ColView const& colView,
                                ElementMatrix& elementMatrix,
116
117
118
                                double fac,
                                Dune::TypeTree::LeafNodeTag,
                                Dune::TypeTree::LeafNodeTag,
119
                                FirstOrderType_t<GRD_PSI>)
120
121
122
123
124
125
    {
      auto geometry = rowView.element().geometry();
      auto const& quad_geometry = Super::getGeometry();

      auto const& rowLocalBasis = rowView.tree().finiteElement().localBasis();
      auto const& colLocalBasis = colView.tree().finiteElement().localBasis();
126
      auto const& quad = Super::getQuadrature().getRule();
127
128
129

      for (std::size_t iq = 0; iq < quad.size(); ++iq) {
        // Position of the current quadrature point in the reference element
130
        decltype(auto) quadPos = Super::map(quad[iq].position());
131
132
133
134
135

        // The transposed inverse Jacobian of the map from the reference element to the element
        const auto jacobian = geometry.jacobianInverseTransposed(quadPos);

        // The multiplicative factor in the integral transformation formula
136
        const double factor = quad_geometry.integrationElement(quad[iq].position()) * quad[iq].weight() * fac;
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160

        // The gradients of the shape functions on the reference element
        std::vector<Dune::FieldMatrix<double,1,dim> > rowReferenceGradients;
        rowLocalBasis.evaluateJacobian(quadPos, rowReferenceGradients);

        // 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;
        colLocalBasis.evaluateFunction(quadPos, colShapeValues);

        for (std::size_t i = 0; i < num_rows(elementMatrix); ++i) {
          for (std::size_t j = 0; j < num_cols(elementMatrix); ++j) {
            const int local_i = rowView.tree().localIndex(i);
            const int local_j = colView.tree().localIndex(j);
            op.evalFot2(elementMatrix[local_i][local_j], iq, factor, rowGradients[i], colShapeValues[j]);
          }
        }
      }
    }

161
162
163
164
165
166
167
168
    template <class Operator, class RowView, class ColView, class RowNodeTag, class ColNodeTag, class FOT>
    void calculateElementMatrix(Operator& op,
                                RowView const& rowView,
                                ColView const& colView,
                                ElementMatrix& elementMatrix,
                                double fac,
                                RowNodeTag, ColNodeTag, FOT)
    {}
169
170
171
172
173

    template <class Operator, class RowView>
    void calculateElementVector(Operator& op,
                                RowView const& rowView,
                                ElementVector& elementVector,
174
175
                                double fac,
                                Dune::TypeTree::LeafNodeTag,
176
                                FirstOrderType_t<GRD_PSI>)
177
178
179
180
181
    {
      auto geometry = rowView.element().geometry();
      auto const& quad_geometry = Super::getGeometry();

      auto const& rowLocalBasis = rowView.tree().finiteElement().localBasis();
182
      auto const& quad = Super::getQuadrature().getRule();
183
184
185

      for (std::size_t iq = 0; iq < quad.size(); ++iq) {
        // Position of the current quadrature point in the reference element
186
        decltype(auto) quadPos = Super::map(quad[iq].position());
187
188
189
190
191

        // The transposed inverse Jacobian of the map from the reference element to the element
        const auto jacobian = geometry.jacobianInverseTransposed(quadPos);

        // The multiplicative factor in the integral transformation formula
192
        const double factor = quad_geometry.integrationElement(quad[iq].position()) * quad[iq].weight() * fac;
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209

        // The gradients of the shape functions on the reference element
        std::vector<Dune::FieldMatrix<double,1,dim> > rowReferenceGradients;
        rowLocalBasis.evaluateJacobian(quadPos, rowReferenceGradients);

        // 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]);

        for (std::size_t i = 0; i < size(elementVector); ++i) {
          const int local_i = rowView.tree().localIndex(i);
          op.evalFot2(elementVector[local_i], iq, factor, rowGradients[i]);
        }
      }
    }
210
211
212
213
214
215
216
217

    template <class Operator, class RowView, class NodeTag, class FOT>
    void calculateElementVector(Operator& op,
                                RowView const& rowView,
                                ElementVector& elementVector,
                                double fac,
                                NodeTag, FOT)
    {}
218
219
220
  };

} // end namespace AMDiS