Add MeshTet.init_ball

LICENSE.md

@@ -34,7 +34,7 @@ later in this document or at the beginning of the respective files.
 Some examples under `docs/examples/` are licensed differently.  In particular,
 they have a GPL-licensed dependency and, therefore, are also GPL-licensed.
 
-- docs/examples/ex{28,32}.py
+- docs/examples/ex{28}.py
 
 > Copyright 2018-2020 scikit-fem developers
 > 

README.md

@@ -192,6 +192,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ### Unreleased
 
+- Added: `MeshTet.init_ball` for meshing a ball.
+
 ### [2.2.3] - 2020-10-16
 
 - Fixed: Remove an unnecessary dependency.

docs/examples/ex32.py

@@ -2,7 +2,7 @@
 
 .. note::
 
-   This examples requires the external package `pygmsh <https://pypi.org/project/pygmsh/>`_ and an implementation of AMG (either `pyamgcl    <https://pypi.org/project/pyamgcl>`_ or `pyamg <https://pypi.org/project/pyamg/>`_).
+   This examples requires an implementation of AMG (either `pyamgcl <https://pypi.org/project/pyamgcl>`_ or `pyamg <https://pypi.org/project/pyamg/>`_).
 
 This example again solves the Stokes problem,
 
@@ -46,27 +46,9 @@
 
 .. [McBAIN] McBain, G. D. (2016). `Creeping convection in a horizontally heated ellipsoid <http://people.eng.unimelb.edu.au/imarusic/proceedings/20/548/%20Paper.pdf>`_. *Proceedings of the Twentieth Australasian Fluid Mechanics Conference*.
 
-License
--------
-
-Copyright 2018-2020 scikit-fem developers
-
-This program is free software: you can redistribute it and/or modify
-it under the terms of the GNU General Public License as published by
-the Free Software Foundation, either version 3 of the License, or
-(at your option) any later version.
-
-This program is distributed in the hope that it will be useful,
-but WITHOUT ANY WARRANTY; without even the implied warranty of
-MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-GNU General Public License for more details.
-
-You should have received a copy of the GNU General Public License
-along with this program.  If not, see <https://www.gnu.org/licenses/>.
 
 """
-from typing import NamedTuple, Optional
-from packaging import version
+from typing import NamedTuple
 
 from skfem import *
 from skfem.io.meshio import from_meshio
@@ -77,19 +59,6 @@
 from scipy.sparse import bmat, spmatrix
 from scipy.sparse.linalg import LinearOperator, minres
 
-import pygmsh
-
-
-if version.parse(pygmsh.__version__) < version.parse('7.0.0'):
-    class NullContextManager():
-        def __enter__(self):
-            return None
-        def __exit__(self, *args):
-            pass
-    geometrycontext = NullContextManager()
-else:
-    geometrycontext = pygmsh.occ.Geometry()
-
 
 try:
     try:
@@ -109,29 +78,10 @@ def build_pc_amg(A: spmatrix, **kwargs) -> LinearOperator:
         return smoothed_aggregation_solver(A, **kwargs).aspreconditioner()
 
 
-class Ellipsoid(NamedTuple):
-
-    a: float
-    b: float
-    c: float
-
-    @property
-    def semiaxes(self) -> np.ndarray:
-        return np.array([self.a, self.b, self.c])
+class Sphere(NamedTuple):
 
-    def mesh(self, geom=None, lcar=.1) -> MeshTet:
-        with geometrycontext as g:
-            if version.parse(pygmsh.__version__) < version.parse('7.0.0'):
-                geom = pygmsh.opencascade.Geometry()
-            else:
-                geom = g
-            geom.add_ellipsoid([0.]*3, self.semiaxes, lcar)
-            if version.parse(pygmsh.__version__) < version.parse('7.0.0'):
-                return from_meshio(
-                    pygmsh.generate_mesh(geom))
-            else:
-                return from_meshio(
-                    geom.generate_mesh())
+    def mesh(self) -> MeshTet:
+        return MeshTet.init_ball()
 
     def pressure(self, x, y, z) -> np.ndarray:
         """Exact pressure at zero Grashof number.
