1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2020 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
minimum()46 ErrorVectorReal ErrorVector::minimum() const
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
mean()69 Real ErrorVector::mean() const
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
median()92 Real ErrorVector::median()
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
102 StatisticsVector<ErrorVectorReal> sv;
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
median()116 Real ErrorVector::median() const
117 {
118 ErrorVector ev = (*this);
119
120 return ev.median();
121 }
122
123
124
125
variance(const Real mean_in)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
cut_below(Real cut)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
cut_above(Real cut)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
is_active_elem(dof_id_type i)197 bool ErrorVector::is_active_elem (dof_id_type i) const
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
plot_error(const std::string & filename,const MeshBase & oldmesh)210 void ErrorVector::plot_error(const std::string & filename,
211 const MeshBase & oldmesh) const
212 {
213 std::unique_ptr<MeshBase> meshptr = oldmesh.clone();
214 MeshBase & mesh = *meshptr;
215
216 // The all_first_order routine will prepare_for_use(), which would
217 // break our ordering if elements get changed.
218 mesh.allow_renumbering(false);
219 mesh.all_first_order();
220
221 #ifdef LIBMESH_ENABLE_AMR
222 // We don't want p elevation when plotting a single constant value
223 // per element
224 for (auto & elem : mesh.element_ptr_range())
225 {
226 elem->set_p_refinement_flag(Elem::DO_NOTHING);
227 elem->set_p_level(0);
228 }
229 #endif // LIBMESH_ENABLE_AMR
230
231 EquationSystems temp_es (mesh);
232 ExplicitSystem & error_system
233 = temp_es.add_system<ExplicitSystem> ("Error");
234 error_system.add_variable("error", CONSTANT, MONOMIAL);
235 temp_es.init();
236
237 const DofMap & error_dof_map = error_system.get_dof_map();
238 std::vector<dof_id_type> dof_indices;
239
240 for (const auto & elem : mesh.active_local_element_ptr_range())
241 {
242 error_dof_map.dof_indices(elem, dof_indices);
243
244 const dof_id_type elem_id = elem->id();
245
246 //0 for the monomial basis
247 const dof_id_type solution_index = dof_indices[0];
248
249 // libMesh::out << "elem_number=" << elem_number << std::endl;
250 libmesh_assert_less (elem_id, (*this).size());
251
252 // We may have zero error values in special circumstances
253 // libmesh_assert_greater ((*this)[elem_id], 0.);
254 error_system.solution->set(solution_index, (*this)[elem_id]);
255 }
256
257 error_system.solution->close();
258
259 // We may have to renumber if the original numbering was not
260 // contiguous. Since this is just a temporary mesh, that's probably
261 // fine.
262 if (mesh.max_elem_id() != mesh.n_elem() ||
263 mesh.max_node_id() != mesh.n_nodes())
264 {
265 mesh.allow_renumbering(true);
266 mesh.renumber_nodes_and_elements();
267 }
268
269 if (filename.rfind(".gmv") < filename.size())
270 {
271 GMVIO(mesh).write_discontinuous_gmv(filename,
272 temp_es, false);
273 }
274 else if (filename.rfind(".plt") < filename.size())
275 {
276 TecplotIO (mesh).write_equation_systems
277 (filename, temp_es);
278 }
279 #if defined(LIBMESH_HAVE_EXODUS_API) && defined(LIBMESH_HAVE_NEMESIS_API)
280 else if ((filename.rfind(".nem") < filename.size()) ||
281 (filename.rfind(".n") < filename.size()))
282 {
283 Nemesis_IO io(mesh);
284 io.write(filename);
285 io.write_element_data(temp_es);
286 }
287 #endif
288 #ifdef LIBMESH_HAVE_EXODUS_API
289 else if ((filename.rfind(".exo") < filename.size()) ||
290 (filename.rfind(".e") < filename.size()))
291 {
292 ExodusII_IO io(mesh);
293 io.write(filename);
294 io.write_element_data(temp_es);
295 }
296 #endif
297 else if (filename.rfind(".xda") < filename.size())
298 {
299 XdrIO(mesh).write("mesh-"+filename);
300 temp_es.write("soln-"+filename,WRITE,
301 EquationSystems::WRITE_DATA |
302 EquationSystems::WRITE_ADDITIONAL_DATA);
303 }
304 else if (filename.rfind(".xdr") < filename.size())
305 {
306 XdrIO(mesh,true).write("mesh-"+filename);
307 temp_es.write("soln-"+filename,ENCODE,
308 EquationSystems::WRITE_DATA |
309 EquationSystems::WRITE_ADDITIONAL_DATA);
310 }
311 else
312 {
313 libmesh_here();
314 libMesh::err << "Warning: ErrorVector::plot_error currently only"
315 << " supports .gmv, .plt, .xdr/.xda, and .exo/.e (if enabled) output;" << std::endl;
316 libMesh::err << "Could not recognize filename: " << filename
317 << std::endl;
318 }
319 }
320
321 } // namespace libMesh
322