Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2026 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 : // Local includes
19 : #include "libmesh/dof_map.h"
20 : #include "libmesh/elem.h"
21 : #include "libmesh/exact_solution.h"
22 : #include "libmesh/equation_systems.h"
23 : #include "libmesh/fe_base.h"
24 : #include "libmesh/function_base.h"
25 : #include "libmesh/mesh_base.h"
26 : #include "libmesh/mesh_function.h"
27 : #include "libmesh/mesh_serializer.h"
28 : #include "libmesh/numeric_vector.h"
29 : #include "libmesh/parallel.h"
30 : #include "libmesh/quadrature.h"
31 : #include "libmesh/wrapped_function.h"
32 : #include "libmesh/wrapped_functor.h"
33 : #include "libmesh/fe_interface.h"
34 : #include "libmesh/raw_accessor.h"
35 : #include "libmesh/tensor_tools.h"
36 : #include "libmesh/enum_norm_type.h"
37 : #include "libmesh/utility.h"
38 :
39 : // C++ Includes
40 : #include <memory>
41 :
42 :
43 : namespace libMesh
44 : {
45 :
46 650 : ExactSolution::ExactSolution(const EquationSystems & es) :
47 290 : _equation_systems(es),
48 290 : _equation_systems_fine(nullptr),
49 1370 : _extra_order(1)
50 : {
51 : // Initialize the _errors data structure which holds all
52 : // the eventual values of the error.
53 1452 : for (auto sys : make_range(_equation_systems.n_systems()))
54 : {
55 : // Reference to the system
56 802 : const System & system = _equation_systems.get_system(sys);
57 :
58 : // The name of the system
59 216 : const std::string & sys_name = system.name();
60 :
61 : // The SystemErrorMap to be inserted
62 432 : ExactSolution::SystemErrorMap sem;
63 :
64 2768 : for (auto var : make_range(system.n_vars()))
65 : {
66 : // The name of this variable
67 1380 : const std::string & var_name = system.variable_name(var);
68 2372 : sem[var_name] = std::vector<Real>(5, 0.);
69 : }
70 :
71 802 : _errors[sys_name] = sem;
72 : }
73 650 : }
74 :
75 :
76 0 : ExactSolution::ExactSolution(ExactSolution &&) = default;
77 580 : ExactSolution::~ExactSolution() = default;
78 :
79 :
80 0 : void ExactSolution::
81 : set_excluded_subdomains(const std::set<subdomain_id_type> & excluded)
82 : {
83 0 : _excluded_subdomains = excluded;
84 0 : }
85 :
86 7 : void ExactSolution::attach_reference_solution (const EquationSystems * es_fine)
87 : {
88 2 : libmesh_assert(es_fine);
89 7 : _equation_systems_fine = es_fine;
90 :
91 : // If we're using a fine grid solution, we're not using exact values
92 5 : _exact_values.clear();
93 5 : _exact_derivs.clear();
94 5 : _exact_hessians.clear();
95 7 : }
96 :
97 :
98 227 : void ExactSolution::attach_exact_value (ValueFunctionPointer fptr)
99 : {
100 66 : libmesh_assert(fptr);
101 :
102 : // Clear out any previous _exact_values entries, then add a new
103 : // entry for each system.
104 227 : _exact_values.clear();
105 :
106 466 : for (auto sys : make_range(_equation_systems.n_systems()))
107 : {
108 239 : const System & system = _equation_systems.get_system(sys);
109 309 : _exact_values.emplace_back(std::make_unique<WrappedFunctor<Number>>(system, fptr, &_equation_systems.parameters));
110 : }
111 :
112 : // If we're using exact values, we're not using a fine grid solution
113 227 : _equation_systems_fine = nullptr;
114 227 : }
115 :
116 :
117 360 : void ExactSolution::attach_exact_values (const std::vector<FunctionBase<Number> *> & f)
118 : {
119 : // Automatically delete any previous _exact_values entries, then add a new
120 : // entry for each system.
121 360 : _exact_values.clear();
122 :
123 720 : for (auto ptr : f)
124 552 : _exact_values.emplace_back(ptr ? std::make_unique<WrappedFunctor<Number>>(*ptr) : nullptr);
125 360 : }
126 :
127 :
128 0 : void ExactSolution::attach_exact_values (const std::vector<FEMFunctionBase<Number> *> & f)
129 : {
130 : // Automatically delete any previous _exact_values entries, then add a new
131 : // entry for each system.
132 0 : _exact_values.clear();
133 :
134 0 : for (auto ptr : f)
135 0 : _exact_values.emplace_back(ptr ? ptr->clone() : nullptr);
136 0 : }
137 :
138 :
139 14 : void ExactSolution::attach_exact_value (unsigned int sys_num,
140 : FunctionBase<Number> * f)
141 : {
142 18 : if (_exact_values.size() <= sys_num)
143 14 : _exact_values.resize(sys_num+1);
144 :
145 14 : if (f)
146 18 : _exact_values[sys_num] = std::make_unique<WrappedFunctor<Number>>(*f);
147 14 : }
148 :
149 :
150 0 : void ExactSolution::attach_exact_value (unsigned int sys_num,
151 : FEMFunctionBase<Number> * f)
152 : {
153 0 : if (_exact_values.size() <= sys_num)
154 0 : _exact_values.resize(sys_num+1);
155 :
156 0 : if (f)
157 0 : _exact_values[sys_num] = f->clone();
158 0 : }
159 :
160 :
161 203 : void ExactSolution::attach_exact_deriv (GradientFunctionPointer gptr)
162 : {
163 58 : libmesh_assert(gptr);
164 :
165 : // Clear out any previous _exact_derivs entries, then add a new
166 : // entry for each system.
167 203 : _exact_derivs.clear();
168 :
169 406 : for (auto sys : make_range(_equation_systems.n_systems()))
170 : {
171 203 : const System & system = _equation_systems.get_system(sys);
172 261 : _exact_derivs.emplace_back(std::make_unique<WrappedFunctor<Gradient>>(system, gptr, &_equation_systems.parameters));
173 : }
174 :
175 : // If we're using exact values, we're not using a fine grid solution
176 203 : _equation_systems_fine = nullptr;
177 203 : }
178 :
179 :
180 360 : void ExactSolution::attach_exact_derivs (const std::vector<FunctionBase<Gradient> *> & g)
181 : {
182 : // Automatically delete any previous _exact_derivs entries, then add a new
183 : // entry for each system.
184 360 : _exact_derivs.clear();
185 :
186 720 : for (auto ptr : g)
187 552 : _exact_derivs.emplace_back(ptr ? std::make_unique<WrappedFunctor<Gradient>>(*ptr) : nullptr);
188 360 : }
189 :
190 :
191 0 : void ExactSolution::attach_exact_derivs (const std::vector<FEMFunctionBase<Gradient> *> & g)
192 : {
193 : // Automatically delete any previous _exact_derivs entries, then add a new
194 : // entry for each system.
195 0 : _exact_derivs.clear();
196 :
197 0 : for (auto ptr : g)
198 0 : _exact_derivs.emplace_back(ptr ? ptr->clone() : nullptr);
199 0 : }
200 :
201 :
202 0 : void ExactSolution::attach_exact_deriv (unsigned int sys_num,
203 : FunctionBase<Gradient> * g)
204 : {
205 0 : if (_exact_derivs.size() <= sys_num)
206 0 : _exact_derivs.resize(sys_num+1);
207 :
208 0 : if (g)
209 0 : _exact_derivs[sys_num] = std::make_unique<WrappedFunctor<Gradient>>(*g);
210 0 : }
211 :
212 :
213 0 : void ExactSolution::attach_exact_deriv (unsigned int sys_num,
214 : FEMFunctionBase<Gradient> * g)
215 : {
216 0 : if (_exact_derivs.size() <= sys_num)
217 0 : _exact_derivs.resize(sys_num+1);
218 :
219 0 : if (g)
220 0 : _exact_derivs[sys_num] = g->clone();
221 0 : }
222 :
223 :
224 42 : void ExactSolution::attach_exact_hessian (HessianFunctionPointer hptr)
225 : {
226 12 : libmesh_assert(hptr);
227 :
228 : // Clear out any previous _exact_hessians entries, then add a new
229 : // entry for each system.
230 42 : _exact_hessians.clear();
231 :
232 84 : for (auto sys : make_range(_equation_systems.n_systems()))
233 : {
234 42 : const System & system = _equation_systems.get_system(sys);
235 54 : _exact_hessians.emplace_back(std::make_unique<WrappedFunctor<Tensor>>(system, hptr, &_equation_systems.parameters));
236 : }
237 :
238 : // If we're using exact values, we're not using a fine grid solution
239 42 : _equation_systems_fine = nullptr;
240 42 : }
241 :
242 :
243 0 : void ExactSolution::attach_exact_hessians (std::vector<FunctionBase<Tensor> *> h)
244 : {
245 : // Automatically delete any previous _exact_hessians entries, then add a new
246 : // entry for each system.
247 0 : _exact_hessians.clear();
248 :
249 0 : for (auto ptr : h)
250 0 : _exact_hessians.emplace_back(ptr ? std::make_unique<WrappedFunctor<Tensor>>(*ptr) : nullptr);
251 0 : }
252 :
253 :
254 0 : void ExactSolution::attach_exact_hessians (std::vector<FEMFunctionBase<Tensor> *> h)
255 : {
256 : // Automatically delete any previous _exact_hessians entries, then add a new
257 : // entry for each system.
258 0 : _exact_hessians.clear();
259 :
260 0 : for (auto ptr : h)
261 0 : _exact_hessians.emplace_back(ptr ? ptr->clone() : nullptr);
262 0 : }
263 :
264 :
265 0 : void ExactSolution::attach_exact_hessian (unsigned int sys_num,
266 : FunctionBase<Tensor> * h)
267 : {
268 0 : if (_exact_hessians.size() <= sys_num)
269 0 : _exact_hessians.resize(sys_num+1);
270 :
271 0 : if (h)
272 0 : _exact_hessians[sys_num] = std::make_unique<WrappedFunctor<Tensor>>(*h);
273 0 : }
274 :
275 :
276 0 : void ExactSolution::attach_exact_hessian (unsigned int sys_num,
277 : FEMFunctionBase<Tensor> * h)
278 : {
279 0 : if (_exact_hessians.size() <= sys_num)
280 0 : _exact_hessians.resize(sys_num+1);
281 :
282 0 : if (h)
283 0 : _exact_hessians[sys_num] = h->clone();
284 0 : }
285 :
286 :
287 12225 : std::vector<Real> & ExactSolution::_check_inputs(std::string_view sys_name,
288 : std::string_view unknown_name)
289 : {
290 : // Return a reference to the proper error entry, or throw an error
291 : // if it doesn't exist.
292 12225 : auto & system_error_map = libmesh_map_find(_errors, sys_name);
293 12225 : return libmesh_map_find(system_error_map, unknown_name);
294 : }
295 :
296 :
297 :
298 2620 : void ExactSolution::compute_error(std::string_view sys_name,
299 : std::string_view unknown_name)
300 : {
301 : // Check the inputs for validity, and get a reference
302 : // to the proper location to store the error
303 1128 : std::vector<Real> & error_vals = this->_check_inputs(sys_name,
304 1492 : unknown_name);
305 :
306 746 : libmesh_assert( _equation_systems.has_system(sys_name) );
307 2620 : const System & sys = _equation_systems.get_system<System>( sys_name );
308 :
309 746 : libmesh_assert( sys.has_variable( unknown_name ) );
310 2620 : switch( FEInterface::field_type(sys.variable_type( unknown_name )) )
311 : {
312 2190 : case TYPE_SCALAR:
313 : {
314 2190 : this->_compute_error<Real>(sys_name,
315 : unknown_name,
316 : error_vals);
317 2190 : break;
318 : }
319 430 : case TYPE_VECTOR:
320 : {
321 430 : this->_compute_error<RealGradient>(sys_name,
322 : unknown_name,
323 : error_vals);
324 430 : break;
325 : }
326 0 : default:
327 0 : libmesh_error_msg("Invalid variable type!");
328 : }
329 2620 : }
330 :
331 :
332 :
333 :
334 :
335 566 : Real ExactSolution::error_norm(std::string_view sys_name,
336 : std::string_view unknown_name,
337 : const FEMNormType & norm)
338 : {
339 : // Check the inputs for validity, and get a reference
340 : // to the proper location to store the error
341 270 : std::vector<Real> & error_vals = this->_check_inputs(sys_name,
342 296 : unknown_name);
343 :
344 148 : libmesh_assert(_equation_systems.has_system(sys_name));
345 148 : libmesh_assert(_equation_systems.get_system(sys_name).has_variable( unknown_name ));
346 566 : const FEType & fe_type = _equation_systems.get_system(sys_name).variable_type(unknown_name);
347 :
348 566 : switch (norm)
349 : {
350 0 : case L2:
351 0 : return std::sqrt(error_vals[0]);
352 0 : case H1:
353 0 : return std::sqrt(error_vals[0] + error_vals[1]);
354 0 : case H2:
355 0 : return std::sqrt(error_vals[0] + error_vals[1] + error_vals[2]);
356 61 : case HCURL:
357 : {
358 61 : libmesh_error_msg_if(FEInterface::field_type(fe_type) == TYPE_SCALAR,
359 : "Cannot compute HCurl error norm of scalar-valued variables!");
360 :
361 61 : return std::sqrt(error_vals[0] + error_vals[5]);
362 : }
363 222 : case HDIV:
364 : {
365 222 : libmesh_error_msg_if(FEInterface::field_type(fe_type) == TYPE_SCALAR,
366 : "Cannot compute HDiv error norm of scalar-valued variables!");
367 :
368 222 : return std::sqrt(error_vals[0] + error_vals[6]);
369 : }
370 0 : case H1_SEMINORM:
371 0 : return std::sqrt(error_vals[1]);
372 0 : case H2_SEMINORM:
373 0 : return std::sqrt(error_vals[2]);
374 61 : case HCURL_SEMINORM:
375 : {
376 61 : libmesh_error_msg_if(FEInterface::field_type(fe_type) == TYPE_SCALAR,
377 : "Cannot compute HCurl error seminorm of scalar-valued variables!");
378 :
379 61 : return std::sqrt(error_vals[5]);
380 : }
381 222 : case HDIV_SEMINORM:
382 : {
383 222 : libmesh_error_msg_if(FEInterface::field_type(fe_type) == TYPE_SCALAR,
384 : "Cannot compute HDiv error seminorm of scalar-valued variables!");
385 :
386 222 : return std::sqrt(error_vals[6]);
387 : }
388 0 : case L1:
389 0 : return error_vals[3];
390 0 : case L_INF:
391 0 : return error_vals[4];
392 :
393 0 : default:
394 0 : libmesh_error_msg("Currently only Sobolev norms/seminorms are supported!");
395 : }
396 : }
397 :
398 :
399 :
400 :
401 :
402 :
403 :
404 4748 : Real ExactSolution::l2_error(std::string_view sys_name,
405 : std::string_view unknown_name)
406 : {
407 :
408 : // Check the inputs for validity, and get a reference
409 : // to the proper location to store the error
410 2040 : std::vector<Real> & error_vals = this->_check_inputs(sys_name,
411 2708 : unknown_name);
412 :
413 : // Return the square root of the first component of the
414 : // computed error.
415 4748 : return std::sqrt(error_vals[0]);
416 : }
417 :
418 :
419 :
420 :
421 :
422 :
423 :
424 0 : Real ExactSolution::l1_error(std::string_view sys_name,
425 : std::string_view unknown_name)
426 : {
427 :
428 : // Check the inputs for validity, and get a reference
429 : // to the proper location to store the error
430 0 : std::vector<Real> & error_vals = this->_check_inputs(sys_name,
431 0 : unknown_name);
432 :
433 : // Return the square root of the first component of the
434 : // computed error.
435 0 : return error_vals[3];
436 : }
437 :
438 :
439 :
440 :
441 :
442 :
443 :
444 0 : Real ExactSolution::l_inf_error(std::string_view sys_name,
445 : std::string_view unknown_name)
446 : {
447 :
448 : // Check the inputs for validity, and get a reference
449 : // to the proper location to store the error
450 0 : std::vector<Real> & error_vals = this->_check_inputs(sys_name,
451 0 : unknown_name);
452 :
453 : // Return the square root of the first component of the
454 : // computed error.
455 0 : return error_vals[4];
456 : }
457 :
458 :
459 :
460 :
461 :
462 :
463 :
464 2947 : Real ExactSolution::h1_error(std::string_view sys_name,
465 : std::string_view unknown_name)
466 : {
467 : // If the user has supplied no exact derivative function, we
468 : // just integrate the H1 norm of the solution; i.e. its
469 : // difference from an "exact solution" of zero.
470 :
471 : // Check the inputs for validity, and get a reference
472 : // to the proper location to store the error
473 1263 : std::vector<Real> & error_vals = this->_check_inputs(sys_name,
474 1684 : unknown_name);
475 :
476 : // Return the square root of the sum of the computed errors.
477 2947 : return std::sqrt(error_vals[0] + error_vals[1]);
478 : }
479 :
480 :
481 61 : Real ExactSolution::hcurl_error(std::string_view sys_name,
482 : std::string_view unknown_name)
483 : {
484 61 : return this->error_norm(sys_name,unknown_name,HCURL);
485 : }
486 :
487 :
488 222 : Real ExactSolution::hdiv_error(std::string_view sys_name,
489 : std::string_view unknown_name)
490 : {
491 222 : return this->error_norm(sys_name,unknown_name,HDIV);
492 : }
493 :
494 :
495 :
496 1344 : Real ExactSolution::h2_error(std::string_view sys_name,
497 : std::string_view unknown_name)
498 : {
499 : // If the user has supplied no exact derivative functions, we
500 : // just integrate the H2 norm of the solution; i.e. its
501 : // difference from an "exact solution" of zero.
502 :
503 : // Check the inputs for validity, and get a reference
504 : // to the proper location to store the error
505 576 : std::vector<Real> & error_vals = this->_check_inputs(sys_name,
506 768 : unknown_name);
507 :
508 : // Return the square root of the sum of the computed errors.
509 1344 : return std::sqrt(error_vals[0] + error_vals[1] + error_vals[2]);
510 : }
511 :
512 :
513 :
514 :
515 :
516 :
517 :
518 : template<typename OutputShape>
519 2620 : void ExactSolution::_compute_error(std::string_view sys_name,
520 : std::string_view unknown_name,
521 : std::vector<Real> & error_vals)
522 : {
523 : // Make sure we aren't "overconfigured"
524 746 : libmesh_assert (!(_exact_values.size() && _equation_systems_fine));
525 :
526 : // We need a communicator.
527 2620 : const Parallel::Communicator & communicator(_equation_systems.comm());
528 :
529 : // This function must be run on all processors at once
530 746 : libmesh_parallel_only(communicator);
531 :
532 : // Get a reference to the system whose error is being computed.
533 : // If we have a fine grid, however, we'll integrate on that instead
534 : // for more accuracy.
535 2620 : const System & computed_system = _equation_systems_fine ?
536 7 : _equation_systems_fine->get_system(sys_name) :
537 2613 : _equation_systems.get_system (sys_name);
538 :
539 4112 : FEMContext context(computed_system);
540 :
541 1492 : const MeshBase & mesh = computed_system.get_mesh();
542 :
543 2620 : const Real time = _equation_systems.get_system(sys_name).time;
544 :
545 1492 : const unsigned int sys_num = computed_system.number();
546 2620 : const unsigned int var = computed_system.variable_number(unknown_name);
547 746 : unsigned int var_component = 0;
548 3564 : for (const auto var_num : make_range(var))
549 : {
550 944 : const auto & var_fe_type = computed_system.variable_type(var_num);
551 944 : const auto var_vec_dim = FEInterface::n_vec_dim(mesh, var_fe_type);
552 944 : var_component += var_vec_dim;
553 : }
554 :
555 : // Prepare a global solution, a serialized mesh, and a MeshFunction
556 : // of the coarse system, if we need them for fine-system integration
557 1874 : std::unique_ptr<MeshFunction> coarse_values;
558 4112 : std::unique_ptr<NumericVector<Number>> comparison_soln =
559 2620 : NumericVector<Number>::build(_equation_systems.comm());
560 : MeshSerializer
561 3366 : serial(const_cast<MeshBase&>(_equation_systems.get_mesh()),
562 2620 : _equation_systems_fine);
563 2620 : if (_equation_systems_fine)
564 : {
565 : const System & comparison_system
566 7 : = _equation_systems.get_system(sys_name);
567 :
568 4 : std::vector<Number> global_soln;
569 7 : comparison_system.update_global_solution(global_soln);
570 9 : comparison_soln->init(comparison_system.solution->size(), true, SERIAL);
571 7 : (*comparison_soln) = global_soln;
572 :
573 10 : coarse_values = std::make_unique<MeshFunction>
574 : (_equation_systems,
575 2 : *comparison_soln,
576 : comparison_system.get_dof_map(),
577 10 : comparison_system.variable_number(unknown_name));
578 7 : coarse_values->init();
579 : }
580 :
581 : // Grab which element dimensions are present in the mesh
582 746 : const std::set<unsigned char> & elem_dims = mesh.elem_dimensions();
583 :
584 : // Initialize any functors we're going to use
585 4933 : for (auto & ev : _exact_values)
586 2313 : if (ev)
587 : {
588 2313 : ev->init();
589 2313 : ev->init_context(context);
590 : }
591 :
592 4811 : for (auto & ed : _exact_derivs)
593 2191 : if (ed)
594 : {
595 2191 : ed->init();
596 2191 : ed->init_context(context);
597 : }
598 :
599 2956 : for (auto & eh : _exact_hessians)
600 336 : if (eh)
601 : {
602 336 : eh->init();
603 336 : eh->init_context(context);
604 : }
605 :
606 : // If we have *no* functors we intend to use (because we're using a
607 : // fine system, or because our exact solution is zero and we're just
608 : // computing norms in an outdated way) then let our FE objects know
609 : // we don't actually need anything from them, so they don't think
610 : // we've just invoked them in a deprecated "compute everything"
611 : // fashion.
612 2718 : if (_exact_values.empty() && _exact_derivs.empty() &&
613 98 : _exact_hessians.empty())
614 686 : for (auto dim : elem_dims)
615 686 : for (auto v : make_range(computed_system.n_vars()))
616 : {
617 : FEAbstract * fe;
618 98 : context.get_element_fe(v, fe, dim);
619 98 : fe->get_nothing();
620 : }
621 :
622 : // Get a reference to the dofmap and mesh for that system
623 746 : const DofMap & computed_dof_map = computed_system.get_dof_map();
624 :
625 : // Zero the error before summation
626 : // 0 - sum of square of function error (L2)
627 : // 1 - sum of square of gradient error (H1 semi)
628 : // 2 - sum of square of Hessian error (H2 semi)
629 : // 3 - sum of sqrt(square of function error) (L1)
630 : // 4 - max of sqrt(square of function error) (Linfty)
631 : // 5 - sum of square of curl error (HCurl semi)
632 : // 6 - sum of square of div error (HDiv semi)
633 4494 : error_vals = std::vector<Real>(7, 0.);
634 :
635 : // Construct Quadrature rule based on default quadrature order
636 746 : const FEType & fe_type = computed_dof_map.variable_type(var);
637 2620 : const auto field_type = FEInterface::field_type(fe_type);
638 :
639 2620 : unsigned int n_vec_dim = FEInterface::n_vec_dim( mesh, fe_type );
640 :
641 : // FIXME: MeshFunction needs to be updated to support vector-valued
642 : // elements before we can use a reference solution.
643 2620 : if ((n_vec_dim > 1) && _equation_systems_fine)
644 : {
645 0 : libMesh::err << "Error calculation using reference solution not yet\n"
646 0 : << "supported for vector-valued elements."
647 0 : << std::endl;
648 0 : libmesh_not_implemented();
649 : }
650 :
651 :
652 : // Allow space for dims 0-3, even if we don't use them all
653 4112 : std::vector<std::unique_ptr<FEGenericBase<OutputShape>>> fe_ptrs(4);
654 4112 : std::vector<std::unique_ptr<QBase>> q_rules(4);
655 :
656 : // Prepare finite elements for each dimension present in the mesh
657 5247 : for (const auto dim : elem_dims)
658 : {
659 : // Build a quadrature rule.
660 2627 : q_rules[dim] = fe_type.default_quadrature_rule (dim, _extra_order);
661 :
662 : // Disallow rules with negative weights. That will use more
663 : // quadrature points, but we're going to be taking square roots
664 : // of element integral results here!
665 2627 : q_rules[dim]->allow_rules_with_negative_weights = false;
666 :
667 : // Construct finite element object
668 3758 : fe_ptrs[dim] = FEGenericBase<OutputShape>::build(dim, fe_type);
669 :
670 : // Attach quadrature rule to FE object
671 4123 : fe_ptrs[dim]->attach_quadrature_rule (q_rules[dim].get());
672 : }
673 :
674 : // The global degree of freedom indices associated
675 : // with the local degrees of freedom.
676 746 : std::vector<dof_id_type> dof_indices;
677 :
678 :
679 : //
680 : // Begin the loop over the elements
681 : //
682 : // TODO: this ought to be threaded (and using subordinate
683 : // MeshFunction objects in each thread rather than a single
684 : // master)
685 787560 : for (const auto & elem : mesh.active_local_element_ptr_range())
686 : {
687 : // Skip this element if it is in a subdomain excluded by the user.
688 516715 : const subdomain_id_type elem_subid = elem->subdomain_id();
689 125746 : if (_excluded_subdomains.count(elem_subid))
690 0 : continue;
691 :
692 : // The spatial dimension of the current Elem. FEs and other data
693 : // are indexed on dim.
694 516715 : const unsigned int dim = elem->dim();
695 :
696 : // If the variable is not active on this subdomain, don't bother
697 516715 : if (!computed_system.variable(var).active_on_subdomain(elem_subid))
698 0 : continue;
699 :
700 : /* If the variable is active, then we're going to restrict the
701 : MeshFunction evaluations to the current element subdomain.
702 : This is for cases such as mixed dimension meshes where we want
703 : to restrict the calculation to one particular domain. */
704 251492 : std::set<subdomain_id_type> subdomain_id;
705 390969 : subdomain_id.insert(elem_subid);
706 :
707 516715 : FEGenericBase<OutputShape> * fe = fe_ptrs[dim].get();
708 251492 : QBase * qrule = q_rules[dim].get();
709 125746 : libmesh_assert(fe);
710 125746 : libmesh_assert(qrule);
711 :
712 : // The Jacobian*weight at the quadrature points.
713 516715 : const std::vector<Real> & JxW = fe->get_JxW();
714 :
715 : // The value of the shape functions at the quadrature points
716 : // i.e. phi(i) = phi_values[i][qp]
717 125746 : const std::vector<std::vector<OutputShape>> & phi_values = fe->get_phi();
718 :
719 : // The value of the shape function gradients at the quadrature points
720 : const std::vector<std::vector<typename FEGenericBase<OutputShape>::OutputGradient>> &
721 125746 : dphi_values = fe->get_dphi();
722 :
723 : // The value of the shape function curls at the quadrature points
724 : // Only computed for vector-valued elements
725 125746 : const std::vector<std::vector<typename FEGenericBase<OutputShape>::OutputShape>> * curl_values = nullptr;
726 :
727 : // The value of the shape function divergences at the quadrature points
728 : // Only computed for vector-valued elements
729 125746 : const std::vector<std::vector<typename FEGenericBase<OutputShape>::OutputDivergence>> * div_values = nullptr;
730 :
731 516715 : if (field_type == TYPE_VECTOR)
732 : {
733 140577 : curl_values = &fe->get_curl_phi();
734 140577 : div_values = &fe->get_div_phi();
735 : }
736 :
737 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
738 : // The value of the shape function second derivatives at the quadrature points
739 : // Not computed for vector-valued elements
740 : const std::vector<std::vector<typename FEGenericBase<OutputShape>::OutputTensor>> *
741 125746 : d2phi_values = nullptr;
742 :
743 516715 : if (field_type != TYPE_VECTOR)
744 78887 : d2phi_values = &fe->get_d2phi();
745 : #endif
746 :
747 : // The XYZ locations (in physical space) of the quadrature points
748 377063 : const std::vector<Point> & q_point = fe->get_xyz();
749 :
750 : // reinitialize the element-specific data
751 : // for the current element
752 516715 : fe->reinit (elem);
753 :
754 516715 : context.pre_fe_reinit(computed_system, elem);
755 516715 : context.elem_fe_reinit();
756 :
757 : // Get the local to global degree of freedom maps
758 516715 : computed_dof_map.dof_indices (elem, dof_indices, var);
759 :
760 : // The number of quadrature points
761 125746 : const unsigned int n_qp = qrule->n_points();
762 :
763 : // The number of shape functions
764 : const unsigned int n_sf =
765 251492 : cast_int<unsigned int>(dof_indices.size());
766 :
767 : //
768 : // Begin the loop over the Quadrature points.
769 : //
770 8650986 : for (unsigned int qp=0; qp<n_qp; qp++)
771 : {
772 : // Real u_h = 0.;
773 : // RealGradient grad_u_h;
774 :
775 1978529 : typename FEGenericBase<OutputShape>::OutputNumber u_h(0.);
776 :
777 2671719 : typename FEGenericBase<OutputShape>::OutputNumberGradient grad_u_h;
778 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
779 3957058 : typename FEGenericBase<OutputShape>::OutputNumberTensor grad2_u_h;
780 : #endif
781 1978529 : typename FEGenericBase<OutputShape>::OutputNumber curl_u_h(0.0);
782 1978529 : typename FEGenericBase<OutputShape>::OutputNumberDivergence div_u_h = 0.0;
783 :
784 : // Compute solution values at the current
785 : // quadrature point. This requires a sum
786 : // over all the shape functions evaluated
787 : // at the quadrature point.
788 78751143 : for (unsigned int i=0; i<n_sf; i++)
789 : {
790 : // Values from current solution.
791 111579127 : u_h += phi_values[i][qp]*computed_system.current_solution (dof_indices[i]);
792 87961917 : grad_u_h += dphi_values[i][qp]*computed_system.current_solution (dof_indices[i]);
793 70616872 : if (field_type == TYPE_VECTOR)
794 : {
795 28766514 : curl_u_h += (*curl_values)[i][qp]*computed_system.current_solution (dof_indices[i]);
796 39839394 : div_u_h += (*div_values)[i][qp]*computed_system.current_solution (dof_indices[i]);
797 : }
798 : else
799 : {
800 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
801 94773596 : grad2_u_h += (*d2phi_values)[i][qp]*computed_system.current_solution (dof_indices[i]);
802 : #endif
803 : }
804 : }
805 :
806 : // Compute the value of the error at this quadrature point
807 5839018 : typename FEGenericBase<OutputShape>::OutputNumber exact_val(0);
808 1978529 : RawAccessor<typename FEGenericBase<OutputShape>::OutputNumber> exact_val_accessor( exact_val, n_vec_dim );
809 10112800 : if (_exact_values.size() > sys_num && _exact_values[sys_num])
810 : {
811 19613358 : for (unsigned int c = 0; c < n_vec_dim; c++)
812 12416147 : exact_val_accessor(c) =
813 : _exact_values[sys_num]->
814 18315535 : component(context, var_component+c, q_point[qp], time);
815 : }
816 937060 : else if (_equation_systems_fine)
817 : {
818 : // FIXME: Needs to be updated for vector-valued elements
819 281880 : DenseVector<Number> output(1);
820 469800 : (*coarse_values)(q_point[qp],time,output,&subdomain_id);
821 375840 : exact_val = output(0);
822 : }
823 4549203 : const typename FEGenericBase<OutputShape>::OutputNumber val_error = u_h - exact_val;
824 :
825 : // Add the squares of the error to each contribution
826 1978529 : Real error_sq = TensorTools::norm_sq(val_error);
827 8134271 : error_vals[0] += JxW[qp]*error_sq;
828 :
829 8134271 : Real norm = sqrt(error_sq);
830 8134271 : error_vals[3] += JxW[qp]*norm;
831 :
832 8134271 : if (error_vals[4]<norm) { error_vals[4] = norm; }
833 :
834 : // Compute the value of the error in the gradient at this
835 : // quadrature point
836 2671719 : typename FEGenericBase<OutputShape>::OutputNumberGradient exact_grad;
837 1978529 : RawAccessor<typename FEGenericBase<OutputShape>::OutputNumberGradient> exact_grad_accessor( exact_grad, LIBMESH_DIM );
838 10112800 : if (_exact_derivs.size() > sys_num && _exact_derivs[sys_num])
839 : {
840 19537590 : for (unsigned int c = 0; c < n_vec_dim; c++)
841 49513052 : for (unsigned int d = 0; d < LIBMESH_DIM; d++)
842 45955026 : exact_grad_accessor(d + c*LIBMESH_DIM) =
843 : _exact_derivs[sys_num]->
844 54775263 : component(context, var_component+c, q_point[qp], time)(d);
845 : }
846 974944 : else if (_equation_systems_fine)
847 : {
848 : // FIXME: Needs to be updated for vector-valued elements
849 375840 : std::vector<Gradient> output(1);
850 469800 : coarse_values->gradient(q_point[qp],time,output,&subdomain_id);
851 375840 : exact_grad = output[0];
852 : }
853 :
854 4966972 : const typename FEGenericBase<OutputShape>::OutputNumberGradient grad_error = grad_u_h - exact_grad;
855 :
856 11714863 : error_vals[1] += JxW[qp]*grad_error.norm_sq();
857 :
858 :
859 8134271 : if (field_type == TYPE_VECTOR)
860 : {
861 : // Compute the value of the error in the curl at this
862 : // quadrature point
863 693190 : typename FEGenericBase<OutputShape>::OutputNumber exact_curl(0.0);
864 3681633 : if (_exact_derivs.size() > sys_num && _exact_derivs[sys_num])
865 : {
866 2988443 : exact_curl = TensorTools::curl_from_grad( exact_grad );
867 : }
868 0 : else if (_equation_systems_fine)
869 : {
870 : // FIXME: Need to implement curl for MeshFunction and support reference
871 : // solution for vector-valued elements
872 : }
873 :
874 693190 : const typename FEGenericBase<OutputShape>::OutputNumber curl_error = curl_u_h - exact_curl;
875 :
876 3681633 : error_vals[5] += JxW[qp]*TensorTools::norm_sq(curl_error);
877 :
878 : // Compute the value of the error in the divergence at this
879 : // quadrature point
880 693190 : typename FEGenericBase<OutputShape>::OutputNumberDivergence exact_div = 0.0;
881 3681633 : if (_exact_derivs.size() > sys_num && _exact_derivs[sys_num])
882 : {
883 2988443 : exact_div = TensorTools::div_from_grad( exact_grad );
884 : }
885 0 : else if (_equation_systems_fine)
886 : {
887 : // FIXME: Need to implement div for MeshFunction and support reference
888 : // solution for vector-valued elements
889 : }
890 :
891 2295253 : const typename FEGenericBase<OutputShape>::OutputNumberDivergence div_error = div_u_h - exact_div;
892 :
893 4374823 : error_vals[6] += JxW[qp]*TensorTools::norm_sq(div_error);
894 : }
895 :
896 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
897 : // Compute the value of the error in the hessian at this
898 : // quadrature point
899 3957058 : typename FEGenericBase<OutputShape>::OutputNumberTensor exact_hess;
900 1978529 : RawAccessor<typename FEGenericBase<OutputShape>::OutputNumberTensor> exact_hess_accessor( exact_hess, dim );
901 10112800 : if (_exact_hessians.size() > sys_num && _exact_hessians[sys_num])
902 : {
903 : //FIXME: This needs to be implemented to support rank 3 tensors
904 : // which can't happen until type_n_tensor is fully implemented
905 : // and a RawAccessor<TypeNTensor> is fully implemented
906 561220 : if (field_type == TYPE_VECTOR)
907 0 : libmesh_not_implemented();
908 :
909 1122440 : for (unsigned int c = 0; c < n_vec_dim; c++)
910 1682892 : for (unsigned int d = 0; d < dim; d++)
911 3364248 : for (unsigned int e =0; e < dim; e++)
912 2797820 : exact_hess_accessor(d + e*dim + c*dim*dim) =
913 : _exact_hessians[sys_num]->
914 3353064 : component(context, var_component+c, q_point[qp], time)(d,e);
915 :
916 : // FIXME: operator- is not currently implemented for TypeNTensor
917 561220 : const typename FEGenericBase<OutputShape>::OutputNumberTensor grad2_error = grad2_u_h - exact_hess;
918 1122440 : error_vals[2] += JxW[qp]*grad2_error.norm_sq();
919 : }
920 7573051 : else if (_equation_systems_fine)
921 : {
922 : // FIXME: Needs to be updated for vector-valued elements
923 469800 : std::vector<Tensor> output(1);
924 469800 : coarse_values->hessian(q_point[qp],time,output,&subdomain_id);
925 375840 : exact_hess = output[0];
926 :
927 : // FIXME: operator- is not currently implemented for TypeNTensor
928 375840 : const typename FEGenericBase<OutputShape>::OutputNumberTensor grad2_error = grad2_u_h - exact_hess;
929 751680 : error_vals[2] += JxW[qp]*grad2_error.norm_sq();
930 : }
931 : #endif
932 :
933 : } // end qp loop
934 : } // end element loop
935 :
936 : // Add up the error values on all processors, except for the L-infty
937 : // norm, for which the maximum is computed.
938 2620 : Real l_infty_norm = error_vals[4];
939 1874 : communicator.max(l_infty_norm);
940 2620 : communicator.sum(error_vals);
941 2620 : error_vals[4] = l_infty_norm;
942 2620 : }
943 :
944 : // Explicit instantiations of templated member functions
945 : template LIBMESH_EXPORT void ExactSolution::_compute_error<Real>(std::string_view, std::string_view, std::vector<Real> &);
946 : template LIBMESH_EXPORT void ExactSolution::_compute_error<RealGradient>(std::string_view, std::string_view, std::vector<Real> &);
947 :
948 : } // namespace libMesh
|
