Some revisions

README.md

@@ -1,6 +1,31 @@
 # ProximalGalerkin
 
-Examples of the proximal Galerkin finite element method.
+This repository contains implementations of the proximal Galerkin finite element method and other proximal numerical methods for variational problems with inequality constraints derived in
+
+```bibtex
+@misc{dokken2025latent,
+      title={The latent variable proximal point algorithm for variational problems with constraints}, 
+      author={Dokken, {J\o rgen} S. and Farrell, Patrick~E. and Keith, Brendan and Papadopoulos, Ioannis~P.A. and Surowiec, Thomas~M.},
+      year={2025},
+}
+```
+
+Please cite the aforementioned manuscript if using the code in this repository.
+
+## Instructions
+
+We encourage using following Docker containers to run the  codes described below:
+
+- DOLFINx: [ghcr.io/methods-group/proximalgalerkin-dolfinx:main](https://github.com/METHODS-Group/ProximalGalerkin/pkgs/container/proximalgalerkin-dolfinx)
+- MFEM: [ghcr.io/methods-group/proximalgalerkin-mfem:main](https://github.com/METHODS-Group/ProximalGalerkin/pkgs/container/proximalgalerkin-mfem)
+- Firedrake: [ghcr.io/methods-group/proximalgalerkin-firedrake:main](https://github.com/METHODS-Group/ProximalGalerkin/pkgs/container/proximalgalerkin-firedrake)
+- Julia/GridAP: [julia:1.10.8](https://hub.docker.com/layers/library/julia/1.10.8/images/sha256-66656909ed7b5e75f4208631b01fc585372f906d68353d97cc06b40a8028c437)
+
+<a name="obstacle"></a>
+
+## Table of Examples and Figures
+
+The following table associates each implementation to the examples and figures in the paper. Further information to run the codes is provided for each specific example.
 
 | Figure |                                        File: examples/                                        |     Backend      | Instructions                     |
 | :----: | :-------------------------------------------------------------------------------------------: | :--------------: | -------------------------------- |
@@ -15,25 +40,14 @@ Examples of the proximal Galerkin finite element method.
 |   8    |                                               ?                                               |    Firedrake     |                                  |
 |   9    |                                               ?                                               |    Firedrake     |                                  |
 |   10   |            [harmonic_maps/harmonic_1d.py](./examples/harmonic_maps/harmonic_1d.py)            |      FEniCS      | [Harmonic map](#harmonic)        |
-|   11   |                        [eikonal/ex40.cpp](./examples/eikonal/ex40.cpp)                        |       MFEM       | [Eikonal](#eikonal)              |
+|   11   |                        [eikonal/ex40.cpp](./examples/eikonal/ex40.cpp)                        |       MFEM/FEniCS       | [Eikonal](#eikonal)              |
 |   12   |                [monge_ampere/cg_cg_dg.py](./examples/monge_ampere/cg_cg_dg.py)                | Firedrake/FEniCS | [Monge-Ampere](#monge)           |
 
-## Installation
+## Example 1 (Figure 2): The Obstacle Problem
 
-We provide several docker container for the various methods used in the paper
+Figures 2 (a) and (b) are generated with `DOLFINx`.
 
-- DOLFINx: [ghcr.io/methods-group/proximalgalerkin-dolfinx:main](https://github.com/METHODS-Group/ProximalGalerkin/pkgs/container/proximalgalerkin-dolfinx)
-- MFEM: [ghcr.io/methods-group/proximalgalerkin-mfem:main](https://github.com/METHODS-Group/ProximalGalerkin/pkgs/container/proximalgalerkin-mfem)
-- Firedrake: [ghcr.io/methods-group/proximalgalerkin-firedrake:main](https://github.com/METHODS-Group/ProximalGalerkin/pkgs/container/proximalgalerkin-firedrake)
-- Julia/GridAP: [julia:1.10.8](https://hub.docker.com/layers/library/julia/1.10.8/images/sha256-66656909ed7b5e75f4208631b01fc585372f906d68353d97cc06b40a8028c437)
-
-<a name="obstacle"></a>
-
-# Obstacle problem
-
-Parts of this table is generated using DOLFINx.
-
-To get the results from Galahad, IPOPT, SNES and LVPP (FEM) use the LVPP docker image for DOLFINx and run (within `examples/obstacle`)
+To reproduce the results in Figures 2 (a) (the comparison between Proximal Galerkin, SNES, Galahad, and IPOPT), first deploy the `DOLFINx` Docker container. Then run the following commands within `examples/obstacle`:
 
 ```bash
 python3 generate_mesh_gmsh.py
@@ -42,44 +56,45 @@ python3 compare_all.py -P ./meshes/disk_2.xdmf -O medium
 python3 compare_all.py -P ./meshes/disk_3.xdmf -O fine
 ```
 
-Julia is used to get results for the finite difference and spectral element method, use the `julia:1.10.8` Docker container and call
+To reproduce the finite difference and spectral element method results in Figure 2 (c), deploy the `julia:1.10.8` Docker container and call
 
 ```bash
 julia finite_difference.jl
 julia spectral.jl
 ```
 
-<a name="signorini"></a>
+within `examples/obstacle`.
 
-# Signorini problem
+<a name="signorini"></a>
 
-Requires DOLFINx and [scifem](https://github.com/scientificcomputing/scifem). These are installed in the docker image.
+## Example 2 (Figure 3): The Signorini Problem
 
-From within `examples/signorini`, call
+Deploy the `DOLFINx` Docker container to reproduce the results in this example.
+Then call
 
 ```bash
 python3 generate_mesh.py
 ```
 
-to generate the mesh file `"meshes/half_sphere.xdmf"`.
-Next run the LVPP algorithm with
+from within `examples/signorini` to generate the mesh file `"meshes/half_sphere.xdmf"`.
+Next, run the proximal Galerkin method with
 
 ```bash
 python3 run_lvpp_problem.py --alpha_0=0.005 --degree=2 --disp=-0.3 --n-max-iterations=250 --alpha_scheme=doubling  --output output_lvpp file --filename=meshes/half_sphere.xdmf
 ```
 
 <a name="fracture"></a>
 
-# Fracture
+## Example 3 (Figure 4): Variational Fracture
 
-Can be simulated with either FEniCS/DOLFINx or Firedrake.
-The DOLFINx code can be executed with
+This example can be run from within `examples/fracture` using both the `DOLFINx` and the `Firedrake` Docker containers.
+The `DOLFINx` code can be executed with
 
 ```bash
 python3 script.py
 ```
 
-while the Firedrake code can be executed with:
+while the `Firedrake` code can be executed with:
 
 > [!WARNING]  
 > Add instructions
@@ -90,65 +105,96 @@ while the Firedrake code can be executed with:
 
 <a name="ch"></a>
 
-# Cahn-Hilliard problem
+## Example 4 (Figure 5): Four-Phase Cahn–Hilliard Gradient Flow
 
-Codes can run from within `examples/cahn-hilliard` with
+To reproduce the results in this example, first deploy the `DOLFINx` Docker container.
+Then run
 
 ```bash
 python3 problem.py
 ```
 
+from within `examples/cahn-hilliard`.
+
 <a name="qvi"></a>
 
-# Thermoforming quasi-variational inequalities
+## Example 5 (Figure 6): Thermoforming Quasi-Variational Inequality
 
-Requires Julia. Can be executed with
+Reproducing the results in this example requires the `julia:1.10.8` Docker container.
+Once this container is deployed, the code can be executed by running
 
 ```bash
 julia theroforming_lvpp.jl
 ```
 
-from the `examples/thermoforming_qvi` folder.
+from within `examples/thermoforming_qvi`.
 
 <a name="gradient"></a>
 
-# Gradient constraint
+## Example 6 (Figure 7): Gradient Norm Constraints
 
-Run the `script.py` with the following input parameters:
+Deploy the `DOLFINx` Docker container to reproduce the results in this example.
+Then run `script.py` within `examples/gradient_constraint` with the following input parameters:
 
 ```bash
 python3 script.py -N 80 -M 80 --alpha_scheme=doubling
 ```
 
 <a name="harmonic"></a>
 
-## Harmonic maps
+## Example 7 (Figure 8): Eigenvalue Constraints
+
+Deploy the `Firedrake` Docker container to reproduce the results in this example.
+Then run the following command within `examples/[TODO]`:
+
+> [!WARNING]  
+> Add instructions 
+
+## Example 8 (Figure 9): Intersections of Constraints
 
-Requires DOLFINx. Run
+Deploy the `Firedrake` Docker container to reproduce the results in this example.
+Then run the following command within `examples/[TODO]`:
+
+> [!WARNING]  
+> Add instructions 
+
+## Example 9 (Figure 10): Harmonic Maps to the Sphere
+
+Deploy the `DOLFINx` Docker container to reproduce the results in this example.
+Then run the following command within `examples/harmonic_maps`:
 
 ```bash
 python3 harmonic_1D.py
 ```
 
-## Eikonal equation
+## Example 10: Linear Equality Constraints
+Note that there is no numerical example for this setting because the derived variational formulation is equivalent to the standard Lagrange multiplier formulation for this class of problems.
+
+## Example 11 (Figure 11): The Eikonal Equation
+
+We have provided code for this example for both the `MFEM` and `DOLFINx` Docker containers.
 
-The MFEM example can be executed by copying
-[./examples/eikonal/ex40.cpp](./examples/eikonal/ex40.cpp) into the `mfem` examples folder
-and call `make ex40`. It can then be executed with:
+To reproduce the Möbius strip solution in Figure 11, you first need to copy [./examples/eikonal/ex40.cpp](./examples/eikonal/ex40.cpp) into the `MFEM` examples folder (`/home/euler/mfem/examples/`) and then calling `make ex40` and `./ex40 -step 10.0 -mi 10`. This following code will execute to entire process:
 
 ```bash
 docker run -it --rm -v ./examples/eikonal:/home/euler/shared -w /home/euler/mfem --rm --entrypoint=/bin/bash ghcr.io/methods-group/proximalgalerkin-mfem:main
 cp /home/euler/shared/ex40.cpp /home/euler/mfem/examples/
 cd examples && make ex40
-./ex40
+./ex40 -step 10.0 -mi 10
 ```
 
-For the non-manifold examples, with the [Star](https://github.com/mfem/mfem/blob/master/data/star.mesh)
-and [Ball](https://github.com/mfem/mfem/blob/master/data/ball-nurbs.mesh) you can compile the [official demo](https://mfem.org/examples/), `ex40.cpp` or `ex40p.cpp` without copying any files from this repository.
-
-The DOLFINx example, in [./examples/eikonal/script.py](./examples/eikonal/script.py) requires to convert the mobius strip mesh from `mfem`, called [mobius-strip.mesh](https://github.com/mfem/mfem/blob/master/data/mobius-strip.mesh)
+To reproduce the results in Figure 11 for the two geometries (i.e., the [Star](https://github.com/mfem/mfem/blob/master/data/star.mesh)
+and [Ball](https://github.com/mfem/mfem/blob/master/data/ball-nurbs.mesh)), you should compile the [official examples](https://mfem.org/examples/) `ex40.cpp` or `ex40p.cpp` without copying any files from this repository
+```bash
+cd examples && make ex40
+# Star Geometry
+./ex40 -step 10.0 -mi 10
+# Ball Geometry
+./ex40 -step 10.0 -mi 10 -m ../data/ball-nurbs.mesh
+```
 
-From the root of the repository, you can call the following commands to compile the code
+The `DOLFINx` implementation, found in [./examples/eikonal/script.py](./examples/eikonal/script.py) requires first converting the `MFEM` Möbius strip mesh [mobius-strip.mesh](https://github.com/mfem/mfem/blob/master/data/mobius-strip.mesh).
+To this end, run the following commands from the root of this repository:
 
 ```bash
 docker run -it --rm -v ./examples/eikonal:/home/euler/shared -w /home/euler/mfem --rm --entrypoint=/bin/bash ghcr.io/methods-group/proximalgalerkin-mfem:main
@@ -158,17 +204,27 @@ cd examples && make convert_mesh
 cp -r  mobius-strip.mesh/ ../../shared/
 ```
 
+The `DOLFINx` code is then executed by calling:
+
+```bash
+python3 script.py
+```
+
+from within `examples/eikonal`.
+
 <a name="monge"></a>
 
-# Monge-Ampere
+## Example 12 (Figure 12): The Monge–Ampere Equation
+
+This example can be run from within `examples/monge_ampere` using both the `DOLFINx` and the `Firedrake` Docker containers.
 
-The firedrake code can be run with
+The `Firedrake` code can be run with the command
 
 ```bash
 python3 cg_cg_dg.py
 ```
 
-The equivalent FEniCS/DOLFINx code can be run with
+The equivalent `DOLFINx` code can be run with
 
 ```bash
 python3 cg_cg_dg_fenics.py

examples/eikonal/ex40.cpp

@@ -314,7 +314,7 @@ int main(int argc, char *argv[])
         }
 
         alpha *= max(growth_rate, 1_r);
-        alpha = min(alpha, 50.0);
+        alpha = min(alpha, 50_r);
         neg_alpha_cf.constant = -alpha;
     }
 

examples/gradient_constraint/script.py

@@ -214,7 +214,7 @@ def main(
     element_options.add_argument(
         "--latent_space",
         type=str,
-        default="RT",
+        default="Lagrange",
         choices=get_args(latent_spaces),
         help="Finite element family for auxiliary variable",
     )

examples/monge_ampere/cg_cg_dg_fenics.py

@@ -29,7 +29,7 @@ def extract_num_dofs(V: dolfinx.fem.FunctionSpace):
 
 
 errors = []
-for i in range(4):
+for j in range(3,15):
     X = (x, y) = ufl.SpatialCoordinate(mesh)
     gdim = mesh.geometry.dim
 
@@ -48,10 +48,10 @@ def hessian(u):
     rho = ufl.det(hessian(u_exact))  # forcing term
     g = u_exact  # boundary data
 
-    k = 6
+    k = j
     el_V = basix.ufl.element("Lagrange", mesh.basix_cell(), k)
     el_U = basix.ufl.element("Lagrange", mesh.basix_cell(), k + 1, shape=(gdim,))
-    el_W = basix.ufl.element("Discontinuous Lagrange", mesh.basix_cell(), k, shape=(3,))
+    el_W = basix.ufl.element("Lagrange", mesh.basix_cell(), k, shape=(3,))
     me = basix.ufl.mixed_element([el_V, el_U, el_W])
     Z = dolfinx.fem.functionspace(mesh, me)
 
@@ -163,21 +163,15 @@ def hessian(u):
     global_error = np.sqrt(mesh.comm.allreduce(local_error, op=MPI.SUM))
     errors.append(global_error)
 
-    skew = 0.5 * (psi - psi.T)
-    skew_psi = ufl.inner(skew, skew) * dx
-    skew_local = dolfinx.fem.assemble_scalar(dolfinx.fem.form(skew_psi))
-    skew_global = np.sqrt(mesh.comm.allreduce(skew_local, op=MPI.SUM))
-    print("||skew(psi)||_L2: ", skew_global, flush=True)
-
-    with dolfinx.io.VTXWriter(mesh.comm, f"output/solution_{i}.bp", [z.sub(0).collapse()]) as vtx:
+    with dolfinx.io.VTXWriter(mesh.comm, f"output/solution_{j}.bp", [z.sub(0).collapse()]) as vtx:
         vtx.write(0.0)
 
-    mesh, _, _ = dolfinx.mesh.refine(mesh)
-
+    # mesh, _, _ = dolfinx.mesh.refine(mesh)
 
-def convergence_orders(x):
-    return np.log2(np.array(x)[:-1] / np.array(x)[1:])
 
+# Not relevant for p-refinement
+# def convergence_orders(x):
+#     return np.log2(np.array(x)[:-1] / np.array(x)[1:])
 
 print("Errors", errors, flush=True)
-print("Convergence orders: ", convergence_orders(errors), flush=True)
+# print("Convergence orders: ", convergence_orders(errors), flush=True)

examples/obstacle/compare_all.py

@@ -1,5 +1,5 @@
 """
-Solve obstacle problem with DOLFINx (LVPP), Galahad and IPOPT and compare the results
+Solve obstacle problem with Proximal Galerkin, Galahad, and IPOPT and compare the results
 
 Author: Jørgen S. Dokken
 SPDX-License-Identifier: MIT

examples/obstacle/lvpp_example.py

@@ -1,20 +1,17 @@
 """
 
-Obstacle problem based of example 3 (P3) in Hintermüller and Kunisch [1].
+Obstacle problem based on experiment 4 in [2].
 
-FEniCSx code solve this problem is based of the paper [2,3]:
+The FEniCSx code solve this problem is based on [3]:
 
 SPXD License: MIT License
 
 Original license file [../../licenses/LICENSE.surowiec](../../licenses/LICENSE.surowiec)
 is included in the repository.
 
-[1] Hintermüller, M. and Kunisch K., Path-following Methods for a Class of
-    Constrained Minimization Problems in Function Space, SIAM Journal on Optimization 2006,
-    https://doi.org/10.1137/040611598
-[2] Keith, B. and Surowiec, T.M., Proximal Galerkin: A Structure-Preserving Finite Element Method
+[1] Keith, B. and Surowiec, T.M., Proximal Galerkin: A Structure-Preserving Finite Element Method
 for Pointwise Bound Constraints. Found Comput Math (2024). https://doi.org/10.1007/s10208-024-09681-8
-[3] Keith, B., Surowiec, T. M., & Dokken, J. S. (2023). Examples for the Proximal Galerkin Method
+[2] Keith, B., Surowiec, T. M., & Dokken, J. S. (2023). Examples for the Proximal Galerkin Method
     (Version 0.1.0) [Computer software]. https://github.com/thomas-surowiec/proximal-galerkin-examples
 """