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

59.74

/ufl/algorithms/change_to_reference.py

"""Algorithm for replacing gradients in an expression."""

# 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.apply_function_pullbacks import apply_function_pullbacks
from ufl.algorithms.apply_geometry_lowering import apply_geometry_lowering
from ufl.checks import is_cellwise_constant
from ufl.classes import Grad, JacobianInverse, ReferenceGrad, ReferenceValue, Restricted
from ufl.core.multiindex import indices
from ufl.corealg.map_dag import map_expr_dag
from ufl.corealg.multifunction import MultiFunction
from ufl.domain import extract_unique_domain
from ufl.tensors import as_tensor

"""
# Some notes:
# Below, let v_i mean physical coordinate of vertex i and V_i mean the
reference cell coordinate of the same vertex.


# Add a type for CellVertices? Note that vertices must be computed in
linear cell cases! triangle_vertices[i,j] = component j of vertex i,
following ufc numbering conventions


# DONE Add a type for CellEdgeLengths? Note that these are only easy to
define in the linear cell case!
triangle_edge_lengths    = [v1v2, v0v2, v0v1] # shape (3,)
tetrahedron_edge_lengths = [v0v1, v0v2, v0v3, v1v2, v1v3, v2v3] # shape (6,)


# DONE Here's how to compute edge lengths from the Jacobian:
J =[ [dx0/dX0, dx0/dX1],
     [dx1/dX0, dx1/dX1] ]
# First compute the edge vector, which is constant for each edge: the
vector from one vertex to the other
reference_edge_vector_0 = V2 - V1 # Example! Add a type ReferenceEdgeVectors?
# Then apply J to it and take the length of the resulting vector, this is generic for affine cells
edge_length_i = || dot(J, reference_edge_vector_i) ||

e2 = || J[:,0] . < 1, 0> || = || J[:,0] || = || dx/dX0 || = edge length of edge 2 (v0-v1)
e1 = || J[:,1] . < 0, 1> || = || J[:,1] || = || dx/dX1 || = edge length of edge 1 (v0-v2)
e0 = || J[:,:] . <-1, 1> || = || < J[0,1]-J[0,0], J[1,1]-J[1,0] > || = || dx/dX <-1,1> ||
    = edge length of edge 0 (v1-v2)

trev = triangle_reference_edge_vector
evec0 = J00 * trev[edge][0] + J01 * trev[edge][1]  =  J*trev[edge]
evec1 = J10 * trev[edge][0] + J11 * trev[edge][1]
elen[edge] = sqrt(evec0*evec0 + evec1*evec1)  = sqrt((J*trev[edge])**2)

trev = triangle_reference_edge_vector
evec0 = J00 * trev[edge][0] + J01 * trev[edge][1]  =  J*trev
evec1 = J10 * trev[edge][0] + J11 * trev[edge][1]
evec2 = J20 * trev[edge][0] + J21 * trev[edge][1] # Manifold: triangle in 3D
elen[edge] = sqrt(evec0*evec0 + evec1*evec1 + evec2*evec2)  = sqrt((J*trev[edge])**2)

trev = tetrahedron_reference_edge_vector
evec0 = sum(J[0,k] * trev[edge][k] for k in range(3))
evec1 = sum(J[1,k] * trev[edge][k] for k in range(3))
evec2 = sum(J[2,k] * trev[edge][k] for k in range(3))
elen[edge] = sqrt(evec0*evec0 + evec1*evec1 + evec2*evec2)  = sqrt((J*trev[edge])**2)


# DONE Here's how to compute min/max facet edge length:
triangle:
  r = facetarea
tetrahedron:
  min(elen[edge] for edge in range(6))
or
  min( min(elen[0], min(elen[1], elen[2])), min(elen[3], min(elen[4], elen[5])) )
or
  min1 = min_value(elen[0], min_value(elen[1], elen[2]))
  min2 = min_value(elen[3], min_value(elen[4], elen[5]))
  r = min_value(min1, min2)
(want proper Min/Max types for this!)


# DONE Here's how to compute circumradius for an interval:
circumradius_interval = cellvolume / 2


# DONE Here's how to compute circumradius for a triangle:
e0 = elen[0]
e1 = elen[1]
e2 = elen[2]
circumradius_triangle = (e0*e1*e2) / (4*cellvolume)


# DONE Here's how to compute circumradius for a tetrahedron:
# v1v2 = edge length between vertex 1 and 2
# la,lb,lc = lengths of the sides of an intermediate triangle
la = v1v2 * v0v3
lb = v0v2 * v1v3
lc = v0v1 * v2v3
# p = perimeter
p = (la + lb + lc)
# s = semiperimeter
s = p / 2
# area of intermediate triangle with Herons formula
tmp_area = sqrt(s * (s - la) * (s - lb) * (s - lc))
circumradius_tetrahedron = tmp_area / (6*cellvolume)

"""


