MFEM  v4.6.0
Finite element discretization library
eval.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 // Internal header, included only by .cpp files.
13 // Template function implementations.
14 
15 #include "../quadinterpolator.hpp"
16 #include "../../general/forall.hpp"
17 #include "../../linalg/dtensor.hpp"
18 #include "../../linalg/kernels.hpp"
19 #include "../kernels.hpp"
20 
21 namespace mfem
22 {
23 
24 namespace internal
25 {
26 
27 namespace quadrature_interpolator
28 {
29 
30 template<QVectorLayout Q_LAYOUT>
31 static void Values1D(const int NE,
32  const double *b_,
33  const double *x_,
34  double *y_,
35  const int vdim,
36  const int d1d,
37  const int q1d)
38 {
39  const auto b = Reshape(b_, q1d, d1d);
40  const auto x = Reshape(x_, d1d, vdim, NE);
41  auto y = Q_LAYOUT == QVectorLayout::byNODES ?
42  Reshape(y_, q1d, vdim, NE):
43  Reshape(y_, vdim, q1d, NE);
44 
45  mfem::forall(NE, [=] MFEM_HOST_DEVICE (int e)
46  {
47  for (int c = 0; c < vdim; c++)
48  {
49  for (int q = 0; q < q1d; q++)
50  {
51  double u = 0.0;
52  for (int d = 0; d < d1d; d++)
53  {
54  u += b(q, d) * x(d, c, e);
55  }
56  if (Q_LAYOUT == QVectorLayout::byVDIM) { y(c, q, e) = u; }
57  if (Q_LAYOUT == QVectorLayout::byNODES) { y(q, c, e) = u; }
58  }
59  }
60  });
61 }
62 
63 // Template compute kernel for Values in 2D: tensor product version.
64 template<QVectorLayout Q_LAYOUT,
65  int T_VDIM = 0, int T_D1D = 0, int T_Q1D = 0,
66  int T_NBZ = 1, int MAX_D1D = 0, int MAX_Q1D = 0>
67 static void Values2D(const int NE,
68  const double *b_,
69  const double *x_,
70  double *y_,
71  const int vdim = 0,
72  const int d1d = 0,
73  const int q1d = 0)
74 {
75  static constexpr int NBZ = T_NBZ ? T_NBZ : 1;
76 
77  const int D1D = T_D1D ? T_D1D : d1d;
78  const int Q1D = T_Q1D ? T_Q1D : q1d;
79  const int VDIM = T_VDIM ? T_VDIM : vdim;
80 
81  const auto b = Reshape(b_, Q1D, D1D);
82  const auto x = Reshape(x_, D1D, D1D, VDIM, NE);
83  auto y = Q_LAYOUT == QVectorLayout::byNODES ?
84  Reshape(y_, Q1D, Q1D, VDIM, NE):
85  Reshape(y_, VDIM, Q1D, Q1D, NE);
86 
87  mfem::forall_2D_batch(NE, Q1D, Q1D, NBZ, [=] MFEM_HOST_DEVICE (int e)
88  {
89  const int D1D = T_D1D ? T_D1D : d1d;
90  const int Q1D = T_Q1D ? T_Q1D : q1d;
91  const int VDIM = T_VDIM ? T_VDIM : vdim;
92  constexpr int MQ1 = T_Q1D ? T_Q1D : DofQuadLimits::MAX_Q1D;
93  constexpr int MD1 = T_D1D ? T_D1D : DofQuadLimits::MAX_D1D;
94  constexpr int MDQ = (MQ1 > MD1) ? MQ1 : MD1;
95  const int tidz = MFEM_THREAD_ID(z);
96 
97  MFEM_SHARED double sB[MQ1*MD1];
98  MFEM_SHARED double sm0[NBZ][MDQ*MDQ];
99  MFEM_SHARED double sm1[NBZ][MDQ*MDQ];
100 
101  kernels::internal::LoadB<MD1,MQ1>(D1D,Q1D,b,sB);
102 
103  ConstDeviceMatrix B(sB, D1D,Q1D);
104  DeviceMatrix DD(sm0[tidz], MD1, MD1);
105  DeviceMatrix DQ(sm1[tidz], MD1, MQ1);
106  DeviceMatrix QQ(sm0[tidz], MQ1, MQ1);
107 
108  for (int c = 0; c < VDIM; c++)
109  {
110  kernels::internal::LoadX(e,D1D,c,x,DD);
111  kernels::internal::EvalX(D1D,Q1D,B,DD,DQ);
112  kernels::internal::EvalY(D1D,Q1D,B,DQ,QQ);
113  MFEM_FOREACH_THREAD(qy,y,Q1D)
114  {
115  MFEM_FOREACH_THREAD(qx,x,Q1D)
116  {
117  double u = QQ(qx,qy);
118  if (Q_LAYOUT == QVectorLayout::byVDIM) { y(c,qx,qy,e) = u; }
119  if (Q_LAYOUT == QVectorLayout::byNODES) { y(qx,qy,c,e) = u; }
120  }
121  }
122  MFEM_SYNC_THREAD;
123  }
124  });
125 }
126 
127 // Template compute kernel for Values in 3D: tensor product version.
128 template<QVectorLayout Q_LAYOUT,
129  int T_VDIM = 0, int T_D1D = 0, int T_Q1D = 0>
130 static void Values3D(const int NE,
131  const double *b_,
132  const double *x_,
133  double *y_,
134  const int vdim = 0,
135  const int d1d = 0,
136  const int q1d = 0)
137 {
138  const int D1D = T_D1D ? T_D1D : d1d;
139  const int Q1D = T_Q1D ? T_Q1D : q1d;
140  const int VDIM = T_VDIM ? T_VDIM : vdim;
141 
142  const auto b = Reshape(b_, Q1D, D1D);
143  const auto x = Reshape(x_, D1D, D1D, D1D, VDIM, NE);
144  auto y = Q_LAYOUT == QVectorLayout:: byNODES ?
145  Reshape(y_, Q1D, Q1D, Q1D, VDIM, NE):
146  Reshape(y_, VDIM, Q1D, Q1D, Q1D, NE);
147 
148  mfem::forall_3D(NE, Q1D, Q1D, Q1D, [=] MFEM_HOST_DEVICE (int e)
149  {
150  const int D1D = T_D1D ? T_D1D : d1d;
151  const int Q1D = T_Q1D ? T_Q1D : q1d;
152  const int VDIM = T_VDIM ? T_VDIM : vdim;
153  constexpr int MQ1 = T_Q1D ? T_Q1D : DofQuadLimits::MAX_INTERP_1D;
154  constexpr int MD1 = T_D1D ? T_D1D : DofQuadLimits::MAX_INTERP_1D;
155  constexpr int MDQ = (MQ1 > MD1) ? MQ1 : MD1;
156 
157  MFEM_SHARED double sB[MQ1*MD1];
158  MFEM_SHARED double sm0[MDQ*MDQ*MDQ];
159  MFEM_SHARED double sm1[MDQ*MDQ*MDQ];
160 
161  kernels::internal::LoadB<MD1,MQ1>(D1D,Q1D,b,sB);
162 
163  ConstDeviceMatrix B(sB, D1D,Q1D);
164  DeviceCube DDD(sm0, MD1,MD1,MD1);
165  DeviceCube DDQ(sm1, MD1,MD1,MQ1);
166  DeviceCube DQQ(sm0, MD1,MQ1,MQ1);
167  DeviceCube QQQ(sm1, MQ1,MQ1,MQ1);
168 
169  for (int c = 0; c < VDIM; c++)
170  {
171  kernels::internal::LoadX(e,D1D,c,x,DDD);
172  kernels::internal::EvalX(D1D,Q1D,B,DDD,DDQ);
173  kernels::internal::EvalY(D1D,Q1D,B,DDQ,DQQ);
174  kernels::internal::EvalZ(D1D,Q1D,B,DQQ,QQQ);
175  MFEM_FOREACH_THREAD(qz,z,Q1D)
176  {
177  MFEM_FOREACH_THREAD(qy,y,Q1D)
178  {
179  MFEM_FOREACH_THREAD(qx,x,Q1D)
180  {
181  const double u = QQQ(qz,qy,qx);
182  if (Q_LAYOUT == QVectorLayout::byVDIM) { y(c,qx,qy,qz,e) = u; }
183  if (Q_LAYOUT == QVectorLayout::byNODES) { y(qx,qy,qz,c,e) = u; }
184  }
185  }
186  }
187  MFEM_SYNC_THREAD;
188  }
189  });
190 }
191 
192 } // namespace quadrature_interpolator
193 
194 } // namespace internal
195 
196 } // namespace mfem
void forall_3D(int N, int X, int Y, int Z, lambda &&body)
Definition: forall.hpp:763
DeviceTensor< 2, const double > ConstDeviceMatrix
Definition: dtensor.hpp:144
DeviceTensor< 2, double > DeviceMatrix
Definition: dtensor.hpp:143
double b
Definition: lissajous.cpp:42
void forall_2D_batch(int N, int X, int Y, int BZ, lambda &&body)
Definition: forall.hpp:757
DeviceTensor< 3, double > DeviceCube
Definition: dtensor.hpp:146
void forall(int N, lambda &&body)
Definition: forall.hpp:742
QVectorLayout
Type describing possible layouts for Q-vectors.
Definition: fespace.hpp:52
NQPT x VDIM x NE (values) / NQPT x VDIM x DIM x NE (grads)
double u(const Vector &xvec)
Definition: lor_mms.hpp:22
MFEM_HOST_DEVICE DeviceTensor< sizeof...(Dims), T > Reshape(T *ptr, Dims... dims)
Wrap a pointer as a DeviceTensor with automatically deduced template parameters.
Definition: dtensor.hpp:131
VDIM x NQPT x NE (values) / VDIM x DIM x NQPT x NE (grads)