libMesh
error_vector.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2017 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/elem.h"
24 #include "libmesh/error_vector.h"
25 #include "libmesh/libmesh_logging.h"
26 
27 #include "libmesh/dof_map.h"
28 #include "libmesh/equation_systems.h"
29 #include "libmesh/explicit_system.h"
30 #include "libmesh/mesh_base.h"
31 #include "libmesh/numeric_vector.h"
32 #include "libmesh/gmv_io.h"
33 #include "libmesh/tecplot_io.h"
34 #include "libmesh/exodusII_io.h"
35 #include "libmesh/xdr_io.h"
36 
37 namespace libMesh
38 {
39 
40 
41 
42 // ------------------------------------------------------------
43 // ErrorVector class member functions
45 {
46  LOG_SCOPE ("minimum()", "ErrorVector");
47 
48  const dof_id_type n = cast_int<dof_id_type>(this->size());
50 
51  for (dof_id_type i=0; i<n; i++)
52  {
53  // Only positive (or zero) values in the error vector
54  libmesh_assert_greater_equal ((*this)[i], 0.);
55  if (this->is_active_elem(i))
56  min = std::min (min, (*this)[i]);
57  }
58 
59  // ErrorVectors are for positive values
60  libmesh_assert_greater_equal (min, 0.);
61 
62  return min;
63 }
64 
65 
66 
68 {
69  LOG_SCOPE ("mean()", "ErrorVector");
70 
71  const dof_id_type n = cast_int<dof_id_type>(this->size());
72 
73  Real the_mean = 0;
74  dof_id_type nnz = 0;
75 
76  for (dof_id_type i=0; i<n; i++)
77  if (this->is_active_elem(i))
78  {
79  the_mean += ( static_cast<Real>((*this)[i]) - the_mean ) / (nnz + 1);
80 
81  nnz++;
82  }
83 
84  return the_mean;
85 }
86 
87 
88 
89 
91 {
92  const dof_id_type n = cast_int<dof_id_type>(this->size());
93 
94  if (n == 0)
95  return 0.;
96 
97 
98  // Build a StatisticsVector<ErrorVectorReal> containing
99  // only our active entries and take its mean
101 
102  sv.reserve (n);
103 
104  for (dof_id_type i=0; i<n; i++)
105  if (this->is_active_elem(i))
106  sv.push_back((*this)[i]);
107 
108  return sv.median();
109 }
110 
111 
112 
113 
115 {
116  ErrorVector ev = (*this);
117 
118  return ev.median();
119 }
120 
121 
122 
123 
124 Real ErrorVector::variance(const Real mean_in) const
125 {
126  const dof_id_type n = cast_int<dof_id_type>(this->size());
127 
128  LOG_SCOPE ("variance()", "ErrorVector");
129 
130  Real the_variance = 0;
131  dof_id_type nnz = 0;
132 
133  for (dof_id_type i=0; i<n; i++)
134  if (this->is_active_elem(i))
135  {
136  const Real delta = ( static_cast<Real>((*this)[i]) - mean_in );
137  the_variance += (delta * delta - the_variance) / (nnz + 1);
138 
139  nnz++;
140  }
141 
142  return the_variance;
143 }
144 
145 
146 
147 
148 std::vector<dof_id_type> ErrorVector::cut_below(Real cut) const
149 {
150  LOG_SCOPE ("cut_below()", "ErrorVector");
151 
152  const dof_id_type n = cast_int<dof_id_type>(this->size());
153 
154  std::vector<dof_id_type> cut_indices;
155  cut_indices.reserve(n/2); // Arbitrary
156 
157  for (dof_id_type i=0; i<n; i++)
158  if (this->is_active_elem(i))
159  {
160  if ((*this)[i] < cut)
161  {
162  cut_indices.push_back(i);
163  }
164  }
165 
166  return cut_indices;
167 }
168 
169 
170 
171 
172 std::vector<dof_id_type> ErrorVector::cut_above(Real cut) const
173 {
174  LOG_SCOPE ("cut_above()", "ErrorVector");
175 
176  const dof_id_type n = cast_int<dof_id_type>(this->size());
177 
178  std::vector<dof_id_type> cut_indices;
179  cut_indices.reserve(n/2); // Arbitrary
180 
181  for (dof_id_type i=0; i<n; i++)
182  if (this->is_active_elem(i))
183  {
184  if ((*this)[i] > cut)
185  {
186  cut_indices.push_back(i);
187  }
188  }
189 
190  return cut_indices;
191 }
192 
193 
194 
196 {
197  libmesh_assert_less (i, this->size());
198 
199  if (_mesh)
200  {
201  return _mesh->elem_ptr(i)->active();
202  }
203  else
204  return ((*this)[i] != 0.);
205 }
206 
207 
208 void ErrorVector::plot_error(const std::string & filename,
209  const MeshBase & oldmesh) const
210 {
211  UniquePtr<MeshBase> meshptr = oldmesh.clone();
212  MeshBase & mesh = *meshptr;
213 
214  // The all_first_order routine will prepare_for_use(), which would
215  // break our ordering if elements get changed.
216  mesh.allow_renumbering(false);
217  mesh.all_first_order();
218 
219 #ifdef LIBMESH_ENABLE_AMR
220  // We don't want p elevation when plotting a single constant value
221  // per element
222  for (auto & elem : mesh.element_ptr_range())
223  {
224  elem->set_p_refinement_flag(Elem::DO_NOTHING);
225  elem->set_p_level(0);
226  }
227 #endif // LIBMESH_ENABLE_AMR
228 
229  EquationSystems temp_es (mesh);
230  ExplicitSystem & error_system
231  = temp_es.add_system<ExplicitSystem> ("Error");
232  error_system.add_variable("error", CONSTANT, MONOMIAL);
233  temp_es.init();
234 
235  const DofMap & error_dof_map = error_system.get_dof_map();
236  std::vector<dof_id_type> dof_indices;
237 
238  for (const auto & elem : mesh.active_local_element_ptr_range())
239  {
240  error_dof_map.dof_indices(elem, dof_indices);
241 
242  const dof_id_type elem_id = elem->id();
243 
244  //0 for the monomial basis
245  const dof_id_type solution_index = dof_indices[0];
246 
247  // libMesh::out << "elem_number=" << elem_number << std::endl;
248  libmesh_assert_less (elem_id, (*this).size());
249 
250  // We may have zero error values in special circumstances
251  // libmesh_assert_greater ((*this)[elem_id], 0.);
252  error_system.solution->set(solution_index, (*this)[elem_id]);
253  }
254 
255  error_system.solution->close();
256 
257  // We may have to renumber if the original numbering was not
258  // contiguous. Since this is just a temporary mesh, that's probably
259  // fine.
260  if (mesh.max_elem_id() != mesh.n_elem() ||
261  mesh.max_node_id() != mesh.n_nodes())
262  {
263  mesh.allow_renumbering(true);
265  }
266 
267  if (filename.rfind(".gmv") < filename.size())
268  {
269  GMVIO(mesh).write_discontinuous_gmv(filename,
270  temp_es, false);
271  }
272  else if (filename.rfind(".plt") < filename.size())
273  {
275  (filename, temp_es);
276  }
277 #ifdef LIBMESH_HAVE_EXODUS_API
278  else if ((filename.rfind(".exo") < filename.size()) ||
279  (filename.rfind(".e") < filename.size()))
280  {
281  ExodusII_IO io(mesh);
282  io.write(filename);
283  io.write_element_data(temp_es);
284  }
285 #endif
286  else if (filename.rfind(".xda") < filename.size())
287  {
288  XdrIO(mesh).write("mesh-"+filename);
289  temp_es.write("soln-"+filename,WRITE,
292  }
293  else if (filename.rfind(".xdr") < filename.size())
294  {
295  XdrIO(mesh,true).write("mesh-"+filename);
296  temp_es.write("soln-"+filename,ENCODE,
299  }
300  else
301  {
302  libmesh_here();
303  libMesh::err << "Warning: ErrorVector::plot_error currently only"
304  << " supports .gmv, .plt, .xdr/.xda, and .exo/.e (if enabled) output;" << std::endl;
305  libMesh::err << "Could not recognize filename: " << filename
306  << std::endl;
307  }
308 }
309 
310 } // namespace libMesh
OStreamProxy err
virtual Real variance() const libmesh_override
Definition: error_vector.h:115
This is the EquationSystems class.
bool active() const
Definition: elem.h:2257
void write(const std::string &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.
void allow_renumbering(bool allow)
If false is passed in then this mesh will no longer be renumbered when being prepared for use...
Definition: mesh_base.h:749
virtual Real median()
Definition: statistics.C:95
MeshBase * _mesh
Pointer to the mesh, which may be used to decide which elements are active.
Definition: error_vector.h:163
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
unsigned int add_variable(const std::string &var, const FEType &type, const std::set< subdomain_id_type > *const active_subdomains=libmesh_nullptr)
Adds the variable var to the list of variables for this system.
Definition: system.C:1101
virtual dof_id_type max_node_id() const =0
MeshBase & mesh
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
virtual void all_first_order()=0
Converts a mesh with higher-order elements into a mesh with linear elements.
long double max(long double a, double b)
MeshIO class used for writing XDR (eXternal Data Representation) and XDA mesh files.
Definition: xdr_io.h:51
DIE A HORRIBLE DEATH HERE typedef float ErrorVectorReal
This is the MeshBase class.
Definition: mesh_base.h:68
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
virtual void write(const std::string &) libmesh_override
This method implements writing a mesh to a specified file.
Definition: xdr_io.C:167
virtual SimpleRange< element_iterator > active_local_element_ptr_range()=0
This class handles the numbering of degrees of freedom on a mesh.
Definition: dof_map.h:167
virtual dof_id_type max_elem_id() const =0
const DofMap & get_dof_map() const
Definition: system.h:2030
virtual SimpleRange< element_iterator > element_ptr_range()=0
void write_discontinuous_gmv(const std::string &name, const EquationSystems &es, const bool write_partitioning, const std::set< std::string > *system_names=libmesh_nullptr) const
Writes a GMV file with discontinuous data.
Definition: gmv_io.C:1556
virtual System & add_system(const std::string &system_type, const std::string &name)
Add the system of type system_type named name to the systems array.
UniquePtr< NumericVector< Number > > solution
Data structure to hold solution values.
Definition: system.h:1523
void plot_error(const std::string &filename, const MeshBase &mesh) 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:208
virtual Real median() libmesh_override
Definition: error_vector.C:90
virtual void write_equation_systems(const std::string &, const EquationSystems &, const std::set< std::string > *system_names=libmesh_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
bool is_active_elem(dof_id_type i) const
Utility function to decide whether element i is active.
Definition: error_vector.C:195
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
virtual void write(const std::string &fname) libmesh_override
This method implements writing a mesh to a specified file.
Definition: exodusII_io.C:789
virtual std::vector< dof_id_type > cut_above(Real cut) const libmesh_override
Definition: error_vector.C:172
void write_element_data(const EquationSystems &es)
Write out element solution.
Definition: exodusII_io.C:535
This class implements writing meshes in the Tecplot format.
Definition: tecplot_io.h:43
virtual void init()
Initialize all the systems.
virtual Real mean() const libmesh_override
Definition: error_vector.C:67
virtual dof_id_type n_nodes() const =0
virtual const Elem * elem_ptr(const dof_id_type i) const =0
virtual dof_id_type n_elem() const =0
virtual UniquePtr< MeshBase > clone() const =0
Virtual "copy constructor".
virtual std::vector< dof_id_type > cut_below(Real cut) const libmesh_override
Definition: error_vector.C:148
long double min(long double a, double b)
virtual ErrorVectorReal minimum() const libmesh_override
Definition: error_vector.C:44
virtual void renumber_nodes_and_elements()=0
After partitioning a mesh it is useful to renumber the nodes and elements so that they lie in contigu...
The ExplicitSystem provides only "right hand side" storage, which should be sufficient for solving mo...
uint8_t dof_id_type
Definition: id_types.h:64
void dof_indices(const Elem *const elem, std::vector< dof_id_type > &di) const
Fills the vector di with the global degree of freedom indices for the element.
Definition: dof_map.C:1917