How to use the underworld.systems.sle function in underworld

def __init__(self, rowVector, colVector, rhs=None, rhs_T=None, assembleOnNodes=False, **kwargs):
        if not isinstance(rowVector, uw.systems.sle.SolutionVector):
            raise TypeError("'rowVector' object passed in must be of type 'SolutionVector'")

        if not isinstance(colVector, uw.systems.sle.SolutionVector):
            raise TypeError("'colVector' object passed in must be of type 'SolutionVector'")

        self._meshVariableRow = rowVector.meshVariable
        self._meshVariableCol = colVector.meshVariable

        if rhs and not isinstance(rhs, _assembledvector.AssembledVector):
            raise TypeError("'rhs' object passed in must be of type 'AssembledVector'")

        if rhs_T and not isinstance(rhs_T, _assembledvector.AssembledVector):
            raise TypeError("'rhs_T' object passed in must be of type 'AssembledVector'")

        self._rhs   = rhs
        self._rhs_T = rhs_T

vnsField     = self._vnsField = self._velocityField.copy()
        vnsEqNum     = uw.systems.sle.EqNumber( vnsField, False )
        self._vnsVec = uw.systems.sle.SolutionVector(vnsField, vnsEqNum) # store on class

        # evaluate vnsField, check compatibility
        if isinstance(mesh, uw.mesh.FeMesh_Annulus): # with FeMesh_Annulus we have the sbr_fn
            self._vnsVec.meshVariable.data[:] = mesh.sbr_fn.evaluate(mesh)

        uw.libUnderworld.StgFEM.SolutionVector_LoadCurrentFeVariableValuesOntoVector(self._vnsVec._cself) # store in petsc vec

        # must be done after vnsField creation
        self._velocityField._cself.nonAABCs = 1

        # build a global 're-rotate' matrix
        self._rot = uw.systems.sle.AssembledMatrix( self._vnsVec, self._vnsVec, rhs=None )
        gaussSwarm = self._constitMatTerm._integrationSwarm
        self._rot._cself.assembleOnNodes = 1 # important doesn't perform FEM integral

        term = self._term = uw.systems.sle.MatrixAssemblyTerm_RotationDof( integrationSwarm = gaussSwarm,
                                                                           assembledObject  = self._rot,
                                                                           fn_basis         = basis_vectors,
                                                                           mesh             = mesh)

        vnsEqNum._cself.removeBCs = True
        uw.libUnderworld.StgFEM.StiffnessMatrix_Assemble(self._rot._cself, None, None);
        vnsEqNum._cself.removeBCs = False

        # add rotation matrix element terms using the following
        uw.libUnderworld.StgFEM.StiffnessMatrix_SetRotationTerm(self._kmatrix._cself, term._cself)
        uw.libUnderworld.StgFEM.StiffnessMatrix_SetRotationTerm(self._gmatrix._cself, term._cself)
        uw.libUnderworld.StgFEM.ForceVector_SetRotationTerm(self._fvector._cself, term._cself)

def addRotationMatrix(self, rMat):
        if rMat:
            if not isinstance( rMat, uw.systems.sle.AssembledMatrix):
                raise TypeError( "Provided 'rMat' must be of or convertible to 'AssembledMatrix' class." )
            # TODO: write a setter for the following 2 lines
            self.rMat = rMat
            self._cself.rMat = rMat._cself
        return

assembledObject=self._kmatrix,
                                                                 fn=_fn_diffusivity)
        self._forceVecTerm = sle.VectorAssemblyTerm_NA__Fn(   integrationSwarm=intswarm,
                                                              assembledObject=self._fvector,
                                                              fn=fn_heating)
        # search for neumann conditions
        for cond in conditions:
            if isinstance( cond, uw.conditions.NeumannCondition ):
                #NOTE many NeumannConditions can be used but the _sufaceFluxTerm only records the last

                ### -VE flux because of the FEM discretisation method of the initial equation
                negativeCond = uw.conditions.NeumannCondition( fn_flux=cond.fn_flux,
                                                               variable=cond.variable,
                                                               indexSetsPerDof=cond.indexSetsPerDof )

                self._surfaceFluxTerm = sle.VectorSurfaceAssemblyTerm_NA__Fn__ni(
                                                                assembledObject  = self._fvector,
                                                                surfaceGaussPoints = 2,
                                                                nbc         = negativeCond )
        super(SteadyStateHeat, self).__init__(**kwargs)

