diff --git a/python/dune/codegen/blockstructured/geometry.py b/python/dune/codegen/blockstructured/geometry.py
index dba4940e1b03236db6bd78062b21acedd5a867e8..3631727497d71e538c1aa3602adf18c11266aab6 100644
--- a/python/dune/codegen/blockstructured/geometry.py
+++ b/python/dune/codegen/blockstructured/geometry.py
@@ -1,10 +1,9 @@
-from dune.codegen.pdelab.restriction import restricted_name
-
from dune.codegen.generation import (geometry_mixin,
get_backend,
temporary_variable,
instruction,
- get_global_context_value)
+ get_global_context_value,
+ domain)
from dune.codegen.tools import get_pymbolic_basename
from dune.codegen.options import (get_form_option,
option_switch, get_option)
@@ -13,7 +12,11 @@ from dune.codegen.pdelab.geometry import (AxiparallelGeometryMixin,
GenericPDELabGeometryMixin,
world_dimension,
local_dimension,
- component_iname)
+ component_iname,
+ enforce_boundary_restriction,
+ restricted_name,
+ name_geometry,
+ )
from dune.codegen.blockstructured.tools import sub_element_inames
import pymbolic.primitives as prim
from loopy.match import Writes
@@ -27,15 +30,49 @@ class BlockStructuredGeometryMixin(GenericPDELabGeometryMixin):
def facet_jacobian_determinant(self, o):
raise NotImplementedError("I think this case was never really implemented")
+ 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)
+
+ 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)})
+ )
+
def jacobian_inverse(self, o):
- raise NotImplementedError("I have not yet ported this one")
+ restriction = enforce_boundary_restriction(self)
+
+ assert(len(self.indices) == 2)
+ i, j = self.indices
+ self.indices = None
+
+ name = restricted_name("jit", restriction)
+ self.define_jacobian_inverse_transposed(name, restriction)
+
+ 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):
def facet_jacobian_determinant(self, o):
- return prim.Quotient(AxiparallelGeometryMixin.facet_jacobian_determinant(self, o),
- prim.Power(get_form_option("number_of_blocks"), local_dimension()))
+ raise NotImplementedError("I think this case was never really implemented")
def jacobian_determinant(self, o):
return prim.Quotient(AxiparallelGeometryMixin.jacobian_determinant(self, o),
@@ -97,8 +134,8 @@ def bilinear_transformation_coefficients():
raise NotImplementedError()
-def compute_jacobian(name):
- pymbolic_pos = get_backend("quad_pos", selector=option_switch("blockstructured"))()
+def compute_jacobian(name, visitor):
+ pymbolic_pos = visitor.quadrature_position_in_macro()
jac_iname = component_iname(context="jac")
coefficients = bilinear_transformation_coefficients()
@@ -137,23 +174,23 @@ def compute_jacobian(name):
for i, expr in enumerate(expr_jac):
instruction(expression=expr,
assignee=prim.Subscript(prim.Variable(name), (prim.Variable(jac_iname), i)),
- within_inames=frozenset((jac_iname, ) + sub_element_inames() + get_backend(interface="quad_inames")()),
+ within_inames=frozenset((jac_iname, ) + sub_element_inames() + visitor.quadrature_inames()),
depends_on=frozenset({Writes(get_pymbolic_basename(pymbolic_pos))} | {Writes(cd) for cd in coefficients})
)
-def define_jacobian_matrix(name):
+def define_jacobian_matrix(name, visitor):
temporary_variable(name, shape=(world_dimension(), world_dimension()), managed=True)
- compute_jacobian(name)
+ compute_jacobian(name, visitor)
-def name_jacobian_matrix():
+def name_jacobian_matrix(visitor):
name = "jacobian"
- define_jacobian_matrix(name)
+ define_jacobian_matrix(name, visitor)
return name
-def compute_determinant(name, matrix):
+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:
@@ -171,52 +208,36 @@ def compute_determinant(name, matrix):
raise NotImplementedError()
instruction(expression=expr_determinant,
assignee=prim.Variable(name),
- within_inames=frozenset(sub_element_inames() + get_backend(interface="quad_inames")()),
+ within_inames=frozenset(sub_element_inames() + visitor.quadrature_inames()),
depends_on=frozenset({Writes(matrix)})
)
-def define_jacobian_determinant(name):
+def define_jacobian_determinant(name, visitor):
temporary_variable(name, shape=(), managed=True)
- jacobian = name_jacobian_matrix()
- compute_determinant(name, jacobian)
+ jacobian = name_jacobian_matrix(visitor)
+ compute_determinant(name, jacobian, visitor)
-def define_jacobian_determinant_inverse(name):
+def define_jacobian_determinant_inverse(name, visitor):
temporary_variable(name, shape=(), managed=True)
- determinant = name_jacobian_determinant()
+ 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() + get_backend(interface="quad_inames")()),
- depends_on=frozenset({Writes(determinant)})
- )
-
-
-def define_jacobian_determinant_abs(name):
- temporary_variable(name, shape=(), managed=True)
- determinant = name_jacobian_determinant()
- return instruction(expression=prim.Call(prim.Variable("abs"), (prim.Variable(determinant),)),
- assignee=prim.Variable(name),
- within_inames=frozenset(sub_element_inames() + get_backend(interface="quad_inames")()),
+ within_inames=frozenset(sub_element_inames() + visitor.quadrature_inames()),
depends_on=frozenset({Writes(determinant)})
)
-def name_jacobian_determinant():
- name = "detjac"
- define_jacobian_determinant(name)
+def name_jacobian_determinant_signed(visitor):
+ name = "detjac_signed"
+ define_jacobian_determinant(name, visitor)
return name
-def name_jacobian_determinant_inverse():
+def name_jacobian_determinant_inverse(visitor):
name = "detjac_inverse"
- define_jacobian_determinant_inverse(name)
- return name
-
-
-def name_jacobian_determinant_abs():
- name = "detjac_abs"
- define_jacobian_determinant_abs(name)
+ define_jacobian_determinant_inverse(name, visitor)
return name
@@ -231,7 +252,7 @@ def pymbolic_jacobian_determinant():
prim.Power(get_form_option("number_of_blocks"), local_dimension()))
-def compute_inverse_transposed(name, det_inv, matrix):
+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)]
@@ -278,20 +299,20 @@ def compute_inverse_transposed(name, det_inv, matrix):
for i in range(dim):
instruction(expression=exprs[i][j],
assignee=assignee[i][j],
- within_inames=frozenset(sub_element_inames() + get_backend(interface="quad_inames")()),
+ within_inames=frozenset(sub_element_inames() + visitor.quadrature_inames()),
depends_on=frozenset({Writes(matrix), Writes(det_inv)}))
-def define_jacobian_inverse_transposed(name):
+def define_jacobian_inverse_transposed(name, visitor):
temporary_variable(name, shape=(world_dimension(), world_dimension()), managed=True)
- jacobian = name_jacobian_matrix()
- det_inv = name_jacobian_determinant_inverse()
- compute_inverse_transposed(name, det_inv, jacobian)
+ 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):
+def name_jacobian_inverse_transposed(restriction, visitor):
name = restricted_name("jit", restriction)
- define_jacobian_inverse_transposed(name)
+ define_jacobian_inverse_transposed(name, visitor)
return name
#
@@ -408,8 +429,8 @@ def apply_constant_to_global_transformation(name, local):
)
-def to_global(local):
- macro = name_point_in_macro(local)
+def to_global(local, visitor):
+ macro = name_point_in_macro(local, visitor)
name = macro + "_global"
it = get_global_context_value("integral_type")
diff --git a/test/blockstructured/poisson/poisson_unstructured.mini b/test/blockstructured/poisson/poisson_unstructured.mini
index c71c381ada38934e5794e2de566ff895e3db4bc9..4da0b9e19b8bdba5cef210281d96e872284279c0 100644
--- a/test/blockstructured/poisson/poisson_unstructured.mini
+++ b/test/blockstructured/poisson/poisson_unstructured.mini
@@ -19,9 +19,10 @@ compare_l2errorsquared = 1e-7
grid_unstructured = 1
[formcompiler.r]
-matrix_free = 1
+matrix_free = 0
blockstructured = 1
number_of_blocks = 16, 8 | expand dimension
+geometry_mixins = blockstructured_multilinear
[formcompiler.ufl_variants]
cell = quadrilateral, hexahedron | expand dimension