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