"""Solver B: 5-point finite differences on the square grid (no mesh, no FEM code shared)."""
import numpy as np, scipy.sparse as sp, scipy.sparse.linalg as spla
from geom import inside, on_boundary
def lap(P,N):
    h=1.0/N; M=3*N+1; xs=np.arange(M)*h
    X,Y=np.meshgrid(xs,xs,indexing='ij'); pts=np.c_[X.ravel(),Y.ravel()]
    inn=inside(P,pts)&~on_boundary(P,pts)
    idx=-np.ones(M*M,int); ids=np.where(inn)[0]; idx[ids]=np.arange(len(ids))
    I,J=np.divmod(ids,M); rows=[np.arange(len(ids))]; cols=[np.arange(len(ids))]; vals=[np.full(len(ids),4.0)]
    for di,dj in ((1,0),(-1,0),(0,1),(0,-1)):
        nb=idx[(I+di)*M+(J+dj)]; m=nb>=0
        rows.append(np.where(m)[0]); cols.append(nb[m]); vals.append(-np.ones(m.sum()))
    A=sp.csr_matrix((np.concatenate(vals),(np.concatenate(rows),np.concatenate(cols))),(len(ids),)*2)/h**2
    return A,pts[ids]
def solve(P,N,k=50,vecs=False):
    A,p=lap(P,N)
    w,V=spla.eigsh(A,k,sigma=0,which='LM'); o=np.argsort(w)
    return (w[o],V[:,o],p) if vecs else w[o]