Functions

Function class objects provide the building blocks for mathematical expression within Underworld2. The primary aim of this class is to enable a natural description of mathematics through the Python syntax so that users may quickly and accurately prototype model behaviour.

Functions are used extensively across the Underworld2 API and provide a unified interface to Underworld2 discrete objects.

Overview:

  1. A simple example.
  2. Usage basics.
  3. Module overview.
  4. The evaluate() method.
  5. The input function.
  6. Branching functions.

Keywords: functions, swarms, meshvariables, materials

A Simple Example

Let us define a function which we might use as a variable heat conductivity for a thermal problem. It will take the following temperature dependent form:

$$ k(\mathbf{x}) = 5 +8\exp({5T(\mathbf{x})}) $$
In [1]:
import underworld as uw
from underworld import function as fn
import glucifer

# first create a mesh and variable
mesh = uw.mesh.FeMesh_Cartesian(elementRes=(64,64))
tempVar = mesh.add_variable(1)
# init the temp variable 
for index,coord in enumerate(mesh.data):
    tempVar.data[index] = coord[1]
    
# and now define the function.
fn_k = 5. + 8.*fn.math.exp(5.*tempVar)            # a

# vis
fig = glucifer.Figure()
fig.append(glucifer.objects.Surface(mesh,fn_k))   # b
fig.show()                                        # c

Let's deconstruct the Function defined at step (a) above:

fn_k = 5.                              # 1
          + 8.*                        # 2
               fn.math.exp(            # 3
                            5.*temp )  # 4

Things to note at the positions above:

  1. You can directly use Python native numerical objects. Under the hood, the native Python float object created here '5.' will be automatically converted to an Underworld2 Constant type function. Note that arithmetic operations only currently support float type objects, and an exception will be thrown for other types. For this reason, you often need to be careful where you use Python natives (for instance, using '5' here instead of '5.' would result in an error).
  2. Again, the native '8.' will be automatically converted. The addition operator here will be automatically converted to an Underworld Addition operation through operator overloading. Likewise for the multiplication operation.
  3. Note that for an exponential function, we need to use the Underworld provided fn.math.exp function, not the Python math module exp function.
  4. Here the argument (5.*temp) is itself a Function, and Function compounding applies. Importantly, note that the MeshVariable is used directly in the arithmetic, and this is possible because it is also a Function class object (more on this soon).

At step (b), we provide the function to the visualisation object (Surface). Function objects are expected in many places across the Underworld2 API. At this stage, tests are performed to ensure that the provided function is valid and compatible.

Finally, at step (c), the actually function evaluation occurs. This will generally be the most compuationally expensive phase, with the provided function being evaluated numerous times.

Usage basics

Underworld data objects

As seen in the example above, Underworld data objects (MeshVariable and SwarmVariable types) may be used directly within functions, as they are indeed themselves Function objects (via Python multiple inheritance). The Function class provides a uniform interface to these objects and is (largely) agnostic to the underlying data discretisation, instead providing mechanisms for evaluation at arbitrary coordinates. Note that true arbitrary coordinate evaluation is not possible for SwarmVariable objects, as they are purely discrete and do not (currently) have supporting interpolation functions. The special case inputs provided by objects of the FunctionInput class may be used for SwarmVariable objects (more on this soon).

The following simple example demonstrates querying a MeshVariable object.

In [2]:
# create a new mesh variable
sinvar = mesh.add_variable(1)

# initialise with a sine perturbation
import numpy as np
sinvar.data[:,0] = np.sin( np.pi*mesh.data[:,0] )

# use the `evaluate()` method to perform query.. more on this below 
result = sinvar.evaluate( (0.25,0.) ) 
print("Evaluation result: {}".format(result))

Evaluation result: [[ 0.70710678]]
In [3]:
# Let's define a test which we will use to ensure we're getting the correct numbers:
def test( result, expected ):
    if not np.allclose(result, expected, rtol=2e-2, atol=2e-2 ): raise RuntimeError("Error! Expected result was not obtained.")
# use test
import math
test(result,math.sqrt(2.)/2.)

Elementary algebraic operations

Use the Python equivalents (+,-,*,/) directly with your Function objects! Function objects are operator overloaded to facilitate this. Note however that only functions which return floating point type values are compatible with elementary operations currently, and an exception will be raised otherwise.

