Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2026 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 :
4 : // This library is free software; you can redistribute it and/or
5 : // modify it under the terms of the GNU Lesser General Public
6 : // License as published by the Free Software Foundation; either
7 : // version 2.1 of the License, or (at your option) any later version.
8 :
9 : // This library is distributed in the hope that it will be useful,
10 : // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 : // Lesser General Public License for more details.
13 :
14 : // You should have received a copy of the GNU Lesser General Public
15 : // License along with this library; if not, write to the Free Software
16 : // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 :
18 :
19 :
20 : #ifndef LIBMESH_FEM_CONTEXT_H
21 : #define LIBMESH_FEM_CONTEXT_H
22 :
23 : // Local Includes
24 : #include "libmesh/diff_context.h"
25 : #include "libmesh/id_types.h"
26 : #include "libmesh/enum_quadrature_type.h"
27 : #include "libmesh/fe_type.h"
28 : #include "libmesh/fe_base.h"
29 : #include "libmesh/vector_value.h"
30 :
31 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
32 : #include "libmesh/tensor_value.h"
33 : #endif
34 :
35 : // C++ includes
36 : #include <map>
37 : #include <set>
38 :
39 : namespace libMesh
40 : {
41 :
42 : // Forward Declarations
43 : class BoundaryInfo;
44 : class Elem;
45 : template <typename T> class FEGenericBase;
46 : typedef FEGenericBase<Real> FEBase;
47 : class QBase;
48 : class Point;
49 : template <typename T> class NumericVector;
50 :
51 : /**
52 : * This class provides all data required for a physics package
53 : * (e.g. an FEMSystem subclass) to perform local element residual
54 : * and jacobian integrations.
55 : *
56 : * This class is part of the new DifferentiableSystem framework,
57 : * which is still experimental. Users of this framework should
58 : * beware of bugs and future API changes.
59 : *
60 : * \author Roy H. Stogner
61 : * \date 2009
62 : */
63 : class FEMContext : public DiffContext
64 : {
65 : public:
66 :
67 : /**
68 : * Constructor. Allocates some but fills no data structures.
69 : *
70 : * Optionally specify a limited number of variables to be "active"
71 : * and thus calculated on. If \p active_vars is null then
72 : * calculations will be prepared for every variable in \p sys.
73 : */
74 : explicit
75 : FEMContext (const System & sys,
76 : const std::vector<unsigned int> * active_vars = nullptr,
77 : bool allocate_local_matrices = true);
78 :
79 : /**
80 : * Constructor. Specify the extra quadrature order instead
81 : * of getting it from \p sys.
82 : *
83 : * Optionally specify a limited number of variables to be "active"
84 : * and thus calculated on. If \p active_vars is null then
85 : * calculations will be prepared for every variable in \p sys.
86 : */
87 : explicit
88 : FEMContext (const System & sys,
89 : int extra_quadrature_order,
90 : const std::vector<unsigned int> * active_vars = nullptr,
91 : bool allocate_local_matrices = true);
92 :
93 :
94 : /**
95 : * Destructor.
96 : */
97 : virtual ~FEMContext ();
98 :
99 : /**
100 : * Use quadrature rules designed to over-integrate a mass matrix,
101 : * plus \p extra_quadrature_order.
102 : */
103 : void use_default_quadrature_rules(int extra_quadrature_order=0);
104 :
105 : /**
106 : * Use quadrature rules designed to exactly integrate unweighted
107 : * undistorted basis functions, plus \p extra_quadrature_order.
108 : */
109 : void use_unweighted_quadrature_rules(int extra_quadrature_order=0);
110 :
111 : /**
112 : * Use a specified quadrature rule type \p qt for every element,
113 : * side, and edge dimension, at the default quadrature order for the
114 : * hardest FE type plus \p extra_quadrature_order.
115 : *
116 : * This is useful when the default Gauss rules (whose points lie in
117 : * the element interior) are inadequate. For example, rules whose
118 : * points include the element vertices (\p QTRAP, \p QSIMPSON, \p
119 : * QNODAL, \p QGAUSS_LOBATTO) can detect integrand behavior that is
120 : * localized at element corners --- a feature interior Gauss points
121 : * can miss entirely. Note that not every quadrature type is defined
122 : * for every element type (e.g. \p QGAUSS_LOBATTO is only available for
123 : * tensor-product elements); it is the caller's responsibility to pick
124 : * a type compatible with the mesh.
125 : */
126 : void use_quadrature_rules(QuadratureType qt, int extra_quadrature_order=0);
127 :
128 : /**
129 : * Reports if the boundary id is found on the current side
130 : */
131 : bool has_side_boundary_id(boundary_id_type id) const;
132 :
133 : /**
134 : * As above, but fills in the std::set provided by the user.
135 : */
136 : void side_boundary_ids(std::vector<boundary_id_type> & vec_to_fill) const;
137 :
138 : /**
139 : * \returns The value of the solution variable \p var at the
140 : * quadrature point \p qp on the current element interior.
141 : *
142 : * \note This API is currently present for backward compatibility.
143 : */
144 : Number interior_value(unsigned int var, unsigned int qp) const;
145 :
146 : /**
147 : * \returns The value of the solution variable \p var at the quadrature
148 : * point \p qp on the current element side.
149 : *
150 : * \note This API currently is present for backward compatibility.
151 : */
152 : Number side_value(unsigned int var, unsigned int qp) const;
153 :
154 : /**
155 : * \returns The value of the solution variable \p var at the physical
156 : * point \p p on the current element.
157 : *
158 : * \note This API is currently present for backward compatibility.
159 : */
160 : Number point_value(unsigned int var, const Point & p) const;
161 :
162 : /**
163 : * \returns The gradient of the solution variable \p var at the quadrature
164 : * point \p qp on the current element interior.
165 : *
166 : * \note This API is currently present for backward compatibility.
167 : */
168 : Gradient interior_gradient(unsigned int var, unsigned int qp) const;
169 :
170 : /**
171 : * \returns The gradient of the solution variable \p var at the quadrature
172 : * point \p qp on the current element side.
173 : *
174 : * \note This API is currently present for backward compatibility.
175 : */
176 : Gradient side_gradient(unsigned int var, unsigned int qp) const;
177 :
178 : /**
179 : * \returns The gradient of the solution variable \p var at the physical
180 : * point \p p on the current element.
181 : *
182 : * \note This API is currently present for backward compatibility.
183 : */
184 : Gradient point_gradient(unsigned int var, const Point & p) const;
185 :
186 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
187 : /**
188 : * \returns The hessian of the solution variable \p var at the quadrature
189 : * point \p qp on the current element interior.
190 : *
191 : * \note This API is currently present for backward compatibility.
192 : */
193 : Tensor interior_hessian(unsigned int var, unsigned int qp) const;
194 :
195 : /**
196 : * \returns The hessian of the solution variable \p var at the quadrature
197 : * point \p qp on the current element side.
198 : *
199 : * \note This API is currently present for backward compatibility.
200 : */
201 : Tensor side_hessian(unsigned int var, unsigned int qp) const;
202 :
203 : /**
204 : * \returns The hessian of the solution variable \p var at the physical
205 : * point \p p on the current element.
206 : *
207 : * \note This API currently present for backward compatibility.
208 : */
209 : Tensor point_hessian(unsigned int var, const Point & p) const;
210 :
211 : #endif // LIBMESH_ENABLE_SECOND_DERIVATIVES
212 :
213 : /**
214 : * \returns The value of the fixed_solution variable \p var at the quadrature
215 : * point \p qp on the current element interior.
216 : *
217 : * \note This API is currently present for backward compatibility.
218 : */
219 : Number fixed_interior_value(unsigned int var, unsigned int qp) const;
220 :
221 : /**
222 : * \returns The value of the fixed_solution variable \p var at the quadrature
223 : * point \p qp on the current element side.
224 : *
225 : * \note This API is currently present for backward compatibility.
226 : */
227 : Number fixed_side_value(unsigned int var, unsigned int qp) const;
228 :
229 : /**
230 : * \returns The value of the fixed_solution variable \p var at the physical
231 : * point \p p on the current element.
232 : *
233 : * \note This API is currently present for backward compatibility.
234 : */
235 : Number fixed_point_value(unsigned int var, const Point & p) const;
236 :
237 : /**
238 : * \returns The gradient of the fixed_solution variable \p var at the quadrature
239 : * point \p qp on the current element interior.
240 : *
241 : * \note This API is currently present for backward compatibility.
242 : */
243 : Gradient fixed_interior_gradient(unsigned int var, unsigned int qp) const;
244 :
245 : /**
246 : * \returns The gradient of the fixed_solution variable \p var at the quadrature
247 : * point \p qp on the current element side.
248 : *
249 : * \note This API is currently present for backward compatibility.
250 : */
251 : Gradient fixed_side_gradient(unsigned int var, unsigned int qp) const;
252 :
253 : /**
254 : * \returns The gradient of the fixed_solution variable \p var at the physical
255 : * point \p p on the current element.
256 : *
257 : * \note This API is currently present for backward compatibility.
258 : */
259 : Gradient fixed_point_gradient(unsigned int var, const Point & p) const;
260 :
261 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
262 : /**
263 : * \returns The hessian of the fixed_solution variable \p var at the quadrature
264 : * point \p qp on the current element interior.
265 : *
266 : * \note This API is currently present for backward compatibility.
267 : */
268 : Tensor fixed_interior_hessian(unsigned int var, unsigned int qp) const;
269 :
270 : /**
271 : * \returns The hessian of the fixed_solution variable \p var at the quadrature
272 : * point \p qp on the current element side.
273 : *
274 : * \note This API is currently present for backward compatibility.
275 : */
276 : Tensor fixed_side_hessian(unsigned int var, unsigned int qp) const;
277 :
278 : /**
279 : * \returns The hessian of the fixed_solution variable \p var at the physical
280 : * point \p p on the current element.
281 : *
282 : * \note This API is currently present for backward compatibility.
283 : */
284 : Tensor fixed_point_hessian (unsigned int var, const Point & p) const;
285 :
286 : #endif // LIBMESH_ENABLE_SECOND_DERIVATIVES
287 :
288 : /**
289 : * Accessor for interior finite element object for variable var for
290 : * the largest dimension in the mesh. We default to the largest mesh dim
291 : * if this method is called before the Elem * is set in the FEMContext,
292 : * e.g. in FEMSystem::init_context (or a subclass).
293 : */
294 : template<typename OutputShape>
295 15589464 : void get_element_fe( unsigned int var, FEGenericBase<OutputShape> *& fe ) const
296 15589464 : { this->get_element_fe<OutputShape>(var,fe,this->get_elem_dim()); }
297 :
298 : /**
299 : * Accessor for interior finite element object for scalar-valued variable var
300 : * for the largest dimension in the mesh. We default to the largest mesh dim
301 : * if this method is called before the Elem * is set in the FEMContext,
302 : * e.g. in FEMSystem::init_context (or a subclass).
303 : */
304 0 : FEBase * get_element_fe( unsigned int var ) const
305 0 : { return this->get_element_fe(var,this->get_elem_dim()); }
306 :
307 : /**
308 : * Accessor for interior finite element object for variable var for
309 : * dimension dim.
310 : */
311 : template<typename OutputShape>
312 : void get_element_fe( unsigned int var, FEGenericBase<OutputShape> *& fe,
313 : unsigned short dim ) const;
314 :
315 : /**
316 : * Accessor for interior finite element object for variable var for
317 : * dimension dim.
318 : */
319 : void get_element_fe( unsigned int var, FEAbstract *& fe,
320 : unsigned short dim ) const;
321 :
322 : /**
323 : * Accessor for interior finite element object for scalar-valued variable var for
324 : * dimension dim.
325 : */
326 : FEBase * get_element_fe( unsigned int var, unsigned short dim ) const;
327 :
328 : /**
329 : * Accessor for edge/face (2D/3D) finite element object for variable var
330 : * for the largest dimension in the mesh. We default to the largest mesh dim
331 : * if this method is called before the Elem * is set in the FEMContext,
332 : * e.g. in FEMSystem::init_context (or a subclass).
333 : */
334 : template<typename OutputShape>
335 569205 : void get_side_fe( unsigned int var, FEGenericBase<OutputShape> *& fe ) const
336 569205 : { this->get_side_fe<OutputShape>(var,fe,this->get_elem_dim()); }
337 :
338 : /**
339 : * Accessor for side finite element object for scalar-valued variable var
340 : * for the largest dimension in the mesh. We default to the largest mesh dim
341 : * if this method is called before the Elem * is set in the FEMContext,
342 : * e.g. in FEMSystem::init_context (or a subclass).
343 : */
344 0 : FEBase * get_side_fe( unsigned int var ) const
345 0 : { return this->get_side_fe(var,this->get_elem_dim()); }
346 :
347 : /**
348 : * Accessor for edge/face (2D/3D) finite element object for variable var
349 : * for dimension dim.
350 : */
351 : template<typename OutputShape>
352 : void get_side_fe( unsigned int var, FEGenericBase<OutputShape> *& fe,
353 : unsigned short dim ) const;
354 :
355 : /**
356 : * Accessor for edge/face (2D/3D) finite element object for variable var
357 : * for dimension dim.
358 : */
359 : void get_side_fe( unsigned int var, FEAbstract *& fe,
360 : unsigned short dim ) const;
361 :
362 : /**
363 : * Accessor for side finite element object for scalar-valued variable var
364 : * for dimension dim.
365 : */
366 : FEBase * get_side_fe( unsigned int var, unsigned short dim ) const;
367 :
368 : /**
369 : * Accessor for edge (3D only!) finite element object for variable var.
370 : */
371 : template<typename OutputShape>
372 : void get_edge_fe( unsigned int var, FEGenericBase<OutputShape> *& fe ) const;
373 :
374 : void get_edge_fe( unsigned int var, FEAbstract *& fe ) const;
375 :
376 : /**
377 : * Accessor for edge (3D only!) finite element object for scalar-valued variable var.
378 : */
379 : FEBase * get_edge_fe( unsigned int var ) const;
380 :
381 : /**
382 : * \returns The value of the solution variable \p var at the quadrature
383 : * point \p qp on the current element interior.
384 : *
385 : * \note This is the preferred API.
386 : */
387 : template<typename OutputType>
388 : void interior_value(unsigned int var,
389 : unsigned int qp,
390 : OutputType & u) const;
391 :
392 : /**
393 : * Fills a vector of values of the _system_vector at the all the quadrature
394 : * points in the current element interior.
395 : */
396 : template<typename OutputType>
397 : void interior_values(unsigned int var,
398 : const NumericVector<Number> & _system_vector,
399 : std::vector<OutputType> & interior_values_vector) const;
400 :
401 : /**
402 : * \returns The value of the solution variable \p var at the quadrature
403 : * point \p qp on the current element side.
404 : *
405 : * \note This is the preferred API.
406 : */
407 : template<typename OutputType>
408 : void side_value(unsigned int var,
409 : unsigned int qp,
410 : OutputType & u) const;
411 :
412 : /**
413 : * Fills a vector of values of the _system_vector at the all the quadrature
414 : * points on the current element side.
415 : */
416 : template<typename OutputType>
417 : void side_values(unsigned int var,
418 : const NumericVector<Number> & _system_vector,
419 : std::vector<OutputType> & side_values_vector) const;
420 :
421 : /**
422 : * \returns The value of the solution variable \p var at the physical
423 : * point \p p on the current element.
424 : *
425 : * \note This is the preferred API.
426 : *
427 : * Allows evaluation of points within a relative tolerance outside
428 : * the element.
429 : */
430 : template<typename OutputType>
431 : void point_value(unsigned int var,
432 : const Point & p,
433 : OutputType & u,
434 : const Real tolerance = TOLERANCE) const;
435 :
436 : /**
437 : * \returns The gradient of the solution variable \p var at the quadrature
438 : * point \p qp on the current element interior.
439 : *
440 : * \note This is the preferred API.
441 : */
442 : template<typename OutputType>
443 : void interior_gradient(unsigned int var, unsigned int qp,
444 : OutputType & du) const;
445 :
446 : /**
447 : * Fills a vector with the gradient of the solution variable \p var at all the quadrature
448 : * points in the current element interior.
449 : *
450 : * \note This is the preferred API.
451 : */
452 : template<typename OutputType>
453 : void interior_gradients(unsigned int var,
454 : const NumericVector<Number> & _system_vector,
455 : std::vector<OutputType> & interior_gradients_vector) const;
456 :
457 : /**
458 : * \returns The gradient of the solution variable \p var at the quadrature
459 : * point \p qp on the current element side.
460 : *
461 : * \note This is the preferred API.
462 : */
463 : template<typename OutputType>
464 : void side_gradient(unsigned int var,
465 : unsigned int qp,
466 : OutputType & du) const;
467 :
468 : /**
469 : * Fills a vector with the gradient of the solution variable \p var at all the quadrature
470 : * points on the current element side.
471 : *
472 : * \note This is the preferred API.
473 : */
474 : template<typename OutputType>
475 : void side_gradients(unsigned int var,
476 : const NumericVector<Number> & _system_vector,
477 : std::vector<OutputType> & side_gradients_vector) const;
478 :
479 : /**
480 : * \returns The gradient of the solution variable \p var at the physical
481 : * point \p p on the current element.
482 : *
483 : * \note This is the preferred API.
484 : *
485 : * Allows evaluation of points within a relative tolerance outside
486 : * the element.
487 : */
488 : template<typename OutputType>
489 : void point_gradient(unsigned int var,
490 : const Point & p,
491 : OutputType & grad_u,
492 : const Real tolerance = TOLERANCE) const;
493 :
494 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
495 : /**
496 : * \returns The hessian of the solution variable \p var at the quadrature
497 : * point \p qp on the current element interior.
498 : *
499 : * \note This is the preferred API.
500 : */
501 : template<typename OutputType>
502 : void interior_hessian(unsigned int var,
503 : unsigned int qp,
504 : OutputType & d2u) const;
505 :
506 : /**
507 : * Fills a vector of hessians of the _system_vector at the all the
508 : * quadrature points in the current element interior. This is the
509 : * preferred API.
510 : */
511 : template<typename OutputType>
512 : void interior_hessians(unsigned int var,
513 : const NumericVector<Number> & _system_vector,
514 : std::vector<OutputType> & d2u_vals) const;
515 :
516 : /**
517 : * \returns The hessian of the solution variable \p var at the quadrature
518 : * point \p qp on the current element side.
519 : *
520 : * \note This is the preferred API.
521 : */
522 : template<typename OutputType>
523 : void side_hessian(unsigned int var,
524 : unsigned int qp,
525 : OutputType & d2u) const;
526 :
527 : /**
528 : * Fills a vector of hessians of the _system_vector at the all the
529 : * quadrature points on the current element side. This is the
530 : * preferred API.
531 : */
532 : template<typename OutputType>
533 : void side_hessians(unsigned int var,
534 : const NumericVector<Number> & _system_vector,
535 : std::vector<OutputType> & d2u_vals) const;
536 :
537 : /**
538 : * \returns The hessian of the solution variable \p var at the physical
539 : * point \p p on the current element.
540 : *
541 : * \note This is the preferred API.
542 : *
543 : * Allows evaluation of points within a relative tolerance outside
544 : * the element.
545 : */
546 : template<typename OutputType>
547 : void point_hessian(unsigned int var,
548 : const Point & p,
549 : OutputType & hess_u,
550 : const Real tolerance = TOLERANCE) const;
551 :
552 : #endif // LIBMESH_ENABLE_SECOND_DERIVATIVES
553 :
554 : /**
555 : * \returns The time derivative (rate) of the solution variable
556 : * \p var at the quadrature point \p qp on the current element
557 : * interior.
558 : */
559 : template<typename OutputType>
560 : void interior_rate(unsigned int var,
561 : unsigned int qp,
562 : OutputType & u) const;
563 :
564 :
565 : /**
566 : * \returns The time derivative (rate) of the solution gradient
567 : * of variable \p var at the quadrature point \p qp on the current
568 : * element interior.
569 : */
570 : template<typename OutputType>
571 : void interior_rate_gradient(unsigned int var,
572 : unsigned int qp,
573 : OutputType & u) const;
574 :
575 :
576 : /**
577 : * \returns The time derivative (rate) of the solution variable
578 : * \p var at the quadrature point \p qp on the current element side.
579 : */
580 : template<typename OutputType>
581 : void side_rate(unsigned int var,
582 : unsigned int qp,
583 : OutputType & u) const;
584 :
585 : /**
586 : * \returns The time derivative (rate) of the solution variable
587 : * \p var at the physical point \p p on the current element.
588 : */
589 : template<typename OutputType>
590 : void point_rate(unsigned int var,
591 : const Point & p,
592 : OutputType & u) const;
593 :
594 : /**
595 : * \returns The second time derivative (acceleration) of the solution variable
596 : * \p var at the quadrature point \p qp on the current element
597 : * interior.
598 : */
599 : template<typename OutputType>
600 : void interior_accel(unsigned int var,
601 : unsigned int qp,
602 : OutputType & u) const;
603 :
604 : /**
605 : * \returns The second time derivative (acceleration) of the solution variable
606 : * \p var at the quadrature point \p qp on the current element side.
607 : */
608 : template<typename OutputType>
609 : void side_accel(unsigned int var,
610 : unsigned int qp,
611 : OutputType & u) const;
612 :
613 : /**
614 : * \returns The second time derivative (acceleration) of the solution variable
615 : * \p var at the physical point \p p on the current element.
616 : */
617 : template<typename OutputType>
618 : void point_accel(unsigned int var,
619 : const Point & p,
620 : OutputType & u) const;
621 :
622 : /**
623 : * \returns The value of the fixed_solution variable \p var at the quadrature
624 : * point \p qp on the current element interior.
625 : *
626 : * \note This is the preferred API.
627 : */
628 : template<typename OutputType>
629 : void fixed_interior_value(unsigned int var,
630 : unsigned int qp,
631 : OutputType & u) const;
632 :
633 : /**
634 : * \returns The value of the fixed_solution variable \p var at the quadrature
635 : * point \p qp on the current element side.
636 : *
637 : * \note This is the preferred API.
638 : */
639 : template<typename OutputType>
640 : void fixed_side_value(unsigned int var,
641 : unsigned int qp,
642 : OutputType & u) const;
643 :
644 : /**
645 : * \returns The value of the fixed_solution variable \p var at the physical
646 : * point \p p on the current element.
647 : *
648 : * \note This is the preferred API.
649 : *
650 : * Allows evaluation of points within a relative tolerance outside
651 : * the element.
652 : */
653 : template<typename OutputType>
654 : void fixed_point_value(unsigned int var,
655 : const Point & p,
656 : OutputType & u,
657 : const Real tolerance = TOLERANCE) const;
658 :
659 : /**
660 : * \returns The gradient of the fixed_solution variable \p var at the quadrature
661 : * point \p qp on the current element interior.
662 : *
663 : * \note This is the preferred API.
664 : */
665 : template<typename OutputType>
666 : void fixed_interior_gradient(unsigned int var,
667 : unsigned int qp,
668 : OutputType & grad_u) const;
669 :
670 : /**
671 : * \returns The gradient of the fixed_solution variable \p var at the quadrature
672 : * point \p qp on the current element side.
673 : *
674 : * \note This is the preferred API.
675 : */
676 : template<typename OutputType>
677 : void fixed_side_gradient(unsigned int var,
678 : unsigned int qp,
679 : OutputType & grad_u) const;
680 :
681 : /**
682 : * \returns The gradient of the fixed_solution variable \p var at the physical
683 : * point \p p on the current element.
684 : *
685 : * \note This is the preferred API.
686 : *
687 : * Allows evaluation of points within a relative tolerance outside
688 : * the element.
689 : */
690 : template<typename OutputType>
691 : void fixed_point_gradient(unsigned int var,
692 : const Point & p,
693 : OutputType & grad_u,
694 : const Real tolerance = TOLERANCE) const;
695 :
696 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
697 : /**
698 : * \returns The hessian of the fixed_solution variable \p var at the quadrature
699 : * point \p qp on the current element interior.
700 : *
701 : * \note This is the preferred API.
702 : */
703 : template<typename OutputType>
704 : void fixed_interior_hessian(unsigned int var,
705 : unsigned int qp,
706 : OutputType & hess_u) const;
707 :
708 : /**
709 : * \returns The hessian of the fixed_solution variable \p var at the quadrature
710 : * point \p qp on the current element side.
711 : *
712 : * \note This is the preferred API.
713 : */
714 : template<typename OutputType>
715 : void fixed_side_hessian(unsigned int var,
716 : unsigned int qp,
717 : OutputType & hess_u) const;
718 :
719 : /**
720 : * \returns The hessian of the fixed_solution variable \p var at the physical
721 : * point \p p on the current element.
722 : *
723 : * \note This is the preferred API.
724 : *
725 : * Allows evaluation of points within a relative tolerance outside
726 : * the element.
727 : */
728 : template<typename OutputType>
729 : void fixed_point_hessian(unsigned int var,
730 : const Point & p,
731 : OutputType & hess_u,
732 : const Real tolerance = TOLERANCE) const;
733 :
734 : #endif // LIBMESH_ENABLE_SECOND_DERIVATIVES
735 :
736 : /**
737 : * \returns The curl of the solution variable \p var at the physical
738 : * point \p p on the current element.
739 : */
740 : template<typename OutputType>
741 : void interior_curl(unsigned int var,
742 : unsigned int qp,
743 : OutputType & curl_u) const;
744 :
745 : /**
746 : * \returns The curl of the solution variable \p var at the physical
747 : * point \p p on the current element.
748 : *
749 : * Allows evaluation of points within a relative tolerance outside
750 : * the element.
751 : */
752 : template<typename OutputType>
753 : void point_curl(unsigned int var,
754 : const Point & p,
755 : OutputType & curl_u,
756 : const Real tolerance = TOLERANCE) const;
757 :
758 : /**
759 : * \returns The divergence of the solution variable \p var at the physical
760 : * point \p p on the current element.
761 : */
762 : template<typename OutputType>
763 : void interior_div(unsigned int var,
764 : unsigned int qp,
765 : OutputType & div_u) const;
766 :
767 : // should be protected:
768 : /**
769 : * Resets the current time in the context. Additionally, reinitialize Elem
770 : * and FE objects if there's a moving mesh present in the system such that
771 : * the mesh is deformed to its position at \f$ t_{\theta} \f$.
772 : */
773 : virtual void elem_reinit(Real theta) override;
774 :
775 : /**
776 : * Resets the current time in the context. Additionally, reinitialize Elem
777 : * and FE objects if there's a moving mesh present in the system such that
778 : * the mesh is deformed to its position at \f$ t_{\theta} \f$.
779 : */
780 : virtual void elem_side_reinit(Real theta) override;
781 :
782 : /**
783 : * Resets the current time in the context. Additionally, reinitialize Elem
784 : * and FE objects if there's a moving mesh present in the system such that
785 : * the mesh is deformed to its position at \f$ t_{\theta} \f$.
786 : */
787 : virtual void elem_edge_reinit(Real theta) override;
788 :
789 : /**
790 : * Gives derived classes the opportunity to reinitialize data needed
791 : * for nonlocal calculations at a new point within a timestep
792 : */
793 : virtual void nonlocal_reinit(Real theta) override;
794 :
795 : /**
796 : * Reinitializes local data vectors/matrices on the current geometric element
797 : */
798 : virtual void pre_fe_reinit(const System &,
799 : const Elem * e);
800 :
801 : /**
802 : * Reinitializes interior FE objects on the current geometric element
803 : */
804 : virtual void elem_fe_reinit(const std::vector<Point> * const pts = nullptr);
805 :
806 : /**
807 : * Reinitializes side FE objects on the current geometric element
808 : */
809 : virtual void side_fe_reinit();
810 :
811 : /**
812 : * Reinitializes edge FE objects on the current geometric element
813 : */
814 : virtual void edge_fe_reinit();
815 :
816 : /**
817 : * Accessor for element interior quadrature rule for the dimension of the
818 : * current _elem.
819 : */
820 0 : const QBase & get_element_qrule() const
821 0 : { return this->get_element_qrule(this->get_elem_dim()); }
822 :
823 : /**
824 : * Accessor for element side quadrature rule for the dimension of the
825 : * current _elem.
826 : */
827 743 : const QBase & get_side_qrule() const
828 743 : { return this->get_side_qrule(this->get_elem_dim()); }
829 :
830 : /**
831 : * Accessor for element interior quadrature rule.
832 : */
833 0 : const QBase & get_element_qrule( unsigned short dim ) const
834 0 : { libmesh_assert(_element_qrule[dim]);
835 16201782 : return *(this->_element_qrule[dim]); }
836 :
837 : /**
838 : * Accessor for element side quadrature rule.
839 : */
840 743 : const QBase & get_side_qrule( unsigned short dim ) const
841 : {
842 743 : libmesh_assert(_side_qrule[dim]);
843 11004 : return *(this->_side_qrule[dim]);
844 : }
845 :
846 : /**
847 : * Accessor for element edge quadrature rule.
848 : */
849 0 : const QBase & get_edge_qrule() const
850 0 : { return *(this->_edge_qrule); }
851 :
852 : /**
853 : * Tells the FEMContext that system \p sys contains the
854 : * isoparametric Lagrangian variables which correspond to the
855 : * coordinates of mesh nodes, in problems where the mesh itself is
856 : * expected to move in time.
857 : *
858 : * This should be set automatically if the FEMPhysics requires it.
859 : */
860 320500 : virtual void set_mesh_system(System * sys)
861 320500 : { this->_mesh_sys = sys; }
862 :
863 : /**
864 : * Accessor for moving mesh System
865 : */
866 : const System * get_mesh_system() const
867 : { return this->_mesh_sys; }
868 :
869 : /**
870 : * Accessor for moving mesh System
871 : */
872 : System * get_mesh_system()
873 : { return this->_mesh_sys; }
874 :
875 : /**
876 : * Accessor for x-variable of moving mesh System
877 : */
878 0 : unsigned int get_mesh_x_var() const
879 23616 : { return _mesh_x_var; }
880 :
881 : /**
882 : * Accessor for x-variable of moving mesh System
883 : *
884 : * This should be set automatically if the FEMPhysics requires it.
885 : */
886 10450 : void set_mesh_x_var(unsigned int x_var)
887 320500 : { _mesh_x_var = x_var; }
888 :
889 : /**
890 : * Accessor for y-variable of moving mesh System
891 : */
892 0 : unsigned int get_mesh_y_var() const
893 23616 : { return _mesh_y_var; }
894 :
895 : /**
896 : * Accessor for y-variable of moving mesh System
897 : *
898 : * This should be set automatically if the FEMPhysics requires it.
899 : */
900 10450 : void set_mesh_y_var(unsigned int y_var)
901 320500 : { _mesh_y_var = y_var; }
902 :
903 : /**
904 : * Accessor for z-variable of moving mesh System
905 : */
906 0 : unsigned int get_mesh_z_var() const
907 23616 : { return _mesh_z_var; }
908 :
909 : /**
910 : * Accessor for z-variable of moving mesh System
911 : *
912 : * This should be set automatically if the FEMPhysics requires it.
913 : */
914 10450 : void set_mesh_z_var(unsigned int z_var)
915 320500 : { _mesh_z_var = z_var; }
916 :
917 : /**
918 : * Test for current Elem object
919 : */
920 18753172 : bool has_elem() const
921 181100200 : { return (this->_elem != nullptr); }
922 :
923 : /**
924 : * Accessor for current Elem object
925 : */
926 5764112 : const Elem & get_elem() const
927 5764112 : { libmesh_assert(this->_elem);
928 16344923 : return *(this->_elem); }
929 :
930 : /**
931 : * Accessor for current Elem object
932 : */
933 25910629 : Elem & get_elem()
934 25910629 : { libmesh_assert(this->_elem);
935 143016860 : return *(const_cast<Elem *>(this->_elem)); }
936 :
937 : /**
938 : * Accessor for current side of Elem object
939 : */
940 1144038 : unsigned char get_side() const
941 9306544 : { return side; }
942 :
943 : /**
944 : * Accessor for current edge of Elem object
945 : */
946 11524 : unsigned char get_edge() const
947 146872 : { return edge; }
948 :
949 : /**
950 : * Accessor for cached mesh dimension. This is the largest dimension
951 : * of the elements in the mesh. For the dimension of this->_elem, use
952 : * get_elem_dim();
953 : */
954 34346 : unsigned char get_dim() const
955 378126 : { return this->_dim; }
956 :
957 : /**
958 : * \returns The dimension of this->_elem. For mixed dimension meshes, this
959 : * may be different from get_dim(). If no element init has happened
960 : * yet, fall back on get_dim().
961 : */
962 1916460 : unsigned char get_elem_dim() const
963 619612553 : { return this->_elem ? this->_elem_dim : this->_dim; }
964 :
965 : /**
966 : * \returns Set of dimensions of elements present in the mesh at
967 : * context initialization.
968 : */
969 96662 : const std::set<unsigned char> & elem_dimensions() const
970 96662 : { return _elem_dims; }
971 :
972 : /**
973 : * Uses the coordinate data specified by mesh_*_position configuration
974 : * to set the geometry of \p elem to the value it would take after a fraction
975 : * \p theta of a timestep.
976 : */
977 : void elem_position_set(Real theta);
978 :
979 : /**
980 : * Uses the geometry of \p elem to set the coordinate data specified
981 : * by mesh_*_position configuration.
982 : */
983 : void elem_position_get();
984 :
985 : /**
986 : * Enum describing what data to use when initializing algebraic
987 : * structures on each element.
988 : */
989 : enum AlgebraicType { NONE = 0, // Do not reinitialize dof_indices
990 : DOFS_ONLY, // Reinitialize dof_indices, not
991 : // algebraic structures
992 : CURRENT, // Use dof_indices, current solution
993 : OLD, // Use old_dof_indices, custom solution
994 : OLD_DOFS_ONLY}; // Reinitialize old_dof_indices, not
995 : // algebraic structures
996 :
997 : /**
998 : * Setting which determines whether to initialize algebraic
999 : * structures (elem_*) on each element and set their values from
1000 : * current_local_solution. Algebraic initialization may be disabled
1001 : * for efficiency in cases where FEMContext is only used as a
1002 : * convenient container of FE objects.
1003 : */
1004 24740 : void set_algebraic_type(const AlgebraicType atype)
1005 330691 : { _atype = atype; }
1006 :
1007 : /*
1008 : * Get the current AlgebraicType setting
1009 : */
1010 418526547 : AlgebraicType algebraic_type() const { return _atype; }
1011 :
1012 : /**
1013 : * Set a NumericVector to be used in place of current_local_solution
1014 : * for calculating elem_solution. Set to nullptr to restore the
1015 : * current_local_solution behavior. Advanced DifferentiableSystem
1016 : * specific capabilities will only be enabled in the
1017 : * current_local_solution case.
1018 : */
1019 13858 : void set_custom_solution(const NumericVector<Number> * custom_sol)
1020 76065 : { _custom_solution = custom_sol; }
1021 :
1022 : /**
1023 : * Calls set_jacobian_tolerance() on all the FE objects controlled
1024 : * by this class. (Actually, it calls this on the underlying)
1025 : */
1026 : void set_jacobian_tolerance(Real tol);
1027 :
1028 : /**
1029 : * Current side for side_* to examine
1030 : */
1031 : unsigned char side;
1032 :
1033 : /**
1034 : * Current edge for edge_* to examine
1035 : */
1036 : unsigned char edge;
1037 :
1038 : /**
1039 : * Helper function for creating quadrature rules
1040 : */
1041 : FEType find_hardest_fe_type();
1042 :
1043 : /**
1044 : * Helper function for attaching quadrature rules
1045 : */
1046 : void attach_quadrature_rules();
1047 :
1048 : /**
1049 : * Return a pointer to the vector of active variables being computed
1050 : * for, or a null pointer if all variables in the system are active.
1051 : */
1052 63698 : const std::vector<unsigned int> * active_vars() const { return _active_vars.get(); }
1053 :
1054 : /**
1055 : * Helper function to reduce some code duplication in the *_point_* methods.
1056 : *
1057 : * get_derivative_level should be -1 to get_ everything, 0 to
1058 : * get_phi, 1 to get_dphi, 2 to get_d2phi, or 3 to get_curl_phi
1059 : */
1060 : template<typename OutputShape>
1061 : FEGenericBase<OutputShape> * build_new_fe(const FEGenericBase<OutputShape> * fe,
1062 : const Point & p,
1063 : const Real tolerance = TOLERANCE,
1064 : const int get_derivative_level = -1) const;
1065 :
1066 : protected:
1067 :
1068 : /**
1069 : * Variables on which to enable calculations, or nullptr if all
1070 : * variables in the System are to be enabled
1071 : */
1072 : std::unique_ptr<const std::vector<unsigned int>> _active_vars;
1073 :
1074 : /**
1075 : * System from which to acquire moving mesh information
1076 : */
1077 : System * _mesh_sys;
1078 :
1079 : /**
1080 : * Variables from which to acquire moving mesh information
1081 : */
1082 : unsigned int _mesh_x_var, _mesh_y_var, _mesh_z_var;
1083 :
1084 : /**
1085 : * Keep track of what type of algebra reinitialization is to be done
1086 : */
1087 : AlgebraicType _atype;
1088 :
1089 : /**
1090 : * Data with which to do algebra reinitialization
1091 : */
1092 : const NumericVector<Number> * _custom_solution;
1093 :
1094 : mutable std::unique_ptr<FEGenericBase<Real>> _real_fe;
1095 : mutable std::unique_ptr<FEGenericBase<RealGradient>> _real_grad_fe;
1096 : mutable int _real_fe_derivative_level;
1097 : mutable int _real_grad_fe_derivative_level;
1098 :
1099 : #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
1100 : mutable bool _real_fe_is_inf;
1101 : mutable bool _real_grad_fe_is_inf;
1102 : #endif
1103 :
1104 : template<typename OutputShape>
1105 : FEGenericBase<OutputShape> * cached_fe( const unsigned int elem_dim,
1106 : const FEType fe_type,
1107 : const int get_derivative_level ) const;
1108 :
1109 : /**
1110 : * Helper function to promote accessor usage
1111 : */
1112 : void set_elem( const Elem * e );
1113 :
1114 : // gcc-3.4, oracle 12.3 require this typedef to be public
1115 : // in order to use it in a return type
1116 : public:
1117 :
1118 : /**
1119 : * Helper typedef to simplify refactoring
1120 : */
1121 : typedef const DenseSubVector<Number> & (DiffContext::*diff_subsolution_getter)(unsigned int) const;
1122 :
1123 : protected:
1124 : /**
1125 : * Helper nested class for C++03-compatible "template typedef"
1126 : */
1127 : template <typename OutputType>
1128 : struct FENeeded
1129 : {
1130 : // Rank decrementer helper types
1131 : typedef typename TensorTools::DecrementRank<OutputType>::type Rank1Decrement;
1132 : typedef typename TensorTools::DecrementRank<Rank1Decrement>::type Rank2Decrement;
1133 :
1134 : // Typedefs for "Value getter" function pointer
1135 : typedef typename TensorTools::MakeReal<OutputType>::type value_shape;
1136 : typedef FEGenericBase<value_shape> value_base;
1137 : typedef void (FEMContext::*value_getter) (unsigned int, value_base *&, unsigned short) const;
1138 :
1139 : // Typedefs for "Grad getter" function pointer
1140 : typedef typename TensorTools::MakeReal<Rank1Decrement>::type grad_shape;
1141 : typedef FEGenericBase<grad_shape> grad_base;
1142 : typedef void (FEMContext::*grad_getter) (unsigned int, grad_base *&, unsigned short) const;
1143 :
1144 : // Typedefs for "Hessian getter" function pointer
1145 : typedef typename TensorTools::MakeReal<Rank2Decrement>::type hess_shape;
1146 : typedef FEGenericBase<hess_shape> hess_base;
1147 : typedef void (FEMContext::*hess_getter) (unsigned int, hess_base *&, unsigned short) const;
1148 : };
1149 :
1150 :
1151 :
1152 : /**
1153 : * Helper function to reduce some code duplication in the
1154 : * *interior_value methods.
1155 : */
1156 : template<typename OutputType,
1157 : typename FENeeded<OutputType>::value_getter fe_getter,
1158 : diff_subsolution_getter subsolution_getter>
1159 : void some_value(unsigned int var, unsigned int qp, OutputType & u) const;
1160 :
1161 : /**
1162 : * Helper function to reduce some code duplication in the
1163 : * *interior_gradient methods.
1164 : */
1165 : template<typename OutputType,
1166 : typename FENeeded<OutputType>::grad_getter fe_getter,
1167 : diff_subsolution_getter subsolution_getter>
1168 : void some_gradient(unsigned int var, unsigned int qp, OutputType & u) const;
1169 :
1170 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1171 : /**
1172 : * Helper function to reduce some code duplication in the
1173 : * *interior_hessian methods.
1174 : */
1175 : template<typename OutputType,
1176 : typename FENeeded<OutputType>::hess_getter fe_getter,
1177 : diff_subsolution_getter subsolution_getter>
1178 : void some_hessian(unsigned int var, unsigned int qp, OutputType & u) const;
1179 : #endif
1180 :
1181 : /**
1182 : * Finite element objects for each variable's interior, sides and edges.
1183 : * We store FE objects for each element dimension present in the mesh,
1184 : * except for edge_fe which only applies to 3D elements.
1185 : */
1186 : std::vector<std::map<FEType, std::unique_ptr<FEAbstract>>> _element_fe;
1187 : std::vector<std::map<FEType, std::unique_ptr<FEAbstract>>> _side_fe;
1188 : std::map<FEType, std::unique_ptr<FEAbstract>> _edge_fe;
1189 :
1190 :
1191 : /**
1192 : * Pointers to the same finite element objects, but indexed
1193 : * by variable number. We store FE objects for each element dimension
1194 : * present in the mesh, except for edge_fe_var which only applies
1195 : * for 3D elements.
1196 : */
1197 : std::vector<std::vector<FEAbstract *>> _element_fe_var;
1198 : std::vector<std::vector<FEAbstract *>> _side_fe_var;
1199 : std::vector<FEAbstract *> _edge_fe_var;
1200 :
1201 : /**
1202 : * Saved reference to BoundaryInfo on the mesh for this System.
1203 : * Used to answer boundary id requests.
1204 : */
1205 : const BoundaryInfo & _boundary_info;
1206 :
1207 : /**
1208 : * Current element for element_* to examine
1209 : */
1210 : const Elem * _elem;
1211 :
1212 : /**
1213 : * Cached dimension of largest dimension element in this mesh
1214 : */
1215 : unsigned char _dim;
1216 :
1217 : /**
1218 : * Cached dimension of this->_elem.
1219 : */
1220 : unsigned char _elem_dim;
1221 :
1222 : /**
1223 : * Cached dimensions of elements in the mesh, plus dimension 0 if
1224 : * SCALAR variables are in use.
1225 : */
1226 : std::set<unsigned char> _elem_dims;
1227 :
1228 : /**
1229 : * Quadrature rule for element interior.
1230 : * The FEM context will try to find a quadrature rule that
1231 : * correctly integrates all variables. We prepare quadrature
1232 : * rules for each element dimension in the mesh.
1233 : */
1234 : std::vector<std::unique_ptr<QBase>> _element_qrule;
1235 :
1236 : /**
1237 : * Quadrature rules for element sides
1238 : * The FEM context will try to find a quadrature rule that
1239 : * correctly integrates all variables. We prepare quadrature
1240 : * rules for each element dimension in the mesh.
1241 : */
1242 : std::vector<std::unique_ptr<QBase>> _side_qrule;
1243 :
1244 : /**
1245 : * Quadrature rules for element edges. If the FEM context is told
1246 : * to prepare for edge integration on 3D elements, it will try to
1247 : * find a quadrature rule that correctly integrates all variables.
1248 : * Because edge rules only apply to 3D elements, we don't need to
1249 : * worry about multiple dimensions
1250 : */
1251 : std::unique_ptr<QBase> _edge_qrule;
1252 :
1253 : /**
1254 : * The extra quadrature order for this context.
1255 : */
1256 : int _extra_quadrature_order;
1257 :
1258 : private:
1259 : /**
1260 : * Helper function used in constructors to set up internal data.
1261 : */
1262 : void init_internal_data(const System & sys);
1263 :
1264 : /**
1265 : * Uses the coordinate data specified by mesh_*_position configuration
1266 : * to set the geometry of \p elem to the value it would take after a fraction
1267 : * \p theta of a timestep.
1268 : *
1269 : * This does the work of elem_position_set, but isn't safe to call
1270 : * without _mesh_sys/etc. defined first.
1271 : */
1272 : void _do_elem_position_set(Real theta);
1273 :
1274 : /**
1275 : * Update the time in the context object for the given value of
1276 : * theta, based on the values of "time" and "deltat" stored in the
1277 : * system which created this context.
1278 : */
1279 : void _update_time_from_system(Real theta);
1280 : };
1281 :
1282 :
1283 :
1284 : // ------------------------------------------------------------
1285 : // FEMContext inline methods
1286 :
1287 : inline
1288 0 : void FEMContext::elem_position_set(Real theta)
1289 : {
1290 23040 : if (_mesh_sys)
1291 23040 : this->_do_elem_position_set(theta);
1292 0 : }
1293 :
1294 : template<typename OutputShape>
1295 : inline
1296 470455 : void FEMContext::get_element_fe( unsigned int var, FEGenericBase<OutputShape> *& fe,
1297 : unsigned short dim ) const
1298 : {
1299 470455 : libmesh_assert( !_element_fe_var[dim].empty() );
1300 470455 : libmesh_assert_less ( var, (_element_fe_var[dim].size() ) );
1301 654385763 : fe = cast_ptr<FEGenericBase<OutputShape> *>( (_element_fe_var[dim][var] ) );
1302 470455 : }
1303 :
1304 : inline
1305 89320 : void FEMContext::get_element_fe( unsigned int var, FEAbstract *& fe,
1306 : unsigned short dim ) const
1307 : {
1308 89320 : libmesh_assert( !_element_fe_var[dim].empty() );
1309 89320 : libmesh_assert_less ( var, (_element_fe_var[dim].size() ) );
1310 3236655 : fe = _element_fe_var[dim][var];
1311 89320 : }
1312 :
1313 : inline
1314 0 : FEBase * FEMContext::get_element_fe( unsigned int var, unsigned short dim ) const
1315 : {
1316 0 : libmesh_assert( !_element_fe_var[dim].empty() );
1317 0 : libmesh_assert_less ( var, (_element_fe_var[dim].size() ) );
1318 421230 : return cast_ptr<FEBase *>( (_element_fe_var[dim][var] ) );
1319 : }
1320 :
1321 : template<typename OutputShape>
1322 : inline
1323 1823638 : void FEMContext::get_side_fe( unsigned int var, FEGenericBase<OutputShape> *& fe,
1324 : unsigned short dim ) const
1325 : {
1326 1823638 : libmesh_assert( !_side_fe_var[dim].empty() );
1327 1823638 : libmesh_assert_less ( var, (_side_fe_var[dim].size() ) );
1328 19146836 : fe = cast_ptr<FEGenericBase<OutputShape> *>( (_side_fe_var[dim][var] ) );
1329 1823638 : }
1330 :
1331 : inline
1332 37638 : void FEMContext::get_side_fe( unsigned int var, FEAbstract *& fe,
1333 : unsigned short dim ) const
1334 : {
1335 37638 : libmesh_assert( !_side_fe_var[dim].empty() );
1336 37638 : libmesh_assert_less ( var, (_side_fe_var[dim].size() ) );
1337 1354276 : fe = _side_fe_var[dim][var];
1338 1241362 : }
1339 :
1340 : inline
1341 0 : FEBase * FEMContext::get_side_fe( unsigned int var, unsigned short dim ) const
1342 : {
1343 0 : libmesh_assert( !_side_fe_var[dim].empty() );
1344 0 : libmesh_assert_less ( var, (_side_fe_var[dim].size() ) );
1345 416 : return cast_ptr<FEBase *>( (_side_fe_var[dim][var] ) );
1346 : }
1347 :
1348 : template<typename OutputShape>
1349 : inline
1350 8748 : void FEMContext::get_edge_fe( unsigned int var, FEGenericBase<OutputShape> *& fe ) const
1351 : {
1352 8748 : libmesh_assert_less ( var, _edge_fe_var.size() );
1353 129722 : fe = cast_ptr<FEGenericBase<OutputShape> *>( _edge_fe_var[var] );
1354 8748 : }
1355 :
1356 : inline
1357 15464 : void FEMContext::get_edge_fe( unsigned int var, FEAbstract *& fe ) const
1358 : {
1359 15464 : libmesh_assert_less ( var, _edge_fe_var.size() );
1360 582831 : fe = _edge_fe_var[var];
1361 551903 : }
1362 :
1363 : inline
1364 : FEBase * FEMContext::get_edge_fe( unsigned int var ) const
1365 : {
1366 : libmesh_assert_less ( var, _edge_fe_var.size() );
1367 : return cast_ptr<FEBase *>( _edge_fe_var[var] );
1368 : }
1369 :
1370 :
1371 : } // namespace libMesh
1372 :
1373 : #endif // LIBMESH_FEM_CONTEXT_H
|
