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 bp3_h 9 #define bp3_h 10 11 #include <ceed.h> 12 13 /// A structure used to pass additional data to f_build_diff and f_apply_diff 14 struct BuildContext { 15 CeedInt dim, space_dim; 16 }; 17 18 /// libCEED Q-function for building quadrature data for a diffusion operator 19 CEED_QFUNCTION(f_build_diff)(void *ctx, const CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 20 BuildContext *bc = (BuildContext *)ctx; 21 // in[0] is Jacobians with shape [dim, nc=dim, Q] 22 // in[1] is quadrature weights, size (Q) 23 // 24 // At every quadrature point, compute w/det(J).adj(J).adj(J)^T and store the symmetric part of the result. 25 const CeedScalar *J = in[0], *w = in[1]; 26 CeedScalar *qdata = out[0]; 27 28 switch (bc->dim + 10 * bc->space_dim) { 29 case 11: 30 // Quadrature Point Loop 31 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { qdata[i] = w[i] / J[i]; } 32 break; 33 case 22: 34 // Quadrature Point Loop 35 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 36 // J: 0 2 qdata: 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 qdata[i + Q * 0] = qw * (J12 * J12 + J22 * J22); 44 qdata[i + Q * 1] = qw * (J11 * J11 + J21 * J21); 45 qdata[i + Q * 2] = -qw * (J11 * J12 + J21 * J22); 46 } 47 break; 48 case 33: 49 // Quadrature Point Loop 50 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 51 // J: 0 3 6 qdata: 0 5 4 52 // 1 4 7 5 1 3 53 // 2 5 8 4 3 2 54 const CeedScalar J11 = J[i + Q * 0]; 55 const CeedScalar J21 = J[i + Q * 1]; 56 const CeedScalar J31 = J[i + Q * 2]; 57 const CeedScalar J12 = J[i + Q * 3]; 58 const CeedScalar J22 = J[i + Q * 4]; 59 const CeedScalar J32 = J[i + Q * 5]; 60 const CeedScalar J13 = J[i + Q * 6]; 61 const CeedScalar J23 = J[i + Q * 7]; 62 const CeedScalar J33 = J[i + Q * 8]; 63 const CeedScalar A11 = J22 * J33 - J23 * J32; 64 const CeedScalar A12 = J13 * J32 - J12 * J33; 65 const CeedScalar A13 = J12 * J23 - J13 * J22; 66 const CeedScalar A21 = J23 * J31 - J21 * J33; 67 const CeedScalar A22 = J11 * J33 - J13 * J31; 68 const CeedScalar A23 = J13 * J21 - J11 * J23; 69 const CeedScalar A31 = J21 * J32 - J22 * J31; 70 const CeedScalar A32 = J12 * J31 - J11 * J32; 71 const CeedScalar A33 = J11 * J22 - J12 * J21; 72 const CeedScalar qw = w[i] / (J11 * A11 + J21 * A12 + J31 * A13); 73 qdata[i + Q * 0] = qw * (A11 * A11 + A12 * A12 + A13 * A13); 74 qdata[i + Q * 1] = qw * (A21 * A21 + A22 * A22 + A23 * A23); 75 qdata[i + Q * 2] = qw * (A31 * A31 + A32 * A32 + A33 * A33); 76 qdata[i + Q * 3] = qw * (A21 * A31 + A22 * A32 + A23 * A33); 77 qdata[i + Q * 4] = qw * (A11 * A31 + A12 * A32 + A13 * A33); 78 qdata[i + Q * 5] = qw * (A11 * A21 + A12 * A22 + A13 * A23); 79 } 80 break; 81 } 82 return 0; 83 } 84 85 /// libCEED Q-function for applying a diff operator 86 CEED_QFUNCTION(f_apply_diff)(void *ctx, const CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 87 BuildContext *bc = (BuildContext *)ctx; 88 // in[0], out[0] have shape [dim, nc=1, Q] 89 const CeedScalar *ug = in[0], *qdata = in[1]; 90 CeedScalar *vg = out[0]; 91 92 switch (bc->dim) { 93 case 1: 94 // Quadrature Point Loop 95 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { vg[i] = ug[i] * qdata[i]; } 96 break; 97 case 2: 98 // Quadrature Point Loop 99 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 100 const CeedScalar ug0 = ug[i + Q * 0]; 101 const CeedScalar ug1 = ug[i + Q * 1]; 102 vg[i + Q * 0] = qdata[i + Q * 0] * ug0 + qdata[i + Q * 2] * ug1; 103 vg[i + Q * 1] = qdata[i + Q * 2] * ug0 + qdata[i + Q * 1] * ug1; 104 } 105 break; 106 case 3: 107 // Quadrature Point Loop 108 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 109 const CeedScalar ug0 = ug[i + Q * 0]; 110 const CeedScalar ug1 = ug[i + Q * 1]; 111 const CeedScalar ug2 = ug[i + Q * 2]; 112 vg[i + Q * 0] = qdata[i + Q * 0] * ug0 + qdata[i + Q * 5] * ug1 + qdata[i + Q * 4] * ug2; 113 vg[i + Q * 1] = qdata[i + Q * 5] * ug0 + qdata[i + Q * 1] * ug1 + qdata[i + Q * 3] * ug2; 114 vg[i + Q * 2] = qdata[i + Q * 4] * ug0 + qdata[i + Q * 3] * ug1 + qdata[i + Q * 2] * ug2; 115 } 116 break; 117 } 118 return 0; 119 } 120 121 #endif // bp3_h 122