DDM-Matlab

Overview

This repository contains MATLAB prototypes of several domain decomposition methods (DDM) for the 2D Poisson equation on triangular meshes produced with Gmsh. Each solver assembles local stiffness matrices for every partition, eliminates interior degrees of freedom to form Schur complements on the interface, and either solves the resulting system directly or with preconditioned conjugate gradients (PCGM). The code is organised by algorithm (direct Schur complement, Neumann preconditioned PCGM, lumped preconditioning, two-level additive Schwarz, and FETI-DP) and shares a common set of mesh-parsing and algebra utilities. Reference solutions are compared with the manufactured solution u(x,y) = sin(pi*x) * sin(2*pi*y) so that every script reports the maximum-norm error and plots the solution, reference solution, and error surface.

Prerequisites

• MATLAB R2018a or newer. The scripts rely only on base MATLAB functions plus pdesurf for plotting (from PDE Toolbox).
• Gmsh 3.x/4.x for generating .geo and .msh files. The sample geometry is in square.geo.
• A mesh partitioner that can inject partition tags into the .msh file. The repository ships the GPL-licensed puffin-0.1.6 source, which can partition a Gmsh mesh into equal-sized subdomains and emit the extended .msh files that the ExtractGmshPart_puff parser expects.
• (Optional) A PDF viewer if you want to build the Puffin manual in puffin-0.1.6/doc/manual.

After cloning the repository, add it to your MATLAB path once per session:

addpath(genpath('/path/to/ddm-matlab'));

Repository Layout

• DirectDDM/ - One-level Schur complement solver that eliminates interior unknowns and solves the reduced interface system directly.
• Neumann/ - PCGM solver with a Neumann-Neumann preconditioner assembled from subdomain solves.
• LumpedPreconditoner/ - PCGM solver that uses a lumped interface preconditioner (diagonal of the Schur complement scaled by the partition-of-unity).
• DDM-TwoLevel/ - Two-level additive Schwarz method that augments the lumped preconditioner with a coarse problem built from corner and remaining interface nodes.
• FETI-DP/ - Dual-primal FETI formulation that enforces continuity with Lagrange multipliers and coarse primal variables at corners and edges.
• ExtractGmsh/ - Standalone scripts to parse standard and partitioned .msh files, dump intermediate point/edge/triangle tables, and debug interface detection.
• sparsity/ - Copy of the two-level solver with AssembleForSparse.m for inspecting the sparsity of the assembled global matrix.
• puffin-0.1.6/ and puffin-0.1.6.zip - Third-party partitioner from Johan Hoffman and Anders Logg (GPLv2). Build it if you need to generate new partitioned Gmsh files.
• ExtractGmshPart_puff.m, GetInterfaceNodes*.m, GetRestrictionMatrix.m, GetScalingMatrix*.m, AssembleMatrix.m, AssembleVector.m, and Poisson.m are replicated inside each solver directory so that the scripts are self-contained when launched from MATLAB.

Most solver folders also include pre-generated meshes such as squarepart4.msh, square2_part4.msh, or mymeshpart10.msh, and the local stiffness/load files A_i.mat / b_i.mat that are regenerated every time you run a script.

Typical Workflow

1. Generate or edit a geometry. Use the provided square.geo as a starting point or author a new .geo file in Gmsh. Export a .msh file with first-order triangles.
2. Partition the mesh. Run Puffin (or another tool) to split the mesh into the desired number of subdomains and write a .msh file that contains, for every triangle, the number of adjacent parts and their IDs. The examples squarepart{4,10,20,...}.msh and square2_part4.msh illustrate the expected format.
3. Place the mesh where the solver expects it. Copy the partitioned .msh file into the solver directory you plan to run and update the filename passed to ExtractGmshPart_puff. Each solver currently references either squarepart4.msh, square2_part4.msh, or mymeshpart4.msh.
4. Launch MATLAB and run the solver script. For example:

   cd('DDM-TwoLevel');
   PoissonSolver_DDM_TwoLevel;

   Every script starts with delete *.mat to remove stale A_i.mat/b_i.mat files, extracts the mesh, assembles local matrices, and executes the solver. Progress is echoed in the Command Window (number of PCGM iterations, residual norms, and the maximum error).
5. Inspect the results. The scripts call pdesurf three times (solution, exact solution, error) and, for iterative solvers, plot residual vs iteration (e.g., Lump.mat, NN_m.mat, 2Level.mat store residual histories).

Solver Directories

DirectDDM

