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:
Calculate Total Face Area
Options for
Grid.calculate_total_face_area
FunctionGetting Area of Individual Faces
Calculate Area of a Single Triangle in Cartesian Coordinates
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])