#!/usr/bin/env python import unittest from math import sqrt, pi from quantity import * #from config import epsilon from Numeric import allclose, array class TestCase(unittest.TestCase): def setUp(self): from domain import Domain a = 0.0 b = 1.0 c = 2.0 d = 2.5 e = 3.1 f = 4.0 self.points = [a, b, c, d, e, f] self.centroids = [(a+b)/2,(b+c)/2,(c+d)/2,(d+e)/2,(e+f)/2] self.vertex_values = [[1.0,2.0],[2.0,3.0],[3.0,4.0],[4.5,5],[5.5,5.6]] self.centroid_values = [[1.5, 2.5, 3.5, 4.75, 5.55]] self.domain1 = Domain(self.points[0:2]) self.domain5 = Domain(self.points) def tearDown(self): pass #print " Tearing down" def test_creation(self): quantity = Quantity(self.domain5, self.vertex_values) assert allclose(quantity.centroid_values, self.centroid_values) try: quantity = Quantity() except: pass else: raise 'Should have raised empty quantity exception' try: quantity = Quantity([1,2]) except AssertionError: pass except: raise 'Should have raised "missing domain object" error' def test_creation_zeros(self): quantity = Quantity(self.domain1) assert allclose(quantity.vertex_values, [[0.,0.]]) quantity = Quantity(self.domain5) assert allclose(quantity.vertex_values, [[0.,0.], [0.,0.], [0.,0.], [0.,0.], [0.,0.]]) def test_interpolation(self): quantity = Quantity(self.domain1, [[1,2]]) assert allclose(quantity.centroid_values, 1.5) #Centroid def test_interpolation2(self): quantity = Quantity(self.domain5, [[1,2], [5,5], [0,9], [-6, 3], [3,4]]) assert allclose(quantity.centroid_values, [1.5, 5., 4.5, -1.5, 3.5 ]) #Centroid ## def test_boundary_allocation(self): ## quantity = Conserved_quantity(self.mesh4, ## [[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]]) ## assert quantity.boundary_values.shape[0] == len(self.mesh4.boundary) def test_set_values(self): quantity = Quantity(self.domain5) quantity.set_values([[1,2], [5,5], [0,0], [-6, 3], [-2,4]], location = 'vertices') assert allclose(quantity.vertex_values, [[1,2], [5,5], [0,0], [-6, 3], [-2,4]]) assert allclose(quantity.centroid_values, [1.5, 5., 0., -1.5, 1.0]) #Centroid #Test default quantity.set_values([[1,2], [5,5], [0,0], [-6, 3], [-2,4]]) assert allclose(quantity.vertex_values, [[1,2], [5,5], [0,0], [-6, 3], [-2,4]]) assert allclose(quantity.centroid_values, [1.5, 5., 0., -1.5, 1.0]) #Centroid #Test centroids quantity.set_values([1,2,3,4,5], location = 'centroids') assert allclose(quantity.centroid_values, [1., 2., 3., 4., 5.]) #Centroid #Test exceptions try: quantity.set_values([[1,2], [5,5], [0,0], [-6, 3], [-2,4]], location = 'bas kamel tuba') except: pass try: quantity.set_values([[1,2], [0,0]]) except AssertionError: pass except: raise 'should have raised AssertionError' def test_set_values_const(self): quantity = Quantity(self.domain5) quantity.set_values(1.0, location = 'vertices') assert allclose(quantity.vertex_values, [[1,1], [1,1], [1,1], [1,1], [1,1]]) assert allclose(quantity.centroid_values, [1, 1, 1, 1, 1]) #Centroid quantity.set_values(2.0, location = 'centroids') assert allclose(quantity.centroid_values, [2, 2, 2, 2, 2]) def test_set_values_func(self): quantity = Quantity(self.domain5) def f(x): return x**2 quantity.set_values(f, location = 'vertices') assert allclose(quantity.vertex_values, [[0,1], [1,4], [4,6.25], [6.25,9.61], [9.61,16]]) assert allclose(quantity.centroid_values, [0.5, 2.5, 5.125, 7.93, 12.805 ]) quantity.set_values(f, location = 'centroids') assert allclose(quantity.centroid_values, [0.25, 1.5**2, 2.25**2, 2.8**2, 3.55**2]) ## def test_gradient(self): ## quantity = Conserved_quantity(self.mesh4) ## #Set up for a gradient of (3,0) at mid triangle ## quantity.set_values([2.0, 4.0, 8.0, 2.0], ## location = 'centroids') ## #Gradients ## a, b = quantity.compute_gradients() ## #gradient bewteen t0 and t1 is undefined as det==0 ## assert a[0] == 0.0 ## assert b[0] == 0.0 ## #The others are OK ## for i in range(1,4): ## assert a[i] == 3.0 ## assert b[i] == 0.0 ## quantity.extrapolate_second_order() ## assert allclose(quantity.vertex_values, [[2., 2., 2.], ## [0., 6., 6.], ## [6., 6., 12.], ## [0., 0., 6.]]) ## def test_second_order_extrapolation2(self): ## quantity = Conserved_quantity(self.mesh4) ## #Set up for a gradient of (3,1), f(x) = 3x+y ## quantity.set_values([2.0+2.0/3, 4.0+4.0/3, 8.0+2.0/3, 2.0+8.0/3], ## location = 'centroids') ## #Gradients ## a, b = quantity.compute_gradients() ## #gradient bewteen t0 and t1 is undefined as det==0 ## assert a[0] == 0.0 ## assert b[0] == 0.0 ## #The others are OK ## 163.968025 for i in range(1,4): ## assert allclose(a[i], 3.0) ## assert allclose(b[i], 1.0) ## quantity.extrapolate_second_order() ## assert allclose(quantity.vertex_values[1,0], 2.0) ## assert allclose(quantity.vertex_values[1,1], 6.0) ## assert allclose(quantity.vertex_values[1,2], 8.0) ## # def test_limiter(self): ## # initialise_consecutive_datastructure(points=6+4, elements=4) ## # a = Point (0.0, 0.0) ## # b = Point (0.0, 2.0) ## # c = Point (2.0, 0.0) ## # d = Point (0.0, 4.0) ## # e = Point (2.0, 2.0) ## # f = Point (4.0, 0.0) ## # #Set up for a gradient of (3,1), f(x) = 3x+y ## # v1 = Volume(b,a,c,array([0.0,0,0])) ## # v2 = Volume(b,c,e,array([1.0,0,0])) ## # v3 = Volume(e,c,f,array([10.0,0,0])) ## # v4 = Volume(d,b,e,array([0.0,0,0])) ## # #Setup neighbour structure ## # domain = Domain([v1,v2,v3,v4]) ## # domain.precompute() ## # #Lets's check first order first, hey ## # domain.order = 1 ## # domain.limiter = None ## # distribute_to_vertices_and_edges(domain) ## # assert allclose(v2.conserved_quantities_vertex0, ## # v2.conserved_quantities_centroid) ## # assert allclose(v2.conserved_quantities_vertex1, ## # v2.conserved_quantities_centroid) ## # assert allclose(v2.conserved_quantities_vertex2, ## # v2.conserved_quantities_centroid) ## # #Gradient of fitted pwl surface ## # a, b = compute_gradient(v2.id) ## # assert abs(a[0] - 5.0) < epsilon ## # assert abs(b[0]) < epsilon ## # #assert qminr[0] == 0.0 ## # #assert qmaxr[0] == 10.0 ## # #And now for the second order stuff ## # # - the full second order extrapolation ## # domain.order = 2 ## # distribute_to_vertices_and_edges(domain) ## # qmin = qmax = v2.conserved_quantities_centroid ## # qmin = minimum(qmin, v1.conserved_quantities_centroid) ## # qmax = maximum(qmax, v1.conserved_quantities_centroid) ## # qmin = minimum(qmin, v3.conserved_quantities_centroid) ## # qmax = maximum(qmax, v3.conserved_quantities_centroid) ## # qmin = minimum(qmin, v4.conserved_quantities_centroid) ## # qmax = maximum(qmax, v4.conserved_quantities_centroid) ## # #assert qminr == qmin ## # #assert qmaxr == qmax ## # assert v2.conserved_quantities_vertex0 <= qmax ## # assert v2.conserved_quantities_vertex0 >= qmin ## # assert v2.conserved_quantities_vertex1 <= qmax ## # assert v2.conserved_quantities_vertex1 >= qmin ## # assert v2.conserved_quantities_vertex2 <= qmax ## # assert v2.conserved_quantities_vertex2 >= qmin 163.968025 ## # #Check that volume has been preserved ## # q = v2.conserved_quantities_centroid[0] ## # w = (v2.conserved_quantities_vertex0[0] + ## # v2.conserved_quantities_vertex1[0] + ## # v2.conserved_quantities_vertex2[0])/3 ## # assert allclose(q, w) ## def test_first_order_extrapolator(self): ## quantity = Conserved_quantity(self.mesh4) ## #Test centroids ## quantity.set_values([1.,2.,3.,4.], location = 'centroids') ## assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid ## #Extrapolate ## quantity.extrapolate_first_order() ## #Check vertices but not edge values ## assert allclose(quantity.vertex_values, ## [[1,1,1], [2,2,2], [3,3,3], [4, 4, 4]]) ## def test_second_order_extrapolator(self): ## quantity = Conserved_quantity(self.mesh4) ## #Set up for a gradient of (3,0) at mid triangle ## quantity.set_values([2.0, 4.0, 8.0, 2.0], ## location = 'centroids') ## quantity.extrapolate_second_order() ## quantity.limit() ## #Assert that central triangle is limited by neighbours ## assert quantity.vertex_values[1,0] >= quantity.vertex_values[0,0] ## assert quantity.vertex_values[1,0] >= quantity.vertex_values[3,1] ## assert quantity.vertex_values[1,1] <= quantity.vertex_values[2,1] ## assert quantity.vertex_values[1,1] >= quantity.vertex_values[0,2] ## assert quantity.vertex_values[1,2] <= quantity.vertex_values[2,0] ## assert quantity.vertex_values[1,2] >= quantity.vertex_values[3,1] ## #Assert that quantities are conserved ## from Numeric import sum ## for k in range(quantity.centroid_values.shape[0]): ## assert allclose (quantity.centroid_values[k], ## sum(quantity.vertex_values[k,:])/3) ## def test_limiter(self): ## quantity = Conserved_quantity(self.mesh4) ## #Create a deliberate overshoot (e.g. from gradient computation) ## quantity.set_values([[3,0,3], [2,2,6], [5,3,8], [8,3,5]]) ## #Limit ## quantity.limit() ## #Assert that central triangle is limited by neighbours ## assert quantity.vertex_values[1,0] >= quantity.vertex_values[0,0] ## assert quantity.vertex_values[1,0] <= quantity.vertex_values[3,1] ## assert quantity.vertex_values[1,1] <= quantity.vertex_values[2,1] ## assert quantity.vertex_values[1,1] >= quantity.vertex_values[0,2] ## assert quantity.vertex_values[1,2] <= quantity.vertex_values[2,0] ## assert quantity.vertex_values[1,2] <= quantity.vertex_values[3,1] ## #Assert that quantities are conserved ## from Numeric import sum ## for k in range(quantity.centroid_values.shape[0]): ## assert allclose (quantity.centroid_values[k], ## sum(quantity.vertex_values[k,:])/3) ## def test_distribute_first_order(self): ## quantity = Conserved_quantity(self.mesh4) ## #Test centroids ## quantity.set_values([1.,2.,3.,4.], location = 'centroids') ## assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid ## #Extrapolate ## quantity.extrapolate_first_order() ## #Interpolate ## quantity.interpolate_from_vertices_to_edges() ## assert allclose(quantity.vertex_values, ## [[1,1,1], [2,2,2], [3,3,3], [4, 4, 4]]) ## assert allclose(quantity.edge_values, [[1,1,1], [2,2,2], ## [3,3,3], [4, 4, 4]]) ## def test_update_explicit(self): ## quantity = Conserved_quantity(self.mesh4) ## #Test centroids ## quantity.set_values([1.,2.,3.,4.], location = 'centroids') ## assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid ## #Set explicit_update ## quantity.explicit_update = array( [1.,1.,1.,1.] ) ## #Update with given timestep ## quantity.update(0.1) ## x = array([1, 2, 3, 4]) + array( [.1,.1,.1,.1] ) ## assert allclose( quantity.centroid_values, x) ## def test_update_semi_implicit(self): ## quantity = Conserved_quantity(self.mesh4) ## #Test centroids ## quantity.set_values([1.,2.,3.,4.], location = 'centroids') ## assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid ## #Set semi implicit update ## quantity.semi_implicit_update = array( [1.,1.,1.,1.] ) ## #Update with given timestep ## quantity.update(0.1) ## x = array([1, 2, 3, 4])/array( [.9,.9,.9,.9] ) ## assert allclose( quantity.centroid_values, x) ## def test_both_updates(self): ## quantity = Conserved_quantity(self.mesh4) ## #Test centroids ## quantity.set_values([1.,2.,3.,4.], location = 'centroids') ## assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid ## #Set explicit_update ## quantity.explicit_update = array( [4.,3.,2.,1.] ) ## #Set semi implicit update ## quantity.semi_implicit_update = array( [1.,1.,1.,1.] ) ## #Update with given timestep ## quantity.update(0.1) ## x = array([1, 2, 3, 4]) + array( [.4,.3,.2,.1] ) ## x /= array( [.9,.9,.9,.9] ) ## assert allclose( quantity.centroid_values, x) #------------------------------------------------------------- if __name__ == "__main__": suite = unittest.makeSuite(TestCase,'test') runner = unittest.TextTestRunner() runner.run(suite)