Merge branch 'csgmsh'

.github/workflows/test.yaml

@@ -38,23 +38,23 @@ jobs:
   test:
     needs: build-python-package
     name: 'Test ${{ matrix.name }}'
-    runs-on: ${{ matrix.os }}
+    runs-on: ${{ matrix.os }}-latest
     strategy:
       matrix:
         include:
-          - {name: "baseline", os: ubuntu-latest, python-version: "3.10", matrix-backend: numpy, nprocs: 1}
-          - {name: "windows", os: windows-latest, python-version: "3.10", matrix-backend: numpy, nprocs: 1}
-          - {name: "macos", os: macos-latest, python-version: "3.10", matrix-backend: numpy, nprocs: 1}
-          - {name: "python 3.8", os: ubuntu-latest, python-version: "3.8", matrix-backend: numpy, nprocs: 1}
-          - {name: "python 3.9", os: ubuntu-latest, python-version: "3.9", matrix-backend: numpy, nprocs: 1}
-          - {name: "scipy matrix", os: ubuntu-latest, python-version: "3.10", matrix-backend: scipy, nprocs: 1}
-          - {name: "mkl linux", os: ubuntu-latest, python-version: "3.10", matrix-backend: mkl, nprocs: 1}
-          - {name: "mkl linux parallel", os: ubuntu-latest, python-version: "3.10", matrix-backend: mkl, nprocs: 2}
-          - {name: "mkl windows", os: windows-latest, python-version: "3.10", matrix-backend: mkl, nprocs: 1}
-          - {name: "mkl macos", os: macos-latest, python-version: "3.10", matrix-backend: mkl, nprocs: 1}
-          - {name: "parallel", os: ubuntu-latest, python-version: "3.10", matrix-backend: numpy, nprocs: 2}
-          - {name: "numpy 1.17", os: ubuntu-latest, python-version: "3.8", matrix-backend: numpy, nprocs: 1, numpy-version: ==1.17.3}
-          - {name: "tensorial", os: ubuntu-latest, python-version: "3.10", matrix-backend: numpy, nprocs: 1, tensorial: test}
+          - {name: "baseline", os: ubuntu, python-version: "3.10", matrix-backend: numpy, nprocs: 1}
+          - {name: "windows", os: windows, python-version: "3.10", matrix-backend: numpy, nprocs: 1}
+          - {name: "macos", os: macos, python-version: "3.10", matrix-backend: numpy, nprocs: 1}
+          - {name: "python 3.8", os: ubuntu, python-version: "3.8", matrix-backend: numpy, nprocs: 1}
+          - {name: "python 3.9", os: ubuntu, python-version: "3.9", matrix-backend: numpy, nprocs: 1}
+          - {name: "scipy matrix", os: ubuntu, python-version: "3.10", matrix-backend: scipy, nprocs: 1}
+          - {name: "mkl linux", os: ubuntu, python-version: "3.10", matrix-backend: mkl, nprocs: 1}
+          - {name: "mkl linux parallel", os: ubuntu, python-version: "3.10", matrix-backend: mkl, nprocs: 2}
+          - {name: "mkl windows", os: windows, python-version: "3.10", matrix-backend: mkl, nprocs: 1}
+          - {name: "mkl macos", os: macos, python-version: "3.10", matrix-backend: mkl, nprocs: 1}
+          - {name: "parallel", os: ubuntu, python-version: "3.10", matrix-backend: numpy, nprocs: 2}
+          - {name: "numpy 1.17", os: ubuntu, python-version: "3.8", matrix-backend: numpy, nprocs: 1, numpy-version: ==1.17.3}
+          - {name: "tensorial", os: ubuntu, python-version: "3.10", matrix-backend: numpy, nprocs: 1, tensorial: test}
       fail-fast: false
     env:
       _wheel: ${{ needs.build-python-package.outputs.wheel }}
@@ -82,8 +82,11 @@ jobs:
           name: python-package
           path: dist/
       - name: Install Graphviz
-        if: ${{ matrix.os == 'ubuntu-latest' }}
+        if: ${{ matrix.os == 'ubuntu' }}
         run: sudo apt install -y graphviz
+      - name: Install LibGLU
+        if: ${{ matrix.os == 'ubuntu' }}
+        run: sudo apt install -y libglu1
       - name: Install Nutils and dependencies
         id: install
         env:
