Pith. sign in

REVIEW 3 major objections 5 minor 38 references

Stencil scaling for vector-valued PDEs on hybrid grids with applications to generalized Newtonian fluids

T0 review · 3 major / 5 minor · reviewed 2026-08-14 · deepseek-v4-flash

Pith's one-line read This paper shows that vector-valued PDE solves on hierarchical hybrid grids can match standard nodal-integration accuracy at roughly half the time, by scaling precomputed reference stencils with local coefficients plus an antisymmetric…

desk verdict Solid vector-valued stencil scaling that actually derives the needed correction term, but the benchmarked operator is nonsymmetric at macro boundaries and that perturbation is asserted, not analyzed. read the letter →

arxiv 1908.08666 v2 pith:AM2NOY4J submitted 2019-08-23 cs.CE

classification cs.CE MSC 65N3065N5565Y0565Y20
keywords matrix-freefiniteelementsstencilscalingvariablecoefficientsvector-valuedPDEshierarchicalhybridgridslinearelastostaticsStokesflowgeneralizedNewtonianfluids
verification ladder T0 review T1 audit T2 compute T3 formal

The pith

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

The reading

This paper aims to show that low-order finite element solves of vector-valued PDEs can be made matrix-free and fast at the same time. On hierarchical hybrid grids, every interior node sees the same element-neighborhood pattern, so the constant-coefficient discretization can be stored once as a reference stencil. The paper's proposal is to scale that stencil by locally varying coefficients instead of integrating the bilinear form element by element, and to add an antisymmetric correction term that restores the exact nodal-integration operator when the coefficient is affine on each element patch. If the correction is omitted, the scaling discretizes a different PDE; with it, numerical errors match standard nodal integration while time-to-solution drops to 45–61% for linear elastostatics and 75–87% for the Stokes variants tested, including a nonlinear shear-thinning fluid.

What carries the argument

The load-bearing object is the reference stencil $\hat S^T_{ij}$, the $d\times d$ matrix of the constant-coefficient bilinear form for strain and divergence operators, paired with the symmetric/antisymmetric split of local stiffness tensors: $a^t_{ij}=a^{s;t}_{ij}+a^{a;t}_{ij}$. Mirror-reflected elements $t$ and $t^m$ around a stencil edge satisfy $a^{s;t}_{ij}=a^{s;t^m}_{ij}$ and $a^{a;t}_{ij}=-a^{a;t^m}_{ij}$, because $\nabla\varphi^t_i=-\nabla\varphi^{t^m}_j$. These relations let the coefficient-weighted element sum in nodal integration be rewritten as a patch-average scaling of the reference stencil plus the antisymmetric correction $R^T(k)_{ij}=\tfrac12\sum_t(\bar k_t-\bar k_{t^m})a^{a;t}_{ij}$. When $k$ is constant the correction vanishes and naive scaling is exact; when $k$ varies, the correction encodes the coefficient asymmetry and must be kept.

What would settle it

On a sequence of meshes with many macro elements, assemble the operator exactly as implemented, compute the norm of $A-A^\top$, and run a symmetric Krylov solver against a version with symmetrized interface contributions; if iteration counts or solution errors diverge as macro-boundary degrees of freedom grow, the higher-order perturbation assumption fails.

Watch

Extended reading notes

Core claim

Under the assumption that the coefficient $k$ is affine linear on each local element patch $\omega_{i,j;T}$, Lemma 2.1 decomposes the nodal-integration bilinear form into a scaled reference stencil plus a correction: $\hat a_h(v_h,w_h)= \sum_T \sum_{i,j} \sum_{l,m} (\hat k^T_{ij}(\hat S^T_{ij})_{lm} + (R^T(k))_{ij})(v^{(j)})_l(w^{(i)})_m$, where $\hat k^T_{ij}$ is the patch average of $k$ and $R^T(k)$ is antisymmetric and vanishes for constant $k$. The antisymmetric part is exactly what prevents the naive scaled stencil from discretizing a different equation: for $-\nabla\cdot(k\varepsilon(u))$, dropping $R$ turns the operator into $-\nabla\cdot A(u)$ plus a first-order term. With $R$ included, the physical scaling matches standard element-wise nodal integration in the experiments, in both error and multigrid iteration counts, while using a fraction of the runtime. In 2D the correction is a fixed $2\times2$ antisymmetric matrix scaled by one coefficient difference; in 3D it requires three additional stored stencils per macro element.

