Revision bcfd72c9bcb903cbda4ed849c958d63446abd881 authored by renchao.lu on 01 May 2021, 12:51:27 UTC, committed by renchao.lu on 01 May 2021, 12:51:27 UTC
[MPL] Add overloaded operator functions at medium, phase, and component levels See merge request ogs/ogs!3590
DOFTableUtil.cpp
/**
* \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 "DOFTableUtil.h"
#include <cassert>
namespace NumLib
{
namespace
{
template <typename CalculateNorm>
double norm(GlobalVector const& x, unsigned const global_component,
LocalToGlobalIndexMap const& dof_table, MeshLib::Mesh const& mesh,
CalculateNorm calculate_norm)
{
#ifdef USE_PETSC
x.setLocalAccessibleVector();
#endif
double res = 0.0;
auto const& ms = dof_table.getMeshSubset(global_component);
assert(ms.getMeshID() == mesh.getID());
for (MeshLib::Node const* node : ms.getNodes())
{
auto const value = getNonGhostNodalValue(
x, mesh, dof_table, node->getID(), global_component);
res = calculate_norm(res, value);
}
return res;
}
double norm1(GlobalVector const& x, unsigned const global_component,
LocalToGlobalIndexMap const& dof_table, MeshLib::Mesh const& mesh)
{
double res =
norm(x, global_component, dof_table, mesh,
[](double res, double value) { return res + std::abs(value); });
#ifdef USE_PETSC
double global_result = 0.0;
MPI_Allreduce(&res, &global_result, 1, MPI_DOUBLE, MPI_SUM,
PETSC_COMM_WORLD);
res = global_result;
#endif
return res;
}
double norm2(GlobalVector const& x, unsigned const global_component,
LocalToGlobalIndexMap const& dof_table, MeshLib::Mesh const& mesh)
{
double res =
norm(x, global_component, dof_table, mesh,
[](double res, double value) { return res + value * value; });
#ifdef USE_PETSC
double global_result = 0.0;
MPI_Allreduce(&res, &global_result, 1, MPI_DOUBLE, MPI_SUM,
PETSC_COMM_WORLD);
res = global_result;
#endif
return std::sqrt(res);
}
double normInfinity(GlobalVector const& x, unsigned const global_component,
LocalToGlobalIndexMap const& dof_table,
MeshLib::Mesh const& mesh)
{
double res = norm(x, global_component, dof_table, mesh,
[](double res, double value) {
return std::max(res, std::abs(value));
});
#ifdef USE_PETSC
double global_result = 0.0;
MPI_Allreduce(&res, &global_result, 1, MPI_DOUBLE, MPI_MAX,
PETSC_COMM_WORLD);
res = global_result;
#endif
return res;
}
} // anonymous namespace
double getNonGhostNodalValue(GlobalVector const& x, MeshLib::Mesh const& mesh,
NumLib::LocalToGlobalIndexMap const& dof_table,
std::size_t const node_id,
std::size_t const global_component_id)
{
MeshLib::Location const l{mesh.getID(), MeshLib::MeshItemType::Node,
node_id};
auto const index = dof_table.getGlobalIndex(l, global_component_id);
assert(index != NumLib::MeshComponentMap::nop);
if (index < 0)
{ // ghost node value
return 0.0;
}
return x.get(index);
}
double getNodalValue(GlobalVector const& x, MeshLib::Mesh const& mesh,
NumLib::LocalToGlobalIndexMap const& dof_table,
std::size_t const node_id,
std::size_t const global_component_id)
{
MeshLib::Location const l{mesh.getID(), MeshLib::MeshItemType::Node,
node_id};
auto const index = dof_table.getGlobalIndex(l, global_component_id);
assert(index != NumLib::MeshComponentMap::nop);
return x.get(index);
}
std::vector<GlobalIndexType> getIndices(
std::size_t const mesh_item_id,
NumLib::LocalToGlobalIndexMap const& dof_table)
{
assert(dof_table.size() > mesh_item_id);
std::vector<GlobalIndexType> indices;
// Local matrices and vectors will always be ordered by component
// no matter what the order of the global matrix is.
for (int c = 0; c < dof_table.getNumberOfGlobalComponents(); ++c)
{
auto const& idcs = dof_table(mesh_item_id, c).rows;
indices.reserve(indices.size() + idcs.size());
indices.insert(indices.end(), idcs.begin(), idcs.end());
}
return indices;
}
NumLib::LocalToGlobalIndexMap::RowColumnIndices getRowColumnIndices(
std::size_t const id,
NumLib::LocalToGlobalIndexMap const& dof_table,
std::vector<GlobalIndexType>& indices)
{
assert(dof_table.size() > id);
indices.clear();
// Local matrices and vectors will always be ordered by component,
// no matter what the order of the global matrix is.
for (int c = 0; c < dof_table.getNumberOfGlobalComponents(); ++c)
{
auto const& idcs = dof_table(id, c).rows;
indices.reserve(indices.size() + idcs.size());
indices.insert(indices.end(), idcs.begin(), idcs.end());
}
return NumLib::LocalToGlobalIndexMap::RowColumnIndices(indices, indices);
}
double norm(GlobalVector const& x, unsigned const global_component,
MathLib::VecNormType norm_type,
LocalToGlobalIndexMap const& dof_table, MeshLib::Mesh const& mesh)
{
switch (norm_type)
{
case MathLib::VecNormType::NORM1:
return norm1(x, global_component, dof_table, mesh);
case MathLib::VecNormType::NORM2:
return norm2(x, global_component, dof_table, mesh);
case MathLib::VecNormType::INFINITY_N:
return normInfinity(x, global_component, dof_table, mesh);
default:
OGS_FATAL("An invalid norm type has been passed.");
}
}
} // namespace NumLib
Computing file changes ...