oxidize transforms

.github/workflows/test.yaml

@@ -12,24 +12,36 @@ env:
   official_container_repository: ghcr.io/evalf/nutils
 jobs:
   build-python-package:
-    name: Build Python package
-    runs-on: ubuntu-20.04
+    name: 'Build Python package on ${{ matrix.os }}'
+    runs-on: ${{ matrix.os }}
+    strategy:
+      matrix:
+        include:
+          - {os: ubuntu-latest}
+          - {os: macos-latest, build-args: --no-sdist}
+          - {os: windows-latest, build-args: --no-sdist}
     steps:
       - name: Checkout
         uses: actions/checkout@v2
-      - name: Install build dependencies
-        run: python3 -m pip install setuptools wheel
-      - name: Build package
-        run: |
-          # To make the wheels reproducible, set the timestamp of the (files in
-          # the) generated wheels to the date of the commit.
-          export SOURCE_DATE_EPOCH=`git show -s --format=%ct`
-          python3 setup.py sdist bdist_wheel
+      - uses: actions/setup-python@v2
+        if: matrix.os == 'windows-latest'
+        with:
+          python-version: '3.10'
+      - name: Install Rust toolchain
+        uses: actions-rs/toolchain@v1
+        with:
+          toolchain: stable
+          profile: minimal
+          default: true
+      - name: Build wheels
+        uses: messense/maturin-action@v1
+        with:
+          args: --release ${{ matrix.build-args }}
       - name: Upload package artifacts
         uses: actions/upload-artifact@v2
         with:
           name: python-package
-          path: dist/
+          path: target/wheels
           if-no-files-found: error
   test:
     needs: build-python-package
@@ -49,7 +61,6 @@ jobs:
           - {name: "mkl matrix parallel", os: ubuntu-latest, python-version: "3.10", matrix-backend: mkl, nprocs: 2}
           - {name: "parallel", os: ubuntu-latest, python-version: "3.10", matrix-backend: numpy, nprocs: 2}
           - {name: "numpy 1.17", os: ubuntu-latest, python-version: "3.7", matrix-backend: numpy, nprocs: 1, numpy-version: ==1.17}
-          - {name: "tensorial", os: ubuntu-latest, python-version: "3.10", matrix-backend: numpy, nprocs: 1, tensorial: test}
       fail-fast: false
     env:
       NUTILS_MATRIX: ${{ matrix.matrix-backend }}
@@ -67,7 +78,7 @@ jobs:
       - name: Move nutils directory
         run: mv nutils _nutils
       - name: Set up Python ${{ matrix.python-version }}
-        uses: actions/setup-python@v1
+        uses: actions/setup-python@v2
         with:
           python-version: ${{ matrix.python-version }}
       - name: Download Python package artifact
@@ -86,7 +97,7 @@ jobs:
           python -um pip install --upgrade --upgrade-strategy eager wheel
           python -um pip install --upgrade --upgrade-strategy eager coverage treelog stringly meshio numpy$_numpy_version
           # Install Nutils from `dist` dir created in job `build-python-package`.
-          python -um pip install --no-index --find-links ./dist nutils
+          python -um pip install --no-index --find-links ./dist --only-binary :all: nutils
       - name: Install Scipy
         if: ${{ matrix.matrix-backend == 'scipy' }}
         run: python -um pip install --upgrade --upgrade-strategy eager scipy
@@ -134,21 +145,38 @@ jobs:
           python -um pip install --upgrade --upgrade-strategy eager wheel
           python -um pip install --upgrade --upgrade-strategy eager treelog stringly matplotlib scipy pillow numpy git+https://github.com/evalf/nutils-SI.git
           # Install Nutils from `dist` dir created in job `build-python-package`.
-          python -um pip install --no-index --find-links ./dist nutils
+          python -um pip install --no-index --find-links ./dist --only-binary :all: nutils
       - name: Test
         run: python -um unittest discover -b -q -t . -s examples
   test-sphinx:
     name: Test building docs
+    needs: build-python-package
     runs-on: ubuntu-20.04
     steps:
       - name: Checkout
         uses: actions/checkout@v2
