arxiv: 2605.14035 · v1 · submitted 2026-05-13 · 🧮 math.NA · cs.NA

Recognition: 2 theorem links

· Lean Theorem

ellFEM: An efficient loop-free Matlab implementation of isoparametric bulk and surface finite elements

Bal\'azs Kov\'acs , Michael Lantelme

Authors on Pith no claims yet

Pith reviewed 2026-05-15 02:16 UTC · model grok-4.3

classification 🧮 math.NA cs.NA

keywords finiteelementsimplementationisoparametricmatlabbulkefficientloop-free

0 comments

The pith

ℓFEM is a loop-free MATLAB package implementing isoparametric bulk and surface finite elements with high-order support, assembly details, and performance tests.

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

This paper introduces a MATLAB software package called ℓFEM that computes finite element solutions for problems inside volumes and on curved surfaces. It focuses on isoparametric elements, which use the same polynomial mapping for both geometry and solution fields, allowing accurate representation of curved boundaries at high order. The key implementation choice is to avoid explicit loops by relying on MATLAB's paged array operations, which vectorize the assembly of stiffness matrices and load vectors. The authors outline the surface element construction, the quadrature rules used, and how the code handles both linear and nonlinear problems. They include runtime comparisons against other approaches and supply a testing suite plus examples to verify the implementation.

Core claim

The ℓFEM MATLAB package provides a simple, efficient, and flexible implementation of isoparametric finite elements in bulk domains and on surfaces with completely loop-free matrix assemblies based on paged operators.

Load-bearing premise

That MATLAB's paged operators deliver competitive performance for the described high-order isoparametric surface elements without hidden overheads that would negate the loop-free advantage.

Figures

Figures reproduced from arXiv: 2605.14035 by Bal\'azs Kov\'acs, Michael Lantelme.

**Figure 2.** Figure 2: The reference element Eˆ and an arbitrary element E, and the transformation maps F and F −1 in between. We will denote the global basis functions ϕj , while on the reference element Eˆ by φˆj , where ˆj is the local label of the globally labelled node j. Recall the Lagrange basis functions satisfy ϕj (xi) = δji. Any element E, with nodes xj , is also given as the image of the map (3.1) F(x) = X 6 ˆj=1 xjφˆ… view at source ↗

**Figure 3.** Figure 3: 2D P2 bulk elements E (Elements) 1 2 3 4 5 6 2 7 3 8 9 5 [PITH_FULL_IMAGE:figures/full_fig_p007_3.png] view at source ↗

**Figure 4.** Figure 4: Coordinate array N, where each page corresponds to some element Ei for all elements (i = 1, . . . , |E|). The vector x Ei j stores all nodal coordinates of the j-th spatial dimension of the nodes of the given Element. O(|E| 2 ) in the number of elements |E|. There is a remedy to the elementwise approach, by storing the contributions of each element in cell arrays—storing the local mass and stiffness matrix… view at source ↗

**Figure 5.** Figure 5: Assembly of mass, stiffness matrices and load vector for 2D surface FEM [PITH_FULL_IMAGE:figures/full_fig_p013_5.png] view at source ↗

**Figure 6.** Figure 6: Assembly of mass, stiffness matrices and load vector for 2D bulk FEM with [PITH_FULL_IMAGE:figures/full_fig_p014_6.png] view at source ↗

**Figure 7.** Figure 7: Comparing the assembly times for the non-vectorized, [PITH_FULL_IMAGE:figures/full_fig_p016_7.png] view at source ↗

read the original abstract

The $\ell$FEM MATLAB package provides a simple, efficient, and flexible implementation of isoparametric finite elements in bulk domains and on surfaces. The finite element matrix assemblies are based on MATLAB's paged operators and therefore completely loop-free. We give a short and conscious description of high-order isoparametric surface finite elements, which is then used to describe the assembly process and the implementation. We report on relevant numerical experiments (runtime comparisons, modifications for non-linear problems, etc.), and on additional functions, examples, and a testing unit which are all part of the $\ell$FEM package.

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.

Desk Editor's Note private letter to a colleague

ℓFEM gives MATLAB users a documented loop-free assembly for high-order isoparametric bulk and surface elements via paged operators, with the main value in the practical implementation rather than new theory.

read the letter

