18629405325

Committed 19 Oct 2025 10:56AM UTC coverage: 77.06% (+0.4%) from 76.622%

Build # 18629405325

Build Type

Pull #401

github

Committed by

schnellerhase

Commit Message

Ruff

Pull Request Pull Request #401: Removal of custom type system

Run Details

494 of 533 new or added lines in 41 files covered. (92.68%)

6 existing lines in 2 files now uncovered.

9325 of 12101 relevant lines covered (77.06%)

0.77 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

79.5

/ufl/algorithms/estimate_degrees.py

"""Algorithms for estimating polynomial degrees of expressions."""

# Copyright (C) 2008-2016 Martin Sandve Alnæs and Anders Logg
#
# This file is part of UFL (https://www.fenicsproject.org)
#
# SPDX-License-Identifier:    LGPL-3.0-or-later
#
# Modified by Anders Logg, 2009-2010
# Modified by Jan Blechta, 2012

import warnings

from ufl.checks import is_cellwise_constant
from ufl.constantvalue import IntValue
from ufl.corealg.map_dag import map_expr_dags
from ufl.corealg.multifunction import MultiFunction
from ufl.domain import extract_domains, extract_unique_domain
from ufl.form import Form
from ufl.integral import Integral


class SumDegreeEstimator(MultiFunction):
    """Sum degree estimator.

    This algorithm is exact for a few operators and heuristic for many.
    """

    def __init__(self, default_degree, element_replace_map):
        """Initialise."""
        MultiFunction.__init__(self)
        self.default_degree = default_degree
        self.element_replace_map = element_replace_map

    def constant_value(self, v):
        """Apply to constant_value.

        Constant values are constant.
        """
        return 0

    def constant(self, v):
        """Apply to constant."""
        return 0

    def geometric_quantity(self, v):
        """Apply to geometric_quantity.

        Some geometric quantities are cellwise constant. Others are
        nonpolynomial and thus hard to estimate.
        """
        if is_cellwise_constant(v):
            return 0
        else:
            # As a heuristic, just returning domain degree to bump up degree somewhat
            return extract_unique_domain(v).ufl_coordinate_element().embedded_superdegree

    def spatial_coordinate(self, v):
        """Apply to spatial_coordinate.

        A coordinate provides additional degrees depending on coordinate field of domain.
        """
        return extract_unique_domain(v).ufl_coordinate_element().embedded_superdegree

    def cell_coordinate(self, v):
        """Apply to cell_coordinate.

        A coordinate provides one additional degree.
        """
        return 1

    def argument(self, v):
        """Apply to argument.

        A form argument provides a degree depending on the element,
        or the default degree if the element has no degree.
        """
        return (
            v.ufl_element().embedded_superdegree
        )  # FIXME: Use component to improve accuracy for mixed elements

    def coefficient(self, v):
        """Apply to coefficient.

        A form argument provides a degree depending on the element,
        or the default degree if the element has no degree.
        """
        e = v.ufl_element()
        e = self.element_replace_map.get(e, e)
        d = e.embedded_superdegree  # FIXME: Use component to improve accuracy for mixed elements
        if d is None:
            d = self.default_degree
        return d

    def _reduce_degree(self, v, f):
        """Reduce the estimated degree by one.

        This is used when derivatives are taken. It does not reduce the degree when
        TensorProduct elements or quadrilateral elements are involved.
        """
        # Can have multiple domains of the same cell type.
        (cell,) = set(d.ufl_cell() for d in extract_domains(v))
        if isinstance(f, int) and cell.cellname() not in ["quadrilateral", "hexahedron"]:
            return max(f - 1, 0)
        else:
            return f

    def _add_degrees(self, v, *ops):
        """Apply to _add_degrees."""
        if any(isinstance(o, tuple) for o in ops):
            # we can add a slight hack here to handle things
            # like adding 0 to (3, 3) [by expanding
            # 0 to (0, 0) when making tempops]
            tempops = [foo if isinstance(foo, tuple) else (foo, foo) for foo in ops]
            return tuple(map(sum, zip(*tempops)))
        else:
            return sum(ops)

    def _max_degrees(self, v, *ops):
        """Apply to _max_degrees."""
        if any(isinstance(o, tuple) for o in ops):
            tempops = [foo if isinstance(foo, tuple) else (foo, foo) for foo in ops]
            return tuple(map(max, zip(*tempops)))
        else:
            return max(ops + (0,))

    def _not_handled(self, v, *args):
        """Apply to _not_handled."""
        raise ValueError(f"Missing degree handler for type {type(v).__name__}")

    def expr(self, v, *ops):
        """Apply to expr.

        For most operators we take the max degree of its operands.
        """
        warnings.warn(f"Missing degree estimation handler for type {type(v).__name__}")
        return self._add_degrees(v, *ops)

    # Utility types with no degree concept
    def multi_index(self, v):
        """Apply to multi_index."""
        return None

    def label(self, v):
        """Apply to label."""
        return None

    # Fall-through, indexing and similar types
    def reference_value(self, rv, f):
        """Apply to reference_value."""
        return f

    def variable(self, v, e, a):
        """Apply to variable."""
        return e

    def transposed(self, v, A):
        """Apply to transposed."""
        return A

    def index_sum(self, v, A, ii):
        """Apply to index_sum."""
        return A

    def indexed(self, v, A, ii):
        """Apply to indexed."""
        return A

    def component_tensor(self, v, A, ii):
        """Apply to component_tensor."""
        return A

    list_tensor = _max_degrees

    def positive_restricted(self, v, a):
        """Apply to positive_restricted."""
        return a

    def negative_restricted(self, v, a):
        """Apply to negative_restricted."""
        return a

    def conj(self, v, a):
        """Apply to conj."""
        return a

    def real(self, v, a):
        """Apply to real."""
        return a

    def imag(self, v, a):
        """Apply to imag."""
        return a

    # A sum takes the max degree of its operands:
    sum = _max_degrees

    # TODO: Need a new algorithm which considers direction of
    # derivatives of form arguments A spatial derivative reduces the
    # degree with one
    grad = _reduce_degree
    reference_grad = _reduce_degree
    # Handling these types although they should not occur... please
    # apply preprocessing before using this algorithm:
    nabla_grad = _reduce_degree
    div = _reduce_degree
    reference_div = _reduce_degree
    nabla_div = _reduce_degree
    curl = _reduce_degree
    reference_curl = _reduce_degree

    def cell_avg(self, v, a):
        """Apply to cell_avg.

        Cell average of a function is always cellwise constant.
        """
        return 0

    def facet_avg(self, v, a):
        """Apply to facet_avg.

        Facet average of a function is always cellwise constant.
        """
        return 0

    # A product accumulates the degrees of its operands:
    product = _add_degrees
    # Handling these types although they should not occur... please
    # apply preprocessing before using this algorithm:
    inner = _add_degrees
    dot = _add_degrees
    outer = _add_degrees
    cross = _add_degrees

    # Explicitly not handling these types, please apply preprocessing
    # before using this algorithm:
    derivative = _not_handled  # base type
    compound_derivative = _not_handled  # base type
    compound_tensor_operator = _not_handled  # base class
    variable_derivative = _not_handled
    trace = _not_handled
    determinant = _not_handled
    cofactor = _not_handled
    inverse = _not_handled
    deviatoric = _not_handled
    skew = _not_handled
    sym = _not_handled

    def abs(self, v, a):
        """Apply to abs.

        This is a heuristic, correct if there is no.
        """
        if a == 0:
            return a
        else:
            return a

    def division(self, v, *ops):
        """Apply to division.

        Using the sum here is a heuristic. Consider e.g. (x+1)/(x-1).
        """
        return self._add_degrees(v, *ops)

    def power(self, v, a, b):
        """Apply to power.

        If b is a positive integer: degree(a**b) == degree(a)*b
        otherwise use the heuristic: degree(a**b) == degree(a) + 2.
        """
        _f, g = v.ufl_operands

        if isinstance(g, IntValue):
            gi = g.value()
            if gi >= 0:
                if isinstance(a, int):
                    return a * gi
                else:
                    return tuple(foo * gi for foo in a)

        # Something to a non-(positive integer) power, e.g. float,
        # negative integer, Coefficient, etc.
        return self._add_degrees(v, a, 2)

    def atan2(self, v, a, b):
        """Apply to atan2.

        Using the heuristic:
        degree(atan2(const,const)) == 0
        degree(atan2(a,b)) == max(degree(a),degree(b))+2
        which can be wildly inaccurate but at least gives a somewhat
        high integration degree.
        """
        if a or b:
            return self._add_degrees(v, self._max_degrees(v, a, b), 2)
        else:
            return self._max_degrees(v, a, b)

    def math_function(self, v, a):
        """Apply to math_function.

        Using the heuristic:
        degree(sin(const)) == 0
        degree(sin(a)) == degree(a)+2
        which can be wildly inaccurate but at least gives a somewhat
        high integration degree.
        """
        if a:
            return self._add_degrees(v, a, 2)
        else:
            return a

    def bessel_function(self, v, nu, x):
        """Apply to bessel_function.

        Using the heuristic
        degree(bessel_*(const)) == 0
        degree(bessel_*(x)) == degree(x)+2
        which can be wildly inaccurate but at least gives a somewhat
        high integration degree.
        """
        if x:
            return self._add_degrees(v, x, 2)
        else:
            return x

    def condition(self, v, *args):
        """Apply to condition."""
        return None

    def conditional(self, v, c, t, f):
        """Apply to conditional.

        Degree of condition does not influence degree of values which conditional takes. So
        heuristicaly taking max of true degree and false degree. This will be exact in cells
        where condition takes single value. For improving accuracy of quadrature near
        condition transition surface quadrature order must be adjusted manually.
        """
        return self._max_degrees(v, t, f)

    def min_value(self, v, a, r):
        """Apply to min_value.

        Same as conditional.
        """
        return self._max_degrees(v, a, r)

    max_value = min_value

    def coordinate_derivative(self, v, integrand_degree, b, direction_degree, d):
        """Apply to coordinate_derivative.

        We use the heuristic that a shape derivative in direction V
        introduces terms V and grad(V) into the integrand. Hence we add the
        degree of the deformation to the estimate.
        """
        return self._add_degrees(v, integrand_degree, direction_degree)

    def expr_list(self, v, *o):
        """Apply to expr_list."""
        return self._max_degrees(v, *o)

    def expr_mapping(self, v, *o):
        """Apply to expr_mapping."""
        return self._max_degrees(v, *o)