-      - name: Install dependencies
+      - name: Move nutils directory
+        run: mv nutils _nutils
+      - name: Set up Python ${{ matrix.python-version }}
+        uses: actions/setup-python@v2
+        with:
+          python-version: '3.10'
+      - name: Download Python package artifact
+        uses: actions/download-artifact@v2
+        with:
+          name: python-package
+          path: dist/
+      - name: Install Nutils and dependencies
+        id: install
+        env:
+          _numpy_version: ${{ matrix.numpy-version }}
         run: |
-          python3 -um pip install setuptools wheel
-          python3 -um pip install --upgrade --upgrade-strategy eager .[docs]
+          python -um pip install --upgrade --upgrade-strategy eager wheel
+          python -um pip install --upgrade --upgrade-strategy eager Sphinx treelog stringly meshio numpy scipy matplotlib
+          # Install Nutils from `dist` dir created in job `build-python-package`.
+          python -um pip install --no-index --find-links ./dist --only-binary :all: nutils
       - name: Build docs
-        run: python3 setup.py build_sphinx --nitpicky --warning-is-error --keep-going
+        run: python3 -m sphinx -n -W --keep-going docs build/sphinx
   build-and-test-container-image:
     name: Build container image
     needs: build-python-package

.gitignore

@@ -4,3 +4,6 @@ __pycache__/
 /dist/
 /nutils.egg-info/
 /.eggs/
+/target
+Cargo.lock
+/nutils/_rust.so

Cargo.toml

@@ -0,0 +1,18 @@
+[package]
+name = "nutils"
+version = "0.1.0"
+authors = ["Evalf <info@evalf.com>"]
+edition = "2021"
+
+[lib]
+name = "nutils"
+crate-type = ["cdylib"]
+
+[package.metadata.maturin]
+name = "nutils._rust"
+
+[dependencies]
+approx = "0.5"
+num = "0.4"
+pyo3 = { version = "0.16", features = ["extension-module", "abi3-py37"] }
+numpy = "0.16"

Makefile

@@ -0,0 +1,25 @@
+develop: build
+	cp --reflink=auto target/debug/libnutils.so nutils/_rust.so
+
+develop-release: build-release
+	cp --reflink=auto target/release/libnutils.so nutils/_rust.so
+
+build:
+	cargo build
+
+build-release:
+	cargo build --release
+
+bench:
+	cargo +nightly bench --features bench
+
+bench-python-compare:
+	python3 -m pytest benches/ --benchmark-compare --benchmark-group-by=name
+
+test-rust:
+	cargo test
+
+docs-python: develop-release
+	python3 -m sphinx -n -W --keep-going -E -D html_theme=sphinx_rtd_theme docs build/sphinx/html
+
+.PHONY: build

examples/adaptivity.py

@@ -39,20 +39,25 @@ def main(etype: str, btype: str, degree: int, nrefine: int):
 
     x, y = geom - .5
     exact = (x**2 + y**2)**(1/3) * numpy.cos(numpy.arctan2(y+x, y-x) * (2/3))
-    domain = domain.trim(exact-1e-15, maxrefine=0)
+    #domain = domain.select(exact, ischeme='gauss0')
     linreg = util.linear_regressor()
 
     for irefine in treelog.iter.fraction('level', range(nrefine+1)):
 
         if irefine:
+            basis = domain.basis(btype, degree=degree)
             refdom = domain.refined
             refbasis = refdom.basis(btype, degree=degree)
             ns.add_field('vref', refbasis)
             res = refdom.integral('∇_k(vref) ∇_k(u) dV' @ ns, degree=degree*2)
             res -= refdom.boundary.integral('vref ∇_k(u) n_k dS' @ ns, degree=degree*2)
             indicator = res.derivative('vref').eval(**args)
             supp = refbasis.get_support(indicator**2 > numpy.mean(indicator**2))
-            domain = domain.refined_by(refdom.transforms[supp])
+            if hasattr(domain, 'coord_system'):
+                domain = domain.refined_by(refdom.coord_system.trans_to(domain.coord_system).apply_indices(supp))
+            else:
+                domain = domain.refined_by(refdom.transforms[supp])
+
 
         ns = Namespace()
         ns.x = geom
