18630567039

Committed 19 Oct 2025 12:42PM UTC coverage: 75.964% (-0.7%) from 76.622%

Build # 18630567039

Build Type

Pull #431

github

Committed by

jorgensd

Commit Message

Test ufl against release branch

Pull Request Pull Request #431: Dokken/test release

Run Details

8960 of 11795 relevant lines covered (75.96%)

0.76 hits per line

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

80.81

/ufl/algorithms/compute_form_data.py

"""This module provides the compute_form_data function.

Form compilers will typically call compute_form_dataprior to code
generation to preprocess/simplify a raw input form given by a user.
"""
# 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 itertools import chain

from ufl.algorithms.analysis import extract_coefficients, extract_sub_elements, unique_tuple
from ufl.algorithms.apply_algebra_lowering import apply_algebra_lowering

# These are the main symbolic processing steps:
from ufl.algorithms.apply_coefficient_split import CoefficientSplitter
from ufl.algorithms.apply_derivatives import apply_coordinate_derivatives, apply_derivatives
from ufl.algorithms.apply_function_pullbacks import apply_function_pullbacks
from ufl.algorithms.apply_geometry_lowering import apply_geometry_lowering
from ufl.algorithms.apply_integral_scaling import apply_integral_scaling
from ufl.algorithms.apply_restrictions import apply_restrictions
from ufl.algorithms.check_arities import check_form_arity
from ufl.algorithms.comparison_checker import do_comparison_check

# See TODOs at the call sites of these below:
from ufl.algorithms.domain_analysis import (
    build_integral_data,
    group_form_integrals,
    reconstruct_form_from_integral_data,
)
from ufl.algorithms.estimate_degrees import estimate_total_polynomial_degree
from ufl.algorithms.formdata import FormData
from ufl.algorithms.formtransformations import compute_form_arities
from ufl.algorithms.remove_complex_nodes import remove_complex_nodes
from ufl.algorithms.remove_component_tensors import remove_component_tensors
from ufl.algorithms.replace import replace
from ufl.classes import Coefficient, Form, FunctionSpace, GeometricFacetQuantity
from ufl.constantvalue import Zero
from ufl.corealg.traversal import traverse_unique_terminals
from ufl.domain import MeshSequence, extract_domains, extract_unique_domain
from ufl.utils.sequences import max_degree


def _auto_select_degree(elements):
    """Automatically select degree for all elements of the form.

    This is be used in cases where the degree has not been specified by the user.
    This feature is used by DOLFIN to allow the specification of Expressions with
    undefined degrees.
    """
    # Use max degree of all elements, at least 1 (to work with
    # Lagrange elements)
    return max_degree({e.embedded_superdegree for e in elements} - {None} | {1})


def _compute_element_mapping(form):
    """Compute element mapping for element replacement."""
    # The element mapping is a slightly messy concept with two use
    # cases:
    # - Expression with missing cell or element TODO: Implement proper
    #   Expression handling in UFL and get rid of this
    # - Constant with missing cell TODO: Fix anything that needs to be
    #   worked around to drop this requirement

    # Extract all elements and include subelements of mixed elements
    elements = [obj.ufl_element() for obj in chain(form.arguments(), form.coefficients())]
    elements = extract_sub_elements(elements)

    # Try to find a common degree for elements
    common_degree = _auto_select_degree(elements)

    # Compute element map
    element_mapping = {}
    for element in elements:
        # Flag for whether element needs to be reconstructed
        reconstruct = False

        # Set cell
        cell = element.cell
        if cell is None:
            domains = form.ufl_domains()
            if not all(domains[0].ufl_cell() == d.ufl_cell() for d in domains):
                raise ValueError(
                    "Cannot replace unknown element cell without unique common cell in form."
                )
            cell = domains[0].ufl_cell()
            reconstruct = True

        # Set degree
        degree = element.embedded_superdegree
        if degree is None:
            degree = common_degree
            reconstruct = True

        # Reconstruct element and add to map
        if reconstruct:
            element_mapping[element] = element.reconstruct(cell=cell, degree=degree)
        else:
            element_mapping[element] = element

    return element_mapping


def _compute_max_subdomain_ids(integral_data):
    """Compute the maximum subdomain ids."""
    max_subdomain_ids = {}
    for itg_data in integral_data:
        it = itg_data.integral_type
        for integral in itg_data.integrals:
            # Convert string for default integral to -1
            sids = (-1 if isinstance(si, str) else si for si in integral.subdomain_id())
            newmax = max(sids) + 1
            prevmax = max_subdomain_ids.get(it, 0)
            max_subdomain_ids[it] = max(prevmax, newmax)
    return max_subdomain_ids