Let's have a play with some mesh variables initialised to constants.

In [4]:
# create some more mesh variables and initialise
two_var = mesh.add_variable(1)
six_var = mesh.add_variable(1)
two_var.data[:] = 2.
six_var.data[:] = 6.

# create some functions via the Python operators
fn_plus  = two_var + six_var
fn_minus = two_var - six_var
fn_div   = two_var / six_var
fn_times = two_var * six_var

# check results.. evaluate anywhere as our mesh variables are constant
coord = (0.1234, 0.5678)
resultplus  =  fn_plus.evaluate( coord )
resultminus = fn_minus.evaluate( coord )
resultdiv   =   fn_div.evaluate( coord )
resulttimes = fn_times.evaluate( coord )

print("Addition result      : {}".format(resultplus))
print("Subtraction result   : {}".format(resultminus))
print("Division result      : {}".format(resultdiv))
print("Multiplication result: {}".format(resulttimes))

# run tests
test(resultplus,  8.)
test(resultminus,-4.)
test(resultdiv,   1./3.)
test(resulttimes, 12.)

Addition result      : [[ 8.]]
Subtraction result   : [[-4.]]
Division result      : [[ 0.33333333]]
Multiplication result: [[ 12.]]

Convenience conversions

Python elementary types (int/floats/etc) may be used directly with Underworld Function objects in algebraic operations. Likewise, Python tuples (or lists) of Underworld Functions are automatically converted into Function objects which return vector results composed of the tuple/list entries. It is often important to remember that the elementary algebraic operations are only compatible with float type objects, and therefore you will often need to write '5.' (which will be converted to a double precision float object) instead of '5' (which will be converted to a integer object).

In [5]:
# create a new mesh variable and init
cosvar = mesh.add_variable(1)
cosvar.data[:,0] = np.cos( np.pi*mesh.data[:,0] )

# create function.. note that `2.` and the `1.` are automatically converted
sin2 = 2.*sinvar*cosvar + 1.

# evaluate somewhere and then test
import random
coord = ( random.random(), 0. )
result = sin2.evaluate(coord)
expected = math.sin(2.*math.pi*coord[0]) + 1.  # via double angle formula
test(result, expected)

# also let's create a vector function on the fly.
# first create a vector mesh variable
vecvar = mesh.add_variable(2)
vecvar.data[:] = (1.,1.)

# now the function definition, evaluation, and test
fn_vec = vecvar + (sinvar, cosvar)
result = fn_vec.evaluate(coord)
expected = (math.sin(math.pi*coord[0]) + 1., math.cos(math.pi*coord[0]) + 1.)
test(result,expected)

Note that conversions can only occur automatically where the Python object (tuple or primary math object) comes in contact with an Underworld function, whereby conversion occurs by virtue of the object overloading. However, occasionally we may wish to explicitly perform conversions, generally where we wish to use the evaluate() method or when we wish to utilise overloading. The convert() static method on the Function class provides this functionality:

In [6]:
vec_as_py_tuple = (sinvar, cosvar)
print("`vec_as_py_tuple` type is: {}".format(type(vec_as_py_tuple)))

# this will not work!
# vec_as_py_tuple.evaluate()  

vec_as_uw_fn = fn.Function.convert(vec_as_py_tuple)
print("`vec_as_uw_fn` type is: {}".format(type(vec_as_uw_fn)))

# this is better
print("evaluate: {}".format(vec_as_uw_fn.evaluate(coord)))

`vec_as_py_tuple` type is: <type 'tuple'>
`vec_as_uw_fn` type is: <class 'underworld.function._function.add'>
evaluate: [[ 0.9401396   0.34053079]]

Basic mathematical functions

Basic functions (such as sin() and exp()) are provided by the underworld.function.math module. Note that the Python math module is not compatible with Underworld2 functions (you must use our math module). Operator overloads are also provided from the indexing operator ([]) and the power operator (**).

We will construct some more double angle formula here. While in the previous example, we used Numpy to initialise a mesh variable object to construct sin and cos like functions, here all mathematical operations will be performed by Underworld Function objects.