def _create_penalty_objects(self):
        # using this function se we don't need to add anything extra to the stokeSLE struct

        velocityField=self._velocityField
        pressureField=self._pressureField
        stokesSLE = self._stokesSLE
        # create junk force vectors -- we provide no assembly terms for these so they are 0 vectors.
        self._vmfvector   = sle.AssembledVector(velocityField, stokesSLE._eqNums[velocityField])
        self._junkfvector = sle.AssembledVector(pressureField, stokesSLE._eqNums[pressureField])

        # and matrices
        self._vmmatrix = sle.AssembledMatrix( stokesSLE._velocitySol, stokesSLE._velocitySol, rhs=self._vmfvector )
        self._mmatrix  = sle.AssembledMatrix( stokesSLE._pressureSol, stokesSLE._pressureSol, rhs=self._junkfvector )

        # create assembly terms
        self._pressMassMatTerm = sle.MatrixAssemblyTerm_NA__NB__Fn( integrationSwarm=uw.swarm.GaussIntegrationSwarm(velocityField.mesh),
                                                             fn=1.0, assembledObject=self._mmatrix,
                                                             mesh = velocityField._mesh)

        # attach terms to live solver struct
        self._cself.vmForceVec = self._vmfvector._cself
        self._cself.vmStiffMat = self._vmmatrix._cself
        self._cself.jForceVec  = self._junkfvector._cself
        self._cself.mStiffMat  = self._mmatrix._cself

def _create_penalty_objects(self):
        # using this function se we don't need to add anything extra to the stokeSLE struct

        velocityField=self._velocityField
        pressureField=self._pressureField
        stokesSLE = self._stokesSLE
        # create junk force vectors -- we provide no assembly terms for these so they are 0 vectors.
        self._vmfvector   = sle.AssembledVector(velocityField, stokesSLE._eqNums[velocityField])
        self._junkfvector = sle.AssembledVector(pressureField, stokesSLE._eqNums[pressureField])

        # and matrices
        self._vmmatrix = sle.AssembledMatrix( stokesSLE._velocitySol, stokesSLE._velocitySol, rhs=self._vmfvector )
        self._mmatrix  = sle.AssembledMatrix( stokesSLE._pressureSol, stokesSLE._pressureSol, rhs=self._junkfvector )

        # create assembly terms
        self._pressMassMatTerm = sle.MatrixAssemblyTerm_NA__NB__Fn( integrationSwarm=uw.swarm.GaussIntegrationSwarm(velocityField.mesh),
                                                             fn=1.0, assembledObject=self._mmatrix,
                                                             mesh = velocityField._mesh)

        # attach terms to live solver struct
        self._cself.vmForceVec = self._vmfvector._cself
        self._cself.vmStiffMat = self._vmmatrix._cself
        self._cself.jForceVec  = self._junkfvector._cself
        self._cself.mStiffMat  = self._mmatrix._cself

fn   = -0.5 * fn_dt * self.fn_diffusivity * self._phiStar.fn_gradient )

        if nbc is not None:
            # -VE flux because of the FEM discretisation method of the initial equation
            negativeCond = uw.conditions.NeumannCondition( 
                            fn_flux  = fn_dt * nbc.fn_flux,
                            variable = nbc.variable,
                            indexSetsPerDof = nbc.indexSetsPerDof )

            # NOTE many NeumannConditions can be used but the _sufaceFluxTerm only records the last
            self._surfaceFluxTerm = sle.VectorSurfaceAssemblyTerm_NA__Fn__ni(
                                        assembledObject    = f,
                                        surfaceGaussPoints = 2,
                                        nbc         = negativeCond )

        self._k_term = uw.systems.sle.MatrixAssemblyTerm_NA_i__NB_i__Fn(
                        assembledObject  = K,
                        integrationSwarm = intSwarm,
                        fn = 0.5 * fn_dt * self.fn_diffusivity)

        self._m_term = uw.systems.sle.MatrixAssemblyTerm_NA__NB__Fn(
                        assembledObject  = K,
                        integrationSwarm = intSwarm,
                        fn   = 1.,
                        mesh = mesh)

        # functions used to calculate the timestep, see function get_max_dt()

        self._maxVsq  = uw.function.view.min_max(velocityField, fn_norm = uw.function.math.dot(velocityField, velocityField) )
        self._maxDiff = uw.function.view.min_max(self.fn_diffusivity)

        # Note that the c level minSep on the mesh is for the local domain

