Expand JOSS paper: add sections, figure, references, impact, AI disclosure

Grow paper.md to a full JOSS-style article (state of the field, software design incl. spectral hierarchy, research impact, AI usage disclosure, acknowledgements/COI) with a rendered figure; expand paper.bib to 16 references with verified DOIs; gitignore the generated paper.pdf.

Author: Sébastien Loisel

.gitignore

@@ -36,3 +36,6 @@ CLAUDE.md
 
 # Development scratch folder
 dev/
+
+# JOSS paper PDF is a build artifact (generated by the Open Journals action)
+paper/paper.pdf

paper/figure.png

@@ -22,6 +22,71 @@ @article{loisel2020efficient
   doi       = {10.1007/s00211-020-01141-z}
 }
 
+@unpublished{loiselmgb,
+  author    = {Loisel, Sébastien},
+  title     = {Algorithm {MGB} to solve highly nonlinear elliptic {PDEs} in $\tilde{O}(n)$ {FLOPS}},
+  year      = {2026},
+  note      = {Submitted}
+}
+
+@unpublished{loiselnonuniform,
+  author    = {Loisel, Sébastien},
+  title     = {The Multi-Grid-Barrier method on nonuniform grids},
+  year      = {2026},
+  note      = {Submitted}
+}
+
+@unpublished{loiselhpc,
+  author    = {Loisel, Sébastien},
+  title     = {A high-performance solver for convex optimization problems in function spaces},
+  year      = {2026},
+  note      = {Submitted}
+}
+
+@inproceedings{loiseldd29,
+  author    = {Loisel, Sébastien},
+  title     = {The Algebraic Multi-Grid-Barrier method for solving $p$-Laplace and other convex optimization problems},
+  booktitle = {Proceedings of the 29th International Conference on Domain Decomposition Methods, Milan, Italy},
+  publisher = {Springer},
+  year      = {2026},
+  note      = {Accepted, to appear}
+}
+
+@article{zhangjiang2025cnn,
+  author    = {Zhang, Xiaoyu and Jiang, Lijian},
+  title     = {Convolutional neural network based reduced order modeling for multiscale problems},
+  journal   = {Journal of Computational Physics},
+  volume    = {524},
+  pages     = {113710},
+  year      = {2025},
+  publisher = {Elsevier},
+  doi       = {10.1016/j.jcp.2024.113710}
+}
+
+@article{muller2021pharmonic,
+  author    = {M{\"u}ller, Peter Marvin and K{\"u}hl, Niklas and Siebenborn, Martin and Deckelnick, Klaus and Hinze, Michael and Rung, Thomas},
+  title     = {A novel p-harmonic descent approach applied to fluid dynamic shape optimization},
+  journal   = {Structural and Multidisciplinary Optimization},
+  volume    = {64},
+  number    = {6},
+  pages     = {3489--3503},
+  year      = {2021},
+  publisher = {Springer},
+  doi       = {10.1007/s00158-021-03030-x}
+}
+
+@article{balci2023kacanov,
+  author    = {Balci, Anna Kh. and Diening, Lars and Storn, Johannes},
+  title     = {Relaxed {Ka{\v{c}}anov} scheme for the $p$-Laplacian with large exponent},
+  journal   = {SIAM Journal on Numerical Analysis},
+  volume    = {61},
+  number    = {6},
+  pages     = {2775--2794},
+  year      = {2023},
+  publisher = {SIAM},
+  doi       = {10.1137/22M1528550}
+}
+
 @article{bezanson2017julia,
   author    = {Bezanson, Jeff and Edelman, Alan and Karpinski, Stefan and Shah, Viral B.},
   title     = {Julia: A fresh approach to numerical computing},
@@ -34,6 +99,35 @@ @article{bezanson2017julia
   doi       = {10.1137/141000671}
 }
 
+@article{gridap,
+  author    = {Badia, Santiago and Verdugo, Francesc},
+  title     = {Gridap: An extensible finite element toolbox in {Julia}},
+  journal   = {Journal of Open Source Software},
+  volume    = {5},
+  number    = {52},
+  pages     = {2520},
+  year      = {2020},
+  doi       = {10.21105/joss.02520}
+}
+
+@article{fenics,
+  author    = {Aln{\ae}s, Martin S. and Blechta, Jan and Hake, Johan and Johansson, August and Kehlet, Benjamin and Logg, Anders and Richardson, Chris and Ring, Johannes and Rognes, Marie E. and Wells, Garth N.},
+  title     = {The {FEniCS} project version 1.5},
+  journal   = {Archive of Numerical Software},
+  volume    = {3},
+  number    = {100},
+  pages     = {9--23},
+  year      = {2015},
+  doi       = {10.11588/ans.2015.100.20553}
+}
+
+@misc{ferrite,
+  author       = {Carlsson, Kristoffer and Ekre, Fredrik and {Ferrite.jl contributors}},
+  title        = {{Ferrite.jl}: A finite element toolbox for {Julia}},
+  year         = {2024},
+  doi          = {10.5281/zenodo.18196628}
+}
+
 @misc{AlgebraicMultigrid,
   author       = {{JuliaLinearAlgebra contributors}},
   title        = {{AlgebraicMultigrid.jl}: Algebraic multigrid (AMG) solvers in {Julia}},

paper/paper.bib

@@ -22,30 +22,28 @@ bibliography: paper.bib
 
 # Summary
 
-`MultiGridBarrier.jl` is a Julia package for solving convex variational problems in
-function spaces. These are the nonlinear partial differential equations (PDEs) and
-boundary-value problems that arise from minimizing a convex functional. Representative examples include
-the $p$-Laplacian for any $p \in [1, \infty]$, total-variation problems, and obstacle
-problems. The most useful of these are *nonsmooth*: the energy is convex but not
-differentiable (e.g. $p = 1$ or total variation), a regime in which Newton-type solvers
-applied naively either fail or require a number of iterations that grows rapidly with the
-mesh resolution.
-
-The package implements the **multigrid barrier method**, which couples an interior-point
-(barrier) method with a multigrid hierarchy. For the problem classes covered by the
-supporting theory the method is *quasi-optimal*: the number of interior-point/Newton
-iterations grows only mildly with the number of degrees of freedom $n$. For instance, this
-count is $O(\sqrt{n}\,\log n)$ for the $p$-Laplacian [@loisel2020efficient], and
-polylogarithmic in the analytic, spectral setting [@loisel2026spectral].
+`MultiGridBarrier.jl` is a Julia [@bezanson2017julia] package for solving convex variational
+problems in function spaces. These are the nonlinear partial differential equations (PDEs) and boundary-value
+problems that arise from minimizing a convex functional. Representative examples include the
+$p$-Laplacian for any $p \in [1, \infty]$, total-variation problems, and obstacle problems.
+The most useful of these are *nonsmooth*: the energy is convex but not differentiable (for
+example $p = 1$, or total variation), a regime in which Newton-type solvers applied naively
+either fail or require a number of iterations that grows rapidly with the mesh resolution.
+
+The package implements the **multigrid barrier method** [@loiselmgb; @loisel2020efficient;
+@loisel2026spectral], which couples an interior-point (barrier) method with a multigrid
+hierarchy. For the problem classes covered by the supporting theory the method is
+*quasi-optimal*: the number of interior-point/Newton iterations grows only mildly with the
+number of degrees of freedom $n$. For instance, this count is $O(\sqrt{n}\,\log n)$ for the
+$p$-Laplacian [@loisel2020efficient], and polylogarithmic in the analytic, spectral setting
+[@loisel2026spectral].
 
 `MultiGridBarrier.jl` provides finite-element discretizations in one, two, and three
 dimensions (simplicial $P_1$/$P_2$ elements and tensor-product $Q_k$ elements), as well as
 Chebyshev spectral discretizations, all with isoparametric element maps. It builds an
-algebraic-multigrid hierarchy automatically (via `AlgebraicMultigrid.jl`
-[@AlgebraicMultigrid] or, optionally, `PyAMG` [@pyamg]), supports user-specified mesh
-connectivity (enabling slit domains, branch cuts, and glued manifolds), solves
-time-dependent problems, and offers optional GPU acceleration through CUDA. A typical solve
-is three lines:
+algebraic-multigrid hierarchy automatically, supports user-specified mesh connectivity
+(enabling slit domains, branch cuts, and glued manifolds), solves time-dependent problems,
+and offers optional GPU acceleration through CUDA. A typical solve is three lines:
 
 ```julia
 using MultiGridBarrier
@@ -58,38 +56,91 @@ sol  = mgb_solve(assemble(amg(geom); p = 1.0))   # a nonsmooth p = 1 problem
 Convex variational problems are ubiquitous in computational science: nonlinear elasticity
 and plasticity, image denoising and segmentation (total variation), contact and obstacle
 problems, and non-Newtonian flow (the $p$-Laplacian). The difficulty is that the most
-interesting cases are nonsmooth (the energy is convex but not differentiable), so
-Newton-type methods applied naively either stagnate or require an iteration count that grows
-rapidly as the mesh is refined.
+interesting cases are nonsmooth (the energy is convex but not differentiable), so Newton-type
+methods applied naively either stagnate or require an iteration count that grows rapidly as
+the mesh is refined.
 
 Interior-point (barrier) methods handle nonsmoothness robustly by following a smooth central
 path, but a single barrier solve still requires solving a sequence of large, increasingly
-ill-conditioned linear systems. The multigrid barrier method addresses both issues at once:
-a multigrid hierarchy preconditions the central-path subproblems so that the *total* cost
-stays close to linear in the number of unknowns, with rigorous bounds for the covered problem
-classes [@loisel2020efficient; @loisel2026spectral].
-
-General-purpose finite-element libraries and algebraic-multigrid libraries in Julia provide
-the building blocks for discretizing PDEs and solving linear systems, but they do not provide
-an out-of-the-box, theoretically grounded solver for nonsmooth convex variational problems.
-`MultiGridBarrier.jl` fills this gap: it packages the discretization, the multigrid hierarchy,
-and the barrier solver behind a small high-level interface (`fem2d_P2`, `amg`, `assemble`,
-`mgb_solve`). Researchers and practitioners can then solve such problems, and reproduce the
-numerical results of the underlying papers, in a few lines on the CPU or the GPU.
-
-# Functionality
-
-- **Discretizations.** Finite elements in 1D/2D/3D: simplicial $P_1$/$P_2$ and tensor-product
-  $Q_k$, plus Chebyshev spectral elements; all isoparametric.
-- **Solver.** An algebraic-multigrid hierarchy (`amg`) drives a barrier (interior-point)
-  method (`mgb_solve`) for user-assembled convex problems (`assemble`).
-- **Convex constraints.** Built-in convex sets for $p$-norm/Euclidian-power, linear, and
-  piecewise constraints, composable via `intersect`.
-- **Topological meshes.** Explicit connectivity (`tensor_dofmap` and the `t=` keyword) lets
-  geometrically coincident nodes remain topologically distinct, supporting slit domains,
-  branch cuts, and glued manifolds.
-- **Time dependence.** `parabolic_solve` for time-dependent problems.
-- **GPU.** Optional CUDA acceleration through a package extension.
-- **Visualization.** Plotting of 1D/2D/3D solutions and animations.
+ill-conditioned linear systems. The multigrid barrier method addresses both issues at once: a
+multigrid hierarchy preconditions the central-path subproblems so that the *total* cost stays
+close to linear in the number of unknowns, with rigorous bounds for the covered problem
+classes [@loiselmgb; @loisel2020efficient; @loisel2026spectral].
+
+# State of the field
+
+General-purpose finite-element libraries in Julia, such as `Gridap.jl` [@gridap] and
+`Ferrite.jl` [@ferrite], and in the wider ecosystem, such as the FEniCS project [@fenics],
+provide flexible tools for discretizing PDEs; algebraic-multigrid libraries such as
+`AlgebraicMultigrid.jl` [@AlgebraicMultigrid] and `PyAMG` [@pyamg] solve the resulting linear
+systems. These are general building blocks, but none provides an out-of-the-box,
+theoretically grounded solver for *nonsmooth convex variational* problems.
+`MultiGridBarrier.jl` fills this gap. It packages a discretization, a multigrid hierarchy, and
+a barrier solver behind a small high-level interface, and builds on (rather than reinvents)
+the Julia ecosystem: it uses `AlgebraicMultigrid.jl`, or optionally `PyAMG`, to coarsen its
+auxiliary problems, and runs on both CPU and GPU.
+
+# Software design
+
+A problem is solved with four composable steps. A *mesh constructor* (`fem1d`, `fem2d`,
+`fem2d_P1`, `fem2d_P2`, `fem3d`, `spectral1d`, `spectral2d`) returns a `Geometry` describing
+the discretization. `amg(geom)` attaches an algebraic-multigrid hierarchy, `assemble(mg; ...)`
+builds the convex problem (the functional, its barrier, and any constraints), and
+`mgb_solve(prob)` runs the barrier method.
+
+Internally, the multigrid barrier method tracks the central path of an interior-point method
+while using a multigrid hierarchy to solve the Newton systems along that path
+[@loiselmgb; @loisel2026spectral]. The finite-element discretizations are isoparametric
+simplicial $P_k$ and tensor-product $Q_k$ elements. For these, the multigrid hierarchy is built
+algebraically: the package coarsens an auxiliary $P_1$/$Q_1$ problem on the element corners and
+lifts the resulting transfer operators to the full high-order basis, so the same machinery
+serves every finite-element family. Mesh topology is represented by an explicit node
+connectivity array, which decouples geometry from topology and lets geometrically coincident
+nodes remain distinct, supporting slit domains, branch cuts, and glued manifolds. The package
+also provides a structured, batched-GEMM assembly of the Newton Hessians that maps efficiently
+onto GPUs [@loiselhpc], a variant for nonuniform grids [@loiselnonuniform], and the algebraic
+formulation used here [@loiseldd29].
+
+The Chebyshev spectral discretizations (`spectral1d`, `spectral2d`) use an intrinsically
+spectral hierarchy instead. The multigrid levels are a sequence of polynomial approximation
+spaces of increasing degree, with exact polynomial interpolation as the inter-level transfer
+and zero-trace boundary conditions imposed by basis construction rather than node masking. This
+is the setting of the spectral barrier method [@loisel2026spectral]: when the solution is
+analytic, the spectral discretization converges geometrically and the overall method is
+quasi-optimal, with the iteration count growing only polylogarithmically in the number of
+degrees of freedom.
+
+Convex constraints (p-norm / Euclidian-power, linear, and piecewise) are built in and
+composable for any discretization, and `parabolic_solve` extends the method to time-dependent
+problems.
+
+![A nonsmooth ($p = 1$) solution of a two-dimensional $p$-Laplace problem, computed with
+`MultiGridBarrier.jl` on a refined $Q_2$ mesh and rendered with the package's plotting
+front-end.](figure.png)
+
+# Research impact
+
+The methods implemented in `MultiGridBarrier.jl` underpin the numerical experiments in the
+papers that introduce them [@loisel2020efficient; @loiselmgb; @loisel2026spectral], and the
+solver and its earlier implementations have been used by other researchers. For example, Zhang
+and Jiang use this algorithm within a convolutional-neural-network reduced-order modeling
+method for multiscale problems [@zhangjiang2025cnn]. More broadly, the underlying $p$-Laplacian
+algorithm [@loisel2020efficient] has been cited around 26 times (Google Scholar, as of June
+2026) across the numerical-PDE and optimization literature, in areas such as computational
+$p$-Laplacian numerics [@balci2023kacanov] and $p$-harmonic shape optimization
+[@muller2021pharmonic].
+
+# AI usage disclosure
+
+Generative AI assistance (Anthropic Claude, via the Claude Code command-line tool) was used in
+preparing this submission: drafting and editing this paper and parts of the documentation and
+README, generating the illustrative figure, and assisting with some software changes
+(refactoring, a correctness fix, and a mesh-connectivity feature). All AI-assisted output was
+reviewed, validated, and edited by the author, who made all core design and research decisions
+and takes full responsibility for the software and the manuscript.
+
+# Acknowledgements
+
+This work received no specific external funding. The author declares no conflicts of interest.
 
 # References