Repository: kinnala/scikit-fem

Issue:
Add doflocs to all Element subclasses
Remaining work from #223 .

Discussion:

Patch:
diff --git a/skfem/assembly/basis/basis.py b/skfem/assembly/basis/basis.py
index cfd8b19c..e998c5b8 100644
--- a/skfem/assembly/basis/basis.py
+++ b/skfem/assembly/basis/basis.py
@@ -48,6 +48,8 @@ class Basis:
         if hasattr(elem, 'doflocs'):
             doflocs = self.mapping.F(elem.doflocs.T)
             self.doflocs = np.zeros((doflocs.shape[0], self.N))
+
+            # match mapped dofs and global dof numbering
             for itr in range(doflocs.shape[0]):
                 for jtr in range(self.element_dofs.shape[0]):
                     self.doflocs[itr, self.element_dofs[jtr]] =\
diff --git a/skfem/element/element_hex/element_hex1.py b/skfem/element/element_hex/element_hex1.py
index 613c4810..0979ddfe 100644
--- a/skfem/element/element_hex/element_hex1.py
+++ b/skfem/element/element_hex/element_hex1.py
@@ -9,14 +9,14 @@ class ElementHex1(ElementH1):
     dim = 3
     maxdeg = 3
     dofnames = ['u']
-    doflocs = np.array([[0., 0., 0.],
-                        [0., 0., 1.],
-                        [0., 1., 0.],
-                        [1., 0., 0.],
-                        [0., 1., 1.],
-                        [1., 0., 1.],
+    doflocs = np.array([[1., 1., 1.],
                         [1., 1., 0.],
-                        [1., 1., 1.]])
+                        [1., 0., 1.],
+                        [0., 1., 1.],
+                        [1., 0., 0.],
+                        [0., 1., 0.],
+                        [0., 0., 1.],
+                        [0., 0., 0.]])
     mesh_type = MeshHex
 
     def lbasis(self, X, i):
diff --git a/skfem/element/element_line/element_line_pp.py b/skfem/element/element_line/element_line_pp.py
index a50adc20..410551bb 100644
--- a/skfem/element/element_line/element_line_pp.py
+++ b/skfem/element/element_line/element_line_pp.py
@@ -21,7 +21,7 @@ class ElementLinePp(ElementH1):
         self.interior_dofs = p - 1
         self.maxdeg = p
         self.dofnames = ['u'] + (p - 1) * ['u']
-        self.doflocs = np.array([[0., 1.] + [.5] * (p - 1)]).T
+        self.doflocs = np.array([[0., 1.] + [np.nan] * (p - 1)]).T
         self.P = np.zeros((0, 0))
         self.dP = np.zeros((0, 0, 1))
         self.p = p
diff --git a/skfem/element/element_quad/element_quadp.py b/skfem/element/element_quad/element_quadp.py
index c6461ae7..691216b8 100644
--- a/skfem/element/element_quad/element_quadp.py
+++ b/skfem/element/element_quad/element_quadp.py
@@ -16,13 +16,12 @@ class ElementQuadP(ElementLinePp):
         self.interior_dofs = (p + 1) ** 2 - 4 * self.facet_dofs - 4
         self.maxdeg = p ** 2
         self.dofnames = (1 + self.facet_dofs + self.interior_dofs) * ['u']
-        flocs = np.linspace(0, 1, self.facet_dofs + 2)[1:-1] if self.facet_dofs > 0 else []
-        self.doflocs = np.array([[0., 0.], [1., 0.], [1., 1.], [0., 1.]] +\
-                                [[x, 0.] for x in flocs] +\
-                                [[1., y] for y in flocs] +\
-                                [[-x + 1., 1.] for x in flocs] +\
-                                [[0., y] for y in flocs] +\
-                                [[.5, .5] for i in range(self.interior_dofs)])
+        N = 4 * self.facet_dofs + self.interior_dofs
+        self.doflocs = np.array([[0., 0.],
+                                 [1., 0.],
+                                 [1., 1.],
+                                 [0., 1.]] +
+                                [[np.nan, np.nan] for i in range(N)])
         self.Px, self.Py = np.zeros((0, 0)), np.zeros((0, 0))
         self.dPx, self.dPy = np.zeros((0, 0, 1)), np.zeros((0, 0, 1))
         self.p = p