# set the bcs on here
            if type(cond) == uw.conditions.NeumannCondition:
                if nbc != None:
                    # check only one nbc condition is given in 'conditions' list
                    RuntimeError( "Provided 'conditions' can only accept one NeumannConditions condition object.")
            elif type(cond) == uw.conditions.DirichletCondition:
                if cond.variable == self._phiField:
                    libUnderworld.StgFEM.FeVariable_SetBC( self._phiField._cself, cond._cself )
                    libUnderworld.StgFEM.FeVariable_SetBC( self._phiDotField._cself, cond._cself )
            else:
                raise RuntimeError("Input condition type not recognised.")
        self._conditions = conditions

        # force removal of BCs as SUPG cannot handle leaving them in
        self._eqNumPhi    = sle.EqNumber( phiField, removeBCs=True )
        self._eqNumPhiDot = sle.EqNumber( phiDotField, removeBCs=True )

        self._phiSolution    = sle.SolutionVector( phiField, self._eqNumPhi )
        self._phiDotSolution = sle.SolutionVector( phiDotField, self._eqNumPhiDot )

        # create force vectors
        self._residualVector = sle.AssembledVector(phiField, self._eqNumPhi )
        self._massVector     = sle.AssembledVector(phiField, self._eqNumPhi )

        # create swarm
        self._gaussSwarm = uw.swarm.GaussIntegrationSwarm(self._phiField.mesh)

        super(_SUPG_AdvectionDiffusion, self).__init__()

        self._cself.phiVector = self._phiSolution._cself
        self._cself.phiDotVector = self._phiDotSolution._cself

for cond in self._conditions:
            if isinstance( cond, uw.conditions.NeumannCondition ):
                #NOTE many NeumannConditions can be used but the _sufaceFluxTerm only records the last
                self._surfaceFluxTerm = sle.VectorSurfaceAssemblyTerm_NA__Fn__ni(
                                                                assembledObject    = self._fvector,
                                                                surfaceGaussPoints = 3, # increase to resolve stress bc fluctuations
                                                                nbc                = cond )
        if self._fn_minus_one_on_lambda != None:
            # add matrix and associated assembly term for compressible stokes formulation
            # a mass matrix goes into the lower right block of the stokes system coeff matrix
            self._mmatrix = sle.AssembledMatrix( self._pressureSol, self._pressureSol, rhs=self._hvector )
            # -1. as per Hughes, The Finite Element Method, 1987, Table 4.3.1, [M]

            self._compressibleTerm = sle.MatrixAssemblyTerm_NA__NB__Fn( integrationSwarm = intswarm,
                                                                        assembledObject  = self._mmatrix,
                                                                        mesh             =  mesh,
                                                                        fn               = self._fn_minus_one_on_lambda )

        if fn_stresshistory != None:
            self._NA_j__Fn_ijTerm    = sle.VectorAssemblyTerm_NA_j__Fn_ij( integrationSwarm = intswarm, 
                                                                           assembledObject  = self._fvector, 
                                                                           fn               = fn_stresshistory )
        # objects used for analysis, dictionary for organisation
        self._aObjects = dict()
        self._aObjects['vdotv_fn'] = uw.function.math.dot(self._velocityField, self._velocityField)

        super(Curvilinear_Stokes, self).__init__(**kwargs)

        for cond in self._conditions:
            if type(cond) in [uw.conditions.RotatedDirichletCondition, uw.conditions.CurvilinearDirichletCondition]:

def _create_penalty_objects(self):
        # using this function se we don't need to add anything extra to the stokeSLE struct

        velocityField=self._velocityField
        pressureField=self._pressureField
        stokesSLE = self._stokesSLE
        # create junk force vectors -- we provide no assembly terms for these so they are 0 vectors.
        self._vmfvector   = sle.AssembledVector(velocityField, stokesSLE._eqNums[velocityField])
        self._junkfvector = sle.AssembledVector(pressureField, stokesSLE._eqNums[pressureField])

        # and matrices
        self._vmmatrix = sle.AssembledMatrix( stokesSLE._velocitySol, stokesSLE._velocitySol, rhs=self._vmfvector )
        self._mmatrix  = sle.AssembledMatrix( stokesSLE._pressureSol, stokesSLE._pressureSol, rhs=self._junkfvector )

        # create assembly terms
        self._pressMassMatTerm = sle.MatrixAssemblyTerm_NA__NB__Fn( integrationSwarm=uw.swarm.GaussIntegrationSwarm(velocityField.mesh),
                                                             fn=1.0, assembledObject=self._mmatrix,
                                                             mesh = velocityField._mesh)

        # attach terms to live solver struct
        self._cself.vmForceVec = self._vmfvector._cself
        self._cself.vmStiffMat = self._vmmatrix._cself
        self._cself.jForceVec  = self._junkfvector._cself
        self._cself.mStiffMat  = self._mmatrix._cself