libMesh
Loading...
Searching...
No Matches
fem_context.C
Go to the documentation of this file.
1// The libMesh Finite Element Library.
2// Copyright (C) 2002-2026 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3
4// This library is free software; you can redistribute it and/or
5// modify it under the terms of the GNU Lesser General Public
6// License as published by the Free Software Foundation; either
7// version 2.1 of the License, or (at your option) any later version.
8
9// This library is distributed in the hope that it will be useful,
10// but WITHOUT ANY WARRANTY; without even the implied warranty of
11// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12// Lesser General Public License for more details.
13
14// You should have received a copy of the GNU Lesser General Public
15// License along with this library; if not, write to the Free Software
16// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17
18
19
20#include "libmesh/fem_context.h"
21
22#include "libmesh/boundary_info.h"
23#include "libmesh/diff_system.h"
24#include "libmesh/dof_map.h"
25#include "libmesh/elem.h"
26#include "libmesh/fe_base.h"
27#include "libmesh/fe_interface.h"
28#include "libmesh/libmesh_logging.h"
29#include "libmesh/mesh_base.h"
30#include "libmesh/numeric_vector.h"
31#include "libmesh/quadrature.h"
32#include "libmesh/system.h"
33#include "libmesh/time_solver.h"
34#include "libmesh/unsteady_solver.h" // For euler_residual
35
36namespace libMesh
37{
38
39FEMContext::FEMContext (const System & sys,
40 const std::vector<unsigned int> * active_vars,
41 bool allocate_local_matrices)
42 : FEMContext(sys, sys.extra_quadrature_order, active_vars,
43 allocate_local_matrices)
44{
45 init_internal_data(sys);
46}
47
48FEMContext::FEMContext (const System & sys,
49 int extra_quadrature_order,
50 const std::vector<unsigned int> * active_vars,
51 bool allocate_local_matrices)
52 : DiffContext(sys, allocate_local_matrices),
53 side(0), edge(0),
54 _mesh_sys(nullptr),
55 _mesh_x_var(0),
56 _mesh_y_var(0),
57 _mesh_z_var(0),
58 _atype(CURRENT),
59 _custom_solution(nullptr),
60 _boundary_info(sys.get_mesh().get_boundary_info()),
61 _elem(nullptr),
62 _dim(cast_int<unsigned char>(sys.get_mesh().mesh_dimension())),
63 _elem_dim(0), /* This will be reset in set_elem(). */
64 _elem_dims(sys.get_mesh().elem_dimensions()),
65 _element_qrule(4),
66 _side_qrule(4),
67 _extra_quadrature_order(extra_quadrature_order)
68{
69 if (active_vars)
70 {
71 libmesh_assert(!active_vars->empty());
72 auto vars_copy =
73 std::make_unique<std::vector<unsigned int>>(*active_vars);
74
75 // We want to do quick binary_search later
76 std::sort(vars_copy->begin(), vars_copy->end());
77
78 _active_vars = std::move(vars_copy);
79 }
80
81 init_internal_data(sys);
82}
83
84
85
86FEType FEMContext::find_hardest_fe_type()
87{
88 const System & sys = this->get_system();
89 FEType hardest_fe_type =
90 sys.variable_type(_active_vars ?
91 (*_active_vars)[0] : 0);
92
93 auto check_var = [&hardest_fe_type, &sys](unsigned int v)
94 {
95 FEType fe_type = sys.variable_type(v);
96
97 // Make sure we find a non-SCALAR FE family, even in the case
98 // where the first variable(s) weren't
99 if (hardest_fe_type.family == SCALAR)
100 {
101 hardest_fe_type.family = fe_type.family;
102 hardest_fe_type.order = fe_type.order;
103 }
104
105 // FIXME - we don't yet handle mixed finite elements from
106 // different families which require different quadrature rules
107 // libmesh_assert_equal_to (fe_type.family, hardest_fe_type.family);
108
109 // We need to detect SCALAR's so we can prepare FE objects for
110 // them, and so we don't mistake high order scalars as a reason
111 // to crank up the quadrature order on other types.
112 if (fe_type.family != SCALAR && fe_type.order > hardest_fe_type.order)
113 hardest_fe_type = fe_type;
114 };
115
116 if (_active_vars)
117 for (auto v : *_active_vars)
118 check_var(v);
119 else
120 for (auto v : make_range(sys.n_vars()))
121 check_var(v);
122
123 return hardest_fe_type;
124}
125
126
127void FEMContext::attach_quadrature_rules()
128{
129 const System & sys = this->get_system();
130
131 auto attach_rules = [this, &sys](unsigned int v)
132 {
133 for (const auto & dim : _elem_dims)
134 {
135 FEType fe_type = sys.variable_type(v);
136
137 _element_fe[dim][fe_type]->attach_quadrature_rule(_element_qrule[dim].get());
138 if (dim)
139 _side_fe[dim][fe_type]->attach_quadrature_rule(_side_qrule[dim].get());
140 if (dim == 3)
141 _edge_fe[fe_type]->attach_quadrature_rule(_edge_qrule.get());
142 };
143 };
144
145 if (_active_vars)
146 for (auto v : *_active_vars)
147 attach_rules(v);
148 else
149 for (auto v : make_range(sys.n_vars()))
150 attach_rules(v);
151}
152
153
154
155void FEMContext::use_default_quadrature_rules(int extra_quadrature_order)
156{
157 _extra_quadrature_order = extra_quadrature_order;
158
159 FEType hardest_fe_type = this->find_hardest_fe_type();
160
161 for (const auto & dim : _elem_dims)
162 {
163 // Create an adequate quadrature rule
164 _element_qrule[dim] =
165 hardest_fe_type.default_quadrature_rule(dim, _extra_quadrature_order);
166 if (dim)
167 _side_qrule[dim] =
168 hardest_fe_type.default_quadrature_rule(dim-1, _extra_quadrature_order);
169 if (dim == 3)
170 _edge_qrule = hardest_fe_type.default_quadrature_rule(1, _extra_quadrature_order);
171 }
172
173 this->attach_quadrature_rules();
174}
175
176
177void FEMContext::use_unweighted_quadrature_rules(int extra_quadrature_order)
178{
179 _extra_quadrature_order = extra_quadrature_order;
180
181 FEType hardest_fe_type = this->find_hardest_fe_type();
182
183 for (const auto & dim : _elem_dims)
184 {
185 // Create an adequate quadrature rule
186 _element_qrule[dim] =
187 hardest_fe_type.unweighted_quadrature_rule(dim, _extra_quadrature_order);
188 _side_qrule[dim] =
189 hardest_fe_type.unweighted_quadrature_rule(dim-1, _extra_quadrature_order);
190 if (dim == 3)
191 _edge_qrule = hardest_fe_type.unweighted_quadrature_rule(1, _extra_quadrature_order);
192 }
193
194 this->attach_quadrature_rules();
195}
196
197
198void FEMContext::init_internal_data(const System & sys)
199{
200 // Reserve space for the FEAbstract and QBase objects for each
201 // element dimension possibility (0,1,2,3)
202
203 // Note: we would simply resize() all four of these vectors, but
204 // some compilers (ICC 19, MSVC) generate a diagnostic about copying
205 // a std::unique_ptr in this case, so the following two lines are a
206 // workaround.
207 _element_fe = std::vector<std::map<FEType, std::unique_ptr<FEAbstract>>>(4);
208 _side_fe = std::vector<std::map<FEType, std::unique_ptr<FEAbstract>>>(4);
209 _element_fe_var.resize(4);
210 _side_fe_var.resize(4);
211
212 // We need to know which of our variables has the hardest
213 // shape functions to numerically integrate.
214
215 unsigned int nv = sys.n_vars();
216 libmesh_assert (nv);
217
218 bool have_scalar = false;
219
220 if (_active_vars)
221 {
222 for (auto v : *_active_vars)
223 if (sys.variable_type(v).family == SCALAR)
224 {
225 have_scalar = true;
226 break;
227 }
228 }
229 else
230 {
231 for (auto v : make_range(sys.n_vars()))
232 if (sys.variable_type(v).family == SCALAR)
233 {
234 have_scalar = true;
235 break;
236 }
237 }
238
239 if (have_scalar)
240 // SCALAR FEs have dimension 0 by assumption
241 _elem_dims.insert(0);
242
243 auto build_var_fe = [this, &sys](unsigned int dim,
244 unsigned int i)
245 {
246 FEType fe_type = sys.variable_type(i);
247 const bool add_p_level = fe_type.p_refinement;
248
249 auto & element_fe = _element_fe[dim][fe_type];
250 auto & side_fe = _side_fe[dim][fe_type];
251 if (!element_fe)
252 {
253 element_fe = FEAbstract::build(dim, fe_type);
254 element_fe->add_p_level_in_reinit(add_p_level);
255 side_fe = FEAbstract::build(dim, fe_type);
256 side_fe->add_p_level_in_reinit(add_p_level);
257
258 if (dim == 3)
259 {
260 auto & edge_fe = _edge_fe[fe_type];
261 edge_fe = FEAbstract::build(dim, fe_type);
262 edge_fe->add_p_level_in_reinit(add_p_level);
263 }
264 }
265
266 _element_fe_var[dim][i] = element_fe.get();
267 _side_fe_var[dim][i] = side_fe.get();
268 if ((dim) == 3)
269 _edge_fe_var[i] = _edge_fe[fe_type].get();
270 };
271
272 for (const auto & dim : _elem_dims)
273 {
274 // Create finite element objects
275 _element_fe_var[dim].resize(nv);
276 _side_fe_var[dim].resize(nv);
277 if (dim == 3)
278 _edge_fe_var.resize(nv);
279
280 if (_active_vars)
281 for (auto v : *_active_vars)
282 build_var_fe(dim, v);
283 else
284 for (auto v : make_range(nv))
285 build_var_fe(dim, v);
286 }
287
288 this->use_default_quadrature_rules(_extra_quadrature_order);
289}
290
294
295
296
301
302
303
304void FEMContext::side_boundary_ids(std::vector<boundary_id_type> & vec_to_fill) const
305{
306 _boundary_info.boundary_ids(&(this->get_elem()), side, vec_to_fill);
307}
308
309
310
311template<typename OutputType,
312 typename FEMContext::FENeeded<OutputType>::value_getter fe_getter,
313 FEMContext::diff_subsolution_getter subsolution_getter>
314void FEMContext::some_value(unsigned int var, unsigned int qp, OutputType & u) const
315{
316 // Get local-to-global dof index lookup
317 const unsigned int n_dofs = cast_int<unsigned int>
318 (this->get_dof_indices(var).size());
319
320 // Get current local coefficients
321 const DenseSubVector<Number> & coef = (this->*subsolution_getter)(var);
322 libmesh_assert_equal_to(coef.size(), n_dofs);
323
324 // Get finite element object
325 typename FENeeded<OutputType>::value_base * fe = nullptr;
326 (this->*fe_getter)( var, fe, this->get_elem_dim() );
327
328 // Get shape function values at quadrature point
329 const std::vector<std::vector
330 <typename FENeeded<OutputType>::value_shape>> & phi = fe->get_phi();
331 libmesh_assert_equal_to(phi.size(), n_dofs);
332
333 // Accumulate solution value
334 u = 0.;
335
336 for (unsigned int l=0; l != n_dofs; l++)
337 {
338 libmesh_assert_less(qp, phi[l].size());
339 u += phi[l][qp] * coef(l);
340 }
341}
342
343
344
345template<typename OutputType,
346 typename FEMContext::FENeeded<OutputType>::grad_getter fe_getter,
347 FEMContext::diff_subsolution_getter subsolution_getter>
348void FEMContext::some_gradient(unsigned int var, unsigned int qp, OutputType & du) const
349{
350 // Get local-to-global dof index lookup
351 const unsigned int n_dofs = cast_int<unsigned int>
352 (this->get_dof_indices(var).size());
353
354 // Get current local coefficients
355 const DenseSubVector<Number> & coef = (this->*subsolution_getter)(var);
356
357 // Get finite element object
358 typename FENeeded<OutputType>::grad_base * fe = nullptr;
359 (this->*fe_getter)( var, fe, this->get_elem_dim() );
360
361 // Get shape function values at quadrature point
362 const std::vector<std::vector
363 <typename FENeeded<OutputType>::grad_base::OutputGradient>>
364 & dphi = fe->get_dphi();
365
366 // Accumulate solution derivatives
367 du = 0;
368
369 for (unsigned int l=0; l != n_dofs; l++)
370 du.add_scaled(dphi[l][qp], coef(l));
371
372 return;
373}
374
375
376
377#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
378template<typename OutputType,
379 typename FEMContext::FENeeded<OutputType>::hess_getter fe_getter,
380 FEMContext::diff_subsolution_getter subsolution_getter>
381void FEMContext::some_hessian(unsigned int var, unsigned int qp, OutputType & d2u) const
382{
383 // Get local-to-global dof index lookup
384 const unsigned int n_dofs = cast_int<unsigned int>
385 (this->get_dof_indices(var).size());
386
387 // Get current local coefficients
388 const DenseSubVector<Number> & coef = (this->*subsolution_getter)(var);
389
390 // Get finite element object
391 typename FENeeded<OutputType>::hess_base * fe = nullptr;
392 (this->*fe_getter)( var, fe, this->get_elem_dim() );
393
394 // Get shape function values at quadrature point
395 const std::vector<std::vector
396 <typename FENeeded<OutputType>::hess_base::OutputTensor>>
397 & d2phi = fe->get_d2phi();
398
399 // Accumulate solution second derivatives
400 d2u = 0.0;
401
402 for (unsigned int l=0; l != n_dofs; l++)
403 d2u.add_scaled(d2phi[l][qp], coef(l));
404
405 return;
406}
407#endif
408
409
410
411Number FEMContext::interior_value(unsigned int var, unsigned int qp) const
412{
413 Number u;
414
415 this->interior_value( var, qp, u );
416
417 return u;
418}
419
420template<typename OutputType>
421void FEMContext::interior_value(unsigned int var, unsigned int qp,
422 OutputType & u) const
423{
424 this->some_value<OutputType,
425 &FEMContext::get_element_fe<typename TensorTools::MakeReal<OutputType>::type>,
426 &DiffContext::get_elem_solution>(var, qp, u);
427}
428
429
430template<typename OutputType>
431void FEMContext::interior_values (unsigned int var,
432 const NumericVector<Number> & _system_vector,
433 std::vector<OutputType> & u_vals) const
434{
435 typedef typename TensorTools::MakeReal<OutputType>::type OutputShape;
436
437 // Get local-to-global dof index lookup
438 const unsigned int n_dofs = cast_int<unsigned int>
439 (this->get_dof_indices(var).size());
440
441 // Get current local coefficients
442 const DenseSubVector<Number> & coef = get_localized_subvector(_system_vector, var);
443
444 // Get the finite element object
445 FEGenericBase<OutputShape> * fe = nullptr;
446 this->get_element_fe<OutputShape>( var, fe, this->get_elem_dim() );
447
448 // Get shape function values at quadrature point
449 const std::vector<std::vector<OutputShape>> & phi = fe->get_phi();
450
451 // Loop over all the q_points on this element
452 for (auto qp : index_range(u_vals))
453 {
454 OutputType & u = u_vals[qp];
455
456 // Compute the value at this q_point
457 u = 0.;
458
459 for (unsigned int l=0; l != n_dofs; l++)
460 u += phi[l][qp] * coef(l);
461 }
462
463 return;
464}
465
466Gradient FEMContext::interior_gradient(unsigned int var,
467 unsigned int qp) const
468{
469 Gradient du;
470
471 this->interior_gradient( var, qp, du );
472
473 return du;
474}
475
476
477
478template<typename OutputType>
479void FEMContext::interior_gradient(unsigned int var,
480 unsigned int qp,
481 OutputType & du) const
482{
483 this->some_gradient<OutputType,
484 &FEMContext::get_element_fe<typename TensorTools::MakeReal
485 <typename TensorTools::DecrementRank
486 <OutputType>::type>::type>,
487 &DiffContext::get_elem_solution>(var, qp, du);
488}
489
490
491
492template<typename OutputType>
493void FEMContext::interior_gradients(unsigned int var,
494 const NumericVector<Number> & _system_vector,
495 std::vector<OutputType> & du_vals) const
496{
497 typedef typename TensorTools::MakeReal
498 <typename TensorTools::DecrementRank<OutputType>::type>::type
499 OutputShape;
500
501 // Get local-to-global dof index lookup
502 const unsigned int n_dofs = cast_int<unsigned int>
503 (this->get_dof_indices(var).size());
504
505 // Get current local coefficients
506 const DenseSubVector<Number> & coef = get_localized_subvector(_system_vector, var);
507
508 // Get finite element object
509 FEGenericBase<OutputShape> * fe = nullptr;
510 this->get_element_fe<OutputShape>( var, fe, this->get_elem_dim() );
511
512 // Get shape function values at quadrature point
513 const std::vector<std::vector<typename FEGenericBase<OutputShape>::OutputGradient>> & dphi = fe->get_dphi();
514
515 // Loop over all the q_points in this finite element
516 for (auto qp : index_range(du_vals))
517 {
518 OutputType & du = du_vals[qp];
519
520 // Compute the gradient at this q_point
521 du = 0;
522
523 for (unsigned int l=0; l != n_dofs; l++)
524 du.add_scaled(dphi[l][qp], coef(l));
525 }
526
527 return;
528}
529
530#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
531Tensor FEMContext::interior_hessian(unsigned int var, unsigned int qp) const
532{
533 Tensor d2u;
534
535 this->interior_hessian( var, qp, d2u );
536
537 return d2u;
538}
539
540template<typename OutputType>
541void FEMContext::interior_hessian(unsigned int var, unsigned int qp,
542 OutputType & d2u) const
543{
544 this->some_hessian<OutputType,
545 &FEMContext::get_element_fe
546 <typename TensorTools::MakeReal
547 <typename TensorTools::DecrementRank
548 <typename TensorTools::DecrementRank
549 <OutputType>::type>::type>::type>,
550 &DiffContext::get_elem_solution>(var, qp, d2u);
551}
552
553
554template<typename OutputType>
555void FEMContext::interior_hessians(unsigned int var,
556 const NumericVector<Number> & _system_vector,
557 std::vector<OutputType> & d2u_vals) const
558{
559 typedef typename TensorTools::DecrementRank<OutputType>::type Rank1Decrement;
560 typedef typename TensorTools::DecrementRank<Rank1Decrement>::type Rank2Decrement;
561 typedef typename TensorTools::MakeReal<Rank2Decrement>::type OutputShape;
562
563 // Get local-to-global dof index lookup
564 const unsigned int n_dofs = cast_int<unsigned int>
565 (this->get_dof_indices(var).size());
566
567 // Get current local coefficients
568 const DenseSubVector<Number> & coef = get_localized_subvector(_system_vector, var);
569
570 // Get finite element object
571 FEGenericBase<OutputShape> * fe = nullptr;
572 this->get_element_fe<OutputShape>( var, fe, this->get_elem_dim() );
573
574 // Get shape function values at quadrature point
575 const std::vector<std::vector<typename FEGenericBase<OutputShape>::OutputTensor>> & d2phi = fe->get_d2phi();
576
577 // Loop over all the q_points in this finite element
578 for (auto qp : index_range(d2u_vals))
579 {
580 OutputType & d2u = d2u_vals[qp];
581
582 // Compute the gradient at this q_point
583 d2u = 0;
584
585 for (unsigned int l=0; l != n_dofs; l++)
586 d2u.add_scaled(d2phi[l][qp], coef(l));
587 }
588
589 return;
590}
591
592
593#endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
594
595
596template<typename OutputType>
597void FEMContext::interior_curl(unsigned int var, unsigned int qp,
598 OutputType & curl_u) const
599{
600 typedef typename TensorTools::MakeReal<OutputType>::type OutputShape;
601
602 // Get local-to-global dof index lookup
603 const unsigned int n_dofs = cast_int<unsigned int>
604 (this->get_dof_indices(var).size());
605
606 // Get current local coefficients
607 libmesh_assert_greater (this->_elem_subsolutions.size(), var);
608 const DenseSubVector<Number> & coef = this->get_elem_solution(var);
609
610 // Get finite element object
611 FEGenericBase<OutputShape> * fe = nullptr;
612 this->get_element_fe<OutputShape>( var, fe, this->get_elem_dim() );
613
614 // Get shape function values at quadrature point
615 const std::vector<std::vector<typename FEGenericBase<OutputShape>::OutputShape>> & curl_phi = fe->get_curl_phi();
616
617 // Accumulate solution curl
618 curl_u = 0.;
619
620 for (unsigned int l=0; l != n_dofs; l++)
621 curl_u.add_scaled(curl_phi[l][qp], coef(l));
622
623 return;
624}
625
626
627template<typename OutputType>
628void FEMContext::interior_div(unsigned int var, unsigned int qp,
629 OutputType & div_u) const
630{
631 typedef typename
632 TensorTools::IncrementRank
633 <typename TensorTools::MakeReal<OutputType>::type>::type OutputShape;
634
635 // Get local-to-global dof index lookup
636 const unsigned int n_dofs = cast_int<unsigned int>
637 (this->get_dof_indices(var).size());
638
639 // Get current local coefficients
640 libmesh_assert_greater (this->_elem_subsolutions.size(), var);
641 const DenseSubVector<Number> & coef = this->get_elem_solution(var);
642
643 // Get finite element object
644 FEGenericBase<OutputShape> * fe = nullptr;
645 this->get_element_fe<OutputShape>( var, fe, this->get_elem_dim() );
646
647 // Get shape function values at quadrature point
648 const std::vector<std::vector<typename FEGenericBase<OutputShape>::OutputDivergence>> & div_phi = fe->get_div_phi();
649
650 // Accumulate solution curl
651 div_u = 0.;
652
653 for (unsigned int l=0; l != n_dofs; l++)
654 div_u += div_phi[l][qp] * coef(l);
655
656 return;
657}
658
659
660Number FEMContext::side_value(unsigned int var,
661 unsigned int qp) const
662{
663 Number u = 0.;
664
665 this->side_value( var, qp, u );
666
667 return u;
668}
669
670
671template<typename OutputType>
672void FEMContext::side_value(unsigned int var,
673 unsigned int qp,
674 OutputType & u) const
675{
676 this->some_value<OutputType,
677 &FEMContext::get_side_fe<typename TensorTools::MakeReal<OutputType>::type>,
678 &DiffContext::get_elem_solution>(var, qp, u);
679}
680
681
682template<typename OutputType>
683void FEMContext::side_values(unsigned int var,
684 const NumericVector<Number> & _system_vector,
685 std::vector<OutputType> & u_vals) const
686{
687 typedef typename TensorTools::MakeReal<OutputType>::type OutputShape;
688
689 // Get local-to-global dof index lookup
690 const unsigned int n_dofs = cast_int<unsigned int>
691 (this->get_dof_indices(var).size());
692
693 // Get current local coefficients
694 const DenseSubVector<Number> & coef = get_localized_subvector(_system_vector, var);
695
696 // Get the finite element object
697 FEGenericBase<OutputShape> * the_side_fe = nullptr;
698 this->get_side_fe<OutputShape>( var, the_side_fe, this->get_elem_dim() );
699
700 // Get shape function values at quadrature point
701 const std::vector<std::vector<OutputShape>> & phi = the_side_fe->get_phi();
702
703 // Loop over all the q_points on this element
704 for (auto qp : index_range(u_vals))
705 {
706 OutputType & u = u_vals[qp];
707
708 // Compute the value at this q_point
709 u = 0.;
710
711 for (unsigned int l=0; l != n_dofs; l++)
712 u += phi[l][qp] * coef(l);
713 }
714
715 return;
716}
717
718Gradient FEMContext::side_gradient(unsigned int var, unsigned int qp) const
719{
720 Gradient du;
721
722 this->side_gradient( var, qp, du );
723
724 return du;
725}
726
727
728template<typename OutputType>
729void FEMContext::side_gradient(unsigned int var, unsigned int qp,
730 OutputType & du) const
731{
732 typedef typename TensorTools::MakeReal
733 <typename TensorTools::DecrementRank<OutputType>::type>::type
734 OutputShape;
735
736 // Get local-to-global dof index lookup
737 const unsigned int n_dofs = cast_int<unsigned int>
738 (this->get_dof_indices(var).size());
739
740 // Get current local coefficients
741 libmesh_assert_greater (this->_elem_subsolutions.size(), var);
742 const DenseSubVector<Number> & coef = this->get_elem_solution(var);
743
744 // Get finite element object
745 FEGenericBase<OutputShape> * the_side_fe = nullptr;
746 this->get_side_fe<OutputShape>( var, the_side_fe, this->get_elem_dim() );
747
748 // Get shape function values at quadrature point
749 const std::vector<std::vector<typename FEGenericBase<OutputShape>::OutputGradient>> & dphi = the_side_fe->get_dphi();
750
751 // Accumulate solution derivatives
752 du = 0.;
753
754 for (unsigned int l=0; l != n_dofs; l++)
755 du.add_scaled(dphi[l][qp], coef(l));
756
757 return;
758}
759
760
761
762template<typename OutputType>
763void FEMContext::side_gradients(unsigned int var,
764 const NumericVector<Number> & _system_vector,
765 std::vector<OutputType> & du_vals) const
766{
767 typedef typename TensorTools::MakeReal
768 <typename TensorTools::DecrementRank<OutputType>::type>::type
769 OutputShape;
770
771 // Get local-to-global dof index lookup
772 const unsigned int n_dofs = cast_int<unsigned int>
773 (this->get_dof_indices(var).size());
774
775 // Get current local coefficients
776 const DenseSubVector<Number> & coef = get_localized_subvector(_system_vector, var);
777
778 // Get finite element object
779 FEGenericBase<OutputShape> * the_side_fe = nullptr;
780 this->get_side_fe<OutputShape>( var, the_side_fe, this->get_elem_dim() );
781
782 // Get shape function values at quadrature point
783 const std::vector<std::vector<typename FEGenericBase<OutputShape>::OutputGradient>> & dphi = the_side_fe->get_dphi();
784
785 // Loop over all the q_points in this finite element
786 for (auto qp : index_range(du_vals))
787 {
788 OutputType & du = du_vals[qp];
789
790 du = 0;
791
792 // Compute the gradient at this q_point
793 for (unsigned int l=0; l != n_dofs; l++)
794 du.add_scaled(dphi[l][qp], coef(l));
795 }
796
797 return;
798}
799
800#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
801Tensor FEMContext::side_hessian(unsigned int var,
802 unsigned int qp) const
803{
804 Tensor d2u;
805
806 this->side_hessian( var, qp, d2u );
807
808 return d2u;
809}
810
811
812
813template<typename OutputType>
814void FEMContext::side_hessian(unsigned int var,
815 unsigned int qp,
816 OutputType & d2u) const
817{
818 this->some_hessian<OutputType,
819 &FEMContext::get_side_fe
820 <typename TensorTools::MakeReal
821 <typename TensorTools::DecrementRank
822 <typename TensorTools::DecrementRank
823 <OutputType>::type>::type>::type>,
824 &DiffContext::get_elem_solution>(var, qp, d2u);
825}
826
827
828
829template<typename OutputType>
830void FEMContext::side_hessians(unsigned int var,
831 const NumericVector<Number> & _system_vector,
832 std::vector<OutputType> & d2u_vals) const
833{
834 typedef typename TensorTools::DecrementRank<OutputType>::type Rank1Decrement;
835 typedef typename TensorTools::DecrementRank<Rank1Decrement>::type Rank2Decrement;
836 typedef typename TensorTools::MakeReal<Rank2Decrement>::type OutputShape;
837
838 // Get local-to-global dof index lookup
839 const unsigned int n_dofs = cast_int<unsigned int>
840 (this->get_dof_indices(var).size());
841
842 // Get current local coefficients
843 const DenseSubVector<Number> & coef = get_localized_subvector(_system_vector, var);
844
845 // Get finite element object
846 FEGenericBase<OutputShape> * the_side_fe = nullptr;
847 this->get_side_fe<OutputShape>( var, the_side_fe, this->get_elem_dim() );
848
849 // Get shape function values at quadrature point
850 const std::vector<std::vector<typename FEGenericBase<OutputShape>::OutputTensor>> & d2phi = the_side_fe->get_d2phi();
851
852 // Loop over all the q_points in this finite element
853 for (auto qp : index_range(d2u_vals))
854 {
855 OutputType & d2u = d2u_vals[qp];
856
857 // Compute the gradient at this q_point
858 d2u = 0;
859
860 for (unsigned int l=0; l != n_dofs; l++)
861 d2u.add_scaled(d2phi[l][qp], coef(l));
862 }
863
864 return;
865}
866
867
868
869#endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
870
871
872
873Number FEMContext::point_value(unsigned int var, const Point & p) const
874{
875 Number u = 0.;
876
877 this->point_value( var, p, u );
878
879 return u;
880}
881
882template<typename OutputType>
883void FEMContext::point_value(unsigned int var,
884 const Point & p,
885 OutputType & u,
886 const Real tolerance) const
887{
888 typedef typename TensorTools::MakeReal<OutputType>::type OutputShape;
889
890 // Get local-to-global dof index lookup
891 const unsigned int n_dofs = cast_int<unsigned int>
892 (this->get_dof_indices(var).size());
893
894 // Get current local coefficients
895 libmesh_assert_greater (this->_elem_subsolutions.size(), var);
896 const DenseSubVector<Number> & coef = this->get_elem_solution(var);
897
898 // Get finite element object
899 FEGenericBase<OutputShape> * fe = nullptr;
900 this->get_element_fe<OutputShape>( var, fe, this->get_elem_dim() );
901
902 // Build a FE for calculating u(p)
903 FEGenericBase<OutputShape> * fe_new =
904 this->build_new_fe( fe, p, tolerance, 0 );
905
906 // Get the values of the shape function derivatives
907 const std::vector<std::vector<OutputShape>> & phi = fe_new->get_phi();
908
909 u = 0.;
910
911 for (unsigned int l=0; l != n_dofs; l++)
912 u += phi[l][0] * coef(l);
913
914 return;
915}
916
917
918
919Gradient FEMContext::point_gradient(unsigned int var, const Point & p) const
920{
921 Gradient grad_u;
922
923 this->point_gradient( var, p, grad_u );
924
925 return grad_u;
926}
927
928
929
930template<typename OutputType>
931void FEMContext::point_gradient(unsigned int var,
932 const Point & p,
933 OutputType & grad_u,
934 const Real tolerance) const
935{
936 typedef typename TensorTools::MakeReal
937 <typename TensorTools::DecrementRank<OutputType>::type>::type
938 OutputShape;
939
940 // Get local-to-global dof index lookup
941 const unsigned int n_dofs = cast_int<unsigned int>
942 (this->get_dof_indices(var).size());
943
944 // Get current local coefficients
945 libmesh_assert_greater (this->_elem_subsolutions.size(), var);
946 const DenseSubVector<Number> & coef = this->get_elem_solution(var);
947
948 // Get finite element object
949 FEGenericBase<OutputShape> * fe = nullptr;
950 this->get_element_fe<OutputShape>( var, fe, this->get_elem_dim() );
951
952 // Build a FE for calculating u(p)
953 FEGenericBase<OutputShape> * fe_new =
954 this->build_new_fe( fe, p, tolerance, 1 );
955
956 // Get the values of the shape function derivatives
957 const std::vector<std::vector<typename FEGenericBase<OutputShape>::OutputGradient>> & dphi = fe_new->get_dphi();
958
959 grad_u = 0.0;
960
961 for (unsigned int l=0; l != n_dofs; l++)
962 grad_u.add_scaled(dphi[l][0], coef(l));
963
964 return;
965}
966
967
968
969#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
970
971Tensor FEMContext::point_hessian(unsigned int var, const Point & p) const
972{
973 Tensor hess_u;
974
975 this->point_hessian( var, p, hess_u );
976
977 return hess_u;
978}
979
980
981template<typename OutputType>
982void FEMContext::point_hessian(unsigned int var,
983 const Point & p,
984 OutputType & hess_u,
985 const Real tolerance) const
986{
987 typedef typename TensorTools::DecrementRank<OutputType>::type Rank1Decrement;
988 typedef typename TensorTools::DecrementRank<Rank1Decrement>::type Rank2Decrement;
989 typedef typename TensorTools::MakeReal<Rank2Decrement>::type OutputShape;
990
991 // Get local-to-global dof index lookup
992 const unsigned int n_dofs = cast_int<unsigned int>
993 (this->get_dof_indices(var).size());
994
995 // Get current local coefficients
996 libmesh_assert_greater (this->_elem_subsolutions.size(), var);
997 const DenseSubVector<Number> & coef = this->get_elem_solution(var);
998
999 // Get finite element object
1000 FEGenericBase<OutputShape> * fe = nullptr;
1001 this->get_element_fe<OutputShape>( var, fe, this->get_elem_dim() );
1002
1003 // Build a FE for calculating u(p)
1004 FEGenericBase<OutputShape> * fe_new =
1005 this->build_new_fe( fe, p, tolerance, 2 );
1006
1007 // Get the values of the shape function derivatives
1008 const std::vector<std::vector<typename FEGenericBase<OutputShape>::OutputTensor>> & d2phi = fe_new->get_d2phi();
1009
1010 hess_u = 0.0;
1011
1012 for (unsigned int l=0; l != n_dofs; l++)
1013 hess_u.add_scaled(d2phi[l][0], coef(l));
1014
1015 return;
1016}
1017
1018#endif // LIBMESH_ENABLE_SECOND_DERIVATIVES
1019
1020
1021template<typename OutputType>
1022void FEMContext::point_curl(unsigned int var,
1023 const Point & p,
1024 OutputType & curl_u,
1025 const Real tolerance) const
1026{
1027 typedef typename TensorTools::MakeReal<OutputType>::type OutputShape;
1028
1029 // Get local-to-global dof index lookup
1030 const unsigned int n_dofs = cast_int<unsigned int>
1031 (this->get_dof_indices(var).size());
1032
1033 // Get current local coefficients
1034 libmesh_assert_greater (this->_elem_subsolutions.size(), var);
1035 const DenseSubVector<Number> & coef = this->get_elem_solution(var);
1036
1037 // Get finite element object
1038 FEGenericBase<OutputShape> * fe = nullptr;
1039 this->get_element_fe<OutputShape>( var, fe, this->get_elem_dim() );
1040
1041 // Build a FE for calculating u(p)
1042 FEGenericBase<OutputShape> * fe_new =
1043 this->build_new_fe( fe, p, tolerance, 3 );
1044
1045 // Get the values of the shape function derivatives
1046 const std::vector<std::vector<typename FEGenericBase<OutputShape>::OutputShape>> & curl_phi = fe_new->get_curl_phi();
1047
1048 curl_u = 0.0;
1049
1050 for (unsigned int l=0; l != n_dofs; l++)
1051 curl_u.add_scaled(curl_phi[l][0], coef(l));
1052
1053 return;
1054}
1055
1056
1057
1058Number FEMContext::fixed_interior_value(unsigned int var, unsigned int qp) const
1059{
1060 Number u = 0.;
1061
1062 this->fixed_interior_value( var, qp, u );
1063
1064 return u;
1065}
1066
1067
1068
1069template<typename OutputType>
1070void FEMContext::fixed_interior_value(unsigned int var, unsigned int qp,
1071 OutputType & u) const
1072{
1073 this->some_value<OutputType,
1074 &FEMContext::get_element_fe
1075 <typename TensorTools::MakeReal<OutputType>::type>,
1076 &DiffContext::get_elem_fixed_solution>(var, qp, u);
1077}
1078
1079
1080
1081Gradient FEMContext::fixed_interior_gradient(unsigned int var, unsigned int qp) const
1082{
1083 Gradient du;
1084
1085 this->fixed_interior_gradient( var, qp, du );
1086
1087 return du;
1088}
1089
1090
1091template<typename OutputType>
1092void FEMContext::fixed_interior_gradient(unsigned int var, unsigned int qp,
1093 OutputType & du) const
1094{
1095 this->some_gradient
1096 <OutputType,
1097 &FEMContext::get_element_fe
1098 <typename TensorTools::MakeReal
1099 <typename TensorTools::DecrementRank
1100 <OutputType>::type>::type>,
1101 &DiffContext::get_elem_fixed_solution>
1102 (var, qp, du);
1103}
1104
1105
1106
1107#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1108Tensor FEMContext::fixed_interior_hessian(unsigned int var, unsigned int qp) const
1109{
1110 Tensor d2u;
1111
1112 this->fixed_interior_hessian( var, qp, d2u );
1113
1114 return d2u;
1115}
1116
1117
1118template<typename OutputType>
1119void FEMContext::fixed_interior_hessian(unsigned int var, unsigned int qp,
1120 OutputType & d2u) const
1121{
1122 this->some_hessian<OutputType,
1123 &FEMContext::get_element_fe
1124 <typename TensorTools::MakeReal
1125 <typename TensorTools::DecrementRank
1126 <typename TensorTools::DecrementRank
1127 <OutputType>::type>::type>::type>,
1128 &DiffContext::get_elem_fixed_solution>(var, qp, d2u);
1129}
1130#endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1131
1132
1133
1134Number FEMContext::fixed_side_value(unsigned int var, unsigned int qp) const
1135{
1136 Number u = 0.;
1137
1138 this->fixed_side_value( var, qp, u );
1139
1140 return u;
1141}
1142
1143
1144template<typename OutputType>
1145void FEMContext::fixed_side_value(unsigned int var, unsigned int qp,
1146 OutputType & u) const
1147{
1148 this->some_value
1149 <OutputType,
1150 &FEMContext::get_side_fe
1151 <typename TensorTools::MakeReal<OutputType>::type>,
1152 &DiffContext::get_elem_fixed_solution>
1153 (var, qp, u);
1154}
1155
1156
1157
1158Gradient FEMContext::fixed_side_gradient(unsigned int var, unsigned int qp) const
1159{
1160 Gradient du;
1161
1162 this->fixed_side_gradient( var, qp, du );
1163
1164 return du;
1165}
1166
1167
1168template<typename OutputType>
1169void FEMContext::fixed_side_gradient(unsigned int var, unsigned int qp,
1170 OutputType & du) const
1171{
1172 this->some_gradient<OutputType,
1173 &FEMContext::get_side_fe
1174 <typename TensorTools::MakeReal
1175 <typename TensorTools::DecrementRank
1176 <OutputType>::type>::type>,
1177 &DiffContext::get_elem_fixed_solution>(var, qp, du);
1178}
1179
1180
1181
1182#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1183Tensor FEMContext::fixed_side_hessian(unsigned int var, unsigned int qp) const
1184{
1185 Tensor d2u;
1186
1187 this->fixed_side_hessian( var, qp, d2u );
1188
1189 return d2u;
1190}
1191
1192template<typename OutputType>
1193void FEMContext::fixed_side_hessian(unsigned int var, unsigned int qp,
1194 OutputType & d2u) const
1195{
1196 this->some_hessian<OutputType,
1197 &FEMContext::get_side_fe
1198 <typename TensorTools::MakeReal
1199 <typename TensorTools::DecrementRank
1200 <typename TensorTools::DecrementRank
1201 <OutputType>::type>::type>::type>,
1202 &DiffContext::get_elem_fixed_solution>(var, qp, d2u);
1203}
1204#endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1205
1206
1207
1208Number FEMContext::fixed_point_value(unsigned int var, const Point & p) const
1209{
1210 Number u = 0.;
1211
1212 this->fixed_point_value( var, p, u );
1213
1214 return u;
1215}
1216
1217template<typename OutputType>
1218void FEMContext::fixed_point_value(unsigned int var,
1219 const Point & p,
1220 OutputType & u,
1221 const Real tolerance) const
1222{
1223 typedef typename TensorTools::MakeReal<OutputType>::type OutputShape;
1224
1225 // Get local-to-global dof index lookup
1226 const unsigned int n_dofs = cast_int<unsigned int>
1227 (this->get_dof_indices(var).size());
1228
1229 // Get current local coefficients
1230 libmesh_assert_greater (_elem_fixed_subsolutions.size(), var);
1231 const DenseSubVector<Number> & coef = this->get_elem_fixed_solution(var);
1232
1233 // Get finite element object
1234 FEGenericBase<OutputShape> * fe = nullptr;
1235 this->get_element_fe<OutputShape>( var, fe, this->get_elem_dim() );
1236
1237 // Build a FE for calculating u(p)
1238 FEGenericBase<OutputShape> * fe_new =
1239 this->build_new_fe( fe, p, tolerance, 0 );
1240
1241 // Get the values of the shape function derivatives
1242 const std::vector<std::vector<OutputShape>> & phi = fe_new->get_phi();
1243
1244 u = 0.;
1245
1246 for (unsigned int l=0; l != n_dofs; l++)
1247 u += phi[l][0] * coef(l);
1248
1249 return;
1250}
1251
1252
1253
1254Gradient FEMContext::fixed_point_gradient(unsigned int var, const Point & p) const
1255{
1256 Gradient grad_u;
1257
1258 this->fixed_point_gradient( var, p, grad_u );
1259
1260 return grad_u;
1261}
1262
1263
1264
1265template<typename OutputType>
1266void FEMContext::fixed_point_gradient(unsigned int var,
1267 const Point & p,
1268 OutputType & grad_u,
1269 const Real tolerance) const
1270{
1271 typedef typename TensorTools::MakeReal
1272 <typename TensorTools::DecrementRank<OutputType>::type>::type
1273 OutputShape;
1274
1275 // Get local-to-global dof index lookup
1276 const unsigned int n_dofs = cast_int<unsigned int>
1277 (this->get_dof_indices(var).size());
1278
1279 // Get current local coefficients
1280 libmesh_assert_greater (_elem_fixed_subsolutions.size(), var);
1281 const DenseSubVector<Number> & coef = this->get_elem_fixed_solution(var);
1282
1283 // Get finite element object
1284 FEGenericBase<OutputShape> * fe = nullptr;
1285 this->get_element_fe<OutputShape>( var, fe, this->get_elem_dim() );
1286
1287 // Build a FE for calculating u(p)
1288 FEGenericBase<OutputShape> * fe_new =
1289 this->build_new_fe( fe, p, tolerance, 1 );
1290
1291 // Get the values of the shape function derivatives
1292 const std::vector<std::vector<typename FEGenericBase<OutputShape>::OutputGradient>> & dphi = fe_new->get_dphi();
1293
1294 grad_u = 0.0;
1295
1296 for (unsigned int l=0; l != n_dofs; l++)
1297 grad_u.add_scaled(dphi[l][0], coef(l));
1298
1299 return;
1300}
1301
1302
1303#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1304
1305Tensor FEMContext::fixed_point_hessian(unsigned int var, const Point & p) const
1306{
1307 Tensor hess_u;
1308
1309 this->fixed_point_hessian( var, p, hess_u );
1310
1311 return hess_u;
1312}
1313
1314
1315
1316template<typename OutputType>
1317void FEMContext::fixed_point_hessian(unsigned int var,
1318 const Point & p,
1319 OutputType & hess_u,
1320 const Real tolerance) const
1321{
1322 typedef typename TensorTools::DecrementRank<OutputType>::type Rank1Decrement;
1323 typedef typename TensorTools::DecrementRank<Rank1Decrement>::type Rank2Decrement;
1324 typedef typename TensorTools::MakeReal<Rank2Decrement>::type OutputShape;
1325
1326 // Get local-to-global dof index lookup
1327 const unsigned int n_dofs = cast_int<unsigned int>
1328 (this->get_dof_indices(var).size());
1329
1330 // Get current local coefficients
1331 libmesh_assert_greater (_elem_fixed_subsolutions.size(), var);
1332 const DenseSubVector<Number> & coef = this->get_elem_fixed_solution(var);
1333
1334 // Get finite element object
1335 FEGenericBase<OutputShape> * fe = nullptr;
1336 this->get_element_fe<OutputShape>( var, fe, this->get_elem_dim() );
1337
1338 // Build a FE for calculating u(p)
1339 FEGenericBase<OutputShape> * fe_new =
1340 this->build_new_fe( fe, p, tolerance, 2 );
1341
1342 // Get the values of the shape function derivatives
1343 const std::vector<std::vector<typename FEGenericBase<OutputShape>::OutputTensor>> & d2phi = fe_new->get_d2phi();
1344
1345 hess_u = 0.0;
1346
1347 for (unsigned int l=0; l != n_dofs; l++)
1348 hess_u.add_scaled(d2phi[l][0], coef(l));
1349
1350 return;
1351}
1352
1353#endif // LIBMESH_ENABLE_SECOND_DERIVATIVES
1354
1355
1356
1357template<typename OutputType>
1358void FEMContext::interior_rate(unsigned int var, unsigned int qp,
1359 OutputType & u) const
1360{
1361 this->some_value<OutputType,
1362 &FEMContext::get_element_fe
1363 <typename TensorTools::MakeReal<OutputType>::type>,
1364 &DiffContext::get_elem_solution_rate>(var, qp, u);
1365}
1366
1367template<typename OutputType>
1368void FEMContext::interior_rate_gradient(unsigned int var, unsigned int qp,
1369 OutputType & dudot) const
1370{
1371 this->some_gradient<OutputType,
1372 &FEMContext::get_element_fe<typename TensorTools::MakeReal
1373 <typename TensorTools::DecrementRank
1374 <OutputType>::type>::type>,
1375 &DiffContext::get_elem_solution_rate>(var, qp, dudot);
1376}
1377
1378template<typename OutputType>
1379void FEMContext::side_rate(unsigned int var, unsigned int qp,
1380 OutputType & u) const
1381{
1382 this->some_value<OutputType,
1383 &FEMContext::get_side_fe
1384 <typename TensorTools::MakeReal<OutputType>::type>,
1385 &DiffContext::get_elem_solution_rate>(var, qp, u);
1386}
1387
1388template<typename OutputType>
1389void FEMContext::interior_accel(unsigned int var, unsigned int qp,
1390 OutputType & u) const
1391{
1392 this->some_value<OutputType,
1393 &FEMContext::get_element_fe
1394 <typename TensorTools::MakeReal<OutputType>::type>,
1395 &DiffContext::get_elem_solution_accel>(var, qp, u);
1396}
1397
1398
1399
1400template<typename OutputType>
1401void FEMContext::side_accel(unsigned int var, unsigned int qp,
1402 OutputType & u) const
1403{
1404 this->some_value<OutputType,
1405 &FEMContext::get_side_fe
1406 <typename TensorTools::MakeReal<OutputType>::type>,
1407 &DiffContext::get_elem_solution_accel>(var, qp, u);
1408}
1409
1410
1411
1412void FEMContext::elem_reinit(Real theta)
1413{
1414 // Update the "time" variable of this context object
1415 this->_update_time_from_system(theta);
1416
1417 // Handle a moving element if necessary.
1418 if (_mesh_sys)
1419 {
1420 // We assume that the ``default'' state
1421 // of the mesh is its final, theta=1.0
1422 // position, so we don't bother with
1423 // mesh motion in that case.
1424 if (theta != 1.0)
1425 {
1426 // FIXME - ALE is not threadsafe yet!
1427 libmesh_assert_equal_to (libMesh::n_threads(), 1);
1428
1429 elem_position_set(theta);
1430 }
1431 elem_fe_reinit();
1432 }
1433}
1434
1435
1436void FEMContext::elem_side_reinit(Real theta)
1437{
1438 // Update the "time" variable of this context object
1439 this->_update_time_from_system(theta);
1440
1441 // Handle a moving element if necessary
1442 if (_mesh_sys)
1443 {
1444 // FIXME - not threadsafe yet!
1445 elem_position_set(theta);
1446 side_fe_reinit();
1447 }
1448}
1449
1450
1451void FEMContext::elem_edge_reinit(Real theta)
1452{
1453 // Update the "time" variable of this context object
1454 this->_update_time_from_system(theta);
1455
1456 // Handle a moving element if necessary
1457 if (_mesh_sys)
1458 {
1459 // FIXME - not threadsafe yet!
1460 elem_position_set(theta);
1461 edge_fe_reinit();
1462 }
1463}
1464
1465
1466void FEMContext::nonlocal_reinit(Real theta)
1467{
1468 // Update the "time" variable of this context object
1469 this->_update_time_from_system(theta);
1470
1471 // We can reuse the Elem FE safely here.
1472 elem_fe_reinit();
1473}
1474
1475
1476void FEMContext::elem_fe_reinit(const std::vector<Point> * const pts)
1477{
1478 // Initialize all the interior FE objects on elem.
1479 // Logging of FE::reinit is done in the FE functions
1480 // We only reinit the FE objects for the current element
1481 // dimension
1482 const unsigned char dim = this->get_elem_dim();
1483
1484 libmesh_assert( !_element_fe[dim].empty() );
1485
1486 for (const auto & pr : _element_fe[dim])
1487 {
1488 if (this->has_elem())
1489 pr.second->reinit(&(this->get_elem()), pts);
1490 // If !this->has_elem(), then still might need to reinit for a
1491 // SCALAR variable; everything else will depend on an elem
1492 else if (pr.first.family == SCALAR)
1493 pr.second->reinit(nullptr);
1494 }
1495}
1496
1497
1498void FEMContext::side_fe_reinit ()
1499{
1500 // Initialize all the side FE objects on elem/side.
1501 // Logging of FE::reinit is done in the FE functions
1502 // We only reinit the FE objects for the current element
1503 // dimension
1504 const unsigned char dim = this->get_elem_dim();
1505
1506 libmesh_assert( !_side_fe[dim].empty() );
1507
1508 for (auto & pr : _side_fe[dim])
1509 pr.second->reinit(&(this->get_elem()), this->get_side());
1510}
1511
1512
1513
1514void FEMContext::edge_fe_reinit ()
1515{
1516 libmesh_assert_equal_to (this->get_elem_dim(), 3);
1517
1518 // Initialize all the interior FE objects on elem/edge.
1519 // Logging of FE::reinit is done in the FE functions
1520 for (auto & pr : _edge_fe)
1521 pr.second->edge_reinit(&(this->get_elem()), this->get_edge());
1522}
1523
1524
1525
1526void FEMContext::elem_position_get()
1527{
1528 // This is too expensive to call unless we've been asked to move the mesh
1529 libmesh_assert (_mesh_sys);
1530
1531 // This will probably break with threading when two contexts are
1532 // operating on elements which share a node
1533 libmesh_assert_equal_to (libMesh::n_threads(), 1);
1534
1535 // If the coordinate data is in our own system, it's already
1536 // been set up for us
1537 // if (_mesh_sys == this->number())
1538 // {
1539 unsigned int n_nodes = this->get_elem().n_nodes();
1540
1541#ifndef NDEBUG
1542 const unsigned char dim = this->get_elem_dim();
1543
1544 // For simplicity we demand that mesh coordinates be stored
1545 // in a format that allows a direct copy
1546 libmesh_assert(this->get_mesh_x_var() == libMesh::invalid_uint ||
1547 (this->get_element_fe(this->get_mesh_x_var(), dim)->get_fe_type().family
1548 == FEMap::map_fe_type(this->get_elem()) &&
1549 this->get_element_fe(this->get_mesh_x_var(), dim)->get_fe_type().order.get_order()
1550 == this->get_elem().default_order()));
1551 libmesh_assert(this->get_mesh_y_var() == libMesh::invalid_uint ||
1552 (this->get_element_fe(this->get_mesh_y_var(), dim)->get_fe_type().family
1553 == FEMap::map_fe_type(this->get_elem()) &&
1554 this->get_element_fe(this->get_mesh_y_var(), dim)->get_fe_type().order.get_order()
1555 == this->get_elem().default_order()));
1556 libmesh_assert(this->get_mesh_z_var() == libMesh::invalid_uint ||
1557 (this->get_element_fe(this->get_mesh_z_var(), dim)->get_fe_type().family
1558 == FEMap::map_fe_type(this->get_elem()) &&
1559 this->get_element_fe(this->get_mesh_z_var(), dim)->get_fe_type().order.get_order()
1560 == this->get_elem().default_order()));
1561#endif
1562
1563 // Get degree of freedom coefficients from point coordinates
1564 if (this->get_mesh_x_var() != libMesh::invalid_uint)
1565 for (unsigned int i=0; i != n_nodes; ++i)
1566 (this->get_elem_solution(this->get_mesh_x_var()))(i) = this->get_elem().point(i)(0);
1567
1568 if (this->get_mesh_y_var() != libMesh::invalid_uint)
1569 for (unsigned int i=0; i != n_nodes; ++i)
1570 (this->get_elem_solution(this->get_mesh_y_var()))(i) = this->get_elem().point(i)(1);
1571
1572 if (this->get_mesh_z_var() != libMesh::invalid_uint)
1573 for (unsigned int i=0; i != n_nodes; ++i)
1574 (this->get_elem_solution(this->get_mesh_z_var()))(i) = this->get_elem().point(i)(2);
1575 // }
1576 // FIXME - If the coordinate data is not in our own system, someone
1577 // had better get around to implementing that... - RHS
1578 // else
1579 // {
1580 // libmesh_not_implemented();
1581 // }
1582}
1583
1584
1585
1586void FEMContext::set_jacobian_tolerance(Real tol)
1587{
1588 for (auto & m : _element_fe)
1589 for (auto & pr : m)
1590 pr.second->get_fe_map().set_jacobian_tolerance(tol);
1591
1592 for (auto & m : _side_fe)
1593 for (auto & pr : m)
1594 pr.second->get_fe_map().set_jacobian_tolerance(tol);
1595
1596 for (auto & pr : _edge_fe)
1597 pr.second->get_fe_map().set_jacobian_tolerance(tol);
1598}
1599
1600
1601
1602// We can ignore the theta argument in the current use of this
1603// function, because elem_subsolutions will already have been set to
1604// the theta value.
1605//
1606// To enable loose mesh movement coupling things will need to change.
1607void FEMContext::_do_elem_position_set(Real)
1608{
1609 // This is too expensive to call unless we've been asked to move the mesh
1610 libmesh_assert (_mesh_sys);
1611
1612 // This will probably break with threading when two contexts are
1613 // operating on elements which share a node
1614 libmesh_assert_equal_to (libMesh::n_threads(), 1);
1615
1616 // If the coordinate data is in our own system, it's already
1617 // been set up for us, and we can ignore our input parameter theta
1618 // if (_mesh_sys == this->number())
1619 // {
1620 unsigned int n_nodes = this->get_elem().n_nodes();
1621
1622#ifndef NDEBUG
1623 const unsigned char dim = this->get_elem_dim();
1624
1625 // For simplicity we demand that mesh coordinates be stored
1626 // in a format that allows a direct copy
1627 libmesh_assert(this->get_mesh_x_var() == libMesh::invalid_uint ||
1628 (this->get_element_fe(this->get_mesh_x_var(), dim)->get_fe_type().family
1629 == FEMap::map_fe_type(this->get_elem()) &&
1630 this->get_elem_solution(this->get_mesh_x_var()).size() == n_nodes));
1631 libmesh_assert(this->get_mesh_y_var() == libMesh::invalid_uint ||
1632 (this->get_element_fe(this->get_mesh_y_var(), dim)->get_fe_type().family
1633 == FEMap::map_fe_type(this->get_elem()) &&
1634 this->get_elem_solution(this->get_mesh_y_var()).size() == n_nodes));
1635 libmesh_assert(this->get_mesh_z_var() == libMesh::invalid_uint ||
1636 (this->get_element_fe(this->get_mesh_z_var(), dim)->get_fe_type().family
1637 == FEMap::map_fe_type(this->get_elem()) &&
1638 this->get_elem_solution(this->get_mesh_z_var()).size() == n_nodes));
1639#endif
1640
1641 // Set the new point coordinates
1642 if (this->get_mesh_x_var() != libMesh::invalid_uint)
1643 for (unsigned int i=0; i != n_nodes; ++i)
1644 const_cast<Elem &>(this->get_elem()).point(i)(0) =
1645 libmesh_real(this->get_elem_solution(this->get_mesh_x_var())(i));
1646
1647 if (this->get_mesh_y_var() != libMesh::invalid_uint)
1648 for (unsigned int i=0; i != n_nodes; ++i)
1649 const_cast<Elem &>(this->get_elem()).point(i)(1) =
1650 libmesh_real(this->get_elem_solution(this->get_mesh_y_var())(i));
1651
1652 if (this->get_mesh_z_var() != libMesh::invalid_uint)
1653 for (unsigned int i=0; i != n_nodes; ++i)
1654 const_cast<Elem &>(this->get_elem()).point(i)(2) =
1655 libmesh_real(this->get_elem_solution(this->get_mesh_z_var())(i));
1656 // }
1657 // FIXME - If the coordinate data is not in our own system, someone
1658 // had better get around to implementing that... - RHS
1659 // else
1660 // {
1661 // libmesh_not_implemented();
1662 // }
1663}
1664
1665
1666
1667
1668
1669/*
1670 void FEMContext::reinit(const FEMSystem & sys, Elem * e)
1671 {
1672 // Initialize our elem pointer, algebraic objects
1673 this->pre_fe_reinit(e);
1674
1675 // Moving the mesh may be necessary
1676 // Reinitializing the FE objects is definitely necessary
1677 this->elem_reinit(1.);
1678 }
1679*/
1680
1681
1682
1683void FEMContext::pre_fe_reinit(const System & sys, const Elem * e)
1684{
1685 this->set_elem(e);
1686
1687 if (algebraic_type() == CURRENT ||
1688 algebraic_type() == DOFS_ONLY)
1689 {
1690 // Initialize the per-element data for elem.
1691 if (this->has_elem())
1692 sys.get_dof_map().dof_indices (&(this->get_elem()), this->get_dof_indices());
1693 else
1694 // If !this->has_elem(), then we assume we are dealing with a SCALAR variable
1695 sys.get_dof_map().dof_indices
1696 (static_cast<Elem*>(nullptr), this->get_dof_indices());
1697 }
1698#ifdef LIBMESH_ENABLE_AMR
1699 else if (algebraic_type() == OLD ||
1700 algebraic_type() == OLD_DOFS_ONLY)
1701 {
1702 // Initialize the per-element data for elem.
1703 if (this->has_elem())
1704 sys.get_dof_map().old_dof_indices (&(this->get_elem()), this->get_dof_indices());
1705 else
1706 // If !this->has_elem(), then we assume we are dealing with a SCALAR variable
1707 sys.get_dof_map().old_dof_indices
1708 (static_cast<Elem*>(nullptr), this->get_dof_indices());
1709 }
1710#endif // LIBMESH_ENABLE_AMR
1711
1712 const unsigned int n_dofs = cast_int<unsigned int>
1713 (this->get_dof_indices().size());
1714 const unsigned int n_qoi = sys.n_qois();
1715
1716 if (this->algebraic_type() != NONE &&
1717 this->algebraic_type() != DOFS_ONLY &&
1718 this->algebraic_type() != OLD_DOFS_ONLY)
1719 {
1720 // This also resizes elem_solution
1721 if (_custom_solution == nullptr)
1722 sys.current_local_solution->get(this->get_dof_indices(), this->get_elem_solution().get_values());
1723 else
1724 _custom_solution->get(this->get_dof_indices(), this->get_elem_solution().get_values());
1725
1726 if (sys.use_fixed_solution)
1727 this->get_elem_fixed_solution().resize(n_dofs);
1728
1729 // Only make space for these if we're using DiffSystem
1730 // This is assuming *only* DiffSystem is using elem_solution_rate/accel
1731 const DifferentiableSystem * diff_system = dynamic_cast<const DifferentiableSystem *>(&sys);
1732 if (diff_system)
1733 {
1734 // Now, we only need these if the solver is unsteady
1735 if (!diff_system->get_time_solver().is_steady())
1736 {
1737 this->get_elem_solution_rate().resize(n_dofs);
1738
1739 // We only need accel space if the TimeSolver is second order
1740 const UnsteadySolver & time_solver = cast_ref<const UnsteadySolver &>(diff_system->get_time_solver());
1741
1742 if (time_solver.time_order() >= 2 || !diff_system->get_second_order_vars().empty())
1743 this->get_elem_solution_accel().resize(n_dofs);
1744 }
1745 }
1746
1747 if (algebraic_type() != OLD)
1748 {
1749 // These resize calls also zero out the residual and jacobian
1750 this->get_elem_residual().resize(n_dofs);
1751 if (this->_have_local_matrices)
1752 this->get_elem_jacobian().resize(n_dofs, n_dofs);
1753
1754 this->get_qoi_derivatives().resize(n_qoi);
1755 this->_elem_qoi_subderivatives.resize(n_qoi);
1756 for (std::size_t q=0; q != n_qoi; ++q)
1757 (this->get_qoi_derivatives())[q].resize(n_dofs);
1758 }
1759 }
1760
1761 // Initialize the per-variable data for elem.
1762 {
1763 unsigned int sub_dofs = 0;
1764 for (auto i : make_range(sys.n_vars()))
1765 {
1766 if (algebraic_type() == CURRENT ||
1767 algebraic_type() == DOFS_ONLY)
1768 {
1769 if (this->has_elem())
1770 sys.get_dof_map().dof_indices (&(this->get_elem()), this->get_dof_indices(i), i);
1771 else
1772 // If !this->has_elem(), then we assume we are dealing with a SCALAR variable
1773 sys.get_dof_map().dof_indices
1774 (static_cast<Elem*>(nullptr), this->get_dof_indices(i), i);
1775 }
1776#ifdef LIBMESH_ENABLE_AMR
1777 else if (algebraic_type() == OLD ||
1778 algebraic_type() == OLD_DOFS_ONLY)
1779 {
1780 if (this->has_elem())
1781 sys.get_dof_map().old_dof_indices (&(this->get_elem()), this->get_dof_indices(i), i);
1782 else
1783 // If !this->has_elem(), then we assume we are dealing with a SCALAR variable
1784 sys.get_dof_map().old_dof_indices
1785 (static_cast<Elem*>(nullptr), this->get_dof_indices(i), i);
1786 }
1787#endif // LIBMESH_ENABLE_AMR
1788
1789 if (this->algebraic_type() != NONE &&
1790 this->algebraic_type() != DOFS_ONLY &&
1791 this->algebraic_type() != OLD_DOFS_ONLY)
1792 {
1793 const unsigned int n_dofs_var = cast_int<unsigned int>
1794 (this->get_dof_indices(i).size());
1795
1796
1797 if (!_active_vars ||
1798 std::binary_search(_active_vars->begin(),
1799 _active_vars->end(), i))
1800 {
1801 this->get_elem_solution(i).reposition
1802 (sub_dofs, n_dofs_var);
1803
1804 // Only make space for these if we're using DiffSystem
1805 // This is assuming *only* DiffSystem is using elem_solution_rate/accel
1806 const DifferentiableSystem * diff_system = dynamic_cast<const DifferentiableSystem *>(&sys);
1807 if (diff_system)
1808 {
1809 // Now, we only need these if the solver is unsteady
1810 if (!diff_system->get_time_solver().is_steady())
1811 {
1812 this->get_elem_solution_rate(i).reposition
1813 (sub_dofs, n_dofs_var);
1814
1815 // We only need accel space if the TimeSolver is second order
1816 const UnsteadySolver & time_solver = cast_ref<const UnsteadySolver &>(diff_system->get_time_solver());
1817
1818 if (time_solver.time_order() >= 2 || !diff_system->get_second_order_vars().empty())
1819 this->get_elem_solution_accel(i).reposition
1820 (sub_dofs, n_dofs_var);
1821 }
1822 }
1823
1824 if (sys.use_fixed_solution)
1825 this->get_elem_fixed_solution(i).reposition
1826 (sub_dofs, n_dofs_var);
1827
1828 if (algebraic_type() != OLD)
1829 {
1830 this->get_elem_residual(i).reposition
1831 (sub_dofs, n_dofs_var);
1832
1833 for (std::size_t q=0; q != n_qoi; ++q)
1834 this->get_qoi_derivatives(q,i).reposition
1835 (sub_dofs, n_dofs_var);
1836
1837 if (this->_have_local_matrices)
1838 {
1839 for (unsigned int j=0; j != i; ++j)
1840 {
1841 const unsigned int n_dofs_var_j =
1842 cast_int<unsigned int>
1843 (this->get_dof_indices(j).size());
1844
1845 this->get_elem_jacobian(i,j).reposition
1846 (sub_dofs, this->get_elem_residual(j).i_off(),
1847 n_dofs_var, n_dofs_var_j);
1848 this->get_elem_jacobian(j,i).reposition
1849 (this->get_elem_residual(j).i_off(), sub_dofs,
1850 n_dofs_var_j, n_dofs_var);
1851 }
1852 this->get_elem_jacobian(i,i).reposition
1853 (sub_dofs, sub_dofs,
1854 n_dofs_var,
1855 n_dofs_var);
1856 }
1857 }
1858 }
1859
1860 sub_dofs += n_dofs_var;
1861 }
1862 }
1863
1864 if (this->algebraic_type() != NONE &&
1865 this->algebraic_type() != DOFS_ONLY &&
1866 this->algebraic_type() != OLD &&
1867 this->algebraic_type() != OLD_DOFS_ONLY)
1868 libmesh_assert_equal_to (sub_dofs, n_dofs);
1869 }
1870
1871 // Now do the localization for the user requested vectors
1872 if (this->algebraic_type() != NONE &&
1873 this->algebraic_type() != DOFS_ONLY &&
1874 this->algebraic_type() != OLD_DOFS_ONLY)
1875 {
1876 DiffContext::localized_vectors_iterator localized_vec_it = this->_localized_vectors.begin();
1877 const DiffContext::localized_vectors_iterator localized_vec_end = this->_localized_vectors.end();
1878
1879 for (; localized_vec_it != localized_vec_end; ++localized_vec_it)
1880 {
1881 const NumericVector<Number> & current_localized_vector = *localized_vec_it->first;
1882 DenseVector<Number> & target_vector = localized_vec_it->second.first;
1883
1884 current_localized_vector.get(this->get_dof_indices(), target_vector.get_values());
1885
1886 // Initialize the per-variable data for elem.
1887 unsigned int sub_dofs = 0;
1888 auto init_localized_var_data = [this, localized_vec_it, &sub_dofs](unsigned int i)
1889 {
1890 const unsigned int n_dofs_var = cast_int<unsigned int>
1891 (this->get_dof_indices(i).size());
1892
1893 // This is redundant with earlier initialization, isn't it? - RHS
1894 // sys.get_dof_map().dof_indices (&(this->get_elem()), this->get_dof_indices(i), i);
1895
1896 localized_vec_it->second.second[i].reposition
1897 (sub_dofs, n_dofs_var);
1898
1899 sub_dofs += n_dofs_var;
1900 };
1901
1902 if (_active_vars)
1903 for (auto v : *_active_vars)
1904 init_localized_var_data(v);
1905 else
1906 for (auto v : make_range(sys.n_vars()))
1907 init_localized_var_data(v);
1908
1909 libmesh_assert_equal_to (sub_dofs, n_dofs);
1910 }
1911 }
1912}
1913
1914void FEMContext::set_elem( const Elem * e )
1915{
1916 this->_elem = e;
1917
1918 // If e is nullptr, we assume it's SCALAR and set _elem_dim to 0.
1919 this->_elem_dim =
1920 cast_int<unsigned char>(this->_elem ? this->_elem->dim() : 0);
1921}
1922
1923void FEMContext::_update_time_from_system(Real theta)
1924{
1925 // Update the "time" variable based on the value of theta. For this
1926 // to work, we need to know the value of deltat, a pointer to which is now
1927 // stored by our parent DiffContext class. Note: get_deltat_value() will
1928 // assert in debug mode if the requested pointer is nullptr.
1929 const Real deltat = this->get_deltat_value();
1930
1931 this->set_time(theta*(this->get_system_time() + deltat) + (1.-theta)*this->get_system_time());
1932}
1933
1934
1935
1936template<>
1937FEGenericBase<Real> *
1938FEMContext::cached_fe( const unsigned int elem_dim,
1939 const FEType fe_type,
1940 const int get_derivative_level ) const
1941{
1942#ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
1943 const bool fe_needs_inf =
1944 this->has_elem() && this->get_elem().infinite();
1945#endif
1946
1947 if (!_real_fe ||
1948 elem_dim != _real_fe->get_dim() ||
1949 fe_type != _real_fe->get_fe_type() ||
1950 get_derivative_level != _real_fe_derivative_level)
1951 {
1952 _real_fe_derivative_level = get_derivative_level;
1953
1954 _real_fe =
1955#ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
1956 fe_needs_inf ?
1957 FEGenericBase<Real>::build_InfFE(elem_dim, fe_type) :
1958#endif
1959 FEGenericBase<Real>::build(elem_dim, fe_type);
1960 }
1961
1962#ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
1963 else if (fe_needs_inf && !_real_fe_is_inf)
1964 _real_fe =
1965 FEGenericBase<Real>::build_InfFE(elem_dim, fe_type);
1966 else if (!fe_needs_inf && _real_fe_is_inf)
1967 _real_fe =
1968 FEGenericBase<Real>::build(elem_dim, fe_type);
1969
1970 _real_fe_is_inf =
1971 (this->has_elem() && this->get_elem().infinite());
1972#endif
1973
1974 return _real_fe.get();
1975}
1976
1977
1978template<>
1979FEGenericBase<RealGradient> *
1980FEMContext::cached_fe( const unsigned int elem_dim,
1981 const FEType fe_type,
1982 const int get_derivative_level ) const
1983{
1984#ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
1985 const bool fe_needs_inf =
1986 this->has_elem() && this->get_elem().infinite();
1987#endif
1988
1989 if (!_real_grad_fe ||
1990 elem_dim != _real_grad_fe->get_dim() ||
1991 fe_type != _real_grad_fe->get_fe_type() ||
1992 get_derivative_level != _real_grad_fe_derivative_level)
1993 {
1994 _real_grad_fe_derivative_level = get_derivative_level;
1995
1996 _real_grad_fe =
1997#ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
1998 fe_needs_inf ?
1999 FEGenericBase<RealGradient>::build_InfFE(elem_dim, fe_type) :
2000#endif
2001 FEGenericBase<RealGradient>::build(elem_dim, fe_type);
2002 }
2003
2004#ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
2005 else if (fe_needs_inf && !_real_grad_fe_is_inf)
2006 _real_grad_fe =
2007 FEGenericBase<RealGradient>::build_InfFE(elem_dim, fe_type);
2008 else if (!fe_needs_inf && _real_grad_fe_is_inf)
2009 _real_grad_fe =
2010 FEGenericBase<RealGradient>::build(elem_dim, fe_type);
2011
2012 _real_grad_fe_is_inf =
2013 (this->has_elem() && this->get_elem().infinite());
2014#endif
2015
2016 return _real_grad_fe.get();
2017}
2018
2019
2020
2021template<typename OutputShape>
2022FEGenericBase<OutputShape> *
2023FEMContext::build_new_fe( const FEGenericBase<OutputShape>* fe,
2024 const Point & p,
2025 const Real tolerance,
2026 const int get_derivative_level) const
2027{
2028 FEType fe_type = fe->get_fe_type();
2029
2030 // If we don't have an Elem to evaluate on, then the only functions
2031 // we can sensibly evaluate are the scalar dofs which are the same
2032 // everywhere.
2033 libmesh_assert(this->has_elem() || fe_type.family == SCALAR);
2034
2035#ifdef LIBMESH_ENABLE_AMR
2036 const int add_p_level = fe->add_p_level_in_reinit();
2037 if ((algebraic_type() == OLD) &&
2038 this->has_elem())
2039 {
2040 if (this->get_elem().p_refinement_flag() == Elem::JUST_REFINED)
2041 fe_type.order -= add_p_level;
2042 else if (this->get_elem().p_refinement_flag() == Elem::JUST_COARSENED)
2043 fe_type.order += add_p_level;
2044 }
2045#endif // LIBMESH_ENABLE_AMR
2046
2047 const unsigned int elem_dim = this->has_elem() ? this->get_elem().dim() : 0;
2048
2049 FEGenericBase<OutputShape>* fe_new =
2050 cached_fe<OutputShape>(elem_dim, fe_type, get_derivative_level);
2051#ifdef LIBMESH_ENABLE_AMR
2052 fe_new->add_p_level_in_reinit(add_p_level);
2053#endif // LIBMESH_ENABLE_AMR
2054
2055 // Map the physical co-ordinates to the master co-ordinates using the inverse_map from fe_interface.h
2056 // Build a vector of point co-ordinates to send to reinit
2057 Point master_point = this->has_elem() ?
2058 FEMap::inverse_map (elem_dim, &this->get_elem(), p, tolerance) :
2059 Point(0);
2060
2061 std::vector<Point> coor(1, master_point);
2062
2063 switch (get_derivative_level)
2064 {
2065 case -1:
2066 fe_new->get_phi();
2067 fe_new->get_dphi();
2068#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
2069 fe_new->get_d2phi();
2070#endif
2071 fe_new->get_curl_phi();
2072 break;
2073 case 0:
2074 fe_new->get_phi();
2075 break;
2076 case 1:
2077 fe_new->get_dphi();
2078 break;
2079 case 2:
2080#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
2081 fe_new->get_d2phi();
2082#else
2083 // here a different configuration is required.
2084 libmesh_not_implemented();
2085#endif
2086 break;
2087 case 3:
2088 fe_new->get_curl_phi();
2089 break;
2090 default:
2091 libmesh_error();
2092 }
2093
2094 // Reinitialize the element and compute the shape function values at coor
2095 if (this->has_elem())
2096 fe_new->reinit (&this->get_elem(), &coor);
2097 else
2098 // If !this->has_elem(), then we assume we are dealing with a SCALAR variable
2099 fe_new->reinit (nullptr, &coor);
2100
2101 return fe_new;
2102}
2103
2104
2105
2106
2107
2108// Instantiate member function templates
2109template LIBMESH_EXPORT void FEMContext::interior_value<Number>(unsigned int, unsigned int, Number &) const;
2110template LIBMESH_EXPORT void FEMContext::interior_values<Number>(unsigned int, const NumericVector<Number> &,
2111 std::vector<Number> &) const;
2112template LIBMESH_EXPORT void FEMContext::interior_value<Gradient>(unsigned int, unsigned int, Gradient &) const;
2113template LIBMESH_EXPORT void FEMContext::interior_values<Gradient>(unsigned int, const NumericVector<Number> &,
2114 std::vector<Gradient> &) const;
2115
2116template LIBMESH_EXPORT void FEMContext::interior_gradient<Gradient>(unsigned int, unsigned int, Gradient &) const;
2117template LIBMESH_EXPORT void FEMContext::interior_gradients<Gradient>(unsigned int, const NumericVector<Number> &,
2118 std::vector<Gradient> &) const;
2119template LIBMESH_EXPORT void FEMContext::interior_gradient<Tensor>(unsigned int, unsigned int, Tensor &) const;
2120template LIBMESH_EXPORT void FEMContext::interior_gradients<Tensor>(unsigned int, const NumericVector<Number> &,
2121 std::vector<Tensor> &) const;
2122
2123#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
2124template LIBMESH_EXPORT void FEMContext::interior_hessian<Tensor>(unsigned int, unsigned int, Tensor &) const;
2125template LIBMESH_EXPORT void FEMContext::interior_hessians<Tensor>(unsigned int, const NumericVector<Number> &,
2126 std::vector<Tensor> &) const;
2127//FIXME: Not everything is implemented yet for second derivatives of RealGradients
2128//template LIBMESH_EXPORT void FEMContext::interior_hessian<??>(unsigned int, unsigned int, ??&) const;
2129//template LIBMESH_EXPORT void FEMContext::interior_hessians<??>(unsigned int, const NumericVector<Number> &,
2130// std::vector<??> &) const;
2131#endif
2132
2133template LIBMESH_EXPORT void FEMContext::interior_curl<Gradient>(unsigned int, unsigned int, Gradient &) const;
2134
2135template LIBMESH_EXPORT void FEMContext::interior_div<Number>(unsigned int, unsigned int, Number &) const;
2136
2137template LIBMESH_EXPORT void FEMContext::side_value<Number>(unsigned int, unsigned int, Number &) const;
2138template LIBMESH_EXPORT void FEMContext::side_value<Gradient>(unsigned int, unsigned int, Gradient &) const;
2139template LIBMESH_EXPORT void FEMContext::side_values<Number>(unsigned int, const NumericVector<Number> &,
2140 std::vector<Number> &) const;
2141template LIBMESH_EXPORT void FEMContext::side_values<Gradient>(unsigned int, const NumericVector<Number> &,
2142 std::vector<Gradient> &) const;
2143
2144template LIBMESH_EXPORT void FEMContext::side_gradient<Gradient>(unsigned int, unsigned int, Gradient &) const;
2145template LIBMESH_EXPORT void FEMContext::side_gradients<Gradient>(unsigned int, const NumericVector<Number> &,
2146 std::vector<Gradient> &) const;
2147template LIBMESH_EXPORT void FEMContext::side_gradient<Tensor>(unsigned int, unsigned int, Tensor &) const;
2148template LIBMESH_EXPORT void FEMContext::side_gradients<Tensor>(unsigned int, const NumericVector<Number> &,
2149 std::vector<Tensor> &) const;
2150
2151
2152#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
2153template LIBMESH_EXPORT void FEMContext::side_hessian<Tensor>(unsigned int, unsigned int, Tensor &) const;
2154template LIBMESH_EXPORT void FEMContext::side_hessians<Tensor>(unsigned int, const NumericVector<Number> &,
2155 std::vector<Tensor> &) const;
2156//FIXME: Not everything is implemented yet for second derivatives of RealGradients
2157//template LIBMESH_EXPORT void FEMContext::side_hessian<??>(unsigned int, unsigned int,
2158// ??&) const;
2159//template LIBMESH_EXPORT void FEMContext::side_hessians<??>(unsigned int, const NumericVector<Number> &,
2160// std::vector<??> &) const;
2161#endif
2162
2163template LIBMESH_EXPORT void FEMContext::point_value<Number>(unsigned int, const Point &, Number &, const Real) const;
2164template LIBMESH_EXPORT void FEMContext::point_value<Gradient>(unsigned int, const Point &, Gradient &, const Real) const;
2165
2166template LIBMESH_EXPORT void FEMContext::point_gradient<Gradient>(unsigned int, const Point &, Gradient &, const Real) const;
2167template LIBMESH_EXPORT void FEMContext::point_gradient<Tensor>(unsigned int, const Point &, Tensor &, const Real) const;
2168
2169#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
2170template LIBMESH_EXPORT void FEMContext::point_hessian<Tensor>(unsigned int, const Point &, Tensor &, const Real) const;
2171//FIXME: Not everything is implemented yet for second derivatives of RealGradients
2172//template LIBMESH_EXPORT void FEMContext::point_hessian<??>(unsigned int, const Point &, ??&) const;
2173#endif
2174
2175template LIBMESH_EXPORT void FEMContext::point_curl<Gradient>(unsigned int, const Point &, Gradient &, const Real) const;
2176
2177template LIBMESH_EXPORT void FEMContext::fixed_interior_value<Number>(unsigned int, unsigned int, Number &) const;
2178template LIBMESH_EXPORT void FEMContext::fixed_interior_value<Gradient>(unsigned int, unsigned int, Gradient &) const;
2179
2180template LIBMESH_EXPORT void FEMContext::fixed_interior_gradient<Gradient>(unsigned int, unsigned int, Gradient &) const;
2181template LIBMESH_EXPORT void FEMContext::fixed_interior_gradient<Tensor>(unsigned int, unsigned int, Tensor &) const;
2182
2183#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
2184template LIBMESH_EXPORT void FEMContext::fixed_interior_hessian<Tensor>(unsigned int, unsigned int, Tensor &) const;
2185//FIXME: Not everything is implemented yet for second derivatives of RealGradients
2186//template LIBMESH_EXPORT void FEMContext::fixed_interior_hessian<??>(unsigned int, unsigned int, ??&) const;
2187#endif
2188
2189template LIBMESH_EXPORT void FEMContext::fixed_side_value<Number>(unsigned int, unsigned int, Number &) const;
2190template LIBMESH_EXPORT void FEMContext::fixed_side_value<Gradient>(unsigned int, unsigned int, Gradient &) const;
2191
2192template LIBMESH_EXPORT void FEMContext::fixed_side_gradient<Gradient>(unsigned int, unsigned int, Gradient &) const;
2193template LIBMESH_EXPORT void FEMContext::fixed_side_gradient<Tensor>(unsigned int, unsigned int, Tensor &) const;
2194
2195#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
2196template LIBMESH_EXPORT void FEMContext::fixed_side_hessian<Tensor>(unsigned int, unsigned int, Tensor &) const;
2197//FIXME: Not everything is implemented yet for second derivatives of RealGradients
2198//template LIBMESH_EXPORT void FEMContext::fixed_side_hessian<??>(unsigned int, unsigned int, ??&) const;
2199#endif
2200
2201template LIBMESH_EXPORT void FEMContext::fixed_point_value<Number>(unsigned int, const Point &, Number &, const Real) const;
2202template LIBMESH_EXPORT void FEMContext::fixed_point_value<Gradient>(unsigned int, const Point &, Gradient &, const Real) const;
2203
2204template LIBMESH_EXPORT void FEMContext::fixed_point_gradient<Gradient>(unsigned int, const Point &, Gradient &, const Real) const;
2205template LIBMESH_EXPORT void FEMContext::fixed_point_gradient<Tensor>(unsigned int, const Point &, Tensor &, const Real) const;
2206
2207#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
2208template LIBMESH_EXPORT void FEMContext::fixed_point_hessian<Tensor>(unsigned int, const Point &, Tensor &, const Real) const;
2209//FIXME: Not everything is implemented yet for second derivatives of RealGradients
2210//template LIBMESH_EXPORT void FEMContext::fixed_point_hessian<??>(unsigned int, const Point &, ??&) const;
2211#endif
2212
2213template LIBMESH_EXPORT void FEMContext::interior_rate<Number>(unsigned int, unsigned int, Number &) const;
2214template LIBMESH_EXPORT void FEMContext::interior_rate<Gradient>(unsigned int, unsigned int, Gradient &) const;
2215
2216template LIBMESH_EXPORT void FEMContext::interior_rate_gradient<Gradient>(unsigned int, unsigned int, Gradient &) const;
2217template LIBMESH_EXPORT void FEMContext::interior_rate_gradient<Tensor>(unsigned int, unsigned int, Tensor &) const;
2218
2219template LIBMESH_EXPORT void FEMContext::side_rate<Number>(unsigned int, unsigned int, Number &) const;
2220template LIBMESH_EXPORT void FEMContext::side_rate<Gradient>(unsigned int, unsigned int, Gradient &) const;
2221
2222template LIBMESH_EXPORT void FEMContext::interior_accel<Number>(unsigned int, unsigned int, Number &) const;
2223template LIBMESH_EXPORT void FEMContext::interior_accel<Gradient>(unsigned int, unsigned int, Gradient &) const;
2224
2225template LIBMESH_EXPORT void FEMContext::side_accel<Number>(unsigned int, unsigned int, Number &) const;
2226template LIBMESH_EXPORT void FEMContext::side_accel<Gradient>(unsigned int, unsigned int, Gradient &) const;
2227
2228template LIBMESH_EXPORT FEGenericBase<Real> *
2229FEMContext::build_new_fe(const FEGenericBase<Real>*,
2230 const Point &,
2231 const Real,
2232 const int) const;
2233
2234template LIBMESH_EXPORT FEGenericBase<RealGradient> *
2235FEMContext::build_new_fe(const FEGenericBase<RealGradient>*,
2236 const Point &,
2237 const Real,
2238 const int) const;
2239
2240} // namespace libMesh
dim
unsigned int dim
Definition adaptivity_ex3.C:114
libMesh::BoundaryInfo::has_boundary_id
bool has_boundary_id(const Node *const node, const boundary_id_type id) const
Definition boundary_info.C:1377
libMesh::BoundaryInfo::boundary_ids
void boundary_ids(const Node *node, std::vector< boundary_id_type > &vec_to_fill) const
Fills a user-provided std::vector with the boundary ids associated with Node node.
Definition boundary_info.C:1389
libMesh::DenseMatrix::resize
void resize(const unsigned int new_m, const unsigned int new_n)
Resizes the matrix to the specified size and calls zero().
Definition dense_matrix.h:895
libMesh::DenseSubVector
Defines a dense subvector for use in finite element computations.
Definition dense_subvector.h:44
libMesh::DenseSubVector::size
virtual unsigned int size() const override final
Definition dense_subvector.h:90
libMesh::DenseVector
Defines a dense vector for use in Finite Element-type computations.
Definition dense_vector.h:59
libMesh::DenseVector::resize
void resize(const unsigned int n)
Resize the vector.
Definition dense_vector.h:396
libMesh::DenseVector::get_values
std::vector< T > & get_values()
Definition dense_vector.h:278
libMesh::DiffContext
This class provides all data required for a physics package (e.g.
Definition diff_context.h:56
libMesh::DiffContext::_localized_vectors
std::map< const NumericVector< Number > *, std::pair< DenseVector< Number >, std::vector< DenseSubVector< Number > > > > _localized_vectors
Contains pointers to vectors the user has asked to be localized, keyed with pairs of element localize...
Definition diff_context.h:571
libMesh::DiffContext::get_elem_solution_rate
const DenseVector< Number > & get_elem_solution_rate() const
Accessor for element solution rate of change w.r.t.
Definition diff_context.h:144
libMesh::DiffContext::_elem_fixed_subsolutions
std::vector< DenseSubVector< Number > > _elem_fixed_subsolutions
Definition diff_context.h:605
libMesh::DiffContext::_elem_subsolutions
std::vector< DenseSubVector< Number > > _elem_subsolutions
Definition diff_context.h:583
libMesh::DiffContext::localized_vectors_iterator
std::map< constNumericVector< Number > *, std::pair< DenseVector< Number >, std::vector< DenseSubVector< Number > > > >::iterator localized_vectors_iterator
Typedef for the localized_vectors iterator.
Definition diff_context.h:540
libMesh::DiffContext::get_dof_indices
const std::vector< dof_id_type > & get_dof_indices() const
Accessor for element dof indices.
Definition diff_context.h:363
libMesh::DiffContext::get_system_time
Real get_system_time() const
Accessor for the time variable stored in the system class.
Definition diff_context.h:411
libMesh::DiffContext::_elem_qoi_subderivatives
std::vector< std::vector< DenseSubVector< Number > > > _elem_qoi_subderivatives
Definition diff_context.h:627
libMesh::DiffContext::get_elem_solution
const DenseVector< Number > & get_elem_solution() const
Accessor for element solution.
Definition diff_context.h:112
libMesh::DiffContext::get_elem_solution_accel
const DenseVector< Number > & get_elem_solution_accel() const
Accessor for element solution accel of change w.r.t.
Definition diff_context.h:177
libMesh::DiffContext::set_time
void set_time(Real time_in)
Set the time for which the current nonlinear_solution is defined.
Definition diff_context.h:423
libMesh::DiffContext::get_localized_subvector
DenseSubVector< Number > & get_localized_subvector(const NumericVector< Number > &localized_vector, unsigned int var)
Return a reference to DenseSubVector localization of localized_vector at variable var contained in th...
Definition diff_context.C:154
libMesh::DiffContext::get_system
const System & get_system() const
Accessor for associated system.
Definition diff_context.h:106
libMesh::DiffContext::get_elem_fixed_solution
const DenseVector< Number > & get_elem_fixed_solution() const
Accessor for element fixed solution.
Definition diff_context.h:210
libMesh::DiffContext::get_qoi_derivatives
const std::vector< DenseVector< Number > > & get_qoi_derivatives() const
Const accessor for QoI derivatives.
Definition diff_context.h:329
libMesh::DiffContext::_have_local_matrices
const bool _have_local_matrices
Whether we have local matrices allocated/initialized.
Definition diff_context.h:576
libMesh::DiffContext::get_elem_residual
const DenseVector< Number > & get_elem_residual() const
Const accessor for element residual.
Definition diff_context.h:242
libMesh::DiffContext::get_elem_jacobian
const DenseMatrix< Number > & get_elem_jacobian() const
Const accessor for element Jacobian.
Definition diff_context.h:274
libMesh::DifferentiablePhysics::get_second_order_vars
const std::set< unsigned int > & get_second_order_vars() const
Definition diff_physics.h:519
libMesh::DifferentiableSystem
This class provides a specific system class.
Definition diff_system.h:57
libMesh::DofMap::dof_indices
void dof_indices(const Elem *const elem, std::vector< dof_id_type > &di) const
Definition dof_map.C:2201
libMesh::DofMap::old_dof_indices
void old_dof_indices(const Elem &elem, unsigned int n, std::vector< dof_id_type > &di, const unsigned int vn) const
Appends to the vector di the old global degree of freedom indices for elem.node_ref(n),...
Definition dof_map.C:2478
libMesh::Elem
This is the base class from which all geometric element types are derived.
Definition elem.h:96
libMesh::Elem::point
const Point & point(const unsigned int i) const
Definition elem.h:2462
libMesh::Elem::n_nodes
virtual unsigned int n_nodes() const =0
libMesh::Elem::JUST_COARSENED
@ JUST_COARSENED
Definition elem.h:1450
libMesh::Elem::dim
virtual unsigned short dim() const =0
libMesh::Elem::infinite
virtual bool infinite() const =0
libMesh::FEAbstract::add_p_level_in_reinit
void add_p_level_in_reinit(bool value)
Indicate whether to add p-refinement levels in init/reinit methods.
Definition fe_abstract.h:631
libMesh::FEAbstract::reinit
virtual void reinit(const Elem *elem, const std::vector< Point > *const pts=nullptr, const std::vector< Real > *const weights=nullptr)=0
This is at the core of this class.
libMesh::FEAbstract::build
static std::unique_ptr< FEAbstract > build(const unsigned int dim, const FEType &type)
Builds a specific finite element type.
Definition fe_abstract.C:78
libMesh::FEAbstract::get_fe_type
FEType get_fe_type() const
Definition fe_abstract.h:520
libMesh::FEGenericBase
This class forms the foundation from which generic finite elements may be derived.
Definition fe_base.h:86
libMesh::FEGenericBase::build
static std::unique_ptr< FEGenericBase > build(const unsigned int dim, const FEType &type)
Builds a specific finite element type.
libMesh::FEGenericBase::get_div_phi
virtual_for_inffe const std::vector< std::vector< OutputDivergence > > & get_div_phi() const
Definition fe_base.h:261
libMesh::FEGenericBase::get_phi
const std::vector< std::vector< OutputShape > > & get_phi() const
Definition fe_base.h:207
libMesh::FEGenericBase::get_dphi
const std::vector< std::vector< OutputGradient > > & get_dphi() const
Definition fe_base.h:230
libMesh::FEGenericBase::build_InfFE
static std::unique_ptr< FEGenericBase > build_InfFE(const unsigned int dim, const FEType &type)
Builds a specific infinite element type.
libMesh::FEGenericBase::get_d2phi
const std::vector< std::vector< OutputTensor > > & get_d2phi() const
Definition fe_base.h:319
libMesh::FEGenericBase::get_curl_phi
virtual_for_inffe const std::vector< std::vector< OutputShape > > & get_curl_phi() const
Definition fe_base.h:252
libMesh::FEMContext
This class provides all data required for a physics package (e.g.
Definition fem_context.h:63
libMesh::FEMContext::_custom_solution
const NumericVector< Number > * _custom_solution
Data with which to do algebra reinitialization.
Definition fem_context.h:1074
libMesh::FEMContext::some_value
void some_value(unsigned int var, unsigned int qp, OutputType &u) const
Helper function to reduce some code duplication in the *interior_value methods.
Definition fem_context.C:314
libMesh::FEMContext::find_hardest_fe_type
FEType find_hardest_fe_type()
Helper function for creating quadrature rules.
Definition fem_context.C:86
libMesh::FEMContext::interior_values
void interior_values(unsigned int var, const NumericVector< Number > &_system_vector, std::vector< OutputType > &interior_values_vector) const
Fills a vector of values of the _system_vector at the all the quadrature points in the current elemen...
Definition fem_context.C:431
libMesh::FEMContext::interior_curl
void interior_curl(unsigned int var, unsigned int qp, OutputType &curl_u) const
Definition fem_context.C:597
libMesh::FEMContext::side_fe_reinit
virtual void side_fe_reinit()
Reinitializes side FE objects on the current geometric element.
Definition fem_context.C:1498
libMesh::FEMContext::get_side_fe
void get_side_fe(unsigned int var, FEGenericBase< OutputShape > *&fe) const
Accessor for edge/face (2D/3D) finite element object for variable var for the largest dimension in th...
Definition fem_context.h:317
libMesh::FEMContext::_do_elem_position_set
void _do_elem_position_set(Real theta)
Uses the coordinate data specified by mesh_*_position configuration to set the geometry of elem to th...
Definition fem_context.C:1607
libMesh::FEMContext::interior_gradient
Gradient interior_gradient(unsigned int var, unsigned int qp) const
Definition fem_context.C:466
libMesh::FEMContext::interior_rate_gradient
void interior_rate_gradient(unsigned int var, unsigned int qp, OutputType &u) const
Definition fem_context.C:1368
libMesh::FEMContext::side_gradients
void side_gradients(unsigned int var, const NumericVector< Number > &_system_vector, std::vector< OutputType > &side_gradients_vector) const
Fills a vector with the gradient of the solution variable var at all the quadrature points on the cur...
Definition fem_context.C:763
libMesh::FEMContext::_edge_fe
std::map< FEType, std::unique_ptr< FEAbstract > > _edge_fe
Definition fem_context.h:1170
libMesh::FEMContext::_mesh_sys
System * _mesh_sys
System from which to acquire moving mesh information.
Definition fem_context.h:1059
libMesh::FEMContext::side_value
Number side_value(unsigned int var, unsigned int qp) const
Definition fem_context.C:660
libMesh::FEMContext::attach_quadrature_rules
void attach_quadrature_rules()
Helper function for attaching quadrature rules.
Definition fem_context.C:127
libMesh::FEMContext::some_hessian
void some_hessian(unsigned int var, unsigned int qp, OutputType &u) const
Helper function to reduce some code duplication in the *interior_hessian methods.
Definition fem_context.C:381
libMesh::FEMContext::interior_rate
void interior_rate(unsigned int var, unsigned int qp, OutputType &u) const
Definition fem_context.C:1358
libMesh::FEMContext::side_hessian
Tensor side_hessian(unsigned int var, unsigned int qp) const
Definition fem_context.C:801
libMesh::FEMContext::get_mesh_y_var
unsigned int get_mesh_y_var() const
Accessor for y-variable of moving mesh System.
Definition fem_context.h:874
libMesh::FEMContext::_element_fe_var
std::vector< std::vector< FEAbstract * > > _element_fe_var
Pointers to the same finite element objects, but indexed by variable number.
Definition fem_context.h:1179
libMesh::FEMContext::point_curl
void point_curl(unsigned int var, const Point &p, OutputType &curl_u, const Real tolerance=TOLERANCE) const
Definition fem_context.C:1022
libMesh::FEMContext::elem_reinit
virtual void elem_reinit(Real theta) override
Resets the current time in the context.
Definition fem_context.C:1412
libMesh::FEMContext::fixed_side_gradient
Gradient fixed_side_gradient(unsigned int var, unsigned int qp) const
Definition fem_context.C:1158
libMesh::FEMContext::fixed_point_value
Number fixed_point_value(unsigned int var, const Point &p) const
Definition fem_context.C:1208
libMesh::FEMContext::_element_fe
std::vector< std::map< FEType, std::unique_ptr< FEAbstract > > > _element_fe
Finite element objects for each variable's interior, sides and edges.
Definition fem_context.h:1168
libMesh::FEMContext::pre_fe_reinit
virtual void pre_fe_reinit(const System &, const Elem *e)
Reinitializes local data vectors/matrices on the current geometric element.
Definition fem_context.C:1683
libMesh::FEMContext::_side_qrule
std::vector< std::unique_ptr< QBase > > _side_qrule
Quadrature rules for element sides The FEM context will try to find a quadrature rule that correctly ...
Definition fem_context.h:1224
libMesh::FEMContext::side_values
void side_values(unsigned int var, const NumericVector< Number > &_system_vector, std::vector< OutputType > &side_values_vector) const
Fills a vector of values of the _system_vector at the all the quadrature points on the current elemen...
Definition fem_context.C:683
libMesh::FEMContext::cached_fe
FEGenericBase< OutputShape > * cached_fe(const unsigned int elem_dim, const FEType fe_type, const int get_derivative_level) const
libMesh::FEMContext::fixed_side_hessian
Tensor fixed_side_hessian(unsigned int var, unsigned int qp) const
Definition fem_context.C:1183
libMesh::FEMContext::side
unsigned char side
Current side for side_* to examine.
Definition fem_context.h:1013
libMesh::FEMContext::side_boundary_ids
void side_boundary_ids(std::vector< boundary_id_type > &vec_to_fill) const
As above, but fills in the std::set provided by the user.
Definition fem_context.C:304
libMesh::FEMContext::side_gradient
Gradient side_gradient(unsigned int var, unsigned int qp) const
Definition fem_context.C:718
libMesh::FEMContext::diff_subsolution_getter
const DenseSubVector< Number > &(DiffContext::* diff_subsolution_getter)(unsigned int) const
Helper typedef to simplify refactoring.
Definition fem_context.h:1103
libMesh::FEMContext::interior_value
Number interior_value(unsigned int var, unsigned int qp) const
Definition fem_context.C:411
libMesh::FEMContext::some_gradient
void some_gradient(unsigned int var, unsigned int qp, OutputType &u) const
Helper function to reduce some code duplication in the *interior_gradient methods.
Definition fem_context.C:348
libMesh::FEMContext::_side_fe_var
std::vector< std::vector< FEAbstract * > > _side_fe_var
Definition fem_context.h:1180
libMesh::FEMContext::_elem_dim
unsigned char _elem_dim
Cached dimension of this->_elem.
Definition fem_context.h:1202
libMesh::FEMContext::side_hessians
void side_hessians(unsigned int var, const NumericVector< Number > &_system_vector, std::vector< OutputType > &d2u_vals) const
Fills a vector of hessians of the _system_vector at the all the quadrature points on the current elem...
Definition fem_context.C:830
libMesh::FEMContext::get_elem
const Elem & get_elem() const
Accessor for current Elem object.
Definition fem_context.h:908
libMesh::FEMContext::elem_side_reinit
virtual void elem_side_reinit(Real theta) override
Resets the current time in the context.
Definition fem_context.C:1436
libMesh::FEMContext::side_accel
void side_accel(unsigned int var, unsigned int qp, OutputType &u) const
Definition fem_context.C:1401
libMesh::FEMContext::set_jacobian_tolerance
void set_jacobian_tolerance(Real tol)
Calls set_jacobian_tolerance() on all the FE objects controlled by this class.
Definition fem_context.C:1586
libMesh::FEMContext::elem_position_get
void elem_position_get()
Uses the geometry of elem to set the coordinate data specified by mesh_*_position configuration.
Definition fem_context.C:1526
libMesh::FEMContext::_real_grad_fe
std::unique_ptr< FEGenericBase< RealGradient > > _real_grad_fe
Definition fem_context.h:1077
libMesh::FEMContext::set_elem
void set_elem(const Elem *e)
Helper function to promote accessor usage.
Definition fem_context.C:1914
libMesh::FEMContext::side_rate
void side_rate(unsigned int var, unsigned int qp, OutputType &u) const
Definition fem_context.C:1379
libMesh::FEMContext::interior_hessians
void interior_hessians(unsigned int var, const NumericVector< Number > &_system_vector, std::vector< OutputType > &d2u_vals) const
Fills a vector of hessians of the _system_vector at the all the quadrature points in the current elem...
Definition fem_context.C:555
libMesh::FEMContext::elem_edge_reinit
virtual void elem_edge_reinit(Real theta) override
Resets the current time in the context.
Definition fem_context.C:1451
libMesh::FEMContext::FEMContext
FEMContext(const System &sys, const std::vector< unsigned int > *active_vars=nullptr, bool allocate_local_matrices=true)
Constructor.
Definition fem_context.C:39
libMesh::FEMContext::use_default_quadrature_rules
void use_default_quadrature_rules(int extra_quadrature_order=0)
Use quadrature rules designed to over-integrate a mass matrix, plus extra_quadrature_order.
Definition fem_context.C:155
libMesh::FEMContext::elem_position_set
void elem_position_set(Real theta)
Uses the coordinate data specified by mesh_*_position configuration to set the geometry of elem to th...
Definition fem_context.h:1270
libMesh::FEMContext::interior_gradients
void interior_gradients(unsigned int var, const NumericVector< Number > &_system_vector, std::vector< OutputType > &interior_gradients_vector) const
Fills a vector with the gradient of the solution variable var at all the quadrature points in the cur...
Definition fem_context.C:493
libMesh::FEMContext::has_elem
bool has_elem() const
Test for current Elem object.
Definition fem_context.h:902
libMesh::FEMContext::interior_div
void interior_div(unsigned int var, unsigned int qp, OutputType &div_u) const
Definition fem_context.C:628
libMesh::FEMContext::build_new_fe
FEGenericBase< OutputShape > * build_new_fe(const FEGenericBase< OutputShape > *fe, const Point &p, const Real tolerance=TOLERANCE, const int get_derivative_level=-1) const
Helper function to reduce some code duplication in the *_point_* methods.
Definition fem_context.C:2023
libMesh::FEMContext::_update_time_from_system
void _update_time_from_system(Real theta)
Update the time in the context object for the given value of theta, based on the values of "time" and...
Definition fem_context.C:1923
libMesh::FEMContext::_extra_quadrature_order
int _extra_quadrature_order
The extra quadrature order for this context.
Definition fem_context.h:1238
libMesh::FEMContext::interior_accel
void interior_accel(unsigned int var, unsigned int qp, OutputType &u) const
Definition fem_context.C:1389
libMesh::FEMContext::_elem
const Elem * _elem
Current element for element_* to examine.
Definition fem_context.h:1192
libMesh::FEMContext::get_elem_dim
unsigned char get_elem_dim() const
Definition fem_context.h:944
libMesh::FEMContext::_active_vars
std::unique_ptr< const std::vector< unsigned int > > _active_vars
Variables on which to enable calculations, or nullptr if all variables in the System are to be enable...
Definition fem_context.h:1054
libMesh::FEMContext::get_mesh_x_var
unsigned int get_mesh_x_var() const
Accessor for x-variable of moving mesh System.
Definition fem_context.h:860
libMesh::FEMContext::fixed_side_value
Number fixed_side_value(unsigned int var, unsigned int qp) const
Definition fem_context.C:1134
libMesh::FEMContext::active_vars
const std::vector< unsigned int > * active_vars() const
Return a pointer to the vector of active variables being computed for, or a null pointer if all varia...
Definition fem_context.h:1034
libMesh::FEMContext::elem_fe_reinit
virtual void elem_fe_reinit(const std::vector< Point > *const pts=nullptr)
Reinitializes interior FE objects on the current geometric element.
Definition fem_context.C:1476
libMesh::FEMContext::_elem_dims
std::set< unsigned char > _elem_dims
Cached dimensions of elements in the mesh, plus dimension 0 if SCALAR variables are in use.
Definition fem_context.h:1208
libMesh::FEMContext::get_side
unsigned char get_side() const
Accessor for current side of Elem object.
Definition fem_context.h:922
libMesh::FEMContext::get_element_fe
void get_element_fe(unsigned int var, FEGenericBase< OutputShape > *&fe) const
Accessor for interior finite element object for variable var for the largest dimension in the mesh.
Definition fem_context.h:277
libMesh::FEMContext::point_value
Number point_value(unsigned int var, const Point &p) const
Definition fem_context.C:873
libMesh::FEMContext::point_hessian
Tensor point_hessian(unsigned int var, const Point &p) const
Definition fem_context.C:971
libMesh::FEMContext::_edge_fe_var
std::vector< FEAbstract * > _edge_fe_var
Definition fem_context.h:1181
libMesh::FEMContext::~FEMContext
virtual ~FEMContext()
Destructor.
Definition fem_context.C:291
libMesh::FEMContext::point_gradient
Gradient point_gradient(unsigned int var, const Point &p) const
Definition fem_context.C:919
libMesh::FEMContext::fixed_interior_hessian
Tensor fixed_interior_hessian(unsigned int var, unsigned int qp) const
Definition fem_context.C:1108
libMesh::FEMContext::_side_fe
std::vector< std::map< FEType, std::unique_ptr< FEAbstract > > > _side_fe
Definition fem_context.h:1169
libMesh::FEMContext::fixed_point_hessian
Tensor fixed_point_hessian(unsigned int var, const Point &p) const
Definition fem_context.C:1305
libMesh::FEMContext::use_unweighted_quadrature_rules
void use_unweighted_quadrature_rules(int extra_quadrature_order=0)
Use quadrature rules designed to exactly integrate unweighted undistorted basis functions,...
Definition fem_context.C:177
libMesh::FEMContext::has_side_boundary_id
bool has_side_boundary_id(boundary_id_type id) const
Reports if the boundary id is found on the current side.
Definition fem_context.C:297
libMesh::FEMContext::get_edge
unsigned char get_edge() const
Accessor for current edge of Elem object.
Definition fem_context.h:928
libMesh::FEMContext::fixed_interior_value
Number fixed_interior_value(unsigned int var, unsigned int qp) const
Definition fem_context.C:1058
libMesh::FEMContext::_boundary_info
const BoundaryInfo & _boundary_info
Saved reference to BoundaryInfo on the mesh for this System.
Definition fem_context.h:1187
libMesh::FEMContext::_edge_qrule
std::unique_ptr< QBase > _edge_qrule
Quadrature rules for element edges.
Definition fem_context.h:1233
libMesh::FEMContext::fixed_point_gradient
Gradient fixed_point_gradient(unsigned int var, const Point &p) const
Definition fem_context.C:1254
libMesh::FEMContext::algebraic_type
AlgebraicType algebraic_type() const
Definition fem_context.h:992
libMesh::FEMContext::fixed_interior_gradient
Gradient fixed_interior_gradient(unsigned int var, unsigned int qp) const
Definition fem_context.C:1081
libMesh::FEMContext::get_mesh_z_var
unsigned int get_mesh_z_var() const
Accessor for z-variable of moving mesh System.
Definition fem_context.h:888
libMesh::FEMContext::nonlocal_reinit
virtual void nonlocal_reinit(Real theta) override
Gives derived classes the opportunity to reinitialize data needed for nonlocal calculations at a new ...
Definition fem_context.C:1466
libMesh::FEMContext::_element_qrule
std::vector< std::unique_ptr< QBase > > _element_qrule
Quadrature rule for element interior.
Definition fem_context.h:1216
libMesh::FEMContext::edge_fe_reinit
virtual void edge_fe_reinit()
Reinitializes edge FE objects on the current geometric element.
Definition fem_context.C:1514
libMesh::FEMContext::_real_fe
std::unique_ptr< FEGenericBase< Real > > _real_fe
Definition fem_context.h:1076
libMesh::FEMContext::init_internal_data
void init_internal_data(const System &sys)
Helper function used in constructors to set up internal data.
Definition fem_context.C:198
libMesh::FEMContext::interior_hessian
Tensor interior_hessian(unsigned int var, unsigned int qp) const
Definition fem_context.C:531
libMesh::FEMap::map_fe_type
static FEFamily map_fe_type(const Elem &elem)
Definition fe_map.C:46
libMesh::FEMap::inverse_map
static Point inverse_map(const unsigned int dim, const Elem *elem, const Point &p, const Real tolerance=TOLERANCE, const bool secure=true, const bool extra_checks=true)
Definition fe_map.C:1512
libMesh::FEType
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition fe_type.h:197
libMesh::FEType::default_quadrature_rule
std::unique_ptr< QBase > default_quadrature_rule(const unsigned int dim, const int extraorder=0) const
Definition fe_type.C:34
libMesh::FEType::unweighted_quadrature_rule
std::unique_ptr< QBase > unweighted_quadrature_rule(const unsigned int dim, const int extraorder=0) const
Definition fe_type.C:53
libMesh::FEType::order
OrderWrapper order
The approximation order of the element (at 0 p-refinement level).
Definition fe_type.h:203
libMesh::FEType::p_refinement
bool p_refinement
Whether or not the finite elements for this type increase their p refinement level on geometric eleme...
Definition fe_type.h:292
libMesh::FEType::family
FEFamily family
The type of finite element.
Definition fe_type.h:228
libMesh::NumericVector
Provides a uniform interface to vector storage schemes for different linear algebra libraries.
Definition numeric_vector.h:70
libMesh::NumericVector::get
virtual void get(const std::vector< numeric_index_type > &index, T *values) const
Access multiple components at once.
Definition numeric_vector.h:967
libMesh::Point
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition point.h:40
libMesh::System
Manages consistently variables, degrees of freedom, and coefficient vectors.
Definition system.h:100
libMesh::System::use_fixed_solution
bool use_fixed_solution
A boolean to be set to true by systems using elem_fixed_solution, for optional use by e....
Definition system.h:1625
libMesh::System::current_local_solution
std::unique_ptr< NumericVector< Number > > current_local_solution
All the values I need to compute my contribution to the simulation at hand.
Definition system.h:1667
libMesh::System::variable_type
const FEType & variable_type(const unsigned int i) const
Definition system.C:2721
libMesh::System::n_qois
unsigned int n_qois() const
Number of currently active quantities of interest.
Definition system.h:2562
libMesh::System::n_vars
unsigned int n_vars() const
Definition system.C:2674
libMesh::System::get_dof_map
const DofMap & get_dof_map() const
Definition system.h:2417
libMesh::TensorValue
This class defines a tensor in LIBMESH_DIM dimensional Real or Complex space.
Definition tensor_value.h:45
libMesh::TimeSolver::is_steady
virtual bool is_steady() const =0
Is this effectively a steady-state solver?
libMesh::UnsteadySolver
This is a generic class that defines a solver to handle time integration of DifferentiableSystems.
Definition unsteady_solver.h:51
libMesh::UnsteadySolver::time_order
virtual unsigned int time_order() const =0
libMesh::VectorValue
This class defines a vector in LIBMESH_DIM dimensional Real or Complex space.
Definition vector_value.h:48
libMesh
The libMesh namespace provides an interface to certain functionality in the library.
libMesh::index_range
auto index_range(const T &sizable)
Helper function that returns an IntRange<std::size_t> representing all the indices of the passed-in v...
Definition int_range.h:153
libMesh::boundary_id_type
int8_t boundary_id_type
Definition id_types.h:51
libMesh::libmesh_real
T libmesh_real(T a)
Definition libmesh_common.h:189
libMesh::libmesh_assert
libmesh_assert(ctx)
libMesh::invalid_uint
const unsigned int invalid_uint
A number which is used quite often to represent an invalid or uninitialized value for an unsigned int...
Definition libmesh.h:303
libMesh::Real
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Definition libmesh_common.h:144
libMesh::make_range
IntRange< T > make_range(T beg, T end)
The 2-parameter make_range() helper function returns an IntRange<T> when both input parameters are of...
Definition int_range.h:176
libMesh::n_threads
unsigned int n_threads()
Definition libmesh_base.h:109
libMesh::FEMContext::FENeeded
Helper nested class for C++03-compatible "template typedef".
Definition fem_context.h:1111
libMesh::FEMContext::FENeeded::value_shape
TensorTools::MakeReal< OutputType >::type value_shape
Definition fem_context.h:1117
n_nodes
const dof_id_type n_nodes
Definition tecplot_io.C:67
Generated on Mon Aug 24 2026 19:00:38 for libMesh by doxygen 1.9.8