remove computation of inverse and determinante from geometry.py

python/dune/codegen/blockstructured/geometry.py

@@ -15,6 +15,7 @@ from dune.codegen.pdelab.geometry import (AxiparallelGeometryMixin,
                                           restricted_name,
                                           name_cell_geometry
                                           )
+from dune.codegen.pdelab.tensors import name_matrix_inverse, name_determinant
 from dune.codegen.blockstructured.tools import (sub_element_inames,
                                                 name_point_in_macro,
                                                 )
@@ -49,21 +50,11 @@ class BlockStructuredGeometryMixin(GenericPDELabGeometryMixin):
                              prim.Power(get_form_option("number_of_blocks"), local_dimension()))
 
     def jacobian_determinant(self, o):
-        name = 'detjac'
-        self.define_jacobian_determinant(name)
-        return prim.Quotient(prim.Variable(name),
-                             prim.Power(get_form_option("number_of_blocks"), local_dimension()))
-
-    def define_jacobian_determinant(self, name):
-        temporary_variable(name, shape=(), managed=True)
-
-        determinant_signed = name_jacobian_determinant_signed(self)
+        jacobian = name_jacobian_matrix(self)
+        name = name_determinant(jacobian, (world_dimension(), world_dimension()), self)
 
-        return instruction(expression=prim.Call(prim.Variable("abs"), (prim.Variable(determinant_signed),)),
-                           assignee=prim.Variable(name),
-                           within_inames=frozenset(sub_element_inames() + self.quadrature_inames()),
-                           depends_on=frozenset({Writes(determinant_signed)})
-                           )
+        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)
@@ -72,20 +63,12 @@ class BlockStructuredGeometryMixin(GenericPDELabGeometryMixin):
         i, j = self.indices
         self.indices = None
 
-        name = restricted_name("jit", restriction)
-        self.define_jacobian_inverse_transposed(name, restriction)
+        jacobian = restricted_name(name_jacobian_matrix(self), restriction)
+        name = name_matrix_inverse(jacobian, (world_dimension(), world_dimension()), self)
 
         return prim.Product((prim.Subscript(prim.Variable(name), (j, i)),
                              get_form_option("number_of_blocks")))
 
-    def define_jacobian_inverse_transposed(self, name, restriction):
-        temporary_variable(name, shape=(world_dimension(), world_dimension()), managed=True)
-
-        jacobian = name_jacobian_matrix(self)
-        det_inv = name_jacobian_determinant_inverse(self)
-
-        compute_inverse_transposed(name, det_inv, jacobian, self)
-
 
 @geometry_mixin("blockstructured_axiparallel")
 class AxiparallelBlockStructuredGeometryMixin(AxiparallelGeometryMixin, BlockStructuredGeometryMixin):
@@ -213,121 +196,6 @@ def name_jacobian_matrix(visitor):
     return name
 
 
