From 02cedcf44d6abe3ebfedb0501806307f6b7f8ddc Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B8rgen=20S=2E=20Dokken?= Date: Wed, 2 Oct 2024 16:02:25 +0000 Subject: [PATCH 1/6] Add handling of case where a split returns a form with no arguments (they have all been reduced to zeros). --- ufl/algorithms/formsplitter.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/ufl/algorithms/formsplitter.py b/ufl/algorithms/formsplitter.py index 00519963d..2273e8bd5 100644 --- a/ufl/algorithms/formsplitter.py +++ b/ufl/algorithms/formsplitter.py @@ -138,6 +138,8 @@ def extract_blocks(form, i: Optional[int] = None, j: Optional[None] = None): f = fs.split(form, pi, pj) if f.empty(): form_i.append(None) + elif f.arguments() == (): + form_i.append(None) else: form_i.append(f) forms.append(tuple(form_i)) From caedf443647484ffdeadc55abd81fce64cefd2b1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B8rgen=20S=2E=20Dokken?= Date: Wed, 2 Oct 2024 16:33:03 +0000 Subject: [PATCH 2/6] Add test case --- test/test_extract_blocks.py | 18 +++++++++++++++++- 1 file changed, 17 insertions(+), 1 deletion(-) diff --git a/test/test_extract_blocks.py b/test/test_extract_blocks.py index 4fdd923e7..f3855fe68 100644 --- a/test/test_extract_blocks.py +++ b/test/test_extract_blocks.py @@ -4,7 +4,6 @@ import ufl.algorithms from ufl.finiteelement import FiniteElement, MixedElement - def epsilon(u): return ufl.sym(ufl.grad(u)) @@ -80,3 +79,20 @@ def lhs(u, p, v, q): for j in range(2): J_ij_ext = ufl.extract_blocks(J, i, j) assert J_sub[2 * i + j].signature() == J_ij_ext.signature() + +def test_postive_restricted_extract_none(): + cell = ufl.triangle + d = cell.topological_dimension() + domain = ufl.Mesh(FiniteElement("Lagrange", cell, 1, (d,), ufl.identity_pullback, ufl.H1)) + el_u = FiniteElement("Lagrange", cell, 2, (d, ), ufl.identity_pullback, ufl.H1) + el_p = FiniteElement( + "Lagrange", cell, 1, (), ufl.identity_pullback, ufl.H1) + V = ufl.FunctionSpace(domain, el_u) + Q = ufl.FunctionSpace(domain, el_p) + W = ufl.MixedFunctionSpace(V, Q) + u, p = ufl.TrialFunctions(W) + v, q = ufl.TestFunctions(W) + a = ufl.inner(ufl.grad(u), ufl.grad(v))*ufl.dx + ufl.div(u)*q*ufl.dx + ufl.div(v)*p*ufl.dx + a += ufl.inner(u("+"),v("+")) * ufl.dS + a_blocks = ufl.extract_blocks(a) + assert a_blocks[1][1] is None \ No newline at end of file From 98ab60885c331fb9873c4542bc1e69f1d67c414c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B8rgen=20S=2E=20Dokken?= Date: Wed, 2 Oct 2024 16:43:38 +0000 Subject: [PATCH 3/6] Add test and some error handling --- test/test_extract_blocks.py | 19 ++++++++++++------- ufl/algorithms/formsplitter.py | 10 +++++++++- 2 files changed, 21 insertions(+), 8 deletions(-) diff --git a/test/test_extract_blocks.py b/test/test_extract_blocks.py index f3855fe68..0ce712418 100644 --- a/test/test_extract_blocks.py +++ b/test/test_extract_blocks.py @@ -4,6 +4,7 @@ import ufl.algorithms from ufl.finiteelement import FiniteElement, MixedElement + def epsilon(u): return ufl.sym(ufl.grad(u)) @@ -80,19 +81,23 @@ def lhs(u, p, v, q): J_ij_ext = ufl.extract_blocks(J, i, j) assert J_sub[2 * i + j].signature() == J_ij_ext.signature() + def test_postive_restricted_extract_none(): cell = ufl.triangle d = cell.topological_dimension() domain = ufl.Mesh(FiniteElement("Lagrange", cell, 1, (d,), ufl.identity_pullback, ufl.H1)) - el_u = FiniteElement("Lagrange", cell, 2, (d, ), ufl.identity_pullback, ufl.H1) - el_p = FiniteElement( - "Lagrange", cell, 1, (), ufl.identity_pullback, ufl.H1) + el_u = FiniteElement("Lagrange", cell, 2, (d,), ufl.identity_pullback, ufl.H1) + el_p = FiniteElement("Lagrange", cell, 1, (), ufl.identity_pullback, ufl.H1) V = ufl.FunctionSpace(domain, el_u) Q = ufl.FunctionSpace(domain, el_p) W = ufl.MixedFunctionSpace(V, Q) u, p = ufl.TrialFunctions(W) - v, q = ufl.TestFunctions(W) - a = ufl.inner(ufl.grad(u), ufl.grad(v))*ufl.dx + ufl.div(u)*q*ufl.dx + ufl.div(v)*p*ufl.dx - a += ufl.inner(u("+"),v("+")) * ufl.dS + v, q = ufl.TestFunctions(W) + a = ( + ufl.inner(ufl.grad(u), ufl.grad(v)) * ufl.dx + + ufl.div(u) * q * ufl.dx + + ufl.div(v) * p * ufl.dx + ) + a += ufl.inner(u("+"), v("+")) * ufl.dS a_blocks = ufl.extract_blocks(a) - assert a_blocks[1][1] is None \ No newline at end of file + assert a_blocks[1][1] is None diff --git a/ufl/algorithms/formsplitter.py b/ufl/algorithms/formsplitter.py index 2273e8bd5..0cb4287ad 100644 --- a/ufl/algorithms/formsplitter.py +++ b/ufl/algorithms/formsplitter.py @@ -10,6 +10,7 @@ from typing import Optional +from ufl.algorithms.apply_restrictions import apply_restrictions from ufl.algorithms.map_integrands import map_integrand_dags from ufl.argument import Argument from ufl.classes import FixedIndex, ListTensor @@ -139,7 +140,14 @@ def extract_blocks(form, i: Optional[int] = None, j: Optional[None] = None): if f.empty(): form_i.append(None) elif f.arguments() == (): - form_i.append(None) + # Restrict interior facet integrals to eliminate zeros + f = apply_restrictions(f) + if f.empty(): + form_i.append(None) + else: + raise RuntimeError( + "Extracted form with no arguments that does not reduce to 0." + ) else: form_i.append(f) forms.append(tuple(form_i)) From cd81e916fa2d766bb12faf755c01e0dabc7383f4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B8rgen=20S=2E=20Dokken?= Date: Wed, 2 Oct 2024 16:57:34 +0000 Subject: [PATCH 4/6] Apply restriction in formsplitter --- ufl/algorithms/formsplitter.py | 23 ++++++++++++----------- 1 file changed, 12 insertions(+), 11 deletions(-) diff --git a/ufl/algorithms/formsplitter.py b/ufl/algorithms/formsplitter.py index 0cb4287ad..a2627ad2e 100644 --- a/ufl/algorithms/formsplitter.py +++ b/ufl/algorithms/formsplitter.py @@ -10,8 +10,7 @@ from typing import Optional -from ufl.algorithms.apply_restrictions import apply_restrictions -from ufl.algorithms.map_integrands import map_integrand_dags +from ufl.algorithms.map_integrands import map_expr_dag, map_integrand_dags from ufl.argument import Argument from ufl.classes import FixedIndex, ListTensor from ufl.constantvalue import Zero @@ -82,6 +81,15 @@ def multi_index(self, obj): """Apply to multi_index.""" return obj + def restricted(self, o): + """Apply to a restricted function""" + # If we hit a restriction first apply form splitter to argument, then check for zero + op_split = map_expr_dag(self, o.ufl_operands[0]) + if isinstance(op_split, Zero): + return op_split + else: + return op_split(o._side) + expr = MultiFunction.reuse_if_untouched @@ -139,16 +147,9 @@ def extract_blocks(form, i: Optional[int] = None, j: Optional[None] = None): f = fs.split(form, pi, pj) if f.empty(): form_i.append(None) - elif f.arguments() == (): - # Restrict interior facet integrals to eliminate zeros - f = apply_restrictions(f) - if f.empty(): - form_i.append(None) - else: - raise RuntimeError( - "Extracted form with no arguments that does not reduce to 0." - ) else: + if (num_args := len(f.arguments())) != 2: + raise RuntimeError(f"Expected 2 arguments, got {num_args}") form_i.append(f) forms.append(tuple(form_i)) else: From 66239732695eb6ed05ce83fbc02aef701c017563 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B8rgen=20S=2E=20Dokken?= Date: Wed, 2 Oct 2024 16:59:48 +0000 Subject: [PATCH 5/6] Fix docstring --- ufl/algorithms/formsplitter.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ufl/algorithms/formsplitter.py b/ufl/algorithms/formsplitter.py index a2627ad2e..e9fefaf68 100644 --- a/ufl/algorithms/formsplitter.py +++ b/ufl/algorithms/formsplitter.py @@ -82,7 +82,7 @@ def multi_index(self, obj): return obj def restricted(self, o): - """Apply to a restricted function""" + """Apply to a restricted function.""" # If we hit a restriction first apply form splitter to argument, then check for zero op_split = map_expr_dag(self, o.ufl_operands[0]) if isinstance(op_split, Zero): From 72948fcb906562a0b9709e842e3878b1f3372bbb Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B8rgen=20S=2E=20Dokken?= Date: Mon, 7 Oct 2024 10:06:53 +0200 Subject: [PATCH 6/6] Fix apply pullback on piola mapped elements for double derivatives --- ufl/algorithms/apply_derivatives.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/ufl/algorithms/apply_derivatives.py b/ufl/algorithms/apply_derivatives.py index bd9624f12..da7b61da1 100644 --- a/ufl/algorithms/apply_derivatives.py +++ b/ufl/algorithms/apply_derivatives.py @@ -30,6 +30,8 @@ Imag, Indexed, IndexSum, + Jacobian, + JacobianDeterminant, JacobianInverse, ListTensor, Product, @@ -672,7 +674,7 @@ def reference_grad(self, o): f = o.ufl_operands[0] valid_operand = f._ufl_is_in_reference_frame_ or isinstance( - f, (JacobianInverse, SpatialCoordinate) + f, (JacobianInverse, SpatialCoordinate, Jacobian, JacobianDeterminant) ) if not valid_operand: raise ValueError("ReferenceGrad can only wrap a reference frame type!")