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