libMesh
error_vector.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2025 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 
4 // This library is free software; you can redistribute it and/or
5 // modify it under the terms of the GNU Lesser General Public
6 // License as published by the Free Software Foundation; either
7 // version 2.1 of the License, or (at your option) any later version.
8 
9 // This library is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 // Lesser General Public License for more details.
13 
14 // You should have received a copy of the GNU Lesser General Public
15 // License along with this library; if not, write to the Free Software
16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 
18 
19 // C++ includes
20 #include <limits>
21 
22 // Local includes
23 #include "libmesh/dof_map.h"
24 #include "libmesh/elem.h"
25 #include "libmesh/enum_xdr_mode.h"
26 #include "libmesh/error_vector.h"
27 #include "libmesh/equation_systems.h"
28 #include "libmesh/explicit_system.h"
29 #include "libmesh/libmesh_logging.h"
30 #include "libmesh/mesh_base.h"
31 #include "libmesh/numeric_vector.h"
32 
33 #include "libmesh/exodusII_io.h"
34 #include "libmesh/gmv_io.h"
35 #include "libmesh/nemesis_io.h"
36 #include "libmesh/tecplot_io.h"
37 #include "libmesh/xdr_io.h"
38 
39 namespace libMesh
40 {
41 
42 
43 
44 // ------------------------------------------------------------
45 // ErrorVector class member functions
47 {
48  LOG_SCOPE ("minimum()", "ErrorVector");
49 
50  const dof_id_type n = cast_int<dof_id_type>(this->size());
51  ErrorVectorReal min = std::numeric_limits<ErrorVectorReal>::max();
52 
53  for (dof_id_type i=0; i<n; i++)
54  {
55  // Only positive (or zero) values in the error vector
56  libmesh_assert_greater_equal ((*this)[i], 0.);
57  if (this->is_active_elem(i))
58  min = std::min (min, (*this)[i]);
59  }
60 
61  // ErrorVectors are for positive values
62  libmesh_assert_greater_equal (min, 0.);
63 
64  return min;
65 }
66 
67 
68 
70 {
71  LOG_SCOPE ("mean()", "ErrorVector");
72 
73  const dof_id_type n = cast_int<dof_id_type>(this->size());
74 
75  Real the_mean = 0;
76  dof_id_type nnz = 0;
77 
78  for (dof_id_type i=0; i<n; i++)
79  if (this->is_active_elem(i))
80  {
81  the_mean += ( static_cast<Real>((*this)[i]) - the_mean ) / (nnz + 1);
82 
83  nnz++;
84  }
85 
86  return the_mean;
87 }
88 
89 
90 
91 
93 {
94  const dof_id_type n = cast_int<dof_id_type>(this->size());
95 
96  if (n == 0)
97  return 0.;
98 
99 
100  // Build a StatisticsVector<ErrorVectorReal> containing
101  // only our active entries and take its mean
103 
104  sv.reserve (n);
105 
106  for (dof_id_type i=0; i<n; i++)
107  if (this->is_active_elem(i))
108  sv.push_back((*this)[i]);
109 
110  return sv.median();
111 }
112 
113 
114 
115 
117 {
118  ErrorVector ev = (*this);
119 
120  return ev.median();
121 }
122 
123 
124 
125 
126 Real ErrorVector::variance(const Real mean_in) const
127 {
128  const dof_id_type n = cast_int<dof_id_type>(this->size());
129 
130  LOG_SCOPE ("variance()", "ErrorVector");
131 
132  Real the_variance = 0;
133  dof_id_type nnz = 0;
134 
135  for (dof_id_type i=0; i<n; i++)
136  if (this->is_active_elem(i))
137  {
138  const Real delta = ( static_cast<Real>((*this)[i]) - mean_in );
139  the_variance += (delta * delta - the_variance) / (nnz + 1);
140 
141  nnz++;
142  }
143 
144  return the_variance;
145 }
146 
147 
148 
149 
150 std::vector<dof_id_type> ErrorVector::cut_below(Real cut) const
151 {
152  LOG_SCOPE ("cut_below()", "ErrorVector");
153 
154  const dof_id_type n = cast_int<dof_id_type>(this->size());
155 
156  std::vector<dof_id_type> cut_indices;
157  cut_indices.reserve(n/2); // Arbitrary
158 
159  for (dof_id_type i=0; i<n; i++)
160  if (this->is_active_elem(i))
161  {
162  if ((*this)[i] < cut)
163  {
164  cut_indices.push_back(i);
165  }
166  }
167 
168  return cut_indices;
169 }
170 
171 
172 
173 
174 std::vector<dof_id_type> ErrorVector::cut_above(Real cut) const
175 {
176  LOG_SCOPE ("cut_above()", "ErrorVector");
177 
178  const dof_id_type n = cast_int<dof_id_type>(this->size());
179 
180  std::vector<dof_id_type> cut_indices;
181  cut_indices.reserve(n/2); // Arbitrary
182 
183  for (dof_id_type i=0; i<n; i++)
184  if (this->is_active_elem(i))
185  {
186  if ((*this)[i] > cut)
187  {
188  cut_indices.push_back(i);
189  }
190  }
191 
192  return cut_indices;
193 }
194 
195 
196 
198 {
199  libmesh_assert_less (i, this->size());
200 
201  if (_mesh)
202  {
203  return _mesh->elem_ptr(i)->active();
204  }
205  else
206  return ((*this)[i] != 0.);
207 }
208 
209 
210 void ErrorVector::plot_error(const std::string & filename,
211  const MeshBase & oldmesh,
212  const std::string & data_type) const
213 {
214  std::unique_ptr<MeshBase> meshptr = oldmesh.clone();
215  MeshBase & mesh = *meshptr;
216 
217  // The all_first_order routine will prepare_for_use(), which would
218  // break our ordering if elements get changed.
219  mesh.allow_renumbering(false);
220  mesh.all_first_order();
221 
222 #ifdef LIBMESH_ENABLE_AMR
223  // We don't want p elevation when plotting a single constant value
224  // per element
225  for (auto & elem : mesh.element_ptr_range())
226  {
227  elem->set_p_refinement_flag(Elem::DO_NOTHING);
228  elem->set_p_level(0);
229  }
230 #endif // LIBMESH_ENABLE_AMR
231 
232  EquationSystems temp_es (mesh);
233  ExplicitSystem & error_system
234  = temp_es.add_system<ExplicitSystem> ("Error");
235  error_system.add_variable(data_type, CONSTANT, MONOMIAL);
236 
237  temp_es.init();
238 
239  const DofMap & error_dof_map = error_system.get_dof_map();
240  std::vector<dof_id_type> dof_indices;
241 
242  for (const auto & elem : mesh.active_local_element_ptr_range())
243  {
244  error_dof_map.dof_indices(elem, dof_indices);
245 
246  const dof_id_type elem_id = elem->id();
247 
248  //0 for the monomial basis
249  const dof_id_type solution_index = dof_indices[0];
250 
251  // libMesh::out << "elem_number=" << elem_number << std::endl;
252  libmesh_assert_less (elem_id, (*this).size());
253 
254  // We may have zero error values in special circumstances
255  // libmesh_assert_greater ((*this)[elem_id], 0.);
256  error_system.solution->set(solution_index, (*this)[elem_id]);
257  }
258  error_system.solution->close();
259  error_system.update();
260 
261  // We may have to renumber if the original numbering was not
262  // contiguous. Since this is just a temporary mesh, that's probably
263  // fine.
264  if (mesh.max_elem_id() != mesh.n_elem() ||
265  mesh.max_node_id() != mesh.n_nodes())
266  {
267  mesh.allow_renumbering(true);
268  mesh.renumber_nodes_and_elements();
269  }
270 
271  if (Utility::contains(filename, ".gmv"))
272  {
274  temp_es, false);
275  }
276  else if (Utility::contains(filename, ".plt"))
277  {
279  (filename, temp_es);
280  }
281 #if defined(LIBMESH_HAVE_EXODUS_API) && defined(LIBMESH_HAVE_NEMESIS_API)
282  else if (Utility::contains(filename, ".nem") ||
283  Utility::contains(filename, ".n"))
284  {
285  Nemesis_IO io(mesh);
286  io.write(filename);
287  io.write_element_data(temp_es);
288  }
289 #endif
290 #ifdef LIBMESH_HAVE_EXODUS_API
291  else if (Utility::contains(filename, ".exo") ||
292  Utility::contains(filename, ".e"))
293  {
294  ExodusII_IO io(mesh);
295  io.write(filename);
296  io.write_element_data(temp_es);
297  }
298 #endif
299  else if (Utility::contains(filename, ".xda"))
300  {
301  XdrIO(mesh).write("mesh-"+filename);
302  temp_es.write("soln-"+filename,WRITE,
305  }
306  else if (Utility::contains(filename, ".xdr"))
307  {
308  XdrIO(mesh,true).write("mesh-"+filename);
309  temp_es.write("soln-"+filename,ENCODE,
312  }
313  else
314  {
315  libmesh_here();
316  libMesh::err << "Warning: ErrorVector::plot_error currently only"
317  << " supports .gmv, .plt, .xdr/.xda, and .exo/.e (if enabled) output;" << std::endl;
318  libMesh::err << "Could not recognize filename: " << filename
319  << std::endl;
320  }
321 }
322 
323 } // namespace libMesh
OStreamProxy err
This is the EquationSystems class.
void write(std::string_view name, const XdrMODE, const unsigned int write_flags=(WRITE_DATA), bool partition_agnostic=true) const
Write the systems to disk using the XDR data format.
MPI_Datatype data_type
virtual Real median()
Definition: statistics.C:96
virtual void write_equation_systems(const std::string &, const EquationSystems &, const std::set< std::string > *system_names=nullptr)
This method implements writing a mesh with data to a specified file where the data is taken from the ...
Definition: mesh_output.C:31
MeshBase * _mesh
Pointer to the mesh, which may be used to decide which elements are active.
Definition: error_vector.h:164
void dof_indices(const Elem *const elem, std::vector< dof_id_type > &di) const
Definition: dof_map.C:2229
The ErrorVector is a specialization of the StatisticsVector for error data computed on a finite eleme...
Definition: error_vector.h:50
The ExodusII_IO class implements reading meshes in the ExodusII file format from Sandia National Labs...
Definition: exodusII_io.h:52
MeshBase & mesh
virtual Real variance() const override
Definition: error_vector.h:115
The libMesh namespace provides an interface to certain functionality in the library.
This class implements writing meshes in the GMV format.
Definition: gmv_io.h:46
MeshIO class used for writing XDR (eXternal Data Representation) and XDA mesh files.
Definition: xdr_io.h:51
virtual std::unique_ptr< MeshBase > clone() const =0
Virtual "copy constructor".
DIE A HORRIBLE DEATH HERE typedef float ErrorVectorReal
virtual std::vector< dof_id_type > cut_above(Real cut) const override
Definition: error_vector.C:174
This is the MeshBase class.
Definition: mesh_base.h:75
This class handles the numbering of degrees of freedom on a mesh.
Definition: dof_map.h:179
void plot_error(const std::string &filename, const MeshBase &mesh, const std::string &data_type="error") const
Plots a data file, of a type determined by looking at the file extension in filename, of the error values on the active elements of mesh.
Definition: error_vector.C:210
virtual Real median() override
Definition: error_vector.C:92
The Nemesis_IO class implements reading parallel meshes in the Nemesis file format from Sandia Nation...
Definition: nemesis_io.h:51
virtual void write(const std::string &) override
This method implements writing a mesh to a specified file.
Definition: xdr_io.C:127
void write_element_data(const EquationSystems &es)
Write out element solution in parallel, without localizing the solution vector.
Definition: nemesis_io.C:1466
std::unique_ptr< NumericVector< Number > > solution
Data structure to hold solution values.
Definition: system.h:1593
unsigned int add_variable(std::string_view var, const FEType &type, const std::set< subdomain_id_type > *const active_subdomains=nullptr)
Adds the variable var to the list of variables for this system.
Definition: system.C:1357
bool contains(std::string_view superstring, std::string_view substring)
Look for a substring within a string.
Definition: utility.C:205
virtual std::vector< dof_id_type > cut_below(Real cut) const override
Definition: error_vector.C:150
virtual ErrorVectorReal minimum() const override
Definition: error_vector.C:46
virtual const Elem * elem_ptr(const dof_id_type i) const =0
virtual void update()
Update the local values to reflect the solution on neighboring processors.
Definition: system.C:510
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
virtual void write(const std::string &fname) override
This method implements writing a mesh to a specified file.
Definition: exodusII_io.C:2180
void write_element_data(const EquationSystems &es)
Write out element solution.
Definition: exodusII_io.C:1359
virtual void write(const std::string &base_filename) override
This method implements writing a mesh to a specified file.
Definition: nemesis_io.C:1228
This class implements writing meshes in the Tecplot format.
Definition: tecplot_io.h:43
virtual void init()
Initialize all the systems.
bool is_active_elem(dof_id_type i) const
Utility function to decide whether element i is active.
Definition: error_vector.C:197
virtual System & add_system(std::string_view system_type, std::string_view name)
Add the system of type system_type named name to the systems array.
bool active() const
Definition: elem.h:2941
const DofMap & get_dof_map() const
Definition: system.h:2374
void write_discontinuous_gmv(const std::string &name, const EquationSystems &es, const bool write_partitioning, const std::set< std::string > *system_names=nullptr) const
Writes a GMV file with discontinuous data.
Definition: gmv_io.C:1539
Manages consistently variables, degrees of freedom, and coefficient vectors for explicit systems...
uint8_t dof_id_type
Definition: id_types.h:67
virtual Real mean() const override
Definition: error_vector.C:69