The paper's real contribution is a clear, loop-free MATLAB implementation of isoparametric finite elements that handles both volume and surface cases. It uses paged operators to avoid explicit loops during matrix assembly, and it includes a concise description of the high-order surface element setup that feeds directly into the code structure. The authors also supply runtime comparisons, notes on nonlinear extensions, examples, and a test suite, which makes the package immediately usable for people already working in MATLAB on surface PDEs. That combination of documented assembly and supporting material is the part that stands out as useful implementation work. The experiments they report are the right way to back the efficiency claim, and if those timings show consistent gains over standard looped code, the approach delivers on its promise. The stress-test worry about memory overheads from strided access in paged tensor contractions for curved high-order elements is reasonable to raise, but it lands only if the paper's own runtime data fails to show net improvement; otherwise the measured results take precedence. No new mathematical theory is claimed, and the citation pattern stays within the expected references for isoparametric surface FEM, so nothing looks circular or invented. This is a tool paper aimed at numerical analysts who need quick, maintainable code for these elements rather than a broad audience. Readers who already code in MATLAB and care about surface problems will get direct value from the package and the assembly description. It is solid enough on its own terms to deserve peer review so that the runtime claims can be checked against the released code and any edge cases in high-order mappings can be discussed.

Referee Report

1 major / 2 minor

Summary. The manuscript presents the ℓFEM MATLAB package for isoparametric finite elements in bulk domains and on surfaces. Matrix assemblies are implemented in a completely loop-free manner via MATLAB paged operators. The paper gives a description of high-order isoparametric surface finite elements, details the assembly process, reports runtime experiments (including comparisons and modifications for nonlinear problems), and describes additional package components such as functions, examples, and a testing unit.

Significance. If the reported runtime advantages hold under scrutiny, the package would be a useful practical contribution to the numerical analysis community. It lowers the barrier to using high-order surface finite elements in MATLAB by replacing explicit loops with native vectorized operations, which aids reproducibility and rapid prototyping for geometric PDEs. The inclusion of a testing unit and examples is a positive feature for usability.

major comments (1)

[Numerical experiments] Numerical experiments section: the runtime comparisons for high-order isoparametric surface elements report speedups but omit hardware specifications, MATLAB version, and any memory profiling (e.g., temporary allocations or cache behavior during paged tensor contractions). This information is required to evaluate whether the paged-operator approach truly avoids the overheads that could negate its advantage over hand-tuned looped code for curved high-order mappings.

minor comments (2)

[Abstract] Abstract: 'conscious description' appears to be a typographical error for 'concise description'.
[Assembly process] The description of the assembly process would benefit from an explicit small-scale example (e.g., a single quadratic surface element) showing the exact paged-operator expressions used for the isoparametric mapping and stiffness matrix contributions.

Simulated Author's Rebuttal

1 responses · 0 unresolved

We thank the referee for the positive evaluation of the manuscript and the constructive suggestion regarding the numerical experiments. We address the comment below and will revise the manuscript accordingly.

read point-by-point responses

Referee: [Numerical experiments] Numerical experiments section: the runtime comparisons for high-order isoparametric surface elements report speedups but omit hardware specifications, MATLAB version, and any memory profiling (e.g., temporary allocations or cache behavior during paged tensor contractions). This information is required to evaluate whether the paged-operator approach truly avoids the overheads that could negate its advantage over hand-tuned looped code for curved high-order mappings.

Authors: We agree that these details are important for a complete assessment of the reported speedups. In the revised manuscript we will add the hardware specifications of the test machine (CPU model, core count, and RAM), the precise MATLAB version used, and basic memory profiling results obtained via MATLAB's built-in memory monitoring (peak memory usage and notes on temporary array allocations during paged tensor operations). A detailed cache-miss analysis lies outside the scope of the current experiments, but the added information will allow readers to judge the practical overhead of the paged-operator approach. revision: yes

Axiom & Free-Parameter Ledger

0 free parameters · 1 axioms · 0 invented entities

The implementation rests on standard finite-element theory without introducing new free parameters, axioms beyond established mathematics, or invented entities.

axioms (1)

standard math Standard isoparametric finite-element mapping and quadrature rules hold for both bulk and surface elements
Invoked when describing the high-order surface element construction and assembly process.

pith-pipeline@v0.9.0 · 5404 in / 1089 out tokens · 39194 ms · 2026-05-15T02:16:46.123243+00:00 · methodology

discussion (0)

Lean theorems connected to this paper

Citations machine-checked in the Pith Canon. Every link opens the source theorem in the public Lean library.

