MFEM v4.7.0
Finite element discretization library
Loading...
Searching...
No Matches
eval.hpp
Go to the documentation of this file.
1// Copyright (c) 2010-2024, 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
19#include "../kernels.hpp"
20
21namespace mfem
22{
23
24namespace internal
25{
26
27namespace quadrature_interpolator
28{
29
30template<QVectorLayout Q_LAYOUT>
31static void Values1D(const int NE,
32 const real_t *b_,
33 const real_t *x_,
34 real_t *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 real_t 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.
64template<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>
67static void Values2D(const int NE,
68 const real_t *b_,
69 const real_t *x_,
70 real_t *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 real_t sB[MQ1*MD1];
98 MFEM_SHARED real_t sm0[NBZ][MDQ*MDQ];
99 MFEM_SHARED real_t 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 real_t 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.
128template<QVectorLayout Q_LAYOUT,
129 int T_VDIM = 0, int T_D1D = 0, int T_Q1D = 0>
130static void Values3D(const int NE,
131 const real_t *b_,
132 const real_t *x_,
133 real_t *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 real_t sB[MQ1*MD1];
158 MFEM_SHARED real_t sm0[MDQ*MDQ*MDQ];
159 MFEM_SHARED real_t 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 real_t 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
real_t b
Definition lissajous.cpp:42
real_t u(const Vector &xvec)
Definition lor_mms.hpp:22
DeviceTensor< 2, real_t > DeviceMatrix
Definition dtensor.hpp:143
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
void forall_2D_batch(int N, int X, int Y, int BZ, lambda &&body)
Definition forall.hpp:769
DeviceTensor< 2, const real_t > ConstDeviceMatrix
Definition dtensor.hpp:144
void forall_3D(int N, int X, int Y, int Z, lambda &&body)
Definition forall.hpp:775
QVectorLayout
Type describing possible layouts for Q-vectors.
Definition fespace.hpp:53
@ byNODES
NQPT x VDIM x NE (values) / NQPT x VDIM x DIM x NE (grads)
@ byVDIM
VDIM x NQPT x NE (values) / VDIM x DIM x NQPT x NE (grads)
float real_t
Definition config.hpp:43
DeviceTensor< 3, real_t > DeviceCube
Definition dtensor.hpp:146
void forall(int N, lambda &&body)
Definition forall.hpp:754