MFEM  v4.1.0
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
quadinterpolator.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2020, 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 "../general/forall.hpp"
14 #include "../linalg/dtensor.hpp"
15 #include "../linalg/kernels.hpp"
16 
17 namespace mfem
18 {
19 
20 QuadratureInterpolator::QuadratureInterpolator(const FiniteElementSpace &fes,
21  const IntegrationRule &ir)
22 {
23  fespace = &fes;
24  qspace = NULL;
25  IntRule = &ir;
26  q_layout = QVectorLayout::byNODES;
27  use_tensor_products = true; // not implemented yet (not used)
28 
29  if (fespace->GetNE() == 0) { return; }
30  const FiniteElement *fe = fespace->GetFE(0);
31  MFEM_VERIFY(dynamic_cast<const ScalarFiniteElement*>(fe) != NULL,
32  "Only scalar finite elements are supported");
33 }
34 
35 QuadratureInterpolator::QuadratureInterpolator(const FiniteElementSpace &fes,
36  const QuadratureSpace &qs)
37 {
38  fespace = &fes;
39  qspace = &qs;
40  IntRule = NULL;
41  q_layout = QVectorLayout::byNODES;
42  use_tensor_products = true; // not implemented yet (not used)
43 
44  if (fespace->GetNE() == 0) { return; }
45  const FiniteElement *fe = fespace->GetFE(0);
46  MFEM_VERIFY(dynamic_cast<const ScalarFiniteElement*>(fe) != NULL,
47  "Only scalar finite elements are supported");
48 }
49 
50 template<const int T_VDIM, const int T_ND, const int T_NQ>
51 void QuadratureInterpolator::Eval2D(
52  const int NE,
53  const int vdim,
54  const DofToQuad &maps,
55  const Vector &e_vec,
56  Vector &q_val,
57  Vector &q_der,
58  Vector &q_det,
59  const int eval_flags)
60 {
61  const int nd = maps.ndof;
62  const int nq = maps.nqpt;
63  const int ND = T_ND ? T_ND : nd;
64  const int NQ = T_NQ ? T_NQ : nq;
65  const int VDIM = T_VDIM ? T_VDIM : vdim;
66  MFEM_VERIFY(ND <= MAX_ND2D, "");
67  MFEM_VERIFY(NQ <= MAX_NQ2D, "");
68  MFEM_VERIFY(VDIM == 2 || !(eval_flags & DETERMINANTS), "");
69  auto B = Reshape(maps.B.Read(), NQ, ND);
70  auto G = Reshape(maps.G.Read(), NQ, 2, ND);
71  auto E = Reshape(e_vec.Read(), ND, VDIM, NE);
72  auto val = Reshape(q_val.Write(), NQ, VDIM, NE);
73  auto der = Reshape(q_der.Write(), NQ, VDIM, 2, NE);
74  auto det = Reshape(q_det.Write(), NQ, NE);
75  MFEM_FORALL(e, NE,
76  {
77  const int ND = T_ND ? T_ND : nd;
78  const int NQ = T_NQ ? T_NQ : nq;
79  const int VDIM = T_VDIM ? T_VDIM : vdim;
80  constexpr int max_ND = T_ND ? T_ND : MAX_ND2D;
81  constexpr int max_VDIM = T_VDIM ? T_VDIM : MAX_VDIM2D;
82  double s_E[max_VDIM*max_ND];
83  for (int d = 0; d < ND; d++)
84  {
85  for (int c = 0; c < VDIM; c++)
86  {
87  s_E[c+d*VDIM] = E(d,c,e);
88  }
89  }
90  for (int q = 0; q < NQ; ++q)
91  {
92  if (eval_flags & VALUES)
93  {
94  double ed[max_VDIM];
95  for (int c = 0; c < VDIM; c++) { ed[c] = 0.0; }
96  for (int d = 0; d < ND; ++d)
97  {
98  const double b = B(q,d);
99  for (int c = 0; c < VDIM; c++) { ed[c] += b*s_E[c+d*VDIM]; }
100  }
101  for (int c = 0; c < VDIM; c++) { val(q,c,e) = ed[c]; }
102  }
103  if ((eval_flags & DERIVATIVES) || (eval_flags & DETERMINANTS))
104  {
105  // use MAX_VDIM2D to avoid "subscript out of range" warnings
106  double D[MAX_VDIM2D*2];
107  for (int i = 0; i < 2*VDIM; i++) { D[i] = 0.0; }
108  for (int d = 0; d < ND; ++d)
109  {
110  const double wx = G(q,0,d);
111  const double wy = G(q,1,d);
112  for (int c = 0; c < VDIM; c++)
113  {
114  double s_e = s_E[c+d*VDIM];
115  D[c+VDIM*0] += s_e * wx;
116  D[c+VDIM*1] += s_e * wy;
117  }
118  }
119  if (eval_flags & DERIVATIVES)
120  {
121  for (int c = 0; c < VDIM; c++)
122  {
123  der(q,c,0,e) = D[c+VDIM*0];
124  der(q,c,1,e) = D[c+VDIM*1];
125  }
126  }
127  if (VDIM == 2 && (eval_flags & DETERMINANTS))
128  {
129  // The check (VDIM == 2) should eliminate this block when VDIM is
130  // known at compile time and (VDIM != 2).
131  det(q,e) = kernels::Det<2>(D);
132  }
133  }
134  }
135  });
136 }
137 
138 template<const int T_VDIM, const int T_ND, const int T_NQ>
139 void QuadratureInterpolator::Eval3D(
140  const int NE,
141  const int vdim,
142  const DofToQuad &maps,
143  const Vector &e_vec,
144  Vector &q_val,
145  Vector &q_der,
146  Vector &q_det,
147  const int eval_flags)
148 {
149  const int nd = maps.ndof;
150  const int nq = maps.nqpt;
151  const int ND = T_ND ? T_ND : nd;
152  const int NQ = T_NQ ? T_NQ : nq;
153  const int VDIM = T_VDIM ? T_VDIM : vdim;
154  MFEM_VERIFY(ND <= MAX_ND3D, "");
155  MFEM_VERIFY(NQ <= MAX_NQ3D, "");
156  MFEM_VERIFY(VDIM == 3 || !(eval_flags & DETERMINANTS), "");
157  auto B = Reshape(maps.B.Read(), NQ, ND);
158  auto G = Reshape(maps.G.Read(), NQ, 3, ND);
159  auto E = Reshape(e_vec.Read(), ND, VDIM, NE);
160  auto val = Reshape(q_val.Write(), NQ, VDIM, NE);
161  auto der = Reshape(q_der.Write(), NQ, VDIM, 3, NE);
162  auto det = Reshape(q_det.Write(), NQ, NE);
163  MFEM_FORALL(e, NE,
164  {
165  const int ND = T_ND ? T_ND : nd;
166  const int NQ = T_NQ ? T_NQ : nq;
167  const int VDIM = T_VDIM ? T_VDIM : vdim;
168  constexpr int max_ND = T_ND ? T_ND : MAX_ND3D;
169  constexpr int max_VDIM = T_VDIM ? T_VDIM : MAX_VDIM3D;
170  double s_E[max_VDIM*max_ND];
171  for (int d = 0; d < ND; d++)
172  {
173  for (int c = 0; c < VDIM; c++)
174  {
175  s_E[c+d*VDIM] = E(d,c,e);
176  }
177  }
178  for (int q = 0; q < NQ; ++q)
179  {
180  if (eval_flags & VALUES)
181  {
182  double ed[max_VDIM];
183  for (int c = 0; c < VDIM; c++) { ed[c] = 0.0; }
184  for (int d = 0; d < ND; ++d)
185  {
186  const double b = B(q,d);
187  for (int c = 0; c < VDIM; c++) { ed[c] += b*s_E[c+d*VDIM]; }
188  }
189  for (int c = 0; c < VDIM; c++) { val(q,c,e) = ed[c]; }
190  }
191  if ((eval_flags & DERIVATIVES) || (eval_flags & DETERMINANTS))
192  {
193  // use MAX_VDIM3D to avoid "subscript out of range" warnings
194  double D[MAX_VDIM3D*3];
195  for (int i = 0; i < 3*VDIM; i++) { D[i] = 0.0; }
196  for (int d = 0; d < ND; ++d)
197  {
198  const double wx = G(q,0,d);
199  const double wy = G(q,1,d);
200  const double wz = G(q,2,d);
201  for (int c = 0; c < VDIM; c++)
202  {
203  double s_e = s_E[c+d*VDIM];
204  D[c+VDIM*0] += s_e * wx;
205  D[c+VDIM*1] += s_e * wy;
206  D[c+VDIM*2] += s_e * wz;
207  }
208  }
209  if (eval_flags & DERIVATIVES)
210  {
211  for (int c = 0; c < VDIM; c++)
212  {
213  der(q,c,0,e) = D[c+VDIM*0];
214  der(q,c,1,e) = D[c+VDIM*1];
215  der(q,c,2,e) = D[c+VDIM*2];
216  }
217  }
218  if (VDIM == 3 && (eval_flags & DETERMINANTS))
219  {
220  // The check (VDIM == 3) should eliminate this block when VDIM is
221  // known at compile time and (VDIM != 3).
222  det(q,e) = kernels::Det<3>(D);
223  }
224  }
225  }
226  });
227 }
228 
229 void QuadratureInterpolator::Mult(
230  const Vector &e_vec, unsigned eval_flags,
231  Vector &q_val, Vector &q_der, Vector &q_det) const
232 {
233  if (q_layout == QVectorLayout::byVDIM)
234  {
235  if (eval_flags & VALUES) { Values(e_vec, q_val); }
236  if (eval_flags & DERIVATIVES) { Derivatives(e_vec, q_der); }
237  if (eval_flags & DETERMINANTS)
238  {
239  MFEM_ABORT("evaluation of determinants with 'byVDIM' output layout"
240  " is not implemented yet!");
241  }
242  return;
243  }
244 
245  // q_layout == QVectorLayout::byNODES
246  const int ne = fespace->GetNE();
247  if (ne == 0) { return; }
248  const int vdim = fespace->GetVDim();
249  const int dim = fespace->GetMesh()->Dimension();
250  const FiniteElement *fe = fespace->GetFE(0);
251  const IntegrationRule *ir =
252  IntRule ? IntRule : &qspace->GetElementIntRule(0);
253  const DofToQuad &maps = fe->GetDofToQuad(*ir, DofToQuad::FULL);
254  const int nd = maps.ndof;
255  const int nq = maps.nqpt;
256  void (*eval_func)(
257  const int NE,
258  const int vdim,
259  const DofToQuad &maps,
260  const Vector &e_vec,
261  Vector &q_val,
262  Vector &q_der,
263  Vector &q_det,
264  const int eval_flags) = NULL;
265  if (vdim == 1)
266  {
267  if (dim == 2)
268  {
269  switch (100*nd + nq)
270  {
271  // Q0
272  case 101: eval_func = &Eval2D<1,1,1>; break;
273  case 104: eval_func = &Eval2D<1,1,4>; break;
274  // Q1
275  case 404: eval_func = &Eval2D<1,4,4>; break;
276  case 409: eval_func = &Eval2D<1,4,9>; break;
277  // Q2
278  case 909: eval_func = &Eval2D<1,9,9>; break;
279  case 916: eval_func = &Eval2D<1,9,16>; break;
280  // Q3
281  case 1616: eval_func = &Eval2D<1,16,16>; break;
282  case 1625: eval_func = &Eval2D<1,16,25>; break;
283  case 1636: eval_func = &Eval2D<1,16,36>; break;
284  // Q4
285  case 2525: eval_func = &Eval2D<1,25,25>; break;
286  case 2536: eval_func = &Eval2D<1,25,36>; break;
287  case 2549: eval_func = &Eval2D<1,25,49>; break;
288  case 2564: eval_func = &Eval2D<1,25,64>; break;
289  }
290  if (nq >= 100 || !eval_func)
291  {
292  eval_func = &Eval2D<1>;
293  }
294  }
295  else if (dim == 3)
296  {
297  switch (1000*nd + nq)
298  {
299  // Q0
300  case 1001: eval_func = &Eval3D<1,1,1>; break;
301  case 1008: eval_func = &Eval3D<1,1,8>; break;
302  // Q1
303  case 8008: eval_func = &Eval3D<1,8,8>; break;
304  case 8027: eval_func = &Eval3D<1,8,27>; break;
305  // Q2
306  case 27027: eval_func = &Eval3D<1,27,27>; break;
307  case 27064: eval_func = &Eval3D<1,27,64>; break;
308  // Q3
309  case 64064: eval_func = &Eval3D<1,64,64>; break;
310  case 64125: eval_func = &Eval3D<1,64,125>; break;
311  case 64216: eval_func = &Eval3D<1,64,216>; break;
312  // Q4
313  case 125125: eval_func = &Eval3D<1,125,125>; break;
314  case 125216: eval_func = &Eval3D<1,125,216>; break;
315  }
316  if (nq >= 1000 || !eval_func)
317  {
318  eval_func = &Eval3D<1>;
319  }
320  }
321  }
322  else if (vdim == dim)
323  {
324  if (dim == 2)
325  {
326  switch (100*nd + nq)
327  {
328  // Q1
329  case 404: eval_func = &Eval2D<2,4,4>; break;
330  case 409: eval_func = &Eval2D<2,4,9>; break;
331  // Q2
332  case 909: eval_func = &Eval2D<2,9,9>; break;
333  case 916: eval_func = &Eval2D<2,9,16>; break;
334  // Q3
335  case 1616: eval_func = &Eval2D<2,16,16>; break;
336  case 1625: eval_func = &Eval2D<2,16,25>; break;
337  case 1636: eval_func = &Eval2D<2,16,36>; break;
338  // Q4
339  case 2525: eval_func = &Eval2D<2,25,25>; break;
340  case 2536: eval_func = &Eval2D<2,25,36>; break;
341  case 2549: eval_func = &Eval2D<2,25,49>; break;
342  case 2564: eval_func = &Eval2D<2,25,64>; break;
343  }
344  if (nq >= 100 || !eval_func)
345  {
346  eval_func = &Eval2D<2>;
347  }
348  }
349  else if (dim == 3)
350  {
351  switch (1000*nd + nq)
352  {
353  // Q1
354  case 8008: eval_func = &Eval3D<3,8,8>; break;
355  case 8027: eval_func = &Eval3D<3,8,27>; break;
356  // Q2
357  case 27027: eval_func = &Eval3D<3,27,27>; break;
358  case 27064: eval_func = &Eval3D<3,27,64>; break;
359  // Q3
360  case 64064: eval_func = &Eval3D<3,64,64>; break;
361  case 64125: eval_func = &Eval3D<3,64,125>; break;
362  case 64216: eval_func = &Eval3D<3,64,216>; break;
363  // Q4
364  case 125125: eval_func = &Eval3D<3,125,125>; break;
365  case 125216: eval_func = &Eval3D<3,125,216>; break;
366  }
367  if (nq >= 1000 || !eval_func)
368  {
369  eval_func = &Eval3D<3>;
370  }
371  }
372  }
373  if (eval_func)
374  {
375  eval_func(ne, vdim, maps, e_vec, q_val, q_der, q_det, eval_flags);
376  }
377  else
378  {
379  MFEM_ABORT("case not supported yet");
380  }
381 }
382 
383 void QuadratureInterpolator::MultTranspose(
384  unsigned eval_flags, const Vector &q_val, const Vector &q_der,
385  Vector &e_vec) const
386 {
387  MFEM_ABORT("this method is not implemented yet");
388 }
389 
390 
391 template<int T_VDIM = 0, int T_D1D = 0, int T_Q1D = 0, int T_NBZ = 0>
392 static void D2QValues2D(const int NE,
393  const Array<double> &b_,
394  const Vector &x_,
395  Vector &y_,
396  const int vdim = 1,
397  const int d1d = 0,
398  const int q1d = 0)
399 {
400  const int D1D = T_D1D ? T_D1D : d1d;
401  const int Q1D = T_Q1D ? T_Q1D : q1d;
402  constexpr int NBZ = T_NBZ ? T_NBZ : 1;
403  const int VDIM = T_VDIM ? T_VDIM : vdim;
404 
405  auto b = Reshape(b_.Read(), Q1D, D1D);
406  auto x = Reshape(x_.Read(), D1D, D1D, VDIM, NE);
407  auto y = Reshape(y_.Write(), VDIM, Q1D, Q1D, NE);
408 
409  MFEM_FORALL_2D(e, NE, Q1D, Q1D, NBZ,
410  {
411  const int D1D = T_D1D ? T_D1D : d1d;
412  const int Q1D = T_Q1D ? T_Q1D : q1d;
413  const int VDIM = T_VDIM ? T_VDIM : vdim;
414  constexpr int MQ1 = T_Q1D ? T_Q1D : MAX_Q1D;
415  constexpr int MD1 = T_D1D ? T_D1D : MAX_D1D;
416  constexpr int NBZ = T_NBZ ? T_NBZ : 1;
417  const int zid = MFEM_THREAD_ID(z);
418  MFEM_SHARED double B[MQ1][MD1];
419 
420  MFEM_SHARED double DDz[NBZ][MD1*MD1];
421  double (*DD)[MD1] = (double (*)[MD1])(DDz + zid);
422 
423  MFEM_SHARED double DQz[NBZ][MD1*MQ1];
424  double (*DQ)[MQ1] = (double (*)[MQ1])(DQz + zid);
425 
426  if (zid == 0)
427  {
428  MFEM_FOREACH_THREAD(d,y,D1D)
429  {
430  MFEM_FOREACH_THREAD(q,x,Q1D)
431  {
432  B[q][d] = b(q,d);
433  }
434  }
435  }
436  MFEM_SYNC_THREAD;
437 
438  for (int c = 0; c < VDIM; c++)
439  {
440  MFEM_FOREACH_THREAD(dy,y,D1D)
441  {
442  MFEM_FOREACH_THREAD(dx,x,D1D)
443  {
444  DD[dy][dx] = x(dx,dy,c,e);
445  }
446  }
447  MFEM_SYNC_THREAD;
448  MFEM_FOREACH_THREAD(dy,y,D1D)
449  {
450  MFEM_FOREACH_THREAD(qx,x,Q1D)
451  {
452  double dq = 0.0;
453  for (int dx = 0; dx < D1D; ++dx)
454  {
455  dq += B[qx][dx] * DD[dy][dx];
456  }
457  DQ[dy][qx] = dq;
458  }
459  }
460  MFEM_SYNC_THREAD;
461  MFEM_FOREACH_THREAD(qy,y,Q1D)
462  {
463  MFEM_FOREACH_THREAD(qx,x,Q1D)
464  {
465  double qq = 0.0;
466  for (int dy = 0; dy < D1D; ++dy)
467  {
468  qq += DQ[dy][qx] * B[qy][dy];
469  }
470  y(c,qx,qy,e) = qq;
471  }
472  }
473  MFEM_SYNC_THREAD;
474  }
475  });
476 }
477 
478 template<int T_VDIM = 0, int T_D1D = 0, int T_Q1D = 0,
479  int MAX_D = 0, int MAX_Q = 0>
480 static void D2QValues3D(const int NE,
481  const Array<double> &b_,
482  const Vector &x_,
483  Vector &y_,
484  const int vdim = 1,
485  const int d1d = 0,
486  const int q1d = 0)
487 {
488  const int D1D = T_D1D ? T_D1D : d1d;
489  const int Q1D = T_Q1D ? T_Q1D : q1d;
490  const int VDIM = T_VDIM ? T_VDIM : vdim;
491 
492  auto b = Reshape(b_.Read(), Q1D, D1D);
493  auto x = Reshape(x_.Read(), D1D, D1D, D1D, VDIM, NE);
494  auto y = Reshape(y_.Write(), VDIM, Q1D, Q1D, Q1D, NE);
495 
496  MFEM_FORALL_3D(e, NE, Q1D, Q1D, Q1D,
497  {
498  const int D1D = T_D1D ? T_D1D : d1d;
499  const int Q1D = T_Q1D ? T_Q1D : q1d;
500  const int VDIM = T_VDIM ? T_VDIM : vdim;
501  constexpr int MQ1 = T_Q1D ? T_Q1D : MAX_Q;
502  constexpr int MD1 = T_D1D ? T_D1D : MAX_D;
503  constexpr int MDQ = (MQ1 > MD1) ? MQ1 : MD1;
504  const int tidz = MFEM_THREAD_ID(z);
505  MFEM_SHARED double B[MQ1][MD1];
506  MFEM_SHARED double sm0[MDQ*MDQ*MDQ];
507  MFEM_SHARED double sm1[MDQ*MDQ*MDQ];
508  double (*X)[MD1][MD1] = (double (*)[MD1][MD1]) sm0;
509  double (*DDQ)[MD1][MQ1] = (double (*)[MD1][MQ1]) sm1;
510  double (*DQQ)[MQ1][MQ1] = (double (*)[MQ1][MQ1]) sm0;
511 
512  if (tidz == 0)
513  {
514  MFEM_FOREACH_THREAD(d,y,D1D)
515  {
516  MFEM_FOREACH_THREAD(q,x,Q1D)
517  {
518  B[q][d] = b(q,d);
519  }
520  }
521  }
522  MFEM_SYNC_THREAD;
523 
524  for (int c = 0; c < VDIM; c++)
525  {
526  MFEM_FOREACH_THREAD(dz,z,D1D)
527  {
528  MFEM_FOREACH_THREAD(dy,y,D1D)
529  {
530  MFEM_FOREACH_THREAD(dx,x,D1D)
531  {
532  X[dz][dy][dx] = x(dx,dy,dz,c,e);
533  }
534  }
535  }
536  MFEM_SYNC_THREAD;
537  MFEM_FOREACH_THREAD(dz,z,D1D)
538  {
539  MFEM_FOREACH_THREAD(dy,y,D1D)
540  {
541  MFEM_FOREACH_THREAD(qx,x,Q1D)
542  {
543  double u = 0.0;
544  for (int dx = 0; dx < D1D; ++dx)
545  {
546  u += B[qx][dx] * X[dz][dy][dx];
547  }
548  DDQ[dz][dy][qx] = u;
549  }
550  }
551  }
552  MFEM_SYNC_THREAD;
553  MFEM_FOREACH_THREAD(dz,z,D1D)
554  {
555  MFEM_FOREACH_THREAD(qy,y,Q1D)
556  {
557  MFEM_FOREACH_THREAD(qx,x,Q1D)
558  {
559  double u = 0.0;
560  for (int dy = 0; dy < D1D; ++dy)
561  {
562  u += DDQ[dz][dy][qx] * B[qy][dy];
563  }
564  DQQ[dz][qy][qx] = u;
565  }
566  }
567  }
568  MFEM_SYNC_THREAD;
569  MFEM_FOREACH_THREAD(qz,z,Q1D)
570  {
571  MFEM_FOREACH_THREAD(qy,y,Q1D)
572  {
573  MFEM_FOREACH_THREAD(qx,x,Q1D)
574  {
575  double u = 0.0;
576  for (int dz = 0; dz < D1D; ++dz)
577  {
578  u += DQQ[dz][qy][qx] * B[qz][dz];
579  }
580  y(c,qx,qy,qz,e) = u;
581  }
582  }
583  }
584  MFEM_SYNC_THREAD;
585  }
586  });
587 }
588 
589 static void D2QValues(const FiniteElementSpace &fes,
590  const DofToQuad *maps,
591  const Vector &e_vec,
592  Vector &q_val)
593 {
594  const int dim = fes.GetMesh()->Dimension();
595  const int vdim = fes.GetVDim();
596  const int NE = fes.GetNE();
597  const int D1D = maps->ndof;
598  const int Q1D = maps->nqpt;
599  const int id = (vdim<<8) | (D1D<<4) | Q1D;
600 
601  if (dim == 2)
602  {
603  switch (id)
604  {
605  case 0x124: return D2QValues2D<1,2,4,8>(NE, maps->B, e_vec, q_val);
606  case 0x136: return D2QValues2D<1,3,6,4>(NE, maps->B, e_vec, q_val);
607  case 0x148: return D2QValues2D<1,4,8,2>(NE, maps->B, e_vec, q_val);
608  case 0x224: return D2QValues2D<2,2,4,8>(NE, maps->B, e_vec, q_val);
609  case 0x236: return D2QValues2D<2,3,6,4>(NE, maps->B, e_vec, q_val);
610  case 0x248: return D2QValues2D<2,4,8,2>(NE, maps->B, e_vec, q_val);
611  default:
612  {
613  MFEM_VERIFY(D1D <= MAX_D1D, "Orders higher than " << MAX_D1D-1
614  << " are not supported!");
615  MFEM_VERIFY(Q1D <= MAX_Q1D, "Quadrature rules with more than "
616  << MAX_Q1D << " 1D points are not supported!");
617  D2QValues2D(NE, maps->B, e_vec, q_val, vdim, D1D, Q1D);
618  return;
619  }
620  }
621  }
622  if (dim == 3)
623  {
624  switch (id)
625  {
626  case 0x124: return D2QValues3D<1,2,4>(NE, maps->B, e_vec, q_val);
627  case 0x136: return D2QValues3D<1,3,6>(NE, maps->B, e_vec, q_val);
628  case 0x148: return D2QValues3D<1,4,8>(NE, maps->B, e_vec, q_val);
629  case 0x324: return D2QValues3D<3,2,4>(NE, maps->B, e_vec, q_val);
630  case 0x336: return D2QValues3D<3,3,6>(NE, maps->B, e_vec, q_val);
631  case 0x348: return D2QValues3D<3,4,8>(NE, maps->B, e_vec, q_val);
632  default:
633  {
634  constexpr int MD = 8;
635  constexpr int MQ = 8;
636  MFEM_VERIFY(D1D <= MD, "Orders higher than " << MD-1
637  << " are not supported!");
638  MFEM_VERIFY(Q1D <= MQ, "Quadrature rules with more than " << MQ
639  << " 1D points are not supported!");
640  D2QValues3D<0,0,0,MD,MQ>(NE, maps->B, e_vec, q_val, vdim, D1D, Q1D);
641  return;
642  }
643  }
644  }
645  mfem::out << "Unknown kernel 0x" << std::hex << id << std::endl;
646  MFEM_ABORT("Unknown kernel");
647 }
648 
649 void QuadratureInterpolator::Values(const Vector &e_vec, Vector &q_val) const
650 {
651  if (q_layout == QVectorLayout::byNODES)
652  {
653  Vector empty;
654  Mult(e_vec, VALUES, q_val, empty, empty);
655  return;
656  }
657 
658  // q_layout == QVectorLayout::byVDIM
659  if (fespace->GetNE() == 0) { return; }
660  const IntegrationRule &ir = *IntRule;
661  const DofToQuad::Mode mode = DofToQuad::TENSOR;
662  const DofToQuad &d2q = fespace->GetFE(0)->GetDofToQuad(ir, mode);
663  D2QValues(*fespace, &d2q, e_vec, q_val);
664 }
665 
666 template<int T_VDIM = 0, int T_D1D = 0, int T_Q1D = 0, int T_NBZ = 0>
667 static void D2QGrad2D(const int NE,
668  const double *b_,
669  const double *g_,
670  const double *x_,
671  double *y_,
672  const int vdim = 1,
673  const int d1d = 0,
674  const int q1d = 0)
675 {
676  const int D1D = T_D1D ? T_D1D : d1d;
677  const int Q1D = T_Q1D ? T_Q1D : q1d;
678  constexpr int NBZ = T_NBZ ? T_NBZ : 1;
679  const int VDIM = T_VDIM ? T_VDIM : vdim;
680 
681  auto b = Reshape(b_, Q1D, D1D);
682  auto g = Reshape(g_, Q1D, D1D);
683  auto x = Reshape(x_, D1D, D1D, VDIM, NE);
684  auto y = Reshape(y_, VDIM, 2, Q1D, Q1D, NE);
685 
686  MFEM_FORALL_2D(e, NE, Q1D, Q1D, NBZ,
687  {
688  const int D1D = T_D1D ? T_D1D : d1d;
689  const int Q1D = T_Q1D ? T_Q1D : q1d;
690  const int VDIM = T_VDIM ? T_VDIM : vdim;
691  constexpr int MQ1 = T_Q1D ? T_Q1D : MAX_Q1D;
692  constexpr int MD1 = T_D1D ? T_D1D : MAX_D1D;
693  constexpr int NBZ = T_NBZ ? T_NBZ : 1;
694  const int tidz = MFEM_THREAD_ID(z);
695  MFEM_SHARED double B[MQ1][MD1];
696  MFEM_SHARED double G[MQ1][MD1];
697 
698  MFEM_SHARED double Xz[NBZ][MD1][MD1];
699  double (*X)[MD1] = (double (*)[MD1])(Xz + tidz);
700 
701  MFEM_SHARED double GD[2][NBZ][MD1][MQ1];
702  double (*DQ0)[MQ1] = (double (*)[MQ1])(GD[0] + tidz);
703  double (*DQ1)[MQ1] = (double (*)[MQ1])(GD[1] + tidz);
704 
705  if (tidz == 0)
706  {
707  MFEM_FOREACH_THREAD(d,y,D1D)
708  {
709  MFEM_FOREACH_THREAD(q,x,Q1D)
710  {
711  B[q][d] = b(q,d);
712  G[q][d] = g(q,d);
713  }
714  }
715  }
716  MFEM_SYNC_THREAD;
717 
718  for (int c = 0; c < VDIM; ++c)
719  {
720  MFEM_FOREACH_THREAD(dx,x,D1D)
721  {
722  MFEM_FOREACH_THREAD(dy,y,D1D)
723  {
724  X[dx][dy] = x(dx,dy,c,e);
725  }
726  }
727  MFEM_SYNC_THREAD;
728  MFEM_FOREACH_THREAD(dy,y,D1D)
729  {
730  MFEM_FOREACH_THREAD(qx,x,Q1D)
731  {
732  double u = 0.0;
733  double v = 0.0;
734  for (int dx = 0; dx < D1D; ++dx)
735  {
736  const double input = X[dx][dy];
737  u += B[qx][dx] * input;
738  v += G[qx][dx] * input;
739  }
740  DQ0[dy][qx] = u;
741  DQ1[dy][qx] = v;
742  }
743  }
744  MFEM_SYNC_THREAD;
745  MFEM_FOREACH_THREAD(qy,y,Q1D)
746  {
747  MFEM_FOREACH_THREAD(qx,x,Q1D)
748  {
749  double u = 0.0;
750  double v = 0.0;
751  for (int dy = 0; dy < D1D; ++dy)
752  {
753  u += DQ1[dy][qx] * B[qy][dy];
754  v += DQ0[dy][qx] * G[qy][dy];
755  }
756  y(c,0,qx,qy,e) = u;
757  y(c,1,qx,qy,e) = v;
758  }
759  }
760  MFEM_SYNC_THREAD;
761  }
762  });
763 }
764 
765 template<int T_VDIM = 0, int T_D1D = 0, int T_Q1D = 0,
766  int MAX_D = 0, int MAX_Q = 0>
767 static void D2QGrad3D(const int NE,
768  const double *b_,
769  const double *g_,
770  const double *x_,
771  double *y_,
772  const int vdim = 1,
773  const int d1d = 0,
774  const int q1d = 0)
775 {
776  const int D1D = T_D1D ? T_D1D : d1d;
777  const int Q1D = T_Q1D ? T_Q1D : q1d;
778  const int VDIM = T_VDIM ? T_VDIM : vdim;
779 
780  auto b = Reshape(b_, Q1D, D1D);
781  auto g = Reshape(g_, Q1D, D1D);
782  auto x = Reshape(x_, D1D, D1D, D1D, VDIM, NE);
783  auto y = Reshape(y_, VDIM, 3, Q1D, Q1D, Q1D, NE);
784 
785  MFEM_FORALL_3D(e, NE, Q1D, Q1D, Q1D,
786  {
787  const int D1D = T_D1D ? T_D1D : d1d;
788  const int Q1D = T_Q1D ? T_Q1D : q1d;
789  const int VDIM = T_VDIM ? T_VDIM : vdim;
790  constexpr int MQ1 = T_Q1D ? T_Q1D : MAX_Q;
791  constexpr int MD1 = T_D1D ? T_D1D : MAX_D;
792  const int tidz = MFEM_THREAD_ID(z);
793  MFEM_SHARED double B[MQ1][MD1];
794  MFEM_SHARED double G[MQ1][MD1];
795 
796  MFEM_SHARED double sm0[3][MQ1*MQ1*MQ1];
797  MFEM_SHARED double sm1[3][MQ1*MQ1*MQ1];
798  double (*X)[MD1][MD1] = (double (*)[MD1][MD1]) (sm0+2);
799  double (*DDQ0)[MD1][MQ1] = (double (*)[MD1][MQ1]) (sm0+0);
800  double (*DDQ1)[MD1][MQ1] = (double (*)[MD1][MQ1]) (sm0+1);
801  double (*DQQ0)[MQ1][MQ1] = (double (*)[MQ1][MQ1]) (sm1+0);
802  double (*DQQ1)[MQ1][MQ1] = (double (*)[MQ1][MQ1]) (sm1+1);
803  double (*DQQ2)[MQ1][MQ1] = (double (*)[MQ1][MQ1]) (sm1+2);
804 
805  if (tidz == 0)
806  {
807  MFEM_FOREACH_THREAD(d,y,D1D)
808  {
809  MFEM_FOREACH_THREAD(q,x,Q1D)
810  {
811  B[q][d] = b(q,d);
812  G[q][d] = g(q,d);
813  }
814  }
815  }
816  MFEM_SYNC_THREAD;
817 
818  for (int c = 0; c < VDIM; ++c)
819  {
820  MFEM_FOREACH_THREAD(dx,x,D1D)
821  {
822  MFEM_FOREACH_THREAD(dy,y,D1D)
823  {
824  MFEM_FOREACH_THREAD(dz,z,D1D)
825  {
826  X[dx][dy][dz] = x(dx,dy,dz,c,e);
827  }
828  }
829  }
830  MFEM_SYNC_THREAD;
831 
832  MFEM_FOREACH_THREAD(dz,z,D1D)
833  {
834  MFEM_FOREACH_THREAD(dy,y,D1D)
835  {
836  MFEM_FOREACH_THREAD(qx,x,Q1D)
837  {
838  double u = 0.0;
839  double v = 0.0;
840  for (int dx = 0; dx < D1D; ++dx)
841  {
842  const double coords = X[dx][dy][dz];
843  u += coords * B[qx][dx];
844  v += coords * G[qx][dx];
845  }
846  DDQ0[dz][dy][qx] = u;
847  DDQ1[dz][dy][qx] = v;
848  }
849  }
850  }
851  MFEM_SYNC_THREAD;
852  MFEM_FOREACH_THREAD(dz,z,D1D)
853  {
854  MFEM_FOREACH_THREAD(qy,y,Q1D)
855  {
856  MFEM_FOREACH_THREAD(qx,x,Q1D)
857  {
858  double u = 0.0;
859  double v = 0.0;
860  double w = 0.0;
861  for (int dy = 0; dy < D1D; ++dy)
862  {
863  u += DDQ1[dz][dy][qx] * B[qy][dy];
864  v += DDQ0[dz][dy][qx] * G[qy][dy];
865  w += DDQ0[dz][dy][qx] * B[qy][dy];
866  }
867  DQQ0[dz][qy][qx] = u;
868  DQQ1[dz][qy][qx] = v;
869  DQQ2[dz][qy][qx] = w;
870  }
871  }
872  }
873  MFEM_SYNC_THREAD;
874  MFEM_FOREACH_THREAD(qz,z,Q1D)
875  {
876  MFEM_FOREACH_THREAD(qy,y,Q1D)
877  {
878  MFEM_FOREACH_THREAD(qx,x,Q1D)
879  {
880  double u = 0.0;
881  double v = 0.0;
882  double w = 0.0;
883  for (int dz = 0; dz < D1D; ++dz)
884  {
885  u += DQQ0[dz][qy][qx] * B[qz][dz];
886  v += DQQ1[dz][qy][qx] * B[qz][dz];
887  w += DQQ2[dz][qy][qx] * G[qz][dz];
888  }
889  y(c,0,qx,qy,qz,e) = u;
890  y(c,1,qx,qy,qz,e) = v;
891  y(c,2,qx,qy,qz,e) = w;
892  }
893  }
894  }
895  MFEM_SYNC_THREAD;
896  }
897  });
898 }
899 
900 static void D2QGrad(const FiniteElementSpace &fes,
901  const DofToQuad *maps,
902  const Vector &e_vec,
903  Vector &q_der)
904 {
905  const int dim = fes.GetMesh()->Dimension();
906  const int vdim = fes.GetVDim();
907  const int NE = fes.GetNE();
908  const int D1D = maps->ndof;
909  const int Q1D = maps->nqpt;
910  const int id = (vdim<<8) | (D1D<<4) | Q1D;
911  const double *B = maps->B.Read();
912  const double *G = maps->G.Read();
913  const double *X = e_vec.Read();
914  double *Y = q_der.Write();
915  if (dim == 2)
916  {
917  switch (id)
918  {
919  case 0x134: return D2QGrad2D<1,3,4,8>(NE, B, G, X, Y);
920  case 0x146: return D2QGrad2D<1,4,6,4>(NE, B, G, X, Y);
921  case 0x158: return D2QGrad2D<1,5,8,2>(NE, B, G, X, Y);
922  case 0x234: return D2QGrad2D<2,3,4,8>(NE, B, G, X, Y);
923  case 0x246: return D2QGrad2D<2,4,6,4>(NE, B, G, X, Y);
924  case 0x258: return D2QGrad2D<2,5,8,2>(NE, B, G, X, Y);
925  default:
926  {
927  MFEM_VERIFY(D1D <= MAX_D1D, "Orders higher than " << MAX_D1D-1
928  << " are not supported!");
929  MFEM_VERIFY(Q1D <= MAX_Q1D, "Quadrature rules with more than "
930  << MAX_Q1D << " 1D points are not supported!");
931  D2QGrad2D(NE, B, G, X, Y, vdim, D1D, Q1D);
932  return;
933  }
934  }
935  }
936  if (dim == 3)
937  {
938  switch (id)
939  {
940  case 0x134: return D2QGrad3D<1,3,4>(NE, B, G, X, Y);
941  case 0x146: return D2QGrad3D<1,4,6>(NE, B, G, X, Y);
942  case 0x158: return D2QGrad3D<1,5,8>(NE, B, G, X, Y);
943  case 0x334: return D2QGrad3D<3,3,4>(NE, B, G, X, Y);
944  case 0x346: return D2QGrad3D<3,4,6>(NE, B, G, X, Y);
945  case 0x358: return D2QGrad3D<3,5,8>(NE, B, G, X, Y);
946  default:
947  {
948  constexpr int MD = 8;
949  constexpr int MQ = 8;
950  MFEM_VERIFY(D1D <= MD, "Orders higher than " << MD-1
951  << " are not supported!");
952  MFEM_VERIFY(Q1D <= MQ, "Quadrature rules with more than " << MQ
953  << " 1D points are not supported!");
954  D2QGrad3D<0,0,0,MD,MQ>(NE, B, G, X, Y, vdim, D1D, Q1D);
955  return;
956  }
957  }
958  }
959  mfem::out << "Unknown kernel 0x" << std::hex << id << std::endl;
960  MFEM_ABORT("Unknown kernel");
961 }
962 
963 template<int T_VDIM = 0, int T_D1D = 0, int T_Q1D = 0, int T_NBZ = 0>
964 static void D2QPhysGrad2D(const int NE,
965  const double *b_,
966  const double *g_,
967  const double *j_,
968  const double *x_,
969  double *y_,
970  const int vdim = 1,
971  const int d1d = 0,
972  const int q1d = 0)
973 {
974  const int D1D = T_D1D ? T_D1D : d1d;
975  const int Q1D = T_Q1D ? T_Q1D : q1d;
976  constexpr int NBZ = T_NBZ ? T_NBZ : 1;
977  const int VDIM = T_VDIM ? T_VDIM : vdim;
978 
979  auto b = Reshape(b_, Q1D, D1D);
980  auto g = Reshape(g_, Q1D, D1D);
981  auto j = Reshape(j_, Q1D, Q1D, 2, 2, NE);
982  auto x = Reshape(x_, D1D, D1D, VDIM, NE);
983  auto y = Reshape(y_, VDIM, 2, Q1D, Q1D, NE);
984 
985  MFEM_FORALL_2D(e, NE, Q1D, Q1D, NBZ,
986  {
987  const int D1D = T_D1D ? T_D1D : d1d;
988  const int Q1D = T_Q1D ? T_Q1D : q1d;
989  const int VDIM = T_VDIM ? T_VDIM : vdim;
990  constexpr int NBZ = T_NBZ ? T_NBZ : 1;
991  constexpr int MQ1 = T_Q1D ? T_Q1D : MAX_Q1D;
992  constexpr int MD1 = T_D1D ? T_D1D : MAX_D1D;
993  const int tidz = MFEM_THREAD_ID(z);
994  MFEM_SHARED double B[MQ1][MD1];
995  MFEM_SHARED double G[MQ1][MD1];
996 
997  MFEM_SHARED double Xz[NBZ][MD1][MD1];
998  double (*X)[MD1] = (double (*)[MD1])(Xz + tidz);
999 
1000  MFEM_SHARED double GD[2][NBZ][MD1][MQ1];
1001  double (*DQ0)[MQ1] = (double (*)[MQ1])(GD[0] + tidz);
1002  double (*DQ1)[MQ1] = (double (*)[MQ1])(GD[1] + tidz);
1003 
1004  if (tidz == 0)
1005  {
1006  MFEM_FOREACH_THREAD(d,y,D1D)
1007  {
1008  MFEM_FOREACH_THREAD(q,x,Q1D)
1009  {
1010  B[q][d] = b(q,d);
1011  G[q][d] = g(q,d);
1012  }
1013  }
1014  }
1015  MFEM_SYNC_THREAD;
1016 
1017  for (int c = 0; c < VDIM; ++c)
1018  {
1019  MFEM_FOREACH_THREAD(dx,x,D1D)
1020  {
1021  MFEM_FOREACH_THREAD(dy,y,D1D)
1022  {
1023  X[dx][dy] = x(dx,dy,c,e);
1024  }
1025  }
1026  MFEM_SYNC_THREAD;
1027  MFEM_FOREACH_THREAD(dy,y,D1D)
1028  {
1029  MFEM_FOREACH_THREAD(qx,x,Q1D)
1030  {
1031  double u = 0.0;
1032  double v = 0.0;
1033  for (int dx = 0; dx < D1D; ++dx)
1034  {
1035  const double input = X[dx][dy];
1036  u += B[qx][dx] * input;
1037  v += G[qx][dx] * input;
1038  }
1039  DQ0[dy][qx] = u;
1040  DQ1[dy][qx] = v;
1041  }
1042  }
1043  MFEM_SYNC_THREAD;
1044 
1045  MFEM_FOREACH_THREAD(qy,y,Q1D)
1046  {
1047  MFEM_FOREACH_THREAD(qx,x,Q1D)
1048  {
1049  double u = 0.0;
1050  double v = 0.0;
1051  for (int dy = 0; dy < D1D; ++dy)
1052  {
1053  u += DQ1[dy][qx] * B[qy][dy];
1054  v += DQ0[dy][qx] * G[qy][dy];
1055  }
1056  double Jloc[4], Jinv[4];
1057  Jloc[0] = j(qx,qy,0,0,e);
1058  Jloc[1] = j(qx,qy,1,0,e);
1059  Jloc[2] = j(qx,qy,0,1,e);
1060  Jloc[3] = j(qx,qy,1,1,e);
1061  kernels::CalcInverse<2>(Jloc, Jinv);
1062  y(c,0,qx,qy,e) = Jinv[0]*u + Jinv[1]*v;
1063  y(c,1,qx,qy,e) = Jinv[2]*u + Jinv[3]*v;
1064  }
1065  }
1066  MFEM_SYNC_THREAD;
1067  }
1068  });
1069 }
1070 
1071 template<int T_VDIM = 0, int T_D1D = 0, int T_Q1D = 0,
1072  int MAX_D = 0, int MAX_Q = 0>
1073 static void D2QPhysGrad3D(const int NE,
1074  const double *b_,
1075  const double *g_,
1076  const double *j_,
1077  const double *x_,
1078  double *y_,
1079  const int vdim = 1,
1080  const int d1d = 0,
1081  const int q1d = 0)
1082 {
1083  const int D1D = T_D1D ? T_D1D : d1d;
1084  const int Q1D = T_Q1D ? T_Q1D : q1d;
1085  const int VDIM = T_VDIM ? T_VDIM : vdim;
1086 
1087  auto b = Reshape(b_, Q1D, D1D);
1088  auto g = Reshape(g_, Q1D, D1D);
1089  auto j = Reshape(j_, Q1D, Q1D, Q1D, 3, 3, NE);
1090  auto x = Reshape(x_, D1D, D1D, D1D, VDIM, NE);
1091  auto y = Reshape(y_, VDIM, 3, Q1D, Q1D, Q1D, NE);
1092 
1093  MFEM_FORALL_3D(e, NE, Q1D, Q1D, Q1D,
1094  {
1095  const int D1D = T_D1D ? T_D1D : d1d;
1096  const int Q1D = T_Q1D ? T_Q1D : q1d;
1097  const int VDIM = T_VDIM ? T_VDIM : vdim;
1098  constexpr int MQ1 = T_Q1D ? T_Q1D : MAX_Q;
1099  constexpr int MD1 = T_D1D ? T_D1D : MAX_D;
1100  const int tidz = MFEM_THREAD_ID(z);
1101  MFEM_SHARED double B[MQ1][MD1];
1102  MFEM_SHARED double G[MQ1][MD1];
1103 
1104  MFEM_SHARED double sm0[3][MQ1*MQ1*MQ1];
1105  MFEM_SHARED double sm1[3][MQ1*MQ1*MQ1];
1106  double (*X)[MD1][MD1] = (double (*)[MD1][MD1]) (sm0+2);
1107  double (*DDQ0)[MD1][MQ1] = (double (*)[MD1][MQ1]) (sm0+0);
1108  double (*DDQ1)[MD1][MQ1] = (double (*)[MD1][MQ1]) (sm0+1);
1109  double (*DQQ0)[MQ1][MQ1] = (double (*)[MQ1][MQ1]) (sm1+0);
1110  double (*DQQ1)[MQ1][MQ1] = (double (*)[MQ1][MQ1]) (sm1+1);
1111  double (*DQQ2)[MQ1][MQ1] = (double (*)[MQ1][MQ1]) (sm1+2);
1112 
1113  if (tidz == 0)
1114  {
1115  MFEM_FOREACH_THREAD(d,y,D1D)
1116  {
1117  MFEM_FOREACH_THREAD(q,x,Q1D)
1118  {
1119  B[q][d] = b(q,d);
1120  G[q][d] = g(q,d);
1121  }
1122  }
1123  }
1124  MFEM_SYNC_THREAD;
1125 
1126  for (int c = 0; c < VDIM; ++c)
1127  {
1128  MFEM_FOREACH_THREAD(dx,x,D1D)
1129  {
1130  MFEM_FOREACH_THREAD(dy,y,D1D)
1131  {
1132  MFEM_FOREACH_THREAD(dz,z,D1D)
1133  {
1134 
1135  X[dx][dy][dz] = x(dx,dy,dz,c,e);
1136  }
1137  }
1138  }
1139  MFEM_SYNC_THREAD;
1140 
1141  MFEM_FOREACH_THREAD(dz,z,D1D)
1142  {
1143  MFEM_FOREACH_THREAD(dy,y,D1D)
1144  {
1145  MFEM_FOREACH_THREAD(qx,x,Q1D)
1146  {
1147  double u = 0.0;
1148  double v = 0.0;
1149  for (int dx = 0; dx < D1D; ++dx)
1150  {
1151  const double coords = X[dx][dy][dz];
1152  u += coords * B[qx][dx];
1153  v += coords * G[qx][dx];
1154  }
1155  DDQ0[dz][dy][qx] = u;
1156  DDQ1[dz][dy][qx] = v;
1157  }
1158  }
1159  }
1160  MFEM_SYNC_THREAD;
1161  MFEM_FOREACH_THREAD(dz,z,D1D)
1162  {
1163  MFEM_FOREACH_THREAD(qy,y,Q1D)
1164  {
1165  MFEM_FOREACH_THREAD(qx,x,Q1D)
1166  {
1167  double u = 0.0;
1168  double v = 0.0;
1169  double w = 0.0;
1170  for (int dy = 0; dy < D1D; ++dy)
1171  {
1172  u += DDQ1[dz][dy][qx] * B[qy][dy];
1173  v += DDQ0[dz][dy][qx] * G[qy][dy];
1174  w += DDQ0[dz][dy][qx] * B[qy][dy];
1175  }
1176  DQQ0[dz][qy][qx] = u;
1177  DQQ1[dz][qy][qx] = v;
1178  DQQ2[dz][qy][qx] = w;
1179  }
1180  }
1181  }
1182  MFEM_SYNC_THREAD;
1183  MFEM_FOREACH_THREAD(qz,z,Q1D)
1184  {
1185  MFEM_FOREACH_THREAD(qy,y,Q1D)
1186  {
1187  MFEM_FOREACH_THREAD(qx,x,Q1D)
1188  {
1189  double u = 0.0;
1190  double v = 0.0;
1191  double w = 0.0;
1192  for (int dz = 0; dz < D1D; ++dz)
1193  {
1194  u += DQQ0[dz][qy][qx] * B[qz][dz];
1195  v += DQQ1[dz][qy][qx] * B[qz][dz];
1196  w += DQQ2[dz][qy][qx] * G[qz][dz];
1197  }
1198  double Jloc[9], Jinv[9];
1199  for (int col = 0; col < 3; col++)
1200  {
1201  for (int row = 0; row < 3; row++)
1202  {
1203  Jloc[row+3*col] = j(qx,qy,qz,row,col,e);
1204  }
1205  }
1206  kernels::CalcInverse<3>(Jloc, Jinv);
1207  y(c,0,qx,qy,qz,e) = Jinv[0]*u + Jinv[1]*v + Jinv[2]*w;
1208  y(c,1,qx,qy,qz,e) = Jinv[3]*u + Jinv[4]*v + Jinv[5]*w;
1209  y(c,2,qx,qy,qz,e) = Jinv[6]*u + Jinv[7]*v + Jinv[8]*w;
1210  }
1211  }
1212  }
1213  MFEM_SYNC_THREAD;
1214  }
1215  });
1216 }
1217 
1218 
1219 static void D2QPhysGrad(const FiniteElementSpace &fes,
1220  const GeometricFactors *geom,
1221  const DofToQuad *maps,
1222  const Vector &e_vec,
1223  Vector &q_der)
1224 {
1225  const int dim = fes.GetMesh()->Dimension();
1226  const int vdim = fes.GetVDim();
1227  const int NE = fes.GetNE();
1228  const int D1D = maps->ndof;
1229  const int Q1D = maps->nqpt;
1230  const int id = (vdim<<8) | (D1D<<4) | Q1D;
1231  const double *B = maps->B.Read();
1232  const double *G = maps->G.Read();
1233  const double *J = geom->J.Read();
1234  const double *X = e_vec.Read();
1235  double *Y = q_der.Write();
1236  if (dim == 2)
1237  {
1238  switch (id)
1239  {
1240  case 0x134: return D2QPhysGrad2D<1,3,4,8>(NE, B, G, J, X, Y);
1241  case 0x146: return D2QPhysGrad2D<1,4,6,4>(NE, B, G, J, X, Y);
1242  case 0x158: return D2QPhysGrad2D<1,5,8,2>(NE, B, G, J, X, Y);
1243  case 0x234: return D2QPhysGrad2D<2,3,4,8>(NE, B, G, J, X, Y);
1244  case 0x246: return D2QPhysGrad2D<2,4,6,4>(NE, B, G, J, X, Y);
1245  case 0x258: return D2QPhysGrad2D<2,5,8,2>(NE, B, G, J, X, Y);
1246  default:
1247  {
1248  MFEM_VERIFY(D1D <= MAX_D1D, "Orders higher than " << MAX_D1D-1
1249  << " are not supported!");
1250  MFEM_VERIFY(Q1D <= MAX_Q1D, "Quadrature rules with more than "
1251  << MAX_Q1D << " 1D points are not supported!");
1252  D2QPhysGrad2D(NE, B, G, J, X, Y, vdim, D1D, Q1D);
1253  return;
1254  }
1255  }
1256  }
1257  if (dim == 3)
1258  {
1259  switch (id)
1260  {
1261  case 0x134: return D2QPhysGrad3D<1,3,4>(NE, B, G, J, X, Y);
1262  case 0x146: return D2QPhysGrad3D<1,4,6>(NE, B, G, J, X, Y);
1263  case 0x158: return D2QPhysGrad3D<1,5,8>(NE, B, G, J, X, Y);
1264  case 0x334: return D2QPhysGrad3D<3,3,4>(NE, B, G, J, X, Y);
1265  case 0x346: return D2QPhysGrad3D<3,4,6>(NE, B, G, J, X, Y);
1266  case 0x358: return D2QPhysGrad3D<3,5,8>(NE, B, G, J, X, Y);
1267  default:
1268  {
1269  constexpr int MD = 8;
1270  constexpr int MQ = 8;
1271  MFEM_VERIFY(D1D <= MD, "Orders higher than " << MD-1
1272  << " are not supported!");
1273  MFEM_VERIFY(Q1D <= MQ, "Quadrature rules with more than " << MQ
1274  << " 1D points are not supported!");
1275  D2QPhysGrad3D<0,0,0,MD,MQ>(NE, B, G, J, X, Y, vdim, D1D, Q1D);
1276  return;
1277  }
1278  }
1279  }
1280  mfem::out << "Unknown kernel 0x" << std::hex << id << std::endl;
1281  MFEM_ABORT("Unknown kernel");
1282 }
1283 
1284 void QuadratureInterpolator::Derivatives(const Vector &e_vec,
1285  Vector &q_der) const
1286 {
1287  if (q_layout == QVectorLayout::byNODES)
1288  {
1289  Vector empty;
1290  Mult(e_vec, DERIVATIVES, empty, q_der, empty);
1291  return;
1292  }
1293 
1294  // q_layout == QVectorLayout::byVDIM
1295  if (fespace->GetNE() == 0) { return; }
1296  const IntegrationRule &ir = *IntRule;
1297  const DofToQuad::Mode mode = DofToQuad::TENSOR;
1298  const DofToQuad &d2q = fespace->GetFE(0)->GetDofToQuad(ir, mode);
1299  D2QGrad(*fespace, &d2q, e_vec, q_der);
1300 }
1301 
1302 void QuadratureInterpolator::PhysDerivatives(const Vector &e_vec,
1303  Vector &q_der) const
1304 {
1305  if (q_layout == QVectorLayout::byNODES)
1306  {
1307  MFEM_ABORT("evaluation of physical derivatives with 'byNODES' output"
1308  " layout is not implemented yet!");
1309  return;
1310  }
1311 
1312  // q_layout == QVectorLayout::byVDIM
1313  Mesh *mesh = fespace->GetMesh();
1314  if (mesh->GetNE() == 0) { return; }
1315  // mesh->DeleteGeometricFactors(); // This should be done outside
1316  const IntegrationRule &ir = *IntRule;
1317  const GeometricFactors *geom =
1318  mesh->GetGeometricFactors(ir, GeometricFactors::JACOBIANS);
1319  const DofToQuad::Mode mode = DofToQuad::TENSOR;
1320  const DofToQuad &d2q = fespace->GetFE(0)->GetDofToQuad(ir, mode);
1321  D2QPhysGrad(*fespace, geom, &d2q, e_vec, q_der);
1322 }
1323 
1324 } // namespace mfem
mfem::FiniteElement
Abstract class for Finite Elements.
Definition: fe.hpp:232
mfem::IntegrationRule
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:90
mfem::DofToQuad::TENSOR
Tensor product representation using 1D matrices/tensors with dimensions using 1D number of quadrature...
Definition: fe.hpp:154
mfem::QuadratureSpace::GetElementIntRule
const IntegrationRule & GetElementIntRule(int idx) const
Get the IntegrationRule associated with mesh element idx.
Definition: fespace.hpp:703
mfem::QuadratureInterpolator::fespace
const FiniteElementSpace * fespace
Not owned.
Definition: quadinterpolator.hpp:41
geom
const Geometry::Type geom
Definition: ex1.cpp:40
mfem::DofToQuad::nqpt
int nqpt
Number of quadrature points. When mode is TENSOR, this is the 1D number.
Definition: fe.hpp:166
mfem::DofToQuad::ndof
int ndof
Number of degrees of freedom = number of basis functions. When mode is TENSOR, this is the 1D number...
Definition: fe.hpp:162
mfem::Mesh::GetNE
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:693
mfem::QuadratureInterpolator::DERIVATIVES
Evaluate the derivatives at quadrature points.
Definition: quadinterpolator.hpp:60
mfem::GeometricFactors
Structure for storing mesh geometric factors: coordinates, Jacobians, and determinants of the Jacobia...
Definition: mesh.hpp:1332
mfem::Vector::Write
double * Write(bool on_dev=true)
Shortcut for mfem::Write(vec.GetMemory(), vec.Size(), on_dev).
Definition: vector.hpp:341
mfem::Reshape
DeviceTensor< sizeof...(Dims), T > Reshape(T *ptr, Dims...dims)
Wrap a pointer as a DeviceTensor with automatically deduced template parameters.
Definition: dtensor.hpp:134
mfem::FiniteElementSpace::GetNE
int GetNE() const
Returns number of elements in the mesh.
Definition: fespace.hpp:408
mfem::Mesh::GetGeometricFactors
const GeometricFactors * GetGeometricFactors(const IntegrationRule &ir, const int flags)
Return the mesh geometric factors corresponding to the given integration rule.
Definition: mesh.cpp:763
mfem::QuadratureInterpolator::q_layout
QVectorLayout q_layout
Output Q-vector layout.
Definition: quadinterpolator.hpp:44
mfem::MAX_Q1D
const int MAX_Q1D
Definition: forall.hpp:36
mfem::FiniteElementSpace::GetMesh
Mesh * GetMesh() const
Returns the mesh.
Definition: fespace.hpp:295
b
double b
Definition: lissajous.cpp:42
mfem::QuadratureInterpolator::IntRule
const IntegrationRule * IntRule
Not owned.
Definition: quadinterpolator.hpp:43
mfem::Array::Read
const T * Read(bool on_dev=true) const
Shortcut for mfem::Read(a.GetMemory(), a.Size(), on_dev).
Definition: array.hpp:273
mfem::FiniteElementSpace::GetVDim
int GetVDim() const
Returns vector dimension.
Definition: fespace.hpp:370
mfem::Mesh::Dimension
int Dimension() const
Definition: mesh.hpp:744
mfem::QuadratureInterpolator::QuadratureInterpolator
QuadratureInterpolator(const FiniteElementSpace &fes, const IntegrationRule &ir)
Definition: quadinterpolator.cpp:20
mfem::FiniteElementSpace
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:87
mfem::QuadratureInterpolator::PhysDerivatives
void PhysDerivatives(const Vector &e_vec, Vector &q_der) const
Interpolate the derivatives in physical space of the E-vector e_vec at quadrature points...
Definition: quadinterpolator.cpp:1302
mfem::QuadratureInterpolator::Eval2D
static void Eval2D(const int NE, const int vdim, const DofToQuad &maps, const Vector &e_vec, Vector &q_val, Vector &q_der, Vector &q_det, const int eval_flags)
Template compute kernel for 2D.
Definition: quadinterpolator.cpp:51
mfem::DofToQuad
Structure representing the matrices/tensors needed to evaluate (in reference space) the values...
Definition: fe.hpp:128
mfem::DofToQuad::B
Array< double > B
Basis functions evaluated at quadrature points.
Definition: fe.hpp:177
mfem::FiniteElement::GetDofToQuad
virtual const DofToQuad & GetDofToQuad(const IntegrationRule &ir, DofToQuad::Mode mode) const
Definition: fe.cpp:370
mfem::QuadratureInterpolator::DETERMINANTS
Assuming the derivative at quadrature points form a matrix, this flag can be used to compute and stor...
Definition: quadinterpolator.hpp:64
mfem::QuadratureInterpolator::Eval3D
static void Eval3D(const int NE, const int vdim, const DofToQuad &maps, const Vector &e_vec, Vector &q_val, Vector &q_der, Vector &q_det, const int eval_flags)
Template compute kernel for 3D.
Definition: quadinterpolator.cpp:139
mfem::DofToQuad::Mode
Mode
Type of data stored in the arrays B, Bt, G, and Gt.
Definition: fe.hpp:141
mfem::DofToQuad::FULL
Full multidimensional representation which does not use tensor product structure. The ordering of the...
Definition: fe.hpp:146
mfem::Vector::Read
const double * Read(bool on_dev=true) const
Shortcut for mfem::Read(vec.GetMemory(), vec.Size(), on_dev).
Definition: vector.hpp:333
mfem::QuadratureInterpolator::Mult
void Mult(const Vector &e_vec, unsigned eval_flags, Vector &q_val, Vector &q_der, Vector &q_det) const
Interpolate the E-vector e_vec to quadrature points.
Definition: quadinterpolator.cpp:229
dim
int dim
Definition: ex24.cpp:43
mfem::MAX_D1D
const int MAX_D1D
Definition: forall.hpp:35
mfem::DofToQuad::G
Array< double > G
Gradients/divergences/curls of basis functions evaluated at quadrature points.
Definition: fe.hpp:198
mfem::FiniteElementSpace::GetFE
const FiniteElement * GetFE(int i) const
Returns pointer to the FiniteElement associated with i&#39;th element.
Definition: fespace.cpp:1642
mfem::QuadratureInterpolator::Derivatives
void Derivatives(const Vector &e_vec, Vector &q_der) const
Interpolate the derivatives of the E-vector e_vec at quadrature points.
Definition: quadinterpolator.cpp:1284
mfem::Vector
Vector data type.
Definition: vector.hpp:48
mfem::QuadratureSpace
Class representing the storage layout of a QuadratureFunction.
Definition: fespace.hpp:671
mfem::QuadratureInterpolator::Values
void Values(const Vector &e_vec, Vector &q_val) const
Interpolate the values of the E-vector e_vec at quadrature points.
Definition: quadinterpolator.cpp:649
mfem::out
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
mfem::QuadratureInterpolator::MultTranspose
void MultTranspose(unsigned eval_flags, const Vector &q_val, const Vector &q_der, Vector &e_vec) const
Perform the transpose operation of Mult(). (TODO)
Definition: quadinterpolator.cpp:383
mfem::QVectorLayout::byVDIM
VDIM x NQPT x NE.
mfem::QuadratureInterpolator::VALUES
Evaluate the values at quadrature points.
Definition: quadinterpolator.hpp:59
mfem::QuadratureInterpolator::qspace
const QuadratureSpace * qspace
Not owned.
Definition: quadinterpolator.hpp:42
Generated on Fri Feb 18 2022 18:54:17 for MFEM by   doxygen 1.8.5