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