MFEM  v4.5.2
Finite element discretization library
grad.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 "../../fem/kernels.hpp"
19 #include "../../linalg/kernels.hpp"
20 
21 namespace mfem
22 {
23 
24 namespace internal
25 {
26 
27 namespace quadrature_interpolator
28 {
29 
30 // Template compute kernel for derivatives in 2D: tensor product version.
31 template<QVectorLayout Q_LAYOUT, bool GRAD_PHYS,
32  int T_VDIM = 0, int T_D1D = 0, int T_Q1D = 0,
33  int T_NBZ = 1, int MAX_D1D = 0, int MAX_Q1D = 0>
34 static void Derivatives2D(const int NE,
35  const double *b_,
36  const double *g_,
37  const double *j_,
38  const double *x_,
39  double *y_,
40  const int vdim = 0,
41  const int d1d = 0,
42  const int q1d = 0)
43 {
44  const int D1D = T_D1D ? T_D1D : d1d;
45  const int Q1D = T_Q1D ? T_Q1D : q1d;
46  const int VDIM = T_VDIM ? T_VDIM : vdim;
47  static constexpr int NBZ = T_NBZ ? T_NBZ : 1;
48 
49  const auto b = Reshape(b_, Q1D, D1D);
50  const auto g = Reshape(g_, Q1D, D1D);
51  const auto j = Reshape(j_, Q1D, Q1D, 2, 2, NE);
52  const auto x = Reshape(x_, D1D, D1D, VDIM, NE);
53  auto y = Q_LAYOUT == QVectorLayout:: byNODES ?
54  Reshape(y_, Q1D, Q1D, VDIM, 2, NE):
55  Reshape(y_, VDIM, 2, Q1D, Q1D, NE);
56 
57  MFEM_FORALL_2D(e, NE, Q1D, Q1D, NBZ,
58  {
59  const int D1D = T_D1D ? T_D1D : d1d;
60  const int Q1D = T_Q1D ? T_Q1D : q1d;
61  const int VDIM = T_VDIM ? T_VDIM : vdim;
62  constexpr int MQ1 = T_Q1D ? T_Q1D : MAX_Q1D;
63  constexpr int MD1 = T_D1D ? T_D1D : MAX_D1D;
64 
65  const int tidz = MFEM_THREAD_ID(z);
66  MFEM_SHARED double BG[2][MQ1*MD1];
67  kernels::internal::LoadBG<MD1,MQ1>(D1D,Q1D,b,g,BG);
68  DeviceMatrix B(BG[0], D1D, Q1D);
69  DeviceMatrix G(BG[1], D1D, Q1D);
70 
71  MFEM_SHARED double XY[NBZ][MD1*MD1];
72  DeviceTensor<2> X((double*)(XY+tidz), D1D, D1D);
73 
74  MFEM_SHARED double s_DQ[2][NBZ][MD1*MQ1];
75  DeviceTensor<2> DQ0(s_DQ[0][tidz], D1D, Q1D);
76  DeviceTensor<2> DQ1(s_DQ[1][tidz], D1D, Q1D);
77 
78  for (int c = 0; c < VDIM; ++c)
79  {
80  kernels::internal::LoadX<MD1,NBZ>(e,D1D,c,x,XY);
81  MFEM_FOREACH_THREAD(dy,y,D1D)
82  {
83  MFEM_FOREACH_THREAD(qx,x,Q1D)
84  {
85  double u = 0.0;
86  double v = 0.0;
87  for (int dx = 0; dx < D1D; ++dx)
88  {
89  const double input = X(dx,dy);
90  u += input * B(dx,qx);
91  v += input * G(dx,qx);
92  }
93  DQ0(dy,qx) = u;
94  DQ1(dy,qx) = v;
95  }
96  }
97  MFEM_SYNC_THREAD;
98  MFEM_FOREACH_THREAD(qy,y,Q1D)
99  {
100  MFEM_FOREACH_THREAD(qx,x,Q1D)
101  {
102  double u = 0.0;
103  double v = 0.0;
104  for (int dy = 0; dy < D1D; ++dy)
105  {
106  u += DQ1(dy,qx) * B(dy,qy);
107  v += DQ0(dy,qx) * G(dy,qy);
108  }
109  if (GRAD_PHYS)
110  {
111  double Jloc[4], Jinv[4];
112  Jloc[0] = j(qx,qy,0,0,e);
113  Jloc[1] = j(qx,qy,1,0,e);
114  Jloc[2] = j(qx,qy,0,1,e);
115  Jloc[3] = j(qx,qy,1,1,e);
116  kernels::CalcInverse<2>(Jloc, Jinv);
117  const double U = Jinv[0]*u + Jinv[1]*v;
118  const double V = Jinv[2]*u + Jinv[3]*v;
119  u = U; v = V;
120  }
121  if (Q_LAYOUT == QVectorLayout::byVDIM)
122  {
123  y(c,0,qx,qy,e) = u;
124  y(c,1,qx,qy,e) = v;
125  }
126  if (Q_LAYOUT == QVectorLayout::byNODES)
127  {
128  y(qx,qy,c,0,e) = u;
129  y(qx,qy,c,1,e) = v;
130  }
131  }
132  }
133  MFEM_SYNC_THREAD;
134  }
135  });
136 }
137 
138 // Template compute kernel for derivatives in 3D: tensor product version.
139 template<QVectorLayout Q_LAYOUT, bool GRAD_PHYS,
140  int T_VDIM = 0, int T_D1D = 0, int T_Q1D = 0,
141  int MAX_D1D = 0, int MAX_Q1D = 0>
142 static void Derivatives3D(const int NE,
143  const double *b_,
144  const double *g_,
145  const double *j_,
146  const double *x_,
147  double *y_,
148  const int vdim = 0,
149  const int d1d = 0,
150  const int q1d = 0)
151 {
152  const int D1D = T_D1D ? T_D1D : d1d;
153  const int Q1D = T_Q1D ? T_Q1D : q1d;
154  const int VDIM = T_VDIM ? T_VDIM : vdim;
155 
156  const auto b = Reshape(b_, Q1D, D1D);
157  const auto g = Reshape(g_, Q1D, D1D);
158  const auto j = Reshape(j_, Q1D, Q1D, Q1D, 3, 3, NE);
159  const auto x = Reshape(x_, D1D, D1D, D1D, VDIM, NE);
160  auto y = Q_LAYOUT == QVectorLayout:: byNODES ?
161  Reshape(y_, Q1D, Q1D, Q1D, VDIM, 3, NE):
162  Reshape(y_, VDIM, 3, Q1D, Q1D, Q1D, NE);
163 
164  MFEM_FORALL_3D(e, NE, Q1D, Q1D, Q1D,
165  {
166  const int D1D = T_D1D ? T_D1D : d1d;
167  const int Q1D = T_Q1D ? T_Q1D : q1d;
168  const int VDIM = T_VDIM ? T_VDIM : vdim;
169  constexpr int MQ1 = T_Q1D ? T_Q1D : MAX_Q1D;
170  constexpr int MD1 = T_D1D ? T_D1D : MAX_D1D;
171 
172  MFEM_SHARED double BG[2][MQ1*MD1];
173  kernels::internal::LoadBG<MD1,MQ1>(D1D,Q1D,b,g,BG);
174  DeviceMatrix B(BG[0], D1D, Q1D);
175  DeviceMatrix G(BG[1], D1D, Q1D);
176 
177  MFEM_SHARED double sm0[3][MQ1*MQ1*MQ1];
178  MFEM_SHARED double sm1[3][MQ1*MQ1*MQ1];
179  DeviceTensor<3> X(sm0[2], D1D, D1D, D1D);
180  DeviceTensor<3> DDQ0(sm0[0], D1D, D1D, Q1D);
181  DeviceTensor<3> DDQ1(sm0[1], D1D, D1D, Q1D);
182  DeviceTensor<3> DQQ0(sm1[0], D1D, Q1D, Q1D);
183  DeviceTensor<3> DQQ1(sm1[1], D1D, Q1D, Q1D);
184  DeviceTensor<3> DQQ2(sm1[2], D1D, Q1D, Q1D);
185 
186  for (int c = 0; c < VDIM; ++c)
187  {
188  kernels::internal::LoadX(e,D1D,c,x,X);
189  MFEM_FOREACH_THREAD(dz,z,D1D)
190  {
191  MFEM_FOREACH_THREAD(dy,y,D1D)
192  {
193  MFEM_FOREACH_THREAD(qx,x,Q1D)
194  {
195  double u = 0.0;
196  double v = 0.0;
197  for (int dx = 0; dx < D1D; ++dx)
198  {
199  const double input = X(dx,dy,dz);
200  u += input * B(dx,qx);
201  v += input * G(dx,qx);
202  }
203  DDQ0(dz,dy,qx) = u;
204  DDQ1(dz,dy,qx) = v;
205  }
206  }
207  }
208  MFEM_SYNC_THREAD;
209  MFEM_FOREACH_THREAD(dz,z,D1D)
210  {
211  MFEM_FOREACH_THREAD(qy,y,Q1D)
212  {
213  MFEM_FOREACH_THREAD(qx,x,Q1D)
214  {
215  double u = 0.0;
216  double v = 0.0;
217  double w = 0.0;
218  for (int dy = 0; dy < D1D; ++dy)
219  {
220  u += DDQ1(dz,dy,qx) * B(dy,qy);
221  v += DDQ0(dz,dy,qx) * G(dy,qy);
222  w += DDQ0(dz,dy,qx) * B(dy,qy);
223  }
224  DQQ0(dz,qy,qx) = u;
225  DQQ1(dz,qy,qx) = v;
226  DQQ2(dz,qy,qx) = w;
227  }
228  }
229  }
230  MFEM_SYNC_THREAD;
231  MFEM_FOREACH_THREAD(qz,z,Q1D)
232  {
233  MFEM_FOREACH_THREAD(qy,y,Q1D)
234  {
235  MFEM_FOREACH_THREAD(qx,x,Q1D)
236  {
237  double u = 0.0;
238  double v = 0.0;
239  double w = 0.0;
240  for (int dz = 0; dz < D1D; ++dz)
241  {
242  u += DQQ0(dz,qy,qx) * B(dz,qz);
243  v += DQQ1(dz,qy,qx) * B(dz,qz);
244  w += DQQ2(dz,qy,qx) * G(dz,qz);
245  }
246  if (GRAD_PHYS)
247  {
248  double Jloc[9], Jinv[9];
249  for (int col = 0; col < 3; col++)
250  {
251  for (int row = 0; row < 3; row++)
252  {
253  Jloc[row+3*col] = j(qx,qy,qz,row,col,e);
254  }
255  }
256  kernels::CalcInverse<3>(Jloc, Jinv);
257  const double U = Jinv[0]*u + Jinv[1]*v + Jinv[2]*w;
258  const double V = Jinv[3]*u + Jinv[4]*v + Jinv[5]*w;
259  const double W = Jinv[6]*u + Jinv[7]*v + Jinv[8]*w;
260  u = U; v = V; w = W;
261  }
262  if (Q_LAYOUT == QVectorLayout::byVDIM)
263  {
264  y(c,0,qx,qy,qz,e) = u;
265  y(c,1,qx,qy,qz,e) = v;
266  y(c,2,qx,qy,qz,e) = w;
267  }
268  if (Q_LAYOUT == QVectorLayout::byNODES)
269  {
270  y(qx,qy,qz,c,0,e) = u;
271  y(qx,qy,qz,c,1,e) = v;
272  y(qx,qy,qz,c,2,e) = w;
273  }
274  }
275  }
276  }
277  MFEM_SYNC_THREAD;
278  }
279  });
280 }
281 
282 } // namespace quadrature_interpolator
283 
284 } // namespace internal
285 
286 } // namespace mfem
mfem::DeviceMatrix
DeviceTensor< 2, double > DeviceMatrix
Definition: dtensor.hpp:143
mfem::MAX_Q1D
const int MAX_Q1D
Definition: forall.hpp:29
b
double b
Definition: lissajous.cpp:42
mfem::QVectorLayout
QVectorLayout
Type describing possible layouts for Q-vectors.
Definition: fespace.hpp:52
mfem::MAX_D1D
const int MAX_D1D
Definition: forall.hpp:28
mfem::QVectorLayout::byNODES
NQPT x VDIM x NE (values) / NQPT x VDIM x DIM x NE (grads)
mfem::u
double u(const Vector &xvec)
Definition: lor_mms.hpp:24
mfem::Reshape
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
mfem::QVectorLayout::byVDIM
VDIM x NQPT x NE (values) / VDIM x DIM x NQPT x NE (grads)
Generated on Thu Mar 23 2023 16:28:21 for MFEM by   doxygen 1.8.14