Welcome to Nutils’s documentation!¶
Nutils: open source numerical utilities for Python, is a collaborative programming effort aimed at the creation of a modern, general purpose programming library for Finite Element applications and related computational methods. Identifying features are a heavily object oriented design, strict separation of topology and geometry, and CAS-like function arithmetic such as found in Maple and Mathematica. Primary design goals are:
- Readability. Finite element scripts built on top of Nutils should focus on work flow and maths, unobscured by Finite Element infrastructure.
- Flexibility. The Nutils are tools; they do not enforce a strict work flow. Missing components can be added locally without loosing interoperability.
- Compatibility. Exposed objects are of native python type or allow for easy conversion to leverage third party tools.
- Speed. Nutils are self-optimizing and support parallel computation. Typical scripting inefficiencies are discouraged by design.
For latest project news and developments visit the project website at nutils.org.
Contents¶
Introduction¶
To get one thing out of the way first, note that Nutils is not your classical Finite Element program. It does not have menus, no buttons to click, nothing to make a screenshot of. To get it to do anything some programming is going to be required.
That said, let’s see what Nutils can be instead.
Design¶
Nutils is a programming library, providing components that are rich enough to handle a wide range of problems by simply linking them together. This blurs the line between classical graphical user interfaces and a programming environment, both of which serve to offer some degree of mixing and matching of available components. The former has a lower entry bar, whereas the latter offers more flexibility, the possibility to extend the toolkit with custom algorithms, and the possibility to pull in third party modules. It is our strong belief that on the edge of science where Nutils strives to be a great degree of extensibility is adamant.
For those so inclined, one of the lesser interesting possibilities this gives is to write a dedicated, Nutils powered GUI application.
What Nutils specifically does not offer are problem specific components, such as, conceivably, a “crack growth” module or “solve navier stokes” function. As a primary design principle we aim for a Nutils application to be closely readable as a high level mathematical problem description; i.e. the weak form, domain, boundary conditions, time stepping of Newton iterations, etc. It is the supporting operations like integrating over a domain or taking gradients of compound functions that are being kept out of sight as much as possible.
Quick demo¶
As a small but representative demonstration of what is involved in setting up a
problem in Nutils we solve the Laplace problem on a unit square, with
zero Dirichlet conditions on the left and bottom boundaries, unit flux at the
top and a natural boundary condition at the right. We begin by creating a
structured nelems
ⅹ nelems
Finite Element mesh using the built-in
generator:
verts = numpy.linspace( 0, 1, nelems+1 )
domain, geom = mesh.rectilinear( [verts,verts] )
Here domain
is topology representing an interconnected set of elements, and
geometry
is a mapping from the topology onto ℝ², representing it placement
in physical space. This strict separation of topological and geometric
information is key design choice in Nutils.
Proceeding to specifying the problem, we create a second order spline basis
funcsp
which doubles as trial and test space (u resp. v). We build a
matrix
by integrating laplace
= ∇v · ∇u over the domain, and a rhs
vector by integrating v over the top boundary. The Dirichlet constraints are
projected over the left and bottom boundaries to find constrained coefficients
cons
. Remaining coefficients are found by solving the system in lhs
.
Finally these are contracted with the basis to form our solution
function:
funcsp = domain.splinefunc( degree=2 )
laplace = function.outer( funcsp.grad(geom) ).sum()
matrix = domain.integrate( laplace, geometry=geom, ischeme='gauss2' )
rhs = domain.boundary['top'].integrate( funcsp, geometry=geom, ischeme='gauss1' )
cons = domain.boundary['left,bottom'].project( 0, ischeme='gauss1', geometry=geom, onto=funcsp )
lhs = matrix.solve( rhs, constrain=cons, tol=1e-8, symmetric=True )
solution = funcsp.dot(lhs)
The solution
function is a mapping from the topology onto ℝ. Sampling this
together with the geometry
generates arrays that we can use for plotting:
points, colors = domain.elem_eval( [ geom, solution ], ischeme='bezier4', separate=True )
with plot.PyPlot( 'solution', index=index ) as plt:
plt.mesh( points, colors, triangulate='bezier' )
plt.colorbar()
Wiki¶
This is a collection of technical notes.
Binary operations on Numpy/Nutils arrays¶
Tensor | Einstein | Nutils | |||
---|---|---|---|---|---|
1 | \(\mathbf{a} \in \mathbb{R}^n\) | \(\mathbf{b} \in \mathbb{R}^n\) | \(c = \mathbf{a} \cdot \mathbf{b} \in \mathbb{R}\) | \(c = a_i b_i\) | c = (a*b).sum(-1) |
2 | \(\mathbf{a} \in \mathbb{R}^n\) | \(\mathbf{b} \in \mathbb{R}^m\) | \(\mathbf{C} = \mathbf{a} \otimes \mathbf{b} \in \mathbb{R}^{n \times m}\) | \(C_{ij} = a_i b_j\) | C = a[:,_]*b[_,:] |
C = function.outer(a,b) |
|||||
3 | \(\mathbf{A} \in \mathbb{R}^{m \times n}\) | \(\mathbf{b} \in \mathbb{R}^n\) | \(\mathbf{c} = \mathbf{A}\mathbf{b} \in \mathbb{R}^{m}\) | \(c_{i} = A_{ij} b_j\) | c = (A[:,:]*b[_,:]).sum(-1) |
4 | \(\mathbf{A} \in \mathbb{R}^{m \times n}\) | \(\mathbf{B} \in \mathbb{R}^{n \times p}\) | \(\mathbf{C} = \mathbf{A} \mathbf{B} \in \mathbb{R}^{m \times p}\) | \(c_{ij} = A_{ik} B_{kj}\) | c = (A[:,:,_]*B[_,:,:]).sum(-2) |
5 | \(\mathbf{A} \in \mathbb{R}^{m \times n}\) | \(\mathbf{B} \in \mathbb{R}^{p \times n}\) | \(\mathbf{C} = \mathbf{A} \mathbf{B}^T \in \mathbb{R}^{m \times p}\) | \(C_{ij} = A_{ik} B_{jk}\) | C = (A[:,_,:]*B[_,:,:]).sum(-1) |
C = function.outer(A,B).sum(-1) |
|||||
6 | \(\mathbf{A} \in \mathbb{R}^{m \times n}\) | \(\mathbf{B} \in \mathbb{R}^{m \times n}\) | \(c = \mathbf{A} : \mathbf{B} \in \mathbb{R}\) | \(c = A_{ij} B_{ij}\) | c = (A*B).sum([-2,-1]) |
Notes:
- In the above table the summation axes are numbered backward. For example,
sum(-1)
is used to sum over the last axis of an array. Although forward numbering is possible in many situations, backward numbering is generally preferred in Nutils code. - When a summation over multiple axes is performed (#6), these axes are to be listed. In the case of single-axis summations listing is optional (for example
sum(-1)
is equivalent tosum([-1])
). The shorter notationsum(-1)
is preferred. - When the numer of dimensions of the two arguments of a binary operation mismatch, singleton axes are automatically prepended to the “shorter” argument. This property can be used to shorten notation. For example, #3 can be written as
(A*b).sum(-1)
. To avoid ambiguities, in general, such abbreviations are discouraged.
Library¶
The Nutils are separated in modules focussing on topics such as mesh generation, function manipulation, debugging, plotting, etc. They are designed to form relatively independent units, though some components such as output logging run through all. Others, such as topology and element, operate in tight connection, but are divided for reasons of scope and scale. A typical Nutils application uses methods from all modules, although, as seen above, very few modules require direct access for standard computations.
What follows is an automatically generated API reference.
Topology¶
The topology module defines the topology objects, notably the
StructuredTopology
and UnstructuredTopology
. Maintaining
strict separation of topological and geometrical information, the topology
represents a set of elements and their interconnectivity, boundaries,
refinements, subtopologies etc, but not their positioning in physical space. The
dimension of the topology represents the dimension of its elements, not that of
the the space they are embedded in.
The primary role of topologies is to form a domain for nutils.function
objects, like the geometry function and function bases for analysis, as well as
provide tools for their construction. It also offers methods for integration and
sampling, thus providing a high level interface to operations otherwise written
out in element loops. For lower level operations topologies can be used as
nutils.element
iterators.
-
class
nutils.topology.
Topology
(ndims)[source]¶ topology base class
-
elem_eval
(funcs, ischeme, separate=False, geometry=None, asfunction=False, edit=<function <lambda>>, *, arguments=None)[source]¶ element-wise evaluation
-
integrate
(funcs, ischeme='gauss', degree=None, geometry=None, force_dense=False, fcache=None, edit=<function <lambda>>, *, arguments=None)[source]¶
-
project
(fun, onto, geometry, tol=0, ischeme='gauss', degree=None, droptol=1e-12, exact_boundaries=False, constrain=None, verify=None, ptype='lsqr', precon='diag', edit=<function <lambda>>, *, arguments=None, **solverargs)[source]¶ L2 projection of function onto function space
-
-
class
nutils.topology.
WithGroupsTopology
(basetopo, vgroups={}, bgroups={}, igroups={}, pgroups={})[source]¶ item topology
-
class
nutils.topology.
StructuredLine
(root, i, j, periodic=False, bnames=None)[source]¶ structured topology
-
class
nutils.topology.
StructuredTopology
(root, axes, nrefine=0, bnames=None)[source]¶ structured topology
-
basis_spline
(degree, knotvalues=None, knotmultiplicities=None, periodic=None, removedofs=None)[source]¶ spline basis
-
refined
¶ refine non-uniformly
-
-
class
nutils.topology.
UnstructuredTopology
(ndims, elements)[source]¶ unstructured topology
-
basis_std
(degree)¶ bernstein shape functions
-
-
class
nutils.topology.
OrientedGroupsTopology
(basetopo, elems)[source]¶ unstructured topology with undirected semi-overlapping basetopology
-
class
nutils.topology.
TrimmedTopologyBoundaryItem
(btopo, trimmed, othertopo)[source]¶ trimmed topology boundary item
-
class
nutils.topology.
HierarchicalTopology
(basetopo, allelements, precise)[source]¶ collection of nested topology elments
-
class
nutils.topology.
RevolutionTopology
[source]¶ topology consisting of a single revolution element
-
class
nutils.topology.
MultipatchTopology
(patches)[source]¶ multipatch topology
-
class
Patch
(topo, verts, boundaries)¶ -
boundaries
¶ Alias for field number 2
-
topo
¶ Alias for field number 0
-
verts
¶ Alias for field number 1
-
-
class
Function¶
The function module defines the Evaluable
class and derived objects,
commonly referred to as nutils functions. They represent mappings from a
nutils.topology
onto Python space. The notabe class of Array
objects map onto the space of Numpy arrays of predefined dimension and shape.
Most functions used in nutils applicatons are of this latter type, including the
geometry and function bases for analysis.
Nutils functions are essentially postponed python functions, stored in a tree
structure of input/output dependencies. Many Array
objects have
directly recognizable numpy equivalents, such as Sin
or
Inverse
. By not evaluating directly but merely stacking operations,
complex operations can be defined prior to entering a quadrature loop, allowing
for a higher level style programming. It also allows for automatic
differentiation and code optimization.
It is important to realize that nutils functions do not map for a physical xy-domain but from a topology, where a point is characterized by the combination of an element and its local coordinate. This is a natural fit for typical finite element operations such as quadrature. Evaluation from physical coordinates is possible only via inverting of the geometry function, which is a fundamentally expensive and currently unsupported operation.
-
exception
nutils.function.
EvaluationError
(etype, evalue, evaluable, values)[source]¶ evaluation error
-
class
nutils.function.
Array
(args: tuple, shape: tuple, dtype: <function <lambda> at 0x7fd93d38de18>)[source]¶ array function
-
class
nutils.function.
Interpolate
(x: <function <lambda> at 0x7fd93d38dea0>, xp: nutils.numeric.const, fp: nutils.numeric.const, left=None, right=None)[source]¶ interpolate uniformly spaced data; stepwise for now
-
class
nutils.function.
BlockAdd
(funcs: nutils.util.frozenmultiset)[source]¶ block addition (used for DG)
-
class
nutils.function.
Sampled
(data: nutils.util.frozendict, trans=<nutils.function.SelectChain object>)[source]¶ sampled
-
class
nutils.function.
Zeros
(shape: tuple, dtype: <function <lambda> at 0x7fd93d38de18>)[source]¶ zero
-
class
nutils.function.
Guard
(fun: <function <lambda> at 0x7fd93d38dea0>)[source]¶ bar all simplifications
-
class
nutils.function.
Find
(where: <function <lambda> at 0x7fd93d38dea0>)[source]¶ indices of boolean index vector
-
class
nutils.function.
DerivativeTargetBase
(args: tuple, shape: tuple, dtype: <function <lambda> at 0x7fd93d38de18>)[source]¶ base class for derivative targets
-
class
nutils.function.
Argument
(name, shape: tuple, nderiv: int = 0)[source]¶ Array argument, to be substituted before evaluation.
The
Argument
is anArray
with a known shape, but whose values are to be defined later, before evaluation, e.g. usingreplace_arguments()
.It is possible to take the derivative of an
Array
to anArgument
:>>> from nutils import function >>> a = function.Argument('x', []) >>> b = function.Argument('y', []) >>> f = a**3 + b**2 >>> function.derivative(f, a).simplified == (3.*a**2).simplified True
Furthermore, derivatives to the local cooardinates are remembered and applied to the replacement when using
replace_arguments()
:>>> from nutils import mesh >>> domain, x = mesh.rectilinear([2,2]) >>> basis = domain.basis('spline', degree=2) >>> c = function.Argument('c', basis.shape) >>> replace_arguments(c.grad(x), dict(c=basis)) == basis.grad(x) True
Parameters:
-
class
nutils.function.
Polyval
(coeffs: <function <lambda> at 0x7fd93d38dea0>, points: <function <lambda> at 0x7fd93d38dea0>, ngrad: int = 0)[source]¶ Computes the \(k\)-dimensional array
\[\begin{split}j_0,\dots,j_{k-1} \mapsto \sum_{\substack{i_0,\dots,i_{n-1}\in\mathbb{N}\\i_0+\cdots+i_{n-1}\le d}} p_0^{i_0} \cdots p_{n-1}^{i_{n-1}} c_{j_0,\dots,j_{k-1},i_0,\dots,i_{n-1}},\end{split}\]where \(p\) are the \(n\)-dimensional local coordinates and \(c\) is the argument
coeffs
and \(d\) is the degree of the polynomial, where \(d\) is the length of the last \(n\) axes ofcoeffs
.Warning
All coefficients with a (combined) degree larger than \(d\) should be zero. Failing to do so won’t raise an
Exception
, but might give incorrect results.
-
nutils.function.
matmat
(arg0, *args)[source]¶ helper function, contracts last axis of arg0 with first axis of arg1, etc
-
nutils.function.
polyfunc
(coeffs, dofs, ndofs, transforms, *, issorted=True)[source]¶ Create an inflated
Polyval
with coefficientscoeffs
and corresponding dofsdofs
. The argumentscoeffs
,dofs
andtransforms
are assumed to have matching order. In addition, ifissorted
is true, thetransforms
argument is assumed to be sorted.
-
nutils.function.
replace_arguments
(value, arguments)[source]¶ Replace
Argument
objects invalue
.Replace
Argument
objects invalue
according to thearguments
map, taking into account derivatives to the local coordinates.Parameters: - value (
Array
) – Array to be edited. - arguments (
collections.abc.Mapping
withArray
s as values) –Argument
s replacements. The key correspond to thename
passed to anArgument
and the value is the replacement.
Returns: The edited
value
.Return type: - value (
-
class
nutils.function.
Namespace
(*, default_geometry_name='x')[source]¶ Namespace for
Array
objects supporting assignments with tensor expressions.The
Namespace
object is used to storeArray
objects.>>> from nutils import function >>> ns = function.Namespace() >>> ns.A = function.zeros([3, 3]) >>> ns.x = function.zeros([3]) >>> ns.c = 2
In addition to the assignment of
Array
objects, it is also possible to specify an array using a tensor expression string — seenutils.expression.parse()
for the syntax. All attributes defined in this namespace are available as variables in the expression. If the array defined by the expression has one or more dimensions the indices of the axes should be appended to the attribute name. Examples:>>> ns.cAx_i = 'c A_ij x_j' >>> ns.xAx = 'x_i A_ij x_j'
It is also possible to simply evaluate an expression without storing its value in the namespace by passing the expression to the method
eval_
suffixed with appropriate indices:>>> ns.eval_('2 c') Array<> >>> ns.eval_i('c A_ij x_j') Array<3> >>> ns.eval_ij('A_ij + A_ji') Array<3,3>
For zero and one dimensional expressions the following shorthand can be used:
>>> '2 c' @ ns Array<> >>> 'A_ij x_j' @ ns Array<3>
When evaluating an expression through this namespace the following functions are available:
opposite
,sin
,cos
,tan
,sinh
,cosh
,tanh
,arcsin
,arccos
,arctan2
,arctanh
,exp
,abs
,ln
,log
,log2
,log10
,sqrt
andsign
.Parameters: default_geometry_name ( str
) – The name of the default geometry. This argument is passed tonutils.expression.parse()
. Default:'x'
.-
arg_shapes
¶ types.MappingProxyType
– A readonly map of argument names and shapes.
-
default_geometry
¶ nutils.function.Array
– The default geometry, shorthand forgetattr(ns, ns.default_geometry_name)
.
-
Expression¶
This module defines the function parse()
, which parses a tensor
expression.
-
nutils.expression.
parse
(expression, variables, functions, indices, arg_shapes={}, default_geometry_name='x')[source]¶ Parse
expression
and return AST.This function parses a tensor expression with Einstein Summation Convection stored in a
str
and returns an Abstract Syntax Tree (AST). The syntax ofexpression
is as follows:Integers or decimal numbers are denoted in the usual way. Examples:
1
,1.2
,.2
. A number may not start with a zero, except when followed by a dot:0.1
is valid, but01
is not.Variables are denoted with a string of alphanumeric characters. The first character may not be a numeral. Unlike Python variables, underscores are not allowed, as they have a special meaning. If the variable is an array with one or more axes, all those axes should be labeled with a latin character, the index, and appended to the variable with an underscore. For example an array
a
with two axes can be denoted witha_ij
. Optionally, a single numeral may be used to select an item at the concerning axis. Example: ina_i0
the first axis ofa
is labeledi
and the first element of the second axis is selected. If the same index occurs twice, the trace is taken along the concerning axes. Example: the trace of the first and third axes ofb
is denoted byb_iji
. It is invalid to specify an index more than twice. The following names cannot be used as variables:n
,δ
,$
. The variable namedx
, or the value of argumentdefault_geometry_name
, has a special meaning, detailed below.A term, the product of two or more arrays or scalars, is denoted by space-separated variables, constants or compound expressions. Example:
a b c
denotes the product of the scalarsa
,b
andc
. A term may start with a number, but a number is not allowed in other parts of the term. Example:2 a
denotes two timesa
;2 2 a
and2 a 2`
are invalid. When two arrays in a term have the same index, this index is summed. Example:a_i b_i
denotes the inner product ofa
andb
andA_ij b_j`
a matrix vector product. It is not allowed to use an index more than twice in a term.The operator
/
denotes a fraction. Example: ina b / c d
a b
is the numerator andc d
the denominator. Both the numerator and the denominator may start with a number. Example:2 a / 3 b
. The denominator must be a scalar. Example:2 / a_i b_i
is valid, but2 a_i / b_i
is not.Warning
This syntax is different from the Python syntax. In Python
a*b / c*d
is mathematically equivalent toa*b*d/c
.The operators
+
and-
denote add and subtract. Both operators should be surrounded by whitespace, e.g.a + b
. Both operands should have the same shape. Example:a_ij + b_i c_j
is a valid, provided that the lengths of the axes with the same indices match, buta_ij + b_i
is invalid. At the beginning of an expression or a compound-
may be used to negate the following term. Example: in-a b + c
the terma b
is negated before addingc
. It is not allowed to negate other terms:a + -b
is invalid, so isa -b
.An expression surrounded by parentheses is a compound expression and can be used as single entity in a term. Example:
(a_i + b_i) c_i
denotes the inner product ofa_i + b_i
withc_i
.Exponentiation is denoted by a
^
, where the left and right operands should be a number, variable or compound expression and the right operand should be a scalar. Example:a^2
denotes the square ofa
,a^-2
denotesa
to the power-2
anda^(1 / 2)
the square root ofa
.An argument is denoted by a name — following the same rules as a variable name — prefixed with a question mark. An argument is a scalar or array with a yet unknown value. Example:
basis_i ?coeffs_i
denotes the inner product of a basis with unknown coefficient vector?coeffs
. If possible the shape of the argument is deduced from the expression. In the previous example the shape of?coeffs
is equal to the shape ofbasis
. If the shape cannot be deduced from the expression the shape should be defined manually (seeparse()
). Arguments and variables live in separate namespaces:?x
andx
are different entities.An argument may be substituted by appending without whitespace
(arg = value)
to a variable of compound expression, wherearg
is an argument andvalue
the substitution. The substitution applies to the variable of compound expression only. The value may be an expression. Example:2 ?x(x = 3 + y)
is equivalent to2 (3 + y)
and2 ?x(x=y) + 3
is equivalent to2 (y) + 3
. It is possible to apply multiple substitutions. Example:(?x + ?y)(x = 1, y = )2
is equivalent to1 + 2
.The gradient of a variable to the default geometry — the default geometry is variable
x
unless overriden by the argumentdefault_geometry_name
— is denoted by an underscore, a comma and an index. If the variable is an array with more than one axis, the underscore is omitted. Example:a_,i
denotes the gradient of the scalara
to the geometry andb_i,j
the gradient of vectorb
. The gradient of a compound expression is denoted by an underscore, a comma and an index. Example:(a_i + b_j)_,k
denotes the gradient ofa_i + b_j
. The usual summation rules apply and it is allowed to use a numeral as index. The surface gradient is denoted with a semicolon instead of a comma, but follows the same rules as the gradient otherwise. Example:a_i;j
is the sufrace gradient ofa_i
to the geometry. It is also possible to take the gradient to another geometry by appending the name of the geometry, which should exist as a variable, and an underscore directly after the comma of semicolon. Example:a_i,altgeom_j
denotes the gradient ofa_i
toaltgeom
and the gradient axis has indexj
. Futhermore, it is possible to take the derivative to an argument by adding the argument with appropriate indices after the comma. Example:(?x^2)_,?x
denotes the derivative of?x^2
to?x
, which is equivalent to2 ?x
, and(?y_i ?y_i),?y_j
is the derivative of?y_i ?y_i
to?y_j
, which is equivalent to2 ?y_j
.The normal of the default geometry is denoted by
n_i
, where the indexi
may be replaced with an index of choice. The normal with respect to different geometry is denoted by appending an underscore with the name of the geometry right aftern
. Example:n_altgeom_j
is the normal with respect to geometryaltgeom
.A dirac is denoted by
δ
or$
and takes two indices. The shape of the dirac is deduced from the expression. Example: letA
be a square matrix with three rows and columns, thenδ_ij
in(A_ij - λ δ_ij) x_j
has three rows and columns as well.An expression surrounded by square brackets or curly braces denotes the jump or mean, respectively, of the enclosed expression. Example:
[ a_i ]
denotes the jump ofa_i
and{ a_i + b_i }
denotes the mean ofa_i + b_i
.A function call is denoted by a name — following the same rules as for a variable name — directly followed by the left parenthesis
(
, without a space. The arguments to the function are separated by a comma and at least one space. The function is applied pointwise to the arguments and all arguments should have the same shape. Example:f(x_i, y_i)
.denotes the call to functionf
with argumentsx_i
andy_i
. Functions and variables share a namespace: defining a variable with the same name as a function renders the function inaccessible.A stack of two or more arrays along an axis is denoted by a
<
followed by comma and space separated arrays followed by>
and an index. If an argument does not have an axis with the specified stack index, the argument is expanded with an axis of length one. Beside the stack axis, all arguments should have the same shape. Example:<1, x_i>_i
, withx
a vector of length three, creates an array with components1
,x_0
,x_1
,x_2
.
Parameters: - expression (
str
) – The expression to parse. Seeexpression
for the expression syntax. - variables (
dict
ofstr
andnutils.function.Array
pairs) – Adict
of variable names and array pairs. All variables used in theexpression
should exist invariables
. - functions (
dict
ofstr
andint
pairs) – Adict
of function names and number of arguments pairs. All functions used in theexpression
should exist infunctions
. - indices (
str
) – The indices used for aligning the resulting array. For example, letexpression
be'a_ij'
. Ifindices
is'ij'
, then the returned array is simplyvariables['a']
, but ifindices
is'ji'
the transpose ofvariables['a']
is returned. All indices of theexpression
should be listed precisely once. - arg_shapes (
dict
ofstr
andtuple
orint
s pairs) – Adict
of argument names and shapes. Ifexpression
contains an argument not present inarg_shapes
the shape will be decuded from the expression and added to a copy ofarg_shapes
. - default_geometry_name (
str
) – The name of the default geometry variable. When computing a gradient or the normal, e.g.'f_,i'
or'n_i'
, this variable is used as the geometry, unless the geometry is explicitly mentioned in the expression. Default:'x'
.
Returns: ast (
tuple
) – The parsedexpression
as an abstract syntax tree (AST). The AST is atuple
of an opcode and arguments. The special opcodeNone
indicates that the single argument is used verbatim. All other opcodes have AST as arguments. The following opcodes exist:(None, const) ('group', group) ('arg', name, *shape) ('substitute', array, arg, value) ('call', func, arg) ('eye', length) ('normal', geom) ('getitem', array, dim, index) ('trace', array, n1, n2) ('sum', array, axis) ('concatenate', *args) ('grad', array, geom) ('surfgrad', array, geom) ('derivative', func, target) ('append_axis', array, length) ('transpose', array, trans) ('jump', array) ('mean', array) ('neg', array) ('add', left, right) ('sub', left, right) ('mul', left, right) ('truediv', left, right) ('pow', left, right)
arg_shapes (
dict
ofstr
andtuple
ofint
s pairs) – A copy ofarg_shapes
updated with shapes of arguments present in thisexpression
.
Core¶
The core module provides a collection of low level constructs that have no dependencies on other nutils modules. Primarily for internal use.
Config¶
This module holds the Nutils global configuration, stored as (immutable)
attributes. To inspect the current configuration, use print()
or
vars()
on this module. The configuration can be changed temporarily by
calling this module with the new settings passed as keyword arguments and
entering the returned context. The old settings are restored as soon as the
context is exited. Example:
>>> from nutils import config
>>> config.verbose
4
>>> with config(verbose=2, nprocs=4):
... # The configuration has been updated.
... config.verbose
2
>>> # Exiting the context reverts the changes:
>>> config.verbose
4
Note
The default entry point for Nutils scripts nutils.cli.run()
(and
nutils.cli.choose()
) will read user configuration from disk.
Important
The configuration is not thread-safe: changing the configuration inside a thread changes the process wide configuration.
Element¶
The element module defines reference elements such as the QuadElement
and TriangularElement
, but also more exotic objects like the
TrimmedElement
. A set of (interconnected) elements together form a
nutils.topology
. Elements have edges and children (for refinement), which
are in turn elements and map onto self by an affine transformation. They also
have a well defined reference coordinate system, and provide pointsets for
purposes of integration and sampling.
-
class
nutils.element.
Element
(reference, trans, opptrans=None, oriented=False)[source]¶ element class
-
class
nutils.element.
RevolutionReference
[source]¶ modify gauss integration to always return a single point
-
class
nutils.element.
TriangleReference
[source]¶ 2D simplex
-
getischeme_gauss
(degree)[source]¶ get integration scheme http://www.cs.rpi.edu/~flaherje/pdf/fea6.pdf
-
-
class
nutils.element.
TetrahedronReference
[source]¶ 3D simplex
-
getischeme_gauss
(degree)[source]¶ get integration scheme http://www.cs.rpi.edu/~flaherje/pdf/fea6.pdf
-
-
class
nutils.element.
NeighborhoodTensorReference
(ref1, ref2, neighborhood, transf)[source]¶ product reference element
Log¶
The log module provides print methods debug
, info
, user
,
warning
, and error
, in increasing order of priority. Output is sent to
stdout as well as to an html formatted log file if so configured.
-
class
nutils.log.
Log
[source]¶ Base class for log objects. A subclass should define a
context()
method that returns a context manager which adds a contextual layer and awrite()
method.
-
class
nutils.log.
ContextLog
[source]¶ Base class for loggers that keep track of the current list of contexts.
The base class implements
context()
which keeps the attribute_context
up-to-date.
-
class
nutils.log.
ContextTreeLog
[source]¶ Base class for loggers that display contexts as a tree.
-
_print_push_context
(title)[source]¶ Push a context to the log.
This method is called just before the first item of this context is added to the log. If no items are added to the log within this context or children of this context this method nor
_print_pop_context()
will be called.Note
This function is abstract.
-
_print_pop_context
()[source]¶ Pop a context from the log.
This method is called whenever a context is exited, but only if
_print_push_context()
has been called before for the same context.Note
This function is abstract.
-
write
(level, text)[source]¶ Write
text
with log levellevel
to the log.This method makes sure the current context is printed and calls
_print_item()
.
-
-
class
nutils.log.
RichOutputLog
(stream=None, *, progressinterval=None)[source]¶ Output rich (colored,unicode) text to stream.
-
class
nutils.log.
HtmlInsertAnchor
[source]¶ Mix-in class for HTML-based loggers that inserts anchor tags for paths.
-
class
nutils.log.
HtmlLog
(file, *, title='nutils', scriptname=None, funcname=None, funcargs=None)[source]¶ Output html nested lists.
Matrix¶
The matrix module defines a number of 2D matrix objects, notably the
ScipyMatrix()
and NumpyMatrix()
. Matrix objects support basic
addition and subtraction operations and provide a consistent insterface for
solving linear systems. Matrices can be converted to numpy arrays via
toarray
or scipy matrices via toscipy
.
Mesh¶
The mesh module provides mesh generators: methods that return a topology and an
accompanying geometry function. Meshes can either be generated on the fly, e.g.
rectilinear()
, or read from external an externally prepared file,
gmsh()
, and converted to nutils format. Note that no mesh writers are
provided at this point; output is handled by the nutils.plot
module.
-
nutils.mesh.
multipatch
(patches, nelems, patchverts=None, name='multipatch')[source]¶ multipatch rectilinear mesh generator
Generator for a
MultipatchTopology
and geometry. TheMultipatchTopology
consists of a set patches, where each patch is aStructuredTopology
and all patches have the same number of dimensions.The
patches
argument, anumpy.ndarray
-like with shape(npatches, 2*ndims)
or(npatches,)+(2,)*ndims
, defines the connectivity by labelling the patch vertices. For example, three one-dimensional patches can be connected at one edge by:# connectivity: 3 # │ # 1──0──2 patches=[[0,1], [0,2], [0,3]]
Or two two-dimensional patches along an edge by:
# connectivity: 3──4──5 # │ │ │ # 0──1──2 patches=[[[0,3],[1,4]], [[1,4],[2,5]]]
The geometry is specified by the
patchverts
argument: anumpy.ndarray
-like with shape(nverts,ngeomdims)
specifying for each vertex a coordinate. Note that the dimension of the geometry may be higher than the dimension of the patches. The created geometry is a patch-wise linear interpolation of the vertex coordinates. If thepatchverts
argument is omitted the geometry describes a unit hypercube per patch.The
nelems
argument is either anint
defining the number of elements per patch per dimension, or adict
with edges (a pair of vertex numbers) as keys and the number of elements (int
) as values, with keyNone
specifying the default number of elements. Example:# connectivity: 3─────4─────5 # │ 4x3 │ 8x3 │ # 0─────1─────2 patches=[[[0,3],[1,4]], [[1,4],[2,5]]] nelems={None: 4, (1,2): 8, (4,5): 8, (0,3): 3, (1,4): 3, (2,5): 3}
Since the patches are structured topologies, the number of elements per patch per dimension should be unambiguous. In above example specifying
nelems={None: 4, (1,2): 8}
will raise an exception because the patch on the right has 8 elements along edge(1,2)
and 4 along(4,5)
.Example
An L-shaped domain can be generated by:
# connectivity: 2──5 # │ | # 1──4─────7 y # │ │ │ │ # 0──3─────6 └──x domain, geom = mesh.multipatch( patches=[[0,1,3,4], [1,2,4,5], [3,4,6,7]], patchverts=[[0,0], [0,1], [0,2], [1,0], [1,1], [1,2], [3,0], [3,1]], nelems={None: 4, (3,6): 8, (4,7): 8})
The number of elements is chosen such that all elements in the domain have the same size.
A topology and geometry describing the surface of a sphere can be generated by creating a multipatch cube surface and inflating the cube to a sphere:
# connectivity: 3────7 # ╱│ ╱│ # 2────6 │ y # │ │ │ │ │ # │ 1──│─5 │ z # │╱ │╱ │╱ # 0────4 *────x topo, cube = multipatch( patches=[ # The order of the vertices is chosen such that normals point outward. [2,3,0,1], [4,5,6,7], [4,6,0,2], [1,3,5,7], [1,5,0,4], [2,6,3,7], ], patchverts=tuple(itertools.product(*([[-1,1]]*3))), nelems=10, ) sphere = cube / function.sqrt((cube**2).sum(0))
Parameters: - patches – A
numpy.ndarray
with shape sequence of patches with each patch being a list of vertex indices. - patchverts – A sequence of coordinates of the vertices.
- nelems – Either an
int
specifying the number of elements per patch per dimension, or adict
with edges (a pair of vertex numbers) as keys and the number of elements (int
) as values, with keyNone
specifying the default number of elements.
Returns: nutils.topology.MultipatchTopology
– The multipatch topology.nutils.function.Array
– The geometry defined by thepatchverts
or a unit hypercube per patch ifpatchverts
is not specified.
- patches – A
-
nutils.mesh.
gmsh
(fname, name=None)[source]¶ Gmsh parser
Parser for Gmsh files in .msh format. Only files with physical groups are supported. See the Gmsh manual for details.
Parameters: Returns: Topology of parsed Gmsh file geom (
nutils.function.Array
): Isoparametric mapReturn type: topo (
nutils.topology.Topology
)
Numeric¶
The numeric module provides methods that are lacking from the numpy module.
-
nutils.numeric.
dot
(A, B, axis=-1)[source]¶ Transform axis of A by contraction with first axis of B and inserting remaining axes. Note: with default axis=-1 this leads to multiplication of vectors and matrices following linear algebra conventions.
-
nutils.numeric.
diagonalize
(arg, axis=-1, newaxis=-1)[source]¶ insert newaxis, place axis on diagonal of axis and newaxis
-
nutils.numeric.
eig
(A)[source]¶ If A has repeated eigenvalues, numpy.linalg.eig sometimes fails to produce the complete eigenbasis. This function aims to fix that by identifying the problem and completing the basis where necessary.
-
nutils.numeric.
ix
(args)[source]¶ version of
numpy.ix_()
that allows for scalars
Parallel¶
The parallel module provides tools aimed at parallel computing. At this point
all parallel solutions use the fork
system call and are supported on limited
platforms, notably excluding Windows. On unsupported platforms parallel features
will disable and a warning is printed.
-
nutils.parallel.
shempty
(shape, dtype=<class 'float'>)[source]¶ create uninitialized array in shared memory
-
nutils.parallel.
shzeros
(shape, dtype=<class 'float'>)[source]¶ create zero-initialized array in shared memory
-
nutils.parallel.
pariter
(iterable, nprocs)[source]¶ iterate in parallel
Fork into
nprocs
subprocesses, then yield items from iterable such that all processes receive a nonoverlapping subset of the total. It is up to the user to prepare shared memory and/or locks for inter-process communication. The following creates a data vector containing the first four quadratics:data = shzeros(shape=[4], dtype=int) for i in pariter(range(4), 2): data[i] = i**2 data
As a safety measure nested pariters are blocked by setting the global
procid
variable; all secundary pariters will be treated like normal serial iterators.Parameters: - iterable (iterable) – The collection of items to be distributed over processors
- nprocs (int) – Maximum number of processers to use
Yields: Items from iterable, distributed over at most nprocs processors.
-
nutils.parallel.
parmap
(func, iterable, nprocs, shape=(), dtype=<class 'float'>)[source]¶ parallel equivalent to builtin map function
Produces an array of
func(item)
values for all items initerable
. Because of shared memory restrictionsfunc
must yield numpy arrays of predetermined shape and type.Parameters: Returns: Return type: Array of shape
len(iterable),+shape
and dtypedtype
Util¶
The util module provides a collection of general purpose methods. Most
importantly it provides the run()
method which is the preferred entry
point of a nutils application, taking care of command line parsing, output dir
creation and initiation of a log file.
-
nutils.util.
single_or_multiple
(f)[source]¶ Method wrapper, converts first positional argument to tuple: tuples/lists are passed on as tuples, other objects are turned into tuple singleton. Return values should match the length of the argument list, and are unpacked if the original argument was not a tuple/list.
>>> class Test: ... @single_or_multiple ... def square(self, args): ... return [v**2 for v in args] ... >>> T = Test() >>> T.square(2) 4 >>> T.square([2,3]) [4, 9]
Parameters: f (method) – Method that expects a tuple as first positional argument, and that returns a list/tuple of the same length. Returns: Return type: Wrapped method.
Plot¶
The plot module aims to provide a consistent interface to various plotting backends. At this point matplotlib and vtk are supported.
-
class
nutils.plot.
BasePlot
(name=None, ndigits=0, index=None)[source]¶ base class for plotting objects
-
class
nutils.plot.
PyPlot
(name=None, imgtype=None, ndigits=3, index=None, **kwargs)[source]¶ matplotlib figure
-
mesh
(points, values=None, edgecolors='k', edgewidth=0.1, mergetol=0, setxylim=True, aspect='equal', tight=True, **kwargs)[source]¶ plot elemtwise mesh
-
slope_marker
(x, y, slope=None, width=0.2, xoffset=0, yoffset=0.2, color='0.5')[source]¶ slope marker
-
slope_triangle
(x, y, fillcolor='0.9', edgecolor='k', xoffset=0, yoffset=0.1, slopefmt='{0:.1f}')[source]¶ Draw slope triangle for supplied y(x) - x, y: coordinates - xoffset, yoffset: distance graph & triangle (points) - fillcolor, edgecolor: triangle style - slopefmt: format string for slope number
-
-
class
nutils.plot.
PyPlotVideo
(name, videotype=None, clearfigure=True, framerate=24)[source]¶ matplotlib based video generator
Video generator based on matplotlib figures. Follows the same syntax as PyPlot.
Parameters: - clearfigure (bool, default: True) – If True clears the matplotlib figure after writing each frame.
- framerate (int, float, default: 24) – Framerate in frames per second of the generated video.
- videotype (str, default: ‘webm’ unless overriden by property
videotype
) – Video type of the generated video. Note that not every video type supports playback before the video has been finalized, i.e. beforeclose
has been called. - properties (Nutils) –
- ----------------- –
- videotype – see parameter with the same name
- videoencoder (str, default: 'ffmpeg') – Name or path of the video encoder. The video encoder should take the same arguments as ‘ffmpeg’.
Examples
Using a
with
-statement:video = PyPlotVideo('video') for timestep in timesteps: ... with video: video.plot(...) video.title('frame {:04d}'.format(video.frame)) video.close()
Using
saveframe
:video = PyPlotVideo('video') for timestep in timesteps: ... video.plot(...) video.title('frame {:04d}'.format(video.frame)) video.saveframe() video.close()
Cache¶
The cache module.
-
nutils.cache.
replace
(func)[source]¶ decorator for deep object replacement
Generates a deep replacement method for Immutable objects based on a callable that is applied (recursively) on individual constructor arguments.
Parameters: func (callable which maps (obj, ..) onto replaced_obj) – Returns: The method that searches the object to perform the replacements. Return type: callable
Cli¶
The cli (command line interface) module provides the cli.run function that can be used set up properties, initiate an output environment, and execute a python function based arguments specified on the command line.
-
nutils.cli.
run
(func, *, skip=1, loaduserconfig=True)[source]¶ parse command line arguments and call function
Solver¶
The solver module defines the Integral
class, which represents an
unevaluated integral. This is useful for fully automated solution procedures
such as Newton, that require functional derivatives of an entire functional.
To demonstrate this consider the following setup:
>>> from nutils import mesh, function, solver
>>> ns = function.Namespace()
>>> domain, ns.x = mesh.rectilinear([4,4])
>>> ns.basis = domain.basis('spline', degree=2)
>>> cons = domain.boundary['left,top'].project(0, onto=ns.basis, geometry=ns.x, ischeme='gauss4')
project > constrained 11/36 dofs, error 0.00e+00/area
>>> ns.u = 'basis_n ?lhs_n'
Function u
represents an element from the discrete space but cannot not
evaluated yet as we did not yet establish values for ?lhs
. It can,
however, be used to construct a residual functional res
. Aiming to solve
the Poisson problem u_,kk = f
we define the residual functional res = v,k
u,k + v f
and solve for res == 0
using solve_linear
:
>>> res = domain.integral('basis_n,i u_,i + basis_n' @ ns, geometry=ns.x, degree=2)
>>> lhs = solver.solve_linear('lhs', residual=res, constrain=cons)
solve > solving system using sparse direct solver
The coefficients lhs
represent the solution to the Poisson problem.
In addition to solve_linear
the solver module defines newton
and
pseudotime
for solving nonlinear problems, as well as impliciteuler
for
time dependent problems.
-
nutils.solver.
solve_linear
(target, residual, constrain=None, *, arguments=None, **solveargs)[source]¶ solve linear problem
Parameters: - target (
str
) – Name of the target: anutils.function.Argument
inresidual
. - residual (Integral) – Residual integral, depends on
target
- constrain (float vector) – Defines the fixed entries of the coefficient vector
- arguments (
collections.abc.Mapping
) – Defines the values fornutils.function.Argument
objects in residual. Thetarget
should not be present inarguments
. Optional.
Returns: Array of
target
values for whichresidual == 0
Return type: vector
- target (
-
nutils.solver.
solve
(gen_lhs_resnorm, tol=1e-10, maxiter=None)[source]¶ execute nonlinear solver
Iterates over nonlinear solver until tolerance is reached. Example:
lhs = solve(newton(target, residual), tol=1e-5)
Parameters: Returns: Coefficient vector that corresponds to a smaller than
tol
residual.Return type: vector
-
nutils.solver.
withsolve
(f)[source]¶ add a .solve method to (lhs,resnorm) iterators
Introduces the convenient form:
newton(target, residual).solve(tol)
Shorthand for:
solve(newton(target, residual), tol)
-
class
nutils.solver.
newton
(*args, **kwargs)[source]¶ iteratively solve nonlinear problem by gradient descent
Generates targets such that residual approaches 0 using Newton procedure with line search based on a residual integral. Suitable to be used inside
solve
.An optimal relaxation value is computed based on the following cubic assumption:
| res(lhs + r * dlhs) |^2 = A + B * r + C * r^2 + D * r^3
where
A
,B
,C
andD
are determined based on the current and updated residual and tangent.Parameters: - target (
str
) – Name of the target: anutils.function.Argument
inresidual
. - residual (Integral) –
- lhs0 (vector) – Coefficient vector, starting point of the iterative procedure.
- constrain (boolean or float vector) – Equal length to
lhs0
, masks the free vector entries asFalse
(boolean) or NaN (float). In the remaining positions the values oflhs0
are returned unchanged (boolean) or overruled by the values in constrain (float). - nrelax (int) – Maximum number of relaxation steps before proceding with the updated coefficient vector (by default unlimited).
- minrelax (float) – Lower bound for the relaxation value, to force re-evaluating the functional in situation where the parabolic assumption would otherwise result in unreasonably small steps.
- maxrelax (float) – Relaxation value below which relaxation continues, unless
nrelax
is reached; should be a value less than or equal to 1. - rebound (float) – Factor by which the relaxation value grows after every update until it reaches unity.
- arguments (
collections.abc.Mapping
) – Defines the values fornutils.function.Argument
objects in residual. Thetarget
should not be present inarguments
. Optional.
Yields: vector – Coefficient vector that approximates residual==0 with increasing accuracy
- target (
-
class
nutils.solver.
pseudotime
(*args, **kwargs)[source]¶ iteratively solve nonlinear problem by pseudo time stepping
Generates targets such that residual approaches 0 using hybrid of Newton and time stepping. Requires an inertia term and initial timestep. Suitable to be used inside
solve
.Parameters: - target (
str
) – Name of the target: anutils.function.Argument
inresidual
. - residual (Integral) –
- inertia (Integral) –
- timestep (float) – Initial time step, will scale up as residual decreases
- lhs0 (vector) – Coefficient vector, starting point of the iterative procedure.
- constrain (boolean or float vector) – Equal length to
lhs0
, masks the free vector entries asFalse
(boolean) or NaN (float). In the remaining positions the values oflhs0
are returned unchanged (boolean) or overruled by the values in constrain (float). - arguments (
collections.abc.Mapping
) – Defines the values fornutils.function.Argument
objects in residual. Thetarget
should not be present inarguments
. Optional.
Yields: vector, float – Tuple of coefficient vector and residual norm
- target (
-
nutils.solver.
thetamethod
(target, residual, inertia, timestep, lhs0, theta, target0='_thetamethod_target0', constrain=None, newtontol=1e-10, *, arguments=None, **newtonargs)[source]¶ solve time dependent problem using the theta method
Parameters: - target (
str
) – Name of the target: anutils.function.Argument
inresidual
. - residual (Integral) –
- inertia (Integral) –
- timestep (float) – Initial time step, will scale up as residual decreases
- lhs0 (vector) – Coefficient vector, starting point of the iterative procedure.
- theta (float) – Theta value (theta=1 for implicit Euler, theta=0.5 for Crank-Nicolson)
- residual0 (Integral) – Optional additional residual component evaluated in previous timestep
- constrain (boolean or float vector) – Equal length to
lhs0
, masks the free vector entries asFalse
(boolean) or NaN (float). In the remaining positions the values oflhs0
are returned unchanged (boolean) or overruled by the values in constrain (float). - newtontol (float) – Residual tolerance of individual timesteps
- arguments (
collections.abc.Mapping
) – Defines the values fornutils.function.Argument
objects in residual. Thetarget
should not be present inarguments
. Optional.
Yields: vector – Coefficient vector for all timesteps after the initial condition.
- target (
-
nutils.solver.
optimize
(target, functional, droptol=None, lhs0=None, constrain=None, newtontol=None, *, arguments=None)[source]¶ find the minimizer of a given functional
Parameters: - target (
str
) – Name of the target: anutils.function.Argument
inresidual
. - functional (scalar Integral) – The functional the should be minimized by varying target
- droptol (
float
) – Threshold for leaving entries in the return value at NaN if they do not contribute to the value of the functional. - lhs0 (vector) – Coefficient vector, starting point of the iterative procedure (if applicable).
- constrain (boolean or float vector) – Equal length to
lhs0
, masks the free vector entries asFalse
(boolean) or NaN (float). In the remaining positions the values oflhs0
are returned unchanged (boolean) or overruled by the values in constrain (float). - newtontol (float) – Residual tolerance of Newton procedure (if applicable)
Yields: vector – Coefficient vector corresponding to the functional optimum
- target (
Transform¶
The transform module.