diff --git a/applications/solvers/basic/laplacianFoam/Make/files b/applications/solvers/basic/laplacianFoam/Make/files
index b1027ce8fd..290d5d6e92 100644
--- a/applications/solvers/basic/laplacianFoam/Make/files
+++ b/applications/solvers/basic/laplacianFoam/Make/files
@@ -1,3 +1,5 @@
+/* lduPrimitiveMesh.C */
+
laplacianFoam.C
EXE = $(FOAM_APPBIN)/laplacianFoam
diff --git a/applications/solvers/basic/laplacianFoam/lduPrimitiveMesh-hack.C b/applications/solvers/basic/laplacianFoam/lduPrimitiveMesh-hack.C
new file mode 100644
index 0000000000..4de262b5d0
--- /dev/null
+++ b/applications/solvers/basic/laplacianFoam/lduPrimitiveMesh-hack.C
@@ -0,0 +1,1414 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | www.openfoam.com
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+ Copyright (C) 2013-2017 OpenFOAM Foundation
+ Copyright (C) 2019,2022 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+ This file is part of OpenFOAM.
+
+ OpenFOAM is free software: you can redistribute it and/or modify it
+ under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+ ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+ FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
+ for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with OpenFOAM. If not, see .
+
+\*---------------------------------------------------------------------------*/
+
+#include "lduPrimitiveMesh.H"
+#include "processorLduInterface.H"
+#include "edgeHashes.H"
+#include "labelPair.H"
+#include "processorGAMGInterface.H"
+#include "cyclicAMIGAMGInterface.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+ defineTypeNameAndDebug(lduPrimitiveMesh, 0);
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
+
+void Foam::lduPrimitiveMesh::checkUpperTriangular
+(
+ const label size,
+ const labelUList& l,
+ const labelUList& u
+)
+{
+ forAll(l, facei)
+ {
+ if (u[facei] < l[facei])
+ {
+ FatalErrorInFunction
+ << "Reversed face. Problem at face " << facei
+ << " l:" << l[facei] << " u:" << u[facei]
+ << abort(FatalError);
+ }
+ if (l[facei] < 0 || u[facei] < 0 || u[facei] >= size)
+ {
+ FatalErrorInFunction
+ << "Illegal cell label. Problem at face " << facei
+ << " l:" << l[facei] << " u:" << u[facei]
+ << abort(FatalError);
+ }
+ }
+
+ for (label facei=1; facei < l.size(); ++facei)
+ {
+ if (l[facei-1] > l[facei])
+ {
+ FatalErrorInFunction
+ << "Lower not in incremental cell order."
+ << " Problem at face " << facei
+ << " l:" << l[facei] << " u:" << u[facei]
+ << " previous l:" << l[facei-1]
+ << abort(FatalError);
+ }
+ else if (l[facei-1] == l[facei])
+ {
+ // Same cell.
+ if (u[facei-1] > u[facei])
+ {
+ FatalErrorInFunction
+ << "Upper not in incremental cell order."
+ << " Problem at face " << facei
+ << " l:" << l[facei] << " u:" << u[facei]
+ << " previous u:" << u[facei-1]
+ << abort(FatalError);
+ }
+ }
+ }
+}
+
+
+Foam::labelListList Foam::lduPrimitiveMesh::globalCellCells
+(
+ const lduMesh& mesh,
+ const globalIndex& globalNumbering
+)
+{
+ const lduAddressing& addr = mesh.lduAddr();
+ lduInterfacePtrsList interfaces = mesh.interfaces();
+
+ const labelList globalIndices
+ (
+ identity
+ (
+ addr.size(),
+ globalNumbering.localStart(UPstream::myProcNo(mesh.comm()))
+ )
+ );
+
+ // Get the interface cells
+ PtrList nbrGlobalCells(interfaces.size());
+ {
+ // Initialise transfer of restrict addressing on the interface
+ forAll(interfaces, inti)
+ {
+ if (interfaces.set(inti))
+ {
+ interfaces[inti].initInternalFieldTransfer
+ (
+ Pstream::commsTypes::nonBlocking,
+ globalIndices
+ );
+ }
+ }
+
+ // Wait for all
+ UPstream::waitRequests();
+
+ forAll(interfaces, inti)
+ {
+ if (interfaces.set(inti))
+ {
+ nbrGlobalCells.set
+ (
+ inti,
+ new labelList
+ (
+ interfaces[inti].internalFieldTransfer
+ (
+ Pstream::commsTypes::nonBlocking,
+ globalIndices
+ )
+ )
+ );
+ }
+ }
+ }
+
+
+ // Scan the neighbour list to find out how many times the cell
+ // appears as a neighbour of the face. Done this way to avoid guessing
+ // and resizing list
+ labelList nNbrs(addr.size(), Zero);
+
+ const labelUList& nbr = addr.upperAddr();
+ const labelUList& own = addr.lowerAddr();
+
+ {
+ forAll(nbr, facei)
+ {
+ nNbrs[nbr[facei]]++;
+ nNbrs[own[facei]]++;
+ }
+
+ forAll(interfaces, inti)
+ {
+ if (interfaces.set(inti))
+ {
+ for (const label celli : interfaces[inti].faceCells())
+ {
+ nNbrs[celli]++;
+ }
+ }
+ }
+ }
+
+
+ // Create cell-cells addressing
+ labelListList cellCells(addr.size());
+
+ forAll(cellCells, celli)
+ {
+ cellCells[celli].setSize(nNbrs[celli], -1);
+ }
+
+ // Reset the list of number of neighbours to zero
+ nNbrs = 0;
+
+ // Scatter the neighbour faces
+ forAll(nbr, facei)
+ {
+ const label c0 = own[facei];
+ const label c1 = nbr[facei];
+
+ cellCells[c0][nNbrs[c0]++] = globalIndices[c1];
+ cellCells[c1][nNbrs[c1]++] = globalIndices[c0];
+ }
+ forAll(interfaces, inti)
+ {
+ if (interfaces.set(inti))
+ {
+ const labelUList& faceCells = interfaces[inti].faceCells();
+
+ forAll(faceCells, facei)
+ {
+ const label c0 = faceCells[facei];
+ cellCells[c0][nNbrs[c0]++] = nbrGlobalCells[inti][facei];
+ }
+ }
+ }
+
+ return cellCells;
+}
+
+
+Foam::label Foam::lduPrimitiveMesh::totalSize
+(
+ const PtrList& meshes
+)
+{
+ label size = 0;
+
+ for (const lduPrimitiveMesh& msh : meshes)
+ {
+ size += msh.lduAddr().size();
+ }
+ return size;
+}
+
+
+Foam::labelList Foam::lduPrimitiveMesh::upperTriOrder
+(
+ const label nCells,
+ const labelUList& lower,
+ const labelUList& upper
+)
+{
+ labelList nNbrs(nCells, Zero);
+
+ // Count number of upper neighbours
+ forAll(lower, facei)
+ {
+ if (upper[facei] < lower[facei])
+ {
+ FatalErrorInFunction
+ << "Problem at face:" << facei
+ << " lower:" << lower[facei]
+ << " upper:" << upper[facei]
+ << exit(FatalError);
+ }
+ nNbrs[lower[facei]]++;
+ }
+
+ // Construct cell-upper cell addressing
+ labelList offsets(nCells+1);
+ offsets[0] = 0;
+ forAll(nNbrs, celli)
+ {
+ offsets[celli+1] = offsets[celli]+nNbrs[celli];
+ }
+
+ nNbrs = offsets;
+
+ labelList cellToFaces(offsets.last());
+ forAll(upper, facei)
+ {
+ const label celli = lower[facei];
+ cellToFaces[nNbrs[celli]++] = facei;
+ }
+
+ // Sort
+
+ labelList oldToNew(lower.size());
+
+ DynamicList