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