error_vector.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2014 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"
26 
27 #include "libmesh/dof_map.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 
36 namespace libMesh
37 {
38 
39 
40 
41 // ------------------------------------------------------------
42 // ErrorVector class member functions
44 {
45  START_LOG ("minimum()", "ErrorVector");
46 
47  const dof_id_type n = this->size();
49 
50  for (dof_id_type i=0; i<n; i++)
51  {
52  // Only positive (or zero) values in the error vector
53  libmesh_assert_greater_equal ((*this)[i], 0.);
54  if (this->is_active_elem(i))
55  min = std::min (min, (*this)[i]);
56  }
57  STOP_LOG ("minimum()", "ErrorVector");
58 
59  // ErrorVectors are for positive values
60  libmesh_assert_greater_equal (min, 0.);
61 
62  return min;
63 }
64 
65 
66 
68 {
69  START_LOG ("mean()", "ErrorVector");
70 
71  const dof_id_type n = 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  STOP_LOG ("mean()", "ErrorVector");
85 
86  return the_mean;
87 }
88 
89 
90 
91 
93 {
94  const dof_id_type n = 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 = this->size();
129 
130  START_LOG ("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  STOP_LOG ("variance()", "ErrorVector");
145 
146  return the_variance;
147 }
148 
149 
150 
151 
152 std::vector<dof_id_type> ErrorVector::cut_below(Real cut) const
153 {
154  START_LOG ("cut_below()", "ErrorVector");
155 
156  const dof_id_type n = this->size();
157 
158  std::vector<dof_id_type> cut_indices;
159  cut_indices.reserve(n/2); // Arbitrary
160 
161  for (dof_id_type i=0; i<n; i++)
162  if (this->is_active_elem(i))
163  {
164  if ((*this)[i] < cut)
165  {
166  cut_indices.push_back(i);
167  }
168  }
169 
170  STOP_LOG ("cut_below()", "ErrorVector");
171 
172  return cut_indices;
173 }
174 
175 
176 
177 
178 std::vector<dof_id_type> ErrorVector::cut_above(Real cut) const
179 {
180  START_LOG ("cut_above()", "ErrorVector");
181 
182  const dof_id_type n = this->size();
183 
184  std::vector<dof_id_type> cut_indices;
185  cut_indices.reserve(n/2); // Arbitrary
186 
187  for (dof_id_type i=0; i<n; i++)
188  if (this->is_active_elem(i))
189  {
190  if ((*this)[i] > cut)
191  {
192  cut_indices.push_back(i);
193  }
194  }
195 
196  STOP_LOG ("cut_above()", "ErrorVector");
197 
198  return cut_indices;
199 }
200 
201 
202 
204 {
205  libmesh_assert_less (i, this->size());
206 
207  if (_mesh)
208  {
210  return _mesh->elem(i)->active();
211  }
212  else
213  return ((*this)[i] != 0.);
214 }
215 
216 
217 void ErrorVector::plot_error(const std::string& filename,
218  const MeshBase& oldmesh) const
219 {
220  AutoPtr<MeshBase> meshptr = oldmesh.clone();
221  MeshBase &mesh = *meshptr;
222  mesh.all_first_order();
223  EquationSystems temp_es (mesh);
224  ExplicitSystem& error_system
225  = temp_es.add_system<ExplicitSystem> ("Error");
226  error_system.add_variable("error", CONSTANT, MONOMIAL);
227  temp_es.init();
228 
229  const DofMap& error_dof_map = error_system.get_dof_map();
230 
233  const MeshBase::const_element_iterator end_el =
235  std::vector<dof_id_type> dof_indices;
236 
237  for ( ; el != end_el; ++el)
238  {
239  const Elem* elem = *el;
240 
241  error_dof_map.dof_indices(elem, dof_indices);
242 
243  const dof_id_type elem_id = elem->id();
244 
245  //0 for the monomial basis
246  const dof_id_type solution_index = dof_indices[0];
247 
248  // libMesh::out << "elem_number=" << elem_number << std::endl;
249  libmesh_assert_less (elem_id, (*this).size());
250 
251  // We may have zero error values in special circumstances
252  // libmesh_assert_greater ((*this)[elem_id], 0.);
253  error_system.solution->set(solution_index, (*this)[elem_id]);
254  }
255 
256  // We may have to renumber if the original numbering was not
257  // contiguous. Since this is just a temporary mesh, that's probably
258  // fine.
259  if (mesh.max_elem_id() != mesh.n_elem() ||
260  mesh.max_node_id() != mesh.n_nodes())
261  {
262  mesh.allow_renumbering(true);
264  }
265 
266  if (filename.rfind(".gmv") < filename.size())
267  {
268  GMVIO(mesh).write_discontinuous_gmv(filename,
269  temp_es, false);
270  }
271  else if (filename.rfind(".plt") < filename.size())
272  {
274  (filename, temp_es);
275  }
276 #ifdef LIBMESH_HAVE_EXODUS_API
277  else if( (filename.rfind(".exo") < filename.size()) ||
278  (filename.rfind(".e") < filename.size()) )
279  {
280  ExodusII_IO io(mesh);
281  io.write(filename);
282  io.write_element_data(temp_es);
283  }
284 #endif
285  else
286  {
287  libmesh_here();
288  libMesh::err << "Warning: ErrorVector::plot_error currently only"
289  << " supports .gmv and .plt and .exo/.e (if enabled) output;" << std::endl;
290  libMesh::err << "Could not recognize filename: " << filename
291  << std::endl;
292  }
293 }
294 
295 } // namespace libMesh

Site Created By: libMesh Developers
Last modified: February 07 2014 16:57:04 UTC

Hosted By:
SourceForge.net Logo