• Main entry point: PoissonSolver_Direct.m.
• What it does: Forms the Schur complement of each subdomain, assembles the global interface system with GetRestrictionMatrix, and solves it directly. It then reconstructs interior unknowns by back substitution and compares the field to the manufactured solution.
• Key helpers: ExtractGmshPart_puff.m, AssembleMatrix.m, AssembleVector.m, GetInterfaceNodes*.m.
• Sample meshes: source/square2_part4.msh, source/square_part8_remote.msh, mac/square2_part8.msh.

Neumann

• Main entry point: PoissonSolver_PCGM_NN.m.
• Algorithm: Preconditioned conjugate gradients on the interface unknowns using a Neumann-Neumann preconditioner assembled in GetPreconditionedNeumann.m. Interface operators are applied through repeated local solves in GetMatrixVectorProduct.m. Scaling matrices D enforce a partition-of-unity on shared nodes.
• Outputs: Stores residual histories in NN_m.mat.

LumpedPreconditoner

• Main entry point: PoissonSolver_PCGM_Lumped.m.
• Algorithm: Identical PCGM loop but the preconditioner uses the diagonal (lumped) interface stiffness scaled by D through GetPreconditionedZLumped.m.
• Outputs: Residual curves saved in Lump.mat.

DDM-TwoLevel

• Main entry point: PoissonSolver_DDM_TwoLevel.m.
• Algorithm: Builds a coarse problem on interface corner nodes (node_corner) and remaining interface nodes (node_remaining). GetCoarseGrid.m, GetRestrictionMatrix.m, and GetMatrixVectorProductTwoLevel.m assemble the coarse operators, while GetPreconditionedTwolevelNeumann.m adds the coarse-space correction. The inner coarse solves reuse DoPCGMLumpedTwoLevel.m.
• Meshes: Includes square.msh, squarepart{2,4,...,80}.msh, mymesh.msh, and mymeshpart*.msh test cases.
• Outputs: Residual curves saved in 2Level.mat and intermediate data for sparsity studies in IterVsNS/.

FETI-DP

• Main entry point: PoissonSolver_FETI_DP.m.
• Algorithm: Implements a dual-primal FETI method. Coarse primal unknowns live on corner/edge nodes via Bc, while continuity on the remaining interface is enforced by the jump operator Br and scaling Dr. Lagrange multipliers (lambda) are iteratively updated with PCGM, using GetMatrixVectorProductFETIDP.m for the action of the dual Schur complement and GetPreconditioner_FETIDP.m for diagonal scaling. GetResidue_FETIDP.m, GetGlobalCornerNode.m, and GetGlobalRemainingNode.m reconstruct primal solutions once the dual problem converges.
• Outputs: Saves residual curves in 2Level.mat (shared naming) and prints the maximum error.

Mesh and Interface Extraction Helpers

All solvers invoke ExtractGmshPart_puff.m, which parses a partitioned .msh file and returns:

• p_f, e_f, t_f: point, edge, and triangle tables per partition ready for AssembleMatrix and AssembleVector.
• GLnode: mapping from renamed local node indices back to the global numbering.
• npart: number of partitions in the mesh.
• point_exact, triangle_exact: domain-wide coordinates for plotting with pdesurf.
• node_gamma, node_interior: renamed interface and interior nodes per partition.
• node_interfacewhole, node_remainingwhole: full lists (with original numbering) used to scatter interface solutions back to U.
• R: restriction matrices that inject interface unknowns into the subdomain order.
• D: partition-of-unity scaling derived by GetScalingMatrix.m.
• node_corner, node_remaining, node_c, node_r: coarse-grid nodes (original and renamed).
• Rc, Rr, Bc: restriction operators for corner-only, remaining-interface-only, and global corner nodes respectively.
• node_*_OG: convenience copies that retain original Gmsh numbering for debugging and plotting.

FETI-DP additionally needs Br (jump operator built by GetJumpOperator.m or GetInterfaceNodes_puff.m), Rrr (restriction from global interface to remaining nodes), and Dr (interface scaling for the dual problem). These are exposed by the extended signature of ExtractGmshPart_puff.m inside FETI-DP/.

The ExtractGmsh folder contains alternative parsers (ExtractGmsh.m, ExtractGmsh_part.m, ExtractGmsh_PartInterface.m) that dump mesh subsets to text files (points.txt, edges.txt, etc.) to help debug new meshes or partition assignments.

Linear Algebra Helpers

