smbinterp/test/baker2dorder.py

#!/usr/bin/env python

import unittest

from interp import baker
from interp.grid import grid
import numpy as np

from interp.grid import contains


def exact_func(point):
    x = point[0]
    y = point[1]
    return 0.5 + x * x + y


def calculate_error_term(self, a, b, c, d, e, f):
    B = np.array([
     self.p1[a] * self.p1[b], self.p1[c] * self.p1[d], self.p1[e] * self.p1[f],
     self.p2[a] * self.p2[b], self.p2[c] * self.p2[d], self.p2[e] * self.p2[f],
     self.p3[a] * self.p3[b], self.p3[c] * self.p3[d], self.p3[e] * self.p3[f],
     self.p4[a] * self.p4[b], self.p4[c] * self.p4[d], self.p4[e] * self.p4[f],
    ])

    B.shape = (4, 3)

    A = np.dot(B.T, B)
    rhs = np.dot(B.T, self.w)
    abc = np.linalg.solve(A, rhs)

    err = \
        abc[0] * self.phis[a] * self.phis[b] + \
        abc[1] * self.phis[c] * self.phis[d] + \
        abc[2] * self.phis[e] * self.phis[f]
    return err


class Test(unittest.TestCase):
    def setUp(self):
        self.verts = [
                         [2, 3],   # 0
                         [7, 4],   # 1
                         [4, 8],   # 2
                         [0, 7],   # 3, 1
                         [5, 0],   # 4, 2
                         [0, 5],   # 5, 3
                         [8, 9],   # 6, 4
                     ]

        self.q = [exact_func(v) for v in self.verts]

        self.g = grid()
        self.g.verts = self.verts
        self.g.q = self.q

        self.R, self.R_q = (self.verts[:3], self.q[:3])
        self.S, self.S_q = (self.verts[3:], self.q[3:])

        self.p1, self.ql1 = baker.qlinear(self.verts[3], self.R, self.q)
        self.p2, self.ql2 = baker.qlinear(self.verts[4], self.R, self.q)
        self.p3, self.ql3 = baker.qlinear(self.verts[5], self.R, self.q)
        self.p4, self.ql4 = baker.qlinear(self.verts[6], self.R, self.q)

        self.q1 = exact_func(self.verts[3])
        self.q2 = exact_func(self.verts[4])
        self.q3 = exact_func(self.verts[5])
        self.q4 = exact_func(self.verts[6])

        self.w = np.array([
            self.q1 - self.ql1,
            self.q2 - self.ql2,
            self.q3 - self.ql3,
            self.q4 - self.ql4,
        ])

        self.X = [4, 5]

        self.phis, self.qlin = baker.qlinear(self.X, self.R, self.q)
        self.exact = exact_func(self.X)
        self.answer = baker.interpolate(self.X, self.R,
                                       self.R_q, self.S, self.S_q)

    def test_R_contains_X(self):
        self.assertTrue(contains(self.X, self.R))

    def test_1(self):
        a, b, c, d, e, f = (0, 1, 1, 2, 2, 0)
        err = calculate_error_term(self, a, b, c, d, e, f)
        self.assertAlmostEqual(err, self.answer.error)

    def test_swap_first_elements(self):
        a, b, c, d, e, f = (1, 0, 1, 2, 2, 0)
        err = calculate_error_term(self, a, b, c, d, e, f)
        self.assertAlmostEqual(err, self.answer.error)

    def test_swap_two_pairs(self):
        a, b, c, d, e, f = (1, 2, 0, 1, 2, 0)
        err = calculate_error_term(self, a, b, c, d, e, f)
        self.assertAlmostEqual(err, self.answer.error)

    def test_swap_all_pairs(self):
        a, b, c, d, e, f = (0, 2, 0, 1, 2, 1)
        err = calculate_error_term(self, a, b, c, d, e, f)
        self.assertAlmostEqual(err, self.answer.error)


if __name__ == '__main__':
    suite = unittest.TestLoader().loadTestsFromTestCase(Test)
    unittest.TextTestRunner(verbosity=3).run(suite)