diff --git a/NumLib/DOF/LocalToGlobalIndexMap.cpp b/NumLib/DOF/LocalToGlobalIndexMap.cpp
index 82ec698b7ca3a6635c7aec953ce872ea6749e3cf..ef7204152ed5fdb346f9d0039e3f1fbdeb91bac1 100644
--- a/NumLib/DOF/LocalToGlobalIndexMap.cpp
+++ b/NumLib/DOF/LocalToGlobalIndexMap.cpp
@@ -29,6 +29,46 @@ std::vector<T> to_cumulative(std::vector<T> const& vec)
} // no named namespace
+template <typename ElementIterator>
+void
+LocalToGlobalIndexMap::findGlobalIndicesWithElementID(
+ ElementIterator first, ElementIterator last,
+ std::vector<MeshLib::Node*> const& nodes, std::size_t const mesh_id,
+ const unsigned comp_id, const unsigned comp_id_write)
+{
+ // _rows should be resized based on an element ID
+ std::size_t max_elem_id = 0;
+ for (ElementIterator e = first; e != last; ++e)
+ max_elem_id = std::max(max_elem_id, (*e)->getID());
+ if (max_elem_id+1 > static_cast<unsigned>(_rows.rows()))
+ _rows.conservativeResize(max_elem_id + 1, _mesh_subsets.size());
+
+ std::vector<MeshLib::Node*> sorted_nodes{nodes};
+ std::sort(std::begin(sorted_nodes), std::end(sorted_nodes));
+
+ // For each element find the global indices for node/element
+ // components.
+ for (ElementIterator e = first; e != last; ++e)
+ {
+ LineIndex indices;
+
+ for (auto* n = (*e)->getNodes();
+ n < (*e)->getNodes()+(*e)->getNumberOfNodes(); ++n)
+ {
+ // Check if the element's node is in the given list of nodes.
+ if (!std::binary_search(std::begin(sorted_nodes),
+ std::end(sorted_nodes), *n))
+ continue;
+ MeshLib::Location l(
+ mesh_id, MeshLib::MeshItemType::Node, (*n)->getID());
+ indices.push_back(_mesh_component_map.getGlobalIndex(l, comp_id));
+ }
+
+ _rows((*e)->getID(), comp_id_write) = std::move(indices);
+ }
+}
+
+
template <typename ElementIterator>
void LocalToGlobalIndexMap::findGlobalIndices(
ElementIterator first, ElementIterator last,
@@ -109,6 +149,50 @@ LocalToGlobalIndexMap::LocalToGlobalIndexMap(
}
}
+
+LocalToGlobalIndexMap::LocalToGlobalIndexMap(
+ std::vector<std::unique_ptr<MeshLib::MeshSubsets>>&& mesh_subsets,
+ std::vector<unsigned> const& vec_var_n_components,
+ std::vector<std::vector<MeshLib::Element*>const*> const& vec_var_elements,
+ NumLib::ComponentOrder const order)
+ : _mesh_subsets(std::move(mesh_subsets)),
+ _mesh_component_map(_mesh_subsets, order),
+ _variable_component_offsets(to_cumulative(vec_var_n_components))
+{
+ assert(vec_var_n_components.size() == vec_var_elements.size());
+
+ // For all MeshSubsets and each of their MeshSubset's and each element
+ // of that MeshSubset save a line of global indices.
+
+
+ std::size_t offset = 0;
+ for (int variable_id = 0; variable_id < static_cast<int>(vec_var_n_components.size());
+ ++variable_id)
+ {
+ std::vector<MeshLib::Element*> const& var_elements = *vec_var_elements[variable_id];
+ for (int component_id = 0; component_id < static_cast<int>(vec_var_n_components[variable_id]);
+ ++component_id)
+ {
+ auto const global_component_id =
+ getGlobalComponent(variable_id, component_id);
+
+ auto const& mss = *_mesh_subsets[global_component_id];
+ for (int subset_id = 0; subset_id < static_cast<int>(mss.size());
+ ++subset_id)
+ {
+ MeshLib::MeshSubset const& ms = mss.getMeshSubset(subset_id);
+ std::size_t const mesh_id = ms.getMeshID();
+
+ findGlobalIndicesWithElementID(var_elements.cbegin(), var_elements.cend(), ms.getNodes(),
+ mesh_id, global_component_id, global_component_id);
+ }
+ // increase by number of components of that variable
+ offset += mss.size();
+ }
+ }
+}
+
+
LocalToGlobalIndexMap::LocalToGlobalIndexMap(
std::vector<std::unique_ptr<MeshLib::MeshSubsets>>&& mesh_subsets,
int const component_id,
diff --git a/NumLib/DOF/LocalToGlobalIndexMap.h b/NumLib/DOF/LocalToGlobalIndexMap.h
index 2ebda14202b8051f64ad208b2257b9ebcf65f738..379da7870f0ec11fee43a1e211aed4d0a52da53e 100644
--- a/NumLib/DOF/LocalToGlobalIndexMap.h
+++ b/NumLib/DOF/LocalToGlobalIndexMap.h
@@ -64,6 +64,21 @@ public:
std::vector<unsigned> const& vec_var_n_components,
NumLib::ComponentOrder const order);
+ /// Creates a MeshComponentMap internally and stores the global indices for
+ /// the given mesh elements
+ ///
+ /// \param mesh_subsets a vector of components
+ /// \param vec_var_n_components a vector of the number of variable components.
+ /// The size of the vector should be equal to the number of variables. Sum of the entries
+ /// should be equal to the size of the mesh_subsets.
+ /// \param vec_var_elements a vector of active mesh elements for each variable.
+ /// \param order type of ordering values in a vector
+ LocalToGlobalIndexMap(
+ std::vector<std::unique_ptr<MeshLib::MeshSubsets>>&& mesh_subsets,
+ std::vector<unsigned> const& vec_var_n_components,
+ std::vector<std::vector<MeshLib::Element*>const*> const& vec_var_elements,
+ NumLib::ComponentOrder const order);
+
/// Derive a LocalToGlobalIndexMap constrained to a set of mesh subsets and
/// elements. A new mesh component map will be constructed using the passed
/// mesh_subsets for the given variable and component ids.
@@ -158,6 +173,13 @@ private:
std::size_t const mesh_id,
const unsigned component_id, const unsigned comp_id_write);
+ template <typename ElementIterator>
+ void
+ findGlobalIndicesWithElementID(ElementIterator first, ElementIterator last,
+ std::vector<MeshLib::Node*> const& nodes,
+ std::size_t const mesh_id,
+ const unsigned component_id, const unsigned comp_id_write);
+
/// The global component id for the specific variable (like velocity) and a
/// component (like x, or y, or z).
std::size_t getGlobalComponent(int const variable_id,
diff --git a/Tests/NumLib/LocalToGlobalIndexMap.cpp b/Tests/NumLib/LocalToGlobalIndexMap.cpp
index 1f4fda2310dfddb3bc902ff5c60a5cffef1428ca..c2cf558a51132ac93a7ede40b9ad46b45d6d8ddc 100644
--- a/Tests/NumLib/LocalToGlobalIndexMap.cpp
+++ b/Tests/NumLib/LocalToGlobalIndexMap.cpp
@@ -174,3 +174,108 @@ TEST_F(NumLibLocalToGlobalIndexMapTest, DISABLED_MultipleVariablesMultipleCompon
ASSERT_EQ(10, dof_map->getGlobalIndex(l_node0, 0, 1));
ASSERT_EQ(20, dof_map->getGlobalIndex(l_node0, 1, 0));
}
+
+
+#ifndef USE_PETSC
+TEST_F(NumLibLocalToGlobalIndexMapTest, MultipleVariablesMultipleComponentsHeterogeneousElements)
+#else
+TEST_F(NumLibLocalToGlobalIndexMapTest, DISABLED_MultipleVariablesMultipleComponentsHeterogeneousElements)
+#endif
+{
+ // test 2 variables
+ // - 1st variable with 2 components for all nodes, elements
+ // - 2nd variable with 1 component for nodes of element id 1
+ std::vector<MeshLib::Node*> var2_nodes{const_cast<MeshLib::Node*>(mesh->getNode(1)), const_cast<MeshLib::Node*>(mesh->getNode(2))};
+ std::unique_ptr<MeshLib::MeshSubset> var2_subset(new MeshLib::MeshSubset(*mesh, &var2_nodes));
+ components.emplace_back(new MeshLib::MeshSubsets{var2_subset.get()});
+
+ std::vector<unsigned> vec_var_n_components{2, 1};
+ std::vector<std::vector<MeshLib::Element*>const*> vec_var_elements;
+ vec_var_elements.push_back(&mesh->getElements());
+ std::vector<MeshLib::Element*> var2_elements{const_cast<MeshLib::Element*>(mesh->getElement(1))};
+ vec_var_elements.push_back(&var2_elements);
+
+ dof_map.reset(new NumLib::LocalToGlobalIndexMap(
+ std::move(components),
+ vec_var_n_components,
+ vec_var_elements,
+ NumLib::ComponentOrder::BY_COMPONENT));
+
+ ASSERT_EQ(22u, dof_map->dofSizeWithGhosts());
+ ASSERT_EQ(3u, dof_map->getNumberOfComponents());
+ ASSERT_EQ(2u, dof_map->getNumberOfVariables());
+ ASSERT_EQ(2u, dof_map->getNumberOfVariableComponents(0));
+ ASSERT_EQ(1u, dof_map->getNumberOfVariableComponents(1));
+
+ MeshLib::Location l_node0(mesh->getID(), MeshLib::MeshItemType::Node, 0);
+
+ ASSERT_EQ(0, dof_map->getGlobalIndex(l_node0, 0, 0));
+ ASSERT_EQ(10, dof_map->getGlobalIndex(l_node0, 0, 1));
+ ASSERT_EQ(std::numeric_limits<GlobalIndexType>::max(), dof_map->getGlobalIndex(l_node0, 1, 0));
+
+ MeshLib::Location l_node1(mesh->getID(), MeshLib::MeshItemType::Node, 1);
+ ASSERT_EQ(1, dof_map->getGlobalIndex(l_node1, 0, 0));
+ ASSERT_EQ(11, dof_map->getGlobalIndex(l_node1, 0, 1));
+ ASSERT_EQ(20, dof_map->getGlobalIndex(l_node1, 1, 0));
+
+ auto ele0_c0_indices = (*dof_map)(0, 0);
+ ASSERT_EQ(2u, ele0_c0_indices.rows.size());
+ auto ele0_c2_indices = (*dof_map)(0, 2);
+ ASSERT_EQ(0u, ele0_c2_indices.rows.size());
+
+ auto ele1_c2_indices = (*dof_map)(1, 2);
+ ASSERT_EQ(2u, ele1_c2_indices.rows.size());
+}
+
+
+#ifndef USE_PETSC
+TEST_F(NumLibLocalToGlobalIndexMapTest, MultipleVariablesMultipleComponentsHeterogeneousWithinElement)
+#else
+TEST_F(NumLibLocalToGlobalIndexMapTest, DISABLED_MultipleVariablesMultipleComponentsHeterogeneousWithinElement)
+#endif
+{
+ // test 2 variables
+ // - 1st variable with 2 components for all nodes, elements
+ // - 2nd variable with 1 component for 1st node of element id 1
+ std::vector<MeshLib::Node*> var2_nodes{const_cast<MeshLib::Node*>(mesh->getNode(1))};
+ std::unique_ptr<MeshLib::MeshSubset> var2_subset(new MeshLib::MeshSubset(*mesh, &var2_nodes));
+ components.emplace_back(new MeshLib::MeshSubsets{var2_subset.get()});
+
+ std::vector<unsigned> vec_var_n_components{2, 1};
+ std::vector<std::vector<MeshLib::Element*>const*> vec_var_elements;
+ vec_var_elements.push_back(&mesh->getElements());
+ std::vector<MeshLib::Element*> var2_elements{const_cast<MeshLib::Element*>(mesh->getElement(1))};
+ vec_var_elements.push_back(&var2_elements);
+
+ dof_map.reset(new NumLib::LocalToGlobalIndexMap(
+ std::move(components),
+ vec_var_n_components,
+ vec_var_elements,
+ NumLib::ComponentOrder::BY_COMPONENT));
+
+ ASSERT_EQ(21u, dof_map->dofSizeWithGhosts());
+ ASSERT_EQ(3u, dof_map->getNumberOfComponents());
+ ASSERT_EQ(2u, dof_map->getNumberOfVariables());
+ ASSERT_EQ(2u, dof_map->getNumberOfVariableComponents(0));
+ ASSERT_EQ(1u, dof_map->getNumberOfVariableComponents(1));
+
+ MeshLib::Location l_node0(mesh->getID(), MeshLib::MeshItemType::Node, 0);
+
+ ASSERT_EQ(0, dof_map->getGlobalIndex(l_node0, 0, 0));
+ ASSERT_EQ(10, dof_map->getGlobalIndex(l_node0, 0, 1));
+ ASSERT_EQ(std::numeric_limits<GlobalIndexType>::max(), dof_map->getGlobalIndex(l_node0, 1, 0));
+
+ MeshLib::Location l_node1(mesh->getID(), MeshLib::MeshItemType::Node, 1);
+ ASSERT_EQ(1, dof_map->getGlobalIndex(l_node1, 0, 0));
+ ASSERT_EQ(11, dof_map->getGlobalIndex(l_node1, 0, 1));
+ ASSERT_EQ(20, dof_map->getGlobalIndex(l_node1, 1, 0));
+
+ auto ele0_c0_indices = (*dof_map)(0, 0);
+ ASSERT_EQ(2u, ele0_c0_indices.rows.size());
+ auto ele0_c2_indices = (*dof_map)(0, 2);
+ ASSERT_EQ(0u, ele0_c2_indices.rows.size());
+
+ auto ele1_c2_indices = (*dof_map)(1, 2);
+ ASSERT_EQ(1u, ele1_c2_indices.rows.size());
+ ASSERT_EQ(20u, ele1_c2_indices.rows[0]);
+}