MFEM  v4.6.0
Finite element discretization library
eval_by_nodes.cpp
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 #include "../quadinterpolator.hpp"
13 #include "dispatch.hpp"
14 #include "eval.hpp"
15 
16 namespace mfem
17 {
18 
19 namespace internal
20 {
21 
22 namespace quadrature_interpolator
23 {
24 
25 // Tensor-product evaluation of quadrature point values: dispatch function.
26 // Instantiation for the case QVectorLayout::byNODES.
27 template<>
28 void TensorValues<QVectorLayout::byNODES>(const int NE,
29  const int vdim,
30  const DofToQuad &maps,
31  const Vector &e_vec,
32  Vector &q_val)
33 {
34  if (NE == 0) { return; }
35  const int dim = maps.FE->GetDim();
36  const int D1D = maps.ndof;
37  const int Q1D = maps.nqpt;
38  const double *B = maps.B.Read();
39  const double *X = e_vec.Read();
40  double *Y = q_val.Write();
41 
43 
44  const int id = (vdim<<8) | (D1D<<4) | Q1D;
45 
46  if (dim == 1)
47  {
48  MFEM_VERIFY(D1D <= DeviceDofQuadLimits::Get().MAX_D1D,
49  "Orders higher than " << DeviceDofQuadLimits::Get().MAX_D1D-1
50  << " are not supported!");
51  MFEM_VERIFY(Q1D <= DeviceDofQuadLimits::Get().MAX_Q1D,
52  "Quadrature rules with more than "
53  << DeviceDofQuadLimits::Get().MAX_Q1D << " 1D points are not supported!");
54  Values1D<L>(NE, B, X, Y, vdim, D1D, Q1D);
55  return;
56  }
57  if (dim == 2)
58  {
59  switch (id)
60  {
61  case 0x133: return Values2D<L,1,3,3>(NE,B,X,Y);
62  case 0x124: return Values2D<L,1,2,4>(NE,B,X,Y);
63  case 0x132: return Values2D<L,1,3,2>(NE,B,X,Y);
64  case 0x134: return Values2D<L,1,3,4>(NE,B,X,Y);
65  case 0x143: return Values2D<L,1,4,3>(NE,B,X,Y);
66  case 0x144: return Values2D<L,1,4,4>(NE,B,X,Y);
67 
68  case 0x222: return Values2D<L,2,2,2>(NE,B,X,Y);
69  case 0x223: return Values2D<L,2,2,3>(NE,B,X,Y);
70  case 0x224: return Values2D<L,2,2,4>(NE,B,X,Y);
71  case 0x225: return Values2D<L,2,2,5>(NE,B,X,Y);
72  case 0x226: return Values2D<L,2,2,6>(NE,B,X,Y);
73 
74  case 0x233: return Values2D<L,2,3,3>(NE,B,X,Y);
75  case 0x234: return Values2D<L,2,3,4>(NE,B,X,Y);
76  case 0x236: return Values2D<L,2,3,6>(NE,B,X,Y);
77 
78  case 0x243: return Values2D<L,2,4,3>(NE,B,X,Y);
79  case 0x244: return Values2D<L,2,4,4>(NE,B,X,Y);
80  case 0x245: return Values2D<L,2,4,5>(NE,B,X,Y);
81  case 0x246: return Values2D<L,2,4,6>(NE,B,X,Y);
82  case 0x247: return Values2D<L,2,4,7>(NE,B,X,Y);
83 
84  case 0x256: return Values2D<L,2,5,6>(NE,B,X,Y);
85 
86  default:
87  {
88  const int MD = DeviceDofQuadLimits::Get().MAX_D1D;
89  const int MQ = DeviceDofQuadLimits::Get().MAX_Q1D;
90  MFEM_VERIFY(D1D <= MD, "Orders higher than " << MD-1
91  << " are not supported!");
92  MFEM_VERIFY(Q1D <= MQ, "Quadrature rules with more than "
93  << MQ << " 1D points are not supported!");
94  Values2D<L>(NE,B,X,Y,vdim,D1D,Q1D);
95  return;
96  }
97  }
98  }
99  if (dim == 3)
100  {
101  switch (id)
102  {
103  case 0x124: return Values3D<L,1,2,4>(NE,B,X,Y);
104  case 0x133: return Values3D<L,1,3,3>(NE,B,X,Y);
105  case 0x134: return Values3D<L,1,3,4>(NE,B,X,Y);
106  case 0x136: return Values3D<L,1,3,6>(NE,B,X,Y);
107  case 0x143: return Values3D<L,1,4,3>(NE,B,X,Y);
108  case 0x144: return Values3D<L,1,4,4>(NE,B,X,Y);
109  case 0x148: return Values3D<L,1,4,8>(NE,B,X,Y);
110 
111  case 0x222: return Values3D<L,2,2,2>(NE,B,X,Y);
112  case 0x223: return Values3D<L,2,2,3>(NE,B,X,Y);
113  case 0x234: return Values3D<L,2,3,4>(NE,B,X,Y);
114 
115  case 0x323: return Values3D<L,3,2,3>(NE,B,X,Y);
116  case 0x324: return Values3D<L,3,2,4>(NE,B,X,Y);
117  case 0x325: return Values3D<L,3,2,5>(NE,B,X,Y);
118  case 0x326: return Values3D<L,3,2,6>(NE,B,X,Y);
119 
120  case 0x333: return Values3D<L,3,3,3>(NE,B,X,Y);
121  case 0x334: return Values3D<L,3,3,4>(NE,B,X,Y);
122  case 0x335: return Values3D<L,3,3,5>(NE,B,X,Y);
123  case 0x336: return Values3D<L,3,3,6>(NE,B,X,Y);
124 
125  case 0x343: return Values3D<L,3,4,3>(NE,B,X,Y);
126  case 0x344: return Values3D<L,3,4,4>(NE,B,X,Y);
127  case 0x346: return Values3D<L,3,4,6>(NE,B,X,Y);
128  case 0x347: return Values3D<L,3,4,7>(NE,B,X,Y);
129  case 0x348: return Values3D<L,3,4,8>(NE,B,X,Y);
130 
131  default:
132  {
133  const int MD = DeviceDofQuadLimits::Get().MAX_INTERP_1D;
134  const int MQ = DeviceDofQuadLimits::Get().MAX_INTERP_1D;
135  MFEM_VERIFY(D1D <= MD, "Orders higher than " << MD-1
136  << " are not supported!");
137  MFEM_VERIFY(Q1D <= MQ, "Quadrature rules with more than "
138  << MQ << " 1D points are not supported!");
139  Values3D<L>(NE,B,X,Y,vdim,D1D,Q1D);
140  return;
141  }
142  }
143  }
144  mfem::out << "Unknown kernel 0x" << std::hex << id << std::endl;
145  MFEM_ABORT("Kernel not supported yet");
146 }
147 
148 } // namespace quadrature_interpolator
149 
150 } // namespace internal
151 
152 } // namespace mfem
int MAX_Q1D
Maximum number of 1D quadrature points.
Definition: forall.hpp:113
int MAX_D1D
Maximum number of 1D nodal points.
Definition: forall.hpp:112
OutStream out(std::cout)
Global stream used by the library for standard output. Initially it uses the same std::streambuf as s...
Definition: globals.hpp:66
static const DeviceDofQuadLimits & Get()
Return a const reference to the DeviceDofQuadLimits singleton.
Definition: forall.hpp:122
QVectorLayout
Type describing possible layouts for Q-vectors.
Definition: fespace.hpp:52
int dim
Definition: ex24.cpp:53
int MAX_INTERP_1D
Maximum number of points for use in QuadratureInterpolator.
Definition: forall.hpp:118
NQPT x VDIM x NE (values) / NQPT x VDIM x DIM x NE (grads)