1 // Copyright (c) 2017-2026, Lawrence Livermore National Security, LLC and other CEED contributors.
2 // All Rights Reserved. See the top-level LICENSE and NOTICE files for details.
3 //
4 // SPDX-License-Identifier: BSD-2-Clause
5 //
6 // This file is part of CEED: http://github.com/ceed
7
8 #include <ceed/types.h>
9
10 /// A structure used to pass additional data to f_build_diff and f_apply_diff
11 struct BuildContext {
12 CeedInt dim, space_dim;
13 };
14
15 /// libCEED Q-function for building quadrature data for a diffusion operator
f_build_diff(void * ctx,const CeedInt Q,const CeedScalar * const * in,CeedScalar * const * out)16 CEED_QFUNCTION(f_build_diff)(void *ctx, const CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
17 BuildContext *bc = (BuildContext *)ctx;
18 // in[0] is Jacobians with shape [dim, nc=dim, Q]
19 // in[1] is quadrature weights, size (Q)
20 //
21 // At every quadrature point, compute w/det(J).adj(J).adj(J)^T and store the symmetric part of the result.
22 const CeedScalar *J = in[0], *w = in[1];
23 CeedScalar *qdata = out[0];
24
25 switch (bc->dim + 10 * bc->space_dim) {
26 case 11:
27 // Quadrature Point Loop
28 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { qdata[i] = w[i] / J[i]; }
29 break;
30 case 22:
31 // Quadrature Point Loop
32 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) {
33 // J: 0 2 qdata: 0 2 adj(J): J22 -J12
34 // 1 3 2 1 -J21 J11
35 const CeedScalar J11 = J[i + Q * 0];
36 const CeedScalar J21 = J[i + Q * 1];
37 const CeedScalar J12 = J[i + Q * 2];
38 const CeedScalar J22 = J[i + Q * 3];
39 const CeedScalar qw = w[i] / (J11 * J22 - J21 * J12);
40 qdata[i + Q * 0] = qw * (J12 * J12 + J22 * J22);
41 qdata[i + Q * 1] = qw * (J11 * J11 + J21 * J21);
42 qdata[i + Q * 2] = -qw * (J11 * J12 + J21 * J22);
43 }
44 break;
45 case 33:
46 // Quadrature Point Loop
47 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) {
48 // J: 0 3 6 qdata: 0 5 4
49 // 1 4 7 5 1 3
50 // 2 5 8 4 3 2
51 const CeedScalar J11 = J[i + Q * 0];
52 const CeedScalar J21 = J[i + Q * 1];
53 const CeedScalar J31 = J[i + Q * 2];
54 const CeedScalar J12 = J[i + Q * 3];
55 const CeedScalar J22 = J[i + Q * 4];
56 const CeedScalar J32 = J[i + Q * 5];
57 const CeedScalar J13 = J[i + Q * 6];
58 const CeedScalar J23 = J[i + Q * 7];
59 const CeedScalar J33 = J[i + Q * 8];
60 const CeedScalar A11 = J22 * J33 - J23 * J32;
61 const CeedScalar A12 = J13 * J32 - J12 * J33;
62 const CeedScalar A13 = J12 * J23 - J13 * J22;
63 const CeedScalar A21 = J23 * J31 - J21 * J33;
64 const CeedScalar A22 = J11 * J33 - J13 * J31;
65 const CeedScalar A23 = J13 * J21 - J11 * J23;
66 const CeedScalar A31 = J21 * J32 - J22 * J31;
67 const CeedScalar A32 = J12 * J31 - J11 * J32;
68 const CeedScalar A33 = J11 * J22 - J12 * J21;
69 const CeedScalar qw = w[i] / (J11 * A11 + J21 * A12 + J31 * A13);
70 qdata[i + Q * 0] = qw * (A11 * A11 + A12 * A12 + A13 * A13);
71 qdata[i + Q * 1] = qw * (A21 * A21 + A22 * A22 + A23 * A23);
72 qdata[i + Q * 2] = qw * (A31 * A31 + A32 * A32 + A33 * A33);
73 qdata[i + Q * 3] = qw * (A21 * A31 + A22 * A32 + A23 * A33);
74 qdata[i + Q * 4] = qw * (A11 * A31 + A12 * A32 + A13 * A33);
75 qdata[i + Q * 5] = qw * (A11 * A21 + A12 * A22 + A13 * A23);
76 }
77 break;
78 }
79 return 0;
80 }
81
82 /// libCEED Q-function for applying a diff operator
f_apply_diff(void * ctx,const CeedInt Q,const CeedScalar * const * in,CeedScalar * const * out)83 CEED_QFUNCTION(f_apply_diff)(void *ctx, const CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
84 BuildContext *bc = (BuildContext *)ctx;
85 // in[0], out[0] have shape [dim, nc=1, Q]
86 const CeedScalar *ug = in[0], *qdata = in[1];
87 CeedScalar *vg = out[0];
88
89 switch (bc->dim) {
90 case 1:
91 // Quadrature Point Loop
92 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { vg[i] = ug[i] * qdata[i]; }
93 break;
94 case 2:
95 // Quadrature Point Loop
96 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) {
97 const CeedScalar ug0 = ug[i + Q * 0];
98 const CeedScalar ug1 = ug[i + Q * 1];
99 vg[i + Q * 0] = qdata[i + Q * 0] * ug0 + qdata[i + Q * 2] * ug1;
100 vg[i + Q * 1] = qdata[i + Q * 2] * ug0 + qdata[i + Q * 1] * ug1;
101 }
102 break;
103 case 3:
104 // Quadrature Point Loop
105 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) {
106 const CeedScalar ug0 = ug[i + Q * 0];
107 const CeedScalar ug1 = ug[i + Q * 1];
108 const CeedScalar ug2 = ug[i + Q * 2];
109 vg[i + Q * 0] = qdata[i + Q * 0] * ug0 + qdata[i + Q * 5] * ug1 + qdata[i + Q * 4] * ug2;
110 vg[i + Q * 1] = qdata[i + Q * 5] * ug0 + qdata[i + Q * 1] * ug1 + qdata[i + Q * 3] * ug2;
111 vg[i + Q * 2] = qdata[i + Q * 4] * ug0 + qdata[i + Q * 3] * ug1 + qdata[i + Q * 2] * ug2;
112 }
113 break;
114 }
115 return 0;
116 }
117