diff --git a/skfem/element/element_tet/element_tet_n0.py b/skfem/element/element_tet/element_tet_n0.py
index 7232a491..a6e3c815 100644
--- a/skfem/element/element_tet/element_tet_n0.py
+++ b/skfem/element/element_tet/element_tet_n0.py
@@ -8,29 +8,47 @@ class ElementTetN0(ElementHcurl):
     dim = 3
     maxdeg = 1
     dofnames = ['u^t']
+    doflocs = np.array([[.5, 0., 0.],
+                        [.5, .5, 0.],
+                        [0., .5, 0.],
+                        [0., 0., .5],
+                        [.5, 0., .5],
+                        [0., .5, .5],])
     mesh_type = MeshTet
 
     def lbasis(self, X, i):
-        x, y, z = X[0, :], X[1, :], X[2, :]
+        x, y, z = X
 
         if i == 0:
-            phi = np.array([1-z-y, x, x])
-            dphi = np.array([0*x, -2 + 0*x, 2 + 0*x])
+            phi = np.array([1 - z - y, x, x])
+            dphi = np.array([0 * x,
+                             -2 + 0 * x,
+                             2 + 0 * x])
         elif i == 1:
-            phi = np.array([-y, x, 0*z])
-            dphi = np.array([0*x, 0*x, 2 + 0*x])
+            phi = np.array([-y, x, 0 * z])
+            dphi = np.array([0 * x,
+                             0 * x,
+                             2 + 0 * x])
         elif i == 2:
-            phi = np.array([y, 1-z-x, y])
-            dphi = np.array([2 + 0*x, 0*x, -2 + 0*x])
+            phi = np.array([y, 1 - z - x, y])
+            dphi = np.array([2 + 0 * x,
+                             0 * x,
+                             -2 + 0 * x])
         elif i == 3:
-            phi = np.array([z, z, 1-x-y])
-            dphi = np.array([-2 + 0*x, 2 + 0*x, 0*x])
+            phi = np.array([z, z, 1 - x - y])
+            dphi = np.array([-2 + 0 * x,
+                             2 + 0 * x,
+                             0 * x])
         elif i == 4:
-            phi = np.array([-z, 0*y, x])
-            dphi = np.array([0*x, -2 + 0*x, 0*x])
+            phi = np.array([-z, 0 * y, x])
+            dphi = np.array([0 * x,
+                             -2 + 0 * x,
+                             0 * x])
         elif i == 5:
-            phi = np.array([0*x, -z, y])
-            dphi = np.array([2 + 0*x, 0*x, 0*x])
+            phi = np.array([0 * x, -z, y])
+            dphi = np.array([2 + 0 * x,
+                             0 * x,
+                             0 * x])
         else:
             self._index_error()
 
diff --git a/skfem/element/element_tet/element_tet_rt0.py b/skfem/element/element_tet/element_tet_rt0.py
index 5a3492c0..e74a6184 100644
--- a/skfem/element/element_tet/element_tet_rt0.py
+++ b/skfem/element/element_tet/element_tet_rt0.py
@@ -9,6 +9,10 @@ class ElementTetRT0(ElementHdiv):
     dim = 3
     maxdeg = 1
     dofnames = ['u^n']
+    doflocs = np.array([[.5, .5, .0],
+                        [.5, .0, .5],
+                        [.0, .5, .5],
+                        [.5, .5, .5],])
     mesh_type = MeshTet
 
     def lbasis(self, X, i):
diff --git a/skfem/element/element_tri/element_tri_mini.py b/skfem/element/element_tri/element_tri_mini.py
index fa611850..d1d1bb13 100644
--- a/skfem/element/element_tri/element_tri_mini.py
+++ b/skfem/element/element_tri/element_tri_mini.py
@@ -13,7 +13,7 @@ class ElementTriMini(ElementH1):
     doflocs = np.array([[0., 0.],
                         [1., 0.],
                         [0., 1.],
-                        [1./3., 1./3.]])
+                        [np.nan, np.nan]])
     mesh_type = MeshTri
 
     def lbasis(self, X, i):
