MFEM  v4.1.0 Finite element discretization library
tintrules.hpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2020, Lawrence Livermore National Security, LLC. Produced
2 // at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3 // LICENSE and NOTICE for details. LLNL-CODE-806117.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability visit https://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the BSD-3 license. We welcome feedback and contributions, see file
10 // CONTRIBUTING.md for details.
11
12 #ifndef MFEM_TEMPLATE_INTEGRATION_RULES
13 #define MFEM_TEMPLATE_INTEGRATION_RULES
14
15 #include "../config/tconfig.hpp"
16 #include "../linalg/ttensor.hpp"
17 #include "geom.hpp"
18
19 namespace mfem
20 {
21
22 // Templated integration rules, cf. intrules.?pp
23
24 template <Geometry::Type G, int Q, int Order, typename real_t>
26 {
27 public:
28  static const Geometry::Type geom = G;
30  static const int qpts = Q;
31  static const int order = Order;
32
33  static const bool tensor_prod = false;
34
35  typedef real_t real_type;
36
37 protected:
39
40 public:
42  {
43  const IntegrationRule &ir = GetIntRule();
44  MFEM_ASSERT(ir.GetNPoints() == qpts, "quadrature rule mismatch");
45  for (int j = 0; j < qpts; j++)
46  {
47  weights[j] = ir.IntPoint(j).weight;
48  }
49  }
50
51  // Default copy constructor
52
53  static const IntegrationRule &GetIntRule()
54  {
55  return IntRules.Get(geom, order);
56  }
57
58  // Multi-component weight assignment. qpt_layout_t must be (qpts x n1 x ...)
59  template <AssignOp::Type Op, typename qpt_layout_t, typename qpt_data_t>
60  void AssignWeights(const qpt_layout_t &qpt_layout,
61  qpt_data_t &qpt_data) const
62  {
63  MFEM_STATIC_ASSERT(qpt_layout_t::rank > 1, "invalid rank");
64  MFEM_STATIC_ASSERT(qpt_layout_t::dim_1 == qpts, "invalid size");
65  for (int j = 0; j < qpts; j++)
66  {
67  TAssign<Op>(qpt_layout.ind1(j), qpt_data, weights.data[j]);
68  }
69  }
70
71  template <typename qpt_data_t>
72  void ApplyWeights(qpt_data_t &qpt_data) const
73  {
74  AssignWeights<AssignOp::Mult>(ColumnMajorLayout2D<qpts,1>(), qpt_data);
75  }
76 };
77
78 template <int Dim, int Q, typename real_t>
80
81 template <int Q, typename real_t>
83 {
84 protected:
86
87 public:
88  // Multi-component weight assignment. qpt_layout_t must be (qpts x n1 x ...)
89  template <AssignOp::Type Op, typename qpt_layout_t, typename qpt_data_t>
90  void AssignWeights(const qpt_layout_t &qpt_layout,
91  qpt_data_t &qpt_data) const
92  {
93  MFEM_STATIC_ASSERT(qpt_layout_t::rank > 1, "invalid rank");
94  MFEM_STATIC_ASSERT(qpt_layout_t::dim_1 == Q, "invalid size");
95  for (int j = 0; j < Q; j++)
96  {
97  TAssign<Op>(qpt_layout.ind1(j), qpt_data, weights_1d.data[j]);
98  }
99  }
100
101  template <typename qpt_data_t>
102  void ApplyWeights(qpt_data_t &qpt_data) const
103  {
104  AssignWeights<AssignOp::Mult>(ColumnMajorLayout2D<Q,1>(), qpt_data);
105  }
106 };
107
108 template <int Q, typename real_t>
110 {
111 protected:
113
114 public:
115  // Multi-component weight assignment. qpt_layout_t must be (qpts x n1 x ...)
116  template <AssignOp::Type Op, typename qpt_layout_t, typename qpt_data_t>
117  void AssignWeights(const qpt_layout_t &qpt_layout,
118  qpt_data_t &qpt_data) const
119  {
120  MFEM_STATIC_ASSERT(qpt_layout_t::rank > 1, "invalid rank");
121  MFEM_STATIC_ASSERT(qpt_layout_t::dim_1 == Q*Q, "invalid size");
122  MFEM_FLOPS_ADD(Q*Q);
123  for (int j2 = 0; j2 < Q; j2++)
124  {
125  for (int j1 = 0; j1 < Q; j1++)
126  {
127  TAssign<Op>(
128  qpt_layout.ind1(TMatrix<Q,Q>::layout.ind(j1,j2)), qpt_data,
129  weights_1d.data[j1]*weights_1d.data[j2]);
130  }
131  }
132  }
133
134  template <typename qpt_data_t>
135  void ApplyWeights(qpt_data_t &qpt_data) const
136  {
137  AssignWeights<AssignOp::Mult>(ColumnMajorLayout2D<Q*Q,1>(), qpt_data);
138  }
139 };
140
141 template <int Q, typename real_t>
143 {
144 protected:
146
147 public:
148  // Multi-component weight assignment. qpt_layout_t must be (qpts x n1 x ...)
149  template <AssignOp::Type Op, typename qpt_layout_t, typename qpt_data_t>
150  void AssignWeights(const qpt_layout_t &qpt_layout,
151  qpt_data_t &qpt_data) const
152  {
153  MFEM_STATIC_ASSERT(qpt_layout_t::rank > 1, "invalid rank");
154  MFEM_STATIC_ASSERT(qpt_layout_t::dim_1 == Q*Q*Q, "invalid size");
155  MFEM_FLOPS_ADD(2*Q*Q*Q);
156  for (int j3 = 0; j3 < Q; j3++)
157  {
158  for (int j2 = 0; j2 < Q; j2++)
159  {
160  for (int j1 = 0; j1 < Q; j1++)
161  {
162  TAssign<Op>(
163  qpt_layout.ind1(TTensor3<Q,Q,Q>::layout.ind(j1,j2,j3)),
164  qpt_data,
165  weights_1d.data[j1]*weights_1d.data[j2]*weights_1d.data[j3]);
166  }
167  }
168  }
169  }
170
171  template <typename qpt_data_t>
172  void ApplyWeights(qpt_data_t &qpt_data) const
173  {
174  AssignWeights<AssignOp::Mult>(ColumnMajorLayout2D<Q*Q*Q,1>(), qpt_data);
175  }
176 };
177
178 template <int Dim, int Q, int Order, typename real_t>
180  : public TProductIntegrationRule_base<Dim,Q,real_t>
181 {
182 public:
183  static const Geometry::Type geom =
184  ((Dim == 1) ? Geometry::SEGMENT :
185  ((Dim == 2) ? Geometry::SQUARE : Geometry::CUBE));
186  static const int dim = Dim;
187  static const int qpts_1d = Q;
188  static const int qpts = (Dim == 1) ? Q : ((Dim == 2) ? (Q*Q) : (Q*Q*Q));
189  static const int order = Order;
190
191  static const bool tensor_prod = true;
192
193  typedef real_t real_type;
194
195 protected:
197
198 public:
199  // default constructor, default copy constructor
200 };
201
202 template <int Dim, int Q, typename real_t>
204  : public TProductIntegrationRule<Dim, Q, 2*Q-1, real_t>
205 {
206 public:
207  typedef TProductIntegrationRule<Dim,Q,2*Q-1,real_t> base_class;
208
209  using base_class::geom;
210  using base_class::order;
211  using base_class::qpts_1d;
212
213 protected:
214  using base_class::weights_1d;
215
216 public:
218  {
219  const IntegrationRule &ir_1d = Get1DIntRule();
220  MFEM_ASSERT(ir_1d.GetNPoints() == qpts_1d, "quadrature rule mismatch");
221  for (int j = 0; j < qpts_1d; j++)
222  {
223  weights_1d.data[j] = ir_1d.IntPoint(j).weight;
224  }
225  }
226
228  {
230  }
231  static const IntegrationRule &GetIntRule()
232  {
233  return IntRules.Get(geom, order);
234  }
235 };
236
237 template <Geometry::Type G, int Order, typename real_t = double>
239
240 template <int Order, typename real_t>
241 class TIntegrationRule<Geometry::SEGMENT, Order, real_t>
242  : public GaussIntegrationRule<1, Order/2+1, real_t> { };
243
244 template <int Order, typename real_t>
245 class TIntegrationRule<Geometry::SQUARE, Order, real_t>
246  : public GaussIntegrationRule<2, Order/2+1, real_t> { };
247
248 template <int Order, typename real_t>
249 class TIntegrationRule<Geometry::CUBE, Order, real_t>
250  : public GaussIntegrationRule<3, Order/2+1, real_t> { };
251
252 // Triangle integration rules (based on intrules.cpp)
253 // These specializations define the number of quadrature points for each rule as
254 // a compile-time constant.
255 // TODO: add higher order rules
256 template <typename real_t>
257 class TIntegrationRule<Geometry::TRIANGLE, 0, real_t>
258  : public GenericIntegrationRule<Geometry::TRIANGLE, 1, 0, real_t> { };
259 template <typename real_t>
260 class TIntegrationRule<Geometry::TRIANGLE, 1, real_t>
261  : public GenericIntegrationRule<Geometry::TRIANGLE, 1, 1, real_t> { };
262 template <typename real_t>
263 class TIntegrationRule<Geometry::TRIANGLE, 2, real_t>
264  : public GenericIntegrationRule<Geometry::TRIANGLE, 3, 2, real_t> { };
265 template <typename real_t>
266 class TIntegrationRule<Geometry::TRIANGLE, 3, real_t>
267  : public GenericIntegrationRule<Geometry::TRIANGLE, 4, 3, real_t> { };
268 template <typename real_t>
269 class TIntegrationRule<Geometry::TRIANGLE, 4, real_t>
270  : public GenericIntegrationRule<Geometry::TRIANGLE, 6, 4, real_t> { };
271 template <typename real_t>
272 class TIntegrationRule<Geometry::TRIANGLE, 5, real_t>
273  : public GenericIntegrationRule<Geometry::TRIANGLE, 7, 5, real_t> { };
274 template <typename real_t>
275 class TIntegrationRule<Geometry::TRIANGLE, 6, real_t>
276  : public GenericIntegrationRule<Geometry::TRIANGLE, 12, 6, real_t> { };
277 template <typename real_t>
278 class TIntegrationRule<Geometry::TRIANGLE, 7, real_t>
279  : public GenericIntegrationRule<Geometry::TRIANGLE, 12, 7, real_t> { };
280
281 // Tetrahedron integration rules (based on intrules.cpp)
282 // These specializations define the number of quadrature points for each rule as
283 // a compile-time constant.
284 // TODO: add higher order rules
285 template <typename real_t>
286 class TIntegrationRule<Geometry::TETRAHEDRON, 0, real_t>
287  : public GenericIntegrationRule<Geometry::TETRAHEDRON, 1, 0, real_t> { };
288 template <typename real_t>
289 class TIntegrationRule<Geometry::TETRAHEDRON, 1, real_t>
290  : public GenericIntegrationRule<Geometry::TETRAHEDRON, 1, 1, real_t> { };
291 template <typename real_t>
292 class TIntegrationRule<Geometry::TETRAHEDRON, 2, real_t>
293  : public GenericIntegrationRule<Geometry::TETRAHEDRON, 4, 2, real_t> { };
294 template <typename real_t>
295 class TIntegrationRule<Geometry::TETRAHEDRON, 3, real_t>
296  : public GenericIntegrationRule<Geometry::TETRAHEDRON, 5, 3, real_t> { };
297 template <typename real_t>
298 class TIntegrationRule<Geometry::TETRAHEDRON, 4, real_t>
299  : public GenericIntegrationRule<Geometry::TETRAHEDRON, 11, 4, real_t> { };
300 template <typename real_t>
301 class TIntegrationRule<Geometry::TETRAHEDRON, 5, real_t>
302  : public GenericIntegrationRule<Geometry::TETRAHEDRON, 14, 5, real_t> { };
303 template <typename real_t>
304 class TIntegrationRule<Geometry::TETRAHEDRON, 6, real_t>
305  : public GenericIntegrationRule<Geometry::TETRAHEDRON, 24, 6, real_t> { };
306 template <typename real_t>
307 class TIntegrationRule<Geometry::TETRAHEDRON, 7, real_t>
308  : public GenericIntegrationRule<Geometry::TETRAHEDRON, 31, 7, real_t> { };
309
310 } // namespace mfem
311
312 #endif // MFEM_TEMPLATE_INTEGRATION_RULES
int GetNPoints() const
Returns the number of the points in the integration rule.
Definition: intrules.hpp:245
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:90
void AssignWeights(const qpt_layout_t &qpt_layout, qpt_data_t &qpt_data) const
Definition: tintrules.hpp:60
void AssignWeights(const qpt_layout_t &qpt_layout, qpt_data_t &qpt_data) const
Definition: tintrules.hpp:117
const IntegrationRule & Get(int GeomType, int Order)
Returns an integration rule for given GeomType and Order.
Definition: intrules.cpp:890
void AssignWeights(const qpt_layout_t &qpt_layout, qpt_data_t &qpt_data) const
Definition: tintrules.hpp:150
static const Geometry::Type geom
Definition: tintrules.hpp:183
static const bool tensor_prod
Definition: tintrules.hpp:33
TProductIntegrationRule< Dim, Q, 2 *Q-1, real_t > base_class
Definition: tintrules.hpp:207
static const bool tensor_prod
Definition: tintrules.hpp:191
static const Geometry::Type geom
Definition: tintrules.hpp:28
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
Definition: intrules.hpp:248
void AssignWeights(const qpt_layout_t &qpt_layout, qpt_data_t &qpt_data) const
Definition: tintrules.hpp:90
TVector< qpts, real_t > weights
Definition: tintrules.hpp:38
static const IntegrationRule & GetIntRule()
Definition: tintrules.hpp:53
static const int order
Definition: tintrules.hpp:31
void ApplyWeights(qpt_data_t &qpt_data) const
Definition: tintrules.hpp:102
static const IntegrationRule & Get1DIntRule()
Definition: tintrules.hpp:227
void ApplyWeights(qpt_data_t &qpt_data) const
Definition: tintrules.hpp:72
data_t data[aligned_size >0?aligned_size:1]
Definition: ttensor.hpp:254
IntegrationRules IntRules(0, Quadrature1D::GaussLegendre)
A global object with all integration rules (defined in intrules.cpp)
Definition: intrules.hpp:376
void ApplyWeights(qpt_data_t &qpt_data) const
Definition: tintrules.hpp:135
void ApplyWeights(qpt_data_t &qpt_data) const
Definition: tintrules.hpp:172
static const IntegrationRule & GetIntRule()
Definition: tintrules.hpp:231