Question #144683 on DOLFIN changed: https://answers.launchpad.net/dolfin/+question/144683
Johan Hake posted a new comment: 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 -- You received this question notification because you are a member of DOLFIN Team, which is an answer contact for DOLFIN. _______________________________________________ Mailing list: https://launchpad.net/~dolfin Post to : [email protected] Unsubscribe : https://launchpad.net/~dolfin More help : https://help.launchpad.net/ListHelp
