Changeset 8553

Timestamp:

Sep 4, 2012, 4:26:13 PM (13 years ago)

Author:

steve

Message:

Committing faster ghost_layer code

Location:

trunk/anuga_core/source/anuga_parallel

Files:

• distribute_mesh.py (modified) (4 diffs)
• parallel_api.py (modified) (1 diff)
• run_parallel_sw_rectangular_cross.py (modified) (2 diffs)

trunk/anuga_core/source/anuga_parallel/distribute_mesh.py

@@ -24 +24 @@
 
 import numpy as num
+import numpy.lib.arraysetops as numset
 
 from anuga.abstract_2d_finite_volumes.neighbour_mesh import Mesh
@@ -498 +499 @@
 #########################################################
 
-def ghost_layer(submesh, mesh, p, tupper, tlower, parameters = None):
+def ghost_layer_old(submesh, mesh, p, tupper, tlower, parameters = None):
 
     ncoord = mesh.number_of_nodes
@@ -590 +591 @@
 
     return subnodes, subtriangles, layer_width
+
+
+
+def ghost_layer(submesh, mesh, p, tupper, tlower, parameters = None):
+
+    ncoord = mesh.number_of_nodes
+    ntriangles = mesh.number_of_triangles
+
+    if parameters is None:
+        layer_width  = 2
+    else:
+        layer_width = parameters['ghost_layer_width']
+
+
+    full_ids = num.arange(tlower, tupper)
+
+    n0 = mesh.neighbours[full_ids, :]
+    n0 = num.unique(n0.flat)
+    n0 = num.extract(n0>=0,n0)
+    n0 = num.extract(num.logical_or(n0<tlower, tupper<= n0), n0)
+
+    layer_cells = {}
+    layer_cells[0] = n0
+
+
+    # Find the subsequent layers of ghost triangles
+    for i in range(layer_width-1):
+
+        # use previous layer as a start
+        n0 = mesh.neighbours[n0, :]
+        n0 = num.unique(n0.flat)
+        n0 = num.extract(n0>=0,n0)
+        n0 = num.extract(num.logical_or(n0<tlower, tupper<= n0), n0)
+
+        for j in xrange(i+1):
+            n0 = numset.setdiff1d(n0,layer_cells[j])
+
+        layer_cells[i+1] = n0
+
+
+    # Build the triangle list and make note of the vertices
+    new_trianglemap = layer_cells[0]
+    for i in range(layer_width-1):
+        new_trianglemap = numset.union1d(new_trianglemap,layer_cells[i+1])
+
+    new_subtriangles = num.concatenate((num.reshape(new_trianglemap, (-1,1)), mesh.triangles[new_trianglemap]), 1)
+
+
+
+
+    fullnodes = submesh["full_nodes"][p]
+    full_nodes_ids = num.array(fullnodes[:,0],num.int)
+
+    new_nodes = num.unique(mesh.triangles[new_trianglemap].flat)
+    new_nodes = numset.setdiff1d(new_nodes,full_nodes_ids)
+
+    new_subnodes = num.concatenate((num.reshape(new_nodes, (-1,1)), mesh.nodes[new_nodes]), 1)
+
+    # Clean up before exiting
+
+    del (new_nodes)
+    del (layer_cells)
+    del (n0)
+    del (new_trianglemap)
+
+    # Return the triangles and vertices sitting on the boundary layer
+
+    return new_subnodes, new_subtriangles, layer_width
 
 #########################################################
@@ -925 +994 @@
 #########################################################
 
-def build_submesh(mesh, quantities,
+def build_submesh(nodes, triangles, boundary, quantities,
                   triangles_per_proc, parameters = None):
 
     # Temporarily build the mesh to find the neighbouring
-    # triangles and true boundary polygon
-
-    #mesh = Mesh(nodes, triangles, boundary)
+    # triangles and true boundary polygon\
+
+    mesh = Mesh(nodes, triangles, boundary)
     boundary_polygon = mesh.get_boundary_polygon()
     

trunk/anuga_core/source/anuga_parallel/parallel_api.py

@@ -266 +266 @@
     #      sp2_map (maps parallel domain triangles to domain triangles)
 
-    mesh = Mesh(new_nodes, new_triangles, new_boundary)
+
 
     # Build the mesh that should be assigned to each processor,
     # this includes ghost nodes and the communication pattern
     if verbose: print 'Build submeshes'    
-    submesh = build_submesh(mesh, quantities, triangles_per_proc, parameters)
+    submesh = build_submesh(new_nodes, new_triangles, new_boundary, quantities, triangles_per_proc, parameters)
 
     if verbose:

trunk/anuga_core/source/anuga_parallel/run_parallel_sw_rectangular_cross.py

@@ -42 +42 @@
     width = 2.0
     dx = dy = 0.005
-    dx = dy = 0.0025
+    dx = dy = 0.005
+    #dx = dy  = 0.5
     domain = rectangular_cross_domain(int(length/dx), int(width/dy),
                                               len1=length, len2=width)
@@ -64 +65 @@
 
 # setup parameters to test using different ghost_layer_widths
-parameters = dict(ghost_layer_width = 2)
+parameters = dict(ghost_layer_width = 4)
 domain = distribute(domain,verbose=verbose, parameters=parameters)