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