diff --git a/python/dune/codegen/blockstructured/geometry.py b/python/dune/codegen/blockstructured/geometry.py index 3631727497d71e538c1aa3602adf18c11266aab6..3bf5acc4f81ef7451c2f6729b9222702febcc829 100644 --- a/python/dune/codegen/blockstructured/geometry.py +++ b/python/dune/codegen/blockstructured/geometry.py @@ -27,6 +27,20 @@ class BlockStructuredGeometryMixin(GenericPDELabGeometryMixin): def spatial_coordinate(self, o): return self.to_global(self.quadrature_position_in_macro()) + def to_global(self, local): + assert isinstance(local, prim.Expression) + name = get_pymbolic_basename(local) + "_global" + + it = get_global_context_value("integral_type") + + if it == 'cell': + compute_multilinear_to_global_transformation(name, local, self) + return prim.Variable(name) + elif it == 'exterior_facet' or it == 'interior_facet': + return GenericPDELabGeometryMixin.to_global(self, local) + else: + raise NotImplementedError + def facet_jacobian_determinant(self, o): raise NotImplementedError("I think this case was never really implemented") @@ -71,6 +85,20 @@ class BlockStructuredGeometryMixin(GenericPDELabGeometryMixin): @geometry_mixin("blockstructured_axiparallel") class AxiparallelBlockStructuredGeometryMixin(AxiparallelGeometryMixin, BlockStructuredGeometryMixin): + def to_global(self, local): + assert isinstance(local, prim.Expression) + name = get_pymbolic_basename(local) + "_global" + + it = get_global_context_value("integral_type") + + if it == 'cell': + compute_axiparallel_to_global_transformation(name, local, self) + return prim.Variable(name) + elif it == 'exterior_facet' or it == 'interior_facet': + return AxiparallelGeometryMixin.to_global(self, local) + else: + raise NotImplementedError + def facet_jacobian_determinant(self, o): raise NotImplementedError("I think this case was never really implemented") @@ -241,17 +269,6 @@ def name_jacobian_determinant_inverse(visitor): return name -# scale determinant according to the order of the blockstructuring -def pymbolic_jacobian_determinant(): - if get_form_option("constant_transformation_matrix"): - from dune.codegen.pdelab.geometry import name_jacobian_determinant - n_det = name_jacobian_determinant() - else: - n_det = name_jacobian_determinant_abs() - return prim.Quotient(prim.Variable(n_det), - prim.Power(get_form_option("number_of_blocks"), local_dimension())) - - 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)] @@ -369,7 +386,7 @@ def name_point_in_macro(point_in_micro, visitor): return name -def apply_default_to_global_transformation(name, local): +def compute_multilinear_to_global_transformation(name, local, visitor): dim = world_dimension() temporary_variable(name, shape=(dim,), managed=True) corners = name_element_corners() @@ -400,13 +417,13 @@ def apply_default_to_global_transformation(name, local): assignee = prim.Subscript(prim.Variable(name), (dim_pym,)) instruction(assignee=assignee, expression=expr, - within_inames=frozenset(sub_element_inames() + get_backend(interface="quad_inames")() + (dim_pym.name,)), + within_inames=frozenset(sub_element_inames() + visitor.quadrature_inames() + (dim_pym.name,)), within_inames_is_final=True, - depends_on=frozenset({Writes(local.name), Writes(corners)}) + depends_on=frozenset({Writes(get_pymbolic_basename(local)), Writes(corners)}) ) -def apply_constant_to_global_transformation(name, local): +def compute_axiparallel_to_global_transformation(name, local, visitor): dim = world_dimension() temporary_variable(name, shape=(dim,), managed=True) corners = name_element_corners() @@ -423,32 +440,12 @@ def apply_constant_to_global_transformation(name, local): assignee = prim.Subscript(prim.Variable(name), (dim_pym,)) instruction(assignee=assignee, expression=expr, - within_inames=frozenset(sub_element_inames() + get_backend(interface="quad_inames")() + (dim_pym.name,)), + within_inames=frozenset(sub_element_inames() + visitor.quadrature_inames() + (dim_pym.name,)), within_inames_is_final=True, - depends_on=frozenset({Writes(local.name), Writes(corners)}) + depends_on=frozenset({Writes(get_pymbolic_basename(local)), Writes(corners)}) ) -def to_global(local, visitor): - macro = name_point_in_macro(local, visitor) - name = macro + "_global" - - it = get_global_context_value("integral_type") - - if it == 'cell': - if get_form_option("constant_transformation_matrix"): - apply_constant_to_global_transformation(name, prim.Variable(macro)) - else: - apply_default_to_global_transformation(name, prim.Variable(macro)) - elif it == 'exterior_facet' or it == 'interior_facet': - from dune.codegen.pdelab.geometry import apply_to_global_transformation - apply_to_global_transformation(name, prim.Variable(macro)) - else: - raise NotImplementedError - - return prim.Variable(name) - - def define_element_corners(name): from dune.codegen.pdelab.driver import get_form n_corners = get_form().ufl_cell().num_vertices()