-
Notifications
You must be signed in to change notification settings - Fork 6
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
Create an `InverseJacobianMatrix` operator. This operator removes the need for the remaining inheritance of jacobian. Closes #410 . See merge request gysela-developpers/gyselalibxx!782 --------------------------------------------
- Loading branch information
1 parent
dffd882
commit b6800fa
Showing
9 changed files
with
219 additions
and
350 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
148 changes: 148 additions & 0 deletions
148
vendor/sll/include/sll/mapping/inverse_jacobian_matrix.hpp
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,148 @@ | ||
| // SPDX-License-Identifier: MIT | ||
| #pragma once | ||
|
|
||
| #include <sll/view.hpp> | ||
|
|
||
| #include "mapping_tools.hpp" | ||
|
|
||
| /** | ||
| * A class to calculate the inverse of the Jacobian matrix. | ||
| * If specialised methods are available then these are used by the class. | ||
| * Otherwise the inverse is calculated from the Jacobian matrix. | ||
| * | ||
| * @tparam Mapping The mapping whose inverse we are interested in. | ||
| * @tparam PositionCoordinate The coordinate system in which the inverse should be calculated. | ||
| */ | ||
| template <class Mapping, class PositionCoordinate = typename Mapping::CoordArg> | ||
| class InverseJacobianMatrix | ||
| { | ||
| private: | ||
| Mapping m_mapping; | ||
|
|
||
| public: | ||
| /** | ||
| * @brief A constructor for the InverseJacobianMatrix. | ||
| * @param[in] mapping The mapping whose inverse we are interested in. | ||
| */ | ||
| KOKKOS_FUNCTION explicit InverseJacobianMatrix(Mapping const& mapping) : m_mapping(mapping) {} | ||
|
|
||
| /** | ||
| * @brief Compute full inverse Jacobian matrix. | ||
| * | ||
| * For some computations, we need the complete inverse Jacobian matrix or just the | ||
| * coefficients. | ||
| * The coefficients can be given indendently with the functions | ||
| * inv_jacobian_11, inv_jacobian_12, inv_jacobian_21 and inv_jacobian_22. | ||
| * | ||
| * @param[in] coord The coordinate where we evaluate the Jacobian matrix. | ||
| * @returns The inverse Jacobian matrix returned. | ||
| * | ||
| * @see inv_jacobian_11 | ||
| * @see inv_jacobian_12 | ||
| * @see inv_jacobian_21 | ||
| * @see inv_jacobian_22 | ||
| */ | ||
| KOKKOS_INLINE_FUNCTION Matrix_2x2 operator()(PositionCoordinate const& coord) const | ||
| { | ||
| Matrix_2x2 matrix; | ||
| if constexpr (has_2d_inv_jacobian_v<Mapping, PositionCoordinate>) { | ||
| m_mapping.inv_jacobian_matrix(coord, matrix); | ||
| } else { | ||
| static_assert(has_2d_jacobian_v<Mapping, PositionCoordinate>); | ||
| double jacob = m_mapping.jacobian(coord); | ||
| assert(fabs(jacob) > 1e-15); | ||
| matrix[0][0] = m_mapping.jacobian_22(coord) / jacob; | ||
| matrix[0][1] = -m_mapping.jacobian_12(coord) / jacob; | ||
| matrix[1][0] = -m_mapping.jacobian_21(coord) / jacob; | ||
| matrix[1][1] = m_mapping.jacobian_11(coord) / jacob; | ||
| } | ||
| return matrix; | ||
| } | ||
|
|
||
| /** | ||
| * @brief Compute the (1,1) coefficient of the inverse Jacobian matrix. | ||
| * | ||
| * Be careful because not all mappings are invertible, especially at the center point. | ||
| * | ||
| * @param[in] coord | ||
| * The coordinate where we evaluate the inverse Jacobian matrix. | ||
| * | ||
| * @return A double with the value of the (1,1) coefficient of the inverse Jacobian matrix. | ||
| */ | ||
| KOKKOS_INLINE_FUNCTION double inv_jacobian_11(PositionCoordinate const& coord) const | ||
| { | ||
| if constexpr (has_2d_inv_jacobian_v<Mapping, PositionCoordinate>) { | ||
| return m_mapping.inv_jacobian_11(coord); | ||
| } else { | ||
| static_assert(has_2d_jacobian_v<Mapping, PositionCoordinate>); | ||
| double jacob = m_mapping.jacobian(coord); | ||
| assert(fabs(jacob) > 1e-15); | ||
| return m_mapping.jacobian_22(coord) / jacob; | ||
| } | ||
| } | ||
|
|
||
| /** | ||
| * @brief Compute the (1,2) coefficient of the inverse Jacobian matrix. | ||
| * | ||
| * Be careful because not all mappings are invertible, especially at the center point. | ||
| * | ||
| * @param[in] coord | ||
| * The coordinate where we evaluate the inverse Jacobian matrix. | ||
| * | ||
| * @return A double with the value of the (1,2) coefficient of the inverse Jacobian matrix. | ||
| */ | ||
| KOKKOS_INLINE_FUNCTION double inv_jacobian_12(PositionCoordinate const& coord) const | ||
| { | ||
| if constexpr (has_2d_inv_jacobian_v<Mapping, PositionCoordinate>) { | ||
| return m_mapping.inv_jacobian_12(coord); | ||
| } else { | ||
| static_assert(has_2d_jacobian_v<Mapping, PositionCoordinate>); | ||
| double jacob = m_mapping.jacobian(coord); | ||
| assert(fabs(jacob) > 1e-15); | ||
| return -m_mapping.jacobian_12(coord) / jacob; | ||
| } | ||
| } | ||
| /** | ||
| * @brief Compute the (2,1) coefficient of the inverse Jacobian matrix. | ||
| * | ||
| * Be careful because not all mappings are invertible, especially at the center point. | ||
| * | ||
| * @param[in] coord | ||
| * The coordinate where we evaluate the inverse Jacobian matrix. | ||
| * | ||
| * @return A double with the value of the (2,1) coefficient of the inverse Jacobian matrix. | ||
| */ | ||
| KOKKOS_INLINE_FUNCTION double inv_jacobian_21(PositionCoordinate const& coord) const | ||
| { | ||
| if constexpr (has_2d_inv_jacobian_v<Mapping, PositionCoordinate>) { | ||
| return m_mapping.inv_jacobian_21(coord); | ||
| } else { | ||
| static_assert(has_2d_jacobian_v<Mapping, PositionCoordinate>); | ||
| double jacob = m_mapping.jacobian(coord); | ||
| assert(fabs(jacob) > 1e-15); | ||
| return -m_mapping.jacobian_21(coord) / jacob; | ||
| } | ||
| } | ||
|
|
||
| /** | ||
| * @brief Compute the (2,2) coefficient of the inverse Jacobian matrix. | ||
| * | ||
| * Be careful because not all mappings are invertible, especially at the center point. | ||
| * | ||
| * @param[in] coord | ||
| * The coordinate where we evaluate the inverse Jacobian matrix. | ||
| * | ||
| * @return A double with the value of the (2,2) coefficient of the inverse Jacobian matrix. | ||
| */ | ||
| KOKKOS_INLINE_FUNCTION double inv_jacobian_22(PositionCoordinate const& coord) const | ||
| { | ||
| if constexpr (has_2d_inv_jacobian_v<Mapping, PositionCoordinate>) { | ||
| return m_mapping.inv_jacobian_22(coord); | ||
| } else { | ||
| static_assert(has_2d_jacobian_v<Mapping, PositionCoordinate>); | ||
| double jacob = m_mapping.jacobian(coord); | ||
| assert(fabs(jacob) > 1e-15); | ||
| return m_mapping.jacobian_11(coord) / jacob; | ||
| } | ||
| } | ||
| }; |
Oops, something went wrong.