Vector check

mint/tests/test_polyline_integral.py

@@ -1,20 +1,32 @@
 from mint import Grid, PolylineIntegral
 import numpy
-import sys
-from pathlib import Path
+import pytest
 
 
 def potentialFunc(p):
     x, y = p[:2]
     return x + 2*y
 
 
-def test_simple():
+def singularPotentialFunc(p):
+    x, y = p[:2]
+    return numpy.arctan2(y, x)/(2.*numpy.pi)
+
+
+@pytest.mark.parametrize("nx", [3])
+@pytest.mark.parametrize("ny", [2])
+@pytest.mark.parametrize("potFunc", [potentialFunc, singularPotentialFunc])
+@pytest.mark.parametrize("xyz", [numpy.array([(1., 0., 0.),
+                                              (0., 1., 0.)]),
+                                 numpy.array([(0., 0., 0.),
+                                              (1., 0., 0.),
+                                              (1., 1., 0.),
+                                              (0., 1., 0.)])])
+def test_simple(nx, ny, potFunc, xyz):
 
     # create the grid and the edge data
     gr = Grid()
 
-    nx, ny = 3, 2
     points = numpy.zeros((nx*ny, 4, 3), numpy.float64)
     data = numpy.zeros((nx*ny, 4))
     dx = 1.0 / float(nx)
@@ -46,40 +58,36 @@ def test_simple():
             #  |     |
             #  +-->--+
             #     0
-            data[k, 0] = potentialFunc(points[k, 1, :]) - potentialFunc(points[k, 0, :])
-            data[k, 1] = potentialFunc(points[k, 2, :]) - potentialFunc(points[k, 1, :])
-            data[k, 2] = potentialFunc(points[k, 2, :]) - potentialFunc(points[k, 3, :])
-            data[k, 3] = potentialFunc(points[k, 3, :]) - potentialFunc(points[k, 0, :])
+            data[k, 0] = potFunc(points[k, 1, :]) - potFunc(points[k, 0, :])
+            data[k, 1] = potFunc(points[k, 2, :]) - potFunc(points[k, 1, :])
+            data[k, 2] = potFunc(points[k, 2, :]) - potFunc(points[k, 3, :])
+            data[k, 3] = potFunc(points[k, 3, :]) - potFunc(points[k, 0, :])
 
             # increment the cell counter
             k += 1
 
     gr.setPoints(points)
     gr.dump('test_polyline_integral.vtk')
 
-
     pli = PolylineIntegral()
-
-    # create the polyline through which the flux will be integrated
-    xyz = numpy.array([(0., 0., 0.),
-                       (1., 0., 0.),
-                       (1., 1., 0.),
-                       (0., 1., 0.)])
-
     # no periodicity in x
     pli.build(gr, xyz, counterclock=False, periodX=0.0)
 
     flux = pli.getIntegral(data)
-    exactFlux = potentialFunc(xyz[-1, :]) - potentialFunc(xyz[0, :])
+    exactFlux = potFunc(xyz[-1, :]) - potFunc(xyz[0, :])
     print(f'total flux: {flux:.3f} exact flux: {exactFlux:.3f}')
     assert abs(flux - exactFlux) < 1.e-10
 
-def test_partially_outside():
 
+@pytest.mark.parametrize("nx", [3])
+@pytest.mark.parametrize("ny", [2])
+@pytest.mark.parametrize("potFunc", [potentialFunc])
+def test_partially_outside(nx, ny, potFunc):
+
+    print('target line is partially outside the domain, expect a warning!')
     # create the grid and the edge data
     gr = Grid()
 
-    nx, ny = 3, 2
     points = numpy.zeros((nx*ny, 4, 3), numpy.float64)
     data = numpy.zeros((nx*ny, 4))
     dx = 1.0 / float(nx)
@@ -111,10 +119,10 @@ def test_partially_outside():
             #  |     |
             #  +-->--+
             #     0
-            data[k, 0] = potentialFunc(points[k, 1, :]) - potentialFunc(points[k, 0, :])
-            data[k, 1] = potentialFunc(points[k, 2, :]) - potentialFunc(points[k, 1, :])
-            data[k, 2] = potentialFunc(points[k, 2, :]) - potentialFunc(points[k, 3, :])
-            data[k, 3] = potentialFunc(points[k, 3, :]) - potentialFunc(points[k, 0, :])
+            data[k, 0] = potFunc(points[k, 1, :]) - potFunc(points[k, 0, :])
+            data[k, 1] = potFunc(points[k, 2, :]) - potFunc(points[k, 1, :])
+            data[k, 2] = potFunc(points[k, 2, :]) - potFunc(points[k, 3, :])
+            data[k, 3] = potFunc(points[k, 3, :]) - potFunc(points[k, 0, :])
 
             # increment the cell counter
             k += 1