In [7]:
# trig funcs
sin = fn.math.sin()
cos = fn.math.cos()
tan = fn.math.tan()

# double angle formula
sin_2theta = 2.*sin*cos
cos_2theta = 1. - 2.*sin**2
tan_2theta = (2.*tan)/(1.-tan**2)

# get somewhere to evaluate
theta = random.random()

# do things check out? 
test( sin_2theta.evaluate(theta), math.sin(2*theta) )
test( cos_2theta.evaluate(theta), math.cos(2*theta) )
test( tan_2theta.evaluate(theta), math.tan(2*theta) )

Relational and logical functions

Relational functions are constructed via the Python relational operators (<,<=,>,>=). Underworld functions for AND, OR and XOR logical operations are also available, and these overload the Python bitwise operators (&,|,^) (though they do not perform bitwise operations). These functions will all return boolean results.

In [8]:
# define a logical function for inside a unit circle.
# we will use the `input()` function as a proxy for the coordinate (see below).
radius = 1.
coord = fn.input()
inside_circle = ( coord[0]**2 + coord[1]**2 < radius**2 )

# test at some locations
test( inside_circle.evaluate( (0. ,0. ) ), True )
test( inside_circle.evaluate( (1. ,1. ) ), False )
test( inside_circle.evaluate( (0.5,0.5) ), True )
test( inside_circle.evaluate( (0.9,0.5) ), False )

# now something a bit more complex.. first, some circles
radius = 0.25
offset = (0.5,0.5)
deltax = radius*math.cos(math.pi/4.)
deltay = radius*math.sin(math.pi/4.)
coord1 = fn.input() - offset 
circle1 = ( fn.math.dot(coord1,coord1) < radius**2 )
coord2 = coord1 - ( deltax, deltay)
circle2 = ( fn.math.dot(coord2,coord2) < radius**2 )
coord3 = coord1 - (-deltax, deltay)
circle3 = ( fn.math.dot(coord3,coord3) < radius**2 )
coord4 = coord1 - ( deltax,-deltay)
circle4 = ( fn.math.dot(coord4,coord4) < radius**2 )
coord5 = coord1 - (-deltax,-deltay)
circle5 = ( fn.math.dot(coord5,coord5) < radius**2 )

# now create a cross.. note the use of the OR operator
cross =  ( (fn.math.abs(coord1[0])<radius/2.) | (fn.math.abs(coord1[1])<radius/2.) )

# visualise the XOR of these shapes
fig = glucifer.Figure()
fig.append(glucifer.objects.Surface(mesh, circle1^circle2^circle3^circle4^circle5^cross, resolution=400, colours=['white','blue']))
fig.show()

Conformal input & output checking

When you define your functions, they are input agnostic in the sense that you describe your function without declaring the type of argument that will be used when the function is eventually evaluated. At evaluation time, checks are performed to ensure the argument is compatible with the provided function. Likewise, the output returned by the evaluated function will be checked to ensure it is of the required form. If a check fails, an exception is thrown.

Some examples:

Input/Output types

In [9]:
# as mentioned already, arithmetic operations are only possible for `float` types
fn_test = vecvar + 1           # note that the definition works fine,
try:
    fn_test.evaluate((0.1,0.2))  # but the failure will occur at evaluation time.
except RuntimeError as e:
    print("RuntimeError: "+str(e))

RuntimeError: Operand in binary function does not appear to return a 'double' type value, as required. Note that where the operand Function you have constructed uses Python numeric objects, those objects must be of 'float' type (so for example '2.' instead of '2').

Cardinality

In [10]:
# the `Surface` object can only visualise scalar objects
fig = glucifer.Figure()
try:
    fig.append(glucifer.objects.Surface(mesh,vecvar))  # failure will occur here
except RuntimeError as e:
    print("RuntimeError: "+str(e))

RuntimeError: Provided function must return a scalar result.

Input/Output classes

In [11]:
# 'SwarmVariable' objects do not accept global coordinate evaluation
swarm = uw.swarm.Swarm(mesh)
svar = swarm.add_variable('double',1)
swarm.populate_using_layout(uw.swarm.layouts.GlobalSpaceFillerLayout(swarm,1))
try:
    svar.evaluate((0.2,0.3))
