Efficient simulation method of cosserat elastic rod based on strain-based admm

Abstract

The application discloses a high-efficiency simulation method of a Cosserat elastic rod based on strain ADMM and belongs to the technical field of physics-based animation and simulation in computer graphics. The application defines auxiliary variables in the strain level, and position-orientation coupling is explicitly exposed to the penalty term of the augmented Lagrangian function. Through two-stage mathematical relaxation of block coordinate descent (BCD) and projected dynamics (PD), the solving is replaced by an inexact proxy solving in the original variable updating of ADMM, and the local-global structure of PD pre-decomposability is restored. The application further augments the orientation block system matrix of bending-torsion constraints and introduces a robust penalty parameter clamping strategy. The application can stably and efficiently simulate slender structures such as trees and hair under real physical parameters without scene-by-scene parameter adjustment. Compared with the existing ADMM, in representative high-stiffness ratio scenes, the total running time is reduced by about 21% under the same number of iterations, and smaller numerical damping is shown under a limited iteration budget.

Description

Technical Field

[0001] This invention belongs to the field of physical-based animation and simulation technology in computer graphics, specifically involving a simulation method for Cosserat elastic rods (micro-polar elastic rods) based on the alternating direction multiplier method (ADMM) of strain. Background Technology

[0002] Slender, rod-like structures are ubiquitous in computer graphics, including hair, ropes, surgical sutures, cables, tree trunks, and various structural components. These slender rods simultaneously undergo stretching, shearing, bending, and torsional deformation, and their interactions produce complex natural phenomena such as trees swaying in the wind, hair flowing with body movement, and the movement of muscles and tendons.

[0003] The Cosserat elastic rod model equips each point on the rod with a material frame, accurately reproducing the aforementioned effects, including torsion and shear. While solving the Cosserat elastic rod system using a direct Newton solver yields accurate results, it requires reassembling and solving a new linear system at each iteration of every time step, resulting in enormous computational overhead. Projective Dynamics (PD), on the other hand, pre-decomposes the system matrix in the initial state and reuses this decomposition in all subsequent time steps, thus achieving fast solution at each time step. Within this framework, PD will... The elastic energy of the constraint is approximately a standard surrogate quadratic form:

[0004]

[0005] in For the local degrees of freedom involved in this constraint, For constant linear operators, To constrain the local projected target on the manifold; the projected target is calculated in the local step during the solution process. The system is solved in one step with a pre-decomposed constant system matrix, which is the local-global iteration of PD. The stiffness coefficient corresponding to the constraint is used to characterize the weight of the constraint in the total elastic energy.

[0006] However, existing Cosserat elastic rod simulators based on PD exhibit poor convergence when the ratio of tensile-shear stiffness to bending-torsional stiffness is large, forcing users to artificially reduce the material stiffness. Since the stiffness ratio is inversely proportional to the square of the cross-sectional radius, this contradiction is particularly pronounced in thinner rods. This parameter manipulation to accommodate the solver inevitably introduces visual defects such as noticeable sag or artificial stretching in real-world physical scenarios.

[0007] On the other hand, existing research has proposed enhancing PD with ADMM to enable it to handle general nonlinear elastic materials. Local updates in this type of method take two forms: projection-based and optimization-based. Directly applying the projection-based form to Cosserat elastic bars severs the information propagation between position and orientation, leading to slower convergence. While the optimization-based form addresses this coupling and accelerates information propagation through local optimization within the near-end operator, it significantly increases the cost of local updates, resulting in limited overall improvement over PD. Furthermore, in interactive applications, non-converged residuals introduce visually perceptible numerical damping. The auxiliary variables in these methods are all defined at the state (position, quaternion) level, which limits the information exchange between position and orientation updates when applied to Cosserat elastic bars. Therefore, the industry urgently needs a Cosserat elastic bar simulation method that is both efficient and robust under real physical parameters. Summary of the Invention

[0008] The purpose of this invention is to address the problem that existing methods cannot efficiently and robustly simulate Cosserat elastic rods under real physical parameters (high stiffness ratio). A strain-based Active Mode Variable Model (ADMM) simulation method for Cosserat elastic rods is proposed. The core of this method is to define the auxiliary variables of the ADMM at the strain level rather than the state level. This directly exposes the position-orientation coupling in the augmented Lagrangian penalty term, allowing the original variable update (x-update) of the ADMM to recover the local-global structure of the dynamic field (PD), while avoiding nonlinear local solutions.