def _compute_form_data_elements(self, arguments, coefficients, domains):
    """Compute form data elements."""
    self.argument_elements = tuple(f.ufl_element() for f in arguments)
    self.coefficient_elements = tuple(f.ufl_element() for f in coefficients)
    self.coordinate_elements = tuple(domain.ufl_coordinate_element() for domain in domains)

    # TODO: Include coordinate elements from argument and coefficient
    # domains as well? Can they differ?

    # Note: Removed self.elements and self.sub_elements to make sure
    #       code that depends on the selection of argument +
    #       coefficient elements blow up, as opposed to silently
    #       almost working, with the introduction of the coordinate
    #       elements here.

    all_elements = self.argument_elements + self.coefficient_elements + self.coordinate_elements
    all_sub_elements = extract_sub_elements(all_elements)

    self.unique_elements = unique_tuple(all_elements)
    self.unique_sub_elements = unique_tuple(all_sub_elements)


def _check_elements(form_data):
    """Check elements."""
    for element in chain(form_data.unique_elements, form_data.unique_sub_elements):
        if element.cell is None:
            raise ValueError(f"Found element with undefined cell: {element}")


def _check_facet_geometry(integral_data):
    """Check facet geometry."""
    for itg_data in integral_data:
        for itg in itg_data.integrals:
            it = itg_data.integral_type
            # Facet geometry is only valid in facet integrals.
            # Allowing custom integrals to pass as well, although
            # that's not really strict enough.
            if not ("facet" in it or "custom" in it or "interface" in it):
                # Not a facet integral
                for expr in traverse_unique_terminals(itg.integrand()):
                    cls = expr._ufl_class_
                    if issubclass(cls, GeometricFacetQuantity):
                        raise ValueError(f"Integral of type {it} cannot contain a {cls.__name__}.")


def _check_form_arity(preprocessed_form):
    """Check that we don't have a mixed linear/bilinear form or anything like that."""
    # FIXME: This is slooow and should be moved to form compiler
    # and/or replaced with something faster
    if 1 != len(compute_form_arities(preprocessed_form)):
        raise ValueError("All terms in form must have same rank.")


def _build_coefficient_replace_map(coefficients, element_mapping=None):
    """Create new Coefficient objects with count starting at 0.

    Returns:
        mapping from old to new objects, and lists of the new objects
    """
    if element_mapping is None:
        element_mapping = {}

    new_coefficients = []
    replace_map = {}
    for i, f in enumerate(coefficients):
        old_e = f.ufl_element()
        new_e = element_mapping.get(old_e, old_e)
        # XXX: This is a hack to ensure that if the original
        # coefficient had a domain, the new one does too.
        # This should be overhauled with requirement that Expressions
        # always have a domain.
        domain = extract_unique_domain(f, expand_mesh_sequence=False)
        if domain is not None:
            new_e = FunctionSpace(domain, new_e)
        new_f = Coefficient(new_e, count=i)
        new_coefficients.append(new_f)
        replace_map[f] = new_f

    return new_coefficients, replace_map


def attach_estimated_degrees(form):
    """Attach estimated polynomial degree to a form's integrals.

    Args:
        form: The Form` to inspect.

    Returns:
        A new Form with estimate degrees attached.
    """
    integrals = form.integrals()

    new_integrals = []
    for integral in integrals:
        md = {}
        md.update(integral.metadata())
        degree = estimate_total_polynomial_degree(integral.integrand())
        md["estimated_polynomial_degree"] = degree
        new_integrals.append(integral.reconstruct(metadata=md))
    return Form(new_integrals)


def preprocess_form(form, complex_mode):
    """Preprocess a form."""
    # Note: Default behaviour here will process form the way that is
    # currently expected by vanilla FFC

    # Check that the form does not try to compare complex quantities:
    # if the quantites being compared are 'provably' real, wrap them
    # with Real, otherwise throw an error.
    if complex_mode:
        form = do_comparison_check(form)

    # Lower abstractions for tensor-algebra types into index notation,
    # reducing the number of operators later algorithms and form
    # compilers need to handle
    form = apply_algebra_lowering(form)

    # After lowering to index notation, remove any complex nodes that
    # have been introduced but are not wanted when working in real mode,
    # allowing for purely real forms to be written
    if not complex_mode:
        form = remove_complex_nodes(form)

    # Apply differentiation before function pullbacks, because for
    # example coefficient derivatives are more complicated to derive
    # after coefficients are rewritten, and in particular for
    # user-defined coefficient relations it just gets too messy
    form = apply_derivatives(form)

    return form


