libMesh
fem_context.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2017 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/boundary_info.h"
21 #include "libmesh/dof_map.h"
22 #include "libmesh/elem.h"
23 #include "libmesh/fe_base.h"
24 #include "libmesh/fe_interface.h"
25 #include "libmesh/fem_context.h"
26 #include "libmesh/libmesh_logging.h"
27 #include "libmesh/mesh_base.h"
28 #include "libmesh/quadrature.h"
29 #include "libmesh/system.h"
30 #include "libmesh/diff_system.h"
31 #include "libmesh/time_solver.h"
32 #include "libmesh/unsteady_solver.h" // For euler_residual
33 
34 namespace libMesh
35 {
36 
38  : DiffContext(sys),
39  _mesh_sys(libmesh_nullptr),
40  _mesh_x_var(0),
41  _mesh_y_var(0),
42  _mesh_z_var(0),
43  side(0), edge(0),
44  _atype(CURRENT),
45  _custom_solution(libmesh_nullptr),
46  _boundary_info(sys.get_mesh().get_boundary_info()),
47  _elem(libmesh_nullptr),
48  _dim(sys.get_mesh().mesh_dimension()),
49  _elem_dim(0), /* This will be reset in set_elem(). */
50  _elem_dims(sys.get_mesh().elem_dimensions()),
51  _element_qrule(4,libmesh_nullptr),
52  _side_qrule(4,libmesh_nullptr),
53  _extra_quadrature_order(sys.extra_quadrature_order)
54 {
55  init_internal_data(sys);
56 }
57 
58 FEMContext::FEMContext (const System & sys, int extra_quadrature_order)
59  : DiffContext(sys),
61  _mesh_x_var(0),
62  _mesh_y_var(0),
63  _mesh_z_var(0),
64  side(0), edge(0),
65  _atype(CURRENT),
67  _boundary_info(sys.get_mesh().get_boundary_info()),
69  _dim(sys.get_mesh().mesh_dimension()),
70  _elem_dim(0), /* This will be reset in set_elem(). */
71  _elem_dims(sys.get_mesh().elem_dimensions()),
74  _extra_quadrature_order(extra_quadrature_order)
75 {
76  init_internal_data(sys);
77 }
78 
80 {
81  // Reserve space for the FEAbstract and QBase objects for each
82  // element dimension possibility (0,1,2,3)
83  _element_fe.resize(4);
84  _side_fe.resize(4);
85  _element_fe_var.resize(4);
86  _side_fe_var.resize(4);
87 
88  // We need to know which of our variables has the hardest
89  // shape functions to numerically integrate.
90 
91  unsigned int nv = sys.n_vars();
92 
93  libmesh_assert (nv);
94  FEType hardest_fe_type = sys.variable_type(0);
95 
96  bool have_scalar = false;
97 
98  for (unsigned int i=0; i != nv; ++i)
99  {
100  FEType fe_type = sys.variable_type(i);
101 
102  // Make sure we find a non-SCALAR FE family, even in the case
103  // where the first variable(s) weren't
104  if (hardest_fe_type.family == SCALAR)
105  {
106  hardest_fe_type.family = fe_type.family;
107  hardest_fe_type.order = fe_type.order;
108  }
109 
110  // FIXME - we don't yet handle mixed finite elements from
111  // different families which require different quadrature rules
112  // libmesh_assert_equal_to (fe_type.family, hardest_fe_type.family);
113 
114  // We need to detect SCALAR's so we can prepare FE objects for
115  // them, and so we don't mistake high order scalars as a reason
116  // to crank up the quadrature order on other types.
117  if (fe_type.family == SCALAR)
118  have_scalar = true;
119  else if (fe_type.order > hardest_fe_type.order)
120  hardest_fe_type = fe_type;
121  }
122 
123  if (have_scalar)
124  // SCALAR FEs have dimension 0 by assumption
125  _elem_dims.insert(0);
126 
127  for (std::set<unsigned char>::const_iterator dim_it = _elem_dims.begin();
128  dim_it != _elem_dims.end(); ++dim_it)
129  {
130  const unsigned char dim = *dim_it;
131 
132  // Create an adequate quadrature rule
133  _element_qrule[dim] = hardest_fe_type.default_quadrature_rule
134  (dim, _extra_quadrature_order).release();
135  _side_qrule[dim] = hardest_fe_type.default_quadrature_rule
136  (dim-1, _extra_quadrature_order).release();
137  if (dim == 3)
138  _edge_qrule.reset(hardest_fe_type.default_quadrature_rule
139  (1, _extra_quadrature_order).release());
140 
141  // Next, create finite element objects
142  _element_fe_var[dim].resize(nv);
143  _side_fe_var[dim].resize(nv);
144  if (dim == 3)
145  _edge_fe_var.resize(nv);
146 
147 
148  for (unsigned int i=0; i != nv; ++i)
149  {
150  FEType fe_type = sys.variable_type(i);
151 
152  if (_element_fe[dim][fe_type] == libmesh_nullptr)
153  {
154  _element_fe[dim][fe_type] = FEAbstract::build(dim, fe_type).release();
155  _element_fe[dim][fe_type]->attach_quadrature_rule(_element_qrule[dim]);
156  _side_fe[dim][fe_type] = FEAbstract::build(dim, fe_type).release();
157  _side_fe[dim][fe_type]->attach_quadrature_rule(_side_qrule[dim]);
158 
159  if (dim == 3)
160  {
161  _edge_fe[fe_type] = FEAbstract::build(dim, fe_type).release();
162  _edge_fe[fe_type]->attach_quadrature_rule(_edge_qrule.get());
163  }
164  }
165 
166  _element_fe_var[dim][i] = _element_fe[dim][fe_type];
167  _side_fe_var[dim][i] = _side_fe[dim][fe_type];
168  if ((dim) == 3)
169  _edge_fe_var[i] = _edge_fe[fe_type];
170 
171  }
172  }
173 }
174 
176 {
177  // We don't want to store UniquePtrs in STL containers, but we don't
178  // want to leak memory either
179  for (std::vector<std::map<FEType, FEAbstract *>>::iterator d = _element_fe.begin();
180  d != _element_fe.end(); ++d)
181  for (std::map<FEType, FEAbstract *>::iterator i = d->begin();
182  i != d->end(); ++i)
183  delete i->second;
184 
185  for (std::vector<std::map<FEType, FEAbstract *>>::iterator d = _side_fe.begin();
186  d != _side_fe.end(); ++d)
187  for (std::map<FEType, FEAbstract *>::iterator i = d->begin();
188  i != d->end(); ++i)
189  delete i->second;
190 
191  for (std::map<FEType, FEAbstract *>::iterator i = _edge_fe.begin();
192  i != _edge_fe.end(); ++i)
193  delete i->second;
194  _edge_fe.clear();
195 
196  for (std::vector<QBase *>::iterator i = _element_qrule.begin();
197  i != _element_qrule.end(); ++i)
198  delete *i;
199  _element_qrule.clear();
200 
201  for (std::vector<QBase *>::iterator i = _side_qrule.begin();
202  i != _side_qrule.end(); ++i)
203  delete *i;
204  _side_qrule.clear();
205 }
206 
207 
208 
210 {
211  return _boundary_info.has_boundary_id(&(this->get_elem()), side, id);
212 }
213 
214 
215 #ifdef LIBMESH_ENABLE_DEPRECATED
216 std::vector<boundary_id_type> FEMContext::side_boundary_ids() const
217 {
218  libmesh_deprecated();
219  return _boundary_info.boundary_ids(&(this->get_elem()), side);
220 }
221 #endif
222 
223 
224 void FEMContext::side_boundary_ids(std::vector<boundary_id_type> & vec_to_fill) const
225 {
226  _boundary_info.boundary_ids(&(this->get_elem()), side, vec_to_fill);
227 }
228 
229 
230 
231 template<typename OutputType,
233  FEMContext::diff_subsolution_getter subsolution_getter>
234 void FEMContext::some_value(unsigned int var, unsigned int qp, OutputType & u) const
235 {
236  // Get local-to-global dof index lookup
237  const unsigned int n_dofs = cast_int<unsigned int>
238  (this->get_dof_indices(var).size());
239 
240  // Get current local coefficients
241  const DenseSubVector<Number> & coef = (this->*subsolution_getter)(var);
242 
243  // Get finite element object
245  (this->*fe_getter)( var, fe, this->get_elem_dim() );
246 
247  // Get shape function values at quadrature point
248  const std::vector<std::vector
249  <typename FENeeded<OutputType>::value_shape>> & phi = fe->get_phi();
250 
251  // Accumulate solution value
252  u = 0.;
253 
254  for (unsigned int l=0; l != n_dofs; l++)
255  u += phi[l][qp] * coef(l);
256 }
257 
258 
259 
260 template<typename OutputType,
262  FEMContext::diff_subsolution_getter subsolution_getter>
263 void FEMContext::some_gradient(unsigned int var, unsigned int qp, OutputType & du) const
264 {
265  // Get local-to-global dof index lookup
266  const unsigned int n_dofs = cast_int<unsigned int>
267  (this->get_dof_indices(var).size());
268 
269  // Get current local coefficients
270  const DenseSubVector<Number> & coef = (this->*subsolution_getter)(var);
271 
272  // Get finite element object
274  (this->*fe_getter)( var, fe, this->get_elem_dim() );
275 
276  // Get shape function values at quadrature point
277  const std::vector<std::vector
279  & dphi = fe->get_dphi();
280 
281  // Accumulate solution derivatives
282  du = 0;
283 
284  for (unsigned int l=0; l != n_dofs; l++)
285  du.add_scaled(dphi[l][qp], coef(l));
286 
287  return;
288 }
289 
290 
291 
292 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
293 template<typename OutputType,
295  FEMContext::diff_subsolution_getter subsolution_getter>
296 void FEMContext::some_hessian(unsigned int var, unsigned int qp, OutputType & d2u) const
297 {
298  // Get local-to-global dof index lookup
299  const unsigned int n_dofs = cast_int<unsigned int>
300  (this->get_dof_indices(var).size());
301 
302  // Get current local coefficients
303  const DenseSubVector<Number> & coef = (this->*subsolution_getter)(var);
304 
305  // Get finite element object
307  (this->*fe_getter)( var, fe, this->get_elem_dim() );
308 
309  // Get shape function values at quadrature point
310  const std::vector<std::vector
312  & d2phi = fe->get_d2phi();
313 
314  // Accumulate solution second derivatives
315  d2u = 0.0;
316 
317  for (unsigned int l=0; l != n_dofs; l++)
318  d2u.add_scaled(d2phi[l][qp], coef(l));
319 
320  return;
321 }
322 #endif
323 
324 
325 
326 Number FEMContext::interior_value(unsigned int var, unsigned int qp) const
327 {
328  Number u;
329 
330  this->interior_value( var, qp, u );
331 
332  return u;
333 }
334 
335 template<typename OutputType>
336 void FEMContext::interior_value(unsigned int var, unsigned int qp,
337  OutputType & u) const
338 {
339  this->some_value<OutputType,
340  &FEMContext::get_element_fe<typename TensorTools::MakeReal<OutputType>::type>,
341  &DiffContext::get_elem_solution>(var, qp, u);
342 }
343 
344 
345 template<typename OutputType>
346 void FEMContext::interior_values (unsigned int var,
347  const NumericVector<Number> & _system_vector,
348  std::vector<OutputType> & u_vals) const
349 {
350  typedef typename TensorTools::MakeReal<OutputType>::type OutputShape;
351 
352  // Get local-to-global dof index lookup
353  const unsigned int n_dofs = cast_int<unsigned int>
354  (this->get_dof_indices(var).size());
355 
356  // Get current local coefficients
357  const DenseSubVector<Number> & coef = get_localized_subvector(_system_vector, var);
358 
359  // Get the finite element object
361  this->get_element_fe<OutputShape>( var, fe, this->get_elem_dim() );
362 
363  // Get shape function values at quadrature point
364  const std::vector<std::vector<OutputShape>> & phi = fe->get_phi();
365 
366  // Loop over all the q_points on this element
367  for (std::size_t qp=0; qp != u_vals.size(); qp++)
368  {
369  OutputType & u = u_vals[qp];
370 
371  // Compute the value at this q_point
372  u = 0.;
373 
374  for (unsigned int l=0; l != n_dofs; l++)
375  u += phi[l][qp] * coef(l);
376  }
377 
378  return;
379 }
380 
382  unsigned int qp) const
383 {
384  Gradient du;
385 
386  this->interior_gradient( var, qp, du );
387 
388  return du;
389 }
390 
391 
392 
393 template<typename OutputType>
394 void FEMContext::interior_gradient(unsigned int var,
395  unsigned int qp,
396  OutputType & du) const
397 {
398  this->some_gradient<OutputType,
401  <OutputType>::type>::type>,
402  &DiffContext::get_elem_solution>(var, qp, du);
403 }
404 
405 
406 
407 template<typename OutputType>
408 void FEMContext::interior_gradients(unsigned int var,
409  const NumericVector<Number> & _system_vector,
410  std::vector<OutputType> & du_vals) const
411 {
412  typedef typename TensorTools::MakeReal
414  OutputShape;
415 
416  // Get local-to-global dof index lookup
417  const unsigned int n_dofs = cast_int<unsigned int>
418  (this->get_dof_indices(var).size());
419 
420  // Get current local coefficients
421  const DenseSubVector<Number> & coef = get_localized_subvector(_system_vector, var);
422 
423  // Get finite element object
425  this->get_element_fe<OutputShape>( var, fe, this->get_elem_dim() );
426 
427  // Get shape function values at quadrature point
428  const std::vector<std::vector<typename FEGenericBase<OutputShape>::OutputGradient>> & dphi = fe->get_dphi();
429 
430  // Loop over all the q_points in this finite element
431  for (std::size_t qp=0; qp != du_vals.size(); qp++)
432  {
433  OutputType & du = du_vals[qp];
434 
435  // Compute the gradient at this q_point
436  du = 0;
437 
438  for (unsigned int l=0; l != n_dofs; l++)
439  du.add_scaled(dphi[l][qp], coef(l));
440  }
441 
442  return;
443 }
444 
445 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
446 Tensor FEMContext::interior_hessian(unsigned int var, unsigned int qp) const
447 {
448  Tensor d2u;
449 
450  this->interior_hessian( var, qp, d2u );
451 
452  return d2u;
453 }
454 
455 template<typename OutputType>
456 void FEMContext::interior_hessian(unsigned int var, unsigned int qp,
457  OutputType & d2u) const
458 {
459  this->some_hessian<OutputType,
461  <typename TensorTools::MakeReal
464  <OutputType>::type>::type>::type>,
465  &DiffContext::get_elem_solution>(var, qp, d2u);
466 }
467 
468 
469 template<typename OutputType>
470 void FEMContext::interior_hessians(unsigned int var,
471  const NumericVector<Number> & _system_vector,
472  std::vector<OutputType> & d2u_vals) const
473 {
474  typedef typename TensorTools::DecrementRank<OutputType>::type Rank1Decrement;
475  typedef typename TensorTools::DecrementRank<Rank1Decrement>::type Rank2Decrement;
476  typedef typename TensorTools::MakeReal<Rank2Decrement>::type OutputShape;
477 
478  // Get local-to-global dof index lookup
479  const unsigned int n_dofs = cast_int<unsigned int>
480  (this->get_dof_indices(var).size());
481 
482  // Get current local coefficients
483  const DenseSubVector<Number> & coef = get_localized_subvector(_system_vector, var);
484 
485  // Get finite element object
487  this->get_element_fe<OutputShape>( var, fe, this->get_elem_dim() );
488 
489  // Get shape function values at quadrature point
490  const std::vector<std::vector<typename FEGenericBase<OutputShape>::OutputTensor>> & d2phi = fe->get_d2phi();
491 
492  // Loop over all the q_points in this finite element
493  for (std::size_t qp=0; qp != d2u_vals.size(); qp++)
494  {
495  OutputType & d2u = d2u_vals[qp];
496 
497  // Compute the gradient at this q_point
498  d2u = 0;
499 
500  for (unsigned int l=0; l != n_dofs; l++)
501  d2u.add_scaled(d2phi[l][qp], coef(l));
502  }
503 
504  return;
505 }
506 
507 
508 #endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
509 
510 
511 template<typename OutputType>
512 void FEMContext::interior_curl(unsigned int var, unsigned int qp,
513  OutputType & curl_u) const
514 {
515  typedef typename TensorTools::MakeReal<OutputType>::type OutputShape;
516 
517  // Get local-to-global dof index lookup
518  const unsigned int n_dofs = cast_int<unsigned int>
519  (this->get_dof_indices(var).size());
520 
521  // Get current local coefficients
522  libmesh_assert_greater (this->_elem_subsolutions.size(), var);
523  const DenseSubVector<Number> & coef = this->get_elem_solution(var);
524 
525  // Get finite element object
527  this->get_element_fe<OutputShape>( var, fe, this->get_elem_dim() );
528 
529  // Get shape function values at quadrature point
530  const std::vector<std::vector<typename FEGenericBase<OutputShape>::OutputShape>> & curl_phi = fe->get_curl_phi();
531 
532  // Accumulate solution curl
533  curl_u = 0.;
534 
535  for (unsigned int l=0; l != n_dofs; l++)
536  curl_u.add_scaled(curl_phi[l][qp], coef(l));
537 
538  return;
539 }
540 
541 
542 template<typename OutputType>
543 void FEMContext::interior_div(unsigned int var, unsigned int qp,
544  OutputType & div_u) const
545 {
546  typedef typename
548  <typename TensorTools::MakeReal<OutputType>::type>::type OutputShape;
549 
550  // Get local-to-global dof index lookup
551  const unsigned int n_dofs = cast_int<unsigned int>
552  (this->get_dof_indices(var).size());
553 
554  // Get current local coefficients
555  libmesh_assert_greater (this->_elem_subsolutions.size(), var);
556  const DenseSubVector<Number> & coef = this->get_elem_solution(var);
557 
558  // Get finite element object
560  this->get_element_fe<OutputShape>( var, fe, this->get_elem_dim() );
561 
562  // Get shape function values at quadrature point
563  const std::vector<std::vector<typename FEGenericBase<OutputShape>::OutputDivergence>> & div_phi = fe->get_div_phi();
564 
565  // Accumulate solution curl
566  div_u = 0.;
567 
568  for (unsigned int l=0; l != n_dofs; l++)
569  div_u += div_phi[l][qp] * coef(l);
570 
571  return;
572 }
573 
574 
576  unsigned int qp) const
577 {
578  Number u = 0.;
579 
580  this->side_value( var, qp, u );
581 
582  return u;
583 }
584 
585 
586 template<typename OutputType>
587 void FEMContext::side_value(unsigned int var,
588  unsigned int qp,
589  OutputType & u) const
590 {
591  this->some_value<OutputType,
592  &FEMContext::get_side_fe<typename TensorTools::MakeReal<OutputType>::type>,
593  &DiffContext::get_elem_solution>(var, qp, u);
594 }
595 
596 
597 template<typename OutputType>
598 void FEMContext::side_values(unsigned int var,
599  const NumericVector<Number> & _system_vector,
600  std::vector<OutputType> & u_vals) const
601 {
602  typedef typename TensorTools::MakeReal<OutputType>::type OutputShape;
603 
604  // Get local-to-global dof index lookup
605  const unsigned int n_dofs = cast_int<unsigned int>
606  (this->get_dof_indices(var).size());
607 
608  // Get current local coefficients
609  const DenseSubVector<Number> & coef = get_localized_subvector(_system_vector, var);
610 
611  // Get the finite element object
613  this->get_side_fe<OutputShape>( var, the_side_fe, this->get_elem_dim() );
614 
615  // Get shape function values at quadrature point
616  const std::vector<std::vector<OutputShape>> & phi = the_side_fe->get_phi();
617 
618  // Loop over all the q_points on this element
619  for (std::size_t qp=0; qp != u_vals.size(); qp++)
620  {
621  OutputType & u = u_vals[qp];
622 
623  // Compute the value at this q_point
624  u = 0.;
625 
626  for (unsigned int l=0; l != n_dofs; l++)
627  u += phi[l][qp] * coef(l);
628  }
629 
630  return;
631 }
632 
633 Gradient FEMContext::side_gradient(unsigned int var, unsigned int qp) const
634 {
635  Gradient du;
636 
637  this->side_gradient( var, qp, du );
638 
639  return du;
640 }
641 
642 
643 template<typename OutputType>
644 void FEMContext::side_gradient(unsigned int var, unsigned int qp,
645  OutputType & du) const
646 {
647  typedef typename TensorTools::MakeReal
649  OutputShape;
650 
651  // Get local-to-global dof index lookup
652  const unsigned int n_dofs = cast_int<unsigned int>
653  (this->get_dof_indices(var).size());
654 
655  // Get current local coefficients
656  libmesh_assert_greater (this->_elem_subsolutions.size(), var);
657  const DenseSubVector<Number> & coef = this->get_elem_solution(var);
658 
659  // Get finite element object
661  this->get_side_fe<OutputShape>( var, the_side_fe, this->get_elem_dim() );
662 
663  // Get shape function values at quadrature point
664  const std::vector<std::vector<typename FEGenericBase<OutputShape>::OutputGradient>> & dphi = the_side_fe->get_dphi();
665 
666  // Accumulate solution derivatives
667  du = 0.;
668 
669  for (unsigned int l=0; l != n_dofs; l++)
670  du.add_scaled(dphi[l][qp], coef(l));
671 
672  return;
673 }
674 
675 
676 
677 template<typename OutputType>
678 void FEMContext::side_gradients(unsigned int var,
679  const NumericVector<Number> & _system_vector,
680  std::vector<OutputType> & du_vals) const
681 {
682  typedef typename TensorTools::MakeReal
684  OutputShape;
685 
686  // Get local-to-global dof index lookup
687  const unsigned int n_dofs = cast_int<unsigned int>
688  (this->get_dof_indices(var).size());
689 
690  // Get current local coefficients
691  const DenseSubVector<Number> & coef = get_localized_subvector(_system_vector, var);
692 
693  // Get finite element object
695  this->get_side_fe<OutputShape>( var, the_side_fe, this->get_elem_dim() );
696 
697  // Get shape function values at quadrature point
698  const std::vector<std::vector<typename FEGenericBase<OutputShape>::OutputGradient>> & dphi = the_side_fe->get_dphi();
699 
700  // Loop over all the q_points in this finite element
701  for (std::size_t qp=0; qp != du_vals.size(); qp++)
702  {
703  OutputType & du = du_vals[qp];
704 
705  du = 0;
706 
707  // Compute the gradient at this q_point
708  for (unsigned int l=0; l != n_dofs; l++)
709  du.add_scaled(dphi[l][qp], coef(l));
710  }
711 
712  return;
713 }
714 
715 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
717  unsigned int qp) const
718 {
719  Tensor d2u;
720 
721  this->side_hessian( var, qp, d2u );
722 
723  return d2u;
724 }
725 
726 
727 
728 template<typename OutputType>
729 void FEMContext::side_hessian(unsigned int var,
730  unsigned int qp,
731  OutputType & d2u) const
732 {
733  this->some_hessian<OutputType,
735  <typename TensorTools::MakeReal
738  <OutputType>::type>::type>::type>,
739  &DiffContext::get_elem_solution>(var, qp, d2u);
740 }
741 
742 
743 
744 template<typename OutputType>
745 void FEMContext::side_hessians(unsigned int var,
746  const NumericVector<Number> & _system_vector,
747  std::vector<OutputType> & d2u_vals) const
748 {
749  typedef typename TensorTools::DecrementRank<OutputType>::type Rank1Decrement;
750  typedef typename TensorTools::DecrementRank<Rank1Decrement>::type Rank2Decrement;
751  typedef typename TensorTools::MakeReal<Rank2Decrement>::type OutputShape;
752 
753  // Get local-to-global dof index lookup
754  const unsigned int n_dofs = cast_int<unsigned int>
755  (this->get_dof_indices(var).size());
756 
757  // Get current local coefficients
758  const DenseSubVector<Number> & coef = get_localized_subvector(_system_vector, var);
759 
760  // Get finite element object
762  this->get_side_fe<OutputShape>( var, the_side_fe, this->get_elem_dim() );
763 
764  // Get shape function values at quadrature point
765  const std::vector<std::vector<typename FEGenericBase<OutputShape>::OutputTensor>> & d2phi = the_side_fe->get_d2phi();
766 
767  // Loop over all the q_points in this finite element
768  for (std::size_t qp=0; qp != d2u_vals.size(); qp++)
769  {
770  OutputType & d2u = d2u_vals[qp];
771 
772  // Compute the gradient at this q_point
773  d2u = 0;
774 
775  for (unsigned int l=0; l != n_dofs; l++)
776  d2u.add_scaled(d2phi[l][qp], coef(l));
777  }
778 
779  return;
780 }
781 
782 
783 
784 #endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
785 
786 
787 
788 Number FEMContext::point_value(unsigned int var, const Point & p) const
789 {
790  Number u = 0.;
791 
792  this->point_value( var, p, u );
793 
794  return u;
795 }
796 
797 template<typename OutputType>
798 void FEMContext::point_value(unsigned int var,
799  const Point & p,
800  OutputType & u,
801  const Real tolerance) const
802 {
803  typedef typename TensorTools::MakeReal<OutputType>::type OutputShape;
804 
805  // Get local-to-global dof index lookup
806  const unsigned int n_dofs = cast_int<unsigned int>
807  (this->get_dof_indices(var).size());
808 
809  // Get current local coefficients
810  libmesh_assert_greater (this->_elem_subsolutions.size(), var);
811  const DenseSubVector<Number> & coef = this->get_elem_solution(var);
812 
813  // Get finite element object
815  this->get_element_fe<OutputShape>( var, fe, this->get_elem_dim() );
816 
817  // Build a FE for calculating u(p)
818  FEGenericBase<OutputShape> * fe_new =
819  this->build_new_fe( fe, p, tolerance );
820 
821  // Get the values of the shape function derivatives
822  const std::vector<std::vector<OutputShape>> & phi = fe_new->get_phi();
823 
824  u = 0.;
825 
826  for (unsigned int l=0; l != n_dofs; l++)
827  u += phi[l][0] * coef(l);
828 
829  return;
830 }
831 
832 
833 
834 Gradient FEMContext::point_gradient(unsigned int var, const Point & p) const
835 {
836  Gradient grad_u;
837 
838  this->point_gradient( var, p, grad_u );
839 
840  return grad_u;
841 }
842 
843 
844 
845 template<typename OutputType>
846 void FEMContext::point_gradient(unsigned int var,
847  const Point & p,
848  OutputType & grad_u,
849  const Real tolerance) const
850 {
851  typedef typename TensorTools::MakeReal
853  OutputShape;
854 
855  // Get local-to-global dof index lookup
856  const unsigned int n_dofs = cast_int<unsigned int>
857  (this->get_dof_indices(var).size());
858 
859  // Get current local coefficients
860  libmesh_assert_greater (this->_elem_subsolutions.size(), var);
861  const DenseSubVector<Number> & coef = this->get_elem_solution(var);
862 
863  // Get finite element object
865  this->get_element_fe<OutputShape>( var, fe, this->get_elem_dim() );
866 
867  // Build a FE for calculating u(p)
868  FEGenericBase<OutputShape> * fe_new =
869  this->build_new_fe( fe, p, tolerance );
870 
871  // Get the values of the shape function derivatives
872  const std::vector<std::vector<typename FEGenericBase<OutputShape>::OutputGradient>> & dphi = fe_new->get_dphi();
873 
874  grad_u = 0.0;
875 
876  for (unsigned int l=0; l != n_dofs; l++)
877  grad_u.add_scaled(dphi[l][0], coef(l));
878 
879  return;
880 }
881 
882 
883 
884 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
885 
886 Tensor FEMContext::point_hessian(unsigned int var, const Point & p) const
887 {
888  Tensor hess_u;
889 
890  this->point_hessian( var, p, hess_u );
891 
892  return hess_u;
893 }
894 
895 
896 template<typename OutputType>
897 void FEMContext::point_hessian(unsigned int var,
898  const Point & p,
899  OutputType & hess_u,
900  const Real tolerance) const
901 {
902  typedef typename TensorTools::DecrementRank<OutputType>::type Rank1Decrement;
903  typedef typename TensorTools::DecrementRank<Rank1Decrement>::type Rank2Decrement;
904  typedef typename TensorTools::MakeReal<Rank2Decrement>::type OutputShape;
905 
906  // Get local-to-global dof index lookup
907  const unsigned int n_dofs = cast_int<unsigned int>
908  (this->get_dof_indices(var).size());
909 
910  // Get current local coefficients
911  libmesh_assert_greater (this->_elem_subsolutions.size(), var);
912  const DenseSubVector<Number> & coef = this->get_elem_solution(var);
913 
914  // Get finite element object
916  this->get_element_fe<OutputShape>( var, fe, this->get_elem_dim() );
917 
918  // Build a FE for calculating u(p)
919  FEGenericBase<OutputShape> * fe_new =
920  this->build_new_fe( fe, p, tolerance );
921 
922  // Get the values of the shape function derivatives
923  const std::vector<std::vector<typename FEGenericBase<OutputShape>::OutputTensor>> & d2phi = fe_new->get_d2phi();
924 
925  hess_u = 0.0;
926 
927  for (unsigned int l=0; l != n_dofs; l++)
928  hess_u.add_scaled(d2phi[l][0], coef(l));
929 
930  return;
931 }
932 
933 #endif // LIBMESH_ENABLE_SECOND_DERIVATIVES
934 
935 
936 template<typename OutputType>
937 void FEMContext::point_curl(unsigned int var,
938  const Point & p,
939  OutputType & curl_u,
940  const Real tolerance) const
941 {
942  typedef typename TensorTools::MakeReal<OutputType>::type OutputShape;
943 
944  // Get local-to-global dof index lookup
945  const unsigned int n_dofs = cast_int<unsigned int>
946  (this->get_dof_indices(var).size());
947 
948  // Get current local coefficients
949  libmesh_assert_greater (this->_elem_subsolutions.size(), var);
950  const DenseSubVector<Number> & coef = this->get_elem_solution(var);
951 
952  // Get finite element object
954  this->get_element_fe<OutputShape>( var, fe, this->get_elem_dim() );
955 
956  // Build a FE for calculating u(p)
957  FEGenericBase<OutputShape> * fe_new =
958  this->build_new_fe( fe, p, tolerance );
959 
960  // Get the values of the shape function derivatives
961  const std::vector<std::vector<typename FEGenericBase<OutputShape>::OutputShape>> & curl_phi = fe_new->get_curl_phi();
962 
963  curl_u = 0.0;
964 
965  for (unsigned int l=0; l != n_dofs; l++)
966  curl_u.add_scaled(curl_phi[l][0], coef(l));
967 
968  return;
969 }
970 
971 
972 
973 Number FEMContext::fixed_interior_value(unsigned int var, unsigned int qp) const
974 {
975  Number u = 0.;
976 
977  this->fixed_interior_value( var, qp, u );
978 
979  return u;
980 }
981 
982 
983 
984 template<typename OutputType>
985 void FEMContext::fixed_interior_value(unsigned int var, unsigned int qp,
986  OutputType & u) const
987 {
988  this->some_value<OutputType,
992 }
993 
994 
995 
996 Gradient FEMContext::fixed_interior_gradient(unsigned int var, unsigned int qp) const
997 {
998  Gradient du;
999 
1000  this->fixed_interior_gradient( var, qp, du );
1001 
1002  return du;
1003 }
1004 
1005 
1006 template<typename OutputType>
1007 void FEMContext::fixed_interior_gradient(unsigned int var, unsigned int qp,
1008  OutputType & du) const
1009 {
1010  this->some_gradient
1011  <OutputType,
1013  <typename TensorTools::MakeReal
1014  <typename TensorTools::DecrementRank
1015  <OutputType>::type>::type>,
1017  (var, qp, du);
1018 }
1019 
1020 
1021 
1022 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1023 Tensor FEMContext::fixed_interior_hessian(unsigned int var, unsigned int qp) const
1024 {
1025  Tensor d2u;
1026 
1027  this->fixed_interior_hessian( var, qp, d2u );
1028 
1029  return d2u;
1030 }
1031 
1032 
1033 template<typename OutputType>
1034 void FEMContext::fixed_interior_hessian(unsigned int var, unsigned int qp,
1035  OutputType & d2u) const
1036 {
1037  this->some_hessian<OutputType,
1039  <typename TensorTools::MakeReal
1040  <typename TensorTools::DecrementRank
1041  <typename TensorTools::DecrementRank
1042  <OutputType>::type>::type>::type>,
1043  &DiffContext::get_elem_fixed_solution>(var, qp, d2u);
1044 }
1045 #endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1046 
1047 
1048 
1049 Number FEMContext::fixed_side_value(unsigned int var, unsigned int qp) const
1050 {
1051  Number u = 0.;
1052 
1053  this->fixed_side_value( var, qp, u );
1054 
1055  return u;
1056 }
1057 
1058 
1059 template<typename OutputType>
1060 void FEMContext::fixed_side_value(unsigned int var, unsigned int qp,
1061  OutputType & u) const
1062 {
1063  this->some_value
1064  <OutputType,
1068  (var, qp, u);
1069 }
1070 
1071 
1072 
1073 Gradient FEMContext::fixed_side_gradient(unsigned int var, unsigned int qp) const
1074 {
1075  Gradient du;
1076 
1077  this->fixed_side_gradient( var, qp, du );
1078 
1079  return du;
1080 }
1081 
1082 
1083 template<typename OutputType>
1084 void FEMContext::fixed_side_gradient(unsigned int var, unsigned int qp,
1085  OutputType & du) const
1086 {
1087  this->some_gradient<OutputType,
1089  <typename TensorTools::MakeReal
1090  <typename TensorTools::DecrementRank
1091  <OutputType>::type>::type>,
1092  &DiffContext::get_elem_fixed_solution>(var, qp, du);
1093 }
1094 
1095 
1096 
1097 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1098 Tensor FEMContext::fixed_side_hessian(unsigned int var, unsigned int qp) const
1099 {
1100  Tensor d2u;
1101 
1102  this->fixed_side_hessian( var, qp, d2u );
1103 
1104  return d2u;
1105 }
1106 
1107 template<typename OutputType>
1108 void FEMContext::fixed_side_hessian(unsigned int var, unsigned int qp,
1109  OutputType & d2u) const
1110 {
1111  this->some_hessian<OutputType,
1113  <typename TensorTools::MakeReal
1114  <typename TensorTools::DecrementRank
1115  <typename TensorTools::DecrementRank
1116  <OutputType>::type>::type>::type>,
1117  &DiffContext::get_elem_fixed_solution>(var, qp, d2u);
1118 }
1119 #endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1120 
1121 
1122 
1123 Number FEMContext::fixed_point_value(unsigned int var, const Point & p) const
1124 {
1125  Number u = 0.;
1126 
1127  this->fixed_point_value( var, p, u );
1128 
1129  return u;
1130 }
1131 
1132 template<typename OutputType>
1133 void FEMContext::fixed_point_value(unsigned int var,
1134  const Point & p,
1135  OutputType & u,
1136  const Real tolerance) const
1137 {
1138  typedef typename TensorTools::MakeReal<OutputType>::type OutputShape;
1139 
1140  // Get local-to-global dof index lookup
1141  const unsigned int n_dofs = cast_int<unsigned int>
1142  (this->get_dof_indices(var).size());
1143 
1144  // Get current local coefficients
1145  libmesh_assert_greater (_elem_fixed_subsolutions.size(), var);
1146  const DenseSubVector<Number> & coef = this->get_elem_fixed_solution(var);
1147 
1148  // Get finite element object
1150  this->get_element_fe<OutputShape>( var, fe, this->get_elem_dim() );
1151 
1152  // Build a FE for calculating u(p)
1153  FEGenericBase<OutputShape> * fe_new =
1154  this->build_new_fe( fe, p, tolerance );
1155 
1156  // Get the values of the shape function derivatives
1157  const std::vector<std::vector<OutputShape>> & phi = fe_new->get_phi();
1158 
1159  u = 0.;
1160 
1161  for (unsigned int l=0; l != n_dofs; l++)
1162  u += phi[l][0] * coef(l);
1163 
1164  return;
1165 }
1166 
1167 
1168 
1169 Gradient FEMContext::fixed_point_gradient(unsigned int var, const Point & p) const
1170 {
1171  Gradient grad_u;
1172 
1173  this->fixed_point_gradient( var, p, grad_u );
1174 
1175  return grad_u;
1176 }
1177 
1178 
1179 
1180 template<typename OutputType>
1181 void FEMContext::fixed_point_gradient(unsigned int var,
1182  const Point & p,
1183  OutputType & grad_u,
1184  const Real tolerance) const
1185 {
1186  typedef typename TensorTools::MakeReal
1188  OutputShape;
1189 
1190  // Get local-to-global dof index lookup
1191  const unsigned int n_dofs = cast_int<unsigned int>
1192  (this->get_dof_indices(var).size());
1193 
1194  // Get current local coefficients
1195  libmesh_assert_greater (_elem_fixed_subsolutions.size(), var);
1196  const DenseSubVector<Number> & coef = this->get_elem_fixed_solution(var);
1197 
1198  // Get finite element object
1200  this->get_element_fe<OutputShape>( var, fe, this->get_elem_dim() );
1201 
1202  // Build a FE for calculating u(p)
1203  FEGenericBase<OutputShape> * fe_new =
1204  this->build_new_fe( fe, p, tolerance );
1205 
1206  // Get the values of the shape function derivatives
1207  const std::vector<std::vector<typename FEGenericBase<OutputShape>::OutputGradient>> & dphi = fe_new->get_dphi();
1208 
1209  grad_u = 0.0;
1210 
1211  for (unsigned int l=0; l != n_dofs; l++)
1212  grad_u.add_scaled(dphi[l][0], coef(l));
1213 
1214  return;
1215 }
1216 
1217 
1218 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1219 
1220 Tensor FEMContext::fixed_point_hessian(unsigned int var, const Point & p) const
1221 {
1222  Tensor hess_u;
1223 
1224  this->fixed_point_hessian( var, p, hess_u );
1225 
1226  return hess_u;
1227 }
1228 
1229 
1230 
1231 template<typename OutputType>
1232 void FEMContext::fixed_point_hessian(unsigned int var,
1233  const Point & p,
1234  OutputType & hess_u,
1235  const Real tolerance) const
1236 {
1237  typedef typename TensorTools::DecrementRank<OutputType>::type Rank1Decrement;
1238  typedef typename TensorTools::DecrementRank<Rank1Decrement>::type Rank2Decrement;
1239  typedef typename TensorTools::MakeReal<Rank2Decrement>::type OutputShape;
1240 
1241  // Get local-to-global dof index lookup
1242  const unsigned int n_dofs = cast_int<unsigned int>
1243  (this->get_dof_indices(var).size());
1244 
1245  // Get current local coefficients
1246  libmesh_assert_greater (_elem_fixed_subsolutions.size(), var);
1247  const DenseSubVector<Number> & coef = this->get_elem_fixed_solution(var);
1248 
1249  // Get finite element object
1251  this->get_element_fe<OutputShape>( var, fe, this->get_elem_dim() );
1252 
1253  // Build a FE for calculating u(p)
1254  FEGenericBase<OutputShape> * fe_new =
1255  this->build_new_fe( fe, p, tolerance );
1256 
1257  // Get the values of the shape function derivatives
1258  const std::vector<std::vector<typename FEGenericBase<OutputShape>::OutputTensor>> & d2phi = fe_new->get_d2phi();
1259 
1260  hess_u = 0.0;
1261 
1262  for (unsigned int l=0; l != n_dofs; l++)
1263  hess_u.add_scaled(d2phi[l][0], coef(l));
1264 
1265  return;
1266 }
1267 
1268 #endif // LIBMESH_ENABLE_SECOND_DERIVATIVES
1269 
1270 
1271 
1272 template<typename OutputType>
1273 void FEMContext::interior_rate(unsigned int var, unsigned int qp,
1274  OutputType & u) const
1275 {
1276  this->some_value<OutputType,
1280 }
1281 
1282 
1283 
1284 template<typename OutputType>
1285 void FEMContext::side_rate(unsigned int var, unsigned int qp,
1286  OutputType & u) const
1287 {
1288  this->some_value<OutputType,
1292 }
1293 
1294 template<typename OutputType>
1295 void FEMContext::interior_accel(unsigned int var, unsigned int qp,
1296  OutputType & u) const
1297 {
1298  this->some_value<OutputType,
1302 }
1303 
1304 
1305 
1306 template<typename OutputType>
1307 void FEMContext::side_accel(unsigned int var, unsigned int qp,
1308  OutputType & u) const
1309 {
1310  this->some_value<OutputType,
1314 }
1315 
1316 
1317 
1319 {
1320  // Update the "time" variable of this context object
1321  this->_update_time_from_system(theta);
1322 
1323  // Handle a moving element if necessary.
1324  if (_mesh_sys)
1325  {
1326  // We assume that the ``default'' state
1327  // of the mesh is its final, theta=1.0
1328  // position, so we don't bother with
1329  // mesh motion in that case.
1330  if (theta != 1.0)
1331  {
1332  // FIXME - ALE is not threadsafe yet!
1333  libmesh_assert_equal_to (libMesh::n_threads(), 1);
1334 
1335  elem_position_set(theta);
1336  }
1337  elem_fe_reinit();
1338  }
1339 }
1340 
1341 
1343 {
1344  // Update the "time" variable of this context object
1345  this->_update_time_from_system(theta);
1346 
1347  // Handle a moving element if necessary
1348  if (_mesh_sys)
1349  {
1350  // FIXME - not threadsafe yet!
1351  elem_position_set(theta);
1352  side_fe_reinit();
1353  }
1354 }
1355 
1356 
1358 {
1359  // Update the "time" variable of this context object
1360  this->_update_time_from_system(theta);
1361 
1362  // Handle a moving element if necessary
1363  if (_mesh_sys)
1364  {
1365  // FIXME - not threadsafe yet!
1366  elem_position_set(theta);
1367  edge_fe_reinit();
1368  }
1369 }
1370 
1371 
1373 {
1374  // Update the "time" variable of this context object
1375  this->_update_time_from_system(theta);
1376 
1377  // We can reuse the Elem FE safely here.
1378  elem_fe_reinit();
1379 }
1380 
1381 
1382 void FEMContext::elem_fe_reinit(const std::vector<Point> * const pts)
1383 {
1384  // Initialize all the interior FE objects on elem.
1385  // Logging of FE::reinit is done in the FE functions
1386  // We only reinit the FE objects for the current element
1387  // dimension
1388  const unsigned char dim = this->get_elem_dim();
1389 
1390  libmesh_assert( !_element_fe[dim].empty() );
1391 
1392  std::map<FEType, FEAbstract *>::iterator local_fe_end = _element_fe[dim].end();
1393  for (std::map<FEType, FEAbstract *>::iterator i = _element_fe[dim].begin();
1394  i != local_fe_end; ++i)
1395  {
1396  if (this->has_elem())
1397  i->second->reinit(&(this->get_elem()), pts);
1398  else
1399  // If !this->has_elem(), then we assume we are dealing with a SCALAR variable
1400  i->second->reinit(libmesh_nullptr);
1401  }
1402 }
1403 
1404 
1406 {
1407  // Initialize all the side FE objects on elem/side.
1408  // Logging of FE::reinit is done in the FE functions
1409  // We only reinit the FE objects for the current element
1410  // dimension
1411  const unsigned char dim = this->get_elem_dim();
1412 
1413  libmesh_assert( !_side_fe[dim].empty() );
1414 
1415  std::map<FEType, FEAbstract *>::iterator local_fe_end = _side_fe[dim].end();
1416  for (std::map<FEType, FEAbstract *>::iterator i = _side_fe[dim].begin();
1417  i != local_fe_end; ++i)
1418  {
1419  i->second->reinit(&(this->get_elem()), this->get_side());
1420  }
1421 }
1422 
1423 
1424 
1426 {
1427  libmesh_assert_equal_to (this->get_elem_dim(), 3);
1428 
1429  // Initialize all the interior FE objects on elem/edge.
1430  // Logging of FE::reinit is done in the FE functions
1431  std::map<FEType, FEAbstract *>::iterator local_fe_end = _edge_fe.end();
1432  for (std::map<FEType, FEAbstract *>::iterator i = _edge_fe.begin();
1433  i != local_fe_end; ++i)
1434  {
1435  i->second->edge_reinit(&(this->get_elem()), this->get_edge());
1436  }
1437 }
1438 
1439 
1440 
1442 {
1443  // This is too expensive to call unless we've been asked to move the mesh
1445 
1446  // This will probably break with threading when two contexts are
1447  // operating on elements which share a node
1448  libmesh_assert_equal_to (libMesh::n_threads(), 1);
1449 
1450  // If the coordinate data is in our own system, it's already
1451  // been set up for us
1452  // if (_mesh_sys == this->number())
1453  // {
1454  unsigned int n_nodes = this->get_elem().n_nodes();
1455 
1456 #ifndef NDEBUG
1457  const unsigned char dim = this->get_elem_dim();
1458 
1459  // For simplicity we demand that mesh coordinates be stored
1460  // in a format that allows a direct copy
1462  (this->get_element_fe(this->get_mesh_x_var(), dim)->get_fe_type().family
1463  == LAGRANGE &&
1464  this->get_element_fe(this->get_mesh_x_var(), dim)->get_fe_type().order.get_order()
1465  == this->get_elem().default_order()));
1467  (this->get_element_fe(this->get_mesh_y_var(), dim)->get_fe_type().family
1468  == LAGRANGE &&
1469  this->get_element_fe(this->get_mesh_y_var(), dim)->get_fe_type().order.get_order()
1470  == this->get_elem().default_order()));
1472  (this->get_element_fe(this->get_mesh_z_var(), dim)->get_fe_type().family
1473  == LAGRANGE &&
1474  this->get_element_fe(this->get_mesh_z_var(), dim)->get_fe_type().order.get_order()
1475  == this->get_elem().default_order()));
1476 #endif
1477 
1478  // Get degree of freedom coefficients from point coordinates
1479  if (this->get_mesh_x_var() != libMesh::invalid_uint)
1480  for (unsigned int i=0; i != n_nodes; ++i)
1481  (this->get_elem_solution(this->get_mesh_x_var()))(i) = this->get_elem().point(i)(0);
1482 
1483  if (this->get_mesh_y_var() != libMesh::invalid_uint)
1484  for (unsigned int i=0; i != n_nodes; ++i)
1485  (this->get_elem_solution(this->get_mesh_y_var()))(i) = this->get_elem().point(i)(1);
1486 
1487  if (this->get_mesh_z_var() != libMesh::invalid_uint)
1488  for (unsigned int i=0; i != n_nodes; ++i)
1489  (this->get_elem_solution(this->get_mesh_z_var()))(i) = this->get_elem().point(i)(2);
1490  // }
1491  // FIXME - If the coordinate data is not in our own system, someone
1492  // had better get around to implementing that... - RHS
1493  // else
1494  // {
1495  // libmesh_not_implemented();
1496  // }
1497 }
1498 
1499 
1500 
1501 // We can ignore the theta argument in the current use of this
1502 // function, because elem_subsolutions will already have been set to
1503 // the theta value.
1504 //
1505 // To enable loose mesh movement coupling things will need to change.
1507 {
1508  // This is too expensive to call unless we've been asked to move the mesh
1510 
1511  // This will probably break with threading when two contexts are
1512  // operating on elements which share a node
1513  libmesh_assert_equal_to (libMesh::n_threads(), 1);
1514 
1515  // If the coordinate data is in our own system, it's already
1516  // been set up for us, and we can ignore our input parameter theta
1517  // if (_mesh_sys == this->number())
1518  // {
1519  unsigned int n_nodes = this->get_elem().n_nodes();
1520 
1521 #ifndef NDEBUG
1522  const unsigned char dim = this->get_elem_dim();
1523 
1524  // For simplicity we demand that mesh coordinates be stored
1525  // in a format that allows a direct copy
1527  (this->get_element_fe(this->get_mesh_x_var(), dim)->get_fe_type().family
1528  == LAGRANGE &&
1529  this->get_elem_solution(this->get_mesh_x_var()).size() == n_nodes));
1531  (this->get_element_fe(this->get_mesh_y_var(), dim)->get_fe_type().family
1532  == LAGRANGE &&
1533  this->get_elem_solution(this->get_mesh_y_var()).size() == n_nodes));
1535  (this->get_element_fe(this->get_mesh_z_var(), dim)->get_fe_type().family
1536  == LAGRANGE &&
1537  this->get_elem_solution(this->get_mesh_z_var()).size() == n_nodes));
1538 #endif
1539 
1540  // Set the new point coordinates
1541  if (this->get_mesh_x_var() != libMesh::invalid_uint)
1542  for (unsigned int i=0; i != n_nodes; ++i)
1543  const_cast<Elem &>(this->get_elem()).point(i)(0) =
1544  libmesh_real(this->get_elem_solution(this->get_mesh_x_var())(i));
1545 
1546  if (this->get_mesh_y_var() != libMesh::invalid_uint)
1547  for (unsigned int i=0; i != n_nodes; ++i)
1548  const_cast<Elem &>(this->get_elem()).point(i)(1) =
1549  libmesh_real(this->get_elem_solution(this->get_mesh_y_var())(i));
1550 
1551  if (this->get_mesh_z_var() != libMesh::invalid_uint)
1552  for (unsigned int i=0; i != n_nodes; ++i)
1553  const_cast<Elem &>(this->get_elem()).point(i)(2) =
1554  libmesh_real(this->get_elem_solution(this->get_mesh_z_var())(i));
1555  // }
1556  // FIXME - If the coordinate data is not in our own system, someone
1557  // had better get around to implementing that... - RHS
1558  // else
1559  // {
1560  // libmesh_not_implemented();
1561  // }
1562 }
1563 
1564 
1565 
1566 
1567 
1568 /*
1569  void FEMContext::reinit(const FEMSystem & sys, Elem * e)
1570  {
1571  // Initialize our elem pointer, algebraic objects
1572  this->pre_fe_reinit(e);
1573 
1574  // Moving the mesh may be necessary
1575  // Reinitializing the FE objects is definitely necessary
1576  this->elem_reinit(1.);
1577  }
1578 */
1579 
1580 
1581 
1582 void FEMContext::pre_fe_reinit(const System & sys, const Elem * e)
1583 {
1584  this->set_elem(e);
1585 
1586  if (algebraic_type() == CURRENT ||
1588  {
1589  // Initialize the per-element data for elem.
1590  if (this->has_elem())
1591  sys.get_dof_map().dof_indices (&(this->get_elem()), this->get_dof_indices());
1592  else
1593  // If !this->has_elem(), then we assume we are dealing with a SCALAR variable
1594  sys.get_dof_map().dof_indices
1595  (static_cast<Elem*>(libmesh_nullptr), this->get_dof_indices());
1596  }
1597 #ifdef LIBMESH_ENABLE_AMR
1598  else if (algebraic_type() == OLD)
1599  {
1600  // Initialize the per-element data for elem.
1601  if (this->has_elem())
1602  sys.get_dof_map().old_dof_indices (&(this->get_elem()), this->get_dof_indices());
1603  else
1604  // If !this->has_elem(), then we assume we are dealing with a SCALAR variable
1606  (static_cast<Elem*>(libmesh_nullptr), this->get_dof_indices());
1607  }
1608 #endif // LIBMESH_ENABLE_AMR
1609 
1610  const unsigned int n_dofs = cast_int<unsigned int>
1611  (this->get_dof_indices().size());
1612  const std::size_t n_qoi = sys.qoi.size();
1613 
1614  if (this->algebraic_type() != NONE &&
1615  this->algebraic_type() != DOFS_ONLY)
1616  {
1617  // This also resizes elem_solution
1619  sys.current_local_solution->get(this->get_dof_indices(), this->get_elem_solution().get_values());
1620  else
1621  _custom_solution->get(this->get_dof_indices(), this->get_elem_solution().get_values());
1622 
1623  if (sys.use_fixed_solution)
1624  this->get_elem_fixed_solution().resize(n_dofs);
1625 
1626  // Only make space for these if we're using DiffSystem
1627  // This is assuming *only* DiffSystem is using elem_solution_rate/accel
1628  const DifferentiableSystem * diff_system = dynamic_cast<const DifferentiableSystem *>(&sys);
1629  if (diff_system)
1630  {
1631  // Now, we only need these if the solver is unsteady
1632  if (!diff_system->get_time_solver().is_steady())
1633  {
1634  this->get_elem_solution_rate().resize(n_dofs);
1635 
1636  // We only need accel space if the TimeSolver is second order
1637  const UnsteadySolver & time_solver = cast_ref<const UnsteadySolver &>(diff_system->get_time_solver());
1638 
1639  if (time_solver.time_order() >= 2 || !diff_system->get_second_order_vars().empty())
1640  this->get_elem_solution_accel().resize(n_dofs);
1641  }
1642  }
1643 
1644  if (algebraic_type() != OLD)
1645  {
1646  // These resize calls also zero out the residual and jacobian
1647  this->get_elem_residual().resize(n_dofs);
1648  this->get_elem_jacobian().resize(n_dofs, n_dofs);
1649 
1650  this->get_qoi_derivatives().resize(n_qoi);
1651  this->_elem_qoi_subderivatives.resize(n_qoi);
1652  for (std::size_t q=0; q != n_qoi; ++q)
1653  (this->get_qoi_derivatives())[q].resize(n_dofs);
1654  }
1655  }
1656 
1657  // Initialize the per-variable data for elem.
1658  {
1659  unsigned int sub_dofs = 0;
1660  for (unsigned int i=0; i != sys.n_vars(); ++i)
1661  {
1662  if (algebraic_type() == CURRENT ||
1664  {
1665  if (this->has_elem())
1666  sys.get_dof_map().dof_indices (&(this->get_elem()), this->get_dof_indices(i), i);
1667  else
1668  // If !this->has_elem(), then we assume we are dealing with a SCALAR variable
1669  sys.get_dof_map().dof_indices
1670  (static_cast<Elem*>(libmesh_nullptr), this->get_dof_indices(i), i);
1671  }
1672 #ifdef LIBMESH_ENABLE_AMR
1673  else if (algebraic_type() == OLD)
1674  {
1675  if (this->has_elem())
1676  sys.get_dof_map().old_dof_indices (&(this->get_elem()), this->get_dof_indices(i), i);
1677  else
1678  // If !this->has_elem(), then we assume we are dealing with a SCALAR variable
1680  (static_cast<Elem*>(libmesh_nullptr), this->get_dof_indices(i), i);
1681  }
1682 #endif // LIBMESH_ENABLE_AMR
1683 
1684  if (this->algebraic_type() != NONE &&
1685  this->algebraic_type() != DOFS_ONLY)
1686  {
1687  const unsigned int n_dofs_var = cast_int<unsigned int>
1688  (this->get_dof_indices(i).size());
1689 
1690  this->get_elem_solution(i).reposition
1691  (sub_dofs, n_dofs_var);
1692 
1693  // Only make space for these if we're using DiffSystem
1694  // This is assuming *only* DiffSystem is using elem_solution_rate/accel
1695  const DifferentiableSystem * diff_system = dynamic_cast<const DifferentiableSystem *>(&sys);
1696  if (diff_system)
1697  {
1698  // Now, we only need these if the solver is unsteady
1699  if (!diff_system->get_time_solver().is_steady())
1700  {
1701  this->get_elem_solution_rate(i).reposition
1702  (sub_dofs, n_dofs_var);
1703 
1704  // We only need accel space if the TimeSolver is second order
1705  const UnsteadySolver & time_solver = cast_ref<const UnsteadySolver &>(diff_system->get_time_solver());
1706 
1707  if (time_solver.time_order() >= 2 || !diff_system->get_second_order_vars().empty())
1708  this->get_elem_solution_accel(i).reposition
1709  (sub_dofs, n_dofs_var);
1710  }
1711  }
1712 
1713  if (sys.use_fixed_solution)
1714  this->get_elem_fixed_solution(i).reposition
1715  (sub_dofs, n_dofs_var);
1716 
1717  if (algebraic_type() != OLD)
1718  {
1719  this->get_elem_residual(i).reposition
1720  (sub_dofs, n_dofs_var);
1721 
1722  for (std::size_t q=0; q != n_qoi; ++q)
1723  this->get_qoi_derivatives(q,i).reposition
1724  (sub_dofs, n_dofs_var);
1725 
1726  for (unsigned int j=0; j != i; ++j)
1727  {
1728  const unsigned int n_dofs_var_j =
1729  cast_int<unsigned int>
1730  (this->get_dof_indices(j).size());
1731 
1732  this->get_elem_jacobian(i,j).reposition
1733  (sub_dofs, this->get_elem_residual(j).i_off(),
1734  n_dofs_var, n_dofs_var_j);
1735  this->get_elem_jacobian(j,i).reposition
1736  (this->get_elem_residual(j).i_off(), sub_dofs,
1737  n_dofs_var_j, n_dofs_var);
1738  }
1739  this->get_elem_jacobian(i,i).reposition
1740  (sub_dofs, sub_dofs,
1741  n_dofs_var,
1742  n_dofs_var);
1743  }
1744 
1745  sub_dofs += n_dofs_var;
1746  }
1747  }
1748 
1749  if (this->algebraic_type() != NONE &&
1750  this->algebraic_type() != DOFS_ONLY &&
1751  this->algebraic_type() != OLD)
1752  libmesh_assert_equal_to (sub_dofs, n_dofs);
1753  }
1754 
1755  // Now do the localization for the user requested vectors
1756  if (this->algebraic_type() != NONE)
1757  {
1758  DiffContext::localized_vectors_iterator localized_vec_it = this->_localized_vectors.begin();
1759  const DiffContext::localized_vectors_iterator localized_vec_end = this->_localized_vectors.end();
1760 
1761  for (; localized_vec_it != localized_vec_end; ++localized_vec_it)
1762  {
1763  const NumericVector<Number> & current_localized_vector = *localized_vec_it->first;
1764  DenseVector<Number> & target_vector = localized_vec_it->second.first;
1765 
1766  current_localized_vector.get(this->get_dof_indices(), target_vector.get_values());
1767 
1768  // Initialize the per-variable data for elem.
1769  unsigned int sub_dofs = 0;
1770  for (unsigned int i=0; i != sys.n_vars(); ++i)
1771  {
1772  const unsigned int n_dofs_var = cast_int<unsigned int>
1773  (this->get_dof_indices(i).size());
1774 
1775  // This is redundant with earlier initialization, isn't it? - RHS
1776  // sys.get_dof_map().dof_indices (&(this->get_elem()), this->get_dof_indices(i), i);
1777 
1778  localized_vec_it->second.second[i]->reposition
1779  (sub_dofs, n_dofs_var);
1780 
1781  sub_dofs += n_dofs_var;
1782  }
1783  libmesh_assert_equal_to (sub_dofs, n_dofs);
1784  }
1785  }
1786 }
1787 
1788 void FEMContext::set_elem( const Elem * e )
1789 {
1790  this->_elem = e;
1791 
1792  // If e is NULL, we assume it's SCALAR and set _elem_dim to 0.
1793  this->_elem_dim = this->_elem ? this->_elem->dim() : 0;
1794 }
1795 
1797 {
1798  // Update the "time" variable based on the value of theta. For this
1799  // to work, we need to know the value of deltat, a pointer to which is now
1800  // stored by our parent DiffContext class. Note: get_deltat_value() will
1801  // assert in debug mode if the requested pointer is NULL.
1802  const Real deltat = this->get_deltat_value();
1803 
1804  this->set_time(theta*(this->get_system_time() + deltat) + (1.-theta)*this->get_system_time());
1805 }
1806 
1807 
1808 
1809 template<>
1811 FEMContext::cached_fe( const unsigned int elem_dim,
1812  const FEType fe_type ) const
1813 {
1814 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
1815  const bool fe_needs_inf =
1816  this->has_elem() && this->get_elem().infinite();
1817 #endif
1818 
1819  if (!_real_fe ||
1820  elem_dim != _real_fe->get_dim() ||
1821  fe_type != _real_fe->get_fe_type())
1822  _real_fe =
1823 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
1824  fe_needs_inf ?
1825  FEGenericBase<Real>::build_InfFE(elem_dim, fe_type) :
1826 #endif
1827  FEGenericBase<Real>::build(elem_dim, fe_type);
1828 
1829 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
1830  else if (fe_needs_inf && !_real_fe_is_inf)
1831  _real_fe =
1832  FEGenericBase<Real>::build_InfFE(elem_dim, fe_type);
1833  else if (!fe_needs_inf && _real_fe_is_inf)
1834  _real_fe =
1835  FEGenericBase<Real>::build(elem_dim, fe_type);
1836 
1837  _real_fe_is_inf =
1838  (this->has_elem() && this->get_elem().infinite());
1839 #endif
1840 
1841  return _real_fe.get();
1842 }
1843 
1844 
1845 template<>
1847 FEMContext::cached_fe( const unsigned int elem_dim,
1848  const FEType fe_type ) const
1849 {
1850 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
1851  const bool fe_needs_inf =
1852  this->has_elem() && this->get_elem().infinite();
1853 #endif
1854 
1855  if (!_real_grad_fe ||
1856  elem_dim != _real_grad_fe->get_dim() ||
1857  fe_type != _real_grad_fe->get_fe_type())
1858  _real_grad_fe =
1859 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
1860  fe_needs_inf ?
1861  FEGenericBase<RealGradient>::build_InfFE(elem_dim, fe_type) :
1862 #endif
1863  FEGenericBase<RealGradient>::build(elem_dim, fe_type);
1864 
1865 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
1866  else if (fe_needs_inf && !_real_grad_fe_is_inf)
1867  _real_grad_fe =
1868  FEGenericBase<RealGradient>::build_InfFE(elem_dim, fe_type);
1869  else if (!fe_needs_inf && _real_grad_fe_is_inf)
1870  _real_grad_fe =
1871  FEGenericBase<RealGradient>::build(elem_dim, fe_type);
1872 
1874  (this->has_elem() && this->get_elem().infinite());
1875 #endif
1876 
1877  return _real_grad_fe.get();
1878 }
1879 
1880 
1881 
1882 template<typename OutputShape>
1885  const Point & p,
1886  const Real tolerance) const
1887 {
1888  FEType fe_type = fe->get_fe_type();
1889 
1890  // If we don't have an Elem to evaluate on, then the only functions
1891  // we can sensibly evaluate are the scalar dofs which are the same
1892  // everywhere.
1893  libmesh_assert(this->has_elem() || fe_type.family == SCALAR);
1894 
1895 #ifdef LIBMESH_ENABLE_AMR
1896  if ((algebraic_type() == OLD) &&
1897  this->has_elem())
1898  {
1899  if (this->get_elem().p_refinement_flag() == Elem::JUST_REFINED)
1900  fe_type.order = static_cast<Order>(fe_type.order - 1);
1901  else if (this->get_elem().p_refinement_flag() == Elem::JUST_COARSENED)
1902  fe_type.order = static_cast<Order>(fe_type.order + 1);
1903  }
1904 #endif // LIBMESH_ENABLE_AMR
1905 
1906  const unsigned int elem_dim = this->has_elem() ? this->get_elem().dim() : 0;
1907 
1908  FEGenericBase<OutputShape>* fe_new = cached_fe<OutputShape>(elem_dim, fe_type);
1909 
1910  // Map the physical co-ordinates to the master co-ordinates using the inverse_map from fe_interface.h
1911  // Build a vector of point co-ordinates to send to reinit
1912  Point master_point = this->has_elem() ?
1913  FEInterface::inverse_map (elem_dim,
1914  fe_type,
1915  &this->get_elem(),
1916  p,
1917  tolerance) : Point(0);
1918 
1919  std::vector<Point> coor(1, master_point);
1920 
1921  // Reinitialize the element and compute the shape function values at coor
1922  if (this->has_elem())
1923  fe_new->reinit (&this->get_elem(), &coor);
1924  else
1925  // If !this->has_elem(), then we assume we are dealing with a SCALAR variable
1926  fe_new->reinit (libmesh_nullptr, &coor);
1927 
1928  return fe_new;
1929 }
1930 
1931 
1932 
1933 
1934 
1935 // Instantiate member function templates
1936 template void FEMContext::interior_value<Number>(unsigned int, unsigned int, Number &) const;
1937 template void FEMContext::interior_values<Number>(unsigned int, const NumericVector<Number> &,
1938  std::vector<Number> &) const;
1939 template void FEMContext::interior_value<Gradient>(unsigned int, unsigned int, Gradient &) const;
1940 template void FEMContext::interior_values<Gradient>(unsigned int, const NumericVector<Number> &,
1941  std::vector<Gradient> &) const;
1942 
1943 template void FEMContext::interior_gradient<Gradient>(unsigned int, unsigned int, Gradient &) const;
1944 template void FEMContext::interior_gradients<Gradient>(unsigned int, const NumericVector<Number> &,
1945  std::vector<Gradient> &) const;
1946 template void FEMContext::interior_gradient<Tensor>(unsigned int, unsigned int, Tensor &) const;
1947 template void FEMContext::interior_gradients<Tensor>(unsigned int, const NumericVector<Number> &,
1948  std::vector<Tensor> &) const;
1949 
1950 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1951 template void FEMContext::interior_hessian<Tensor>(unsigned int, unsigned int, Tensor &) const;
1952 template void FEMContext::interior_hessians<Tensor>(unsigned int, const NumericVector<Number> &,
1953  std::vector<Tensor> &) const;
1954 //FIXME: Not everything is implemented yet for second derivatives of RealGradients
1955 //template void FEMContext::interior_hessian<??>(unsigned int, unsigned int, ??&) const;
1956 //template void FEMContext::interior_hessians<??>(unsigned int, const NumericVector<Number> &,
1957 // std::vector<??> &) const;
1958 #endif
1959 
1960 template void FEMContext::interior_curl<Gradient>(unsigned int, unsigned int, Gradient &) const;
1961 
1962 template void FEMContext::interior_div<Number>(unsigned int, unsigned int, Number &) const;
1963 
1964 template void FEMContext::side_value<Number>(unsigned int, unsigned int, Number &) const;
1965 template void FEMContext::side_value<Gradient>(unsigned int, unsigned int, Gradient &) const;
1966 template void FEMContext::side_values<Number>(unsigned int, const NumericVector<Number> &,
1967  std::vector<Number> &) const;
1968 template void FEMContext::side_values<Gradient>(unsigned int, const NumericVector<Number> &,
1969  std::vector<Gradient> &) const;
1970 
1971 template void FEMContext::side_gradient<Gradient>(unsigned int, unsigned int, Gradient &) const;
1972 template void FEMContext::side_gradients<Gradient>(unsigned int, const NumericVector<Number> &,
1973  std::vector<Gradient> &) const;
1974 template void FEMContext::side_gradient<Tensor>(unsigned int, unsigned int, Tensor &) const;
1975 template void FEMContext::side_gradients<Tensor>(unsigned int, const NumericVector<Number> &,
1976  std::vector<Tensor> &) const;
1977 
1978 
1979 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1980 template void FEMContext::side_hessian<Tensor>(unsigned int, unsigned int, Tensor &) const;
1981 template void FEMContext::side_hessians<Tensor>(unsigned int, const NumericVector<Number> &,
1982  std::vector<Tensor> &) const;
1983 //FIXME: Not everything is implemented yet for second derivatives of RealGradients
1984 //template void FEMContext::side_hessian<??>(unsigned int, unsigned int,
1985 // ??&) const;
1986 //template void FEMContext::side_hessians<??>(unsigned int, const NumericVector<Number> &,
1987 // std::vector<??> &) const;
1988 #endif
1989 
1990 template void FEMContext::point_value<Number>(unsigned int, const Point &, Number &, const Real) const;
1991 template void FEMContext::point_value<Gradient>(unsigned int, const Point &, Gradient &, const Real) const;
1992 
1993 template void FEMContext::point_gradient<Gradient>(unsigned int, const Point &, Gradient &, const Real) const;
1994 template void FEMContext::point_gradient<Tensor>(unsigned int, const Point &, Tensor &, const Real) const;
1995 
1996 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1997 template void FEMContext::point_hessian<Tensor>(unsigned int, const Point &, Tensor &, const Real) const;
1998 //FIXME: Not everything is implemented yet for second derivatives of RealGradients
1999 //template void FEMContext::point_hessian<??>(unsigned int, const Point &, ??&) const;
2000 #endif
2001 
2002 template void FEMContext::point_curl<Gradient>(unsigned int, const Point &, Gradient &, const Real) const;
2003 
2004 template void FEMContext::fixed_interior_value<Number>(unsigned int, unsigned int, Number &) const;
2005 template void FEMContext::fixed_interior_value<Gradient>(unsigned int, unsigned int, Gradient &) const;
2006 
2007 template void FEMContext::fixed_interior_gradient<Gradient>(unsigned int, unsigned int, Gradient &) const;
2008 template void FEMContext::fixed_interior_gradient<Tensor>(unsigned int, unsigned int, Tensor &) const;
2009 
2010 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
2011 template void FEMContext::fixed_interior_hessian<Tensor>(unsigned int, unsigned int, Tensor &) const;
2012 //FIXME: Not everything is implemented yet for second derivatives of RealGradients
2013 //template void FEMContext::fixed_interior_hessian<??>(unsigned int, unsigned int, ??&) const;
2014 #endif
2015 
2016 template void FEMContext::fixed_side_value<Number>(unsigned int, unsigned int, Number &) const;
2017 template void FEMContext::fixed_side_value<Gradient>(unsigned int, unsigned int, Gradient &) const;
2018 
2019 template void FEMContext::fixed_side_gradient<Gradient>(unsigned int, unsigned int, Gradient &) const;
2020 template void FEMContext::fixed_side_gradient<Tensor>(unsigned int, unsigned int, Tensor &) const;
2021 
2022 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
2023 template void FEMContext::fixed_side_hessian<Tensor>(unsigned int, unsigned int, Tensor &) const;
2024 //FIXME: Not everything is implemented yet for second derivatives of RealGradients
2025 //template void FEMContext::fixed_side_hessian<??>(unsigned int, unsigned int, ??&) const;
2026 #endif
2027 
2028 template void FEMContext::fixed_point_value<Number>(unsigned int, const Point &, Number &, const Real) const;
2029 template void FEMContext::fixed_point_value<Gradient>(unsigned int, const Point &, Gradient &, const Real) const;
2030 
2031 template void FEMContext::fixed_point_gradient<Gradient>(unsigned int, const Point &, Gradient &, const Real) const;
2032 template void FEMContext::fixed_point_gradient<Tensor>(unsigned int, const Point &, Tensor &, const Real) const;
2033 
2034 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
2035 template void FEMContext::fixed_point_hessian<Tensor>(unsigned int, const Point &, Tensor &, const Real) const;
2036 //FIXME: Not everything is implemented yet for second derivatives of RealGradients
2037 //template void FEMContext::fixed_point_hessian<??>(unsigned int, const Point &, ??&) const;
2038 #endif
2039 
2040 template void FEMContext::interior_rate<Number>(unsigned int, unsigned int, Number &) const;
2041 template void FEMContext::interior_rate<Gradient>(unsigned int, unsigned int, Gradient &) const;
2042 
2043 template void FEMContext::side_rate<Number>(unsigned int, unsigned int, Number &) const;
2044 template void FEMContext::side_rate<Gradient>(unsigned int, unsigned int, Gradient &) const;
2045 
2046 template void FEMContext::interior_accel<Number>(unsigned int, unsigned int, Number &) const;
2047 template void FEMContext::interior_accel<Gradient>(unsigned int, unsigned int, Gradient &) const;
2048 
2049 template void FEMContext::side_accel<Number>(unsigned int, unsigned int, Number &) const;
2050 template void FEMContext::side_accel<Gradient>(unsigned int, unsigned int, Gradient &) const;
2051 
2052 template FEGenericBase<Real> *
2054  const Point &,
2055  const Real) const;
2056 
2057 template FEGenericBase<RealGradient> *
2059  const Point &,
2060  const Real) const;
2061 
2062 } // namespace libMesh
T libmesh_real(T a)
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition: fe_type.h:178
FEFamily family
The type of finite element.
Definition: fe_type.h:203
Number fixed_interior_value(unsigned int var, unsigned int qp) const
Definition: fem_context.C:973
unsigned char get_side() const
Accessor for current side of Elem object.
Definition: fem_context.h:885
const std::set< unsigned int > & get_second_order_vars() const
Definition: diff_physics.h:529
void interior_curl(unsigned int var, unsigned int qp, OutputType &curl_u) const
Definition: fem_context.C:512
unsigned int n_threads()
Definition: libmesh_base.h:125
const DenseVector< Number > & get_elem_solution() const
Accessor for element solution.
Definition: diff_context.h:110
void elem_position_set(Real theta)
Uses the coordinate data specified by mesh_*_position configuration to set the geometry of elem to th...
Definition: fem_context.h:1195
This class provides all data required for a physics package (e.g.
Definition: diff_context.h:54
virtual void pre_fe_reinit(const System &, const Elem *e)
Reinitializes local data vectors/matrices on the current geometric element.
Definition: fem_context.C:1582
const unsigned int invalid_uint
A number which is used quite often to represent an invalid or uninitialized value.
Definition: libmesh.h:184
unsigned int get_mesh_x_var() const
Accessor for x-variable of moving mesh System.
Definition: fem_context.h:823
std::map< const NumericVector< Number > *, std::pair< DenseVector< Number >, std::vector< DenseSubVector< Number > * > > > _localized_vectors
Contains pointers to vectors the user has asked to be localized, keyed with pairs of element localize...
Definition: diff_context.h:558
const NumericVector< Number > * _custom_solution
Data with which to do algebra reinitialization.
Definition: fem_context.h:994
Tensor interior_hessian(unsigned int var, unsigned int qp) const
Definition: fem_context.C:446
void set_elem(const Elem *e)
Helper function to promote accessor usage.
Definition: fem_context.C:1788
unsigned int _mesh_z_var
Definition: fem_context.h:972
void side_values(unsigned int var, const NumericVector< Number > &_system_vector, std::vector< OutputType > &side_values_vector) const
Fills a vector of values of the _system_vector at the all the quadrature points on the current elemen...
Definition: fem_context.C:598
void point_curl(unsigned int var, const Point &p, OutputType &curl_u, const Real tolerance=TOLERANCE) const
Definition: fem_context.C:937
Number point_value(unsigned int var, const Point &p) const
Definition: fem_context.C:788
unsigned int dim
Tensor fixed_point_hessian(unsigned int var, const Point &p) const
Definition: fem_context.C:1220
void side_hessians(unsigned int var, const NumericVector< Number > &_system_vector, std::vector< OutputType > &d2u_vals) const
Fills a vector of hessians of the _system_vector at the all the quadrature points on the current elem...
Definition: fem_context.C:745
ImplicitSystem & sys
bool has_elem() const
Test for current Elem object.
Definition: fem_context.h:865
std::vector< boundary_id_type > side_boundary_ids() const
Lists the boundary ids found on the current side.
Definition: fem_context.C:216
TensorTools::MakeReal< OutputType >::type value_shape
Definition: fem_context.h:1042
std::vector< std::vector< DenseSubVector< Number > * > > _elem_qoi_subderivatives
Definition: diff_context.h:609
const DenseVector< Number > & get_elem_solution_accel() const
Accessor for element solution accel of change w.r.t.
Definition: diff_context.h:179
void interior_hessians(unsigned int var, const NumericVector< Number > &_system_vector, std::vector< OutputType > &d2u_vals) const
Fills a vector of hessians of the _system_vector at the all the quadrature points in the current elem...
Definition: fem_context.C:470
UniquePtr< FEGenericBase< Real > > _real_fe
Definition: fem_context.h:1004
void resize(const unsigned int n)
Resize the vector.
Definition: dense_vector.h:350
std::vector< T > & get_values()
Definition: dense_vector.h:251
unsigned short int side
Definition: xdr_io.C:49
const std::vector< dof_id_type > & get_dof_indices() const
Accessor for element dof indices.
Definition: diff_context.h:366
This is the base class from which all geometric element types are derived.
Definition: elem.h:89
void some_value(unsigned int var, unsigned int qp, OutputType &u) const
Helper function to reduce some code duplication in the *interior_value methods.
Definition: fem_context.C:234
const DenseVector< Number > & get_elem_fixed_solution() const
Accessor for element fixed solution.
Definition: diff_context.h:214
UniquePtr< NumericVector< Number > > current_local_solution
All the values I need to compute my contribution to the simulation at hand.
Definition: system.h:1535
const class libmesh_nullptr_t libmesh_nullptr
std::vector< QBase * > _side_qrule
Quadrature rules for element sides The FEM context will try to find a quadrature rule that correctly ...
Definition: fem_context.h:1149
OrderWrapper order
The approximation order of the element.
Definition: fe_type.h:197
void get_element_fe(unsigned int var, FEGenericBase< OutputShape > *&fe) const
Accessor for interior finite element object for variable var for the largest dimension in the mesh...
Definition: fem_context.h:262
int _extra_quadrature_order
The extra quadrature order for this context.
Definition: fem_context.h:1163
This class defines a vector in LIBMESH_DIM dimensional Real or Complex space.
void old_dof_indices(const Elem *const elem, std::vector< dof_id_type > &di, const unsigned int vn=libMesh::invalid_uint) const
After a mesh is refined and repartitioned it is possible that the _send_list will need to be augmente...
Definition: dof_map.C:2378
void elem_position_get()
Uses the geometry of elem to set the coordinate data specified by mesh_*_position configuration...
Definition: fem_context.C:1441
The libMesh namespace provides an interface to certain functionality in the library.
void some_hessian(unsigned int var, unsigned int qp, OutputType &u) const
Helper function to reduce some code duplication in the *interior_hessian methods. ...
Definition: fem_context.C:296
const DenseVector< Number > & get_elem_residual() const
Const accessor for element residual.
Definition: diff_context.h:248
const FEType & variable_type(const unsigned int i) const
Definition: system.h:2164
virtual ~FEMContext()
Destructor.
Definition: fem_context.C:175
FEType get_fe_type() const
Definition: fe_abstract.h:421
Number fixed_side_value(unsigned int var, unsigned int qp) const
Definition: fem_context.C:1049
const Elem & get_elem() const
Accessor for current Elem object.
Definition: fem_context.h:871
Number interior_value(unsigned int var, unsigned int qp) const
Definition: fem_context.C:326
void side_gradients(unsigned int var, const NumericVector< Number > &_system_vector, std::vector< OutputType > &side_gradients_vector) const
Fills a vector with the gradient of the solution variable var at all the quadrature points on the cur...
Definition: fem_context.C:678
void interior_rate(unsigned int var, unsigned int qp, OutputType &u) const
Definition: fem_context.C:1273
libmesh_assert(j)
virtual void nonlocal_reinit(Real theta) libmesh_override
Gives derived classes the opportunity to reinitialize data needed for nonlocal calculations at a new ...
Definition: fem_context.C:1372
virtual void elem_side_reinit(Real theta) libmesh_override
Resets the current time in the context.
Definition: fem_context.C:1342
Defines a dense subvector for use in finite element computations.
Gradient fixed_point_gradient(unsigned int var, const Point &p) const
Definition: fem_context.C:1169
UniquePtr< FEGenericBase< RealGradient > > _real_grad_fe
Definition: fem_context.h:1005
virtual bool infinite() const =0
virtual unsigned int n_nodes() const =0
This class provides a specific system class.
Definition: diff_system.h:53
std::vector< boundary_id_type > boundary_ids(const Node *node) const
unsigned char edge
Current edge for edge_* to examine.
Definition: fem_context.h:982
const DenseMatrix< Number > & get_elem_jacobian() const
Const accessor for element Jacobian.
Definition: diff_context.h:282
const std::set< unsigned char > & elem_dimensions() const
Definition: fem_context.h:913
std::set< unsigned char > _elem_dims
Cached dimensions of elements in the mesh, plus dimension 0 if SCALAR variables are in use...
Definition: fem_context.h:1133
const dof_id_type n_nodes
Definition: tecplot_io.C:67
virtual bool is_steady() const =0
Is this effectively a steady-state solver?
const std::vector< DenseVector< Number > > & get_qoi_derivatives() const
Const accessor for QoI derivatives.
Definition: diff_context.h:330
const DofMap & get_dof_map() const
Definition: system.h:2030
int8_t boundary_id_type
Definition: id_types.h:51
virtual unsigned int time_order() const =0
unsigned char _elem_dim
Cached dimension of this->_elem.
Definition: fem_context.h:1127
std::map< const NumericVector< Number > *, std::pair< DenseVector< Number >, std::vector< DenseSubVector< Number > * > > >::iterator localized_vectors_iterator
Typedef for the localized_vectors iterator.
Definition: diff_context.h:527
virtual void side_fe_reinit()
Reinitializes side FE objects on the current geometric element.
Definition: fem_context.C:1405
std::vector< DenseSubVector< Number > * > _elem_fixed_subsolutions
Definition: diff_context.h:587
std::vector< Number > qoi
Values of the quantities of interest.
Definition: system.h:1553
void _update_time_from_system(Real theta)
Update the time in the context object for the given value of theta, based on the values of "time" and...
Definition: fem_context.C:1796
void init_internal_data(const System &sys)
Helper function used in constructors to set up internal data.
Definition: fem_context.C:79
FEGenericBase< OutputShape > * cached_fe(const unsigned int elem_dim, const FEType fe_type) const
bool use_fixed_solution
A boolean to be set to true by systems using elem_fixed_solution, for optional use by e...
Definition: system.h:1493
This is the base class for classes which contain information related to any physical process that mig...
Definition: system.h:76
System * _mesh_sys
System from which to acquire moving mesh information.
Definition: fem_context.h:967
Tensor fixed_side_hessian(unsigned int var, unsigned int qp) const
Definition: fem_context.C:1098
const std::vector< std::vector< OutputShape > > & get_phi() const
Definition: fe_base.h:208
void interior_accel(unsigned int var, unsigned int qp, OutputType &u) const
Definition: fem_context.C:1295
unsigned int _mesh_x_var
Variables from which to acquire moving mesh information.
Definition: fem_context.h:972
AlgebraicType _atype
Keep track of what type of algebra reinitialization is to be done.
Definition: fem_context.h:989
static Point inverse_map(const unsigned int dim, const FEType &fe_t, const Elem *elem, const Point &p, const Real tolerance=TOLERANCE, const bool secure=true)
Definition: fe_interface.C:569
Tensor point_hessian(unsigned int var, const Point &p) const
Definition: fem_context.C:886
RefinementState p_refinement_flag() const
Definition: elem.h:2521
std::map< FEType, FEAbstract * > _edge_fe
Definition: fem_context.h:1095
virtual Order default_order() const =0
const BoundaryInfo & _boundary_info
Saved reference to BoundaryInfo on the mesh for this System.
Definition: fem_context.h:1112
void set_time(Real time_in)
Set the time for which the current nonlinear_solution is defined.
Definition: diff_context.h:410
std::vector< std::vector< FEAbstract * > > _element_fe_var
Pointers to the same finite element objects, but indexed by variable number.
Definition: fem_context.h:1104
Gradient interior_gradient(unsigned int var, unsigned int qp) const
Definition: fem_context.C:381
static UniquePtr< FEGenericBase > build_InfFE(const unsigned int dim, const FEType &type)
Builds a specific infinite element type.
virtual void elem_fe_reinit(const std::vector< Point > *const pts=libmesh_nullptr)
Reinitializes interior FE objects on the current geometric element.
Definition: fem_context.C:1382
unsigned char get_edge() const
Accessor for current edge of Elem object.
Definition: fem_context.h:891
DenseSubVector< Number > & get_localized_subvector(const NumericVector< Number > &localized_vector, unsigned int var)
Return a reference to DenseSubVector localization of localized_vector at variable var contained in th...
Definition: diff_context.C:184
void interior_gradients(unsigned int var, const NumericVector< Number > &_system_vector, std::vector< OutputType > &interior_gradients_vector) const
Fills a vector with the gradient of the solution variable var at all the quadrature points in the cur...
Definition: fem_context.C:408
const DenseVector< Number > & get_elem_solution_rate() const
Accessor for element solution rate of change w.r.t.
Definition: diff_context.h:144
unsigned char _dim
Cached dimension of largest dimension element in this mesh.
Definition: fem_context.h:1122
virtual void edge_fe_reinit()
Reinitializes edge FE objects on the current geometric element.
Definition: fem_context.C:1425
std::vector< std::map< FEType, FEAbstract * > > _side_fe
Definition: fem_context.h:1094
std::vector< FEAbstract * > _edge_fe_var
Definition: fem_context.h:1106
AlgebraicType algebraic_type() const
Definition: fem_context.h:952
void side_rate(unsigned int var, unsigned int qp, OutputType &u) const
Definition: fem_context.C:1285
std::vector< std::map< FEType, FEAbstract * > > _element_fe
Finite element objects for each variable&#39;s interior, sides and edges.
Definition: fem_context.h:1093
bool has_side_boundary_id(boundary_id_type id) const
Reports if the boundary id is found on the current side.
Definition: fem_context.C:209
std::vector< QBase * > _element_qrule
Quadrature rule for element interior.
Definition: fem_context.h:1141
This is a generic class that defines a solver to handle time integration of DifferentiableSystems.
virtual void elem_reinit(Real theta) libmesh_override
Resets the current time in the context.
Definition: fem_context.C:1318
const std::vector< std::vector< OutputGradient > > & get_dphi() const
Definition: fe_base.h:216
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
const Point & point(const unsigned int i) const
Definition: elem.h:1809
unsigned char side
Current side for side_* to examine.
Definition: fem_context.h:977
const std::vector< std::vector< OutputTensor > > & get_d2phi() const
Definition: fe_base.h:290
Number side_value(unsigned int var, unsigned int qp) const
Definition: fem_context.C:575
UniquePtr< QBase > default_quadrature_rule(const unsigned int dim, const int extraorder=0) const
Definition: fe_type.C:30
virtual void get(const std::vector< numeric_index_type > &index, T *values) const
Access multiple components at once.
bool has_boundary_id(const Node *const node, const boundary_id_type id) const
Gradient side_gradient(unsigned int var, unsigned int qp) const
Definition: fem_context.C:633
void resize(const unsigned int new_m, const unsigned int new_n)
Resize the matrix.
Definition: dense_matrix.h:776
Gradient point_gradient(unsigned int var, const Point &p) const
Definition: fem_context.C:834
virtual unsigned int dim() const =0
Order
defines an enum for polynomial orders.
Definition: enum_order.h:32
FEMContext(const System &sys)
Constructor.
Definition: fem_context.C:37
Helper nested class for C++03-compatible "template typedef".
Definition: fem_context.h:1035
void some_gradient(unsigned int var, unsigned int qp, OutputType &u) const
Helper function to reduce some code duplication in the *interior_gradient methods.
Definition: fem_context.C:263
unsigned int _mesh_y_var
Definition: fem_context.h:972
Tensor side_hessian(unsigned int var, unsigned int qp) const
Definition: fem_context.C:716
unsigned int n_vars() const
Definition: system.h:2086
void side_accel(unsigned int var, unsigned int qp, OutputType &u) const
Definition: fem_context.C:1307
FEGenericBase< OutputShape > * build_new_fe(const FEGenericBase< OutputShape > *fe, const Point &p, const Real tolerance=TOLERANCE) const
Helper function to reduce some code duplication in the *_point_* methods.
Definition: fem_context.C:1884
unsigned int get_mesh_z_var() const
Accessor for z-variable of moving mesh System.
Definition: fem_context.h:851
static UniquePtr< FEAbstract > build(const unsigned int dim, const FEType &type)
Builds a specific finite element type.
Definition: fe_abstract.C:44
std::vector< DenseSubVector< Number > * > _elem_subsolutions
Definition: diff_context.h:565
std::vector< std::vector< FEAbstract * > > _side_fe_var
Definition: fem_context.h:1105
const DenseSubVector< Number > &(DiffContext::* diff_subsolution_getter)(unsigned int) const
Helper typedef to simplify refactoring.
Definition: fem_context.h:1028
virtual void elem_edge_reinit(Real theta) libmesh_override
Resets the current time in the context.
Definition: fem_context.C:1357
Gradient fixed_interior_gradient(unsigned int var, unsigned int qp) const
Definition: fem_context.C:996
const Elem * _elem
Current element for element_* to examine.
Definition: fem_context.h:1117
void interior_values(unsigned int var, const NumericVector< Number > &_system_vector, std::vector< OutputType > &interior_values_vector) const
Fills a vector of values of the _system_vector at the all the quadrature points in the current elemen...
Definition: fem_context.C:346
void _do_elem_position_set(Real theta)
Uses the coordinate data specified by mesh_*_position configuration to set the geometry of elem to th...
Definition: fem_context.C:1506
void interior_div(unsigned int var, unsigned int qp, OutputType &div_u) const
Definition: fem_context.C:543
void get_side_fe(unsigned int var, FEGenericBase< OutputShape > *&fe) const
Accessor for edge/face (2D/3D) finite element object for variable var for the largest dimension in th...
Definition: fem_context.h:299
Number fixed_point_value(unsigned int var, const Point &p) const
Definition: fem_context.C:1123
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:38
Gradient fixed_side_gradient(unsigned int var, unsigned int qp) const
Definition: fem_context.C:1073
Tensor fixed_interior_hessian(unsigned int var, unsigned int qp) const
Definition: fem_context.C:1023
UniquePtr< QBase > _edge_qrule
Quadrature rules for element edges.
Definition: fem_context.h:1158
void ErrorVector unsigned int
Definition: adjoints_ex3.C:360
This class forms the foundation from which generic finite elements may be derived.
unsigned char get_elem_dim() const
Definition: fem_context.h:906
unsigned int get_mesh_y_var() const
Accessor for y-variable of moving mesh System.
Definition: fem_context.h:837
Real get_system_time() const
Accessor for the time variable stored in the system class.
Definition: diff_context.h:398
This class defines a tensor in LIBMESH_DIM dimensional Real or Complex space.
virtual void reinit(const Elem *elem, const std::vector< Point > *const pts=libmesh_nullptr, const std::vector< Real > *const weights=libmesh_nullptr)=0
This is at the core of this class.
TimeSolver & get_time_solver()
Definition: diff_system.h:402
void dof_indices(const Elem *const elem, std::vector< dof_id_type > &di) const
Fills the vector di with the global degree of freedom indices for the element.
Definition: dof_map.C:1917
static UniquePtr< FEGenericBase > build(const unsigned int dim, const FEType &type)
Builds a specific finite element type.