except RuntimeError as e:
    print("RuntimeError: "+str(e))

RuntimeError: Provided input to SwarmVariableFn does not appear to be of supported type. Supported types are 'FEMCoordinate' and 'ParticleCoordinate'. From Python, these are provided by respectively 'VoronoiIntegrationSwarm' and 'Swarm' type `FunctionInput` objects.

Function definitions are templates

When you are defining your functions, you are actually constructing a tree structure which templates the actual function definition. At evaluation time, this tree is walked starting with the provided input data, evaluating all subfunctions and passing through the results. So functions are dynamic in the sense that they do not capture a static representation of their constituent functions at definition time, but instead simply refer to these functions at evaluation time. This detail is mainly of importance when using data functions such as a MeshVariable objects, and means that any changes to the mesh variable will be reflected immediately in any functions which utilise it.

In [12]:
# create mesh var and init to 1.
meshvar = mesh.add_variable(1)
meshvar.data[:] = 1.

# now create a function which utilises the mesh var and evaluate
fn_test = meshvar * 5.
result_before = fn_test.evaluate((0.1,0.2)) 

# now modify the meshvar and evaluate the fn again
meshvar.data[:] = -1.
result_after = fn_test.evaluate((0.1,0.2)) 

print(result_before,result_after)

(array([[ 5.]]), array([[-5.]]))

Module Overview

The user is encouraged to drill down interactively into submodules to discover all available functionality. Documentation and examples are provided via class docstrings. A general overview of the top level is provided here:

Submodules:

analytic: Analytic solution functions.
branching: Functions which provide branching type behaviour.
exception: Functions which raise exceptions under certain conditions.
math: Basic mathematics functions.
misc: Miscellaneous other functions.
rheology: Functions which perform rheology specific behaviours.
shape: Functions which perform geometric operations.
tensor: Functions which perform operations on tensor functions.
view: Functions which only observe function evaluations without modifying results.

Classes:

coord: This function is an alias to the input() function.
input: This is the identity function. It simply returns its input.
Function: The Function base class.
FunctionInput: The FunctionInput base class.

The evaluate() Method

Once you have created functions, you will pass these into various objects within the Underworld2 API, however you will often also wish to directly evaluate the functions you have created for analytic or testing purposes. The evaluate() method provides this capability. Here are some basic examples using the function created previously:

In [13]:
# evaluate at single location, provide as a coordinate tuple or list.
fn_k.evaluate( (0.3,0.2) )
Out[13]:
array([[ 26.74625463]])
In [14]:
# evaluate at a set of locations, provide these as a numpy array.
import numpy as np
count = 10
# create an empty array
locations = np.zeros( (count,2) )
# specify evaluation coodinates
locations[:,0] = 0.5
locations[:,1] = np.linspace(0.,1.,count)
# evaluate
fn_k.evaluate(locations)
Out[14]:
array([[   13.        ],
       [   18.94327199],
       [   29.30185422],
       [   47.3559204 ],
       [   78.82251482],
       [  133.66592538],
       [  229.25299916],
       [  395.8525702 ],
       [  686.22046174],
       [ 1192.30527282]])

The FunctionInput class

A further type of input to function evaluation are FunctionInput class objects. These are shortcuts to their underlying data, but leverage the fundamental object data for higher efficiency evaluations, and sometimes as a necessity to facilitate the evaluation. Currently provided FunctionInput classes include:

  • underworld.mesh.FeMesh: Evaluation at all mesh vertex coordinates.
  • underworld.mesh.FeMesh_IndexSet: Evaluation at the mesh vertex coordinates within the set.
  • underworld.swarm.Swarm: Evaluation at all swarm particle coordinates.
  • underworld.swarm.VoronoiIntegrationSwarm: Evaluation at all voronoi swarm particle coordinates.

The above behave as FunctionInput classes by way of the multiple inheritence mechanisms of Python. For example, mesh objects (ie, objects of class FeMesh) are also FunctionInput objects.

Note that you will get identical evaluation results using the entire mesh vertex Numpy array as evaluation input, or using the mesh directly as an input:

In [15]:
# first lets confirm that mesh is indeed a `FunctionInput` object
if not isinstance(mesh, fn.FunctionInput):
    raise RuntimeError("Error! The mesh does not appear to be an instance of the `FunctionInput` class.")
    
