Add simple support for quadratic refine (#796)

skfem/mesh/mesh_2d_2.py

@@ -1,22 +1,6 @@
 class Mesh2D2:
     """Mixin for quadratic 2D meshes."""
 
-    def plot(self, *args, **kwargs):
-        """Convenience wrapper for :func:`skfem.visuals.matplotlib.plot`."""
-        from skfem.visuals.matplotlib import plot, show
-        from skfem.assembly import CellBasis
-        ax = plot(CellBasis(self, self.elem()), *args, **kwargs)
-        ax.show = show
-        return ax
-
-    def draw(self, *args, **kwargs):
-        """Convenience wrapper for :func:`skfem.visuals.matplotlib.draw`."""
-        from skfem.visuals.matplotlib import draw, show
-        from skfem.assembly import CellBasis
-        ax = draw(CellBasis(self, self.elem()), *args, **kwargs)
-        ax.show = show
-        return ax
-
     def _repr_svg_(self) -> str:
         from skfem.visuals.svg import draw
         return draw(self, nrefs=2, boundaries_only=True).svg

skfem/mesh/mesh_hex_2.py

@@ -39,3 +39,6 @@ class MeshHex2(MeshHex1):
                                  [1., .5, .5],
                                  [.5, .5, .5]], dtype=np.float64).T
     elem: Type[Element] = ElementHex2
+
+    def _uniform(self):
+        return MeshHex2.from_mesh(MeshHex1.from_mesh(self).refined())

skfem/mesh/mesh_quad_2.py

@@ -23,3 +23,6 @@ class MeshQuad2(Mesh2D2, MeshQuad1):
                                  [.5, 1.],
                                  [.5, .5]], dtype=np.float64).T
     elem: Type[Element] = ElementQuad2
+
+    def _uniform(self):
+        return MeshQuad2.from_mesh(MeshQuad1.from_mesh(self).refined())

skfem/mesh/mesh_tet_2.py

@@ -39,6 +39,7 @@ class MeshTet2(MeshTet1):
                                  [1., .5, 1.],
                                  [1., 1., .5]], dtype=np.float64).T
     elem: Type[Element] = ElementTetP2
+    affine: bool = False
 
     @classmethod
     def init_ball(cls: Type, nrefs: int = 3) -> 'MeshTet2':
@@ -48,3 +49,9 @@ def init_ball(cls: Type, nrefs: int = 3) -> 'MeshTet2':
         doflocs = M.doflocs.copy()
         doflocs[:, D] /= np.linalg.norm(doflocs[:, D], axis=0)
         return replace(M, doflocs=doflocs)
+
+    def _uniform(self):
+        return MeshTet2.from_mesh(MeshTet1.from_mesh(self).refined())
+
+    def _adaptive(self, marked):
+        return MeshTet2.from_mesh(MeshTet1.from_mesh(self).refined(marked))

skfem/mesh/mesh_tri_2.py

@@ -35,3 +35,9 @@ def init_circle(cls: Type,
         doflocs = M.doflocs.copy()
         doflocs[:, D] /= np.linalg.norm(doflocs[:, D], axis=0)
         return replace(M, doflocs=doflocs)
+
+    def _uniform(self):
+        return MeshTri2.from_mesh(MeshTri1.from_mesh(self).refined())
+
+    def _adaptive(self, marked):
+        return MeshTri2.from_mesh(MeshTri1.from_mesh(self).refined(marked))

tests/test_convergence.py

@@ -31,7 +31,7 @@
 )
 from skfem.assembly import FacetBasis, Basis
 from skfem.mesh import (MeshHex, MeshLine, MeshQuad, MeshTet, MeshTri,
-                        MeshWedge1)
+                        MeshWedge1, MeshQuad2, MeshTri2, MeshTet2, MeshHex2)
 
 
 class ConvergenceQ1(unittest.TestCase):
@@ -170,6 +170,12 @@ def setUp(self):
         self.mesh = MeshQuad().refined(2)
 
 
+class ConvergenceQ1QuadraticMesh(ConvergenceQ1):
+
+    def setUp(self):
+        self.mesh = MeshQuad2().refined(2)
+
+
 class ConvergenceQ2(ConvergenceQ1):
 
     rateL2 = 3.0
@@ -200,6 +206,12 @@ def setUp(self):
         self.mesh = MeshTri.init_sqsymmetric().refined(2)
 
 
+class ConvergenceTriP1QuadraticMesh(ConvergenceTriP1):
+
+    def setUp(self):
+        self.mesh = MeshTri2.from_mesh(MeshTri.init_sqsymmetric().refined(2))
+
+
 class ConvergenceTriP2(ConvergenceTriP1):
 
     rateL2 = 3.0
@@ -311,6 +323,12 @@ def setUp(self):
         self.mesh = MeshHex().refined(2)
 
 
+class ConvergenceHex1QuadraticMesh(ConvergenceHex1):
+
+    def setUp(self):
+        self.mesh = MeshHex2().refined(1)
+
+
 class ConvergenceHexSplitTet1(ConvergenceQ1):
 
     rateL2 = 2.0

tests/test_manufactured.py

@@ -295,7 +295,7 @@ class SolveCirclePoissonTet2(SolveCirclePoisson):
     mesh_type = MeshTet2
     element_type = ElementTetP2
     filename = "quadratic_sphere_tet.msh"
-    maxval = 0.0405901240018571
+    maxval = 0.041655619826791175
 
 
 @pytest.mark.parametrize(