@@ -92,7 +95,7 @@ jobs:
           python -um pip install --upgrade --upgrade-strategy eager wheel
           python -um pip install --upgrade --upgrade-strategy eager coverage numpy$_numpy_version
           # Install Nutils from `dist` dir created in job `build-python-package`.
-          python -um pip install "$_wheel[import_gmsh,export_mpl]"
+          python -um pip install "$_wheel[import_gmsh,export_mpl,cs]"
       - name: Install Scipy
         if: ${{ matrix.matrix-backend == 'scipy' }}
         run: python -um pip install --upgrade scipy
@@ -135,12 +138,15 @@ jobs:
         with:
           name: python-package
           path: dist/
+      - name: Install LibGLU
+        if: ${{ matrix.os == 'ubuntu' }}
+        run: sudo apt install -y libglu1
       - name: Install Nutils and dependencies
         id: install
         run: |
           python -um pip install --upgrade --upgrade-strategy eager wheel
           # Install Nutils from `dist` dir created in job `build-python-package`.
-          python -um pip install "$_wheel[matrix_scipy,export_mpl]"
+          python -um pip install "$_wheel[matrix_scipy,export_mpl,cs]"
       - name: Test
         run: python -um unittest discover -b -q -t . -s examples
   test-sphinx:

examples/adaptivity.py

@@ -1,10 +1,43 @@
-from nutils import mesh, function, solver, export, testing
+from nutils import mesh, function, solver, export, testing, cs
 from nutils.expression_v2 import Namespace
+from dataclasses import dataclass
+from typing import Union
 import numpy
 import treelog
 
 
