Merge branch 'firedrake/complex'

Add support for complex-valued elements

ChangeLog.rst

@@ -1,10 +1,18 @@
 Changelog
 =========
 
+2018.2.0
+--------
+
+- Add support for complex valued elements; complex mode
+  is chosen by ``compute_form_data(form, complex_mode=True)``
+  typically by a form compiler; otherwise UFL language is
+  agnostic to the choice of real/complex domain
+
 2018.1.0 (2018-06-14)
--------------
+---------------------
 
-- Remove python2 support. 
+- Remove python2 support
 
 2017.2.0 (2017-12-05)
 ---------------------

doc/sphinx/source/manual/form_language.rst

@@ -871,9 +871,9 @@ orthonormal basis for a Euclidean space):
 
   \mathbf{i}_i \cdot \mathbf{i}_j = \delta_{ij}
 
-where :math:`\delta_{ij}` is the Kronecker delta function.
-The dot product of higher order tensors follow from this, as illustrated
-with the following examples.
+where :math:`\delta_{ij}` is the Kronecker delta function.  The dot
+product of higher order tensors follow from this, as illustrated with
+the following examples.
 
 An example with two vectors
 
@@ -882,6 +882,7 @@ An example with two vectors
    \mathbf{v} \cdot \mathbf{u} = (v_i \mathbf{i}_i) \cdot (u_j \mathbf{i}_j)
         = v_i u_j (\mathbf{i}_i \cdot \mathbf{i}_j) = v_i u_j \delta_{ij} = v_i u_i
 
+
 An example with a tensor of rank two
 
 .. math::
@@ -913,20 +914,23 @@ The tensor rank of the product is rank(a)+rank(b)-2.
 ``inner``
 ---------
 
-The inner product is a contraction over all axes of a and b, that is
-the sum of all component-wise products.  The operands must have exactly the
+The inner product is a contraction over all axes of a and b, that is the
+sum of all component-wise products.  The operands must have exactly the
 same dimensions.  For two vectors it is equivalent to the dot product.
+Complex values are supported by UFL taking the complex conjugate of the
+second operand.  This has no impact if the values are real.
 
 If :math:`\mathbf{A}` and :math:`\mathbf{B}` are rank two tensors and
 :math:`\mathcal{C}` and :math:`\mathcal{D}` are rank 3 tensors
 their inner products are
 
 .. math::
-   \mathbf{A} : \mathbf{B}   &= A_{ij} B_{ij}
+   \mathbf{A} : \mathbf{B}   &= A_{ij} B^*_{ij}
    \\
-   \mathcal{C} : \mathcal{D} &= C_{ijk} D_{ijk}
+   \mathcal{C} : \mathcal{D} &= C_{ijk} D^*_{ijk}
 
-Using UFL notation, the following sets of declarations are equivalent::
+Using UFL notation, for real values, the following sets of declarations are
+equivalent::
 
   # Vectors
   f = dot(a, b)
@@ -941,6 +945,8 @@ Using UFL notation, the following sets of declarations are equivalent::
   f = inner(C, D)
   f = C[i,j,k]*D[i,j,k]
 
+Note that, in the UFL notation, `dot` and `inner` products are not equivalent
+for complex values.
 
 ``outer``
 ---------
@@ -957,18 +963,20 @@ The general definition of the outer product of two tensors
 
    \mathcal{C} \otimes \mathcal{D}
     =
-    C_{\iota^a_0 \ldots \iota^a_{r-1}} D_{\iota^b_0 \ldots\iota^b_{s-1}}
+    C^*_{\iota^a_0 \ldots \iota^a_{r-1}} D_{\iota^b_0 \ldots\iota^b_{s-1}}
     \mathbf{i}_{\iota^a_0}\otimes\cdots\otimes\mathbf{i}_{\iota^a_{r-2}}
     \otimes
     \mathbf{i}_{\iota^b_1} \otimes \cdots \otimes \mathbf{i}_{\iota^b_{s-1}}
 
+For consistency with the inner product, the complex conjugate is taken of the first operand.
+
 Some examples with vectors and matrices are easier to understand:
 
 .. math::
 
-   \mathbf{v} \otimes \mathbf{u} = v_i u_j \mathbf{i}_i \mathbf{i}_j, \\
-   \mathbf{v} \otimes \mathbf{B} = v_i B_{kl} \mathbf{i}_i \mathbf{i}_k \mathbf{i}_l, \\
-   \mathbf{A} \otimes \mathbf{B} = A_{ij} B_{kl} \mathbf{i}_i \mathbf{i}_j \mathbf{i}_k \mathbf{i}_l .
+   \mathbf{v} \otimes \mathbf{u} = v^*_i u_j \mathbf{i}_i \mathbf{i}_j, \\
+   \mathbf{v} \otimes \mathbf{B} = v^*_i B_{kl} \mathbf{i}_i \mathbf{i}_k \mathbf{i}_l, \\
+   \mathbf{A} \otimes \mathbf{B} = A^*_{ij} B_{kl} \mathbf{i}_i \mathbf{i}_j \mathbf{i}_k \mathbf{i}_l .
 
 The outer product of vectors is often written simply as
 
