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 25 // the symmetric part of the result. 26 const CeedScalar *J = in[0], *w = in[1]; 27 CeedScalar *qdata = out[0]; 28 29 switch (bc->dim + 10 * bc->space_dim) { 30 case 11: 31 // Quadrature Point Loop 32 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { qdata[i] = w[i] / J[i]; } 33 break; 34 case 22: 35 // Quadrature Point Loop 36 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 37 // J: 0 2 qdata: 0 2 adj(J): J22 -J12 38 // 1 3 2 1 -J21 J11 39 const CeedScalar J11 = J[i + Q * 0]; 40 const CeedScalar J21 = J[i + Q * 1]; 41 const CeedScalar J12 = J[i + Q * 2]; 42 const CeedScalar J22 = J[i + Q * 3]; 43 const CeedScalar qw = w[i] / (J11 * J22 - J21 * J12); 44 qdata[i + Q * 0] = qw * (J12 * J12 + J22 * J22); 45 qdata[i + Q * 1] = qw * (J11 * J11 + J21 * J21); 46 qdata[i + Q * 2] = -qw * (J11 * J12 + J21 * J22); 47 } 48 break; 49 case 33: 50 // Quadrature Point Loop 51 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 52 // J: 0 3 6 qdata: 0 5 4 53 // 1 4 7 5 1 3 54 // 2 5 8 4 3 2 55 const CeedScalar J11 = J[i + Q * 0]; 56 const CeedScalar J21 = J[i + Q * 1]; 57 const CeedScalar J31 = J[i + Q * 2]; 58 const CeedScalar J12 = J[i + Q * 3]; 59 const CeedScalar J22 = J[i + Q * 4]; 60 const CeedScalar J32 = J[i + Q * 5]; 61 const CeedScalar J13 = J[i + Q * 6]; 62 const CeedScalar J23 = J[i + Q * 7]; 63 const CeedScalar J33 = J[i + Q * 8]; 64 const CeedScalar A11 = J22 * J33 - J23 * J32; 65 const CeedScalar A12 = J13 * J32 - J12 * J33; 66 const CeedScalar A13 = J12 * J23 - J13 * J22; 67 const CeedScalar A21 = J23 * J31 - J21 * J33; 68 const CeedScalar A22 = J11 * J33 - J13 * J31; 69 const CeedScalar A23 = J13 * J21 - J11 * J23; 70 const CeedScalar A31 = J21 * J32 - J22 * J31; 71 const CeedScalar A32 = J12 * J31 - J11 * J32; 72 const CeedScalar A33 = J11 * J22 - J12 * J21; 73 const CeedScalar qw = w[i] / (J11 * A11 + J21 * A12 + J31 * A13); 74 qdata[i + Q * 0] = qw * (A11 * A11 + A12 * A12 + A13 * A13); 75 qdata[i + Q * 1] = qw * (A21 * A21 + A22 * A22 + A23 * A23); 76 qdata[i + Q * 2] = qw * (A31 * A31 + A32 * A32 + A33 * A33); 77 qdata[i + Q * 3] = qw * (A21 * A31 + A22 * A32 + A23 * A33); 78 qdata[i + Q * 4] = qw * (A11 * A31 + A12 * A32 + A13 * A33); 79 qdata[i + Q * 5] = qw * (A11 * A21 + A12 * A22 + A13 * A23); 80 } 81 break; 82 } 83 return 0; 84 } 85 86 /// libCEED Q-function for applying a diff operator 87 CEED_QFUNCTION(f_apply_diff)(void *ctx, const CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 88 BuildContext *bc = (BuildContext *)ctx; 89 // in[0], out[0] have shape [dim, nc=1, Q] 90 const CeedScalar *ug = in[0], *qdata = in[1]; 91 CeedScalar *vg = out[0]; 92 93 switch (bc->dim) { 94 case 1: 95 // Quadrature Point Loop 96 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { vg[i] = ug[i] * qdata[i]; } 97 break; 98 case 2: 99 // Quadrature Point Loop 100 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 101 const CeedScalar ug0 = ug[i + Q * 0]; 102 const CeedScalar ug1 = ug[i + Q * 1]; 103 vg[i + Q * 0] = qdata[i + Q * 0] * ug0 + qdata[i + Q * 2] * ug1; 104 vg[i + Q * 1] = qdata[i + Q * 2] * ug0 + qdata[i + Q * 1] * ug1; 105 } 106 break; 107 case 3: 108 // Quadrature Point Loop 109 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 110 const CeedScalar ug0 = ug[i + Q * 0]; 111 const CeedScalar ug1 = ug[i + Q * 1]; 112 const CeedScalar ug2 = ug[i + Q * 2]; 113 vg[i + Q * 0] = qdata[i + Q * 0] * ug0 + qdata[i + Q * 5] * ug1 + qdata[i + Q * 4] * ug2; 114 vg[i + Q * 1] = qdata[i + Q * 5] * ug0 + qdata[i + Q * 1] * ug1 + qdata[i + Q * 3] * ug2; 115 vg[i + Q * 2] = qdata[i + Q * 4] * ug0 + qdata[i + Q * 3] * ug1 + qdata[i + Q * 2] * ug2; 116 } 117 break; 118 } 119 return 0; 120 } 121 122 #endif // bp3_h 123