@@ -135,19 +143,23 @@ def test_partially_outside():
     flux = pli.getIntegral(data)
 
     # because the first point is outside the domain, only the contribution
-    # stemming from the path inside the domain will be computed. Let's 
+    # stemming from the path inside the domain will be computed. Let's
     # correct for this by moving the first point inwards
     xyz[0, 0] = 0.
     exactFlux = potentialFunc(xyz[-1, :]) - potentialFunc(xyz[0, :])
     print(f'total flux: {flux:.3f} exact flux: {exactFlux:.3f}')
     assert abs(flux - exactFlux) < 1.e-10
 
-def test_completely_outside():
 
+@pytest.mark.parametrize("nx", [3])
+@pytest.mark.parametrize("ny", [2])
+@pytest.mark.parametrize("potFunc", [potentialFunc])
+def test_completely_outside(nx, ny, potFunc):
+
+    print('target line is outside the domain, expect warnings!')
     # create the grid and the edge data
     gr = Grid()
 
-    nx, ny = 3, 2
     points = numpy.zeros((nx*ny, 4, 3), numpy.float64)
     data = numpy.zeros((nx*ny, 4))
     dx = 1.0 / float(nx)
@@ -179,10 +191,10 @@ def test_completely_outside():
             #  |     |
             #  +-->--+
             #     0
-            data[k, 0] = potentialFunc(points[k, 1, :]) - potentialFunc(points[k, 0, :])
-            data[k, 1] = potentialFunc(points[k, 2, :]) - potentialFunc(points[k, 1, :])
-            data[k, 2] = potentialFunc(points[k, 2, :]) - potentialFunc(points[k, 3, :])
-            data[k, 3] = potentialFunc(points[k, 3, :]) - potentialFunc(points[k, 0, :])
+            data[k, 0] = potFunc(points[k, 1, :]) - potFunc(points[k, 0, :])
+            data[k, 1] = potFunc(points[k, 2, :]) - potFunc(points[k, 1, :])
+            data[k, 2] = potFunc(points[k, 2, :]) - potFunc(points[k, 3, :])
+            data[k, 3] = potFunc(points[k, 3, :]) - potFunc(points[k, 0, :])
 
             # increment the cell counter
             k += 1
@@ -208,8 +220,16 @@ def test_completely_outside():
 
 if __name__ == '__main__':
 
-    test_simple()
-    test_partially_outside()
-    test_completely_outside()
+    # polyline through which the line integral will be computed
+    xyz = numpy.array([(1., 0., 0.),
+                       (0., 1., 0.)])
+    test_simple(3, 2, singularPotentialFunc, xyz)
 
+    xyz = numpy.array([(0., 0., 0.),
+                       (1., 0., 0.),
+                       (1., 1., 0.),
+                       (0., 1., 0.)])
+    test_simple(3, 2, potentialFunc, xyz)
 
+    test_partially_outside(2, 3, potentialFunc)
+    test_completely_outside(2, 3, potentialFunc)

mint/tests/test_vector_interp.py

@@ -128,11 +128,17 @@ def test_rectilinear():
             # set one edge to 1, all other edges to zero
             data[cellId, edgeIndex] = 1.0
 
-            # get the interpolated vectors
-            vectorData = vi.getVectors(numpy.array(data))
+            # get the edge interpolated vectors
+            vectorData = vi.getEdgeVectors(data)
             assert(abs(vectorData.max() - 1.) < 1.e-12)
             assert(abs(vectorData.min() - 0.) < 1.e-12)
 
+            # get the lateral flux interpolated vectors
+            vectorData = vi.getFaceVectors(data)
+            # face vectors take the sign of the area vector, negative if pointing down
+            assert(abs(numpy.fabs(vectorData).max() - 1.) < 1.e-12)
+            assert(abs(numpy.fabs(vectorData).min() - 0.) < 1.e-12)            
+
             # reset this edge's value back to its original
             data[cellId, edgeIndex] = 0.0
 
@@ -174,15 +180,19 @@ def test_slanted():
             # set one edge to 1, all other edges to zero
             data[cellId, edgeIndex] = 1.0
 