def compute_form_data(
    form,
    do_apply_function_pullbacks=False,
    do_apply_integral_scaling=False,
    do_apply_geometry_lowering=False,
    preserve_geometry_types=(),
    do_apply_default_restrictions=True,
    do_apply_restrictions=True,
    do_estimate_degrees=True,
    do_append_everywhere_integrals=True,
    do_replace_functions=False,
    coefficients_to_split=None,
    complex_mode=False,
    do_remove_component_tensors=False,
):
    """Compute form data.

    The default arguments configured to behave the way old FFC expects.
    """
    # Currently, only integral_type="cell" can be used with MeshSequence.
    for integral in form.integrals():
        if integral.integral_type() != "cell":
            all_domains = extract_domains(integral.integrand(), expand_mesh_sequence=False)
            if any(isinstance(m, MeshSequence) for m in all_domains):
                raise NotImplementedError(f"""
                    Only integral_type="cell" can be used with MeshSequence;
                    got integral_type={integral.integral_type()}
                """)

    # TODO: Move this to the constructor instead
    self = FormData()

    # --- Store untouched form for reference.
    # The user of FormData may get original arguments,
    # original coefficients, and form signature from this object.
    # But be aware that the set of original coefficients are not
    # the same as the ones used in the final UFC form.
    # See 'reduced_coefficients' below.
    self.original_form = form

    # --- Pass form integrands through some symbolic manipulation

    form = preprocess_form(form, complex_mode)

    # --- Group form integrals
    # TODO: Refactor this, it's rather opaque what this does
    # TODO: Is self.original_form.ufl_domains() right here?
    #       It will matter when we start including 'num_domains' in ufc form.
    form = group_form_integrals(
        form,
        self.original_form.ufl_domains(),
        do_append_everywhere_integrals=do_append_everywhere_integrals,
    )

    # Estimate polynomial degree of integrands now, before applying
    # any pullbacks and geometric lowering.  Otherwise quad degrees
    # blow up horrifically.
    if do_estimate_degrees:
        form = attach_estimated_degrees(form)

    if do_apply_function_pullbacks:
        # Rewrite coefficients and arguments in terms of their
        # reference cell values with Piola transforms and symmetry
        # transforms injected where needed.
        # Decision: Not supporting grad(dolfin.Expression) without a
        #           Domain.  Current dolfin works if Expression has a
        #           cell but this should be changed to a mesh.
        form = apply_function_pullbacks(form)

    # Scale integrals to reference cell frames
    if do_apply_integral_scaling:
        form = apply_integral_scaling(form)

    # Lower abstractions for geometric quantities into a smaller set
    # of quantities, allowing the form compiler to deal with a smaller
    # set of types and treating geometric quantities like any other
    # expressions w.r.t. loop-invariant code motion etc.
    if do_apply_geometry_lowering:
        form = apply_geometry_lowering(form, preserve_geometry_types)

    # Apply differentiation again, because the algorithms above can
    # generate new derivatives or rewrite expressions inside
    # derivatives
    if do_apply_function_pullbacks or do_apply_geometry_lowering:
        form = apply_derivatives(form)

        # Neverending story: apply_derivatives introduces new Jinvs,
        # which needs more geometry lowering
        if do_apply_geometry_lowering:
            form = apply_geometry_lowering(form, preserve_geometry_types)
            # Lower derivatives that may have appeared
            form = apply_derivatives(form)

    form = apply_coordinate_derivatives(form)

    # Propagate restrictions to terminals
    if do_apply_restrictions:
        form = apply_restrictions(form, apply_default=do_apply_default_restrictions)

    # If in real mode, remove any complex nodes introduced during form processing.
    if not complex_mode:
        form = remove_complex_nodes(form)

    # Remove component tensors
    if do_remove_component_tensors:
        form = remove_component_tensors(form)

    # --- Group integrals into IntegralData objects
    # Most of the heavy lifting is done above in group_form_integrals.
    self.integral_data = build_integral_data(form.integrals())

    # --- Create replacements for arguments and coefficients

    # Figure out which form coefficients each integral should enable
    for itg_data in self.integral_data:
        itg_coeffs = set()
        # Get all coefficients in integrand
        for itg in itg_data.integrals:
            itg_coeffs.update(extract_coefficients(itg.integrand()))
        # Store with IntegralData object
        itg_data.integral_coefficients = itg_coeffs

    # Figure out which coefficients from the original form are
    # actually used in any integral (Differentiation may reduce the
    # set of coefficients w.r.t. the original form)
    reduced_coefficients_set = set()
    for itg_data in self.integral_data:
        reduced_coefficients_set.update(itg_data.integral_coefficients)
    self.reduced_coefficients = sorted(reduced_coefficients_set, key=lambda c: c.count())
    self.num_coefficients = len(self.reduced_coefficients)
    self.original_coefficient_positions = [
        i for i, c in enumerate(self.original_form.coefficients()) if c in self.reduced_coefficients
    ]

    # Store back into integral data which form coefficients are used
    # by each integral
    for itg_data in self.integral_data:
        itg_data.enabled_coefficients = [
            bool(coeff in itg_data.integral_coefficients) for coeff in self.reduced_coefficients
        ]

    # --- Collect some trivial data

    # Get rank of form from argument list (assuming not a mixed arity form)
    self.rank = len(self.original_form.arguments())

    # Extract common geometric dimension (topological is not common!)
    self.geometric_dimension = self.original_form.integrals()[0].ufl_domain().geometric_dimension()

    # --- Build mapping from old incomplete element objects to new
    # well defined elements.  This is to support the Expression
    # construct in dolfin which subclasses Coefficient but doesn't
    # provide an element, and the Constant construct that doesn't
    # provide the domain that a Coefficient is supposed to have. A
    # future design iteration in UFL/UFC/FFC/DOLFIN may allow removal
    # of this mapping with the introduction of UFL types for
    # Expression-like functions that can be evaluated in quadrature
    # points.
    self.element_replace_map = _compute_element_mapping(self.original_form)

    # Mappings from elements and coefficients that reside in form to
    # objects with canonical numbering as well as completed cells and
    # elements
    renumbered_coefficients, function_replace_map = _build_coefficient_replace_map(
        self.reduced_coefficients, self.element_replace_map
    )
    self.function_replace_map = function_replace_map

    # --- Store various lists of elements and sub elements (adds
    #     members to self)
    _compute_form_data_elements(
        self,
        self.original_form.arguments(),
        renumbered_coefficients,
        self.original_form.ufl_domains(),
    )

    # --- Store number of domains for integral types
    # TODO: Group this by domain first. For now keep a backwards
    # compatible data structure.
    self.max_subdomain_ids = _compute_max_subdomain_ids(self.integral_data)

    # --- Apply replace(integrand, self.function_replace_map)
    if do_replace_functions:
        for itg_data in self.integral_data:
            new_integrals = []
            for integral in itg_data.integrals:
                integrand = replace(integral.integrand(), self.function_replace_map)
                new_integrals.append(integral.reconstruct(integrand=integrand))
            itg_data.integrals = new_integrals

    # --- Split mixed coefficients with their components
    if coefficients_to_split is None:
        self.coefficient_split = {}
    else:
        if not do_replace_functions:
            raise ValueError("Must call with do_replace_functions=True")
        # Split coefficients that are contained in ``coefficients_to_split``
        # into components, and store a dict in ``self`` that maps
        # each coefficient to its components.
        coefficient_split = {}
        for o in self.reduced_coefficients:
            if o in coefficients_to_split:
                c = self.function_replace_map[o]
                mesh = extract_unique_domain(c, expand_mesh_sequence=False)
                elem = c.ufl_element()
                coefficient_split[c] = [
                    Coefficient(FunctionSpace(m, e))
                    for m, e in zip(mesh.iterable_like(elem), elem.sub_elements)
                ]
        self.coefficient_split = coefficient_split
        coeff_splitter = CoefficientSplitter(self.coefficient_split)
        for itg_data in self.integral_data:
            new_integrals = []
            for integral in itg_data.integrals:
                integrand = coeff_splitter(integral.integrand())
                if not isinstance(integrand, Zero):
                    new_integrals.append(integral.reconstruct(integrand=integrand))
            itg_data.integrals = new_integrals

    # --- Checks
    _check_elements(self)
    _check_facet_geometry(self.integral_data)

    # TODO: This is a very expensive check... Replace with something
    # faster!
    preprocessed_form = reconstruct_form_from_integral_data(self.integral_data)

    # TODO: Test how fast this is
    check_form_arity(preprocessed_form, self.original_form.arguments(), complex_mode)

    # TODO: This member is used by unit tests, change the tests to
    # remove this!
    self.preprocessed_form = preprocessed_form

    return self

