diff --git a/python/dune/codegen/__init__.py b/python/dune/codegen/__init__.py
index 243911af33e57531a147b9b4814b14ad65fd408b..e892baf44c4c9fb301536c15c40059a5acc5d602 100644
--- a/python/dune/codegen/__init__.py
+++ b/python/dune/codegen/__init__.py
@@ -6,8 +6,7 @@ os.environ["OMP_NUM_THREADS"] = "1"
# Trigger imports that involve monkey patching!
import dune.codegen.loopy.symbolic # noqa
-# Trigger some imports that are needed to have all backend implementations visible
-# to the selection mechanisms
+# Trigger some imports that are needed to have all mixin implementations visible
import dune.codegen.pdelab # noqa
import dune.codegen.sumfact # noqa
import dune.codegen.blockstructured # noqa
diff --git a/python/dune/codegen/blockstructured/__init__.py b/python/dune/codegen/blockstructured/__init__.py
index 138b11463a554ee58ac236d4439120172aabb640..dc9db5e5a7204099e6e8345a45dd5889de66f589 100644
--- a/python/dune/codegen/blockstructured/__init__.py
+++ b/python/dune/codegen/blockstructured/__init__.py
@@ -1,62 +1,7 @@
+import dune.codegen.blockstructured.accumulation
import dune.codegen.blockstructured.quadrature
import dune.codegen.blockstructured.argument
import dune.codegen.blockstructured.geometry
import dune.codegen.blockstructured.spaces
import dune.codegen.blockstructured.basis
import dune.codegen.blockstructured.transformations
-from dune.codegen.options import get_form_option
-from dune.codegen.pdelab.quadrature import pymbolic_quadrature_position
-from dune.codegen.blockstructured.spaces import lfs_inames
-from dune.codegen.blockstructured.basis import (pymbolic_reference_gradient,
- pymbolic_basis)
-from dune.codegen.blockstructured.geometry import (pymbolic_jacobian_inverse,
- pymbolic_jacobian_determinant,
- pymbolic_facet_jacobian_determinant,
- to_global)
-from dune.codegen.blockstructured.tools import sub_element_inames
-from dune.codegen.pdelab import PDELabInterface
-
-
-class BlockStructuredInterface(PDELabInterface):
- def __init__(self):
- PDELabInterface.__init__(self)
-
- def generate_accumulation_instruction(self, expr, visitor):
- if get_form_option('vectorization_blockstructured'):
- from dune.codegen.blockstructured.accumulation import generate_accumulation_instruction
- return generate_accumulation_instruction(expr, visitor)
- else:
- return PDELabInterface.generate_accumulation_instruction(self, expr, visitor)
- #
- # Local function space related generator functions
- #
- # TODO current way to squeeze subelem iname in, not really ideal
-
- def lfs_inames(self, element, restriction, number=None, context=''):
- return sub_element_inames() + lfs_inames(element, restriction, number, context)
-
- #
- # Test and trial function related generator functions
- #
-
- def pymbolic_basis(self, element, restriction, number, context=''):
- return pymbolic_basis(element, restriction, number, context)
-
- def pymbolic_reference_gradient(self, element, restriction, number, context=''):
- return pymbolic_reference_gradient(element, restriction, number, context)
-
- #
- # Geometry related generator functions
- #
-
- def pymbolic_spatial_coordinate(self):
- return to_global(pymbolic_quadrature_position())
-
- def pymbolic_jacobian_determinant(self):
- return pymbolic_jacobian_determinant()
-
- def pymbolic_jacobian_inverse(self, i, j, restriction):
- return pymbolic_jacobian_inverse(i, j, restriction)
-
- def pymbolic_facet_jacobian_determinant(self):
- return pymbolic_facet_jacobian_determinant()
diff --git a/python/dune/codegen/blockstructured/accumulation.py b/python/dune/codegen/blockstructured/accumulation.py
index 469e69b4bba506ed50d2800c8bdd98c71338ca03..9af9759d5e3426df12a10db39f94048e48006867 100644
--- a/python/dune/codegen/blockstructured/accumulation.py
+++ b/python/dune/codegen/blockstructured/accumulation.py
@@ -1,17 +1,27 @@
-from dune.codegen.generation import instruction
+from dune.codegen.blockstructured.tools import sub_element_inames
+from dune.codegen.generation import accumulation_mixin, instruction
from dune.codegen.loopy.target import dtype_floatingpoint
from dune.codegen.options import get_form_option
-from dune.codegen.pdelab.geometry import world_dimension
-from dune.codegen.pdelab.localoperator import determine_accumulation_space
+from dune.codegen.pdelab.geometry import world_dimension, name_intersection_geometry_wrapper
+from dune.codegen.pdelab.localoperator import determine_accumulation_space, GenericAccumulationMixin
from dune.codegen.pdelab.argument import name_accumulation_variable
from dune.codegen.pdelab.localoperator import boundary_predicates
-from dune.codegen.generation.loopy import function_mangler, globalarg
+from dune.codegen.generation.loopy import function_mangler, globalarg, temporary_variable
import loopy as lp
import pymbolic.primitives as prim
from loopy.match import Writes
+@accumulation_mixin("blockstructured")
+class BlockStructuredAccumulationMixin(GenericAccumulationMixin):
+ def generate_accumulation_instruction(self, expr):
+ if get_form_option('vectorization_blockstructured'):
+ return generate_accumulation_instruction_vectorized(expr, self)
+ else:
+ return generate_accumulation_instruction(expr, self)
+
+
def name_accumulation_alias(container, accumspace):
name = container + "_" + accumspace.lfs.name + "_alias"
name_tail = container + "_" + accumspace.lfs.name + "_alias_tail"
@@ -40,24 +50,90 @@ def residual_weight_mangler(knl, func, arg_dtypes):
return lp.CallMangleInfo(func, (lp.types.NumpyType(dtype_floatingpoint()),), ())
+def blockstructured_boundary_predicated(measure, subdomain_id):
+ predicates = []
+
+ if subdomain_id not in ['everywhere', 'otherwise']:
+ subelem_inames = sub_element_inames()
+
+ face_id = "face_id"
+ temporary_variable(face_id, shape=())
+ instruction(code="{} = {}.indexInInside();".format(face_id, name_intersection_geometry_wrapper()),
+ assignees=(face_id,))
+
+ def face_id_equals(id):
+ return prim.Comparison(prim.Variable(face_id), "==", id)
+
+ def iname_equals(iname, i):
+ return prim.Comparison(prim.Variable(iname), "==", i)
+
+ k = get_form_option("number_of_blocks")
+
+ if world_dimension() >= 2:
+ predicates.append(prim.If(face_id_equals(0), iname_equals(subelem_inames[0], 0), True))
+ predicates.append(prim.If(face_id_equals(1), iname_equals(subelem_inames[0], k - 1), True))
+ predicates.append(prim.If(face_id_equals(2), iname_equals(subelem_inames[1], 0), True))
+ predicates.append(prim.If(face_id_equals(3), iname_equals(subelem_inames[1], k - 1), True))
+ if world_dimension() == 3:
+ predicates.append(prim.If(face_id_equals(4), iname_equals(subelem_inames[2], 0), True))
+ predicates.append(prim.If(face_id_equals(5), iname_equals(subelem_inames[2], k - 1), True))
+ return frozenset(predicates)
+
+
def generate_accumulation_instruction(expr, visitor):
# Collect the lfs and lfs indices for the accumulate call
test_lfs = determine_accumulation_space(visitor.test_info, 0)
# In the jacobian case, also determine the space for the ansatz space
ansatz_lfs = determine_accumulation_space(visitor.trial_info, 1)
+ # Collect the lfs and lfs indices for the accumulate call
+ accumvar = name_accumulation_variable(test_lfs.get_restriction() + ansatz_lfs.get_restriction())
+
+ predicates = boundary_predicates(visitor.measure, visitor.subdomain_id)
+ predicates = predicates.union(blockstructured_boundary_predicated(visitor.measure, visitor.subdomain_id))
+
+ rank = 1 if ansatz_lfs.lfs is None else 2
+
+ from dune.codegen.pdelab.argument import PDELabAccumulationFunction
+ from pymbolic.primitives import Call
+ accexpr = Call(PDELabAccumulationFunction(accumvar, rank),
+ (test_lfs.get_args() + ansatz_lfs.get_args() + (expr,))
+ )
+
+ from dune.codegen.generation import instruction
+ quad_inames = visitor.quadrature_inames()
+ lfs_inames = frozenset(visitor.test_info.inames)
+ if visitor.trial_info:
+ lfs_inames = lfs_inames.union(visitor.trial_info.inames)
+
+ instruction(assignees=(),
+ expression=accexpr,
+ forced_iname_deps=lfs_inames.union(frozenset(quad_inames)),
+ forced_iname_deps_is_final=True,
+ predicates=predicates
+ )
+
+
+def generate_accumulation_instruction_vectorized(expr, visitor):
+ # Collect the lfs and lfs indices for the accumulate call
+ test_lfs = determine_accumulation_space(visitor.test_info, 0)
+ # In the jacobian case, also determine the space for the ansatz space
+ ansatz_lfs = determine_accumulation_space(visitor.trial_info, 1)
+
# Collect the lfs and lfs indices for the accumulate call
accumvar = name_accumulation_variable(test_lfs.get_restriction() + ansatz_lfs.get_restriction())
accumvar_alias = name_accumulation_alias(accumvar, test_lfs)
predicates = boundary_predicates(visitor.measure, visitor.subdomain_id)
+ predicates = predicates.union(blockstructured_boundary_predicated(visitor.measure, visitor.subdomain_id))
- quad_inames = visitor.interface.quadrature_inames()
+ quad_inames = visitor.quadrature_inames()
lfs_inames = visitor.test_info.inames
if visitor.trial_info:
lfs_inames = lfs_inames + visitor.trial_info.inames
- assignee = prim.Subscript(prim.Variable(accumvar_alias), tuple(prim.Variable(i) for i in lfs_inames))
+ assignee = prim.Subscript(prim.Variable(accumvar_alias),
+ tuple(prim.Variable(i) for i in sub_element_inames() + lfs_inames))
expr_with_weight = prim.Product((expr, prim.Call(prim.Variable(accumvar + '.weight'), ())))
diff --git a/python/dune/codegen/blockstructured/argument.py b/python/dune/codegen/blockstructured/argument.py
index afe120a8ec142820cf8ea3944538802b721f4469..420773e85bea93ee55cb310255da7fa60d55d9de 100644
--- a/python/dune/codegen/blockstructured/argument.py
+++ b/python/dune/codegen/blockstructured/argument.py
@@ -1,5 +1,4 @@
-from dune.codegen.generation import (backend,
- kernel_cached,
+from dune.codegen.generation import (kernel_cached,
valuearg, instruction, globalarg)
from dune.codegen.options import get_form_option
from dune.codegen.pdelab.argument import CoefficientAccess
@@ -27,7 +26,6 @@ def name_alias(container, lfs, element):
# TODO remove the need for element
@kernel_cached
-@backend(interface="pymbolic_coefficient", name="blockstructured")
def pymbolic_coefficient(container, lfs, element, index):
# TODO introduce a proper type for local function spaces!
if isinstance(lfs, str):
diff --git a/python/dune/codegen/blockstructured/basis.py b/python/dune/codegen/blockstructured/basis.py
index dcaf258b17e0dc4128860d3c58b9801a454f4471..5af8bb48f2bb625a0b336c5cae6d42863f0063c5 100644
--- a/python/dune/codegen/blockstructured/basis.py
+++ b/python/dune/codegen/blockstructured/basis.py
@@ -1,35 +1,142 @@
-from dune.codegen.generation import (backend,
+from loopy import Reduction
+
+from dune.codegen.generation import (basis_mixin,
kernel_cached,
- get_backend,
instruction,
temporary_variable,
globalarg,
class_member,
initializer_list,
include_file,)
-from dune.codegen.tools import get_pymbolic_basename
+from dune.codegen.tools import get_pymbolic_basename, get_pymbolic_indices
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)
from dune.codegen.pdelab.driver import (isPk,
isQk,
)
-from dune.codegen.pdelab.geometry import world_dimension
-from dune.codegen.pdelab.quadrature import pymbolic_quadrature_position_in_cell
-from dune.codegen.pdelab.spaces import type_leaf_gfs
+from dune.codegen.pdelab.geometry import world_dimension, component_iname
+from dune.codegen.pdelab.spaces import type_leaf_gfs, name_lfs
from dune.codegen.pdelab.restriction import restricted_name
from dune.codegen.blockstructured.spaces import lfs_inames
-from dune.codegen.blockstructured.tools import tensor_index_to_sequential_index
+from dune.codegen.blockstructured.tools import tensor_index_to_sequential_index, sub_element_inames
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_in_micro())
+ 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_in_micro())
+ 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)}),
+ )
+
+ @kernel_cached
+ def evaluate_coefficient_gradient(self, element, name, container, restriction, index):
+ sub_element = element
+ if isinstance(element, MixedElement):
+ sub_element = element.extract_component(index)[1]
+ from ufl import FiniteElement
+ assert isinstance(sub_element, FiniteElement)
+
+ temporary_variable(name, shape=(element.cell().geometric_dimension(),), managed=True)
+
+ dimindex = component_iname(count=0)
+
+ lfs = name_lfs(element, restriction, index)
+ basis = self.implement_reference_gradient(sub_element, restriction, 0, context='trialgrad')
+ basisindex = get_pymbolic_indices(basis)[:-1]
+ from dune.codegen.tools import maybe_wrap_subscript
+ basis = maybe_wrap_subscript(basis, prim.Variable(dimindex))
+
+ from dune.codegen.blockstructured.argument import pymbolic_coefficient
+ coeff = pymbolic_coefficient(container, lfs, sub_element, basisindex)
+
+ assignee = prim.Subscript(prim.Variable(name), (prim.Variable(dimindex),))
+ reduction_expr = prim.Product((coeff, basis))
+
+ instruction(expression=Reduction("sum", basisindex, reduction_expr, allow_simultaneous=True),
+ assignee=assignee,
+ forced_iname_deps=frozenset(self.quadrature_inames() + (dimindex,) + sub_element_inames()),
+ forced_iname_deps_is_final=True,
+ )
+
+ @kernel_cached
+ def evaluate_coefficient(self, element, name, container, restriction, index):
+ sub_element = element
+ if isinstance(element, MixedElement):
+ sub_element = element.extract_component(index)[1]
+
+ from ufl import FiniteElement
+ assert isinstance(sub_element, FiniteElement)
+
+ temporary_variable(name, shape=(), managed=True)
+
+ lfs = name_lfs(element, restriction, index)
+ basis = self.implement_basis(sub_element, restriction, 0, context='trial')
+ basisindex = get_pymbolic_indices(basis)[:-1]
+
+ from dune.codegen.blockstructured.argument import pymbolic_coefficient
+ coeff = pymbolic_coefficient(container, lfs, sub_element, basisindex)
+
+ assignee = prim.Variable(name)
+ reduction_expr = prim.Product((coeff, basis))
+
+ instruction(expression=Reduction("sum", basisindex, reduction_expr, allow_simultaneous=True),
+ assignee=assignee,
+ forced_iname_deps=frozenset(self.quadrature_inames() + sub_element_inames()),
+ forced_iname_deps_is_final=True,
+ )
+
+
# define FE basis explicitly in localoperator
-@backend(interface="typedef_localbasis", name="blockstructured")
@class_member(classtag="operator")
def typedef_localbasis(element, name):
df = "typename {}::Traits::GridView::ctype".format(type_leaf_gfs(element))
@@ -67,51 +174,3 @@ def name_localbasis(leaf_element):
globalarg(name)
define_localbasis(leaf_element, name)
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))
diff --git a/python/dune/codegen/blockstructured/geometry.py b/python/dune/codegen/blockstructured/geometry.py
index 4f1a71d3e3674ea768ef63842ca5aa4def21db61..3ad3a9026d8a5287e54e3950aef0aa199d08cc4e 100644
--- a/python/dune/codegen/blockstructured/geometry.py
+++ b/python/dune/codegen/blockstructured/geometry.py
@@ -1,22 +1,116 @@
-from dune.codegen.pdelab.restriction import restricted_name
-
-from dune.codegen.generation import (get_backend,
+from dune.codegen.generation import (geometry_mixin,
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)
-from dune.codegen.pdelab.geometry import (world_dimension,
+from dune.codegen.options import get_form_option
+from dune.codegen.pdelab.geometry import (AxiparallelGeometryMixin,
+ EquidistantGeometryMixin,
+ GenericPDELabGeometryMixin,
+ world_dimension,
local_dimension,
- name_facet_jacobian_determinant,
- name_element_corners,
- component_iname)
-from dune.codegen.blockstructured.tools import sub_element_inames
+ component_iname,
+ enforce_boundary_restriction,
+ restricted_name,
+ name_cell_geometry
+ )
+from dune.codegen.blockstructured.tools import (sub_element_inames,
+ name_point_in_macro,
+ )
+from dune.codegen.ufl.modified_terminals import Restriction
import pymbolic.primitives as prim
from loopy.match import Writes
+@geometry_mixin("blockstructured_multilinear")
+class BlockStructuredGeometryMixin(GenericPDELabGeometryMixin):
+ def spatial_coordinate(self, o):
+ if self.measure == 'cell':
+ return self.to_global(self.quadrature_position_in_macro())
+ else:
+ assert self.measure == 'exterior_facet' or self.measure == 'interior_facet'
+ micro_qp = self.to_cell(self.quadrature_position_in_micro())
+
+ macro_qp = prim.Variable(name_point_in_macro(micro_qp, self), )
+
+ return self.to_global(macro_qp)
+
+ def to_global(self, local):
+ assert isinstance(local, prim.Expression)
+ name = get_pymbolic_basename(local) + "_global"
+
+ # TODO need to assert somehow that local is of codim 0
+ compute_multilinear_to_global_transformation(name, local, self)
+ return prim.Variable(name)
+
+ def facet_jacobian_determinant(self, o):
+ return prim.Quotient(GenericPDELabGeometryMixin.facet_jacobian_determinant(self, o),
+ 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)
+
+ 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):
+ 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 to_global(self, local):
+ assert isinstance(local, prim.Expression)
+ name = get_pymbolic_basename(local) + "_global"
+
+ # TODO need to assert somehow that local is of codim 0
+ compute_axiparallel_to_global_transformation(name, local, self)
+ return prim.Variable(name)
+
+ def jacobian_determinant(self, o):
+ return prim.Quotient(AxiparallelGeometryMixin.jacobian_determinant(self, o),
+ prim.Power(get_form_option("number_of_blocks"), local_dimension()))
+
+ def jacobian_inverse(self, o):
+ return prim.Product((AxiparallelGeometryMixin.jacobian_inverse(self, o),
+ get_form_option("number_of_blocks")))
+
+
+@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()
@@ -63,8 +157,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()
@@ -103,23 +197,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:
@@ -137,67 +231,40 @@ 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
-# 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):
+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)]
@@ -244,76 +311,24 @@ 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)
- 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():
- 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
-def define_point_in_macro(name, point_in_micro):
- dim = local_dimension()
- if get_form_option('vectorization_blockstructured'):
- temporary_variable(name, shape=(dim,), managed=True)
- else:
- temporary_variable(name, shape=(dim,), shape_impl=('fv',))
-
- # point_macro = (point_micro + index_micro) / number_of_blocks
- # where index_micro = tensor index of the micro element
- subelem_inames = sub_element_inames()
- for i in range(dim):
- if isinstance(point_in_micro, prim.Subscript):
- expr = prim.Subscript(point_in_micro.aggregate, point_in_micro.index + (i,))
- else:
- expr = prim.Subscript(point_in_micro, (i,))
- expr = prim.Sum((expr, prim.Variable(subelem_inames[i]),))
- expr = prim.Quotient(expr, get_form_option('number_of_blocks'))
- # TODO relax within inames
- instruction(assignee=prim.Subscript(prim.Variable(name), (i,)),
- expression=expr,
- within_inames=frozenset(subelem_inames + get_backend(interface="quad_inames")()),
- tags=frozenset({subelem_inames[i]})
- )
-
-
-# TODO add subelem inames if this function gets called
-# TODO change input parameter to string
-def name_point_in_macro(point_in_micro):
- assert isinstance(point_in_micro, prim.Expression)
- name = get_pymbolic_basename(point_in_micro) + "_macro"
- define_point_in_macro(name, point_in_micro)
+ define_jacobian_inverse_transposed(name, 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()
@@ -344,13 +359,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()
@@ -367,27 +382,34 @@ 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):
- macro = name_point_in_macro(local)
- name = macro + "_global"
+def define_element_corners(name):
+ from dune.codegen.pdelab.driver import get_form
+ n_corners = get_form().ufl_cell().num_vertices()
+ temporary_variable(name, shape_impl=('fv', 'fv'), shape=(n_corners, world_dimension()))
- it = get_global_context_value("integral_type")
+ iname = "i_corner"
+ domain(iname, n_corners)
+ 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))
+ restriction = Restriction.NONE
+ elif it == 'exterior_facet':
+ restriction = Restriction.POSITIVE
else:
- raise NotImplementedError
+ raise NotImplementedError()
- return prim.Variable(name)
+ instruction(code="{0}[{1}] = {2}.corner({3});".format(name, iname, name_cell_geometry(restriction), iname),
+ assignees=frozenset({name}),
+ within_inames=frozenset({iname}), within_inames_is_final=True)
+
+
+def name_element_corners():
+ name = "corners"
+ define_element_corners(name)
+ return name
diff --git a/python/dune/codegen/blockstructured/quadrature.py b/python/dune/codegen/blockstructured/quadrature.py
index 4c78d733866b79817f1d7f43189e157d8967fa76..d3ba4950e8034232ebd41258f571b887caf198a4 100644
--- a/python/dune/codegen/blockstructured/quadrature.py
+++ b/python/dune/codegen/blockstructured/quadrature.py
@@ -1,22 +1,22 @@
-from dune.codegen.generation import (backend)
-from dune.codegen.pdelab.quadrature import (name_quadrature_points,
- quadrature_iname)
-from dune.codegen.blockstructured.geometry import name_point_in_macro
+from dune.codegen.error import CodegenError
+from dune.codegen.generation import quadrature_mixin
+from dune.codegen.pdelab.quadrature import GenericQuadratureMixin
+from dune.codegen.blockstructured.tools import name_point_in_macro
import pymbolic.primitives as prim
-@backend(interface="quad_pos", name='blockstructured')
-def pymbolic_quadrature_position():
- quad_points = name_quadrature_points()
- quad_iname = quadrature_iname()
+@quadrature_mixin("blockstructured")
+class BlockstructuredQuadratureMixin(GenericQuadratureMixin):
+ def quadrature_position_in_micro(self, index=None):
+ return GenericQuadratureMixin.quadrature_position(self, index)
- qp_in_micro = prim.Subscript(prim.Variable(quad_points), (prim.Variable(quad_iname),))
+ def quadrature_position_in_macro(self, index=None):
+ original = self.quadrature_position_in_micro(index)
+ qp = prim.Variable(name_point_in_macro(original, self), )
+ if index is not None:
+ return prim.Subscript(qp, (index,))
+ else:
+ return qp
- 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, index=None):
+ raise CodegenError('Decide if the quadrature point is in the macro or micro element')
diff --git a/python/dune/codegen/blockstructured/spaces.py b/python/dune/codegen/blockstructured/spaces.py
index 688b1984a7cac29922de43d2dae5ad35313cebc5..4602789a77a12bd342adbf4ed8bebc59cb989d0c 100644
--- a/python/dune/codegen/blockstructured/spaces.py
+++ b/python/dune/codegen/blockstructured/spaces.py
@@ -1,5 +1,4 @@
-from dune.codegen.generation import (backend,
- domain)
+from dune.codegen.generation import domain
from dune.codegen.pdelab.geometry import world_dimension
from dune.codegen.pdelab.spaces import name_leaf_lfs
diff --git a/python/dune/codegen/blockstructured/tools.py b/python/dune/codegen/blockstructured/tools.py
index 3eafcae5914f21538766960d75dbdb2fa2e728ad..802b819f23754aa0862e122f4127a3db18593117 100644
--- a/python/dune/codegen/blockstructured/tools.py
+++ b/python/dune/codegen/blockstructured/tools.py
@@ -1,16 +1,9 @@
-from dune.codegen.ufl.modified_terminals import Restriction
+from dune.codegen.tools import get_pymbolic_basename
from dune.codegen.generation import (iname,
domain,
- get_global_context_value,
temporary_variable,
instruction)
-from dune.codegen.pdelab.geometry import (local_dimension,
- world_dimension,
- name_localcenter,
- pymbolic_in_cell_coordinates)
-
-from dune.codegen.pdelab.quadrature import quadrature_inames
-from dune.codegen.generation.counter import get_counted_variable
+from dune.codegen.pdelab.geometry import world_dimension
from dune.codegen.options import get_form_option
import pymbolic.primitives as prim
@@ -20,7 +13,7 @@ import pymbolic.primitives as prim
@iname
def sub_element_inames():
name = "subel"
- dim = local_dimension()
+ dim = world_dimension()
dim_names = ["x", "y", "z"] + [str(i) for i in range(4, dim + 1)]
inames = tuple()
for i in range(dim):
@@ -29,72 +22,6 @@ def sub_element_inames():
return inames
-# define inames for boundary integration
-# In the case of boundary integrationsince we have only d-1 loops,
-# but we need always d inames to compute the macro index.
-# Therefore we must find out which iname is constant.
-# Since loo.py cannot handle if-else blocks very well, this whole
-# computation seems a bit cumbersome.
-def sub_facet_inames():
- subelem_inames = sub_element_inames()
-
- center = pymbolic_in_cell_coordinates(prim.Variable(name_localcenter()), Restriction.POSITIVE)
-
- # check if iname[index] must be constant or not
- def predicate(index):
- return prim.Comparison(prim.Subscript(center, index), "==", 0.5)
-
- def conditional_instruction(index):
- # instruction for not constant iname
- # special case for the third iname
- instruction(assignee=prim.Variable(inames[index]),
- expression=prim.Variable(subelem_inames[1 if index == 2 else 0]),
- within_inames=frozenset(subelem_inames).union(frozenset(quadrature_inames())),
- predicates=frozenset([predicate(index)])
- )
- # instruction for constant iname
- instruction(assignee=prim.Variable(inames[index]),
- expression=prim.Product(((k - 1), prim.Subscript(center, (index,)))),
- within_inames=frozenset(subelem_inames).union(frozenset(quadrature_inames())),
- predicates=frozenset([prim.LogicalNot(predicate(index))])
- )
-
- k = get_form_option("number_of_blocks")
-
- inames = ("x",)
- temporary_variable(inames[0])
- conditional_instruction(0)
- if world_dimension() < 3:
- inames = inames + ("y",)
- temporary_variable(inames[1])
- conditional_instruction(1)
- else:
- inames = inames + ("y",)
- temporary_variable(inames[1])
- # one additional condition is needed in 3d for the second iname
- instruction(assignee=prim.Variable(inames[1]),
- expression=prim.Variable(subelem_inames[0]),
- within_inames=frozenset(subelem_inames).union(frozenset(quadrature_inames())),
- predicates=frozenset([prim.LogicalAnd((predicate(1), predicate(2)))])
- )
- instruction(assignee=prim.Variable(inames[1]),
- expression=prim.Variable(subelem_inames[1]),
- within_inames=frozenset(subelem_inames).union(frozenset(quadrature_inames())),
- predicates=frozenset([prim.LogicalAnd((predicate(1), predicate(0)))])
- )
- instruction(assignee=prim.Variable(inames[1]),
- expression=prim.Product(((k - 1), prim.Subscript(center, (1,)))),
- within_inames=frozenset(subelem_inames).union(frozenset(quadrature_inames())),
- predicates=frozenset([prim.LogicalNot(predicate(1))])
- )
-
- inames = inames + ("z",)
- temporary_variable(inames[2])
- conditional_instruction(2)
-
- return inames
-
-
# compute sequential index for given tensor index, the same as index in base-k to base-10
def tensor_index_to_sequential_index(indices, k):
return prim.Sum(tuple(prim.Variable(index) * k ** i for i, index in enumerate(indices)))
@@ -107,13 +34,43 @@ def sequential_index_to_tensor_index(iname, k):
# compute index for higher order FEM for a given Q1 index
def micro_index_to_macro_index(element, inames):
- it = get_global_context_value("integral_type")
- if it == "cell":
- subelem_inames = sub_element_inames()
- elif it == "exterior_facet" or it == "interior_facet":
- subelem_inames = sub_facet_inames()
+ subelem_inames = sub_element_inames()
k = get_form_option("number_of_blocks")
p = element.degree()
return prim.Sum(tuple((p * prim.Variable(si) + prim.Variable(bi)) * (p * k + 1) ** i
for i, (si, bi) in enumerate(zip(subelem_inames, inames))))
+
+
+# translate a point in the micro element into macro coordinates
+def define_point_in_macro(name, point_in_micro, visitor):
+ # TODO this won't work for 2d mannifolds
+ dim = world_dimension()
+ if get_form_option('vectorization_blockstructured'):
+ temporary_variable(name, shape=(dim,), managed=True)
+ else:
+ temporary_variable(name, shape=(dim,), shape_impl=('fv',))
+
+ # point_macro = (point_micro + index_micro) / number_of_blocks
+ # where index_micro = tensor index of the micro element
+ subelem_inames = sub_element_inames()
+ for i in range(dim):
+ if isinstance(point_in_micro, prim.Subscript):
+ expr = prim.Subscript(point_in_micro.aggregate, point_in_micro.index + (i,))
+ else:
+ expr = prim.Subscript(point_in_micro, (i,))
+ expr = prim.Sum((expr, prim.Variable(subelem_inames[i]),))
+ expr = prim.Quotient(expr, get_form_option('number_of_blocks'))
+ # TODO relax within inames
+ instruction(assignee=prim.Subscript(prim.Variable(name), (i,)),
+ expression=expr,
+ within_inames=frozenset(subelem_inames + visitor.quadrature_inames()),
+ tags=frozenset({subelem_inames[i]})
+ )
+
+
+def name_point_in_macro(point_in_micro, visitor):
+ assert isinstance(point_in_micro, prim.Expression)
+ name = get_pymbolic_basename(point_in_micro) + "_macro"
+ define_point_in_macro(name, point_in_micro, visitor)
+ return name
diff --git a/python/dune/codegen/generation/__init__.py b/python/dune/codegen/generation/__init__.py
index 3103bac657648779921eb66c7170f5a55059e209..d0cf1d4dc8b6880db13bfdec4458d815ea8208c2 100644
--- a/python/dune/codegen/generation/__init__.py
+++ b/python/dune/codegen/generation/__init__.py
@@ -1,9 +1,3 @@
-from __future__ import absolute_import
-
-from dune.codegen.generation.backend import (backend,
- get_backend,
- )
-
from dune.codegen.generation.counter import (get_counter,
get_counted_variable,
)
@@ -58,3 +52,10 @@ from dune.codegen.generation.context import (cache_restoring,
global_context,
get_global_context_value,
)
+
+from dune.codegen.generation.mixins import (accumulation_mixin,
+ basis_mixin,
+ construct_from_mixins,
+ geometry_mixin,
+ quadrature_mixin,
+ )
diff --git a/python/dune/codegen/generation/backend.py b/python/dune/codegen/generation/backend.py
deleted file mode 100644
index a45f7ac8be21f3c358e3473a4bc9a51e4deeb9f2..0000000000000000000000000000000000000000
--- a/python/dune/codegen/generation/backend.py
+++ /dev/null
@@ -1,44 +0,0 @@
-from dune.codegen.generation.cache import _RegisteredFunction
-from dune.codegen.options import option_switch
-from pytools import Record
-
-
-_backend_mapping = {}
-
-
-class FuncProxy(Record):
- def __init__(self, interface, name, func):
- Record.__init__(self, interface=interface, name=name, func=func)
-
- def __call__(self, *args, **kwargs):
- return self.func(*args, **kwargs)
-
-
-def register_backend(interface, name, func):
- _backend_mapping.setdefault(interface, {})
- _backend_mapping[interface][name] = func
-
-
-def backend(interface=None, name='default'):
- assert interface
-
- def _dec(func):
- if not isinstance(func, _RegisteredFunction):
- # Allow order independence of the generator decorators
- # and the backend decorator by delaying the registration
- func = FuncProxy(interface, name, func)
-
- register_backend(interface, name, func)
-
- return func
-
- return _dec
-
-
-def get_backend(interface=None, selector=option_switch("sumfact"), **kwargs):
- assert interface and selector
-
- select = selector(**kwargs)
- assert select in _backend_mapping[interface], "Implementation '{}' for interface '{}' missing!".format(select, interface)
-
- return _backend_mapping[interface][select]
diff --git a/python/dune/codegen/generation/cache.py b/python/dune/codegen/generation/cache.py
index 223ab1bf6012bab7e19043a6163b936221735fe2..2e555ae99e73b5917d082d6c631d675c70cfc3de 100644
--- a/python/dune/codegen/generation/cache.py
+++ b/python/dune/codegen/generation/cache.py
@@ -8,9 +8,10 @@ from dune.codegen.generation.context import (get_global_context_value,
)
from dune.codegen.generation.counter import get_counter
from dune.codegen.options import get_option
+from pytools import ImmutableRecord
# Store a global list of generator functions
-_generators = []
+_generators = {}
def _freeze(data):
@@ -85,16 +86,6 @@ class _RegisteredFunction(object):
# Initialize the memoization cache
self._memoize_cache = {}
- # Register this generator function
- _generators.append(self)
-
- # Allow order independence of the backend and the generator decorators.
- # If backend was applied first, we resolve the issued FuncProxy object
- from dune.codegen.generation.backend import FuncProxy, register_backend
- if isinstance(self.func, FuncProxy):
- register_backend(self.func.interface, self.func.name, self)
- self.func = self.func.func
-
def _get_content(self, key):
return self._memoize_cache[key].value
@@ -180,11 +171,32 @@ def generator_factory(**factory_kwargs):
# Modify the kwargs according to the factorys kwargs
for k in factory_kwargs:
kwargs[k] = factory_kwargs[k]
+ # If there args, this function is used as a decorator, as in this example
+ #
+ # @decorator
+ # def foo():
+ # pass
+ #
+ # If there are no args, this is used as a decorator factory:
+ #
+ # @decorator(bar=42)
+ # def foo():
+ # pass
+ #
if args:
assert len(args) == 1
- return _RegisteredFunction(args[0], **kwargs)
+ key = args[0]
+
+ if hasattr(key, "__name__") and key.__name__ == '<lambda>':
+ key = str(kwargs)
+
+ funcobj = _generators.setdefault(key, _RegisteredFunction(args[0], **kwargs))
+ return lambda *a, **ka: funcobj(*a, **ka)
else:
- return lambda f: _RegisteredFunction(f, **kwargs)
+ def __dec(f):
+ funcobj = _generators.setdefault(f, _RegisteredFunction(f, **kwargs))
+ return lambda *a, **ka: funcobj(*a, **ka)
+ return __dec
if no_deco:
return _dec(lambda x: x)
@@ -241,13 +253,13 @@ def retrieve_cache_items(condition=True, make_generable=False):
return content
# First yield all those items that are not sorted
- for gen in filter(lambda g: not g.counted, _generators):
+ for gen in filter(lambda g: not g.counted, _generators.values()):
for item in _filter_cache_items(gen, condition).values():
yield as_generable(item.value)
# And now the sorted ones
counted_ones = []
- for gen in filter(lambda g: g.counted, _generators):
+ for gen in filter(lambda g: g.counted, _generators.values()):
counted_ones.extend(_filter_cache_items(gen, condition).values())
for item in sorted(counted_ones, key=lambda i: i.count):
@@ -264,12 +276,12 @@ def delete_cache_items(condition=True, keep=False):
if not keep:
condition = "not ({})".format(condition)
- for gen in _generators:
+ for gen in _generators.values():
gen._memoize_cache = _filter_cache_items(gen, condition)
def retrieve_cache_functions(condition="True"):
- return [g.func for g in _generators if eval(condition, _ConditionDict(g.item_tags))]
+ return [g.func for g in _generators.values() if eval(condition, _ConditionDict(g.item_tags))]
def inspect_generator(gen):
diff --git a/python/dune/codegen/generation/context.py b/python/dune/codegen/generation/context.py
index cf22ddd4ef4a0dd731c2d3f6cffd805b8b47248c..b90e46a5f6c4acd09b8200ce45b00f14a565beec 100644
--- a/python/dune/codegen/generation/context.py
+++ b/python/dune/codegen/generation/context.py
@@ -37,16 +37,17 @@ class _CacheRestoringContext(object):
def __enter__(self):
from dune.codegen.generation.cache import _generators as g
self.cache = {}
- for i in g:
- self.cache[i] = {}
- for k, v in i._memoize_cache.items():
- self.cache[i][k] = v
+ for original_func, cache_func in g.items():
+ self.cache[original_func] = {}
+ for k, v in cache_func._memoize_cache.items():
+ self.cache[original_func][k] = v
def __exit__(self, exc_type, exc_value, traceback):
+ from dune.codegen.generation.cache import _generators as g
for i, c in self.cache.items():
- i._memoize_cache = {}
+ g[i]._memoize_cache = {}
for k, v in c.items():
- i._memoize_cache[k] = v
+ g[i]._memoize_cache[k] = v
def cache_restoring():
diff --git a/python/dune/codegen/generation/mixins.py b/python/dune/codegen/generation/mixins.py
new file mode 100644
index 0000000000000000000000000000000000000000..06141dbe873f53d4c4644ee7863bb3fca9fad90b
--- /dev/null
+++ b/python/dune/codegen/generation/mixins.py
@@ -0,0 +1,28 @@
+""" The infrastructure for registered mixin classes """
+
+from functools import partial
+
+
+_mixin_registry = {}
+
+
+def mixin_base(mixintype, name):
+ _mixin_registry.setdefault(mixintype, {})
+
+ def _dec(cls):
+ _mixin_registry[mixintype][name] = cls
+ return cls
+
+ return _dec
+
+
+def construct_from_mixins(base=object, mixins=[], mixintype="geometry", name="GeometryInterface"):
+ mixins = tuple(_mixin_registry[mixintype][m] for m in mixins)
+ return type(name, mixins + (base,), {})
+
+
+# A list of specific mixins that we keep around explicitly
+geometry_mixin = partial(mixin_base, "geometry")
+quadrature_mixin = partial(mixin_base, "quadrature")
+basis_mixin = partial(mixin_base, "basis")
+accumulation_mixin = partial(mixin_base, "accumulation")
diff --git a/python/dune/codegen/loopy/transformations/vectorize_quad.py b/python/dune/codegen/loopy/transformations/vectorize_quad.py
index 2ff6361f445ce44e73c659aa34bf7e4e7f02646b..c3989d8011f74a662ed82d73551a337109c9a776 100644
--- a/python/dune/codegen/loopy/transformations/vectorize_quad.py
+++ b/python/dune/codegen/loopy/transformations/vectorize_quad.py
@@ -11,7 +11,6 @@ from dune.codegen.loopy.transformations.vectorview import (add_vector_view,
get_vector_view_name,
)
from dune.codegen.loopy.symbolic import substitute
-from dune.codegen.sumfact.quadrature import quadrature_inames
from dune.codegen.tools import get_pymbolic_basename, get_pymbolic_tag, ceildiv
from dune.codegen.options import get_option
@@ -324,8 +323,10 @@ def _vectorize_quadrature_loop(knl, inames, suffix):
def vectorize_quadrature_loop(knl):
# Loop over the quadrature loops that exist in the kernel.
- # This is implemented a bit hacky right now...
- for key, inames in quadrature_inames._memoize_cache.items():
+ # This is implemented very hacky right now...
+ from dune.codegen.generation.cache import _generators as _g
+ gen = list(filter(lambda i: hasattr(i[0], "__name__") and i[0].__name__ == "_quadrature_inames", _g.items()))[0][1]
+ for key, inames in gen._memoize_cache.items():
element = key[0][0]
if element is None:
suffix = ''
diff --git a/python/dune/codegen/options.py b/python/dune/codegen/options.py
index 40f2d1a520ca5506604d8f5ec98709d692b75fc3..bb6bbafaf09dd0b0ae7330b174dfa5805769d09f 100644
--- a/python/dune/codegen/options.py
+++ b/python/dune/codegen/options.py
@@ -98,7 +98,6 @@ class CodegenFormOptionsArray(ImmutableRecord):
generate_jacobian_apply = CodegenOption(default=False, helpstr="Wether jacobian_allpy_* methods should be generated.")
generate_residuals = CodegenOption(default=True, helpstr="Whether alpha_* methods should be generated.")
unroll_dimension_loops = CodegenOption(default=False, helpstr="whether loops over the geometric dimension should be unrolled")
- precompute_quadrature_info = CodegenOption(default=True, helpstr="compute quadrature points and weights in the constructor of the local operator")
blockstructured = CodegenOption(default=False, helpstr="Use block structure")
number_of_blocks = CodegenOption(default=1, helpstr="Number of sub blocks in one direction")
vectorization_blockstructured = CodegenOption(default=False, helpstr="Vectorize block structuring")
@@ -108,6 +107,10 @@ class CodegenFormOptionsArray(ImmutableRecord):
control_variable = CodegenOption(default=None, helpstr="Name of control variable in UFL file")
block_preconditioner_diagonal = CodegenOption(default=False, helpstr="Whether this operator should implement the diagonal part of a block preconditioner")
block_preconditioner_offdiagonal = CodegenOption(default=False, helpstr="Whether this operator should implement the off-diagonal part of a block preconditioner")
+ geometry_mixins = CodegenOption(default="generic", helpstr="A comma separated list of mixin identifiers to use for geometries. Currently implemented mixins: generic, axiparallel, equidistant, sumfact_multilinear, sumfact_axiparallel, sumfact_equidistant")
+ quadrature_mixins = CodegenOption(default="generic", helpstr="A comma separated list of mixin identifiers to use for quadrature. Currently implemented: generic, sumfact")
+ basis_mixins = CodegenOption(default="generic", helpstr="A comma separated list of mixin identifiers to use for basis function evaluation. Currently implemented: generic, sumfact")
+ accumulation_mixins = CodegenOption(default="generic", helpstr="A comma separated list of mixin identifiers to use for accumulation. Currently implemented: generic, sumfact, control, blockstructured")
enable_volume = CodegenOption(default=True, helpstr="Whether to assemble volume integrals")
enable_skeleton = CodegenOption(default=True, helpstr="Whether to assemble skeleton integrals")
enable_boundary = CodegenOption(default=True, helpstr="Whether to assemble boundary integrals")
@@ -175,7 +178,20 @@ def process_global_options(opt):
@memoize
def process_form_options(opt, form):
if opt.sumfact:
- opt = opt.copy(unroll_dimension_loops=True)
+ opt = opt.copy(unroll_dimension_loops=True,
+ quadrature_mixins="sumfact",
+ basis_mixins="sumfact",
+ accumulation_mixins="sumfact",
+ )
+
+ if opt.blockstructured:
+ opt = opt.copy(accumulation_mixins="blockstructured",
+ quadrature_mixins="blockstructured",
+ basis_mixins="blockstructured"
+ )
+
+ if opt.control:
+ opt = opt.copy(accumulation_mixins="control")
if opt.numerical_jacobian:
opt = opt.copy(generate_jacobians=False, generate_jacobian_apply=False)
diff --git a/python/dune/codegen/pdelab/__init__.py b/python/dune/codegen/pdelab/__init__.py
index affa0725c89539fa88e470ac143f5d76fdb8a18b..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391 100644
--- a/python/dune/codegen/pdelab/__init__.py
+++ b/python/dune/codegen/pdelab/__init__.py
@@ -1,158 +0,0 @@
-""" The pdelab specific parts of the code generation process """
-
-# Trigger some imports that are needed to have all backend implementations visible
-# to the selection mechanisms
-from dune.codegen.generation import (get_backend)
-from dune.codegen.options import option_switch
-from dune.codegen.pdelab.argument import (pymbolic_apply_function,
- pymbolic_apply_function_gradient,
- pymbolic_trialfunction,
- pymbolic_trialfunction_gradient,
- pymbolic_coefficient
- )
-from dune.codegen.pdelab.basis import (pymbolic_basis,
- pymbolic_reference_gradient,
- )
-from dune.codegen.pdelab.function import pymbolic_gridfunction
-from dune.codegen.pdelab.geometry import (component_iname,
- pymbolic_cell_volume,
- pymbolic_facet_area,
- pymbolic_facet_jacobian_determinant,
- pymbolic_jacobian_determinant,
- pymbolic_jacobian_inverse,
- pymbolic_unit_inner_normal,
- pymbolic_unit_outer_normal,
- to_global,
- )
-from dune.codegen.pdelab.index import (name_index,
- )
-from dune.codegen.pdelab.quadrature import (pymbolic_quadrature_weight,
- pymbolic_quadrature_position,
- quadrature_inames,
- )
-from dune.codegen.pdelab.spaces import (lfs_inames,
- )
-from dune.codegen.pdelab.tensors import (pymbolic_list_tensor,
- pymbolic_identity,
- pymbolic_matrix_inverse,
- )
-
-
-class PDELabInterface(object):
- def __init__(self):
- # The visitor instance will be registered by its init method
- self.visitor = None
-
- # Accumulation interfaces
- def get_accumulation_info(self, expr, visitor):
- from dune.codegen.pdelab.localoperator import get_accumulation_info
- return get_accumulation_info(expr, visitor)
-
- def list_accumulation_infos(self, expr, visitor):
- from dune.codegen.pdelab.localoperator import list_accumulation_infos
- return list_accumulation_infos(expr, visitor)
-
- def generate_accumulation_instruction(self, expr, visitor):
- from dune.codegen.pdelab.localoperator import generate_accumulation_instruction
- return generate_accumulation_instruction(expr, visitor)
-
- #
- # TODO: The following ones are actually entirely PDELab independent!
- # They should be placed elsewhere and be used directly in the visitor.
- #
-
- def component_iname(self, context=None, count=None):
- return component_iname(context=context, count=count)
-
- def name_index(self, ind):
- return name_index(ind)
-
- #
- # Local function space related generator functions
- #
-
- def lfs_inames(self, element, restriction, number=None, context=''):
- return lfs_inames(element, restriction, number, context)
-
- def initialize_function_spaces(self, expr, visitor):
- from dune.codegen.pdelab.spaces import initialize_function_spaces
- return initialize_function_spaces(expr, visitor)
-
- #
- # Test and trial function related generator functions
- #
-
- # TODO: should pymbolic_basis also pass the interface?
- def pymbolic_basis(self, element, restriction, number, context=''):
- return pymbolic_basis(element, restriction, number, context)
-
- # TODO: should pymbolic_reference_gradient also pass the interface?
- def pymbolic_reference_gradient(self, element, restriction, number, context=''):
- return pymbolic_reference_gradient(element, restriction, number, context)
-
- def pymbolic_trialfunction_gradient(self, element, restriction, index):
- return pymbolic_trialfunction_gradient(self.visitor, element, restriction, index)
-
- def pymbolic_apply_function_gradient(self, element, restriction, index):
- return pymbolic_apply_function_gradient(self.visitor, element, restriction, index)
-
- def pymbolic_trialfunction(self, element, restriction, index):
- return pymbolic_trialfunction(self.visitor, element, restriction, index)
-
- def pymbolic_apply_function(self, element, restriction, index):
- return pymbolic_apply_function(self.visitor, element, restriction, index)
-
- def pymbolic_gridfunction(self, coeff, restriction, grad):
- return pymbolic_gridfunction(coeff, restriction, grad)
-
- #
- # Tensor expression related generator functions
- #
-
- def pymbolic_list_tensor(self, o):
- return pymbolic_list_tensor(o, self.visitor)
-
- def pymbolic_identity(self, o):
- return pymbolic_identity(o)
-
- def pymbolic_matrix_inverse(self, o):
- return pymbolic_matrix_inverse(o, self.visitor)
-
- #
- # Geometry related generator functions
- #
-
- def pymbolic_spatial_coordinate(self):
- return to_global(pymbolic_quadrature_position())
-
- def pymbolic_facet_jacobian_determinant(self):
- return pymbolic_facet_jacobian_determinant()
-
- def pymbolic_jacobian_determinant(self):
- return pymbolic_jacobian_determinant()
-
- def pymbolic_jacobian_inverse(self, i, j, restriction):
- return pymbolic_jacobian_inverse(i, j, restriction)
-
- def pymbolic_unit_inner_normal(self):
- return pymbolic_unit_inner_normal()
-
- def pymbolic_unit_outer_normal(self):
- return pymbolic_unit_outer_normal()
-
- def pymbolic_cell_volume(self, restriction):
- return pymbolic_cell_volume(restriction)
-
- def pymbolic_facet_area(self):
- return pymbolic_facet_area()
-
- #
- # Quadrature related generator functions
- #
-
- def pymbolic_quadrature_weight(self):
- return pymbolic_quadrature_weight()
-
- # TODO Should this be part of interface or not?
- def quadrature_inames(self):
- return quadrature_inames()
diff --git a/python/dune/codegen/pdelab/adjoint.py b/python/dune/codegen/pdelab/adjoint.py
index 2681b00391d46532a37846ab2c550f98ff5a3d96..d1a986165671817f8de808dd3ea77d3e0536e3f7 100644
--- a/python/dune/codegen/pdelab/adjoint.py
+++ b/python/dune/codegen/pdelab/adjoint.py
@@ -8,7 +8,8 @@ from loopy.types import NumpyType
import pymbolic.primitives as prim
-from dune.codegen.generation import (class_member,
+from dune.codegen.generation import (accumulation_mixin,
+ class_member,
constructor_parameter,
function_mangler,
get_global_context_value,
@@ -20,10 +21,11 @@ from dune.codegen.generation import (class_member,
from dune.codegen.options import (get_form_option,
)
from dune.codegen.loopy.target import dtype_floatingpoint
-from dune.codegen.pdelab import PDELabInterface
from dune.codegen.pdelab.localoperator import (boundary_predicates,
determine_accumulation_space,
extract_kernel_from_cache,
+ GenericAccumulationMixin,
+ get_visitor,
)
@@ -61,7 +63,7 @@ def generate_accumulation_instruction(expr, visitor, accumulation_index, number_
expr = prim.Sum((assignee, expr))
from dune.codegen.generation import instruction
- quad_inames = visitor.interface.quadrature_inames()
+ quad_inames = visitor.quadrature_inames()
instruction(assignee=assignee,
expression=expr,
forced_iname_deps=frozenset(quad_inames),
@@ -69,44 +71,22 @@ def generate_accumulation_instruction(expr, visitor, accumulation_index, number_
)
-def list_accumulation_infos(expr, visitor):
- return ["control", ]
-
-
-class ControlInterface(PDELabInterface):
- """Interface for generating the control localoperator
-
- In this case we will not accumulate in the residual vector but use
- a class member representing dJdm instead.
-
- """
- def __init__(self, accumulation_index, number_of_controls):
- """Create ControlInterface
-
- Arguments:
- ----------
- accumulation_index: In which component of the dJdm should be accumulated.
- number_of_controls: Number of components of dJdm. Needed for creating the member variable.
- """
+@accumulation_mixin("control")
+class AdjointAccumulationMixin(GenericAccumulationMixin):
+ def set_adjoint_information(self, accumulation_index, number_of_controls):
self.accumulation_index = accumulation_index
self.number_of_controls = number_of_controls
- def list_accumulation_infos(self, expr, visitor):
- return list_accumulation_infos(expr, visitor)
+ def list_accumulation_infos(self, expr):
+ return ["control"]
- def generate_accumulation_instruction(self, expr, visitor):
+ def generate_accumulation_instruction(self, expr):
return generate_accumulation_instruction(expr,
- visitor,
+ self,
self.accumulation_index,
self.number_of_controls)
-def get_visitor(measure, subdomain_id, accumulation_index, number_of_controls):
- interface = ControlInterface(accumulation_index, number_of_controls)
- from dune.codegen.ufl.visitor import UFL2LoopyVisitor
- return UFL2LoopyVisitor(interface, measure, subdomain_id)
-
-
def visit_integral(integral, accumulation_index, number_of_controls):
integrand = integral.integrand()
measure = integral.integral_type()
@@ -114,7 +94,8 @@ def visit_integral(integral, accumulation_index, number_of_controls):
# The visitor needs to know about the current index and the number
# of controls in order to generate the accumulation instruction
- visitor = get_visitor(measure, subdomain_id, accumulation_index, number_of_controls)
+ visitor = get_visitor(measure, subdomain_id)
+ visitor.set_adjoint_information(accumulation_index, number_of_controls)
# Start the visiting process!
visitor.accumulate(integrand)
@@ -162,7 +143,6 @@ def generate_kernel(forms):
return knl
-# @backend(interface="generate_kernels_per_integral")
def control_generate_kernels_per_integral(forms):
"""For the control problem forms will have one form for every
measure. Every form will only contain integrals of one type.
diff --git a/python/dune/codegen/pdelab/argument.py b/python/dune/codegen/pdelab/argument.py
index 6bf800d59cc6ddeb683873793a4e72f9d0b11f0b..c3cc48298c1b2efb70c6c7716b29b221486044e1 100644
--- a/python/dune/codegen/pdelab/argument.py
+++ b/python/dune/codegen/pdelab/argument.py
@@ -11,13 +11,9 @@ from dune.codegen.generation import (domain,
valuearg,
get_global_context_value,
kernel_cached,
- backend
)
from dune.codegen.loopy.target import dtype_floatingpoint
from dune.codegen.pdelab.index import name_index
-from dune.codegen.pdelab.basis import (evaluate_coefficient,
- evaluate_coefficient_gradient,
- )
from dune.codegen.pdelab.spaces import (lfs_iname,
name_lfs_bound,
)
@@ -89,38 +85,6 @@ def accumulation_mangler(target, func, dtypes):
)
-def pymbolic_trialfunction_gradient(visitor, element, restriction, index):
- rawname = "gradu_{}".format(index)
- name = restricted_name(rawname, restriction)
- container = name_coefficientcontainer(restriction)
- evaluate_coefficient_gradient(visitor, element, name, container, restriction, index)
- return Variable(name)
-
-
-def pymbolic_trialfunction(visitor, element, restriction, index):
- rawname = "u_{}".format(index)
- name = restricted_name(rawname, restriction)
- container = name_coefficientcontainer(restriction)
- evaluate_coefficient(visitor, element, name, container, restriction, index)
- return Variable(name)
-
-
-def pymbolic_apply_function_gradient(visitor, element, restriction, index):
- rawname = "gradz_func_{}".format(index)
- name = restricted_name(rawname, restriction)
- container = name_applycontainer(restriction)
- evaluate_coefficient_gradient(visitor, element, name, container, restriction, index)
- return Variable(name)
-
-
-def pymbolic_apply_function(visitor, element, restriction, index):
- rawname = "z_func_{}".format(index)
- name = restricted_name(rawname, restriction)
- container = name_applycontainer(restriction)
- evaluate_coefficient(visitor, element, name, container, restriction, index)
- return Variable(name)
-
-
def name_coefficientcontainer(restriction):
name = restricted_name("x", restriction)
return name
@@ -132,7 +96,6 @@ def name_applycontainer(restriction):
@kernel_cached
-@backend(interface="pymbolic_coefficient")
def pymbolic_coefficient(container, lfs, index):
# TODO introduce a proper type for local function spaces!
if isinstance(lfs, str):
diff --git a/python/dune/codegen/pdelab/basis.py b/python/dune/codegen/pdelab/basis.py
index c65e71ebfc9e3a38023033e27dc65d445734edbf..03c069ac518f01efaec8304356fddcf34915668f 100644
--- a/python/dune/codegen/pdelab/basis.py
+++ b/python/dune/codegen/pdelab/basis.py
@@ -1,8 +1,7 @@
""" Generators for basis evaluations """
-from dune.codegen.generation import (backend,
+from dune.codegen.generation import (basis_mixin,
class_member,
- get_backend,
include_file,
instruction,
kernel_cached,
@@ -12,6 +11,9 @@ from dune.codegen.generation import (backend,
from dune.codegen.options import (option_switch,
get_form_option,
)
+from dune.codegen.pdelab.argument import (name_applycontainer,
+ name_coefficientcontainer,
+ )
from dune.codegen.pdelab.spaces import (lfs_iname,
lfs_inames,
name_leaf_lfs,
@@ -19,15 +21,15 @@ from dune.codegen.pdelab.spaces import (lfs_iname,
name_lfs_bound,
type_gfs,
type_leaf_gfs,
+ initialize_function_spaces,
)
from dune.codegen.pdelab.geometry import (component_iname,
world_dimension,
- name_jacobian_inverse_transposed,
- to_cell_coordinates,
name_cell,
)
from dune.codegen.pdelab.localoperator import (lop_template_ansatz_gfs,
lop_template_test_gfs,
+ name_gridfunction_member,
)
from dune.codegen.tools import (get_pymbolic_basename,
get_pymbolic_indices,
@@ -39,13 +41,241 @@ from dune.codegen.pdelab.driver import (isPk,
isDG)
from pymbolic.primitives import Product, Subscript, Variable
+import pymbolic.primitives as prim
from ufl import MixedElement
from loopy import Reduction
-@backend(interface="typedef_localbasis")
+@basis_mixin("base")
+class BasisMixinBase(object):
+ def initialize_function_spaces(self, expr):
+ pass
+
+ def lfs_inames(self, element, restriction, number):
+ raise NotImplementedError("Basis Mixins should implement local function space inames")
+
+ def implement_basis(self, element, restriction, number):
+ raise NotImplementedError("Basis Mixins should implement the basis")
+
+ def implement_reference_gradient(self, element, restriction, number, index):
+ raise NotImplementedError("Basis Mixins should implement the basis gradient")
+
+ def implement_trialfunction(self, element, restriction, index):
+ raise NotImplementedError("Basis Mixins should implement trial function evaluation")
+
+ def implement_trialfunction_gradient(self, element, restriction, index):
+ raise NotImplementedError("Basis Mixins should implement trial function gradient evaluation")
+
+ def implement_apply_function(self, element, restriction, index):
+ raise NotImplementedError("Basis Mixins should implement linearization point evaluation")
+
+ def implement_apply_function_gradient(self, element, restriction, index):
+ raise NotImplementedError("Basis Mixins should implement linearization point gradient evaluation")
+
+ def implement_grid_function(self, coeff, restriction, grad):
+ raise NotImplementedError("Basis Mixins should implement grid function evaluation")
+
+
+@basis_mixin("generic")
+class GenericBasisMixin(BasisMixinBase):
+ def initialize_function_spaces(self, expr):
+ return initialize_function_spaces(expr, self)
+
+ def lfs_inames(self, element, restriction, number, context=""):
+ return (lfs_iname(element, restriction, number, 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)
+ iname, = self.lfs_inames(element, restriction, number, context=context)
+
+ return prim.Subscript(prim.Variable(name), (prim.Variable(iname), 0))
+
+ @kernel_cached
+ def evaluate_basis(self, element, name, restriction):
+ lfs = name_leaf_lfs(element, restriction)
+ temporary_variable(name,
+ shape=(name_lfs_bound(lfs), 1),
+ decl_method=declare_cache_temporary(element, restriction, 'Function'),
+ )
+ cache = name_localbasis_cache(element)
+ qp = self.to_cell(self.quadrature_position())
+
+ instruction(inames=self.quadrature_inames(),
+ code='{} = {}.evaluateFunction({}, {}.finiteElement().localBasis());'.format(name,
+ cache,
+ str(qp),
+ lfs,
+ ),
+ 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)
+ iname, = self.lfs_inames(element, restriction, number, context=context)
+
+ return prim.Subscript(prim.Variable(name), (prim.Variable(iname), 0))
+
+ @kernel_cached
+ def evaluate_reference_gradient(self, element, name, restriction):
+ lfs = name_leaf_lfs(element, restriction)
+ temporary_variable(name,
+ shape=(name_lfs_bound(lfs), 1, world_dimension()),
+ decl_method=declare_cache_temporary(element, restriction, 'Jacobian'),
+ )
+ cache = name_localbasis_cache(element)
+ qp = self.to_cell(self.quadrature_position())
+ instruction(inames=self.quadrature_inames(),
+ code='{} = {}.evaluateJacobian({}, {}.finiteElement().localBasis());'.format(name,
+ cache,
+ str(qp),
+ lfs,
+ ),
+ assignees=frozenset({name}),
+ read_variables=frozenset({get_pymbolic_basename(qp)}),
+ )
+
+ def implement_trialfunction(self, element, restriction, index):
+ rawname = "u_{}".format(index)
+ name = restricted_name(rawname, restriction)
+ container = name_coefficientcontainer(restriction)
+ self.evaluate_coefficient(element, name, container, restriction, index)
+ return prim.Variable(name)
+
+ def implement_trialfunction_gradient(self, element, restriction, index):
+ rawname = "gradu_{}".format(index)
+ name = restricted_name(rawname, restriction)
+ container = name_coefficientcontainer(restriction)
+ self.evaluate_coefficient_gradient(element, name, container, restriction, index)
+ return prim.Variable(name)
+
+ def implement_apply_function(self, element, restriction, index):
+ rawname = "z_func_{}".format(index)
+ name = restricted_name(rawname, restriction)
+ container = name_applycontainer(restriction)
+ self.evaluate_coefficient(element, name, container, restriction, index)
+ return prim.Variable(name)
+
+ def implement_apply_function_gradient(self, element, restriction, index):
+ rawname = "gradz_func_{}".format(index)
+ name = restricted_name(rawname, restriction)
+ container = name_applycontainer(restriction)
+ self.evaluate_coefficient_gradient(element, name, container, restriction, index)
+ return prim.Variable(name)
+
+ @kernel_cached
+ def evaluate_coefficient(self, element, name, container, restriction, index):
+ sub_element = element
+ if isinstance(element, MixedElement):
+ sub_element = element.extract_component(index)[1]
+
+ from ufl import FiniteElement
+ assert isinstance(sub_element, FiniteElement)
+
+ temporary_variable(name, shape=(), managed=True)
+
+ lfs = name_lfs(element, restriction, index)
+ basis = self.implement_basis(sub_element, restriction, 0, context='trial')
+ basisindex = get_pymbolic_indices(basis)[:-1]
+
+ # TODO get rid ot this!
+ if get_form_option("blockstructured"):
+ from dune.codegen.blockstructured.argument import pymbolic_coefficient
+ coeff = pymbolic_coefficient(container, lfs, sub_element, basisindex)
+ else:
+ from dune.codegen.pdelab.argument import pymbolic_coefficient
+ coeff = pymbolic_coefficient(container, lfs, basisindex[0])
+
+ assignee = prim.Variable(name)
+ reduction_expr = prim.Product((coeff, basis))
+
+ instruction(expression=Reduction("sum", basisindex, reduction_expr, allow_simultaneous=True),
+ assignee=assignee,
+ forced_iname_deps=frozenset(self.quadrature_inames()),
+ forced_iname_deps_is_final=True,
+ )
+
+ @kernel_cached
+ def evaluate_coefficient_gradient(self, element, name, container, restriction, index):
+ sub_element = element
+ if isinstance(element, MixedElement):
+ sub_element = element.extract_component(index)[1]
+ from ufl import FiniteElement
+ assert isinstance(sub_element, FiniteElement)
+
+ temporary_variable(name, shape=(element.cell().geometric_dimension(),), managed=True)
+
+ dimindex = component_iname(count=0)
+
+ lfs = name_lfs(element, restriction, index)
+ basis = self.implement_reference_gradient(sub_element, restriction, 0, context='trialgrad')
+ basisindex = get_pymbolic_indices(basis)[:-1]
+ from dune.codegen.tools import maybe_wrap_subscript
+ basis = maybe_wrap_subscript(basis, prim.Variable(dimindex))
+
+ # TODO get rid of this
+ if get_form_option("blockstructured"):
+ from dune.codegen.blockstructured.argument import pymbolic_coefficient
+ coeff = pymbolic_coefficient(container, lfs, sub_element, basisindex)
+ else:
+ from dune.codegen.pdelab.argument import pymbolic_coefficient
+ coeff = pymbolic_coefficient(container, lfs, basisindex[0])
+
+ assignee = prim.Subscript(prim.Variable(name), (prim.Variable(dimindex),))
+ reduction_expr = prim.Product((coeff, basis))
+
+ instruction(expression=Reduction("sum", basisindex, reduction_expr, allow_simultaneous=True),
+ assignee=assignee,
+ forced_iname_deps=frozenset(self.quadrature_inames()).union(frozenset({dimindex})),
+ forced_iname_deps_is_final=True,
+ )
+
+ def implement_gridfunction(self, coeff, restriction, grad):
+ name = "coeff{}{}".format(coeff.count(), "_grad" if grad else "")
+ self.define_grid_function(name, coeff, restriction, grad)
+ return prim.Subscript(prim.Variable(name), (0,))
+
+ def define_grid_function(self, name, coeff, restriction, grad):
+ diffOrder = 1 if grad else 0
+
+ gridfunction = name_gridfunction_member(coeff, restriction, diffOrder)
+ bind_gridfunction_to_element(gridfunction, restriction)
+
+ temporary_variable(name,
+ shape=(1,) + (world_dimension(),) * diffOrder,
+ decl_method=declare_grid_function_range(gridfunction),
+ managed=False,
+ )
+
+ quadpos = self.to_cell(self.quadrature_position())
+ instruction(code="{} = {}({});".format(name, gridfunction, quadpos),
+ assignees=frozenset({name}),
+ within_inames=frozenset(self.quadrature_inames()),
+ within_inames_is_final=True,
+ )
+
+
+@preamble
+def bind_gridfunction_to_element(gf, restriction):
+ element = name_cell(restriction)
+ return "{}.bind({});".format(gf, element)
+
+
+def declare_grid_function_range(gridfunction):
+ def _decl(name, kernel, decl_info):
+ return "typename decltype({})::Range {};".format(gridfunction, name)
+
+ return _decl
+
+
@class_member(classtag="operator")
def typedef_localbasis(element, name):
basis_type = "{}::Traits::FiniteElementMap::Traits::FiniteElementType::Traits::LocalBasisType".format(type_leaf_gfs(element))
@@ -61,7 +291,14 @@ def type_localbasis(element):
name = "DG{}_LocalBasis".format(element._degree)
else:
raise NotImplementedError("Element type not known in code generation")
- get_backend("typedef_localbasis", selector=option_switch("blockstructured"))(element, name)
+
+ # TODO get rid of this
+ if get_form_option("blockstructured"):
+ from dune.codegen.blockstructured.basis import typedef_localbasis as bs_typedef_localbasis
+ bs_typedef_localbasis(element, name)
+ else:
+ typedef_localbasis(element, name)
+
return name
@@ -94,62 +331,6 @@ def declare_cache_temporary(element, restriction, which):
return decl
-@backend(interface="evaluate_basis")
-@kernel_cached
-def evaluate_basis(leaf_element, name, restriction):
- lfs = name_leaf_lfs(leaf_element, restriction)
- temporary_variable(name, shape=(name_lfs_bound(lfs), 1), decl_method=declare_cache_temporary(leaf_element, restriction, 'Function'))
- cache = name_localbasis_cache(leaf_element)
- qp = get_backend("qp_in_cell")(restriction)
- instruction(inames=get_backend("quad_inames")(),
- code='{} = {}.evaluateFunction({}, {}.finiteElement().localBasis());'.format(name,
- cache,
- str(qp),
- lfs,
- ),
- 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)
- iname, = lfs_inames(leaf_element, restriction, number, context=context)
-
- return Subscript(Variable(name), (Variable(iname), 0))
-
-
-@backend(interface="evaluate_grad")
-@kernel_cached
-def evaluate_reference_gradient(leaf_element, name, restriction):
- lfs = name_leaf_lfs(leaf_element, restriction)
- temporary_variable(name, shape=(name_lfs_bound(lfs), 1, world_dimension()), decl_method=declare_cache_temporary(leaf_element, restriction, 'Jacobian'))
- cache = name_localbasis_cache(leaf_element)
- qp = get_backend("qp_in_cell")(restriction)
- instruction(inames=get_backend("quad_inames")(),
- code='{} = {}.evaluateJacobian({}, {}.finiteElement().localBasis());'.format(name,
- cache,
- str(qp),
- lfs,
- ),
- 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)
- iname, = lfs_inames(leaf_element, restriction, number, context=context)
-
- return Subscript(Variable(name), (Variable(iname), 0))
-
-
def shape_as_pymbolic(shape):
def _shape_as_pymbolic(s):
if isinstance(s, str):
@@ -157,71 +338,3 @@ def shape_as_pymbolic(shape):
else:
return s
return tuple(_shape_as_pymbolic(s) for s in shape)
-
-
-@kernel_cached
-def evaluate_coefficient(visitor, element, name, container, restriction, index):
- sub_element = element
- if isinstance(element, MixedElement):
- sub_element = element.extract_component(index)[1]
-
- from ufl import FiniteElement
- assert isinstance(sub_element, FiniteElement)
-
- temporary_variable(name, shape=(), managed=True)
-
- lfs = name_lfs(element, restriction, index)
- basis = visitor.interface.pymbolic_basis(sub_element, restriction, 0, context='trial')
- basisindex = get_pymbolic_indices(basis)[:-1]
- if get_form_option("blockstructured"):
- from dune.codegen.blockstructured.argument import pymbolic_coefficient
- coeff = pymbolic_coefficient(container, lfs, sub_element, basisindex)
- else:
- from dune.codegen.pdelab.argument import pymbolic_coefficient
- coeff = pymbolic_coefficient(container, lfs, basisindex[0])
-
- assignee = Variable(name)
-
- reduction_expr = Product((coeff, basis))
-
- instruction(expression=Reduction("sum", basisindex, reduction_expr, allow_simultaneous=True),
- assignee=assignee,
- forced_iname_deps=frozenset(get_backend("quad_inames")()),
- forced_iname_deps_is_final=True,
- )
-
-
-@kernel_cached
-def evaluate_coefficient_gradient(visitor, element, name, container, restriction, index):
- sub_element = element
- if isinstance(element, MixedElement):
- sub_element = element.extract_component(index)[1]
- from ufl import FiniteElement
- assert isinstance(sub_element, FiniteElement)
-
- temporary_variable(name, shape=(element.cell().geometric_dimension(),), managed=True)
-
- dimindex = component_iname(count=0)
-
- lfs = name_lfs(element, restriction, index)
- basis = visitor.interface.pymbolic_reference_gradient(sub_element, restriction, 0, context='trialgrad')
- basisindex = get_pymbolic_indices(basis)[:-1]
- from dune.codegen.tools import maybe_wrap_subscript
- basis = maybe_wrap_subscript(basis, Variable(dimindex))
-
- if get_form_option("blockstructured"):
- from dune.codegen.blockstructured.argument import pymbolic_coefficient
- coeff = pymbolic_coefficient(container, lfs, sub_element, basisindex)
- else:
- from dune.codegen.pdelab.argument import pymbolic_coefficient
- coeff = pymbolic_coefficient(container, lfs, basisindex[0])
-
- assignee = Subscript(Variable(name), (Variable(dimindex),))
-
- reduction_expr = Product((coeff, basis))
-
- instruction(expression=Reduction("sum", basisindex, reduction_expr, allow_simultaneous=True),
- assignee=assignee,
- forced_iname_deps=frozenset(get_backend("quad_inames")()).union(frozenset({dimindex})),
- forced_iname_deps_is_final=True,
- )
diff --git a/python/dune/codegen/pdelab/driver/visitor.py b/python/dune/codegen/pdelab/driver/visitor.py
index 0c45871f5d940a62f85386e7aaccd41d210d1abb..6fab5574d1af6c97841535264677a511d96a727f 100644
--- a/python/dune/codegen/pdelab/driver/visitor.py
+++ b/python/dune/codegen/pdelab/driver/visitor.py
@@ -14,8 +14,7 @@ def set_lop_to_starting_time():
class DriverUFL2PymbolicVisitor(UFL2LoopyVisitor):
def __init__(self):
- from dune.codegen.pdelab import PDELabInterface
- UFL2LoopyVisitor.__init__(self, PDELabInterface(), "exterior_facet", {})
+ UFL2LoopyVisitor.__init__(self, "exterior_facet", {})
def __call__(self, expr):
self.preambles = []
diff --git a/python/dune/codegen/pdelab/function.py b/python/dune/codegen/pdelab/function.py
deleted file mode 100644
index 74e6b6b6ad08cdf8b6b9c1f0a1c5aeeae5f3f3c9..0000000000000000000000000000000000000000
--- a/python/dune/codegen/pdelab/function.py
+++ /dev/null
@@ -1,52 +0,0 @@
-from dune.codegen.generation import (get_backend,
- instruction,
- kernel_cached,
- preamble,
- temporary_variable,
- )
-from dune.codegen.pdelab.geometry import (name_cell,
- world_dimension,
- )
-from dune.codegen.pdelab.localoperator import name_gridfunction_member
-
-import pymbolic.primitives as prim
-
-
-@preamble
-def bind_gridfunction_to_element(gf, restriction):
- element = name_cell(restriction)
- return "{}.bind({});".format(gf, element)
-
-
-def declare_grid_function_range(gridfunction):
- def _decl(name, kernel, decl_info):
- return "typename decltype({})::Range {};".format(gridfunction, name)
-
- return _decl
-
-
-@kernel_cached
-def pymbolic_evaluate_gridfunction(name, coeff, restriction, grad):
- diffOrder = 1 if grad else 0
-
- gridfunction = name_gridfunction_member(coeff, restriction, diffOrder)
- bind_gridfunction_to_element(gridfunction, restriction)
-
- temporary_variable(name,
- shape=(1,) + (world_dimension(),) * diffOrder,
- decl_method=declare_grid_function_range(gridfunction),
- managed=False,
- )
-
- quadpos = get_backend(interface="qp_in_cell")(restriction)
- instruction(code="{} = {}({});".format(name, gridfunction, quadpos),
- assignees=frozenset({name}),
- within_inames=frozenset(get_backend(interface="quad_inames")()),
- within_inames_is_final=True,
- )
-
-
-def pymbolic_gridfunction(coeff, restriction, grad):
- name = "coeff{}{}".format(coeff.count(), "_grad" if grad else "")
- pymbolic_evaluate_gridfunction(name, coeff, restriction, grad)
- return prim.Subscript(prim.Variable(name), (0,))
diff --git a/python/dune/codegen/pdelab/geometry.py b/python/dune/codegen/pdelab/geometry.py
index 7d1e722ad8d92143d44747cbc2622b7eec0431b2..d37a315f2926ffe2f914589925eb77eee9352d5e 100644
--- a/python/dune/codegen/pdelab/geometry.py
+++ b/python/dune/codegen/pdelab/geometry.py
@@ -1,9 +1,8 @@
from dune.codegen.ufl.modified_terminals import Restriction
from dune.codegen.pdelab.restriction import restricted_name
-from dune.codegen.generation import (backend,
- class_member,
+from dune.codegen.generation import (class_member,
domain,
- get_backend,
+ geometry_mixin,
get_global_context_value,
globalarg,
iname,
@@ -17,11 +16,11 @@ from dune.codegen.options import (get_form_option,
option_switch,
)
from dune.codegen.loopy.target import dtype_floatingpoint, type_floatingpoint
-from dune.codegen.pdelab.quadrature import quadrature_preamble
+from dune.codegen.pdelab.quadrature import (quadrature_preamble,
+ )
from dune.codegen.tools import get_pymbolic_basename
from ufl.algorithms import MultiFunction
from pymbolic.primitives import Variable
-from pymbolic.primitives import Expression as PymbolicExpression
from loopy.match import Writes
@@ -29,6 +28,278 @@ import numpy as np
import pymbolic.primitives as prim
+@geometry_mixin("base")
+class GeometryMixinBase(object):
+ """
+ This mixin will be in by default by the mixin magic,
+ so it should only define an interface and throw exceptions.
+ """
+ def geo_unhandled(self, o):
+ raise NotImplementedError("Geometry Mixins should handle {}".format(type(o)))
+
+ def jacobian(self, o):
+ raise NotImplementedError("How did you get Jacobian into your form? We only support JacobianInverse right now. Report!")
+
+ spatial_coordinate = geo_unhandled
+ facet_normal = geo_unhandled
+ jacobian_determinant = geo_unhandled
+ jacobian_inverse = geo_unhandled
+ facet_jacobian_determinant = geo_unhandled
+ cell_volume = geo_unhandled
+ facet_area = geo_unhandled
+
+ def to_global(self, local):
+ raise NotImplementedError("Geometry Mixins should implement a to_global mapping")
+
+ def to_cell(self, local):
+ raise NotImplementedError("Geometry Mixins should implement a to_cell mapping")
+
+
+@geometry_mixin("generic")
+class GenericPDELabGeometryMixin(GeometryMixinBase):
+ def spatial_coordinate(self, o):
+ return self.to_global(self.quadrature_position())
+
+ def to_global(self, local):
+ assert isinstance(local, prim.Expression)
+ name = get_pymbolic_basename(local) + "_global"
+
+ temporary_variable(name, shape=(world_dimension(),), shape_impl=("fv",))
+ geo = name_geometry()
+ code = "{} = {}.global({});".format(name,
+ geo,
+ local,
+ )
+
+ quadrature_preamble(self,
+ code,
+ assignees=frozenset({name}),
+ read_variables=frozenset({get_pymbolic_basename(local)}),
+ depends_on=frozenset({Writes(get_pymbolic_basename(local))})
+ )
+
+ return prim.Variable(name)
+
+ def to_cell(self, local):
+ assert isinstance(local, prim.Expression)
+
+ restriction = enforce_boundary_restriction(self)
+ if restriction == Restriction.NONE:
+ return local
+
+ basename = get_pymbolic_basename(local)
+ name = "{}_in_{}side".format(basename, "in" if restriction is Restriction.POSITIVE else "out")
+ temporary_variable(name, shape=(world_dimension(),), shape_impl=("fv",))
+ geo = name_in_cell_geometry(restriction)
+ quadrature_preamble(self,
+ "{} = {}.global({});".format(name,
+ geo,
+ str(local),
+ ),
+ assignees=frozenset({name}),
+ read_variables=frozenset({get_pymbolic_basename(local)}),
+ depends_on=frozenset({Writes(get_pymbolic_basename(local))}),
+ )
+
+ return prim.Variable(name)
+
+ def facet_jacobian_determinant(self, o):
+ name = "fdetjac"
+ self.define_facet_jacobian_determinant(name)
+ return prim.Variable(name)
+
+ def define_facet_jacobian_determinant(self, name):
+ # NB: This can reuse the same implementation as the cell jacobian
+ # determinant, because the obtain geometry will be a face geometry
+ return self.define_jacobian_determinant(name)
+
+ def jacobian_determinant(self, o):
+ name = 'detjac'
+ self.define_jacobian_determinant(name)
+ return prim.Variable(name)
+
+ def define_jacobian_determinant(self, name):
+ temporary_variable(name, shape=())
+ geo = name_geometry()
+ pos = self.quadrature_position()
+ code = "{} = {}.integrationElement({});".format(name,
+ geo,
+ str(pos),
+ )
+
+ return quadrature_preamble(self,
+ code,
+ assignees=frozenset({name}),
+ read_variables=frozenset({get_pymbolic_basename(pos)}),
+ depends_on=frozenset({Writes(get_pymbolic_basename(pos))}),
+ )
+
+ def jacobian_inverse(self, o):
+ 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.Subscript(prim.Variable(name), (j, i))
+
+ def define_jacobian_inverse_transposed(self, name, restriction):
+ dim = world_dimension()
+ temporary_variable(name, decl_method=define_jacobian_inverse_transposed_temporary(restriction), shape=(dim, dim))
+ geo = name_cell_geometry(restriction)
+ pos = self.to_cell(self.quadrature_position())
+
+ return quadrature_preamble(self,
+ "{} = {}.jacobianInverseTransposed({});".format(name,
+ geo,
+ str(pos),
+ ),
+ assignees=frozenset({name}),
+ read_variables=frozenset({get_pymbolic_basename(pos)}),
+ depends_on=frozenset({Writes(get_pymbolic_basename(pos))}),
+ )
+
+ def facet_normal(self, o):
+ if self.restriction == Restriction.NEGATIVE:
+ raise NotImplementedError("Inner Normals not implemented!")
+
+ name = "unit_outer_normal"
+ self.define_unit_outer_normal(name)
+ return prim.Variable(name)
+
+ def define_unit_outer_normal(self, name):
+ temporary_variable(name, shape=(world_dimension(),), decl_method=declare_normal)
+
+ ig = name_intersection_geometry_wrapper()
+ qp = self.quadrature_position()
+ return quadrature_preamble(self,
+ "{} = {}.unitOuterNormal({});".format(name, ig, qp),
+ assignees=frozenset({name}),
+ read_variables=frozenset({get_pymbolic_basename(qp)}),
+ depends_on=frozenset({Writes(get_pymbolic_basename(qp))}),
+ )
+
+ def cell_volume(self, o):
+ restriction = enforce_boundary_restriction(self)
+ name = restricted_name("volume", restriction)
+ self.define_cell_volume(name, restriction)
+ return prim.Variable(name)
+
+ @preamble
+ def define_cell_volume(self, name, restriction):
+ geo = name_cell_geometry(restriction)
+ valuearg(name)
+ return "auto {} = {}.volume();".format(name, geo)
+
+ def facet_area(self, o):
+ name = "area"
+ self.define_facet_area(name)
+ return prim.Variable(name)
+
+ @preamble
+ def define_facet_area(self, name):
+ geo = name_intersection_geometry()
+ valuearg(name)
+ return "auto {} = {}.volume();".format(name, geo)
+
+
+@geometry_mixin("axiparallel")
+class AxiparallelGeometryMixin(GenericPDELabGeometryMixin):
+ def define_unit_outer_normal(self, name):
+ # NB: Strictly speaking, this implementation holds for any non-curved intersection.
+ declare_normal(name, None, None)
+ globalarg(name, shape=(world_dimension(),))
+
+ def jacobian_inverse(self, o):
+ i, j = self.indices
+ assert isinstance(i, int) and isinstance(j, int)
+
+ if i != j:
+ self.indices = None
+ return 0
+
+ jac = GenericPDELabGeometryMixin.jacobian_inverse(self, o)
+
+ return jac
+
+ def facet_area(self, o):
+ # This is not 100% correct, but in practice structured grid implementations are not
+ # embedded into higher dimensional world space.
+ return self.facet_jacobian_determinant(o)
+
+ def cell_volume(self, o):
+ # This is not 100% correct, but in practice structured grid implementations are not
+ # embedded into higher dimensional world space.
+ return self.jacobian_determinant(o)
+
+ def define_facet_jacobian_determinant(self, name):
+ valuearg(name)
+ self._define_facet_jacobian_determinant(name)
+
+ @preamble
+ def _define_facet_jacobian_determinant(self, name):
+ # NB: In the equidistant case, this might be optimized to store *d* values on the operator level
+ # We don't do that right now for laziness.
+ geo = name_geometry()
+ pos = name_localcenter()
+
+ return "auto {} = {}.integrationElement({});".format(name,
+ geo,
+ pos,
+ )
+
+
+@geometry_mixin("equidistant")
+class EquidistantGeometryMixin(AxiparallelGeometryMixin):
+ def define_jacobian_determinant(self, name):
+ valuearg(name)
+ self._define_jacobian_determinant(name)
+
+ @class_member(classtag="operator")
+ def _define_jacobian_determinant(self, name):
+ from dune.codegen.pdelab.localoperator import lop_template_ansatz_gfs
+ gfst = lop_template_ansatz_gfs()
+ self._define_jacobian_determinant_eval(name)
+ return "typename {}::Traits::GridView::template Codim<0>::Geometry::ctype {};".format(gfst, name)
+
+ @preamble(kernel="operator")
+ def _define_jacobian_determinant_eval(self, name):
+ from dune.codegen.pdelab.localoperator import name_ansatz_gfs_constructor_param, lop_template_range_field
+ gfs = name_ansatz_gfs_constructor_param()
+ rft = lop_template_range_field()
+ return "{} = {}.gridView().template begin<0>()->geometry().integrationElement(Dune::FieldVector<{}, {}>());".format(name, gfs, rft, world_dimension())
+
+ def define_jacobian_inverse_transposed(self, name, restriction):
+ dim = world_dimension()
+ globalarg(name, shape=(dim, dim), managed=False)
+ self._define_jacobian_inverse_transposed(name, restriction)
+
+ @class_member(classtag="operator")
+ def _define_jacobian_inverse_transposed(self, name, restriction):
+ dim = world_dimension()
+ self._define_jacobian_inverse_transposed_eval(name)
+ from dune.codegen.pdelab.localoperator import lop_template_ansatz_gfs
+ gfst = lop_template_ansatz_gfs()
+ return "typename {}::Traits::GridView::template Codim<0>::Geometry::JacobianInverseTransposed {};".format(gfst, name)
+
+ @preamble(kernel="operator")
+ def _define_jacobian_inverse_transposed_eval(self, name):
+ from dune.codegen.pdelab.localoperator import name_ansatz_gfs_constructor_param, lop_template_range_field
+ gfs = name_ansatz_gfs_constructor_param()
+ rft = lop_template_range_field()
+ return "{} = {}.gridView().template begin<0>()->geometry().jacobianInverseTransposed(Dune::FieldVector<{}, {}>());".format(name, gfs, rft, world_dimension())
+
+
+def enforce_boundary_restriction(visitor):
+ if visitor.measure == 'exterior_facet' and visitor.restriction == Restriction.NONE:
+ return Restriction.POSITIVE
+ else:
+ return visitor.restriction
+
+
@preamble
def define_reference_element(name):
geo = name_geometry()
@@ -201,35 +472,6 @@ def name_in_cell_geometry(restriction):
return name
-# TODO is it always necessary to add quadrature inames?
-def apply_in_cell_transformation(name, local, restriction):
- geo = name_in_cell_geometry(restriction)
- return quadrature_preamble("{} = {}.global({});".format(name,
- geo,
- str(local),
- ),
- assignees=frozenset({name}),
- read_variables=frozenset({get_pymbolic_basename(local)}),
- depends_on=frozenset({Writes(get_pymbolic_basename(local))}),
- )
-
-
-def pymbolic_in_cell_coordinates(local, restriction):
- basename = get_pymbolic_basename(local)
- name = "{}_in_{}side".format(basename, "in" if restriction is Restriction.POSITIVE else "out")
- temporary_variable(name, shape=(world_dimension(),), shape_impl=("fv",))
- apply_in_cell_transformation(name, local, restriction)
- return Variable(name)
-
-
-def to_cell_coordinates(local, restriction):
- assert isinstance(local, PymbolicExpression)
- if restriction == Restriction.NONE:
- return local
- else:
- return pymbolic_in_cell_coordinates(local, restriction)
-
-
def world_dimension():
data = get_global_context_value("data")
form = data.object_by_name[get_form_option("form")]
@@ -250,48 +492,12 @@ def local_dimension():
return intersection_dimension()
-def evaluate_unit_outer_normal(name):
- ig = name_intersection_geometry_wrapper()
- qp = get_backend("quad_pos")()
- return quadrature_preamble("{} = {}.unitOuterNormal({});".format(name, ig, qp),
- assignees=frozenset({name}),
- read_variables=frozenset({get_pymbolic_basename(qp)}),
- depends_on=frozenset({Writes(get_pymbolic_basename(qp))}),
- )
-
-
@preamble
def declare_normal(name, kernel, decl_info):
ig = name_intersection_geometry_wrapper()
return "auto {} = {}.centerUnitOuterNormal();".format(name, ig)
-def pymbolic_unit_outer_normal():
- name = "unit_outer_normal"
- if not get_form_option("diagonal_transformation_matrix"):
- temporary_variable(name, shape=(world_dimension(),), decl_method=declare_normal)
- evaluate_unit_outer_normal(name)
- else:
- declare_normal(name, None, None)
- globalarg(name, shape=(world_dimension(),))
- return prim.Variable(name)
-
-
-def evaluate_unit_inner_normal(name):
- outer = name_unit_outer_normal()
- return quadrature_preamble("auto {} = -1. * {};".format(name, outer),
- assignees=frozenset({name}),
- read_variables=frozenset({outer}),
- )
-
-
-def pymbolic_unit_inner_normal():
- name = "inner_normal"
- temporary_variable(name, shape=(world_dimension(),), decl_method=declare_normal)
- evaluate_unit_inner_normal(name)
- return prim.Variable(name)
-
-
def type_jacobian_inverse_transposed(restriction):
geo = type_cell_geometry(restriction)
return "typename {}::JacobianInverseTransposed".format(geo)
@@ -306,213 +512,3 @@ def define_jacobian_inverse_transposed_temporary(restriction):
name,
)
return _define_jacobian_inverse_transposed_temporary
-
-
-@preamble(kernel="operator")
-def define_constant_jacobian_inverse_transposed_eval(name):
- from dune.codegen.pdelab.localoperator import name_ansatz_gfs_constructor_param, lop_template_range_field
- gfs = name_ansatz_gfs_constructor_param()
- rft = lop_template_range_field()
-
- return "{} = {}.gridView().template begin<0>()->geometry().jacobianInverseTransposed(Dune::FieldVector<{}, {}>());".format(name, gfs, rft, world_dimension())
-
-
-@class_member(classtag="operator")
-def define_constant_jacobian_inverse_transposed(name):
- define_constant_jacobian_inverse_transposed_eval(name)
- from dune.codegen.pdelab.localoperator import lop_template_ansatz_gfs
- gfst = lop_template_ansatz_gfs()
- return "typename {}::Traits::GridView::template Codim<0>::Geometry::JacobianInverseTransposed {};".format(gfst, name)
-
-
-@backend(interface="name_jit", name="constant_transformation_matrix")
-def name_constant_jacobian_inverse_transposed(restriction):
- name = "jit"
- dim = world_dimension()
- globalarg(name, shape=(dim, dim), managed=False)
- define_constant_jacobian_inverse_transposed(name)
- return name
-
-
-def define_jacobian_inverse_transposed(name, restriction):
- dim = world_dimension()
- temporary_variable(name, decl_method=define_jacobian_inverse_transposed_temporary(restriction), shape=(dim, dim))
- geo = name_cell_geometry(restriction)
- pos = get_backend("qp_in_cell", selector=option_switch("blockstructured"))(restriction)
-
- return quadrature_preamble("{} = {}.jacobianInverseTransposed({});".format(name,
- geo,
- str(pos),
- ),
- assignees=frozenset({name}),
- read_variables=frozenset({get_pymbolic_basename(pos)}),
- depends_on=frozenset({Writes(get_pymbolic_basename(pos))}),
- )
-
-
-@backend(interface="name_jit", name="default")
-def name_jacobian_inverse_transposed(restriction):
- name = restricted_name("jit", restriction)
- define_jacobian_inverse_transposed(name, restriction)
- return name
-
-
-def pymbolic_jacobian_inverse(i, j, restriction):
- # Calculate jacobian_inverse_transposed
- name = get_backend(interface="name_jit", selector=option_switch("constant_transformation_matrix"))(restriction)
-
- # Return jacobian inverse -> (j, i) instead of (i, j)
- return prim.Subscript(prim.Variable(name), (j, i))
-
-
-@preamble(kernel="operator")
-def define_constant_jacobian_determinant_eval(name):
- from dune.codegen.pdelab.localoperator import name_ansatz_gfs_constructor_param, lop_template_range_field
- gfs = name_ansatz_gfs_constructor_param()
- rft = lop_template_range_field()
-
- return "{} = {}.gridView().template begin<0>()->geometry().integrationElement(Dune::FieldVector<{}, {}>());".format(name, gfs, rft, world_dimension())
-
-
-@class_member(classtag="operator")
-def _define_constant_jacobian_determinant(name):
- define_constant_jacobian_determinant_eval(name)
- from dune.codegen.pdelab.localoperator import lop_template_ansatz_gfs
- gfst = lop_template_ansatz_gfs()
- return "typename {}::Traits::GridView::template Codim<0>::Geometry::ctype {};".format(gfst, name)
-
-
-@backend(interface="detjac", name="constant_transformation_matrix")
-def define_constant_jacobian_determinant(name):
- valuearg(name)
- _define_constant_jacobian_determinant(name)
-
-
-@backend(interface="detjac", name="default")
-def define_jacobian_determinant(name):
- temporary_variable(name, shape=())
- geo = name_geometry()
- pos = get_backend("quad_pos", selector=option_switch("blockstructured"))()
- code = "{} = {}.integrationElement({});".format(name,
- geo,
- str(pos),
- )
-
- return quadrature_preamble(code,
- assignees=frozenset({name}),
- read_variables=frozenset({get_pymbolic_basename(pos)}),
- depends_on=frozenset({Writes(get_pymbolic_basename(pos))}),
- )
-
-
-@backend(interface="fdetjac", name="constant_transformation_matrix")
-@preamble
-def define_facet_jacobian_determinant(name):
- # NB: This might be optimized to store *d* values on the operator level
- # We don't do that right now for laziness.
- geo = name_geometry()
- pos = name_localcenter()
-
- valuearg(name)
-
- return "auto {} = {}.integrationElement({});".format(name,
- geo,
- pos,
- )
-
-
-@backend(interface="fdetjac", name="default")
-def define_facet_jacobian_determinant(name):
- return define_jacobian_determinant(name)
-
-
-def name_jacobian_determinant():
- name = 'detjac'
- get_backend(interface="detjac", selector=option_switch("constant_transformation_matrix"))(name)
- return name
-
-
-def pymbolic_jacobian_determinant():
- return prim.Variable(name_jacobian_determinant())
-
-
-def name_facet_jacobian_determinant():
- name = 'fdetjac'
- get_backend(interface="fdetjac", selector=option_switch("constant_transformation_matrix"))(name)
- return name
-
-
-def pymbolic_facet_jacobian_determinant():
- return prim.Variable(name_facet_jacobian_determinant())
-
-
-def apply_to_global_transformation(name, local):
- temporary_variable(name, shape=(world_dimension(),), shape_impl=("fv",))
- geo = name_geometry()
- code = "{} = {}.global({});".format(name,
- geo,
- local,
- )
- return quadrature_preamble(code,
- assignees=frozenset({name}),
- read_variables=frozenset({get_pymbolic_basename(local)}),
- depends_on=frozenset({Writes(get_pymbolic_basename(local))})
- )
-
-
-def to_global(local):
- assert isinstance(local, prim.Expression)
- name = get_pymbolic_basename(local) + "_global"
- apply_to_global_transformation(name, local)
- return prim.Variable(name)
-
-
-@preamble
-def define_cell_volume(name, restriction):
- geo = name_cell_geometry(restriction)
- valuearg(name)
- return "auto {} = {}.volume();".format(name, geo)
-
-
-def pymbolic_cell_volume(restriction):
- if get_form_option("constant_transformation_matrix"):
- return pymbolic_jacobian_determinant()
- else:
- name = restricted_name("volume", restriction)
- define_cell_volume(name, restriction)
- return prim.Variable(name)
-
-
-@preamble
-def define_facet_area(name):
- geo = name_intersection_geometry()
- valuearg(name)
- return "auto {} = {}.volume();".format(name, geo)
-
-
-def pymbolic_facet_area():
- if get_form_option("constant_transformation_matrix"):
- return pymbolic_facet_jacobian_determinant()
- else:
- name = "area"
- define_facet_area(name)
- 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()
- temporary_variable(name, shape_impl=('fv', 'fv'), shape=(n_corners, local_dimension()))
-
- iname = "i_corner"
- domain(iname, n_corners)
- from dune.codegen.generation import instruction
- instruction(code="{0}[{1}] = {2}.corner({3});".format(name, iname, name_geometry(), iname),
- assignees=frozenset({name}),
- within_inames=frozenset({iname}), within_inames_is_final=True)
-
-
-def name_element_corners():
- name = "corners"
- define_element_corners(name)
- return name
diff --git a/python/dune/codegen/pdelab/localoperator.py b/python/dune/codegen/pdelab/localoperator.py
index 12793f96977e46ccb252285454f11344810892bc..70ab9f4793a6163559d070e33e3600f9cd8c5578 100644
--- a/python/dune/codegen/pdelab/localoperator.py
+++ b/python/dune/codegen/pdelab/localoperator.py
@@ -9,17 +9,17 @@ from dune.codegen.options import (get_form_option,
get_option,
option_switch,
set_form_option)
-from dune.codegen.generation import (backend,
+from dune.codegen.generation import (accumulation_mixin,
base_class,
class_basename,
class_member,
+ construct_from_mixins,
constructor_parameter,
domain,
dump_accumulate_timer,
end_of_file,
function_mangler,
generator_factory,
- get_backend,
get_global_context_value,
global_context,
iname,
@@ -328,6 +328,30 @@ def boundary_predicates(measure, subdomain_id):
return frozenset(predicates)
+@accumulation_mixin("base")
+class AccumulationMixinBase(object):
+ def get_accumulation_info(self, expr):
+ raise NotImplementedError
+
+ def list_accumulation_infos(self, expr):
+ raise NotImplementedError
+
+ def generate_accumulation_instruction(self, expr):
+ raise NotImplementedError
+
+
+@accumulation_mixin("generic")
+class GenericAccumulationMixin(AccumulationMixinBase):
+ def get_accumulation_info(self, expr):
+ return get_accumulation_info(expr, self)
+
+ def list_accumulation_infos(self, expr):
+ return list_accumulation_infos(expr, self)
+
+ def generate_accumulation_instruction(self, expr):
+ return generate_accumulation_instruction(expr, self)
+
+
class PDELabAccumulationInfo(ImmutableRecord):
def __init__(self,
element=None,
@@ -390,10 +414,10 @@ def get_accumulation_info(expr, visitor):
if visitor.measure == 'exterior_facet':
restriction = Restriction.POSITIVE
- inames = visitor.interface.lfs_inames(leaf_element,
- restriction,
- expr.number()
- )
+ inames = visitor.lfs_inames(leaf_element,
+ restriction,
+ expr.number()
+ )
return PDELabAccumulationInfo(element=expr.ufl_element(),
element_index=element_index,
@@ -424,7 +448,7 @@ def generate_accumulation_instruction(expr, visitor):
from dune.codegen.generation import instruction
from dune.codegen.options import option_switch
- quad_inames = visitor.interface.quadrature_inames()
+ quad_inames = visitor.quadrature_inames()
lfs_inames = frozenset(visitor.test_info.inames)
if visitor.trial_info:
lfs_inames = lfs_inames.union(visitor.trial_info.inames)
@@ -438,19 +462,24 @@ def generate_accumulation_instruction(expr, visitor):
def get_visitor(measure, subdomain_id):
- # Get a transformer instance for this kernel
- if get_form_option('sumfact'):
- from dune.codegen.sumfact import SumFactInterface
- interface = SumFactInterface()
- elif get_form_option('blockstructured'):
- from dune.codegen.blockstructured import BlockStructuredInterface
- interface = BlockStructuredInterface()
- else:
- from dune.codegen.pdelab import PDELabInterface
- interface = PDELabInterface()
-
+ # Construct the visitor taking into account geometry mixins
from dune.codegen.ufl.visitor import UFL2LoopyVisitor
- return UFL2LoopyVisitor(interface, measure, subdomain_id)
+ mixins = get_form_option("geometry_mixins").split(",")
+ VisitorType = construct_from_mixins(base=UFL2LoopyVisitor, mixins=mixins, mixintype="geometry", name="UFLVisitor")
+
+ # Mix quadrature mixins in
+ mixins = get_form_option("quadrature_mixins").split(",")
+ VisitorType = construct_from_mixins(base=VisitorType, mixins=mixins, mixintype="quadrature", name="UFLVisitor")
+
+ # Mix basis mixins in
+ mixins = get_form_option("basis_mixins").split(",")
+ VisitorType = construct_from_mixins(base=VisitorType, mixins=mixins, mixintype="basis", name="UFLVisitor")
+
+ # Mix accumulation mixins in
+ mixins = get_form_option("accumulation_mixins").split(",")
+ VisitorType = construct_from_mixins(base=VisitorType, mixins=mixins, mixintype="accumulation", name="UFLVisitor")
+
+ return VisitorType(measure, subdomain_id)
def visit_integral(integral):
@@ -514,9 +543,13 @@ def generate_kernel(integrals):
return knl
-@backend(interface="generate_kernels_per_integral")
def generate_kernels_per_integral(integrals):
- yield generate_kernel(integrals)
+ if get_form_option("sumfact"):
+ from dune.codegen.sumfact.switch import sumfact_generate_kernels_per_integral
+ for k in sumfact_generate_kernels_per_integral(integrals):
+ yield k
+ else:
+ yield generate_kernel(integrals)
def extract_kernel_from_cache(tag, name, signature, wrap_in_cgen=True, add_timings=True):
@@ -759,14 +792,6 @@ def local_operator_default_settings(operator, form):
for c in sorted(filter(lambda c: c.count() > 2, form.coefficients()), key=lambda c: c.count()):
name_gridfunction_constructor_argument(c)
- # Set some options!
- from dune.codegen.pdelab.driver import isQuadrilateral
- if isQuadrilateral(form.arguments()[0].ufl_element().cell()) and not get_option("grid_unstructured"):
- from dune.codegen.options import set_form_option
- # For Yasp Grids the jacobian of the transformation is diagonal and constant on each cell
- set_form_option('diagonal_transformation_matrix', True)
- set_form_option('constant_transformation_matrix', True)
-
# Add right base classes for stationary/instationary operators
base_class('Dune::PDELab::LocalOperatorDefaultFlags', classtag="operator")
from dune.codegen.pdelab.driver import is_stationary
@@ -803,7 +828,7 @@ def generate_residual_kernels(form, original_form):
with global_context(integral_type=measure):
from dune.codegen.pdelab.signatures import assembler_routine_name
with global_context(kernel=assembler_routine_name()):
- kernel = [k for k in get_backend(interface="generate_kernels_per_integral")(form.integrals_by_type(measure))]
+ kernel = [k for k in generate_kernels_per_integral(form.integrals_by_type(measure))]
# The integrals might vanish due to unfulfillable boundary conditions.
# We only generate the local operator enums/base classes if they did not.
@@ -890,7 +915,7 @@ def generate_jacobian_kernels(form, original_form):
with global_context(integral_type=measure):
from dune.codegen.pdelab.signatures import assembler_routine_name
with global_context(kernel=assembler_routine_name()):
- kernel = [k for k in get_backend(interface="generate_kernels_per_integral")(jac_apply_form.integrals_by_type(measure))]
+ kernel = [k for k in generate_kernels_per_integral(jac_apply_form.integrals_by_type(measure))]
operator_kernels[(measure, 'jacobian_apply')] = kernel
# Generate dummy functions for those kernels, that vanished in the differentiation process
@@ -908,6 +933,8 @@ def generate_jacobian_kernels(form, original_form):
if get_form_option("sumfact"):
was_sumfact = True
if get_form_option("sumfact_regular_jacobians"):
+ old_geometry_mixins = get_form_option("geometry_mixins")
+ set_form_option("geometry_mixins", "generic")
set_form_option("sumfact", False)
for measure in set(i.integral_type() for i in jacform.integrals()):
if not measure_is_enabled(measure):
@@ -917,7 +944,7 @@ def generate_jacobian_kernels(form, original_form):
with global_context(integral_type=measure):
from dune.codegen.pdelab.signatures import assembler_routine_name
with global_context(kernel=assembler_routine_name()):
- kernel = [k for k in get_backend(interface="generate_kernels_per_integral")(jacform.integrals_by_type(measure))]
+ kernel = [k for k in generate_kernels_per_integral(jacform.integrals_by_type(measure))]
operator_kernels[(measure, 'jacobian')] = kernel
# Generate dummy functions for those kernels, that vanished in the differentiation process
@@ -932,6 +959,7 @@ def generate_jacobian_kernels(form, original_form):
if get_form_option("sumfact_regular_jacobians"):
if was_sumfact:
set_form_option("sumfact", True)
+ set_form_option("geometry_mixins", old_geometry_mixins)
return operator_kernels
diff --git a/python/dune/codegen/pdelab/quadrature.py b/python/dune/codegen/pdelab/quadrature.py
index 8c3e2b36df63b3120883fd51737fc5ded051b157..ae6a7e2db7212b254e3857bc86813ea420ed0251 100644
--- a/python/dune/codegen/pdelab/quadrature.py
+++ b/python/dune/codegen/pdelab/quadrature.py
@@ -1,22 +1,109 @@
import numpy
-from dune.codegen.generation import (backend,
- class_member,
+from dune.codegen.generation import (class_member,
domain,
- get_backend,
get_global_context_value,
globalarg,
iname,
include_file,
instruction,
preamble,
+ quadrature_mixin,
temporary_variable,
valuearg,
)
-from dune.codegen.pdelab.localoperator import lop_template_range_field
+from dune.codegen.pdelab.localoperator import lop_template_range_field, name_ansatz_gfs_constructor_param
from dune.codegen.options import get_form_option
from pymbolic.primitives import Variable, Subscript
+import pymbolic.primitives as prim
+
+
+@quadrature_mixin("base")
+class QuadratureMixinBase(object):
+ def quad_unhandled(self, o):
+ raise NotImplementedError("Quadrature Mixins should handle {}".format(type(o)))
+
+ quadrature_weight = quad_unhandled
+
+ def quadrature_inames(self):
+ raise NotImplementedError("Quadrature mixins should provide the quadrature inames!")
+
+ def quadrature_position(self, index=None):
+ raise NotImplementedError("Quadrature mixins should provide the quadrature position")
+
+
+@quadrature_mixin("generic")
+class GenericQuadratureMixin(QuadratureMixinBase):
+ def quadrature_weight(self, o):
+ # Create a name that has all the necessary quantities mangled in
+ from dune.codegen.pdelab.geometry import local_dimension
+ dim = local_dimension()
+ order = quadrature_order()
+ name = "qw_dim{}_order{}".format(dim, order)
+
+ shape = name_quadrature_bound()
+ globalarg(name, shape=(shape,))
+ self.define_quadrature_weights(name)
+
+ return prim.Subscript(prim.Variable(name),
+ tuple(prim.Variable(i) for i in self.quadrature_inames()))
+
+ @class_member(classtag="operator")
+ def define_quadrature_weights(self, name):
+ rf = lop_template_range_field()
+ from dune.codegen.pdelab.geometry import local_dimension
+ dim = local_dimension()
+ self.eval_quadrature_weights(name)
+ return "mutable std::vector<typename Dune::QuadraturePoint<{}, {}>::Field> {};".format(rf, dim, name)
+
+ @preamble(kernel="operator")
+ def eval_quadrature_weights(self, name):
+ gfs = name_ansatz_gfs_constructor_param()
+ quad_order = quadrature_order()
+ include_file("dune/codegen/localbasiscache.hh", filetag='operatorfile')
+ entity = "{}.gridView().template begin<0>()".format(gfs)
+ if self.measure != "cell":
+ entity = "{}.gridView().ibegin(*({}))".format(gfs, entity)
+ return "fillQuadratureWeightsCache({}->geometry(), {}, {});".format(entity, quad_order, name)
+
+ def quadrature_inames(self):
+ return (quadrature_iname(),)
+
+ def quadrature_position(self, index=None):
+ from dune.codegen.pdelab.geometry import local_dimension
+ dim = local_dimension()
+ order = quadrature_order()
+ name = "qp_dim{}_order{}".format(dim, order)
+
+ shape = (name_quadrature_bound(), dim)
+ globalarg(name, shape=shape, managed=False)
+ self.define_quadrature_points(name)
+
+ qp = prim.Subscript(prim.Variable(name), tuple(prim.Variable(i) for i in self.quadrature_inames()))
+
+ if index is not None:
+ qp = prim.Subscript(qp, (index,))
+
+ return qp
+
+ @class_member(classtag="operator")
+ def define_quadrature_points(self, name):
+ rf = lop_template_range_field()
+ from dune.codegen.pdelab.geometry import local_dimension
+ dim = local_dimension()
+ self.eval_quadrature_points(name)
+ return "mutable std::vector<typename Dune::QuadraturePoint<{}, {}>::Vector> {};".format(rf, dim, name)
+
+ @preamble(kernel="operator")
+ def eval_quadrature_points(self, name):
+ gfs = name_ansatz_gfs_constructor_param()
+ quad_order = quadrature_order()
+ include_file("dune/codegen/localbasiscache.hh", filetag='operatorfile')
+ entity = "{}.gridView().template begin<0>()".format(gfs)
+ if self.measure != "cell":
+ entity = "{}.gridView().ibegin(*({}))".format(gfs, entity)
+ return "fillQuadraturePointsCache({}->geometry(), {}, {});".format(entity, quad_order, name)
@preamble
@@ -56,14 +143,9 @@ def quadrature_iname():
return "q"
-@backend(interface="quad_inames")
-def quadrature_inames():
- return (quadrature_iname(),)
-
-
-def quadrature_preamble(code, **kw):
+def quadrature_preamble(visitor, code, **kw):
kw['tags'] = kw.get('tags', frozenset({})).union(frozenset({"quad"}))
- return instruction(inames=get_backend(interface="quad_inames")(), code=code, **kw)
+ return instruction(inames=visitor.quadrature_inames(), code=code, **kw)
def name_quadrature_point():
@@ -71,114 +153,6 @@ def name_quadrature_point():
return "qp"
-@preamble
-def fill_quadrature_points_cache(name):
- from dune.codegen.pdelab.geometry import name_geometry
- geo = name_geometry()
- quad_order = quadrature_order()
- include_file("dune/codegen/localbasiscache.hh", filetag='operatorfile')
- return "fillQuadraturePointsCache({}, {}, {});".format(geo, quad_order, name)
-
-
-@class_member(classtag="operator")
-def typedef_quadrature_points(name):
- range_field = lop_template_range_field()
- from dune.codegen.pdelab.geometry import local_dimension
- dim = local_dimension()
- return "using {} = typename Dune::QuadraturePoint<{}, {}>::Vector;".format(name, range_field, dim)
-
-
-def type_quadrature_points(name):
- name = name.upper()
- typedef_quadrature_points(name)
- return name
-
-
-@class_member(classtag="operator")
-def define_quadrature_points(name):
- qp_type = type_quadrature_points(name)
- return "mutable std::vector<{}> {};".format(qp_type, name)
-
-
-def name_quadrature_points():
- """Name of vector storing quadrature points as class member"""
- from dune.codegen.pdelab.geometry import local_dimension
- dim = local_dimension()
- order = quadrature_order()
- name = "qp_dim{}_order{}".format(dim, order)
- shape = (name_quadrature_bound(), dim)
- globalarg(name, shape=shape, managed=False)
- define_quadrature_points(name)
- fill_quadrature_points_cache(name)
- return name
-
-
-@backend(interface="quad_pos")
-def pymbolic_quadrature_position():
- quad_points = name_quadrature_points()
- quad_iname = quadrature_iname()
- from pymbolic.primitives import Subscript, Variable
- return Subscript(Variable(quad_points), (Variable(quad_iname),))
-
-
-@backend(interface="qp_in_cell")
-def pymbolic_quadrature_position_in_cell(restriction):
- from dune.codegen.pdelab.geometry import to_cell_coordinates
- quad_pos = get_backend("quad_pos")()
- return to_cell_coordinates(quad_pos, restriction)
-
-
-@preamble
-def fill_quadrature_weights_cache(name):
- from dune.codegen.pdelab.geometry import name_geometry
- geo = name_geometry()
- quad_order = quadrature_order()
- include_file("dune/codegen/localbasiscache.hh", filetag='operatorfile')
- return "fillQuadratureWeightsCache({}, {}, {});".format(geo, quad_order, name)
-
-
-@class_member(classtag="operator")
-def typedef_quadrature_weights(name):
- range_field = lop_template_range_field()
- from dune.codegen.pdelab.geometry import local_dimension
- dim = local_dimension()
- return "using {} = typename Dune::QuadraturePoint<{}, {}>::Field;".format(name, range_field, dim)
-
-
-def pymbolic_quadrature_weight():
- vec = name_quadrature_weights()
- return Subscript(Variable(vec),
- tuple(Variable(i) for i in quadrature_inames()))
-
-
-def type_quadrature_weights(name):
- name = name.upper()
- typedef_quadrature_weights(name)
- return name
-
-
-@class_member(classtag="operator")
-def define_quadrature_weights(name):
- qw_type = type_quadrature_weights(name)
- return "mutable std::vector<{}> {};".format(qw_type, name)
-
-
-def name_quadrature_weights():
- """"Name of vector storing quadrature weights as class member"""
- from dune.codegen.pdelab.geometry import local_dimension
- dim = local_dimension()
- order = quadrature_order()
- name = "qw_dim{}_order{}".format(dim, order)
- define_quadrature_weights(name)
- fill_quadrature_weights_cache(name)
-
- # Quadrature weighs is a globar argument for loopy
- shape = name_quadrature_bound()
- globalarg(name, shape=(shape,))
-
- return name
-
-
def _estimate_quadrature_order():
"""Estimate quadrature order using polynomial degree estimation from UFL"""
# According to UFL documentation estimate_total_polynomial_degree
diff --git a/python/dune/codegen/pdelab/tensors.py b/python/dune/codegen/pdelab/tensors.py
index 754e7b24a8f639a3a73a3532ef56c888368dd2af..a924a39ae1bba19a94c47e5542567d8717cedf57 100644
--- a/python/dune/codegen/pdelab/tensors.py
+++ b/python/dune/codegen/pdelab/tensors.py
@@ -14,61 +14,6 @@ import loopy as lp
import itertools as it
-def define_list_tensor(name, expr, visitor, stack=()):
- for i, child in enumerate(expr.ufl_operands):
- from ufl.classes import ListTensor
- if isinstance(child, ListTensor):
- define_list_tensor(name, child, visitor, stack=stack + (i,))
- else:
- visexpr = visitor.call(child)
- from loopy.symbolic import DependencyMapper
- deps = DependencyMapper(include_subscripts=False, include_lookups=False, include_calls=False)(visexpr)
- instruction(assignee=prim.Subscript(prim.Variable(name), stack + (i,)),
- expression=visitor.call(child),
- forced_iname_deps=frozenset(visitor.interface.quadrature_inames()),
- depends_on=frozenset({lp.match.Tagged("sumfact_stage1")}),
- tags=frozenset({"quad"}),
- )
-
-
-@kernel_cached
-def pymbolic_list_tensor(expr, visitor):
- name = get_counted_variable("listtensor")
- temporary_variable(name,
- shape=expr.ufl_shape,
- managed=True,
- )
- define_list_tensor(name, expr, visitor)
- return prim.Variable(name)
-
-
-@iname
-def identity_iname(name, bound):
- name = "id_{}_{}".format(name, bound)
- domain(name, bound)
- return name
-
-
-def define_identity(name, expr):
- i = identity_iname("i", expr.ufl_shape[0])
- j = identity_iname("j", expr.ufl_shape[1])
- instruction(assignee=prim.Subscript(prim.Variable(name), (prim.Variable(i), prim.Variable(j))),
- expression=prim.If(prim.Comparison(prim.Variable(i), "==", prim.Variable(j)), 1, 0),
- forced_iname_deps_is_final=True,
- )
-
-
-@kernel_cached
-def pymbolic_identity(expr):
- name = "identity_{}_{}".format(expr.ufl_shape[0], expr.ufl_shape[1])
- temporary_variable(name,
- shape=expr.ufl_shape,
- shape_impl=('fm',),
- )
- define_identity(name, expr)
- return prim.Variable(name)
-
-
def define_assembled_tensor(name, expr, visitor):
temporary_variable(name,
shape=expr.ufl_shape,
diff --git a/python/dune/codegen/sumfact/__init__.py b/python/dune/codegen/sumfact/__init__.py
index 22a6049614784ced07c1e2a68739a094ec328516..0b5fe7330a7cd5821646aec853e21fc6df35fe2d 100644
--- a/python/dune/codegen/sumfact/__init__.py
+++ b/python/dune/codegen/sumfact/__init__.py
@@ -1,102 +1,3 @@
-from dune.codegen.generation import get_backend
-from dune.codegen.options import option_switch
-from dune.codegen.pdelab.argument import (name_applycontainer,
- name_coefficientcontainer,
- )
-from dune.codegen.sumfact.quadrature import (quadrature_inames,
- quadrature_weight,
- )
-
-from dune.codegen.sumfact.basis import (lfs_inames,
- pymbolic_basis,
- pymbolic_reference_gradient,
- pymbolic_coefficient,
- pymbolic_coefficient_gradient,
- )
-
+import dune.codegen.sumfact.geometry
import dune.codegen.sumfact.accumulation
import dune.codegen.sumfact.switch
-
-from dune.codegen.pdelab import PDELabInterface
-
-
-class SumFactInterface(PDELabInterface):
- def get_accumulation_info(self, expr, visitor):
- from dune.codegen.sumfact.accumulation import get_accumulation_info
- return get_accumulation_info(expr, visitor)
-
- def list_accumulation_infos(self, expr, visitor):
- from dune.codegen.sumfact.accumulation import list_accumulation_infos
- return list_accumulation_infos(expr, visitor)
-
- def generate_accumulation_instruction(self, expr, visitor):
- from dune.codegen.sumfact.accumulation import generate_accumulation_instruction
- return generate_accumulation_instruction(expr, visitor)
-
- def lfs_inames(self, element, restriction, number=None, context=''):
- return lfs_inames(element, restriction, number, context)
-
- def pymbolic_basis(self, element, restriction, number):
- return pymbolic_basis(element, restriction, number)
-
- def pymbolic_reference_gradient(self, element, restriction, number):
- ret, indices = pymbolic_reference_gradient(element, restriction, number, self.visitor.indices)
- self.visitor.indices = indices
- return ret
-
- def pymbolic_trialfunction_gradient(self, element, restriction, index):
- ret = pymbolic_coefficient_gradient(element, restriction, index, name_coefficientcontainer, self.visitor)
- return ret
-
- def pymbolic_trialfunction(self, element, restriction, index):
- ret = pymbolic_coefficient(element, restriction, index, name_coefficientcontainer, self.visitor)
- return ret
-
- def pymbolic_apply_function_gradient(self, element, restriction, index):
- ret = pymbolic_coefficient_gradient(element, restriction, index, name_applycontainer, self.visitor)
- return ret
-
- def pymbolic_apply_function(self, element, restriction, index):
- ret = pymbolic_coefficient(element, restriction, index, name_applycontainer, self.visitor)
- return ret
-
- def quadrature_inames(self):
- # TODO Is broken at the moment. Remove from interface? See also pdelab interface.
- assert False
- return quadrature_inames()
-
- def pymbolic_quadrature_weight(self):
- return quadrature_weight(self.visitor)
-
- #
- # Geometry related generator functions
- #
-
- def pymbolic_spatial_coordinate(self):
- import dune.codegen.sumfact.geometry
- ret = get_backend(interface="spatial_coordinate", selector=option_switch("diagonal_transformation_matrix"))(self.visitor.do_predicates, self.visitor)
- return ret
-
- def pymbolic_unit_outer_normal(self):
- from dune.codegen.sumfact.geometry import pymbolic_unit_outer_normal
- return pymbolic_unit_outer_normal(self.visitor)
-
- def pymbolic_unit_inner_normal(self):
- from dune.codegen.sumfact.geometry import pymbolic_unit_inner_normal
- return pymbolic_unit_inner_normal(self.visitor)
-
- def pymbolic_facet_jacobian_determinant(self):
- from dune.codegen.sumfact.geometry import pymbolic_facet_jacobian_determinant
- return pymbolic_facet_jacobian_determinant(self.visitor)
-
- def pymbolic_facet_area(self):
- from dune.codegen.sumfact.geometry import pymbolic_facet_area
- return pymbolic_facet_area(self.visitor)
-
- def pymbolic_jacobian_inverse(self, i, j, restriction):
- from dune.codegen.sumfact.geometry import pymbolic_jacobian_inverse
- return pymbolic_jacobian_inverse(i, j, restriction, self.visitor)
-
- def pymbolic_jacobian_determinant(self):
- from dune.codegen.sumfact.geometry import pymbolic_jacobian_determinant
- return pymbolic_jacobian_determinant(self.visitor)
diff --git a/python/dune/codegen/sumfact/accumulation.py b/python/dune/codegen/sumfact/accumulation.py
index 516da6f08b6c2825f616c925448175b61d454a43..eeb383da8fcec4d7799a15befa6f38698376fdb5 100644
--- a/python/dune/codegen/sumfact/accumulation.py
+++ b/python/dune/codegen/sumfact/accumulation.py
@@ -3,7 +3,7 @@ import itertools
from dune.codegen.pdelab.argument import (name_accumulation_variable,
PDELabAccumulationFunction,
)
-from dune.codegen.generation import (backend,
+from dune.codegen.generation import (accumulation_mixin,
domain,
dump_accumulate_timer,
generator_factory,
@@ -26,7 +26,7 @@ from dune.codegen.options import (get_form_option,
from dune.codegen.loopy.flatten import flatten_index
from dune.codegen.loopy.target import type_floatingpoint
from dune.codegen.pdelab.driver import FEM_name_mangling
-from dune.codegen.pdelab.localoperator import determine_accumulation_space
+from dune.codegen.pdelab.localoperator import determine_accumulation_space, AccumulationMixinBase
from dune.codegen.pdelab.restriction import restricted_name
from dune.codegen.pdelab.signatures import assembler_routine_name
from dune.codegen.pdelab.geometry import world_dimension
@@ -46,7 +46,6 @@ from dune.codegen.sumfact.symbolic import SumfactKernel, SumfactKernelInterfaceB
from dune.codegen.ufl.modified_terminals import extract_modified_arguments
from dune.codegen.tools import get_pymbolic_basename, get_leaf
from dune.codegen.error import CodegenError
-from dune.codegen.sumfact.quadrature import quadrature_inames
from pytools import ImmutableRecord, product
@@ -168,8 +167,8 @@ class AccumulationOutput(SumfactKernelInterfaceBase, ImmutableRecord):
if self.trial_element is None:
return ()
else:
- from dune.codegen.sumfact.basis import lfs_inames
- return lfs_inames(get_leaf(self.trial_element, self.trial_element_index), self.restriction)
+ from dune.codegen.sumfact.basis import SumfactBasisMixin
+ return SumfactBasisMixin.lfs_inames(None, get_leaf(self.trial_element, self.trial_element_index), self.restriction)
def realize_input(self, inames, shape, vec_iname, vec_shape, buf, ftags):
# The result of stage 2 has the correct quadrature permutation but no
@@ -351,6 +350,18 @@ def _dof_offset(element, component):
return sum(sizes[0:component])
+@accumulation_mixin("sumfact")
+class SumfactAccumulationMixin(AccumulationMixinBase):
+ def get_accumulation_info(self, expr):
+ return get_accumulation_info(expr, self)
+
+ def list_accumulation_infos(self, expr):
+ return list_accumulation_infos(expr, self)
+
+ def generate_accumulation_instruction(self, expr):
+ return generate_accumulation_instruction(expr, self)
+
+
class SumfactAccumulationInfo(ImmutableRecord):
def __init__(self,
element=None,
@@ -388,10 +399,10 @@ def get_accumulation_info(expr, visitor):
from dune.codegen.pdelab.restriction import Restriction
restriction = Restriction.POSITIVE
- inames = visitor.interface.lfs_inames(leaf_element,
- restriction,
- expr.number()
- )
+ inames = visitor.lfs_inames(leaf_element,
+ restriction,
+ expr.number()
+ )
grad_index = None
if visitor.reference_grad and expr.number() == 0:
@@ -564,7 +575,7 @@ def generate_accumulation_instruction(expr, visitor):
assignee = prim.Subscript(lp.TaggedVariable(temp, vsf.tag), pad)
instruction(assignee=assignee,
expression=0,
- forced_iname_deps=frozenset(quadrature_inames(trial_leaf_element) + jacobian_inames),
+ forced_iname_deps=frozenset(visitor.quadrature_inames() + jacobian_inames),
forced_iname_deps_is_final=True,
tags=frozenset(["quadvec", "gradvec"]),
)
@@ -579,9 +590,10 @@ def generate_accumulation_instruction(expr, visitor):
# with the evaluation of the contribution at all quadrature points
assignee = prim.Subscript(lp.TaggedVariable(temp, vsf.tag),
vsf.quadrature_index(sf, visitor))
+
contrib_dep = instruction(assignee=assignee,
expression=expr,
- forced_iname_deps=frozenset(quadrature_inames(trial_leaf_element) + jacobian_inames),
+ forced_iname_deps=frozenset(visitor.quadrature_inames() + jacobian_inames),
forced_iname_deps_is_final=True,
tags=frozenset({"quadvec", "sumfact_stage2"}).union(vectag),
depends_on=frozenset({deps}).union(frozenset({lp.match.Tagged("sumfact_stage1")})),
diff --git a/python/dune/codegen/sumfact/basis.py b/python/dune/codegen/sumfact/basis.py
index 049e0be9f5adfae7d48ef9adc40a951adcc55924..d8194af77400fc26148263d367fde33f25ed2480 100644
--- a/python/dune/codegen/sumfact/basis.py
+++ b/python/dune/codegen/sumfact/basis.py
@@ -5,7 +5,7 @@ multiplication with the test function is part of the sum factorization kernel.
"""
import itertools
-from dune.codegen.generation import (backend,
+from dune.codegen.generation import (basis_mixin,
domain,
get_counted_variable,
get_counter,
@@ -31,11 +31,11 @@ from dune.codegen.sumfact.permutation import (permute_backward,
sumfact_cost_permutation_strategy,
sumfact_quadrature_permutation_strategy,
)
-from dune.codegen.sumfact.quadrature import quadrature_inames
from dune.codegen.sumfact.switch import (get_facedir,
get_facemod,
)
from dune.codegen.pdelab.argument import name_coefficientcontainer
+from dune.codegen.pdelab.basis import GenericBasisMixin
from dune.codegen.pdelab.geometry import (local_dimension,
world_dimension,
)
@@ -57,6 +57,184 @@ from loopy.match import Writes
import pymbolic.primitives as prim
+@basis_mixin("sumfact")
+class SumfactBasisMixin(GenericBasisMixin):
+ def lfs_inames(self, element, restriction, number=1):
+ from ufl import FiniteElement, TensorProductElement
+ assert isinstance(element, (FiniteElement, TensorProductElement))
+ if number == 0:
+ return ()
+ else:
+ dim = world_dimension()
+ return tuple(sumfact_lfs_iname(element, _basis_functions_per_direction(element)[d], d) for d in range(dim))
+
+ def implement_basis(self, element, restriction, number):
+ # If this is a test function we omit it!
+ if number == 0:
+ return 1
+
+ assert not isinstance(element, MixedElement)
+
+ name = "phi_{}".format(FEM_name_mangling(element))
+ name = restricted_name(name, restriction)
+ self.evaluate_basis(element, name, restriction)
+
+ return prim.Variable(name)
+
+ @kernel_cached
+ def evaluate_basis(self, element, name, restriction):
+ temporary_variable(name, shape=())
+ quad_inames = self.quadrature_inames()
+ inames = self.lfs_inames(element, restriction)
+ facedir = get_facedir(restriction)
+
+ # Collect the pairs of lfs/quad inames that are in use
+ # On facets, the normal direction of the facet is excluded
+ #
+ # If facedir is not none, the length of inames and quad_inames is
+ # different. For inames we want to skip the facedir direction, for
+ # quad_inames we need all entries. Thats the reason for the
+ # help_index.
+ basis_per_dir = _basis_functions_per_direction(element)
+ prod = ()
+ help_index = 0
+ for direction in range(len(inames)):
+ if direction != facedir:
+ prod = prod + (BasisTabulationMatrix(basis_size=basis_per_dir[direction],
+ direction=direction).pymbolic((prim.Variable(quad_inames[help_index]), prim.Variable(inames[direction]))),)
+ help_index += 1
+
+ # Add the missing direction on facedirs by evaluating at either 0 or 1
+ if facedir is not None:
+ facemod = get_facemod(restriction)
+ prod = prod + (prim.Call(PolynomialLookup(name_polynomials(element.degree()), False),
+ (prim.Variable(inames[facedir]), facemod)),)
+
+ # Issue the product
+ instruction(expression=prim.Product(prod),
+ assignee=prim.Variable(name),
+ forced_iname_deps=frozenset(quad_inames + inames),
+ forced_iname_deps_is_final=True,
+ )
+
+ def implement_reference_gradient(self, element, restriction, number):
+ # If this is a test function, we omit it.
+ if number == 0:
+ self.indices = None
+ return 1
+
+ assert len(self.indices) == 1
+ index, = self.indices
+
+ # TODO: Change name?
+ name = "js_{}".format(FEM_name_mangling(element))
+ name = restricted_name(name, restriction)
+ name = "{}_{}".format(name, index)
+ self.evaluate_reference_gradient(element, name, restriction, index)
+
+ self.indices = None
+ return prim.Variable(name)
+
+ @kernel_cached
+ def evaluate_reference_gradient(self, element, name, restriction, index):
+ dim = world_dimension()
+ temporary_variable(name, shape=())
+ quad_inames = self.quadrature_inames()
+ inames = self.lfs_inames(element, restriction)
+ facedir = get_facedir(restriction)
+
+ # Map the direction to a quadrature iname
+ quadinamemapping = {}
+ i = 0
+ for d in range(local_dimension()):
+ if d == facedir:
+ i = i + 1
+ quadinamemapping[i] = quad_inames[d]
+ i = i + 1
+
+ prod = []
+ for i in range(dim):
+ if i != facedir:
+ tab = BasisTabulationMatrix(basis_size=_basis_functions_per_direction(element)[i],
+ direction=i,
+ derivative=index == i)
+ prod.append(tab.pymbolic((prim.Variable(quadinamemapping[i]), prim.Variable(inames[i]))))
+
+ if facedir is not None:
+ facemod = get_facemod(restriction)
+ prod.append(prim.Call(PolynomialLookup(name_polynomials(element.degree()), index == facedir),
+ (prim.Variable(inames[facedir]), facemod)),)
+
+ instruction(assignee=prim.Variable(name),
+ expression=prim.Product(tuple(prod)),
+ forced_iname_deps=frozenset(quad_inames + inames),
+ forced_iname_deps_is_final=True,
+ )
+
+ def implement_trialfunction(self, element, restriction, index):
+ return self.implement_coefficient(element, restriction, index, name_coefficientcontainer)
+
+ def implement_trialfunction_gradient(self, element, restriction, index):
+ derivative = self.indices[0]
+ return self.implement_coefficient(element, restriction, index, name_coefficientcontainer, derivative=derivative)
+
+ def implement_apply_function(self, element, restriction, index):
+ return self.implement_coefficient(element, restriction, index, name_applycontainer)
+
+ def implement_apply_function_gradient(self, element, restriction, index):
+ derivative = self.indices[0]
+ return self.implement_coefficient(element, restriction, index, name_applycontainer, derivative=derivative)
+
+ def implement_coefficient(self, element, restriction, index, coeff_func, derivative=None):
+ sub_element = element
+ if isinstance(element, MixedElement):
+ sub_element = element.extract_component(index)[1]
+ from ufl import FiniteElement
+ assert isinstance(sub_element, (FiniteElement, TensorProductElement))
+
+ # Basis functions per direction
+ basis_size = _basis_functions_per_direction(sub_element)
+
+ # Construct the matrix sequence for this sum factorization
+ matrix_sequence = construct_basis_matrix_sequence(derivative=derivative,
+ facedir=get_facedir(restriction),
+ facemod=get_facemod(restriction),
+ basis_size=basis_size)
+
+ inp = LFSSumfactKernelInput(matrix_sequence=matrix_sequence,
+ coeff_func=coeff_func,
+ element=element,
+ element_index=index,
+ restriction=restriction,
+ )
+
+ position_priority = derivative
+ if position_priority is None:
+ position_priority = 3
+
+ sf = SumfactKernel(matrix_sequence=matrix_sequence,
+ interface=inp,
+ position_priority=position_priority,
+ )
+
+ from dune.codegen.sumfact.vectorization import attach_vectorization_info
+ vsf = attach_vectorization_info(sf)
+
+ self.indices = None
+
+ # If this sum factorization kernel was not used in the dry run we
+ # just return 0
+ if vsf == 0:
+ return 0
+
+ # Add a sum factorization kernel that implements the evaluation of
+ # the basis functions at quadrature points (stage 1)
+ from dune.codegen.sumfact.realization import realize_sum_factorization_kernel
+ var, _ = realize_sum_factorization_kernel(vsf)
+
+ return prim.Subscript(var, vsf.quadrature_index(sf, self))
+
+
class LFSSumfactKernelInput(SumfactKernelInterfaceBase, ImmutableRecord):
def __init__(self,
matrix_sequence=None,
@@ -260,99 +438,6 @@ def _basis_functions_per_direction(element):
return basis_size
-def pymbolic_coefficient_gradient(element, restriction, index, coeff_func, visitor):
- sub_element = element
- grad_index = visitor.indices[0]
- if isinstance(element, MixedElement):
- sub_element = element.extract_component(index)[1]
-
- from ufl import FiniteElement
- assert isinstance(sub_element, (FiniteElement, TensorProductElement))
-
- # Number of basis functions per direction
- basis_size = _basis_functions_per_direction(sub_element)
-
- # Construct the matrix sequence for this sum factorization
- matrix_sequence = construct_basis_matrix_sequence(derivative=grad_index,
- facedir=get_facedir(restriction),
- facemod=get_facemod(restriction),
- basis_size=basis_size,
- )
-
- inp = LFSSumfactKernelInput(matrix_sequence=matrix_sequence,
- coeff_func=coeff_func,
- element=element,
- element_index=index,
- restriction=restriction,
- )
-
- # The sum factorization kernel object gathering all relevant information
- sf = SumfactKernel(matrix_sequence=matrix_sequence,
- position_priority=grad_index,
- interface=inp,
- )
-
- from dune.codegen.sumfact.vectorization import attach_vectorization_info
- vsf = attach_vectorization_info(sf)
-
- # If this sum factorization kernel was not used in the dry run we
- # just return 0
- if vsf == 0:
- visitor.indices = None
- return 0
-
- from dune.codegen.sumfact.realization import realize_sum_factorization_kernel
- var, insn_dep = realize_sum_factorization_kernel(vsf)
-
- visitor.indices = None
- return prim.Subscript(var, vsf.quadrature_index(sf, visitor))
-
-
-def pymbolic_coefficient(element, restriction, index, coeff_func, visitor):
- sub_element = element
- if isinstance(element, MixedElement):
- sub_element = element.extract_component(index)[1]
- from ufl import FiniteElement
- assert isinstance(sub_element, (FiniteElement, TensorProductElement))
-
- # Basis functions per direction
- basis_size = _basis_functions_per_direction(sub_element)
-
- # Construct the matrix sequence for this sum factorization
- matrix_sequence = construct_basis_matrix_sequence(facedir=get_facedir(restriction),
- facemod=get_facemod(restriction),
- basis_size=basis_size)
-
- inp = LFSSumfactKernelInput(matrix_sequence=matrix_sequence,
- coeff_func=coeff_func,
- element=element,
- element_index=index,
- restriction=restriction,
- )
-
- sf = SumfactKernel(matrix_sequence=matrix_sequence,
- interface=inp,
- position_priority=3,
- )
-
- from dune.codegen.sumfact.vectorization import attach_vectorization_info
- vsf = attach_vectorization_info(sf)
-
- # If this sum factorization kernel was not used in the dry run we
- # just return 0
- if vsf == 0:
- visitor.indices = None
- return 0
-
- # Add a sum factorization kernel that implements the evaluation of
- # the basis functions at quadrature points (stage 1)
- from dune.codegen.sumfact.realization import realize_sum_factorization_kernel
- var, _ = realize_sum_factorization_kernel(vsf)
-
- visitor.indices = None
- return prim.Subscript(var, vsf.quadrature_index(sf, visitor))
-
-
@iname
def sumfact_lfs_iname(element, bound, dim):
assert(isinstance(bound, int))
@@ -360,121 +445,3 @@ def sumfact_lfs_iname(element, bound, dim):
name = "sumfac_lfs_{}_{}".format(FEM_name_mangling(element), dim)
domain(name, bound)
return name
-
-
-@backend(interface="lfs_inames", name="sumfact")
-def lfs_inames(element, restriction, number=1, context=''):
- from ufl import FiniteElement, TensorProductElement
- assert isinstance(element, (FiniteElement, TensorProductElement))
- if number == 0:
- return ()
- else:
- dim = world_dimension()
- return tuple(sumfact_lfs_iname(element, _basis_functions_per_direction(element)[d], d) for d in range(dim))
-
-
-@backend(interface="evaluate_basis", name="sumfact")
-@kernel_cached
-def evaluate_basis(element, name, restriction):
- temporary_variable(name, shape=())
- quad_inames = quadrature_inames(element)
- inames = lfs_inames(element, restriction)
- facedir = get_facedir(restriction)
-
- # Collect the pairs of lfs/quad inames that are in use
- # On facets, the normal direction of the facet is excluded
- #
- # If facedir is not none, the length of inames and quad_inames is
- # different. For inames we want to skip the facedir direction, for
- # quad_inames we need all entries. Thats the reason for the
- # help_index.
- basis_per_dir = _basis_functions_per_direction(element)
- prod = ()
- help_index = 0
- for direction in range(len(inames)):
- if direction != facedir:
- prod = prod + (BasisTabulationMatrix(basis_size=basis_per_dir[direction],
- direction=direction).pymbolic((prim.Variable(quad_inames[help_index]), prim.Variable(inames[direction]))),)
- help_index += 1
-
- # Add the missing direction on facedirs by evaluating at either 0 or 1
- if facedir is not None:
- facemod = get_facemod(restriction)
- prod = prod + (prim.Call(PolynomialLookup(name_polynomials(element.degree()), False),
- (prim.Variable(inames[facedir]), facemod)),)
-
- # Issue the product
- instruction(expression=prim.Product(prod),
- assignee=prim.Variable(name),
- forced_iname_deps=frozenset(quad_inames + inames),
- forced_iname_deps_is_final=True,
- )
-
-
-def pymbolic_basis(element, restriction, number):
- # If this is a test function we omit it!
- if number == 0:
- return 1
-
- assert not isinstance(element, MixedElement)
-
- name = "phi_{}".format(FEM_name_mangling(element))
- name = restricted_name(name, restriction)
- evaluate_basis(element, name, restriction)
-
- return prim.Variable(name)
-
-
-@backend(interface="evaluate_grad", name="sumfact")
-@kernel_cached
-def evaluate_reference_gradient(element, name, restriction, index):
- dim = world_dimension()
- temporary_variable(name, shape=())
- quad_inames = quadrature_inames(element)
- inames = lfs_inames(element, restriction)
- facedir = get_facedir(restriction)
-
- # Map the direction to a quadrature iname
- quadinamemapping = {}
- i = 0
- for d in range(local_dimension()):
- if d == facedir:
- i = i + 1
- quadinamemapping[i] = quad_inames[d]
- i = i + 1
-
- prod = []
- for i in range(dim):
- if i != facedir:
- tab = BasisTabulationMatrix(basis_size=_basis_functions_per_direction(element)[i],
- direction=i,
- derivative=index == i)
- prod.append(tab.pymbolic((prim.Variable(quadinamemapping[i]), prim.Variable(inames[i]))))
-
- if facedir is not None:
- facemod = get_facemod(restriction)
- prod.append(prim.Call(PolynomialLookup(name_polynomials(element.degree()), index == facedir),
- (prim.Variable(inames[facedir]), facemod)),)
-
- instruction(assignee=prim.Variable(name),
- expression=prim.Product(tuple(prod)),
- forced_iname_deps=frozenset(quad_inames + inames),
- forced_iname_deps_is_final=True,
- )
-
-
-def pymbolic_reference_gradient(element, restriction, number, indices):
- # If this is a test function, we omit it.
- if number == 0:
- return 1, None
-
- assert len(indices) == 1
- index, = indices
-
- # TODO: Change name?
- name = "js_{}".format(FEM_name_mangling(element))
- name = restricted_name(name, restriction)
- name = "{}_{}".format(name, index)
- evaluate_reference_gradient(element, name, restriction, index)
-
- return prim.Variable(name), None
diff --git a/python/dune/codegen/sumfact/geometry.py b/python/dune/codegen/sumfact/geometry.py
index b95b348c43d5820aea9e7e6b2d32491f1fc03023..cef07354adb7dbfbbb85d42e72fc5e7e14d57a08 100644
--- a/python/dune/codegen/sumfact/geometry.py
+++ b/python/dune/codegen/sumfact/geometry.py
@@ -1,8 +1,8 @@
""" Sum factorized geometry evaluations """
-from dune.codegen.generation import (backend,
- class_member,
+from dune.codegen.generation import (class_member,
domain,
+ geometry_mixin,
get_counted_variable,
iname,
include_file,
@@ -17,10 +17,14 @@ from dune.codegen.generation import (backend,
)
from dune.codegen.loopy.flatten import flatten_index
from dune.codegen.options import get_option
-from dune.codegen.pdelab.geometry import (local_dimension,
+from dune.codegen.pdelab.geometry import (enforce_boundary_restriction,
+ local_dimension,
world_dimension,
name_cell_geometry,
name_geometry,
+ GenericPDELabGeometryMixin,
+ AxiparallelGeometryMixin,
+ EquidistantGeometryMixin,
)
from dune.codegen.pdelab.localoperator import (name_ansatz_gfs_constructor_param,
lop_template_ansatz_gfs,
@@ -34,14 +38,9 @@ from dune.codegen.sumfact.permutation import (permute_backward,
sumfact_cost_permutation_strategy,
sumfact_quadrature_permutation_strategy,
)
-from dune.codegen.sumfact.quadrature import (additional_inames,
- default_quadrature_inames)
-from dune.codegen.sumfact.realization import (name_buffer_storage,
- realize_sum_factorization_kernel,
- )
+from dune.codegen.sumfact.quadrature import additional_inames
from dune.codegen.sumfact.switch import get_facedir, get_facemod
from dune.codegen.sumfact.symbolic import SumfactKernelInterfaceBase, SumfactKernel
-from dune.codegen.sumfact.vectorization import attach_vectorization_info
from dune.codegen.tools import get_pymbolic_basename
from dune.codegen.options import get_form_option, option_switch
from dune.codegen.ufl.modified_terminals import Restriction
@@ -52,6 +51,274 @@ from loopy.match import Writes
import pymbolic.primitives as prim
import numpy as np
+import loopy as lp
+
+
+@geometry_mixin("sumfact_multilinear")
+class SumfactMultiLinearGeometryMixin(GenericPDELabGeometryMixin):
+ def _jacobian_inverse(self):
+ return "jit"
+
+ def jacobian_inverse(self, o):
+ # This differs from the generic one by not falling back to the
+ # transpose of the jacobian inverse.
+ restriction = enforce_boundary_restriction(self)
+
+ assert(len(self.indices) == 2)
+ i, j = self.indices
+ self.indices = None
+
+ name = restricted_name(self._jacobian_inverse(), restriction)
+ self.define_jacobian_inverse(name, restriction)
+
+ return prim.Subscript(prim.Variable(name), (i, j))
+
+ def _jacobian_determinant(self):
+ return "detjac"
+
+ def define_jacobian_determinant(self, o):
+ restriction = Restriction.NONE if self.measure == "cell" else Restriction.POSITIVE
+ self.define_jacobian_inverse(self._jacobian_inverse(), restriction)
+ return prim.Variable(self._jacobian_determinant())
+
+ def define_jacobian_inverse(self, name, restriction):
+ """Return jacobian inverse of the geometry mapping (of a cell)
+
+ At the moment this only works for geometry mappings of cells and not for
+ intersection. We only consider this case as it greatly simplifies code
+ generation for unstructured grids by making the grid edge consistent and
+ avoiding the face-twist problem.
+ """
+ # Calculate the jacobian of the geometry mapping using sum factorization
+ dim = world_dimension()
+ names_jacobian = []
+ for j in range(dim):
+ for i in range(dim):
+ name_jacobian = restricted_name("jacobian_geometry_{}{}".format(i, j), restriction)
+ temporary_variable(name_jacobian)
+ assignee = prim.Variable(name_jacobian)
+ expression = _name_jacobian(i, j, restriction, self)
+ inames = self.quadrature_inames() + additional_inames(self)
+ instruction(assignee=assignee,
+ expression=expression,
+ within_inames=frozenset(inames),
+ depends_on=frozenset({lp.match.Tagged("sumfact_stage1")})
+ )
+ names_jacobian.append(name_jacobian)
+
+ # The sum factorization kernels from the geometry evaluation of the
+ # jacobians will never appear in the expression for the input of
+ # stage 3. This way the SumfactCollectMapper will never see them
+ # and they will be marked as inactive. Here we explicitly mark the
+ # as used.
+ basis_sf_kernels(expression.aggregate)
+
+ # Calculate the inverse of the jacobian of the geometry mapping and the
+ # determinant by calling a c++ function. Note: The result will be column
+ # major -> fortran style.
+ name_detjac = self._jacobian_determinant()
+ temporary_variable(name_detjac, shape=())
+ ftags = "f,f"
+ temporary_variable(name, shape=(dim, dim), dim_tags=ftags, managed=True)
+ include_file('dune/codegen/sumfact/invertgeometry.hh', filetag='operatorfile')
+ code = "{} = invert_and_return_determinant({}, {});".format(name_detjac,
+ ", ".join(names_jacobian),
+ name)
+ silenced_warning("read_no_write({})".format(name_detjac))
+ inames = self.quadrature_inames() + additional_inames(self)
+ return instruction(code=code,
+ assignees=frozenset([name, name_detjac]),
+ read_variables=frozenset(names_jacobian),
+ inames=frozenset(inames),
+ )
+
+ def facet_jacobian_determinant(self, o):
+ """Absolute value of determinant of jacobian of facet geometry mapping
+
+ Calculate the absolute value of the determinant of the jacobian of the
+ geometry mapping from the reference cell of the intersection to the
+ intersection:
+
+ |\det(\nabla \mu_F)|
+
+ This is done using the surface metric tensor lemma from the lecture notes
+ of Jean-Luc Guermond:
+
+ |\det(\nabla \mu_F)| = \| \tilde{n} \| |\det(\nabla \mu_{T_s})|
+
+ Here \tilde{n} is the outer normal defined by:
+
+ \tilde{n} = (\nabla \mu_{T_s})^{-T} \hat{n}_s
+
+ where \hat{n}_s is the unit outer normal of the corresponding face of the
+ reference cell.
+ """
+ detjac = self.jacobian_determinant(None)
+ norm_of_outer_normal = normalize(self.outer_normal(), world_dimension())
+
+ return prim.Product((detjac, norm_of_outer_normal))
+
+ def outer_normal(self):
+ """ This is the *unnormalized* outer normal """
+ name = "outer_normal"
+ facedir_s = get_facedir(Restriction.POSITIVE)
+ facemod_s = get_facemod(Restriction.POSITIVE)
+
+ temporary_variable(name, shape=(world_dimension(),))
+ for i in range(world_dimension()):
+ assignee = prim.Subscript(prim.Variable(name), i)
+ self.indices = (facedir_s, i)
+ ji = self.jacobian_inverse(None)
+
+ # Note: 2*facemod_s-1 because of
+ # -1 if facemod_s = 0
+ # +1 if facemod_s = 1
+ expression = ji * (2 * facemod_s - 1)
+
+ inames = self.quadrature_inames() + additional_inames(self)
+ instruction(assignee=assignee,
+ expression=expression,
+ within_inames=frozenset(inames),
+ )
+
+ return prim.Variable(name)
+
+ def facet_normal(self, o):
+ index, = self.indices
+ self.indices = None
+
+ outer_normal = self.outer_normal()
+ norm = normalize(outer_normal, world_dimension())
+
+ return prim.Subscript(outer_normal, (index,)) / norm
+
+ def spatial_coordinate(self, o):
+ """Calculate global coordinates of quadrature points for multilinear geometry mappings
+
+ The evalualation is done using sum factorization techniques.
+
+ On facets we use the geometry mapping of the self cell to get a global
+ evaluation of the quadrature points. Avoiding the geometry mapping of the
+ face itself greatly simplifies the generation of code for unstructured
+ grids. Instead of looking at all the possible embeddings of the reference
+ element of the face into the reference element of the cell we can make the
+ grid edge consistent and choose a suitable convention for the order of
+ directions in our sum factorization kernels.
+ """
+ assert len(self.indices) == 1
+ restriction = enforce_boundary_restriction(self)
+
+ # Generate sum factorization kernel and add vectorization info
+ matrix_sequence = construct_basis_matrix_sequence(facedir=get_facedir(restriction),
+ facemod=get_facemod(restriction),
+ basis_size=(2,) * world_dimension())
+ inp = GeoCornersInput(matrix_sequence=matrix_sequence,
+ direction=self.indices[0],
+ restriction=restriction,
+ )
+ sf = SumfactKernel(matrix_sequence=matrix_sequence,
+ interface=inp,
+ )
+
+ from dune.codegen.sumfact.vectorization import attach_vectorization_info
+ vsf = attach_vectorization_info(sf)
+
+ # If this sum factorization kernel was not used in the dry run we
+ # just return 0
+ if vsf == 0:
+ self.indices = None
+ return 0
+
+ # Add a sum factorization kernel that implements the evaluation of
+ # the basis functions at quadrature points (stage 1)
+ from dune.codegen.sumfact.realization import realize_sum_factorization_kernel
+ var, _ = realize_sum_factorization_kernel(vsf)
+
+ # The results of this function is already the right component of the
+ # spatial coordinate and we don't need to index this in the visitor
+ self.indices = None
+
+ return prim.Subscript(var, vsf.quadrature_index(sf, self))
+
+
+@geometry_mixin("sumfact_axiparallel")
+class SumfactAxiParallelGeometryMixin(AxiparallelGeometryMixin):
+ def facet_normal(self, o):
+ i, = self.indices
+ self.indices = None
+ assert isinstance(i, int)
+
+ # Use facemod_s and facedir_s
+ if i == get_facedir(Restriction.POSITIVE):
+ if get_facemod(Restriction.POSITIVE):
+ return 1
+ else:
+ return -1
+ else:
+ return 0
+
+
+@geometry_mixin("sumfact_equidistant")
+class SumfactEqudistantGeometryMixin(EquidistantGeometryMixin, SumfactAxiParallelGeometryMixin):
+ def facet_jacobian_determinant(self, o):
+ name = "fdetjac"
+ self.define_facet_jacobian_determinant(name)
+ facedir = get_facedir(Restriction.POSITIVE)
+ globalarg(name, shape=(world_dimension(),))
+ return prim.Subscript(prim.Variable(name), (facedir,))
+
+ @class_member(classtag="operator")
+ def define_facet_jacobian_determinant(self, name):
+ self._define_facet_jacobian_determinant_eval(name)
+ gfst = lop_template_ansatz_gfs()
+ return "std::array<typename {}::Traits::GridView::template Codim<0>::Geometry::ctype, {}> {};".format(gfst, world_dimension(), name)
+
+ @kernel_cached(kernel="operator")
+ def _define_facet_jacobian_determinant_eval(self, name):
+ gfs = name_ansatz_gfs_constructor_param()
+ rft = lop_template_range_field()
+ code = ["{",
+ " auto e = *({}.gridView().template begin<0>());".format(gfs),
+ " int dir=0;",
+ " for(auto is = {0}.gridView().ibegin(e); is != {0}.gridView().iend(e); ++(++is), ++dir)".format(gfs),
+ " {}[dir] = is->geometry().integrationElement(Dune::FieldVector<{}, {}>());".format(name, rft, world_dimension() - 1),
+ "}",
+ ]
+ instruction(code="\n".join(code),
+ kernel="operator",
+ )
+
+ def spatial_coordinate(self, o):
+ index, = self.indices
+
+ # Urgh: *SOMEHOW* construct a face direction. This is not breaking in the unstructured
+ # case, because we do not enter this code path...
+ from dune.codegen.pdelab.restriction import Restriction
+ restriction = Restriction.NONE
+ if self.measure == "interior_facet":
+ restriction = Restriction.POSITIVE
+ from dune.codegen.sumfact.switch import get_facedir
+ face = get_facedir(restriction)
+
+ lowcorner = name_lowerleft_corner()
+ meshwidth = name_meshwidth()
+
+ # If we have to decide which boundary condition to take for this
+ # intersection we always take the boundary condition of the center
+ # of the intersection. We assume that there are no cells with more
+ # than one boundary condition.
+ if self.do_predicates:
+ x = 0.5
+ elif index == face:
+ x = 0
+ else:
+ iindex = index
+ if face is not None and index > face:
+ iindex = iindex - 1
+ x = self.quadrature_position(iindex)
+
+ self.indices = None
+ return prim.Subscript(prim.Variable(lowcorner), (index,)) + x * prim.Subscript(prim.Variable(meshwidth), (index,))
@iname
@@ -154,6 +421,7 @@ class GeoCornersInput(SumfactKernelInterfaceBase, ImmutableRecord):
# The array that will be passed to the sum factorization kernel
# function should contain the coefficients in the cost permuted order!
+ from dune.codegen.sumfact.realization import name_buffer_storage
name = "input_{}".format(sf.buffer)
ftags = ",".join(["f"] * (sf.length + 1))
temporary_variable(name,
@@ -205,56 +473,6 @@ class GeoCornersInput(SumfactKernelInterfaceBase, ImmutableRecord):
return prim.Subscript(prim.Variable(inp), inames + vec_iname)
-@backend(interface="spatial_coordinate", name="default")
-def pymbolic_spatial_coordinate_multilinear(do_predicates, visitor):
- """Calculate global coordinates of quadrature points for multilinear geometry mappings
-
- The evalualation is done using sum factorization techniques.
-
- On facets we use the geometry mapping of the self cell to get a global
- evaluation of the quadrature points. Avoiding the geometry mapping of the
- face itself greatly simplifies the generation of code for unstructured
- grids. Instead of looking at all the possible embeddings of the reference
- element of the face into the reference element of the cell we can make the
- grid edge consistent and choose a suitable convention for the order of
- directions in our sum factorization kernels.
- """
- assert len(visitor.indices) == 1
-
- # Adjust restriction for boundary facets
- restriction = visitor.restriction
- if get_global_context_value("integral_type") == "exterior_facet":
- restriction = 1
-
- # Generate sum factorization kernel and add vectorization info
- matrix_sequence = construct_basis_matrix_sequence(facedir=get_facedir(restriction),
- facemod=get_facemod(restriction),
- basis_size=(2,) * world_dimension())
- inp = GeoCornersInput(matrix_sequence=matrix_sequence,
- direction=visitor.indices[0],
- restriction=restriction)
- sf = SumfactKernel(matrix_sequence=matrix_sequence,
- interface=inp,
- )
- vsf = attach_vectorization_info(sf)
-
- # If this sum factorization kernel was not used in the dry run we
- # just return 0
- if vsf == 0:
- visitor.indices = None
- return 0
-
- # Add a sum factorization kernel that implements the evaluation of
- # the basis functions at quadrature points (stage 1)
- var, _ = realize_sum_factorization_kernel(vsf)
-
- # The results of this function is already the right component of the
- # spatial coordinate and we don't need to index this in the visitor
- visitor.indices = None
-
- return prim.Subscript(var, vsf.quadrature_index(sf, visitor))
-
-
@preamble
def define_corner(name, low):
geo = name_geometry()
@@ -298,299 +516,6 @@ def name_meshwidth():
return name
-@backend(interface="spatial_coordinate", name="diagonal_transformation_matrix")
-def pymbolic_spatial_coordinate_axiparallel(do_predicates, visitor):
- assert len(visitor.indices) == 1
- index, = visitor.indices
-
- # Urgh: *SOMEHOW* construct a face direction
- from dune.codegen.pdelab.restriction import Restriction
- restriction = Restriction.NONE
- if get_global_context_value("integral_type") == "interior_facet":
- restriction = Restriction.POSITIVE
- from dune.codegen.sumfact.switch import get_facedir
- face = get_facedir(restriction)
-
- lowcorner = name_lowerleft_corner()
- meshwidth = name_meshwidth()
-
- # If we have to decide which boundary condition to take for this
- # intersection we always take the boundary condition of the center
- # of the intersection. We assume that there are no cells with more
- # than one boundary condition.
- if do_predicates:
- x = 0.5
- elif index == face:
- x = 0
- else:
- iindex = index
- if face is not None and index > face:
- iindex = iindex - 1
- from dune.codegen.sumfact.quadrature import pymbolic_indexed_quadrature_position
- x = pymbolic_indexed_quadrature_position(iindex, visitor)
-
- visitor.indices = None
- return prim.Subscript(prim.Variable(lowcorner), (index,)) + x * prim.Subscript(prim.Variable(meshwidth), (index,))
-
-
-def define_outer_normal(name, visitor):
- """Calculate outer normal (length not necessarily 1)
-
- This will be used to calculate the unit outer normal and the determinant of
- the Jacobian of the geometry mapping of the face.
-
- The outer normal is calculated by multiplying the Jacobian inverse
- transposed of the geometry mapping of the 'self'-cell with the unit outer
- normal \hat{n_s} of the face corresponding to the intersetcion in the
- reference element of the 'self'-cell:
-
- n = (\nabla \mu_{T_s})^{-T} \hat{n_s}
-
- Instead of setting up \hat{n_s} we just look at facedir_s and facemod_s.
- """
- facedir_s = get_facedir(Restriction.POSITIVE)
- facemod_s = get_facemod(Restriction.POSITIVE)
- temporary_variable(name, shape=(world_dimension(),))
- for i in range(world_dimension()):
- assignee = prim.Subscript(prim.Variable(name), i)
- ji = pymbolic_jacobian_inverse(facedir_s, i, Restriction.POSITIVE, visitor)
-
- # Note: 2*facemod_s-1 because of
- # -1 if facemod_s = 0
- # +1 if facemod_s = 1
- expression = ji * (2 * facemod_s - 1)
-
- inames = default_quadrature_inames(visitor) + additional_inames(visitor)
- instruction(assignee=assignee,
- expression=expression,
- within_inames=frozenset(inames),
- )
-
-
-def name_outer_normal(visitor):
- name = "outer_normal"
- define_outer_normal(name, visitor)
- return name
-
-
-def pymbolic_outer_normal(visitor):
- name = name_outer_normal(visitor)
- return prim.Variable(name)
-
-
-def define_norm_of_outer_normal(name, visitor):
- outer_normal = pymbolic_outer_normal(visitor)
-
- # Calculate the norm of outer_normal
- temporary_variable(name, shape=())
- expression = 0
- for i in range(world_dimension()):
- expression = prim.Sum((expression, prim.Subscript(outer_normal, i) * prim.Subscript(outer_normal, i)))
- expression = prim.Call(prim.Variable("sqrt"), (expression,))
- inames = default_quadrature_inames(visitor) + additional_inames(visitor)
-
- # Apparently the call of sqrt shadows the dependency on outer normal. Add manually.
- instruction(assignee=prim.Variable(name),
- expression=expression,
- depends_on=frozenset({Writes(get_pymbolic_basename(outer_normal))}),
- within_inames=frozenset(inames),
- )
-
-
-def name_norm_of_outer_normal(visitor):
- name = "norm_of_outer_normal"
- define_norm_of_outer_normal(name, visitor)
- return name
-
-
-def pymbolic_norm_of_outer_normal(visitor):
- name = name_norm_of_outer_normal(visitor)
- return prim.Variable(name)
-
-
-def pymbolic_unit_outer_normal(visitor):
- assert (len(visitor.indices) is 1)
- index = visitor.indices[0]
- assert isinstance(index, int)
- if get_form_option("diagonal_transformation_matrix"):
- # In this case we just return -1, 0 or 1 depending on the current
- # facemod and facedir. We don't need to (and cannot) index these ints
- # so we set the indices of the visitor to None.
- visitor.indices = None
-
- # Use facemod_s and facedir_s
- from dune.codegen.sumfact.switch import get_facedir, get_facemod
- if index == get_facedir(Restriction.POSITIVE):
- if get_facemod(Restriction.POSITIVE):
- return 1
- else:
- return -1
- else:
- return 0
- else:
- # In pymbolic_outer_normal we need to create the
- # jacobian_inverse_transposed which doesn't play nicely with the
- # visitor having indices. Maybe this should be improved. For now let's
- # just set them to None and restore them afterwards.
- visitor.indices = None
- outer_normal = pymbolic_outer_normal(visitor)
- norm = pymbolic_norm_of_outer_normal(visitor)
- visitor.indices = (index,)
-
- # Create unit_outer_normal by scaling outer_normal to length 1
- name = "unit_outer_normal"
- temporary_variable(name, shape=(world_dimension(),))
- inames = default_quadrature_inames(visitor) + additional_inames(visitor)
- for i in range(world_dimension()):
- instruction(assignee=prim.Subscript(prim.Variable(name), i),
- expression=prim.Subscript(outer_normal, i) / norm,
- within_inames=frozenset(inames),
- )
-
- return prim.Variable(name)
-
-
-def pymbolic_unit_inner_normal(visitor_indices):
- # Note: The implementation of this was adjusted for unstructured grids but
- # it was not properly tested. If you need the unit inner normal just remove
- # the assert(False) statement but make sure to verify that this does the
- # right thing.
- assert(False)
-
- assert (len(visitor.indices) is 1)
- index = visitor.indices[0]
- assert isinstance(index, int)
- if get_form_option("diagonal_transformation_matrix"):
- # In this case we just return -1, 0 or 1 depending on the current
- # facemod and facedir. We don't need to (and cannot) index these ints
- # so we set the indices of the visitor to None.
- visitor.indices = None
- if index == get_facedir(Restriction.POSITIVE):
- if get_facemod(Restriction.POSITIVE):
- return -1
- else:
- return 1
- else:
- return 0
- else:
- # In pymbolic_outer_normal we need to create the
- # jacobian_inverse_transposed which doesn't play nicely with the
- # visitor having indices. Maybe this should be improved. For now let's
- # just set them to None and restore them afterwards.
- visitor.indices = None
- outer_normal = pymbolic_outer_normal(visitor)
- norm = pymbolic_norm_of_outer_normal(visitor)
- visitor.indices = (index,)
-
- # Create unit_outer_normal by scaling outer_normal to length 1
- name = "unit_outer_normal"
- temporary_variable(name, shape=(world_dimension(),))
-
- inames = default_quadrature_inames(visitor) + additional_inames(visitor)
- for i in range(world_dimension()):
- instruction(assignee=prim.Subscript(prim.Variable(name), i),
- expression=-1 * prim.Subscript(outer_normal, i) / norm,
- within_inames=frozenset(inames),
- )
-
- return prim.Variable(name)
-
-
-@kernel_cached(kernel="operator")
-def define_constant_facet_jacobian_determinant_eval(name):
- gfs = name_ansatz_gfs_constructor_param()
- rft = lop_template_range_field()
- code = ["{",
- " auto e = *({}.gridView().template begin<0>());".format(gfs),
- " int dir=0;",
- " for(auto is = {0}.gridView().ibegin(e); is != {0}.gridView().iend(e); ++(++is), ++dir)".format(gfs),
- " {}[dir] = is->geometry().integrationElement(Dune::FieldVector<{}, {}>());".format(name, rft, world_dimension() - 1),
- "}",
- ]
- instruction(code="\n".join(code),
- kernel="operator",
- )
-
-
-@class_member(classtag="operator")
-def define_constant_facet_jacobian_determinant(name):
- define_constant_facet_jacobian_determinant_eval(name)
- gfst = lop_template_ansatz_gfs()
- return "std::array<typename {}::Traits::GridView::template Codim<0>::Geometry::ctype, {}> {};".format(gfst, world_dimension(), name)
-
-
-def pymbolic_constant_facet_jacobian_determinant():
- facedir = get_facedir(Restriction.POSITIVE)
- assert isinstance(facedir, int)
-
- name = "fdetjac"
- define_constant_facet_jacobian_determinant(name)
- globalarg(name, shape=(world_dimension(),))
- return prim.Subscript(prim.Variable(name), (facedir,))
-
-
-def define_facet_jacobian_determinant(name, visitor):
- """Absolute value of determinant of jacobian of facet geometry mapping
-
- Calculate the absolute value of the determinant of the jacobian of the
- geometry mapping from the reference cell of the intersection to the
- intersection:
-
- |\det(\nabla \mu_F)|
-
- This is done using the surface metric tensor lemma from the lecture notes
- of Jean-Luc Guermond:
-
- |\det(\nabla \mu_F)| = \| \tilde{n} \| |\det(\nabla \mu_{T_s})|
-
- Here \tilde{n} is the outer normal defined by:
-
- \tilde{n} = (\nabla \mu_{T_s})^{-T} \hat{n}_s
-
- where \hat{n}_s is the unit outer normal of the corresponding face of the
- reference cell.
- """
- # Calculate norm of outer normal and determinant of jacobian of geometry
- # mapping of the inside cell -> set restriction to 1 for boundary facets
- if get_global_context_value("integral_type") == "exterior_facet":
- assert visitor.restriction is 0
- visitor.restriction = 1
- norm_of_outer_normal = pymbolic_norm_of_outer_normal(visitor)
- detjac = pymbolic_jacobian_determinant(visitor)
- if get_global_context_value("integral_type") == "exterior_facet":
- visitor.restriction = 0
-
- # Multiply detjac and norm_of_outer_normal
- temporary_variable(name, shape=())
- inames = default_quadrature_inames(visitor) + additional_inames(visitor)
- instruction(assignee=prim.Variable(name),
- expression=detjac * norm_of_outer_normal,
- within_inames=frozenset(inames),
- )
-
-
-def name_facet_jacobian_determinant(visitor):
- name = "fdetjac"
- define_facet_jacobian_determinant(name, visitor)
- return name
-
-
-def pymbolic_facet_jacobian_determinant(visitor):
- if get_form_option("constant_transformation_matrix"):
- return pymbolic_constant_facet_jacobian_determinant()
- else:
- name = name_facet_jacobian_determinant(visitor)
- return prim.Variable(name)
-
-
-def pymbolic_facet_area(visitor):
- if get_form_option("constant_transformation_matrix"):
- return pymbolic_facet_jacobian_determinant(visitor)
- else:
- from dune.codegen.pdelab.geometry import pymbolic_facet_area as nonconstant_facet_area
- return nonconstant_facet_area()
-
-
def _name_jacobian(i, j, restriction, visitor):
"""Return the (i, j) component of the jacobian of the geometry mapping
@@ -615,6 +540,7 @@ def _name_jacobian(i, j, restriction, visitor):
sf = SumfactKernel(matrix_sequence=matrix_sequence,
interface=inp,
)
+ from dune.codegen.sumfact.vectorization import attach_vectorization_info
vsf = attach_vectorization_info(sf)
# If this sum factorization kernel was not used in the dry run we
@@ -625,99 +551,12 @@ def _name_jacobian(i, j, restriction, visitor):
# Add a sum factorization kernel that implements the evaluation of
# the basis functions at quadrature points (stage 1)
+ from dune.codegen.sumfact.realization import realize_sum_factorization_kernel
var, _ = realize_sum_factorization_kernel(vsf)
assert(visitor.indices is None)
return prim.Subscript(var, vsf.quadrature_index(sf, visitor))
-def define_jacobian_inverse(name, restriction, visitor):
- """Return jacobian inverse of the geometry mapping (of a cell)
-
- At the moment this only works for geometry mappings of cells and not for
- intersection. We only consider this case as it greatly simplifies code
- generation for unstructured grids by making the grid edge consistent and
- avoiding the face-twist problem.
- """
- # Calculate the jacobian of the geometry mapping using sum factorization
- dim = world_dimension()
- names_jacobian = []
- for j in range(dim):
- for i in range(dim):
- name_jacobian = restricted_name("jacobian_geometry_{}{}".format(i, j), restriction)
- temporary_variable(name_jacobian)
- assignee = prim.Variable(name_jacobian)
- expression = _name_jacobian(i, j, restriction, visitor)
- inames = default_quadrature_inames(visitor) + additional_inames(visitor)
- instruction(assignee=assignee,
- expression=expression,
- within_inames=frozenset(inames),
- )
- names_jacobian.append(name_jacobian)
-
- # The sum factorization kernels from the geometry evaluation of the
- # jacobians will never appear in the expression for the input of
- # stage 3. This way the SumfactCollectMapper will never see them
- # and they will be marked as inactive. Here we explicitly mark the
- # as used.
- basis_sf_kernels(expression.aggregate)
-
- # Calculate the inverse of the jacobian of the geometry mapping and the
- # determinant by calling a c++ function. Note: The result will be column
- # major -> fortran style.
- name_detjac = name_jacobian_determinant(visitor)
- temporary_variable(name_detjac, shape=())
- ftags = ",".join(["f"] * 2)
- temporary_variable(name, shape=(dim, dim), dim_tags=ftags, managed=True)
- include_file('dune/codegen/sumfact/invertgeometry.hh', filetag='operatorfile')
- code = "{} = invert_and_return_determinant({}, {});".format(name_detjac,
- ", ".join(names_jacobian),
- name)
- silenced_warning("read_no_write({})".format(name_detjac))
- inames = default_quadrature_inames(visitor) + additional_inames(visitor)
- return instruction(code=code,
- assignees=frozenset([name, name_detjac]),
- read_variables=frozenset(names_jacobian),
- inames=frozenset(inames),
- )
-
-
-def name_jacobian_inverse(restriction, visitor):
- name = restricted_name("jacobian_inverse", restriction)
- define_jacobian_inverse(name, restriction, visitor)
- return name
-
-
-def pymbolic_jacobian_inverse(i, j, restriction, visitor):
- # Switch between implementations without using the backend mechanism. In
- # the case of constant_transformation_matrix we want to use the pdelab
- # branch, otherwise we want to go into the sumfact implementation instead
- # of pdelab.
- if get_form_option("constant_transformation_matrix"):
- from dune.codegen.pdelab.geometry import name_constant_jacobian_inverse_transposed
- name = name_constant_jacobian_inverse_transposed(restriction)
- # Return jacobian inverse -> (j, i) instead of (i, j)
- return prim.Subscript(prim.Variable(name), (j, i))
- else:
- name = name_jacobian_inverse(restriction, visitor)
- return prim.Subscript(prim.Variable(name), (i, j))
-
-
-def name_jacobian_determinant(visitor):
- name = "detjac"
- return name
-
-
-def pymbolic_jacobian_determinant(visitor):
- if get_form_option("constant_transformation_matrix"):
- from dune.codegen.pdelab.geometry import pymbolic_jacobian_determinant
- return pymbolic_jacobian_determinant()
- else:
- # The calculation of the jacobian determinant happens in this function
- if visitor.measure == 'cell':
- restriction = 0
- else:
- restriction = 1
- name_jacobian_inverse(restriction, visitor)
-
- return prim.Variable(name_jacobian_determinant(visitor))
+def normalize(expr, dim):
+ return prim.Call(prim.Variable("sqrt"), (prim.Sum(tuple(expr[i] * expr[i] for i in range(dim))),))
diff --git a/python/dune/codegen/sumfact/quadrature.py b/python/dune/codegen/sumfact/quadrature.py
index eedd85da732c60fbe3729451ddbe87ee6319f39c..91e99c4c3c0cd333d6e64353ed06ac0c129a7e05 100644
--- a/python/dune/codegen/sumfact/quadrature.py
+++ b/python/dune/codegen/sumfact/quadrature.py
@@ -1,5 +1,4 @@
-from dune.codegen.generation import (backend,
- domain,
+from dune.codegen.generation import (domain,
function_mangler,
get_global_context_value,
globalarg,
@@ -8,6 +7,7 @@ from dune.codegen.generation import (backend,
kernel_cached,
loopy_class_member,
preamble,
+ quadrature_mixin,
temporary_variable,
)
from dune.codegen.sumfact.switch import get_facedir
@@ -20,6 +20,7 @@ from dune.codegen.pdelab.argument import name_accumulation_variable
from dune.codegen.pdelab.geometry import (local_dimension,
world_dimension,
)
+from dune.codegen.pdelab.quadrature import GenericQuadratureMixin
from dune.codegen.options import get_form_option
from dune.codegen.sumfact.switch import get_facedir
from dune.codegen.loopy.target import dtype_floatingpoint
@@ -39,6 +40,85 @@ import loopy as lp
import numpy as np
+@quadrature_mixin("sumfact")
+class SumfactQuadratureMixin(GenericQuadratureMixin):
+ def quadrature_inames(self):
+ info = self.current_info[1]
+ if info is not None:
+ info = info.element
+
+ return _quadrature_inames(info)
+
+ def quadrature_weight(self, o):
+ # Quadrature points per (local) direction
+ dim = local_dimension()
+ local_qps_per_dir = local_quadrature_points_per_direction()
+ local_qps_per_dir_str = '_'.join(map(str, local_qps_per_dir))
+ name = "quad_weights_dim{}_num{}".format(dim, local_qps_per_dir_str)
+
+ # Add a class member
+ loopy_class_member(name,
+ shape=local_qps_per_dir,
+ classtag="operator",
+ dim_tags=",".join(["f"] * dim),
+ managed=True,
+ potentially_vectorized=True,
+ )
+
+ # Precompute it in the constructor
+ inames = constructor_quadrature_inames(dim)
+ assignee = prim.Subscript(prim.Variable(name), tuple(prim.Variable(i) for i in inames))
+ expression = prim.Product(tuple(Subscript(Variable(name_oned_quadrature_weights(local_qps_per_dir[index])),
+ (prim.Variable(iname),)) for index, iname in enumerate(inames)))
+ instruction(assignee=assignee,
+ expression=expression,
+ within_inames=frozenset(inames),
+ kernel="operator",
+ )
+
+ self.quadrature_inames()
+ ret = prim.Subscript(lp.symbolic.TaggedVariable(name, "operator_precomputed"),
+ tuple(prim.Variable(i) for i in self.quadrature_inames()))
+
+ # Multiply the accumulation weight
+ baseweight = 1.0
+ if get_form_option("fastdg"):
+ baseweight = prim.Call(BaseWeight(name_accumulation_variable()), ())
+
+ return prim.Product((baseweight, ret))
+
+ def quadrature_position(self, index=None):
+ assert index is not None
+ dim = local_dimension()
+ qps_per_dir = quadrature_points_per_direction()
+ local_qps_per_dir = local_quadrature_points_per_direction()
+ local_qps_per_dir_str = '_'.join(map(str, local_qps_per_dir))
+ name = "quad_points_dim{}_num{}_localdir{}".format(dim, local_qps_per_dir_str, index)
+
+ loopy_class_member(name,
+ shape=local_qps_per_dir,
+ classtag="operator",
+ dim_tags=",".join(["f"] * dim),
+ managed=True,
+ potentially_vectorized=True,
+ )
+
+ # Precompute it in the constructor
+ inames = constructor_quadrature_inames(dim)
+ assignee = prim.Subscript(prim.Variable(name), tuple(prim.Variable(i) for i in inames))
+ expression = Subscript(Variable(name_oned_quadrature_points(local_qps_per_dir[index])),
+ (prim.Variable(inames[index])))
+ instruction(assignee=assignee,
+ expression=expression,
+ within_inames=frozenset(inames),
+ kernel="operator",
+ )
+
+ return prim.Subscript(lp.symbolic.TaggedVariable(name, "operator_precomputed"),
+ tuple(prim.Variable(i) for i in self.quadrature_inames())
+ )
+
+
class BaseWeight(FunctionIdentifier):
def __init__(self, accumobj):
self.accumobj = accumobj
@@ -60,28 +140,8 @@ def base_weight_function_mangler(target, func, dtypes):
return CallMangleInfo(func.name, (NumpyType(dtype_floatingpoint()),), ())
-def pymbolic_base_weight():
- """ This is the base weight that should be multiplied to the quadrature
- weight. With the fast DG assembler this will handle the weighting of the
- time discretization scheme.
- """
- if get_form_option("fastdg"):
- return prim.Call(BaseWeight(name_accumulation_variable()), ())
- else:
- return 1.0
-
-
-def default_quadrature_inames(visitor):
- info = visitor.current_info[1]
- if info is None:
- element = None
- else:
- element = info.element
- return quadrature_inames(element)
-
-
@kernel_cached
-def quadrature_inames(element):
+def _quadrature_inames(element):
if element is None:
names = tuple("quad_{}".format(d) for d in range(local_dimension()))
else:
@@ -109,157 +169,6 @@ def constructor_quadrature_inames(dim):
return tuple(constructor_quad_iname(name, d, local_qps_per_dir[d]) for d in range(local_dimension()))
-def define_recursive_quadrature_weight(visitor, name, direction):
- info = visitor.current_info[1]
- if info is None:
- element = None
- else:
- element = info.element
- iname = quadrature_inames(element)[direction]
- temporary_variable(name, shape=(), shape_impl=())
-
- qps_per_dir = quadrature_points_per_direction()
- qp_bound = qps_per_dir[direction]
- instruction(expression=Product((recursive_quadrature_weight(direction + 1),
- Subscript(Variable(name_oned_quadrature_weights(qp_bound)),
- (Variable(iname),)
- ))
- ),
- assignee=Variable(name),
- forced_iname_deps=frozenset(quadrature_inames()[direction:]),
- forced_iname_deps_is_final=True,
- tags=frozenset({"quad"}),
- )
-
-
-def recursive_quadrature_weight(visitor, direction=0):
- if direction == local_dimension():
- return pymbolic_base_weight()
- else:
- name = 'weight_{}'.format(direction)
- define_recursive_quadrature_weight(visitor, name, direction)
- return Variable(name)
-
-
-def quadrature_weight(visitor):
- # Return non-precomputed version
- if not get_form_option("precompute_quadrature_info"):
- return recursive_quadrature_weight(visitor)
-
- # Quadrature points per (local) direction
- dim = local_dimension()
- local_qps_per_dir = local_quadrature_points_per_direction()
- local_qps_per_dir_str = '_'.join(map(str, local_qps_per_dir))
- name = "quad_weights_dim{}_num{}".format(dim, local_qps_per_dir_str)
-
- # Add a class member
- loopy_class_member(name,
- shape=local_qps_per_dir,
- classtag="operator",
- dim_tags=",".join(["f"] * dim),
- managed=True,
- potentially_vectorized=True,
- )
-
- # Precompute it in the constructor
- inames = constructor_quadrature_inames(dim)
- assignee = prim.Subscript(prim.Variable(name), tuple(prim.Variable(i) for i in inames))
- expression = prim.Product(tuple(Subscript(Variable(name_oned_quadrature_weights(local_qps_per_dir[index])),
- (prim.Variable(iname),)) for index, iname in enumerate(inames)))
- instruction(assignee=assignee,
- expression=expression,
- within_inames=frozenset(inames),
- kernel="operator",
- )
-
- info = visitor.current_info[1]
- if info is None:
- element = None
- else:
- element = info.element
-
- ret = prim.Subscript(lp.symbolic.TaggedVariable(name, "operator_precomputed"),
- tuple(prim.Variable(i) for i in quadrature_inames(element)))
-
- return prim.Product((pymbolic_base_weight(), ret))
-
-
-def define_quadrature_position(name, local_index):
- # TODO: This whole function looks quite suspicious:
- # - qps_per_dir not used
- # - probably rename index to local_index
- #
- # This whole function is probably not used anymore since we use sum
- # factorization for evaluation of global coordinate, jacobian inverse
- # transposed and determinant of the jacobian of the geometry mapping.
- #
- # Needs to be fixed. Or removed. Added to my to do list ;)
- assert False
-
- qps_per_dir = quadrature_points_per_direction()
- local_qps_per_dir = local_quadrature_points_per_direction()
- qp_bound = local_qps_per_dir[index]
- expression = Subscript(Variable(name_oned_quadrature_points(qp_bound)),
- (Variable(quadrature_inames()[index]),))
- instruction(expression=expression,
- assignee=Subscript(Variable(name), (index,)),
- forced_iname_deps=frozenset(quadrature_inames(element)),
- forced_iname_deps_is_final=True,
- tags=frozenset({"quad"}),
- )
-
-
-def pymbolic_indexed_quadrature_position(local_index, visitor):
- """Return the value of the indexed local quadrature points
-
- Arguments:
- ----------
-
- local_index: Local index to the quadrature point. Local means that on the
- face of a 3D Element the index always will be 0 or 1 but never 2.
- visitor: UFL tree visitor
- """
- # Return the non-precomputed version
- if not get_form_option("precompute_quadrature_info"):
- name = 'pos'
- temporary_variable(name, shape=(local_dimension(),), shape_impl=("fv",))
- define_quadrature_position(name, local_index)
- return prim.Subscript(prim.Variable(name), (local_index,))
-
- assert isinstance(local_index, int)
- dim = local_dimension()
- qps_per_dir = quadrature_points_per_direction()
- local_qps_per_dir = local_quadrature_points_per_direction()
- local_qps_per_dir_str = '_'.join(map(str, local_qps_per_dir))
- name = "quad_points_dim{}_num{}_localdir{}".format(dim, local_qps_per_dir_str, local_index)
-
- loopy_class_member(name,
- shape=local_qps_per_dir,
- classtag="operator",
- dim_tags=",".join(["f"] * dim),
- managed=True,
- potentially_vectorized=True,
- )
-
- # Precompute it in the constructor
- inames = constructor_quadrature_inames(dim)
- assignee = prim.Subscript(prim.Variable(name), tuple(prim.Variable(i) for i in inames))
- expression = Subscript(Variable(name_oned_quadrature_points(local_qps_per_dir[local_index])),
- (prim.Variable(inames[local_index])))
- instruction(assignee=assignee,
- expression=expression,
- within_inames=frozenset(inames),
- kernel="operator",
- )
-
- info = visitor.current_info[1]
- if info is None:
- element = None
- else:
- element = info.element
- return prim.Subscript(lp.symbolic.TaggedVariable(name, "operator_precomputed"), tuple(prim.Variable(i) for i in quadrature_inames(element)))
-
-
def additional_inames(visitor):
"""Return inames for iterating over ansatz space in the case of jacobians
"""
@@ -274,8 +183,7 @@ def additional_inames(visitor):
restriction = Restriction.NONE
if get_global_context_value("integral_type") == "interior_facet":
restriction = Restriction.POSITIVE
- from dune.codegen.sumfact.basis import lfs_inames
- lfs_inames = lfs_inames(element, restriction)
+ lfs_inames = visitor.lfs_inames(element, restriction)
return lfs_inames
else:
return ()
diff --git a/python/dune/codegen/sumfact/realization.py b/python/dune/codegen/sumfact/realization.py
index f14599e5e7009997da267bfd16402d1371e6a6ec..8b585d4ce9a91c8565cb9e6b2b9b495b3734b2e6 100644
--- a/python/dune/codegen/sumfact/realization.py
+++ b/python/dune/codegen/sumfact/realization.py
@@ -17,7 +17,6 @@ from dune.codegen.generation import (barrier,
transform,
)
from dune.codegen.loopy.flatten import flatten_index
-from dune.codegen.pdelab.argument import pymbolic_coefficient
from dune.codegen.pdelab.basis import shape_as_pymbolic
from dune.codegen.pdelab.geometry import world_dimension
from dune.codegen.options import (get_form_option,
diff --git a/python/dune/codegen/sumfact/switch.py b/python/dune/codegen/sumfact/switch.py
index d18ab344b9c1c78d08c85185219a6b2facf95fc7..67cb9df358112e456b4ef9250c3602f1f9aacdec 100644
--- a/python/dune/codegen/sumfact/switch.py
+++ b/python/dune/codegen/sumfact/switch.py
@@ -2,13 +2,11 @@
import csv
-from dune.codegen.generation import (backend,
- get_backend,
- get_global_context_value,
+from dune.codegen.generation import (get_global_context_value,
global_context,
)
from dune.codegen.pdelab.geometry import world_dimension
-from dune.codegen.pdelab.localoperator import generate_kernel
+from dune.codegen.pdelab.localoperator import generate_kernel, generate_kernels_per_integral
from dune.codegen.pdelab.signatures import (assembly_routine_args,
assembly_routine_signature,
kernel_name,
@@ -17,8 +15,7 @@ from dune.codegen.options import get_form_option, get_option, set_form_option
from dune.codegen.cgen.clazz import ClassMember
-@backend(interface="generate_kernels_per_integral", name="sumfact")
-def generate_kernels_per_integral(integrals):
+def sumfact_generate_kernels_per_integral(integrals):
dim = world_dimension()
measure = get_global_context_value("integral_type")
@@ -29,7 +26,7 @@ def generate_kernels_per_integral(integrals):
# Maybe skip sum factorization on boundary integrals
if not get_form_option("sumfact_on_boundary"):
set_form_option("sumfact", False)
- for k in get_backend(interface="generate_kernels_per_integral")(integrals):
+ for k in generate_kernels_per_integral(integrals):
yield k
set_form_option("sumfact", True)
return
@@ -64,7 +61,7 @@ def get_kernel_name(facedir_s=None, facemod_s=None, facedir_n=None, facemod_n=No
def decide_if_kernel_is_necessary(facedir_s, facemod_s, facedir_n, facemod_n):
# If we are not using YaspGrid, all variants need to be realized
- if not get_form_option("diagonal_transformation_matrix"):
+ if get_form_option("geometry_mixins") == "sumfact_multilinear":
# Reduce the variability according to grid info file
if get_option("grid_info") is not None:
filename = get_option("grid_info")
diff --git a/python/dune/codegen/sumfact/symbolic.py b/python/dune/codegen/sumfact/symbolic.py
index b8b5b857728ea063edec3812b2132d210ffe93b3..2768ece8993b01565a2eeacb4311fa976174c236 100644
--- a/python/dune/codegen/sumfact/symbolic.py
+++ b/python/dune/codegen/sumfact/symbolic.py
@@ -14,7 +14,6 @@ from dune.codegen.sumfact.permutation import (flop_cost,
sumfact_cost_permutation_strategy,
sumfact_quadrature_permutation_strategy,
)
-from dune.codegen.sumfact.quadrature import quadrature_inames
from dune.codegen.sumfact.tabulation import BasisTabulationMatrixBase, BasisTabulationMatrixArray
from dune.codegen.loopy.target import dtype_floatingpoint, type_floatingpoint
from dune.codegen.loopy.vcl import ExplicitVCLCast, VCLLowerUpperLoad
@@ -713,16 +712,7 @@ class SumfactKernel(SumfactKernelBase, ImmutableRecord, prim.Variable):
return tuple(mat.quadrature_size for mat in self.matrix_sequence_quadrature_permuted)
def quadrature_index(self, sf, visitor):
- if visitor.current_info[1] is None:
- element = None
- element_index = 0
- else:
- element = visitor.current_info[1].element
- element_index = visitor.current_info[1].element_index
- if isinstance(element, MixedElement):
- element = element.extract_component(element_index)[1]
-
- quad_inames = quadrature_inames(element)
+ quad_inames = visitor.quadrature_inames()
if len(self.matrix_sequence_quadrature_permuted) == local_dimension():
return tuple(prim.Variable(i) for i in quad_inames)
@@ -1039,16 +1029,7 @@ class VectorizedSumfactKernel(SumfactKernelBase, ImmutableRecord, prim.Variable)
return 0
def _quadrature_index(self, sf, visitor):
- if visitor.current_info[1] is None:
- element = None
- element_index = 0
- else:
- element = visitor.current_info[1].element
- element_index = visitor.current_info[1].element_index
- if isinstance(element, MixedElement):
- element = element.extract_component(element_index)[1]
-
- quad_inames = quadrature_inames(element)
+ quad_inames = visitor.quadrature_inames()
index = []
if len(self.matrix_sequence_quadrature_permuted) == local_dimension():
@@ -1077,16 +1058,8 @@ class VectorizedSumfactKernel(SumfactKernelBase, ImmutableRecord, prim.Variable)
return tuple(index)
def vec_index(self, sf, visitor):
- if visitor.current_info[1] is None:
- element = None
- element_index = 0
- else:
- element = visitor.current_info[1].element
- element_index = visitor.current_info[1].element_index
- if isinstance(element, MixedElement):
- element = element.extract_component(element_index)[1]
+ quad_inames = visitor.quadrature_inames()
- quad_inames = quadrature_inames(element)
sliced = 0
if len(self.matrix_sequence_quadrature_permuted) == local_dimension():
for d in range(local_dimension()):
diff --git a/python/dune/codegen/sumfact/vectorization.py b/python/dune/codegen/sumfact/vectorization.py
index e7f36a8f0a741c8e4d262762249c2d005a33ffac..e753652b10b3ac9b1765ee071bce1ccebd15f6b1 100644
--- a/python/dune/codegen/sumfact/vectorization.py
+++ b/python/dune/codegen/sumfact/vectorization.py
@@ -7,12 +7,11 @@ import logging
from dune.codegen.loopy.target import dtype_floatingpoint
from dune.codegen.loopy.vcl import get_vcl_type_size
from dune.codegen.loopy.symbolic import SumfactKernel, VectorizedSumfactKernel
-from dune.codegen.generation import (backend,
- generator_factory,
- get_backend,
+from dune.codegen.generation import (generator_factory,
get_counted_variable,
get_global_context_value,
kernel_cached,
+ retrieve_cache_items,
)
from dune.codegen.pdelab.restriction import (Restriction,
restricted_name,
@@ -151,7 +150,7 @@ class PrimitiveApproximateOpcounter(FlopCounter):
raise NotImplementedError("The class {} should implement a symbolic flopcounter.".format(type(expr)))
-@kernel_cached
+@generator_factory(item_tags=("opcounts",), context_tags="kernel")
def store_operation_count(expr, count):
return count
@@ -185,7 +184,7 @@ def quadrature_penalized_strategy_cost(strat_tuple):
# Determine the number of floating point operations per quadrature point.
# This flop counting is a very crude approximation, but that is totally sufficient here.
- ops_per_qp = sum(i.value for i in store_operation_count._memoize_cache.values())
+ ops_per_qp = sum(i for i in retrieve_cache_items("opcounts"))
# Do the actual scaling.
return float((sf_flops + ops_per_qp * num_qp_new) / (sf_flops + ops_per_qp * num_qp_old)) * cost
diff --git a/python/dune/codegen/ufl/visitor.py b/python/dune/codegen/ufl/visitor.py
index 996024cea6ec966cbb8f5463949c7bb84b2560b2..e55232d26839abb2a8866dc2a57af715156237a7 100644
--- a/python/dune/codegen/ufl/visitor.py
+++ b/python/dune/codegen/ufl/visitor.py
@@ -42,9 +42,7 @@ import numpy as np
class UFL2LoopyVisitor(ModifiedTerminalTracker):
- def __init__(self, interface, measure, subdomain_id):
- self.interface = interface
- self.interface.visitor = self
+ def __init__(self, measure, subdomain_id):
self.measure = measure
self.subdomain_id = subdomain_id
@@ -56,14 +54,14 @@ class UFL2LoopyVisitor(ModifiedTerminalTracker):
return self._call(o, do_predicates)
def accumulate(self, o):
- for info in self.interface.list_accumulation_infos(o, self):
+ for info in self.list_accumulation_infos(o):
self.current_info = info
expr = self._call(o, False)
if expr != 0:
if get_form_option("simplify"):
from dune.codegen.sympy import simplify_pymbolic_expression
expr = simplify_pymbolic_expression(expr)
- self.interface.generate_accumulation_instruction(expr, self)
+ self.generate_accumulation_instruction(expr)
def _call(self, o, do_predicates):
# Reset state variables
@@ -94,9 +92,9 @@ class UFL2LoopyVisitor(ModifiedTerminalTracker):
#
def argument(self, o):
- self.interface.initialize_function_spaces(o, self)
+ self.initialize_function_spaces(o)
# Update the information on where to accumulate this
- info = self.interface.get_accumulation_info(o, self)
+ info = self.get_accumulation_info(o)
if o.number() == 0:
if info != self.current_info[0]:
self.indices = None
@@ -131,9 +129,9 @@ class UFL2LoopyVisitor(ModifiedTerminalTracker):
raise CodegenUFLError("Gradients should have been transformed to reference gradients!!!")
if self.reference_grad:
- return self.interface.pymbolic_reference_gradient(leaf_element, restriction, o.number())
+ return self.implement_reference_gradient(leaf_element, restriction, o.number())
else:
- return self.interface.pymbolic_basis(leaf_element, restriction, o.number())
+ return self.implement_basis(leaf_element, restriction, o.number())
def coefficient(self, o):
# Correct the restriction on boundary integrals
@@ -143,7 +141,7 @@ class UFL2LoopyVisitor(ModifiedTerminalTracker):
# Do something different for trial function and coefficients from jacobian apply
if o.count() == 0 or o.count() == 1:
- self.interface.initialize_function_spaces(o, self)
+ self.initialize_function_spaces(o)
index = None
if isinstance(o.ufl_element(), MixedElement):
@@ -158,14 +156,14 @@ class UFL2LoopyVisitor(ModifiedTerminalTracker):
if self.reference_grad:
if o.count() == 0:
- return self.interface.pymbolic_trialfunction_gradient(o.ufl_element(), restriction, index)
+ return self.implement_trialfunction_gradient(o.ufl_element(), restriction, index)
else:
- return self.interface.pymbolic_apply_function_gradient(o.ufl_element(), restriction, index)
+ return self.implement_apply_function_gradient(o.ufl_element(), restriction, index)
else:
if o.count() == 0:
- return self.interface.pymbolic_trialfunction(o.ufl_element(), restriction, index)
+ return self.implement_trialfunction(o.ufl_element(), restriction, index)
else:
- return self.interface.pymbolic_apply_function(o.ufl_element(), restriction, index)
+ return self.implement_apply_function(o.ufl_element(), restriction, index)
# In this case it represents the time variable
elif o.count() == 2:
@@ -176,7 +174,7 @@ class UFL2LoopyVisitor(ModifiedTerminalTracker):
if self.reference_grad:
raise CodegenUFLError("Coefficient gradients should not be transformed to reference element")
- return self.interface.pymbolic_gridfunction(o, restriction, self.grad)
+ return self.implement_gridfunction(o, restriction, self.grad)
def variable(self, o):
# Right now only scalar varibables are supported
@@ -232,7 +230,7 @@ class UFL2LoopyVisitor(ModifiedTerminalTracker):
if index in self.indexmap:
return self.indexmap[index]
else:
- return Variable(self.interface.name_index(index))
+ raise CodegenUFLError("Index should have been unrolled!")
def multi_index(self, o):
return tuple(self._index_or_fixed_index(i) for i in o)
@@ -248,7 +246,7 @@ class UFL2LoopyVisitor(ModifiedTerminalTracker):
self.indices = None
return self.call(o.ufl_operands[index])
else:
- return self.interface.pymbolic_list_tensor(o)
+ raise CodegenUFLError("Index should have been unrolled!")
def component_tensor(self, o):
assert len(self.indices) == len(o.ufl_operands[1])
@@ -265,10 +263,13 @@ class UFL2LoopyVisitor(ModifiedTerminalTracker):
return ret
def identity(self, o):
- return self.interface.pymbolic_identity(o)
+ i, j = self.indices
+ assert isinstance(i, int) and isinstance(j, int)
+ return 1 if i == j else 0
def inverse(self, o):
- return self.interface.pymbolic_matrix_inverse(o)
+ from dune.codegen.pdelab.tensors import pymbolic_matrix_inverse
+ return pymbolic_matrix_inverse(o, self)
#
# Handlers for arithmetic operators and functions
@@ -421,68 +422,9 @@ class UFL2LoopyVisitor(ModifiedTerminalTracker):
return prim.LogicalNot(self.call(o.ufl_operands[0]))
#
- # Handlers for geometric quantities
+ # NB: Geometry and quadrature related handlers have been moved into configurable mixin classes!
#
- def spatial_coordinate(self, o):
- return self.interface.pymbolic_spatial_coordinate()
-
- def facet_normal(self, o):
- # The normal must be restricted to be well-defined
- assert self.restriction is not Restriction.NONE
-
- # Note: In UFL the jump is defined as: jump(v) = v('+') -
- # v('-'). The corresponding outer unit normal is
- # n=FacetNormal(cell)('+'). In order to be consisten with UFL
- # we need to create the outer unit normal if the restriction
- # is positive.
- if self.restriction == Restriction.POSITIVE:
- return self.interface.pymbolic_unit_outer_normal()
- if self.restriction == Restriction.NEGATIVE:
- # It is highly unnatural to have this generator function,
- # but I do run into subtle trouble with return -1*outer
- # as the indexing into the normal happens only later.
- # Not investing more time into this cornercase right now.
- return self.interface.pymbolic_unit_inner_normal()
-
- def quadrature_weight(self, o):
- return self.interface.pymbolic_quadrature_weight()
-
- def jacobian_determinant(self, o):
- return self.interface.pymbolic_jacobian_determinant()
-
- def jacobian_inverse(self, o):
- restriction = self.restriction
- if self.measure == 'exterior_facet':
- restriction = Restriction.POSITIVE
-
- assert(len(self.indices) == 2)
- i, j = self.indices
- self.indices = None
-
- # Implement diagonal jacobians for unrolled matrices!
- if get_form_option("diagonal_transformation_matrix"):
- if isinstance(i, int) and isinstance(j, int) and i != j:
- return 0
-
- return self.interface.pymbolic_jacobian_inverse(i, j, restriction)
-
- def jacobian(self, o):
- raise NotImplementedError("How did you get Jacobian into your form? We only support JacobianInverse right now. Report!")
-
- def facet_jacobian_determinant(self, o):
- return self.interface.pymbolic_facet_jacobian_determinant()
-
- def cell_volume(self, o):
- restriction = self.restriction
- if self.measure == 'exterior_facet':
- restriction = Restriction.POSITIVE
-
- return self.interface.pymbolic_cell_volume(restriction)
-
- def facet_area(self, o):
- return self.interface.pymbolic_facet_area()
-
#
# Equality/Hashability of the visitor
# In order to pass this visitor into (cached) generator functions, it is beneficial
@@ -491,7 +433,7 @@ class UFL2LoopyVisitor(ModifiedTerminalTracker):
#
def __eq__(self, other):
- return type(self) == type(other)
+ return isinstance(other, UFL2LoopyVisitor)
def __hash__(self):
return 0
diff --git a/python/test/dune/codegen/generation/test_backend.py b/python/test/dune/codegen/generation/test_backend.py
deleted file mode 100644
index 76bda979394f0293f5b9d1a9314a80ce58ba8a66..0000000000000000000000000000000000000000
--- a/python/test/dune/codegen/generation/test_backend.py
+++ /dev/null
@@ -1,35 +0,0 @@
-from dune.codegen.generation import (backend,
- get_backend,
- generator_factory,
- )
-
-
-@backend(interface="foo", name="f1")
-@generator_factory()
-def f1():
- return 1
-
-
-@generator_factory()
-@backend(interface="foo", name="f2")
-def f2():
- return 2
-
-
-@backend(interface="bar", name="f3")
-@generator_factory()
-def f3():
- return 3
-
-
-@generator_factory()
-@backend(interface="bar", name="f4")
-def f4():
- return 4
-
-
-def test_backend():
- assert get_backend(interface="foo", selector=lambda: "f1")() == 1
- assert get_backend(interface="foo", selector=lambda: "f2")() == 2
- assert get_backend(interface="bar", selector=lambda: "f3")() == 3
- assert get_backend(interface="bar", selector=lambda: "f4")() == 4
diff --git a/test/blockstructured/nonlinear/nonlinear.mini b/test/blockstructured/nonlinear/nonlinear.mini
index 5e9835a0311a98e815e7c54e51a02714b3bef9ee..29d8d531c8a15967f2673ff5a1b9326ab832fae9 100644
--- a/test/blockstructured/nonlinear/nonlinear.mini
+++ b/test/blockstructured/nonlinear/nonlinear.mini
@@ -14,3 +14,4 @@ compare_l2errorsquared = 6e-4
[formcompiler.r]
blockstructured = 1
number_of_blocks = 5
+geometry_mixins = blockstructured_equidistant
diff --git a/test/blockstructured/poisson/3d/poisson.mini b/test/blockstructured/poisson/3d/poisson.mini
index e9e34187a447198666d8de5a11e77e54aaa4ebb1..796c8652795d4e6a9715df7d2ac364d12b1638b9 100644
--- a/test/blockstructured/poisson/3d/poisson.mini
+++ b/test/blockstructured/poisson/3d/poisson.mini
@@ -17,4 +17,5 @@ compare_l2errorsquared = 1e-7
numerical_jacobian = 1, 0 | expand num
blockstructured = 1
number_of_blocks = 3
+geometry_mixins = blockstructured_equidistant
diff --git a/test/blockstructured/poisson/poisson.mini b/test/blockstructured/poisson/poisson.mini
index a3a1ff0545f378fe1bed3efefb8911b3645d269f..1d7df30a14d0436c1337734f263a5894d528820f 100644
--- a/test/blockstructured/poisson/poisson.mini
+++ b/test/blockstructured/poisson/poisson.mini
@@ -16,6 +16,7 @@ compare_l2errorsquared = 1e-7
numerical_jacobian = 1, 0 | expand num
blockstructured = 1
number_of_blocks = 3
+geometry_mixins = blockstructured_equidistant
[formcompiler.ufl_variants]
cell = quadrilateral
diff --git a/test/blockstructured/poisson/poisson_grid.mini b/test/blockstructured/poisson/poisson_grid.mini
index c3c256486cb78e7fd74135abc978a5010f8f08db..e7c979d2e8f4da3595367a81ce60e7e27cc04d59 100644
--- a/test/blockstructured/poisson/poisson_grid.mini
+++ b/test/blockstructured/poisson/poisson_grid.mini
@@ -1,9 +1,12 @@
__name = blockstructured_poisson_grid_{__exec_suffix}
-__exec_suffix = vec, nonvec | expand vectorize
+__exec_suffix = {vec_suffix}_{dim_suffix}
-dim = 2
+dim = 2, 3 | expand dimension
-gmshFile = square.msh
+vec_suffix = nonvec, vec | expand vectorized
+dim_suffix = 2d, 3d | expand dimension
+
+gmshFile = square.msh, cube.msh | expand dimension
[wrapper.vtkcompare]
name = {__name}
@@ -11,16 +14,17 @@ reference = poisson_ref
extension = vtu
[formcompiler]
-compare_l2errorsquared = 1e-7
+compare_l2errorsquared = 1e-7, 1e-4 | expand dimension
grid_unstructured = 1
[formcompiler.r]
matrix_free = 1
generate_jacobians = 0
blockstructured = 1
-number_of_blocks = 16
-vectorization_blockstructured = 1, 0 | expand vectorize
+number_of_blocks = 16, 4 | expand dimension
+vectorization_blockstructured = 1, 0 | expand vectorized
+geometry_mixins = blockstructured_multilinear
[formcompiler.ufl_variants]
-cell = quadrilateral
+cell = quadrilateral, hexahedron | expand dimension
degree = 1
\ No newline at end of file
diff --git a/test/blockstructured/poisson/poisson_matrix_free.mini b/test/blockstructured/poisson/poisson_matrix_free.mini
index 3a76151cdc4d6efe9c069e87c67417063e655c57..514ab95b79ea5febb48da6d4a610c9abbb4d5ba4 100644
--- a/test/blockstructured/poisson/poisson_matrix_free.mini
+++ b/test/blockstructured/poisson/poisson_matrix_free.mini
@@ -15,6 +15,7 @@ compare_l2errorsquared = 1e-7
matrix_free = 1
blockstructured = 1
number_of_blocks = 4
+geometry_mixins = blockstructured_equidistant
[formcompiler.ufl_variants]
cell = quadrilateral
diff --git a/test/blockstructured/poisson/poisson_neumann.mini b/test/blockstructured/poisson/poisson_neumann.mini
index 1512f88bc1ba25141e121de948da3aed06933b97..26e75dd9cd63ec92179434d4bec5d9e288cc9566 100644
--- a/test/blockstructured/poisson/poisson_neumann.mini
+++ b/test/blockstructured/poisson/poisson_neumann.mini
@@ -16,3 +16,4 @@ compare_l2errorsquared = 1e-8
numerical_jacobian = 1, 0 | expand num
blockstructured = 1
number_of_blocks = 4
+geometry_mixins = blockstructured_equidistant
diff --git a/test/blockstructured/poisson/poisson_tensor.mini b/test/blockstructured/poisson/poisson_tensor.mini
index 2923263610d8260066e2809090b64ce2d10948c9..6f18a924c4deabebaeb8447abd77076b0e2e127f 100644
--- a/test/blockstructured/poisson/poisson_tensor.mini
+++ b/test/blockstructured/poisson/poisson_tensor.mini
@@ -27,6 +27,7 @@ vectorization_blockstructured = 0, 1 | expand vectorized
generate_jacobians = 0
blockstructured = 1
number_of_blocks = 8, 4 | expand dimension
+geometry_mixins = blockstructured_equidistant, blockstructured_multilinear | expand unstructured
[formcompiler.ufl_variants]
cell = quadrilateral, hexahedron | expand dimension
diff --git a/test/blockstructured/poisson/poisson_unstructured.mini b/test/blockstructured/poisson/poisson_unstructured.mini
index c71c381ada38934e5794e2de566ff895e3db4bc9..73ee2072de494acb4af02e3514a03fd058e5e4ed 100644
--- a/test/blockstructured/poisson/poisson_unstructured.mini
+++ b/test/blockstructured/poisson/poisson_unstructured.mini
@@ -4,7 +4,6 @@ __exec_suffix = {dimname}
dim = 2, 3 | expand dimension
dimname = 2d, 3d | expand dimension
-
lowerleft = 0.0 | repeat {dim}
upperright = 1.0 | repeat {dim}
elements = 8, 2 | expand dimension | repeat {dim}
@@ -19,9 +18,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
diff --git a/test/blockstructured/poisson/poisson_unstructured_vec.mini b/test/blockstructured/poisson/poisson_unstructured_vec.mini
index ddae7714c65011f131ca3a31e85b9197b10901c3..681addc52358037b8daf740c7c4e9aa79a9327e3 100644
--- a/test/blockstructured/poisson/poisson_unstructured_vec.mini
+++ b/test/blockstructured/poisson/poisson_unstructured_vec.mini
@@ -9,7 +9,6 @@ upperright = 1.0 | repeat {dim}
elements = 8, 2 | expand dimension | repeat {dim}
elementType = quadrilateral
-
[wrapper.vtkcompare]
name = {__name}
reference = poisson_ref
@@ -25,6 +24,7 @@ generate_jacobians = 0
blockstructured = 1
number_of_blocks = 16, 8 | expand dimension
vectorization_blockstructured = 1
+geometry_mixins = blockstructured_multilinear
[formcompiler.ufl_variants]
cell = quadrilateral, hexahedron | expand dimension
diff --git a/test/blockstructured/poisson/poisson_vec.mini b/test/blockstructured/poisson/poisson_vec.mini
index 9ec5cb244d9e56ca31eb3fe67ecf1c45d6e88f24..7b132230b5ab54469d227389eea0b9ba94c78f4a 100644
--- a/test/blockstructured/poisson/poisson_vec.mini
+++ b/test/blockstructured/poisson/poisson_vec.mini
@@ -21,6 +21,7 @@ generate_jacobians = 0
blockstructured = 1
number_of_blocks = 16, 8 | expand dimension
vectorization_blockstructured = 1
+geometry_mixins = blockstructured_equidistant
[formcompiler.ufl_variants]
cell = quadrilateral, hexahedron | expand dimension
diff --git a/test/blockstructured/stokes/stokes.mini b/test/blockstructured/stokes/stokes.mini
index 532a4159b6b6019525acdea50da3c18f14e22f9a..251f8a842ab4f6940c4a7015982d00ea3947f2ba 100644
--- a/test/blockstructured/stokes/stokes.mini
+++ b/test/blockstructured/stokes/stokes.mini
@@ -16,3 +16,4 @@ compare_l2errorsquared = 1e-9
numerical_jacobian = 0, 1 | expand num
blockstructured = 1
number_of_blocks = 3
+geometry_mixins = blockstructured_equidistant
diff --git a/test/hyperbolic/linearacoustics.mini b/test/hyperbolic/linearacoustics.mini
index ad1cc95b405ebd4a7624ef8057b99834fb52ccbc..249d3637136504341bccc3f951271c1c68db6e3b 100644
--- a/test/hyperbolic/linearacoustics.mini
+++ b/test/hyperbolic/linearacoustics.mini
@@ -17,7 +17,9 @@ explicit_time_stepping = 1
operators = mass, r
[formcompiler.mass]
+geometry_mixins = equidistant
numerical_jacobian = 1
[formcompiler.r]
+geometry_mixins = equidistant
numerical_jacobian = 1
diff --git a/test/navier-stokes/navierstokes_2d_dg_quadrilateral.mini b/test/navier-stokes/navierstokes_2d_dg_quadrilateral.mini
index 5700107d03d62655e87b50b129aabc5b534134d2..787b3408ce17d074c5c3114ce11cca3e3fddae6f 100644
--- a/test/navier-stokes/navierstokes_2d_dg_quadrilateral.mini
+++ b/test/navier-stokes/navierstokes_2d_dg_quadrilateral.mini
@@ -26,9 +26,11 @@ overlapping = 1
[formcompiler.mass]
numerical_jacobian = 0, 1 | expand num
+geometry_mixins = equidistant
[formcompiler.r]
numerical_jacobian = 0, 1 | expand num
+geometry_mixins = equidistant
[instat]
T = 1e-2
diff --git a/test/navier-stokes/navierstokes_3d_dg_quadrilateral.mini b/test/navier-stokes/navierstokes_3d_dg_quadrilateral.mini
index 8823cb1ba26196a5f0aa01cb78d68ff60625f88b..c0bf23eadb1b5238d94424a424217a1b5ffc616f 100644
--- a/test/navier-stokes/navierstokes_3d_dg_quadrilateral.mini
+++ b/test/navier-stokes/navierstokes_3d_dg_quadrilateral.mini
@@ -21,9 +21,11 @@ operators = mass, r
[formcompiler.mass]
numerical_jacobian = 0, 1 | expand num
+geometry_mixins = equidistant
[formcompiler.r]
numerical_jacobian = 0, 1 | expand num
+geometry_mixins = equidistant
[instat]
T = 5e-2
diff --git a/test/nonlinear/diffusivewave.mini b/test/nonlinear/diffusivewave.mini
index c82f8ed23330982d6caf40f21d71dbb5d96f0588..e02ec95646d5838b95d3be5793adeebc4e060a72 100644
--- a/test/nonlinear/diffusivewave.mini
+++ b/test/nonlinear/diffusivewave.mini
@@ -17,9 +17,11 @@ T = 0.01
operators = mass, poisson
[formcompiler.mass]
+geometry_mixins = equidistant, sumfact_equidistant | expand sf
sumfact = 0, 1 | expand sf
fastdg = 0, 0 | expand sf
[formcompiler.operator]
+geometry_mixins = equidistant, sumfact_equidistant | expand sf
sumfact = 0, 1 | expand sf
fastdg = 0, 0 | expand sf
diff --git a/test/poisson/CMakeLists.txt b/test/poisson/CMakeLists.txt
index 0a4f9897176e1aaa2a2864c23bf97135eb67261b..5eeae1a434b4051470540c44d766360d904b4fdd 100644
--- a/test/poisson/CMakeLists.txt
+++ b/test/poisson/CMakeLists.txt
@@ -75,6 +75,12 @@ dune_add_formcompiler_system_test(UFLFILE poisson_customfunction.ufl
INIFILE poisson_customfunction.mini
)
+# 13. Poisson Test Case: DG quadrilaterals with axiparallel mixin
+dune_add_formcompiler_system_test(UFLFILE poisson_dg_quadrilateral.ufl
+ BASENAME poisson_dg_axiparallel
+ INIFILE poisson_dg_axiparallel.mini
+ )
+
# the reference vtk file
add_executable(poisson_dg_ref reference_main.cc)
set_target_properties(poisson_dg_ref PROPERTIES EXCLUDE_FROM_ALL 1)
diff --git a/test/poisson/dimension-grid-variations/poisson_2d_cg_quadrilateral.mini b/test/poisson/dimension-grid-variations/poisson_2d_cg_quadrilateral.mini
index c5ec6e72a870b19eee8d9b8f881c02c770c65da9..7cc6619beb1761050f15e9a71b5f78078838e36f 100644
--- a/test/poisson/dimension-grid-variations/poisson_2d_cg_quadrilateral.mini
+++ b/test/poisson/dimension-grid-variations/poisson_2d_cg_quadrilateral.mini
@@ -9,3 +9,6 @@ extension = vtu
[formcompiler]
compare_l2errorsquared = 9e-8
+
+[formcompiler.r]
+geometry_mixins = equidistant
diff --git a/test/poisson/dimension-grid-variations/poisson_2d_dg_quadrilateral.mini b/test/poisson/dimension-grid-variations/poisson_2d_dg_quadrilateral.mini
index cce6d681b76f031b3a239a04b7d3cb371bd99a85..a4478b53e9f399fab159dead0010286d0c9cc6a3 100644
--- a/test/poisson/dimension-grid-variations/poisson_2d_dg_quadrilateral.mini
+++ b/test/poisson/dimension-grid-variations/poisson_2d_dg_quadrilateral.mini
@@ -9,3 +9,6 @@ extension = vtu
[formcompiler]
compare_l2errorsquared = 6e-7
+
+[formcompiler.r]
+geometry_mixins = equidistant
diff --git a/test/poisson/dimension-grid-variations/poisson_3d_cg_hexahedron.mini b/test/poisson/dimension-grid-variations/poisson_3d_cg_hexahedron.mini
index 1d4c4e9660cb9f285b3dc51db9c3c525181dd1ea..23bc16ba4b87ef9381aedbf8e3fb773e02dc5ae0 100644
--- a/test/poisson/dimension-grid-variations/poisson_3d_cg_hexahedron.mini
+++ b/test/poisson/dimension-grid-variations/poisson_3d_cg_hexahedron.mini
@@ -9,3 +9,6 @@ extension = vtu
[formcompiler]
compare_l2errorsquared = 9e-5
+
+[formcompiler.r]
+geometry_mixins = equidistant
diff --git a/test/poisson/dimension-grid-variations/poisson_3d_dg_hexahedron.mini b/test/poisson/dimension-grid-variations/poisson_3d_dg_hexahedron.mini
index 1442988525b7b3080bcc0efd51e69ccca4e84f63..3c8ebda08525392d88eb1a900572616cd40052e4 100644
--- a/test/poisson/dimension-grid-variations/poisson_3d_dg_hexahedron.mini
+++ b/test/poisson/dimension-grid-variations/poisson_3d_dg_hexahedron.mini
@@ -9,3 +9,6 @@ extension = vtu
[formcompiler]
compare_l2errorsquared = 1e-4
+
+[formcompiler.r]
+geometry_mixins = equidistant
diff --git a/test/poisson/opcount_poisson_dg_symdiff.mini b/test/poisson/opcount_poisson_dg_symdiff.mini
index d98044df779fa8be1d3c2ee09d4e2023f0d95da9..186c4344f9fed7512cde913433145259a76bf7b5 100644
--- a/test/poisson/opcount_poisson_dg_symdiff.mini
+++ b/test/poisson/opcount_poisson_dg_symdiff.mini
@@ -10,3 +10,4 @@ extension = vtu
[formcompiler]
opcounter = 1
instrumentation_level = 3
+geometry_mixins = equidistant
diff --git a/test/poisson/poisson_dg_axiparallel.mini b/test/poisson/poisson_dg_axiparallel.mini
new file mode 100644
index 0000000000000000000000000000000000000000..1f2490eba8b1f2ce06fc150b5fb71972338c6502
--- /dev/null
+++ b/test/poisson/poisson_dg_axiparallel.mini
@@ -0,0 +1,17 @@
+__name = poisson_dg_axiparallel_{__exec_suffix}
+__exec_suffix = {__exec_symdiff}
+__exec_symdiff = numdiff, symdiff |expand num
+
+extension = 1.0 1.0
+cells = 32 32
+
+[wrapper.vtkcompare]
+name = {__name}
+extension = vtu
+
+[formcompiler]
+compare_l2errorsquared = 7e-7
+
+[formcompiler.r]
+numerical_jacobian = 1, 0 | expand num
+geometry_mixins = axiparallel
diff --git a/test/poisson/poisson_dg_quadrilateral.mini b/test/poisson/poisson_dg_quadrilateral.mini
index dde2e495ab819ed22be6b848518ed71d41cfe30f..3715da8457913ecac94772fd7d045f7d9326ed2b 100644
--- a/test/poisson/poisson_dg_quadrilateral.mini
+++ b/test/poisson/poisson_dg_quadrilateral.mini
@@ -14,3 +14,4 @@ compare_l2errorsquared = 7e-7
[formcompiler.r]
numerical_jacobian = 1, 0 | expand num
+geometry_mixins = equidistant
diff --git a/test/poisson/poisson_dg_tensor.mini b/test/poisson/poisson_dg_tensor.mini
index d696cebc5196bd711ef4e8e2e1371a977be3a31f..e73c1ebff8b536644c7aefce95451854ea97b9c1 100644
--- a/test/poisson/poisson_dg_tensor.mini
+++ b/test/poisson/poisson_dg_tensor.mini
@@ -15,4 +15,5 @@ extension = vtu
compare_l2errorsquared = 4e-6
[formcompiler.r]
-numerical_jacobian = 1, 0 | expand num
\ No newline at end of file
+numerical_jacobian = 1, 0 | expand num
+geometry_mixins = equidistant
\ No newline at end of file
diff --git a/test/stokes/stokes_3d_dg_quadrilateral.mini b/test/stokes/stokes_3d_dg_quadrilateral.mini
index 59396277ede15f6a563d04fb448c4d3b5a445b3b..18bd43c75cb65131496e7ecd2beaff69999612c3 100644
--- a/test/stokes/stokes_3d_dg_quadrilateral.mini
+++ b/test/stokes/stokes_3d_dg_quadrilateral.mini
@@ -13,4 +13,5 @@ extension = vtu
compare_l2errorsquared = 6e-8
[formcompiler.r]
+geometry_mixins = equidistant
numerical_jacobian = 0, 1 | expand num
diff --git a/test/stokes/stokes_3d_quadrilateral.mini b/test/stokes/stokes_3d_quadrilateral.mini
index 17f3d9f510054ff9014830e258edd7840180248e..d008c6363ccd2381848c09976450d6db2e810001 100644
--- a/test/stokes/stokes_3d_quadrilateral.mini
+++ b/test/stokes/stokes_3d_quadrilateral.mini
@@ -14,4 +14,5 @@ extension = vtu
compare_l2errorsquared = 1e-10
[formcompiler.r]
+geometry_mixins = equidistant
numerical_jacobian = 1, 0 | expand num
diff --git a/test/stokes/stokes_dg_quadrilateral.mini b/test/stokes/stokes_dg_quadrilateral.mini
index 78954b12873569589c2874d15858c6121a242eb5..5a1fb4127984eb1b60027247e23efa03e6eacc92 100644
--- a/test/stokes/stokes_dg_quadrilateral.mini
+++ b/test/stokes/stokes_dg_quadrilateral.mini
@@ -13,4 +13,5 @@ extension = vtu
compare_l2errorsquared = 1e-8
[formcompiler.r]
+geometry_mixins = equidistant
numerical_jacobian = 0, 1 | expand num
diff --git a/test/stokes/stokes_quadrilateral.mini b/test/stokes/stokes_quadrilateral.mini
index 6ee36e8220463cbb75764a91b0ae1d2970f28eb3..cdaa55d544dd29c6354190291e12704e610ee064 100644
--- a/test/stokes/stokes_quadrilateral.mini
+++ b/test/stokes/stokes_quadrilateral.mini
@@ -14,4 +14,5 @@ extension = vtu
compare_l2errorsquared = 1e-10
[formcompiler.r]
+geometry_mixins = equidistant
numerical_jacobian = 1, 0 | expand num
diff --git a/test/sumfact/hyperbolic/linearacoustics.mini b/test/sumfact/hyperbolic/linearacoustics.mini
index 1f7b1fa542172fd28e45321d36ab82dc3cbdc1ef..df4cc0521d7019dde13bb60464206045b3ccc3a3 100644
--- a/test/sumfact/hyperbolic/linearacoustics.mini
+++ b/test/sumfact/hyperbolic/linearacoustics.mini
@@ -19,7 +19,9 @@ operators = mass, r
[formcompiler.mass]
numerical_jacobian = 1
sumfact = 1
+geometry_mixins = sumfact_equidistant
[formcompiler.r]
numerical_jacobian = 1
sumfact = 1
+geometry_mixins = sumfact_equidistant
diff --git a/test/sumfact/hyperbolic/lineartransport.mini b/test/sumfact/hyperbolic/lineartransport.mini
index 4ea0034f57922c41dda48dbac96d8ef89ec426b0..07f39ed3a8828bfd9b88f535a2c3dcfb6832cb5e 100644
--- a/test/sumfact/hyperbolic/lineartransport.mini
+++ b/test/sumfact/hyperbolic/lineartransport.mini
@@ -22,7 +22,9 @@ operators = mass, r
[formcompiler.mass]
sumfact = 1
numerical_jacobian = 1
+geometry_mixins = sumfact_equidistant
[formcompiler.r]
sumfact = 1
numerical_jacobian = 1
+geometry_mixins = sumfact_equidistant
diff --git a/test/sumfact/hyperbolic/shallowwater.mini b/test/sumfact/hyperbolic/shallowwater.mini
index 5f63f48bd0bd55537d1055ff685d0998ac28ad15..b9c18a34fc8e4759c8ca08e811c8437368a7f899 100644
--- a/test/sumfact/hyperbolic/shallowwater.mini
+++ b/test/sumfact/hyperbolic/shallowwater.mini
@@ -18,6 +18,8 @@ operators = mass, r
[formcompiler.mass]
sumfact = 1
+geometry_mixins = sumfact_equidistant
[formcompiler.r]
sumfact = 1
+geometry_mixins = sumfact_equidistant
diff --git a/test/sumfact/mass/mass.mini b/test/sumfact/mass/mass.mini
index 6b0e9db8144fe18f7cf5e89f016d228d74ae9173..6b2f8c1b92674b9d62a7c50eb69b0362665a3897 100644
--- a/test/sumfact/mass/mass.mini
+++ b/test/sumfact/mass/mass.mini
@@ -16,3 +16,4 @@ extension = vtu
numerical_jacobian = 1, 0 | expand num
vectorization_quadloop = 1, 0 | expand vec
sumfact = 1
+geometry_mixins = sumfact_equidistant
diff --git a/test/sumfact/mass/mass_3d.mini b/test/sumfact/mass/mass_3d.mini
index fff87d11bb272dc5d5bd57e8f12569581ddc2485..88795d74e1bc1a84faf73412210ca0a24bc5c517 100644
--- a/test/sumfact/mass/mass_3d.mini
+++ b/test/sumfact/mass/mass_3d.mini
@@ -17,6 +17,7 @@ extension = vtu
numerical_jacobian = 1, 0 | expand num
vectorization_quadloop = 1, 0 | expand vec
sumfact = 1
+geometry_mixins = sumfact_equidistant
[formcompiler.ufl_variants]
degree = 3
diff --git a/test/sumfact/mass/sliced.mini b/test/sumfact/mass/sliced.mini
index 17d331901a999c5d2d90a0453dfc524183f6e132..c4f1c56eac78046fbca4a75ff93136e9b68a6d1b 100644
--- a/test/sumfact/mass/sliced.mini
+++ b/test/sumfact/mass/sliced.mini
@@ -15,6 +15,7 @@ vectorization_strategy = explicit
vectorization_horizontal = 1
vectorization_vertical = 4
sumfact = 1
+geometry_mixins = sumfact_equidistant
[formcompiler.ufl_variants]
degree = 3
diff --git a/test/sumfact/poisson/diagonal.mini b/test/sumfact/poisson/diagonal.mini
index 298fadba9554771a7bdf4810b15ad7e925861cbf..067a8076ab0922e42f34e5fd65c0565b6d00bc0c 100644
--- a/test/sumfact/poisson/diagonal.mini
+++ b/test/sumfact/poisson/diagonal.mini
@@ -18,6 +18,7 @@ vectorization_horizontal = 2
vectorization_vertical = 2
quadrature_order = 6, 6, 6
fastdg = 1
+geometry_mixins = sumfact_equidistant
[formcompiler.ufl_variants]
degree = 3
diff --git a/test/sumfact/poisson/facedir-facemod-variation/poisson_dg_3d_facedir_facemod_variation.mini b/test/sumfact/poisson/facedir-facemod-variation/poisson_dg_3d_facedir_facemod_variation.mini
index 45927298669692f8d3675e1d7896d565b922ae53..6384d72b60f76fa53dea540fff49f32e281dd91f 100644
--- a/test/sumfact/poisson/facedir-facemod-variation/poisson_dg_3d_facedir_facemod_variation.mini
+++ b/test/sumfact/poisson/facedir-facemod-variation/poisson_dg_3d_facedir_facemod_variation.mini
@@ -36,6 +36,7 @@ numerical_jacobian = 1, 0 | expand num
sumfact = 0
vectorization_quadloop = 1, 0 | expand quad
vectorization_strategy = model, none | expand grad
+geometry_mixins = sumfact_multilinear
[formcompiler.ufl_variants]
degree = 1, 2 | expand deg
diff --git a/test/sumfact/poisson/facedir-facemod-variation/sumfact_poisson_dg_3d_facedir_facemod_variation.mini b/test/sumfact/poisson/facedir-facemod-variation/sumfact_poisson_dg_3d_facedir_facemod_variation.mini
index 664421e4f91a3893c42783c209068e306b110c7e..58090aeb5cedde40572f0371eb839688645c0af9 100644
--- a/test/sumfact/poisson/facedir-facemod-variation/sumfact_poisson_dg_3d_facedir_facemod_variation.mini
+++ b/test/sumfact/poisson/facedir-facemod-variation/sumfact_poisson_dg_3d_facedir_facemod_variation.mini
@@ -44,6 +44,7 @@ sumfact = 1
sumfact_regular_jacobians = 1
vectorization_quadloop = 1, 0 | expand quad
vectorization_strategy = model, none | expand grad
+geometry_mixins = sumfact_multilinear
[formcompiler.ufl_variants]
degree = 1, 2 | expand deg
diff --git a/test/sumfact/poisson/facedir-facemod-variation/sumfact_poisson_fastdg_3d_facedir_facemod_variation.mini b/test/sumfact/poisson/facedir-facemod-variation/sumfact_poisson_fastdg_3d_facedir_facemod_variation.mini
index f814a9aac356cf416a358ff08816d346138d3d9a..5e0b4c29b2d6f89913c9254a7c5b1070855ca8d0 100644
--- a/test/sumfact/poisson/facedir-facemod-variation/sumfact_poisson_fastdg_3d_facedir_facemod_variation.mini
+++ b/test/sumfact/poisson/facedir-facemod-variation/sumfact_poisson_fastdg_3d_facedir_facemod_variation.mini
@@ -45,6 +45,7 @@ sumfact_regular_jacobians = 1
vectorization_quadloop = 1, 0 | expand quad
vectorization_strategy = model, none | expand grad
fastdg = 1
+geometry_mixins = sumfact_multilinear
[formcompiler.ufl_variants]
degree = 1, 2 | expand deg
diff --git a/test/sumfact/poisson/opcount_poisson_2d_order2.mini b/test/sumfact/poisson/opcount_poisson_2d_order2.mini
index 538189b2e5c14bd745c748fc43e3161c4013cc77..4d76eac6dc2fbf91edf7fe924a7cc50bc30be5df 100644
--- a/test/sumfact/poisson/opcount_poisson_2d_order2.mini
+++ b/test/sumfact/poisson/opcount_poisson_2d_order2.mini
@@ -18,6 +18,7 @@ instrumentation_level = 4
[formcompiler.r]
sumfact = 1
+geometry_mixins = sumfact_equidistant
[formcompiler.ufl_variants]
degree = 2
diff --git a/test/sumfact/poisson/opcount_sumfact_poisson_dg_2d_vec.mini b/test/sumfact/poisson/opcount_sumfact_poisson_dg_2d_vec.mini
index b657c1a2c731e0606093ab3c0e00044afaadd3ae..d12c7896ed219317f8ba895ddf61365011d9d2f6 100644
--- a/test/sumfact/poisson/opcount_sumfact_poisson_dg_2d_vec.mini
+++ b/test/sumfact/poisson/opcount_sumfact_poisson_dg_2d_vec.mini
@@ -15,6 +15,7 @@ instrumentation_level = 4
[formcompiler.r]
sumfact = 1
+geometry_mixins = sumfact_equidistant
[formcompiler.ufl_variants]
degree = 1
diff --git a/test/sumfact/poisson/poisson_2d.mini b/test/sumfact/poisson/poisson_2d.mini
index d9ce1773e575fab9ec9e05057207093c8ea8b747..90cb97db1229c502ada8fa02d190e30394cad4f5 100644
--- a/test/sumfact/poisson/poisson_2d.mini
+++ b/test/sumfact/poisson/poisson_2d.mini
@@ -21,6 +21,7 @@ numerical_jacobian = 1, 0 | expand num
sumfact = 1
vectorization_strategy = explicit, none | expand grad
quadrature_order = 2, 4
+geometry_mixins = sumfact_equidistant
[formcompiler.ufl_variants]
degree = 1, 2 | expand deg
diff --git a/test/sumfact/poisson/poisson_2d_unstructured.mini b/test/sumfact/poisson/poisson_2d_unstructured.mini
index 81a8246ea330ca157edcc61e0ceaf0839ffb4465..ce12ec78f77d940a8b8cec3081a5a9708ece0f56 100644
--- a/test/sumfact/poisson/poisson_2d_unstructured.mini
+++ b/test/sumfact/poisson/poisson_2d_unstructured.mini
@@ -30,6 +30,7 @@ numerical_jacobian = 1, 0 | expand num
sumfact = 1
vectorization_quadloop = 1, 0 | expand quad
vectorization_strategy = model, none | expand grad
+geometry_mixins = sumfact_multilinear
[formcompiler.ufl_variants]
degree = 1, 2 | expand deg
diff --git a/test/sumfact/poisson/poisson_3d.mini b/test/sumfact/poisson/poisson_3d.mini
index bf7b0f3bfdc7823e970694ecead48e6ecb4d7f6a..27ef4608c64394b6bf626b423a6f9d1eab1b7de3 100644
--- a/test/sumfact/poisson/poisson_3d.mini
+++ b/test/sumfact/poisson/poisson_3d.mini
@@ -22,6 +22,7 @@ numerical_jacobian = 1, 0 | expand num
sumfact = 1
vectorization_quadloop = 1, 0 | expand quad
vectorization_strategy = explicit, autotune, none | expand grad
+geometry_mixins = sumfact_equidistant
[formcompiler.ufl_variants]
degree = 1, 2 | expand deg
diff --git a/test/sumfact/poisson/poisson_3d_unstructured.mini b/test/sumfact/poisson/poisson_3d_unstructured.mini
index 8c1ad02cf98fabd8d69e5df48fcd88e529611e03..7d5710fa68d588a96d4ed23aac0ea72dff6276c6 100644
--- a/test/sumfact/poisson/poisson_3d_unstructured.mini
+++ b/test/sumfact/poisson/poisson_3d_unstructured.mini
@@ -30,6 +30,7 @@ numerical_jacobian = 1, 0 | expand num
sumfact = 1
vectorization_quadloop = 1, 0 | expand quad
vectorization_strategy = model, none | expand grad
+geometry_mixins = sumfact_multilinear
[formcompiler.ufl_variants]
degree = 1, 2 | expand deg
diff --git a/test/sumfact/poisson/poisson_dg_2d.mini b/test/sumfact/poisson/poisson_dg_2d.mini
index d6799eac4600300f86303ea63853365698efdf1f..8903d6597f9d71cecdc027cabaae894b8d3f5820 100644
--- a/test/sumfact/poisson/poisson_dg_2d.mini
+++ b/test/sumfact/poisson/poisson_dg_2d.mini
@@ -21,6 +21,7 @@ numerical_jacobian = 1, 0 | expand num
sumfact = 1
vectorization_quadloop = 1, 0 | expand quad
vectorization_strategy = explicit, none | expand grad
+geometry_mixins = sumfact_equidistant
[formcompiler.ufl_variants]
degree = 1, 2 | expand deg
diff --git a/test/sumfact/poisson/poisson_dg_2d_gmsh.mini b/test/sumfact/poisson/poisson_dg_2d_gmsh.mini
index 3437a746a086965c33243c354ca768c944cc9427..8d622b06ff3e6a83b140a6c6db00929fd8fb8c86 100644
--- a/test/sumfact/poisson/poisson_dg_2d_gmsh.mini
+++ b/test/sumfact/poisson/poisson_dg_2d_gmsh.mini
@@ -28,6 +28,7 @@ sumfact = 1
sumfact_regular_jacobians = 1
vectorization_quadloop = 1, 0 | expand quad
vectorization_strategy = model, none | expand grad
+geometry_mixins = sumfact_multilinear
[formcompiler.ufl_variants]
degree = 1, 2 | expand deg
diff --git a/test/sumfact/poisson/poisson_dg_2d_unstructured.mini b/test/sumfact/poisson/poisson_dg_2d_unstructured.mini
index 0d9423ccbc65a91dfac244006f93331703c913a3..eca089ff89316487fda300c4b45947b5b2d0b83b 100644
--- a/test/sumfact/poisson/poisson_dg_2d_unstructured.mini
+++ b/test/sumfact/poisson/poisson_dg_2d_unstructured.mini
@@ -9,7 +9,6 @@ deg_suffix = deg{formcompiler.ufl_variants.degree}
{deg_suffix} == deg1 | exclude
{diff_suffix} == numdiff | exclude
{quadvec_suffix} == quadvec | exclude
-# {gradvec_suffix} == gradvec | exclude
lowerleft = 0.0 0.0
upperright = 1.0 1.0
@@ -30,6 +29,7 @@ numerical_jacobian = 1, 0 | expand num
sumfact = 1
vectorization_quadloop = 1, 0 | expand quad
vectorization_strategy = model, none | expand grad
+geometry_mixins = sumfact_multilinear
[formcompiler.ufl_variants]
degree = 1, 2 | expand deg
diff --git a/test/sumfact/poisson/poisson_dg_3d.mini b/test/sumfact/poisson/poisson_dg_3d.mini
index 91556ff0ca2f70d49e9ac89cfaad206040ebda03..3f0485672102ab0d7fffd7ce4f490c0f86a66cff 100644
--- a/test/sumfact/poisson/poisson_dg_3d.mini
+++ b/test/sumfact/poisson/poisson_dg_3d.mini
@@ -22,6 +22,7 @@ numerical_jacobian = 1, 0 | expand num
sumfact = 1
vectorization_quadloop = 1, 0 | expand quad
vectorization_strategy = explicit, none | expand grad
+geometry_mixins = sumfact_equidistant
[formcompiler.ufl_variants]
degree = 1, 2 | expand deg
diff --git a/test/sumfact/poisson/poisson_dg_3d_gmsh.mini b/test/sumfact/poisson/poisson_dg_3d_gmsh.mini
index 3a78bb50786f23711f67165afd987005f267b4e6..fb9b26801ac504ecac69f58d4c0a1f305ce16ccd 100644
--- a/test/sumfact/poisson/poisson_dg_3d_gmsh.mini
+++ b/test/sumfact/poisson/poisson_dg_3d_gmsh.mini
@@ -29,6 +29,7 @@ sumfact = 1
sumfact_regular_jacobians = 1
vectorization_quadloop = 1, 0 | expand quad
vectorization_strategy = model, none | expand grad
+geometry_mixins = sumfact_multilinear
[formcompiler.ufl_variants]
degree = 1, 2 | expand deg
diff --git a/test/sumfact/poisson/poisson_dg_3d_unstructured.mini b/test/sumfact/poisson/poisson_dg_3d_unstructured.mini
index 7303abfedd817a5df78380c53bc3547ae2b02dc2..ed80ba617c1d91b66d1df9e7b5ecfc508d71fc77 100644
--- a/test/sumfact/poisson/poisson_dg_3d_unstructured.mini
+++ b/test/sumfact/poisson/poisson_dg_3d_unstructured.mini
@@ -9,7 +9,6 @@ deg_suffix = deg{formcompiler.ufl_variants.degree}
{deg_suffix} == deg1 | exclude
{diff_suffix} == numdiff | exclude
{quadvec_suffix} == quadvec | exclude
-# {gradvec_suffix} == gradvec | exclude
lowerleft = 0.0 0.0 0.0
upperright = 1.0 1.0 1.0
@@ -31,6 +30,7 @@ sumfact = 1
sumfact_regular_jacobians = 1
vectorization_quadloop = 1, 0 | expand quad
vectorization_strategy = model, none | expand grad
+geometry_mixins = sumfact_multilinear
[formcompiler.ufl_variants]
degree = 1, 2 | expand deg
diff --git a/test/sumfact/poisson/poisson_dg_tensor.mini b/test/sumfact/poisson/poisson_dg_tensor.mini
index f6884f965eac0a47416af97b26aa849a6be688ad..d0b164bd7f25831da4d2dbc755590d6ee6b5aca8 100644
--- a/test/sumfact/poisson/poisson_dg_tensor.mini
+++ b/test/sumfact/poisson/poisson_dg_tensor.mini
@@ -18,6 +18,7 @@ compare_l2errorsquared = 3e-4
sumfact = 1
vectorization_quadloop = 1, 0 | expand quad
vectorization_strategy = explicit, none | expand grad
+geometry_mixins = sumfact_equidistant
[formcompiler.ufl_variants]
degree = 2
diff --git a/test/sumfact/poisson/poisson_fastdg_2d.mini b/test/sumfact/poisson/poisson_fastdg_2d.mini
index 53012e325a57387003dc42f2af4268f48cbdaa98..3f33180cf2e0e57ffbe360bd92375ac2c67655ee 100644
--- a/test/sumfact/poisson/poisson_fastdg_2d.mini
+++ b/test/sumfact/poisson/poisson_fastdg_2d.mini
@@ -20,6 +20,7 @@ sumfact = 1
vectorization_quadloop = 1, 0 | expand quadvec
vectorization_strategy = explicit, none | expand gradvec
fastdg = 1
+geometry_mixins = sumfact_equidistant
[formcompiler.ufl_variants]
degree = 2
\ No newline at end of file
diff --git a/test/sumfact/poisson/poisson_fastdg_2d_gmsh.mini b/test/sumfact/poisson/poisson_fastdg_2d_gmsh.mini
index 56b69c773b1ce95fa6387f5313af47a94980f8c0..651c14b14555f35fd52429ad8f56bdbd04acce64 100644
--- a/test/sumfact/poisson/poisson_fastdg_2d_gmsh.mini
+++ b/test/sumfact/poisson/poisson_fastdg_2d_gmsh.mini
@@ -27,6 +27,7 @@ fastdg = 1
sumfact_regular_jacobians = 1
vectorization_quadloop = 1, 0 | expand quad
vectorization_strategy = model, none | expand grad
+geometry_mixins = sumfact_multilinear
[formcompiler.ufl_variants]
degree = 1, 2 | expand deg
diff --git a/test/sumfact/poisson/poisson_fastdg_3d.mini b/test/sumfact/poisson/poisson_fastdg_3d.mini
index 46552ce9e960e8a53715016e3de0bd6e770b5425..cc2bd596d77cb6623a1b246f5d1f55ea12d51d76 100644
--- a/test/sumfact/poisson/poisson_fastdg_3d.mini
+++ b/test/sumfact/poisson/poisson_fastdg_3d.mini
@@ -20,6 +20,7 @@ sumfact = 1
vectorization_quadloop = 1, 0 | expand quadvec
vectorization_strategy = explicit, none | expand gradvec
fastdg = 1
+geometry_mixins = sumfact_equidistant
[formcompiler.ufl_variants]
degree = 2
diff --git a/test/sumfact/poisson/poisson_fastdg_3d_gmsh.mini b/test/sumfact/poisson/poisson_fastdg_3d_gmsh.mini
index 438b8d009c56e42e9f65bd4b708e8f6acd1f6e06..ada6f4aed8cc1769b4793d8f1c6f7ee55ba2cbb3 100644
--- a/test/sumfact/poisson/poisson_fastdg_3d_gmsh.mini
+++ b/test/sumfact/poisson/poisson_fastdg_3d_gmsh.mini
@@ -7,7 +7,6 @@ deg_suffix = deg{formcompiler.ufl_variants.degree}
{deg_suffix} == deg1 | exclude
{quadvec_suffix} == quadvec | exclude
-# {gradvec_suffix} == gradvec | exclude
gmshFile = cube_hexa.msh
@@ -27,6 +26,7 @@ fastdg = 1
sumfact_regular_jacobians = 1
vectorization_quadloop = 1, 0 | expand quad
vectorization_strategy = model, none | expand grad
+geometry_mixins = sumfact_multilinear
[formcompiler.ufl_variants]
degree = 1, 2 | expand deg
diff --git a/test/sumfact/poisson/sliced.mini b/test/sumfact/poisson/sliced.mini
index 55b6fcf7df3b0ee33f06bb97da0e8b555104c161..29a8a6bdb02d2173bd0f327f49af32a22062cb90 100644
--- a/test/sumfact/poisson/sliced.mini
+++ b/test/sumfact/poisson/sliced.mini
@@ -17,6 +17,7 @@ vectorization_strategy = explicit
vectorization_horizontal = 1
vectorization_vertical = 4
quadrature_order = 6
+geometry_mixins = sumfact_equidistant
[formcompiler.ufl_variants]
degree = 3
diff --git a/test/sumfact/stokes/stokes.mini b/test/sumfact/stokes/stokes.mini
index cfb89ec57d5504fca215ddc08d862e5e2484fc19..91b80c9971051768bed35342e47973380a0b4949 100644
--- a/test/sumfact/stokes/stokes.mini
+++ b/test/sumfact/stokes/stokes.mini
@@ -18,3 +18,4 @@ compare_l2errorsquared = 1e-12
numerical_jacobian = 1, 0 | expand num
vectorization_quadloop = 1, 0 | expand quad
sumfact = 1
+geometry_mixins = sumfact_equidistant
diff --git a/test/sumfact/stokes/stokes_3d_dg.mini b/test/sumfact/stokes/stokes_3d_dg.mini
index 7fb1b22ed1de823afff15f519beb4aa474a3ec11..2f19c49b5309d6e2ea28d966096b169a64d2693f 100644
--- a/test/sumfact/stokes/stokes_3d_dg.mini
+++ b/test/sumfact/stokes/stokes_3d_dg.mini
@@ -17,4 +17,5 @@ compare_l2errorsquared = 1e-10
[formcompiler.r]
numerical_jacobian = 0
sumfact = 1
-fastdg = 1, 0 | expand fastdg
\ No newline at end of file
+fastdg = 1, 0 | expand fastdg
+geometry_mixins = sumfact_equidistant
\ No newline at end of file
diff --git a/test/sumfact/stokes/stokes_dg.mini b/test/sumfact/stokes/stokes_dg.mini
index f34f23422ae888f7ea5d3085b41b87c2ac929346..c58cbc91e7207fba28a45d86de1aa0d892a3dea5 100644
--- a/test/sumfact/stokes/stokes_dg.mini
+++ b/test/sumfact/stokes/stokes_dg.mini
@@ -19,5 +19,6 @@ compare_l2errorsquared = 1e-8
numerical_jacobian = 0, 1 | expand num
sumfact = 1
fastdg = 1, 0 | expand fastdg
+geometry_mixins = sumfact_equidistant
{formcompiler.r.fastdg} == 1 and {formcompiler.r.numerical_jacobian} == 1 | exclude
\ No newline at end of file