Add Interface object for multipatch.

Add Interface object for multipatch geometry and coordinate transformation on the edge from one patch to another.

(From @pvidal and @ahoffmann )

See merge request gysela-developpers/gyselalibxx!485

--------------------------------------------

Co-authored-by: Pauline Vidal <[email protected]>
Co-authored-by: Emily Bourne <[email protected]>
Co-authored-by: alex-m-h <[email protected]>

src/CMakeLists.txt

@@ -11,6 +11,7 @@ add_subdirectory(geometryXVx)
 add_subdirectory(geometryRTheta)
 add_subdirectory(geometryXYVxVy)
 add_subdirectory(interpolation)
+add_subdirectory(multipatch)
 add_subdirectory(advection)
 add_subdirectory(timestepper)
 add_subdirectory(pde_solvers)

src/README.md

@@ -8,6 +8,7 @@ The `src/` folder contains all the code necessary to build a gyrokinetic semi-La
 - [geometryXYVxVy](./geometryXYVxVy/README.md) - Code describing methods which are specific to a simulation with 2 spatial dimensions and 2 velocity dimension.
 - [geometry5D](./geometry5D/README.md) - Code describing methods which are specific to a simulation with 3 spatial dimensions and 2 velocity dimension. This will be the general geometry used for GyselaX++ development.
 - [interpolation](./interpolation/README.md) - Code describing interpolation methods.
