Changeset 6902 for branches/numpy/anuga/abstract_2d_finite_volumes
- Timestamp:
- Apr 24, 2009, 5:22:14 PM (16 years ago)
- Location:
- branches/numpy/anuga/abstract_2d_finite_volumes
- Files:
-
- 5 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/numpy/anuga/abstract_2d_finite_volumes/domain.py
r6553 r6902 92 92 """ 93 93 94 number_of_full_nodes=None 95 number_of_full_triangles=None 96 94 97 # Determine whether source is a mesh filename or coordinates 95 98 if type(source) == types.StringType: … … 117 120 118 121 # Expose Mesh attributes (FIXME: Maybe turn into methods) 122 self.triangles = self.mesh.triangles 119 123 self.centroid_coordinates = self.mesh.centroid_coordinates 120 124 self.vertex_coordinates = self.mesh.vertex_coordinates … … 205 209 # the ghost triangles. 206 210 if not num.allclose(self.tri_full_flag[:self.number_of_full_nodes], 1): 207 print ('WARNING: ' 208 'Not all full triangles are store before ghost triangles') 211 if self.numproc>1: 212 print ('WARNING: ' 213 'Not all full triangles are store before ghost triangles') 209 214 210 215 # Defaults … … 306 311 def get_normal(self, *args, **kwargs): 307 312 return self.mesh.get_normal(*args, **kwargs) 313 def get_triangle_containing_point(self, *args, **kwargs): 314 return self.mesh.get_triangle_containing_point(*args, **kwargs) 308 315 309 316 def get_intersecting_segments(self, *args, **kwargs): … … 1745 1752 1746 1753 ## 1747 # @brief ??1754 # @brief Sequential update of ghost cells 1748 1755 def update_ghosts(self): 1749 pass 1750 1756 # We must send the information from the full cells and 1757 # receive the information for the ghost cells 1758 # We have a list with ghosts expecting updates 1759 1760 from Numeric import take,put 1761 1762 1763 #Update of ghost cells 1764 iproc = self.processor 1765 if self.full_send_dict.has_key(iproc): 1766 1767 # now store full as local id, global id, value 1768 Idf = self.full_send_dict[iproc][0] 1769 1770 # now store ghost as local id, global id, value 1771 Idg = self.ghost_recv_dict[iproc][0] 1772 1773 for i, q in enumerate(self.conserved_quantities): 1774 Q_cv = self.quantities[q].centroid_values 1775 put(Q_cv, Idg, take(Q_cv, Idf)) 1776 1777 1751 1778 ## 1752 1779 # @brief Extrapolate conserved quantities from centroid to vertices -
branches/numpy/anuga/abstract_2d_finite_volumes/generic_boundary_conditions.py
r6553 r6902 82 82 """Time dependent boundary returns values for the 83 83 conserved quantities as a function of time. 84 Must specify domain to get access to model time and a function f(t) 85 which must return conserved quantities as a function time 84 Must specify domain to get access to model time and a function of t 85 which must return conserved quantities as a function time. 86 87 Example: 88 B = Time_boundary(domain, 89 function=lambda t: [(60<t<3660)*2, 0, 0]) 90 91 This will produce a boundary condition with is a 2m high square wave 92 starting 60 seconds into the simulation and lasting one hour. 93 Momentum applied will be 0 at all times. 94 86 95 """ 87 96 … … 98 107 self.verbose = verbose 99 108 109 110 # FIXME: Temporary code to deal with both f and function 111 if function is not None and f is not None: 112 raise Exception, 'Specify either function or f to Time_boundary' 113 114 if function is None and f is None: 115 raise Exception, 'You must specify a function to Time_boundary' 116 117 if f is None: 118 f = function 119 ##### 120 100 121 try: 101 122 q = f(0.0) -
branches/numpy/anuga/abstract_2d_finite_volumes/mesh_factory.py
r6304 r6902 201 201 202 202 return points, elements, boundary 203 204 205 def rectangular_periodic(m_g, n_g, len1_g=1.0, len2_g=1.0, origin_g = (0.0, 0.0)): 206 207 208 """Setup a rectangular grid of triangles 209 with m+1 by n+1 grid points 210 and side lengths len1, len2. If side lengths are omitted 211 the mesh defaults to the unit square, divided between all the 212 processors 213 214 len1: x direction (left to right) 215 len2: y direction (bottom to top) 216 217 """ 218 219 processor = 0 220 numproc = 1 221 222 223 n = n_g 224 m_low = -1 225 m_high = m_g +1 226 227 m = m_high - m_low 228 229 delta1 = float(len1_g)/m_g 230 delta2 = float(len2_g)/n_g 231 232 len1 = len1_g*float(m)/float(m_g) 233 len2 = len2_g 234 origin = ( origin_g[0]+float(m_low)/float(m_g)*len1_g, origin_g[1] ) 235 236 #Calculate number of points 237 Np = (m+1)*(n+1) 238 239 class VIndex: 240 241 def __init__(self, n,m): 242 self.n = n 243 self.m = m 244 245 def __call__(self, i,j): 246 return j+i*(self.n+1) 247 248 class EIndex: 249 250 def __init__(self, n,m): 251 self.n = n 252 self.m = m 253 254 def __call__(self, i,j): 255 return 2*(j+i*self.n) 256 257 258 I = VIndex(n,m) 259 E = EIndex(n,m) 260 261 points = num.zeros( (Np,2), num.Float) 262 263 for i in range(m+1): 264 for j in range(n+1): 265 266 points[I(i,j),:] = [i*delta1 + origin[0], j*delta2 + origin[1]] 267 268 #Construct 2 triangles per rectangular element and assign tags to boundary 269 #Calculate number of triangles 270 Nt = 2*m*n 271 272 273 elements = num.zeros( (Nt,3), num.Int) 274 boundary = {} 275 Idgl = [] 276 Idfl = [] 277 Idgr = [] 278 Idfr = [] 279 280 full_send_dict = {} 281 ghost_recv_dict = {} 282 nt = -1 283 for i in range(m): 284 for j in range(n): 285 286 i1 = I(i,j+1) 287 i2 = I(i,j) 288 i3 = I(i+1,j+1) 289 i4 = I(i+1,j) 290 291 #Lower Element 292 nt = E(i,j) 293 if i == 0: 294 Idgl.append(nt) 295 296 if i == 1: 297 Idfl.append(nt) 298 299 if i == m-2: 300 Idfr.append(nt) 301 302 if i == m-1: 303 Idgr.append(nt) 304 305 if i == m-1: 306 if processor == numproc-1: 307 boundary[nt, 2] = 'right' 308 else: 309 boundary[nt, 2] = 'ghost' 310 311 if j == 0: 312 boundary[nt, 1] = 'bottom' 313 elements[nt,:] = [i4,i3,i2] 314 315 #Upper Element 316 nt = E(i,j)+1 317 if i == 0: 318 Idgl.append(nt) 319 320 if i == 1: 321 Idfl.append(nt) 322 323 if i == m-2: 324 Idfr.append(nt) 325 326 if i == m-1: 327 Idgr.append(nt) 328 329 if i == 0: 330 if processor == 0: 331 boundary[nt, 2] = 'left' 332 else: 333 boundary[nt, 2] = 'ghost' 334 if j == n-1: 335 boundary[nt, 1] = 'top' 336 elements[nt,:] = [i1,i2,i3] 337 338 Idfl.extend(Idfr) 339 Idgr.extend(Idgl) 340 341 Idfl = num.array(Idfl, num.Int) 342 Idgr = num.array(Idgr, num.Int) 343 344 full_send_dict[processor] = [Idfl, Idfl] 345 ghost_recv_dict[processor] = [Idgr, Idgr] 346 347 348 return points, elements, boundary, full_send_dict, ghost_recv_dict 349 203 350 204 351 def oblique(m, n, lenx = 1.0, leny = 1.0, theta = 8.95, origin = (0.0, 0.0)): -
branches/numpy/anuga/abstract_2d_finite_volumes/quantity.py
r6689 r6902 1178 1178 # Edges have already been deprecated in set_values, see changeset:5521, 1179 1179 # but *might* be useful in get_values. Any thoughts anyone? 1180 # YES (Ole): Edge values are necessary for volumetric balance check 1180 1181 1181 1182 if location not in ['vertices', 'centroids', -
branches/numpy/anuga/abstract_2d_finite_volumes/test_domain.py
r6553 r6902 397 397 398 398 def test_setting_timestepping_method(self): 399 """test_set _quanitities_to_be_monitored399 """test_setting_timestepping_method 400 400 """ 401 401 … … 803 803 [ 11.0, 11.0, 11.0]]) 804 804 805 def test_rectangular_periodic_and_ghosts(self): 806 807 from mesh_factory import rectangular_periodic 808 809 810 M=5 811 N=2 812 points, vertices, boundary, full_send_dict, ghost_recv_dict = rectangular_periodic(M, N) 813 814 assert num.allclose(ghost_recv_dict[0][0], [24, 25, 26, 27, 0, 1, 2, 3]) 815 assert num.allclose(full_send_dict[0][0] , [ 4, 5, 6, 7, 20, 21, 22, 23]) 816 817 #print 'HERE' 818 819 conserved_quantities = ['quant1', 'quant2'] 820 domain = Domain(points, vertices, boundary, conserved_quantities, 821 full_send_dict=full_send_dict, 822 ghost_recv_dict=ghost_recv_dict) 823 824 825 826 827 assert num.allclose(ghost_recv_dict[0][0], [24, 25, 26, 27, 0, 1, 2, 3]) 828 assert num.allclose(full_send_dict[0][0] , [ 4, 5, 6, 7, 20, 21, 22, 23]) 829 830 def xylocation(x,y): 831 return 15*x + 9*y 832 833 834 domain.set_quantity('quant1',xylocation,location='centroids') 835 domain.set_quantity('quant2',xylocation,location='centroids') 836 837 838 assert num.allclose(domain.quantities['quant1'].centroid_values, 839 [ 0.5, 1., 5., 5.5, 3.5, 4., 8., 8.5, 6.5, 7., 11., 11.5, 9.5, 840 10., 14., 14.5, 12.5, 13., 17., 17.5, 15.5, 16., 20., 20.5, 841 18.5, 19., 23., 23.5]) 842 843 844 845 assert num.allclose(domain.quantities['quant2'].centroid_values, 846 [ 0.5, 1., 5., 5.5, 3.5, 4., 8., 8.5, 6.5, 7., 11., 11.5, 9.5, 847 10., 14., 14.5, 12.5, 13., 17., 17.5, 15.5, 16., 20., 20.5, 848 18.5, 19., 23., 23.5]) 849 850 domain.update_ghosts() 851 852 853 assert num.allclose(domain.quantities['quant1'].centroid_values, 854 [ 15.5, 16., 20., 20.5, 3.5, 4., 8., 8.5, 6.5, 7., 11., 11.5, 9.5, 855 10., 14., 14.5, 12.5, 13., 17., 17.5, 15.5, 16., 20., 20.5, 856 3.5, 4., 8., 8.5]) 857 858 859 860 assert num.allclose(domain.quantities['quant2'].centroid_values, 861 [ 15.5, 16., 20., 20.5, 3.5, 4., 8., 8.5, 6.5, 7., 11., 11.5, 9.5, 862 10., 14., 14.5, 12.5, 13., 17., 17.5, 15.5, 16., 20., 20.5, 863 3.5, 4., 8., 8.5]) 864 865 866 assert num.allclose(domain.tri_full_flag, [0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 867 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0]) 868 869 #Test that points are arranged in a counter clock wise order 870 domain.check_integrity() 871 872 805 873 def test_that_mesh_methods_exist(self): 806 874 """test_that_mesh_methods_exist … … 828 896 domain.get_number_of_nodes() 829 897 domain.get_normal(0,0) 898 domain.get_triangle_containing_point([0.4,0.5]) 830 899 domain.get_intersecting_segments([[0.0, 0.0], [0.0, 1.0]]) 831 900 domain.get_disconnected_triangles()
Note: See TracChangeset
for help on using the changeset viewer.