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