1 // Copyright (c) 2017-2022, 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 #ifndef ex2_surface_h 9 #define ex2_surface_h 10 11 /// A structure used to pass additional data to f_build_diff 12 struct BuildContext { CeedInt dim, space_dim; }; 13 14 /// libCEED Q-function for building quadrature data for a diffusion operator 15 CEED_QFUNCTION(f_build_diff)(void *ctx, const CeedInt Q, 16 const CeedScalar *const *in, CeedScalar *const *out) { 17 struct BuildContext *bc = (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 (bc->dim + 10*bc->space_dim) { 27 case 11: 28 CeedPragmaSIMD 29 for (CeedInt i=0; i<Q; i++) { 30 q_data[i] = w[i] / J[i]; 31 } // End of Quadrature Point Loop 32 break; 33 case 22: 34 CeedPragmaSIMD 35 for (CeedInt i=0; i<Q; i++) { 36 // J: 0 2 q_data: 0 2 adj(J): J22 -J12 37 // 1 3 2 1 -J21 J11 38 const CeedScalar J11 = J[i+Q*0]; 39 const CeedScalar J21 = J[i+Q*1]; 40 const CeedScalar J12 = J[i+Q*2]; 41 const CeedScalar J22 = J[i+Q*3]; 42 const CeedScalar qw = w[i] / (J11*J22 - J21*J12); 43 q_data[i+Q*0] = qw * (J12*J12 + J22*J22); 44 q_data[i+Q*1] = qw * (J11*J11 + J21*J21); 45 q_data[i+Q*2] = - qw * (J11*J12 + J21*J22); 46 } // End of Quadrature Point Loop 47 break; 48 case 33: 49 CeedPragmaSIMD 50 for (CeedInt i=0; i<Q; i++) { 51 // Compute the adjoint 52 CeedScalar A[3][3]; 53 for (CeedInt j=0; j<3; j++) 54 for (CeedInt k=0; k<3; k++) 55 // Equivalent code with J as a VLA and no mod operations: 56 // A[k][j] = J[j+1][k+1]*J[j+2][k+2] - J[j+1][k+2]*J[j+2][k+1] 57 A[k][j] = J[i+Q*((j+1)%3+3*((k+1)%3))]*J[i+Q*((j+2)%3+3*((k+2)%3))] - 58 J[i+Q*((j+1)%3+3*((k+2)%3))]*J[i+Q*((j+2)%3+3*((k+1)%3))]; 59 60 // Compute quadrature weight / det(J) 61 const CeedScalar qw = w[i] / (J[i+Q*0]*A[0][0] + J[i+Q*1]*A[0][1] + 62 J[i+Q*2]*A[0][2]); 63 64 // Compute geometric factors 65 // Stored in Voigt convention 66 // 0 5 4 67 // 5 1 3 68 // 4 3 2 69 q_data[i+Q*0] = qw * (A[0][0]*A[0][0] + A[0][1]*A[0][1] + A[0][2]*A[0][2]); 70 q_data[i+Q*1] = qw * (A[1][0]*A[1][0] + A[1][1]*A[1][1] + A[1][2]*A[1][2]); 71 q_data[i+Q*2] = qw * (A[2][0]*A[2][0] + A[2][1]*A[2][1] + A[2][2]*A[2][2]); 72 q_data[i+Q*3] = qw * (A[1][0]*A[2][0] + A[1][1]*A[2][1] + A[1][2]*A[2][2]); 73 q_data[i+Q*4] = qw * (A[0][0]*A[2][0] + A[0][1]*A[2][1] + A[0][2]*A[2][2]); 74 q_data[i+Q*5] = qw * (A[0][0]*A[1][0] + A[0][1]*A[1][1] + A[0][2]*A[1][2]); 75 } // End of Quadrature Point Loop 76 break; 77 } 78 return 0; 79 } 80 81 /// libCEED Q-function for applying a diff operator 82 CEED_QFUNCTION(f_apply_diff)(void *ctx, const CeedInt Q, 83 const CeedScalar *const *in, CeedScalar *const *out) { 84 struct BuildContext *bc = (struct BuildContext *)ctx; 85 // in[0], out[0] have shape [dim, nc=1, Q] 86 const CeedScalar *ug = in[0], *q_data = in[1]; 87 CeedScalar *vg = out[0]; 88 89 switch (bc->dim) { 90 case 1: 91 CeedPragmaSIMD 92 for (CeedInt i=0; i<Q; i++) { 93 vg[i] = ug[i] * q_data[i]; 94 } // End of Quadrature Point Loop 95 break; 96 case 2: 97 CeedPragmaSIMD 98 for (CeedInt i=0; i<Q; i++) { 99 // Read spatial derivatives of u 100 const CeedScalar du[2] = {ug[i+Q*0], 101 ug[i+Q*1] 102 }; 103 104 // Read q_data (dXdxdXdx_T symmetric matrix) 105 // Stored in Voigt convention 106 // 0 2 107 // 2 1 108 // *INDENT-OFF* 109 const CeedScalar dXdxdXdx_T[2][2] = {{q_data[i+0*Q], 110 q_data[i+2*Q]}, 111 {q_data[i+2*Q], 112 q_data[i+1*Q]}}; 113 // *INDENT-ON* 114 // j = direction of vg 115 for (int j=0; j<2; j++) 116 vg[i+j*Q] = (du[0] * dXdxdXdx_T[0][j] + 117 du[1] * dXdxdXdx_T[1][j]); 118 } // End of Quadrature Point Loop 119 break; 120 case 3: 121 CeedPragmaSIMD 122 for (CeedInt i=0; i<Q; i++) { 123 // Read spatial derivatives of u 124 const CeedScalar du[3] = {ug[i+Q*0], 125 ug[i+Q*1], 126 ug[i+Q*2] 127 }; 128 129 // Read q_data (dXdxdXdx_T symmetric matrix) 130 // Stored in Voigt convention 131 // 0 5 4 132 // 5 1 3 133 // 4 3 2 134 // *INDENT-OFF* 135 const CeedScalar dXdxdXdx_T[3][3] = {{q_data[i+0*Q], 136 q_data[i+5*Q], 137 q_data[i+4*Q]}, 138 {q_data[i+5*Q], 139 q_data[i+1*Q], 140 q_data[i+3*Q]}, 141 {q_data[i+4*Q], 142 q_data[i+3*Q], 143 q_data[i+2*Q]} 144 }; 145 // *INDENT-ON* 146 // j = direction of vg 147 for (int j=0; j<3; j++) 148 vg[i+j*Q] = (du[0] * dXdxdXdx_T[0][j] + 149 du[1] * dXdxdXdx_T[1][j] + 150 du[2] * dXdxdXdx_T[2][j]); 151 } // End of Quadrature Point Loop 152 break; 153 } 154 return 0; 155 } 156 157 #endif // ex2_surface_h 158