Commit bab35660 authored by GILLES Sebastien's avatar GILLES Sebastien
Browse files

#1292 Use Xtensor for the mass operator. Currently I'm introducing the...

#1292 Use Xtensor for the mass operator. Currently I'm introducing the possibility to use Seldon or Xtensor; it is only temporary and the goal is to use Xtensor everywhere at the end of the day. Some asserts are deactivated and should be adapted and reactivated before closing this ticket.
parent 2fdea851
......@@ -105,7 +105,7 @@ install(DIRECTORY ${CMAKE_CURRENT_LIST_DIR}/ DESTINATION ${MOREFEM_INSTALL_DIR_I
#
add_library(morefem_cmake INTERFACE)
target_include_directories(morefem_cmake INTERFACE $<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}> $<INSTALL_INTERFACE:${PROJECT_NAME}/include> ${OPEN_MPI_INCL_DIR} ${PETSC_INCL_DIR} ${PARMETIS_INCL_DIR} ${LUA_INCL_DIR} ${BOOST_INCL_DIR} ${LIBMESH_INCL_DIR} $<INSTALL_INTERFACE:${PROJECT_NAME}/include/ThirdParty/Source/Tclap/include> $<INSTALL_INTERFACE:${PROJECT_NAME}/include/ThirdParty/Source/Seldon>)
target_include_directories(morefem_cmake INTERFACE $<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}> $<INSTALL_INTERFACE:${PROJECT_NAME}/include> ${OPEN_MPI_INCL_DIR} ${PETSC_INCL_DIR} ${PARMETIS_INCL_DIR} ${LUA_INCL_DIR} ${BOOST_INCL_DIR} ${LIBMESH_INCL_DIR} ${XTENSOR_INC_DIR} $<INSTALL_INTERFACE:${PROJECT_NAME}/include/ThirdParty/Source/Tclap/include> $<INSTALL_INTERFACE:${PROJECT_NAME}/include/ThirdParty/Source/Seldon>)
target_link_libraries(morefem_cmake INTERFACE ${MOREFEM_POST_PROCESSING} ${MOREFEM_TEST_TOOLS})
......
......@@ -36,7 +36,7 @@ namespace MoReFEM
/*!
* \brief Implementation of mass operator.
*/
class Mass final : public BilinearLocalVariationalOperator<LocalMatrix>
class Mass final : public BilinearLocalVariationalOperator<XtensorMatrix>
{
public:
......
......@@ -173,14 +173,14 @@ namespace MoReFEM
const auto& row_local_2_global = local_felt_space.GetLocal2Global<MpiScale::program_wise>(row_test_extended_unknown_list);
assert(static_cast<int>(row_local_2_global.size()) == local_matrix.GetM());
// \todo #1292 assert(static_cast<int>(row_local_2_global.size()) == local_matrix.GetM());
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);
assert(static_cast<int>(col_local_2_global.size()) == local_matrix.GetN());
// \todo #1292 assert(static_cast<int>(col_local_2_global.size()) == local_matrix.GetN());
int row_local_matrix = 0;
int col_local_matrix = 0;
......
......@@ -76,8 +76,17 @@ namespace MoReFEM
{
assert(Ndof_row > 0);
assert(Ndof_col > 0);
matrix_.Resize(static_cast<int>(Ndof_row), static_cast<int>(Ndof_col));
matrix_.Zero();
if constexpr(std::is_same<MatrixTypeT, XtensorMatrix>())
{
matrix_ = XtensorMatrix({ Ndof_row, Ndof_col });
matrix_.fill(0.);
}
else
{
matrix_.Resize(static_cast<int>(Ndof_row), static_cast<int>(Ndof_col));
matrix_.Zero();
}
}
......
......@@ -140,7 +140,10 @@ namespace MoReFEM
>
void ElementaryData<OperatorNS::Nature::bilinear, MatrixTypeT, std::false_type>::Zero()
{
this->GetNonCstMatrixResult().Zero();
if constexpr(std::is_same<MatrixTypeT, XtensorMatrix>())
this->GetNonCstMatrixResult().fill(0.);
else
this->GetNonCstMatrixResult().Zero();
}
......
......@@ -8,8 +8,8 @@
*/
#ifndef MOREFEM_x_TEST_x_THIRD_PARTY_x_P_E_T_SC_x_VECTOR_I_O_x_INPUT_DATA_HPP_
# define MOREFEM_x_TEST_x_THIRD_PARTY_x_P_E_T_SC_x_VECTOR_I_O_x_INPUT_DATA_HPP_
#ifndef MOREFEM_x_TEST_x_THIRD_PARTY_x_P_E_T_SC_x_MATRIX_OPERATIONS_x_INPUT_DATA_HPP_
# define MOREFEM_x_TEST_x_THIRD_PARTY_x_P_E_T_SC_x_MATRIX_OPERATIONS_x_INPUT_DATA_HPP_
# include "Utilities/Containers/EnumClass.hpp"
......@@ -61,4 +61,4 @@ namespace MoReFEM::TestNS::PetscNS::MatrixOperationsNS
} // namespace MoReFEM::TestNS::PetscNS::MatrixOperationsNS
#endif // MOREFEM_x_TEST_x_THIRD_PARTY_x_P_E_T_SC_x_VECTOR_I_O_x_INPUT_DATA_HPP_
#endif // MOREFEM_x_TEST_x_THIRD_PARTY_x_P_E_T_SC_x_MATRIX_OPERATIONS_x_INPUT_DATA_HPP_
......@@ -8,8 +8,8 @@
*/
#ifndef MOREFEM_x_TEST_x_THIRD_PARTY_x_P_E_T_SC_x_VECTOR_I_O_x_MODEL_HPP_
# define MOREFEM_x_TEST_x_THIRD_PARTY_x_P_E_T_SC_x_VECTOR_I_O_x_MODEL_HPP_
#ifndef MOREFEM_x_TEST_x_THIRD_PARTY_x_P_E_T_SC_x_MATRIX_OPERATIONS_x_MODEL_HPP_
# define MOREFEM_x_TEST_x_THIRD_PARTY_x_P_E_T_SC_x_MATRIX_OPERATIONS_x_MODEL_HPP_
# include <memory>
# include <vector>
......@@ -125,10 +125,10 @@ namespace MoReFEM
//! Accessor to the matrix A used for tests.
const GlobalMatrix& GetMatrixA() const noexcept;
//! Accessor to the matrix B used for tests.
const GlobalMatrix& GetMatrixB() const noexcept;
//! Accessor to the matrix C used for tests.
const GlobalMatrix& GetMatrixC() const noexcept;
......@@ -156,14 +156,14 @@ namespace MoReFEM
//! Global matrix A which is used in the tests.
GlobalMatrix::unique_ptr matrix_a_ = nullptr;
//! Global matrix B which is used in the tests.
GlobalMatrix::unique_ptr matrix_b_ = nullptr;
//! Global matrix C which is used in the tests.
GlobalMatrix::unique_ptr matrix_c_ = nullptr;
};
......@@ -176,4 +176,4 @@ namespace MoReFEM
#include "Test/ThirdParty/PETSc/MatrixOperations/Model.hxx"
#endif // MOREFEM_x_TEST_x_THIRD_PARTY_x_P_E_T_SC_x_VECTOR_I_O_x_MODEL_HPP_
#endif // MOREFEM_x_TEST_x_THIRD_PARTY_x_P_E_T_SC_x_MATRIX_OPERATIONS_x_MODEL_HPP_
......@@ -8,8 +8,8 @@
*/
#ifndef MOREFEM_x_TEST_x_THIRD_PARTY_x_P_E_T_SC_x_VECTOR_I_O_x_MODEL_HXX_
# define MOREFEM_x_TEST_x_THIRD_PARTY_x_P_E_T_SC_x_VECTOR_I_O_x_MODEL_HXX_
#ifndef MOREFEM_x_TEST_x_THIRD_PARTY_x_P_E_T_SC_x_MATRIX_OPERATIONS_x_MODEL_HXX_
# define MOREFEM_x_TEST_x_THIRD_PARTY_x_P_E_T_SC_x_MATRIX_OPERATIONS_x_MODEL_HXX_
namespace MoReFEM
......@@ -43,27 +43,27 @@ namespace MoReFEM
return *vector_;
}
inline const GlobalMatrix& Model::GetMatrixA() const noexcept
{
assert(!(!matrix_a_));
return *matrix_a_;
}
inline const GlobalMatrix& Model::GetMatrixB() const noexcept
{
assert(!(!matrix_b_));
return *matrix_b_;
}
inline const GlobalMatrix& Model::GetMatrixC() const noexcept
{
assert(!(!matrix_c_));
return *matrix_c_;
}
} // namespace TestNS::PetscNS::MatrixOperationsNS
......@@ -71,4 +71,4 @@ namespace MoReFEM
} // namespace MoReFEM
#endif // MOREFEM_x_TEST_x_THIRD_PARTY_x_P_E_T_SC_x_VECTOR_I_O_x_MODEL_HXX_
#endif // MOREFEM_x_TEST_x_THIRD_PARTY_x_P_E_T_SC_x_MATRIX_OPERATIONS_x_MODEL_HXX_
......@@ -46,7 +46,7 @@ namespace MoReFEM
* };
* \endcode
*/
template<class DerivedT, std::size_t NlocalMatricesT>
template<class DerivedT, std::size_t NlocalMatricesT, class LocalMatrixT = LocalMatrix>
class LocalMatrixStorage
{
public:
......@@ -120,13 +120,13 @@ namespace MoReFEM
* \endinternal
*/
template<std::size_t IndexT>
LocalMatrix& GetLocalMatrix() const;
LocalMatrixT& GetLocalMatrix() const;
private:
//! Local matrices stored.
mutable std::array<LocalMatrix, NlocalMatricesT> matrix_list_;
mutable std::array<LocalMatrixT, NlocalMatricesT> matrix_list_;
};
......
......@@ -24,30 +24,38 @@ namespace MoReFEM
{
template<class DerivedT, std::size_t NlocalMatricesT>
constexpr std::size_t LocalMatrixStorage<DerivedT, NlocalMatricesT>::N()
template<class DerivedT, std::size_t NlocalMatricesT, class LocalMatrixT>
constexpr std::size_t LocalMatrixStorage<DerivedT, NlocalMatricesT, LocalMatrixT>::N()
{
return NlocalMatricesT;
}
template<class DerivedT, std::size_t NlocalMatricesT>
void LocalMatrixStorage<DerivedT, NlocalMatricesT>
template<class DerivedT, std::size_t NlocalMatricesT, class LocalMatrixT>
void LocalMatrixStorage<DerivedT, NlocalMatricesT, LocalMatrixT>
::InitLocalMatrixStorage(const std::array<std::pair<unsigned int, unsigned int>, NlocalMatricesT>& matrices_dimension)
{
for (std::size_t i = 0ul; i < NlocalMatricesT; ++i)
{
auto& matrix = matrix_list_[i];
matrix.Resize(static_cast<int>(matrices_dimension[i].first),
static_cast<int>(matrices_dimension[i].second));
matrix.Zero();
if constexpr(std::is_same<LocalMatrixT, LocalMatrix>())
{
matrix.Resize(static_cast<int>(matrices_dimension[i].first),
static_cast<int>(matrices_dimension[i].second));
matrix.Zero();
}
else
{
}
}
}
template<class DerivedT, std::size_t NlocalMatricesT>
template<class DerivedT, std::size_t NlocalMatricesT, class LocalMatrixT>
template<std::size_t IndexT>
inline LocalMatrix& LocalMatrixStorage<DerivedT, NlocalMatricesT>::GetLocalMatrix() const
inline LocalMatrixT& LocalMatrixStorage<DerivedT, NlocalMatricesT, LocalMatrixT>::GetLocalMatrix() const
{
static_assert(IndexT < NlocalMatricesT, "Check index is within bounds!");
return matrix_list_[IndexT];
......
......@@ -18,9 +18,11 @@
# include <memory>
# include <vector>
# include "ThirdParty/Wrappers/Petsc/Vector/Vector.hpp"
# include "ThirdParty/Wrappers/Petsc/Matrix/Matrix.hpp"
# include "ThirdParty/Wrappers/Seldon/SubVector.hpp"
# include "ThirdParty/IncludeWithoutWarning/Xtensor/Xtensor.hpp"
namespace MoReFEM
......@@ -33,6 +35,8 @@ namespace MoReFEM
//! Convenient alias for the type used most of the time in MoReFEM for local matrices.
using LocalMatrix = ::Seldon::Matrix<double, Seldon::General, Seldon::RowMajor, Seldon::MallocAlloc<double>>;
using XtensorMatrix = xt::xtensor<double, 2, xt::layout_type::row_major>;
//! Convenient alias for the local symmetrix matrices, used in some specific operators.
using LocalSymmetricMatrix = ::Seldon::Matrix<double, Seldon::General, Seldon::RowSym, Seldon::MallocAlloc<double>>;
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment