From a49736c7e4372b9bac61d728c223c8b6784b25a3 Mon Sep 17 00:00:00 2001 From: Matthew Scroggs Date: Wed, 9 Aug 2023 16:13:24 +0100 Subject: [PATCH 01/11] avoid error if Label not in renumbering (#187) --- ufl/variable.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/ufl/variable.py b/ufl/variable.py index b5a29dc26..935c2f8c8 100644 --- a/ufl/variable.py +++ b/ufl/variable.py @@ -51,6 +51,8 @@ def ufl_domains(self): return () def _ufl_signature_data_(self, renumbering): + if self not in renumbering: + return ("Label", self._count) return ("Label", renumbering[self]) From 9e120c6ace16f13c7b8bf9398d1dd1a71d695f7f Mon Sep 17 00:00:00 2001 From: Matthew Scroggs Date: Thu, 10 Aug 2023 07:23:43 +0100 Subject: [PATCH 02/11] Added 4d element types (#188) * [feature][ElementTypes] Added 4d types pentatope (4d simplex) and tesseract (4d cube). (#95) Co-authored-by: Jack S. Hale Co-authored-by: David A. Ham Co-authored-by: Matthew Scroggs * fix merge, add test --------- Co-authored-by: dr-robertk Co-authored-by: Jack S. Hale Co-authored-by: David A. Ham --- test/test_cell.py | 23 ++++++++++++++++++++++- ufl/__init__.py | 6 ++++-- ufl/cell.py | 8 ++++++++ ufl/finiteelement/elementlist.py | 7 ++++--- ufl/objects.py | 2 ++ 5 files changed, 40 insertions(+), 6 deletions(-) diff --git a/test/test_cell.py b/test/test_cell.py index 21a5240fa..8ec754622 100644 --- a/test/test_cell.py +++ b/test/test_cell.py @@ -37,6 +37,20 @@ def test_hexahedron(): assert cell.num_faces() == 6 +def test_pentatope(): + cell = ufl.pentatope + assert cell.num_vertices() == 5 + assert cell.num_edges() == 10 + assert cell.num_faces() == 10 + + +def test_tesseract(): + cell = ufl.tesseract + assert cell.num_vertices() == 16 + assert cell.num_edges() == 32 + assert cell.num_faces() == 24 + + @pytest.mark.parametrize("cell", [ufl.interval]) def test_cells_1d(cell): @@ -53,12 +67,19 @@ def test_cells_2d(cell): @pytest.mark.parametrize("cell", [ufl.tetrahedron, ufl.hexahedron, ufl.prism, ufl.pyramid]) -def test_cells_2d(cell): +def test_cells_3d(cell): assert cell.num_facets() == cell.num_faces() assert cell.num_ridges() == cell.num_edges() assert cell.num_peaks() == cell.num_vertices() +@pytest.mark.parametrize("cell", [ufl.tesseract, ufl.pentatope]) +def test_cells_4d(cell): + assert cell.num_facets() == cell.num_sub_entities(3) + assert cell.num_ridges() == cell.num_faces() + assert cell.num_peaks() == cell.num_edges() + + def test_tensorproductcell(): orig = ufl.TensorProductCell(ufl.interval, ufl.interval) cell = orig.reconstruct() diff --git a/ufl/__init__.py b/ufl/__init__.py index 7e7ca3cdb..7cb9cf0e0 100644 --- a/ufl/__init__.py +++ b/ufl/__init__.py @@ -53,6 +53,8 @@ - hexahedron - prism - pyramid + - pentatope + - tesseract * Domains:: @@ -339,7 +341,7 @@ # Predefined convenience objects from ufl.objects import ( - vertex, interval, triangle, tetrahedron, + vertex, interval, triangle, tetrahedron, pentatope, tesseract, quadrilateral, hexahedron, prism, pyramid, facet, i, j, k, l, p, q, r, s, dx, ds, dS, dP, @@ -400,7 +402,7 @@ 'dc', 'dC', 'dO', 'dI', 'dX', 'ds_b', 'ds_t', 'ds_tb', 'ds_v', 'dS_h', 'dS_v', 'vertex', 'interval', 'triangle', 'tetrahedron', - 'prism', 'pyramid', + 'prism', 'pyramid', 'pentatope', 'tesseract', 'quadrilateral', 'hexahedron', 'facet', 'i', 'j', 'k', 'l', 'p', 'q', 'r', 's', 'e', 'pi', diff --git a/ufl/cell.py b/ufl/cell.py index 9ef8e9975..2d50f6a77 100644 --- a/ufl/cell.py +++ b/ufl/cell.py @@ -196,6 +196,10 @@ def peak_types(self) -> typing.Tuple[AbstractCell, ...]: ("triangle", "quadrilateral", "quadrilateral", "quadrilateral", "triangle"), ("prism", )], "pyramid": [tuple("vertex" for i in range(5)), tuple("interval" for i in range(8)), ("quadrilateral", "triangle", "triangle", "triangle", "triangle"), ("pyramid", )], + "pentatope": [tuple("vertex" for i in range(5)), tuple("interval" for i in range(10)), + tuple("triangle" for i in range(10)), tuple("tetrahedron" for i in range(5)), ("pentatope", )], + "tesseract": [tuple("vertex" for i in range(16)), tuple("interval" for i in range(32)), + tuple("quadrilateral" for i in range(24)), tuple("hexahedron" for i in range(8)), ("tesseract", )], } @@ -422,6 +426,8 @@ def simplex(topological_dimension: int, geometric_dimension: typing.Optional[int return Cell("triangle", geometric_dimension) if topological_dimension == 3: return Cell("tetrahedron", geometric_dimension) + if topological_dimension == 4: + return Cell("pentatope", geometric_dimension) raise ValueError(f"Unsupported topological dimension for simplex: {topological_dimension}") @@ -435,6 +441,8 @@ def hypercube(topological_dimension, geometric_dimension=None): return Cell("quadrilateral", geometric_dimension) if topological_dimension == 3: return Cell("hexahedron", geometric_dimension) + if topological_dimension == 4: + return Cell("tesseract", geometric_dimension) raise ValueError(f"Unsupported topological dimension for hypercube: {topological_dimension}") diff --git a/ufl/finiteelement/elementlist.py b/ufl/finiteelement/elementlist.py index da3855641..6ea6a6fa2 100644 --- a/ufl/finiteelement/elementlist.py +++ b/ufl/finiteelement/elementlist.py @@ -12,6 +12,7 @@ # Modified by Marie E. Rognes , 2010 # Modified by Lizao Li , 2015, 2016 # Modified by Massimiliano Leoni, 2016 +# Modified by Robert Kloefkorn, 2022 import warnings from numpy import asarray @@ -82,12 +83,12 @@ def show_elements(): # the future, add mapping name as another element property. # Cell groups -simplices = ("interval", "triangle", "tetrahedron") -cubes = ("interval", "quadrilateral", "hexahedron") +simplices = ("interval", "triangle", "tetrahedron", "pentatope") +cubes = ("interval", "quadrilateral", "hexahedron", "tesseract") any_cell = (None, "vertex", "interval", "triangle", "tetrahedron", "prism", - "pyramid", "quadrilateral", "hexahedron") + "pyramid", "quadrilateral", "hexahedron", "pentatope", "tesseract") # Elements in the periodic table # TODO: Register these as aliases of # periodic table element description instead of the other way around diff --git a/ufl/objects.py b/ufl/objects.py index ab414091d..84f3f360c 100644 --- a/ufl/objects.py +++ b/ufl/objects.py @@ -37,6 +37,8 @@ pyramid = Cell("pyramid", 3) quadrilateral = Cell("quadrilateral", 2) hexahedron = Cell("hexahedron", 3) +tesseract = Cell("tesseract", 4) +pentatope = Cell("pentatope", 4) # Facet is just a dummy declaration for RestrictedElement facet = "facet" From 5d6df80985dc52f009fe91c3b46991a263fdcef5 Mon Sep 17 00:00:00 2001 From: mrambausek Date: Tue, 12 Sep 2023 11:02:43 +0200 Subject: [PATCH 03/11] Fix issue #36 (#37) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit * removed probably unneeded troublesome lines in VectorElement The removed code overwrites some properties for whatever reason, in particular, quad_scheme. In TensorElement these two lines were not there at all. * added back FiniteElementBase constructor but (re)use sub-element data properly Note: the call to FiniteElementBase.__init__ overrides some properties set by the MixedElement constructor. Essentially the same seems to be done in TensorElement, but more explicitly. --------- Co-authored-by: rambausek Co-authored-by: Jørgen Schartum Dokken Co-authored-by: Matthew Scroggs --- ufl/finiteelement/mixedelement.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/ufl/finiteelement/mixedelement.py b/ufl/finiteelement/mixedelement.py index b67fc9f48..f49662122 100644 --- a/ufl/finiteelement/mixedelement.py +++ b/ufl/finiteelement/mixedelement.py @@ -311,8 +311,10 @@ def __init__(self, family, cell=None, degree=None, dim=None, # Initialize element data MixedElement.__init__(self, sub_elements, value_shape=value_shape, reference_value_shape=reference_value_shape) - FiniteElementBase.__init__(self, sub_element.family(), cell, sub_element.degree(), quad_scheme, - value_shape, reference_value_shape) + + FiniteElementBase.__init__(self, sub_element.family(), sub_element.cell(), sub_element.degree(), + sub_element.quadrature_scheme(), value_shape, reference_value_shape) + self._sub_element = sub_element if variant is None: From 551f6735fa238a6edc70a4de48178872190e6074 Mon Sep 17 00:00:00 2001 From: Yang Zongze Date: Tue, 12 Sep 2023 17:03:57 +0800 Subject: [PATCH 04/11] Add function determinant_expr_nxn (#101) * Add function determinant_expr_nxn * Change the order of the assert expression --------- Co-authored-by: Matthew Scroggs --- test/test_apply_algebra_lowering.py | 2 +- ufl/compound_expressions.py | 12 ++++++++++-- 2 files changed, 11 insertions(+), 3 deletions(-) diff --git a/test/test_apply_algebra_lowering.py b/test/test_apply_algebra_lowering.py index a9a097b15..a66332313 100755 --- a/test/test_apply_algebra_lowering.py +++ b/test/test_apply_algebra_lowering.py @@ -55,7 +55,7 @@ def test_determinant2(A2): def test_determinant3(A3): assert determinant_expr(A3) == (A3[0, 0]*(A3[1, 1]*A3[2, 2] - A3[1, 2]*A3[2, 1]) - + A3[0, 1]*(A3[1, 2]*A3[2, 0] - A3[1, 0]*A3[2, 2]) + + (A3[1, 0]*A3[2, 2] - A3[1, 2]*A3[2, 0])*(-A3[0, 1]) + A3[0, 2]*(A3[1, 0]*A3[2, 1] - A3[1, 1]*A3[2, 0])) diff --git a/ufl/compound_expressions.py b/ufl/compound_expressions.py index a25307298..af3e05fd9 100644 --- a/ufl/compound_expressions.py +++ b/ufl/compound_expressions.py @@ -93,6 +93,8 @@ def determinant_expr(A): return determinant_expr_2x2(A) elif sh[0] == 3: return determinant_expr_3x3(A) + else: + return determinant_expr_nxn(A) else: return pseudo_determinant_expr(A) @@ -116,6 +118,12 @@ def determinant_expr_3x3(A): return codeterminant_expr_nxn(A, [0, 1, 2], [0, 1, 2]) +def determinant_expr_nxn(A): + nrow, ncol = A.ufl_shape + assert nrow == ncol + return codeterminant_expr_nxn(A, list(range(nrow)), list(range(ncol))) + + def codeterminant_expr_nxn(A, rows, cols): if len(rows) == 2: return _det_2x2(A, rows[0], rows[1], cols[0], cols[1]) @@ -123,8 +131,8 @@ def codeterminant_expr_nxn(A, rows, cols): r = rows[0] subrows = rows[1:] for i, c in enumerate(cols): - subcols = cols[i + 1:] + cols[:i] - codet += A[r, c] * codeterminant_expr_nxn(A, subrows, subcols) + subcols = cols[:i] + cols[i + 1:] + codet += (-1)**i * A[r, c] * codeterminant_expr_nxn(A, subrows, subcols) return codet From 37ae3531c8e411fef275cfb74712b284a231fc8a Mon Sep 17 00:00:00 2001 From: Cian Wilson Date: Tue, 12 Sep 2023 05:04:23 -0400 Subject: [PATCH 05/11] Possible fix for taking scalar derivatives of vector coefficients using coefficient derivative maps (#124) * A possible fix for taking the derivative of vector coefficients w.r.t. a scalar. Also adding a test. * A slightly neater version that handles scalar derivatives of vectors. --------- Co-authored-by: Matthew Scroggs --- test/test_derivative.py | 24 ++++++++++++++++++++++++ ufl/algorithms/apply_derivatives.py | 6 +++--- 2 files changed, 27 insertions(+), 3 deletions(-) diff --git a/test/test_derivative.py b/test/test_derivative.py index 0259f066e..2a023b9ae 100755 --- a/test/test_derivative.py +++ b/test/test_derivative.py @@ -483,6 +483,30 @@ def test_coefficient_derivatives(self): self.assertEqual(replace(actual, fd.function_replace_map), expected) +def test_vector_coefficient_scalar_derivatives(self): + V = FiniteElement("Lagrange", triangle, 1) + VV = VectorElement("Lagrange", triangle, 1) + + dv = TestFunction(V) + + df = Coefficient(VV, count=0) + g = Coefficient(VV, count=1) + f = Coefficient(VV, count=2) + u = Coefficient(V, count=3) + cd = {f: df} + + integrand = inner(f, g) + + i0, i1, i2, i3, i4 = [Index(count=c) for c in range(5)] + expected = as_tensor(df[i1]*dv, (i1,))[i0]*g[i0] + + F = integrand*dx + J = derivative(F, u, dv, cd) + fd = compute_form_data(J) + actual = fd.preprocessed_form.integrals()[0].integrand() + assert (actual*dx).signature() == (expected*dx).signature() + + def test_vector_coefficient_derivatives(self): V = VectorElement("Lagrange", triangle, 1) VV = TensorElement("Lagrange", triangle, 1) diff --git a/ufl/algorithms/apply_derivatives.py b/ufl/algorithms/apply_derivatives.py index 01c79b8cd..680969426 100644 --- a/ufl/algorithms/apply_derivatives.py +++ b/ufl/algorithms/apply_derivatives.py @@ -835,9 +835,9 @@ def coefficient(self, o): dosum = Zero(o.ufl_shape) for do, v in zip(dos, self._v): so, oi = as_scalar(do) - rv = len(v.ufl_shape) - oi1 = oi[:-rv] - oi2 = oi[-rv:] + rv = len(oi) - len(v.ufl_shape) + oi1 = oi[:rv] + oi2 = oi[rv:] prod = so * v[oi2] if oi1: dosum += as_tensor(prod, oi1) From f08d7c3b738abdbf9eb93f0d57ed1d6747a65fde Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Tue, 12 Sep 2023 10:04:37 +0100 Subject: [PATCH 06/11] Return correct Sobolev space for RestrictedElement (#128) * return correct Sobolev space for FacetElement * InteriorElement in L2 --------- Co-authored-by: Matthew Scroggs --- ufl/finiteelement/restrictedelement.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/ufl/finiteelement/restrictedelement.py b/ufl/finiteelement/restrictedelement.py index 82cd581f4..c1d10b65c 100644 --- a/ufl/finiteelement/restrictedelement.py +++ b/ufl/finiteelement/restrictedelement.py @@ -39,7 +39,10 @@ def __repr__(self): return f"RestrictedElement({repr(self._element)}, {repr(self._restriction_domain)})" def sobolev_space(self): - return L2 + if self._restriction_domain == "interior": + return L2 + else: + return self._element.sobolev_space() def is_cellwise_constant(self): """Return whether the basis functions of this element is spatially From 5803fcb6dab627a462dd378bcc17f7f058aefd3e Mon Sep 17 00:00:00 2001 From: tommbendall Date: Tue, 12 Sep 2023 10:05:01 +0100 Subject: [PATCH 07/11] Implements the perp operator as a class (#184) * Change perp to being a named operator * PR #184: fix behaviour when perp is Zero and expand testing --------- Co-authored-by: Matthew Scroggs --- test/test_tensoralgebra.py | 21 ++++++++++++++++++ ufl/algorithms/apply_algebra_lowering.py | 3 +++ ufl/operators.py | 4 ++-- ufl/tensoralgebra.py | 28 ++++++++++++++++++++++++ 4 files changed, 54 insertions(+), 2 deletions(-) diff --git a/test/test_tensoralgebra.py b/test/test_tensoralgebra.py index 0217311dc..e8831a1d8 100755 --- a/test/test_tensoralgebra.py +++ b/test/test_tensoralgebra.py @@ -105,6 +105,27 @@ def test_cross(self): self.assertEqualValues(C, D) +def test_perp(self): + # Test perp is generally doing the correct thing + u = as_vector([3, 1]) + v = perp(u) + w = as_vector([-1, 3]) + self.assertEqualValues(v, w) + + # Test that a perp does the correct thing to Zero + u = zero(2) + v = perp(u) + self.assertEqualValues(u, v) + + # Test that perp throws an error if the wrong thing is provided + u = as_vector([3, 1, -1]) # 3D vector instead of 2D + with pytest.raises(ValueError): + v = perp(u) + u = as_matrix([[1, 3], [0, 4]]) # Matrix instead of vector + with pytest.raises(ValueError): + v = perp(u) + + def xtest_dev(self, A): C = dev(A) D = 0*C # FIXME: Add expected value here diff --git a/ufl/algorithms/apply_algebra_lowering.py b/ufl/algorithms/apply_algebra_lowering.py index 1aec94467..3938ed01f 100644 --- a/ufl/algorithms/apply_algebra_lowering.py +++ b/ufl/algorithms/apply_algebra_lowering.py @@ -55,6 +55,9 @@ def c(i, j): return Product(a[i], b[j]) - Product(a[j], b[i]) return as_vector((c(1, 2), c(2, 0), c(0, 1))) + def perp(self, o, a): + return as_vector([-a[1], a[0]]) + def dot(self, o, a, b): ai = indices(len(a.ufl_shape) - 1) bi = indices(len(b.ufl_shape) - 1) diff --git a/ufl/operators.py b/ufl/operators.py index 44dfea0c4..baf66f2bf 100644 --- a/ufl/operators.py +++ b/ufl/operators.py @@ -20,7 +20,7 @@ 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, + Transposed, Inner, Outer, Dot, Cross, Perp, Determinant, Inverse, Cofactor, Trace, Deviatoric, Skew, Sym) from ufl.coefficient import Coefficient from ufl.variable import Variable @@ -175,7 +175,7 @@ def perp(v): v = as_ufl(v) if v.ufl_shape != (2,): raise ValueError("Expecting a 2D vector expression.") - return as_vector((-v[1], v[0])) + return Perp(v) def cross(a, b): diff --git a/ufl/tensoralgebra.py b/ufl/tensoralgebra.py index 374931dc2..e1f40d8e5 100644 --- a/ufl/tensoralgebra.py +++ b/ufl/tensoralgebra.py @@ -217,6 +217,34 @@ def __str__(self): parstr(self.ufl_operands[1], self)) +@ufl_type(is_index_free=True, num_ops=1) +class Perp(CompoundTensorOperator): + __slots__ = () + + def __new__(cls, A): + sh = A.ufl_shape + + # Checks + if not len(sh) == 1: + raise ValueError(f"Perp requires arguments of rank 1, got {ufl_err_str(A)}") + if not sh[0] == 2: + raise ValueError(f"Perp can only work on 2D vectors, got {ufl_err_str(A)}") + + # Simplification + if isinstance(A, Zero): + return Zero(sh, A.ufl_free_indices, A.ufl_index_dimensions) + + return CompoundTensorOperator.__new__(cls) + + def __init__(self, A): + CompoundTensorOperator.__init__(self, (A,)) + + ufl_shape = (2,) + + def __str__(self): + return "perp(%s)" % self.ufl_operands[0] + + @ufl_type(num_ops=2) class Cross(CompoundTensorOperator): __slots__ = ("ufl_free_indices", "ufl_index_dimensions") From 0bafaf5bba781d31c2a32f3d8b48ab8e5fb59386 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Tue, 12 Sep 2023 10:05:20 +0100 Subject: [PATCH 08/11] Fix coefficients optional kwarg when calling a Form (#193) Co-authored-by: Matthew Scroggs From c01759f1f78b42e217a058c1c2bc8b7a4209bfaf Mon Sep 17 00:00:00 2001 From: "David A. Ham" Date: Tue, 12 Sep 2023 10:59:09 +0100 Subject: [PATCH 09/11] Add trimmed serendipity to element list. (#195) * Add trimmed serendipity * Fix value shape for trimmed serendipity * ufl plumbing update for trimmed serendipity. * Plumbing for SminusDiv.py * Adding in element stuff for SminusCurl. * Fix typo * remove spurioius names --------- Co-authored-by: Rob Kirby Co-authored-by: Justincrum Co-authored-by: nbouziani Co-authored-by: ksagiyam Co-authored-by: ksagiyam <46749170+ksagiyam@users.noreply.github.com> Co-authored-by: Connor Ward --- ufl/finiteelement/elementlist.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/ufl/finiteelement/elementlist.py b/ufl/finiteelement/elementlist.py index 6ea6a6fa2..73ca9fe90 100644 --- a/ufl/finiteelement/elementlist.py +++ b/ufl/finiteelement/elementlist.py @@ -193,6 +193,10 @@ def show_elements(): ("quadrilateral",)) register_element("BDMCF", None, 1, HDiv, "contravariant Piola", (1, None), ("quadrilateral",)) +register_element("SminusE", "SminusE", 1, HCurl, "covariant Piola", (1, None), cubes[1:3]) +register_element("SminusF", "SminusF", 1, HDiv, "contravariant Piola", (1, None), cubes[1:2]) +register_element("SminusDiv", "SminusDiv", 1, HDiv, "contravariant Piola", (1, None), cubes[1:3]) +register_element("SminusCurl", "SminusCurl", 1, HCurl, "covariant Piola", (1, None), cubes[1:3]) register_element("AAE", None, 1, HCurl, "covariant Piola", (1, None), ("hexahedron",)) register_element("AAF", None, 1, HDiv, "contravariant Piola", (1, None), From 8f44275dcf88f8c853487fa832c41098f9d9506c Mon Sep 17 00:00:00 2001 From: Ignacia Fierro-Piccardo Date: Tue, 12 Sep 2023 12:14:25 +0100 Subject: [PATCH 10/11] atan_2 changed to atan2 (#196) --- ufl/__init__.py | 6 +++--- ufl/algorithms/apply_derivatives.py | 2 +- ufl/algorithms/estimate_degrees.py | 2 +- ufl/mathfunctions.py | 4 ++-- ufl/operators.py | 4 ++-- 5 files changed, 9 insertions(+), 9 deletions(-) diff --git a/ufl/__init__.py b/ufl/__init__.py index 7cb9cf0e0..a31d341ee 100644 --- a/ufl/__init__.py +++ b/ufl/__init__.py @@ -185,7 +185,7 @@ - sqrt - exp, ln, erf - cos, sin, tan - - acos, asin, atan, atan_2 + - acos, asin, atan, atan2 - cosh, sinh, tanh - bessel_J, bessel_Y, bessel_I, bessel_K @@ -320,7 +320,7 @@ from ufl.operators import ( rank, shape, 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, + sqrt, exp, ln, erf, cos, sin, tan, acos, asin, atan, atan2, 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, variable, diff, Dx, grad, div, curl, rot, nabla_grad, nabla_div, Dn, exterior_derivative, @@ -385,7 +385,7 @@ 'transpose', 'tr', 'diag', 'diag_vector', 'dev', 'skew', 'sym', 'sqrt', 'exp', 'ln', 'erf', 'cos', 'sin', 'tan', - 'acos', 'asin', 'atan', 'atan_2', + 'acos', 'asin', 'atan', 'atan2', 'cosh', 'sinh', 'tanh', 'bessel_J', 'bessel_Y', 'bessel_I', 'bessel_K', 'eq', 'ne', 'le', 'ge', 'lt', 'gt', 'And', 'Or', 'Not', diff --git a/ufl/algorithms/apply_derivatives.py b/ufl/algorithms/apply_derivatives.py index 680969426..befb4067b 100644 --- a/ufl/algorithms/apply_derivatives.py +++ b/ufl/algorithms/apply_derivatives.py @@ -354,7 +354,7 @@ def atan(self, o, fp): f, = o.ufl_operands return fp / (1.0 + f**2) - def atan_2(self, o, fp, gp): + def atan2(self, o, fp, gp): f, g = o.ufl_operands return (g * fp - f * gp) / (f**2 + g**2) diff --git a/ufl/algorithms/estimate_degrees.py b/ufl/algorithms/estimate_degrees.py index 53dae9354..de39d26b3 100644 --- a/ufl/algorithms/estimate_degrees.py +++ b/ufl/algorithms/estimate_degrees.py @@ -242,7 +242,7 @@ def power(self, v, a, b): # negative integer, Coefficient, etc. return self._add_degrees(v, a, 2) - def atan_2(self, v, a, b): + def atan2(self, v, a, b): """Using the heuristic degree(atan2(const,const)) == 0 degree(atan2(a,b)) == max(degree(a),degree(b))+2 diff --git a/ufl/mathfunctions.py b/ufl/mathfunctions.py index c4bc9aca4..fc0619cd6 100644 --- a/ufl/mathfunctions.py +++ b/ufl/mathfunctions.py @@ -290,12 +290,12 @@ def evaluate(self, x, mapping, component, index_values): except TypeError: raise ValueError('Atan2 does not support complex numbers.') except ValueError: - warnings.warn('Value error in evaluation of function atan_2 with arguments %s, %s.' % (a, b)) + warnings.warn('Value error in evaluation of function atan2 with arguments %s, %s.' % (a, b)) raise return res def __str__(self): - return "atan_2(%s,%s)" % (self.ufl_operands[0], self.ufl_operands[1]) + return "atan2(%s,%s)" % (self.ufl_operands[0], self.ufl_operands[1]) @ufl_type() diff --git a/ufl/operators.py b/ufl/operators.py index baf66f2bf..dfdc01956 100644 --- a/ufl/operators.py +++ b/ufl/operators.py @@ -605,12 +605,12 @@ def atan(f): return _mathfunction(f, Atan) -def atan_2(f1, f2): +def atan2(f1, f2): "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.') + raise TypeError('atan2 is incompatible with complex numbers.') r = Atan2(f1, f2) if isinstance(r, (RealValue, Zero, int, float)): return float(r) From e2632367d14ca25e0246d93c2aaf253fc2fe3ed4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B8rgen=20Schartum=20Dokken?= Date: Tue, 12 Sep 2023 14:54:59 +0200 Subject: [PATCH 11/11] Group integrals with common integrand in IntegralData (#92) --------- Co-authored-by: Igor Baratta Co-authored-by: Matthew Scroggs Co-authored-by: David A. Ham --- .github/workflows/fenicsx-tests.yml | 4 +- test/test_domains.py | 57 +++++++++++++++++++++++++---- ufl/algorithms/compute_form_data.py | 13 +++---- ufl/algorithms/domain_analysis.py | 48 ++++++++++++++++-------- ufl/form.py | 12 +++++- 5 files changed, 102 insertions(+), 32 deletions(-) diff --git a/.github/workflows/fenicsx-tests.yml b/.github/workflows/fenicsx-tests.yml index 30d8094c2..39a2e8191 100644 --- a/.github/workflows/fenicsx-tests.yml +++ b/.github/workflows/fenicsx-tests.yml @@ -41,6 +41,8 @@ jobs: with: path: ./ffcx repository: FEniCS/ffcx + ref: dokken/group_integrals_v2 + - name: Install FFCx run: | cd ffcx @@ -77,7 +79,7 @@ jobs: - name: Install Basix and FFCx run: | python3 -m pip install git+https://github.com/FEniCS/basix.git - python3 -m pip install git+https://github.com/FEniCS/ffcx.git + python3 -m pip install git+https://github.com/FEniCS/ffcx.git@dokken/group_integrals_v2 - name: Clone DOLFINx uses: actions/checkout@v3 diff --git a/test/test_domains.py b/test/test_domains.py index b072274a6..2422a8822 100755 --- a/test/test_domains.py +++ b/test/test_domains.py @@ -4,6 +4,7 @@ Tests of domain language and attaching domains to forms. """ +from mockobjects import MockMesh, MockMeshFunction import pytest from ufl import * @@ -13,8 +14,6 @@ all_cells = (interval, triangle, tetrahedron, quadrilateral, hexahedron) -from mockobjects import MockMesh, MockMeshFunction - def test_construct_domains_from_cells(): for cell in all_cells: @@ -198,7 +197,8 @@ def test_everywhere_integrals_with_backwards_compatibility(): # Check some integral data assert ida.integral_type == "cell" - assert ida.subdomain_id == "otherwise" + assert(len(ida.subdomain_id) == 1) + assert ida.subdomain_id[0] == "otherwise" assert ida.metadata == {} # Integrands are not equal because of renumbering @@ -208,6 +208,47 @@ def test_everywhere_integrals_with_backwards_compatibility(): assert itg1.ufl_element() == itg2.ufl_element() +def test_merge_sort_integral_data(): + D = Mesh(triangle) + + V = FunctionSpace(D, FiniteElement("CG", triangle, 1)) + + u = Coefficient(V) + c = Constant(D) + a = c * dS((2, 4)) + u * dx + u * ds + 2 * u * dx(3) + 2 * c * dS + 2 * u * dx((1, 4)) + form_data = compute_form_data(a, do_append_everywhere_integrals=False).integral_data + assert(len(form_data) == 5) + + # Check some integral data + assert form_data[0].integral_type == "cell" + assert(len(form_data[0].subdomain_id) == 1) + assert form_data[0].subdomain_id[0] == "otherwise" + assert form_data[0].metadata == {} + + assert form_data[1].integral_type == "cell" + assert(len(form_data[1].subdomain_id) == 3) + assert form_data[1].subdomain_id[0] == 1 + assert form_data[1].subdomain_id[1] == 3 + assert form_data[1].subdomain_id[2] == 4 + assert form_data[1].metadata == {} + + assert form_data[2].integral_type == "exterior_facet" + assert(len(form_data[2].subdomain_id) == 1) + assert form_data[2].subdomain_id[0] == "otherwise" + assert form_data[2].metadata == {} + + assert form_data[3].integral_type == "interior_facet" + assert(len(form_data[3].subdomain_id) == 1) + assert form_data[3].subdomain_id[0] == "otherwise" + assert form_data[3].metadata == {} + + assert form_data[4].integral_type == "interior_facet" + assert(len(form_data[4].subdomain_id) == 2) + assert form_data[4].subdomain_id[0] == 2 + assert form_data[4].subdomain_id[1] == 4 + assert form_data[4].metadata == {} + + def xtest_mixed_elements_on_overlapping_regions(): # Old sketch, not working # Create domain and both disjoint and overlapping regions @@ -366,10 +407,10 @@ def xtest_form_domain_model(): # Old sketch, not working # domains and regions to be part of their integrands dxb = dx('DB') # Get Mesh by name dxbl = dx(Region(DB, (1, 4), 'DBL2')) - # Provide a region with different name but same subdomain ids as - # DBL + # Provide a region with different name but same subdomain ids as + # DBL dxbr = dx((1, 4)) - # Assume unique Mesh and provide subdomain ids explicitly + # Assume unique Mesh and provide subdomain ids explicitly # Not checking measure objects in detail, as long as # they carry information to construct integrals below @@ -466,8 +507,8 @@ def sub_elements_on_subdomains(W): # Disjunctified by UFL: alonly = dot(ul, vl) * dx(D1) - # integral_1 knows that only subelement VL is active + # integral_1 knows that only subelement VL is active am = (dot(ul, vl) + ur * vr) * dx(D2) - # integral_2 knows that both subelements are active + # integral_2 knows that both subelements are active aronly = ur * vr * \ dx(D3) # integral_3 knows that only subelement VR is active diff --git a/ufl/algorithms/compute_form_data.py b/ufl/algorithms/compute_form_data.py index 842d43d9f..348e27e17 100644 --- a/ufl/algorithms/compute_form_data.py +++ b/ufl/algorithms/compute_form_data.py @@ -102,13 +102,12 @@ def _compute_max_subdomain_ids(integral_data): max_subdomain_ids = {} for itg_data in integral_data: it = itg_data.integral_type - si = itg_data.subdomain_id - if isinstance(si, int): - newmax = si + 1 - else: - newmax = 0 - prevmax = max_subdomain_ids.get(it, 0) - max_subdomain_ids[it] = max(prevmax, newmax) + for integral in itg_data.integrals: + # Convert string for default integral to -1 + sids = (-1 if isinstance(si, str) else si for si in integral.subdomain_id()) + newmax = max(sids) + 1 + prevmax = max_subdomain_ids.get(it, 0) + max_subdomain_ids[it] = max(prevmax, newmax) return max_subdomain_ids diff --git a/ufl/algorithms/domain_analysis.py b/ufl/algorithms/domain_analysis.py index e12668ac9..5afec433a 100644 --- a/ufl/algorithms/domain_analysis.py +++ b/ufl/algorithms/domain_analysis.py @@ -12,11 +12,13 @@ import ufl from ufl.integral import Integral from ufl.form import Form +from ufl.protocols import id_or_none 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 +import typing class IntegralData(object): @@ -134,14 +136,18 @@ def integral_subdomain_ids(integral): raise ValueError(f"Invalid domain id {did}.") -def rearrange_integrals_by_single_subdomains(integrals, do_append_everywhere_integrals): +def rearrange_integrals_by_single_subdomains( + integrals: typing.List[Integral], + do_append_everywhere_integrals: bool) -> typing.Dict[int, typing.List[Integral]]: """Rearrange integrals over multiple subdomains to single subdomain integrals. Input: - integrals: list(Integral) + integrals: List of integrals + do_append_everywhere_integrals: Boolean indicating if integrals defined on the whole domain should + just be restricted to the set of input subdomain ids. Output: - integrals: dict: subdomain_id -> list(Integral) (reconstructed with single subdomain_id) + integrals: The integrals reconstructed with single subdomain_id """ # Split integrals into lists of everywhere and subdomain integrals everywhere_integrals = [] @@ -231,23 +237,37 @@ def build_integral_data(integrals): """ itgs = defaultdict(list) + # --- Merge integral data that has the same integrals, + unique_integrals = defaultdict(tuple) + metadata_table = defaultdict(dict) for integral in integrals: - domain = integral.ufl_domain() + integrand = integral.integrand() integral_type = integral.integral_type() + ufl_domain = integral.ufl_domain() + metadata = integral.metadata() + meta_hash = hash(canonicalize_metadata(metadata)) subdomain_id = integral.subdomain_id() + subdomain_data = id_or_none(integral.subdomain_data()) if subdomain_id == "everywhere": raise ValueError("'everywhere' not a valid subdomain id. Did you forget to call group_form_integrals?") + unique_integrals[(integral_type, ufl_domain, meta_hash, integrand, subdomain_data)] += (subdomain_id,) + metadata_table[(integral_type, ufl_domain, meta_hash, integrand, subdomain_data)] = metadata + + for integral_data, subdomain_ids in unique_integrals.items(): + (integral_type, ufl_domain, metadata, integrand, subdomain_data) = integral_data + + integral = Integral(integrand, integral_type, ufl_domain, subdomain_ids, + metadata_table[integral_data], subdomain_data) # Group for integral data (One integral data object for all - # integrals with same domain, itype, subdomain_id (but - # possibly different metadata). - itgs[(domain, integral_type, subdomain_id)].append(integral) + # integrals with same domain, itype, (but possibly different metadata). + itgs[(ufl_domain, integral_type, subdomain_ids)].append(integral) # Build list with canonical ordering, iteration over dicts # is not deterministic across python versions def keyfunc(item): (d, itype, sid), integrals = item - return (d._ufl_sort_key_(), itype, (type(sid).__name__, sid)) - + sid_int = tuple(-1 if i == "otherwise" else i for i in sid) + return (d._ufl_sort_key_(), itype, (type(sid).__name__, ), sid_int) integral_datas = [] for (d, itype, sid), integrals in sorted(itgs.items(), key=keyfunc): integral_datas.append(IntegralData(d, itype, sid, integrals, {})) @@ -262,8 +282,7 @@ def group_form_integrals(form, domains, do_append_everywhere_integrals=True): :returns: A new :class:`~.Form` with gathered integrands. """ # Group integrals by domain and type - integrals_by_domain_and_type = \ - group_integrals_by_domain_and_type(form.integrals(), domains) + integrals_by_domain_and_type = group_integrals_by_domain_and_type(form.integrals(), domains) integrals = [] for domain in domains: @@ -277,8 +296,8 @@ def group_form_integrals(form, domains, do_append_everywhere_integrals=True): # f*dx((1,2)) + g*dx((2,3)) -> f*dx(1) + (f+g)*dx(2) + g*dx(3) # (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, do_append_everywhere_integrals) + single_subdomain_integrals = rearrange_integrals_by_single_subdomains( + ddt_integrals, do_append_everywhere_integrals) for subdomain_id, ss_integrals in sorted_by_key(single_subdomain_integrals): @@ -307,8 +326,7 @@ def calc_hash(cd): # Accumulate integrands of integrals that share the # same compiler data - integrands_and_cds = \ - accumulate_integrands_with_same_metadata(samecd_integrals[1]) + integrands_and_cds = accumulate_integrands_with_same_metadata(samecd_integrals[1]) for integrand, metadata in integrands_and_cds: integral = Integral(integrand, integral_type, domain, diff --git a/ufl/form.py b/ufl/form.py index eeff85020..d4370b65a 100644 --- a/ufl/form.py +++ b/ufl/form.py @@ -10,7 +10,9 @@ # Modified by Anders Logg, 2009-2011. # Modified by Massimiliano Leoni, 2016. # Modified by Cecile Daversin-Catty, 2018. +# Modified by Jørgen S. Dokken 2023. +import numbers import warnings from collections import defaultdict from itertools import chain @@ -51,10 +53,18 @@ def _sorted_integrals(integrals): all_integrals = [] + def keyfunc(item): + if isinstance(item, numbers.Integral): + sid_int = item + else: + # As subdomain ids can be either int or tuples, we need to compare them + sid_int = tuple(-1 if i == "otherwise" else i for i in item) + return (type(item).__name__, sid_int) + # 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=keyfunc): unsorted_integrals = integrals_dict[d][it][si] # TODO: At this point we could order integrals by # metadata and integrand, or even add the