Commit 0200cd3f authored by GILLES Sebastien's avatar GILLES Sebastien
Browse files

#1443 Modify the code to make the test concerning matrix pattern work (in the...

#1443 Modify the code to make the test concerning matrix pattern work (in the end we had to keep it, as NodeBearer that are neither processor-wise nor ghost ones are involved in its computation).
parent 5950e8cd
......@@ -4207,6 +4207,7 @@
BE93B31C1F0536DA004F84CF /* OutputDeformedMesh.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = OutputDeformedMesh.cpp; sourceTree = "<group>"; };
BE93B31D1F0536DA004F84CF /* OutputDeformedMesh.hpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.h; path = OutputDeformedMesh.hpp; sourceTree = "<group>"; };
BE93B31E1F0536DA004F84CF /* OutputDeformedMesh.hxx */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.h; path = OutputDeformedMesh.hxx; sourceTree = "<group>"; };
BE940F672481390D0095C5BB /* demo_initial.lua */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = text; path = demo_initial.lua; sourceTree = "<group>"; };
BE9453992100BCA500953B84 /* CMakeLists.txt */ = {isa = PBXFileReference; lastKnownFileType = text; path = CMakeLists.txt; sourceTree = "<group>"; };
BE94539B2100BCA500953B84 /* CMakeLists.txt */ = {isa = PBXFileReference; lastKnownFileType = text; path = CMakeLists.txt; sourceTree = "<group>"; };
BE94539C2100BCA500953B84 /* demo_morefem_data.lua */ = {isa = PBXFileReference; lastKnownFileType = text; path = demo_morefem_data.lua; sourceTree = "<group>"; };
......@@ -5562,6 +5563,7 @@
137F2B081E38B22900BD6083 /* Hyperelasticity */ = {
isa = PBXGroup;
children = (
BE940F672481390D0095C5BB /* demo_initial.lua */,
BE0AFB7220750C570089FD9D /* CMakeLists.txt */,
BEB0C6171EC3054500D62905 /* README */,
BE6EA6201EC330820085E651 /* demo_input_hyperelasticity.lua */,
......@@ -75,6 +75,10 @@
argument = "-i ${HOME}/Codes/MoReFEM/CoreLibrary/Sources/ModelInstances/Elasticity/demo_2d.lua"
isEnabled = "YES">
</CommandLineArgument>
<CommandLineArgument
argument = "--overwrite_directory"
isEnabled = "YES">
</CommandLineArgument>
</CommandLineArguments>
</LaunchAction>
<ProfileAction
......
......@@ -85,11 +85,11 @@
</CommandLineArgument>
<CommandLineArgument
argument = "-i ${HOME}/Codes/MoReFEM/CoreLibrary/Sources/ModelInstances/Hyperelasticity/demo_initial.lua"
isEnabled = "NO">
isEnabled = "YES">
</CommandLineArgument>
<CommandLineArgument
argument = "-i ${HOME}/Codes/MoReFEM/CoreLibrary/Sources/ModelInstances/Hyperelasticity/demo_rerun.lua"
isEnabled = "YES">
isEnabled = "NO">
</CommandLineArgument>
</CommandLineArguments>
</LaunchAction>
......
......@@ -83,8 +83,7 @@ namespace MoReFEM
{
return domain.IsGeometricEltInside(local_felt_space.GetGeometricElt());
}
void LocalFEltSpace
::InitLocal2Global(const NumberingSubset::vector_const_shared_ptr& numbering_subset_list,
......
......@@ -446,7 +446,6 @@ namespace MoReFEM
}
void GodOfDof::OndomaticLikePrint(const NumberingSubset& numbering_subset,
const std::string& output_directory) const
{
......
......@@ -70,7 +70,8 @@ namespace MoReFEM
col_numbering_subset_ptr->GetUniqueId());
auto it = std::find(already_seen.cbegin(), already_seen.cend(), new_pair);
assert(it == already_seen.cend() && "One given numbering subset should be only once in the list!");
assert(it == already_seen.cend()
&& "One given numbering subset should be only once in the list!");
already_seen.push_back(new_pair);
#endif // NDEBUG
......
......@@ -87,15 +87,17 @@ namespace MoReFEM
public:
/*!
* \brief Compute the pattern(s) that have to be covered within a GodOfDof.
* \brief Compute the pattern(s) that have to be covered within a \a GodOfDof.
*
* \param[in] mpi_rank Rank of current processor.
* \param[in] felt_space_list The list of finite elements spaces considered.
* \param[in] node_bearer_list The list of node bearers. In output, some informations
* have been provided, such as the processor and the index of each node bearer.
* \param[in,out] node_bearer_list The list of node bearers that must encompass all the ones relevant for the current processor,
* including the ghost ones). This argument is tagged as an output parameter because the \a NodeBearer objects in the list gets new
* informations in the process.
* \param[in] Ndof_holder Object that knows the number of dofs per numbering subset.
* \param[in] numbering_subset_list List of \a NumberingSubset to consider.
*
*
* \return Pattern of the global matrix.
*/
static MatrixPattern::vector_const_unique_ptr
......
......@@ -35,7 +35,7 @@ namespace MoReFEM
void ComputeNodeBearerConnectivityInFEltSpace(const FEltSpace& felt_space,
connectivity_type& connectivity,
const unsigned int mpi_rank);
void CleanUpConnectivity(connectivity_type& connectivity,
KeepSelfConnexion do_keep_self_connexion);
......@@ -94,36 +94,37 @@ namespace MoReFEM
// Iterate through all node_bearers and determine the list of connected node_bearers.
// In this first approach, a same node_bearer might be listed twice in the list.
for (const auto& pair : felt_list_per_type)
for (const auto& [ref_local_felt_space, local_felt_space_list_per_geom_index ] : felt_list_per_type)
{
const auto& local_felt_space_list = pair.second;
for (const auto& local_felt_space_pair : local_felt_space_list)
static_cast<void>(ref_local_felt_space);
for (const auto& [geom_elt_index, local_felt_space_ptr] : local_felt_space_list_per_geom_index)
{
const auto& local_felt_space_ptr = local_felt_space_pair.second;
static_cast<void>(geom_elt_index);
assert(!(!local_felt_space_ptr));
const auto& node_bearers_in_current_finite_element =
local_felt_space_ptr->GetNodeBearerList();
for (const auto& node_bearer_ptr : node_bearers_in_current_finite_element)
const auto& local_felt_space = *local_felt_space_ptr;
const auto& node_bearers_in_current_local_felt_space =
local_felt_space.GetNodeBearerList();
decltype(auto) geometric_element = local_felt_space.GetGeometricElt();
assert(geometric_element.GetIndex() == geom_elt_index);
for (const auto& node_bearer_ptr : node_bearers_in_current_local_felt_space)
{
if (!is_on_local_processor(node_bearer_ptr))
continue;
// The node_bearer might already exist; check that here!
auto it = connectivity.find(node_bearer_ptr);
if (it == connectivity.cend())
connectivity.insert({node_bearer_ptr, node_bearers_in_current_finite_element});
else
{
// If already existing, add the new node_bearers to the existing vector.
auto& second_member = it->second;
for (const auto& connected_node_bearer_ptr : node_bearers_in_current_finite_element)
second_member.push_back(connected_node_bearer_ptr);
}
// Use of C++ 17 "structure binding".
// The following code works whether the \a node_bearer_ptr is already in the
// unordered_map or not.
auto [it, is_new_element] = connectivity.insert({node_bearer_ptr, {}});
static_cast<void>(is_new_element);
auto& second_member = it->second;
std::copy(node_bearers_in_current_local_felt_space.cbegin(),
node_bearers_in_current_local_felt_space.cend(),
std::back_inserter(second_member));
}
}
}
......
......@@ -143,7 +143,7 @@ namespace MoReFEM
ghost_dof_list,
it,
end);
assert(it == end && "The data in the prepartition data file should cover exactly the number of "
"reconstructed dofs in the numbering subset.");
}
......
......@@ -32,7 +32,9 @@ namespace MoReFEM
const NumberingSubset::const_shared_ptr& column_numbering_subset_ptr,
const connectivity_type& connectivity,
const NodeBearer::vector_shared_ptr& processor_wise_node_bearer_list,
const NdofHolder& Ndof_holder)
const NdofHolder& Ndof_holder,
std::optional<std::reference_wrapper<NodeBearer::vector_shared_ptr>> ghost_like_node_bearer_list,
bool do_print)
: row_numbering_subset_(row_numbering_subset_ptr),
column_numbering_subset_(column_numbering_subset_ptr)
{
......@@ -94,6 +96,9 @@ namespace MoReFEM
if (!node.DoBelongToNumberingSubset(column_numbering_subset))
continue;
if (ghost_like_node_bearer_list.has_value())
ghost_like_node_bearer_list.value().get().push_back(connected_node_bearer_ptr);
const auto& dof_list = node.GetDofList();
......@@ -140,9 +145,69 @@ namespace MoReFEM
matrix_pattern_ =
std::make_unique<::MoReFEM::Wrappers::Petsc::MatrixPattern>(content_for_each_local_row);
if (ghost_like_node_bearer_list.has_value())
{
auto& vector = ghost_like_node_bearer_list.value().get();
Utilities::EliminateDuplicate(vector,
Utilities::PointerComparison::Less<NodeBearer::shared_ptr>(),
Utilities::PointerComparison::Equal<NodeBearer::shared_ptr>());
NodeBearer::vector_shared_ptr ret;
std::set_difference(vector.cbegin(),
vector.cend(),
processor_wise_node_bearer_list.cbegin(),
processor_wise_node_bearer_list.cend(),
std::back_inserter(ret),
Utilities::PointerComparison::Less<NodeBearer::shared_ptr>());
#ifndef NDEBUG
{
if (!processor_wise_node_bearer_list.empty())
{
// std::cout << "SIZE = " << processor_wise_node_bearer_list.size() << std::endl;
decltype(auto) first_node_bearer_ptr = *processor_wise_node_bearer_list.cbegin();
assert(!(!first_node_bearer_ptr));
const auto mpi_rank = first_node_bearer_ptr->GetProcessor();
auto belong_to_rank = [mpi_rank](const auto& node_bearer_ptr)
{
assert(!(!node_bearer_ptr));
return node_bearer_ptr->GetProcessor() == mpi_rank;
};
assert(std::all_of(processor_wise_node_bearer_list.cbegin(),
processor_wise_node_bearer_list.cend(),
belong_to_rank));
assert(std::none_of(ret.cbegin(),
ret.cend(),
belong_to_rank));
}
}
#endif // NDEBUG
std::swap(ret, vector);
//std::cout << "MATRIX PATTERN Nghost_like = " << ret.size() << std::endl;
}
}
MatrixPattern::MatrixPattern(const NumberingSubset::const_shared_ptr& row_numbering_subset,
const NumberingSubset::const_shared_ptr& column_numbering_subset,
std::vector<PetscInt>&& iCSR,
std::vector<PetscInt>&& jCSR)
: row_numbering_subset_(row_numbering_subset),
column_numbering_subset_(column_numbering_subset),
matrix_pattern_(std::make_unique<const ::MoReFEM::Wrappers::Petsc::MatrixPattern>(std::move(iCSR),
std::move(jCSR)))
{ }
bool operator==(const MatrixPattern& lhs, const MatrixPattern& rhs)
{
if (lhs.GetRowNumberingSubset() != rhs.GetRowNumberingSubset())
......
......@@ -111,7 +111,28 @@ namespace MoReFEM
const NumberingSubset::const_shared_ptr& column_numbering_subset,
const connectivity_type& connectivity,
const NodeBearer::vector_shared_ptr& processor_wise_node_bearer_list,
const NdofHolder& Ndof_holder);
const NdofHolder& Ndof_holder,
std::optional<std::reference_wrapper<NodeBearer::vector_shared_ptr>> ghost_like_node_bearer_list = std::nullopt,
bool do_print = false); // TMP #1443
/*!
* \brief Computes the number of diagonal and off-diagonal terms for each row of a dof matrix.
*
*
* Beware: The indexes of the row are processor-wise, whereas the values for each vector are program-wise!
* (not my fault: it is the way Petsc is designed...)
*
* \param[in] row_numbering_subset Numbering subset over which all dofs considered in the rows of the
* matrix are defined.
* \param[in] column_numbering_subset Numbering subset over which all dofs considered in the rows of the
* matrix are defined.
* \param[in] iCSR i-part of the CSR pattern.
* \param[in] jCSR j-part of the CSR pattern.
*/
explicit MatrixPattern(const NumberingSubset::const_shared_ptr& row_numbering_subset,
const NumberingSubset::const_shared_ptr& column_numbering_subset,
std::vector<PetscInt>&& iCSR,
std::vector<PetscInt>&& jCSR);
//! Destructor.
~MatrixPattern() = default;
......
......@@ -68,7 +68,7 @@ namespace MoReFEM
void ComputeGhostNodeBearerList(const unsigned int mpi_rank,
const LocalFEltSpacePerRefLocalFEltSpace& processor_wise_felt_list_per_type,
NodeBearer::vector_shared_ptr& ghost_node_bearer_list);
} // namespace anonymous
......@@ -106,10 +106,25 @@ namespace MoReFEM
Mesh& mesh)
{
// if (mpi.GetRank<int>() == 1)
// {
// std::cout << "//////////// ReduceToProcessorWise::Perform ///////////////" << std::endl;
//
// auto fii = std::find_if(processor_wise_node_bearer_list.cbegin(),
// processor_wise_node_bearer_list.cend(),
// [](const auto& ptr)
// {
// return ptr->GetNature() == ::MoReFEM::InterfaceNS::Nature::vertex
// && ptr->GetInterface().GetVertexCoordsIndexList().back() == 18;
// });
//
// if (fii != processor_wise_node_bearer_list.cend())
// std::cout << "\t COORDS 18 STILL FOUND in ReduceToProcessorWise::Perform() ON PROC 1!" << std::endl;
//
// }
const unsigned int mpi_rank = mpi.GetRank<unsigned int>();
assert(ghost_node_bearer_list.empty());
std::unordered_map<unsigned int, unsigned int> processor_for_each_geom_elt;
processor_for_each_geom_elt.max_load_factor(Utilities::DefaultMaxLoadFactor());
......@@ -131,12 +146,44 @@ namespace MoReFEM
Utilities::EliminateDuplicate(ghost_node_bearer_list,
Utilities::PointerComparison::Less<NodeBearer::shared_ptr>(),
Utilities::PointerComparison::Equal<NodeBearer::shared_ptr>());
// std::cout << "IS SUBSET? = " << std::includes(ancillary_ghost_node_bearer_list.cbegin(),
// ancillary_ghost_node_bearer_list.cend(),
// ghost_node_bearer_list.cbegin(),
// ghost_node_bearer_list.cend(),
// Utilities::PointerComparison::Less<NodeBearer::shared_ptr>()) << std::endl;
#ifndef NDEBUG
{
assert(std::is_sorted(processor_wise_node_bearer_list.cbegin(),
processor_wise_node_bearer_list.cend(),
Utilities::PointerComparison::Less<NodeBearer::shared_ptr>()));
assert(std::is_sorted(ghost_node_bearer_list.cbegin(),
ghost_node_bearer_list.cend(),
Utilities::PointerComparison::Less<NodeBearer::shared_ptr>()));
NodeBearer::vector_shared_ptr intersection;
std::set_intersection(processor_wise_node_bearer_list.cbegin(),
processor_wise_node_bearer_list.cend(),
ghost_node_bearer_list.cbegin(),
ghost_node_bearer_list.cend(),
std::back_inserter(intersection),
Utilities::PointerComparison::Less<NodeBearer::shared_ptr>());
assert(intersection.empty());
}
#endif // NDEBUG
ReduceMesh::Perform(mpi,
felt_space_list,
processor_wise_node_bearer_list,
ghost_node_bearer_list,
mesh);
// if (mpi.GetRank<int>() == 1)
// std::cout << "//////////// ENDReduceToProcessorWise::Perform ///////////////" << std::endl;
}
......@@ -205,7 +252,12 @@ namespace MoReFEM
CoordsListHelper(processor_wise_node_bearer_list, processor_wise_coords_list);
CoordsListHelper(ghost_node_bearer_list, ghosted_coords_list);
// std::cout << "NGHOST COORDS LIST = " << ghosted_coords_list.size() << std::endl;
// LOOK FOR NODE BEARR 109 HERE!
// It is more than weith than Nghost at point 3 is always 0...
//
Utilities::EliminateDuplicate(processor_wise_coords_list,
Utilities::PointerComparison::Less<Coords*>(),
Utilities::PointerComparison::Equal<Coords*>());
......@@ -213,6 +265,8 @@ namespace MoReFEM
Utilities::EliminateDuplicate(ghosted_coords_list,
Utilities::PointerComparison::Less<Coords*>(),
Utilities::PointerComparison::Equal<Coords*>());
// std::cout << "2 - NGHOST COORDS LIST = " << ghosted_coords_list.size() << std::endl;
const auto proc_wise_coords_list_begin = processor_wise_coords_list.cbegin();
const auto proc_wise_coords_list_end = processor_wise_coords_list.cend();
......@@ -227,8 +281,10 @@ namespace MoReFEM
coords,
Utilities::PointerComparison::Less<const Coords*>());
});
ghosted_coords_list.erase(it_partition, ghosted_coords_list.end());
// std::cout << "3 - NGHOST COORDS LIST = " << ghosted_coords_list.size() << std::endl;
std::vector<unsigned int> vertex_index_list;
std::vector<unsigned int> edge_index_list;
......@@ -331,17 +387,10 @@ namespace MoReFEM
// Determine to which processor the local finite element space should be attributed.
const auto chosen_processor_for_local_felt_space =
attribute_processor_helper.ProcessorForCurrentLocalFEltSpace(local_felt_space);
attribute_processor_helper.ProcessorForCurrentLocalFEltSpace(local_felt_space);
if (chosen_processor_for_local_felt_space == mpi_rank)
buf.insert({ local_felt_space.GetGeometricElt().GetIndex(), local_felt_space_ptr });
// IF NOT BUT IF THERE ARE SOME NODE BEARERS DOF SHOULD BE ADDED TO THE DOF LIST OF THE
// FINITE ELEMENT SPACE (WHERE IT COULD IN FACT BE GHOST OR NOT!)
// IF THIS DOES NOT WORK, another SOLUTION (perhaps even better!) IS TO GIVE YET ANOTHER INDEX TO DOFS, THAT
// ACTS AS A unique_id. We could then store the indexes in the reduction, not getting in the way
// of the reduction, and then retrieve any relevant dof.
}
#ifndef NDEBUG
......@@ -358,7 +407,7 @@ namespace MoReFEM
}
void
ComputeGhostNodeBearerList(const unsigned int mpi_rank,
const LocalFEltSpacePerRefLocalFEltSpace& processor_wise_felt_list_per_type,
......@@ -367,26 +416,26 @@ namespace MoReFEM
for (const auto& pair : processor_wise_felt_list_per_type)
{
const auto& local_felt_space_list = pair.second;
for (const auto& local_felt_space_pair : local_felt_space_list)
{
const auto& local_felt_space_ptr = local_felt_space_pair.second;
assert(!(!local_felt_space_ptr));
const auto& node_bearer_list_in_geometric_element = local_felt_space_ptr->GetNodeBearerList();
// All node_bearers that are not on the current processor must be ghosted here.
for (const auto& node_bearer_ptr : node_bearer_list_in_geometric_element)
{
assert(!(!node_bearer_ptr));
if (node_bearer_ptr->GetProcessor() != mpi_rank)
ghost_node_bearer_list.push_back(node_bearer_ptr);
}
}
}
}
} // namespace anonymous
......
......@@ -87,7 +87,6 @@ namespace MoReFEM
* processor-wise finite elements but not in processor-wise node bearer list. It can also be an input
* parameter: if several finite element spaces are involved the same list us kept for all of them.
*/
static void Perform(const ::MoReFEM::Wrappers::Mpi& mpi,
const FEltSpace::vector_unique_ptr& felt_space_list,
const NodeBearer::vector_shared_ptr& processor_wise_node_bearer_list,
......
......@@ -122,6 +122,40 @@ namespace MoReFEM::Internal::GodOfDofNS
out, ", ", "{ ", " }\n");
out << std::endl;
}
out << "\n-- For each numbering subset pair, the CSR pattern of the matrix." << std::endl;
out << "-- Unfortunately this can't be rebuild from local partitioned data: we need MORE NodeBearer "
"than those already covered by processor-wise and ghost, and they aren't otherwise needed in the code."
<< std::endl;
for (const auto& row_numbering_subset_ptr : numbering_subset_list)
{
assert(!(!row_numbering_subset_ptr));
const auto& row_numbering_subset = *row_numbering_subset_ptr;
for (const auto& col_numbering_subset_ptr : numbering_subset_list)
{
assert(!(!col_numbering_subset_ptr));
const auto& col_numbering_subset = *col_numbering_subset_ptr;
decltype(auto) matrix_pattern = god_of_dof.GetMatrixPattern(row_numbering_subset,
col_numbering_subset);
out << "iCSR_matrix_pattern_for_numbering_subset_pair_" << row_numbering_subset.GetUniqueId()
<< "_" << col_numbering_subset.GetUniqueId() << " = ";
Utilities::PrintContainer<>::Do(matrix_pattern.GetICsr(),
out, ", ", "{ ", " }\n");
out << "jCSR_matrix_pattern_for_numbering_subset_pair_" << row_numbering_subset.GetUniqueId()
<< "_" << col_numbering_subset.GetUniqueId() << " = ";
Utilities::PrintContainer<>::Do(matrix_pattern.GetJCsr(),
out, ", ", "{ ", " }\n");
out << std::endl;
}
}
out << std::endl;
}
}
......
......@@ -58,53 +58,6 @@ namespace MoReFEM
const ::MoReFEM::FilesystemNS::Directory& mesh_directory);
// class WritePrepartitionedData
// {
//
// public:
//
// //! \copydoc doxygen_hide_alias_self
// // \TODO This might seem a bit dumb but is actually very convenient for template classes.
// using self = WritePrepartitionedData;
//
// //! Alias to unique pointer.
// using unique_ptr = std::unique_ptr<self>;
//
// //! Alias to vector of unique pointers.
// using vector_unique_ptr = std::vector<unique_ptr>;
//
// public:
//
// /// \name Special members.
// ///@{
//
// //! Constructor.
// explicit WritePrepartitionedData() = default;
//
// //! Destructor.
// ~WritePrepartitionedData() = default;
//
// //! \copydoc doxygen_hide_copy_constructor
// WritePrepartitionedData(const WritePrepartitionedData& rhs) = delete;
//
// //! \copydoc doxygen_hide_move_constructor
// WritePrepartitionedData(WritePrepartitionedData&& rhs) = delete;
//
// //! \copydoc doxygen_hide_copy_affectation
// WritePrepartitionedData& operator=(const WritePrepartitionedData& rhs) = delete;
//
// //! \copydoc doxygen_hide_move_affectation
// WritePrepartitionedData& operator=(WritePrepartitionedData&& rhs) = delete;
//
// ///@}