IndisputableMonolith/Cost/FunctionalEquation.lean washburn_uniqueness_aczel unclear

?

unclear
Relation between the paper passage and the cited Recognition theorem.

The finite element matrix assemblies are based on MATLAB's paged operators and therefore completely loop-free... assembly functions can readily be extended to higher dimension, or to finite elements of polynomial degree k≥3.
IndisputableMonolith/Foundation/AbsoluteFloorClosure.lean reality_from_one_distinction unclear

?

unclear
Relation between the paper passage and the cited Recognition theorem.

We use the classical reference element technique to assemble finite element matrices, but it is completely loop-free... operations are automatically parallel, and allow for computations on a GPU by a simple command.

What do these tags mean?

matches: The paper's claim is directly supported by a theorem in the formal canon.
supports: The theorem supports part of the paper's argument, but the paper may add assumptions or extra steps.
extends: The paper goes beyond the formal theorem; the theorem is a base layer rather than the whole result.
uses: The paper appears to rely on the theorem as machinery.
contradicts: The paper's claim conflicts with a theorem or certificate in the canon.
unclear: Pith found a possible connection, but the passage is too broad, indirect, or ambiguous to say the theorem truly supports the claim.

Reference graph

Works this paper leans on

40 extracted references · 40 canonical work pages

[1]

Alberty, C

J. Alberty, C. Carstensen, and S. A. Funken , Remarks around 50 lines of Matlab: short finite element implementation, Numer. Algorithms, 20 (1999), pp. 117–137, https://doi.or g/10.1023/A:1019155918070

work page doi:10.1023/a:1019155918070 1999
[2]

Anderson, J

R. Anderson, J. Andrej, A. Barker, J. Bramwell, J.-S. Camier, J. Cerveny, V. Dobrev, Y. Dudouit, A. Fisher, T. Kolev, W. Pazner, M. Stowell, V. Tomov, I. Akkerman, J. Dahm, D. Medina, and S. Zampini , Mfem: A modular finite element methods library, Computers & Mathematics with Applications, 81 (2021), pp. 42–74, https://doi.org/https: //doi.org/10.1016/j....

work page doi:10.1016/j.camwa.2020.06.009 2021
[3]

J. W. Barrett, H. Garcke, and R. Nürnberg , Parametric finite element approximations of curvature-driven interface evolutions, in Geometric partial differential equations. Part I, vol. 21 of Handb. Numer. Anal., Elsevier/North-Holland, Amsterdam, [2020]©2020, pp. 275–423

work page 2020
[4]

Bartels, C

S. Bartels, C. Carstensen, and A. Hecht , P2Q2Iso2D = 2D isoparametric FEM in Matlab, J. Comput. Appl. Math., 192 (2006), pp. 219–250

work page 2006
[5]

Bastian, M

P. Bastian, M. Blatt, A. Dedner, C. Engwer, R. Klöfkorn, R. Kornhuber, M. Ohlberger, and O. Sander , A generic grid interface for parallel and adaptive scientific computing. part ii: Implementation and tests in dune, Computing, 82 (2008), pp. 121–138, https://doi.org/10.1007/s00607-008-0004-9

work page doi:10.1007/s00607-008-0004-9 2008
[6]

Beuter and S

S. Beuter and S. A. Funken , Efficient P1-FEM for any space dimension in Matlab, Comput. Methods Appl. Math., 24 (2024), pp. 283–315, https://doi.org/10.1515/cmam-2022-0239

work page doi:10.1515/cmam-2022-0239 2024
[7]

Bonito, R

A. Bonito, R. A. DeVore, and R. H. Nochetto , Adaptive finite element methods for elliptic problems with discontinuous coefficients, SIAM J. Numer. Anal., 51 (2013), pp. 3106–3134, https://doi.org/10.1137/130905757

work page doi:10.1137/130905757 2013
[8]

Braess , Finite Elements: Theory, Fast Solvers, and Applications in Solid Mechanics, Cambridge University Press, 3 ed., 2007

D. Braess , Finite Elements: Theory, Fast Solvers, and Applications in Solid Mechanics, Cambridge University Press, 3 ed., 2007

work page 2007
[9]

Bullerjahn and B

N. Bullerjahn and B. Kovács , Error estimates for full discretization of Cahn–Hilliard equation with dynamic boundary conditions, IMA J. Numer. Anal., 46 (2026), pp. 713–757, https://doi.org/10.1093/imanum/draf009

