Zhenhao Zhou, Ruochun Zhang, Ganesh Arivoli et al.
We present the numerical methods and GPU-accelerated implementation underlying a Total Lagrangian finite element framework for finite-deformation flexible multibody dynamics, introduced in the companion paper [1]. The framework supports 10-node quadratic tetrahedral (T10) elements and ANCF beam and shell elements, with quadrature-based hyperelastic response (St. Venant-Kirchhoff and Mooney-Rivlin) and an optional Kelvin-Voigt viscous stress contribution. Time stepping employs a velocity-based implicit backward-Euler scheme, yielding a nonlinear residual in velocity that couples inertia, internal and external forces, and bilateral constraints. Constraints are enforced via an augmented Lagrangian method (ALM), structured as an outer loop alternating an inner velocity solve with a dual-ascent multiplier update. We introduce a two-stage GPU parallelization strategy for internal force and tangent stiffness evaluation, and provide two inner solvers: a first-order AdamW optimizer and a second-order Newton solver that assembles and factorizes a sparse global Hessian on the GPU using cuDSS. A fixed-sparsity matrix strategy eliminates repeated symbolic analysis and enables efficient numerical refactorization across Newton iterations. For collision detection, we present a GPU-native two-thread asynchronous algorithm operating on triangle soups, avoiding bounding-volume hierarchies entirely. Systematic scaling benchmarks across all three supported element types and six mesh resolutions show that the Newton solver achieves approximately one order of magnitude reduction in real-time factor relative to CPU baselines at the largest resolutions tested. The frictional contact model is validated against closed-form rigid-body predictions through quasi-static and dynamic impact unit tests.