Putting this into a dev branch. We'll need to check against v0.99

.gitignore

@@ -13,6 +13,12 @@ Jupyterbook/**/*.png
 Jupyterbook/**/*.vtk
 Jupyterbook/**/*.h5
 
+**/WIP/**/*.png
+**/BatPumpkin/**/*.png
+
+**/?????_debug.txt
+**/Untitled*.ipynb
+
 **/.jupyter/**
 
 # Byte-compiled / optimized / DLL files

Documentation/Ex_Anisotropy.ipynb

@@ -757,7 +757,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.11.5"
+   "version": "3.12.8"
   },
   "title": "Constitutive Models in Underworld 3"
  },

Notebooks/Developers/Timing/Timing_NonLinear_ShearBand.py

@@ -5,7 +5,7 @@
 #       extension: .py
 #       format_name: light
 #       format_version: '1.5'
-#       jupytext_version: 1.14.4
+#       jupytext_version: 1.16.6
 #   kernelspec:
 #     display_name: Python 3 (ipykernel)
 #     language: python

Notebooks/Developers/Timing/Timing_StokesSinker.py

@@ -5,7 +5,7 @@
 #       extension: .py
 #       format_name: light
 #       format_version: '1.5'
-#       jupytext_version: 1.14.1
+#       jupytext_version: 1.16.6
 #   kernelspec:
 #     display_name: Python 3 (ipykernel)
 #     language: python
@@ -200,7 +200,5 @@
     time += dt
 
 
-# + tags=[]
 timing.print_table(display_fraction=0.999)
-# -
 mesh.write_checkpoint("stokesSinkerTiming", meshUpdates=False, meshVars=[p, v])

Notebooks/Developers/WIP/ParallelEvaluate-viz.py

@@ -0,0 +1,138 @@
+# ---
+# jupyter:
+#   jupytext:
+#     formats: py:percent
+#     text_representation:
+#       extension: .py
+#       format_name: percent
+#       format_version: '1.3'
+#       jupytext_version: 1.16.6
+#   kernelspec:
+#     display_name: Python (Pixi)
+#     language: python
+#     name: pixi-kernel-python3
+# ---
+
+# %%
+# This is required to fix pyvista
+# (visualisation) crashes in interactive notebooks (including on binder)
+
+import nest_asyncio
+nest_asyncio.apply()
+
+# %%
+import numpy as np
+import sympy
+import underworld3 as uw
+
+# %%
+index = 12
+basename = "output/Parallel_Evaluate"
+
+# %%
+meshball = uw.discretisation.Mesh(f"{basename}.mesh.{index:05d}.h5")
+
+# meshball.view()
+
+x, y = meshball.CoordinateSystem.X
+r, th = meshball.CoordinateSystem.xR
+unit_rvec = meshball.CoordinateSystem.unit_e_0
+
+
+# Orientation of surface normals
+Gamma_N = meshball.Gamma
+
+
+# %%
+scalar_var = uw.discretisation.MeshVariable(
+    varname="Radius", mesh=meshball, vtype=uw.VarType.SCALAR, varsymbol=r"r"
+)
+
+rank_var = uw.discretisation.MeshVariable(
+    varname="Rank",
+    mesh=meshball,
+    vtype=uw.VarType.SCALAR,
+    varsymbol=r"Ra",
+    degree=0,
+)
+
+# %%
+scalar_var.read_timestep(f"{basename}", "Radius", index, verbose=True)
+rank_var.read_timestep(f"{basename}", "Rank", index, verbose=True)
+
+# %%
+
+# %%
+swarm = uw.swarm.Swarm(meshball)
+
+error_var = uw.swarm.SwarmVariable(
+    "Error",
+    swarm,
+    vtype=uw.VarType.SCALAR,
+    _proxy=False,
+)
+
+s_rank_var = uw.swarm.SwarmVariable(
+    "FRank",
+    swarm,
+    vtype=uw.VarType.SCALAR,
+    _proxy=False,
+)
+
+
+swarm.read_timestep(
+    basename,
+    "data_swarm",
+    index,
+    migrate=False,
+)
+
+error_var.read_timestep(
+    basename,
+    "data_swarm",
+    "Error",
+    index,
+)
+
+s_rank_var.read_timestep(
+    basename,
+    "data_swarm",
+    "FRank",
+    index,
+)
+
+
+
+# %%
+in_or_out = meshball.points_in_domain(swarm.particle_coordinates.array[...].reshape(-1,2))
+
+# %%
+
+# %%
+import pyvista as pv
+
+pvmesh = uw.visualisation.mesh_to_pv_mesh(meshball)
+
+pvmesh.cell_data["rank"] = uw.visualisation.scalar_fn_to_pv_points(
+    pvmesh.cell_centers(), rank_var.sym
+)
+
+error_swarm = uw.visualisation.swarm_to_pv_cloud(swarm)
+error_swarm.point_data["E"] = np.sqrt(error_var.array[:, 0, 0]**2)
+error_swarm.point_data["E2"] = np.sqrt(error_var.array[:, 0, 0]**2) * in_or_out.astype(float)
+error_swarm.point_data["R"] = np.sqrt(s_rank_var.array[:, 0, 0]**2) * in_or_out.astype(float)
+
+
+
+plotter = pv.Plotter()
+
+plotter.add_mesh(pvmesh, scalars="rank", cmap="rainbow", show_edges=True, opacity=0.2)
+plotter.add_points(error_swarm, cmap="rainbow", scalars="R", point_size=4)
+# plotter.add_points(error_swarm,  scalars="E", point_size=3)
+
+plotter.show()
+
+# %%
+swarm.dm.getLocalSize()
+
+# %%

Notebooks/Developers/WIP/ParallelEvaluate.py

@@ -0,0 +1,142 @@
+# ---
+# jupyter:
+#   jupytext:
+#     formats: py:percent
+#     text_representation:
+#       extension: .py
+#       format_name: percent
+#       format_version: '1.3'
+#       jupytext_version: 1.16.6
+#   kernelspec:
+#     display_name: Python (Pixi)
+#     language: python
+#     name: pixi-kernel-python3
+# ---
+
+# %%
+import numpy as np
+import sympy
+import underworld3 as uw
+
+# %%
+res = 0.05
+r_o = 1.0
+r_int = 0.825
+r_i = 0.55
+
+meshball = uw.meshing.AnnulusInternalBoundary(radiusOuter=r_o, 
+                                              radiusInternal=r_int, 
+                                              radiusInner=r_i, 
+                                              cellSize_Inner=res,
+                                              cellSize_Internal=res*0.5,
+                                              cellSize_Outer=res,
+                                              centre=False,)
+
+
+meshbox = uw.meshing.UnstructuredSimplexBox(
+    cellSize=res,
+    minCoords=(-1.0, -1.0),
+    maxCoords=(+1.0, +1.0),
+    qdegree=3,
+)
+
+mesh = meshball
+
+x, y = mesh.CoordinateSystem.X
+r, th = mesh.CoordinateSystem.xR
+unit_rvec = mesh.CoordinateSystem.unit_e_0
+
+
+# Orientation of surface normals
+Gamma_N = mesh.Gamma
+
+
+# %%
+rank_var = uw.discretisation.MeshVariable(
+    varname="Rank",
+    mesh=mesh, 
+    vtype = uw.VarType.SCALAR,
+    varsymbol=r"Ra",
+    continuous=False,
+    degree=0,
+)
+
+
+scalar_var = uw.discretisation.MeshVariable(
+    varname="Radius",
+    mesh=mesh, 
+    vtype = uw.VarType.SCALAR,
+    varsymbol=r"r"
+)
+
+scalar_var.array[...] = uw.function.evaluate(r, scalar_var.coords)
+rank_var.array[...] = uw.mpi.rank
+
+# %%
+## Functions that we can evaluate anywhere
+scalar_fn = scalar_var.sym[0,0] * sympy.sin(x) * sympy.cos(x)
+
+# %%
+# Random points in the bounding box of the domain (on each processor)
+coords = 1 - 2 * np.random.random(size=(1000,2)) 
+
+# %%
+results, is_extrapolated =  uw.function.global_evaluate(scalar_fn, coords, rbf=False, check_extrapolated=True)
+
+# %%
+rank, is_extrapolated =  uw.function.global_evaluate(rank_var.sym[0], coords, rbf=False, check_extrapolated=True)
+
+# %%
+true_results = uw.function.global_evaluate(r * sympy.sin(x) * sympy.cos(x)  , coords, rbf=False, check_extrapolated=False)
+
+# %%
+inside_domain = mesh.points_in_domain(coords)
+
+print(f"{uw.mpi.rank} Max difference              - {(results - true_results).max()}", flush=True)
+print(f"{uw.mpi.rank} Mean difference             - {(results - true_results).mean()}", flush=True)
+
+print(f"{uw.mpi.rank} Max difference  [in domain] - {(results - true_results)[inside_domain].max()}", flush=True)
+print(f"{uw.mpi.rank} Mean difference [in domain] - {(results - true_results)[inside_domain].mean()}", flush=True)
+
+# %%
+swarm = uw.swarm.Swarm(mesh)
+error_var = uw.swarm.SwarmVariable("Error", swarm, vtype=uw.VarType.SCALAR, _proxy=False)
+s_rank_var = uw.swarm.SwarmVariable("FRank", swarm, vtype=uw.VarType.SCALAR, _proxy=False)
+
+# %%
+swarm.add_particles_with_global_coordinates(coords, migrate=False)
+
+# %%
+error_var.array[...] = (results - true_results)
+s_rank_var.array[...] = rank[...]
+uw.mpi.barrier()
+
+# %%
+mesh.write_timestep(
+    "Parallel_Evaluate",
+    meshUpdates=True,
+    meshVars=[scalar_var, rank_var],
+    outputPath="output",
+    index=uw.mpi.size,
+
+)
+
+# %%
+swarm.write_timestep(
+    "Parallel_Evaluate",
+    "data_swarm",
+    swarmVars=[error_var, s_rank_var],
+    outputPath="output",
+    index=uw.mpi.size,
+    force_sequential=True,
+)
+
+
+# %%
+if uw.mpi.rank == 0:
+    print("Complete", flush=True)
+
+# %%
+exit(0)
+
+# %%