class ChangeToReferenceGrad(MultiFunction):
    """Change to reference grad."""

    def __init__(self):
        """Initalise."""
        MultiFunction.__init__(self)

    expr = MultiFunction.reuse_if_untouched

    def terminal(self, o):
        """Apply to terminal."""
        return o

    def grad(self, o):
        """Apply to grad."""
        # Peel off the Grads and count them, and get restriction if
        # it's between the grad and the terminal
        ngrads = 0
        restricted = ""
        rv = False
        while not o._ufl_is_terminal_:
            if isinstance(o, Grad):
                (o,) = o.ufl_operands
                ngrads += 1
            elif isinstance(o, Restricted):
                restricted = o.side()
                (o,) = o.ufl_operands
            elif isinstance(o, ReferenceValue):
                rv = True
                (o,) = o.ufl_operands
            else:
                raise ValueError(f"Invalid type {type(o).__name__}")
        f = o
        if rv:
            f = ReferenceValue(f)

        # Get domain and create Jacobian inverse object
        domain = extract_unique_domain(o)
        Jinv = JacobianInverse(domain)

        if is_cellwise_constant(Jinv):
            # Optimise slightly by turning Grad(Grad(...)) into
            # J^(-T)J^(-T)RefGrad(RefGrad(...))
            # rather than J^(-T)RefGrad(J^(-T)RefGrad(...))

            # Create some new indices
            ii = indices(len(f.ufl_shape))  # Indices to get to the scalar component of f
            jj = indices(
                ngrads
            )  # Indices to sum over the local gradient axes with the inverse Jacobian
            kk = indices(ngrads)  # Indices for the leftover inverse Jacobian axes

            # Preserve restricted property
            if restricted:
                Jinv = Jinv(restricted)
                f = f(restricted)

            # Apply the same number of ReferenceGrad without mappings
            lgrad = f
            for i in range(ngrads):
                lgrad = ReferenceGrad(lgrad)

            # Apply mappings with scalar indexing operations (assumes
            # ReferenceGrad(Jinv) is zero)
            jinv_lgrad_f = lgrad[ii + jj]
            for j, k in zip(jj, kk):
                jinv_lgrad_f = Jinv[j, k] * jinv_lgrad_f

            # Wrap back in tensor shape, derivative axes at the end
            jinv_lgrad_f = as_tensor(jinv_lgrad_f, ii + kk)

        else:
            # J^(-T)RefGrad(J^(-T)RefGrad(...))

            # Preserve restricted property
            if restricted:
                Jinv = Jinv(restricted)
                f = f(restricted)

            jinv_lgrad_f = f
            for foo in range(ngrads):
                ii = indices(
                    len(jinv_lgrad_f.ufl_shape)
                )  # Indices to get to the scalar component of f
                j, k = indices(2)

                lgrad = ReferenceGrad(jinv_lgrad_f)
                jinv_lgrad_f = Jinv[j, k] * lgrad[ii + (j,)]

                # Wrap back in tensor shape, derivative axes at the end
                jinv_lgrad_f = as_tensor(jinv_lgrad_f, ii + (k,))

        return jinv_lgrad_f

    def reference_grad(self, o):
        """Apply to reference_grad."""
        raise ValueError("Not expecting reference grad while applying change to reference grad.")

    def coefficient_derivative(self, o):
        """Apply to coefficient_derivative."""
        raise ValueError(
            "Coefficient derivatives should be expanded before applying change to reference grad."
        )