diff --git a/skfem/element/element_tri/element_tri_rt0.py b/skfem/element/element_tri/element_tri_rt0.py
index 5ecc72ff..dc7d81d3 100644
--- a/skfem/element/element_tri/element_tri_rt0.py
+++ b/skfem/element/element_tri/element_tri_rt0.py
@@ -9,9 +9,9 @@ class ElementTriRT0(ElementHdiv):
     dim = 2
     maxdeg = 1
     dofnames = ['u^n']
-    doflocs = np.array([[.5, 0.],
+    doflocs = np.array([[.5, .0],
                         [.5, .5],
-                        [0., .5]])
+                        [.0, .5]])
     mesh_type = MeshTri
 
     def lbasis(self, X, i):
diff --git a/skfem/mapping/mapping_affine.py b/skfem/mapping/mapping_affine.py
index b3aa5d9f..659de44a 100644
--- a/skfem/mapping/mapping_affine.py
+++ b/skfem/mapping/mapping_affine.py
@@ -16,10 +16,10 @@ class MappingAffine(Mapping):
             self.b = np.empty((dim, nt))
 
             for i in range(dim):
-                self.b[i] = mesh.p[i, mesh.t[0, :]]
+                self.b[i] = mesh.p[i, mesh.t[0]]
                 for j in range(dim):
-                    self.A[i, j] = (mesh.p[i, mesh.t[j + 1, :]] -
-                                    mesh.p[i, mesh.t[0, :]])
+                    self.A[i, j] = (mesh.p[i, mesh.t[j + 1]] -
+                                    mesh.p[i, mesh.t[0]])
 
             # determinants
             if dim == 1:
diff --git a/skfem/mesh/mesh2d/mesh_quad.py b/skfem/mesh/mesh2d/mesh_quad.py
index 5e43fc59..7d720b21 100644
--- a/skfem/mesh/mesh2d/mesh_quad.py
+++ b/skfem/mesh/mesh2d/mesh_quad.py
@@ -161,7 +161,7 @@ class MeshQuad(Mesh2D):
 
         # get unique facets and build quad-to-facet mapping: 4 (edges) x Nquads
         tmp = np.ascontiguousarray(self.facets.T)
