1 // Copyright (c) 2010-2021, 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 /// A structure used to pass additional data to f_build_conv and f_apply_conv
13 struct ConvectionContext {
14 CeedInt dim, space_dim, vdim;
15 CeedScalar coeff[3];
16 CeedScalar alpha;
17 };
18
19 /// libCEED Q-function for building quadrature data for a convection operator
20 /// with a constant coefficient
CEED_QFUNCTION(f_build_conv_const)21 CEED_QFUNCTION(f_build_conv_const)(void *ctx, CeedInt Q,
22 const CeedScalar *const *in,
23 CeedScalar *const *out)
24 {
25 ConvectionContext *bc = (ConvectionContext*)ctx;
26 // in[0] is Jacobians with shape [dim, nc=dim, Q]
27 // in[1] is quadrature weights, size (Q)
28 //
29 // At every quadrature point, compute and store qw * adj(J).
30 const CeedScalar coeff0 = bc->coeff[0];
31 const CeedScalar coeff1 = bc->coeff[1];
32 const CeedScalar coeff2 = bc->coeff[2];
33 const CeedScalar alpha = bc->alpha;
34 const CeedScalar *J = in[0], *qw = in[1];
35 CeedScalar *qd = out[0];
36 switch (bc->dim + 10 * bc->space_dim)
37 {
38 case 11:
39 for (CeedInt i = 0; i < Q; i++)
40 {
41 qd[i] = alpha * coeff0 * qw[i] * J[i];
42 }
43 break;
44 case 22:
45 for (CeedInt i = 0; i < Q; i++)
46 {
47 // J: 0 2 qd: 0 1 adj(J): J22 -J12
48 // 1 3 1 2 -J21 J11
49 const CeedScalar J11 = J[i + Q * 0];
50 const CeedScalar J21 = J[i + Q * 1];
51 const CeedScalar J12 = J[i + Q * 2];
52 const CeedScalar J22 = J[i + Q * 3];
53 const CeedScalar w = alpha * qw[i];
54 const CeedScalar wx = w * coeff0;
55 const CeedScalar wy = w * coeff1;
56 qd[i + Q * 0] = wx * J22 - wy * J12;
57 qd[i + Q * 1] = -wx * J21 + wy * J11;
58 }
59 break;
60 case 33:
61 for (CeedInt i = 0; i < Q; i++)
62 {
63 // J: 0 3 6 qd: 0 1 2
64 // 1 4 7 1 3 4
65 // 2 5 8 2 4 5
66 const CeedScalar J11 = J[i + Q * 0];
67 const CeedScalar J21 = J[i + Q * 1];
68 const CeedScalar J31 = J[i + Q * 2];
69 const CeedScalar J12 = J[i + Q * 3];
70 const CeedScalar J22 = J[i + Q * 4];
71 const CeedScalar J32 = J[i + Q * 5];
72 const CeedScalar J13 = J[i + Q * 6];
73 const CeedScalar J23 = J[i + Q * 7];
74 const CeedScalar J33 = J[i + Q * 8];
75 const CeedScalar A11 = J22 * J33 - J23 * J32;
76 const CeedScalar A12 = J13 * J32 - J12 * J33;
77 const CeedScalar A13 = J12 * J23 - J13 * J22;
78 const CeedScalar A21 = J23 * J31 - J21 * J33;
79 const CeedScalar A22 = J11 * J33 - J13 * J31;
80 const CeedScalar A23 = J13 * J21 - J11 * J23;
81 const CeedScalar A31 = J21 * J32 - J22 * J31;
82 const CeedScalar A32 = J12 * J31 - J11 * J32;
83 const CeedScalar A33 = J11 * J22 - J12 * J21;
84 const CeedScalar w = alpha * qw[i];
85 const CeedScalar wx = w * coeff0;
86 const CeedScalar wy = w * coeff1;
87 const CeedScalar wz = w * coeff2;
88 qd[i + Q * 0] = wx * A11 + wy * A12 + wz * A13;
89 qd[i + Q * 1] = wx * A21 + wy * A22 + wz * A23;
90 qd[i + Q * 2] = wx * A31 + wy * A32 + wz * A33;
91 }
92 break;
93 }
94 return 0;
95 }
96
97 /// libCEED Q-function for building quadrature data for a convection operator
98 /// coefficient evaluated at quadrature points.
CEED_QFUNCTION(f_build_conv_quad)99 CEED_QFUNCTION(f_build_conv_quad)(void *ctx, CeedInt Q,
100 const CeedScalar *const *in,
101 CeedScalar *const *out)
102 {
103 ConvectionContext *bc = (ConvectionContext *)ctx;
104 // in[1] is Jacobians with shape [dim, nc=dim, Q]
105 // in[2] is quadrature weights, size (Q)
106 //
107 // At every quadrature point, compute and store qw * adj(J).
108 const CeedScalar *c = in[0], *J = in[1], *qw = in[2];
109 const CeedScalar alpha = bc->alpha;
110 CeedScalar *qd = out[0];
111 switch (bc->dim + 10 * bc->space_dim)
112 {
113 case 11:
114 for (CeedInt i = 0; i < Q; i++)
115 {
116 const CeedScalar coeff = c[i];
117 qd[i] = alpha * coeff * qw[i] * J[i];
118 }
119 break;
120 case 22:
121 for (CeedInt i = 0; i < Q; i++)
122 {
123 // J: 0 2 qd: 0 1 adj(J): J22 -J12
124 // 1 3 1 2 -J21 J11
125 const CeedScalar J11 = J[i + Q * 0];
126 const CeedScalar J21 = J[i + Q * 1];
127 const CeedScalar J12 = J[i + Q * 2];
128 const CeedScalar J22 = J[i + Q * 3];
129 const CeedScalar w = alpha * qw[i];
130 const CeedScalar wx = w * c[i + Q * 0];
131 const CeedScalar wy = w * c[i + Q * 1];
132 qd[i + Q * 0] = wx * J22 - wy * J12;
133 qd[i + Q * 1] = -wx * J21 + wy * J11;
134 }
135 break;
136 case 33:
137 for (CeedInt i = 0; i < Q; i++)
138 {
139 // J: 0 3 6 qd: 0 1 2
140 // 1 4 7 1 3 4
141 // 2 5 8 2 4 5
142 const CeedScalar J11 = J[i + Q * 0];
143 const CeedScalar J21 = J[i + Q * 1];
144 const CeedScalar J31 = J[i + Q * 2];
145 const CeedScalar J12 = J[i + Q * 3];
146 const CeedScalar J22 = J[i + Q * 4];
147 const CeedScalar J32 = J[i + Q * 5];
148 const CeedScalar J13 = J[i + Q * 6];
149 const CeedScalar J23 = J[i + Q * 7];
150 const CeedScalar J33 = J[i + Q * 8];
151 const CeedScalar A11 = J22 * J33 - J23 * J32;
152 const CeedScalar A12 = J13 * J32 - J12 * J33;
153 const CeedScalar A13 = J12 * J23 - J13 * J22;
154 const CeedScalar A21 = J23 * J31 - J21 * J33;
155 const CeedScalar A22 = J11 * J33 - J13 * J31;
156 const CeedScalar A23 = J13 * J21 - J11 * J23;
157 const CeedScalar A31 = J21 * J32 - J22 * J31;
158 const CeedScalar A32 = J12 * J31 - J11 * J32;
159 const CeedScalar A33 = J11 * J22 - J12 * J21;
160 const CeedScalar w = alpha * qw[i];
161 const CeedScalar wx = w * c[i + Q * 0];
162 const CeedScalar wy = w * c[i + Q * 1];
163 const CeedScalar wz = w * c[i + Q * 2];
164 qd[i + Q * 0] = wx * A11 + wy * A12 + wz * A13;
165 qd[i + Q * 1] = wx * A21 + wy * A22 + wz * A23;
166 qd[i + Q * 2] = wx * A31 + wy * A32 + wz * A33;
167 }
168 break;
169 }
170 return 0;
171 }
172
173 /// libCEED Q-function for applying a conv operator
CEED_QFUNCTION(f_apply_conv)174 CEED_QFUNCTION(f_apply_conv)(void *ctx, CeedInt Q,
175 const CeedScalar *const *in,
176 CeedScalar *const *out)
177 {
178 ConvectionContext *bc = (ConvectionContext *)ctx;
179 // in[0], out[0] have shape [dim, nc=1, Q]
180 const CeedScalar *ug = in[0], *qd = in[1];
181 CeedScalar *vg = out[0];
182 switch (10*bc->dim + bc->vdim)
183 {
184 case 11:
185 for (CeedInt i = 0; i < Q; i++)
186 {
187 vg[i] = ug[i] * qd[i];
188 }
189 break;
190 case 21:
191 for (CeedInt i = 0; i < Q; i++)
192 {
193 const CeedScalar ug0 = ug[i + Q * 0];
194 const CeedScalar ug1 = ug[i + Q * 1];
195 vg[i] = qd[i + Q * 0] * ug0 + qd[i + Q * 1] * ug1;
196 }
197 break;
198 case 22:
199 for (CeedInt i = 0; i < Q; i++)
200 {
201 const CeedScalar qd0 = qd[i + Q * 0];
202 const CeedScalar qd1 = qd[i + Q * 1];
203 for (CeedInt c = 0; c < 2; c++)
204 {
205 const CeedScalar ug0 = ug[i + Q * (c+2*0)];
206 const CeedScalar ug1 = ug[i + Q * (c+2*1)];
207 vg[i + Q * c] = qd0 * ug0 + qd1 * ug1;
208 }
209 }
210 break;
211 case 31:
212 for (CeedInt i = 0; i < Q; i++)
213 {
214 const CeedScalar ug0 = ug[i + Q * 0];
215 const CeedScalar ug1 = ug[i + Q * 1];
216 const CeedScalar ug2 = ug[i + Q * 2];
217 vg[i] = qd[i + Q * 0] * ug0 + qd[i + Q * 1] * ug1 + qd[i + Q * 2] * ug2;
218 }
219 break;
220 case 33:
221 for (CeedInt i = 0; i < Q; i++)
222 {
223 const CeedScalar qd0 = qd[i + Q * 0];
224 const CeedScalar qd1 = qd[i + Q * 1];
225 const CeedScalar qd2 = qd[i + Q * 2];
226 for (CeedInt c = 0; c < 3; c++)
227 {
228 const CeedScalar ug0 = ug[i + Q * (c+3*0)];
229 const CeedScalar ug1 = ug[i + Q * (c+3*1)];
230 const CeedScalar ug2 = ug[i + Q * (c+3*2)];
231 vg[i + Q * c] = qd0 * ug0 + qd1 * ug1 + qd2 * ug2;
232 }
233 }
234 break;
235 }
236 return 0;
237 }
238
239 /// libCEED Q-function for applying a conv operator
CEED_QFUNCTION(f_apply_conv_mf_const)240 CEED_QFUNCTION(f_apply_conv_mf_const)(void *ctx, CeedInt Q,
241 const CeedScalar *const *in,
242 CeedScalar *const *out)
243 {
244 ConvectionContext *bc = (ConvectionContext*)ctx;
245 // in[0], out[0] have shape [dim, nc=1, Q]
246 // in[1] is Jacobians with shape [dim, nc=dim, Q]
247 // in[2] is quadrature weights, size (Q)
248 //
249 // At every quadrature point, compute qw * adj(J).
250 const CeedScalar coeff0 = bc->coeff[0];
251 const CeedScalar coeff1 = bc->coeff[1];
252 const CeedScalar coeff2 = bc->coeff[2];
253 const CeedScalar alpha = bc->alpha;
254 const CeedScalar *ug = in[0], *J = in[1], *qw = in[2];
255 CeedScalar *vg = out[0];
256 switch (10 * bc->dim + bc->vdim)
257 {
258 case 11:
259 for (CeedInt i = 0; i < Q; i++)
260 {
261 const CeedScalar qd = alpha * coeff0 * qw[i] * J[i];
262 vg[i] = ug[i] * qd;
263 }
264 break;
265 case 21:
266 for (CeedInt i = 0; i < Q; i++)
267 {
268 // J: 0 2 qd: 0 1 adj(J): J22 -J12
269 // 1 3 1 2 -J21 J11
270 const CeedScalar J11 = J[i + Q * 0];
271 const CeedScalar J21 = J[i + Q * 1];
272 const CeedScalar J12 = J[i + Q * 2];
273 const CeedScalar J22 = J[i + Q * 3];
274 const CeedScalar w = alpha * qw[i];
275 const CeedScalar wx = w * coeff0;
276 const CeedScalar wy = w * coeff1;
277 const CeedScalar qd0 = wx * J22 - wy * J12;
278 const CeedScalar qd1 = -wx * J21 + wy * J11;
279 const CeedScalar ug0 = ug[i + Q * 0];
280 const CeedScalar ug1 = ug[i + Q * 1];
281 vg[i] = qd0 * ug0 + qd1 * ug1;
282 }
283 break;
284 case 22:
285 for (CeedInt i = 0; i < Q; i++)
286 {
287 // J: 0 2 qd: 0 1 adj(J): J22 -J12
288 // 1 3 1 2 -J21 J11
289 const CeedScalar J11 = J[i + Q * 0];
290 const CeedScalar J21 = J[i + Q * 1];
291 const CeedScalar J12 = J[i + Q * 2];
292 const CeedScalar J22 = J[i + Q * 3];
293 const CeedScalar w = alpha * qw[i];
294 const CeedScalar wx = w * coeff0;
295 const CeedScalar wy = w * coeff1;
296 const CeedScalar qd0 = wx * J22 - wy * J12;
297 const CeedScalar qd1 = -wx * J21 + wy * J11;
298 for (CeedInt c = 0; c < 2; c++)
299 {
300 const CeedScalar ug0 = ug[i + Q * (c+2*0)];
301 const CeedScalar ug1 = ug[i + Q * (c+2*1)];
302 vg[i + Q * c] = qd0 * ug0 + qd1 * ug1;
303 }
304 }
305 break;
306 case 31:
307 for (CeedInt i = 0; i < Q; i++)
308 {
309 // J: 0 3 6 qd: 0 1 2
310 // 1 4 7 1 3 4
311 // 2 5 8 2 4 5
312 const CeedScalar J11 = J[i + Q * 0];
313 const CeedScalar J21 = J[i + Q * 1];
314 const CeedScalar J31 = J[i + Q * 2];
315 const CeedScalar J12 = J[i + Q * 3];
316 const CeedScalar J22 = J[i + Q * 4];
317 const CeedScalar J32 = J[i + Q * 5];
318 const CeedScalar J13 = J[i + Q * 6];
319 const CeedScalar J23 = J[i + Q * 7];
320 const CeedScalar J33 = J[i + Q * 8];
321 const CeedScalar A11 = J22 * J33 - J23 * J32;
322 const CeedScalar A12 = J13 * J32 - J12 * J33;
323 const CeedScalar A13 = J12 * J23 - J13 * J22;
324 const CeedScalar A21 = J23 * J31 - J21 * J33;
325 const CeedScalar A22 = J11 * J33 - J13 * J31;
326 const CeedScalar A23 = J13 * J21 - J11 * J23;
327 const CeedScalar A31 = J21 * J32 - J22 * J31;
328 const CeedScalar A32 = J12 * J31 - J11 * J32;
329 const CeedScalar A33 = J11 * J22 - J12 * J21;
330 const CeedScalar w = alpha * qw[i];
331 const CeedScalar wx = w * coeff0;
332 const CeedScalar wy = w * coeff1;
333 const CeedScalar wz = w * coeff2;
334 const CeedScalar qd0 = wx * A11 + wy * A12 + wz * A13;
335 const CeedScalar qd1 = wx * A21 + wy * A22 + wz * A23;
336 const CeedScalar qd2 = wx * A31 + wy * A32 + wz * A33;
337 const CeedScalar ug0 = ug[i + Q * 0];
338 const CeedScalar ug1 = ug[i + Q * 1];
339 const CeedScalar ug2 = ug[i + Q * 2];
340 vg[i] = qd0 * ug0 + qd1 * ug1 + qd2 * ug2;
341 }
342 break;
343 case 33:
344 for (CeedInt i = 0; i < Q; i++)
345 {
346 // J: 0 3 6 qd: 0 1 2
347 // 1 4 7 1 3 4
348 // 2 5 8 2 4 5
349 const CeedScalar J11 = J[i + Q * 0];
350 const CeedScalar J21 = J[i + Q * 1];
351 const CeedScalar J31 = J[i + Q * 2];
352 const CeedScalar J12 = J[i + Q * 3];
353 const CeedScalar J22 = J[i + Q * 4];
354 const CeedScalar J32 = J[i + Q * 5];
355 const CeedScalar J13 = J[i + Q * 6];
356 const CeedScalar J23 = J[i + Q * 7];
357 const CeedScalar J33 = J[i + Q * 8];
358 const CeedScalar A11 = J22 * J33 - J23 * J32;
359 const CeedScalar A12 = J13 * J32 - J12 * J33;
360 const CeedScalar A13 = J12 * J23 - J13 * J22;
361 const CeedScalar A21 = J23 * J31 - J21 * J33;
362 const CeedScalar A22 = J11 * J33 - J13 * J31;
363 const CeedScalar A23 = J13 * J21 - J11 * J23;
364 const CeedScalar A31 = J21 * J32 - J22 * J31;
365 const CeedScalar A32 = J12 * J31 - J11 * J32;
366 const CeedScalar A33 = J11 * J22 - J12 * J21;
367 const CeedScalar w = alpha * qw[i];
368 const CeedScalar wx = w * coeff0;
369 const CeedScalar wy = w * coeff1;
370 const CeedScalar wz = w * coeff2;
371 const CeedScalar qd0 = wx * A11 + wy * A12 + wz * A13;
372 const CeedScalar qd1 = wx * A21 + wy * A22 + wz * A23;
373 const CeedScalar qd2 = wx * A31 + wy * A32 + wz * A33;
374 for (CeedInt c = 0; c < 3; c++)
375 {
376 const CeedScalar ug0 = ug[i + Q * (c+3*0)];
377 const CeedScalar ug1 = ug[i + Q * (c+3*1)];
378 const CeedScalar ug2 = ug[i + Q * (c+3*2)];
379 vg[i + Q * c] = qd0 * ug0 + qd1 * ug1 + qd2 * ug2;
380 }
381 }
382 break;
383 }
384 return 0;
385 }
386
CEED_QFUNCTION(f_apply_conv_mf_quad)387 CEED_QFUNCTION(f_apply_conv_mf_quad)(void *ctx, CeedInt Q,
388 const CeedScalar *const *in,
389 CeedScalar *const *out)
390 {
391 ConvectionContext *bc = (ConvectionContext*)ctx;
392 // in[0], out[0] have shape [dim, nc=1, Q]
393 // in[1] is Jacobians with shape [dim, nc=dim, Q]
394 // in[2] is quadrature weights, size (Q)
395 //
396 // At every quadrature point, compute qw * adj(J).
397 const CeedScalar *c = in[0], *ug = in[1], *J = in[2], *qw = in[3];
398 const CeedScalar alpha = bc->alpha;
399 CeedScalar *vg = out[0];
400 switch (10 * bc->dim + bc->vdim)
401 {
402 case 11:
403 for (CeedInt i = 0; i < Q; i++)
404 {
405 const CeedScalar qd = alpha * c[i] * qw[i] * J[i];
406 vg[i] = ug[i] * qd;
407 }
408 break;
409 case 21:
410 for (CeedInt i = 0; i < Q; i++)
411 {
412 // J: 0 2 qd: 0 1 adj(J): J22 -J12
413 // 1 3 1 2 -J21 J11
414 const CeedScalar J11 = J[i + Q * 0];
415 const CeedScalar J21 = J[i + Q * 1];
416 const CeedScalar J12 = J[i + Q * 2];
417 const CeedScalar J22 = J[i + Q * 3];
418 const CeedScalar w = alpha * qw[i];
419 const CeedScalar wx = w * c[i + Q * 0];
420 const CeedScalar wy = w * c[i + Q * 1];
421 const CeedScalar qd0 = wx * J22 - wy * J12;
422 const CeedScalar qd1 = -wx * J21 + wy * J11;
423 const CeedScalar ug0 = ug[i + Q * 0];
424 const CeedScalar ug1 = ug[i + Q * 1];
425 vg[i] = qd0 * ug0 + qd1 * ug1;
426 }
427 break;
428 case 22:
429 for (CeedInt i = 0; i < Q; i++)
430 {
431 // J: 0 2 qd: 0 1 adj(J): J22 -J12
432 // 1 3 1 2 -J21 J11
433 const CeedScalar J11 = J[i + Q * 0];
434 const CeedScalar J21 = J[i + Q * 1];
435 const CeedScalar J12 = J[i + Q * 2];
436 const CeedScalar J22 = J[i + Q * 3];
437 const CeedScalar w = alpha * qw[i];
438 const CeedScalar wx = w * c[i + Q * 0];
439 const CeedScalar wy = w * c[i + Q * 1];
440 const CeedScalar qd0 = wx * J22 - wy * J12;
441 const CeedScalar qd1 = -wx * J21 + wy * J11;
442 for (CeedInt c = 0; c < 2; c++)
443 {
444 const CeedScalar ug0 = ug[i + Q * (c+2*0)];
445 const CeedScalar ug1 = ug[i + Q * (c+2*1)];
446 vg[i + Q * c] = qd0 * ug0 + qd1 * ug1;
447 }
448 }
449 break;
450 case 31:
451 for (CeedInt i = 0; i < Q; i++)
452 {
453 // J: 0 3 6 qd: 0 1 2
454 // 1 4 7 1 3 4
455 // 2 5 8 2 4 5
456 const CeedScalar J11 = J[i + Q * 0];
457 const CeedScalar J21 = J[i + Q * 1];
458 const CeedScalar J31 = J[i + Q * 2];
459 const CeedScalar J12 = J[i + Q * 3];
460 const CeedScalar J22 = J[i + Q * 4];
461 const CeedScalar J32 = J[i + Q * 5];
462 const CeedScalar J13 = J[i + Q * 6];
463 const CeedScalar J23 = J[i + Q * 7];
464 const CeedScalar J33 = J[i + Q * 8];
465 const CeedScalar A11 = J22 * J33 - J23 * J32;
466 const CeedScalar A12 = J13 * J32 - J12 * J33;
467 const CeedScalar A13 = J12 * J23 - J13 * J22;
468 const CeedScalar A21 = J23 * J31 - J21 * J33;
469 const CeedScalar A22 = J11 * J33 - J13 * J31;
470 const CeedScalar A23 = J13 * J21 - J11 * J23;
471 const CeedScalar A31 = J21 * J32 - J22 * J31;
472 const CeedScalar A32 = J12 * J31 - J11 * J32;
473 const CeedScalar A33 = J11 * J22 - J12 * J21;
474 const CeedScalar w = alpha * qw[i];
475 const CeedScalar wx = w * c[i + Q * 0];
476 const CeedScalar wy = w * c[i + Q * 1];
477 const CeedScalar wz = w * c[i + Q * 2];
478 const CeedScalar qd0 = wx * A11 + wy * A12 + wz * A13;
479 const CeedScalar qd1 = wx * A21 + wy * A22 + wz * A23;
480 const CeedScalar qd2 = wx * A31 + wy * A32 + wz * A33;
481 const CeedScalar ug0 = ug[i + Q * 0];
482 const CeedScalar ug1 = ug[i + Q * 1];
483 const CeedScalar ug2 = ug[i + Q * 2];
484 vg[i] = qd0 * ug0 + qd1 * ug1 + qd2 * ug2;
485 }
486 break;
487 case 33:
488 for (CeedInt i = 0; i < Q; i++)
489 {
490 // J: 0 3 6 qd: 0 1 2
491 // 1 4 7 1 3 4
492 // 2 5 8 2 4 5
493 const CeedScalar J11 = J[i + Q * 0];
494 const CeedScalar J21 = J[i + Q * 1];
495 const CeedScalar J31 = J[i + Q * 2];
496 const CeedScalar J12 = J[i + Q * 3];
497 const CeedScalar J22 = J[i + Q * 4];
498 const CeedScalar J32 = J[i + Q * 5];
499 const CeedScalar J13 = J[i + Q * 6];
500 const CeedScalar J23 = J[i + Q * 7];
501 const CeedScalar J33 = J[i + Q * 8];
502 const CeedScalar A11 = J22 * J33 - J23 * J32;
503 const CeedScalar A12 = J13 * J32 - J12 * J33;
504 const CeedScalar A13 = J12 * J23 - J13 * J22;
505 const CeedScalar A21 = J23 * J31 - J21 * J33;
506 const CeedScalar A22 = J11 * J33 - J13 * J31;
507 const CeedScalar A23 = J13 * J21 - J11 * J23;
508 const CeedScalar A31 = J21 * J32 - J22 * J31;
509 const CeedScalar A32 = J12 * J31 - J11 * J32;
510 const CeedScalar A33 = J11 * J22 - J12 * J21;
511 const CeedScalar w = alpha * qw[i];
512 const CeedScalar wx = w * c[i + Q * 0];
513 const CeedScalar wy = w * c[i + Q * 1];
514 const CeedScalar wz = w * c[i + Q * 2];
515 const CeedScalar qd0 = wx * A11 + wy * A12 + wz * A13;
516 const CeedScalar qd1 = wx * A21 + wy * A22 + wz * A23;
517 const CeedScalar qd2 = wx * A31 + wy * A32 + wz * A33;
518 for (CeedInt c = 0; c < 3; c++)
519 {
520 const CeedScalar ug0 = ug[i + Q * (c+3*0)];
521 const CeedScalar ug1 = ug[i + Q * (c+3*1)];
522 const CeedScalar ug2 = ug[i + Q * (c+3*2)];
523 vg[i + Q * c] = qd0 * ug0 + qd1 * ug1 + qd2 * ug2;
524 }
525 }
526 break;
527 }
528 return 0;
529 }
530