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