Load-bearing premise

The practical implementation applies the scaled stencil only in macro-element interiors and falls back to the standard bilinear form for nodes on macro-element boundaries, and the claim that this loss of global symmetry is a harmless higher-order effect is asserted without proof.

Editorial extensions

If this is right

  • Linear elastostatics with a spatially varying Young's modulus can be solved with essentially the same discrete errors and multigrid iteration counts as nodal integration, while using 45–61% of the time-to-solution across the tested coefficient frequencies.
  • For incompressible Stokes flow the velocity block enjoys the same savings, but total solve time drops only to 75–87% because the coefficient-independent divergence and stabilization blocks cost the same in both approaches.
  • Shear-thinning generalized Newtonian fluids fit the same pipeline: a local least-squares projection of the piecewise-constant strain rate supplies node-wise viscosities, and the fixed-point iterations converge to solutions matching nodal integration at about 87% of its time-to-solution on the finest level.
  • The storage needed for the operator itself is constant per macro element in the matrix-free variants, whereas stored-matrix approaches require memory proportional to the number of degrees of freedom.
  • The method relies on the regular element-neighborhood pattern created by uniform refinement of an unstructured coarse triangulation; it does not transfer directly to general unstructured meshes.

Reading between the lines

Editorial extensions of the paper, not claims the author makes directly.

  • Enforcing exact global symmetry across macro-element boundaries by averaging the two one-sided interface bilinear forms is a natural hardening step; it would remove the unproved symmetry-loss caveat at the price of a little extra interface work, and would let symmetric solvers run without qualification.
  • The same symmetric/antisymmetric splitting should apply to other vector-valued operators whose local tensors flip sign under element reflection, such as anisotropic elasticity or couple-stress terms, provided the coefficient is smooth enough on patch scale.
  • In 3D the three stored correction stencils partially offset the memory savings; compressing them via the color-class symmetries of the refined tetrahedra is an obvious optimization the paper leaves open.
  • Accuracy for non-smooth coefficients is likely controlled by the local least-squares recovery that maps element-wise viscosities to nodes; comparing alternative recovery rules on a discontinuous-viscosity Stokes problem would reveal how far the method extends beyond smooth test cases.
Share X Bluesky LinkedIn Reddit HN

Signed reviews

No signed human review yet.

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, and a circularity audit.

Referee Report

3 major / 5 minor

Summary. The paper proposes a matrix-free stencil scaling technique for vector-valued second-order elliptic PDEs on hierarchical hybrid grids (HHG). The central theoretical result, Lemma 2.1, shows that under the assumption that the coefficient is affine linear on each local element patch, the nodal-integration bilinear form (2.5) decomposes exactly into a scaled constant reference stencil plus an antisymmetric correction term, yielding the physical-scaling form (2.6). Closed-form correction terms are derived for 2D and 3D, and a cost/memory analysis compares physical scaling with on-the-fly nodal integration and stored-stencil approaches. Numerical experiments cover linear elastostatics, a heterogeneous metal-foam elasticity problem, an incompressible Stokes problem with a viscosity jump, and a non-linear Carreau generalized Newtonian fluid. The reported results show that physical scaling achieves errors matching nodal integration while reducing time-to-solution by roughly 35% to 55% in the elasticity benchmarks and by smaller amounts in the Stokes variants.

Significance. If the method holds up, it is significant for large-scale low-order finite element simulations on HHG-type grids: it enables matrix-free operator evaluation with variable coefficients at a fraction of the runtime cost of standard on-the-fly integration, and it substantially reduces memory traffic compared with stored matrices. The paper's strengths include a first-principles derivation of the correction terms, a detailed performance model validated by hardware counters and roofline analysis, and large-scale runs up to about 1e11 degrees of freedom on SuperMUC-NG. The linear elasticity benchmarks provide quantitative error and convergence tables that support the accuracy claim in the coefficient regimes tested. However, the practical implementation deviates from the analyzed symmetric operator at macro-element boundaries, and the nonlinear/discontinuous-coefficient examples lack quantitative error comparisons; these are the main weaknesses.

