 Timestamp:
 Apr 24, 2009, 5:22:14 PM (14 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

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 == m2: 300 Idfr.append(nt) 301 302 if i == m1: 303 Idgr.append(nt) 304 305 if i == m1: 306 if processor == numproc1: 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 == m2: 324 Idfr.append(nt) 325 326 if i == m1: 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 == n1: 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)):
Note: See TracChangeset
for help on using the changeset viewer.