Maintenance (#1082)

* Fix MeshDG.draw

* Deprecate skfem.visuals.glvis

* skip ex32 test for python 3.12 due to pyamg assuming setuptools

* Fix FacetBasis for quadratic meshes

* Add Mesh.remove_unused_nodes

* add Mesh.remove_duplicate_nodes

* reorient boundary normals to outward normals in ElementGlobal

* colorbar=False will now disable colorbar

Author: kinnala

README.md

@@ -214,6 +214,8 @@ with respect to documented and/or tested features.
 - Removed: Python 3.7 support
 - Removed: `MappingMortar` and `MortarFacetBasis` in favor of
   `skfem.supermeshing`
+- Deprecated: `skfem.visuals.glvis`; current version is broken and no
+  replacement is being planned
 - Added: Python 3.12 support
 - Added: `Mesh.load` supports new keyword arguments
   `ignore_orientation=True` and `ignore_interior_facets=True` which
@@ -222,7 +224,18 @@ with respect to documented and/or tested features.
   respectively.
 - Added: `skfem.supermeshing` (requires `shapely>=2`) for creating quadrature
   rules for interpolating between two 1D or 2D meshes.
+- Added: `Mesh.remove_unused_nodes`
+- Added: `Mesh.remove_duplicate_nodes`
+- Added: `Mesh.remove_elements` now supports passing any subdomain
+  reference through `Mesh.normalize_elements`; subdomains and
+  boundaries are also properly preserved
 - Fixed: `MeshTet` uniform refine was reindexing subdomains incorrectly
+- Fixed: `MeshDG.draw` did not work; now calls `Basis.draw` which
+  works for any mesh topology
+- Fixed: `FacetBasis` now works with `MeshTri2`, `MeshQuad2`,
+  `MeshTet2` and `MeshHex2`
+- Fixed: `ElementGlobal` now uses outward normals to initialize DOFs
+  on boundary factes
 
 ## [8.1.0] - 2023-06-16
 

docs/examples/ex35.py

@@ -221,10 +221,9 @@
 from packaging import version
 from pathlib import Path
 
-from skfem.mesh import MeshTri
-from skfem.assembly import Basis, FacetBasis
-from skfem.utils import solve, asm, condense, projection
-from skfem.element import ElementTriP1
+from skfem import (MeshTri, Basis, FacetBasis,
+                   solve, asm, condense, projection,
+                   ElementTriP1)
 from skfem.models.poisson import laplace, unit_load, mass
 from skfem.io.json import from_file
 

skfem/element/element_global.py

@@ -166,6 +166,20 @@ def _eval_dofs(self, mesh, tind=None):
                                         -w['n'][itr, 0, :]])
                 w['n'][itr] /= np.linalg.norm(w['n'][itr], axis=0)
 
+            # swap
+            from skfem.assembly import FacetBasis
+            fb = FacetBasis(mesh, mesh.elem(), intorder=0)
+            ix = np.isin(mesh.t2f, mesh.boundary_facets())
+            sortix = np.argsort(mesh.t2f[ix])
+            norms1 = np.vstack((w['n'][:, 0, :][ix][sortix],
+                                w['n'][:, 1, :][ix][sortix]))
+            norms2 = fb.normals[:, :, 0]
+            signs = np.diag(norms1.T @ norms2)
+            signs1 = signs.copy()
+            signs1[sortix] = signs
+            w['n'][:, 0, :][ix] *= signs1
+            w['n'][:, 1, :][ix] *= signs1
+
         # evaluate dofs, gdof implemented in subclasses
         for itr in range(N):
             for jtr in range(N):

skfem/mapping/mapping_isoparametric.py

@@ -66,6 +66,13 @@ def Fmap(self, i, X, tind=None):
     def bndmap(self, i, X, find=None):
         p = self.mesh.doflocs
         facets = self.mesh.facets
+        if len(self.mesh.dofs.edge_dofs) > 0:
+            facets = np.vstack((facets,
+                                self.mesh.dofs.edge_dofs[0, self.mesh.f2e]))
+            # TODO currently supports only one DOF per edge (slice 0 idx)
+        if len(self.mesh.dofs.facet_dofs) > 0:
+            facets = np.vstack((facets,
+                                self.mesh.dofs.facet_dofs))
         if find is None:
             out = np.zeros((facets.shape[1], X.shape[1]))
             for itr in range(facets.shape[0]):
@@ -82,6 +89,13 @@ def bndmap(self, i, X, find=None):
     def bndJ(self, i, j, X, find=None):
         p = self.mesh.doflocs
         facets = self.mesh.facets
+        if len(self.mesh.dofs.edge_dofs) > 0:
+            facets = np.vstack((facets,
+                                self.mesh.dofs.edge_dofs[0, self.mesh.f2e]))
+            # TODO currently supports only one DOF per edge (slice 0 idx)
+        if len(self.mesh.dofs.facet_dofs) > 0:
+            facets = np.vstack((facets,
+                                self.mesh.dofs.facet_dofs))
         if find is None:
             out = np.zeros((facets.shape[1], X.shape[1]))
             for itr in range(facets.shape[0]):