-def main(etype: str = 'square',
+@dataclass
+class structured:
+
+    def generate(self):
+        topo, geom = mesh.rectilinear([2,2])
+        quadrant = topo[1:,:1].withboundary(corner='left,top')
+        return topo - quadrant, (geom - 1) / 2
+
+
+@dataclass
+class unstructured:
+
+    elemsize: float = .15
+
+    def generate(self):
+        square = cs.Rectangle(cs.Interval(-.5, .5), cs.Interval(-.5, .5))
+        quadrant = cs.Rectangle(cs.Interval(0, .5), cs.Interval(-.5, 0))
+        shapes = dict(dom=square-quadrant, corner=quadrant.boundary)
+        return mesh.csgmsh(shapes, elemsize=self.elemsize)
+
+
+@dataclass
+class multipatch:
+
+    def generate(self):
+        vertices = [0,-.5], [-.5,-.5], [0,0], [-.5,.5], [.5,0], [.5,.5]
+        patches = [0,1,2,3], [2,3,4,5]
+        topo, geom = mesh.multipatch(patchverts=vertices, patches=patches, nelems=1)
+        return topo.withboundary(corner='patch0-bottom,patch1-bottom'), geom
+
+
+def main(basemesh: Union[structured, unstructured, multipatch] = structured(),
          btype: str = 'h-std',
          degree: int = 2,
          nrefine: int = 5):
@@ -28,24 +61,22 @@ def main(etype: str = 'square',
 
     Parameters
     ----------
-    etype
-        Type of elements (square/triangle/mixed).
+    basemesh
+        Initial mesh: structured, unstructured, or multipatch.
     btype
         Type of basis function (h/th-std/spline), with availability depending
         on the configured element type.
     degree
         Polynomial degree.
     nrefine
         Number of refinement steps to perform.
+    unstructured
+        Generate triangulated domain if True, or structured domain if False.
     '''
 
-    domain, geom = mesh.unitsquare(2, etype)
-    geom -= .5 # shift domain center to origin
-
+    domain, geom = basemesh.generate()
     x, y = geom
     exact = (x**2 + y**2)**(1/3) * numpy.cos(numpy.arctan2(y+x, y-x) * (2/3))
-    selection = domain.select(exact, ischeme='gauss1')
-    domain = domain.subset(selection, newboundary='corner')
     linreg = LinearRegressor(bias=1)
 
     for irefine in treelog.iter.fraction('level', range(nrefine+1)):
@@ -135,8 +166,8 @@ def offset(self):
 
 class test(testing.TestCase):
 
-    def test_square_quadratic(self):
-        error, u = main(nrefine=2)
+    def test_structured(self):
+        error, u = main(nrefine=2, basemesh=structured())
         with self.subTest('degrees of freedom'):
             self.assertEqual(len(u), 149)
         with self.subTest('L2-error'):
@@ -151,32 +182,28 @@ def test_square_quadratic(self):
                 7AYLvMPpsqGkCTPumzWf+qV92kKevjK36ozDP/FSnh1iteWiqWuf+oMaKuyKaC1i52rKPokiF2WLA/20
                 bya+ZCPbWKRPpvgFaedebw==''')
 
-    def test_triangle_quadratic(self):
-        error, u = main(nrefine=2, etype='triangle')
+    def test_unstructured(self):
+        error, u = main(nrefine=2, basemesh=unstructured(elemsize=.5))
         with self.subTest('degrees of freedom'):
-            self.assertEqual(len(u), 98)
+            self.assertEqual(len(u), 97)
         with self.subTest('L2-error'):
-            self.assertAlmostEqual(error[0], 0.00138, places=5)
+            self.assertAlmostEqual(error[0], 0.00095, places=5)
         with self.subTest('H1-error'):
-            self.assertAlmostEqual(error[1], 0.05326, places=5)
-        with self.subTest('left-hand side'):
-            self.assertAlmostEqual64(u, '''
-                eNprMV1oesqU2VTO1Nbko6myWbhpq+kckwST90avjRgYzptYm+YYMwBBk3GQWavZb1NXs2+mm83um1WY
-                bQbyXYEiQWbKZjNM7wJVzjBlYICoPW8CMiXH+LXRR9NwoPkg82xN5IB2MZu2mGabSBnnAbGscYEJj3GV
-                YQAQg/TVGfaA7RI0BsErRjeNeowDgDQPmF9gkmciaJxtArGjzrAKCGWNpYAQAL0kOBE=''')
+            self.assertAlmostEqual(error[1], 0.04006, places=5)
 
-    def test_mixed_linear(self):
-        error, u = main(nrefine=2, etype='mixed', degree=1)
+    def test_multipatch(self):
+        error, u = main(nrefine=2, basemesh=multipatch())
         with self.subTest('degrees of freedom'):
-            self.assertEqual(len(u), 34)
+            self.assertEqual(len(u), 93)
         with self.subTest('L2-error'):
-            self.assertAlmostEqual(error[0], 0.00450, places=5)
+            self.assertAlmostEqual(error[0], 0.00128, places=5)
         with self.subTest('H1-error'):
-            self.assertAlmostEqual(error[1], 0.11692, places=5)
+            self.assertAlmostEqual(error[1], 0.05662, places=5)
         with self.subTest('left-hand side'):
             self.assertAlmostEqual64(u, '''
-                eNprMT1u6mQyxUTRzMCUAQhazL6b3jNrMYPxp5iA5FtMD+lcMgDxHa4aXzS+6HDV+fKO85cMnC8zMBzS
-                AQDBThbY''')
+                eNpjYHhnZGfCwCBgnA4ky4wfmTCZtpqamO4yXWWqbAYRgbAgKkCqQbIgVQwM5wxqjRgYTIy2GicZfzKe
+                CNQbb7LaJNF0qelbEyXTh6ZCQH2tRtOB+opNuExdTTeZzjX9Zspl5m3mYlZn1mK2GQjvw9WAzAGpAZkD
+                UgMyB6QGYg7ILpAtIBtApgMA3a820A==''')
 
 
 if __name__ == '__main__':

examples/platewithhole.py

@@ -1,4 +1,4 @@
-from nutils import mesh, function, solver, export, testing
+from nutils import mesh, function, solver, export, testing, cs
 from nutils.expression_v2 import Namespace
 from dataclasses import dataclass
 from typing import Union
@@ -42,8 +42,8 @@ def generate(self, radius):
 
 
 @dataclass
-class NURBS:
-    '''Non-Uniform Radional B-Splines
+class IGA:
+    '''Isogeometric Analysis
 
     Generate a 1x2 structured topology, map it using quadratic NURBS to a
     square domain with circular cut-out, and refine several times before
@@ -80,7 +80,33 @@ def generate(self, radius):
         return topo.withboundary(hole='left', sym='top,bottom', far='right'), geom, nurbsbasis, 5
 
 
-def main(mode: Union[FCM, NURBS] = NURBS(),
+@dataclass
+class TRI:
+    '''Triangulation
+
+    Generate a triangulated mesh with second order geometry using gmsh.
+
+    Parameters
+    ----------
+    elemsize
+        Target element size.
+    degree
+        Polynomial degree for the basis.
+    '''
+
+    elemsize: int = .1
+    degree: int = 2
+
+    def generate(self, radius):
+        rect = cs.Rectangle()
+        circ = cs.Circle(radius=radius)
+        shapes = dict(dom=rect-circ, sym=rect.boundary['left|bottom'], far=rect.boundary['right|top'], hole=circ.boundary)
+        topo, geom = mesh.csgmsh(shapes, elemsize=self.elemsize, order=2)
+        basis = topo.basis('std', self.degree)
+        return topo, geom, basis, self.degree
+
+
+def main(mode: Union[FCM, IGA, TRI] = IGA(),
          radius: float = .5,
          traction: float = .1,
          poisson: float = .3):
@@ -95,15 +121,16 @@ def main(mode: Union[FCM, NURBS] = NURBS(),
 
     The script can be run in two modes: by specifying `mode=FCM`, the circular
     hole is cut out of a regular finite element mesh by means of the Finite Cell
-    Method; by specifying `mode=NURBS` a Non-Uniform Rational BSpline geometry is
+    Method; by specifying `mode=IGA` a Non-Uniform Rational BSpline geometry is
     created to map a regular domain onto the desired shape. Either mode supports
     sub-parameters which can be specified from the command-line by attaching them
     in curly braces (e.g. `FCM{nelems=20,degree=1}`).
 
     Parameters
     ----------
     mode
-        Discretization strategy: FCM (Finite Cell Method) or NURBS.
+        Discretization strategy: FCM (Finite Cell Method), IGA (Isogeometric
+        Analysis) or TRI (triangulation).
     radius
         Cut-out radius.
     traction
@@ -136,10 +163,10 @@ def main(mode: Union[FCM, NURBS] = NURBS(),
     log.info('hole radius exact up to L2 error {:.2e}'.format(radiuserr))
 
     sqr = topo.boundary['sym'].integral('(u_i n_i)^2 dS' @ ns, degree=degree*2)
-    cons = solver.optimize(('u',), sqr, droptol=1e-15)
+    cons = solver.optimize(('u',), sqr, droptol=1e-10)
 
     sqr = topo.boundary['far'].integral('du_k du_k dS' @ ns, degree=20)
-    cons = solver.optimize(('u',), sqr, droptol=1e-15, constrain=cons)
+    cons = solver.optimize(('u',), sqr, droptol=1e-10, constrain=cons)
 
     res = topo.integral('∇_j(v_i) σ_ij dV' @ ns, degree=degree*2)
     args = solver.solve_linear('u:v', res, constrain=cons)
@@ -192,7 +219,7 @@ def test_mixed(self):
                 fyhZkYI=''')
 
     def test_nurbs0(self):
