diff --git a/.github/workflows/test.yaml b/.github/workflows/test.yaml index a8175fa97..003423841 100644 --- a/.github/workflows/test.yaml +++ b/.github/workflows/test.yaml @@ -134,7 +134,7 @@ jobs: id: install run: | python -um pip install --upgrade wheel - python -um pip install --upgrade treelog stringly matplotlib scipy pillow numpy + python -um pip install --upgrade 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 - name: Test diff --git a/devtools/container/build.py b/devtools/container/build.py index 04a9c62f7..49a58902f 100644 --- a/devtools/container/build.py +++ b/devtools/container/build.py @@ -57,9 +57,10 @@ if not examples.exists(): raise SystemExit(f'cannot find example dir: {examples}') - container = stack.enter_context(Container.new_from(base, network='none', mounts=[Mount(src=dist, dst='/mnt')])) + container = stack.enter_context(Container.new_from(base, mounts=[Mount(src=dist, dst='/mnt')])) container.run('pip', 'install', '--no-cache-dir', '--no-index', '--find-links', 'file:///mnt/', 'nutils', env=dict(PYTHONHASHSEED='0')) + container.run('pip', 'install', '--no-cache-dir', 'https://github.com/evalf/nutils-SI/archive/main.tar.gz', env=dict(PYTHONHASHSEED='0')) container.add_label('org.opencontainers.image.url', 'https://github.com/evalf/nutils') container.add_label('org.opencontainers.image.source', 'https://github.com/evalf/nutils') container.add_label('org.opencontainers.image.authors', 'Evalf') diff --git a/examples/cahnhilliard.py b/examples/cahnhilliard.py index e653cf5fa..58ad3ac4b 100644 --- a/examples/cahnhilliard.py +++ b/examples/cahnhilliard.py @@ -1,101 +1,203 @@ #! /usr/bin/env python3 # -# In this script we solve the Cahn-Hiilliard equation, which models the -# unmixing of two phases under the effect of surface tension. +# In this script we solve the Cahn-Hilliard equation, which models the unmixing +# of two phases (φ=+1 and φ=-1) under the influence of surface tension. It is a +# mixed equation of two scalar equations for phase φ and chemical potential η:: +# +# dφ/dt = -div(J(η)) +# η = ψ'(φ) σ / ε - σ ε Δ(φ) +# +# along with constitutive relations for the flux vector J = -M ∇η and the +# double well potential ψ = .25 (φ² - 1)², and subject to boundary conditions +# ∂ₙφ = -σd / σ ε and ∂ₙη = 0. Parameters are the interface thickness ε, fluid +# surface tension σ, differential wall surface tension σd, and mobility M. +# +# Cahn-Hilliard is a diffuse interface model, which means that phases do not +# separate sharply, but instead form a transition zone between the phases. The +# transition zone has a thickness proportional to ε, as is readily confirmed in +# one dimension, where a steady state solution on an infinite domain is formed +# by η(x) = 0, φ(x) = tanh(x / √2 ε). +# +# The main driver of the unmixing process is the double well potential ψ that +# is proportional to the mixing energy, taking its minima at the pure phases +# φ=+1 and φ=-1. The interface between the phases adds an energy contribution +# proportional to its length. At the wall we have a phase-dependent fluid-solid +# energy. Over time, the system minimizes the total energy:: +# +# E(φ) = ∫_Ω ψ(φ) σ / ε + ∫_Ω .5 σ ε ‖∇φ‖² + ∫_Γ (σm + φ σd) +# \ \ \ +# mixing energy interface energy wall energy +# +# Proof: the time derivative of E followed by substitution of the strong form +# and boundary conditions yields dE/dt = ∫_Ω η dφ/dt = -∫_Ω M ‖∇η‖² ≤ 0. □ +# +# Switching to discrete time we set dφ/dt = (φ - φ0) / dt and add a stabilizing +# perturbation term δψ(φ, φ0) to the doube well potential for reasons outlined +# below. This yields the following time discrete system:: +# +# φ = φ0 - dt div(J(η)) +# η = (ψ'(φ) + δψ'(φ, ψ0)) σ / ε - σ ε Δ(φ) +# +# with the equivalent weak formulation:: +# +# ∂/∂η ∫_Ω [ η (φ - φ0) + .5 dt J(η)·∇η ] = 0 +# ∂/∂φ ∫_Ω [ E(φ) + δψ(φ, φ0) σ / ε - η φ ] = 0 +# +# For stability we wish for the perturbation δψ to be such that the time +# discrete system preserves the energy dissipation property E(φ) ≤ E(φ0) for +# any timestep dt. To derive suitable perturbation terms to this effect, we +# define without loss of generality δψ'(φ, φ0) = .5 (φ - φ0) f(φ, φ0) and +# derive the following condition for unconditional stability:: +# +# E(φ) - E(φ0) = ∫_Ω .5 (1 - φ² - .5 (φ + φ0)² - f(φ, φ0)) (φ - φ0)² σ / ε +# - ∫_Ω (.5 σ ε ‖∇φ - ∇φ0‖² + dt M ‖∇η‖²) ≤ 0 +# +# The inequality holds true if the perturbation f is bounded from below such +# that f(φ, φ0) ≥ 1 - φ² - .5 (φ + φ0)². To keep the energy minima at the pure +# phases we additionally impose that f(±1, φ0) = 0, and select 1 - φ² as a +# suitable upper bound which we will call "nonlinear". +# +# We next observe that η is linear in φ if f(φ, φ0) = g(φ0) - φ² - (φ + φ0)² +# for any function g, which dominates if g(φ0) ≥ 1 + .5 (φ + φ0)². While this +# cannot be made to hold true for all φ, we make it hold for -√2 ≤ φ, φ0 ≤ √2 +# by defining g(φ0) = 2 + 2 ‖φ0‖ + φ0², which we will call "linear". This +# scheme further satisfies a weak minima preservation f(±1, ±‖φ0‖) = 0. +# +# We have thus arrived at the three stabilization schemes implemented here: +# +# - nonlinear: f(φ, φ0) = 1 - φ² +# - linear: f(φ, φ0) = 2 + 2 ‖φ0‖ - 2 φ (φ + φ0) +# - none: f(φ, φ0) = 0 (not unconditionally stable) +# +# The stab enum in this script defines the schemes in terms of δψ to allow +# Nutils to construct the residuals through automatic differentiation. +# +# NOTE: This script uses dimensional quantities and requires the nutils.SI +# module, which is installed separate from the the core nutils. -from nutils import mesh, function, solver, export, cli, testing +from nutils import mesh, function, solver, numeric, export, cli, testing from nutils.expression_v2 import Namespace -import numpy, treelog, itertools, enum, typing +from nutils.SI import Length, Time, Density, Tension, Energy, Pressure, Velocity, parse +import numpy, itertools, enum +import treelog as log class stab(enum.Enum): - none = '0' # for educational purposes only - linear = '.5 dc^2 (6 - 6 c^2 + 8 c dc - 3 dc^2) / epsilon^2' - optimal = '.5 dc^2 (1 - dc^2 / 12) / epsilon^2' - -def main(nelems:int, etype:str, btype:str, degree:int, epsilon:typing.Optional[float], - contactangle:float, timestep:float, mtol:float, seed:int, circle:bool, stab:stab): + nonlinear = '.25 dφ^2 (1 - φ^2 + φ dφ (2 / 3) - dφ^2 / 6)' + linear = '.25 dφ^2 (2 + 2 abs(φ0) - (φ + φ0)^2)' + none = '0' # not unconditionally stable + +def main(size:Length, epsilon:Length, mobility:Time/Density, stens:Tension, + wtensn:Tension, wtensp:Tension, nelems:int, etype:str, degree:int, + timestep:Time, tol:Energy/Length, endtime:Time, seed:int, circle:bool, + stab:stab, showflux:bool): ''' Cahn-Hilliard equation on a unit square/circle. .. arguments:: - nelems [20] - Number of elements along domain edge. + size [10cm] + Domain size. + epsilon [2mm] + Interface thickness; defaults to an automatic value based on the + configured mesh density if left unspecified. + mobility [1mL*s/kg] + Mobility. + stens [50mN/m] + Surface tension. + wtensn [30mN/m] + Wall surface tension for phase -1. + wtensp [20mN/m] + Wall surface tension for phase +1. + nelems [0] + Number of elements along domain edge. When set to zero a value is set + automatically based on the configured domain size and epsilon. etype [square] Type of elements (square/triangle/mixed). - btype [std] - Type of basis function (std/spline), with availability depending on the - configured element type. degree [2] Polynomial degree. - epsilon [] - Interface thickness; defaults to an automatic value based on the - configured mesh density if left unspecified. - contactangle [90] - Wall contact angle in degrees. - timestep [.01] + timestep [.5s] Time step. - mtol [.01] - Threshold value for chemical potential peak to peak difference, used as - a stop criterion. + tol [1nJ/m] + Newton tolerance. + endtime [1min] + End of the simulation. seed [0] Random seed for the initial condition. circle [no] Select circular domain as opposed to a unit square. stab [linear] - Stabilization method (linear/optimal/none). + Stabilization method (linear/nonlinear/none). + showflux [no] + Overlay flux vectors on phase plot ''' - mineps = 1./nelems - if epsilon is None: - treelog.info('setting epsilon={}'.format(mineps)) - epsilon = mineps - elif epsilon < mineps: - treelog.warning('epsilon under crititical threshold: {} < {}'.format(epsilon, mineps)) + nmin = round(size / epsilon) + if nelems <= 0: + nelems = nmin + log.info('setting nelems to {}'.format(nelems)) + elif nelems < nmin: + log.warning('mesh is too coarse, consider increasing nelems to {:.0f}'.format(nmin)) + + log.info('contact angle: {:.0f}°'.format(numpy.arccos((wtensn - wtensp) / stens) * 180 / numpy.pi)) domain, geom = mesh.unitsquare(nelems, etype) - bezier = domain.sample('bezier', 5) # sample for plotting + if circle: + angle = (geom-.5) * (numpy.pi/2) + geom = .5 + function.sin(angle) * function.cos(angle)[[1,0]] / numpy.sqrt(2) + + bezier = domain.sample('bezier', 5) # sample for surface plots + grid = domain.locate(geom, numeric.simplex_grid([1,1], 1/40), maxdist=1/nelems, skip_missing=True, tol=1e-5) # sample for quivers + + φbasis = ηbasis = domain.basis('std', degree=degree) + ηbasis *= stens / epsilon # basis scaling to give η the required unit ns = Namespace() - if not circle: - ns.x = geom - else: - angle = (geom-.5) * (numpy.pi/2) - ns.x = function.sin(angle) * function.cos(angle)[[1,0]] / numpy.sqrt(2) + ns.x = size * geom ns.define_for('x', gradient='∇', normal='n', jacobians=('dV', 'dS')) - ns.epsilon = epsilon - ns.ewall = .5 * numpy.cos(contactangle * numpy.pi / 180) - ns.cbasis = ns.mbasis = domain.basis('std', degree=degree) - ns.c = function.dotarg('c', ns.cbasis) - ns.dc = ns.c - function.dotarg('c0', ns.cbasis) - ns.m = function.dotarg('m', ns.mbasis) - ns.F = '.5 (c^2 - 1)^2 / epsilon^2' - ns.dF = stab.value + ns.ε = epsilon + ns.σ = stens + ns.φ = function.dotarg('φ', φbasis) + ns.σmean = (wtensp + wtensn) / 2 + ns.σdiff = (wtensp - wtensn) / 2 + ns.σwall = 'σmean + φ σdiff' + ns.φ0 = function.dotarg('φ0', φbasis) + ns.dφ = 'φ - φ0' + ns.η = function.dotarg('η', ηbasis) + ns.ψ = '.25 (φ^2 - 1)^2' + ns.δψ = stab.value + ns.M = mobility + ns.J_i = '-M ∇_i(η)' ns.dt = timestep - nrg_mix = domain.integral('F dV' @ ns, degree=7) - nrg_iface = domain.integral('.5 ∇_k(c) ∇_k(c) dV' @ ns, degree=7) - nrg_wall = domain.boundary.integral('(abs(ewall) + c ewall) dS' @ ns, degree=7) - nrg = nrg_mix + nrg_iface + nrg_wall + domain.integral('(dF - m dc - .5 dt epsilon^2 ∇_k(m) ∇_k(m)) dV' @ ns, degree=7) + nrg_mix = domain.integral('(ψ σ / ε) dV' @ ns, degree=7) + nrg_iface = domain.integral('.5 σ ε ∇_k(φ) ∇_k(φ) dV' @ ns, degree=7) + nrg_wall = domain.boundary.integral('σwall dS' @ ns, degree=7) + nrg = nrg_mix + nrg_iface + nrg_wall + domain.integral('(δψ σ / ε - η dφ + .5 dt J_k ∇_k(η)) dV' @ ns, degree=7) numpy.random.seed(seed) - state = dict(c=numpy.random.normal(0,.5,ns.cbasis.shape), m=numpy.random.normal(0,.5,ns.mbasis.shape)) # initial condition + state = dict(φ=numpy.random.normal(0, .5, φbasis.shape)) # initial condition - with treelog.iter.plain('timestep', itertools.count()) as steps: + with log.iter.fraction('timestep', range(round(endtime / timestep))) as steps: for istep in steps: - E = function.eval([nrg_mix, nrg_iface, nrg_wall], **state) - treelog.user('energy: {0:.3f} ({1[0]:.0f}% mixture, {1[1]:.0f}% interface, {1[2]:.0f}% wall)'.format(sum(E), 100*numpy.array(E)/sum(E))) - - x, c, m = bezier.eval(['x_i', 'c', 'm'] @ ns, **state) - export.triplot('phase.png', x, c, tri=bezier.tri, clim=(-1,1)) - export.triplot('chempot.png', x, m, tri=bezier.tri) - - if numpy.ptp(m) < mtol: - break - - state['c0'] = state['c'] - state = solver.optimize(['c', 'm'], nrg, arguments=state, tol=1e-10) + E = numpy.stack(function.eval([nrg_mix, nrg_iface, nrg_wall], **state)) + log.user('energy: {0:,.0μJ/m} ({1[0]:.0f}% mixture, {1[1]:.0f}% interface, {1[2]:.0f}% wall)'.format(numpy.sum(E), 100*E/numpy.sum(E))) + + state['φ0'] = state['φ'] + state = solver.optimize(['φ', 'η'], nrg / tol, arguments=state, tol=1) + + with export.mplfigure('phase.png') as fig: + ax = fig.add_subplot(aspect='equal', xlabel='[mm]', ylabel='[mm]') + x, φ = bezier.eval(['x_i', 'φ'] @ ns, **state) + im = ax.tripcolor(*(x/'mm').T, bezier.tri, φ, shading='gouraud', cmap='bwr') + im.set_clim(-1, 1) + fig.colorbar(im) + if showflux: + x, J = grid.eval(['x_i', 'J_i'] @ ns, **state) + log.info('largest flux: {:.1mm/h}'.format(numpy.max(numpy.hypot(J[:,0], J[:,1])))) + ax.quiver(*(x/'mm').T, *(J/'m/s').T, color='r') + ax.quiver(*(x/'mm').T, *-(J/'m/s').T, color='b') + ax.autoscale(enable=True, axis='both', tight=True) return state @@ -114,37 +216,34 @@ def main(nelems:int, etype:str, btype:str, degree:int, epsilon:typing.Optional[f class test(testing.TestCase): def test_initial(self): - state = main(nelems=3, etype='square', btype='std', degree=2, epsilon=None, contactangle=90, timestep=1, mtol=float('inf'), seed=0, circle=False, stab=stab.linear) - with self.subTest('concentration'): self.assertAlmostEqual64(state['c'], ''' + state = main(size=parse('10cm'), epsilon=parse('5cm'), mobility=parse('1μL*s/kg'), + stens=parse('50mN/m'), wtensn=parse('30mN/m'), wtensp=parse('20mN/m'), nelems=3, + etype='square', degree=2, timestep=parse('1h'), tol=parse('1nJ/m'), + endtime=parse('1h'), seed=0, circle=False, stab=stab.linear, showflux=True) + with self.subTest('concentration'): self.assertAlmostEqual64(state['φ0'], ''' eNoBYgCd/xM3LjTtNYs3MDcUyt41uc14zjo0LzKzNm812jFhNNMzwDYgzbMzV8o0yCM1rzWeypE3Tcnx L07NzTa4NlMyETREyrPIGMxYMl82VDbjy1/M8clZyf3IRjday6XLmMl6NRnJMF4tqQ==''') - with self.subTest('chemical-potential'): self.assertAlmostEqual64(state['m'], ''' - eNoBYgCd/w7NQModNFjLtckB0VA0rDCiM+zKBMzPygjMcMr4yJcy0MsUyXY0OcoxMFo19zE5Np/JMDTG - yk7KGstQzFgwvMnDNXk0Msm+Njc3SjZizebJEjbPy1w25zLsNfQzSjURLRk3Qt4uBQ==''') def test_square(self): - state = main(nelems=3, etype='square', btype='std', degree=2, epsilon=None, contactangle=90, timestep=1, mtol=.1, seed=0, circle=False, stab=stab.linear) - with self.subTest('concentration'): self.assertAlmostEqual64(state['c'], ''' - eNoBYgCd/5o2nDZANsQ1dTQO1QXQmTaaNjw2vDVVNOrPqM5mNmY2xDWwNPDRNMsZyys2KjYdNQoyh8v9 - yffJ1zXTNSA0Ks2JyorJicluNWQ1IDJXy+/JNck3yWQ1WjX0MVHL78k3yTnJbYgt7Q==''') - with self.subTest('chemical-potential'): self.assertAlmostEqual64(state['m'], ''' - eNpVyU0KgCAABeGrK9SiQCiQCsrwB8GjzLF8uJNvNYzBsuP5yGKmWlhxXKMSmxzchPGceF4iVU55+Ck0 - masD28JDNQ==''') - - def test_contactangle(self): - state = main(nelems=3, etype='square', btype='std', degree=2, epsilon=None, contactangle=45, timestep=1, mtol=.1, seed=0, circle=False, stab=stab.linear) - with self.subTest('concentration'): self.assertAlmostEqual64(state['c'], ''' - eNoBYgCd/zs2bjYbNqs19zNqzRXMazaZNkw25zVxNLXOncxDNnU28zUqNeIwTsvKyhc2TDaANdozSMwp - yuXJojXiNYo05c/Oyp3Jb8n7NE81iDK2ywfKO8kayXI03jQjLyrLzskZyfrIkb8vTg==''') - with self.subTest('chemical-potential'): self.assertAlmostEqual64(state['m'], ''' - eNoNzCEOggAYgNFLvbPguILeQRpWgyQPoAYCuBGhWHVzGshOLegMSvffS1/5EqmF2sUnTKJylbO3l6mZ - pb2TZ5jLrDWObmGlUDroDWFjq43H3dcvaqdz9TCGP1tYLVU=''') + state = main(size=parse('10cm'), epsilon=parse('5cm'), mobility=parse('1μL*s/kg'), + stens=parse('50mN/m'), wtensn=parse('30mN/m'), wtensp=parse('20mN/m'), nelems=3, + etype='square', degree=2, timestep=parse('1h'), tol=parse('1nJ/m'), + endtime=parse('2h'), seed=0, circle=False, stab=stab.linear, showflux=True) + with self.subTest('concentration'): self.assertAlmostEqual64(state['φ'], ''' + eNoBYgCd/zE1EjX1NaA2+TXiMxkz0TS9NL01ajaRNZoxYNElNRM1LDUlNZQw0cqgysI1nTWcNN4xLsuk + ybDJvDWaNTQ07s7nysnJ6ckPNQY1CzNozKjK58kOysQ0zTQKM73M3coVyjfKR9cuPg==''') + with self.subTest('chemical-potential'): self.assertAlmostEqual64(state['η'], ''' + eNoBYgCd/3TIkccBNkQ6IDqIN3/MF8cSx9Y02TmdOVHLMcecxxLIEjQUOAHOa8a1xWw3izb2M9UzPMc0 + xmnGpzibODY34tETyJHHp8hbyWU2xzZTydfIOsrNyo3Gi8jCyyXIm8hkzD3K1IAxtQ==''') def test_mixedcircle(self): - state = main(nelems=3, etype='mixed', btype='std', degree=2, epsilon=None, contactangle=90, timestep=1, mtol=.1, seed=0, circle=True, stab=stab.linear) - with self.subTest('concentration'): self.assertAlmostEqual64(state['c'], ''' - eNoBYgCd/y415zQQNUw1dzVTNB01fTS4My81dDQpNWU03DJENDo0rjEVMCfMj8vjysUzozMmz3E0VCyg - 1I00MDMWM3bQU8t4ym7Kd8uKzoTPgc7PzmfOYs3ZzJvM6jIozunL7sqezgXMYsUuJg==''') - with self.subTest('chemical-potential'): self.assertAlmostEqual64(state['m'], ''' - eNoljEkKgEAQxP7/B/EkiCCO27gOuJxy8FEGpSEUoaoTkYWTgyAjHRuTZpGzHGgpKchMMz2JRnN/l6j1 - uctge2W08W8fdlPF5ccXGqBI7Q==''') + state = main(size=parse('10cm'), epsilon=parse('5cm'), mobility=parse('1μL*s/kg'), + stens=parse('50mN/m'), wtensn=parse('30mN/m'), wtensp=parse('20mN/m'), nelems=3, + etype='mixed', degree=2, timestep=parse('1h'), tol=parse('1nJ/m'), + endtime=parse('2h'), seed=0, circle=True, stab=stab.linear, showflux=True) + with self.subTest('concentration'): self.assertAlmostEqual64(state['φ'], ''' + eNoBYgCd/w01AjX+NAw1IjXTNMw0iTRPNDI1vDQcNTk0uzJ9NFM0HS4P0SbMcssOy0wzZjNw0b0zljHK + z6U0ps8zM/LPjspVypDKUsuLzk3MgM3OzYnN7s/61KfP2zH4MADNhst3z7DMoBcvyQ==''') + with self.subTest('chemical-potential'): self.assertAlmostEqual64(state['η'], ''' + eNoBYgCd/+s1Bcp4ztI3gjYFyZk4YzVjyfA2AzdAMj032zfLNTE4fMm7yLnGisbqxZPJ2MsfyD81csiv + x+E5xDhjOJA3msZ1xZTFa8ddx/fG88eCx73H1MieM/c0WDihMUrLvMYZNpvIrWQ0sw==''') diff --git a/nutils/export.py b/nutils/export.py index 474dbe65f..a7894f0b9 100644 --- a/nutils/export.py +++ b/nutils/export.py @@ -54,15 +54,20 @@ def mplfigure(name, kwargs=...): finally: fig.set_canvas(None) # break circular reference -def triplot(name, points, values=None, *, tri=None, hull=None, cmap=None, clim=None, linewidth=.1, linecolor='k'): +def plotlines_(ax, xy, lines, **kwargs): + from matplotlib import collections + lc = collections.LineCollection(numpy.asarray(xy).T[lines], **kwargs) + ax.add_collection(lc) + return lc + +def triplot(name, points, values=None, *, tri=None, hull=None, cmap=None, clim=None, linewidth=.1, linecolor='k', plabel=None, vlabel=None): if (tri is None) != (values is None): raise Exception('tri and values can only be specified jointly') with mplfigure(name) as fig: - ax = fig.add_subplot(111) if points.shape[1] == 1: + ax = fig.add_subplot(111, xlabel=plabel, ylabel=vlabel) if tri is not None: - import matplotlib.collections - ax.add_collection(matplotlib.collections.LineCollection(numpy.array([points[:,0], values]).T[tri])) + plotlines_(ax, [points[:,0], values], tri) if hull is not None: for x in points[hull[:,0],0]: ax.axvline(x, color=linecolor, linewidth=linewidth) @@ -72,15 +77,14 @@ def triplot(name, points, values=None, *, tri=None, hull=None, cmap=None, clim=N else: ax.set_ylim(clim) elif points.shape[1] == 2: - ax.set_aspect('equal') + ax = fig.add_subplot(111, xlabel=plabel, ylabel=plabel, aspect='equal') if tri is not None: - im = ax.tripcolor(points[:,0], points[:,1], tri, values, shading='gouraud', cmap=cmap, rasterized=True) + im = ax.tripcolor(*points.T, tri, values, shading='gouraud', cmap=cmap, rasterized=True) if clim is not None: im.set_clim(clim) - fig.colorbar(im) + fig.colorbar(im, label=vlabel) if hull is not None: - import matplotlib.collections - ax.add_collection(matplotlib.collections.LineCollection(points[hull], colors=linecolor, linewidths=linewidth, alpha=1 if tri is None else .5)) + plotlines_(ax, points.T, hull, colors=linecolor, linewidths=linewidth, alpha=1 if tri is None else .5) ax.autoscale(enable=True, axis='both', tight=True) else: raise Exception('invalid spatial dimension: {}'.format(points.shape[1])) diff --git a/nutils/expression_v2.py b/nutils/expression_v2.py index b901b820e..7fdddb8c8 100644 --- a/nutils/expression_v2.py +++ b/nutils/expression_v2.py @@ -595,33 +595,33 @@ class Namespace: def __init__(self) -> None: self.opposite = function.opposite - self.sin = function.sin - self.cos = function.cos - self.tan = function.tan - self.sinh = function.sinh - self.cosh = function.cosh - self.tanh = function.tanh - self.arcsin = function.arcsin - self.arccos = function.arccos - self.arctan = function.arctan - self.arctanh = function.arctanh - self.exp = function.exp - self.abs = function.abs - self.ln = function.ln - self.log = function.ln - self.log2 = function.log2 - self.log10 = function.log10 - self.sqrt = function.sqrt - self.sign = function.sign - self.conj = function.conj - self.real = function.real - self.imag = function.imag + self.sin = numpy.sin + self.cos = numpy.cos + self.tan = numpy.tan + self.sinh = numpy.sinh + self.cosh = numpy.cosh + self.tanh = numpy.tanh + self.arcsin = numpy.arcsin + self.arccos = numpy.arccos + self.arctan = numpy.arctan + self.arctanh = numpy.arctanh + self.exp = numpy.exp + self.abs = numpy.abs + self.ln = numpy.log + self.log = numpy.log + self.log2 = numpy.log2 + self.log10 = numpy.log10 + self.sqrt = numpy.sqrt + self.sign = numpy.sign + self.conj = numpy.conj + self.real = numpy.real + self.imag = numpy.imag def __setattr__(self, attr: str, value: Union[function.Array, str]) -> None: name, underscore, indices = attr.partition('_') if isinstance(value, (int, float, complex, numpy.ndarray)): value = function.Array.cast(value) - if isinstance(value, function.Array): + if hasattr(value, '__array_ufunc__') and hasattr(value, '__array_function__'): if underscore: raise AttributeError('Cannot assign an array to an attribute with an underscore.') super().__setattr__(name, value) @@ -632,7 +632,7 @@ def __setattr__(self, attr: str, value: Union[function.Array, str]) -> None: raise AttributeError('All indices must be unique.') ops = _FunctionArrayOps(self) array, shape, expression_indices, summed = _Parser(ops).parse_expression(_Substring(value)) - assert array.shape == shape + assert numpy.shape(array) == shape if expression_indices != indices: for index in sorted(set(indices) - set(expression_indices)): raise AttributeError('Index {} of the namespace attribute is missing in the expression.'.format(index)) @@ -640,7 +640,7 @@ def __setattr__(self, attr: str, value: Union[function.Array, str]) -> None: raise AttributeError('Index {} of the expression is missing in the namespace attribute.'.format(index)) array = ops.align(array, expression_indices, indices) super().__setattr__(name, array) - elif isinstance(value, Callable): + elif callable(value): if underscore: raise AttributeError('Cannot assign a function to an attribute with an underscore.') super().__setattr__(name, value) @@ -652,7 +652,7 @@ def __rmatmul__(self, expression): parser = _Parser(ops) if isinstance(expression, str): array, shape, indices, summed = parser.parse_expression(_Substring(expression)) - assert array.shape == shape + assert numpy.shape(array) == shape array = ops.align(array, indices, ''.join(sorted(indices))) return array elif isinstance(expression, tuple): @@ -702,17 +702,17 @@ def define_for(self, __name: str, *, gradient: Optional[str] = None, curl: Optio if gradient: setattr(self, gradient, lambda arg: function.grad(arg, geom)) if curl: - if geom.shape != (3,): - raise ValueError('The curl can only be defined for a geometry with shape (3,) but got {}.'.format(geom.shape)) + if numpy.shape(geom) != (3,): + raise ValueError('The curl can only be defined for a geometry with shape (3,) but got {}.'.format(numpy.shape(geom))) # Definition: `curl_ki(u_...)` := `ε_kji ∇_j(u_...)`. Should be used as # `curl_ki(u_i)`, which is equivalent to `ε_kji ∇_j(u_i)`. setattr(self, curl, lambda arg: (function.levicivita(3) * function.grad(arg, geom)[...,numpy.newaxis,:,numpy.newaxis]).sum(-2)) if normal: setattr(self, normal, function.normal(geom)) for i, jacobian in enumerate(jacobians): - if i > geom.size: + if i > numpy.size(geom): raise ValueError('Cannot define the jacobian {!r}: dimension is negative.'.format(jacobian)) - setattr(self, jacobian, function.jacobian(geom, geom.size - i)) + setattr(self, jacobian, function.jacobian(geom, numpy.size(geom) - i)) def copy_(self, **replacements: Mapping[str, function.Array]) -> 'Namespace': '''Return a copy of this namespace. @@ -755,12 +755,12 @@ def get_variable(self, name: str, ndim: int) -> Optional[Union[Tuple[function.Ar array = getattr(self.namespace, name) except AttributeError: return None - if not isinstance(array, function.Array): + if callable(array): return None - if array.ndim == ndim: - return array, array.shape + elif numpy.ndim(array) == ndim: + return array, numpy.shape(array) else: - return _InvalidDimension(array.ndim) + return _InvalidDimension(numpy.ndim(array)) def call(self, name: str, ngenerates: int, arg: function.Array) -> Optional[Union[Tuple[function.Array, _Shape],_InvalidDimension]]: try: @@ -768,23 +768,22 @@ def call(self, name: str, ngenerates: int, arg: function.Array) -> Optional[Unio except AttributeError: return None array = func(arg) - assert isinstance(array, function.Array) - assert array.shape[:arg.ndim] == arg.shape - if array.ndim == arg.ndim + ngenerates: - return array, array.shape[arg.ndim:] + assert numpy.shape(array)[:numpy.ndim(arg)] == numpy.shape(arg) + if numpy.ndim(array) == numpy.ndim(arg) + ngenerates: + return array, numpy.shape(array)[numpy.ndim(arg):] else: - return _InvalidDimension(array.ndim - arg.ndim) + return _InvalidDimension(numpy.ndim(array) - numpy.ndim(arg)) def get_element(self, array: function.Array, axis: int, index: int) -> function.Array: - assert 0 <= axis < array.ndim and 0 <= index < array.shape[axis] - return function.get(array, axis, index) + assert 0 <= axis < numpy.ndim(array) and 0 <= index < numpy.shape(array)[axis] + return numpy.take(array, index, axis) def transpose(self, array: function.Array, axes: Tuple[int, ...]) -> function.Array: - assert array.ndim == len(axes) - return function.transpose(array, axes) + assert numpy.ndim(array) == len(axes) + return numpy.transpose(array, axes) def trace(self, array: function.Array, axis1: int, axis2: int) -> function.Array: - return function.trace(array, axis1, axis2) + return numpy.trace(array, axis1, axis2) def scope(self, array: function.Array) -> function.Array: return array @@ -796,20 +795,24 @@ def jump(self, array: function.Array) -> function.Array: return function.jump(array) def add(self, *args: Tuple[bool, function.Array]) -> function.Array: - assert all(arg.shape == args[0][1].shape for neg, arg in args[1:]) + assert all(numpy.shape(arg) == numpy.shape(args[0][1]) for neg, arg in args[1:]) negated = (-arg if neg else arg for neg, arg in args) - return functools.reduce(function.add, negated) + return functools.reduce(numpy.add, negated) + + def append_axes(self, array, shape): + shuffle = numpy.concatenate([len(shape) + numpy.arange(numpy.ndim(array)), numpy.arange(len(shape))]) + return numpy.transpose(numpy.broadcast_to(array, shape + numpy.shape(array)), shuffle) def multiply(self, *args: function.Array) -> function.Array: result = args[0] for arg in args[1:]: - result = function.multiply(function._append_axes(result, arg.shape), function._prepend_axes(arg, result.shape)) + result = numpy.multiply(self.append_axes(result, numpy.shape(arg)), arg) return result def divide(self, numerator: function.Array, denominator: function.Array) -> function.Array: - assert denominator.ndim == 0 - return function.divide(numerator, function._append_axes(denominator, numerator.shape)) + assert numpy.ndim(denominator) == 0 + return numpy.true_divide(numerator, denominator) def power(self, base: function.Array, exponent: function.Array) -> function.Array: - assert exponent.ndim == 0 - return function.power(base, function._append_axes(exponent, base.shape)) + assert numpy.ndim(exponent) == 0 + return numpy.power(base, exponent) diff --git a/nutils/function.py b/nutils/function.py index d37f4fea6..6dfee58a6 100644 --- a/nutils/function.py +++ b/nutils/function.py @@ -141,7 +141,7 @@ def __array_ufunc__(self, ufunc, method, *inputs, out=None, **kwargs): if method != '__call__' or ufunc not in HANDLED_FUNCTIONS: return NotImplemented try: - arrays = [Array.cast(v) for v in inputs] + arrays = [v if isinstance(v, (Array, bool, int, float, complex, numpy.ndarray)) else Array.cast(v) for v in inputs] except ValueError: return NotImplemented return HANDLED_FUNCTIONS[ufunc](*arrays, **kwargs) @@ -177,6 +177,17 @@ def cast(cls, __value: IntoArray, dtype: Optional[DType] = None, ndim: Optional[ raise ValueError('expected an array with dimension `{}` but got `{}`'.format(ndim, value.ndim)) return value + @classmethod + def cast_withscale(cls, __value: IntoArray, dtype: Optional[DType] = None, ndim: Optional[int] = None): + try: + scale = type(__value).reference_quantity + except AttributeError: + value = cls.cast(__value, dtype=dtype, ndim=ndim) + scale = value.dtype(1) + else: + value = cls.cast(__value / scale, dtype=dtype, ndim=ndim) + return value, scale + def __init__(self, shape: Shape, dtype: DType, spaces: FrozenSet[str], arguments: Mapping[str, Tuple[Shape, DType]]) -> None: self.shape = tuple(sh.__index__() for sh in shape) self.dtype = dtype @@ -1074,6 +1085,18 @@ def asarray(__arg: IntoArray) -> Array: return Array.cast(__arg) +@implements(numpy.shape) +def shape(__arg: IntoArray) -> Tuple[int, ...]: + return asarray(__arg).shape + +@implements(numpy.ndim) +def ndim(__arg: IntoArray) -> int: + return asarray(__arg).ndim + +@implements(numpy.size) +def size(__arg: IntoArray) -> int: + return asarray(__arg).size + def zeros(shape: Shape, dtype: DType = float) -> Array: '''Create a new :class:`Array` of given shape and dtype, filled with zeros. @@ -1217,7 +1240,10 @@ def multiply(__left: IntoArray, __right: IntoArray) -> Array: :class:`Array` ''' - return _Wrapper.broadcasted_arrays(evaluable.multiply, __left, __right) + left, right = typecast_arrays(__left, __right) + return right if _is_unit_scalar(__left) \ + else left if _is_unit_scalar(__right) \ + else _Wrapper.broadcasted_arrays(evaluable.multiply, left, right) @implements(numpy.true_divide) def divide(__dividend: IntoArray, __divisor: IntoArray) -> Array: @@ -1232,7 +1258,10 @@ def divide(__dividend: IntoArray, __divisor: IntoArray) -> Array: :class:`Array` ''' - return multiply(__dividend, reciprocal(__divisor)) + dividend, divisor = typecast_arrays(__dividend, __divisor, min_dtype=float) + return dividend if _is_unit_scalar(__divisor) \ + else reciprocal(divisor) if _is_unit_scalar(__dividend) \ + else _Wrapper.broadcasted_arrays(evaluable.divide, dividend, divisor) @implements(numpy.floor_divide) def floor_divide(__dividend: IntoArray, __divisor: IntoArray) -> Array: @@ -2714,7 +2743,7 @@ def derivative(__arg: IntoArray, __var: Union[str, 'Argument']) -> Array: :class:`Array` ''' - arg = Array.cast(__arg) + arg, argscale = Array.cast_withscale(__arg) if isinstance(__var, str): if __var not in arg.arguments: raise ValueError('no such argument: {}'.format(__var)) @@ -2728,7 +2757,7 @@ def derivative(__arg: IntoArray, __var: Union[str, 'Argument']) -> Array: raise ValueError('Argument {!r} has shape {} in the function, but the derivative to {!r} with shape {} was requested.'.format(__var.name, shape, __var.name, __var.shape)) if __var.dtype != dtype: raise ValueError('Argument {!r} has dtype {} in the function, but the derivative to {!r} with dtype {} was requested.'.format(__var.name, dtype.__name__ if dtype in _dtypes else dtype, __var.name, __var.dtype.__name__ if __var.dtype in _dtypes else __var.dtype)) - return _Derivative(arg, __var) + return _Derivative(arg, __var) * argscale def grad(__arg: IntoArray, __geom: IntoArray, ndims: int = 0) -> Array: '''Return the gradient of the argument to the given geometry. @@ -2744,21 +2773,22 @@ def grad(__arg: IntoArray, __geom: IntoArray, ndims: int = 0) -> Array: :class:`Array` ''' - arg = Array.cast(__arg) - geom = Array.cast(__geom) + arg, argscale = Array.cast_withscale(__arg) + geom, geomscale = Array.cast_withscale(__geom) if geom.dtype != float: raise ValueError('The geometry must be real-valued.') if geom.ndim == 0: - return grad(arg, _append_axes(geom, (1,)))[...,0] + ret = grad(arg, _append_axes(geom, (1,)))[...,0] elif geom.ndim > 1: sh = geom.shape[-2:] - return unravel(grad(arg, ravel(geom, geom.ndim-2)), arg.ndim+geom.ndim-2, sh) + ret = unravel(grad(arg, ravel(geom, geom.ndim-2)), arg.ndim+geom.ndim-2, sh) elif ndims == 0 or ndims == geom.shape[0]: - return _Gradient(arg, geom) + ret = _Gradient(arg, geom) elif ndims == -1 or ndims == geom.shape[0] - 1: - return _SurfaceGradient(arg, geom) + ret = _SurfaceGradient(arg, geom) else: raise NotImplementedError + return ret * (argscale / geomscale) def curl(__arg: IntoArray, __geom: IntoArray) -> Array: '''Return the curl of the argument w.r.t. the given geometry. @@ -2801,7 +2831,7 @@ def normal(__geom: IntoArray, refgeom: Optional[Array] = None) -> Array: :class:`Array` ''' - geom = Array.cast(__geom) + geom, geomscale_ = Array.cast_withscale(__geom) if geom.dtype != float: raise ValueError('The geometry must be real-valued.') if geom.ndim == 0: @@ -2870,15 +2900,22 @@ def jacobian(__geom: IntoArray, __ndims: Optional[int] = None) -> Array: :class:`Array` ''' - geom = Array.cast(__geom) + geom, geomscale = Array.cast_withscale(__geom) if geom.dtype != float: raise ValueError('The geometry must be real-valued.') + if __ndims is not None: + scale = geomscale**__ndims + elif _is_unit_scalar(geomscale): + scale = 1. + else: + raise ValueError('scaled arrays require an explicit dimension') if geom.ndim == 0: - return jacobian(insertaxis(geom, 0, 1), __ndims) + J = jacobian(insertaxis(geom, 0, 1), __ndims) elif geom.ndim > 1: - return jacobian(ravel(geom, geom.ndim-2), __ndims) + J = jacobian(ravel(geom, geom.ndim-2), __ndims) else: - return _Jacobian(geom, __ndims) + J = _Jacobian(geom, __ndims) + return J * scale def J(__geom: IntoArray, __ndims: Optional[int] = None) -> Array: '''Return the absolute value of the determinant of the Jacobian matrix of the given geometry. @@ -3026,7 +3063,8 @@ def eval(funcs: evaluable.AsEvaluableArray, **arguments: Mapping[str, numpy.ndar results : :class:`tuple` of arrays ''' - return map(sparse.toarray, evaluable.eval_sparse(funcs, **arguments)) + funcs, funcscales = zip(*map(Array.cast_withscale, funcs)) + return tuple(sparse.toarray(data) * scale for data, scale in zip(evaluable.eval_sparse(funcs, **arguments), funcscales)) def isarray(__arg: Any) -> bool: 'Test if the argument is an instance of :class:`Array`.' @@ -3636,3 +3674,7 @@ def f_dofs_coeffs(self, index: evaluable.Array) -> Tuple[evaluable.Array,evaluab def Namespace(*args, **kwargs): from .expression_v1 import Namespace return Namespace(*args, **kwargs) + +def _is_unit_scalar(v): + T = type(v) + return T in _dtypes and v == T(1) diff --git a/nutils/numeric.py b/nutils/numeric.py index da55e0698..a9b72f0d0 100644 --- a/nutils/numeric.py +++ b/nutils/numeric.py @@ -145,6 +145,54 @@ def meshgrid(*args, dtype=None): assert n == 0 return grid +def _simplex_grid_helper(bounds): + if bounds.ndim != 1 or len(bounds) == 0: + raise ValueError + nd = len(bounds) + spacing = [numpy.sqrt((1+i/2)/(1+i)) for i in range(nd)] # simplex height orthogonal to lower dimension + grid = meshgrid(*[numpy.arange(bound, step=step) for step, bound in zip(spacing, bounds)]) + out_of_bounds = [] + for idim in range(nd-1): + stripes = grid[(idim,)+(slice(None),)*(idim+1)+(slice(1,None,2),)] + stripes += spacing[idim] * (idim+1) / (idim+2) + if stripes.size and stripes.flat[-1] >= bounds[idim]: + out_of_bounds.append(idim) + if out_of_bounds: + select = numpy.ones(grid.shape[1:], dtype=bool) + for idim in out_of_bounds: + select[(slice(None),)*(idim)+(-1,)+(slice(1,None,2),)] = False + points = grid[:,select].T + else: + points = grid.reshape(nd, -1).T + d = numpy.subtract(bounds, points.max(axis=0)) + assert (d > 0).all() + points += d / 2 + return points + +def simplex_grid(shape, spacing): + '''Multi-dimensional generator for equilateral simplex grids. + + Simplex_grid generates a point cloud within an n-dimensional orthotope, which + ranges from zero to a specified shape. The point coordinates are spaced in + such a way that the nearest neighbours are at distance `spacing`, thus + forming vertices of regular simplices. The returned array is two-dimensional, + with the first axis being the spatial dimension (matching `shape`) and the + second a stacking of the generated points. + + Parameters + ---------- + shape : :class:`tuple` + list or tuple of dimensions of the orthotope to be filled. + spacing : :class:`float` + minimum spacing in the generated point cloud. + + Returns + ------- + :class:`numpy.ndarray` + ''' + + return _simplex_grid_helper(numpy.divide(shape, spacing)) * spacing + def takediag(A, axis=-2, rmaxis=-1): axis = normdim(A.ndim, axis) rmaxis = normdim(A.ndim, rmaxis) diff --git a/nutils/sample.py b/nutils/sample.py index 469762cd7..0e155a6b5 100644 --- a/nutils/sample.py +++ b/nutils/sample.py @@ -189,9 +189,10 @@ def integrate(self, funcs: Iterable[function.IntoArray], arguments: Mapping[str, Optional arguments for function evaluation. ''' + funcs, funcscales = zip(*map(function.Array.cast_withscale, funcs)) datas = self.integrate_sparse(funcs, argdict(arguments)) with log.iter.fraction('assembling', datas) as items: - return tuple(_convert(data, inplace=True) for data in items) + return tuple(_convert(data, inplace=True) * scale for data, scale in zip(items, funcscales)) @util.single_or_multiple def integrate_sparse(self, funcs: Iterable[function.IntoArray], arguments: Optional[Mapping[str, numpy.ndarray]] = None) -> Tuple[numpy.ndarray, ...]: @@ -216,7 +217,8 @@ def integral(self, __func: function.IntoArray) -> function.Array: Integrand. ''' - return _Integral(function.Array.cast(__func), self) + func, funcscale = function.Array.cast_withscale(__func) + return _Integral(func, self) * funcscale @util.positional_only @util.single_or_multiple @@ -231,9 +233,10 @@ def eval(self, funcs: Iterable[function.IntoArray], arguments: Mapping[str, nump Optional arguments for function evaluation. ''' + funcs, funcscales = zip(*map(function.Array.cast_withscale, funcs)) datas = self.eval_sparse(funcs, arguments) with log.iter.fraction('assembling', datas) as items: - return tuple(map(sparse.toarray, items)) + return tuple(sparse.toarray(data) * funcscale for data, funcscale in zip(datas, funcscales)) @util.positional_only @util.single_or_multiple @@ -251,10 +254,10 @@ def eval_sparse(self, funcs: Iterable[function.IntoArray], arguments: Optional[M return evaluable.eval_sparse(map(self, funcs), **(arguments or {})) def __call__(self, __func: function.IntoArray) -> function.Array: - func = _ConcatenatePoints(function.Array.cast(__func), self) + func, funcscale = function.Array.cast_withscale(__func) ielem = evaluable.loop_index('_sample_' + '_'.join(self.spaces), self.nelems) indices = evaluable.loop_concatenate(evaluable._flat(self.get_evaluable_indices(ielem)), ielem) - return _ReorderPoints(func, indices) + return _ReorderPoints(_ConcatenatePoints(func, self), indices) * funcscale def basis(self) -> function.Array: '''Basis-like function that for every point in the sample evaluates to the diff --git a/nutils/topology.py b/nutils/topology.py index 66b3fac5a..5463500a8 100644 --- a/nutils/topology.py +++ b/nutils/topology.py @@ -697,7 +697,7 @@ def select(self, indicator: function.Array, ischeme: str = 'bezier2', **kwargs: selected = types.frozenarray(numpy.unique(ielem[isactive])) return self[selected] - def locate(self, geom, coords, *, tol=0, eps=0, maxiter=0, arguments=None, weights=None, maxdist=None, ischeme=None, scale=None) -> Sample: + def locate(self, geom, coords, *, tol=0, eps=0, maxiter=0, arguments=None, weights=None, maxdist=None, ischeme=None, scale=None, skip_missing=False) -> Sample: '''Create a sample based on physical coordinates. In a finite element application, functions are commonly evaluated in points @@ -745,6 +745,10 @@ def locate(self, geom, coords, *, tol=0, eps=0, maxiter=0, arguments=None, weigh element centroid above which the element is rejected immediately. If all points are expected to be located then this can safely be left unspecified. + skip_missing : :class:`bool` (default: False) + When set to true, skip points that are not found (for instance because + they fall outside the domain) in the returned sample. When set to false + (the default) missing points raise a ``LocateError``. Returns ------- @@ -950,7 +954,7 @@ def withgroups(self, vgroups: Mapping[str, Union[str, Topology]] = {}, bgroups: def indicator(self, subtopo: Union[str, Topology]) -> Topology: raise SkipTest('`{}` does not implement `Topology.indicator`'.format(type(self).__qualname__)) - def locate(self, geom, coords, *, tol=0, eps=0, maxiter=0, arguments=None, weights=None, maxdist=None, ischeme=None, scale=None) -> Sample: + def locate(self, geom, coords, *, tol=0, eps=0, maxiter=0, arguments=None, weights=None, maxdist=None, ischeme=None, scale=None, skip_missing=False) -> Sample: raise SkipTest('`{}` does not implement `Topology.locate`'.format(type(self).__qualname__)) else: @@ -1269,8 +1273,8 @@ def indicator(self, subtopo: Union[str, Topology]) -> Topology: else: return super().indicator(subtopo) - def locate(self, geom, coords, *, tol=0, eps=0, maxiter=0, arguments=None, weights=None, maxdist=None, ischeme=None, scale=None) -> Sample: - return self.parent.locate(geom, coords, tol=tol, eps=eps, maxiter=maxiter, arguments=arguments, weights=weights, maxdist=maxdist, ischeme=ischeme, scale=scale) + def locate(self, geom, coords, *, tol=0, eps=0, maxiter=0, arguments=None, weights=None, maxdist=None, ischeme=None, scale=None, skip_missing=False) -> Sample: + return self.parent.locate(geom, coords, tol=tol, eps=eps, maxiter=maxiter, arguments=arguments, weights=weights, maxdist=maxdist, ischeme=ischeme, scale=scale, skip_missing=skip_missing) def boundary_spaces_unchecked(self, spaces: FrozenSet[str]) -> Topology: return _WithGroupAliases(self.parent.boundary_spaces_unchecked(spaces), self.bgroups, types.frozendict({}), types.frozendict({})) @@ -1473,7 +1477,7 @@ def indicator(self, subtopo): return function.get(values, 0, self.f_index) @log.withcontext - def locate(self, geom, coords, *, tol=0, eps=0, maxiter=0, arguments=None, weights=None, maxdist=None, ischeme=None, scale=None): + def locate(self, geom, coords, *, tol=0, eps=0, maxiter=0, arguments=None, weights=None, maxdist=None, ischeme=None, scale=None, skip_missing=False): if ischeme is not None: warnings.deprecation('the ischeme argument is deprecated and will be removed in future') if scale is not None: @@ -1495,6 +1499,10 @@ def locate(self, geom, coords, *, tol=0, eps=0, maxiter=0, arguments=None, weigh egeom = geom.lower((), {self.space: (self.transforms.get_evaluable(_ielem), self.opposites.get_evaluable(_ielem))}, {self.space: _point}) xJ = evaluable.Tuple((egeom, evaluable.derivative(egeom, _point))).simplified arguments = dict(arguments or ()) + if skip_missing: + if weights is not None: + raise ValueError('weights and skip_missing are mutually exclusive') + missing = parallel.shzeros(len(coords), dtype=bool) with parallel.ctxrange('locating', len(coords)) as ipoints: for ipoint in ipoints: xt = coords[ipoint] # target @@ -1528,7 +1536,13 @@ def locate(self, geom, coords, *, tol=0, eps=0, maxiter=0, arguments=None, weigh points[ipoint] = p break else: - raise LocateError('failed to locate point: {}'.format(xt)) + if skip_missing: + missing[ipoint] = True + else: + raise LocateError('failed to locate point: {}'.format(xt)) + if skip_missing: + ielems = ielems[~missing] + points = points[~missing] return self._sample(ielems, points, weights) def _sample(self, ielems, coords, weights=None): @@ -2215,7 +2229,7 @@ def refined(self): axes = [axis.refined for axis in self.axes] return StructuredTopology(self.space, self.root, axes, self.nrefine+1, bnames=self._bnames) - def locate(self, geom, coords, *, tol, eps=0, weights=None, **kwargs): + def locate(self, geom, coords, *, tol, eps=0, weights=None, skip_missing=False, **kwargs): coords = numpy.asarray(coords, dtype=float) if geom.ndim == 0: geom = geom[_] @@ -2225,9 +2239,9 @@ def locate(self, geom, coords, *, tol, eps=0, weights=None, **kwargs): geom0, scale, index = self._asaffine(geom) e = self.sample('uniform', 2).eval(function.norm2(geom0 + index * scale - geom)).max() # inf-norm on non-gauss sample if e > tol: - return super().locate(geom, coords, eps=eps, tol=tol, weights=weights, **kwargs) + return super().locate(geom, coords, eps=eps, tol=tol, weights=weights, skip_missing=skip_missing, **kwargs) log.info('locate detected linear geometry: x = {} + {} xi ~{:+.1e}'.format(geom0, scale, e)) - return self._locate(geom0, scale, coords, eps=eps, weights=weights) + return self._locate(geom0, scale, coords, eps=eps, weights=weights, skip_missing=skip_missing) def _asaffine(self, geom): p0 = p1 = self @@ -2243,14 +2257,19 @@ def _asaffine(self, geom): scale = (geom1 - geom0) / numpy.array(self.shape) return geom0, scale, index - def _locate(self, geom0, scale, coords, *, eps=0, weights=None): + def _locate(self, geom0, scale, coords, *, eps=0, weights=None, skip_missing=False): mincoords, maxcoords = numpy.sort([geom0, geom0 + scale * self.shape], axis=0) - outofbounds = numpy.less(coords, mincoords - eps) | numpy.greater(coords, maxcoords + eps) - if outofbounds.any(): - raise LocateError('failed to locate {}/{} points'.format(outofbounds.any(axis=1).sum(), len(coords))) + missing = numpy.any(numpy.less(coords, mincoords - eps) | numpy.greater(coords, maxcoords + eps), axis=1) + if not skip_missing and missing.any(): + raise LocateError('failed to locate {}/{} points'.format(missing.sum(), len(coords))) xi = (coords - geom0) / scale ielem = numpy.minimum(numpy.maximum(xi.astype(int), 0), numpy.array(self.shape)-1) - return self._sample(numpy.ravel_multi_index(ielem.T, self.shape), xi - ielem, weights) + ielems = numpy.ravel_multi_index(ielem.T, self.shape) + points = xi - ielem + if skip_missing: + ielems = ielems[~missing] + points = points[~missing] + return self._sample(ielems, points, weights) def __str__(self): 'string representation' diff --git a/tests/test_function.py b/tests/test_function.py index 94ef4e309..31c05b577 100644 --- a/tests/test_function.py +++ b/tests/test_function.py @@ -1310,3 +1310,21 @@ def test_multiple(self): g = function.dotarg('v', numpy.array([3,2,1])) retvals = function.eval([f, g], v=numpy.array([4,5,6])) self.assertEqual(retvals, (4+10+18, 12+10+6)) + +class simplifications(TestCase): + + def test_multiply(self): + f = function.Argument('test', shape=(2,3), dtype=float) + self.assertIs(f * 1, f) + self.assertIs(f * 1., f) + f = function.Argument('test', shape=(2,3), dtype=int) + self.assertIs(f * 1, f) + self.assertIsNot(f * 1., f) + + def test_divide(self): + f = function.Argument('test', shape=(2,3), dtype=float) + self.assertIs(f / 1, f) + self.assertIs(f / 1., f) + f = function.Argument('test', shape=(2,3), dtype=int) + self.assertIsNot(f / 1, f) + self.assertIsNot(f / 1., f) diff --git a/tests/test_numeric.py b/tests/test_numeric.py index d18999c06..b24f725a7 100644 --- a/tests/test_numeric.py +++ b/tests/test_numeric.py @@ -226,6 +226,31 @@ def test_dtype(self): self.assertEqual(m.dtype, float) self.assertAllEqual(m, numpy.ones((1,))) +class simplex_grid(TestCase): + + def simplex_grid(self, shape, spacing): + coords = numeric.simplex_grid(shape, spacing) + self.assertEqual(coords.ndim, 2) + self.assertEqual(coords.shape[1], len(shape)) + self.assertTrue((coords > 0).all()) + self.assertTrue((coords < shape).all()) + mindist = min(numpy.linalg.norm(c1 - c2) for i, c1 in enumerate(coords) for c2 in coords[:i]) + self.assertAlmostEqual(mindist, spacing) + return coords + + def test_1d(self): + coords = self.simplex_grid([2], .8) + self.assertEqual(len(coords), 3) + self.assertAllAlmostEqual(coords[:,0], [.2,1,1.8]) + + def test_2d(self): + coords = self.simplex_grid([2,3], .8) + self.assertEqual(len(coords), 13) + + def test_3d(self): + coords = self.simplex_grid([2,3,4], .8) + self.assertEqual(len(coords), 82) + class overlapping(TestCase): def test_pairwise(self):