pub fn quad4_shape(xi: f64, eta: f64) -> [f64; 4]
N_i(xi, eta) = (1 + xi_i xi)(1 + eta_i eta) / 4.
N_i(xi, eta) = (1 + xi_i xi)(1 + eta_i eta) / 4