1 // Copyright (c) 2017-2024, 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 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(build_diff)(void *ctx, const CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 17 struct BuildContext *build_data = (struct 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 22 // the symmetric part of the result. 23 const CeedScalar *J = in[0], *w = in[1]; 24 CeedScalar *q_data = out[0]; 25 26 switch (build_data->dim + 10 * build_data->space_dim) { 27 case 11: 28 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { q_data[i] = w[i] / J[i]; } // End of Quadrature Point Loop 29 break; 30 case 22: 31 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 32 // J: 0 2 q_data: 0 2 adj(J): J22 -J12 33 // 1 3 2 1 -J21 J11 34 const CeedScalar J11 = J[i + Q * 0]; 35 const CeedScalar J21 = J[i + Q * 1]; 36 const CeedScalar J12 = J[i + Q * 2]; 37 const CeedScalar J22 = J[i + Q * 3]; 38 const CeedScalar qw = w[i] / (J11 * J22 - J21 * J12); 39 q_data[i + Q * 0] = qw * (J12 * J12 + J22 * J22); 40 q_data[i + Q * 1] = qw * (J11 * J11 + J21 * J21); 41 q_data[i + Q * 2] = -qw * (J11 * J12 + J21 * J22); 42 } // End of Quadrature Point Loop 43 break; 44 case 33: 45 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 46 // Compute the adjoint 47 CeedScalar A[3][3]; 48 for (CeedInt j = 0; j < 3; j++) 49 for (CeedInt k = 0; k < 3; k++) 50 // Equivalent code with J as a VLA and no mod operations: 51 // A[k][j] = J[j+1][k+1]*J[j+2][k+2] - J[j+1][k+2]*J[j+2][k+1] 52 A[k][j] = J[i + Q * ((j + 1) % 3 + 3 * ((k + 1) % 3))] * J[i + Q * ((j + 2) % 3 + 3 * ((k + 2) % 3))] - 53 J[i + Q * ((j + 1) % 3 + 3 * ((k + 2) % 3))] * J[i + Q * ((j + 2) % 3 + 3 * ((k + 1) % 3))]; 54 55 // Compute quadrature weight / det(J) 56 const CeedScalar qw = w[i] / (J[i + Q * 0] * A[0][0] + J[i + Q * 1] * A[0][1] + J[i + Q * 2] * A[0][2]); 57 58 // Compute geometric factors 59 // Stored in Voigt convention 60 // 0 5 4 61 // 5 1 3 62 // 4 3 2 63 q_data[i + Q * 0] = qw * (A[0][0] * A[0][0] + A[0][1] * A[0][1] + A[0][2] * A[0][2]); 64 q_data[i + Q * 1] = qw * (A[1][0] * A[1][0] + A[1][1] * A[1][1] + A[1][2] * A[1][2]); 65 q_data[i + Q * 2] = qw * (A[2][0] * A[2][0] + A[2][1] * A[2][1] + A[2][2] * A[2][2]); 66 q_data[i + Q * 3] = qw * (A[1][0] * A[2][0] + A[1][1] * A[2][1] + A[1][2] * A[2][2]); 67 q_data[i + Q * 4] = qw * (A[0][0] * A[2][0] + A[0][1] * A[2][1] + A[0][2] * A[2][2]); 68 q_data[i + Q * 5] = qw * (A[0][0] * A[1][0] + A[0][1] * A[1][1] + A[0][2] * A[1][2]); 69 } // End of Quadrature Point Loop 70 break; 71 } 72 return 0; 73 } 74 75 /// libCEED Q-function for applying a diff operator 76 CEED_QFUNCTION(apply_diff)(void *ctx, const CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 77 struct BuildContext *build_data = (struct BuildContext *)ctx; 78 // in[0], out[0] have shape [dim, nc=1, Q] 79 const CeedScalar *ug = in[0], *q_data = in[1]; 80 CeedScalar *vg = out[0]; 81 82 switch (build_data->dim) { 83 case 1: 84 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { vg[i] = ug[i] * q_data[i]; } // End of Quadrature Point Loop 85 break; 86 case 2: 87 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 88 // Read spatial derivatives of u 89 const CeedScalar du[2] = {ug[i + Q * 0], ug[i + Q * 1]}; 90 91 // Read q_data (dXdxdXdx_T symmetric matrix) 92 // Stored in Voigt convention 93 // 0 2 94 // 2 1 95 const CeedScalar dXdxdXdx_T[2][2] = { 96 {q_data[i + 0 * Q], q_data[i + 2 * Q]}, 97 {q_data[i + 2 * Q], q_data[i + 1 * Q]} 98 }; 99 // j = direction of vg 100 for (int j = 0; j < 2; j++) vg[i + j * Q] = (du[0] * dXdxdXdx_T[0][j] + du[1] * dXdxdXdx_T[1][j]); 101 } // End of Quadrature Point Loop 102 break; 103 case 3: 104 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 105 // Read spatial derivatives of u 106 const CeedScalar du[3] = {ug[i + Q * 0], ug[i + Q * 1], ug[i + Q * 2]}; 107 108 // Read q_data (dXdxdXdx_T symmetric matrix) 109 // Stored in Voigt convention 110 // 0 5 4 111 // 5 1 3 112 // 4 3 2 113 const CeedScalar dXdxdXdx_T[3][3] = { 114 {q_data[i + 0 * Q], q_data[i + 5 * Q], q_data[i + 4 * Q]}, 115 {q_data[i + 5 * Q], q_data[i + 1 * Q], q_data[i + 3 * Q]}, 116 {q_data[i + 4 * Q], q_data[i + 3 * Q], q_data[i + 2 * Q]} 117 }; 118 // j = direction of vg 119 for (int j = 0; j < 3; j++) vg[i + j * Q] = (du[0] * dXdxdXdx_T[0][j] + du[1] * dXdxdXdx_T[1][j] + du[2] * dXdxdXdx_T[2][j]); 120 } // End of Quadrature Point Loop 121 break; 122 } 123 return 0; 124 } 125