OGS
DOFTableUtil.cpp
Go to the documentation of this file.
1 
11 #include "DOFTableUtil.h"
12 
13 #include <cassert>
14 
15 namespace NumLib
16 {
17 namespace
18 {
19 template <typename CalculateNorm>
20 double norm(GlobalVector const& x, unsigned const global_component,
21  LocalToGlobalIndexMap const& dof_table, MeshLib::Mesh const& mesh,
22  CalculateNorm calculate_norm)
23 {
24 #ifdef USE_PETSC
25  x.setLocalAccessibleVector();
26 #endif
27 
28  double res = 0.0;
29  auto const& ms = dof_table.getMeshSubset(global_component);
30 
31  assert(ms.getMeshID() == mesh.getID());
32  for (MeshLib::Node const* node : ms.getNodes())
33  {
34  auto const value = getNonGhostNodalValue(
35  x, mesh, dof_table, node->getID(), global_component);
36  res = calculate_norm(res, value);
37  }
38  return res;
39 }
40 
41 double norm1(GlobalVector const& x, unsigned const global_component,
42  LocalToGlobalIndexMap const& dof_table, MeshLib::Mesh const& mesh)
43 {
44  double res =
45  norm(x, global_component, dof_table, mesh,
46  [](double res, double value) { return res + std::abs(value); });
47 
48 #ifdef USE_PETSC
49  double global_result = 0.0;
50  MPI_Allreduce(&res, &global_result, 1, MPI_DOUBLE, MPI_SUM,
51  PETSC_COMM_WORLD);
52  res = global_result;
53 #endif
54  return res;
55 }
56 
57 double norm2(GlobalVector const& x, unsigned const global_component,
58  LocalToGlobalIndexMap const& dof_table, MeshLib::Mesh const& mesh)
59 {
60  double res =
61  norm(x, global_component, dof_table, mesh,
62  [](double res, double value) { return res + value * value; });
63 
64 #ifdef USE_PETSC
65  double global_result = 0.0;
66  MPI_Allreduce(&res, &global_result, 1, MPI_DOUBLE, MPI_SUM,
67  PETSC_COMM_WORLD);
68  res = global_result;
69 #endif
70  return std::sqrt(res);
71 }
72 
73 double normInfinity(GlobalVector const& x, unsigned const global_component,
74  LocalToGlobalIndexMap const& dof_table,
75  MeshLib::Mesh const& mesh)
76 {
77  double res = norm(x, global_component, dof_table, mesh,
78  [](double res, double value)
79  { return std::max(res, std::abs(value)); });
80 
81 #ifdef USE_PETSC
82  double global_result = 0.0;
83  MPI_Allreduce(&res, &global_result, 1, MPI_DOUBLE, MPI_MAX,
84  PETSC_COMM_WORLD);
85  res = global_result;
86 #endif
87  return res;
88 }
89 } // anonymous namespace
90 
91 double getNonGhostNodalValue(GlobalVector const& x, MeshLib::Mesh const& mesh,
92  NumLib::LocalToGlobalIndexMap const& dof_table,
93  std::size_t const node_id,
94  std::size_t const global_component_id)
95 {
97  node_id};
98 
99  auto const index = dof_table.getGlobalIndex(l, global_component_id);
100  assert(index != NumLib::MeshComponentMap::nop);
101 
102  if (index < 0)
103  { // ghost node value
104  return 0.0;
105  }
106 
107  return x.get(index);
108 }
109 
110 double getNodalValue(GlobalVector const& x, MeshLib::Mesh const& mesh,
111  NumLib::LocalToGlobalIndexMap const& dof_table,
112  std::size_t const node_id,
113  std::size_t const global_component_id)
114 {
116  node_id};
117 
118  auto const index = dof_table.getGlobalIndex(l, global_component_id);
119  assert(index != NumLib::MeshComponentMap::nop);
120 
121  return x.get(index);
122 }
123 
124 std::vector<GlobalIndexType> getIndices(
125  std::size_t const mesh_item_id,
126  NumLib::LocalToGlobalIndexMap const& dof_table)
127 {
128  assert(dof_table.size() > mesh_item_id);
129  std::vector<GlobalIndexType> indices;
130 
131  // Local matrices and vectors will always be ordered by component
132  // no matter what the order of the global matrix is.
133  for (int c = 0; c < dof_table.getNumberOfGlobalComponents(); ++c)
134  {
135  auto const& idcs = dof_table(mesh_item_id, c).rows;
136  indices.reserve(indices.size() + idcs.size());
137  indices.insert(indices.end(), idcs.begin(), idcs.end());
138  }
139 
140  return indices;
141 }
142 
144  std::size_t const id,
145  NumLib::LocalToGlobalIndexMap const& dof_table,
146  std::vector<GlobalIndexType>& indices)
147 {
148  assert(dof_table.size() > id);
149  indices.clear();
150 
151  // Local matrices and vectors will always be ordered by component,
152  // no matter what the order of the global matrix is.
153  for (int c = 0; c < dof_table.getNumberOfGlobalComponents(); ++c)
154  {
155  auto const& idcs = dof_table(id, c).rows;
156  indices.reserve(indices.size() + idcs.size());
157  indices.insert(indices.end(), idcs.begin(), idcs.end());
158  }
159 
160  return NumLib::LocalToGlobalIndexMap::RowColumnIndices(indices, indices);
161 }
162 
163 double norm(GlobalVector const& x, unsigned const global_component,
164  MathLib::VecNormType norm_type,
165  LocalToGlobalIndexMap const& dof_table, MeshLib::Mesh const& mesh)
166 {
167  switch (norm_type)
168  {
170  return norm1(x, global_component, dof_table, mesh);
172  return norm2(x, global_component, dof_table, mesh);
174  return normInfinity(x, global_component, dof_table, mesh);
175  default:
176  OGS_FATAL("An invalid norm type has been passed.");
177  }
178 }
179 
180 } // namespace NumLib
#define OGS_FATAL(...)
Definition: Error.h:26
Global vector based on Eigen vector.
Definition: EigenVector.h:26
double get(IndexType rowId) const
get entry
Definition: EigenVector.h:62
std::size_t getID() const
Get id of the mesh.
Definition: Mesh.h:110
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 NUMLIB_EXPORT GlobalIndexType const nop
VecNormType
Norm type. Not declared as class type in order to use the members as integers.
Definition: LinAlgEnums.h:22
double norm1(GlobalVector const &x, unsigned const global_component, LocalToGlobalIndexMap const &dof_table, MeshLib::Mesh const &mesh)
double normInfinity(GlobalVector const &x, unsigned const global_component, LocalToGlobalIndexMap const &dof_table, MeshLib::Mesh const &mesh)
double norm2(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)