[0009] The technical solution of the present invention is as follows:

[0010] An efficient simulation method for Cosserat elastic rods based on strain-dependent ADMM is proposed. This method establishes a discrete model of the Cosserat elastic rod and arranges position variables at the vertices of the rod. Arrange unit quaternions representing the material frame on the line segments between adjacent vertices. Define stretch-shear (SS) constraints for each line segment and bend-torsion (BT) constraints for each vertex, and iteratively update the position and orientation of the elastic rod at each time step using implicit Euler integrals; the iterative update at each time step includes:

[0011] (1) For tension-shear constraints only, introduce auxiliary variables at the strain level. ,in For the direction of the line segment, For the third axis of the material frame, Let the resting length of the line segment be denoted by , and introduce a scaling dual variable. The position-orientation coupling of the elastic rod is explicitly placed inside the penalty term of the augmented Lagrangian function, and the augmented Lagrangian function is rewritten.

[0012] In some embodiments of the present invention, the augmented Lagrangian function places the nonlinear position-orientation coupling term entirely within the ADMM penalty term, such that the penalty term simultaneously includes both position and orientation variables and explicitly couples them.

[0013] In some embodiments of the present invention, the auxiliary variable Dual variable with scaling Each has closed-loop update and The updated value is used for subsequent elastic rod orientation and position updates.

[0014] (2) Inaccurate solution and local-global recovery of original variable update: Inaccurate solution is performed by two-level mathematical relaxation of block coordinate descent and projection dynamics, decoupling position and orientation into two subproblems and reducing them to local-global iteration;

[0015] In some embodiments of the present invention, the block coordinate descent and the two-level mathematical relaxation of the projection dynamics are both in the original variables of the ADMM. The update was completed within the specified timeframe; the block coordinate descent decouples position and orientation into two subproblems, which are then solved separately using the local-global form of projection dynamics: the position subproblem recovers the global step of projection dynamics and sets the current material frame orientation. By explicitly embedding the right-hand side of the linear system, the target unit quaternion is obtained by rotating the subproblem towards the shortest arc in the local step, and the nonlinear direction penalty is replaced by the quaternion quadratic potential energy in the global step. As a result, the two global system matrices are constant and state-independent, and are pre-decomposed only once during initialization and reused in all time steps. The local steps of each line segment and each vertex are decoupled from each other and can be computed in parallel.

[0016] (3) Reconstruct the bending-torsion constraint to obtain the explicit linear constraint form with respect to quaternions. Bake the rotational coupling between adjacent material frames into the orientation block of the global system matrix. The bending-torsion constraint is processed by augmenting the system matrix.

[0017] In some embodiments of the present invention, the reconstruction of the bending-torsion constraint involves multiplying both ends of the constraint by a unit quaternion. , get about , Linear constraints; assumption of symbolic terms when initializing the global system matrix. To construct constant linear operators and state-related symbol items Correcting the migration to a local step, through the projection vector Dynamically correct the right-hand side terms of the global linear system.

[0018] (4) Adopt a penalty parameter clamping strategy By constraining the tensile-shear penalty parameters to a softer local average bending-torsional stiffness, the benign state of the global linear system is ensured.

[0019] In some embodiments of the present invention, in the penalty parameter clamping strategy, the tension-shear constraint is defined for each line segment, and the bending-torsion constraint is defined for each vertex. The average local bending-torsion stiffness is first calculated for each line segment. Take again This allows the global step-linear system to maintain a good state and fast, stable convergence without manual parameter adjustment at extreme stiffness ratios.

[0020] Advantages and beneficial effects of the present invention:

[0021] (1) It can perform robust simulations without scene-by-scene parameter tuning under real physical parameters (high stiffness ratio), which helps to reduce visual artifacts such as sagging, overstretching, jelly-like oscillation and branch twisting caused by artificially softening stiffness; (2) It retains the advantages of PD one-time pre-decomposition and local-global fast iteration, avoiding nonlinear local solutions; (3) Compared with the existing state hierarchy ADMM, under the representative high stiffness ratio scene (dense tree), the total running time for simulating 100 frames with the same number of iterations is reduced by about 21%, of which the time spent in the local solution stage is about 1 / 8 of that; (4) It has smaller numerical damping under a limited iteration budget, which helps to retain a more vivid high-frequency dynamic response. Attached Figure Description