-        err, cons, args = main(mode=NURBS(nrefine=0))
+        err, cons, args = main(mode=IGA(nrefine=0))
         with self.subTest('l2-error'):
             self.assertAlmostEqual(err[0], .00200, places=5)
         with self.subTest('h1-error'):
@@ -205,7 +232,7 @@ def test_nurbs0(self):
                 eNpjYJh07qLhhnOTjb0vTDdmAAKVcy/1u85lGYforQDzFc6pGSedlzd+eP4ykA8AvkQRaA==''')
 
     def test_nurbs2(self):
-        err, cons, args = main(mode=NURBS(nrefine=2))
+        err, cons, args = main(mode=IGA(nrefine=2))
         with self.subTest('l2-error'):
             self.assertAlmostEqual(err[0], .00009, places=5)
         with self.subTest('h1-error'):
@@ -221,10 +248,17 @@ def test_nurbs2(self):
                 n2c23nZe3djqQqpx88XNxrOv7gOr0zwXZeBxztro/bmnRp7nVY1zgTjvvIXxSaBfnl3YYbzmygmgOgDU
                 Imlr''')
 
+    def test_tri(self):
+        err, cons, args = main(mode=TRI(elemsize=.5))
+        with self.subTest('l2-error'):
+            self.assertAlmostEqual(err[0], .00053, places=5)
+        with self.subTest('h1-error'):
+            self.assertAlmostEqual(err[1], .0131, places=4)
+
 
 if __name__ == '__main__':
     from nutils import cli
     cli.run(main)
 
 
-# example:tags=elasticity,FCM,NURBS
+# example:tags=elasticity,FCM,IGA

nutils/__init__.py

@@ -1,4 +1,4 @@
 'Numerical Utilities for Finite Element Analysis'
 
-__version__ = version = '9a4'
+__version__ = version = '9a5'
 version_name = 'jook-sing'