17702620171

Committed 13 Sep 2025 10:07PM UTC coverage: 75.985% (+0.07%) from 75.917%

Build # 17702620171

Build Type

Pull #385

github

Committed by

schnellerhase

Commit Message

One more

Pull Request Pull Request #385: Change `AbstractCell` members to properties

Run Details

80 of 121 new or added lines in 11 files covered. (66.12%)

5 existing lines in 1 file now uncovered.

8970 of 11805 relevant lines covered (75.98%)

0.76 hits per line

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

36.52

/ufl/algorithms/apply_geometry_lowering.py

"""Algorithm for lowering abstractions of geometric types.

This means replacing high-level types with expressions
of mostly the Jacobian and reference cell data.
"""

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

import warnings
from functools import reduce
from itertools import combinations

from ufl.classes import (
    CellCoordinate,
    CellEdgeVectors,
    CellFacetJacobian,
    CellOrientation,
    CellOrigin,
    CellRidgeJacobian,
    CellVertices,
    CellVolume,
    Expr,
    FacetEdgeVectors,
    FacetJacobian,
    FacetJacobianDeterminant,
    FloatValue,
    Form,
    Integral,
    Jacobian,
    JacobianDeterminant,
    JacobianInverse,
    MaxCellEdgeLength,
    ReferenceCellVolume,
    ReferenceFacetVolume,
    ReferenceGrad,
    ReferenceNormal,
    RidgeJacobian,
    SpatialCoordinate,
)
from ufl.compound_expressions import cross_expr, determinant_expr, inverse_expr
from ufl.core.multiindex import Index, indices
from ufl.corealg.map_dag import map_expr_dag
from ufl.corealg.multifunction import MultiFunction, memoized_handler
from ufl.domain import extract_unique_domain
from ufl.measure import custom_integral_types, point_integral_types
from ufl.operators import conj, max_value, min_value, real, sqrt
from ufl.tensors import as_tensor, as_vector