[0022] Figure 1 This is a schematic diagram illustrating the simulation effect of the present invention applied to dense trees during the application and removal of wind loads;

[0023] Figure 2 This is a schematic diagram of the effect of the present invention applied to hair simulation. From left to right, the diagram shows the initial configuration, the configuration after gravity descent, and several subsequent frames showing the hair roots rotating around the central axis and causing the hair to move.

[0024] Figure 3 This is a schematic diagram comparing the simulation results of the present invention and the traditional PD method on sparse trees under wind load under real physical parameters.

[0025] Figure 4 This diagram illustrates the convergence comparison between the present invention and methods such as Newton's reference solution and PD in the case of a cantilever beam with a high stiffness ratio. Detailed Implementation

[0026] The Cosserat elastic rod simulation method based on strain ADMM provided by this invention will be further described in detail below with reference to the accompanying drawings and specific embodiments. This invention can be implemented on a general-purpose computer using C++ combined with a sparse linear algebra library (such as Eigen), and local steps can be parallelized using OpenMP through multi-threading.

[0027] The method of this invention establishes a discrete model of the Cosserat elastic rod and arranges position variables at the vertices of the elastic rod. Arrange unit quaternions representing the material frame on the line segments between adjacent vertices. (satisfy The four-dimensional vector is used to obtain the material frame of the point by rotating the basis vector of the world coordinate system (used to represent the local rotation orientation). The tension-shear (SS) constraint of each line segment and the bending-torsion (BT) constraint of each vertex are defined, and the position and orientation of the elastic rod are iteratively updated in each time step by implicit Euler integral.

[0028] Based on the above technical concept, the iterative update within each time step of this invention includes:

[0029] I. Augmented Lagrangian Structure. Auxiliary variables are introduced only for the tension-shear (SS) constraint. ,in ( ) represents the direction of the line segment. For the third axis of the material frame, Let the resting length of the line segment be denoted by , and introduce a scaling dual variable. The system is rewritten as the following augmented Lagrange function. :

[0030]

[0031] in ADMM penalty parameters are used to control the coupling strength during global updates. It is the sum of the bending and torsional energies at each vertex.

[0032] Unlike the state-level ADMM, this invention places the nonlinear coupling term within the penalty term, making the position and orientation variables explicitly coupled in the penalty term, thereby directly exchanging information during global updates.

[0033] II. Inaccurate solution and local-global recovery of the original variable update (i.e., x-update of ADMM, corresponding to primal update): Inaccurate solution is performed through two levels of mathematical relaxation of block coordinate descent and projection dynamics, decoupling position and orientation into two subproblems, and equalizing them into local-global iteration.

[0034] To preserve the constant and predecomposable system matrix, this invention does not perform exact joint minimization in the original variable update. Instead, it introduces two levels of mathematical relaxation: block coordinate descent (BCD) relaxation for joint minimization to prevent coupling terms from contaminating the constant linear system; and (ii) introducing a second level of imprecision in the orientation subproblem using standard PD (first local manifold projection, then optimization of quadratic potential in quaternion space). Thus, the original variable update is constructed as an imprecise proxy solution that restores the local-global architecture of PD. Unlike state-level ADMM, which completes PD derivation during the Lagrange assembly stage, this invention postpones PD derivation to the original variable update stage. This is the key difference between this invention and existing ADMM methods.

[0035] Location update: Frozen orientation Auxiliary variables With dual variables The energy of the position subproblem is

[0036]

[0037] The linear system obtained by taking its gradient and setting it to zero precisely recovers the standard global step of PD and sets the current material frame orientation. The right-hand side of the explicitly embedded positional linear system. This positional subproblem corresponds to the PD standard surrogate form described above. Among them, the local degrees of freedom acting on the position The constant finite difference operator is

[0038]

[0039] in It is a 3×3 identity matrix. The local projected target is... .

[0040] Orientation Update: Freeze position, collect direction tracking target The optimal target unit quaternion is obtained through local step rotation using the shortest arc. (That is, rotate and align the current material frame's third axis to the normalized value) ),in Rotate to The shortest arc rotation quaternion is , The global step replaces the nonlinear direction penalty with the quadratic potential energy in the quaternion space.

[0041]

[0042] Its contribution to the global Hessian is a constant diagonal block, thus preserving the system matrix. Similarly, the orientation subproblem corresponds to the standard PD surrogate form, whose constant linear operator is... .

[0043] Auxiliary variables Dual variable with scaling Update: Fixed position and orientation, with auxiliary variables Minimization has a closed-form solution, followed by scaling the dual variable. Update according to standard duality:

[0044]

[0045] Updated auxiliary variables Dual variable with scaling The value is immediately used for subsequent elastic rod orientation and position updates.

[0046] Third, the bending-torsion constraints are reconstructed to obtain an explicit linear constraint form with respect to quaternions, and the rotational coupling between adjacent material frames is baked into the orientation block of the global system matrix. The bending-torsion constraints are handled by system matrix augmentation; the augmented Lagrangian function places the nonlinear position-orientation coupling term entirely within the ADMM penalty term, so that the penalty term simultaneously includes position and orientation variables and explicitly couples them.

[0047] In this invention, both ends of the BT constraint are simultaneously multiplied by a unit quaternion on the left. , get about , Explicit linear constraints During initialization, it is assumed that... Construct a constant linear operator that acts on adjacent oriented blocks.

[0048]

[0049] in for of Right multiply by the matrix. The contribution of this constraint to the Hessian of the system is... Its off-diagonal blocks capture rotational coupling between adjacent material frames; state-dependent sign correction. The vector is transferred to a local step and then projected. Dynamically correct the right-hand side of the global system. In the above expression, the superscript T denotes matrix transpose. It is a 4×4 identity matrix.

[0050] IV. Employing a penalty parameter clamping strategy By constraining the tensile-shear penalty parameters to a softer local average bending-torsional stiffness, the benign state of the global linear system is ensured.

[0051] In this invention, due to the thin rod Compare Larger by several orders of magnitude, if we directly let This can cause the entire system to become highly ill and stagnate at extreme stiffness ratios. This invention first calculates the average local bending-torsional stiffness for each segment. (Take the average bending-torsional stiffness at both ends of the segment, i.e.) ), then take By constraining the penalty to a relatively soft BT stiffness, the well-state of the globally linear system is ensured, allowing for rapid and stable convergence without manual parameter tuning.

[0052] In this embodiment of the invention, the core of the Cosserat elastic rod discrete model and time integration construction is:

[0053] The Cosserat elastic rod is modeled as having an arc length The parameterized 3D curve has its centerline located at... Unit quaternion Encode the material frame at each point. Use a standard staggered layout: Location Located at the vertex, unit quaternion lie in and On the line segment between them. The stretch-shear constraint for each segment is... The bending-torsion constraint at each vertex is: The total elastic energy is the sum of the SS energy of each segment and the BT energy of each vertex, and the stiffness... , ,in , By Young's modulus shear modulus Cross-sectional area With polar moment of inertia Physical quantities are derived. Using implicit Euler integrals, the endpoint states are obtained from the following equation:

[0054]

[0055] in For time steps, For inertial targets, To concentrate mass, a quasi-static assumption is adopted for the material frame.

[0056] The meanings of the symbols in the above formula are as follows: For the stretch-shear constraint of line segment i, Let i be the resting length of line segment i. The direction vector of the third axis of the material frame; For the bending-torsion constraint at vertex i, The sign correction term, used to guide constraints to the nearest quaternion polarity, is determined by... Sure, The relative rotational quaternion in the resting state; total elastic energy The tensile-shear energy of each segment Bending-torsional energy at each vertex The sum constitutes; , The tensile-shear stiffness and bending-torsional stiffness of line segment i are given by the resting length. and material stiffness coefficient , Confirmed, among which For Young's modulus, Shear modulus For cross-sectional area, It is the polar moment of inertia; For speed, This refers to external acceleration.

[0057] In an embodiment of the present invention, the single-step iterative process is as follows:

[0058] In the following iterative process, the superscript t represents the state at the start of the current time step, i.e. , These represent the position and orientation at the start of the current time step.

[0059] Execute the following iterations within each time step to initialize , , , ;Execute the following sequentially for each ADMM iteration:

[0060] Step 1: Update auxiliary variables ;

[0061] Step 2, update the dual variable ;

[0062] Step 3: Assemble and solve for the incremental values ​​of the system. and normalization ;

[0063] Step 4: Assemble and solve the position system to obtain the increment. ,renew ;

[0064] Update speed after iteration and order , To reduce the floating-point truncation error in absolute coordinate evaluation under high stiffness, numerical implementation is unified for incremental... , Solve for the system matrices at the two left ends of the position and orientation under BCD. , All values ​​are constants. During initialization, a sparse LDLT decomposition is performed once and reused across all time steps; each orientation is solved and normalized. It is immediately used to evaluate the position projection in order to propagate the kinematic coupling during iteration.

[0065] Among them, the location system matrix It is a weighted quadratic sum of the inertial mass term and the contributions of each tension-shear constraint; oriented system matrix It is formed by accumulating the constant diagonal blocks contributed by each tension-shear constraint and the rotational coupling blocks contributed by each bending-torsion constraint; the above accumulation refers to accumulating the weighted matrix blocks corresponding to each constraint to the corresponding positions of the global sparse matrix according to the mapping from their local degrees of freedom to global degrees of freedom.

[0066] In the embodiments of the present invention, the key implementation details of the penalty parameter and the unit quaternion constraint are:

[0067] (a) Penalty parameter: take ;when hour Restore PD without loss of speed, when The clamping effect is timely, preventing pathological changes. Ablation experiments show that in a cantilever beam scenario with a stiffness ratio of approximately 596, as... The relative energy error increases monotonically (from) The time of the appointment Rise to The time of the appointment ), among the values ​​tested, will Take as The convergence error is minimized at this time.

