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