ufz · norihiro-w · Oct 7, 2016 · Sep 30, 2016
diff --git 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
@@ -64,6 +64,21 @@ class LocalToGlobalIndexMap final
         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 @@ class LocalToGlobalIndexMap final
         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
@@ -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]);
+}