diff --git a/python/dune/codegen/loopy/transformations/vectorize_quad.py b/python/dune/codegen/loopy/transformations/vectorize_quad.py
index 77686ee2355fc4828322437aeb3edcbaeb66770b..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
@@ -326,7 +325,7 @@ def vectorize_quadrature_loop(knl):
# Loop over the quadrature loops that exist in the kernel.
# 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]
+ 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:
diff --git a/python/dune/codegen/options.py b/python/dune/codegen/options.py
index f89281b888cb951ddf688c684f779c87dccc8643..e850df54c5da547f53c2798ba322025d0f4b635c 100644
--- a/python/dune/codegen/options.py
+++ b/python/dune/codegen/options.py
@@ -101,7 +101,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")
@@ -111,6 +110,9 @@ 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")
+ quadrature_mixins = CodegenOption(default="generic", helpstr="A comma separated list of mixin identifiers to use for quadrature. Currently implemented: generic")
+ basis_mixins = CodegenOption(default="generic", helpstr="A comma separated list of mixin identifiers to use for basis function evaluation. Currently implemented: generic")
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")
@@ -178,7 +180,13 @@ 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",
+ constant_transformation_matrix=True,
+ diagonal_transformation_matrix=True,
+ )
+ #TODO Remove the trafo matrix ones!
if opt.numerical_jacobian:
opt = opt.copy(generate_jacobians=False, generate_jacobian_apply=False)
diff --git a/python/dune/codegen/pdelab/basis.py b/python/dune/codegen/pdelab/basis.py
index c88ab605e0672cd96e43b0ed175c9697826b8db1..ec03a182ecd476908259c9e289292391c59b019d 100644
--- a/python/dune/codegen/pdelab/basis.py
+++ b/python/dune/codegen/pdelab/basis.py
@@ -48,25 +48,25 @@ from loopy import Reduction
@basis_mixin("base")
class BasisMixinBase(object):
- def lfs_inames(self, *args):
+ def lfs_inames(self, element, restriction, number):
raise NotImplementedError("Basis Mixins should implement local function space inames")
- def implement_basis(self, *args):
+ def implement_basis(self, element, restriction, number):
raise NotImplementedError("Basis Mixins should implement the basis")
- def implement_reference_gradient(self, *args):
+ def implement_reference_gradient(self, element, restriction, number, index):
raise NotImplementedError("Basis Mixins should implement the basis gradient")
- def implement_trialfunction(self, *args):
+ def implement_trialfunction(self, element, restriction, index):
raise NotImplementedError("Basis Mixins should implement trial function evaluation")
- def implement_trialfunction_gradient(self, *args):
+ def implement_trialfunction_gradient(self, element, restriction, index):
raise NotImplementedError("Basis Mixins should implement trial function gradient evaluation")
- def implement_apply_function(self, *args):
+ def implement_apply_function(self, element, restriction, index):
raise NotImplementedError("Basis Mixins should implement linearization point evaluation")
- def implement_apply_function_gradient(self, *args):
+ def implement_apply_function_gradient(self, element, restriction, index):
raise NotImplementedError("Basis Mixins should implement linearization point gradient evaluation")
diff --git a/python/dune/codegen/pdelab/geometry.py b/python/dune/codegen/pdelab/geometry.py
index e1f2bd4b328e32ec4750d8a5b178472ba8fe357c..5bdc2fac640e5c6eb3d2363ad2d972c6560cafd2 100644
--- a/python/dune/codegen/pdelab/geometry.py
+++ b/python/dune/codegen/pdelab/geometry.py
@@ -217,13 +217,20 @@ class AxiparallelGeometryMixin(GenericPDELabGeometryMixin):
def jacobian_inverse(self, o):
i, j = self.indices
- self.indices = None
assert isinstance(i, int) and isinstance(j, int)
if i != j:
+ self.indices = None
return 0
- return GenericPDELabGeometryMixin.jacobian_determinant(self, o)
+ jac = GenericPDELabGeometryMixin.jacobian_inverse(self, o)
+
+ name = get_pymbolic_basename(jac)
+ dim = world_dimension()
+ globalarg(name, shape=(dim, dim), managed=False)
+
+ return jac
+
def facet_area(self, o):
# This is not 100% correct, but in practice structured grid implementations are not
@@ -262,6 +269,7 @@ class EquidistantGeometryMixin(AxiparallelGeometryMixin):
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")
@@ -274,7 +282,6 @@ class EquidistantGeometryMixin(AxiparallelGeometryMixin):
@class_member(classtag="operator")
def define_jacobian_inverse_transposed(self, name, restriction):
dim = world_dimension()
- globalarg(name, shape=(dim, dim), managed=False)
self._define_jacobian_inverse_transposed_eval(name)
from dune.codegen.pdelab.localoperator import lop_template_ansatz_gfs
gfst = lop_template_ansatz_gfs()
diff --git a/python/dune/codegen/pdelab/localoperator.py b/python/dune/codegen/pdelab/localoperator.py
index ef6687033f3a12f7bbbb574daea10b75d34c4495..59a669209edb9ec39767ac4fbfff7e92debc61aa 100644
--- a/python/dune/codegen/pdelab/localoperator.py
+++ b/python/dune/codegen/pdelab/localoperator.py
@@ -452,15 +452,15 @@ def get_visitor(measure, subdomain_id):
# Construct the visitor taking into account geometry mixins
from dune.codegen.ufl.visitor import UFL2LoopyVisitor
- mixins = get_option("geometry_mixins").split(",")
+ 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_option("quadrature_mixins").split(",")
+ 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_option("basis_mixins").split(",")
+ mixins = get_form_option("basis_mixins").split(",")
VisitorType = construct_from_mixins(base=VisitorType, mixins=mixins, mixintype="basis", name="UFLVisitor")
return VisitorType(interface, measure, subdomain_id)
diff --git a/python/dune/codegen/pdelab/quadrature.py b/python/dune/codegen/pdelab/quadrature.py
index 8066423b641c7ca32911d1ead28b459382b228b2..71972c2c92cdb374045e3f8ce00355c7fdaab335 100644
--- a/python/dune/codegen/pdelab/quadrature.py
+++ b/python/dune/codegen/pdelab/quadrature.py
@@ -31,12 +31,9 @@ class QuadratureMixinBase(object):
def quadrature_inames(self):
raise NotImplementedError("Quadrature mixins should provide the quadrature inames!")
- def quadrature_position(self):
+ def quadrature_position(self, index=None):
raise NotImplementedError("Quadrature mixins should provide the quadrature position")
- def quadrature_position_in_cell(self, restriction):
- raise NotImplementedError("Quadrature mixins hsould provide the embedding of the quadrature rule into neighboring cells")
-
@quadrature_mixin("generic")
class GenericQuadratureMixin(QuadratureMixinBase):
@@ -75,7 +72,7 @@ class GenericQuadratureMixin(QuadratureMixinBase):
def quadrature_inames(self):
return (quadrature_iname(),)
- def quadrature_position(self):
+ def quadrature_position(self, index=None):
from dune.codegen.pdelab.geometry import local_dimension
dim = local_dimension()
order = quadrature_order()
@@ -85,7 +82,12 @@ class GenericQuadratureMixin(QuadratureMixinBase):
globalarg(name, shape=shape, managed=False)
self.define_quadrature_points(name)
- return prim.Subscript(prim.Variable(name), tuple(prim.Variable(i) for i in self.quadrature_inames()))
+ 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):
diff --git a/python/dune/codegen/sumfact/__init__.py b/python/dune/codegen/sumfact/__init__.py
index 5dc7f3b159bebec7d9a667cca1966174cac84fd8..bfe02a45337c049b6297146163cd77c1b34db877 100644
--- a/python/dune/codegen/sumfact/__init__.py
+++ b/python/dune/codegen/sumfact/__init__.py
@@ -1,18 +1,10 @@
+import dune.codegen.sumfact.geometry
+
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.accumulation
import dune.codegen.sumfact.switch
@@ -32,42 +24,3 @@ class SumFactInterface(PDELabInterface):
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 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
diff --git a/python/dune/codegen/sumfact/accumulation.py b/python/dune/codegen/sumfact/accumulation.py
index 5bf6703f3133bc4cdae6716480228c651e927976..7a3caf3e1958ab345b23a1348660ad046f25caa3 100644
--- a/python/dune/codegen/sumfact/accumulation.py
+++ b/python/dune/codegen/sumfact/accumulation.py
@@ -40,7 +40,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
@@ -120,8 +119,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(self, sf, result, insn_dep, inames=None, additional_inames=()):
trial_leaf_element = get_leaf(self.trial_element, self.trial_element_index) if self.trial_element is not None else None
@@ -305,10 +304,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:
@@ -480,7 +479,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"]),
)
@@ -495,9 +494,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 ce29589a1e3a8d7739519341578813fd0946603e..a6fdbb301f52f9ebbbc57bdff6a40e1cda511800 100644
--- a/python/dune/codegen/sumfact/basis.py
+++ b/python/dune/codegen/sumfact/basis.py
@@ -6,6 +6,7 @@ multiplication with the test function is part of the sum factorization kernel.
import itertools
from dune.codegen.generation import (backend,
+ basis_mixin,
domain,
get_counted_variable,
get_counter,
@@ -24,11 +25,11 @@ from dune.codegen.sumfact.tabulation import (basis_functions_per_direction,
name_polynomials,
polynomial_degree,
)
-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,
)
@@ -50,6 +51,183 @@ 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(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,
coeff_func=None,
@@ -161,97 +339,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(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(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))
@@ -259,121 +346,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 f7f1fcdc92854195e1b8536b4d8728283a033036..e4c23caf215adc0b07f76c112c35ae0c347827a1 100644
--- a/python/dune/codegen/sumfact/geometry.py
+++ b/python/dune/codegen/sumfact/geometry.py
@@ -34,12 +34,8 @@ from dune.codegen.sumfact.accumulation import basis_sf_kernels
from dune.codegen.sumfact.basis import construct_basis_matrix_sequence
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.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
@@ -65,6 +61,7 @@ class SumfactMultiLinearGeometryMixin(GenericPDELabGeometryMixin):
class SumfactAxiParallelGeometryMixin(AxiparallelGeometryMixin):
def facet_normal(self, o):
i, = self.indices
+ self.indices = None
assert isinstance(i, int)
# Use facemod_s and facedir_s
@@ -78,7 +75,7 @@ class SumfactAxiParallelGeometryMixin(AxiparallelGeometryMixin):
@geometry_mixin("sumfact_equidistant")
-class SumfactAxiParallelGeometryMixin(EquidistantGeometryMixin, SumfactAxiParallelGeometryMixin):
+class SumfactEqudistantGeometryMixin(EquidistantGeometryMixin, SumfactAxiParallelGeometryMixin):
def facet_jacobian_determinant(self, o):
name = "fdetjac"
self.define_facet_jacobian_determinant(name)
@@ -107,6 +104,37 @@ class SumfactAxiParallelGeometryMixin(EquidistantGeometryMixin, SumfactAxiParall
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
+ 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 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
def global_corner_iname(restriction):
@@ -145,6 +173,7 @@ class GeoCornersInput(SumfactKernelInterfaceBase, ImmutableRecord):
def realize(self, sf, insn_dep, index=0):
# Note: world_dimension, since we only do evaluation of cell geometry mappings
+ from dune.codegen.sumfact.realization import name_buffer_storage
name = "input_{}".format(sf.buffer)
temporary_variable(name,
shape=(2 ** world_dimension(), sf.vector_width),
@@ -209,6 +238,8 @@ def pymbolic_spatial_coordinate_multilinear(do_predicates, 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
@@ -219,6 +250,7 @@ def pymbolic_spatial_coordinate_multilinear(do_predicates, 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)
# The results of this function is already the right component of the
@@ -548,6 +580,8 @@ 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
@@ -558,6 +592,7 @@ 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)
diff --git a/python/dune/codegen/sumfact/quadrature.py b/python/dune/codegen/sumfact/quadrature.py
index d5181a282dea1ed14054069e0e0193dbd8a11308..b68da208360fe12aa7b85894759b9dfa53f4136f 100644
--- a/python/dune/codegen/sumfact/quadrature.py
+++ b/python/dune/codegen/sumfact/quadrature.py
@@ -44,11 +44,80 @@ import numpy as np
@quadrature_mixin("sumfact")
class SumfactQuadratureMixin(GenericQuadratureMixin):
def quadrature_inames(self):
- element = self.trial_info
- if self.trial_info is not None:
- element = element.element
-
- return quadrature_inames(element)
+ 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):
@@ -72,17 +141,6 @@ 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:
@@ -93,7 +151,7 @@ def default_quadrature_inames(visitor):
@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:
@@ -153,49 +211,6 @@ def recursive_quadrature_weight(visitor, direction=0):
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
@@ -221,57 +236,6 @@ def define_quadrature_position(name, local_index):
)
-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
"""
diff --git a/python/dune/codegen/sumfact/symbolic.py b/python/dune/codegen/sumfact/symbolic.py
index 24f924555b0573d895ed7f574aa6d744de44f6f3..05dc74b12b296f6f627aea56a3832110c292ddd2 100644
--- a/python/dune/codegen/sumfact/symbolic.py
+++ b/python/dune/codegen/sumfact/symbolic.py
@@ -7,7 +7,6 @@ from dune.codegen.generation import (get_counted_variable,
)
from dune.codegen.pdelab.geometry import local_dimension, world_dimension
from dune.codegen.sumfact.permutation import permute_forward, 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
@@ -462,16 +461,7 @@ class SumfactKernel(SumfactKernelBase, ImmutableRecord, prim.Variable):
return tuple(mat.quadrature_size for mat in self.permuted_matrix_sequence)
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.permuted_matrix_sequence) == local_dimension():
return tuple(prim.Variable(i) for i in quad_inames)
@@ -753,16 +743,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) == local_dimension():
@@ -791,16 +772,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(sf.permuted_matrix_sequence) == local_dimension():
for d in range(local_dimension()):
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/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 cb10c5bf36232f295b1e4cf665f6b4ccf556ef75..90cb97db1229c502ada8fa02d190e30394cad4f5 100644
--- a/test/sumfact/poisson/poisson_2d.mini
+++ b/test/sumfact/poisson/poisson_2d.mini
@@ -21,7 +21,7 @@ numerical_jacobian = 1, 0 | expand num
sumfact = 1
vectorization_strategy = explicit, none | expand grad
quadrature_order = 2, 4
-geometry = sumfact_axiparallel
+geometry_mixins = sumfact_equidistant
[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_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_3d.mini b/test/sumfact/poisson/poisson_dg_3d.mini
index f0b4ef26f73509e9dea47a1cca366c83a51bfb93..db9dc2e3c2f719c17d0b619abc946664103fe8b5 100644
--- a/test/sumfact/poisson/poisson_dg_3d.mini
+++ b/test/sumfact/poisson/poisson_dg_3d.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_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_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