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 bp4_h 12 #define bp4_h 13 14 #include <ceed.h> 15 #include <math.h> 16 17 // ----------------------------------------------------------------------------- 18 // This QFunction sets up the rhs and true solution for the problem 19 // ----------------------------------------------------------------------------- 20 CEED_QFUNCTION(SetupDiffRhs3)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 21 #ifndef M_PI 22 #define M_PI 3.14159265358979323846 23 #endif 24 const CeedScalar *x = in[0], *w = in[1]; 25 CeedScalar *true_soln = out[0], *rhs = out[1]; 26 27 // Quadrature Point Loop 28 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 29 const CeedScalar c[3] = {0, 1., 2.}; 30 const CeedScalar k[3] = {1., 2., 3.}; 31 32 // Component 1 33 true_soln[i + 0 * Q] = 34 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])); 35 // Component 2 36 true_soln[i + 1 * Q] = 2 * true_soln[i + 0 * Q]; 37 // Component 3 38 true_soln[i + 2 * Q] = 3 * true_soln[i + 0 * Q]; 39 40 // Component 1 41 rhs[i + 0 * Q] = w[i + Q * 0] * M_PI * M_PI * (k[0] * k[0] + k[1] * k[1] + k[2] * k[2]) * true_soln[i + 0 * Q]; 42 // Component 2 43 rhs[i + 1 * Q] = 2 * rhs[i + 0 * Q]; 44 // Component 3 45 rhs[i + 2 * Q] = 3 * rhs[i + 0 * Q]; 46 } // End of Quadrature Point Loop 47 48 return 0; 49 } 50 51 // ----------------------------------------------------------------------------- 52 // This QFunction applies the diffusion operator for a vector field of 3 components. 53 // 54 // Inputs: 55 // ug - Input vector Jacobian at quadrature points 56 // q_data - Geometric factors 57 // 58 // Output: 59 // vJ - Output vector (test functions) Jacobian at quadrature points 60 // ----------------------------------------------------------------------------- 61 CEED_QFUNCTION(Diff3)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 62 const CeedScalar *ug = in[0], *qd = in[1]; 63 CeedScalar *vg = out[0]; 64 65 // Quadrature Point Loop 66 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 67 // Read spatial derivatives of u components 68 const CeedScalar uJ[3][3] = { 69 {ug[i + (0 + 0 * 3) * Q], ug[i + (0 + 1 * 3) * Q], ug[i + (0 + 2 * 3) * Q]}, 70 {ug[i + (1 + 0 * 3) * Q], ug[i + (1 + 1 * 3) * Q], ug[i + (1 + 2 * 3) * Q]}, 71 {ug[i + (2 + 0 * 3) * Q], ug[i + (2 + 1 * 3) * Q], ug[i + (2 + 2 * 3) * Q]} 72 }; 73 // Read q_data (dXdxdXdx_T symmetric matrix) 74 const CeedScalar dXdxdXdx_T[3][3] = { 75 {qd[i + 1 * Q], qd[i + 2 * Q], qd[i + 3 * Q]}, 76 {qd[i + 2 * Q], qd[i + 4 * Q], qd[i + 5 * Q]}, 77 {qd[i + 3 * Q], qd[i + 5 * Q], qd[i + 6 * Q]} 78 }; 79 80 for (int k = 0; k < 3; k++) { // k = component 81 for (int j = 0; j < 3; j++) { // j = direction of vg 82 vg[i + (k + j * 3) * Q] = (uJ[k][0] * dXdxdXdx_T[0][j] + uJ[k][1] * dXdxdXdx_T[1][j] + uJ[k][2] * dXdxdXdx_T[2][j]); 83 } 84 } 85 } // End of Quadrature Point Loop 86 87 return 0; 88 } 89 // ----------------------------------------------------------------------------- 90 91 #endif // bp4_h 92