Skip to content
Snippets Groups Projects
Commit e0e60cb8 authored by René Heß's avatar René Heß
Browse files

Use ufl polynomial degree estimation

Note: needs ufl patch
parent b3d02d24
No related branches found
No related tags found
No related merge requests found
from ufl.algorithms import estimate_total_polynomial_degree
from dune.perftool import Restriction from dune.perftool import Restriction
from dune.perftool.generation import (cached, from dune.perftool.generation import (cached,
domain, domain,
get_global_context_value,
iname, iname,
include_file, include_file,
instruction, instruction,
...@@ -66,11 +69,36 @@ def name_quadrature_weight(): ...@@ -66,11 +69,36 @@ def name_quadrature_weight():
return name return name
def estimate_quadrature_order():
"""Estimate quadrature order using polynomial degree estimation from UFL"""
# According to UFL documentation estimate_total_polynomial_degree
# should only be called on preprocessed forms.
form = get_global_context_value("formdata").preprocessed_form
# 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
for i in integrals:
polynomial_degree = max(polynomial_degree, estimate_total_polynomial_degree(i))
# Note: In PDELab quadrature order m means that integration of
# polynomials of degree m is exact.
return polynomial_degree
def name_order(): def name_order():
# TODO how to determine integration order in a generic fashion here? # TODO: It should be possible to set the quadrature order as a
# Or maybe just read it from the ini file? Or define some mapping from the # user. A global option in the ini file is not really what we want
# finite elements? # since it would not be possible to take different quadrature
return "2" # order for different integral types. Maybe the following would
# work:
#
# - One quadrature order that can be used everywhere
# - [formcompiler.quadrature]: Give orders for specific integral
# types.
quadrature_order = estimate_quadrature_order()
return str(quadrature_order)
def quadrature_loop_statement(): def quadrature_loop_statement():
......
...@@ -10,6 +10,7 @@ from dune.perftool.generation import (class_member, ...@@ -10,6 +10,7 @@ from dune.perftool.generation import (class_member,
from dune.perftool.pdelab.localoperator import (name_domain_field, from dune.perftool.pdelab.localoperator import (name_domain_field,
lop_template_range_field, lop_template_range_field,
) )
from dune.perftool.pdelab.quadrature import estimate_quadrature_order
@class_member("operator") @class_member("operator")
...@@ -29,7 +30,7 @@ def quadrature_points_per_direction(): ...@@ -29,7 +30,7 @@ def quadrature_points_per_direction():
# as soon as we have a generic implementation # as soon as we have a generic implementation
# Quadrature order # Quadrature order
q = 2 q = estimate_quadrature_order()
# Quadrature points in per direction # Quadrature points in per direction
nb_qp = q // 2 + 1 nb_qp = q // 2 + 1
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment