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 :
20 : // Local includes
21 : #include "libmesh/dof_map.h"
22 : #include "libmesh/equation_systems.h"
23 : #include "libmesh/int_range.h"
24 : #include "libmesh/libmesh_logging.h"
25 : #include "libmesh/mesh_base.h"
26 : #include "libmesh/numeric_vector.h"
27 : #include "libmesh/parameter_vector.h"
28 : #include "libmesh/point.h" // For point_value
29 : #include "libmesh/point_locator_base.h" // For point_value
30 : #include "libmesh/qoi_set.h"
31 : #include "libmesh/enum_to_string.h"
32 : #include "libmesh/sparse_matrix.h"
33 : #include "libmesh/system.h"
34 : #include "libmesh/system_norm.h"
35 : #include "libmesh/utility.h"
36 : #include "libmesh/elem.h"
37 : #include "libmesh/fe_type.h"
38 : #include "libmesh/parallel_fe_type.h"
39 : #include "libmesh/fe_interface.h"
40 : #include "libmesh/fe_compute_data.h"
41 : #include "libmesh/static_condensation.h"
42 :
43 : // includes for calculate_norm, point_*
44 : #include "libmesh/fe_base.h"
45 : #include "libmesh/fe_interface.h"
46 : #include "libmesh/parallel.h"
47 : #include "libmesh/parallel_algebra.h"
48 : #include "libmesh/quadrature.h"
49 : #include "libmesh/tensor_value.h"
50 : #include "libmesh/vector_value.h"
51 : #include "libmesh/tensor_tools.h"
52 : #include "libmesh/enum_norm_type.h"
53 : #include "libmesh/enum_fe_family.h"
54 :
55 : // C++ includes
56 : #include <sstream> // for std::ostringstream
57 :
58 : namespace libMesh
59 : {
60 :
61 :
62 : // ------------------------------------------------------------
63 : // System implementation
64 237355 : System::System (EquationSystems & es,
65 : const std::string & name_in,
66 237355 : const unsigned int number_in) :
67 :
68 : ParallelObject (es),
69 223807 : assemble_before_solve (true),
70 223807 : use_fixed_solution (false),
71 223807 : extra_quadrature_order (0),
72 223807 : solution (NumericVector<Number>::build(this->comm())),
73 223807 : current_local_solution (NumericVector<Number>::build(this->comm())),
74 223807 : time (0.),
75 223807 : qoi (0),
76 223807 : qoi_error_estimates (0),
77 223807 : _init_system_function (nullptr),
78 223807 : _init_system_object (nullptr),
79 223807 : _assemble_system_function (nullptr),
80 223807 : _assemble_system_object (nullptr),
81 223807 : _constrain_system_function (nullptr),
82 223807 : _constrain_system_object (nullptr),
83 223807 : _qoi_evaluate_function (nullptr),
84 223807 : _qoi_evaluate_object (nullptr),
85 223807 : _qoi_evaluate_derivative_function (nullptr),
86 223807 : _qoi_evaluate_derivative_object (nullptr),
87 223807 : _dof_map (std::make_unique<DofMap>(number_in, es.get_mesh())),
88 223807 : _equation_systems (es),
89 244129 : _mesh (es.get_mesh()),
90 223807 : _sys_name (name_in),
91 223807 : _sys_number (number_in),
92 223807 : _active (true),
93 223807 : _matrices_initialized (false),
94 223807 : _solution_projection (true),
95 223807 : _basic_system_only (false),
96 223807 : _is_initialized (false),
97 223807 : _identify_variable_groups (true),
98 223807 : _additional_data_written (false),
99 223807 : adjoint_already_solved (false),
100 223807 : _hide_output (false),
101 223807 : project_with_constraints (true),
102 223807 : _prefer_hash_table_matrix_assembly(false),
103 223807 : _require_sparsity_pattern (false),
104 467936 : _prefix_with_name (false)
105 : {
106 467936 : if (libMesh::on_command_line("--solver-system-names"))
107 0 : this->prefix_with_name(true);
108 691743 : if (libMesh::on_command_line("--" + name_in + "-static-condensation"))
109 280 : this->create_static_condensation();
110 237355 : }
111 :
112 :
113 :
114 672621 : System::~System ()
115 : {
116 6774 : libmesh_exceptionless_assert (!libMesh::closed());
117 889654 : }
118 :
119 :
120 :
121 1377591 : dof_id_type System::n_dofs() const
122 : {
123 1377591 : return _dof_map->n_dofs();
124 : }
125 :
126 :
127 :
128 40540 : dof_id_type System::n_constrained_dofs() const
129 : {
130 : #ifdef LIBMESH_ENABLE_CONSTRAINTS
131 :
132 40540 : return _dof_map->n_constrained_dofs();
133 :
134 : #else
135 :
136 : return 0;
137 :
138 : #endif
139 : }
140 :
141 :
142 :
143 35042 : dof_id_type System::n_local_constrained_dofs() const
144 : {
145 : #ifdef LIBMESH_ENABLE_CONSTRAINTS
146 :
147 35042 : return _dof_map->n_local_constrained_dofs();
148 :
149 : #else
150 :
151 : return 0;
152 :
153 : #endif
154 : }
155 :
156 :
157 :
158 1365958 : dof_id_type System::n_local_dofs() const
159 : {
160 1365958 : return _dof_map->n_local_dofs();
161 : }
162 :
163 :
164 :
165 1083326068 : Number System::current_solution (const dof_id_type global_dof_number) const
166 : {
167 : // Check the sizes
168 100501773 : libmesh_assert_less (global_dof_number, _dof_map->n_dofs());
169 100501773 : libmesh_assert_less (global_dof_number, current_local_solution->size());
170 :
171 1083326068 : return (*current_local_solution)(global_dof_number);
172 : }
173 :
174 :
175 :
176 849 : void System::clear ()
177 : {
178 823 : _variables.clear();
179 26 : _variable_numbers.clear();
180 849 : _dof_map->clear ();
181 849 : solution->clear ();
182 849 : current_local_solution->clear ();
183 :
184 : // clear any user-added vectors
185 26 : _vectors.clear();
186 26 : _vector_projections.clear();
187 26 : _vector_is_adjoint.clear();
188 849 : _is_initialized = false;
189 :
190 : // clear any user-added matrices
191 26 : _matrices.clear();
192 849 : _matrices_initialized = false;
193 849 : }
194 :
195 :
196 :
197 493 : void System::init ()
198 : {
199 : // Calling init() twice on the same system currently works evil
200 : // magic, whether done directly or via EquationSystems::read()
201 14 : libmesh_assert(!this->is_initialized());
202 :
203 493 : this->reinit_mesh();
204 493 : }
205 :
206 :
207 :
208 237380 : void System::init_data ()
209 : {
210 6776 : parallel_object_only();
211 :
212 13552 : MeshBase & mesh = this->get_mesh();
213 :
214 : // Add all variable groups to our underlying DofMap
215 6776 : unsigned int n_dof_map_vg = _dof_map->n_variable_groups();
216 485020 : for (auto vg : make_range(this->n_variable_groups()))
217 : {
218 7066 : const VariableGroup & group = this->variable_group(vg);
219 247640 : if (vg < n_dof_map_vg)
220 12 : libmesh_assert(group == _dof_map->variable_group(vg));
221 : else
222 487374 : _dof_map->add_variable_group(group);
223 : }
224 :
225 : // Distribute the degrees of freedom on the mesh
226 237380 : auto total_dofs = _dof_map->distribute_dofs (mesh);
227 :
228 : // Throw an error if the total number of DOFs is not capable of
229 : // being indexed by our solution vector.
230 237357 : auto max_allowed_id = solution->max_allowed_id();
231 237357 : libmesh_error_msg_if(total_dofs > max_allowed_id,
232 : "Cannot allocate a NumericVector with " << total_dofs << " degrees of freedom. "
233 : "The vector can only index up to " << max_allowed_id << " entries.");
234 :
235 : // Recreate any user or internal constraints
236 237357 : this->reinit_constraints();
237 :
238 : // Even if there weren't any constraint changes,
239 : // reinit_constraints() did prepare_send_list() for us.
240 :
241 : // Now finally after dof distribution and construction of any
242 : // possible constraints, we may init any static condensation
243 : // data
244 237286 : _dof_map->reinit_static_condensation();
245 :
246 : // Resize the solution conformal to the current mesh
247 237286 : solution->init (this->n_dofs(), this->n_local_dofs(), false, PARALLEL);
248 :
249 : // Resize the current_local_solution for the current mesh
250 : #ifdef LIBMESH_ENABLE_GHOSTED
251 244058 : current_local_solution->init (this->n_dofs(), this->n_local_dofs(),
252 : _dof_map->get_send_list(), /*fast=*/false,
253 13544 : GHOSTED);
254 : #else
255 : current_local_solution->init (this->n_dofs(), false, SERIAL);
256 : #endif
257 :
258 : // from now on, adding additional vectors or variables can't be done
259 : // without immediately initializing them
260 237286 : _is_initialized = true;
261 :
262 : // initialize & zero other vectors, if necessary
263 279793 : for (auto & [vec_name, vec] : _vectors)
264 : {
265 1292 : libmesh_ignore(vec_name); // spurious warning from old gcc
266 42507 : const ParallelType type = vec->type();
267 :
268 42507 : if (type == GHOSTED)
269 : {
270 : #ifdef LIBMESH_ENABLE_GHOSTED
271 5824 : vec->init (this->n_dofs(), this->n_local_dofs(),
272 : _dof_map->get_send_list(), /*fast=*/false,
273 328 : GHOSTED);
274 : #else
275 : libmesh_error_msg("Cannot initialize ghosted vectors when they are not enabled.");
276 : #endif
277 : }
278 36847 : else if (type == SERIAL)
279 : {
280 0 : vec->init (this->n_dofs(), false, type);
281 : }
282 : else
283 : {
284 1128 : libmesh_assert_equal_to(type, PARALLEL);
285 36847 : vec->init (this->n_dofs(), this->n_local_dofs(), false, type);
286 : }
287 : }
288 :
289 : // Add matrices
290 237286 : this->add_matrices();
291 :
292 : // Clear any existing matrices
293 256295 : for (auto & pr : _matrices)
294 19009 : pr.second->clear();
295 :
296 : // Initialize the matrices for the system
297 237286 : if (!_basic_system_only)
298 237286 : this->init_matrices();
299 237286 : }
300 :
301 237380 : void System::reinit_mesh ()
302 : {
303 6776 : parallel_object_only();
304 :
305 : // First initialize any required data:
306 : // either only the basic System data
307 237380 : if (_basic_system_only)
308 0 : System::init_data();
309 : // or all the derived class' data too
310 : else
311 237380 : this->init_data();
312 :
313 : // If no variables have been added to this system
314 : // don't do anything
315 237286 : if (!this->n_vars())
316 24 : return;
317 :
318 : // Then call the user-provided initialization function
319 236443 : this->user_initialization();
320 :
321 : }
322 :
323 237286 : void System::init_matrices ()
324 : {
325 6772 : parallel_object_only();
326 :
327 : // No matrices to init
328 237286 : if (_matrices.empty())
329 : {
330 : // any future matrices to be added will need their own
331 : // initialization
332 218916 : _matrices_initialized = true;
333 :
334 218916 : return;
335 : }
336 :
337 : // Check for quick return in case the first matrix
338 : // (and by extension all the matrices) has already
339 : // been initialized
340 18370 : if (_matrices.begin()->second->initialized())
341 : {
342 0 : libmesh_assert(_matrices_initialized);
343 0 : return;
344 : }
345 :
346 18370 : _matrices_initialized = true;
347 :
348 : // Tell the matrices about the dof map, and vice versa
349 37379 : for (auto & pr : _matrices)
350 : {
351 550 : SparseMatrix<Number> & m = *(pr.second);
352 550 : libmesh_assert (!m.initialized());
353 :
354 : // We want to allow repeated init() on systems, but we don't
355 : // want to attach the same matrix to the DofMap twice
356 19009 : if (!this->get_dof_map().is_attached(m))
357 19009 : this->get_dof_map().attach_matrix(m);
358 :
359 : // If the user has already explicitly requested that this matrix use a hash table, then we
360 : // always honor that
361 : const bool use_hash =
362 19559 : pr.second->use_hash_table() ||
363 18943 : (this->_prefer_hash_table_matrix_assembly && pr.second->supports_hash_table());
364 19009 : pr.second->use_hash_table(use_hash);
365 : // Make this call after we've determined whether the matrix is using a hash table
366 19009 : if (pr.second->require_sparsity_pattern())
367 18527 : this->_require_sparsity_pattern = true;
368 : }
369 :
370 : // Compute the sparsity pattern for the current
371 : // mesh and DOF distribution. This also updates
372 : // additional matrices, \p DofMap now knows them
373 18370 : if (this->_require_sparsity_pattern)
374 17954 : this->get_dof_map().compute_sparsity(this->get_mesh());
375 :
376 : // Initialize matrices and set to zero
377 37379 : for (auto & [name, mat] : _matrices)
378 : {
379 19009 : mat->init(_matrix_types[name]);
380 19009 : mat->zero();
381 : }
382 : }
383 :
384 :
385 :
386 24142 : void System::restrict_vectors ()
387 : {
388 800 : parallel_object_only();
389 :
390 : #ifdef LIBMESH_ENABLE_AMR
391 : // Restrict the _vectors on the coarsened cells
392 73475 : for (auto & [vec_name, vec] : _vectors)
393 : {
394 1830 : NumericVector<Number> * v = vec.get();
395 :
396 49333 : if (_vector_projections[vec_name])
397 : {
398 20375 : this->project_vector (*v, this->vector_is_adjoint(vec_name));
399 : }
400 : else
401 : {
402 28958 : const ParallelType type = vec->type();
403 :
404 28958 : if (type == GHOSTED)
405 : {
406 : #ifdef LIBMESH_ENABLE_GHOSTED
407 0 : vec->init (this->n_dofs(), this->n_local_dofs(),
408 : _dof_map->get_send_list(), /*fast=*/false,
409 0 : GHOSTED);
410 : #else
411 : libmesh_error_msg("Cannot initialize ghosted vectors when they are not enabled.");
412 : #endif
413 : }
414 : else
415 28958 : vec->init (this->n_dofs(), this->n_local_dofs(), false, type);
416 : }
417 : }
418 :
419 800 : const std::vector<dof_id_type> & send_list = _dof_map->get_send_list ();
420 :
421 : // Restrict the solution on the coarsened cells
422 24142 : if (_solution_projection)
423 23290 : this->project_vector (*solution);
424 : // Or at least make sure the solution vector is the correct size
425 : else
426 852 : solution->init (this->n_dofs(), this->n_local_dofs(), true, PARALLEL);
427 :
428 : #ifdef LIBMESH_ENABLE_GHOSTED
429 24942 : current_local_solution->init(this->n_dofs(),
430 : this->n_local_dofs(), send_list,
431 1600 : false, GHOSTED);
432 : #else
433 : current_local_solution->init(this->n_dofs());
434 : #endif
435 :
436 24142 : if (_solution_projection)
437 24066 : solution->localize (*current_local_solution, send_list);
438 :
439 : #endif // LIBMESH_ENABLE_AMR
440 24142 : }
441 :
442 :
443 :
444 24142 : void System::prolong_vectors ()
445 : {
446 : #ifdef LIBMESH_ENABLE_AMR
447 : // Currently project_vector handles both restriction and prolongation
448 24142 : this->restrict_vectors();
449 : #endif
450 24142 : }
451 :
452 :
453 :
454 23000 : void System::reinit ()
455 : {
456 766 : parallel_object_only();
457 :
458 : // project_vector handles vector initialization now
459 766 : libmesh_assert_equal_to (solution->size(), current_local_solution->size());
460 :
461 : // Make sure our static condensation dof map is up-to-date before we init any
462 : // static condensation matrices
463 23000 : this->get_dof_map().reinit_static_condensation();
464 :
465 23000 : if (!_matrices.empty() && !_basic_system_only)
466 : {
467 : // Clear the matrices
468 39814 : for (auto & pr : _matrices)
469 : {
470 19977 : pr.second->clear();
471 19977 : pr.second->attach_dof_map(this->get_dof_map());
472 : }
473 :
474 19837 : if (this->_require_sparsity_pattern)
475 : {
476 : // Clear the sparsity pattern
477 19697 : this->get_dof_map().clear_sparsity();
478 :
479 : // Compute the sparsity pattern for the current
480 : // mesh and DOF distribution. This also updates
481 : // additional matrices, \p DofMap now knows them
482 19697 : this->get_dof_map().compute_sparsity (this->get_mesh());
483 : }
484 :
485 : // Initialize matrices and set to zero
486 39814 : for (auto & pr : _matrices)
487 : {
488 19977 : pr.second->init();
489 19977 : pr.second->zero();
490 : }
491 : }
492 23000 : }
493 :
494 :
495 261779 : void System::reinit_constraints()
496 : {
497 7582 : parallel_object_only();
498 :
499 : #ifdef LIBMESH_ENABLE_CONSTRAINTS
500 261779 : get_dof_map().create_dof_constraints(_mesh, this->time);
501 261779 : user_constrain();
502 261779 : get_dof_map().process_constraints(_mesh);
503 515836 : if (libMesh::on_command_line ("--print-constraints"))
504 0 : get_dof_map().print_dof_constraints(libMesh::out);
505 : #endif
506 261708 : get_dof_map().prepare_send_list();
507 261708 : }
508 :
509 :
510 1434109 : void System::update ()
511 : {
512 36816 : parallel_object_only();
513 :
514 36816 : libmesh_assert(solution->closed());
515 :
516 36816 : const std::vector<dof_id_type> & send_list = _dof_map->get_send_list ();
517 :
518 : // Check sizes
519 36816 : libmesh_assert_equal_to (current_local_solution->size(), solution->size());
520 : // More processors than elements => empty send_list
521 : // libmesh_assert (!send_list.empty());
522 36816 : libmesh_assert_less_equal (send_list.size(), solution->size());
523 :
524 : // Create current_local_solution from solution. This will
525 : // put a local copy of solution into current_local_solution.
526 : // Only the necessary values (specified by the send_list)
527 : // are copied to minimize communication
528 1470109 : solution->localize (*current_local_solution, send_list);
529 1434109 : }
530 :
531 :
532 :
533 22935 : void System::re_update ()
534 : {
535 766 : parallel_object_only();
536 :
537 : // If this system is empty... don't do anything!
538 22935 : if (!this->n_vars())
539 8 : return;
540 :
541 758 : const std::vector<dof_id_type> & send_list = this->get_dof_map().get_send_list ();
542 :
543 : // Check sizes
544 758 : libmesh_assert_equal_to (current_local_solution->size(), solution->size());
545 : // Not true with ghosted vectors
546 : // libmesh_assert_equal_to (current_local_solution->local_size(), solution->size());
547 : // libmesh_assert (!send_list.empty());
548 758 : libmesh_assert_less_equal (send_list.size(), solution->size());
549 :
550 : // Create current_local_solution from solution. This will
551 : // put a local copy of solution into current_local_solution.
552 23411 : solution->localize (*current_local_solution, send_list);
553 : }
554 :
555 :
556 :
557 0 : void System::restrict_solve_to (const SystemSubset * subset,
558 : const SubsetSolveMode /*subset_solve_mode*/)
559 : {
560 0 : if (subset != nullptr)
561 0 : libmesh_not_implemented();
562 0 : }
563 :
564 :
565 :
566 39219 : void System::assemble ()
567 : {
568 : // Log how long the user's assembly code takes
569 2704 : LOG_SCOPE("assemble()", "System");
570 :
571 : // Call the user-specified assembly function
572 39219 : this->user_assembly();
573 39219 : }
574 :
575 :
576 :
577 0 : void System::assemble_qoi (const QoISet & qoi_indices)
578 : {
579 : // Log how long the user's assembly code takes
580 0 : LOG_SCOPE("assemble_qoi()", "System");
581 :
582 : // Call the user-specified quantity of interest function
583 0 : this->user_QOI(qoi_indices);
584 0 : }
585 :
586 :
587 :
588 0 : void System::assemble_qoi_derivative(const QoISet & qoi_indices,
589 : bool include_liftfunc,
590 : bool apply_constraints)
591 : {
592 : // Log how long the user's assembly code takes
593 0 : LOG_SCOPE("assemble_qoi_derivative()", "System");
594 :
595 : // Call the user-specified quantity of interest function
596 0 : this->user_QOI_derivative(qoi_indices, include_liftfunc,
597 0 : apply_constraints);
598 0 : }
599 :
600 :
601 :
602 0 : void System::qoi_parameter_sensitivity (const QoISet & qoi_indices,
603 : const ParameterVector & parameters_vec,
604 : SensitivityData & sensitivities)
605 : {
606 : // Forward sensitivities are more efficient for Nq > Np
607 0 : if (qoi_indices.size(*this) > parameters_vec.size())
608 0 : forward_qoi_parameter_sensitivity(qoi_indices, parameters_vec, sensitivities);
609 : // Adjoint sensitivities are more efficient for Np > Nq,
610 : // and an adjoint may be more reusable than a forward
611 : // solution sensitivity in the Np == Nq case.
612 : else
613 0 : adjoint_qoi_parameter_sensitivity(qoi_indices, parameters_vec, sensitivities);
614 0 : }
615 :
616 :
617 :
618 0 : bool System::compare (const System & other_system,
619 : const Real threshold,
620 : const bool verbose) const
621 : {
622 : // we do not care for matrices, but for vectors
623 0 : libmesh_assert (_is_initialized);
624 0 : libmesh_assert (other_system._is_initialized);
625 :
626 0 : if (verbose)
627 : {
628 0 : libMesh::out << " Systems \"" << _sys_name << "\"" << std::endl;
629 0 : libMesh::out << " comparing matrices not supported." << std::endl;
630 0 : libMesh::out << " comparing names...";
631 : }
632 :
633 : // compare the name: 0 means identical
634 0 : const int name_result = _sys_name.compare(other_system.name());
635 0 : if (verbose)
636 : {
637 0 : if (name_result == 0)
638 0 : libMesh::out << " identical." << std::endl;
639 : else
640 0 : libMesh::out << " names not identical." << std::endl;
641 0 : libMesh::out << " comparing solution vector...";
642 : }
643 :
644 :
645 : // compare the solution: -1 means identical
646 0 : const int solu_result = solution->compare (*other_system.solution.get(),
647 0 : threshold);
648 :
649 0 : if (verbose)
650 : {
651 0 : if (solu_result == -1)
652 0 : libMesh::out << " identical up to threshold." << std::endl;
653 : else
654 0 : libMesh::out << " first difference occurred at index = "
655 0 : << solu_result << "." << std::endl;
656 : }
657 :
658 :
659 : // safety check, whether we handle at least the same number
660 : // of vectors
661 0 : std::vector<int> ov_result;
662 :
663 0 : if (this->n_vectors() != other_system.n_vectors())
664 : {
665 0 : if (verbose)
666 : {
667 0 : libMesh::out << " Fatal difference. This system handles "
668 0 : << this->n_vectors() << " add'l vectors," << std::endl
669 0 : << " while the other system handles "
670 0 : << other_system.n_vectors()
671 0 : << " add'l vectors." << std::endl
672 0 : << " Aborting comparison." << std::endl;
673 : }
674 0 : return false;
675 : }
676 0 : else if (this->n_vectors() == 0)
677 : {
678 : // there are no additional vectors...
679 0 : ov_result.clear ();
680 : }
681 : else
682 : {
683 : // compare other vectors
684 0 : for (auto & [vec_name, vec] : _vectors)
685 : {
686 0 : if (verbose)
687 0 : libMesh::out << " comparing vector \""
688 0 : << vec_name << "\" ...";
689 :
690 : // assume they have the same name
691 : const NumericVector<Number> & other_system_vector =
692 0 : other_system.get_vector(vec_name);
693 :
694 0 : ov_result.push_back(vec->compare(other_system_vector, threshold));
695 :
696 0 : if (verbose)
697 : {
698 0 : if (ov_result[ov_result.size()-1] == -1)
699 0 : libMesh::out << " identical up to threshold." << std::endl;
700 : else
701 0 : libMesh::out << " first difference occurred at" << std::endl
702 0 : << " index = " << ov_result[ov_result.size()-1] << "." << std::endl;
703 : }
704 : }
705 : } // finished comparing additional vectors
706 :
707 :
708 : bool overall_result;
709 :
710 : // sum up the results
711 0 : if ((name_result==0) && (solu_result==-1))
712 : {
713 0 : if (ov_result.size()==0)
714 0 : overall_result = true;
715 : else
716 : {
717 : bool ov_identical;
718 0 : unsigned int n = 0;
719 0 : do
720 : {
721 0 : ov_identical = (ov_result[n]==-1);
722 0 : n++;
723 : }
724 0 : while (ov_identical && n<ov_result.size());
725 0 : overall_result = ov_identical;
726 0 : }
727 : }
728 : else
729 0 : overall_result = false;
730 :
731 0 : if (verbose)
732 : {
733 0 : libMesh::out << " finished comparisons, ";
734 0 : if (overall_result)
735 0 : libMesh::out << "found no differences." << std::endl << std::endl;
736 : else
737 0 : libMesh::out << "found differences." << std::endl << std::endl;
738 : }
739 :
740 0 : return overall_result;
741 : }
742 :
743 :
744 :
745 71 : void System::update_global_solution (std::vector<Number> & global_soln) const
746 : {
747 2 : parallel_object_only();
748 :
749 71 : global_soln.resize (solution->size());
750 :
751 71 : solution->localize (global_soln);
752 71 : }
753 :
754 :
755 :
756 5578 : void System::update_global_solution (std::vector<Number> & global_soln,
757 : const processor_id_type dest_proc) const
758 : {
759 232 : parallel_object_only();
760 :
761 5578 : global_soln.resize (solution->size());
762 :
763 5578 : solution->localize_to_one (global_soln, dest_proc);
764 5578 : }
765 :
766 :
767 :
768 208446 : NumericVector<Number> & System::add_vector (std::string_view vec_name,
769 : const bool projections,
770 : const ParallelType type)
771 : {
772 6220 : parallel_object_only();
773 :
774 6220 : libmesh_assert(this->comm().verify(std::string(vec_name)));
775 6220 : libmesh_assert(this->comm().verify(int(type)));
776 6220 : libmesh_assert(this->comm().verify(projections));
777 :
778 : // Return the vector if it is already there.
779 208446 : if (auto it = this->_vectors.find(vec_name);
780 6220 : it != this->_vectors.end())
781 : {
782 : // If the projection setting has *upgraded*, change it.
783 119165 : if (projections) // only do expensive lookup if needed
784 46734 : libmesh_map_find(_vector_projections, vec_name) = projections;
785 :
786 3550 : NumericVector<Number> & vec = *it->second;
787 :
788 : // If we're in serial, our vectors are effectively SERIAL, so
789 : // we'll ignore any type setting. If we're in parallel, we
790 : // might have a type change to deal with.
791 :
792 122715 : if (this->n_processors() > 1)
793 : {
794 : // If the type setting has changed in a way we can't
795 : // perceive as an upgrade or a downgrade, scream.
796 3550 : libmesh_assert_equal_to(type == SERIAL,
797 : vec.type() == SERIAL);
798 :
799 : // If the type setting has *upgraded*, change it.
800 117339 : if (type == GHOSTED && vec.type() == PARALLEL)
801 : {
802 : // A *really* late upgrade is expensive, but better not
803 : // to risk zeroing data.
804 138 : if (vec.initialized())
805 : {
806 69 : if (!vec.closed())
807 0 : vec.close();
808 :
809 : // Ideally we'd move parallel coefficients and then
810 : // add ghosted coefficients, but copy and swap is
811 : // simpler. If anyone actually ever uses this case
812 : // for real we can look into optimizing it.
813 71 : auto new_vec = NumericVector<Number>::build(this->comm());
814 : #ifdef LIBMESH_ENABLE_GHOSTED
815 71 : new_vec->init (this->n_dofs(), this->n_local_dofs(),
816 : _dof_map->get_send_list(), /*fast=*/false,
817 4 : GHOSTED);
818 : #else
819 : libmesh_error_msg("Cannot initialize ghosted vectors when they are not enabled.");
820 : #endif
821 :
822 69 : *new_vec = vec;
823 71 : vec.swap(*new_vec);
824 65 : }
825 : else
826 : // The PARALLEL vec is not yet initialized, so we can
827 : // just "upgrade" it to GHOSTED.
828 69 : vec.set_type(type);
829 : }
830 : }
831 :
832 : // Any upgrades are done; we're happy here.
833 3550 : return vec;
834 : }
835 :
836 : // Otherwise, build the vector. The following emplace() is
837 : // guaranteed to succeed because, if we made it here, we don't
838 : // already have a vector named "vec_name". We pass the user's
839 : // requested ParallelType directly to NumericVector::build() so
840 : // that, even if the vector is not initialized now, it will get the
841 : // right type when it is initialized later.
842 : auto pr =
843 : _vectors.emplace(vec_name,
844 178562 : NumericVector<Number>::build(this->comm(),
845 : libMesh::default_solver_package(),
846 5340 : type));
847 2670 : auto buf = pr.first->second.get();
848 2670 : _vector_projections.emplace(vec_name, projections);
849 :
850 : // Vectors are primal by default
851 89281 : _vector_is_adjoint.emplace(vec_name, -1);
852 :
853 : // Initialize it if necessary
854 89281 : if (_is_initialized)
855 : {
856 45702 : if (type == GHOSTED)
857 : {
858 : #ifdef LIBMESH_ENABLE_GHOSTED
859 4200 : buf->init (this->n_dofs(), this->n_local_dofs(),
860 : _dof_map->get_send_list(), /*fast=*/false,
861 240 : GHOSTED);
862 : #else
863 : libmesh_error_msg("Cannot initialize ghosted vectors when they are not enabled.");
864 : #endif
865 : }
866 : else
867 41502 : buf->init (this->n_dofs(), this->n_local_dofs(), false, type);
868 : }
869 :
870 2670 : return *buf;
871 : }
872 :
873 39413 : void System::remove_vector (std::string_view vec_name)
874 : {
875 1126 : parallel_object_only(); // Not strictly needed, but the only safe way to keep in sync
876 :
877 39413 : if (const auto pos = _vectors.find(vec_name);
878 1126 : pos != _vectors.end())
879 : {
880 38287 : _vectors.erase(pos);
881 1126 : auto proj_it = _vector_projections.find(vec_name);
882 1126 : libmesh_assert(proj_it != _vector_projections.end());
883 38287 : _vector_projections.erase(proj_it);
884 :
885 1126 : auto adj_it = _vector_is_adjoint.find(vec_name);
886 1126 : libmesh_assert(adj_it != _vector_is_adjoint.end());
887 38287 : _vector_is_adjoint.erase(adj_it);
888 : }
889 39413 : }
890 :
891 0 : const NumericVector<Number> * System::request_vector (std::string_view vec_name) const
892 : {
893 0 : if (const auto pos = _vectors.find(vec_name);
894 0 : pos != _vectors.end())
895 0 : return pos->second.get();
896 :
897 : // Otherwise, vec_name was not found
898 0 : return nullptr;
899 : }
900 :
901 :
902 :
903 0 : NumericVector<Number> * System::request_vector (std::string_view vec_name)
904 : {
905 0 : if (auto pos = _vectors.find(vec_name);
906 0 : pos != _vectors.end())
907 0 : return pos->second.get();
908 :
909 : // Otherwise, vec_name was not found
910 0 : return nullptr;
911 : }
912 :
913 :
914 :
915 0 : const NumericVector<Number> * System::request_vector (const unsigned int vec_num) const
916 : {
917 : // If we don't have that many vectors, return nullptr
918 0 : if (vec_num >= _vectors.size())
919 0 : return nullptr;
920 :
921 : // Otherwise return a pointer to the vec_num'th vector
922 0 : auto it = vectors_begin();
923 0 : std::advance(it, vec_num);
924 0 : return it->second.get();
925 : }
926 :
927 :
928 :
929 0 : NumericVector<Number> * System::request_vector (const unsigned int vec_num)
930 : {
931 : // If we don't have that many vectors, return nullptr
932 0 : if (vec_num >= _vectors.size())
933 0 : return nullptr;
934 :
935 : // Otherwise return a pointer to the vec_num'th vector
936 0 : auto it = vectors_begin();
937 0 : std::advance(it, vec_num);
938 0 : return it->second.get();
939 : }
940 :
941 :
942 :
943 5792 : const NumericVector<Number> & System::get_vector (std::string_view vec_name) const
944 : {
945 5792 : return *(libmesh_map_find(_vectors, vec_name));
946 : }
947 :
948 :
949 :
950 6382860 : NumericVector<Number> & System::get_vector (std::string_view vec_name)
951 : {
952 6382860 : return *(libmesh_map_find(_vectors, vec_name));
953 : }
954 :
955 :
956 :
957 0 : const NumericVector<Number> & System::get_vector (const unsigned int vec_num) const
958 : {
959 : // If we don't have that many vectors, throw an error
960 0 : libmesh_assert_less(vec_num, _vectors.size());
961 :
962 : // Otherwise return a reference to the vec_num'th vector
963 0 : auto it = vectors_begin();
964 0 : std::advance(it, vec_num);
965 0 : return *(it->second);
966 : }
967 :
968 :
969 :
970 0 : NumericVector<Number> & System::get_vector (const unsigned int vec_num)
971 : {
972 : // If we don't have that many vectors, throw an error
973 0 : libmesh_assert_less(vec_num, _vectors.size());
974 :
975 : // Otherwise return a reference to the vec_num'th vector
976 0 : auto it = vectors_begin();
977 0 : std::advance(it, vec_num);
978 0 : return *(it->second);
979 : }
980 :
981 :
982 :
983 0 : const std::string & System::vector_name (const unsigned int vec_num) const
984 : {
985 : // If we don't have that many vectors, throw an error
986 0 : libmesh_assert_less(vec_num, _vectors.size());
987 :
988 : // Otherwise return a reference to the vec_num'th vector name
989 0 : auto it = vectors_begin();
990 0 : std::advance(it, vec_num);
991 0 : return it->first;
992 : }
993 :
994 66238 : const std::string & System::vector_name (const NumericVector<Number> & vec_reference) const
995 : {
996 : // Linear search for a vector whose pointer matches vec_reference
997 62454 : auto it = std::find_if(vectors_begin(), vectors_end(),
998 31852 : [&vec_reference](const decltype(_vectors)::value_type & pr)
999 17818 : { return &vec_reference == pr.second.get(); });
1000 :
1001 : // Before returning, make sure we didn't loop till the end and not find any match
1002 1892 : libmesh_assert (it != vectors_end());
1003 :
1004 : // Return the string associated with the current vector
1005 66238 : return it->first;
1006 : }
1007 :
1008 :
1009 :
1010 18659 : SparseMatrix<Number> & System::add_matrix (std::string_view mat_name,
1011 : const ParallelType type,
1012 : const MatrixBuildType mat_build_type)
1013 : {
1014 540 : parallel_object_only();
1015 :
1016 540 : libmesh_assert(this->comm().verify(std::string(mat_name)));
1017 540 : libmesh_assert(this->comm().verify(int(type)));
1018 540 : libmesh_assert(this->comm().verify(int(mat_build_type)));
1019 :
1020 : // Return the matrix if it is already there.
1021 18659 : if (auto it = this->_matrices.find(mat_name);
1022 540 : it != this->_matrices.end())
1023 0 : return *it->second;
1024 :
1025 : // Otherwise build the matrix to return.
1026 18119 : std::unique_ptr<SparseMatrix<Number>> matrix;
1027 18659 : if (this->has_static_condensation())
1028 : {
1029 0 : if (mat_build_type == MatrixBuildType::DIAGONAL)
1030 0 : libmesh_error_msg(
1031 : "We do not currently support static condensation of the diagonal matrix type");
1032 0 : matrix = std::make_unique<StaticCondensation>(this->get_mesh(),
1033 : *this,
1034 : this->get_dof_map(),
1035 0 : this->get_dof_map().get_static_condensation());
1036 : }
1037 : else
1038 36778 : matrix = SparseMatrix<Number>::build(this->comm(), libMesh::default_solver_package());
1039 540 : auto & mat = *matrix;
1040 :
1041 540 : _matrices.emplace(mat_name, std::move(matrix));
1042 :
1043 540 : _matrix_types.emplace(mat_name, type);
1044 :
1045 : // Initialize it first if we've already initialized the others.
1046 18659 : this->late_matrix_init(mat, type);
1047 :
1048 540 : return mat;
1049 17579 : }
1050 :
1051 :
1052 :
1053 350 : SparseMatrix<Number> & System::add_matrix (std::string_view mat_name,
1054 : std::unique_ptr<SparseMatrix<Number>> matrix,
1055 : const ParallelType type)
1056 : {
1057 10 : parallel_object_only();
1058 :
1059 10 : libmesh_assert(this->comm().verify(std::string(mat_name)));
1060 10 : libmesh_assert(this->comm().verify(int(type)));
1061 :
1062 10 : auto [it, inserted] = _matrices.emplace(mat_name, std::move(matrix));
1063 350 : libmesh_error_msg_if(!inserted,
1064 : "Tried to add '" << mat_name << "' but the matrix already exists");
1065 :
1066 10 : _matrix_types.emplace(mat_name, type);
1067 :
1068 10 : SparseMatrix<Number> & mat = *(it->second);
1069 :
1070 : // Initialize it first if we've already initialized the others.
1071 350 : this->late_matrix_init(mat, type);
1072 :
1073 350 : return mat;
1074 : }
1075 :
1076 19009 : void System::late_matrix_init(SparseMatrix<Number> & mat,
1077 : ParallelType type)
1078 : {
1079 19009 : if (_matrices_initialized)
1080 : {
1081 0 : this->get_dof_map().attach_matrix(mat);
1082 0 : mat.init(type);
1083 : }
1084 19009 : }
1085 :
1086 :
1087 :
1088 :
1089 0 : void System::remove_matrix (std::string_view mat_name)
1090 : {
1091 0 : parallel_object_only(); // Not strictly needed, but the only safe way to keep in sync
1092 :
1093 0 : if (const auto pos = _matrices.find(mat_name);
1094 0 : pos != _matrices.end())
1095 0 : _matrices.erase(pos); // erase()'d entries are destroyed
1096 0 : }
1097 :
1098 :
1099 :
1100 14 : const SparseMatrix<Number> * System::request_matrix (std::string_view mat_name) const
1101 : {
1102 14 : if (const auto pos = _matrices.find(mat_name);
1103 14 : pos != _matrices.end())
1104 14 : return pos->second.get();
1105 :
1106 : // Otherwise, mat_name does not exist
1107 0 : return nullptr;
1108 : }
1109 :
1110 :
1111 :
1112 244176 : SparseMatrix<Number> * System::request_matrix (std::string_view mat_name)
1113 : {
1114 244176 : if (auto pos = _matrices.find(mat_name);
1115 7160 : pos != _matrices.end())
1116 2 : return pos->second.get();
1117 :
1118 : // Otherwise, mat_name does not exist
1119 7158 : return nullptr;
1120 : }
1121 :
1122 :
1123 :
1124 0 : const SparseMatrix<Number> & System::get_matrix (std::string_view mat_name) const
1125 : {
1126 0 : return *libmesh_map_find(_matrices, mat_name);
1127 : }
1128 :
1129 :
1130 :
1131 659231 : SparseMatrix<Number> & System::get_matrix (std::string_view mat_name)
1132 : {
1133 659231 : return *libmesh_map_find(_matrices, mat_name);
1134 : }
1135 :
1136 :
1137 :
1138 69038 : void System::set_vector_preservation (const std::string & vec_name,
1139 : bool preserve)
1140 : {
1141 1972 : parallel_object_only(); // Not strictly needed, but the only safe way to keep in sync
1142 :
1143 69038 : _vector_projections[vec_name] = preserve;
1144 69038 : }
1145 :
1146 :
1147 :
1148 159257 : bool System::vector_preservation (std::string_view vec_name) const
1149 : {
1150 159257 : if (auto it = _vector_projections.find(vec_name);
1151 4550 : it != _vector_projections.end())
1152 159257 : return it->second;
1153 :
1154 : // vec_name was not in the map, return false
1155 0 : return false;
1156 : }
1157 :
1158 :
1159 :
1160 24746 : void System::set_vector_as_adjoint (const std::string & vec_name,
1161 : int qoi_num)
1162 : {
1163 798 : parallel_object_only(); // Not strictly needed, but the only safe way to keep in sync
1164 :
1165 : // We reserve -1 for vectors which get primal constraints, -2 for
1166 : // vectors which get no constraints
1167 798 : libmesh_assert_greater_equal(qoi_num, -2);
1168 24746 : _vector_is_adjoint[vec_name] = qoi_num;
1169 24746 : }
1170 :
1171 :
1172 :
1173 20375 : int System::vector_is_adjoint (std::string_view vec_name) const
1174 : {
1175 820 : const auto it = _vector_is_adjoint.find(vec_name);
1176 820 : libmesh_assert(it != _vector_is_adjoint.end());
1177 20375 : return it->second;
1178 : }
1179 :
1180 :
1181 :
1182 142 : NumericVector<Number> & System::add_sensitivity_solution (unsigned int i)
1183 : {
1184 142 : std::ostringstream sensitivity_name;
1185 138 : sensitivity_name << "sensitivity_solution" << i;
1186 :
1187 284 : return this->add_vector(sensitivity_name.str());
1188 134 : }
1189 :
1190 :
1191 :
1192 284 : NumericVector<Number> & System::get_sensitivity_solution (unsigned int i)
1193 : {
1194 284 : std::ostringstream sensitivity_name;
1195 276 : sensitivity_name << "sensitivity_solution" << i;
1196 :
1197 568 : return this->get_vector(sensitivity_name.str());
1198 268 : }
1199 :
1200 :
1201 :
1202 0 : const NumericVector<Number> & System::get_sensitivity_solution (unsigned int i) const
1203 : {
1204 0 : std::ostringstream sensitivity_name;
1205 0 : sensitivity_name << "sensitivity_solution" << i;
1206 :
1207 0 : return this->get_vector(sensitivity_name.str());
1208 0 : }
1209 :
1210 :
1211 :
1212 0 : NumericVector<Number> & System::add_weighted_sensitivity_solution ()
1213 : {
1214 0 : return this->add_vector("weighted_sensitivity_solution");
1215 : }
1216 :
1217 :
1218 :
1219 0 : NumericVector<Number> & System::get_weighted_sensitivity_solution ()
1220 : {
1221 0 : return this->get_vector("weighted_sensitivity_solution");
1222 : }
1223 :
1224 :
1225 :
1226 0 : const NumericVector<Number> & System::get_weighted_sensitivity_solution () const
1227 : {
1228 0 : return this->get_vector("weighted_sensitivity_solution");
1229 : }
1230 :
1231 :
1232 :
1233 24746 : NumericVector<Number> & System::add_adjoint_solution (unsigned int i)
1234 : {
1235 24746 : std::ostringstream adjoint_name;
1236 23948 : adjoint_name << "adjoint_solution" << i;
1237 :
1238 24746 : NumericVector<Number> & returnval = this->add_vector(adjoint_name.str());
1239 25544 : this->set_vector_as_adjoint(adjoint_name.str(), i);
1240 25544 : return returnval;
1241 23150 : }
1242 :
1243 :
1244 :
1245 2885957 : NumericVector<Number> & System::get_adjoint_solution (unsigned int i)
1246 : {
1247 2885957 : std::ostringstream adjoint_name;
1248 2663375 : adjoint_name << "adjoint_solution" << i;
1249 :
1250 5771914 : return this->get_vector(adjoint_name.str());
1251 2440605 : }
1252 :
1253 :
1254 :
1255 5792 : const NumericVector<Number> & System::get_adjoint_solution (unsigned int i) const
1256 : {
1257 5792 : std::ostringstream adjoint_name;
1258 5480 : adjoint_name << "adjoint_solution" << i;
1259 :
1260 11584 : return this->get_vector(adjoint_name.str());
1261 5168 : }
1262 :
1263 :
1264 :
1265 0 : NumericVector<Number> & System::add_weighted_sensitivity_adjoint_solution (unsigned int i)
1266 : {
1267 0 : std::ostringstream adjoint_name;
1268 0 : adjoint_name << "weighted_sensitivity_adjoint_solution" << i;
1269 :
1270 0 : NumericVector<Number> & returnval = this->add_vector(adjoint_name.str());
1271 0 : this->set_vector_as_adjoint(adjoint_name.str(), i);
1272 0 : return returnval;
1273 0 : }
1274 :
1275 :
1276 :
1277 0 : NumericVector<Number> & System::get_weighted_sensitivity_adjoint_solution (unsigned int i)
1278 : {
1279 0 : std::ostringstream adjoint_name;
1280 0 : adjoint_name << "weighted_sensitivity_adjoint_solution" << i;
1281 :
1282 0 : return this->get_vector(adjoint_name.str());
1283 0 : }
1284 :
1285 :
1286 :
1287 0 : const NumericVector<Number> & System::get_weighted_sensitivity_adjoint_solution (unsigned int i) const
1288 : {
1289 0 : std::ostringstream adjoint_name;
1290 0 : adjoint_name << "weighted_sensitivity_adjoint_solution" << i;
1291 :
1292 0 : return this->get_vector(adjoint_name.str());
1293 0 : }
1294 :
1295 :
1296 :
1297 24817 : NumericVector<Number> & System::add_adjoint_rhs (unsigned int i)
1298 : {
1299 24817 : std::ostringstream adjoint_rhs_name;
1300 24017 : adjoint_rhs_name << "adjoint_rhs" << i;
1301 :
1302 49634 : return this->add_vector(adjoint_rhs_name.str(), false);
1303 23217 : }
1304 :
1305 :
1306 :
1307 2563069 : NumericVector<Number> & System::get_adjoint_rhs (unsigned int i)
1308 : {
1309 2563069 : std::ostringstream adjoint_rhs_name;
1310 2332935 : adjoint_rhs_name << "adjoint_rhs" << i;
1311 :
1312 5126138 : return this->get_vector(adjoint_rhs_name.str());
1313 2102801 : }
1314 :
1315 :
1316 :
1317 0 : const NumericVector<Number> & System::get_adjoint_rhs (unsigned int i) const
1318 : {
1319 0 : std::ostringstream adjoint_rhs_name;
1320 0 : adjoint_rhs_name << "adjoint_rhs" << i;
1321 :
1322 0 : return this->get_vector(adjoint_rhs_name.str());
1323 0 : }
1324 :
1325 :
1326 :
1327 7430 : NumericVector<Number> & System::add_sensitivity_rhs (unsigned int i)
1328 : {
1329 7430 : std::ostringstream sensitivity_rhs_name;
1330 7218 : sensitivity_rhs_name << "sensitivity_rhs" << i;
1331 :
1332 14860 : return this->add_vector(sensitivity_rhs_name.str(), false);
1333 7006 : }
1334 :
1335 :
1336 :
1337 7430 : NumericVector<Number> & System::get_sensitivity_rhs (unsigned int i)
1338 : {
1339 7430 : std::ostringstream sensitivity_rhs_name;
1340 7218 : sensitivity_rhs_name << "sensitivity_rhs" << i;
1341 :
1342 14860 : return this->get_vector(sensitivity_rhs_name.str());
1343 7006 : }
1344 :
1345 :
1346 :
1347 0 : const NumericVector<Number> & System::get_sensitivity_rhs (unsigned int i) const
1348 : {
1349 0 : std::ostringstream sensitivity_rhs_name;
1350 0 : sensitivity_rhs_name << "sensitivity_rhs" << i;
1351 :
1352 0 : return this->get_vector(sensitivity_rhs_name.str());
1353 0 : }
1354 :
1355 :
1356 :
1357 265515 : unsigned int System::add_variable (std::string_view var,
1358 : const FEType & type,
1359 : const std::set<subdomain_id_type> * const active_subdomains)
1360 : {
1361 7562 : parallel_object_only(); // Not strictly needed, but the only safe way to keep in sync
1362 :
1363 7562 : libmesh_assert(this->comm().verify(std::string(var)));
1364 7562 : libmesh_assert(this->comm().verify(type));
1365 7562 : libmesh_assert(this->comm().verify((active_subdomains == nullptr)));
1366 :
1367 7562 : if (active_subdomains)
1368 148 : libmesh_assert(this->comm().verify(active_subdomains->size()));
1369 :
1370 : // Make sure the variable isn't there already
1371 : // or if it is, that it's the type we want
1372 337794 : for (auto v : make_range(this->n_vars()))
1373 72279 : if (this->variable_name(v) == var)
1374 : {
1375 0 : if (this->variable_type(v) == type)
1376 : {
1377 : // Check whether the existing variable's active subdomains also matches
1378 : // the incoming variable's active subdomains. If they don't match, then
1379 : // either it is an error by the user or the user is trying to change the
1380 : // subdomain restriction after the variable has already been added, which
1381 : // is not supported.
1382 0 : const Variable & existing_var = this->variable(v);
1383 :
1384 : // Check whether active_subdomains is not provided/empty and the existing_var is implicitly_active()
1385 0 : bool check1 = (!active_subdomains || active_subdomains->empty()) && existing_var.implicitly_active();
1386 :
1387 : // Check if the provided active_subdomains is equal to the existing_var's active_subdomains
1388 0 : bool check2 = (active_subdomains && (*active_subdomains == existing_var.active_subdomains()));
1389 :
1390 : // If either of these checks passed, then we already have this variable
1391 0 : if (check1 || check2)
1392 0 : return _variables[v].number();
1393 : }
1394 :
1395 0 : libmesh_error_msg("ERROR: incompatible variable " << var << " has already been added for this system!");
1396 : }
1397 :
1398 7562 : libmesh_assert(!this->is_initialized());
1399 :
1400 265515 : if (this->n_variable_groups())
1401 : {
1402 : // Optimize for VariableGroups here - if the user is adding multiple
1403 : // variables of the same FEType and subdomain restriction, catch
1404 : // that here and add them as members of the same VariableGroup.
1405 : //
1406 : // start by setting this flag to whatever the user has requested
1407 : // and then consider the conditions which should negate it.
1408 1640 : bool should_be_in_vg = this->identify_variable_groups();
1409 :
1410 820 : VariableGroup & vg(_variable_groups.back());
1411 :
1412 : // get a pointer to their subdomain restriction, if any.
1413 : const std::set<subdomain_id_type> * const
1414 29287 : their_active_subdomains (vg.implicitly_active() ?
1415 104 : nullptr : &vg.active_subdomains());
1416 :
1417 : // Different types?
1418 1151 : if (vg.type() != type)
1419 280 : should_be_in_vg = false;
1420 :
1421 : // they are restricted, we aren't?
1422 29391 : if (their_active_subdomains &&
1423 3674 : (!active_subdomains || (active_subdomains && active_subdomains->empty())))
1424 0 : should_be_in_vg = false;
1425 :
1426 : // they aren't restricted, we are?
1427 29287 : if (!their_active_subdomains && (active_subdomains && !active_subdomains->empty()))
1428 2 : should_be_in_vg = false;
1429 :
1430 29287 : if (their_active_subdomains && active_subdomains)
1431 : // restricted to different sets?
1432 3674 : if (*their_active_subdomains != *active_subdomains)
1433 52 : should_be_in_vg = false;
1434 :
1435 : // OK, after all that, append the variable to the vg if none of the conditions
1436 : // were violated
1437 27499 : if (should_be_in_vg)
1438 : {
1439 506 : const unsigned int curr_n_vars = this->n_vars();
1440 :
1441 17678 : std::string varstr(var);
1442 :
1443 18184 : _variable_numbers[varstr] = curr_n_vars;
1444 18184 : vg.append (std::move(varstr));
1445 36368 : _variables.push_back(vg(vg.n_variables()-1));
1446 :
1447 506 : return curr_n_vars;
1448 : }
1449 : }
1450 :
1451 : // otherwise, fall back to adding a single variable group
1452 487606 : return this->add_variables (std::vector<std::string>(1, std::string(var)),
1453 : type,
1454 7056 : active_subdomains);
1455 : }
1456 :
1457 :
1458 :
1459 260395 : unsigned int System::add_variable (std::string_view var,
1460 : const Order order,
1461 : const FEFamily family,
1462 : const std::set<subdomain_id_type> * const active_subdomains)
1463 : {
1464 260395 : return this->add_variable(var,
1465 267807 : FEType(order, family),
1466 267807 : active_subdomains);
1467 : }
1468 :
1469 :
1470 :
1471 247402 : unsigned int System::add_variables (const std::vector<std::string> & vars,
1472 : const FEType & type,
1473 : const std::set<subdomain_id_type> * const active_subdomains)
1474 : {
1475 7058 : parallel_object_only(); // Not strictly needed, but the only safe way to keep in sync
1476 :
1477 7058 : libmesh_assert(!this->is_initialized());
1478 :
1479 7058 : libmesh_assert(this->comm().verify(vars.size()));
1480 7058 : libmesh_assert(this->comm().verify(type));
1481 7058 : libmesh_assert(this->comm().verify((active_subdomains == nullptr)));
1482 :
1483 7058 : if (active_subdomains)
1484 96 : libmesh_assert(this->comm().verify(active_subdomains->size()));
1485 :
1486 : // Make sure the variable isn't there already
1487 : // or if it is, that it's the type we want
1488 7594733 : for (auto ovar : vars)
1489 : {
1490 207056 : libmesh_assert(this->comm().verify(ovar));
1491 :
1492 7369683 : for (auto v : make_range(this->n_vars()))
1493 21718 : if (this->variable_name(v) == ovar)
1494 : {
1495 0 : if (this->variable_type(v) == type)
1496 0 : return _variables[v].number();
1497 :
1498 0 : libmesh_error_msg("ERROR: incompatible variable " << ovar << " has already been added for this system!");
1499 : }
1500 : }
1501 :
1502 247402 : if (this->n_variable_groups())
1503 : {
1504 : // Optimize for VariableGroups here - if the user is adding multiple
1505 : // variables of the same FEType and subdomain restriction, catch
1506 : // that here and add them as members of the same VariableGroup.
1507 : //
1508 : // start by setting this flag to whatever the user has requested
1509 : // and then consider the conditions which should negate it.
1510 628 : bool should_be_in_vg = this->identify_variable_groups();
1511 :
1512 314 : VariableGroup & vg(_variable_groups.back());
1513 :
1514 : // get a pointer to their subdomain restriction, if any.
1515 : const std::set<subdomain_id_type> * const
1516 11103 : their_active_subdomains (vg.implicitly_active() ?
1517 52 : nullptr : &vg.active_subdomains());
1518 :
1519 : // Different types?
1520 427 : if (vg.type() != type)
1521 280 : should_be_in_vg = false;
1522 :
1523 : // they are restricted, we aren't?
1524 11155 : if (their_active_subdomains &&
1525 1840 : (!active_subdomains || (active_subdomains && active_subdomains->empty())))
1526 0 : should_be_in_vg = false;
1527 :
1528 : // they aren't restricted, we are?
1529 11103 : if (!their_active_subdomains && (active_subdomains && !active_subdomains->empty()))
1530 2 : should_be_in_vg = false;
1531 :
1532 11103 : if (their_active_subdomains && active_subdomains)
1533 : // restricted to different sets?
1534 1840 : if (*their_active_subdomains != *active_subdomains)
1535 52 : should_be_in_vg = false;
1536 :
1537 : // If after all that none of the conditions were violated,
1538 : // append the variables to the vg and we're done
1539 9315 : if (should_be_in_vg)
1540 : {
1541 0 : unsigned int curr_n_vars = this->n_vars();
1542 :
1543 0 : for (auto ovar : vars)
1544 : {
1545 0 : curr_n_vars = this->n_vars();
1546 :
1547 0 : vg.append (ovar);
1548 :
1549 0 : _variables.push_back(vg(vg.n_variables()-1));
1550 0 : _variable_numbers[ovar] = curr_n_vars;
1551 : }
1552 0 : return curr_n_vars;
1553 : }
1554 : }
1555 :
1556 7058 : const unsigned int curr_n_vars = this->n_vars();
1557 :
1558 247402 : const unsigned int next_first_component = this->n_components();
1559 :
1560 : // We weren't able to add to an existing variable group, so
1561 : // add a new variable group to the list
1562 494708 : _variable_groups.push_back((active_subdomains == nullptr) ?
1563 : VariableGroup(this, vars, curr_n_vars,
1564 : next_first_component, type) :
1565 : VariableGroup(this, vars, curr_n_vars,
1566 : next_first_component, type, *active_subdomains));
1567 :
1568 7058 : const VariableGroup & vg (_variable_groups.back());
1569 :
1570 : // Add each component of the group individually
1571 7594733 : for (auto v : make_range(vars.size()))
1572 : {
1573 14487606 : _variables.push_back (vg(v));
1574 7554387 : _variable_numbers[vars[v]] = curr_n_vars+v;
1575 : }
1576 :
1577 7058 : libmesh_assert_equal_to ((curr_n_vars+vars.size()), this->n_vars());
1578 :
1579 : // BSK - Defer this now to System::init_data() so we can detect
1580 : // VariableGroups 12/28/2012
1581 : // // Add the variable group to the _dof_map
1582 : // _dof_map->add_variable_group (vg);
1583 :
1584 : // Return the number of the new variable
1585 261518 : return cast_int<unsigned int>(curr_n_vars+vars.size()-1);
1586 : }
1587 :
1588 :
1589 :
1590 71 : unsigned int System::add_variables (const std::vector<std::string> & vars,
1591 : const Order order,
1592 : const FEFamily family,
1593 : const std::set<subdomain_id_type> * const active_subdomains)
1594 : {
1595 71 : return this->add_variables(vars,
1596 73 : FEType(order, family),
1597 73 : active_subdomains);
1598 : }
1599 :
1600 :
1601 :
1602 960 : bool System::has_variable (std::string_view var) const
1603 : {
1604 960 : return _variable_numbers.count(var);
1605 : }
1606 :
1607 :
1608 :
1609 8865541 : unsigned int System::variable_number (std::string_view var) const
1610 : {
1611 8865541 : auto var_num = libmesh_map_find(_variable_numbers, var);
1612 385695 : libmesh_assert_equal_to (_variables[var_num].name(), var);
1613 8865541 : return var_num;
1614 : }
1615 :
1616 :
1617 493 : void System::get_all_variable_numbers(std::vector<unsigned int> & all_variable_numbers) const
1618 : {
1619 493 : all_variable_numbers.resize(n_vars());
1620 :
1621 14 : unsigned int count = 0;
1622 1266 : for (auto vn : _variable_numbers)
1623 795 : all_variable_numbers[count++] = vn.second;
1624 493 : }
1625 :
1626 :
1627 142 : void System::local_dof_indices(const unsigned int var,
1628 : std::set<dof_id_type> & var_indices) const
1629 : {
1630 : // Make sure the set is clear
1631 4 : var_indices.clear();
1632 :
1633 8 : std::vector<dof_id_type> dof_indices;
1634 :
1635 : const dof_id_type
1636 4 : first_local = this->get_dof_map().first_dof(),
1637 4 : end_local = this->get_dof_map().end_dof();
1638 :
1639 : // Begin the loop over the elements
1640 984 : for (const auto & elem : this->get_mesh().active_local_element_ptr_range())
1641 : {
1642 384 : this->get_dof_map().dof_indices (elem, dof_indices, var);
1643 :
1644 1152 : for (dof_id_type dof : dof_indices)
1645 : //If the dof is owned by the local processor
1646 768 : if (first_local <= dof && dof < end_local)
1647 503 : var_indices.insert(dof);
1648 134 : }
1649 :
1650 : // we may have missed assigning DOFs to nodes that we own
1651 : // but to which we have no connected elements matching our
1652 : // variable restriction criterion. this will happen, for example,
1653 : // if variable V is restricted to subdomain S. We may not own
1654 : // any elements which live in S, but we may own nodes which are
1655 : // *connected* to elements which do.
1656 1380 : for (const auto & node : this->get_mesh().local_node_ptr_range())
1657 : {
1658 50 : libmesh_assert(node);
1659 600 : this->get_dof_map().dof_indices (node, dof_indices, var);
1660 960 : for (auto dof : dof_indices)
1661 360 : if (first_local <= dof && dof < end_local)
1662 330 : var_indices.insert(dof);
1663 134 : }
1664 142 : }
1665 :
1666 :
1667 :
1668 0 : void System::zero_variable (NumericVector<Number> & v,
1669 : unsigned int var_num) const
1670 : {
1671 : /* Make sure the call makes sense. */
1672 0 : libmesh_assert_less (var_num, this->n_vars());
1673 :
1674 : /* Get a reference to the mesh. */
1675 0 : const MeshBase & mesh = this->get_mesh();
1676 :
1677 : /* Check which system we are. */
1678 0 : const unsigned int sys_num = this->number();
1679 :
1680 : // Loop over nodes.
1681 0 : for (const auto & node : mesh.local_node_ptr_range())
1682 : {
1683 0 : unsigned int n_comp = node->n_comp(sys_num,var_num);
1684 0 : for (unsigned int i=0; i<n_comp; i++)
1685 : {
1686 0 : const dof_id_type index = node->dof_number(sys_num,var_num,i);
1687 0 : v.set(index,0.0);
1688 : }
1689 0 : }
1690 :
1691 : // Loop over elements.
1692 0 : for (const auto & elem : mesh.active_local_element_ptr_range())
1693 : {
1694 0 : unsigned int n_comp = elem->n_comp(sys_num,var_num);
1695 0 : for (unsigned int i=0; i<n_comp; i++)
1696 : {
1697 0 : const dof_id_type index = elem->dof_number(sys_num,var_num,i);
1698 0 : v.set(index,0.0);
1699 : }
1700 0 : }
1701 0 : }
1702 :
1703 :
1704 :
1705 0 : Real System::discrete_var_norm(const NumericVector<Number> & v,
1706 : unsigned int var,
1707 : FEMNormType norm_type) const
1708 : {
1709 0 : std::set<dof_id_type> var_indices;
1710 0 : local_dof_indices(var, var_indices);
1711 :
1712 0 : if (norm_type == DISCRETE_L1)
1713 0 : return v.subset_l1_norm(var_indices);
1714 0 : if (norm_type == DISCRETE_L2)
1715 0 : return v.subset_l2_norm(var_indices);
1716 0 : if (norm_type == DISCRETE_L_INF)
1717 0 : return v.subset_linfty_norm(var_indices);
1718 : else
1719 0 : libmesh_error_msg("Invalid norm_type = " << Utility::enum_to_string(norm_type));
1720 : }
1721 :
1722 :
1723 :
1724 106884 : Real System::calculate_norm(const NumericVector<Number> & v,
1725 : unsigned int var,
1726 : FEMNormType norm_type,
1727 : std::set<unsigned int> * skip_dimensions) const
1728 : {
1729 : //short circuit to save time
1730 106884 : if (norm_type == DISCRETE_L1 ||
1731 106884 : norm_type == DISCRETE_L2 ||
1732 : norm_type == DISCRETE_L_INF)
1733 0 : return discrete_var_norm(v,var,norm_type);
1734 :
1735 : // Not a discrete norm
1736 113164 : std::vector<FEMNormType> norms(this->n_vars(), L2);
1737 109930 : std::vector<Real> weights(this->n_vars(), 0.0);
1738 106884 : norms[var] = norm_type;
1739 106884 : weights[var] = 1.0;
1740 207488 : Real val = this->calculate_norm(v, SystemNorm(norms, weights), skip_dimensions);
1741 3234 : return val;
1742 : }
1743 :
1744 :
1745 :
1746 125981 : Real System::calculate_norm(const NumericVector<Number> & v,
1747 : const SystemNorm & norm,
1748 : std::set<unsigned int> * skip_dimensions) const
1749 : {
1750 : // This function must be run on all processors at once
1751 3968 : parallel_object_only();
1752 :
1753 7936 : LOG_SCOPE ("calculate_norm()", "System");
1754 :
1755 : // Zero the norm before summation
1756 125981 : Real v_norm = 0.;
1757 :
1758 125981 : if (norm.is_discrete())
1759 : {
1760 : //Check to see if all weights are 1.0 and all types are equal
1761 11340 : FEMNormType norm_type0 = norm.type(0);
1762 324 : unsigned int check_var = 0, check_end = this->n_vars();
1763 22680 : for (; check_var != check_end; ++check_var)
1764 11340 : if ((norm.weight(check_var) != 1.0) || (norm.type(check_var) != norm_type0))
1765 0 : break;
1766 :
1767 : //All weights were 1.0 so just do the full vector discrete norm
1768 11340 : if (check_var == this->n_vars())
1769 : {
1770 11340 : if (norm_type0 == DISCRETE_L1)
1771 0 : return v.l1_norm();
1772 11340 : if (norm_type0 == DISCRETE_L2)
1773 11340 : return v.l2_norm();
1774 0 : if (norm_type0 == DISCRETE_L_INF)
1775 0 : return v.linfty_norm();
1776 : else
1777 0 : libmesh_error_msg("Invalid norm_type0 = " << Utility::enum_to_string(norm_type0));
1778 : }
1779 :
1780 0 : for (auto var : make_range(this->n_vars()))
1781 : {
1782 : // Skip any variables we don't need to integrate
1783 0 : if (norm.weight(var) == 0.0)
1784 0 : continue;
1785 :
1786 0 : v_norm += norm.weight(var) * discrete_var_norm(v, var, norm.type(var));
1787 : }
1788 :
1789 0 : return v_norm;
1790 : }
1791 :
1792 : // Localize the potentially parallel vector
1793 114641 : std::unique_ptr<NumericVector<Number>> local_v = NumericVector<Number>::build(this->comm());
1794 118285 : local_v->init(v.size(), v.local_size(), _dof_map->get_send_list(),
1795 7100 : true, GHOSTED);
1796 118097 : v.localize (*local_v, _dof_map->get_send_list());
1797 :
1798 : // I'm not sure how best to mix Hilbert norms on some variables (for
1799 : // which we'll want to square then sum then square root) with norms
1800 : // like L_inf (for which we'll just want to take an absolute value
1801 : // and then sum).
1802 3644 : bool using_hilbert_norm = true,
1803 3644 : using_nonhilbert_norm = true;
1804 :
1805 : // Loop over all variables
1806 229282 : for (auto var : make_range(this->n_vars()))
1807 : {
1808 : // Skip any variables we don't need to integrate
1809 114641 : Real norm_weight_sq = norm.weight_sq(var);
1810 114641 : if (norm_weight_sq == 0.0)
1811 0 : continue;
1812 114641 : Real norm_weight = norm.weight(var);
1813 :
1814 : // Check for unimplemented norms (rather than just returning 0).
1815 114641 : FEMNormType norm_type = norm.type(var);
1816 114641 : if ((norm_type==H1) ||
1817 111041 : (norm_type==H2) ||
1818 111021 : (norm_type==L2) ||
1819 111013 : (norm_type==H1_SEMINORM) ||
1820 : (norm_type==H2_SEMINORM))
1821 : {
1822 114041 : if (!using_hilbert_norm)
1823 0 : libmesh_not_implemented();
1824 3628 : using_nonhilbert_norm = false;
1825 : }
1826 616 : else if ((norm_type==L1) ||
1827 592 : (norm_type==L_INF) ||
1828 584 : (norm_type==W1_INF_SEMINORM) ||
1829 : (norm_type==W2_INF_SEMINORM))
1830 : {
1831 600 : if (!using_nonhilbert_norm)
1832 0 : libmesh_not_implemented();
1833 16 : using_hilbert_norm = false;
1834 : }
1835 : else
1836 0 : libmesh_not_implemented();
1837 :
1838 3644 : const FEType & fe_type = this->get_dof_map().variable_type(var);
1839 :
1840 : // Allow space for dims 0-3, and for both scalar and vector
1841 : // elements, even if we don't use them all
1842 121741 : std::vector<std::unique_ptr<FEBase>> fe_ptrs(4);
1843 121741 : std::vector<std::unique_ptr<FEVectorBase>> vec_fe_ptrs(4);
1844 121741 : std::vector<std::unique_ptr<QBase>> q_rules(4);
1845 :
1846 114641 : const std::set<unsigned char> & elem_dims = _mesh.elem_dimensions();
1847 :
1848 : // Prepare finite elements for each dimension present in the mesh
1849 230692 : for (const auto & dim : elem_dims)
1850 : {
1851 116127 : if (skip_dimensions && skip_dimensions->find(dim) != skip_dimensions->end())
1852 1372 : continue;
1853 :
1854 : // Construct quadrature and finite element objects
1855 114641 : q_rules[dim] = fe_type.default_quadrature_rule (dim);
1856 :
1857 114641 : const FEFieldType field_type = FEInterface::field_type(fe_type);
1858 114641 : if (field_type == TYPE_SCALAR)
1859 : {
1860 114215 : fe_ptrs[dim] = FEBase::build(dim, fe_type);
1861 121103 : fe_ptrs[dim]->attach_quadrature_rule (q_rules[dim].get());
1862 : }
1863 : else
1864 : {
1865 426 : vec_fe_ptrs[dim] = FEVectorBase::build(dim, fe_type);
1866 450 : vec_fe_ptrs[dim]->attach_quadrature_rule (q_rules[dim].get());
1867 12 : libmesh_assert_equal_to(field_type, TYPE_VECTOR);
1868 : }
1869 :
1870 : }
1871 :
1872 7288 : std::vector<dof_id_type> dof_indices;
1873 :
1874 : // Begin the loop over the elements
1875 31243676 : for (const auto & elem : this->get_mesh().active_local_element_ptr_range())
1876 : {
1877 17135144 : const unsigned int dim = elem->dim();
1878 :
1879 : // One way for implementing this would be to exchange the fe with the FEInterface- class.
1880 : // However, it needs to be discussed whether integral-norms make sense for infinite elements.
1881 : // or in which sense they could make sense.
1882 3459437 : if (elem->infinite() )
1883 0 : libmesh_not_implemented();
1884 :
1885 17138942 : if (skip_dimensions && skip_dimensions->find(dim) != skip_dimensions->end())
1886 25062 : continue;
1887 :
1888 18705353 : QBase * qrule = q_rules[dim].get();
1889 1652935 : libmesh_assert(qrule);
1890 :
1891 17110082 : this->get_dof_map().dof_indices (elem, dof_indices, var);
1892 :
1893 13861876 : auto element_calculation = [&dof_indices, &elem,
1894 : norm_type, norm_weight, norm_weight_sq, &qrule,
1895 792377197 : &local_v, &v_norm](auto & fe) {
1896 : typedef typename std::remove_reference<decltype(fe)>::type::OutputShape OutputShape;
1897 : typedef typename TensorTools::MakeNumber<OutputShape>::type OutputNumberShape;
1898 : typedef typename std::remove_reference<decltype(fe)>::type::OutputGradient OutputGradient;
1899 : typedef typename TensorTools::MakeNumber<OutputGradient>::type OutputNumberGradient;
1900 :
1901 17110082 : const std::vector<Real> & JxW = fe.get_JxW();
1902 1652935 : const std::vector<std::vector<OutputShape>> * phi = nullptr;
1903 17110082 : if (norm_type == H1 ||
1904 15459356 : norm_type == H2 ||
1905 15457975 : norm_type == L2 ||
1906 15457975 : norm_type == L1 ||
1907 : norm_type == L_INF)
1908 1652383 : phi = &(fe.get_phi());
1909 :
1910 1652935 : const std::vector<std::vector<OutputGradient>> * dphi = nullptr;
1911 18705353 : if (norm_type == H1 ||
1912 15459356 : norm_type == H2 ||
1913 15459080 : norm_type == H1_SEMINORM ||
1914 : norm_type == W1_INF_SEMINORM)
1915 1651278 : dphi = &(fe.get_dphi());
1916 :
1917 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1918 : typedef typename std::remove_reference<decltype(fe)>::type::OutputTensor OutputTensor;
1919 :
1920 1652935 : const std::vector<std::vector<OutputTensor>> * d2phi = nullptr;
1921 18705353 : if (norm_type == H2 ||
1922 17110082 : norm_type == H2_SEMINORM ||
1923 : norm_type == W2_INF_SEMINORM)
1924 0 : d2phi = &(fe.get_d2phi());
1925 : #endif
1926 :
1927 17110082 : fe.reinit (elem);
1928 :
1929 17110082 : const unsigned int n_qp = qrule->n_points();
1930 :
1931 : const unsigned int n_sf = cast_int<unsigned int>
1932 3248206 : (dof_indices.size());
1933 :
1934 : // Begin the loop over the Quadrature points.
1935 85755590 : for (unsigned int qp=0; qp<n_qp; qp++)
1936 : {
1937 68645508 : if (norm_type == L1)
1938 : {
1939 0 : OutputNumberShape u_h = 0.;
1940 0 : for (unsigned int i=0; i != n_sf; ++i)
1941 0 : u_h += (*phi)[i][qp] * (*local_v)(dof_indices[i]);
1942 0 : v_norm += norm_weight *
1943 0 : JxW[qp] * TensorTools::norm(u_h);
1944 : }
1945 :
1946 68645508 : if (norm_type == L_INF)
1947 : {
1948 2772 : OutputNumberShape u_h = 0.;
1949 440028 : for (unsigned int i=0; i != n_sf; ++i)
1950 537030 : u_h += (*phi)[i][qp] * (*local_v)(dof_indices[i]);
1951 34529 : v_norm = std::max(v_norm, norm_weight * TensorTools::norm(u_h));
1952 : }
1953 :
1954 68645508 : if (norm_type == H1 ||
1955 62041932 : norm_type == H2 ||
1956 : norm_type == L2)
1957 : {
1958 6620419 : OutputNumberShape u_h = 0.;
1959 344614584 : for (unsigned int i=0; i != n_sf; ++i)
1960 377435308 : u_h += (*phi)[i][qp] * (*local_v)(dof_indices[i]);
1961 74933589 : v_norm += norm_weight_sq *
1962 68543826 : JxW[qp] * TensorTools::norm_sq(u_h);
1963 : }
1964 :
1965 75043587 : if (norm_type == H1 ||
1966 62041932 : norm_type == H2 ||
1967 62016773 : norm_type == H1_SEMINORM)
1968 : {
1969 6606348 : OutputNumberGradient grad_u_h;
1970 338447388 : for (unsigned int i=0; i != n_sf; ++i)
1971 295288586 : grad_u_h.add_scaled((*dphi)[i][qp], (*local_v)(dof_indices[i]));
1972 74750666 : v_norm += norm_weight_sq *
1973 68374974 : JxW[qp] * grad_u_h.norm_sq();
1974 : }
1975 :
1976 68645508 : if (norm_type == W1_INF_SEMINORM)
1977 : {
1978 2772 : OutputNumberGradient grad_u_h;
1979 440028 : for (unsigned int i=0; i != n_sf; ++i)
1980 438858 : grad_u_h.add_scaled((*dphi)[i][qp], (*local_v)(dof_indices[i]));
1981 36441 : v_norm = std::max(v_norm, norm_weight * grad_u_h.norm());
1982 : }
1983 :
1984 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1985 : typedef typename TensorTools::MakeNumber<OutputTensor>::type OutputNumberTensor;
1986 :
1987 75043587 : if (norm_type == H2 ||
1988 62016773 : norm_type == H2_SEMINORM)
1989 : {
1990 0 : OutputNumberTensor hess_u_h;
1991 0 : for (unsigned int i=0; i != n_sf; ++i)
1992 0 : hess_u_h.add_scaled((*d2phi)[i][qp], (*local_v)(dof_indices[i]));
1993 0 : v_norm += norm_weight_sq *
1994 0 : JxW[qp] * hess_u_h.norm_sq();
1995 : }
1996 :
1997 68645508 : if (norm_type == W2_INF_SEMINORM)
1998 : {
1999 0 : OutputNumberTensor hess_u_h;
2000 0 : for (unsigned int i=0; i != n_sf; ++i)
2001 0 : hess_u_h.add_scaled((*d2phi)[i][qp], (*local_v)(dof_indices[i]));
2002 0 : v_norm = std::max(v_norm, norm_weight * hess_u_h.norm());
2003 : }
2004 : #endif
2005 : }
2006 34220164 : };
2007 :
2008 3248206 : FEBase * scalar_fe = fe_ptrs[dim].get();
2009 3248206 : FEVectorBase * vec_fe = vec_fe_ptrs[dim].get();
2010 :
2011 17110082 : if (scalar_fe)
2012 : {
2013 1651830 : libmesh_assert(!vec_fe);
2014 17096822 : element_calculation(*scalar_fe);
2015 : }
2016 :
2017 17110082 : if (vec_fe)
2018 : {
2019 1105 : libmesh_assert(!scalar_fe);
2020 13260 : element_calculation(*vec_fe);
2021 : }
2022 107541 : }
2023 107541 : }
2024 :
2025 114641 : if (using_hilbert_norm)
2026 : {
2027 114041 : this->comm().sum(v_norm);
2028 114041 : v_norm = std::sqrt(v_norm);
2029 : }
2030 : else
2031 : {
2032 600 : this->comm().max(v_norm);
2033 : }
2034 :
2035 114641 : return v_norm;
2036 107541 : }
2037 :
2038 :
2039 :
2040 17521 : std::string System::get_info() const
2041 : {
2042 18537 : std::ostringstream oss;
2043 :
2044 :
2045 508 : const std::string & sys_name = this->name();
2046 :
2047 17521 : oss << " System #" << this->number() << ", \"" << sys_name << "\"\n"
2048 34534 : << " Type \"" << this->system_type() << "\"\n"
2049 50531 : << " Variables=";
2050 :
2051 39914 : for (auto vg : make_range(this->n_variable_groups()))
2052 : {
2053 646 : const VariableGroup & vg_description (this->variable_group(vg));
2054 :
2055 22393 : if (vg_description.n_variables() > 1) oss << "{ ";
2056 55889 : for (auto vn : make_range(vg_description.n_variables()))
2057 65072 : oss << "\"" << vg_description.name(vn) << "\" ";
2058 22393 : if (vg_description.n_variables() > 1) oss << "} ";
2059 : }
2060 :
2061 17521 : oss << '\n';
2062 :
2063 17521 : oss << " Finite Element Types=";
2064 : #ifndef LIBMESH_ENABLE_INFINITE_ELEMENTS
2065 37176 : for (auto vg : make_range(this->n_variable_groups()))
2066 : oss << "\""
2067 41730 : << Utility::enum_to_string<FEFamily>(this->get_dof_map().variable_group(vg).type().family)
2068 41730 : << "\" ";
2069 : #else
2070 2738 : for (auto vg : make_range(this->n_variable_groups()))
2071 : {
2072 : oss << "\""
2073 3056 : << Utility::enum_to_string<FEFamily>(this->get_dof_map().variable_group(vg).type().family)
2074 : << "\", \""
2075 3056 : << Utility::enum_to_string<FEFamily>(this->get_dof_map().variable_group(vg).type().radial_family)
2076 3292 : << "\" ";
2077 : }
2078 :
2079 1210 : oss << '\n' << " Infinite Element Mapping=";
2080 2738 : for (auto vg : make_range(this->n_variable_groups()))
2081 : oss << "\""
2082 3056 : << Utility::enum_to_string<InfMapType>(this->get_dof_map().variable_group(vg).type().inf_map)
2083 2410 : << "\" ";
2084 : #endif
2085 :
2086 17521 : oss << '\n';
2087 :
2088 17521 : oss << " Approximation Orders=";
2089 39914 : for (auto vg : make_range(this->n_variable_groups()))
2090 : {
2091 : #ifndef LIBMESH_ENABLE_INFINITE_ELEMENTS
2092 : oss << "\""
2093 41730 : << Utility::enum_to_string<Order>(this->get_dof_map().variable_group(vg).type().order)
2094 41730 : << "\" ";
2095 : #else
2096 : oss << "\""
2097 3056 : << Utility::enum_to_string<Order>(this->get_dof_map().variable_group(vg).type().order)
2098 : << "\", \""
2099 3056 : << Utility::enum_to_string<Order>(this->get_dof_map().variable_group(vg).type().radial_order)
2100 3292 : << "\" ";
2101 : #endif
2102 : }
2103 :
2104 17521 : oss << '\n';
2105 :
2106 34026 : oss << " n_dofs()=" << this->n_dofs() << '\n';
2107 17521 : dof_id_type local_dofs = this->n_local_dofs();
2108 17521 : oss << " n_local_dofs()=" << local_dofs << '\n';
2109 17521 : this->comm().max(local_dofs);
2110 18029 : oss << " max(n_local_dofs())=" << local_dofs << '\n';
2111 : #ifdef LIBMESH_ENABLE_CONSTRAINTS
2112 34534 : oss << " n_constrained_dofs()=" << this->n_constrained_dofs() << '\n';
2113 34026 : oss << " n_local_constrained_dofs()=" << this->n_local_constrained_dofs() << '\n';
2114 17521 : dof_id_type local_unconstrained_dofs = this->n_local_dofs() - this->n_local_constrained_dofs();
2115 17521 : this->comm().max(local_unconstrained_dofs);
2116 18029 : oss << " max(local unconstrained dofs)=" << local_unconstrained_dofs << '\n';
2117 : #endif
2118 :
2119 34026 : oss << " " << "n_vectors()=" << this->n_vectors() << '\n';
2120 34026 : oss << " " << "n_matrices()=" << this->n_matrices() << '\n';
2121 : // oss << " " << "n_additional_matrices()=" << this->n_additional_matrices() << '\n';
2122 :
2123 34534 : oss << this->get_dof_map().get_info();
2124 :
2125 18029 : return oss.str();
2126 16505 : }
2127 :
2128 :
2129 :
2130 360 : void System::attach_init_function (void fptr(EquationSystems & es,
2131 : const std::string & name))
2132 : {
2133 12 : libmesh_assert(fptr);
2134 :
2135 360 : if (_init_system_object != nullptr)
2136 : {
2137 : libmesh_warning("WARNING: Cannot specify both initialization function and object!");
2138 :
2139 0 : _init_system_object = nullptr;
2140 : }
2141 :
2142 360 : _init_system_function = fptr;
2143 360 : }
2144 :
2145 :
2146 :
2147 280 : void System::attach_init_object (System::Initialization & init_in)
2148 : {
2149 280 : if (_init_system_function != nullptr)
2150 : {
2151 : libmesh_warning("WARNING: Cannot specify both initialization object and function!");
2152 :
2153 0 : _init_system_function = nullptr;
2154 : }
2155 :
2156 280 : _init_system_object = &init_in;
2157 280 : }
2158 :
2159 :
2160 :
2161 9734 : void System::attach_assemble_function (void fptr(EquationSystems & es,
2162 : const std::string & name))
2163 : {
2164 282 : libmesh_assert(fptr);
2165 :
2166 9734 : if (_assemble_system_object != nullptr)
2167 : {
2168 : libmesh_warning("WARNING: Cannot specify both assembly function and object!");
2169 :
2170 0 : _assemble_system_object = nullptr;
2171 : }
2172 :
2173 9734 : _assemble_system_function = fptr;
2174 9734 : }
2175 :
2176 :
2177 :
2178 142 : void System::attach_assemble_object (System::Assembly & assemble_in)
2179 : {
2180 142 : if (_assemble_system_function != nullptr)
2181 : {
2182 : libmesh_warning("WARNING: Cannot specify both assembly object and function!");
2183 :
2184 0 : _assemble_system_function = nullptr;
2185 : }
2186 :
2187 142 : _assemble_system_object = &assemble_in;
2188 142 : }
2189 :
2190 :
2191 :
2192 0 : void System::attach_constraint_function(void fptr(EquationSystems & es,
2193 : const std::string & name))
2194 : {
2195 0 : libmesh_assert(fptr);
2196 :
2197 0 : if (_constrain_system_object != nullptr)
2198 : {
2199 : libmesh_warning("WARNING: Cannot specify both constraint function and object!");
2200 :
2201 0 : _constrain_system_object = nullptr;
2202 : }
2203 :
2204 0 : _constrain_system_function = fptr;
2205 0 : }
2206 :
2207 :
2208 :
2209 923 : void System::attach_constraint_object (System::Constraint & constrain)
2210 : {
2211 923 : if (_constrain_system_function != nullptr)
2212 : {
2213 : libmesh_warning("WARNING: Cannot specify both constraint object and function!");
2214 :
2215 0 : _constrain_system_function = nullptr;
2216 : }
2217 :
2218 923 : _constrain_system_object = &constrain;
2219 923 : }
2220 :
2221 0 : bool System::has_constraint_object () const
2222 : {
2223 0 : return _constrain_system_object != nullptr;
2224 : }
2225 :
2226 0 : System::Constraint& System::get_constraint_object ()
2227 : {
2228 0 : libmesh_assert_msg(_constrain_system_object,"No constraint object available.");
2229 0 : return *_constrain_system_object;
2230 : }
2231 :
2232 :
2233 :
2234 0 : void System::attach_QOI_function(void fptr(EquationSystems &,
2235 : const std::string &,
2236 : const QoISet &))
2237 : {
2238 0 : libmesh_assert(fptr);
2239 :
2240 0 : if (_qoi_evaluate_object != nullptr)
2241 : {
2242 : libmesh_warning("WARNING: Cannot specify both QOI function and object!");
2243 :
2244 0 : _qoi_evaluate_object = nullptr;
2245 : }
2246 :
2247 0 : _qoi_evaluate_function = fptr;
2248 0 : }
2249 :
2250 :
2251 :
2252 0 : void System::attach_QOI_object (QOI & qoi_in)
2253 : {
2254 0 : if (_qoi_evaluate_function != nullptr)
2255 : {
2256 : libmesh_warning("WARNING: Cannot specify both QOI object and function!");
2257 :
2258 0 : _qoi_evaluate_function = nullptr;
2259 : }
2260 :
2261 0 : _qoi_evaluate_object = &qoi_in;
2262 0 : }
2263 :
2264 :
2265 :
2266 0 : void System::attach_QOI_derivative(void fptr(EquationSystems &, const std::string &,
2267 : const QoISet &, bool, bool))
2268 : {
2269 0 : libmesh_assert(fptr);
2270 :
2271 0 : if (_qoi_evaluate_derivative_object != nullptr)
2272 : {
2273 : libmesh_warning("WARNING: Cannot specify both QOI derivative function and object!");
2274 :
2275 0 : _qoi_evaluate_derivative_object = nullptr;
2276 : }
2277 :
2278 0 : _qoi_evaluate_derivative_function = fptr;
2279 0 : }
2280 :
2281 :
2282 :
2283 0 : void System::attach_QOI_derivative_object (QOIDerivative & qoi_derivative)
2284 : {
2285 0 : if (_qoi_evaluate_derivative_function != nullptr)
2286 : {
2287 : libmesh_warning("WARNING: Cannot specify both QOI derivative object and function!");
2288 :
2289 0 : _qoi_evaluate_derivative_function = nullptr;
2290 : }
2291 :
2292 0 : _qoi_evaluate_derivative_object = &qoi_derivative;
2293 0 : }
2294 :
2295 :
2296 :
2297 236514 : void System::user_initialization ()
2298 : {
2299 : // Call the user-provided initialization function,
2300 : // if it was provided
2301 236514 : if (_init_system_function != nullptr)
2302 431 : this->_init_system_function (_equation_systems, this->name());
2303 :
2304 : // ...or the user-provided initialization object.
2305 236083 : else if (_init_system_object != nullptr)
2306 280 : this->_init_system_object->initialize();
2307 236514 : }
2308 :
2309 :
2310 :
2311 39219 : void System::user_assembly ()
2312 : {
2313 : // Call the user-provided assembly function,
2314 : // if it was provided
2315 39219 : if (_assemble_system_function != nullptr)
2316 39077 : this->_assemble_system_function (_equation_systems, this->name());
2317 :
2318 : // ...or the user-provided assembly object.
2319 142 : else if (_assemble_system_object != nullptr)
2320 142 : this->_assemble_system_object->assemble();
2321 39219 : }
2322 :
2323 :
2324 :
2325 261779 : void System::user_constrain ()
2326 : {
2327 : // Call the user-provided constraint function,
2328 : // if it was provided
2329 261779 : if (_constrain_system_function!= nullptr)
2330 0 : this->_constrain_system_function(_equation_systems, this->name());
2331 :
2332 : // ...or the user-provided constraint object.
2333 261779 : else if (_constrain_system_object != nullptr)
2334 923 : this->_constrain_system_object->constrain();
2335 261779 : }
2336 :
2337 :
2338 :
2339 0 : void System::user_QOI (const QoISet & qoi_indices)
2340 : {
2341 : // Call the user-provided quantity of interest function,
2342 : // if it was provided
2343 0 : if (_qoi_evaluate_function != nullptr)
2344 0 : this->_qoi_evaluate_function(_equation_systems, this->name(), qoi_indices);
2345 :
2346 : // ...or the user-provided QOI function object.
2347 0 : else if (_qoi_evaluate_object != nullptr)
2348 0 : this->_qoi_evaluate_object->qoi(qoi_indices);
2349 0 : }
2350 :
2351 :
2352 :
2353 0 : void System::user_QOI_derivative(const QoISet & qoi_indices,
2354 : bool include_liftfunc,
2355 : bool apply_constraints)
2356 : {
2357 : // Call the user-provided quantity of interest derivative,
2358 : // if it was provided
2359 0 : if (_qoi_evaluate_derivative_function != nullptr)
2360 0 : this->_qoi_evaluate_derivative_function
2361 0 : (_equation_systems, this->name(), qoi_indices, include_liftfunc,
2362 : apply_constraints);
2363 :
2364 : // ...or the user-provided QOI derivative function object.
2365 0 : else if (_qoi_evaluate_derivative_object != nullptr)
2366 0 : this->_qoi_evaluate_derivative_object->qoi_derivative
2367 0 : (qoi_indices, include_liftfunc, apply_constraints);
2368 0 : }
2369 :
2370 :
2371 2272 : void System::init_qois(unsigned int n_qois)
2372 : {
2373 2272 : qoi.resize(n_qois);
2374 2272 : qoi_error_estimates.resize(n_qois);
2375 2272 : }
2376 :
2377 :
2378 117555 : void System::set_qoi(unsigned int qoi_index, Number qoi_value)
2379 : {
2380 3358 : libmesh_assert(qoi_index < qoi.size());
2381 :
2382 120913 : qoi[qoi_index] = qoi_value;
2383 117555 : }
2384 :
2385 :
2386 51380 : Number System::get_qoi_value(unsigned int qoi_index) const
2387 : {
2388 1468 : libmesh_assert(qoi_index < qoi.size());
2389 52848 : return qoi[qoi_index];
2390 : }
2391 :
2392 :
2393 54015 : std::vector<Number> System::get_qoi_values() const
2394 : {
2395 54015 : return this->qoi;
2396 : }
2397 :
2398 :
2399 39155 : void System::set_qoi(std::vector<Number> new_qoi)
2400 : {
2401 1118 : libmesh_assert_equal_to(this->qoi.size(), new_qoi.size());
2402 39155 : this->qoi = std::move(new_qoi);
2403 39155 : }
2404 :
2405 :
2406 22400 : void System::set_qoi_error_estimate(unsigned int qoi_index, Number qoi_error_estimate)
2407 : {
2408 640 : libmesh_assert(qoi_index < qoi_error_estimates.size());
2409 :
2410 23040 : qoi_error_estimates[qoi_index] = qoi_error_estimate;
2411 22400 : }
2412 :
2413 22400 : Number System::get_qoi_error_estimate_value(unsigned int qoi_index) const
2414 : {
2415 640 : libmesh_assert(qoi_index < qoi_error_estimates.size());
2416 23040 : return qoi_error_estimates[qoi_index];
2417 : }
2418 :
2419 :
2420 :
2421 285117 : Number System::point_value(unsigned int var,
2422 : const Point & p,
2423 : const bool insist_on_success,
2424 : const NumericVector<Number> *sol) const
2425 : {
2426 : // This function must be called on every processor; there's no
2427 : // telling where in the partition p falls.
2428 7996 : parallel_object_only();
2429 :
2430 : // And every processor had better agree about which point we're
2431 : // looking for
2432 : #ifndef NDEBUG
2433 7996 : libmesh_assert(this->comm().verify(p(0)));
2434 : #if LIBMESH_DIM > 1
2435 7996 : libmesh_assert(this->comm().verify(p(1)));
2436 : #endif
2437 : #if LIBMESH_DIM > 2
2438 7996 : libmesh_assert(this->comm().verify(p(2)));
2439 : #endif
2440 : #endif // NDEBUG
2441 :
2442 : // Get a reference to the mesh object associated with the system object that calls this function
2443 15992 : const MeshBase & mesh = this->get_mesh();
2444 :
2445 : // Use an existing PointLocator or create a new one
2446 285117 : std::unique_ptr<PointLocatorBase> locator_ptr = mesh.sub_point_locator();
2447 7996 : PointLocatorBase & locator = *locator_ptr;
2448 :
2449 285117 : if (!insist_on_success || !mesh.is_serial())
2450 244858 : locator.enable_out_of_mesh_mode();
2451 :
2452 : // Get a pointer to an element that contains p and allows us to
2453 : // evaluate var
2454 : const std::set<subdomain_id_type> & raw_subdomains =
2455 7996 : this->variable(var).active_subdomains();
2456 : const std::set<subdomain_id_type> * implicit_subdomains =
2457 285117 : raw_subdomains.empty() ? nullptr : &raw_subdomains;
2458 285117 : const Elem * e = locator(p, implicit_subdomains);
2459 :
2460 285117 : Number u = 0;
2461 :
2462 285117 : if (e && this->get_dof_map().is_evaluable(*e, var))
2463 100848 : u = point_value(var, p, *e, sol);
2464 :
2465 : // If I have an element containing p, then let's let everyone know
2466 : processor_id_type lowest_owner =
2467 285117 : (e && (e->processor_id() == this->processor_id())) ?
2468 281119 : this->processor_id() : this->n_processors();
2469 285117 : this->comm().min(lowest_owner);
2470 :
2471 : // Everybody should get their value from a processor that was able
2472 : // to compute it.
2473 : // If nobody admits owning the point, we have a problem.
2474 293113 : if (lowest_owner != this->n_processors())
2475 277121 : this->comm().broadcast(u, lowest_owner);
2476 : else
2477 0 : libmesh_assert(!insist_on_success);
2478 :
2479 301109 : return u;
2480 269125 : }
2481 :
2482 119229 : Number System::point_value(unsigned int var,
2483 : const Point & p,
2484 : const Elem & e,
2485 : const NumericVector<Number> *sol) const
2486 : {
2487 : // Ensuring that the given point is really in the element is an
2488 : // expensive assert, but as long as debugging is turned on we might
2489 : // as well try to catch a particularly nasty potential error
2490 6122 : libmesh_assert (e.contains_point(p));
2491 :
2492 119229 : if (!sol)
2493 3460 : sol = this->current_local_solution.get();
2494 :
2495 : // Get the dof map to get the proper indices for our computation
2496 6122 : const DofMap & dof_map = this->get_dof_map();
2497 :
2498 : // Make sure we can evaluate on this element.
2499 6122 : libmesh_assert (dof_map.is_evaluable(e, var));
2500 :
2501 : // Need dof_indices for phi[i][j]
2502 12244 : std::vector<dof_id_type> dof_indices;
2503 :
2504 : // Fill in the dof_indices for our element
2505 119229 : dof_map.dof_indices (&e, dof_indices, var);
2506 :
2507 : // Get the no of dofs associated with this point
2508 : const unsigned int num_dofs = cast_int<unsigned int>
2509 12244 : (dof_indices.size());
2510 :
2511 119229 : FEType fe_type = dof_map.variable_type(var);
2512 :
2513 : // Map the physical co-ordinates to the master co-ordinates
2514 119229 : Point coor = FEMap::inverse_map(e.dim(), &e, p);
2515 :
2516 : // get the shape function value via the FEInterface to also handle the case
2517 : // of infinite elements correctly, the shape function is not fe->phi().
2518 125351 : FEComputeData fe_data(this->get_equation_systems(), coor);
2519 119229 : FEInterface::compute_data(e.dim(), fe_type, &e, fe_data);
2520 :
2521 : // Get ready to accumulate a value
2522 6122 : Number u = 0;
2523 :
2524 2442313 : for (unsigned int l=0; l<num_dofs; l++)
2525 : {
2526 2521150 : u += fe_data.shape[l] * (*sol)(dof_indices[l]);
2527 : }
2528 :
2529 125351 : return u;
2530 106985 : }
2531 :
2532 :
2533 :
2534 18381 : Number System::point_value(unsigned int var, const Point & p, const Elem * e) const
2535 : {
2536 337 : libmesh_assert(e);
2537 18381 : return this->point_value(var, p, *e);
2538 : }
2539 :
2540 :
2541 :
2542 123722 : Number System::point_value(unsigned int var, const Point & p, const NumericVector<Number> * sol) const
2543 : {
2544 123722 : return this->point_value(var, p, true, sol);
2545 : }
2546 :
2547 :
2548 :
2549 :
2550 0 : Gradient System::point_gradient(unsigned int var,
2551 : const Point & p,
2552 : const bool insist_on_success,
2553 : const NumericVector<Number> *sol) const
2554 : {
2555 : // This function must be called on every processor; there's no
2556 : // telling where in the partition p falls.
2557 0 : parallel_object_only();
2558 :
2559 : // And every processor had better agree about which point we're
2560 : // looking for
2561 : #ifndef NDEBUG
2562 0 : libmesh_assert(this->comm().verify(p(0)));
2563 : #if LIBMESH_DIM > 1
2564 0 : libmesh_assert(this->comm().verify(p(1)));
2565 : #endif
2566 : #if LIBMESH_DIM > 2
2567 0 : libmesh_assert(this->comm().verify(p(2)));
2568 : #endif
2569 : #endif // NDEBUG
2570 :
2571 : // Get a reference to the mesh object associated with the system object that calls this function
2572 0 : const MeshBase & mesh = this->get_mesh();
2573 :
2574 : // Use an existing PointLocator or create a new one
2575 0 : std::unique_ptr<PointLocatorBase> locator_ptr = mesh.sub_point_locator();
2576 0 : PointLocatorBase & locator = *locator_ptr;
2577 :
2578 0 : if (!insist_on_success || !mesh.is_serial())
2579 0 : locator.enable_out_of_mesh_mode();
2580 :
2581 : // Get a pointer to an element that contains p and allows us to
2582 : // evaluate var
2583 : const std::set<subdomain_id_type> & raw_subdomains =
2584 0 : this->variable(var).active_subdomains();
2585 : const std::set<subdomain_id_type> * implicit_subdomains =
2586 0 : raw_subdomains.empty() ? nullptr : &raw_subdomains;
2587 0 : const Elem * e = locator(p, implicit_subdomains);
2588 :
2589 0 : Gradient grad_u;
2590 :
2591 0 : if (e && this->get_dof_map().is_evaluable(*e, var))
2592 0 : grad_u = point_gradient(var, p, *e, sol);
2593 :
2594 : // If I have an element containing p, then let's let everyone know
2595 : processor_id_type lowest_owner =
2596 0 : (e && (e->processor_id() == this->processor_id())) ?
2597 0 : this->processor_id() : this->n_processors();
2598 0 : this->comm().min(lowest_owner);
2599 :
2600 : // Everybody should get their value from a processor that was able
2601 : // to compute it.
2602 : // If nobody admits owning the point, we may have a problem.
2603 0 : if (lowest_owner != this->n_processors())
2604 0 : this->comm().broadcast(grad_u, lowest_owner);
2605 : else
2606 0 : libmesh_assert(!insist_on_success);
2607 :
2608 0 : return grad_u;
2609 0 : }
2610 :
2611 :
2612 0 : Gradient System::point_gradient(unsigned int var,
2613 : const Point & p,
2614 : const Elem & e,
2615 : const NumericVector<Number> *sol) const
2616 : {
2617 : // Ensuring that the given point is really in the element is an
2618 : // expensive assert, but as long as debugging is turned on we might
2619 : // as well try to catch a particularly nasty potential error
2620 0 : libmesh_assert (e.contains_point(p));
2621 :
2622 0 : if (!sol)
2623 0 : sol = this->current_local_solution.get();
2624 :
2625 : // Get the dof map to get the proper indices for our computation
2626 0 : const DofMap & dof_map = this->get_dof_map();
2627 :
2628 : // write the element dimension into a separate variable.
2629 0 : const unsigned int dim = e.dim();
2630 :
2631 : // Make sure we can evaluate on this element.
2632 0 : libmesh_assert (dof_map.is_evaluable(e, var));
2633 :
2634 : // Need dof_indices for phi[i][j]
2635 0 : std::vector<dof_id_type> dof_indices;
2636 :
2637 : // Fill in the dof_indices for our element
2638 0 : dof_map.dof_indices (&e, dof_indices, var);
2639 :
2640 : // Get the no of dofs associated with this point
2641 : const unsigned int num_dofs = cast_int<unsigned int>
2642 0 : (dof_indices.size());
2643 :
2644 0 : FEType fe_type = dof_map.variable_type(var);
2645 :
2646 : // Map the physical co-ordinates to the master co-ordinates
2647 0 : Point coor = FEMap::inverse_map(dim, &e, p);
2648 :
2649 : // get the shape function value via the FEInterface to also handle the case
2650 : // of infinite elements correctly, the shape function is not fe->phi().
2651 0 : FEComputeData fe_data(this->get_equation_systems(), coor);
2652 0 : fe_data.enable_derivative();
2653 0 : FEInterface::compute_data(dim, fe_type, &e, fe_data);
2654 :
2655 : // Get ready to accumulate a gradient
2656 0 : Gradient grad_u;
2657 :
2658 0 : for (unsigned int l=0; l<num_dofs; l++)
2659 : {
2660 : // Chartesian coordinates have always LIBMESH_DIM entries,
2661 : // local coordinates have as many coordinates as the element has.
2662 0 : for (std::size_t v=0; v<dim; v++)
2663 0 : for (std::size_t xyz=0; xyz<LIBMESH_DIM; xyz++)
2664 : {
2665 : // FIXME: this needs better syntax: It is matrix-vector multiplication.
2666 0 : grad_u(xyz) += fe_data.local_transform[v][xyz]
2667 0 : * fe_data.dshape[l](v)
2668 0 : * (*sol)(dof_indices[l]);
2669 : }
2670 : }
2671 :
2672 0 : return grad_u;
2673 0 : }
2674 :
2675 :
2676 :
2677 0 : Gradient System::point_gradient(unsigned int var, const Point & p, const Elem * e) const
2678 : {
2679 0 : libmesh_assert(e);
2680 0 : return this->point_gradient(var, p, *e);
2681 : }
2682 :
2683 :
2684 :
2685 0 : Gradient System::point_gradient(unsigned int var, const Point & p, const NumericVector<Number> * sol) const
2686 : {
2687 0 : return this->point_gradient(var, p, true, sol);
2688 : }
2689 :
2690 :
2691 :
2692 : // We can only accumulate a hessian with --enable-second
2693 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
2694 0 : Tensor System::point_hessian(unsigned int var,
2695 : const Point & p,
2696 : const bool insist_on_success,
2697 : const NumericVector<Number> *sol) const
2698 : {
2699 : // This function must be called on every processor; there's no
2700 : // telling where in the partition p falls.
2701 0 : parallel_object_only();
2702 :
2703 : // And every processor had better agree about which point we're
2704 : // looking for
2705 : #ifndef NDEBUG
2706 0 : libmesh_assert(this->comm().verify(p(0)));
2707 : #if LIBMESH_DIM > 1
2708 0 : libmesh_assert(this->comm().verify(p(1)));
2709 : #endif
2710 : #if LIBMESH_DIM > 2
2711 0 : libmesh_assert(this->comm().verify(p(2)));
2712 : #endif
2713 : #endif // NDEBUG
2714 :
2715 : // Get a reference to the mesh object associated with the system object that calls this function
2716 0 : const MeshBase & mesh = this->get_mesh();
2717 :
2718 : // Use an existing PointLocator or create a new one
2719 0 : std::unique_ptr<PointLocatorBase> locator_ptr = mesh.sub_point_locator();
2720 0 : PointLocatorBase & locator = *locator_ptr;
2721 :
2722 0 : if (!insist_on_success || !mesh.is_serial())
2723 0 : locator.enable_out_of_mesh_mode();
2724 :
2725 : // Get a pointer to an element that contains p and allows us to
2726 : // evaluate var
2727 : const std::set<subdomain_id_type> & raw_subdomains =
2728 0 : this->variable(var).active_subdomains();
2729 : const std::set<subdomain_id_type> * implicit_subdomains =
2730 0 : raw_subdomains.empty() ? nullptr : &raw_subdomains;
2731 0 : const Elem * e = locator(p, implicit_subdomains);
2732 :
2733 0 : Tensor hess_u;
2734 :
2735 0 : if (e && this->get_dof_map().is_evaluable(*e, var))
2736 0 : hess_u = point_hessian(var, p, *e, sol);
2737 :
2738 : // If I have an element containing p, then let's let everyone know
2739 : processor_id_type lowest_owner =
2740 0 : (e && (e->processor_id() == this->processor_id())) ?
2741 0 : this->processor_id() : this->n_processors();
2742 0 : this->comm().min(lowest_owner);
2743 :
2744 : // Everybody should get their value from a processor that was able
2745 : // to compute it.
2746 : // If nobody admits owning the point, we may have a problem.
2747 0 : if (lowest_owner != this->n_processors())
2748 0 : this->comm().broadcast(hess_u, lowest_owner);
2749 : else
2750 0 : libmesh_assert(!insist_on_success);
2751 :
2752 0 : return hess_u;
2753 0 : }
2754 :
2755 0 : Tensor System::point_hessian(unsigned int var,
2756 : const Point & p,
2757 : const Elem & e,
2758 : const NumericVector<Number> *sol) const
2759 : {
2760 : // Ensuring that the given point is really in the element is an
2761 : // expensive assert, but as long as debugging is turned on we might
2762 : // as well try to catch a particularly nasty potential error
2763 0 : libmesh_assert (e.contains_point(p));
2764 :
2765 0 : if (!sol)
2766 0 : sol = this->current_local_solution.get();
2767 :
2768 0 : if (e.infinite())
2769 0 : libmesh_not_implemented();
2770 :
2771 : // Get the dof map to get the proper indices for our computation
2772 0 : const DofMap & dof_map = this->get_dof_map();
2773 :
2774 : // Make sure we can evaluate on this element.
2775 0 : libmesh_assert (dof_map.is_evaluable(e, var));
2776 :
2777 : // Need dof_indices for phi[i][j]
2778 0 : std::vector<dof_id_type> dof_indices;
2779 :
2780 : // Fill in the dof_indices for our element
2781 0 : dof_map.dof_indices (&e, dof_indices, var);
2782 :
2783 : // Get the no of dofs associated with this point
2784 : const unsigned int num_dofs = cast_int<unsigned int>
2785 0 : (dof_indices.size());
2786 :
2787 0 : FEType fe_type = dof_map.variable_type(var);
2788 :
2789 : // Build a FE again so we can calculate u(p)
2790 0 : std::unique_ptr<FEBase> fe (FEBase::build(e.dim(), fe_type));
2791 :
2792 : // Map the physical co-ordinates to the master co-ordinates
2793 : // Build a vector of point co-ordinates to send to reinit
2794 0 : std::vector<Point> coor(1, FEMap::inverse_map(e.dim(), &e, p));
2795 :
2796 : // Get the values of the shape function derivatives
2797 0 : const std::vector<std::vector<RealTensor>> & d2phi = fe->get_d2phi();
2798 :
2799 : // Reinitialize the element and compute the shape function values at coor
2800 0 : fe->reinit (&e, &coor);
2801 :
2802 : // Get ready to accumulate a hessian
2803 0 : Tensor hess_u;
2804 :
2805 0 : for (unsigned int l=0; l<num_dofs; l++)
2806 : {
2807 0 : hess_u.add_scaled (d2phi[l][0], (*sol)(dof_indices[l]));
2808 : }
2809 :
2810 0 : return hess_u;
2811 0 : }
2812 :
2813 :
2814 :
2815 0 : Tensor System::point_hessian(unsigned int var, const Point & p, const Elem * e) const
2816 : {
2817 0 : libmesh_assert(e);
2818 0 : return this->point_hessian(var, p, *e);
2819 : }
2820 :
2821 :
2822 :
2823 0 : Tensor System::point_hessian(unsigned int var, const Point & p, const NumericVector<Number> * sol) const
2824 : {
2825 0 : return this->point_hessian(var, p, true, sol);
2826 : }
2827 :
2828 : #else
2829 :
2830 : Tensor System::point_hessian(unsigned int, const Point &, const bool,
2831 : const NumericVector<Number> *) const
2832 : {
2833 : libmesh_error_msg("We can only accumulate a hessian with --enable-second");
2834 :
2835 : // Avoid compiler warnings
2836 : return Tensor();
2837 : }
2838 :
2839 : Tensor System::point_hessian(unsigned int, const Point &, const Elem &,
2840 : const NumericVector<Number> *) const
2841 : {
2842 : libmesh_error_msg("We can only accumulate a hessian with --enable-second");
2843 :
2844 : // Avoid compiler warnings
2845 : return Tensor();
2846 : }
2847 :
2848 : Tensor System::point_hessian(unsigned int, const Point &, const Elem *) const
2849 : {
2850 : libmesh_error_msg("We can only accumulate a hessian with --enable-second");
2851 :
2852 : // Avoid compiler warnings
2853 : return Tensor();
2854 : }
2855 :
2856 : Tensor System::point_hessian(unsigned int, const Point &, const NumericVector<Number> *) const
2857 : {
2858 : libmesh_error_msg("We can only accumulate a hessian with --enable-second");
2859 :
2860 : // Avoid compiler warnings
2861 : return Tensor();
2862 : }
2863 :
2864 : #endif // LIBMESH_ENABLE_SECOND_DERIVATIVES
2865 :
2866 350 : void System::create_static_condensation()
2867 : {
2868 350 : this->get_dof_map().create_static_condensation(this->get_mesh(), *this);
2869 350 : }
2870 :
2871 50568 : bool System::has_static_condensation() const
2872 : {
2873 50568 : return this->get_dof_map().has_static_condensation();
2874 : }
2875 :
2876 : } // namespace libMesh
|