# now evaluate at all mesh vertices:
results_using_functioninput = fn_k.evaluate(mesh)

# likewise, let's do the numpy equivalent
results_using_numpy = fn_k.evaluate(mesh.data)

# confirm identical results
if not np.allclose( results_using_functioninput, results_using_numpy ):
    raise RuntimeError("Error! Results differ where they should be the same.")

Note that the FeMesh_IndexSet which contains all mesh vertices should also return identical results:

In [16]:
# first create the 'empty' set
allindices = mesh.specialSets['Empty']

# now invert to obtain the 'full' set
allindices.invert()

# evaluate
results_using_indexset = fn_k.evaluate(allindices)

# again confirm identical results
if not np.allclose( results_using_functioninput, results_using_indexset ):
    raise RuntimeError("Error! Results differ where they should be the same.")

While all the above methods yield the same numerical results, note that there are efficiency differences. When you provide a Numpy array as input to a MeshVariable object evaluation, it is first necessary to determine which element the evaluation coordinate belongs to, which can be an expensive operation (in particular for deformed mesh). Once the owning element is determined, the interpolation itself needs to be calculated. However, when using the mesh object directly as a function input, the evaluation leverages higher level object information to directly extract the nodal values (ie, owning element and interpolation are not required). Consider the following timing results:

In [17]:
# first evaluate directly using mesh
%timeit tempVar.evaluate(mesh)

1000 loops, best of 3: 426 µs per loop
In [18]:
# now evaluate via numpy array
%timeit tempVar.evaluate(mesh.data)

100 loops, best of 3: 4.2 ms per loop
In [19]:
# now flag the mesh as deformed. note that we are not
# moving any mesh vertices, but simply flagging the 
# mesh as deformed.
with mesh.deform_mesh():
    pass
%timeit fn_k.evaluate(mesh.data)
mesh.reset()       # reset the mesh for future operations

10 loops, best of 3: 96.9 ms per loop

So we see that there isn't a dramatic difference until the mesh is deformed. Note also that as we haven't truly deformed the mesh, the observed result will generally be the best case scenario.

The other important FunctionInput objects are the swarm based items. In the case of functions which utilise SwarmVariable objects, these are the only option for evaluation from within Python. This is because particles are inherently discrete, and therefore evaluation at arbitrary locations is not directly possible. In some case where you supply SwarmVariable functions to be used within Underworld operations (such as a viscosity for the Stokes system), a nearest neighbour calculation is implicitly utilised, however these mechanisms are not (yet) exposed for calls to the evaluate() method. Future releases will provide greater functionality for operating over SwarmVariable functions. Again, simply pass the swarm to the evaluate function to utilse the swarm FunctionInput behaviour:

In [20]:
# first add a swarm to the existing mesh and populate
swarm = uw.swarm.Swarm(mesh)
swarm.populate_using_layout( uw.swarm.layouts.GlobalSpaceFillerLayout(swarm,20) )

# add variable to use for conductivity
k_var = swarm.add_variable('double',1)

# now initialise these particles to have values which correspond to those
# provided by fn_k.  we can use either `swarm` or `swarm.data` as inputs 
# here and the results in this case will be identical.
k_var.data[:] = fn_k.evaluate(swarm)

# note that k_var is itself of the `Function` class
if not isinstance(k_var, fn.Function):
    raise RuntimeError("Error! Swarm variable does not appear to be of `Function` class")
    
# now try to evaluate at arbitrary location
# k_var.evaluate( (0.3,0.2) )   # this will raise a RuntimeError!

# now evaluate using swarm 
results_using_swarm_k_var = k_var.evaluate( swarm )
# also evaluate fn_k using swarm
results_using_swarm_fn_k = fn_k.evaluate( swarm )
if not np.allclose( results_using_swarm_fn_k, results_using_swarm_fn_k ):
    raise RuntimeError("Error! These arrays should have identical values.")

The input function

The function.input class simply provides the identity function. It returns whatever values are passed to it!

In [21]:
# create the input class/function
fn_input = fn.input()

