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