class GeometryLoweringApplier(MultiFunction):
    """Geometry lowering."""

    def __init__(self, preserve_types=()):
        """Initialise."""
        MultiFunction.__init__(self)
        # Store preserve_types as boolean lookup table
        self._preserve_types = [False] * Expr._ufl_num_typecodes_
        for cls in preserve_types:
            self._preserve_types[cls._ufl_typecode_] = True

    expr = MultiFunction.reuse_if_untouched

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

    @memoized_handler
    def jacobian(self, o):
        """Apply to jacobian."""
        if self._preserve_types[o._ufl_typecode_]:
            return o
        domain = extract_unique_domain(o)
        if not domain.ufl_coordinate_element().pullback.is_identity:
            raise ValueError("Piola mapped coordinates are not implemented.")
        # Note: No longer supporting domain.coordinates(), always
        # preserving SpatialCoordinate object.  However if Jacobians
        # are not preserved, using
        # ReferenceGrad(SpatialCoordinate(domain)) to represent them.
        x = self.spatial_coordinate(SpatialCoordinate(domain))
        return ReferenceGrad(x)

    @memoized_handler
    def _future_jacobian(self, o):
        """Apply to _future_jacobian."""
        # If we're not using Coefficient to represent the spatial
        # coordinate, we can just as well just return o here too
        # unless we add representation of basis functions and dofs to
        # the ufl layer (which is nice to avoid).
        return o

    @memoized_handler
    def jacobian_inverse(self, o):
        """Apply to jacobian_inverse."""
        if self._preserve_types[o._ufl_typecode_]:
            return o

        domain = extract_unique_domain(o)
        J = self.jacobian(Jacobian(domain))
        # TODO: This could in principle use
        # preserve_types[JacobianDeterminant] with minor refactoring:
        K = inverse_expr(J)
        return K

    @memoized_handler
    def jacobian_determinant(self, o):
        """Apply to jacobian_determinant."""
        if self._preserve_types[o._ufl_typecode_]:
            return o

        domain = extract_unique_domain(o)
        J = self.jacobian(Jacobian(domain))
        detJ = determinant_expr(J)

        # TODO: Is "signing" the determinant for manifolds the
        #       cleanest approach?  The alternative is to have a
        #       specific type for the unsigned pseudo-determinant.
        if domain.topological_dimension < domain.geometric_dimension:
            co = CellOrientation(domain)
            detJ = co * detJ
        return detJ

    @memoized_handler
    def facet_jacobian(self, o):
        """Apply to facet_jacobian."""
        if self._preserve_types[o._ufl_typecode_]:
            return o

        domain = extract_unique_domain(o)
        J = self.jacobian(Jacobian(domain))
        RFJ = CellFacetJacobian(domain)
        i, j, k = indices(3)
        return as_tensor(J[i, k] * RFJ[k, j], (i, j))

    @memoized_handler
    def facet_jacobian_inverse(self, o):
        """Apply to facet_jacobian_inverse."""
        if self._preserve_types[o._ufl_typecode_]:
            return o

        domain = extract_unique_domain(o)
        FJ = self.facet_jacobian(FacetJacobian(domain))
        # This could in principle use
        # preserve_types[JacobianDeterminant] with minor refactoring:
        return inverse_expr(FJ)

    @memoized_handler
    def facet_jacobian_determinant(self, o):
        """Apply to facet_jacobian_determinant."""
        if self._preserve_types[o._ufl_typecode_]:
            return o

        domain = extract_unique_domain(o)
        FJ = self.facet_jacobian(FacetJacobian(domain))
        detFJ = determinant_expr(FJ)

        # TODO: Should we "sign" the facet jacobian determinant for
        #       manifolds?  It's currently used unsigned in
        #       apply_integral_scaling.
        # if domain.topological_dimension() < domain.geometric_dimension():
        #     co = CellOrientation(domain)
        #     detFJ = co*detFJ

        return detFJ

    @memoized_handler
    def ridge_jacobian(self, o):
        """Apply to ridge_jacobian."""
        if self._preserve_types[o._ufl_typecode_]:
            return o

        domain = extract_unique_domain(o)
        J = self.jacobian(Jacobian(domain))
        REJ = CellRidgeJacobian(domain)
        i, j, k = indices(3)
        return as_tensor(J[i, k] * REJ[k, j], (i, j))

    @memoized_handler
    def ridge_jacobian_inverse(self, o):
        """Apply to edge_jacobian_inverse."""
        if self._preserve_types[o._ufl_typecode_]:
            return o

        domain = extract_unique_domain(o)
        EJ = self.ridge_jacobian(RidgeJacobian(domain))
        return inverse_expr(EJ)

    @memoized_handler
    def ridge_jacobian_determinant(self, o):
        """Apply to edge_jacobian_determinant."""
        if self._preserve_types[o._ufl_typecode_]:
            return o

        domain = extract_unique_domain(o)
        EJ = self.ridge_jacobian(RidgeJacobian(domain))
        detEJ = determinant_expr(EJ)
        return detEJ

    @memoized_handler
    def spatial_coordinate(self, o):
        """Apply to spatial_coordinate.

        Fall through to coordinate field of domain if it exists.
        """
        if self._preserve_types[o._ufl_typecode_]:
            return o
        if not extract_unique_domain(o).ufl_coordinate_element().pullback.is_identity:
            raise ValueError("Piola mapped coordinates are not implemented.")
        # No longer supporting domain.coordinates(), always preserving
        # SpatialCoordinate object.
        return o

    @memoized_handler
    def cell_coordinate(self, o):
        """Apply to cell_coordinate.

        Compute from physical coordinates if they are known, using the appropriate mappings.
        """
        if self._preserve_types[o._ufl_typecode_]:
            return o

        domain = extract_unique_domain(o)
        K = self.jacobian_inverse(JacobianInverse(domain))
        x = self.spatial_coordinate(SpatialCoordinate(domain))
        x0 = CellOrigin(domain)
        i, j = indices(2)
        X = as_tensor(K[i, j] * (x[j] - x0[j]), (i,))
        return X

    @memoized_handler
    def facet_cell_coordinate(self, o):
        """Apply to facet_cell_coordinate."""
        if self._preserve_types[o._ufl_typecode_]:
            return o

        raise ValueError(
            "Missing computation of facet reference coordinates "
            "from physical coordinates via mappings."
        )

    @memoized_handler
    def cell_volume(self, o):
        """Apply to cell_volume."""
        if self._preserve_types[o._ufl_typecode_]:
            return o

        domain = extract_unique_domain(o)
        if not domain.is_piecewise_linear_simplex_domain():
            # Don't lower for non-affine cells, instead leave it to
            # form compiler
            warnings.warn("Only know how to compute the cell volume of an affine cell.")
            return o

        r = self.jacobian_determinant(JacobianDeterminant(domain))
        r0 = ReferenceCellVolume(domain)
        return abs(r * r0)

    @memoized_handler
    def facet_area(self, o):
        """Apply to facet_area."""
        if self._preserve_types[o._ufl_typecode_]:
            return o

        domain = extract_unique_domain(o)
        tdim = domain.topological_dimension
        if not domain.is_piecewise_linear_simplex_domain():
            # Don't lower for non-affine cells, instead leave it to
            # form compiler
            warnings.warn("Only know how to compute the facet area of an affine cell.")
            return o

        # Area of "facet" of interval (i.e. "area" of a vertex) is defined as 1.0
        if tdim == 1:
            return FloatValue(1.0)

        r = self.facet_jacobian_determinant(FacetJacobianDeterminant(domain))
        r0 = ReferenceFacetVolume(domain)
        return abs(r * r0)

    @memoized_handler
    def circumradius(self, o):
        """Apply to circumradius."""
        if self._preserve_types[o._ufl_typecode_]:
            return o

        domain = extract_unique_domain(o)

        if not domain.is_piecewise_linear_simplex_domain():
            raise ValueError("Circumradius only makes sense for affine simplex cells")

        cellname = domain.ufl_cell().cellname
        cellvolume = self.cell_volume(CellVolume(domain))

        if cellname == "interval":
            # Optimization for square interval; no square root needed
            return 0.5 * cellvolume

        # Compute lengths of cell edges
        edges = CellEdgeVectors(domain)
        num_edges = edges.ufl_shape[0]
        j = Index()
        elen = [real(sqrt(real(edges[e, j] * conj(edges[e, j])))) for e in range(num_edges)]

        if cellname == "triangle":
            return (elen[0] * elen[1] * elen[2]) / (4.0 * cellvolume)

        elif cellname == "tetrahedron":
            # la, lb, lc = lengths of the sides of an intermediate triangle
            # NOTE: Is here some hidden numbering assumption?
            la = elen[3] * elen[2]
            lb = elen[4] * elen[1]
            lc = elen[5] * elen[0]
            # p = perimeter
            p = la + lb + lc
            # s = semiperimeter
            s = p / 2
            # area of intermediate triangle with Herons formula
            triangle_area = sqrt(s * (s - la) * (s - lb) * (s - lc))
            return triangle_area / (6.0 * cellvolume)

    @memoized_handler
    def max_cell_edge_length(self, o):
        """Apply to max_cell_edge_length."""
        return self._reduce_cell_edge_length(o, max_value)

    @memoized_handler
    def min_cell_edge_length(self, o):
        """Apply to min_cell_edge_length."""
        return self._reduce_cell_edge_length(o, min_value)

    def _reduce_cell_edge_length(self, o, reduction_op):
        """Apply to _reduce_cell_edge_length."""
        if self._preserve_types[o._ufl_typecode_]:
            return o

        domain = extract_unique_domain(o)

        if domain.ufl_coordinate_element().embedded_subdegree > 1:
            # Don't lower bendy cells, instead leave it to form compiler
            warnings.warn("Only know how to compute cell edge lengths of P1 or Q1 cell.")
            return o

        elif domain.ufl_cell().cellname == "interval":
            # Interval optimization, square root not needed
            return self.cell_volume(CellVolume(domain))

        else:
            # Other P1 or Q1 cells
            edges = CellEdgeVectors(domain)
            num_edges = edges.ufl_shape[0]
            j = Index()
            elen2 = [real(edges[e, j] * conj(edges[e, j])) for e in range(num_edges)]
            return real(sqrt(reduce(reduction_op, elen2)))

    @memoized_handler
    def cell_diameter(self, o):
        """Apply to cell_diameter."""
        if self._preserve_types[o._ufl_typecode_]:
            return o

        domain = extract_unique_domain(o)

        if domain.ufl_coordinate_element().embedded_subdegree > 1:
            # Don't lower bendy cells, instead leave it to form compiler
            warnings.warn("Only know how to compute cell diameter of P1 or Q1 cell.")
            return o

        elif domain.is_piecewise_linear_simplex_domain():
            # Simplices
            return self.max_cell_edge_length(MaxCellEdgeLength(domain))

        else:
            # Q1 cells, maximal distance between any two vertices
            verts = CellVertices(domain)
            verts = [verts[v, ...] for v in range(verts.ufl_shape[0])]
            j = Index()
            elen2 = (real((v0 - v1)[j] * conj((v0 - v1)[j])) for v0, v1 in combinations(verts, 2))
            return real(sqrt(reduce(max_value, elen2)))

    @memoized_handler
    def max_facet_edge_length(self, o):
        """Apply to max_facet_edge_length."""
        return self._reduce_facet_edge_length(o, max_value)

    @memoized_handler
    def min_facet_edge_length(self, o):
        """Apply to min_facet_edge_length."""
        return self._reduce_facet_edge_length(o, min_value)

    def _reduce_facet_edge_length(self, o, reduction_op):
        """Apply to _reduce_facet_edge_length."""
        if self._preserve_types[o._ufl_typecode_]:
            return o

        domain = extract_unique_domain(o)

        if domain.ufl_cell().topological_dimension < 3:
            raise ValueError("Facet edge lengths only make sense for topological dimension >= 3.")

        elif domain.ufl_coordinate_element().embedded_subdegree > 1:
            # Don't lower bendy cells, instead leave it to form compiler
            warnings.warn("Only know how to compute facet edge lengths of P1 or Q1 cell.")
            return o

        else:
            # P1 tetrahedron or Q1 hexahedron
            edges = FacetEdgeVectors(domain)
            num_edges = edges.ufl_shape[0]
            j = Index()
            elen2 = [real(edges[e, j] * conj(edges[e, j])) for e in range(num_edges)]
            return real(sqrt(reduce(reduction_op, elen2)))

    @memoized_handler
    def cell_normal(self, o):
        """Apply to cell_normal."""
        if self._preserve_types[o._ufl_typecode_]:
            return o

        domain = extract_unique_domain(o)
        gdim = domain.geometric_dimension()
        tdim = domain.topological_dimension()

        if tdim == gdim - 1:  # n-manifold embedded in n-1 space
            i = Index()
            J = self.jacobian(Jacobian(domain))

            if tdim == 2:
                # Surface in 3D
                t0 = as_vector(J[i, 0], i)
                t1 = as_vector(J[i, 1], i)
                cell_normal = cross_expr(t0, t1)
            elif tdim == 1:
                # Line in 2D (cell normal is 'up' for a line pointing
                # to the 'right')
                cell_normal = as_vector((-J[1, 0], J[0, 0]))
            else:
                raise ValueError(f"Cell normal not implemented for tdim {tdim}, gdim {gdim}")

            # Return normalized vector, sign corrected by cell
            # orientation
            co = CellOrientation(domain)
            return co * cell_normal / sqrt(cell_normal[i] * cell_normal[i])
        else:
            raise ValueError(f"Cell normal undefined for tdim {tdim}, gdim {gdim}")

    @memoized_handler
    def facet_normal(self, o):
        """Apply to facet_normal."""
        if self._preserve_types[o._ufl_typecode_]:
            return o

        domain = extract_unique_domain(o)
        tdim = domain.topological_dimension

        if tdim == 1:
            # Special-case 1D (possibly immersed), for which we say
            # that n is just in the direction of J.
            J = self.jacobian(Jacobian(domain))  # dx/dX
            ndir = J[:, 0]

            gdim = domain.geometric_dimension
            if gdim == 1:
                nlen = abs(ndir[0])
            else:
                i = Index()
                nlen = sqrt(ndir[i] * ndir[i])

            rn = ReferenceNormal(domain)  # +/- 1.0 here
            n = rn[0] * ndir / nlen
            r = n
        else:
            # Recall that the covariant Piola transform u -> J^(-T)*u
            # preserves tangential components. The normal vector is
            # characterised by having zero tangential component in
            # reference and physical space.
            Jinv = self.jacobian_inverse(JacobianInverse(domain))
            i, j = indices(2)

            rn = ReferenceNormal(domain)
            # compute signed, unnormalised normal; note transpose
            ndir = as_vector(Jinv[j, i] * rn[j], i)

            # normalise
            i = Index()
            n = ndir / sqrt(ndir[i] * ndir[i])
            r = n

        if r.ufl_shape != o.ufl_shape:
            raise ValueError(
                f"Inconsistent dimensions (in={o.ufl_shape[0]}, out={r.ufl_shape[0]})."
            )
        return r


def apply_geometry_lowering(form, preserve_types=()):
    """Change GeometricQuantity objects in expression to the lowest level GeometricQuantity objects.

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

    Args:
        form: An Expr or Form.
        preserve_types: Preserved types
    """
    if isinstance(form, Form):
        newintegrals = [
            apply_geometry_lowering(integral, preserve_types) for integral in form.integrals()
        ]
        return Form(newintegrals)

    elif isinstance(form, Integral):
        integral = form
        if integral.integral_type() in (custom_integral_types + point_integral_types):
            automatic_preserve_types = [SpatialCoordinate, Jacobian]
        else:
            automatic_preserve_types = [CellCoordinate]
        preserve_types = set(preserve_types) | set(automatic_preserve_types)

        mf = GeometryLoweringApplier(preserve_types)
        newintegrand = map_expr_dag(mf, integral.integrand())
        return integral.reconstruct(integrand=newintegrand)

    elif isinstance(form, Expr):
        expr = form
        mf = GeometryLoweringApplier(preserve_types)
        return map_expr_dag(mf, expr)

    else:
        raise ValueError(f"Invalid type {form.__class__.__name__}")

1	"""Algorithm for lowering abstractions of geometric types.
2
3	This means replacing high-level types with expressions
4	of mostly the Jacobian and reference cell data.
5	"""
6
7	# Copyright (C) 2013-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	import warnings	1✔
14	from functools import reduce	1✔
15	from itertools import combinations	1✔
16
17	from ufl.classes import (	1✔
18	CellCoordinate,
19	CellEdgeVectors,
20	CellFacetJacobian,
21	CellOrientation,
22	CellOrigin,
23	CellRidgeJacobian,
24	CellVertices,
25	CellVolume,
26	Expr,
27	FacetEdgeVectors,
28	FacetJacobian,
29	FacetJacobianDeterminant,
30	FloatValue,
31	Form,
32	Integral,
33	Jacobian,
34	JacobianDeterminant,
35	JacobianInverse,
36	MaxCellEdgeLength,
37	ReferenceCellVolume,
38	ReferenceFacetVolume,
39	ReferenceGrad,
40	ReferenceNormal,
41	RidgeJacobian,
42	SpatialCoordinate,
43	)
44	from ufl.compound_expressions import cross_expr, determinant_expr, inverse_expr	1✔
45	from ufl.core.multiindex import Index, indices	1✔
46	from ufl.corealg.map_dag import map_expr_dag	1✔
47	from ufl.corealg.multifunction import MultiFunction, memoized_handler	1✔
48	from ufl.domain import extract_unique_domain	1✔
49	from ufl.measure import custom_integral_types, point_integral_types	1✔
50	from ufl.operators import conj, max_value, min_value, real, sqrt	1✔
51	from ufl.tensors import as_tensor, as_vector	1✔
52
53
54	class GeometryLoweringApplier(MultiFunction):	1✔
55	"""Geometry lowering."""
56
57	def __init__(self, preserve_types=()):	1✔
58	"""Initialise."""
59	MultiFunction.__init__(self)	1✔
60	# Store preserve_types as boolean lookup table
61	self._preserve_types = [False] * Expr._ufl_num_typecodes_	1✔
62	for cls in preserve_types:	1✔
63	self._preserve_types[cls._ufl_typecode_] = True	1✔
64
65	expr = MultiFunction.reuse_if_untouched	1✔
66
67	def terminal(self, t):	1✔
68	"""Apply to terminal."""
69	return t	1✔
70
71	@memoized_handler	1✔
72	def jacobian(self, o):	1✔
73	"""Apply to jacobian."""
74	if self._preserve_types[o._ufl_typecode_]:	1✔
75	return o	1✔
76	domain = extract_unique_domain(o)	1✔
77	if not domain.ufl_coordinate_element().pullback.is_identity:	1✔
78	raise ValueError("Piola mapped coordinates are not implemented.")	×
79	# Note: No longer supporting domain.coordinates(), always
80	# preserving SpatialCoordinate object. However if Jacobians
81	# are not preserved, using
82	# ReferenceGrad(SpatialCoordinate(domain)) to represent them.
83	x = self.spatial_coordinate(SpatialCoordinate(domain))	1✔
84	return ReferenceGrad(x)	1✔
85
86	@memoized_handler	1✔
87	def _future_jacobian(self, o):	1✔
88	"""Apply to _future_jacobian."""
89	# If we're not using Coefficient to represent the spatial
90	# coordinate, we can just as well just return o here too
91	# unless we add representation of basis functions and dofs to
92	# the ufl layer (which is nice to avoid).
93	return o	×
94
95	@memoized_handler	1✔
96	def jacobian_inverse(self, o):	1✔
97	"""Apply to jacobian_inverse."""
98	if self._preserve_types[o._ufl_typecode_]:	1✔
99	return o	×
100
101	domain = extract_unique_domain(o)	1✔
102	J = self.jacobian(Jacobian(domain))	1✔
103	# TODO: This could in principle use
104	# preserve_types[JacobianDeterminant] with minor refactoring:
105	K = inverse_expr(J)	1✔
106	return K	1✔
107
108	@memoized_handler	1✔
109	def jacobian_determinant(self, o):	1✔
110	"""Apply to jacobian_determinant."""
111	if self._preserve_types[o._ufl_typecode_]:	1✔
112	return o	×
113
114	domain = extract_unique_domain(o)	1✔
115	J = self.jacobian(Jacobian(domain))	1✔
116	detJ = determinant_expr(J)	1✔
117
118	# TODO: Is "signing" the determinant for manifolds the
119	# cleanest approach? The alternative is to have a
120	# specific type for the unsigned pseudo-determinant.
121	if domain.topological_dimension < domain.geometric_dimension:	1✔
122	co = CellOrientation(domain)	1✔
123	detJ = co * detJ	1✔
124	return detJ	1✔
125
126	@memoized_handler	1✔
127	def facet_jacobian(self, o):	1✔
128	"""Apply to facet_jacobian."""
129	if self._preserve_types[o._ufl_typecode_]:	×
130	return o	×
131
132	domain = extract_unique_domain(o)	×
133	J = self.jacobian(Jacobian(domain))	×
134	RFJ = CellFacetJacobian(domain)	×
135	i, j, k = indices(3)	×
136	return as_tensor(J[i, k] * RFJ[k, j], (i, j))	×
137
138	@memoized_handler	1✔
139	def facet_jacobian_inverse(self, o):	1✔
140	"""Apply to facet_jacobian_inverse."""
141	if self._preserve_types[o._ufl_typecode_]:	×
142	return o	×
143
144	domain = extract_unique_domain(o)	×
145	FJ = self.facet_jacobian(FacetJacobian(domain))	×
146	# This could in principle use
147	# preserve_types[JacobianDeterminant] with minor refactoring:
148	return inverse_expr(FJ)	×
149
150	@memoized_handler	1✔
151	def facet_jacobian_determinant(self, o):	1✔
152	"""Apply to facet_jacobian_determinant."""
153	if self._preserve_types[o._ufl_typecode_]:	×
154	return o	×
155
156	domain = extract_unique_domain(o)	×
157	FJ = self.facet_jacobian(FacetJacobian(domain))	×
158	detFJ = determinant_expr(FJ)	×
159
160	# TODO: Should we "sign" the facet jacobian determinant for
161	# manifolds? It's currently used unsigned in
162	# apply_integral_scaling.
163	# if domain.topological_dimension() < domain.geometric_dimension():
164	# co = CellOrientation(domain)
165	# detFJ = co*detFJ
166
167	return detFJ	×
168
169	@memoized_handler	1✔
170	def ridge_jacobian(self, o):	1✔
171	"""Apply to ridge_jacobian."""
172	if self._preserve_types[o._ufl_typecode_]:	×
173	return o	×
174
175	domain = extract_unique_domain(o)	×
176	J = self.jacobian(Jacobian(domain))	×
177	REJ = CellRidgeJacobian(domain)	×
178	i, j, k = indices(3)	×
179	return as_tensor(J[i, k] * REJ[k, j], (i, j))	×
180
181	@memoized_handler	1✔
182	def ridge_jacobian_inverse(self, o):	1✔
183	"""Apply to edge_jacobian_inverse."""
184	if self._preserve_types[o._ufl_typecode_]:	×
185	return o	×
186
187	domain = extract_unique_domain(o)	×
188	EJ = self.ridge_jacobian(RidgeJacobian(domain))	×
189	return inverse_expr(EJ)	×
190
191	@memoized_handler	1✔
192	def ridge_jacobian_determinant(self, o):	1✔
193	"""Apply to edge_jacobian_determinant."""
194	if self._preserve_types[o._ufl_typecode_]:	×
195	return o	×
196
197	domain = extract_unique_domain(o)	×
198	EJ = self.ridge_jacobian(RidgeJacobian(domain))	×
199	detEJ = determinant_expr(EJ)	×
200	return detEJ	×
201
202	@memoized_handler	1✔
203	def spatial_coordinate(self, o):	1✔
204	"""Apply to spatial_coordinate.
205
206	Fall through to coordinate field of domain if it exists.
207	"""
208	if self._preserve_types[o._ufl_typecode_]:	1✔
209	return o	×
210	if not extract_unique_domain(o).ufl_coordinate_element().pullback.is_identity:	1✔
211	raise ValueError("Piola mapped coordinates are not implemented.")	×
212	# No longer supporting domain.coordinates(), always preserving
213	# SpatialCoordinate object.
214	return o	1✔
215
216	@memoized_handler	1✔
217	def cell_coordinate(self, o):	1✔
218	"""Apply to cell_coordinate.
219
220	Compute from physical coordinates if they are known, using the appropriate mappings.
221	"""
222	if self._preserve_types[o._ufl_typecode_]:	×
223	return o	×
224
225	domain = extract_unique_domain(o)	×
226	K = self.jacobian_inverse(JacobianInverse(domain))	×
227	x = self.spatial_coordinate(SpatialCoordinate(domain))	×
228	x0 = CellOrigin(domain)	×
229	i, j = indices(2)	×
230	X = as_tensor(K[i, j] * (x[j] - x0[j]), (i,))	×
231	return X	×
232
233	@memoized_handler	1✔
234	def facet_cell_coordinate(self, o):	1✔
235	"""Apply to facet_cell_coordinate."""
236	if self._preserve_types[o._ufl_typecode_]:	×
237	return o	×
238
239	raise ValueError(	×
240	"Missing computation of facet reference coordinates "
241	"from physical coordinates via mappings."
242	)
243
244	@memoized_handler	1✔
245	def cell_volume(self, o):	1✔
246	"""Apply to cell_volume."""
247	if self._preserve_types[o._ufl_typecode_]:	×
248	return o	×
249
250	domain = extract_unique_domain(o)	×
251	if not domain.is_piecewise_linear_simplex_domain():	×
252	# Don't lower for non-affine cells, instead leave it to
253	# form compiler
254	warnings.warn("Only know how to compute the cell volume of an affine cell.")	×
255	return o	×
256
257	r = self.jacobian_determinant(JacobianDeterminant(domain))	×
258	r0 = ReferenceCellVolume(domain)	×
259	return abs(r * r0)	×
260
261	@memoized_handler	1✔
262	def facet_area(self, o):	1✔
263	"""Apply to facet_area."""
264	if self._preserve_types[o._ufl_typecode_]:	×
265	return o	×
266
267	domain = extract_unique_domain(o)	×
NEW 268	tdim = domain.topological_dimension	×
269	if not domain.is_piecewise_linear_simplex_domain():	×
270	# Don't lower for non-affine cells, instead leave it to
271	# form compiler
272	warnings.warn("Only know how to compute the facet area of an affine cell.")	×
273	return o	×
274
275	# Area of "facet" of interval (i.e. "area" of a vertex) is defined as 1.0
276	if tdim == 1:	×
277	return FloatValue(1.0)	×
278
279	r = self.facet_jacobian_determinant(FacetJacobianDeterminant(domain))	×
280	r0 = ReferenceFacetVolume(domain)	×
281	return abs(r * r0)	×
282
283	@memoized_handler	1✔
284	def circumradius(self, o):	1✔
285	"""Apply to circumradius."""
286	if self._preserve_types[o._ufl_typecode_]:	×
287	return o	×
288
289	domain = extract_unique_domain(o)	×
290
291	if not domain.is_piecewise_linear_simplex_domain():	×
292	raise ValueError("Circumradius only makes sense for affine simplex cells")	×
293
NEW 294	cellname = domain.ufl_cell().cellname	×
295	cellvolume = self.cell_volume(CellVolume(domain))	×
296
297	if cellname == "interval":	×
298	# Optimization for square interval; no square root needed
299	return 0.5 * cellvolume	×
300
301	# Compute lengths of cell edges
302	edges = CellEdgeVectors(domain)	×
303	num_edges = edges.ufl_shape[0]	×
304	j = Index()	×
305	elen = [real(sqrt(real(edges[e, j] * conj(edges[e, j])))) for e in range(num_edges)]	×
306
307	if cellname == "triangle":	×
308	return (elen[0] * elen[1] * elen[2]) / (4.0 * cellvolume)	×
309
310	elif cellname == "tetrahedron":	×
311	# la, lb, lc = lengths of the sides of an intermediate triangle
312	# NOTE: Is here some hidden numbering assumption?
313	la = elen[3] * elen[2]	×
314	lb = elen[4] * elen[1]	×
315	lc = elen[5] * elen[0]	×
316	# p = perimeter
317	p = la + lb + lc	×
318	# s = semiperimeter
319	s = p / 2	×
320	# area of intermediate triangle with Herons formula
321	triangle_area = sqrt(s * (s - la) * (s - lb) * (s - lc))	×
322	return triangle_area / (6.0 * cellvolume)	×
323
324	@memoized_handler	1✔
325	def max_cell_edge_length(self, o):	1✔
326	"""Apply to max_cell_edge_length."""
327	return self._reduce_cell_edge_length(o, max_value)	×
328
329	@memoized_handler	1✔
330	def min_cell_edge_length(self, o):	1✔
331	"""Apply to min_cell_edge_length."""
332	return self._reduce_cell_edge_length(o, min_value)	×
333
334	def _reduce_cell_edge_length(self, o, reduction_op):	1✔
335	"""Apply to _reduce_cell_edge_length."""
336	if self._preserve_types[o._ufl_typecode_]:	×
337	return o	×
338
339	domain = extract_unique_domain(o)	×
340
341	if domain.ufl_coordinate_element().embedded_subdegree > 1:	×
342	# Don't lower bendy cells, instead leave it to form compiler
343	warnings.warn("Only know how to compute cell edge lengths of P1 or Q1 cell.")	×
344	return o	×
345
NEW 346	elif domain.ufl_cell().cellname == "interval":	×
347	# Interval optimization, square root not needed
348	return self.cell_volume(CellVolume(domain))	×
349
350	else:
351	# Other P1 or Q1 cells
352	edges = CellEdgeVectors(domain)	×
353	num_edges = edges.ufl_shape[0]	×
354	j = Index()	×
355	elen2 = [real(edges[e, j] * conj(edges[e, j])) for e in range(num_edges)]	×
356	return real(sqrt(reduce(reduction_op, elen2)))	×
357
358	@memoized_handler	1✔
359	def cell_diameter(self, o):	1✔
360	"""Apply to cell_diameter."""
361	if self._preserve_types[o._ufl_typecode_]:	×
362	return o	×
363
364	domain = extract_unique_domain(o)	×
365
366	if domain.ufl_coordinate_element().embedded_subdegree > 1:	×
367	# Don't lower bendy cells, instead leave it to form compiler
368	warnings.warn("Only know how to compute cell diameter of P1 or Q1 cell.")	×
369	return o	×
370
371	elif domain.is_piecewise_linear_simplex_domain():	×
372	# Simplices
373	return self.max_cell_edge_length(MaxCellEdgeLength(domain))	×
374
375	else:
376	# Q1 cells, maximal distance between any two vertices
377	verts = CellVertices(domain)	×
378	verts = [verts[v, ...] for v in range(verts.ufl_shape[0])]	×
379	j = Index()	×
380	elen2 = (real((v0 - v1)[j] * conj((v0 - v1)[j])) for v0, v1 in combinations(verts, 2))	×
381	return real(sqrt(reduce(max_value, elen2)))	×
382
383	@memoized_handler	1✔
384	def max_facet_edge_length(self, o):	1✔
385	"""Apply to max_facet_edge_length."""
386	return self._reduce_facet_edge_length(o, max_value)	×
387
388	@memoized_handler	1✔
389	def min_facet_edge_length(self, o):	1✔
390	"""Apply to min_facet_edge_length."""
391	return self._reduce_facet_edge_length(o, min_value)	×
392
393	def _reduce_facet_edge_length(self, o, reduction_op):	1✔
394	"""Apply to _reduce_facet_edge_length."""
395	if self._preserve_types[o._ufl_typecode_]:	×
396	return o	×
397
398	domain = extract_unique_domain(o)	×
399
NEW 400	if domain.ufl_cell().topological_dimension < 3:	×
401	raise ValueError("Facet edge lengths only make sense for topological dimension >= 3.")	×
402
403	elif domain.ufl_coordinate_element().embedded_subdegree > 1:	×
404	# Don't lower bendy cells, instead leave it to form compiler
405	warnings.warn("Only know how to compute facet edge lengths of P1 or Q1 cell.")	×
406	return o	×
407
408	else:
409	# P1 tetrahedron or Q1 hexahedron
410	edges = FacetEdgeVectors(domain)	×
411	num_edges = edges.ufl_shape[0]	×
412	j = Index()	×
413	elen2 = [real(edges[e, j] * conj(edges[e, j])) for e in range(num_edges)]	×
414	return real(sqrt(reduce(reduction_op, elen2)))	×
415
416	@memoized_handler	1✔
417	def cell_normal(self, o):	1✔
418	"""Apply to cell_normal."""
419	if self._preserve_types[o._ufl_typecode_]:	×
420	return o	×
421
422	domain = extract_unique_domain(o)	×
423	gdim = domain.geometric_dimension()	×
424	tdim = domain.topological_dimension()	×
425
426	if tdim == gdim - 1: # n-manifold embedded in n-1 space	×
427	i = Index()	×
428	J = self.jacobian(Jacobian(domain))	×
429
430	if tdim == 2:	×
431	# Surface in 3D
432	t0 = as_vector(J[i, 0], i)	×
433	t1 = as_vector(J[i, 1], i)	×
434	cell_normal = cross_expr(t0, t1)	×
435	elif tdim == 1:	×
436	# Line in 2D (cell normal is 'up' for a line pointing
437	# to the 'right')
438	cell_normal = as_vector((-J[1, 0], J[0, 0]))	×
439	else:
440	raise ValueError(f"Cell normal not implemented for tdim {tdim}, gdim {gdim}")	×
441
442	# Return normalized vector, sign corrected by cell
443	# orientation
444	co = CellOrientation(domain)	×
445	return co * cell_normal / sqrt(cell_normal[i] * cell_normal[i])	×
446	else:
447	raise ValueError(f"Cell normal undefined for tdim {tdim}, gdim {gdim}")	×
448
449	@memoized_handler	1✔
450	def facet_normal(self, o):	1✔
451	"""Apply to facet_normal."""
452	if self._preserve_types[o._ufl_typecode_]:	×
453	return o	×
454
455	domain = extract_unique_domain(o)	×
NEW 456	tdim = domain.topological_dimension	×
457
458	if tdim == 1:	×
459	# Special-case 1D (possibly immersed), for which we say
460	# that n is just in the direction of J.
461	J = self.jacobian(Jacobian(domain)) # dx/dX	×
462	ndir = J[:, 0]	×
463
NEW 464	gdim = domain.geometric_dimension	×
465	if gdim == 1:	×
466	nlen = abs(ndir[0])	×
467	else:
468	i = Index()	×
469	nlen = sqrt(ndir[i] * ndir[i])	×
470
471	rn = ReferenceNormal(domain) # +/- 1.0 here	×
472	n = rn[0] * ndir / nlen	×
473	r = n	×
474	else:
475	# Recall that the covariant Piola transform u -> J^(-T)*u
476	# preserves tangential components. The normal vector is
477	# characterised by having zero tangential component in
478	# reference and physical space.
479	Jinv = self.jacobian_inverse(JacobianInverse(domain))	×
480	i, j = indices(2)	×
481
482	rn = ReferenceNormal(domain)	×
483	# compute signed, unnormalised normal; note transpose
484	ndir = as_vector(Jinv[j, i] * rn[j], i)	×
485
486	# normalise
487	i = Index()	×
488	n = ndir / sqrt(ndir[i] * ndir[i])	×
489	r = n	×
490
491	if r.ufl_shape != o.ufl_shape:	×
492	raise ValueError(	×
493	f"Inconsistent dimensions (in={o.ufl_shape[0]}, out={r.ufl_shape[0]})."
494	)
495	return r	×
496
497
498	def apply_geometry_lowering(form, preserve_types=()):	1✔
499	"""Change GeometricQuantity objects in expression to the lowest level GeometricQuantity objects.
500
501	Assumes the expression is preprocessed or at least that derivatives have been expanded.
502
503	Args:
504	form: An Expr or Form.
505	preserve_types: Preserved types
506	"""
507	if isinstance(form, Form):	1✔
508	newintegrals = [	1✔
509	apply_geometry_lowering(integral, preserve_types) for integral in form.integrals()
510	]
511	return Form(newintegrals)	1✔
512
513	elif isinstance(form, Integral):	1✔
514	integral = form	1✔
515	if integral.integral_type() in (custom_integral_types + point_integral_types):	1✔
516	automatic_preserve_types = [SpatialCoordinate, Jacobian]	×
517	else:
518	automatic_preserve_types = [CellCoordinate]	1✔
519	preserve_types = set(preserve_types) \| set(automatic_preserve_types)	1✔
520
521	mf = GeometryLoweringApplier(preserve_types)	1✔
522	newintegrand = map_expr_dag(mf, integral.integrand())	1✔
523	return integral.reconstruct(integrand=newintegrand)	1✔
524
525	elif isinstance(form, Expr):	1✔
526	expr = form	1✔
527	mf = GeometryLoweringApplier(preserve_types)	1✔
528	return map_expr_dag(mf, expr)	1✔
529
530	else:
531	raise ValueError(f"Invalid type {form.__class__.__name__}")	×

FEniCS / ufl / 17702620171

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