Hypergeometric2F1

ufl/__init__.py

@@ -380,6 +380,8 @@
     elem_mult,
     elem_op,
     elem_pow,
+    elliptic_E,
+    elliptic_K,
     eq,
     erf,
     exp,
@@ -388,6 +390,7 @@
     ge,
     grad,
     gt,
+    hyp2f1,
     imag,
     inner,
     inv,
@@ -595,6 +598,8 @@
     "elem_mult",
     "elem_op",
     "elem_pow",
+    "elliptic_E",
+    "elliptic_K",
     "energy_norm",
     "eq",
     "erf",
@@ -608,6 +613,7 @@
     "grad",
     "gt",
     "hexahedron",
+    "hyp2f1",
     "i",
     "identity_pullback",
     "imag",

ufl/algorithms/apply_derivatives.py

@@ -90,6 +90,7 @@
     Cosh,
     Erf,
     Exp,
+    Hypergeometric2F1,
     Ln,
     MathFunction,
     Sin,
@@ -112,6 +113,7 @@
     cosh,
     exp,
     facet_avg,
+    hyp2f1,
     ln,
     sign,
     sin,
@@ -583,6 +585,19 @@ def _(self, o: Expr, fp: Expr) -> Expr:
         (f,) = o.ufl_operands
         return fp * (2.0 / sqrt(pi) * exp(-(f**2)))
 
+    @process.register(Hypergeometric2F1)
+    @DAGTraverser.postorder
+    def _(self, o: Expr, ap: Expr, bp: Expr, cp: Expr, fp: Expr) -> Expr:
+        """Differentiate a hyp2f1."""
+        a, b, c, f = o.ufl_operands
+        if not all(nup is None or isinstance(nup, Zero) for nup in (ap, bp, cp)):
+            raise NotImplementedError(
+                "Differentiation of hypergeometric function w.r.t. (a, b, c) is not supported."
+            )
+
+        op = (a * b / c) * hyp2f1(a + 1, b + 1, c + 1, f)
+        return op * fp
+
     # --- Bessel functions
 
     @process.register(BesselJ)

ufl/mathfunctions.py

@@ -388,6 +388,42 @@ def evaluate(self, x, mapping, component, index_values):
         return math.erf(a)
 
 
+@ufl_type(is_scalar=True, num_ops=4)
+class Hypergeometric2F1(Operator):
+    """The Gaussian hypergeometric function 2F1."""
+
+    __slots__ = "_name"
+
+    def __init__(self, a, b, c, argument):
+        """Initialise."""
+        if not all(is_true_ufl_scalar(f) for f in (a, b, c)):
+            raise ValueError("Expecting scalar (a, b, c).")
+        if not is_true_ufl_scalar(argument):
+            raise ValueError("Expecting scalar argument.")
+
+        Operator.__init__(self, (a, b, c, argument))
+
+        self._name = "hyp2f1"
+
+    def evaluate(self, x, mapping, component, index_values):
+        """Evaluate."""
+        try:
+            import scipy.special  # type: ignore
+        except ImportError:
+            raise ValueError(
+                "You must have scipy installed to evaluate hypergeometric functions in python."
+            )
+        func = getattr(scipy.special, self._name)
+        return func(
+            *(arg.evaluate(x, mapping, component, index_values) for arg in self.ufl_operands)
+        )
+
+    def __str__(self):
+        """Format as a string."""
+        a, b, c, arg = self.ufl_operands
+        return f"{self._name}({a}, {b}, {c}, {arg})"
+
+
 @ufl_type(is_abstract=True, is_scalar=True, num_ops=2)
 class BesselFunction(Operator):
     """Base class for all bessel functions."""

ufl/operators.py

@@ -16,6 +16,7 @@
 
 import operator
 import warnings
+from math import pi
 
 from ufl import sobolevspace
 from ufl.algebra import Conj, Imag, Real
@@ -51,6 +52,7 @@
     Cosh,
     Erf,
     Exp,
+    Hypergeometric2F1,
     Ln,
     Sin,
     Sinh,
@@ -700,6 +702,35 @@ def erf(f):
     return _mathfunction(f, Erf)
 
 
+def hyp2f1(a, b, c, f):
+    """The Gaussian hypergeometric function 2F1."""
+    a = as_ufl(a)
+    b = as_ufl(b)
+    c = as_ufl(c)
+    f = as_ufl(f)
+    return Hypergeometric2F1(a, b, c, f)
+
+
+def elliptic_K(k):
+    r"""Complete elliptic integral of the first kind.
+
+    This function is defined as
+
+    .. math:: K(k) = \int_0^{\pi/2} [1 - k^2 sin(t)^2]^{-1/2} dt
+    """
+    return (pi / 2) * hyp2f1(1 / 2, 1 / 2, 1, k**2)
+
+
+def elliptic_E(k):
+    r"""Complete elliptic integral of the second kind.
+
+    This function is defined as
+
+    .. math:: E(k) = \int_0^{\pi/2} [1 - k^2 sin(t)^2]^{1/2} dt
+    """
+    return (pi / 2) * hyp2f1(1 / 2, -1 / 2, 1, k**2)
+
+
 def bessel_J(nu, f):
     """Cylindrical Bessel function of the first kind."""
     nu = as_ufl(nu)