@@ -1423,6 +1431,56 @@ strain operator may be defined as follows::
   epsilon = lambda v: 0.5*(grad(v) + grad(v).T)
 
 
+Complex values
+==============
+
+UFL supports the definition of forms over either the real or the
+complex field.  Indeed, UFL does not explicitly define whether
+``Coefficient`` or ``Constant`` are real or complex. This is instead a
+matter for the form compiler to define. The complex-valued finite
+element spaces supported by UFL always have a real basis but complex
+coefficients. This means that ``Constant`` are ``Coefficient`` are
+complex-valued, but ``Argument`` is real-valued.
+
+Complex operators
+-----------------
+
+* ``conj(f)`` :: complex conjugate of ``f``.
+* ``imag(f)`` :: imaginary part of ``f``.
+* ``real(f)`` :: real part of ``f``.
+
+Sesquilinearity
+---------------
+
+``inner`` and ``outer`` are sesquilinear rather than linear
+when applied to complex values. Consequently, forms with two arguments
+are also sesquilinear in this case. UFL adopts the convention that
+inner products take the complex conjugate of the second operand. This
+is the usual convention in complex analysis but the reverse of the
+usual convention in physics.
+
+Complex values and conditionals
+-------------------------------
+
+Since the field of complex numbers does not admit a well order,
+complex expressions are not permissable as operands to ``lt``, ``gt``,
+``le``, or ``ge``. When compiling complex forms, the preprocessing
+stage of a compiler will attempt to prove that the relevant operands
+are real and will raise an exception if it is unable to do so. The
+user may always explicitly use ``real`` (or ``imag``) in order to
+ensure that the operand is real.
+
+
+Compiling real forms
+--------------------
+
+When the compiler treats a form as real, the preprocessing stage will
+discard all instances of ``conj`` and ``real`` in the form. Any
+instances of ``imag`` or complex literal constants will cause an
+exception.
+
+
+
 Form Transformations
 ====================
 
@@ -1435,9 +1493,9 @@ Replacing arguments of a Form
 -----------------------------
 
 The function ``replace`` lets you replace terminal objects with
-other values, using a mapping defined by a Python dicaionaryt. This can be
+other values, using a mapping defined by a Python dictionary. This can be
 used for example to replace a ``Coefficient`` with a fixed value for
-optimized runtime evaluation.
+optimized run-time evaluation.
 
 Example::
 
@@ -1619,7 +1677,7 @@ Assume in the following examples that
 
 The stiffness matrix can be computed from the functional
 :math:`\int_\Omega \nabla w : \nabla w \, dx`, by
-                              
+
 ::
 
   f = inner(grad(w), grad(w))/2 * dx

doc/sphinx/source/releases/next.rst

@@ -6,6 +6,8 @@ Changes in the next release
 Summary of changes
 ==================
 
+- Add support for complex valued elements
+
 .. note:: Developers should use this page to track and list changes
           during development. At the time of release, this page should
           be published (and renamed) to list the most important

scripts/ufl-convert

@@ -5,7 +5,8 @@ import subprocess
 import io
 import os
 import optparse
-from pprint import pprint
+
+from six import string_types
 
 from ufl.algorithms import tree_format, compute_form_data
 from ufl.formatting.ufl2dot import ufl2dot
@@ -37,6 +38,7 @@ def runcmd(cmd):
         print(output)
         sys.exit(-1)
 
+
 def write_file(filename, text):
     "Write text to a file and close it."
     with io.open(filename, "w", encoding="utf-8") as f:
@@ -45,22 +47,27 @@ def write_file(filename, text):
 
 # --- Option parsing
 
+
 usage = """Convert a .ufl file to some other format.
 
 Examples:
 
   ufl-convert -omydir -iyourdir -c -f -tpdf -s mass.ufl"""
 
+
 def opt(long, short, t, default, help):
     return optparse.make_option("--%s" % long, "-%s" % short, action="store", type=t, dest=long, default=default, help=help)
 
