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