Face Area Calculations#

Computing and working with face areas is important for the analysis of unstructured grids, with many algorithms and workflows requiring them. This section will showcase the different face area calculation options provided with UXarray:

  1. Calculate Total Face Area

  2. Options for Grid.calculate_total_face_area Function

  3. Getting Area of Individual Faces

  4. Calculate Area of a Single Triangle in Cartesian Coordinates

  5. Calculate Area from Multiple Faces in Spherical Coordinates

import uxarray as ux
import numpy as np

We will be using the outCSne30.ug grid file, which is encoded in the UGRID convention.

base_path = "../../test/meshfiles/"
grid_path = base_path + "/ugrid/outCSne30/outCSne30.ug"

ugrid = ux.open_grid(grid_path)
ugrid
<uxarray.Grid>
Original Grid Type: UGRID
Grid Dimensions:
  * n_node: 5402
  * n_face: 5400
  * n_max_face_nodes: 4
  * n_nodes_per_face: (5400,)
Grid Coordinates (Spherical):
  * node_lon: (5402,)
  * node_lat: (5402,)
Grid Coordinates (Cartesian):
Grid Connectivity Variables:
  * face_node_connectivity: (5400, 4)
Grid Descriptor Variables:
  * n_nodes_per_face: (5400,)

1. Calculate Total Face Area#

We can calculate the total face area by calling the function Grid.calculate_total_face_area(). Since our dataset lies on the unit sphere, our expected area is 4π, which is approximately 12.56

t4_area = ugrid.calculate_total_face_area()
t4_area
np.float64(12.566370614678554)

2. Options for Grid.calculate_total_face_area Function#

By default, Grid.calculate_total_face_area uses a Triangular Quadrature Rule and an Order of 4. However, you can specify the Quadrature Rule and Order as follows:

Order:

   1 to 10              for gaussian

   1, 4, 8, 10 and 12   for triangular
t1_area = ugrid.calculate_total_face_area(quadrature_rule="triangular", order=1)
t1_area
np.float64(12.571403993719983)

For the result above, notice that the area is slightly different than the first calculation we made.

Now we use Triangular Quadrature Rule and Order 1

Using a lower order is faster, but at the sacrifice of accuracy.

Generally, gaussian quadrature rule is more accurate than the triangular quadrature rule. Additionally, a higher order comes at the cost of computation time, but produces a more accurate result. See uxarray/grid/area.py file and function get_gauss_quadratureDG for details on quadrature points and weights.

3. Getting Area of Individual Faces#

We can access the Grid attribute Grid.face_area to access the individual face areas. If you have not run a face area calculation yet, it will run the Grid.compute_face_areas and cache the value.

ugrid.face_areas
<xarray.DataArray 'face_areas' (n_face: 5400)> Size: 43kB
array([0.00211174, 0.00211221, 0.00210723, ..., 0.00210723, 0.00211221,
       0.00211174])
Dimensions without coordinates: n_face
Attributes:
    cf_role:    face_areas
    long_name:  Area of each face.

Now calculate the area again with the Grid.compute_face_areas function using arguments: quadrature_rule “gaussian” and order 4

all_face_areas, all_face_jacobians = ugrid.compute_face_areas(
    quadrature_rule="gaussian", order=4
)
g4_area = all_face_areas.sum()
g4_area
np.float64(12.566370614359112)

Now we compare the values with actual know value and report error for each of the three cases above.

actual_area = 4 * np.pi
diff_t4_area = np.abs(t4_area - actual_area)
diff_t1_area = np.abs(t1_area - actual_area)
diff_g4_area = np.abs(g4_area - actual_area)

diff_t1_area, diff_t4_area, diff_g4_area
(np.float64(0.005033379360810386),
 np.float64(3.1938185429680743e-10),
 np.float64(6.039613253960852e-14))

As we can see, it is clear that the Gaussian Quadrature Rule with Order 4 is the most accurate, and the Triangular Quadrature Rule with Order 1 is the least accurate.

4. Calculate Area of a Single Triangle in Cartesian Coordinates#

For this section, we create a single triangle face with 3 vertices. By default, in uxarray, we assume that the coordinate system is spherical (lat / lon), however if you want to use cartesian coordinates, you must pass through latlon = False into the Grid constructor.

Assume the units in meters - this is a big triangle!

verts = [
    [
        [0.57735027, -5.77350269e-01, -0.57735027],
        [0.57735027, 5.77350269e-01, -0.57735027],
        [-0.57735027, 5.77350269e-01, -0.57735027],
    ]
]

# load our vertices into a UXarray Grid object
vgrid = ux.open_grid(verts, latlon=False)

vgrid
<uxarray.Grid>
Original Grid Type: Face Vertices
Grid Dimensions:
  * n_node: 3
  * n_face: 1
  * n_max_face_nodes: 3
  * n_nodes_per_face: (1,)
Grid Coordinates (Spherical):
Grid Coordinates (Cartesian):
  * node_x: (3,)
  * node_y: (3,)
  * node_z: (3,)
Grid Connectivity Variables:
  * face_node_connectivity: (1, 3)
Grid Descriptor Variables:
  * n_nodes_per_face: (1,)
vgrid.calculate_total_face_area()
np.float64(1.0719419938548218)

Additionally, if you are using a unit other than meters, you can update the units as follows

# Calculate the area of the triangle
area_gaussian = vgrid.calculate_total_face_area(quadrature_rule="gaussian", order=5)
area_gaussian
np.float64(1.0475702709086991)

5. Calculate Area from Multiple Faces in Spherical Coordinates#

Similar to above, we can construct a Grid object with multiple faces by passing through a set of vertices. Here we define 3 six-sided faces.

faces_verts_ndarray = np.array(
    [
        np.array(
            [
                [150, 10, 0],
                [160, 20, 0],
                [150, 30, 0],
                [135, 30, 0],
                [125, 20, 0],
                [135, 10, 0],
            ]
        ),
        np.array(
            [
                [125, 20, 0],
                [135, 30, 0],
                [125, 60, 0],
                [110, 60, 0],
                [100, 30, 0],
                [105, 20, 0],
            ]
        ),
        np.array(
            [
                [95, 10, 0],
                [105, 20, 0],
                [100, 30, 0],
                [85, 30, 0],
                [75, 20, 0],
                [85, 10, 0],
            ]
        ),
    ]
)

We want our units to be spherical, so we pass through latlon=True. Additionally, if latlon is not passed through, it will default to spherical coordinates.

verts_grid = ux.open_grid(faces_verts_ndarray, latlon=True)

verts_grid
<uxarray.Grid>
Original Grid Type: Face Vertices
Grid Dimensions:
  * n_node: 14
  * n_face: 3
  * n_max_face_nodes: 6
  * n_nodes_per_face: (3,)
Grid Coordinates (Spherical):
  * node_lon: (14,)
  * node_lat: (14,)
Grid Coordinates (Cartesian):
Grid Connectivity Variables:
  * face_node_connectivity: (3, 6)
Grid Descriptor Variables:
  * n_nodes_per_face: (3,)
area, jacobian = verts_grid.compute_face_areas()
area
array([0.14323746, 0.25118746, 0.12141312])