+
 option_list = [ \
     # Directories:
     opt("outputdir", "o", "str", "", "Output directory."),
     opt("inputdir",  "i", "str", "", "Input directory."),
     # Expression transformations:
     opt("labeling",  "l", "str", "repr", "Set to 'repr' or 'compact' for different naming of graph nodes."),
     opt("compile",   "c", "int", 0, "'Compile' forms: apply expression transformations like in a quadrature based form compilation. Only used for latex formatting."),
+    opt("preprocess", "p", "int", 0, "Set to 1 to return the preprocessed form. Only used for dot formatting."),
+    opt("complex", "z", "int", 0, "Set to 1 to for output containing complex nodes. Otherwise all nodes are assumed real. Only used for dot formatting."),
     # Output formats:
     opt("format",    "f", "str", "", "Rendering format (str, repr, tree, dot, latex, unicode)."),
     opt("filetype",  "t", "str", "", "Output file type (txt, py, dot, tex, ps, pdf, png)."),
@@ -92,10 +99,10 @@ for arg in args:
     # 1) Load forms
     ufl_data = load_ufl_file(uflfilename)
     forms = ufl_data.forms
-    #expressions = ufl_data.expressions # TODO: Allow rendering expressions without form stuff!
+    # expressions = ufl_data.expressions # TODO: Allow rendering expressions without form stuff!
 
     # Preprocess forms
-    form_datas = [compute_form_data(f) for f in forms]
+    form_datas = [compute_form_data(f, do_apply_geometry_lowering=True, complex_mode=options.complex) for f in forms]
 
     # 3) Render result
     format = options.format
@@ -131,7 +138,10 @@ for arg in args:
         data = []
         nodeoffset = 0
         for i, fd in enumerate(form_datas):
-            f = fd.original_form
+            if options.preprocess:
+                f = fd.preprocessed_form
+            else:
+                f = fd.original_form
             name = ufl_data.object_names.get(f, "form")
 
             begin = (i == 0)
@@ -175,7 +185,7 @@ for arg in args:
 
     # Conversions from tex:
     elif format == "tex":
-        texfile = os.path.join(options.outputdir, basename + ".tex") # TODO: Use a proper temp file?
+        texfile = os.path.join(options.outputdir, basename + ".tex")  # TODO: Use a proper temp file?
         write_file(texfile, rendered)
         if filetype == "pdf":
             flags = "-file-line-error-style -interaction=nonstopmode"
@@ -202,7 +212,7 @@ for arg in args:
 
     # That's all we know!
     else:
-        print("*** Error: Sorry, don't know how to render format '%s' for file type '%s'." \
-            % (format, filetype))
+        print("*** Error: Sorry, don't know how to render format '%s' for file type '%s'."
+              % (format, filetype))
         print("Please try another combination, perhaps -fdot -tpdf?")
         sys.exit(-1)

setup.cfg

@@ -1,3 +1,3 @@
 [flake8]
-ignore = E501,E226
+ignore = E501,E226,E741
 exclude = .git,__pycache__,doc/sphinx/source/conf.py,build,dist,test

test/test_arithmetic.py

@@ -4,16 +4,19 @@
 import pytest
 
 from ufl import *
-from ufl.classes import Division, FloatValue, IntValue
+from ufl.classes import Division, FloatValue, IntValue, ComplexValue
 
 
 def test_scalar_casting(self):
     f = as_ufl(2.0)
     r = as_ufl(4)
+    c = as_ufl(1 + 2j)
     self.assertIsInstance(f, FloatValue)
     self.assertIsInstance(r, IntValue)
+    self.assertIsInstance(c, ComplexValue)
     assert float(f) == 2.0
     assert int(r) == 4
+    assert complex(c) == 1.0 + 2.0j
 
 
 def test_ufl_float_division(self):

test/test_check_arities.py

@@ -3,6 +3,7 @@
 import pytest
 from ufl import *
 from ufl.algorithms.compute_form_data import compute_form_data
+from ufl.algorithms.check_arities import ArityMismatch
 
 
 def test_check_arities():
@@ -30,3 +31,23 @@ def test_check_arities():
     fd = compute_form_data(a)
 
     assert True
+
+
+def test_complex_arities():
+    cell = tetrahedron
+    D = Mesh(cell)
+    V = FunctionSpace(D, VectorElement("P", cell, 2))
+    v = TestFunction(V)
+    u = TrialFunction(V)
+
+    # Valid form.
+    F = inner(u, v) * dx
+    compute_form_data(F, complex_mode=True)
+    # Check that adjoint conjugates correctly
+    compute_form_data(adjoint(F), complex_mode=True)
+
+    with pytest.raises(ArityMismatch):
+        compute_form_data(inner(v, u) * dx, complex_mode=True)
+
+    with pytest.raises(ArityMismatch):
+        compute_form_data(inner(conj(v), u) * dx, complex_mode=True)