diff -Nru ufl-2018.1.0/bitbucket-pipelines.yml ufl-2019.1.0/bitbucket-pipelines.yml --- ufl-2018.1.0/bitbucket-pipelines.yml 2018-06-14 14:20:02.000000000 +0000 +++ ufl-2019.1.0/bitbucket-pipelines.yml 1970-01-01 00:00:00.000000000 +0000 @@ -1,9 +0,0 @@ -image: quay.io/fenicsproject/pipelines - -pipelines: - default: - - step: - script: - - pip3 install . - - python3 -m flake8 ufl/ - - python3 -m pytest -v test/ diff -Nru ufl-2018.1.0/ChangeLog.rst ufl-2019.1.0/ChangeLog.rst --- ufl-2018.1.0/ChangeLog.rst 2018-06-14 14:20:02.000000000 +0000 +++ ufl-2019.1.0/ChangeLog.rst 2019-04-17 18:51:30.000000000 +0000 @@ -1,10 +1,20 @@ Changelog ========= +2019.1.0 (2019-04-17) +--------------------- + +- Remove scripts +- Remove LaTeX support (not functional) +- 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) --------------------- diff -Nru ufl-2018.1.0/debian/changelog ufl-2019.1.0/debian/changelog --- ufl-2018.1.0/debian/changelog 2019-01-30 08:22:32.000000000 +0000 +++ ufl-2019.1.0/debian/changelog 2019-07-31 04:00:45.000000000 +0000 @@ -1,3 +1,26 @@ +ufl (2019.1.0-2) unstable; urgency=medium + + * Standards-Version: 4.4.0 + * drop python-ufl + + -- Drew Parsons Wed, 31 Jul 2019 12:00:45 +0800 + +ufl (2019.1.0-1) experimental; urgency=medium + + * New upstream release (official). + * debhelper 12: Build-Depends: debhelper-compat (= 12) + + -- Drew Parsons Thu, 25 Apr 2019 10:41:16 +0800 + +ufl (2019.1.0~git20190220-1) experimental; urgency=medium + + * New upstream test release. + - applies patches inspect_API_update_09d5ec8.patch + inspect_API_update_94c0873.patch + and test_indexing_fixture_9d11804.patch + + -- Drew Parsons Wed, 06 Mar 2019 02:59:10 +0800 + ufl (2018.1.0-5) unstable; urgency=medium * add debian patch inspect_API_update_94c0873.patch diff -Nru ufl-2018.1.0/debian/compat ufl-2019.1.0/debian/compat --- ufl-2018.1.0/debian/compat 2019-01-30 08:22:32.000000000 +0000 +++ ufl-2019.1.0/debian/compat 1970-01-01 00:00:00.000000000 +0000 @@ -1 +0,0 @@ -12 diff -Nru ufl-2018.1.0/debian/control ufl-2019.1.0/debian/control --- ufl-2018.1.0/debian/control 2019-01-30 08:22:32.000000000 +0000 +++ ufl-2019.1.0/debian/control 2019-07-31 04:00:45.000000000 +0000 @@ -4,13 +4,13 @@ Maintainer: Debian Science Team Uploaders: Johannes Ring , Drew Parsons Build-Depends: - debhelper (>= 12), + debhelper-compat (= 12), dh-python, python3-all (>= 3.4), python3-setuptools, python3-numpy, python3-pytest -Standards-Version: 4.3.0 +Standards-Version: 4.4.0 Homepage: http://fenicsproject.org Vcs-Git: https://salsa.debian.org/science-team/fenics/ufl.git Vcs-Browser: https://salsa.debian.org/science-team/fenics/ufl @@ -36,23 +36,6 @@ . This package installs the library for Python 3. -Package: python-ufl -Architecture: all -Depends: - python3-ufl (>= 2018.1), ${misc:Depends} -Suggests: python-ufl-doc -Description: unified language for form-compilers - UFL (Unified Form Language) is a unified language for definition of - variational forms intended for finite element discretization. More - precisely, it defines a fixed interface for choosing finite element - spaces and defining expressions for weak forms in a notation close to - mathematical notation. The form compilers FFC and SyFi use UFL as - their end-user interface, producing UFC implementations as their - output. - . - This is a dummy package that depends on python3-ufl. - (UFL is no longer available for Python 2). - Package: python-ufl-doc Section: doc Architecture: all diff -Nru ufl-2018.1.0/debian/patches/inspect_API_update_09d5ec8.patch ufl-2019.1.0/debian/patches/inspect_API_update_09d5ec8.patch --- ufl-2018.1.0/debian/patches/inspect_API_update_09d5ec8.patch 2019-01-30 08:22:32.000000000 +0000 +++ ufl-2019.1.0/debian/patches/inspect_API_update_09d5ec8.patch 1970-01-01 00:00:00.000000000 +0000 @@ -1,36 +0,0 @@ -From 09d5ec8feac3dcdb0a3207d62ac9d8d17eb1d59c Mon Sep 17 00:00:00 2001 -From: Chris Richardson -Date: Fri, 14 Sep 2018 12:04:28 +0100 -Subject: [PATCH] Update syntax to suppress warning - ---- - ufl/corealg/multifunction.py | 6 +++--- - 1 file changed, 3 insertions(+), 3 deletions(-) - -diff --git a/ufl/corealg/multifunction.py b/ufl/corealg/multifunction.py -index 8fa615a..f901d8e 100644 ---- a/ufl/corealg/multifunction.py -+++ b/ufl/corealg/multifunction.py -@@ -20,7 +20,7 @@ - # - # Modified by Massimiliano Leoni, 2016 - --from inspect import getargspec -+import inspect - - from ufl.log import error - from ufl.core.expr import Expr -@@ -28,8 +28,8 @@ from ufl.core.expr import Expr - - def get_num_args(function): - "Return the number of arguments accepted by *function*." -- insp = getargspec(function) -- return len(insp[0]) + int(insp[1] is not None) -+ sig = inspect.signature(function) -+ return len(sig.parameters) + 1 - - - def memoized_handler(handler): --- -2.10.5 - diff -Nru ufl-2018.1.0/debian/patches/inspect_API_update_94c0873.patch ufl-2019.1.0/debian/patches/inspect_API_update_94c0873.patch --- ufl-2018.1.0/debian/patches/inspect_API_update_94c0873.patch 2019-01-30 08:22:32.000000000 +0000 +++ ufl-2019.1.0/debian/patches/inspect_API_update_94c0873.patch 1970-01-01 00:00:00.000000000 +0000 @@ -1,34 +0,0 @@ -From 94c087314152bdebae53584677c7c4f8392aa7aa Mon Sep 17 00:00:00 2001 -From: "Garth N. Wells" -Date: Fri, 28 Sep 2018 19:33:01 +0100 -Subject: [PATCH] Update inspect function. Fixes #99. - ---- - ufl/algorithms/transformer.py | 4 ++-- - 1 file changed, 2 insertions(+), 2 deletions(-) - -diff --git a/ufl/algorithms/transformer.py b/ufl/algorithms/transformer.py -index 7efbd0e..3c59ef4 100644 ---- a/ufl/algorithms/transformer.py -+++ b/ufl/algorithms/transformer.py -@@ -23,7 +23,7 @@ algorithms.""" - # - # Modified by Anders Logg, 2009-2010 - --from inspect import getargspec -+import inspect - from ufl.log import error - from ufl.classes import Variable, all_ufl_classes - from ufl.algorithms.map_integrands import map_integrands -@@ -31,7 +31,7 @@ from ufl.algorithms.map_integrands import map_integrands - - def is_post_handler(function): - "Is this a handler that expects transformed children as input?" -- insp = getargspec(function) -+ insp = inspect.getfullargspec(function) - num_args = len(insp[0]) + int(insp[1] is not None) - visit_children_first = num_args > 2 - return visit_children_first --- -2.10.5 - diff -Nru ufl-2018.1.0/debian/patches/series ufl-2019.1.0/debian/patches/series --- ufl-2018.1.0/debian/patches/series 2019-01-30 08:22:32.000000000 +0000 +++ ufl-2019.1.0/debian/patches/series 2019-07-31 04:00:45.000000000 +0000 @@ -1,3 +0,0 @@ -inspect_API_update_09d5ec8.patch -inspect_API_update_94c0873.patch -test_indexing_fixture_9d11804.patch diff -Nru ufl-2018.1.0/debian/patches/test_indexing_fixture_9d11804.patch ufl-2019.1.0/debian/patches/test_indexing_fixture_9d11804.patch --- ufl-2018.1.0/debian/patches/test_indexing_fixture_9d11804.patch 2019-01-30 08:22:32.000000000 +0000 +++ ufl-2019.1.0/debian/patches/test_indexing_fixture_9d11804.patch 1970-01-01 00:00:00.000000000 +0000 @@ -1,32 +0,0 @@ -From 9d11804abc2442619138d2f112a0e8a1e7e878b3 Mon Sep 17 00:00:00 2001 -From: Chris Richardson -Date: Thu, 15 Nov 2018 11:32:59 +0000 -Subject: [PATCH] Fix test fixture - ---- - test/test_indexing.py | 4 ++-- - 1 file changed, 2 insertions(+), 2 deletions(-) - -diff --git a/test/test_indexing.py b/test/test_indexing.py -index a53cdab..95d1029 100755 ---- a/test/test_indexing.py -+++ b/test/test_indexing.py -@@ -13,13 +13,13 @@ def x1(): - - @pytest.fixture - def x2(): -- x = x1() -+ x = SpatialCoordinate(triangle) - return outer(x, x) - - - @pytest.fixture - def x3(): -- x = x1() -+ x = SpatialCoordinate(triangle) - return outer(outer(x, x), x) - - --- -2.10.5 - diff -Nru ufl-2018.1.0/doc/sphinx/source/manual/command_line_utils.rst ufl-2019.1.0/doc/sphinx/source/manual/command_line_utils.rst --- ufl-2018.1.0/doc/sphinx/source/manual/command_line_utils.rst 2018-06-14 14:20:02.000000000 +0000 +++ ufl-2019.1.0/doc/sphinx/source/manual/command_line_utils.rst 1970-01-01 00:00:00.000000000 +0000 @@ -1,67 +0,0 @@ -********************* -Commandline utilities -********************* - -See installed version: ``ufl-analyse`` -====================================== - -Run - -:: - - # ufl-version - -to see the currently installed version of UFL printed to the terminal. - - -Validation and debugging: ``ufl-analyse`` -========================================= - -The command ``ufl-analyse`` loads all forms found in a ``.ufl`` -file, tries to discover any errors in them, and prints various kinds of -information about each form. Basic usage is - -:: - - # ufl-analyse myform.ufl - -For more information, type - -:: - - # ufl-analyse --help - -Note: This script is not well maintained, you will likely get more -useful information from your form compiler. - - -Formatting and visualization: ``ufl-convert`` -============================================= - -The command ``ufl-convert`` loads all forms found in a ``.ufl`` -file, compiles them into a different form or extracts some information -from them, and writes the result in a suitable file format. - -To try this tool, go to the ``demo/`` directory of the UFL source -tree. Some of the features to try are basic printing of ``str`` and -``repr`` string representations of each form:: - - # ufl-convert --format=str stiffness.ufl - # ufl-convert --format=repr stiffness.ufl - -compilation of forms to mathematical notation in LaTeX:: - - # ufl-convert --filetype=pdf --format=tex --show=1 stiffness.ufl - -LaTeX output of forms after processing with UFL compiler utilities:: - - # ufl-convert -tpdf -ftex -s1 --compile=1 stiffness.ufl - -and visualization of expression trees using graphviz via compilation of -forms to the dot format:: - - # ufl-convert -tpdf -fdot -s1 stiffness.ufl - -Type ``ufl-convert --help`` for more details. - - diff -Nru ufl-2018.1.0/doc/sphinx/source/manual/form_language.rst ufl-2019.1.0/doc/sphinx/source/manual/form_language.rst --- ufl-2018.1.0/doc/sphinx/source/manual/form_language.rst 2018-06-14 14:20:02.000000000 +0000 +++ ufl-2019.1.0/doc/sphinx/source/manual/form_language.rst 2019-04-17 18:51:30.000000000 +0000 @@ -871,9 +871,9 @@ \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 @@ \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 @@ ``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 @@ 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 @@ \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 @@ 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 @@ ----------------------------- 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 @@ 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 diff -Nru ufl-2018.1.0/doc/sphinx/source/manual/internal_representation.rst ufl-2019.1.0/doc/sphinx/source/manual/internal_representation.rst --- ufl-2018.1.0/doc/sphinx/source/manual/internal_representation.rst 2018-06-14 14:20:02.000000000 +0000 +++ ufl-2019.1.0/doc/sphinx/source/manual/internal_representation.rst 2019-04-17 18:51:30.000000000 +0000 @@ -178,10 +178,10 @@ you can get away with implementing a new operator as a combination of existing ones, that is the easiest route. The reason is that only some of the properties of an operator is represented by the ``Expr`` -subclass. Other properties are part of the various algorithms in -UFL. One example is derivatives, which are defined in the differentiation -algorithm, and how to render a type to the ``LaTeX`` or dot formats. These -properties could be merged into the class hierarchy, but other properties -like how to map a UFL type to some ``ffc`` or ``dolfin`` type -cannot be part of UFL. So before adding a new class, consider that doing -so may require changes in multiple algorithms and even other projects. +subclass. Other properties are part of the various algorithms in UFL. +One example is derivatives, which are defined in the differentiation +algorithm, and how to render a type to the dot formats. These properties +could be merged into the class hierarchy, but other properties like how +to map a UFL type to some ``ffc`` or ``dolfin`` type cannot be part of +UFL. So before adding a new class, consider that doing so may require +changes in multiple algorithms and even other projects. diff -Nru ufl-2018.1.0/doc/sphinx/source/manual/user_manual.rst ufl-2019.1.0/doc/sphinx/source/manual/user_manual.rst --- ufl-2018.1.0/doc/sphinx/source/manual/user_manual.rst 2018-06-14 14:20:02.000000000 +0000 +++ ufl-2019.1.0/doc/sphinx/source/manual/user_manual.rst 2019-04-17 18:51:30.000000000 +0000 @@ -14,4 +14,3 @@ examples internal_representation algorithms - command_line_utils diff -Nru ufl-2018.1.0/doc/sphinx/source/releases/next.rst ufl-2019.1.0/doc/sphinx/source/releases/next.rst --- ufl-2018.1.0/doc/sphinx/source/releases/next.rst 2018-06-14 14:20:02.000000000 +0000 +++ ufl-2019.1.0/doc/sphinx/source/releases/next.rst 2019-04-17 18:51:30.000000000 +0000 @@ -6,7 +6,7 @@ Summary of changes ================== -- Remove Python 2 support. Now requires Python 3.5 or above. +- 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 diff -Nru ufl-2018.1.0/doc/sphinx/source/releases/v2019.1.0.rst ufl-2019.1.0/doc/sphinx/source/releases/v2019.1.0.rst --- ufl-2018.1.0/doc/sphinx/source/releases/v2019.1.0.rst 1970-01-01 00:00:00.000000000 +0000 +++ ufl-2019.1.0/doc/sphinx/source/releases/v2019.1.0.rst 2019-04-17 18:51:30.000000000 +0000 @@ -0,0 +1,18 @@ +=========================== +Changes in version 2019.1.0 +=========================== + +Summary of changes +================== + +- Add support for complex valued elements +- Remove LaTeX support (not functional) +- Remove scripts + +Detailed changes +================ + +- 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 diff -Nru ufl-2018.1.0/doc/sphinx/source/releases.rst ufl-2019.1.0/doc/sphinx/source/releases.rst --- ufl-2018.1.0/doc/sphinx/source/releases.rst 2018-06-14 14:20:02.000000000 +0000 +++ ufl-2019.1.0/doc/sphinx/source/releases.rst 2019-04-17 18:51:30.000000000 +0000 @@ -9,6 +9,7 @@ :maxdepth: 2 releases/next + releases/v2019.1.0 releases/v2018.1.0 releases/v2017.2.0 releases/v2017.1.0.post1 diff -Nru ufl-2018.1.0/.mailmap ufl-2019.1.0/.mailmap --- ufl-2018.1.0/.mailmap 2018-06-14 14:20:02.000000000 +0000 +++ ufl-2019.1.0/.mailmap 1970-01-01 00:00:00.000000000 +0000 @@ -1,132 +0,0 @@ -Anders Logg -Anders Logg -Anders Logg -Anders Logg -Anders Logg logg -Anders Logg fenics -Anders Logg hg -Anders Logg root -Anders Logg -Anders Logg -Anders Logg -Anders Logg -Anders Logg -Anders Logg anders -Anders Logg logg -Anders Logg logg -Anders Logg logg -Anders Logg logg -Andy R. Terrel -Andy R. Terrel -Andy R. Terrel aterrel -Andy R. Terrel -Benjamin Kehlet -Benjamin Kehlet -Benjamin Kehlet -Benjamin Kehlet -Benjamin Kehlet -Chris Richardson -Chris Richardson chris -Chris Richardson Chris Richardson -Chris Richardson root -Corrado Maurini -Corrado Maurini -Dag Lindbo -Dag Lindbo dag -Dag Lindbo dag -David Ham -David Ham david.ham@imperial.ac.uk <> -Evan Lezar -Evan Lezar -Evan Lezar elezar -Evan Lezar elezar -Fredrik Valdmanis -Fredrik Valdmanis Fredrik Valdmanis -Garth N. Wells -Garth N. Wells root -Garth N. Wells garth -Garth N. Wells -Garth N. Wells -Garth N. Wells -Garth N. Wells -Garth N. Wells -Garth N. Wells -Garth N. Wells -Garth N. Wells -Garth N. Wells -Garth N. Wells -Garth N. Wells gnw20@cam.ac.uk <> -gideonsimpson -Gustav Magnus Vikström -Gustav Magnus Vikström -Gustav Magnus Vikström -Gustav Magnus Vikström -Harish Narayanan -Harish Narayanan Harish Narayanan -Johan Hoffman -Johan Hoffman hoffman -Johan Hoffman -Johan Hoffman -Ilmar Wilbers -Ilmar Wilbers -Ilmar Wilbers -Jack S. Hale -Jack S. Hale -Johan Hake -Johan Hake -Johan Hake -Johan Jansson -Johan Jansson johan -Johan Jansson johan -Johan Jansson johanjan -Johan Jansson johanjan -Johan Jansson Johan Jansson -Johannes Ring -Johannes Ring -Kent-Andre Mardal -Kent-Andre Mardal -Kristian B. Ølgaard -Kristian B. Ølgaard -Kristian B. Ølgaard -Magnus Vikstrøm -Marco Morandini -Marco Morandini -Marie E. Rognes -Marie E. Rognes -Marie E. Rognes -Marie E. Rognes Marie E. Rognes (meg@simula.no) -Marie E. Rognes meg@simula.no <> -Martin Sandve Alnæs -Martin Sandve Alnæs -Martin Sandve Alnæs -Michele Zaffalon -Mikael Mortensen -Mikael Mortensen -Nate Sime -Nate Sime -Nuno Lopes -Nuno Lopes N.Lopes -Patrick Farrell -Patrick Farrell -Patrick Farrell -Quang Ha -Kristoffer Selim -Kristoffer Selim -Simon Funke -Simon Funke -Simon Funke -Solveig Bruvoll -Solveig Masvie -Steven Vandekerckhove -Steven Vandekerckhove -stockli -stockli -Tianyi Li -Steffen Müthing -Steffen Müthing Steffen Müthing steffen.muething@ipvs.uni-stuttgart.de <> -Miklós Homolya -Åsmund Ødegård -Åsmund Ødegård -Ola Skavhaug -Andre Massing -Andrew McRae diff -Nru ufl-2018.1.0/release.conf ufl-2019.1.0/release.conf --- ufl-2018.1.0/release.conf 2018-06-14 14:20:02.000000000 +0000 +++ ufl-2019.1.0/release.conf 1970-01-01 00:00:00.000000000 +0000 @@ -1,8 +0,0 @@ -# Configuration file for fenics-release - -PACKAGE="ufl" -BRANCH="master" -FILES="ChangeLog.rst \ - setup.py \ - doc/sphinx/source/releases/next.rst \ - doc/sphinx/source/releases.rst" diff -Nru ufl-2018.1.0/scripts/ufl2py ufl-2019.1.0/scripts/ufl2py --- ufl-2018.1.0/scripts/ufl2py 2018-06-14 14:20:02.000000000 +0000 +++ ufl-2019.1.0/scripts/ufl2py 1970-01-01 00:00:00.000000000 +0000 @@ -1,56 +0,0 @@ -#!/usr/bin/env python - -import os, sys, optparse -from ufl.algorithms import validate_form, ufl2latex, tree_format -from ufl.algorithms.formfiles import read_ufl_file - -# Get commandline options - -usage = """Convert a .ufl file to an executable .py file for debugging. - -Example: - - ufl2py Poisson.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 = [] -# opt("quiet", "q", "int", 1, "Do not print form information to screen."), -# opt("write", "w", "int", 0, "Write form information to file."), -# ] - -parser = optparse.OptionParser(usage=usage, option_list=option_list) -args = sys.argv[1:] -(options, args) = parser.parse_args(args=args) - -if not args: - print("Missing files!") - print() - parser.print_usage() - sys.exit(-1) -filenames = args - -header = """#!/usr/bin/env python -# -*- coding: utf-8 -*- -from ufl import * -set_level(DEBUG) -""" - -footer = "" - -# Handle each form file separately -for filename in filenames: - if not filename.endswith(".ufl"): - print("Warning: Filename '%s' doesn't end with .ufl." % filename) - - # Read code - fcode = read_ufl_file(filename) - code = header + fcode + footer - - # Dump code to python file - basename = os.path.splitext(os.path.basename(filename))[0] - basename = "%s_debug" % basename - pyname = "%s.py" % basename - with file(pyname, "w") as f: - f.write(code) diff -Nru ufl-2018.1.0/scripts/ufl-analyse ufl-2019.1.0/scripts/ufl-analyse --- ufl-2018.1.0/scripts/ufl-analyse 2018-06-14 14:20:02.000000000 +0000 +++ ufl-2019.1.0/scripts/ufl-analyse 1970-01-01 00:00:00.000000000 +0000 @@ -1,101 +0,0 @@ -#!/usr/bin/env python - -__authors__ = "Martin Sandve Alnes" -__date__ = "2008-05-09" - -# Modified by Anders Logg, 2009. -# Last changed: 2015-01-05 - -import io -import sys, optparse -from ufl.log import warning -from ufl.algorithms import load_ufl_file, validate_form, ufl2latex, tree_format - -# Get commandline options - -usage = """Analyse a .ufl file to find errors. -Optionally write information about -the forms for further inspection. - -Examples: - - ufl-analyse --quiet=0 --write=1 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 = [ \ - opt("quiet", "q", "int", 1, "Do not print form information to screen."), - opt("write", "w", "int", 0, "Write form information to file."), - ] - -parser = optparse.OptionParser(usage=usage, option_list=option_list) -args = sys.argv[1:] -(options, args) = parser.parse_args(args=args) - -if not args: - print("Missing files!") - print() - parser.print_usage() - sys.exit(-1) -filenames = args - -write_file = options.write -quiet = options.quiet - -# Handle each form file separately -for filename in filenames: - - # Check file suffix - if not filename.endswith(".ufl"): - warning("Filename '%s' does not end with .ufl." % filename) - - # Load form file, which triggers many consistency - # checks while the form is being built - print("Loading form file '%s'" % filename) - try: - # TODO: Forms that fail will usually fail inside this, - # which doesn't produce any log... - # Perhaps we should pass a log file to load_forms? - data = load_ufl_file(filename) - forms = data.forms - except: - print("Failed to load form file.") - raise - - outputfilename = filename + ".log" - if write_file: - outputfile = io.open(outputfilename, "w", encoding="utf-8") - - def write(*items): - text = " ".join(str(s) for s in items) - if write_file: - outputfile.write(text) - outputfile.flush() - if not quiet: - print(text) - - # Analyse each form separately - for form in forms: - - # Validate form - validate_form(form) - - # Compute form metadata and extract preprocessed form - form_data = form.compute_form_data() - preprocessed_form = form_data.preprocessed_form - - # Print form data - write("\nForm data:\n", str(form_data)) - - # Print different representations - write("\n\nForm pretty-print (original):\n", str(form)) - write("\n\nForm pretty-print (preprocessed):\n", str(preprocessed_form)) - write("\n\nForm representation (original):\n", repr(form)) - write("\n\nForm representation (preprocessed):\n", repr(preprocessed_form)) - write("\n\nForm tree formatting (original):\n", tree_format(form)) - write("\n\nForm tree formatting (preprocessed):\n", tree_format(preprocessed_form)) - write("\n\nForm LaTeX code (preprocessed):\n", ufl2latex(form)) - - if write_file: - outputfile.close() diff -Nru ufl-2018.1.0/scripts/ufl-convert ufl-2019.1.0/scripts/ufl-convert --- ufl-2018.1.0/scripts/ufl-convert 2018-06-14 14:20:02.000000000 +0000 +++ ufl-2019.1.0/scripts/ufl-convert 1970-01-01 00:00:00.000000000 +0000 @@ -1,208 +0,0 @@ -#!/usr/bin/env python - -import sys -import subprocess -import io -import os -import optparse -from pprint import pprint - -from ufl.algorithms import tree_format, compute_form_data -from ufl.formatting.ufl2dot import ufl2dot -from ufl.formatting.ufl2latex import forms2latexdocument -from ufl.formatting.ufl2unicode import form2unicode -from ufl.algorithms.formfiles import load_ufl_file - -# --- Utilities - - -def get_status_output(cmd, input=None, cwd=None, env=None): - """Replacement for commands.getstatusoutput which does not work on Windows (or Python 3).""" - if isinstance(cmd, string_types): - cmd = cmd.strip().split() - pipe = subprocess.Popen(cmd, shell=False, cwd=cwd, env=env, - stdout=subprocess.PIPE, stderr=subprocess.STDOUT) - (output, errout) = pipe.communicate(input=input) - assert not errout - status = pipe.returncode - if isinstance(output, bytes): - output = output.decode('utf-8') - return (status, output) - - -def runcmd(cmd): - status, output = get_status_output(cmd) - if status != 0: - print("*** Error:") - 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: - f.write(text) - print("Wrote file %s" % filename) - -# --- 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."), - # 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)."), - ] - -parser = optparse.OptionParser(usage=usage, option_list=option_list) -args = sys.argv[1:] -(options, args) = parser.parse_args(args=args) - -if not args: - print("Missing files!") - print() - parser.print_usage() - sys.exit(-1) - - -# --- Handle each file - -for arg in args: - - # 0) Get and check filename - uflfilename = os.path.join(options.inputdir, arg) - path, name = os.path.split(uflfilename) - basename, ext = os.path.splitext(name) - if ext != ".ufl": - print("Expecting a .ufl file, not ", uflfilename) - sys.exit(-1) - - # 1) Load forms - ufl_data = load_ufl_file(uflfilename) - forms = ufl_data.forms - #expressions = ufl_data.expressions # TODO: Allow rendering expressions without form stuff! - - # Preprocess forms - form_datas = [compute_form_data(f) for f in forms] - - # 3) Render result - format = options.format - - # Make format string conform - if format == "latex": - format = "tex" - if format == "tex": - rendered = forms2latexdocument(forms, uflfilename, compile=options.compile) - elif format == "unicode": - data = [] - for form, form_data in zip(forms, form_datas): - tmp = form2unicode(form, form_data) - data.append(tmp) - rendered = "\n\n".join(data) - elif format in ("str", "repr", "tree"): - data = [] - for i, fd in enumerate(form_datas): - f = fd.original_form - name = ufl_data.object_names.get(f, "form") - - if format == "str": - s = str(f) - elif format == "repr": - s = repr(f) - elif format == "tree": - s = tree_format(f) - tmp = "Form %s:\n%s\n" % (name, s) - - data.append(tmp) - rendered = "\n\n".join(data) - elif format == "dot": - data = [] - nodeoffset = 0 - for i, fd in enumerate(form_datas): - f = fd.original_form - name = ufl_data.object_names.get(f, "form") - - begin = (i == 0) - end = (i == len(forms) - 1) - dot, nodeoffset = ufl2dot(f, name, nodeoffset, begin, end, - options.labeling, ufl_data.object_names) - tmp = "/* Form %s: */\n%s\n" % (name, dot) - data.append(tmp) - rendered = "\n\n".join(data) - else: - print("Unknown rendering format ", format) - sys.exit(-1) - - # 4) Convert file format - filetype = options.filetype - - # Default filetypes: - if not filetype: - if format == "str": - filetype = "str" - elif format == "repr": - filetype = "repr" - elif format == "tree": - filetype = "tree" - elif format == "dot": - filetype = "dot" - elif format == "tex": - filetype = "tex" - elif format == "unicode": - filetype = "txt" - - # Guess that the filetype is the ext, usually the case - ext = filetype - if ext and not ext.startswith("."): - ext = "." + ext - outputfilename = os.path.join(options.outputdir, basename + ext) - - # Pure text files: - if filetype == "txt" or filetype == format: - write_file(outputfilename, rendered) - - # Conversions from tex: - elif format == "tex": - 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" - cmd = "pdflatex %s '%s'" % (flags, texfile) - runcmd(cmd) - else: - print("Unknown format and filetype combination:", format, filetype) - sys.exit(-1) - - # Conversions from dot: - elif format == "dot": - tempfile = os.path.join(options.outputdir, basename + ".dot") # TODO: Use a proper temp file? - write_file(tempfile, rendered) - if filetype in ("png", "ps", "svg", "gif", "dia", "imap", "cmapx"): # taken from "man dot" - runcmd("dot -T%s -o'%s' '%s'" % (filetype, outputfilename, tempfile)) - elif filetype == "pdf": - psfilename = os.path.join(options.outputdir, basename + ".ps") - pdffilename = os.path.join(options.outputdir, basename + ".pdf") - runcmd("dot -T%s -o'%s' '%s'" % (filetype, psfilename, tempfile)) - runcmd("ps2pdf '%s' '%s'" % (psfilename, pdffilename)) - else: - print("Unknown format and filetype combination:", format, filetype) - sys.exit(-1) - - # That's all we know! - else: - 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) diff -Nru ufl-2018.1.0/scripts/ufl-version ufl-2019.1.0/scripts/ufl-version --- ufl-2018.1.0/scripts/ufl-version 2018-06-14 14:20:02.000000000 +0000 +++ ufl-2019.1.0/scripts/ufl-version 1970-01-01 00:00:00.000000000 +0000 @@ -1,23 +0,0 @@ -#!/usr/bin/env python -# -*- coding: utf-8 -*- -# Copyright (C) 2011-2016 Anders Logg. -# -# This file is part of UFL. -# -# UFL is free software: you can redistribute it and/or modify -# it under the terms of the GNU Lesser General Public License as published by -# the Free Software Foundation, either version 3 of the License, or -# (at your option) any later version. -# -# UFL is distributed in the hope that it will be useful, -# but WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -# GNU Lesser General Public License for more details. -# -# You should have received a copy of the GNU Lesser General Public License -# along with UFL. If not, see . - -"""This is a simple script that just prints the UFL version number.""" - -from ufl import __version__ -print(__version__) diff -Nru ufl-2018.1.0/setup.cfg ufl-2019.1.0/setup.cfg --- ufl-2018.1.0/setup.cfg 2018-06-14 14:20:02.000000000 +0000 +++ ufl-2019.1.0/setup.cfg 2019-04-17 18:51:30.000000000 +0000 @@ -1,3 +1,3 @@ [flake8] -ignore = E501,E226 +ignore = E501, W504 exclude = .git,__pycache__,doc/sphinx/source/conf.py,build,dist,test diff -Nru ufl-2018.1.0/setup.py ufl-2019.1.0/setup.py --- ufl-2018.1.0/setup.py 2018-06-14 14:20:02.000000000 +0000 +++ ufl-2019.1.0/setup.py 2019-04-17 18:51:30.000000000 +0000 @@ -1,9 +1,7 @@ # -*- coding: utf-8 -*- from setuptools import setup -from os.path import join, split import sys -import platform module_name = "ufl" @@ -11,30 +9,12 @@ print("Python 3.5 or higher required, please upgrade.") sys.exit(1) -version = "2018.1.0" +version = "2019.1.0" -url = "https://bitbucket.org/fenics-project/%s/" % module_name +url = "https://bitbucket.org/fenics-project/{}/".format(module_name) tarball = None if 'dev' not in version: - tarball = url + "downloads/fenics-%s-%s.tar.gz" % (module_name, version) - -script_names = ("ufl-analyse", "ufl-convert", "ufl-version", "ufl2py") - -scripts = [join("scripts", script) for script in script_names] -man_files = [join("doc", "man", "man1", "%s.1.gz" % (script,)) for script in script_names] -data_files = [(join("share", "man", "man1"), man_files)] - -if platform.system() == "Windows" or "bdist_wininst" in sys.argv: - # In the Windows command prompt we can't execute Python scripts - # without a .py extension. A solution is to create batch files - # that runs the different scripts. - batch_files = [] - for script in scripts: - batch_file = script + ".bat" - with open(batch_file, "w") as f: - f.write(sys.executable + ' "%%~dp0\%s" %%*' % split(script)[1]) - batch_files.append(batch_file) - scripts.extend(batch_files) + tarball = url + "downloads/fenics-{}-{}.tar.gz".format(module_name, version) CLASSIFIERS = """\ Development Status :: 5 - Production/Stable @@ -53,25 +33,23 @@ Topic :: Software Development :: Libraries :: Python Modules """ -setup(name="fenics-ufl", - version=version, - description="Unified Form Language", - author="Martin Sandve Alnæs, Anders Logg", - author_email="fenics-dev@googlegroups.com", - url=url, - download_url=tarball, - classifiers=[_f for _f in CLASSIFIERS.split('\n') if _f], - scripts=scripts, - packages=[ - "ufl", - "ufl.utils", - "ufl.finiteelement", - "ufl.core", - "ufl.corealg", - "ufl.algorithms", - "ufl.formatting", - ], - package_dir={"ufl": "ufl"}, - install_requires=["numpy"], - data_files=data_files - ) +setup( + name="fenics-ufl", + version=version, + description="Unified Form Language", + author="Martin Sandve Alnæs, Anders Logg", + author_email="fenics-dev@googlegroups.com", + url=url, + download_url=tarball, + classifiers=[_f for _f in CLASSIFIERS.split('\n') if _f], + packages=[ + "ufl", + "ufl.utils", + "ufl.finiteelement", + "ufl.core", + "ufl.corealg", + "ufl.algorithms", + "ufl.formatting", + ], + package_dir={"ufl": "ufl"}, + install_requires=["numpy"]) diff -Nru ufl-2018.1.0/test/clean.sh ufl-2019.1.0/test/clean.sh --- ufl-2018.1.0/test/clean.sh 2018-06-14 14:20:02.000000000 +0000 +++ ufl-2019.1.0/test/clean.sh 1970-01-01 00:00:00.000000000 +0000 @@ -1,18 +0,0 @@ -#!/bin/bash - -# ufl-analyse output -rm -f *_debug.py - -# python compiled files -rm -f *.pyc - -# logs -rm -f *.log - -# temp files from emacs -rm -f *~ - -# temp files from py.test -rm -rf __pycache__ - - diff -Nru ufl-2018.1.0/test/pytest.ini ufl-2019.1.0/test/pytest.ini --- ufl-2018.1.0/test/pytest.ini 2018-06-14 14:20:02.000000000 +0000 +++ ufl-2019.1.0/test/pytest.ini 1970-01-01 00:00:00.000000000 +0000 @@ -1,3 +0,0 @@ -[pytest] -minversion = 2.4 -looponfailroots = . ../ufl diff -Nru ufl-2018.1.0/test/README ufl-2019.1.0/test/README --- ufl-2018.1.0/test/README 2018-06-14 14:20:02.000000000 +0000 +++ ufl-2019.1.0/test/README 1970-01-01 00:00:00.000000000 +0000 @@ -1,20 +0,0 @@ -To run tests, you need the py.test module. - -Just run - - cd /test - py.test - -Or to run a single test file - - py.test test_literals.py - -To run these tests from within the source tree without -needing to install the UFL Python module, update your -PYTHONPATH and PATH by running - - source sourceme.sh - -in a bash shell or equivalent. If on other OSes, you -must set the paths whichever way your OS requires. - diff -Nru ufl-2018.1.0/test/test_apply_function_pullbacks.py ufl-2019.1.0/test/test_apply_function_pullbacks.py --- ufl-2018.1.0/test/test_apply_function_pullbacks.py 2018-06-14 14:20:02.000000000 +0000 +++ ufl-2019.1.0/test/test_apply_function_pullbacks.py 2019-04-17 18:51:30.000000000 +0000 @@ -41,6 +41,8 @@ Vc = FiniteElement("N1curl", cell, 1) T = TensorElement("CG", cell, 1) S = TensorElement("CG", cell, 1, symmetry=True) + COV2T = FiniteElement("Regge", cell, 0) # (0, 2)-symmetric tensors + CONTRA2T = FiniteElement("HHJ", cell, 0) # (2, 0)-symmetric tensors Um = U*U Vm = U*V @@ -59,6 +61,8 @@ vc = Coefficient(Vc) t = Coefficient(T) s = Coefficient(S) + cov2t = Coefficient(COV2T) + contra2t = Coefficient(CONTRA2T) um = Coefficient(Um) vm = Coefficient(Vm) @@ -77,6 +81,8 @@ rvc = ReferenceValue(vc) rt = ReferenceValue(t) rs = ReferenceValue(s) + rcov2t = ReferenceValue(cov2t) + rcontra2t = ReferenceValue(contra2t) rum = ReferenceValue(um) rvm = ReferenceValue(vm) @@ -117,6 +123,9 @@ s: as_tensor([[rs[0], rs[1], rs[2]], [rs[1], rs[3], rs[4]], [rs[2], rs[4], rs[5]]]), + cov2t: as_tensor(Jinv[k, i] * rcov2t[k, l] * Jinv[l, j], (i, j)), + contra2t: as_tensor((1.0 / detJ) * (1.0 / detJ) + * J[i, k] * rcontra2t[k, l] * J[j, l], (i, j)), # Mixed elements become a bit more complicated um: rum, vm: rvm, @@ -201,6 +210,8 @@ check_single_function_pullback(vc, mappings) check_single_function_pullback(t, mappings) check_single_function_pullback(s, mappings) + check_single_function_pullback(cov2t, mappings) + check_single_function_pullback(contra2t, mappings) # Check functions of various elements inside a mixed context check_single_function_pullback(um, mappings) diff -Nru ufl-2018.1.0/test/test_arithmetic.py ufl-2019.1.0/test/test_arithmetic.py --- ufl-2018.1.0/test/test_arithmetic.py 2018-06-14 14:20:02.000000000 +0000 +++ ufl-2019.1.0/test/test_arithmetic.py 2019-04-17 18:51:30.000000000 +0000 @@ -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): diff -Nru ufl-2018.1.0/test/test_change_to_reference_frame.py ufl-2019.1.0/test/test_change_to_reference_frame.py --- ufl-2018.1.0/test/test_change_to_reference_frame.py 2018-06-14 14:20:02.000000000 +0000 +++ ufl-2019.1.0/test/test_change_to_reference_frame.py 2019-04-17 18:51:30.000000000 +0000 @@ -51,8 +51,8 @@ integrand = change_to_reference_frame(integral.integrand()) # Compute and apply integration scaling factor - scale = compute_integrand_scaling_factor(integral.ufl_domain(), - integral.integral_type()) + scale, degree = compute_integrand_scaling_factor(integral.ufl_domain(), + integral.integral_type()) return integral.reconstruct(integrand * scale) # TODO: , reference=True) diff -Nru ufl-2018.1.0/test/test_check_arities.py ufl-2019.1.0/test/test_check_arities.py --- ufl-2018.1.0/test/test_check_arities.py 2018-06-14 14:20:02.000000000 +0000 +++ ufl-2019.1.0/test/test_check_arities.py 2019-04-17 18:51:30.000000000 +0000 @@ -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 @@ 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) diff -Nru ufl-2018.1.0/test/test_complex.py ufl-2019.1.0/test/test_complex.py --- ufl-2018.1.0/test/test_complex.py 1970-01-01 00:00:00.000000000 +0000 +++ ufl-2019.1.0/test/test_complex.py 2019-04-17 18:51:30.000000000 +0000 @@ -0,0 +1,171 @@ +#!/usr/bin/env py.test +# -*- coding: utf-8 -*- +import pytest +import cmath +import ufl +from ufl.constantvalue import Zero, ComplexValue +from ufl.algebra import Conj, Real, Imag +from ufl.algorithms.apply_algebra_lowering import apply_algebra_lowering +from ufl.algorithms.remove_complex_nodes import remove_complex_nodes +from ufl.algorithms import estimate_total_polynomial_degree +from ufl.algorithms.comparison_checker import do_comparison_check, ComplexComparisonError +from ufl.algorithms.formtransformations import compute_form_adjoint +from ufl import TestFunction, TrialFunction, triangle, FiniteElement, \ + as_ufl, inner, grad, dx, dot, outer, conj, sqrt, sin, cosh, \ + atan, ln, exp, as_tensor, real, imag, conditional, \ + min_value, max_value, gt, lt, cos, ge, le, Coefficient + + +def test_conj(self): + z1 = ComplexValue(1+2j) + z2 = ComplexValue(1-2j) + + assert z1 == Conj(z2) + assert z2 == Conj(z1) + + +def test_real(self): + z0 = Zero() + z1 = as_ufl(1.0) + z2 = ComplexValue(1j) + z3 = ComplexValue(1+1j) + assert Real(z1) == z1 + assert Real(z3) == z1 + assert Real(z2) == z0 + + +def test_imag(self): + z0 = Zero() + z1 = as_ufl(1.0) + z2 = as_ufl(1j) + z3 = ComplexValue(1+1j) + + assert Imag(z2) == z1 + assert Imag(z3) == z1 + assert Imag(z1) == z0 + + +def test_compute_form_adjoint(self): + cell = triangle + element = FiniteElement('Lagrange', cell, 1) + + u = TrialFunction(element) + v = TestFunction(element) + + a = inner(grad(u), grad(v)) * dx + + assert compute_form_adjoint(a) == conj(inner(grad(v), grad(u))) * dx + + +def test_complex_algebra(self): + z1 = ComplexValue(1j) + z2 = ComplexValue(1+1j) + + # Remember that ufl.algebra functions return ComplexValues, but ufl.mathfunctions return complex Python scalar + # Any operations with a ComplexValue and a complex Python scalar promote to ComplexValue + assert z1*z2 == ComplexValue(-1+1j) + assert z2/z1 == ComplexValue(1-1j) + assert pow(z2, z1) == ComplexValue((1+1j)**1j) + assert sqrt(z2) * as_ufl(1) == ComplexValue(cmath.sqrt(1+1j)) + assert ((sin(z2) + cosh(z2) - atan(z2)) * z1) == ComplexValue((cmath.sin(1+1j) + cmath.cosh(1+1j) - cmath.atan(1+1j))*1j) + assert (abs(z2) - ln(z2))/exp(z1) == ComplexValue((abs(1+1j) - cmath.log(1+1j))/cmath.exp(1j)) + + +def test_automatic_simplification(self): + cell = triangle + element = FiniteElement("Lagrange", cell, 1) + + v = TestFunction(element) + u = TrialFunction(element) + + assert inner(u, v) == u * conj(v) + assert dot(u, v) == u * v + assert outer(u, v) == conj(u) * v + + +def test_apply_algebra_lowering_complex(self): + cell = triangle + element = FiniteElement("Lagrange", cell, 1) + + v = TestFunction(element) + u = TrialFunction(element) + + gv = grad(v) + gu = grad(u) + + a = dot(gu, gv) + b = inner(gv, gu) + c = outer(gu, gv) + + lowered_a = apply_algebra_lowering(a) + lowered_b = apply_algebra_lowering(b) + lowered_c = apply_algebra_lowering(c) + lowered_a_index = lowered_a.index() + lowered_b_index = lowered_b.index() + lowered_c_indices = lowered_c.indices() + + assert lowered_a == gu[lowered_a_index] * gv[lowered_a_index] + assert lowered_b == gv[lowered_b_index] * conj(gu[lowered_b_index]) + assert lowered_c == as_tensor(conj(gu[lowered_c_indices[0]]) * gv[lowered_c_indices[1]], (lowered_c_indices[0],) + (lowered_c_indices[1],)) + + +def test_remove_complex_nodes(self): + cell = triangle + element = FiniteElement("Lagrange", cell, 1) + + u = TrialFunction(element) + v = TestFunction(element) + f = Coefficient(element) + + a = conj(v) + b = real(u) + c = imag(f) + d = conj(real(v))*imag(conj(u)) + + assert remove_complex_nodes(a) == v + assert remove_complex_nodes(b) == u + with pytest.raises(ufl.log.UFLException): + remove_complex_nodes(c) + with pytest.raises(ufl.log.UFLException): + remove_complex_nodes(d) + + +def test_comparison_checker(self): + cell = triangle + element = FiniteElement("Lagrange", cell, 1) + + u = TrialFunction(element) + v = TestFunction(element) + + a = conditional(ge(abs(u), imag(v)), u, v) + b = conditional(le(sqrt(abs(u)), imag(v)), as_ufl(1), as_ufl(1j)) + c = conditional(gt(abs(u), pow(imag(v), 0.5)), sin(u), cos(v)) + d = conditional(lt(as_ufl(-1), as_ufl(1)), u, v) + e = max_value(as_ufl(0), real(u)) + f = min_value(sin(u), cos(v)) + g = min_value(sin(pow(u, 3)), cos(abs(v))) + + assert do_comparison_check(a) == conditional(ge(real(abs(u)), real(imag(v))), u, v) + with pytest.raises(ComplexComparisonError): + b = do_comparison_check(b) + with pytest.raises(ComplexComparisonError): + c = do_comparison_check(c) + assert do_comparison_check(d) == conditional(lt(real(as_ufl(-1)), real(as_ufl(1))), u, v) + assert do_comparison_check(e) == max_value(real(as_ufl(0)), real(real(u))) + assert do_comparison_check(f) == min_value(real(sin(u)), real(cos(v))) + assert do_comparison_check(g) == min_value(real(sin(pow(u, 3))), real(cos(abs(v)))) + + +def test_complex_degree_handling(self): + cell = triangle + element = FiniteElement("Lagrange", cell, 3) + + v = TestFunction(element) + + a = conj(v) + b = imag(v) + c = real(v) + + assert estimate_total_polynomial_degree(a) == 3 + assert estimate_total_polynomial_degree(b) == 3 + assert estimate_total_polynomial_degree(c) == 3 diff -Nru ufl-2018.1.0/test/test_form.py ufl-2019.1.0/test/test_form.py --- ufl-2018.1.0/test/test_form.py 2018-06-14 14:20:02.000000000 +0000 +++ ufl-2019.1.0/test/test_form.py 2019-04-17 18:51:30.000000000 +0000 @@ -125,7 +125,7 @@ u = TrialFunction(V) f = Coefficient(V) g = Coefficient(V) - a = g*inner(grad(u), grad(v))*dx + a = g*inner(grad(v), grad(u))*dx M = a(f, f, coefficients={ g: 1 }) assert M == grad(f)**2*dx diff -Nru ufl-2018.1.0/test/test_indexing.py ufl-2019.1.0/test/test_indexing.py --- ufl-2018.1.0/test/test_indexing.py 2018-06-14 14:20:02.000000000 +0000 +++ ufl-2019.1.0/test/test_indexing.py 2019-04-17 18:51:30.000000000 +0000 @@ -13,13 +13,13 @@ @pytest.fixture def x2(): - x = x1() + x = SpatialCoordinate(triangle) return outer(x, x) @pytest.fixture def x3(): - x = x1() + x = SpatialCoordinate(triangle) return outer(outer(x, x), x) diff -Nru ufl-2018.1.0/test/test_literals.py ufl-2019.1.0/test/test_literals.py --- ufl-2018.1.0/test/test_literals.py 2018-06-14 14:20:02.000000000 +0000 +++ ufl-2019.1.0/test/test_literals.py 2019-04-17 18:51:30.000000000 +0000 @@ -8,7 +8,7 @@ from ufl import * from ufl.classes import Indexed -from ufl.constantvalue import Zero, FloatValue, IntValue, as_ufl +from ufl.constantvalue import Zero, FloatValue, IntValue, ComplexValue, as_ufl def test_zero(self): @@ -27,6 +27,7 @@ assert z1 == z1 assert int(z1) == 0 assert float(z1) == 0.0 + assert complex(z1) == 0.0 + 0.0j self.assertNotEqual(z1, 1.0) self.assertFalse(z1) @@ -67,6 +68,25 @@ assert f2 == f6 # Division produces a FloatValue +def test_complex(self): + f1 = as_ufl(1 + 1j) + f2 = as_ufl(1) + f3 = as_ufl(1j) + f4 = ComplexValue(1 + 1j) + f5 = ComplexValue(1.0 + 1.0j) + f6 = as_ufl(1.0) + f7 = as_ufl(1.0j) + + assert f1 == f1 + assert f1 == f4 + assert f1 == f5 # ComplexValue uses floats + assert f1 == f2 + f3 # Type promotion of IntValue to ComplexValue with arithmetic + assert f4 == f2 + f3 + assert f5 == f2 + f3 + assert f4 == f5 + assert f6 + f7 == f2 + f3 + + def test_scalar_sums(self): n = 10 s = [as_ufl(i) for i in range(n)] diff -Nru ufl-2018.1.0/test/test_sobolevspace.py ufl-2019.1.0/test/test_sobolevspace.py --- ufl-2018.1.0/test/test_sobolevspace.py 2018-06-14 14:20:02.000000000 +0000 +++ ufl-2019.1.0/test/test_sobolevspace.py 2019-04-17 18:51:30.000000000 +0000 @@ -119,8 +119,8 @@ def test_contains_h2(): h2_elements = [ - FiniteElement("ARG", triangle, 1), - FiniteElement("MOR", triangle), + FiniteElement("ARG", triangle, 5), + FiniteElement("MOR", triangle, 2), ] for h2_element in h2_elements: assert h2_element in H2 @@ -235,7 +235,7 @@ P2 = FiniteElement("CG", interval, 2) P3 = FiniteElement("CG", interval, 3) RT1 = FiniteElement("RT", triangle, 1) - ARG = FiniteElement("ARG", triangle, 1) + ARG = FiniteElement("ARG", triangle, 5) # Tensor product elements P1DGP2 = TensorProductElement(P1DG_t, P2) diff -Nru ufl-2018.1.0/test/test_tensoralgebra.py ufl-2019.1.0/test/test_tensoralgebra.py --- ufl-2018.1.0/test/test_tensoralgebra.py 2018-06-14 14:20:02.000000000 +0000 +++ ufl-2019.1.0/test/test_tensoralgebra.py 2019-04-17 18:51:30.000000000 +0000 @@ -5,6 +5,7 @@ import pytest from ufl import * +from ufl.algorithms.remove_complex_nodes import remove_complex_nodes @pytest.fixture(scope="module") @@ -65,13 +66,13 @@ def test_pow2_inner(self, A, u): f = FacetNormal(triangle)[0] f2 = f*f - assert f2 == inner(f, f) + assert f2 == remove_complex_nodes(inner(f, f)) u2 = u**2 - assert u2 == inner(u, u) + assert u2 == remove_complex_nodes(inner(u, u)) A2 = A**2 - assert A2 == inner(A, A) + assert A2 == remove_complex_nodes(inner(A, A)) # Only tensor**2 notation is supported: self.assertRaises(UFLException, lambda: A**3) diff -Nru ufl-2018.1.0/.travis.yml ufl-2019.1.0/.travis.yml --- ufl-2018.1.0/.travis.yml 2018-06-14 14:20:02.000000000 +0000 +++ ufl-2019.1.0/.travis.yml 1970-01-01 00:00:00.000000000 +0000 @@ -1,18 +0,0 @@ -notifications: - slack: - secure: r9t4/J1Y85Tv/D/bJBVPGaW153mSU+I5r+EZrq7E6HnYqFRQ0X0AJw0pUX1ap/MXGtSrHFmQRYmIN1TTUMLPAHlzF4YA4rFA1LxJxxlPQ1TuhmkLbu5DEf5qkSsVeOjIe8+ptx8E8rSUSVfzv1g/RpdgO6mOD6qge94gQhQ4XekrWkZE1OudkqPvHfznCvxtFs3T5dxZWX3suQsGA+l1Gn3gp2YFITTGEfCXFek+rRWL3aqn6Oq+nvI6jfrooei7NSmNQOpGawAW868uvOlWQ2NiH1nsUkOH9U9c4YpTPD7RLLc2r3pAisSzdhZ5LsT+9o+lycqJleGizcVqtF7BP0tBzc1G4uJQwOCCcypLG62VwuwF6HD3zx4bDXD2MTqEneUaLVNlT2JiAsXmg2RKn/z8+5OKrPNe23pDvD3h+NdtHzJxPSnSRynpKD1FBXh1tfiGSVvgY+Mm7Bdu1unuK2tbpb6OLhcq50kYzqW1sbH8GEceb2v10YWV1tWDQ/JORFUtSuM3DFLlTO5XFjTCHqcGMkPQnpmXo+KF4D4uMHmkqVQu1aHdWn99uxoJzyoQrYwa7+94Z/UJuhwSCc+HOr28vXfd5O9ouo3E4FlEMMaAeRiCIGh88Ota6x+7BHmihijroKb9NdB2J2E0HpbLtEQha1yWIsuF0UbWywCojo4= -language: python -python: - - "3.5" - - "3.6" - -before_install: - - pip install flake8 - - pip install pytest - -install: - - pip install . - -script: - - flake8 ufl/ - - py.test test/ diff -Nru ufl-2018.1.0/ufl/algebra.py ufl-2019.1.0/ufl/algebra.py --- ufl-2018.1.0/ufl/algebra.py 2018-06-14 14:20:02.000000000 +0000 +++ ufl-2019.1.0/ufl/algebra.py 2019-04-17 18:51:30.000000000 +0000 @@ -23,9 +23,9 @@ from ufl.log import error from ufl.utils.str import as_native_strings from ufl.core.ufl_type import ufl_type -from ufl.core.expr import Expr, ufl_err_str +from ufl.core.expr import ufl_err_str from ufl.core.operator import Operator -from ufl.constantvalue import Zero, zero, ScalarValue, IntValue, as_ufl +from ufl.constantvalue import Zero, zero, ScalarValue, IntValue, ComplexValue, as_ufl from ufl.checks import is_ufl_scalar, is_true_ufl_scalar from ufl.index_combination_utils import merge_unique_indices from ufl.sorting import sorted_expr @@ -110,7 +110,7 @@ n = len(s) for o in ops[1:]: m = len(o) - if n+m > limit: + if n + m > limit: s += delimop n = m else: @@ -242,7 +242,10 @@ # Simplification "literal a / literal b" -> "literal value of # a/b". Avoiding integer division by casting to float if isinstance(a, ScalarValue) and isinstance(b, ScalarValue): - return as_ufl(float(a._value) / float(b._value)) + try: + return as_ufl(float(a._value) / float(b._value)) + except TypeError: + return as_ufl(complex(a._value) / complex(b._value)) # Simplification "a / a" -> "1" # if not a.ufl_free_indices and not a.ufl_shape and a == b: # return as_ufl(1) @@ -265,7 +268,11 @@ a = a.evaluate(x, mapping, component, index_values) b = b.evaluate(x, mapping, component, index_values) # Avoiding integer division by casting to float - return float(a) / float(b) + try: + e = float(a) / float(b) + except TypeError: + e = complex(a) / complex(b) + return e def __str__(self): return "%s / %s" % (parstr(self.ufl_operands[0], self), @@ -292,7 +299,11 @@ # Simplification if isinstance(a, ScalarValue) and isinstance(b, ScalarValue): return as_ufl(a._value ** b._value) + if isinstance(b, Zero): + return IntValue(1) if isinstance(a, Zero) and isinstance(b, ScalarValue): + if isinstance(b, ComplexValue): + error("Cannot raise zero to a complex power.") bf = float(b) if bf < 0: error("Division by zero, cannot raise 0 to a negative power.") @@ -300,8 +311,6 @@ return zero() if isinstance(b, ScalarValue) and b._value == 1: return a - if isinstance(b, Zero): - return IntValue(1) # Construction self = Operator.__new__(cls) @@ -333,10 +342,21 @@ class Abs(Operator): __slots__ = () + def __new__(cls, a): + a = as_ufl(a) + + # Simplification + if isinstance(a, (Zero, Abs)): + return a + if isinstance(a, Conj): + return Abs(a.ufl_operands[0]) + if isinstance(a, ScalarValue): + return as_ufl(abs(a._value)) + + return Operator.__new__(cls) + def __init__(self, a): Operator.__init__(self, (a,)) - if not isinstance(a, Expr): - error("Expecting Expr instance, not %s." % ufl_err_str(a)) def evaluate(self, x, mapping, component, index_values): a = self.ufl_operands[0].evaluate(x, mapping, component, index_values) @@ -345,3 +365,95 @@ def __str__(self): a, = self.ufl_operands return "|%s|" % (parstr(a, self),) + + +@ufl_type(num_ops=1, + inherit_shape_from_operand=0, inherit_indices_from_operand=0) +class Conj(Operator): + __slots__ = () + + def __new__(cls, a): + a = as_ufl(a) + + # Simplification + if isinstance(a, (Abs, Real, Imag, Zero)): + return a + if isinstance(a, Conj): + return a.ufl_operands[0] + if isinstance(a, ScalarValue): + return as_ufl(a._value.conjugate()) + + return Operator.__new__(cls) + + def __init__(self, a): + Operator.__init__(self, (a,)) + + def evaluate(self, x, mapping, component, index_values): + a = self.ufl_operands[0].evaluate(x, mapping, component, index_values) + return a.conjugate() + + def __str__(self): + a, = self.ufl_operands + return "conj(%s)" % (parstr(a, self),) + + +@ufl_type(num_ops=1, + inherit_shape_from_operand=0, inherit_indices_from_operand=0) +class Real(Operator): + __slots__ = () + + def __new__(cls, a): + a = as_ufl(a) + + # Simplification + if isinstance(a, Conj): + a = a.ufl_operands[0] + if isinstance(a, Zero): + return a + if isinstance(a, ScalarValue): + return as_ufl(a.real()) + if isinstance(a, Real): + a = a.ufl_operands[0] + + return Operator.__new__(cls) + + def __init__(self, a): + Operator.__init__(self, (a,)) + + def evaluate(self, x, mapping, component, index_values): + a = self.ufl_operands[0].evaluate(x, mapping, component, index_values) + return a.real + + def __str__(self): + a, = self.ufl_operands + return "Re[%s]" % (parstr(a, self),) + + +@ufl_type(num_ops=1, + inherit_shape_from_operand=0, inherit_indices_from_operand=0) +class Imag(Operator): + __slots__ = () + + def __new__(cls, a): + a = as_ufl(a) + + # Simplification + if isinstance(a, Zero): + return a + if isinstance(a, (Real, Imag, Abs)): + return Zero(a.ufl_shape, a.ufl_free_indices, a.ufl_index_dimensions) + if isinstance(a, ScalarValue): + return as_ufl(a.imag()) + + return Operator.__new__(cls) + + def __init__(self, a): + Operator.__init__(self, (a,)) + + def evaluate(self, x, mapping, component, index_values): + a = self.ufl_operands[0].evaluate(x, mapping, component, index_values) + return a.imag + + def __str__(self): + a, = self.ufl_operands + return "Im[%s]" % (parstr(a, self),) diff -Nru ufl-2018.1.0/ufl/algorithms/apply_algebra_lowering.py ufl-2019.1.0/ufl/algorithms/apply_algebra_lowering.py --- ufl-2018.1.0/ufl/algorithms/apply_algebra_lowering.py 2018-06-14 14:20:02.000000000 +0000 +++ ufl-2019.1.0/ufl/algorithms/apply_algebra_lowering.py 2019-04-17 18:51:30.000000000 +0000 @@ -23,7 +23,7 @@ from ufl.log import error -from ufl.classes import Product, Grad +from ufl.classes import Product, Grad, Conj from ufl.core.multiindex import indices, Index, FixedIndex from ufl.tensors import as_tensor, as_matrix, as_vector @@ -36,6 +36,7 @@ class LowerCompoundAlgebra(MultiFunction): """Expands high level compound operators (e.g. inner) to equivalent representations using basic operators (e.g. index notation).""" + def __init__(self): MultiFunction.__init__(self) @@ -87,16 +88,16 @@ else: k = (Index(),) # Potentially creates a single IndexSum over a Product - s = a[ai+k]*b[k+bi] - return as_tensor(s, ai+bi) + s = a[ai + k] * b[k + bi] + return as_tensor(s, ai + bi) def dot(self, o, a, b): - ai = indices(len(a.ufl_shape)-1) - bi = indices(len(b.ufl_shape)-1) + ai = indices(len(a.ufl_shape) - 1) + bi = indices(len(b.ufl_shape) - 1) k = (Index(),) # Creates a single IndexSum over a Product - s = a[ai+k]*b[k+bi] - return as_tensor(s, ai+bi) + s = a[ai + k] * b[k + bi] + return as_tensor(s, ai + bi) def alternative_inner(self, o, a, b): # TODO: Test this ash = a.ufl_shape @@ -115,7 +116,7 @@ ii = tuple(ii) # ii = indices(len(a.ufl_shape)) # Potentially creates multiple IndexSums over a Product - s = a[ii]*b[ii] + s = a[ii] * Conj(b[ii]) return s def inner(self, o, a, b): @@ -125,15 +126,15 @@ error("Nonmatching shapes.") ii = indices(len(ash)) # Creates multiple IndexSums over a Product - s = a[ii]*b[ii] + s = a[ii] * Conj(b[ii]) return s def outer(self, o, a, b): ii = indices(len(a.ufl_shape)) jj = indices(len(b.ufl_shape)) # Create a Product with no shared indices - s = a[ii]*b[jj] - return as_tensor(s, ii+jj) + s = Conj(a[ii]) * b[jj] + return as_tensor(s, ii + jj) def determinant(self, o, A): return determinant_expr(A) diff -Nru ufl-2018.1.0/ufl/algorithms/apply_derivatives.py ufl-2019.1.0/ufl/algorithms/apply_derivatives.py --- ufl-2018.1.0/ufl/algorithms/apply_derivatives.py 2018-06-14 14:20:02.000000000 +0000 +++ ufl-2019.1.0/ufl/algorithms/apply_derivatives.py 2019-04-17 18:51:30.000000000 +0000 @@ -33,6 +33,7 @@ from ufl.classes import Indexed, ListTensor, ComponentTensor from ufl.classes import ExprList, ExprMapping from ufl.classes import Product, Sum, IndexSum +from ufl.classes import Conj, Real, Imag from ufl.classes import JacobianInverse from ufl.classes import SpatialCoordinate @@ -49,7 +50,7 @@ from ufl.algorithms.map_integrands import map_integrand_dags from ufl.checks import is_cellwise_constant - +from ufl.differentiation import CoordinateDerivative # TODO: Add more rulesets? # - DivRuleset # - CurlRuleset @@ -276,11 +277,11 @@ if isinstance(gp, Zero): # This probably produces better results for the common # case of f**constant - op = fp * g * f**(g-1) + op = fp * g * f**(g - 1) else: # Note: This produces expressions like (1/w)*w**5 instead of w**4 # op = o * (fp * g / f + gp * ln(f)) # This reuses o - op = f**(g-1) * (g*fp + f*ln(f)*gp) # This gives better accuracy in dolfin integration test + op = f**(g - 1) * (g * fp + f * ln(f) * gp) # This gives better accuracy in dolfin integration test # Example: d/dx[x**(x**3)]: # f = x @@ -299,6 +300,17 @@ # return conditional(eq(f, 0), 0, Product(sign(f), df)) return sign(f) * df + # --- Complex algebra + + def conj(self, o, df): + return Conj(df) + + def real(self, o, df): + return Real(df) + + def imag(self, o, df): + return Imag(df) + # --- Mathfunctions def math_function(self, o, df): @@ -310,7 +322,7 @@ error("Unknown math function.") def sqrt(self, o, fp): - return fp / (2*o) + return fp / (2 * o) def exp(self, o, fp): return fp * o @@ -331,7 +343,7 @@ def tan(self, o, fp): f, = o.ufl_operands - return 2.0*fp / (cos(2.0*f) + 1.0) + return 2.0 * fp / (cos(2.0 * f) + 1.0) def cosh(self, o, fp): f, = o.ufl_operands @@ -345,7 +357,7 @@ f, = o.ufl_operands def sech(y): - return (2.0*cosh(y)) / (cosh(2.0*y) + 1.0) + return (2.0 * cosh(y)) / (cosh(2.0 * y) + 1.0) return fp * sech(f)**2 def acos(self, o, fp): @@ -362,7 +374,7 @@ def atan_2(self, o, fp, gp): f, g = o.ufl_operands - return (g*fp - f*gp) / (f**2 + g**2) + return (g * fp - f * gp) / (f**2 + g**2) def erf(self, o, fp): f, = o.ufl_operands @@ -378,7 +390,7 @@ if isinstance(nu, Zero): op = -bessel_J(1, f) else: - op = 0.5 * (bessel_J(nu-1, f) - bessel_J(nu+1, f)) + op = 0.5 * (bessel_J(nu - 1, f) - bessel_J(nu + 1, f)) return op * fp def bessel_y(self, o, nup, fp): @@ -389,7 +401,7 @@ if isinstance(nu, Zero): op = -bessel_Y(1, f) else: - op = 0.5 * (bessel_Y(nu-1, f) - bessel_Y(nu+1, f)) + op = 0.5 * (bessel_Y(nu - 1, f) - bessel_Y(nu + 1, f)) return op * fp def bessel_i(self, o, nup, fp): @@ -400,7 +412,7 @@ if isinstance(nu, Zero): op = bessel_I(1, f) else: - op = 0.5 * (bessel_I(nu-1, f) + bessel_I(nu+1, f)) + op = 0.5 * (bessel_I(nu - 1, f) + bessel_I(nu + 1, f)) return op * fp def bessel_k(self, o, nup, fp): @@ -411,7 +423,7 @@ if isinstance(nu, Zero): op = -bessel_K(1, f) else: - op = -0.5 * (bessel_K(nu-1, f) + bessel_K(nu+1, f)) + op = -0.5 * (bessel_K(nu - 1, f) + bessel_K(nu + 1, f)) return op * fp # --- Restrictions @@ -445,7 +457,7 @@ # or df become NaN or Inf in floating point computations! c = o.ufl_operands[0] dc = conditional(c, 1, 0) - return dc*dt + (1.0 - dc)*df + return dc * dt + (1.0 - dc) * df else: # Not placing t[1],f[1] outside, allowing arguments inside # conditionals. This will make legacy ffc fail, but @@ -461,7 +473,7 @@ # conditionals f, g = o.ufl_operands dc = conditional(f > g, 1, 0) - return dc*df + (1.0 - dc)*dg + return dc * df + (1.0 - dc) * dg def min_value(self, o, df, dg): # d/dx min(f, g) = @@ -471,7 +483,7 @@ # inside conditionals f, g = o.ufl_operands dc = conditional(f < g, 1, 0) - return dc*df + (1.0 - dc)*dg + return dc * df + (1.0 - dc) * dg class GradRuleset(GenericDerivativeRuleset): @@ -499,7 +511,7 @@ error("ReferenceValue can only wrap a terminal") r = indices(len(o.ufl_shape)) i, j = indices(2) - Do = as_tensor(o[j, i]*ReferenceGrad(o)[r + (j,)], r + (i,)) + Do = as_tensor(o[j, i] * ReferenceGrad(o)[r + (j,)], r + (i,)) return Do # TODO: Add more explicit geometry type handlers here, with @@ -541,7 +553,7 @@ K = JacobianInverse(domain) r = indices(len(o.ufl_shape)) i, j = indices(2) - Do = as_tensor(K[j, i]*ReferenceGrad(o)[r + (j,)], r + (i,)) + Do = as_tensor(K[j, i] * ReferenceGrad(o)[r + (j,)], r + (i,)) return Do def reference_grad(self, o): @@ -554,7 +566,7 @@ K = JacobianInverse(domain) r = indices(len(o.ufl_shape)) i, j = indices(2) - Do = as_tensor(K[j, i]*ReferenceGrad(o)[r + (j,)], r + (i,)) + Do = as_tensor(K[j, i] * ReferenceGrad(o)[r + (j,)], r + (i,)) return Do # --- Nesting of gradients @@ -754,6 +766,7 @@ D_w[v](e) = d/dtau e(w+tau v)|tau=0 """ + def __init__(self, coefficients, arguments, coefficient_derivatives): GenericDerivativeRuleset.__init__(self, var_shape=()) @@ -774,7 +787,7 @@ # Build more convenient dict {f: df/dw} for each coefficient f # where df/dw is nonzero cd = coefficient_derivatives.ufl_operands - self._cd = {cd[2*i]: cd[2*i+1] for i in range(len(cd)//2)} + self._cd = {cd[2 * i]: cd[2 * i + 1] for i in range(len(cd) // 2)} # Explicitly defining dg/dw == 0 geometric_quantity = GenericDerivativeRuleset.independent_terminal @@ -829,7 +842,7 @@ rv = len(v.ufl_shape) oi1 = oi[:-rv] oi2 = oi[-rv:] - prod = so*v[oi2] + prod = so * v[oi2] if oi1: dosum += as_tensor(prod, oi1) else: @@ -912,14 +925,14 @@ # Apply gradients directly to argument vval, and get the # right indexed scalar component(s) kk = indices(ngrads) - Dvkk = apply_grads(vval)[vcomp+kk] + Dvkk = apply_grads(vval)[vcomp + kk] # Place scalar component(s) Dvkk into the right tensor # positions if wshape: Ejj, jj = unit_indexed_tensor(wshape, wcomp) else: Ejj, jj = 1, () - gprimeterm = as_tensor(Ejj*Dvkk, jj+kk) + gprimeterm = as_tensor(Ejj * Dvkk, jj + kk) return gprimeterm # Accumulate contributions from variations in different @@ -1001,7 +1014,7 @@ rv = len(v.ufl_shape) oi1 = oi[:-rv] oi2 = oi[-rv:] - prod = so*v[oi2] + prod = so * v[oi2] if oi1: gprimesum += as_tensor(prod, oi1) else: @@ -1009,6 +1022,10 @@ return gprimesum + def coordinate_derivative(self, o): + o = o.ufl_operands + return CoordinateDerivative(map_expr_dag(self, o[0]), o[1], o[2], o[3]) + class DerivativeRuleDispatcher(MultiFunction): def __init__(self): @@ -1039,6 +1056,10 @@ rules = GateauxDerivativeRuleset(w, v, cd) return map_expr_dag(rules, f) + def coordinate_derivative(self, o, f, dummy_w, dummy_v, dummy_cd): + o_ = o.ufl_operands + return CoordinateDerivative(map_expr_dag(self, o_[0]), o_[1], o_[2], o_[3]) + def indexed(self, o, Ap, ii): # TODO: (Partially) duplicated in generic rules # Reuse if untouched if Ap is o.ufl_operands[0]: @@ -1072,3 +1093,133 @@ def apply_derivatives(expression): rules = DerivativeRuleDispatcher() return map_integrand_dags(rules, expression) + + +class CoordinateDerivativeRuleset(GenericDerivativeRuleset): + """Apply AFD (Automatic Functional Differentiation) to expression. + + Implements rules for the Gateaux derivative D_w[v](...) defined as + + D_w[v](e) = d/dtau e(w+tau v)|tau=0 + + where 'e' is a ufl form after pullback and w is a SpatialCoordinate. + + """ + + def __init__(self, coefficients, arguments, coefficient_derivatives): + GenericDerivativeRuleset.__init__(self, var_shape=()) + + # Type checking + if not isinstance(coefficients, ExprList): + error("Expecting a ExprList of coefficients.") + if not isinstance(arguments, ExprList): + error("Expecting a ExprList of arguments.") + if not isinstance(coefficient_derivatives, ExprMapping): + error("Expecting a coefficient-coefficient ExprMapping.") + + # The coefficient(s) to differentiate w.r.t. and the + # argument(s) s.t. D_w[v](e) = d/dtau e(w+tau v)|tau=0 + self._w = coefficients.ufl_operands + self._v = arguments.ufl_operands + self._w2v = {w: v for w, v in zip(self._w, self._v)} + + # Build more convenient dict {f: df/dw} for each coefficient f + # where df/dw is nonzero + cd = coefficient_derivatives.ufl_operands + self._cd = {cd[2 * i]: cd[2 * i + 1] for i in range(len(cd) // 2)} + + # Explicitly defining dg/dw == 0 + geometric_quantity = GenericDerivativeRuleset.independent_terminal + + # Explicitly defining da/dw == 0 + argument = GenericDerivativeRuleset.independent_terminal + + def coefficient(self, o): + error("CoordinateDerivative of coefficient in physical space is not implemented.") + + def grad(self, o): + error("CoordinateDerivative grad in physical space is not implemented.") + + def spatial_coordinate(self, o): + do = self._w2v.get(o) + # d x /d x => Argument(x.function_space()) + if do is not None: + return do + else: + error("Not implemented: CoordinateDerivative found a SpatialCoordinate that is different from the one being differentiated.") + + def reference_value(self, o): + do = self._cd.get(o) + if do is not None: + return do + else: + return self.independent_terminal(o) + + def reference_grad(self, g): + # d (grad_X(...(x)) / dx => grad_X(...(Argument(x.function_space())) + o = g + ngrads = 0 + while isinstance(o, ReferenceGrad): + o, = o.ufl_operands + ngrads += 1 + if not (isinstance(o, SpatialCoordinate) or isinstance(o.ufl_operands[0], FormArgument)): + error("Expecting gradient of a FormArgument, not %s" % ufl_err_str(o)) + + def apply_grads(f): + for i in range(ngrads): + f = ReferenceGrad(f) + return f + + # Find o among all w without any indexing, which makes this + # easy + for (w, v) in zip(self._w, self._v): + if o == w and isinstance(v, ReferenceValue) and isinstance(v.ufl_operands[0], FormArgument): + # Case: d/dt [w + t v] + return apply_grads(v) + return self.independent_terminal(o) + + def jacobian(self, o): + # d (grad_X(x))/d x => grad_X(Argument(x.function_space()) + for (w, v) in zip(self._w, self._v): + if o.ufl_domain() == w.ufl_domain() and isinstance(v.ufl_operands[0], FormArgument): + return ReferenceGrad(v) + return self.independent_terminal(o) + + +class CoordinateDerivativeRuleDispatcher(MultiFunction): + def __init__(self): + MultiFunction.__init__(self) + + def terminal(self, o): + return o + + def derivative(self, o): + error("Missing derivative handler for {0}.".format(type(o).__name__)) + + expr = MultiFunction.reuse_if_untouched + + def grad(self, o): + return o + + def reference_grad(self, o): + return o + + def coefficient_derivative(self, o): + return o + + def coordinate_derivative(self, o): + from ufl.algorithms import extract_unique_elements + spaces = set(c.family() for c in extract_unique_elements(o)) + unsupported_spaces = {"Argyris", "Bell", "Hermite", "Morley"} + if spaces & unsupported_spaces: + error("CoordinateDerivative is not supported for elements of type %s. " + "This is because their pullback is not implemented in UFL." % unsupported_spaces) + f, w, v, cd = o.ufl_operands + f = self(f) # transform f + rules = CoordinateDerivativeRuleset(w, v, cd) + return map_expr_dag(rules, f) + + +def apply_coordinate_derivatives(expression): + rules = CoordinateDerivativeRuleDispatcher() + return map_integrand_dags(rules, expression) diff -Nru ufl-2018.1.0/ufl/algorithms/apply_function_pullbacks.py ufl-2019.1.0/ufl/algorithms/apply_function_pullbacks.py --- ufl-2018.1.0/ufl/algorithms/apply_function_pullbacks.py 2018-06-14 14:20:02.000000000 +0000 +++ ufl-2019.1.0/ufl/algorithms/apply_function_pullbacks.py 2019-04-17 18:51:30.000000000 +0000 @@ -52,7 +52,7 @@ if len(shape) == 0: return [None] elif len(shape) == 1: - return [None]*shape[0] + return [None] * shape[0] else: return [create_nested_lists(shape[1:]) for i in range(shape[0])] @@ -66,7 +66,7 @@ return components else: n = product(shape[1:]) - return [reshape_to_nested_list(components[n*i:n*(i+1)], shape[1:]) for i in range(shape[0])] + return [reshape_to_nested_list(components[n * i:n * (i + 1)], shape[1:]) for i in range(shape[0])] def apply_single_function_pullbacks(g): @@ -96,7 +96,7 @@ # Create contravariant transform for reuse (note that detJ is the # _signed_ (pseudo-)determinant) - transform_hdiv = (1.0/detJ) * J + transform_hdiv = (1.0 / detJ) * J # Shortcut simple cases for a more efficient representation, # including directly Piola-mapped elements and mixed elements of @@ -114,25 +114,25 @@ elif mapping == "contravariant Piola": assert transform_hdiv.ufl_shape == (gsize, rsize) i, j = indices(2) - f = as_vector(transform_hdiv[i, j]*r[j], i) + f = as_vector(transform_hdiv[i, j] * r[j], i) # f = as_tensor(transform_hdiv[i, j]*r[k,j], (k,i)) # FIXME: Handle Vector(Piola) here? assert f.ufl_shape == g.ufl_shape return f elif mapping == "covariant Piola": assert Jinv.ufl_shape == (rsize, gsize) i, j = indices(2) - f = as_vector(Jinv[j, i]*r[j], i) + f = as_vector(Jinv[j, i] * r[j], i) # f = as_tensor(Jinv[j, i]*r[k,j], (k,i)) # FIXME: Handle Vector(Piola) here? assert f.ufl_shape == g.ufl_shape return f elif mapping == "double covariant Piola": i, j, m, n = indices(4) - f = as_tensor(Jinv[m, i]*r[m, n]*Jinv[n, j], (i, j)) + f = as_tensor(Jinv[m, i] * r[m, n] * Jinv[n, j], (i, j)) assert f.ufl_shape == g.ufl_shape return f elif mapping == "double contravariant Piola": i, j, m, n = indices(4) - f = as_tensor((1.0/detJ)*(1.0/detJ)*J[i, m]*r[m, n]*J[j, n], (i, j)) + f = as_tensor((1.0 / detJ) * (1.0 / detJ) * J[i, m] * r[m, n] * J[j, n], (i, j)) assert f.ufl_shape == g.ufl_shape return f @@ -151,7 +151,7 @@ assert len(gsh) == 1 assert len(rsh) == 1 - g_components = [None]*gsize + g_components = [None] * gsize gpos = 0 rpos = 0 for subelm in sub_elements_with_mappings(element): @@ -182,46 +182,49 @@ elif mp == "contravariant Piola": assert transform_hdiv.ufl_shape == (gm, rm) # Get reference value vector corresponding to this subelement: - rv = as_vector([r[rpos+k] for k in range(rm)]) + rv = as_vector([r[rpos + k] for k in range(rm)]) # Apply transform with IndexSum over j for each row j = Index() for i in range(gm): - g_components[gpos + i] = transform_hdiv[i, j]*rv[j] + g_components[gpos + i] = transform_hdiv[i, j] * rv[j] elif mp == "covariant Piola": assert Jinv.ufl_shape == (rm, gm) # Get reference value vector corresponding to this subelement: - rv = as_vector([r[rpos+k] for k in range(rm)]) + rv = as_vector([r[rpos + k] for k in range(rm)]) # Apply transform with IndexSum over j for each row j = Index() for i in range(gm): - g_components[gpos + i] = Jinv[j, i]*rv[j] + g_components[gpos + i] = Jinv[j, i] * rv[j] elif mp == "double covariant Piola": # components are flatten, map accordingly - rv = as_vector([r[rpos+k] for k in range(rm)]) - dim = subelm.value_shape()[0] - for i in range(dim): - for j in range(dim): + rv = as_vector([r[rpos + k] for k in range(rm)]) + (gdim, _) = subelm.value_shape() + (rdim, _) = subelm.reference_value_shape() + for i in range(gdim): + for j in range(gdim): gv = 0 # int times Index is not allowed. so sum by hand - for m in range(dim): - for n in range(dim): - gv += Jinv[m, i]*rv[m*dim+n]*Jinv[n, j] - g_components[gpos + i * dim + j] = gv + for m in range(rdim): + for n in range(rdim): + gv += Jinv[m, i] * rv[m * rdim + n] * Jinv[n, j] + g_components[gpos + i * gdim + j] = gv elif mp == "double contravariant Piola": # components are flatten, map accordingly - rv = as_vector([r[rpos+k] for k in range(rm)]) - dim = subelm.value_shape()[0] - for i in range(dim): - for j in range(dim): + rv = as_vector([r[rpos + k] for k in range(rm)]) + (gdim, _) = subelm.value_shape() + (rdim, _) = subelm.reference_value_shape() + for i in range(gdim): + for j in range(gdim): gv = 0 # int times Index is not allowed. so sum by hand - for m in range(dim): - for n in range(dim): - gv += (1.0/detJ)*(1.0/detJ)*J[i, m]*rv[m*dim+n]*J[j, n] - g_components[gpos + i * dim + j] = gv + for m in range(rdim): + for n in range(rdim): + gv += ((1.0 / detJ) * (1.0 / detJ) * + J[i, m] * rv[m * rdim + n] * J[j, n]) + g_components[gpos + i * gdim + j] = gv else: error("Unknown subelement mapping type %s for element %s." % (mp, str(subelm))) diff -Nru ufl-2018.1.0/ufl/algorithms/apply_geometry_lowering.py ufl-2019.1.0/ufl/algorithms/apply_geometry_lowering.py --- ufl-2018.1.0/ufl/algorithms/apply_geometry_lowering.py 2018-06-14 14:20:02.000000000 +0000 +++ ufl-2019.1.0/ufl/algorithms/apply_geometry_lowering.py 2019-04-17 18:51:30.000000000 +0000 @@ -57,7 +57,7 @@ def __init__(self, preserve_types=()): MultiFunction.__init__(self) # Store preserve_types as boolean lookup table - self._preserve_types = [False]*Expr._ufl_num_typecodes_ + self._preserve_types = [False] * Expr._ufl_num_typecodes_ for cls in preserve_types: self._preserve_types[cls._ufl_typecode_] = True @@ -114,7 +114,7 @@ # specific type for the unsigned pseudo-determinant. if domain.topological_dimension() < domain.geometric_dimension(): co = CellOrientation(domain) - detJ = co*detJ + detJ = co * detJ return detJ @@ -127,7 +127,7 @@ J = self.jacobian(Jacobian(domain)) RFJ = CellFacetJacobian(domain) i, j, k = indices(3) - return as_tensor(J[i, k]*RFJ[k, j], (i, j)) + return as_tensor(J[i, k] * RFJ[k, j], (i, j)) @memoized_handler def facet_jacobian_inverse(self, o): @@ -249,7 +249,7 @@ edges = CellEdgeVectors(domain) num_edges = edges.ufl_shape[0] j = Index() - elen = [sqrt(edges[e, j]*edges[e, j]) for e in range(num_edges)] + elen = [sqrt(edges[e, j] * edges[e, j]) for e in range(num_edges)] if cellname == "triangle": return (elen[0] * elen[1] * elen[2]) / (4.0 * cellvolume) @@ -296,7 +296,7 @@ edges = CellEdgeVectors(domain) num_edges = edges.ufl_shape[0] j = Index() - elen2 = [edges[e, j]*edges[e, j] for e in range(num_edges)] + elen2 = [edges[e, j] * edges[e, j] for e in range(num_edges)] return sqrt(reduce(reduction_op, elen2)) @memoized_handler @@ -320,7 +320,7 @@ verts = CellVertices(domain) verts = [verts[v, ...] for v in range(verts.ufl_shape[0])] j = Index() - elen2 = ((v0-v1)[j]*(v0-v1)[j] for v0, v1 in combinations(verts, 2)) + elen2 = ((v0 - v1)[j] * (v0 - v1)[j] for v0, v1 in combinations(verts, 2)) return sqrt(reduce(max_value, elen2)) @memoized_handler @@ -350,7 +350,7 @@ edges = FacetEdgeVectors(domain) num_edges = edges.ufl_shape[0] j = Index() - elen2 = [edges[e, j]*edges[e, j] for e in range(num_edges)] + elen2 = [edges[e, j] * edges[e, j] for e in range(num_edges)] return sqrt(reduce(reduction_op, elen2)) @memoized_handler @@ -381,7 +381,7 @@ # Return normalized vector, sign corrected by cell # orientation co = CellOrientation(domain) - return co * cell_normal / sqrt(cell_normal[i]*cell_normal[i]) + return co * cell_normal / sqrt(cell_normal[i] * cell_normal[i]) else: error("What do you want cell normal in gdim={0}, tdim={1} to be?".format(gdim, tdim)) @@ -404,7 +404,7 @@ nlen = abs(ndir[0]) else: i = Index() - nlen = sqrt(ndir[i]*ndir[i]) + nlen = sqrt(ndir[i] * ndir[i]) rn = ReferenceNormal(domain) # +/- 1.0 here n = rn[0] * ndir / nlen @@ -423,7 +423,7 @@ # normalise i = Index() - n = ndir / sqrt(ndir[i]*ndir[i]) + n = ndir / sqrt(ndir[i] * ndir[i]) r = n if r.ufl_shape != o.ufl_shape: diff -Nru ufl-2018.1.0/ufl/algorithms/apply_integral_scaling.py ufl-2019.1.0/ufl/algorithms/apply_integral_scaling.py --- ufl-2018.1.0/ufl/algorithms/apply_integral_scaling.py 2018-06-14 14:20:02.000000000 +0000 +++ ufl-2019.1.0/ufl/algorithms/apply_integral_scaling.py 2019-04-17 18:51:30.000000000 +0000 @@ -21,6 +21,9 @@ from ufl.log import error from ufl.classes import JacobianDeterminant, FacetJacobianDeterminant, QuadratureWeight, Form, Integral from ufl.measure import custom_integral_types, point_integral_types +from ufl.differentiation import CoordinateDerivative +from ufl.algorithms.apply_geometry_lowering import apply_geometry_lowering +from ufl.algorithms.estimate_degrees import estimate_total_polynomial_degree def compute_integrand_scaling_factor(integral): @@ -33,14 +36,22 @@ tdim = domain.topological_dimension() # gdim = domain.geometric_dimension() + # Polynomial degree of integrand scaling + degree = 0 if integral_type == "cell": - scale = abs(JacobianDeterminant(domain)) * weight + detJ = JacobianDeterminant(domain) + degree = estimate_total_polynomial_degree(apply_geometry_lowering(detJ)) + # Despite the abs, |detJ| is polynomial except for + # self-intersecting cells, where we have other problems. + scale = abs(detJ) * weight elif integral_type.startswith("exterior_facet"): if tdim > 1: # Scaling integral by facet jacobian determinant and # quadrature weight - scale = FacetJacobianDeterminant(domain) * weight + detFJ = FacetJacobianDeterminant(domain) + degree = estimate_total_polynomial_degree(apply_geometry_lowering(detFJ)) + scale = detFJ * weight else: # No need to scale 'integral' over a vertex scale = 1 @@ -49,7 +60,9 @@ if tdim > 1: # Scaling integral by facet jacobian determinant from one # side and quadrature weight - scale = FacetJacobianDeterminant(domain)('+') * weight + detFJ = FacetJacobianDeterminant(domain) + degree = estimate_total_polynomial_degree(apply_geometry_lowering(detFJ)) + scale = detFJ('+') * weight else: # No need to scale 'integral' over a vertex scale = 1 @@ -66,7 +79,7 @@ else: error("Unknown integral type {}, don't know how to scale.".format(integral_type)) - return scale + return scale, degree def apply_integral_scaling(form): @@ -80,10 +93,34 @@ elif isinstance(form, Integral): integral = form - # Compute and apply integration scaling factor - scale = compute_integrand_scaling_factor(integral) - newintegrand = integral.integrand() * scale - return integral.reconstruct(integrand=newintegrand) + integrand = integral.integrand() + # Compute and apply integration scaling factor since we want to compute + # coordinate derivatives at the end, the scaling factor has to be moved + # inside those + scale, degree = compute_integrand_scaling_factor(integral) + md = {} + md.update(integral.metadata()) + new_degree = degree + cur_degree = md.get("estimated_polynomial_degree") + if cur_degree is not None: + if isinstance(cur_degree, tuple) and isinstance(degree, tuple): + new_degree = tuple(d[0] + d[1] for d in zip(cur_degree, degree)) + elif isinstance(cur_degree, tuple): + new_degree = tuple(d + degree for d in cur_degree) + elif isinstance(degree, tuple): + new_degree = tuple(cur_degree + d for d in degree) + else: + new_degree = cur_degree + degree + md["estimated_polynomial_degree"] = new_degree + + def scale_coordinate_derivative(o, scale): + o_ = o.ufl_operands + if isinstance(o, CoordinateDerivative): + return CoordinateDerivative(scale_coordinate_derivative(o_[0], scale), o_[1], o_[2], o_[3]) + else: + return scale * o + newintegrand = scale_coordinate_derivative(integrand, scale) + return integral.reconstruct(integrand=newintegrand, metadata=md) else: error("Invalid type %s" % (form.__class__.__name__,)) diff -Nru ufl-2018.1.0/ufl/algorithms/balancing.py ufl-2019.1.0/ufl/algorithms/balancing.py --- ufl-2018.1.0/ufl/algorithms/balancing.py 1970-01-01 00:00:00.000000000 +0000 +++ ufl-2019.1.0/ufl/algorithms/balancing.py 2019-04-17 18:51:30.000000000 +0000 @@ -0,0 +1,76 @@ +# -*- coding: utf-8 -*- +# Copyright (C) 2011-2017 Martin Sandve Alnæs +# +# This file is part of UFL (https://www.fenicsproject.org) +# +# SPDX-License-Identifier: LGPL-3.0-or-later + +from ufl.classes import (CellAvg, FacetAvg, Grad, Indexed, NegativeRestricted, + PositiveRestricted, ReferenceGrad, ReferenceValue) +from ufl.corealg.map_dag import map_expr_dag +from ufl.corealg.multifunction import MultiFunction + +modifier_precedence = [ + ReferenceValue, ReferenceGrad, Grad, CellAvg, FacetAvg, PositiveRestricted, + NegativeRestricted, Indexed +] + +modifier_precedence = { + m._ufl_handler_name_: i + for i, m in enumerate(modifier_precedence) +} + + +def balance_modified_terminal(expr): + # NB! Assuming e.g. grad(cell_avg(expr)) does not occur, + # i.e. it is simplified to 0 immediately. + + if expr._ufl_is_terminal_: + return expr + + assert expr._ufl_is_terminal_modifier_ + + orig = expr + + # Build list of modifier layers + layers = [expr] + while not expr._ufl_is_terminal_: + assert expr._ufl_is_terminal_modifier_ + expr = expr.ufl_operands[0] + layers.append(expr) + assert layers[-1] is expr + assert expr._ufl_is_terminal_ + + # Apply modifiers in order + layers = sorted( + layers[:-1], key=lambda e: modifier_precedence[e._ufl_handler_name_]) + for op in layers: + ops = (expr, ) + op.ufl_operands[1:] + expr = op._ufl_expr_reconstruct_(*ops) + + # Preserve id if nothing has changed + return orig if expr == orig else expr + + +class BalanceModifiers(MultiFunction): + def expr(self, expr, *ops): + return expr._ufl_expr_reconstruct_(*ops) + + def terminal(self, expr): + return expr + + def _modifier(self, expr, *ops): + return balance_modified_terminal(expr) + + reference_value = _modifier + reference_grad = _modifier + grad = _modifier + cell_avg = _modifier + facet_avg = _modifier + positive_restricted = _modifier + negative_restricted = _modifier + + +def balance_modifiers(expr): + mf = BalanceModifiers() + return map_expr_dag(mf, expr) diff -Nru ufl-2018.1.0/ufl/algorithms/change_to_reference.py ufl-2019.1.0/ufl/algorithms/change_to_reference.py --- ufl-2018.1.0/ufl/algorithms/change_to_reference.py 2018-06-14 14:20:02.000000000 +0000 +++ ufl-2019.1.0/ufl/algorithms/change_to_reference.py 2019-04-17 18:51:30.000000000 +0000 @@ -226,7 +226,7 @@ K = K(restricted) # Create Hdiv mapping from possibly restricted geometry objects - Mdiv = (1.0/detJ) * J + Mdiv = (1.0 / detJ) * J # Get component indices of global and reference terminal objects gtsh = g.ufl_shape @@ -261,7 +261,7 @@ if len(shape) == 0: return [None] elif len(shape) == 1: - return [None]*shape[-1] + return [None] * shape[-1] else: return [ndarray(shape[1:]) for i in range(shape[0])] global_components = ndarray(gsh) @@ -421,12 +421,12 @@ # Apply mappings with scalar indexing operations (assumes # ReferenceGrad(Jinv) is zero) - jinv_lgrad_f = lgrad[ii+jj] + jinv_lgrad_f = lgrad[ii + jj] for j, k in zip(jj, kk): - jinv_lgrad_f = Jinv[j, k]*jinv_lgrad_f + jinv_lgrad_f = Jinv[j, k] * jinv_lgrad_f # Wrap back in tensor shape, derivative axes at the end - jinv_lgrad_f = as_tensor(jinv_lgrad_f, ii+kk) + jinv_lgrad_f = as_tensor(jinv_lgrad_f, ii + kk) else: # J^(-T)RefGrad(J^(-T)RefGrad(...)) @@ -442,10 +442,10 @@ j, k = indices(2) lgrad = ReferenceGrad(jinv_lgrad_f) - jinv_lgrad_f = Jinv[j, k]*lgrad[ii+(j,)] + jinv_lgrad_f = Jinv[j, k] * lgrad[ii + (j,)] # Wrap back in tensor shape, derivative axes at the end - jinv_lgrad_f = as_tensor(jinv_lgrad_f, ii+(k,)) + jinv_lgrad_f = as_tensor(jinv_lgrad_f, ii + (k,)) return jinv_lgrad_f diff -Nru ufl-2018.1.0/ufl/algorithms/check_arities.py ufl-2019.1.0/ufl/algorithms/check_arities.py --- ufl-2018.1.0/ufl/algorithms/check_arities.py 2018-06-14 14:20:02.000000000 +0000 +++ ufl-2019.1.0/ufl/algorithms/check_arities.py 2019-04-17 18:51:30.000000000 +0000 @@ -14,6 +14,12 @@ pass +# String representation of an arity tuple: +def _afmt(atuple): + return tuple("conj({0})".format(arg) if conj else str(arg) + for arg, conj in atuple) + + class ArityChecker(MultiFunction): def __init__(self, arguments): MultiFunction.__init__(self) @@ -24,7 +30,7 @@ return self._et def argument(self, o): - return (o,) + return ((o, False),) def nonlinear_operator(self, o): # Cutoff traversal by not having *ops in argument list of this @@ -39,7 +45,7 @@ def sum(self, o, a, b): if a != b: - raise ArityMismatch("Adding expressions with non-matching form arguments {0} vs {1}.".format(a, b)) + raise ArityMismatch("Adding expressions with non-matching form arguments {0} vs {1}.".format(_afmt(a), _afmt(b))) return a def division(self, o, a, b): @@ -51,16 +57,16 @@ if a and b: # Check that we don't have test*test, trial*trial, even # for different parts in a block system - anumbers = set(x.number() for x in a) + anumbers = set(x[0].number() for x in a) for x in b: - if x.number() in anumbers: - raise ArityMismatch("Multiplying expressions with overlapping form argument number {0}, argument is {1}.".format(x.number(), x)) + if x[0].number() in anumbers: + raise ArityMismatch("Multiplying expressions with overlapping form argument number {0}, argument is {1}.".format(x[0].number(), _afmt(x))) # Combine argument lists - c = tuple(sorted(set(a + b), key=lambda x: (x.number(), x.part()))) + c = tuple(sorted(set(a + b), key=lambda x: (x[0].number(), x[0].part()))) # Check that we don't have any arguments shared between a # and b - if len(c) != len(a) + len(b): - raise ArityMismatch("Multiplying expressions with overlapping form arguments {0} vs {1}.".format(a, b)) + if len(c) != len(a) + len(b) or len(c) != len({x[0] for x in c}): + raise ArityMismatch("Multiplying expressions with overlapping form arguments {0} vs {1}.".format(_afmt(a), _afmt(b))) # It's fine for argument parts to overlap return c elif a: @@ -68,10 +74,14 @@ else: return b - # inner, outer and dot all behave as product - inner = product - outer = product - dot = product + # inner, outer and dot all behave as product but for conjugates + def inner(self, o, a, b): + return self.product(o, a, self.conj(None, b)) + + dot = inner + + def outer(self, o, a, b): + return self.product(o, self.conj(None, a), b) def linear_operator(self, o, a): return a @@ -89,6 +99,10 @@ reference_grad = linear_operator reference_value = linear_operator + # Conj, is a sesquilinear operator + def conj(self, o, a): + return tuple((a_[0], not a_[1]) for a_ in a) + # Does it make sense to have a Variable(Argument)? I see no # problem. def variable(self, o, f, l): @@ -97,7 +111,7 @@ # Conditional is linear on each side of the condition def conditional(self, o, c, a, b): if c: - raise ArityMismatch("Condition cannot depend on form arguments ({0}).".format(a)) + raise ArityMismatch("Condition cannot depend on form arguments ({0}).".format(_afmt(a))) if a and isinstance(o.ufl_operands[2], Zero): # Allow conditional(c, arg, 0) return a @@ -110,7 +124,7 @@ else: # Do not allow e.g. conditional(c, test, trial), # conditional(c, test, nonzeroconstant) - raise ArityMismatch("Conditional subexpressions with non-matching form arguments {0} vs {1}.".format(a, b)) + raise ArityMismatch("Conditional subexpressions with non-matching form arguments {0} vs {1}.".format(_afmt(a), _afmt(b))) def linear_indexed_type(self, o, a, i): return a @@ -125,28 +139,39 @@ if args: # Check that each list tensor component has the same # argument numbers (ignoring parts) - numbers = set(tuple(sorted(set(arg.number() for arg in op))) for op in ops) + numbers = set(tuple(sorted(set(arg[0].number() for arg in op))) for op in ops) if () in numbers: # Allow e.g. but not numbers.remove(()) if len(numbers) > 1: raise ArityMismatch("Listtensor components must depend on the same argument numbers, found {0}.".format(numbers)) # Allow different parts with the same number - return tuple(sorted(args, key=lambda x: (x.number(), x.part()))) + return tuple(sorted(args, key=lambda x: (x[0].number(), x[0].part()))) else: # No argument dependencies return self._et -def check_integrand_arity(expr, arguments): +def check_integrand_arity(expr, arguments, complex_mode=False): arguments = tuple(sorted(set(arguments), key=lambda x: (x.number(), x.part()))) rules = ArityChecker(arguments) - args = map_expr_dag(rules, expr, compress=False) + arg_tuples = map_expr_dag(rules, expr, compress=False) + args = tuple(a[0] for a in arg_tuples) if args != arguments: raise ArityMismatch("Integrand arguments {0} differ from form arguments {1}.".format(args, arguments)) + if complex_mode: + # Check that the test function is conjugated and that any + # trial function is not conjugated. Further arguments are + # treated as trial funtions (i.e. no conjugation) but this + # might not be correct. + for arg, conj in arg_tuples: + if arg.number() == 0 and not conj: + raise ArityMismatch("Failure to conjugate test function in complex Form") + elif arg.number() > 0 and conj: + raise ArityMismatch("Argument {0} is spuriously conjugated in complex Form".format(arg)) -def check_form_arity(form, arguments): +def check_form_arity(form, arguments, complex_mode=False): for itg in form.integrals(): - check_integrand_arity(itg.integrand(), arguments) + check_integrand_arity(itg.integrand(), arguments, complex_mode) diff -Nru ufl-2018.1.0/ufl/algorithms/comparison_checker.py ufl-2019.1.0/ufl/algorithms/comparison_checker.py --- ufl-2018.1.0/ufl/algorithms/comparison_checker.py 1970-01-01 00:00:00.000000000 +0000 +++ ufl-2019.1.0/ufl/algorithms/comparison_checker.py 2019-04-17 18:51:30.000000000 +0000 @@ -0,0 +1,132 @@ +# -*- coding: utf-8 -*- +"""Algorithm to check for 'comparison' nodes +in a form when the user is in 'complex mode'""" + +from ufl.corealg.multifunction import MultiFunction +from ufl.algorithms.map_integrands import map_integrand_dags +from ufl.algebra import Real +from ufl.constantvalue import RealValue, Zero +from ufl.argument import Argument +from ufl.geometry import GeometricQuantity + + +class CheckComparisons(MultiFunction): + """Raises an error if comparisons are done with complex quantities. + + If quantities are real, adds the Real operator to the compared quantities. + + Terminals that are real are RealValue, Zero, and Argument + (even in complex FEM, the basis functions are real) + Operations that produce reals are Abs, Real, Imag. + Terminals default to complex, and Sqrt, Pow (defensively) imply complex. + Otherwise, operators preserve the type of their operands. + """ + + def __init__(self): + MultiFunction.__init__(self) + self.nodetype = {} + + def expr(self, o, *ops): + """Defaults expressions to complex unless they only + act on real quantities. Overridden for specific operators. + + Rebuilds objects if necessary. + """ + + types = {self.nodetype[op] for op in ops} + + if types: + t = "complex" if "complex" in types else "real" + else: + t = "complex" + + o = self.reuse_if_untouched(o, *ops) + self.nodetype[o] = t + return o + + def compare(self, o, *ops): + types = {self.nodetype[op] for op in ops} + + if "complex" in types: + raise ComplexComparisonError("Ordering undefined for complex values.") + else: + o = o._ufl_expr_reconstruct_(*map(Real, ops)) + self.nodetype[o] = "bool" + return o + + gt = compare + lt = compare + ge = compare + le = compare + sign = compare + + def max_value(self, o, *ops): + types = {self.nodetype[op] for op in ops} + + if "complex" in types: + raise ComplexComparisonError("You can't compare complex numbers with max.") + else: + o = o._ufl_expr_reconstruct_(*map(Real, ops)) + self.nodetype[o] = "bool" + return o + + def min_value(self, o, *ops): + types = {self.nodetype[op] for op in ops} + + if "complex" in types: + raise ComplexComparisonError("You can't compare complex numbers with min.") + else: + o = o._ufl_expr_reconstruct_(*map(Real, ops)) + self.nodetype[o] = "bool" + return o + + def real(self, o, *ops): + o = self.reuse_if_untouched(o, *ops) + self.nodetype[o] = 'real' + return o + + def imag(self, o, *ops): + o = self.reuse_if_untouched(o, *ops) + self.nodetype[o] = 'real' + return o + + def sqrt(self, o, *ops): + o = self.reuse_if_untouched(o, *ops) + self.nodetype[o] = 'complex' + return o + + def power(self, o, base, exponent): + o = self.reuse_if_untouched(o, base, exponent) + try: + # Attempt to diagnose circumstances in which the result must be real. + exponent = float(exponent) + if self.nodetype[base] == 'real' and int(exponent) == exponent: + self.nodetype[o] = 'real' + return o + except TypeError: + pass + + self.nodetype[o] = 'complex' + return o + + def abs(self, o, *ops): + o = self.reuse_if_untouched(o, *ops) + self.nodetype[o] = 'real' + return o + + def terminal(self, term, *ops): + # default terminals to complex, except the ones we *know* are real + if isinstance(term, (RealValue, Zero, Argument, GeometricQuantity)): + self.nodetype[term] = 'real' + else: + self.nodetype[term] = 'complex' + return term + + +def do_comparison_check(form): + """Raises an error if invalid comparison nodes exist""" + return map_integrand_dags(CheckComparisons(), form) + + +class ComplexComparisonError(Exception): + pass diff -Nru ufl-2018.1.0/ufl/algorithms/compute_form_data.py ufl-2019.1.0/ufl/algorithms/compute_form_data.py --- ufl-2018.1.0/ufl/algorithms/compute_form_data.py 2018-06-14 14:20:02.000000000 +0000 +++ ufl-2019.1.0/ufl/algorithms/compute_form_data.py 2019-04-17 18:51:30.000000000 +0000 @@ -35,11 +35,13 @@ # These are the main symbolic processing steps: from ufl.algorithms.apply_function_pullbacks import apply_function_pullbacks from ufl.algorithms.apply_algebra_lowering import apply_algebra_lowering -from ufl.algorithms.apply_derivatives import apply_derivatives +from ufl.algorithms.apply_derivatives import apply_derivatives, apply_coordinate_derivatives from ufl.algorithms.apply_integral_scaling import apply_integral_scaling from ufl.algorithms.apply_geometry_lowering import apply_geometry_lowering from ufl.algorithms.apply_restrictions import apply_restrictions, apply_default_restrictions from ufl.algorithms.estimate_degrees import estimate_total_polynomial_degree +from ufl.algorithms.remove_complex_nodes import remove_complex_nodes +from ufl.algorithms.comparison_checker import do_comparison_check # See TODOs at the call sites of these below: from ufl.algorithms.domain_analysis import build_integral_data @@ -230,6 +232,8 @@ do_apply_default_restrictions=True, do_apply_restrictions=True, do_estimate_degrees=True, + do_append_everywhere_integrals=True, + complex_mode=False, ): # TODO: Move this to the constructor instead @@ -248,11 +252,23 @@ # Note: Default behaviour here will process form the way that is # currently expected by vanilla FFC + # Check that the form does not try to compare complex quantities: + # if the quantites being compared are 'provably' real, wrap them + # with Real, otherwise throw an error. + if complex_mode: + form = do_comparison_check(form) + # Lower abstractions for tensor-algebra types into index notation, # reducing the number of operators later algorithms and form # compilers need to handle form = apply_algebra_lowering(form) + # After lowering to index notation, remove any complex nodes that + # have been introduced but are not wanted when working in real mode, + # allowing for purely real forms to be written + if not complex_mode: + form = remove_complex_nodes(form) + # Apply differentiation before function pullbacks, because for # example coefficient derivatives are more complicated to derive # after coefficients are rewritten, and in particular for @@ -263,7 +279,8 @@ # TODO: Refactor this, it's rather opaque what this does # TODO: Is self.original_form.ufl_domains() right here? # It will matter when we start including 'num_domains' in ufc form. - form = group_form_integrals(form, self.original_form.ufl_domains()) + form = group_form_integrals(form, self.original_form.ufl_domains(), + do_append_everywhere_integrals=do_append_everywhere_integrals) # Estimate polynomial degree of integrands now, before applying # any pullbacks and geometric lowering. Otherwise quad degrees @@ -308,6 +325,8 @@ # Lower derivatives that may have appeared form = apply_derivatives(form) + form = apply_coordinate_derivatives(form) + # Propagate restrictions to terminals if do_apply_restrictions: form = apply_restrictions(form) @@ -391,7 +410,12 @@ # TODO: This is a very expensive check... Replace with something # faster! preprocessed_form = reconstruct_form_from_integral_data(self.integral_data) - check_form_arity(preprocessed_form, self.original_form.arguments()) # Currently testing how fast this is + + # If in real mode, remove complex nodes entirely. + if not complex_mode: + preprocessed_form = remove_complex_nodes(preprocessed_form) + + check_form_arity(preprocessed_form, self.original_form.arguments(), complex_mode) # Currently testing how fast this is # TODO: This member is used by unit tests, change the tests to # remove this! diff -Nru ufl-2018.1.0/ufl/algorithms/coordinate_derivative_helpers.py ufl-2019.1.0/ufl/algorithms/coordinate_derivative_helpers.py --- ufl-2018.1.0/ufl/algorithms/coordinate_derivative_helpers.py 1970-01-01 00:00:00.000000000 +0000 +++ ufl-2019.1.0/ufl/algorithms/coordinate_derivative_helpers.py 2019-04-17 18:51:30.000000000 +0000 @@ -0,0 +1,99 @@ +# -*- coding: utf-8 -*- +"""This module provides the necessary tools to strip away and then reattach the +coordinate derivatives at the right time point in compute_form_data.""" + +# Copyright (C) 2018 Florian Wechsung +# +# This file is part of UFL. +# +# UFL is free software: you can redistribute it and/or modify +# it under the terms of the GNU Lesser General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# UFL is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU Lesser General Public License for more details. +# +# You should have received a copy of the GNU Lesser General Public License +# along with UFL. If not, see . + +from ufl.log import error +from ufl.differentiation import CoordinateDerivative +from ufl.algorithms.multifunction import MultiFunction +from ufl.corealg.map_dag import map_expr_dags +from ufl.classes import Integral + + +class CoordinateDerivativeIsOutermostChecker(MultiFunction): + + """ Traverses the tree to make sure that CoordinateDerivatives are only on + the outside. The visitor returns False as long as no CoordinateDerivative + has been seen. """ + + def multi_index(self, o): + return False + + def terminal(self, o): + return False + + def expr(self, o, *operands): + """ If we have already seen a CoordinateDerivative, then no other + expressions apart from more CoordinateDerivatives are allowed to wrap + around it. """ + if any(operands): + error("CoordinateDerivative(s) must be outermost") + return False + + def coordinate_derivative(self, o, expr, *_): + return True + + +def strip_coordinate_derivatives(integrals): + + if isinstance(integrals, list): + if len(integrals) == 0: + return integrals, None + stripped_integrals_and_cds = [] + for integral in integrals: + (si, cd) = strip_coordinate_derivatives(integral) + stripped_integrals_and_cds.append((si, cd)) + return stripped_integrals_and_cds + + elif isinstance(integrals, Integral): + integral = integrals + integrand = integral.integrand() + checker = CoordinateDerivativeIsOutermostChecker() + map_expr_dags(checker, [integrand]) + coordinate_derivatives = [] + + # grab all coordinate derivatives and store them, so that we can apply + # them later again + def take_top_coordinate_derivatives(o): + o_ = o.ufl_operands + if isinstance(o, CoordinateDerivative): + coordinate_derivatives.append((o_[1], o_[2], o_[3])) + return take_top_coordinate_derivatives(o_[0]) + else: + return o + + newintegrand = take_top_coordinate_derivatives(integrand) + return (integral.reconstruct(integrand=newintegrand), coordinate_derivatives) + + else: + error("Invalid type %s" % (integrals.__class__.__name__,)) + + +def attach_coordinate_derivatives(integral, coordinate_derivatives): + if coordinate_derivatives is None: + return integral + + if isinstance(integral, Integral): + integrand = integral.integrand() + # apply the stored coordinate derivatives back onto the integrand + for tup in reversed(coordinate_derivatives): + integrand = CoordinateDerivative(integrand, tup[0], tup[1], tup[2]) + return integral.reconstruct(integrand=integrand) + else: + error("Invalid type %s" % (integral.__class__.__name__,)) diff -Nru ufl-2018.1.0/ufl/algorithms/domain_analysis.py ufl-2019.1.0/ufl/algorithms/domain_analysis.py --- ufl-2018.1.0/ufl/algorithms/domain_analysis.py 2018-06-14 14:20:02.000000000 +0000 +++ ufl-2019.1.0/ufl/algorithms/domain_analysis.py 2019-04-17 18:51:30.000000000 +0000 @@ -27,6 +27,7 @@ from ufl.form import Form from ufl.sorting import cmp_expr, sorted_expr from ufl.utils.sorting import canonicalize_metadata, sorted_by_key +from ufl.algorithms.coordinate_derivative_helpers import attach_coordinate_derivatives, strip_coordinate_derivatives import numbers @@ -160,7 +161,7 @@ error("Invalid domain id %s." % did) -def rearrange_integrals_by_single_subdomains(integrals): +def rearrange_integrals_by_single_subdomains(integrals, do_append_everywhere_integrals): """Rearrange integrals over multiple subdomains to single subdomain integrals. Input: @@ -196,14 +197,14 @@ otherwise_integrals = [] for ev_itg in everywhere_integrals: # Restrict everywhere integral to 'otherwise' - otherwise_integrals.append( - ev_itg.reconstruct(subdomain_id="otherwise")) + otherwise_integrals.append(ev_itg.reconstruct(subdomain_id="otherwise")) # Restrict everywhere integral to each subdomain # and append to each integral list - for subdomain_id in sorted(single_subdomain_integrals.keys()): - single_subdomain_integrals[subdomain_id].append( - ev_itg.reconstruct(subdomain_id=subdomain_id)) + if do_append_everywhere_integrals: + for subdomain_id in sorted(single_subdomain_integrals.keys()): + single_subdomain_integrals[subdomain_id].append( + ev_itg.reconstruct(subdomain_id=subdomain_id)) if otherwise_integrals: single_subdomain_integrals["otherwise"] = otherwise_integrals @@ -280,11 +281,11 @@ return integral_datas -def group_form_integrals(form, domains): +def group_form_integrals(form, domains, do_append_everywhere_integrals=True): """Group integrals by domain and type, performing canonical simplification. :arg form: the :class:`~.Form` to group the integrals of. - :arg domains: an iterable of :class:`~.Domain`\s. + :arg domains: an iterable of :class:`~.Domain`s. :returns: A new :class:`~.Form` with gathered integrands. """ # Group integrals by domain and type @@ -304,17 +305,43 @@ # (note: before this call, 'everywhere' is a valid subdomain_id, # and after this call, 'otherwise' is a valid subdomain_id) single_subdomain_integrals = \ - rearrange_integrals_by_single_subdomains(ddt_integrals) + rearrange_integrals_by_single_subdomains(ddt_integrals, do_append_everywhere_integrals) for subdomain_id, ss_integrals in sorted_by_key(single_subdomain_integrals): - # Accumulate integrands of integrals that share the - # same compiler data - integrands_and_cds = \ - accumulate_integrands_with_same_metadata(ss_integrals) - - for integrand, metadata in integrands_and_cds: - integrals.append(Integral(integrand, integral_type, domain, - subdomain_id, metadata, None)) + + # strip the coordinate derivatives from all integrals + # this yields a list of the form [(coordinate derivative, integral), ...] + stripped_integrals_and_coordderivs = strip_coordinate_derivatives(ss_integrals) + + # now group the integrals by the coordinate derivative + def calc_hash(cd): + return sum(sum(tuple_elem._ufl_compute_hash_() + for tuple_elem in tuple_) for tuple_ in cd) + coordderiv_integrals_dict = {} + for integral, coordderiv in stripped_integrals_and_coordderivs: + coordderivhash = calc_hash(coordderiv) + if coordderivhash in coordderiv_integrals_dict: + coordderiv_integrals_dict[coordderivhash][1].append(integral) + else: + coordderiv_integrals_dict[coordderivhash] = (coordderiv, [integral]) + + # cd_integrals_dict is now a dict of the form + # { hash: (CoordinateDerivative, [integral, integral, ...]), ... } + # we can now put the integrals back together and then afterwards + # apply the CoordinateDerivative again + + for cdhash, samecd_integrals in sorted_by_key(coordderiv_integrals_dict): + + # Accumulate integrands of integrals that share the + # same compiler data + integrands_and_cds = \ + accumulate_integrands_with_same_metadata(samecd_integrals[1]) + + for integrand, metadata in integrands_and_cds: + integral = Integral(integrand, integral_type, domain, + subdomain_id, metadata, None) + integral = attach_coordinate_derivatives(integral, samecd_integrals[0]) + integrals.append(integral) return Form(integrals) diff -Nru ufl-2018.1.0/ufl/algorithms/estimate_degrees.py ufl-2019.1.0/ufl/algorithms/estimate_degrees.py --- ufl-2018.1.0/ufl/algorithms/estimate_degrees.py 2018-06-14 14:20:02.000000000 +0000 +++ ufl-2019.1.0/ufl/algorithms/estimate_degrees.py 2019-04-17 18:51:30.000000000 +0000 @@ -86,7 +86,7 @@ are taken. Does not reduce the degree when TensorProduct elements or quadrilateral elements are involved.""" if isinstance(f, int) and not isinstance(f, IrreducibleInt): - return max(f-1, 0) + return max(f - 1, 0) else: # if tuple, do not reduce return f @@ -162,6 +162,15 @@ def negative_restricted(self, v, a): return a + def conj(self, v, a): + return a + + def real(self, v, a): + return a + + def imag(self, v, a): + return a + # A sum takes the max degree of its operands: sum = _max_degrees @@ -232,9 +241,9 @@ gi = g.value() if gi >= 0: if isinstance(a, int): - return a*gi + return a * gi else: - return tuple(foo*gi for foo in a) + return tuple(foo * gi for foo in a) # Something to a non-(positive integer) power, e.g. float, # negative integer, Coefficient, etc. @@ -296,6 +305,18 @@ return self._max_degrees(v, l, r) max_value = min_value + def coordinate_derivative(self, v, integrand_degree, b, direction_degree, d): + """ We use the heuristic that a shape derivative in direction V + introduces terms V and grad(V) into the integrand. Hence we add the + degree of the deformation to the estimate. """ + return self._add_degrees(v, integrand_degree, direction_degree) + + def expr_list(self, v, *o): + return self._max_degrees(v, *o) + + def expr_mapping(self, v, *o): + return self._max_degrees(v, *o) + def estimate_total_polynomial_degree(e, default_degree=1, element_replace_map={}): diff -Nru ufl-2018.1.0/ufl/algorithms/expand_indices.py ufl-2019.1.0/ufl/algorithms/expand_indices.py --- ufl-2018.1.0/ufl/algorithms/expand_indices.py 2018-06-14 14:20:02.000000000 +0000 +++ ufl-2019.1.0/ufl/algorithms/expand_indices.py 2019-04-17 18:51:30.000000000 +0000 @@ -34,6 +34,7 @@ class IndexExpander(ReuseTransformer): """...""" + def __init__(self): ReuseTransformer.__init__(self) self._components = Stack() diff -Nru ufl-2018.1.0/ufl/algorithms/formfiles.py ufl-2019.1.0/ufl/algorithms/formfiles.py --- ufl-2018.1.0/ufl/algorithms/formfiles.py 2018-06-14 14:20:02.000000000 +0000 +++ ufl-2019.1.0/ufl/algorithms/formfiles.py 2019-04-17 18:51:30.000000000 +0000 @@ -66,7 +66,7 @@ if m: encoding, = m.groups() # Drop encoding line - lines = lines[:i] + lines[i+1:] + lines = lines[:i] + lines[i + 1:] break else: # Default to utf-8 (works for ascii files diff -Nru ufl-2018.1.0/ufl/algorithms/formtransformations.py ufl-2019.1.0/ufl/algorithms/formtransformations.py --- ufl-2018.1.0/ufl/algorithms/formtransformations.py 2018-06-14 14:20:02.000000000 +0000 +++ ufl-2019.1.0/ufl/algorithms/formtransformations.py 2019-04-17 18:51:30.000000000 +0000 @@ -31,6 +31,7 @@ from ufl.argument import Argument from ufl.coefficient import Coefficient from ufl.constantvalue import Zero +from ufl.algebra import Conj # Other algorithms: from ufl.algorithms.map_integrands import map_integrands @@ -261,6 +262,11 @@ # Grad is a linear operator grad = linear_operator + # Conj, Real, Imag are linear operators + conj = linear_operator + real = linear_operator + imag = linear_operator + def linear_indexed_type(self, x): """Return parts of expression belonging to this indexed expression.""" @@ -322,7 +328,7 @@ if len(arguments) < arity: warning("Form has no parts with arity %d." % arity) - return 0*form + return 0 * form # Assuming that the form is not a sum of terms # that depend on different arguments, e.g. (u+v)*dx @@ -350,7 +356,7 @@ error("compute_form_arities cannot handle parts.") arities = set() - for arity in range(len(arguments)+1): + for arity in range(len(arguments) + 1): # Compute parts with arity "arity" parts = compute_form_with_arity(form, arity, arguments) @@ -492,4 +498,4 @@ if reordered_v.ufl_function_space() != v.ufl_function_space(): error("Element mismatch between new and old arguments (test functions).") - return replace(form, {v: reordered_v, u: reordered_u}) + return map_integrands(Conj, replace(form, {v: reordered_v, u: reordered_u})) diff -Nru ufl-2018.1.0/ufl/algorithms/__init__.py ufl-2019.1.0/ufl/algorithms/__init__.py --- ufl-2018.1.0/ufl/algorithms/__init__.py 2018-06-14 14:20:02.000000000 +0000 +++ ufl-2019.1.0/ufl/algorithms/__init__.py 2019-04-17 18:51:30.000000000 +0000 @@ -53,7 +53,6 @@ "change_to_reference_grad", "expand_compounds", "validate_form", - "ufl2latex", "FormSplitter", "extract_arguments", "compute_form_adjoint", @@ -82,7 +81,7 @@ extract_type, extract_arguments, extract_coefficients, - #extract_arguments_and_coefficients, + # extract_arguments_and_coefficients, extract_elements, extract_unique_elements, extract_sub_elements, @@ -136,4 +135,3 @@ # Utilities for UFL object printing # from ufl.formatting.printing import integral_info, form_info from ufl.formatting.printing import tree_format -from ufl.formatting.ufl2latex import ufl2latex diff -Nru ufl-2018.1.0/ufl/algorithms/multifunction.py ufl-2019.1.0/ufl/algorithms/multifunction.py --- ufl-2018.1.0/ufl/algorithms/multifunction.py 2018-06-14 14:20:02.000000000 +0000 +++ ufl-2019.1.0/ufl/algorithms/multifunction.py 2019-04-17 18:51:30.000000000 +0000 @@ -2,4 +2,4 @@ # Moved here to be usable in ufl.* files without depending on # ufl.algorithms.*... -from ufl.corealg.multifunction import MultiFunction # flake8: noqa +from ufl.corealg.multifunction import MultiFunction # noqa: F401 diff -Nru ufl-2018.1.0/ufl/algorithms/remove_complex_nodes.py ufl-2019.1.0/ufl/algorithms/remove_complex_nodes.py --- ufl-2018.1.0/ufl/algorithms/remove_complex_nodes.py 1970-01-01 00:00:00.000000000 +0000 +++ ufl-2019.1.0/ufl/algorithms/remove_complex_nodes.py 2019-04-17 18:51:30.000000000 +0000 @@ -0,0 +1,37 @@ +# -*- coding: utf-8 -*- +"""Algorithm for removing conj, real, and imag nodes +from a form for when the user is in 'real mode'""" + +from ufl.corealg.multifunction import MultiFunction +from ufl.constantvalue import ComplexValue +from ufl.algorithms.map_integrands import map_integrand_dags +from ufl.log import error + + +class ComplexNodeRemoval(MultiFunction): + """Replaces complex operator nodes with their children""" + expr = MultiFunction.reuse_if_untouched + + def conj(self, o, a): + return a + + def real(self, o, a): + return a + + def imag(self, o, a): + error("Unexpected imag in real expression.") + + def terminal(self, t, *ops): + if isinstance(t, ComplexValue): + error('Unexpected complex value in real expression.') + else: + return t + + +def remove_complex_nodes(expr): + """Replaces complex operator nodes with their children. This is called + during compute_form_data if the compiler wishes to compile + real-valued forms. In essence this strips all trace of complex + support from the preprocessed form. + """ + return map_integrand_dags(ComplexNodeRemoval(), expr) diff -Nru ufl-2018.1.0/ufl/algorithms/renumbering.py ufl-2019.1.0/ufl/algorithms/renumbering.py --- ufl-2018.1.0/ufl/algorithms/renumbering.py 2018-06-14 14:20:02.000000000 +0000 +++ ufl-2019.1.0/ufl/algorithms/renumbering.py 2019-04-17 18:51:30.000000000 +0000 @@ -44,6 +44,7 @@ class IndexRenumberingTransformer(VariableRenumberingTransformer): "This is a poorly designed algorithm. It is used in some tests, please do not use for anything else." + def __init__(self): VariableRenumberingTransformer.__init__(self) self.index_map = {} diff -Nru ufl-2018.1.0/ufl/algorithms/signature.py ufl-2019.1.0/ufl/algorithms/signature.py --- ufl-2018.1.0/ufl/algorithms/signature.py 2018-06-14 14:20:02.000000000 +0000 +++ ufl-2019.1.0/ufl/algorithms/signature.py 2019-04-17 18:51:30.000000000 +0000 @@ -36,7 +36,7 @@ j = index_numbering.get(i) if j is None: # Use negative ints for Index - j = -(len(index_numbering)+1) + j = -(len(index_numbering) + 1) index_numbering[i] = j data.append(j) else: diff -Nru ufl-2018.1.0/ufl/algorithms/transformer.py ufl-2019.1.0/ufl/algorithms/transformer.py --- ufl-2018.1.0/ufl/algorithms/transformer.py 2018-06-14 14:20:02.000000000 +0000 +++ ufl-2019.1.0/ufl/algorithms/transformer.py 2019-04-17 18:51:30.000000000 +0000 @@ -23,15 +23,16 @@ # # Modified by Anders Logg, 2009-2010 -from inspect import getargspec -from ufl.log import error -from ufl.classes import Variable, all_ufl_classes +import inspect + from ufl.algorithms.map_integrands import map_integrands +from ufl.classes import Variable, all_ufl_classes +from ufl.log import error def is_post_handler(function): "Is this a handler that expects transformed children as input?" - insp = getargspec(function) + insp = inspect.getfullargspec(function) num_args = len(insp[0]) + int(insp[1] is not None) visit_children_first = num_args > 2 return visit_children_first @@ -51,7 +52,7 @@ # first time this is run for a particular class cache_data = Transformer._handlers_cache.get(type(self)) if not cache_data: - cache_data = [None]*len(all_ufl_classes) + cache_data = [None] * len(all_ufl_classes) # For all UFL classes for classobject in all_ufl_classes: # Iterate over the inheritance chain @@ -63,27 +64,32 @@ handler_name = c._ufl_handler_name_ function = getattr(self, handler_name, None) if function: - cache_data[classobject._ufl_typecode_] = handler_name, is_post_handler(function) + cache_data[ + classobject. + _ufl_typecode_] = handler_name, is_post_handler( + function) break Transformer._handlers_cache[type(self)] = cache_data # Build handler list for this particular class (get functions # bound to self) - self._handlers = [(getattr(self, name), post) for (name, post) in cache_data] + self._handlers = [(getattr(self, name), post) + for (name, post) in cache_data] # Keep a stack of objects visit is called on, to ease # backtracking self._visit_stack = [] def print_visit_stack(self): - print("/"*80) + print("/" * 80) print("Visit stack in Transformer:") def sstr(s): ss = str(type(s)) + " ; " n = 160 - len(ss) return ss + str(s)[:n] + print("\n".join(sstr(s) for s in self._visit_stack)) - print("\\"*80) + print("\\" * 80) def visit(self, o): # Update stack @@ -224,7 +230,8 @@ def apply_transformer(e, transformer, integral_type=None): """Apply transformer.visit(expression) to each integrand expression in form, or to form if it is an Expr.""" - return map_integrands(lambda expr: transformer.visit(expr), e, integral_type) + return map_integrands(lambda expr: transformer.visit(expr), e, + integral_type) def ufl2ufl(e): diff -Nru ufl-2018.1.0/ufl/cell.py ufl-2019.1.0/ufl/cell.py --- ufl-2018.1.0/ufl/cell.py 2018-06-14 14:20:02.000000000 +0000 +++ ufl-2019.1.0/ufl/cell.py 2019-04-17 18:51:30.000000000 +0000 @@ -101,7 +101,7 @@ "hexahedron": (8, 12, 6, 1)} # Mapping from cell name to topological dimension -cellname2dim = dict((k, len(v)-1) for k, v in num_cell_entities.items()) +cellname2dim = dict((k, len(v) - 1) for k, v in num_cell_entities.items()) # Mapping from cell name to facet name # Note: This is not generalizable to product elements but it's still @@ -171,7 +171,7 @@ def num_facets(self): "The number of cell facets." tdim = self.topological_dimension() - return num_cell_entities[self.cellname()][tdim-1] + return num_cell_entities[self.cellname()][tdim - 1] # --- Facet properties --- diff -Nru ufl-2018.1.0/ufl/checks.py ufl-2019.1.0/ufl/checks.py --- ufl-2018.1.0/ufl/checks.py 2018-06-14 14:20:02.000000000 +0000 +++ ufl-2019.1.0/ufl/checks.py 2019-04-17 18:51:30.000000000 +0000 @@ -26,7 +26,7 @@ def is_python_scalar(expression): "Return True iff expression is of a Python scalar type." - return isinstance(expression, (int, float)) + return isinstance(expression, (int, float, complex)) def is_ufl_scalar(expression): diff -Nru ufl-2018.1.0/ufl/compound_expressions.py ufl-2019.1.0/ufl/compound_expressions.py --- ufl-2018.1.0/ufl/compound_expressions.py 2018-06-14 14:20:02.000000000 +0000 +++ ufl-2019.1.0/ufl/compound_expressions.py 2019-04-17 18:51:30.000000000 +0000 @@ -41,14 +41,14 @@ assert len(b) == 3 def c(i, j): - return a[i]*b[j] - a[j]*b[i] + return a[i] * b[j] - a[j] * b[i] return as_vector((c(1, 2), c(2, 0), c(0, 1))) def generic_pseudo_determinant_expr(A): """Compute the pseudo-determinant of A: sqrt(det(A.T*A)).""" i, j, k = indices(3) - ATA = as_tensor(A[k, i]*A[k, j], (i, j)) + ATA = as_tensor(A[k, i] * A[k, j], (i, j)) return sqrt(determinant_expr(ATA)) @@ -58,12 +58,12 @@ if n == 1: # Special case 1xm for simpler expression i = Index() - return sqrt(A[i, 0]*A[i, 0]) + return sqrt(A[i, 0] * A[i, 0]) elif n == 2 and m == 3: # Special case 2x3 for simpler expression c = cross_expr(A[:, 0], A[:, 1]) i = Index() - return sqrt(c[i]*c[i]) + return sqrt(c[i] * c[i]) else: # Generic formulation based on A.T*A return generic_pseudo_determinant_expr(A) @@ -72,7 +72,7 @@ def generic_pseudo_inverse_expr(A): """Compute the Penrose-Moore pseudo-inverse of A: (A.T*A)^-1 * A.T.""" i, j, k = indices(3) - ATA = as_tensor(A[k, i]*A[k, j], (i, j)) + ATA = as_tensor(A[k, i] * A[k, j], (i, j)) ATAinv = inverse_expr(ATA) q, r, s = indices(3) return as_tensor(ATAinv[r, q] * A[s, q], (r, s)) @@ -84,7 +84,7 @@ if n == 1: # Simpler special case for 1d i, j, k = indices(3) - return as_tensor(A[i, j], (j, i)) / (A[k, 0]*A[k, 0]) + return as_tensor(A[i, j], (j, i)) / (A[k, 0] * A[k, 0]) else: # Generic formulation return generic_pseudo_inverse_expr(A) @@ -110,7 +110,7 @@ def _det_2x2(B, i, j, k, l): - return B[i, k]*B[j, l] - B[i, l]*B[j, k] + return B[i, k] * B[j, l] - B[i, l] * B[j, k] def determinant_expr_2x2(B): @@ -118,9 +118,9 @@ def old_determinant_expr_3x3(A): - return (A[0, 0]*_det_2x2(A, 1, 2, 1, 2) + - A[0, 1]*_det_2x2(A, 1, 2, 2, 0) + - A[0, 2]*_det_2x2(A, 1, 2, 0, 1)) + return (A[0, 0] * _det_2x2(A, 1, 2, 1, 2) + + A[0, 1] * _det_2x2(A, 1, 2, 2, 0) + + A[0, 2] * _det_2x2(A, 1, 2, 0, 1)) def determinant_expr_3x3(A): @@ -134,7 +134,7 @@ r = rows[0] subrows = rows[1:] for i, c in enumerate(cols): - subcols = cols[i+1:] + cols[:i] + subcols = cols[i + 1:] + cols[:i] codet += A[r, c] * codeterminant_expr_nxn(A, subrows, subcols) return codet @@ -175,30 +175,30 @@ def adj_expr_3x3(A): return as_matrix([ - [A[2, 2]*A[1, 1] - A[1, 2]*A[2, 1], -A[0, 1]*A[2, 2] + A[0, 2]*A[2, 1], A[0, 1]*A[1, 2] - A[0, 2]*A[1, 1]], - [-A[2, 2]*A[1, 0] + A[1, 2]*A[2, 0], -A[0, 2]*A[2, 0] + A[2, 2]*A[0, 0], A[0, 2]*A[1, 0] - A[1, 2]*A[0, 0]], - [A[1, 0]*A[2, 1] - A[2, 0]*A[1, 1], A[0, 1]*A[2, 0] - A[0, 0]*A[2, 1], A[0, 0]*A[1, 1] - A[0, 1]*A[1, 0]], + [A[2, 2] * A[1, 1] - A[1, 2] * A[2, 1], -A[0, 1] * A[2, 2] + A[0, 2] * A[2, 1], A[0, 1] * A[1, 2] - A[0, 2] * A[1, 1]], + [-A[2, 2] * A[1, 0] + A[1, 2] * A[2, 0], -A[0, 2] * A[2, 0] + A[2, 2] * A[0, 0], A[0, 2] * A[1, 0] - A[1, 2] * A[0, 0]], + [A[1, 0] * A[2, 1] - A[2, 0] * A[1, 1], A[0, 1] * A[2, 0] - A[0, 0] * A[2, 1], A[0, 0] * A[1, 1] - A[0, 1] * A[1, 0]], ]) def adj_expr_4x4(A): return as_matrix([ - [-A[3, 3]*A[2, 1]*A[1, 2] + A[1, 2]*A[3, 1]*A[2, 3] + A[1, 1]*A[3, 3]*A[2, 2] - A[3, 1]*A[2, 2]*A[1, 3] + A[2, 1]*A[1, 3]*A[3, 2] - A[1, 1]*A[3, 2]*A[2, 3], - -A[3, 1]*A[0, 2]*A[2, 3] + A[0, 1]*A[3, 2]*A[2, 3] - A[0, 3]*A[2, 1]*A[3, 2] + A[3, 3]*A[2, 1]*A[0, 2] - A[3, 3]*A[0, 1]*A[2, 2] + A[0, 3]*A[3, 1]*A[2, 2], - A[3, 1]*A[1, 3]*A[0, 2] + A[1, 1]*A[0, 3]*A[3, 2] - A[0, 3]*A[1, 2]*A[3, 1] - A[0, 1]*A[1, 3]*A[3, 2] + A[3, 3]*A[1, 2]*A[0, 1] - A[1, 1]*A[3, 3]*A[0, 2], - A[1, 1]*A[0, 2]*A[2, 3] - A[2, 1]*A[1, 3]*A[0, 2] + A[0, 3]*A[2, 1]*A[1, 2] - A[1, 2]*A[0, 1]*A[2, 3] - A[1, 1]*A[0, 3]*A[2, 2] + A[0, 1]*A[2, 2]*A[1, 3]], - [A[3, 3]*A[1, 2]*A[2, 0] - A[3, 0]*A[1, 2]*A[2, 3] + A[1, 0]*A[3, 2]*A[2, 3] - A[3, 3]*A[1, 0]*A[2, 2] - A[1, 3]*A[3, 2]*A[2, 0] + A[3, 0]*A[2, 2]*A[1, 3], - A[0, 3]*A[3, 2]*A[2, 0] - A[0, 3]*A[3, 0]*A[2, 2] + A[3, 3]*A[0, 0]*A[2, 2] + A[3, 0]*A[0, 2]*A[2, 3] - A[0, 0]*A[3, 2]*A[2, 3] - A[3, 3]*A[0, 2]*A[2, 0], - -A[3, 3]*A[0, 0]*A[1, 2] + A[0, 0]*A[1, 3]*A[3, 2] - A[3, 0]*A[1, 3]*A[0, 2] + A[3, 3]*A[1, 0]*A[0, 2] + A[0, 3]*A[3, 0]*A[1, 2] - A[0, 3]*A[1, 0]*A[3, 2], - A[0, 3]*A[1, 0]*A[2, 2] + A[1, 3]*A[0, 2]*A[2, 0] - A[0, 0]*A[2, 2]*A[1, 3] - A[0, 3]*A[1, 2]*A[2, 0] + A[0, 0]*A[1, 2]*A[2, 3] - A[1, 0]*A[0, 2]*A[2, 3]], - [A[3, 1]*A[1, 3]*A[2, 0] + A[3, 3]*A[2, 1]*A[1, 0] + A[1, 1]*A[3, 0]*A[2, 3] - A[1, 0]*A[3, 1]*A[2, 3] - A[3, 0]*A[2, 1]*A[1, 3] - A[1, 1]*A[3, 3]*A[2, 0], - A[3, 3]*A[0, 1]*A[2, 0] - A[3, 3]*A[0, 0]*A[2, 1] - A[0, 3]*A[3, 1]*A[2, 0] - A[3, 0]*A[0, 1]*A[2, 3] + A[0, 0]*A[3, 1]*A[2, 3] + A[0, 3]*A[3, 0]*A[2, 1], - -A[0, 0]*A[3, 1]*A[1, 3] + A[0, 3]*A[1, 0]*A[3, 1] - A[3, 3]*A[1, 0]*A[0, 1] + A[1, 1]*A[3, 3]*A[0, 0] - A[1, 1]*A[0, 3]*A[3, 0] + A[3, 0]*A[0, 1]*A[1, 3], - A[0, 0]*A[2, 1]*A[1, 3] + A[1, 0]*A[0, 1]*A[2, 3] - A[0, 3]*A[2, 1]*A[1, 0] + A[1, 1]*A[0, 3]*A[2, 0] - A[1, 1]*A[0, 0]*A[2, 3] - A[0, 1]*A[1, 3]*A[2, 0]], - [-A[1, 2]*A[3, 1]*A[2, 0] - A[2, 1]*A[1, 0]*A[3, 2] + A[3, 0]*A[2, 1]*A[1, 2] - A[1, 1]*A[3, 0]*A[2, 2] + A[1, 0]*A[3, 1]*A[2, 2] + A[1, 1]*A[3, 2]*A[2, 0], - -A[3, 0]*A[2, 1]*A[0, 2] - A[0, 1]*A[3, 2]*A[2, 0] + A[3, 1]*A[0, 2]*A[2, 0] - A[0, 0]*A[3, 1]*A[2, 2] + A[3, 0]*A[0, 1]*A[2, 2] + A[0, 0]*A[2, 1]*A[3, 2], - A[0, 0]*A[1, 2]*A[3, 1] - A[1, 0]*A[3, 1]*A[0, 2] + A[1, 1]*A[3, 0]*A[0, 2] + A[1, 0]*A[0, 1]*A[3, 2] - A[3, 0]*A[1, 2]*A[0, 1] - A[1, 1]*A[0, 0]*A[3, 2], - -A[1, 1]*A[0, 2]*A[2, 0] + A[2, 1]*A[1, 0]*A[0, 2] + A[1, 2]*A[0, 1]*A[2, 0] + A[1, 1]*A[0, 0]*A[2, 2] - A[1, 0]*A[0, 1]*A[2, 2] - A[0, 0]*A[2, 1]*A[1, 2]], + [-A[3, 3] * A[2, 1] * A[1, 2] + A[1, 2] * A[3, 1] * A[2, 3] + A[1, 1] * A[3, 3] * A[2, 2] - A[3, 1] * A[2, 2] * A[1, 3] + A[2, 1] * A[1, 3] * A[3, 2] - A[1, 1] * A[3, 2] * A[2, 3], + -A[3, 1] * A[0, 2] * A[2, 3] + A[0, 1] * A[3, 2] * A[2, 3] - A[0, 3] * A[2, 1] * A[3, 2] + A[3, 3] * A[2, 1] * A[0, 2] - A[3, 3] * A[0, 1] * A[2, 2] + A[0, 3] * A[3, 1] * A[2, 2], + A[3, 1] * A[1, 3] * A[0, 2] + A[1, 1] * A[0, 3] * A[3, 2] - A[0, 3] * A[1, 2] * A[3, 1] - A[0, 1] * A[1, 3] * A[3, 2] + A[3, 3] * A[1, 2] * A[0, 1] - A[1, 1] * A[3, 3] * A[0, 2], + A[1, 1] * A[0, 2] * A[2, 3] - A[2, 1] * A[1, 3] * A[0, 2] + A[0, 3] * A[2, 1] * A[1, 2] - A[1, 2] * A[0, 1] * A[2, 3] - A[1, 1] * A[0, 3] * A[2, 2] + A[0, 1] * A[2, 2] * A[1, 3]], + [A[3, 3] * A[1, 2] * A[2, 0] - A[3, 0] * A[1, 2] * A[2, 3] + A[1, 0] * A[3, 2] * A[2, 3] - A[3, 3] * A[1, 0] * A[2, 2] - A[1, 3] * A[3, 2] * A[2, 0] + A[3, 0] * A[2, 2] * A[1, 3], + A[0, 3] * A[3, 2] * A[2, 0] - A[0, 3] * A[3, 0] * A[2, 2] + A[3, 3] * A[0, 0] * A[2, 2] + A[3, 0] * A[0, 2] * A[2, 3] - A[0, 0] * A[3, 2] * A[2, 3] - A[3, 3] * A[0, 2] * A[2, 0], + -A[3, 3] * A[0, 0] * A[1, 2] + A[0, 0] * A[1, 3] * A[3, 2] - A[3, 0] * A[1, 3] * A[0, 2] + A[3, 3] * A[1, 0] * A[0, 2] + A[0, 3] * A[3, 0] * A[1, 2] - A[0, 3] * A[1, 0] * A[3, 2], + A[0, 3] * A[1, 0] * A[2, 2] + A[1, 3] * A[0, 2] * A[2, 0] - A[0, 0] * A[2, 2] * A[1, 3] - A[0, 3] * A[1, 2] * A[2, 0] + A[0, 0] * A[1, 2] * A[2, 3] - A[1, 0] * A[0, 2] * A[2, 3]], + [A[3, 1] * A[1, 3] * A[2, 0] + A[3, 3] * A[2, 1] * A[1, 0] + A[1, 1] * A[3, 0] * A[2, 3] - A[1, 0] * A[3, 1] * A[2, 3] - A[3, 0] * A[2, 1] * A[1, 3] - A[1, 1] * A[3, 3] * A[2, 0], + A[3, 3] * A[0, 1] * A[2, 0] - A[3, 3] * A[0, 0] * A[2, 1] - A[0, 3] * A[3, 1] * A[2, 0] - A[3, 0] * A[0, 1] * A[2, 3] + A[0, 0] * A[3, 1] * A[2, 3] + A[0, 3] * A[3, 0] * A[2, 1], + -A[0, 0] * A[3, 1] * A[1, 3] + A[0, 3] * A[1, 0] * A[3, 1] - A[3, 3] * A[1, 0] * A[0, 1] + A[1, 1] * A[3, 3] * A[0, 0] - A[1, 1] * A[0, 3] * A[3, 0] + A[3, 0] * A[0, 1] * A[1, 3], + A[0, 0] * A[2, 1] * A[1, 3] + A[1, 0] * A[0, 1] * A[2, 3] - A[0, 3] * A[2, 1] * A[1, 0] + A[1, 1] * A[0, 3] * A[2, 0] - A[1, 1] * A[0, 0] * A[2, 3] - A[0, 1] * A[1, 3] * A[2, 0]], + [-A[1, 2] * A[3, 1] * A[2, 0] - A[2, 1] * A[1, 0] * A[3, 2] + A[3, 0] * A[2, 1] * A[1, 2] - A[1, 1] * A[3, 0] * A[2, 2] + A[1, 0] * A[3, 1] * A[2, 2] + A[1, 1] * A[3, 2] * A[2, 0], + -A[3, 0] * A[2, 1] * A[0, 2] - A[0, 1] * A[3, 2] * A[2, 0] + A[3, 1] * A[0, 2] * A[2, 0] - A[0, 0] * A[3, 1] * A[2, 2] + A[3, 0] * A[0, 1] * A[2, 2] + A[0, 0] * A[2, 1] * A[3, 2], + A[0, 0] * A[1, 2] * A[3, 1] - A[1, 0] * A[3, 1] * A[0, 2] + A[1, 1] * A[3, 0] * A[0, 2] + A[1, 0] * A[0, 1] * A[3, 2] - A[3, 0] * A[1, 2] * A[0, 1] - A[1, 1] * A[0, 0] * A[3, 2], + -A[1, 1] * A[0, 2] * A[2, 0] + A[2, 1] * A[1, 0] * A[0, 2] + A[1, 2] * A[0, 1] * A[2, 0] + A[1, 1] * A[0, 0] * A[2, 2] - A[1, 0] * A[0, 1] * A[2, 2] - A[0, 0] * A[2, 1] * A[1, 2]], ]) @@ -224,30 +224,30 @@ def cofactor_expr_3x3(A): return as_matrix([ - [A[1, 1]*A[2, 2] - A[2, 1]*A[1, 2], A[2, 0]*A[1, 2] - A[1, 0]*A[2, 2], - A[2, 0]*A[1, 1] + A[1, 0]*A[2, 1]], - [A[2, 1]*A[0, 2] - A[0, 1]*A[2, 2], A[0, 0]*A[2, 2] - A[2, 0]*A[0, 2], - A[0, 0]*A[2, 1] + A[2, 0]*A[0, 1]], - [A[0, 1]*A[1, 2] - A[1, 1]*A[0, 2], A[1, 0]*A[0, 2] - A[0, 0]*A[1, 2], - A[1, 0]*A[0, 1] + A[0, 0]*A[1, 1]], + [A[1, 1] * A[2, 2] - A[2, 1] * A[1, 2], A[2, 0] * A[1, 2] - A[1, 0] * A[2, 2], - A[2, 0] * A[1, 1] + A[1, 0] * A[2, 1]], + [A[2, 1] * A[0, 2] - A[0, 1] * A[2, 2], A[0, 0] * A[2, 2] - A[2, 0] * A[0, 2], - A[0, 0] * A[2, 1] + A[2, 0] * A[0, 1]], + [A[0, 1] * A[1, 2] - A[1, 1] * A[0, 2], A[1, 0] * A[0, 2] - A[0, 0] * A[1, 2], - A[1, 0] * A[0, 1] + A[0, 0] * A[1, 1]], ]) def cofactor_expr_4x4(A): return as_matrix([ - [-A[3, 1]*A[2, 2]*A[1, 3] - A[3, 2]*A[2, 3]*A[1, 1] + A[1, 3]*A[3, 2]*A[2, 1] + A[3, 1]*A[2, 3]*A[1, 2] + A[2, 2]*A[1, 1]*A[3, 3] - A[3, 3]*A[2, 1]*A[1, 2], - -A[1, 0]*A[2, 2]*A[3, 3] + A[2, 0]*A[3, 3]*A[1, 2] + A[2, 2]*A[1, 3]*A[3, 0] - A[2, 3]*A[3, 0]*A[1, 2] + A[1, 0]*A[3, 2]*A[2, 3] - A[1, 3]*A[3, 2]*A[2, 0], - A[1, 0]*A[3, 3]*A[2, 1] + A[2, 3]*A[1, 1]*A[3, 0] - A[2, 0]*A[1, 1]*A[3, 3] - A[1, 3]*A[3, 0]*A[2, 1] - A[1, 0]*A[3, 1]*A[2, 3] + A[3, 1]*A[1, 3]*A[2, 0], - A[3, 0]*A[2, 1]*A[1, 2] + A[1, 0]*A[3, 1]*A[2, 2] + A[3, 2]*A[2, 0]*A[1, 1] - A[2, 2]*A[1, 1]*A[3, 0] - A[3, 1]*A[2, 0]*A[1, 2] - A[1, 0]*A[3, 2]*A[2, 1]], - [A[3, 1]*A[2, 2]*A[0, 3] + A[0, 2]*A[3, 3]*A[2, 1] + A[0, 1]*A[3, 2]*A[2, 3] - A[3, 1]*A[0, 2]*A[2, 3] - A[0, 1]*A[2, 2]*A[3, 3] - A[3, 2]*A[0, 3]*A[2, 1], - -A[2, 2]*A[0, 3]*A[3, 0] - A[0, 2]*A[2, 0]*A[3, 3] - A[3, 2]*A[2, 3]*A[0, 0] + A[2, 2]*A[3, 3]*A[0, 0] + A[0, 2]*A[2, 3]*A[3, 0] + A[3, 2]*A[2, 0]*A[0, 3], - A[3, 1]*A[2, 3]*A[0, 0] - A[0, 1]*A[2, 3]*A[3, 0] - A[3, 1]*A[2, 0]*A[0, 3] - A[3, 3]*A[0, 0]*A[2, 1] + A[0, 3]*A[3, 0]*A[2, 1] + A[0, 1]*A[2, 0]*A[3, 3], - A[3, 2]*A[0, 0]*A[2, 1] - A[0, 2]*A[3, 0]*A[2, 1] + A[0, 1]*A[2, 2]*A[3, 0] + A[3, 1]*A[0, 2]*A[2, 0] - A[0, 1]*A[3, 2]*A[2, 0] - A[3, 1]*A[2, 2]*A[0, 0]], - [A[3, 1]*A[1, 3]*A[0, 2] - A[0, 2]*A[1, 1]*A[3, 3] - A[3, 1]*A[0, 3]*A[1, 2] + A[3, 2]*A[1, 1]*A[0, 3] + A[0, 1]*A[3, 3]*A[1, 2] - A[0, 1]*A[1, 3]*A[3, 2], - A[1, 3]*A[3, 2]*A[0, 0] - A[1, 0]*A[3, 2]*A[0, 3] - A[1, 3]*A[0, 2]*A[3, 0] + A[0, 3]*A[3, 0]*A[1, 2] + A[1, 0]*A[0, 2]*A[3, 3] - A[3, 3]*A[0, 0]*A[1, 2], - -A[1, 0]*A[0, 1]*A[3, 3] + A[0, 1]*A[1, 3]*A[3, 0] - A[3, 1]*A[1, 3]*A[0, 0] - A[1, 1]*A[0, 3]*A[3, 0] + A[1, 0]*A[3, 1]*A[0, 3] + A[1, 1]*A[3, 3]*A[0, 0], - A[0, 2]*A[1, 1]*A[3, 0] - A[3, 2]*A[1, 1]*A[0, 0] - A[0, 1]*A[3, 0]*A[1, 2] - A[1, 0]*A[3, 1]*A[0, 2] + A[3, 1]*A[0, 0]*A[1, 2] + A[1, 0]*A[0, 1]*A[3, 2]], - [A[0, 3]*A[2, 1]*A[1, 2] + A[0, 2]*A[2, 3]*A[1, 1] + A[0, 1]*A[2, 2]*A[1, 3] - A[2, 2]*A[1, 1]*A[0, 3] - A[1, 3]*A[0, 2]*A[2, 1] - A[0, 1]*A[2, 3]*A[1, 2], - A[1, 0]*A[2, 2]*A[0, 3] + A[1, 3]*A[0, 2]*A[2, 0] - A[1, 0]*A[0, 2]*A[2, 3] - A[2, 0]*A[0, 3]*A[1, 2] - A[2, 2]*A[1, 3]*A[0, 0] + A[2, 3]*A[0, 0]*A[1, 2], - -A[0, 1]*A[1, 3]*A[2, 0] + A[2, 0]*A[1, 1]*A[0, 3] + A[1, 3]*A[0, 0]*A[2, 1] - A[1, 0]*A[0, 3]*A[2, 1] + A[1, 0]*A[0, 1]*A[2, 3] - A[2, 3]*A[1, 1]*A[0, 0], - A[1, 0]*A[0, 2]*A[2, 1] - A[0, 2]*A[2, 0]*A[1, 1] + A[0, 1]*A[2, 0]*A[1, 2] + A[2, 2]*A[1, 1]*A[0, 0] - A[1, 0]*A[0, 1]*A[2, 2] - A[0, 0]*A[2, 1]*A[1, 2]] + [-A[3, 1] * A[2, 2] * A[1, 3] - A[3, 2] * A[2, 3] * A[1, 1] + A[1, 3] * A[3, 2] * A[2, 1] + A[3, 1] * A[2, 3] * A[1, 2] + A[2, 2] * A[1, 1] * A[3, 3] - A[3, 3] * A[2, 1] * A[1, 2], + -A[1, 0] * A[2, 2] * A[3, 3] + A[2, 0] * A[3, 3] * A[1, 2] + A[2, 2] * A[1, 3] * A[3, 0] - A[2, 3] * A[3, 0] * A[1, 2] + A[1, 0] * A[3, 2] * A[2, 3] - A[1, 3] * A[3, 2] * A[2, 0], + A[1, 0] * A[3, 3] * A[2, 1] + A[2, 3] * A[1, 1] * A[3, 0] - A[2, 0] * A[1, 1] * A[3, 3] - A[1, 3] * A[3, 0] * A[2, 1] - A[1, 0] * A[3, 1] * A[2, 3] + A[3, 1] * A[1, 3] * A[2, 0], + A[3, 0] * A[2, 1] * A[1, 2] + A[1, 0] * A[3, 1] * A[2, 2] + A[3, 2] * A[2, 0] * A[1, 1] - A[2, 2] * A[1, 1] * A[3, 0] - A[3, 1] * A[2, 0] * A[1, 2] - A[1, 0] * A[3, 2] * A[2, 1]], + [A[3, 1] * A[2, 2] * A[0, 3] + A[0, 2] * A[3, 3] * A[2, 1] + A[0, 1] * A[3, 2] * A[2, 3] - A[3, 1] * A[0, 2] * A[2, 3] - A[0, 1] * A[2, 2] * A[3, 3] - A[3, 2] * A[0, 3] * A[2, 1], + -A[2, 2] * A[0, 3] * A[3, 0] - A[0, 2] * A[2, 0] * A[3, 3] - A[3, 2] * A[2, 3] * A[0, 0] + A[2, 2] * A[3, 3] * A[0, 0] + A[0, 2] * A[2, 3] * A[3, 0] + A[3, 2] * A[2, 0] * A[0, 3], + A[3, 1] * A[2, 3] * A[0, 0] - A[0, 1] * A[2, 3] * A[3, 0] - A[3, 1] * A[2, 0] * A[0, 3] - A[3, 3] * A[0, 0] * A[2, 1] + A[0, 3] * A[3, 0] * A[2, 1] + A[0, 1] * A[2, 0] * A[3, 3], + A[3, 2] * A[0, 0] * A[2, 1] - A[0, 2] * A[3, 0] * A[2, 1] + A[0, 1] * A[2, 2] * A[3, 0] + A[3, 1] * A[0, 2] * A[2, 0] - A[0, 1] * A[3, 2] * A[2, 0] - A[3, 1] * A[2, 2] * A[0, 0]], + [A[3, 1] * A[1, 3] * A[0, 2] - A[0, 2] * A[1, 1] * A[3, 3] - A[3, 1] * A[0, 3] * A[1, 2] + A[3, 2] * A[1, 1] * A[0, 3] + A[0, 1] * A[3, 3] * A[1, 2] - A[0, 1] * A[1, 3] * A[3, 2], + A[1, 3] * A[3, 2] * A[0, 0] - A[1, 0] * A[3, 2] * A[0, 3] - A[1, 3] * A[0, 2] * A[3, 0] + A[0, 3] * A[3, 0] * A[1, 2] + A[1, 0] * A[0, 2] * A[3, 3] - A[3, 3] * A[0, 0] * A[1, 2], + -A[1, 0] * A[0, 1] * A[3, 3] + A[0, 1] * A[1, 3] * A[3, 0] - A[3, 1] * A[1, 3] * A[0, 0] - A[1, 1] * A[0, 3] * A[3, 0] + A[1, 0] * A[3, 1] * A[0, 3] + A[1, 1] * A[3, 3] * A[0, 0], + A[0, 2] * A[1, 1] * A[3, 0] - A[3, 2] * A[1, 1] * A[0, 0] - A[0, 1] * A[3, 0] * A[1, 2] - A[1, 0] * A[3, 1] * A[0, 2] + A[3, 1] * A[0, 0] * A[1, 2] + A[1, 0] * A[0, 1] * A[3, 2]], + [A[0, 3] * A[2, 1] * A[1, 2] + A[0, 2] * A[2, 3] * A[1, 1] + A[0, 1] * A[2, 2] * A[1, 3] - A[2, 2] * A[1, 1] * A[0, 3] - A[1, 3] * A[0, 2] * A[2, 1] - A[0, 1] * A[2, 3] * A[1, 2], + A[1, 0] * A[2, 2] * A[0, 3] + A[1, 3] * A[0, 2] * A[2, 0] - A[1, 0] * A[0, 2] * A[2, 3] - A[2, 0] * A[0, 3] * A[1, 2] - A[2, 2] * A[1, 3] * A[0, 0] + A[2, 3] * A[0, 0] * A[1, 2], + -A[0, 1] * A[1, 3] * A[2, 0] + A[2, 0] * A[1, 1] * A[0, 3] + A[1, 3] * A[0, 0] * A[2, 1] - A[1, 0] * A[0, 3] * A[2, 1] + A[1, 0] * A[0, 1] * A[2, 3] - A[2, 3] * A[1, 1] * A[0, 0], + A[1, 0] * A[0, 2] * A[2, 1] - A[0, 2] * A[2, 0] * A[1, 1] + A[0, 1] * A[2, 0] * A[1, 2] + A[2, 2] * A[1, 1] * A[0, 0] - A[1, 0] * A[0, 1] * A[2, 2] - A[0, 0] * A[2, 1] * A[1, 2]] ]) @@ -265,11 +265,11 @@ def deviatoric_expr_2x2(A): - return as_matrix([[-1./2*A[1, 1]+1./2*A[0, 0], A[0, 1]], - [A[1, 0], 1./2*A[1, 1]-1./2*A[0, 0]]]) + return as_matrix([[-1. / 2 * A[1, 1] + 1. / 2 * A[0, 0], A[0, 1]], + [A[1, 0], 1. / 2 * A[1, 1] - 1. / 2 * A[0, 0]]]) def deviatoric_expr_3x3(A): - return as_matrix([[-1./3*A[1, 1]-1./3*A[2, 2]+2./3*A[0, 0], A[0, 1], A[0, 2]], - [A[1, 0], 2./3*A[1, 1]-1./3*A[2, 2]-1./3*A[0, 0], A[1, 2]], - [A[2, 0], A[2, 1], -1./3*A[1, 1]+2./3*A[2, 2]-1./3*A[0, 0]]]) + return as_matrix([[-1. / 3 * A[1, 1] - 1. / 3 * A[2, 2] + 2. / 3 * A[0, 0], A[0, 1], A[0, 2]], + [A[1, 0], 2. / 3 * A[1, 1] - 1. / 3 * A[2, 2] - 1. / 3 * A[0, 0], A[1, 2]], + [A[2, 0], A[2, 1], -1. / 3 * A[1, 1] + 2. / 3 * A[2, 2] - 1. / 3 * A[0, 0]]]) diff -Nru ufl-2018.1.0/ufl/constantvalue.py ufl-2019.1.0/ufl/constantvalue.py --- ufl-2018.1.0/ufl/constantvalue.py 2018-06-14 14:20:02.000000000 +0000 +++ ufl-2019.1.0/ufl/constantvalue.py 2019-04-17 18:51:30.000000000 +0000 @@ -21,6 +21,8 @@ # Modified by Anders Logg, 2011. # Modified by Massimiliano Leoni, 2016. +from math import atan2 + from ufl.utils.str import as_native_str from ufl.utils.str import as_native_strings from ufl.log import error, UFLValueError @@ -192,6 +194,9 @@ def __int__(self): return 0 + def __complex__(self): + return 0 + 0j + def zero(*shape): "UFL literal constant: Return a zero tensor with the given shape." @@ -246,15 +251,67 @@ def __int__(self): return int(self._value) + def __complex__(self): + return complex(self._value) + def __neg__(self): return type(self)(-self._value) def __abs__(self): return type(self)(abs(self._value)) + def real(self): + return self._value.real + + def imag(self): + return self._value.imag + + +@ufl_type(wraps_type=complex, is_literal=True) +class ComplexValue(ScalarValue): + "UFL literal type: Representation of a constant, complex scalar" + __slots__ = () + + def __getnewargs__(self): + return (self._value,) + + def __new__(cls, value): + if value.imag == 0: + if value.real == 0: + return Zero() + else: + return FloatValue(value.real) + else: + return ConstantValue.__new__(cls) + + def __init__(self, value): + ScalarValue.__init__(self, complex(value)) + + def modulus(self): + return abs(self.value()) + + def argument(self): + return atan2(self.value().imag, self.value().real) + + def __repr__(self): + r = "%s(%s)" % (type(self).__name__, repr(self._value)) + return as_native_str(r) + + def __float__(self): + raise TypeError("ComplexValues cannot be cast to float") + + def __int__(self): + raise TypeError("ComplexValues cannot be cast to int") + + +@ufl_type(is_abstract=True, is_scalar=True) +class RealValue(ScalarValue): + "Abstract class used to differentiate real values from complex ones" + __slots__ = () + @ufl_type(wraps_type=float, is_literal=True) -class FloatValue(ScalarValue): +class FloatValue(RealValue): "UFL literal type: Representation of a constant scalar floating point value." __slots__ = () @@ -268,7 +325,7 @@ return ConstantValue.__new__(cls) def __init__(self, value): - ScalarValue.__init__(self, float(value)) + super(FloatValue, self).__init__(float(value)) def __repr__(self): r = "%s(%s)" % (type(self).__name__, format_float(self._value)) @@ -276,7 +333,7 @@ @ufl_type(wraps_type=int, is_literal=True) -class IntValue(ScalarValue): +class IntValue(RealValue): "UFL literal type: Representation of a constant scalar integer value." __slots__ = () @@ -295,15 +352,15 @@ self = IntValue._cache.get(value) if self is not None: return self - self = ScalarValue.__new__(cls) + self = RealValue.__new__(cls) IntValue._cache[value] = self else: - self = ScalarValue.__new__(cls) + self = RealValue.__new__(cls) self._init(value) return self def _init(self, value): - ScalarValue.__init__(self, int(value)) + super(IntValue, self).__init__(int(value)) def __init__(self, value): pass @@ -361,7 +418,7 @@ def __init__(self, dim): ConstantValue.__init__(self) self._dim = dim - self.ufl_shape = (dim,)*dim + self.ufl_shape = (dim,) * dim def evaluate(self, x, mapping, component, index_values): "Evaluates the permutation symbol." @@ -404,6 +461,8 @@ "Converts expression to an Expr if possible." if isinstance(expression, Expr): return expression + elif isinstance(expression, complex): + return ComplexValue(expression) elif isinstance(expression, float): return FloatValue(expression) elif isinstance(expression, int): diff -Nru ufl-2018.1.0/ufl/core/compute_expr_hash.py ufl-2019.1.0/ufl/core/compute_expr_hash.py --- ufl-2018.1.0/ufl/core/compute_expr_hash.py 2018-06-14 14:20:02.000000000 +0000 +++ ufl-2019.1.0/ufl/core/compute_expr_hash.py 2019-04-17 18:51:30.000000000 +0000 @@ -33,7 +33,7 @@ if expr._hash is not None: return expr._hash - stack = [None]*_recursion_limit_ + stack = [None] * _recursion_limit_ stacksize = 0 ops = expr.ufl_operands diff -Nru ufl-2018.1.0/ufl/core/expr.py ufl-2019.1.0/ufl/core/expr.py --- ufl-2018.1.0/ufl/core/expr.py 2018-06-14 14:20:02.000000000 +0000 +++ ufl-2019.1.0/ufl/core/expr.py 2019-04-17 18:51:30.000000000 +0000 @@ -336,6 +336,14 @@ v = NotImplemented return v + def __complex__(self): + "Try to evaluate as scalar and cast to complex." + try: + v = complex(self._ufl_evaluate_scalar_()) + except TypeError: + v = NotImplemented + return v + def __bool__(self): "By default, all Expr are nonzero/False." return True @@ -382,14 +390,6 @@ "Return a short string to represent this Expr in an error message." return "<%s id=%d>" % (self._ufl_class_.__name__, id(self)) - def _repr_latex_(self): - from ufl.algorithms import ufl2latex - return "$%s$" % ufl2latex(self) - - def _repr_png_(self): - from IPython.lib.latextools import latex_to_png - return latex_to_png(self._repr_latex_()) - # --- Special functions used for processing expressions --- def __eq__(self, other): @@ -420,7 +420,15 @@ def __round__(self, n=None): "Round to nearest integer or to nearest nth decimal." - return round(float(self), n) + try: + val = float(self._ufl_evaluate_scalar_()) + val = round(val, n) + except TypeError: + val = complex(self._ufl_evaluate_scalar_()) + val = round(val.real, n) + round(val.imag, n) * 1j + except TypeError: + val = NotImplemented + return val # --- Deprecated functions diff -Nru ufl-2018.1.0/ufl/core/terminal.py ufl-2019.1.0/ufl/core/terminal.py --- ufl-2018.1.0/ufl/core/terminal.py 2018-06-14 14:20:02.000000000 +0000 +++ ufl-2019.1.0/ufl/core/terminal.py 2019-04-17 18:51:30.000000000 +0000 @@ -57,7 +57,10 @@ # No mapping, trying to evaluate self as a constant if f is None: try: - f = float(self) + try: + f = float(self) + except TypeError: + f = complex(self) if derivatives: f = 0.0 return f diff -Nru ufl-2018.1.0/ufl/core/ufl_id.py ufl-2019.1.0/ufl/core/ufl_id.py --- ufl-2018.1.0/ufl/core/ufl_id.py 2018-06-14 14:20:02.000000000 +0000 +++ ufl-2019.1.0/ufl/core/ufl_id.py 2019-04-17 18:51:30.000000000 +0000 @@ -59,6 +59,7 @@ def _init_ufl_id(cls): "Initialize new ufl_id for the object under construction." # Bind cls with closure here + def init_ufl_id(self, ufl_id): if ufl_id is None: ufl_id = cls._ufl_global_id diff -Nru ufl-2018.1.0/ufl/corealg/map_dag.py ufl-2019.1.0/ufl/corealg/map_dag.py --- ufl-2018.1.0/ufl/corealg/map_dag.py 2018-06-14 14:20:02.000000000 +0000 +++ ufl-2019.1.0/ufl/corealg/map_dag.py 2019-04-17 18:51:30.000000000 +0000 @@ -58,8 +58,8 @@ handlers = function._handlers # Optimization else: # Regular function: no skipping supported - cutoff_types = [False]*Expr._ufl_num_typecodes_ - handlers = [function]*Expr._ufl_num_typecodes_ + cutoff_types = [False] * Expr._ufl_num_typecodes_ + handlers = [function] * Expr._ufl_num_typecodes_ # Create visited set here to share between traversal calls visited = set() diff -Nru ufl-2018.1.0/ufl/corealg/multifunction.py ufl-2019.1.0/ufl/corealg/multifunction.py --- ufl-2018.1.0/ufl/corealg/multifunction.py 2018-06-14 14:20:02.000000000 +0000 +++ ufl-2019.1.0/ufl/corealg/multifunction.py 2019-04-17 18:51:30.000000000 +0000 @@ -20,7 +20,7 @@ # # Modified by Massimiliano Leoni, 2016 -from inspect import getargspec +import inspect from ufl.log import error from ufl.core.expr import Expr @@ -28,12 +28,13 @@ def get_num_args(function): "Return the number of arguments accepted by *function*." - insp = getargspec(function) - return len(insp[0]) + int(insp[1] is not None) + sig = inspect.signature(function) + return len(sig.parameters) + 1 def memoized_handler(handler): "Function decorator to memoize ``MultiFunction`` handlers." + def _memoized_handler(self, o): c = getattr(self, "_memoized_handler_cache") r = c.get(o) @@ -66,7 +67,7 @@ algorithm_class = type(self) cache_data = MultiFunction._handlers_cache.get(algorithm_class) if not cache_data: - handler_names = [None]*len(Expr._ufl_all_classes_) + handler_names = [None] * len(Expr._ufl_all_classes_) # Iterate over the inheritance chain for each Expr # subclass (NB! This assumes that all UFL classes inherits diff -Nru ufl-2018.1.0/ufl/corealg/traversal.py ufl-2019.1.0/ufl/corealg/traversal.py --- ufl-2018.1.0/ufl/corealg/traversal.py 2018-06-14 14:20:02.000000000 +0000 +++ ufl-2019.1.0/ufl/corealg/traversal.py 2019-04-17 18:51:30.000000000 +0000 @@ -30,7 +30,7 @@ def pre_traversal(expr): """Yield ``o`` for each tree node ``o`` in *expr*, parent before child.""" - stack = [None]*_recursion_limit_ + stack = [None] * _recursion_limit_ stack[0] = expr stacksize = 1 while stacksize > 0: @@ -44,7 +44,7 @@ def post_traversal(expr): """Yield ``o`` for each node ``o`` in *expr*, child before parent.""" - stack = [None]*_recursion_limit_ + stack = [None] * _recursion_limit_ stacksize = 0 ops = expr.ufl_operands @@ -67,7 +67,7 @@ def cutoff_post_traversal(expr, cutofftypes): """Yield ``o`` for each node ``o`` in *expr*, child before parent, but skipping subtrees of the cutofftypes.""" - stack = [None]*_recursion_limit_ + stack = [None] * _recursion_limit_ stacksize = 0 ops = expr.ufl_operands @@ -96,7 +96,7 @@ This version only visits each node once. """ - stack = [None]*_recursion_limit_ + stack = [None] * _recursion_limit_ stack[0] = expr stacksize = 1 if visited is None: @@ -116,7 +116,7 @@ """Yield ``o`` for each node ``o`` in *expr*, child before parent. Never visit a node twice.""" - stack = [None]*_recursion_limit_ + stack = [None] * _recursion_limit_ stack[0] = (expr, list(expr.ufl_operands)) stacksize = 1 if visited is None: @@ -139,7 +139,7 @@ """Yield ``o`` for each node ``o`` in *expr*, child before parent. Never visit a node twice.""" - stack = [None]*_recursion_limit_ + stack = [None] * _recursion_limit_ stack[0] = (expr, () if cutofftypes[expr._ufl_typecode_] else list(expr.ufl_operands)) stacksize = 1 if visited is None: @@ -160,7 +160,7 @@ def traverse_terminals(expr): "Iterate over all terminal objects in *expr*, including duplicates." - stack = [None]*_recursion_limit_ + stack = [None] * _recursion_limit_ stack[0] = expr stacksize = 1 while stacksize > 0: @@ -176,7 +176,7 @@ def traverse_unique_terminals(expr, visited=None): "Iterate over all terminal objects in *expr*, not including duplicates." - stack = [None]*_recursion_limit_ + stack = [None] * _recursion_limit_ stack[0] = expr stacksize = 1 if visited is None: diff -Nru ufl-2018.1.0/ufl/differentiation.py ufl-2019.1.0/ufl/differentiation.py --- ufl-2018.1.0/ufl/differentiation.py 2018-06-14 14:20:02.000000000 +0000 +++ ufl-2019.1.0/ufl/differentiation.py 2019-04-17 18:51:30.000000000 +0000 @@ -80,6 +80,18 @@ self.ufl_operands[2], self.ufl_operands[3]) +@ufl_type(num_ops=4, inherit_shape_from_operand=0, + inherit_indices_from_operand=0) +class CoordinateDerivative(CoefficientDerivative): + """Derivative of the integrand of a form w.r.t. the SpatialCoordinates.""" + __slots__ = () + + def __str__(self): + return "d/dfj { %s }, with fh=%s, dfh/dfj = %s, and coordinate derivatives %s"\ + % (self.ufl_operands[0], self.ufl_operands[1], + self.ufl_operands[2], self.ufl_operands[3]) + + @ufl_type(num_ops=2) class VariableDerivative(Derivative): __slots__ = as_native_strings(( diff -Nru ufl-2018.1.0/ufl/domain.py ufl-2019.1.0/ufl/domain.py --- ufl-2018.1.0/ufl/domain.py 2018-06-14 14:20:02.000000000 +0000 +++ ufl-2019.1.0/ufl/domain.py 2019-04-17 18:51:30.000000000 +0000 @@ -43,6 +43,7 @@ and topological dimension. """ + def __init__(self, topological_dimension, geometric_dimension): # Validate dimensions if not isinstance(geometric_dimension, numbers.Integral): @@ -76,6 +77,7 @@ @attach_ufl_id class Mesh(AbstractDomain): """Symbolic representation of a mesh.""" + def __init__(self, coordinate_element, ufl_id=None, cargo=None): self._ufl_id = self._init_ufl_id(ufl_id) @@ -142,6 +144,7 @@ @attach_ufl_id class MeshView(AbstractDomain): """Symbolic representation of a mesh.""" + def __init__(self, mesh, topological_dimension, ufl_id=None): self._ufl_id = self._init_ufl_id(ufl_id) @@ -191,6 +194,7 @@ @attach_ufl_id class TensorProductMesh(AbstractDomain): """Symbolic representation of a mesh.""" + def __init__(self, meshes, ufl_id=None): self._ufl_id = self._init_ufl_id(ufl_id) @@ -268,7 +272,7 @@ if domain is None: # Create one and only one affine Mesh with a negative ufl_id # to avoid id collision - ufl_id = -(len(_default_domains)+1) + ufl_id = -(len(_default_domains) + 1) domain = affine_mesh(cell, ufl_id=ufl_id) _default_domains[cell] = domain return domain diff -Nru ufl-2018.1.0/ufl/exprequals.py ufl-2019.1.0/ufl/exprequals.py --- ufl-2018.1.0/ufl/exprequals.py 2018-06-14 14:20:02.000000000 +0000 +++ ufl-2019.1.0/ufl/exprequals.py 2019-04-17 18:51:30.000000000 +0000 @@ -28,9 +28,9 @@ # Skip those that are all not equal if sn != on and ne == tot: continue - print(fmt % (k, eq, int(100.0*eq/tot), - ne, int(100.0*ne/tot), - co, int(100.0*co/tot), + print(fmt % (k, eq, int(100.0 * eq / tot), + ne, int(100.0 * ne / tot), + co, int(100.0 * co / tot), tot)) diff -Nru ufl-2018.1.0/ufl/exproperators.py ufl-2019.1.0/ufl/exproperators.py --- ufl-2018.1.0/ufl/exproperators.py 2018-06-14 14:20:02.000000000 +0000 +++ ufl-2019.1.0/ufl/exproperators.py 2019-04-17 18:51:30.000000000 +0000 @@ -183,7 +183,7 @@ # --- Extend Expr with algebraic operators --- -_valid_types = (Expr, numbers.Real, numbers.Integral) +_valid_types = (Expr, numbers.Real, numbers.Integral, numbers.Complex) def _mul(self, o): @@ -293,7 +293,7 @@ # TODO: Add Negated class for this? Might simplify reductions in Add. def _neg(self): - return -1*self + return -1 * self Expr.__neg__ = _neg diff -Nru ufl-2018.1.0/ufl/finiteelement/elementlist.py ufl-2019.1.0/ufl/finiteelement/elementlist.py --- ufl-2018.1.0/ufl/finiteelement/elementlist.py 2018-06-14 14:20:02.000000000 +0000 +++ ufl-2019.1.0/ufl/finiteelement/elementlist.py 2019-04-17 18:51:30.000000000 +0000 @@ -125,8 +125,9 @@ (1, None), simplices[1:]) # "RTF" (2d), "N1F" (3d) # Elements not in the periodic table -register_element("Argyris", "ARG", 0, H2, "identity", (1, None), simplices[1:]) +register_element("Argyris", "ARG", 0, H2, "identity", (5, 5), ("triangle",)) register_element("Arnold-Winther", "AW", 0, H1, "identity", None, ("triangle",)) +register_element("Bell", "BELL", 0, H2, "identity", (5, 5), ("triangle",)) register_element("Brezzi-Douglas-Fortin-Marini", "BDFM", 1, HDiv, "contravariant Piola", (1, None), simplices[1:]) register_element("Crouzeix-Raviart", "CR", 0, L2, "identity", (1, 1), @@ -134,15 +135,16 @@ # TODO: Implement generic Tear operator for elements instead of this: register_element("Discontinuous Raviart-Thomas", "DRT", 1, L2, "contravariant Piola", (1, None), simplices[1:]) -register_element("Hermite", "HER", 0, H1, "identity", None, simplices[1:]) +register_element("Hermite", "HER", 0, H1, "identity", (3, 3), simplices) register_element("Mardal-Tai-Winther", "MTW", 0, H1, "identity", None, ("triangle",)) -register_element("Morley", "MOR", 0, H2, "identity", None, ("triangle",)) +register_element("Morley", "MOR", 0, H2, "identity", (2, 2), ("triangle",)) # Special elements register_element("Boundary Quadrature", "BQ", 0, L2, "identity", (0, None), any_cell) register_element("Bubble", "B", 0, H1, "identity", (2, None), simplices) +register_element("FacetBubble", "FB", 0, H1, "identity", (2, None), simplices) register_element("Quadrature", "Quadrature", 0, L2, "identity", (0, None), any_cell) register_element("Real", "R", 0, L2, "identity", (0, 0), @@ -164,6 +166,8 @@ register_alias("Lob", lambda family, dim, order, degree: ("Gauss-Lobatto-Legendre", order)) +register_element2("Bernstein", 0, H1, "identity", (1, None), simplices) + # Let Nedelec H(div) elements be aliases to BDMs/RTs register_alias("Nedelec 1st kind H(div)", diff -Nru ufl-2018.1.0/ufl/finiteelement/mixedelement.py ufl-2019.1.0/ufl/finiteelement/mixedelement.py --- ufl-2018.1.0/ufl/finiteelement/mixedelement.py 2018-06-14 14:20:02.000000000 +0000 +++ ufl-2019.1.0/ufl/finiteelement/mixedelement.py 2019-04-17 18:51:30.000000000 +0000 @@ -295,7 +295,7 @@ dim = cell.geometric_dimension() # Create list of sub elements for mixed element constructor - sub_elements = [sub_element]*dim + sub_elements = [sub_element] * dim # Compute value shapes value_shape = (dim,) + sub_element.value_shape() @@ -417,7 +417,7 @@ # Compute reference value shape based on symmetries if symmetry: # Flatten and subtract symmetries - reference_value_shape = (product(shape)-len(symmetry),) + reference_value_shape = (product(shape) - len(symmetry),) self._mapping = "symmetries" else: # Do not flatten if there are no symmetries diff -Nru ufl-2018.1.0/ufl/formatting/graph.py ufl-2019.1.0/ufl/formatting/graph.py --- ufl-2018.1.0/ufl/formatting/graph.py 2018-06-14 14:20:02.000000000 +0000 +++ ufl-2019.1.0/ufl/formatting/graph.py 2019-04-17 18:51:30.000000000 +0000 @@ -115,6 +115,7 @@ class Graph: "Graph class which computes connectivity on demand." + def __init__(self, expression): self._V, self._E = build_graph(expression) self._Ein = None @@ -214,6 +215,7 @@ "x" - depends on local coordinates "v%d" % i - depends on argument i, for i in [0,rank) """ + def __init__(self, argument_deps=None, coefficient_deps=None): MultiFunction.__init__(self) if argument_deps is None: @@ -273,7 +275,7 @@ Vout = G.Vout() partitions = defaultdict(list) - keys = [None]*n + keys = [None] * n for iv, v in enumerate(V): # Get keys from all outgoing edges edge_keys = [keys[ii] for ii in Vout[iv]] diff -Nru ufl-2018.1.0/ufl/formatting/latextools.py ufl-2019.1.0/ufl/formatting/latextools.py --- ufl-2018.1.0/ufl/formatting/latextools.py 2018-06-14 14:20:02.000000000 +0000 +++ ufl-2019.1.0/ufl/formatting/latextools.py 1970-01-01 00:00:00.000000000 +0000 @@ -1,151 +0,0 @@ -# -*- coding: utf-8 -*- -"This module defines basic utilities for stitching together LaTeX documents." - -# Copyright (C) 2008-2016 Martin Sandve Alnæs -# -# This file is part of UFL. -# -# UFL is free software: you can redistribute it and/or modify -# it under the terms of the GNU Lesser General Public License as published by -# the Free Software Foundation, either version 3 of the License, or -# (at your option) any later version. -# -# UFL is distributed in the hope that it will be useful, -# but WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -# GNU Lesser General Public License for more details. -# -# You should have received a copy of the GNU Lesser General Public License -# along with UFL. If not, see . - - -# --- Basic LaTeX tools --- - -documenttemplate = """\ -\\documentclass{article} -\\usepackage{amsmath} - -\\title{%s} - -\\begin{document} -\\maketitle{} - -%s -\\end{document} -""" - -sectiontemplate = """\ -\\section{%s} -%s -""" - -subsectiontemplate = """\ -\\subsection{%s} -%s -""" - -subsubsectiontemplate = """\ -\\subsubsection{%s} -%s -""" - -itemizetemplate = """\ -\\begin{itemize} -%s -\\end{itemize} -""" - -aligntemplate = """\ -\\begin{align} -%s -\\end{align} -""" - -verbatimtemplate = """\ -\\begin{verbatim} -%s -\\end{verbatim} -""" - - -def verbatim(string): - return verbatimtemplate % string - - -def align(lines): - # Calculate column lengths - if isinstance(lines[0], str): - body = " \\\\\n".join(l for l in lines) - else: - n = len(lines[0]) - collengths = [0]*n - for l in lines: - for i, s in enumerate(l): - collengths[i] = max(collengths[i], len(s)) - - def coljoin(cols): - return " & ".join(c.ljust(collengths[i]) for (i, c) in enumerate(cols)) - body = " \\\\\n".join(coljoin(l) for l in lines) - return aligntemplate % body - - -def itemize(items): - body = "\n".join(items) - return itemizetemplate % body - - -def subsubsection(s): - if isinstance(s, str): - return s - if isinstance(s, tuple): - title, body = s - if isinstance(body, list): - body = itemize(list(map(str, body))) - return subsubsectiontemplate % (title, body) - - -def subsection(s): - if isinstance(s, str): - return s - if isinstance(s, tuple): - title, body = s - if isinstance(body, list): - body = "\n".join(subsubsection(ss) for ss in body) - return subsectiontemplate % (title, body) - - -def section(s): - if isinstance(s, str): - return s - if isinstance(s, tuple): - title, body = s - if isinstance(body, list): - body = "\n".join(subsection(ss) for ss in body) - return sectiontemplate % (title, body) - - -def document(title, sections): - body = "\n".join(section(s) for s in sections) - return documenttemplate % (title, body) - - -def testdocument(): - title = "Test title 1" - sections = ["sec1", "sec2"] - print(document(title, sections)) - - title = "Test title 2" - sections = [("sec1", "secbody1"), ("sec2", "secbody2")] - print(document(title, sections)) - - title = "Test title 3" - section1 = [("subsec1", "subsecbody1"), ("subsec2", "subsecbody2")] - section2 = [("subsec1", "subsecbody1"), ("subsec2", "subsecbody2"), ("subsec3", "subsecbody3"), ] - section3 = "\\section{manual sec}\ntestelest" - sections = [("sec1", section1), ("sec2", section2), ("sec3", section3)] - print(document(title, sections)) - - matrix = [("a(...) ", "= \\int_\\Omega foo dx0"), - ("", "+ \\int_\\Omega foo dx1"), - ("", "+ \\int_\\Omega foo dx1"), ] - print(align(matrix)) diff -Nru ufl-2018.1.0/ufl/formatting/printing.py ufl-2019.1.0/ufl/formatting/printing.py --- ufl-2018.1.0/ufl/formatting/printing.py 2018-06-14 14:20:02.000000000 +0000 +++ ufl-2019.1.0/ufl/formatting/printing.py 2019-04-17 18:51:30.000000000 +0000 @@ -81,7 +81,7 @@ def _indent_string(n): - return " "*n + return " " * n def _tree_format_expression(expression, indentation, parentheses): @@ -89,7 +89,7 @@ if expression._ufl_is_terminal_: s = "%s%s" % (ind, repr(expression)) else: - sops = [_tree_format_expression(o, indentation+1, parentheses) for o in expression.ufl_operands] + sops = [_tree_format_expression(o, indentation + 1, parentheses) for o in expression.ufl_operands] s = "%s%s\n" % (ind, expression._ufl_class_.__name__) if parentheses and len(sops) > 1: s += "%s(\n" % (ind,) @@ -112,16 +112,16 @@ ind = _indent_string(indentation) s += ind + "Form:\n" - s += "\n".join(tree_format(itg, indentation+1, parentheses) for itg in itgs) + s += "\n".join(tree_format(itg, indentation + 1, parentheses) for itg in itgs) elif isinstance(expression, Integral): ind = _indent_string(indentation) s += ind + "Integral:\n" - ind = _indent_string(indentation+1) + ind = _indent_string(indentation + 1) s += ind + "integral type: %s\n" % expression.integral_type() s += ind + "subdomain id: %s\n" % expression.subdomain_id() s += ind + "integrand:\n" - s += tree_format(expression._integrand, indentation+2, parentheses) + s += tree_format(expression._integrand, indentation + 2, parentheses) elif isinstance(expression, Expr): s += _tree_format_expression(expression, indentation, parentheses) diff -Nru ufl-2018.1.0/ufl/formatting/ufl2dot.py ufl-2019.1.0/ufl/formatting/ufl2dot.py --- ufl-2018.1.0/ufl/formatting/ufl2dot.py 2018-06-14 14:20:02.000000000 +0000 +++ ufl-2019.1.0/ufl/formatting/ufl2dot.py 2019-04-17 18:51:30.000000000 +0000 @@ -99,6 +99,15 @@ def facet_avg(self, e): # TODO: Understandable short notation for this? return "_F_" + def conj(self, e): + return "conj" + + def real(self, e): + return "real" + + def imag(self, e): + return "imag" + def inner(self, e): return "inner" @@ -192,7 +201,7 @@ elif n > 2: oplabels = ["op%d" % i for i in range(n)] else: - oplabels = [None]*n + oplabels = [None] * n for i, o in enumerate(ops): # Handle entire subtree for expression o diff -Nru ufl-2018.1.0/ufl/formatting/ufl2latex.py ufl-2019.1.0/ufl/formatting/ufl2latex.py --- ufl-2018.1.0/ufl/formatting/ufl2latex.py 2018-06-14 14:20:02.000000000 +0000 +++ ufl-2019.1.0/ufl/formatting/ufl2latex.py 1970-01-01 00:00:00.000000000 +0000 @@ -1,972 +0,0 @@ -# -*- coding: utf-8 -*- -"""This module defines expression transformation utilities, -either converting UFL expressions to new UFL expressions or -converting UFL expressions to other representations.""" - -# Copyright (C) 2008-2016 Martin Sandve Alnæs -# -# This file is part of UFL. -# -# UFL is free software: you can redistribute it and/or modify -# it under the terms of the GNU Lesser General Public License as published by -# the Free Software Foundation, either version 3 of the License, or -# (at your option) any later version. -# -# UFL is distributed in the hope that it will be useful, -# but WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -# GNU Lesser General Public License for more details. -# -# You should have received a copy of the GNU Lesser General Public License -# along with UFL. If not, see . -# -# Modified by Anders Logg, 2008-2009. -# Modified by Kristian B. Oelgaard, 2011 - -import ufl -from ufl.log import error -from ufl.permutation import compute_indices -from ufl.algorithms.traversal import iter_expressions - -# All classes: -from ufl.variable import Variable -from ufl.core.multiindex import Index, FixedIndex -from ufl.indexed import Indexed -from ufl.tensors import ListTensor, ComponentTensor -from ufl.algebra import Sum, Product, Division, Power, Abs -from ufl.indexsum import IndexSum -from ufl.tensoralgebra import Transposed, Outer, Inner, Dot, Cross, Trace, Determinant, Inverse, Deviatoric, Cofactor -from ufl.mathfunctions import Sqrt, Exp, Ln, Cos, Sin, Tan, Cosh, Sinh, Tanh, Acos, Asin, Atan, Atan2, Erf, BesselJ, BesselY, BesselI, BesselK -from ufl.restriction import PositiveRestricted, NegativeRestricted -from ufl.averaging import CellAvg -from ufl.differentiation import VariableDerivative, Grad, Div, Curl, NablaGrad, NablaDiv -from ufl.conditional import EQ, NE, LE, GE, LT, GT, Conditional -from ufl.form import Form -from ufl.classes import terminal_classes - -# Other algorithms: -from ufl.algorithms.compute_form_data import compute_form_data - -from ufl.corealg.multifunction import MultiFunction -from ufl.corealg.map_dag import map_expr_dag -from ufl.corealg.traversal import unique_post_traversal - -from ufl.formatting.graph import build_graph, partition, extract_outgoing_vertex_connections -from ufl.formatting.latextools import align, document, verbatim - - -# TODO: Maybe this can be cleaner written using the graph utilities - -def _extract_variables(a): - """Build a list of all Variable objects in a, - which can be a Form, Integral or Expr. - The ordering in the list obeys dependency order.""" - handled = set() - variables = [] - for e in iter_expressions(a): - for o in unique_post_traversal(e): - if isinstance(o, Variable): - expr, label = o.ufl_operands - if label not in handled: - variables.append(o) - handled.add(label) - return variables - - -# --- Tools for LaTeX rendering of UFL expressions --- - -# TODO: Finish precedence mapping -def build_precedence_map(): - precedence_list = [] # TODO: Review this list very carefully! - - precedence_list.append((Sum,)) - precedence_list.append((IndexSum,)) - - # TODO: What to do with these? - precedence_list.append((ListTensor, ComponentTensor)) - precedence_list.append((CellAvg, )) - precedence_list.append((NegativeRestricted, PositiveRestricted)) - precedence_list.append((Conditional,)) - precedence_list.append((LE, GT, GE, NE, EQ, LT)) - - precedence_list.append((Div, Grad, NablaGrad, NablaDiv, Curl, - VariableDerivative, - Determinant, Trace, Cofactor, Inverse, Deviatoric)) - precedence_list.append((Product, Division, Cross, Dot, Outer, Inner)) - precedence_list.append((Indexed, Transposed, Power)) - precedence_list.append((Abs, Cos, Cosh, Exp, Ln, Sin, Sinh, Sqrt, Tan, - Tanh, Acos, Asin, Atan, Atan2, Erf, BesselJ, - BesselY, BesselI, BesselK)) - precedence_list.append((Variable,)) - precedence_list.append(terminal_classes) - - precedence_map = {} - k = 0 - for p in precedence_list: - for c in p: - precedence_map[c] = k - k += 1 - return precedence_map - - -# Utility for parentesizing string -def par(s, condition=True): # TODO: Finish precedence handling by adding condition argument to calls to this function! - if condition: - return r"\left(%s\right)" % s - return str(s) - - -def format_index(ii): - if isinstance(ii, FixedIndex): - s = "%d" % ii._value - elif isinstance(ii, Index): - s = "i_{%d}" % ii._count - else: - error("Invalid index type %s." % type(ii)) - return s - - -def format_multi_index(ii, formatstring="%s"): - return "".join(formatstring % format_index(i) for i in ii) - - -def bfname(i, p): - s = "" if p is None else (",%d" % (p,)) - return "{v_h^{%d%s}}" % (i, s) - - -def cfname(i): - return "{w_h^%d}" % i - - -# TODO: Handle line wrapping -# TODO: Handle ListTensors of rank > 1 correctly -class Expression2LatexHandler(MultiFunction): - def __init__(self, argument_names=None, coefficient_names=None): - MultiFunction.__init__(self) - self.argument_names = argument_names - self.coefficient_names = coefficient_names - - # --- Terminal objects --- - - def scalar_value(self, o): - if o.ufl_shape: - return r"{\mathbf %s}" % o._value - return "{%s}" % o._value - - def zero(self, o): - return "0" if not o.ufl_shape else r"{\mathbf 0}" - - def identity(self, o): - return r"{\mathbf I}" - - def permutation_symbol(self, o): - return r"{\mathbf \varepsilon}" - - def facet_normal(self, o): - return r"{\mathbf n}" - - def argument(self, o): - # Using ^ for argument numbering and _ for indexing since - # indexing is more common than exponentiation - if self.argument_names is None: - return bfname(o.number(), o.part()) - return self.argument_names[(o.number(), o.part())] - - def coefficient(self, o): - # Using ^ for coefficient numbering and _ for indexing since - # indexing is more common than exponentiation - if self.coefficient_names is None: - return cfname(o.count()) - return self.coefficient_names[o.count()] - constant = coefficient - - def multi_index(self, o): - return format_multi_index(o, formatstring="{%s}") - - def variable(self, o): - # TODO: Ensure variable has been handled - e, l = o.ufl_operands # noqa: E741 - return "s_{%d}" % l._count - - # --- Non-terminal objects --- - - def index_sum(self, o, f, i): - return r"\sum_{%s}%s" % (i, par(f)) - - def sum(self, o, *ops): - return " + ".join(par(op) for op in ops) - - def product(self, o, *ops): - return " ".join(par(op) for op in ops) - - def division(self, o, a, b): - return r"\frac{%s}{%s}" % (a, b) - - def abs(self, o, a): - return r"\|%s\|" % a - - def transposed(self, o, a): - return "{%s}^T" % par(a) - - def indexed(self, o, a, b): - return "{%s}_{%s}" % (par(a), b) - - def variable_derivative(self, o, f, v): - nom = r"\partial%s" % par(f) - denom = r"\partial%s" % par(v) - return r"\frac{%s}{%s}" % (nom, denom) - - def coefficient_derivative(self, o, f, w, v): - nom = r"\partial%s" % par(f) - denom = r"\partial%s" % par(w) - return r"\frac{%s}{%s}[%s]" % (nom, denom, v) # TODO: Fix this syntax... - - def grad(self, o, f): - return r"\mathbf{grad}{%s}" % par(f) - - def div(self, o, f): - return r"\mathbf{grad}{%s}" % par(f) - - def nabla_grad(self, o, f): - return r"\nabla{\otimes %s}" % par(f) - - def nabla_div(self, o, f): - return r"\nabla{\cdot %s}" % par(f) - - def curl(self, o, f): - return r"\nabla{\times %s}" % par(f) - - def sqrt(self, o, f): - return r"%s^{\frac 1 2}" % par(f) - - def exp(self, o, f): - return "e^{%s}" % f - - def ln(self, o, f): - return r"\ln{%s}" % par(f) - - def cos(self, o, f): - return r"\cos{%s}" % par(f) - - def sin(self, o, f): - return r"\sin{%s}" % par(f) - - def tan(self, o, f): - return r"\tan{%s}" % par(f) - - def cosh(self, o, f): - return r"\cosh{%s}" % par(f) - - def sinh(self, o, f): - return r"\sinh{%s}" % par(f) - - def tanh(self, o, f): - return r"\tanh{%s}" % par(f) - - def acos(self, o, f): - return r"\arccos{%s}" % par(f) - - def asin(self, o, f): - return r"\arcsin{%s}" % par(f) - - def atan(self, o, f): - return r"\arctan{%s}" % par(f) - - def atan2(self, o, f1, f2): - return r"\arctan_2{%s,%s}" % (par(f1), par(f2)) - - def erf(self, o, f): - return r"\erf{%s}" % par(f) - - def bessel_j(self, o, nu, f): - return r"J_{%s}{%s}" % (nu, par(f)) - - def bessel_y(self, o, nu, f): - return r"Y_{%s}{%s}" % (nu, par(f)) - - def bessel_i(self, o, nu, f): - return r"I_{%s}{%s}" % (nu, par(f)) - - def bessel_K(self, o, nu, f): - return r"K_{%s}{%s}" % (nu, par(f)) - - def power(self, o, a, b): - return "{%s}^{%s}" % (par(a), par(b)) - - def outer(self, o, a, b): - return r"{%s}\otimes{%s}" % (par(a), par(b)) - - def inner(self, o, a, b): - return "{%s}:{%s}" % (par(a), par(b)) - - def dot(self, o, a, b): - return r"{%s}\cdot{%s}" % (par(a), par(b)) - - def cross(self, o, a, b): - return r"{%s}\times{%s}" % (par(a), par(b)) - - def trace(self, o, A): - return "tr{%s}" % par(A) # TODO: Get built-in function syntax like \sin for this - - def determinant(self, o, A): - return "det{%s}" % par(A) # TODO: Get built-in function syntax like \sin for this - - def inverse(self, o, A): - return "{%s}^{-1}" % par(A) - - def deviatoric(self, o, A): - return "dev{%s}" % par(A) # TODO: Get built-in function syntax like \sin for this - - def cofactor(self, o, A): - return "cofac{%s}" % par(A) # TODO: Get built-in function syntax like \sin for this - - def skew(self, o, A): - return "skew{%s}" % par(A) # TODO: Get built-in function syntax like \sin for this - - def sym(self, o, A): - return "sym{%s}" % par(A) # TODO: Get built-in function syntax like \sin for this - - def list_tensor(self, o): - shape = o.ufl_shape - if len(shape) == 1: - ops = [self.visit(op) for op in o.ufl_operands] - l = " \\\\ \n ".join(ops) # noqa: E741 - elif len(shape) == 2: - rows = [] - for row in o.ufl_operands: - cols = (self.visit(op) for op in row.ufl_operands) - rows.append(" & \n ".join(cols)) - l = " \\\\ \n ".join(rows) # noqa: E741 - else: - error("TODO: LaTeX handler for list tensor of rank 3 or higher not implemented!") - return "\\left[\\begin{matrix}{%s}\\end{matrix}\\right]^T" % l - - def component_tensor(self, o, *ops): - A, ii = ops - return "\\left[A \\quad | \\quad A_{%s} = {%s} \\quad \\forall {%s} \\right]" % (ii, A, ii) - - def positive_restricted(self, o, f): - return "{%s}^+" % par(f) - - def negative_restricted(self, o, f): - return "{%s}^-" % par(f) - - def cell_avg(self, o, f): - return "{%s}_K" % par(f) - - def eq(self, o, a, b): - return "(%s = %s)" % (a, b) - - def ne(self, o, a, b): - return r"(%s \ne %s)" % (a, b) - - def le(self, o, a, b): - return r"(%s \le %s)" % (a, b) - - def ge(self, o, a, b): - return r"(%s \ge %s)" % (a, b) - - def lt(self, o, a, b): - return "(%s < %s)" % (a, b) - - def gt(self, o, a, b): - return "(%s > %s)" % (a, b) - - def and_condition(self, o, a, b): - return "(%s && %s)" % (a, b) - - def or_condition(self, o, a, b): - return "(%s || %s)" % (a, b) - - def not_condition(self, o, a): - return "!(%s)" % (a,) - - def conditional(self, o, c, t, f): - l = "\\begin{cases}\n" # noqa: E741 - l += "%s, &\text{if }\quad %s, \\\\\n" % (t, c) # noqa: E741 - l += "%s, &\text{otherwise.}\n" % f # noqa: E741 - l += "\\end{cases}" # noqa: E741 - return l - - def min_value(self, o, a, b): - return "min(%s, %s)" % (a, b) - - def max_value(self, o, a, b): - return "max(%s, %s)" % (a, b) - - def expr(self, o): - error("Missing handler for type %s" % str(type(o))) - - -def expression2latex(expression, argument_names=None, coefficient_names=None): - rules = Expression2LatexHandler(argument_names, coefficient_names) - return map_expr_dag(rules, expression) - - -def element2latex(element): - e = str(element) - e = e.replace("<", "") - e = e.replace(">", "") - e = "fixme" - return r"{\mbox{%s}}" % e - - -domain_strings = {"cell": r"\Omega", - "exterior_facet": r"\Gamma^{ext}", - "exterior_facet_bottom": r"\Gamma_{bottom}^{ext}", - "exterior_facet_top": r"\Gamma_{top}^{ext}", - "exterior_facet_vert": r"\Gamma_{vert}^{ext}", - "interior_facet": r"\Gamma^{int}", - "interior_facet_horiz": r"\Gamma_{horiz}^{int}", - "interior_facet_vert": r"\Gamma_{vert}^{int}", - "vertex": r"\Gamma^{vertex}", - "custom": r"\Gamma^{custom}", } -default_domain_string = "d(?)" - - -def form2latex(form, formdata): - - formname = formdata.name - argument_names = formdata.argument_names - coefficient_names = formdata.coefficient_names - - # List of sections to make latex document from - sections = [] - - # Define elements - lines = [] - for i, f in enumerate(formdata.original_arguments): - lines.append(r"\mathcal{P}_{%d} = \{%s\} " % (i, element2latex(f.ufl_element()))) - for i, f in enumerate(formdata.original_coefficients): - lines.append(r"\mathcal{Q}_{%d} = \{%s\} " % (i, element2latex(f.ufl_element()))) - if lines: - sections.append(("Finite elements", align(lines))) - - # Define function spaces - lines = [] - for i, f in enumerate(formdata.original_arguments): - lines.append("V_h^{%d} = \\{v : v \\vert_K \\in \\mathcal{P}_{%d}(K) \\quad \\forall K \\in \\mathcal{T}\\} " % (i, i)) - for i, f in enumerate(formdata.original_coefficients): - lines.append("W_h^{%d} = \\{v : v \\vert_K \\in \\mathcal{Q}_{%d}(K) \\quad \\forall K \\in \\mathcal{T}\\} " % (i, i)) - if lines: - sections.append(("Function spaces", align(lines))) - - # Define arguments and coefficients - lines = [] - for f in formdata.original_arguments: - i = f.number() - p = f.part() - lines.append("%s = %s \\in V_h^{%d} " % (argument_names[(i, p)], bfname(i, p), i)) # FIXME: Handle part in V_h - for i, f in enumerate(formdata.original_coefficients): - lines.append("%s = %s \\in W_h^{%d} " % (coefficient_names[i], cfname(i), i)) - if lines: - sections.append(("Form arguments", align(lines))) - - # TODO: Wrap ListTensors, ComponentTensor and Conditionals in - # expression as variables before transformation - - # Define variables - handled_variables = set() - integrals = form.integrals() - lines = [] - for itg in integrals: - variables = _extract_variables(itg.integrand()) - for v in variables: - l = v._label # noqa: E741 - if l not in handled_variables: - handled_variables.add(l) - exprlatex = expression2latex(v._expression, - formdata.argument_names, - formdata.coefficient_names) - lines.append(("s_{%d}" % l._count, "= %s" % exprlatex)) - if lines: - sections.append(("Variables", align(lines))) - - # Join form arguments for signature "a(...) =" - b = ", ".join(formdata.argument_names) - c = ", ".join(formdata.coefficient_names) - arguments = "; ".join((b, c)) - signature = "%s(%s) = " % (formname, arguments, ) - - # Define form as sum of integrals - lines = [] - a = signature - p = "" - for itg in integrals: - # TODO: Get list of expression strings instead of single - # expression! - integrand_string = expression2latex(itg.integrand(), - formdata.argument_names, - formdata.coefficient_names) - - integral_type = itg.integral_type() - dstr = domain_strings[integral_type] - - # domain = itg.ufl_domain() - # TODO: Render domain description - - subdomain_id = itg.subdomain_id() - if isinstance(subdomain_id, int): - dstr += "_{%d}" % subdomain_id - elif subdomain_id == "everywhere": - pass - elif subdomain_id == "otherwise": - dstr += "_{\text{oth}}" - elif isinstance(subdomain_id, tuple): - dstr += "_{%s}" % subdomain_id - - b = p + "\\int_{%s}" % (dstr,) - dxstr = ufl.measure.integral_type_to_measure_name[integral_type] - c = "{ %s } \\,%s" % (integrand_string, dxstr) - lines.append((a, b, c)) - a = "{}" - p = "{}+ " - - sections.append(("Form", align(lines))) - - return sections - - -def ufl2latex(expression): - "Generate LaTeX code for a UFL expression or form (wrapper for form2latex and expression2latex)." - if isinstance(expression, Form): - form_data = compute_form_data(expression) - preprocessed_form = form_data.preprocessed_form - return form2latex(preprocessed_form, form_data) - else: - return expression2latex(expression) - - -# --- LaTeX rendering of composite UFL objects --- - -def deps2latex(deps): - return "Dependencies: ${ %s }$." % ", ".join(sorted(deps)) - - -def dependency_sorting(deplist, rank): - - def split(deps, state): - left = [] - todo = [] - for dep in deps: - if dep - state: - left.append(dep) - else: - todo.append(dep) - return todo, left - - deplistlist = [] - state = set() - left = deplist - - # --- Initialization time - precompute, left = split(left, state) - deplistlist.append(precompute) - - state.add("x") - precompute_quad, left = split(left, state) - deplistlist.append(precompute_quad) - - # Permutations of 0/1 dependence of arguments - indices = compute_indices((2,)*rank) - for bfs in indices[1:]: # skip (0,...,0), already handled that - for i, bf in reversed(enumerate(bfs)): - n = "v%d" % i - if bf: - if n in state: - state.remove(n) - else: - state.add(n) - next, left = split(left, state) - deplistlist.append(next) - - # --- Runtime - state.add("c") - state.add("w") - - state.remove("x") - runtime, left = split(left, state) - deplistlist.append(runtime) - - state.add("x") - runtime_quad, left = split(left, state) - deplistlist.append(runtime_quad) - - indices = compute_indices((2,)*rank) - for bfs in indices[1:]: # skip (0,...,0), already handled that - for i, bf in reversed(enumerate(bfs)): - n = "v%d" % i - if bf: - state.add(n) - else: - if n in state: - state.remove(n) - next, left = split(left, state) - deplistlist.append(next) - - if left: - error("Shouldn't have anything left!") - - return deplistlist - - -def code2latex(G, partitions, formdata): - "TODO: Document me" - bfn = formdata.argument_names - - V, E = G - Vout = extract_outgoing_vertex_connections(G) - - # Sort dependency sets in a sensible way (preclude to a good - # quadrature code generator) - deplistlist = dependency_sorting(list(partitions.keys()), len(bfn)) - - def format_v(i): - return "s_{%d}" % i - - pieces = [] - for deplist in deplistlist: - for dep in deplist: - lines = [] - for iv in partitions[dep]: - v = V[iv] - vout = Vout[iv] - vl = format_v(iv) - args = ", ".join(format_v(i) for i in vout) - if args: - el = r"{\mbox{%s}}(%s)" % (v._ufl_class_.__name__, args) - else: # terminal - el = r"{\mbox{%s}}" % (repr(v),) - lines.append((vl, "= " + el)) - pieces.extend(("\n", deps2latex(dep), align(lines))) - - # Add final variable representing integrand - vl = format_v(len(V)-1) - pieces.append("\n") - pieces.append("Variable representing integrand: %s" % vl) - - # Could also return list of (title, body) parts for subsections if - # wanted - body = "\n".join(pieces) - return body - - -def integrand2code(integrand, formdata): - G = build_graph(integrand) - partitions, keys = partition(G) - return G, partitions - - -def formdata2latex(formdata): # TODO: Format better - return verbatim(str(formdata)) - - -def form2code2latex(formdata): - # Render introductory sections - title = "Form data" - body = formdata2latex(formdata) - sections = [(title, body)] - - # Render each integral as a separate section - for itg in formdata.form.integrals(): - title = "%s integral over domain %d" % (itg.integral_type(), - itg.subdomain_id()) - - G, partitions = integrand2code(itg.integrand(), formdata) - - body = code2latex(G, partitions, formdata) - - sections.append((title, body)) - - return sections - - -# --- Creating complete documents --- - -def forms2latexdocument(forms, uflfilename, compile=False): - "Render forms from a .ufl file as a LaTeX document." - - # Render one section for each form - sections = [] - for form in forms: - - # Compute form data - form_data = compute_form_data(form) - - # Generate LaTex code - title = "Form %s" % form_data.name - if compile: - body = form2code2latex(form, form_data) - else: - body = form2latex(form, form_data) - sections.append((title, body)) - - # Render title - suffix = "from UFL file %s" % uflfilename.replace("_", "\\_") - if compile: - title = "Compiled forms " + suffix - else: - title = "Forms " + suffix - return document(title, sections) - - -"""# Code from uflacs: - -from ffc.log import error -from ffc.log import ffc_assert - -import ufl - -from ufl.corealg.multifunction import MultiFunction - -# TODO: Assuming in this code that preprocessed expressions -# are formatted, so no compounds etc. are included here. -# Would be nice to format e.g. dot(u, v) -> u \cdot v. - - -class LatexFormattingRules(object): - - # === Error rules catching groups of missing types by their superclasses === - - # Generic fallback error messages for missing rules: - def expr(self, o): - error("Missing LaTeX formatting rule for expr type %s." % o._ufl_class_) - - def terminal(self, o): - error("Missing LaTeX formatting rule for terminal type %s." % o._ufl_class_) - - def constant_value(self, o, component=(), derivatives=(), restriction=None): - error("Missing LaTeX rule for constant value type %s." % o._ufl_class_) - - def geometric_quantity(self, o, component=(), derivatives=()): - error("Missing LaTeX formatting rule for geometric quantity type %s." % o._ufl_class_) - - # Unexcepted type checks: - def variable(self, o): - error("Should strip away variables before formatting LaTeX code.") - return o # or just do this if necessary - - def invalid_request(self, o, *ops): - error("Invalid request for LaTeX formatting of a %s." % o._ufl_class_) - wrapper_type = invalid_request - index_sum = invalid_request - indexed = invalid_request - derivative = invalid_request - restricted = invalid_request - - # === Formatting rules for literal constants === - - def zero(self, o, component=(), derivatives=(), restriction=None): - return "0" if not o.ufl_shape else r"{\mathbf 0}" - - def int_value(self, o, component=(), derivatives=(), restriction=None): - if derivatives: - return self.zero(0 * o) - else: - return "%d" % int(o) - - def float_value(self, o, component=(), derivatives=(), restriction=None): - # Using configurable precision parameter from ufl - if derivatives: - return self.zero(0 * o) - else: - return ufl.constantvalue.format_float(float(o)) - - # ... The compound literals below are removed during preprocessing - - def identity(self, o): - return r"{\mathbf I}" - - def permutation_symbol(self, o): - return r"{\mathbf \varepsilon}" - - # === Formatting rules for geometric quantities === - - # TODO: Add all geometric quantities here, use restriction - - def spatial_coordinate(self, o, component=(), derivatives=(), restriction=None): - if component: - i, = component - else: - i = 0 - if derivatives: - return "x_{%d, %s}" % (i, ' '.join('%d' % d for d in derivatives)) - else: - return "x_%d" % i - - def facet_normal(self, o, component=(), derivatives=(), restriction=None): - if component: - i, = component - else: - i = 0 - if derivatives: - return "n_{%d, %s}" % (i, ' '.join('%d' % d for d in derivatives)) - else: - return "n_%d" % i - - def cell_volume(self, o, component=(), derivatives=(), restriction=None): - ffc_assert(not component, "Expecting no component for scalar value.") - if derivatives: - return "0" - else: - return r"K_{\text{vol}}" - - def circumradius(self, o, component=(), derivatives=(), restriction=None): - ffc_assert(not component, "Expecting no component for scalar value.") - if derivatives: - return "0" - else: - return r"K_{\text{rad}}" - - # === Formatting rules for functions === - - def coefficient(self, o, component=(), derivatives=(), restriction=None): - common_name = "w" - c = o.count() - - ffc_assert(c >= 0, "Expecting positive count, have you preprocessed the expression?") - - name = r"\overset{%d}{%s}" % (c, common_name) - - # TODO: Use restriction - - if component: - cstr = ' '.join('%d' % d for d in component) - else: - cstr = '' - - if derivatives: - dstr = ' '.join('%d' % d for d in derivatives) - return "%s_{%s, %s}" % (name, cstr, dstr) - elif not component: - return name - else: - return "%s_{%s}" % (name, cstr) - - def argument(self, o, component=(), derivatives=(), restriction=None): - common_name = "v" - c = o.number() - - name = r"\overset{%d}{%s}" % (c, common_name) - - # TODO: Use restriction - - if component: - cstr = ' '.join('%d' % d for d in component) - else: - cstr = '' - - if derivatives: - dstr = ' '.join('%d' % d for d in derivatives) - return "%s_{%s, %s}" % (name, cstr, dstr) - elif not component: - return name - else: - return "%s_{%s}" % (name, cstr) - - # === Formatting rules for arithmetic operations === - - def sum(self, o, *ops): - return " + ".join(ops) - - def product(self, o, *ops): - return " ".join(ops) - - def division(self, o, a, b): - return r"\frac{%s}{%s}" % (a, b) - - # === Formatting rules for cmath functions === - - def power(self, o, a, b): - return "{%s}^{%s}" % (a, b) - - def sqrt(self, o, op): - return "\sqrt{%s}" % (op,) - - def ln(self, o, op): - return r"\ln(%s)" % (op,) - - def exp(self, o, op): - return "e^{%s}" % (op,) - - def abs(self, o, op): - return r"\|%s\|" % (op,) - - def cos(self, o, op): - return r"\cos(%s)" % (op,) - - def sin(self, o, op): - return r"\sin(%s)" % (op,) - - def tan(self, o, op): - return r"\tan(%s)" % (op,) - - def cosh(self, o, op): - return r"\cosh(%s)" % (op,) - - def sinh(self, o, op): - return r"\sinh(%s)" % (op,) - - def tanh(self, o, op): - return r"\tanh(%s)" % (op,) - - def acos(self, o, op): - return r"\arccos(%s)" % (op,) - - def asin(self, o, op): - return r"\arcsin(%s)" % (op,) - - def atan(self, o, op): - return r"\arctan(%s)" % (op,) - - # === Formatting rules for bessel functions === - - # TODO: Bessel functions, erf - - # === Formatting rules for conditional expressions === - - def conditional(self, o, c, t, f): - return r"\left{{%s} \text{if} {%s} \text{else} {%s}\right}" % (t, c, f) - - def eq(self, o, a, b): - return r" = ".join((a, b)) - - def ne(self, o, a, b): - return r" \ne ".join((a, b)) - - def le(self, o, a, b): - return r" \le ".join((a, b)) - - def ge(self, o, a, b): - return r" \ge ".join((a, b)) - - def lt(self, o, a, b): - return r" \lt ".join((a, b)) - - def gt(self, o, a, b): - return r" \gt ".join((a, b)) - - def and_condition(self, o, a, b): - return r" \land ".join((a, b)) - - def or_condition(self, o, a, b): - return r" \lor ".join((a, b)) - - def not_condition(self, o, a): - return r" \lnot %s" % (a,) - - # === Formatting rules for restrictions === - - def positive_restricted(self, o, a): - return r"%s^{[+]}" % (a,) # TODO - - def negative_restricted(self, o, a): - return r"%s^{[-]}" % (a,) # TODO - - -class LatexFormatter(MultiFunction, LatexFormattingRules): - def __init__(self): - MultiFunction.__init__(self) - -""" diff -Nru ufl-2018.1.0/ufl/formatting/ufl2unicode.py ufl-2019.1.0/ufl/formatting/ufl2unicode.py --- ufl-2018.1.0/ufl/formatting/ufl2unicode.py 2018-06-14 14:20:02.000000000 +0000 +++ ufl-2019.1.0/ufl/formatting/ufl2unicode.py 2019-04-17 18:51:30.000000000 +0000 @@ -13,6 +13,7 @@ class PrecedenceRules(MultiFunction): "An enum-like class for C operator precedence levels." + def __init__(self): MultiFunction.__init__(self) @@ -617,6 +618,16 @@ def sym(self, o, A): return mathop(o, A, "sym") + def conj(self, o, a): + # Overbar is already taken for average, and there is no superscript asterix in unicode. + return mathop(o, a, "conj") + + def real(self, o, a): + return mathop(o, a, "Re") + + def imag(self, o, a): + return mathop(o, a, "Im") + def list_tensor(self, o, *ops): l = ", ".join(ops) # noqa: E741 return "%s%s%s" % ("[", l, "]") diff -Nru ufl-2018.1.0/ufl/formoperators.py ufl-2019.1.0/ufl/formoperators.py --- ufl-2018.1.0/ufl/formoperators.py 2018-06-14 14:20:02.000000000 +0000 +++ ufl-2019.1.0/ufl/formoperators.py 2019-04-17 18:51:30.000000000 +0000 @@ -30,13 +30,14 @@ from ufl.finiteelement import MixedElement from ufl.argument import Argument from ufl.coefficient import Coefficient -from ufl.differentiation import CoefficientDerivative +from ufl.differentiation import CoefficientDerivative, CoordinateDerivative from ufl.constantvalue import is_true_ufl_scalar, as_ufl from ufl.indexed import Indexed from ufl.core.multiindex import FixedIndex, MultiIndex from ufl.tensors import as_tensor, ListTensor from ufl.sorting import sorted_expr from ufl.functionspace import FunctionSpace +from ufl.geometry import SpatialCoordinate # An exception to the rule that ufl.* does not depend on ufl.algorithms.* ... from ufl.algorithms import compute_form_adjoint, compute_form_action @@ -129,7 +130,8 @@ def adjoint(form, reordered_arguments=None): """UFL form operator: Given a combined bilinear form, compute the adjoint form by - changing the ordering (count) of the test and trial functions. + changing the ordering (count) of the test and trial functions, and + taking the complex conjugate of the result. By default, new ``Argument`` objects will be created with opposite ordering. However, if the adjoint form is to @@ -145,7 +147,7 @@ if len(shape) == 0: error("Invalid shape.") elif len(shape) == 1: - return [0]*shape[0] + return [0] * shape[0] else: return [zero_lists(shape[1:]) for i in range(shape[0])] @@ -221,7 +223,7 @@ for (c, a) in zip(coefficients, arguments): if c.ufl_shape != a.ufl_shape: error("Coefficient and argument shapes do not match!") - if isinstance(c, Coefficient): + if isinstance(c, Coefficient) or isinstance(c, SpatialCoordinate): m[c] = a else: if not isinstance(c, Indexed): @@ -290,15 +292,23 @@ if isinstance(form, Form): integrals = [] for itg in form.integrals(): - fd = CoefficientDerivative(itg.integrand(), coefficients, - arguments, coefficient_derivatives) + if not isinstance(coefficient, SpatialCoordinate): + fd = CoefficientDerivative(itg.integrand(), coefficients, + arguments, coefficient_derivatives) + else: + fd = CoordinateDerivative(itg.integrand(), coefficients, + arguments, coefficient_derivatives) integrals.append(itg.reconstruct(fd)) return Form(integrals) elif isinstance(form, Expr): # What we got was in fact an integrand - return CoefficientDerivative(form, coefficients, arguments, - coefficient_derivatives) + if not isinstance(coefficient, SpatialCoordinate): + return CoefficientDerivative(form, coefficients, + arguments, coefficient_derivatives) + else: + return CoordinateDerivative(form, coefficients, + arguments, coefficient_derivatives) error("Invalid argument type %s." % str(type(form))) diff -Nru ufl-2018.1.0/ufl/form.py ufl-2019.1.0/ufl/form.py --- ufl-2018.1.0/ufl/form.py 2018-06-14 14:20:02.000000000 +0000 +++ ufl-2019.1.0/ufl/form.py 2019-04-17 18:51:30.000000000 +0000 @@ -46,11 +46,14 @@ # Group integrals in multilevel dict by keys # [domain][integral_type][subdomain_id] - integrals_dict = defaultdict(lambda: defaultdict(lambda: defaultdict(list))) + integrals_dict = defaultdict( + lambda: defaultdict(lambda: defaultdict(list))) for integral in integrals: d = integral.ufl_domain() if d is None: - error("Each integral in a form must have a uniquely defined integration domain.") + error( + "Each integral in a form must have a uniquely defined integration domain." + ) it = integral.integral_type() si = integral.subdomain_id() integrals_dict[d][it][si] += [integral] @@ -60,8 +63,9 @@ # Order integrals canonically to increase signature stability for d in sort_domains(integrals_dict): for it in sorted(integrals_dict[d]): # str is sortable - for si in sorted(integrals_dict[d][it], - key=lambda x: (type(x).__name__, x)): # int/str are sortable + for si in sorted( + integrals_dict[d][it], key=lambda x: (type(x).__name__, x) + ): # int/str are sortable unsorted_integrals = integrals_dict[d][it][si] # TODO: At this point we could order integrals by # metadata and integrand, or even add the @@ -175,17 +179,20 @@ # Check that all are equal TODO: don't return more than one if # all are equal? if not all(domain == domains[0] for domain in domains): - error("Calling Form.ufl_domain() is only valid if all integrals share domain.") + error( + "Calling Form.ufl_domain() is only valid if all integrals share domain." + ) # Return the one and only domain return domains[0] def geometric_dimension(self): "Return the geometric dimension shared by all domains and functions in this form." - gdims = tuple(set(domain.geometric_dimension() - for domain in self.ufl_domains())) + gdims = tuple( + set(domain.geometric_dimension() for domain in self.ufl_domains())) if len(gdims) != 1: error("Expecting all domains and functions in a form " - "to share geometric dimension, got %s." % str(tuple(sorted(gdims)))) + "to share geometric dimension, got %s." % str( + tuple(sorted(gdims)))) return gdims[0] def domain_numbering(self): @@ -278,7 +285,9 @@ # Allow adding 0 or 0.0 as a no-op, needed for sum([a,b]) return self - elif isinstance(other, Zero) and not (other.ufl_shape or other.ufl_free_indices): + elif isinstance( + other, + Zero) and not (other.ufl_shape or other.ufl_free_indices): # Allow adding ufl Zero as a no-op, needed for sum([a,b]) return self @@ -305,7 +314,7 @@ "Multiply all integrals in form with constant scalar value." # This enables the handy "0*form" or "dt*form" syntax if is_scalar_constant_expression(scalar): - return Form([scalar*itg for itg in self.integrals()]) + return Form([scalar * itg for itg in self.integrals()]) return NotImplemented def __mul__(self, coefficient): @@ -344,7 +353,8 @@ if args: arguments = self.arguments() if len(arguments) != len(args): - error("Need %d arguments to form(), got %d." % (len(arguments), len(args))) + error("Need %d arguments to form(), got %d." % (len(arguments), + len(args))) repdict.update(zip(arguments, args)) coefficients = kwargs.pop("coefficients") @@ -386,28 +396,22 @@ r = "Form([" + itgs + "])" return as_native_str(r) - def x_repr_latex_(self): # TODO: This works, but enable when form latex rendering is fixed - from ufl.algorithms import ufl2latex - return "$$%s$$" % ufl2latex(self) - - def x_repr_png_(self): # TODO: This works, but enable when form latex rendering is fixed - from IPython.lib.latextools import latex_to_png - return latex_to_png(self._repr_latex_()) - # --- Analysis functions, precomputation and caching of various quantities def _analyze_domains(self): from ufl.domain import join_domains, sort_domains # Collect unique integration domains - integration_domains = join_domains([itg.ufl_domain() for itg in self._integrals]) + integration_domains = join_domains( + [itg.ufl_domain() for itg in self._integrals]) # Make canonically ordered list of the domains self._integration_domains = sort_domains(integration_domains) # TODO: Not including domains from coefficients and arguments # here, may need that later - self._domain_numbering = dict((d, i) for i, d in enumerate(self._integration_domains)) + self._domain_numbering = dict( + (d, i) for i, d in enumerate(self._integration_domains)) def _analyze_subdomain_data(self): integration_domains = self.ufl_domains() @@ -430,7 +434,9 @@ subdomain_data[domain][it] = sd elif sd is not None: if data.ufl_id() != sd.ufl_id(): - error("Integrals in form have different subdomain_data objects.") + error( + "Integrals in form have different subdomain_data objects." + ) self._subdomain_data = subdomain_data def _analyze_form_arguments(self): @@ -439,12 +445,12 @@ arguments, coefficients = extract_arguments_and_coefficients(self) # Define canonical numbering of arguments and coefficients - self._arguments = tuple(sorted(set(arguments), - key=lambda x: x.number())) - self._coefficients = tuple(sorted(set(coefficients), - key=lambda x: x.count())) - self._coefficient_numbering = dict((c, i) for i, - c in enumerate(self._coefficients)) + self._arguments = tuple( + sorted(set(arguments), key=lambda x: x.number())) + self._coefficients = tuple( + sorted(set(coefficients), key=lambda x: x.count())) + self._coefficient_numbering = dict( + (c, i) for i, c in enumerate(self._coefficients)) def _compute_renumbering(self): # Include integration domains and coefficients in renumbering diff -Nru ufl-2018.1.0/ufl/geometry.py ufl-2019.1.0/ufl/geometry.py --- ufl-2018.1.0/ufl/geometry.py 2018-06-14 14:20:02.000000000 +0000 +++ ufl-2019.1.0/ufl/geometry.py 2019-04-17 18:51:30.000000000 +0000 @@ -182,6 +182,11 @@ else: return float(x[component[0]]) + def count(self): + # FIXME: Hack to make SpatialCoordinate behave like a coefficient. + # When calling `derivative`, the count is used to sort over. + return -1 + @ufl_type() class CellCoordinate(GeometricCellQuantity): @@ -230,7 +235,7 @@ @property def ufl_shape(self): t = self._domain.topological_dimension() - return (t-1,) + return (t - 1,) def is_cellwise_constant(self): "Return whether this expression is spatially constant over each cell." @@ -330,7 +335,7 @@ def ufl_shape(self): g = self._domain.geometric_dimension() t = self._domain.topological_dimension() - return (g, t-1) + return (g, t - 1) def is_cellwise_constant(self): "Return whether this expression is spatially constant over each cell." @@ -357,7 +362,7 @@ @property def ufl_shape(self): t = self._domain.topological_dimension() - return (t, t-1) + return (t, t - 1) def is_cellwise_constant(self): "Return whether this expression is spatially constant over each cell." @@ -573,7 +578,7 @@ def ufl_shape(self): g = self._domain.geometric_dimension() t = self._domain.topological_dimension() - return (t-1, g) + return (t - 1, g) def is_cellwise_constant(self): "Return whether this expression is spatially constant over each cell." @@ -597,7 +602,7 @@ @property def ufl_shape(self): t = self._domain.topological_dimension() - return (t-1, t) + return (t - 1, t) def is_cellwise_constant(self): "Return whether this expression is spatially constant over each cell." diff -Nru ufl-2018.1.0/ufl/index_combination_utils.py ufl-2019.1.0/ufl/index_combination_utils.py --- ufl-2018.1.0/ufl/index_combination_utils.py 2018-06-14 14:20:02.000000000 +0000 +++ ufl-2019.1.0/ufl/index_combination_utils.py 2019-04-17 18:51:30.000000000 +0000 @@ -104,7 +104,7 @@ nrfi = len(rfi) nfi = len(fi) - shape = [None]*nrfi + shape = [None] * nrfi k = 0 pos = 0 newfiid = [] @@ -246,7 +246,7 @@ # Consistency checks if len(set(free_indices)) != len(free_indices): error("Not expecting repeated indices left.") - if len(free_indices) + 2*len(repeated_indices) != an + bn: + if len(free_indices) + 2 * len(repeated_indices) != an + bn: error("Expecting only twice repeated indices.") return tuple(free_indices), tuple(index_dimensions), tuple(repeated_indices), tuple(repeated_index_dimensions) diff -Nru ufl-2018.1.0/ufl/indexsum.py ufl-2019.1.0/ufl/indexsum.py --- ufl-2018.1.0/ufl/indexsum.py 2018-06-14 14:20:02.000000000 +0000 +++ ufl-2019.1.0/ufl/indexsum.py 2019-04-17 18:51:30.000000000 +0000 @@ -46,7 +46,7 @@ if not isinstance(index, MultiIndex): error("Expecting MultiIndex instance, got %s" % ufl_err_str(index)) if len(index) != 1: - error("Expecting a single Index onlym got %d." % len(index)) + error("Expecting a single Index but got %d." % len(index)) # Simplification to zero if isinstance(summand, Zero): @@ -55,8 +55,8 @@ fi = summand.ufl_free_indices fid = summand.ufl_index_dimensions pos = fi.index(j.count()) - fi = fi[:pos] + fi[pos+1:] - fid = fid[:pos] + fid[pos+1:] + fi = fi[:pos] + fi[pos + 1:] + fid = fid[:pos] + fid[pos + 1:] return Zero(sh, fi, fid) return Operator.__new__(cls) @@ -67,8 +67,8 @@ fid = summand.ufl_index_dimensions pos = fi.index(j.count()) self._dimension = fid[pos] - self.ufl_free_indices = fi[:pos] + fi[pos+1:] - self.ufl_index_dimensions = fid[:pos] + fid[pos+1:] + self.ufl_free_indices = fi[:pos] + fi[pos + 1:] + self.ufl_index_dimensions = fid[:pos] + fid[pos + 1:] Operator.__init__(self, (summand, index)) def index(self): diff -Nru ufl-2018.1.0/ufl/__init__.py ufl-2019.1.0/ufl/__init__.py --- ufl-2018.1.0/ufl/__init__.py 2018-06-14 14:20:02.000000000 +0000 +++ ufl-2019.1.0/ufl/__init__.py 2019-04-17 18:51:30.000000000 +0000 @@ -185,6 +185,11 @@ - cosh, sinh, tanh - bessel_J, bessel_Y, bessel_I, bessel_K +* Complex operations:: + + - conj, real, imag + conjugate is an alias for conj + * Discontinuous Galerkin operators:: - v('+'), v('-') @@ -246,7 +251,7 @@ __version__ = pkg_resources.get_distribution("fenics-ufl").version -########## README +# README # Imports here should be what the user sees when doing "from ufl import *", # which means we should _not_ import f.ex. "Grad", but "grad". # This way we expose the language, the operation "grad", but less @@ -271,7 +276,7 @@ CellVolume, CellDiameter, Circumradius, MinCellEdgeLength, MaxCellEdgeLength, FacetArea, MinFacetEdgeLength, MaxFacetEdgeLength, Jacobian, JacobianDeterminant, JacobianInverse - ) +) # Sobolev spaces from ufl.sobolevspace import L2, H1, H2, HDiv, HCurl @@ -284,18 +289,18 @@ FacetElement, InteriorElement # Hook to extend predefined element families -from ufl.finiteelement.elementlist import register_element, show_elements #, ufl_elements +from ufl.finiteelement.elementlist import register_element, show_elements # , ufl_elements # Function spaces from ufl.functionspace import FunctionSpace, MixedFunctionSpace # Arguments from ufl.argument import Argument, TestFunction, TrialFunction, \ - Arguments, TestFunctions, TrialFunctions + Arguments, TestFunctions, TrialFunctions # Coefficients from ufl.coefficient import Coefficient, Coefficients, \ - Constant, VectorConstant, TensorConstant + Constant, VectorConstant, TensorConstant # Split function from ufl.split_functions import split @@ -316,21 +321,22 @@ # Operators from ufl.operators import rank, shape, \ - outer, inner, dot, cross, perp, \ - det, inv, cofac, \ - transpose, tr, diag, diag_vector, \ - dev, skew, sym, \ - sqrt, exp, ln, erf, \ - cos, sin, tan, \ - acos, asin, atan, atan_2, \ - cosh, sinh, tanh, \ - bessel_J, bessel_Y, bessel_I, bessel_K, \ - eq, ne, le, ge, lt, gt, And, Or, Not, \ - conditional, sign, max_value, min_value, Max, Min, \ - variable, diff, \ - Dx, grad, div, curl, rot, nabla_grad, nabla_div, Dn, exterior_derivative, \ - jump, avg, cell_avg, facet_avg, \ - elem_mult, elem_div, elem_pow, elem_op + conj, real, imag, \ + outer, inner, dot, cross, perp, \ + det, inv, cofac, \ + transpose, tr, diag, diag_vector, \ + dev, skew, sym, \ + sqrt, exp, ln, erf, \ + cos, sin, tan, \ + acos, asin, atan, atan_2, \ + cosh, sinh, tanh, \ + bessel_J, bessel_Y, bessel_I, bessel_K, \ + eq, ne, le, ge, lt, gt, And, Or, Not, \ + conditional, sign, max_value, min_value, Max, Min, \ + variable, diff, \ + Dx, grad, div, curl, rot, nabla_grad, nabla_div, Dn, exterior_derivative, \ + jump, avg, cell_avg, facet_avg, \ + elem_mult, elem_div, elem_pow, elem_op # Measure classes from ufl.measure import Measure, register_integral_type, integral_types, custom_integral_types @@ -347,7 +353,7 @@ # Representations of transformed forms from ufl.formoperators import replace, derivative, action, energy_norm, rhs, lhs,\ - system, functional, adjoint, sensitivity_rhs, block_split #, dirichlet_functional + system, functional, adjoint, sensitivity_rhs, block_split # , dirichlet_functional # Predefined convenience objects from ufl.objects import ( @@ -357,7 +363,7 @@ dx, ds, dS, dP, dc, dC, dO, dI, dX, ds_b, ds_t, ds_tb, ds_v, dS_h, dS_v - ) +) # Useful constants from math import e, pi @@ -394,7 +400,7 @@ 'Index', 'indices', 'as_tensor', 'as_vector', 'as_matrix', 'relabel', 'unit_vector', 'unit_vectors', 'unit_matrix', 'unit_matrices', - 'rank', 'shape', + 'rank', 'shape', 'conj', 'real', 'imag', 'outer', 'inner', 'dot', 'cross', 'perp', 'det', 'inv', 'cofac', 'transpose', 'tr', 'diag', 'diag_vector', 'dev', 'skew', 'sym', diff -Nru ufl-2018.1.0/ufl/integral.py ufl-2019.1.0/ufl/integral.py --- ufl-2018.1.0/ufl/integral.py 2018-06-14 14:20:02.000000000 +0000 +++ ufl-2019.1.0/ufl/integral.py 2019-04-17 18:51:30.000000000 +0000 @@ -112,13 +112,13 @@ def __mul__(self, scalar): if not is_python_scalar(scalar): error("Cannot multiply an integral with non-constant values.") - return self.reconstruct(scalar*self._integrand) + return self.reconstruct(scalar * self._integrand) def __rmul__(self, scalar): if not is_scalar_constant_expression(scalar): error("An integral can only be multiplied by a " "globally constant scalar expression.") - return self.reconstruct(scalar*self._integrand) + return self.reconstruct(scalar * self._integrand) def __str__(self): fmt = "{ %s } * %s(%s[%s], %s)" diff -Nru ufl-2018.1.0/ufl/log.py ufl-2019.1.0/ufl/log.py --- ufl-2018.1.0/ufl/log.py 2018-06-14 14:20:02.000000000 +0000 +++ ufl-2019.1.0/ufl/log.py 2019-04-17 18:51:30.000000000 +0000 @@ -174,7 +174,7 @@ def begin(self, *message): "Begin task: write message and increase indentation level." self.info(*message) - self.info("-"*len(self._format_raw(*message))) + self.info("-" * len(self._format_raw(*message))) self.add_indent() def end(self): @@ -231,7 +231,7 @@ def _format(self, *message): "Format message including indentation." - indent = self._prefix + 2*self._indent_level*" " + indent = self._prefix + 2 * self._indent_level * " " return "\n".join([indent + line for line in (message[0] % message[1:]).split("\n")]) def _format_raw(self, *message): diff -Nru ufl-2018.1.0/ufl/mathfunctions.py ufl-2019.1.0/ufl/mathfunctions.py --- ufl-2018.1.0/ufl/mathfunctions.py 2018-06-14 14:20:02.000000000 +0000 +++ ufl-2019.1.0/ufl/mathfunctions.py 2019-04-17 18:51:30.000000000 +0000 @@ -22,11 +22,14 @@ # Modified by Kristian B. Oelgaard, 2011 import math +import cmath +import numbers + from ufl.utils.str import as_native_strings 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, ScalarValue, Zero, FloatValue, IntValue, as_ufl +from ufl.constantvalue import is_true_ufl_scalar, Zero, RealValue, FloatValue, IntValue, ComplexValue, ConstantValue, as_ufl """ TODO: Include additional functions available in (need derivatives as well): @@ -67,7 +70,10 @@ def evaluate(self, x, mapping, component, index_values): a = self.ufl_operands[0].evaluate(x, mapping, component, index_values) try: - res = getattr(math, self._name)(a) + 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 @@ -82,8 +88,13 @@ __slots__ = () def __new__(cls, argument): - if isinstance(argument, (ScalarValue, Zero)): - return FloatValue(math.sqrt(float(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): @@ -95,8 +106,10 @@ __slots__ = () def __new__(cls, argument): - if isinstance(argument, (ScalarValue, Zero)): + 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): @@ -108,8 +121,10 @@ __slots__ = () def __new__(cls, argument): - if isinstance(argument, (ScalarValue, Zero)): + 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): @@ -117,7 +132,10 @@ def evaluate(self, x, mapping, component, index_values): a = self.ufl_operands[0].evaluate(x, mapping, component, index_values) - return math.log(a) + try: + return math.log(a) + except TypeError: + return cmath.log(a) @ufl_type() @@ -125,8 +143,10 @@ __slots__ = () def __new__(cls, argument): - if isinstance(argument, (ScalarValue, Zero)): + 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): @@ -138,8 +158,10 @@ __slots__ = () def __new__(cls, argument): - if isinstance(argument, (ScalarValue, Zero)): + 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): @@ -151,8 +173,10 @@ __slots__ = () def __new__(cls, argument): - if isinstance(argument, (ScalarValue, Zero)): + 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): @@ -164,8 +188,10 @@ __slots__ = () def __new__(cls, argument): - if isinstance(argument, (ScalarValue, Zero)): + 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): @@ -177,8 +203,10 @@ __slots__ = () def __new__(cls, argument): - if isinstance(argument, (ScalarValue, Zero)): + 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): @@ -190,8 +218,10 @@ __slots__ = () def __new__(cls, argument): - if isinstance(argument, (ScalarValue, Zero)): + 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): @@ -203,8 +233,10 @@ __slots__ = () def __new__(cls, argument): - if isinstance(argument, (ScalarValue, Zero)): + 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): @@ -216,8 +248,10 @@ __slots__ = () def __new__(cls, argument): - if isinstance(argument, (ScalarValue, Zero)): + 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): @@ -229,8 +263,10 @@ __slots__ = () def __new__(cls, argument): - if isinstance(argument, (ScalarValue, Zero)): + 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): @@ -242,12 +278,16 @@ __slots__ = () def __new__(cls, arg1, arg2): - if isinstance(arg1, (ScalarValue, Zero)) and isinstance(arg2, (ScalarValue, Zero)): + 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): @@ -258,6 +298,8 @@ 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 @@ -282,10 +324,14 @@ __slots__ = () def __new__(cls, argument): - if isinstance(argument, (ScalarValue, Zero)): + 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): diff -Nru ufl-2018.1.0/ufl/measure.py ufl-2019.1.0/ufl/measure.py --- ufl-2018.1.0/ufl/measure.py 2018-06-14 14:20:02.000000000 +0000 +++ ufl-2019.1.0/ufl/measure.py 2019-04-17 18:51:30.000000000 +0000 @@ -423,7 +423,7 @@ # Integral.__add__: subdomain_id = self.subdomain_id() if isinstance(subdomain_id, tuple): - return sum(integrand*self.reconstruct(subdomain_id=d) for d in subdomain_id) + return sum(integrand * self.reconstruct(subdomain_id=d) for d in subdomain_id) # Check that we have an integer subdomain or a string # ("everywhere" or "otherwise", any more?) @@ -469,7 +469,7 @@ self._measures = measures def __rmul__(self, other): - integrals = [other*m for m in self._measures] + integrals = [other * m for m in self._measures] return sum(integrals) def __add__(self, other): diff -Nru ufl-2018.1.0/ufl/operators.py ufl-2019.1.0/ufl/operators.py --- ufl-2018.1.0/ufl/operators.py 2018-06-14 14:20:02.000000000 +0000 +++ ufl-2019.1.0/ufl/operators.py 2019-04-17 18:51:30.000000000 +0000 @@ -28,7 +28,7 @@ from ufl.log import error, warning from ufl.form import Form -from ufl.constantvalue import Zero, ScalarValue, as_ufl +from ufl.constantvalue import Zero, RealValue, ComplexValue, as_ufl from ufl.differentiation import VariableDerivative, Grad, Div, Curl, NablaGrad, NablaDiv from ufl.tensoralgebra import Transposed, Inner, Outer, Dot, Cross, \ Determinant, Inverse, Cofactor, Trace, Deviatoric, Skew, Sym @@ -37,6 +37,7 @@ from ufl.tensors import as_tensor, as_matrix, as_vector, ListTensor from ufl.conditional import EQ, NE, \ AndCondition, OrCondition, NotCondition, Conditional, MaxValue, MinValue +from ufl.algebra import Conj, Real, Imag from ufl.mathfunctions import Sqrt, Exp, Ln, Erf,\ Cos, Sin, Tan, Cosh, Sinh, Tanh, Acos, Asin, Atan, Atan2,\ BesselJ, BesselY, BesselI, BesselK @@ -47,9 +48,9 @@ from ufl.checks import is_cellwise_constant from ufl.domain import extract_domains - # --- Basic operators --- + def rank(f): "UFL operator: The rank of *f*." f = as_ufl(f) @@ -62,6 +63,30 @@ return f.ufl_shape +# --- Complex operators --- + +def conj(f): + "UFL operator: The complex conjugate of *f*" + f = as_ufl(f) + return Conj(f) + + +# Alias because both conj and conjugate are in numpy and we wish to be consistent. +conjugate = conj + + +def real(f): + "UFL operator: The real part of *f*" + f = as_ufl(f) + return Real(f) + + +def imag(f): + "UFL operator: The imaginary part of *f*" + f = as_ufl(f) + return Imag(f) + + # --- Elementwise tensor operators --- def elem_op_items(op_ind, indices, *args): @@ -72,7 +97,7 @@ def extind(ii): return indices + (ii,) - if len(sh) == len(indices)+1: + if len(sh) == len(indices) + 1: return [op_ind(extind(i), *args) for i in range(n)] else: return [elem_op_items(op_ind, extind(i), *args) for i in range(n)] @@ -119,7 +144,7 @@ def outer(*operands): - "UFL operator: Take the outer product of two or more operands." + "UFL operator: Take the outer product of two or more operands. The complex conjugate of the first argument is taken." n = len(operands) if n == 1: return operands[0] @@ -131,16 +156,16 @@ a = as_ufl(a) b = as_ufl(b) if a.ufl_shape == () and b.ufl_shape == (): - return a*b + return Conj(a) * b return Outer(a, b) def inner(a, b): - "UFL operator: Take the inner product of *a* and *b*." + "UFL operator: Take the inner product of *a* and *b*. The complex conjugate of the second argument is taken." a = as_ufl(a) b = as_ufl(b) if a.ufl_shape == () and b.ufl_shape == (): - return a*b + return a * Conj(b) return Inner(a, b) @@ -150,15 +175,15 @@ "UFL operator: Take the partial inner product of a and b." ar, br = len(a.ufl_shape), len(b.ufl_shape) n = min(ar, br) - return contraction(a, list(range(n-ar, n-ar+n)), b, list(range(n))) + return contraction(a, list(range(n - ar, n - ar + n)), b, list(range(n))) def dot(a, b): - "UFL operator: Take the dot product of *a* and *b*." + "UFL operator: Take the dot product of *a* and *b*. The complex conjugate of the second argument is taken." a = as_ufl(a) b = as_ufl(b) if a.ufl_shape == () and b.ufl_shape == (): - return a*b + return a * b return Dot(a, b) @@ -172,7 +197,7 @@ aii = indices(len(a.ufl_shape)) bii = indices(len(b.ufl_shape)) cii = indices(len(ai)) - shape = [None]*len(ai) + shape = [None] * len(ai) for i, j in enumerate(ai): aii[j] = cii[i] shape[i] = ash[j] @@ -180,7 +205,7 @@ bii[j] = cii[i] if shape[i] != bsh[j]: error("Shape mismatch in contraction.") - s = a[aii]*b[bii] + s = a[aii] * b[bii] cii = set(cii) ii = tuple(i for i in (aii + bii) if i not in cii) return as_tensor(s, ii) @@ -251,7 +276,7 @@ # Build matrix row by row rows = [] for i in range(n): - row = [0]*n + row = [0] * n row[i] = A[i] if r == 1 else A[i, i] rows.append(row) return as_matrix(rows) @@ -432,7 +457,7 @@ return v('+') - v('-') r = len(v.ufl_shape) if r == 0: - return v('+')*n('+') + v('-')*n('-') + return v('+') * n('+') + v('-') * n('-') else: return dot(v('+'), n('+')) + dot(v('-'), n('-')) else: @@ -449,7 +474,7 @@ def avg(v): "UFL operator: Take the average of *v* across a facet." v = as_ufl(v) - return 0.5*(v('+') + v('-')) + return 0.5 * (v('+') + v('-')) def cell_avg(f): @@ -568,8 +593,10 @@ def _mathfunction(f, cls): f = as_ufl(f) r = cls(f) - if isinstance(r, (ScalarValue, Zero, int, float)): + if isinstance(r, (RealValue, Zero, int, float)): return float(r) + if isinstance(r, (ComplexValue, complex)): + return complex(r) return r @@ -637,9 +664,13 @@ "UFL operator: Take the inverse tangent with two the arguments *f1* and *f2*." f1 = as_ufl(f1) f2 = as_ufl(f2) + if isinstance(f1, (ComplexValue, complex)) or isinstance(f2, (ComplexValue, complex)): + raise TypeError('atan_2 is incompatible with complex numbers.') r = Atan2(f1, f2) - if isinstance(r, (ScalarValue, Zero, int, float)): + if isinstance(r, (RealValue, Zero, int, float)): return float(r) + if isinstance(r, (ComplexValue, complex)): + return complex(r) return r diff -Nru ufl-2018.1.0/ufl/sobolevspace.py ufl-2019.1.0/ufl/sobolevspace.py --- ufl-2018.1.0/ufl/sobolevspace.py 2018-06-14 14:20:02.000000000 +0000 +++ ufl-2019.1.0/ufl/sobolevspace.py 2019-04-17 18:51:30.000000000 +0000 @@ -25,9 +25,10 @@ # Modified by Lizao Li 2015 # Modified by Thomas Gibson 2017 -from ufl.utils.str import as_native_str from functools import total_ordering +from ufl.utils.str import as_native_str + @total_ordering class SobolevSpace(object): @@ -47,30 +48,27 @@ p = frozenset(parents or []) # Ensure that the inclusion operations are transitive. self.parents = p.union(*[p_.parents for p_ in p]) - self._order = {"L2": 0, - "H1": 1, - "H2": 2, - # Order for the elements below is taken from - # its parent Sobolev space - "HDiv": 0, - "HCurl": 0, - "HEin": 0, - "HDivDiv": 0, - "DirectionalH": 0}[self.name] + self._order = { + "L2": 0, + "H1": 1, + "H2": 2, + # Order for the elements below is taken from + # its parent Sobolev space + "HDiv": 0, + "HCurl": 0, + "HEin": 0, + "HDivDiv": 0, + "DirectionalH": 0 + }[self.name] def __str__(self): return self.name def __repr__(self): - r = "SobolevSpace(%s, %s)" % (repr(self.name), repr(list(self.parents))) + r = "SobolevSpace(%s, %s)" % (repr(self.name), repr( + list(self.parents))) return as_native_str(r) - def _repr_latex_(self): - if len(self.name) == 2: - return "$%s^%s$" % tuple(self.name) - else: - return "$%s(%s)$" % (self.name[0], self.name[1:].lower()) - def __eq__(self, other): return isinstance(other, SobolevSpace) and self.name == other.name @@ -110,7 +108,9 @@ elif self.name == "HCurl": from ufl.finiteelement import HCurlElement return HCurlElement(element) - raise NotImplementedError("SobolevSpace has no call operator (only the specific HDiv and HCurl instances).") + raise NotImplementedError( + "SobolevSpace has no call operator (only the specific HDiv and HCurl instances)." + ) @total_ordering @@ -127,12 +127,11 @@ the position denotes in what spatial variable the smoothness requirement is enforced. """ - assert all(isinstance(x, int) for x in orders), ( - "Order must be an integer." - ) - assert all(x < 3 for x in orders), ( - "Not implemented for orders greater than 2" - ) + assert all( + isinstance(x, int) for x in orders), ("Order must be an integer.") + assert all( + x < 3 + for x in orders), ("Not implemented for orders greater than 2") name = "DirectionalH" parents = [L2] super(DirectionalSobolevSpace, self).__init__(name, parents) @@ -145,9 +144,7 @@ """ if spatial_index not in range(len(self._orders)): raise IndexError("Spatial index out of range.") - spaces = {0: L2, - 1: H1, - 2: H2} + spaces = {0: L2, 1: H1, 2: H2} return spaces[self._orders[spatial_index]] def __contains__(self, other): @@ -181,18 +178,14 @@ elif other.name in ["HDivDiv", "HEin"]: # Don't know how these spaces compare return NotImplementedError( - "Don't know how to compare with %s" % other.name - ) + "Don't know how to compare with %s" % other.name) else: - return any(self._orders[i] > other._order - for i in self._spatial_indices) + return any( + self._orders[i] > other._order for i in self._spatial_indices) def __str__(self): return self.name + "(%s)" % ", ".join(map(str, self._orders)) - def _repr_latex_(self): - return "H(%s)" % ", ".join(map(str, self._orders)) - L2 = SobolevSpace("L2") HDiv = SobolevSpace("HDiv", [L2]) diff -Nru ufl-2018.1.0/ufl/sorting.py ufl-2019.1.0/ufl/sorting.py --- ufl-2018.1.0/ufl/sorting.py 2018-06-14 14:20:02.000000000 +0000 +++ ufl-2019.1.0/ufl/sorting.py 2019-04-17 18:51:30.000000000 +0000 @@ -104,7 +104,7 @@ # Hack up a MultiFunction-like type dispatch for terminal comparisons -_terminal_cmps = [_cmp_terminal_by_repr]*Expr._ufl_num_typecodes_ +_terminal_cmps = [_cmp_terminal_by_repr] * Expr._ufl_num_typecodes_ _terminal_cmps[MultiIndex._ufl_typecode_] = _cmp_multi_index _terminal_cmps[Argument._ufl_typecode_] = _cmp_argument _terminal_cmps[Coefficient._ufl_typecode_] = _cmp_coefficient diff -Nru ufl-2018.1.0/ufl/split_functions.py ufl-2019.1.0/ufl/split_functions.py --- ufl-2018.1.0/ufl/split_functions.py 2018-06-14 14:20:02.000000000 +0000 +++ ufl-2019.1.0/ufl/split_functions.py 2019-04-17 18:51:30.000000000 +0000 @@ -108,7 +108,7 @@ elif rank <= 1: subv = as_vector(components) elif rank == 2: - subv = as_matrix([components[i*shape[1]: (i+1)*shape[1]] + subv = as_matrix([components[i * shape[1]: (i + 1) * shape[1]] for i in range(shape[0])]) else: error("Don't know how to split functions with sub functions of rank %d." % rank) diff -Nru ufl-2018.1.0/ufl/tensoralgebra.py ufl-2019.1.0/ufl/tensoralgebra.py --- ufl-2018.1.0/ufl/tensoralgebra.py 2018-06-14 14:20:02.000000000 +0000 +++ ufl-2019.1.0/ufl/tensoralgebra.py 2019-04-17 18:51:30.000000000 +0000 @@ -23,7 +23,7 @@ from ufl.core.expr import ufl_err_str from ufl.core.ufl_type import ufl_type from ufl.constantvalue import Zero -from ufl.algebra import Operator +from ufl.algebra import Operator, Conj from ufl.precedence import parstr from ufl.sorting import sorted_expr from ufl.index_combination_utils import merge_nonoverlapping_indices @@ -131,7 +131,7 @@ fi, fid = merge_nonoverlapping_indices(a, b) return Zero(ash + bsh, fi, fid) if ash == () or bsh == (): - return a * b + return Conj(a) * b return CompoundTensorOperator.__new__(cls) def __init__(self, a, b): @@ -165,14 +165,16 @@ fi, fid = merge_nonoverlapping_indices(a, b) return Zero((), fi, fid) elif ash == (): - return a*b + return a * Conj(b) + # sort operands for unique representation, + # must be independent of various counts etc. + # as explained in cmp_expr + if (a, b) != tuple(sorted_expr((a, b))): + return Conj(Inner(b, a)) return CompoundTensorOperator.__new__(cls) def __init__(self, a, b): - # sort operands for unique representation, must be independent - # of various counts etc. as explained in cmp_expr - a, b = sorted_expr((a, b)) CompoundTensorOperator.__init__(self, (a, b)) diff -Nru ufl-2018.1.0/ufl/tensors.py ufl-2019.1.0/ufl/tensors.py --- ufl-2018.1.0/ufl/tensors.py 2018-06-14 14:20:02.000000000 +0000 +++ ufl-2019.1.0/ufl/tensors.py 2019-04-17 18:51:30.000000000 +0000 @@ -104,12 +104,12 @@ def __str__(self): def substring(expressions, indent): - ind = " "*indent + ind = " " * indent if any(isinstance(e, ListTensor) for e in expressions): substrings = [] for e in expressions: if isinstance(e, ListTensor): - substrings.append(substring(e.ufl_operands, indent+2)) + substrings.append(substring(e.ufl_operands, indent + 2)) else: substrings.append(str(e)) s = (",\n" + ind).join(substrings) @@ -432,5 +432,5 @@ else: for s, sub in enumerate(subs): for c, v in unwrap_list_tensor(sub): - components.append(((s,)+c, v)) + components.append(((s,) + c, v)) return components diff -Nru ufl-2018.1.0/ufl/utils/derivativetuples.py ufl-2019.1.0/ufl/utils/derivativetuples.py --- ufl-2018.1.0/ufl/utils/derivativetuples.py 2018-06-14 14:20:02.000000000 +0000 +++ ufl-2019.1.0/ufl/utils/derivativetuples.py 2019-04-17 18:51:30.000000000 +0000 @@ -30,7 +30,7 @@ """ derivatives = [] # = 1 for i, d in enumerate(derivative_counts): - derivatives.extend((i,)*d) # *= d/dx_i^d + derivatives.extend((i,) * d) # *= d/dx_i^d return tuple(derivatives) @@ -41,7 +41,7 @@ in counting form as (0, 2, 1) meaning (dx^0, dy^2, dz^1) and in listing form as (1, 1, 2) meaning (dy, dy, dz). """ - derivative_counts = [0]*gdim + derivative_counts = [0] * gdim for d in derivatives: derivative_counts[d] += 1 return tuple(derivative_counts) @@ -71,7 +71,7 @@ """ # Create list of derivatives (note that we have d^n derivatives) - deriv_tuples = [d for d in itertools.product(*(n*[range(0, gdim)]))] + deriv_tuples = [d for d in itertools.product(*(n * [range(0, gdim)]))] # Translate from list of derivative tuples to list of tuples # expressing the number of derivatives in each dimension... diff -Nru ufl-2018.1.0/ufl/utils/formatting.py ufl-2019.1.0/ufl/utils/formatting.py --- ufl-2018.1.0/ufl/utils/formatting.py 2018-06-14 14:20:02.000000000 +0000 +++ ufl-2019.1.0/ufl/utils/formatting.py 2019-04-17 18:51:30.000000000 +0000 @@ -70,13 +70,13 @@ value = "'%s'" % value else: value = str(value) - s += key + ":" + " "*(keylen - len(key) + 1) + s += key + ":" + " " * (keylen - len(key) + 1) space = "" while len(value) > 0: end = min(len(value), colsize - keylen) s += space + value[:end] + "\n" value = value[end:] - space = " "*(keylen + 2) + space = " " * (keylen + 2) return s diff -Nru ufl-2018.1.0/ufl/utils/indexflattening.py ufl-2019.1.0/ufl/utils/indexflattening.py --- ufl-2018.1.0/ufl/utils/indexflattening.py 2018-06-14 14:20:02.000000000 +0000 +++ ufl-2019.1.0/ufl/utils/indexflattening.py 2019-04-17 18:51:30.000000000 +0000 @@ -24,10 +24,10 @@ n = len(sh) if not n: return () - strides = [None]*n - strides[n-1] = 1 - for i in range(n-1, 0, -1): - strides[i-1] = strides[i]*sh[i] + strides = [None] * n + strides[n - 1] = 1 + for i in range(n - 1, 0, -1): + strides[i - 1] = strides[i] * sh[i] return tuple(strides) diff -Nru ufl-2018.1.0/ufl/utils/stacks.py ufl-2019.1.0/ufl/utils/stacks.py --- ufl-2018.1.0/ufl/utils/stacks.py 2018-06-14 14:20:02.000000000 +0000 +++ ufl-2019.1.0/ufl/utils/stacks.py 2019-04-17 18:51:30.000000000 +0000 @@ -21,6 +21,7 @@ class Stack(list): "A stack datastructure." + def __init__(self, *args): list.__init__(self, *args) @@ -33,6 +34,7 @@ class StackDict(dict): "A dict that can be changed incrementally with 'd.push(k,v)' and have changes rolled back with 'k,v = d.pop()'." + def __init__(self, *args, **kwargs): dict.__init__(self, *args, **kwargs) self._l = []