- Timestamp:
- Apr 3, 2009, 3:03:07 PM (15 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
anuga_core/source/anuga/abstract_2d_finite_volumes/mesh_factory.py
r6145 r6717 201 201 202 202 return points, elements, boundary 203 204 205 206 def rectangular_periodic(m_g, n_g, len1_g=1.0, len2_g=1.0, origin_g = (0.0, 0.0)): 207 208 209 """Setup a rectangular grid of triangles 210 with m+1 by n+1 grid points 211 and side lengths len1, len2. If side lengths are omitted 212 the mesh defaults to the unit square, divided between all the 213 processors 214 215 len1: x direction (left to right) 216 len2: y direction (bottom to top) 217 218 """ 219 220 processor = 0 221 numproc = 1 222 223 224 n = n_g 225 m_low = -1 226 m_high = m_g +1 227 228 m = m_high - m_low 229 230 delta1 = float(len1_g)/m_g 231 delta2 = float(len2_g)/n_g 232 233 len1 = len1_g*float(m)/float(m_g) 234 len2 = len2_g 235 origin = ( origin_g[0]+float(m_low)/float(m_g)*len1_g, origin_g[1] ) 236 237 #Calculate number of points 238 Np = (m+1)*(n+1) 239 240 class VIndex: 241 242 def __init__(self, n,m): 243 self.n = n 244 self.m = m 245 246 def __call__(self, i,j): 247 return j+i*(self.n+1) 248 249 class EIndex: 250 251 def __init__(self, n,m): 252 self.n = n 253 self.m = m 254 255 def __call__(self, i,j): 256 return 2*(j+i*self.n) 257 258 259 I = VIndex(n,m) 260 E = EIndex(n,m) 261 262 points = num.zeros( (Np,2), num.Float) 263 264 for i in range(m+1): 265 for j in range(n+1): 266 267 points[I(i,j),:] = [i*delta1 + origin[0], j*delta2 + origin[1]] 268 269 #Construct 2 triangles per rectangular element and assign tags to boundary 270 #Calculate number of triangles 271 Nt = 2*m*n 272 273 274 elements = num.zeros( (Nt,3), num.Int) 275 boundary = {} 276 Idgl = [] 277 Idfl = [] 278 Idgr = [] 279 Idfr = [] 280 281 full_send_dict = {} 282 ghost_recv_dict = {} 283 nt = -1 284 for i in range(m): 285 for j in range(n): 286 287 i1 = I(i,j+1) 288 i2 = I(i,j) 289 i3 = I(i+1,j+1) 290 i4 = I(i+1,j) 291 292 #Lower Element 293 nt = E(i,j) 294 if i == 0: 295 Idgl.append(nt) 296 297 if i == 1: 298 Idfl.append(nt) 299 300 if i == m-2: 301 Idfr.append(nt) 302 303 if i == m-1: 304 Idgr.append(nt) 305 306 if i == m-1: 307 if processor == numproc-1: 308 boundary[nt, 2] = 'right' 309 else: 310 boundary[nt, 2] = 'ghost' 311 312 if j == 0: 313 boundary[nt, 1] = 'bottom' 314 elements[nt,:] = [i4,i3,i2] 315 316 #Upper Element 317 nt = E(i,j)+1 318 if i == 0: 319 Idgl.append(nt) 320 321 if i == 1: 322 Idfl.append(nt) 323 324 if i == m-2: 325 Idfr.append(nt) 326 327 if i == m-1: 328 Idgr.append(nt) 329 330 if i == 0: 331 if processor == 0: 332 boundary[nt, 2] = 'left' 333 else: 334 boundary[nt, 2] = 'ghost' 335 if j == n-1: 336 boundary[nt, 1] = 'top' 337 elements[nt,:] = [i1,i2,i3] 338 339 Idfl.extend(Idfr) 340 Idgr.extend(Idgl) 341 342 Idfl = num.array(Idfl, num.Int) 343 Idgr = num.array(Idgr, num.Int) 344 345 full_send_dict[processor] = [Idfl, Idfl] 346 ghost_recv_dict[processor] = [Idgr, Idgr] 347 348 349 return points, elements, boundary, full_send_dict, ghost_recv_dict 350 203 351 204 352 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.