AMDiS  0.1
The Adaptive Multi-Dimensional Simulation Toolbox
SecondOrderGradTestGradTrial.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 #include <amdis/common/ValueCategory.hpp>
9 
10 namespace AMDiS
11 {
17  namespace tag
18  {
19  struct gradtest_gradtrial {};
20  }
21 
22 
24  template <class LC, class GridFct>
25  class GridFunctionOperator<tag::gradtest_gradtrial, LC, GridFct>
26  : public GridFunctionOperatorBase<GridFunctionOperator<tag::gradtest_gradtrial, LC, GridFct>, LC, GridFct>
27  {
28  static const int dow = ContextGeometry<LC>::dow;
29  using Self = GridFunctionOperator;
31 
32  using expr_value_type = typename GridFct::Range;
33  static_assert( static_size_v<expr_value_type> == 1 || (static_num_rows_v<expr_value_type> == dow && static_num_cols_v<expr_value_type> == dow),
34  "Expression must be of scalar or matrix type." );
35 
36  public:
38  : Super(expr, 2)
39  {}
40 
41  template <class CG, class RN, class CN, class Mat>
42  void getElementMatrix(CG const& contextGeo, RN const& rowNode, CN const& colNode, Mat& elementMatrix)
43  {
44  static_assert(RN::isLeaf && CN::isLeaf,
45  "Operator can be applied to Leaf-Nodes only.");
46 
47  const bool sameFE = std::is_same_v<FiniteElementType_t<RN>, FiniteElementType_t<CN>>;
48  const bool sameNode = rowNode.treeIndex() == colNode.treeIndex();
49  using Category = ValueCategory_t<typename GridFct::Range>;
50 
51  auto const& quad = this->getQuadratureRule(contextGeo.type(), rowNode, colNode);
52  if (sameFE && sameNode)
53  getElementMatrixOptimized(contextGeo, quad, rowNode, colNode, elementMatrix, Category{});
54  else if (sameFE)
55  getElementMatrixStandard(contextGeo, quad, rowNode, colNode, elementMatrix);
56  else
57  error_exit("Not implemented: currently only the implementation for equal fespaces available");
58  }
59 
60  protected:
61 
62  template <class CG, class QR, class RN, class CN, class Mat>
63  void getElementMatrixStandard(CG const& contextGeo, QR const& quad,
64  RN const& rowNode, CN const& colNode,
65  Mat& elementMatrix)
66  {
67  std::size_t size = rowNode.size();
68 
69  using RangeFieldType = typename RN::LocalBasis::Traits::RangeFieldType;
70  using WorldVector = Dune::FieldVector<RangeFieldType,CG::dow>;
71  std::vector<WorldVector> gradients;
72 
73  for (auto const& qp : quad) {
74  // Position of the current quadrature point in the reference element
75  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 = contextGeo.integrationElement(qp.position()) * qp.weight();
82  const auto exprValue = Super::coefficient(local);
83 
84  // The gradients of the shape functions on the reference element
85  auto const& shapeGradients = rowNode.localBasisJacobiansAt(local);
86 
87  // Compute the shape function gradients on the real element
88  gradients.resize(shapeGradients.size());
89 
90  for (std::size_t i = 0; i < gradients.size(); ++i)
91  jacobian.mv(shapeGradients[i][0], gradients[i]);
92 
93  for (std::size_t i = 0; i < size; ++i) {
94  const auto local_i = rowNode.localIndex(i);
95  for (std::size_t j = 0; j < size; ++j) {
96  const auto local_j = colNode.localIndex(j);
97  elementMatrix[local_i][local_j] += eval(exprValue, factor, gradients[i], gradients[j]);
98  }
99  }
100  }
101  }
102 
103  template <class CG, class QR, class RN, class CN, class Mat>
104  void getElementMatrixOptimized(CG const& contextGeo, QR const& quad,
105  RN const& node, CN const& /*colNode*/,
106  Mat& elementMatrix, tag::scalar)
107  {
108  std::size_t size = node.size();
109 
110  using RangeFieldType = typename RN::LocalBasis::Traits::RangeFieldType;
111  using WorldVector = Dune::FieldVector<RangeFieldType,CG::dow>;
112  std::vector<WorldVector> gradients;
113 
114  for (auto const& qp : quad) {
115  // Position of the current quadrature point in the reference element
116  auto&& local = contextGeo.local(qp.position());
117 
118  // The transposed inverse Jacobian of the map from the reference element to the element
119  const auto jacobian = contextGeo.geometry().jacobianInverseTransposed(local);
120 
121  // The multiplicative factor in the integral transformation formula
122  const auto factor = Super::coefficient(local) * contextGeo.integrationElement(qp.position())
123  * qp.weight();
124 
125  // The gradients of the shape functions on the reference element
126  auto const& shapeGradients = node.localBasisJacobiansAt(local);
127 
128  // Compute the shape function gradients on the real element
129  gradients.resize(shapeGradients.size());
130  for (std::size_t i = 0; i < gradients.size(); ++i)
131  jacobian.mv(shapeGradients[i][0], gradients[i]);
132 
133  for (std::size_t i = 0; i < size; ++i) {
134  const auto local_i = node.localIndex(i);
135 
136  elementMatrix[local_i][local_i] += factor * (gradients[i] * gradients[i]);
137 
138  for (std::size_t j = i+1; j < size; ++j) {
139  const auto local_j = node.localIndex(j);
140  const auto value = factor * (gradients[i] * gradients[j]);
141 
142  elementMatrix[local_i][local_j] += value;
143  elementMatrix[local_j][local_i] += value;
144  }
145  }
146  }
147  }
148 
149  template <class CG, class QR, class RN, class CN, class Mat>
150  void getElementMatrixOptimized(CG const& contextGeo, QR const& quad,
151  RN const& node, CN const& /*colNode*/,
152  Mat& elementMatrix, tag::matrix)
153  {
154  std::size_t size = node.size();
155 
156  using RangeFieldType = typename RN::LocalBasis::Traits::RangeFieldType;
157  using WorldVector = Dune::FieldVector<RangeFieldType,CG::dow>;
158  std::vector<WorldVector> gradients;
159 
160  for (auto const& qp : quad) {
161  // Position of the current quadrature point in the reference element
162  auto&& local = contextGeo.local(qp.position());
163 
164  // The transposed inverse Jacobian of the map from the reference element to the element
165  const auto jacobian = contextGeo.geometry().jacobianInverseTransposed(local);
166 
167  // The multiplicative factor in the integral transformation formula
168  const auto factor = contextGeo.integrationElement(qp.position()) * qp.weight();
169  const auto exprValue = Super::coefficient(local);
170 
171  // The gradients of the shape functions on the reference element
172  auto const& shapeGradients = node.localBasisJacobiansAt(local);
173 
174  // Compute the shape function gradients on the real element
175  gradients.resize(shapeGradients.size());
176  for (std::size_t i = 0; i < gradients.size(); ++i)
177  jacobian.mv(shapeGradients[i][0], gradients[i]);
178 
179  for (std::size_t i = 0; i < size; ++i) {
180  const auto local_i = node.localIndex(i);
181  for (std::size_t j = 0; j < size; ++j) {
182  const auto local_j = node.localIndex(j);
183  elementMatrix[local_i][local_j] += eval(exprValue, factor, gradients[i], gradients[j]);
184  }
185  }
186  }
187  }
188 
189  protected:
190 
191  template <class S, class F, class T, int dow,
192  std::enable_if_t<Category::Scalar<S>,int> = 0>
193  T eval(S const& scalar, F factor,
194  Dune::FieldVector<T,dow> const& grad_test,
195  Dune::FieldVector<T,dow> const& grad_trial) const
196  {
197  return (scalar * factor) * (grad_test * grad_trial);
198  }
199 
200  template <class M, class F, class T, int dow,
201  std::enable_if_t<Category::Matrix<M>,int> = 0>
202  T eval(M const& mat, F factor,
203  Dune::FieldVector<T,dow> const& grad_test,
204  Dune::FieldVector<T,dow> const& grad_trial) const
205  {
206  return factor * (grad_test * (mat * grad_trial));
207  }
208  };
209 
211  template <class Expr, class... QuadratureArgs>
212  auto sot(Expr&& expr, QuadratureArgs&&... args)
213  {
214  return makeOperator(tag::gradtest_gradtrial{}, FWD(expr), FWD(args)...);
215  }
216 
219 } // end namespace AMDiS
Definition: Tags.hpp:11
Definition: SecondOrderGradTestGradTrial.hpp:19
auto makeOperator(Tag tag, Expr &&expr, QuadratureArgs &&... args)
Store tag and expression into a PreGridFunctionOperator to create a GridFunctionOperator.
Definition: GridFunctionOperator.hpp:220
The base-template for GridFunctionOperators.
Definition: GridFunctionOperator.hpp:242
void error_exit(std::string const &str, Args &&... args)
print a message and exit
Definition: Output.hpp:142
Contains all classes needed for solving linear and non linear equation systems.
Definition: AdaptBase.hpp:6
Definition: Tags.hpp:9
auto sot(Expr &&expr, QuadratureArgs &&... args)
Create a second-order term.
Definition: SecondOrderGradTestGradTrial.hpp:212
Wrapper class for element and geometry.
Definition: ContextGeometry.hpp:43
The main implementation of an CRTP-base class for operators using a grid-function coefficient to be u...
Definition: GridFunctionOperator.hpp:39