AMDiS  0.3
The Adaptive Multi-Dimensional Simulation Toolbox
FirstOrderTestDivTrialvec.hpp
1 #pragma once
2 
3 #include <type_traits>
4 
5 #include <amdis/GridFunctionOperator.hpp>
6 #include <amdis/Output.hpp>
7 #include <amdis/common/StaticSize.hpp>
8 
9 namespace AMDiS
10 {
16  namespace tag
17  {
18  struct test_divtrialvec {};
19  }
20 
21 
23  template <class LC, class GridFct>
24  class GridFunctionOperator<tag::test_divtrialvec, LC, GridFct>
25  : public GridFunctionOperatorBase<GridFunctionOperator<tag::test_divtrialvec, LC, GridFct>, LC, GridFct>
26  {
27  using Self = GridFunctionOperator;
29 
30  static_assert( static_size_v<typename GridFct::Range> == 1, "Expression must be of scalar type." );
31 
32  public:
34  : Super(expr, 1)
35  {}
36 
37  template <class CG, class RN, class CN, class Mat>
38  void getElementMatrix(CG const& contextGeo, RN const& rowNode, CN const& colNode, Mat& elementMatrix)
39  {
40  static_assert(RN::isLeaf && CN::isPower,
41  "row-node must be Leaf-Node and col-node must be a Power-Node.");
42 
43  static const std::size_t CHILDREN = CN::CHILDREN;
44  static_assert( CHILDREN == CG::dow, "" );
45 
46  auto const& quad = this->getQuadratureRule(contextGeo.type(), rowNode, colNode);
47  std::size_t rowSize = rowNode.size();
48  std::size_t colSize = colNode.child(0).size();
49 
50  using RangeFieldType = typename CN::ChildType::LocalBasis::Traits::RangeFieldType;
51  using WorldVector = FieldVector<RangeFieldType,CG::dow>;
52  std::vector<WorldVector> colGradients;
53 
54  for (auto const& qp : quad) {
55  // Position of the current quadrature point in the reference element
56  auto&& local = contextGeo.local(qp.position());
57 
58  // The transposed inverse Jacobian of the map from the reference element to the element
59  const auto jacobian = contextGeo.geometry().jacobianInverseTransposed(local);
60 
61  // The multiplicative factor in the integral transformation formula
62  const auto factor = Super::coefficient(local) * contextGeo.integrationElement(qp.position()) * qp.weight();
63 
64  // the values of the shape functions on the reference element at the quadrature point
65  auto const& shapeValues = rowNode.localBasisValuesAt(local);
66 
67  // The gradients of the shape functions on the reference element
68  auto const& shapeGradients = colNode.child(0).localBasisJacobiansAt(local);
69 
70  // Compute the shape function gradients on the real element
71  colGradients.resize(shapeGradients.size());
72 
73  for (std::size_t i = 0; i < colGradients.size(); ++i)
74  jacobian.mv(shapeGradients[i][0], colGradients[i]);
75 
76  for (std::size_t i = 0; i < rowSize; ++i) {
77  const auto local_i = rowNode.localIndex(i);
78  const auto value = factor * shapeValues[i];
79  for (std::size_t j = 0; j < colSize; ++j) {
80  for (std::size_t k = 0; k < CHILDREN; ++k) {
81  const auto local_kj = colNode.child(k).localIndex(j);
82  elementMatrix[local_i][local_kj] += value * colGradients[j][k];
83  }
84  }
85  }
86  }
87 
88  }
89  };
90 
93 } // end namespace AMDiS
The base-template for GridFunctionOperators.
Definition: GridFunctionOperator.hpp:242
Contains all classes needed for solving linear and non linear equation systems.
Definition: AdaptBase.hpp:6
Definition: FirstOrderTestDivTrialvec.hpp:18
The main implementation of an CRTP-base class for operators using a grid-function coefficient to be u...
Definition: GridFunctionOperator.hpp:39