13d8e8822SJeremy L Thompson // Copyright (c) 2017-2022, Lawrence Livermore National Security, LLC and other CEED contributors. 23d8e8822SJeremy L Thompson // All Rights Reserved. See the top-level LICENSE and NOTICE files for details. 35754ecacSJeremy L Thompson // 43d8e8822SJeremy L Thompson // SPDX-License-Identifier: BSD-2-Clause 55754ecacSJeremy L Thompson // 63d8e8822SJeremy L Thompson // This file is part of CEED: http://github.com/ceed 75754ecacSJeremy L Thompson 85754ecacSJeremy L Thompson /// @file 95754ecacSJeremy L Thompson /// Geometric factors for solid mechanics example using PETSc 105754ecacSJeremy L Thompson 115754ecacSJeremy L Thompson #ifndef TRACTION_BOUNDARY_H 125754ecacSJeremy L Thompson #define TRACTION_BOUNDARY_H 135754ecacSJeremy L Thompson 14*c9c2c079SJeremy L Thompson #include <ceed.h> 15*c9c2c079SJeremy L Thompson 165754ecacSJeremy L Thompson // ----------------------------------------------------------------------------- 175754ecacSJeremy L Thompson // This QFunction computes the surface integral of the user traction vector on 185754ecacSJeremy L Thompson // the constrained faces. 195754ecacSJeremy L Thompson // 205754ecacSJeremy L Thompson // Reference (parent) 2D coordinates: X 215754ecacSJeremy L Thompson // Physical (current) 3D coordinates: x 225754ecacSJeremy L Thompson // Change of coordinate matrix: 235754ecacSJeremy L Thompson // dxdX_{i,j} = dx_i/dX_j (indicial notation) [3 * 2] 245754ecacSJeremy L Thompson // 255754ecacSJeremy L Thompson // (J1,J2,J3) is given by the cross product of the columns of dxdX_{i,j} 265754ecacSJeremy L Thompson // 275754ecacSJeremy L Thompson // detJb is the magnitude of (J1,J2,J3) 285754ecacSJeremy L Thompson // 295754ecacSJeremy L Thompson // Computed: 305754ecacSJeremy L Thompson // t * (w detJb) 315754ecacSJeremy L Thompson // 325754ecacSJeremy L Thompson // ----------------------------------------------------------------------------- 335754ecacSJeremy L Thompson CEED_QFUNCTION(SetupTractionBCs)(void *ctx, CeedInt Q, 345754ecacSJeremy L Thompson const CeedScalar *const *in, CeedScalar *const *out) { 355754ecacSJeremy L Thompson // *INDENT-OFF* 365754ecacSJeremy L Thompson // Inputs 375754ecacSJeremy L Thompson const CeedScalar(*J)[3][CEED_Q_VLA] = (const CeedScalar(*)[3][CEED_Q_VLA])in[0], 385754ecacSJeremy L Thompson (*w) = in[1]; 395754ecacSJeremy L Thompson // Outputs 405754ecacSJeremy L Thompson CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 415754ecacSJeremy L Thompson // *INDENT-ON* 425754ecacSJeremy L Thompson 435754ecacSJeremy L Thompson // User stress tensor 445754ecacSJeremy L Thompson const CeedScalar (*traction) = (const CeedScalar(*))ctx; 455754ecacSJeremy L Thompson 465754ecacSJeremy L Thompson CeedPragmaSIMD 475754ecacSJeremy L Thompson // Quadrature Point Loop 485754ecacSJeremy L Thompson for (CeedInt i = 0; i < Q; i++) { 495754ecacSJeremy L Thompson // Setup 505754ecacSJeremy L Thompson // *INDENT-OFF* 515754ecacSJeremy L Thompson const CeedScalar dxdX[3][2] = {{J[0][0][i], 525754ecacSJeremy L Thompson J[1][0][i]}, 535754ecacSJeremy L Thompson {J[0][1][i], 545754ecacSJeremy L Thompson J[1][1][i]}, 555754ecacSJeremy L Thompson {J[0][2][i], 565754ecacSJeremy L Thompson J[1][2][i]}}; 575754ecacSJeremy L Thompson // *INDENT-ON* 585754ecacSJeremy L Thompson // J1, J2, and J3 are given by the cross product of the columns of dxdX 595754ecacSJeremy L Thompson const CeedScalar J1 = dxdX[1][0] * dxdX[2][1] - dxdX[2][0] * dxdX[1][1]; 605754ecacSJeremy L Thompson const CeedScalar J2 = dxdX[2][0] * dxdX[0][1] - dxdX[0][0] * dxdX[2][1]; 615754ecacSJeremy L Thompson const CeedScalar J3 = dxdX[0][0] * dxdX[1][1] - dxdX[1][0] * dxdX[0][1]; 625754ecacSJeremy L Thompson 635754ecacSJeremy L Thompson // Qdata 645754ecacSJeremy L Thompson // -- Interp-to-Interp q_data 655754ecacSJeremy L Thompson CeedScalar wdetJb = w[i] * sqrt(J1 * J1 + J2 * J2 + J3 * J3); 665754ecacSJeremy L Thompson 675754ecacSJeremy L Thompson // Traction surface integral 685754ecacSJeremy L Thompson for (CeedInt j = 0; j < 3; j++) 695754ecacSJeremy L Thompson v[j][i] = traction[j] * wdetJb; 705754ecacSJeremy L Thompson 715754ecacSJeremy L Thompson } // End of Quadrature Point Loop 725754ecacSJeremy L Thompson 735754ecacSJeremy L Thompson // Return 745754ecacSJeremy L Thompson return 0; 755754ecacSJeremy L Thompson } 765754ecacSJeremy L Thompson // ----------------------------------------------------------------------------- 775754ecacSJeremy L Thompson 785754ecacSJeremy L Thompson #endif // End of TRACTION_BOUNDARY_H 79