OGS 6.2.2-330-gf48c72f61.dirty.20200225212913
DOFTableUtil.cpp
Go to the documentation of this file.
1 
11 #include "DOFTableUtil.h"
12 #include <cassert>
13 
14 namespace NumLib
15 {
16 namespace
17 {
18 template <typename CalculateNorm>
19 double norm(GlobalVector const& x, unsigned const global_component,
20  LocalToGlobalIndexMap const& dof_table, MeshLib::Mesh const& mesh,
21  CalculateNorm calculate_norm)
22 {
23 #ifdef USE_PETSC
24  x.setLocalAccessibleVector();
25 #endif
26 
27  double res = 0.0;
28  auto const& ms = dof_table.getMeshSubset(global_component);
29 
30  assert(ms.getMeshID() == mesh.getID());
31  for (MeshLib::Node const* node : ms.getNodes())
32  {
33  auto const value = getNonGhostNodalValue(
34  x, mesh, dof_table, node->getID(), global_component);
35  res = calculate_norm(res, value);
36  }
37  return res;
38 }
39 
40 double norm1(GlobalVector const& x, unsigned const global_component,
41  LocalToGlobalIndexMap const& dof_table, MeshLib::Mesh const& mesh)
42 {
43  double res =
44  norm(x, global_component, dof_table, mesh,
45  [](double res, double value) { return res + std::abs(value); });
46 
47 #ifdef USE_PETSC
48  double global_result = 0.0;
49  MPI_Allreduce(&res, &global_result, 1, MPI_DOUBLE, MPI_SUM,
50  PETSC_COMM_WORLD);
51  res = global_result;
52 #endif
53  return res;
54 }
55 
56 double norm2(GlobalVector const& x, unsigned const global_component,
57  LocalToGlobalIndexMap const& dof_table, MeshLib::Mesh const& mesh)
58 {
59  double res =
60  norm(x, global_component, dof_table, mesh,
61  [](double res, double value) { return res + value * value; });
62 
63 #ifdef USE_PETSC
64  double global_result = 0.0;
65  MPI_Allreduce(&res, &global_result, 1, MPI_DOUBLE, MPI_SUM,
66  PETSC_COMM_WORLD);
67  res = global_result;
68 #endif
69  return std::sqrt(res);
70 }
71 
72 double normInfinity(GlobalVector const& x, unsigned const global_component,
73  LocalToGlobalIndexMap const& dof_table,
74  MeshLib::Mesh const& mesh)
75 {
76  double res = norm(x, global_component, dof_table, mesh,
77  [](double res, double value) {
78  return std::max(res, std::abs(value));
79  });
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.getNumberOfComponents(); ++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.getNumberOfComponents(); ++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
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 norm1(GlobalVector const &x, unsigned const global_component, LocalToGlobalIndexMap const &dof_table, MeshLib::Mesh const &mesh)
MathLib::RowColumnIndices< GlobalIndexType > RowColumnIndices
std::size_t getID() const
Get id of the mesh.
Definition: Mesh.h:120
VecNormType
Norm type. Not declared as class type in order to use the members as integers.
Definition: LinAlgEnums.h:21
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)
double norm2(GlobalVector const &x, unsigned const global_component, LocalToGlobalIndexMap const &dof_table, MeshLib::Mesh const &mesh)
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
#define OGS_FATAL(fmt,...)
Definition: Error.h:64
NumLib::LocalToGlobalIndexMap::RowColumnIndices getRowColumnIndices(std::size_t const id, NumLib::LocalToGlobalIndexMap const &dof_table, std::vector< GlobalIndexType > &indices)
double norm(GlobalVector const &x, unsigned const global_component, MathLib::VecNormType norm_type, LocalToGlobalIndexMap const &dof_table, MeshLib::Mesh const &mesh)
static NUMLIB_EXPORT GlobalIndexType const nop