Source code for ufl.mathfunctions

# -*- coding: utf-8 -*-
"""This module provides basic mathematical functions."""

# Copyright (C) 2008-2016 Martin Sandve Aln├Žs
#
# This file is part of UFL (https://www.fenicsproject.org)
#
# SPDX-License-Identifier:    LGPL-3.0-or-later
#
# Modified by Anders Logg, 2008
# Modified by Kristian B. Oelgaard, 2011

import math
import cmath
import numbers

from ufl.log import warning, error
from ufl.core.operator import Operator
from ufl.core.ufl_type import ufl_type
from ufl.constantvalue import is_true_ufl_scalar, Zero, RealValue, FloatValue, IntValue, ComplexValue, ConstantValue, as_ufl

"""
TODO: Include additional functions available in <cmath> (need derivatives as well):

Exponential and logarithmic functions:
log10    Compute common logarithm (function)

TODO: Any other useful special functions?

About bessel functions:
http://en.wikipedia.org/wiki/Bessel_function

Portable implementations of bessel functions:
http://www.boost.org/doc/libs/1_47_0/libs/math/doc/sf_and_dist/html/math_toolkit/main_overview/tr1.html

Implementation in C++ std::tr1:: or boost::math::tr1::
- BesselK: cyl_bessel_k(nu, x)
- BesselI: cyl_bessel_i(nu, x)
- BesselJ: cyl_bessel_j(nu, x)
- BesselY: cyl_neumann(nu, x)
"""


# --- Function representations ---

