#1292

- Replace the computation of determinant by the native Xtensor function.
- Reactivate some asserts that were deactivated during the transition.
- Update some ticket numbers in the code, including fixing some that were incorrect.

Sources/FiniteElement/RefFiniteElement/Advanced/RefFEltInLocalOperator.cpp

@@ -72,7 +72,7 @@ namespace MoReFEM
         }
         
         
-        // \todo #1292 Probably might be done with native Xtensor functions
+        // \todo #1491 Probably might be done with native Xtensor functions
         const LocalMatrix& ExtractSubMatrix(const LocalMatrix& full_matrix,
                                             const RefFEltInLocalOperator& ref_felt)
         {

Sources/FiniteElement/RefFiniteElement/Advanced/RefFEltInLocalOperator.hpp

@@ -18,8 +18,6 @@
 # include <memory>
 # include <vector>
 
-// \todo #1292 # include ""
-
 # include "Utilities/LinearAlgebra/Storage/Local/LocalMatrixStorage.hpp"
 # include "Utilities/LinearAlgebra/Storage/Local/LocalVectorStorage.hpp"
 

Sources/Geometry/Mesh/Mesh.hpp

@@ -25,7 +25,6 @@
 # include "Utilities/Mpi/Mpi.hpp"
 # include "Utilities/Filesystem/Directory.hpp"
 
-// \todo #1292 # include ""
 # include "ThirdParty/IncludeWithoutWarning/Petsc/PetscSys.hpp"
 
 # include "Core/InputData/InputData.hpp"

Sources/OperatorInstances/ConformInterpolator/Local/P1_to_P1b.cpp

@@ -86,7 +86,7 @@ namespace MoReFEM
                 LocalMatrix block({ Nnode_in_row, Nnode_in_col });
                 block.fill(0.);
 
-				// \todo #1492 Xtensor probably has something better for that...
+				// \todo #1491 Xtensor probably has something better for that...
                 for (auto i = 0ul; i < Nnode_in_col; ++i)
                     block(i, i) = 1.;
 

Sources/Operators/GlobalVariationalOperator/Internal/Helper.hxx

@@ -173,21 +173,21 @@ namespace MoReFEM
 
                 const auto& row_local_2_global = local_felt_space.GetLocal2Global<MpiScale::program_wise>(row_test_extended_unknown_list);
 
-              // \todo #1292  assert(static_cast<int>(row_local_2_global.size()) == local_matrix.shape(0));
+				assert(row_local_2_global.size() == local_matrix.shape(0));
 
                 decltype(auto) col_extended_unknown_list =
                     global_variational_operator.GetExtendedUnknownList(col_numbering_subset);
                 
                 const auto& col_local_2_global = local_felt_space.GetLocal2Global<MpiScale::program_wise>(col_extended_unknown_list);
 
-              // \todo #1292  assert(static_cast<int>(col_local_2_global.size()) == local_matrix.shape(1));
+              	assert(col_local_2_global.size() == local_matrix.shape(1));
 
-                int row_local_matrix = 0;
-                int col_local_matrix = 0;
+                auto row_local_matrix = 0ul;
+                auto col_local_matrix = 0ul;
 
                 for (auto& row : row_local_2_global)
                 {
-                    col_local_matrix = 0;
+                    col_local_matrix = 0ul;
 
                     for (auto& col : col_local_2_global)
                     {
@@ -247,7 +247,7 @@ namespace MoReFEM
                 decltype(auto) local_2_global =
                     local_felt_space.GetLocal2Global<MpiScale::program_wise>(test_unknown_list);
 
-                // \todo #1292 assert(static_cast<int>(local_2_global.size()) == local_vector.GetSize());
+                assert(local_2_global.size() == local_vector.size());
                 if constexpr(std::is_same<std::decay_t<decltype(local_vector)>, LocalVector>())
                     global_vector.SetValues(local_2_global,
                                             local_vector.data(),

Sources/Operators/LocalVariationalOperator/Advanced/GradientDisplacementMatrix.cpp

@@ -32,7 +32,7 @@ namespace MoReFEM
             {
                 
                 
-                // \todo #1292 Evaluate if Xtensor can't do that much better...
+                // \todo #1491 Evaluate if Xtensor can't do that much better...
                 /*!
                  * \brief Perform matrix-vector product on a given subset of the input data.
                  *

Sources/Operators/LocalVariationalOperator/Advanced/InformationsAtQuadraturePoint/ForUnknownList.cpp

@@ -201,7 +201,7 @@ namespace MoReFEM
                         xt::noalias(GetNonCstGradientFEltPhi()) = xt::linalg::dot(GetGradientRefFEltPhi(), inverse_jacobian);
                     }
                     else
-                        GetNonCstDeterminant() = Wrappers::Xtensor::ComputeDeterminant(jacobian);
+                        GetNonCstDeterminant() = xt::linalg::det(jacobian);
                 }
                 else // Case in which a geometric with lower dimension than the mesh is considered.
                 {

Sources/ThirdParty/Wrappers/Xtensor/Functions.cpp

@@ -76,35 +76,6 @@ namespace MoReFEM::Wrappers::Xtensor
 	}
 
 
-	// \todo #1292 See to use native function for this!
-	double ComputeDeterminant(const LocalMatrix& local_matrix)
-	{
-		assert(local_matrix.shape(0) == local_matrix.shape(1) && "Matrix must be square!");
-		assert(local_matrix.shape(0) < 4 && "At the moment, only determinants for 2x2 or 3x3 matrices are "
-			   "computed.");
-
-		const auto dimension = local_matrix.shape(0);
-
-		switch(dimension)
-		{
-			case 1:
-				return local_matrix(0, 0);
-			case 2:
-				return local_matrix(0, 0) * local_matrix(1, 1) - local_matrix(0, 1) * local_matrix(1, 0);
-			case 3:
-				return local_matrix(0, 0) * (local_matrix(1, 1) * local_matrix(2, 2)
-											 - local_matrix(1, 2) * local_matrix(2, 1))
-				- local_matrix(0, 1) * (local_matrix(1, 0) * local_matrix(2, 2)
-										- local_matrix(1, 2) * local_matrix(2, 0))
-				+ local_matrix(0, 2) * (local_matrix(1, 0) * local_matrix(2, 1)
-										- local_matrix(1, 1) * local_matrix(2, 0));
-			default:
-				assert(false && "Not handled by this function!");
-				exit(EXIT_FAILURE);
-		}
-	}
-
-
 	void ComputeInverseSquareMatrix(const LocalMatrix& local_matrix, LocalMatrix& inverse, double& determinant)
 	{
 		const auto Nrow = local_matrix.shape(0);
@@ -112,7 +83,7 @@ namespace MoReFEM::Wrappers::Xtensor
 		assert(Nrow == inverse.shape(0));
 		assert(Nrow == inverse.shape(1));
 
-		determinant = ComputeDeterminant(local_matrix);
+		determinant = xt::linalg::det(local_matrix);
 
 		if (NumericNS::IsZero(determinant))
 			throw Exception("Determinant is zero!", __FILE__, __LINE__);

Sources/ThirdParty/Wrappers/Xtensor/Functions.hpp

@@ -77,18 +77,6 @@ namespace MoReFEM
              */
             bool IsZeroVector(const LocalVector& vector);
 
-
-			/*!
-			 * \brief Compute the determinant of \a local_matrix.
-			 *
-			 * Works only for dimensions 1 to 3.
-			 *
-			 * \param[in] local_matrix Matrix for which determinant is computed.
-			 *
-			 * \return Determinant of the matrix.
-			 */
-			double ComputeDeterminant(const LocalMatrix& local_matrix);
-
 			
 			/*!
 			 * \brief Check whether there is a nan or inf value in a \a LocalMatrix.