-        tmp, ixa, ixb = np.unique(tmp.view([('', tmp.dtype)]*tmp.shape[1]),
+        tmp, ixa, ixb = np.unique(tmp.view([('', tmp.dtype)] * tmp.shape[1]),
                                   return_index=True, return_inverse=True)
         self.facets = self.facets[:, ixa]
         self.t2f = ixb.reshape((4, self.t.shape[1]))
diff --git a/skfem/mesh/mesh3d/mesh_hex.py b/skfem/mesh/mesh3d/mesh_hex.py
index cb03e523..e1893d78 100644
--- a/skfem/mesh/mesh3d/mesh_hex.py
+++ b/skfem/mesh/mesh3d/mesh_hex.py
@@ -84,6 +84,11 @@ class MeshHex(Mesh3D):
             self._validate()
         self._build_mappings()
 
+    @classmethod
+    def init_refdom(cls: Type[MeshType]):
+        """Initialise a mesh of the reference domain."""
+        return cls()
+
     @classmethod
     def init_tensor(cls: Type[MeshType],
                     x: ndarray,
diff --git a/skfem/mesh/mesh3d/mesh_tet.py b/skfem/mesh/mesh3d/mesh_tet.py
index 628e0e7f..871c68ce 100644
--- a/skfem/mesh/mesh3d/mesh_tet.py
+++ b/skfem/mesh/mesh3d/mesh_tet.py
@@ -78,6 +78,16 @@ class MeshTet(Mesh3D):
         self.enable_facets = True
         self._build_mappings()
 
+    @classmethod
+    def init_refdom(cls: Type[MeshType]):
+        """Initialise a mesh of the reference domain."""
+        p = np.array([[0., 0., 0.],
+                      [1., 0., 0.],
+                      [0., 1., 0.],
+                      [0., 0., 1.]]).T
+        t = np.array([[0, 1, 2, 3]]).T
+        return cls(p, t)
+
     @classmethod
     def init_tensor(cls: Type[MeshType],
                     x: ndarray,
@@ -378,22 +388,22 @@ class MeshTet(Mesh3D):
         sz = p.shape[1]
         t2e = self.t2e + sz
         # new vertices are the midpoints of edges
-        newp = 0.5*np.vstack((p[0, e[0, :]] + p[0, e[1, :]],
-                              p[1, e[0, :]] + p[1, e[1, :]],
-                              p[2, e[0, :]] + p[2, e[1, :]]))
+        newp = 0.5*np.vstack((p[0, e[0]] + p[0, e[1]],
+                              p[1, e[0]] + p[1, e[1]],
+                              p[2, e[0]] + p[2, e[1]]))
         newp = np.hstack((p, newp))
         # new tets
-        newt = np.vstack((t[0, :], t2e[0, :], t2e[2, :], t2e[3, :]))
-        newt = np.hstack((newt, np.vstack((t[1, :], t2e[0, :], t2e[1, :], t2e[4, :]))))
-        newt = np.hstack((newt, np.vstack((t[2, :], t2e[1, :], t2e[2, :], t2e[5, :]))))
-        newt = np.hstack((newt, np.vstack((t[3, :], t2e[3, :], t2e[4, :], t2e[5, :]))))
+        newt = np.vstack((t[0], t2e[0], t2e[2], t2e[3]))
+        newt = np.hstack((newt, np.vstack((t[1], t2e[0], t2e[1], t2e[4]))))
+        newt = np.hstack((newt, np.vstack((t[2], t2e[1], t2e[2], t2e[5]))))
+        newt = np.hstack((newt, np.vstack((t[3], t2e[3], t2e[4], t2e[5]))))
         # compute middle pyramid diagonal lengths and choose shortest
-        d1 = ((newp[0, t2e[2, :]] - newp[0, t2e[4, :]])**2 +
-              (newp[1, t2e[2, :]] - newp[1, t2e[4, :]])**2)
-        d2 = ((newp[0, t2e[1, :]] - newp[0, t2e[3, :]])**2 +
-              (newp[1, t2e[1, :]] - newp[1, t2e[3, :]])**2)
-        d3 = ((newp[0, t2e[0, :]] - newp[0, t2e[5, :]])**2 +
-              (newp[1, t2e[0, :]] - newp[1, t2e[5, :]])**2)
+        d1 = ((newp[0, t2e[2]] - newp[0, t2e[4]])**2 +
+              (newp[1, t2e[2]] - newp[1, t2e[4]])**2)
+        d2 = ((newp[0, t2e[1]] - newp[0, t2e[3]])**2 +
+              (newp[1, t2e[1]] - newp[1, t2e[3]])**2)
+        d3 = ((newp[0, t2e[0]] - newp[0, t2e[5]])**2 +
+              (newp[1, t2e[0]] - newp[1, t2e[5]])**2)
         I1 = d1 < d2
         I2 = d1 < d3
         I3 = d2 < d3
@@ -438,11 +448,11 @@ class MeshTet(Mesh3D):
     def shapereg(self):
         """Return the largest shape-regularity constant."""
         def edgelen(n):
-            return np.sqrt(np.sum((self.p[:, self.edges[0, self.t2e[n, :]]] -
-                                   self.p[:, self.edges[1, self.t2e[n, :]]])**2,
+            return np.sqrt(np.sum((self.p[:, self.edges[0, self.t2e[n]]] -
+                                   self.p[:, self.edges[1, self.t2e[n]]])**2,
                                   axis=0))
         edgelenmat = np.vstack(tuple(edgelen(i) for i in range(6)))
-        return np.max(np.max(edgelenmat, axis=0)/np.min(edgelenmat, axis=0))
+        return np.max(np.max(edgelenmat, axis=0) / np.min(edgelenmat, axis=0))
 
     def mapping(self):
         from skfem.mapping import MappingAffine


Test patch:
diff --git a/tests/test_elements.py b/tests/test_elements.py
index da22d6b6..8a928d17 100644
--- a/tests/test_elements.py
+++ b/tests/test_elements.py
@@ -7,63 +7,72 @@ from skfem import *
 
 
 class TestNodality(TestCase):
-    elem = ElementTriP2()
-    N = 6
+    """Test for Element.doflocs."""
+
+    elems = [
+        ElementLineP1(),
+        ElementLineP2(),
+        ElementLinePp(1),
+        ElementLinePp(3),
+        ElementTriP0(),
+        ElementTriP1(),
+        ElementTriP2(),
+        ElementTriMini(),
+        ElementQuad0(),
+        ElementQuad1(),
+        ElementQuad2(),
+        ElementQuadS2(),
+        ElementQuadP(1),
+        ElementQuadP(3),
+        ElementTetP0(),
+        ElementTetP1(),
+        ElementTetP2(),
+        ElementHex1(),
+    ]
 
     def runTest(self):
-        Ih = np.zeros((self.N, self.N))
-        for itr in range(self.N):
-            for jtr in range(self.N):
-                Ih[itr, jtr] = self.elem.lbasis(
-                    self.elem.doflocs[itr, None].T,
-                    jtr
-                )[0][0]
-        self.assertTrue(np.sum(Ih - np.eye(self.N)) < 1e-17)
-
-
-class TestQuad2Nodality(TestNodality):
-    elem = ElementQuad2()
-    N = 9
-
-
-class TestQuadS2Nodality(TestNodality):
-    elem = ElementQuadS2()
-    N = 8
-
-
-class TestTriP1Nodality(TestNodality):
-    elem = ElementTriP1()
-    N = 3
+        for e in self.elems:
+            N = e.doflocs.shape[0]
+            Ih = np.zeros((N, N))
+            for itr in range(N):
+                Ih[itr] = e.lbasis(e.doflocs.T, itr)[0]
 
+            # Remove nan-rows: test nodality only on non-nan doflocs.
+            #
+            # Some elements, such as ElementTriMini might have a combination
+            # of nodal dofs and non-nodal dofs.
+            #
+            # Nodal dof is defined so that there exists a point where the
+            # corresponding basis function is one, and other basis functions are
+            # zero. Non-nodal dof does not satisfy this property.
+            ix = np.isnan(np.sum(Ih, axis=1))
+            Nnan = np.sum(ix)
+            ixs = np.nonzero(~ix)[0]
+            Ih = Ih[ixs].T[ixs].T
 
-class TestTetP1Nodality(TestNodality):
-    elem = ElementTetP1()
-    N = 4
+            assert_array_equal(Ih, np.eye(N - Nnan),
+                               err_msg = "{}".format(type(e)))
 
 
-class TestTetP2Nodality(TestNodality):
-    elem = ElementTetP2()
-    N = 10
 
+class TestNodalityTriRT0(TestCase):
 
-class TestTetP0Nodality(TestNodality):
-    elem = ElementTetP0()
-    N = 1
+    elem = ElementTriRT0()
 
-
-class TestTri0Nodality(TestNodality):
-    elem = ElementTriP2()
-    N = 1
-
-
-class TestLineP1Nodality(TestNodality):
-    elem = ElementLineP1()
-    N = 2
-
-
-class TestLineP2Nodality(TestNodality):
-    elem = ElementLineP2()
-    N = 3
+    def runTest(self):
+        e = self.elem
+        N = e.doflocs.shape[0]
+        Ih = np.zeros((N, N))
+        normals = np.array([[0., -1.],
+                            [1 / np.sqrt(2), 1 / np.sqrt(2)],
+                            [-1., 0.]]).T
+        for itr in range(N):
+            # calculate integral of normal component over edge
+            Ih[itr] = np.sum(e.lbasis(e.doflocs.T, itr)[0] * normals, axis=0) *\
+                np.array([1., np.sqrt(2), 1.])
+
+        assert_array_equal(Ih, np.eye(N),
+                       err_msg = "{}".format(type(e)))
 
 
 class TestComposite(TestCase):