OGS
DOFTableUtil.cpp
Go to the documentation of this file.
1
11#include "DOFTableUtil.h"
12
13#include <algorithm>
14#include <cassert>
15#include <functional>
16
17#include "BaseLib/MPI.h"
19namespace NumLib
20{
21namespace
22{
23template <typename CalculateNorm>
24double norm(GlobalVector const& x, unsigned const global_component,
25 LocalToGlobalIndexMap const& dof_table, MeshLib::Mesh const& mesh,
26 CalculateNorm calculate_norm)
27{
28#ifdef USE_PETSC
29 x.setLocalAccessibleVector();
30#endif
31
32 double res = 0.0;
33 auto const& ms = dof_table.getMeshSubset(global_component);
34
35 assert(ms.getMeshID() == mesh.getID());
36 for (MeshLib::Node const* node : ms.getNodes())
37 {
38 auto const value = getNonGhostNodalValue(
39 x, mesh, dof_table, node->getID(), global_component);
40 res = calculate_norm(res, value);
41 }
42 return res;
43}
44
45double norm1(GlobalVector const& x, unsigned const global_component,
46 LocalToGlobalIndexMap const& dof_table, MeshLib::Mesh const& mesh)
47{
48 double const res =
49 norm(x, global_component, dof_table, mesh,
50 [](double res, double value) { return res + std::abs(value); });
51
52#ifdef USE_PETSC
53 return BaseLib::MPI::allreduce(res, MPI_SUM, BaseLib::MPI::Mpi{});
54#endif
55 return res;
56}
57
58double norm2(GlobalVector const& x, unsigned const global_component,
59 LocalToGlobalIndexMap const& dof_table, MeshLib::Mesh const& mesh)
60{
61 double const res =
62 norm(x, global_component, dof_table, mesh,
63 [](double res, double value) { return res + value * value; });
64
65#ifdef USE_PETSC
66 return std::sqrt(
68#endif
69 return std::sqrt(res);
70}
71
72double normInfinity(GlobalVector const& x, unsigned const global_component,
73 LocalToGlobalIndexMap const& dof_table,
74 MeshLib::Mesh const& mesh)
75{
76 double const res =
77 norm(x, global_component, dof_table, mesh, [](double res, double value)
78 { return std::max(res, std::abs(value)); });
79
80#ifdef USE_PETSC
81 return BaseLib::MPI::allreduce(res, MPI_MAX, BaseLib::MPI::Mpi{});
82#endif
83 return res;
84}
85} // anonymous namespace
86
88 NumLib::LocalToGlobalIndexMap const& dof_table,
89 std::size_t const node_id,
90 std::size_t const global_component_id)
91{
93 node_id};
94
95 auto const index = dof_table.getGlobalIndex(l, global_component_id);
96 assert(index != NumLib::MeshComponentMap::nop);
97
98 if (index < 0)
99 { // ghost node value
100 return 0.0;
101 }
102
103 return x.get(index);
104}
105
106double getNodalValue(GlobalVector const& x, MeshLib::Mesh const& mesh,
107 NumLib::LocalToGlobalIndexMap const& dof_table,
108 std::size_t const node_id,
109 std::size_t const global_component_id)
110{
112 node_id};
113
114 auto const index = dof_table.getGlobalIndex(l, global_component_id);
115 assert(index != NumLib::MeshComponentMap::nop);
116
117 return x.get(index);
118}
119
120std::vector<GlobalIndexType> getIndices(
121 std::size_t const mesh_item_id,
122 NumLib::LocalToGlobalIndexMap const& dof_table)
123{
124 assert(dof_table.size() > mesh_item_id);
125 std::vector<GlobalIndexType> indices;
126
127 // Local matrices and vectors will always be ordered by component
128 // no matter what the order of the global matrix is.
129 for (int c = 0; c < dof_table.getNumberOfGlobalComponents(); ++c)
130 {
131 auto const& idcs = dof_table(mesh_item_id, c).rows;
132 indices.reserve(indices.size() + idcs.size());
133 indices.insert(indices.end(), idcs.begin(), idcs.end());
134 }
135
136 return indices;
137}
138
140 std::size_t const id,
141 NumLib::LocalToGlobalIndexMap const& dof_table,
142 std::vector<GlobalIndexType>& indices)
143{
144 assert(dof_table.size() > id);
145 indices.clear();
146
147 // Local matrices and vectors will always be ordered by component,
148 // no matter what the order of the global matrix is.
149 for (int c = 0; c < dof_table.getNumberOfGlobalComponents(); ++c)
150 {
151 auto const& idcs = dof_table(id, c).rows;
152 indices.reserve(indices.size() + idcs.size());
153 indices.insert(indices.end(), idcs.begin(), idcs.end());
154 }
155
157}
158
159double norm(GlobalVector const& x, unsigned const global_component,
160 MathLib::VecNormType norm_type,
161 LocalToGlobalIndexMap const& dof_table, MeshLib::Mesh const& mesh)
162{
163 switch (norm_type)
164 {
166 return norm1(x, global_component, dof_table, mesh);
168 return norm2(x, global_component, dof_table, mesh);
170 return normInfinity(x, global_component, dof_table, mesh);
171 default:
172 OGS_FATAL("An invalid norm type has been passed.");
173 }
174}
175
176Eigen::VectorXd getLocalX(
177 std::size_t const mesh_item_id,
178 std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_tables,
179 std::vector<GlobalVector*> const& x)
180{
181 Eigen::VectorXd local_x_vec;
182
183 auto const n_processes = x.size();
184 for (std::size_t process_id = 0; process_id < n_processes; ++process_id)
185 {
186 auto const indices =
187 NumLib::getIndices(mesh_item_id, *dof_tables[process_id]);
188 assert(!indices.empty());
189 auto const last = local_x_vec.size();
190 local_x_vec.conservativeResize(last + indices.size());
191 auto const local_solution = x[process_id]->get(indices);
192 assert(indices.size() == local_solution.size());
193 local_x_vec.tail(local_solution.size()).noalias() =
194 MathLib::toVector(local_solution);
195 }
196 return local_x_vec;
197}
198
199} // namespace NumLib
#define OGS_FATAL(...)
Definition Error.h:26
Global vector based on Eigen vector.
Definition EigenVector.h:25
double get(IndexType rowId) const
get entry
Definition EigenVector.h:58
std::size_t getID() const
Get id of the mesh.
Definition Mesh.h:121
MathLib::RowColumnIndices< GlobalIndexType > RowColumnIndices
GlobalIndexType getGlobalIndex(MeshLib::Location const &l, int const variable_id, int const component_id) const
MeshLib::MeshSubset const & getMeshSubset(int const variable_id, int const component_id) const
static constexpr NUMLIB_EXPORT GlobalIndexType const nop
static T allreduce(T const &value, MPI_Op const &mpi_op, Mpi const &mpi)
Definition MPI.h:128
Eigen::Map< const Vector > toVector(std::vector< double > const &data, Eigen::VectorXd::Index size)
Creates an Eigen mapped vector from the given data vector.
double normInfinity(GlobalVector const &x, unsigned const global_component, LocalToGlobalIndexMap const &dof_table, MeshLib::Mesh const &mesh)
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)
double norm(GlobalVector const &x, unsigned const global_component, MathLib::VecNormType norm_type, LocalToGlobalIndexMap const &dof_table, MeshLib::Mesh const &mesh)
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)
std::vector< GlobalIndexType > getIndices(std::size_t const mesh_item_id, NumLib::LocalToGlobalIndexMap const &dof_table)
NumLib::LocalToGlobalIndexMap::RowColumnIndices getRowColumnIndices(std::size_t const id, NumLib::LocalToGlobalIndexMap const &dof_table, std::vector< GlobalIndexType > &indices)
Eigen::VectorXd getLocalX(std::size_t const mesh_item_id, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_tables, std::vector< GlobalVector * > const &x)