 Timestamp:
 Jun 18, 2010, 5:49:57 PM (13 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

anuga_work/development/2010projects/anuga_1d/base/generic_domain.py
r7842 r7860 40 40 self.set_timestepping_method(timestepping_method) 41 41 42 self.wet_nodes = numpy.zeros((N,2), numpy.int) # should this be here 42 43 # work array for special vertices (like wet or dry cells 44 self.special_vertices = numpy.zeros((N,2), numpy.int) # should this be here 43 45 44 46 #Allocate space for neighbour and boundary structures … … 58 60 59 61 self.normals = numpy.zeros((N, 2), numpy.float) 60 61 for i in range(N): 62 xl = self.coordinates[i] 63 xr = self.coordinates[i+1] 64 self.vertices[i,0] = xl 65 self.vertices[i,1] = xr 66 67 centroid = (xl+xr)/2.0 68 self.centroids[i] = centroid 69 70 msg = 'Coordinates should be ordered, smallest to largest' 71 assert xr>xl, msg 72 73 #The normal vectors 74 #  point outward from each edge 75 #  are orthogonal to the edge 76 #  have unit length 77 #  Are enumerated by left vertex then right vertex normals 78 79 nl = 1.0 80 nr = 1.0 81 self.normals[i,:] = [nl, nr] 82 83 self.areas[i] = (xrxl) 84 85 # # print 'N', N 86 # # print 'Centroid', self.centroids 87 # # print 'Areas', self.areas 88 # # print 'Vertex_Coordinates', self.vertices 89 90 #Initialise Neighbours (1 means that it is a boundary neighbour) 91 self.neighbours[i, :] = [1, 1] 92 #Initialise edge ids of neighbours 93 #Initialise vertex ids of neighbours 94 #In case of boundaries this slot is not used 95 #self.neighbour_edges[i, :] = [1, 1] 96 self.neighbour_vertices[i, :] = [1, 1] 62 63 64 xl = self.coordinates[0:1] 65 xr = self.coordinates[1: ] 66 67 self.vertices[:,0] = xl 68 self.vertices[:,1] = xr 69 70 self.centroids[:] = (xl+xr)/2.0 71 72 msg = 'Coordinates should be ordered, smallest to largest' 73 assert (xr>xl).all(), msg 74 75 #The normal vectors 76 #  point outward from each edge 77 #  are orthogonal to the edge 78 #  have unit length 79 #  Are enumerated by left vertex then right vertex normals 80 81 82 self.normals[:,0] = 1.0 83 self.normals[:,1] = 1.0 84 85 self.areas[:] = (xrxl) 86 87 #Initialise Neighbours (1 means that it is a boundary neighbour) 88 89 self.neighbours[:, :] = 1 90 91 #Initialise edge ids of neighbours 92 self.neighbour_vertices[:, :] = 1 93 94 95 # for i in range(N): 96 # xl = self.coordinates[i] 97 # xr = self.coordinates[i+1] 98 # self.vertices[i,0] = xl 99 # self.vertices[i,1] = xr 100 # 101 # centroid = (xl+xr)/2.0 102 # self.centroids[i] = centroid 103 # 104 # msg = 'Coordinates should be ordered, smallest to largest' 105 # assert xr>xl, msg 106 # 107 # #The normal vectors 108 # #  point outward from each edge 109 # #  are orthogonal to the edge 110 # #  have unit length 111 # #  Are enumerated by left vertex then right vertex normals 112 # 113 # nl = 1.0 114 # nr = 1.0 115 # self.normals[i,:] = [nl, nr] 116 # 117 # self.areas[i] = (xrxl) 118 # 119 ## # print 'N', N 120 ## # print 'Centroid', self.centroids 121 ## # print 'Areas', self.areas 122 ## # print 'Vertex_Coordinates', self.vertices 123 # 124 # #Initialise Neighbours (1 means that it is a boundary neighbour) 125 # self.neighbours[i, :] = [1, 1] 126 # #Initialise edge ids of neighbours 127 # #Initialise vertex ids of neighbours 128 # #In case of boundaries this slot is not used 129 # #self.neighbour_edges[i, :] = [1, 1] 130 # self.neighbour_vertices[i, :] = [1, 1] 97 131 98 132 self.build_vertexlist() … … 630 664 return 631 665 666 if timestepping_method in [1, 2, 3]: 667 self.timestepping_method = ['euler', 'rk2', 'rk3'][timestepping_method1] 668 return 669 632 670 msg = '%s is an incorrect timestepping type'% timestepping_method 633 671 raise Exception, msg … … 814 852 return self.filename 815 853 854 855 def get_spatial_order(self): 856 return self.order 857 858 def set_spatial_order(self, order): 859 860 possible_orders = [1,2] 861 862 if order in possible_orders: 863 self.order = order 864 return 865 866 msg = '%s is an incorrect spatial order.\n'% limiter 867 msg += 'Possible orders are: '+ ", ".join(["%s" % el for el in possible_orders]) 868 raise Exception, msg 869 870 816 871 def get_limiter(self): 817 872 return self.limiter
Note: See TracChangeset
for help on using the changeset viewer.