diff --git a/python/dune/codegen/blockstructured/__init__.py b/python/dune/codegen/blockstructured/__init__.py
index 8810750ff4f0d1122429917c48f304b05d101667..dc9db5e5a7204099e6e8345a45dd5889de66f589 100644
--- a/python/dune/codegen/blockstructured/__init__.py
+++ b/python/dune/codegen/blockstructured/__init__.py
@@ -1,3 +1,4 @@
+import dune.codegen.blockstructured.accumulation
import dune.codegen.blockstructured.quadrature
import dune.codegen.blockstructured.argument
import dune.codegen.blockstructured.geometry
diff --git a/python/dune/codegen/blockstructured/basis.py b/python/dune/codegen/blockstructured/basis.py
index 3477a18295fdae14a324a8c309635bb6df565537..e5aa4d20ebef0eb320fff559c2d2fb8b54fe909b 100644
--- a/python/dune/codegen/blockstructured/basis.py
+++ b/python/dune/codegen/blockstructured/basis.py
@@ -1,4 +1,5 @@
from dune.codegen.generation import (backend,
+ basis_mixin,
kernel_cached,
get_backend,
instruction,
@@ -9,7 +10,8 @@ from dune.codegen.generation import (backend,
include_file,)
from dune.codegen.tools import get_pymbolic_basename
from dune.codegen.loopy.target import type_floatingpoint
-from dune.codegen.pdelab.basis import (declare_cache_temporary,
+from dune.codegen.pdelab.basis import (GenericBasisMixin,
+ declare_cache_temporary,
name_localbasis_cache,
type_localbasis,
FEM_name_mangling)
@@ -27,6 +29,56 @@ from ufl import MixedElement
import pymbolic.primitives as prim
+@basis_mixin("blockstructured")
+class BlockStructuredBasisMixin(GenericBasisMixin):
+ def lfs_inames(self, element, restriction, number, context=""):
+ return lfs_inames(element, restriction, number, context=context)
+
+ def implement_basis(self, element, restriction, number, context=''):
+ assert not isinstance(element, MixedElement)
+ name = "phi_{}".format(FEM_name_mangling(element))
+ name = restricted_name(name, restriction)
+ self.evaluate_basis(element, name, restriction)
+ inames = self.lfs_inames(element, restriction, number, context=context)
+
+ return prim.Subscript(prim.Variable(name), (tensor_index_to_sequential_index(inames, element.degree() + 1), 0))
+
+ @kernel_cached
+ def evaluate_basis(self, element, name, restriction):
+ temporary_variable(name, shape=((element.degree() + 1)**world_dimension(), 1),
+ decl_method=declare_cache_temporary(element, restriction, 'Function'))
+ cache = name_localbasis_cache(element)
+ qp = self.to_cell(self.quadrature_position())
+ localbasis = name_localbasis(element)
+ instruction(inames=self.quadrature_inames(),
+ code='{} = {}.evaluateFunction({}, {});'.format(name, cache, str(qp), localbasis),
+ assignees=frozenset({name}),
+ read_variables=frozenset({get_pymbolic_basename(qp)}),
+ )
+
+ def implement_reference_gradient(self, element, restriction, number, context=''):
+ assert not isinstance(element, MixedElement)
+ name = "js_{}".format(FEM_name_mangling(element))
+ name = restricted_name(name, restriction)
+ self.evaluate_reference_gradient(element, name, restriction)
+ inames = self.lfs_inames(element, restriction, number, context=context)
+
+ return prim.Subscript(prim.Variable(name), (tensor_index_to_sequential_index(inames, element.degree() + 1), 0))
+
+ @kernel_cached
+ def evaluate_reference_gradient(self, element, name, restriction):
+ temporary_variable(name, shape=((element.degree() + 1)**world_dimension(), 1, world_dimension()),
+ decl_method=declare_cache_temporary(element, restriction, 'Jacobian'))
+ cache = name_localbasis_cache(element)
+ qp = self.to_cell(self.quadrature_position())
+ localbasis = name_localbasis(element)
+ instruction(inames=self.quadrature_inames(),
+ code='{} = {}.evaluateJacobian({}, {});'.format(name, cache, str(qp), localbasis),
+ assignees=frozenset({name}),
+ read_variables=frozenset({get_pymbolic_basename(qp)}),
+ )
+
+
# define FE basis explicitly in localoperator
@backend(interface="typedef_localbasis", name="blockstructured")
@class_member(classtag="operator")
@@ -68,49 +120,50 @@ def name_localbasis(leaf_element):
return name
-@kernel_cached
-def evaluate_basis(leaf_element, name, restriction):
- temporary_variable(name, shape=((leaf_element.degree() + 1)**world_dimension(), 1),
- decl_method=declare_cache_temporary(leaf_element, restriction, 'Function'))
- cache = name_localbasis_cache(leaf_element)
- qp = pymbolic_quadrature_position_in_cell(restriction)
- localbasis = name_localbasis(leaf_element)
- instruction(inames=get_backend("quad_inames")(),
- code='{} = {}.evaluateFunction({}, {});'.format(name, cache, str(qp), localbasis),
- assignees=frozenset({name}),
- read_variables=frozenset({get_pymbolic_basename(qp)}),
- )
-
-
-def pymbolic_basis(leaf_element, restriction, number, context=''):
- assert not isinstance(leaf_element, MixedElement)
- name = "phi_{}".format(FEM_name_mangling(leaf_element))
- name = restricted_name(name, restriction)
- evaluate_basis(leaf_element, name, restriction)
- inames = lfs_inames(leaf_element, restriction, number, context=context)
-
- return prim.Subscript(prim.Variable(name), (tensor_index_to_sequential_index(inames, leaf_element.degree() + 1), 0))
-
-
-@kernel_cached
-def evaluate_reference_gradient(leaf_element, name, restriction):
- temporary_variable(name, shape=((leaf_element.degree() + 1)**world_dimension(), 1, world_dimension()),
- decl_method=declare_cache_temporary(leaf_element, restriction, 'Jacobian'))
- cache = name_localbasis_cache(leaf_element)
- qp = pymbolic_quadrature_position_in_cell(restriction)
- localbasis = name_localbasis(leaf_element)
- instruction(inames=get_backend("quad_inames")(),
- code='{} = {}.evaluateJacobian({}, {});'.format(name, cache, str(qp), localbasis),
- assignees=frozenset({name}),
- read_variables=frozenset({get_pymbolic_basename(qp)}),
- )
-
-
-def pymbolic_reference_gradient(leaf_element, restriction, number, context=''):
- assert not isinstance(leaf_element, MixedElement)
- name = "js_{}".format(FEM_name_mangling(leaf_element))
- name = restricted_name(name, restriction)
- evaluate_reference_gradient(leaf_element, name, restriction)
- inames = lfs_inames(leaf_element, restriction, number, context=context)
-
- return prim.Subscript(prim.Variable(name), (tensor_index_to_sequential_index(inames, leaf_element.degree() + 1), 0))
+# @kernel_cached
+# def evaluate_basis(leaf_element, name, restriction):
+# temporary_variable(name, shape=((leaf_element.degree() + 1)**world_dimension(), 1),
+# decl_method=declare_cache_temporary(leaf_element, restriction, 'Function'))
+# cache = name_localbasis_cache(leaf_element)
+# qp = pymbolic_quadrature_position_in_cell(restriction)
+# localbasis = name_localbasis(leaf_element)
+# instruction(inames=get_backend("quad_inames")(),
+# code='{} = {}.evaluateFunction({}, {});'.format(name, cache, str(qp), localbasis),
+# assignees=frozenset({name}),
+# read_variables=frozenset({get_pymbolic_basename(qp)}),
+# )
+#
+#
+# def pymbolic_basis(leaf_element, restriction, number, context=''):
+# assert not isinstance(leaf_element, MixedElement)
+# name = "phi_{}".format(FEM_name_mangling(leaf_element))
+# name = restricted_name(name, restriction)
+# evaluate_basis(leaf_element, name, restriction)
+# inames = lfs_inames(leaf_element, restriction, number, context=context)
+#
+# return prim.Subscript(prim.Variable(name), (tensor_index_to_sequential_index(inames, leaf_element.degree() + 1), 0))
+
+#
+# @kernel_cached
+# def evaluate_reference_gradient(leaf_element, name, restriction):
+# from pudb import set_trace; set_trace()
+# temporary_variable(name, shape=((leaf_element.degree() + 1)**world_dimension(), 1, world_dimension()),
+# decl_method=declare_cache_temporary(leaf_element, restriction, 'Jacobian'))
+# cache = name_localbasis_cache(leaf_element)
+# qp = pymbolic_quadrature_position_in_cell(restriction)
+# localbasis = name_localbasis(leaf_element)
+# instruction(inames=get_backend("quad_inames")(),
+# code='{} = {}.evaluateJacobian({}, {});'.format(name, cache, str(qp), localbasis),
+# assignees=frozenset({name}),
+# read_variables=frozenset({get_pymbolic_basename(qp)}),
+# )
+#
+#
+# def pymbolic_reference_gradient(leaf_element, restriction, number, context=''):
+# assert not isinstance(leaf_element, MixedElement)
+# name = "js_{}".format(FEM_name_mangling(leaf_element))
+# name = restricted_name(name, restriction)
+# evaluate_reference_gradient(leaf_element, name, restriction)
+# inames = lfs_inames(leaf_element, restriction, number, context=context)
+#
+# return prim.Subscript(prim.Variable(name), (tensor_index_to_sequential_index(inames, leaf_element.degree() + 1), 0))
diff --git a/python/dune/codegen/blockstructured/geometry.py b/python/dune/codegen/blockstructured/geometry.py
index 4bc0e6131ccfd192f14cc6e2ab15f37b99a7841b..f19202ec8ed3a5b527302edab638d841b336663e 100644
--- a/python/dune/codegen/blockstructured/geometry.py
+++ b/python/dune/codegen/blockstructured/geometry.py
@@ -1,13 +1,17 @@
from dune.codegen.pdelab.restriction import restricted_name
-from dune.codegen.generation import (get_backend,
+from dune.codegen.generation import (geometry_mixin,
+ get_backend,
temporary_variable,
instruction,
get_global_context_value)
from dune.codegen.tools import get_pymbolic_basename
from dune.codegen.options import (get_form_option,
option_switch, get_option)
-from dune.codegen.pdelab.geometry import (world_dimension,
+from dune.codegen.pdelab.geometry import (AxiparallelGeometryMixin,
+ EquidistantGeometryMixin,
+ GenericPDELabGeometryMixin,
+ world_dimension,
local_dimension,
component_iname)
from dune.codegen.blockstructured.tools import sub_element_inames
@@ -15,6 +19,31 @@ import pymbolic.primitives as prim
from loopy.match import Writes
+@geometry_mixin("blockstructured_multilinear")
+class BlockStructuredGeometryMixin(GenericPDELabGeometryMixin):
+ def facet_jacobian_determinant(self, o):
+ raise NotImplementedError("I think this case was never really implemented")
+
+ def jacobian_inverse(self, o):
+ raise NotImplementedError("I have not yet ported this one")
+
+
+@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()))
+
+ def jacobian_inverse(self, o):
+ return prim.Quotient(AxiparallelGeometryMixin.facet_jacobian_determinant(self, o),
+ prim.Power(get_form_option("number_of_blocks"), local_dimension()))
+
+
+@geometry_mixin("blockstructured_equidistant")
+class EquidistantBlockStructuredGeometryMixin(EquidistantGeometryMixin, AxiparallelBlockStructuredGeometryMixin):
+ pass
+
+
# TODO warum muss within_inames_is_final=True gesetzt werden?
def compute_corner_diff(first, second, additional_terms=tuple()):
corners = name_element_corners()
@@ -258,23 +287,23 @@ def name_jacobian_inverse_transposed(restriction):
define_jacobian_inverse_transposed(name)
return name
-
-# scale Jacobian according to the order of the blockstructure
-def pymbolic_jacobian_inverse(i, j, restriction):
- if get_form_option("constant_transformation_matrix"):
- from dune.codegen.pdelab.geometry import name_constant_jacobian_inverse_transposed
- name_jit = name_constant_jacobian_inverse_transposed(restriction)
- else:
- name_jit = name_jacobian_inverse_transposed(restriction)
- return prim.Product((get_form_option("number_of_blocks"),
- prim.Subscript(prim.Variable(name_jit), (j, i))))
-
-
-# scale determinant according to the order of the blockstructure
-def pymbolic_facet_jacobian_determinant():
- from dune.codegen.pdelab.geometry import name_facet_jacobian_determinant
- return prim.Quotient(prim.Variable(name_facet_jacobian_determinant()),
- prim.Power(get_form_option("number_of_blocks"), local_dimension()))
+#
+# # scale Jacobian according to the order of the blockstructure
+# def pymbolic_jacobian_inverse(i, j, restriction):
+# if get_form_option("constant_transformation_matrix"):
+# from dune.codegen.pdelab.geometry import name_constant_jacobian_inverse_transposed
+# name_jit = name_constant_jacobian_inverse_transposed(restriction)
+# else:
+# name_jit = name_jacobian_inverse_transposed(restriction)
+# return prim.Product((get_form_option("number_of_blocks"),
+# prim.Subscript(prim.Variable(name_jit), (j, i))))
+#
+#
+# # scale determinant according to the order of the blockstructure
+# def pymbolic_facet_jacobian_determinant():
+# from dune.codegen.pdelab.geometry import name_facet_jacobian_determinant
+# return prim.Quotient(prim.Variable(name_facet_jacobian_determinant()),
+# prim.Power(get_form_option("number_of_blocks"), local_dimension()))
# translate a point in the micro element into macro coordinates
diff --git a/python/dune/codegen/blockstructured/quadrature.py b/python/dune/codegen/blockstructured/quadrature.py
index 64282c687bc4d41a2b5c5fe28303cbc555a5ddae..d2875435fb72345a3bd87ee7e23983edbd9343fe 100644
--- a/python/dune/codegen/blockstructured/quadrature.py
+++ b/python/dune/codegen/blockstructured/quadrature.py
@@ -1,29 +1,30 @@
from dune.codegen.generation import (backend,
quadrature_mixin,
)
-from dune.codegen.pdelab.quadrature import (GenericQuadratureMixin,
- quadrature_iname)
+from dune.codegen.pdelab.quadrature import GenericQuadratureMixin
from dune.codegen.blockstructured.geometry import name_point_in_macro
import pymbolic.primitives as prim
@quadrature_mixin("blockstructured")
class BlockstructuredQuadratureMixin(GenericQuadratureMixin):
- pass
-
-
-@backend(interface="quad_pos", name='blockstructured')
-def pymbolic_quadrature_position():
- quad_points = name_quadrature_points()
- quad_iname = quadrature_iname()
-
- qp_in_micro = prim.Subscript(prim.Variable(quad_points), (prim.Variable(quad_iname),))
-
- return prim.Variable(name_point_in_macro(qp_in_micro))
-
-
-@backend(interface="qp_in_cell", name="blockstructured")
-def pymbolic_quadrature_position_in_cell(restriction):
- from dune.codegen.pdelab.geometry import to_cell_coordinates
- quad_pos = pymbolic_quadrature_position()
- return to_cell_coordinates(quad_pos, restriction)
+ def quadrature_position(self):
+ original = GenericQuadratureMixin.quadrature_position(self)
+ return prim.Variable(name_point_in_macro(original))
+
+#
+# @backend(interface="quad_pos", name='blockstructured')
+# def pymbolic_quadrature_position():
+# quad_points = name_quadrature_points()
+# quad_iname = quadrature_iname()
+#
+# qp_in_micro = prim.Subscript(prim.Variable(quad_points), (prim.Variable(quad_iname),))
+#
+# return prim.Variable(name_point_in_macro(qp_in_micro))
+#
+#
+# @backend(interface="qp_in_cell", name="blockstructured")
+# def pymbolic_quadrature_position_in_cell(restriction):
+# from dune.codegen.pdelab.geometry import to_cell_coordinates
+# quad_pos = pymbolic_quadrature_position()
+# return to_cell_coordinates(quad_pos, restriction)
diff --git a/python/dune/codegen/options.py b/python/dune/codegen/options.py
index 3e41ef98368ac2cf75517dec5824a4a3440ee6b0..83b57cf48b036302d0ff64e196ffe36062cfc0fa 100644
--- a/python/dune/codegen/options.py
+++ b/python/dune/codegen/options.py
@@ -184,6 +184,12 @@ def process_form_options(opt, form):
accumulation_mixins="sumfact",
)
+ if opt.blockstructured:
+ opt = opt.copy(accumulation_mixins="blockstructured",
+ geometry_mixins="blockstructured",
+ quadrature_mixins="blockstructured",
+ )
+
if opt.control:
opt = opt.copy(accumulation_mixins="control")
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index 7f2fef85e537fbe264712de625dc4f4b61208ad8..4be79371bc51e40c4334f0e91343c522d00d2147 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -7,5 +7,5 @@ add_subdirectory(stokes)
add_subdirectory(sumfact)
add_subdirectory(coeffeval)
-#add_subdirectory(blockstructured)
+add_subdirectory(blockstructured)
add_subdirectory(adjoint)