def estimate_total_polynomial_degree(e, default_degree=1, element_replace_map={}):
    """Estimate total polynomial degree of integrand.

    NB: Although some compound types are supported here,
    some derivatives and compounds must be preprocessed
    prior to degree estimation. In generic code, this algorithm
    should only be applied after preprocessing.

    For coefficients defined on an element with unspecified degree
    (None), the degree is set to the given default degree.
    """
    de = SumDegreeEstimator(default_degree, element_replace_map)
    if isinstance(e, Form):
        if not e.integrals():
            raise ValueError("Form has no integrals.")
        degrees = map_expr_dags(de, [it.integrand() for it in e.integrals()])
    elif isinstance(e, Integral):
        degrees = map_expr_dags(de, [e.integrand()])
    else:
        degrees = map_expr_dags(de, [e])
    degree = max(degrees) if degrees else default_degree
    return degree

1	"""Algorithms for estimating polynomial degrees of expressions."""
2
3	# Copyright (C) 2008-2016 Martin Sandve Alnæs and Anders Logg
4	#
5	# This file is part of UFL (https://www.fenicsproject.org)
6	#
7	# SPDX-License-Identifier: LGPL-3.0-or-later
8	#
9	# Modified by Anders Logg, 2009-2010
10	# Modified by Jan Blechta, 2012
11
12	import warnings	1✔
13
14	from ufl.checks import is_cellwise_constant	1✔
15	from ufl.constantvalue import IntValue	1✔
16	from ufl.corealg.map_dag import map_expr_dags	1✔
17	from ufl.corealg.multifunction import MultiFunction	1✔
18	from ufl.domain import extract_domains, extract_unique_domain	1✔
19	from ufl.form import Form	1✔
20	from ufl.integral import Integral	1✔
21
22
23	class SumDegreeEstimator(MultiFunction):	1✔
24	"""Sum degree estimator.
25
26	This algorithm is exact for a few operators and heuristic for many.
27	"""
28
29	def __init__(self, default_degree, element_replace_map):	1✔
30	"""Initialise."""
31	MultiFunction.__init__(self)	1✔
32	self.default_degree = default_degree	1✔
33	self.element_replace_map = element_replace_map	1✔
34
35	def constant_value(self, v):	1✔
36	"""Apply to constant_value.
37
38	Constant values are constant.
39	"""
40	return 0	1✔
41
42	def constant(self, v):	1✔
43	"""Apply to constant."""
44	return 0	1✔
45
46	def geometric_quantity(self, v):	1✔
47	"""Apply to geometric_quantity.
48
49	Some geometric quantities are cellwise constant. Others are
50	nonpolynomial and thus hard to estimate.
51	"""
52	if is_cellwise_constant(v):	1✔
53	return 0	1✔
54	else:
55	# As a heuristic, just returning domain degree to bump up degree somewhat
56	return extract_unique_domain(v).ufl_coordinate_element().embedded_superdegree	1✔
57
58	def spatial_coordinate(self, v):	1✔
59	"""Apply to spatial_coordinate.
60
61	A coordinate provides additional degrees depending on coordinate field of domain.
62	"""
63	return extract_unique_domain(v).ufl_coordinate_element().embedded_superdegree	1✔
64
65	def cell_coordinate(self, v):	1✔
66	"""Apply to cell_coordinate.
67
68	A coordinate provides one additional degree.
69	"""
70	return 1	×
71
72	def argument(self, v):	1✔
73	"""Apply to argument.
74
75	A form argument provides a degree depending on the element,
76	or the default degree if the element has no degree.
77	"""
78	return (	1✔
79	v.ufl_element().embedded_superdegree
80	) # FIXME: Use component to improve accuracy for mixed elements
81
82	def coefficient(self, v):	1✔
83	"""Apply to coefficient.
84
85	A form argument provides a degree depending on the element,
86	or the default degree if the element has no degree.
87	"""
88	e = v.ufl_element()	1✔
89	e = self.element_replace_map.get(e, e)	1✔
90	d = e.embedded_superdegree # FIXME: Use component to improve accuracy for mixed elements	1✔
91	if d is None:	1✔
92	d = self.default_degree	×
93	return d	1✔
94
95	def _reduce_degree(self, v, f):	1✔
96	"""Reduce the estimated degree by one.
97
98	This is used when derivatives are taken. It does not reduce the degree when
99	TensorProduct elements or quadrilateral elements are involved.
100	"""
101	# Can have multiple domains of the same cell type.
102	(cell,) = set(d.ufl_cell() for d in extract_domains(v))	1✔
103	if isinstance(f, int) and cell.cellname() not in ["quadrilateral", "hexahedron"]:	1✔
104	return max(f - 1, 0)	1✔
105	else:
106	return f	×
107
108	def _add_degrees(self, v, *ops):	1✔
109	"""Apply to _add_degrees."""
110	if any(isinstance(o, tuple) for o in ops):	1✔
111	# we can add a slight hack here to handle things
112	# like adding 0 to (3, 3) [by expanding
113	# 0 to (0, 0) when making tempops]
114	tempops = [foo if isinstance(foo, tuple) else (foo, foo) for foo in ops]	×
115	return tuple(map(sum, zip(*tempops)))	×
116	else:
117	return sum(ops)	1✔
118
119	def _max_degrees(self, v, *ops):	1✔
120	"""Apply to _max_degrees."""
121	if any(isinstance(o, tuple) for o in ops):	1✔
122	tempops = [foo if isinstance(foo, tuple) else (foo, foo) for foo in ops]	×
123	return tuple(map(max, zip(*tempops)))	×
124	else:
125	return max(ops + (0,))	1✔
126
127	def _not_handled(self, v, *args):	1✔
128	"""Apply to _not_handled."""
NEW 129	raise ValueError(f"Missing degree handler for type {type(v).__name__}")	×
130
131	def expr(self, v, *ops):	1✔
132	"""Apply to expr.
133
134	For most operators we take the max degree of its operands.
135	"""
NEW 136	warnings.warn(f"Missing degree estimation handler for type {type(v).__name__}")	×
137	return self._add_degrees(v, *ops)	×
138
139	# Utility types with no degree concept
140	def multi_index(self, v):	1✔
141	"""Apply to multi_index."""
142	return None	1✔
143
144	def label(self, v):	1✔
145	"""Apply to label."""
146	return None	1✔
147
148	# Fall-through, indexing and similar types
149	def reference_value(self, rv, f):	1✔
150	"""Apply to reference_value."""
151	return f	×
152
153	def variable(self, v, e, a):	1✔
154	"""Apply to variable."""
155	return e	1✔
156
157	def transposed(self, v, A):	1✔
158	"""Apply to transposed."""
159	return A	×
160
161	def index_sum(self, v, A, ii):	1✔
162	"""Apply to index_sum."""
163	return A	1✔
164
165	def indexed(self, v, A, ii):	1✔
166	"""Apply to indexed."""
167	return A	1✔
168
169	def component_tensor(self, v, A, ii):	1✔
170	"""Apply to component_tensor."""
171	return A	1✔
172
173	list_tensor = _max_degrees	1✔
174
175	def positive_restricted(self, v, a):	1✔
176	"""Apply to positive_restricted."""
177	return a	1✔
178
179	def negative_restricted(self, v, a):	1✔
180	"""Apply to negative_restricted."""
181	return a	1✔
182
183	def conj(self, v, a):	1✔
184	"""Apply to conj."""
185	return a	1✔
186
187	def real(self, v, a):	1✔
188	"""Apply to real."""
189	return a	1✔
190
191	def imag(self, v, a):	1✔
192	"""Apply to imag."""
193	return a	1✔
194
195	# A sum takes the max degree of its operands:
196	sum = _max_degrees	1✔
197
198	# TODO: Need a new algorithm which considers direction of
199	# derivatives of form arguments A spatial derivative reduces the
200	# degree with one
201	grad = _reduce_degree	1✔
202	reference_grad = _reduce_degree	1✔
203	# Handling these types although they should not occur... please
204	# apply preprocessing before using this algorithm:
205	nabla_grad = _reduce_degree	1✔
206	div = _reduce_degree	1✔
207	reference_div = _reduce_degree	1✔
208	nabla_div = _reduce_degree	1✔
209	curl = _reduce_degree	1✔
210	reference_curl = _reduce_degree	1✔
211
212	def cell_avg(self, v, a):	1✔
213	"""Apply to cell_avg.
214
215	Cell average of a function is always cellwise constant.
216	"""
217	return 0	×
218
219	def facet_avg(self, v, a):	1✔
220	"""Apply to facet_avg.
221
222	Facet average of a function is always cellwise constant.
223	"""
224	return 0	×
225
226	# A product accumulates the degrees of its operands:
227	product = _add_degrees	1✔
228	# Handling these types although they should not occur... please
229	# apply preprocessing before using this algorithm:
230	inner = _add_degrees	1✔
231	dot = _add_degrees	1✔
232	outer = _add_degrees	1✔
233	cross = _add_degrees	1✔
234
235	# Explicitly not handling these types, please apply preprocessing
236	# before using this algorithm:
237	derivative = _not_handled # base type	1✔
238	compound_derivative = _not_handled # base type	1✔
239	compound_tensor_operator = _not_handled # base class	1✔
240	variable_derivative = _not_handled	1✔
241	trace = _not_handled	1✔
242	determinant = _not_handled	1✔
243	cofactor = _not_handled	1✔
244	inverse = _not_handled	1✔
245	deviatoric = _not_handled	1✔
246	skew = _not_handled	1✔
247	sym = _not_handled	1✔
248
249	def abs(self, v, a):	1✔
250	"""Apply to abs.
251
252	This is a heuristic, correct if there is no.
253	"""
254	if a == 0:	1✔
255	return a	1✔
256	else:
257	return a	×
258
259	def division(self, v, *ops):	1✔
260	"""Apply to division.
261
262	Using the sum here is a heuristic. Consider e.g. (x+1)/(x-1).
263	"""
264	return self._add_degrees(v, *ops)	1✔
265
266	def power(self, v, a, b):	1✔
267	"""Apply to power.
268
269	If b is a positive integer: degree(a*b) == degree(a)b
270	otherwise use the heuristic: degree(a**b) == degree(a) + 2.
271	"""
272	_f, g = v.ufl_operands	1✔
273
274	if isinstance(g, IntValue):	1✔
275	gi = g.value()	1✔
276	if gi >= 0:	1✔
277	if isinstance(a, int):	1✔
278	return a * gi	1✔
279	else:
280	return tuple(foo * gi for foo in a)	×
281
282	# Something to a non-(positive integer) power, e.g. float,
283	# negative integer, Coefficient, etc.
284	return self._add_degrees(v, a, 2)	×
285
286	def atan2(self, v, a, b):	1✔
287	"""Apply to atan2.
288
289	Using the heuristic:
290	degree(atan2(const,const)) == 0
291	degree(atan2(a,b)) == max(degree(a),degree(b))+2
292	which can be wildly inaccurate but at least gives a somewhat
293	high integration degree.
294	"""
295	if a or b:	×
296	return self._add_degrees(v, self._max_degrees(v, a, b), 2)	×
297	else:
298	return self._max_degrees(v, a, b)	×
299
300	def math_function(self, v, a):	1✔
301	"""Apply to math_function.
302
303	Using the heuristic:
304	degree(sin(const)) == 0
305	degree(sin(a)) == degree(a)+2
306	which can be wildly inaccurate but at least gives a somewhat
307	high integration degree.
308	"""
309	if a:	1✔
310	return self._add_degrees(v, a, 2)	1✔
311	else:
312	return a	1✔
313
314	def bessel_function(self, v, nu, x):	1✔
315	"""Apply to bessel_function.
316
317	Using the heuristic
318	degree(bessel_*(const)) == 0
319	degree(bessel_*(x)) == degree(x)+2
320	which can be wildly inaccurate but at least gives a somewhat
321	high integration degree.
322	"""
323	if x:	×
324	return self._add_degrees(v, x, 2)	×
325	else:
326	return x	×
327
328	def condition(self, v, *args):	1✔
329	"""Apply to condition."""
330	return None	×
331
332	def conditional(self, v, c, t, f):	1✔
333	"""Apply to conditional.
334
335	Degree of condition does not influence degree of values which conditional takes. So
336	heuristicaly taking max of true degree and false degree. This will be exact in cells
337	where condition takes single value. For improving accuracy of quadrature near
338	condition transition surface quadrature order must be adjusted manually.
339	"""
340	return self._max_degrees(v, t, f)	×
341
342	def min_value(self, v, a, r):	1✔
343	"""Apply to min_value.
344
345	Same as conditional.
346	"""
347	return self._max_degrees(v, a, r)	×
348
349	max_value = min_value	1✔
350
351	def coordinate_derivative(self, v, integrand_degree, b, direction_degree, d):	1✔
352	"""Apply to coordinate_derivative.
353
354	We use the heuristic that a shape derivative in direction V
355	introduces terms V and grad(V) into the integrand. Hence we add the
356	degree of the deformation to the estimate.
357	"""
358	return self._add_degrees(v, integrand_degree, direction_degree)	×
359
360	def expr_list(self, v, *o):	1✔
361	"""Apply to expr_list."""
362	return self._max_degrees(v, *o)	×
363
364	def expr_mapping(self, v, *o):	1✔
365	"""Apply to expr_mapping."""
366	return self._max_degrees(v, *o)	×
367
368
369	def estimate_total_polynomial_degree(e, default_degree=1, element_replace_map={}):	1✔
370	"""Estimate total polynomial degree of integrand.
371
372	NB: Although some compound types are supported here,
373	some derivatives and compounds must be preprocessed
374	prior to degree estimation. In generic code, this algorithm
375	should only be applied after preprocessing.
376
377	For coefficients defined on an element with unspecified degree
378	(None), the degree is set to the given default degree.
379	"""
380	de = SumDegreeEstimator(default_degree, element_replace_map)	1✔
381	if isinstance(e, Form):	1✔
382	if not e.integrals():	×
383	raise ValueError("Form has no integrals.")	×
384	degrees = map_expr_dags(de, [it.integrand() for it in e.integrals()])	×
385	elif isinstance(e, Integral):	1✔
386	degrees = map_expr_dags(de, [e.integrand()])	×
387	else:
388	degrees = map_expr_dags(de, [e])	1✔
389	degree = max(degrees) if degrees else default_degree	1✔
390	return degree	1✔

FEniCS / ufl / 18629405325

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous