diff --git a/python/dune/codegen/blockstructured/geometry.py b/python/dune/codegen/blockstructured/geometry.py
index 3ad3a9026d8a5287e54e3950aef0aa199d08cc4e..08c19399314f3b7f9c2923a68882998319c23dbb 100644
--- a/python/dune/codegen/blockstructured/geometry.py
+++ b/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)