major comments (3)
  1. [Section 2.2 (paragraph after Remark 2.2)] The operator actually implemented and benchmarked is not the symmetric bilinear form \hat a_h analyzed in Lemma 2.1; the paper states that in practice the standard bilinear form is used whenever the center node lies on a macro-element boundary, which makes the global operator nonsymmetric. The text calls this symmetry loss a "higher order perturbation" but provides no estimate of the norm of the nonsymmetric part, no consistency or convergence analysis, and no solver-specific perturbation bound for the hybrid Gauss-Seidel multigrid, inexact Uzawa, or MINRES iterations. Since all numerical results in Section 5 are obtained with this hybrid operator, the central claim that physical scaling reproduces nodal integration at reduced cost rests on an unverified assumption. The authors should either implement and benchmark the symmetric variant they describe, or supply a rigorous bound on the perturbation and demonstrate that the solver behavior is unaffected across a wider range of coefficient contrasts and mesh sizes.
  2. [Section 5.2.1 (geophysics example)] The viscosity field in this example has a jump across a horizontal plane, so the coefficient is not affine on the patches covered by Lemma 2.1. The paper says the jump is resolved via the macro mesh and that standard on-the-fly integration is performed across these interfaces, which again invokes the hybrid heuristic from Section 2.2. No quantitative error table is reported for this example; only iteration counts, time-to-solution, and a visual comparison of streamlines are given. Because this is a regime outside the lemma's assumptions, a quantitative accuracy comparison against nodal integration is needed to support the claim that physical scaling remains accurate for discontinuous coefficients.
  3. [Section 5.2.2 (nonlinear Carreau example)] For the generalized Newtonian fluid, the viscosity is piecewise constant per element and is only mapped to nodal values via the local least-squares procedure of Section 3, so the assumptions of Lemma 2.1 do not apply. The comparison in Figure 12 is qualitative (curve coincidence), and Table 9 reports only time-to-solution and iteration counts, not an error norm. The central claim that physical scaling reproduces nodal integration for generalized Newtonian fluids would be substantially strengthened by reporting discrete L2 errors of velocity and pressure between the two approaches on each refinement level.
minor comments (5)
  1. [Section 2.2] In the redefinition of \hat a_h, "if xi∈∂T and xj∈∂T of at least one T∈T_H" reads more clearly as "for at least one T∈T_H"; the same applies to the following sentence describing the alternative practical definition.
  2. [Lemma 2.1 proof] The displayed equation rewriting the sum over elements as a half-sum over pairs t and t^m slightly abuses the summation index because both t and t^m lie in the same set; a sentence clarifying the pairing over reflected elements would improve readability.
  3. [Table 5] For m=8 at level 6, the physical scaling error is 6.90e-6 versus 5.54e-6 for nodal integration, about 25% larger; a brief interpretation of this difference (for example, whether it is consistent with the local patch-averaging of the coefficient) would be helpful.
  4. [References] Reference [13] contains a typo in the author name "Eisenst a t", which should be "Eisenstat".
  5. [Section 5.2.1] The sentence "it is possible to resolve the jump via the macro mesh, since the standard on-the-fly integration is performed across these interfaces" is slightly confusing; rephrasing to state explicitly that the macro mesh is aligned with the jump so that the jump lies on macro boundaries where the standard form is used would improve clarity.

Circularity Check

0 steps flagged · score 0.0 of 10

No circularity: physical scaling is an exact algebraic decomposition validated against an external nodal-integration baseline.

full rationale

The central derivation is self-contained. Lemma 2.1 proves, under an affine-coefficient patch assumption, that the nodal-integration bilinear form (2.5) algebraically equals a scaled constant reference stencil plus an antisymmetric correction term, yielding (2.6). No parameter is fitted to the quantities later reported as predictions, and the physical-scaling operator is not defined in terms of the benchmark outcomes. The numerical sections compare the scaled operator against the standard nodal-integration operator as an external baseline, and the Carreau parameters are taken from the experimental literature (reference [16]). Citations to the authors' own HHG framework and prior scalar stencil scaling work [4] provide infrastructure and context; the vector-valued decomposition and its correction term are proven in this paper rather than imported through self-citation. The only flagged limitation is the practical macro-boundary fallback described in Section 2.2 after Remark 2.2: the implemented operator uses the standard form on boundary rows, losing global symmetry, and the paper calls this a higher-order perturbation without a proof. That is a correctness and robustness caveat, not a circularity, because it does not reduce the claimed result to its inputs. No circular step is established.