# confirm behaviour
print(fn_input.evaluate( 2. ) )
print(fn_input.evaluate( ((0.2), (0.3)) ))
print(fn_input.evaluate( ((0.2), (0.3), (0.4)) ))
print(fn_input.evaluate( mesh.data[0:5] ))

if not np.allclose(mesh.data, fn_input.evaluate(mesh.data)):
    raise RuntimeError("Error! These arrays should have identical values.")

[[ 2.]]
[[ 0.2  0.3]]
[[ 0.2  0.3  0.4]]
[[ 0.        0.      ]
 [ 0.015625  0.      ]
 [ 0.03125   0.      ]
 [ 0.046875  0.      ]
 [ 0.0625    0.      ]]

The input function is most often used when you wish to construct functions which operate over a particular coordinate axis:

In [22]:
# Create sinusoidal function in the horizontal direction.
# Note that the square brackets operator is overloaded to extract 
# a certain axis' value from vector functions.
fn_sin = fn.math.sin(8.*np.pi*fn_input[0])
# take a look
fig = glucifer.Figure()
fig.append(glucifer.objects.Surface(mesh,fn_sin))
fig.show()

We also provide an alias to this class in function.coord. This is strictly a Python alias, and is identical in all ways (except name!) to function.input. It is provided because we often use the input function to operate on coordinate values (as above), but it is worth noting that function inputs are certainly not restricted to coordinates or vectors, and in these cases it would be misleading to use the coord function.

Branching Functions

Branching functions are functions which nest other functions, and select which function to execute based on some condition (also expressed as a function!). The most common use case for branching functions are for defining materials within your model which exhibit different behaviours (based on their associated function). Note that Underworld2 does not include an explict 'material' construct (unlike Underworld1), however the equivalent functionality may be obtained through branching functions.

Two branching functions are currently provided:

function.branching.map: Key/value type mapping of behavior.
function.branching.conditional: if/elif type behaviour.

Refer to respective API documentation for further information on these function classes. The map function is usually used to construct material type behaviours. Note that the map function provides a subset of the conditional function behaviour, though it has performance advantages.

Let us construct a basic model with material type behaviours:

In [23]:
# setup a mesh and swarm
mesh = uw.mesh.FeMesh_Cartesian(elementRes=(8,8),minCoord=(-1.0, -1.0), maxCoord=(1.0, 1.0))
swarm = uw.swarm.Swarm(mesh)
swarm.populate_using_layout(uw.swarm.layouts.GlobalSpaceFillerLayout(swarm,80))

# add a variable which will act as the key function for the map function.
material_index = swarm.add_variable("int",1)

# add the following for convenience and clarity.
outside_circle = 0
inside_circle = 1

# initialise material_index to be outside_circle everywhere
material_index.data[:] = outside_circle
# now set to inside_circle where inside circle!
r2 = 0.8
for index, position in enumerate(swarm.particleCoordinates.data):
    if position[0]**2 + position[1]**2 < r2:
        material_index.data[index] = inside_circle

# create mapped behaviour dictionary
fn_z = fn.input()[1]
map_dict = { outside_circle:fn_z,
              inside_circle:fn_sin }
# create function
fn_map = fn.branching.map( fn_key = material_index, 
                          mapping = map_dict )

# viz
fig = glucifer.Figure()
fig.append( glucifer.objects.Points(swarm,fn_map, pointSize=10.) )
fig.show()

Alternatively, we can achieve an identical result using the conditional function, and its evaluate() method:

In [24]:
# the 'position' variable (created above) retains a reference to the 
# particleCoordinates numpy array. this inteferes with adding a new
# swarm variable, so we must delete it. 
del position  

# add a new variable first (we will test later for equality)
material_index2 = swarm.add_variable("int",1)

# first define a function which returns true if inside the circle
coord = fn.input()
fn_in_circle = fn.math.dot(coord,coord) < r2  # note the overload of the '<' operator

# now create the conditional 
fn_conditional = fn.branching.conditional(  ( (fn_in_circle,  inside_circle),
                                              (        True, outside_circle)   ) )

# use the evaluate, writing the results out to the material_index2 
material_index2.data[:] = fn_conditional.evaluate(swarm)
# check to ensure that we have identical results
if not (material_index2.data == material_index.data).all():
    raise RuntimeError("Error! These arrays should have identical values.")