AMDiS  0.3
The Adaptive Multi-Dimensional Simulation Toolbox
SecondOrderDivTestvecDivTrialvec.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  {
19  }
20 
21 
23  template <class LC, class GridFct>
24  class GridFunctionOperator<tag::divtestvec_divtrialvec, LC, GridFct>
25  : public GridFunctionOperatorBase<GridFunctionOperator<tag::divtestvec_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, 2)
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::isPower && CN::isPower,
41  "Operator can be applied to Power-Nodes only.");
42 
43  const bool sameFE = std::is_same_v<FiniteElementType_t<RN>, FiniteElementType_t<CN>>;
44  const bool sameNode = rowNode.treeIndex() == colNode.treeIndex();
45 
46  auto const& quad = this->getQuadratureRule(contextGeo.type(), rowNode, colNode);
47  if (sameFE && sameNode)
48  getElementMatrixOptimized(contextGeo, quad, rowNode, colNode, elementMatrix);
49  else
50  getElementMatrixStandard(contextGeo, quad, rowNode, colNode, elementMatrix);
51  }
52 
53  protected:
54  template <class CG, class QR, class RN, class CN, class Mat>
55  void getElementMatrixStandard(CG const& contextGeo, QR const& quad,
56  RN const& rowNode, CN const& colNode,
57  Mat& elementMatrix)
58  {
59  static const std::size_t CHILDREN = RN::CHILDREN;
60  static_assert( RN::CHILDREN == CN::CHILDREN, "" );
61 
62  std::size_t rowSize = rowNode.child(0).size();
63  std::size_t colSize = colNode.child(0).size();
64 
65  using RowFieldType = typename RN::ChildType::LocalBasis::Traits::RangeFieldType;
66  using RowWorldVector = FieldVector<RowFieldType,CG::dow>;
67  std::vector<RowWorldVector> rowGradients;
68 
69  using ColFieldType = typename CN::ChildType::LocalBasis::Traits::RangeFieldType;
70  using ColWorldVector = FieldVector<ColFieldType,CG::dow>;
71  std::vector<ColWorldVector> colGradients;
72 
73  for (auto const& qp : quad) {
74  // Position of the current quadrature point in the reference element
75  decltype(auto) local = contextGeo.local(qp.position());
76 
77  // The transposed inverse Jacobian of the map from the reference element to the element
78  const auto jacobian = contextGeo.geometry().jacobianInverseTransposed(local);
79 
80  // The multiplicative factor in the integral transformation formula
81  const auto factor = Super::coefficient(local) * contextGeo.integrationElement(qp.position()) * qp.weight();
82 
83  // The gradients of the shape functions on the reference element
84  auto const& rowShapeGradients = rowNode.child(0).localBasisJacobiansAt(local);
85  auto const& colShapeGradients = colNode.child(0).localBasisJacobiansAt(local);
86 
87  // Compute the shape function gradients on the real element
88  rowGradients.resize(rowShapeGradients.size());
89 
90  for (std::size_t i = 0; i < rowGradients.size(); ++i)
91  jacobian.mv(rowShapeGradients[i][0], rowGradients[i]);
92 
93  colGradients.resize(colShapeGradients.size());
94 
95  for (std::size_t i = 0; i < colGradients.size(); ++i)
96  jacobian.mv(colShapeGradients[i][0], colGradients[i]);
97 
98  for (std::size_t i = 0; i < rowSize; ++i) {
99  for (std::size_t j = 0; j < colSize; ++j) {
100  for (std::size_t k = 0; k < CHILDREN; ++k) {
101  const auto local_ki = rowNode.child(k).localIndex(i);
102  const auto local_kj = colNode.child(k).localIndex(j);
103  elementMatrix[local_ki][local_kj] += factor * rowGradients[i][k] * colGradients[j][k];
104  }
105  }
106  }
107  }
108  }
109 
110 
111  template <class CG, class QR, class Node, class Mat>
112  void getElementMatrixOptimized(CG const& contextGeo, QR const& quad,
113  Node const& node, Node const& /*colNode*/,
114  Mat& elementMatrix)
115  {
116  static const std::size_t CHILDREN = Node::CHILDREN;
117 
118  std::size_t size = node.child(0).size();
119 
120  using RangeFieldType = typename Node::ChildType::LocalBasis::Traits::RangeFieldType;
121  using WorldVector = FieldVector<RangeFieldType,CG::dow>;
122  std::vector<WorldVector> gradients;
123 
124  for (auto const& qp : quad) {
125  // Position of the current quadrature point in the reference element
126  auto&& local = contextGeo.local(qp.position());
127 
128  // The transposed inverse Jacobian of the map from the reference element to the element
129  const auto jacobian = contextGeo.geometry().jacobianInverseTransposed(local);
130 
131  // The multiplicative factor in the integral transformation formula
132  const auto factor = Super::coefficient(local) * contextGeo.integrationElement(qp.position()) * qp.weight();
133 
134  // The gradients of the shape functions on the reference element
135  auto const& shapeGradients = node.child(0).localBasisJacobiansAt(local);
136 
137  // Compute the shape function gradients on the real element
138  gradients.resize(shapeGradients.size());
139 
140  for (std::size_t i = 0; i < gradients.size(); ++i)
141  jacobian.mv(shapeGradients[i][0], gradients[i]);
142 
143  for (std::size_t k = 0; k < CHILDREN; ++k) {
144  for (std::size_t i = 0; i < size; ++i) {
145  const auto local_ki = node.child(k).localIndex(i);
146  const auto value = factor * gradients[i][k];
147  elementMatrix[local_ki][local_ki] += value * gradients[i][k];
148 
149  for (std::size_t j = i+1; j < size; ++j) {
150  const auto local_kj = node.child(k).localIndex(j);
151  elementMatrix[local_ki][local_kj] += value * gradients[j][k];
152  elementMatrix[local_kj][local_ki] += value * gradients[j][k];
153  }
154  }
155  }
156  }
157  }
158  };
159 
162 } // 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: SecondOrderDivTestvecDivTrialvec.hpp:18
The main implementation of an CRTP-base class for operators using a grid-function coefficient to be u...
Definition: GridFunctionOperator.hpp:39