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 /// Geometric factors for solid mechanics example using PETSc 10 11 #ifndef TRACTION_BOUNDARY_H 12 #define TRACTION_BOUNDARY_H 13 14 #include <ceed.h> 15 16 // ----------------------------------------------------------------------------- 17 // This QFunction computes the surface integral of the user traction vector on 18 // the constrained faces. 19 // 20 // Reference (parent) 2D coordinates: X 21 // Physical (current) 3D coordinates: x 22 // Change of coordinate matrix: 23 // dxdX_{i,j} = dx_i/dX_j (indicial notation) [3 * 2] 24 // 25 // (J1,J2,J3) is given by the cross product of the columns of dxdX_{i,j} 26 // 27 // detJb is the magnitude of (J1,J2,J3) 28 // 29 // Computed: 30 // t * (w detJb) 31 // 32 // ----------------------------------------------------------------------------- 33 CEED_QFUNCTION(SetupTractionBCs)(void *ctx, CeedInt Q, 34 const CeedScalar *const *in, CeedScalar *const *out) { 35 // *INDENT-OFF* 36 // Inputs 37 const CeedScalar(*J)[3][CEED_Q_VLA] = (const CeedScalar(*)[3][CEED_Q_VLA])in[0], 38 (*w) = in[1]; 39 // Outputs 40 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 41 // *INDENT-ON* 42 43 // User stress tensor 44 const CeedScalar (*traction) = (const CeedScalar(*))ctx; 45 46 CeedPragmaSIMD 47 // Quadrature Point Loop 48 for (CeedInt i = 0; i < Q; i++) { 49 // Setup 50 // *INDENT-OFF* 51 const CeedScalar dxdX[3][2] = {{J[0][0][i], 52 J[1][0][i]}, 53 {J[0][1][i], 54 J[1][1][i]}, 55 {J[0][2][i], 56 J[1][2][i]}}; 57 // *INDENT-ON* 58 // J1, J2, and J3 are given by the cross product of the columns of dxdX 59 const CeedScalar J1 = dxdX[1][0] * dxdX[2][1] - dxdX[2][0] * dxdX[1][1]; 60 const CeedScalar J2 = dxdX[2][0] * dxdX[0][1] - dxdX[0][0] * dxdX[2][1]; 61 const CeedScalar J3 = dxdX[0][0] * dxdX[1][1] - dxdX[1][0] * dxdX[0][1]; 62 63 // Qdata 64 // -- Interp-to-Interp q_data 65 CeedScalar wdetJb = w[i] * sqrt(J1 * J1 + J2 * J2 + J3 * J3); 66 67 // Traction surface integral 68 for (CeedInt j = 0; j < 3; j++) 69 v[j][i] = traction[j] * wdetJb; 70 71 } // End of Quadrature Point Loop 72 73 // Return 74 return 0; 75 } 76 // ----------------------------------------------------------------------------- 77 78 #endif // End of TRACTION_BOUNDARY_H 79