pith. sign in

arxiv: 2604.07307 · v1 · submitted 2026-04-08 · 💻 cs.CE · cond-mat.mtrl-sci· cs.NA· math.NA

Improved Implementation of Approximate Full Mass Matrix Inverse Methods into Material Point Method Simulations

John A. Nairn This is my paper

Pith reviewed 2026-05-10 17:18 UTC · model grok-4.3

classification 💻 cs.CE cond-mat.mtrl-scics.NAmath.NA
keywords material point methodFMPMfull mass matrixgrid velocitiesMPM implementationcontact algorithmsboundary conditionsnumerical stability
0
0 comments X

The pith

A revised FMPM loop lets approximate full mass matrix methods work with standard MPM features like contact and boundary conditions.

A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.

The paper derives a simplified FMPM(k) implementation that clarifies the loop for inserting approximate full mass matrix inversion into any MPM code. It then modifies that loop so FMPM(k) can coexist with other common MPM features that depend on lumped-mass calculations, such as velocity boundary conditions, multimaterial contact, crack contact, and imperfect interfaces. The revised approach performs the needed calculations only once per time step, just before updating particle positions and velocities. A sympathetic reader would care because these changes make the accuracy gains from full mass matrix methods available in the kinds of complex simulations where MPM is typically used.

Core claim

The paper first derives a revised FMPM(k) implementation that both simplifies and clarifies the FMPM Loop that can be added to MPM codes. Next, that loop is modified to allow FMPM(k) to work well even in simulations that need other MPM features that previously caused conflicts. The implementation requires these calculations only once per time step just before updating particle positions and velocities. Additional discussion covers apparent loss of stability at very high order k and inherent computational cost.

What carries the argument

The revised FMPM Loop, a calculation sequence inserted into MPM codes to compute grid velocities via approximate full mass matrix inversion of order k while handling lumped-mass features.

If this is right

  • FMPM(k) calculations need to occur only once per time step before particle updates.
  • Simulations using velocity boundary conditions, multimaterial contact, crack contact, or imperfect interfaces can now use the accuracy benefits of FMPM(k).
  • Stability of the method varies with order k and may degrade at very high orders.
  • Options exist to mitigate the computational cost while retaining the revised loop.

Where Pith is reading between the lines

These are editorial extensions of the paper, not claims the author makes directly.

  • The simplification of the loop could make it easier to add FMPM(k) to existing open-source or custom MPM codes without extensive refactoring.
  • Similar loop modifications might address mass-matrix conflicts in other particle or meshfree methods that mix inversion techniques with contact algorithms.
  • Systematic testing across increasing values of k in production-scale runs would help identify the practical optimum order for accuracy versus stability.

Load-bearing premise

The modified FMPM loop preserves the accuracy gains of full mass matrix inversion while fully resolving conflicts with lumped-mass features such as contact and boundary conditions, without introducing new instabilities or errors.

What would settle it

A simulation run with multimaterial contact and grid-based velocity boundary conditions that compares particle velocities, stability, and accuracy between the original FMPM(k) and the modified loop version.

Figures

Figures reproduced from arXiv: 2604.07307 by John A. Nairn.

