Skip to content
Snippets Groups Projects
Select Git revision
  • 9ab2f3e9543175c3d8cfdb6ddb7a6bdf683e798f
  • develop default protected
  • feature/gmsh-reader
  • origin/stage/bouguettaia
  • feature/kinetic-schemes
  • feature/reconstruction
  • feature/local-dt-fsi
  • feature/composite-scheme-sources
  • feature/composite-scheme-other-fluxes
  • feature/serraille
  • feature/variational-hydro
  • feature/composite-scheme
  • hyperplastic
  • feature/polynomials
  • feature/gks
  • feature/implicit-solver-o2
  • feature/coupling_module
  • feature/implicit-solver
  • feature/merge-local-dt-fsi
  • master protected
  • feature/escobar-smoother
  • v0.5.0 protected
  • v0.4.1 protected
  • v0.4.0 protected
  • v0.3.0 protected
  • v0.2.0 protected
  • v0.1.0 protected
  • Kidder
  • v0.0.4 protected
  • v0.0.3 protected
  • v0.0.2 protected
  • v0 protected
  • v0.0.1 protected
33 results

IBoundaryConditionDescriptor.hpp

Blame
    • Stéphane Del Pino's avatar
      9ab2f3e9
      Reorganize boundary conditions · 9ab2f3e9
      Stéphane Del Pino authored
      Now BC types are
      - Dirichlet
      - Fourier
      - Neumann
      - Symmetry
      
      One can set tags for Dirichlet, Neumann and Fourier in order to
      distinguish variables concerned by these conditions.
      
      In the future, it is planned to allow freefem like conditions such
      as `u = g on XMIN` or `alpha * u + dnu (u) = g on 3`,...
      
      However this will only make sense when "algorithms" will be set in
      pugs' language.
      9ab2f3e9
      History
      Reorganize boundary conditions
      Stéphane Del Pino authored
      Now BC types are
      - Dirichlet
      - Fourier
      - Neumann
      - Symmetry
      
      One can set tags for Dirichlet, Neumann and Fourier in order to
      distinguish variables concerned by these conditions.
      
      In the future, it is planned to allow freefem like conditions such
      as `u = g on XMIN` or `alpha * u + dnu (u) = g on 3`,...
      
      However this will only make sense when "algorithms" will be set in
      pugs' language.
    test_TinyMatrix.cpp 12.14 KiB
    #include <catch2/catch_approx.hpp>
    #include <catch2/catch_test_macros.hpp>
    #include <catch2/matchers/catch_matchers_all.hpp>
    
    #include <Kokkos_Core.hpp>
    
    #include <utils/PugsAssert.hpp>
    #include <utils/Types.hpp>
    
    #include <algebra/TinyMatrix.hpp>
    
    #include <sstream>
    
    // Instantiate to ensure full coverage is performed
    template class TinyMatrix<1, 1, int>;
    template class TinyMatrix<2, 2, int>;
    template class TinyMatrix<3, 4, int>;
    template class TinyMatrix<4, 4, double>;
    
    // clazy:excludeall=non-pod-global-static
    
    TEST_CASE("TinyMatrix", "[algebra]")
    {
      REQUIRE(TinyMatrix<1, 1, int>::Dimension == 1);
      REQUIRE(TinyMatrix<1, 1, int>::NumberOfRows == 1);
      REQUIRE(TinyMatrix<1, 1, int>::NumberOfColumns == 1);
      REQUIRE(TinyMatrix<2, 3, int>::Dimension == 6);
      REQUIRE(TinyMatrix<2, 3, int>::NumberOfRows == 2);
      REQUIRE(TinyMatrix<2, 3, int>::NumberOfColumns == 3);
      REQUIRE(TinyMatrix<5, 4, int>::Dimension == 20);
      REQUIRE(TinyMatrix<5, 4, int>::NumberOfRows == 5);
      REQUIRE(TinyMatrix<5, 4, int>::NumberOfColumns == 4);
    
      TinyMatrix<3, 4, int> A(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12);
      TinyVector<3, int> x1(1, 5, 9);
      TinyVector<3, int> x2(2, 6, 10);
      TinyVector<3, int> x3(3, 7, 11);
      TinyVector<3, int> x4(4, 8, 12);
      REQUIRE(((A(0, 0) == 1) and (A(0, 1) == 2) and (A(0, 2) == 3) and (A(0, 3) == 4) and   //
               (A(1, 0) == 5) and (A(1, 1) == 6) and (A(1, 2) == 7) and (A(1, 3) == 8) and   //
               (A(2, 0) == 9) and (A(2, 1) == 10) and (A(2, 2) == 11) and (A(2, 3) == 12)));
      REQUIRE(A == TinyMatrix<3, 4, int>(x1, x2, x3, x4));
      TinyMatrix<3, 4, int> B(6, 5, 3, 8, 34, 6, 35, 6, 7, 1, 3, 6);
    
      SECTION("checking for opposed matrix")
      {
        const TinyMatrix minus_A = -A;
        REQUIRE(
          ((minus_A(0, 0) == -1) and (minus_A(0, 1) == -2) and (minus_A(0, 2) == -3) and (minus_A(0, 3) == -4) and   //
           (minus_A(1, 0) == -5) and (minus_A(1, 1) == -6) and (minus_A(1, 2) == -7) and (minus_A(1, 3) == -8) and   //
           (minus_A(2, 0) == -9) and (minus_A(2, 1) == -10) and (minus_A(2, 2) == -11) and (minus_A(2, 3) == -12)));
      }
    
      SECTION("checking for equality and difference tests")
      {
        const TinyMatrix copy_A = A;
        REQUIRE(((copy_A(0, 0) == 1) and (copy_A(0, 1) == 2) and (copy_A(0, 2) == 3) and (copy_A(0, 3) == 4) and   //
                 (copy_A(1, 0) == 5) and (copy_A(1, 1) == 6) and (copy_A(1, 2) == 7) and (copy_A(1, 3) == 8) and   //
                 (copy_A(2, 0) == 9) and (copy_A(2, 1) == 10) and (copy_A(2, 2) == 11) and (copy_A(2, 3) == 12)));
        REQUIRE(copy_A == A);
        REQUIRE_FALSE(copy_A != A);
    
        TinyMatrix<3, 4, int> affected_A;
        affected_A = A;
        REQUIRE(affected_A == A);
        REQUIRE_FALSE(affected_A != A);
    
        REQUIRE(A != B);
        REQUIRE_FALSE(A == B);
      }
    
      SECTION("checking for scalar left product")
      {
        const int a                    = 2;
        const TinyMatrix<3, 4, int> aA = a * A;
    
        REQUIRE(aA == TinyMatrix<3, 4, int>(2, 4, 6, 8, 10, 12, 14, 16, 18, 20, 22, 24));
      }
    
      SECTION("checking for scalar seft product")
      {
        const int a = 2;
    
        TinyMatrix copy_A = A;
    
        REQUIRE((copy_A *= a) == TinyMatrix<3, 4, int>(2, 4, 6, 8, 10, 12, 14, 16, 18, 20, 22, 24));
      }
    
      SECTION("checking for null matrix management")
      {
        TinyMatrix<4, 3, int> Z = zero;
        REQUIRE(Z == TinyMatrix<4, 3, int>(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0));
    
        TinyMatrix<4, 3, int> affected_Z;
        affected_Z = zero;
        REQUIRE(affected_Z == Z);
      }
    
      SECTION("checking for identity management")
      {
        TinyMatrix<3, 3, int> I = identity;
        REQUIRE(I == TinyMatrix<3, 3, int>(1, 0, 0, 0, 1, 0, 0, 0, 1));
    
        TinyMatrix<3, 3, int> affected_I;
        affected_I = identity;
        REQUIRE(affected_I == I);
      }
    
      SECTION("checking for matrices sum")
      {
        REQUIRE(A + B == TinyMatrix<3, 4, int>(7, 7, 6, 12, 39, 12, 42, 14, 16, 11, 14, 18));
    
        TinyMatrix<3, 4, int> ApB = A;
        ApB += B;
        REQUIRE(ApB == A + B);
    
        TinyMatrix<3, 4, int> Ap2B = A + 2 * B;
        REQUIRE(Ap2B == ApB + B);
      }
    
      SECTION("checking for matrices difference ")
      {
        REQUIRE(A - B == TinyMatrix<3, 4, int>(-5, -3, 0, -4, -29, 0, -28, 2, 2, 9, 8, 6));
    
        TinyMatrix<3, 4, int> AmB = A;
        AmB -= B;
        REQUIRE(AmB == A - B);
    
        TinyMatrix<3, 4, int> Am2B = A - 2 * B;
        REQUIRE(Am2B == AmB - B);
      }
    
      SECTION("checking for matrices product")
      {
        TinyMatrix<4, 2, int> C{3, -2, 2, 6, -2, 5, 7, 2};
        REQUIRE(A * C == TinyMatrix<3, 2, int>(29, 33, 69, 77, 109, 121));
    
        TinyMatrix<2, 3, int> D{-3, 2, 3, 2, 5, -7};
        REQUIRE(D * A == TinyMatrix<2, 4, int>(34, 36, 38, 40, -36, -36, -36, -36));
      }
    
      SECTION("checking for matrix-vector product")
      {
        REQUIRE(A * TinyVector<4, int>(2, -3, 5, 2) == TinyVector<3, int>(19, 43, 67));
      }
    
      SECTION("checking for tensor product")
      {
        const TinyVector<4, int> u(1, 3, 7, 5);
        const TinyVector<3, int> v(6, 2, -3);
    
        REQUIRE(tensorProduct(u, v) == TinyMatrix<4, 3, int>(6, 2, -3, 18, 6, -9, 42, 14, -21, 30, 10, -15));
      }
    
      SECTION("checking for minor calculation")
      {
        REQUIRE(getMinor(A, 0, 0) == TinyMatrix<2, 3, int>(6, 7, 8, 10, 11, 12));
        REQUIRE(getMinor(A, 1, 0) == TinyMatrix<2, 3, int>(2, 3, 4, 10, 11, 12));
        REQUIRE(getMinor(A, 2, 0) == TinyMatrix<2, 3, int>(2, 3, 4, 6, 7, 8));
    
        REQUIRE(getMinor(A, 0, 1) == TinyMatrix<2, 3, int>(5, 7, 8, 9, 11, 12));
        REQUIRE(getMinor(A, 1, 1) == TinyMatrix<2, 3, int>(1, 3, 4, 9, 11, 12));
        REQUIRE(getMinor(A, 2, 1) == TinyMatrix<2, 3, int>(1, 3, 4, 5, 7, 8));
    
        REQUIRE(getMinor(A, 0, 2) == TinyMatrix<2, 3, int>(5, 6, 8, 9, 10, 12));
        REQUIRE(getMinor(A, 1, 2) == TinyMatrix<2, 3, int>(1, 2, 4, 9, 10, 12));
        REQUIRE(getMinor(A, 2, 2) == TinyMatrix<2, 3, int>(1, 2, 4, 5, 6, 8));
    
        REQUIRE(getMinor(A, 0, 3) == TinyMatrix<2, 3, int>(5, 6, 7, 9, 10, 11));
        REQUIRE(getMinor(A, 1, 3) == TinyMatrix<2, 3, int>(1, 2, 3, 9, 10, 11));
        REQUIRE(getMinor(A, 2, 3) == TinyMatrix<2, 3, int>(1, 2, 3, 5, 6, 7));
      }
    
      SECTION("checking for cofactors")
      {
        TinyMatrix<3, 3, int> A(1, 2, 3, 4, 5, 6, 7, 8, 9);
        REQUIRE(cofactor(A, 0, 0) == (5 * 9 - 8 * 6));
        REQUIRE(cofactor(A, 1, 0) == -(2 * 9 - 8 * 3));
        REQUIRE(cofactor(A, 2, 0) == (2 * 6 - 5 * 3));
    
        REQUIRE(cofactor(A, 0, 1) == -(4 * 9 - 7 * 6));
        REQUIRE(cofactor(A, 1, 1) == (1 * 9 - 7 * 3));
        REQUIRE(cofactor(A, 2, 1) == -(1 * 6 - 4 * 3));
    
        REQUIRE(cofactor(A, 0, 2) == (4 * 8 - 5 * 7));
        REQUIRE(cofactor(A, 1, 2) == -(1 * 8 - 7 * 2));
        REQUIRE(cofactor(A, 2, 2) == (1 * 5 - 4 * 2));
      }
    
      SECTION("checking for determinant calculations")
      {
        REQUIRE(det(TinyMatrix<1, 1, int>(6)) == 6);
        REQUIRE(det(TinyMatrix<2, 2, int>(3, 1, -3, 6)) == 21);
        REQUIRE(det(TinyMatrix<3, 3, int>(1, 1, 1, 1, 2, 1, 2, 1, 3)) == 1);
        REQUIRE(det(TinyMatrix<3, 3, int>(6, 5, 3, 8, 34, 6, 35, 6, 7)) == -1444);
        REQUIRE(det(TinyMatrix<4>(1, 2.3, 7, -6.2, 3, 4, 9, 1, 4.1, 5, 2, -3, 2, 27, 3, 17.5)) ==
                Catch::Approx(6661.455).epsilon(1E-14));
        REQUIRE(det(TinyMatrix<4>(1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 2, 0, 0, 2, 2)) == 0);
      }
    
      SECTION("checking for trace calculations")
      {
        REQUIRE(trace(TinyMatrix<1, 1, int>(6)) == 6);
        REQUIRE(trace(TinyMatrix<2, 2, int>(5, 1, -3, 6)) == 5 + 6);
        REQUIRE(trace(TinyMatrix<3, 3, int>(1, 1, 1, 1, 2, 1, 2, 1, 3)) == 1 + 2 + 3);
        REQUIRE(trace(TinyMatrix<3, 3, int>(6, 5, 3, 8, 34, 6, 35, 6, 7)) == 6 + 34 + 7);
        REQUIRE(trace(TinyMatrix<4>(1, 2.3, 7, -6.2, 3, 4, 9, 1, 4.1, 5, 2, -3, 2, 27, 3, 17.5)) ==
                Catch::Approx(1 + 4 + 2 + 17.5).epsilon(1E-14));
        REQUIRE(trace(TinyMatrix<4>(1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 2, 0, 0, 2, 2)) == 1 + 0 + 1 + 2);
      }
    
      SECTION("checking for inverse calculations")
      {
        {
          const TinyMatrix<1> A1(2);
          REQUIRE(inverse(A1)(0, 0) == Catch::Approx(0.5).epsilon(1E-14));
        }
    
        {
          const TinyMatrix<2> A2(2, 3, 4, 1);
          const TinyMatrix<2> I = inverse(A2) * A2;
    
          REQUIRE(I(0, 0) == Catch::Approx(1).epsilon(1E-14));
          REQUIRE(I(0, 1) == Catch::Approx(0).margin(1E-14));
          REQUIRE(I(1, 0) == Catch::Approx(0).margin(1E-14));
          REQUIRE(I(1, 1) == Catch::Approx(1).epsilon(1E-14));
        }
    
        {
          const TinyMatrix<3> A3(2, 3, 1, 4, -1, 5, -2, 3, 4);
          const TinyMatrix<3> I = inverse(A3) * A3;
    
          REQUIRE(I(0, 0) == Catch::Approx(1).epsilon(1E-14));
          REQUIRE(I(0, 1) == Catch::Approx(0).margin(1E-14));
          REQUIRE(I(0, 2) == Catch::Approx(0).margin(1E-14));
          REQUIRE(I(1, 0) == Catch::Approx(0).margin(1E-14));
          REQUIRE(I(1, 1) == Catch::Approx(1).epsilon(1E-14));
          REQUIRE(I(1, 2) == Catch::Approx(0).margin(1E-14));
          REQUIRE(I(2, 0) == Catch::Approx(0).margin(1E-14));
          REQUIRE(I(2, 1) == Catch::Approx(0).margin(1E-14));
          REQUIRE(I(2, 2) == Catch::Approx(1).epsilon(1E-14));
        }
      }
    
      SECTION("checking for sizes")
      {
        REQUIRE(TinyMatrix<1>{}.numberOfRows() == 1);
        REQUIRE(TinyMatrix<1>{}.numberOfColumns() == 1);
        REQUIRE(TinyMatrix<1>{}.dimension() == 1);
        REQUIRE(TinyMatrix<1>{}.numberOfValues() == 1);
    
        REQUIRE(TinyMatrix<2>{}.numberOfRows() == 2);
        REQUIRE(TinyMatrix<2>{}.numberOfColumns() == 2);
        REQUIRE(TinyMatrix<2>{}.numberOfValues() == 4);
    
        REQUIRE(TinyMatrix<3>{}.numberOfRows() == 3);
        REQUIRE(TinyMatrix<3>{}.numberOfColumns() == 3);
        REQUIRE(TinyMatrix<3>{}.dimension() == 9);
        REQUIRE(TinyMatrix<3>{}.numberOfValues() == 9);
    
        REQUIRE(TinyMatrix<3, 4>{}.numberOfRows() == 3);
        REQUIRE(TinyMatrix<3, 4>{}.numberOfColumns() == 4);
        REQUIRE(TinyMatrix<3, 4>{}.dimension() == 12);
        REQUIRE(TinyMatrix<3, 4>{}.numberOfValues() == 12);
      }
    
      SECTION("is square matrix")
      {
        REQUIRE(TinyMatrix<1>{}.isSquare());
        REQUIRE(TinyMatrix<2>{}.isSquare());
        REQUIRE(TinyMatrix<3>{}.isSquare());
    
        REQUIRE_FALSE(TinyMatrix<3, 4>{}.isSquare());
      }
    
      SECTION("transpose")
      {
        TinyMatrix tA = transpose(A);
    
        REQUIRE(((tA(0, 0) == 1) and (tA(1, 0) == 2) and (tA(2, 0) == 3) and (tA(3, 0) == 4) and   //
                 (tA(0, 1) == 5) and (tA(1, 1) == 6) and (tA(2, 1) == 7) and (tA(3, 1) == 8) and   //
                 (tA(0, 2) == 9) and (tA(1, 2) == 10) and (tA(2, 2) == 11) and (tA(3, 2) == 12)));
    
        TinyMatrix ttA = transpose(tA);
        REQUIRE(ttA == A);
      }
    
      SECTION("checking for matrices output")
      {
        REQUIRE(Catch::Detail::stringify(A) == "[[1,2,3,4],[5,6,7,8],[9,10,11,12]]");
        REQUIRE(Catch::Detail::stringify(TinyMatrix<1, 1, int>(7)) == "[[7]]");
      }
    
      SECTION("checking scalarProduct")
      {
        TinyMatrix<3, 4, int> B(0, -1, 1, -1, 1, -1, 1, -1, 1, -1, 1, -1);
        //  TinyMatrix<3, 4, int> A(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12);
    
        REQUIRE(scalarProduct(A, B) == -7);
      }
      SECTION("checking norm")
      {
        TinyMatrix<3, 4, int> B(0, 0, 0, -1, 1, -1, 1, -1, 1, -1, 1, -1);
        //  TinyMatrix<3, 4, int> A(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12);
    
        REQUIRE(norm(B) == 3);
      }
    #ifndef NDEBUG
      SECTION("output with signaling NaN")
      {
        TinyMatrix<2, 3> A;
        A(0, 0) = 1;
        A(0, 2) = 3;
        A(1, 0) = 2;
        std::ostringstream A_ost;
        A_ost << A;
    
        std::ostringstream ref_ost;
        ref_ost << "[[1,nan,3],[2,nan,nan]]";
    
        REQUIRE(A_ost.str() == ref_ost.str());
      }
    
      SECTION("checking for bounds violation")
      {
        REQUIRE_THROWS_AS(A(3, 0), AssertError);
        REQUIRE_THROWS_AS(A(0, 4), AssertError);
    
        REQUIRE_THROWS_AS(getMinor(A, 3, 0), AssertError);
        REQUIRE_THROWS_AS(getMinor(A, 0, 4), AssertError);
    
        const TinyMatrix<3, 4, int>& constA = A;
        REQUIRE_THROWS_AS(constA(3, 0), AssertError);
        REQUIRE_THROWS_AS(constA(0, 4), AssertError);
      }
    
      SECTION("checking for nan initialization")
      {
        TinyMatrix<3, 4, double> B;
    
        for (size_t i = 0; i < B.numberOfRows(); ++i) {
          for (size_t j = 0; j < B.numberOfColumns(); ++j) {
            REQUIRE(std::isnan(B(i, j)));
          }
        }
      }
    
      SECTION("checking for bad initialization")
      {
        TinyMatrix<3, 4, int> B;
    
        for (size_t i = 0; i < B.numberOfRows(); ++i) {
          for (size_t j = 0; j < B.numberOfColumns(); ++j) {
            REQUIRE(B(i, j) == std::numeric_limits<int>::max() / 2);
          }
        }
      }
      SECTION("checking for bad initialization")
      {
        TinyMatrix<3, 4, int> B;
    
        for (size_t i = 0; i < B.numberOfRows(); ++i) {
          for (size_t j = 0; j < B.numberOfColumns(); ++j) {
            REQUIRE(B(i, j) == std::numeric_limits<int>::max() / 2);
          }
        }
      }
    #endif   // NDEBUG
    }