1	"""This module provides the compute_form_data function.
2
3	Form compilers will typically call compute_form_dataprior to code
4	generation to preprocess/simplify a raw input form given by a user.
5	"""
6	# Copyright (C) 2008-2016 Martin Sandve Alnæs
7	#
8	# This file is part of UFL (https://www.fenicsproject.org)
9	#
10	# SPDX-License-Identifier: LGPL-3.0-or-later
11
12	from itertools import chain	1✔
13
14	from ufl.algorithms.analysis import extract_coefficients, extract_sub_elements, unique_tuple	1✔
15	from ufl.algorithms.apply_algebra_lowering import apply_algebra_lowering	1✔
16
17	# These are the main symbolic processing steps:
18	from ufl.algorithms.apply_coefficient_split import CoefficientSplitter	1✔
19	from ufl.algorithms.apply_derivatives import apply_coordinate_derivatives, apply_derivatives	1✔
20	from ufl.algorithms.apply_function_pullbacks import apply_function_pullbacks	1✔
21	from ufl.algorithms.apply_geometry_lowering import apply_geometry_lowering	1✔
22	from ufl.algorithms.apply_integral_scaling import apply_integral_scaling	1✔
23	from ufl.algorithms.apply_restrictions import apply_restrictions	1✔
24	from ufl.algorithms.check_arities import check_form_arity	1✔
25	from ufl.algorithms.comparison_checker import do_comparison_check	1✔
26
27	# See TODOs at the call sites of these below:
28	from ufl.algorithms.domain_analysis import (	1✔
29	build_integral_data,
30	group_form_integrals,
31	reconstruct_form_from_integral_data,
32	)
33	from ufl.algorithms.estimate_degrees import estimate_total_polynomial_degree	1✔
34	from ufl.algorithms.formdata import FormData	1✔
35	from ufl.algorithms.formtransformations import compute_form_arities	1✔
36	from ufl.algorithms.remove_complex_nodes import remove_complex_nodes	1✔
37	from ufl.algorithms.remove_component_tensors import remove_component_tensors	1✔
38	from ufl.algorithms.replace import replace	1✔
39	from ufl.classes import Coefficient, Form, FunctionSpace, GeometricFacetQuantity	1✔
40	from ufl.constantvalue import Zero	1✔
41	from ufl.corealg.traversal import traverse_unique_terminals	1✔
42	from ufl.domain import MeshSequence, extract_domains, extract_unique_domain	1✔
43	from ufl.utils.sequences import max_degree	1✔
44
45
46	def _auto_select_degree(elements):	1✔
47	"""Automatically select degree for all elements of the form.
48
49	This is be used in cases where the degree has not been specified by the user.
50	This feature is used by DOLFIN to allow the specification of Expressions with
51	undefined degrees.
52	"""
53	# Use max degree of all elements, at least 1 (to work with
54	# Lagrange elements)
55	return max_degree({e.embedded_superdegree for e in elements} - {None} \| {1})	1✔
56
57
58	def _compute_element_mapping(form):	1✔
59	"""Compute element mapping for element replacement."""
60	# The element mapping is a slightly messy concept with two use
61	# cases:
62	# - Expression with missing cell or element TODO: Implement proper
63	# Expression handling in UFL and get rid of this
64	# - Constant with missing cell TODO: Fix anything that needs to be
65	# worked around to drop this requirement
66
67	# Extract all elements and include subelements of mixed elements
68	elements = [obj.ufl_element() for obj in chain(form.arguments(), form.coefficients())]	1✔
69	elements = extract_sub_elements(elements)	1✔
70
71	# Try to find a common degree for elements
72	common_degree = _auto_select_degree(elements)	1✔
73
74	# Compute element map
75	element_mapping = {}	1✔
76	for element in elements:	1✔
77	# Flag for whether element needs to be reconstructed
78	reconstruct = False	1✔
79
80	# Set cell
81	cell = element.cell	1✔
82	if cell is None:	1✔
83	domains = form.ufl_domains()	×
84	if not all(domains[0].ufl_cell() == d.ufl_cell() for d in domains):	×
85	raise ValueError(	×
86	"Cannot replace unknown element cell without unique common cell in form."
87	)
88	cell = domains[0].ufl_cell()	×
89	reconstruct = True	×
90
91	# Set degree
92	degree = element.embedded_superdegree	1✔
93	if degree is None:	1✔
94	degree = common_degree	×
95	reconstruct = True	×
96
97	# Reconstruct element and add to map
98	if reconstruct:	1✔
99	element_mapping[element] = element.reconstruct(cell=cell, degree=degree)	×
100	else:
101	element_mapping[element] = element	1✔
102
103	return element_mapping	1✔
104
105
106	def _compute_max_subdomain_ids(integral_data):	1✔
107	"""Compute the maximum subdomain ids."""
108	max_subdomain_ids = {}	1✔
109	for itg_data in integral_data:	1✔
110	it = itg_data.integral_type	1✔
111	for integral in itg_data.integrals:	1✔
112	# Convert string for default integral to -1
113	sids = (-1 if isinstance(si, str) else si for si in integral.subdomain_id())	1✔
114	newmax = max(sids) + 1	1✔
115	prevmax = max_subdomain_ids.get(it, 0)	1✔
116	max_subdomain_ids[it] = max(prevmax, newmax)	1✔
117	return max_subdomain_ids	1✔
118
119
120	def _compute_form_data_elements(self, arguments, coefficients, domains):	1✔
121	"""Compute form data elements."""
122	self.argument_elements = tuple(f.ufl_element() for f in arguments)	1✔
123	self.coefficient_elements = tuple(f.ufl_element() for f in coefficients)	1✔
124	self.coordinate_elements = tuple(domain.ufl_coordinate_element() for domain in domains)	1✔
125
126	# TODO: Include coordinate elements from argument and coefficient
127	# domains as well? Can they differ?
128
129	# Note: Removed self.elements and self.sub_elements to make sure
130	# code that depends on the selection of argument +
131	# coefficient elements blow up, as opposed to silently
132	# almost working, with the introduction of the coordinate
133	# elements here.
134
135	all_elements = self.argument_elements + self.coefficient_elements + self.coordinate_elements	1✔
136	all_sub_elements = extract_sub_elements(all_elements)	1✔
137
138	self.unique_elements = unique_tuple(all_elements)	1✔
139	self.unique_sub_elements = unique_tuple(all_sub_elements)	1✔
140
141
142	def _check_elements(form_data):	1✔
143	"""Check elements."""
144	for element in chain(form_data.unique_elements, form_data.unique_sub_elements):	1✔
145	if element.cell is None:	1✔
146	raise ValueError(f"Found element with undefined cell: {element}")	×
147
148
149	def _check_facet_geometry(integral_data):	1✔
150	"""Check facet geometry."""
151	for itg_data in integral_data:	1✔
152	for itg in itg_data.integrals:	1✔
153	it = itg_data.integral_type	1✔
154	# Facet geometry is only valid in facet integrals.
155	# Allowing custom integrals to pass as well, although
156	# that's not really strict enough.
157	if not ("facet" in it or "custom" in it or "interface" in it):	1✔
158	# Not a facet integral
159	for expr in traverse_unique_terminals(itg.integrand()):	1✔
160	cls = expr._ufl_class_	1✔
161	if issubclass(cls, GeometricFacetQuantity):	1✔
162	raise ValueError(f"Integral of type {it} cannot contain a {cls.__name__}.")	×
163
164
165	def _check_form_arity(preprocessed_form):	1✔
166	"""Check that we don't have a mixed linear/bilinear form or anything like that."""
167	# FIXME: This is slooow and should be moved to form compiler
168	# and/or replaced with something faster
169	if 1 != len(compute_form_arities(preprocessed_form)):	×
170	raise ValueError("All terms in form must have same rank.")	×
171
172
173	def _build_coefficient_replace_map(coefficients, element_mapping=None):	1✔
174	"""Create new Coefficient objects with count starting at 0.
175
176	Returns:
177	mapping from old to new objects, and lists of the new objects
178	"""
179	if element_mapping is None:	1✔
180	element_mapping = {}	×
181
182	new_coefficients = []	1✔
183	replace_map = {}	1✔
184	for i, f in enumerate(coefficients):	1✔
185	old_e = f.ufl_element()	1✔
186	new_e = element_mapping.get(old_e, old_e)	1✔
187	# XXX: This is a hack to ensure that if the original
188	# coefficient had a domain, the new one does too.
189	# This should be overhauled with requirement that Expressions
190	# always have a domain.
191	domain = extract_unique_domain(f, expand_mesh_sequence=False)	1✔
192	if domain is not None:	1✔
193	new_e = FunctionSpace(domain, new_e)	1✔
194	new_f = Coefficient(new_e, count=i)	1✔
195	new_coefficients.append(new_f)	1✔
196	replace_map[f] = new_f	1✔
197
198	return new_coefficients, replace_map	1✔
199
200
201	def attach_estimated_degrees(form):	1✔
202	"""Attach estimated polynomial degree to a form's integrals.
203
204	Args:
205	form: The Form` to inspect.
206
207	Returns:
208	A new Form with estimate degrees attached.
209	"""
210	integrals = form.integrals()	1✔
211
212	new_integrals = []	1✔
213	for integral in integrals:	1✔
214	md = {}	1✔
215	md.update(integral.metadata())	1✔
216	degree = estimate_total_polynomial_degree(integral.integrand())	1✔
217	md["estimated_polynomial_degree"] = degree	1✔
218	new_integrals.append(integral.reconstruct(metadata=md))	1✔
219	return Form(new_integrals)	1✔
220
221
222	def preprocess_form(form, complex_mode):	1✔
223	"""Preprocess a form."""
224	# Note: Default behaviour here will process form the way that is
225	# currently expected by vanilla FFC
226
227	# Check that the form does not try to compare complex quantities:
228	# if the quantites being compared are 'provably' real, wrap them
229	# with Real, otherwise throw an error.
230	if complex_mode:	1✔
231	form = do_comparison_check(form)	1✔
232
233	# Lower abstractions for tensor-algebra types into index notation,
234	# reducing the number of operators later algorithms and form
235	# compilers need to handle
236	form = apply_algebra_lowering(form)	1✔
237
238	# After lowering to index notation, remove any complex nodes that
239	# have been introduced but are not wanted when working in real mode,
240	# allowing for purely real forms to be written
241	if not complex_mode:	1✔
242	form = remove_complex_nodes(form)	1✔
243
244	# Apply differentiation before function pullbacks, because for
245	# example coefficient derivatives are more complicated to derive
246	# after coefficients are rewritten, and in particular for
247	# user-defined coefficient relations it just gets too messy
248	form = apply_derivatives(form)	1✔
249
250	return form	1✔
251
252
253	def compute_form_data(	1✔
254	form,
255	do_apply_function_pullbacks=False,
256	do_apply_integral_scaling=False,
257	do_apply_geometry_lowering=False,
258	preserve_geometry_types=(),
259	do_apply_default_restrictions=True,
260	do_apply_restrictions=True,
261	do_estimate_degrees=True,
262	do_append_everywhere_integrals=True,
263	do_replace_functions=False,
264	coefficients_to_split=None,
265	complex_mode=False,
266	do_remove_component_tensors=False,
267	):
268	"""Compute form data.
269
270	The default arguments configured to behave the way old FFC expects.
271	"""
272	# Currently, only integral_type="cell" can be used with MeshSequence.
273	for integral in form.integrals():	1✔
274	if integral.integral_type() != "cell":	1✔
275	all_domains = extract_domains(integral.integrand(), expand_mesh_sequence=False)	1✔
276	if any(isinstance(m, MeshSequence) for m in all_domains):	1✔
277	raise NotImplementedError(f"""	×
278	Only integral_type="cell" can be used with MeshSequence;
279	got integral_type={integral.integral_type()}
280	""")
281
282	# TODO: Move this to the constructor instead
283	self = FormData()	1✔
284
285	# --- Store untouched form for reference.
286	# The user of FormData may get original arguments,
287	# original coefficients, and form signature from this object.
288	# But be aware that the set of original coefficients are not
289	# the same as the ones used in the final UFC form.
290	# See 'reduced_coefficients' below.
291	self.original_form = form	1✔
292
293	# --- Pass form integrands through some symbolic manipulation
294
295	form = preprocess_form(form, complex_mode)	1✔
296
297	# --- Group form integrals
298	# TODO: Refactor this, it's rather opaque what this does
299	# TODO: Is self.original_form.ufl_domains() right here?
300	# It will matter when we start including 'num_domains' in ufc form.
301	form = group_form_integrals(	1✔
302	form,
303	self.original_form.ufl_domains(),
304	do_append_everywhere_integrals=do_append_everywhere_integrals,
305	)
306
307	# Estimate polynomial degree of integrands now, before applying
308	# any pullbacks and geometric lowering. Otherwise quad degrees
309	# blow up horrifically.
310	if do_estimate_degrees:	1✔
311	form = attach_estimated_degrees(form)	1✔
312
313	if do_apply_function_pullbacks:	1✔
314	# Rewrite coefficients and arguments in terms of their
315	# reference cell values with Piola transforms and symmetry
316	# transforms injected where needed.
317	# Decision: Not supporting grad(dolfin.Expression) without a
318	# Domain. Current dolfin works if Expression has a
319	# cell but this should be changed to a mesh.
320	form = apply_function_pullbacks(form)	1✔
321
322	# Scale integrals to reference cell frames
323	if do_apply_integral_scaling:	1✔
324	form = apply_integral_scaling(form)	1✔
325
326	# Lower abstractions for geometric quantities into a smaller set
327	# of quantities, allowing the form compiler to deal with a smaller
328	# set of types and treating geometric quantities like any other
329	# expressions w.r.t. loop-invariant code motion etc.
330	if do_apply_geometry_lowering:	1✔
331	form = apply_geometry_lowering(form, preserve_geometry_types)	1✔
332
333	# Apply differentiation again, because the algorithms above can
334	# generate new derivatives or rewrite expressions inside
335	# derivatives
336	if do_apply_function_pullbacks or do_apply_geometry_lowering:	1✔
337	form = apply_derivatives(form)	1✔
338
339	# Neverending story: apply_derivatives introduces new Jinvs,
340	# which needs more geometry lowering
341	if do_apply_geometry_lowering:	1✔
342	form = apply_geometry_lowering(form, preserve_geometry_types)	1✔
343	# Lower derivatives that may have appeared
344	form = apply_derivatives(form)	1✔
345
346	form = apply_coordinate_derivatives(form)	1✔
347
348	# Propagate restrictions to terminals
349	if do_apply_restrictions:	1✔
350	form = apply_restrictions(form, apply_default=do_apply_default_restrictions)	1✔
351
352	# If in real mode, remove any complex nodes introduced during form processing.
353	if not complex_mode:	1✔
354	form = remove_complex_nodes(form)	1✔
355
356	# Remove component tensors
357	if do_remove_component_tensors:	1✔
358	form = remove_component_tensors(form)	1✔
359
360	# --- Group integrals into IntegralData objects
361	# Most of the heavy lifting is done above in group_form_integrals.
362	self.integral_data = build_integral_data(form.integrals())	1✔
363
364	# --- Create replacements for arguments and coefficients
365
366	# Figure out which form coefficients each integral should enable
367	for itg_data in self.integral_data:	1✔
368	itg_coeffs = set()	1✔
369	# Get all coefficients in integrand
370	for itg in itg_data.integrals:	1✔
371	itg_coeffs.update(extract_coefficients(itg.integrand()))	1✔
372	# Store with IntegralData object
373	itg_data.integral_coefficients = itg_coeffs	1✔
374
375	# Figure out which coefficients from the original form are
376	# actually used in any integral (Differentiation may reduce the
377	# set of coefficients w.r.t. the original form)
378	reduced_coefficients_set = set()	1✔
379	for itg_data in self.integral_data:	1✔
380	reduced_coefficients_set.update(itg_data.integral_coefficients)	1✔
381	self.reduced_coefficients = sorted(reduced_coefficients_set, key=lambda c: c.count())	1✔
382	self.num_coefficients = len(self.reduced_coefficients)	1✔
383	self.original_coefficient_positions = [	1✔
384	i for i, c in enumerate(self.original_form.coefficients()) if c in self.reduced_coefficients
385	]
386
387	# Store back into integral data which form coefficients are used
388	# by each integral
389	for itg_data in self.integral_data:	1✔
390	itg_data.enabled_coefficients = [	1✔
391	bool(coeff in itg_data.integral_coefficients) for coeff in self.reduced_coefficients
392	]
393
394	# --- Collect some trivial data
395
396	# Get rank of form from argument list (assuming not a mixed arity form)
397	self.rank = len(self.original_form.arguments())	1✔
398
399	# Extract common geometric dimension (topological is not common!)
400	self.geometric_dimension = self.original_form.integrals()[0].ufl_domain().geometric_dimension()	1✔
401
402	# --- Build mapping from old incomplete element objects to new
403	# well defined elements. This is to support the Expression
404	# construct in dolfin which subclasses Coefficient but doesn't
405	# provide an element, and the Constant construct that doesn't
406	# provide the domain that a Coefficient is supposed to have. A
407	# future design iteration in UFL/UFC/FFC/DOLFIN may allow removal
408	# of this mapping with the introduction of UFL types for
409	# Expression-like functions that can be evaluated in quadrature
410	# points.
411	self.element_replace_map = _compute_element_mapping(self.original_form)	1✔
412
413	# Mappings from elements and coefficients that reside in form to
414	# objects with canonical numbering as well as completed cells and
415	# elements
416	renumbered_coefficients, function_replace_map = _build_coefficient_replace_map(	1✔
417	self.reduced_coefficients, self.element_replace_map
418	)
419	self.function_replace_map = function_replace_map	1✔
420
421	# --- Store various lists of elements and sub elements (adds
422	# members to self)
423	_compute_form_data_elements(	1✔
424	self,
425	self.original_form.arguments(),
426	renumbered_coefficients,
427	self.original_form.ufl_domains(),
428	)
429
430	# --- Store number of domains for integral types
431	# TODO: Group this by domain first. For now keep a backwards
432	# compatible data structure.
433	self.max_subdomain_ids = _compute_max_subdomain_ids(self.integral_data)	1✔
434
435	# --- Apply replace(integrand, self.function_replace_map)
436	if do_replace_functions:	1✔
437	for itg_data in self.integral_data:	×
438	new_integrals = []	×
439	for integral in itg_data.integrals:	×
440	integrand = replace(integral.integrand(), self.function_replace_map)	×
441	new_integrals.append(integral.reconstruct(integrand=integrand))	×
442	itg_data.integrals = new_integrals	×
443
444	# --- Split mixed coefficients with their components
445	if coefficients_to_split is None:	1✔
446	self.coefficient_split = {}	1✔
447	else:
448	if not do_replace_functions:	×
449	raise ValueError("Must call with do_replace_functions=True")	×
450	# Split coefficients that are contained in ``coefficients_to_split``
451	# into components, and store a dict in ``self`` that maps
452	# each coefficient to its components.
453	coefficient_split = {}	×
454	for o in self.reduced_coefficients:	×
455	if o in coefficients_to_split:	×
456	c = self.function_replace_map[o]	×
457	mesh = extract_unique_domain(c, expand_mesh_sequence=False)	×
458	elem = c.ufl_element()	×
459	coefficient_split[c] = [	×
460	Coefficient(FunctionSpace(m, e))
461	for m, e in zip(mesh.iterable_like(elem), elem.sub_elements)
462	]
463	self.coefficient_split = coefficient_split	×
464	coeff_splitter = CoefficientSplitter(self.coefficient_split)	×
465	for itg_data in self.integral_data:	×
466	new_integrals = []	×
467	for integral in itg_data.integrals:	×
468	integrand = coeff_splitter(integral.integrand())	×
469	if not isinstance(integrand, Zero):	×
470	new_integrals.append(integral.reconstruct(integrand=integrand))	×
471	itg_data.integrals = new_integrals	×
472
473	# --- Checks
474	_check_elements(self)	1✔
475	_check_facet_geometry(self.integral_data)	1✔
476
477	# TODO: This is a very expensive check... Replace with something
478	# faster!
479	preprocessed_form = reconstruct_form_from_integral_data(self.integral_data)	1✔
480
481	# TODO: Test how fast this is
482	check_form_arity(preprocessed_form, self.original_form.arguments(), complex_mode)	1✔
483
484	# TODO: This member is used by unit tests, change the tests to
485	# remove this!
486	self.preprocessed_form = preprocessed_form	1✔
487
488	return self	1✔

FEniCS / ufl / 18630567039

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