This repository provides a Python implementation of the Meshfree Geometric Multilevel (MGM) Method, originally developed by Dr. Grady Wright, Professor at Boise State University, and his collaborators. The MGM method is designed for solving elliptic equations on surfaces, particularly suited for problems arising in computational science and engineering.
The MGM method offers a powerful approach to solve and precondition linear systems resulting from meshfree discretizations of elliptic equations. It encompasses algorithms for both 2D surfaces embedded in 3D space and 2D/3D Euclidean domains with boundaries. This Python implementation extends the capabilities of the original MGM method for 2D surfaces, allowing users to leverage its functionality within the Python ecosystem. For quick reference on the current features and limitations of pyMGM, please see the Quick Reference section below.
The code makes use of the weighted sample elmination (WSE) method from the cyCodeBase Package. This method is implemented in C++ and the mgm
package includes pybind interface files in the +util
directory to access this method. However, before you can use pyMGM, the shared object (so) files must
be compiled on your machine. Below are the instructions to build these files, which only has to be done once after cloning the repository.
Using terminal enter the following command to compile the shared object files and link to the python interface files:
g++ -O3 -Wall -shared -std=c++11 -fPIC `python3 -m pybind11 --includes` src/pymgm_test/utils/bindings.cpp src/pymgm_test/utils/PcCoarsen2D.cpp -o src/pymgm_test/utils/PcCoarsen`python3-config --extension-suffix`
This should create a .so file in the src/pymgm_test/utils/
directory.
Note that currently the files from the cyCodeBase Package are included in this repository. However, once this repository is made public, users will be required to clone the cyCodeBase Package to their machine. Depending on the location where this package is cloned, the build instructions above may change.
To use the pyMGM library you will need to add the top level folder to the Python path. This can be done as follows:
export PYTHONPATH=$PYTHONPATH:$(pwd)
. This will add the top level folder to the Python path. set PYTHONPATH=%PYTHONPATH%;%cd%
MGM was initially developed as a multilevel method for solving/preconditioning linear systems that arise for meshfree discretizations of elliptic equations on 2D surfaces embedded in 3D space as described in [1]. However, we have found that it also appears to work for other discretizations and for Euclidean domains with boundaries. In this sense, the method can be used as a black-box solver for systems $L_h u_h = f_h$ corresponding to the discretization of an elliptic PDE over a domain $\Omega$.
The inputs for MGM are simple:
Lh
Lh
, stored as a $N\times d$ MATLAB matrix X
, where $d$ is the dimension of (embedding) space. For nodal FE discretizations, $X$ would be the quadrature nodes over the elements of the mesh (without repeats). domainVol
With these three components, one first sets up the MGM method, which consists of computing the point sets and discrete operators on the various coarser levels
After importing the mgm2D class from src.pymgm_test.mgm2d.mgm2d within the project an instance of the class can be created by calling:
# For a 2D Euclidean domain:
mgm = mgm2D(Lh,X,domainVol, False, 1)
Note:
Lh
.The following code is then used to solve the system for a given right hand side fh
using the default acceleration method of 'none'
uh,flag,relres,iters,resvec = mgmobj.solve(mgm.obj,fh,tol,'none');
Note:
uh,flag,relres,iters,resvec = mgmobj.solve(mgm.obj,fh,tol,'none', True, 100);
See the help text of mgm.solve
for a description of the outputs and other optional input arguments.
For PDEs such as the surface Poisson problem or Poisson problems on Euclidean domains with pure Neumann boundary conditions the matrix $L_h$ will have a one dimensional nullspace corresponding to constant vectors. In these cases where there is a one-dimensional nullspace, it is necessary to alter the above construction procedure. For example, for the surface problem one would instead do:
% For a 2D surface embedded in 3D space where Lh has a one dimensional nullspace (e.g. surface Poisson problem)
mgmobj = mgm2D(Lh,X,domainVol,True, 1)
Again see the help text of mgm2ds
, mgm2d
, or mgm3d
for a description of other optional input arguments.
util.squarepoisson
The util.squarepoisson
function generates a square poisson distribution problem with matrices, points, right hand sides, and exact solutions already computed
from a given n for an $N\times N$ matrix. For the current implementation of pyMGM n ≥ 50 is recommended.
Below is an example using this with pyMGM
Lh, x, vol, fh, uexact = squarepoissond(50)
# Set-up MGM2D for this problem
mgm = mgm2D(Lh,x,vol,False,1)
# Solve the linear system using pyMGM with default acceleration of 'none'
uh,flag,relres,iters,resvec = mgm.solve(mgm.obj,fh,1e-10,'none',100)
To optionally create a plot of the solution, the following code can be used:
import matplotlib.pyplot as plt
# Scatter plot of the solution
plt.figure()
plt.scatter(x[:, 0], x[:, 1], c=uh, cmap='viridis', s=20)
plt.colorbar(label='Solution')
plt.title('Scatter plot of the solution')
plt.xlabel('x')
plt.ylabel('y')
plt.show()
hasConstNullSpace
parameter for a value of True is not yet implemented. The MGM method is designed to handle systems with a one-dimensional nullspace, but this feature is not yet implemented in the current version of pyMGM.[1] G. B. Wright, A. M. Jones, and V. Shankar . MGM: A meshfree geometric multilevel method for systems arising from elliptic equations on point cloud surfaces. SIAM J. Sci. Comput., 45, A312-A337 (2023) (arxiv, journal)
This website/app was created for a
Boise State University
Computer Science Senior Design Project by
Sean Calkins
Tanner Frost
Paul Vanderveen
For information about sponsoring a project go to
https://www.boisestate.edu/coen-cs/community/cs481-senior-design-project/