From f3a45e3db83d803f2486576e4911d2fc4882cf4a Mon Sep 17 00:00:00 2001 From: Joe Dean Date: Mon, 5 Jul 2021 21:49:55 +0100 Subject: [PATCH 01/10] Refactor duplicate code --- ufl/algorithms/apply_derivatives.py | 19 ++++++++++--------- 1 file changed, 10 insertions(+), 9 deletions(-) diff --git a/ufl/algorithms/apply_derivatives.py b/ufl/algorithms/apply_derivatives.py index 013618503..b42953598 100644 --- a/ufl/algorithms/apply_derivatives.py +++ b/ufl/algorithms/apply_derivatives.py @@ -506,9 +506,7 @@ def jacobian_inverse(self, o): return self.independent_terminal(o) if not o._ufl_is_terminal_: 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 = self._grad_to_reference_grad(o, o) return Do # TODO: Add more explicit geometry type handlers here, with @@ -552,9 +550,7 @@ def reference_value(self, o): error("ReferenceValue can only wrap a terminal") domain = f.ufl_domain() 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 = self._grad_to_reference_grad(o, K) return Do def reference_grad(self, o): @@ -566,9 +562,7 @@ def reference_grad(self, o): error("ReferenceGrad can only wrap a reference frame type!") domain = f.ufl_domain() 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 = self._grad_to_reference_grad(o, K) return Do # --- Nesting of gradients @@ -593,6 +587,13 @@ def _grad(self, o): # 1) n = count number of Grads, get f # 2) if not f.has_derivatives(n): return zero(...) + def _grad_to_reference_grad(self, o, K): + "Transforms grad to reference_grad" + r = indices(len(o.ufl_shape)) + i, j = indices(2) + Do = as_tensor(K[j, i] * ReferenceGrad(o)[r + (j,)], r + (i,)) + return Do + cell_avg = GenericDerivativeRuleset.independent_operator facet_avg = GenericDerivativeRuleset.independent_operator From e3e5964514ca760d8b2f459be4144497c4056717 Mon Sep 17 00:00:00 2001 From: Joe Dean Date: Mon, 5 Jul 2021 21:55:53 +0100 Subject: [PATCH 02/10] Transform grad(J) to reference_grad(J) --- ufl/algorithms/apply_derivatives.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/ufl/algorithms/apply_derivatives.py b/ufl/algorithms/apply_derivatives.py index b42953598..a17286d62 100644 --- a/ufl/algorithms/apply_derivatives.py +++ b/ufl/algorithms/apply_derivatives.py @@ -498,7 +498,10 @@ def geometric_quantity(self, o): else: # TODO: Which types does this involve? I don't think the # form compilers will handle this. - return Grad(o) + domain = o.ufl_domain() + K = JacobianInverse(domain) + Do = self._grad_to_reference_grad(o, K) + return Do def jacobian_inverse(self, o): # grad(K) == K_ji rgrad(K)_rj From bb2b76441b81a96ba97bc7a707d9964bc2cbe818 Mon Sep 17 00:00:00 2001 From: Joe Dean Date: Mon, 5 Jul 2021 22:00:00 +0100 Subject: [PATCH 03/10] Add comment --- ufl/algorithms/apply_derivatives.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/ufl/algorithms/apply_derivatives.py b/ufl/algorithms/apply_derivatives.py index a17286d62..9a85dd300 100644 --- a/ufl/algorithms/apply_derivatives.py +++ b/ufl/algorithms/apply_derivatives.py @@ -591,9 +591,10 @@ def _grad(self, o): # 2) if not f.has_derivatives(n): return zero(...) def _grad_to_reference_grad(self, o, K): - "Transforms grad to reference_grad" + "Transforms grad(o) to reference_grad(o)." r = indices(len(o.ufl_shape)) i, j = indices(2) + # grad(o) == K_ji rgrad(o)_rj Do = as_tensor(K[j, i] * ReferenceGrad(o)[r + (j,)], r + (i,)) return Do From d81c3b3c881d38a9c2804ca38c30ff48b90d4a78 Mon Sep 17 00:00:00 2001 From: Joe Dean Date: Mon, 5 Jul 2021 22:03:17 +0100 Subject: [PATCH 04/10] Change docstring --- ufl/algorithms/apply_derivatives.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/ufl/algorithms/apply_derivatives.py b/ufl/algorithms/apply_derivatives.py index 9a85dd300..75c67ed3a 100644 --- a/ufl/algorithms/apply_derivatives.py +++ b/ufl/algorithms/apply_derivatives.py @@ -491,7 +491,8 @@ def __init__(self, geometric_dimension): # --- Specialized rules for geometric quantities def geometric_quantity(self, o): - """Default for geometric quantities is dg/dx = 0 if piecewise constant, otherwise keep Grad(g). + """Default for geometric quantities is do/dx = 0 if piecewise constant, + otherwise transform grad(o) to reference_grad(o). Override for specific types if other behaviour is needed.""" if is_cellwise_constant(o): return self.independent_terminal(o) From 287ecb0f25363b8440a67401260958b37815228f Mon Sep 17 00:00:00 2001 From: Joe Dean Date: Mon, 5 Jul 2021 22:18:37 +0100 Subject: [PATCH 05/10] Improve docstring --- ufl/algorithms/apply_derivatives.py | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/ufl/algorithms/apply_derivatives.py b/ufl/algorithms/apply_derivatives.py index 75c67ed3a..a97de3432 100644 --- a/ufl/algorithms/apply_derivatives.py +++ b/ufl/algorithms/apply_derivatives.py @@ -592,7 +592,13 @@ def _grad(self, o): # 2) if not f.has_derivatives(n): return zero(...) def _grad_to_reference_grad(self, o, K): - "Transforms grad(o) to reference_grad(o)." + """Relates grad(o) to reference_grad(o) using the Jacobian inverse. + Parameters: + o: operand + K: Jacobian inverse + + Returns: + grad(o) written in terms of reference_grad(o) and K""" r = indices(len(o.ufl_shape)) i, j = indices(2) # grad(o) == K_ji rgrad(o)_rj From 4b82ef6cda37ea1ae3bd1ea97a25f53a75bebad0 Mon Sep 17 00:00:00 2001 From: Joe Dean Date: Mon, 5 Jul 2021 22:21:23 +0100 Subject: [PATCH 06/10] Fix docstring style --- ufl/algorithms/apply_derivatives.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/ufl/algorithms/apply_derivatives.py b/ufl/algorithms/apply_derivatives.py index a97de3432..7913d2850 100644 --- a/ufl/algorithms/apply_derivatives.py +++ b/ufl/algorithms/apply_derivatives.py @@ -593,11 +593,11 @@ def _grad(self, o): def _grad_to_reference_grad(self, o, K): """Relates grad(o) to reference_grad(o) using the Jacobian inverse. - Parameters: + Parameters o: operand K: Jacobian inverse - Returns: + Returns grad(o) written in terms of reference_grad(o) and K""" r = indices(len(o.ufl_shape)) i, j = indices(2) From ea7cbf950f5ba03495d44fbdc01daa7a6efee121 Mon Sep 17 00:00:00 2001 From: Joe Dean Date: Mon, 5 Jul 2021 22:23:24 +0100 Subject: [PATCH 07/10] Modify docstring --- ufl/algorithms/apply_derivatives.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ufl/algorithms/apply_derivatives.py b/ufl/algorithms/apply_derivatives.py index 7913d2850..c637d8c4c 100644 --- a/ufl/algorithms/apply_derivatives.py +++ b/ufl/algorithms/apply_derivatives.py @@ -492,7 +492,7 @@ def __init__(self, geometric_dimension): def geometric_quantity(self, o): """Default for geometric quantities is do/dx = 0 if piecewise constant, - otherwise transform grad(o) to reference_grad(o). + otherwise transform derivatives to reference derivatives. Override for specific types if other behaviour is needed.""" if is_cellwise_constant(o): return self.independent_terminal(o) From ae1dd488abacbe89690c2fb6b5a4e8aa68934c67 Mon Sep 17 00:00:00 2001 From: Joe Dean Date: Mon, 5 Jul 2021 22:31:38 +0100 Subject: [PATCH 08/10] pydocstyle fixes --- ufl/algorithms/apply_derivatives.py | 15 ++++++++++----- 1 file changed, 10 insertions(+), 5 deletions(-) diff --git a/ufl/algorithms/apply_derivatives.py b/ufl/algorithms/apply_derivatives.py index c637d8c4c..8b48be3f4 100644 --- a/ufl/algorithms/apply_derivatives.py +++ b/ufl/algorithms/apply_derivatives.py @@ -593,12 +593,17 @@ def _grad(self, o): def _grad_to_reference_grad(self, o, K): """Relates grad(o) to reference_grad(o) using the Jacobian inverse. - Parameters - o: operand - K: Jacobian inverse - Returns - grad(o) written in terms of reference_grad(o) and K""" + Args + ---- + o: operand + K: Jacobian inverse + + Returns + ------- + Do: grad(o) written in terms of reference_grad(o) and K + + """ r = indices(len(o.ufl_shape)) i, j = indices(2) # grad(o) == K_ji rgrad(o)_rj From 39cec450a0b4fb81b036281cb8a72ee1ba5dc8a9 Mon Sep 17 00:00:00 2001 From: Joe Dean Date: Wed, 7 Jul 2021 15:45:49 +0100 Subject: [PATCH 09/10] Move function --- ufl/algorithms/apply_derivatives.py | 43 +++++++++++++++-------------- 1 file changed, 22 insertions(+), 21 deletions(-) diff --git a/ufl/algorithms/apply_derivatives.py b/ufl/algorithms/apply_derivatives.py index 8b48be3f4..d7932df07 100644 --- a/ufl/algorithms/apply_derivatives.py +++ b/ufl/algorithms/apply_derivatives.py @@ -501,7 +501,7 @@ def geometric_quantity(self, o): # form compilers will handle this. domain = o.ufl_domain() K = JacobianInverse(domain) - Do = self._grad_to_reference_grad(o, K) + Do = grad_to_reference_grad(o, K) return Do def jacobian_inverse(self, o): @@ -510,7 +510,7 @@ def jacobian_inverse(self, o): return self.independent_terminal(o) if not o._ufl_is_terminal_: error("ReferenceValue can only wrap a terminal") - Do = self._grad_to_reference_grad(o, o) + Do = grad_to_reference_grad(o, o) return Do # TODO: Add more explicit geometry type handlers here, with @@ -554,7 +554,7 @@ def reference_value(self, o): error("ReferenceValue can only wrap a terminal") domain = f.ufl_domain() K = JacobianInverse(domain) - Do = self._grad_to_reference_grad(o, K) + Do = grad_to_reference_grad(o, K) return Do def reference_grad(self, o): @@ -566,7 +566,7 @@ def reference_grad(self, o): error("ReferenceGrad can only wrap a reference frame type!") domain = f.ufl_domain() K = JacobianInverse(domain) - Do = self._grad_to_reference_grad(o, K) + Do = grad_to_reference_grad(o, K) return Do # --- Nesting of gradients @@ -591,27 +591,28 @@ def _grad(self, o): # 1) n = count number of Grads, get f # 2) if not f.has_derivatives(n): return zero(...) - def _grad_to_reference_grad(self, o, K): - """Relates grad(o) to reference_grad(o) using the Jacobian inverse. + cell_avg = GenericDerivativeRuleset.independent_operator + facet_avg = GenericDerivativeRuleset.independent_operator - Args - ---- - o: operand - K: Jacobian inverse - Returns - ------- - Do: grad(o) written in terms of reference_grad(o) and K +def grad_to_reference_grad(o, K): + """Relates grad(o) to reference_grad(o) using the Jacobian inverse. - """ - r = indices(len(o.ufl_shape)) - i, j = indices(2) - # grad(o) == K_ji rgrad(o)_rj - Do = as_tensor(K[j, i] * ReferenceGrad(o)[r + (j,)], r + (i,)) - return Do + Args + ---- + o: Operand + K: Jacobian inverse - cell_avg = GenericDerivativeRuleset.independent_operator - facet_avg = GenericDerivativeRuleset.independent_operator + Returns + ------- + Do: grad(o) written in terms of reference_grad(o) and K + + """ + r = indices(len(o.ufl_shape)) + i, j = indices(2) + # grad(o) == K_ji rgrad(o)_rj + Do = as_tensor(K[j, i] * ReferenceGrad(o)[r + (j,)], r + (i,)) + return Do class ReferenceGradRuleset(GenericDerivativeRuleset): From 9855a098e0a0900536f24ab65e6cc558b81eb5a6 Mon Sep 17 00:00:00 2001 From: Joe Dean Date: Wed, 7 Jul 2021 15:47:45 +0100 Subject: [PATCH 10/10] Remove TODO --- ufl/algorithms/apply_derivatives.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/ufl/algorithms/apply_derivatives.py b/ufl/algorithms/apply_derivatives.py index d7932df07..d1fc39f99 100644 --- a/ufl/algorithms/apply_derivatives.py +++ b/ufl/algorithms/apply_derivatives.py @@ -497,8 +497,6 @@ def geometric_quantity(self, o): if is_cellwise_constant(o): return self.independent_terminal(o) else: - # TODO: Which types does this involve? I don't think the - # form compilers will handle this. domain = o.ufl_domain() K = JacobianInverse(domain) Do = grad_to_reference_grad(o, K)