diff --git a/LICENSE.md b/LICENSE.md
index 37e9356f9..103c0a8bb 100644
--- a/LICENSE.md
+++ b/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
>
diff --git a/README.md b/README.md
index ac2b12c9c..71da5f89d 100644
--- a/README.md
+++ b/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.
diff --git a/docs/examples/ex32.py b/docs/examples/ex32.py
index 36b155864..ac017428a 100644
--- a/docs/examples/ex32.py
+++ b/docs/examples/ex32.py
@@ -2,7 +2,7 @@
.. note::
- This examples requires the external package `pygmsh `_ and an implementation of AMG (either `pyamgcl `_ or `pyamg `_).
+ This examples requires an implementation of AMG (either `pyamgcl `_ or `pyamg `_).
This example again solves the Stokes problem,
@@ -46,27 +46,9 @@
.. [McBAIN] McBain, G. D. (2016). `Creeping convection in a horizontally heated ellipsoid `_. *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 .
"""
-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)
diff --git a/docs/listofexamples.rst b/docs/listofexamples.rst
index dd1f50253..d6b329253 100644
--- a/docs/listofexamples.rst
+++ b/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 `__ and an implementation of algebraic multigrid (either `pyamgcl `_ or `pyamg `_).
+ This examples requires an implementation of algebraic multigrid (either `pyamgcl `_ or `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 `_.
+ The velocity and pressure fields of Example 32, clipped in the plane of spanwise symmetry, *z* = 0.
+ The figure was created using `ParaView `_ 5.8.1.
See the `source code of Example 32 `_ for more information.
diff --git a/skfem/mesh/mesh3d/mesh_tet.py b/skfem/mesh/mesh3d/mesh_tet.py
index d920a20c1..1735d8b26 100644
--- a/skfem/mesh/mesh3d/mesh_tet.py
+++ b/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)
diff --git a/tests/test_examples.py b/tests/test_examples.py
index 561cc589b..9d6bdde81 100644
--- a/tests/test_examples.py
+++ b/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):