[0068] (b) Unit Quaternion Constraint: This continues the practice of not explicitly applying unit norm constraints in the global step, but applying them after each orientation update. The strategy of renormalization: Since the information propagation of BCD and local updates is relatively slow, successive normalization is sufficient to avoid artifacts such as single-segment spurious stretching, without the need to introduce additional penalty terms or Lie algebra parameterization on the rotating manifold.

[0069] The effects of the present invention will be further illustrated below through specific examples.

[0070] Example 1 (Cantilever Beam, see reference) Figure 4):

[0071] A cantilever beam clamped at one end deflects under its own weight, and its stiffness is directly derived from the cross-sectional radius, exhibiting a high stiffness ratio. Using the results of two Newton iterations as a reference solution, the present invention only requires about 10 iterations per time step to reach the equilibrium configuration, while PD (regardless of whether or not the augmented BT constraint of the present invention is used) still fails to converge to the equilibrium configuration even with 200 iterations per step; in this scenario, the present invention takes about 0.12ms per frame, while the Newton reference takes about 0.93ms per frame.

[0072] Example 2 (Sparse tree / Dense tree, see reference) Figure 1 , Figure 3 ):

[0073] The radius of a tree from its trunk to its terminal branch varies greatly (from approximately 0.1m to 10m). -3 (m), the stiffness ratio has a wide spatial distribution range. With a fixed 10 iterations per step, traditional PD needs to... Scaling down the entire tree to force convergence can over-soften the trunk due to the same scaling factor, causing it to tilt excessively in the wind and exhibit jelly-like elastic oscillations after the wind subsides. Insufficient scaling results in visible bending and twisting artifacts on the branches. This invention, under realistic physical parameters and the same iteration budget, can robustly and efficiently calculate the true elastic deformation and rebound of trees under wind, without jelly-like oscillations or twisting. In a dense tree scene, simulating 100 frames with 20 iterations per step, the total time consumption of this invention is reduced by approximately 21% compared to the optimization-based state hierarchy ADMM, and the time consumption in the local solution stage is approximately 1 / 8 of that of ADMM.