work page doi:10.1093/imanum/draf009 2026
[10]

Chen , ifem: an integrated finite element methods package in MATLAB, tech

L. Chen , ifem: an integrated finite element methods package in MATLAB, tech. report, University of California at Irvine, 2009, https://github.com/lyc102/ifem. Technical Report

work page 2009
[11]

Croucher, Paged matrix functions in matlab (2024 edition), May 29 2024, https://blog s.mathworks.com/matlab/2024/05/29/paged-matrix-functions-in-matlab-2024-edition/

M. Croucher, Paged matrix functions in matlab (2024 edition), May 29 2024, https://blog s.mathworks.com/matlab/2024/05/29/paged-matrix-functions-in-matlab-2024-edition/. MATLAB & Simulink - MathWorks Blogs

work page 2024
[12]

Cuvelier, C

F. Cuvelier, C. Japhet, and G. Scarella , An efficient way to perform the assembly of finite element matrices in matlab and octave, (2013)

work page 2013
[13]

Deckelnick, G

K. Deckelnick, G. Dziuk, and C. M. Elliott , Computation of geometric partial differential equations and mean curvature flow, Acta Numerica, 14 (2005), pp. 139–232

work page 2005
[14]

Deckelnick, H

K. Deckelnick, H. Garcke, and B. Kovács , Error estimates for a surface finite element method for anisotropic mean curvature flow, (2025). arxiv:2508.01078

work page arXiv 2025
[15]

Demlow, Higher–order finite element methods and pointwise error estimates for elliptic problems on surfaces, SIAM J

A. Demlow, Higher–order finite element methods and pointwise error estimates for elliptic problems on surfaces, SIAM J. Numer. Anal., 47 (2009), pp. 805–807, https://doi.org/10.1 137/070708135

work page 2009
[16]

Demlow and G

A. Demlow and G. Dziuk , An adaptive finite element method for the laplace–beltrami operator on implicitly defined surfaces, SIAM Journal on Numerical Analysis, 45 (2007), pp. 421–442, https://doi.org/10.1137/050642873

work page doi:10.1137/050642873 2007
[17]

Dunavant, High degree efficient symmetrical Gaussian quadrature rules for the triangle, International journal for numerical methods in engineering, 21 (1985), pp

D. Dunavant, High degree efficient symmetrical Gaussian quadrature rules for the triangle, International journal for numerical methods in engineering, 21 (1985), pp. 1129–1148

work page 1985
[18]

G. Dziuk, Finite elements for the Beltrami operator on arbitrary surfaces, Partial differential equations and calculus of variations, Lecture Notes in Math., 1357, Springer, Berlin, (1988), pp. 142–155

work page 1988
[19]

Dziuk, An algorithm for evolutionary surfaces, Numer

G. Dziuk, An algorithm for evolutionary surfaces, Numer. Math., 58 (1990), pp. 603–611

work page 1990
[20]

Dziuk and C

G. Dziuk and C. M. Elliott , Finite element methods for surface PDEs, Acta Numerica, 22 (2013), pp. 289–396., https://doi.org/10.1017/S0962492913000056

work page doi:10.1017/s0962492913000056 2013
[21]

Edelmann, B

D. Edelmann, B. Kovács, and C. Lubich , Numerical analysis of an evolving bulk–surface model of tumour growth, IMA Journal of Numerical Analysis, (2024), https://doi.org/10.1 093/imanum/drae077

work page 2024
[22]

C. M. Elliott, H. Garcke, and B. Kovács , Numerical analysis for the interaction of mean curvature flow and diffusion on closed surfaces, Numerische Mathematik, 151 (2022), pp. 873–925. ℓFEM: LOOP-FREE BULK AND SURFACE FEM ASSEMBLY 19

work page 2022
[23]

C. M. Elliott and T. Ranner , Finite element analysis for a coupled bulk-surface partial differential equation, IMA J. Numer. Anal., 33 (2013), pp. 377–402

work page 2013
[24]

C. M. Elliott and T. Ranner , A unified theory for continuous-in-time evolving finite element space approximations to partial differential equations in evolving domains, IMA Journal of Numerical Analysis, 41 (2021), pp. 1696–1845

work page 2021
[25]

Ern and J.-L

A. Ern and J.-L. Guermond , Theory and Practice of Finite Elements, vol. 159 of Applied Mathematical Sciences, Springer, New York, 2004, https://doi.org/10.1007/978-1-4757-435 5-5

