diff --git a/python/dune/codegen/blockstructured/geometry.py b/python/dune/codegen/blockstructured/geometry.py
index 4f1a71d3e3674ea768ef63842ca5aa4def21db61..4bc0e6131ccfd192f14cc6e2ab15f37b99a7841b 100644
--- a/python/dune/codegen/blockstructured/geometry.py
+++ b/python/dune/codegen/blockstructured/geometry.py
@@ -9,8 +9,6 @@ from dune.codegen.options import (get_form_option,
option_switch, get_option)
from dune.codegen.pdelab.geometry import (world_dimension,
local_dimension,
- name_facet_jacobian_determinant,
- name_element_corners,
component_iname)
from dune.codegen.blockstructured.tools import sub_element_inames
import pymbolic.primitives as prim
@@ -274,6 +272,7 @@ def pymbolic_jacobian_inverse(i, j, restriction):
# scale determinant according to the order of the blockstructure
def pymbolic_facet_jacobian_determinant():
+ from dune.codegen.pdelab.geometry import name_facet_jacobian_determinant
return prim.Quotient(prim.Variable(name_facet_jacobian_determinant()),
prim.Power(get_form_option("number_of_blocks"), local_dimension()))
@@ -391,3 +390,22 @@ def to_global(local):
raise NotImplementedError
return prim.Variable(name)
+
+
+def define_element_corners(name):
+ from dune.codegen.pdelab.driver import get_form
+ n_corners = get_form().ufl_cell().num_vertices()
+ temporary_variable(name, shape_impl=('fv', 'fv'), shape=(n_corners, local_dimension()))
+
+ iname = "i_corner"
+ domain(iname, n_corners)
+ from dune.codegen.generation import instruction
+ instruction(code="{0}[{1}] = {2}.corner({3});".format(name, iname, name_geometry(), iname),
+ assignees=frozenset({name}),
+ within_inames=frozenset({iname}), within_inames_is_final=True)
+
+
+def name_element_corners():
+ name = "corners"
+ define_element_corners(name)
+ return name
diff --git a/python/dune/codegen/generation/__init__.py b/python/dune/codegen/generation/__init__.py
index 3103bac657648779921eb66c7170f5a55059e209..cf729760ff5751967e28c0547645b2616512d0f5 100644
--- a/python/dune/codegen/generation/__init__.py
+++ b/python/dune/codegen/generation/__init__.py
@@ -58,3 +58,7 @@ from dune.codegen.generation.context import (cache_restoring,
global_context,
get_global_context_value,
)
+
+from dune.codegen.generation.mixins import (construct_from_mixins,
+ geometry_mixin,
+ )
diff --git a/python/dune/codegen/generation/cache.py b/python/dune/codegen/generation/cache.py
index c1622cb1528f78fe4425334343484d963e52d333..c727fe0516e02e1197dc37bd6168103668492a78 100644
--- a/python/dune/codegen/generation/cache.py
+++ b/python/dune/codegen/generation/cache.py
@@ -8,6 +8,7 @@ from dune.codegen.generation.context import (get_global_context_value,
)
from dune.codegen.generation.counter import get_counter
from dune.codegen.options import get_option
+from pytools import ImmutableRecord
# Store a global list of generator functions
_generators = {}
diff --git a/python/dune/codegen/generation/mixins.py b/python/dune/codegen/generation/mixins.py
new file mode 100644
index 0000000000000000000000000000000000000000..a8a7e72b91c11d12eb24a71e8a3b50fcd8a0b1b4
--- /dev/null
+++ b/python/dune/codegen/generation/mixins.py
@@ -0,0 +1,22 @@
+""" The infrastructure for registered mixin classes """
+
+_mixin_registry = {}
+
+
+def mixin_base(name, mixintype):
+ _mixin_registry.setdefault(mixintype, {})
+
+ def _dec(cls):
+ _mixin_registry[mixintype][name] = cls
+ return cls
+
+ 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,), {})
diff --git a/python/dune/codegen/options.py b/python/dune/codegen/options.py
index 40f2d1a520ca5506604d8f5ec98709d692b75fc3..b75d7b41dcc75e964fdd1a90a37c899a6ca26f6d 100644
--- a/python/dune/codegen/options.py
+++ b/python/dune/codegen/options.py
@@ -55,6 +55,7 @@ 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")
# 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 affa0725c89539fa88e470ac143f5d76fdb8a18b..42915eac2c5dbac4f3af71d18c3000357986e231 100644
--- a/python/dune/codegen/pdelab/__init__.py
+++ b/python/dune/codegen/pdelab/__init__.py
@@ -14,18 +14,7 @@ from dune.codegen.pdelab.basis import (pymbolic_basis,
pymbolic_reference_gradient,
)
from dune.codegen.pdelab.function import pymbolic_gridfunction
-from dune.codegen.pdelab.geometry import (component_iname,
- pymbolic_cell_volume,
- pymbolic_facet_area,
- pymbolic_facet_jacobian_determinant,
- pymbolic_jacobian_determinant,
- pymbolic_jacobian_inverse,
- pymbolic_unit_inner_normal,
- pymbolic_unit_outer_normal,
- to_global,
- )
-from dune.codegen.pdelab.index import (name_index,
- )
+from dune.codegen.pdelab.index import name_index
from dune.codegen.pdelab.quadrature import (pymbolic_quadrature_weight,
pymbolic_quadrature_position,
quadrature_inames,
@@ -118,34 +107,6 @@ class PDELabInterface(object):
def pymbolic_matrix_inverse(self, o):
return pymbolic_matrix_inverse(o, self.visitor)
- #
- # Geometry related generator functions
- #
-
- def pymbolic_spatial_coordinate(self):
- return to_global(pymbolic_quadrature_position())
-
- def pymbolic_facet_jacobian_determinant(self):
- return pymbolic_facet_jacobian_determinant()
-
- def pymbolic_jacobian_determinant(self):
- return pymbolic_jacobian_determinant()
-
- def pymbolic_jacobian_inverse(self, i, j, restriction):
- return pymbolic_jacobian_inverse(i, j, restriction)
-
- def pymbolic_unit_inner_normal(self):
- return pymbolic_unit_inner_normal()
-
- def pymbolic_unit_outer_normal(self):
- return pymbolic_unit_outer_normal()
-
- def pymbolic_cell_volume(self, restriction):
- return pymbolic_cell_volume(restriction)
-
- def pymbolic_facet_area(self):
- return pymbolic_facet_area()
-
#
# Quadrature related generator functions
#
diff --git a/python/dune/codegen/pdelab/basis.py b/python/dune/codegen/pdelab/basis.py
index c65e71ebfc9e3a38023033e27dc65d445734edbf..6b6afb3b95dd65fcfa7ef5bb27167da674a090a9 100644
--- a/python/dune/codegen/pdelab/basis.py
+++ b/python/dune/codegen/pdelab/basis.py
@@ -22,7 +22,6 @@ from dune.codegen.pdelab.spaces import (lfs_iname,
)
from dune.codegen.pdelab.geometry import (component_iname,
world_dimension,
- name_jacobian_inverse_transposed,
to_cell_coordinates,
name_cell,
)
diff --git a/python/dune/codegen/pdelab/geometry.py b/python/dune/codegen/pdelab/geometry.py
index 7d1e722ad8d92143d44747cbc2622b7eec0431b2..07b2f0b88af3751c9d9bea94d893f3aa24036429 100644
--- a/python/dune/codegen/pdelab/geometry.py
+++ b/python/dune/codegen/pdelab/geometry.py
@@ -3,6 +3,7 @@ from dune.codegen.pdelab.restriction import restricted_name
from dune.codegen.generation import (backend,
class_member,
domain,
+ geometry_mixin,
get_backend,
get_global_context_value,
globalarg,
@@ -17,7 +18,9 @@ from dune.codegen.options import (get_form_option,
option_switch,
)
from dune.codegen.loopy.target import dtype_floatingpoint, type_floatingpoint
-from dune.codegen.pdelab.quadrature import quadrature_preamble
+from dune.codegen.pdelab.quadrature import (pymbolic_quadrature_position,
+ quadrature_preamble,
+ )
from dune.codegen.tools import get_pymbolic_basename
from ufl.algorithms import MultiFunction
from pymbolic.primitives import Variable
@@ -29,6 +32,214 @@ import numpy as np
import pymbolic.primitives as prim
+@geometry_mixin("base")
+class GeometryMixinBase(object):
+ """
+ This mixin will be in by default by the mixin magic,
+ so it should only define an interface and throw exceptions.
+ """
+ def 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
+
+
+@geometry_mixin("generic")
+class GenericPDELabGeometry(object):
+ def spatial_coordinate(self, o):
+ return to_global(pymbolic_quadrature_position())
+
+ def facet_jacobian_determinant(self, o):
+ name = "fdetjac"
+ self.define_facet_jacobian_determinant(name)
+ return prim.Variable(name)
+
+ def define_facet_jacobian_determinant(self, name):
+ # NB: This can reuse the same implementation as the cell jacobian
+ # determinant, because the obtain geometry will be a face geometry
+ return self.define_jacobian_determinant(name)
+
+ def jacobian_determinant(self, o):
+ name = 'detjac'
+ self.define_jacobian_determinant(name)
+ return prim.Variable(name)
+
+ def define_jacobian_determinant(self, name):
+ temporary_variable(name, shape=())
+ geo = name_geometry()
+ pos = get_backend("quad_pos", selector=option_switch("blockstructured"))()
+ code = "{} = {}.integrationElement({});".format(name,
+ geo,
+ str(pos),
+ )
+
+ return quadrature_preamble(code,
+ assignees=frozenset({name}),
+ read_variables=frozenset({get_pymbolic_basename(pos)}),
+ depends_on=frozenset({Writes(get_pymbolic_basename(pos))}),
+ )
+
+ def jacobian_inverse(self, o):
+ restriction = enforce_boundary_restriction(self)
+
+ assert(len(self.indices) == 2)
+ i, j = self.indices
+ self.indices = None
+
+ name = restricted_name("jit", restriction)
+ self.define_jacobian_inverse_transposed(name, restriction)
+
+ return prim.Subscript(prim.Variable(name), (j,i))
+
+ def define_jacobian_inverse_transposed(self, name, restriction):
+ dim = world_dimension()
+ temporary_variable(name, decl_method=define_jacobian_inverse_transposed_temporary(restriction), shape=(dim, dim))
+ geo = name_cell_geometry(restriction)
+ pos = get_backend("qp_in_cell", selector=option_switch("blockstructured"))(restriction)
+
+ return quadrature_preamble("{} = {}.jacobianInverseTransposed({});".format(name,
+ geo,
+ str(pos),
+ ),
+ assignees=frozenset({name}),
+ read_variables=frozenset({get_pymbolic_basename(pos)}),
+ depends_on=frozenset({Writes(get_pymbolic_basename(pos))}),
+ )
+
+ def facet_normal(self, o):
+ if self.restriction == Restriction.NEGATIVE:
+ raise NotImplementedError("Inner Normals not implemented!")
+
+ name = "unit_outer_normal"
+ self.define_unit_outer_normal(name)
+ return prim.Variable(name)
+
+ def define_unit_outer_normal(self, name):
+ temporary_variable(name, shape=(world_dimension(),), decl_method=declare_normal)
+
+ ig = name_intersection_geometry_wrapper()
+ qp = get_backend("quad_pos")()
+ return quadrature_preamble("{} = {}.unitOuterNormal({});".format(name, ig, qp),
+ assignees=frozenset({name}),
+ read_variables=frozenset({get_pymbolic_basename(qp)}),
+ depends_on=frozenset({Writes(get_pymbolic_basename(qp))}),
+ )
+
+ def cell_volume(self, o):
+ restriction = enforce_boundary_restriction(self)
+ name = restricted_name("volume", restriction)
+ self.define_cell_volume(name, restriction)
+ return prim.Variable(name)
+
+ @preamble
+ def define_cell_volume(self, name, restriction):
+ geo = name_cell_geometry(restriction)
+ valuearg(name)
+ return "auto {} = {}.volume();".format(name, geo)
+
+ def facet_area(self, o):
+ name = "area"
+ self.define_facet_area(name)
+ return prim.Variable(name)
+
+ @preamble
+ def define_facet_area(self, name):
+ geo = name_intersection_geometry()
+ valuearg(name)
+ return "auto {} = {}.volume();".format(name, geo)
+
+
+@geometry_mixin("axiparallel")
+class AxiparallelGeometry(GenericPDELabGeometry):
+ def define_unit_outer_normal(self, name):
+ # NB: Strictly speaking, this implementation holds for any non-curved intersection.
+ declare_normal(name, None, None)
+ globalarg(name, shape=(world_dimension(),))
+
+ def jacobian_inverse(self, o):
+ i, j = self.indices
+ self.indices = None
+ assert isinstance(i, int) and isinstance(j, int)
+
+ if i != j:
+ return 0
+
+ return GenericPDELabGeometry.jacobian_determinant(self, o)
+
+ def facet_area(self, o):
+ # This is not 100% correct, but in practice structured grid implementations are not
+ # embedded into higher dimensional world space.
+ return self.facet_jacobian_determinant(o)
+
+ def cell_volume(self, o):
+ # This is not 100% correct, but in practice structured grid implementations are not
+ # embedded into higher dimensional world space.
+ return self.jacobian_determinant(o)
+
+ @preamble
+ def define_facet_jacobian_determinant(self, name):
+ # NB: In the equidistant case, this might be optimized to store *d* values on the operator level
+ # We don't do that right now for laziness.
+ geo = name_geometry()
+ pos = name_localcenter()
+ valuearg(name)
+
+ return "auto {} = {}.integrationElement({});".format(name,
+ geo,
+ pos,
+ )
+
+
+@geometry_mixin("equidistant")
+class EquidistantGeometry(AxiparallelGeometry):
+ @class_member(classtag="operator")
+ def define_jacobian_determinant(self, name):
+ valuearg(name)
+ self._define_jacobian_determinant_eval(name)
+ from dune.codegen.pdelab.localoperator import lop_template_ansatz_gfs
+ gfst = lop_template_ansatz_gfs()
+ return "typename {}::Traits::GridView::template Codim<0>::Geometry::ctype {};".format(gfst, name)
+
+ @preamble(kernel="operator")
+ def _define_jacobian_determinant_eval(name):
+ from dune.codegen.pdelab.localoperator import name_ansatz_gfs_constructor_param, lop_template_range_field
+ gfs = name_ansatz_gfs_constructor_param()
+ rft = lop_template_range_field()
+ return "{} = {}.gridView().template begin<0>()->geometry().integrationElement(Dune::FieldVector<{}, {}>());".format(name, gfs, rft, world_dimension())
+
+ @class_member(classtag="operator")
+ def define_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()
+ return "typename {}::Traits::GridView::template Codim<0>::Geometry::JacobianInverseTransposed {};".format(gfst, name)
+
+ @preamble(kernel="operator")
+ def _define_jacobian_inverse_transposed_eval(name):
+ from dune.codegen.pdelab.localoperator import name_ansatz_gfs_constructor_param, lop_template_range_field
+ gfs = name_ansatz_gfs_constructor_param()
+ rft = lop_template_range_field()
+ return "{} = {}.gridView().template begin<0>()->geometry().jacobianInverseTransposed(Dune::FieldVector<{}, {}>());".format(name, gfs, rft, world_dimension())
+
+
+def enforce_boundary_restriction(visitor):
+ if visitor.measure == 'exterior_facet' and visitor.restriction == Restriction.NONE:
+ return Restriction.POSITIVE
+ else:
+ return visitor.restriction
+
+
@preamble
def define_reference_element(name):
geo = name_geometry()
@@ -250,48 +461,12 @@ def local_dimension():
return intersection_dimension()
-def evaluate_unit_outer_normal(name):
- ig = name_intersection_geometry_wrapper()
- qp = get_backend("quad_pos")()
- return quadrature_preamble("{} = {}.unitOuterNormal({});".format(name, ig, qp),
- assignees=frozenset({name}),
- read_variables=frozenset({get_pymbolic_basename(qp)}),
- depends_on=frozenset({Writes(get_pymbolic_basename(qp))}),
- )
-
-
@preamble
def declare_normal(name, kernel, decl_info):
ig = name_intersection_geometry_wrapper()
return "auto {} = {}.centerUnitOuterNormal();".format(name, ig)
-def pymbolic_unit_outer_normal():
- name = "unit_outer_normal"
- if not get_form_option("diagonal_transformation_matrix"):
- temporary_variable(name, shape=(world_dimension(),), decl_method=declare_normal)
- evaluate_unit_outer_normal(name)
- else:
- declare_normal(name, None, None)
- globalarg(name, shape=(world_dimension(),))
- return prim.Variable(name)
-
-
-def evaluate_unit_inner_normal(name):
- outer = name_unit_outer_normal()
- return quadrature_preamble("auto {} = -1. * {};".format(name, outer),
- assignees=frozenset({name}),
- read_variables=frozenset({outer}),
- )
-
-
-def pymbolic_unit_inner_normal():
- name = "inner_normal"
- temporary_variable(name, shape=(world_dimension(),), decl_method=declare_normal)
- evaluate_unit_inner_normal(name)
- return prim.Variable(name)
-
-
def type_jacobian_inverse_transposed(restriction):
geo = type_cell_geometry(restriction)
return "typename {}::JacobianInverseTransposed".format(geo)
@@ -308,144 +483,6 @@ def define_jacobian_inverse_transposed_temporary(restriction):
return _define_jacobian_inverse_transposed_temporary
-@preamble(kernel="operator")
-def define_constant_jacobian_inverse_transposed_eval(name):
- from dune.codegen.pdelab.localoperator import name_ansatz_gfs_constructor_param, lop_template_range_field
- gfs = name_ansatz_gfs_constructor_param()
- rft = lop_template_range_field()
-
- return "{} = {}.gridView().template begin<0>()->geometry().jacobianInverseTransposed(Dune::FieldVector<{}, {}>());".format(name, gfs, rft, world_dimension())
-
-
-@class_member(classtag="operator")
-def define_constant_jacobian_inverse_transposed(name):
- define_constant_jacobian_inverse_transposed_eval(name)
- from dune.codegen.pdelab.localoperator import lop_template_ansatz_gfs
- gfst = lop_template_ansatz_gfs()
- return "typename {}::Traits::GridView::template Codim<0>::Geometry::JacobianInverseTransposed {};".format(gfst, name)
-
-
-@backend(interface="name_jit", name="constant_transformation_matrix")
-def name_constant_jacobian_inverse_transposed(restriction):
- name = "jit"
- dim = world_dimension()
- globalarg(name, shape=(dim, dim), managed=False)
- define_constant_jacobian_inverse_transposed(name)
- return name
-
-
-def define_jacobian_inverse_transposed(name, restriction):
- dim = world_dimension()
- temporary_variable(name, decl_method=define_jacobian_inverse_transposed_temporary(restriction), shape=(dim, dim))
- geo = name_cell_geometry(restriction)
- pos = get_backend("qp_in_cell", selector=option_switch("blockstructured"))(restriction)
-
- return quadrature_preamble("{} = {}.jacobianInverseTransposed({});".format(name,
- geo,
- str(pos),
- ),
- assignees=frozenset({name}),
- read_variables=frozenset({get_pymbolic_basename(pos)}),
- depends_on=frozenset({Writes(get_pymbolic_basename(pos))}),
- )
-
-
-@backend(interface="name_jit", name="default")
-def name_jacobian_inverse_transposed(restriction):
- name = restricted_name("jit", restriction)
- define_jacobian_inverse_transposed(name, restriction)
- return name
-
-
-def pymbolic_jacobian_inverse(i, j, restriction):
- # Calculate jacobian_inverse_transposed
- name = get_backend(interface="name_jit", selector=option_switch("constant_transformation_matrix"))(restriction)
-
- # Return jacobian inverse -> (j, i) instead of (i, j)
- return prim.Subscript(prim.Variable(name), (j, i))
-
-
-@preamble(kernel="operator")
-def define_constant_jacobian_determinant_eval(name):
- from dune.codegen.pdelab.localoperator import name_ansatz_gfs_constructor_param, lop_template_range_field
- gfs = name_ansatz_gfs_constructor_param()
- rft = lop_template_range_field()
-
- return "{} = {}.gridView().template begin<0>()->geometry().integrationElement(Dune::FieldVector<{}, {}>());".format(name, gfs, rft, world_dimension())
-
-
-@class_member(classtag="operator")
-def _define_constant_jacobian_determinant(name):
- define_constant_jacobian_determinant_eval(name)
- from dune.codegen.pdelab.localoperator import lop_template_ansatz_gfs
- gfst = lop_template_ansatz_gfs()
- return "typename {}::Traits::GridView::template Codim<0>::Geometry::ctype {};".format(gfst, name)
-
-
-@backend(interface="detjac", name="constant_transformation_matrix")
-def define_constant_jacobian_determinant(name):
- valuearg(name)
- _define_constant_jacobian_determinant(name)
-
-
-@backend(interface="detjac", name="default")
-def define_jacobian_determinant(name):
- temporary_variable(name, shape=())
- geo = name_geometry()
- pos = get_backend("quad_pos", selector=option_switch("blockstructured"))()
- code = "{} = {}.integrationElement({});".format(name,
- geo,
- str(pos),
- )
-
- return quadrature_preamble(code,
- assignees=frozenset({name}),
- read_variables=frozenset({get_pymbolic_basename(pos)}),
- depends_on=frozenset({Writes(get_pymbolic_basename(pos))}),
- )
-
-
-@backend(interface="fdetjac", name="constant_transformation_matrix")
-@preamble
-def define_facet_jacobian_determinant(name):
- # NB: This might be optimized to store *d* values on the operator level
- # We don't do that right now for laziness.
- geo = name_geometry()
- pos = name_localcenter()
-
- valuearg(name)
-
- return "auto {} = {}.integrationElement({});".format(name,
- geo,
- pos,
- )
-
-
-@backend(interface="fdetjac", name="default")
-def define_facet_jacobian_determinant(name):
- return define_jacobian_determinant(name)
-
-
-def name_jacobian_determinant():
- name = 'detjac'
- get_backend(interface="detjac", selector=option_switch("constant_transformation_matrix"))(name)
- return name
-
-
-def pymbolic_jacobian_determinant():
- return prim.Variable(name_jacobian_determinant())
-
-
-def name_facet_jacobian_determinant():
- name = 'fdetjac'
- get_backend(interface="fdetjac", selector=option_switch("constant_transformation_matrix"))(name)
- return name
-
-
-def pymbolic_facet_jacobian_determinant():
- return prim.Variable(name_facet_jacobian_determinant())
-
-
def apply_to_global_transformation(name, local):
temporary_variable(name, shape=(world_dimension(),), shape_impl=("fv",))
geo = name_geometry()
@@ -465,54 +502,3 @@ def to_global(local):
name = get_pymbolic_basename(local) + "_global"
apply_to_global_transformation(name, local)
return prim.Variable(name)
-
-
-@preamble
-def define_cell_volume(name, restriction):
- geo = name_cell_geometry(restriction)
- valuearg(name)
- return "auto {} = {}.volume();".format(name, geo)
-
-
-def pymbolic_cell_volume(restriction):
- if get_form_option("constant_transformation_matrix"):
- return pymbolic_jacobian_determinant()
- else:
- name = restricted_name("volume", restriction)
- define_cell_volume(name, restriction)
- return prim.Variable(name)
-
-
-@preamble
-def define_facet_area(name):
- geo = name_intersection_geometry()
- valuearg(name)
- return "auto {} = {}.volume();".format(name, geo)
-
-
-def pymbolic_facet_area():
- if get_form_option("constant_transformation_matrix"):
- return pymbolic_facet_jacobian_determinant()
- else:
- name = "area"
- define_facet_area(name)
- return prim.Variable(name)
-
-
-def define_element_corners(name):
- from dune.codegen.pdelab.driver import get_form
- n_corners = get_form().ufl_cell().num_vertices()
- temporary_variable(name, shape_impl=('fv', 'fv'), shape=(n_corners, local_dimension()))
-
- iname = "i_corner"
- domain(iname, n_corners)
- from dune.codegen.generation import instruction
- instruction(code="{0}[{1}] = {2}.corner({3});".format(name, iname, name_geometry(), iname),
- assignees=frozenset({name}),
- within_inames=frozenset({iname}), within_inames_is_final=True)
-
-
-def name_element_corners():
- name = "corners"
- define_element_corners(name)
- return name
diff --git a/python/dune/codegen/pdelab/localoperator.py b/python/dune/codegen/pdelab/localoperator.py
index 12793f96977e46ccb252285454f11344810892bc..2845b6893f8dd58485320b7455d236257933657d 100644
--- a/python/dune/codegen/pdelab/localoperator.py
+++ b/python/dune/codegen/pdelab/localoperator.py
@@ -13,6 +13,7 @@ from dune.codegen.generation import (backend,
base_class,
class_basename,
class_member,
+ construct_from_mixins,
constructor_parameter,
domain,
dump_accumulate_timer,
@@ -449,8 +450,12 @@ def get_visitor(measure, subdomain_id):
from dune.codegen.pdelab import PDELabInterface
interface = PDELabInterface()
+ # Construct the visitor taking into account geometry mixins
from dune.codegen.ufl.visitor import UFL2LoopyVisitor
- return UFL2LoopyVisitor(interface, measure, subdomain_id)
+ mixins = get_option("geometry").split(",")
+ VisitorType = construct_from_mixins(base=UFL2LoopyVisitor, mixins=mixins, name="UFLVisitor")
+
+ return VisitorType(interface, measure, subdomain_id)
def visit_integral(integral):
@@ -759,14 +764,6 @@ def local_operator_default_settings(operator, form):
for c in sorted(filter(lambda c: c.count() > 2, form.coefficients()), key=lambda c: c.count()):
name_gridfunction_constructor_argument(c)
- # Set some options!
- from dune.codegen.pdelab.driver import isQuadrilateral
- if isQuadrilateral(form.arguments()[0].ufl_element().cell()) and not get_option("grid_unstructured"):
- from dune.codegen.options import set_form_option
- # For Yasp Grids the jacobian of the transformation is diagonal and constant on each cell
- set_form_option('diagonal_transformation_matrix', True)
- set_form_option('constant_transformation_matrix', True)
-
# Add right base classes for stationary/instationary operators
base_class('Dune::PDELab::LocalOperatorDefaultFlags', classtag="operator")
from dune.codegen.pdelab.driver import is_stationary
diff --git a/python/dune/codegen/ufl/visitor.py b/python/dune/codegen/ufl/visitor.py
index 996024cea6ec966cbb8f5463949c7bb84b2560b2..9c54018eae15debabb462e49804acb20f0bb619f 100644
--- a/python/dune/codegen/ufl/visitor.py
+++ b/python/dune/codegen/ufl/visitor.py
@@ -421,67 +421,15 @@ class UFL2LoopyVisitor(ModifiedTerminalTracker):
return prim.LogicalNot(self.call(o.ufl_operands[0]))
#
- # Handlers for geometric quantities
+ # Handlers for quadrature related stuff
#
- def spatial_coordinate(self, o):
- return self.interface.pymbolic_spatial_coordinate()
-
- def facet_normal(self, o):
- # The normal must be restricted to be well-defined
- assert self.restriction is not Restriction.NONE
-
- # Note: In UFL the jump is defined as: jump(v) = v('+') -
- # v('-'). The corresponding outer unit normal is
- # n=FacetNormal(cell)('+'). In order to be consisten with UFL
- # we need to create the outer unit normal if the restriction
- # is positive.
- if self.restriction == Restriction.POSITIVE:
- return self.interface.pymbolic_unit_outer_normal()
- if self.restriction == Restriction.NEGATIVE:
- # It is highly unnatural to have this generator function,
- # but I do run into subtle trouble with return -1*outer
- # as the indexing into the normal happens only later.
- # Not investing more time into this cornercase right now.
- return self.interface.pymbolic_unit_inner_normal()
-
def quadrature_weight(self, o):
return self.interface.pymbolic_quadrature_weight()
- def jacobian_determinant(self, o):
- return self.interface.pymbolic_jacobian_determinant()
-
- def jacobian_inverse(self, o):
- restriction = self.restriction
- if self.measure == 'exterior_facet':
- restriction = Restriction.POSITIVE
-
- assert(len(self.indices) == 2)
- i, j = self.indices
- self.indices = None
-
- # Implement diagonal jacobians for unrolled matrices!
- if get_form_option("diagonal_transformation_matrix"):
- if isinstance(i, int) and isinstance(j, int) and i != j:
- return 0
-
- return self.interface.pymbolic_jacobian_inverse(i, j, restriction)
-
- def jacobian(self, o):
- raise NotImplementedError("How did you get Jacobian into your form? We only support JacobianInverse right now. Report!")
-
- def facet_jacobian_determinant(self, o):
- return self.interface.pymbolic_facet_jacobian_determinant()
-
- def cell_volume(self, o):
- restriction = self.restriction
- if self.measure == 'exterior_facet':
- restriction = Restriction.POSITIVE
-
- return self.interface.pymbolic_cell_volume(restriction)
-
- def facet_area(self, o):
- return self.interface.pymbolic_facet_area()
+ #
+ # NB: Geometry related handlers have been moved into configurable mixin classes!
+ #
#
# Equality/Hashability of the visitor