diff --git a/patches/apply_patches.sh b/patches/apply_patches.sh
index f95f1f6a28289d3ef0b62fdef5e818d5b87530c4..0ac07e6614e86c9a439584a71d6f17b0a1a8a90d 100755
--- a/patches/apply_patches.sh
+++ b/patches/apply_patches.sh
@@ -12,3 +12,7 @@ popd
pushd python/ufl
git apply ../../patches/ufl/conditional-uflid.patch
popd
+
+pushd python/ufl
+git apply ../../patches/ufl/tensor-product-element.patch
+popd
diff --git a/patches/ufl/tensor-product-element.patch b/patches/ufl/tensor-product-element.patch
new file mode 100644
index 0000000000000000000000000000000000000000..9fc64f124e95bcbf28391ed5be4c87e1040a4e28
--- /dev/null
+++ b/patches/ufl/tensor-product-element.patch
@@ -0,0 +1,19 @@
+commit f87dcd18d765b0200808b79b2e7374f82a0c6199
+Author: René Heß <rene.hess@iwr.uni-heidelberg.de>
+Date: Tue Aug 29 14:56:17 2017 +0200
+
+ Patch for TensorProductElements
+
+diff --git a/ufl/algorithms/compute_form_data.py b/ufl/algorithms/compute_form_data.py
+index 3388bbfc..1cef3924 100644
+--- a/ufl/algorithms/compute_form_data.py
++++ b/ufl/algorithms/compute_form_data.py
+@@ -56,7 +56,7 @@ def _auto_select_degree(elements):
+ """
+ # Use max degree of all elements, at least 1 (to work with
+ # Lagrange elements)
+- return max({e.degree() for e in elements} - {None} | {1})
++ return max({e.degree() if not isinstance(e.degree(), tuple) else max(e.degree()) for e in elements} - {None} | {1})
+
+
+ def _compute_element_mapping(form):
diff --git a/python/dune/perftool/generation/loopy.py b/python/dune/perftool/generation/loopy.py
index 171ba6709c27fa8e88876da6180ef1b9debd26f4..fb3dbaba5ab97dbc7f3f40b1bae02bf15777695e 100644
--- a/python/dune/perftool/generation/loopy.py
+++ b/python/dune/perftool/generation/loopy.py
@@ -50,10 +50,21 @@ def valuearg(name, **kw):
@generator_factory(item_tags=("domain",), context_tags="kernel")
def domain(iname, shape):
- if isinstance(iname, tuple):
+ if isinstance(iname, tuple) and isinstance(shape, tuple):
+ assert(len(iname) == len(shape))
+ condition = ""
+ for index, (i, s) in enumerate(zip(iname, shape)):
+ if index > 0:
+ condition += " and "
+ else:
+ condition += " "
+ condition += "0<={}<{}".format(i, s)
iname = ",".join(iname)
- if isinstance(shape, str):
- valuearg(shape)
+ return "{{ [{}] : {} }}".format(iname, condition)
+ elif isinstance(iname, tuple):
+ iname = ",".join(iname)
+ if isinstance(shape, str):
+ valuearg(shape)
return "{{ [{0}] : 0<={0}<{1} }}".format(iname, shape)
diff --git a/python/dune/perftool/loopy/transformations/vectorize_quad.py b/python/dune/perftool/loopy/transformations/vectorize_quad.py
index 149002c4b872ba28f21e790a4df6c2a044ae5b40..abc7294086f52f9281090fcbe36de59bfd162ce4 100644
--- a/python/dune/perftool/loopy/transformations/vectorize_quad.py
+++ b/python/dune/perftool/loopy/transformations/vectorize_quad.py
@@ -136,6 +136,7 @@ def _vectorize_quadrature_loop(knl, inames, suffix):
# Do precomputation of the quantity
prec_quantity = "{}_precomputed".format(quantity)
prec_quantities.append(prec_quantity)
+
knl = lp.precompute(knl, subst_name, inames,
temporary_name=prec_quantity,
)
diff --git a/python/dune/perftool/options.py b/python/dune/perftool/options.py
index 49e6941a443493ac7d82294bc461a2a79b5c6c1d..ed37897178ab6fa5b74b987e38ce81f341eefa9a 100644
--- a/python/dune/perftool/options.py
+++ b/python/dune/perftool/options.py
@@ -63,8 +63,8 @@ class PerftoolOptionsArray(ImmutableRecord):
# Arguments that are mainly to be set by logic depending on other options
max_vector_width = PerftoolOption(default=256, helpstr=None)
- unroll_dimension_loops = PerftoolOption(default=False, helpstr="whether loops over the gemetric dimension should be unrolled.")
- precompute_quadrature_info = PerftoolOption(default=True, helpstr="whether loops over the gemetric dimension should be unrolled.")
+ unroll_dimension_loops = PerftoolOption(default=False, helpstr="whether loops over the geometric dimension should be unrolled")
+ precompute_quadrature_info = PerftoolOption(default=True, helpstr="compute quadrature points and weights in the constructor of the local operator")
blockstructured = PerftoolOption(default=False, helpstr="Use block structure")
number_of_blocks = PerftoolOption(default=1, helpstr="Number of sub blocks in one direction")
diff --git a/python/dune/perftool/pdelab/driver/__init__.py b/python/dune/perftool/pdelab/driver/__init__.py
index 81ebce1ea48efc25b88d9e3863df15535afc96db..2c38881f4a7351a199bc43f1515d8fc55c9914d2 100644
--- a/python/dune/perftool/pdelab/driver/__init__.py
+++ b/python/dune/perftool/pdelab/driver/__init__.py
@@ -136,7 +136,7 @@ def isDG(fem):
def FEM_name_mangling(fem):
- from ufl import MixedElement, VectorElement, FiniteElement, TensorElement
+ from ufl import MixedElement, VectorElement, FiniteElement, TensorElement, TensorProductElement
if isinstance(fem, VectorElement):
return FEM_name_mangling(fem.sub_elements()[0]) + "_" + str(fem.num_sub_elements())
if isinstance(fem, TensorElement):
@@ -155,6 +155,14 @@ def FEM_name_mangling(fem):
return "Q" + str(fem._degree)
if isDG(fem):
return "DG" + str(fem._degree)
+ if isinstance(fem, TensorProductElement):
+ assert(len(set(subel._short_name for subel in fem.sub_elements())) == 1)
+
+ if isLagrange(fem.sub_elements()[0]):
+ return "TensorQ" + '_'.join(map(str, fem._degree))
+ if isDG(fem.sub_elements()[0]):
+ return "TensorDG" + '_'.join(map(str, fem._degree))
+ raise NotImplementedError("fem name mangling")
raise NotImplementedError("FEM NAME MANGLING")
diff --git a/python/dune/perftool/pdelab/driver/constraints.py b/python/dune/perftool/pdelab/driver/constraints.py
index 2052532851c119d49bcb47284cb0e40531f3b573..9158fc6be52d5c95c976eb08041fc3ad5b4a5da6 100644
--- a/python/dune/perftool/pdelab/driver/constraints.py
+++ b/python/dune/perftool/pdelab/driver/constraints.py
@@ -14,7 +14,7 @@ from dune.perftool.pdelab.driver.gridfunctionspace import (name_gfs,
preprocess_leaf_data,
)
-from ufl import FiniteElement, MixedElement, TensorElement, VectorElement
+from ufl import FiniteElement, MixedElement, TensorElement, VectorElement, TensorProductElement
def name_assembled_constraints():
@@ -53,7 +53,7 @@ def name_bctype_function(element, is_dirichlet):
define_composite_bctype_function(element, is_dirichlet, name, tuple(childs))
return name
else:
- assert isinstance(element, FiniteElement)
+ assert isinstance(element, (FiniteElement, TensorProductElement))
name = "{}_bctype".format(FEM_name_mangling(element).lower())
define_bctype_function(element, is_dirichlet[0], name)
return name
diff --git a/python/dune/perftool/pdelab/driver/gridfunctionspace.py b/python/dune/perftool/pdelab/driver/gridfunctionspace.py
index e442c08ed9fc7b5f7b12f7505fdbd86287918bd3..bf47920a94513d0075b92d3b65bc8357081e371e 100644
--- a/python/dune/perftool/pdelab/driver/gridfunctionspace.py
+++ b/python/dune/perftool/pdelab/driver/gridfunctionspace.py
@@ -7,6 +7,7 @@ from dune.perftool.pdelab.driver import (FEM_name_mangling,
get_dimension,
get_test_element,
get_trial_element,
+ isLagrange,
isDG,
isPk,
isQk,
@@ -16,7 +17,7 @@ from dune.perftool.pdelab.driver import (FEM_name_mangling,
preprocess_leaf_data,
)
-from ufl import FiniteElement, MixedElement, TensorElement, VectorElement
+from ufl import FiniteElement, MixedElement, TensorElement, VectorElement, TensorProductElement
@preamble
@@ -124,24 +125,52 @@ def typedef_fem(element, name):
df = type_domainfield()
r = type_range()
dim = get_dimension()
+
if get_option("blockstructured"):
include_file("dune/perftool/blockstructured/blockstructuredqkfem.hh", filetag="driver")
degree = element.degree() * get_option("number_of_blocks")
- return "using {} = Dune::PDELab::BlockstructuredQkLocalFiniteElementMap<{}, {}, {}, {}>;".format(name, gv, df, r, degree)
- if isQk(element):
+ return "using {} = Dune::PDELab::BlockstructuredQkLocalFiniteElementMap<{}, {}, {}, {}>;" \
+ .format(name, gv, df, r, degree)
+
+ if isinstance(element, TensorProductElement):
+ # Only allow TensorProductElements where all subelements are
+ # of the same type ('CG' or 'DG')
+ assert(len(set(subel._short_name for subel in element.sub_elements())) == 1)
+
+ # Anisotropic degree is not yet supported in Dune
+ degrees = element.degree()
+ for deg in degrees:
+ assert (deg == degrees[0])
+
+ # TensorProductElements have Qk structure -> no Pk
+ if isLagrange(element.sub_elements()[0]):
+ include_file("dune/pdelab/finiteelementmap/qkfem.hh", filetag="driver")
+ return "using {} = Dune::PDELab::QkLocalFiniteElementMap<{}, {}, {}, {}>;" \
+ .format(name, gv, df, r, degrees[0])
+ elif isDG(element.sub_elements()[0]):
+ include_file("dune/pdelab/finiteelementmap/qkdg.hh", filetag="driver")
+ # TODO allow switching the basis here!
+ return "using {} = Dune::PDELab::QkDGLocalFiniteElementMap<{}, {}, {}, {}>;" \
+ .format(name, df, r, degrees[0], dim)
+ raise NotImplementedError("FEM not implemented in dune-perftool")
+ elif isQk(element):
include_file("dune/pdelab/finiteelementmap/qkfem.hh", filetag="driver")
- return "using {} = Dune::PDELab::QkLocalFiniteElementMap<{}, {}, {}, {}>;".format(name, gv, df, r, element.degree())
- if isPk(element):
+ return "using {} = Dune::PDELab::QkLocalFiniteElementMap<{}, {}, {}, {}>;" \
+ .format(name, gv, df, r, element.degree())
+ elif isPk(element):
include_file("dune/pdelab/finiteelementmap/pkfem.hh", filetag="driver")
- return "using {} = Dune::PDELab::PkLocalFiniteElementMap<{}, {}, {}, {}>;".format(name, gv, df, r, element.degree())
- if isDG(element):
+ return "using {} = Dune::PDELab::PkLocalFiniteElementMap<{}, {}, {}, {}>;" \
+ .format(name, gv, df, r, element.degree())
+ elif isDG(element):
if isQuadrilateral(element.cell()):
include_file("dune/pdelab/finiteelementmap/qkdg.hh", filetag="driver")
# TODO allow switching the basis here!
- return "using {} = Dune::PDELab::QkDGLocalFiniteElementMap<{}, {}, {}, {}>;".format(name, df, r, element.degree(), dim)
+ return "using {} = Dune::PDELab::QkDGLocalFiniteElementMap<{}, {}, {}, {}>;" \
+ .format(name, df, r, element.degree(), dim)
if isSimplical(element.cell()):
include_file("dune/pdelab/finiteelementmap/opbfem.hh", filetag="driver")
- return "using {} = Dune::PDELab::OPBLocalFiniteElementMap<{}, {}, {}, {}, Dune::GeometryType::simplex>;".format(name, df, r, element.degree(), dim)
+ return "using {} = Dune::PDELab::OPBLocalFiniteElementMap<{}, {}, {}, {}, Dune::GeometryType::simplex>;" \
+ .format(name, df, r, element.degree(), dim)
raise NotImplementedError("Geometry type not known in code generation")
raise NotImplementedError("FEM not implemented in dune-perftool")
@@ -157,15 +186,26 @@ def type_fem(element):
def define_fem(element, name):
femtype = type_fem(element)
from dune.perftool.pdelab.driver import isDG
- if isDG(element):
+ if isinstance(element, TensorProductElement):
+ # Only allow TensorProductElements where all subelements are
+ # of the same type ('CG' or 'DG')
+ assert(len(set(subel._short_name for subel in element.sub_elements())) == 1)
+ if isDG(element.sub_elements()[0]):
+ return "{} {};".format(femtype, name)
+ else:
+ assert(isLagrange(element.sub_elements()[0]))
+ gv = name_leafview()
+ return "{} {}({});".format(femtype, name, gv)
+ elif isDG(element):
return "{} {};".format(femtype, name)
else:
+ assert(isLagrange(element))
gv = name_leafview()
return "{} {}({});".format(femtype, name, gv)
def name_fem(element):
- assert isinstance(element, FiniteElement)
+ assert isinstance(element, (FiniteElement, TensorProductElement))
name = "{}_fem".format(FEM_name_mangling(element).lower())
define_fem(element, name)
return name
@@ -192,7 +232,7 @@ def name_gfs(element, is_dirichlet, treepath=()):
"_".join(str(t) for t in treepath))
define_power_gfs(element, is_dirichlet, name, subgfs)
return name
- if isinstance(element, MixedElement):
+ elif isinstance(element, MixedElement):
k = 0
subgfs = []
for i, subel in enumerate(element.sub_elements()):
@@ -203,7 +243,7 @@ def name_gfs(element, is_dirichlet, treepath=()):
define_composite_gfs(element, is_dirichlet, name, tuple(subgfs))
return name
else:
- assert isinstance(element, FiniteElement)
+ assert isinstance(element, (FiniteElement, TensorProductElement))
name = "{}{}_gfs_{}".format(FEM_name_mangling(element).lower(),
"_dirichlet" if is_dirichlet[0] else "",
"_".join(str(t) for t in treepath))
@@ -230,7 +270,7 @@ def type_gfs(element, is_dirichlet):
name = "{}_POW{}GFS".format(subgfs, element.num_sub_elements())
typedef_power_gfs(element, is_dirichlet, name, subgfs)
return name
- if isinstance(element, MixedElement):
+ elif isinstance(element, MixedElement):
k = 0
subgfs = []
for subel in element.sub_elements():
@@ -240,7 +280,7 @@ def type_gfs(element, is_dirichlet):
typedef_composite_gfs(element, name, tuple(subgfs))
return name
else:
- assert isinstance(element, FiniteElement)
+ assert isinstance(element, (FiniteElement, TensorProductElement))
name = "{}{}_GFS".format(FEM_name_mangling(element).upper(),
"_dirichlet" if is_dirichlet[0] else "",
)
diff --git a/python/dune/perftool/pdelab/driver/interpolate.py b/python/dune/perftool/pdelab/driver/interpolate.py
index cacb75aa1642180a9c8a014380e46b842faed09a..f3cc2afbeb3f083e41ed0c66a4060c08a2b0e5a1 100644
--- a/python/dune/perftool/pdelab/driver/interpolate.py
+++ b/python/dune/perftool/pdelab/driver/interpolate.py
@@ -15,7 +15,7 @@ from dune.perftool.pdelab.driver.gridfunctionspace import (name_trial_gfs,
)
from dune.perftool.pdelab.driver.gridoperator import (name_parameters,)
-from ufl import FiniteElement, MixedElement, TensorElement, VectorElement
+from ufl import FiniteElement, MixedElement, TensorElement, VectorElement, TensorProductElement
def _do_interpolate(dirichlet):
@@ -55,7 +55,7 @@ def name_boundary_function(element, func):
define_compositegfs_parameterfunction(name, tuple(childs))
return name
else:
- assert isinstance(element, FiniteElement)
+ assert isinstance(element, (FiniteElement, TensorProductElement))
name = get_counted_variable("func")
define_boundary_function(name, func[0])
return name
diff --git a/python/dune/perftool/pdelab/driver/vtk.py b/python/dune/perftool/pdelab/driver/vtk.py
index 7edf92f357533393c9fb81370743904808885133..824f24ded82a193caaa924e6bfee5c4fa552faf3 100644
--- a/python/dune/perftool/pdelab/driver/vtk.py
+++ b/python/dune/perftool/pdelab/driver/vtk.py
@@ -41,7 +41,10 @@ def type_vtkwriter():
@preamble
def define_subsamplinglevel(name):
ini = name_initree()
- return "int {} = {}.get<int>(\"vtk.subsamplinglevel\", {});".format(name, ini, max(get_trial_element().degree() - 1, 0))
+ degree = get_trial_element().degree()
+ if isinstance(degree, tuple):
+ degree = max(degree)
+ return "int {} = {}.get<int>(\"vtk.subsamplinglevel\", {});".format(name, ini, max(degree - 1, 0))
def name_subsamplinglevel():
diff --git a/python/dune/perftool/pdelab/localoperator.py b/python/dune/perftool/pdelab/localoperator.py
index a1f8e04417e2caae95135090578dc0fb9f69e01b..90b320db0cfb13fabfc2ab31f7389042145ea46e 100644
--- a/python/dune/perftool/pdelab/localoperator.py
+++ b/python/dune/perftool/pdelab/localoperator.py
@@ -775,9 +775,9 @@ def generate_localoperator_kernels(formdata, data):
operator_kernels[(measure, 'residual')] = kernel
- logger.info("generate_localoperator_kernels: create jacobian methods")
# Generate the necessary jacobian methods
if not get_option("numerical_jacobian"):
+ logger.info("generate_localoperator_kernels: create jacobian methods")
from ufl import derivative
jacform = derivative(formdata.original_form, formdata.original_form.coefficients()[0])
diff --git a/python/dune/perftool/pdelab/quadrature.py b/python/dune/perftool/pdelab/quadrature.py
index 8356ecd22e28809825f334e03f6cf4e58734dbdb..4f28b9c52c8a2fe51b0b31ce7c34145ea63a7985 100644
--- a/python/dune/perftool/pdelab/quadrature.py
+++ b/python/dune/perftool/pdelab/quadrature.py
@@ -189,9 +189,25 @@ def _estimate_quadrature_order():
# Estimate polynomial degree of integrals of current type (eg 'Cell')
integral_type = get_global_context_value("integral_type")
integrals = form.integrals_by_type(integral_type)
- polynomial_degree = 0
+
+ # Degree could be a tuple (for TensorProductElements)
+ degree = integrals[0].metadata()['estimated_polynomial_degree']
+ if isinstance(degree, int):
+ degree = (degree,)
+ polynomial_degree = [0, ] * len(degree)
+
for i in integrals:
- polynomial_degree = max(polynomial_degree, i.metadata()['estimated_polynomial_degree'])
+ degree = i.metadata()['estimated_polynomial_degree']
+ if isinstance(degree, int):
+ degree = [degree, ]
+ assert(len(polynomial_degree) == len(degree))
+ for i in range(len(polynomial_degree)):
+ polynomial_degree[i] = max(polynomial_degree[i], degree[i])
+
+ # Return either a tuple or an int
+ polynomial_degree = tuple(polynomial_degree)
+ if len(polynomial_degree) == 1:
+ polynomial_degree = polynomial_degree[0]
return polynomial_degree
@@ -199,14 +215,31 @@ def _estimate_quadrature_order():
def quadrature_order():
"""Get quadrature order
- Note: In PDELab quadrature order m means that integration of
- polynomials of degree m is excat.
+ Notes:
+
+ - In PDELab quadrature order m means that integration of
+ polynomials of degree m is exact.
+
+ - If you use sum factorization and TensorProductElement it is
+ possible to use a different quadrature_order per direction.
"""
if get_option("quadrature_order"):
- quadrature_order = int(get_option("quadrature_order"))
+ quadrature_order = tuple(map(int, get_option("quadrature_order").split(',')))
else:
quadrature_order = _estimate_quadrature_order()
+ # TensorProductElements can have different quadrature order for
+ # every direction so quadrature_order may be a tuple. If we do not
+ # use TensorProductElements we want to return an int.
+ if isinstance(quadrature_order, tuple):
+ if len(quadrature_order) == 1:
+ quadrature_order = quadrature_order[0]
+ if isinstance(quadrature_order, tuple):
+ if not get_option('sumfact'):
+ raise NotImplementedError("Different quadrature order per direction is only implemented for kernels using sum factorization.")
+ from dune.perftool.pdelab.geometry import world_dimension
+ assert(len(quadrature_order) == world_dimension())
+
return quadrature_order
diff --git a/python/dune/perftool/pdelab/spaces.py b/python/dune/perftool/pdelab/spaces.py
index c9883dd67763fdccd37a25859bdf1804423ce6d1..706be5b91a186b990719227fe61f1336d01313a0 100644
--- a/python/dune/perftool/pdelab/spaces.py
+++ b/python/dune/perftool/pdelab/spaces.py
@@ -6,6 +6,7 @@ from dune.perftool.generation import (class_member,
generator_factory,
include_file,
preamble,
+ valuearg,
)
from dune.perftool.pdelab.restriction import restricted_name
from dune.perftool.ufl.modified_terminals import Restriction
@@ -73,6 +74,7 @@ def define_lfs(name, father, child):
@preamble
def define_lfs_size(lfs, element, restriction):
name = name_lfs_bound(lfs)
+ valuearg(name, dtype=numpy.int32)
return "auto {} = {}.size();".format(name, lfs)
diff --git a/python/dune/perftool/sumfact/accumulation.py b/python/dune/perftool/sumfact/accumulation.py
index fa48689f1034952acf21fcab15ca88775d24ea8f..50f4d169935cae6ec421da62ba432506d055f298 100644
--- a/python/dune/perftool/sumfact/accumulation.py
+++ b/python/dune/perftool/sumfact/accumulation.py
@@ -41,6 +41,7 @@ import loopy as lp
import numpy as np
import pymbolic.primitives as prim
import ufl.classes as uc
+from ufl import FiniteElement, MixedElement, TensorProductElement
@iname
@@ -138,8 +139,7 @@ def get_accumulation_info(expr, visitor):
def _get_childs(element):
- from ufl import FiniteElement
- if isinstance(element, FiniteElement):
+ if isinstance(element, (FiniteElement, TensorProductElement)):
yield (0, element)
else:
for i in range(element.value_size()):
@@ -199,10 +199,18 @@ def generate_accumulation_instruction(expr, visitor):
test_info = visitor.test_info
trial_info = visitor.trial_info
+ # Number of basis functions per direction
leaf_element = test_info.element
- if leaf_element.num_sub_elements() > 0:
+ from ufl import MixedElement
+ if isinstance(leaf_element, MixedElement):
leaf_element = leaf_element.extract_component(test_info.element_index)[1]
- basis_size = leaf_element.degree() + 1
+ degree = leaf_element._degree
+ if isinstance(degree, int):
+ degree = (degree,) * dim
+ basis_size = tuple(deg + 1 for deg in degree)
+
+ # Anisotropic finite elements are not (yet) supported by Dune
+ assert(size == basis_size[0] for size in basis_size)
from dune.perftool.pdelab.localoperator import boundary_predicates
predicates = boundary_predicates(expr,
@@ -316,7 +324,7 @@ def generate_accumulation_instruction(expr, visitor):
# Collect the lfs and lfs indices for the accumulate call
test_lfs.index = flatten_index(tuple(prim.Variable(i) for i in inames),
- (basis_size,) * dim,
+ basis_size,
order="f"
)
@@ -327,7 +335,7 @@ def generate_accumulation_instruction(expr, visitor):
from dune.perftool.sumfact.basis import _basis_functions_per_direction
ansatz_lfs.index = flatten_index(tuple(prim.Variable(jacobian_inames[i])
for i in range(world_dimension())),
- (_basis_functions_per_direction(trial_leaf_element),) * dim,
+ _basis_functions_per_direction(trial_leaf_element),
order="f"
)
diff --git a/python/dune/perftool/sumfact/basis.py b/python/dune/perftool/sumfact/basis.py
index a58a7bb7df83123d5756d2ba4b9bb8f7a3d0be33..8e745c7daf65294db5e9f34b03c7a19736918f8c 100644
--- a/python/dune/perftool/sumfact/basis.py
+++ b/python/dune/perftool/sumfact/basis.py
@@ -42,7 +42,7 @@ from dune.perftool.tools import maybe_wrap_subscript
from dune.perftool.pdelab.basis import shape_as_pymbolic
from dune.perftool.sumfact.accumulation import sumfact_iname
-from ufl import VectorElement, TensorElement
+from ufl import MixedElement, VectorElement, TensorElement, TensorProductElement
from pytools import product, ImmutableRecord
@@ -96,20 +96,30 @@ class LFSSumfactKernelInput(SumfactKernelInputBase, ImmutableRecord):
def _basis_functions_per_direction(element):
"""Number of basis functions per direction """
- from ufl import FiniteElement
- assert isinstance(element, FiniteElement)
- return element.degree() + 1
+ from ufl import FiniteElement, TensorProductElement
+ assert isinstance(element, (FiniteElement, TensorProductElement))
+ degree = element.degree()
+ if isinstance(degree, int):
+ degree = (degree,) * world_dimension()
+
+ basis_size = tuple(deg + 1 for deg in degree)
+
+ # Anisotropic finite elements are not (yet) supported by Dune
+ for size in basis_size:
+ assert(size == basis_size[0])
+
+ return basis_size
@kernel_cached
def pymbolic_coefficient_gradient(element, restriction, index, coeff_func, visitor_indices):
sub_element = element
grad_index = visitor_indices[0]
- if element.num_sub_elements() > 0:
+ if isinstance(element, MixedElement):
sub_element = element.extract_component(index)[1]
from ufl import FiniteElement
- assert isinstance(sub_element, FiniteElement)
+ assert isinstance(sub_element, (FiniteElement, TensorProductElement))
# Number of basis functions per direction
basis_size = _basis_functions_per_direction(sub_element)
@@ -182,6 +192,7 @@ def pymbolic_coefficient(element, restriction, index, coeff_func, visitor_indice
@iname
def sumfact_lfs_iname(element, bound, dim):
+ assert(isinstance(bound, int))
from dune.perftool.pdelab.driver import FEM_name_mangling
name = "sumfac_lfs_{}_{}".format(FEM_name_mangling(element), dim)
domain(name, bound)
@@ -190,13 +201,13 @@ def sumfact_lfs_iname(element, bound, dim):
@backend(interface="lfs_inames", name="sumfact")
def lfs_inames(element, restriction, number=1, context=''):
- from ufl import FiniteElement
- assert isinstance(element, FiniteElement)
+ from ufl import FiniteElement, TensorProductElement
+ assert isinstance(element, (FiniteElement, TensorProductElement))
if number == 0:
return ()
else:
dim = world_dimension()
- return tuple(sumfact_lfs_iname(element, _basis_functions_per_direction(element), d) for d in range(dim))
+ return tuple(sumfact_lfs_iname(element, _basis_functions_per_direction(element)[d], d) for d in range(dim))
@backend(interface="evaluate_basis", name="sumfact")
@@ -209,10 +220,18 @@ def evaluate_basis(element, name, restriction):
# Collect the pairs of lfs/quad inames that are in use
# On facets, the normal direction of the facet is excluded
- tab = BasisTabulationMatrix(basis_size=_basis_functions_per_direction(element))
- prod = tuple(tab.pymbolic((prim.Variable(i), prim.Variable(j)))
- for (i, j) in zip(quad_inames, tuple(iname for i, iname in enumerate(inames) if i != facedir))
- )
+ #
+ # If facedir is not none, the length of inames and quad_inames is
+ # different. For inames we want to skip the facedir direction, for
+ # quad_inames we need all entries. Thats the reason for the
+ # help_index.
+ basis_per_dir = _basis_functions_per_direction(element)
+ prod = ()
+ help_index = 0
+ for direction in range(len(inames)):
+ if direction != facedir:
+ prod = prod + (BasisTabulationMatrix(basis_size=basis_per_dir[direction]).pymbolic((prim.Variable(quad_inames[help_index]), prim.Variable(inames[direction]))),)
+ help_index += 1
# Add the missing direction on facedirs by evaluating at either 0 or 1
if facedir is not None:
@@ -263,7 +282,7 @@ def evaluate_reference_gradient(element, name, restriction, index):
prod = []
for i in range(dim):
if i != facedir:
- tab = BasisTabulationMatrix(basis_size=_basis_functions_per_direction(element),
+ tab = BasisTabulationMatrix(basis_size=_basis_functions_per_direction(element)[i],
derivative=index == i)
prod.append(tab.pymbolic((prim.Variable(quadinamemapping[i]), prim.Variable(inames[i]))))
diff --git a/python/dune/perftool/sumfact/quadrature.py b/python/dune/perftool/sumfact/quadrature.py
index a5ede367812b2c0c50153538d7bb893f8c077b37..d42296f6aa533049402a6d8f3f45695cddf90339 100644
--- a/python/dune/perftool/sumfact/quadrature.py
+++ b/python/dune/perftool/sumfact/quadrature.py
@@ -10,13 +10,16 @@ from dune.perftool.generation import (backend,
)
from dune.perftool.sumfact.tabulation import (quadrature_points_per_direction,
+ local_quadrature_points_per_direction,
name_oned_quadrature_points,
name_oned_quadrature_weights,
)
from dune.perftool.pdelab.argument import name_accumulation_variable
from dune.perftool.pdelab.geometry import (local_dimension,
+ world_dimension,
)
from dune.perftool.options import get_option
+from dune.perftool.sumfact.switch import get_facedir
from loopy import CallMangleInfo
from loopy.symbolic import FunctionIdentifier
@@ -80,52 +83,59 @@ def quadrature_inames(element):
if element is None:
names = tuple("quad_{}".format(d) for d in range(local_dimension()))
else:
- from ufl import FiniteElement
- assert isinstance(element, FiniteElement)
+ from ufl import FiniteElement, TensorProductElement
+ assert isinstance(element, (FiniteElement, TensorProductElement))
from dune.perftool.pdelab.driver import FEM_name_mangling
names = tuple("quad_{}_{}".format(FEM_name_mangling(element), d) for d in range(local_dimension()))
- domain(names, quadrature_points_per_direction())
+
+ local_qps_per_dir = local_quadrature_points_per_direction()
+ domain(names, local_qps_per_dir)
return names
@iname(kernel="operator")
def constructor_quad_iname(name, d, bound):
- name = "{}_{}".format(name, d)
- domain(name, quadrature_points_per_direction(), kernel="operator")
+ name = "{}_localdim{}".format(name, d)
+ domain(name, bound, kernel="operator")
return name
-def constructor_quadrature_inames(dim, num1d):
- name = "quadiname_dim{}_num{}".format(dim, num1d)
- return tuple(constructor_quad_iname(name, d, quadrature_points_per_direction()) for d in range(local_dimension()))
+def constructor_quadrature_inames(dim):
+ local_qps_per_dir = local_quadrature_points_per_direction()
+ local_qps_per_dir_str = '_'.join(map(str, local_qps_per_dir))
+ name = "quadiname_dim{}_num{}".format(dim, local_qps_per_dir_str)
+ return tuple(constructor_quad_iname(name, d, local_qps_per_dir[d]) for d in range(local_dimension()))
-def define_recursive_quadrature_weight(visitor, name, dir):
+def define_recursive_quadrature_weight(visitor, name, direction):
info = visitor.current_info[1]
if info is None:
element = None
else:
element = info.element
- iname = quadrature_inames(element)[dir]
+ iname = quadrature_inames(element)[direction]
temporary_variable(name, shape=(), shape_impl=())
- instruction(expression=Product((recursive_quadrature_weight(dir + 1),
- Subscript(Variable(name_oned_quadrature_weights()),
+
+ qps_per_dir = quadrature_points_per_direction()
+ qp_bound = qps_per_dir[direction]
+ instruction(expression=Product((recursive_quadrature_weight(direction + 1),
+ Subscript(Variable(name_oned_quadrature_weights(qp_bound)),
(Variable(iname),)
))
),
assignee=Variable(name),
- forced_iname_deps=frozenset(quadrature_inames()[dir:]),
+ forced_iname_deps=frozenset(quadrature_inames()[direction:]),
forced_iname_deps_is_final=True,
tags=frozenset({"quad"}),
)
-def recursive_quadrature_weight(visitor, dir=0):
- if dir == local_dimension():
+def recursive_quadrature_weight(visitor, direction=0):
+ if direction == local_dimension():
return pymbolic_base_weight()
else:
- name = 'weight_{}'.format(dir)
- define_recursive_quadrature_weight(visitor, name, dir)
+ name = 'weight_{}'.format(direction)
+ define_recursive_quadrature_weight(visitor, name, direction)
return Variable(name)
@@ -134,14 +144,16 @@ def quadrature_weight(visitor):
if not get_option("precompute_quadrature_info"):
return recursive_quadrature_weight(visitor)
+ # Quadrature points per (local) direction
dim = local_dimension()
- num1d = quadrature_points_per_direction()
- name = "quad_weights_dim{}_num{}".format(dim, num1d)
+ local_qps_per_dir = local_quadrature_points_per_direction()
+ local_qps_per_dir_str = '_'.join(map(str, local_qps_per_dir))
+ name = "quad_weights_dim{}_num{}".format(dim, local_qps_per_dir_str)
# Add a class member
loopy_class_member(name,
dtype=np.float64,
- shape=(num1d,) * dim,
+ shape=local_qps_per_dir,
classtag="operator",
dim_tags=",".join(["f"] * dim),
managed=True,
@@ -149,9 +161,12 @@ def quadrature_weight(visitor):
)
# Precompute it in the constructor
- inames = constructor_quadrature_inames(dim, num1d)
- instruction(assignee=prim.Subscript(prim.Variable(name), tuple(prim.Variable(i) for i in inames)),
- expression=prim.Product(tuple(Subscript(Variable(name_oned_quadrature_weights()), (prim.Variable(i),)) for i in inames)),
+ inames = constructor_quadrature_inames(dim)
+ assignee = prim.Subscript(prim.Variable(name), tuple(prim.Variable(i) for i in inames))
+ expression = prim.Product(tuple(Subscript(Variable(name_oned_quadrature_weights(local_qps_per_dir[index])),
+ (prim.Variable(iname),)) for index, iname in enumerate(inames)))
+ instruction(assignee=assignee,
+ expression=expression,
within_inames=frozenset(inames),
kernel="operator",
)
@@ -165,7 +180,11 @@ def quadrature_weight(visitor):
def define_quadrature_position(name, index):
- instruction(expression=Subscript(Variable(name_oned_quadrature_points()), (Variable(quadrature_inames()[index]),)),
+ qps_per_dir = quadrature_points_per_direction()
+ qp_bound = qps_per_dir[index]
+ expression = Subscript(Variable(name_oned_quadrature_points(qp_bound)),
+ (Variable(quadrature_inames()[index]),))
+ instruction(expression=expression,
assignee=Subscript(Variable(name), (index,)),
forced_iname_deps=frozenset(quadrature_inames(element)),
forced_iname_deps_is_final=True,
@@ -184,12 +203,14 @@ def pymbolic_quadrature_position(index, visitor):
assert isinstance(index, int)
dim = local_dimension()
- num1d = quadrature_points_per_direction()
- name = "quad_points_dim{}_num{}_dir{}".format(dim, num1d, index)
+ 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)
loopy_class_member(name,
dtype=np.float64,
- shape=(num1d,) * dim,
+ shape=local_qps_per_dir,
classtag="operator",
dim_tags=",".join(["f"] * dim),
managed=True,
@@ -197,9 +218,12 @@ def pymbolic_quadrature_position(index, visitor):
)
# Precompute it in the constructor
- inames = constructor_quadrature_inames(dim, num1d)
- instruction(assignee=prim.Subscript(prim.Variable(name), tuple(prim.Variable(i) for i in inames)),
- expression=Subscript(Variable(name_oned_quadrature_points()), (prim.Variable(inames[index]))),
+ 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(qps_per_dir[index])),
+ (prim.Variable(inames[index])))
+ instruction(assignee=assignee,
+ expression=expression,
within_inames=frozenset(inames),
kernel="operator",
)
diff --git a/python/dune/perftool/sumfact/symbolic.py b/python/dune/perftool/sumfact/symbolic.py
index 288374daeebaa7780862bb33e2565486e2f7c3f1..a0cc1054afac0a67f312667f7fa4243700c1ccac 100644
--- a/python/dune/perftool/sumfact/symbolic.py
+++ b/python/dune/perftool/sumfact/symbolic.py
@@ -4,7 +4,7 @@ from dune.perftool.options import get_option
from dune.perftool.generation import get_counted_variable
from dune.perftool.pdelab.geometry import local_dimension, world_dimension
from dune.perftool.sumfact.quadrature import quadrature_inames
-from dune.perftool.sumfact.tabulation import BasisTabulationMatrixArray
+from dune.perftool.sumfact.tabulation import BasisTabulationMatrixBase, BasisTabulationMatrixArray
from pytools import ImmutableRecord, product
@@ -107,7 +107,6 @@ class SumfactKernel(SumfactKernelBase, ImmutableRecord, prim.Variable):
"""
# Assert the inputs!
assert isinstance(matrix_sequence, tuple)
- from dune.perftool.sumfact.tabulation import BasisTabulationMatrixBase
assert all(isinstance(m, BasisTabulationMatrixBase) for m in matrix_sequence)
assert stage in (1, 3)
diff --git a/python/dune/perftool/sumfact/tabulation.py b/python/dune/perftool/sumfact/tabulation.py
index 736f64b282a0bce3a06c4d682e55213ebac9f1ea..e8518668c227c7528ab72712bbd9e2cba92c147f 100644
--- a/python/dune/perftool/sumfact/tabulation.py
+++ b/python/dune/perftool/sumfact/tabulation.py
@@ -3,7 +3,7 @@ from dune.perftool.ufl.modified_terminals import Restriction
from dune.perftool.options import get_option
from dune.perftool.pdelab.argument import name_coefficientcontainer
-from dune.perftool.pdelab.geometry import world_dimension
+from dune.perftool.pdelab.geometry import world_dimension, local_dimension
from dune.perftool.generation import (class_member,
domain,
function_mangler,
@@ -51,11 +51,13 @@ class BasisTabulationMatrix(BasisTabulationMatrixBase, ImmutableRecord):
slice_size=None,
slice_index=None,
):
+ assert(isinstance(basis_size, int))
if quadrature_size is None:
quadrature_size = quadrature_points_per_direction()
+ assert(qp == quadrature_size[0] for qp in quadrature_size)
+ quadrature_size = quadrature_size[0]
if slice_size is not None:
quadrature_size = ceildiv(quadrature_size, slice_size)
-
ImmutableRecord.__init__(self,
quadrature_size=quadrature_size,
basis_size=basis_size,
@@ -81,13 +83,14 @@ class BasisTabulationMatrix(BasisTabulationMatrixBase, ImmutableRecord):
return self.basis_size
def pymbolic(self, indices):
- name = "{}{}Theta{}{}_qp{}_dof{}".format("face{}_".format(self.face) if self.face is not None else "",
- "d" if self.derivative else "",
- "T" if self.transpose else "",
- "_slice{}".format(self.slice_index) if self.slice_size is not None else "",
- self.quadrature_size,
- self.basis_size,
- )
+ name = "{}{}Theta{}{}_qp{}_dof{}" \
+ .format("face{}_".format(self.face) if self.face is not None else "",
+ "d" if self.derivative else "",
+ "T" if self.transpose else "",
+ "_slice{}".format(self.slice_index) if self.slice_size is not None else "",
+ self.quadrature_size,
+ self.basis_size,
+ )
define_theta(name, self)
return prim.Subscript(prim.Variable(name), indices)
@@ -199,51 +202,68 @@ class BasisTabulationMatrixArray(BasisTabulationMatrixBase):
def quadrature_points_per_direction():
- # Quadrature order
+ # Quadrature order per direction
q = quadrature_order()
+ if isinstance(q, int):
+ q = (q,) * world_dimension()
# Quadrature points in per direction
- nb_qp = q // 2 + 1
+ nb_qp = tuple(order // 2 + 1 for order in q)
return nb_qp
+def local_quadrature_points_per_direction():
+ """On a volume this simply gives the number of quadrature points per
+ direction. On a face it only returns the worlddim - 1 number of
+ quadrature points belonging to this face. (Delete normal 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')
+ qps_per_dir = qps_per_dir[:facedir] + qps_per_dir[facedir + 1:]
+ return qps_per_dir
+
+
def polynomial_degree():
form = get_global_context_value("formdata").preprocessed_form
- return form.coefficients()[0].ufl_element()._degree
+ degree = form.coefficients()[0].ufl_element().degree()
+ if isinstance(degree, int):
+ degree = (degree,) * world_dimension()
+ return degree
def basis_functions_per_direction():
- return polynomial_degree() + 1
+ return tuple(degree + 1 for degree in polynomial_degree())
-def define_oned_quadrature_weights(name):
+def define_oned_quadrature_weights(name, bound):
loopy_class_member(name,
dtype=np.float64,
classtag="operator",
- shape=(quadrature_points_per_direction(),),
+ shape=(bound,),
)
-def name_oned_quadrature_weights():
- num = quadrature_points_per_direction()
- name = "qw_num{}".format(num)
- define_oned_quadrature_weights(name)
+def name_oned_quadrature_weights(bound):
+ name = "qw_num{}".format(bound)
+ define_oned_quadrature_weights(name, bound)
return name
-def define_oned_quadrature_points(name):
+def define_oned_quadrature_points(name, bound):
loopy_class_member(name,
dtype=np.float64,
classtag="operator",
- shape=(quadrature_points_per_direction(),),
+ shape=(bound,),
)
-def name_oned_quadrature_points():
- num = quadrature_points_per_direction()
- name = "qp_num{}".format(num)
- define_oned_quadrature_points(name)
+def name_oned_quadrature_points(bound):
+ name = "qp_num{}".format(bound)
+ define_oned_quadrature_points(name, bound)
return name
@@ -284,12 +304,12 @@ def name_polynomials(degree):
@preamble(kernel="operator")
-def sort_quadrature_points_weights(qp, qw):
+def sort_quadrature_points_weights(qp, qw, bound):
range_field = lop_template_range_field()
domain_field = name_domain_field()
- number_qp = quadrature_points_per_direction()
include_file("dune/perftool/sumfact/onedquadrature.hh", filetag="operatorfile")
- return "onedQuadraturePointsWeights<{}, {}, {}>({}, {});".format(range_field, domain_field, number_qp, qp, qw)
+ return "onedQuadraturePointsWeights<{}, {}, {}>({}, {});" \
+ .format(range_field, domain_field, bound, qp, qw)
@iname(kernel="operator")
@@ -325,9 +345,12 @@ def polynomial_lookup_mangler(target, func, dtypes):
def define_theta(name, tabmat, additional_indices=(), width=None):
assert isinstance(tabmat, BasisTabulationMatrix)
- qp = name_oned_quadrature_points()
- qw = name_oned_quadrature_weights()
- sort_quadrature_points_weights(qp, qw)
+ bound = tabmat.quadrature_size
+ if tabmat.slice_size is not None:
+ bound *= tabmat.slice_size
+ qp = name_oned_quadrature_points(bound)
+ qw = name_oned_quadrature_weights(bound)
+ sort_quadrature_points_weights(qp, qw, bound)
degree = tabmat.basis_size - 1
polynomials = name_polynomials(degree)
@@ -370,17 +393,18 @@ def define_theta(name, tabmat, additional_indices=(), width=None):
def construct_basis_matrix_sequence(transpose=False, derivative=None, facedir=None, facemod=None, basis_size=None):
dim = world_dimension()
result = [None] * dim
+ quadrature_size = quadrature_points_per_direction()
+ assert (basis_size is not None)
+ if facedir is not None:
+ quadrature_size = quadrature_size[:facedir] + (1,) + quadrature_size[facedir:]
for i in range(dim):
- quadrature_size = quadrature_points_per_direction()
-
onface = None
if facedir == i:
- quadrature_size = 1
onface = facemod
- result[i] = BasisTabulationMatrix(quadrature_size=quadrature_size,
- basis_size=basis_size,
+ result[i] = BasisTabulationMatrix(quadrature_size=quadrature_size[i],
+ basis_size=basis_size[i],
transpose=transpose,
derivative=derivative == i,
face=onface)
diff --git a/python/dune/perftool/sumfact/vectorization.py b/python/dune/perftool/sumfact/vectorization.py
index 790fa59a1f052078800785cf10ed35fbf988da03..d5e855291816e27cc6627591a54bfa32b43bc1a6 100644
--- a/python/dune/perftool/sumfact/vectorization.py
+++ b/python/dune/perftool/sumfact/vectorization.py
@@ -9,7 +9,6 @@ from dune.perftool.generation import (generator_factory,
from dune.perftool.pdelab.restriction import (Restriction,
restricted_name,
)
-from dune.perftool.sumfact.tabulation import BasisTabulationMatrixArray
from dune.perftool.error import PerftoolError
from dune.perftool.options import get_option
diff --git a/python/dune/perftool/ufl/visitor.py b/python/dune/perftool/ufl/visitor.py
index b4283031662228c785235c4e5f69fd21f032c0e9..7e6187e0bf30b79a6d770a80c44492ee2d6835be 100644
--- a/python/dune/perftool/ufl/visitor.py
+++ b/python/dune/perftool/ufl/visitor.py
@@ -23,7 +23,9 @@ from pymbolic.primitives import (Call,
from ufl.algorithms import MultiFunction
from ufl.checks import is_cellwise_constant
from ufl import (VectorElement,
+ MixedElement,
TensorElement,
+ TensorProductElement,
)
from ufl.classes import (FixedIndex,
IndexSum,
@@ -98,7 +100,7 @@ class UFL2LoopyVisitor(ModifiedTerminalTracker):
leaf_element = o.ufl_element()
# Select the correct leaf element in the case of this being a mixed finite element
- if o.ufl_element().num_sub_elements() > 0:
+ if isinstance(o.ufl_element(), MixedElement):
index = self.indices[0]
assert isinstance(index, int)
self.indices = self.indices[1:]
@@ -127,7 +129,7 @@ class UFL2LoopyVisitor(ModifiedTerminalTracker):
self.interface.initialize_function_spaces(o, self)
index = None
- if o.ufl_element().num_sub_elements() > 0:
+ if isinstance(o.ufl_element(), MixedElement):
index = self.indices[0]
assert isinstance(index, int)
self.indices = self.indices[1:]
diff --git a/test/sumfact/poisson/poisson_2d.mini b/test/sumfact/poisson/poisson_2d.mini
index ea68946824a53387f3f6ce6f19f0349ef309af05..129c2ff373937cd9fecd86dcdf38b78723fee339 100644
--- a/test/sumfact/poisson/poisson_2d.mini
+++ b/test/sumfact/poisson/poisson_2d.mini
@@ -18,6 +18,7 @@ numerical_jacobian = 1, 0 | expand num
compare_l2errorsquared = 4e-5, 4e-9 | expand deg
sumfact = 1
vectorize_grads = 1, 0 | expand grad
+quadrature_order = 2, 4
[formcompiler.ufl_variants]
degree = 1, 2 | expand deg
diff --git a/test/sumfact/poisson/poisson_2d.ufl b/test/sumfact/poisson/poisson_2d.ufl
index 0d494b70453fa4cc15bae46080540c39d550c520..d2c78a8d2a1a928ae1f941db4e1a337c8d308bbb 100644
--- a/test/sumfact/poisson/poisson_2d.ufl
+++ b/test/sumfact/poisson/poisson_2d.ufl
@@ -5,7 +5,10 @@ c = (0.5-x[0])**2 + (0.5-x[1])**2
g = exp(-1.*c)
f = 2*(2.-2*c)*g
-V = FiniteElement("CG", cell, degree)
+V_0 = FiniteElement("CG", interval, degree)
+V_1 = FiniteElement("CG", interval, degree)
+V = TensorProductElement(V_0, V_1, cell=cell)
+
u = TrialFunction(V)
v = TestFunction(V)