WIP stuff on sumfact geometry mixins

python/dune/codegen/pdelab/geometry.py

@@ -159,7 +159,7 @@ class GenericPDELabGeometry(object):
 
 
 @geometry_mixin("axiparallel")
-class AxiparallelGeometry(GenericPDELabGeometry):
+class AxiparallelGeometryMixin(GenericPDELabGeometryMixin):
     def define_unit_outer_normal(self, name):
         # NB: Strictly speaking, this implementation holds for any non-curved intersection.
         declare_normal(name, None, None)
@@ -173,7 +173,7 @@ class AxiparallelGeometry(GenericPDELabGeometry):
         if i != j:
             return 0
 
-        return GenericPDELabGeometry.jacobian_determinant(self, o)
+        return GenericPDELabGeometryMixin.jacobian_determinant(self, o)
 
     def facet_area(self, o):
         # This is not 100% correct, but in practice structured grid implementations are not
@@ -200,7 +200,7 @@ class AxiparallelGeometry(GenericPDELabGeometry):
 
 
 @geometry_mixin("equidistant")
-class EquidistantGeometry(AxiparallelGeometry):
+class EquidistantGeometryMixin(AxiparallelGeometryMixin):
     @class_member(classtag="operator")
     def define_jacobian_determinant(self, name):
         valuearg(name)

python/dune/codegen/sumfact/__init__.py

@@ -76,27 +76,3 @@ class SumFactInterface(PDELabInterface):
         import dune.codegen.sumfact.geometry
         ret = get_backend(interface="spatial_coordinate", selector=option_switch("diagonal_transformation_matrix"))(self.visitor.do_predicates, self.visitor)
         return ret
-
-    def pymbolic_unit_outer_normal(self):
-        from dune.codegen.sumfact.geometry import pymbolic_unit_outer_normal
-        return pymbolic_unit_outer_normal(self.visitor)
-
-    def pymbolic_unit_inner_normal(self):
-        from dune.codegen.sumfact.geometry import pymbolic_unit_inner_normal
-        return pymbolic_unit_inner_normal(self.visitor)
-
-    def pymbolic_facet_jacobian_determinant(self):
-        from dune.codegen.sumfact.geometry import pymbolic_facet_jacobian_determinant
-        return pymbolic_facet_jacobian_determinant(self.visitor)
-
-    def pymbolic_facet_area(self):
-        from dune.codegen.sumfact.geometry import pymbolic_facet_area
-        return pymbolic_facet_area(self.visitor)
-
-    def pymbolic_jacobian_inverse(self, i, j, restriction):
-        from dune.codegen.sumfact.geometry import pymbolic_jacobian_inverse
-        return pymbolic_jacobian_inverse(i, j, restriction, self.visitor)
-
-    def pymbolic_jacobian_determinant(self):
-        from dune.codegen.sumfact.geometry import pymbolic_jacobian_determinant
-        return pymbolic_jacobian_determinant(self.visitor)

python/dune/codegen/sumfact/geometry.py

@@ -3,6 +3,7 @@
 from dune.codegen.generation import (backend,
                                      class_member,
                                      domain,
+                                     geometry_mixin,
                                      get_counted_variable,
                                      iname,
                                      include_file,
@@ -20,6 +21,9 @@ from dune.codegen.pdelab.geometry import (local_dimension,
                                           world_dimension,
                                           name_cell_geometry,
                                           name_geometry,
+                                          GenericPDELabGeometryMixin,
+                                          AxiparallelGeometryMixin,
+                                          EquidistantGeometryMixin,
                                           )
 from dune.codegen.pdelab.localoperator import (name_ansatz_gfs_constructor_param,
                                                lop_template_ansatz_gfs,
@@ -48,6 +52,62 @@ import pymbolic.primitives as prim
 import numpy as np
 
 
+@geometry_mixin("sumfact_multilinear")
+class SumfactMultiLinearGeometryMixin(GenericPDELabGeometryMixin):
+    def pymbolic_unit_outer_normal(self):
+        return pymbolic_unit_outer_normal(self)
+
+    def pymbolic_jacobian_inverse(self, i, j, restriction):
+        return pymbolic_jacobian_inverse(i, j, restriction, self.visitor)
+
+
+@geometry_mixin("sumfact_axiparallel")
+class SumfactAxiParallelGeometryMixin(AxiparallelGeometryMixin):
+    def facet_normal(self, o):
+        i, = self.indices
+        assert isinstance(i, int)
+
+        # Use facemod_s and facedir_s
+        if i == get_facedir(Restriction.POSITIVE):
+            if get_facemod(Restriction.POSITIVE):
+                return 1
+            else:
+                return -1
+        else:
+            return 0
+
+
+@geometry_mixin("sumfact_equidistant")
+class SumfactAxiParallelGeometryMixin(EquidistantGeometryMixin, SumfactAxiParallelGeometryMixin):
+    def facet_jacobian_determinant(self, o):
+        name = "fdetjac"
+        self.define_facet_jacobian_determinant(name)
+        facedir = get_facedir(Restriction.POSITIVE)
+        globalarg(name, shape=(world_dimension(),))
+        return prim.Subscript(prim.Variable(name), (facedir,))
+
+    @class_member(classtag="operator")
+    def define_facet_jacobian_determinant(self, name):
+        self._define_facet_jacobian_determinant_eval(name)
+        gfst = lop_template_ansatz_gfs()
+        return "std::array<typename {}::Traits::GridView::template Codim<0>::Geometry::ctype, {}> {};".format(gfst, world_dimension(), name)
+
+    @kernel_cached(kernel="operator")
+    def _define_facet_jacobian_determinant_eval(self, name):
+        gfs = name_ansatz_gfs_constructor_param()
+        rft = lop_template_range_field()
+        code = ["{",
+                "  auto e = *({}.gridView().template begin<0>());".format(gfs),
+                "  int dir=0;",
+                "  for(auto is = {0}.gridView().ibegin(e); is != {0}.gridView().iend(e); ++(++is), ++dir)".format(gfs),
+                "    {}[dir] = is->geometry().integrationElement(Dune::FieldVector<{}, {}>());".format(name, rft, world_dimension() - 1),
+                "}",
+                ]
+        instruction(code="\n".join(code),
+                    kernel="operator",
+                    )
+
+
 @iname
 def global_corner_iname(restriction):
     name = get_counted_variable(restricted_name("global_corneriname", restriction))
@@ -409,39 +469,6 @@ def pymbolic_unit_inner_normal(visitor_indices):
         return prim.Variable(name)
 
 
-@kernel_cached(kernel="operator")
-def define_constant_facet_jacobian_determinant_eval(name):
-    gfs = name_ansatz_gfs_constructor_param()
-    rft = lop_template_range_field()
-    code = ["{",
-            "  auto e = *({}.gridView().template begin<0>());".format(gfs),
-            "  int dir=0;",
-            "  for(auto is = {0}.gridView().ibegin(e); is != {0}.gridView().iend(e); ++(++is), ++dir)".format(gfs),
-            "    {}[dir] = is->geometry().integrationElement(Dune::FieldVector<{}, {}>());".format(name, rft, world_dimension() - 1),
-            "}",
-            ]
-    instruction(code="\n".join(code),
-                kernel="operator",
-                )
-
-
-@class_member(classtag="operator")
-def define_constant_facet_jacobian_determinant(name):
-    define_constant_facet_jacobian_determinant_eval(name)
-    gfst = lop_template_ansatz_gfs()
-    return "std::array<typename {}::Traits::GridView::template Codim<0>::Geometry::ctype, {}> {};".format(gfst, world_dimension(), name)
-
-
-def pymbolic_constant_facet_jacobian_determinant():
-    facedir = get_facedir(Restriction.POSITIVE)
-    assert isinstance(facedir, int)
-
-    name = "fdetjac"
-    define_constant_facet_jacobian_determinant(name)
-    globalarg(name, shape=(world_dimension(),))
-    return prim.Subscript(prim.Variable(name), (facedir,))
-
-
 def define_facet_jacobian_determinant(name, visitor):
     """Absolute value of determinant of jacobian of facet geometry mapping
 
@@ -489,19 +516,13 @@ def name_facet_jacobian_determinant(visitor):
 
 
 def pymbolic_facet_jacobian_determinant(visitor):
-    if get_form_option("constant_transformation_matrix"):
-        return pymbolic_constant_facet_jacobian_determinant()
-    else:
-        name = name_facet_jacobian_determinant(visitor)
-        return prim.Variable(name)
+    name = name_facet_jacobian_determinant(visitor)
+    return prim.Variable(name)
 
 
 def pymbolic_facet_area(visitor):
-    if get_form_option("constant_transformation_matrix"):
-        return pymbolic_facet_jacobian_determinant(visitor)
-    else:
-        from dune.codegen.pdelab.geometry import pymbolic_facet_area as nonconstant_facet_area
-        return nonconstant_facet_area()
+    from dune.codegen.pdelab.geometry import pymbolic_facet_area as nonconstant_facet_area
+    return nonconstant_facet_area()
 
 
 def _name_jacobian(i, j, restriction, visitor):
@@ -601,18 +622,8 @@ def name_jacobian_inverse(restriction, visitor):
 
 
 def pymbolic_jacobian_inverse(i, j, restriction, visitor):
-    # Switch between implementations without using the backend mechanism. In
-    # the case of constant_transformation_matrix we want to use the pdelab
-    # branch, otherwise we want to go into the sumfact implementation instead
-    # of pdelab.
-    if get_form_option("constant_transformation_matrix"):
-        from dune.codegen.pdelab.geometry import name_constant_jacobian_inverse_transposed
-        name = name_constant_jacobian_inverse_transposed(restriction)
-        # Return jacobian inverse -> (j, i) instead of (i, j)
-        return prim.Subscript(prim.Variable(name), (j, i))
-    else:
-        name = name_jacobian_inverse(restriction, visitor)
-        return prim.Subscript(prim.Variable(name), (i, j))
+    name = name_jacobian_inverse(restriction, visitor)
+    return prim.Subscript(prim.Variable(name), (i, j))
 
 
 def name_jacobian_determinant(visitor):
@@ -621,15 +632,11 @@ def name_jacobian_determinant(visitor):
 
 
 def pymbolic_jacobian_determinant(visitor):
-    if get_form_option("constant_transformation_matrix"):
-        from dune.codegen.pdelab.geometry import pymbolic_jacobian_determinant
-        return pymbolic_jacobian_determinant()
+    # The calculation of the jacobian determinant happens in this function
+    if visitor.measure == 'cell':
+        restriction = 0
     else:
-        # The calculation of the jacobian determinant happens in this function
-        if visitor.measure == 'cell':
-            restriction = 0
-        else:
-            restriction = 1
-        name_jacobian_inverse(restriction, visitor)
+        restriction = 1
+    name_jacobian_inverse(restriction, visitor)
 
-        return prim.Variable(name_jacobian_determinant(visitor))
+    return prim.Variable(name_jacobian_determinant(visitor))

test/sumfact/poisson/poisson_2d.mini

@@ -21,6 +21,7 @@ numerical_jacobian = 1, 0 | expand num
 sumfact = 1
 vectorization_strategy = explicit, none | expand grad
 quadrature_order = 2, 4
+geometry = sumfact_axiparallel
 
 [formcompiler.ufl_variants]
 degree = 1, 2 | expand deg