• AssembleMatrix.m and AssembleVector.m assemble the finite-element stiffness matrix and load vector for Poisson's equation by quadrature over each triangle and boundary edge (sourced from Hoffman & Logg, GPLv2).
• Poisson.m specifies the bilinear and linear forms, boundary conditions, and the manufactured forcing f(x,y) = 5*pi^2*sin(pi*x)*sin(2*pi*y). Modify f, gd, or gn to solve a different problem.
• GetMatrixVectorProduct*.m files evaluate Schur complement actions by solving interior blocks and combining with interface data.
• GetPreconditioned*.m files assemble Neumann, lumped, two-level, and FETI-DP preconditioners by inverting (or partially inverting) local Schur complements.
• DoPCGMLumped*.m implement the PCGM loop used both as an outer solver (Neumann/Lumped) and as an inner coarse-grid solver (TwoLevel/FETI).
• AssembleForSparse.m (in DDM-TwoLevel/ and sparsity/) is a diagnostic utility that reconstructs the global stiffness pattern for plotting sparsity.

Understanding these helpers is useful when adapting the code to a different PDE or when adding diagnostics such as energy norms, iteration logs, or stopping criteria.

Puffin Partitioner

The puffin-0.1.6 directory contains the source, scripts, and documentation for Puffin (a light-weight partitioner originally bundled with the FEniCS project). Typical usage:

1. Build Puffin following the instructions in puffin-0.1.6/README.
2. Run Puffin on a .msh file to produce squarepart4.msh, squarepart10.msh, etc. The resulting mesh has extended tags describing how many partitions share each element and lists the adjacent part IDs. Those tags are what ExtractGmshPart_puff.m consumes.
3. Copy the new .msh file into the solver directory and point ExtractGmshPart_puff to it.

You can keep the zipped archive (puffin-0.1.6.zip) for reproducibility.

Saved Data and Experiment Scripts

• A_i.mat / b_i.mat - Local stiffness and load matrices per subdomain written during every solver run. They are re-used inside GetMatrixVectorProduct* and the PCGM loops. The scripts delete them on startup to avoid mixing meshes.
• Lump.mat, NN_m.mat, 2Level.mat - Residual histories (iter, nres) saved by the respective solvers. These files drive the "residual vs iteration" plots.
• IterVsNS/ - Placeholder directories used to store parametric studies (iterations vs number of subdomains). Drop your experiment logs here to keep the solver folders tidy.

Customisation Tips

• Change the PDE or forcing term: edit Poisson.m (inside the solver you are running) and update f, gd, gn, or the penalty g. Keep the function signature unchanged because it is called by AssembleMatrix/AssembleVector.
• Use a different mesh: point ExtractGmshPart_puff to your .msh file. The parser assumes linear triangles and that Gmsh stores elements in the standard order (edges before triangles). Verify that the partition tags are present; otherwise node_gamma will be empty.
• Vary tolerances: the Neumann and Lumped solvers stop when norm(r_new)/norm(r_0) drops below 1e-5. The two-level and FETI solvers accept a tolc or use 1e-6 by default - modify the value defined near the top of those scripts.
• Reuse assemblies: comment out delete *.mat if you want to inspect A_i.mat/b_i.mat across runs, but be mindful that changing the mesh without deleting the files will corrupt subsequent runs.
• Plot sparsity: run AssembleForSparse.m (from sparsity/) after building A_i.mat to capture the global stiffness pattern for debugging or publications.
• Batch execution: launch MATLAB with -nodesktop/-nodisplay and run PoissonSolver_*; exit to execute solvers from a shell script. This is how parameter sweeps in IterVsNS were produced.

Troubleshooting

• manpath warnings on macOS are benign and originate from the shell environment when running MATLAB from the terminal.
• If ExtractGmshPart_puff reports missing interface nodes, double-check that the mesh contains partition IDs for every element and that interior nodes were not duplicated when editing the .msh.
• Singular local Schur complements typically indicate a disconnected partition or a mesh with hanging nodes. Regenerate the partition or slightly perturb the geometry.
• The plots require MATLAB's PDE Toolbox. If pdesurf is unavailable, replace the plotting section with trisurf or patch.

Licensing and Credits

• The finite element assembly routines (AssembleMatrix.m, AssembleVector.m, Poisson.m) and the Puffin partitioner are copyrighted by Johan Hoffman and Anders Logg and distributed under the GNU GPL v2. This repository preserves that license.
• The domain decomposition scripts, mesh utilities, and experiments were authored by Sudhi Sharma P V for research and teaching purposes. Cite accordingly if you use the code in publications.

Please reach out via the repository issue tracker with questions or improvements.