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