Loading [MathJax]/extensions/tex2jax.js
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,
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 for (auto const l :
37 {
38 auto const value =
39 getNonGhostNodalValue(x, l, dof_table, 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)
47{
48 double const res =
49 norm(x, global_component, dof_table,
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)
60{
61 double const res =
62 norm(x, global_component, dof_table,
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{
75 double const res =
76 norm(x, global_component, dof_table, [](double res, double value)
77 { return std::max(res, std::abs(value)); });
78
79#ifdef USE_PETSC
80 return BaseLib::MPI::allreduce(res, MPI_MAX, BaseLib::MPI::Mpi{});
81#endif
82 return res;
83}
84} // anonymous namespace
85
87 MeshLib::Location const& location,
88 NumLib::LocalToGlobalIndexMap const& dof_table,
89 std::size_t const global_component_id)
90{
91 auto const index = dof_table.getGlobalIndex(location, global_component_id);
92 assert(index != NumLib::MeshComponentMap::nop);
93
94 if (index < 0)
95 { // ghost node value
96 return 0.0;
97 }
98
99 return x.get(index);
100}
101
102double getNodalValue(GlobalVector const& x, MeshLib::Location const& location,
103 NumLib::LocalToGlobalIndexMap const& dof_table,
104 std::size_t const global_component_id)
105{
106 auto const index = dof_table.getGlobalIndex(location, global_component_id);
107 assert(index != NumLib::MeshComponentMap::nop);
108
109 return x.get(index);
110}
111
112std::vector<GlobalIndexType> getIndices(
113 std::size_t const mesh_item_id,
114 NumLib::LocalToGlobalIndexMap const& dof_table)
115{
116 assert(dof_table.size() > mesh_item_id);
117 std::vector<GlobalIndexType> indices;
118
119 // Local matrices and vectors will always be ordered by component
120 // no matter what the order of the global matrix is.
121 for (int c = 0; c < dof_table.getNumberOfGlobalComponents(); ++c)
122 {
123 auto const& idcs = dof_table(mesh_item_id, c).rows;
124 indices.reserve(indices.size() + idcs.size());
125 indices.insert(indices.end(), idcs.begin(), idcs.end());
126 }
127
128 return indices;
129}
130
132 std::size_t const id,
133 NumLib::LocalToGlobalIndexMap const& dof_table,
134 std::vector<GlobalIndexType>& indices)
135{
136 assert(dof_table.size() > id);
137 indices.clear();
138
139 // Local matrices and vectors will always be ordered by component,
140 // no matter what the order of the global matrix is.
141 for (int c = 0; c < dof_table.getNumberOfGlobalComponents(); ++c)
142 {
143 auto const& idcs = dof_table(id, c).rows;
144 indices.reserve(indices.size() + idcs.size());
145 indices.insert(indices.end(), idcs.begin(), idcs.end());
146 }
147
149}
150
151double norm(GlobalVector const& x, unsigned const global_component,
152 MathLib::VecNormType norm_type,
153 LocalToGlobalIndexMap const& dof_table)
154{
155 switch (norm_type)
156 {
158 return norm1(x, global_component, dof_table);
160 return norm2(x, global_component, dof_table);
162 return normInfinity(x, global_component, dof_table);
163 default:
164 OGS_FATAL("An invalid norm type has been passed.");
165 }
166}
167
168Eigen::VectorXd getLocalX(
169 std::size_t const mesh_item_id,
170 std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_tables,
171 std::vector<GlobalVector*> const& x)
172{
173 Eigen::VectorXd local_x_vec;
174
175 auto const n_processes = x.size();
176 for (std::size_t process_id = 0; process_id < n_processes; ++process_id)
177 {
178 auto const indices =
179 NumLib::getIndices(mesh_item_id, *dof_tables[process_id]);
180 assert(!indices.empty());
181 auto const last = local_x_vec.size();
182 local_x_vec.conservativeResize(last + indices.size());
183 auto const local_solution = x[process_id]->get(indices);
184 assert(indices.size() == local_solution.size());
185 local_x_vec.tail(local_solution.size()).noalias() =
186 MathLib::toVector(local_solution);
187 }
188 return local_x_vec;
189}
190
191} // 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
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:124
Eigen::Map< const Vector > toVector(std::vector< double > const &data, Eigen::VectorXd::Index size)
Creates an Eigen mapped vector from the given data vector.
auto meshLocations(Mesh const &mesh, MeshItemType const item_type)
Definition Mesh.h:238
double normInfinity(GlobalVector const &x, unsigned const global_component, LocalToGlobalIndexMap const &dof_table)
double getNodalValue(GlobalVector const &x, MeshLib::Location const &location, NumLib::LocalToGlobalIndexMap const &dof_table, 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)
double getNonGhostNodalValue(GlobalVector const &x, MeshLib::Location const &location, NumLib::LocalToGlobalIndexMap const &dof_table, std::size_t const global_component_id)
double norm(GlobalVector const &x, unsigned const global_component, MathLib::VecNormType norm_type, LocalToGlobalIndexMap const &dof_table)