Assumptions & free parameters 0 free parameters · 3 assumptions · 0 invented entities

No free parameters are fitted in this paper; Carreau model parameters are taken from the experimental literature (Section 5.2.2). The only invented objects are algorithmic, namely reference and correction stencils, not physical entities.

assumptions (3)
  • domain assumption The HHG mesh family provides identical stencil neighborhoods for every interior node of a macro element, with reflected element pairs across each interior edge.
    This structural property is stated in Section 2 and is the basis for the constant reference stencil and for the cancellation of symmetric and antisymmetric parts in Lemma 2.1. It holds for the grids used but is not universal for other mesh families.
  • domain assumption In Lemma 2.1 the coefficient k is affine linear on each local element patch.
    The exactness of the decomposition (2.6) requires this; for non-affine k the paper relies on nodal averaging and reports convergence, but the exact identity holds only under this assumption.
  • ad hoc to paper The loss of global symmetry when switching to the standard bilinear form at macro interfaces is a higher-order perturbation that does not degrade solver convergence.
    Asserted in Section 2.2 without proof; only empirical evidence from the experiments is provided. This is load-bearing for the practical equivalence of the implemented operator with the symmetric form.

how reviews work

0 comments
Cite this review

Pith. "Pith review of Stencil scaling for vector-valued PDEs on hybrid grids with applications to generalized Newtonian fluids." pith.science (2026). https://pith.science/paper/AM2NOY4J

