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 :
43 : // includes for calculate_norm, point_*
44 : #include "libmesh/fe_base.h"
45 : #include "libmesh/fe_interface.h"
46 : #include "libmesh/parallel.h"
47 : #include "libmesh/parallel_algebra.h"
48 : #include "libmesh/quadrature.h"
49 : #include "libmesh/tensor_value.h"
50 : #include "libmesh/vector_value.h"
51 : #include "libmesh/tensor_tools.h"
52 : #include "libmesh/enum_norm_type.h"
53 : #include "libmesh/enum_fe_family.h"
54 :
55 : // C++ includes
56 : #include <sstream> // for std::ostringstream
57 :
58 : namespace libMesh
59 : {
60 :
61 :
62 : // ------------------------------------------------------------
63 : // System implementation
64 242602 : System::System (EquationSystems & es,
65 : const std::string & name_in,
66 242602 : const unsigned int number_in) :
67 :
68 : ParallelObject (es),
69 228826 : assemble_before_solve (true),
70 228826 : use_fixed_solution (false),
71 228826 : extra_quadrature_order (0),
72 228826 : solution (NumericVector<Number>::build(this->comm())),
73 228826 : current_local_solution (NumericVector<Number>::build(this->comm())),
74 228826 : time (0.),
75 228826 : _init_system_function (nullptr),
76 228826 : _init_system_object (nullptr),
77 228826 : _assemble_system_function (nullptr),
78 228826 : _assemble_system_object (nullptr),
79 228826 : _constrain_system_function (nullptr),
80 228826 : _constrain_system_object (nullptr),
81 228826 : _qoi_evaluate_function (nullptr),
82 228826 : _qoi_evaluate_object (nullptr),
83 228826 : _qoi_evaluate_derivative_function (nullptr),
84 228826 : _qoi_evaluate_derivative_object (nullptr),
85 228826 : _dof_map (std::make_unique<DofMap>(number_in, es.get_mesh())),
86 228826 : _equation_systems (es),
87 249490 : _mesh (es.get_mesh()),
88 228826 : _sys_name (name_in),
89 228826 : _sys_number (number_in),
90 228826 : _active (true),
91 228826 : _matrices_initialized (false),
92 228826 : _solution_projection (true),
93 228826 : _basic_system_only (false),
94 228826 : _is_initialized (false),
95 228826 : _additional_data_written (false),
96 228826 : adjoint_already_solved (false),
97 228826 : _hide_output (false),
98 228826 : project_with_constraints (true),
99 228826 : _prefer_hash_table_matrix_assembly(false),
100 228826 : _require_sparsity_pattern (false),
101 249490 : _prefix_with_name (false)
102 : {
103 478316 : if (libMesh::on_command_line("--solver-system-names"))
104 0 : this->prefix_with_name(true);
105 707142 : if (libMesh::on_command_line("--" + name_in + "-static-condensation"))
106 420 : this->create_static_condensation();
107 242602 : }
108 :
109 :
110 :
111 679434 : System::~System ()
112 : {
113 6888 : libmesh_exceptionless_assert (!libMesh::closed());
114 450608 : }
115 :
116 :
117 :
118 1423402 : dof_id_type System::n_dofs() const
119 : {
120 1423402 : return _dof_map->n_dofs();
121 : }
122 :
123 :
124 :
125 50957 : dof_id_type System::n_constrained_dofs() const
126 : {
127 : #ifdef LIBMESH_ENABLE_CONSTRAINTS
128 :
129 50957 : return _dof_map->n_constrained_dofs();
130 :
131 : #else
132 :
133 : return 0;
134 :
135 : #endif
136 : }
137 :
138 :
139 :
140 20102 : dof_id_type System::n_local_constrained_dofs() const
141 : {
142 : #ifdef LIBMESH_ENABLE_CONSTRAINTS
143 :
144 20102 : return _dof_map->n_local_constrained_dofs();
145 :
146 : #else
147 :
148 : return 0;
149 :
150 : #endif
151 : }
152 :
153 :
154 :
155 1404299 : dof_id_type System::n_local_dofs() const
156 : {
157 1404299 : return _dof_map->n_local_dofs();
158 : }
159 :
160 :
161 :
162 1076433086 : Number System::current_solution (const dof_id_type global_dof_number) const
163 : {
164 : // Check the sizes
165 96669372 : libmesh_assert_less (global_dof_number, _dof_map->n_dofs());
166 96669372 : libmesh_assert_less (global_dof_number, current_local_solution->size());
167 :
168 1076433086 : return (*current_local_solution)(global_dof_number);
169 : }
170 :
171 :
172 :
173 919 : void System::clear ()
174 : {
175 919 : _dof_map->clear ();
176 919 : solution->clear ();
177 919 : current_local_solution->clear ();
178 :
179 : // clear any user-added vectors
180 28 : _vectors.clear();
181 28 : _vector_projections.clear();
182 28 : _vector_is_adjoint.clear();
183 919 : _is_initialized = false;
184 :
185 : // clear any user-added matrices
186 28 : _matrices.clear();
187 919 : _matrices_initialized = false;
188 :
189 : // But our "basic"/"null" state may still have a StaticCondensation
190 2673 : if (libMesh::on_command_line("--" + _sys_name + "-static-condensation"))
191 70 : this->create_static_condensation();
192 919 : }
193 :
194 :
195 :
196 633 : void System::init ()
197 : {
198 : // Calling init() twice on the same system currently works evil
199 : // magic, whether done directly or via EquationSystems::read()
200 18 : libmesh_assert(!this->is_initialized());
201 :
202 633 : this->reinit_mesh();
203 633 : }
204 :
205 :
206 :
207 242881 : void System::init_data ()
208 : {
209 6886 : parallel_object_only();
210 :
211 13772 : MeshBase & mesh = this->get_mesh();
212 :
213 : // Distribute the degrees of freedom on the mesh
214 242881 : auto total_dofs = _dof_map->distribute_dofs (mesh);
215 :
216 : // Throw an error if the total number of DOFs is not capable of
217 : // being indexed by our solution vector.
218 242858 : auto max_allowed_id = solution->max_allowed_id();
219 242858 : libmesh_error_msg_if(total_dofs > max_allowed_id,
220 : "Cannot allocate a NumericVector with " << total_dofs << " degrees of freedom. "
221 : "The vector can only index up to " << max_allowed_id << " entries.");
222 :
223 : // Recreate any user or internal constraints
224 242858 : this->reinit_constraints();
225 :
226 : // Even if there weren't any constraint changes,
227 : // reinit_constraints() did prepare_send_list() for us.
228 :
229 : // Now finally after dof distribution and construction of any
230 : // possible constraints, we may init any static condensation
231 : // data
232 242787 : _dof_map->reinit_static_condensation();
233 :
234 : // Resize the solution conformal to the current mesh
235 242787 : solution->init (this->n_dofs(), this->n_local_dofs(), false, PARALLEL);
236 :
237 : // Resize the current_local_solution for the current mesh
238 : #ifdef LIBMESH_ENABLE_GHOSTED
239 249669 : current_local_solution->init (this->n_dofs(), this->n_local_dofs(),
240 : _dof_map->get_send_list(), /*fast=*/false,
241 13764 : GHOSTED);
242 : #else
243 : current_local_solution->init (this->n_dofs(), false, SERIAL);
244 : #endif
245 :
246 : // from now on, adding additional vectors or variables can't be done
247 : // without immediately initializing them
248 242787 : _is_initialized = true;
249 :
250 : // initialize & zero other vectors, if necessary
251 287429 : for (auto & [vec_name, vec] : _vectors)
252 : {
253 1330 : libmesh_ignore(vec_name); // spurious warning from old gcc
254 44642 : const ParallelType type = vec->type();
255 :
256 44642 : if (type == GHOSTED)
257 : {
258 : #ifdef LIBMESH_ENABLE_GHOSTED
259 6112 : vec->init (this->n_dofs(), this->n_local_dofs(),
260 : _dof_map->get_send_list(), /*fast=*/false,
261 344 : GHOSTED);
262 : #else
263 : libmesh_error_msg("Cannot initialize ghosted vectors when they are not enabled.");
264 : #endif
265 : }
266 38702 : else if (type == SERIAL)
267 : {
268 0 : vec->init (this->n_dofs(), false, type);
269 : }
270 : else
271 : {
272 1158 : libmesh_assert_equal_to(type, PARALLEL);
273 38702 : vec->init (this->n_dofs(), this->n_local_dofs(), false, type);
274 : }
275 : }
276 :
277 : // Add matrices
278 242787 : this->add_matrices();
279 :
280 : // Clear any existing matrices
281 263651 : for (auto & pr : _matrices)
282 20864 : pr.second->clear();
283 :
284 : // Initialize the matrices for the system
285 242787 : if (!_basic_system_only)
286 242787 : this->init_matrices();
287 242787 : }
288 :
289 242881 : void System::reinit_mesh ()
290 : {
291 6886 : parallel_object_only();
292 :
293 : // First initialize any required data:
294 : // either only the basic System data
295 242881 : if (_basic_system_only)
296 0 : System::init_data();
297 : // or all the derived class' data too
298 : else
299 242881 : this->init_data();
300 :
301 : // If no variables have been added to this system
302 : // don't do anything
303 242787 : if (!this->n_vars())
304 28 : return;
305 :
306 : // Then call the user-provided initialization function
307 241804 : this->user_initialization();
308 :
309 : }
310 :
311 242787 : void System::init_matrices ()
312 : {
313 6882 : parallel_object_only();
314 :
315 : // No matrices to init
316 242787 : if (_matrices.empty())
317 : {
318 : // any future matrices to be added will need their own
319 : // initialization
320 222562 : _matrices_initialized = true;
321 :
322 222562 : return;
323 : }
324 :
325 : // Check for quick return in case the first matrix
326 : // (and by extension all the matrices) has already
327 : // been initialized
328 20225 : if (_matrices.begin()->second->initialized())
329 : {
330 0 : libmesh_assert(_matrices_initialized);
331 0 : return;
332 : }
333 :
334 20225 : _matrices_initialized = true;
335 :
336 : // Tell the matrices about the dof map, and vice versa
337 41089 : for (auto & pr : _matrices)
338 : {
339 580 : SparseMatrix<Number> & m = *(pr.second);
340 580 : libmesh_assert (!m.initialized());
341 :
342 : // We want to allow repeated init() on systems, but we don't
343 : // want to attach the same matrix to the DofMap twice
344 20864 : if (!this->get_dof_map().is_attached(m))
345 20864 : this->get_dof_map().attach_matrix(m);
346 :
347 : // If the user has already explicitly requested that this matrix use a hash table, then we
348 : // always honor that
349 : const bool use_hash =
350 21444 : pr.second->use_hash_table() ||
351 20798 : (this->_prefer_hash_table_matrix_assembly && pr.second->supports_hash_table());
352 20864 : pr.second->use_hash_table(use_hash);
353 : // Make this call after we've determined whether the matrix is using a hash table
354 20864 : if (pr.second->require_sparsity_pattern())
355 20242 : this->_require_sparsity_pattern = true;
356 : }
357 :
358 : // Compute the sparsity pattern for the current
359 : // mesh and DOF distribution. This also updates
360 : // additional matrices, \p DofMap now knows them
361 20225 : if (this->_require_sparsity_pattern)
362 19669 : this->get_dof_map().compute_sparsity(this->get_mesh());
363 :
364 : // Initialize matrices and set to zero
365 41089 : for (auto & [name, mat] : _matrices)
366 : {
367 20864 : mat->init(_matrix_types[name]);
368 20864 : mat->zero();
369 : }
370 : }
371 :
372 :
373 :
374 31142 : void System::restrict_vectors ()
375 : {
376 1000 : parallel_object_only();
377 :
378 : #ifdef LIBMESH_ENABLE_AMR
379 : // Restrict the _vectors on the coarsened cells
380 90975 : for (auto & [vec_name, vec] : _vectors)
381 : {
382 2130 : NumericVector<Number> * v = vec.get();
383 :
384 59833 : if (_vector_projections[vec_name])
385 : {
386 53730 : this->project_vector (*v, this->vector_is_adjoint(vec_name));
387 : }
388 : else
389 : {
390 32458 : const ParallelType type = vec->type();
391 :
392 32458 : if (type == GHOSTED)
393 : {
394 : #ifdef LIBMESH_ENABLE_GHOSTED
395 0 : vec->init (this->n_dofs(), this->n_local_dofs(),
396 : _dof_map->get_send_list(), /*fast=*/false,
397 0 : GHOSTED);
398 : #else
399 : libmesh_error_msg("Cannot initialize ghosted vectors when they are not enabled.");
400 : #endif
401 : }
402 : else
403 32458 : vec->init (this->n_dofs(), this->n_local_dofs(), false, type);
404 : }
405 : }
406 :
407 1000 : const std::vector<dof_id_type> & send_list = _dof_map->get_send_list ();
408 :
409 : // Restrict the solution on the coarsened cells
410 31142 : if (_solution_projection)
411 59604 : this->project_vector (*solution);
412 : // Or at least make sure the solution vector is the correct size
413 : else
414 852 : solution->init (this->n_dofs(), this->n_local_dofs(), true, PARALLEL);
415 :
416 : #ifdef LIBMESH_ENABLE_GHOSTED
417 32142 : current_local_solution->init(this->n_dofs(),
418 : this->n_local_dofs(), send_list,
419 2000 : false, GHOSTED);
420 : #else
421 : current_local_solution->init(this->n_dofs());
422 : #endif
423 :
424 31142 : if (_solution_projection)
425 31266 : solution->localize (*current_local_solution, send_list);
426 :
427 : #endif // LIBMESH_ENABLE_AMR
428 31142 : }
429 :
430 :
431 :
432 31142 : void System::prolong_vectors ()
433 : {
434 : #ifdef LIBMESH_ENABLE_AMR
435 : // Currently project_vector handles both restriction and prolongation
436 31142 : this->restrict_vectors();
437 : #endif
438 31142 : }
439 :
440 :
441 :
442 30000 : void System::reinit ()
443 : {
444 966 : parallel_object_only();
445 :
446 : // project_vector handles vector initialization now
447 966 : libmesh_assert_equal_to (solution->size(), current_local_solution->size());
448 :
449 : // Make sure our static condensation dof map is up-to-date before we init any
450 : // static condensation matrices
451 30000 : this->get_dof_map().reinit_static_condensation();
452 :
453 30000 : if (!_matrices.empty() && !_basic_system_only)
454 : {
455 : // Clear the matrices
456 46814 : for (auto & pr : _matrices)
457 : {
458 23477 : pr.second->clear();
459 23477 : pr.second->attach_dof_map(this->get_dof_map());
460 : }
461 :
462 23337 : if (this->_require_sparsity_pattern)
463 : {
464 : // Clear the sparsity pattern
465 19697 : this->get_dof_map().clear_sparsity();
466 :
467 : // Compute the sparsity pattern for the current
468 : // mesh and DOF distribution. This also updates
469 : // additional matrices, \p DofMap now knows them
470 19697 : this->get_dof_map().compute_sparsity (this->get_mesh());
471 : }
472 :
473 : // Initialize matrices and set to zero
474 46814 : for (auto & pr : _matrices)
475 : {
476 23477 : pr.second->init();
477 23477 : pr.second->zero();
478 : }
479 : }
480 30000 : }
481 :
482 :
483 274280 : void System::reinit_constraints()
484 : {
485 7892 : parallel_object_only();
486 :
487 : #ifdef LIBMESH_ENABLE_CONSTRAINTS
488 274280 : get_dof_map().create_dof_constraints(_mesh, this->time);
489 274280 : user_constrain();
490 274280 : get_dof_map().process_constraints(_mesh);
491 540528 : if (libMesh::on_command_line ("--print-constraints"))
492 0 : get_dof_map().print_dof_constraints(libMesh::out);
493 : #endif
494 274209 : get_dof_map().prepare_send_list();
495 274209 : }
496 :
497 :
498 1463490 : void System::update ()
499 : {
500 37556 : parallel_object_only();
501 :
502 37556 : libmesh_assert(solution->closed());
503 :
504 37556 : const std::vector<dof_id_type> & send_list = _dof_map->get_send_list ();
505 :
506 : // Check sizes
507 37556 : libmesh_assert_equal_to (current_local_solution->size(), solution->size());
508 : // More processors than elements => empty send_list
509 : // libmesh_assert (!send_list.empty());
510 37556 : libmesh_assert_less_equal (send_list.size(), solution->size());
511 :
512 : // Create current_local_solution from solution. This will
513 : // put a local copy of solution into current_local_solution.
514 : // Only the necessary values (specified by the send_list)
515 : // are copied to minimize communication
516 1500230 : solution->localize (*current_local_solution, send_list);
517 1463490 : }
518 :
519 :
520 :
521 29935 : void System::re_update ()
522 : {
523 966 : parallel_object_only();
524 :
525 : // If this system is empty... don't do anything!
526 29935 : if (!this->n_vars())
527 108 : return;
528 :
529 858 : const std::vector<dof_id_type> & send_list = this->get_dof_map().get_send_list ();
530 :
531 : // Check sizes
532 858 : libmesh_assert_equal_to (current_local_solution->size(), solution->size());
533 : // Not true with ghosted vectors
534 : // libmesh_assert_equal_to (current_local_solution->local_size(), solution->size());
535 : // libmesh_assert (!send_list.empty());
536 858 : libmesh_assert_less_equal (send_list.size(), solution->size());
537 :
538 : // Create current_local_solution from solution. This will
539 : // put a local copy of solution into current_local_solution.
540 27011 : solution->localize (*current_local_solution, send_list);
541 : }
542 :
543 :
544 :
545 0 : void System::restrict_solve_to (const SystemSubset * subset,
546 : const SubsetSolveMode /*subset_solve_mode*/)
547 : {
548 0 : if (subset != nullptr)
549 0 : libmesh_not_implemented();
550 0 : }
551 :
552 :
553 :
554 46331 : void System::assemble ()
555 : {
556 : // Log how long the user's assembly code takes
557 3088 : LOG_SCOPE("assemble()", "System");
558 :
559 : // Call the user-specified assembly function
560 46331 : this->user_assembly();
561 46331 : }
562 :
563 :
564 :
565 0 : void System::assemble_qoi (const QoISet & qoi_indices)
566 : {
567 : // Log how long the user's assembly code takes
568 0 : LOG_SCOPE("assemble_qoi()", "System");
569 :
570 : // Call the user-specified quantity of interest function
571 0 : this->user_QOI(qoi_indices);
572 0 : }
573 :
574 :
575 :
576 0 : void System::assemble_qoi_derivative(const QoISet & qoi_indices,
577 : bool include_liftfunc,
578 : bool apply_constraints)
579 : {
580 : // Log how long the user's assembly code takes
581 0 : LOG_SCOPE("assemble_qoi_derivative()", "System");
582 :
583 : // Call the user-specified quantity of interest function
584 0 : this->user_QOI_derivative(qoi_indices, include_liftfunc,
585 0 : apply_constraints);
586 0 : }
587 :
588 :
589 :
590 0 : void System::qoi_parameter_sensitivity (const QoISet & qoi_indices,
591 : const ParameterVector & parameters_vec,
592 : SensitivityData & sensitivities)
593 : {
594 : // Forward sensitivities are more efficient for Nq > Np
595 0 : if (qoi_indices.size(*this) > parameters_vec.size())
596 0 : forward_qoi_parameter_sensitivity(qoi_indices, parameters_vec, sensitivities);
597 : // Adjoint sensitivities are more efficient for Np > Nq,
598 : // and an adjoint may be more reusable than a forward
599 : // solution sensitivity in the Np == Nq case.
600 : else
601 0 : adjoint_qoi_parameter_sensitivity(qoi_indices, parameters_vec, sensitivities);
602 0 : }
603 :
604 :
605 :
606 0 : bool System::compare (const System & other_system,
607 : const Real threshold,
608 : const bool verbose) const
609 : {
610 : // we do not care for matrices, but for vectors
611 0 : libmesh_assert (_is_initialized);
612 0 : libmesh_assert (other_system._is_initialized);
613 :
614 0 : if (verbose)
615 : {
616 0 : libMesh::out << " Systems \"" << _sys_name << "\"" << std::endl;
617 0 : libMesh::out << " comparing matrices not supported." << std::endl;
618 0 : libMesh::out << " comparing names...";
619 : }
620 :
621 : // compare the name: 0 means identical
622 0 : const int name_result = _sys_name.compare(other_system.name());
623 0 : if (verbose)
624 : {
625 0 : if (name_result == 0)
626 0 : libMesh::out << " identical." << std::endl;
627 : else
628 0 : libMesh::out << " names not identical." << std::endl;
629 0 : libMesh::out << " comparing solution vector...";
630 : }
631 :
632 :
633 : // compare the solution: -1 means identical
634 0 : const int solu_result = solution->compare (*other_system.solution.get(),
635 0 : threshold);
636 :
637 0 : if (verbose)
638 : {
639 0 : if (solu_result == -1)
640 0 : libMesh::out << " identical up to threshold." << std::endl;
641 : else
642 0 : libMesh::out << " first difference occurred at index = "
643 0 : << solu_result << "." << std::endl;
644 : }
645 :
646 :
647 : // safety check, whether we handle at least the same number
648 : // of vectors
649 0 : std::vector<int> ov_result;
650 :
651 0 : if (this->n_vectors() != other_system.n_vectors())
652 : {
653 0 : if (verbose)
654 : {
655 0 : libMesh::out << " Fatal difference. This system handles "
656 0 : << this->n_vectors() << " add'l vectors," << std::endl
657 0 : << " while the other system handles "
658 0 : << other_system.n_vectors()
659 0 : << " add'l vectors." << std::endl
660 0 : << " Aborting comparison." << std::endl;
661 : }
662 0 : return false;
663 : }
664 0 : else if (this->n_vectors() == 0)
665 : {
666 : // there are no additional vectors...
667 0 : ov_result.clear ();
668 : }
669 : else
670 : {
671 : // compare other vectors
672 0 : for (auto & [vec_name, vec] : _vectors)
673 : {
674 0 : if (verbose)
675 0 : libMesh::out << " comparing vector \""
676 0 : << vec_name << "\" ...";
677 :
678 : // assume they have the same name
679 : const NumericVector<Number> & other_system_vector =
680 0 : other_system.get_vector(vec_name);
681 :
682 0 : ov_result.push_back(vec->compare(other_system_vector, threshold));
683 :
684 0 : if (verbose)
685 : {
686 0 : if (ov_result[ov_result.size()-1] == -1)
687 0 : libMesh::out << " identical up to threshold." << std::endl;
688 : else
689 0 : libMesh::out << " first difference occurred at" << std::endl
690 0 : << " index = " << ov_result[ov_result.size()-1] << "." << std::endl;
691 : }
692 : }
693 : } // finished comparing additional vectors
694 :
695 :
696 : bool overall_result;
697 :
698 : // sum up the results
699 0 : if ((name_result==0) && (solu_result==-1))
700 : {
701 0 : if (ov_result.size()==0)
702 0 : overall_result = true;
703 : else
704 : {
705 : bool ov_identical;
706 0 : unsigned int n = 0;
707 0 : do
708 : {
709 0 : ov_identical = (ov_result[n]==-1);
710 0 : n++;
711 : }
712 0 : while (ov_identical && n<ov_result.size());
713 0 : overall_result = ov_identical;
714 0 : }
715 : }
716 : else
717 0 : overall_result = false;
718 :
719 0 : if (verbose)
720 : {
721 0 : libMesh::out << " finished comparisons, ";
722 0 : if (overall_result)
723 0 : libMesh::out << "found no differences." << std::endl << std::endl;
724 : else
725 0 : libMesh::out << "found differences." << std::endl << std::endl;
726 : }
727 :
728 0 : return overall_result;
729 : }
730 :
731 :
732 :
733 71 : void System::update_global_solution (std::vector<Number> & global_soln) const
734 : {
735 2 : parallel_object_only();
736 :
737 71 : global_soln.resize (solution->size());
738 :
739 71 : solution->localize (global_soln);
740 71 : }
741 :
742 :
743 :
744 5578 : void System::update_global_solution (std::vector<Number> & global_soln,
745 : const processor_id_type dest_proc) const
746 : {
747 232 : parallel_object_only();
748 :
749 5578 : global_soln.resize (solution->size());
750 :
751 5578 : solution->localize_to_one (global_soln, dest_proc);
752 5578 : }
753 :
754 :
755 :
756 210862 : NumericVector<Number> & System::add_vector (std::string_view vec_name,
757 : const bool projections,
758 : const ParallelType type)
759 : {
760 6266 : parallel_object_only();
761 :
762 6266 : libmesh_assert(this->comm().verify(std::string(vec_name)));
763 6266 : libmesh_assert(this->comm().verify(int(type)));
764 6266 : libmesh_assert(this->comm().verify(projections));
765 :
766 : // Return the vector if it is already there.
767 210862 : if (auto it = this->_vectors.find(vec_name);
768 6266 : it != this->_vectors.end())
769 : {
770 : // If the projection setting has *upgraded*, change it.
771 119165 : if (projections) // only do expensive lookup if needed
772 46734 : libmesh_map_find(_vector_projections, vec_name) = projections;
773 :
774 3550 : NumericVector<Number> & vec = *it->second;
775 :
776 : // If we're in serial, our vectors are effectively SERIAL, so
777 : // we'll ignore any type setting. If we're in parallel, we
778 : // might have a type change to deal with.
779 :
780 122715 : if (this->n_processors() > 1)
781 : {
782 : // If the type setting has changed in a way we can't
783 : // perceive as an upgrade or a downgrade, scream.
784 3550 : libmesh_assert_equal_to(type == SERIAL,
785 : vec.type() == SERIAL);
786 :
787 : // If the type setting has *upgraded*, change it.
788 117339 : if (type == GHOSTED && vec.type() == PARALLEL)
789 : {
790 : // A *really* late upgrade is expensive, but better not
791 : // to risk zeroing data.
792 138 : if (vec.initialized())
793 : {
794 69 : if (!vec.closed())
795 0 : vec.close();
796 :
797 : // Ideally we'd move parallel coefficients and then
798 : // add ghosted coefficients, but copy and swap is
799 : // simpler. If anyone actually ever uses this case
800 : // for real we can look into optimizing it.
801 71 : auto new_vec = NumericVector<Number>::build(this->comm());
802 : #ifdef LIBMESH_ENABLE_GHOSTED
803 71 : new_vec->init (this->n_dofs(), this->n_local_dofs(),
804 : _dof_map->get_send_list(), /*fast=*/false,
805 4 : GHOSTED);
806 : #else
807 : libmesh_error_msg("Cannot initialize ghosted vectors when they are not enabled.");
808 : #endif
809 :
810 69 : *new_vec = vec;
811 71 : vec.swap(*new_vec);
812 65 : }
813 : else
814 : // The PARALLEL vec is not yet initialized, so we can
815 : // just "upgrade" it to GHOSTED.
816 69 : vec.set_type(type);
817 : }
818 : }
819 :
820 : // Any upgrades are done; we're happy here.
821 3550 : return vec;
822 : }
823 :
824 : // Otherwise, build the vector. The following emplace() is
825 : // guaranteed to succeed because, if we made it here, we don't
826 : // already have a vector named "vec_name". We pass the user's
827 : // requested ParallelType directly to NumericVector::build() so
828 : // that, even if the vector is not initialized now, it will get the
829 : // right type when it is initialized later.
830 : auto pr =
831 : _vectors.emplace(vec_name,
832 183394 : NumericVector<Number>::build(this->comm(),
833 : libMesh::default_solver_package(),
834 5432 : type));
835 2716 : auto buf = pr.first->second.get();
836 2716 : _vector_projections.emplace(vec_name, projections);
837 :
838 : // Vectors are primal by default
839 91697 : _vector_is_adjoint.emplace(vec_name, -1);
840 :
841 : // Initialize it if necessary
842 91697 : if (_is_initialized)
843 : {
844 45702 : if (type == GHOSTED)
845 : {
846 : #ifdef LIBMESH_ENABLE_GHOSTED
847 4200 : buf->init (this->n_dofs(), this->n_local_dofs(),
848 : _dof_map->get_send_list(), /*fast=*/false,
849 240 : GHOSTED);
850 : #else
851 : libmesh_error_msg("Cannot initialize ghosted vectors when they are not enabled.");
852 : #endif
853 : }
854 : else
855 41502 : buf->init (this->n_dofs(), this->n_local_dofs(), false, type);
856 : }
857 :
858 2716 : return *buf;
859 : }
860 :
861 39413 : void System::remove_vector (std::string_view vec_name)
862 : {
863 1126 : parallel_object_only(); // Not strictly needed, but the only safe way to keep in sync
864 :
865 39413 : if (const auto pos = _vectors.find(vec_name);
866 1126 : pos != _vectors.end())
867 : {
868 38287 : _vectors.erase(pos);
869 1126 : auto proj_it = _vector_projections.find(vec_name);
870 1126 : libmesh_assert(proj_it != _vector_projections.end());
871 38287 : _vector_projections.erase(proj_it);
872 :
873 1126 : auto adj_it = _vector_is_adjoint.find(vec_name);
874 1126 : libmesh_assert(adj_it != _vector_is_adjoint.end());
875 38287 : _vector_is_adjoint.erase(adj_it);
876 : }
877 39413 : }
878 :
879 0 : const NumericVector<Number> * System::request_vector (std::string_view vec_name) const
880 : {
881 0 : if (const auto pos = _vectors.find(vec_name);
882 0 : pos != _vectors.end())
883 0 : return pos->second.get();
884 :
885 : // Otherwise, vec_name was not found
886 0 : return nullptr;
887 : }
888 :
889 :
890 :
891 0 : NumericVector<Number> * System::request_vector (std::string_view vec_name)
892 : {
893 0 : if (auto pos = _vectors.find(vec_name);
894 0 : pos != _vectors.end())
895 0 : return pos->second.get();
896 :
897 : // Otherwise, vec_name was not found
898 0 : return nullptr;
899 : }
900 :
901 :
902 :
903 0 : const NumericVector<Number> * System::request_vector (const unsigned int vec_num) const
904 : {
905 : // If we don't have that many vectors, return nullptr
906 0 : if (vec_num >= _vectors.size())
907 0 : return nullptr;
908 :
909 : // Otherwise return a pointer to the vec_num'th vector
910 0 : auto it = vectors_begin();
911 0 : std::advance(it, vec_num);
912 0 : return it->second.get();
913 : }
914 :
915 :
916 :
917 0 : NumericVector<Number> * System::request_vector (const unsigned int vec_num)
918 : {
919 : // If we don't have that many vectors, return nullptr
920 0 : if (vec_num >= _vectors.size())
921 0 : return nullptr;
922 :
923 : // Otherwise return a pointer to the vec_num'th vector
924 0 : auto it = vectors_begin();
925 0 : std::advance(it, vec_num);
926 0 : return it->second.get();
927 : }
928 :
929 :
930 :
931 5792 : const NumericVector<Number> & System::get_vector (std::string_view vec_name) const
932 : {
933 5792 : return *(libmesh_map_find(_vectors, vec_name));
934 : }
935 :
936 :
937 :
938 6382872 : NumericVector<Number> & System::get_vector (std::string_view vec_name)
939 : {
940 6382872 : return *(libmesh_map_find(_vectors, vec_name));
941 : }
942 :
943 :
944 :
945 0 : const NumericVector<Number> & System::get_vector (const unsigned int vec_num) const
946 : {
947 : // If we don't have that many vectors, throw an error
948 0 : libmesh_assert_less(vec_num, _vectors.size());
949 :
950 : // Otherwise return a reference to the vec_num'th vector
951 0 : auto it = vectors_begin();
952 0 : std::advance(it, vec_num);
953 0 : return *(it->second);
954 : }
955 :
956 :
957 :
958 0 : NumericVector<Number> & System::get_vector (const unsigned int vec_num)
959 : {
960 : // If we don't have that many vectors, throw an error
961 0 : libmesh_assert_less(vec_num, _vectors.size());
962 :
963 : // Otherwise return a reference to the vec_num'th vector
964 0 : auto it = vectors_begin();
965 0 : std::advance(it, vec_num);
966 0 : return *(it->second);
967 : }
968 :
969 :
970 :
971 0 : const std::string & System::vector_name (const unsigned int vec_num) const
972 : {
973 : // If we don't have that many vectors, throw an error
974 0 : libmesh_assert_less(vec_num, _vectors.size());
975 :
976 : // Otherwise return a reference to the vec_num'th vector name
977 0 : auto it = vectors_begin();
978 0 : std::advance(it, vec_num);
979 0 : return it->first;
980 : }
981 :
982 66238 : const std::string & System::vector_name (const NumericVector<Number> & vec_reference) const
983 : {
984 : // Linear search for a vector whose pointer matches vec_reference
985 62454 : auto it = std::find_if(vectors_begin(), vectors_end(),
986 31852 : [&vec_reference](const decltype(_vectors)::value_type & pr)
987 17818 : { return &vec_reference == pr.second.get(); });
988 :
989 : // Before returning, make sure we didn't loop till the end and not find any match
990 1892 : libmesh_assert (it != vectors_end());
991 :
992 : // Return the string associated with the current vector
993 66238 : return it->first;
994 : }
995 :
996 :
997 :
998 20374 : SparseMatrix<Number> & System::add_matrix (std::string_view mat_name,
999 : const ParallelType type,
1000 : const MatrixBuildType mat_build_type)
1001 : {
1002 566 : parallel_object_only();
1003 :
1004 566 : libmesh_assert(this->comm().verify(std::string(mat_name)));
1005 566 : libmesh_assert(this->comm().verify(int(type)));
1006 566 : libmesh_assert(this->comm().verify(int(mat_build_type)));
1007 :
1008 : // Return the matrix if it is already there.
1009 20374 : if (auto it = this->_matrices.find(mat_name);
1010 566 : it != this->_matrices.end())
1011 0 : return *it->second;
1012 :
1013 : // Otherwise build the matrix to return.
1014 19808 : std::unique_ptr<SparseMatrix<Number>> matrix;
1015 20374 : if (this->has_static_condensation())
1016 : {
1017 0 : if (mat_build_type == MatrixBuildType::DIAGONAL)
1018 0 : libmesh_error_msg(
1019 : "We do not currently support static condensation of the diagonal matrix type");
1020 0 : matrix = std::make_unique<StaticCondensation>(this->get_mesh(),
1021 : *this,
1022 : this->get_dof_map(),
1023 0 : this->get_dof_map().get_static_condensation());
1024 : }
1025 : else
1026 40182 : matrix = SparseMatrix<Number>::build(this->comm(), libMesh::default_solver_package());
1027 566 : auto & mat = *matrix;
1028 :
1029 566 : _matrices.emplace(mat_name, std::move(matrix));
1030 :
1031 566 : _matrix_types.emplace(mat_name, type);
1032 :
1033 : // Initialize it first if we've already initialized the others.
1034 20374 : this->late_matrix_init(mat, type);
1035 :
1036 566 : return mat;
1037 19242 : }
1038 :
1039 :
1040 :
1041 630 : SparseMatrix<Number> & System::add_matrix (std::string_view mat_name,
1042 : std::unique_ptr<SparseMatrix<Number>> matrix,
1043 : const ParallelType type)
1044 : {
1045 18 : parallel_object_only();
1046 :
1047 612 : const std::string namestr{mat_name};
1048 :
1049 18 : libmesh_assert(this->comm().verify(namestr));
1050 18 : libmesh_assert(this->comm().verify(int(type)));
1051 :
1052 18 : SparseMatrix<Number> & mat = *matrix;
1053 :
1054 630 : _matrices[namestr] = std::move(matrix);
1055 630 : _matrix_types[namestr] = type;
1056 :
1057 : // Initialize it first if we've already initialized the others.
1058 630 : this->late_matrix_init(mat, type);
1059 :
1060 648 : return mat;
1061 : }
1062 :
1063 21004 : void System::late_matrix_init(SparseMatrix<Number> & mat,
1064 : ParallelType type)
1065 : {
1066 21004 : if (_matrices_initialized)
1067 : {
1068 0 : this->get_dof_map().attach_matrix(mat);
1069 0 : mat.init(type);
1070 : }
1071 21004 : }
1072 :
1073 :
1074 :
1075 :
1076 0 : void System::remove_matrix (std::string_view mat_name)
1077 : {
1078 0 : parallel_object_only(); // Not strictly needed, but the only safe way to keep in sync
1079 :
1080 0 : if (const auto pos = _matrices.find(mat_name);
1081 0 : pos != _matrices.end())
1082 0 : _matrices.erase(pos); // erase()'d entries are destroyed
1083 0 : }
1084 :
1085 :
1086 :
1087 14 : const SparseMatrix<Number> * System::request_matrix (std::string_view mat_name) const
1088 : {
1089 14 : if (const auto pos = _matrices.find(mat_name);
1090 14 : pos != _matrices.end())
1091 14 : return pos->second.get();
1092 :
1093 : // Otherwise, mat_name does not exist
1094 0 : return nullptr;
1095 : }
1096 :
1097 :
1098 :
1099 257233 : SparseMatrix<Number> * System::request_matrix (std::string_view mat_name)
1100 : {
1101 257233 : if (auto pos = _matrices.find(mat_name);
1102 7508 : pos != _matrices.end())
1103 2 : return pos->second.get();
1104 :
1105 : // Otherwise, mat_name does not exist
1106 7506 : return nullptr;
1107 : }
1108 :
1109 :
1110 :
1111 0 : const SparseMatrix<Number> & System::get_matrix (std::string_view mat_name) const
1112 : {
1113 0 : return *libmesh_map_find(_matrices, mat_name);
1114 : }
1115 :
1116 :
1117 :
1118 665417 : SparseMatrix<Number> & System::get_matrix (std::string_view mat_name)
1119 : {
1120 665417 : return *libmesh_map_find(_matrices, mat_name);
1121 : }
1122 :
1123 :
1124 :
1125 69038 : void System::set_vector_preservation (const std::string & vec_name,
1126 : bool preserve)
1127 : {
1128 1972 : parallel_object_only(); // Not strictly needed, but the only safe way to keep in sync
1129 :
1130 69038 : _vector_projections[vec_name] = preserve;
1131 69038 : }
1132 :
1133 :
1134 :
1135 159257 : bool System::vector_preservation (std::string_view vec_name) const
1136 : {
1137 159257 : if (auto it = _vector_projections.find(vec_name);
1138 4550 : it != _vector_projections.end())
1139 159257 : return it->second;
1140 :
1141 : // vec_name was not in the map, return false
1142 0 : return false;
1143 : }
1144 :
1145 :
1146 :
1147 24746 : void System::set_vector_as_adjoint (const std::string & vec_name,
1148 : int qoi_num)
1149 : {
1150 798 : parallel_object_only(); // Not strictly needed, but the only safe way to keep in sync
1151 :
1152 : // We reserve -1 for vectors which get primal constraints, -2 for
1153 : // vectors which get no constraints
1154 798 : libmesh_assert_greater_equal(qoi_num, -2);
1155 24746 : _vector_is_adjoint[vec_name] = qoi_num;
1156 24746 : }
1157 :
1158 :
1159 :
1160 27375 : int System::vector_is_adjoint (std::string_view vec_name) const
1161 : {
1162 1020 : const auto it = _vector_is_adjoint.find(vec_name);
1163 1020 : libmesh_assert(it != _vector_is_adjoint.end());
1164 27375 : return it->second;
1165 : }
1166 :
1167 :
1168 :
1169 142 : NumericVector<Number> & System::add_sensitivity_solution (unsigned int i)
1170 : {
1171 142 : std::ostringstream sensitivity_name;
1172 138 : sensitivity_name << "sensitivity_solution" << i;
1173 :
1174 284 : return this->add_vector(sensitivity_name.str());
1175 134 : }
1176 :
1177 :
1178 :
1179 284 : NumericVector<Number> & System::get_sensitivity_solution (unsigned int i)
1180 : {
1181 284 : std::ostringstream sensitivity_name;
1182 276 : sensitivity_name << "sensitivity_solution" << i;
1183 :
1184 568 : return this->get_vector(sensitivity_name.str());
1185 268 : }
1186 :
1187 :
1188 :
1189 0 : const NumericVector<Number> & System::get_sensitivity_solution (unsigned int i) const
1190 : {
1191 0 : std::ostringstream sensitivity_name;
1192 0 : sensitivity_name << "sensitivity_solution" << i;
1193 :
1194 0 : return this->get_vector(sensitivity_name.str());
1195 0 : }
1196 :
1197 :
1198 :
1199 0 : NumericVector<Number> & System::add_weighted_sensitivity_solution ()
1200 : {
1201 0 : return this->add_vector("weighted_sensitivity_solution");
1202 : }
1203 :
1204 :
1205 :
1206 0 : NumericVector<Number> & System::get_weighted_sensitivity_solution ()
1207 : {
1208 0 : return this->get_vector("weighted_sensitivity_solution");
1209 : }
1210 :
1211 :
1212 :
1213 0 : const NumericVector<Number> & System::get_weighted_sensitivity_solution () const
1214 : {
1215 0 : return this->get_vector("weighted_sensitivity_solution");
1216 : }
1217 :
1218 :
1219 :
1220 24746 : NumericVector<Number> & System::add_adjoint_solution (unsigned int i)
1221 : {
1222 24746 : std::ostringstream adjoint_name;
1223 23948 : adjoint_name << "adjoint_solution" << i;
1224 :
1225 24746 : NumericVector<Number> & returnval = this->add_vector(adjoint_name.str());
1226 25544 : this->set_vector_as_adjoint(adjoint_name.str(), i);
1227 25544 : return returnval;
1228 23150 : }
1229 :
1230 :
1231 :
1232 2885957 : NumericVector<Number> & System::get_adjoint_solution (unsigned int i)
1233 : {
1234 2885957 : std::ostringstream adjoint_name;
1235 2663375 : adjoint_name << "adjoint_solution" << i;
1236 :
1237 5771914 : return this->get_vector(adjoint_name.str());
1238 2440605 : }
1239 :
1240 :
1241 :
1242 5792 : const NumericVector<Number> & System::get_adjoint_solution (unsigned int i) const
1243 : {
1244 5792 : std::ostringstream adjoint_name;
1245 5480 : adjoint_name << "adjoint_solution" << i;
1246 :
1247 11584 : return this->get_vector(adjoint_name.str());
1248 5168 : }
1249 :
1250 :
1251 :
1252 0 : NumericVector<Number> & System::add_weighted_sensitivity_adjoint_solution (unsigned int i)
1253 : {
1254 0 : std::ostringstream adjoint_name;
1255 0 : adjoint_name << "weighted_sensitivity_adjoint_solution" << i;
1256 :
1257 0 : NumericVector<Number> & returnval = this->add_vector(adjoint_name.str());
1258 0 : this->set_vector_as_adjoint(adjoint_name.str(), i);
1259 0 : return returnval;
1260 0 : }
1261 :
1262 :
1263 :
1264 0 : NumericVector<Number> & System::get_weighted_sensitivity_adjoint_solution (unsigned int i)
1265 : {
1266 0 : std::ostringstream adjoint_name;
1267 0 : adjoint_name << "weighted_sensitivity_adjoint_solution" << i;
1268 :
1269 0 : return this->get_vector(adjoint_name.str());
1270 0 : }
1271 :
1272 :
1273 :
1274 0 : const NumericVector<Number> & System::get_weighted_sensitivity_adjoint_solution (unsigned int i) const
1275 : {
1276 0 : std::ostringstream adjoint_name;
1277 0 : adjoint_name << "weighted_sensitivity_adjoint_solution" << i;
1278 :
1279 0 : return this->get_vector(adjoint_name.str());
1280 0 : }
1281 :
1282 :
1283 :
1284 24817 : NumericVector<Number> & System::add_adjoint_rhs (unsigned int i)
1285 : {
1286 24817 : std::ostringstream adjoint_rhs_name;
1287 24017 : adjoint_rhs_name << "adjoint_rhs" << i;
1288 :
1289 49634 : return this->add_vector(adjoint_rhs_name.str(), false);
1290 23217 : }
1291 :
1292 :
1293 :
1294 2563081 : NumericVector<Number> & System::get_adjoint_rhs (unsigned int i)
1295 : {
1296 2563081 : std::ostringstream adjoint_rhs_name;
1297 2332959 : adjoint_rhs_name << "adjoint_rhs" << i;
1298 :
1299 5126162 : return this->get_vector(adjoint_rhs_name.str());
1300 2102825 : }
1301 :
1302 :
1303 :
1304 0 : const NumericVector<Number> & System::get_adjoint_rhs (unsigned int i) const
1305 : {
1306 0 : std::ostringstream adjoint_rhs_name;
1307 0 : adjoint_rhs_name << "adjoint_rhs" << i;
1308 :
1309 0 : return this->get_vector(adjoint_rhs_name.str());
1310 0 : }
1311 :
1312 :
1313 :
1314 7430 : NumericVector<Number> & System::add_sensitivity_rhs (unsigned int i)
1315 : {
1316 7430 : std::ostringstream sensitivity_rhs_name;
1317 7218 : sensitivity_rhs_name << "sensitivity_rhs" << i;
1318 :
1319 14860 : return this->add_vector(sensitivity_rhs_name.str(), false);
1320 7006 : }
1321 :
1322 :
1323 :
1324 7430 : NumericVector<Number> & System::get_sensitivity_rhs (unsigned int i)
1325 : {
1326 7430 : std::ostringstream sensitivity_rhs_name;
1327 7218 : sensitivity_rhs_name << "sensitivity_rhs" << i;
1328 :
1329 14860 : return this->get_vector(sensitivity_rhs_name.str());
1330 7006 : }
1331 :
1332 :
1333 :
1334 0 : const NumericVector<Number> & System::get_sensitivity_rhs (unsigned int i) const
1335 : {
1336 0 : std::ostringstream sensitivity_rhs_name;
1337 0 : sensitivity_rhs_name << "sensitivity_rhs" << i;
1338 :
1339 0 : return this->get_vector(sensitivity_rhs_name.str());
1340 0 : }
1341 :
1342 :
1343 :
1344 273574 : unsigned int System::add_variable (std::string_view var,
1345 : const FEType & type,
1346 : const std::set<subdomain_id_type> * const active_subdomains)
1347 : {
1348 273574 : return this->get_dof_map().add_variable(*this, var, type, active_subdomains);
1349 : }
1350 :
1351 :
1352 :
1353 268242 : unsigned int System::add_variable (std::string_view var,
1354 : const Order order,
1355 : const FEFamily family,
1356 : const std::set<subdomain_id_type> * const active_subdomains,
1357 : const bool p_refinement)
1358 : {
1359 268242 : return this->add_variable(var,
1360 275830 : FEType(order, family).set_p_refinement(p_refinement),
1361 275830 : active_subdomains);
1362 : }
1363 :
1364 :
1365 :
1366 71 : unsigned int System::add_variables (const std::vector<std::string> & vars,
1367 : const FEType & type,
1368 : const std::set<subdomain_id_type> * const active_subdomains)
1369 : {
1370 71 : return this->get_dof_map().add_variables(*this, vars, type, active_subdomains);
1371 : }
1372 :
1373 :
1374 :
1375 71 : unsigned int System::add_variables (const std::vector<std::string> & vars,
1376 : const Order order,
1377 : const FEFamily family,
1378 : const std::set<subdomain_id_type> * const active_subdomains,
1379 : const bool p_refinement)
1380 : {
1381 71 : return this->add_variables(vars,
1382 73 : FEType(order, family).set_p_refinement(p_refinement),
1383 73 : active_subdomains);
1384 : }
1385 :
1386 568 : unsigned int System::add_variable_array (const std::vector<std::string> & vars,
1387 : const FEType & type,
1388 : const std::set<subdomain_id_type> * const active_subdomains)
1389 : {
1390 568 : return this->get_dof_map().add_variable_array(*this, vars, type, active_subdomains);
1391 : }
1392 :
1393 894 : bool System::has_variable (std::string_view var) const
1394 : {
1395 894 : return this->get_dof_map().has_variable(var);
1396 : }
1397 :
1398 8903978 : unsigned int System::variable_number (std::string_view var) const
1399 : {
1400 8903978 : return this->get_dof_map().variable_number(var);
1401 : }
1402 :
1403 493 : void System::get_all_variable_numbers(std::vector<unsigned int> & all_variable_numbers) const
1404 : {
1405 493 : this->get_dof_map().get_all_variable_numbers(all_variable_numbers);
1406 493 : }
1407 :
1408 :
1409 142 : void System::local_dof_indices(const unsigned int var,
1410 : std::set<dof_id_type> & var_indices) const
1411 : {
1412 : // Make sure the set is clear
1413 4 : var_indices.clear();
1414 :
1415 8 : std::vector<dof_id_type> dof_indices;
1416 :
1417 : const dof_id_type
1418 4 : first_local = this->get_dof_map().first_dof(),
1419 4 : end_local = this->get_dof_map().end_dof();
1420 :
1421 : // Begin the loop over the elements
1422 984 : for (const auto & elem : this->get_mesh().active_local_element_ptr_range())
1423 : {
1424 384 : this->get_dof_map().dof_indices (elem, dof_indices, var);
1425 :
1426 1152 : for (dof_id_type dof : dof_indices)
1427 : //If the dof is owned by the local processor
1428 768 : if (first_local <= dof && dof < end_local)
1429 503 : var_indices.insert(dof);
1430 134 : }
1431 :
1432 : // we may have missed assigning DOFs to nodes that we own
1433 : // but to which we have no connected elements matching our
1434 : // variable restriction criterion. this will happen, for example,
1435 : // if variable V is restricted to subdomain S. We may not own
1436 : // any elements which live in S, but we may own nodes which are
1437 : // *connected* to elements which do.
1438 1380 : for (const auto & node : this->get_mesh().local_node_ptr_range())
1439 : {
1440 50 : libmesh_assert(node);
1441 600 : this->get_dof_map().dof_indices (node, dof_indices, var);
1442 960 : for (auto dof : dof_indices)
1443 360 : if (first_local <= dof && dof < end_local)
1444 330 : var_indices.insert(dof);
1445 134 : }
1446 142 : }
1447 :
1448 :
1449 :
1450 0 : void System::zero_variable (NumericVector<Number> & v,
1451 : unsigned int var_num) const
1452 : {
1453 : /* Make sure the call makes sense. */
1454 0 : libmesh_assert_less (var_num, this->n_vars());
1455 :
1456 : /* Get a reference to the mesh. */
1457 0 : const MeshBase & mesh = this->get_mesh();
1458 :
1459 : /* Check which system we are. */
1460 0 : const unsigned int sys_num = this->number();
1461 :
1462 : // Loop over nodes.
1463 0 : for (const auto & node : mesh.local_node_ptr_range())
1464 : {
1465 0 : unsigned int n_comp = node->n_comp(sys_num,var_num);
1466 0 : for (unsigned int i=0; i<n_comp; i++)
1467 : {
1468 0 : const dof_id_type index = node->dof_number(sys_num,var_num,i);
1469 0 : v.set(index,0.0);
1470 : }
1471 0 : }
1472 :
1473 : // Loop over elements.
1474 : Threads::parallel_for
1475 0 : (mesh.active_local_element_stored_range(),
1476 0 : [sys_num, var_num, &v](const ConstElemRange & range)
1477 : {
1478 0 : for (const Elem * elem : range)
1479 : {
1480 0 : unsigned int n_comp = elem->n_comp(sys_num,var_num);
1481 0 : for (unsigned int i=0; i<n_comp; i++)
1482 : {
1483 0 : const dof_id_type index = elem->dof_number(sys_num,var_num,i);
1484 0 : v.set(index,0.0);
1485 : }
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 8344 : 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 31907456 : for (const auto & elem : this->get_mesh().active_local_element_ptr_range())
1663 : {
1664 17489993 : 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 3519431 : if (elem->infinite() )
1670 0 : libmesh_not_implemented();
1671 :
1672 17493791 : if (skip_dimensions && skip_dimensions->find(dim) != skip_dimensions->end())
1673 25062 : continue;
1674 :
1675 19090199 : QBase * qrule = q_rules[dim].get();
1676 1682932 : libmesh_assert(qrule);
1677 :
1678 17464931 : this->get_dof_map().dof_indices (elem, dof_indices, var);
1679 :
1680 14156731 : auto element_calculation = [&dof_indices, &elem,
1681 : norm_type, norm_weight, norm_weight_sq, &qrule,
1682 953826418 : &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 17464931 : const std::vector<Real> & JxW = fe.get_JxW();
1689 1682932 : const std::vector<std::vector<OutputShape>> * phi = nullptr;
1690 17464931 : if (norm_type == H1 ||
1691 15784208 : norm_type == H2 ||
1692 15782827 : norm_type == L2 ||
1693 15782827 : norm_type == L1 ||
1694 : norm_type == L_INF)
1695 1682380 : phi = &(fe.get_phi());
1696 :
1697 1682932 : const std::vector<std::vector<OutputGradient>> * dphi = nullptr;
1698 19090199 : if (norm_type == H1 ||
1699 15784208 : norm_type == H2 ||
1700 15783932 : norm_type == H1_SEMINORM ||
1701 : norm_type == W1_INF_SEMINORM)
1702 1681275 : dphi = &(fe.get_dphi());
1703 :
1704 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1705 : typedef typename std::remove_reference<decltype(fe)>::type::OutputTensor OutputTensor;
1706 :
1707 1682932 : const std::vector<std::vector<OutputTensor>> * d2phi = nullptr;
1708 19090199 : if (norm_type == H2 ||
1709 17464931 : norm_type == H2_SEMINORM ||
1710 : norm_type == W2_INF_SEMINORM)
1711 0 : d2phi = &(fe.get_d2phi());
1712 : #endif
1713 :
1714 17464931 : fe.reinit (elem);
1715 :
1716 17464931 : const unsigned int n_qp = qrule->n_points();
1717 :
1718 : const unsigned int n_sf = cast_int<unsigned int>
1719 3308200 : (dof_indices.size());
1720 :
1721 : // Begin the loop over the Quadrature points.
1722 91185407 : for (unsigned int qp=0; qp<n_qp; qp++)
1723 : {
1724 73720476 : 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 73720476 : 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 73720476 : if (norm_type == H1 ||
1742 66685908 : norm_type == H2 ||
1743 : norm_type == L2)
1744 : {
1745 7051411 : OutputNumberShape u_h = 0.;
1746 420041952 : for (unsigned int i=0; i != n_sf; ++i)
1747 471846076 : u_h += (*phi)[i][qp] * (*local_v)(dof_indices[i]);
1748 80439549 : v_norm += norm_weight_sq *
1749 73618794 : JxW[qp] * TensorTools::norm_sq(u_h);
1750 : }
1751 :
1752 80549547 : if (norm_type == H1 ||
1753 66685908 : norm_type == H2 ||
1754 66660749 : norm_type == H1_SEMINORM)
1755 : {
1756 7037340 : OutputNumberGradient grad_u_h;
1757 413874756 : for (unsigned int i=0; i != n_sf; ++i)
1758 371655578 : grad_u_h.add_scaled((*dphi)[i][qp], (*local_v)(dof_indices[i]));
1759 80256626 : v_norm += norm_weight_sq *
1760 73449942 : JxW[qp] * grad_u_h.norm_sq();
1761 : }
1762 :
1763 73720476 : 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 80549547 : if (norm_type == H2 ||
1775 66660749 : 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 73720476 : 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 34929862 : };
1794 :
1795 3308200 : FEBase * scalar_fe = fe_ptrs[dim].get();
1796 3308200 : FEVectorBase * vec_fe = vec_fe_ptrs[dim].get();
1797 :
1798 17464931 : if (scalar_fe)
1799 : {
1800 1681827 : libmesh_assert(!vec_fe);
1801 17451671 : element_calculation(*scalar_fe);
1802 : }
1803 :
1804 17464931 : 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 17958 : std::string System::get_info() const
1828 : {
1829 18954 : std::ostringstream oss;
1830 :
1831 :
1832 498 : const std::string & sys_name = this->name();
1833 :
1834 17958 : oss << " System #" << this->number() << ", \"" << sys_name << "\"\n"
1835 17958 : << " Type \"" << this->system_type() << "\"\n"
1836 51882 : << " Variables=";
1837 :
1838 58240 : for (auto vg : make_range(this->n_variable_groups()))
1839 : {
1840 23320 : const VariableGroup & vg_description (this->variable_group(vg));
1841 :
1842 23320 : if (vg_description.n_variables() > 1) oss << "{ ";
1843 57743 : for (auto vn : make_range(vg_description.n_variables()))
1844 66918 : oss << "\"" << vg_description.name(vn) << "\" ";
1845 23320 : if (vg_description.n_variables() > 1) oss << "} ";
1846 : }
1847 :
1848 17958 : oss << '\n';
1849 :
1850 17958 : oss << " Finite Element Types=";
1851 : #ifndef LIBMESH_ENABLE_INFINITE_ELEMENTS
1852 55335 : for (auto vg : make_range(this->n_variable_groups()))
1853 : oss << "\""
1854 43578 : << Utility::enum_to_string<FEFamily>(this->get_dof_map().variable_group(vg).type().family)
1855 43578 : << "\" ";
1856 : #else
1857 2905 : for (auto vg : make_range(this->n_variable_groups()))
1858 : {
1859 : oss << "\""
1860 3062 : << Utility::enum_to_string<FEFamily>(this->get_dof_map().variable_group(vg).type().family)
1861 : << "\", \""
1862 3062 : << Utility::enum_to_string<FEFamily>(this->get_dof_map().variable_group(vg).type().radial_family)
1863 3293 : << "\" ";
1864 : }
1865 :
1866 1185 : oss << '\n' << " Infinite Element Mapping=";
1867 2905 : for (auto vg : make_range(this->n_variable_groups()))
1868 : oss << "\""
1869 3062 : << Utility::enum_to_string<InfMapType>(this->get_dof_map().variable_group(vg).type().inf_map)
1870 2412 : << "\" ";
1871 : #endif
1872 :
1873 17958 : oss << '\n';
1874 :
1875 17958 : oss << " Approximation Orders=";
1876 58240 : for (auto vg : make_range(this->n_variable_groups()))
1877 : {
1878 : #ifndef LIBMESH_ENABLE_INFINITE_ELEMENTS
1879 : oss << "\""
1880 43578 : << Utility::enum_to_string<Order>(this->get_dof_map().variable_group(vg).type().order)
1881 43578 : << "\" ";
1882 : #else
1883 : oss << "\""
1884 3062 : << Utility::enum_to_string<Order>(this->get_dof_map().variable_group(vg).type().order)
1885 : << "\", \""
1886 3062 : << Utility::enum_to_string<Order>(this->get_dof_map().variable_group(vg).type().radial_order)
1887 3293 : << "\" ";
1888 : #endif
1889 : }
1890 :
1891 17958 : oss << '\n';
1892 :
1893 17958 : if (this->is_initialized())
1894 : {
1895 34782 : oss << " n_dofs()=" << this->n_dofs() << '\n';
1896 17887 : dof_id_type local_dofs = this->n_local_dofs();
1897 17887 : oss << " n_local_dofs()=" << local_dofs << '\n';
1898 17887 : this->comm().max(local_dofs);
1899 17887 : oss << " max(n_local_dofs())=" << local_dofs << '\n';
1900 : #ifdef LIBMESH_ENABLE_CONSTRAINTS
1901 17887 : if (this->n_constrained_dofs())
1902 : {
1903 19816 : oss << " n_constrained_dofs()=" << this->n_constrained_dofs() << '\n';
1904 19530 : oss << " n_local_constrained_dofs()=" << this->n_local_constrained_dofs() << '\n';
1905 10051 : dof_id_type local_unconstrained_dofs = this->n_local_dofs() - this->n_local_constrained_dofs();
1906 10051 : this->comm().max(local_unconstrained_dofs);
1907 10337 : oss << " max(local unconstrained dofs)=" << local_unconstrained_dofs << '\n';
1908 : }
1909 : #endif
1910 17887 : 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 34920 : oss << " " << "n_vectors()=" << this->n_vectors() << '\n';
1918 34920 : oss << " " << "n_matrices()=" << this->n_matrices() << '\n';
1919 : // oss << " " << "n_additional_matrices()=" << this->n_additional_matrices() << '\n';
1920 :
1921 35418 : oss << this->get_dof_map().get_info();
1922 :
1923 18456 : return oss.str();
1924 16962 : }
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 9986 : void System::attach_assemble_function (void fptr(EquationSystems & es,
1960 : const std::string & name))
1961 : {
1962 278 : libmesh_assert(fptr);
1963 :
1964 9986 : 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 9986 : _assemble_system_function = fptr;
1972 9986 : }
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 2485 : void System::attach_constraint_object (System::Constraint & constrain)
2008 : {
2009 2485 : 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 2485 : _constrain_system_object = &constrain;
2017 2485 : }
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 241875 : void System::user_initialization ()
2096 : {
2097 : // Call the user-provided initialization function,
2098 : // if it was provided
2099 241875 : if (_init_system_function != nullptr)
2100 501 : this->_init_system_function (_equation_systems, this->name());
2101 :
2102 : // ...or the user-provided initialization object.
2103 241374 : else if (_init_system_object != nullptr)
2104 280 : this->_init_system_object->initialize();
2105 241875 : }
2106 :
2107 :
2108 :
2109 46331 : void System::user_assembly ()
2110 : {
2111 : // Call the user-provided assembly function,
2112 : // if it was provided
2113 46331 : if (_assemble_system_function != nullptr)
2114 46189 : 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 46331 : }
2120 :
2121 :
2122 :
2123 274280 : void System::user_constrain ()
2124 : {
2125 : // Call the user-provided constraint function,
2126 : // if it was provided
2127 274280 : if (_constrain_system_function!= nullptr)
2128 0 : this->_constrain_system_function(_equation_systems, this->name());
2129 :
2130 : // ...or the user-provided constraint object.
2131 274280 : else if (_constrain_system_object != nullptr)
2132 2485 : this->_constrain_system_object->constrain();
2133 274280 : }
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 339077 : 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 9516 : parallel_object_only();
2227 :
2228 : // And every processor had better agree about which point we're
2229 : // looking for
2230 : #ifndef NDEBUG
2231 9516 : libmesh_assert(this->comm().verify(p(0)));
2232 : #if LIBMESH_DIM > 1
2233 9516 : libmesh_assert(this->comm().verify(p(1)));
2234 : #endif
2235 : #if LIBMESH_DIM > 2
2236 9516 : 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 19032 : const MeshBase & mesh = this->get_mesh();
2242 :
2243 : // Use an existing PointLocator or create a new one
2244 339077 : std::unique_ptr<PointLocatorBase> locator_ptr = mesh.sub_point_locator();
2245 9516 : PointLocatorBase & locator = *locator_ptr;
2246 :
2247 339077 : if (!insist_on_success || !mesh.is_serial())
2248 293498 : 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 339077 : this->variable(var).active_subdomains();
2254 : const std::set<subdomain_id_type> * implicit_subdomains =
2255 339077 : raw_subdomains.empty() ? nullptr : &raw_subdomains;
2256 339077 : const Elem * e = locator(p, implicit_subdomains);
2257 :
2258 339077 : Number u = 0;
2259 :
2260 339077 : if (e && this->get_dof_map().is_evaluable(*e, var))
2261 120171 : 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 339077 : (e && (e->processor_id() == this->processor_id())) ?
2266 334319 : this->processor_id() : this->n_processors();
2267 339077 : 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 348593 : if (lowest_owner != this->n_processors())
2273 329561 : this->comm().broadcast(u, lowest_owner);
2274 : else
2275 0 : libmesh_assert(!insist_on_success);
2276 :
2277 358109 : return u;
2278 320045 : }
2279 :
2280 138552 : 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 7642 : libmesh_assert (e.contains_point(p));
2289 :
2290 138552 : if (!sol)
2291 4980 : sol = this->current_local_solution.get();
2292 :
2293 : // Get the dof map to get the proper indices for our computation
2294 7642 : const DofMap & dof_map = this->get_dof_map();
2295 :
2296 : // Make sure we can evaluate on this element.
2297 7642 : libmesh_assert (dof_map.is_evaluable(e, var));
2298 :
2299 : // Need dof_indices for phi[i][j]
2300 15284 : std::vector<dof_id_type> dof_indices;
2301 :
2302 : // Fill in the dof_indices for our element
2303 138552 : 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 15284 : (dof_indices.size());
2308 :
2309 138552 : FEType fe_type = dof_map.variable_type(var);
2310 :
2311 : // Map the physical co-ordinates to the master co-ordinates
2312 138552 : 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 146194 : FEComputeData fe_data(this->get_equation_systems(), coor);
2317 138552 : FEInterface::compute_data(e.dim(), fe_type, &e, fe_data);
2318 :
2319 : // Get ready to accumulate a value
2320 7642 : Number u = 0;
2321 :
2322 2538928 : for (unsigned int l=0; l<num_dofs; l++)
2323 : {
2324 2610602 : u += fe_data.shape[l] * (*sol)(dof_indices[l]);
2325 : }
2326 :
2327 146194 : return u;
2328 123268 : }
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 87371 : bool System::has_static_condensation() const
2670 : {
2671 87371 : return this->get_dof_map().has_static_condensation();
2672 : }
2673 :
2674 125798812 : unsigned int System::n_vars() const
2675 : {
2676 125798812 : return this->get_dof_map().n_vars();
2677 : }
2678 :
2679 30610369 : const std::string & System::variable_name (const unsigned int i) const
2680 : {
2681 30610369 : 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 273913 : unsigned int System::n_variable_groups() const
2700 : {
2701 273913 : return this->get_dof_map().n_variable_groups();
2702 : }
2703 :
2704 51609426 : const Variable & System::variable (const unsigned int i) const
2705 : {
2706 51609426 : return this->get_dof_map().variable(i);
2707 : }
2708 :
2709 285704 : const VariableGroup & System::variable_group (const unsigned int vg) const
2710 : {
2711 285704 : return this->get_dof_map().variable_group(vg);
2712 : }
2713 :
2714 : unsigned int
2715 6651847 : System::variable_scalar_number (unsigned int var_num,
2716 : unsigned int component) const
2717 : {
2718 6651847 : return this->get_dof_map().variable_scalar_number(var_num, component);
2719 : }
2720 :
2721 61847204 : const FEType & System::variable_type (const unsigned int i) const
2722 : {
2723 61847204 : return this->get_dof_map().variable_type(i);
2724 : }
2725 :
2726 34524 : const FEType & System::variable_type (std::string_view var) const
2727 : {
2728 34524 : return this->get_dof_map().variable_type(var);
2729 : }
2730 :
2731 : } // namespace libMesh
|
