diff --git a/MeshLib/MeshSubset.h b/MeshLib/MeshSubset.h
index 6323599928fbb07ff50b609b6a1e2072cd1dd301..49f7291d718cf0c6fe56692178eb2e118cc475ec 100644
--- a/MeshLib/MeshSubset.h
+++ b/MeshLib/MeshSubset.h
@@ -125,6 +125,12 @@ public:
return _msh.getElements().cend();
}
+ std::vector<Node*> const& getNodes() const
+ {
+ assert(_nodes);
+ return *_nodes;
+ }
+
/// Constructs a new mesh subset which is a set intersection of the current
/// nodes and the provided vector of nodes.
/// An empty mesh subset may be returned, not a nullptr, in case of empty
diff --git a/NumLib/DOF/LocalToGlobalIndexMap.cpp b/NumLib/DOF/LocalToGlobalIndexMap.cpp
index aed289d168388ef050ed63b2acef9b032778d6e0..6dcfd18a79f6b64b0c31a32d948d45afdf220132 100644
--- a/NumLib/DOF/LocalToGlobalIndexMap.cpp
+++ b/NumLib/DOF/LocalToGlobalIndexMap.cpp
@@ -16,6 +16,7 @@ template <typename ElementIterator>
void
LocalToGlobalIndexMap::findGlobalIndices(
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)
{
@@ -26,16 +27,16 @@ LocalToGlobalIndexMap::findGlobalIndices(
std::size_t elem_id = 0;
for (ElementIterator e = first; e != last; ++e, ++elem_id)
{
- std::size_t const nnodes = (*e)->getNumberOfNodes();
-
LineIndex indices;
- indices.reserve(nnodes);
- for (unsigned n = 0; n < nnodes; n++)
+ 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::find(std::begin(nodes), std::end(nodes), *n) == std::end(nodes))
+ continue;
MeshLib::Location l(mesh_id,
MeshLib::MeshItemType::Node,
- (*e)->getNode(n)->getID());
+ (*n)->getID());
indices.push_back(_mesh_component_map.getGlobalIndex(l, comp_id));
}
@@ -70,8 +71,8 @@ LocalToGlobalIndexMap::LocalToGlobalIndexMap(
auto const global_component_id =
getGlobalComponent(variable_id, component_id);
- findGlobalIndices(ms.elementsBegin(), ms.elementsEnd(), mesh_id,
- global_component_id, global_component_id);
+ findGlobalIndices(ms.elementsBegin(), ms.elementsEnd(), ms.getNodes(),
+ mesh_id, global_component_id, global_component_id);
}
// increase by number of components of that variable
@@ -98,7 +99,7 @@ LocalToGlobalIndexMap::LocalToGlobalIndexMap(
{
std::size_t const mesh_id = ms->getMeshID();
- findGlobalIndices(elements.cbegin(), elements.cend(), mesh_id,
+ findGlobalIndices(elements.cbegin(), elements.cend(), ms->getNodes(), mesh_id,
component_id, 0); // There is only one component to
// write out, therefore the zero
// parameter.
diff --git a/NumLib/DOF/LocalToGlobalIndexMap.h b/NumLib/DOF/LocalToGlobalIndexMap.h
index 6fab44ea25fd646b114cca7b4085082bbb193ab1..2b71d267357190beb5d4e6380e5bfac9929c17d1 100644
--- a/NumLib/DOF/LocalToGlobalIndexMap.h
+++ b/NumLib/DOF/LocalToGlobalIndexMap.h
@@ -130,6 +130,7 @@ private:
template <typename ElementIterator>
void
findGlobalIndices(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);
diff --git a/ProcessLib/BoundaryCondition/BoundaryCondition.cpp b/ProcessLib/BoundaryCondition/BoundaryCondition.cpp
index 4365f48f8328824344c95310da63293da6755fa7..270433a5725e896e8ded14746dc265291edfbf13 100644
--- a/ProcessLib/BoundaryCondition/BoundaryCondition.cpp
+++ b/ProcessLib/BoundaryCondition/BoundaryCondition.cpp
@@ -56,6 +56,26 @@ std::unique_ptr<BoundaryCondition> createBoundaryCondition(
// is defined.
std::vector<std::size_t> ids =
mesh_node_searcher.getMeshNodeIDs(config.geometry);
+ // Filter out ids, which are not part of mesh subsets corresponding to
+ // the variable_id and component_id.
+ auto const& mesh_subsets =
+ dof_table.getMeshSubsets(variable_id, config.component_id);
+ auto ids_new_end_iterator = std::end(ids);
+ for (auto const& mesh_subset : mesh_subsets)
+ {
+ auto const& nodes = mesh_subset->getNodes();
+ ids_new_end_iterator = std::remove_if(std::begin(ids), ids_new_end_iterator,
+ [&nodes](std::size_t const node_id)
+ {
+ return std::end(nodes) == std::find_if(
+ std::begin(nodes), std::end(nodes),
+ [&node_id](MeshLib::Node* const node)
+ {
+ return node->getID() == node_id;
+ });
+ });
+ }
+ ids.erase(ids_new_end_iterator, std::end(ids));
return createDirichletBoundaryCondition(
config.config, std::move(ids), dof_table, mesh.getID(), variable_id,
diff --git a/ProcessLib/Process.cpp b/ProcessLib/Process.cpp
index 407b6e12fb2f4415df17de06cdf9dd55f0eec4b0..9feda0c3284979bb8c5123fa92db36d330fbe516 100644
--- a/ProcessLib/Process.cpp
+++ b/ProcessLib/Process.cpp
@@ -60,7 +60,6 @@ void Process::initialize()
void Process::setInitialConditions(double const t, GlobalVector& x)
{
DBUG("Set initial conditions.");
- std::size_t const n_nodes = _mesh.getNumberOfNodes();
SpatialPosition pos;
@@ -73,31 +72,38 @@ void Process::setInitialConditions(double const t, GlobalVector& x)
auto const num_comp = pv.getNumberOfComponents();
- for (std::size_t node_id = 0; node_id < n_nodes; ++node_id)
+ for (int component_id = 0; component_id < num_comp; ++component_id)
{
- MeshLib::Location const l(_mesh.getID(),
- MeshLib::MeshItemType::Node, node_id);
-
- pos.setNodeID(node_id);
- auto const& ic_value = ic(t, pos);
-
- for (int comp_id = 0; comp_id < num_comp; ++comp_id)
+ auto const& mesh_subsets =
+ _local_to_global_index_map->getMeshSubsets(variable_id,
+ component_id);
+ for (auto const& mesh_subset : mesh_subsets)
{
- auto global_index =
- std::abs(_local_to_global_index_map->getGlobalIndex(
- l, variable_id, comp_id));
+ auto const mesh_id = mesh_subset->getMeshID();
+ for (auto const* node : mesh_subset->getNodes())
+ {
+ MeshLib::Location const l(
+ mesh_id, MeshLib::MeshItemType::Node, node->getID());
+
+ pos.setNodeID(node->getID());
+ auto const& ic_value = ic(t, pos);
+
+ auto global_index =
+ std::abs(_local_to_global_index_map->getGlobalIndex(
+ l, variable_id, component_id));
#ifdef USE_PETSC
- // The global indices of the ghost entries of the global matrix
- // or the global vectors need to be set as negative values for
- // equation assembly, however the global indices start from
- // zero. Therefore, any ghost entry with zero index is assigned
- // an negative value of the vector size or the matrix dimension.
- // To assign the initial value for the ghost entries, the
- // negative indices of the ghost entries are restored to zero.
- if (global_index == x.size())
- global_index = 0;
+ // The global indices of the ghost entries of the global matrix
+ // or the global vectors need to be set as negative values for
+ // equation assembly, however the global indices start from
+ // zero. Therefore, any ghost entry with zero index is assigned
+ // an negative value of the vector size or the matrix dimension.
+ // To assign the initial value for the ghost entries, the
+ // negative indices of the ghost entries are restored to zero.
+ if (global_index == x.size())
+ global_index = 0;
#endif
- x.set(global_index, ic_value[comp_id]);
+ x.set(global_index, ic_value[component_id]);
+ }
}
}
}
diff --git a/ProcessLib/ProcessOutput.cpp b/ProcessLib/ProcessOutput.cpp
index f57906d55806ae39537fccc2dcbba194eb6296b9..0bac1d2701bea3abd2c252b3b593e95a630a1e96 100644
--- a/ProcessLib/ProcessOutput.cpp
+++ b/ProcessLib/ProcessOutput.cpp
@@ -67,13 +67,15 @@ void doProcessOutput(
auto const& output_variables = process_output.output_variables;
std::set<std::string> already_output;
- std::size_t const n_nodes = mesh.getNumberOfNodes();
int global_component_offset = 0;
int global_component_offset_next = 0;
// primary variables
- for (ProcessVariable& pv : process_variables)
+ for (int variable_id = 0;
+ variable_id < static_cast<int>(process_variables.size());
+ ++variable_id)
{
+ ProcessVariable& pv = process_variables[variable_id];
int const n_components = pv.getNumberOfComponents();
global_component_offset = global_component_offset_next;
global_component_offset_next += n_components;
@@ -87,22 +89,30 @@ void doProcessOutput(
auto& output_data = pv.getOrCreateMeshProperty();
- for (std::size_t node_id = 0; node_id < n_nodes; ++node_id)
+ auto const num_comp = pv.getNumberOfComponents();
+
+ for (int component_id = 0; component_id < num_comp; ++component_id)
{
- MeshLib::Location const l(mesh.getID(),
- MeshLib::MeshItemType::Node, node_id);
- // TODO extend component ids to multiple process variables.
- for (int component_id = 0; component_id < n_components;
- ++component_id)
+ auto const& mesh_subsets =
+ dof_table.getMeshSubsets(variable_id,
+ component_id);
+ for (auto const& mesh_subset : mesh_subsets)
{
+ auto const mesh_id = mesh_subset->getMeshID();
+ for (auto const* node : mesh_subset->getNodes())
+ {
+ MeshLib::Location const l(
+ mesh_id, MeshLib::MeshItemType::Node, node->getID());
+
auto const global_component_id = global_component_offset + component_id;
auto const index =
dof_table.getLocalIndex(
l, global_component_id, x.getRangeBegin(),
x.getRangeEnd());
- output_data[node_id * n_components + component_id] =
+ output_data[node->getID() * n_components + component_id] =
x_copy[index];
+ }
}
}
}