#!/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 Test_Quantity(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_integral(self): quantity = Quantity(self.mesh4) #Try constants first const = 5 quantity.set_values(const, location = 'vertices') #print 'Q', quantity.get_integral() assert allclose(quantity.get_integral(), self.mesh4.get_area() * const) #Try with a linear function def f(x, y): return x+y quantity.set_values(f, location = 'vertices') ref_integral = (4.0/3 + 8.0/3 + 10.0/3 + 10.0/3) * 2 assert allclose (quantity.get_integral(), ref_integral) 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_subset(self): quantity = Quantity(self.mesh4) quantity.set_vertex_values([0,1,2,3,4,5]) quantity.set_vertex_values([0,20,30,50], indexes = [0,2,3,5]) assert allclose(quantity.vertex_values, [[1,0,20], [1,20,4], [4,20,50], [30,1,4]]) 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]]) #Centroid assert allclose(quantity.centroid_values, [1., 7./3, 11./3, 8./3]) 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_limiter2(self): """Taken from test_shallow_water """ quantity = Conserved_quantity(self.mesh4) #Test centroids quantity.set_values([2.,4.,8.,2.], location = 'centroids') assert allclose(quantity.centroid_values, [2, 4, 8, 2]) #Centroid #Extrapolate quantity.extrapolate_second_order() assert allclose(quantity.vertex_values[1,:], [0.0, 6, 6]) #Limit quantity.limit() assert allclose(quantity.vertex_values[1,:], [2.2, 4.9, 4.9]) #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_distribute_second_order(self): quantity = Conserved_quantity(self.mesh4) #Test centroids quantity.set_values([2.,4.,8.,2.], location = 'centroids') assert allclose(quantity.centroid_values, [2, 4, 8, 2]) #Centroid #Extrapolate quantity.extrapolate_second_order() assert allclose(quantity.vertex_values[1,:], [0.0, 6, 6]) 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 stage = zeros(bed.shape, Float) h = 0.03 for i in range(stage.shape[0]): if i % 2 == 0: stage[i,:] = bed[i,:] + h else: stage[i,:] = bed[i,:] domain.set_quantity('stage', stage) stage = domain.quantities['stage'] #Get smoothed stage A, V = stage.get_vertex_values(xy=False, smooth=True) Q = stage.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 stage with XY X, Y, A1, V1 = stage.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 stage = zeros(bed.shape, Float) h = 0.03 for i in range(stage.shape[0]): if i % 2 == 0: stage[i,:] = bed[i,:] + h else: stage[i,:] = bed[i,:] domain.set_quantity('stage', stage) #Get stage stage = domain.quantities['stage'] A, V = stage.get_vertex_values(xy=False, smooth=False) Q = stage.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 = stage.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) #print "vertices",vertices #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]) #print "quantity",quantity.vertex_values 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]) quantity = Quantity(domain,[[1,1,1],[2,2,2],[3,3,3], [4,4,4],[5,5,5],[6,6,6]]) values = [10,100,50] #this will be per unique vertex, indexing the vertices #print "quantity.vertex_values",quantity.vertex_values quantity.set_values(values, indexes = [0,1,5]) #print "quantity.vertex_values",quantity.vertex_values assert allclose(quantity.vertex_values[0], [1,50,10]) assert allclose(quantity.vertex_values[5], [6,6,6]) assert allclose(quantity.vertex_values[1], [100,10,50]) quantity = Quantity(domain,[[1,1,1],[2,2,2],[3,3,3], [4,4,4],[5,5,5],[6,6,6]]) 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, [1,2,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_setting_unique_vertex_values(self): """ set values based on unique_vertex 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 "vertices",vertices #Create shallow water domain domain = Domain(points, vertices, boundary) #print "domain.number_of_elements ",domain.number_of_elements quantity = Quantity(domain,[[0,0,0],[1,1,1],[2,2,2],[3,3,3], [4,4,4],[5,5,5]]) value = 7 indexes = [1,5] quantity.set_values(value, location = 'unique vertices', indexes = indexes) #print "quantity.centroid_values",quantity.centroid_values assert allclose(quantity.vertex_values[0], [0,7,0]) assert allclose(quantity.vertex_values[1], [7,1,7]) assert allclose(quantity.vertex_values[2], [7,2,7]) def test_get_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,[[0,0,0],[1,1,1],[2,2,2],[3,3,3], [4,4,4],[5,5,5]]) #print "quantity.get_values(location = 'unique vertices')", \ # quantity.get_values(location = 'unique vertices') #print "quantity.get_values(location = 'unique vertices')", \ # quantity.get_values(indexes=[0,1,2,3,4,5,6,7], \ # location = 'unique vertices') answer = [0.5,2,4,5,0,1,3,4.5] assert allclose(answer, quantity.get_values(location = 'unique vertices')) indexes = [0,5,3] answer = [0.5,1,5] assert allclose(answer, quantity.get_values(indexes=indexes, \ location = 'unique vertices')) #print "quantity.centroid_values",quantity.centroid_values #print "quantity.get_values(location = 'centroids') ",\ # quantity.get_values(location = 'centroids') 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(Test_Quantity,'test') #print "restricted test" #suite = unittest.makeSuite(Test_Quantity,'test_set_vertex_values_subset') runner = unittest.TextTestRunner() runner.run(suite)