@@ -142,7 +92,7 @@ def pressure(self, x, y, z) -> np.ndarray:
 
         """
 
-        a, b, c = self.semiaxes
+        a, b, _ = np.ones(3)  # semiaxes of ellipsoid
         return (a**2 * (3 * a**2 + b**2) * x * y
                 / (3 * a**4 + 2 * a**2 * b**2 + 3 * b**4))
 
@@ -154,8 +104,8 @@ def form(v, w):
         return LinearForm(form)
 
 
-ellipsoid = Ellipsoid(.5, .3, .2)
-mesh = ellipsoid.mesh(lcar=.1)
+ball = Sphere()
+mesh = ball.mesh()
 
 element = {'u': ElementVectorH1(ElementTetP2()),
            'p': ElementTetP1()}
@@ -198,7 +148,7 @@ def precondition(uvp):
           solver=solver_iter_krylov(minres, verbose=True, M=M)),
     [basis['u'].N])
 
-error_p = asm(ellipsoid.pressure_error(), basis['p'],
+error_p = asm(ball.pressure_error(), basis['p'],
               p=basis['p'].interpolate(pressure))
 l2error_p = np.sqrt(error_p.T @ Q @ error_p)
 

docs/listofexamples.rst

@@ -443,12 +443,12 @@ algorithm that scales to reasonably fine meshes (a million tetrahedra in a few
 minutes).
 
 .. note::
-   This examples requires the external package `pygmsh <https://pypi.org/project/pygmsh/>`__ and an implementation of algebraic multigrid (either `pyamgcl    <https://pypi.org/project/pyamgcl>`_ or `pyamg <https://pypi.org/project/pyamg/>`_).
+   This examples requires an implementation of algebraic multigrid (either `pyamgcl <https://pypi.org/project/pyamgcl>`_ or `pyamg <https://pypi.org/project/pyamg/>`_).
 
-.. figure:: https://user-images.githubusercontent.com/973268/87859195-fcd08980-c93b-11ea-930e-ddcd26aabdb4.png
+.. figure:: https://user-images.githubusercontent.com/1588947/96520786-8a18d680-12bb-11eb-981a-c3388f2c8e35.png
 
-   The pressure field of Example 32.
-   The figure was created using `ParaView <https://www.paraview.org/>`_.
+   The velocity and pressure fields of Example 32, clipped in the plane of spanwise symmetry, *z* = 0.
+   The figure was created using `ParaView <https://www.paraview.org/>`_ 5.8.1.
 
 See the `source code of Example 32 <https://github.com/kinnala/scikit-fem/blob/master/docs/examples/ex32.py>`_ for more information.
 

skfem/mesh/mesh3d/mesh_tet.py

@@ -173,6 +173,39 @@ def init_tensor(cls: Type,
 
         return cls(p, T.astype(np.int64))
 
+    @classmethod
+    def init_ball(cls: Type,
+                  Nrefs: int = 3) -> Mesh3D:
+        r"""Initialize a ball mesh.
+
+        Parameters
+        ----------
+        Nrefs
+            Number of refinements, by default 3.
+
+        """
+        p = np.array([[0., 0., 0.],
+                      [1., 0., 0.],
+                      [0., 1., 0.],
+                      [0., 0., 1.],
+                      [-1., 0., 0.],
+                      [0., -1., 0.],
+                      [0., 0., -1.]]).T
+        t = np.array([[0, 1, 2, 3],
+                      [0, 4, 5, 6],
+                      [0, 1, 2, 6],
+                      [0, 1, 3, 5],
+                      [0, 2, 3, 4],
+                      [0, 4, 5, 3],
+                      [0, 4, 6, 2],
+                      [0, 5, 6, 1]], dtype=np.intp).T
+        m = cls(p, t)
+        for _ in range(Nrefs):
+            m.refine()
+            D = m.boundary_nodes()
+            m.p[:, D] = m.p[:, D] / np.linalg.norm(m.p[:, D], axis=0)
+        return m
+
     def _build_mappings(self):
         """Build element-to-facet, element-to-edges, etc. mappings."""
         # define edges: in the order (0,1) (1,2) (0,2) (0,3) (1,3) (2,3)

tests/test_examples.py

@@ -251,7 +251,7 @@ class TestEx32(TestCase):
 
     def runTest(self):
         from docs.examples.ex32 import l2error_p
-        self.assertLess(l2error_p, 3e-7)
+        self.assertLess(l2error_p, 1e-5)
 
 
 class TestEx33(TestCase):