Figure 1
Figure 1. Figure 1: Flow chart for an MPM time step. The branches at the two “?” decision points refer to calculations using USF, USL, USL+, USAVG, and USAVG+ methods for updating particle stresses and strains. All FMPM(k) calculations are in tasks 1.c, 3.b, and 4.c. When using the recommended USL method, the only FMPM(k) task in each time step is 3.b. that updates particle position and velocity with the same grid velocities … view at source ↗
Figure 2
Figure 2. Figure 2: A moving wall on the right edge of a particle domain should set velocity boundary conditions on all circled nodes. d is the wall depth that determines which nodes inside the particle domain get velocity boundary conditions. solely by moving wall velocity BCs on each end that set spatial velocity and gradient in the x direction to v (b) = εx˙ wall 1 + ˙εt and ∇v (b) = ε˙ 1 + ˙εt =⇒ v (b) i = εx˙ i 1 + ˙εt B… view at source ↗
Figure 3
Figure 3. Figure 3: A. Error of FMPM(k) calculations as function of order k using the new methods for velocity BCs and the old combined method from Ref. [2]. B. Error in calculations as a function of depth of the wall in the moving-wall velocity BCs for FLIP and FMPM(k) for k = 1, 2, 4, and 8. Solid lines used CPDI shape function and dashed lines using B2CPDI shape functions. The solid, red lines are the old combined method f… view at source ↗
Figure 4
Figure 4. Figure 4: Pressure shock wave passing a material interface using “stick” contact conditions by FLIP, FMPM(1), FMPM(2), and FMPM(k ≥ 4); FMPM(k ≥ 2) used either “Evolving” or “Net” methods with identical results. All results were also identical to simulations run in single material mode. The pressure wave is moving towards the left with leading edge a 34 mm and trailing edge near the wall at 46 mm. The dashed lines s… view at source ↗
Figure 5
Figure 5. Figure 5: Shock wave passing a material interface using frictional contact conditions by FMPM(k ≥ 4) using the “Net” method and FMPM(4) or FMPM(8) using the “Evolving” method. The pressure wave is moving towards the left with leading edge a 34 mm and trailing edge near the wall at 46 mm. The dashed lines show the position of the material interface and the moving wall when the simulations was stopped. MPM results wer… view at source ↗
Figure 6
Figure 6. Figure 6: A plot of ∥fc∥/Ac as a function of Nc for three values of Tc and µ = 0.8. Note that the plot for Tc = 0 is also the plot for frictionless contact for all values of Tc. This contact model is plotted in [PITH_FULL_IMAGE:figures/full_fig_p019_6.png] view at source ↗
Figure 7
Figure 7. Figure 7: Snapshots of two disks in off-center impact. The color is plotting temperature (from black for 300 K to red for 456 K). the cell size was 1/3 mm, coefficient of friction was µ = 0.6, and calculations used FMPM(4). Arrows indicate direction of motion [PITH_FULL_IMAGE:figures/full_fig_p020_7.png] view at source ↗
Figure 8
Figure 8. Figure 8: The trajectories of the center of mass for left and right disks from 0 to 16 ms using FLIP (solid lines), FMPM(4) (red dots only for left disk) or FMPM(1) (blue dots only for right disk) for µ = 0, 0.3, and 0.6. Simulations used cell size = 1/3 mm (60 cells along diameter of the disks) and C = 0.2. For clarity, the y positions are relative to their starting positions and x positions for left and right disk… view at source ↗
Figure 9
Figure 9. Figure 9: Internal plus kinetic energy (U + K) energy loss at the end of the impact event (recorded at 16 ms) using cell size 1/3 mm and C = 0.2. A. FLIP (solid lines) and FMPM(4) (dashed lines) results for µ = 0, 0.3, or 0.6. B. FLIP results (horizontal line) compared to FMPM(k) as a function of k using µ = 0.3. 21 [PITH_FULL_IMAGE:figures/full_fig_p021_9.png] view at source ↗
Figure 10
Figure 10. Figure 10: Total Helmholz free energy as a function of time for cell size = 1/3 mm and C = 0.2. The peak is caused by the impact event. The post-impact value returns close to zero with oscillations likely caused by remaining stress waves after the disks have fully separated (i.e., not an equilibrium state). separated. The presence of oscillating stress waves means the final state has not yet returned to a new equili… view at source ↗
Figure 11
Figure 11. Figure 11: Internal plus kinetic energy (U + K) energy loss at the end of the impact event (recorded at 16 ms) using µ = 0.3. A. Energy loss as function of cell size for C = 0.2 for FLIP or FMPM(4) with update methods USL or USL+. B. Energy loss as a function of Courant number for cell size 1/3 mm by FLIP, FMPM(4), or periodic FMPM(4). simulation has more time steps. If the FMPM loop causes dissipation, even if only… view at source ↗
Figure 12
Figure 12. Figure 12: A. Stability limit for FMPM(k), FMPM(k) blended with FMPM(1) using α = 0.8, FMPM(k) blended with FMPM(2) using α = 0.8, and periodic FMPM(k) using CX = 2 compared to stability limit for FLIP (dashed line at 0.84). B. Energy dissipation at the end of five vibrations using FMPM(4) as a function of fraction FMPM(4) (or α) for fixed C = 0.539. The only change needed to the FMPM loop in [PITH_FULL_IMAGE:figur… view at source ↗
Figure 13
Figure 13. Figure 13: A. Sample convergence of “Means” and “Changes” metrics for time step in the middle of the block impact or shock wave simulations. B Energy loss and calculation time relative to FLIP calculations as a function of critical metric to find FMPM(k) convergences for block impact simulations. Solid lines used “Means” metric. Dashed lines used “Changes” metric. Horizontal lines plot results for non-dynamic FMPM(2… view at source ↗
Figure 14
Figure 14. Figure 14: The average FMPM(k) order for convergence as a function of the convergence criterion for the two-block impact and shock wave examples. Solid lines are when using the “Means” metric while dashed lines are when using the “Changes” metric. another test of dyamic FMPM(k), the method was tried on the shock wave example. Figure 13A shows convergence of the two metrics for an arbitrary time step in the middle of… view at source ↗
read the original abstract

Approximate full mass matrix methods for the material point method, known as FMPM(k) of order k, can improve the calculation of grid velocities from grid momentum. It can be implemented in any MPM code by inserting a new calculation task whenever grid velocities are needed. The implementation recommended in this paper only needs these calculations once per time step just before when updating particle positions and velocities. FMPM implementation issues arise, however, when its methods are mixed with other MPM feature that rely on lumped mass calculations. Some common lumped-mass MPM features are grid-based, velocity boundary condition, multimaterial contact calculations, crack contact calculations, and imperfect interfaces. This paper first derives a revised FMPM(k) implementation that both simplifies and clarifies the "FMPM Loop" that can be added to MPM codes. Next, that loop is modified to allow FMPM(k) to work well even in simulations that need other MPM features that previously caused conflicts. Two other FMPM(k) issues are apparent loss of stability at very higher order k and inherent computational cost. These issues are discussed in an analysis of temporal stability as a function of order k and in consideration of options to improve efficiency.

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, simulated authors' rebuttal, and a circularity audit. Tearing a paper down is the easy half of reading it; the pith above is the substance, this is the friction.

Referee Report

2 major / 2 minor

Summary. The paper claims to derive a revised FMPM(k) implementation for the Material Point Method that simplifies the 'FMPM Loop' for insertion into existing MPM codes, requiring the full-mass calculations only once per timestep just before particle position/velocity updates. It then modifies this loop to resolve prior conflicts with lumped-mass-dependent features including grid-based velocity boundary conditions, multimaterial contact, crack contact, and imperfect interfaces. The work additionally analyzes temporal stability as a function of order k and discusses options for improving computational efficiency of the approximate full-mass inversion.

Significance. If the modified loop preserves the documented accuracy gains of FMPM(k) while enabling seamless use with standard MPM contact and boundary features, the contribution would be significant for practical adoption in computational mechanics codes. The derivation of a clarified, once-per-step loop and the stability analysis as a function of k provide concrete implementation guidance and insight into high-order behavior that could reduce barriers to using full-mass methods in complex engineering simulations.

major comments (2)
  1. [Modified FMPM loop] Modified FMPM loop (as described following the initial derivation): The central claim that the loop can be isolated or reordered to coexist with velocity BCs, multimaterial contact, crack contact, and imperfect interfaces without corrupting lumped-mass quantities or introducing new truncation errors, altered contact normals, or stability shifts is load-bearing but unsupported by explicit verification. The stability analysis is stated only for plain FMPM(k); no corresponding checks or error metrics are provided for the combined setting with contact/BC cases.
  2. [Abstract and implementation sections] Abstract and implementation sections: The paper asserts that the modifications 'allow FMPM(k) to work well' with previously conflicting features, yet provides no quantitative demonstration (e.g., convergence rates, contact force errors, or stability thresholds) comparing the modified loop against both standard lumped-mass MPM and unmodified FMPM(k) on benchmark problems that exercise those features.
minor comments (2)
  1. [Abstract] Abstract: The summary outlines the contributions but contains no equations, pseudocode for the revised loop, or numerical results, which is acceptable for an abstract but reduces immediate clarity on the exact form of the simplification and modification.
  2. [Efficiency discussion] Efficiency discussion: Options for improving computational cost are mentioned but lack concrete benchmarks, flop counts, or scaling comparisons relative to the original FMPM(k) implementation or standard MPM.

Simulated Author's Rebuttal

2 responses · 0 unresolved

We thank the referee for the careful review and constructive comments on our manuscript. The feedback highlights the need for stronger verification of the proposed modifications. We address each major comment below and commit to revisions that will incorporate the requested demonstrations and checks.

read point-by-point responses

  1. Referee: [Modified FMPM loop] Modified FMPM loop (as described following the initial derivation): The central claim that the loop can be isolated or reordered to coexist with velocity BCs, multimaterial contact, crack contact, and imperfect interfaces without corrupting lumped-mass quantities or introducing new truncation errors, altered contact normals, or stability shifts is load-bearing but unsupported by explicit verification. The stability analysis is stated only for plain FMPM(k); no corresponding checks or error metrics are provided for the combined setting with contact/BC cases.

    Authors: We agree that the manuscript's derivation of the reordered loop is intended to ensure compatibility without corrupting lumped-mass quantities or introducing the listed errors, but explicit numerical verification for the combined settings is not provided. The stability analysis focuses on the base FMPM(k) method. In the revised manuscript, we will add benchmark simulations exercising velocity boundary conditions, multimaterial contact, crack contact, and imperfect interfaces. These will include error metrics, contact force comparisons, and stability checks to confirm the modified loop preserves accuracy and does not shift stability thresholds relative to the plain case. revision: yes

  2. Referee: [Abstract and implementation sections] Abstract and implementation sections: The paper asserts that the modifications 'allow FMPM(k) to work well' with previously conflicting features, yet provides no quantitative demonstration (e.g., convergence rates, contact force errors, or stability thresholds) comparing the modified loop against both standard lumped-mass MPM and unmodified FMPM(k) on benchmark problems that exercise those features.

    Authors: The abstract and implementation sections summarize the derivation and modifications to resolve prior conflicts. We acknowledge that no quantitative comparisons (such as convergence rates or contact force errors) against standard MPM and unmodified FMPM(k) are included for the relevant benchmarks. We will revise these sections and add a dedicated results subsection with the requested quantitative demonstrations on benchmark problems that exercise the lumped-mass features. revision: yes

Circularity Check

0 steps flagged

No circularity: algorithmic derivation of revised FMPM loop is self-contained

full rationale

The paper's central contribution is an explicit derivation of a simplified FMPM(k) loop (inserted once per timestep before particle updates) followed by targeted modifications to isolate it from lumped-mass features such as velocity BCs, multimaterial contact, and crack contact. These steps are presented as implementation changes with accompanying stability analysis as a function of k; no equations reduce by construction to fitted parameters, self-citations, or renamed inputs. The abstract and description contain no load-bearing self-referential definitions or uniqueness theorems imported from prior author work. The derivation chain is therefore independent of its own outputs.

Axiom & Free-Parameter Ledger

0 free parameters · 0 axioms · 0 invented entities

Only the abstract is available, so no specific free parameters, axioms, or invented entities can be identified from the text. The work appears to rest on standard MPM framework assumptions without introducing new entities.

pith-pipeline@v0.9.0 · 5523 in / 1131 out tokens · 40153 ms · 2026-05-10T17:18:00.486534+00:00 · methodology

discussion (0)

Sign in with ORCID, Apple, or X to comment. Anyone can read and Pith papers without signing in.

Reference graph

Works this paper leans on

32 extracted references · 32 canonical work pages

  1. [1]

    Hammerquist and John A

    Chad C. Hammerquist and John A. Nairn. A new method for material point method particle updates that reduces noise and enhances stability.Computer Methods in Applied Mechanics and Engineering, 318:724–738, 2017. 30

    work page 2017

  2. [2]

    Nairn and Chad C

    John A. Nairn and Chad C. Hammerquist. Material point method simulations using an ap- proximate full mass matrix inverse.Computer Methods in Applied Mechanics and Engineering, 337:113667, 2021

    work page 2021

  3. [3]

    J. A. Nairn. Coupling transport equations to mechanics in the material point method using an approximate full capacity matrix inverse.Computer Methods in Applied Mechanics and Engineering, 420:116757, 2024

    work page 2024

  4. [4]

    S. G. Bardenhagen, J. E. Guilkey, K. M. Roessig, J. U. Brackbill, W. M. Witzel, and J. C. Foster. An improved contact algorithm for the material point method and application to stress propagation in granular material.Computer Modeling in Engineering & Sciences, 2:509–522, 2001

    work page 2001

  5. [5]

    Nairn, S

    John A. Nairn, S. G. Bardenhagen, and G. S. Smith. Generalized contact and improved frictional heating in the material point method.Computational Particle Mechanics, 5(3):285– 296, 2018

    work page 2018

  6. [6]

    Nairn, Chad C

    John A. Nairn, Chad C. Hammerquist, and Grant Smith. New material point method contact algorithms for improved accuracy, large-deformation problems, and proper null-space filtering. Computer Methods in Applied Mechanics and Engineering 362, 362:112859, 2020

    work page 2020

  7. [7]

    Material point method calculations with explicit cracks.Computer Modeling in Engineering and Sciences, 4(6):649–664, 2003

    John A Nairn. Material point method calculations with explicit cracks.Computer Modeling in Engineering and Sciences, 4(6):649–664, 2003

    work page 2003

  8. [8]

    J. A. Nairn and Y. E. Aimene. Modeling multiple crack propagation in the material point method by j-integral methods accounting for other cracks intersecting the j contour.Engineer- ing Fracture Mechanics, 322:111143, 2025

    work page 2025

  9. [9]

    John A. Nairn. Numerical implementation of imperfect interfaces.Computational Materials Science, 40:525–536, 2007

    work page 2007

  10. [10]

    John A. Nairn. Modeling of imperfect interfaces in the material point method using multima- terial methods.Computer Modeling in Eng. & Sci., 92(3):271–299, 2013

    work page 2013

  11. [11]

    A hybrid penalty and grid based con- tact method for the material point method.Computer Methods in Applied Mechanics and Engineering, 379:113739, 2021

    James Guilkey, Robert Lander, and Linda Bonnell. A hybrid penalty and grid based con- tact method for the material point method.Computer Methods in Applied Mechanics and Engineering, 379:113739, 2021

    work page 2021

  12. [12]

    John A. Nairn. Material point method (OSparticulas and Nairn- MPM) and finite element analysis (NairnFEA) open-source software. http://osupdocs.forestry.oregonstate.edu/index.php/Main Page, 2026

    work page 2026

  13. [13]

    John A. Nairn. Source code repository for open-source material point method (NairnMPM) and finite element analysis (NairnFEA) software. https://github.com/nairnj/nairn-mpm-fea, March 2026

    work page 2026

  14. [14]

    Sulsky, Z

    D. Sulsky, Z. Chen, and H. L. Schreyer. A particle method for history-dependent materials. Comput. Methods Appl. Mech. Engrg., 118:179–186, 1994

    work page 1994

  15. [15]

    Sulsky, S.-J

    D. Sulsky, S.-J. Zhou, and H. L. Schreyer. Application of a particle-in-cell method to solid mechanics.Comput. Phys. Commun., 87:236–252, 1995. 31

    work page 1995

  16. [16]

    S. G. Bardenhagen and E. M. Kober. The generalized interpolation material point method. Computer Modeling in Engineering & Sciences, 5:477–496, 2004

    work page 2004

  17. [17]

    Sadeghirad, R

    A. Sadeghirad, R. M. Brannon, and J. Burghardt. A convected particle domain interpolation technique to extend applicability of the material point method for problems involving massive deformations.Int. J. Num. Meth. Engng., 86(12):1435–1456, 2011

    work page 2011

  18. [18]

    Nairn and J

    John A. Nairn and J. E. Guilkey. Axisymmetric form of the generalized interpolation material point method.Int. J. for Numerical Methods in Engineering, 101:127–147, 2015

    work page 2015

  19. [19]

    Steffen, P

    M. Steffen, P. C. Wallstedt, J. E. Guilkey, R.M. Kirby, and M. Berzins. Examination and anal- ysis of implementation choices within the material point method (MPM).Computer Modeling in Engineering & Sciences, 31(2):107–127, 2008

    work page 2008

  20. [20]

    A generalized inverse for matrices.Proceedings of the Cambridge Philosophical Society, 51(3):406–413, 1955

    Roger Penrose. A generalized inverse for matrices.Proceedings of the Cambridge Philosophical Society, 51(3):406–413, 1955

    work page 1955

  21. [21]

    J. U. Brackbill and H. M. Ruppel. FLIP: A method for adaptively zoned, particle-in-cell calculations of fluid flows in two dimensions.Journal of Computational Physics, 65(2):314 – 343, 1986

    work page 1986

  22. [22]

    S. G. Bardenhagen. Energy conservation error in the material point method.J. Comp. Phys., 180:383–403, 2002

    work page 2002

  23. [23]

    Zhou.The Numerical Prediction of Material Failure Based on the Material Point Method

    S. Zhou.The Numerical Prediction of Material Failure Based on the Material Point Method. PhD thesis, University of New Mexico, 1998

    work page 1998

  24. [24]

    Kamojjala, R

    K. Kamojjala, R. Brannon, A. Sadeghirad, and J. Guilkey. Verification tests in solid mechanics. Engineering with Computers, 31(2):193–213, 2015

    work page 2015

  25. [25]

    R. W. Ogden.Non-Linear Elastic Deformations. Ellis-Harwood, New York, 1984

    work page 1984

  26. [26]

    O. C. Zienkiewicz and R. L. Taylor.The Finite Element Methods for Solid and Structural Mechanics. Elsevier Butterworth-Heinemann, Oxford, UK, 2000

    work page 2000

  27. [27]

    Courant, K

    R. Courant, K. Friedrichs, and H. Lewy. On the partial difference equations of mathematical physics.IBM J. Res. Dev., 11(2):215–234, March 1967

    work page 1967

  28. [28]

    Von Neumann and R

    J. Von Neumann and R. D. Richtmyer. A method for the numerical calculation of hydrody- namic shocks.J. Appl. Phys., 21:232–237, 1950

    work page 1950

  29. [29]

    Stability analyses and instability mitigation for the material point method.Computer Methods in Applied Mechanics and Engineering, 452:118784, 2026

    Wen-Chia Yang and Deborah Sulsky. Stability analyses and instability mitigation for the material point method.Computer Methods in Applied Mechanics and Engineering, 452:118784, 2026

    work page 2026

  30. [30]

    P. C. Wallstedt and J. E. Guilkey. Improved velocity projection for the material point method. CMES, 19(3):223–232, 2007

    work page 2007

  31. [31]

    Chris Gritton, Martin Berzins, and Robert M. Kirby. Improving accuracy in particle methods using null spaces and filters. InProceedings of the 4th International Conference on Particle- Based Methods - Fundamentals and Applications, PARTICLES 2015, pages 202–213. Interna- tional Center for Numerical Methods in Engineering, 2015

    work page 2015

  32. [32]

    J. A Nairn. Improved modeling of imperfect interfaces in material point method simulations. in preparation, 2026. 32

    work page 2026