dolfin team mailing list archive

[Chaffra <question144683@xxxxxxxxxxxxxxxxxxxxx>]

New question #144683 on DOLFIN:
https://answers.launchpad.net/dolfin/+question/144683

Dear all:

Is it possible to create a custom ufl mathfunction to use in a nonlinear problem? In the script below I can create a custom math function "fermi_integral_ufl" but I run into problem when I try to solve let's div(grad(u)) = fermi_integral_ufl(u). I get:

/usr/lib/python2.6/dist-packages/ufl/algorithms/forward_ad.pyc in math_function(self, o, a)
    371 
    372     def math_function(self, o, a):
--> 373         error("Unknown math function.")
    374 
    375     def sqrt(self, o, a):

/usr/lib/python2.6/dist-packages/ufl/log.pyc in error(self, *message)
    122         "Write error message and raise an exception."
    123         self._log.error(*message)
--> 124         raise UFLException(self._format_raw(*message))
    125 
    126     def begin(self, *message):

UFLException: Unknown math function.

when solving the problem. I think, I am running into problems when computing the deriravative of fermi_integral_ufl. Any thoughts?

Thanks,
Chaffra

my script ------------------

'''
Created on Feb 8, 2011

@author: Chaffra Affouda
'''

from scipy.special import gamma
from scipy import exp, log as ln, Inf
from scipy.integrate import quad
from numpy import vectorize
from ufl.mathfunctions import MathFunction
from ufl.constantvalue import ScalarValue, Zero, FloatValue, as_ufl
 

def fermi_function(t,j, x): return pow(t,j)/(exp(t-x) + 1)

def fermi_integral(order,x, 
                   full_output=0, 
                   epsabs=1e-13, 
                   epsrel=1e-13, 
                   limit=50,):
    if order == 0: 
        return ln(1+exp(x))
    else:
        val = quad(fermi_function, 0, Inf, args=(order,x),
                   full_output=full_output,
                   epsabs=epsabs,
                   epsrel=epsrel,
                   limit=limit,)
        return gamma(order+1)*val[0]
fermi_integral_vector = vectorize(fermi_integral)

class FermiIntegralUFL(MathFunction):
    
    def __new__(cls, order, argument):
        if isinstance(argument, (ScalarValue, Zero)):
            return FloatValue(fermi_integral(order,float(argument)))
        return MathFunction.__new__(cls)
    
    def evaluate(self, x, mapping, component, index_values):    
        a = self._argument.evaluate(x, mapping, component, index_values)
        return fermi_integral(self.order, a)
    
    def __init__(self, order, argument):
        MathFunction.__init__(self, "fermi_integral_"+str(order), argument)
        self.order = order
        
def fermi_integral_ufl(order,f):
    "The Fermi integral of f."
    f = as_ufl(f)
    r = FermiIntegralUFL(order,f)
    if isinstance(r, (ScalarValue, Zero, int, float)):
        return float(r)
    return r


-- 
You received this question notification because you are a member of
DOLFIN Team, which is an answer contact for DOLFIN.