evs.py

from nutils import mesh, function, plot, _
import numpy as np
import scipy as sp


nelems = 40
xs = np.linspace(0, 1, nelems + 1)
domain, geom = mesh.rectilinear([xs])
basis = domain.basis('spline', degree=2)

# Eigenvalue decomposition and formation of reduced system
integrand = function.outer(basis.grad(geom)).sum(-1)
matrix = domain.integrate(integrand, geometry=geom, ischeme='gauss4')

cons = domain.boundary.project(1, onto=basis, geometry=geom, ischeme='gauss4')
I, = np.where(np.isnan(cons))

mx = matrix.toarray()
mx = mx[np.ix_(I, I)]
w, v = np.linalg.eigh(mx)
vv = np.zeros((len(basis), v.shape[1]))
vv[I,:] = v

newbasis = (basis[:,_] * vv).sum(0)

# # Plot reduced basis functions
# for i in range(len(I)):
#     dd = np.zeros((len(newbasis),))
#     dd[i] = 1
#     bfun = newbasis.dot(dd)
#     points, colors = domain.elem_eval([geom, bfun], ischeme='bezier9', separate=True)
#     with plot.PyPlot('bfun', index=i) as plt:
#         plt.mesh(points, colors)
#         plt.ylim(-1, 1)
#         plt.colorbar()

# Solve with reduced system
integrand = function.outer(newbasis.grad(geom)).sum(-1)
matrix = domain.integrate(integrand, geometry=geom, ischeme='gauss4')
rhs = domain.integrate(newbasis, geometry=geom, ischeme='gauss4')
lhs = matrix.solve(rhs)
red_sol = newbasis.dot(lhs)

# Solve with original system
integrand = function.outer(basis.grad(geom)).sum(-1)
matrix = domain.integrate(integrand, geometry=geom, ischeme='gauss4')
rhs = domain.integrate(basis, geometry=geom, ischeme='gauss4')
cons = domain.boundary.project(0, onto=basis, geometry=geom, ischeme='gauss4')
lhs = matrix.solve(rhs, constrain=cons)
orig_sol = basis.dot(lhs)

points, red, orig = domain.elem_eval([geom, red_sol, orig_sol], ischeme='bezier9', separate=True)
with plot.PyPlot('redsol', index=0) as plt:
    plt.mesh(points, red)
    plt.colorbar()
with plot.PyPlot('origsol', index=0) as plt:
    plt.mesh(points, orig)
    plt.colorbar()