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

#include <vector>

#include <dune/amdis/Assembler.hpp>
#include <dune/amdis/common/Mpl.hpp>
#include <dune/amdis/common/TypeDefs.hpp>
#include <dune/amdis/utility/GetEntity.hpp>

namespace AMDiS
{
12
  template <class GridView, class LocalContext, FirstOrderType type>
13
  class FirstOrderAssembler
14
      : public Assembler<GridView, LocalContext>
15
  {
16
    using Super = Assembler<GridView, LocalContext>;
17
18
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;

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


36
    // tag dispatching for FirstOrderType...
37
    template <class Operator, class RowView, class ColView>
38
    void calculateElementMatrix(Operator& op, RowView const& rowView, ColView const& colView, ElementMatrix& elementMatrix)
39
    {
40
      calculateElementMatrix(op, rowView, colView, elementMatrix, FirstOrderType_<type>);
41
42
43
44
45
    }

    template <class Operator, class RowView>
    void calculateElementVector(Operator& op, RowView const& rowView, ElementVector& elementVector)
    {
46
      calculateElementVector(op, rowView, elementVector, FirstOrderType_<type>);
47
48
49
50
51
52
53
54
    }


    template <class Operator, class RowView, class ColView>
    void calculateElementMatrix(Operator& op,
                                RowView const& rowView,
                                ColView const& colView,
                                ElementMatrix& elementMatrix,
55
                                FirstOrderType_t<GRD_PHI>)
56
57
58
59
60
61
    {
      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();
62
      auto const& quad = Super::getQuadrature().getRule();
63
64
65

      for (std::size_t iq = 0; iq < quad.size(); ++iq) {
        // Position of the current quadrature point in the reference element
66
        decltype(auto) quadPos = Super::map(quad[iq].position());
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102

        // 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
        const double factor = quad_geometry.integrationElement(quad[iq].position()) * quad[iq].weight();

        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,
103
                                FirstOrderType_t<GRD_PSI>)
104
105
106
107
108
109
    {
      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();
110
      auto const& quad = Super::getQuadrature().getRule();
111
112
113

      for (std::size_t iq = 0; iq < quad.size(); ++iq) {
        // Position of the current quadrature point in the reference element
114
        decltype(auto) quadPos = Super::map(quad[iq].position());
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149

        // 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
        const double factor = quad_geometry.integrationElement(quad[iq].position()) * quad[iq].weight();

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


    template <class Operator, class RowView>
    void calculateElementVector(Operator& op,
                                RowView const& rowView,
                                ElementVector& elementVector,
150
                                FirstOrderType_t<GRD_PSI>)
151
152
153
154
155
    {
      auto geometry = rowView.element().geometry();
      auto const& quad_geometry = Super::getGeometry();

      auto const& rowLocalBasis = rowView.tree().finiteElement().localBasis();
156
      auto const& quad = Super::getQuadrature().getRule();
157
158
159

      for (std::size_t iq = 0; iq < quad.size(); ++iq) {
        // Position of the current quadrature point in the reference element
160
        decltype(auto) quadPos = Super::map(quad[iq].position());
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186

        // 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
        const double factor = quad_geometry.integrationElement(quad[iq].position()) * quad[iq].weight();

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

} // end namespace AMDiS