[0074] Example 3 (hair, reference) Figure 2 ):

[0075] The cross-sectional radius of a hair is extremely small (approximately With a high stiffness ratio (m), the solver of this invention remains stable and can represent the natural movement of hair under gravity and the rotation of the hair roots around the central axis.

[0076] It should be further noted that the above embodiments and examples are merely for understanding the technical solution of the present invention and are not intended to limit the scope of protection of the present invention. The present invention recovers the fast, predecomposable local-global architecture of PD through strain-level auxiliary variables and two-level relaxation, while faithfully handling the position-orientation coupling that causes classical PD to stagnate at high stiffness ratios. Combined with an augmented BT system matrix and a penalty clamping strategy, it can achieve stable convergence in realistic physical parameter scenarios such as cantilever beams, trees, and hair, and maintain a relatively vivid dynamic response under a finite iteration budget. The present invention can also be combined with Lie algebra-based PD solvers, and is applicable to scenarios involving contact and extensions such as local-global multigrids.

Claims

1. An efficient simulation method for Cosserat elastic rods based on strain-dependent ADMM, wherein the method establishes a discrete model of the Cosserat elastic rod and arranges position variables at the vertices of the elastic rod. Arrange unit quaternions representing the material frame on the line segments between adjacent vertices. Define the tension-shear constraint for each line segment and the bending-torsion constraint for each vertex, and iteratively update the position and orientation of the elastic rod at each time step using implicit Euler integrals; characterized in that... The iterative updates within each time step include: (1) For tension-shear constraints only, introduce auxiliary variables at the strain level. ,in For the direction of the line segment, For the third axis of the material frame, Let the resting length of the line segment be denoted by , and introduce a scaling dual variable. The position-orientation coupling of the elastic rod is explicitly placed inside the penalty term of the augmented Lagrangian function, and the augmented Lagrangian function is rewritten. (2) Inaccurate solution and local-global recovery of original variable update: Inaccurate solution is performed by two-level mathematical relaxation of block coordinate descent and projection dynamics, decoupling position and orientation into two subproblems and reducing them to local-global iteration; (3) Reconstruct the bending-torsion constraint to obtain the explicit linear constraint form with respect to quaternions. Bake the rotational coupling between adjacent material frames into the orientation block of the global system matrix. The bending-torsion constraint is processed by augmenting the system matrix. (4) Adopt a penalty parameter clamping strategy By constraining the tensile-shear penalty parameters to a softer local average bending-torsional stiffness, the benign state of the global linear system is ensured.

2. The method according to claim 1, characterized in that: The augmented Lagrangian function in step (1) places the nonlinear position-orientation coupling term as a whole within the ADMM penalty term, so that the penalty term simultaneously includes both position and orientation variables and explicitly couples them.

3. The method according to claim 1, characterized in that: The auxiliary variable mentioned in step (1) Dual variable with scaling Each has closed-loop update and The updated value is used for subsequent elastic rod orientation and position updates.

4. The method according to claim 1, characterized in that: The two-stage mathematical relaxation of block coordinate descent and projection dynamics in step (2) are both in the original variables of the ADMM. The update was completed within the specified timeframe; the block coordinate descent decouples position and orientation into two subproblems, which are then solved separately using the local-global form of projection dynamics: the position subproblem recovers the global step of projection dynamics and sets the current material frame orientation. By explicitly embedding the right-hand side of the linear system, the target unit quaternion is obtained by rotating the subproblem towards the shortest arc in the local step, and the nonlinear direction penalty is replaced by the quaternion quadratic potential energy in the global step. As a result, the two global system matrices are constant and state-independent, and are pre-decomposed only once during initialization and reused in all time steps. The local steps of each line segment and each vertex are decoupled from each other and can be computed in parallel.

5. The method according to claim 1, characterized in that: The reconstruction of the bending-torsion constraint in step (3) involves multiplying both ends of the constraint by a unit quaternion. , get about , Linear constraints; assumption of symbolic terms when initializing the global system matrix. To construct constant linear operators and state-related symbol items Correcting the migration to a local step, through the projection vector Dynamically correct the right-hand side terms of the global linear system.

6. The method according to claim 1, characterized in that: In the penalty parameter clamping strategy described in step (4), the tension-shear constraint is defined for each line segment, and the bending-torsion constraint is defined for each vertex. First, the average local bending-torsion stiffness is calculated for each line segment. Take again This allows the global step-linear system to maintain a good state and fast, stable convergence without manual parameter adjustment at extreme stiffness ratios.

7. The application of the method according to any one of claims 1 to 6, characterized in that: The method is used for physical simulation and animation generation of any slender elastic structure such as hair, rope, suture, cable, and tree branches under real physical parameters, and can be applied to film and television special effects, virtual reality / augmented reality and interactive simulation systems.

Images