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