From 784f9044d5ce2680f493e9b81d40e210c86811df Mon Sep 17 00:00:00 2001
From: filoferra <66474175+filoferra@users.noreply.github.com>
Date: Tue, 27 Oct 2026 15:53:16 +0100
Subject: [PATCH] Pass the library to the pointwise macros

@simple_fn, @residual_fn, @jacobian_fn and @bd_fn take the library or its type first, @residual_fn(petsclib, f0, Nf), and build the C signature from its PetscInt, PetscReal and PetscScalar.

The old form, which read those three types from consts the caller had to define and so allowed one library per module, still works and warns.
---
 CHANGELOG.md              |  1 +
 docs/src/man/dmplex.md    | 13 +++---
 examples/dmplex/ex12.jl   | 26 +++++------
 examples/dmplex/ex17.jl   | 40 ++++++++---------
 examples/dmplex/ex62.jl   | 30 ++++++-------
 examples/dmplex/ex62b.jl  | 57 ++++++++++++-----------
 examples/dmplex/ex69.jl   | 36 +++++++--------
 scripts/renames.jl        | 34 ++++++++++++++
 src/deprecations.jl       |  5 +++
 src/dmplex.jl             | 95 ++++++++++++++++++++++++++++-----------
 test/core/deprecations.jl | 18 ++++++++
 test/dm/dmplex.jl         | 37 ++++++++++-----
 12 files changed, 250 insertions(+), 142 deletions(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 3c658075c..71207cd9b 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -6,6 +6,7 @@
 
 - `local_arrays` and `restore_local_arrays!` warn: use `with_local_arrays!`. An exception between the two calls left both Vecs checked out.
 - `local_coordinate_array(dm)` warns: use `with_local_coordinate_array(f, dm)`. It now returns a copy. It used to leave the coordinate vector checked out until the garbage collector freed the array, and if the DM was destroyed first, that late hand-back crashed Julia.
+- `@simple_fn`, `@residual_fn`, `@jacobian_fn` and `@bd_fn` take the library first, `@residual_fn(petsclib, f0, Nf)`, and use its `PetscInt`, `PetscReal` and `PetscScalar` for the C signature. The library value or its type both work; it is evaluated when the macro expands, so it must be a global defined before the call. Without the library, the macros take `PetscInt`, `PetscScalar` and `PetscReal` from the calling module's globals as before, and warn. That form needed those three consts and allowed only one library per module.
 
 ### Added
 
diff --git a/docs/src/man/dmplex.md b/docs/src/man/dmplex.md
index c881ab571..5b2a00471 100644
--- a/docs/src/man/dmplex.md
+++ b/docs/src/man/dmplex.md
@@ -98,8 +98,11 @@ PETSc.set_name!(petsclib, dm, "stokes")
 ## Pointwise Callbacks
 
 PETSc's FEM assembly calls user-supplied C function pointers at each
-quadrature point.  Three macros generate the required `@cfunction` wrappers
-from plain Julia functions.
+quadrature point.  Four macros (`@simple_fn`, `@residual_fn`, `@jacobian_fn`,
+`@bd_fn`) generate the required `@cfunction` wrappers from plain Julia
+functions.  Their first argument is the library (or its type), which sets the
+C types; it is evaluated when the macro expands, so it must be a global
+defined before the macro call.
 
 ### `@simple_fn` — boundary / exact-solution functions
 
@@ -111,7 +114,7 @@ function bg_vel(t, x, u, ctx)
     u[1] = x[1] * 0.1    # Vx = exx * x
     u[2] = x[2] * (-0.1) # Vy = eyy * y
 end
-const bg_vel_ptr = PETSc.@simple_fn(bg_vel)
+const bg_vel_ptr = PETSc.@simple_fn(petsclib, bg_vel)
 ```
 
 ### `@residual_fn` — f0 / f1 residual and post-processing functions
@@ -127,7 +130,7 @@ function f0_p(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
     f0[1] = 0.0
     for d in 0:dim_-1; f0[1] -= u_x[d*dim_+d+1]; end
 end
-const f0_p_ptr = PETSc.@residual_fn(f0_p, 1)
+const f0_p_ptr = PETSc.@residual_fn(petsclib, f0_p, 1)
 ```
 
 ### `@jacobian_fn` — g0/g1/g2/g3 Jacobian functions
@@ -140,7 +143,7 @@ function g1_pu(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
                aOff, aOff_x, a, a_t, a_x, t, utShift, x, nC, cst, g1)
     for d in 0:dim_-1; g1[d*dim_+d+1] = -1.0; end
 end
-const g1_pu_ptr = PETSc.@jacobian_fn(g1_pu, dim_*dim_)
+const g1_pu_ptr = PETSc.@jacobian_fn(petsclib, g1_pu, dim_*dim_)
 ```
 
 !!! note "Automatic differentiation"
diff --git a/examples/dmplex/ex12.jl b/examples/dmplex/ex12.jl
index cde7f0801..eee9e8e0b 100644
--- a/examples/dmplex/ex12.jl
+++ b/examples/dmplex/ex12.jl
@@ -130,8 +130,6 @@ MPI.Initialized() || MPI.Init()
 PETSc.initialize(petsclib)
 
 const PetscInt    = petsclib.PetscInt
-const PetscScalar = Float64
-const PetscReal   = Float64
 
 comm = MPI.COMM_WORLD
 
@@ -148,13 +146,13 @@ let _b = get(NamedTuple(pairs(opts)), :bc,    nothing); global bc    = _b === no
 function quadratic_u(t, x, u, ctx)
     u[1] = sum(xi^2 for xi in x)
 end
-const quadratic_u_ptr = PETSc.@simple_fn(quadratic_u)
+const quadratic_u_ptr = PETSc.@simple_fn(petsclib, quadratic_u)
 
 # Zero initial guess for the interior DOFs.
 function zero_u(t, x, u, ctx)
     fill!(u, 0)
 end
-const zero_u_ptr = PETSc.@simple_fn(zero_u)
+const zero_u_ptr = PETSc.@simple_fn(petsclib, zero_u)
 
 # ── COEFF_FIELD:  κ(x) = x₂ + 1 ─────────────────────────────────────────────
 #
@@ -166,14 +164,14 @@ function f0_field(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
                   aOff, aOff_x, a, a_t, a_x, t, x, numConstants, constants, f0)
     f0[1] = (2*dim_+2)*x[2] + 2*dim_
 end
-const f0_field_ptr = PETSc.@residual_fn(f0_field, Nf)
+const f0_field_ptr = PETSc.@residual_fn(petsclib, f0_field, Nf)
 
 function f1_field(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
                   aOff, aOff_x, a, a_t, a_x, t, x, numConstants, constants, f1)
     kap = x[2] + 1.0
     @inbounds for d in 1:dim_; f1[d] = kap * u_x[d]; end
 end
-const f1_field_ptr = PETSc.@residual_fn(f1_field, dim_)
+const f1_field_ptr = PETSc.@residual_fn(petsclib, f1_field, dim_)
 
 function g3_field(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
                   aOff, aOff_x, a, a_t, a_x, t, u_tShift, x, numConstants, constants, g3)
@@ -181,7 +179,7 @@ function g3_field(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
     fill!(g3, 0)
     @inbounds for d in 1:dim_; g3[(d-1)*dim_ + d] = kap; end
 end
-const g3_field_ptr = PETSc.@jacobian_fn(g3_field, dim_*dim_)
+const g3_field_ptr = PETSc.@jacobian_fn(petsclib, g3_field, dim_*dim_)
 
 # Neumann flux for COEFF_FIELD:  h = κ(x)·∇u·n = (x₂+1)·2(x·n)
 # Convention: f0_bd = −h  (PETSc boundary residual sign)
@@ -192,7 +190,7 @@ function f0_bd_field(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
     @inbounds for d in 1:dim_; xdotn += x[d]*n[d]; end
     f0[1] = -kap * 2.0 * xdotn
 end
-const f0_bd_field_ptr = PETSc.@bd_fn(f0_bd_field, Nf)
+const f0_bd_field_ptr = PETSc.@bd_fn(petsclib, f0_bd_field, Nf)
 
 # ── COEFF_NONLINEAR:  κ(u) = u + 1 ──────────────────────────────────────────
 #
@@ -214,28 +212,28 @@ function f0_nl(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
     u_val = u[Int(uOff[1])+1]
     f0[1] = (2*dim_+4)*u_val + 2*dim_
 end
-const f0_nl_ptr = PETSc.@residual_fn(f0_nl, Nf)
+const f0_nl_ptr = PETSc.@residual_fn(petsclib, f0_nl, Nf)
 
 function f1_nl(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
                aOff, aOff_x, a, a_t, a_x, t, x, numConstants, constants, f1)
     kap = u[Int(uOff[1])+1] + 1.0
     @inbounds for d in 1:dim_; f1[d] = kap * u_x[d]; end
 end
-const f1_nl_ptr = PETSc.@residual_fn(f1_nl, dim_)
+const f1_nl_ptr = PETSc.@residual_fn(petsclib, f1_nl, dim_)
 
 # g0: ∂f0/∂u = 2·dim+4
 function g0_nl(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
                aOff, aOff_x, a, a_t, a_x, t, u_tShift, x, numConstants, constants, g0)
     g0[1] = 2*dim_ + 4
 end
-const g0_nl_ptr = PETSc.@jacobian_fn(g0_nl, Nf*Nf)
+const g0_nl_ptr = PETSc.@jacobian_fn(petsclib, g0_nl, Nf*Nf)
 
 # g2: ∂f1[d]/∂u = u_x[d]   (κ′(u)=1, so ∂((u+1)·u_x[d])/∂u = u_x[d])
 function g2_nl(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
                aOff, aOff_x, a, a_t, a_x, t, u_tShift, x, numConstants, constants, g2)
     @inbounds for d in 1:dim_; g2[d] = u_x[d]; end
 end
-const g2_nl_ptr = PETSc.@jacobian_fn(g2_nl, Nf*Nf*dim_)
+const g2_nl_ptr = PETSc.@jacobian_fn(petsclib, g2_nl, Nf*Nf*dim_)
 
 # g3: ∂f1[d]/∂(∂u/∂xd′) = (u+1)·δ_{dd′}
 function g3_nl(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
@@ -244,7 +242,7 @@ function g3_nl(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
     fill!(g3, 0)
     @inbounds for d in 1:dim_; g3[(d-1)*dim_ + d] = kap; end
 end
-const g3_nl_ptr = PETSc.@jacobian_fn(g3_nl, dim_*dim_)
+const g3_nl_ptr = PETSc.@jacobian_fn(petsclib, g3_nl, dim_*dim_)
 
 # Neumann flux for COEFF_NONLINEAR:  h = (u+1)·2(x·n)
 # f0_bd = −h
@@ -255,7 +253,7 @@ function f0_bd_nl(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
     @inbounds for d in 1:dim_; xdotn += x[d]*n[d]; end
     f0[1] = -kap * 2.0 * xdotn
 end
-const f0_bd_nl_ptr = PETSc.@bd_fn(f0_bd_nl, Nf)
+const f0_bd_nl_ptr = PETSc.@bd_fn(petsclib, f0_bd_nl, Nf)
 
 # ── Build the mesh ────────────────────────────────────────────────────────────
 # 3-D simplex (TetGen) is broken in the current PETSc_jll; use hexahedra for 3D.
diff --git a/examples/dmplex/ex17.jl b/examples/dmplex/ex17.jl
index 946ce7508..8868242a9 100644
--- a/examples/dmplex/ex17.jl
+++ b/examples/dmplex/ex17.jl
@@ -123,8 +123,6 @@ MPI.Initialized() || MPI.Init()
 PETSc.initialize(petsclib)
 
 const PetscInt    = petsclib.PetscInt
-const PetscScalar = Float64
-const PetscReal   = Float64
 
 comm = MPI.COMM_WORLD
 
@@ -154,39 +152,39 @@ VALID_SOL = ("vlap_quad", "elas_quad", "vlap_trig", "elas_trig",
 function zero_u(t, x, u, ctx)
     fill!(u, 0)
 end
-const zero_u_ptr = PETSc.@simple_fn(zero_u)
+const zero_u_ptr = PETSc.@simple_fn(petsclib, zero_u)
 
 function ge_shift_u(t, x, u, ctx)
     u[1] = 0.1
     for d in 2:length(u); u[d] = 0.0; end
 end
-const ge_shift_u_ptr = PETSc.@simple_fn(ge_shift_u)
+const ge_shift_u_ptr = PETSc.@simple_fn(petsclib, ge_shift_u)
 
 function quadratic_2d_u(t, x, u, ctx)
     u[1] = x[1]^2
     u[2] = x[2]^2 - 2*x[1]*x[2]
 end
-const quadratic_2d_u_ptr = PETSc.@simple_fn(quadratic_2d_u)
+const quadratic_2d_u_ptr = PETSc.@simple_fn(petsclib, quadratic_2d_u)
 
 function quadratic_3d_u(t, x, u, ctx)
     u[1] = x[1]^2
     u[2] = x[2]^2 - 2*x[1]*x[2]
     u[3] = x[3]^2 - 2*x[2]*x[3]
 end
-const quadratic_3d_u_ptr = PETSc.@simple_fn(quadratic_3d_u)
+const quadratic_3d_u_ptr = PETSc.@simple_fn(petsclib, quadratic_3d_u)
 
 function trig_2d_u(t, x, u, ctx)
     u[1] = sin(2π * x[1])
     u[2] = sin(2π * x[2]) - 2*x[1]*x[2]
 end
-const trig_2d_u_ptr = PETSc.@simple_fn(trig_2d_u)
+const trig_2d_u_ptr = PETSc.@simple_fn(petsclib, trig_2d_u)
 
 function trig_3d_u(t, x, u, ctx)
     u[1] = sin(2π * x[1])
     u[2] = sin(2π * x[2]) - 2*x[1]*x[2]
     u[3] = sin(2π * x[3]) - 2*x[2]*x[3]
 end
-const trig_3d_u_ptr = PETSc.@simple_fn(trig_3d_u)
+const trig_3d_u_ptr = PETSc.@simple_fn(petsclib, trig_3d_u)
 
 function axial_disp_u(t, x, u, ctx)
     mu = mu_val; lambda = lambda_val; N = N_val
@@ -195,7 +193,7 @@ function axial_disp_u(t, x, u, ctx)
     u[2] = 0.25*lambda / (mu*(lambda+mu)) * N * x[2]
     for d in 3:length(u); u[d] = 0.0; end
 end
-const axial_disp_u_ptr = PETSc.@simple_fn(axial_disp_u)
+const axial_disp_u_ptr = PETSc.@simple_fn(petsclib, axial_disp_u)
 
 function uniform_strain_u(t, x, u, ctx)
     eps_xx = 0.1; eps_xy = 0.3; eps_yy = 0.25
@@ -203,7 +201,7 @@ function uniform_strain_u(t, x, u, ctx)
     u[2] = eps_xy * x[1] + eps_yy * x[2]
     for d in 3:length(u); u[d] = 0.0; end
 end
-const uniform_strain_u_ptr = PETSc.@simple_fn(uniform_strain_u)
+const uniform_strain_u_ptr = PETSc.@simple_fn(petsclib, uniform_strain_u)
 
 # ── Residual functions f0 (zeroth-order) ───────────────────────────────────────
 # f0 output size = Nc = dim (one entry per displacement component).
@@ -214,7 +212,7 @@ function f0_vlap_quadratic(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
                            aOff, aOff_x, a, a_t, a_x, t, x, numConstants, constants, f0)
     for d in 1:dim_; f0[d] += 2.0; end
 end
-const f0_vlap_quadratic_ptr = PETSc.@residual_fn(f0_vlap_quadratic, dim_)
+const f0_vlap_quadratic_ptr = PETSc.@residual_fn(petsclib, f0_vlap_quadratic, dim_)
 
 # Linear elasticity, quadratic solution: body force from −∇·σ.
 function f0_elas_quadratic(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
@@ -223,14 +221,14 @@ function f0_elas_quadratic(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
     for d in 1:dim_-1; f0[d] += 2.0*mu; end
     f0[dim_] += 2.0*lambda + 4.0*mu
 end
-const f0_elas_quadratic_ptr = PETSc.@residual_fn(f0_elas_quadratic, dim_)
+const f0_elas_quadratic_ptr = PETSc.@residual_fn(petsclib, f0_elas_quadratic, dim_)
 
 # Vector Laplacian, trig solution: body force = 4π² sin(2πx_d) per component.
 function f0_vlap_trig(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
                       aOff, aOff_x, a, a_t, a_x, t, x, numConstants, constants, f0)
     for d in 1:dim_; f0[d] += -4.0*π^2 * sin(2π * x[d]); end
 end
-const f0_vlap_trig_ptr = PETSc.@residual_fn(f0_vlap_trig, dim_)
+const f0_vlap_trig_ptr = PETSc.@residual_fn(petsclib, f0_vlap_trig, dim_)
 
 # Linear elasticity, trig solution: body force from −∇·σ.
 function f0_elas_trig(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
@@ -242,7 +240,7 @@ function f0_elas_trig(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
                  (d < dim_ ? 2.0*(mu + lambda) : 0.0)
     end
 end
-const f0_elas_trig_ptr = PETSc.@residual_fn(f0_elas_trig, dim_)
+const f0_elas_trig_ptr = PETSc.@residual_fn(petsclib, f0_elas_trig, dim_)
 
 # ── Residual functions f1 (first-order, diffusion flux) ────────────────────────
 # f1 output size = Nc * dim = dim * dim (gradient of Nc-component field).
@@ -254,7 +252,7 @@ function f1_vlap(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
         f1[c*dim_ + d + 1] += u_x[c*dim_ + d + 1]
     end
 end
-const f1_vlap_ptr = PETSc.@residual_fn(f1_vlap, dim_*dim_)
+const f1_vlap_ptr = PETSc.@residual_fn(petsclib, f1_vlap, dim_*dim_)
 
 function f1_elas(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
                  aOff, aOff_x, a, a_t, a_x, t, x, numConstants, constants, f1)
@@ -264,7 +262,7 @@ function f1_elas(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
         f1[c*dim_ + c + 1] += lambda * u_x[d*dim_ + d + 1]
     end
 end
-const f1_elas_ptr = PETSc.@residual_fn(f1_elas, dim_*dim_)
+const f1_elas_ptr = PETSc.@residual_fn(petsclib, f1_elas, dim_*dim_)
 
 # ── Jacobian functions ─────────────────────────────────────────────────────────
 # g3 output size = Nc*Nc*dim*dim = dim^4 (full stiffness tensor).
@@ -277,7 +275,7 @@ function g3_vlap(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
         g3[((c*dim_ + c)*dim_ + d)*dim_ + d + 1] = 1.0
     end
 end
-const g3_vlap_ptr = PETSc.@jacobian_fn(g3_vlap, dim_*dim_*dim_*dim_)
+const g3_vlap_ptr = PETSc.@jacobian_fn(petsclib, g3_vlap, dim_*dim_*dim_*dim_)
 
 function g3_elas(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
                  aOff, aOff_x, a, a_t, a_x, t, u_tShift, x, numConstants, constants, g3)
@@ -289,7 +287,7 @@ function g3_elas(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
         g3[((c*dim_ + d)*dim_ + c)*dim_ + d + 1] += lambda
     end
 end
-const g3_elas_ptr = PETSc.@jacobian_fn(g3_elas, dim_*dim_*dim_*dim_)
+const g3_elas_ptr = PETSc.@jacobian_fn(petsclib, g3_elas, dim_*dim_*dim_*dim_)
 
 # ── VTK post-processing: displacement (3-vector) and stress (3×3 tensor) ──────
 #
@@ -307,7 +305,7 @@ function copy_displacement_3(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
     for c in 1:dim_; out[c] = u[c]; end
     for c in dim_+1:3; out[c] = 0.0; end
 end
-const copy_displacement_3_ptr = PETSc.@residual_fn(copy_displacement_3, 3)
+const copy_displacement_3_ptr = PETSc.@residual_fn(petsclib, copy_displacement_3, 3)
 
 # Compute Cauchy stress σ_{cd} = μ(∂u_c/∂x_d + ∂u_d/∂x_c) + λ tr(ε) δ_{cd}
 # and place it in a 3×3 output (9 components, row-major).
@@ -323,7 +321,7 @@ function compute_stress_3x3(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
     end
     for c in 0:dim_-1; out[c*3+c+1] += lambda * div_u; end
 end
-const compute_stress_3x3_ptr = PETSc.@residual_fn(compute_stress_3x3, 9)
+const compute_stress_3x3_ptr = PETSc.@residual_fn(petsclib, compute_stress_3x3, 9)
 
 # ── Neumann BC for elas_axial_disp: traction N on right wall ──────────────────
 # f0_bd output size = Nc = dim.
@@ -332,7 +330,7 @@ function f0_elas_axial_disp_bd(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
                                 numConstants, constants, f0)
     f0[1] = real(constants[3])   # N — tension force, constants[3] = N (1-indexed)
 end
-const f0_elas_axial_disp_bd_ptr = PETSc.@bd_fn(f0_elas_axial_disp_bd, dim_)
+const f0_elas_axial_disp_bd_ptr = PETSc.@bd_fn(petsclib, f0_elas_axial_disp_bd, dim_)
 
 # ── Near-null space constructor for AMG ────────────────────────────────────────
 # Called by PETSc for each level when setting up the preconditioner.
diff --git a/examples/dmplex/ex62.jl b/examples/dmplex/ex62.jl
index bb4c5ef64..ebd658f87 100644
--- a/examples/dmplex/ex62.jl
+++ b/examples/dmplex/ex62.jl
@@ -104,8 +104,6 @@ MPI.Initialized() || MPI.Init()
 PETSc.initialize(petsclib)
 
 const PetscInt    = petsclib.PetscInt
-const PetscScalar = Float64
-const PetscReal   = Float64
 
 comm = MPI.COMM_WORLD
 
@@ -135,7 +133,7 @@ function f1_u(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
         f1[c*dim_+c+1] -= u[dim_+1]   # subtract pressure (field 1)
     end
 end
-const f1_u_ptr = PETSc.@residual_fn(f1_u, dim_*dim_)
+const f1_u_ptr = PETSc.@residual_fn(petsclib, f1_u, dim_*dim_)
 
 # Continuity: f0 = −∇·u
 function f0_p(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
@@ -145,7 +143,7 @@ function f0_p(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
         f0[1] -= u_x[d*dim_+d+1]
     end
 end
-const f0_p_ptr = PETSc.@residual_fn(f0_p, 1)
+const f0_p_ptr = PETSc.@residual_fn(petsclib, f0_p, 1)
 
 # ── Body-force functions (one per MMS solution type) ──────────────────────────
 #
@@ -163,7 +161,7 @@ function f0_quadratic_u(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
         f0[d] = 4mu - 1.0
     end
 end
-const f0_quadratic_u_ptr = PETSc.@residual_fn(f0_quadratic_u, dim_)
+const f0_quadratic_u_ptr = PETSc.@residual_fn(petsclib, f0_quadratic_u, dim_)
 
 # Trigonometric MMS body force:
 #   f_x = 2π cos(2πx) + μ(dim−1)π² sin(πx) + μπ² Σ_{d>0} sin(πx_d)
@@ -178,7 +176,7 @@ function f0_trig_u(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
         f0[d]  = -2π*cos(2π*x[d]) + mu*π^3*cos(π*x[1])*x[d]
     end
 end
-const f0_trig_u_ptr = PETSc.@residual_fn(f0_trig_u, dim_)
+const f0_trig_u_ptr = PETSc.@residual_fn(petsclib, f0_trig_u, dim_)
 
 # ── Jacobians ─────────────────────────────────────────────────────────────────
 
@@ -191,28 +189,28 @@ function g3_uu(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
         g3[((c*dim_+d)*dim_+d)*dim_+c+1] += mu
     end
 end
-const g3_uu_ptr = PETSc.@jacobian_fn(g3_uu, dim_*dim_*dim_*dim_)
+const g3_uu_ptr = PETSc.@jacobian_fn(petsclib, g3_uu, dim_*dim_*dim_*dim_)
 
 # J_up (g2): −⟨∇·v, p⟩
 function g2_up(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
                aOff, aOff_x, a, a_t, a_x, t, utShift, x, nC, cst, g2)
     for d in 0:dim_-1; g2[d*dim_+d+1] = -1.0; end
 end
-const g2_up_ptr = PETSc.@jacobian_fn(g2_up, dim_*dim_)
+const g2_up_ptr = PETSc.@jacobian_fn(petsclib, g2_up, dim_*dim_)
 
 # J_pu (g1): ⟨q, −∇·u⟩
 function g1_pu(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
                aOff, aOff_x, a, a_t, a_x, t, utShift, x, nC, cst, g1)
     for d in 0:dim_-1; g1[d*dim_+d+1] = -1.0; end
 end
-const g1_pu_ptr = PETSc.@jacobian_fn(g1_pu, dim_*dim_)
+const g1_pu_ptr = PETSc.@jacobian_fn(petsclib, g1_pu, dim_*dim_)
 
 # J_pp preconditioner (g0): (1/μ)⟨q, p⟩
 function g0_pp(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
                aOff, aOff_x, a, a_t, a_x, t, utShift, x, nC, cst, g0)
     g0[1] = 1.0 / real(cst[1])
 end
-const g0_pp_ptr = PETSc.@jacobian_fn(g0_pp, 1)
+const g0_pp_ptr = PETSc.@jacobian_fn(petsclib, g0_pp, 1)
 
 # ── Exact solutions ───────────────────────────────────────────────────────────
 #
@@ -228,12 +226,12 @@ function quadratic_vel(t, x, u, ctx)
         u[c]  = 2*x[1]^2 - 2*x[1]*x[c]
     end
 end
-const quadratic_vel_ptr = PETSc.@simple_fn(quadratic_vel)
+const quadratic_vel_ptr = PETSc.@simple_fn(petsclib, quadratic_vel)
 
 function quadratic_pres(t, x, u, ctx)
     u[1] = sum(x) - 0.5*length(x)
 end
-const quadratic_pres_ptr = PETSc.@simple_fn(quadratic_pres)
+const quadratic_pres_ptr = PETSc.@simple_fn(petsclib, quadratic_pres)
 
 #
 # Trigonometric MMS:
@@ -249,23 +247,23 @@ function trig_vel(t, x, u, ctx)
         u[c]  = -π*cos(π*x[1])*x[c]
     end
 end
-const trig_vel_ptr = PETSc.@simple_fn(trig_vel)
+const trig_vel_ptr = PETSc.@simple_fn(petsclib, trig_vel)
 
 function trig_pres(t, x, u, ctx)
     u[1] = sum(sin(2π*xi) for xi in x)
 end
-const trig_pres_ptr = PETSc.@simple_fn(trig_pres)
+const trig_pres_ptr = PETSc.@simple_fn(petsclib, trig_pres)
 
 # Helper functions for the pressure null space constructor
 function zero_vel(t, x, u, ctx)
     for c in 1:length(u); u[c] = 0.0; end
 end
-const zero_vel_ptr = PETSc.@simple_fn(zero_vel)
+const zero_vel_ptr = PETSc.@simple_fn(petsclib, zero_vel)
 
 function one_pres(t, x, u, ctx)
     u[1] = 1.0
 end
-const one_pres_ptr = PETSc.@simple_fn(one_pres)
+const one_pres_ptr = PETSc.@simple_fn(petsclib, one_pres)
 
 # ── Select exact solution functions based on sol_type ─────────────────────────
 exact_vel_ptr  = sol_type == "quadratic" ? quadratic_vel_ptr  : trig_vel_ptr
diff --git a/examples/dmplex/ex62b.jl b/examples/dmplex/ex62b.jl
index 1a4a8f199..28007dd3f 100644
--- a/examples/dmplex/ex62b.jl
+++ b/examples/dmplex/ex62b.jl
@@ -350,7 +350,6 @@ PETSc.initialize(petsclib)
 
 @isdefined(PetscInt)    || const PetscInt    = petsclib.PetscInt
 @isdefined(PetscScalar) || const PetscScalar = Float64
-@isdefined(PetscReal)   || const PetscReal   = Float64
 
 # Module-level globals for background strain rates — read by @cfunction-based BCs.
 bg_exx = 0.0
@@ -518,7 +517,7 @@ function f1_u(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
     compute_strainrate!(f1, dim_, u_x)
     deviatoric_stress_pressure!(f1, dim_, f1, phase, cst, p, tau_old)
 end
-const f1_u_ptr = PETSc.@residual_fn(f1_u, dim_*dim_)
+const f1_u_ptr = PETSc.@residual_fn(petsclib, f1_u, dim_*dim_)
 
 # Continuity: f0 = −∇·u
 function f0_p(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
@@ -528,7 +527,7 @@ function f0_p(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
         f0[1] -= u_x[d*dim_+d+1]
     end
 end
-const f0_p_ptr = PETSc.@residual_fn(f0_p, 1)
+const f0_p_ptr = PETSc.@residual_fn(petsclib, f0_p, 1)
 
 # ── Body-force functions (one per MMS solution type) ──────────────────────────
 #
@@ -546,7 +545,7 @@ function f0_quadratic_u(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
         f0[d] = 4mu - 1.0
     end
 end
-const f0_quadratic_u_ptr = PETSc.@residual_fn(f0_quadratic_u, dim_)
+const f0_quadratic_u_ptr = PETSc.@residual_fn(petsclib, f0_quadratic_u, dim_)
 
 # Trigonometric MMS body force:
 #   f_x = 2π cos(2πx) + μ(dim−1)π² sin(πx) + μπ² Σ_{d>0} sin(πx_d)
@@ -561,14 +560,14 @@ function f0_trig_u(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
         f0[d]  = -2π*cos(2π*x[d]) + mu*π^3*cos(π*x[1])*x[d]
     end
 end
-const f0_trig_u_ptr = PETSc.@residual_fn(f0_trig_u, dim_)
+const f0_trig_u_ptr = PETSc.@residual_fn(petsclib, f0_trig_u, dim_)
 
 # Zero body force for background deformation mode (no MMS forcing)
 function f0_zero_u(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
                    aOff, aOff_x, a, a_t, a_x, t, x, nC, cst, f0)
     for d in 1:dim_; f0[d] = 0.0; end
 end
-const f0_zero_u_ptr = PETSc.@residual_fn(f0_zero_u, dim_)
+const f0_zero_u_ptr = PETSc.@residual_fn(petsclib, f0_zero_u, dim_)
 
 # Gravity body force: f_body = ρ g (−ŷ direction).
 # In PETSc weak-form sign convention, f0 = −f_body, so:
@@ -581,7 +580,7 @@ function f0_grav_u(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
     for d in 1:dim_; f0[d] = 0.0; end
     f0[dim_] = rho * g
 end
-const f0_grav_u_ptr = PETSc.@residual_fn(f0_grav_u, dim_)
+const f0_grav_u_ptr = PETSc.@residual_fn(petsclib, f0_grav_u, dim_)
 
 # ── Jacobians ─────────────────────────────────────────────────────────────────
 
@@ -593,7 +592,7 @@ function g3_uu(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
     p         = real(u[dim_+1])
     deviatoric_stress_tangent!(g3, dim_, u_x, phase, cst, p, tau_old)
 end
-const g3_uu_ptr = PETSc.@jacobian_fn(g3_uu, dim_*dim_*dim_*dim_)
+const g3_uu_ptr = PETSc.@jacobian_fn(petsclib, g3_uu, dim_*dim_*dim_*dim_)
 
 # J_up (g2): ∂σ_{fc,df}/∂p via AD of deviatoric_stress_pressure! w.r.t. p.
 # Layout: g2[fc*dim+df + 1] = ∂σ_{fc,df}/∂p.
@@ -615,14 +614,14 @@ function g2_up(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
     end
     for i in 1:n; g2[i] = dσ_dp[i]; end
 end
-const g2_up_ptr = PETSc.@jacobian_fn(g2_up, dim_*dim_)
+const g2_up_ptr = PETSc.@jacobian_fn(petsclib, g2_up, dim_*dim_)
 
 # J_pu (g1): ⟨q, −∇·u⟩ — kinematic incompressibility, independent of rheology.
 function g1_pu(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
                aOff, aOff_x, a, a_t, a_x, t, utShift, x, nC, cst, g1)
     for d in 0:dim_-1; g1[d*dim_+d+1] = -1.0; end
 end
-const g1_pu_ptr = PETSc.@jacobian_fn(g1_pu, dim_*dim_)
+const g1_pu_ptr = PETSc.@jacobian_fn(petsclib, g1_pu, dim_*dim_)
 
 # J_pp preconditioner (g0): (1/η_eff)⟨q, p⟩ using the Maxwell effective viscosity
 function g0_pp(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
@@ -640,7 +639,7 @@ function g0_pp(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
     eta_eff = eta * G_mod * dt / (eta + G_mod * dt)
     g0[1]   = 1.0 / eta_eff
 end
-const g0_pp_ptr = PETSc.@jacobian_fn(g0_pp, 1)
+const g0_pp_ptr = PETSc.@jacobian_fn(petsclib, g0_pp, 1)
 
 # ── Deviatoric stress computation (for VTK and viscoelastic history) ─────────
 #
@@ -668,7 +667,7 @@ function compute_tau(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
         out[i] = 2 * eta_eff * eps_buf[i] + relax_fac * real(tau_old[i])
     end
 end
-const compute_tau_ptr = PETSc.@residual_fn(compute_tau, dim_*dim_)
+const compute_tau_ptr = PETSc.@residual_fn(petsclib, compute_tau, dim_*dim_)
 
 # Same as compute_tau but zero-padded to 3×3 row-major for VTK tensor output.
 function compute_tau_3x3(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
@@ -691,7 +690,7 @@ function compute_tau_3x3(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
         out[c*3+d+1] = 2 * eta_eff * eps_buf[c*dim_+d+1] + relax_fac * real(tau_old[c*dim_+d+1])
     end
 end
-const compute_tau_3x3_ptr = PETSc.@residual_fn(compute_tau_3x3, 9)
+const compute_tau_3x3_ptr = PETSc.@residual_fn(petsclib, compute_tau_3x3, 9)
 
 # Second invariants (DG-P0, 1 DOF/cell) — used for domain averages and VTK output.
 # ε_II = sqrt(½ ε:ε),  τ_II = sqrt(½ τ:τ)
@@ -703,7 +702,7 @@ function compute_eps_II(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
     compute_strainrate!(eps_buf, dim_, eps_buf)
     out[1] = real(second_invariant(eps_buf, dim_))
 end
-const compute_eps_II_ptr = PETSc.@residual_fn(compute_eps_II, 1)
+const compute_eps_II_ptr = PETSc.@residual_fn(petsclib, compute_eps_II, 1)
 
 function compute_tau_II(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
                         aOff, aOff_x, a, a_t, a_x, t, x, nC, cst, out)
@@ -727,7 +726,7 @@ function compute_tau_II(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
     end
     out[1] = sqrt(0.5 * s)
 end
-const compute_tau_II_ptr = PETSc.@residual_fn(compute_tau_II, 1)
+const compute_tau_II_ptr = PETSc.@residual_fn(petsclib, compute_tau_II, 1)
 
 # ── VTK post-processing: strain rate tensor ─────────────────────────────────
 #
@@ -741,7 +740,7 @@ function compute_strainrate_3x3(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
         out[c*3+d+1] = 0.5*(u_x[c*dim_+d+1] + u_x[d*dim_+c+1])
     end
 end
-const compute_strainrate_3x3_ptr = PETSc.@residual_fn(compute_strainrate_3x3, 9)
+const compute_strainrate_3x3_ptr = PETSc.@residual_fn(petsclib, compute_strainrate_3x3, 9)
 
 # VTK copy kernels: extract velocity (field 0) and pressure (field 1) from the
 # combined solution vector.  Must be at top-level scope for @cfunction.
@@ -749,13 +748,13 @@ function copy_vel(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
                   aOff, aOff_x, a, a_t, a_x, t, x, nC, cst, out_)
     for c in 1:dim_; out_[c] = u[c]; end
 end
-const copy_vel_ptr = PETSc.@residual_fn(copy_vel, dim_)
+const copy_vel_ptr = PETSc.@residual_fn(petsclib, copy_vel, dim_)
 
 function copy_pres(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
                    aOff, aOff_x, a, a_t, a_x, t, x, nC, cst, out_)
     out_[1] = u[dim_+1]
 end
-const copy_pres_ptr = PETSc.@residual_fn(copy_pres, 1)
+const copy_pres_ptr = PETSc.@residual_fn(petsclib, copy_pres, 1)
 
 # copy_mu reports the Maxwell effective viscosity η_eff = η·G·dt/(η+G·dt).
 function copy_mu(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
@@ -772,13 +771,13 @@ function copy_mu(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
     dt      = real(cst[17])
     out_[1] = eta * G_mod * dt / (eta + G_mod * dt)
 end
-const copy_mu_ptr = PETSc.@residual_fn(copy_mu, 1)
+const copy_mu_ptr = PETSc.@residual_fn(petsclib, copy_mu, 1)
 
 function copy_phase(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
                     aOff, aOff_x, a, a_t, a_x, t, x, nC, cst, out_)
     out_[1] = real(a[1])   # phase from aux field — correct after Lagrangian advection
 end
-const copy_phase_ptr = PETSc.@residual_fn(copy_phase, 1)
+const copy_phase_ptr = PETSc.@residual_fn(petsclib, copy_phase, 1)
 
 # Pass-through kernel: copies dim×dim values from the source field to the output.
 # Used by advance_tau! to project tau_vec (on dm_tau) into field 1 of aux_vec.
@@ -786,7 +785,7 @@ function identity_tau(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
                       aOff, aOff_x, a, a_t, a_x, t, x, nC, cst, out)
     for i in 1:dim_*dim_; out[i] = u[i]; end
 end
-const identity_tau_ptr = PETSc.@residual_fn(identity_tau, dim_*dim_)
+const identity_tau_ptr = PETSc.@residual_fn(petsclib, identity_tau, dim_*dim_)
 
 # ── Time-stepping helpers ─────────────────────────────────────────────────────
 
@@ -986,12 +985,12 @@ function quadratic_vel(t, x, u, ctx)
         u[c]  = 2*x[1]^2 - 2*x[1]*x[c]
     end
 end
-const quadratic_vel_ptr = PETSc.@simple_fn(quadratic_vel)
+const quadratic_vel_ptr = PETSc.@simple_fn(petsclib, quadratic_vel)
 
 function quadratic_pres(t, x, u, ctx)
     u[1] = sum(x) - 0.5*length(x)
 end
-const quadratic_pres_ptr = PETSc.@simple_fn(quadratic_pres)
+const quadratic_pres_ptr = PETSc.@simple_fn(petsclib, quadratic_pres)
 
 #
 # Trigonometric MMS:
@@ -1007,23 +1006,23 @@ function trig_vel(t, x, u, ctx)
         u[c]  = -π*cos(π*x[1])*x[c]
     end
 end
-const trig_vel_ptr = PETSc.@simple_fn(trig_vel)
+const trig_vel_ptr = PETSc.@simple_fn(petsclib, trig_vel)
 
 function trig_pres(t, x, u, ctx)
     u[1] = sum(sin(2π*xi) for xi in x)
 end
-const trig_pres_ptr = PETSc.@simple_fn(trig_pres)
+const trig_pres_ptr = PETSc.@simple_fn(petsclib, trig_pres)
 
 # Helper functions for the pressure null space constructor
 function zero_vel(t, x, u, ctx)
     for c in 1:length(u); u[c] = 0.0; end
 end
-const zero_vel_ptr = PETSc.@simple_fn(zero_vel)
+const zero_vel_ptr = PETSc.@simple_fn(petsclib, zero_vel)
 
 function one_pres(t, x, u, ctx)
     u[1] = 1.0
 end
-const one_pres_ptr = PETSc.@simple_fn(one_pres)
+const one_pres_ptr = PETSc.@simple_fn(petsclib, one_pres)
 
 # Background deformation BC: Vx = x*exx_bg, Vy = y*eyy_bg, Vz = z*ezz_bg.
 # ezz_bg defaults to -(exx_bg+eyy_bg) to satisfy incompressibility.
@@ -1032,7 +1031,7 @@ function bg_vel(t, x, u, ctx)
     u[2] = x[2] * bg_eyy
     length(u) >= 3 && (u[3] = x[3] * bg_ezz)
 end
-const bg_vel_ptr = PETSc.@simple_fn(bg_vel)
+const bg_vel_ptr = PETSc.@simple_fn(petsclib, bg_vel)
 
 """
     update_tau!(petsclib, dm_tau, u, tau_vec, t=0.0)
@@ -1055,7 +1054,7 @@ end
 function phase_fn(t, x, u, ctx)
     u[1] = sum(xi^2 for xi in x) < inclusion_R^2 ? 2.0 : 1.0
 end
-const phase_fn_ptr = PETSc.@simple_fn(phase_fn)
+const phase_fn_ptr = PETSc.@simple_fn(petsclib, phase_fn)
 
 exact_vel_ptr  = sol_type == "quadratic" ? quadratic_vel_ptr  :
                  sol_type == "trig"      ? trig_vel_ptr        :
diff --git a/examples/dmplex/ex69.jl b/examples/dmplex/ex69.jl
index fa4168c73..725e1f29e 100644
--- a/examples/dmplex/ex69.jl
+++ b/examples/dmplex/ex69.jl
@@ -144,8 +144,6 @@ MPI.Initialized() || MPI.Init()
 PETSc.initialize(petsclib)
 
 const PetscInt    = petsclib.PetscInt
-const PetscScalar = Float64
-const PetscReal   = Float64
 
 comm = MPI.COMM_WORLD
 
@@ -176,7 +174,7 @@ function f0_stokes_u(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
     f0[1] = 0.0
     f0[2] = -sin(real(cst[2]) * pi * x[2]) * cos(real(cst[1]) * pi * x[1])
 end
-const f0_stokes_u_ptr = PETSc.@residual_fn(f0_stokes_u, dim_)
+const f0_stokes_u_ptr = PETSc.@residual_fn(petsclib, f0_stokes_u, dim_)
 
 # Momentum flux σ = 2μ ε(u) − pI  for SolKx (μ = exp(2Bx), cst[3]=B)
 function stokes_momentum_kx(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
@@ -189,7 +187,7 @@ function stokes_momentum_kx(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
         f1[c*dim_+c+1] -= u[dim_+1]      # subtract pressure
     end
 end
-const stokes_momentum_kx_ptr = PETSc.@residual_fn(stokes_momentum_kx, dim_*dim_)
+const stokes_momentum_kx_ptr = PETSc.@residual_fn(petsclib, stokes_momentum_kx, dim_*dim_)
 
 # Momentum flux for SolCx (μ = etaA or etaB, cst[3]=etaA, cst[4]=etaB, cst[5]=xc)
 function stokes_momentum_cx(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
@@ -202,7 +200,7 @@ function stokes_momentum_cx(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
         f1[c*dim_+c+1] -= u[dim_+1]
     end
 end
-const stokes_momentum_cx_ptr = PETSc.@residual_fn(stokes_momentum_cx, dim_*dim_)
+const stokes_momentum_cx_ptr = PETSc.@residual_fn(petsclib, stokes_momentum_cx, dim_*dim_)
 
 # Continuity: f0 = −∇·u
 function stokes_mass(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
@@ -212,14 +210,14 @@ function stokes_mass(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
         f0[1] -= u_x[d*dim_+d+1]
     end
 end
-const stokes_mass_ptr = PETSc.@residual_fn(stokes_mass, 1)
+const stokes_mass_ptr = PETSc.@residual_fn(petsclib, stokes_mass, 1)
 
 # Zero f1 for pressure field
 function f1_zero_p(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
                    aOff, aOff_x, a, a_t, a_x, t, x, nC, cst, f1)
     for d in 1:dim_; f1[d] = 0.0; end
 end
-const f1_zero_p_ptr = PETSc.@residual_fn(f1_zero_p, dim_)
+const f1_zero_p_ptr = PETSc.@residual_fn(petsclib, f1_zero_p, dim_)
 
 # ── Pointwise functions: Jacobians ────────────────────────────────────────────
 
@@ -232,7 +230,7 @@ function stokes_vel_J_kx(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
         g3[((cI*dim_+d)*dim_+d)*dim_+cI+1] += mu
     end
 end
-const stokes_vel_J_kx_ptr = PETSc.@jacobian_fn(stokes_vel_J_kx, dim_*dim_*dim_*dim_)
+const stokes_vel_J_kx_ptr = PETSc.@jacobian_fn(petsclib, stokes_vel_J_kx, dim_*dim_*dim_*dim_)
 
 function stokes_vel_J_cx(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
                           aOff, aOff_x, a, a_t, a_x, t, utShift, x, nC, cst, g3)
@@ -242,35 +240,35 @@ function stokes_vel_J_cx(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
         g3[((cI*dim_+d)*dim_+d)*dim_+cI+1] += mu
     end
 end
-const stokes_vel_J_cx_ptr = PETSc.@jacobian_fn(stokes_vel_J_cx, dim_*dim_*dim_*dim_)
+const stokes_vel_J_cx_ptr = PETSc.@jacobian_fn(petsclib, stokes_vel_J_cx, dim_*dim_*dim_*dim_)
 
 # J_up (g2): −⟨∇·v, p⟩  →  g2[d*dim+d] = -1
 function stokes_pres_J(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
                         aOff, aOff_x, a, a_t, a_x, t, utShift, x, nC, cst, g2)
     for d in 0:dim_-1; g2[d*dim_+d+1] = -1.0; end
 end
-const stokes_pres_J_ptr = PETSc.@jacobian_fn(stokes_pres_J, dim_*dim_)
+const stokes_pres_J_ptr = PETSc.@jacobian_fn(petsclib, stokes_pres_J, dim_*dim_)
 
 # J_pu (g1): ⟨q, ∇·u⟩  →  g1[d*dim+d] = -1
 function stokes_mass_J(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
                         aOff, aOff_x, a, a_t, a_x, t, utShift, x, nC, cst, g1)
     for d in 0:dim_-1; g1[d*dim_+d+1] = -1.0; end
 end
-const stokes_mass_J_ptr = PETSc.@jacobian_fn(stokes_mass_J, dim_*dim_)
+const stokes_mass_J_ptr = PETSc.@jacobian_fn(petsclib, stokes_mass_J, dim_*dim_)
 
 # J_pp preconditioner (g0): (1/μ)⟨q, p⟩
 function stokes_id_J_kx(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
                           aOff, aOff_x, a, a_t, a_x, t, utShift, x, nC, cst, g0)
     g0[1] = 1.0 / exp(2.0 * real(cst[3]) * x[1])
 end
-const stokes_id_J_kx_ptr = PETSc.@jacobian_fn(stokes_id_J_kx, 1)
+const stokes_id_J_kx_ptr = PETSc.@jacobian_fn(petsclib, stokes_id_J_kx, 1)
 
 function stokes_id_J_cx(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
                           aOff, aOff_x, a, a_t, a_x, t, utShift, x, nC, cst, g0)
     mu = x[1] < real(cst[5]) ? real(cst[3]) : real(cst[4])
     g0[1] = 1.0 / mu
 end
-const stokes_id_J_cx_ptr = PETSc.@jacobian_fn(stokes_id_J_cx, 1)
+const stokes_id_J_cx_ptr = PETSc.@jacobian_fn(petsclib, stokes_id_J_cx, 1)
 
 # ── Exact solution / BC callbacks ─────────────────────────────────────────────
 # PetscSimplePointFn: f(t, x, u, ctx) → nothing
@@ -278,38 +276,38 @@ const stokes_id_J_cx_ptr = PETSc.@jacobian_fn(stokes_id_J_cx, 1)
 function zero_vel(t, x, u, ctx)
     for c in 1:length(u); u[c] = 0.0; end
 end
-const zero_vel_ptr = PETSc.@simple_fn(zero_vel)
+const zero_vel_ptr = PETSc.@simple_fn(petsclib, zero_vel)
 
 function one_pres(t, x, u, ctx)
     u[1] = 1.0
 end
-const one_pres_ptr = PETSc.@simple_fn(one_pres)
+const one_pres_ptr = PETSc.@simple_fn(petsclib, one_pres)
 
 # SolKx exact solution (captures parameters via module-level globals)
 function exact_vel_kx(t, x, u, ctx)
     vx, vz, _ = SolKxSolution(x, m_val, n_val, B_val)
     u[1] = vx; u[2] = vz
 end
-const exact_vel_kx_ptr = PETSc.@simple_fn(exact_vel_kx)
+const exact_vel_kx_ptr = PETSc.@simple_fn(petsclib, exact_vel_kx)
 
 function exact_pres_kx(t, x, u, ctx)
     _, _, p = SolKxSolution(x, m_val, n_val, B_val)
     u[1] = p
 end
-const exact_pres_kx_ptr = PETSc.@simple_fn(exact_pres_kx)
+const exact_pres_kx_ptr = PETSc.@simple_fn(petsclib, exact_pres_kx)
 
 # SolCx exact solution
 function exact_vel_cx(t, x, u, ctx)
     vx, vz, _ = SolCxSolution(x, m_val, n_val, xc_val, etaA_val, etaB_val)
     u[1] = vx; u[2] = vz
 end
-const exact_vel_cx_ptr = PETSc.@simple_fn(exact_vel_cx)
+const exact_vel_cx_ptr = PETSc.@simple_fn(petsclib, exact_vel_cx)
 
 function exact_pres_cx(t, x, u, ctx)
     _, _, p = SolCxSolution(x, m_val, n_val, xc_val, etaA_val, etaB_val)
     u[1] = p
 end
-const exact_pres_cx_ptr = PETSc.@simple_fn(exact_pres_cx)
+const exact_pres_cx_ptr = PETSc.@simple_fn(petsclib, exact_pres_cx)
 
 exact_vel_ptr  = sol_type == "solkx" ? exact_vel_kx_ptr  : exact_vel_cx_ptr
 exact_pres_ptr = sol_type == "solkx" ? exact_pres_kx_ptr : exact_pres_cx_ptr
diff --git a/scripts/renames.jl b/scripts/renames.jl
index 2c239d566..a7cba9363 100644
--- a/scripts/renames.jl
+++ b/scripts/renames.jl
@@ -741,6 +741,33 @@ const EXTRA_SHIMS = NamedTuple{(:name, :code, :test), Tuple{String, String, Stri
         end
         """,
     ),
+    (
+        name = "pointwise macros without a library",
+        code = """
+        # Warns in src/dmplex.jl: the one-argument-shorter forms of @simple_fn,
+        # @residual_fn, @jacobian_fn and @bd_fn take PetscInt, PetscScalar and
+        # PetscReal from the calling module.
+        """,
+        test = """
+        @testset "pointwise macros without a library" begin
+            m = Module()
+            Core.eval(m, quote
+                using PETSc
+                const PetscInt = Int64
+                const PetscScalar = Float64
+                const PetscReal = Float64
+                one_u(t, x, u, ctx) = (u[1] = x[1] + 1; nothing)
+            end)
+            ptr = warns(() -> Core.eval(m, :(PETSc.@simple_fn(one_u))))
+            x = [2.0]
+            u = [0.0]
+            err = ccall(ptr, Cint, (Int64, Float64, Ptr{Float64}, Int64, Ptr{Float64}, Ptr{Cvoid}),
+                        1, 0.0, x, 1, u, C_NULL)
+            @test err == 0
+            @test u[1] == 3
+        end
+        """,
+    ),
 ]
 
 # ---------------------------------------------------------------------------
@@ -752,6 +779,13 @@ const INTERNAL = Set{Symbol}([
     :checkout_local_arrays,
     # the per-library body of the MatShell constructor
     :create_matshell,
+    # the expansions behind the pointwise macros and their type lookup
+    :simple_fn_expr,
+    :residual_fn_expr,
+    :jacobian_fn_expr,
+    :bd_fn_expr,
+    :pointwise_types,
+    :pointwise_types_from_scope,
     # leading underscore dropped (naming.md "Internals")
     :petsc_link,
     :petsc_subst,
diff --git a/src/deprecations.jl b/src/deprecations.jl
index 04127261d..ea44a7b81 100644
--- a/src/deprecations.jl
+++ b/src/deprecations.jl
@@ -432,6 +432,11 @@ end
 # Warns in src/dm.jl and returns a copy; `with_local_coordinate_array` hands
 # the array back when its block returns.
 
+# pointwise macros without a library
+# Warns in src/dmplex.jl: the one-argument-shorter forms of @simple_fn,
+# @residual_fn, @jacobian_fn and @bd_fn take PetscInt, PetscScalar and
+# PetscReal from the calling module.
+
 
 # Renamed types: a `const` alias, so `isa` and dispatch keep working
 
diff --git a/src/dmplex.jl b/src/dmplex.jl
index f7ee17a91..962d0c922 100644
--- a/src/dmplex.jl
+++ b/src/dmplex.jl
@@ -13,20 +13,42 @@ Base.unsafe_convert(::Type{LibPETSc.PetscDS}, v::PetscDS) = v.ptr
 # These macros generate a thin C-callable wrapper that unwraps all C pointer
 # arguments into Julia Vectors, then calls the user's pure-Julia function.
 #
-# PetscInt, PetscReal, PetscScalar must be concrete types in the calling scope
-# (e.g. `const PetscInt = petsclib.PetscInt`).
+# The first argument names the library (a PetscLibType value or type). It is
+# evaluated when the macro expands, so it must be a global defined before the
+# macro call; the C signature uses that library's PetscInt, PetscReal and
+# PetscScalar.
 #
 # Usage:
+#   const petsclib = PETSc.getlib(; PetscScalar = Float64)
 #   function my_fn(t, x, u, ctx);  u[1] = x[1]^2 + x[2]^2;  end
-#   const my_fn_ptr = PETSc.@simple_fn(my_fn)
+#   const my_fn_ptr = PETSc.@simple_fn(petsclib, my_fn)
 #
 #   function f0(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
 #               aOff, aOff_x, a, a_t, a_x, t, x, numConstants, constants, out)
 #       out .= 4
 #   end
-#   const f0_ptr = PETSc.@residual_fn(f0, Nf)        # outsz = Nf
-#   const f1_ptr = PETSc.@residual_fn(f1, dim_*Nf)   # outsz = dim*Nf
-#   const g3_ptr = PETSc.@jacobian_fn(g3, dim_*dim_) # outsz = dim² (Nf=1)
+#   const f0_ptr = PETSc.@residual_fn(petsclib, f0, Nf)        # outsz = Nf
+#   const f1_ptr = PETSc.@residual_fn(petsclib, f1, dim_*Nf)   # outsz = dim*Nf
+#   const g3_ptr = PETSc.@jacobian_fn(petsclib, g3, dim_*dim_) # outsz = dim² (Nf=1)
+#
+# The wrapper's C pointer comes from an `@eval` at run time: within a single
+# top-level block the `@cfunction` is compiled before the wrapper exists.
+
+# Look up the PetscInt, PetscScalar and PetscReal of the library `lib` names
+# in module `mod`.
+function pointwise_types(mod::Module, lib)
+    L = Core.eval(mod, lib)
+    L isa LibPETSc.PetscLibType || (L isa Type && L <: LibPETSc.PetscLibType && isconcretetype(L)) ||
+        throw(ArgumentError("the first argument of a pointwise macro must name a PETSc library or its type, got `$lib` = $L"))
+    return LibPETSc.inttype(L), LibPETSc.scalartype(L), LibPETSc.realtype(L)
+end
+
+# Old form: the types are the globals PetscInt, PetscScalar and PetscReal of `mod`.
+# TODO: remove in v0.6 with the other 0.5 deprecations.
+function pointwise_types_from_scope(mod::Module, name::Symbol, f)
+    @warn "$name without a library is deprecated, use $name(petsclib, $f, …)" maxlog = 1 _id = name
+    return Core.eval(mod, :PetscInt), Core.eval(mod, :PetscScalar), Core.eval(mod, :PetscReal)
+end
 
 # Helper: substitute gensymmed param names into an output-size expression.
 petsc_subst(ex, m) =
@@ -35,7 +57,7 @@ petsc_subst(ex, m) =
     ex
 
 """
-    @simple_fn(f)
+    @simple_fn(petsclib, f)
 
 Generate a C-callable wrapper for `f(t, x, u, ctx)` and return its C pointer,
 matching the `PetscSimplePointFn` signature.  Arguments seen by `f`:
@@ -45,14 +67,20 @@ matching the `PetscSimplePointFn` signature.  Arguments seen by `f`:
   - `ctx` — user context (`Ptr{Cvoid}`)
 
 Used by `project_function!`, `add_boundary!`, `set_exact_solution!`.
-`PetscInt`, `PetscReal`, `PetscScalar` must be in scope.
+`petsclib` is a PETSc library or its type, and must be a global defined before
+the macro call: it sets the C types of the arguments.
 """
+macro simple_fn(lib, f)
+    return simple_fn_expr(f, pointwise_types(__module__, lib)...)
+end
+
 macro simple_fn(f)
+    return simple_fn_expr(f, pointwise_types_from_scope(__module__, Symbol("@simple_fn"), f)...)
+end
+
+function simple_fn_expr(f, PI, PS, PR)
     cfn = Symbol(f, :__cfn__)
     gdim = gensym("dim"); gNc = gensym("Nc")
-    PI = Core.eval(__module__, :PetscInt)
-    PS = Core.eval(__module__, :PetscScalar)
-    PR = Core.eval(__module__, :PetscReal)
     cfn_ptr_expr = :(Base.@cfunction($cfn, Cint,
         ($PI, $PR, Ptr{$PR}, $PI, Ptr{$PS}, Ptr{Cvoid})))
     quote
@@ -70,7 +98,7 @@ macro simple_fn(f)
 end
 
 """
-    @residual_fn(f, outsz)
+    @residual_fn(petsclib, f, outsz)
 
 Generate a C-callable wrapper for the pure-Julia residual function `f` and
 return its C pointer, matching the `PetscPointFn` signature.  All pointer
@@ -85,16 +113,21 @@ User function signature:
 f(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
   aOff, aOff_x, a, a_t, a_x, t, x, numConstants, constants, out)
 ```
-`PetscInt`, `PetscReal`, `PetscScalar` must be in scope.
+`petsclib` is as in `@simple_fn`.
 """
+macro residual_fn(lib, f, outsz)
+    return residual_fn_expr(f, outsz, pointwise_types(__module__, lib)...)
+end
+
 macro residual_fn(f, outsz)
+    return residual_fn_expr(f, outsz, pointwise_types_from_scope(__module__, Symbol("@residual_fn"), f)...)
+end
+
+function residual_fn_expr(f, outsz, PI, PS, PR)
     cfn = Symbol(f, :__cfn__)
     gdim = gensym("dim"); gNf = gensym("Nf"); gNfAux = gensym("NfAux"); gnC = gensym("nC")
     m = Dict{Symbol,Any}(:dim_ => gdim, :Nf => gNf, :NfAux => gNfAux, :numConstants => gnC)
     outsz_s = petsc_subst(outsz, m)
-    PI = Core.eval(__module__, :PetscInt)
-    PS = Core.eval(__module__, :PetscScalar)
-    PR = Core.eval(__module__, :PetscReal)
     cfn_ptr_expr = :(Base.@cfunction($cfn, Cvoid,
         ($PI, $PI, $PI, Ptr{$PI}, Ptr{$PI},
          Ptr{$PS}, Ptr{$PS}, Ptr{$PS}, Ptr{$PI}, Ptr{$PI},
@@ -132,7 +165,7 @@ macro residual_fn(f, outsz)
 end
 
 """
-    @jacobian_fn(f, outsz)
+    @jacobian_fn(petsclib, f, outsz)
 
 Generate a C-callable wrapper for the pure-Julia Jacobian function `f` and
 return its C pointer, matching the `PetscPointJacFn` signature.  Identical to
@@ -143,16 +176,21 @@ User function signature:
 f(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
   aOff, aOff_x, a, a_t, a_x, t, u_tShift, x, numConstants, constants, out)
 ```
-`PetscInt`, `PetscReal`, `PetscScalar` must be in scope.
+`petsclib` is as in `@simple_fn`.
 """
+macro jacobian_fn(lib, f, outsz)
+    return jacobian_fn_expr(f, outsz, pointwise_types(__module__, lib)...)
+end
+
 macro jacobian_fn(f, outsz)
+    return jacobian_fn_expr(f, outsz, pointwise_types_from_scope(__module__, Symbol("@jacobian_fn"), f)...)
+end
+
+function jacobian_fn_expr(f, outsz, PI, PS, PR)
     cfn = Symbol(f, :__cfn__)
     gdim = gensym("dim"); gNf = gensym("Nf"); gNfAux = gensym("NfAux"); gnC = gensym("nC")
     m = Dict{Symbol,Any}(:dim_ => gdim, :Nf => gNf, :NfAux => gNfAux, :numConstants => gnC)
     outsz_s = petsc_subst(outsz, m)
-    PI = Core.eval(__module__, :PetscInt)
-    PS = Core.eval(__module__, :PetscScalar)
-    PR = Core.eval(__module__, :PetscReal)
     cfn_ptr_expr = :(Base.@cfunction($cfn, Cvoid,
         ($PI, $PI, $PI, Ptr{$PI}, Ptr{$PI},
          Ptr{$PS}, Ptr{$PS}, Ptr{$PS}, Ptr{$PI}, Ptr{$PI},
@@ -191,7 +229,7 @@ macro jacobian_fn(f, outsz)
 end
 
 """
-    @bd_fn(f, outsz)
+    @bd_fn(petsclib, f, outsz)
 
 Generate a C-callable wrapper for the pure-Julia boundary pointwise function `f`
 and return its C pointer, matching the `PetscBdPointFn` signature.  Identical to
@@ -203,16 +241,21 @@ User function signature:
 f(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
   aOff, aOff_x, a, a_t, a_x, t, x, n, numConstants, constants, out)
 ```
-`PetscInt`, `PetscReal`, `PetscScalar` must be in scope.
+`petsclib` is as in `@simple_fn`.
 """
+macro bd_fn(lib, f, outsz)
+    return bd_fn_expr(f, outsz, pointwise_types(__module__, lib)...)
+end
+
 macro bd_fn(f, outsz)
+    return bd_fn_expr(f, outsz, pointwise_types_from_scope(__module__, Symbol("@bd_fn"), f)...)
+end
+
+function bd_fn_expr(f, outsz, PI, PS, PR)
     cfn = Symbol(f, :__cfn__)
     gdim = gensym("dim"); gNf = gensym("Nf"); gNfAux = gensym("NfAux"); gnC = gensym("nC")
     m = Dict{Symbol,Any}(:dim_ => gdim, :Nf => gNf, :NfAux => gNfAux, :numConstants => gnC)
     outsz_s = petsc_subst(outsz, m)
-    PI = Core.eval(__module__, :PetscInt)
-    PS = Core.eval(__module__, :PetscScalar)
-    PR = Core.eval(__module__, :PetscReal)
     cfn_ptr_expr = :(Base.@cfunction($cfn, Cvoid,
         ($PI, $PI, $PI, Ptr{$PI}, Ptr{$PI},
          Ptr{$PS}, Ptr{$PS}, Ptr{$PS}, Ptr{$PI}, Ptr{$PI},
diff --git a/test/core/deprecations.jl b/test/core/deprecations.jl
index 69e418bde..fac6caa63 100644
--- a/test/core/deprecations.jl
+++ b/test/core/deprecations.jl
@@ -380,4 +380,22 @@ const SHIMS = [
         @test coord[1, 3] ≈ 1
     end
 
+    @testset "pointwise macros without a library" begin
+        m = Module()
+        Core.eval(m, quote
+            using PETSc
+            const PetscInt = Int64
+            const PetscScalar = Float64
+            const PetscReal = Float64
+            one_u(t, x, u, ctx) = (u[1] = x[1] + 1; nothing)
+        end)
+        ptr = warns(() -> Core.eval(m, :(PETSc.@simple_fn(one_u))))
+        x = [2.0]
+        u = [0.0]
+        err = ccall(ptr, Cint, (Int64, Float64, Ptr{Float64}, Int64, Ptr{Float64}, Ptr{Cvoid}),
+                    1, 0.0, x, 1, u, C_NULL)
+        @test err == 0
+        @test u[1] == 3
+    end
+
 end
diff --git a/test/dm/dmplex.jl b/test/dm/dmplex.jl
index f47a29d9e..e9fc1814c 100644
--- a/test/dm/dmplex.jl
+++ b/test/dm/dmplex.jl
@@ -10,18 +10,15 @@ using MPI
 
 MPI.Initialized() || MPI.Init()
 
-# ── PETSc library and scalar types ───────────────────────────────────────────
-# @simple_fn / @residual_fn / @jacobian_fn resolve
-# PetscInt / PetscScalar / PetscReal from the calling module at macro-expansion
-# time, so they MUST be module-level names before the macros are invoked.
-# We define them once for the Float64/Int64 lib that owns the FEM callbacks.
+# ── PETSc library ────────────────────────────────────────────────────────────
+# @simple_fn / @residual_fn / @jacobian_fn evaluate their library argument at
+# macro-expansion time, so it must be a module-level name before the macros
+# are invoked. The FEM callbacks use the Float64/Int64 lib.
 
 const _dmplex_petsclib = PETSc.petsclibs[1]
 PETSc.initialize(_dmplex_petsclib)
 
 const PetscInt    = _dmplex_petsclib.PetscInt
-const PetscScalar = Float64
-const PetscReal   = Float64
 
 # Serial-safe communicator (shared across all lib loops)
 const _TC = MPI.COMM_SELF
@@ -46,7 +43,7 @@ function _dm_f0(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
                 aOff, aOff_x, a, a_t, a_x, t, x, nC, cst, f0)
     f0[1] = 0.0
 end
-const _dm_f0_ptr = PETSc.@residual_fn(_dm_f0, 1)
+const _dm_f0_ptr = PETSc.@residual_fn(_dmplex_petsclib, _dm_f0, 1)
 
 # f1: diffusion flux ∇u
 function _dm_f1(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
@@ -55,7 +52,7 @@ function _dm_f1(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
         f1[d] = u_x[d]
     end
 end
-const _dm_f1_ptr = PETSc.@residual_fn(_dm_f1, dim_)
+const _dm_f1_ptr = PETSc.@residual_fn(_dmplex_petsclib, _dm_f1, dim_)
 
 # g3: identity stiffness tensor (Laplacian Jacobian)
 function _dm_g3(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
@@ -65,7 +62,7 @@ function _dm_g3(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
         g3[d * dim_ + d + 1] = 1.0
     end
 end
-const _dm_g3_ptr = PETSc.@jacobian_fn(_dm_g3, dim_ * dim_)
+const _dm_g3_ptr = PETSc.@jacobian_fn(_dmplex_petsclib, _dm_g3, dim_ * dim_)
 
 # exact solution: u(t, x) = x₁ + x₂  (linear → exact in P1/Q1)
 function _dm_exact(t, x, u, ctx)
@@ -74,13 +71,17 @@ function _dm_exact(t, x, u, ctx)
         u[1] += x[d]
     end
 end
-const _dm_exact_ptr = PETSc.@simple_fn(_dm_exact)
+const _dm_exact_ptr = PETSc.@simple_fn(_dmplex_petsclib, _dm_exact)
 
 # zero BC function
 function _dm_zero(t, x, u, ctx)
     fill!(u, 0.0)
 end
-const _dm_zero_ptr = PETSc.@simple_fn(_dm_zero)
+const _dm_zero_ptr = PETSc.@simple_fn(_dmplex_petsclib, _dm_zero)
+
+# The library's type names the C types as well as the library does.
+_dm_one(t, x, u, ctx) = (u[1] = x[1] + 1; nothing)
+const _dm_one_ptr = PETSc.@simple_fn(typeof(_dmplex_petsclib), _dm_one)
 
 # ─────────────────────────────────────────────────────────────────────────────
 @testset "DMPlex" begin
@@ -89,6 +90,18 @@ if isintelmac
     @test_skip "DMPlex+PETSc tests skipped on macOS Intel (known SIGSEGV in PETSc_jll)"
 else
 
+@testset "pointwise macro library argument" begin
+    @test _dmplex_petsclib.PetscScalar === Float64
+    x = [2.0]
+    u = [0.0]
+    err = ccall(_dm_one_ptr, Cint, (PetscInt, Float64, Ptr{Float64}, PetscInt, Ptr{Float64}, Ptr{Cvoid}),
+                1, 0, x, 1, u, C_NULL)
+    @test err == 0
+    @test u[1] == 3
+    @test_throws "must name a PETSc library" macroexpand(@__MODULE__, :(PETSc.@simple_fn(1, _dm_one)))
+end
+
+
 for petsclib in PETSc.petsclibs
     PETSc.initialize(petsclib)
     PetscInt_t    = petsclib.PetscInt