-            # get the interpolated vectors
-            vectorData = vi.getVectors(numpy.array(data))
+            # get the edge interpolated vectors
+            vectorData = vi.getEdgeVectors(data)
+            fileName = f'slanted_edgeVectors_cellId{cellId}edgeIndex{edgeIndex}.vtk'
+            saveVectorsVTKFile(targetPoints, vectorData, fileName)
+
+            # get the lateral face interpolated vectors
+            vectorData = vi.getFaceVectors(data)
+            fileName = f'slanted_faceVectors_cellId{cellId}edgeIndex{edgeIndex}.vtk'
+            saveVectorsVTKFile(targetPoints, vectorData, fileName)
 
             # reset this edge's value back to its original
             data[cellId, edgeIndex] = 0.0
 
-            fileName = f'slanted_cellId{cellId}edgeIndex{edgeIndex}.vtk'
-            saveVectorsVTKFile(targetPoints, vectorData, fileName)
-
 
 def test_degenerate():
 
@@ -221,8 +231,11 @@ def test_degenerate():
             # set one edge to 1, all other edges to zero
             data[cellId, edgeIndex] = 1.0
 
-            # get the interpolated vectors
-            vectorData = vi.getVectors(numpy.array(data))
+            # get the edge interpolated vectors
+            vectorData = vi.getEdgeVectors(data)
+
+            # get the face interpolated vectors
+            vectorData = vi.getFaceVectors(data)
 
             # reset this edge's value back to its original
             data[cellId, edgeIndex] = 0.0

mint/vector_interp.py

@@ -5,14 +5,20 @@
 import logging
 
 
+FILE = 'vectorinterp.py'
+DOUBLE_ARRAY_PTR = numpy.ctypeslib.ndpointer(dtype=numpy.float64)
+
+
 def error_handler(filename, methodname, ier):
     msg = f'ier={ier} after calling {methodname} in {filename}!'
     logging.error(msg)
     raise RuntimeError(msg)
 
 
-FILE = 'vectorinterp.py'
-DOUBLE_ARRAY_PTR = numpy.ctypeslib.ndpointer(dtype=numpy.float64)
+def warning_handler(filename, methodname, ier, detailedMsg=''):
+    msg = f'ier={ier} after calling {methodname} in {filename}!\n'
+    msg += detailedMsg
+    logging.warning(msg)
 
 
 class VectorInterp(object):
@@ -95,19 +101,38 @@ def findPoints(self, targetPoints, tol2):
                                                  targetPoints, tol2)
         return numBad
 
-    def getVectors(self, data):
+    def getEdgeVectors(self, data):
         """
-        Get the vectors at the target points.
+        Get the edge vectors at given target points.
 
         :param data: edge data array of size numCells times 4
         :returns vector array of size numTargetPoints times 3
-        :note: call this after invoking findPoints
+        :note: call this after invoking findPoints.
         """
-        LIB.mnt_vectorinterp_getVectors.argtypes = [POINTER(c_void_p),
-                                                    DOUBLE_ARRAY_PTR,
-                                                    DOUBLE_ARRAY_PTR]
+        LIB.mnt_vectorinterp_getEdgeVectors.argtypes = [POINTER(c_void_p),
+                                                        DOUBLE_ARRAY_PTR,
+                                                        DOUBLE_ARRAY_PTR]
+        res = numpy.zeros((self.numTargetPoints, 3), numpy.float64)
+        ier = LIB.mnt_vectorinterp_getEdgeVectors(self.obj, data, res)
+        if ier:
+            msg = "Some target lines fall outside the grid."
+            warning_handler(FILE, 'getEdgeVectors', ier, detailedMsg=msg)
+        return res
+
+    def getFaceVectors(self, data):
+        """
+        Get the lateral face vectors at given target points.
+
+        :param data: edge data array of size numCells times 4
+        :returns vector array of size numTargetPoints times 3
+        :note: call this after invoking findPoints.
+        """
+        LIB.mnt_vectorinterp_getFaceVectors.argtypes = [POINTER(c_void_p),
+                                                        DOUBLE_ARRAY_PTR,
+                                                        DOUBLE_ARRAY_PTR]
         res = numpy.zeros((self.numTargetPoints, 3), numpy.float64)
-        ier = LIB.mnt_vectorinterp_getVectors(self.obj, data, res)
+        ier = LIB.mnt_vectorinterp_getFaceVectors(self.obj, data, res)
         if ier:
-            error_handler(FILE, 'getVectors', ier)
+            msg = "Some target lines fall outside the grid."
+            warning_handler(FILE, 'getFaceVectors', ier, detailedMsg=msg)
         return res