diff --git a/python/dune/perftool/pdelab/quadrature.py b/python/dune/perftool/pdelab/quadrature.py index 91a25b3318a047ccfb85a2be19dd5311d80cad1f..b87c79dcb5874710fd600a3d5f36f9f4f3dce2e3 100644 --- a/python/dune/perftool/pdelab/quadrature.py +++ b/python/dune/perftool/pdelab/quadrature.py @@ -1,6 +1,9 @@ +from ufl.algorithms import estimate_total_polynomial_degree + from dune.perftool import Restriction from dune.perftool.generation import (cached, domain, + get_global_context_value, iname, include_file, instruction, @@ -66,11 +69,36 @@ def name_quadrature_weight(): 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(): - # TODO how to determine integration order in a generic fashion here? - # Or maybe just read it from the ini file? Or define some mapping from the - # finite elements? - return "2" + # TODO: It should be possible to set the quadrature order as a + # user. A global option in the ini file is not really what we want + # since it would not be possible to take different quadrature + # 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(): diff --git a/python/dune/perftool/sumfact/amatrix.py b/python/dune/perftool/sumfact/amatrix.py index f1b8a6209591b08b543b3a906dec037bbecf1da4..8755a56abe2b8a77efd7840c33b13337c5a102eb 100644 --- a/python/dune/perftool/sumfact/amatrix.py +++ b/python/dune/perftool/sumfact/amatrix.py @@ -10,6 +10,7 @@ from dune.perftool.generation import (class_member, from dune.perftool.pdelab.localoperator import (name_domain_field, lop_template_range_field, ) +from dune.perftool.pdelab.quadrature import estimate_quadrature_order @class_member("operator") @@ -29,7 +30,7 @@ def quadrature_points_per_direction(): # as soon as we have a generic implementation # Quadrature order - q = 2 + q = estimate_quadrature_order() # Quadrature points in per direction nb_qp = q // 2 + 1