From dfb82f862dd9ee0169dc51ca48efb19a488cdfb0 Mon Sep 17 00:00:00 2001
From: Dominic Kempf <dominic.kempf@iwr.uni-heidelberg.de>
Date: Tue, 4 Dec 2018 10:40:49 +0100
Subject: [PATCH] WIP stuff on sumfact geometry mixins
---
python/dune/codegen/pdelab/geometry.py | 6 +-
python/dune/codegen/sumfact/__init__.py | 24 -----
python/dune/codegen/sumfact/geometry.py | 137 +++++++++++++-----------
test/sumfact/poisson/poisson_2d.mini | 1 +
4 files changed, 76 insertions(+), 92 deletions(-)
diff --git a/python/dune/codegen/pdelab/geometry.py b/python/dune/codegen/pdelab/geometry.py
index 07b2f0b8..ba59e267 100644
--- a/python/dune/codegen/pdelab/geometry.py
+++ b/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)
diff --git a/python/dune/codegen/sumfact/__init__.py b/python/dune/codegen/sumfact/__init__.py
index 22a60496..f887ceda 100644
--- a/python/dune/codegen/sumfact/__init__.py
+++ b/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)
diff --git a/python/dune/codegen/sumfact/geometry.py b/python/dune/codegen/sumfact/geometry.py
index 3d7df550..f7f1fcdc 100644
--- a/python/dune/codegen/sumfact/geometry.py
+++ b/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))
diff --git a/test/sumfact/poisson/poisson_2d.mini b/test/sumfact/poisson/poisson_2d.mini
index d9ce1773..cb10c5bf 100644
--- a/test/sumfact/poisson/poisson_2d.mini
+++ b/test/sumfact/poisson/poisson_2d.mini
@@ -21,6 +21,7 @@ numerical_jacobian = 1, 0 | expand num
sumfact = 1
vectorization_strategy = explicit, none | expand grad
quadrature_order = 2, 4
+geometry = sumfact_axiparallel
[formcompiler.ufl_variants]
degree = 1, 2 | expand deg
--
GitLab