diff --git a/hog/forms.py b/hog/forms.py index 439e19a7cd59790807e94ad434c91dc896406d6d..28b030a61de8590d81d0d91ab78391e0487d7660 100644 --- a/hog/forms.py +++ b/hog/forms.py @@ -71,8 +71,15 @@ Weak formulation integrand as integrand_affine, ) + # mypy type checking supposedly cannot figure out ternaries. + # https://stackoverflow.com/a/70832391 + if blending == IdentityMap(): + integr = integrand_affine + else: + integr = integrand + return process_integrand( - (integrand_affine if blending == IdentityMap() else integrand), + integr, trial, test, geometry, @@ -144,8 +151,15 @@ Weak formulation integrand as integrand_affine, ) + # mypy type checking supposedly cannot figure out ternaries. + # https://stackoverflow.com/a/70832391 + if blending == IdentityMap(): + integr = integrand_affine + else: + integr = integrand + return process_integrand( - integrand_affine if blending == IdentityMap() else integrand, + integr, trial, test, geometry, @@ -336,8 +350,15 @@ where integrand as integrand_affine, ) + # mypy type checking supposedly cannot figure out ternaries. + # https://stackoverflow.com/a/70832391 + if blending == IdentityMap(): + integr = integrand_affine + else: + integr = integrand + return process_integrand( - (integrand_affine if blending == IdentityMap() else integrand), + integr, trial, test, geometry, diff --git a/hog/forms_vectorial.py b/hog/forms_vectorial.py index 237b3619db6ce8e13bba2d0084e8ca581953391c..a49ce7058925404dea4270ff20986d3ebf162c8a 100644 --- a/hog/forms_vectorial.py +++ b/hog/forms_vectorial.py @@ -30,7 +30,7 @@ from hog.fem_helpers import ( element_matrix_iterator, scalar_space_dependent_coefficient, ) -from hog.forms import Form +from hog.integrand import Form from hog.function_space import FunctionSpace, EnrichedGalerkinFunctionSpace, N1E1Space from hog.math_helpers import inv, abs, det, double_contraction, dot, curl from hog.quadrature import Quadrature, Tabulation diff --git a/hog/function_space.py b/hog/function_space.py index b240a53cbdb1f5dc5f65bdde5b98742836b21328..0f4029c13b287d177d4177679520cc773669583c 100644 --- a/hog/function_space.py +++ b/hog/function_space.py @@ -370,7 +370,9 @@ class TensorialVectorFunctionSpace(FunctionSpace): """ - def __init__(self, function_space: FunctionSpace, single_component: int = None): + def __init__( + self, function_space: FunctionSpace, single_component: Union[None, int] = None + ): """ Initializes a tensorial vector function space from a scalar function space. diff --git a/hog/integrand.py b/hog/integrand.py index 875d38bcaf4ea4de701726fe61b2b57cc111030a..e45907ecb4f151d012fb6b36e6c0571e9394c05d 100644 --- a/hog/integrand.py +++ b/hog/integrand.py @@ -46,7 +46,7 @@ This module contains various functions and data structures to formulate the inte The actual integration and code generation is handled somewhere else. """ -from typing import Callable, List, Union, Tuple +from typing import Callable, List, Union, Tuple, Any from dataclasses import dataclass, asdict, fields import sympy as sp @@ -114,66 +114,68 @@ class IntegrandSymbols: return [f.name for f in fields(IntegrandSymbols) if not f.name.startswith("_")] # Jacobian from reference to affine space. - jac_a: sp.Matrix = None + jac_a: sp.Matrix | None = None # Its inverse. - jac_a_inv: sp.Matrix = None + jac_a_inv: sp.Matrix | None = None # The absolute of its determinant. - jac_a_abs_det: sp.Symbol = None + jac_a_abs_det: sp.Symbol | None = None # Jacobian from affine to physical space. - jac_b: sp.Matrix = None + jac_b: sp.Matrix | None = None # Its inverse. - jac_b_inv: sp.Matrix = None + jac_b_inv: sp.Matrix | None = None # The absolute of its determinant. - jac_b_abs_det: sp.Symbol = None + jac_b_abs_det: sp.Symbol | None = None # The trial shape function (reference space!). - u: sp.Expr = None + u: sp.Expr | None = None # The gradient of the trial shape function (reference space!). - grad_u: sp.Matrix = None + grad_u: sp.Matrix | None = None # The test shape function (reference space!). - v: sp.Expr = None + v: sp.Expr | None = None # The gradient of the test shape function (reference space!). - grad_v: sp.Matrix = None + grad_v: sp.Matrix | None = None # A list of scalar constants. - c: List[sp.Symbol] = None + c: List[sp.Symbol] | None = None # A list of finite element functions that can be used as function parameters. - k: List[sp.Symbol] = None + k: List[sp.Symbol] | None = None + # A list of the gradients of the parameter finite element functions. + grad_k: List[sp.Matrix] | None = None # The geometry of the volume element. - volume_geometry: ElementGeometry = None + volume_geometry: ElementGeometry | None = None # The geometry of the boundary element. - boundary_geometry: ElementGeometry = None + boundary_geometry: ElementGeometry | None = None # If a boundary geometry is available, this is populated with the Jacobian of the affine mapping from the reference # space of the boundary element to the computational (affine) space. # The reference space has the dimensions of the boundary element. # The affine space has the space dimension (aka the dimension of the space it is embedded in) of the boundary # element. - jac_a_boundary: sp.Matrix = None + jac_a_boundary: sp.Matrix | None = None # A callback to tabulate (aka precompute) terms that are identical on all elements of the same type. # Use at your own risk, you may get wrong code if used on terms that are not element-invariant! - tabulate: Callable = None + tabulate: Union[Callable, None] = None # You can also give the tabulated variable a name. That has no effect other than the generated code to be more # readable. So not encouraged. But nice for debugging. - _tabulate_named: Callable = None + _tabulate_named: Union[Callable, None] = None def process_integrand( - integrand: Callable, + integrand: Callable[..., Any], trial: FunctionSpace, test: FunctionSpace, volume_geometry: ElementGeometry, symbolizer: Symbolizer, blending: GeometryMap = IdentityMap(), - boundary_geometry: ElementGeometry = None, - scalar_coefficients: List[str] = None, - fe_coefficients: List[Tuple[str, Union[FunctionSpace, None]]] = None, + boundary_geometry: ElementGeometry | None = None, + scalar_coefficients: List[str] | None = None, + fe_coefficients: List[Tuple[str, Union[FunctionSpace, None]]] | None = None, is_symmetric: bool = False, docstring: str = "", ) -> Form: @@ -186,32 +188,32 @@ def process_integrand( Integrands are passed in as a callable (aka function). For instance: - ```python - # The arguments of the function must begin with an asterisk (*), followed by keyword arguments, followed by the - # unused keyword arguments (**_). All keyword arguments must be members of the IntegrandSymbols class. - # The function must return the integrand. You can use functions from the module hog.math_helpers module. - # Many integrands are already implemented under hog/recipes/integrands/. - - def my_diffusion_integrand( - *, - jac_a_inv, - jac_a_abs_det, - jac_b_inv, - jac_b_abs_det, - grad_u, - grad_v, - tabulate, - **_, - ): - return ( - double_contraction( - jac_b_inv.T * tabulate(jac_a_inv.T * grad_u), - jac_b_inv.T * tabulate(jac_a_inv.T * grad_v), + .. code-block:: python + + # The arguments of the function must begin with an asterisk (*), followed by keyword arguments, followed by the + # unused keyword arguments (**_). All keyword arguments must be members of the IntegrandSymbols class. + # The function must return the integrand. You can use functions from the module hog.math_helpers module. + # Many integrands are already implemented under hog/recipes/integrands/. + + def my_diffusion_integrand( + *, + jac_a_inv, + jac_a_abs_det, + jac_b_inv, + jac_b_abs_det, + grad_u, + grad_v, + tabulate, + **_, + ): + return ( + double_contraction( + jac_b_inv.T * tabulate(jac_a_inv.T * grad_u), + jac_b_inv.T * tabulate(jac_a_inv.T * grad_v), + ) + * tabulate(jac_a_abs_det) + * jac_b_abs_det ) - * tabulate(jac_a_abs_det) - * jac_b_abs_det - ) - ``` The callable (here `my_diffusion_integrand`, not `my_diffusion_integrand()`) is then passed to this function. @@ -242,7 +244,7 @@ def process_integrand( tabulation = Tabulation(symbolizer) - def _tabulate(factor: Union[sp.Expr, sp.Matrix]): + def _tabulate(factor: Union[sp.Expr, sp.Matrix]) -> sp.Matrix: if isinstance(factor, sp.Expr): factor = sp.Matrix([factor]) @@ -250,7 +252,7 @@ def process_integrand( f"tabulated_factor_{symbolizer.get_next_running_integer()}", factor ) - def _tabulate_named(factor_name: str, factor: sp.Matrix): + def _tabulate_named(factor_name: str, factor: sp.Matrix) -> sp.Matrix: return tabulation.register_factor(factor_name, factor) s.tabulate = _tabulate diff --git a/hog/manifold_forms.py b/hog/manifold_forms.py index 81f92c5b6941735829caab8151b67bbf980ad2c8..d6ccfc20a39d9b95932d4ce03e77aeddfc96ece0 100644 --- a/hog/manifold_forms.py +++ b/hog/manifold_forms.py @@ -33,7 +33,7 @@ from hog.quadrature import Tabulation from hog.symbolizer import Symbolizer from hog.logger import TimedLogger from hog.blending import GeometryMap, ExternalMap, IdentityMap -from hog.forms import Form +from hog.integrand import Form from hog.manifold_helpers import face_projection, embedded_normal diff --git a/hog/operator_generation/operators.py b/hog/operator_generation/operators.py index d8fd1167c70af2feb3ab02e54db1104b5eeaf20c..1052857a20166b5bb241f62f5508a886b51b5b27 100644 --- a/hog/operator_generation/operators.py +++ b/hog/operator_generation/operators.py @@ -70,7 +70,7 @@ from pystencils.typing.transformations import add_types from pystencils.backends.cbackend import CustomCodeNode from pystencils.typing.typed_sympy import FieldPointerSymbol -from hog.forms import Form +from hog.integrand import Form from hog.ast import Operations, count_operations from hog.blending import GeometryMap import hog.code_generation diff --git a/hog_tests/test_forms_new.py b/hog_tests/test_forms_new.py deleted file mode 100644 index fd96275e4d64583133c81b2277126bbe0e45692a..0000000000000000000000000000000000000000 --- a/hog_tests/test_forms_new.py +++ /dev/null @@ -1,81 +0,0 @@ -# HyTeG Operator Generator -# Copyright (C) 2024 HyTeG Team -# -# This program is free software: you can redistribute it and/or modify -# it under the terms of the GNU General Public License as published by -# the Free Software Foundation, either version 3 of the License, or -# (at your option) any later version. -# -# This program is distributed in the hope that it will be useful, -# but WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -# GNU General Public License for more details. -# -# You should have received a copy of the GNU General Public License -# along with this program. If not, see <https://www.gnu.org/licenses/>. - -from hog.forms_new import process_integrand -from hog.symbolizer import Symbolizer -from hog.element_geometry import TriangleElement -from hog.blending import AnnulusMap -from hog.function_space import LagrangianFunctionSpace -from hog.math_helpers import dot - -import hog.forms - - -def test_my_form() -> None: - - symbolizer = Symbolizer() - - trial = LagrangianFunctionSpace(2, symbolizer) - test = LagrangianFunctionSpace(2, symbolizer) - - geometry = TriangleElement() - blending = AnnulusMap() - - # This is what we have right now. - old_form_diffusion = hog.forms.diffusion( - trial, test, geometry, symbolizer, blending=blending - ) - old_form_mass = hog.forms.mass(trial, test, geometry, symbolizer, blending=blending) - - # Proposal: - # - # Have short (one-liner) functions that only describe the integrand. - # The arguments are automatically handed in from a populated IntegrandSymbols object. - # They may already be tabulated. - # Also, we can apply the chain rule already depending on the function spaces, such that writing forms becomes much - # easier. - def my_new_diffusion(*, grad_u, grad_v, dx, **kwargs): - return dot(grad_u, grad_v) * dx - - # We can add all arguments here that are given by the IntegrandSymbols object! - def my_new_mass(*, u, v, dx, **kwargs): - return u * v * dx - - integrands = [my_new_diffusion, my_new_mass] - - # Where the magic happens. - # The only variable argument is the form! - new_forms = [ - process_integrand( - integrand, - trial, - test, - geometry, - symbolizer, - blending, - is_symmetric=True, - ) - for integrand in integrands - ] - - # Prints "True" on my machine :) - print(new_forms[0].integrand == old_form_diffusion.integrand) - - # Mass is not tabulated correctly, but you get the point. - - -if __name__ == "__main__": - test_my_form()