Changeset 7484 for anuga_core/source/anuga

Timestamp:

Sep 4, 2009, 6:22:15 PM (16 years ago)

Author:

steve

Message:

Added edge_midpoint_coordinates together with get_edge_midpoint_coordinates
and get_edge_midpoint_coordinate into the general mesh. This is so that
boundary classes can get access to the location of the edge midpoint.

I will use this to allow Dirichlet_boundary to take a function of
time and space.

Location:

anuga_core/source/anuga/abstract_2d_finite_volumes

Files:

• general_mesh.py (modified) (3 diffs)
• test_general_mesh.py (modified) (2 diffs)

anuga_core/source/anuga/abstract_2d_finite_volumes/general_mesh.py

@@ -144 +144 @@
         self.edgelengths = num.zeros((N, 3), num.float)
 
-        # Get x,y coordinates for all triangles and store
+        # Get x,y coordinates for all triangle vertices and store
         self.vertex_coordinates = V = self.compute_vertex_coordinates()
+
+        # Get x,y coordinates for all triangle edge midpoints and store
+        self.edge_midpoint_coordinates  = self.compute_edge_midpoint_coordinates()
 
         # Initialise each triangle
@@ -327 +330 @@
         return V[j,:]
 
+
+
     def compute_vertex_coordinates(self):
         """Return all vertex coordinates for all triangles as a 3*M x 2 array
@@ -344 +349 @@
 
         return vertex_coordinates
+
+
+    def get_edge_midpoint_coordinates(self, triangle_id=None, absolute=False):
+        """Return edge midpoint coordinates for all triangles or from particular triangle.
+
+        Return all edge midpoint coordinates for all triangles as a 3*M x 2 array
+        where the jth midpoint of the ith triangle is located in row 3*i+j and
+        M the number of triangles in the mesh.
+
+        if triangle_id is specified (an integer) the 3 midpoint coordinates
+        for triangle_id are returned.
+
+        Boolean keyword argument absolute determines whether coordinates
+        are to be made absolute by taking georeference into account
+        Default is False as many parts of ANUGA expects relative coordinates.
+        """
+
+        E = self.edge_midpoint_coordinates
+        
+        if triangle_id is None:
+            if absolute is True:
+                if not self.geo_reference.is_absolute():
+                    E = self.geo_reference.get_absolute(E)
+            return E
+        else:
+            i = triangle_id
+            msg = 'triangle_id must be an integer'
+            assert int(i) == i, msg
+            assert 0 <= i < self.number_of_triangles
+
+            i3 = 3*i
+            if absolute is True and not self.geo_reference.is_absolute():
+                offset=num.array([self.geo_reference.get_xllcorner(),
+                                  self.geo_reference.get_yllcorner()], num.float)
+                return num.array([E[i3,:]+offset,
+                                  E[i3+1,:]+offset,
+                                  E[i3+2,:]+offset], num.float)
+            else:
+                return num.array([E[i3,:], E[i3+1,:], E[i3+2,:]], num.float)    
+
+
+    def get_edge_midpoint_coordinate(self, i, j, absolute=False):
+        """Return coordinates for edge midpoint j of the i'th triangle.
+        Return value is the numeric array slice [x, y]
+        """
+
+        msg = 'edge midpoint id j must be an integer in [0,1,2]'
+        assert j in [0,1,2], msg
+
+        E = self.get_edge_midpoint_coordinates(triangle_id=i, absolute=absolute)
+        return E[j,:]
+
+    
+    def compute_edge_midpoint_coordinates(self):
+        """Return all edge midpoint coordinates for all triangles as a 3*M x 2 array
+        where the jth edge midpoint of the ith triangle is located in row 3*i+j.
+
+        This function is used to precompute this important structure. Use
+        get_edge_midpoint_coordinates to retrieve the points.
+
+        Assumes that vertex_coordinates have been computed
+        """
+
+        M = self.number_of_triangles
+        E = num.zeros((3*M, 2), num.float)
+
+        V = self.vertex_coordinates
+
+        V0 = V[0:3*M:3, :]
+        V1 = V[1:3*M:3, :]
+        V2 = V[2:3*M:3, :]
+
+        
+        #print V.shape, V0.shape, V1.shape, V2.shape
+
+        #print E.shape, E[0:3*M:3, :].shape, E[1:3*M:3, :].shape, E[2:3*M:3, :].shape
+        E[0:3*M:3, :] = 0.5*(V1+V2)
+        E[1:3*M:3, :] = 0.5*(V2+V0)
+        E[2:3*M:3, :] = 0.5*(V0+V1)
+
+        return E
+
+
 
     def get_triangles(self, indices=None):

anuga_core/source/anuga/abstract_2d_finite_volumes/test_general_mesh.py

@@ -118 +118 @@
                 k = triangles[i,j]  #Index of vertex j in triangle i
                 assert num.allclose(V[j,:], nodes[k])
+
+
+    def test_get_edge_midpoint_coordinates(self):
+        from mesh_factory import rectangular
+
+        #Create basic mesh
+        nodes, triangles, _ = rectangular(1, 3)
+        domain = General_mesh(nodes, triangles)
+
+
+        assert num.allclose(domain.get_nodes(), nodes)
+
+
+        M = domain.number_of_triangles        
+
+        E = domain.get_edge_midpoint_coordinates()
+        assert E.shape[0] == 3*M
+
+        for i in range(M):
+            k0 = triangles[i,0]  #Index of vertex 0 in triangle i
+            k1 = triangles[i,1]  #Index of vertex 1 in triangle i
+            k2 = triangles[i,2]  #Index of vertex 2 in triangle i
+
+            assert num.allclose(E[3*i+0,:], 0.5*(nodes[k1]+nodes[k2]) )
+            assert num.allclose(E[3*i+1,:], 0.5*(nodes[k0]+nodes[k2]) )
+            assert num.allclose(E[3*i+2,:], 0.5*(nodes[k1]+nodes[k0]) )
+
+    def test_get_edge_midpoint_coordinates_with_geo_ref(self):
+        x0 = 314036.58727982
+        y0 = 6224951.2960092
+        geo = Geo_reference(56, x0, y0)
+        
+        a = num.array([0.0, 0.0])
+        b = num.array([0.0, 2.0])
+        c = num.array([2.0, 0.0])
+        d = num.array([0.0, 4.0])
+        e = num.array([2.0, 2.0])
+        f = num.array([4.0, 0.0])
+        nodes = num.array([a, b, c, d, e, f])
+
+        nodes_absolute = geo.get_absolute(nodes)
+        
+        #                        bac,     bce,     ecf,     dbe
+        triangles = num.array([[1,0,2], [1,2,4], [4,2,5], [3,1,4]], num.int)
+
+        domain = General_mesh(nodes, triangles, geo_reference=geo)
+
+        verts = domain.get_edge_midpoint_coordinates(triangle_id=0)    # bac
+        msg = ("num.array(1/2[a+c,b+c,a+b])=\n%s\nshould be close to 'verts'=\n%s"
+               % (str(num.array([0.5*(a+c),0.5*(b+c),0.5*(a+b)])), str(verts)))
+        self.failUnless(num.allclose(num.array([0.5*(a+c),0.5*(b+c),0.5*(a+b)]), verts), msg)
+
+
+        verts = domain.get_edge_midpoint_coordinates(triangle_id=0, absolute=True)
+        msg = ("num.array([...])=\n%s\nshould be close to 'verts'=\n%s"
+               % (str(0.5*num.array([nodes_absolute[0]+nodes_absolute[2],
+                                     nodes_absolute[1]+nodes_absolute[2],
+                                     nodes_absolute[1]+nodes_absolute[0]])),
+                  str(verts)))
+        self.assert_(num.allclose(0.5*num.array([nodes_absolute[0]+nodes_absolute[2],
+                                                 nodes_absolute[1]+nodes_absolute[2],
+                                                 nodes_absolute[1]+nodes_absolute[0]]),
+                                  verts), msg)
+
+    def test_get_edge_midpoint_coordinates_triangle_id(self):
+        """test_get_vertex_coordinates_triangle_id
+        Test that vertices for one triangle can be returned.
+        """
+        from mesh_factory import rectangular
+
+        #Create basic mesh
+        nodes, triangles, _ = rectangular(1, 3)
+        domain = General_mesh(nodes, triangles)
+
+
+        assert num.allclose(domain.get_nodes(), nodes)
+
+
+        M = domain.number_of_triangles        
+
+        for i in range(M):
+            E = domain.get_edge_midpoint_coordinates(triangle_id=i)
+            assert E.shape[0] == 3
+
+
+            k0 = triangles[i,0]  #Index of vertex 0 in triangle i
+            k1 = triangles[i,1]  #Index of vertex 0 in triangle i
+            k2 = triangles[i,2]  #Index of vertex 0 in triangle i
+
+            assert num.allclose(E[0,:], 0.5*(nodes[k1]+nodes[k2]))
+            assert num.allclose(E[1,:], 0.5*(nodes[k0]+nodes[k2]))
+            assert num.allclose(E[2,:], 0.5*(nodes[k1]+nodes[k0]))
+
+            E0 = domain.get_edge_midpoint_coordinate(i, 0 )
+            E1 = domain.get_edge_midpoint_coordinate(i, 1 )
+            E2 = domain.get_edge_midpoint_coordinate(i, 2 )
+
+            assert num.allclose(E0, 0.5*(nodes[k1]+nodes[k2]))
+            assert num.allclose(E1, 0.5*(nodes[k0]+nodes[k2]))
+            assert num.allclose(E2, 0.5*(nodes[k1]+nodes[k0]))
+            
 
 
@@ -356 +457 @@
 
 if __name__ == "__main__":
-    suite = unittest.makeSuite(Test_General_Mesh, 'test') 
+    #suite = unittest.makeSuite(Test_General_Mesh, 'test')
+    suite = unittest.makeSuite(Test_General_Mesh, 'test_get_edge_midpoint_coordinates')     
     runner = unittest.TextTestRunner()
     runner.run(suite)