1 // Copyright (c) 2017-2025, 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 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 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