test: updating test for unstructured

Author: Lachlan Grose

tests/unit_tests/interpolator/test_unstructured_supports.py

@@ -15,27 +15,37 @@ def test_get_elements():
 
     vertices = nodes[elements,:]
     pos = points[:,:]
-    vap = pos[None,:, :] - vertices[:, None, 0, :]
-    vbp = pos[None,:, :] - vertices[:, None, 1, :]
-
-    vab = vertices[:, None, 1, :] - vertices[:,None,  0, :]
-    vac = vertices[:, None, 2, :] - vertices[:,None,  0, :]
-    vad = vertices[:, None, 3, :] - vertices[:,None,  0, :]
-    vbc = vertices[:, None, 2, :] - vertices[:,None,  1, :]
-    vbd = vertices[:, None, 3, :] - vertices[:, None, 1, :]
+    vap = pos[:,None, :] - vertices[None, :, 0, :]
+    vbp = pos[:,None, :] - vertices[None, :, 1, :]
+    #         # vcp = p - points[:, 2, :]
+    #         # vdp = p - points[:, 3, :]
+    vab = vertices[None, :, 1, :] - vertices[None, :,  0, :]
+    vac = vertices[None, :, 2, :] - vertices[None, :,  0, :]
+    vad = vertices[None, :, 3, :] - vertices[None, :,  0, :]
+    vbc = vertices[None, :, 2, :] - vertices[None, :,  1, :]
+    vbd = vertices[None, :, 3, :] - vertices[None, :, 1, :]
 
     va = np.einsum('ikj, ikj->ik', vbp, np.cross(vbd, vbc, axisa=2, axisb=2)) / 6.
     vb = np.einsum('ikj, ikj->ik', vap, np.cross(vac, vad, axisa=2, axisb=2)) / 6.
     vc = np.einsum('ikj, ikj->ik', vap, np.cross(vad, vab, axisa=2, axisb=2)) / 6.
     vd = np.einsum('ikj, ikj->ik', vap, np.cross(vab, vac, axisa=2, axisb=2)) / 6.
     v = np.einsum('ikj, ikj->ik', vab, np.cross(vac, vad, axisa=2, axisb=2)) / 6.
-
-    c = np.zeros((va.shape[0], pos.shape[0], 4))
+    # c = np.zeros((va.shape[0], 4))
+    # c[:,  0] = va / v
+    # c[:,  1] = vb / v
+    # c[:,  2] = vc / v
+    # c[:,  3] = vd / v
+    c = np.zeros((pos.shape[0],va.shape[1],  4))
     c[:,:,  0] = va / v
     c[:,:,  1] = vb / v
     c[:,:,  2] = vc / v
+
     c[:,:,  3] = vd / v
     row,col = np.where(np.all(c>=0,axis=2))
-    tetra_idx = np.sort(row)
 
-    assert np.all(elements[tetra_idx]-tetra == 0)
+    tetra_idx = col
+    
+    # check if the calculated tetra from the mesh method using aabb
+    # is the same as using the barycentric coordinates on all elelemts for 
+    # all points
+    assert np.all(elements[tetra_idx]-tetra==0)