swh:1:snp:6088ab52ef49920e01e3f334cdf4d5d6c8a822b9
Tip revision: ccd7e321dfa91b05eddac5864df55263b944c2dd authored by wenqing on 12 May 2021, 10:04:07 UTC
Merge branch 'fix_MPL' into 'master'
Merge branch 'fix_MPL' into 'master'
Tip revision: ccd7e32
MeshComponentMap.cpp
/**
* \author Norihiro Watanabe
* \date 2013-04-16
*
* \file
* \copyright
* Copyright (c) 2012-2021, OpenGeoSys Community (http://www.opengeosys.org)
* Distributed under a Modified BSD License.
* See accompanying file LICENSE.txt or
* http://www.opengeosys.org/project/license
*/
#include "MeshComponentMap.h"
#include "BaseLib/Error.h"
#include "MeshLib/MeshSubset.h"
#ifdef USE_PETSC
#include "MeshLib/NodePartitionedMesh.h"
#endif
namespace NumLib
{
using namespace detail;
GlobalIndexType const MeshComponentMap::nop =
std::numeric_limits<GlobalIndexType>::max();
#ifdef USE_PETSC
MeshComponentMap::MeshComponentMap(
std::vector<MeshLib::MeshSubset> const& components, ComponentOrder order)
{
// Use PETSc with single thread
const MeshLib::NodePartitionedMesh& partitioned_mesh =
static_cast<const MeshLib::NodePartitionedMesh&>(
components[0].getMesh());
if (partitioned_mesh.isForSingleThread())
{
createSerialMeshComponentMap(components, order);
return;
}
else
{
createParallelMeshComponentMap(components, order);
}
}
#else
MeshComponentMap::MeshComponentMap(
std::vector<MeshLib::MeshSubset> const& components, ComponentOrder order)
{
createSerialMeshComponentMap(components, order);
}
#endif // end of USE_PETSC
MeshComponentMap MeshComponentMap::getSubset(
std::vector<MeshLib::MeshSubset> const& bulk_mesh_subsets,
MeshLib::MeshSubset const& new_mesh_subset,
std::vector<int> const& new_global_component_ids) const
{
{ // Testing first an assumption met later in the code that the meshes for
// all the bulk_mesh_subsets are equal.
auto const first_mismatch =
std::adjacent_find(begin(bulk_mesh_subsets), end(bulk_mesh_subsets),
[](auto const& a, auto const& b) {
return a.getMeshID() != b.getMeshID();
});
if (first_mismatch != end(bulk_mesh_subsets))
{
OGS_FATAL(
"Assumption in the MeshComponentMap violated. Expecting all of "
"mesh ids to be the same, but it is not true for "
"the mesh '{:s}' with id {:d}.",
first_mismatch->getMesh().getName(),
first_mismatch->getMeshID());
}
}
// Mapping of the nodes in the new_mesh_subset to the bulk mesh nodes
auto const& new_mesh_properties = new_mesh_subset.getMesh().getProperties();
if (!new_mesh_properties.template existsPropertyVector<std::size_t>(
"bulk_node_ids"))
{
OGS_FATAL(
"Bulk node ids map expected in the construction of the mesh "
"subset.");
}
auto const& bulk_node_ids_map =
*new_mesh_properties.template getPropertyVector<std::size_t>(
"bulk_node_ids", MeshLib::MeshItemType::Node, 1);
// New dictionary for the subset.
ComponentGlobalIndexDict subset_dict;
std::size_t const new_mesh_id = new_mesh_subset.getMeshID();
// Lookup the locations in the current mesh component map and
// insert the full lines into the new subset dictionary.
for (auto* const node : new_mesh_subset.getNodes())
{
auto const node_id = node->getID();
bool const is_base_node = isBaseNode(*node);
MeshLib::Location const new_location{
new_mesh_id, MeshLib::MeshItemType::Node, node_id};
// Assuming the meshes for all the bulk_mesh_subsets are equal.
MeshLib::Location const bulk_location{
bulk_mesh_subsets.front().getMeshID(), MeshLib::MeshItemType::Node,
bulk_node_ids_map[node_id]};
for (auto component_id : new_global_component_ids)
{
auto const global_index =
getGlobalIndex(bulk_location, component_id);
if (global_index == nop)
{
if (is_base_node)
{
OGS_FATAL(
"Could not find a global index for global component "
"{:d} for the mesh '{:s}', node {:d}, in the "
"corresponding bulk mesh '{:s}' and node {:d}. This "
"happens because the boundary mesh is larger then the "
"definition region of the bulk component, usually "
"because the geometry for the boundary condition is "
"too large.",
component_id,
new_mesh_subset.getMesh().getName(),
node_id,
bulk_mesh_subsets.front().getMesh().getName(),
bulk_node_ids_map[node_id]);
}
continue;
}
subset_dict.insert({new_location, component_id, global_index});
}
}
return MeshComponentMap(subset_dict);
}
void MeshComponentMap::renumberByLocation(GlobalIndexType offset)
{
GlobalIndexType global_index = offset;
auto& m = _dict.get<ByLocation>(); // view as sorted by mesh item
for (auto itr_mesh_item = m.begin(); itr_mesh_item != m.end();
++itr_mesh_item)
{
Line pos = *itr_mesh_item;
pos.global_index = global_index++;
m.replace(itr_mesh_item, pos);
}
}
std::vector<int> MeshComponentMap::getComponentIDs(const Location& l) const
{
auto const& m = _dict.get<ByLocation>();
auto const p = m.equal_range(Line(l));
std::vector<int> vec_compID;
for (auto itr = p.first; itr != p.second; ++itr)
{
vec_compID.push_back(itr->comp_id);
}
return vec_compID;
}
Line MeshComponentMap::getLine(Location const& l, int const comp_id) const
{
auto const& m = _dict.get<ByLocationAndComponent>();
auto const itr = m.find(Line(l, comp_id));
assert(itr != m.end()); // The line must exist in the current dictionary.
return *itr;
}
GlobalIndexType MeshComponentMap::getGlobalIndex(Location const& l,
int const comp_id) const
{
auto const& m = _dict.get<ByLocationAndComponent>();
auto const itr = m.find(Line(l, comp_id));
return itr != m.end() ? itr->global_index : nop;
}
std::vector<GlobalIndexType> MeshComponentMap::getGlobalIndices(
const Location& l) const
{
auto const& m = _dict.get<ByLocation>();
auto const p = m.equal_range(Line(l));
std::vector<GlobalIndexType> global_indices;
for (auto itr = p.first; itr != p.second; ++itr)
{
global_indices.push_back(itr->global_index);
}
return global_indices;
}
std::vector<GlobalIndexType> MeshComponentMap::getGlobalIndicesByLocation(
std::vector<Location> const& ls) const
{
// Create vector of global indices sorted by location containing all
// locations given in ls parameter.
std::vector<GlobalIndexType> global_indices;
global_indices.reserve(ls.size());
auto const& m = _dict.get<ByLocation>();
for (const auto& l : ls)
{
auto const p = m.equal_range(Line(l));
for (auto itr = p.first; itr != p.second; ++itr)
{
global_indices.push_back(itr->global_index);
}
}
return global_indices;
}
std::vector<GlobalIndexType> MeshComponentMap::getGlobalIndicesByComponent(
std::vector<Location> const& ls) const
{
// vector of (Component, global Index) pairs.
using CIPair = std::pair<int, GlobalIndexType>;
std::vector<CIPair> pairs;
pairs.reserve(ls.size());
// Create a sub dictionary containing all lines with location from ls.
auto const& m = _dict.get<ByLocation>();
for (const auto& l : ls)
{
auto const p = m.equal_range(Line(l));
for (auto itr = p.first; itr != p.second; ++itr)
{
pairs.emplace_back(itr->comp_id, itr->global_index);
}
}
auto CIPairLess = [](CIPair const& a, CIPair const& b) {
return a.first < b.first;
};
// Create vector of global indices from sub dictionary sorting by component.
if (!std::is_sorted(pairs.begin(), pairs.end(), CIPairLess))
{
std::stable_sort(pairs.begin(), pairs.end(), CIPairLess);
}
std::vector<GlobalIndexType> global_indices;
global_indices.reserve(pairs.size());
transform(cbegin(pairs), cend(pairs), back_inserter(global_indices),
[&](const auto& pair) { return pair.second; });
return global_indices;
}
GlobalIndexType MeshComponentMap::getLocalIndex(
Location const& l,
int const comp_id,
std::size_t const range_begin,
std::size_t const range_end) const
{
GlobalIndexType const global_index = getGlobalIndex(l, comp_id);
#ifndef USE_PETSC
(void)range_begin;
(void)range_end;
return global_index;
#else
if (global_index >= 0) // non-ghost location.
return global_index - range_begin;
//
// For a ghost location look up the global index in ghost indices.
//
// A special case for a ghost location with global index equal to the size
// of the local vector:
GlobalIndexType const real_global_index =
(-global_index == static_cast<GlobalIndexType>(_num_global_dof))
? 0
: -global_index;
// TODO Find in ghost indices is O(n^2/2) for n being the length of
// _ghosts_indices. Providing an inverted table would be faster.
auto const ghost_index_it = std::find(
_ghosts_indices.begin(), _ghosts_indices.end(), real_global_index);
if (ghost_index_it == _ghosts_indices.end())
{
OGS_FATAL("index {:d} not found in ghost_indices", real_global_index);
}
// Using std::distance on a std::vector is O(1). As long as _ghost_indices
// remains of std::vector type, this shall be fast.
return range_end - range_begin +
std::distance(_ghosts_indices.begin(), ghost_index_it);
#endif
}
void MeshComponentMap::createSerialMeshComponentMap(
std::vector<MeshLib::MeshSubset> const& components, ComponentOrder order)
{
// construct dict (and here we number global_index by component type)
GlobalIndexType global_index = 0;
int comp_id = 0;
for (auto const& c : components)
{
std::size_t const mesh_id = c.getMeshID();
// mesh items are ordered first by node, cell, ....
for (std::size_t j = 0; j < c.getNumberOfNodes(); j++)
{
_dict.insert(Line(
Location(mesh_id, MeshLib::MeshItemType::Node, c.getNodeID(j)),
comp_id, global_index++));
}
comp_id++;
}
_num_local_dof = _dict.size();
if (order == ComponentOrder::BY_LOCATION)
{
renumberByLocation();
}
}
#ifdef USE_PETSC
void MeshComponentMap::createParallelMeshComponentMap(
std::vector<MeshLib::MeshSubset> const& components, ComponentOrder order)
{
if (order != ComponentOrder::BY_LOCATION)
{
// Not allowed in parallel case since this is not suitable to
// arrange non ghost entries of a partition within
// a rank in the parallel computing.
OGS_FATAL(
"Global index in the system of equations can only be numbered by "
"the order type of ComponentOrder::BY_LOCATION");
}
// get number of unknowns
GlobalIndexType num_unknowns = 0;
for (auto const& c : components)
{
const MeshLib::NodePartitionedMesh& partitioned_mesh =
static_cast<const MeshLib::NodePartitionedMesh&>(c.getMesh());
num_unknowns += partitioned_mesh.getNumberOfGlobalNodes();
}
// construct dict (and here we number global_index by component type)
int comp_id = 0;
_num_global_dof = 0;
_num_local_dof = 0;
for (auto const& c : components)
{
assert(dynamic_cast<MeshLib::NodePartitionedMesh const*>(
&c.getMesh()) != nullptr);
std::size_t const mesh_id = c.getMeshID();
const MeshLib::NodePartitionedMesh& partitioned_mesh =
static_cast<const MeshLib::NodePartitionedMesh&>(c.getMesh());
// mesh items are ordered first by node, cell, ....
for (std::size_t j = 0; j < c.getNumberOfNodes(); j++)
{
GlobalIndexType global_index = static_cast<GlobalIndexType>(
components.size() * partitioned_mesh.getGlobalNodeID(j) +
comp_id);
const bool is_ghost = partitioned_mesh.isGhostNode(
partitioned_mesh.getNode(j)->getID());
if (is_ghost)
{
_ghosts_indices.push_back(global_index);
global_index = -global_index;
// If the ghost entry has an index of 0,
// its index is set to the negative value of unknowns.
if (global_index == 0)
{
global_index = -num_unknowns;
}
}
else
{
_num_local_dof++;
}
_dict.insert(Line(Location(mesh_id, MeshLib::MeshItemType::Node, j),
comp_id, global_index));
}
_num_global_dof += partitioned_mesh.getNumberOfGlobalNodes();
comp_id++;
}
}
#endif
} // namespace NumLib