1 // Copyright (c) 2017, Lawrence Livermore National Security, LLC. Produced at 2 // the Lawrence Livermore National Laboratory. LLNL-CODE-734707. All Rights 3 // reserved. See files LICENSE and NOTICE for details. 4 // 5 // This file is part of CEED, a collection of benchmarks, miniapps, software 6 // libraries and APIs for efficient high-order finite element and spectral 7 // element discretizations for exascale applications. For more information and 8 // source code availability see http://github.com/ceed. 9 // 10 // The CEED research is supported by the Exascale Computing Project 17-SC-20-SC, 11 // a collaborative effort of two U.S. Department of Energy organizations (Office 12 // of Science and the National Nuclear Security Administration) responsible for 13 // the planning and preparation of a capable exascale ecosystem, including 14 // software, applications, hardware, advanced system engineering and early 15 // testbed platforms, in support of the nation's exascale computing imperative. 16 17 /// @file 18 /// Geometric factors for solid mechanics example using PETSc 19 20 #ifndef TRACTION_BOUNDARY_H 21 #define TRACTION_BOUNDARY_H 22 23 // ----------------------------------------------------------------------------- 24 // This QFunction computes the surface integral of the user traction vector on 25 // the constrained faces. 26 // 27 // Reference (parent) 2D coordinates: X 28 // Physical (current) 3D coordinates: x 29 // Change of coordinate matrix: 30 // dxdX_{i,j} = dx_i/dX_j (indicial notation) [3 * 2] 31 // 32 // (J1,J2,J3) is given by the cross product of the columns of dxdX_{i,j} 33 // 34 // detJb is the magnitude of (J1,J2,J3) 35 // 36 // Computed: 37 // t * (w detJb) 38 // 39 // ----------------------------------------------------------------------------- 40 CEED_QFUNCTION(SetupTractionBCs)(void *ctx, CeedInt Q, 41 const CeedScalar *const *in, CeedScalar *const *out) { 42 // *INDENT-OFF* 43 // Inputs 44 const CeedScalar(*J)[3][CEED_Q_VLA] = (const CeedScalar(*)[3][CEED_Q_VLA])in[0], 45 (*w) = in[1]; 46 // Outputs 47 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 48 // *INDENT-ON* 49 50 // User stress tensor 51 const CeedScalar (*traction) = (const CeedScalar(*))ctx; 52 53 CeedPragmaSIMD 54 // Quadrature Point Loop 55 for (CeedInt i = 0; i < Q; i++) { 56 // Setup 57 // *INDENT-OFF* 58 const CeedScalar dxdX[3][2] = {{J[0][0][i], 59 J[1][0][i]}, 60 {J[0][1][i], 61 J[1][1][i]}, 62 {J[0][2][i], 63 J[1][2][i]}}; 64 // *INDENT-ON* 65 // J1, J2, and J3 are given by the cross product of the columns of dxdX 66 const CeedScalar J1 = dxdX[1][0] * dxdX[2][1] - dxdX[2][0] * dxdX[1][1]; 67 const CeedScalar J2 = dxdX[2][0] * dxdX[0][1] - dxdX[0][0] * dxdX[2][1]; 68 const CeedScalar J3 = dxdX[0][0] * dxdX[1][1] - dxdX[1][0] * dxdX[0][1]; 69 70 // Qdata 71 // -- Interp-to-Interp q_data 72 CeedScalar wdetJb = w[i] * sqrt(J1 * J1 + J2 * J2 + J3 * J3); 73 74 // Traction surface integral 75 for (CeedInt j = 0; j < 3; j++) 76 v[j][i] = traction[j] * wdetJb; 77 78 } // End of Quadrature Point Loop 79 80 // Return 81 return 0; 82 } 83 // ----------------------------------------------------------------------------- 84 85 #endif // End of TRACTION_BOUNDARY_H 86