From f22c5b4e96a87f08f80b279fdb57f314ae5dbdb5 Mon Sep 17 00:00:00 2001
From: Dominic Kempf <dominic.kempf@iwr.uni-heidelberg.de>
Date: Mon, 12 Dec 2016 11:35:49 +0100
Subject: [PATCH] Implement basis evaluation for symdiffs on facets
---
python/dune/perftool/sumfact/basis.py | 44 ++++++++++++++++++---------
test/sumfact/poisson/poisson_dg.mini | 1 -
2 files changed, 29 insertions(+), 16 deletions(-)
diff --git a/python/dune/perftool/sumfact/basis.py b/python/dune/perftool/sumfact/basis.py
index 794c3467..b346a2f9 100644
--- a/python/dune/perftool/sumfact/basis.py
+++ b/python/dune/perftool/sumfact/basis.py
@@ -29,7 +29,9 @@ from dune.perftool.sumfact.quadrature import quadrature_inames
from dune.perftool.sumfact.switch import (get_facedir,
get_facemod,
)
-from dune.perftool.pdelab.geometry import world_dimension
+from dune.perftool.pdelab.geometry import (local_dimension,
+ world_dimension,
+ )
from dune.perftool.loopy.buffer import initialize_buffer
from dune.perftool.pdelab.driver import FEM_name_mangling
from dune.perftool.pdelab.restriction import restricted_name
@@ -206,7 +208,7 @@ def evaluate_basis(element, name, restriction):
)
# Add the missing direction on facedirs by evaluating at either 0 or 1
- if facedir:
+ if facedir is not None:
facemod = get_facemod(restriction)
from dune.perftool.sumfact.amatrix import PolynomialLookup, name_polynomials
prod = prod + (prim.Call(PolynomialLookup(name_polynomials(), False),
@@ -238,20 +240,32 @@ def evaluate_reference_gradient(element, name, restriction):
temporary_variable(name, shape=(dim,))
quad_inames = quadrature_inames()
inames = lfs_inames(element, restriction)
- assert(len(quad_inames) == len(inames))
-
- # Matrices for sumfactorization
- theta = name_theta()
- dtheta = name_theta(derivative=True)
-
- for i in range(dim):
- calls = [prim.Subscript(prim.Variable(theta), (prim.Variable(m), prim.Variable(n)))
- for (m, n) in zip(quad_inames, inames)]
- calls[i] = prim.Subscript(prim.Variable(dtheta), (prim.Variable(quad_inames[i]), prim.Variable(inames[i])))
- calls = tuple(calls)
+ facedir = get_facedir(restriction)
- assignee = prim.Subscript(prim.Variable(name), (i,))
- expression = prim.Product(calls)
+ # Map the direction to a quadrature iname
+ quadinamemapping = {}
+ i = 0
+ for d in range(local_dimension()):
+ if d == facedir:
+ i = i+1
+ quadinamemapping[i] = quad_inames[d]
+ i = i+1
+
+ for d in range(dim):
+ prod = []
+ for i in range(dim):
+ if i != facedir:
+ prod.append(prim.Subscript(prim.Variable(name_theta(derivative=d == i)),
+ (prim.Variable(quadinamemapping[i]), prim.Variable(inames[i]))
+ ))
+ if facedir is not None:
+ facemod = get_facemod(restriction)
+ from dune.perftool.sumfact.amatrix import PolynomialLookup, name_polynomials
+ prod.append(prim.Call(PolynomialLookup(name_polynomials(), facedir==d),
+ (prim.Variable(inames[facedir]), facemod)),)
+
+ assignee = prim.Subscript(prim.Variable(name), (d,))
+ expression = prim.Product(tuple(prod))
instruction(assignee=assignee,
expression=expression,
diff --git a/test/sumfact/poisson/poisson_dg.mini b/test/sumfact/poisson/poisson_dg.mini
index c59d3c05..eeca1791 100644
--- a/test/sumfact/poisson/poisson_dg.mini
+++ b/test/sumfact/poisson/poisson_dg.mini
@@ -1,6 +1,5 @@
__name = sumfact_poisson_dg_{__exec_suffix}
__exec_suffix = numdiff, symdiff | expand num
-{__exec_suffix} == symdiff | exclude
cells = 1 1
extension = 1. 1.
--
GitLab