@@ -61,7 +66,9 @@ def main(etype: str, btype: str, degree: int, nrefine: int):
         ns.uexact = exact
         ns.du = 'u - uexact'
 
-        sqr = domain.boundary['trimmed'].integral('u^2 dS' @ ns, degree=degree*2)
+        #sqr = domain.boundary['trimmed'].integral('u^2 dS' @ ns, degree=degree*2)
+        sqr = domain.boundary['right,bottom'].integral('u^2 dS' @ ns, degree=degree*2)
+        #sqr = domain.boundary.select(x - y, 'gauss0').integral('u^2 dS' @ ns, degree=degree*2)
         cons = solver.optimize('u,', sqr, droptol=1e-15)
 
         sqr = domain.boundary.integral('du^2 dS' @ ns, degree=7)

examples/coil.py

@@ -65,7 +65,7 @@ def main(nelems: int = 50, degree: int = 3, freq: float = 0., nturns: int = 1, r
     # domain is mapped from [0,1] to [0,inf) using an arctanh transform. Finally,
     # a Neumann boundary condition is used at z=0 to obtain symmetry in z=0.
 
-    RZ, ns.rz0 = mesh.rectilinear([numpy.linspace(0, 1, nelems)]*2, space='RZ')
+    RZ, ns.rz0 = mesh.rectilinear([numpy.linspace(0, 1, nelems)]*2, spaces=('R', 'Z'))
     REV, ns.θ = mesh.line([-numpy.pi, numpy.pi], bnames=['start', 'end'], space='Θ')
     REV0 = REV.refined[:1].boundary['end'].sample('bezier', 2)
     ns.rz = numpy.arctanh(ns.rz0) * 2 * rcoil

nutils/debug_flags.py

@@ -7,6 +7,7 @@
 sparse = _env.pop('sparse', _all or __debug__)  # check sparse chunks in evaluable
 lower = _env.pop('lower', _all or __debug__)  # check lowered shape, dtype in function
 evalf = _env.pop('evalf', _all)  # check evaluated arrays in evaluable
+hierarchical = _env.pop('hierarchical', _all)  # check hierarchical topologies
 
 if _env:
     warnings.warn('unused debug flags: {}'.format(', '.join(_env)))

nutils/element.py

@@ -8,7 +8,8 @@
 system, and provide pointsets for purposes of integration and sampling.
 """
 
-from . import util, numeric, cache, transform, warnings, types, points
+from . import util, numeric, cache, transform, warnings, types, points, _rust
+from ._rust import CoordSystem
 import numpy
 import re
 import math
@@ -353,6 +354,15 @@ def edge_transforms(self):
         assert self.ndims > 0
         return tuple(transform.SimplexEdge(self.ndims, i) for i in range(self.ndims+1))
 
+    def uniform_edges_coord_system(self, coord_system: CoordSystem, offset: int) -> CoordSystem:
+        if self.ndims == 1:
+            s = _rust.Simplex.line
+        elif self.ndims == 2:
+            s = _rust.Simplex.triangle
+        else:
+            raise NotImplementedError
+        return coord_system.edges(s, offset)
+
     @property
     def child_refs(self):
         return tuple([self] * (2**self.ndims))
@@ -361,6 +371,15 @@ def child_refs(self):
     def child_transforms(self):
         return tuple(transform.SimplexChild(self.ndims, ichild) for ichild in range(2**self.ndims))
 
+    def uniform_children_coord_system(self, coord_system: CoordSystem, offset: int) -> CoordSystem:
+        if self.ndims == 1:
+            s = _rust.Simplex.line
+        elif self.ndims == 2:
+            s = nutils._rust.Simplex.triangle
+        else:
+            raise NotImplementedError
+        return coord_system.children(s, offset)
+
     @property
     def ribbons(self):
         return tuple(((iedge1, iedge2), (iedge2+1, iedge1)) for iedge1 in range(self.ndims+1) for iedge2 in range(iedge1, self.ndims))

nutils/elementseq.py

@@ -3,6 +3,7 @@
 from . import types, numeric, util
 from .element import Reference
 from .pointsseq import PointsSequence
+from ._rust import CoordSystem
 from typing import Tuple, Sequence, Iterable, Iterator, Optional, Union, overload
 import abc
 import itertools
@@ -122,10 +123,7 @@ def __getitem__(self, index):
         if numeric.isint(index):
             return self.get(index)
         elif isinstance(index, slice):
-            index = range(len(self))[index]
-            if index == range(len(self)):
-                return self
-            return self.take(numpy.arange(index.start, index.stop, index.step))
+            return self.slice(index)
         elif numeric.isintarray(index):
             return self.take(index)
         elif numeric.isboolarray(index):
@@ -175,6 +173,25 @@ def get(self, index: int) -> Reference:
 
         raise NotImplementedError
 
+    def slice(self, __s: slice) -> 'References':
+        '''Return a slice of this sequence.
+
+        Parameters
+        ----------
+        s : :class:`slice`
+            The slice.
+
+        Returns
+        -------
+        :class:`References`
+            The slice.
+        '''
+
+        start, stop, step = __s.indices(len(self))
+        if start == 0 and stop == len(self) and step == 1:
+            return self
+        return self.take(numpy.arange(start, stop, step))
+
     def take(self, indices: numpy.ndarray) -> 'References':
         '''Return a selection of this sequence.
 
@@ -295,6 +312,9 @@ def children(self) -> 'References':
 
         return _Derived(self, 'child_refs', self.ndims)
 
+    def children_coord_system(self, coord_system: CoordSystem) -> CoordSystem:
+        raise NotImplementedError
+
     @property
     def edges(self) -> 'References':
         '''Return the sequence of edge references.
@@ -309,6 +329,9 @@ def edges(self) -> 'References':
 
         return _Derived(self, 'edge_refs', self.ndims-1)
 
+    def edges_coord_system(self, coord_system: CoordSystem) -> CoordSystem:
+        raise NotImplementedError
+
     @property
     def isuniform(self) -> 'bool':
         '''``True`` if all reference in this sequence are equal.'''
@@ -411,10 +434,16 @@ def product(self, other: References) -> References:
     def children(self) -> References:
         return References.from_iter(self.item.child_refs, self.ndims).repeat(len(self))
 
+    def children_coord_system(self, coord_system: CoordSystem) -> CoordSystem:
+        return self.item.uniform_children_coord_system(coord_system, 0)
+
     @property
     def edges(self) -> References:
         return References.from_iter(self.item.edge_refs, self.ndims-1).repeat(len(self))
 
+    def edges_coord_system(self, coord_system: CoordSystem) -> CoordSystem:
+        return self.item.uniform_edges_coord_system(coord_system, 0)
+
     @property
     def isuniform(self) -> bool:
         return True
@@ -562,10 +591,20 @@ def compress(self, mask: numpy.ndarray) -> References:
     def children(self) -> References:
         return self.sequence1.children.chain(self.sequence2.children)
 
+    def children_coord_system(self, coord_system: CoordSystem) -> CoordSystem:
+        coord_system1 = self.sequence1.children_coord_system(coord_system.slice(slice(0, len(self.sequence1))))
+        coord_system2 = self.sequence2.children_coord_system(coord_system.slice(slice(len(self.sequence1), None)))
+        return coord_system1.concat(coord_system2)
+
     @property
     def edges(self) -> References:
         return self.sequence1.edges.chain(self.sequence2.edges)
 
+    def edges_coord_system(self, coord_system: CoordSystem) -> CoordSystem:
+        coord_system1 = self.sequence1.edges_coord_system(coord_system.slice(slice(0, len(self.sequence1))))
+        coord_system2 = self.sequence2.edges_coord_system(coord_system.slice(slice(len(self.sequence1), None)))
+        return coord_system1.concat(coord_system2)
+
     def getpoints(self, ischeme: str, degree: int) -> PointsSequence:
         return self.sequence1.getpoints(ischeme, degree).chain(self.sequence2.getpoints(ischeme, degree))