1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2020 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 // C++ includes
21 #include <sstream> // for std::ostringstream
22
23 // Local includes
24 #include "libmesh/dof_map.h"
25 #include "libmesh/equation_systems.h"
26 #include "libmesh/int_range.h"
27 #include "libmesh/libmesh_logging.h"
28 #include "libmesh/mesh_base.h"
29 #include "libmesh/numeric_vector.h"
30 #include "libmesh/parameter_vector.h"
31 #include "libmesh/point.h" // For point_value
32 #include "libmesh/point_locator_base.h" // For point_value
33 #include "libmesh/qoi_set.h"
34 #include "libmesh/enum_to_string.h"
35 #include "libmesh/system.h"
36 #include "libmesh/system_norm.h"
37 #include "libmesh/utility.h"
38 #include "libmesh/elem.h"
39 #include "libmesh/fe_type.h"
40 #include "libmesh/fe_interface.h"
41 #include "libmesh/fe_compute_data.h"
42 #include "libmesh/auto_ptr.h" // libmesh_make_unique
43
44 // includes for calculate_norm, point_*
45 #include "libmesh/fe_base.h"
46 #include "libmesh/fe_interface.h"
47 #include "libmesh/parallel.h"
48 #include "libmesh/parallel_algebra.h"
49 #include "libmesh/quadrature.h"
50 #include "libmesh/tensor_value.h"
51 #include "libmesh/vector_value.h"
52 #include "libmesh/tensor_tools.h"
53 #include "libmesh/enum_norm_type.h"
54
55 namespace libMesh
56 {
57
58
59 // ------------------------------------------------------------
60 // System implementation
System(EquationSystems & es,const std::string & name_in,const unsigned int number_in)61 System::System (EquationSystems & es,
62 const std::string & name_in,
63 const unsigned int number_in) :
64
65 ParallelObject (es),
66 assemble_before_solve (true),
67 use_fixed_solution (false),
68 extra_quadrature_order (0),
69 solution (NumericVector<Number>::build(this->comm())),
70 current_local_solution (NumericVector<Number>::build(this->comm())),
71 time (0.),
72 qoi (0),
73 _init_system_function (nullptr),
74 _init_system_object (nullptr),
75 _assemble_system_function (nullptr),
76 _assemble_system_object (nullptr),
77 _constrain_system_function (nullptr),
78 _constrain_system_object (nullptr),
79 _qoi_evaluate_function (nullptr),
80 _qoi_evaluate_object (nullptr),
81 _qoi_evaluate_derivative_function (nullptr),
82 _qoi_evaluate_derivative_object (nullptr),
83 _dof_map (libmesh_make_unique<DofMap>(number_in, es.get_mesh())),
84 _equation_systems (es),
85 _mesh (es.get_mesh()),
86 _sys_name (name_in),
87 _sys_number (number_in),
88 _active (true),
89 _solution_projection (true),
90 _basic_system_only (false),
91 _is_initialized (false),
92 _identify_variable_groups (true),
93 _additional_data_written (false),
94 adjoint_already_solved (false),
95 _hide_output (false),
96 project_with_constraints (true)
97 {
98 }
99
100
101
102 // No copy construction of System objects!
System(const System & other)103 System::System (const System & other) :
104 ReferenceCountedObject<System>(),
105 ParallelObject(other),
106 _equation_systems(other._equation_systems),
107 _mesh(other._mesh),
108 _sys_number(other._sys_number)
109 {
110 libmesh_not_implemented();
111 }
112
113
114
115 System & System::operator= (const System &)
116 {
117 libmesh_not_implemented();
118 }
119
120
~System()121 System::~System ()
122 {
123 // Null-out the function pointers. Since this
124 // class is getting destructed it is pointless,
125 // but a good habit.
126 _init_system_function =
127 _assemble_system_function =
128 _constrain_system_function = nullptr;
129
130 _qoi_evaluate_function = nullptr;
131 _qoi_evaluate_derivative_function = nullptr;
132
133 // nullptr-out user-provided objects.
134 _init_system_object = nullptr;
135 _assemble_system_object = nullptr;
136 _constrain_system_object = nullptr;
137 _qoi_evaluate_object = nullptr;
138 _qoi_evaluate_derivative_object = nullptr;
139
140 // Clear data
141 // Note: although clear() is virtual, C++ only calls
142 // the clear() of the base class in the destructor.
143 // Thus we add System namespace to make it clear.
144 System::clear ();
145
146 libmesh_exceptionless_assert (!libMesh::closed());
147 }
148
149
150
n_dofs()151 dof_id_type System::n_dofs() const
152 {
153 return _dof_map->n_dofs();
154 }
155
156
157
n_constrained_dofs()158 dof_id_type System::n_constrained_dofs() const
159 {
160 #ifdef LIBMESH_ENABLE_CONSTRAINTS
161
162 return _dof_map->n_constrained_dofs();
163
164 #else
165
166 return 0;
167
168 #endif
169 }
170
171
172
n_local_constrained_dofs()173 dof_id_type System::n_local_constrained_dofs() const
174 {
175 #ifdef LIBMESH_ENABLE_CONSTRAINTS
176
177 return _dof_map->n_local_constrained_dofs();
178
179 #else
180
181 return 0;
182
183 #endif
184 }
185
186
187
n_local_dofs()188 dof_id_type System::n_local_dofs() const
189 {
190 return _dof_map->n_dofs_on_processor (this->processor_id());
191 }
192
193
194
current_solution(const dof_id_type global_dof_number)195 Number System::current_solution (const dof_id_type global_dof_number) const
196 {
197 // Check the sizes
198 libmesh_assert_less (global_dof_number, _dof_map->n_dofs());
199 libmesh_assert_less (global_dof_number, current_local_solution->size());
200
201 return (*current_local_solution)(global_dof_number);
202 }
203
204
205
clear()206 void System::clear ()
207 {
208 _variables.clear();
209
210 _variable_numbers.clear();
211
212 _dof_map->clear ();
213
214 solution->clear ();
215
216 current_local_solution->clear ();
217
218 // clear any user-added vectors
219 {
220 for (auto & pr : _vectors)
221 {
222 pr.second->clear ();
223 delete pr.second;
224 pr.second = nullptr;
225 }
226
227 _vectors.clear();
228 _vector_projections.clear();
229 _vector_is_adjoint.clear();
230 _vector_types.clear();
231 _is_initialized = false;
232 }
233
234 }
235
236
237
init()238 void System::init ()
239 {
240 // Calling init() twice on the same system currently works evil
241 // magic, whether done directly or via EquationSystems::read()
242 libmesh_assert(!this->is_initialized());
243
244 // First initialize any required data:
245 // either only the basic System data
246 if (_basic_system_only)
247 System::init_data();
248 // or all the derived class' data too
249 else
250 this->init_data();
251
252 // If no variables have been added to this system
253 // don't do anything
254 if (!this->n_vars())
255 return;
256
257 // Then call the user-provided initialization function
258 this->user_initialization();
259 }
260
261
262
init_data()263 void System::init_data ()
264 {
265 MeshBase & mesh = this->get_mesh();
266
267 // Add all variable groups to our underlying DofMap
268 for (auto vg : make_range(this->n_variable_groups()))
269 _dof_map->add_variable_group(this->variable_group(vg));
270
271 // Distribute the degrees of freedom on the mesh
272 _dof_map->distribute_dofs (mesh);
273
274 // Recreate any user or internal constraints
275 this->reinit_constraints();
276
277 // And clean up the send_list before we first use it
278 _dof_map->prepare_send_list();
279
280 // Resize the solution conformal to the current mesh
281 solution->init (this->n_dofs(), this->n_local_dofs(), false, PARALLEL);
282
283 // Resize the current_local_solution for the current mesh
284 #ifdef LIBMESH_ENABLE_GHOSTED
285 current_local_solution->init (this->n_dofs(), this->n_local_dofs(),
286 _dof_map->get_send_list(), false,
287 GHOSTED);
288 #else
289 current_local_solution->init (this->n_dofs(), false, SERIAL);
290 #endif
291
292 // from now on, adding additional vectors or variables can't be done
293 // without immediately initializing them
294 _is_initialized = true;
295
296 // initialize & zero other vectors, if necessary
297 for (auto & pr : _vectors)
298 {
299 ParallelType type = _vector_types[pr.first];
300
301 if (type == GHOSTED)
302 {
303 #ifdef LIBMESH_ENABLE_GHOSTED
304 pr.second->init (this->n_dofs(), this->n_local_dofs(),
305 _dof_map->get_send_list(), false,
306 GHOSTED);
307 #else
308 libmesh_error_msg("Cannot initialize ghosted vectors when they are not enabled.");
309 #endif
310 }
311 else if (type == SERIAL)
312 {
313 pr.second->init (this->n_dofs(), false, type);
314 }
315 else
316 {
317 libmesh_assert_equal_to(type, PARALLEL);
318 pr.second->init (this->n_dofs(), this->n_local_dofs(), false, type);
319 }
320 }
321 }
322
323
324
restrict_vectors()325 void System::restrict_vectors ()
326 {
327 #ifdef LIBMESH_ENABLE_AMR
328 // Restrict the _vectors on the coarsened cells
329 for (auto & pr : _vectors)
330 {
331 NumericVector<Number> * v = pr.second;
332
333 if (_vector_projections[pr.first])
334 {
335 this->project_vector (*v, this->vector_is_adjoint(pr.first));
336 }
337 else
338 {
339 ParallelType type = _vector_types[pr.first];
340
341 if (type == GHOSTED)
342 {
343 #ifdef LIBMESH_ENABLE_GHOSTED
344 pr.second->init (this->n_dofs(), this->n_local_dofs(),
345 _dof_map->get_send_list(), false,
346 GHOSTED);
347 #else
348 libmesh_error_msg("Cannot initialize ghosted vectors when they are not enabled.");
349 #endif
350 }
351 else
352 pr.second->init (this->n_dofs(), this->n_local_dofs(), false, type);
353 }
354 }
355
356 const std::vector<dof_id_type> & send_list = _dof_map->get_send_list ();
357
358 // Restrict the solution on the coarsened cells
359 if (_solution_projection)
360 this->project_vector (*solution);
361 // Or at least make sure the solution vector is the correct size
362 else
363 solution->init (this->n_dofs(), this->n_local_dofs(), true, PARALLEL);
364
365 #ifdef LIBMESH_ENABLE_GHOSTED
366 current_local_solution->init(this->n_dofs(),
367 this->n_local_dofs(), send_list,
368 false, GHOSTED);
369 #else
370 current_local_solution->init(this->n_dofs());
371 #endif
372
373 if (_solution_projection)
374 solution->localize (*current_local_solution, send_list);
375
376 #endif // LIBMESH_ENABLE_AMR
377 }
378
379
380
prolong_vectors()381 void System::prolong_vectors ()
382 {
383 #ifdef LIBMESH_ENABLE_AMR
384 // Currently project_vector handles both restriction and prolongation
385 this->restrict_vectors();
386 #endif
387 }
388
389
390
reinit()391 void System::reinit ()
392 {
393 // project_vector handles vector initialization now
394 libmesh_assert_equal_to (solution->size(), current_local_solution->size());
395 }
396
397
reinit_constraints()398 void System::reinit_constraints()
399 {
400 #ifdef LIBMESH_ENABLE_CONSTRAINTS
401 get_dof_map().create_dof_constraints(_mesh, this->time);
402 user_constrain();
403 get_dof_map().process_constraints(_mesh);
404 #endif
405 get_dof_map().prepare_send_list();
406 }
407
408
update()409 void System::update ()
410 {
411 libmesh_assert(solution->closed());
412
413 const std::vector<dof_id_type> & send_list = _dof_map->get_send_list ();
414
415 // Check sizes
416 libmesh_assert_equal_to (current_local_solution->size(), solution->size());
417 // More processors than elements => empty send_list
418 // libmesh_assert (!send_list.empty());
419 libmesh_assert_less_equal (send_list.size(), solution->size());
420
421 // Create current_local_solution from solution. This will
422 // put a local copy of solution into current_local_solution.
423 // Only the necessary values (specified by the send_list)
424 // are copied to minimize communication
425 solution->localize (*current_local_solution, send_list);
426 }
427
428
429
re_update()430 void System::re_update ()
431 {
432 parallel_object_only();
433
434 // If this system is empty... don't do anything!
435 if (!this->n_vars())
436 return;
437
438 const std::vector<dof_id_type> & send_list = this->get_dof_map().get_send_list ();
439
440 // Check sizes
441 libmesh_assert_equal_to (current_local_solution->size(), solution->size());
442 // Not true with ghosted vectors
443 // libmesh_assert_equal_to (current_local_solution->local_size(), solution->size());
444 // libmesh_assert (!send_list.empty());
445 libmesh_assert_less_equal (send_list.size(), solution->size());
446
447 // Create current_local_solution from solution. This will
448 // put a local copy of solution into current_local_solution.
449 solution->localize (*current_local_solution, send_list);
450 }
451
452
453
restrict_solve_to(const SystemSubset * subset,const SubsetSolveMode)454 void System::restrict_solve_to (const SystemSubset * subset,
455 const SubsetSolveMode /*subset_solve_mode*/)
456 {
457 if (subset != nullptr)
458 libmesh_not_implemented();
459 }
460
461
462
assemble()463 void System::assemble ()
464 {
465 // Log how long the user's assembly code takes
466 LOG_SCOPE("assemble()", "System");
467
468 // Call the user-specified assembly function
469 this->user_assembly();
470 }
471
472
473
assemble_qoi(const QoISet & qoi_indices)474 void System::assemble_qoi (const QoISet & qoi_indices)
475 {
476 // Log how long the user's assembly code takes
477 LOG_SCOPE("assemble_qoi()", "System");
478
479 // Call the user-specified quantity of interest function
480 this->user_QOI(qoi_indices);
481 }
482
483
484
assemble_qoi_derivative(const QoISet & qoi_indices,bool include_liftfunc,bool apply_constraints)485 void System::assemble_qoi_derivative(const QoISet & qoi_indices,
486 bool include_liftfunc,
487 bool apply_constraints)
488 {
489 // Log how long the user's assembly code takes
490 LOG_SCOPE("assemble_qoi_derivative()", "System");
491
492 // Call the user-specified quantity of interest function
493 this->user_QOI_derivative(qoi_indices, include_liftfunc,
494 apply_constraints);
495 }
496
497
498
qoi_parameter_sensitivity(const QoISet & qoi_indices,const ParameterVector & parameters,SensitivityData & sensitivities)499 void System::qoi_parameter_sensitivity (const QoISet & qoi_indices,
500 const ParameterVector & parameters,
501 SensitivityData & sensitivities)
502 {
503 // Forward sensitivities are more efficient for Nq > Np
504 if (qoi_indices.size(*this) > parameters.size())
505 forward_qoi_parameter_sensitivity(qoi_indices, parameters, sensitivities);
506 // Adjoint sensitivities are more efficient for Np > Nq,
507 // and an adjoint may be more reusable than a forward
508 // solution sensitivity in the Np == Nq case.
509 else
510 adjoint_qoi_parameter_sensitivity(qoi_indices, parameters, sensitivities);
511 }
512
513
514
compare(const System & other_system,const Real threshold,const bool verbose)515 bool System::compare (const System & other_system,
516 const Real threshold,
517 const bool verbose) const
518 {
519 // we do not care for matrices, but for vectors
520 libmesh_assert (_is_initialized);
521 libmesh_assert (other_system._is_initialized);
522
523 if (verbose)
524 {
525 libMesh::out << " Systems \"" << _sys_name << "\"" << std::endl;
526 libMesh::out << " comparing matrices not supported." << std::endl;
527 libMesh::out << " comparing names...";
528 }
529
530 // compare the name: 0 means identical
531 const int name_result = _sys_name.compare(other_system.name());
532 if (verbose)
533 {
534 if (name_result == 0)
535 libMesh::out << " identical." << std::endl;
536 else
537 libMesh::out << " names not identical." << std::endl;
538 libMesh::out << " comparing solution vector...";
539 }
540
541
542 // compare the solution: -1 means identical
543 const int solu_result = solution->compare (*other_system.solution.get(),
544 threshold);
545
546 if (verbose)
547 {
548 if (solu_result == -1)
549 libMesh::out << " identical up to threshold." << std::endl;
550 else
551 libMesh::out << " first difference occurred at index = "
552 << solu_result << "." << std::endl;
553 }
554
555
556 // safety check, whether we handle at least the same number
557 // of vectors
558 std::vector<int> ov_result;
559
560 if (this->n_vectors() != other_system.n_vectors())
561 {
562 if (verbose)
563 {
564 libMesh::out << " Fatal difference. This system handles "
565 << this->n_vectors() << " add'l vectors," << std::endl
566 << " while the other system handles "
567 << other_system.n_vectors()
568 << " add'l vectors." << std::endl
569 << " Aborting comparison." << std::endl;
570 }
571 return false;
572 }
573 else if (this->n_vectors() == 0)
574 {
575 // there are no additional vectors...
576 ov_result.clear ();
577 }
578 else
579 {
580 // compare other vectors
581 for (auto & pr : _vectors)
582 {
583 if (verbose)
584 libMesh::out << " comparing vector \""
585 << pr.first << "\" ...";
586
587 // assume they have the same name
588 const NumericVector<Number> & other_system_vector =
589 other_system.get_vector(pr.first);
590
591 ov_result.push_back(pr.second->compare (other_system_vector,
592 threshold));
593
594 if (verbose)
595 {
596 if (ov_result[ov_result.size()-1] == -1)
597 libMesh::out << " identical up to threshold." << std::endl;
598 else
599 libMesh::out << " first difference occurred at" << std::endl
600 << " index = " << ov_result[ov_result.size()-1] << "." << std::endl;
601 }
602 }
603 } // finished comparing additional vectors
604
605
606 bool overall_result;
607
608 // sum up the results
609 if ((name_result==0) && (solu_result==-1))
610 {
611 if (ov_result.size()==0)
612 overall_result = true;
613 else
614 {
615 bool ov_identical;
616 unsigned int n = 0;
617 do
618 {
619 ov_identical = (ov_result[n]==-1);
620 n++;
621 }
622 while (ov_identical && n<ov_result.size());
623 overall_result = ov_identical;
624 }
625 }
626 else
627 overall_result = false;
628
629 if (verbose)
630 {
631 libMesh::out << " finished comparisons, ";
632 if (overall_result)
633 libMesh::out << "found no differences." << std::endl << std::endl;
634 else
635 libMesh::out << "found differences." << std::endl << std::endl;
636 }
637
638 return overall_result;
639 }
640
641
642
update_global_solution(std::vector<Number> & global_soln)643 void System::update_global_solution (std::vector<Number> & global_soln) const
644 {
645 global_soln.resize (solution->size());
646
647 solution->localize (global_soln);
648 }
649
650
651
update_global_solution(std::vector<Number> & global_soln,const processor_id_type dest_proc)652 void System::update_global_solution (std::vector<Number> & global_soln,
653 const processor_id_type dest_proc) const
654 {
655 global_soln.resize (solution->size());
656
657 solution->localize_to_one (global_soln, dest_proc);
658 }
659
660
661
add_vector(const std::string & vec_name,const bool projections,const ParallelType type)662 NumericVector<Number> & System::add_vector (const std::string & vec_name,
663 const bool projections,
664 const ParallelType type)
665 {
666 // Return the vector if it is already there.
667 if (this->have_vector(vec_name))
668 return *(_vectors[vec_name]);
669
670 // Otherwise build the vector
671 NumericVector<Number> * buf = NumericVector<Number>::build(this->comm()).release();
672 _vectors.emplace(vec_name, buf);
673 _vector_projections.emplace(vec_name, projections);
674 _vector_types.emplace(vec_name, type);
675
676 // Vectors are primal by default
677 _vector_is_adjoint.emplace(vec_name, -1);
678
679 // Initialize it if necessary
680 if (_is_initialized)
681 {
682 if (type == GHOSTED)
683 {
684 #ifdef LIBMESH_ENABLE_GHOSTED
685 buf->init (this->n_dofs(), this->n_local_dofs(),
686 _dof_map->get_send_list(), false,
687 GHOSTED);
688 #else
689 libmesh_error_msg("Cannot initialize ghosted vectors when they are not enabled.");
690 #endif
691 }
692 else
693 buf->init (this->n_dofs(), this->n_local_dofs(), false, type);
694 }
695
696 return *buf;
697 }
698
remove_vector(const std::string & vec_name)699 void System::remove_vector (const std::string & vec_name)
700 {
701 vectors_iterator pos = _vectors.find(vec_name);
702
703 //Return if the vector does not exist
704 if (pos == _vectors.end())
705 return;
706
707 delete pos->second;
708
709 _vectors.erase(pos);
710
711 _vector_projections.erase(vec_name);
712 _vector_is_adjoint.erase(vec_name);
713 _vector_types.erase(vec_name);
714 }
715
request_vector(const std::string & vec_name)716 const NumericVector<Number> * System::request_vector (const std::string & vec_name) const
717 {
718 const_vectors_iterator pos = _vectors.find(vec_name);
719
720 if (pos == _vectors.end())
721 return nullptr;
722
723 return pos->second;
724 }
725
726
727
request_vector(const std::string & vec_name)728 NumericVector<Number> * System::request_vector (const std::string & vec_name)
729 {
730 vectors_iterator pos = _vectors.find(vec_name);
731
732 if (pos == _vectors.end())
733 return nullptr;
734
735 return pos->second;
736 }
737
738
739
request_vector(const unsigned int vec_num)740 const NumericVector<Number> * System::request_vector (const unsigned int vec_num) const
741 {
742 const_vectors_iterator v = vectors_begin();
743 const_vectors_iterator v_end = vectors_end();
744 unsigned int num = 0;
745 while ((num<vec_num) && (v!=v_end))
746 {
747 num++;
748 ++v;
749 }
750 if (v==v_end)
751 return nullptr;
752 return v->second;
753 }
754
755
756
request_vector(const unsigned int vec_num)757 NumericVector<Number> * System::request_vector (const unsigned int vec_num)
758 {
759 vectors_iterator v = vectors_begin();
760 vectors_iterator v_end = vectors_end();
761 unsigned int num = 0;
762 while ((num<vec_num) && (v!=v_end))
763 {
764 num++;
765 ++v;
766 }
767 if (v==v_end)
768 return nullptr;
769 return v->second;
770 }
771
772
773
get_vector(const std::string & vec_name)774 const NumericVector<Number> & System::get_vector (const std::string & vec_name) const
775 {
776 return *(libmesh_map_find(_vectors, vec_name));
777 }
778
779
780
get_vector(const std::string & vec_name)781 NumericVector<Number> & System::get_vector (const std::string & vec_name)
782 {
783 return *(libmesh_map_find(_vectors, vec_name));
784 }
785
786
787
get_vector(const unsigned int vec_num)788 const NumericVector<Number> & System::get_vector (const unsigned int vec_num) const
789 {
790 const_vectors_iterator v = vectors_begin();
791 const_vectors_iterator v_end = vectors_end();
792 unsigned int num = 0;
793 while ((num<vec_num) && (v!=v_end))
794 {
795 num++;
796 ++v;
797 }
798 libmesh_assert (v != v_end);
799 return *(v->second);
800 }
801
802
803
get_vector(const unsigned int vec_num)804 NumericVector<Number> & System::get_vector (const unsigned int vec_num)
805 {
806 vectors_iterator v = vectors_begin();
807 vectors_iterator v_end = vectors_end();
808 unsigned int num = 0;
809 while ((num<vec_num) && (v!=v_end))
810 {
811 num++;
812 ++v;
813 }
814 libmesh_assert (v != v_end);
815 return *(v->second);
816 }
817
818
819
vector_name(const unsigned int vec_num)820 const std::string & System::vector_name (const unsigned int vec_num) const
821 {
822 const_vectors_iterator v = vectors_begin();
823 const_vectors_iterator v_end = vectors_end();
824 unsigned int num = 0;
825 while ((num<vec_num) && (v!=v_end))
826 {
827 num++;
828 ++v;
829 }
830 libmesh_assert (v != v_end);
831 return v->first;
832 }
833
vector_name(const NumericVector<Number> & vec_reference)834 const std::string & System::vector_name (const NumericVector<Number> & vec_reference) const
835 {
836 const_vectors_iterator v = vectors_begin();
837 const_vectors_iterator v_end = vectors_end();
838
839 for (; v != v_end; ++v)
840 {
841 // Check if the current vector is the one whose name we want
842 if (&vec_reference == v->second)
843 break; // exit loop if it is
844 }
845
846 // Before returning, make sure we didnt loop till the end and not find any match
847 libmesh_assert (v != v_end);
848
849 // Return the string associated with the current vector
850 return v->first;
851 }
852
853
854
set_vector_preservation(const std::string & vec_name,bool preserve)855 void System::set_vector_preservation (const std::string & vec_name,
856 bool preserve)
857 {
858 _vector_projections[vec_name] = preserve;
859 }
860
861
862
vector_preservation(const std::string & vec_name)863 bool System::vector_preservation (const std::string & vec_name) const
864 {
865 if (_vector_projections.find(vec_name) == _vector_projections.end())
866 return false;
867
868 return _vector_projections.find(vec_name)->second;
869 }
870
871
872
set_vector_as_adjoint(const std::string & vec_name,int qoi_num)873 void System::set_vector_as_adjoint (const std::string & vec_name,
874 int qoi_num)
875 {
876 // We reserve -1 for vectors which get primal constraints, -2 for
877 // vectors which get no constraints
878 libmesh_assert_greater_equal(qoi_num, -2);
879 _vector_is_adjoint[vec_name] = qoi_num;
880 }
881
882
883
vector_is_adjoint(const std::string & vec_name)884 int System::vector_is_adjoint (const std::string & vec_name) const
885 {
886 libmesh_assert(_vector_is_adjoint.find(vec_name) !=
887 _vector_is_adjoint.end());
888
889 return _vector_is_adjoint.find(vec_name)->second;
890 }
891
892
893
add_sensitivity_solution(unsigned int i)894 NumericVector<Number> & System::add_sensitivity_solution (unsigned int i)
895 {
896 std::ostringstream sensitivity_name;
897 sensitivity_name << "sensitivity_solution" << i;
898
899 return this->add_vector(sensitivity_name.str());
900 }
901
902
903
get_sensitivity_solution(unsigned int i)904 NumericVector<Number> & System::get_sensitivity_solution (unsigned int i)
905 {
906 std::ostringstream sensitivity_name;
907 sensitivity_name << "sensitivity_solution" << i;
908
909 return this->get_vector(sensitivity_name.str());
910 }
911
912
913
get_sensitivity_solution(unsigned int i)914 const NumericVector<Number> & System::get_sensitivity_solution (unsigned int i) const
915 {
916 std::ostringstream sensitivity_name;
917 sensitivity_name << "sensitivity_solution" << i;
918
919 return this->get_vector(sensitivity_name.str());
920 }
921
922
923
add_weighted_sensitivity_solution()924 NumericVector<Number> & System::add_weighted_sensitivity_solution ()
925 {
926 return this->add_vector("weighted_sensitivity_solution");
927 }
928
929
930
get_weighted_sensitivity_solution()931 NumericVector<Number> & System::get_weighted_sensitivity_solution ()
932 {
933 return this->get_vector("weighted_sensitivity_solution");
934 }
935
936
937
get_weighted_sensitivity_solution()938 const NumericVector<Number> & System::get_weighted_sensitivity_solution () const
939 {
940 return this->get_vector("weighted_sensitivity_solution");
941 }
942
943
944
add_adjoint_solution(unsigned int i)945 NumericVector<Number> & System::add_adjoint_solution (unsigned int i)
946 {
947 std::ostringstream adjoint_name;
948 adjoint_name << "adjoint_solution" << i;
949
950 NumericVector<Number> & returnval = this->add_vector(adjoint_name.str());
951 this->set_vector_as_adjoint(adjoint_name.str(), i);
952 return returnval;
953 }
954
955
956
get_adjoint_solution(unsigned int i)957 NumericVector<Number> & System::get_adjoint_solution (unsigned int i)
958 {
959 std::ostringstream adjoint_name;
960 adjoint_name << "adjoint_solution" << i;
961
962 return this->get_vector(adjoint_name.str());
963 }
964
965
966
get_adjoint_solution(unsigned int i)967 const NumericVector<Number> & System::get_adjoint_solution (unsigned int i) const
968 {
969 std::ostringstream adjoint_name;
970 adjoint_name << "adjoint_solution" << i;
971
972 return this->get_vector(adjoint_name.str());
973 }
974
975
976
add_weighted_sensitivity_adjoint_solution(unsigned int i)977 NumericVector<Number> & System::add_weighted_sensitivity_adjoint_solution (unsigned int i)
978 {
979 std::ostringstream adjoint_name;
980 adjoint_name << "weighted_sensitivity_adjoint_solution" << i;
981
982 NumericVector<Number> & returnval = this->add_vector(adjoint_name.str());
983 this->set_vector_as_adjoint(adjoint_name.str(), i);
984 return returnval;
985 }
986
987
988
get_weighted_sensitivity_adjoint_solution(unsigned int i)989 NumericVector<Number> & System::get_weighted_sensitivity_adjoint_solution (unsigned int i)
990 {
991 std::ostringstream adjoint_name;
992 adjoint_name << "weighted_sensitivity_adjoint_solution" << i;
993
994 return this->get_vector(adjoint_name.str());
995 }
996
997
998
get_weighted_sensitivity_adjoint_solution(unsigned int i)999 const NumericVector<Number> & System::get_weighted_sensitivity_adjoint_solution (unsigned int i) const
1000 {
1001 std::ostringstream adjoint_name;
1002 adjoint_name << "weighted_sensitivity_adjoint_solution" << i;
1003
1004 return this->get_vector(adjoint_name.str());
1005 }
1006
1007
1008
add_adjoint_rhs(unsigned int i)1009 NumericVector<Number> & System::add_adjoint_rhs (unsigned int i)
1010 {
1011 std::ostringstream adjoint_rhs_name;
1012 adjoint_rhs_name << "adjoint_rhs" << i;
1013
1014 return this->add_vector(adjoint_rhs_name.str(), false);
1015 }
1016
1017
1018
get_adjoint_rhs(unsigned int i)1019 NumericVector<Number> & System::get_adjoint_rhs (unsigned int i)
1020 {
1021 std::ostringstream adjoint_rhs_name;
1022 adjoint_rhs_name << "adjoint_rhs" << i;
1023
1024 return this->get_vector(adjoint_rhs_name.str());
1025 }
1026
1027
1028
get_adjoint_rhs(unsigned int i)1029 const NumericVector<Number> & System::get_adjoint_rhs (unsigned int i) const
1030 {
1031 std::ostringstream adjoint_rhs_name;
1032 adjoint_rhs_name << "adjoint_rhs" << i;
1033
1034 return this->get_vector(adjoint_rhs_name.str());
1035 }
1036
1037
1038
add_sensitivity_rhs(unsigned int i)1039 NumericVector<Number> & System::add_sensitivity_rhs (unsigned int i)
1040 {
1041 std::ostringstream sensitivity_rhs_name;
1042 sensitivity_rhs_name << "sensitivity_rhs" << i;
1043
1044 return this->add_vector(sensitivity_rhs_name.str(), false);
1045 }
1046
1047
1048
get_sensitivity_rhs(unsigned int i)1049 NumericVector<Number> & System::get_sensitivity_rhs (unsigned int i)
1050 {
1051 std::ostringstream sensitivity_rhs_name;
1052 sensitivity_rhs_name << "sensitivity_rhs" << i;
1053
1054 return this->get_vector(sensitivity_rhs_name.str());
1055 }
1056
1057
1058
get_sensitivity_rhs(unsigned int i)1059 const NumericVector<Number> & System::get_sensitivity_rhs (unsigned int i) const
1060 {
1061 std::ostringstream sensitivity_rhs_name;
1062 sensitivity_rhs_name << "sensitivity_rhs" << i;
1063
1064 return this->get_vector(sensitivity_rhs_name.str());
1065 }
1066
1067
1068
add_variable(const std::string & var,const FEType & type,const std::set<subdomain_id_type> * const active_subdomains)1069 unsigned int System::add_variable (const std::string & var,
1070 const FEType & type,
1071 const std::set<subdomain_id_type> * const active_subdomains)
1072 {
1073 libmesh_assert(!this->is_initialized());
1074
1075 // Make sure the variable isn't there already
1076 // or if it is, that it's the type we want
1077 for (auto v : make_range(this->n_vars()))
1078 if (this->variable_name(v) == var)
1079 {
1080 if (this->variable_type(v) == type)
1081 return _variables[v].number();
1082
1083 libmesh_error_msg("ERROR: incompatible variable " << var << " has already been added for this system!");
1084 }
1085
1086 // Optimize for VariableGroups here - if the user is adding multiple
1087 // variables of the same FEType and subdomain restriction, catch
1088 // that here and add them as members of the same VariableGroup.
1089 //
1090 // start by setting this flag to whatever the user has requested
1091 // and then consider the conditions which should negate it.
1092 bool should_be_in_vg = this->identify_variable_groups();
1093
1094 // No variable groups, nothing to add to
1095 if (!this->n_variable_groups())
1096 should_be_in_vg = false;
1097
1098 else
1099 {
1100 VariableGroup & vg(_variable_groups.back());
1101
1102 // get a pointer to their subdomain restriction, if any.
1103 const std::set<subdomain_id_type> * const
1104 their_active_subdomains (vg.implicitly_active() ?
1105 nullptr : &vg.active_subdomains());
1106
1107 // Different types?
1108 if (vg.type() != type)
1109 should_be_in_vg = false;
1110
1111 // they are restricted, we aren't?
1112 if (their_active_subdomains &&
1113 (!active_subdomains || (active_subdomains && active_subdomains->empty())))
1114 should_be_in_vg = false;
1115
1116 // they aren't restricted, we are?
1117 if (!their_active_subdomains && (active_subdomains && !active_subdomains->empty()))
1118 should_be_in_vg = false;
1119
1120 if (their_active_subdomains && active_subdomains)
1121 // restricted to different sets?
1122 if (*their_active_subdomains != *active_subdomains)
1123 should_be_in_vg = false;
1124
1125 // OK, after all that, append the variable to the vg if none of the conditions
1126 // were violated
1127 if (should_be_in_vg)
1128 {
1129 const unsigned short curr_n_vars = cast_int<unsigned short>
1130 (this->n_vars());
1131
1132 vg.append (var);
1133
1134 _variables.push_back(vg(vg.n_variables()-1));
1135 _variable_numbers[var] = curr_n_vars;
1136 return curr_n_vars;
1137 }
1138 }
1139
1140 // otherwise, fall back to adding a single variable group
1141 return this->add_variables (std::vector<std::string>(1, var),
1142 type,
1143 active_subdomains);
1144 }
1145
1146
1147
add_variable(const std::string & var,const Order order,const FEFamily family,const std::set<subdomain_id_type> * const active_subdomains)1148 unsigned int System::add_variable (const std::string & var,
1149 const Order order,
1150 const FEFamily family,
1151 const std::set<subdomain_id_type> * const active_subdomains)
1152 {
1153 return this->add_variable(var,
1154 FEType(order, family),
1155 active_subdomains);
1156 }
1157
1158
1159
add_variables(const std::vector<std::string> & vars,const FEType & type,const std::set<subdomain_id_type> * const active_subdomains)1160 unsigned int System::add_variables (const std::vector<std::string> & vars,
1161 const FEType & type,
1162 const std::set<subdomain_id_type> * const active_subdomains)
1163 {
1164 libmesh_assert(!this->is_initialized());
1165
1166 // Make sure the variable isn't there already
1167 // or if it is, that it's the type we want
1168 for (auto ovar : vars)
1169 for (auto v : make_range(this->n_vars()))
1170 if (this->variable_name(v) == ovar)
1171 {
1172 if (this->variable_type(v) == type)
1173 return _variables[v].number();
1174
1175 libmesh_error_msg("ERROR: incompatible variable " << ovar << " has already been added for this system!");
1176 }
1177
1178 // Optimize for VariableGroups here - if the user is adding multiple
1179 // variables of the same FEType and subdomain restriction, catch
1180 // that here and add them as members of the same VariableGroup.
1181 //
1182 // start by setting this flag to whatever the user has requested
1183 // and then consider the conditions which should negate it.
1184 bool should_be_in_vg = this->identify_variable_groups();
1185
1186 // No variable groups, nothing to add to
1187 if (!this->n_variable_groups())
1188 should_be_in_vg = false;
1189 else
1190 {
1191 VariableGroup & vg(_variable_groups.back());
1192
1193 // get a pointer to their subdomain restriction, if any.
1194 const std::set<subdomain_id_type> * const
1195 their_active_subdomains (vg.implicitly_active() ?
1196 nullptr : &vg.active_subdomains());
1197
1198 // Different types?
1199 if (vg.type() != type)
1200 should_be_in_vg = false;
1201
1202 // they are restricted, we aren't?
1203 if (their_active_subdomains &&
1204 (!active_subdomains || (active_subdomains && active_subdomains->empty())))
1205 should_be_in_vg = false;
1206
1207 // they aren't restricted, we are?
1208 if (!their_active_subdomains && (active_subdomains && !active_subdomains->empty()))
1209 should_be_in_vg = false;
1210
1211 if (their_active_subdomains && active_subdomains)
1212 // restricted to different sets?
1213 if (*their_active_subdomains != *active_subdomains)
1214 should_be_in_vg = false;
1215
1216 // If after all that none of the conditions were violated,
1217 // append the variables to the vg and we're done
1218 if (should_be_in_vg)
1219 {
1220 unsigned short curr_n_vars = cast_int<unsigned short>
1221 (this->n_vars());
1222
1223 for (auto ovar : vars)
1224 {
1225 curr_n_vars = cast_int<unsigned short> (this->n_vars());
1226
1227 vg.append (ovar);
1228
1229 _variables.push_back(vg(vg.n_variables()-1));
1230 _variable_numbers[ovar] = curr_n_vars;
1231 }
1232 return curr_n_vars;
1233 }
1234 }
1235
1236 const unsigned short curr_n_vars = cast_int<unsigned short>
1237 (this->n_vars());
1238
1239 const unsigned int next_first_component = this->n_components();
1240
1241 // We weren't able to add to an existing variable group, so
1242 // add a new variable group to the list
1243 _variable_groups.push_back((active_subdomains == nullptr) ?
1244 VariableGroup(this, vars, curr_n_vars,
1245 next_first_component, type) :
1246 VariableGroup(this, vars, curr_n_vars,
1247 next_first_component, type, *active_subdomains));
1248
1249 const VariableGroup & vg (_variable_groups.back());
1250
1251 // Add each component of the group individually
1252 for (auto v : make_range(vars.size()))
1253 {
1254 _variables.push_back (vg(v));
1255 _variable_numbers[vars[v]] = cast_int<unsigned short>
1256 (curr_n_vars+v);
1257 }
1258
1259 libmesh_assert_equal_to ((curr_n_vars+vars.size()), this->n_vars());
1260
1261 // BSK - Defer this now to System::init_data() so we can detect
1262 // VariableGroups 12/28/2012
1263 // // Add the variable group to the _dof_map
1264 // _dof_map->add_variable_group (vg);
1265
1266 // Return the number of the new variable
1267 return cast_int<unsigned int>(curr_n_vars+vars.size()-1);
1268 }
1269
1270
1271
add_variables(const std::vector<std::string> & vars,const Order order,const FEFamily family,const std::set<subdomain_id_type> * const active_subdomains)1272 unsigned int System::add_variables (const std::vector<std::string> & vars,
1273 const Order order,
1274 const FEFamily family,
1275 const std::set<subdomain_id_type> * const active_subdomains)
1276 {
1277 return this->add_variables(vars,
1278 FEType(order, family),
1279 active_subdomains);
1280 }
1281
1282
1283
has_variable(const std::string & var)1284 bool System::has_variable (const std::string & var) const
1285 {
1286 return _variable_numbers.count(var);
1287 }
1288
1289
1290
variable_number(const std::string & var)1291 unsigned short int System::variable_number (const std::string & var) const
1292 {
1293 auto var_num = libmesh_map_find(_variable_numbers, var);
1294 libmesh_assert_equal_to (_variables[var_num].name(), var);
1295 return var_num;
1296 }
1297
1298
get_all_variable_numbers(std::vector<unsigned int> & all_variable_numbers)1299 void System::get_all_variable_numbers(std::vector<unsigned int> & all_variable_numbers) const
1300 {
1301 all_variable_numbers.resize(n_vars());
1302
1303 // Make sure the variable exists
1304 std::map<std::string, unsigned short int>::const_iterator
1305 it = _variable_numbers.begin();
1306 std::map<std::string, unsigned short int>::const_iterator
1307 it_end = _variable_numbers.end();
1308
1309 unsigned int count = 0;
1310 for ( ; it != it_end; ++it)
1311 {
1312 all_variable_numbers[count] = it->second;
1313 count++;
1314 }
1315 }
1316
1317
local_dof_indices(const unsigned int var,std::set<dof_id_type> & var_indices)1318 void System::local_dof_indices(const unsigned int var,
1319 std::set<dof_id_type> & var_indices) const
1320 {
1321 // Make sure the set is clear
1322 var_indices.clear();
1323
1324 std::vector<dof_id_type> dof_indices;
1325
1326 const dof_id_type
1327 first_local = this->get_dof_map().first_dof(),
1328 end_local = this->get_dof_map().end_dof();
1329
1330 // Begin the loop over the elements
1331 for (const auto & elem : this->get_mesh().active_local_element_ptr_range())
1332 {
1333 this->get_dof_map().dof_indices (elem, dof_indices, var);
1334
1335 for (dof_id_type dof : dof_indices)
1336 //If the dof is owned by the local processor
1337 if (first_local <= dof && dof < end_local)
1338 var_indices.insert(dof);
1339 }
1340
1341 // we may have missed assigning DOFs to nodes that we own
1342 // but to which we have no connected elements matching our
1343 // variable restriction criterion. this will happen, for example,
1344 // if variable V is restricted to subdomain S. We may not own
1345 // any elements which live in S, but we may own nodes which are
1346 // *connected* to elements which do.
1347 for (const auto & node : this->get_mesh().local_node_ptr_range())
1348 {
1349 libmesh_assert(node);
1350 this->get_dof_map().dof_indices (node, dof_indices, var);
1351 for (auto dof : dof_indices)
1352 if (first_local <= dof && dof < end_local)
1353 var_indices.insert(dof);
1354 }
1355 }
1356
1357
1358
zero_variable(NumericVector<Number> & v,unsigned int var_num)1359 void System::zero_variable (NumericVector<Number> & v,
1360 unsigned int var_num) const
1361 {
1362 /* Make sure the call makes sense. */
1363 libmesh_assert_less (var_num, this->n_vars());
1364
1365 /* Get a reference to the mesh. */
1366 const MeshBase & mesh = this->get_mesh();
1367
1368 /* Check which system we are. */
1369 const unsigned int sys_num = this->number();
1370
1371 // Loop over nodes.
1372 for (const auto & node : mesh.local_node_ptr_range())
1373 {
1374 unsigned int n_comp = node->n_comp(sys_num,var_num);
1375 for (unsigned int i=0; i<n_comp; i++)
1376 {
1377 const dof_id_type index = node->dof_number(sys_num,var_num,i);
1378 v.set(index,0.0);
1379 }
1380 }
1381
1382 // Loop over elements.
1383 for (const auto & elem : mesh.active_local_element_ptr_range())
1384 {
1385 unsigned int n_comp = elem->n_comp(sys_num,var_num);
1386 for (unsigned int i=0; i<n_comp; i++)
1387 {
1388 const dof_id_type index = elem->dof_number(sys_num,var_num,i);
1389 v.set(index,0.0);
1390 }
1391 }
1392 }
1393
1394
1395
discrete_var_norm(const NumericVector<Number> & v,unsigned int var,FEMNormType norm_type)1396 Real System::discrete_var_norm(const NumericVector<Number> & v,
1397 unsigned int var,
1398 FEMNormType norm_type) const
1399 {
1400 std::set<dof_id_type> var_indices;
1401 local_dof_indices(var, var_indices);
1402
1403 if (norm_type == DISCRETE_L1)
1404 return v.subset_l1_norm(var_indices);
1405 if (norm_type == DISCRETE_L2)
1406 return v.subset_l2_norm(var_indices);
1407 if (norm_type == DISCRETE_L_INF)
1408 return v.subset_linfty_norm(var_indices);
1409 else
1410 libmesh_error_msg("Invalid norm_type = " << norm_type);
1411 }
1412
1413
1414
calculate_norm(const NumericVector<Number> & v,unsigned int var,FEMNormType norm_type,std::set<unsigned int> * skip_dimensions)1415 Real System::calculate_norm(const NumericVector<Number> & v,
1416 unsigned int var,
1417 FEMNormType norm_type,
1418 std::set<unsigned int> * skip_dimensions) const
1419 {
1420 //short circuit to save time
1421 if (norm_type == DISCRETE_L1 ||
1422 norm_type == DISCRETE_L2 ||
1423 norm_type == DISCRETE_L_INF)
1424 return discrete_var_norm(v,var,norm_type);
1425
1426 // Not a discrete norm
1427 std::vector<FEMNormType> norms(this->n_vars(), L2);
1428 std::vector<Real> weights(this->n_vars(), 0.0);
1429 norms[var] = norm_type;
1430 weights[var] = 1.0;
1431 Real val = this->calculate_norm(v, SystemNorm(norms, weights), skip_dimensions);
1432 return val;
1433 }
1434
1435
1436
calculate_norm(const NumericVector<Number> & v,const SystemNorm & norm,std::set<unsigned int> * skip_dimensions)1437 Real System::calculate_norm(const NumericVector<Number> & v,
1438 const SystemNorm & norm,
1439 std::set<unsigned int> * skip_dimensions) const
1440 {
1441 // This function must be run on all processors at once
1442 parallel_object_only();
1443
1444 LOG_SCOPE ("calculate_norm()", "System");
1445
1446 // Zero the norm before summation
1447 Real v_norm = 0.;
1448
1449 if (norm.is_discrete())
1450 {
1451 //Check to see if all weights are 1.0 and all types are equal
1452 FEMNormType norm_type0 = norm.type(0);
1453 unsigned int check_var = 0, check_end = this->n_vars();
1454 for (; check_var != check_end; ++check_var)
1455 if ((norm.weight(check_var) != 1.0) || (norm.type(check_var) != norm_type0))
1456 break;
1457
1458 //All weights were 1.0 so just do the full vector discrete norm
1459 if (check_var == this->n_vars())
1460 {
1461 if (norm_type0 == DISCRETE_L1)
1462 return v.l1_norm();
1463 if (norm_type0 == DISCRETE_L2)
1464 return v.l2_norm();
1465 if (norm_type0 == DISCRETE_L_INF)
1466 return v.linfty_norm();
1467 else
1468 libmesh_error_msg("Invalid norm_type0 = " << norm_type0);
1469 }
1470
1471 for (auto var : make_range(this->n_vars()))
1472 {
1473 // Skip any variables we don't need to integrate
1474 if (norm.weight(var) == 0.0)
1475 continue;
1476
1477 v_norm += norm.weight(var) * discrete_var_norm(v, var, norm.type(var));
1478 }
1479
1480 return v_norm;
1481 }
1482
1483 // Localize the potentially parallel vector
1484 std::unique_ptr<NumericVector<Number>> local_v = NumericVector<Number>::build(this->comm());
1485 local_v->init(v.size(), v.local_size(), _dof_map->get_send_list(),
1486 true, GHOSTED);
1487 v.localize (*local_v, _dof_map->get_send_list());
1488
1489 // I'm not sure how best to mix Hilbert norms on some variables (for
1490 // which we'll want to square then sum then square root) with norms
1491 // like L_inf (for which we'll just want to take an absolute value
1492 // and then sum).
1493 bool using_hilbert_norm = true,
1494 using_nonhilbert_norm = true;
1495
1496 // Loop over all variables
1497 for (auto var : make_range(this->n_vars()))
1498 {
1499 // Skip any variables we don't need to integrate
1500 Real norm_weight_sq = norm.weight_sq(var);
1501 if (norm_weight_sq == 0.0)
1502 continue;
1503 Real norm_weight = norm.weight(var);
1504
1505 // Check for unimplemented norms (rather than just returning 0).
1506 FEMNormType norm_type = norm.type(var);
1507 if ((norm_type==H1) ||
1508 (norm_type==H2) ||
1509 (norm_type==L2) ||
1510 (norm_type==H1_SEMINORM) ||
1511 (norm_type==H2_SEMINORM))
1512 {
1513 if (!using_hilbert_norm)
1514 libmesh_not_implemented();
1515 using_nonhilbert_norm = false;
1516 }
1517 else if ((norm_type==L1) ||
1518 (norm_type==L_INF) ||
1519 (norm_type==W1_INF_SEMINORM) ||
1520 (norm_type==W2_INF_SEMINORM))
1521 {
1522 if (!using_nonhilbert_norm)
1523 libmesh_not_implemented();
1524 using_hilbert_norm = false;
1525 }
1526 else
1527 libmesh_not_implemented();
1528
1529 const FEType & fe_type = this->get_dof_map().variable_type(var);
1530
1531 // Allow space for dims 0-3, even if we don't use them all
1532 std::vector<std::unique_ptr<FEBase>> fe_ptrs(4);
1533 std::vector<std::unique_ptr<QBase>> q_rules(4);
1534
1535 const std::set<unsigned char> & elem_dims = _mesh.elem_dimensions();
1536
1537 // Prepare finite elements for each dimension present in the mesh
1538 for (const auto & dim : elem_dims)
1539 {
1540 if (skip_dimensions && skip_dimensions->find(dim) != skip_dimensions->end())
1541 continue;
1542
1543 // Construct quadrature and finite element objects
1544 q_rules[dim] = fe_type.default_quadrature_rule (dim);
1545 fe_ptrs[dim] = FEBase::build(dim, fe_type);
1546
1547 // Attach quadrature rule to FE object
1548 fe_ptrs[dim]->attach_quadrature_rule (q_rules[dim].get());
1549 }
1550
1551 std::vector<dof_id_type> dof_indices;
1552
1553 // Begin the loop over the elements
1554 for (const auto & elem : this->get_mesh().active_local_element_ptr_range())
1555 {
1556 const unsigned int dim = elem->dim();
1557
1558 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
1559
1560 // One way for implementing this would be to exchange the fe with the FEInterface- class.
1561 // However, it needs to be discussed whether integral-norms make sense for infinite elements.
1562 // or in which sense they could make sense.
1563 if (elem->infinite() )
1564 libmesh_not_implemented();
1565
1566 #endif
1567
1568 if (skip_dimensions && skip_dimensions->find(dim) != skip_dimensions->end())
1569 continue;
1570
1571 FEBase * fe = fe_ptrs[dim].get();
1572 QBase * qrule = q_rules[dim].get();
1573 libmesh_assert(fe);
1574 libmesh_assert(qrule);
1575
1576 const std::vector<Real> & JxW = fe->get_JxW();
1577 const std::vector<std::vector<Real>> * phi = nullptr;
1578 if (norm_type == H1 ||
1579 norm_type == H2 ||
1580 norm_type == L2 ||
1581 norm_type == L1 ||
1582 norm_type == L_INF)
1583 phi = &(fe->get_phi());
1584
1585 const std::vector<std::vector<RealGradient>> * dphi = nullptr;
1586 if (norm_type == H1 ||
1587 norm_type == H2 ||
1588 norm_type == H1_SEMINORM ||
1589 norm_type == W1_INF_SEMINORM)
1590 dphi = &(fe->get_dphi());
1591 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1592 const std::vector<std::vector<RealTensor>> * d2phi = nullptr;
1593 if (norm_type == H2 ||
1594 norm_type == H2_SEMINORM ||
1595 norm_type == W2_INF_SEMINORM)
1596 d2phi = &(fe->get_d2phi());
1597 #endif
1598
1599 fe->reinit (elem);
1600
1601 this->get_dof_map().dof_indices (elem, dof_indices, var);
1602
1603 const unsigned int n_qp = qrule->n_points();
1604
1605 const unsigned int n_sf = cast_int<unsigned int>
1606 (dof_indices.size());
1607
1608 // Begin the loop over the Quadrature points.
1609 for (unsigned int qp=0; qp<n_qp; qp++)
1610 {
1611 if (norm_type == L1)
1612 {
1613 Number u_h = 0.;
1614 for (unsigned int i=0; i != n_sf; ++i)
1615 u_h += (*phi)[i][qp] * (*local_v)(dof_indices[i]);
1616 v_norm += norm_weight *
1617 JxW[qp] * std::abs(u_h);
1618 }
1619
1620 if (norm_type == L_INF)
1621 {
1622 Number u_h = 0.;
1623 for (unsigned int i=0; i != n_sf; ++i)
1624 u_h += (*phi)[i][qp] * (*local_v)(dof_indices[i]);
1625 v_norm = std::max(v_norm, norm_weight * std::abs(u_h));
1626 }
1627
1628 if (norm_type == H1 ||
1629 norm_type == H2 ||
1630 norm_type == L2)
1631 {
1632 Number u_h = 0.;
1633 for (unsigned int i=0; i != n_sf; ++i)
1634 u_h += (*phi)[i][qp] * (*local_v)(dof_indices[i]);
1635 v_norm += norm_weight_sq *
1636 JxW[qp] * TensorTools::norm_sq(u_h);
1637 }
1638
1639 if (norm_type == H1 ||
1640 norm_type == H2 ||
1641 norm_type == H1_SEMINORM)
1642 {
1643 Gradient grad_u_h;
1644 for (unsigned int i=0; i != n_sf; ++i)
1645 grad_u_h.add_scaled((*dphi)[i][qp], (*local_v)(dof_indices[i]));
1646 v_norm += norm_weight_sq *
1647 JxW[qp] * grad_u_h.norm_sq();
1648 }
1649
1650 if (norm_type == W1_INF_SEMINORM)
1651 {
1652 Gradient grad_u_h;
1653 for (unsigned int i=0; i != n_sf; ++i)
1654 grad_u_h.add_scaled((*dphi)[i][qp], (*local_v)(dof_indices[i]));
1655 v_norm = std::max(v_norm, norm_weight * grad_u_h.norm());
1656 }
1657
1658 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1659 if (norm_type == H2 ||
1660 norm_type == H2_SEMINORM)
1661 {
1662 Tensor hess_u_h;
1663 for (unsigned int i=0; i != n_sf; ++i)
1664 hess_u_h.add_scaled((*d2phi)[i][qp], (*local_v)(dof_indices[i]));
1665 v_norm += norm_weight_sq *
1666 JxW[qp] * hess_u_h.norm_sq();
1667 }
1668
1669 if (norm_type == W2_INF_SEMINORM)
1670 {
1671 Tensor hess_u_h;
1672 for (unsigned int i=0; i != n_sf; ++i)
1673 hess_u_h.add_scaled((*d2phi)[i][qp], (*local_v)(dof_indices[i]));
1674 v_norm = std::max(v_norm, norm_weight * hess_u_h.norm());
1675 }
1676 #endif
1677 }
1678 }
1679 }
1680
1681 if (using_hilbert_norm)
1682 {
1683 this->comm().sum(v_norm);
1684 v_norm = std::sqrt(v_norm);
1685 }
1686 else
1687 {
1688 this->comm().max(v_norm);
1689 }
1690
1691 return v_norm;
1692 }
1693
1694
1695
get_info()1696 std::string System::get_info() const
1697 {
1698 std::ostringstream oss;
1699
1700
1701 const std::string & sys_name = this->name();
1702
1703 oss << " System #" << this->number() << ", \"" << sys_name << "\"\n"
1704 << " Type \"" << this->system_type() << "\"\n"
1705 << " Variables=";
1706
1707 for (auto vg : make_range(this->n_variable_groups()))
1708 {
1709 const VariableGroup & vg_description (this->variable_group(vg));
1710
1711 if (vg_description.n_variables() > 1) oss << "{ ";
1712 for (auto vn : make_range(vg_description.n_variables()))
1713 oss << "\"" << vg_description.name(vn) << "\" ";
1714 if (vg_description.n_variables() > 1) oss << "} ";
1715 }
1716
1717 oss << '\n';
1718
1719 oss << " Finite Element Types=";
1720 #ifndef LIBMESH_ENABLE_INFINITE_ELEMENTS
1721 for (auto vg : make_range(this->n_variable_groups()))
1722 oss << "\""
1723 << Utility::enum_to_string<FEFamily>(this->get_dof_map().variable_group(vg).type().family)
1724 << "\" ";
1725 #else
1726 for (auto vg : make_range(this->n_variable_groups()))
1727 {
1728 oss << "\""
1729 << Utility::enum_to_string<FEFamily>(this->get_dof_map().variable_group(vg).type().family)
1730 << "\", \""
1731 << Utility::enum_to_string<FEFamily>(this->get_dof_map().variable_group(vg).type().radial_family)
1732 << "\" ";
1733 }
1734
1735 oss << '\n' << " Infinite Element Mapping=";
1736 for (auto vg : make_range(this->n_variable_groups()))
1737 oss << "\""
1738 << Utility::enum_to_string<InfMapType>(this->get_dof_map().variable_group(vg).type().inf_map)
1739 << "\" ";
1740 #endif
1741
1742 oss << '\n';
1743
1744 oss << " Approximation Orders=";
1745 for (auto vg : make_range(this->n_variable_groups()))
1746 {
1747 #ifndef LIBMESH_ENABLE_INFINITE_ELEMENTS
1748 oss << "\""
1749 << Utility::enum_to_string<Order>(this->get_dof_map().variable_group(vg).type().order)
1750 << "\" ";
1751 #else
1752 oss << "\""
1753 << Utility::enum_to_string<Order>(this->get_dof_map().variable_group(vg).type().order)
1754 << "\", \""
1755 << Utility::enum_to_string<Order>(this->get_dof_map().variable_group(vg).type().radial_order)
1756 << "\" ";
1757 #endif
1758 }
1759
1760 oss << '\n';
1761
1762 oss << " n_dofs()=" << this->n_dofs() << '\n';
1763 oss << " n_local_dofs()=" << this->n_local_dofs() << '\n';
1764 #ifdef LIBMESH_ENABLE_CONSTRAINTS
1765 oss << " n_constrained_dofs()=" << this->n_constrained_dofs() << '\n';
1766 oss << " n_local_constrained_dofs()=" << this->n_local_constrained_dofs() << '\n';
1767 #endif
1768
1769 oss << " " << "n_vectors()=" << this->n_vectors() << '\n';
1770 oss << " " << "n_matrices()=" << this->n_matrices() << '\n';
1771 // oss << " " << "n_additional_matrices()=" << this->n_additional_matrices() << '\n';
1772
1773 oss << this->get_dof_map().get_info();
1774
1775 return oss.str();
1776 }
1777
1778
1779
attach_init_function(void fptr (EquationSystems & es,const std::string & name))1780 void System::attach_init_function (void fptr(EquationSystems & es,
1781 const std::string & name))
1782 {
1783 libmesh_assert(fptr);
1784
1785 if (_init_system_object != nullptr)
1786 {
1787 libmesh_here();
1788 libMesh::out << "WARNING: Cannot specify both initialization function and object!"
1789 << std::endl;
1790
1791 _init_system_object = nullptr;
1792 }
1793
1794 _init_system_function = fptr;
1795 }
1796
1797
1798
attach_init_object(System::Initialization & init_in)1799 void System::attach_init_object (System::Initialization & init_in)
1800 {
1801 if (_init_system_function != nullptr)
1802 {
1803 libmesh_here();
1804 libMesh::out << "WARNING: Cannot specify both initialization object and function!"
1805 << std::endl;
1806
1807 _init_system_function = nullptr;
1808 }
1809
1810 _init_system_object = &init_in;
1811 }
1812
1813
1814
attach_assemble_function(void fptr (EquationSystems & es,const std::string & name))1815 void System::attach_assemble_function (void fptr(EquationSystems & es,
1816 const std::string & name))
1817 {
1818 libmesh_assert(fptr);
1819
1820 if (_assemble_system_object != nullptr)
1821 {
1822 libmesh_here();
1823 libMesh::out << "WARNING: Cannot specify both assembly function and object!"
1824 << std::endl;
1825
1826 _assemble_system_object = nullptr;
1827 }
1828
1829 _assemble_system_function = fptr;
1830 }
1831
1832
1833
attach_assemble_object(System::Assembly & assemble_in)1834 void System::attach_assemble_object (System::Assembly & assemble_in)
1835 {
1836 if (_assemble_system_function != nullptr)
1837 {
1838 libmesh_here();
1839 libMesh::out << "WARNING: Cannot specify both assembly object and function!"
1840 << std::endl;
1841
1842 _assemble_system_function = nullptr;
1843 }
1844
1845 _assemble_system_object = &assemble_in;
1846 }
1847
1848
1849
attach_constraint_function(void fptr (EquationSystems & es,const std::string & name))1850 void System::attach_constraint_function(void fptr(EquationSystems & es,
1851 const std::string & name))
1852 {
1853 libmesh_assert(fptr);
1854
1855 if (_constrain_system_object != nullptr)
1856 {
1857 libmesh_here();
1858 libMesh::out << "WARNING: Cannot specify both constraint function and object!"
1859 << std::endl;
1860
1861 _constrain_system_object = nullptr;
1862 }
1863
1864 _constrain_system_function = fptr;
1865 }
1866
1867
1868
attach_constraint_object(System::Constraint & constrain)1869 void System::attach_constraint_object (System::Constraint & constrain)
1870 {
1871 if (_constrain_system_function != nullptr)
1872 {
1873 libmesh_here();
1874 libMesh::out << "WARNING: Cannot specify both constraint object and function!"
1875 << std::endl;
1876
1877 _constrain_system_function = nullptr;
1878 }
1879
1880 _constrain_system_object = &constrain;
1881 }
1882
1883
1884
attach_QOI_function(void fptr (EquationSystems &,const std::string &,const QoISet &))1885 void System::attach_QOI_function(void fptr(EquationSystems &,
1886 const std::string &,
1887 const QoISet &))
1888 {
1889 libmesh_assert(fptr);
1890
1891 if (_qoi_evaluate_object != nullptr)
1892 {
1893 libmesh_here();
1894 libMesh::out << "WARNING: Cannot specify both QOI function and object!"
1895 << std::endl;
1896
1897 _qoi_evaluate_object = nullptr;
1898 }
1899
1900 _qoi_evaluate_function = fptr;
1901 }
1902
1903
1904
attach_QOI_object(QOI & qoi_in)1905 void System::attach_QOI_object (QOI & qoi_in)
1906 {
1907 if (_qoi_evaluate_function != nullptr)
1908 {
1909 libmesh_here();
1910 libMesh::out << "WARNING: Cannot specify both QOI object and function!"
1911 << std::endl;
1912
1913 _qoi_evaluate_function = nullptr;
1914 }
1915
1916 _qoi_evaluate_object = &qoi_in;
1917 }
1918
1919
1920
attach_QOI_derivative(void fptr (EquationSystems &,const std::string &,const QoISet &,bool,bool))1921 void System::attach_QOI_derivative(void fptr(EquationSystems &, const std::string &,
1922 const QoISet &, bool, bool))
1923 {
1924 libmesh_assert(fptr);
1925
1926 if (_qoi_evaluate_derivative_object != nullptr)
1927 {
1928 libmesh_here();
1929 libMesh::out << "WARNING: Cannot specify both QOI derivative function and object!"
1930 << std::endl;
1931
1932 _qoi_evaluate_derivative_object = nullptr;
1933 }
1934
1935 _qoi_evaluate_derivative_function = fptr;
1936 }
1937
1938
1939
attach_QOI_derivative_object(QOIDerivative & qoi_derivative)1940 void System::attach_QOI_derivative_object (QOIDerivative & qoi_derivative)
1941 {
1942 if (_qoi_evaluate_derivative_function != nullptr)
1943 {
1944 libmesh_here();
1945 libMesh::out << "WARNING: Cannot specify both QOI derivative object and function!"
1946 << std::endl;
1947
1948 _qoi_evaluate_derivative_function = nullptr;
1949 }
1950
1951 _qoi_evaluate_derivative_object = &qoi_derivative;
1952 }
1953
1954
1955
user_initialization()1956 void System::user_initialization ()
1957 {
1958 // Call the user-provided initialization function,
1959 // if it was provided
1960 if (_init_system_function != nullptr)
1961 this->_init_system_function (_equation_systems, this->name());
1962
1963 // ...or the user-provided initialization object.
1964 else if (_init_system_object != nullptr)
1965 this->_init_system_object->initialize();
1966 }
1967
1968
1969
user_assembly()1970 void System::user_assembly ()
1971 {
1972 // Call the user-provided assembly function,
1973 // if it was provided
1974 if (_assemble_system_function != nullptr)
1975 this->_assemble_system_function (_equation_systems, this->name());
1976
1977 // ...or the user-provided assembly object.
1978 else if (_assemble_system_object != nullptr)
1979 this->_assemble_system_object->assemble();
1980 }
1981
1982
1983
user_constrain()1984 void System::user_constrain ()
1985 {
1986 // Call the user-provided constraint function,
1987 // if it was provided
1988 if (_constrain_system_function!= nullptr)
1989 this->_constrain_system_function(_equation_systems, this->name());
1990
1991 // ...or the user-provided constraint object.
1992 else if (_constrain_system_object != nullptr)
1993 this->_constrain_system_object->constrain();
1994 }
1995
1996
1997
user_QOI(const QoISet & qoi_indices)1998 void System::user_QOI (const QoISet & qoi_indices)
1999 {
2000 // Call the user-provided quantity of interest function,
2001 // if it was provided
2002 if (_qoi_evaluate_function != nullptr)
2003 this->_qoi_evaluate_function(_equation_systems, this->name(), qoi_indices);
2004
2005 // ...or the user-provided QOI function object.
2006 else if (_qoi_evaluate_object != nullptr)
2007 this->_qoi_evaluate_object->qoi(qoi_indices);
2008 }
2009
2010
2011
user_QOI_derivative(const QoISet & qoi_indices,bool include_liftfunc,bool apply_constraints)2012 void System::user_QOI_derivative(const QoISet & qoi_indices,
2013 bool include_liftfunc,
2014 bool apply_constraints)
2015 {
2016 // Call the user-provided quantity of interest derivative,
2017 // if it was provided
2018 if (_qoi_evaluate_derivative_function != nullptr)
2019 this->_qoi_evaluate_derivative_function
2020 (_equation_systems, this->name(), qoi_indices, include_liftfunc,
2021 apply_constraints);
2022
2023 // ...or the user-provided QOI derivative function object.
2024 else if (_qoi_evaluate_derivative_object != nullptr)
2025 this->_qoi_evaluate_derivative_object->qoi_derivative
2026 (qoi_indices, include_liftfunc, apply_constraints);
2027 }
2028
2029
2030
point_value(unsigned int var,const Point & p,const bool insist_on_success,const NumericVector<Number> * sol)2031 Number System::point_value(unsigned int var,
2032 const Point & p,
2033 const bool insist_on_success,
2034 const NumericVector<Number> *sol) const
2035 {
2036 // This function must be called on every processor; there's no
2037 // telling where in the partition p falls.
2038 parallel_object_only();
2039
2040 // And every processor had better agree about which point we're
2041 // looking for
2042 #ifndef NDEBUG
2043 libmesh_assert(this->comm().verify(p(0)));
2044 #if LIBMESH_DIM > 1
2045 libmesh_assert(this->comm().verify(p(1)));
2046 #endif
2047 #if LIBMESH_DIM > 2
2048 libmesh_assert(this->comm().verify(p(2)));
2049 #endif
2050 #endif // NDEBUG
2051
2052 // Get a reference to the mesh object associated with the system object that calls this function
2053 const MeshBase & mesh = this->get_mesh();
2054
2055 // Use an existing PointLocator or create a new one
2056 std::unique_ptr<PointLocatorBase> locator_ptr = mesh.sub_point_locator();
2057 PointLocatorBase & locator = *locator_ptr;
2058
2059 if (!insist_on_success || !mesh.is_serial())
2060 locator.enable_out_of_mesh_mode();
2061
2062 // Get a pointer to an element that contains p and allows us to
2063 // evaluate var
2064 const std::set<subdomain_id_type> & raw_subdomains =
2065 this->variable(var).active_subdomains();
2066 const std::set<subdomain_id_type> * implicit_subdomains =
2067 raw_subdomains.empty() ? nullptr : &raw_subdomains;
2068 const Elem * e = locator(p, implicit_subdomains);
2069
2070 Number u = 0;
2071
2072 if (e && this->get_dof_map().is_evaluable(*e, var))
2073 u = point_value(var, p, *e, sol);
2074
2075 // If I have an element containing p, then let's let everyone know
2076 processor_id_type lowest_owner =
2077 (e && (e->processor_id() == this->processor_id())) ?
2078 this->processor_id() : this->n_processors();
2079 this->comm().min(lowest_owner);
2080
2081 // Everybody should get their value from a processor that was able
2082 // to compute it.
2083 // If nobody admits owning the point, we have a problem.
2084 if (lowest_owner != this->n_processors())
2085 this->comm().broadcast(u, lowest_owner);
2086 else
2087 libmesh_assert(!insist_on_success);
2088
2089 return u;
2090 }
2091
point_value(unsigned int var,const Point & p,const Elem & e,const NumericVector<Number> * sol)2092 Number System::point_value(unsigned int var,
2093 const Point & p,
2094 const Elem & e,
2095 const NumericVector<Number> *sol) const
2096 {
2097 // Ensuring that the given point is really in the element is an
2098 // expensive assert, but as long as debugging is turned on we might
2099 // as well try to catch a particularly nasty potential error
2100 libmesh_assert (e.contains_point(p));
2101
2102 if (!sol)
2103 sol = this->current_local_solution.get();
2104
2105 // Get the dof map to get the proper indices for our computation
2106 const DofMap & dof_map = this->get_dof_map();
2107
2108 // Make sure we can evaluate on this element.
2109 libmesh_assert (dof_map.is_evaluable(e, var));
2110
2111 // Need dof_indices for phi[i][j]
2112 std::vector<dof_id_type> dof_indices;
2113
2114 // Fill in the dof_indices for our element
2115 dof_map.dof_indices (&e, dof_indices, var);
2116
2117 // Get the no of dofs associated with this point
2118 const unsigned int num_dofs = cast_int<unsigned int>
2119 (dof_indices.size());
2120
2121 FEType fe_type = dof_map.variable_type(var);
2122
2123 // Map the physical co-ordinates to the master co-ordinates
2124 Point coor = FEMap::inverse_map(e.dim(), &e, p);
2125
2126 // get the shape function value via the FEInterface to also handle the case
2127 // of infinite elements correcly, the shape function is not fe->phi().
2128 FEComputeData fe_data(this->get_equation_systems(), coor);
2129 FEInterface::compute_data(e.dim(), fe_type, &e, fe_data);
2130
2131 // Get ready to accumulate a value
2132 Number u = 0;
2133
2134 for (unsigned int l=0; l<num_dofs; l++)
2135 {
2136 u += fe_data.shape[l] * (*sol)(dof_indices[l]);
2137 }
2138
2139 return u;
2140 }
2141
2142
2143
point_value(unsigned int var,const Point & p,const Elem * e)2144 Number System::point_value(unsigned int var, const Point & p, const Elem * e) const
2145 {
2146 libmesh_assert(e);
2147 return this->point_value(var, p, *e);
2148 }
2149
2150
2151
point_value(unsigned int var,const Point & p,const NumericVector<Number> * sol)2152 Number System::point_value(unsigned int var, const Point & p, const NumericVector<Number> * sol) const
2153 {
2154 return this->point_value(var, p, true, sol);
2155 }
2156
2157
2158
2159
point_gradient(unsigned int var,const Point & p,const bool insist_on_success,const NumericVector<Number> * sol)2160 Gradient System::point_gradient(unsigned int var,
2161 const Point & p,
2162 const bool insist_on_success,
2163 const NumericVector<Number> *sol) const
2164 {
2165 // This function must be called on every processor; there's no
2166 // telling where in the partition p falls.
2167 parallel_object_only();
2168
2169 // And every processor had better agree about which point we're
2170 // looking for
2171 #ifndef NDEBUG
2172 libmesh_assert(this->comm().verify(p(0)));
2173 #if LIBMESH_DIM > 1
2174 libmesh_assert(this->comm().verify(p(1)));
2175 #endif
2176 #if LIBMESH_DIM > 2
2177 libmesh_assert(this->comm().verify(p(2)));
2178 #endif
2179 #endif // NDEBUG
2180
2181 // Get a reference to the mesh object associated with the system object that calls this function
2182 const MeshBase & mesh = this->get_mesh();
2183
2184 // Use an existing PointLocator or create a new one
2185 std::unique_ptr<PointLocatorBase> locator_ptr = mesh.sub_point_locator();
2186 PointLocatorBase & locator = *locator_ptr;
2187
2188 if (!insist_on_success || !mesh.is_serial())
2189 locator.enable_out_of_mesh_mode();
2190
2191 // Get a pointer to an element that contains p and allows us to
2192 // evaluate var
2193 const std::set<subdomain_id_type> & raw_subdomains =
2194 this->variable(var).active_subdomains();
2195 const std::set<subdomain_id_type> * implicit_subdomains =
2196 raw_subdomains.empty() ? nullptr : &raw_subdomains;
2197 const Elem * e = locator(p, implicit_subdomains);
2198
2199 Gradient grad_u;
2200
2201 if (e && this->get_dof_map().is_evaluable(*e, var))
2202 grad_u = point_gradient(var, p, *e, sol);
2203
2204 // If I have an element containing p, then let's let everyone know
2205 processor_id_type lowest_owner =
2206 (e && (e->processor_id() == this->processor_id())) ?
2207 this->processor_id() : this->n_processors();
2208 this->comm().min(lowest_owner);
2209
2210 // Everybody should get their value from a processor that was able
2211 // to compute it.
2212 // If nobody admits owning the point, we may have a problem.
2213 if (lowest_owner != this->n_processors())
2214 this->comm().broadcast(grad_u, lowest_owner);
2215 else
2216 libmesh_assert(!insist_on_success);
2217
2218 return grad_u;
2219 }
2220
2221
point_gradient(unsigned int var,const Point & p,const Elem & e,const NumericVector<Number> * sol)2222 Gradient System::point_gradient(unsigned int var,
2223 const Point & p,
2224 const Elem & e,
2225 const NumericVector<Number> *sol) const
2226 {
2227 // Ensuring that the given point is really in the element is an
2228 // expensive assert, but as long as debugging is turned on we might
2229 // as well try to catch a particularly nasty potential error
2230 libmesh_assert (e.contains_point(p));
2231
2232 if (!sol)
2233 sol = this->current_local_solution.get();
2234
2235 // Get the dof map to get the proper indices for our computation
2236 const DofMap & dof_map = this->get_dof_map();
2237
2238 // write the element dimension into a separate variable.
2239 const unsigned int dim = e.dim();
2240
2241 // Make sure we can evaluate on this element.
2242 libmesh_assert (dof_map.is_evaluable(e, var));
2243
2244 // Need dof_indices for phi[i][j]
2245 std::vector<dof_id_type> dof_indices;
2246
2247 // Fill in the dof_indices for our element
2248 dof_map.dof_indices (&e, dof_indices, var);
2249
2250 // Get the no of dofs associated with this point
2251 const unsigned int num_dofs = cast_int<unsigned int>
2252 (dof_indices.size());
2253
2254 FEType fe_type = dof_map.variable_type(var);
2255
2256 // Map the physical co-ordinates to the master co-ordinates
2257 Point coor = FEMap::inverse_map(dim, &e, p);
2258
2259 // get the shape function value via the FEInterface to also handle the case
2260 // of infinite elements correcly, the shape function is not fe->phi().
2261 FEComputeData fe_data(this->get_equation_systems(), coor);
2262 fe_data.enable_derivative();
2263 FEInterface::compute_data(dim, fe_type, &e, fe_data);
2264
2265 // Get ready to accumulate a gradient
2266 Gradient grad_u;
2267
2268 for (unsigned int l=0; l<num_dofs; l++)
2269 {
2270 // Chartesian coordinates have allways LIBMESH_DIM entries,
2271 // local coordinates have as many coordinates as the element has.
2272 for (std::size_t v=0; v<dim; v++)
2273 for (std::size_t xyz=0; xyz<LIBMESH_DIM; xyz++)
2274 {
2275 // FIXME: this needs better syntax: It is matrix-vector multiplication.
2276 grad_u(xyz) += fe_data.local_transform[v][xyz]
2277 * fe_data.dshape[l](v)
2278 * (*sol)(dof_indices[l]);
2279 }
2280 }
2281
2282 return grad_u;
2283 }
2284
2285
2286
point_gradient(unsigned int var,const Point & p,const Elem * e)2287 Gradient System::point_gradient(unsigned int var, const Point & p, const Elem * e) const
2288 {
2289 libmesh_assert(e);
2290 return this->point_gradient(var, p, *e);
2291 }
2292
2293
2294
point_gradient(unsigned int var,const Point & p,const NumericVector<Number> * sol)2295 Gradient System::point_gradient(unsigned int var, const Point & p, const NumericVector<Number> * sol) const
2296 {
2297 return this->point_gradient(var, p, true, sol);
2298 }
2299
2300
2301
2302 // We can only accumulate a hessian with --enable-second
2303 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
point_hessian(unsigned int var,const Point & p,const bool insist_on_success,const NumericVector<Number> * sol)2304 Tensor System::point_hessian(unsigned int var,
2305 const Point & p,
2306 const bool insist_on_success,
2307 const NumericVector<Number> *sol) const
2308 {
2309 // This function must be called on every processor; there's no
2310 // telling where in the partition p falls.
2311 parallel_object_only();
2312
2313 // And every processor had better agree about which point we're
2314 // looking for
2315 #ifndef NDEBUG
2316 libmesh_assert(this->comm().verify(p(0)));
2317 #if LIBMESH_DIM > 1
2318 libmesh_assert(this->comm().verify(p(1)));
2319 #endif
2320 #if LIBMESH_DIM > 2
2321 libmesh_assert(this->comm().verify(p(2)));
2322 #endif
2323 #endif // NDEBUG
2324
2325 // Get a reference to the mesh object associated with the system object that calls this function
2326 const MeshBase & mesh = this->get_mesh();
2327
2328 // Use an existing PointLocator or create a new one
2329 std::unique_ptr<PointLocatorBase> locator_ptr = mesh.sub_point_locator();
2330 PointLocatorBase & locator = *locator_ptr;
2331
2332 if (!insist_on_success || !mesh.is_serial())
2333 locator.enable_out_of_mesh_mode();
2334
2335 // Get a pointer to an element that contains p and allows us to
2336 // evaluate var
2337 const std::set<subdomain_id_type> & raw_subdomains =
2338 this->variable(var).active_subdomains();
2339 const std::set<subdomain_id_type> * implicit_subdomains =
2340 raw_subdomains.empty() ? nullptr : &raw_subdomains;
2341 const Elem * e = locator(p, implicit_subdomains);
2342
2343 Tensor hess_u;
2344
2345 if (e && this->get_dof_map().is_evaluable(*e, var))
2346 hess_u = point_hessian(var, p, *e, sol);
2347
2348 // If I have an element containing p, then let's let everyone know
2349 processor_id_type lowest_owner =
2350 (e && (e->processor_id() == this->processor_id())) ?
2351 this->processor_id() : this->n_processors();
2352 this->comm().min(lowest_owner);
2353
2354 // Everybody should get their value from a processor that was able
2355 // to compute it.
2356 // If nobody admits owning the point, we may have a problem.
2357 if (lowest_owner != this->n_processors())
2358 this->comm().broadcast(hess_u, lowest_owner);
2359 else
2360 libmesh_assert(!insist_on_success);
2361
2362 return hess_u;
2363 }
2364
point_hessian(unsigned int var,const Point & p,const Elem & e,const NumericVector<Number> * sol)2365 Tensor System::point_hessian(unsigned int var,
2366 const Point & p,
2367 const Elem & e,
2368 const NumericVector<Number> *sol) const
2369 {
2370 // Ensuring that the given point is really in the element is an
2371 // expensive assert, but as long as debugging is turned on we might
2372 // as well try to catch a particularly nasty potential error
2373 libmesh_assert (e.contains_point(p));
2374
2375 if (!sol)
2376 sol = this->current_local_solution.get();
2377
2378 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
2379 if (e.infinite())
2380 libmesh_not_implemented();
2381 #endif
2382
2383 // Get the dof map to get the proper indices for our computation
2384 const DofMap & dof_map = this->get_dof_map();
2385
2386 // Make sure we can evaluate on this element.
2387 libmesh_assert (dof_map.is_evaluable(e, var));
2388
2389 // Need dof_indices for phi[i][j]
2390 std::vector<dof_id_type> dof_indices;
2391
2392 // Fill in the dof_indices for our element
2393 dof_map.dof_indices (&e, dof_indices, var);
2394
2395 // Get the no of dofs associated with this point
2396 const unsigned int num_dofs = cast_int<unsigned int>
2397 (dof_indices.size());
2398
2399 FEType fe_type = dof_map.variable_type(var);
2400
2401 // Build a FE again so we can calculate u(p)
2402 std::unique_ptr<FEBase> fe (FEBase::build(e.dim(), fe_type));
2403
2404 // Map the physical co-ordinates to the master co-ordinates
2405 // Build a vector of point co-ordinates to send to reinit
2406 std::vector<Point> coor(1, FEMap::inverse_map(e.dim(), &e, p));
2407
2408 // Get the values of the shape function derivatives
2409 const std::vector<std::vector<RealTensor>> & d2phi = fe->get_d2phi();
2410
2411 // Reinitialize the element and compute the shape function values at coor
2412 fe->reinit (&e, &coor);
2413
2414 // Get ready to accumulate a hessian
2415 Tensor hess_u;
2416
2417 for (unsigned int l=0; l<num_dofs; l++)
2418 {
2419 hess_u.add_scaled (d2phi[l][0], (*sol)(dof_indices[l]));
2420 }
2421
2422 return hess_u;
2423 }
2424
2425
2426
point_hessian(unsigned int var,const Point & p,const Elem * e)2427 Tensor System::point_hessian(unsigned int var, const Point & p, const Elem * e) const
2428 {
2429 libmesh_assert(e);
2430 return this->point_hessian(var, p, *e);
2431 }
2432
2433
2434
point_hessian(unsigned int var,const Point & p,const NumericVector<Number> * sol)2435 Tensor System::point_hessian(unsigned int var, const Point & p, const NumericVector<Number> * sol) const
2436 {
2437 return this->point_hessian(var, p, true, sol);
2438 }
2439
2440 #else
2441
point_hessian(unsigned int,const Point &,const bool,const NumericVector<Number> *)2442 Tensor System::point_hessian(unsigned int, const Point &, const bool,
2443 const NumericVector<Number> *) const
2444 {
2445 libmesh_error_msg("We can only accumulate a hessian with --enable-second");
2446
2447 // Avoid compiler warnings
2448 return Tensor();
2449 }
2450
point_hessian(unsigned int,const Point &,const Elem &,const NumericVector<Number> *)2451 Tensor System::point_hessian(unsigned int, const Point &, const Elem &,
2452 const NumericVector<Number> *) const
2453 {
2454 libmesh_error_msg("We can only accumulate a hessian with --enable-second");
2455
2456 // Avoid compiler warnings
2457 return Tensor();
2458 }
2459
point_hessian(unsigned int,const Point &,const Elem *)2460 Tensor System::point_hessian(unsigned int, const Point &, const Elem *) const
2461 {
2462 libmesh_error_msg("We can only accumulate a hessian with --enable-second");
2463
2464 // Avoid compiler warnings
2465 return Tensor();
2466 }
2467
point_hessian(unsigned int,const Point &,const NumericVector<Number> *)2468 Tensor System::point_hessian(unsigned int, const Point &, const NumericVector<Number> *) const
2469 {
2470 libmesh_error_msg("We can only accumulate a hessian with --enable-second");
2471
2472 // Avoid compiler warnings
2473 return Tensor();
2474 }
2475
2476 #endif // LIBMESH_ENABLE_SECOND_DERIVATIVES
2477
2478 } // namespace libMesh
2479