-def compute_determinant(name, matrix, visitor):
-    dim = world_dimension()
-    matrix_entry = [[prim.Subscript(prim.Variable(matrix), (i, j)) for j in range(dim)] for i in range(dim)]
-    if dim == 2:
-        expr_determinant = prim.Sum((prim.Product((matrix_entry[0][0], matrix_entry[1][1])),
-                                     -1 * prim.Product((matrix_entry[1][0], matrix_entry[0][1]))))
-    elif dim == 3:
-        expr_determinant = prim.Sum((prim.Product((matrix_entry[0][0], matrix_entry[1][1], matrix_entry[2][2])),
-                                     prim.Product((matrix_entry[0][1], matrix_entry[1][2], matrix_entry[2][0])),
-                                     prim.Product((matrix_entry[0][2], matrix_entry[1][0], matrix_entry[2][1])),
-                                     -1 * prim.Product((matrix_entry[0][2], matrix_entry[1][1], matrix_entry[2][0])),
-                                     -1 * prim.Product((matrix_entry[0][0], matrix_entry[1][2], matrix_entry[2][1])),
-                                     -1 * prim.Product((matrix_entry[0][1], matrix_entry[1][0], matrix_entry[2][2]))
-                                     ))
-    else:
-        raise NotImplementedError()
-    instruction(expression=expr_determinant,
-                assignee=prim.Variable(name),
-                within_inames=frozenset(sub_element_inames() + visitor.quadrature_inames()),
-                depends_on=frozenset({Writes(matrix)})
-                )
-
-
-def define_jacobian_determinant(name, visitor):
-    temporary_variable(name, shape=(), managed=True)
-    jacobian = name_jacobian_matrix(visitor)
-    compute_determinant(name, jacobian, visitor)
-
-
-def define_jacobian_determinant_inverse(name, visitor):
-    temporary_variable(name, shape=(), managed=True)
-    determinant = name_jacobian_determinant_signed(visitor)
-    return instruction(expression=prim.Quotient(1., prim.Variable(determinant)),
-                       assignee=prim.Variable(name),
-                       within_inames=frozenset(sub_element_inames() + visitor.quadrature_inames()),
-                       depends_on=frozenset({Writes(determinant)})
-                       )
-
-
-def name_jacobian_determinant_signed(visitor):
-    name = "detjac_signed"
-    define_jacobian_determinant(name, visitor)
-    return name
-
-
-def name_jacobian_determinant_inverse(visitor):
-    name = "detjac_inverse"
-    define_jacobian_determinant_inverse(name, visitor)
-    return name
-
-
-def compute_inverse_transposed(name, det_inv, matrix, visitor):
-    dim = world_dimension()
-    matrix_entry = [[prim.Subscript(prim.Variable(matrix), (i, j)) for j in range(dim)] for i in range(dim)]
-    assignee = [[prim.Subscript(prim.Variable(name), (i, j)) for j in range(dim)] for i in range(dim)]
-    exprs = [[None for _ in range(dim)] for _ in range(dim)]
-
-    if dim == 2:
-        for i in range(2):
-            for j in range(2):
-                sign = 1. if i == j else -1.
-                exprs[i][j] = prim.Product((sign, prim.Variable(det_inv), matrix_entry[1 - i][1 - j]))
-    elif dim == 3:
-        exprs[0][0] = prim.Product((1., prim.Variable(det_inv),
-                                    prim.Sum((prim.Product((matrix_entry[1][1], matrix_entry[2][2])),
-                                              -1 * prim.Product((matrix_entry[1][2], matrix_entry[2][1]))))))
-        exprs[1][0] = prim.Product((-1., prim.Variable(det_inv),
-                                    prim.Sum((prim.Product((matrix_entry[0][1], matrix_entry[2][2])),
-                                              -1 * prim.Product((matrix_entry[0][2], matrix_entry[2][1]))))))
-        exprs[2][0] = prim.Product((1., prim.Variable(det_inv),
-                                    prim.Sum((prim.Product((matrix_entry[0][1], matrix_entry[1][2])),
-                                              -1 * prim.Product((matrix_entry[0][2], matrix_entry[1][1]))))))
-
-        exprs[0][1] = prim.Product((-1., prim.Variable(det_inv),
-                                    prim.Sum((prim.Product((matrix_entry[1][0], matrix_entry[2][2])),
-                                              -1 * prim.Product((matrix_entry[1][2], matrix_entry[2][0]))))))
-        exprs[1][1] = prim.Product((1., prim.Variable(det_inv),
-                                    prim.Sum((prim.Product((matrix_entry[0][0], matrix_entry[2][2])),
-                                              -1 * prim.Product((matrix_entry[0][2], matrix_entry[2][0]))))))
-        exprs[2][1] = prim.Product((-1., prim.Variable(det_inv),
-                                    prim.Sum((prim.Product((matrix_entry[0][0], matrix_entry[1][2])),
-                                              -1 * prim.Product((matrix_entry[0][2], matrix_entry[1][0]))))))
-
-        exprs[0][2] = prim.Product((1., prim.Variable(det_inv),
-                                    prim.Sum((prim.Product((matrix_entry[1][0], matrix_entry[2][1])),
-                                              -1 * prim.Product((matrix_entry[1][1], matrix_entry[2][0]))))))
-        exprs[1][2] = prim.Product((-1., prim.Variable(det_inv),
-                                    prim.Sum((prim.Product((matrix_entry[0][0], matrix_entry[2][1])),
-                                              -1 * prim.Product((matrix_entry[0][1], matrix_entry[2][0]))))))
-        exprs[2][2] = prim.Product((1., prim.Variable(det_inv),
-                                    prim.Sum((prim.Product((matrix_entry[0][0], matrix_entry[1][1])),
-                                              -1 * prim.Product((matrix_entry[0][1], matrix_entry[1][0]))))))
-    else:
-        raise NotImplementedError
-    for j in range(dim):
-        for i in range(dim):
-            instruction(expression=exprs[i][j],
-                        assignee=assignee[i][j],
-                        within_inames=frozenset(sub_element_inames() + visitor.quadrature_inames()),
-                        depends_on=frozenset({Writes(matrix), Writes(det_inv)}))
-
-
-def define_jacobian_inverse_transposed(name, visitor):
-    temporary_variable(name, shape=(world_dimension(), world_dimension()), managed=True)
-    jacobian = name_jacobian_matrix(visitor)
-    det_inv = name_jacobian_determinant_inverse(visitor)
-    compute_inverse_transposed(name, det_inv, jacobian, visitor)
-
-
-def name_jacobian_inverse_transposed(restriction, visitor):
-    name = restricted_name("jit", restriction)
-    define_jacobian_inverse_transposed(name, visitor)
-    return name
-
-
 def compute_multilinear_to_global_transformation(name, local, visitor):
     dim = world_dimension()
     temporary_variable(name, shape=(dim,), managed=True)