I'm looking at solving systems with the FEM discretization $$ -\int_\Omega (\Delta u) v = \int_\Omega \nabla u \cdot \nabla v - \int_\Gamma (n\cdot\nabla u) v. $$ without applying Dirichlet- or Neumann-type boundary conditions. The resulting matrix is generally not self-adjoint (except in for 1D meshes).
The kernel consists of all linear functions over the domain, so there are always $n+1$ eigenvalues 0 (where $n$ is the dimensionality of the mesh).
The right-hand side is such that the system is consistent and I would like to find a solution of the system. The idea is to use GMRES, starting off with an all-zero initial guess.
Is there a good preconditioner for this known in literature? I played around with Dirichlet- and Neumann-Laplace without much success.
To get a feel, here's how to create the matrix and plot the spectrum in sckit-fem:
import matplotlib.pyplot as plt import meshzoo import numpy as np import skfem from skfem.helpers import dot, grad @skfem.BilinearForm def laplace(u, v, _): return dot(grad(u), grad(v)) @skfem.BilinearForm def flux(u, v, w): return dot(w.n, u.grad) * v points, cells = meshzoo.disk(6, 20) pT = np.ascontiguousarray(points.T) cT = np.ascontiguousarray(cells.T) mesh = skfem.MeshTri(pT, cT) element = skfem.ElementTriP1() basis = skfem.CellBasis(mesh, element) facet_basis = skfem.FacetBasis(basis.mesh, basis.elem) lap = skfem.asm(laplace, basis) boundary_terms = skfem.asm(flux, facet_basis) A = lap - boundary_terms out = np.linalg.eigvals(A.toarray()) plt.plot(out.real, out.imag, "o") plt.savefig("out.png", bbox_inches="tight") plt.show()