-
Notifications
You must be signed in to change notification settings - Fork 95
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
Merge pull request #3018 from boutproject/move-invert3x3
Move `invert3x3` out of general purpose `utils.hxx` header
- Loading branch information
Showing
9 changed files
with
164 additions
and
159 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
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
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,77 @@ | ||
| /*!************************************************************************* | ||
| * \file invert3x3.hxx | ||
| * | ||
| * A mix of short utilities for memory management, strings, and some | ||
| * simple but common calculations | ||
| * | ||
| ************************************************************************** | ||
| * Copyright 2010-2024 B.D.Dudson, BOUT++ Team | ||
| * | ||
| * Contact: Ben Dudson, [email protected] | ||
| * | ||
| * This file is part of BOUT++. | ||
| * | ||
| * BOUT++ is free software: you can redistribute it and/or modify | ||
| * it under the terms of the GNU Lesser General Public License as published by | ||
| * the Free Software Foundation, either version 3 of the License, or | ||
| * (at your option) any later version. | ||
| * | ||
| * BOUT++ is distributed in the hope that it will be useful, | ||
| * but WITHOUT ANY WARRANTY; without even the implied warranty of | ||
| * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the | ||
| * GNU Lesser General Public License for more details. | ||
| * | ||
| * You should have received a copy of the GNU Lesser General Public License | ||
| * along with BOUT++. If not, see <http://www.gnu.org/licenses/>. | ||
| * | ||
| **************************************************************************/ | ||
|
|
||
| #pragma once | ||
|
|
||
| #include <bout/utils.hxx> | ||
| #include <optional> | ||
|
|
||
| /// Explicit inversion of a 3x3 matrix \p a | ||
| /// | ||
| /// If the matrix is singular (ill conditioned), the determinant is | ||
| /// return. Otherwise, an empty `std::optional` is return | ||
| namespace bout { | ||
| inline std::optional<BoutReal> invert3x3(Matrix<BoutReal>& a) { | ||
| TRACE("invert3x3"); | ||
|
|
||
| // Calculate the first co-factors | ||
| const BoutReal A = a(1, 1) * a(2, 2) - a(1, 2) * a(2, 1); | ||
| const BoutReal B = a(1, 2) * a(2, 0) - a(1, 0) * a(2, 2); | ||
| const BoutReal C = a(1, 0) * a(2, 1) - a(1, 1) * a(2, 0); | ||
|
|
||
| // Calculate the determinant | ||
| const BoutReal det = a(0, 0) * A + a(0, 1) * B + a(0, 2) * C; | ||
| constexpr BoutReal small = 1.0e-15; | ||
| if (std::abs(det) < std::abs(small)) { | ||
| return det; | ||
| } | ||
|
|
||
| // Calculate the rest of the co-factors | ||
| const BoutReal D = a(0, 2) * a(2, 1) - a(0, 1) * a(2, 2); | ||
| const BoutReal E = a(0, 0) * a(2, 2) - a(0, 2) * a(2, 0); | ||
| const BoutReal F = a(0, 1) * a(2, 0) - a(0, 0) * a(2, 1); | ||
| const BoutReal G = a(0, 1) * a(1, 2) - a(0, 2) * a(1, 1); | ||
| const BoutReal H = a(0, 2) * a(1, 0) - a(0, 0) * a(1, 2); | ||
| const BoutReal I = a(0, 0) * a(1, 1) - a(0, 1) * a(1, 0); | ||
|
|
||
| // Now construct the output, overwrites input | ||
| const BoutReal detinv = 1.0 / det; | ||
|
|
||
| a(0, 0) = A * detinv; | ||
| a(0, 1) = D * detinv; | ||
| a(0, 2) = G * detinv; | ||
| a(1, 0) = B * detinv; | ||
| a(1, 1) = E * detinv; | ||
| a(1, 2) = H * detinv; | ||
| a(2, 0) = C * detinv; | ||
| a(2, 1) = F * detinv; | ||
| a(2, 2) = I * detinv; | ||
|
|
||
| return std::nullopt; | ||
| } | ||
| } // namespace bout |
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
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,74 @@ | ||
| #include "../../src/mesh/invert3x3.hxx" | ||
|
|
||
| #include "gtest/gtest.h" | ||
|
|
||
| TEST(Invert3x3Test, Identity) { | ||
| Matrix<BoutReal> input(3, 3); | ||
| input = 0; | ||
| for (int i = 0; i < 3; i++) { | ||
| input(i, i) = 1.0; | ||
| } | ||
| auto expected = input; | ||
| bout::invert3x3(input); | ||
|
|
||
| for (int j = 0; j < 3; j++) { | ||
| for (int i = 0; i < 3; i++) { | ||
| EXPECT_EQ(input(i, j), expected(i, j)); | ||
| } | ||
| } | ||
| } | ||
|
|
||
| TEST(Invert3x3Test, InvertTwice) { | ||
| std::vector<BoutReal> rawDataMat = {0.05567105, 0.92458227, 0.19954631, | ||
| 0.28581972, 0.54009039, 0.13234403, | ||
| 0.8841194, 0.161224, 0.74853209}; | ||
| std::vector<BoutReal> rawDataInv = {-2.48021781, 4.27410022, -0.09449605, | ||
| 0.6278449, 0.87275842, -0.32168092, | ||
| 2.79424897, -5.23628123, 1.51684677}; | ||
|
|
||
| Matrix<BoutReal> input(3, 3); | ||
| Matrix<BoutReal> expected(3, 3); | ||
|
|
||
| int counter = 0; | ||
| for (int j = 0; j < 3; j++) { | ||
| for (int i = 0; i < 3; i++) { | ||
| input(i, j) = rawDataMat[counter]; | ||
| expected(i, j) = rawDataInv[counter]; | ||
| counter++; | ||
| } | ||
| } | ||
|
|
||
| // Invert twice to check if we get back to where we started | ||
| bout::invert3x3(input); | ||
|
|
||
| for (int j = 0; j < 3; j++) { | ||
| for (int i = 0; i < 3; i++) { | ||
| // Note we only check to single tolerance here | ||
| EXPECT_FLOAT_EQ(input(i, j), expected(i, j)); | ||
| } | ||
| } | ||
| } | ||
|
|
||
| TEST(Invert3x3Test, Singular) { | ||
| Matrix<BoutReal> input(3, 3); | ||
| input = 0; | ||
| auto result = bout::invert3x3(input); | ||
| EXPECT_TRUE(result.has_value()); | ||
| } | ||
|
|
||
| TEST(Invert3x3Test, BadCondition) { | ||
| Matrix<BoutReal> input(3, 3); | ||
|
|
||
| input = 0.; | ||
| input(0, 0) = 1.0e-16; | ||
| input(1, 1) = 1.0; | ||
| input(2, 2) = 1.0; | ||
| EXPECT_TRUE(bout::invert3x3(input).has_value()); | ||
|
|
||
| // not quite bad enough condition | ||
| input = 0.; | ||
| input(0, 0) = 1.0e-12; | ||
| input(1, 1) = 1.0; | ||
| input(2, 2) = 1.0; | ||
| EXPECT_FALSE(bout::invert3x3(input).has_value()); | ||
| } |
Oops, something went wrong.