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 FEniCS/SciPy:
import matplotlib.pyplot as plt import numpy as np from dolfin import ( EigenMatrix, FacetNormal, FunctionSpace, TestFunction, TrialFunction, UnitIntervalMesh, UnitSquareMesh, assemble, dot, ds, dx, grad, ) def _assemble_eigen(form): L = EigenMatrix() assemble(form, tensor=L) return L k = 20 # mesh = UnitIntervalMesh(n) mesh = UnitSquareMesh(k, k) degree = 1 V = FunctionSpace(mesh, "CG", degree) u = TrialFunction(V) v = TestFunction(V) mesh = V.mesh() n = FacetNormal(mesh) A = _assemble_eigen(dot(grad(u), grad(v)) * dx - dot(n, grad(u)) * v * ds).sparray() out = np.linalg.eigvals(A.toarray()) plt.plot(out.real, out.imag, "o") plt.show()