Line data Source code
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
46 0 : ErrorVectorReal ErrorVector::minimum() const
47 : {
48 0 : LOG_SCOPE ("minimum()", "ErrorVector");
49 :
50 0 : const dof_id_type n = cast_int<dof_id_type>(this->size());
51 0 : ErrorVectorReal min = std::numeric_limits<ErrorVectorReal>::max();
52 :
53 0 : for (dof_id_type i=0; i<n; i++)
54 : {
55 : // Only positive (or zero) values in the error vector
56 0 : libmesh_assert_greater_equal ((*this)[i], 0.);
57 0 : if (this->is_active_elem(i))
58 0 : min = std::min (min, (*this)[i]);
59 : }
60 :
61 : // ErrorVectors are for positive values
62 0 : libmesh_assert_greater_equal (min, 0.);
63 :
64 0 : return min;
65 : }
66 :
67 :
68 :
69 14409 : Real ErrorVector::mean() const
70 : {
71 457 : LOG_SCOPE ("mean()", "ErrorVector");
72 :
73 914 : const dof_id_type n = cast_int<dof_id_type>(this->size());
74 :
75 457 : Real the_mean = 0;
76 457 : dof_id_type nnz = 0;
77 :
78 762634 : for (dof_id_type i=0; i<n; i++)
79 748225 : if (this->is_active_elem(i))
80 : {
81 550350 : the_mean += ( static_cast<Real>((*this)[i]) - the_mean ) / (nnz + 1);
82 :
83 15837 : nnz++;
84 : }
85 :
86 14866 : return the_mean;
87 : }
88 :
89 :
90 :
91 :
92 0 : Real ErrorVector::median()
93 : {
94 0 : const dof_id_type n = cast_int<dof_id_type>(this->size());
95 :
96 0 : if (n == 0)
97 0 : return 0.;
98 :
99 :
100 : // Build a StatisticsVector<ErrorVectorReal> containing
101 : // only our active entries and take its mean
102 0 : StatisticsVector<ErrorVectorReal> sv;
103 :
104 0 : sv.reserve (n);
105 :
106 0 : for (dof_id_type i=0; i<n; i++)
107 0 : if (this->is_active_elem(i))
108 0 : sv.push_back((*this)[i]);
109 :
110 0 : return sv.median();
111 : }
112 :
113 :
114 :
115 :
116 0 : Real ErrorVector::median() const
117 : {
118 0 : ErrorVector ev = (*this);
119 :
120 0 : return ev.median();
121 : }
122 :
123 :
124 :
125 :
126 2982 : Real ErrorVector::variance(const Real mean_in) const
127 : {
128 168 : const dof_id_type n = cast_int<dof_id_type>(this->size());
129 :
130 84 : LOG_SCOPE ("variance()", "ErrorVector");
131 :
132 84 : Real the_variance = 0;
133 84 : dof_id_type nnz = 0;
134 :
135 178658 : for (dof_id_type i=0; i<n; i++)
136 175676 : if (this->is_active_elem(i))
137 : {
138 135236 : const Real delta = ( static_cast<Real>((*this)[i]) - mean_in );
139 135236 : the_variance += (delta * delta - the_variance) / (nnz + 1);
140 :
141 3812 : nnz++;
142 : }
143 :
144 3066 : return the_variance;
145 : }
146 :
147 :
148 :
149 :
150 0 : std::vector<dof_id_type> ErrorVector::cut_below(Real cut) const
151 : {
152 0 : LOG_SCOPE ("cut_below()", "ErrorVector");
153 :
154 0 : const dof_id_type n = cast_int<dof_id_type>(this->size());
155 :
156 0 : std::vector<dof_id_type> cut_indices;
157 0 : cut_indices.reserve(n/2); // Arbitrary
158 :
159 0 : for (dof_id_type i=0; i<n; i++)
160 0 : if (this->is_active_elem(i))
161 : {
162 0 : if ((*this)[i] < cut)
163 : {
164 0 : cut_indices.push_back(i);
165 : }
166 : }
167 :
168 0 : return cut_indices;
169 : }
170 :
171 :
172 :
173 :
174 0 : std::vector<dof_id_type> ErrorVector::cut_above(Real cut) const
175 : {
176 0 : LOG_SCOPE ("cut_above()", "ErrorVector");
177 :
178 0 : const dof_id_type n = cast_int<dof_id_type>(this->size());
179 :
180 0 : std::vector<dof_id_type> cut_indices;
181 0 : cut_indices.reserve(n/2); // Arbitrary
182 :
183 0 : for (dof_id_type i=0; i<n; i++)
184 0 : if (this->is_active_elem(i))
185 : {
186 0 : if ((*this)[i] > cut)
187 : {
188 0 : cut_indices.push_back(i);
189 : }
190 : }
191 :
192 0 : return cut_indices;
193 : }
194 :
195 :
196 :
197 923901 : bool ErrorVector::is_active_elem (dof_id_type i) const
198 : {
199 27383 : libmesh_assert_less (i, this->size());
200 :
201 923901 : if (_mesh)
202 : {
203 0 : return _mesh->elem_ptr(i)->active();
204 : }
205 : else
206 951288 : return ((*this)[i] != 0.);
207 : }
208 :
209 :
210 11117 : void ErrorVector::plot_error(const std::string & filename,
211 : const MeshBase & oldmesh,
212 : const std::string & data_type) const
213 : {
214 11505 : std::unique_ptr<MeshBase> meshptr = oldmesh.clone();
215 388 : 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 388 : mesh.allow_renumbering(false);
220 11117 : 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 480398 : for (auto & elem : mesh.element_ptr_range())
226 : {
227 255194 : elem->set_p_refinement_flag(Elem::DO_NOTHING);
228 255194 : elem->set_p_level(0);
229 10341 : }
230 : #endif // LIBMESH_ENABLE_AMR
231 :
232 11893 : EquationSystems temp_es (mesh);
233 : ExplicitSystem & error_system
234 11117 : = temp_es.add_system<ExplicitSystem> ("Error");
235 11505 : error_system.add_variable(data_type, CONSTANT, MONOMIAL);
236 :
237 11117 : temp_es.init();
238 :
239 388 : const DofMap & error_dof_map = error_system.get_dof_map();
240 776 : std::vector<dof_id_type> dof_indices;
241 :
242 178194 : for (const auto & elem : mesh.active_local_element_ptr_range())
243 : {
244 91417 : error_dof_map.dof_indices(elem, dof_indices);
245 :
246 91417 : const dof_id_type elem_id = elem->id();
247 :
248 : //0 for the monomial basis
249 91417 : const dof_id_type solution_index = dof_indices[0];
250 :
251 : // libMesh::out << "elem_number=" << elem_number << std::endl;
252 9608 : 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 101025 : error_system.solution->set(solution_index, (*this)[elem_id]);
257 10341 : }
258 11117 : error_system.solution->close();
259 11117 : 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 22234 : if (mesh.max_elem_id() != mesh.n_elem() ||
265 11117 : mesh.max_node_id() != mesh.n_nodes())
266 : {
267 384 : mesh.allow_renumbering(true);
268 9644 : mesh.renumber_nodes_and_elements();
269 : }
270 :
271 11117 : if (Utility::contains(filename, ".gmv"))
272 : {
273 3946 : GMVIO(mesh).write_discontinuous_gmv(filename,
274 : temp_es, false);
275 : }
276 7171 : else if (Utility::contains(filename, ".plt"))
277 : {
278 0 : TecplotIO (mesh).write_equation_systems
279 0 : (filename, temp_es);
280 : }
281 : #if defined(LIBMESH_HAVE_EXODUS_API) && defined(LIBMESH_HAVE_NEMESIS_API)
282 14342 : else if (Utility::contains(filename, ".nem") ||
283 7373 : Utility::contains(filename, ".n"))
284 : {
285 7575 : Nemesis_IO io(mesh);
286 7171 : io.write(filename);
287 7171 : io.write_element_data(temp_es);
288 6767 : }
289 : #endif
290 : #ifdef LIBMESH_HAVE_EXODUS_API
291 0 : else if (Utility::contains(filename, ".exo") ||
292 0 : Utility::contains(filename, ".e"))
293 : {
294 0 : ExodusII_IO io(mesh);
295 0 : io.write(filename);
296 0 : io.write_element_data(temp_es);
297 0 : }
298 : #endif
299 0 : else if (Utility::contains(filename, ".xda"))
300 : {
301 0 : XdrIO(mesh).write("mesh-"+filename);
302 0 : temp_es.write("soln-"+filename,WRITE,
303 : EquationSystems::WRITE_DATA |
304 : EquationSystems::WRITE_ADDITIONAL_DATA);
305 : }
306 0 : else if (Utility::contains(filename, ".xdr"))
307 : {
308 0 : XdrIO(mesh,true).write("mesh-"+filename);
309 0 : temp_es.write("soln-"+filename,ENCODE,
310 : EquationSystems::WRITE_DATA |
311 : EquationSystems::WRITE_ADDITIONAL_DATA);
312 : }
313 : else
314 : {
315 0 : libmesh_here();
316 0 : libMesh::err << "Warning: ErrorVector::plot_error currently only"
317 0 : << " supports .gmv, .plt, .xdr/.xda, and .exo/.e (if enabled) output;" << std::endl;
318 0 : libMesh::err << "Could not recognize filename: " << filename
319 0 : << std::endl;
320 : }
321 11117 : }
322 :
323 : } // namespace libMesh
|
