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