skfem/mesh/mesh.py

@@ -112,6 +112,15 @@ def f2t(self):
             self._f2t = self.build_inverse(self.t, self.t2f)
         return self._f2t
 
+    @property
+    def f2e(self):
+        if not hasattr(self, '_f2e'):
+            _, self._f2e = self.build_entities(
+                self.facets,
+                self.bndelem.refdom.facets,
+            )
+        return self._f2e
+
     @property
     def edges(self):
         if not hasattr(self, '_edges'):
@@ -573,6 +582,13 @@ def is_valid(self, raise_=False) -> bool:
         logger.debug("Mesh validation completed with no warnings.")
         return True
 
+    @staticmethod
+    def _remove_duplicate_nodes(p, t):
+        tmp = np.ascontiguousarray(p.T)
+        tmp, ixa, ixb = np.unique(tmp.view([('', tmp.dtype)] * tmp.shape[1]),
+                                  return_index=True, return_inverse=True)
+        return p[:, ixa], ixb[t]
+
     def __iter__(self):
         return iter((self.doflocs, self.t))
 
@@ -606,12 +622,7 @@ def __add__(self, other):
         p = np.hstack((self.p.round(decimals=8),
                        other.p.round(decimals=8)))
         t = np.hstack((self.t, other.t + self.p.shape[1]))
-        tmp = np.ascontiguousarray(p.T)
-        tmp, ixa, ixb = np.unique(tmp.view([('', tmp.dtype)] * tmp.shape[1]),
-                                  return_index=True, return_inverse=True)
-        p = p[:, ixa]
-        t = ixb[t]
-        return cls(p, t)
+        return cls(*self._remove_duplicate_nodes(p, t))
 
     def __repr__(self):
         rep = ""
@@ -1038,27 +1049,42 @@ def restrict(self, elements):
                                         newf[self.boundaries[k]])
                           for k in self.boundaries}
                          if self.boundaries is not None else None),
-            _subdomains=({k: newt[np.intersect1d(self.subdomains[k], elements)]
+            _subdomains=({k: newt[np.intersect1d(self.subdomains[k],
+                                                 elements).astype(np.int64)]
                           for k in self.subdomains}
                          if self.subdomains is not None else None),
         )
 
-    def remove_elements(self, element_indices: ndarray):
+    def remove_elements(self, elements):
         """Construct a new mesh by removing elements.
 
         Parameters
         ----------
-        element_indices
-            List of element indices to remove.
+        elements
+            Criteria of which elements to include.  This input is normalized
+            using ``self.normalize_elements``.
 
         """
-        p, t, _ = self._reix(np.delete(self.t, element_indices, axis=1))
+        elements = self.normalize_elements(elements)
+        return self.restrict(np.setdiff1d(np.arange(self.t.shape[1],
+                                                    dtype=np.int64),
+                                          elements))
+
+    def remove_unused_nodes(self):
+        p, t, _ = self._reix(self.t)
+        return replace(
+            self,
+            doflocs=p,
+            t=t,
+        )
+
+    def remove_duplicate_nodes(self):
+        p, t = self._remove_duplicate_nodes(self.doflocs,
+                                            self.t)
         return replace(
             self,
             doflocs=p,
             t=t,
-            _boundaries=None,
-            _subdomains=None,
         )
 
     def element_finder(self, mapping=None):

skfem/mesh/mesh_2d_2.py

@@ -11,3 +11,7 @@ def element_finder(self, *args, **kwargs):
     @classmethod
     def init_refdom(cls):
         return cls.__bases__[-1].init_refdom()
+
+    def draw(self, *args, **kwargs):
+        from ..assembly import CellBasis
+        return CellBasis(self, self.elem()).draw(*args, **kwargs)

skfem/mesh/mesh_dg.py

@@ -109,3 +109,7 @@ def load(cls, *args, **kwargs):
 
     def element_finder(self, *args, **kwargs):
         raise NotImplementedError
+
+    def draw(self, *args, **kwargs):
+        from ..assembly import CellBasis
+        return CellBasis(self, self.elem()).draw(*args, **kwargs)

skfem/mesh/mesh_tri_1_dg.py

@@ -16,16 +16,25 @@ class MeshTri1DG(MeshDG, MeshTri1):
     doflocs: ndarray = field(
         default_factory=lambda: np.array(
             [
-                [0.0, 0.0],
-                [1.0, 0.0],
-                [0.0, 1.0],
-                [1.0, 0.0],
-                [0.0, 1.0],
-                [1.0, 1.0],
+                [0., 0.],
+                [1., 0.],
+                [1., 1.],
+                [0., 0.],
+                [0., 1.],
+                [1., 1.],
             ],
             dtype=np.float64,
         ).T
     )
+    t: ndarray = field(
+        default_factory=lambda: np.array(
+            [
+                [0, 1, 3],
+                [0, 2, 3],
+            ],
+            dtype=np.int64,
+        ).T
+    )
     elem: Type[Element] = ElementTriP1DG
     affine: bool = False
     sort_t: bool = False

skfem/refdom.py

@@ -104,6 +104,10 @@ class RefTet(Refdom):
              [0, 3],
              [1, 3],
              [2, 3]]
