MFEM  v4.6.0
Finite element discretization library
tevaluator.hpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2023, 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_EVALUATOR
13 #define MFEM_TEMPLATE_EVALUATOR
14 
15 #include "../config/tconfig.hpp"
16 #include "../linalg/ttensor.hpp"
17 #include "../general/error.hpp"
18 #include "fespace.hpp"
19 
20 namespace mfem
21 {
22 
23 // Templated classes for transitioning between degrees of freedom and quadrature
24 // points values.
25 
26 /** @brief Shape evaluators -- values of basis functions on the reference element
27  @tparam FE some form of TFiniteElement, probably got from TMesh::FE_type
28  @tparam IR some form of TIntegrationRule
29  @tparam TP tensor product or not
30  @tparam real_t data type for mesh nodes, solution basis, mesh basis
31 */
32 template <class FE, class IR, bool TP, typename real_t>
34 
35 /// ShapeEvaluator without tensor-product structure
36 template <class FE, class IR, typename real_t>
37 class ShapeEvaluator_base<FE, IR, false, real_t>
38 {
39 public:
40  static const int DOF = FE::dofs;
41  static const int NIP = IR::qpts;
42  static const int DIM = FE::dim;
43 
44 protected:
49 
50 public:
51  ShapeEvaluator_base(const FE &fe)
52  {
53  fe.CalcShapes(IR::GetIntRule(), B.data, G.data);
54  TAssign<AssignOp::Set>(Bt.layout, Bt, B.layout.transpose_12(), B);
55  TAssign<AssignOp::Set>(Gt.layout.merge_23(), Gt,
56  G.layout.merge_12().transpose_12(), G);
57  }
58 
59  // default copy constructor
60 
61  /** @brief Multi-component shape evaluation from DOFs to quadrature points.
62  dof_layout is (DOF x NumComp) and qpt_layout is (NIP x NumComp). */
63  template <typename dof_layout_t, typename dof_data_t,
64  typename qpt_layout_t, typename qpt_data_t>
65  inline MFEM_ALWAYS_INLINE
66  void Calc(const dof_layout_t &dof_layout, const dof_data_t &dof_data,
67  const qpt_layout_t &qpt_layout, qpt_data_t &qpt_data) const
68  {
69  MFEM_STATIC_ASSERT(dof_layout_t::rank == 2 &&
70  dof_layout_t::dim_1 == DOF,
71  "invalid dof_layout_t.");
72  MFEM_STATIC_ASSERT(qpt_layout_t::rank == 2 &&
73  qpt_layout_t::dim_1 == NIP,
74  "invalid qpt_layout_t.");
75  MFEM_STATIC_ASSERT(dof_layout_t::dim_2 == qpt_layout_t::dim_2,
76  "incompatible dof- and qpt- layouts.");
77 
78  Mult_AB<false>(B.layout, B,
79  dof_layout, dof_data,
80  qpt_layout, qpt_data);
81  }
82 
83  /** @brief Multi-component shape evaluation transpose from quadrature points to
84  DOFs. qpt_layout is (NIP x NumComp) and dof_layout is (DOF x NumComp). */
85  template <bool Add,
86  typename qpt_layout_t, typename qpt_data_t,
87  typename dof_layout_t, typename dof_data_t>
88  inline MFEM_ALWAYS_INLINE
89  void CalcT(const qpt_layout_t &qpt_layout, const qpt_data_t &qpt_data,
90  const dof_layout_t &dof_layout, dof_data_t &dof_data) const
91  {
92  MFEM_STATIC_ASSERT(dof_layout_t::rank == 2 &&
93  dof_layout_t::dim_1 == DOF,
94  "invalid dof_layout_t.");
95  MFEM_STATIC_ASSERT(qpt_layout_t::rank == 2 &&
96  qpt_layout_t::dim_1 == NIP,
97  "invalid qpt_layout_t.");
98  MFEM_STATIC_ASSERT(dof_layout_t::dim_2 == qpt_layout_t::dim_2,
99  "incompatible dof- and qpt- layouts.");
100 
101  Mult_AB<Add>(Bt.layout, Bt,
102  qpt_layout, qpt_data,
103  dof_layout, dof_data);
104  }
105 
106  /** @brief Multi-component gradient evaluation from DOFs to quadrature points.
107  dof_layout is (DOF x NumComp) and grad_layout is (NIP x DIM x NumComp). */
108  template <typename dof_layout_t, typename dof_data_t,
109  typename grad_layout_t, typename grad_data_t>
110  inline MFEM_ALWAYS_INLINE
111  void CalcGrad(const dof_layout_t &dof_layout,
112  const dof_data_t &dof_data,
113  const grad_layout_t &grad_layout,
114  grad_data_t &grad_data) const
115  {
116  MFEM_STATIC_ASSERT(dof_layout_t::rank == 2 &&
117  dof_layout_t::dim_1 == DOF,
118  "invalid dof_layout_t.");
119  MFEM_STATIC_ASSERT(grad_layout_t::rank == 3 &&
120  grad_layout_t::dim_1 == NIP &&
121  grad_layout_t::dim_2 == DIM,
122  "invalid grad_layout_t.");
123  MFEM_STATIC_ASSERT(dof_layout_t::dim_2 == grad_layout_t::dim_3,
124  "incompatible dof- and grad- layouts.");
125 
126  Mult_AB<false>(G.layout.merge_12(), G,
127  dof_layout, dof_data,
128  grad_layout.merge_12(), grad_data);
129  }
130 
131  /** @brief Multi-component gradient evaluation transpose from quadrature points to
132  DOFs. grad_layout is (NIP x DIM x NumComp), dof_layout is (DOF x NumComp). */
133  template <bool Add,
134  typename grad_layout_t, typename grad_data_t,
135  typename dof_layout_t, typename dof_data_t>
136  inline MFEM_ALWAYS_INLINE
137  void CalcGradT(const grad_layout_t &grad_layout,
138  const grad_data_t &grad_data,
139  const dof_layout_t &dof_layout,
140  dof_data_t &dof_data) const
141  {
142  MFEM_STATIC_ASSERT(dof_layout_t::rank == 2 &&
143  dof_layout_t::dim_1 == DOF,
144  "invalid dof_layout_t.");
145  MFEM_STATIC_ASSERT(grad_layout_t::rank == 3 &&
146  grad_layout_t::dim_1 == NIP &&
147  grad_layout_t::dim_2 == DIM,
148  "invalid grad_layout_t.");
149  MFEM_STATIC_ASSERT(dof_layout_t::dim_2 == grad_layout_t::dim_3,
150  "incompatible dof- and grad- layouts.");
151 
152  Mult_AB<Add>(Gt.layout.merge_23(), Gt,
153  grad_layout.merge_12(), grad_data,
154  dof_layout, dof_data);
155  }
156 
157  /** @brief Multi-component assemble.
158  qpt_layout is (NIP x NumComp),
159  M_layout is (DOF x DOF x NumComp) */
160  template <typename qpt_layout_t, typename qpt_data_t,
161  typename M_layout_t, typename M_data_t>
162  inline MFEM_ALWAYS_INLINE
163  void Assemble(const qpt_layout_t &qpt_layout, const qpt_data_t &qpt_data,
164  const M_layout_t &M_layout, M_data_t &M_data) const
165  {
166  // M_{i,j,k} = \sum_{s} B_{s,i} B_{s,j} qpt_data_{s,k}
167  // Using TensorAssemble: <1,NIP,NC> --> <DOF,1,DOF,NC>
168 #if 0
169  TensorAssemble<false>(
170  B.layout, B,
171  qpt_layout.template split_1<1,NIP>(), qpt_data,
172  M_layout.template split_1<DOF,1>(), M_data);
173 #else
174  TensorAssemble<false>(
175  Bt.layout, Bt, B.layout, B,
176  qpt_layout.template split_1<1,NIP>(), qpt_data,
177  M_layout.template split_1<DOF,1>(), M_data);
178 #endif
179  }
180 
181  /** @brief Multi-component assemble of grad-grad element matrices.
182  qpt_layout is (NIP x DIM x DIM x NumComp), and
183  D_layout is (DOF x DOF x NumComp). */
184  template <typename qpt_layout_t, typename qpt_data_t,
185  typename D_layout_t, typename D_data_t>
186  inline MFEM_ALWAYS_INLINE
187  void AssembleGradGrad(const qpt_layout_t &qpt_layout,
188  const qpt_data_t &qpt_data,
189  const D_layout_t &D_layout,
190  D_data_t &D_data) const
191  {
192  const int NC = qpt_layout_t::dim_4;
193  typedef typename qpt_data_t::data_type entry_type;
195  for (int k = 0; k < NC; k++)
196  {
197  // Next loop performs a batch of matrix-matrix products of size
198  // (DIM x DIM) x (DIM x DOF) --> (DIM x DOF)
199  for (int j = 0; j < NIP; j++)
200  {
201  Mult_AB<false>(qpt_layout.ind14(j,k), qpt_data,
202  G.layout.ind1(j), G,
203  F.layout.ind14(j,k), F);
204  }
205  }
206  // (DOF x (NIP x DIM)) x ((NIP x DIM) x DOF x NC) --> (DOF x DOF x NC)
207  Mult_2_1<false>(Gt.layout.merge_23(), Gt,
208  F.layout.merge_12(), F,
209  D_layout, D_data);
210  }
211 };
212 
213 template <int Dim, int DOF, int NIP, typename real_t>
215 
216 /// ShapeEvaluator with 1D tensor-product structure
217 template <int DOF, int NIP, typename real_t>
218 class TProductShapeEvaluator<1, DOF, NIP, real_t>
219 {
220 protected:
221  static const int TDOF = DOF; // total dofs
222 
225 
226 public:
228 
229  /** @brief Multi-component shape evaluation from DOFs to quadrature points.
230  dof_layout is (DOF x NumComp) and qpt_layout is (NIP x NumComp). */
231  template <typename dof_layout_t, typename dof_data_t,
232  typename qpt_layout_t, typename qpt_data_t>
233  inline MFEM_ALWAYS_INLINE
234  void Calc(const dof_layout_t &dof_layout, const dof_data_t &dof_data,
235  const qpt_layout_t &qpt_layout, qpt_data_t &qpt_data) const
236  {
237  Mult_AB<false>(B_1d.layout, B_1d,
238  dof_layout, dof_data,
239  qpt_layout, qpt_data);
240  }
241 
242  /** @brief Multi-component shape evaluation transpose from quadrature points
243  to DOFs. qpt_layout is (NIP x NumComp) and dof_layout is (DOF x NumComp). */
244  template <bool Add,
245  typename qpt_layout_t, typename qpt_data_t,
246  typename dof_layout_t, typename dof_data_t>
247  inline MFEM_ALWAYS_INLINE
248  void CalcT(const qpt_layout_t &qpt_layout, const qpt_data_t &qpt_data,
249  const dof_layout_t &dof_layout, dof_data_t &dof_data) const
250  {
251  Mult_AB<Add>(Bt_1d.layout, Bt_1d,
252  qpt_layout, qpt_data,
253  dof_layout, dof_data);
254  }
255 
256  /** @brief Multi-component gradient evaluation from DOFs to quadrature points.
257  dof_layout is (DOF x NumComp) and grad_layout is (NIP x DIM x NumComp). */
258  template <typename dof_layout_t, typename dof_data_t,
259  typename grad_layout_t, typename grad_data_t>
260  inline MFEM_ALWAYS_INLINE
261  void CalcGrad(const dof_layout_t &dof_layout,
262  const dof_data_t &dof_data,
263  const grad_layout_t &grad_layout,
264  grad_data_t &grad_data) const
265  {
266  // grad_data(nip,dim,comp) = sum(dof) G(nip,dim,dof) * dof_data(dof,comp)
267  Mult_AB<false>(G_1d.layout, G_1d,
268  dof_layout, dof_data,
269  grad_layout.merge_12(), grad_data);
270  }
271 
272  /** @brief Multi-component gradient evaluation transpose from quadrature points to
273  DOFs. grad_layout is (NIP x DIM x NumComp), dof_layout is (DOF x NumComp). */
274  template <bool Add,
275  typename grad_layout_t, typename grad_data_t,
276  typename dof_layout_t, typename dof_data_t>
277  inline MFEM_ALWAYS_INLINE
278  void CalcGradT(const grad_layout_t &grad_layout,
279  const grad_data_t &grad_data,
280  const dof_layout_t &dof_layout,
281  dof_data_t &dof_data) const
282  {
283  // dof_data(dof,comp) +=
284  // sum(nip,dim) G(nip,dim,dof) * grad_data(nip,dim,comp)
285  Mult_AB<Add>(Gt_1d.layout, Gt_1d,
286  grad_layout.merge_12(), grad_data,
287  dof_layout, dof_data);
288  }
289 
290  /** @brief Multi-component assemble.
291  qpt_layout is (NIP x NumComp), M_layout is (DOF x DOF x NumComp) */
292  template <typename qpt_layout_t, typename qpt_data_t,
293  typename M_layout_t, typename M_data_t>
294  inline MFEM_ALWAYS_INLINE
295  void Assemble(const qpt_layout_t &qpt_layout, const qpt_data_t &qpt_data,
296  const M_layout_t &M_layout, M_data_t &M_data) const
297  {
298  // M_{i,j,k} = \sum_{s} B_1d_{s,i} B_{s,j} qpt_data_{s,k}
299  // Using TensorAssemble: <1,NIP,NC> --> <DOF,1,DOF,NC>
300 #if 0
301  TensorAssemble<false>(
302  B_1d.layout, B_1d,
303  qpt_layout.template split_1<1,NIP>(), qpt_data,
304  M_layout.template split_1<DOF,1>(), M_data);
305 #else
306  TensorAssemble<false>(
307  Bt_1d.layout, Bt_1d, B_1d.layout, B_1d,
308  qpt_layout.template split_1<1,NIP>(), qpt_data,
309  M_layout.template split_1<DOF,1>(), M_data);
310 #endif
311  }
312 
313  /** @brief Multi-component assemble of grad-grad element matrices.
314  qpt_layout is (NIP x DIM x DIM x NumComp), and
315  D_layout is (DOF x DOF x NumComp). */
316  template <typename qpt_layout_t, typename qpt_data_t,
317  typename D_layout_t, typename D_data_t>
318  inline MFEM_ALWAYS_INLINE
319  void AssembleGradGrad(const qpt_layout_t &qpt_layout,
320  const qpt_data_t &qpt_data,
321  const D_layout_t &D_layout,
322  D_data_t &D_data) const
323  {
324  // D_{i,j,k} = \sum_{s} G_1d_{s,i} G_{s,j} qpt_data_{s,k}
325  // Using TensorAssemble: <1,NIP,NC> --> <DOF,1,DOF,NC>
326 #if 0
327  TensorAssemble<false>(
328  G_1d.layout, G_1d,
329  qpt_layout.merge_12().merge_23().template split_1<1,NIP>(), qpt_data,
330  D_layout.template split_1<DOF,1>(), D_data);
331 #else
332  TensorAssemble<false>(
333  Gt_1d.layout, Gt_1d, G_1d.layout, G_1d,
334  qpt_layout.merge_12().merge_23().template split_1<1,NIP>(), qpt_data,
335  D_layout.template split_1<DOF,1>(), D_data);
336 #endif
337  }
338 };
339 
340 /// ShapeEvaluator with 2D tensor-product structure
341 template <int DOF, int NIP, typename real_t>
342 class TProductShapeEvaluator<2, DOF, NIP, real_t>
343 {
344 protected:
347 
348 public:
349  static const int TDOF = DOF*DOF; // total dofs
350  static const int TNIP = NIP*NIP; // total qpts
351 
353 
354  template <bool Dx, bool Dy,
355  typename dof_layout_t, typename dof_data_t,
356  typename qpt_layout_t, typename qpt_data_t>
357  inline MFEM_ALWAYS_INLINE
358  void Calc(const dof_layout_t &dof_layout, const dof_data_t &dof_data,
359  const qpt_layout_t &qpt_layout, qpt_data_t &qpt_data) const
360  {
361  const int NC = dof_layout_t::dim_2;
362  typedef typename qpt_data_t::data_type entry_type;
363  // DOF x DOF x NC --> NIP x DOF x NC --> NIP x NIP x NC
365 
366  // (1) A_{i,j,k} = \sum_s B_1d_{i,s} dof_data_{s,j,k}
367  Mult_2_1<false>(B_1d.layout, Dx ? G_1d : B_1d,
368  dof_layout. template split_1<DOF,DOF>(), dof_data,
369  A.layout, A);
370  // (2) qpt_data_{i,j,k} = \sum_s B_1d_{j,s} A_{i,s,k}
371  Mult_1_2<false>(Bt_1d.layout, Dy ? Gt_1d : Bt_1d,
372  A.layout, A,
373  qpt_layout.template split_1<NIP,NIP>(), qpt_data);
374  }
375 
376  /** @brief Multi-component shape evaluation from DOFs to quadrature points.
377  dof_layout is (TDOF x NumComp) and qpt_layout is (TNIP x NumComp). */
378  template <typename dof_layout_t, typename dof_data_t,
379  typename qpt_layout_t, typename qpt_data_t>
380  inline MFEM_ALWAYS_INLINE
381  void Calc(const dof_layout_t &dof_layout, const dof_data_t &dof_data,
382  const qpt_layout_t &qpt_layout, qpt_data_t &qpt_data) const
383  {
384  Calc<false,false>(dof_layout, dof_data, qpt_layout, qpt_data);
385  }
386 
387  template <bool Dx, bool Dy, bool Add,
388  typename qpt_layout_t, typename qpt_data_t,
389  typename dof_layout_t, typename dof_data_t>
390  inline MFEM_ALWAYS_INLINE
391  void CalcT(const qpt_layout_t &qpt_layout, const qpt_data_t &qpt_data,
392  const dof_layout_t &dof_layout, dof_data_t &dof_data) const
393  {
394  const int NC = dof_layout_t::dim_2;
395  typedef typename qpt_data_t::data_type entry_type;
396  // NIP x NIP X NC --> NIP x DOF x NC --> DOF x DOF x NC
398 
399  // (1) A_{i,j,k} = \sum_s B_1d_{s,j} qpt_data_{i,s,k}
400  Mult_1_2<false>(B_1d.layout, Dy ? G_1d : B_1d,
401  qpt_layout.template split_1<NIP,NIP>(), qpt_data,
402  A.layout, A);
403  // (2) dof_data_{i,j,k} = \sum_s B_1d_{s,i} A_{s,j,k}
404  Mult_2_1<Add>(Bt_1d.layout, Dx ? Gt_1d : Bt_1d,
405  A.layout, A,
406  dof_layout.template split_1<DOF,DOF>(), dof_data);
407  }
408 
409  /** @brief Multi-component shape evaluation transpose from quadrature points to DOFs.
410  qpt_layout is (TNIP x NumComp) and dof_layout is (TDOF x NumComp). */
411  template <bool Add,
412  typename qpt_layout_t, typename qpt_data_t,
413  typename dof_layout_t, typename dof_data_t>
414  inline MFEM_ALWAYS_INLINE
415  void CalcT(const qpt_layout_t &qpt_layout, const qpt_data_t &qpt_data,
416  const dof_layout_t &dof_layout, dof_data_t &dof_data) const
417  {
418  CalcT<false,false,Add>(qpt_layout, qpt_data, dof_layout, dof_data);
419  }
420 
421  /** @brief Multi-component gradient evaluation from DOFs to quadrature points.
422  dof_layout is (TDOF x NumComp) and grad_layout is (TNIP x DIM x NumComp). */
423  template <typename dof_layout_t, typename dof_data_t,
424  typename grad_layout_t, typename grad_data_t>
425  inline MFEM_ALWAYS_INLINE
426  void CalcGrad(const dof_layout_t &dof_layout,
427  const dof_data_t &dof_data,
428  const grad_layout_t &grad_layout,
429  grad_data_t &grad_data) const
430  {
431  Calc<true,false>(dof_layout, dof_data,
432  grad_layout.ind2(0), grad_data);
433  Calc<false,true>(dof_layout, dof_data,
434  grad_layout.ind2(1), grad_data);
435  }
436 
437  /** @brief Multi-component gradient evaluation transpose from quadrature points to
438  DOFs. grad_layout is (TNIP x DIM x NumComp), dof_layout is
439  (TDOF x NumComp). */
440  template <bool Add,
441  typename grad_layout_t, typename grad_data_t,
442  typename dof_layout_t, typename dof_data_t>
443  inline MFEM_ALWAYS_INLINE
444  void CalcGradT(const grad_layout_t &grad_layout,
445  const grad_data_t &grad_data,
446  const dof_layout_t &dof_layout,
447  dof_data_t &dof_data) const
448  {
449  CalcT<true,false, Add>(grad_layout.ind2(0), grad_data,
450  dof_layout, dof_data);
451  CalcT<false,true,true>(grad_layout.ind2(1), grad_data,
452  dof_layout, dof_data);
453  }
454 
455  /** @brief Multi-component assemble.
456  qpt_layout is (TNIP x NumComp), M_layout is (TDOF x TDOF x NumComp) */
457  template <typename qpt_layout_t, typename qpt_data_t,
458  typename M_layout_t, typename M_data_t>
459  inline MFEM_ALWAYS_INLINE
460  void Assemble(const qpt_layout_t &qpt_layout, const qpt_data_t &qpt_data,
461  const M_layout_t &M_layout, M_data_t &M_data) const
462  {
463  const int NC = qpt_layout_t::dim_2;
464  typedef typename qpt_data_t::data_type entry_type;
465 
466  // Using TensorAssemble: <I,NIP,J> --> <DOF,I,DOF,J>
467 
468 #if 0
470  // qpt_data<NIP1,NIP2,NC> --> A<DOF2,NIP1,DOF2,NC>
471  TensorAssemble<false>(
472  B_1d.layout, B_1d,
473  qpt_layout.template split_1<NIP,NIP>(), qpt_data,
474  A.layout, A);
475  // A<DOF2,NIP1,DOF2*NC> --> M<DOF1,DOF2,DOF1,DOF2*NC>
476  TensorAssemble<false>(
477  B_1d.layout, B_1d,
479  M_layout.merge_23().template split_12<DOF,DOF,DOF,DOF*NC>(), M_data);
480 #elif 1
482  // qpt_data<NIP1,NIP2,NC> --> A<DOF2,NIP1,DOF2,NC>
483  TensorAssemble<false>(
484  Bt_1d.layout, Bt_1d, B_1d.layout, B_1d,
485  qpt_layout.template split_1<NIP,NIP>(), qpt_data,
486  A.layout, A);
487  // A<DOF2,NIP1,DOF2*NC> --> M<DOF1,DOF2,DOF1,DOF2*NC>
488  TensorAssemble<false>(
489  Bt_1d.layout, Bt_1d, B_1d.layout, B_1d,
490  A.layout.merge_34(), A,
491  M_layout.merge_23().template split_12<DOF,DOF,DOF,DOF*NC>(), M_data);
492 #else
496  for (int k = 0; k < NC; k++)
497  {
498  // <1,NIP1,NIP2> --> <1,NIP1,NIP2,DOF1>
499  TensorProduct<AssignOp::Set>(
500  qpt_layout.ind2(k).template split_1<NIP,NIP>().
501  template split_1<1,NIP>(), qpt_data,
502  B_1d.layout, B_1d, F3.layout.template split_1<1,NIP>(), F3);
503  // <NIP1,NIP2,DOF1> --> <NIP1,NIP2,DOF1,DOF2>
504  TensorProduct<AssignOp::Set>(
505  F3.layout, F3, B_1d.layout, B_1d, F4.layout, F4);
506  // <NIP1,NIP2,DOF1,DOF2> --> <NIP1,DOF2,DOF1,DOF2>
507  Mult_1_2<false>(B_1d.layout, B_1d,
508  F4.layout.merge_34(), F4,
509  H3.layout, H3);
510  // <NIP1,DOF2,DOF1,DOF2> --> <DOF1,DOF2,DOF1,DOF2>
511  Mult_2_1<false>(Bt_1d.layout, Bt_1d,
512  H3.layout, H3,
513  M_layout.ind3(k).template split_1<DOF,DOF>(),
514  M_data);
515  }
516 #endif
517  }
518 
519  template <int D1, int D2, bool Add,
520  typename qpt_layout_t, typename qpt_data_t,
521  typename D_layout_t, typename D_data_t>
522  inline MFEM_ALWAYS_INLINE
523  void Assemble(const qpt_layout_t &qpt_layout,
524  const qpt_data_t &qpt_data,
525  const D_layout_t &D_layout,
526  D_data_t &D_data) const
527  {
528  const int NC = qpt_layout_t::dim_2;
529  typedef typename qpt_data_t::data_type entry_type;
531 
532  // Using TensorAssemble: <I,NIP,J> --> <DOF,I,DOF,J>
533 
534  // qpt_data<NIP1,NIP2,NC> --> A<DOF2,NIP1,DOF2,NC>
535  TensorAssemble<false>(
536  Bt_1d.layout, D1 == 0 ? Bt_1d : Gt_1d,
537  B_1d.layout, D2 == 0 ? B_1d : G_1d,
538  qpt_layout.template split_1<NIP,NIP>(), qpt_data,
539  A.layout, A);
540  // A<DOF2,NIP1,DOF2*NC> --> M<DOF1,DOF2,DOF1,DOF2*NC>
541  TensorAssemble<Add>(
542  Bt_1d.layout, D1 == 1 ? Bt_1d : Gt_1d,
543  B_1d.layout, D2 == 1 ? B_1d : G_1d,
544  A.layout.merge_34(), A,
545  D_layout.merge_23().template split_12<DOF,DOF,DOF,DOF*NC>(), D_data);
546  }
547 
548  /** @brief Multi-component assemble of grad-grad element matrices.
549  qpt_layout is (TNIP x DIM x DIM x NumComp), and
550  D_layout is (TDOF x TDOF x NumComp). */
551  template <typename qpt_layout_t, typename qpt_data_t,
552  typename D_layout_t, typename D_data_t>
553  inline MFEM_ALWAYS_INLINE
554  void AssembleGradGrad(const qpt_layout_t &qpt_layout,
555  const qpt_data_t &qpt_data,
556  const D_layout_t &D_layout,
557  D_data_t &D_data) const
558  {
559 #if 1
560  Assemble<0,0,false>(qpt_layout.ind23(0,0), qpt_data, D_layout, D_data);
561  Assemble<1,0,true >(qpt_layout.ind23(1,0), qpt_data, D_layout, D_data);
562  Assemble<0,1,true >(qpt_layout.ind23(0,1), qpt_data, D_layout, D_data);
563  Assemble<1,1,true >(qpt_layout.ind23(1,1), qpt_data, D_layout, D_data);
564 #else
565  const int NC = qpt_layout_t::dim_4;
569 
570  for (int k = 0; k < NC; k++)
571  {
572  for (int d1 = 0; d1 < 2; d1++)
573  {
576  // <1,NIP1,NIP2> --> <1,NIP1,NIP2,DOF1>
577  TensorProduct<Set>(qpt_layout.ind23(d1,0).ind2(k).
578  template split_1<NIP,NIP>().
579  template split_1<1,NIP>(), qpt_data,
580  G_1d.layout, G_1d,
581  F3.layout.template split_1<1,NIP>(), F3);
582  // <NIP1,NIP2,DOF1> --> <NIP1,NIP2,DOF1,DOF2>
583  TensorProduct<Set>(F3.layout, F3,
584  B_1d.layout, B_1d,
585  F4.layout, F4);
586  // <1,NIP1,NIP2> --> <1,NIP1,NIP2,DOF1>
587  TensorProduct<Set>(qpt_layout.ind23(d1,1).ind2(k).
588  template split_1<NIP,NIP>().
589  template split_1<1,NIP>(), qpt_data,
590  B_1d.layout, B_1d,
591  F3.layout.template split_1<1,NIP>(), F3);
592  // <NIP1,NIP2,DOF1> --> <NIP1,NIP2,DOF1,DOF2>
593  TensorProduct<Add>(F3.layout, F3,
594  G_1d.layout, G_1d,
595  F4.layout, F4);
596 
597  Mult_1_2<false>(B_1d.layout, d1 == 0 ? B_1d : G_1d,
598  F4.layout.merge_34(), F4,
599  H3.layout, H3);
600  if (d1 == 0)
601  {
602  Mult_2_1<false>(Bt_1d.layout, Gt_1d,
603  H3.layout, H3,
604  D_layout.ind3(k).template split_1<DOF,DOF>(),
605  D_data);
606  }
607  else
608  {
609  Mult_2_1<true>(Bt_1d.layout, Bt_1d,
610  H3.layout, H3,
611  D_layout.ind3(k).template split_1<DOF,DOF>(),
612  D_data);
613  }
614  }
615  }
616 #endif
617  }
618 };
619 
620 /// ShapeEvaluator with 3D tensor-product structure
621 template <int DOF, int NIP, typename real_t>
623 {
624 protected:
627 
628 public:
629  static const int TDOF = DOF*DOF*DOF; // total dofs
630  static const int TNIP = NIP*NIP*NIP; // total qpts
631 
633 
634  template <bool Dx, bool Dy, bool Dz,
635  typename dof_layout_t, typename dof_data_t,
636  typename qpt_layout_t, typename qpt_data_t>
637  inline MFEM_ALWAYS_INLINE
638  void Calc(const dof_layout_t &dof_layout, const dof_data_t &dof_data,
639  const qpt_layout_t &qpt_layout, qpt_data_t &qpt_data) const
640  {
641  const int NC = dof_layout_t::dim_2;
642  typedef typename qpt_data_t::data_type entry_type;
645 
646  // QDD_{i,jj,k} = \sum_s B_1d_{i,s} dof_data_{s,jj,k}
647  Mult_2_1<false>(B_1d.layout, Dx ? G_1d : B_1d,
648  dof_layout.template split_1<DOF,DOF*DOF>(), dof_data,
650  // QQD_{i,j,kk} = \sum_s B_1d_{j,s} QDD_{i,s,kk}
651  Mult_1_2<false>(Bt_1d.layout, Dy ? Gt_1d : Bt_1d,
654  // qpt_data_{ii,j,k} = \sum_s B_1d_{j,s} QQD_{ii,s,k}
655  Mult_1_2<false>(Bt_1d.layout, Dz ? Gt_1d : Bt_1d,
657  qpt_layout.template split_1<NIP*NIP,NIP>(), qpt_data);
658  }
659 
660  /** @brief Multi-component shape evaluation from DOFs to quadrature points.
661  dof_layout is (TDOF x NumComp) and qpt_layout is (TNIP x NumComp). */
662  template <typename dof_layout_t, typename dof_data_t,
663  typename qpt_layout_t, typename qpt_data_t>
664  inline MFEM_ALWAYS_INLINE
665  void Calc(const dof_layout_t &dof_layout, const dof_data_t &dof_data,
666  const qpt_layout_t &qpt_layout, qpt_data_t &qpt_data) const
667  {
668  Calc<false,false,false>(dof_layout, dof_data, qpt_layout, qpt_data);
669  }
670 
671  template <bool Dx, bool Dy, bool Dz, bool Add,
672  typename qpt_layout_t, typename qpt_data_t,
673  typename dof_layout_t, typename dof_data_t>
674  inline MFEM_ALWAYS_INLINE
675  void CalcT(const qpt_layout_t &qpt_layout, const qpt_data_t &qpt_data,
676  const dof_layout_t &dof_layout, dof_data_t &dof_data) const
677  {
678  const int NC = dof_layout_t::dim_2;
679  typedef typename qpt_data_t::data_type entry_type;
682 
683  // QQD_{ii,j,k} = \sum_s B_1d_{s,j} qpt_data_{ii,s,k}
684  Mult_1_2<false>(B_1d.layout, Dz ? G_1d : B_1d,
685  qpt_layout.template split_1<NIP*NIP,NIP>(), qpt_data,
687  // QDD_{i,j,kk} = \sum_s B_1d_{s,j} QQD_{i,s,kk}
688  Mult_1_2<false>(B_1d.layout, Dy ? G_1d : B_1d,
691  // dof_data_{i,jj,k} = \sum_s B_1d_{s,i} QDD_{s,jj,k}
692  Mult_2_1<Add>(Bt_1d.layout, Dx ? Gt_1d : Bt_1d,
694  dof_layout.template split_1<DOF,DOF*DOF>(), dof_data);
695  }
696 
697  /** @brief Multi-component shape evaluation transpose from quadrature points to DOFs.
698  qpt_layout is (TNIP x NumComp) and dof_layout is (TDOF x NumComp). */
699  template <bool Add,
700  typename qpt_layout_t, typename qpt_data_t,
701  typename dof_layout_t, typename dof_data_t>
702  inline MFEM_ALWAYS_INLINE
703  void CalcT(const qpt_layout_t &qpt_layout, const qpt_data_t &qpt_data,
704  const dof_layout_t &dof_layout, dof_data_t &dof_data) const
705  {
706  CalcT<false,false,false,Add>(qpt_layout, qpt_data, dof_layout, dof_data);
707  }
708 
709  /** @brief Multi-component gradient evaluation from DOFs to quadrature points.
710  dof_layout is (TDOF x NumComp) and grad_layout is (TNIP x DIM x NumComp). */
711  template <typename dof_layout_t, typename dof_data_t,
712  typename grad_layout_t, typename grad_data_t>
713  inline MFEM_ALWAYS_INLINE
714  void CalcGrad(const dof_layout_t &dof_layout,
715  const dof_data_t &dof_data,
716  const grad_layout_t &grad_layout,
717  grad_data_t &grad_data) const
718  {
719  Calc<true,false,false>(dof_layout, dof_data,
720  grad_layout.ind2(0), grad_data);
721  Calc<false,true,false>(dof_layout, dof_data,
722  grad_layout.ind2(1), grad_data);
723  Calc<false,false,true>(dof_layout, dof_data,
724  grad_layout.ind2(2), grad_data);
725  // optimization: the x-transition (dof->nip) is done twice -- once for the
726  // y-derivatives and second time for the z-derivatives.
727  }
728 
729  /** @brief Multi-component gradient evaluation transpose from quadrature points to
730  DOFs. grad_layout is (TNIP x DIM x NumComp), dof_layout is
731  (TDOF x NumComp). */
732  template <bool Add,
733  typename grad_layout_t, typename grad_data_t,
734  typename dof_layout_t, typename dof_data_t>
735  inline MFEM_ALWAYS_INLINE
736  void CalcGradT(const grad_layout_t &grad_layout,
737  const grad_data_t &grad_data,
738  const dof_layout_t &dof_layout,
739  dof_data_t &dof_data) const
740  {
741  CalcT<true,false,false, Add>(grad_layout.ind2(0), grad_data,
742  dof_layout, dof_data);
743  CalcT<false,true,false,true>(grad_layout.ind2(1), grad_data,
744  dof_layout, dof_data);
745  CalcT<false,false,true,true>(grad_layout.ind2(2), grad_data,
746  dof_layout, dof_data);
747  }
748 
749  /** @brief Multi-component assemble.
750  qpt_layout is (TNIP x NumComp), M_layout is (TDOF x TDOF x NumComp) */
751  template <typename qpt_layout_t, typename qpt_data_t,
752  typename M_layout_t, typename M_data_t>
753  inline MFEM_ALWAYS_INLINE
754  void Assemble(const qpt_layout_t &qpt_layout, const qpt_data_t &qpt_data,
755  const M_layout_t &M_layout, M_data_t &M_data) const
756  {
757  const int NC = qpt_layout_t::dim_2;
758  typedef typename qpt_data_t::data_type entry_type;
761 
762  // Using TensorAssemble: <I,NIP,J> --> <DOF,I,DOF,J>
763 
764 #if 0
765  // qpt_data<NIP1*NIP2,NIP3,NC> --> A1<DOF3,NIP1*NIP2,DOF3,NC>
766  TensorAssemble<false>(
767  B_1d.layout, B_1d,
768  qpt_layout.template split_1<NIP*NIP,NIP>(), qpt_data,
769  A1.layout, A1);
770  // A1<DOF3*NIP1,NIP2,DOF3*NC> --> A2<DOF2,DOF3*NIP1,DOF2,DOF3*NC>
771  TensorAssemble<false>(
772  B_1d.layout, B_1d,
774  A2.layout, A2);
775  // A2<DOF2*DOF3,NIP1,DOF2*DOF3*NC> --> M<DOF1,DOF2*DOF3,DOF1,DOF2*DOF3*NC>
776  TensorAssemble<false>(
777  B_1d.layout, B_1d,
779  M_layout.merge_23().template split_12<DOF,DOF*DOF,DOF,DOF*DOF*NC>(),
780  M_data);
781 #else
782  // qpt_data<NIP1*NIP2,NIP3,NC> --> A1<DOF3,NIP1*NIP2,DOF3,NC>
783  TensorAssemble<false>(
784  Bt_1d.layout, Bt_1d, B_1d.layout, B_1d,
785  qpt_layout.template split_1<NIP*NIP,NIP>(), qpt_data,
786  A1.layout, A1);
787  // A1<DOF3*NIP1,NIP2,DOF3*NC> --> A2<DOF2,DOF3*NIP1,DOF2,DOF3*NC>
788  TensorAssemble<false>(
789  Bt_1d.layout, Bt_1d, B_1d.layout, B_1d,
791  A2.layout, A2);
792  // A2<DOF2*DOF3,NIP1,DOF2*DOF3*NC> --> M<DOF1,DOF2*DOF3,DOF1,DOF2*DOF3*NC>
793  TensorAssemble<false>(
794  Bt_1d.layout, Bt_1d, B_1d.layout, B_1d,
796  M_layout.merge_23().template split_12<DOF,DOF*DOF,DOF,DOF*DOF*NC>(),
797  M_data);
798 #endif
799  }
800 
801  template <int D1, int D2, bool Add,
802  typename qpt_layout_t, typename qpt_data_t,
803  typename D_layout_t, typename D_data_t>
804  inline MFEM_ALWAYS_INLINE
805  void Assemble(const qpt_layout_t &qpt_layout,
806  const qpt_data_t &qpt_data,
807  const D_layout_t &D_layout,
808  D_data_t &D_data) const
809  {
810  const int NC = qpt_layout_t::dim_2;
811  typedef typename qpt_data_t::data_type entry_type;
814 
815  // Using TensorAssemble: <I,NIP,J> --> <DOF,I,DOF,J>
816 
817  // qpt_data<NIP1*NIP2,NIP3,NC> --> A1<DOF3,NIP1*NIP2,DOF3,NC>
818  TensorAssemble<false>(
819  Bt_1d.layout, D1 != 2 ? Bt_1d : Gt_1d,
820  B_1d.layout, D2 != 2 ? B_1d : G_1d,
821  qpt_layout.template split_1<NIP*NIP,NIP>(), qpt_data,
822  A1.layout, A1);
823  // A1<DOF3*NIP1,NIP2,DOF3*NC> --> A2<DOF2,DOF3*NIP1,DOF2,DOF3*NC>
824  TensorAssemble<false>(
825  Bt_1d.layout, D1 != 1 ? Bt_1d : Gt_1d,
826  B_1d.layout, D2 != 1 ? B_1d : G_1d,
828  A2.layout, A2);
829  // A2<DOF2*DOF3,NIP1,DOF2*DOF3*NC> --> M<DOF1,DOF2*DOF3,DOF1,DOF2*DOF3*NC>
830  TensorAssemble<Add>(
831  Bt_1d.layout, D1 != 0 ? Bt_1d : Gt_1d,
832  B_1d.layout, D2 != 0 ? B_1d : G_1d,
834  D_layout.merge_23().template split_12<DOF,DOF*DOF,DOF,DOF*DOF*NC>(),
835  D_data);
836  }
837 
838 #if 0
839  template <typename qpt_layout_t, typename qpt_data_t,
840  typename D_layout_t, typename D_data_t>
841  inline MFEM_ALWAYS_INLINE
842  void Assemble(int D1, int D2,
843  const qpt_layout_t &qpt_layout,
844  const qpt_data_t &qpt_data,
845  const D_layout_t &D_layout,
846  D_data_t &D_data) const
847  {
848  const int NC = qpt_layout_t::dim_2;
851 
852  // Using TensorAssemble: <I,NIP,J> --> <DOF,I,DOF,J>
853 
854  // qpt_data<NIP1*NIP2,NIP3,NC> --> A1<DOF3,NIP1*NIP2,DOF3,NC>
855  TensorAssemble<false>(
856  Bt_1d.layout, D1 != 2 ? Bt_1d : Gt_1d,
857  B_1d.layout, D2 != 2 ? B_1d : G_1d,
858  qpt_layout.template split_1<NIP*NIP,NIP>(), qpt_data,
859  A1.layout, A1);
860  // A1<DOF3*NIP1,NIP2,DOF3*NC> --> A2<DOF2,DOF3*NIP1,DOF2,DOF3*NC>
861  TensorAssemble<false>(
862  Bt_1d.layout, D1 != 1 ? Bt_1d : Gt_1d,
863  B_1d.layout, D2 != 1 ? B_1d : G_1d,
865  A2.layout, A2);
866  // A2<DOF2*DOF3,NIP1,DOF2*DOF3*NC> --> M<DOF1,DOF2*DOF3,DOF1,DOF2*DOF3*NC>
867  TensorAssemble<true>(
868  Bt_1d.layout, D1 != 0 ? Bt_1d : Gt_1d,
869  B_1d.layout, D2 != 0 ? B_1d : G_1d,
871  D_layout.merge_23().template split_12<DOF,DOF*DOF,DOF,DOF*DOF*NC>(),
872  D_data);
873  }
874 #endif
875 
876  /** @brief Multi-component assemble of grad-grad element matrices.
877  qpt_layout is (TNIP x DIM x DIM x NumComp), and
878  D_layout is (TDOF x TDOF x NumComp). */
879  template <typename qpt_layout_t, typename qpt_data_t,
880  typename D_layout_t, typename D_data_t>
881  inline MFEM_ALWAYS_INLINE
882  void AssembleGradGrad(const qpt_layout_t &qpt_layout,
883  const qpt_data_t &qpt_data,
884  const D_layout_t &D_layout,
885  D_data_t &D_data) const
886  {
887 #if 1
888  // NOTE: This function compiles into a large chunk of machine code
889  Assemble<0,0,false>(qpt_layout.ind23(0,0), qpt_data, D_layout, D_data);
890  Assemble<1,0,true >(qpt_layout.ind23(1,0), qpt_data, D_layout, D_data);
891  Assemble<2,0,true >(qpt_layout.ind23(2,0), qpt_data, D_layout, D_data);
892  Assemble<0,1,true >(qpt_layout.ind23(0,1), qpt_data, D_layout, D_data);
893  Assemble<1,1,true >(qpt_layout.ind23(1,1), qpt_data, D_layout, D_data);
894  Assemble<2,1,true >(qpt_layout.ind23(2,1), qpt_data, D_layout, D_data);
895  Assemble<0,2,true >(qpt_layout.ind23(0,2), qpt_data, D_layout, D_data);
896  Assemble<1,2,true >(qpt_layout.ind23(1,2), qpt_data, D_layout, D_data);
897  Assemble<2,2,true >(qpt_layout.ind23(2,2), qpt_data, D_layout, D_data);
898 #else
899  TAssign<AssignOp::Set>(D_layout, D_data, 0.0);
900  for (int d2 = 0; d2 < 3; d2++)
901  {
902  for (int d1 = 0; d1 < 3; d1++)
903  {
904  Assemble(d1, d2, qpt_layout.ind23(d1,d2), qpt_data,
905  D_layout, D_data);
906  }
907  }
908 #endif
909  }
910 };
911 
912 /// ShapeEvaluator with tensor-product structure in any dimension
913 template <class FE, class IR, typename real_t>
914 class ShapeEvaluator_base<FE, IR, true, real_t>
915  : public TProductShapeEvaluator<FE::dim, FE::dofs_1d, IR::qpts_1d, real_t>
916 {
917 protected:
918  typedef TProductShapeEvaluator<FE::dim,FE::dofs_1d,
919  IR::qpts_1d,real_t> base_class;
920  using base_class::B_1d;
921  using base_class::Bt_1d;
922  using base_class::G_1d;
923  using base_class::Gt_1d;
924 
925 public:
926  ShapeEvaluator_base(const FE &fe)
927  {
928  fe.Calc1DShapes(IR::Get1DIntRule(), B_1d.data, G_1d.data);
929  TAssign<AssignOp::Set>(Bt_1d.layout, Bt_1d,
930  B_1d.layout.transpose_12(), B_1d);
931  TAssign<AssignOp::Set>(Gt_1d.layout, Gt_1d,
932  G_1d.layout.transpose_12(), G_1d);
933  }
934 
935  // default copy constructor
936 };
937 
938 /// General ShapeEvaluator for any scalar FE type (L2 or H1)
939 template <class FE, class IR, typename real_t>
941  : public ShapeEvaluator_base<FE,IR,FE::tensor_prod && IR::tensor_prod,real_t>
942 {
943 public:
944  typedef real_t real_type;
945  static const int dim = FE::dim;
946  static const int qpts = IR::qpts;
947  static const bool tensor_prod = FE::tensor_prod && IR::tensor_prod;
948  typedef FE FE_type;
949  typedef IR IR_type;
951 
952  using base_class::Calc;
953  using base_class::CalcT;
954  using base_class::CalcGrad;
955  using base_class::CalcGradT;
956 
957  ShapeEvaluator(const FE &fe) : base_class(fe) { }
958 
959  // default copy constructor
960 };
961 
962 
963 /** @brief Field evaluators -- values of a given global FE grid function
964  This is roughly speaking a templated version of GridFunction
965 */
966 template <typename FESpace_t, typename VecLayout_t, typename IR,
967  typename complex_t, typename real_t>
969 {
970 protected:
971  typedef typename FESpace_t::FE_type FE_type;
973 
974  FESpace_t fespace;
976  VecLayout_t vec_layout;
977 
978  /// With this constructor, fespace is a shallow copy.
979  inline MFEM_ALWAYS_INLINE
980  FieldEvaluator_base(const FESpace_t &tfes, const ShapeEval_type &shape_eval,
981  const VecLayout_t &vec_layout)
982  : fespace(tfes),
983  shapeEval(shape_eval),
985  { }
986 
987  /// This constructor creates new fespace, not a shallow copy.
988  inline MFEM_ALWAYS_INLINE
990  : fespace(fe, fes), shapeEval(fe), vec_layout(fes)
991  { }
992 };
993 
994 /// complex_t - dof/qpt data type, real_t - ShapeEvaluator (FE basis) data type
995 template <typename FESpace_t, typename VecLayout_t, typename IR,
996  typename complex_t = double, typename real_t = double>
998  : public FieldEvaluator_base<FESpace_t,VecLayout_t,IR,complex_t,real_t>
999 {
1000 public:
1001  typedef complex_t complex_type;
1002  typedef FESpace_t FESpace_type;
1003  typedef typename FESpace_t::FE_type FE_type;
1005  typedef VecLayout_t VecLayout_type;
1006 
1007  // this type
1009 
1010  static const int dofs = FE_type::dofs;
1011  static const int dim = FE_type::dim;
1012  static const int qpts = IR::qpts;
1013  static const int vdim = VecLayout_t::vec_dim;
1014 
1015 protected:
1016 
1019 
1020  using base_class::fespace;
1021  using base_class::shapeEval;
1022  using base_class::vec_layout;
1023  const complex_t *data_in;
1024  complex_t *data_out;
1025 
1026 public:
1027  /// With this constructor, fespace is a shallow copy of tfes.
1028  inline MFEM_ALWAYS_INLINE
1029  FieldEvaluator(const FESpace_t &tfes, const ShapeEval_type &shape_eval,
1030  const VecLayout_type &vec_layout,
1031  const complex_t *global_data_in, complex_t *global_data_out)
1032  : base_class(tfes, shape_eval, vec_layout),
1033  data_in(global_data_in),
1034  data_out(global_data_out)
1035  { }
1036 
1037  /// With this constructor, fespace is a shallow copy of f.fespace.
1038  inline MFEM_ALWAYS_INLINE
1040  const complex_t *global_data_in, complex_t *global_data_out)
1042  data_in(global_data_in),
1043  data_out(global_data_out)
1044  { }
1045 
1046  /// This constructor creates a new fespace, not a shallow copy.
1047  inline MFEM_ALWAYS_INLINE
1049  const complex_t *global_data_in, complex_t *global_data_out)
1050  : base_class(FE_type(*fes.FEColl()), fes),
1051  data_in(global_data_in),
1052  data_out(global_data_out)
1053  { }
1054 
1055  // Default copy constructor
1056 
1057  inline MFEM_ALWAYS_INLINE FESpace_type &FESpace() { return fespace; }
1058  inline MFEM_ALWAYS_INLINE ShapeEval_type &ShapeEval() { return shapeEval; }
1059  inline MFEM_ALWAYS_INLINE VecLayout_type &VecLayout() { return vec_layout; }
1060 
1061  inline MFEM_ALWAYS_INLINE
1062  void SetElement(int el)
1063  {
1064  fespace.SetElement(el);
1065  }
1066 
1067  /// val_layout_t is (qpts x vdim x NE)
1068  template <typename val_layout_t, typename val_data_t>
1069  inline MFEM_ALWAYS_INLINE
1070  void GetValues(int el, const val_layout_t &l, val_data_t &vals)
1071  {
1072  const int ne = val_layout_t::dim_3;
1074  SetElement(el);
1075  fespace.VectorExtract(vec_layout, data_in, val_dofs.layout, val_dofs);
1076  shapeEval.Calc(val_dofs.layout.merge_23(), val_dofs, l.merge_23(), vals);
1077  }
1078 
1079  /// grad_layout_t is (qpts x dim x vdim x NE)
1080  template <typename grad_layout_t, typename grad_data_t>
1081  inline MFEM_ALWAYS_INLINE
1082  void GetGradients(int el, const grad_layout_t &l, grad_data_t &grad)
1083  {
1084  const int ne = grad_layout_t::dim_4;
1086  SetElement(el);
1087  fespace.VectorExtract(vec_layout, data_in, val_dofs.layout, val_dofs);
1088  shapeEval.CalcGrad(val_dofs.layout.merge_23(), val_dofs,
1089  l.merge_34(), grad);
1090  }
1091 
1092  // TODO: add method GetValuesAndGradients()
1093 
1094  template <typename DataType>
1095  inline MFEM_ALWAYS_INLINE
1096  void Eval(DataType &F)
1097  {
1098  // T.SetElement() must be called outside
1100  }
1101 
1102  template <typename DataType>
1103  inline MFEM_ALWAYS_INLINE
1104  void Eval(int el, DataType &F)
1105  {
1106  SetElement(el);
1107  Eval(F);
1108  }
1109 
1110  template <bool Add, typename DataType>
1111  inline MFEM_ALWAYS_INLINE
1112  void Assemble(DataType &F)
1113  {
1114  // T.SetElement() must be called outside
1116  template Assemble<Add>(vec_layout, *this, F);
1117  }
1118 
1119  template <bool Add, typename DataType>
1120  inline MFEM_ALWAYS_INLINE
1121  void Assemble(int el, DataType &F)
1122  {
1123  SetElement(el);
1124  Assemble<Add>(F);
1125  }
1126 
1127 #ifdef MFEM_TEMPLATE_ENABLE_SERIALIZE
1128  template <typename DataType>
1129  inline MFEM_ALWAYS_INLINE
1130  void EvalSerialized(const typename DataType::vcomplex_t *loc_dofs,
1131  DataType &F)
1132  {
1134  }
1135 
1136  template <bool Add, typename DataType>
1137  inline MFEM_ALWAYS_INLINE
1138  void AssembleSerialized(const DataType &F,
1139  typename DataType::vcomplex_t *loc_dofs)
1140  {
1142  template AssembleSerialized<Add>(*this, F, loc_dofs);
1143  }
1144 #endif
1145 
1146  /** @brief Enumeration for the data type used by the Eval() and Assemble() methods.
1147  The types can be obtained by summing constants from this enumeration and used
1148  as a template parameter in struct Data. */
1150  {
1151  None = 0,
1152  Values = 1,
1154  };
1155 
1156  /** @brief Auxiliary templated struct AData, used by the Eval() and Assemble()
1157  methods.
1158 
1159  The template parameter IOData is "bitwise or" of constants from
1160  the enum InOutData. The parameter NE is the number of elements to be
1161  processed in the Eval() and Assemble() methods. */
1162  template<int IOData, typename impl_traits_t> struct AData;
1163 
1164  template <typename it_t> struct AData<0,it_t> // 0 = None
1165  {
1166  // Do we need this?
1167  };
1168 
1169  template <typename it_t> struct AData<1,it_t> // 1 = Values
1170  {
1171  static const int ne = it_t::batch_size;
1172  typedef typename it_t::vcomplex_t vcomplex_t;
1173 #ifdef MFEM_TEMPLATE_FIELD_EVAL_DATA_HAS_DOFS
1176 #else
1178 #endif
1180  };
1181 
1182  template <typename it_t> struct AData<2,it_t> // 2 = Gradients
1183  {
1184  static const int ne = it_t::batch_size;
1185  typedef typename it_t::vcomplex_t vcomplex_t;
1186 #ifdef MFEM_TEMPLATE_FIELD_EVAL_DATA_HAS_DOFS
1189 #else
1191 #endif
1193  };
1194 
1195  template <typename it_t> struct AData<3,it_t> // 3 = Values+Gradients
1196  {
1197  static const int ne = it_t::batch_size;
1198  typedef typename it_t::vcomplex_t vcomplex_t;
1199 #ifdef MFEM_TEMPLATE_FIELD_EVAL_DATA_HAS_DOFS
1202 #else
1204 #endif
1207  };
1208 
1209  /** @brief This struct is similar to struct AData, adding separate static data
1210  members for the input (InData) and output (OutData) data types. */
1211  template <int IData, int OData, typename it_t>
1212  struct BData : public AData<IData|OData,it_t>
1213  {
1215  static const int InData = IData;
1216  static const int OutData = OData;
1217  };
1218 
1219  /** @brief This struct implements the input (Eval, EvalSerialized) and output
1220  (Assemble, AssembleSerialized) operations for the given Ops.
1221  Ops is "bitwise or" of constants from the enum InOutData. */
1222  template <int Ops, bool dummy> struct Action;
1223 
1224  template <bool dummy> struct Action<0,dummy> // 0 = None
1225  {
1226  // Do we need this?
1227  };
1228 
1229  template <bool dummy> struct Action<1,dummy> // 1 = Values
1230  {
1231  template <typename vec_layout_t, typename AData_t>
1232  static inline MFEM_ALWAYS_INLINE
1233  void Eval(const vec_layout_t &l, T_type &T, AData_t &D)
1234  {
1235 #ifdef MFEM_TEMPLATE_FIELD_EVAL_DATA_HAS_DOFS
1236  typename AData_t::val_dofs_t &val_dofs = D.val_dofs;
1237 #else
1238  typename AData_t::val_dofs_t val_dofs;
1239 #endif
1240  T.fespace.VectorExtract(l, T.data_in, val_dofs.layout, val_dofs);
1241  T.shapeEval.Calc(val_dofs.layout.merge_23(), val_dofs,
1242  D.val_qpts.layout.merge_23(), D.val_qpts);
1243  }
1244 
1245  template <bool Add, typename vec_layout_t, typename AData_t>
1246  static inline MFEM_ALWAYS_INLINE
1247  void Assemble(const vec_layout_t &l, T_type &T, AData_t &D)
1248  {
1250 #ifdef MFEM_TEMPLATE_FIELD_EVAL_DATA_HAS_DOFS
1251  typename AData_t::val_dofs_t &val_dofs = D.val_dofs;
1252 #else
1253  typename AData_t::val_dofs_t val_dofs;
1254 #endif
1255  T.shapeEval.template CalcT<false>(
1256  D.val_qpts.layout.merge_23(), D.val_qpts,
1257  val_dofs.layout.merge_23(), val_dofs);
1258  T.fespace.template VectorAssemble<Op>(
1259  val_dofs.layout, val_dofs, l, T.data_out);
1260  }
1261 
1262 #ifdef MFEM_TEMPLATE_ENABLE_SERIALIZE
1263  template <typename AData_t>
1264  static inline MFEM_ALWAYS_INLINE
1266  const typename AData_t::vcomplex_t *loc_dofs,
1267  AData_t &D)
1268  {
1269  T.shapeEval.Calc(AData_t::val_dofs_t::layout.merge_23(), loc_dofs,
1270  D.val_qpts.layout.merge_23(), D.val_qpts);
1271  }
1272 
1273  template <bool Add, typename AData_t>
1274  static inline MFEM_ALWAYS_INLINE
1275  void AssembleSerialized(T_type &T, const AData_t &D,
1276  typename AData_t::vcomplex_t *loc_dofs)
1277  {
1278  T.shapeEval.template CalcT<Add>(
1279  D.val_qpts.layout.merge_23(), D.val_qpts,
1280  AData_t::val_dofs_t::layout.merge_23(), loc_dofs);
1281  }
1282 #endif
1283  };
1284 
1285  template <bool dummy> struct Action<2,dummy> // 2 = Gradients
1286  {
1287  template <typename vec_layout_t, typename AData_t>
1288  static inline MFEM_ALWAYS_INLINE
1289  void Eval(const vec_layout_t &l, T_type &T, AData_t &D)
1290  {
1291 #ifdef MFEM_TEMPLATE_FIELD_EVAL_DATA_HAS_DOFS
1292  typename AData_t::val_dofs_t &val_dofs = D.val_dofs;
1293 #else
1294  typename AData_t::val_dofs_t val_dofs;
1295 #endif
1296  T.fespace.VectorExtract(l, T.data_in, val_dofs.layout, val_dofs);
1297  T.shapeEval.CalcGrad(val_dofs.layout.merge_23(), val_dofs,
1298  D.grad_qpts.layout.merge_34(), D.grad_qpts);
1299  }
1300 
1301  template <bool Add, typename vec_layout_t, typename AData_t>
1302  static inline MFEM_ALWAYS_INLINE
1303  void Assemble(const vec_layout_t &l, T_type &T, AData_t &D)
1304  {
1306 #ifdef MFEM_TEMPLATE_FIELD_EVAL_DATA_HAS_DOFS
1307  typename AData_t::val_dofs_t &val_dofs = D.val_dofs;
1308 #else
1309  typename AData_t::val_dofs_t val_dofs;
1310 #endif
1311  T.shapeEval.template CalcGradT<false>(
1312  D.grad_qpts.layout.merge_34(), D.grad_qpts,
1313  val_dofs.layout.merge_23(), val_dofs);
1314  T.fespace.template VectorAssemble<Op>(
1315  val_dofs.layout, val_dofs, l, T.data_out);
1316  }
1317 
1318 #ifdef MFEM_TEMPLATE_ENABLE_SERIALIZE
1319  template <typename AData_t>
1320  static inline MFEM_ALWAYS_INLINE
1322  const typename AData_t::vcomplex_t *loc_dofs,
1323  AData_t &D)
1324  {
1325  T.shapeEval.CalcGrad(AData_t::val_dofs_t::layout.merge_23(), loc_dofs,
1326  D.grad_qpts.layout.merge_34(), D.grad_qpts);
1327  }
1328 
1329  template <bool Add, typename AData_t>
1330  static inline MFEM_ALWAYS_INLINE
1331  void AssembleSerialized(T_type &T, const AData_t &D,
1332  typename AData_t::vcomplex_t *loc_dofs)
1333  {
1334  T.shapeEval.template CalcGradT<Add>(
1335  D.grad_qpts.layout.merge_34(), D.grad_qpts,
1336  AData_t::val_dofs_t::layout.merge_23(), loc_dofs);
1337  }
1338 #endif
1339  };
1340 
1341  template <bool dummy> struct Action<3,dummy> // 3 = Values+Gradients
1342  {
1343  template <typename vec_layout_t, typename AData_t>
1344  static inline MFEM_ALWAYS_INLINE
1345  void Eval(const vec_layout_t &l, T_type &T, AData_t &D)
1346  {
1347 #ifdef MFEM_TEMPLATE_FIELD_EVAL_DATA_HAS_DOFS
1348  typename AData_t::val_dofs_t &val_dofs = D.val_dofs;
1349 #else
1350  typename AData_t::val_dofs_t val_dofs;
1351 #endif
1352  T.fespace.VectorExtract(l, T.data_in, val_dofs.layout, val_dofs);
1353  T.shapeEval.Calc(val_dofs.layout.merge_23(), val_dofs,
1354  D.val_qpts.layout.merge_23(), D.val_qpts);
1355  T.shapeEval.CalcGrad(val_dofs.layout.merge_23(), val_dofs,
1356  D.grad_qpts.layout.merge_34(), D.grad_qpts);
1357  }
1358 
1359  template <bool Add, typename vec_layout_t, typename AData_t>
1360  static inline MFEM_ALWAYS_INLINE
1361  void Assemble(const vec_layout_t &l, T_type &T, AData_t &D)
1362  {
1364 #ifdef MFEM_TEMPLATE_FIELD_EVAL_DATA_HAS_DOFS
1365  typename AData_t::val_dofs_t &val_dofs = D.val_dofs;
1366 #else
1367  typename AData_t::val_dofs_t val_dofs;
1368 #endif
1369  T.shapeEval.template CalcT<false>(
1370  D.val_qpts.layout.merge_23(), D.val_qpts,
1371  val_dofs.layout.merge_23(), val_dofs);
1372  T.shapeEval.template CalcGradT<true>(
1373  D.grad_qpts.layout.merge_34(), D.grad_qpts,
1374  val_dofs.layout.merge_23(), val_dofs);
1375  T.fespace.template VectorAssemble<Op>(
1376  val_dofs.layout, val_dofs, l, T.data_out);
1377  }
1378 
1379 #ifdef MFEM_TEMPLATE_ENABLE_SERIALIZE
1380  template <typename AData_t>
1381  static inline MFEM_ALWAYS_INLINE
1383  const typename AData_t::vcomplex_t *loc_dofs,
1384  AData_t &D)
1385  {
1386  T.shapeEval.Calc(AData_t::val_dofs_t::layout.merge_23(), loc_dofs,
1387  D.val_qpts.layout.merge_23(), D.val_qpts);
1388  T.shapeEval.CalcGrad(AData_t::val_dofs_t::layout.merge_23(), loc_dofs,
1389  D.grad_qpts.layout.merge_34(), D.grad_qpts);
1390  }
1391 
1392  template <bool Add, typename AData_t>
1393  static inline MFEM_ALWAYS_INLINE
1394  void AssembleSerialized(T_type &T, const AData_t &D,
1395  typename AData_t::vcomplex_t *loc_dofs)
1396  {
1397  T.shapeEval.template CalcT<Add>(
1398  D.val_qpts.layout.merge_23(), D.val_qpts,
1399  AData_t::val_dofs_t::layout.merge_23(), loc_dofs);
1400  T.shapeEval.template CalcGradT<true>(
1401  D.grad_qpts.layout.merge_34(), D.grad_qpts,
1402  AData_t::val_dofs_t::layout.merge_23(), loc_dofs);
1403  }
1404 #endif
1405  };
1406 
1407  /** @brief This struct implements element matrix computation for some combinations
1408  of input (InOps) and output (OutOps) operations. */
1409  template <int InOps, int OutOps, typename it_t> struct TElementMatrix;
1410 
1411  // Case 1,1 = Values,Values
1412  template <typename it_t> struct TElementMatrix<1,1,it_t>
1413  {
1414  // qpt_layout_t is (nip), M_layout_t is (dof x dof)
1415  // it_t::batch_size = 1 is assumed
1416  template <typename qpt_layout_t, typename qpt_data_t,
1417  typename M_layout_t, typename M_data_t>
1418  static inline MFEM_ALWAYS_INLINE
1419  void Compute(const qpt_layout_t &a, const qpt_data_t &A,
1420  const M_layout_t &m, M_data_t &M, ShapeEval_type &ev)
1421  {
1422  ev.Assemble(a.template split_1<qpts,1>(), A,
1423  m.template split_2<dofs,1>(), M);
1424  }
1425  };
1426 
1427  // Case 2,2 = Gradients,Gradients
1428  template <typename it_t> struct TElementMatrix<2,2,it_t>
1429  {
1430  /** @brief Assemble element mass matrix
1431  @param a the layout for the quadrature point data
1432  @param A given quadrature point data for element (incl. coefficient,
1433  geometry)
1434  @param m the layout for the resulting element mass matrix
1435  @param M the resulting element mass matrix
1436  @param ev the shape evaluator
1437  qpt_layout_t is (nip), M_layout_t is (dof x dof)
1438  NE = 1 is assumed */
1439  template <typename qpt_layout_t, typename qpt_data_t,
1440  typename M_layout_t, typename M_data_t>
1441  static inline MFEM_ALWAYS_INLINE
1442  void Compute(const qpt_layout_t &a, const qpt_data_t &A,
1443  const M_layout_t &m, M_data_t &M, ShapeEval_type &ev)
1444  {
1445  ev.AssembleGradGrad(a.template split_3<dim,1>(), A,
1446  m.template split_2<dofs,1>(), M);
1447  }
1448  };
1449 
1450  template <typename kernel_t, typename impl_traits_t> struct Spec
1451  {
1452  static const int InData =
1453  Values*kernel_t::in_values + Gradients*kernel_t::in_gradients;
1454  static const int OutData =
1455  Values*kernel_t::out_values + Gradients*kernel_t::out_gradients;
1456 
1459  };
1460 };
1461 
1462 } // namespace mfem
1463 
1464 #endif // MFEM_TEMPLATE_EVALUATOR
MFEM_ALWAYS_INLINE void Assemble(const qpt_layout_t &qpt_layout, const qpt_data_t &qpt_data, const D_layout_t &D_layout, D_data_t &D_data) const
Definition: tevaluator.hpp:523
TMatrix< NIP, DOF, real_t, true > G_1d
Definition: tevaluator.hpp:223
MFEM_ALWAYS_INLINE void Assemble(int D1, int D2, const qpt_layout_t &qpt_layout, const qpt_data_t &qpt_data, const D_layout_t &D_layout, D_data_t &D_data) const
Definition: tevaluator.hpp:842
static const bool tensor_prod
Definition: tevaluator.hpp:947
MFEM_ALWAYS_INLINE void CalcT(const qpt_layout_t &qpt_layout, const qpt_data_t &qpt_data, const dof_layout_t &dof_layout, dof_data_t &dof_data) const
Multi-component shape evaluation transpose from quadrature points to DOFs. qpt_layout is (NIP x NumCo...
Definition: tevaluator.hpp:248
MFEM_ALWAYS_INLINE void CalcGradT(const grad_layout_t &grad_layout, const grad_data_t &grad_data, const dof_layout_t &dof_layout, dof_data_t &dof_data) const
Multi-component gradient evaluation transpose from quadrature points to DOFs. grad_layout is (NIP x D...
Definition: tevaluator.hpp:137
TTensor4< qpts, dim, vdim, ne, vcomplex_t > grad_qpts
static MFEM_ALWAYS_INLINE void AssembleSerialized(T_type &T, const AData_t &D, typename AData_t::vcomplex_t *loc_dofs)
MFEM_ALWAYS_INLINE void AssembleSerialized(const DataType &F, typename DataType::vcomplex_t *loc_dofs)
MFEM_ALWAYS_INLINE void CalcGrad(const dof_layout_t &dof_layout, const dof_data_t &dof_data, const grad_layout_t &grad_layout, grad_data_t &grad_data) const
Multi-component gradient evaluation from DOFs to quadrature points. dof_layout is (DOF x NumComp) and...
Definition: tevaluator.hpp:261
MFEM_ALWAYS_INLINE VecLayout_type & VecLayout()
MFEM_ALWAYS_INLINE void CalcT(const qpt_layout_t &qpt_layout, const qpt_data_t &qpt_data, const dof_layout_t &dof_layout, dof_data_t &dof_data) const
Multi-component shape evaluation transpose from quadrature points to DOFs. qpt_layout is (NIP x NumCo...
Definition: tevaluator.hpp:89
ShapeEvaluator< FE_type, IR, real_t > ShapeEval_type
Definition: tevaluator.hpp:972
This struct implements element matrix computation for some combinations of input (InOps) and output (...
TMatrix< DOF, NIP, real_t, true > Bt
Definition: tevaluator.hpp:46
VecLayout_t VecLayout_type
MFEM_ALWAYS_INLINE void Calc(const dof_layout_t &dof_layout, const dof_data_t &dof_data, const qpt_layout_t &qpt_layout, qpt_data_t &qpt_data) const
Multi-component shape evaluation from DOFs to quadrature points. dof_layout is (DOF x NumComp) and qp...
Definition: tevaluator.hpp:66
static MFEM_ALWAYS_INLINE void Eval(const vec_layout_t &l, T_type &T, AData_t &D)
TMatrix< DOF, NIP, real_t, true > Gt_1d
Definition: tevaluator.hpp:346
static MFEM_ALWAYS_INLINE void Compute(const qpt_layout_t &a, const qpt_data_t &A, const M_layout_t &m, M_data_t &M, ShapeEval_type &ev)
MFEM_ALWAYS_INLINE void Calc(const dof_layout_t &dof_layout, const dof_data_t &dof_data, const qpt_layout_t &qpt_layout, qpt_data_t &qpt_data) const
Multi-component shape evaluation from DOFs to quadrature points. dof_layout is (TDOF x NumComp) and q...
Definition: tevaluator.hpp:381
static MFEM_ALWAYS_INLINE void Compute(const qpt_layout_t &a, const qpt_data_t &A, const M_layout_t &m, M_data_t &M, ShapeEval_type &ev)
Assemble element mass matrix.
TProductShapeEvaluator< FE::dim, FE::dofs_1d, IR::qpts_1d, real_t > base_class
Definition: tevaluator.hpp:919
MFEM_ALWAYS_INLINE void Calc(const dof_layout_t &dof_layout, const dof_data_t &dof_data, const qpt_layout_t &qpt_layout, qpt_data_t &qpt_data) const
Multi-component shape evaluation from DOFs to quadrature points. dof_layout is (TDOF x NumComp) and q...
Definition: tevaluator.hpp:665
static const layout_type layout
Definition: ttensor.hpp:356
TTensor3< qpts, vdim, ne, vcomplex_t > val_qpts
MFEM_ALWAYS_INLINE void CalcGrad(const dof_layout_t &dof_layout, const dof_data_t &dof_data, const grad_layout_t &grad_layout, grad_data_t &grad_data) const
Multi-component gradient evaluation from DOFs to quadrature points. dof_layout is (DOF x NumComp) and...
Definition: tevaluator.hpp:111
MFEM_ALWAYS_INLINE void Eval(DataType &F)
MFEM_ALWAYS_INLINE void CalcGrad(const dof_layout_t &dof_layout, const dof_data_t &dof_data, const grad_layout_t &grad_layout, grad_data_t &grad_data) const
Multi-component gradient evaluation from DOFs to quadrature points. dof_layout is (TDOF x NumComp) an...
Definition: tevaluator.hpp:426
MFEM_ALWAYS_INLINE void CalcGradT(const grad_layout_t &grad_layout, const grad_data_t &grad_data, const dof_layout_t &dof_layout, dof_data_t &dof_data) const
Multi-component gradient evaluation transpose from quadrature points to DOFs. grad_layout is (TNIP x ...
Definition: tevaluator.hpp:736
TMatrix< DOF, NIP, real_t, true > Gt_1d
Definition: tevaluator.hpp:224
static const int qpts
MFEM_ALWAYS_INLINE void Assemble(DataType &F)
constexpr int DIM
MFEM_ALWAYS_INLINE FieldEvaluator_base(const FE_type &fe, const FiniteElementSpace &fes)
This constructor creates new fespace, not a shallow copy.
Definition: tevaluator.hpp:989
MFEM_ALWAYS_INLINE FieldEvaluator(const FiniteElementSpace &fes, const complex_t *global_data_in, complex_t *global_data_out)
This constructor creates a new fespace, not a shallow copy.
MFEM_ALWAYS_INLINE void EvalSerialized(const typename DataType::vcomplex_t *loc_dofs, DataType &F)
MFEM_ALWAYS_INLINE void AssembleGradGrad(const qpt_layout_t &qpt_layout, const qpt_data_t &qpt_data, const D_layout_t &D_layout, D_data_t &D_data) const
Multi-component assemble of grad-grad element matrices. qpt_layout is (NIP x DIM x DIM x NumComp)...
Definition: tevaluator.hpp:187
static const layout_type layout
Definition: ttensor.hpp:392
MFEM_ALWAYS_INLINE void AssembleGradGrad(const qpt_layout_t &qpt_layout, const qpt_data_t &qpt_data, const D_layout_t &D_layout, D_data_t &D_data) const
Multi-component assemble of grad-grad element matrices. qpt_layout is (TNIP x DIM x DIM x NumComp)...
Definition: tevaluator.hpp:882
This struct implements the input (Eval, EvalSerialized) and output (Assemble, AssembleSerialized) ope...
MFEM_ALWAYS_INLINE void Assemble(const qpt_layout_t &qpt_layout, const qpt_data_t &qpt_data, const M_layout_t &M_layout, M_data_t &M_data) const
Multi-component assemble. qpt_layout is (NIP x NumComp), M_layout is (DOF x DOF x NumComp) ...
Definition: tevaluator.hpp:295
MFEM_ALWAYS_INLINE void Assemble(const qpt_layout_t &qpt_layout, const qpt_data_t &qpt_data, const M_layout_t &M_layout, M_data_t &M_data) const
Multi-component assemble. qpt_layout is (TNIP x NumComp), M_layout is (TDOF x TDOF x NumComp) ...
Definition: tevaluator.hpp:460
TTensor3< qpts, vdim, ne, vcomplex_t, true > val_qpts
ShapeEval_type shapeEval
Definition: tevaluator.hpp:975
TTensor3< DOF, NIP, DIM, real_t > Gt
Definition: tevaluator.hpp:48
MFEM_ALWAYS_INLINE FieldEvaluator(const FieldEvaluator &f, const complex_t *global_data_in, complex_t *global_data_out)
With this constructor, fespace is a shallow copy of f.fespace.
FESpace_t::FE_type FE_type
Definition: tevaluator.hpp:971
static MFEM_ALWAYS_INLINE void Assemble(const vec_layout_t &l, T_type &T, AData_t &D)
static MFEM_ALWAYS_INLINE void EvalSerialized(T_type &T, const typename AData_t::vcomplex_t *loc_dofs, AData_t &D)
std::function< double(const Vector &)> f(double mass_coeff)
Definition: lor_mms.hpp:30
ShapeEvaluator(const FE &fe)
Definition: tevaluator.hpp:957
MFEM_ALWAYS_INLINE void Calc(const dof_layout_t &dof_layout, const dof_data_t &dof_data, const qpt_layout_t &qpt_layout, qpt_data_t &qpt_data) const
Definition: tevaluator.hpp:358
MFEM_ALWAYS_INLINE void Assemble(const qpt_layout_t &qpt_layout, const qpt_data_t &qpt_data, const D_layout_t &D_layout, D_data_t &D_data) const
Definition: tevaluator.hpp:805
static const int OutData
void Add(const DenseMatrix &A, const DenseMatrix &B, double alpha, DenseMatrix &C)
C = A + alpha*B.
Definition: densemat.cpp:2320
TTensor3< dofs, vdim, ne, vcomplex_t, true > val_dofs_t
static const layout_type layout
Definition: ttensor.hpp:414
MFEM_ALWAYS_INLINE void Assemble(const qpt_layout_t &qpt_layout, const qpt_data_t &qpt_data, const M_layout_t &M_layout, M_data_t &M_data) const
Multi-component assemble. qpt_layout is (NIP x NumComp), M_layout is (DOF x DOF x NumComp) ...
Definition: tevaluator.hpp:163
MFEM_ALWAYS_INLINE FieldEvaluator(const FESpace_t &tfes, const ShapeEval_type &shape_eval, const VecLayout_type &vec_layout, const complex_t *global_data_in, complex_t *global_data_out)
With this constructor, fespace is a shallow copy of tfes.
MFEM_ALWAYS_INLINE void AssembleGradGrad(const qpt_layout_t &qpt_layout, const qpt_data_t &qpt_data, const D_layout_t &D_layout, D_data_t &D_data) const
Multi-component assemble of grad-grad element matrices. qpt_layout is (NIP x DIM x DIM x NumComp)...
Definition: tevaluator.hpp:319
MFEM_ALWAYS_INLINE FESpace_type & FESpace()
TTensor3< dofs, vdim, ne, vcomplex_t, true > val_dofs_t
TMatrix< DOF, NIP, real_t, true > Gt_1d
Definition: tevaluator.hpp:626
ShapeEvaluator_base< FE, IR, tensor_prod, real_t > base_class
Definition: tevaluator.hpp:950
static MFEM_ALWAYS_INLINE void Eval(const vec_layout_t &l, T_type &T, AData_t &D)
MFEM_ALWAYS_INLINE void CalcGrad(const dof_layout_t &dof_layout, const dof_data_t &dof_data, const grad_layout_t &grad_layout, grad_data_t &grad_data) const
Multi-component gradient evaluation from DOFs to quadrature points. dof_layout is (TDOF x NumComp) an...
Definition: tevaluator.hpp:714
static MFEM_ALWAYS_INLINE void AssembleSerialized(T_type &T, const AData_t &D, typename AData_t::vcomplex_t *loc_dofs)
TTensor3< NIP, DIM, DOF, real_t, true > G
Definition: tevaluator.hpp:47
MFEM_ALWAYS_INLINE void Calc(const dof_layout_t &dof_layout, const dof_data_t &dof_data, const qpt_layout_t &qpt_layout, qpt_data_t &qpt_data) const
Multi-component shape evaluation from DOFs to quadrature points. dof_layout is (DOF x NumComp) and qp...
Definition: tevaluator.hpp:234
TTensor3< dofs, vdim, ne, vcomplex_t, true > val_dofs_t
MFEM_ALWAYS_INLINE ShapeEval_type & ShapeEval()
static MFEM_ALWAYS_INLINE void Assemble(const vec_layout_t &l, T_type &T, AData_t &D)
MFEM_ALWAYS_INLINE void Eval(int el, DataType &F)
static MFEM_ALWAYS_INLINE void EvalSerialized(T_type &T, const typename AData_t::vcomplex_t *loc_dofs, AData_t &D)
This struct is similar to struct AData, adding separate static data members for the input (InData) an...
static const int dim
MFEM_ALWAYS_INLINE void CalcT(const qpt_layout_t &qpt_layout, const qpt_data_t &qpt_data, const dof_layout_t &dof_layout, dof_data_t &dof_data) const
Definition: tevaluator.hpp:391
TTensor4< qpts, dim, vdim, ne, vcomplex_t > grad_qpts
static MFEM_ALWAYS_INLINE void AssembleSerialized(T_type &T, const AData_t &D, typename AData_t::vcomplex_t *loc_dofs)
static const int dofs
ShapeEvaluator< FE_type, IR, real_t > ShapeEval_type
static const int dim
Definition: tevaluator.hpp:945
MFEM_ALWAYS_INLINE void CalcT(const qpt_layout_t &qpt_layout, const qpt_data_t &qpt_data, const dof_layout_t &dof_layout, dof_data_t &dof_data) const
Definition: tevaluator.hpp:675
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:219
General ShapeEvaluator for any scalar FE type (L2 or H1)
Definition: tevaluator.hpp:940
MFEM_ALWAYS_INLINE void CalcGradT(const grad_layout_t &grad_layout, const grad_data_t &grad_data, const dof_layout_t &dof_layout, dof_data_t &dof_data) const
Multi-component gradient evaluation transpose from quadrature points to DOFs. grad_layout is (TNIP x ...
Definition: tevaluator.hpp:444
const complex_t * data_in
MFEM_ALWAYS_INLINE void GetValues(int el, const val_layout_t &l, val_data_t &vals)
val_layout_t is (qpts x vdim x NE)
complex_t - dof/qpt data type, real_t - ShapeEvaluator (FE basis) data type
Definition: tevaluator.hpp:997
static MFEM_ALWAYS_INLINE void Assemble(const vec_layout_t &l, T_type &T, AData_t &D)
double a
Definition: lissajous.cpp:41
Auxiliary templated struct AData, used by the Eval() and Assemble() methods.
MFEM_ALWAYS_INLINE void SetElement(int el)
TMatrix< NIP, DOF, real_t, true > G_1d
Definition: tevaluator.hpp:345
int dim
Definition: ex24.cpp:53
FieldEvaluator_base< FESpace_t, VecLayout_t, IR, complex_t, real_t > base_class
TTensor3< dofs, vdim, ne, vcomplex_t > val_dofs_t
MFEM_ALWAYS_INLINE void CalcT(const qpt_layout_t &qpt_layout, const qpt_data_t &qpt_data, const dof_layout_t &dof_layout, dof_data_t &dof_data) const
Multi-component shape evaluation transpose from quadrature points to DOFs. qpt_layout is (TNIP x NumC...
Definition: tevaluator.hpp:415
MFEM_ALWAYS_INLINE void GetGradients(int el, const grad_layout_t &l, grad_data_t &grad)
grad_layout_t is (qpts x dim x vdim x NE)
FieldEvaluator< FESpace_t, VecLayout_t, IR, complex_t, real_t > T_type
MFEM_ALWAYS_INLINE void CalcGradT(const grad_layout_t &grad_layout, const grad_data_t &grad_data, const dof_layout_t &dof_layout, dof_data_t &dof_data) const
Multi-component gradient evaluation transpose from quadrature points to DOFs. grad_layout is (NIP x D...
Definition: tevaluator.hpp:278
MFEM_ALWAYS_INLINE void AssembleGradGrad(const qpt_layout_t &qpt_layout, const qpt_data_t &qpt_data, const D_layout_t &D_layout, D_data_t &D_data) const
Multi-component assemble of grad-grad element matrices. qpt_layout is (TNIP x DIM x DIM x NumComp)...
Definition: tevaluator.hpp:554
TMatrix< NIP, DOF, real_t, true > B
Definition: tevaluator.hpp:45
MFEM_ALWAYS_INLINE void Calc(const dof_layout_t &dof_layout, const dof_data_t &dof_data, const qpt_layout_t &qpt_layout, qpt_data_t &qpt_data) const
Definition: tevaluator.hpp:638
static const int qpts
Definition: tevaluator.hpp:946
TTensor3< dofs, vdim, ne, vcomplex_t > val_dofs_t
MFEM_ALWAYS_INLINE void Assemble(const qpt_layout_t &qpt_layout, const qpt_data_t &qpt_data, const M_layout_t &M_layout, M_data_t &M_data) const
Multi-component assemble. qpt_layout is (TNIP x NumComp), M_layout is (TDOF x TDOF x NumComp) ...
Definition: tevaluator.hpp:754
MFEM_HOST_DEVICE void Set(const int height, const int width, const double alpha, const TA *Adata, TB *Bdata)
Compute B = alpha*A, where the matrices A and B are of size height x width with data Adata and Bdata...
Definition: kernels.hpp:326
ShapeEvaluator with 3D tensor-product structure.
Definition: tevaluator.hpp:622
FESpace_t::FE_type FE_type
MFEM_ALWAYS_INLINE void CalcT(const qpt_layout_t &qpt_layout, const qpt_data_t &qpt_data, const dof_layout_t &dof_layout, dof_data_t &dof_data) const
Multi-component shape evaluation transpose from quadrature points to DOFs. qpt_layout is (TNIP x NumC...
Definition: tevaluator.hpp:703
static MFEM_ALWAYS_INLINE void Eval(const vec_layout_t &l, T_type &T, AData_t &D)
data_t data[aligned_size >0?aligned_size:1]
Definition: ttensor.hpp:294
TTensor3< dofs, vdim, ne, vcomplex_t > val_dofs_t
Shape evaluators – values of basis functions on the reference element.
Definition: tevaluator.hpp:33
MFEM_ALWAYS_INLINE void Assemble(int el, DataType &F)
BData< InData, OutData, impl_traits_t > DataType
MFEM_ALWAYS_INLINE FieldEvaluator_base(const FESpace_t &tfes, const ShapeEval_type &shape_eval, const VecLayout_t &vec_layout)
With this constructor, fespace is a shallow copy.
Definition: tevaluator.hpp:980
InOutData
Enumeration for the data type used by the Eval() and Assemble() methods. The types can be obtained by...
TMatrix< NIP, DOF, real_t, true > G_1d
Definition: tevaluator.hpp:625
Field evaluators – values of a given global FE grid function This is roughly speaking a templated ve...
Definition: tevaluator.hpp:968
static const int vdim
static MFEM_ALWAYS_INLINE void EvalSerialized(T_type &T, const typename AData_t::vcomplex_t *loc_dofs, AData_t &D)
TElementMatrix< InData, OutData, impl_traits_t > ElementMatrix