work page doi:10.1007/978-1-4757-435 2004
[26]

Eyles, J

J. Eyles, J. R. King, and V. Styles , A tractable mathematical model for tissue growth, Interfaces Free Bound., 21 (2019), pp. 463–493, https://doi.org/10.4171/ifb/428

work page doi:10.4171/ifb/428 2019
[27]

Feifel and S

S. Feifel and S. A. Funken ,Efficient p1-fem for any space dimension in matlab, Computational Methods in Applied Mathematics, 24 (2024), pp. 283–324, https://doi.org/10.1515/cmam-2 022-0239

work page doi:10.1515/cmam-2 2024
[28]

Funken, D

S. Funken, D. Praetorius, and P. Wissgott , Efficient implementation of adaptive P1-FEM in Matlab, Comput. Methods Appl. Math., 11 (2011), pp. 460–490, https://doi.org/10.247 8/cmam-2011-0026

work page 2011
[29]

D. Hipp, M. Hochbruck, and C. Stohrer , Unified error analysis for nonconforming space discretizations of wave-type equations, IMA J. Numer. Anal., 39 (2019), pp. 1206–1245, https://doi.org/10.1093/imanum/dry036

work page doi:10.1093/imanum/dry036 2019
[30]

Jaśkowiec and N

J. Jaśkowiec and N. Sukumar , High-order symmetric cubature rules for tetrahedra and pyramids, Internat. J. Numer. Methods Engrg., 122 (2021), pp. 148–171, https://doi.org/10 .1002/nme.6528

work page 2021
[31]

Kovács, Numerical surgery for mean curvature flow of surfaces, SIAM Journal on Scientific Computing, 46 (2024), pp

B. Kovács, Numerical surgery for mean curvature flow of surfaces, SIAM Journal on Scientific Computing, 46 (2024), pp. A645–A669, https://doi.org/10.1137/22M1531919

work page doi:10.1137/22m1531919 2024
[32]

Kovács and M

B. Kovács and M. Lantelme , A posteriori error estimates and space–time adaptivity for parabolic partial differential equations on stationary surfaces, to appear in SIAM J. Numer. Anal., (2026). arXiv:2407.02101

work page arXiv 2026
[33]

Kovács, B

B. Kovács, B. Li, and C. Lubich , A convergent evolving finite element algorithm for mean cur- vature flow of closed surfaces, Numer. Math., 143 (2019), pp. 797–853. DOI:10.1007/s00211- 019-01074-2

work page doi:10.1007/s00211- 2019
[34]

Kovács, B

B. Kovács, B. Li, and C. Lubich , A convergent evolving finite element algorithm for mean curvature flow of closed surfaces, Numer. Math., 143 (2019), pp. 797–853, https://doi.org/ 10.1007/s00211-019-01074-2

work page doi:10.1007/s00211-019-01074-2 2019
[35]

Kovács and C

B. Kovács and C. Lubich , Numerical analysis of parabolic problems with dynamic boundary conditions, IMA J. Numer. Anal., 37 (2017), pp. 1–39, https://doi.org/10.1093/imanum/d rw015

work page doi:10.1093/imanum/d 2017
[36]

H. N. Mallesham, S. Gaddam, J. V aldman, and S. K. Acharya , Vectorized 3d mesh refinement and implementation of primal hybrid fem in matlab, 2025, https://arxiv.org/ab s/2509.11133

work page arXiv 2025
[37]

Persson and G

P.-O. Persson and G. Strang , A simple mesh generator in MATLAB, SIAM Review, 46 (2004), pp. 329–345

work page 2004
[38]

Rahman and J

T. Rahman and J. V aldman , Fast MATLAB assembly of FEM matrices in 2D and 3D: nodal elements, Appl. Math. Comput., 219 (2013), pp. 7151–7158, https://doi.org/10.1016/j.amc. 2011.08.043

work page doi:10.1016/j.amc 2013
[39]

Schmidt and K

A. Schmidt and K. Siebert , Design of adaptive finite element software: The finite element toolbox alberta, (2005)

work page 2005
[40]

Styles and J

V. Styles and J. V an Yperen , Numerical analysis for a system coupling curve evolution attached orthogonally to a fixed boundary, to a reaction-diffusion equation on the curve, Num. Meth. PDE, (2021)

work page 2021