diff --git a/python/dune/perftool/pdelab/geometry.py b/python/dune/perftool/pdelab/geometry.py
index e3026745c6f20aaabc3254f21839f23e4cd93bf9..39a46f211fb85e48384943ee46fd13880514b518 100644
--- a/python/dune/perftool/pdelab/geometry.py
+++ b/python/dune/perftool/pdelab/geometry.py
@@ -253,11 +253,11 @@ def local_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))}),
- )
+ return get_backend("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
diff --git a/python/dune/perftool/sumfact/__init__.py b/python/dune/perftool/sumfact/__init__.py
index aee933137a1508c7e14984d04debcd3f085ee003..d5eed5544f52c76be02cce4e59455d5a53770ba1 100644
--- a/python/dune/perftool/sumfact/__init__.py
+++ b/python/dune/perftool/sumfact/__init__.py
@@ -78,8 +78,10 @@ class SumFactInterface(PDELabInterface):
return ret
def pymbolic_unit_outer_normal(self):
- from dune.perftool.sumfact.geometry import pymbolic_unit_outer_normal
- ret, indices = pymbolic_unit_outer_normal(self.visitor.indices)
+ from dune.perftool.generation import global_context
+ with global_context(visitor=self.visitor):
+ from dune.perftool.sumfact.geometry import pymbolic_unit_outer_normal
+ ret, indices = pymbolic_unit_outer_normal(self.visitor.indices)
self.visitor.indices = indices
return ret
@@ -90,9 +92,23 @@ class SumFactInterface(PDELabInterface):
return ret
def pymbolic_facet_jacobian_determinant(self):
- from dune.perftool.sumfact.geometry import pymbolic_facet_jacobian_determinant
- return pymbolic_facet_jacobian_determinant()
+ from dune.perftool.generation import global_context
+ with global_context(visitor=self.visitor):
+ from dune.perftool.sumfact.geometry import pymbolic_facet_jacobian_determinant
+ return pymbolic_facet_jacobian_determinant()
def pymbolic_facet_area(self):
from dune.perftool.sumfact.geometry import pymbolic_facet_area
return pymbolic_facet_area()
+
+ def pymbolic_jacobian_inverse_transposed(self, i, j, restriction):
+ from dune.perftool.pdelab.geometry import pymbolic_jacobian_inverse_transposed
+ from dune.perftool.generation import global_context
+ with global_context(visitor=self.visitor):
+ return pymbolic_jacobian_inverse_transposed(i, j, restriction)
+
+ def pymbolic_jacobian_determinant(self):
+ from dune.perftool.pdelab.geometry import pymbolic_jacobian_determinant
+ from dune.perftool.generation import global_context
+ with global_context(visitor=self.visitor):
+ return pymbolic_jacobian_determinant()
diff --git a/python/dune/perftool/sumfact/geometry.py b/python/dune/perftool/sumfact/geometry.py
index 0f0f25e8635fd2e6cd89d6e21d811bc7c2f47f9d..f3e5b9fde10587ffb02f802b6984b3309f9b52b5 100644
--- a/python/dune/perftool/sumfact/geometry.py
+++ b/python/dune/perftool/sumfact/geometry.py
@@ -39,7 +39,17 @@ class GeoCornersInput(SumfactKernelInterfaceBase, ImmutableRecord):
def __init__(self, dir):
ImmutableRecord.__init__(self, dir=dir)
- def realize(self, sf, index, insn_dep):
+ @property
+ def stage(self):
+ return 1
+
+ def __repr__(self):
+ return ImmutableRecord.__repr__(self)
+
+ def realize(self, sf, insn_dep):
+ # TODO: When we use vectorization this should be the offset
+ index = 0
+
from dune.perftool.sumfact.realization import name_buffer_storage
name = "input_{}".format(sf.buffer)
temporary_variable(name,
@@ -64,11 +74,13 @@ class GeoCornersInput(SumfactKernelInterfaceBase, ImmutableRecord):
self.dir,
)
- instruction(code=code,
- within_inames=frozenset({ciname}),
- assignees=(name,),
- tags=frozenset({"sumfact_stage{}".format(sf.stage)}),
- )
+ insn = instruction(code=code,
+ within_inames=frozenset({ciname}),
+ assignees=(name,),
+ tags=frozenset({"sumfact_stage{}".format(sf.stage)}),
+ )
+
+ return insn_dep.union(frozenset({insn}))
@backend(interface="spatial_coordinate", name="default")
diff --git a/python/dune/perftool/sumfact/quadrature.py b/python/dune/perftool/sumfact/quadrature.py
index fdd359578cec007e0bfe0ac9fa52fe7783afa0fc..a48d9887f30c3fea1d58b97d7aacf4c5b7d39dde 100644
--- a/python/dune/perftool/sumfact/quadrature.py
+++ b/python/dune/perftool/sumfact/quadrature.py
@@ -2,13 +2,15 @@ from dune.perftool.generation import (backend,
domain,
function_mangler,
get_global_context_value,
+ globalarg,
iname,
instruction,
kernel_cached,
loopy_class_member,
+ preamble,
temporary_variable,
)
-
+from dune.perftool.sumfact.switch import get_facedir
from dune.perftool.sumfact.tabulation import (quadrature_points_per_direction,
local_quadrature_points_per_direction,
name_oned_quadrature_points,
@@ -78,7 +80,6 @@ def pymbolic_base_weight():
return 1.0
-@backend(interface="quad_inames", name="sumfact")
@kernel_cached
def quadrature_inames(element):
if element is None:
@@ -179,9 +180,10 @@ def quadrature_weight(visitor):
return prim.Subscript(lp.symbolic.TaggedVariable(name, "operator_precomputed"), tuple(prim.Variable(i) for i in quadrature_inames(element)))
-def define_quadrature_position(name, index):
+def define_quadrature_position(name, local_index):
qps_per_dir = quadrature_points_per_direction()
- qp_bound = qps_per_dir[index]
+ local_qps_per_dir = local_quadrature_points_per_direction()
+ qp_bound = local_qps_per_dir[index]
expression = Subscript(Variable(name_oned_quadrature_points(qp_bound)),
(Variable(quadrature_inames()[index]),))
instruction(expression=expression,
@@ -192,20 +194,29 @@ def define_quadrature_position(name, index):
)
-def pymbolic_indexed_quadrature_position(index, visitor):
+def pymbolic_indexed_quadrature_position(local_index, visitor):
+ """Return the value of the indexed local quadrature points
+
+ Arguments:
+ ----------
+
+ local_index: Local index to the quadrature point. Local means that on the
+ face of a 3D Element the index always will be 0 or 1 but never 2.
+ visitor: UFL tree visitor
+ """
# Return the non-precomputed version
if not get_form_option("precompute_quadrature_info"):
name = 'pos'
temporary_variable(name, shape=(local_dimension(),), shape_impl=("fv",))
- define_quadrature_position(name, index)
- return prim.Subscript(prim.Variable(name), (index,))
+ define_quadrature_position(name, local_index)
+ return prim.Subscript(prim.Variable(name), (local_index,))
- assert isinstance(index, int)
+ assert isinstance(local_index, int)
dim = local_dimension()
qps_per_dir = quadrature_points_per_direction()
local_qps_per_dir = local_quadrature_points_per_direction()
local_qps_per_dir_str = '_'.join(map(str, local_qps_per_dir))
- name = "quad_points_dim{}_num{}_dir{}".format(dim, local_qps_per_dir_str, index)
+ name = "quad_points_dim{}_num{}_localdir{}".format(dim, local_qps_per_dir_str, local_index)
loopy_class_member(name,
shape=local_qps_per_dir,
@@ -218,8 +229,8 @@ def pymbolic_indexed_quadrature_position(index, visitor):
# Precompute it in the constructor
inames = constructor_quadrature_inames(dim)
assignee = prim.Subscript(prim.Variable(name), tuple(prim.Variable(i) for i in inames))
- expression = Subscript(Variable(name_oned_quadrature_points(local_qps_per_dir[index])),
- (prim.Variable(inames[index])))
+ expression = Subscript(Variable(name_oned_quadrature_points(local_qps_per_dir[local_index])),
+ (prim.Variable(inames[local_index])))
instruction(assignee=assignee,
expression=expression,
within_inames=frozenset(inames),
@@ -234,10 +245,79 @@ def pymbolic_indexed_quadrature_position(index, visitor):
return prim.Subscript(lp.symbolic.TaggedVariable(name, "operator_precomputed"), tuple(prim.Variable(i) for i in quadrature_inames(element)))
+def _additional_inames(visitor):
+ info = visitor.current_info[1]
+ if info is None:
+ element = None
+ else:
+ element = info.element
+
+ if element is not None:
+ from dune.perftool.pdelab.restriction import Restriction
+ restriction = Restriction.NONE
+ if get_global_context_value("integral_type") == "interior_facet":
+ restriction = Restriction.POSITIVE
+ from dune.perftool.sumfact.basis import lfs_inames
+ lfs_inames = lfs_inames(element, restriction)
+ return lfs_inames
+ else:
+ return ()
+
+
+@backend(interface="quadrature_preamble", name="sumfact")
+def quadrature_preamble(code, **kw):
+ visitor = get_global_context_value('visitor')
+
+ info = visitor.current_info[1]
+ if info is None:
+ element = None
+ else:
+ element = info.element
+
+ inames = quadrature_inames(element)
+ inames = inames + _additional_inames(visitor)
+
+ return instruction(inames=inames, code=code, **kw)
+
+
+@preamble
+def define_quadrature_point(name):
+ local_dim = local_dimension()
+ return "Dune::FieldVector<double, {}> {};".format(local_dim, name)
+
+
+def name_quadrature_point():
+ name = "quadrature_point"
+ define_quadrature_point(name)
+ return name
+
+
+@backend(interface="quad_pos", name="sumfact")
+def pymbolic_quadrature_position():
+ visitor = get_global_context_value('visitor')
+
+ try:
+ info = visitor.current_info[1]
+ except:
+ from pudb import set_trace; set_trace()
+ if info is None:
+ element = None
+ else:
+ element = info.element
+
+ quad_point = name_quadrature_point()
+ inames = quadrature_inames(element) + _additional_inames(visitor)
+ globalarg(quad_point, shape=(local_dimension(),))
+ facedir = get_facedir(visitor.restriction)
+ for i in range(local_dimension()):
+ expression = pymbolic_indexed_quadrature_position(i, visitor)
+ instruction(assignee=prim.Subscript(prim.Variable(quad_point), i),
+ expression=expression,
+ within_inames=frozenset(inames))
+ return Variable(quad_point)
+
+
@backend(interface="qp_in_cell", name="sumfact")
def pymbolic_quadrature_position_in_cell(restriction):
- # TODO: This code path is broken at the moment.
- assert False
-
- from dune.perftool.pdelab.geometry import to_cell_coordinates
- return to_cell_coordinates(pymbolic_quadrature_position(), restriction)
+ from dune.perftool.pdelab.quadrature import pymbolic_quadrature_position_in_cell
+ return pymbolic_quadrature_position_in_cell(restriction)
diff --git a/python/dune/perftool/sumfact/switch.py b/python/dune/perftool/sumfact/switch.py
index 8c6a0f13ace27b6e3081030607149b7b003aa82c..aee3c5a1bc70e87754bbc7eb5f1c687424a401fa 100644
--- a/python/dune/perftool/sumfact/switch.py
+++ b/python/dune/perftool/sumfact/switch.py
@@ -10,7 +10,7 @@ from dune.perftool.pdelab.signatures import (assembly_routine_args,
assembly_routine_signature,
kernel_name,
)
-from dune.perftool.options import get_form_option
+from dune.perftool.options import get_form_option, get_option
from dune.perftool.cgen.clazz import ClassMember
diff --git a/python/dune/perftool/sumfact/tabulation.py b/python/dune/perftool/sumfact/tabulation.py
index 99107c9cd8f6a4429965985fb8fbbdcbd3e898e9..51e232581cb3bf0e140088b5335af876f532898d 100644
--- a/python/dune/perftool/sumfact/tabulation.py
+++ b/python/dune/perftool/sumfact/tabulation.py
@@ -20,6 +20,7 @@ from dune.perftool.generation import (class_member,
)
from dune.perftool.loopy.target import dtype_floatingpoint
from dune.perftool.loopy.vcl import ExplicitVCLCast, get_vcl_type_size
+from dune.perftool.options import get_option
from dune.perftool.pdelab.localoperator import (name_domain_field,
lop_template_range_field,
)
@@ -273,8 +274,13 @@ def local_quadrature_points_per_direction():
qps_per_dir = quadrature_points_per_direction()
if local_dimension() != world_dimension():
facedir = get_global_context_value('facedir_s')
- assert(get_global_context_value('facedir_s') == get_global_context_value('facedir_n') or
- get_global_context_value('integral_type') == 'exterior_facet')
+ if not get_option("grid_unstructured"):
+ assert(get_global_context_value('facedir_s') == get_global_context_value('facedir_n') or
+ get_global_context_value('integral_type') == 'exterior_facet')
+ if get_option("grid_unstructured"):
+ # For unstructured grids the amount of quadrature points could be different for
+ # self and neighbor. For now assert that this case is not happining.
+ assert len(set(qps_per_dir))==1
qps_per_dir = qps_per_dir[:facedir] + qps_per_dir[facedir + 1:]
return qps_per_dir
diff --git a/test/sumfact/poisson/poisson_2d_unstructured.mini b/test/sumfact/poisson/poisson_2d_unstructured.mini
index b2df475e90ee5b11b2f511d2b84df347bc31c58e..e043425520d84bf6db2cc6c0ce5d7ad631dc155b 100644
--- a/test/sumfact/poisson/poisson_2d_unstructured.mini
+++ b/test/sumfact/poisson/poisson_2d_unstructured.mini
@@ -6,6 +6,11 @@ quadvec_suffix = quadvec, nonquadvec | expand quad
gradvec_suffix = gradvec, nongradvec | expand grad
deg_suffix = deg{formcompiler.ufl_variants.degree}
+{deg_suffix} == deg1 | exclude
+{diff_suffix} == numdiff | exclude
+{quadvec_suffix} == quadvec | exclude
+
+
lowerleft = 0.0 0.0
upperright = 1.0 1.0
elements = 16 16
@@ -25,9 +30,5 @@ sumfact = 1
vectorization_quadloop = 1, 0 | expand quad
vectorization_strategy = explicit, none | expand grad
-# palpo TODO
-# diagonal_transformation_matrix = 1
-# constant_transformation_matrix = 1
-
[formcompiler.ufl_variants]
degree = 1, 2 | expand deg
diff --git a/test/sumfact/poisson/poisson_dg_2d_unstructured.mini b/test/sumfact/poisson/poisson_dg_2d_unstructured.mini
index f5c76c1352b7a269598a6289d754720a40cb411a..bd85b05834d1a87a301f37b93e65a592876afb37 100644
--- a/test/sumfact/poisson/poisson_dg_2d_unstructured.mini
+++ b/test/sumfact/poisson/poisson_dg_2d_unstructured.mini
@@ -6,6 +6,11 @@ quadvec_suffix = quadvec, nonquadvec | expand quad
gradvec_suffix = gradvec, nongradvec | expand grad
deg_suffix = deg{formcompiler.ufl_variants.degree}
+{deg_suffix} == deg1 | exclude
+{diff_suffix} == numdiff | exclude
+{quadvec_suffix} == quadvec | exclude
+
+
lowerleft = 0.0 0.0
upperright = 1.0 1.0
elements = 16 16
@@ -24,9 +29,6 @@ numerical_jacobian = 1, 0 | expand num
sumfact = 1
vectorization_quadloop = 1, 0 | expand quad
vectorization_strategy = explicit, none | expand grad
-# palpo TODO
-diagonal_transformation_matrix = 1
-constant_transformation_matrix = 1
[formcompiler.ufl_variants]
degree = 1, 2 | expand deg