[docs]@ufl_type(is_abstract=True, is_scalar=True, num_ops=1) class MathFunction(Operator): "Base class for all unary scalar math functions." # Freeze member variables for objects in this class __slots__ = ("_name",) def __init__(self, name, argument): Operator.__init__(self, (argument,)) if not is_true_ufl_scalar(argument): error("Expecting scalar argument.") self._name = name
[docs] def evaluate(self, x, mapping, component, index_values): a = self.ufl_operands[0].evaluate(x, mapping, component, index_values) try: if isinstance(a, numbers.Real): res = getattr(math, self._name)(a) else: res = getattr(cmath, self._name)(a) except ValueError: warning('Value error in evaluation of function %s with argument %s.' % (self._name, a)) raise return res
def __str__(self): return "%s(%s)" % (self._name, self.ufl_operands[0])
[docs]@ufl_type() class Sqrt(MathFunction): __slots__ = () def __new__(cls, argument): if isinstance(argument, (RealValue, Zero, numbers.Real)): if float(argument) < 0: return ComplexValue(cmath.sqrt(complex(argument))) else: return FloatValue(math.sqrt(float(argument))) if isinstance(argument, (ComplexValue, complex)): return ComplexValue(cmath.sqrt(complex(argument))) return MathFunction.__new__(cls) def __init__(self, argument): MathFunction.__init__(self, "sqrt", argument)
[docs]@ufl_type() class Exp(MathFunction): __slots__ = () def __new__(cls, argument): if isinstance(argument, (RealValue, Zero)): return FloatValue(math.exp(float(argument))) if isinstance(argument, (ComplexValue)): return ComplexValue(cmath.exp(complex(argument))) return MathFunction.__new__(cls) def __init__(self, argument): MathFunction.__init__(self, "exp", argument)
[docs]@ufl_type() class Ln(MathFunction): __slots__ = () def __new__(cls, argument): if isinstance(argument, (RealValue, Zero)): return FloatValue(math.log(float(argument))) if isinstance(argument, (ComplexValue)): return ComplexValue(cmath.log(complex(argument))) return MathFunction.__new__(cls) def __init__(self, argument): MathFunction.__init__(self, "ln", argument)
[docs] def evaluate(self, x, mapping, component, index_values): a = self.ufl_operands[0].evaluate(x, mapping, component, index_values) try: return math.log(a) except TypeError: return cmath.log(a)
[docs]@ufl_type() class Cos(MathFunction): __slots__ = () def __new__(cls, argument): if isinstance(argument, (RealValue, Zero)): return FloatValue(math.cos(float(argument))) if isinstance(argument, (ComplexValue)): return ComplexValue(cmath.cos(complex(argument))) return MathFunction.__new__(cls) def __init__(self, argument): MathFunction.__init__(self, "cos", argument)
[docs]@ufl_type() class Sin(MathFunction): __slots__ = () def __new__(cls, argument): if isinstance(argument, (RealValue, Zero)): return FloatValue(math.sin(float(argument))) if isinstance(argument, (ComplexValue)): return ComplexValue(cmath.sin(complex(argument))) return MathFunction.__new__(cls) def __init__(self, argument): MathFunction.__init__(self, "sin", argument)
[docs]@ufl_type() class Tan(MathFunction): __slots__ = () def __new__(cls, argument): if isinstance(argument, (RealValue, Zero)): return FloatValue(math.tan(float(argument))) if isinstance(argument, (ComplexValue)): return ComplexValue(cmath.tan(complex(argument))) return MathFunction.__new__(cls) def __init__(self, argument): MathFunction.__init__(self, "tan", argument)
[docs]@ufl_type() class Cosh(MathFunction): __slots__ = () def __new__(cls, argument): if isinstance(argument, (RealValue, Zero)): return FloatValue(math.cosh(float(argument))) if isinstance(argument, (ComplexValue)): return ComplexValue(cmath.cosh(complex(argument))) return MathFunction.__new__(cls) def __init__(self, argument): MathFunction.__init__(self, "cosh", argument)
[docs]@ufl_type() class Sinh(MathFunction): __slots__ = () def __new__(cls, argument): if isinstance(argument, (RealValue, Zero)): return FloatValue(math.sinh(float(argument))) if isinstance(argument, (ComplexValue)): return ComplexValue(cmath.sinh(complex(argument))) return MathFunction.__new__(cls) def __init__(self, argument): MathFunction.__init__(self, "sinh", argument)
[docs]@ufl_type() class Tanh(MathFunction): __slots__ = () def __new__(cls, argument): if isinstance(argument, (RealValue, Zero)): return FloatValue(math.tanh(float(argument))) if isinstance(argument, (ComplexValue)): return ComplexValue(cmath.tanh(complex(argument))) return MathFunction.__new__(cls) def __init__(self, argument): MathFunction.__init__(self, "tanh", argument)
[docs]@ufl_type() class Acos(MathFunction): __slots__ = () def __new__(cls, argument): if isinstance(argument, (RealValue, Zero)): return FloatValue(math.acos(float(argument))) if isinstance(argument, (ComplexValue)): return ComplexValue(cmath.acos(complex(argument))) return MathFunction.__new__(cls) def __init__(self, argument): MathFunction.__init__(self, "acos", argument)
[docs]@ufl_type() class Asin(MathFunction): __slots__ = () def __new__(cls, argument): if isinstance(argument, (RealValue, Zero)): return FloatValue(math.asin(float(argument))) if isinstance(argument, (ComplexValue)): return ComplexValue(cmath.asin(complex(argument))) return MathFunction.__new__(cls) def __init__(self, argument): MathFunction.__init__(self, "asin", argument)
[docs]@ufl_type() class Atan(MathFunction): __slots__ = () def __new__(cls, argument): if isinstance(argument, (RealValue, Zero)): return FloatValue(math.atan(float(argument))) if isinstance(argument, (ComplexValue)): return ComplexValue(cmath.atan(complex(argument))) return MathFunction.__new__(cls) def __init__(self, argument): MathFunction.__init__(self, "atan", argument)
[docs]@ufl_type(is_scalar=True, num_ops=2) class Atan2(Operator): __slots__ = () def __new__(cls, arg1, arg2): if isinstance(arg1, (RealValue, Zero)) and isinstance(arg2, (RealValue, Zero)): return FloatValue(math.atan2(float(arg1), float(arg2))) if isinstance(arg1, (ComplexValue)) or isinstance(arg2, (ComplexValue)): raise TypeError("Atan2 does not support complex numbers.") return Operator.__new__(cls) def __init__(self, arg1, arg2): Operator.__init__(self, (arg1, arg2)) if isinstance(arg1, (ComplexValue, complex)) or isinstance(arg2, (ComplexValue, complex)): raise TypeError("Atan2 does not support complex numbers.") if not is_true_ufl_scalar(arg1): error("Expecting scalar argument 1.") if not is_true_ufl_scalar(arg2): error("Expecting scalar argument 2.")
[docs] def evaluate(self, x, mapping, component, index_values): a = self.ufl_operands[0].evaluate(x, mapping, component, index_values) b = self.ufl_operands[1].evaluate(x, mapping, component, index_values) try: res = math.atan2(a, b) except TypeError: error('Atan2 does not support complex numbers.') except ValueError: warning('Value error in evaluation of function atan_2 with arguments %s, %s.' % (a, b)) raise return res
def __str__(self): return "atan_2(%s,%s)" % (self.ufl_operands[0], self.ufl_operands[1])
def _find_erf(): import math if hasattr(math, 'erf'): return math.erf import scipy.special if hasattr(scipy.special, 'erf'): return scipy.special.erf return None
[docs]@ufl_type() class Erf(MathFunction): __slots__ = () def __new__(cls, argument): if isinstance(argument, (RealValue, Zero)): erf = _find_erf() if erf is not None: return FloatValue(erf(float(argument))) if isinstance(argument, (ConstantValue)): erf = _find_erf() if erf is not None: return ComplexValue(erf(complex(argument))) return MathFunction.__new__(cls) def __init__(self, argument): MathFunction.__init__(self, "erf", argument)
[docs] def evaluate(self, x, mapping, component, index_values): a = self.ufl_operands[0].evaluate(x, mapping, component, index_values) erf = _find_erf() if erf is None: error("No python implementation of erf available on this system, cannot evaluate. Upgrade python or install scipy.") return erf(a)
[docs]@ufl_type(is_abstract=True, is_scalar=True, num_ops=2) class BesselFunction(Operator): "Base class for all bessel functions" __slots__ = ("_name", "_classname") def __init__(self, name, classname, nu, argument): if not is_true_ufl_scalar(nu): error("Expecting scalar nu.") if not is_true_ufl_scalar(argument): error("Expecting scalar argument.") # Use integer representation if suitable fnu = float(nu) inu = int(nu) if fnu == inu: nu = as_ufl(inu) else: nu = as_ufl(fnu) Operator.__init__(self, (nu, argument)) self._classname = classname self._name = name
[docs] def evaluate(self, x, mapping, component, index_values): a = self.ufl_operands[1].evaluate(x, mapping, component, index_values) try: import scipy.special except ImportError: error("You must have scipy installed to evaluate bessel functions in python.") name = self._name[-1] if isinstance(self.ufl_operands[0], IntValue): nu = int(self.ufl_operands[0]) functype = 'n' if name != 'i' else 'v' else: nu = self.ufl_operands[0].evaluate(x, mapping, component, index_values) functype = 'v' func = getattr(scipy.special, name + functype) return func(nu, a)
def __str__(self): return "%s(%s, %s)" % (self._name, self.ufl_operands[0], self.ufl_operands[1])
[docs]@ufl_type() class BesselJ(BesselFunction): __slots__ = () def __init__(self, nu, argument): BesselFunction.__init__(self, "cyl_bessel_j", "BesselJ", nu, argument)
[docs]@ufl_type() class BesselY(BesselFunction): __slots__ = () def __init__(self, nu, argument): BesselFunction.__init__(self, "cyl_bessel_y", "BesselY", nu, argument)
[docs]@ufl_type() class BesselI(BesselFunction): __slots__ = () def __init__(self, nu, argument): BesselFunction.__init__(self, "cyl_bessel_i", "BesselI", nu, argument)
[docs]@ufl_type() class BesselK(BesselFunction): __slots__ = () def __init__(self, nu, argument): BesselFunction.__init__(self, "cyl_bessel_k", "BesselK", nu, argument)