From 5871149ae3941d656cd44d743425a7da1a3f1cb5 Mon Sep 17 00:00:00 2001
From: Dominic Kempf <dominic.kempf@iwr.uni-heidelberg.de>
Date: Wed, 5 Dec 2018 09:56:44 +0100
Subject: [PATCH] Generic PDELab stuff is working again
---
.../dune/codegen/blockstructured/__init__.py | 34 --
.../codegen/blockstructured/accumulation.py | 2 +-
python/dune/codegen/blockstructured/basis.py | 1 -
.../codegen/blockstructured/quadrature.py | 11 +-
python/dune/codegen/blockstructured/tools.py | 4 +-
python/dune/codegen/generation/__init__.py | 4 +-
python/dune/codegen/generation/mixins.py | 15 +-
python/dune/codegen/options.py | 4 +-
python/dune/codegen/pdelab/__init__.py | 46 +--
python/dune/codegen/pdelab/argument.py | 35 --
python/dune/codegen/pdelab/basis.py | 310 +++++++++++-------
python/dune/codegen/pdelab/geometry.py | 138 ++++----
python/dune/codegen/pdelab/localoperator.py | 14 +-
python/dune/codegen/pdelab/quadrature.py | 206 +++++-------
python/dune/codegen/sumfact/__init__.py | 5 -
python/dune/codegen/sumfact/quadrature.py | 12 +
python/dune/codegen/ufl/visitor.py | 23 +-
test/poisson/poisson_dg_quadrilateral.mini | 1 +
18 files changed, 404 insertions(+), 461 deletions(-)
diff --git a/python/dune/codegen/blockstructured/__init__.py b/python/dune/codegen/blockstructured/__init__.py
index 138b1146..fe3eb1f7 100644
--- a/python/dune/codegen/blockstructured/__init__.py
+++ b/python/dune/codegen/blockstructured/__init__.py
@@ -5,7 +5,6 @@ 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)
@@ -27,36 +26,3 @@ class BlockStructuredInterface(PDELabInterface):
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 469e69b4..0f54fbe9 100644
--- a/python/dune/codegen/blockstructured/accumulation.py
+++ b/python/dune/codegen/blockstructured/accumulation.py
@@ -52,7 +52,7 @@ def generate_accumulation_instruction(expr, visitor):
predicates = boundary_predicates(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
diff --git a/python/dune/codegen/blockstructured/basis.py b/python/dune/codegen/blockstructured/basis.py
index dcaf258b..3477a182 100644
--- a/python/dune/codegen/blockstructured/basis.py
+++ b/python/dune/codegen/blockstructured/basis.py
@@ -17,7 +17,6 @@ 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.restriction import restricted_name
from dune.codegen.blockstructured.spaces import lfs_inames
diff --git a/python/dune/codegen/blockstructured/quadrature.py b/python/dune/codegen/blockstructured/quadrature.py
index 4c78d733..64282c68 100644
--- a/python/dune/codegen/blockstructured/quadrature.py
+++ b/python/dune/codegen/blockstructured/quadrature.py
@@ -1,10 +1,17 @@
-from dune.codegen.generation import (backend)
-from dune.codegen.pdelab.quadrature import (name_quadrature_points,
+from dune.codegen.generation import (backend,
+ quadrature_mixin,
+ )
+from dune.codegen.pdelab.quadrature import (GenericQuadratureMixin,
quadrature_iname)
from dune.codegen.blockstructured.geometry import name_point_in_macro
import pymbolic.primitives as prim
+@quadrature_mixin("blockstructured")
+class BlockstructuredQuadratureMixin(GenericQuadratureMixin):
+ pass
+
+
@backend(interface="quad_pos", name='blockstructured')
def pymbolic_quadrature_position():
quad_points = name_quadrature_points()
diff --git a/python/dune/codegen/blockstructured/tools.py b/python/dune/codegen/blockstructured/tools.py
index 3eafcae5..e3150804 100644
--- a/python/dune/codegen/blockstructured/tools.py
+++ b/python/dune/codegen/blockstructured/tools.py
@@ -6,10 +6,8 @@ from dune.codegen.generation import (iname,
instruction)
from dune.codegen.pdelab.geometry import (local_dimension,
world_dimension,
- name_localcenter,
- pymbolic_in_cell_coordinates)
+ name_localcenter,)
-from dune.codegen.pdelab.quadrature import quadrature_inames
from dune.codegen.generation.counter import get_counted_variable
from dune.codegen.options import get_form_option
import pymbolic.primitives as prim
diff --git a/python/dune/codegen/generation/__init__.py b/python/dune/codegen/generation/__init__.py
index cf729760..ac8de9f4 100644
--- a/python/dune/codegen/generation/__init__.py
+++ b/python/dune/codegen/generation/__init__.py
@@ -59,6 +59,8 @@ from dune.codegen.generation.context import (cache_restoring,
get_global_context_value,
)
-from dune.codegen.generation.mixins import (construct_from_mixins,
+from dune.codegen.generation.mixins import (basis_mixin,
+ construct_from_mixins,
geometry_mixin,
+ quadrature_mixin,
)
diff --git a/python/dune/codegen/generation/mixins.py b/python/dune/codegen/generation/mixins.py
index a8a7e72b..ce267972 100644
--- a/python/dune/codegen/generation/mixins.py
+++ b/python/dune/codegen/generation/mixins.py
@@ -1,9 +1,12 @@
""" The infrastructure for registered mixin classes """
+from functools import partial
+
+
_mixin_registry = {}
-def mixin_base(name, mixintype):
+def mixin_base(mixintype, name):
_mixin_registry.setdefault(mixintype, {})
def _dec(cls):
@@ -13,10 +16,12 @@ def mixin_base(name, mixintype):
return _dec
-def geometry_mixin(name):
- return mixin_base(name, "geometry")
-
-
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")
\ No newline at end of file
diff --git a/python/dune/codegen/options.py b/python/dune/codegen/options.py
index fcb5d9ed..8526339b 100644
--- a/python/dune/codegen/options.py
+++ b/python/dune/codegen/options.py
@@ -55,7 +55,9 @@ class CodegenGlobalOptionsArray(ImmutableRecord):
target_name = CodegenOption(default=None, helpstr="The target name from CMake")
operator_to_build = CodegenOption(default=None, helpstr="The operators from the list that is about to be build now. CMake sets this one!!!")
debug_interpolate_input = CodegenOption(default=False, helpstr="Should the input for printresidual and printmatix be interpolated (instead of random input).")
- geometry = CodegenOption(default="generic", helpstr="A comma separated list of mixin identifiers to use for geometries. Currently implemented mixins: generic, axiparallel, equidistant")
+ 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")
# Arguments that are mainly to be set by logic depending on other options
max_vector_width = CodegenOption(default=256, helpstr=None)
diff --git a/python/dune/codegen/pdelab/__init__.py b/python/dune/codegen/pdelab/__init__.py
index 42915eac..604cacf0 100644
--- a/python/dune/codegen/pdelab/__init__.py
+++ b/python/dune/codegen/pdelab/__init__.py
@@ -4,21 +4,9 @@
# 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.argument import pymbolic_coefficient
from dune.codegen.pdelab.function import pymbolic_gridfunction
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,
@@ -70,27 +58,6 @@ class PDELabInterface(object):
#
# 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)
@@ -106,14 +73,3 @@ class PDELabInterface(object):
def pymbolic_matrix_inverse(self, o):
return pymbolic_matrix_inverse(o, self.visitor)
-
- #
- # 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/argument.py b/python/dune/codegen/pdelab/argument.py
index 8972c32c..f9c6c56b 100644
--- a/python/dune/codegen/pdelab/argument.py
+++ b/python/dune/codegen/pdelab/argument.py
@@ -15,9 +15,6 @@ from dune.codegen.generation import (domain,
)
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 +86,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
diff --git a/python/dune/codegen/pdelab/basis.py b/python/dune/codegen/pdelab/basis.py
index 6b6afb3b..c88ab605 100644
--- a/python/dune/codegen/pdelab/basis.py
+++ b/python/dune/codegen/pdelab/basis.py
@@ -1,6 +1,7 @@
""" Generators for basis evaluations """
from dune.codegen.generation import (backend,
+ basis_mixin,
class_member,
get_backend,
include_file,
@@ -12,6 +13,7 @@ from dune.codegen.generation import (backend,
from dune.codegen.options import (option_switch,
get_form_option,
)
+from dune.codegen.pdelab.argument import name_coefficientcontainer
from dune.codegen.pdelab.spaces import (lfs_iname,
lfs_inames,
name_leaf_lfs,
@@ -22,7 +24,6 @@ from dune.codegen.pdelab.spaces import (lfs_iname,
)
from dune.codegen.pdelab.geometry import (component_iname,
world_dimension,
- to_cell_coordinates,
name_cell,
)
from dune.codegen.pdelab.localoperator import (lop_template_ansatz_gfs,
@@ -38,12 +39,195 @@ 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
+@basis_mixin("base")
+class BasisMixinBase(object):
+ def lfs_inames(self, *args):
+ raise NotImplementedError("Basis Mixins should implement local function space inames")
+
+ def implement_basis(self, *args):
+ raise NotImplementedError("Basis Mixins should implement the basis")
+
+ def implement_reference_gradient(self, *args):
+ raise NotImplementedError("Basis Mixins should implement the basis gradient")
+
+ def implement_trialfunction(self, *args):
+ raise NotImplementedError("Basis Mixins should implement trial function evaluation")
+
+ def implement_trialfunction_gradient(self, *args):
+ raise NotImplementedError("Basis Mixins should implement trial function gradient evaluation")
+
+ def implement_apply_function(self, *args):
+ raise NotImplementedError("Basis Mixins should implement linearization point evaluation")
+
+ def implement_apply_function_gradient(self, *args):
+ raise NotImplementedError("Basis Mixins should implement linearization point gradient evaluation")
+
+
+@basis_mixin("generic")
+class GenericBasisMixin(BasisMixinBase):
+ 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 rif 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,
+ )
+
+
@backend(interface="typedef_localbasis")
@class_member(classtag="operator")
def typedef_localbasis(element, name):
@@ -93,62 +277,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):
@@ -156,71 +284,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/geometry.py b/python/dune/codegen/pdelab/geometry.py
index ba59e267..0efd5a27 100644
--- a/python/dune/codegen/pdelab/geometry.py
+++ b/python/dune/codegen/pdelab/geometry.py
@@ -18,13 +18,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 (pymbolic_quadrature_position,
- 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
@@ -38,25 +36,74 @@ 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 unhandled(self, o):
+ 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 = unhandled
- facet_normal = unhandled
- jacobian_determinant = unhandled
- jacobian_inverse = unhandled
- facet_jacobian_determinant = unhandled
- cell_volume = unhandled
- facet_area = unhandled
+ 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 GenericPDELabGeometry(object):
+class GenericPDELabGeometryMixin(GeometryMixinBase):
def spatial_coordinate(self, o):
- return to_global(pymbolic_quadrature_position())
+ 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"
@@ -76,13 +123,14 @@ class GenericPDELabGeometry(object):
def define_jacobian_determinant(self, name):
temporary_variable(name, shape=())
geo = name_geometry()
- pos = get_backend("quad_pos", selector=option_switch("blockstructured"))()
+ pos = self.quadrature_position()
code = "{} = {}.integrationElement({});".format(name,
geo,
str(pos),
)
- return quadrature_preamble(code,
+ return quadrature_preamble(self,
+ code,
assignees=frozenset({name}),
read_variables=frozenset({get_pymbolic_basename(pos)}),
depends_on=frozenset({Writes(get_pymbolic_basename(pos))}),
@@ -104,9 +152,10 @@ class GenericPDELabGeometry(object):
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)
+ pos = self.to_cell(self.quadrature_position())
- return quadrature_preamble("{} = {}.jacobianInverseTransposed({});".format(name,
+ return quadrature_preamble(self,
+ "{} = {}.jacobianInverseTransposed({});".format(name,
geo,
str(pos),
),
@@ -127,8 +176,9 @@ class GenericPDELabGeometry(object):
temporary_variable(name, shape=(world_dimension(),), decl_method=declare_normal)
ig = name_intersection_geometry_wrapper()
- qp = get_backend("quad_pos")()
- return quadrature_preamble("{} = {}.unitOuterNormal({});".format(name, ig, qp),
+ 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))}),
@@ -412,35 +462,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")]
@@ -481,24 +502,3 @@ def define_jacobian_inverse_transposed_temporary(restriction):
name,
)
return _define_jacobian_inverse_transposed_temporary
-
-
-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)
diff --git a/python/dune/codegen/pdelab/localoperator.py b/python/dune/codegen/pdelab/localoperator.py
index 44446de5..252b7f8b 100644
--- a/python/dune/codegen/pdelab/localoperator.py
+++ b/python/dune/codegen/pdelab/localoperator.py
@@ -425,7 +425,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)
@@ -452,8 +452,16 @@ 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").split(",")
- VisitorType = construct_from_mixins(base=UFL2LoopyVisitor, mixins=mixins, name="UFLVisitor")
+ mixins = get_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(",")
+ VisitorType = construct_from_mixins(base=VisitorType, mixins=mixins, mixintype="quadrature", name="UFLVisitor")
+
+ # Mix basis mixins in
+ mixins = get_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 8c3e2b36..8066423b 100644
--- a/python/dune/codegen/pdelab/quadrature.py
+++ b/python/dune/codegen/pdelab/quadrature.py
@@ -10,13 +10,100 @@ from dune.codegen.generation import (backend,
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):
+ 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):
+ 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):
+ 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)
+
+ return prim.Subscript(prim.Variable(name), tuple(prim.Variable(i) for i in self.quadrature_inames()))
+
+ @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/sumfact/__init__.py b/python/dune/codegen/sumfact/__init__.py
index f887ceda..5dc7f3b1 100644
--- a/python/dune/codegen/sumfact/__init__.py
+++ b/python/dune/codegen/sumfact/__init__.py
@@ -60,11 +60,6 @@ class SumFactInterface(PDELabInterface):
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)
diff --git a/python/dune/codegen/sumfact/quadrature.py b/python/dune/codegen/sumfact/quadrature.py
index eedd85da..d5181a28 100644
--- a/python/dune/codegen/sumfact/quadrature.py
+++ b/python/dune/codegen/sumfact/quadrature.py
@@ -8,6 +8,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 +21,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 +41,16 @@ import loopy as lp
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)
+
+
class BaseWeight(FunctionIdentifier):
def __init__(self, accumobj):
self.accumobj = accumobj
diff --git a/python/dune/codegen/ufl/visitor.py b/python/dune/codegen/ufl/visitor.py
index 9c54018e..10afd7bf 100644
--- a/python/dune/codegen/ufl/visitor.py
+++ b/python/dune/codegen/ufl/visitor.py
@@ -131,9 +131,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
@@ -158,14 +158,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:
@@ -421,14 +421,7 @@ class UFL2LoopyVisitor(ModifiedTerminalTracker):
return prim.LogicalNot(self.call(o.ufl_operands[0]))
#
- # Handlers for quadrature related stuff
- #
-
- def quadrature_weight(self, o):
- return self.interface.pymbolic_quadrature_weight()
-
- #
- # NB: Geometry related handlers have been moved into configurable mixin classes!
+ # NB: Geometry and quadrature related handlers have been moved into configurable mixin classes!
#
#
@@ -439,7 +432,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/test/poisson/poisson_dg_quadrilateral.mini b/test/poisson/poisson_dg_quadrilateral.mini
index dde2e495..162e2479 100644
--- a/test/poisson/poisson_dg_quadrilateral.mini
+++ b/test/poisson/poisson_dg_quadrilateral.mini
@@ -11,6 +11,7 @@ extension = vtu
[formcompiler]
compare_l2errorsquared = 7e-7
+geometry = equidistant
[formcompiler.r]
numerical_jacobian = 1, 0 | expand num
--
GitLab