17557791029

Committed 08 Sep 2025 04:38PM UTC coverage: 76.363% (+0.4%) from 75.917%

Build # 17557791029

Build Type

Pull #401

github

Committed by

web-flow

Commit Message

Merge branch 'main' into schnellerhase/remove-type-system

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

Run Details

495 of 534 new or added lines in 42 files covered. (92.7%)

6 existing lines in 2 files now uncovered.

9133 of 11960 relevant lines covered (76.36%)

0.76 hits per line

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

86.81

/ufl/algorithms/apply_restrictions.py

"""Apply restrictions.

This module contains the apply_restrictions algorithm which propagates
restrictions in a form towards the terminals.
"""

# Copyright (C) 2008-2016 Martin Sandve Alnæs
#
# This file is part of UFL (https://www.fenicsproject.org)
#
# SPDX-License-Identifier:    LGPL-3.0-or-later

from ufl.algorithms.map_integrands import map_integrand_dags
from ufl.classes import Restricted
from ufl.corealg.map_dag import map_expr_dag
from ufl.corealg.multifunction import MultiFunction
from ufl.domain import extract_unique_domain
from ufl.measure import integral_type_to_measure_name
from ufl.sobolevspace import H1


class RestrictionPropagator(MultiFunction):
    """Restriction propagator."""

    def __init__(self, side=None, apply_default=True):
        """Initialise."""
        MultiFunction.__init__(self)
        self.current_restriction = side
        self.default_restriction = "+"
        self.apply_default = apply_default
        # Caches for propagating the restriction with map_expr_dag
        self.vcaches = {"+": {}, "-": {}}
        self.rcaches = {"+": {}, "-": {}}
        if self.current_restriction is None:
            self._rp = {
                "+": RestrictionPropagator(side="+", apply_default=apply_default),
                "-": RestrictionPropagator(side="-", apply_default=apply_default),
            }

    def restricted(self, o):
        """When hitting a restricted quantity, visit child with a separate restriction algorithm."""
        # Assure that we have only two levels here, inside or outside
        # the Restricted node
        if self.current_restriction is not None:
            raise ValueError("Cannot restrict an expression twice.")
        # Configure a propagator for this side and apply to subtree
        side = o.side()
        return map_expr_dag(
            self._rp[side], o.ufl_operands[0], vcache=self.vcaches[side], rcache=self.rcaches[side]
        )

    # --- Reusable rules

    def _ignore_restriction(self, o):
        """Ignore current restriction.

        Quantity is independent of side also from a computational point
        of view.
        """
        return o

    def _require_restriction(self, o):
        """Restrict a discontinuous quantity to current side, require a side to be set."""
        if self.current_restriction is None:
            raise ValueError(f"Discontinuous type {type(o).__name__} must be restricted.")
        return o(self.current_restriction)

    def _default_restricted(self, o):
        """Restrict a continuous quantity to default side if no current restriction is set."""
        r = self.current_restriction
        if r is not None:
            return o(r)
        if self.apply_default:
            return o(self.default_restriction)
        else:
            return o

    def _opposite(self, o):
        """Restrict a quantity to default side.

        If the current restriction is different swap the sign, require a side to be set.
        """
        if self.current_restriction is None:
            raise ValueError(f"Discontinuous type {type(o).__name__} must be restricted.")
        elif self.current_restriction == self.default_restriction:
            return o(self.default_restriction)
        else:
            return -o(self.default_restriction)

    def _missing_rule(self, o):
        """Raise an error."""
        raise ValueError(f"Missing rule for {type(o).__name__}")

    # --- Rules for operators

    # Default: Operators should reconstruct only if subtrees are not touched
    operator = MultiFunction.reuse_if_untouched

    # Assuming apply_derivatives has been called,
    # propagating Grad inside the Restricted nodes.
    # Considering all grads to be discontinuous, may
    # want something else for facet functions in future.
    grad = _require_restriction

    def variable(self, o, op, label):
        """Strip variable."""
        return op

    def reference_value(self, o):
        """Reference value of something follows same restriction rule as the underlying object."""
        (f,) = o.ufl_operands
        assert f._ufl_is_terminal_
        g = self(f)
        if isinstance(g, Restricted):
            side = g.side()
            return o(side)
        else:
            return o

    # --- Rules for terminals

    # Require handlers to be specified for all terminals
    terminal = _missing_rule

    multi_index = _ignore_restriction
    label = _ignore_restriction

    # Default: Literals should ignore restriction
    constant_value = _ignore_restriction
    constant = _ignore_restriction

    # Even arguments with continuous elements such as Lagrange must be
    # restricted to associate with the right part of the element
    # matrix
    argument = _require_restriction

    # Defaults for geometric quantities
    geometric_cell_quantity = _require_restriction
    geometric_facet_quantity = _require_restriction

    # Only a few geometric quantities are independent on the restriction:
    facet_coordinate = _ignore_restriction
    quadrature_weight = _ignore_restriction

    # Assuming homogeoneous mesh
    reference_cell_volume = _ignore_restriction
    reference_facet_volume = _ignore_restriction

    # These are the same from either side but to compute them
    # cell (or facet) data from one side must be selected:
    spatial_coordinate = _default_restricted
    # Depends on cell only to get to the facet:
    facet_jacobian = _default_restricted
    facet_jacobian_determinant = _default_restricted
    facet_jacobian_inverse = _default_restricted
    facet_area = _default_restricted
    min_facet_edge_length = _default_restricted
    max_facet_edge_length = _default_restricted
    facet_origin = _default_restricted  # FIXME: Is this valid for quads?

    def coefficient(self, o):
        """Restrict a coefficient.

        Allow coefficients to be unrestricted (apply default if so) if
        the values are fully continuous across the facet.
        """
        if o.ufl_element() in H1:
            # If the coefficient _value_ is _fully_ continuous
            # It must still be computed from one of the sides, we just don't care which
            return self._default_restricted(o)
        else:
            return self._require_restriction(o)

    def facet_normal(self, o):
        """Restrict a facet_normal."""
        D = extract_unique_domain(o)
        e = D.ufl_coordinate_element()
        gd = D.geometric_dimension()
        td = D.topological_dimension()

        if e.embedded_superdegree <= 1 and e in H1 and gd == td:
            # For meshes with a continuous linear non-manifold
            # coordinate field, the facet normal from side - points in
            # the opposite direction of the one from side +.  We must
            # still require a side to be chosen by the user but
            # rewrite n- -> n+.  This is an optimization, possibly
            # premature, however it's more difficult to do at a later
            # stage.
            return self._opposite(o)
        else:
            # For other meshes, we require a side to be
            # chosen by the user and respect that
            return self._require_restriction(o)


