Changeset 7484


Ignore:
Timestamp:
Sep 4, 2009, 6:22:15 PM (11 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:
2 edited

Legend:

Unmodified
Added
Removed
  • anuga_core/source/anuga/abstract_2d_finite_volumes/general_mesh.py

    r7317 r7484  
    144144        self.edgelengths = num.zeros((N, 3), num.float)
    145145
    146         # Get x,y coordinates for all triangles and store
     146        # Get x,y coordinates for all triangle vertices and store
    147147        self.vertex_coordinates = V = self.compute_vertex_coordinates()
     148
     149        # Get x,y coordinates for all triangle edge midpoints and store
     150        self.edge_midpoint_coordinates  = self.compute_edge_midpoint_coordinates()
    148151
    149152        # Initialise each triangle
     
    327330        return V[j,:]
    328331
     332
     333
    329334    def compute_vertex_coordinates(self):
    330335        """Return all vertex coordinates for all triangles as a 3*M x 2 array
     
    344349
    345350        return vertex_coordinates
     351
     352
     353    def get_edge_midpoint_coordinates(self, triangle_id=None, absolute=False):
     354        """Return edge midpoint coordinates for all triangles or from particular triangle.
     355
     356        Return all edge midpoint coordinates for all triangles as a 3*M x 2 array
     357        where the jth midpoint of the ith triangle is located in row 3*i+j and
     358        M the number of triangles in the mesh.
     359
     360        if triangle_id is specified (an integer) the 3 midpoint coordinates
     361        for triangle_id are returned.
     362
     363        Boolean keyword argument absolute determines whether coordinates
     364        are to be made absolute by taking georeference into account
     365        Default is False as many parts of ANUGA expects relative coordinates.
     366        """
     367
     368        E = self.edge_midpoint_coordinates
     369       
     370        if triangle_id is None:
     371            if absolute is True:
     372                if not self.geo_reference.is_absolute():
     373                    E = self.geo_reference.get_absolute(E)
     374            return E
     375        else:
     376            i = triangle_id
     377            msg = 'triangle_id must be an integer'
     378            assert int(i) == i, msg
     379            assert 0 <= i < self.number_of_triangles
     380
     381            i3 = 3*i
     382            if absolute is True and not self.geo_reference.is_absolute():
     383                offset=num.array([self.geo_reference.get_xllcorner(),
     384                                  self.geo_reference.get_yllcorner()], num.float)
     385                return num.array([E[i3,:]+offset,
     386                                  E[i3+1,:]+offset,
     387                                  E[i3+2,:]+offset], num.float)
     388            else:
     389                return num.array([E[i3,:], E[i3+1,:], E[i3+2,:]], num.float)   
     390
     391
     392    def get_edge_midpoint_coordinate(self, i, j, absolute=False):
     393        """Return coordinates for edge midpoint j of the i'th triangle.
     394        Return value is the numeric array slice [x, y]
     395        """
     396
     397        msg = 'edge midpoint id j must be an integer in [0,1,2]'
     398        assert j in [0,1,2], msg
     399
     400        E = self.get_edge_midpoint_coordinates(triangle_id=i, absolute=absolute)
     401        return E[j,:]
     402
     403   
     404    def compute_edge_midpoint_coordinates(self):
     405        """Return all edge midpoint coordinates for all triangles as a 3*M x 2 array
     406        where the jth edge midpoint of the ith triangle is located in row 3*i+j.
     407
     408        This function is used to precompute this important structure. Use
     409        get_edge_midpoint_coordinates to retrieve the points.
     410
     411        Assumes that vertex_coordinates have been computed
     412        """
     413
     414        M = self.number_of_triangles
     415        E = num.zeros((3*M, 2), num.float)
     416
     417        V = self.vertex_coordinates
     418
     419        V0 = V[0:3*M:3, :]
     420        V1 = V[1:3*M:3, :]
     421        V2 = V[2:3*M:3, :]
     422
     423       
     424        #print V.shape, V0.shape, V1.shape, V2.shape
     425
     426        #print E.shape, E[0:3*M:3, :].shape, E[1:3*M:3, :].shape, E[2:3*M:3, :].shape
     427        E[0:3*M:3, :] = 0.5*(V1+V2)
     428        E[1:3*M:3, :] = 0.5*(V2+V0)
     429        E[2:3*M:3, :] = 0.5*(V0+V1)
     430
     431        return E
     432
     433
    346434
    347435    def get_triangles(self, indices=None):
  • anuga_core/source/anuga/abstract_2d_finite_volumes/test_general_mesh.py

    r7276 r7484  
    118118                k = triangles[i,j]  #Index of vertex j in triangle i
    119119                assert num.allclose(V[j,:], nodes[k])
     120
     121
     122    def test_get_edge_midpoint_coordinates(self):
     123        from mesh_factory import rectangular
     124
     125        #Create basic mesh
     126        nodes, triangles, _ = rectangular(1, 3)
     127        domain = General_mesh(nodes, triangles)
     128
     129
     130        assert num.allclose(domain.get_nodes(), nodes)
     131
     132
     133        M = domain.number_of_triangles       
     134
     135        E = domain.get_edge_midpoint_coordinates()
     136        assert E.shape[0] == 3*M
     137
     138        for i in range(M):
     139            k0 = triangles[i,0]  #Index of vertex 0 in triangle i
     140            k1 = triangles[i,1]  #Index of vertex 1 in triangle i
     141            k2 = triangles[i,2]  #Index of vertex 2 in triangle i
     142
     143            assert num.allclose(E[3*i+0,:], 0.5*(nodes[k1]+nodes[k2]) )
     144            assert num.allclose(E[3*i+1,:], 0.5*(nodes[k0]+nodes[k2]) )
     145            assert num.allclose(E[3*i+2,:], 0.5*(nodes[k1]+nodes[k0]) )
     146
     147    def test_get_edge_midpoint_coordinates_with_geo_ref(self):
     148        x0 = 314036.58727982
     149        y0 = 6224951.2960092
     150        geo = Geo_reference(56, x0, y0)
     151       
     152        a = num.array([0.0, 0.0])
     153        b = num.array([0.0, 2.0])
     154        c = num.array([2.0, 0.0])
     155        d = num.array([0.0, 4.0])
     156        e = num.array([2.0, 2.0])
     157        f = num.array([4.0, 0.0])
     158        nodes = num.array([a, b, c, d, e, f])
     159
     160        nodes_absolute = geo.get_absolute(nodes)
     161       
     162        #                        bac,     bce,     ecf,     dbe
     163        triangles = num.array([[1,0,2], [1,2,4], [4,2,5], [3,1,4]], num.int)
     164
     165        domain = General_mesh(nodes, triangles, geo_reference=geo)
     166
     167        verts = domain.get_edge_midpoint_coordinates(triangle_id=0)    # bac
     168        msg = ("num.array(1/2[a+c,b+c,a+b])=\n%s\nshould be close to 'verts'=\n%s"
     169               % (str(num.array([0.5*(a+c),0.5*(b+c),0.5*(a+b)])), str(verts)))
     170        self.failUnless(num.allclose(num.array([0.5*(a+c),0.5*(b+c),0.5*(a+b)]), verts), msg)
     171
     172
     173        verts = domain.get_edge_midpoint_coordinates(triangle_id=0, absolute=True)
     174        msg = ("num.array([...])=\n%s\nshould be close to 'verts'=\n%s"
     175               % (str(0.5*num.array([nodes_absolute[0]+nodes_absolute[2],
     176                                     nodes_absolute[1]+nodes_absolute[2],
     177                                     nodes_absolute[1]+nodes_absolute[0]])),
     178                  str(verts)))
     179        self.assert_(num.allclose(0.5*num.array([nodes_absolute[0]+nodes_absolute[2],
     180                                                 nodes_absolute[1]+nodes_absolute[2],
     181                                                 nodes_absolute[1]+nodes_absolute[0]]),
     182                                  verts), msg)
     183
     184    def test_get_edge_midpoint_coordinates_triangle_id(self):
     185        """test_get_vertex_coordinates_triangle_id
     186        Test that vertices for one triangle can be returned.
     187        """
     188        from mesh_factory import rectangular
     189
     190        #Create basic mesh
     191        nodes, triangles, _ = rectangular(1, 3)
     192        domain = General_mesh(nodes, triangles)
     193
     194
     195        assert num.allclose(domain.get_nodes(), nodes)
     196
     197
     198        M = domain.number_of_triangles       
     199
     200        for i in range(M):
     201            E = domain.get_edge_midpoint_coordinates(triangle_id=i)
     202            assert E.shape[0] == 3
     203
     204
     205            k0 = triangles[i,0]  #Index of vertex 0 in triangle i
     206            k1 = triangles[i,1]  #Index of vertex 0 in triangle i
     207            k2 = triangles[i,2]  #Index of vertex 0 in triangle i
     208
     209            assert num.allclose(E[0,:], 0.5*(nodes[k1]+nodes[k2]))
     210            assert num.allclose(E[1,:], 0.5*(nodes[k0]+nodes[k2]))
     211            assert num.allclose(E[2,:], 0.5*(nodes[k1]+nodes[k0]))
     212
     213            E0 = domain.get_edge_midpoint_coordinate(i, 0 )
     214            E1 = domain.get_edge_midpoint_coordinate(i, 1 )
     215            E2 = domain.get_edge_midpoint_coordinate(i, 2 )
     216
     217            assert num.allclose(E0, 0.5*(nodes[k1]+nodes[k2]))
     218            assert num.allclose(E1, 0.5*(nodes[k0]+nodes[k2]))
     219            assert num.allclose(E2, 0.5*(nodes[k1]+nodes[k0]))
     220           
    120221
    121222
     
    356457
    357458if __name__ == "__main__":
    358     suite = unittest.makeSuite(Test_General_Mesh, 'test')
     459    #suite = unittest.makeSuite(Test_General_Mesh, 'test')
     460    suite = unittest.makeSuite(Test_General_Mesh, 'test_get_edge_midpoint_coordinates')     
    359461    runner = unittest.TextTestRunner()
    360462    runner.run(suite)
Note: See TracChangeset for help on using the changeset viewer.