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 : // Local Includes
20 : #include "libmesh/dof_map.h"
21 : #include "libmesh/elem.h"
22 : #include "libmesh/equation_systems.h"
23 : #include "libmesh/error_vector.h"
24 : #include "libmesh/fe.h"
25 : #include "libmesh/libmesh_common.h"
26 : #include "libmesh/libmesh_logging.h"
27 : #include "libmesh/mesh_base.h"
28 : #include "libmesh/mesh_refinement.h"
29 : #include "libmesh/numeric_vector.h"
30 : #include "libmesh/quadrature.h"
31 : #include "libmesh/system.h"
32 : #include "libmesh/uniform_refinement_estimator.h"
33 : #include "libmesh/partitioner.h"
34 : #include "libmesh/tensor_tools.h"
35 : #include "libmesh/enum_error_estimator_type.h"
36 : #include "libmesh/enum_norm_type.h"
37 : #include "libmesh/int_range.h"
38 :
39 : // C++ includes
40 : #include <algorithm> // for std::fill
41 : #include <sstream>
42 : #include <cstdlib> // *must* precede <cmath> for proper std:abs() on PGI, Sun Studio CC
43 : #include <cmath> // for sqrt
44 : #include <memory>
45 :
46 : #ifdef LIBMESH_ENABLE_AMR
47 :
48 : namespace libMesh
49 : {
50 :
51 : //-----------------------------------------------------------------
52 : // ErrorEstimator implementations
53 :
54 0 : UniformRefinementEstimator::UniformRefinementEstimator() :
55 : ErrorEstimator(),
56 0 : number_h_refinements(1),
57 0 : number_p_refinements(0),
58 0 : _extra_order(1)
59 : {
60 0 : error_norm = H1;
61 0 : }
62 :
63 :
64 :
65 0 : ErrorEstimatorType UniformRefinementEstimator::type() const
66 : {
67 0 : return UNIFORM_REFINEMENT;
68 : }
69 :
70 :
71 0 : void UniformRefinementEstimator::estimate_error (const System & _system,
72 : ErrorVector & error_per_cell,
73 : const NumericVector<Number> * solution_vector,
74 : bool estimate_parent_error)
75 : {
76 0 : LOG_SCOPE("estimate_error()", "UniformRefinementEstimator");
77 0 : std::map<const System *, const NumericVector<Number> *> solution_vectors;
78 0 : solution_vectors[&_system] = solution_vector;
79 0 : this->_estimate_error (nullptr,
80 : &_system,
81 : &error_per_cell,
82 : nullptr,
83 : nullptr,
84 : &solution_vectors,
85 0 : estimate_parent_error);
86 0 : }
87 :
88 0 : void UniformRefinementEstimator::estimate_errors (const EquationSystems & _es,
89 : ErrorVector & error_per_cell,
90 : const std::map<const System *, SystemNorm> & error_norms,
91 : const std::map<const System *, const NumericVector<Number> *> * solution_vectors,
92 : bool estimate_parent_error)
93 : {
94 0 : LOG_SCOPE("estimate_errors()", "UniformRefinementEstimator");
95 0 : this->_estimate_error (&_es,
96 : nullptr,
97 : &error_per_cell,
98 : nullptr,
99 : &error_norms,
100 : solution_vectors,
101 0 : estimate_parent_error);
102 0 : }
103 :
104 0 : void UniformRefinementEstimator::estimate_errors (const EquationSystems & _es,
105 : ErrorMap & errors_per_cell,
106 : const std::map<const System *, const NumericVector<Number> *> * solution_vectors,
107 : bool estimate_parent_error)
108 : {
109 0 : LOG_SCOPE("estimate_errors()", "UniformRefinementEstimator");
110 0 : this->_estimate_error (&_es,
111 : nullptr,
112 : nullptr,
113 : &errors_per_cell,
114 : nullptr,
115 : solution_vectors,
116 0 : estimate_parent_error);
117 0 : }
118 :
119 0 : void UniformRefinementEstimator::_estimate_error (const EquationSystems * _es,
120 : const System * _system,
121 : ErrorVector * error_per_cell,
122 : ErrorMap * errors_per_cell,
123 : const std::map<const System *, SystemNorm> * _error_norms,
124 : const std::map<const System *, const NumericVector<Number> *> * solution_vectors,
125 : bool)
126 : {
127 : // Get a vector of the Systems we're going to work on,
128 : // and set up a error_norms map if necessary
129 0 : std::vector<System *> system_list;
130 0 : auto error_norms = std::make_unique<std::map<const System *, SystemNorm>>();
131 :
132 0 : if (_es)
133 : {
134 0 : libmesh_assert(!_system);
135 0 : libmesh_assert(_es->n_systems());
136 0 : _system = &(_es->get_system(0));
137 0 : libmesh_assert_equal_to (&(_system->get_equation_systems()), _es);
138 :
139 0 : libmesh_assert(_es->n_systems());
140 0 : for (auto i : make_range(_es->n_systems()))
141 : // We have to break the rules here, because we can't refine a const System
142 0 : system_list.push_back(const_cast<System *>(&(_es->get_system(i))));
143 :
144 : // If we're computing one vector, we need to know how to scale
145 : // each variable's contributions to it.
146 0 : if (_error_norms)
147 : {
148 0 : libmesh_assert(!errors_per_cell);
149 : }
150 : else
151 : // If we're computing many vectors, we just need to know which
152 : // variables to skip
153 : {
154 0 : libmesh_assert (errors_per_cell);
155 :
156 0 : _error_norms = error_norms.get();
157 :
158 0 : for (auto i : make_range(_es->n_systems()))
159 : {
160 0 : const System & sys = _es->get_system(i);
161 0 : unsigned int n_vars = sys.n_vars();
162 :
163 0 : std::vector<Real> weights(n_vars, 0.0);
164 0 : for (unsigned int v = 0; v != n_vars; ++v)
165 : {
166 0 : if (!errors_per_cell->count(std::make_pair(&sys, v)))
167 0 : continue;
168 :
169 0 : weights[v] = 1.0;
170 : }
171 0 : (*error_norms)[&sys] =
172 0 : SystemNorm(std::vector<FEMNormType>(n_vars, error_norm.type(0)),
173 0 : weights);
174 : }
175 : }
176 : }
177 : else
178 : {
179 0 : libmesh_assert(_system);
180 : // We have to break the rules here, because we can't refine a const System
181 0 : system_list.push_back(const_cast<System *>(_system));
182 :
183 0 : libmesh_assert(!_error_norms);
184 0 : (*error_norms)[_system] = error_norm;
185 0 : _error_norms = error_norms.get();
186 : }
187 :
188 : // An EquationSystems reference will be convenient.
189 : // We have to break the rules here, because we can't refine a const System
190 : EquationSystems & es =
191 0 : const_cast<EquationSystems &>(_system->get_equation_systems());
192 :
193 : // The current mesh
194 0 : MeshBase & mesh = es.get_mesh();
195 :
196 : // The dimensionality of the mesh
197 0 : const unsigned int dim = mesh.mesh_dimension();
198 :
199 : // Resize the error_per_cell vectors to be
200 : // the number of elements, initialize them to 0.
201 0 : if (error_per_cell)
202 : {
203 0 : error_per_cell->clear();
204 0 : error_per_cell->resize (mesh.max_elem_id(), 0.);
205 : }
206 : else
207 : {
208 0 : libmesh_assert(errors_per_cell);
209 0 : for (const auto & pr : *errors_per_cell)
210 : {
211 0 : ErrorVector * e = pr.second.get();
212 0 : e->clear();
213 0 : e->resize(mesh.max_elem_id(), 0.);
214 : }
215 : }
216 :
217 : // We'll want to back up all coarse grid vectors
218 0 : std::vector<std::map<std::string, std::unique_ptr<NumericVector<Number>>>> coarse_vectors(system_list.size());
219 0 : std::vector<std::unique_ptr<NumericVector<Number>>> coarse_solutions(system_list.size());
220 0 : std::vector<std::unique_ptr<NumericVector<Number>>> coarse_local_solutions(system_list.size());
221 : // And make copies of projected solutions
222 0 : std::vector<std::unique_ptr<NumericVector<Number>>> projected_solutions(system_list.size());
223 :
224 : // And we'll need to temporarily change solution projection settings
225 0 : std::vector<bool> old_projection_settings(system_list.size());
226 :
227 : // And it'll be best to avoid any repartitioning
228 0 : std::unique_ptr<Partitioner> old_partitioner(mesh.partitioner().release());
229 :
230 : // And we can't allow any renumbering
231 0 : const bool old_renumbering_setting = mesh.allow_renumbering();
232 0 : mesh.allow_renumbering(false);
233 :
234 0 : for (auto i : index_range(system_list))
235 : {
236 0 : System & system = *system_list[i];
237 :
238 : // Check for valid error_norms
239 0 : libmesh_assert (_error_norms->find(&system) !=
240 : _error_norms->end());
241 :
242 : // Back up the solution vector
243 0 : coarse_solutions[i] = system.solution->clone();
244 0 : coarse_local_solutions[i] = system.current_local_solution->clone();
245 :
246 : // Back up all other coarse grid vectors
247 0 : for (System::vectors_iterator vec = system.vectors_begin(); vec !=
248 0 : system.vectors_end(); ++vec)
249 : {
250 : // The (string) name of this vector
251 0 : const std::string & var_name = vec->first;
252 :
253 0 : coarse_vectors[i][var_name] = vec->second->clone();
254 : }
255 :
256 : // Use a non-standard solution vector if necessary
257 0 : if (solution_vectors &&
258 0 : solution_vectors->find(&system) != solution_vectors->end() &&
259 0 : solution_vectors->find(&system)->second &&
260 0 : solution_vectors->find(&system)->second != system.solution.get())
261 : {
262 : NumericVector<Number> * newsol =
263 : const_cast<NumericVector<Number> *>
264 0 : (solution_vectors->find(&system)->second);
265 0 : newsol->swap(*system.solution);
266 0 : system.update();
267 : }
268 :
269 : // Make sure the solution is projected when we refine the mesh
270 0 : old_projection_settings[i] = system.project_solution_on_reinit();
271 0 : system.project_solution_on_reinit() = true;
272 : }
273 :
274 : // Find the number of coarse mesh elements, to make it possible
275 : // to find correct coarse elem ids later
276 0 : const dof_id_type max_coarse_elem_id = mesh.max_elem_id();
277 : #ifndef NDEBUG
278 : // n_coarse_elem is only used in an assertion later so
279 : // avoid declaring it unless asserts are active.
280 0 : const dof_id_type n_coarse_elem = mesh.n_elem();
281 : #endif
282 :
283 : // Uniformly refine the mesh
284 0 : MeshRefinement mesh_refinement(mesh);
285 :
286 0 : libmesh_assert (number_h_refinements > 0 || number_p_refinements > 0);
287 :
288 : // FIXME: this may break if there is more than one System
289 : // on this mesh but estimate_error was still called instead of
290 : // estimate_errors
291 0 : for (unsigned int i = 0; i != number_h_refinements; ++i)
292 : {
293 0 : mesh_refinement.uniformly_refine(1);
294 0 : es.reinit();
295 : }
296 :
297 0 : for (unsigned int i = 0; i != number_p_refinements; ++i)
298 : {
299 0 : mesh_refinement.uniformly_p_refine(1);
300 0 : es.reinit();
301 : }
302 :
303 0 : for (auto i : index_range(system_list))
304 : {
305 0 : System & system = *system_list[i];
306 :
307 : // Copy the projected coarse grid solutions, which will be
308 : // overwritten by solve()
309 0 : projected_solutions[i] = NumericVector<Number>::build(system.comm());
310 0 : projected_solutions[i]->init(system.solution->size(),
311 0 : system.solution->local_size(),
312 0 : system.get_dof_map().get_send_list(),
313 0 : true, GHOSTED);
314 0 : system.solution->localize(*projected_solutions[i],
315 0 : system.get_dof_map().get_send_list());
316 : }
317 :
318 : // Are we doing a forward or an adjoint solve?
319 0 : bool solve_adjoint = false;
320 0 : if (solution_vectors)
321 : {
322 0 : System * sys = system_list[0];
323 0 : libmesh_assert (solution_vectors->find(sys) !=
324 : solution_vectors->end());
325 0 : const NumericVector<Number> * vec = solution_vectors->find(sys)->second;
326 0 : for (auto j : make_range(sys->n_qois()))
327 : {
328 0 : std::ostringstream adjoint_name;
329 0 : adjoint_name << "adjoint_solution" << j;
330 :
331 0 : if (vec == sys->request_vector(adjoint_name.str()))
332 : {
333 0 : solve_adjoint = true;
334 0 : break;
335 : }
336 0 : }
337 : }
338 :
339 : // Get the uniformly refined solution.
340 :
341 0 : if (_es)
342 : {
343 : // Even if we had a decent preconditioner, valid matrix etc. before
344 : // refinement, we don't any more.
345 0 : for (auto i : make_range(_es->n_systems()))
346 0 : es.get_system(i).disable_cache();
347 :
348 : // No specified vectors == forward solve
349 0 : if (!solution_vectors)
350 0 : es.solve();
351 : else
352 : {
353 0 : libmesh_assert_equal_to (solution_vectors->size(), es.n_systems());
354 0 : libmesh_assert (solution_vectors->find(system_list[0]) !=
355 : solution_vectors->end());
356 0 : libmesh_assert(solve_adjoint ||
357 : (solution_vectors->find(system_list[0])->second ==
358 : system_list[0]->solution.get()) ||
359 : !solution_vectors->find(system_list[0])->second);
360 :
361 : #ifdef DEBUG
362 0 : for (const auto & sys : system_list)
363 : {
364 0 : libmesh_assert (solution_vectors->find(sys) !=
365 : solution_vectors->end());
366 0 : const NumericVector<Number> * vec = solution_vectors->find(sys)->second;
367 0 : if (solve_adjoint)
368 : {
369 0 : bool found_vec = false;
370 0 : for (auto j : make_range(sys->n_qois()))
371 : {
372 0 : std::ostringstream adjoint_name;
373 0 : adjoint_name << "adjoint_solution" << j;
374 :
375 0 : if (vec == sys->request_vector(adjoint_name.str()))
376 : {
377 0 : found_vec = true;
378 0 : break;
379 : }
380 : }
381 0 : libmesh_assert(found_vec);
382 : }
383 : else
384 0 : libmesh_assert(vec == sys->solution.get() || !vec);
385 : }
386 : #endif
387 :
388 0 : if (solve_adjoint)
389 : {
390 : std::vector<unsigned int> adjs(system_list.size(),
391 0 : libMesh::invalid_uint);
392 : // Set up proper initial guesses
393 0 : for (auto i : index_range(system_list))
394 : {
395 0 : System * sys = system_list[i];
396 0 : libmesh_assert (solution_vectors->find(sys) !=
397 : solution_vectors->end());
398 0 : const NumericVector<Number> * vec = solution_vectors->find(sys)->second;
399 0 : for (auto j : make_range(sys->n_qois()))
400 : {
401 0 : std::ostringstream adjoint_name;
402 0 : adjoint_name << "adjoint_solution" << j;
403 :
404 0 : if (vec == sys->request_vector(adjoint_name.str()))
405 : {
406 0 : adjs[i] = j;
407 0 : break;
408 : }
409 0 : }
410 0 : libmesh_assert_not_equal_to (adjs[i], libMesh::invalid_uint);
411 0 : sys->get_adjoint_solution(adjs[i]) = *sys->solution;
412 : }
413 :
414 0 : es.adjoint_solve();
415 :
416 : // Put the adjoint_solution into solution for
417 : // comparisons
418 0 : for (auto i : index_range(system_list))
419 : {
420 0 : system_list[i]->get_adjoint_solution(adjs[i]).swap(*system_list[i]->solution);
421 0 : system_list[i]->update();
422 : }
423 : }
424 : else
425 0 : es.solve();
426 : }
427 : }
428 : else
429 : {
430 0 : System * sys = system_list[0];
431 :
432 : // Even if we had a decent preconditioner, valid matrix etc. before
433 : // refinement, we don't any more.
434 0 : sys->disable_cache();
435 :
436 : // No specified vectors == forward solve
437 0 : if (!solution_vectors)
438 0 : sys->solve();
439 : else
440 : {
441 0 : libmesh_assert (solution_vectors->find(sys) !=
442 : solution_vectors->end());
443 :
444 0 : const NumericVector<Number> * vec = solution_vectors->find(sys)->second;
445 :
446 0 : libmesh_assert(solve_adjoint ||
447 : (solution_vectors->find(sys)->second ==
448 : sys->solution.get()) ||
449 : !solution_vectors->find(sys)->second);
450 :
451 0 : if (solve_adjoint)
452 : {
453 0 : unsigned int adj = libMesh::invalid_uint;
454 0 : for (unsigned int j=0, n_qois = sys->n_qois();
455 0 : j != n_qois; ++j)
456 : {
457 0 : std::ostringstream adjoint_name;
458 0 : adjoint_name << "adjoint_solution" << j;
459 :
460 0 : if (vec == sys->request_vector(adjoint_name.str()))
461 : {
462 0 : adj = j;
463 0 : break;
464 : }
465 0 : }
466 0 : libmesh_assert_not_equal_to (adj, libMesh::invalid_uint);
467 :
468 : // Set up proper initial guess
469 0 : sys->get_adjoint_solution(adj) = *sys->solution;
470 0 : sys->adjoint_solve();
471 : // Put the adjoint_solution into solution for
472 : // comparisons
473 0 : sys->get_adjoint_solution(adj).swap(*sys->solution);
474 0 : sys->update();
475 : }
476 : else
477 0 : sys->solve();
478 : }
479 : }
480 :
481 : // Get the error in the uniformly refined solution(s).
482 0 : for (auto sysnum : index_range(system_list))
483 : {
484 0 : System & system = *system_list[sysnum];
485 :
486 0 : unsigned int n_vars = system.n_vars();
487 :
488 0 : DofMap & dof_map = system.get_dof_map();
489 :
490 : const SystemNorm & system_i_norm =
491 0 : _error_norms->find(&system)->second;
492 :
493 0 : NumericVector<Number> * projected_solution = projected_solutions[sysnum].get();
494 :
495 : // Loop over all the variables in the system
496 0 : for (unsigned int var=0; var<n_vars; var++)
497 : {
498 : // Get the error vector to fill for this system and variable
499 0 : ErrorVector * err_vec = error_per_cell;
500 0 : if (!err_vec)
501 : {
502 0 : libmesh_assert(errors_per_cell);
503 : err_vec =
504 0 : (*errors_per_cell)[std::make_pair(&system,var)].get();
505 : }
506 :
507 : // The type of finite element to use for this variable
508 0 : const FEType & fe_type = dof_map.variable_type (var);
509 :
510 : // Finite element object for each fine element
511 0 : std::unique_ptr<FEBase> fe (FEBase::build (dim, fe_type));
512 :
513 : // Build and attach an appropriate quadrature rule
514 0 : std::unique_ptr<QBase> qrule = fe_type.default_quadrature_rule(dim, _extra_order);
515 0 : fe->attach_quadrature_rule (qrule.get());
516 :
517 0 : const std::vector<Real> & JxW = fe->get_JxW();
518 0 : const std::vector<std::vector<Real>> & phi = fe->get_phi();
519 : const std::vector<std::vector<RealGradient>> & dphi =
520 0 : fe->get_dphi();
521 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
522 : const std::vector<std::vector<RealTensor>> & d2phi =
523 0 : fe->get_d2phi();
524 : #endif
525 :
526 : // The global DOF indices for the fine element
527 0 : std::vector<dof_id_type> dof_indices;
528 :
529 : // Iterate over all the active elements in the fine mesh
530 : // that live on this processor.
531 0 : for (const auto & elem : mesh.active_local_element_ptr_range())
532 : {
533 : // Find the element id for the corresponding coarse grid element
534 0 : const Elem * coarse = elem;
535 0 : dof_id_type e_id = coarse->id();
536 0 : while (e_id >= max_coarse_elem_id)
537 : {
538 0 : libmesh_assert (coarse->parent());
539 0 : coarse = coarse->parent();
540 0 : e_id = coarse->id();
541 : }
542 :
543 0 : Real L2normsq = 0., H1seminormsq = 0.;
544 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
545 0 : Real H2seminormsq = 0.;
546 : #endif
547 :
548 : // reinitialize the element-specific data
549 : // for the current element
550 0 : fe->reinit (elem);
551 :
552 : // Get the local to global degree of freedom maps
553 0 : dof_map.dof_indices (elem, dof_indices, var);
554 :
555 : // The number of quadrature points
556 0 : const unsigned int n_qp = qrule->n_points();
557 :
558 : // The number of shape functions
559 : const unsigned int n_sf =
560 0 : cast_int<unsigned int>(dof_indices.size());
561 :
562 : //
563 : // Begin the loop over the Quadrature points.
564 : //
565 0 : for (unsigned int qp=0; qp<n_qp; qp++)
566 : {
567 0 : Number u_fine = 0., u_coarse = 0.;
568 :
569 0 : Gradient grad_u_fine, grad_u_coarse;
570 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
571 0 : Tensor grad2_u_fine, grad2_u_coarse;
572 : #endif
573 :
574 : // Compute solution values at the current
575 : // quadrature point. This requires a sum
576 : // over all the shape functions evaluated
577 : // at the quadrature point.
578 0 : for (unsigned int i=0; i<n_sf; i++)
579 : {
580 0 : u_fine += phi[i][qp]*system.current_solution (dof_indices[i]);
581 0 : u_coarse += phi[i][qp]*(*projected_solution) (dof_indices[i]);
582 0 : grad_u_fine += dphi[i][qp]*system.current_solution (dof_indices[i]);
583 0 : grad_u_coarse += dphi[i][qp]*(*projected_solution) (dof_indices[i]);
584 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
585 0 : grad2_u_fine += d2phi[i][qp]*system.current_solution (dof_indices[i]);
586 0 : grad2_u_coarse += d2phi[i][qp]*(*projected_solution) (dof_indices[i]);
587 : #endif
588 : }
589 :
590 : // Compute the value of the error at this quadrature point
591 0 : const Number val_error = u_fine - u_coarse;
592 :
593 : // Add the squares of the error to each contribution
594 0 : if (system_i_norm.type(var) == L2 ||
595 0 : system_i_norm.type(var) == H1 ||
596 0 : system_i_norm.type(var) == H2)
597 : {
598 0 : L2normsq += JxW[qp] * system_i_norm.weight_sq(var) *
599 0 : TensorTools::norm_sq(val_error);
600 0 : libmesh_assert_greater_equal (L2normsq, 0.);
601 : }
602 :
603 :
604 : // Compute the value of the error in the gradient at this
605 : // quadrature point
606 0 : if (system_i_norm.type(var) == H1 ||
607 0 : system_i_norm.type(var) == H2 ||
608 0 : system_i_norm.type(var) == H1_SEMINORM)
609 : {
610 0 : Gradient grad_error = grad_u_fine - grad_u_coarse;
611 :
612 0 : H1seminormsq += JxW[qp] * system_i_norm.weight_sq(var) *
613 0 : grad_error.norm_sq();
614 0 : libmesh_assert_greater_equal (H1seminormsq, 0.);
615 : }
616 :
617 : // Compute the value of the error in the hessian at this
618 : // quadrature point
619 0 : if (system_i_norm.type(var) == H2 ||
620 0 : system_i_norm.type(var) == H2_SEMINORM)
621 : {
622 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
623 0 : Tensor grad2_error = grad2_u_fine - grad2_u_coarse;
624 :
625 0 : H2seminormsq += JxW[qp] * system_i_norm.weight_sq(var) *
626 0 : grad2_error.norm_sq();
627 0 : libmesh_assert_greater_equal (H2seminormsq, 0.);
628 : #else
629 : libmesh_error_msg
630 : ("libMesh was not configured with --enable-second");
631 : #endif
632 : }
633 : } // end qp loop
634 :
635 0 : if (system_i_norm.type(var) == L2 ||
636 0 : system_i_norm.type(var) == H1 ||
637 0 : system_i_norm.type(var) == H2)
638 0 : (*err_vec)[e_id] +=
639 0 : static_cast<ErrorVectorReal>(L2normsq);
640 0 : if (system_i_norm.type(var) == H1 ||
641 0 : system_i_norm.type(var) == H2 ||
642 0 : system_i_norm.type(var) == H1_SEMINORM)
643 0 : (*err_vec)[e_id] +=
644 0 : static_cast<ErrorVectorReal>(H1seminormsq);
645 :
646 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
647 0 : if (system_i_norm.type(var) == H2 ||
648 0 : system_i_norm.type(var) == H2_SEMINORM)
649 0 : (*err_vec)[e_id] +=
650 0 : static_cast<ErrorVectorReal>(H2seminormsq);
651 : #endif
652 0 : } // End loop over active local elements
653 0 : } // End loop over variables
654 :
655 : // Don't bother projecting the solution; we'll restore from backup
656 : // after coarsening
657 0 : system.project_solution_on_reinit() = false;
658 : }
659 :
660 :
661 : // Uniformly coarsen the mesh, without projecting the solution
662 0 : libmesh_assert (number_h_refinements > 0 || number_p_refinements > 0);
663 :
664 0 : for (unsigned int i = 0; i != number_h_refinements; ++i)
665 : {
666 0 : mesh_refinement.uniformly_coarsen(1);
667 : // FIXME - should the reinits here be necessary? - RHS
668 0 : es.reinit();
669 : }
670 :
671 0 : for (unsigned int i = 0; i != number_p_refinements; ++i)
672 : {
673 0 : mesh_refinement.uniformly_p_coarsen(1);
674 0 : es.reinit();
675 : }
676 :
677 : // We should be back where we started
678 0 : libmesh_assert_equal_to (n_coarse_elem, mesh.n_elem());
679 :
680 : // Each processor has now computed the error contributions
681 : // for its local elements. We need to sum the vector
682 : // and then take the square-root of each component. Note
683 : // that we only need to sum if we are running on multiple
684 : // processors, and we only need to take the square-root
685 : // if the value is nonzero. There will in general be many
686 : // zeros for the inactive elements.
687 :
688 0 : if (error_per_cell)
689 : {
690 : // First sum the vector of estimated error values
691 0 : this->reduce_error(*error_per_cell, es.comm());
692 :
693 : // Compute the square-root of each component.
694 0 : LOG_SCOPE("std::sqrt()", "UniformRefinementEstimator");
695 0 : for (auto & val : *error_per_cell)
696 0 : if (val != 0.)
697 0 : val = std::sqrt(val);
698 : }
699 : else
700 : {
701 0 : for (const auto & pr : *errors_per_cell)
702 : {
703 0 : ErrorVector & e = *(pr.second);
704 : // First sum the vector of estimated error values
705 0 : this->reduce_error(e, es.comm());
706 :
707 : // Compute the square-root of each component.
708 0 : LOG_SCOPE("std::sqrt()", "UniformRefinementEstimator");
709 0 : for (auto & val : e)
710 0 : if (val != 0.)
711 0 : val = std::sqrt(val);
712 : }
713 : }
714 :
715 : // Restore old solutions and clean up the heap
716 0 : for (auto i : index_range(system_list))
717 : {
718 0 : System & system = *system_list[i];
719 :
720 0 : system.project_solution_on_reinit() = old_projection_settings[i];
721 :
722 : // Restore the coarse solution vectors and delete their copies
723 0 : *system.solution = *coarse_solutions[i];
724 0 : *system.current_local_solution = *coarse_local_solutions[i];
725 :
726 0 : for (System::vectors_iterator vec = system.vectors_begin(); vec !=
727 0 : system.vectors_end(); ++vec)
728 : {
729 : // The (string) name of this vector
730 0 : const std::string & var_name = vec->first;
731 :
732 0 : system.get_vector(var_name) = *coarse_vectors[i][var_name];
733 :
734 0 : coarse_vectors[i][var_name]->clear();
735 : }
736 : }
737 :
738 : // Restore old partitioner and renumbering settings
739 0 : mesh.partitioner().reset(old_partitioner.release());
740 0 : mesh.allow_renumbering(old_renumbering_setting);
741 0 : }
742 :
743 : } // namespace libMesh
744 :
745 : #endif // #ifdef LIBMESH_ENABLE_AMR
|
