diff --git a/python/dune/codegen/blockstructured/geometry.py b/python/dune/codegen/blockstructured/geometry.py index 08c19399314f3b7f9c2923a68882998319c23dbb..753dba0f12ef227ab8f6db21b5bc64735e03988b 100644 --- a/python/dune/codegen/blockstructured/geometry.py +++ b/python/dune/codegen/blockstructured/geometry.py @@ -50,22 +50,32 @@ class BlockStructuredGeometryMixin(GenericPDELabGeometryMixin): prim.Power(get_form_option("number_of_blocks"), local_dimension())) def jacobian_determinant(self, o): + inames = self.inames + self.inames = inames + sub_element_inames() + jacobian = name_jacobian_matrix(self) name = name_determinant(jacobian, (world_dimension(), world_dimension()), self) + self.inames = inames + return prim.Quotient(prim.Call(prim.Variable("abs"), (prim.Variable(name),)), prim.Power(get_form_option("number_of_blocks"), local_dimension())) def jacobian_inverse(self, o): - restriction = enforce_boundary_restriction(self) - assert(len(self.indices) == 2) i, j = self.indices self.indices = None + restriction = enforce_boundary_restriction(self) + + inames = self.inames + self.inames = inames + sub_element_inames() + jacobian = restricted_name(name_jacobian_matrix(self), restriction) name = name_matrix_inverse(jacobian, (world_dimension(), world_dimension()), self) + self.inames = inames + return prim.Product((prim.Subscript(prim.Variable(name), (j, i)), get_form_option("number_of_blocks"))) diff --git a/python/dune/codegen/pdelab/tensors.py b/python/dune/codegen/pdelab/tensors.py index 5a7f48c23aeb8b5451be30a493a5a2312d273dee..fa9d52e6f4bd83aa8ca93f03b4167d5754a08a9e 100644 --- a/python/dune/codegen/pdelab/tensors.py +++ b/python/dune/codegen/pdelab/tensors.py @@ -38,7 +38,7 @@ def define_determinant(name, matrix, shape, visitor): raise NotImplementedError() instruction(expression=expr_determinant, assignee=prim.Variable(name), - within_inames=frozenset(visitor.quadrature_inames()), + within_inames=frozenset(visitor.inames), depends_on=frozenset({Writes(matrix)}) ) @@ -50,8 +50,8 @@ def define_determinant_inverse(name, matrix, shape, visitor): instruction(expression=prim.Quotient(1, prim.Variable(det)), assignee=prim.Variable(name), - within_inames=frozenset(visitor.quadrature_inames()), - depends_on=frozenset({Writes(matrix)}) + within_inames=frozenset(visitor.inames), + depends_on=frozenset({Writes(matrix), Writes(det)}) ) @@ -108,7 +108,7 @@ def define_matrix_inverse(name, name_inv, shape, visitor): for i in range(dim): instruction(expression=exprs[i][j], assignee=assignee[i][j], - within_inames=frozenset(visitor.quadrature_inames()), + within_inames=frozenset(visitor.inames), depends_on=frozenset({Writes(name), Writes(det_inv)})) @@ -150,7 +150,7 @@ def define_assembled_tensor(name, expr, visitor): visitor.indices = indices instruction(assignee=prim.Subscript(prim.Variable(name), indices), expression=visitor.call(expr), - forced_iname_deps=frozenset(visitor.quadrature_inames()), + forced_iname_deps=frozenset(visitor.inames), depends_on=frozenset({lp.match.Tagged("sumfact_stage1")}), tags=frozenset({"quad"}), ) diff --git a/python/dune/codegen/ufl/visitor.py b/python/dune/codegen/ufl/visitor.py index e55232d26839abb2a8866dc2a57af715156237a7..ed4635dc26995b2d80632d8f0e161b6c85bb6caa 100644 --- a/python/dune/codegen/ufl/visitor.py +++ b/python/dune/codegen/ufl/visitor.py @@ -56,7 +56,9 @@ class UFL2LoopyVisitor(ModifiedTerminalTracker): def accumulate(self, o): for info in self.list_accumulation_infos(o): self.current_info = info + self.inames = self.quadrature_inames() expr = self._call(o, False) + self.inames = None if expr != 0: if get_form_option("simplify"): from dune.codegen.sympy import simplify_pymbolic_expression @@ -70,7 +72,6 @@ class UFL2LoopyVisitor(ModifiedTerminalTracker): self._indices_backup = [] self.test_info = None self.trial_info = None - self.inames = () self.do_predicates = do_predicates return self.call(o)