@misc{pith2026190808666,
  author       = {Pith},
  title        = {Pith review of: Stencil scaling for vector-valued PDEs on hybrid grids with applications to generalized Newtonian fluids},
  year         = {2026},
  howpublished = {\url{https://pith.science/paper/AM2NOY4J}},
  note         = {Machine review of arXiv:1908.08666}
}
read the original abstract

Matrix-free finite element implementations for large applications provide an attractive alternative to standard sparse matrix data formats due to the significantly reduced memory consumption. Here, we show that they are also competitive with respect to the run time in the low order case if combined with suitable stencil scaling techniques. We focus on variable coefficient vector-valued partial differential equations as they arise in many physical applications. The presented method is based on scaling constant reference stencils originating from a linear finite element discretization instead of evaluating the bilinear forms on-the-fly. This method assumes the usage of hierarchical hybrid grids, and it may be applied to vector-valued second-order elliptic partial differential equations directly or as a part of more complicated problems. We provide theoretical and experimental performance estimates showing the advantages of this new approach compared to the traditional on-the-fly integration and stored matrix approaches. In our numerical experiments, we consider two specific mathematical models. Namely, linear elastostatics and incompressible Stokes flow. The final example considers a non-linear shear-thinning generalized Newtonian fluid. For this type of non-linearity, we present an efficient approach to compute a regularized strain rate which is then used to define the node-wise viscosity. Depending on the compute architecture, we could observe maximum speedups of 64% and 122% compared to the on-the-fly integration. The largest considered example involved solving a Stokes problem with 12288 compute cores on the state of the art supercomputer SuperMUC-NG.

Discussion (0). Continue with ORCID to comment.

Reference graph

Works this paper leans on

38 extracted references · 38 canonical work pages

  1. [1]

    A. Afzal, The cost of computation: Metrics and models for modern multicore-based systems in scientific computing, Master’s thesis, Department Informatik, Friedrich Alexander Universität Erlangen-Nürnberg, 2015

  2. [2]

    C. L. Alappat, J. Hofmann, G. Hager, H. Fehske, A. R. Bishop, and G. Wellein , Under- standing hpc benchmark performance on intel broadwell and cascade lake processors, arXiv preprint arXiv:2002.03344, (2020)

  3. [3]

    Arbenz, G

    P. Arbenz, G. H. v an Lenthe, U. Mennel, R. Müller, and M. Sala , A scalable multi-level preconditioner for matrix-freeµ-finite element analysis of human bone structures, International Journal for Numerical Methods in Engineering, 73 (2008), pp. 927–947

  4. [4]

    Bauer, D

    S. Bauer, D. Drzisga, M. Mohr, U. Rüde, C. W aluga, and B. Wohlmuth , A stencil scaling approach for accelerating matrix-free finite element implementations, SIAM Journal on Scientific Computing, 40 (2018), pp. C748–C778

  5. [5]

    Bergen, Hierarchical Hybrid Grids: Data Structures and Core Algorithms for Efficient Finite Element Simulations on Supercomputers, SCS Publishing House, Erlangen, 2005

    B. Bergen, Hierarchical Hybrid Grids: Data Structures and Core Algorithms for Efficient Finite Element Simulations on Supercomputers, SCS Publishing House, Erlangen, 2005

  6. [6]

    Bergen and F

    B. Bergen and F. Hülsemann , Hierarchical hybrid grids: data structures and core algorithms for multigrid, Numer. Lin. Alg. Appl., 11 (2004), pp. 279–291

  7. [7]

    Bergen, G

    B. Bergen, G. Wellein, F. Hülsemann, and U. Rüde , Hierarchical hybrid grids: Achieving TER- AFLOP performance on large scale finite element simulations, International Journal of Parallel, Emergent and Distributed Systems, 22 (2007), pp. 311–329

  8. [8]

    Bey, Tetrahedral grid refinement, Computing, 55 (1995), pp

    J. Bey, Tetrahedral grid refinement, Computing, 55 (1995), pp. 355–378

Show all 38 references
  1. [9]

    Bielak, O

    J. Bielak, O. Ghattas, and E.-J. Kim , Parallel Octree-Based Finite Element Method for Large-Scale Earthquake Ground Motion Simulation, Computer Modeling in Engineering & Sciences, 10 (2005), pp. 99–112

  2. [10]

    Brezzi and J

    F. Brezzi and J. Pitkäranta , On the stabilization of finite element approximations of the Stokes equations, in Efficient solutions of elliptic systems, Springer, 1984, pp. 11–19

  3. [11]

    Brown, Efficient Nonlinear Solvers for Nodal High-Order Finite Elements in 3D, J

    J. Brown, Efficient Nonlinear Solvers for Nodal High-Order Finite Elements in 3D, J. Scientific Computing, 45 (2010), pp. 48–63

  4. [12]

    G. F. Carey and B.-N. Jiang , Element-by-element linear and nonlinear solution schemes, Communica- tions in Applied Numerical Methods, 2 (1986), pp. 145–153

  5. [13]

    R. S. Dembo, S. C. Eisenst a t, and T. Steihaug , Inexact Newton methods, SIAM Journal on Numerical analysis, 19 (1982), pp. 400–408

  6. [14]

    Drzisga, L

    D. Drzisga, L. John, U. Rüde, B. Wohlmuth, and W. Zulehner , On the analysis of block smoothers 29 for saddle point problems, SIAM Journal on Matrix Analysis and Applications, 39 (2018), pp. 932–960

  7. [15]

    Flaig and P

    C. Flaig and P. Arbenz , A Highly Scalable Matrix-Free Multigrid Solver forµFE Analysis Based on a Pointer-Less Octree, in Large-Scale Scientific Computing: 8th International Conference, LSSC 2011, Sozopol, Bulgaria, June 6-10, 2011, Revised Selected Papers, I. Lirkov, S. Marge...

  8. [16]

    G. P. Galdi, R. Rannacher, A. M. Robertson, and S. Turek , Hemodynamical flows, Delhi Book Store, (2008)

  9. [17]

    L. J. Gibson and M. F. Ashby , Cellular solids: structure and properties, Cambridge university press, 1999

  10. [18]

    Gmeiner, T

    B. Gmeiner, T. Gradl, H. Köstler, and U. Rüde , Highly parallel geometric multigrid algorithm for hierarchical hybrid grids, in NIC Symposium, vol. 45, 2012, pp. 323–330

  11. [19]

    Hager and G

    G. Hager and G. Wellein , Introduction to high performance computing for scientists and engineers, CRC Press, 2010

  12. [20]

    A. Ilic, F. Pratas, and L. Sousa , Cache-aware Roofline model: Upgrading the loft, IEEE Computer Architecture Letters, 13 (2013), pp. 21–24

  13. [21]

    https://software.intel.com/en-us/intel-advisor-xe, 2019

    Intel Corp., Intel Advisor. https://software.intel.com/en-us/intel-advisor-xe, 2019

  14. [22]

    https://software.intel.com/en-us/vtune, 2019

    , Intel VTune Profiler. https://software.intel.com/en-us/vtune, 2019

  15. [23]

    N. Kohl, D. Thönnes, D. Drzisga, D. Bar tuscha t, and U. Rüde , The HyTeG finite-element software framework for scalable multigrid solvers, International Journal of Parallel, Emergent and Distributed Systems, (2018), pp. 1–20

  16. [24]

    Kronbichler and K

    M. Kronbichler and K. Kormann , A generic interface for parallel cell-based finite element operator application, Computers and Fluids, 63 (2012), pp. 135–147

  17. [25]

    K. Ljungkvist, Matrix-free Finite-element Computations on Graphics Processors with Adaptively Refined Unstructured Meshes, in Proceedings of the 25th High Performance Computing Symposium, HPC ’17, Society for Computer Simulation International, 2017, pp. 1:1–1:12

  18. [26]

    Ljungkvist and M

    K. Ljungkvist and M. Kronbichler , Multigrid for Matrix-Free Finite Element Computations on Graphics Processors, Tech. Rep. 2017-006, Department of Information Technology, Uppsala University, 2017

  19. [27]

    Loffeld and J

    J. Loffeld and J. Hittinger , On the arithmetic intensity of high-order finite-volume discretizations for hyperbolic systems of conservation laws, The International Journal of High Performance Computing Applications, (2017)

  20. [28]

    https://doku.lrz.de/display/PUBLIC/Hardware+of+SuperMUC-NG (retrieved on 25 February 2020)

    LRZ, Hardware of SuperMUC-NG. https://doku.lrz.de/display/PUBLIC/Hardware+of+SuperMUC-NG (retrieved on 25 February 2020)

  21. [29]

    https://www.lrz.de/services/compute/supermuc/systemdescription/ (retrieved on 29 November 2018)

    , SuperMUC Petascale System. https://www.lrz.de/services/compute/supermuc/systemdescription/ (retrieved on 29 November 2018)

  22. [30]

    D. A. May, J. Brown, and L. L. Pourhiet , A scalable, matrix-free multigrid preconditioner for finite element discretizations of heterogeneous Stokes flow, Computer Methods in Applied Mechanics and Engineering, 290 (2015), pp. 496–523

  23. [31]

    Ricard, Physics of mantle convection, Treatise on Geophysics, 7 (2007), pp

    Y. Ricard, Physics of mantle convection, Treatise on Geophysics, 7 (2007), pp. 31–81

  24. [32]

    U. Rüde, C. W aluga, and B. Wohlmuth , Mass-corrections for the conservative coupling of flow and transport on collocated meshes, Journal of Computational Physics, 305 (2016), p. 319âĂŞ332

  25. [33]

    J. Rudi, G. Stadler, and O. Ghattas , Weighted BFBT preconditioner for Stokes flow problems with highly heterogeneous viscosity, SIAM Journal on Scientific Computing, 39 (2017), pp. S272–S297

  26. [34]

    Stengel, J

    H. Stengel, J. Treibig, G. Hager, and G. Wellein , Quantifying performance bottlenecks of stencil computations using the execution-cache-memory model, inProceedingsofthe29thACMonInternational Conference on Supercomputing, ACM, 2015, pp. 207–216

  27. [35]

    v an Rietbergen, H

    B. v an Rietbergen, H. Weinans, R. Huiskes, and B. Polman , Computational strategies for iterative solutions of large FEM applications employing voxel data, International Journal for Numerical Methods in Engineering, 39 (1996), pp. 2743–2767

  28. [36]

    Williams, A

    S. Williams, A. W a terman, and D. P a tterson, Roofline: An Insightful Visual Performance Model for Multicore Architectures, Commun. ACM, 52 (2009), pp. 65–76

  29. [37]

    Zhang, D

    Y. Zhang, D. Rodrigue, and A. Ait-Kadi , High density polyethylene foams. ii. elastic modulus, Journal of applied polymer science, 90 (2003), pp. 2120–2129

  30. [38]

    O. C. Zienkiewicz and J. Z. Zhu , The superconvergent patch recovery and a posteriori error estimates. part 1: The recovery technique, International Journal for Numerical Methods in Engineering, 33 (1992), pp. 1331–1364. 30

Pith tools

Reviewed August 14, 2026 · model on record in the stance chip above.