ninepointsvd.py

#!/usr/bin/env python3                                                                                                                         
from numpy.linalg import pinv
from numpy import diag
from numpy import transpose
from numpy import dot
import numpy as np

# a 3x3 stencil with centroids between (-1.5,1) and (0.5,1)
# to interp onto face at (0,0)
# upwind cell first, then downwind cell, then top row L-R, bottom row L-R, then left-most cell in middle row                                   
X = [ \
    [-0.5, 0], \
    [0.5, 0], \
    [-1.5, 1], \
    [-0.5, 1], \
    [0.5, 1], \
    [-1.5, -1], \
    [-0.5, -1], \
    [0.5, -1], \
    [-1.5, 0] \
]

# unweighted matrix with one row per cell and one column per polynomial term
# in local coordinate system
Btwiddle = [ \
    [1, X[0][0], X[0][1], X[0][0]**2, X[0][0]*X[0][1], X[0][1]**2, X[0][0]**3, X[0][0]**2*X[0][1], X[0][0]*X[0][1]**2], \
    [1, X[1][0], X[1][1], X[1][0]**2, X[1][0]*X[1][1], X[1][1]**2, X[1][0]**3, X[1][0]**2*X[1][1], X[1][0]*X[1][1]**2], \
    [1, X[2][0], X[2][1], X[2][0]**2, X[2][0]*X[2][1], X[2][1]**2, X[2][0]**3, X[2][0]**2*X[2][1], X[2][0]*X[2][1]**2], \
    [1, X[3][0], X[3][1], X[3][0]**2, X[3][0]*X[3][1], X[3][1]**2, X[3][0]**3, X[3][0]**2*X[3][1], X[3][0]*X[3][1]**2], \
    [1, X[4][0], X[4][1], X[4][0]**2, X[4][0]*X[4][1], X[4][1]**2, X[4][0]**3, X[4][0]**2*X[4][1], X[4][0]*X[4][1]**2], \
    [1, X[5][0], X[5][1], X[5][0]**2, X[5][0]*X[5][1], X[5][1]**2, X[5][0]**3, X[5][0]**2*X[5][1], X[5][0]*X[5][1]**2], \
    [1, X[6][0], X[6][1], X[6][0]**2, X[6][0]*X[6][1], X[6][1]**2, X[6][0]**3, X[6][0]**2*X[6][1], X[6][0]*X[6][1]**2], \
    [1, X[7][0], X[7][1], X[7][0]**2, X[7][0]*X[7][1], X[7][1]**2, X[7][0]**3, X[7][0]**2*X[7][1], X[7][0]*X[7][1]**2], \
    [1, X[8][0], X[8][1], X[8][0]**2, X[8][0]*X[8][1], X[8][1]**2, X[8][0]**3, X[8][0]**2*X[8][1], X[8][0]*X[8][1]**2] \
]

wp = [10, 10, 1, 1, 1, 1, 1, 1, 1]
wc = wp
WP = diag(wp)
WC = diag(wc)

print("unweighted ", np.array(Btwiddle))
B = dot(dot(WC, Btwiddle), WP)
print("weighted ", np.array(B))

Binv = pinv(B)
c = dot(dot(Binv[0], wc), wp)
print("coeffsi ", c)
print("Sum to ", sum(c))