How to use the underworld.function function in underworld

variable=tField,
                                                    indexSetsPerDof=(inner+outer) )

        self.dirichletBCs['velocity_freeSlip'] =  uw.conditions.RotatedDirichletCondition(
                                    variable        = vField,
                                    indexSetsPerDof = (inner+outer, None),
                                    basis_vectors   = (annulus.bnd_vec_normal, annulus.bnd_vec_tangent))

        self.dirichletBCs['velocity_noSlip']   =  uw.conditions.RotatedDirichletCondition(
                                    variable        = vField,
                                    indexSetsPerDof = (inner+outer, None),
                                    basis_vectors   = (annulus.bnd_vec_normal, annulus.bnd_vec_tangent))

        # set up analytics logging
        self.f     = radialLengths[0]/radialLengths[1]
        self.dT_dr = fn.math.dot( tField.fn_gradient, annulus._fn_unitvec_radial() )

        self.dT_dr_outer_integral  = uw.utils.Integral( mesh=annulus, fn=self.dT_dr,
                                                   integrationType="surface", surfaceIndexSet=outer )

        self.dT_dr_inner_integral  = uw.utils.Integral( mesh=annulus, fn=self.dT_dr,
                                                   integrationType="surface", surfaceIndexSet=inner )

        # start the log file
        self.log_titles += ['','Nu_u','Nu_b']

        # setup visualisation with lavavu
        animation = self.animation
        self.view = glucifer.Figure(store=animation, name="scene1")
        # self.view.append(glucifer.objects.Mesh(annulus))
        self.view.append(glucifer.objects.Surface(mesh=annulus, fn=self.fields['temperature'],
                                              onMesh=True, name='temperature'))

def layer(top, bottom, minX, maxX):
    if bottom == None or top == None:
        return None
    vertex_array = np.array( [(minX, bottom),(minX, top),(maxX, top),(maxX, bottom)] )
    return uw.function.shape.Polygon(vertex_array)

import underworld as uw
from underworld import function as fn
import numpy as np

nSpheres = 1
fn_conds = []
for s_i in xrange(nSpheres):
    s,lati,longi, radius =  (7, 0, 0, 1.3)
    #s,lati,longi, radius =  (3, 180, 360, 2.5)*np.random.rand(4) + (6., -90, 0, 0.5)
    print s, lati, longi, radius
    lati,longi = np.radians([lati, longi])
    center = (s * np.cos(lati) * np.cos(longi),
              s * np.cos(lati) * np.sin(longi),
              s * np.sin(lati)                 )
    fn_pos = fn.coord() - center
    fn_conds.append( (fn.math.dot( fn_pos, fn_pos ) < radius**2, 1.) )
    
fn_conds.append( (True, 0.0 ) ) # default condition
fn_density = fn.branching.conditional( fn_conds )

fn_r = fn.math.sqrt( fn.math.dot( fn.coord(), fn.coord() ) )
fn_gravity =  fn.misc.constant(-9.8) / fn_r * (fn.coord()[0], fn.coord()[1], fn.coord()[2])

'''
sphereRadius = 0.2         # define radius
spheres = [ (.5,.5,0.7),   # define list of sphere centres
            (.25,.25,0.6),
            (.25,.75, 0.5),
            (.75,.75,0.6),
            (.75,.25,0.4) ]

>>> fn_1 = uw.function.misc.constant(2.0)
        >>> np.allclose( mesh.integrate( fn_1 )[0], 4 )
        True

        >>> fn_2 = uw.function.misc.constant(2.0) * (0.5, 1.0)
        >>> np.allclose( mesh.integrate( fn_2 ), [2,4] )
        True

        """
        # julian, i've dumbed this down for now as i'm not sure how to handle the
        # circular dependency.  i think the performance hit will be generally
        # negligible
        _volume_integral = uw.utils.Integral(mesh=self, fn=fn)
        return _volume_integral.evaluate()

class FeMesh_IndexSet(uw.container.ObjectifiedIndexSet, function.FunctionInput):
    """
    This class ties the FeMesh instance to an index set, and stores other
    metadata relevant to the set.

    Parameters
    ----------
    object: underworld.mesh.FeMesh
        The FeMesh instance from which the IndexSet was extracted.
    topologicalIndex: int
        Mesh topological index for which the IndexSet relates. See
        docstring for further info.

    Example
    -------

    >>> amesh = uw.mesh.FeMesh_Cartesian( elementType='Q1/dQ0', elementRes=(4,4), minCoord=(0.,0.), maxCoord=(1.,1.) )

def inCircleFnGenerator(centre, radius):
    coord = fn.input()
    offsetFn = coord - centre
    return fn.math.dot( offsetFn, offsetFn ) < radius**2

# **Create variables required for plasticity calculations**

# In[ ]:

secinv = fn.tensor.second_invariant( fn.tensor.symmetric( velocityField.gradientFn ) )
coordinate = fn.coord()


# **Setup viscosity functions**
# 
# Remember to use floats everywhere when setting up functions

# In[ ]:

viscosityl1 = fn.math.exp(math.log(ETA_T)*-1.*temperatureField)
viscosityl2 = fn.math.exp((math.log(ETA_T)*-1.*temperatureField) + (1.-coordinate[1])*math.log(ETA_Y))

#Von Mises effective viscosity
viscosityp = ETA0 + YSTRESS/(secinv/math.sqrt(0.5)) #extra factor to account for underworld second invariant form

if CASE == 1:
    viscosityFn = viscosityl1
elif CASE == 2:
    viscosityFn = 2./(1./viscosityl1 + 1./viscosityp)
elif CASE == 3:
    viscosityFn = viscosityl2
else:
    viscosityFn = 2./(1./viscosityl2 + 1./viscosityp)


# **Add functions for density and buoyancy**

# ### Material distribution in the domain.
# 
# 

# In[11]:

# Initialise the 'materialVariable' data to represent different materials. 
materialV = 0 # viscoplastic
materialW = 1 # weak
materialA = 2 # accommodation layer a.k.a. Sticky Air

# The particle coordinates will be the input to the function evaluate (see final line in this cell).
# We get proxy for this now using the input() function.
coord = fn.input()

# Setup the conditions list for the following conditional function. Where the
# z coordinate (coordinate[1]) is less than the perturbation, set to lightIndex.
conditions = [ (                                  coord[1] > thicknessV , materialA ),
               ( ((coord[1] < dWeak) & (coord[0]**2. < (dWeak**2.)/4.)) , materialW ),
               (                                                   True , materialV ) ]

# The actual function evaluation. Here the conditional function is evaluated at the location
# of each swarm particle. The results are then written to the materialVariable swarm variable.
materialVariable.data[:] = fn.branching.conditional( conditions ).evaluate(swarm)


# Define the density function
# ---

# In[12]:

def _effectiveViscosity(self):
        return fn.misc.constant(nd(self._viscosity))