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