1 // Copyright (c) 2017-2018, Lawrence Livermore National Security, LLC. 2 // Produced at the Lawrence Livermore National Laboratory. LLNL-CODE-734707. 3 // All Rights reserved. See files LICENSE and NOTICE for details. 4 // 5 // This file is part of CEED, a collection of benchmarks, miniapps, software 6 // libraries and APIs for efficient high-order finite element and spectral 7 // element discretizations for exascale applications. For more information and 8 // source code availability see http://github.com/ceed. 9 // 10 // The CEED research is supported by the Exascale Computing Project 17-SC-20-SC, 11 // a collaborative effort of two U.S. Department of Energy organizations (Office 12 // of Science and the National Nuclear Security Administration) responsible for 13 // the planning and preparation of a capable exascale ecosystem, including 14 // software, applications, hardware, advanced system engineering and early 15 // testbed platforms, in support of the nation's exascale computing imperative. 16 17 #ifndef ex2_surface_h 18 #define ex2_surface_h 19 #include <ceed.h> 20 21 /// A structure used to pass additional data to f_build_diff 22 struct BuildContext { CeedInt dim, space_dim; }; 23 24 /// libCEED Q-function for building quadrature data for a diffusion operator 25 CEED_QFUNCTION(f_build_diff)(void *ctx, const CeedInt Q, 26 const CeedScalar *const *in, CeedScalar *const *out) { 27 struct BuildContext *bc = (struct BuildContext *)ctx; 28 // in[0] is Jacobians with shape [dim, nc=dim, Q] 29 // in[1] is quadrature weights, size (Q) 30 // 31 // At every quadrature point, compute w/det(J).adj(J).adj(J)^T and store 32 // the symmetric part of the result. 33 const CeedScalar *J = in[0], *w = in[1]; 34 CeedScalar *qdata = out[0]; 35 36 switch (bc->dim + 10*bc->space_dim) { 37 case 11: 38 CeedPragmaSIMD 39 for (CeedInt i=0; i<Q; i++) { 40 qdata[i] = w[i] / J[i]; 41 } // End of Quadrature Point Loop 42 break; 43 case 22: 44 CeedPragmaSIMD 45 for (CeedInt i=0; i<Q; i++) { 46 // J: 0 2 qdata: 0 2 adj(J): J22 -J12 47 // 1 3 2 1 -J21 J11 48 const CeedScalar J11 = J[i+Q*0]; 49 const CeedScalar J21 = J[i+Q*1]; 50 const CeedScalar J12 = J[i+Q*2]; 51 const CeedScalar J22 = J[i+Q*3]; 52 const CeedScalar qw = w[i] / (J11*J22 - J21*J12); 53 qdata[i+Q*0] = qw * (J12*J12 + J22*J22); 54 qdata[i+Q*1] = qw * (J11*J11 + J21*J21); 55 qdata[i+Q*2] = - qw * (J11*J12 + J21*J22); 56 } // End of Quadrature Point Loop 57 break; 58 case 33: 59 CeedPragmaSIMD 60 for (CeedInt i=0; i<Q; i++) { 61 // Compute the adjoint 62 CeedScalar A[3][3]; 63 for (CeedInt j=0; j<3; j++) 64 for (CeedInt k=0; k<3; k++) 65 // Equivalent code with J as a VLA and no mod operations: 66 // A[k][j] = J[j+1][k+1]*J[j+2][k+2] - J[j+1][k+2]*J[j+2][k+1] 67 A[k][j] = J[i+Q*((j+1)%3+3*((k+1)%3))]*J[i+Q*((j+2)%3+3*((k+2)%3))] - 68 J[i+Q*((j+1)%3+3*((k+2)%3))]*J[i+Q*((j+2)%3+3*((k+1)%3))]; 69 70 // Compute quadrature weight / det(J) 71 const CeedScalar qw = w[i] / (J[i+Q*0]*A[0][0] + J[i+Q*1]*A[1][1] + 72 J[i+Q*2]*A[2][2]); 73 74 // Compute geometric factors 75 // Stored in Voigt convention 76 // 0 5 4 77 // 5 1 3 78 // 4 3 2 79 qdata[i+Q*0] = qw * (A[0][0]*A[0][0] + A[0][1]*A[0][1] + A[0][2]*A[0][2]); 80 qdata[i+Q*1] = qw * (A[1][0]*A[1][0] + A[1][1]*A[1][1] + A[1][2]*A[1][2]); 81 qdata[i+Q*2] = qw * (A[2][0]*A[2][0] + A[2][1]*A[2][1] + A[2][2]*A[2][2]); 82 qdata[i+Q*3] = qw * (A[1][0]*A[2][0] + A[1][1]*A[2][1] + A[1][2]*A[2][2]); 83 qdata[i+Q*4] = qw * (A[0][0]*A[2][0] + A[0][1]*A[2][1] + A[0][2]*A[2][2]); 84 qdata[i+Q*5] = qw * (A[0][0]*A[1][0] + A[0][1]*A[1][1] + A[0][2]*A[1][2]); 85 } // End of Quadrature Point Loop 86 break; 87 } 88 return 0; 89 } 90 91 /// libCEED Q-function for applying a diff operator 92 CEED_QFUNCTION(f_apply_diff)(void *ctx, const CeedInt Q, 93 const CeedScalar *const *in, CeedScalar *const *out) { 94 struct BuildContext *bc = (struct BuildContext *)ctx; 95 // in[0], out[0] have shape [dim, nc=1, Q] 96 const CeedScalar *ug = in[0], *qdata = in[1]; 97 CeedScalar *vg = out[0]; 98 99 switch (bc->dim) { 100 case 1: 101 CeedPragmaSIMD 102 for (CeedInt i=0; i<Q; i++) { 103 vg[i] = ug[i] * qdata[i]; 104 } // End of Quadrature Point Loop 105 break; 106 case 2: 107 CeedPragmaSIMD 108 for (CeedInt i=0; i<Q; i++) { 109 // Read spatial derivatives of u 110 const CeedScalar du[2] = {ug[i+Q*0], 111 ug[i+Q*1] 112 }; 113 114 // Read qdata (dXdxdXdxT symmetric matrix) 115 // Stored in Voigt convention 116 // 0 2 117 // 2 1 118 // *INDENT-OFF* 119 const CeedScalar dXdxdXdxT[2][2] = {{qdata[i+0*Q], 120 qdata[i+2*Q]}, 121 {qdata[i+2*Q], 122 qdata[i+1*Q]}}; 123 // *INDENT-ON* 124 // j = direction of vg 125 for (int j=0; j<2; j++) 126 vg[i+j*Q] = (du[0] * dXdxdXdxT[0][j] + 127 du[1] * dXdxdXdxT[1][j]); 128 } // End of Quadrature Point Loop 129 break; 130 case 3: 131 CeedPragmaSIMD 132 for (CeedInt i=0; i<Q; i++) { 133 // Read spatial derivatives of u 134 const CeedScalar du[3] = {ug[i+Q*0], 135 ug[i+Q*1], 136 ug[i+Q*2] 137 }; 138 139 // Read qdata (dXdxdXdxT symmetric matrix) 140 // Stored in Voigt convention 141 // 0 5 4 142 // 5 1 3 143 // 4 3 2 144 // *INDENT-OFF* 145 const CeedScalar dXdxdXdxT[3][3] = {{qdata[i+0*Q], 146 qdata[i+5*Q], 147 qdata[i+4*Q]}, 148 {qdata[i+5*Q], 149 qdata[i+1*Q], 150 qdata[i+3*Q]}, 151 {qdata[i+4*Q], 152 qdata[i+3*Q], 153 qdata[i+2*Q]} 154 }; 155 // *INDENT-ON* 156 // j = direction of vg 157 for (int j=0; j<3; j++) 158 vg[i+j*Q] = (du[0] * dXdxdXdxT[0][j] + 159 du[1] * dXdxdXdxT[1][j] + 160 du[2] * dXdxdXdxT[2][j]); 161 } // End of Quadrature Point Loop 162 break; 163 } 164 return 0; 165 } 166 167 #endif // ex2_surface_h 168