def change_to_reference_grad(e):
    """Change Grad objects in expression to products of JacobianInverse and ReferenceGrad.

    Assumes the expression is preprocessed or at least that derivatives
    have been expanded.

    Args:
        e: An Expr or Form.
    """
    mf = ChangeToReferenceGrad()
    return map_expr_dag(mf, e)


def change_integrand_geometry_representation(integrand, scale, integral_type):
    """Change integrand geometry to the right representations."""
    integrand = apply_function_pullbacks(integrand)

    integrand = change_to_reference_grad(integrand)

    integrand = integrand * scale

    if integral_type == "quadrature":
        physical_coordinates_known = True
    else:
        physical_coordinates_known = False
    integrand = apply_geometry_lowering(integrand, physical_coordinates_known)

    return integrand

1	"""Algorithm for replacing gradients in an expression."""
2
3	# Copyright (C) 2008-2016 Martin Sandve Alnæs
4	#
5	# This file is part of UFL (https://www.fenicsproject.org)
6	#
7	# SPDX-License-Identifier: LGPL-3.0-or-later
8
9	from ufl.algorithms.apply_function_pullbacks import apply_function_pullbacks	1✔
10	from ufl.algorithms.apply_geometry_lowering import apply_geometry_lowering	1✔
11	from ufl.checks import is_cellwise_constant	1✔
12	from ufl.classes import Grad, JacobianInverse, ReferenceGrad, ReferenceValue, Restricted	1✔
13	from ufl.core.multiindex import indices	1✔
14	from ufl.corealg.map_dag import map_expr_dag	1✔
15	from ufl.corealg.multifunction import MultiFunction	1✔
16	from ufl.domain import extract_unique_domain	1✔
17	from ufl.tensors import as_tensor	1✔
18
19	"""	1✔
20	# Some notes:
21	# Below, let v_i mean physical coordinate of vertex i and V_i mean the
22	reference cell coordinate of the same vertex.
23
24
25	# Add a type for CellVertices? Note that vertices must be computed in
26	linear cell cases! triangle_vertices[i,j] = component j of vertex i,
27	following ufc numbering conventions
28
29
30	# DONE Add a type for CellEdgeLengths? Note that these are only easy to
31	define in the linear cell case!
32	triangle_edge_lengths = [v1v2, v0v2, v0v1] # shape (3,)
33	tetrahedron_edge_lengths = [v0v1, v0v2, v0v3, v1v2, v1v3, v2v3] # shape (6,)
34
35
36	# DONE Here's how to compute edge lengths from the Jacobian:
37	J =[ [dx0/dX0, dx0/dX1],
38	[dx1/dX0, dx1/dX1] ]
39	# First compute the edge vector, which is constant for each edge: the
40	vector from one vertex to the other
41	reference_edge_vector_0 = V2 - V1 # Example! Add a type ReferenceEdgeVectors?
42	# Then apply J to it and take the length of the resulting vector, this is generic for affine cells
43	edge_length_i = \|\| dot(J, reference_edge_vector_i) \|\|
44
45	e2 = \|\| J[:,0] . < 1, 0> \|\| = \|\| J[:,0] \|\| = \|\| dx/dX0 \|\| = edge length of edge 2 (v0-v1)
46	e1 = \|\| J[:,1] . < 0, 1> \|\| = \|\| J[:,1] \|\| = \|\| dx/dX1 \|\| = edge length of edge 1 (v0-v2)
47	e0 = \|\| J[:,:] . <-1, 1> \|\| = \|\| < J[0,1]-J[0,0], J[1,1]-J[1,0] > \|\| = \|\| dx/dX <-1,1> \|\|
48	= edge length of edge 0 (v1-v2)
49
50	trev = triangle_reference_edge_vector
51	evec0 = J00 * trev[edge][0] + J01 * trev[edge][1] = J*trev[edge]
52	evec1 = J10 * trev[edge][0] + J11 * trev[edge][1]
53	elen[edge] = sqrt(evec0evec0 + evec1evec1) = sqrt((Jtrev[edge])*2)
54
55	trev = triangle_reference_edge_vector
56	evec0 = J00 * trev[edge][0] + J01 * trev[edge][1] = J*trev
57	evec1 = J10 * trev[edge][0] + J11 * trev[edge][1]
58	evec2 = J20 * trev[edge][0] + J21 * trev[edge][1] # Manifold: triangle in 3D
59	elen[edge] = sqrt(evec0evec0 + evec1evec1 + evec2evec2) = sqrt((Jtrev[edge])**2)
60
61	trev = tetrahedron_reference_edge_vector
62	evec0 = sum(J[0,k] * trev[edge][k] for k in range(3))
63	evec1 = sum(J[1,k] * trev[edge][k] for k in range(3))
64	evec2 = sum(J[2,k] * trev[edge][k] for k in range(3))
65	elen[edge] = sqrt(evec0evec0 + evec1evec1 + evec2evec2) = sqrt((Jtrev[edge])**2)
66
67
68	# DONE Here's how to compute min/max facet edge length:
69	triangle:
70	r = facetarea
71	tetrahedron:
72	min(elen[edge] for edge in range(6))
73	or
74	min( min(elen[0], min(elen[1], elen[2])), min(elen[3], min(elen[4], elen[5])) )
75	or
76	min1 = min_value(elen[0], min_value(elen[1], elen[2]))
77	min2 = min_value(elen[3], min_value(elen[4], elen[5]))
78	r = min_value(min1, min2)
79	(want proper Min/Max types for this!)
80
81
82	# DONE Here's how to compute circumradius for an interval:
83	circumradius_interval = cellvolume / 2
84
85
86	# DONE Here's how to compute circumradius for a triangle:
87	e0 = elen[0]
88	e1 = elen[1]
89	e2 = elen[2]
90	circumradius_triangle = (e0e1e2) / (4*cellvolume)
91
92
93	# DONE Here's how to compute circumradius for a tetrahedron:
94	# v1v2 = edge length between vertex 1 and 2
95	# la,lb,lc = lengths of the sides of an intermediate triangle
96	la = v1v2 * v0v3
97	lb = v0v2 * v1v3
98	lc = v0v1 * v2v3
99	# p = perimeter
100	p = (la + lb + lc)
101	# s = semiperimeter
102	s = p / 2
103	# area of intermediate triangle with Herons formula
104	tmp_area = sqrt(s * (s - la) * (s - lb) * (s - lc))
105	circumradius_tetrahedron = tmp_area / (6*cellvolume)
106
107	"""
108
109
110	class ChangeToReferenceGrad(MultiFunction):	1✔
111	"""Change to reference grad."""
112
113	def __init__(self):	1✔
114	"""Initalise."""
115	MultiFunction.__init__(self)	1✔
116
117	expr = MultiFunction.reuse_if_untouched	1✔
118
119	def terminal(self, o):	1✔
120	"""Apply to terminal."""
121	return o	×
122
123	def grad(self, o):	1✔
124	"""Apply to grad."""
125	# Peel off the Grads and count them, and get restriction if
126	# it's between the grad and the terminal
127	ngrads = 0	1✔
128	restricted = ""	1✔
129	rv = False	1✔
130	while not o._ufl_is_terminal_:	1✔
131	if isinstance(o, Grad):	1✔
132	(o,) = o.ufl_operands	1✔
133	ngrads += 1	1✔
134	elif isinstance(o, Restricted):	×
135	restricted = o.side()	×
136	(o,) = o.ufl_operands	×
137	elif isinstance(o, ReferenceValue):	×
138	rv = True	×
139	(o,) = o.ufl_operands	×
140	else:
NEW 141	raise ValueError(f"Invalid type {type(o).__name__}")	×
142	f = o	1✔
143	if rv:	1✔
144	f = ReferenceValue(f)	×
145
146	# Get domain and create Jacobian inverse object
147	domain = extract_unique_domain(o)	1✔
148	Jinv = JacobianInverse(domain)	1✔
149
150	if is_cellwise_constant(Jinv):	1✔
151	# Optimise slightly by turning Grad(Grad(...)) into
152	# J^(-T)J^(-T)RefGrad(RefGrad(...))
153	# rather than J^(-T)RefGrad(J^(-T)RefGrad(...))
154
155	# Create some new indices
156	ii = indices(len(f.ufl_shape)) # Indices to get to the scalar component of f	1✔
157	jj = indices(	1✔
158	ngrads
159	) # Indices to sum over the local gradient axes with the inverse Jacobian
160	kk = indices(ngrads) # Indices for the leftover inverse Jacobian axes	1✔
161
162	# Preserve restricted property
163	if restricted:	1✔
164	Jinv = Jinv(restricted)	×
165	f = f(restricted)	×
166
167	# Apply the same number of ReferenceGrad without mappings
168	lgrad = f	1✔
169	for i in range(ngrads):	1✔
170	lgrad = ReferenceGrad(lgrad)	1✔
171
172	# Apply mappings with scalar indexing operations (assumes
173	# ReferenceGrad(Jinv) is zero)
174	jinv_lgrad_f = lgrad[ii + jj]	1✔
175	for j, k in zip(jj, kk):	1✔
176	jinv_lgrad_f = Jinv[j, k] * jinv_lgrad_f	1✔
177
178	# Wrap back in tensor shape, derivative axes at the end
179	jinv_lgrad_f = as_tensor(jinv_lgrad_f, ii + kk)	1✔
180
181	else:
182	# J^(-T)RefGrad(J^(-T)RefGrad(...))
183
184	# Preserve restricted property
185	if restricted:	×
186	Jinv = Jinv(restricted)	×
187	f = f(restricted)	×
188
189	jinv_lgrad_f = f	×
190	for foo in range(ngrads):	×
191	ii = indices(	×
192	len(jinv_lgrad_f.ufl_shape)
193	) # Indices to get to the scalar component of f
194	j, k = indices(2)	×
195
196	lgrad = ReferenceGrad(jinv_lgrad_f)	×
197	jinv_lgrad_f = Jinv[j, k] * lgrad[ii + (j,)]	×
198
199	# Wrap back in tensor shape, derivative axes at the end
200	jinv_lgrad_f = as_tensor(jinv_lgrad_f, ii + (k,))	×
201
202	return jinv_lgrad_f	1✔
203
204	def reference_grad(self, o):	1✔
205	"""Apply to reference_grad."""
206	raise ValueError("Not expecting reference grad while applying change to reference grad.")	×
207
208	def coefficient_derivative(self, o):	1✔
209	"""Apply to coefficient_derivative."""
210	raise ValueError(	×
211	"Coefficient derivatives should be expanded before applying change to reference grad."
212	)
213
214
215	def change_to_reference_grad(e):	1✔
216	"""Change Grad objects in expression to products of JacobianInverse and ReferenceGrad.
217
218	Assumes the expression is preprocessed or at least that derivatives
219	have been expanded.
220
221	Args:
222	e: An Expr or Form.
223	"""
224	mf = ChangeToReferenceGrad()	1✔
225	return map_expr_dag(mf, e)	1✔
226
227
228	def change_integrand_geometry_representation(integrand, scale, integral_type):	1✔
229	"""Change integrand geometry to the right representations."""
230	integrand = apply_function_pullbacks(integrand)	×
231
232	integrand = change_to_reference_grad(integrand)	×
233
234	integrand = integrand * scale	×
235
236	if integral_type == "quadrature":	×
237	physical_coordinates_known = True	×
238	else:
239	physical_coordinates_known = False	×
240	integrand = apply_geometry_lowering(integrand, physical_coordinates_known)	×
241
242	return integrand	×

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