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