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