+- [multipatch](./multipatch/README.md) - Code describing multipatch geometry.
 <!-- - [paraconfpp](./paraconfpp/README.md) - Paraconf utility functions. -->
 - [pde\_solvers](./pde_solvers/README.md) - Code describing different methods for solving PDEs (e.g. Poisson's equation).
 - [quadrature](./quadrature/README.md) - Code describing different quadrature methods.

src/multipatch/CMakeLists.txt

@@ -0,0 +1,4 @@
+
+
+
+add_subdirectory(interfaces)

src/multipatch/README.md

@@ -0,0 +1,5 @@
+# Multipatch 
+
+The `multipatch` folder contains all the code describing methods and classes which are specific to a multipatch geometry. It is broken up into the following sub-folders:
+
+- [interface](./interfaces/README.md) - Defines structures and methods for sticking patches together.

src/multipatch/interfaces/CMakeLists.txt

@@ -0,0 +1,13 @@
+cmake_minimum_required(VERSION 3.15)
+
+
+add_library("multipatch_interface" INTERFACE)
+target_include_directories("multipatch_interface"
+    INTERFACE
+        "$<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}>"
+)
+target_link_libraries("multipatch_interface" INTERFACE
+    DDC::DDC
+    gslx::utils
+)
+add_library("gslx::multipatch_interface" ALIAS "multipatch_interface")

src/multipatch/interfaces/README.md

@@ -0,0 +1,106 @@
+# Interfaces 
+
+This directory defines structures and methods for the sticking of
+patches of a multi-patch domain.
+
+## Sticking of Two Edges
+
+We follow the idea to represent the topology of a multipatch domain using
+the idea of *domain manifolds* (see [1]). This means that we have several
+independent tensor-product logical domains and define their sticking via coordinate 
+transformations.
+
+
+**Remark:** In the referenced paper, the authors use affine isometries of the entire 
+patch instead of coordinate transformations of the faces. This is equivalent to our approach in 
+the sense that, when the scaling of the domains is removed and we consider only 
+unit squares as logical domains, the coordinate transformations of the faces and the information 
+about which edges are identified determine the 
+isometries used in the paper and vice versa. This is easy to show.
+
+### Multi-patch domain
+
+For simplicity, we will constrain ourselves to the 2D case, but this approach
+could be generalized to arbitrary dimensions.
+
+Let $\Omega$ be the domain of interest and assume that we have patches 
+$\Omega^{(i)}$, $i=1,...,K$, which are disjoint, s.t.
+$$
+\Omega = \dot{\bigcup}_{i=1}^K \Omega^{(i)}.
+$$
+For every patch $\Omega^{(i)}$ we have a mapping 
+$$
+F^{(i)}: [a_x^{(i)}, b_x^{(i)}]\times[a_y^{(i)} b_y^{(i)}] \rightarrow \Omega^{(i)}, 
+$$
+where the rectangular domain $`[a_x^{(i)}, b_x^{(i)}]\times[a_y^{(i)} b_y^{(i)}]`$ defines 
+the logical coordinates for the patch. We call this logical coordinate domain 
+the *logical patch*.
+
+### Edges
+
+So we obtain four logical edges
+$$
+[a_x^{(i)}, b_x^{(i)}] \times \{ a_y^{(i)} \}, \quad
+[a_x^{(i)}, b_x^{(i)}] \times \{ b_y^{(i)} \}, \quad
+\{ a_x^{(i)} \} \times [a_y^{(i)}, b_y^{(i)}], \quad
+\{ b_x^{(i)} \} \times [a_y^{(i)}, b_y^{(i)}] 
+$$
+
+In the code, we define edges as follows. Every edge of a 
+logical patch is identified via the patch it belongs to, the dimension 
+and whether it is at the front or the back of the domain. So e.g.
+the edge $`[a_x^{(i)}, b_x^{(i)}] \times \{ a_y^{(i)} \}`$ would be
+identified with patch $i$, dimensions `RDimXi` and `FRONT`.
+$`[a_x^{(i)}, b_x^{(i)}] \times \{ b_y^{(i)} \}`$ would be 
+identified with patch $i$, dimensions `RDimXi` and `BACK` and 
+$`\{ b_x^{(i)} \} \times [a_y^{(i)}, b_y^{(i)}]`$ would be 
+identified with patch $i$, dimensions `RDimYi` and `BACK`.
+
+### Sticking and Coordinate Transformation
+
+Any edge can then be identified with an edge on another patch. This corresponds
+to the 'sticking'. So for a different patch $`\Omega^{(j)}`$, we have the logical coordinate domain
+$`[a_x^{(j)}, b_x^{(j)}] \times [a_y^{(j)} b_y^{(j)}]`$.
+If we want to stick patch $i$ to patch $j$, we have to determine the edges that
+are identified and how they are identified. The way that they are identified
+is mathematically determined via the coordinate transformation from one edge
+to the other. Since the transformations are supposed to be affine and bijective, 
+there are only two options: The transformation can be order preserving or 
+order reversing (this corresponds to the orientation of the phyiscal edge where 
+two parametrizations coming from the two patches can have either the same or the 
+opposite orientation respectively).
+
+So for example, we want to stick the edge $`\{ a_x^{(i)} \} \times [a_y^{(i)}, b_y^{(i)}]`$
+on patch $i$ to the edge $`[a_x^{(j)}, b_x^{(j)}] \times \{ b_y^{(j)} \}`$ on patch 
+$j$. If the transformation is order-preserving (i.e. the orientations of the parametrizations 
+of the physical edge agree), then the transformation from the first edge to the second is 
+$$
+t \mapsto a_x^{(j)} + \frac{t - a_y^{(i)}}{b_y^{(i)} - a_y^{(i)}} \, (b_x^{(j)} - a_x^{(j)}).
+$$
+If the transformation is order-reversing (i.e. the orientations of the parametrizations 
+of the physical edge are opposite), then it is
+$$
+t \mapsto b_x^{(j)} - \frac{t - a_y^{(i)}}{b_y^{(i)} - a_y^{(i)}} \, (b_x^{(j)} - a_x^{(j)}).
+$$
+
+
+In the code, the sticking of two edges is represented by an `Interface` structure
+which contains references to 
+the first patch and the second patch as well as the boolean member
+`orientations_agree`. The transformation from one edge to the other
+is done using the 
+`EdgeCoordinatesTransformation` operator.
+
+## References
+[1] Buchegger, F., Jüttler, B., Mantzaflaris, A., 
+*Adaptively refined multi-patch B-splines with enhanced smoothness*.
+Applied Mathematics and Computation,
+Volume 272, Part 1
+(2016).
+https://www.sciencedirect.com/science/article/pii/S009630031500836X.
+
+
+## Contents
+- `coord_transformation.hpp`: Defines the `EdgeCoordinatesTransformation` operator to transform the coordinate from one edge to the other (see above).
+- `edge.hpp`: Defines the `Edge` structure.
+- `interface.hpp`: Defines the `Interface` structure.

src/multipatch/interfaces/coord_transformation.hpp

@@ -0,0 +1,114 @@
+// SPDX-License-Identifier: MIT
+
+#pragma once
+
+#include <ddc/ddc.hpp>
+
+#include <ddc_helper.hpp>
+
+#include "interface.hpp"
+
+/**
+ * @brief Transform a coordinate from one edge to a coordinate on the other edge.
+ * 
+ * According to the orientation stored in the interface, we compute 
+ * * if True, 
+ * 
+ *  @f$ t \mapsto min_2 + \frac{t - min_1}{max_1 - min_1}(max_2 - min_2) @f$
+ * 
+ * * if False, 
+ * 
+ *  @f$ t \mapsto max_2 - \frac{t - min_1}{max_1 - min_1}(max_2 - min_2) @f$
+ * 
+ * where @f$ min_i @f$ and @f$ max_i @f$ are the minimum and maximum 
+ * coordinates of the edge @f$ i @f$. 
+ * 
+ * @tparam IDim1
+ *          The discrete dimension of the first edge. 
+ * @tparam IDim2
+ *          The discrete dimension of the second edge. 
+ * 
+*/
+template <class IDim1, class IDim2>
+class EdgeCoordinatesTransformation
+{
+    using RDim1 = typename IDim1::continuous_dimension_type;
+    using RDim2 = typename IDim2::continuous_dimension_type;
+
+    Interface<IDim1, IDim2> const m_interface;
+
+public:
+    /**
+     * @brief Instantiate an EdgeCoordinatesTransformation.
+     * 
+     * @param interface 
+     *      An Interface object between the two patches. 
+    */
+    EdgeCoordinatesTransformation(Interface<IDim1, IDim2> const interface)
+        : m_interface(interface) {};
+    ~EdgeCoordinatesTransformation() = default;
+
+
+    /**
+     * @brief Transform a coordinate on the edge in the dimension 
+     * of the current patch to the analogous coordinate on the target patch. 
+     * 
+     * @param current_coord
+     *      A coordinate on the edge of the current patch.
+     * 
+     * @tparam CurrentRDim 
+     *      The current continuous dimension of the given coordinate coord. 
+     * 
+     * @return The analogous coordinate on the target patch. 
+    */
+    template <class CurrentRDim>
+    ddc::Coordinate<std::conditional_t<std::is_same_v<CurrentRDim, RDim1>, RDim2, RDim1>>
+    operator()(ddc::Coordinate<CurrentRDim> current_coord) const
+    {
+        // The discrete dimension of CurrentRDim
+        using CurrentIDim = std::conditional_t<std::is_same_v<CurrentRDim, RDim1>, IDim1, IDim2>;
+        // The other continuous dimension
+        using ORDim = std::conditional_t<std::is_same_v<CurrentRDim, RDim1>, RDim2, RDim1>;
+        // The other discrete dimension
+        using OIDim = std::conditional_t<std::is_same_v<CurrentIDim, IDim1>, IDim2, IDim1>;
+
+
+        bool const orientations_agree = m_interface.orientations_agree;
+
+        ddc::Coordinate<CurrentRDim> current_min;
+        double current_length;
+        get_min_and_length<CurrentIDim>(current_min, current_length);
+
+        ddc::Coordinate<ORDim> target_min;
+        double target_length;
+        get_min_and_length<OIDim>(target_min, target_length);
+
+        double rescale_x = (current_coord - current_min) / current_length * target_length;
+
+        if (!orientations_agree) {
+            rescale_x = target_length - rescale_x;
+        }
+        return target_min + rescale_x; // This has type ddc::Coordinate<ODim>.
+    };
+
+
+private:
+    /**
+     * @brief Get the minimum coordinate of a patch on the edge and its length.
+     * 
+     * @tparam IDim 
+     *      The discrete dimension of the edge of the current patch. 
+    */
+    template <class IDim>
+    void get_min_and_length(
+            ddc::Coordinate<typename IDim::continuous_dimension_type>& min,
+            double& length) const
+    {
+        Edge<IDim> const edge = m_interface.template get_edge<IDim>();
+
+        ddc::DiscreteDomain<IDim> const dom = edge.domain;
+
+        min = ddc::coordinate(dom.front());
+        length = ddcHelper::total_interval_length(dom);
+    };
+};

src/multipatch/interfaces/edge.hpp

@@ -0,0 +1,52 @@
+// SPDX-License-Identifier: MIT
+
+#pragma once
+
+#include <ddc/ddc.hpp>
+
+
+enum Extremity { FRONT, BACK };
+
+
+
+/**
+ * @brief Define an edge of a given patch.
+ * 
+ * An edge is defined by a patch index, a dimension
+ * and an extremity.
+ * For example, in the patch defined on logical domain 
+ * @f$ [a_x, b_x]\times[a_y, b_y] @f$, 
+ * 
+ * * the edge IDimX, BACK refers to the set @f$ [a_x, b_x]\times\{b_y\} @f$, 
+ * * and the edge IDimX, FRONT refers to the set @f$ [a_x, b_x]\times\{a_y\} @f$. 
+ * 
+ * Here, the domain given as input corresponds to @f$ [a_x, b_x] @f$
+ * in both case. 
+ * 
+ * 
+ * @tparam IDim
+ *      The discrete dimension of the edge. 
+*/
+template <class IDim>
+struct Edge
+{
+public:
+    /**
+     * The index of the patch. 
+    */
+    std::size_t patch_index;
+    /**
+     * The discrete domain of the edge. 
+     * See example with @f$ [a_x, b_x] @f$. 
+    */
+    ddc::DiscreteDomain<IDim> domain;
+    /**
+     * A Extremity type containing "BACK" or "FRONT" tag. 
+    */
+    Extremity extremity;
+
+    /**
+     * The discrete dimension of the edge. 
+    */
+    using discrete_dimension = IDim;
+};

src/multipatch/interfaces/interface.hpp

@@ -0,0 +1,56 @@
+// SPDX-License-Identifier: MIT
+
+#pragma once
+
+#include <ddc/ddc.hpp>
+
+#include "edge.hpp"
+
+/**
+ * @brief Represent a simple sticking of two edges.
+ * 
+ * @tparam IDim1 
+ *      The discrete dimension of the first edge. 
+ * @tparam IDim2
+ *      The discrete dimension of the second edge. 
+ *  
+ * @see EdgeCoordinatesTransformation
+*/
+template <class IDim1, class IDim2>
+struct Interface
+{
+    /**
+     * The Edge of the first patch. 
+    */
+    Edge<IDim1> edge_1;
+    /**
+     * The Edge of the second patch. 
+    */
+    Edge<IDim2> edge_2;
+    /**
+     * If True, the parametrisations of the physical edge have the same orientation.
+     * If False, the parametrisations of the physical edge have the opposite orientation.
+     * (See EdgeCoordinatesTransformation).
+    */
+    bool orientations_agree;
+
+    /**
+     * @brief Get the edge on the given dimension. 
+     * 
+     * @tparam IDim
+     *      The discrete dimension of the edge. 
+     * 
+     * @return The edge of the interface defined on the given dimension. 
+    */
+    template <class IDim>
+    constexpr Edge<IDim> get_edge() const
+    {
+        static_assert(((std::is_same_v<IDim, IDim1>) || (std::is_same_v<IDim, IDim2>)));
+
+        if constexpr (std::is_same_v<IDim, IDim1>) {
+            return edge_1;
+        } else {
+            return edge_2;
+        }
+    };
+};

tests/CMakeLists.txt

@@ -18,6 +18,7 @@ gtest_discover_tests(unit_tests_common DISCOVERY_MODE PRE_TEST)
 add_subdirectory(geometryXVx)
 add_subdirectory(geometryXYVxVy)
 add_subdirectory(geometryRTheta)
+add_subdirectory(multipatch)
 add_subdirectory(pde_solvers)
 add_subdirectory(timestepper)
 add_subdirectory(utils)