Merge pull request #267 from gdmcbain/q0

Q0

docs/examples/ex25.py

@@ -1,5 +1,5 @@
 from skfem import *
-from skfem.models.poisson import laplace
+from skfem.models.poisson import laplace, mass
 
 from math import ceil
 
@@ -10,9 +10,10 @@
 length = 10.
 peclet = 1e2
 
-mesh = (MeshLine(np.linspace(0, length, ceil(mesh_inlet_n / height * length)))
-        * MeshLine(np.linspace(0, height / 2, mesh_inlet_n)))._splitquads()
-basis = InteriorBasis(mesh, ElementTriP2())
+mesh = MeshQuad.init_tensor(
+    np.linspace(0, length, ceil(mesh_inlet_n / height * length)),
+    np.linspace(0, height / 2, mesh_inlet_n))
+basis = InteriorBasis(mesh, ElementQuad2())
 
 
 @bilinear_form
@@ -31,11 +32,15 @@ def advection(u, du, v, dv, w):
 t[dofs['floor'].all()] = 1.
 t = solve(*condense(A, np.zeros_like(t), t, I=interior))
 
+basis0 = InteriorBasis(mesh, ElementQuad0(), intorder=basis.intorder)
+t0 = solve(asm(mass, basis0),
+           asm(mass, basis, basis0) @ t)
+
 
 if __name__ == '__main__':
 
     from pathlib import Path
 
-    mesh.plot(t[basis.nodal_dofs.flatten()], edgecolors='none')
+    mesh.plot(t0, edgecolors='none')
     mesh.savefig(Path(__file__).with_suffix('.png'),
                  bbox_inches='tight', pad_inches=0)

skfem/__init__.py

@@ -42,6 +42,7 @@
            'ElementHdiv',
            'ElementHex1',
            'ElementTriMorley',
+           'ElementQuad0',
            'ElementQuad1',
            'ElementQuad2',
            'ElementTetN0',

skfem/element/__init__.py

@@ -21,7 +21,7 @@
 from .element_tri import ElementTriP1, ElementTriP2, ElementTriDG,\
     ElementTriP0, ElementTriRT0, ElementTriMorley,\
     ElementTriArgyris
-from .element_quad import ElementQuad1, ElementQuad2
+from .element_quad import ElementQuad0, ElementQuad1, ElementQuad2
 from .element_tet import ElementTetP0, ElementTetP1, ElementTetP2,\
     ElementTetRT0, ElementTetN0
 from .element_hex import ElementHex1

skfem/element/element_quad/__init__.py

@@ -1,2 +1,3 @@
+from .element_quad0 import ElementQuad0
 from .element_quad1 import ElementQuad1
 from .element_quad2 import ElementQuad2

skfem/element/element_quad/element_quad0.py

@@ -0,0 +1,13 @@
+import numpy as np
+from ..element_h1 import ElementH1
+
+
+class ElementQuad0(ElementH1):
+    interior_dofs = 1
+    dim = 2
+    maxdeg = 0
+    dofnames = ['u']
+    doflocs = np.zeros((1, 2))
+
+    def lbasis(self, X, i):
+        return np.ones(X.shape[1]), np.zeros_like(X)

tests/test_examples.py

@@ -187,7 +187,9 @@ def runTest(self):
 class TestEx25(unittest.TestCase):
     def runTest(self):
         import docs.examples.ex25 as ex
-        self.assertAlmostEqual(np.mean(ex.t), 0.4642600944590631, places=5)
+        mu = np.mean(ex.t)
+        self.assertAlmostEqual(mu, 0.4642600944590631, places=5)
+        self.assertAlmostEqual(np.mean(ex.t0), mu, places=2)
 
 
 class TestEx26(unittest.TestCase):