Merge pull request #33 from UXARRAY/integration

rajeeja · web-flow · commit b003f30134d2 · 2022-07-22T13:24:36.000-05:00
Integration. Thanks team!
diff --git a/.gitignore b/.gitignore
@@ -142,7 +142,6 @@ test/dask-worker-space/
 .vscode/
 generated/
 
-
 # Test grid files.
 test/meshfiles/*.g
 test/meshfiles/*.ug
diff --git a/ci/environment.yml b/ci/environment.yml
@@ -4,9 +4,13 @@ channels:
 dependencies:
   - dask
   - netcdf4
+  - numba
+  # pinning the numpy version for numba compatibility
+  - numpy>=1.2, <1.22
   - pathlib
   - pre_commit
   - pytest
   - pytest-cov
   - python
+  - scipy
   - xarray
diff --git a/docs/internal_api/index.rst b/docs/internal_api/index.rst
@@ -34,4 +34,5 @@ Grid Helper Modules
    _ugrid._read_ugrid
    _scrip._read_scrip
    _scrip._to_ugrid
+   helpers._spherical_to_cartesian_unit_
    helpers._is_ugrid
diff --git a/docs/user_api/index.rst b/docs/user_api/index.rst
@@ -21,10 +21,19 @@ Grid Methods
    :toctree: _autosummary
 
    grid.Grid.write
+   grid.Grid.calculate_total_face_area
+   grid.Grid.integrate
+   grid.Grid.compute_face_areas
 
 Helper Functions
 ----------------
 .. autosummary::
    :toctree: _autosummary
 
-   parse_grid_type
+   helpers.get_all_face_area_from_coords
+   helpers.calculate_face_area
+   helpers.calculate_spherical_triangle_jacobian
+   helpers.calculate_spherical_triangle_jacobian_barycentric
+   get_gauss_quadratureDG
+   get_tri_quadratureDG
+   helpers.parse_grid_type
diff --git a/requirements.txt b/requirements.txt
@@ -1 +1,2 @@
 xarray
+numpy>=1.2,<=1.22
diff --git a/test/__init__.py b/test/__init__.py
@@ -0,0 +1 @@
+from . import constants
diff --git a/test/constants.py b/test/constants.py
@@ -4,3 +4,8 @@
 NNODES_outCSne30 = 5402
 NNODES_outRLL1deg = 64442
 DATAVARS_outCSne30 = 2
+TRI_AREA = 1.047
+# 4*Pi is 12.56
+MESH30_AREA = 12.566
+PSI_INTG = 12.566
+VAR2_INTG = 12.566
diff --git a/test/test_dataset.py b/test/test_dataset.py
@@ -3,7 +3,11 @@
 from pathlib import Path
 
 import uxarray as ux
-from . import constants
+
+try:
+    import constants
+except ImportError:
+    from . import constants
 
 current_path = Path(os.path.dirname(os.path.realpath(__file__)))
 
diff --git a/test/test_grid.py b/test/test_grid.py
@@ -75,7 +75,7 @@ def test_read_shpfile(self):
         """Reads a shape file and write ugrid file."""
         with self.assertRaises(RuntimeError):
             shp_filename = current_path / "meshfiles" / "grid_fire.shp"
-            tgrid = ux.Grid(str(shp_filename))
+            tgrid = ux.open_dataset(str(shp_filename))
 
     def test_read_scrip(self):
         """Reads a scrip file and write ugrid file."""
diff --git a/test/test_helpers.py b/test/test_helpers.py
@@ -0,0 +1,62 @@
+import os
+import numpy as np
+import numpy.testing as nt
+
+from unittest import TestCase
+from pathlib import Path
+
+import uxarray as ux
+
+try:
+    import constants
+except ImportError:
+    from . import constants
+
+# Data files
+current_path = Path(os.path.dirname(os.path.realpath(__file__)))
+
+exodus = current_path / "meshfiles" / "outCSne8.g"
+
+
+class TestIntegrate(TestCase):
+
+    def test_face_area_coords(self):
+        """Test function for helper function get_all_face_area_from_coords."""
+        # Note: currently only testing one face, but this can be used to get area of multiple faces
+        x = np.array([0.57735027, 0.57735027, -0.57735027])
+        y = np.array([-5.77350269e-01, 5.77350269e-01, 5.77350269e-01])
+        z = np.array([-0.57735027, -0.57735027, -0.57735027])
+        face_nodes = np.array([[0, 1, 2]])
+        area = ux.get_all_face_area_from_coords(x,
+                                                y,
+                                                z,
+                                                face_nodes,
+                                                3,
+                                                coords_type="cartesian")
+        nt.assert_almost_equal(area, constants.TRI_AREA, decimal=1)
+
+    def test_calculate_face_area(self):
+        """Test function for helper function calculate_face_area - only one face."""
+        # Note: currently only testing one face, but this can be used to get area of multiple faces
+        # Also note, this does not need face_nodes, assumes nodes are in counterclockwise orientation
+        x = np.array([0.57735027, 0.57735027, -0.57735027])
+        y = np.array([-5.77350269e-01, 5.77350269e-01, 5.77350269e-01])
+        z = np.array([-0.57735027, -0.57735027, -0.57735027])
+        area = ux.calculate_face_area(x, y, z, "gaussian", 5, "cartesian")
+        nt.assert_almost_equal(area, constants.TRI_AREA, decimal=3)
+
+    def test_quadrature(self):
+        order = 1
+        dG, dW = ux.get_tri_quadratureDG(order)
+        G = np.array([[0.33333333, 0.33333333, 0.33333333]])
+        W = np.array([1.0])
+
+        np.testing.assert_array_almost_equal(G, dG)
+        np.testing.assert_array_almost_equal(W, dW)
+
+        dG, dW = ux.get_gauss_quadratureDG(order)
+        G = np.array([[0.5]])
+        W = np.array([1.0])
+
+        np.testing.assert_array_almost_equal(G, dG)
+        np.testing.assert_array_almost_equal(W, dW)
diff --git a/test/test_integrate.py b/test/test_integrate.py
@@ -0,0 +1,65 @@
+import os
+import numpy as np
+import numpy.testing as nt
+
+from unittest import TestCase
+from pathlib import Path
+
+import uxarray as ux
+
+try:
+    import constants
+except ImportError:
+    from . import constants
+
+# Data files
+current_path = Path(os.path.dirname(os.path.realpath(__file__)))
+
+mesh_file30 = current_path / "meshfiles" / "outCSne30.ug"
+data_file30 = current_path / "meshfiles" / "outCSne30_vortex.nc"
+data_file30_v2 = current_path / "meshfiles" / "outCSne30_var2.ug"
+
+
+class TestIntegrate(TestCase):
+
+    def test_calculate_total_face_area_triangle(self):
+        """Create a uxarray grid from vertices and saves an exodus file."""
+        verts = np.array([[0.57735027, -5.77350269e-01, -0.57735027],
+                          [0.57735027, 5.77350269e-01, -0.57735027],
+                          [-0.57735027, 5.77350269e-01, -0.57735027]])
+        vgrid = ux.Grid(verts)
+
+        # get node names for each grid object
+        x_var = vgrid.ds_var_names["Mesh2_node_x"]
+        y_var = vgrid.ds_var_names["Mesh2_node_y"]
+        z_var = vgrid.ds_var_names["Mesh2_node_z"]
+
+        vgrid.ds[x_var].attrs["units"] = "m"
+        vgrid.ds[y_var].attrs["units"] = "m"
+        vgrid.ds[z_var].attrs["units"] = "m"
+
+        area_gaussian = vgrid.calculate_total_face_area(
+            quadrature_rule="gaussian", order=5)
+        nt.assert_almost_equal(area_gaussian, constants.TRI_AREA, decimal=3)
+
+        area_triangular = vgrid.calculate_total_face_area(
+            quadrature_rule="triangular", order=4)
+        nt.assert_almost_equal(area_triangular, constants.TRI_AREA, decimal=1)
+
+    def test_calculate_total_face_area_file(self):
+        """Create a uxarray grid from vertices and saves an exodus file."""
+
+        grid = ux.open_dataset(str(mesh_file30))
+
+        area = grid.calculate_total_face_area()
+
+        nt.assert_almost_equal(area, constants.MESH30_AREA, decimal=3)
+
+    def test_integrate(self):
+        uds = ux.open_dataset(mesh_file30, data_file30, data_file30_v2)
+
+        integral_psi = uds.integrate("psi")
+        integral_var2 = uds.integrate("var2")
+
+        nt.assert_almost_equal(integral_psi, constants.PSI_INTG, decimal=3)
+        nt.assert_almost_equal(integral_var2, constants.VAR2_INTG, decimal=3)
diff --git a/test/test_ugrid.py b/test/test_ugrid.py
@@ -5,7 +5,11 @@
 from pathlib import Path
 
 import uxarray as ux
-from . import constants
+
+try:
+    import constants
+except ImportError:
+    from . import constants
 
 current_path = Path(os.path.dirname(os.path.realpath(__file__)))
 
diff --git a/uxarray/__init__.py b/uxarray/__init__.py
@@ -1,2 +1,3 @@
 from .grid import *
 from .dataset import *
+from .helpers import *
diff --git a/uxarray/_exodus.py b/uxarray/_exodus.py
@@ -53,7 +53,7 @@ def _read_exodus(ext_ds, ds_var_names):
                 attrs={
                     "standard_name": "longitude",
                     "long_name": "longitude of mesh nodes",
-                    "units": "degress_east",
+                    "units": "degrees_east",
                 })
             ds["Mesh2_node_y"] = xr.DataArray(
                 data=ext_ds.coord[1],
@@ -79,7 +79,7 @@ def _read_exodus(ext_ds, ds_var_names):
                 attrs={
                     "standard_name": "longitude",
                     "long_name": "longitude of mesh nodes",
-                    "units": "degress_east",
+                    "units": "degrees_east",
                 })
         elif key == "coordy":
             ds["Mesh2_node_y"] = xr.DataArray(
@@ -102,21 +102,14 @@ def _read_exodus(ext_ds, ds_var_names):
                     })
         elif "connect" in key:
             # check if num face nodes is less than max.
-            if value.data.shape[1] < max_face_nodes:
-                # create a temporary array to store connectivity
-                tmp_conn = np.empty((value.data.shape[0], max_face_nodes))
-                tmp_conn.fill(
-                    0
-                )  # exodus in 1-based, fill with zeros here; during assignment subtract 1
-                tmp_conn[:value.data.shape[0], :value.data.
-                         shape[1]] = value.data
-            elif value.data.shape[1] == max_face_nodes:
-                tmp_conn = value.data
+            if value.data.shape[1] <= max_face_nodes:
+                conn = np.full((value.data.shape[1], max_face_nodes),
+                               0,
+                               dtype=conn.dtype)
+                conn = value.data
             else:
                 raise "found face_nodes_dim greater than nMaxMesh2_face_nodes"
 
-            # concatenate to the previous blk
-            conn = np.concatenate((conn, tmp_conn))
             # find the elem_type as etype for this element
             for k, v in value.attrs.items():
                 if k == "elem_type":
diff --git a/uxarray/dataset.py b/uxarray/dataset.py
@@ -19,7 +19,7 @@ def open_dataset(grid_filename, *args, **kw):
     -------
 
     object : xarray Dataset
-        Contains the grid and corresponding data
+        uxarray grid object: Contains the grid and corresponding data
 
     Examples
     --------
diff --git a/uxarray/get_quadratureDG.py b/uxarray/get_quadratureDG.py
diff --git a/uxarray/grid.py b/uxarray/grid.py
diff --git a/uxarray/helpers.py b/uxarray/helpers.py
