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
131std::vector<GlobalIndexType> getIndices(
132 std::size_t const variable_id,
133 std::size_t const component_id,
134 std::size_t const bulk_element_id,
135 NumLib::LocalToGlobalIndexMap const& dof_table)
136{
137 auto const& ms =
138 dof_table.getMeshSubset(variable_id, component_id).getMesh();
139 std::vector<GlobalIndexType> indices;
140 auto const& bulk_element = ms.getElement(bulk_element_id);
141 auto const global_component =
142 dof_table.getGlobalComponent(variable_id, component_id);
143 auto const mesh_id = ms.getID();
144
145 for (auto const node_id : bulk_element->nodes() | MeshLib::views::ids)
146 {
148 node_id};
149 // TODO: possible issues on higher order meshes with first order
150 // variables (e.g. pressure in HM).
151 indices.push_back(dof_table.getGlobalIndex(loc, global_component));
152 }
153 return indices;
154}
155
157 std::size_t const id,
158 NumLib::LocalToGlobalIndexMap const& dof_table,
159 std::vector<GlobalIndexType>& indices)
160{
161 assert(dof_table.size() > id);
162 indices.clear();
163
164 // Local matrices and vectors will always be ordered by component,
165 // no matter what the order of the global matrix is.
166 for (int c = 0; c < dof_table.getNumberOfGlobalComponents(); ++c)
167 {
168 auto const& idcs = dof_table(id, c).rows;
169 indices.reserve(indices.size() + idcs.size());
170 indices.insert(indices.end(), idcs.begin(), idcs.end());
171 }
172
174}
175
176double norm(GlobalVector const& x, unsigned const global_component,
177 MathLib::VecNormType norm_type,
178 LocalToGlobalIndexMap const& dof_table)
179{
180 switch (norm_type)
181 {
183 return norm1(x, global_component, dof_table);
185 return norm2(x, global_component, dof_table);
187 return normInfinity(x, global_component, dof_table);
188 default:
189 OGS_FATAL("An invalid norm type has been passed.");
190 }
191}
192
193Eigen::VectorXd getLocalX(
194 std::size_t const mesh_item_id,
195 std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_tables,
196 std::vector<GlobalVector*> const& x)
197{
198 Eigen::VectorXd local_x_vec;
199
200 auto const n_processes = x.size();
201 for (std::size_t process_id = 0; process_id < n_processes; ++process_id)
202 {
203 auto const indices =
204 NumLib::getIndices(mesh_item_id, *dof_tables[process_id]);
205 assert(!indices.empty());
206 auto const last = local_x_vec.size();
207 local_x_vec.conservativeResize(last + indices.size());
208 auto const local_solution = x[process_id]->get(indices);
209 assert(indices.size() == local_solution.size());
210 local_x_vec.tail(local_solution.size()).noalias() =
211 MathLib::toVector(local_solution);
212 }
213 return local_x_vec;
214}
215
216} // 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
Mesh const & getMesh() const
Definition MeshSubset.h:92
const Element * getElement(std::size_t idx) const
Get the element with the given index.
Definition Mesh.h:96
int getGlobalComponent(int const variable_id, int const component_id) const
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.
constexpr ranges::views::view_closure ids
For an element of a range view return its id.
Definition Mesh.h:227
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)