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