@@ -1011,6 +1011,41 @@ def face_area(self):
10111011 area_j = np .pi * (r0 * np .sin (t0 ) + r1 * np .sin (t0 )) * np .sqrt (np .square (r1 * np .sin (t0 ) - r0 * np .sin (t0 )) + np .square (r1 * np .cos (t0 ) - r0 * np .cos (t0 )))
10121012 return tr (np .array ([area_i , area_j ]))
10131013
1014+ def area_x (self ):
1015+ """
1016+ Return an array of the face areas.
1017+ The shape of the returned array is (ni, nj).
1018+ """
1019+ tr = lambda arr : arr .transpose (1 , 2 , 0 )
1020+ x = self .cell_vertices ()[:,0 ]
1021+ y = self .cell_vertices ()[0 ,:]
1022+ r0 = x [:- 1 , :- 1 ]
1023+ r1 = x [+ 1 :, :- 1 ]
1024+ t0 = y [:- 1 , :- 1 ]
1025+ t1 = y [+ 1 :, + 1 :]
1026+
1027+ area = r1 - r0
1028+ return area
1029+
1030+ def area_y (self ):
1031+ """
1032+ Return an array of the face areas.
1033+ The shape of the returned array is (ni, nj).
1034+ """
1035+ tr = lambda arr : arr .transpose (1 , 2 , 0 )
1036+ x = self .cell_vertices ()[:,0 ]
1037+ y = self .cell_vertices ()[0 ,:]
1038+ r0 = x [:- 1 , :- 1 ]
1039+ r1 = x [+ 1 :, :- 1 ]
1040+ t0 = y [:- 1 , :- 1 ]
1041+ t1 = y [+ 1 :, + 1 :]
1042+
1043+ # ** the area of a part of an annulus
1044+
1045+ area = 0.5 * (r1 ** 2 - r0 ** 2 ) * (t1 - t0 )
1046+ return area
1047+
1048+
10141049 def cell_volumes (self ):
10151050 """
10161051 Return an array of the cell volume data for the given coordinate box
@@ -1025,8 +1060,7 @@ def cell_volumes(self):
10251060 t0 = y [:- 1 , :- 1 ]
10261061 t1 = y [+ 1 :, + 1 :]
10271062
1028- return (r1 ** 3 - r0 ** 3 ) * (np .cos (t1 ) - np .cos (t0 )) * (- 2.0 * np .pi ) / 3.0
1029- # return r1
1063+ return 0.5 * (r1 ** 2 - r0 ** 2 ) * (t1 - t0 ) * (r1 - r0 )
10301064
10311065 def cell_vertices (self ):
10321066 """
0 commit comments