diff --git a/ProcessLib/AssemblyMixin.cpp b/ProcessLib/AssemblyMixin.cpp index 2b955a2364459b3e3e100eacb14fdb2a77d4232d..bca2c4a68599056432ed315d1cb17331d1c2f0e7 100644 --- a/ProcessLib/AssemblyMixin.cpp +++ b/ProcessLib/AssemblyMixin.cpp @@ -10,6 +10,10 @@ #include "AssemblyMixin.h" +#include <range/v3/view/enumerate.hpp> +#include <range/v3/view/for_each.hpp> +#include <range/v3/view/transform.hpp> +#include <range/v3/view/zip.hpp> #include <unordered_set> #include "MeshLib/Utils/getOrCreateMeshProperty.h" @@ -18,30 +22,79 @@ namespace { -std::vector<std::reference_wrapper<MeshLib::PropertyVector<double>>> -createResiduumVectors( - MeshLib::Mesh& mesh, - std::vector<std::string> const& residuum_names, - std::vector<std::reference_wrapper<ProcessLib::ProcessVariable>> - pvs) -{ - auto const num_residua = residuum_names.size(); - std::vector<std::reference_wrapper<MeshLib::PropertyVector<double>>> - residuum_vectors; - residuum_vectors.reserve(num_residua); +void checkResiduumNamesVsProcessVariables( + std::vector<std::vector<std::string>> const& per_process_residuum_names, + std::vector< + std::vector<std::reference_wrapper<ProcessLib::ProcessVariable>>> const& + per_process_pvs) +{ + if (per_process_pvs.size() != per_process_residuum_names.size()) + { + OGS_FATAL( + "The number of passed residuum names ({}) does not match the " + "number of processes ({}).", + per_process_residuum_names.size(), per_process_pvs.size()); + } - for (std::size_t i = 0; i < num_residua; ++i) + auto check_sizes = [](std::size_t const process_id, auto const& rns_pvs) { - auto const& residuum_name = residuum_names[i]; - auto const& pv = pvs[i].get(); + auto const& [rns, pvs] = rns_pvs; - residuum_vectors.emplace_back(*MeshLib::getOrCreateMeshProperty<double>( - mesh, residuum_name, MeshLib::MeshItemType::Node, - pv.getNumberOfGlobalComponents())); + if (rns.size() != pvs.size()) + { + OGS_FATAL( + "The number of passed residuum names ({}) does not match the " + "number of process variables ({}) for process {}.", + rns.size(), pvs.size(), process_id); + } + }; + for (auto const& [process_id, rns_pvs] : + ranges::views::zip(per_process_residuum_names, per_process_pvs) | + ranges::views::enumerate) + { + check_sizes(process_id, rns_pvs); } +} - return residuum_vectors; +std::vector< + std::vector<std::reference_wrapper<MeshLib::PropertyVector<double>>>> +createResiduumVectors( + MeshLib::Mesh& mesh, + std::vector<std::vector<std::string>> const& per_process_residuum_names, + std::vector< + std::vector<std::reference_wrapper<ProcessLib::ProcessVariable>>> const& + per_process_pvs) +{ + checkResiduumNamesVsProcessVariables(per_process_residuum_names, + per_process_pvs); + + auto create_mesh_property_for_residuum = + [&mesh](std::pair<std::string const&, + ProcessLib::ProcessVariable const&> const& + residuum_name_process_variable) + -> std::reference_wrapper<MeshLib::PropertyVector<double>> + { + auto const& [rn, pv] = residuum_name_process_variable; + return *MeshLib::getOrCreateMeshProperty<double>( + mesh, rn, MeshLib::MeshItemType::Node, + pv.getNumberOfGlobalComponents()); + }; + + auto create_mesh_properties = + [&create_mesh_property_for_residuum](auto const& rns_pvs) + { + auto const& [rns, pvs] = rns_pvs; + return ranges::views::zip(rns, pvs) | + ranges::views::transform(create_mesh_property_for_residuum) | + ranges::to<std::vector< + std::reference_wrapper<MeshLib::PropertyVector<double>>>>; + }; + + return ranges::views::zip(per_process_residuum_names, per_process_pvs) | + ranges::views::transform(create_mesh_properties) | + ranges::to<std::vector<std::vector< + std::reference_wrapper<MeshLib::PropertyVector<double>>>>>; } } // namespace @@ -49,7 +102,8 @@ namespace ProcessLib { SubmeshAssemblyData::SubmeshAssemblyData( MeshLib::Mesh const& mesh, - std::vector<std::reference_wrapper<MeshLib::PropertyVector<double>>>&& + std::vector< + std::vector<std::reference_wrapper<MeshLib::PropertyVector<double>>>>&& residuum_vectors) : bulk_element_ids{*MeshLib::bulkElementIDs(mesh)}, bulk_node_ids{*MeshLib::bulkNodeIDs(mesh)}, @@ -58,24 +112,14 @@ SubmeshAssemblyData::SubmeshAssemblyData( } void AssemblyMixinBase::initializeAssemblyOnSubmeshes( - const int process_id, MeshLib::Mesh& bulk_mesh, std::vector<std::reference_wrapper<MeshLib::Mesh>> const& submeshes, - std::vector<std::string> const& residuum_names, - std::vector<std::reference_wrapper<ProcessVariable>> const& pvs) + std::vector<std::vector<std::string>> const& residuum_names, + std::vector<std::vector<std::reference_wrapper<ProcessVariable>>> const& + pvs) { DBUG("AssemblyMixinBase initializeSubmeshOutput()."); - auto const num_residua = residuum_names.size(); - - if (pvs.size() != num_residua) - { - OGS_FATAL( - "The number of passed residuum names does not match the number " - "of process variables for process id {}: {} != {}", - process_id, num_residua, pvs.size()); - } - submesh_assembly_data_.reserve(submeshes.size()); for (auto& mesh_ref : submeshes) { @@ -184,16 +228,18 @@ void AssemblyMixinBase::copyResiduumVectorsToBulkMesh( } void AssemblyMixinBase::copyResiduumVectorsToSubmesh( + int const process_id, GlobalVector const& rhs, NumLib::LocalToGlobalIndexMap const& local_to_global_index_map, SubmeshAssemblyData const& sad) { - for (std::size_t variable_id = 0; variable_id < sad.residuum_vectors.size(); + auto const& residuum_vectors = sad.residuum_vectors[process_id]; + for (std::size_t variable_id = 0; variable_id < residuum_vectors.size(); ++variable_id) { transformVariableFromGlobalVector( rhs, variable_id, local_to_global_index_map, - sad.residuum_vectors[variable_id].get(), sad.bulk_node_ids, + residuum_vectors[variable_id].get(), sad.bulk_node_ids, std::negate<double>()); } } diff --git a/ProcessLib/AssemblyMixin.h b/ProcessLib/AssemblyMixin.h index c68fb644e36603d067886217d60964c5e7b84279..6b4bc09511fc86d3da2973553ad7f3104e82efe8 100644 --- a/ProcessLib/AssemblyMixin.h +++ b/ProcessLib/AssemblyMixin.h @@ -24,13 +24,15 @@ struct SubmeshAssemblyData { explicit SubmeshAssemblyData( MeshLib::Mesh const& mesh, - std::vector<std::reference_wrapper<MeshLib::PropertyVector<double>>>&& + std::vector<std::vector< + std::reference_wrapper<MeshLib::PropertyVector<double>>>>&& residuum_vectors); MeshLib::PropertyVector<std::size_t> const& bulk_element_ids; MeshLib::PropertyVector<std::size_t> const& bulk_node_ids; std::vector<std::size_t> active_element_ids; - std::vector<std::reference_wrapper<MeshLib::PropertyVector<double>>> + std::vector< + std::vector<std::reference_wrapper<MeshLib::PropertyVector<double>>>> residuum_vectors; }; @@ -50,11 +52,11 @@ protected: : pvma_{jacobian_assembler} {}; void initializeAssemblyOnSubmeshes( - const int process_id, MeshLib::Mesh& bulk_mesh, std::vector<std::reference_wrapper<MeshLib::Mesh>> const& submeshes, - std::vector<std::string> const& residuum_names, - std::vector<std::reference_wrapper<ProcessVariable>> const& pvs); + std::vector<std::vector<std::string>> const& residuum_names, + std::vector<std::vector<std::reference_wrapper<ProcessVariable>>> const& + pvs); void updateActiveElements(ProcessLib::Process const& process); @@ -65,6 +67,7 @@ protected: residuum_vectors); static void copyResiduumVectorsToSubmesh( + int const process_id, GlobalVector const& rhs, NumLib::LocalToGlobalIndexMap const& local_to_global_index_map, SubmeshAssemblyData const& sad); @@ -74,7 +77,8 @@ private: protected: std::vector<SubmeshAssemblyData> submesh_assembly_data_; - std::vector<std::reference_wrapper<MeshLib::PropertyVector<double>>> + std::vector< + std::vector<std::reference_wrapper<MeshLib::PropertyVector<double>>>> residuum_vectors_bulk_; /// ID of the b vector on submeshes, cf. NumLib::VectorProvider. @@ -122,13 +126,12 @@ public: * and staggered coupling. */ void initializeAssemblyOnSubmeshes( - const int process_id, std::vector<std::reference_wrapper<MeshLib::Mesh>> const& submeshes, - std::vector<std::string> const& residuum_names) + std::vector<std::vector<std::string>> const& residuum_names) { AssemblyMixinBase::initializeAssemblyOnSubmeshes( - process_id, derived().getMesh(), submeshes, residuum_names, - derived().getProcessVariables(process_id)); + derived().getMesh(), submeshes, residuum_names, + derived().getProcessVariables()); } void updateActiveElements() @@ -182,7 +185,7 @@ public: MathLib::LinAlg::axpy(b, 1.0, b_submesh); AssemblyMixinBase::copyResiduumVectorsToSubmesh( - b_submesh, *(dof_tables.front()), sad); + process_id, b_submesh, *(dof_tables[process_id]), sad); } NumLib::GlobalVectorProvider::provider.releaseVector(b_submesh); @@ -195,7 +198,7 @@ public: } AssemblyMixinBase::copyResiduumVectorsToBulkMesh( - b, *(dof_tables.front()), residuum_vectors_bulk_); + b, *(dof_tables[process_id]), residuum_vectors_bulk_[process_id]); } private: diff --git a/ProcessLib/CreateTimeLoop.cpp b/ProcessLib/CreateTimeLoop.cpp index 0b22415787a7e8e2df50f860c3480879ea5d8b83..49b5bc64e85b28a65b6cbc82b5b5ee24add67179 100644 --- a/ProcessLib/CreateTimeLoop.cpp +++ b/ProcessLib/CreateTimeLoop.cpp @@ -12,6 +12,7 @@ #include <algorithm> #include <range/v3/algorithm/any_of.hpp> +#include <range/v3/view/join.hpp> #include <string> #include "BaseLib/ConfigTree.h" @@ -61,7 +62,7 @@ std::unique_ptr<TimeLoop> createTimeLoop( auto const& residuum_vector_names = process->initializeAssemblyOnSubmeshes(smroc.meshes); - for (auto const& name : residuum_vector_names) + for (auto const& name : residuum_vector_names | ranges::views::join) { smroc.output.doNotProjectFromBulkMeshToSubmeshes( name, MeshLib::MeshItemType::Node); diff --git a/ProcessLib/Process.h b/ProcessLib/Process.h index b8c4fd9ed64b77eb76548a8296a3b2392bd1f436..954a71821f086a5e45d265f86eff1bf23926886d 100644 --- a/ProcessLib/Process.h +++ b/ProcessLib/Process.h @@ -151,6 +151,13 @@ public: } MeshLib::Mesh& getMesh() const { return _mesh; } + + std::vector<std::vector<std::reference_wrapper<ProcessVariable>>> const& + getProcessVariables() const + { + return _process_variables; + } + std::vector<std::reference_wrapper<ProcessVariable>> const& getProcessVariables(const int process_id) const { diff --git a/ProcessLib/SubmeshAssemblySupport.h b/ProcessLib/SubmeshAssemblySupport.h index 669b3029c3418b7a57173d0119094ca9e4c0b719..e188674ddf661dd87153663c52432ad6ca3201dc 100644 --- a/ProcessLib/SubmeshAssemblySupport.h +++ b/ProcessLib/SubmeshAssemblySupport.h @@ -32,8 +32,10 @@ public: /// \attention \c meshes must be a must be a non-overlapping cover of the /// entire simulation domain (bulk mesh)! /// - /// \return The names of the residuum vectors that will be assembled. - virtual std::vector<std::string> initializeAssemblyOnSubmeshes( + /// \return The names of the residuum vectors that will be assembled for + /// each process: outer vector of size 1 for monolithic schemes and greater + /// for staggered schemes. + virtual std::vector<std::vector<std::string>> initializeAssemblyOnSubmeshes( std::vector<std::reference_wrapper<MeshLib::Mesh>> const& meshes) { DBUG( diff --git a/ProcessLib/TH2M/TH2MProcess.cpp b/ProcessLib/TH2M/TH2MProcess.cpp index c7f66f32522cb681da24cf5ce46d044c069c9b1a..3f7e384cadeb18af92e19f754b698e8eeaf29ba9 100644 --- a/ProcessLib/TH2M/TH2MProcess.cpp +++ b/ProcessLib/TH2M/TH2MProcess.cpp @@ -316,17 +316,17 @@ void TH2MProcess<DisplacementDim>::computeSecondaryVariableConcrete( } template <int DisplacementDim> -std::vector<std::string> +std::vector<std::vector<std::string>> TH2MProcess<DisplacementDim>::initializeAssemblyOnSubmeshes( std::vector<std::reference_wrapper<MeshLib::Mesh>> const& meshes) { INFO("TH2M process initializeSubmeshOutput()."); - const int process_id = 0; - std::vector<std::string> residuum_names{ - "GasMassFlowRate", "LiquidMassFlowRate", "HeatFlowRate", "NodalForces"}; + std::vector<std::vector<std::string>> residuum_names{ + {"GasMassFlowRate", "LiquidMassFlowRate", "HeatFlowRate", + "NodalForces"}}; AssemblyMixin<TH2MProcess<DisplacementDim>>::initializeAssemblyOnSubmeshes( - process_id, meshes, residuum_names); + meshes, residuum_names); return residuum_names; } diff --git a/ProcessLib/TH2M/TH2MProcess.h b/ProcessLib/TH2M/TH2MProcess.h index 6b8942efd4411dcd75074fe61452df056788ba4a..96ff5ddb1a2e36b10ff9bc015f13c40c60486682 100644 --- a/ProcessLib/TH2M/TH2MProcess.h +++ b/ProcessLib/TH2M/TH2MProcess.h @@ -105,7 +105,7 @@ private: NumLib::LocalToGlobalIndexMap const& getDOFTable( const int process_id) const override; - std::vector<std::string> initializeAssemblyOnSubmeshes( + std::vector<std::vector<std::string>> initializeAssemblyOnSubmeshes( std::vector<std::reference_wrapper<MeshLib::Mesh>> const& meshes) override; diff --git a/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsProcess.cpp b/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsProcess.cpp index 3813235dd1bfc89445f697f364d6ab01ef2ed13a..c763ac46da494aa8d16c2d3e152c0d49e21603a2 100644 --- a/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsProcess.cpp +++ b/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsProcess.cpp @@ -239,19 +239,18 @@ void ThermoRichardsMechanicsProcess<DisplacementDim, ConstitutiveTraits>:: } template <int DisplacementDim, typename ConstitutiveTraits> -std::vector<std::string> +std::vector<std::vector<std::string>> ThermoRichardsMechanicsProcess<DisplacementDim, ConstitutiveTraits>:: initializeAssemblyOnSubmeshes( std::vector<std::reference_wrapper<MeshLib::Mesh>> const& meshes) { INFO("TRM process initializeSubmeshOutput()."); - const int process_id = 0; - std::vector<std::string> residuum_names{"HeatFlowRate", "MassFlowRate", - "NodalForces"}; + std::vector<std::vector<std::string>> residuum_names{ + {"HeatFlowRate", "MassFlowRate", "NodalForces"}}; AssemblyMixin< ThermoRichardsMechanicsProcess<DisplacementDim, ConstitutiveTraits>>:: - initializeAssemblyOnSubmeshes(process_id, meshes, residuum_names); + initializeAssemblyOnSubmeshes(meshes, residuum_names); return residuum_names; } diff --git a/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsProcess.h b/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsProcess.h index a13c2b0d33896f9cb9a390e6d55d541bc53e9925..6ae284ce890c20b636d34beadc8eb81c730c6307 100644 --- a/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsProcess.h +++ b/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsProcess.h @@ -199,7 +199,7 @@ private: NumLib::LocalToGlobalIndexMap const& getDOFTable( const int process_id) const override; - std::vector<std::string> initializeAssemblyOnSubmeshes( + std::vector<std::vector<std::string>> initializeAssemblyOnSubmeshes( std::vector<std::reference_wrapper<MeshLib::Mesh>> const& meshes) override;