Changeset 6717

Timestamp:

Apr 3, 2009, 3:03:07 PM (16 years ago)

Author:

steve

Message:

Added a unit test for ghost cells in sequential code

Location:

anuga_core/source/anuga

Files:

• abstract_2d_finite_volumes/domain.py (modified) (5 diffs)
• abstract_2d_finite_volumes/mesh_factory.py (modified) (1 diff)
• abstract_2d_finite_volumes/test_domain.py (modified) (2 diffs)
• abstract_2d_finite_volumes/test_general_mesh.py (modified) (1 diff)
• abstract_2d_finite_volumes/test_neighbour_mesh.py (modified) (4 diffs)
• utilities/util_ext.h (modified) (2 diffs)

anuga_core/source/anuga/abstract_2d_finite_volumes/domain.py

@@ -92 +92 @@
         """
 
+
+        number_of_full_nodes=None
+        number_of_full_triangles=None
+        
         # Determine whether source is a mesh filename or coordinates
         if type(source) == types.StringType:
@@ -117 +121 @@
 
         # Expose Mesh attributes (FIXME: Maybe turn into methods)
+        self.triangles = self.mesh.triangles        
         self.centroid_coordinates = self.mesh.centroid_coordinates
         self.vertex_coordinates = self.mesh.vertex_coordinates
@@ -205 +210 @@
         # the ghost triangles.
         if not num.allclose(self.tri_full_flag[:self.number_of_full_nodes], 1):
-            print ('WARNING:  '
-                   'Not all full triangles are store before ghost triangles')
+            if self.numproc>1:
+                print ('WARNING:  '
+                       'Not all full triangles are store before ghost triangles')
 
         # Defaults
@@ -306 +312 @@
     def get_normal(self, *args, **kwargs):
         return self.mesh.get_normal(*args, **kwargs)
+
+    def get_triangle_containing_point(self, *args, **kwargs):
+        return self.mesh.get_triangle_containing_point(*args, **kwargs)
 
     def get_intersecting_segments(self, *args, **kwargs):
@@ -1745 +1754 @@
 
     ##
-    # @brief ??
+    # @brief Sequential update of ghost cells 
     def update_ghosts(self):
-        pass
-
+        # We must send the information from the full cells and
+        # receive the information for the ghost cells
+        # We have a list with ghosts expecting updates
+
+
+        from Numeric import take,put
+
+
+        #Update of ghost cells
+        iproc = self.processor
+        if self.full_send_dict.has_key(iproc):
+
+            # now store full as local id, global id, value
+            Idf  = self.full_send_dict[iproc][0]
+
+            # now store ghost as local id, global id, value
+            Idg = self.ghost_recv_dict[iproc][0]
+
+            for i, q in enumerate(self.conserved_quantities):
+                Q_cv =  self.quantities[q].centroid_values
+                put(Q_cv,     Idg, take(Q_cv,     Idf))
+
+ 
     ##
     # @brief Extrapolate conserved quantities from centroid to vertices

anuga_core/source/anuga/abstract_2d_finite_volumes/mesh_factory.py

@@ -201 +201 @@
 
     return points, elements, boundary
+
+
+
+def rectangular_periodic(m_g, n_g, len1_g=1.0, len2_g=1.0, origin_g = (0.0, 0.0)):
+
+
+    """Setup a rectangular grid of triangles
+    with m+1 by n+1 grid points
+    and side lengths len1, len2. If side lengths are omitted
+    the mesh defaults to the unit square, divided between all the
+    processors
+
+    len1: x direction (left to right)
+    len2: y direction (bottom to top)
+
+    """
+
+    processor = 0
+    numproc   = 1
+
+    
+    n = n_g
+    m_low  = -1
+    m_high = m_g +1
+
+    m = m_high - m_low
+
+    delta1 = float(len1_g)/m_g
+    delta2 = float(len2_g)/n_g
+
+    len1 = len1_g*float(m)/float(m_g)
+    len2 = len2_g
+    origin = ( origin_g[0]+float(m_low)/float(m_g)*len1_g, origin_g[1] )
+
+    #Calculate number of points
+    Np = (m+1)*(n+1)
+
+    class VIndex:
+
+        def __init__(self, n,m):
+            self.n = n
+            self.m = m
+
+        def __call__(self, i,j):
+            return j+i*(self.n+1)
+
+    class EIndex:
+
+        def __init__(self, n,m):
+            self.n = n
+            self.m = m
+
+        def __call__(self, i,j):
+            return 2*(j+i*self.n)
+
+
+    I = VIndex(n,m)
+    E = EIndex(n,m)
+
+    points = num.zeros( (Np,2), num.Float)
+
+    for i in range(m+1):
+        for j in range(n+1):
+
+            points[I(i,j),:] = [i*delta1 + origin[0], j*delta2 + origin[1]]
+
+    #Construct 2 triangles per rectangular element and assign tags to boundary
+    #Calculate number of triangles
+    Nt = 2*m*n
+
+
+    elements = num.zeros( (Nt,3), num.Int)
+    boundary = {}
+    Idgl = []
+    Idfl = []
+    Idgr = []
+    Idfr = []
+
+    full_send_dict = {}
+    ghost_recv_dict = {}
+    nt = -1
+    for i in range(m):
+        for j in range(n):
+
+            i1 = I(i,j+1)
+            i2 = I(i,j)
+            i3 = I(i+1,j+1)
+            i4 = I(i+1,j)
+
+            #Lower Element
+            nt = E(i,j)
+            if i == 0:
+                Idgl.append(nt)
+
+            if i == 1:
+                Idfl.append(nt)
+
+            if i == m-2:
+                Idfr.append(nt)
+
+            if i == m-1:
+                Idgr.append(nt)
+
+            if i == m-1:
+                if processor == numproc-1:
+                    boundary[nt, 2] = 'right'
+                else:
+                    boundary[nt, 2] = 'ghost'
+        
+            if j == 0:
+                boundary[nt, 1] = 'bottom'
+            elements[nt,:] = [i4,i3,i2]
+
+            #Upper Element
+            nt = E(i,j)+1
+            if i == 0:
+                Idgl.append(nt)
+
+            if i == 1:
+                Idfl.append(nt)
+
+            if i == m-2:
+                Idfr.append(nt)
+
+            if i == m-1:
+                Idgr.append(nt)
+
+            if i == 0:
+                if processor == 0:
+                    boundary[nt, 2] = 'left'
+                else:
+                    boundary[nt, 2] = 'ghost'
+            if j == n-1:
+                boundary[nt, 1] = 'top'
+            elements[nt,:] = [i1,i2,i3]
+
+    Idfl.extend(Idfr)
+    Idgr.extend(Idgl)
+
+    Idfl = num.array(Idfl, num.Int)
+    Idgr = num.array(Idgr, num.Int)
+    
+    full_send_dict[processor]  = [Idfl, Idfl]
+    ghost_recv_dict[processor] = [Idgr, Idgr]
+
+
+    return  points, elements, boundary, full_send_dict, ghost_recv_dict
+
 
 def oblique(m, n, lenx = 1.0, leny = 1.0, theta = 8.95, origin = (0.0, 0.0)):

anuga_core/source/anuga/abstract_2d_finite_volumes/test_domain.py

@@ -794 +794 @@
                              [ 11.0,  11.0,  11.0]])
                              
+
+    def test_rectangular_periodic_and_ghosts(self):
+
+        from mesh_factory import rectangular_periodic
+        
+
+        M=5
+        N=2
+        points, vertices, boundary, full_send_dict, ghost_recv_dict = rectangular_periodic(M, N)
+
+        assert num.allclose(ghost_recv_dict[0][0], [24, 25, 26, 27,  0,  1,  2,  3])
+        assert num.allclose(full_send_dict[0][0] , [ 4,  5,  6,  7, 20, 21, 22, 23])
+
+        #print 'HERE'
+        
+        conserved_quantities = ['quant1', 'quant2']
+        domain = Domain(points, vertices, boundary, conserved_quantities,
+                        full_send_dict=full_send_dict,
+                        ghost_recv_dict=ghost_recv_dict)
+
+
+
+
+        assert num.allclose(ghost_recv_dict[0][0], [24, 25, 26, 27,  0,  1,  2,  3])
+        assert num.allclose(full_send_dict[0][0] , [ 4,  5,  6,  7, 20, 21, 22, 23])
+
+        def xylocation(x,y):
+            return 15*x + 9*y
+
+        
+        domain.set_quantity('quant1',xylocation,location='centroids')
+        domain.set_quantity('quant2',xylocation,location='centroids')
+
+
+        assert num.allclose(domain.quantities['quant1'].centroid_values,
+                            [  0.5,   1.,   5.,    5.5,   3.5,   4.,    8.,    8.5,   6.5,  7.,   11.,   11.5,   9.5,
+                               10.,   14.,   14.5,  12.5,  13.,   17.,   17.5,  15.5,  16.,   20.,   20.5,
+                               18.5,  19.,   23.,   23.5])
+
+
+
+        assert num.allclose(domain.quantities['quant2'].centroid_values,
+                            [  0.5,   1.,   5.,    5.5,   3.5,   4.,    8.,    8.5,   6.5,  7.,   11.,   11.5,   9.5,
+                               10.,   14.,   14.5,  12.5,  13.,   17.,   17.5,  15.5,  16.,   20.,   20.5,
+                               18.5,  19.,   23.,   23.5])
+
+        domain.update_ghosts()
+
+
+        assert num.allclose(domain.quantities['quant1'].centroid_values,
+                            [  15.5,  16.,   20.,   20.5,   3.5,   4.,    8.,    8.5,   6.5,  7.,   11.,   11.5,   9.5,
+                               10.,   14.,   14.5,  12.5,  13.,   17.,   17.5,  15.5,  16.,   20.,   20.5,
+                                3.5,   4.,    8.,    8.5])
+
+
+
+        assert num.allclose(domain.quantities['quant2'].centroid_values,
+                            [  15.5,  16.,   20.,   20.5,   3.5,   4.,    8.,    8.5,   6.5,  7.,   11.,   11.5,   9.5,
+                               10.,   14.,   14.5,  12.5,  13.,   17.,   17.5,  15.5,  16.,   20.,   20.5,
+                                3.5,   4.,    8.,    8.5])
+
+        
+        assert num.allclose(domain.tri_full_flag, [0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
+                                                   1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0])
+
+        #Test that points are arranged in a counter clock wise order
+        domain.check_integrity()
+
+
     def test_that_mesh_methods_exist(self):
         """test_that_mesh_methods_exist
@@ -819 +888 @@
         domain.get_number_of_nodes()
         domain.get_normal(0,0)
+        domain.get_triangle_containing_point([0.4,0.5])
         domain.get_intersecting_segments([[0.0, 0.0], [0.0, 1.0]])
         domain.get_disconnected_triangles()

anuga_core/source/anuga/abstract_2d_finite_volumes/test_general_mesh.py

@@ -322 +322 @@
                               nodes, triangles, geo_reference=geo)
         
-
+ 
 
 #-------------------------------------------------------------

anuga_core/source/anuga/abstract_2d_finite_volumes/test_neighbour_mesh.py

@@ -12 +12 @@
 from neighbour_mesh import *
 from mesh_factory import rectangular
+from mesh_factory import rectangular_periodic
 from anuga.config import epsilon
 
@@ -397 +398 @@
             raise "triangle edge duplicates not caught"
 
+
     def test_rectangular_mesh_basic(self):
         M=1
@@ -418 +420 @@
         assert mesh.boundary[(7,1)] == 'top' # top
         assert mesh.boundary[(3,1)] == 'top' # top
+
+
+
+
 
 
@@ -475 +481 @@
         points, vertices, boundary = rectangular(2*N, N, len1=10, len2=10)
         mesh = Mesh(points, vertices, boundary)
-
 
 

anuga_core/source/anuga/utilities/util_ext.h

@@ -261 +261 @@
     
   B = (PyArrayObject*) PyObject_GetAttrString(O, name);
+
+  //printf("B = %p\n",(void*)B);
   if (!B) {
-    PyErr_SetString(PyExc_RuntimeError, "util_ext.h: get_consecutive_array could not obtain array");
+    printf("util_ext.h: get_consecutive_array could not obtain python object");
+    printf(" %s\n",name);
+    fflush(stdout);
+    PyErr_SetString(PyExc_RuntimeError, "util_ext.h: get_consecutive_array could not obtain python object");
     return NULL;
   }     
@@ -273 +278 @@
   
   if (!A) {
+    printf("util_ext.h: get_consecutive_array could not obtain array object");
+    printf(" %s \n",name);
+    fflush(stdout);
     PyErr_SetString(PyExc_RuntimeError, "util_ext.h: get_consecutive_array could not obtain array");
     return NULL;
   }
+
+
   return A;
 }