+    normals = np.array([[0., 0., -1.],
+                        [0., -1., 0.],
+                        [-1., 0., 0.],
+                        [1., 1., 1.]])
     brefdom = RefTri
     nnodes = 4
     nfacets = 4

skfem/utils.py

@@ -717,29 +717,6 @@ def projection(fun,
                diff: Optional[int] = None,
                I: Optional[ndarray] = None,
                expand: bool = False) -> ndarray:
-    """Perform projections onto a finite element basis.
-
-    Parameters
-    ----------
-    fun
-        A solution vector or a function handle.
-    basis_to
-        The finite element basis to project to.
-    basis_from
-        The finite element basis to project from.
-    diff
-        Differentiate with respect to the given dimension.
-    I
-        Index set for limiting the projection to a subset.
-    expand
-        Passed to :func:`skfem.utils.condense`.
-
-    Returns
-    -------
-    ndarray
-        The projected solution vector.
-
-    """
 
     @BilinearForm
     def mass(u, v, w):
@@ -783,7 +760,7 @@ def deriv(u, v, w):
     return solve_linear(M, f)
 
 
-@deprecated("Basis.project")
+@deprecated("Basis.project (will be removed in the next release)")
 def project(fun,
             basis_from: Optional[AbstractBasis] = None,
             basis_to: Optional[AbstractBasis] = None,

skfem/visuals/glvis.py

@@ -5,6 +5,8 @@
 from typing import Optional
 from glvis import glvis
 
+from skfem.generic_utils import deprecated
+
 
 MESH_TYPE_MAPPING = {
     skfem.MeshTet1: 4,
@@ -64,6 +66,7 @@ def _to_float_string(arr):
     return s.getvalue().decode()
 
 
+@deprecated("skfem.visuals.matplotlib.plot (no replacement)")
 def plot(basis, x, keys: Optional[str] = None):
     m = basis.mesh
     if isinstance(basis.elem, skfem.ElementVector):
@@ -99,5 +102,6 @@ def plot(basis, x, keys: Optional[str] = None):
     ))
 
 
+@deprecated("skfem.visuals.matplotlib.draw (no replacement)")
 def draw(basis, keys: Optional[str] = None):
     return plot(basis, basis.zeros(), keys=keys)

skfem/visuals/matplotlib.py

@@ -291,7 +291,7 @@ def plot_meshtri(m: MeshTri1, z: ndarray, **kwargs) -> Axes:
                       levels=kwargs["levels"],
                       **{**{'colors': 'k'}, **plot_kwargs})
 
-    if "colorbar" in kwargs:
+    if "colorbar" in kwargs and kwargs["colorbar"] is not False:
         if isinstance(kwargs["colorbar"], str):
             plt.colorbar(im, ax=ax, label=kwargs["colorbar"])
         elif isinstance(kwargs["colorbar"], dict):

tests/test_assembly.py

@@ -747,3 +747,24 @@ def test_hcurl_projections_3d():
 
     assert_almost_equal(curly, curlx)
     assert_almost_equal(curly, curla)
+
+
+def test_element_global_boundary_normal():
+
+    mesh = MeshTri.init_sqsymmetric().refined(3)
+    basis = Basis(mesh, ElementTriMorley())
+
+    D = basis.get_dofs().all('u_n')
+    x = basis.zeros()
+    x[D] = 1
+
+    @BilinearForm
+    def bilinf(u, v, w):
+        return ddot(u.hess, v.hess)
+
+    A = bilinf.assemble(basis)
+    f = 0 * unit_load.assemble(basis)
+
+    y = solve(*condense(A, f, D=basis.get_dofs(), x=x))
+
+    assert (y[basis.get_dofs(elements=True).all('u')] <= 0).all()

tests/test_examples.py

@@ -1,4 +1,6 @@
 from unittest import TestCase, main
+import pytest
+import sys
 
 import numpy as np
 
@@ -14,7 +16,8 @@ class TestEx02(TestCase):
 
     def runTest(self):
         import docs.examples.ex02 as ex02
-        self.assertAlmostEqual(np.max(ex02.x), 0.001217973811129439)
+        self.assertAlmostEqual(np.max(ex02.x[ex02.ib.nodal_dofs[0]]),
+                               0.00033840961095522285)
 
 
 class TestEx03(TestCase):
@@ -254,6 +257,8 @@ def runTest(self):
         self.assertAlmostEqual(L[0], 22.597202568397734, delta=1e-6)
 
 
+@pytest.mark.skipif(sys.version_info > (3,12),
+                    reason="Python 3.12 has no setuptools; assumed by pyamg")
 class TestEx32(TestCase):
 
     def runTest(self):