#!/usr/bin/env python import unittest from math import sqrt, pi from quantity import * from config import epsilon from Numeric import allclose, array, ones, Float class TestCase(unittest.TestCase): def setUp(self): from domain import Domain a = [0.0, 0.0] b = [0.0, 2.0] c = [2.0,0.0] d = [0.0, 4.0] e = [2.0, 2.0] f = [4.0,0.0] points = [a, b, c, d, e, f] #bac, bce, ecf, dbe elements = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ] self.mesh1 = Domain(points[:3], [elements[0]]) self.mesh1.check_integrity() self.mesh4 = Domain(points, elements) self.mesh4.check_integrity() def tearDown(self): pass #print " Tearing down" def test_creation(self): quantity = Quantity(self.mesh1, [[1,2,3]]) assert allclose(quantity.vertex_values, [[1.,2.,3.]]) try: quantity = Quantity() except: pass else: raise 'Should have raised empty quantity exception' try: quantity = Quantity([1,2,3]) except AssertionError: pass except: raise 'Should have raised "mising mesh object" error' def test_creation_zeros(self): quantity = Quantity(self.mesh1) assert allclose(quantity.vertex_values, [[0.,0.,0.]]) quantity = Quantity(self.mesh4) assert allclose(quantity.vertex_values, [[0.,0.,0.], [0.,0.,0.], [0.,0.,0.], [0.,0.,0.]]) def test_interpolation(self): quantity = Quantity(self.mesh1, [[1,2,3]]) assert allclose(quantity.centroid_values, [2.0]) #Centroid assert allclose(quantity.edge_values, [[2.5, 2.0, 1.5]]) def test_interpolation2(self): quantity = Conserved_quantity(self.mesh4, [[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]]) assert allclose(quantity.centroid_values, [2., 5., 3., 0.]) #Centroid quantity.extrapolate_second_order() assert allclose(quantity.vertex_values, [[2., 2., 2.], [3.+2./3, 6.+2./3, 4.+2./3], [7.5, 0.5, 1.], [-5, -2.5, 7.5]]) #print quantity.edge_values assert allclose(quantity.edge_values, [[2.5, 2.0, 1.5], [5., 5., 5.], [4.5, 4.5, 0.], [3.0, -1.5, -1.5]]) 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.mesh4) quantity.set_values([[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]], location = 'vertices') assert allclose(quantity.vertex_values, [[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]]) assert allclose(quantity.centroid_values, [2., 5., 3., 0.]) #Centroid assert allclose(quantity.edge_values, [[2.5, 2.0, 1.5], [5., 5., 5.], [4.5, 4.5, 0.], [3.0, -1.5, -1.5]]) #Test default quantity.set_values([[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]]) assert allclose(quantity.vertex_values, [[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]]) assert allclose(quantity.centroid_values, [2., 5., 3., 0.]) #Centroid assert allclose(quantity.edge_values, [[2.5, 2.0, 1.5], [5., 5., 5.], [4.5, 4.5, 0.], [3.0, -1.5, -1.5]]) #Test centroids quantity.set_values([1,2,3,4], location = 'centroids') assert allclose(quantity.centroid_values, [1., 2., 3., 4.]) #Centroid #Test edges quantity.set_values([[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]], location = 'edges') assert allclose(quantity.edge_values, [[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]]) #Test exceptions try: quantity.set_values([[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]], location = 'bas kamel tuba') except: pass try: quantity.set_values([[1,2,3], [0,0,9]]) except AssertionError: pass except: raise 'should have raised Assertionerror' def test_set_values_const(self): quantity = Quantity(self.mesh4) quantity.set_values(1.0, location = 'vertices') assert allclose(quantity.vertex_values, [[1,1,1], [1,1,1], [1,1,1], [1, 1, 1]]) assert allclose(quantity.centroid_values, [1, 1, 1, 1]) #Centroid assert allclose(quantity.edge_values, [[1, 1, 1], [1, 1, 1], [1, 1, 1], [1, 1, 1]]) quantity.set_values(2.0, location = 'centroids') assert allclose(quantity.centroid_values, [2, 2, 2, 2]) quantity.set_values(3.0, location = 'edges') assert allclose(quantity.edge_values, [[3, 3, 3], [3, 3, 3], [3, 3, 3], [3, 3, 3]]) def test_set_values_func(self): quantity = Quantity(self.mesh4) def f(x, y): return x+y quantity.set_values(f, location = 'vertices') #print "quantity.vertex_values",quantity.vertex_values assert allclose(quantity.vertex_values, [[2,0,2], [2,2,4], [4,2,4], [4,2,4]]) assert allclose(quantity.centroid_values, [4.0/3, 8.0/3, 10.0/3, 10.0/3]) assert allclose(quantity.edge_values, [[1,2,1], [3,3,2], [3,4,3], [3,4,3]]) quantity.set_values(f, location = 'centroids') assert allclose(quantity.centroid_values, [4.0/3, 8.0/3, 10.0/3, 10.0/3]) def test_set_vertex_values(self): quantity = Quantity(self.mesh4) quantity.set_vertex_values([0,1,2,3,4,5]) assert allclose(quantity.vertex_values, [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]) assert allclose(quantity.centroid_values, [1., 7./3, 11./3, 8./3]) #Centroid assert allclose(quantity.edge_values, [[1., 1.5, 0.5], [3., 2.5, 1.5], [3.5, 4.5, 3.], [2.5, 3.5, 2]]) def test_set_vertex_values_using_general_interface(self): quantity = Quantity(self.mesh4) quantity.set_values([0,1,2,3,4,5]) assert allclose(quantity.vertex_values, [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]) assert allclose(quantity.centroid_values, [1., 7./3, 11./3, 8./3]) #Centroid assert allclose(quantity.edge_values, [[1., 1.5, 0.5], [3., 2.5, 1.5], [3.5, 4.5, 3.], [2.5, 3.5, 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() #print quantity.vertex_values 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 for i in range(1,4): assert allclose(a[i], 3.0) assert allclose(b[i], 1.0) quantity.extrapolate_second_order() #print quantity.vertex_values 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_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 timestep = 0.1 quantity.update(timestep) sem = array([1.,1.,1.,1.])/array([1, 2, 3, 4]) denom = ones(4, Float)-timestep*sem x = array([1, 2, 3, 4])/denom 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 timestep = 0.1 quantity.update(0.1) sem = array([1.,1.,1.,1.])/array([1, 2, 3, 4]) denom = ones(4, Float)-timestep*sem x = array([1, 2, 3, 4]) + array( [.4,.3,.2,.1] ) x /= denom assert allclose( quantity.centroid_values, x) #Test smoothing def test_smoothing(self): from mesh_factory import rectangular from shallow_water import Domain, Transmissive_boundary from Numeric import zeros, Float from util import mean #Create basic mesh points, vertices, boundary = rectangular(2, 2) #Create shallow water domain domain = Domain(points, vertices, boundary) domain.default_order=2 domain.reduction = mean #Set some field values domain.set_quantity('elevation', lambda x,y: x) domain.set_quantity('friction', 0.03) ###################### # Boundary conditions B = Transmissive_boundary(domain) domain.set_boundary( {'left': B, 'right': B, 'top': B, 'bottom': B}) ###################### #Initial condition - with jumps bed = domain.quantities['elevation'].vertex_values level = zeros(bed.shape, Float) h = 0.03 for i in range(level.shape[0]): if i % 2 == 0: level[i,:] = bed[i,:] + h else: level[i,:] = bed[i,:] domain.set_quantity('level', level) level = domain.quantities['level'] #Get smoothed level A, V = level.get_vertex_values(xy=False, smooth=True) Q = level.vertex_values assert A.shape[0] == 9 assert V.shape[0] == 8 assert V.shape[1] == 3 #First four points assert allclose(A[0], (Q[0,2] + Q[1,1])/2) assert allclose(A[1], (Q[1,0] + Q[3,1] + Q[2,2])/3) assert allclose(A[2], Q[3,0]) assert allclose(A[3], (Q[0,0] + Q[5,1] + Q[4,2])/3) #Center point assert allclose(A[4], (Q[0,1] + Q[1,2] + Q[2,0] +\ Q[5,0] + Q[6,2] + Q[7,1])/6) #Check V assert allclose(V[0,:], [3,4,0]) assert allclose(V[1,:], [1,0,4]) assert allclose(V[2,:], [4,5,1]) assert allclose(V[3,:], [2,1,5]) assert allclose(V[4,:], [6,7,3]) assert allclose(V[5,:], [4,3,7]) assert allclose(V[6,:], [7,8,4]) assert allclose(V[7,:], [5,4,8]) #Get smoothed level with XY X, Y, A1, V1 = level.get_vertex_values(xy=True, smooth=True) assert allclose(A, A1) assert allclose(V, V1) #Check XY assert allclose(X[4], 0.5) assert allclose(Y[4], 0.5) assert allclose(X[7], 1.0) assert allclose(Y[7], 0.5) def test_vertex_values_no_smoothing(self): from mesh_factory import rectangular from shallow_water import Domain, Transmissive_boundary from Numeric import zeros, Float from util import mean #Create basic mesh points, vertices, boundary = rectangular(2, 2) #Create shallow water domain domain = Domain(points, vertices, boundary) domain.default_order=2 domain.reduction = mean #Set some field values domain.set_quantity('elevation', lambda x,y: x) domain.set_quantity('friction', 0.03) ###################### #Initial condition - with jumps bed = domain.quantities['elevation'].vertex_values level = zeros(bed.shape, Float) h = 0.03 for i in range(level.shape[0]): if i % 2 == 0: level[i,:] = bed[i,:] + h else: level[i,:] = bed[i,:] domain.set_quantity('level', level) #Get level level = domain.quantities['level'] A, V = level.get_vertex_values(xy=False, smooth=False) Q = level.vertex_values.flat for k in range(8): assert allclose(A[k], Q[k]) for k in range(8): assert V[k, 0] == 3*k assert V[k, 1] == 3*k+1 assert V[k, 2] == 3*k+2 X, Y, A1, V1 = level.get_vertex_values(xy=True, smooth=False) assert allclose(A, A1) assert allclose(V, V1) #Check XY assert allclose(X[1], 0.5) assert allclose(Y[1], 0.5) assert allclose(X[4], 0.0) assert allclose(Y[4], 0.0) assert allclose(X[12], 1.0) assert allclose(Y[12], 0.0) def set_array_values_by_index(self): from mesh_factory import rectangular from shallow_water import Domain from Numeric import zeros, Float #Create basic mesh points, vertices, boundary = rectangular(1, 1) #Create shallow water domain domain = Domain(points, vertices, boundary) #print "domain.number_of_elements ",domain.number_of_elements quantity = Quantity(domain,[[1,1,1],[2,2,2]]) value = [7] indexes = [1] quantity.set_array_values_by_index(value, location = 'centroids', indexes = indexes) #print "quantity.centroid_values",quantity.centroid_values assert allclose(quantity.centroid_values, [1,7]) quantity.set_array_values([15,20,25], indexes = indexes) assert allclose(quantity.centroid_values, [1,20]) quantity.set_array_values([15,20,25], indexes = indexes) assert allclose(quantity.centroid_values, [1,20]) def test_setting_some_vertex_values(self): """ set values based on triangle lists. """ from mesh_factory import rectangular from shallow_water import Domain from Numeric import zeros, Float #Create basic mesh points, vertices, boundary = rectangular(1, 3) #Create shallow water domain domain = Domain(points, vertices, boundary) #print "domain.number_of_elements ",domain.number_of_elements quantity = Quantity(domain,[[1,1,1],[2,2,2],[3,3,3], [4,4,4],[5,5,5],[6,6,6]]) value = [7] indexes = [1] quantity.set_values(value, location = 'centroids', indexes = indexes) #print "quantity.centroid_values",quantity.centroid_values assert allclose(quantity.centroid_values, [1,7,3,4,5,6]) value = [[15,20,25]] quantity.set_values(value, indexes = indexes) #print "1 quantity.vertex_values",quantity.vertex_values assert allclose(quantity.vertex_values[1], value[0]) # FIXMEDSG this seems wrong -DSG # It's right. The values is an array, indexed by the element id's # given in indexes. It sets values per triangle. values = [10,100,50] quantity.set_values(values, indexes = [0,1,5]) # , location = 'centroids' #print "2 quantity.vertex_values",quantity.vertex_values assert allclose(quantity.vertex_values[0], [10,10,10]) assert allclose(quantity.vertex_values[5], [50,50,50]) #quantity.interpolate() #print "quantity.centroid_values",quantity.centroid_values assert allclose(quantity.centroid_values, [10,100,3,4,5,50]) values = [[31,30,29],[400,400,400],[1000,999,998]] quantity.set_values(values, indexes = [3,3,5]) quantity.interpolate() assert allclose(quantity.centroid_values, [10,100,3,400,5,999]) values = [[1,1,1],[2,2,2],[3,3,3], [4,4,4],[5,5,5],[6,6,6]] quantity.set_values(values) # testing the standard set values by vertex # indexed by vertex_id in general_mesh.coordinates values = [0,1,2,3,4,5,6,7] quantity.set_values(values) #print "1 quantity.vertex_values",quantity.vertex_values assert allclose(quantity.vertex_values,[[ 4., 5., 0.], [ 1., 0., 5.], [ 5., 6., 1.], [ 2., 1., 6.], [ 6., 7., 2.], [ 3., 2., 7.]]) def test_getting_some_vertex_values(self): """ get values based on triangle lists. """ from mesh_factory import rectangular from shallow_water import Domain from Numeric import zeros, Float #Create basic mesh points, vertices, boundary = rectangular(1, 3) #print "points",points #print "vertices",vertices #print "boundary",boundary #Create shallow water domain domain = Domain(points, vertices, boundary) #print "domain.number_of_elements ",domain.number_of_elements quantity = Quantity(domain,[[1,1,1],[2,2,2],[3,3,3], [4,4,4],[5,5,5],[6,6,6]]) value = [7] indexes = [1] quantity.set_values(value, location = 'centroids', indexes = indexes) #print "quantity.centroid_values",quantity.centroid_values #print "quantity.get_values(location = 'centroids') ",\ # quantity.get_values(location = 'centroids') assert allclose(quantity.centroid_values, quantity.get_values(location = 'centroids')) value = [[15,20,25]] quantity.set_values(value, indexes = indexes) #print "1 quantity.vertex_values",quantity.vertex_values assert allclose(quantity.vertex_values, quantity.get_values()) assert allclose(quantity.edge_values, quantity.get_values(location = 'edges')) # get a subset of elements subset = quantity.get_values(location='centroids', indexes=[0,5]) answer = [quantity.centroid_values[0],quantity.centroid_values[5]] assert allclose(subset, answer) subset = quantity.get_values(location='edges', indexes=[0,5]) answer = [quantity.edge_values[0],quantity.edge_values[5]] #print "subset",subset #print "answer",answer assert allclose(subset, answer) subset = quantity.get_values( indexes=[1,5]) answer = [quantity.vertex_values[1],quantity.vertex_values[5]] #print "subset",subset #print "answer",answer assert allclose(subset, answer) #------------------------------------------------------------- if __name__ == "__main__": suite = unittest.makeSuite(TestCase,'test') runner = unittest.TextTestRunner() runner.run(suite)