Merge pull request #471 from kinnala/add-hermite-convergence-test

Add convergence test for ElementLineHermite

tests/test_convergence_h2.py

@@ -1,14 +1,15 @@
-import unittest
+from unittest import TestCase
 import numpy as np
 
 from skfem import BilinearForm, LinearForm, Functional, asm, condense, solve
 from skfem.helpers import dd, ddot
-from skfem.mesh import MeshQuad, MeshTri
-from skfem.element import ElementQuadBFS, ElementTriArgyris, ElementTriMorley
+from skfem.mesh import MeshQuad, MeshTri, MeshLine
+from skfem.element import (ElementQuadBFS, ElementTriArgyris,
+                           ElementTriMorley, ElementLineHermite)
 from skfem.assembly import InteriorBasis
 
 
-class ConvergenceMorley(unittest.TestCase):
+class ConvergenceMorley(TestCase):
 
     case = (MeshTri, ElementTriMorley)
     prerefs = 3
@@ -94,3 +95,57 @@ class ConvergenceBFS(ConvergenceMorley):
     preref = 1
     limits = (3.9, 4.5)
     abs_limit = 5e-9
+
+
+class ConvergenceHermite(TestCase):
+
+    case = (MeshLine, ElementLineHermite)
+    prerefs = 3
+    limits = (3.9, 4.1)
+    abs_limit = 8e-5
+
+    def runTest(self):
+        m = self.case[0]()
+        m.refine(self.prerefs)
+
+        hs = []
+        L2s = []
+
+        for itr in range(3):
+            e = self.case[1]()
+            ib = InteriorBasis(m, e)
+
+            @BilinearForm
+            def bilinf(u, v, w):
+                return ddot(dd(u), dd(v))
+
+            @LinearForm
+            def linf(v, w):
+                return 1. * v
+
+            K = asm(bilinf, ib)
+            f = asm(linf, ib)
+
+            x = solve(*condense(K, f, D=ib.get_dofs().all()))
+
+            X = ib.interpolate(x)
+
+            def exact(x):
+                return (x ** 2 - 2. * x ** 3 + x ** 4) / 24.
+
+            @Functional
+            def error(w):
+                return (w.w - exact(w.x)) ** 2
+
+            L2 = np.sqrt(error.assemble(ib, w=X))
+
+            L2s.append(L2)
+            hs.append(m.param())
+            m.refine()
+
+        hs = np.array(hs)
+        L2s = np.array(L2s)
+        pfit = np.polyfit(np.log10(hs), np.log10(L2s), 1)
+        self.assertGreater(pfit[0], self.limits[0])
+        self.assertLess(pfit[0], self.limits[1])
+        self.assertLess(L2s[-1], self.abs_limit)