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