diff --git a/CHANGELOG.md b/CHANGELOG.md index 1332d2f6..6ec77526 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -6,6 +6,11 @@ All notable changes to this project will be documented in this file. ### Added +- #188: Support multipatch interfaces between faces with different normal axes and with any edge or face orientation, several interfaces between the same pair of patches, and periodic and self-interfaces. +- #188: Add `Domain.patches`, `Domain.interface_map`, `Domain.exterior_sides`, `Domain.shared_vertices`, `Boundary.patch`, `Boundary.normal_axis` and `Interface.axis_map` to inspect the topology of multipatch domains, and the aliases `PatchVertex` and `SharedVertex` of `CornerBoundary` and `CornerInterface`. +- #188: Add the module `sympde.utilities.plotting`, with `plot_domain()` for plotting 2D and 3D domains with optional topology annotations, and `collect_topology_vertices()`, `format_topology()` and `print_topology()` for inspecting their shared vertices. `plot_domain()`, `collect_topology_vertices()` and `print_topology()` are available in `sympde.api`. +- #188: Add the module `sympde.topology.multipatch_gallery`, which builds a collection of multipatch domains (previously in PSYDAC), and the analytical mapping `TransposedPolarMapping`. +- #188: Document the multipatch topology in a new page of the user guide. - #194: Add `sympde.api` as the supported flat public API. It contains the objects used by PSYDAC, IGA-Python, GeLaTo and the SymPDE documentation. - #194: Add example scripts to the documentation, which are run by the documentation workflow. - #194: Draw expression trees in the documentation with Graphviz, generated from the current code at build time. @@ -18,6 +23,8 @@ All notable changes to this project will be documented in this file. ### Changed +- #143: Replace the `ornt` argument of `Interface` and `Boundary.join()` with a required `orientation` argument: `None` in 1D, `+1` or `-1` in 2D, and `(flag, sign1, sign2)` in 3D. Use `interface.orientation` instead of `interface.ornt`. +- #188: Accept the connectivity of `Domain.join()` only as `(minus, plus, orientation[, name])` tuples, and name repeated interfaces between the same pair of patches `A|B`, `A|B#2`, and so on. For mapped patches, the names are built from the logical patch names, and a mapped interface and its logical interface have the same name. - #200: Store `ldim` and `pdim` instead of `dim` in the topology of domain files. This is a breaking change: `Domain.from_file()` rejects files with `dim`. - #190: Convert the README to Markdown and expand the installation, testing, and documentation-build instructions. - #194: Use explicit internal imports and load optional dependencies only when needed. @@ -26,7 +33,7 @@ All notable changes to this project will be documented in this file. - #106: Require SymPy 1.12 or newer and support maintained SymPy releases below 1.15. - #106: Update SymPDE's symbolic objects for modern SymPy's structural arguments, traversal, substitution, and reconstruction requirements. - #106: Order the parameters of `Domain` like its `args` (`name`, `interiors`, `boundaries`, `mapping`, `logical_domain`, `connectivity`) and allow passing them positionally. `dim` is now keyword-only. -- #106: Allow passing `ornt`, `mapping` and `logical_domain` of an `Interface` positionally, in this order. +- #106: Allow passing `mapping` and `logical_domain` of an `Interface` positionally, in this order. - #106: Store only `axis` and `weight` in the `args` of the 1D forms (`Mass`, `Stiffness`, `Advection`, `AdvectionT`, `Bilaplacian`). Their name is now a class attribute. - [DEVELOPER] Write the HDF5 files of the unit tests to `tmp_path`, so that the tests can run in parallel without the `xdist_group` marker. - [DEVELOPER] Modernize the package metadata and consolidate the pytest and coverage configuration in `pyproject.toml`. @@ -38,8 +45,12 @@ All notable changes to this project will be documented in this file. ### Fixed -- #200: Keep the physical dimension of the mappings in `Domain.export()` and `Domain.from_file()`. +- #188: Reject interfaces in `Domain.join()` that join a face with itself or use a face that is already part of another interface. +- #188: Reject patches in `Domain.join()` of which only some have a mapping, instead of dropping the mapping of the domain. +- #140: Reject patches with different physical dimensions in `Domain.join()`, instead of creating an inconsistent domain. +- #144: Restore the interface names and exterior boundaries of a domain read with `Domain.from_file()` as they were exported, and read domains whose patches are `InteriorDomain` objects without a `dtype` instead of failing with an `AssertionError`. - #140: Fix the rendering of the `Domain.join()` API documentation. +- #200: Keep the physical dimension of the mappings in `Domain.export()` and `Domain.from_file()`. - #166: Restore warning-free documentation builds, including bibliography and mathematical notation support. - #106: Preserve symbol assumptions and free-symbol information with modern SymPy. - #106: Keep same-named functions from different function spaces distinct during symbolic simplification. @@ -54,6 +65,8 @@ All notable changes to this project will be documented in this file. ### Removed +- #73: Remove the `Edge` class and the `Interface.axis` property. Use `Boundary`, and `interface.minus.axis` and `interface.plus.axis`. +- #188: Remove `plot_domain()`, `plot_2d()`, `plot_3d()`, `plot_2d_single_patch()` and `plot_3d_single_patch()` from `sympde.utilities.utils`. Use `plot_domain()` from `sympde.api` instead. - #194: Remove the re-exports of public objects from the package and subpackage initializers. Replace imports such as `from sympde.topology import Square` with `from sympde.api import Square`. `from sympde import __version__` keeps working. - #194: Remove the `sympde.old_sympy_utilities` module. Use `sympy.utilities.iterables.is_sequence` and Python's `metaclass` keyword instead of its `is_sequence` and `with_metaclass`. - #194: Remove the outdated Jupyter notebooks and the Binder badge. The example scripts in the documentation replace them. diff --git a/doc/conf.py b/doc/conf.py index 370977ce..49c42977 100644 --- a/doc/conf.py +++ b/doc/conf.py @@ -51,6 +51,46 @@ 'sphinxcontrib.bibtex', ] +# Keep the custom commands used by ``math/feec.rst`` available in HTML +# output. The corresponding LaTeX definitions live in ``latex_macros.sty`` +# and are loaded separately for PDF builds below. +mathjax4_config = { + 'tex': { + 'macros': { + 'Alt': r'\mathrm{Alt}\,', + 'Curl': r'\nabla \times', + 'Div': r'\nabla \cdot', + 'Grad': r'\boldsymbol{\nabla}', + 'Hcurl': r'{H}(\mbox{curl}, \Omega)', + 'HcurlLogical': r'{H}(\mbox{curl}, \hat{\Omega})', + 'Hdiv': r'{H}(\mbox{div}, \Omega)', + 'HdivLogical': r'{H}(\mbox{div}, \hat{\Omega})', + 'Hgrad': r'H^1(\Omega)', + 'HgradLogical': r'H^1(\hat{\Omega})', + 'Ltwo': r'L^2(\Omega)', + 'LtwoLogical': r'L^2(\hat{\Omega})', + 'Vcurl': r'V_h(\mbox{curl}, \Omega)', + 'Vdiv': r'V_h(\mbox{div}, \Omega)', + 'Vgrad': r'V_h(\mbox{grad}, \Omega)', + 'Vltwo': r'V_h(L^2, \Omega)', + 'diff': r'\mathrm{d}\,', + 'divs': r'\mathrm{div}\,', + 'icurl': r'\imath^{1}', + 'id': r'\mathrm{id}\,', + 'idiv': r'\imath^{2}', + 'igrad': r'\imath^{0}', + 'iltwo': r'\imath^{3}', + 'nn': r'\mathbf{n}', + 'rots': r'\mathrm{rot}\,', + 'tr': r'\operatorname{tr}', + 'trace': r'\mathrm{Tr}\,', + 'uu': r'\mathbf{u}', + 'volume': r'\mathrm{vol}\,', + 'xx': r'\mathbf{x}', + }, + }, +} + # Add any paths that contain templates here, relative to this directory. templates_path = ['_templates'] diff --git a/doc/index.rst b/doc/index.rst index 6f19fd85..2a91e150 100644 --- a/doc/index.rst +++ b/doc/index.rst @@ -11,6 +11,7 @@ Welcome to sympde's documentation! :caption: Contents: quick + multipatch advanced examples/index diff --git a/doc/multipatch.rst b/doc/multipatch.rst new file mode 100644 index 00000000..65e99bc0 --- /dev/null +++ b/doc/multipatch.rst @@ -0,0 +1,669 @@ +Multipatch domains +****************** + +This chapter describes the symbolic multipatch topology in SymPDE. No mesh, +spline space, degrees of freedom, or MPI decomposition is created here. A +SymPDE multipatch domain answers topological questions such as: + +* which tensor-product patches constitute the domain; +* which boundary faces are identified; +* how the tangential coordinates of the two faces correspond; +* which faces remain on the exterior boundary; and +* which patch-local corners represent the same topological vertex. + +Psydac consumes this information later when it creates a discrete geometry. +Keeping the topology independent of a particular discretization is important: +the same domain may be discretized with different cell counts, degrees, knot +vectors, and process counts. + +The stable collection views described below are additive, and ordinary +single-patch construction with ``Line``, ``Square``, ``Cube``, mappings, +``domain.interior``, and ``domain.boundary`` is unchanged. The multipatch +interface API and serialized interface format intentionally replace older, +ambiguous representations; see `Migration from the previous multipatch API`_. + +The object model +================ + +A multipatch domain has four important levels:: + + logical NCube patches Square('A'), Cube('B'), ... + | + | optional patch mappings + v + physical patches F_A(A), F_B(B), ... + | + | Domain.join(...) + v + Domain + |-- patches stable tuple of patch interiors + |-- boundary faces not used by an interface + |-- interface_map name -> Interface + |-- exterior_sides stable tuple of exterior faces + `-- shared_vertices stable tuple of vertex classes + +These descriptive properties are views of the existing symbolic objects. +``domain.interior``, ``domain.subdomains``, ``domain.connectivity``, +``domain.boundary``, and ``domain.corners`` retain their original behavior. + +The logical dimension is denoted by ``dim`` below. A patch face has dimension +``dim - 1``. Logical axes are zero-based: a square has axes 0 and 1, and a +cube has axes 0, 1, and 2. + +Creating patches +================ + +The usual logical patches are ``Line``, ``Square``, and ``Cube``: + +.. code-block:: python + + from sympde.api import Line, Square, Cube + + interval = Line('I', bounds=(-1, 2)) + square = Square('A', bounds1=(0, 2), bounds2=(-1, 1)) + cube = Cube( + 'B', + bounds1=(0, 1), + bounds2=(0, 2), + bounds3=(-1, 1), + ) + +The parameters have the following meanings. + +``name`` + A unique patch name. Names are used for inspection, serialization, and + patch-specific discretization parameters. + +``bounds``, ``bounds1``, ``bounds2``, ``bounds3`` + The lower and upper logical coordinate of each tensor-product direction. + Every lower bound must be strictly smaller than its upper bound. Bounds + describe the parameter domain; they do not place a mapped patch in physical + space. + +``NCube(name, dim, min_coords, max_coords)`` is the dimension-independent +constructor. ``Line``, ``Square``, and ``Cube`` are its convenient 1D, 2D, +and 3D forms. Multipatch geometry is developed and tested primarily for these +three dimensions. + +Logical and physical patches +---------------------------- + +A mapping is applied to a logical patch by calling it: + +.. code-block:: python + + from sympde.api import AffineMapping, IdentityMapping, Square + + logical_a = Square('A') + logical_b = Square('B') + + patch_a = IdentityMapping('F_A', dim=2)(logical_a) + patch_b = AffineMapping( + 'F_B', dim=2, + c1=1, c2=0, + a11=1, a12=0, + a21=0, a22=1, + )(logical_b) + +For ``Mapping(name, dim=dim)``, ``dim`` sets both the logical and physical +dimension. More generally, ``ldim`` and ``pdim`` may be supplied by mapping +classes which support embedded geometries. Parameters such as ``c1`` and +``a11`` are specific to the mapping class. An ``AffineMapping`` in 2D is + +.. math:: + + F(\eta_1,\eta_2) = + \begin{pmatrix}c_1\\c_2\end{pmatrix} + + \begin{pmatrix}a_{11}&a_{12}\\a_{21}&a_{22}\end{pmatrix} + \begin{pmatrix}\eta_1\\\eta_2\end{pmatrix}. + +If every patch passed to ``Domain.join`` is mapped, the joined domain receives +a ``MultiPatchMapping`` and a corresponding ``logical_domain``. Mixing mapped +and unmapped patches is representable, but no aggregate multipatch mapping is +created. Downstream geometry and variational code generally expects the +patches to be consistently all mapped or all logical. + +A mapping may also be applied to an already joined logical domain. In that +case its interfaces are mapped together with its interiors and boundaries. +Each mapped interface retains the original interface as its ``logical_domain`` +and receives an ``InterfaceMapping`` for interface integrals. + +Domain gallery +============== + +Reusable symbolic examples live in ``sympde.topology.multipatch_gallery``. +They have no dependency on Psydac or on a discretization: + +.. code-block:: python + + from sympde.topology.multipatch_gallery import ( + build_cartesian_multipatch_domain_2d, + build_multipatch_domain_2d, + ) + + annulus = build_multipatch_domain_2d('annulus_4') + square_with_a_hole = build_multipatch_domain_2d('square_8') + +``build_multipatch_domain_2d`` also provides ``square_2``, ``square_4``, +``square_6``, ``square_9``, ``annulus_3``, ``curved_L_shape``, ``pretzel``, +``pretzel_f``, ``pretzel_annulus``, and ``pretzel_debug``. +Each named geometry has a dedicated public builder in the same module; the +name-based function is only a convenience dispatcher. + +``build_cartesian_multipatch_domain_2d`` creates a rectangular grid from a 2D +layout. Every non-``None`` entry creates one patch and a ``None`` entry leaves +a hole. The values themselves are not interpreted by SymPDE, so the same array +may later carry patch-specific cell counts for a discretizer. Both identity +and polar patch mappings are available. + +The gallery also contains 3D domains. ``build_two_cube_3d`` constructs two +sheared cubes with a cross-axis interface, and ``build_torus_3d`` constructs +the four-patch hollow or solid toroidal volume. In the solid torus, the faces +at minor radius zero collapse onto the centreline of the tube. They are still +exterior sides of the domain, so they must be excluded when boundary +conditions are imposed on ``domain.boundary``. Additional diagnostic domains +exercise oriented cube faces and three patches sharing an edge. +``build_multipatch_domain_3d`` selects these builders by name, including +``two_cube`` and ``torus``, and forwards optional arguments to the selected +builder. + +The gallery module can also plot any registered 2D or 3D domain directly. +For example, the following commands open an ordinary annulus plot and save an +annotated two-cube topology plot, respectively: + +.. code-block:: console + + python -m sympde.topology.multipatch_gallery annulus_4 + python -m sympde.topology.multipatch_gallery two_cube \ + --topology --no-show --output two_cube.png + +The corresponding Python entry point is ``plot_multipatch_domain``. Its +``builder_options`` argument forwards domain-specific parameters to the +selected 2D or 3D builder. + +Joining patches +=============== + +The public constructor is + +.. code-block:: python + + Domain.join(patches, interfaces=interfaces, name=name) + +``patches`` + A non-empty ordered list or tuple of ``Domain`` objects. All patches must + have the same logical dimension and the same physical dimension, which is + the logical dimension for a patch without a mapping. Either all patches or + none have a mapping. The order defines the integer patch indices accepted + in the interface specifications. Patch objects, physical interior names, + and logical names used by serialization must all be unique. + +``interfaces`` + A list or tuple of interface descriptions. Each description is the tuple + ``(minus, plus, orientation)`` or + ``(minus, plus, orientation, name)``. The optional explicit name is + preserved exactly and must be unique. The whole collection is a sequence + because several distinct interfaces may connect the same pair of patches. + ``connectivity`` is an equivalent positional/keyword spelling. + Omitting both arguments is equivalent to providing an empty tuple. + +``name`` + The name of the resulting domain. + +If one patch and an empty connectivity are supplied, that patch is returned +unchanged. A one-patch domain with non-empty connectivity follows the normal +join path and can therefore contain self-interfaces. + +Describing an interface side +---------------------------- + +Each side is the triple ``(patch, axis, ext)``: + +``patch`` + Either the patch object itself or its integer position in ``patches``. + +``axis`` + The logical direction normal to the face, in ``range(dim)``. + +``ext`` + ``-1`` selects the minimum-coordinate face and ``+1`` the + maximum-coordinate face. + +For example: + +.. code-block:: python + + left = (A, 0, -1) + top = (A, 1, +1) + +An interface combines two side triples with their orientation: + +.. code-block:: python + + interface = ((A, 0, +1), (B, 1, -1), -1) + +``Domain.join`` resolves each side to an existing symbolic ``Boundary`` +object. + +``minus`` and ``plus`` are ordered labels. They do not mean that one patch is +geometrically smaller, coarser, or located on a particular side in physical +space. The ordering establishes the direction in which the orientation is +read and remains the ordering stored by the resulting ``Interface``. + +Interface orientation +===================== + +The public orientation value depends on the patch dimension: + +``1D`` + ``None``. An endpoint has no tangential coordinate. + +``2D`` + One sign, ``+1`` or ``-1``. It says whether the coordinates along the two + edges increase in the same or opposite directions. + +``3D`` + ``(flag, sign1, sign2)``, where all three entries are ``+1`` or ``-1``. + ``flag=+1`` preserves the order of the two tangential coordinates; + ``flag=-1`` exchanges their order. ``sign1`` and ``sign2`` then give the + directions of the first and second resulting correspondences. + +Thus the common orientations are: + +.. code-block:: python + + identity_2d = +1 + reversed_2d = -1 + + identity_3d = (+1, +1, +1) + swap_3d = (-1, +1, +1) + swap_and_flip_3d = (-1, +1, -1) + +The ordered tangential axes are obtained by removing the face-normal axis from +``range(dim)`` without otherwise changing the order. The most convenient +way to inspect the decoded result is ``interface.axis_map``. It returns +``(minus_axis, plus_axis, direction)`` for every tangential direction. + +For example, consider a 3D interface whose minus normal is axis 0 and whose +plus normal is axis 1: + +.. code-block:: python + + from sympde.api import Cube, Domain + + patch_a = Cube('A') + patch_b = Cube('B') + domain = Domain.join( + patches=[patch_a, patch_b], + interfaces=[( + (0, 0, +1), + (1, 1, -1), + (-1, +1, -1), + )], + name='Omega', + ) + interface, = domain.interface_map.values() + assert interface.orientation == (-1, +1, -1) + assert interface.axis_map == ((1, 2, 1), (2, 0, -1)) + +The first tuple says minus axis 1 maps to plus axis 2 in the same direction. +The second says minus axis 2 maps to plus axis 0 in the opposite direction. +Normal axes are intentionally absent from the orientation value: they already +belong to the two side triples. The orientation is decoded against each +face's own normal axis. Therefore the two normal axes need not be equal. + +An ``Interface`` deliberately has no single ``axis`` property because its two +faces may have different normal axes. Read ``interface.minus.axis`` and +``interface.plus.axis`` explicitly. New multipatch code may use the clearer +aliases ``interface.minus_side.normal_axis``, +``interface.plus_side.normal_axis``, and ``interface.axis_map``. + +Identity orientation means that tangential positions retain their order and +all signs are positive. It does not mean that the two normal axes must be +equal, their extremities must differ, or the physical mappings must be +identical. + +A complete 2D example +===================== + +This example joins two mapped squares twice. The touching seam preserves its +tangential coordinate, while the outer seam reverses it: + +.. code-block:: python + + from sympde.api import ( + AffineMapping, Domain, IdentityMapping, Square, + ) + + A = IdentityMapping('F_A', dim=2)(Square('A')) + B = AffineMapping( + 'F_B', dim=2, + c1=1, c2=0, + a11=1, a12=0, + a21=0, a22=1, + )(Square('B')) + + omega = Domain.join( + patches=[A, B], + interfaces=[ + ((A, 0, +1), (B, 0, -1), +1), + ((A, 0, -1), (B, 0, +1), -1), + ], + name='twisted_strip', + ) + +This example illustrates two independent facts: + +* the same ordered patch pair can have more than one interface; and +* the orientation belongs to each interface, not to the patch pair. + +Interface names are generated from the side-domain names; for mapped patches, +from the names of their logical patches, so that ``F(A)`` and ``G(B)`` give +``A|B``. If that base name is repeated, the interfaces receive ``#2``, +``#3``, and so on, in the order of their faces (patch, axis and extremity of +the minus side, then of the plus side), so that the names do not depend on +the order of the connectivity. The names are checked against all previously generated +names, including names which already contain ``#`` or ``|``, so one interface +cannot overwrite another. A mapped interface and its logical interface have +the same name. The names are stable keys in ``omega.interface_map``, which +is ordered by name. ``omega.connectivity`` is +the unchanged symbolic storage object behind that descriptive view. + +Self-interfaces +=============== + +Two distinct faces of one patch may be identified: + +.. code-block:: python + + P = Square('P') + periodic_x = Domain.join( + patches=[P], + interfaces=[((P, 0, -1), (P, 0, +1), +1)], + name='periodic_x', + ) + +Both ``interface.minus.domain`` and ``interface.plus.domain`` are ``P``. The +two boundary objects remain distinct because their extremities differ. A +self-interface is part of the topology; whether a particular discretization +consumer realizes it as periodic basis functions or as an explicit constraint +is a separate question. + +Shared vertices and corners +=========================== + +Joining faces also identifies their corners. SymPDE computes this transitively: +if an interface identifies a corner of A with a corner of B and another +identifies that B corner with a corner of C, all three belong to one +shared vertex. ``SharedVertex`` is the descriptive alias of the historical +``CornerInterface`` class, while ``PatchVertex`` aliases the patch-local +``CornerBoundary`` class. The historical names remain canonical for SymPy +printing and serialization. + +For each interface, SymPDE enumerates the ``2**(dim - 1)`` corners of its minus +face, maps their tangential extremities with +``interface.axis_map``, and unions the resulting minus and plus +``PatchVertex`` objects. This supports, for +example, three or more 2D patches meeting at one vertex and 3D faces with +permuted axes. + +Use ``domain.interface_map`` and ``domain.shared_vertices`` for stable +multipatch collection views: + +.. code-block:: python + + for name, interface in omega.interface_map.items(): + print(name) + print(interface.minus_side.patch, + interface.minus_side.normal_axis, + interface.minus_side.ext) + print(interface.plus_side.patch, + interface.plus_side.normal_axis, + interface.plus_side.ext) + print(interface.orientation) + print(interface.axis_map) + + for vertex in omega.shared_vertices: + for local_vertex in vertex.corners: + print(local_vertex.domain, local_vertex.coordinates) + +``PatchVertex.coordinates`` is a tuple of zeros and ones, one per logical +axis. It identifies the minimum or maximum vertex in patch-local coordinates. +``domain.shared_vertices`` contains shared vertices induced by connectivity; +it is not a list of every unshared exterior vertex. The symbolic +``domain.corners`` property remains available with its original +``None``/object/``Union`` cardinality behavior. + +Inspecting and plotting 2D and 3D topology +========================================== + +The plotting utilities can turn the symbolic connectivity into structured +vertex data, a textual summary, or an annotated figure. No mesh or Psydac +discretization is created: + +.. code-block:: python + + from sympde.api import ( + collect_topology_vertices, + plot_domain, + print_topology, + ) + + vertices = collect_topology_vertices(omega) + for vertex in vertices: + print(vertex.index, vertex.patches, vertex.is_boundary) + for incidence in vertex.incidences: + print( + incidence.patch_index, + incidence.patch_name, + incidence.logical_corner, + ) + + print_topology(omega) + + figure = plot_domain( + omega, + draw=False, + refinement=60, + isolines=True, + topology=True, + vertex_labels=True, + ) + figure.savefig('omega-topology.png', dpi=180, bbox_inches='tight') + +``collect_topology_vertices`` differs deliberately from ``domain.corners``: +it begins with every corner of every patch and merges the shared-corner +equivalence classes. Its result therefore includes unshared exterior +vertices as well as interior or multiply represented vertices. Patch names +and integer patch indices are both recorded. + +In a 2D annotated topology plot, each side of an interface is drawn slightly +inside its own patch: the minus side solid, the plus side dashed, in a colour +shared by both sides. The arrow of the minus side points along its increasing +tangential coordinate, and the arrow of the plus side shows where that +coordinate goes on the plus side: it points against the plus side's own +coordinate for a reversing interface, such as the twisted strip. If identified sides are separated in the drawing, they are labelled +``I0-`` and ``I0+`` and a dotted curve connects their midpoints; this makes +cut-open periodic or quotient topologies visible without moving the patches. +The exterior boundary is drawn in black, and the vertices are marked white on +the boundary and black in the interior. Patches get their own colour from a +colour-blind safe palette, or, when there are more patches than colours, +neighbouring patches get different colours. + +The 3D view is intentionally lighter. It draws the exterior patch faces as +translucent surfaces, fills the minus face of every interface in the +interface colour, and places arrows along the two tangential coordinates of the +minus face, solid, and along their images on the plus face, dashed. All arrows +are labelled with the coordinates of the minus face, so the arrows of the two +faces coincide when the declared orientation matches the geometry. The result is +intended for connectivity inspection rather than publication-quality +rendering of a complex volume decomposition. :func:`print_topology` reports +the same topology as text, including the signed ``axis_map`` of every +interface. + +Complete-vertex inspection and annotated plotting support 2D and 3D domains +whose physical dimension equals their logical dimension. The ordinary +geometric ``plot_domain`` paths remain available when ``topology=False``. + +Selecting patch subdomains +========================== + +``domain.get_subdomain(names)`` returns the topology induced by the requested +patch names. An interface is retained when both of its sides belong to the +selection. This includes repeated interfaces and interfaces from a patch to +itself. When exactly one side is selected, that side becomes part of the new +exterior boundary. + +For mapped patches, the returned domain receives the corresponding patch +mapping or ``MultiPatchMapping`` and a logical subdomain with matching +connectivity. Interface objects and their names are retained, so orientation, +``axis_map``, ``InterfaceMapping``, and logical-interface identity are not +reconstructed or lost. + +``Domain.join`` accepts atomic patch domains without existing connectivity. +Use ``get_subdomain`` to select from an existing joined domain; nested joined +domains are not implicitly flattened by another call to ``Domain.join``. + +Exterior boundary construction +============================== + +``Domain.join`` starts from every boundary face of every patch and removes all +faces used as an interface side. The remaining faces form ``domain.boundary``. +Consequently a weak form integrated on ``domain.boundary`` sees only the +exterior boundary, not internal seams. + +The join operation does not infer connectivity by comparing physical +coordinates. Every interface must be declared. Conversely, a declared face +is removed from the exterior boundary even if its two physical mappings do not +actually coincide. + +What SymPDE validates +===================== + +``Domain.join`` validates the structural contract: + +* ``patches`` and ``connectivity`` have the expected container types; +* every input to ``Domain.join`` is an atomic, unconnected patch domain; +* patch objects, interior names, and serialized logical names are unique; +* all patch objects have a common logical dimension and a common physical + dimension; +* either all patches or none have a mapping; +* patch indices are in range and patch objects belong to ``patches``; +* axes are valid and extremities are ``-1`` or ``+1``; +* every face is part of at most one interface, and no face is joined with + itself; +* every interface has an explicit dimension-appropriate orientation; and +* every orientation flag and sign belongs to ``{-1, +1}``. + +Axes, extremities, and orientation flags must be integer values. Floating +point values are rejected even when they are numerically equal to an integer, +and booleans are not interpreted as ``0`` or ``1``. + +It does *not* establish geometric or discretization compatibility. The caller +is responsible for ensuring that: + +* the mapped physical faces coincide to the desired tolerance; +* patch maps are regular and have the intended physical dimension; +* the face parameterizations satisfy the declared orientation; and +* later discrete trace meshes and spaces meet the requirements of the chosen + Psydac operator. + +This separation permits general quotient topologies, repeated patch-pair +interfaces, self-identifications, cross-axis gluing, reversed coordinates, and +3D face-axis permutations without baking a mesh policy into the symbolic +geometry. + +Migration from the previous multipatch API +========================================== + +The explicit orientation contract removes compatibility paths whose meaning +was ambiguous for cross-axis or three-dimensional interfaces. Existing +multipatch callers should apply the following replacements: + +.. list-table:: + :header-rows: 1 + + * - Previous API or data + - Current replacement + * - ``Boundary.join(other, ornt=...)`` + - ``Boundary.join(other, orientation)`` + * - ``Interface(..., ornt=...)`` + - ``Interface(..., orientation)`` + * - ``interface.ornt`` + - ``interface.orientation`` + * - ``interface.axis`` + - ``interface.minus.axis`` and ``interface.plus.axis`` + * - ``Edge`` + - ``Boundary`` + * - interface dictionaries + - ``(minus, plus, orientation[, name])`` tuples + * - omitted interface orientation + - an explicit dimension-appropriate orientation + * - ``sympde.utilities.utils.plot_domain`` + - ``sympde.api.plot_domain`` + +The topology file-format migration is described in `Serialization`_. In +particular, legacy interface entries without an explicit orientation are not +loaded implicitly. + +Serialization +============= + +``domain.todict()`` contains the interior, exterior boundary, and connectivity. +``domain.export('topology.h5')`` writes this dictionary as ``topology.yml`` in +an HDF5 file. The exterior boundary is always represented as a list: an empty +list for a fully closed domain, a one-element list for one exterior side, and a +longer list otherwise. The reader also accepts the historical single-side +dictionary representation. + +Every serialized interface contains three entries: its minus side, plus side, +and explicit orientation. Its connectivity dictionary key is the interface +name and is preserved when the file is loaded. For example: + +.. code-block:: yaml + + connectivity: + A|B: + - {patch: A, axis: '0', ext: '1', mapping: None, name: A_0} + - {patch: B, axis: '1', ext: '-1', mapping: None, name: B_1} + - [-1, 1, -1] + +Load topology with ``Domain.from_file(filename)``. Files written by the old +orientation-less format are deliberately rejected: assigning an implicit +orientation would be ambiguous and could silently change the topology. Such +files must be regenerated or migrated. Files containing the experimental +``permutation``/``directions`` orientation mapping are likewise obsolete; +current files store the same compact values accepted by ``Domain.join``. +The loader derives the exterior boundary from the connectivity and rejects a +file when that result disagrees with boundary metadata that is present. Missing +or null legacy metadata is treated as unavailable and replaced by the derived +boundary. + +A domain whose patches are ``InteriorDomain`` objects without a ``dtype`` has +no parametric box from which to derive its boundary. The loader restores its +exterior boundaries and interfaces as they were stored instead, including +their names and any boundary without an axis or extremity. Such patches are +only supported without a mapping. + +SymPDE's ``Domain.export`` stores symbolic topology only. A Psydac geometry +file additionally stores spline or NURBS knots, control points, weights, and +parallelizable mapping data; see Psydac's multipatch discretization tutorial. + +Current scope +============= + +The symbolic representation covers the following independently: + +* 1D endpoint, 2D edge, and 3D face connections; +* unequal normal-axis indices on the two patches; +* all two 2D and all eight 3D face orientations; +* several interfaces between the same two patches; +* interfaces from a patch to itself; and +* arbitrary transitive corner incidence generated by those faces. + +That is the scope of the *topology model*. A downstream consumer may support +only a subset. Psydac's standard variational and parallel interface paths +support signed cross-axis orientations for matching traces, while conforming +projectors and nested unequal trace couplings retain additional restrictions. +Those restrictions belong to Psydac and should not be encoded as assumptions +in a SymPDE domain. diff --git a/sympde/api/__init__.py b/sympde/api/__init__.py index 4cebf564..e6f98852 100644 --- a/sympde/api/__init__.py +++ b/sympde/api/__init__.py @@ -153,7 +153,11 @@ from sympde.exterior.form import DifferentialForm from sympde.printing.latex import latex -from sympde.utilities.utils import plot_domain +from sympde.utilities.plotting import ( + collect_topology_vertices, + plot_domain, + print_topology, +) __all__ = ( 'Constant', @@ -284,6 +288,8 @@ 'ld', 'wedge', # Output and utilities + 'collect_topology_vertices', 'latex', 'plot_domain', + 'print_topology', ) diff --git a/sympde/expr/evaluation.py b/sympde/expr/evaluation.py index f5b50f34..919c7b8a 100644 --- a/sympde/expr/evaluation.py +++ b/sympde/expr/evaluation.py @@ -482,7 +482,7 @@ class TerminalExpr(CalculusFunction): Parameters ---------- - expr : Expr or Matrix or ImmutableDenseMatrix or LogicalExpr + expr : Expr | Matrix | ImmutableDenseMatrix | LogicalExpr The mathematical expression. domain : BasicDomain (from sympde.topology.basic) diff --git a/sympde/expr/tests/test_expr.py b/sympde/expr/tests/test_expr.py index 8d0114cd..dbd01b92 100644 --- a/sympde/expr/tests/test_expr.py +++ b/sympde/expr/tests/test_expr.py @@ -1332,44 +1332,13 @@ def test_linearity_bilinear_form_2d_1(): #============================================================================== def test_interface_2d_1(): -# # ... -# def two_patches(): -# -# from sympde.topology import Connectivity, Interface -# -# A = Square('A') -# B = Square('B') -# -# A = A.interior -# B = B.interior -# -# bnd_A_1 = Boundary(r'\Gamma_1', A, axis=0, ext=-1) -# bnd_A_2 = Boundary(r'\Gamma_2', A, axis=0, ext=1) -# bnd_A_3 = Boundary(r'\Gamma_3', A, axis=1, ext=-1) -# bnd_A_4 = Boundary(r'\Gamma_4', A, axis=1, ext=1) -# -# bnd_B_1 = Boundary(r'\Gamma_1', B, axis=0, ext=-1) -# bnd_B_2 = Boundary(r'\Gamma_2', B, axis=0, ext=1) -# bnd_B_3 = Boundary(r'\Gamma_3', B, axis=1, ext=-1) -# bnd_B_4 = Boundary(r'\Gamma_4', B, axis=1, ext=1) -# -# connectivity = Connectivity() -# connectivity['I'] = Interface('I', bnd_A_2, bnd_B_1, ornt=1) -# Omega = Domain('Omega', -# interiors=[A, B], -# boundaries=[bnd_A_1, bnd_A_2, bnd_A_3, bnd_A_4, bnd_B_1, bnd_B_2, bnd_B_3, bnd_B_4], -# connectivity=connectivity) -# -# return Omega -# # ... - def two_patches(): A = Square('A') B = Square('B') patches = [A, B] - connectivity = [((A, 0, 1), (B, 0, -1), 1)] + connectivity = [((A, 0, 1), (B, 0, -1), +1)] return Domain.join(patches, connectivity, 'Omega') # create a domain with an interface @@ -1381,7 +1350,7 @@ def two_patches(): assert isinstance(interface, Interface) assert interface.minus is domain.subdomains[0].get_boundary(axis=0, ext= 1) assert interface.plus is domain.subdomains[1].get_boundary(axis=0, ext=-1) - assert interface.ornt == 1 + assert interface.orientation == +1 assert interface.dim == 2 assert interface.mapping is None assert interface.logical_domain is None @@ -1405,7 +1374,7 @@ def test_interface_integral_1(): B = Square('B') domains = [A, B] - connectivity = [((0, 0, 1), (1, 0, -1), 1)] + connectivity = [((0, 0, 1), (1, 0, -1), +1)] domain = Domain.join(domains, connectivity, 'domain') x,y = domain.coordinates @@ -1466,7 +1435,7 @@ def test_interface_integral_2(): patches = [A, B] - connectivity = [((0, 0, 1), (1, 0, -1), 1)] + connectivity = [((0, 0, 1), (1, 0, -1), +1)] domain = Domain.join(patches, connectivity, 'domain') x,y = domain.coordinates @@ -1504,8 +1473,8 @@ def test_interface_integral_3(): patches = [A, B, C] - connectivity = [((0, 0, 1), (1, 0, -1), 1), - ((1, 0, 1), (2, 0, -1), 1)] + connectivity = [((0, 0, 1), (1, 0, -1), +1), + ((1, 0, 1), (2, 0, -1), +1)] domain = Domain.join(patches, connectivity, 'domain') x,y = domain.coordinates @@ -1543,7 +1512,7 @@ def test_interface_integral_4(): B = Square('B') patches = [A, B] - connectivity = [((0, 0, 1), (1, 0, -1), 1)] + connectivity = [((0, 0, 1), (1, 0, -1), +1)] domain = Domain.join(patches, connectivity, 'AB') x,y = domain.coordinates diff --git a/sympde/topology/analytical_mapping.py b/sympde/topology/analytical_mapping.py index 05b9760b..f976057a 100644 --- a/sympde/topology/analytical_mapping.py +++ b/sympde/topology/analytical_mapping.py @@ -37,6 +37,19 @@ class PolarMapping(Mapping): _ldim = 2 _pdim = 2 +#============================================================================== +class TransposedPolarMapping(Mapping): + """Polar 2D mapping with angle and radius coordinates transposed. + + The first logical coordinate is the angle and the second one is the + radius. This is the coordinate ordering opposite to :class:`PolarMapping`. + """ + _expressions = {'x': 'c1 + (rmin*(1-x2)+rmax*x2)*cos(x1)', + 'y': 'c2 + (rmin*(1-x2)+rmax*x2)*sin(x1)'} + + _ldim = 2 + _pdim = 2 + #============================================================================== class TargetMapping(Mapping): """ diff --git a/sympde/topology/basic.py b/sympde/topology/basic.py index 522e0862..7a1ceabc 100644 --- a/sympde/topology/basic.py +++ b/sympde/topology/basic.py @@ -3,8 +3,9 @@ from collections import abc +from numbers import Integral -from sympy.core import Basic, Symbol +from sympy.core import Basic, Symbol, Integer from sympy.core.containers import Tuple from sympde.core.basic import ( @@ -15,6 +16,13 @@ restore_string_argument, ) + +def _as_integer(value, name): + """Return an integer input without silently truncating other scalars.""" + if isinstance(value, bool) or not isinstance(value, (Integral, Integer)): + raise TypeError(f'{name} must be an integer') + return int(value) + #============================================================================== class BasicDomain(Basic): _dim = None @@ -51,6 +59,7 @@ def _sympystr(self, printer): sstr = printer.doprint return '{}'.format(sstr(self.name)) + #============================================================================== class InteriorDomain(BasicDomain): """ @@ -290,10 +299,15 @@ def __new__(cls, name, domain, axis=None, ext=None, mapping=None, logical_domain logical_domain = restore_optional_argument(logical_domain) if axis is not None: - assert isinstance(axis, int) + axis = _as_integer(axis, 'boundary axis') + if not 0 <= axis < domain.dim: + raise ValueError( + f'boundary axis must be between 0 and {domain.dim - 1}') if ext is not None: - assert isinstance(ext, int) + ext = _as_integer(ext, 'boundary extremity') + if ext not in (-1, 1): + raise ValueError('boundary extremity must be either -1 or 1') obj = new_basic(cls, name, domain, axis, ext, mapping, logical_domain) @@ -314,10 +328,33 @@ def domain(self): def axis(self): return restore_int_argument(self.args[2]) + @property + def normal_axis(self): + """Logical axis normal to this patch side. + + This is the descriptive multipatch alias of :attr:`axis`. + """ + return self.axis + @property def ext(self): return restore_int_argument(self.args[3]) + @property + def patch(self): + """Patch interior owning this side; an alias of :attr:`domain`.""" + return self.domain + + @property + def patch_dim(self): + """Logical dimension of the patch owning this side.""" + return self.dim + + @property + def intrinsic_dim(self): + """Intrinsic dimension of this codimension-one side.""" + return self.dim - 1 + @property def mapping(self): return self._mapping @@ -349,22 +386,27 @@ def rotate(self, *directions): raise NotImplementedError('only 2d case is available') - def join(self, boundary, ornt=None): - from sympde.topology.mapping import InterfaceMapping - # TODO be careful with '|' in psydac + def join(self, boundary, orientation, *, name=None, logical_name=None): + if not isinstance(boundary, Boundary): + raise TypeError('boundary must be a Boundary') + if self.mapping and boundary.mapping: + # Imported lazily to keep the basic topology layer independent of + # the mapping/domain import cycle. + from sympde.topology.mapping import InterfaceMapping int_map = InterfaceMapping(self.mapping , boundary.mapping) a,b = self.logical_domain, boundary.logical_domain - l_name = '{l}|{r}'.format(l=a.domain.name, r=b.domain.name) - int_logical_domain = Interface(l_name, a,b, ornt=ornt) + l_name = logical_name or '{l}|{r}'.format(l=a.domain.name, r=b.domain.name) + int_logical_domain = Interface(l_name, a, b, orientation=orientation) else: int_map = None int_logical_domain = None - name = '{l}|{r}'.format(l=self.domain.name, r=boundary.domain.name) + name = name or '{l}|{r}'.format(l=self.domain.name, r=boundary.domain.name) interface = Interface(name, self, boundary, mapping=int_map, - logical_domain=int_logical_domain, ornt=ornt) + logical_domain=int_logical_domain, + orientation=orientation) return interface def _sympystr(self, printer): @@ -389,7 +431,7 @@ def todict(self): #============================================================================== class CornerBoundary(BasicDomain): """ - Represents an undefined corner over a domain in 2D. + Represents a vertex as the intersection of a patch's boundary faces. """ def __new__(cls, *boundaries): @@ -432,7 +474,7 @@ def _sympystr(self, printer): #============================================================================== class CornerInterface(BasicDomain): """ - Represents a shared corner over multiple patches in 2D. + Represents a vertex shared by multiple patches. """ def __new__(cls, *corners): @@ -460,6 +502,34 @@ def _sympystr(self, printer): corners = ', '.join(sstr(b) for b in self.corners) return 'CornerInterface({})'.format(corners) + +# Descriptive multipatch aliases. The original class names remain canonical +# SymPy and serialization names for backwards compatibility. +PatchVertex = CornerBoundary +SharedVertex = CornerInterface + + +#============================================================================== +def _interface_axis_map(dim, minus_axis, plus_axis, orientation): + """Decode a validated compact orientation into a signed axis map.""" + dim = int(dim) + minus_tangents = tuple(axis for axis in range(dim) if axis != minus_axis) + plus_tangents = tuple(axis for axis in range(dim) if axis != plus_axis) + + if dim == 1: + return () + if dim == 2: + return ((minus_tangents[0], plus_tangents[0], orientation),) + + flag, sign1, sign2 = orientation + plus_positions = (0, 1) if flag == 1 else (1, 0) + return tuple( + (minus_axis, plus_tangents[plus_position], direction) + for minus_axis, plus_position, direction + in zip(minus_tangents, plus_positions, (sign1, sign2)) + ) + + #============================================================================== class Interface(BasicDomain): """ @@ -476,6 +546,19 @@ class Interface(BasicDomain): bnd_plus : Boundary Boundary on the "plus" side of the interface. + orientation : None, int, or tuple[int, int, int] + Compact orientation from the minus face to the plus face. Let ``M`` + and ``P`` be the ordered tangential axes obtained by removing the + normal axis of the minus and plus face, respectively, from + ``range(dim)``. + + In 1D the orientation is ``None``. In 2D it is one sign, and maps + ``M[0]`` to ``P[0]`` with that direction. In 3D it is + ``(flag, sign1, sign2)``. If ``flag`` is ``+1``, ``M[0]`` and ``M[1]`` + map to ``P[0]`` and ``P[1]``. If ``flag`` is ``-1``, they map to + ``P[1]`` and ``P[0]``. ``sign1`` and ``sign2`` give the respective + directions. Every flag and sign is either ``-1`` or ``+1``. + mapping : Mapping, optional Mapping from the logical domain to the physical domain, if available. @@ -483,24 +566,48 @@ class Interface(BasicDomain): Logical domain associated with the interface, if available. It should be consistent with the mapping if provided. - ornt : int | Iterable[int], optional - Orientation of the interface. For 1D interfaces, this is not needed and - should be set to None. For 2D interfaces, this should be either -1 or 1. - For 3D interfaces, this should be a tuple of three integers, each being - either -1 or 1. - Notes ----- - The orientations are specified in the same manner as in GeoPDES, see e.g. - - and - T. Dokken, E. Quak, V. Skytt. Requirements from Isogeometric Analysis for changes in product design ontologies, 2010. + ``minus`` and ``plus`` are ordered labels: they define the direction in + which ``orientation`` is read. They do not describe geometric position or + mesh size. The two faces may have different normal axes; those axes are + part of ``bnd_minus`` and ``bnd_plus`` and are deliberately not repeated + in the orientation. + + This compact orientation convention follows the multipatch convention + used by GeoPDEs; see its `multipatch geometry specification + `_ + and T. Dokken, E. Quak, V. Skytt, *Requirements from Isogeometric + Analysis for Changes in Product Design Ontologies* (2010). + + :attr:`axis_map` exposes the decoded convention as + ``(minus_axis, plus_axis, direction)`` triples. + + Examples + -------- + A cube interface may exchange its two tangential axes and reverse one of + them, even when the faces have different normal axes: + + >>> from sympde.api import Cube, Interface + >>> A = Cube('A') + >>> B = Cube('B') + >>> interface = Interface( + ... 'A|B', + ... A.get_boundary(axis=0, ext=+1), + ... B.get_boundary(axis=1, ext=-1), + ... orientation=(-1, +1, -1), + ... ) + >>> interface.orientation + (-1, 1, -1) + >>> interface.axis_map + ((1, 2, 1), (2, 0, -1)) + """ - def __new__(cls, name, bnd_minus, bnd_plus, ornt=None, - mapping=None, logical_domain=None): + def __new__(cls, name, bnd_minus, bnd_plus, orientation, mapping=None, + logical_domain=None): name = restore_string_argument(name) - ornt = restore_container_argument(ornt) + orientation = restore_container_argument(orientation) mapping = restore_optional_argument(mapping) logical_domain = restore_optional_argument(logical_domain) @@ -534,23 +641,45 @@ def __new__(cls, name, bnd_minus, bnd_plus, ornt=None, if logical_domain.dim != ldim: raise ValueError(f'Logical domain should have dimension = {ldim}, got {logical_domain.dim} instead') - # Check that orientation is provided in the correct format depending on the dimension + # Validate the public input and convert it to the canonical SymPy + # representation stored in ``args``. ``axis_map`` can consequently + # decode ``self.orientation`` without validating it a second time. if ldim == 1: - assert ornt is None, 'ornt is not needed for 1D interfaces' + if orientation is not None: + raise ValueError('a 1D interface orientation must be None') + orientation_arg = None elif ldim == 2: - assert ornt in (-1, 1), 'ornt must be either -1 or 1 for 2D interfaces' + try: + orientation = _as_integer( + orientation, 'a 2D interface orientation') + except TypeError as error: + raise TypeError( + 'a 2D interface orientation must be an integer') \ + from error + if orientation not in (-1, 1): + raise ValueError('a 2D interface orientation must be +1 or -1') + orientation_arg = Integer(orientation) elif ldim == 3: - ornt = tuple(ornt) - assert len(ornt) == 3, 'ornt must be a tuple of length 3 for 3D interfaces' - assert all(o in (-1, 1) for o in ornt), 'each element of ornt must be either -1 or 1 for 3D interfaces' + if (not isinstance(orientation, (tuple, list, Tuple)) or + len(orientation) != 3): + raise TypeError( + 'a 3D interface orientation must be ' + '(flag, sign1, sign2)') + try: + orientation = tuple( + _as_integer(value, 'a 3D interface orientation value') + for value in orientation) + except TypeError as error: + message = 'a 3D interface orientation must contain integers' + raise TypeError(message) from error + if any(value not in (-1, 1) for value in orientation): + raise ValueError( + 'each 3D interface orientation value must be +1 or -1') + orientation_arg = Tuple(*orientation) else: - raise ValueError(f'Unsupported dimension: {ldim}') - - # Strong requirement: the two boundaries must be defined over the same axis - # TODO [YG 10.02.2026]: relax this requirement ASAP - assert bnd_minus.axis == bnd_plus.axis + raise ValueError(f'unsupported interface dimension: {ldim}') - obj = new_basic(cls, name, bnd_minus, bnd_plus, ornt, + obj = new_basic(cls, name, bnd_minus, bnd_plus, orientation_arg, mapping, logical_domain) obj._mapping = mapping obj._logical_domain = logical_domain @@ -568,17 +697,44 @@ def name(self): def minus(self): return self.args[1] + @property + def minus_side(self): + """Descriptive alias of :attr:`minus`.""" + return self.minus + @property def plus(self): return self.args[2] @property - def ornt(self): + def plus_side(self): + """Descriptive alias of :attr:`plus`.""" + return self.plus + + @property + def orientation(self): + """Compact dimension-specific orientation of this interface.""" + if self.dim == 1: + return None + if self.dim == 2: + return int(self.args[3]) return restore_container_argument(self.args[3]) @property - def axis(self): - return self.plus.axis + def axis_map(self): + """Tangential-axis correspondence from the minus to the plus side.""" + return _interface_axis_map( + self.dim, self.minus.axis, self.plus.axis, self.orientation) + + @property + def patch_dim(self): + """Logical dimension of each adjacent patch.""" + return self.dim + + @property + def intrinsic_dim(self): + """Intrinsic dimension of the interface.""" + return self.dim - 1 @property def mapping(self): @@ -592,20 +748,6 @@ def _sympystr(self, printer): sstr = printer.doprint return '{}'.format(sstr(self.name)) -#============================================================================== -# TODO [MCP 16.02.2026]: remove this class as it is not used anywhere, neither -# in SymPDE nor in PSYDAC. It appears to have been replaced by Boundary, which -# is more general as it can represent geometrical objects in any dimension. -class Edge(object): - def __init__(self, name): - self._name = name - - @property - def name(self): - return self._name - - def __lt__(self, other): - return self.name.__lt__(other.name) #============================================================================== class Connectivity(abc.Mapping): @@ -638,7 +780,11 @@ def todict(self): connectivity = {} data = dict(sorted(self._data.items())) for name, v in data.items(): - connectivity[name] = [v.minus.todict(), v.plus.todict(), v.ornt] + orientation = v.orientation + if isinstance(orientation, tuple): + orientation = list(orientation) + connectivity[name] = [ + v.minus.todict(), v.plus.todict(), orientation] connectivity = dict(sorted(connectivity.items())) return connectivity diff --git a/sympde/topology/domain.py b/sympde/topology/domain.py index 6ff0affe..622f9f9e 100644 --- a/sympde/topology/domain.py +++ b/sympde/topology/domain.py @@ -4,6 +4,7 @@ import os from collections import OrderedDict +from itertools import product from typing import Union as TypeUnion, Optional, List, Iterable, TYPE_CHECKING # Union clashes with core.basic.Union @@ -25,8 +26,10 @@ restore_string_argument, ) from .basic import BasicDomain, InteriorDomain, Boundary, Union, Connectivity -from .basic import Interval, Interface, CornerBoundary, CornerInterface +from .basic import Interval, Interface +from .basic import CornerBoundary, CornerInterface from .basic import ProductDomain +from .basic import _as_integer # TODO fix circular dependency between domain and mapping if TYPE_CHECKING: @@ -34,6 +37,16 @@ iterable_types = (tuple, list, Tuple, Union) +#============================================================================== +def _serialized_integer(value, name): + """Read legacy decimal strings without accepting numeric floats.""" + if isinstance(value, str): + try: + return int(value) + except ValueError as error: + raise TypeError(f'{name} must be an integer') from error + return _as_integer(value, name) + #============================================================================== class Domain(BasicDomain): """ @@ -141,7 +154,14 @@ def __new__(cls, name : str, if not all([isinstance(i, Interface) for i in connectivity]): raise TypeError('> all interfaces must be of type Interface') - connectivity = Connectivity({str(i.name): i for i in connectivity}) + names = [str(i.name) for i in connectivity] + if len(names) != len(set(names)): + raise ValueError('> interface names must be unique') + connectivity = Connectivity(dict(zip(names, connectivity))) + + # Order the interfaces by name, so that the domain does not depend on + # the order in which they are given + connectivity = Connectivity(dict(sorted(connectivity.items()))) # TODO check that patches appearing in connectivity are in interiors @@ -197,13 +217,25 @@ def mapping(self) -> Optional[Mapping]: @property def subdomains(self) -> tuple: - """returns subdomains as tuple of Domains""" + """Return top-dimensional interior regions as a tuple.""" if isinstance( self.interior, iterable_types): subs = self.interior else: subs = [self.interior] return tuple(subs) + @property + def patches(self) -> tuple: + """Top-dimensional patch regions as a stable tuple. + + This additive multipatch view leaves the symbolic ``interior`` and + ``subdomains`` APIs unchanged. For a single-patch domain it contains + exactly that domain's interior. + """ + if isinstance(self.interior, InteriorDomain): + return (self.interior,) + return self.interior.as_tuple() + @property def mappings(self) -> OrderedDict: return OrderedDict([(P.logical_domain, P.mapping) @@ -219,6 +251,25 @@ def connectivity(self) -> Connectivity: """Contains information about the interfaces""" return self._connectivity + @property + def interface_map(self) -> Connectivity: + """Interfaces keyed by their unique names. + + Unlike the symbolic :attr:`interfaces` property, this view always has + mapping semantics, including when the domain has zero or one interface. + """ + return self.connectivity + + @property + def exterior_sides(self) -> tuple: + """Exterior boundary sides as a stable tuple.""" + boundary = self.boundary + if boundary is None: + return () + if isinstance(boundary, Boundary): + return (boundary,) + return boundary.as_tuple() + @property def dim(self) -> int: """Dimension of the space""" @@ -248,6 +299,16 @@ def corners(self): self._corners = corners return corners + @property + def shared_vertices(self) -> tuple: + """Shared patch-vertex equivalence classes as a stable tuple.""" + corners = self.corners + if corners is None: + return () + if isinstance(corners, CornerInterface): + return (corners,) + return corners.as_tuple() + def __len__(self): if isinstance(self.interior, InteriorDomain): return 1 @@ -291,8 +352,18 @@ def get_boundary(self, axis, ext): The domain boundary of interest. """ if axis is None: - assert(self.interior.dim == 1) + if self.interior.dim != 1: + raise ValueError('axis may be None only for a 1D domain') axis = 0 + else: + axis = _as_integer(axis, 'boundary axis') + ext = _as_integer(ext, 'boundary extremity') + + if not 0 <= axis < self.dim: + raise ValueError( + f'boundary axis must be between 0 and {self.dim - 1}') + if ext not in (-1, 1): + raise ValueError('boundary extremity must be either -1 or 1') if isinstance(self.boundary, Union): x = [i for i in self.boundary.args if i.ext == ext and i.axis == axis] @@ -341,7 +412,7 @@ def todict(self): raise ValueError(f'all mappings must have the same physical dimension, got {sorted(pdims)}') pdim = str(pdims.pop()) if pdims else ldim interior = self.interior.todict() - boundary = self.boundary.todict() + boundary = [side.todict() for side in self.exterior_sides] connectivity = self.connectivity.todict() dtype = self.dtype @@ -366,17 +437,12 @@ def export( self, filename ): yml = self.todict() # Dump metadata to string in YAML file format - geo = yaml.dump( data = yml, - sort_keys = None) + geo = yaml.safe_dump(data=yml, sort_keys=None) # Create HDF5 file (in parallel mode if MPI communicator size > 1) - h5 = h5py.File( filename, mode='w' ) - - # Write geometry metadata as fixed-length array of ASCII characters - h5['topology.yml'] = np.array( geo, dtype='S' ) - - # Close HDF5 file - h5.close() + with h5py.File(filename, mode='w') as h5: + # Write geometry metadata as fixed-length array of ASCII characters + h5['topology.yml'] = np.array(geo, dtype='S') @classmethod def from_file(cls, filename): @@ -405,8 +471,8 @@ def from_file(cls, filename): from sympde.topology.mapping import Mapping - h5 = h5py.File( filename, mode='r' ) - yml = yaml.load( h5['topology.yml'][()], Loader=yaml.SafeLoader ) + with h5py.File(filename, mode='r') as h5: + yml = yaml.safe_load(h5['topology.yml'][()]) if 'dim' in yml: raise ValueError(f"'{filename}' was written by an older version of SymPDE:" @@ -417,17 +483,42 @@ def from_file(cls, filename): pdim = int(yml['pdim']) dtype = yml['dtype'] d_interior = yml['interior'] - d_boundary = yml['boundary'] - d_connectivity = yml['connectivity'] + has_boundary_metadata = ( + 'boundary' in yml and yml['boundary'] is not None) + d_boundary = yml.get('boundary', []) + d_connectivity = yml.get('connectivity', {}) + + # Boundary metadata was historically cardinality-dependent: a single + # side was stored as a dictionary and several sides as a list. Accept + # that representation while making the canonical format always a list. + if d_boundary is None: + d_boundary = [] + elif isinstance(d_boundary, dict): + d_boundary = [d_boundary] + elif not isinstance(d_boundary, list): + raise TypeError('boundary topology must be a list or a dictionary') + + if d_connectivity is None: + d_connectivity = {} + elif not isinstance(d_connectivity, dict): + raise TypeError('connectivity topology must be a dictionary') if dtype == 'None': dtype = None - - assert dtype is not None - assert all(dtype) if isinstance(d_interior, dict): d_interior = [d_interior] dtype = [dtype] + if dtype is None: + dtype = [None] * len(d_interior) + + # Patches without a dtype are topological InteriorDomain objects, + # whose boundaries cannot be derived from a parametric box + if not any(dtype): + return cls._from_topological_dict( + domain_name, ldim, d_interior, d_boundary, d_connectivity) + if not all(dtype): + raise ValueError( + 'either all or none of the serialized patches must have a dtype') constructors = [globals()[dt['type']] for dt in dtype] interiors = [cs(i['name'], **dt['parameters']) for cs,i,dt in zip(constructors, d_interior, dtype)] @@ -435,51 +526,177 @@ def from_file(cls, filename): domains = [mapping(i) if mapping else i for i,mapping in zip(interiors, mappings)] patch_index = {I.name:ind for ind,I in enumerate(interiors)} - boundaries = [] - for bd in d_boundary: - name = bd['patch'] - axis = bd['axis'] - ext = bd['ext'] - i = patch_index[name] - bd = domains[i].get_boundary(axis=int(axis), ext=int(ext)) - boundaries.append(bd) + stored_boundary_sides = [] + for boundary_data in d_boundary: + if not isinstance(boundary_data, dict): + raise TypeError('each serialized boundary must be a dictionary') + name = boundary_data['patch'] + axis = _serialized_integer( + boundary_data['axis'], 'serialized boundary axis') + ext = _serialized_integer( + boundary_data['ext'], 'serialized boundary extremity') + if name not in patch_index: + raise ValueError( + f'serialized boundary references unknown patch: {name}') + domains[patch_index[name]].get_boundary(axis=axis, ext=ext) + stored_boundary_sides.append((str(name), axis, ext)) + + if len(stored_boundary_sides) != len(set(stored_boundary_sides)): + raise ValueError('serialized boundary contains duplicate sides') connectivity = [] - for _,(minus, plus, ornt) in d_connectivity.items(): + for interface_name, interface_data in d_connectivity.items(): + if not isinstance(interface_name, str): + raise TypeError('serialized interface names must be strings') + if len(interface_data) != 3: + raise ValueError( + 'multipatch topology files must store two interface sides ' + 'and an explicit orientation; regenerate legacy files ' + 'that contain only two interface sides') + minus, plus, orientation = interface_data + minus_name = minus['patch'] - minus_axis = int(minus['axis']) - minus_ext = int(minus['ext']) + minus_axis = _serialized_integer( + minus['axis'], 'serialized interface axis') + minus_ext = _serialized_integer( + minus['ext'], 'serialized interface extremity') minus_patch_i = patch_index[minus_name] plus_name = plus['patch'] - plus_axis = int(plus['axis']) - plus_ext = int(plus['ext']) + plus_axis = _serialized_integer( + plus['axis'], 'serialized interface axis') + plus_ext = _serialized_integer( + plus['ext'], 'serialized interface extremity') plus_patch_i = patch_index[plus_name] interface = ((minus_patch_i, minus_axis, minus_ext), - ( plus_patch_i, plus_axis, plus_ext), ornt) + (plus_patch_i, plus_axis, plus_ext), orientation, + interface_name) connectivity.append(interface) - if len(domains) == 1: - return domains[0] + if len(domains) == 1 and not connectivity: + domain = domains[0] + else: + domain = Domain.join(domains, connectivity, domain_name) + + derived_boundary_sides = [] + for boundary in domain.exterior_sides: + boundary_data = boundary.todict() + derived_boundary_sides.append(( + str(boundary_data['patch']), + int(boundary_data['axis']), + int(boundary_data['ext']))) + + if has_boundary_metadata and \ + set(stored_boundary_sides) != set(derived_boundary_sides): + raise ValueError( + 'serialized boundary does not match the exterior boundary ' + 'derived from connectivity') + + return domain + + @classmethod + def _from_topological_dict(cls, name, dim, d_interior, d_boundary, + d_connectivity): + """ + Create a domain whose patches are `InteriorDomain` objects without a + dtype, as read by `from_file`. + + The exterior boundaries cannot be derived from the patches, so they are + restored as they were stored, with their names, and with their axes + and extremities when these are known. + """ + if any(I.get('mapping', 'None') != 'None' for I in d_interior): + raise NotImplementedError( + 'reading mapped patches without a dtype is not supported') + + interiors = {I['name']: InteriorDomain(I['name'], dim=dim) + for I in d_interior} + + def boundary(data, kind): + if not isinstance(data, dict): + raise TypeError(f'each serialized {kind} side must be a dictionary') + patch = data['patch'] + if patch not in interiors: + raise ValueError( + f'serialized {kind} references unknown patch: {patch}') + + def optional_integer(key, label): + value = data[key] + if value in (None, 'None'): + return None + return _serialized_integer(value, f'serialized {kind} {label}') - return Domain.join(domains, connectivity, domain_name) + axis = optional_integer('axis', 'axis') + ext = optional_integer('ext', 'extremity') + return Boundary(data['name'], interiors[patch], axis=axis, ext=ext) + + boundaries = [boundary(data, 'boundary') for data in d_boundary] + + interfaces = [] + for interface_name, interface_data in d_connectivity.items(): + if not isinstance(interface_name, str): + raise TypeError('serialized interface names must be strings') + if len(interface_data) != 3: + raise ValueError( + 'multipatch topology files must store two interface sides ' + 'and an explicit orientation; regenerate legacy files ' + 'that contain only two interface sides') + minus, plus, orientation = interface_data + if isinstance(orientation, list): + orientation = tuple(orientation) + interfaces.append(Interface( + interface_name, boundary(minus, 'interface'), + boundary(plus, 'interface'), orientation)) + + return Domain(name, interiors=list(interiors.values()), + boundaries=boundaries, connectivity=interfaces) @classmethod - def join(cls, patches, connectivity, name): - """Create a multipatch domain by joining patches in 2D or 3D. + def join(cls, patches, connectivity=None, name=None, *, interfaces=None): + """ + Create a domain by joining one or more patches in 1D, 2D, or 3D. + + A single patch is returned unchanged when ``connectivity`` is empty. + With non-empty connectivity, its distinct boundaries may be joined to + create self-interfaces such as periodic identifications. Parameters ---------- patches : sequence of Domain - Atomic patches in the joined domain. - - connectivity : sequence of tuple - Interface descriptions of the form `(minus, plus, orientation)`. - Each side is `(patch, axis, ext)`, where `patch` is a patch - object or its index in `patches` and `ext` is `-1` or `1`. - A 2D orientation is `-1` or `1`. A 3D orientation is a tuple - of three values, each equal to `-1` or `1`. + Ordered non-empty collection of atomic, unconnected patches. All + patches must have the same logical dimension and the same + physical dimension, which is the logical dimension for a patch + without a mapping. Either all patches or none have a mapping. + Their positions define the integer patch indices accepted by an + interface side. + Patch objects, interior names, and serialized logical names must be + unique. Joined multipatch domains must be flattened before being + passed to this method. + + connectivity : sequence, optional + Interface descriptions. Each description is the tuple + `(minus, plus, orientation)` or + `(minus, plus, orientation, name)`. The optional name is + preserved exactly and must be unique among all interfaces. + A face may be part of at most one interface, and an interface + cannot join a face with itself. + If both this argument and `interfaces` are omitted, the + connectivity is empty. + + Each side is `(patch, axis, ext)`. `patch` is either a patch + object or its integer position in `patches`. `axis` is the + zero-based logical axis normal to the face. `ext=-1` selects the + lower-coordinate face and `ext=+1` the upper-coordinate face. + + `minus` and `plus` define the direction in which orientation + is read. In 1D, orientation is `None`. In 2D, it is `+1` or + `-1`. In 3D, it is `(flag, sign1, sign2)` as documented by + :class:`Interface`. + + interfaces : sequence, optional + Keyword alias for `connectivity`. Provide exactly one of + `connectivity` and `interfaces`. name : str Name of the domain. @@ -491,106 +708,241 @@ def join(cls, patches, connectivity, name): Notes ----- - The orientations are specified in the same manner as in GeoPDES, see e.g. - - and - T. Dokken, E. Quak, V. Skytt. Requirements from Isogeometric Analysis for changes in product design ontologies, 2010. - - Example - ------- - .. code-block:: python - - # list of patches (mapped domains) - Omega_0 = F0(A) - Omega_1 = F1(A) - Omega_2 = F2(A) - Omega_3 = F3(A) - - patches = [Omega_0, Omega_1, Omega_2, Omega_3] - - # integers representing the axes - axis_0 = 0 - axis_1 = 1 - axis_2 = 2 - - # integers representing the extremities: left (-1) or right (+1) - ext_0 = -1 - ext_1 = +1 - - # A connectivity list in 2D - connectivity = [((Omega_0, axis_0, ext_0), (Omega_1, axis_0, ext_1), 1), - ((Omega_1, axis_1, ext_0), (Omega_3, axis_1, ext_1), -1), - ((Omega_0, axis_1, ext_0), (Omega_2, axis_1, ext_1), 1), - ((Omega_2, axis_0, ext_0), (Omega_3, axis_0, ext_1), -1)] - - # alternative option (passing interface patches by their indices in the patches list): - connectivity = [((0, axis_0, ext_0), (1, axis_0, ext_1), 1), - ((1, axis_1, ext_0), (3, axis_1, ext_1), -1), - ((0, axis_1, ext_0), (2, axis_1, ext_1), 1), - ((2, axis_0, ext_0), (3, axis_0, ext_1), -1)] - - # A connectivity list in 3D - connectivity = [((Omega_0, axis_0, ext_1), (Omega_1, axis_0, ext_0), ( 1, 1, 1)), - ((Omega_0, axis_1, ext_1), (Omega_2, axis_1, ext_0), ( 1, -1, 1)), - ((Omega_1, axis_1, ext_1), (Omega_3, axis_1, ext_0), (-1, 1, -1)), - ((Omega_2, axis_0, ext_1), (Omega_3, axis_0, ext_0), (-1, 1, 1))] - - # alternative option (passing interface patches by their indices in the patches list): - connectivity = [((0, axis_0, ext_1), (1, axis_0, ext_0), ( 1, 1, 1)), - ((0, axis_1, ext_1), (2, axis_1, ext_0), ( 1, -1, 1)), - ((1, axis_1, ext_1), (3, axis_1, ext_0), (-1, 1, -1)), - ((2, axis_0, ext_1), (3, axis_0, ext_0), (-1, 1, 1))] - - # the multi-patch domain - Omega = Domain.join(patches=patches, connectivity=connectivity, name='Omega') + The compact orientation convention follows the convention used by + GeoPDEs; see its `multipatch geometry specification + `_ + and T. Dokken, E. Quak, V. Skytt, *Requirements from Isogeometric + Analysis for Changes in Product Design Ontologies* (2010). + + Connectivity is symbolic. This method does not infer interfaces from + physical coordinates or verify that mapped faces coincide. It also + does not impose a mesh, spline space, or trace-coupling policy. + + Examples + -------- + A four-patch 2D connectivity may reference patch objects directly: + + >>> from sympde.api import Cube, Domain, Square + >>> patches = [Square(f'P{i}') for i in range(4)] + >>> connectivity = [ + ... ((patches[0], 0, -1), (patches[1], 0, +1), +1), + ... ((patches[1], 1, -1), (patches[3], 1, +1), -1), + ... ((patches[0], 1, -1), (patches[2], 1, +1), +1), + ... ((patches[2], 0, -1), (patches[3], 0, +1), -1), + ... ] + >>> omega = Domain.join(patches, connectivity, name='Omega') + + The equivalent connectivity may use positions in ``patches``: + + >>> connectivity = [ + ... ((0, 0, -1), (1, 0, +1), +1), + ... ((1, 1, -1), (3, 1, +1), -1), + ... ((0, 1, -1), (2, 1, +1), +1), + ... ((2, 0, -1), (3, 0, +1), -1), + ... ] + >>> omega = Domain.join(patches, connectivity, name='Omega') + + A tuple also describes a cross-axis 3D interface: + + >>> A = Cube('A') + >>> B = Cube('B') + >>> omega = Domain.join( + ... patches=[A, B], + ... interfaces=[( + ... (A, 0, +1), + ... (B, 1, -1), + ... (-1, +1, -1), + ... )], + ... name='Omega3D', + ... ) + >>> interface, = omega.interface_map.values() + >>> interface.axis_map + ((1, 2, 1), (2, 0, -1)) """ - assert isinstance(patches, (tuple, list)) - assert isinstance(connectivity, (tuple, list)) - assert isinstance(name, str) + if not isinstance(patches, (tuple, list)): + raise TypeError('patches must be a list or tuple') + if not patches or not all(isinstance(patch, Domain) for patch in patches): + raise TypeError('patches must contain Domain objects') + + if len({id(patch) for patch in patches}) != len(patches): + raise ValueError('patches must not contain the same object twice') + + if any(not isinstance(patch.interior, InteriorDomain) or + patch.connectivity for patch in patches): + raise ValueError( + 'Domain.join expects atomic patch domains without existing ' + 'connectivity; flatten joined domains before joining them') + + patch_names = [str(patch.interior.name) for patch in patches] + if len(set(patch_names)) != len(patch_names): + raise ValueError('patch interior names must be unique') + + serialized_patch_names = [ + str(patch.interior.logical_domain.name) + if patch.interior.logical_domain is not None + else str(patch.interior.name) + for patch in patches + ] + if len(set(serialized_patch_names)) != len(serialized_patch_names): + raise ValueError( + 'logical patch names used for serialization must be unique') + + if interfaces is not None: + if connectivity is not None: + raise TypeError( + 'provide either connectivity or interfaces, not both') + connectivity = interfaces + elif connectivity is None: + connectivity = () + if not isinstance(connectivity, (tuple, list)): + raise TypeError('connectivity/interfaces must be a list or tuple') + if not isinstance(name, str): + raise TypeError('name must be a string') - if len(patches) == 1: - # single patch domain: return the patch - assert len(connectivity) == 0 + if len(patches) == 1 and not connectivity: + # Preserve the original single-patch domain when there is no + # topology to add. A non-empty connectivity may contain a + # self-interface and must follow the normal joining path. return patches[0] - assert all(p.dim==patches[0].dim for p in patches) + if not all(p.dim == patches[0].dim for p in patches): + raise ValueError('all patches must have the same logical dimension') + pdims = {p.mapping.pdim if p.mapping is not None else p.dim for p in patches} + if len(pdims) > 1: + raise ValueError('all patches must have the same physical dimension') + if len({p.mapping is None for p in patches}) > 1: + raise ValueError('either all patches or none must have a mapping') ldim = int(patches[0].dim) - patch_given_by_indices = (len(connectivity) > 0 and isinstance(connectivity[0][0][0], int)) + normalized_connectivity = [] + explicit_interface_names = set() + for interface_data in connectivity: + if not isinstance(interface_data, tuple) or \ + len(interface_data) not in (3, 4): + raise TypeError( + 'an interface must be the tuple ' + '(minus, plus, orientation) or ' + '(minus, plus, orientation, name)') + + minus_spec, plus_spec, orientation = interface_data[:3] + interface_name = ( + interface_data[3] if len(interface_data) == 4 else None) + if interface_name is not None: + if not isinstance(interface_name, str): + raise TypeError('an explicit interface name must be a string') + if not interface_name: + raise ValueError('an explicit interface name cannot be empty') + if interface_name in explicit_interface_names: + raise ValueError( + f'duplicate explicit interface name: {interface_name}') + explicit_interface_names.add(interface_name) + + normalized_connectivity.append( + (minus_spec, plus_spec, orientation, interface_name)) from sympde.topology.mapping import MultiPatchMapping # ... connectivity interfaces = {} boundaries = [] - for cn in connectivity: - if patch_given_by_indices: - patch_minus = patches[cn[0][0]] - patch_plus = patches[cn[1][0]] + # Reserve every explicit name before allocating generated ones. This + # makes explicit identity independent of the interface input order. + interface_names = set(explicit_interface_names) + + def get_unique_interface_name(base, used_names): + """Return an unused name while preserving the historical pattern.""" + if base not in used_names: + used_names.add(base) + return base + + occurrence = 2 + while f'{base}#{occurrence}' in used_names: + occurrence += 1 + name = f'{base}#{occurrence}' + used_names.add(name) + return name + + def logical_patch_name(boundary): + patch = boundary.logical_domain or boundary + return patch.domain.name + + def get_boundary(boundary_spec): + if not isinstance(boundary_spec, tuple) or len(boundary_spec) != 3: + raise TypeError( + 'an interface side must be the tuple ' + '(patch, axis, extremity)') + patch_ref, axis, ext = boundary_spec + if isinstance(patch_ref, bool): + raise TypeError( + 'a patch reference must be a patch or integer index') + if isinstance(patch_ref, Domain): + matching_patches = [ + patch for patch in patches if patch_ref is patch] + if not matching_patches: + raise ValueError( + 'an interface references a patch not present in patches') + patch = matching_patches[0] else: - patch_minus = cn[0][0] - patch_plus = cn[1][0] - bnd_minus = patch_minus.get_boundary(axis=cn[0][1], ext=cn[0][2]) - bnd_plus = patch_plus.get_boundary(axis=cn[1][1], ext=cn[1][2]) - - # Check that orientation is provided in the correct format depending on the dimension - ornt = cn[2] if len(cn) == 3 else None - if ldim == 1: - assert ornt is None, 'ornt is not needed for 1D interfaces' - elif ldim == 2: - assert ornt in (-1, 1), 'ornt must be either -1 or 1 for 2D interfaces' - elif ldim == 3: - ornt = tuple(ornt) - assert len(ornt) == 3, 'ornt must be a tuple of length 3 for 3D interfaces' - assert all(o in (-1, 1) for o in ornt), 'each element of ornt must be either -1 or 1 for 3D interfaces' + try: + patch_index = _as_integer(patch_ref, 'patch index') + except TypeError as error: + raise TypeError( + 'a patch reference must be a patch or integer index') \ + from error + if not 0 <= patch_index < len(patches): + raise IndexError(f'patch index {patch_index} is out of range') + patch = patches[patch_index] + axis = _as_integer(axis, 'interface side axis') + ext = _as_integer(ext, 'interface side extremity') + if not 0 <= axis < ldim: + raise ValueError(f'axis must be between 0 and {ldim - 1}') + if ext not in (-1, 1): + raise ValueError('boundary extremity must be either -1 or 1') + return patch.get_boundary(axis=axis, ext=ext) + + sides = [] + joined_boundaries = set() + for minus_spec, plus_spec, orientation, explicit_name \ + in normalized_connectivity: + + bnd_minus = get_boundary(minus_spec) + bnd_plus = get_boundary(plus_spec) + # A face is part of at most one interface, and it is not joined + # with itself + if bnd_minus == bnd_plus: + raise ValueError( + f'an interface joins the face {bnd_minus} with itself') + for bnd in (bnd_minus, bnd_plus): + if bnd in joined_boundaries: + raise ValueError( + f'the face {bnd} is part of more than one interface') + joined_boundaries.add(bnd) + sides.append((bnd_minus, bnd_plus, orientation, explicit_name)) + + def side_key(side): + bnd_minus, bnd_plus = side[:2] + return (logical_patch_name(bnd_minus), bnd_minus.axis, bnd_minus.ext, + logical_patch_name(bnd_plus), bnd_plus.axis, bnd_plus.ext) + + # Generate the missing names in the order of the faces, so that they + # do not depend on the order of the connectivity. A mapped interface + # and its logical interface share one name, generated from the + # logical patch names. This keeps the name of the logical interface + # when the domain is read from a file. + for bnd_minus, bnd_plus, orientation, explicit_name \ + in sorted(sides, key=side_key): + if explicit_name is None: + base = (f'{logical_patch_name(bnd_minus)}|' + f'{logical_patch_name(bnd_plus)}') + interface_name = get_unique_interface_name( + base, interface_names) else: - raise ValueError(f'Unsupported dimension: {ldim}') + interface_name = explicit_name - # Create a new Interface object using the join method of the Boundary objects - interface = bnd_minus.join(bnd_plus, ornt=ornt) - if interface.name in interfaces: - interface = bnd_plus.join(bnd_minus, ornt=ornt) + interface = bnd_minus.join( + bnd_plus, orientation, name=interface_name, + logical_name=interface_name) - interfaces[interface.name] = interface + interfaces[str(interface.name)] = interface boundaries.append(bnd_minus) boundaries.append(bnd_plus) @@ -602,20 +954,27 @@ def join(cls, patches, connectivity, name): boundaries = Union(*[b for p in patches for b in p.boundary]).complement(Union(*boundaries)) if boundaries is None: boundaries = () - else : + elif isinstance(boundaries, Boundary): + boundaries = (boundaries,) + else: boundaries = boundaries.as_tuple() # ... interiors - interiors = Union(*[p.interior for p in patches]) + interior_patches = [p.interior for p in patches] + interiors = Union(*interior_patches) - if all(e.mapping for e in interiors): - logical_interiors = Union(*[e.logical_domain for e in interiors]) + if all(e.mapping for e in interior_patches): + logical_interiors = Union(*[e.logical_domain for e in interior_patches]) logical_boundaries = [e.logical_domain for e in boundaries] logical_connectivity = Connectivity() for k,v in connectivity.items(): logical_connectivity[v.logical_domain.name] = v.logical_domain - mapping = MultiPatchMapping({e.logical_domain: e.mapping for e in interiors}) + patch_mappings = { + e.logical_domain: e.mapping for e in interior_patches} + mapping = ( + interior_patches[0].mapping if len(interior_patches) == 1 + else MultiPatchMapping(patch_mappings)) logical_domain = Domain(name, interiors=logical_interiors, boundaries=logical_boundaries, @@ -632,162 +991,198 @@ def join(cls, patches, connectivity, name): logical_domain=logical_domain) def get_shared_corners(self): - """ Compute the corners shared by multiple patches in 2D """ - - interfaces = self.interfaces - interfaces = (interfaces,) if isinstance(interfaces, Interface) else interfaces - - directions = {i.plus:i.ornt for i in interfaces} - directions.update({i.minus:i.ornt for i in interfaces}) - - boundaries = {i.minus:i.plus for i in interfaces} - boundaries.update({value:key for key, value in boundaries.items()}) + """Return equivalence classes of patch vertices joined by interfaces. - not_treated_corners = set([tuple(set((b, n))) for b in boundaries for n in b.adjacent_boundaries]) - grouped_corners = [] - - while not_treated_corners: - corner = not_treated_corners.pop() - grouped_corners.append([corner]) - if not ( corner[0] in boundaries and corner[1] in boundaries): - while corner[1] in boundaries: - bd1 = boundaries[corner[1]] - bd2 = bd1.domain.get_boundary(axis=corner[0].axis, ext=corner[0].ext) - corner = (bd1, bd2.rotate(directions[bd1])) - grouped_corners[-1].append(corner) - - corner = grouped_corners[-1][0] - while corner[0] in boundaries: - bd2 = boundaries[corner[0]] - bd1 = bd2.domain.get_boundary(axis=corner[1].axis, ext=corner[1].ext) - corner = (bd1.rotate(directions[bd2]), bd2) - grouped_corners[-1].insert(0, corner) - - else: - while corner[1] in boundaries: - bd1 = boundaries[corner[1]] - bd2 = bd1.domain.get_boundary(axis=corner[0].axis, ext=corner[0].ext) - corner = (bd1, bd2.rotate(directions[bd1])) - if corner == grouped_corners[-1][0]: - break - grouped_corners[-1].append(corner) - else: - corner = grouped_corners[-1][0] - while corner[0] in boundaries: - bd2 = boundaries[corner[0]] - bd1 = bd2.domain.get_boundary(axis=corner[1].axis, ext=corner[1].ext) - corner = (bd1.rotate(directions[bd2]), bd2) - grouped_corners[-1].insert(0, corner) - - grouped_corners[-1] = tuple(tuple(set(c)) for c in grouped_corners[-1]) - not_treated_corners = not_treated_corners.difference(grouped_corners[-1]) - - grouped_corners = set(tuple(grouped_corners)) - grouped_corners = Union(*[CornerInterface(*[CornerBoundary(*e) for e in cs]) for cs in grouped_corners]) - return grouped_corners + Face vertices are related using each interface's signed tangential-axis + permutation. Connected components then handle any patch incidence in + both 2D and 3D, without walking around interfaces in a prescribed order. + """ + interfaces = self.interfaces + if interfaces is None: + return None + interfaces = (interfaces,) if isinstance(interfaces, Interface) else tuple(interfaces) + + # The vertices of a patch are the corners of its parametric box. A + # topological InteriorDomain, as read from a file without dtype, has + # no such box. + for interface in interfaces: + for side in (interface.minus, interface.plus): + if not hasattr(side.domain, 'get_boundary'): + raise ValueError( + 'shared vertices are only defined for patches with a ' + f'parametric box, such as Square or Cube; patch ' + f'{side.domain.name} is a topological InteriorDomain') + + parent = {} + + def find(vertex): + parent.setdefault(vertex, vertex) + while parent[vertex] != vertex: + parent[vertex] = parent[parent[vertex]] + vertex = parent[vertex] + return vertex + + def union(vertex_1, vertex_2): + root_1 = find(vertex_1) + root_2 = find(vertex_2) + if root_1 != root_2: + if root_2.sort_key() < root_1.sort_key(): + root_1, root_2 = root_2, root_1 + parent[root_2] = root_1 + + def make_vertex(face, extents): + boundaries = [ + face.domain.get_boundary(axis=axis, ext=ext) + for axis, ext in enumerate(extents) + ] + return CornerBoundary(*boundaries) + + for interface in interfaces: + minus_axis = interface.minus.axis + plus_axis = interface.plus.axis + + for tangent_extents in product((-1, 1), repeat=self.dim - 1): + minus_extents = [None] * self.dim + plus_extents = [None] * self.dim + minus_extents[minus_axis] = interface.minus.ext + plus_extents[plus_axis] = interface.plus.ext + for source_position, (source_axis, target_axis, direction) \ + in enumerate(interface.axis_map): + source_ext = tangent_extents[source_position] + minus_extents[source_axis] = source_ext + plus_extents[target_axis] = direction * source_ext + + union( + make_vertex(interface.minus, minus_extents), + make_vertex(interface.plus, plus_extents)) + + groups = {} + for vertex in parent: + groups.setdefault(find(vertex), []).append(vertex) + + shared_corners = [ + CornerInterface(*vertices) + for _, vertices in sorted(groups.items(), key=lambda item: item[0].sort_key()) + if len(vertices) > 1 + ] + return Union(*shared_corners) def get_subdomain(self, names): """ - Returns an individual patch or a Union of patches of a multipatch domain. + Return the subdomain induced by a selection of patch names. Parameters ---------- names : tuple of str or str - Names of the patches to join. - If a string is given, the corresponding patch will be returned. - If a tuple of strings is given, the Union of the corresponding subdomains will be returned. + Names of the patches to retain. If a string is given, a one-patch + subdomain is returned. If a tuple is given, its order determines + the generated subdomain name. Notes ----- - The subdomain is returned as it was before being joined, which means that its boundary includes the - boundaries that are part of an interface in the multipatch domain. + Interfaces whose two sides belong to selected patches are retained, + including repeated and self-interfaces. If exactly one side is + selected, that side becomes part of the subdomain's exterior boundary. """ if names == (): return None if isinstance(names, str): names = (names,) - assert names[0] in self.interior_names - - elif isinstance(names, tuple): - assert all(isinstance(name, str) for name in names) - assert len(set(names)) == len(names) - assert all(name in self.interior_names or name == self.name for name in names) - - # Check trivial case of single patch domain - if isinstance(self.interior, InteriorDomain): - assert names[0] == self.interior.name + elif not isinstance(names, tuple): + raise TypeError('names must be a string or a tuple of strings') + + if not all(isinstance(name, str) for name in names): + raise TypeError('every subdomain name must be a string') + if len(set(names)) != len(names): + raise ValueError('subdomain names must be unique') + + interior_names = tuple(self.interior_names) + unknown_names = [ + name for name in names + if name not in interior_names and name != self.name + ] + if unknown_names: + unknown = ', '.join(unknown_names) + raise ValueError(f'unknown subdomain name(s): {unknown}') + + # Selecting the containing domain or every patch is an identity + # operation, preserving its name and object identity. + if self.name in names or set(names) == set(interior_names): return self - # If all patches are joined we get the full domain - # Same if the full domain is part of the union - if len(names) == len(self.interior_names) or self.name in names: - return self - - # Build dictionary of interiors accessed by names - interior_dict = {i.name: i for i in self.interior.as_tuple()} - - # Build dictionary of boundaries - if self.boundary is not None: - if isinstance(self.boundary, Union): - boundary_dict = {(b.domain.name, b.axis, b.ext): b for b in self.boundary.as_tuple()} + selected_names = set(names) + interior_dict = {interior.name: interior for interior in self.patches} + interiors = [interior_dict[name] for name in names] + + # Existing exterior sides remain exterior when their patch is selected. + boundaries = [ + boundary for boundary in self.exterior_sides + if boundary.domain.name in selected_names + ] + interfaces = OrderedDict() + + # Connectivity induces both the retained interfaces and the new cut + # boundary. Iterating the original mapping preserves interface keys and + # distinguishes repeated interfaces between the same patch pair. + for key, interface in self.connectivity.items(): + minus_selected = interface.minus.domain.name in selected_names + plus_selected = interface.plus.domain.name in selected_names + + if minus_selected and plus_selected: + interfaces[key] = interface + elif minus_selected: + boundaries.append(interface.minus) + elif plus_selected: + boundaries.append(interface.plus) + + connectivity = Connectivity(interfaces) + name = '|'.join(names) + mapping = None + logical_domain = None + + if all(interior.mapping is not None for interior in interiors): + logical_interiors = [ + interior.logical_domain for interior in interiors] + logical_boundaries = [ + boundary.logical_domain for boundary in boundaries] + logical_interfaces = OrderedDict() + + if any(interior is None for interior in logical_interiors) or \ + any(boundary is None for boundary in logical_boundaries): + raise ValueError( + 'mapped patch topology is missing logical-domain metadata') + + for interface in interfaces.values(): + logical_interface = interface.logical_domain + if logical_interface is None: + raise ValueError( + 'mapped interface is missing logical-domain metadata') + logical_interfaces[str(logical_interface.name)] = \ + logical_interface + + logical_domain = Domain( + name, + interiors=logical_interiors, + boundaries=logical_boundaries, + connectivity=Connectivity(logical_interfaces)) + + if len(interiors) == 1: + mapping = interiors[0].mapping else: - b = self.boundary - boundary_dict = {(b.domain.name, b.axis, b.ext): b} - else: - boundary_dict = {} - - # Build dictionary of interfaces - if self.interfaces is not None: - if isinstance(self.interfaces, Union): - interfaces_dict = {(i.minus.domain.name, i.plus.domain.name): i for i in self.interfaces.as_tuple()} - elif isinstance(self.interfaces, Interface): - i = self.interfaces - interfaces_dict = {(i.minus.domain.name, i.plus.domain.name): i} - else: - interfaces_dict = {} - - interfaces = [] - - for name in names: - if name == self.name: - return self - interior = interior_dict[name] - boundaries = [boundary_dict.get((name, axis, ext)) for axis in range(self.dim) for ext in [-1, 1]] - - boundaries = [b for b in boundaries if b is not None] - - # Extract boundaries and interfaces from interfaces_dict - for other_name in self.interior_names: - if other_name != name: - i_minus = interfaces_dict.pop((name, other_name), None) - i_plus = interfaces_dict.pop((other_name, name), None) - - if other_name not in names: - if i_minus is not None: - boundaries.append(i_minus.minus) - if i_plus is not None: - boundaries.append(i_plus.plus) - else: - if i_plus is not None: - interfaces.append((i_plus.name, i_plus)) - if i_minus is not None: - interfaces.append((i_minus.name, i_minus)) + from sympde.topology.mapping import MultiPatchMapping + mapping = MultiPatchMapping({ + interior.logical_domain: interior.mapping + for interior in interiors + }) - # Create domain with name, interior and boundaries - new_domain = Domain(name=name, interiors=interior, boundaries=boundaries, - mapping=interior.mapping, logical_domain=interior.logical_domain) - try: - previous_domain = Domain.join([previous_domain,new_domain], [], name=f"{previous_domain.name}|{new_domain.name}") - except NameError: - previous_domain = new_domain - - # Add interfaces - joined_domain = previous_domain - for k,v in interfaces: - joined_domain.connectivity[k] = v - - return joined_domain + return Domain( + name, + interiors=interiors, + boundaries=boundaries, + connectivity=connectivity, + mapping=mapping, + logical_domain=logical_domain) #============================================================================== @@ -897,12 +1292,22 @@ def boundary(self): def get_boundary(self, axis=None, ext=None): """return boundary by (axis, ext).""" - - assert(not( ext is None )) + if ext is None: + raise ValueError('boundary extremity must be provided') if axis is None: - assert(self.dim == 1) + if self.dim != 1: + raise ValueError('axis may be None only for a 1D domain') axis = 0 + else: + axis = _as_integer(axis, 'boundary axis') + ext = _as_integer(ext, 'boundary extremity') + + if not 0 <= axis < self.dim: + raise ValueError( + f'boundary axis must be between 0 and {self.dim - 1}') + if ext not in (-1, 1): + raise ValueError('boundary extremity must be either -1 or 1') if isinstance(self.boundary, Union): x = [i for i in self.boundary.args if i.ext == ext and i.axis==axis] diff --git a/sympde/topology/mapping.py b/sympde/topology/mapping.py index 103d9cd4..e8a4bc7c 100644 --- a/sympde/topology/mapping.py +++ b/sympde/topology/mapping.py @@ -628,20 +628,26 @@ def __new__(cls, mapping, logical_domain): elif boundaries: kwargs['boundaries'] = mapping(boundaries) - interfaces = logical_domain.connectivity.interfaces - if interfaces: - if isinstance(interfaces, Union): - interfaces = interfaces.args - else: - interfaces = [interfaces] - connectivity = {} - for e in interfaces: - connectivity[e.name] = Interface(e.name, mapping(e.minus), mapping(e.plus)) + if logical_domain.connectivity: + connectivity = { + name: mapping(interface) + for name, interface in logical_domain.connectivity.items() + } kwargs['connectivity'] = Connectivity(connectivity) name = '{}({})'.format(str(mapping.name), str(logical_domain.name)) return Domain(name, **kwargs) + elif isinstance(logical_domain, Interface): + minus = mapping(logical_domain.minus) + plus = mapping(logical_domain.plus) + interface_mapping = InterfaceMapping( + minus.mapping, plus.mapping) + return Interface( + str(logical_domain.name), minus, plus, + logical_domain.orientation, + mapping=interface_mapping, + logical_domain=logical_domain) elif isinstance(logical_domain, NCubeInterior): name = logical_domain.name dim = logical_domain.dim @@ -758,12 +764,12 @@ def eval(cls, F): Parameters ---------- - F: Mapping + F : Mapping mapping object Returns ------- - expr : ImmutableDenseMatrix + expr : ImmutableDenseMatrix the jacobian matrix """ @@ -806,15 +812,15 @@ def eval(cls, F, v): Parameters ---------- - F: Mapping + F : Mapping mapping object - v: + v : tuple | list | Tuple | ImmutableDenseMatrix | Matrix the basis function Returns ------- - expr : Tuple + expr : Tuple the covariant transformation """ @@ -855,15 +861,15 @@ def eval(cls, F, v): Parameters ---------- - F: Mapping + F : Mapping mapping object - v: + v : tuple | list | Tuple | ImmutableDenseMatrix | Matrix the basis function Returns ------- - expr : Tuple + expr : Tuple the contravariant transformation """ @@ -1228,7 +1234,11 @@ def eval(cls, expr, domain, **options): J = mapping.jacobian det = sqrt((J.T*J).det()) else: - axis = domain.axis + # Interface integration uses the minus side as its canonical + # parametrization. The plus face may have a different normal + # axis, so an Interface has deliberately no single ``axis``. + axis = (domain.minus.axis + if isinstance(domain, Interface) else domain.axis) J = JacobianSymbol(mapping, axis=axis) det = sqrt((J.T*J).det()) diff --git a/sympde/topology/multipatch_gallery.py b/sympde/topology/multipatch_gallery.py new file mode 100644 index 00000000..234620d8 --- /dev/null +++ b/sympde/topology/multipatch_gallery.py @@ -0,0 +1,1480 @@ +# --------------------------------------------------------------------------- # +# This file is part of SymPDE. # +# --------------------------------------------------------------------------- # +"""Gallery of symbolic 2D and 3D multipatch domains.""" + +from argparse import ArgumentParser +from pathlib import Path + +import numpy as np + +from .analytical_mapping import ( + AffineMapping, + IdentityMapping, + PolarMapping, + TorusMapping, + TransposedPolarMapping, +) +from .domain import Cube, Domain, Square + +__all__ = ( + 'build_annulus_3', + 'build_annulus_4', + 'build_cartesian_multipatch_domain_2d', + 'build_curved_l_shape', + 'build_multipatch_domain_2d', + 'build_multipatch_domain_3d', + 'build_oriented_two_cube', + 'build_pretzel', + 'build_pretzel_annulus', + 'build_pretzel_debug', + 'build_pretzel_f', + 'build_self_connected_patch', + 'build_square_2', + 'build_square_4', + 'build_square_6', + 'build_square_8', + 'build_square_9', + 'build_three_patch_shared_edge', + 'build_three_patch_shared_vertex', + 'build_torus_3d', + 'build_two_cube_3d', + 'build_two_patch_twisted_strip', + 'plot_multipatch_domain', +) + + +# ============================================================================= +# 2D domains +# ============================================================================= + +def build_square_2(): + """Build a square decomposed into two horizontal patches.""" + logical_patch_1 = Square( + 'P1', + bounds1=(0., np.pi), + bounds2=(0., np.pi / 2), + ) + logical_patch_2 = Square( + 'P2', + bounds1=(0., np.pi), + bounds2=(np.pi / 2, np.pi), + ) + patch_1 = IdentityMapping('M1', dim=2)(logical_patch_1) + patch_2 = IdentityMapping('M2', dim=2)(logical_patch_2) + + patches = (patch_1, patch_2) + connectivity = ( + ((patch_1, 1, +1), (patch_2, 1, -1), +1), + ) + return Domain.join(patches, connectivity, name='square_2') + + +def build_square_4(): + """Build a square decomposed into a two-by-two patch grid.""" + logical_patch_a = Square( + 'A', + bounds1=(0., np.pi / 2), + bounds2=(0., np.pi / 2), + ) + logical_patch_b = Square( + 'B', + bounds1=(np.pi / 2, np.pi), + bounds2=(0., np.pi / 2), + ) + logical_patch_c = Square( + 'C', + bounds1=(0., np.pi / 2), + bounds2=(np.pi / 2, np.pi), + ) + logical_patch_d = Square( + 'D', + bounds1=(np.pi / 2, np.pi), + bounds2=(np.pi / 2, np.pi), + ) + patch_a = IdentityMapping('M1', dim=2)(logical_patch_a) + patch_b = IdentityMapping('M2', dim=2)(logical_patch_b) + patch_c = IdentityMapping('M3', dim=2)(logical_patch_c) + patch_d = IdentityMapping('M4', dim=2)(logical_patch_d) + + patches = (patch_a, patch_b, patch_c, patch_d) + connectivity = ( + ((patch_a, 0, +1), (patch_b, 0, -1), +1), + ((patch_a, 1, +1), (patch_c, 1, -1), +1), + ((patch_c, 0, +1), (patch_d, 0, -1), +1), + ((patch_b, 1, +1), (patch_d, 1, -1), +1), + ) + return Domain.join(patches, connectivity, name='square_4') + + +def build_square_6(): + """Build a square decomposed into a two-by-three patch grid.""" + logical_patch_1 = Square( + 'P1', + bounds1=(0., np.pi / 2), + bounds2=(0., np.pi / 3), + ) + logical_patch_2 = Square( + 'P2', + bounds1=(np.pi / 2, np.pi), + bounds2=(0., np.pi / 3), + ) + logical_patch_3 = Square( + 'P3', + bounds1=(0., np.pi / 2), + bounds2=(np.pi / 3, 2 * np.pi / 3), + ) + logical_patch_4 = Square( + 'P4', + bounds1=(np.pi / 2, np.pi), + bounds2=(np.pi / 3, 2 * np.pi / 3), + ) + logical_patch_5 = Square( + 'P5', + bounds1=(0., np.pi / 2), + bounds2=(2 * np.pi / 3, np.pi), + ) + logical_patch_6 = Square( + 'P6', + bounds1=(np.pi / 2, np.pi), + bounds2=(2 * np.pi / 3, np.pi), + ) + patch_1 = IdentityMapping('M1', dim=2)(logical_patch_1) + patch_2 = IdentityMapping('M2', dim=2)(logical_patch_2) + patch_3 = IdentityMapping('M3', dim=2)(logical_patch_3) + patch_4 = IdentityMapping('M4', dim=2)(logical_patch_4) + patch_5 = IdentityMapping('M5', dim=2)(logical_patch_5) + patch_6 = IdentityMapping('M6', dim=2)(logical_patch_6) + + patches = (patch_1, patch_2, patch_3, patch_4, patch_5, patch_6) + connectivity = ( + ((patch_1, 0, +1), (patch_2, 0, -1), +1), + ((patch_3, 0, +1), (patch_4, 0, -1), +1), + ((patch_5, 0, +1), (patch_6, 0, -1), +1), + ((patch_1, 1, +1), (patch_3, 1, -1), +1), + ((patch_3, 1, +1), (patch_5, 1, -1), +1), + ((patch_2, 1, +1), (patch_4, 1, -1), +1), + ((patch_4, 1, +1), (patch_6, 1, -1), +1), + ) + return Domain.join(patches, connectivity, name='square_6') + + +def build_square_8(): + """Build a three-by-three square grid without its center patch.""" + logical_patch_1 = Square( + 'P1', + bounds1=(0., np.pi / 3), + bounds2=(0., np.pi / 3), + ) + logical_patch_2 = Square( + 'P2', + bounds1=(np.pi / 3, 2 * np.pi / 3), + bounds2=(0., np.pi / 3), + ) + logical_patch_3 = Square( + 'P3', + bounds1=(2 * np.pi / 3, np.pi), + bounds2=(0., np.pi / 3), + ) + logical_patch_4 = Square( + 'P4', + bounds1=(0., np.pi / 3), + bounds2=(np.pi / 3, 2 * np.pi / 3), + ) + logical_patch_5 = Square( + 'P5', + bounds1=(2 * np.pi / 3, np.pi), + bounds2=(np.pi / 3, 2 * np.pi / 3), + ) + logical_patch_6 = Square( + 'P6', + bounds1=(0., np.pi / 3), + bounds2=(2 * np.pi / 3, np.pi), + ) + logical_patch_7 = Square( + 'P7', + bounds1=(np.pi / 3, 2 * np.pi / 3), + bounds2=(2 * np.pi / 3, np.pi), + ) + logical_patch_8 = Square( + 'P8', + bounds1=(2 * np.pi / 3, np.pi), + bounds2=(2 * np.pi / 3, np.pi), + ) + patch_1 = IdentityMapping('M1', dim=2)(logical_patch_1) + patch_2 = IdentityMapping('M2', dim=2)(logical_patch_2) + patch_3 = IdentityMapping('M3', dim=2)(logical_patch_3) + patch_4 = IdentityMapping('M4', dim=2)(logical_patch_4) + patch_5 = IdentityMapping('M5', dim=2)(logical_patch_5) + patch_6 = IdentityMapping('M6', dim=2)(logical_patch_6) + patch_7 = IdentityMapping('M7', dim=2)(logical_patch_7) + patch_8 = IdentityMapping('M8', dim=2)(logical_patch_8) + + patches = ( + patch_1, patch_2, patch_3, patch_4, + patch_5, patch_6, patch_7, patch_8, + ) + connectivity = ( + ((patch_1, 0, +1), (patch_2, 0, -1), +1), + ((patch_2, 0, +1), (patch_3, 0, -1), +1), + ((patch_6, 0, +1), (patch_7, 0, -1), +1), + ((patch_7, 0, +1), (patch_8, 0, -1), +1), + ((patch_1, 1, +1), (patch_4, 1, -1), +1), + ((patch_4, 1, +1), (patch_6, 1, -1), +1), + ((patch_3, 1, +1), (patch_5, 1, -1), +1), + ((patch_5, 1, +1), (patch_8, 1, -1), +1), + ) + return Domain.join(patches, connectivity, name='square_8') + + +def build_square_9(): + """Build a square decomposed into a three-by-three patch grid.""" + logical_patch_1 = Square( + 'P1', + bounds1=(0., np.pi / 3), + bounds2=(0., np.pi / 3), + ) + logical_patch_2 = Square( + 'P2', + bounds1=(np.pi / 3, 2 * np.pi / 3), + bounds2=(0., np.pi / 3), + ) + logical_patch_3 = Square( + 'P3', + bounds1=(2 * np.pi / 3, np.pi), + bounds2=(0., np.pi / 3), + ) + logical_patch_4 = Square( + 'P4', + bounds1=(0., np.pi / 3), + bounds2=(np.pi / 3, 2 * np.pi / 3), + ) + logical_patch_5 = Square( + 'P5', + bounds1=(2 * np.pi / 3, np.pi), + bounds2=(np.pi / 3, 2 * np.pi / 3), + ) + logical_patch_6 = Square( + 'P6', + bounds1=(0., np.pi / 3), + bounds2=(2 * np.pi / 3, np.pi), + ) + logical_patch_7 = Square( + 'P7', + bounds1=(np.pi / 3, 2 * np.pi / 3), + bounds2=(2 * np.pi / 3, np.pi), + ) + logical_patch_8 = Square( + 'P8', + bounds1=(2 * np.pi / 3, np.pi), + bounds2=(2 * np.pi / 3, np.pi), + ) + logical_patch_9 = Square( + 'P9', + bounds1=(np.pi / 3, 2 * np.pi / 3), + bounds2=(np.pi / 3, 2 * np.pi / 3), + ) + patch_1 = IdentityMapping('M1', dim=2)(logical_patch_1) + patch_2 = IdentityMapping('M2', dim=2)(logical_patch_2) + patch_3 = IdentityMapping('M3', dim=2)(logical_patch_3) + patch_4 = IdentityMapping('M4', dim=2)(logical_patch_4) + patch_5 = IdentityMapping('M5', dim=2)(logical_patch_5) + patch_6 = IdentityMapping('M6', dim=2)(logical_patch_6) + patch_7 = IdentityMapping('M7', dim=2)(logical_patch_7) + patch_8 = IdentityMapping('M8', dim=2)(logical_patch_8) + patch_9 = IdentityMapping('M9', dim=2)(logical_patch_9) + + patches = ( + patch_1, patch_2, patch_3, patch_4, patch_5, + patch_6, patch_7, patch_8, patch_9, + ) + connectivity = ( + ((patch_1, 0, +1), (patch_2, 0, -1), +1), + ((patch_2, 0, +1), (patch_3, 0, -1), +1), + ((patch_4, 0, +1), (patch_9, 0, -1), +1), + ((patch_9, 0, +1), (patch_5, 0, -1), +1), + ((patch_6, 0, +1), (patch_7, 0, -1), +1), + ((patch_7, 0, +1), (patch_8, 0, -1), +1), + ((patch_1, 1, +1), (patch_4, 1, -1), +1), + ((patch_4, 1, +1), (patch_6, 1, -1), +1), + ((patch_2, 1, +1), (patch_9, 1, -1), +1), + ((patch_9, 1, +1), (patch_7, 1, -1), +1), + ((patch_3, 1, +1), (patch_5, 1, -1), +1), + ((patch_5, 1, +1), (patch_8, 1, -1), +1), + ) + return Domain.join(patches, connectivity, name='square_9') + + +def build_annulus_3(r_min=None, r_max=None): + """Build a complete annulus from three polar patches.""" + r_min = 0.5 if r_min is None else r_min + r_max = 1.0 if r_max is None else r_max + if not 0 < r_min < r_max: + raise ValueError('annulus radii must satisfy 0 < r_min < r_max') + + logical_patch_1 = Square( + 'P1', + bounds1=(r_min, r_max), + bounds2=(0., np.pi / 2), + ) + logical_patch_2 = Square( + 'P2', + bounds1=(r_min, r_max), + bounds2=(np.pi / 2, np.pi), + ) + logical_patch_3 = Square( + 'P3', + bounds1=(r_min, r_max), + bounds2=(np.pi, 2 * np.pi), + ) + mapping_1 = PolarMapping( + 'M1', dim=2, c1=0., c2=0., rmin=0., rmax=1.) + mapping_2 = PolarMapping( + 'M2', dim=2, c1=0., c2=0., rmin=0., rmax=1.) + mapping_3 = PolarMapping( + 'M3', dim=2, c1=0., c2=0., rmin=0., rmax=1.) + patch_1 = mapping_1(logical_patch_1) + patch_2 = mapping_2(logical_patch_2) + patch_3 = mapping_3(logical_patch_3) + + patches = (patch_1, patch_2, patch_3) + connectivity = ( + ((patch_1, 1, +1), (patch_2, 1, -1), +1), + ((patch_2, 1, +1), (patch_3, 1, -1), +1), + ((patch_3, 1, +1), (patch_1, 1, -1), +1), + ) + return Domain.join(patches, connectivity, name='annulus_3') + + +def build_annulus_4(r_min=None, r_max=None): + """Build a complete annulus from four quarter-annulus patches.""" + r_min = 0.5 if r_min is None else r_min + r_max = 1.0 if r_max is None else r_max + if not 0 < r_min < r_max: + raise ValueError('annulus radii must satisfy 0 < r_min < r_max') + + logical_patch_1 = Square( + 'P1', + bounds1=(r_min, r_max), + bounds2=(0., np.pi / 2), + ) + logical_patch_2 = Square( + 'P2', + bounds1=(r_min, r_max), + bounds2=(np.pi / 2, np.pi), + ) + logical_patch_3 = Square( + 'P3', + bounds1=(r_min, r_max), + bounds2=(np.pi, 3 * np.pi / 2), + ) + logical_patch_4 = Square( + 'P4', + bounds1=(r_min, r_max), + bounds2=(3 * np.pi / 2, 2 * np.pi), + ) + mapping_1 = PolarMapping( + 'M1', dim=2, c1=0., c2=0., rmin=0., rmax=1.) + mapping_2 = PolarMapping( + 'M2', dim=2, c1=0., c2=0., rmin=0., rmax=1.) + mapping_3 = PolarMapping( + 'M3', dim=2, c1=0., c2=0., rmin=0., rmax=1.) + mapping_4 = PolarMapping( + 'M4', dim=2, c1=0., c2=0., rmin=0., rmax=1.) + patch_1 = mapping_1(logical_patch_1) + patch_2 = mapping_2(logical_patch_2) + patch_3 = mapping_3(logical_patch_3) + patch_4 = mapping_4(logical_patch_4) + + patches = (patch_1, patch_2, patch_3, patch_4) + connectivity = ( + ((patch_1, 1, +1), (patch_2, 1, -1), +1), + ((patch_2, 1, +1), (patch_3, 1, -1), +1), + ((patch_3, 1, +1), (patch_4, 1, -1), +1), + ((patch_4, 1, +1), (patch_1, 1, -1), +1), + ) + return Domain.join(patches, connectivity, name='annulus_4') + + +def build_curved_l_shape(): + """Build the three-patch curved L-shaped benchmark domain.""" + logical_patch_1 = Square( + 'P1', + bounds1=(2, 3), + bounds2=(0., np.pi / 8), + ) + logical_patch_2 = Square( + 'P2', + bounds1=(2, 3), + bounds2=(np.pi / 8, np.pi / 4), + ) + logical_patch_3 = Square( + 'P3', + bounds1=(1, 2), + bounds2=(np.pi / 8, np.pi / 4), + ) + mapping_1 = PolarMapping( + 'M1', dim=2, c1=0., c2=0., rmin=0., rmax=1.) + mapping_2 = PolarMapping( + 'M2', dim=2, c1=0., c2=0., rmin=0., rmax=1.) + mapping_3 = PolarMapping( + 'M3', dim=2, c1=0., c2=0., rmin=0., rmax=1.) + patch_1 = mapping_1(logical_patch_1) + patch_2 = mapping_2(logical_patch_2) + patch_3 = mapping_3(logical_patch_3) + + patches = (patch_1, patch_2, patch_3) + connectivity = ( + ((patch_1, 1, +1), (patch_2, 1, -1), +1), + ((patch_3, 0, +1), (patch_2, 0, -1), +1), + ) + return Domain.join(patches, connectivity, name='curved_L_shape') + + +def build_pretzel(r_min=None, r_max=None): + """Build the coarse eleven-patch pretzel domain.""" + r_min = 1 if r_min is None else r_min + r_max = 2 if r_max is None else r_max + if not 0 < r_min < r_max: + raise ValueError('pretzel radii must satisfy 0 < r_min < r_max') + + h = r_max - r_min + half_h = h / 2 + center_radius = h + (r_max + r_min) / 2 + + logical_patch_1 = Square( + 'P1', + bounds1=(r_min, r_max), + bounds2=(0., np.pi / 2), + ) + logical_patch_2 = Square( + 'P2', + bounds1=(r_min, r_max), + bounds2=(np.pi / 2, np.pi), + ) + logical_patch_3 = Square( + 'P3', + bounds1=(r_min, r_max), + bounds2=(np.pi, 3 * np.pi / 2), + ) + logical_patch_4 = Square( + 'P4', + bounds1=(r_min, r_max), + bounds2=(3 * np.pi / 2, 2 * np.pi), + ) + mapping_1 = PolarMapping( + 'M1', dim=2, c1=h, c2=h, rmin=0., rmax=1.) + mapping_2 = PolarMapping( + 'M2', dim=2, c1=-h, c2=h, rmin=0., rmax=1.) + mapping_3 = PolarMapping( + 'M3', dim=2, c1=-h, c2=0, rmin=0., rmax=1.) + mapping_4 = PolarMapping( + 'M4', dim=2, c1=h, c2=0, rmin=0., rmax=1.) + patch_1 = mapping_1(logical_patch_1) + patch_2 = mapping_2(logical_patch_2) + patch_3 = mapping_3(logical_patch_3) + patch_4 = mapping_4(logical_patch_4) + + logical_patch_5 = Square( + 'P5', + bounds1=(-half_h, half_h), + bounds2=(-h / 2, h / 2), + ) + logical_patch_6 = Square( + 'P6', + bounds1=(-half_h, half_h), + bounds2=(-h / 2, h / 2), + ) + logical_patch_7 = Square( + 'P7', + bounds1=(-half_h, half_h), + bounds2=(-h / 2, h / 2), + ) + logical_patch_9 = Square( + 'P9', + bounds1=(-half_h, half_h), + bounds2=(-h, h), + ) + logical_patch_12 = Square( + 'P12', + bounds1=(-half_h, half_h), + bounds2=(-h / 2, h / 2), + ) + mapping_5 = AffineMapping( + 'M5', dim=2, + c1=h / 2, c2=center_radius, + a11=np.cos(np.pi / 2), a12=-np.sin(np.pi / 2), + a21=np.sin(np.pi / 2), a22=np.cos(np.pi / 2), + ) + mapping_6 = AffineMapping( + 'M6', dim=2, + c1=-h / 2, c2=center_radius, + a11=0, a12=1, a21=1, a22=0, + ) + mapping_7 = AffineMapping( + 'M7', dim=2, + c1=-center_radius, c2=h / 2, + a11=np.cos(np.pi), a12=-np.sin(np.pi), + a21=np.sin(np.pi), a22=np.cos(np.pi), + ) + mapping_9 = AffineMapping( + 'M9', dim=2, + c1=0, c2=h - center_radius, + a11=np.cos(3 * np.pi / 2), a12=-np.sin(3 * np.pi / 2), + a21=np.sin(3 * np.pi / 2), a22=np.cos(3 * np.pi / 2), + ) + mapping_12 = AffineMapping( + 'M12', dim=2, + c1=center_radius, c2=h / 2, + a11=1, a12=0, a21=0, a22=-1, + ) + patch_5 = mapping_5(logical_patch_5) + patch_6 = mapping_6(logical_patch_6) + patch_7 = mapping_7(logical_patch_7) + patch_9 = mapping_9(logical_patch_9) + patch_12 = mapping_12(logical_patch_12) + + logical_patch_13 = Square( + 'P13', + bounds1=(3 * np.pi / 2, 2 * np.pi), + bounds2=(r_min, r_max), + ) + logical_patch_14 = Square( + 'P14', + bounds1=(np.pi, 3 * np.pi / 2), + bounds2=(r_min, r_max), + ) + mapping_13 = TransposedPolarMapping( + 'M13', dim=2, + c1=-r_min - h, c2=r_min + h, rmin=0., rmax=1., + ) + mapping_14 = TransposedPolarMapping( + 'M14', dim=2, + c1=r_min + h, c2=r_min + h, rmin=0., rmax=1., + ) + patch_13 = mapping_13(logical_patch_13) + patch_14 = mapping_14(logical_patch_14) + + patches = ( + patch_1, patch_2, patch_3, patch_4, patch_5, patch_6, + patch_7, patch_9, patch_12, patch_13, patch_14, + ) + connectivity = ( + ((patch_1, 1, +1), (patch_5, 1, -1), +1), + ((patch_5, 1, +1), (patch_6, 1, +1), +1), + ((patch_6, 1, -1), (patch_2, 1, -1), +1), + ((patch_2, 1, +1), (patch_7, 1, -1), +1), + ((patch_7, 1, +1), (patch_3, 1, -1), +1), + ((patch_3, 1, +1), (patch_9, 1, -1), +1), + ((patch_9, 1, +1), (patch_4, 1, -1), +1), + ((patch_4, 1, +1), (patch_12, 1, +1), +1), + ((patch_12, 1, -1), (patch_1, 1, -1), +1), + ((patch_6, 0, -1), (patch_13, 0, +1), +1), + ((patch_7, 0, -1), (patch_13, 0, -1), +1), + ((patch_5, 0, -1), (patch_14, 0, -1), +1), + ((patch_12, 0, -1), (patch_14, 0, +1), +1), + ) + return Domain.join(patches, connectivity, name='pretzel') + + +def build_pretzel_annulus(r_min=None, r_max=None): + """Build the nine-patch outer annular part of the pretzel domain.""" + r_min = 1 if r_min is None else r_min + r_max = 2 if r_max is None else r_max + if not 0 < r_min < r_max: + raise ValueError('pretzel radii must satisfy 0 < r_min < r_max') + + h = r_max - r_min + half_h = h / 2 + center_radius = h + (r_max + r_min) / 2 + + logical_patch_1 = Square( + 'P1', bounds1=(r_min, r_max), bounds2=(0., np.pi / 2)) + logical_patch_2 = Square( + 'P2', bounds1=(r_min, r_max), bounds2=(np.pi / 2, np.pi)) + logical_patch_3 = Square( + 'P3', bounds1=(r_min, r_max), + bounds2=(np.pi, 3 * np.pi / 2)) + logical_patch_4 = Square( + 'P4', bounds1=(r_min, r_max), + bounds2=(3 * np.pi / 2, 2 * np.pi)) + mapping_1 = PolarMapping( + 'M1', dim=2, c1=h, c2=h, rmin=0., rmax=1.) + mapping_2 = PolarMapping( + 'M2', dim=2, c1=-h, c2=h, rmin=0., rmax=1.) + mapping_3 = PolarMapping( + 'M3', dim=2, c1=-h, c2=0, rmin=0., rmax=1.) + mapping_4 = PolarMapping( + 'M4', dim=2, c1=h, c2=0, rmin=0., rmax=1.) + patch_1 = mapping_1(logical_patch_1) + patch_2 = mapping_2(logical_patch_2) + patch_3 = mapping_3(logical_patch_3) + patch_4 = mapping_4(logical_patch_4) + + logical_patch_5 = Square( + 'P5', bounds1=(-half_h, half_h), bounds2=(-h / 2, h / 2)) + logical_patch_6 = Square( + 'P6', bounds1=(-half_h, half_h), bounds2=(-h / 2, h / 2)) + logical_patch_7 = Square( + 'P7', bounds1=(-half_h, half_h), bounds2=(-h / 2, h / 2)) + logical_patch_9 = Square( + 'P9', bounds1=(-half_h, half_h), bounds2=(-h, h)) + logical_patch_12 = Square( + 'P12', bounds1=(-half_h, half_h), bounds2=(-h / 2, h / 2)) + mapping_5 = AffineMapping( + 'M5', dim=2, + c1=h / 2, c2=center_radius, + a11=np.cos(np.pi / 2), a12=-np.sin(np.pi / 2), + a21=np.sin(np.pi / 2), a22=np.cos(np.pi / 2), + ) + mapping_6 = AffineMapping( + 'M6', dim=2, + c1=-h / 2, c2=center_radius, + a11=0, a12=1, a21=1, a22=0, + ) + mapping_7 = AffineMapping( + 'M7', dim=2, + c1=-center_radius, c2=h / 2, + a11=np.cos(np.pi), a12=-np.sin(np.pi), + a21=np.sin(np.pi), a22=np.cos(np.pi), + ) + mapping_9 = AffineMapping( + 'M9', dim=2, + c1=0, c2=h - center_radius, + a11=np.cos(3 * np.pi / 2), a12=-np.sin(3 * np.pi / 2), + a21=np.sin(3 * np.pi / 2), a22=np.cos(3 * np.pi / 2), + ) + mapping_12 = AffineMapping( + 'M12', dim=2, + c1=center_radius, c2=h / 2, + a11=1, a12=0, a21=0, a22=-1, + ) + patch_5 = mapping_5(logical_patch_5) + patch_6 = mapping_6(logical_patch_6) + patch_7 = mapping_7(logical_patch_7) + patch_9 = mapping_9(logical_patch_9) + patch_12 = mapping_12(logical_patch_12) + + patches = ( + patch_1, patch_5, patch_6, patch_2, patch_7, + patch_3, patch_9, patch_4, patch_12, + ) + connectivity = ( + ((patch_1, 1, +1), (patch_5, 1, -1), +1), + ((patch_5, 1, +1), (patch_6, 1, +1), +1), + ((patch_6, 1, -1), (patch_2, 1, -1), +1), + ((patch_2, 1, +1), (patch_7, 1, -1), +1), + ((patch_7, 1, +1), (patch_3, 1, -1), +1), + ((patch_3, 1, +1), (patch_9, 1, -1), +1), + ((patch_9, 1, +1), (patch_4, 1, -1), +1), + ((patch_4, 1, +1), (patch_12, 1, +1), +1), + ((patch_12, 1, -1), (patch_1, 1, -1), +1), + ) + return Domain.join(patches, connectivity, name='pretzel_annulus') + + +def build_pretzel_debug(r_min=None, r_max=None): + """Build the two-patch pretzel geometry used for interface debugging.""" + r_min = 1 if r_min is None else r_min + r_max = 2 if r_max is None else r_max + if not 0 < r_min < r_max: + raise ValueError('pretzel radii must satisfy 0 < r_min < r_max') + + h = r_max - r_min + logical_patch_1 = Square( + 'P1', + bounds1=(r_min, r_max), + bounds2=(0., np.pi / 2), + ) + logical_patch_10 = Square( + 'P10', + bounds1=(r_min, r_max), + bounds2=(np.pi / 2, np.pi), + ) + mapping_1 = PolarMapping( + 'M1', dim=2, c1=h, c2=h, rmin=0., rmax=1.) + mapping_10 = PolarMapping( + 'M10', dim=2, c1=h, c2=h, rmin=0., rmax=1.) + patch_1 = mapping_1(logical_patch_1) + patch_10 = mapping_10(logical_patch_10) + + patches = (patch_1, patch_10) + connectivity = ( + ((patch_1, 1, +1), (patch_10, 1, -1), +1), + ) + return Domain.join(patches, connectivity, name='pretzel_debug') + + +def build_pretzel_f(r_min=None, r_max=None): + """Build the refined eighteen-patch pretzel domain.""" + r_min = 1 if r_min is None else r_min + r_max = 2 if r_max is None else r_max + if not 0 < r_min < r_max: + raise ValueError('pretzel radii must satisfy 0 < r_min < r_max') + + h = r_max - r_min + half_h = h / 2 + center_radius = h + (r_max + r_min) / 2 + + logical_patch_1_1 = Square( + 'P1_1', + bounds1=(r_min, r_max), + bounds2=(0., np.pi / 4), + ) + logical_patch_1_2 = Square( + 'P1_2', + bounds1=(r_min, r_max), + bounds2=(np.pi / 4, np.pi / 2), + ) + logical_patch_2_1 = Square( + 'P2_1', + bounds1=(r_min, r_max), + bounds2=(np.pi / 2, 3 * np.pi / 4), + ) + logical_patch_2_2 = Square( + 'P2_2', + bounds1=(r_min, r_max), + bounds2=(3 * np.pi / 4, np.pi), + ) + logical_patch_3_1 = Square( + 'P3_1', + bounds1=(r_min, r_max), + bounds2=(np.pi, 5 * np.pi / 4), + ) + logical_patch_3_2 = Square( + 'P3_2', + bounds1=(r_min, r_max), + bounds2=(5 * np.pi / 4, 3 * np.pi / 2), + ) + logical_patch_4_1 = Square( + 'P4_1', + bounds1=(r_min, r_max), + bounds2=(3 * np.pi / 2, 7 * np.pi / 4), + ) + logical_patch_4_2 = Square( + 'P4_2', + bounds1=(r_min, r_max), + bounds2=(7 * np.pi / 4, 2 * np.pi), + ) + mapping_1_1 = PolarMapping( + 'M1_1', dim=2, c1=h, c2=h, rmin=0., rmax=1.) + mapping_1_2 = PolarMapping( + 'M1_2', dim=2, c1=h, c2=h, rmin=0., rmax=1.) + mapping_2_1 = PolarMapping( + 'M2_1', dim=2, c1=-h, c2=h, rmin=0., rmax=1.) + mapping_2_2 = PolarMapping( + 'M2_2', dim=2, c1=-h, c2=h, rmin=0., rmax=1.) + mapping_3_1 = PolarMapping( + 'M3_1', dim=2, c1=-h, c2=0, rmin=0., rmax=1.) + mapping_3_2 = PolarMapping( + 'M3_2', dim=2, c1=-h, c2=0, rmin=0., rmax=1.) + mapping_4_1 = PolarMapping( + 'M4_1', dim=2, c1=h, c2=0, rmin=0., rmax=1.) + mapping_4_2 = PolarMapping( + 'M4_2', dim=2, c1=h, c2=0, rmin=0., rmax=1.) + patch_1_1 = mapping_1_1(logical_patch_1_1) + patch_1_2 = mapping_1_2(logical_patch_1_2) + patch_2_1 = mapping_2_1(logical_patch_2_1) + patch_2_2 = mapping_2_2(logical_patch_2_2) + patch_3_1 = mapping_3_1(logical_patch_3_1) + patch_3_2 = mapping_3_2(logical_patch_3_2) + patch_4_1 = mapping_4_1(logical_patch_4_1) + patch_4_2 = mapping_4_2(logical_patch_4_2) + + logical_patch_5 = Square( + 'P5', + bounds1=(-half_h, half_h), + bounds2=(-h / 2, h / 2), + ) + logical_patch_6 = Square( + 'P6', + bounds1=(-half_h, half_h), + bounds2=(-h / 2, h / 2), + ) + logical_patch_7 = Square( + 'P7', + bounds1=(-half_h, half_h), + bounds2=(-h / 2, h / 2), + ) + logical_patch_9_1 = Square( + 'P9_1', + bounds1=(-half_h, half_h), + bounds2=(-h, 0), + ) + logical_patch_9_2 = Square( + 'P9_2', + bounds1=(-half_h, half_h), + bounds2=(0, h), + ) + logical_patch_12 = Square( + 'P12', + bounds1=(-half_h, half_h), + bounds2=(-h / 2, h / 2), + ) + mapping_5 = AffineMapping( + 'M5', dim=2, + c1=h / 2, c2=center_radius, + a11=np.cos(np.pi / 2), a12=-np.sin(np.pi / 2), + a21=np.sin(np.pi / 2), a22=np.cos(np.pi / 2), + ) + mapping_6 = AffineMapping( + 'M6', dim=2, + c1=-h / 2, c2=center_radius, + a11=0, a12=1, a21=1, a22=0, + ) + mapping_7 = AffineMapping( + 'M7', dim=2, + c1=-center_radius, c2=h / 2, + a11=np.cos(np.pi), a12=-np.sin(np.pi), + a21=np.sin(np.pi), a22=np.cos(np.pi), + ) + mapping_9_1 = AffineMapping( + 'M9_1', dim=2, + c1=0, c2=h - center_radius, + a11=np.cos(3 * np.pi / 2), a12=-np.sin(3 * np.pi / 2), + a21=np.sin(3 * np.pi / 2), a22=np.cos(3 * np.pi / 2), + ) + mapping_9_2 = AffineMapping( + 'M9_2', dim=2, + c1=0, c2=h - center_radius, + a11=np.cos(3 * np.pi / 2), a12=-np.sin(3 * np.pi / 2), + a21=np.sin(3 * np.pi / 2), a22=np.cos(3 * np.pi / 2), + ) + mapping_12 = AffineMapping( + 'M12', dim=2, + c1=center_radius, c2=h / 2, + a11=1, a12=0, a21=0, a22=-1, + ) + patch_5 = mapping_5(logical_patch_5) + patch_6 = mapping_6(logical_patch_6) + patch_7 = mapping_7(logical_patch_7) + patch_9_1 = mapping_9_1(logical_patch_9_1) + patch_9_2 = mapping_9_2(logical_patch_9_2) + patch_12 = mapping_12(logical_patch_12) + + logical_patch_13_1 = Square( + 'P13_1', + bounds1=(3 * np.pi / 2, 7 * np.pi / 4), + bounds2=(r_min, r_max), + ) + logical_patch_13_2 = Square( + 'P13_2', + bounds1=(7 * np.pi / 4, 2 * np.pi), + bounds2=(r_min, r_max), + ) + logical_patch_14_1 = Square( + 'P14_1', + bounds1=(np.pi, 5 * np.pi / 4), + bounds2=(r_min, r_max), + ) + logical_patch_14_2 = Square( + 'P14_2', + bounds1=(5 * np.pi / 4, 3 * np.pi / 2), + bounds2=(r_min, r_max), + ) + mapping_13_1 = TransposedPolarMapping( + 'M13_1', dim=2, + c1=-r_min - h, c2=r_min + h, rmin=0., rmax=1., + ) + mapping_13_2 = TransposedPolarMapping( + 'M13_2', dim=2, + c1=-r_min - h, c2=r_min + h, rmin=0., rmax=1., + ) + mapping_14_1 = TransposedPolarMapping( + 'M14_1', dim=2, + c1=r_min + h, c2=r_min + h, rmin=0., rmax=1., + ) + mapping_14_2 = TransposedPolarMapping( + 'M14_2', dim=2, + c1=r_min + h, c2=r_min + h, rmin=0., rmax=1., + ) + patch_13_1 = mapping_13_1(logical_patch_13_1) + patch_13_2 = mapping_13_2(logical_patch_13_2) + patch_14_1 = mapping_14_1(logical_patch_14_1) + patch_14_2 = mapping_14_2(logical_patch_14_2) + + patches = ( + patch_1_1, patch_1_2, patch_2_1, patch_2_2, + patch_3_1, patch_3_2, patch_4_1, patch_4_2, + patch_5, patch_6, patch_7, patch_9_1, patch_9_2, patch_12, + patch_13_1, patch_13_2, patch_14_1, patch_14_2, + ) + connectivity = ( + ((patch_1_1, 1, +1), (patch_1_2, 1, -1), +1), + ((patch_1_2, 1, +1), (patch_5, 1, -1), +1), + ((patch_5, 1, +1), (patch_6, 1, +1), +1), + ((patch_6, 1, -1), (patch_2_1, 1, -1), +1), + ((patch_2_1, 1, +1), (patch_2_2, 1, -1), +1), + ((patch_2_2, 1, +1), (patch_7, 1, -1), +1), + ((patch_7, 1, +1), (patch_3_1, 1, -1), +1), + ((patch_3_1, 1, +1), (patch_3_2, 1, -1), +1), + ((patch_3_2, 1, +1), (patch_9_1, 1, -1), +1), + ((patch_9_1, 1, +1), (patch_9_2, 1, -1), +1), + ((patch_9_2, 1, +1), (patch_4_1, 1, -1), +1), + ((patch_4_1, 1, +1), (patch_4_2, 1, -1), +1), + ((patch_4_2, 1, +1), (patch_12, 1, +1), +1), + ((patch_12, 1, -1), (patch_1_1, 1, -1), +1), + ((patch_6, 0, -1), (patch_13_2, 0, +1), +1), + ((patch_13_2, 0, -1), (patch_13_1, 0, +1), +1), + ((patch_7, 0, -1), (patch_13_1, 0, -1), +1), + ((patch_5, 0, -1), (patch_14_1, 0, -1), +1), + ((patch_14_1, 0, +1), (patch_14_2, 0, -1), +1), + ((patch_12, 0, -1), (patch_14_2, 0, +1), +1), + ) + return Domain.join(patches, connectivity, name='pretzel_f') + + +def build_two_patch_twisted_strip(): + """Return two patches connected by interfaces of both 2D orientations.""" + logical_patch_a = Square('A', bounds1=(0, 1), bounds2=(0, 1)) + logical_patch_b = Square('B', bounds1=(0, 1), bounds2=(0, 1)) + + patch_a = IdentityMapping('F_A', dim=2)(logical_patch_a) + patch_b = AffineMapping( + 'F_B', dim=2, + c1=1, c2=0, + a11=1, a12=0, + a21=0, a22=1, + )(logical_patch_b) + + return Domain.join( + patches=[patch_a, patch_b], + interfaces=[ + # Touching sides preserve their tangential coordinates. + ((patch_a, 0, +1), (patch_b, 0, -1), +1), + # Outer sides are identified with reversed coordinates. + ((patch_a, 0, -1), (patch_b, 0, +1), -1), + ], + name='two_patch_twisted_strip', + ) + + +def build_self_connected_patch(): + """Return a single square with a left-to-right self-interface.""" + logical_patch = Square('P', bounds1=(0, 1), bounds2=(0, 1)) + patch = IdentityMapping('F_P', dim=2)(logical_patch) + + return Domain.join( + patches=[patch], + interfaces=[((patch, 0, -1), (patch, 0, +1), +1)], + name='self_connected_patch', + ) + + +def build_three_patch_shared_vertex(): + """Return three affine patches meeting at one interior vertex.""" + root_three_over_two = np.sqrt(3.0) / 2.0 + directions = ( + np.array((1.0, 0.0)), + np.array((-0.5, root_three_over_two)), + np.array((-0.5, -root_three_over_two)), + ) + patches = [] + for patch_name, first, second in ( + ('A', directions[0], directions[1]), + ('B', directions[1], directions[2]), + ('C', directions[2], directions[0]), + ): + logical_patch = Square( + patch_name, bounds1=(0, 1), bounds2=(0, 1)) + mapping = AffineMapping( + f'F_{patch_name}', dim=2, + c1=0, c2=0, + a11=first[0], a12=second[0], + a21=first[1], a22=second[1], + ) + patches.append(mapping(logical_patch)) + + patch_a, patch_b, patch_c = patches + return Domain.join( + patches=patches, + interfaces=[ + ((patch_a, 0, -1), (patch_b, 1, -1), +1), + ((patch_b, 0, -1), (patch_c, 1, -1), +1), + ((patch_c, 0, -1), (patch_a, 1, -1), +1), + ], + name='three_patch_shared_vertex', + ) + + +_DOMAIN_BUILDERS_2D = { + 'square_2': build_square_2, + 'square_4': build_square_4, + 'square_6': build_square_6, + 'square_8': build_square_8, + 'square_9': build_square_9, + 'annulus_3': build_annulus_3, + 'annulus_4': build_annulus_4, + 'curved_L_shape': build_curved_l_shape, + 'pretzel': build_pretzel, + 'pretzel_f': build_pretzel_f, + 'pretzel_annulus': build_pretzel_annulus, + 'pretzel_debug': build_pretzel_debug, + 'two_patch_twisted_strip': build_two_patch_twisted_strip, + 'self_connected_patch': build_self_connected_patch, + 'three_patch_shared_vertex': build_three_patch_shared_vertex, +} + +_RADIAL_DOMAINS_2D = { + 'annulus_3', + 'annulus_4', + 'pretzel', + 'pretzel_f', + 'pretzel_annulus', + 'pretzel_debug', +} + + +def build_multipatch_domain_2d(domain_name='square_2', r_min=None, r_max=None): + """Build a named 2D multipatch domain from the gallery. + + Parameters + ---------- + domain_name : str + Name registered in the 2D multipatch gallery. + r_min, r_max : float, optional + Inner and outer radii for annulus and pretzel domains. + """ + try: + builder = _DOMAIN_BUILDERS_2D[domain_name] + except (KeyError, TypeError) as error: + choices = ', '.join(_DOMAIN_BUILDERS_2D) + raise ValueError( + f'unknown 2D multipatch domain {domain_name!r}; ' + f'choose from {choices}') from error + + if domain_name in _RADIAL_DOMAINS_2D: + return builder(r_min=r_min, r_max=r_max) + return builder() + + +def build_cartesian_multipatch_domain_2d( + ncells, log_interval_x, log_interval_y, mapping='identity'): + """Create a 2D multipatch domain from a rectangular patch layout. + + Parameters + ---------- + ncells : array-like + Two-dimensional patch layout. Non-``None`` entries create patches; + their numerical values are ignored by this symbolic builder and may + subsequently be used as per-patch cell counts. Row zero is the top + row of the domain. + log_interval_x : tuple + The interval in the x direction in the logical domain. + log_interval_y : tuple + The interval in the y direction in the logical domain. + mapping : str + The type of mapping to use. Can be ``identity`` or ``polar``. + + Returns + ------- + domain : Domain + The symbolic multipatch domain. + """ + layout = np.asarray(ncells, dtype=object) + if layout.ndim != 2: + raise ValueError('ncells must be a two-dimensional patch layout') + + nrows, ncols = layout.shape + if nrows == 0 or ncols == 0: + raise ValueError('ncells must contain at least one layout position') + + try: + ax, bx = log_interval_x + ay, by = log_interval_y + except (TypeError, ValueError) as error: + raise ValueError( + 'logical intervals must each contain exactly two bounds') from error + + mapping_types = { + 'identity': IdentityMapping, + 'polar': PolarMapping, + } + try: + mapping_type = mapping_types[mapping] + except (KeyError, TypeError) as error: + raise ValueError( + "mapping must be either 'identity' or 'polar'") from error + + patches_by_position = {} + patches = [] + for row in range(nrows): + bounds2 = ( + by - (row + 1) / nrows * (by - ay), + by - row / nrows * (by - ay), + ) + for column in range(ncols): + if layout[row, column] is None: + continue + + bounds1 = ( + ax + column / ncols * (bx - ax), + ax + (column + 1) / ncols * (bx - ax), + ) + logical_patch = Square( + f'Log_{row}_{column}', + bounds1=bounds1, + bounds2=bounds2, + ) + mapping_kwargs = {} + if mapping == 'polar': + mapping_kwargs = dict(c1=0., c2=0., rmin=0., rmax=1.) + patch_mapping = mapping_type( + f'M_{row}_{column}', dim=2, **mapping_kwargs) + patch = patch_mapping(logical_patch) + patches_by_position[row, column] = patch + patches.append(patch) + + if not patches: + raise ValueError('ncells must contain at least one non-None entry') + + connectivity = [] + for row in range(nrows): + for column in range(ncols - 1): + left = patches_by_position.get((row, column)) + right = patches_by_position.get((row, column + 1)) + if left is not None and right is not None: + connectivity.append( + ((left, 0, +1), (right, 0, -1), +1)) + + for row in range(nrows - 1): + for column in range(ncols): + top = patches_by_position.get((row, column)) + bottom = patches_by_position.get((row + 1, column)) + if top is not None and bottom is not None: + connectivity.append( + ((top, 1, -1), (bottom, 1, +1), +1)) + + return Domain.join(patches, connectivity, name='domain') + + +# ============================================================================= +# 3D domains +# ============================================================================= + +def build_two_cube_3d(): + """Build two mapped cubes joined across differently parameterized faces.""" + logical_patch_a = Cube('A') + logical_patch_b = Cube('B') + + mapping_a = AffineMapping( + 'F_A', dim=3, + c1=0.0, c2=0.0, c3=0.0, + a11=1.0, a12=0.25, a13=0.0, + a21=0.2, a22=1.00, a23=0.15, + a31=0.0, a32=0.10, a33=1.0, + ) + mapping_b = AffineMapping( + 'F_B', dim=3, + c1=1.0, c2=0.35, c3=1.0, + a11=0.00, a12=1.0, a13=0.25, + a21=-0.15, a22=0.2, a23=1.00, + a31=-1.00, a32=0.0, a33=0.10, + ) + patch_a = mapping_a(logical_patch_a) + patch_b = mapping_b(logical_patch_b) + + patches = (patch_a, patch_b) + connectivity = ( + ( + (patch_a, 0, +1), + (patch_b, 1, -1), + (-1, +1, -1), + ), + ) + return Domain.join( + patches, + connectivity, + name='mapped_two_cube_3d', + ) + + +def build_torus_3d( + major_radius=2.0, minor_bounds=(0.35, 0.80), *, hollow=True, + toroidal_angle=2.0 * np.pi, close_torus=None): + """Build a four-patch torus with two angular cuts per direction. + + Parameters + ---------- + major_radius : float, default=2.0 + Distance from the centre of the tube to the axis of revolution. + minor_bounds : tuple[float, float], default=(0.35, 0.80) + Inner and outer minor radii used by the regular hollow variant. The + lower value must be positive and the upper value must be smaller than + ``major_radius``. + hollow : bool, default=True + Keep the positive inner minor radius. If false, replace it by zero to + fill the tube. The resulting mapping is singular on the centreline: + the faces at minor radius zero collapse onto it, but they remain + exterior sides of the domain and part of ``domain.boundary``. + Exclude them when imposing boundary conditions on the solid torus. + toroidal_angle : float, default=2*pi + Total angle swept around the axis of revolution, in radians. It must + be greater than zero and at most ``2*pi``. The two toroidal patch + sectors each cover half of this angle. + close_torus : bool or None, default=None + Join the two toroidal end faces. ``None`` closes a full ``2*pi`` sweep + and leaves every smaller sweep open. Explicit closure is valid only + for a full sweep; ``False`` may also leave a full sweep cut open. + """ + inner_radius, outer_radius = map(float, minor_bounds) + major_radius = float(major_radius) + if not 0.0 < inner_radius < outer_radius < major_radius: + raise ValueError( + 'expected 0 < inner minor radius < outer minor radius ' + '< major radius') + if not isinstance(hollow, bool): + raise TypeError('hollow must be a bool') + if not hollow: + inner_radius = 0.0 + + full_angle = 2.0 * np.pi + toroidal_angle = float(toroidal_angle) + if not 0.0 < toroidal_angle <= full_angle or not np.isfinite( + toroidal_angle): + raise ValueError('toroidal_angle must satisfy 0 < angle <= 2*pi') + is_full_sweep = bool(np.isclose(toroidal_angle, full_angle)) + if is_full_sweep: + toroidal_angle = full_angle + if close_torus is None: + close_torus = is_full_sweep + elif not isinstance(close_torus, bool): + raise TypeError('close_torus must be a bool or None') + if close_torus and not is_full_sweep: + raise ValueError('only a full 2*pi sweep can be closed') + + theta_bounds = (0.0, np.pi, full_angle) + phi_bounds = (0.0, 0.5 * toroidal_angle, toroidal_angle) + patch_grid = {} + patches = [] + + for theta_sector in range(2): + for phi_sector in range(2): + patch_name = f'P{theta_sector}{phi_sector}' + logical_patch = Cube( + patch_name, + bounds1=(inner_radius, outer_radius), + bounds2=theta_bounds[theta_sector:theta_sector + 2], + bounds3=phi_bounds[phi_sector:phi_sector + 2], + ) + mapping = TorusMapping( + f'F_{patch_name}', dim=3, R0=major_radius) + patch = mapping(logical_patch) + patch_grid[theta_sector, phi_sector] = patch + patches.append(patch) + + identity = (+1, +1, +1) + connectivity = [] + for phi_sector in range(2): + lower = patch_grid[0, phi_sector] + upper = patch_grid[1, phi_sector] + connectivity.extend(( + ((lower, 1, +1), (upper, 1, -1), identity), + ((upper, 1, +1), (lower, 1, -1), identity), + )) + + for theta_sector in range(2): + lower = patch_grid[theta_sector, 0] + upper = patch_grid[theta_sector, 1] + connectivity.append( + ((lower, 2, +1), (upper, 2, -1), identity)) + if close_torus: + connectivity.append( + ((upper, 2, +1), (lower, 2, -1), identity)) + + shape = 'hollow' if hollow else 'solid' + closure = '' if close_torus else 'open_' + return Domain.join( + patches, + connectivity, + name=f'{closure}{shape}_torus_3d', + ) + + +def build_oriented_two_cube(): + """Return two cubes with a swapped and reversed face parameterization.""" + logical_patch_a = Cube('A') + logical_patch_b = Cube('B') + + mapping_a = IdentityMapping('F_A', dim=3) + mapping_b = AffineMapping( + 'F_B', dim=3, + c1=1, c2=0, c3=1, + a11=0, a12=1, a13=0, + a21=0, a22=0, a23=1, + a31=-1, a32=0, a33=0, + ) + patch_a = mapping_a(logical_patch_a) + patch_b = mapping_b(logical_patch_b) + + return Domain.join( + patches=[patch_a, patch_b], + interfaces=[( + (patch_a, 0, +1), + (patch_b, 1, -1), + (-1, +1, -1), + )], + name='oriented_two_cube', + ) + + +def build_three_patch_shared_edge(): + """Return three affine hexahedra meeting at one interior edge.""" + root_three_over_two = np.sqrt(3.0) / 2.0 + directions = ( + np.array((1.0, 0.0)), + np.array((-0.5, root_three_over_two)), + np.array((-0.5, -root_three_over_two)), + ) + + patches = [] + for patch_name, first, second in ( + ('H0', directions[0], directions[1]), + ('H1', directions[1], directions[2]), + ('H2', directions[2], directions[0]), + ): + logical_patch = Cube(patch_name) + mapping = AffineMapping( + f'G_{patch_name}', dim=3, + c1=0, c2=0, c3=0, + a11=first[0], a12=second[0], a13=0, + a21=first[1], a22=second[1], a23=0, + a31=0, a32=0, a33=1, + ) + patches.append(mapping(logical_patch)) + + patch_0, patch_1, patch_2 = patches + return Domain.join( + patches=patches, + interfaces=[ + ((patch_0, 0, -1), (patch_1, 1, -1), (+1, +1, +1)), + ((patch_1, 0, -1), (patch_2, 1, -1), (+1, +1, +1)), + ((patch_2, 0, -1), (patch_0, 1, -1), (+1, +1, +1)), + ], + name='three_patch_shared_edge', + ) + + +_DOMAIN_BUILDERS_3D = { + 'two_cube': build_two_cube_3d, + 'torus': build_torus_3d, + 'oriented_two_cube': build_oriented_two_cube, + 'three_patch_shared_edge': build_three_patch_shared_edge, +} + + +def build_multipatch_domain_3d(domain_name='two_cube', **kwargs): + """Build a named 3D multipatch domain from the gallery. + + Parameters + ---------- + domain_name : str + Name registered in the 3D multipatch gallery. + kwargs : object + Optional arguments passed to the selected domain builder. + """ + try: + builder = _DOMAIN_BUILDERS_3D[domain_name] + except (KeyError, TypeError) as error: + choices = ', '.join(_DOMAIN_BUILDERS_3D) + raise ValueError( + f'unknown 3D multipatch domain {domain_name!r}; ' + f'choose from {choices}') from error + + return builder(**kwargs) + + +# ============================================================================= +# Plotting +# ============================================================================= + +def plot_multipatch_domain( + domain_name='torus', *, output=None, show=True, topology=False, + builder_options=None): + """Build, print, and plot a named 2D or 3D gallery domain. + + Parameters + ---------- + domain_name : str + Name registered in the 2D or 3D multipatch gallery. + output : path-like or None, default=None + Optional path of the generated PNG or PDF figure. + show : bool, default=True + Open an interactive Matplotlib window. + topology : bool, default=False + Add patch, interface, orientation, and vertex annotations. + builder_options : dict or None, default=None + Optional arguments passed to the selected domain builder. + + Returns + ------- + domain : Domain + The selected symbolic multipatch domain. + figure : matplotlib.figure.Figure + The generated Matplotlib figure. + """ + options = {} if builder_options is None else dict(builder_options) + if domain_name in _DOMAIN_BUILDERS_2D: + domain = build_multipatch_domain_2d(domain_name, **options) + elif domain_name in _DOMAIN_BUILDERS_3D: + domain = build_multipatch_domain_3d(domain_name, **options) + else: + choices = ', '.join((*_DOMAIN_BUILDERS_2D, *_DOMAIN_BUILDERS_3D)) + raise ValueError( + f'unknown multipatch domain {domain_name!r}; ' + f'choose from {choices}') + + import matplotlib.pyplot as plt + from sympde.utilities.plotting import plot_domain, print_topology + + print_topology(domain) + figure = plot_domain( + domain, + draw=False, + isolines=not topology, + topology=topology, + interface_labels=topology, + vertex_labels=topology, + ) + if output is not None: + output = Path(output) + output.parent.mkdir(parents=True, exist_ok=True) + figure.savefig(output, dpi=180, bbox_inches='tight') + print(f'Saved plot to {output}') + + if show: + plt.show() + else: + plt.close(figure) + + return domain, figure + + +def main(argv=None): + """Plot a gallery domain selected from the command line.""" + domain_names = (*_DOMAIN_BUILDERS_2D, *_DOMAIN_BUILDERS_3D) + parser = ArgumentParser(description=__doc__) + parser.add_argument( + 'domain_name', nargs='?', default='torus', choices=domain_names, + help='gallery domain to build and plot (default: torus)') + parser.add_argument( + '--output', type=Path, + help='optional path of the generated PNG or PDF figure') + parser.add_argument( + '--no-show', action='store_true', + help='create the plot without opening an interactive window') + parser.add_argument( + '--topology', action='store_true', + help='add patch, interface, orientation, and vertex annotations') + arguments = parser.parse_args(argv) + + return plot_multipatch_domain( + arguments.domain_name, + output=arguments.output, + show=not arguments.no_show, + topology=arguments.topology, + ) + + +if __name__ == '__main__': + main() diff --git a/sympde/topology/tests/test_gallery.py b/sympde/topology/tests/test_gallery.py index ffca551c..05dbb483 100644 --- a/sympde/topology/tests/test_gallery.py +++ b/sympde/topology/tests/test_gallery.py @@ -1,8 +1,15 @@ +import numpy as np +import h5py +import pytest +import yaml + from sympy.abc import x,y,z from sympy import Tuple from sympy import symbols +from sympde.topology.analytical_mapping import IdentityMapping from sympde.topology.basic import Interval, ProductDomain, InteriorDomain +from sympde.topology.basic import Boundary from sympde.topology.domain import Domain from sympde.topology.domain import Line, Square, Cube, NCubeInterior @@ -226,6 +233,241 @@ def test_orthogonal_hexahedron(tmp_path): D = Domain.from_file(filename) assert D == domain + +#============================================================================== +def test_multipatch_orientation_roundtrip(tmp_path): + domain = Domain.join( + [Cube('A'), Cube('B')], + [((0, 0, 1), (1, 1, -1), (-1, +1, -1))], + 'AB') + filename = tmp_path / 'multipatch.h5' + + domain.export(str(filename)) + restored = Domain.from_file(str(filename)) + + assert restored.interfaces.orientation == (-1, +1, -1) + assert restored.interfaces.minus.axis == 0 + assert restored.interfaces.plus.axis == 1 + assert restored.interfaces.axis_map == ((1, 2, 1), (2, 0, -1)) + + +#============================================================================== +def test_multipatch_file_requires_explicit_orientation(tmp_path): + domain = Domain.join( + [Square('A'), Square('B')], + [((0, 0, 1), (1, 0, -1), +1)], + 'AB') + topology = domain.todict() + for interface_data in topology['connectivity'].values(): + interface_data.pop() + + filename = tmp_path / 'missing_orientation.h5' + with h5py.File(filename, mode='w') as h5: + h5['topology.yml'] = np.array( + yaml.safe_dump(topology, sort_keys=None), dtype='S') + + with pytest.raises(ValueError, match='explicit orientation'): + Domain.from_file(str(filename)) + + +#============================================================================== +def test_multipatch_file_rejects_legacy_orientation_mapping(tmp_path): + domain = Domain.join( + [Square('A'), Square('B')], + [((0, 0, 1), (1, 0, -1), +1)], + 'AB') + topology = domain.todict() + for interface_data in topology['connectivity'].values(): + interface_data[2] = { + 'permutation': [0], + 'directions': [1], + } + + filename = tmp_path / 'legacy_orientation.h5' + with h5py.File(filename, mode='w') as h5: + h5['topology.yml'] = np.array( + yaml.safe_dump(topology, sort_keys=None), dtype='S') + + with pytest.raises(TypeError, match='2D interface orientation'): + Domain.from_file(str(filename)) + + +#============================================================================== +@pytest.mark.parametrize('location, message', [ + ('boundary-axis', 'serialized boundary axis must be an integer'), + ('interface-extremity', + 'serialized interface extremity must be an integer'), +]) +def test_multipatch_file_rejects_non_integer_side_metadata( + tmp_path, location, message): + domain = Domain.join( + [Square('A'), Square('B')], + [((0, 0, +1), (1, 0, -1), +1)], + 'AB') + topology = domain.todict() + + if location == 'boundary-axis': + topology['boundary'][0]['axis'] = 0.0 + else: + interface_data = next(iter(topology['connectivity'].values())) + interface_data[0]['ext'] = 1.0 + + filename = tmp_path / f'{location}.h5' + with h5py.File(filename, mode='w') as h5: + h5['topology.yml'] = np.array( + yaml.safe_dump(topology, sort_keys=None), dtype='S') + + with pytest.raises(TypeError, match=message): + Domain.from_file(str(filename)) + + +#============================================================================== +def test_self_interface_roundtrip(tmp_path): + patch = Square('A') + domain = Domain.join( + [patch], + [((0, 0, -1), (0, 0, 1), +1)], + 'periodic-A') + filename = tmp_path / 'self_interface.h5' + + domain.export(str(filename)) + restored = Domain.from_file(str(filename)) + + assert restored.todict() == domain.todict() + assert restored.interfaces.minus.domain == restored.interfaces.plus.domain + + +#============================================================================== +@pytest.mark.parametrize('mapped', [False, True]) +def test_boundary_free_cube_roundtrip(tmp_path, mapped): + logical_patch = Cube('A') + patch = (IdentityMapping('F', dim=3)(logical_patch) + if mapped else logical_patch) + domain = Domain.join( + [patch], + [ + ((0, 0, -1), (0, 0, +1), (+1, +1, +1)), + ((0, 1, -1), (0, 1, +1), (+1, +1, +1)), + ((0, 2, -1), (0, 2, +1), (+1, +1, +1)), + ], + 'periodic-cube') + filename = tmp_path / f'periodic_cube_{mapped}.h5' + + assert domain.boundary is None + assert domain.todict()['boundary'] == [] + + domain.export(str(filename)) + restored = Domain.from_file(str(filename)) + + assert restored.boundary is None + assert restored.todict() == domain.todict() + assert tuple(restored.interface_map) == tuple(domain.interface_map) + assert all((interface.mapping is not None) == mapped + for interface in restored.interface_map.values()) + assert all((interface.logical_domain is not None) == mapped + for interface in restored.interface_map.values()) + + for legacy_value in ('missing', None): + topology = domain.todict() + if legacy_value == 'missing': + topology.pop('boundary') + else: + topology['boundary'] = None + legacy_filename = tmp_path / f'periodic_cube_{legacy_value}.h5' + with h5py.File(legacy_filename, mode='w') as h5: + h5['topology.yml'] = np.array( + yaml.safe_dump(topology, sort_keys=None), dtype='S') + legacy = Domain.from_file(str(legacy_filename)) + assert legacy.boundary is None + + +#============================================================================== +def test_single_boundary_and_legacy_dictionary_roundtrip(tmp_path): + # Every interface of a joined domain uses two faces, so a domain with a + # single exterior side is a topological one + patch = InteriorDomain('A', dim=1) + domain = Domain('half-open-line', interiors=[patch], + boundaries=[Boundary('Gamma_1', patch, axis=0, ext=-1)]) + canonical_filename = tmp_path / 'single_boundary.h5' + legacy_filename = tmp_path / 'legacy_single_boundary.h5' + topology = domain.todict() + + assert isinstance(topology['boundary'], list) + assert len(topology['boundary']) == 1 + + domain.export(str(canonical_filename)) + canonical = Domain.from_file(str(canonical_filename)) + assert canonical.todict() == domain.todict() + + topology['boundary'] = topology['boundary'][0] + with h5py.File(legacy_filename, mode='w') as h5: + h5['topology.yml'] = np.array( + yaml.safe_dump(topology, sort_keys=None), dtype='S') + + legacy = Domain.from_file(str(legacy_filename)) + assert legacy.todict() == domain.todict() + + +#============================================================================== +def test_missing_boundary_is_derived_from_connectivity(tmp_path): + domain = Domain.join( + [Line('A'), Line('B')], + [((0, 0, +1), (1, 0, -1), None)], + 'line') + + for legacy_value in ('missing', None): + topology = domain.todict() + if legacy_value == 'missing': + topology.pop('boundary') + else: + topology['boundary'] = None + incomplete_filename = tmp_path / f'boundary_{legacy_value}.h5' + with h5py.File(incomplete_filename, mode='w') as h5: + h5['topology.yml'] = np.array( + yaml.safe_dump(topology, sort_keys=None), dtype='S') + incomplete = Domain.from_file(str(incomplete_filename)) + assert incomplete.todict() == domain.todict() + + +#============================================================================== +def test_explicit_interface_name_roundtrip(tmp_path): + patch_a = IdentityMapping('FA', dim=2)(Square('A')) + patch_b = IdentityMapping('FB', dim=2)(Square('B')) + domain = Domain.join( + [patch_a, patch_b], + [ + ((0, 0, -1), (1, 0, -1), +1, 'lower-seam'), + ((0, 0, +1), (1, 0, +1), -1, 'upper-seam'), + ], + 'named-seams') + filename = tmp_path / 'named_interfaces.h5' + + domain.export(str(filename)) + restored = Domain.from_file(str(filename)) + + assert set(restored.interface_map) == {'lower-seam', 'upper-seam'} + assert restored.todict() == domain.todict() + for name, interface in restored.interface_map.items(): + assert str(interface.logical_domain.name) == name + + +#============================================================================== +def test_file_boundary_must_match_connectivity(tmp_path): + domain = Domain.join( + [Square('A'), Square('B')], + [((0, 0, +1), (1, 0, -1), +1)], + 'AB') + topology = domain.todict() + topology['boundary'].pop() + filename = tmp_path / 'inconsistent_boundary.h5' + + with h5py.File(filename, mode='w') as h5: + h5['topology.yml'] = np.array( + yaml.safe_dump(topology, sort_keys=None), dtype='S') + + with pytest.raises(ValueError, match='does not match'): + Domain.from_file(str(filename)) + #============================================================================== # CLEAN UP SYMPY NAMESPACE #============================================================================== diff --git a/sympde/topology/tests/test_logical_expr.py b/sympde/topology/tests/test_logical_expr.py index f3aefebe..f7484a75 100644 --- a/sympde/topology/tests/test_logical_expr.py +++ b/sympde/topology/tests/test_logical_expr.py @@ -193,7 +193,7 @@ def test_logical_expr_2d_2(): D2 = M2(B) patches = [D1, D2] - connectivity = [((0, 0, 1), (1, 0, -1), 1)] + connectivity = [((0, 0, 1), (1, 0, -1), +1)] domain = Domain.join(patches, connectivity, 'domain') @@ -233,6 +233,24 @@ def test_logical_expr_2d_2(): assert str(expr) == 'Integral(A|B, Dot((Jacobian(M2)**(-1)).T * Grad(PlusInterfaceOperator(u2)), (Jacobian(M2)**(-1)).T * Grad(PlusInterfaceOperator(v2)))*sqrt(det(Jacobian(M1|M2).T * Jacobian(M1|M2))))' +#============================================================================== +def test_logical_expr_on_mapped_multipatch_interface(): + logical_domain = Domain.join( + [Square('A'), Square('B')], + [((0, 0, +1), (1, 1, -1), -1)], + 'AB') + mapped_domain = Mapping('M', dim=2)(logical_domain) + V = ScalarFunctionSpace('V', mapped_domain, kind='h1') + u = element_of(V, name='u') + + expr = LogicalExpr(integral(mapped_domain.interfaces, u**2), mapped_domain) + + assert expr.domain is logical_domain.interfaces + assert str(expr) == ( + 'Integral(A|B, u**2*sqrt(det(' + 'Jacobian(M|M).T * Jacobian(M|M))))') + + def test_logical_expr_2d_3(): dim = 2 @@ -247,7 +265,7 @@ def test_logical_expr_2d_3(): patches = [D1, D2] - connectivity = [((0, 0, 1), (1, 0, -1), 1)] + connectivity = [((0, 0, 1), (1, 0, -1), +1)] domain = Domain.join(patches, connectivity, 'domain') V = VectorFunctionSpace('V', domain, kind='hcurl') diff --git a/sympde/topology/tests/test_multipatch_gallery.py b/sympde/topology/tests/test_multipatch_gallery.py new file mode 100644 index 00000000..5e63861f --- /dev/null +++ b/sympde/topology/tests/test_multipatch_gallery.py @@ -0,0 +1,307 @@ +import numpy as np +import pytest + +from sympde.topology.analytical_mapping import TransposedPolarMapping +from sympde.topology.multipatch_gallery import ( + build_annulus_3, + build_annulus_4, + build_cartesian_multipatch_domain_2d, + build_curved_l_shape, + build_multipatch_domain_2d, + build_multipatch_domain_3d, + build_oriented_two_cube, + build_pretzel, + build_pretzel_annulus, + build_pretzel_debug, + build_pretzel_f, + build_self_connected_patch, + build_square_2, + build_square_4, + build_square_6, + build_square_8, + build_square_9, + build_three_patch_shared_edge, + build_three_patch_shared_vertex, + build_torus_3d, + build_two_cube_3d, + build_two_patch_twisted_strip, + main, + plot_multipatch_domain, +) + + +def _physical_point(side, parameter): + patch = side.domain + logical_patch = patch.logical_domain + logical_point = [ + minimum + parameter * (maximum - minimum) + for minimum, maximum in zip( + logical_patch.min_coords, logical_patch.max_coords) + ] + logical_point[side.axis] = ( + logical_patch.min_coords[side.axis] + if side.ext == -1 else logical_patch.max_coords[side.axis] + ) + substitutions = dict(zip( + patch.mapping.logical_coordinates, logical_point)) + return np.array([ + float(expression.subs(substitutions)) + for expression in patch.mapping.expressions + ]) + + +@pytest.mark.parametrize( + 'builder, name, npatches, ninterfaces, nboundaries', + [ + (build_square_2, 'square_2', 2, 1, 6), + (build_square_4, 'square_4', 4, 4, 8), + (build_square_6, 'square_6', 6, 7, 10), + (build_square_8, 'square_8', 8, 8, 16), + (build_square_9, 'square_9', 9, 12, 12), + (build_annulus_3, 'annulus_3', 3, 3, 6), + (build_annulus_4, 'annulus_4', 4, 4, 8), + (build_curved_l_shape, 'curved_L_shape', 3, 2, 8), + (build_pretzel, 'pretzel', 11, 13, 18), + (build_pretzel_f, 'pretzel_f', 18, 20, 32), + (build_pretzel_annulus, 'pretzel_annulus', 9, 9, 18), + (build_pretzel_debug, 'pretzel_debug', 2, 1, 6), + ], +) +def test_multipatch_domain_builder( + builder, name, npatches, ninterfaces, nboundaries): + domain = builder() + + assert str(domain.name) == name + assert len(domain.patches) == npatches + assert len(domain.interface_map) == ninterfaces + assert len(domain.exterior_sides) == nboundaries + assert all(interface.orientation == +1 + for interface in domain.interface_map.values()) + for interface in domain.interface_map.values(): + for parameter in (0, 0.5, 1): + np.testing.assert_allclose( + _physical_point(interface.minus, parameter), + _physical_point(interface.plus, parameter), + atol=1e-14, + ) + + +@pytest.mark.parametrize('builder', [ + build_square_2, + build_square_6, + build_square_8, + build_square_9, + build_annulus_3, + build_annulus_4, + build_curved_l_shape, + build_pretzel, + build_pretzel_f, + build_pretzel_annulus, + build_pretzel_debug, +]) +def test_numbered_gallery_domains_use_patch_names(builder): + names = tuple(str(patch.logical_domain.name) for patch in builder().patches) + + assert len(names) == len(set(names)) + assert all(name.startswith('P') for name in names) + + +@pytest.mark.parametrize('mapping', ['identity', 'polar']) +def test_build_rectangular_cartesian_multipatch_domain(mapping): + layout = np.array([ + [1, None, 5], + [2, 3, 4], + ], dtype=object) + + domain = build_cartesian_multipatch_domain_2d( + layout, (0, 3), (0, 2), mapping=mapping) + + assert tuple(str(p.logical_domain.name) for p in domain.patches) == ( + 'Log_0_0', 'Log_0_2', 'Log_1_0', 'Log_1_1', 'Log_1_2') + assert len(domain.interface_map) == 4 + assert len(domain.exterior_sides) == 12 + + patches = {str(p.logical_domain.name): p for p in domain.patches} + top_left = patches['Log_0_0'].logical_domain + bottom_right = patches['Log_1_2'].logical_domain + assert top_left.min_coords == (0, 1) + assert top_left.max_coords == (1, 2) + assert bottom_right.min_coords == (2, 0) + assert bottom_right.max_coords == (3, 1) + + +@pytest.mark.parametrize('layout, message', [ + ([1, 2], 'two-dimensional'), + (np.empty((0, 2)), 'at least one layout position'), + ([[None]], 'at least one non-None entry'), +]) +def test_build_cartesian_multipatch_domain_2d_rejects_invalid_layout( + layout, message): + with pytest.raises(ValueError, match=message): + build_cartesian_multipatch_domain_2d(layout, (0, 1), (0, 1)) + + +def test_build_cartesian_multipatch_domain_2d_rejects_unknown_mapping(): + with pytest.raises(ValueError, match="'identity' or 'polar'"): + build_cartesian_multipatch_domain_2d( + [[1]], (0, 1), (0, 1), mapping='unknown') + + +def test_transposed_polar_mapping_is_public(): + mapping = TransposedPolarMapping( + 'F', dim=2, c1=0, c2=0, rmin=0, rmax=1) + + angle, radius = mapping.logical_coordinates + x, y = mapping.expressions + assert x.subs({angle: 0, radius: 2}) == 2 + assert y.subs({angle: 0, radius: 2}) == 0 + + +@pytest.mark.parametrize('name', ['unknown', None]) +def test_build_multipatch_domain_2d_rejects_unknown_name(name): + with pytest.raises(ValueError, match='unknown 2D multipatch domain'): + build_multipatch_domain_2d(name) + + +@pytest.mark.parametrize('name', ['annulus_4', 'pretzel']) +def test_build_multipatch_domain_2d_rejects_invalid_radii(name): + with pytest.raises(ValueError, match='0 < r_min < r_max'): + build_multipatch_domain_2d(name, r_min=2, r_max=1) + + +def test_build_multipatch_domain_2d_dispatches_names_and_radii(): + square = build_multipatch_domain_2d('square_2') + annulus = build_multipatch_domain_2d( + 'annulus_3', r_min=1, r_max=3) + + assert str(square.name) == 'square_2' + assert str(annulus.name) == 'annulus_3' + assert annulus.patches[0].logical_domain.min_coords[0] == 1 + assert annulus.patches[0].logical_domain.max_coords[0] == 3 + + +def test_build_two_cube_3d(): + domain = build_two_cube_3d() + + assert str(domain.name) == 'mapped_two_cube_3d' + assert len(domain.patches) == 2 + assert len(domain.interface_map) == 1 + assert len(domain.exterior_sides) == 10 + + interface, = domain.interface_map.values() + assert interface.orientation == (-1, +1, -1) + assert interface.axis_map == ((1, 2, +1), (2, 0, -1)) + + +@pytest.mark.parametrize( + 'arguments, name, ninterfaces, nboundaries', + [ + ({}, 'hollow_torus_3d', 8, 8), + ( + {'toroidal_angle': 1.5 * np.pi}, + 'open_hollow_torus_3d', 6, 12, + ), + ( + {'hollow': False, 'close_torus': False}, + 'open_solid_torus_3d', 6, 12, + ), + ], +) +def test_build_torus_3d( + arguments, name, ninterfaces, nboundaries): + domain = build_torus_3d(**arguments) + + assert str(domain.name) == name + assert len(domain.patches) == 4 + assert len(domain.interface_map) == ninterfaces + assert len(domain.exterior_sides) == nboundaries + assert all(interface.orientation == (+1, +1, +1) + for interface in domain.interface_map.values()) + + +@pytest.mark.parametrize('arguments, error, message', [ + ({'minor_bounds': (0.8, 0.35)}, ValueError, 'minor radius'), + ({'hollow': 1}, TypeError, 'hollow must be a bool'), + ({'toroidal_angle': 0}, ValueError, '0 < angle <= 2\\*pi'), + ( + {'toroidal_angle': np.pi, 'close_torus': True}, + ValueError, + 'only a full 2\\*pi sweep can be closed', + ), +]) +def test_build_torus_3d_rejects_invalid_parameters( + arguments, error, message): + with pytest.raises(error, match=message): + build_torus_3d(**arguments) + + +@pytest.mark.parametrize( + 'builder, name, npatches, ninterfaces, nboundaries', + [ + ( + build_two_patch_twisted_strip, + 'two_patch_twisted_strip', 2, 2, 4, + ), + (build_self_connected_patch, 'self_connected_patch', 1, 1, 2), + ( + build_three_patch_shared_vertex, + 'three_patch_shared_vertex', 3, 3, 6, + ), + (build_oriented_two_cube, 'oriented_two_cube', 2, 1, 10), + ( + build_three_patch_shared_edge, + 'three_patch_shared_edge', 3, 3, 12, + ), + ], +) +def test_diagnostic_domain_builder( + builder, name, npatches, ninterfaces, nboundaries): + domain = builder() + + assert str(domain.name) == name + assert len(domain.patches) == npatches + assert len(domain.interface_map) == ninterfaces + assert len(domain.exterior_sides) == nboundaries + + +@pytest.mark.parametrize('name', ['unknown', None]) +def test_build_multipatch_domain_3d_rejects_unknown_name(name): + with pytest.raises(ValueError, match='unknown 3D multipatch domain'): + build_multipatch_domain_3d(name) + + +def test_build_multipatch_domain_3d_dispatches_names_and_arguments(): + two_cube = build_multipatch_domain_3d('two_cube') + torus = build_multipatch_domain_3d( + 'torus', hollow=False, close_torus=False) + + assert str(two_cube.name) == 'mapped_two_cube_3d' + assert str(torus.name) == 'open_solid_torus_3d' + + +@pytest.mark.parametrize('domain_name, dimension', [ + ('square_2', 2), + ('two_cube', 3), +]) +def test_plot_multipatch_domain(domain_name, dimension, tmp_path): + output = tmp_path / f'{domain_name}.png' + domain, figure = plot_multipatch_domain( + domain_name, output=output, show=False, topology=True) + + assert int(domain.dim) == dimension + assert figure.axes + assert output.is_file() + + +def test_multipatch_gallery_main(tmp_path): + output = tmp_path / 'square_2.png' + + main([ + 'square_2', + '--topology', + '--no-show', + '--output', + str(output), + ]) + + assert output.is_file() diff --git a/sympde/topology/tests/test_topology.py b/sympde/topology/tests/test_topology.py index ddbf725a..4358cf4c 100644 --- a/sympde/topology/tests/test_topology.py +++ b/sympde/topology/tests/test_topology.py @@ -1,4 +1,7 @@ import os +import copy +import numpy as np +import pytest import h5py import numpy as np @@ -7,14 +10,14 @@ from sympy import preorder_traversal from sympy.core import Basic +from sympde.topology.analytical_mapping import IdentityMapping from sympde.topology.basic import InteriorDomain, Union from sympde.topology.basic import Boundary -from sympde.topology.domain import Domain, ElementDomain -from sympde.topology.domain import Area -from sympde.topology.mapping import Mapping -from sympde.topology.basic import Interface +from sympde.topology.basic import Interface, PatchVertex, SharedVertex +from sympde.topology.domain import Area, Domain, ElementDomain from sympde.topology.domain import Line, Square, Cube, NCube -from sympde.topology.analytical_mapping import IdentityMapping +from sympde.topology.mapping import InterfaceMapping, Mapping, MultiPatchMapping + base_dir = os.path.dirname(os.path.realpath(__file__)) topo_dir = os.path.join(base_dir, 'data') @@ -23,13 +26,13 @@ def test_structural_arguments_are_sympy_objects(): domain = Domain('Omega', dim=2) boundary = Boundary(r'\Gamma_1', domain) - interface = Interface('I', boundary, boundary, ornt=1) + interface = Interface('I', boundary, boundary, 1) square = Square('S', bounds1=(-1, 2), bounds2=(3, 4)) assert domain.name == 'Omega' assert boundary.name == r'\Gamma_1' assert boundary.axis is None - assert interface.ornt == 1 + assert interface.orientation == 1 for expr in (domain, domain.interior, boundary, interface, square, square.interior, square.boundary.args[0]): @@ -103,7 +106,7 @@ def test_topology_1(tmp_path): D2 = M2(B) patches = [D1, D2] - connectivity = [((0, 0, 1), (1, 0, -1), 1)] + connectivity = [((0, 0, 1), (1, 0, -1), +1)] Omega = Domain.join(patches, connectivity, 'domain') interfaces = Omega.interfaces @@ -118,6 +121,47 @@ def test_topology_1(tmp_path): assert D.todict() == Omega.todict() +#============================================================================== +def test_topological_domain_export_and_from_file(tmp_path): + # Regression test for #144: the patches are InteriorDomain objects without + # a dtype, and their boundaries have no axis or extremity + A = InteriorDomain('A', dim=2) + B = InteriorDomain('B', dim=2) + A1, A2, A3 = [Boundary(f'Gamma_{i}', A) for i in (1, 2, 3)] + B1, B2, B3 = [Boundary(f'Gamma_{i}', B) for i in (1, 2, 3)] + + Omega = Domain('Omega', + interiors=[A, B], + boundaries=[A2, A3, B1, B3], + connectivity=[Interface('I', A1, B2, 1)]) + + filename = str(tmp_path / 'omega.h5') + Omega.export(filename) + D = Domain.from_file(filename) + + assert D == Omega + assert D.todict() == Omega.todict() + + # 3D patches with known boundary axes and extremities + C = InteriorDomain('C', dim=3) + E = InteriorDomain('E', dim=3) + interface = Interface('I', + Boundary('Gamma_1', C, axis=0, ext=1), + Boundary('Gamma_2', E, axis=1, ext=-1), + (-1, 1, -1)) + Omega = Domain('Omega', interiors=[C, E], boundaries=[], + connectivity=[interface]) + + Omega.export(filename) + D = Domain.from_file(filename) + + assert D == Omega + assert D.interface_map['I'].orientation == (-1, 1, -1) + + # Topological patches have no parametric box, hence no vertices + with pytest.raises(ValueError, match='parametric box'): + _ = D.shared_vertices + #============================================================================== def test_export_surface_mapping(tmp_path): domain = Mapping('F', ldim=2, pdim=3)(Square('Omega')) @@ -244,7 +288,7 @@ def test_domain_join_line(): domains = [A, B] - connectivity = [((0, 0, 1), (1, 0, -1))] + connectivity = [((0, 0, 1), (1, 0, -1), None)] AB = Domain.join(domains, connectivity, 'AB') AB_bnd_minus = A.get_boundary(axis=0, ext=1) @@ -254,17 +298,20 @@ def test_domain_join_line(): BC_bnd_plus = C.get_boundary(axis=0, ext=-1) assert AB.interior == Union(A.interior, B.interior) - assert AB.interfaces == Interface('A|B', AB_bnd_minus, AB_bnd_plus) + assert AB.interfaces == Interface( + 'A|B', AB_bnd_minus, AB_bnd_plus, None) print(AB.connectivity) print('') domains = [A, B, C] - connectivity = [((0, 0, 1), (1, 0, -1)), ((1, 0, 1), (2, 0, -1))] + connectivity = [((0, 0, 1), (1, 0, -1), None), + ((1, 0, 1), (2, 0, -1), None)] ABC = Domain.join(domains, connectivity, 'ABC') assert ABC.interior == Union(A.interior, B.interior, C.interior) - assert ABC.interfaces == Union(Interface('A|B', AB_bnd_minus, AB_bnd_plus), - Interface('B|C', BC_bnd_minus, BC_bnd_plus)) + assert ABC.interfaces == Union( + Interface('A|B', AB_bnd_minus, AB_bnd_plus, None), + Interface('B|C', BC_bnd_minus, BC_bnd_plus, None)) print(list(ABC.connectivity.items())) print('') @@ -278,7 +325,7 @@ def test_domain_join_square(): C = Square('C') patches = [A, B] - connectivity = [((0, 0, 1), (1, 0, -1), 1)] + connectivity = [((0, 0, 1), (1, 0, -1), +1)] AB = Domain.join(patches, connectivity, 'AB') AB_bnd_minus = A.get_boundary(axis=0, ext= 1) @@ -289,19 +336,21 @@ def test_domain_join_square(): print(AB) assert AB.interior == Union(A.interior, B.interior) - assert AB.interfaces == Interface('A|B', AB_bnd_minus, AB_bnd_plus, ornt=1) + assert AB.interfaces == Interface( + 'A|B', AB_bnd_minus, AB_bnd_plus, +1) print(AB.connectivity) patches = [A, B, C] - connectivity = [((0, 0, 1), (1, 0, -1), 1), - ((1, 0, 1), (2, 0, -1), 1)] + connectivity = [((0, 0, 1), (1, 0, -1), +1), + ((1, 0, 1), (2, 0, -1), +1)] ABC = Domain.join(patches, connectivity, 'ABC') print(ABC) assert ABC.interior == Union(A.interior, B.interior, C.interior) - assert ABC.interfaces == Union(Interface('A|B', AB_bnd_minus, AB_bnd_plus, ornt=1), - Interface('B|C', BC_bnd_minus, BC_bnd_plus, ornt=1)) + assert ABC.interfaces == Union( + Interface('A|B', AB_bnd_minus, AB_bnd_plus, +1), + Interface('B|C', BC_bnd_minus, BC_bnd_plus, +1)) print(list(ABC.connectivity.items())) print('') @@ -313,11 +362,12 @@ def test_get_subdomain(): C = Square('C') patches = [A, B] - connectivity = [((0, 0, 1), (1, 0, -1), 1)] + connectivity = [((0, 0, 1), (1, 0, -1), +1)] AB = Domain.join(patches, connectivity, 'AB') patches = [A, B, C] - connectivity = [((0, 0, 1), (1, 0, -1), 1), ((1, 0, 1), (2, 0, -1), 1)] + connectivity = [((0, 0, 1), (1, 0, -1), +1), + ((1, 0, 1), (2, 0, -1), +1)] ABC = Domain.join(patches, connectivity, 'ABC') A_1 = AB.get_subdomain('A') @@ -338,6 +388,147 @@ def test_get_subdomain(): assert A_pipe_B_pipe_C_2 is ABC + +#============================================================================== +def test_get_subdomain_of_three_from_four_patches(): + patches = [Square(name) for name in 'ABCD'] + connectivity = [ + ((0, 0, +1), (1, 0, -1), +1), + ((1, 0, +1), (2, 0, -1), +1), + ((2, 0, +1), (3, 0, -1), +1), + ] + domain = Domain.join(patches, connectivity, 'ABCD') + subdomain = domain.get_subdomain(('A', 'B', 'C')) + expected = Domain.join(patches[:3], connectivity[:2], 'expected') + + assert subdomain.name == 'A|B|C' + assert subdomain.patches == expected.patches + assert subdomain.boundary == expected.boundary + assert tuple(subdomain.interface_map) == ('A|B', 'B|C') + assert subdomain.interface_map['A|B'] is domain.interface_map['A|B'] + assert subdomain.interface_map['B|C'] is domain.interface_map['B|C'] + + +#============================================================================== +def test_get_subdomain_preserves_repeated_and_self_interfaces(): + patch_a = Square('A') + patch_b = Square('B') + patch_c = Square('C') + domain = Domain.join( + [patch_a, patch_b, patch_c], + [ + ((0, 1, -1), (0, 1, +1), +1), + ((0, 0, -1), (1, 0, -1), +1), + ((0, 0, +1), (1, 0, +1), +1), + ((1, 1, +1), (2, 1, -1), +1), + ], + 'ABC') + + patch_subdomain = domain.get_subdomain('A') + pair_subdomain = domain.get_subdomain(('A', 'B')) + + assert tuple(patch_subdomain.interface_map) == ('A|A',) + assert patch_subdomain.interface_map['A|A'] is \ + domain.interface_map['A|A'] + assert { + (side.domain.name, side.axis, side.ext) + for side in patch_subdomain.exterior_sides + } == {('A', 0, -1), ('A', 0, +1)} + + assert tuple(pair_subdomain.interface_map) == ( + 'A|A', 'A|B', 'A|B#2') + assert all( + pair_subdomain.interface_map[name] is domain.interface_map[name] + for name in pair_subdomain.interface_map) + assert ('B', 1, +1) in { + (side.domain.name, side.axis, side.ext) + for side in pair_subdomain.exterior_sides + } + + +#============================================================================== +def test_get_subdomain_preserves_mapped_3d_topology(): + logical_patches = [Cube(name) for name in 'ABCD'] + patches = [ + IdentityMapping(f'F{name}', dim=3)(patch) + for name, patch in zip('ABCD', logical_patches) + ] + domain = Domain.join( + patches, + [ + ((0, 0, +1), (1, 1, -1), (-1, +1, -1)), + ((1, 2, +1), (2, 0, -1), (+1, -1, +1)), + ((2, 1, +1), (3, 2, -1), (-1, -1, +1)), + ], + 'mapped-ABCD') + selected_names = tuple(patch.interior.name for patch in patches[:3]) + subdomain = domain.get_subdomain(selected_names) + retained_interfaces = tuple(domain.interface_map.values())[:2] + cut_side = tuple(domain.interface_map.values())[2].minus + + assert isinstance(subdomain.mapping, MultiPatchMapping) + assert subdomain.logical_domain is not None + assert subdomain.patches == tuple(patch.interior for patch in patches[:3]) + assert tuple(subdomain.interface_map.values()) == retained_interfaces + assert cut_side in subdomain.exterior_sides + assert tuple( + interface.axis_map for interface in subdomain.interface_map.values() + ) == tuple(interface.axis_map for interface in retained_interfaces) + assert set(subdomain.mapping.mappings) == { + patch.interior.logical_domain for patch in patches[:3]} + + for interface in subdomain.interface_map.values(): + logical_interface = interface.logical_domain + assert logical_interface is not None + assert subdomain.logical_domain.interface_map[ + str(logical_interface.name)] is logical_interface + + +#============================================================================== +def test_domain_join_rejects_nested_multipatch_inputs(): + patch_a = Square('A') + patch_b = Square('B') + patch_c = Square('C') + joined = Domain.join( + [patch_a, patch_b], + [((0, 0, +1), (1, 0, -1), +1)], + 'AB') + + with pytest.raises(ValueError, match='atomic patch domains'): + Domain.join([joined, patch_c], [], 'ABC') + +#============================================================================== +def test_domain_join_rejects_different_physical_dimensions(): + # Regression test for #140: a planar patch joined with a surface in 3D + A = Mapping('FA', ldim=2, pdim=2)(Square('A')) + B = Mapping('FB', ldim=2, pdim=3)(Square('B')) + C = Square('C') + + with pytest.raises(ValueError, match='same physical dimension'): + Domain.join([A, B], [((0, 0, 1), (1, 0, -1), 1)], 'Omega') + with pytest.raises(ValueError, match='same physical dimension'): + Domain.join([B, C], [((0, 0, 1), (1, 0, -1), 1)], 'Omega') + +#============================================================================== +def test_domain_join_rejects_mapped_and_unmapped_patches(): + A = Mapping('F', dim=2)(Square('A')) + B = Square('B') + + with pytest.raises(ValueError, match='all patches or none'): + Domain.join([A, B], [((0, 0, 1), (1, 0, -1), 1)], 'Omega') + +#============================================================================== +def test_domain_join_rejects_degenerate_connectivity(): + A = Square('A') + B = Square('B') + + with pytest.raises(ValueError, match='with itself'): + Domain.join([A], [((0, 0, 1), (0, 0, 1), 1)], 'Omega') + + with pytest.raises(ValueError, match='more than one interface'): + Domain.join([A, B], [((0, 0, 1), (1, 0, -1), 1), + ((0, 0, 1), (1, 1, -1), 1)], 'Omega') + #============================================================================== def test_2d_domain_without_bnd(): @@ -355,14 +546,14 @@ def test_2d_domain_without_bnd(): domain_4 = mapping_4(OmegaLog4) patches = [domain_1, domain_2, domain_3, domain_4] - connectivity = [((0, 0, 1), (2, 0,-1), 1), - ((1, 0, 1), (3, 0,-1), 1), - ((2, 0, 1), (0, 0,-1), 1), - ((3, 0, 1), (1, 0,-1), 1), - ((0, 1, 1), (1, 1,-1), 1), - ((2, 1, 1), (3, 1,-1), 1), - ((1, 1, 1), (0, 1,-1), 1), - ((3, 1, 1), (2, 1,-1), 1)] + connectivity = [((0, 0, 1), (2, 0,-1), +1), + ((1, 0, 1), (3, 0,-1), +1), + ((2, 0, 1), (0, 0,-1), +1), + ((3, 0, 1), (1, 0,-1), +1), + ((0, 1, 1), (1, 1,-1), +1), + ((2, 1, 1), (3, 1,-1), +1), + ((1, 1, 1), (0, 1,-1), +1), + ((3, 1, 1), (2, 1,-1), +1)] domain = Domain.join(patches, connectivity, 'domain') assert len(domain.interior) == 4 @@ -387,14 +578,14 @@ def test_3d_domain_with_doubly_connected_patches(): domain_4 = mapping_4(OmegaLog4) patches = [domain_1, domain_2, domain_3, domain_4] - connectivity = [((0, 0, 1), (2, 0,-1), (1, 1, 1)), - ((1, 0, 1), (3, 0,-1), (1, 1, 1)), - ((2, 0, 1), (0, 0,-1), (1, 1, 1)), - ((3, 0, 1), (1, 0,-1), (1, 1, 1)), - ((0, 1, 1), (1, 1,-1), (1, 1, 1)), - ((2, 1, 1), (3, 1,-1), (1, 1, 1)), - ((1, 1, 1), (0, 1,-1), (1, 1, 1)), - ((3, 1, 1), (2, 1,-1), (1, 1, 1))] + connectivity = [((0, 0, 1), (2, 0,-1), (+1, +1, +1)), + ((1, 0, 1), (3, 0,-1), (+1, +1, +1)), + ((2, 0, 1), (0, 0,-1), (+1, +1, +1)), + ((3, 0, 1), (1, 0,-1), (+1, +1, +1)), + ((0, 1, 1), (1, 1,-1), (+1, +1, +1)), + ((2, 1, 1), (3, 1,-1), (+1, +1, +1)), + ((1, 1, 1), (0, 1,-1), (+1, +1, +1)), + ((3, 1, 1), (2, 1,-1), (+1, +1, +1))] domain = Domain.join(patches, connectivity, 'domain') assert len(domain.interior) == 4 @@ -404,6 +595,435 @@ def test_3d_domain_with_doubly_connected_patches(): assert all(isinstance(b, Boundary) for b in domain.boundary) assert len(domain.boundary) == 8 +#============================================================================== +def test_interface_orientation_axis_map(): + domain = Domain.join( + [Cube('A'), Cube('B')], + [((0, 0, +1), (1, 1, -1), (-1, +1, -1))], + 'AB') + interface = domain.interfaces + + assert interface.orientation == (-1, +1, -1) + assert interface.axis_map == ((1, 2, +1), (2, 0, -1)) + + +#============================================================================== +@pytest.mark.parametrize('data, permutation, directions', [ + ((+1, +1, +1), (0, 1), (+1, +1)), + ((+1, +1, -1), (0, 1), (+1, -1)), + ((+1, -1, +1), (0, 1), (-1, +1)), + ((+1, -1, -1), (0, 1), (-1, -1)), + ((-1, +1, +1), (1, 0), (+1, +1)), + ((-1, +1, -1), (1, 0), (+1, -1)), + ((-1, -1, +1), (1, 0), (-1, +1)), + ((-1, -1, -1), (1, 0), (-1, -1)), +]) +def test_public_interface_orientation_codec_3d(data, permutation, directions): + domain = Domain.join( + [Cube('A'), Cube('B')], + [((0, 0, +1), (1, 1, -1), data)], + 'AB') + interface = domain.interfaces + plus_tangents = (0, 2) + + assert interface.orientation == data + assert interface.axis_map == tuple( + (minus_axis, plus_tangents[plus_position], direction) + for minus_axis, plus_position, direction + in zip((1, 2), permutation, directions) + ) + + +#============================================================================== +def test_public_interface_orientation_codec_1d_and_2d(): + line_domain = Domain.join( + [Line('A'), Line('B')], + [((0, 0, +1), (1, 0, -1), None)], + 'lines') + forward = Domain.join( + [Square('C'), Square('D')], + [((0, 0, +1), (1, 0, -1), +1)], + 'forward') + reversed_domain = Domain.join( + [Square('E'), Square('F')], + [((0, 0, +1), (1, 0, -1), -1)], + 'reversed') + + assert line_domain.interfaces.orientation is None + assert line_domain.interfaces.axis_map == () + assert forward.interfaces.orientation == +1 + assert reversed_domain.interfaces.orientation == -1 + + with pytest.raises(ValueError, match='1D'): + Domain.join( + [Line('G'), Line('H')], + [((0, 0, +1), (1, 0, -1), +1)], 'invalid') + with pytest.raises(ValueError, match='2D'): + Domain.join( + [Square('G'), Square('H')], + [((0, 0, +1), (1, 0, -1), 0)], 'invalid') + with pytest.raises(TypeError, match='3D'): + Domain.join( + [Cube('G'), Cube('H')], + [((0, 0, +1), (1, 0, -1), (+1, +1))], 'invalid') + with pytest.raises(ValueError, match='each 3D'): + Domain.join( + [Cube('I'), Cube('J')], + [((0, 0, +1), (1, 0, -1), (+1, 0, +1))], 'invalid') + + +#============================================================================== +@pytest.mark.parametrize('orientation', [1.0, True]) +def test_2d_interface_orientation_rejects_non_integer_values(orientation): + with pytest.raises(TypeError, match='2D interface orientation.*integer'): + Domain.join( + [Square('A'), Square('B')], + [((0, 0, +1), (1, 0, -1), orientation)], + 'invalid') + + +#============================================================================== +@pytest.mark.parametrize('orientation', [ + (1.9, -1.2, 1.0), + (+1, -1, True), +]) +def test_3d_interface_orientation_rejects_non_integer_values(orientation): + with pytest.raises(TypeError, match='3D.*contain integers'): + Domain.join( + [Cube('A'), Cube('B')], + [((0, 0, +1), (1, 0, -1), orientation)], + 'invalid') + + +#============================================================================== +@pytest.mark.parametrize('side, expected', [ + ((0, 0.0, +1), 'axis'), + ((0, False, +1), 'axis'), + ((0, 0, 1.0), 'extremity'), + ((0, 0, True), 'extremity'), +]) +def test_interface_sides_reject_non_integer_axis_and_extremity(side, expected): + with pytest.raises(TypeError, match=rf'{expected} must be an integer'): + Domain.join( + [Square('A'), Square('B')], + [(side, (1, 0, -1), +1)], + 'invalid') + + +#============================================================================== +def test_boundary_lookup_rejects_non_integer_axis_and_extremity(): + patch = Square('A') + + with pytest.raises(TypeError, match='axis must be an integer'): + patch.get_boundary(axis=0.0, ext=+1) + with pytest.raises(TypeError, match='extremity must be an integer'): + patch.get_boundary(axis=0, ext=1.0) + with pytest.raises(TypeError, match='axis must be an integer'): + Boundary('invalid', patch.interior, axis=False, ext=+1) + + +#============================================================================== +def test_domain_join_accepts_omitted_connectivity(): + patch_a = Square('A') + patch_b = Square('B') + + assert Domain.join([patch_a], name='unused') is patch_a + + domain = Domain.join([patch_a, patch_b], name='AB') + assert domain.patches == (patch_a.interior, patch_b.interior) + assert domain.interfaces is None + assert len(domain.exterior_sides) == 8 + + +#============================================================================== +def test_domain_join_accepts_numpy_integer_patch_indices(): + domain = Domain.join( + [Square('A'), Square('B')], + [((np.int64(0), 0, +1), (np.int32(1), 0, -1), +1)], + 'AB') + + interface = domain.interfaces + assert interface.minus.domain.name == 'A' + assert interface.plus.domain.name == 'B' + + +#============================================================================== +def test_domain_join_rejects_duplicate_patch_identity_and_names(): + patch = Square('A') + with pytest.raises(ValueError, match='same object twice'): + Domain.join([patch, patch], [], 'duplicate-object') + + with pytest.raises(ValueError, match='interior names must be unique'): + Domain.join( + [Square('A', bounds1=(0, 1)), + Square('A', bounds1=(1, 2))], + [], + 'duplicate-name') + + logical_patch = Square('logical') + mapped_patches = [ + IdentityMapping(name, dim=2)(logical_patch) + for name in ('F', 'G') + ] + with pytest.raises(ValueError, match='serialization must be unique'): + Domain.join(mapped_patches, [], 'duplicate-logical-name') + + +#============================================================================== +def test_named_multipatch_api_preserves_single_patch_domain_api(): + single_patch = Square('single') + + assert Domain.join([single_patch], [], 'unused') is single_patch + assert single_patch.patches == (single_patch.interior,) + assert single_patch.interface_map is single_patch.connectivity + assert tuple(single_patch.interface_map.values()) == () + assert single_patch.exterior_sides == single_patch.boundary.as_tuple() + assert single_patch.shared_vertices == () + + patch_a = Square('A') + patch_b = Square('B') + domain = Domain.join( + [patch_a, patch_b], + interfaces=[( + (patch_a, 0, +1), + (patch_b, 1, -1), + -1, + )], + name='AB') + + interface, = domain.interface_map.values() + assert domain.patches == (patch_a.interior, patch_b.interior) + assert interface is domain.interfaces + assert interface.minus_side is interface.minus + assert interface.plus_side is interface.plus + assert interface.minus_side.patch is patch_a.interior + assert interface.minus_side.normal_axis == interface.minus.axis == 0 + assert interface.minus_side.ext == +1 + assert interface.orientation == -1 + assert interface.axis_map == ((1, 0, -1),) + assert interface.patch_dim == 2 + assert interface.intrinsic_dim == 1 + assert len(domain.exterior_sides) == 6 + assert len(domain.shared_vertices) == 2 + assert all(isinstance(vertex, SharedVertex) + for vertex in domain.shared_vertices) + assert all(isinstance(corner, PatchVertex) + for vertex in domain.shared_vertices + for corner in vertex.corners) + + with pytest.raises(TypeError, match='either connectivity or interfaces'): + Domain.join( + [patch_a, patch_b], connectivity=[], interfaces=[], name='AB') + + +#============================================================================== +def test_interface_can_be_reconstructed_from_sympy_args(): + A = Cube('A') + B = Cube('B') + interface = Interface( + 'A|B', A.get_boundary(0, 1), B.get_boundary(1, -1), + (-1, +1, -1)) + + reconstructed = interface.func(*interface.args) + copied = copy.copy(interface) + + assert reconstructed == interface + assert reconstructed.orientation == interface.orientation + assert copied.orientation == interface.orientation + + +#============================================================================== +def test_3d_interface_with_different_face_axes_and_shared_corners(): + A = Cube('A') + B = Cube('B') + domain = Domain.join( + [A, B], + [((0, 0, +1), (1, 1, -1), (-1, +1, -1))], + 'AB') + + interface = domain.interfaces + assert interface.minus.axis == 0 + assert interface.plus.axis == 1 + assert interface.orientation == (-1, +1, -1) + assert interface.axis_map == ((1, 2, 1), (2, 0, -1)) + + corners = domain.corners.as_tuple() + assert len(corners) == 4 + coordinate_pairs = { + tuple(corner.coordinates for corner in shared_corner.corners) + for shared_corner in corners + } + assert ((1, 0, 1), (0, 0, 0)) in coordinate_pairs + + +#============================================================================== +def test_mapping_multipatch_domain_preserves_interface_orientation(): + logical_domain = Domain.join( + [Cube('A'), Cube('B')], + [((0, 0, 1), (1, 1, -1), (-1, +1, -1))], + 'AB') + mapping = Mapping('F', dim=3) + mapped_domain = mapping(logical_domain) + logical_interface = logical_domain.interfaces + mapped_interface = mapped_domain.interfaces + + assert mapped_interface.orientation == logical_interface.orientation + assert mapped_interface.axis_map == logical_interface.axis_map + assert isinstance(mapped_interface.mapping, InterfaceMapping) + assert mapped_interface.logical_domain is logical_interface + assert mapped_interface.minus.logical_domain == logical_interface.minus + assert mapped_interface.plus.logical_domain == logical_interface.plus + assert mapping(logical_interface) is mapped_interface + + +#============================================================================== +def test_generated_interface_names_do_not_collide_with_patch_names(): + patch_a = Square('A') + patch_b = Square('B') + patch_b2 = Square('B#2') + domain = Domain.join( + [patch_a, patch_b, patch_b2], + [ + ((patch_a, 0, -1), (patch_b, 0, -1), +1), + ((patch_a, 0, +1), (patch_b, 0, +1), +1), + ((patch_a, 1, -1), (patch_b2, 1, -1), +1), + ], + 'collision-safe') + + assert tuple(domain.interface_map) == ( + 'A|B', 'A|B#2', 'A|B#2#2') + assert len(domain.interface_map) == 3 + assert domain.interface_map['A|B#2'].plus.domain == patch_b.interior + assert domain.interface_map['A|B#2#2'].plus.domain == patch_b2.interior + + +#============================================================================== +def test_domain_does_not_depend_on_interface_order(tmp_path): + A = Square('A') + B = Square('B') + first = ((0, 0, +1), (1, 0, -1), +1, 'zz') + second = ((0, 1, +1), (1, 1, -1), +1) + + domain = Domain.join([A, B], [first, second], 'Omega') + + assert domain == Domain.join([A, B], [second, first], 'Omega') + assert tuple(domain.interface_map) == ('A|B', 'zz') + + # Generated names do not depend on the order of the interfaces either + unnamed = first[:3] + assert Domain.join([A, B], [unnamed, second], 'Omega') == \ + Domain.join([A, B], [second, unnamed], 'Omega') + + filename = str(tmp_path / 'omega.h5') + domain.export(filename) + assert Domain.from_file(filename) == domain + + # A mapped interface and its logical interface have the same name, also + # after an export/from_file round trip + F = Mapping('F', dim=2) + G = Mapping('G', dim=2) + mapped = Domain.join([F(A), G(B)], [first, second], 'Omega') + assert tuple(mapped.interface_map) == ('A|B', 'zz') + assert mapped.logical_domain == domain + + mapped.export(filename) + restored = Domain.from_file(filename) + assert restored == mapped + for name, interface in restored.interface_map.items(): + assert str(interface.logical_domain.name) == name + + interface = domain.interface_map['zz'] + with pytest.raises(ValueError, match='interface names must be unique'): + Domain('Omega', interiors=[A.interior, B.interior], + connectivity=[interface, interface]) + +#============================================================================== +def test_explicit_interface_names_are_reserved_and_unique(): + patch_a = Square('A') + patch_b = Square('B') + domain = Domain.join( + [patch_a, patch_b], + [ + ((0, 0, -1), (1, 0, -1), +1), + ((0, 0, +1), (1, 0, +1), +1, 'A|B'), + ], + 'named') + + assert tuple(domain.interface_map) == ('A|B', 'A|B#2') + assert str(domain.interface_map['A|B'].name) == 'A|B' + assert domain.interface_map['A|B'].minus.ext == +1 + + with pytest.raises(ValueError, match='duplicate explicit interface name'): + Domain.join( + [patch_a, patch_b], + [ + ((0, 0, -1), (1, 0, -1), +1, 'seam'), + ((0, 0, +1), (1, 0, +1), +1, 'seam'), + ], + 'duplicate') + + with pytest.raises(TypeError, match='name must be a string'): + Domain.join( + [patch_a, patch_b], + [((0, 0, -1), (1, 0, -1), +1, 1)], + 'invalid-name') + + +#============================================================================== +def test_interface_requires_orientation(): + with pytest.raises(TypeError, match='minus, plus, orientation'): + Domain.join( + [Square('A'), Square('B')], + [((0, 0, 1), (1, 0, -1))], + 'AB') + + with pytest.raises(TypeError, match='minus, plus, orientation'): + Domain.join( + [Square('C'), Square('D')], + [{'minus': (0, 0, 1), + 'plus': (1, 0, -1), + 'orientation': +1}], + 'CD') + + with pytest.raises(ValueError, match='2D'): + Domain.join( + [Square('A'), Square('B')], + [((0, 0, 1), (1, 0, -1), 0)], + 'AB') + + +#============================================================================== +def test_single_patch_self_interface(): + patch = Square('A') + domain = Domain.join( + [patch], + [((0, 0, -1), (0, 0, 1), +1)], + 'periodic-A') + + interface = domain.interfaces + assert interface.minus.domain == patch.interior + assert interface.plus.domain == patch.interior + assert interface.minus.axis == interface.plus.axis == 0 + assert interface.orientation == +1 + assert len(domain.corners) == 2 + + +#============================================================================== +def test_mapped_single_patch_self_interface_uses_single_mapping(): + logical_patch = Square('A') + mapping = IdentityMapping('F', dim=2) + patch = mapping(logical_patch) + domain = Domain.join( + [patch], + [((0, 0, -1), (0, 0, 1), +1)], + 'periodic-A') + + assert domain.mapping == mapping + assert domain.logical_domain.interfaces.orientation == +1 + assert domain.interfaces.minus.domain == patch.interior + assert domain.interfaces.plus.domain == patch.interior + + #============================================================================== def test_hash(): A = Square('A', bounds1=(0, 1), bounds2=(0, 1)) diff --git a/sympde/utilities/plotting.py b/sympde/utilities/plotting.py new file mode 100644 index 00000000..a858ee33 --- /dev/null +++ b/sympde/utilities/plotting.py @@ -0,0 +1,821 @@ +#---------------------------------------------------------------------------# +# This file is part of SymPDE. # +#---------------------------------------------------------------------------# +"""Inspection and plotting helpers for symbolic multipatch topology.""" + +from dataclasses import dataclass +from itertools import product +import numbers +import sys + +import numpy as np + +__all__ = ( + 'TopologyVertexIncidence', + 'TopologyVertex', + 'collect_topology_vertices', + 'format_topology', + 'print_topology', + 'plot_domain', +) + + +@dataclass(frozen=True) +class TopologyVertexIncidence: + """One patch-local representative of a topological vertex.""" + + patch_index: int + patch_name: str + logical_corner: tuple + + +@dataclass(frozen=True) +class TopologyVertex: + """A topological vertex and all patch-local corners identified with it.""" + + index: int + incidences: tuple + is_boundary: bool + + @property + def patch_indices(self): + """Indices of incident patches, without duplicates.""" + return tuple(dict.fromkeys(i.patch_index for i in self.incidences)) + + @property + def patches(self): + """Names of incident patches, without duplicates.""" + return tuple(dict.fromkeys(i.patch_name for i in self.incidences)) + + +def _patch_name(patch): + logical_patch = patch.logical_domain + return str(logical_patch.name if logical_patch is not None else patch.name) + + +def _interfaces(domain): + return tuple(domain.interface_map.values()) + + +def _shared_vertices(domain): + return domain.shared_vertices + + +def _check_topology_dimension(domain): + if int(domain.dim) not in (2, 3): + raise ValueError( + 'topology vertex collection and annotated topology plotting ' + 'support only two- and three-dimensional domains') + + +def collect_topology_vertices(domain): + """Return all vertices of a 2D or 3D domain and their patch incidences. + + Unlike ``domain.corners``, which stores only equivalence classes created + by interface identifications, this function also includes unshared patch + corners. Each returned vertex records all of its patch-local + representatives and whether it lies on the exterior boundary. + + Parameters + ---------- + domain : sympde.topology.domain.Domain + A two- or three-dimensional logical or mapped domain. + + Returns + ------- + tuple[TopologyVertex, ...] + Topological vertices in deterministic patch/corner order. + """ + _check_topology_dimension(domain) + + ndim = int(domain.dim) + patches = domain.patches + patch_indices = {patch: index for index, patch in enumerate(patches)} + patch_names = tuple(_patch_name(patch) for patch in patches) + keys = [ + (patch_index, tuple(corner)) + for patch_index in range(len(patches)) + for corner in product((0, 1), repeat=ndim) + ] + parent = {key: key for key in keys} + + def find(key): + while parent[key] != key: + parent[key] = parent[parent[key]] + key = parent[key] + return key + + def union(first, second): + first_root = find(first) + second_root = find(second) + if first_root != second_root: + parent[second_root] = first_root + + for shared_vertex in _shared_vertices(domain): + members = [ + (patch_indices[corner.domain], tuple(map(int, corner.coordinates))) + for corner in shared_vertex.corners + ] + for member in members[1:]: + union(members[0], member) + + groups = {} + for key in keys: + groups.setdefault(find(key), []).append(key) + + glued_faces = set() + for interface in _interfaces(domain): + for face in (interface.minus, interface.plus): + glued_faces.add( + (patch_indices[face.domain], int(face.axis), int(face.ext))) + + ordered_groups = sorted( + (sorted(group) for group in groups.values()), + key=lambda group: group[0], + ) + + vertices = [] + for index, group in enumerate(ordered_groups): + incidences = tuple( + TopologyVertexIncidence( + patch_index=patch_index, + patch_name=patch_names[patch_index], + logical_corner=corner, + ) + for patch_index, corner in group + ) + local_faces = ( + ( + incidence.patch_index, + axis, + -1 if incidence.logical_corner[axis] == 0 else 1, + ) + for incidence in incidences + for axis in range(ndim) + ) + is_boundary = any(face not in glued_faces for face in local_faces) + vertices.append(TopologyVertex(index, incidences, is_boundary)) + + return tuple(vertices) + + +def format_topology(domain): + """Return a human-readable summary of patches, interfaces, and vertices.""" + _check_topology_dimension(domain) + + patches = domain.patches + lines = [ + f'=== {domain.name} ===', + 'patches: ' + ', '.join(_patch_name(patch) for patch in patches), + 'interfaces:', + ] + + interfaces = _interfaces(domain) + if not interfaces: + lines.append(' (none)') + for index, interface in enumerate(interfaces): + minus = interface.minus + plus = interface.plus + lines.extend(( + ( + f' I{index} ({interface.name}): ' + f'{_patch_name(minus.domain)}' + f'(axis={minus.axis}, ext={int(minus.ext):+d}) <-> ' + f'{_patch_name(plus.domain)}' + f'(axis={plus.axis}, ext={int(plus.ext):+d})' + ), + ( + f' orientation={interface.orientation}, ' + f'axis_map={interface.axis_map}' + ), + )) + + lines.append('vertices:') + for vertex in collect_topology_vertices(domain): + incidences = ', '.join( + f'{incidence.patch_name}{incidence.logical_corner}' + for incidence in vertex.incidences + ) + location = 'boundary' if vertex.is_boundary else 'interior' + lines.append( + f' V{vertex.index}: {location}; patches={vertex.patches}; ' + f'local corners=[{incidences}]') + + return '\n'.join(lines) + + +def print_topology(domain, file=None): + """Print :func:`format_topology` to ``file`` or standard output.""" + print(format_topology(domain), file=sys.stdout if file is None else file) + + +#============================================================================== +# Plotting style +#============================================================================== +# Okabe-Ito palette, which is safe for colour-blind readers +_PATCH_COLORS = ( + '#E69F00', '#56B4E9', '#009E73', '#F0E442', + '#0072B2', '#D55E00', '#CC79A7', '#999999', +) +_INTERFACE_COLORMAP = 'tab10' +_PLUS_LINESTYLE = (0, (3, 2)) + +_STYLE = { + # Patches + 'patch_alpha' : 0.3, # 2D patch fill + 'surface_alpha' : 0.6, # 3D patch faces, geometry plot + 'face_alpha' : 0.08, # 3D exterior faces, topology plot + 'interface_alpha' : 0.3, # 3D minus faces of the interfaces + 'boundary_width' : 1.8, # exterior boundary + 'edge_width' : 0.8, # other patch edges + 'edge_color' : '0.45', + 'isoline_width' : 0.5, + 'isoline_color' : '0.55', + 'isoline_count' : 7, # interior isolines per direction + # Interface sides + 'side_width' : 2.2, + 'side_inset' : 0.06, # fraction of the patch width, in 2D + 'arrow_size' : 14, + 'connector_width' : 1.2, + # Labels and vertices + 'patch_fontsize' : 11, + 'label_fontsize' : 9, + 'vertex_fontsize' : 7, + 'vertex_size' : 5, + 'halo_width' : 2.5, + # Figures + 'figsize_2d' : (7.5, 6.5), + 'figsize_3d' : (8.0, 7.0), + 'view' : (24, -58), # elevation and azimuth in 3D + 'min_box_aspect' : 0.3, # smallest 3D box side, relative +} + + +def _halo(): + import matplotlib.patheffects as path_effects + return [path_effects.withStroke( + linewidth=_STYLE['halo_width'], foreground='white')] + + +#============================================================================== +# Sampling of the geometry +#============================================================================== +def _logical_limits(patch): + logical_patch = patch.logical_domain + parameter_domain = logical_patch if logical_patch is not None else patch + return ( + np.asarray(parameter_domain.min_coords, dtype=float), + np.asarray(parameter_domain.max_coords, dtype=float), + ) + + +def _map_points(patch, logical_points): + """Map an array of logical points, of shape (n, dim), to the physical + space.""" + logical_points = np.asarray(logical_points, dtype=float) + if patch.mapping is None: + return logical_points + mapping = patch.mapping.get_callable_mapping() + components = mapping(*logical_points.T) + return np.column_stack([ + np.broadcast_to(np.asarray(c, dtype=float), len(logical_points)) + for c in components]) + + +def _physical_dimension(domain): + dimensions = { + int(patch.mapping.pdim) if patch.mapping is not None + else int(patch.dim) + for patch in domain.patches + } + if len(dimensions) != 1: + raise ValueError('all plotted patches must have one physical dimension') + return dimensions.pop() + + +def _face_points(patch, axis, ext, samples, inset=0.0): + """Logical points on a face of a patch, of shape (samples, ..., dim). + + The face is moved into the patch by the fraction `inset` of the patch + width along its normal axis. + """ + lower, upper = _logical_limits(patch) + dim = len(lower) + tangents = [a for a in range(dim) if a != axis] + grids = np.meshgrid( + *[np.linspace(lower[a], upper[a], samples) for a in tangents], + indexing='ij') + width = upper[axis] - lower[axis] + position = lower[axis] + inset * width if ext == -1 else \ + upper[axis] - inset * width + points = np.empty((*grids[0].shape, dim)) + points[..., axis] = position + for a, grid in zip(tangents, grids): + points[..., a] = grid + return points + + +def _map_face(patch, axis, ext, samples, inset=0.0): + points = _face_points(patch, axis, ext, samples, inset) + shape = points.shape + return _map_points(patch, points.reshape(-1, shape[-1])).reshape(shape) + + +def _patch_center(patch): + lower, upper = _logical_limits(patch) + return _map_points(patch, [(lower + upper) / 2])[0] + + +def _face_center(face): + lower, upper = _logical_limits(face.domain) + point = (lower + upper) / 2 + point[face.axis] = lower[face.axis] if face.ext == -1 else upper[face.axis] + return _map_points(face.domain, [point])[0] + + +def _corner_point(patch, logical_corner): + lower, upper = _logical_limits(patch) + point = np.where(np.asarray(logical_corner) == 0, lower, upper) + return _map_points(patch, [point])[0] + + +def _patch_polygon(patch, samples): + """Closed boundary of a 2D patch, counter-clockwise in logical space.""" + return np.vstack(( + _map_face(patch, 1, -1, samples), + _map_face(patch, 0, +1, samples), + _map_face(patch, 1, +1, samples)[::-1], + _map_face(patch, 0, -1, samples)[::-1], + )) + + +#============================================================================== +# Plot model: what is drawn, independently of Matplotlib +#============================================================================== +@dataclass(frozen=True) +class _PlotInterface: + """An interface with the points needed to draw it.""" + + index: int + interface: object + color: tuple + minus_center: np.ndarray + plus_center: np.ndarray + separated: bool # the two sides are drawn at different places + + +@dataclass(frozen=True) +class _PlotModel: + """The sampled geometry and topology of a domain.""" + + domain: object + dim: int + refinement: int + patches: tuple + patch_colors: tuple + faces: tuple # per patch, {(axis, ext): sampled face} + exterior_faces: frozenset # (patch index, axis, ext) + interfaces: tuple # _PlotInterface + vertices: tuple # (TopologyVertex, tuple of physical points) + extent: np.ndarray # size of the bounding box + + @property + def scale(self): + return max(float(np.max(self.extent)), 1.0) + + +def _patch_neighbours(domain, patches): + """Pairs of distinct patches that share an interface or a vertex.""" + index = {patch: i for i, patch in enumerate(patches)} + neighbours = {i: set() for i in range(len(patches))} + + def connect(group): + for i in group: + neighbours[i].update(j for j in group if j != i) + + for interface in _interfaces(domain): + connect({index[interface.minus.domain], index[interface.plus.domain]}) + for vertex in collect_topology_vertices(domain): + connect(set(vertex.patch_indices)) + return neighbours + + +def _color_patches(domain, patches): + """Give every patch its own colour if the palette is large enough, and + otherwise give neighbouring patches different colours, when possible.""" + if len(patches) <= len(_PATCH_COLORS): + return _PATCH_COLORS[:len(patches)] + neighbours = _patch_neighbours(domain, patches) + colors = {} + for i in range(len(patches)): + used = {colors[j] for j in neighbours[i] if j in colors} + free = [c for c in range(len(_PATCH_COLORS)) if c not in used] + colors[i] = free[0] if free else i % len(_PATCH_COLORS) + return tuple(_PATCH_COLORS[colors[i]] for i in range(len(patches))) + + +def _build_model(domain, refinement, topology): + import matplotlib.pyplot as plt + + dim = int(domain.dim) + patches = domain.patches + faces = tuple( + {(axis, ext): _map_face(patch, axis, ext, refinement) + for axis in range(dim) for ext in (-1, 1)} + for patch in patches) + points = np.vstack([ + face.reshape(-1, dim) for patch_faces in faces + for face in patch_faces.values()]) + extent = np.ptp(points, axis=0) + + patch_index = {patch: i for i, patch in enumerate(patches)} + glued_faces = { + (patch_index[face.domain], int(face.axis), int(face.ext)) + for interface in _interfaces(domain) + for face in (interface.minus, interface.plus)} + exterior_faces = frozenset( + (i, axis, ext) for i in range(len(patches)) + for axis in range(dim) for ext in (-1, 1) + if (i, axis, ext) not in glued_faces) + + if not topology: + return _PlotModel( + domain, dim, refinement, patches, + _color_patches(domain, patches), + faces, exterior_faces, (), (), extent) + + scale = max(float(np.max(extent)), 1.0) + interface_colors = plt.get_cmap(_INTERFACE_COLORMAP).colors + interfaces = [] + for index, interface in enumerate(_interfaces(domain)): + minus_center = _face_center(interface.minus) + plus_center = _face_center(interface.plus) + separated = ( + np.linalg.norm(minus_center - plus_center) > 1e-10 * scale) + interfaces.append(_PlotInterface( + index, interface, + interface_colors[index % len(interface_colors)], + minus_center, plus_center, separated)) + + vertices = [] + for vertex in collect_topology_vertices(domain): + locations = {} + for incidence in vertex.incidences: + point = _corner_point( + patches[incidence.patch_index], incidence.logical_corner) + locations.setdefault(tuple(np.round(point, decimals=12)), point) + vertices.append((vertex, tuple(locations.values()))) + + return _PlotModel( + domain, dim, refinement, patches, _color_patches(domain, patches), + faces, exterior_faces, tuple(interfaces), tuple(vertices), extent) + + +#============================================================================== +# Drawing helpers shared by 2D and 3D +#============================================================================== +def _new_axes(dim, ax, topology): + import matplotlib.pyplot as plt + + if ax is None: + figsize = _STYLE['figsize_2d'] if dim == 2 else _STYLE['figsize_3d'] + figure = plt.figure( + figsize=figsize if topology else None, layout='constrained') + ax = figure.add_subplot(projection='3d' if dim == 3 else None) + return figure, ax + + if hasattr(ax, 'zaxis') != (dim == 3): + raise ValueError( + f'a {dim}D Matplotlib axes is required for a {dim}D domain') + return ax.figure, ax + + +def _text(ax, point, text, **options): + """Write a label with a white halo, in 2D or 3D.""" + options.setdefault('ha', 'center') + options.setdefault('va', 'center') + ax.text(*point, text, path_effects=_halo(), **options) + + +def _draw_interface_labels(ax, model): + for item in model.interfaces: + style = dict( + color=item.color, fontsize=_STYLE['label_fontsize'], + fontweight='bold', zorder=7, + bbox=dict(boxstyle='round,pad=0.15', facecolor='white', + edgecolor=item.color, linewidth=0.6)) + if not item.separated: + _text(ax, item.minus_center, f'I{item.index}', **style) + continue + _text(ax, item.minus_center, f'I{item.index}-', **style) + _text(ax, item.plus_center, f'I{item.index}+', **style) + + +def _draw_vertices(ax, model, labels): + for vertex, points in model.vertices: + face_color = 'white' if vertex.is_boundary else 'black' + for point in points: + ax.plot(*([coordinate] for coordinate in point), marker='o', + linestyle='', markersize=_STYLE['vertex_size'], + markerfacecolor=face_color, markeredgecolor='black', + zorder=8) + if labels: + _text(ax, point, f'V{vertex.index}', ha='left', va='bottom', + fontsize=_STYLE['vertex_fontsize'], color='0.2', + zorder=9) + + +def _draw_legend(ax, dim): + from matplotlib.lines import Line2D + + side = 'side' if dim == 2 else 'face arrows' + handles = [ + Line2D([], [], color='0.3', linewidth=_STYLE['side_width'], + label=f'minus {side}'), + Line2D([], [], color='0.3', linewidth=_STYLE['side_width'], + linestyle=_PLUS_LINESTYLE, label=f'plus {side}'), + ] + if dim == 2: + handles.append(Line2D( + [], [], color='black', linewidth=_STYLE['boundary_width'], + label='exterior boundary')) + handles += [ + Line2D([], [], marker='o', linestyle='', markerfacecolor='white', + markeredgecolor='black', label='boundary vertex'), + Line2D([], [], marker='o', linestyle='', markerfacecolor='black', + markeredgecolor='black', label='interior vertex'), + ] + ax.legend(handles=handles, loc='upper center', bbox_to_anchor=(0.5, 0.0), + ncol=3, fontsize=8, frameon=False) + + +#============================================================================== +# 2D plots +#============================================================================== +def _draw_patches_2d(ax, model, isolines, topology, labels): + for i, patch in enumerate(model.patches): + polygon = _patch_polygon(patch, model.refinement) + ax.fill(*polygon.T, facecolor=model.patch_colors[i], + alpha=_STYLE['patch_alpha'], linewidth=0, zorder=1) + if isolines: + _draw_isolines_2d(ax, patch, model.refinement) + for (axis, ext), curve in model.faces[i].items(): + exterior = (i, axis, ext) in model.exterior_faces + ax.plot(*curve.T, zorder=3, **( + dict(color='black', linewidth=_STYLE['boundary_width']) + if exterior or not topology else + dict(color=_STYLE['edge_color'], + linewidth=_STYLE['edge_width']))) + if labels: + _text(ax, _patch_center(patch), _patch_name(patch), + fontsize=_STYLE['patch_fontsize'], fontweight='bold', + zorder=6) + + +def _draw_isolines_2d(ax, patch, samples): + lower, upper = _logical_limits(patch) + fractions = np.linspace(0, 1, _STYLE['isoline_count'] + 2)[1:-1] + for axis in range(2): + tangent = 1 - axis + points = np.empty((samples, 2)) + points[:, tangent] = np.linspace(lower[tangent], upper[tangent], samples) + for fraction in fractions: + points[:, axis] = lower[axis] + fraction * (upper[axis] - lower[axis]) + ax.plot(*_map_points(patch, points).T, + color=_STYLE['isoline_color'], + linewidth=_STYLE['isoline_width'], zorder=2) + + +def _draw_side_2d(ax, face, color, linestyle, samples, direction=1): + """Draw an interface side inside its patch, with an arrow along its + tangential coordinate, reversed if `direction` is -1.""" + from matplotlib.patches import FancyArrowPatch + + curve = _map_face(face.domain, int(face.axis), int(face.ext), samples, + inset=_STYLE['side_inset']) + ax.plot(*curve.T, color=color, linewidth=_STYLE['side_width'], + linestyle=linestyle, solid_capstyle='butt', zorder=4) + start, end = curve[int(0.45 * samples)], curve[int(0.62 * samples)] + if direction == -1: + start, end = end, start + ax.add_patch(FancyArrowPatch( + start, end, arrowstyle='-|>', mutation_scale=_STYLE['arrow_size'], + color=color, linewidth=0, zorder=5)) + + +def _draw_interfaces_2d(ax, model): + from matplotlib.patches import FancyArrowPatch + + for item in model.interfaces: + # The arrow of the plus side shows where the coordinate of the minus + # side goes, so that it is reversed for a reversing interface + (_, _, direction), = item.interface.axis_map + _draw_side_2d(ax, item.interface.minus, item.color, '-', + model.refinement) + _draw_side_2d(ax, item.interface.plus, item.color, _PLUS_LINESTYLE, + model.refinement, direction) + if item.separated: + ax.add_patch(FancyArrowPatch( + item.minus_center, item.plus_center, arrowstyle='-', + connectionstyle='arc3,rad=0.25', linestyle=':', + linewidth=_STYLE['connector_width'], color=item.color, + zorder=2)) + + +def _plot_domain_2d(model, ax, isolines, topology, patch_labels, + interface_labels, vertex_labels, legend): + _draw_patches_2d(ax, model, isolines, topology, topology and patch_labels) + if topology: + _draw_interfaces_2d(ax, model) + if interface_labels: + _draw_interface_labels(ax, model) + _draw_vertices(ax, model, vertex_labels) + if legend: + _draw_legend(ax, 2) + ax.set_axis_off() + ax.margins(0.03) + else: + ax.set_xlabel('x') + ax.set_ylabel('y') + ax.set_aspect('equal', adjustable='datalim') + + +#============================================================================== +# 3D plots +#============================================================================== +def _surface_edges(surface): + return (surface[0], surface[-1], surface[:, 0], surface[:, -1]) + + +def _draw_patches_3d(ax, model, isolines, topology, labels): + for i, patch in enumerate(model.patches): + color = model.patch_colors[i] + for (axis, ext), surface in model.faces[i].items(): + if not topology: + ax.plot_surface(*np.moveaxis(surface, -1, 0), color=color, + alpha=_STYLE['surface_alpha'], linewidth=0, + shade=True) + elif (i, axis, ext) in model.exterior_faces: + ax.plot_surface(*np.moveaxis(surface, -1, 0), color=color, + alpha=_STYLE['face_alpha'], linewidth=0, + shade=False) + for edge in _surface_edges(surface): + ax.plot(*edge.T, color=_STYLE['edge_color'], + linewidth=_STYLE['edge_width']) + if isolines: + stride = max(1, model.refinement // 5) + ax.plot_wireframe(*np.moveaxis(surface, -1, 0), + rstride=stride, cstride=stride, + color=_STYLE['isoline_color'], + linewidth=_STYLE['isoline_width'], + alpha=0.5) + if labels: + _text(ax, _patch_center(patch), _patch_name(patch), + fontsize=_STYLE['patch_fontsize'], fontweight='bold') + + +def _draw_side_arrows_3d(ax, face, color, linestyle, alpha, arrows): + """Draw arrows on a face along its tangential axes. + + Each arrow is given by `(axis, direction, label)`: it follows the logical + `axis` of the face, reversed if `direction` is -1, and is labelled with + the coordinate eta of index `label`. + """ + lower, upper = _logical_limits(face.domain) + outward = 0.10 * (_face_center(face) - _patch_center(face.domain)) + for tangent, direction, label in arrows: + start = (lower + upper) / 2 + start[face.axis] = (lower if face.ext == -1 else upper)[face.axis] + end = start.copy() + width = upper[tangent] - lower[tangent] + start[tangent] = lower[tangent] + 0.25 * width + end[tangent] = lower[tangent] + 0.72 * width + start, end = _map_points(face.domain, [start, end]) + if direction == -1: + start, end = end, start + delta = end - start + ax.quiver(*start, *delta, color=color, linestyles=linestyle, + linewidth=1.5, arrow_length_ratio=0.18, alpha=alpha) + _text(ax, end + 0.08 * delta + outward, rf'$\eta_{label + 1}$', + color=color, fontsize=_STYLE['label_fontsize']) + + +def _draw_interfaces_3d(ax, model): + for item in model.interfaces: + interface = item.interface + surface = model.faces[model.patches.index(interface.minus.domain)][ + (int(interface.minus.axis), int(interface.minus.ext))] + ax.plot_surface(*np.moveaxis(surface, -1, 0), color=item.color, + alpha=_STYLE['interface_alpha'], linewidth=0, + shade=False) + # The arrows of the plus face show where the coordinates of the + # minus face go, with the labels of the minus coordinates + axis_map = interface.axis_map + _draw_side_arrows_3d( + ax, interface.minus, item.color, '-', 1.0, + [(minus_axis, 1, minus_axis) for minus_axis, _, _ in axis_map]) + _draw_side_arrows_3d( + ax, interface.plus, item.color, '--', 0.7, + [(plus_axis, direction, minus_axis) + for minus_axis, plus_axis, direction in axis_map]) + if item.separated: + ax.plot(*np.column_stack((item.minus_center, item.plus_center)), + color=item.color, linestyle=':', + linewidth=_STYLE['connector_width']) + + +def _plot_domain_3d(model, ax, isolines, topology, patch_labels, + interface_labels, vertex_labels, legend): + _draw_patches_3d(ax, model, isolines, topology, topology and patch_labels) + if topology: + _draw_interfaces_3d(ax, model) + if interface_labels: + _draw_interface_labels(ax, model) + _draw_vertices(ax, model, vertex_labels) + if legend: + _draw_legend(ax, 3) + ax.view_init(*_STYLE['view']) + ax.set_xlabel('x') + ax.set_ylabel('y') + ax.set_zlabel('z') + extent = np.maximum(model.extent, 1e-12) + ax.set_box_aspect(np.maximum(extent, + _STYLE['min_box_aspect'] * extent.max())) + + +#============================================================================== +# Entry point +#============================================================================== +def plot_domain( + domain, draw=True, refinement=None, isolines=False, ax=None, + topology=False, + patch_labels=True, interface_labels=True, vertex_labels=True, + legend=True): + """Plot a 2D or 3D domain, optionally with topology annotations. + + This is the common plotting entry point for ordinary geometry plots and + annotated multipatch topology plots. Both modes use the same sampling of + the mapped patch edges and faces, and colour neighbouring patches + differently. + + With ``topology=True``, each side of an interface is drawn inside its own + patch, the minus side solid and the plus side dashed, in a colour shared + by both sides. The arrows of the minus side show its increasing + tangential coordinates, and those of the plus side where these go on the + plus side: they point the other way for a reversing interface. In 3D, the + arrows are labelled with the coordinates of the minus side. + A dotted connector joins the two sides of an interface when they are drawn + at different places, as for periodic or self-interfaces. The exterior + boundary is drawn in black, and the vertices are marked: white for + boundary vertices, black for interior ones. :func:`print_topology` gives + the same information as text, including the axis map of every interface. + + Parameters + ---------- + domain : sympde.topology.domain.Domain + Two- or three-dimensional logical or mapped domain with equal logical + and physical dimensions. + draw : bool, default=True + Call ``matplotlib.pyplot.show`` after constructing the figure. + refinement : int or None + Samples per edge. Defaults to 40 in 2D and 15 in 3D. + isolines : bool, default=False + Draw logical isolines or face wireframes. + ax : matplotlib.axes.Axes, optional + Existing 2D or 3D axes on which to draw. + topology : bool, default=False + Add symbolic multipatch topology annotations. + patch_labels, interface_labels, vertex_labels, legend : bool + Control annotations when ``topology=True``. They are ignored in the + ordinary geometry mode. + + Returns + ------- + matplotlib.figure.Figure + The containing figure. + """ + import matplotlib.pyplot as plt + + dim = int(domain.dim) + if dim not in (2, 3): + raise ValueError('domain plotting supports only dimensions 2 and 3') + if _physical_dimension(domain) != dim: + raise ValueError( + 'domain plotting currently requires equal logical and physical ' + 'dimensions') + if refinement is None: + refinement = 40 if dim == 2 else 15 + if isinstance(refinement, bool) or not isinstance( + refinement, numbers.Integral) or refinement < 2: + raise ValueError('refinement must be an integer greater than one') + + model = _build_model(domain, int(refinement), topology) + figure, ax = _new_axes(dim, ax, topology) + plot = _plot_domain_2d if dim == 2 else _plot_domain_3d + plot(model, ax, isolines, topology, patch_labels, interface_labels, + vertex_labels, legend) + ax.set_title(str(domain.name).replace('_', ' ')) + + if draw: + plt.show() + return figure diff --git a/sympde/utilities/tests/test_utilities.py b/sympde/utilities/tests/test_utilities.py index 2b2e8805..9eec36d4 100644 --- a/sympde/utilities/tests/test_utilities.py +++ b/sympde/utilities/tests/test_utilities.py @@ -1,13 +1,21 @@ -import pytest +from io import StringIO + +import matplotlib.pyplot as plt import numpy as np +import pytest from sympy import Matrix, symbols, Array from sympy import S +import sympde.utilities.utils as utilities_utils +from sympde.topology.analytical_mapping import AffineMapping, IdentityMapping +from sympde.topology.domain import Cube, Domain, Square +from sympde.utilities.plotting import ( + collect_topology_vertices, + format_topology, + plot_domain, + print_topology, +) from sympde.utilities.utils import lambdify_sympde -from sympde.utilities.utils import plot_domain -from sympde.topology.domain import Square - - def test_lambdify_sympde_1d(): @@ -157,3 +165,174 @@ def test_lambdify_sympde_2d(): assert np.array_equal(sparse_input_ia, expected_ia) assert np.array_equal(sparse_input_da, expected_da) assert np.array_equal(sparse_input_sa, expected_sa) + + +def _make_twisted_strip(): + patch_a = IdentityMapping('F_A', dim=2)(Square('A')) + patch_b = AffineMapping( + 'F_B', dim=2, + c1=1, c2=0, + a11=1, a12=0, + a21=0, a22=1, + )(Square('B')) + return Domain.join( + [patch_a, patch_b], + [ + ((0, 0, 1), (1, 0, -1), +1), + ((0, 0, -1), (1, 0, 1), -1), + ], + 'twisted_strip', + ) + + +def _make_three_patch_vertex(): + patches = [Square(name) for name in ('A', 'B', 'C')] + return Domain.join( + patches, + [ + ((0, 0, -1), (1, 1, -1), +1), + ((1, 0, -1), (2, 1, -1), +1), + ((2, 0, -1), (0, 1, -1), +1), + ], + 'three_patch_vertex', + ) + + +def _make_oriented_cubes(): + patch_a = IdentityMapping('F_A3', dim=3)(Cube('A3')) + patch_b = AffineMapping( + 'F_B3', dim=3, + c1=1, c2=0, c3=1, + a11=0, a12=1, a13=0, + a21=0, a22=0, a23=1, + a31=-1, a32=0, a33=0, + )(Cube('B3')) + return Domain.join( + [patch_a, patch_b], + [((0, 0, 1), (1, 1, -1), (-1, +1, -1))], + 'oriented_cubes', + ) + + +def test_collect_topology_vertices_shared_and_self_interfaces(): + vertices = collect_topology_vertices(_make_three_patch_vertex()) + interior_vertices = [vertex for vertex in vertices if not vertex.is_boundary] + + assert len(interior_vertices) == 1 + center = interior_vertices[0] + assert center.patch_indices == (0, 1, 2) + assert center.patches == ('A', 'B', 'C') + assert all( + incidence.logical_corner == (0, 0) + for incidence in center.incidences) + + patch = Square('P') + self_connected = Domain.join( + [patch], + [((0, 0, -1), (0, 0, 1), +1)], + 'self_connected', + ) + self_vertices = collect_topology_vertices(self_connected) + assert len(self_vertices) == 2 + assert all(len(vertex.incidences) == 2 for vertex in self_vertices) + assert all(vertex.patches == ('P',) for vertex in self_vertices) + assert all(vertex.is_boundary for vertex in self_vertices) + + +def test_collect_and_format_3d_topology(): + domain = _make_oriented_cubes() + vertices = collect_topology_vertices(domain) + shared_vertices = [ + vertex for vertex in vertices if len(vertex.incidences) == 2] + + assert len(vertices) == 12 + assert len(shared_vertices) == 4 + assert all(vertex.patches == ('A3', 'B3') for vertex in shared_vertices) + + text = format_topology(domain) + assert 'patches: A3, B3' in text + assert 'orientation=(-1, 1, -1)' in text + assert 'axis_map=((1, 2, 1), (2, 0, -1))' in text + + +def test_format_and_print_topology(): + domain = _make_twisted_strip() + text = format_topology(domain) + + assert 'patches: A, B' in text + assert 'I0' in text and 'I1' in text + assert 'orientation=1' in text + assert 'orientation=-1' in text + assert 'axis_map=((1, 1, -1),)' in text + assert 'vertices:' in text + + stream = StringIO() + print_topology(domain, file=stream) + assert stream.getvalue() == text + '\n' + + +def test_plot_domain_2d_integration(): + assert not hasattr(utilities_utils, 'plot_domain') + + domain = _make_twisted_strip() + figure = plot_domain( + domain, draw=False, refinement=8, isolines=True, topology=True) + labels = {text.get_text() for text in figure.axes[0].texts} + # The sides of I0 are drawn at different places, those of I1 coincide + assert {'A', 'B', 'I0-', 'I0+', 'I1'}.issubset(labels) + assert any(label.startswith('V') for label in labels) + assert figure.axes[0].get_title() == 'twisted strip' + plt.close(figure) + + figure = plot_domain(Square('single'), draw=False, refinement=8) + assert len(figure.axes[0].lines) == 4 + plt.close(figure) + + +def test_plot_domain_3d_integration(): + domain = _make_oriented_cubes() + figure = plot_domain( + domain, draw=False, refinement=4, isolines=True, topology=True) + axis = figure.axes[0] + labels = {text.get_text() for text in axis.texts} + assert {'A3', 'B3', 'I0'}.issubset(labels) + assert any(label.startswith('V') for label in labels) + assert axis.get_zlabel() == 'z' + assert figure.axes[0].get_title() == 'oriented cubes' + plt.close(figure) + + figure = plot_domain(Cube('single3d'), draw=False, refinement=4) + axis = figure.axes[0] + assert axis.get_zlabel() == 'z' + assert len(axis.collections) == 6 + plt.close(figure) + + +def test_plot_domain_on_existing_axes_with_numpy_refinement(): + figure, (axis_2d, axis_3d) = plt.subplots(1, 2) + + returned = plot_domain( + _make_three_patch_vertex(), draw=False, refinement=np.int64(6), + ax=axis_2d, topology=True) + assert returned is figure + + with pytest.raises(ValueError, match='3D Matplotlib axes'): + plot_domain(Cube('C'), draw=False, ax=axis_3d) + with pytest.raises(ValueError, match='refinement'): + plot_domain(Square('S'), draw=False, refinement=True, ax=axis_2d) + plt.close(figure) + + +def test_neighbouring_patches_have_different_colors(): + from sympde.topology.multipatch_gallery import build_multipatch_domain_2d + from sympde.utilities.plotting import _color_patches, _patch_neighbours + + # More patches than colours in the palette + domain = build_multipatch_domain_2d('pretzel_f') + patches = domain.patches + colors = _color_patches(domain, patches) + neighbours = _patch_neighbours(domain, patches) + + assert len(patches) > len(set(colors)) + assert all(colors[i] != colors[j] + for i in neighbours for j in neighbours[i]) diff --git a/sympde/utilities/utils.py b/sympde/utilities/utils.py index 5236fffa..858fd69a 100644 --- a/sympde/utilities/utils.py +++ b/sympde/utilities/utils.py @@ -1,10 +1,8 @@ -import numpy as np import itertools as it + +import numpy as np from sympy import lambdify -from sympde.topology.analytical_mapping import IdentityMapping -from sympde.topology.basic import InteriorDomain -from sympde.topology.mapping import MultiPatchMapping def lambdify_sympde(variables, expr): """ @@ -58,6 +56,7 @@ def lambdify_sympde(variables, expr): scalar_shape = array_expr.shape if scalar_shape == (): f = lambdify(variables, expr, 'numpy') + def f_vec_sc(*XYZ): b = np.broadcast(*XYZ) if b.ndim == 0: @@ -71,232 +70,15 @@ def f_vec_sc(*XYZ): return result return f_vec_sc - else: - scalar_functions = {} - for multi_index in it.product(*tuple(range(s) for s in scalar_shape)): - scalar_functions[multi_index] = lambdify(variables, array_expr[multi_index], 'numpy') - - def f_vec_v(*XYZ): - b = np.broadcast(*XYZ) - result = np.zeros(scalar_shape + b.shape) - for multi_index in it.product(*tuple(range(s) for s in scalar_shape)): - result[multi_index] = scalar_functions[multi_index](*XYZ) - return result - return f_vec_v - - -def plot_domain(domain, draw=True, isolines=False, refinement=None): - """ - Plots a 2D or 3D domain using matplotlib - - Parameters - ---------- - domain : sympde.topology.domain.Domain - Domain to plot - - draw : bool, default=True - If true, plt.show() will be called. - - isolines : bool, default=False - If true and the domain is 2D, also plots iso-lines. - - refinement : int or None - Number of straight line segments used to approximate each boundary edge. - If None, uses 15 for 3D domains and 40 for 2D domains - """ - pdim = domain.dim if domain.mapping is None else domain.mapping.pdim - if pdim == 2: - if refinement is None: - plot_2d(domain, draw=draw, isolines=isolines) - else: - plot_2d(domain, draw=draw, isolines=isolines, refinement=refinement) - elif pdim ==3: - if refinement is None: - plot_3d(domain, draw=draw) - else: - plot_3d(domain, draw=draw, refinement=refinement) - - -def plot_2d(domain, draw=True, isolines=False, refinement=40): - """ - Plot a 2D domain - - Parameters - ---------- - domain : sympde.topology.domain.Domain - Domain to plot - - draw : bool - if true, plt.show() will be called. - - refinement : int - Number of straight line segments used to approximate each boundary edge. - """ - import matplotlib.pyplot as plt - - fig = plt.figure() - ax = fig.add_subplot(111) - - if isinstance(domain.interior, InteriorDomain): - plot_2d_single_patch(domain.interior, domain.mapping, ax, isolines=isolines, refinement=refinement) - else: - if isinstance(domain.mapping, MultiPatchMapping): - for patch, mapping in domain.mapping.mappings.items(): - plot_2d_single_patch(patch, mapping, ax, isolines=isolines, refinement=refinement) - else: - for interior in domain.interior.as_tuple(): - plot_2d_single_patch(interior, interior.mapping, ax, isolines=isolines, refinement=refinement) - - ax.set_aspect('equal', adjustable='box') - ax.set_xlabel('X') - ax.set_ylabel('Y', rotation='horizontal') - if draw: - plt.show() - -def plot_3d(domain, draw=True, refinement=15): - """ - Plot a 3D domain - - Parameters - ---------- - domain : sympde.topology.domain.Domain - Domain to plot - - draw : bool - if true, plt.show() will be called. - - refinement : int - Number of straight line segments used to approximate each boundary edge. - """ - import matplotlib.pyplot as plt - - mapping = domain.mapping - - fig = plt.figure() - ax = fig.add_subplot(111, projection="3d") + scalar_functions = {} + for multi_index in it.product(*tuple(range(s) for s in scalar_shape)): + scalar_functions[multi_index] = lambdify( + variables, array_expr[multi_index], 'numpy') - if isinstance(domain.interior, InteriorDomain): - plot_3d_single_patch(domain.interior, domain.mapping, ax, refinement=refinement) - else: - if isinstance(domain.mapping, MultiPatchMapping): - for patch, mapping in domain.mapping.mappings.items(): - plot_3d_single_patch(patch, mapping, ax, refinement=refinement) - else: - for interior in domain.interior.as_tuple(): - plot_3d_single_patch(interior, interior.mapping, ax, refinement=refinement) - - ax.set_xlabel('X') - ax.set_ylabel('Y', rotation='horizontal') - ax.set_zlabel('Z') - if draw: - plt.show() - -def plot_3d_single_patch(patch, mapping, ax, refinement=15): - """ - Plot a singe patch in a 3D domain - - Parameters - ---------- - patch : sympde.topology.basic.InteriorDomain - - mapping : sympde.topology.mapping - - ax : mpl_toolkits.mplot3d.axes3d.Axes3D - Axes object on which the patch is drawn. - - refinement : int, default=15 - Number of straight line segments used to approximate each boundary edge. - """ - if mapping is None: - mapping = IdentityMapping('Id', dim=3) - - map_call = mapping.get_callable_mapping() - refinement += 1 - - linspace_0 = np.linspace(patch.min_coords[0], patch.max_coords[0], refinement, endpoint=True) - linspace_1 = np.linspace(patch.min_coords[1], patch.max_coords[1], refinement, endpoint=True) - linspace_2 = np.linspace(patch.min_coords[2], patch.max_coords[2], refinement, endpoint=True) - - grid_01 = np.meshgrid(linspace_0, linspace_1, indexing='ij', sparse=True) - grid_02 = np.meshgrid(linspace_0, linspace_2, indexing='ij', sparse=True) - grid_12 = np.meshgrid(linspace_1, linspace_2, indexing='ij', sparse=True) - - full_00 = np.full((refinement, refinement), linspace_0[0]) - full_01 = np.full((refinement, refinement), linspace_0[-1]) - full_10 = np.full((refinement, refinement), linspace_1[0]) - full_11 = np.full((refinement, refinement), linspace_1[-1]) - full_20 = np.full((refinement, refinement), linspace_2[0]) - full_21 = np.full((refinement, refinement), linspace_2[-1]) - - mesh_01_0 = map_call(*grid_01, full_20) - mesh_01_1 = map_call(*grid_01, full_21) - - mesh_02_0 = map_call(grid_02[0], full_10, grid_02[1]) - mesh_02_1 = map_call(grid_02[0], full_11, grid_02[1]) - - mesh_12_0 = map_call(full_00, *grid_12) - mesh_12_1 = map_call(full_01, *grid_12) - - kwargs_plot = {'color': 'c', 'alpha': 0.7} - - ax.plot_surface(*mesh_01_0, **kwargs_plot) - ax.plot_surface(*mesh_01_1, **kwargs_plot) - ax.plot_surface(*mesh_02_0, **kwargs_plot) - ax.plot_surface(*mesh_02_1, **kwargs_plot) - ax.plot_surface(*mesh_12_0, **kwargs_plot) - ax.plot_surface(*mesh_12_1, **kwargs_plot) - - -def plot_2d_single_patch(patch, mapping, ax, isolines=False, refinement=40): - """ - Plots a singe patch in a 2D domain - - Parameters - ---------- - patch : sympde.topology.basic.InteriorDomain - - mapping : sympde.topology.mapping - - ax : matplotlib.axes.Axes - Axes object on which the patch is drawn. - - isolines : bool, default=False - If true also plots some iso-lines - - refinement : int, default=40 - Number of straight line segments used to approximate each boundary edge. - """ - if mapping is None: - mapping = IdentityMapping('Id', dim=2) - - refinement+=1 - map_call = mapping.get_callable_mapping() - linspace_0 = np.linspace(patch.min_coords[0], patch.max_coords[0], refinement, endpoint=True) - linspace_1 = np.linspace(patch.min_coords[1], patch.max_coords[1], refinement, endpoint=True) - - if isolines: - mesh_grid = np.meshgrid(linspace_0, linspace_1, indexing='ij') - - XX, YY = map_call(*mesh_grid) - - ax.plot(XX[:, ::5], YY[:, ::5], color='darkgrey') - ax.plot(XX[::5, :].T, YY[::5, :].T, color='darkgrey') - - X_00, Y_00 = map_call(linspace_0, np.full(refinement, linspace_1[0])) - X_01, Y_01 = map_call(linspace_0, np.full(refinement, linspace_1[-1])) - X_10, Y_10 = map_call(np.full(refinement, linspace_0[0]), linspace_1) - X_11, Y_11 = map_call(np.full(refinement, linspace_0[-1]), linspace_1) - - ax.plot(X_00, Y_00, 'k') - ax.plot(X_01, Y_01, 'k') - ax.plot(X_10, Y_10, 'k') - ax.plot(X_11, Y_11, 'k') - -if __name__ == '__main__': - from sympde.topology.analytical_mapping import PolarMapping - from sympde.topology.domain import Square - A = Square('A', bounds1=(0, 1), bounds2=(0, np.pi/2)) - F = PolarMapping('F', c1=0, c2=0, rmin=0.5, rmax=1) - Omega = F(A) - - plot_domain(Omega, draw=True, isolines=True) + def f_vec_v(*XYZ): + b = np.broadcast(*XYZ) + result = np.zeros(scalar_shape + b.shape) + for multi_index in it.product(*tuple(range(s) for s in scalar_shape)): + result[multi_index] = scalar_functions[multi_index](*XYZ) + return result + return f_vec_v