dolfin team mailing list archive

[Johan Hake <johan.hake@xxxxxxxxx>]

Chaffra!

If this still is a problem maybe you can re-post the question at:

  https://answers.launchpad.net/ufl

Johan


On Tuesday February 8 2011 19:15:47 Chaffra wrote:
> 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