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 /// @file 9 /// libCEED QFunctions for diffusion operator example using PETSc 10 11 #ifndef bp3_h 12 #define bp3_h 13 14 #include <ceed.h> 15 #include <math.h> 16 17 // ----------------------------------------------------------------------------- 18 // This QFunction sets up the geometric factors required to apply the 19 // diffusion operator 20 // 21 // We require the product of the inverse of the Jacobian and its transpose to 22 // properly compute integrals of the form: int( gradv gradu) 23 // 24 // Determinant of Jacobian: 25 // detJ = J11*A11 + J21*A12 + J31*A13 26 // Jij = Jacobian entry ij 27 // Aij = Adjoint ij 28 // 29 // Inverse of Jacobian: 30 // Bij = Aij / detJ 31 // 32 // Product of Inverse and Transpose: 33 // BBij = sum( Bik Bkj ) 34 // 35 // Stored: w B^T B detJ = w A^T A / detJ 36 // Note: This matrix is symmetric, so we only store 6 distinct entries 37 // qd: 0 3 6 38 // 1 4 7 39 // 2 5 8 40 // ----------------------------------------------------------------------------- 41 CEED_QFUNCTION(SetupDiffGeo)(void *ctx, CeedInt Q, 42 const CeedScalar *const *in, 43 CeedScalar *const *out) { 44 // Inputs 45 const CeedScalar(*J)[3][CEED_Q_VLA] = (const CeedScalar(*)[3][CEED_Q_VLA])in[1]; 46 const CeedScalar(*w) = in[2]; // Note: *X = in[0] 47 // Outputs 48 CeedScalar(*qd) = out[0]; 49 50 const CeedInt dim = 3; 51 // Quadrature Point Loop 52 CeedPragmaSIMD 53 for (CeedInt i=0; i<Q; i++) { 54 // Setup 55 CeedScalar A[3][3]; 56 for (CeedInt j = 0; j < dim; j++) { 57 for (CeedInt k = 0; k < dim; k++) { 58 // Equivalent code with no mod operations: 59 // A[k][j] = J[k+1][j+1]*J[k+2][j+2] - J[k+1][j+2]*J[k+2][j+1] 60 A[k][j] = J[(k + 1) % dim][(j + 1) % dim][i] * J[(k + 2) % dim][(j + 2) % dim][i] - 61 J[(k + 1) % dim][(j + 2) % dim][i] * J[(k + 2) % dim][(j + 1) % dim][i]; 62 } 63 } 64 const CeedScalar detJ = J[0][0][i] * A[0][0] + J[0][1][i] * A[0][1] + J[0][2][i] * A[0][2]; 65 66 const CeedScalar qw = w[i] / detJ; 67 qd[i+Q*0] = qw * (A[0][0]*A[0][0] + A[0][1]*A[0][1] + A[0][2]*A[0][2]); 68 qd[i+Q*1] = qw * (A[0][0]*A[1][0] + A[0][1]*A[1][1] + A[0][2]*A[1][2]); 69 qd[i+Q*2] = qw * (A[0][0]*A[2][0] + A[0][1]*A[2][1] + A[0][2]*A[2][2]); 70 qd[i+Q*3] = qw * (A[1][0]*A[1][0] + A[1][1]*A[1][1] + A[1][2]*A[1][2]); 71 qd[i+Q*4] = qw * (A[1][0]*A[2][0] + A[1][1]*A[2][1] + A[1][2]*A[2][2]); 72 qd[i+Q*5] = qw * (A[2][0]*A[2][0] + A[2][1]*A[2][1] + A[2][2]*A[2][2]); 73 qd[i+Q*6] = w[i] * detJ; 74 } // End of Quadrature Point Loop 75 76 return 0; 77 } 78 79 // ----------------------------------------------------------------------------- 80 // This QFunction sets up the rhs and true solution for the problem 81 // ----------------------------------------------------------------------------- 82 CEED_QFUNCTION(SetupDiffRhs)(void *ctx, CeedInt Q, 83 const CeedScalar *const *in, 84 CeedScalar *const *out) { 85 #ifndef M_PI 86 # define M_PI 3.14159265358979323846 87 #endif 88 const CeedScalar *x = in[0], *w = in[1]; 89 CeedScalar *true_soln = out[0], *rhs = out[1]; 90 91 // Quadrature Point Loop 92 CeedPragmaSIMD 93 for (CeedInt i=0; i<Q; i++) { 94 const CeedScalar c[3] = { 0, 1., 2. }; 95 const CeedScalar k[3] = { 1., 2., 3. }; 96 97 true_soln[i] = sin(M_PI*(c[0] + k[0]*x[i+Q*0])) * 98 sin(M_PI*(c[1] + k[1]*x[i+Q*1])) * 99 sin(M_PI*(c[2] + k[2]*x[i+Q*2])); 100 101 rhs[i] = w[i+Q*6] * M_PI*M_PI * (k[0]*k[0] + k[1]*k[1] + k[2]*k[2]) * 102 true_soln[i]; 103 } // End of Quadrature Point Loop 104 105 return 0; 106 } 107 108 // ----------------------------------------------------------------------------- 109 // This QFunction applies the diffusion operator for a scalar field. 110 // 111 // Inputs: 112 // ug - Input vector gradient at quadrature points 113 // q_data - Geometric factors 114 // 115 // Output: 116 // vg - Output vector (test functions) gradient at quadrature points 117 // 118 // ----------------------------------------------------------------------------- 119 CEED_QFUNCTION(Diff)(void *ctx, CeedInt Q, 120 const CeedScalar *const *in, CeedScalar *const *out) { 121 const CeedScalar *ug = in[0], *q_data = in[1]; 122 CeedScalar *vg = out[0]; 123 124 // Quadrature Point Loop 125 CeedPragmaSIMD 126 for (CeedInt i=0; i<Q; i++) { 127 // Read spatial derivatives of u 128 const CeedScalar du[3] = {ug[i+Q*0], 129 ug[i+Q*1], 130 ug[i+Q*2] 131 }; 132 // Read q_data (dXdxdXdx_T symmetric matrix) 133 const CeedScalar dXdxdXdx_T[3][3] = {{q_data[i+0*Q], 134 q_data[i+1*Q], 135 q_data[i+2*Q]}, 136 {q_data[i+1*Q], 137 q_data[i+3*Q], 138 q_data[i+4*Q]}, 139 {q_data[i+2*Q], 140 q_data[i+4*Q], 141 q_data[i+5*Q]} 142 }; 143 144 for (int j=0; j<3; j++) // j = direction of vg 145 vg[i+j*Q] = (du[0] * dXdxdXdx_T[0][j] + 146 du[1] * dXdxdXdx_T[1][j] + 147 du[2] * dXdxdXdx_T[2][j]); 148 149 } // End of Quadrature Point Loop 150 return 0; 151 } 152 // ----------------------------------------------------------------------------- 153 154 #endif // bp3_h 155