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