OGS 6.1.0-1699-ge946d4c5f
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) // ghost node value
102  return 0.0;
103 
104  return x.get(index);
105 }
106 
107 double getNodalValue(GlobalVector const& x, MeshLib::Mesh const& mesh,
108  NumLib::LocalToGlobalIndexMap const& dof_table,
109  std::size_t const node_id,
110  std::size_t const global_component_id)
111 {
113  node_id};
114 
115  auto const index = dof_table.getGlobalIndex(l, global_component_id);
116  assert (index != NumLib::MeshComponentMap::nop);
117 
118  return x.get(index);
119 }
120 
121 std::vector<GlobalIndexType> getIndices(
122  std::size_t const mesh_item_id,
123  NumLib::LocalToGlobalIndexMap const& dof_table)
124 {
125  assert(dof_table.size() > mesh_item_id);
126  std::vector<GlobalIndexType> indices;
127 
128  // Local matrices and vectors will always be ordered by component
129  // no matter what the order of the global matrix is.
130  for (int c = 0; c < dof_table.getNumberOfComponents(); ++c)
131  {
132  auto const& idcs = dof_table(mesh_item_id, c).rows;
133  indices.reserve(indices.size() + idcs.size());
134  indices.insert(indices.end(), idcs.begin(), idcs.end());
135  }
136 
137  return indices;
138 }
139 
141  std::size_t const id,
142  NumLib::LocalToGlobalIndexMap const& dof_table,
143  std::vector<GlobalIndexType>& indices)
144 {
145  assert(dof_table.size() > id);
146  indices.clear();
147 
148  // Local matrices and vectors will always be ordered by component,
149  // no matter what the order of the global matrix is.
150  for (int c = 0; c < dof_table.getNumberOfComponents(); ++c)
151  {
152  auto const& idcs = dof_table(id, c).rows;
153  indices.reserve(indices.size() + idcs.size());
154  indices.insert(indices.end(), idcs.begin(), idcs.end());
155  }
156 
157  return NumLib::LocalToGlobalIndexMap::RowColumnIndices(indices, indices);
158 }
159 
160 double norm(GlobalVector const& x, unsigned const global_component,
161  MathLib::VecNormType norm_type,
162  LocalToGlobalIndexMap const& dof_table, MeshLib::Mesh const& mesh)
163 {
164  switch (norm_type)
165  {
167  return norm1(x, global_component, dof_table, mesh);
169  return norm2(x, global_component, dof_table, mesh);
171  return normInfinity(x, global_component, dof_table, mesh);
172  default:
173  OGS_FATAL("An invalid norm type has been passed.");
174  }
175 }
176 
177 } // 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: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:71
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