Changeset 7793 for anuga_work/development
- Timestamp:
- Jun 6, 2010, 10:48:41 PM (15 years ago)
- Location:
- anuga_work/development/anuga_1d
- Files:
-
- 3 edited
Legend:
- Unmodified
- Added
- Removed
-
anuga_work/development/anuga_1d/domain.py
r6694 r7793 112 112 self.build_tagged_elements_dictionary(tagged_elements) 113 113 114 from quantity import Quantity , Conserved_quantity115 #from quantity_domain import Quantity, Conserved_quantity 114 from quantity import Quantity 115 116 116 117 117 #List of quantity names entering … … 592 592 return self.vertices[elem_id,vertex] 593 593 594 def get_area(self, elem_id): 595 """Return area of element id 596 """ 597 598 return self.areas[elem_id] 594 def get_area(self, elem_id=None): 595 """Return total domain area or area of element id 596 """ 597 598 if elem_id is None: 599 return sum(self.areas) 600 else: 601 return self.areas[elem_id] 599 602 600 603 def get_quantity(self, name, location='vertices', indices = None): -
anuga_work/development/anuga_1d/quantity.py
r6453 r7793 1 """Class Quantity - Implements values at each 1d element 1 """ 2 Class Quantity - Implements values at each 1d element 2 3 3 4 To create: … … 85 86 """ 86 87 return self.centroid_values.shape[0] 87 88 89 def __neg__(self): 90 """Negate all values in this quantity giving meaning to the 91 expression -Q where Q is an instance of class Quantity 92 """ 93 94 Q = Quantity(self.domain) 95 Q.set_values_from_numeric(-self.vertex_values) 96 return Q 97 98 99 100 def __add__(self, other): 101 """Add to self anything that could populate a quantity 102 103 E.g other can be a constant, an array, a function, another quantity 104 (except for a filename or points, attributes (for now)) 105 - see set_values_from_numeric for details 106 """ 107 108 Q = Quantity(self.domain) 109 Q.set_values_from_numeric(other) 110 111 result = Quantity(self.domain) 112 result.set_values_from_numeric(self.vertex_values + Q.vertex_values) 113 return result 114 115 def __radd__(self, other): 116 """Handle cases like 7+Q, where Q is an instance of class Quantity 117 """ 118 119 return self + other 120 121 def __sub__(self, other): 122 return self + -other # Invoke self.__neg__() 123 124 def __mul__(self, other): 125 """Multiply self with anything that could populate a quantity 126 127 E.g other can be a constant, an array, a function, another quantity 128 (except for a filename or points, attributes (for now)) 129 - see set_values_from_numeric for details 130 """ 131 132 if isinstance(other, Quantity): 133 Q = other 134 else: 135 Q = Quantity(self.domain) 136 Q.set_values_from_numeric(other) 137 138 result = Quantity(self.domain) 139 140 # The product of vertex_values, edge_values and centroid_values 141 # are calculated and assigned directly without using 142 # set_values_from_numeric (which calls interpolate). Otherwise 143 # centroid values wouldn't be products from q1 and q2 144 result.vertex_values = self.vertex_values * Q.vertex_values 145 result.centroid_values = self.centroid_values * Q.centroid_values 146 147 return result 148 149 def __rmul__(self, other): 150 """Handle cases like 3*Q, where Q is an instance of class Quantity 151 """ 152 153 return self * other 154 155 def __div__(self, other): 156 """Divide self with anything that could populate a quantity 157 158 E.g other can be a constant, an array, a function, another quantity 159 (except for a filename or points, attributes (for now)) 160 - see set_values_from_numeric for details 161 162 Zero division is dealt with by adding an epsilon to the divisore 163 FIXME (Ole): Replace this with native INF once we migrate to NumPy 164 """ 165 166 if isinstance(other, Quantity): 167 Q = other 168 else: 169 Q = Quantity(self.domain) 170 Q.set_values_from_numeric(other) 171 172 result = Quantity(self.domain) 173 174 # The quotient of vertex_values, edge_values and centroid_values 175 # are calculated and assigned directly without using 176 # set_values_from_numeric (which calls interpolate). Otherwise 177 # centroid values wouldn't be quotient of q1 and q2 178 result.vertex_values = self.vertex_values/(Q.vertex_values + epsilon) 179 result.centroid_values = self.centroid_values/(Q.centroid_values + epsilon) 180 181 return result 182 183 def __rdiv__(self, other): 184 """Handle cases like 3/Q, where Q is an instance of class Quantity 185 """ 186 187 return self / other 188 189 def __pow__(self, other): 190 """Raise quantity to (numerical) power 191 192 As with __mul__ vertex values are processed entry by entry 193 while centroid and edge values are re-interpolated. 194 195 Example using __pow__: 196 Q = (Q1**2 + Q2**2)**0.5 197 """ 198 199 if isinstance(other, Quantity): 200 Q = other 201 else: 202 Q = Quantity(self.domain) 203 Q.set_values_from_numeric(other) 204 205 result = Quantity(self.domain) 206 207 # The power of vertex_values, edge_values and centroid_values 208 # are calculated and assigned directly without using 209 # set_values_from_numeric (which calls interpolate). Otherwise 210 # centroid values wouldn't be correct 211 result.vertex_values = self.vertex_values ** other 212 result.centroid_values = self.centroid_values ** other 213 214 return result 215 216 def set_values_from_numeric(self, numeric): 217 import Numeric 218 x = Numeric.array([1.0,2.0]) 219 y = [1.0,2.0] 220 221 if type(numeric) == type(y): 222 self.set_values_from_array(numeric) 223 elif type(numeric) == type(x): 224 self.set_values_from_array(numeric) 225 elif callable(numeric): 226 self.set_values_from_function(numeric) 227 elif isinstance(numeric, Quantity): 228 self.set_values_from_quantity(numeric) 229 else: # see if it's coercible to a float (float, int or long, etc) 230 try: 231 numeric = float(numeric) 232 except ValueError: 233 msg = ("Illegal type for variable 'numeric': %s" % type(numeric)) 234 raise Exception, msg 235 self.set_values_from_constant(numeric) 236 237 def set_values_from_constant(self,numeric): 238 239 self.vertex_values[:] = numeric 240 self.centroid_values[:] = numeric 241 242 243 def set_values_from_array(self,numeric): 244 245 self.vertex_values[:] = numeric 246 self.interpolate() 247 248 249 def set_values_from_quantity(self,quantity): 250 251 self.vertex_values[:] = quantity.vertex_values 252 self.centroid_values[:] = quantity.centroid_values 253 254 def set_values_from_function(self,function): 255 256 self.vertex_values[:] = map(function, self.domain.vertices) 257 self.centroid_values[:] = map(function, self.domain.centroids) 258 259 88 260 def interpolate(self): 89 261 """ … … 218 390 assert values.shape[1] == 2, msg 219 391 220 self.vertex_values = values392 self.vertex_values[:] = values 221 393 222 394 … … 428 600 """ 429 601 430 #print 'compute_gradient'431 602 432 603 from Numeric import array, zeros, Float … … 908 1079 self.centroid_values[:] = (a*self.centroid_values + b*self.centroid_backup_values).astype('f') 909 1080 910 class Conserved_quantity(Quantity): 911 """Class conserved quantity adds to Quantity: 912 913 storage and method for updating, and 914 methods for extrapolation from centropid to vertices inluding 915 gradients and limiters 916 """ 917 918 def __init__(self, domain, vertex_values=None): 919 Quantity.__init__(self, domain, vertex_values) 920 921 print "Use Quantity instead of Conserved_quantity" 922 923 924 ## 925 ##def newLinePlot(title='Simple Plot'): 926 ## import Gnuplot 927 ## g = Gnuplot.Gnuplot() 928 ## g.title(title) 929 ## g('set data style linespoints') 930 ## g.xlabel('x') 931 ## g.ylabel('y') 932 ## return g 933 ## 934 ##def linePlot(g,x,y): 935 ## import Gnuplot 936 ## g.plot(Gnuplot.PlotItems.Data(x.flat,y.flat)) 1081 937 1082 938 1083 def newLinePlot(title='Simple Plot'): -
anuga_work/development/anuga_1d/test_quantity.py
r5858 r7793 1 import quantity 1 2 #!/usr/bin/env python 2 3 3 4 import unittest 4 5 from math import sqrt, pi 5 import tempfile 6 6 7 7 8 from quantity import * 8 from config import epsilon 9 from Numeric import allclose, array, ones, Float 10 11 from domain import Domain 12 13 14 #Aux for fit_interpolate.fit example 15 def linear_function(point): 16 point = array(point) 17 return point[:,0]+point[:,1] 9 from domain import * 10 11 12 from Numeric import allclose, array, ones, Float, maximum, zeros 18 13 19 14 20 15 class Test_Quantity(unittest.TestCase): 21 16 def setUp(self): 22 from domain import Domain 23 24 a = [0.0, 0.0] 25 b = [0.0, 2.0] 26 c = [2.0, 0.0] 27 d = [0.0, 4.0] 28 e = [2.0, 2.0] 29 f = [4.0, 0.0] 30 31 points = [a, b, c, d, e, f] 32 33 #bac, bce, ecf, dbe 34 elements = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ] 35 36 self.mesh1 = Domain(points[:3], [elements[0]]) 37 self.mesh1.check_integrity() 38 39 self.mesh4 = Domain(points, elements) 40 self.mesh4.check_integrity() 41 42 # UTM round Onslow 43 a = [240000, 7620000] 44 b = [240000, 7680000] 45 c = [300000, 7620000] 46 47 points = [a, b, c] 48 elements = [[0,2,1]] 49 50 self.mesh_onslow = Domain(points, elements) 51 self.mesh_onslow.check_integrity() 52 17 self.points3 = [0.0, 1.0, 2.0, 3.0] 18 self.vertex_values3 = [[1.0,2.0],[4.0,5.0],[-1.0,2.0]] 19 self.domain3 = Domain(self.points3) 20 21 22 23 self.points4 = [0.0, 1.0, 2.5, 3.0, 5.0] 24 self.vertex_values4 = [[1.0,2.0],[4.0,5.0],[-1.0,2.0],[3.0,6.0]] 25 self.centroid_values4 = [1.5, 4.5, 0.5, 4.5] 26 self.boundary4 = {(0, 0): 'left', (3, 1): 'right'} 27 self.domain4 = Domain(self.points4,self.boundary4) 28 29 self.points10 = [0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0] 30 self.domain10 = Domain(self.points10) 31 32 self.points6 = [0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0] 33 self.domain6 = Domain(self.points6) 34 35 53 36 def tearDown(self): 54 37 pass … … 56 39 57 40 41 def test_creat_with_boundary(self): 42 43 assert self.domain4.boundary == {(0, 0): 'left', (3, 1): 'right'} 44 58 45 def test_creation(self): 59 46 60 quantity = Quantity(self.mesh1, [[1,2,3]]) 61 assert allclose(quantity.vertex_values, [[1.,2.,3.]]) 62 47 quantity = Quantity(self.domain3) 48 assert allclose(quantity.vertex_values, [[0.0,0.0],[0.0,0.0],[0.0,0.0]]) 49 50 63 51 try: 64 52 quantity = Quantity() … … 74 62 pass 75 63 except: 76 raise 'Should have raised "mising meshobject" error'64 raise 'Should have raised "mising domain object" error' 77 65 78 66 79 67 def test_creation_zeros(self): 80 68 81 quantity = Quantity(self. mesh1)82 assert allclose(quantity. vertex_values, [[0.,0.,0.]])83 84 85 quantity = Quantity(self. mesh4)86 assert allclose(quantity.vertex_values, [[0.,0. ,0.], [0.,0.,0.],87 [0.,0. ,0.], [0.,0.,0.]])69 quantity = Quantity(self.domain3) 70 assert allclose(quantity.centroid_values, [[0.,0.,0.]]) 71 72 73 quantity = Quantity(self.domain4) 74 assert allclose(quantity.vertex_values, [[0.,0.], [0.,0.], 75 [0.,0.], [0.,0.]]) 88 76 89 77 90 78 def test_interpolation(self): 91 quantity = Quantity(self.mesh1, [[1,2,3]]) 92 assert allclose(quantity.centroid_values, [2.0]) #Centroid 93 94 assert allclose(quantity.edge_values, [[2.5, 2.0, 1.5]]) 79 quantity = Quantity(self.domain4, self.vertex_values4) 80 assert allclose(quantity.centroid_values, self.centroid_values4) #Centroid 81 95 82 96 83 97 84 def test_interpolation2(self): 98 quantity = Quantity(self.mesh4, 99 [[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]]) 100 assert allclose(quantity.centroid_values, [2., 5., 3., 0.]) #Centroid 101 85 quantity = Quantity(self.domain4, self.vertex_values4) 86 assert allclose(quantity.centroid_values, self.centroid_values4) #Centroid 102 87 103 88 quantity.extrapolate_second_order() 104 89 105 90 #print quantity.vertex_values 106 assert allclose(quantity.vertex_values, [[3.5, -1.0, 3.5], 107 [3.+2./3, 6.+2./3, 4.+2./3], 108 [4.6, 3.4, 1.], 109 [-5.0, 1.0, 4.0]]) 110 111 #print quantity.edge_values 112 assert allclose(quantity.edge_values, [[1.25, 3.5, 1.25], 113 [5. + 2/3.0, 4.0 + 1.0/6, 5.0 + 1.0/6], 114 [2.2, 2.8, 4.0], 115 [2.5, -0.5, -2.0]]) 116 117 118 #assert allclose(quantity.edge_values, [[2.5, 2.0, 1.5], 119 # [5., 5., 5.], 120 # [4.5, 4.5, 0.], 121 # [3.0, -1.5, -1.5]]) 122 123 def test_get_extrema_1(self): 124 quantity = Quantity(self.mesh4, 125 [[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]]) 126 assert allclose(quantity.centroid_values, [2., 5., 3., 0.]) #Centroids 127 128 v = quantity.get_maximum_value() 129 assert v == 5 130 131 v = quantity.get_minimum_value() 132 assert v == 0 133 134 i = quantity.get_maximum_index() 135 assert i == 1 136 137 i = quantity.get_minimum_index() 138 assert i == 3 139 140 x,y = quantity.get_maximum_location() 141 xref, yref = 4.0/3, 4.0/3 142 assert x == xref 143 assert y == yref 144 145 v = quantity.get_values(interpolation_points = [[x,y]]) 146 assert allclose(v, 5) 147 148 149 x,y = quantity.get_minimum_location() 150 v = quantity.get_values(interpolation_points = [[x,y]]) 151 assert allclose(v, 0) 152 153 154 def test_get_maximum_2(self): 155 156 a = [0.0, 0.0] 157 b = [0.0, 2.0] 158 c = [2.0,0.0] 159 d = [0.0, 4.0] 160 e = [2.0, 2.0] 161 f = [4.0,0.0] 162 163 points = [a, b, c, d, e, f] 164 #bac, bce, ecf, dbe 165 vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]] 166 167 domain = Domain(points, vertices) 168 169 quantity = Quantity(domain) 170 quantity.set_values(lambda x, y: x+2*y) #2 4 4 6 171 172 v = quantity.get_maximum_value() 173 assert v == 6 174 175 v = quantity.get_minimum_value() 176 assert v == 2 177 178 i = quantity.get_maximum_index() 179 assert i == 3 180 181 i = quantity.get_minimum_index() 182 assert i == 0 183 184 x,y = quantity.get_maximum_location() 185 xref, yref = 2.0/3, 8.0/3 186 assert x == xref 187 assert y == yref 188 189 v = quantity.get_values(interpolation_points = [[x,y]]) 190 assert allclose(v, 6) 191 192 x,y = quantity.get_minimum_location() 193 v = quantity.get_values(interpolation_points = [[x,y]]) 194 assert allclose(v, 2) 195 196 #Multiple locations for maximum - 197 #Test that the algorithm picks the first occurrence 198 v = quantity.get_maximum_value(indices=[0,1,2]) 199 assert allclose(v, 4) 200 201 i = quantity.get_maximum_index(indices=[0,1,2]) 202 assert i == 1 203 204 x,y = quantity.get_maximum_location(indices=[0,1,2]) 205 xref, yref = 4.0/3, 4.0/3 206 assert x == xref 207 assert y == yref 208 209 v = quantity.get_values(interpolation_points = [[x,y]]) 210 assert allclose(v, 4) 211 212 # More test of indices...... 213 v = quantity.get_maximum_value(indices=[2,3]) 214 assert allclose(v, 6) 215 216 i = quantity.get_maximum_index(indices=[2,3]) 217 assert i == 3 218 219 x,y = quantity.get_maximum_location(indices=[2,3]) 220 xref, yref = 2.0/3, 8.0/3 221 assert x == xref 222 assert y == yref 223 224 v = quantity.get_values(interpolation_points = [[x,y]]) 225 assert allclose(v, 6) 226 227 228 91 assert allclose(quantity.vertex_values,[[ 0.3, 2.7], 92 [ 4.5, 4.5], 93 [ 0.5, 0.5], 94 [ 1.3, 7.7]]) 95 96 97 98 229 99 def test_boundary_allocation(self): 230 quantity = Quantity(self.mesh4, 231 [[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]]) 232 233 assert quantity.boundary_values.shape[0] == len(self.mesh4.boundary) 100 quantity = Quantity(self.domain4, 101 [[1,2], [5,5], [0,9], [-6, 3]]) 102 103 104 assert quantity.boundary_values.shape[0] == len(self.domain4.boundary) 234 105 235 106 236 107 def test_set_values(self): 237 quantity = Quantity(self. mesh4)108 quantity = Quantity(self.domain4) 238 109 239 110 … … 281 152 282 153 def test_set_values_const(self): 283 quantity = Quantity(self. mesh4)154 quantity = Quantity(self.domain4) 284 155 285 156 quantity.set_values(1.0, location = 'vertices') … … 299 170 300 171 def test_set_values_func(self): 301 quantity = Quantity(self. mesh4)172 quantity = Quantity(self.domain4) 302 173 303 174 def f(x, y): … … 320 191 321 192 def test_integral(self): 322 quantity = Quantity(self. mesh4)193 quantity = Quantity(self.domain4) 323 194 324 195 #Try constants first … … 327 198 #print 'Q', quantity.get_integral() 328 199 329 assert allclose(quantity.get_integral(), self. mesh4.get_area() * const)200 assert allclose(quantity.get_integral(), self.domain4.get_area() * const) 330 201 331 202 #Try with a linear function 332 def f(x , y):333 return x +y203 def f(x): 204 return x 334 205 335 206 quantity.set_values(f, location = 'vertices') 336 207 337 338 ref_integral = (4.0/3 + 8.0/3 + 10.0/3 + 10.0/3) * 2 208 209 assert allclose (quantity.centroid_values, 210 [ 0.5, 1.75, 2.75, 4. ]) 211 212 assert allclose (quantity.vertex_values, [[ 0., 1. ], 213 [ 1., 2.5], 214 [ 2.5, 3. ], 215 [ 3., 5. ]]) 216 217 218 ref_integral = 0.5 + 1.5*1.75 + 0.5*2.75 + 2.0*4.0 339 219 340 220 assert allclose (quantity.get_integral(), ref_integral) … … 343 223 344 224 def test_set_vertex_values(self): 345 quantity = Quantity(self. mesh4)225 quantity = Quantity(self.domain4) 346 226 quantity.set_vertex_values([0,1,2,3,4,5]) 347 227 … … 357 237 358 238 def test_set_vertex_values_subset(self): 359 quantity = Quantity(self. mesh4)239 quantity = Quantity(self.domain4) 360 240 quantity.set_vertex_values([0,1,2,3,4,5]) 361 241 quantity.set_vertex_values([0,20,30,50], indices = [0,2,3,5]) … … 366 246 367 247 def test_set_vertex_values_using_general_interface(self): 368 quantity = Quantity(self. mesh4)248 quantity = Quantity(self.domain4) 369 249 370 250 … … 391 271 """ 392 272 393 quantity = Quantity(self. mesh4)273 quantity = Quantity(self.domain4) 394 274 395 275 … … 462 342 """ 463 343 464 quantity = Quantity(self. mesh4)344 quantity = Quantity(self.domain4) 465 345 G = Geo_reference(56, 10, 100) 466 346 quantity.domain.geo_reference = G … … 507 387 508 388 509 def test_set_values_using_fit(self):510 511 512 quantity = Quantity(self.mesh4)513 514 #Get (enough) datapoints515 data_points = [[ 0.66666667, 0.66666667],516 [ 1.33333333, 1.33333333],517 [ 2.66666667, 0.66666667],518 [ 0.66666667, 2.66666667],519 [ 0.0, 1.0],520 [ 0.0, 3.0],521 [ 1.0, 0.0],522 [ 1.0, 1.0],523 [ 1.0, 2.0],524 [ 1.0, 3.0],525 [ 2.0, 1.0],526 [ 3.0, 0.0],527 [ 3.0, 1.0]]528 529 z = linear_function(data_points)530 531 #Use built-in fit_interpolate.fit532 quantity.set_values( Geospatial_data(data_points, z), alpha = 0 )533 #quantity.set_values(points = data_points, values = z, alpha = 0)534 535 536 answer = linear_function(quantity.domain.get_vertex_coordinates())537 #print quantity.vertex_values, answer538 assert allclose(quantity.vertex_values.flat, answer)539 540 541 #Now try by setting the same values directly542 vertex_attributes = fit_to_mesh(data_points,543 quantity.domain.get_nodes(),544 quantity.domain.triangles, #FIXME545 point_attributes=z,546 alpha = 0,547 verbose=False)548 549 #print vertex_attributes550 quantity.set_values(vertex_attributes)551 assert allclose(quantity.vertex_values.flat, answer)552 553 554 555 556 557 def test_test_set_values_using_fit_w_geo(self):558 559 560 #Mesh561 vertex_coordinates = [[0.76, 0.76],562 [0.76, 5.76],563 [5.76, 0.76]]564 triangles = [[0,2,1]]565 566 mesh_georef = Geo_reference(56,-0.76,-0.76)567 mesh1 = Domain(vertex_coordinates, triangles,568 geo_reference = mesh_georef)569 mesh1.check_integrity()570 571 #Quantity572 quantity = Quantity(mesh1)573 574 #Data575 data_points = [[ 201.0, 401.0],576 [ 201.0, 403.0],577 [ 203.0, 401.0]]578 579 z = [2, 4, 4]580 581 data_georef = Geo_reference(56,-200,-400)582 583 584 #Reference585 ref = fit_to_mesh(data_points, vertex_coordinates, triangles,586 point_attributes=z,587 data_origin = data_georef.get_origin(),588 mesh_origin = mesh_georef.get_origin(),589 alpha = 0)590 591 assert allclose( ref, [0,5,5] )592 593 594 #Test set_values595 596 quantity.set_values( Geospatial_data(data_points, z, data_georef), alpha = 0 )597 598 #quantity.set_values(points = data_points,599 # values = z,600 # data_georef = data_georef,601 # alpha = 0)602 603 604 #quantity.set_values(points = data_points,605 # values = z,606 # data_georef = data_georef,607 # alpha = 0)608 assert allclose(quantity.vertex_values.flat, ref)609 610 611 612 #Test set_values using geospatial data object613 quantity.vertex_values[:] = 0.0614 615 geo = Geospatial_data(data_points, z, data_georef)616 617 618 quantity.set_values(geospatial_data = geo, alpha = 0)619 assert allclose(quantity.vertex_values.flat, ref)620 621 622 623 def test_set_values_from_file1(self):624 quantity = Quantity(self.mesh4)625 626 #Get (enough) datapoints627 data_points = [[ 0.66666667, 0.66666667],628 [ 1.33333333, 1.33333333],629 [ 2.66666667, 0.66666667],630 [ 0.66666667, 2.66666667],631 [ 0.0, 1.0],632 [ 0.0, 3.0],633 [ 1.0, 0.0],634 [ 1.0, 1.0],635 [ 1.0, 2.0],636 [ 1.0, 3.0],637 [ 2.0, 1.0],638 [ 3.0, 0.0],639 [ 3.0, 1.0]]640 641 data_geo_spatial = Geospatial_data(data_points,642 geo_reference = Geo_reference(56, 0, 0))643 data_points_absolute = data_geo_spatial.get_data_points(absolute=True)644 attributes = linear_function(data_points_absolute)645 att = 'spam_and_eggs'646 647 #Create .txt file648 ptsfile = tempfile.mktemp(".txt")649 file = open(ptsfile,"w")650 file.write(" x,y," + att + " \n")651 for data_point, attribute in map(None, data_points_absolute652 ,attributes):653 row = str(data_point[0]) + ',' + str(data_point[1]) \654 + ',' + str(attribute)655 file.write(row + "\n")656 file.close()657 658 659 #Check that values can be set from file660 quantity.set_values(filename = ptsfile,661 attribute_name = att, alpha = 0)662 answer = linear_function(quantity.domain.get_vertex_coordinates())663 664 #print quantity.vertex_values.flat665 #print answer666 667 668 assert allclose(quantity.vertex_values.flat, answer)669 670 671 #Check that values can be set from file using default attribute672 quantity.set_values(filename = ptsfile, alpha = 0)673 assert allclose(quantity.vertex_values.flat, answer)674 675 #Cleanup676 import os677 os.remove(ptsfile)678 679 680 681 def Xtest_set_values_from_file_using_polygon(self):682 """test_set_values_from_file_using_polygon(self):683 684 Test that polygon restriction works for general points data685 """686 687 quantity = Quantity(self.mesh4)688 689 #Get (enough) datapoints690 data_points = [[ 0.66666667, 0.66666667],691 [ 1.33333333, 1.33333333],692 [ 2.66666667, 0.66666667],693 [ 0.66666667, 2.66666667],694 [ 0.0, 1.0],695 [ 0.0, 3.0],696 [ 1.0, 0.0],697 [ 1.0, 1.0],698 [ 1.0, 2.0],699 [ 1.0, 3.0],700 [ 2.0, 1.0],701 [ 3.0, 0.0],702 [ 3.0, 1.0]]703 704 data_geo_spatial = Geospatial_data(data_points,705 geo_reference = Geo_reference(56, 0, 0))706 data_points_absolute = data_geo_spatial.get_data_points(absolute=True)707 attributes = linear_function(data_points_absolute)708 att = 'spam_and_eggs'709 710 #Create .txt file711 ptsfile = tempfile.mktemp(".txt")712 file = open(ptsfile,"w")713 file.write(" x,y," + att + " \n")714 for data_point, attribute in map(None, data_points_absolute715 ,attributes):716 row = str(data_point[0]) + ',' + str(data_point[1]) \717 + ',' + str(attribute)718 file.write(row + "\n")719 file.close()720 721 # Create restricting polygon (containing node #4 (2,2) and722 # centroid of triangle #1 (bce)723 polygon = [[1.0, 1.0], [4.0, 1.0],724 [4.0, 4.0], [1.0, 4.0]]725 726 #print self.mesh4.nodes727 #print inside_polygon(self.mesh4.nodes, polygon)728 assert allclose(inside_polygon(self.mesh4.nodes, polygon), 4)729 730 #print quantity.domain.get_vertex_coordinates()731 #print quantity.domain.get_nodes()732 733 # Check that values can be set from file734 quantity.set_values(filename=ptsfile,735 polygon=polygon,736 location='unique vertices',737 alpha=0)738 739 # Get indices for vertex coordinates in polygon740 indices = inside_polygon(quantity.domain.get_vertex_coordinates(),741 polygon)742 points = take(quantity.domain.get_vertex_coordinates(), indices)743 744 answer = linear_function(points)745 746 #print quantity.vertex_values.flat747 #print answer748 749 # Check vertices in polygon have been set750 assert allclose(take(quantity.vertex_values.flat, indices),751 answer)752 753 # Check vertices outside polygon are zero754 indices = outside_polygon(quantity.domain.get_vertex_coordinates(),755 polygon)756 assert allclose(take(quantity.vertex_values.flat, indices),757 0.0)758 759 #Cleanup760 import os761 os.remove(ptsfile)762 763 764 765 766 def test_cache_test_set_values_from_file(self):767 # FIXME (Ole): What is this about?768 # I don't think it checks anything new769 quantity = Quantity(self.mesh4)770 771 #Get (enough) datapoints772 data_points = [[ 0.66666667, 0.66666667],773 [ 1.33333333, 1.33333333],774 [ 2.66666667, 0.66666667],775 [ 0.66666667, 2.66666667],776 [ 0.0, 1.0],777 [ 0.0, 3.0],778 [ 1.0, 0.0],779 [ 1.0, 1.0],780 [ 1.0, 2.0],781 [ 1.0, 3.0],782 [ 2.0, 1.0],783 [ 3.0, 0.0],784 [ 3.0, 1.0]]785 786 georef = Geo_reference(56, 0, 0)787 data_geo_spatial = Geospatial_data(data_points,788 geo_reference=georef)789 790 data_points_absolute = data_geo_spatial.get_data_points(absolute=True)791 attributes = linear_function(data_points_absolute)792 att = 'spam_and_eggs'793 794 # Create .txt file795 ptsfile = tempfile.mktemp(".txt")796 file = open(ptsfile,"w")797 file.write(" x,y," + att + " \n")798 for data_point, attribute in map(None, data_points_absolute799 ,attributes):800 row = str(data_point[0]) + ',' + str(data_point[1]) \801 + ',' + str(attribute)802 file.write(row + "\n")803 file.close()804 805 806 # Check that values can be set from file807 quantity.set_values(filename=ptsfile,808 attribute_name=att,809 alpha=0,810 use_cache=True,811 verbose=False)812 answer = linear_function(quantity.domain.get_vertex_coordinates())813 assert allclose(quantity.vertex_values.flat, answer)814 815 816 # Check that values can be set from file using default attribute817 quantity.set_values(filename=ptsfile,818 alpha=0)819 assert allclose(quantity.vertex_values.flat, answer)820 821 # Check cache822 quantity.set_values(filename=ptsfile,823 attribute_name=att,824 alpha=0,825 use_cache=True,826 verbose=False)827 828 829 #Cleanup830 import os831 os.remove(ptsfile)832 833 def test_set_values_from_lat_long(self):834 quantity = Quantity(self.mesh_onslow)835 836 #Get (enough) datapoints837 data_points = [[-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65],838 [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6]]839 840 data_geo_spatial = Geospatial_data(data_points,841 points_are_lats_longs=True)842 points_UTM = data_geo_spatial.get_data_points(absolute=True)843 attributes = linear_function(points_UTM)844 att = 'elevation'845 846 #Create .txt file847 txt_file = tempfile.mktemp(".txt")848 file = open(txt_file,"w")849 file.write(" lat,long," + att + " \n")850 for data_point, attribute in map(None, data_points, attributes):851 row = str(data_point[0]) + ',' + str(data_point[1]) \852 + ',' + str(attribute)853 #print "row", row854 file.write(row + "\n")855 file.close()856 857 858 #Check that values can be set from file859 quantity.set_values(filename=txt_file,860 attribute_name=att,861 alpha=0)862 answer = linear_function(quantity.domain.get_vertex_coordinates())863 864 #print "quantity.vertex_values.flat", quantity.vertex_values.flat865 #print "answer",answer866 867 assert allclose(quantity.vertex_values.flat, answer)868 869 870 #Check that values can be set from file using default attribute871 quantity.set_values(filename=txt_file, alpha=0)872 assert allclose(quantity.vertex_values.flat, answer)873 874 #Cleanup875 import os876 os.remove(txt_file)877 878 def test_set_values_from_lat_long(self):879 quantity = Quantity(self.mesh_onslow)880 881 #Get (enough) datapoints882 data_points = [[-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65],883 [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6]]884 885 data_geo_spatial = Geospatial_data(data_points,886 points_are_lats_longs=True)887 points_UTM = data_geo_spatial.get_data_points(absolute=True)888 attributes = linear_function(points_UTM)889 att = 'elevation'890 891 #Create .txt file892 txt_file = tempfile.mktemp(".txt")893 file = open(txt_file,"w")894 file.write(" lat,long," + att + " \n")895 for data_point, attribute in map(None, data_points, attributes):896 row = str(data_point[0]) + ',' + str(data_point[1]) \897 + ',' + str(attribute)898 #print "row", row899 file.write(row + "\n")900 file.close()901 902 903 #Check that values can be set from file904 quantity.set_values(filename=txt_file,905 attribute_name=att, alpha=0)906 answer = linear_function(quantity.domain.get_vertex_coordinates())907 908 #print "quantity.vertex_values.flat", quantity.vertex_values.flat909 #print "answer",answer910 911 assert allclose(quantity.vertex_values.flat, answer)912 913 914 #Check that values can be set from file using default attribute915 quantity.set_values(filename=txt_file, alpha=0)916 assert allclose(quantity.vertex_values.flat, answer)917 918 #Cleanup919 import os920 os.remove(txt_file)921 922 def test_set_values_from_UTM_pts(self):923 quantity = Quantity(self.mesh_onslow)924 925 #Get (enough) datapoints926 data_points = [[-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65],927 [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6]]928 929 data_geo_spatial = Geospatial_data(data_points,930 points_are_lats_longs=True)931 points_UTM = data_geo_spatial.get_data_points(absolute=True)932 attributes = linear_function(points_UTM)933 att = 'elevation'934 935 #Create .txt file936 txt_file = tempfile.mktemp(".txt")937 file = open(txt_file,"w")938 file.write(" x,y," + att + " \n")939 for data_point, attribute in map(None, points_UTM, attributes):940 row = str(data_point[0]) + ',' + str(data_point[1]) \941 + ',' + str(attribute)942 #print "row", row943 file.write(row + "\n")944 file.close()945 946 947 pts_file = tempfile.mktemp(".pts")948 convert = Geospatial_data(txt_file)949 convert.export_points_file(pts_file)950 951 #Check that values can be set from file952 quantity.set_values_from_file(pts_file, att, 0,953 'vertices', None)954 answer = linear_function(quantity.domain.get_vertex_coordinates())955 #print "quantity.vertex_values.flat", quantity.vertex_values.flat956 #print "answer",answer957 assert allclose(quantity.vertex_values.flat, answer)958 959 #Check that values can be set from file960 quantity.set_values(filename=pts_file,961 attribute_name=att, alpha=0)962 answer = linear_function(quantity.domain.get_vertex_coordinates())963 #print "quantity.vertex_values.flat", quantity.vertex_values.flat964 #print "answer",answer965 assert allclose(quantity.vertex_values.flat, answer)966 967 968 #Check that values can be set from file using default attribute969 quantity.set_values(filename=txt_file, alpha=0)970 assert allclose(quantity.vertex_values.flat, answer)971 972 #Cleanup973 import os974 os.remove(txt_file)975 os.remove(pts_file)976 977 def verbose_test_set_values_from_UTM_pts(self):978 quantity = Quantity(self.mesh_onslow)979 980 #Get (enough) datapoints981 data_points = [[-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65],982 [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6],983 [-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65],984 [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6],985 [-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65],986 [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6],987 [-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65],988 [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6],989 [-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65],990 [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6],991 [-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65],992 [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6],993 [-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65],994 [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6],995 [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6],996 [-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65],997 [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6],998 [-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65],999 [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6],1000 [-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65],1001 [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6],1002 [-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65],1003 [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6],1004 [-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65],1005 [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6],1006 [-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65],1007 [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6],1008 ]1009 1010 data_geo_spatial = Geospatial_data(data_points,1011 points_are_lats_longs=True)1012 points_UTM = data_geo_spatial.get_data_points(absolute=True)1013 attributes = linear_function(points_UTM)1014 att = 'elevation'1015 1016 #Create .txt file1017 txt_file = tempfile.mktemp(".txt")1018 file = open(txt_file,"w")1019 file.write(" x,y," + att + " \n")1020 for data_point, attribute in map(None, points_UTM, attributes):1021 row = str(data_point[0]) + ',' + str(data_point[1]) \1022 + ',' + str(attribute)1023 #print "row", row1024 file.write(row + "\n")1025 file.close()1026 1027 1028 pts_file = tempfile.mktemp(".pts")1029 convert = Geospatial_data(txt_file)1030 convert.export_points_file(pts_file)1031 1032 #Check that values can be set from file1033 quantity.set_values_from_file(pts_file, att, 0,1034 'vertices', None, verbose = True,1035 max_read_lines=2)1036 answer = linear_function(quantity.domain.get_vertex_coordinates())1037 #print "quantity.vertex_values.flat", quantity.vertex_values.flat1038 #print "answer",answer1039 assert allclose(quantity.vertex_values.flat, answer)1040 1041 #Check that values can be set from file1042 quantity.set_values(filename=pts_file,1043 attribute_name=att, alpha=0)1044 answer = linear_function(quantity.domain.get_vertex_coordinates())1045 #print "quantity.vertex_values.flat", quantity.vertex_values.flat1046 #print "answer",answer1047 assert allclose(quantity.vertex_values.flat, answer)1048 1049 1050 #Check that values can be set from file using default attribute1051 quantity.set_values(filename=txt_file, alpha=0)1052 assert allclose(quantity.vertex_values.flat, answer)1053 1054 #Cleanup1055 import os1056 os.remove(txt_file)1057 os.remove(pts_file)1058 1059 def test_set_values_from_file_with_georef1(self):1060 1061 #Mesh in zone 56 (absolute coords)1062 1063 x0 = 314036.587279821064 y0 = 6224951.29600921065 1066 a = [x0+0.0, y0+0.0]1067 b = [x0+0.0, y0+2.0]1068 c = [x0+2.0, y0+0.0]1069 d = [x0+0.0, y0+4.0]1070 e = [x0+2.0, y0+2.0]1071 f = [x0+4.0, y0+0.0]1072 1073 points = [a, b, c, d, e, f]1074 1075 #bac, bce, ecf, dbe1076 elements = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]1077 1078 #absolute going in ..1079 mesh4 = Domain(points, elements,1080 geo_reference = Geo_reference(56, 0, 0))1081 mesh4.check_integrity()1082 quantity = Quantity(mesh4)1083 1084 #Get (enough) datapoints (relative to georef)1085 data_points_rel = [[ 0.66666667, 0.66666667],1086 [ 1.33333333, 1.33333333],1087 [ 2.66666667, 0.66666667],1088 [ 0.66666667, 2.66666667],1089 [ 0.0, 1.0],1090 [ 0.0, 3.0],1091 [ 1.0, 0.0],1092 [ 1.0, 1.0],1093 [ 1.0, 2.0],1094 [ 1.0, 3.0],1095 [ 2.0, 1.0],1096 [ 3.0, 0.0],1097 [ 3.0, 1.0]]1098 1099 data_geo_spatial = Geospatial_data(data_points_rel,1100 geo_reference = Geo_reference(56, x0, y0))1101 data_points_absolute = data_geo_spatial.get_data_points(absolute=True)1102 attributes = linear_function(data_points_absolute)1103 att = 'spam_and_eggs'1104 1105 #Create .txt file1106 ptsfile = tempfile.mktemp(".txt")1107 file = open(ptsfile,"w")1108 file.write(" x,y," + att + " \n")1109 for data_point, attribute in map(None, data_points_absolute1110 ,attributes):1111 row = str(data_point[0]) + ',' + str(data_point[1]) \1112 + ',' + str(attribute)1113 file.write(row + "\n")1114 file.close()1115 1116 #file = open(ptsfile, 'r')1117 #lines = file.readlines()1118 #file.close()1119 1120 1121 #Check that values can be set from file1122 quantity.set_values(filename=ptsfile,1123 attribute_name=att, alpha=0)1124 answer = linear_function(quantity.domain.get_vertex_coordinates())1125 1126 assert allclose(quantity.vertex_values.flat, answer)1127 1128 1129 #Check that values can be set from file using default attribute1130 quantity.set_values(filename=ptsfile, alpha=0)1131 assert allclose(quantity.vertex_values.flat, answer)1132 1133 #Cleanup1134 import os1135 os.remove(ptsfile)1136 1137 1138 def test_set_values_from_file_with_georef2(self):1139 1140 #Mesh in zone 56 (relative coords)1141 1142 x0 = 314036.587279821143 y0 = 6224951.29600921144 #x0 = 0.01145 #y0 = 0.01146 1147 a = [0.0, 0.0]1148 b = [0.0, 2.0]1149 c = [2.0, 0.0]1150 d = [0.0, 4.0]1151 e = [2.0, 2.0]1152 f = [4.0, 0.0]1153 1154 points = [a, b, c, d, e, f]1155 1156 #bac, bce, ecf, dbe1157 elements = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]1158 1159 mesh4 = Domain(points, elements,1160 geo_reference = Geo_reference(56, x0, y0))1161 mesh4.check_integrity()1162 quantity = Quantity(mesh4)1163 1164 #Get (enough) datapoints1165 data_points = [[ x0+0.66666667, y0+0.66666667],1166 [ x0+1.33333333, y0+1.33333333],1167 [ x0+2.66666667, y0+0.66666667],1168 [ x0+0.66666667, y0+2.66666667],1169 [ x0+0.0, y0+1.0],1170 [ x0+0.0, y0+3.0],1171 [ x0+1.0, y0+0.0],1172 [ x0+1.0, y0+1.0],1173 [ x0+1.0, y0+2.0],1174 [ x0+1.0, y0+3.0],1175 [ x0+2.0, y0+1.0],1176 [ x0+3.0, y0+0.0],1177 [ x0+3.0, y0+1.0]]1178 1179 1180 data_geo_spatial = Geospatial_data(data_points,1181 geo_reference = Geo_reference(56, 0, 0))1182 data_points_absolute = data_geo_spatial.get_data_points(absolute=True)1183 attributes = linear_function(data_points_absolute)1184 att = 'spam_and_eggs'1185 1186 #Create .txt file1187 ptsfile = tempfile.mktemp(".txt")1188 file = open(ptsfile,"w")1189 file.write(" x,y," + att + " \n")1190 for data_point, attribute in map(None, data_points_absolute1191 ,attributes):1192 row = str(data_point[0]) + ',' + str(data_point[1]) \1193 + ',' + str(attribute)1194 file.write(row + "\n")1195 file.close()1196 1197 1198 #Check that values can be set from file1199 quantity.set_values(filename=ptsfile,1200 attribute_name=att, alpha=0)1201 answer = linear_function(quantity.domain. \1202 get_vertex_coordinates(absolute=True))1203 1204 1205 assert allclose(quantity.vertex_values.flat, answer)1206 1207 1208 #Check that values can be set from file using default attribute1209 quantity.set_values(filename=ptsfile, alpha=0)1210 assert allclose(quantity.vertex_values.flat, answer)1211 1212 #Cleanup1213 import os1214 os.remove(ptsfile)1215 1216 1217 1218 389 1219 390 def test_set_values_from_quantity(self): 1220 391 1221 quantity1 = Quantity(self. mesh4)392 quantity1 = Quantity(self.domain4) 1222 393 quantity1.set_vertex_values([0,1,2,3,4,5]) 1223 394 … … 1226 397 1227 398 1228 quantity2 = Quantity(self. mesh4)399 quantity2 = Quantity(self.domain4) 1229 400 quantity2.set_values(quantity=quantity1) 1230 401 assert allclose(quantity2.vertex_values, … … 1247 418 1248 419 1249 def Xtest_set_values_from_quantity_using_polygon(self): 1250 """test_set_values_from_quantity_using_polygon(self): 1251 1252 Check that polygon can be used to restrict set_values when 1253 using another quantity as argument. 1254 """ 1255 1256 # Create restricting polygon (containing node #4 (2,2) and 1257 # centroid of triangle #1 (bce) 1258 polygon = [[1.0, 1.0], [4.0, 1.0], 1259 [4.0, 4.0], [1.0, 4.0]] 1260 assert allclose(inside_polygon(self.mesh4.nodes, polygon), 4) 1261 1262 quantity1 = Quantity(self.mesh4) 1263 quantity1.set_vertex_values([0,1,2,3,4,5]) 420 421 422 def test_overloading(self): 423 424 quantity1 = Quantity(self.domain4) 425 quantity1.set_values( [[0,1],[1,2],[2,3],[3,4]], 426 location = 'vertices') 1264 427 1265 428 assert allclose(quantity1.vertex_values, 1266 [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]) 1267 1268 1269 quantity2 = Quantity(self.mesh4) 1270 quantity2.set_values(quantity=quantity1, 1271 polygon=polygon) 1272 1273 msg = 'Only node #4(e) at (2,2) should have values applied ' 1274 assert allclose(quantity2.vertex_values, 1275 [[0,0,0], [0,0,4], [4,0,0], [0,0,4]]), msg 1276 #bac, bce, ecf, dbe 1277 1278 1279 1280 def test_overloading(self): 1281 1282 quantity1 = Quantity(self.mesh4) 1283 quantity1.set_vertex_values([0,1,2,3,4,5]) 1284 1285 assert allclose(quantity1.vertex_values, 1286 [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]) 1287 1288 1289 quantity2 = Quantity(self.mesh4) 1290 quantity2.set_values([[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]], 429 [[0,1], [1,2], [2,3], [3,4]]) 430 431 432 quantity2 = Quantity(self.domain4) 433 quantity2.set_values([[1,2], [5,5], [0,9], [-6, 3]], 1291 434 location = 'vertices') 1292 435 1293 436 1294 437 1295 quantity3 = Quantity(self. mesh4)1296 quantity3.set_values([[2,2 ,2], [7,8,9], [7,6,3], [3, 8, -8]],438 quantity3 = Quantity(self.domain4) 439 quantity3.set_values([[2,2], [7,8], [7,6], [3, -8]], 1297 440 location = 'vertices') 1298 441 … … 1302 445 assert allclose(Q.vertex_values, -quantity1.vertex_values) 1303 446 assert allclose(Q.centroid_values, -quantity1.centroid_values) 1304 assert allclose(Q.edge_values, -quantity1.edge_values) 447 1305 448 1306 449 # Addition … … 1308 451 assert allclose(Q.vertex_values, quantity1.vertex_values + 7) 1309 452 assert allclose(Q.centroid_values, quantity1.centroid_values + 7) 1310 assert allclose(Q.edge_values, quantity1.edge_values + 7)1311 453 1312 454 Q = 7 + quantity1 1313 455 assert allclose(Q.vertex_values, quantity1.vertex_values + 7) 1314 456 assert allclose(Q.centroid_values, quantity1.centroid_values + 7) 1315 assert allclose(Q.edge_values, quantity1.edge_values + 7)1316 457 1317 458 Q = quantity1 + quantity2 … … 1320 461 assert allclose(Q.centroid_values, 1321 462 quantity1.centroid_values + quantity2.centroid_values) 1322 assert allclose(Q.edge_values,1323 quantity1.edge_values + quantity2.edge_values)1324 1325 463 1326 464 Q = quantity1 + quantity2 - 3 … … 1336 474 assert allclose(Q.vertex_values, quantity1.vertex_values*3) 1337 475 assert allclose(Q.centroid_values, quantity1.centroid_values*3) 1338 assert allclose(Q.edge_values, quantity1.edge_values*3) 476 1339 477 Q = 3*quantity1 1340 478 assert allclose(Q.vertex_values, quantity1.vertex_values*3) … … 1342 480 #Multiplication 1343 481 Q = quantity1 * quantity2 1344 #print Q.vertex_values1345 #print Q.centroid_values1346 #print quantity1.centroid_values1347 #print quantity2.centroid_values1348 1349 482 assert allclose(Q.vertex_values, 1350 483 quantity1.vertex_values * quantity2.vertex_values) … … 1403 536 quantity2.vertex_values**2)**0.5) 1404 537 1405 1406 1407 1408 1409 1410 1411 538 def test_compute_gradient(self): 1412 quantity = Quantity(self. mesh4)539 quantity = Quantity(self.domain6) 1413 540 1414 541 #Set up for a gradient of (2,0) at mid triangle 1415 quantity.set_values([2.0, 4.0, 6.0,2.0],542 quantity.set_values([2.0, 4.0, 4.0, 5.0, 10.0, 12.0], 1416 543 location = 'centroids') 1417 544 … … 1420 547 quantity.compute_gradients() 1421 548 1422 a = quantity.x_gradient 1423 b = quantity.y_gradient 1424 #print self.mesh4.centroid_coordinates 1425 #print a, b 1426 1427 #The central triangle (1) 1428 #(using standard gradient based on neigbours controid values) 1429 assert allclose(a[1], 2.0) 1430 assert allclose(b[1], 0.0) 1431 1432 1433 #Left triangle (0) using two point gradient 1434 #q0 = q1 + a*(x0-x1) + b*(y0-y1) <=> 1435 #2 = 4 + a*(-2/3) + b*(-2/3) 1436 assert allclose(a[0] + b[0], 3) 1437 #From orthogonality (a*(y0-y1) + b*(x0-x1) == 0) 1438 assert allclose(a[0] - b[0], 0) 1439 1440 1441 #Right triangle (2) using two point gradient 1442 #q2 = q1 + a*(x2-x1) + b*(y2-y1) <=> 1443 #6 = 4 + a*(4/3) + b*(-2/3) 1444 assert allclose(2*a[2] - b[2], 3) 1445 #From orthogonality (a*(y1-y2) + b*(x2-x1) == 0) 1446 assert allclose(a[2] + 2*b[2], 0) 1447 1448 1449 #Top triangle (3) using two point gradient 1450 #q3 = q1 + a*(x3-x1) + b*(y3-y1) <=> 1451 #2 = 4 + a*(-2/3) + b*(4/3) 1452 assert allclose(a[3] - 2*b[3], 3) 1453 #From orthogonality (a*(y1-y3) + b*(x3-x1) == 0) 1454 assert allclose(2*a[3] + b[3], 0) 1455 1456 1457 1458 #print a, b 549 a = quantity.gradients 550 551 assert allclose(a, [ 3., 1., 0.5, 3., 3.5, 0.5]) 552 1459 553 quantity.extrapolate_second_order() 1460 554 1461 #Apply q(x,y) = qc + a*(x-xc) + b*(y-yc) 1462 assert allclose(quantity.vertex_values[0,:], [3., 0., 3.]) 1463 assert allclose(quantity.vertex_values[1,:], [4./3, 16./3, 16./3]) 1464 1465 1466 #a = 1.2, b=-0.6 1467 #q(4,0) = 6 + a*(4 - 8/3) + b*(-2/3) 1468 assert allclose(quantity.vertex_values[2,2], 8) 1469 1470 def test_get_gradients(self): 1471 quantity = Quantity(self.mesh4) 1472 1473 #Set up for a gradient of (2,0) at mid triangle 1474 quantity.set_values([2.0, 4.0, 6.0, 2.0], 1475 location = 'centroids') 1476 1477 1478 #Gradients 1479 quantity.compute_gradients() 1480 1481 a, b = quantity.get_gradients() 1482 #print self.mesh4.centroid_coordinates 1483 #print a, b 1484 1485 #The central triangle (1) 1486 #(using standard gradient based on neigbours controid values) 1487 assert allclose(a[1], 2.0) 1488 assert allclose(b[1], 0.0) 1489 1490 1491 #Left triangle (0) using two point gradient 1492 #q0 = q1 + a*(x0-x1) + b*(y0-y1) <=> 1493 #2 = 4 + a*(-2/3) + b*(-2/3) 1494 assert allclose(a[0] + b[0], 3) 1495 #From orthogonality (a*(y0-y1) + b*(x0-x1) == 0) 1496 assert allclose(a[0] - b[0], 0) 1497 1498 1499 #Right triangle (2) using two point gradient 1500 #q2 = q1 + a*(x2-x1) + b*(y2-y1) <=> 1501 #6 = 4 + a*(4/3) + b*(-2/3) 1502 assert allclose(2*a[2] - b[2], 3) 1503 #From orthogonality (a*(y1-y2) + b*(x2-x1) == 0) 1504 assert allclose(a[2] + 2*b[2], 0) 1505 1506 1507 #Top triangle (3) using two point gradient 1508 #q3 = q1 + a*(x3-x1) + b*(y3-y1) <=> 1509 #2 = 4 + a*(-2/3) + b*(4/3) 1510 assert allclose(a[3] - 2*b[3], 3) 1511 #From orthogonality (a*(y1-y3) + b*(x3-x1) == 0) 1512 assert allclose(2*a[3] + b[3], 0) 1513 555 556 assert allclose(quantity.vertex_values, [[ 1., 3. ], 557 [ 4., 4. ], 558 [ 4., 4. ], 559 [ 4., 6.], 560 [ 8.25, 11.75], 561 [ 11., 13. ]]) 562 563 1514 564 1515 565 def test_second_order_extrapolation2(self): 1516 quantity = Quantity(self. mesh4)566 quantity = Quantity(self.domain4) 1517 567 1518 568 #Set up for a gradient of (3,1), f(x) = 3x+y … … 1543 593 1544 594 def test_backup_saxpy_centroid_values(self): 1545 quantity = Quantity(self. mesh4)595 quantity = Quantity(self.domain4) 1546 596 1547 597 #Set up for a gradient of (3,1), f(x) = 3x+y … … 1566 616 1567 617 def test_first_order_extrapolator(self): 1568 quantity = Quantity(self. mesh4)1569 1570 #Test centroids1571 quantity.set_values( [1.,2.,3.,4.], location = 'centroids')1572 assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid618 quantity = Quantity(self.domain4) 619 620 centroid_values = array([1.,2.,3.,4.]) 621 quantity.set_values(centroid_values, location = 'centroids') 622 assert allclose(quantity.centroid_values, centroid_values) #Centroid 1573 623 1574 624 #Extrapolate … … 1576 626 1577 627 #Check that gradient is zero 1578 a ,b = quantity.get_gradients()628 a = quantity.gradients 1579 629 assert allclose(a, [0,0,0,0]) 1580 assert allclose(b, [0,0,0,0])630 1581 631 1582 632 #Check vertices but not edge values 1583 633 assert allclose(quantity.vertex_values, 1584 [[1,1 ,1], [2,2,2], [3,3,3], [4, 4,4]])634 [[1,1], [2,2], [3,3], [4,4]]) 1585 635 1586 636 1587 637 def test_second_order_extrapolator(self): 1588 quantity = Quantity(self. mesh4)638 quantity = Quantity(self.domain4) 1589 639 1590 640 #Set up for a gradient of (3,0) at mid triangle … … 1616 666 1617 667 1618 1619 1620 1621 def test_limit_vertices_by_all_neighbours(self): 1622 quantity = Quantity(self.mesh4) 668 def test_limit(self): 669 quantity = Quantity(self.domain4) 1623 670 1624 671 #Create a deliberate overshoot (e.g. from gradient computation) 1625 quantity.set_values([[3,0,3], [2,2,6], [5,3,8], [8,3,5]]) 1626 672 quantity.set_values([[0,0], [2,20], [-20,3], [8,3]]) 1627 673 1628 674 #Limit 1629 quantity.limit_vertices_by_all_neighbours() 1630 1631 #Assert that central triangle is limited by neighbours 1632 assert quantity.vertex_values[1,0] >= quantity.vertex_values[0,0] 1633 assert quantity.vertex_values[1,0] <= quantity.vertex_values[3,1] 1634 1635 assert quantity.vertex_values[1,1] <= quantity.vertex_values[2,1] 1636 assert quantity.vertex_values[1,1] >= quantity.vertex_values[0,2] 1637 1638 assert quantity.vertex_values[1,2] <= quantity.vertex_values[2,0] 1639 assert quantity.vertex_values[1,2] <= quantity.vertex_values[3,1] 1640 1641 1642 1643 #Assert that quantities are conserved 1644 from Numeric import sum 1645 for k in range(quantity.centroid_values.shape[0]): 1646 assert allclose (quantity.centroid_values[k], 1647 sum(quantity.vertex_values[k,:])/3) 1648 1649 1650 1651 def test_limit_edges_by_all_neighbours(self): 1652 quantity = Quantity(self.mesh4) 1653 1654 #Create a deliberate overshoot (e.g. from gradient computation) 1655 quantity.set_values([[3,0,3], [2,2,6], [5,3,8], [8,3,5]]) 1656 1657 1658 #Limit 1659 quantity.limit_edges_by_all_neighbours() 1660 1661 #Assert that central triangle is limited by neighbours 1662 assert quantity.edge_values[1,0] <= quantity.centroid_values[2] 1663 assert quantity.edge_values[1,0] >= quantity.centroid_values[0] 1664 1665 assert quantity.edge_values[1,1] <= quantity.centroid_values[2] 1666 assert quantity.edge_values[1,1] >= quantity.centroid_values[0] 1667 1668 assert quantity.edge_values[1,2] <= quantity.centroid_values[2] 1669 assert quantity.edge_values[1,2] >= quantity.centroid_values[0] 1670 1671 1672 1673 #Assert that quantities are conserved 1674 from Numeric import sum 1675 for k in range(quantity.centroid_values.shape[0]): 1676 assert allclose (quantity.centroid_values[k], 1677 sum(quantity.vertex_values[k,:])/3) 1678 1679 1680 def test_limit_edges_by_neighbour(self): 1681 quantity = Quantity(self.mesh4) 1682 1683 #Create a deliberate overshoot (e.g. from gradient computation) 1684 quantity.set_values([[3,0,3], [2,2,6], [5,3,8], [8,3,5]]) 1685 1686 1687 #Limit 1688 quantity.limit_edges_by_neighbour() 1689 1690 #Assert that central triangle is limited by neighbours 1691 assert quantity.edge_values[1,0] <= quantity.centroid_values[3] 1692 assert quantity.edge_values[1,0] >= quantity.centroid_values[1] 1693 1694 assert quantity.edge_values[1,1] <= quantity.centroid_values[2] 1695 assert quantity.edge_values[1,1] >= quantity.centroid_values[1] 1696 1697 assert quantity.edge_values[1,2] <= quantity.centroid_values[1] 1698 assert quantity.edge_values[1,2] >= quantity.centroid_values[0] 1699 1700 1701 1702 #Assert that quantities are conserved 1703 from Numeric import sum 1704 for k in range(quantity.centroid_values.shape[0]): 1705 assert allclose (quantity.centroid_values[k], 1706 sum(quantity.vertex_values[k,:])/3) 1707 1708 def test_limiter2(self): 1709 """Taken from test_shallow_water 1710 """ 1711 quantity = Quantity(self.mesh4) 1712 quantity.domain.beta_w = 0.9 1713 675 quantity.limit_minmod() 676 677 678 #cells 1 and 2 are local max and min 679 assert quantity.vertex_values[1][0] == quantity.centroid_values[1] 680 assert quantity.vertex_values[1][1] == quantity.centroid_values[1] 681 682 assert quantity.vertex_values[2][0] == quantity.centroid_values[2] 683 assert quantity.vertex_values[2][1] == quantity.centroid_values[2] 684 685 686 687 def test_distribute_first_order(self): 688 quantity = Quantity(self.domain4) 689 1714 690 #Test centroids 1715 quantity.set_values([2.,4.,8.,2.], location = 'centroids') 1716 assert allclose(quantity.centroid_values, [2, 4, 8, 2]) #Centroid 691 centroid_values = array([1.,2.,3.,4.]) 692 quantity.set_values(centroid_values, location = 'centroids') 693 assert allclose(quantity.centroid_values, centroid_values) #Centroid 694 695 696 #Extrapolate from centroid to vertices and edges 697 quantity.extrapolate_first_order() 698 699 assert allclose(quantity.vertex_values,[[ 1., 1.], 700 [ 2., 2.], 701 [ 3., 3.], 702 [ 4., 4.]]) 703 704 705 706 def test_distribute_second_order(self): 707 quantity = Quantity(self.domain4) 708 709 #Test centroids 710 centroid_values = array([2.,4.,8.,2.]) 711 quantity.set_values(centroid_values, location = 'centroids') 712 assert allclose(quantity.centroid_values, centroid_values) #Centroid 1717 713 1718 714 … … 1720 716 quantity.extrapolate_second_order() 1721 717 1722 assert allclose(quantity.vertex_values[1,:], [0.0, 6, 6]) 1723 1724 #Limit 1725 quantity.limit() 1726 1727 # limited value for beta_w = 0.9 1728 1729 assert allclose(quantity.vertex_values[1,:], [2.2, 4.9, 4.9]) 1730 # limited values for beta_w = 0.5 1731 #assert allclose(quantity.vertex_values[1,:], [3.0, 4.5, 4.5]) 1732 1733 1734 #Assert that quantities are conserved 1735 from Numeric import sum 1736 for k in range(quantity.centroid_values.shape[0]): 1737 assert allclose (quantity.centroid_values[k], 1738 sum(quantity.vertex_values[k,:])/3) 1739 1740 1741 1742 1743 1744 def test_distribute_first_order(self): 1745 quantity = Quantity(self.mesh4) 718 assert allclose(quantity.vertex_values, [[ 1.2, 2.8], 719 [ 2., 6. ], 720 [ 8., 8. ], 721 [ 6.8, -2.8]]) 722 723 724 def test_update_explicit(self): 725 quantity = Quantity(self.domain4) 1746 726 1747 727 #Test centroids … … 1749 729 assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid 1750 730 1751 1752 #Extrapolate from centroid to vertices and edges 1753 quantity.extrapolate_first_order() 1754 1755 #Interpolate 1756 #quantity.interpolate_from_vertices_to_edges() 1757 1758 assert allclose(quantity.vertex_values, 1759 [[1,1,1], [2,2,2], [3,3,3], [4, 4, 4]]) 1760 assert allclose(quantity.edge_values, [[1,1,1], [2,2,2], 1761 [3,3,3], [4, 4, 4]]) 1762 1763 1764 def test_interpolate_from_vertices_to_edges(self): 1765 quantity = Quantity(self.mesh4) 1766 1767 quantity.vertex_values = array([[1,0,2], [1,2,4], [4,2,5], [3,1,4]],Float) 1768 1769 quantity.interpolate_from_vertices_to_edges() 1770 1771 assert allclose(quantity.edge_values, [[1., 1.5, 0.5], 1772 [3., 2.5, 1.5], 1773 [3.5, 4.5, 3.], 1774 [2.5, 3.5, 2]]) 1775 1776 1777 def test_interpolate_from_edges_to_vertices(self): 1778 quantity = Quantity(self.mesh4) 1779 1780 quantity.edge_values = array([[1., 1.5, 0.5], 1781 [3., 2.5, 1.5], 1782 [3.5, 4.5, 3.], 1783 [2.5, 3.5, 2]],Float) 1784 1785 quantity.interpolate_from_edges_to_vertices() 1786 1787 assert allclose(quantity.vertex_values, 1788 [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]) 1789 1790 1791 1792 def test_distribute_second_order(self): 1793 quantity = Quantity(self.mesh4) 1794 1795 #Test centroids 1796 quantity.set_values([2.,4.,8.,2.], location = 'centroids') 1797 assert allclose(quantity.centroid_values, [2, 4, 8, 2]) #Centroid 1798 1799 1800 #Extrapolate 1801 quantity.extrapolate_second_order() 1802 1803 assert allclose(quantity.vertex_values[1,:], [0.0, 6, 6]) 1804 1805 1806 def test_update_explicit(self): 1807 quantity = Quantity(self.mesh4) 731 #Set explicit_update 732 quantity.explicit_update = array( [1.,1.,1.,1.] ) 733 734 #Update with given timestep 735 quantity.update(0.1) 736 737 x = array([1, 2, 3, 4]) + array( [.1,.1,.1,.1] ) 738 assert allclose( quantity.centroid_values, x) 739 740 def test_update_semi_implicit(self): 741 quantity = Quantity(self.domain4) 1808 742 1809 743 #Test centroids … … 1811 745 assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid 1812 746 1813 #Set explicit_update1814 quantity.explicit_update = array( [1.,1.,1.,1.] )1815 1816 #Update with given timestep1817 quantity.update(0.1)1818 1819 x = array([1, 2, 3, 4]) + array( [.1,.1,.1,.1] )1820 assert allclose( quantity.centroid_values, x)1821 1822 def test_update_semi_implicit(self):1823 quantity = Quantity(self.mesh4)1824 1825 #Test centroids1826 quantity.set_values([1.,2.,3.,4.], location = 'centroids')1827 assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid1828 1829 747 #Set semi implicit update 1830 748 quantity.semi_implicit_update = array([1.,1.,1.,1.]) … … 1842 760 1843 761 def test_both_updates(self): 1844 quantity = Quantity(self. mesh4)762 quantity = Quantity(self.domain4) 1845 763 1846 764 #Test centroids 1847 quantity.set_values([1.,2.,3.,4.], location = 'centroids') 1848 assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid 765 centroid_values = array( [1, 2, 3, 4] ) 766 quantity.set_values(centroid_values, location = 'centroids') 767 assert allclose(quantity.centroid_values, centroid_values) #Centroid 1849 768 1850 769 #Set explicit_update 1851 quantity.explicit_update = array( [4.,3.,2.,1.] ) 770 explicit_update = array( [4.,3.,2.,1.] ) 771 quantity.explicit_update[:] = explicit_update 1852 772 1853 773 #Set semi implicit update 1854 quantity.semi_implicit_update = array( [1.,1.,1.,1.] ) 774 semi_implicit_update = array( [1.,1.,1.,1.] ) 775 quantity.semi_implicit_update[:] = semi_implicit_update 1855 776 1856 777 #Update with given timestep … … 1858 779 quantity.update(0.1) 1859 780 1860 sem = array([1.,1.,1.,1.])/array([1, 2, 3, 4]) 1861 denom = ones(4, Float)-timestep*sem 1862 1863 x = array([1., 2., 3., 4.]) 1864 x /= denom 1865 x += timestep*array( [4.0, 3.0, 2.0, 1.0] ) 1866 781 denom = 1.0-timestep*semi_implicit_update 782 x = (centroid_values + timestep*explicit_update)/denom 783 1867 784 assert allclose( quantity.centroid_values, x) 1868 1869 1870 1871 785 1872 786 #Test smoothing 1873 787 def test_smoothing(self): 1874 788 1875 from mesh_factory import rectangular789 1876 790 from shallow_water import Domain, Transmissive_boundary 1877 791 from Numeric import zeros, Float 1878 792 from anuga.utilities.numerical_tools import mean 1879 793 1880 #Create basic mesh1881 points, vertices, boundary = rectangular(2, 2)1882 794 1883 795 #Create shallow water domain 1884 domain = Domain(points , vertices, boundary)796 domain = Domain(points10) 1885 797 domain.default_order=2 1886 798 domain.reduction = mean … … 1888 800 1889 801 #Set some field values 1890 domain.set_quantity('elevation', lambda x ,y: x)802 domain.set_quantity('elevation', lambda x: x) 1891 803 domain.set_quantity('friction', 0.03) 1892 804 … … 2030 942 2031 943 2032 def set_array_values_by_index(self): 2033 2034 from mesh_factory import rectangular 2035 from shallow_water import Domain 2036 from Numeric import zeros, Float 2037 2038 #Create basic mesh 2039 points, vertices, boundary = rectangular(1, 1) 2040 2041 #Create shallow water domain 2042 domain = Domain(points, vertices, boundary) 2043 #print "domain.number_of_elements ",domain.number_of_elements 2044 quantity = Quantity(domain,[[1,1,1],[2,2,2]]) 2045 value = [7] 2046 indices = [1] 2047 quantity.set_array_values_by_index(value, 2048 location = 'centroids', 2049 indices = indices) 2050 #print "quantity.centroid_values",quantity.centroid_values 2051 2052 assert allclose(quantity.centroid_values, [1,7]) 2053 2054 quantity.set_array_values([15,20,25], indices = indices) 2055 assert allclose(quantity.centroid_values, [1,20]) 2056 2057 quantity.set_array_values([15,20,25], indices = indices) 2058 assert allclose(quantity.centroid_values, [1,20]) 2059 2060 def test_setting_some_vertex_values(self): 2061 """ 2062 set values based on triangle lists. 2063 """ 2064 from mesh_factory import rectangular 2065 from shallow_water import Domain 2066 from Numeric import zeros, Float 2067 2068 #Create basic mesh 2069 points, vertices, boundary = rectangular(1, 3) 2070 #print "vertices",vertices 2071 #Create shallow water domain 2072 domain = Domain(points, vertices, boundary) 2073 #print "domain.number_of_elements ",domain.number_of_elements 2074 quantity = Quantity(domain,[[1,1,1],[2,2,2],[3,3,3], 2075 [4,4,4],[5,5,5],[6,6,6]]) 2076 2077 2078 # Check that constants work 2079 value = 7 2080 indices = [1] 2081 quantity.set_values(value, 2082 location = 'centroids', 2083 indices = indices) 2084 #print "quantity.centroid_values",quantity.centroid_values 2085 assert allclose(quantity.centroid_values, [1,7,3,4,5,6]) 2086 2087 value = [7] 2088 indices = [1] 2089 quantity.set_values(value, 2090 location = 'centroids', 2091 indices = indices) 2092 #print "quantity.centroid_values",quantity.centroid_values 2093 assert allclose(quantity.centroid_values, [1,7,3,4,5,6]) 2094 2095 value = [[15,20,25]] 2096 quantity.set_values(value, indices = indices) 2097 #print "1 quantity.vertex_values",quantity.vertex_values 2098 assert allclose(quantity.vertex_values[1], value[0]) 2099 2100 2101 #print "quantity",quantity.vertex_values 2102 values = [10,100,50] 2103 quantity.set_values(values, indices = [0,1,5], location = 'centroids') 2104 #print "2 quantity.vertex_values",quantity.vertex_values 2105 assert allclose(quantity.vertex_values[0], [10,10,10]) 2106 assert allclose(quantity.vertex_values[5], [50,50,50]) 2107 #quantity.interpolate() 2108 #print "quantity.centroid_values",quantity.centroid_values 2109 assert allclose(quantity.centroid_values, [10,100,3,4,5,50]) 2110 2111 2112 quantity = Quantity(domain,[[1,1,1],[2,2,2],[3,3,3], 2113 [4,4,4],[5,5,5],[6,6,6]]) 2114 values = [10,100,50] 2115 #this will be per unique vertex, indexing the vertices 2116 #print "quantity.vertex_values",quantity.vertex_values 2117 quantity.set_values(values, indices = [0,1,5]) 2118 #print "quantity.vertex_values",quantity.vertex_values 2119 assert allclose(quantity.vertex_values[0], [1,50,10]) 2120 assert allclose(quantity.vertex_values[5], [6,6,6]) 2121 assert allclose(quantity.vertex_values[1], [100,10,50]) 2122 2123 quantity = Quantity(domain,[[1,1,1],[2,2,2],[3,3,3], 2124 [4,4,4],[5,5,5],[6,6,6]]) 2125 values = [[31,30,29],[400,400,400],[1000,999,998]] 2126 quantity.set_values(values, indices = [3,3,5]) 2127 quantity.interpolate() 2128 assert allclose(quantity.centroid_values, [1,2,3,400,5,999]) 2129 2130 values = [[1,1,1],[2,2,2],[3,3,3], 2131 [4,4,4],[5,5,5],[6,6,6]] 2132 quantity.set_values(values) 2133 2134 # testing the standard set values by vertex 2135 # indexed by vertex_id in general_mesh.coordinates 2136 values = [0,1,2,3,4,5,6,7] 2137 2138 quantity.set_values(values) 2139 #print "1 quantity.vertex_values",quantity.vertex_values 2140 assert allclose(quantity.vertex_values,[[ 4., 5., 0.], 2141 [ 1., 0., 5.], 2142 [ 5., 6., 1.], 2143 [ 2., 1., 6.], 2144 [ 6., 7., 2.], 2145 [ 3., 2., 7.]]) 2146 2147 def test_setting_unique_vertex_values(self): 2148 """ 2149 set values based on unique_vertex lists. 2150 """ 2151 from mesh_factory import rectangular 2152 from shallow_water import Domain 2153 from Numeric import zeros, Float 2154 2155 #Create basic mesh 2156 points, vertices, boundary = rectangular(1, 3) 2157 #print "vertices",vertices 2158 #Create shallow water domain 2159 domain = Domain(points, vertices, boundary) 2160 #print "domain.number_of_elements ",domain.number_of_elements 2161 quantity = Quantity(domain,[[0,0,0],[1,1,1],[2,2,2],[3,3,3], 2162 [4,4,4],[5,5,5]]) 2163 value = 7 2164 indices = [1,5] 2165 quantity.set_values(value, 2166 location = 'unique vertices', 2167 indices = indices) 2168 #print "quantity.centroid_values",quantity.centroid_values 2169 assert allclose(quantity.vertex_values[0], [0,7,0]) 2170 assert allclose(quantity.vertex_values[1], [7,1,7]) 2171 assert allclose(quantity.vertex_values[2], [7,2,7]) 2172 2173 2174 def test_get_values(self): 2175 """ 2176 get values based on triangle lists. 2177 """ 2178 from mesh_factory import rectangular 2179 from shallow_water import Domain 2180 from Numeric import zeros, Float 2181 2182 #Create basic mesh 2183 points, vertices, boundary = rectangular(1, 3) 2184 2185 #print "points",points 2186 #print "vertices",vertices 2187 #print "boundary",boundary 2188 2189 #Create shallow water domain 2190 domain = Domain(points, vertices, boundary) 2191 #print "domain.number_of_elements ",domain.number_of_elements 2192 quantity = Quantity(domain,[[0,0,0],[1,1,1],[2,2,2],[3,3,3], 2193 [4,4,4],[5,5,5]]) 2194 2195 #print "quantity.get_values(location = 'unique vertices')", \ 2196 # quantity.get_values(location = 'unique vertices') 2197 2198 #print "quantity.get_values(location = 'unique vertices')", \ 2199 # quantity.get_values(indices=[0,1,2,3,4,5,6,7], \ 2200 # location = 'unique vertices') 2201 2202 answer = [0.5,2,4,5,0,1,3,4.5] 2203 assert allclose(answer, 2204 quantity.get_values(location = 'unique vertices')) 2205 2206 indices = [0,5,3] 2207 answer = [0.5,1,5] 2208 assert allclose(answer, 2209 quantity.get_values(indices=indices, \ 2210 location = 'unique vertices')) 2211 #print "quantity.centroid_values",quantity.centroid_values 2212 #print "quantity.get_values(location = 'centroids') ",\ 2213 # quantity.get_values(location = 'centroids') 2214 2215 2216 2217 2218 def test_get_values_2(self): 2219 """Different mesh (working with domain object) - also check centroids. 2220 """ 2221 2222 2223 a = [0.0, 0.0] 2224 b = [0.0, 2.0] 2225 c = [2.0,0.0] 2226 d = [0.0, 4.0] 2227 e = [2.0, 2.0] 2228 f = [4.0,0.0] 2229 2230 points = [a, b, c, d, e, f] 2231 #bac, bce, ecf, dbe 2232 vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]] 2233 2234 domain = Domain(points, vertices) 2235 2236 quantity = Quantity(domain) 2237 quantity.set_values(lambda x, y: x+2*y) #2 4 4 6 2238 2239 assert allclose(quantity.get_values(location='centroids'), [2,4,4,6]) 2240 assert allclose(quantity.get_values(location='centroids', indices=[1,3]), [4,6]) 2241 2242 2243 assert allclose(quantity.get_values(location='vertices'), [[4,0,2], 2244 [4,2,6], 2245 [6,2,4], 2246 [8,4,6]]) 2247 2248 assert allclose(quantity.get_values(location='vertices', indices=[1,3]), [[4,2,6], 2249 [8,4,6]]) 2250 2251 2252 assert allclose(quantity.get_values(location='edges'), [[1,3,2], 2253 [4,5,3], 2254 [3,5,4], 2255 [5,7,6]]) 2256 assert allclose(quantity.get_values(location='edges', indices=[1,3]), 2257 [[4,5,3], 2258 [5,7,6]]) 2259 2260 # Check averaging over vertices 2261 #a: 0 2262 #b: (4+4+4)/3 2263 #c: (2+2+2)/3 2264 #d: 8 2265 #e: (6+6+6)/3 2266 #f: 4 2267 assert(quantity.get_values(location='unique vertices'), [0, 4, 2, 8, 6, 4]) 2268 2269 2270 2271 2272 2273 2274 def test_get_interpolated_values(self): 2275 2276 from mesh_factory import rectangular 2277 from shallow_water import Domain 2278 from Numeric import zeros, Float 2279 2280 #Create basic mesh 2281 points, vertices, boundary = rectangular(1, 3) 2282 domain = Domain(points, vertices, boundary) 2283 2284 #Constant values 2285 quantity = Quantity(domain,[[0,0,0],[1,1,1],[2,2,2],[3,3,3], 2286 [4,4,4],[5,5,5]]) 2287 2288 2289 2290 # Get interpolated values at centroids 2291 interpolation_points = domain.get_centroid_coordinates() 2292 answer = quantity.get_values(location='centroids') 2293 2294 2295 #print quantity.get_values(points=interpolation_points) 2296 assert allclose(answer, quantity.get_values(interpolation_points=interpolation_points)) 2297 2298 2299 #Arbitrary values 2300 quantity = Quantity(domain,[[0,1,2],[3,1,7],[2,1,2],[3,3,7], 2301 [1,4,-9],[2,5,0]]) 2302 2303 2304 # Get interpolated values at centroids 2305 interpolation_points = domain.get_centroid_coordinates() 2306 answer = quantity.get_values(location='centroids') 2307 #print answer 2308 #print quantity.get_values(points=interpolation_points) 2309 assert allclose(answer, quantity.get_values(interpolation_points=interpolation_points)) 2310 2311 2312 #FIXME TODO 2313 #indices = [0,5,3] 2314 #answer = [0.5,1,5] 2315 #assert allclose(answer, 2316 # quantity.get_values(indices=indices, \ 2317 # location = 'unique vertices')) 2318 2319 2320 2321 2322 def test_get_interpolated_values_2(self): 2323 a = [0.0, 0.0] 2324 b = [0.0, 2.0] 2325 c = [2.0,0.0] 2326 d = [0.0, 4.0] 2327 e = [2.0, 2.0] 2328 f = [4.0,0.0] 2329 2330 points = [a, b, c, d, e, f] 2331 #bac, bce, ecf, dbe 2332 vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]] 2333 2334 domain = Domain(points, vertices) 2335 2336 quantity = Quantity(domain) 2337 quantity.set_values(lambda x, y: x+2*y) #2 4 4 6 2338 2339 #First pick one point 2340 x, y = 2.0/3, 8.0/3 2341 v = quantity.get_values(interpolation_points = [[x,y]]) 2342 assert allclose(v, 6) 2343 2344 # Then another to test that algorithm won't blindly 2345 # reuse interpolation matrix 2346 x, y = 4.0/3, 4.0/3 2347 v = quantity.get_values(interpolation_points = [[x,y]]) 2348 assert allclose(v, 4) 2349 2350 2351 2352 2353 def test_getting_some_vertex_values(self): 2354 """ 2355 get values based on triangle lists. 2356 """ 2357 from mesh_factory import rectangular 2358 from shallow_water import Domain 2359 from Numeric import zeros, Float 2360 2361 #Create basic mesh 2362 points, vertices, boundary = rectangular(1, 3) 2363 2364 #print "points",points 2365 #print "vertices",vertices 2366 #print "boundary",boundary 2367 2368 #Create shallow water domain 2369 domain = Domain(points, vertices, boundary) 2370 #print "domain.number_of_elements ",domain.number_of_elements 2371 quantity = Quantity(domain,[[1,1,1],[2,2,2],[3,3,3], 2372 [4,4,4],[5,5,5],[6,6,6]]) 2373 value = [7] 2374 indices = [1] 2375 quantity.set_values(value, 2376 location = 'centroids', 2377 indices = indices) 2378 #print "quantity.centroid_values",quantity.centroid_values 2379 #print "quantity.get_values(location = 'centroids') ",\ 2380 # quantity.get_values(location = 'centroids') 2381 assert allclose(quantity.centroid_values, 2382 quantity.get_values(location = 'centroids')) 2383 2384 2385 value = [[15,20,25]] 2386 quantity.set_values(value, indices = indices) 2387 #print "1 quantity.vertex_values",quantity.vertex_values 2388 assert allclose(quantity.vertex_values, quantity.get_values()) 2389 2390 assert allclose(quantity.edge_values, 2391 quantity.get_values(location = 'edges')) 2392 2393 # get a subset of elements 2394 subset = quantity.get_values(location='centroids', indices=[0,5]) 2395 answer = [quantity.centroid_values[0],quantity.centroid_values[5]] 2396 assert allclose(subset, answer) 2397 2398 2399 subset = quantity.get_values(location='edges', indices=[0,5]) 2400 answer = [quantity.edge_values[0],quantity.edge_values[5]] 2401 #print "subset",subset 2402 #print "answer",answer 2403 assert allclose(subset, answer) 2404 2405 subset = quantity.get_values( indices=[1,5]) 2406 answer = [quantity.vertex_values[1],quantity.vertex_values[5]] 2407 #print "subset",subset 2408 #print "answer",answer 2409 assert allclose(subset, answer) 2410 2411 def test_smooth_vertex_values(self): 2412 """ 2413 get values based on triangle lists. 2414 """ 2415 from mesh_factory import rectangular 2416 from shallow_water import Domain 2417 from Numeric import zeros, Float 2418 2419 #Create basic mesh 2420 points, vertices, boundary = rectangular(2, 2) 2421 2422 #print "points",points 2423 #print "vertices",vertices 2424 #print "boundary",boundary 2425 2426 #Create shallow water domain 2427 domain = Domain(points, vertices, boundary) 2428 #print "domain.number_of_elements ",domain.number_of_elements 2429 quantity = Quantity(domain,[[0,0,0],[1,1,1],[2,2,2],[3,3,3], 2430 [4,4,4],[5,5,5],[6,6,6],[7,7,7]]) 2431 2432 #print "quantity.get_values(location = 'unique vertices')", \ 2433 # quantity.get_values(location = 'unique vertices') 2434 2435 #print "quantity.get_values(location = 'unique vertices')", \ 2436 # quantity.get_values(indices=[0,1,2,3,4,5,6,7], \ 2437 # location = 'unique vertices') 2438 2439 #print quantity.get_values(location = 'unique vertices') 2440 #print quantity.domain.number_of_triangles_per_node 2441 #print quantity.vertex_values 2442 2443 #answer = [0.5, 2, 3, 3, 3.5, 4, 4, 5, 6.5] 2444 #assert allclose(answer, 2445 # quantity.get_values(location = 'unique vertices')) 2446 2447 quantity.smooth_vertex_values() 2448 2449 #print quantity.vertex_values 2450 2451 2452 answer_vertex_values = [[3,3.5,0.5],[2,0.5,3.5],[3.5,4,2],[3,2,4], 2453 [4,5,3],[3.5,3,5],[5,6.5,3.5],[4,3.5,6.5]] 2454 2455 assert allclose(answer_vertex_values, 2456 quantity.vertex_values) 2457 #print "quantity.centroid_values",quantity.centroid_values 2458 #print "quantity.get_values(location = 'centroids') ",\ 2459 # quantity.get_values(location = 'centroids') 2460 2461 2462 944 2463 945 #------------------------------------------------------------- 2464 946 if __name__ == "__main__": 2465 suite = unittest.makeSuite(Test_Quantity, 'test') 2466 #suite = unittest.makeSuite(Test_Quantity, 'test_set_values_from_file_using_polygon') 2467 2468 #suite = unittest.makeSuite(Test_Quantity, 'test_set_vertex_values_using_general_interface_with_subset') 2469 #print "restricted test" 2470 #suite = unittest.makeSuite(Test_Quantity,'verbose_test_set_values_from_UTM_pts') 947 suite = unittest.makeSuite(Test_Quantity, 'test') 2471 948 runner = unittest.TextTestRunner() 2472 949 runner.run(suite)
Note: See TracChangeset
for help on using the changeset viewer.