def apply_restrictions(expression, apply_default=True):
    """Propagate restriction nodes to wrap differential terminals directly."""
    integral_types = [
        k for k in integral_type_to_measure_name.keys() if k.startswith("interior_facet")
    ]
    rules = RestrictionPropagator(apply_default=apply_default)
    return map_integrand_dags(rules, expression, only_integral_type=integral_types)

1	"""Apply restrictions.
2
3	This module contains the apply_restrictions algorithm which propagates
4	restrictions in a form towards the terminals.
5	"""
6
7	# Copyright (C) 2008-2016 Martin Sandve Alnæs
8	#
9	# This file is part of UFL (https://www.fenicsproject.org)
10	#
11	# SPDX-License-Identifier: LGPL-3.0-or-later
12
13	from ufl.algorithms.map_integrands import map_integrand_dags	1✔
14	from ufl.classes import Restricted	1✔
15	from ufl.corealg.map_dag import map_expr_dag	1✔
16	from ufl.corealg.multifunction import MultiFunction	1✔
17	from ufl.domain import extract_unique_domain	1✔
18	from ufl.measure import integral_type_to_measure_name	1✔
19	from ufl.sobolevspace import H1	1✔
20
21
22	class RestrictionPropagator(MultiFunction):	1✔
23	"""Restriction propagator."""
24
25	def __init__(self, side=None, apply_default=True):	1✔
26	"""Initialise."""
27	MultiFunction.__init__(self)	1✔
28	self.current_restriction = side	1✔
29	self.default_restriction = "+"	1✔
30	self.apply_default = apply_default	1✔
31	# Caches for propagating the restriction with map_expr_dag
32	self.vcaches = {"+": {}, "-": {}}	1✔
33	self.rcaches = {"+": {}, "-": {}}	1✔
34	if self.current_restriction is None:	1✔
35	self._rp = {	1✔
36	"+": RestrictionPropagator(side="+", apply_default=apply_default),
37	"-": RestrictionPropagator(side="-", apply_default=apply_default),
38	}
39
40	def restricted(self, o):	1✔
41	"""When hitting a restricted quantity, visit child with a separate restriction algorithm."""
42	# Assure that we have only two levels here, inside or outside
43	# the Restricted node
44	if self.current_restriction is not None:	1✔
45	raise ValueError("Cannot restrict an expression twice.")	×
46	# Configure a propagator for this side and apply to subtree
47	side = o.side()	1✔
48	return map_expr_dag(	1✔
49	self._rp[side], o.ufl_operands[0], vcache=self.vcaches[side], rcache=self.rcaches[side]
50	)
51
52	# --- Reusable rules
53
54	def _ignore_restriction(self, o):	1✔
55	"""Ignore current restriction.
56
57	Quantity is independent of side also from a computational point
58	of view.
59	"""
60	return o	1✔
61
62	def _require_restriction(self, o):	1✔
63	"""Restrict a discontinuous quantity to current side, require a side to be set."""
64	if self.current_restriction is None:	1✔
65	raise ValueError(f"Discontinuous type {type(o).__name__} must be restricted.")	1✔
66	return o(self.current_restriction)	1✔
67
68	def _default_restricted(self, o):	1✔
69	"""Restrict a continuous quantity to default side if no current restriction is set."""
70	r = self.current_restriction	1✔
71	if r is not None:	1✔
72	return o(r)	1✔
73	if self.apply_default:	1✔
74	return o(self.default_restriction)	1✔
75	else:
76	return o	×
77
78	def _opposite(self, o):	1✔
79	"""Restrict a quantity to default side.
80
81	If the current restriction is different swap the sign, require a side to be set.
82	"""
83	if self.current_restriction is None:	1✔
84	raise ValueError(f"Discontinuous type {type(o).__name__} must be restricted.")	1✔
85	elif self.current_restriction == self.default_restriction:	1✔
86	return o(self.default_restriction)	1✔
87	else:
88	return -o(self.default_restriction)	1✔
89
90	def _missing_rule(self, o):	1✔
91	"""Raise an error."""
NEW 92	raise ValueError(f"Missing rule for {type(o).__name__}")	×
93
94	# --- Rules for operators
95
96	# Default: Operators should reconstruct only if subtrees are not touched
97	operator = MultiFunction.reuse_if_untouched	1✔
98
99	# Assuming apply_derivatives has been called,
100	# propagating Grad inside the Restricted nodes.
101	# Considering all grads to be discontinuous, may
102	# want something else for facet functions in future.
103	grad = _require_restriction	1✔
104
105	def variable(self, o, op, label):	1✔
106	"""Strip variable."""
107	return op	×
108
109	def reference_value(self, o):	1✔
110	"""Reference value of something follows same restriction rule as the underlying object."""
111	(f,) = o.ufl_operands	×
112	assert f._ufl_is_terminal_	×
113	g = self(f)	×
114	if isinstance(g, Restricted):	×
115	side = g.side()	×
116	return o(side)	×
117	else:
118	return o	×
119
120	# --- Rules for terminals
121
122	# Require handlers to be specified for all terminals
123	terminal = _missing_rule	1✔
124
125	multi_index = _ignore_restriction	1✔
126	label = _ignore_restriction	1✔
127
128	# Default: Literals should ignore restriction
129	constant_value = _ignore_restriction	1✔
130	constant = _ignore_restriction	1✔
131
132	# Even arguments with continuous elements such as Lagrange must be
133	# restricted to associate with the right part of the element
134	# matrix
135	argument = _require_restriction	1✔
136
137	# Defaults for geometric quantities
138	geometric_cell_quantity = _require_restriction	1✔
139	geometric_facet_quantity = _require_restriction	1✔
140
141	# Only a few geometric quantities are independent on the restriction:
142	facet_coordinate = _ignore_restriction	1✔
143	quadrature_weight = _ignore_restriction	1✔
144
145	# Assuming homogeoneous mesh
146	reference_cell_volume = _ignore_restriction	1✔
147	reference_facet_volume = _ignore_restriction	1✔
148
149	# These are the same from either side but to compute them
150	# cell (or facet) data from one side must be selected:
151	spatial_coordinate = _default_restricted	1✔
152	# Depends on cell only to get to the facet:
153	facet_jacobian = _default_restricted	1✔
154	facet_jacobian_determinant = _default_restricted	1✔
155	facet_jacobian_inverse = _default_restricted	1✔
156	facet_area = _default_restricted	1✔
157	min_facet_edge_length = _default_restricted	1✔
158	max_facet_edge_length = _default_restricted	1✔
159	facet_origin = _default_restricted # FIXME: Is this valid for quads?	1✔
160
161	def coefficient(self, o):	1✔
162	"""Restrict a coefficient.
163
164	Allow coefficients to be unrestricted (apply default if so) if
165	the values are fully continuous across the facet.
166	"""
167	if o.ufl_element() in H1:	1✔
168	# If the coefficient _value_ is _fully_ continuous
169	# It must still be computed from one of the sides, we just don't care which
170	return self._default_restricted(o)	1✔
171	else:
172	return self._require_restriction(o)	1✔
173
174	def facet_normal(self, o):	1✔
175	"""Restrict a facet_normal."""
176	D = extract_unique_domain(o)	1✔
177	e = D.ufl_coordinate_element()	1✔
178	gd = D.geometric_dimension()	1✔
179	td = D.topological_dimension()	1✔
180
181	if e.embedded_superdegree <= 1 and e in H1 and gd == td:	1✔
182	# For meshes with a continuous linear non-manifold
183	# coordinate field, the facet normal from side - points in
184	# the opposite direction of the one from side +. We must
185	# still require a side to be chosen by the user but
186	# rewrite n- -> n+. This is an optimization, possibly
187	# premature, however it's more difficult to do at a later
188	# stage.
189	return self._opposite(o)	1✔
190	else:
191	# For other meshes, we require a side to be
192	# chosen by the user and respect that
193	return self._require_restriction(o)	×
194
195
196	def apply_restrictions(expression, apply_default=True):	1✔
197	"""Propagate restriction nodes to wrap differential terminals directly."""
198	integral_types = [	1✔
199	k for k in integral_type_to_measure_name.keys() if k.startswith("interior_facet")
200	]
201	rules = RestrictionPropagator(apply_default=apply_default)	1✔
202	return map_integrand_dags(rules, expression, only_integral_type=integral_types)	1✔

FEniCS / ufl / 17557791029

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