<commit_before>// The libMesh Finite Element Library.
// Copyright (C) 2002-2014 <NAME>, <NAME>, <NAME>

// This library is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.

// This library is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// Lesser General Public License for more details.

// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA



#include "libmesh/boundary_info.h"
#include "libmesh/dof_map.h"
#include "libmesh/elem.h"
#include "libmesh/fe_base.h"
#include "libmesh/fe_interface.h"
#include "libmesh/fem_context.h"
#include "libmesh/libmesh_logging.h"
#include "libmesh/mesh_base.h"
#include "libmesh/quadrature.h"
#include "libmesh/system.h"
#include "libmesh/time_solver.h"
#include "libmesh/unsteady_solver.h" // For euler_residual

namespace libMesh
{

FEMContext::FEMContext (const System &sys)
  : DiffContext(sys),
    _mesh_sys(NULL),
    _mesh_x_var(0),
    _mesh_y_var(0),
    _mesh_z_var(0),
    side(0), edge(0),
    _boundary_info(sys.get_mesh().get_boundary_info()),
    _elem(NULL),
    _dim(sys.get_mesh().mesh_dimension()),
    _elem_dim(0), /* This will be reset in set_elem(). */
    _elem_dims(sys.get_mesh().elem_dimensions()),
    _element_qrule(4,NULL),
    _side_qrule(4,NULL),
    _edge_qrule(NULL)
{
  // Reserve space for the FEAbstract and QBase objects for each
  // element dimension possiblity (0,1,2,3)
  _element_fe.resize(4);
  _side_fe.resize(4);
  _element_fe_var.resize(4);
  _side_fe_var.resize(4);

  // We need to know which of our variables has the hardest
  // shape functions to numerically integrate.

  unsigned int nv = sys.n_vars();

  libmesh_assert (nv);
  FEType hardest_fe_type = sys.variable_type(0);

  bool have_scalar = false;

  for (unsigned int i=0; i != nv; ++i)
    {
      FEType fe_type = sys.variable_type(i);

      // FIXME - we don't yet handle mixed finite elements from
      // different families which require different quadrature rules
      // libmesh_assert_equal_to (fe_type.family, hardest_fe_type.family);

      if (fe_type.order > hardest_fe_type.order)
        hardest_fe_type = fe_type;

      // We need to detect SCALAR's so we can prepare FE objects for them.
      if( fe_type.family == SCALAR )
        have_scalar = true;
    }

  if(have_scalar)
    // SCALAR FEs have dimension 0 by assumption
    _elem_dims.insert(0);

  for( std::set<unsigned char>::const_iterator dim_it = _elem_dims.begin();
       dim_it != _elem_dims.end(); ++dim_it )
    {
      const unsigned char dim = *dim_it;

      // Create an adequate quadrature rule
      _element_qrule[dim] = hardest_fe_type.default_quadrature_rule
        (dim, sys.extra_quadrature_order).release();
      _side_qrule[dim] = hardest_fe_type.default_quadrature_rule
        (dim-1, sys.extra_quadrature_order).release();
      if (dim == 3)
        _edge_qrule = hardest_fe_type.default_quadrature_rule
          (1, sys.extra_quadrature_order).release();

      // Next, create finite element objects
      _element_fe_var[dim].resize(nv);
      _side_fe_var[dim].resize(nv);
      if (dim == 3)
        _edge_fe_var.resize(nv);


      for (unsigned int i=0; i != nv; ++i)
        {
          FEType fe_type = sys.variable_type(i);

          if ( _element_fe[dim][fe_type] == NULL )
            {
              _element_fe[dim][fe_type] = FEAbstract::build(dim, fe_type).release();
              _element_fe[dim][fe_type]->attach_quadrature_rule(_element_qrule[dim]);
              _side_fe[dim][fe_type] = FEAbstract::build(dim, fe_type).release();
              _side_fe[dim][fe_type]->attach_quadrature_rule(_side_qrule[dim]);

              if (this->_dim == 3)
                {
                  _edge_fe[fe_type] = FEAbstract::build(dim, fe_type).release();
                  _edge_fe[fe_type]->attach_quadrature_rule(_edge_qrule);
                }
            }

          _element_fe_var[dim][i] = _element_fe[dim][fe_type];
          _side_fe_var[dim][i] = _side_fe[dim][fe_type];
          if ((dim) == 3)
            _edge_fe_var[i] = _edge_fe[fe_type];

        }
    }
}

FEMContext::~FEMContext()
{
  // We don't want to store UniquePtrs in STL containers, but we don't
  // want to leak memory either
  for (std::vector<std::map<FEType, FEAbstract *> >::iterator d = _element_fe.begin();
       d != _element_fe.end(); ++d)
    for (std::map<FEType, FEAbstract *>::iterator i = d->begin();
         i != d->end(); ++i)
      delete i->second;

  for (std::vector<std::map<FEType, FEAbstract *> >::iterator d = _side_fe.begin();
       d != _side_fe.end(); ++d)
    for (std::map<FEType, FEAbstract *>::iterator i = d->begin();
         i != d->end(); ++i)
      delete i->second;

  for (std::map<FEType, FEAbstract *>::iterator i = _edge_fe.begin();
       i != _edge_fe.end(); ++i)
    delete i->second;
  _edge_fe.clear();

  for (std::vector<QBase*>::iterator i = _element_qrule.begin();
       i != _element_qrule.end(); ++i)
    delete *i;
  _element_qrule.clear();

  for (std::vector<QBase*>::iterator i = _side_qrule.begin();
       i != _side_qrule.end(); ++i)
    delete *i;
  _side_qrule.clear();

  delete _edge_qrule;
  _edge_qrule = NULL;
}



bool FEMContext::has_side_boundary_id(boundary_id_type id) const
{
  return _boundary_info.has_boundary_id(&(this->get_elem()), side, id);
}


std::vector<boundary_id_type> FEMContext::side_boundary_ids() const
{
  return _boundary_info.boundary_ids(&(this->get_elem()), side);
}



template<typename OutputType,
         typename FEMContext::FENeeded<OutputType>::value_getter fe_getter,
         FEMContext::diff_subsolution_getter subsolution_getter>
void FEMContext::some_value(unsigned int var, unsigned int qp, OutputType& u) const
{
  // Get local-to-global dof index lookup
  libmesh_assert_greater (this->get_dof_indices().size(), var);
  const unsigned int n_dofs = cast_int<unsigned int>
    (this->get_dof_indices(var).size());

  // Get current local coefficients
  const DenseSubVector<Number> &coef = (this->*subsolution_getter)(var);

  // Get finite element object
  typename FENeeded<OutputType>::value_base* fe = NULL;
  (this->*fe_getter)( var, fe );

  // Get shape function values at quadrature point
  const std::vector<std::vector
    <typename FENeeded<OutputType>::value_shape> > &phi = fe->get_phi();

  // Accumulate solution value
  u = 0.;

  for (unsigned int l=0; l != n_dofs; l++)
    u += phi[l][qp] * coef(l);
}



template<typename OutputType,
         typename FEMContext::FENeeded<OutputType>::grad_getter fe_getter,
         FEMContext::diff_subsolution_getter subsolution_getter>
void FEMContext::some_gradient(unsigned int var, unsigned int qp, OutputType& du) const
{
  // Get local-to-global dof index lookup
  libmesh_assert_greater (this->get_dof_indices().size(), var);
  const unsigned int n_dofs = cast_int<unsigned int>
    (this->get_dof_indices(var).size());

  // Get current local coefficients
  const DenseSubVector<Number> &coef = (this->*subsolution_getter)(var);

  // Get finite element object
  typename FENeeded<OutputType>::grad_base* fe = NULL;
  (this->*fe_getter)( var, fe );

  // Get shape function values at quadrature point
  const std::vector<std::vector
    <typename FENeeded<OutputType>::grad_base::OutputGradient> >
    &dphi = fe->get_dphi();

  // Accumulate solution derivatives
  du = 0;

  for (unsigned int l=0; l != n_dofs; l++)
    du.add_scaled(dphi[l][qp], coef(l));

  return;
}



#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
template<typename OutputType,
         typename FEMContext::FENeeded<OutputType>::hess_getter fe_getter,
         FEMContext::diff_subsolution_getter subsolution_getter>
void FEMContext::some_hessian(unsigned int var, unsigned int qp, OutputType& d2u) const
{
  // Get local-to-global dof index lookup
  libmesh_assert_greater (this->get_dof_indices().size(), var);
  const unsigned int n_dofs = cast_int<unsigned int>
    (this->get_dof_indices(var).size());

  // Get current local coefficients
  const DenseSubVector<Number> &coef = (this->*subsolution_getter)(var);

  // Get finite element object
  typename FENeeded<OutputType>::hess_base* fe = NULL;
  (this->*fe_getter)( var, fe );

  // Get shape function values at quadrature point
  const std::vector<std::vector
    <typename FENeeded<OutputType>::hess_base::OutputTensor> >
    &d2phi = fe->get_d2phi();

  // Accumulate solution second derivatives
  d2u = 0.0;

  for (unsigned int l=0; l != n_dofs; l++)
    d2u.add_scaled(d2phi[l][qp], coef(l));

  return;
}
#endif



Number FEMContext::interior_value(unsigned int var, unsigned int qp) const
{
  Number u;

  this->interior_value( var, qp, u );

  return u;
}

template<typename OutputType>
void FEMContext::interior_value(unsigned int var, unsigned int qp,
                                OutputType& u) const
{
  this->some_value
    <OutputType,
     &FEMContext::get_element_fe
       <typename TensorTools::MakeReal<OutputType>::type>,
     &DiffContext::get_elem_solution>
    (var, qp, u);
}


template<typename OutputType>
void FEMContext::interior_values (unsigned int var,
                                  const NumericVector<Number> & _system_vector,
                                  std::vector<OutputType>& u_vals) const
{
  typedef typename TensorTools::MakeReal<OutputType>::type OutputShape;

  // Get local-to-global dof index lookup
  libmesh_assert_greater (this->get_dof_indices().size(), var);
  const unsigned int n_dofs = cast_int<unsigned int>
    (this->get_dof_indices(var).size());

  // Get current local coefficients
  const DenseSubVector<Number> &coef = get_localized_subvector(_system_vector, var);

  // Get the finite element object
  FEGenericBase<OutputShape>* fe = NULL;
  this->get_element_fe<OutputShape>( var, fe );

  // Get shape function values at quadrature point
  const std::vector<std::vector<OutputShape> > &phi = fe->get_phi();

  // Loop over all the q_points on this element
  for (unsigned int qp=0; qp != u_vals.size(); qp++)
    {
      OutputType &u = u_vals[qp];

      // Compute the value at this q_point
      u = 0.;

      for (unsigned int l=0; l != n_dofs; l++)
        u += phi[l][qp] * coef(l);
    }

  return;
}

Gradient FEMContext::interior_gradient(unsigned int var, unsigned int qp) const
{
  Gradient du;

  this->interior_gradient( var, qp, du );

  return du;
}



template<typename OutputType>
void FEMContext::interior_gradient(unsigned int var, unsigned int qp,
                                   OutputType& du) const
{
  this->some_gradient
    <OutputType,
     &FEMContext::get_element_fe
       <typename TensorTools::MakeReal
         <typename TensorTools::DecrementRank
           <OutputType>::type>::type>,
     &DiffContext::get_elem_solution>
    (var, qp, du);
}



template<typename OutputType>
void FEMContext::interior_gradients
(unsigned int var,
 const NumericVector<Number> & _system_vector,
 std::vector<OutputType>& du_vals) const
{
  typedef typename TensorTools::MakeReal
    <typename TensorTools::DecrementRank<OutputType>::type>::type
    OutputShape;

  // Get local-to-global dof index lookup
  libmesh_assert_greater (this->get_dof_indices().size(), var);
  const unsigned int n_dofs = cast_int<unsigned int>
    (this->get_dof_indices(var).size());

  // Get current local coefficients
  const DenseSubVector<Number> &coef = get_localized_subvector(_system_vector, var);

  // Get finite element object
  FEGenericBase<OutputShape>* fe = NULL;
  this->get_element_fe<OutputShape>( var, fe );

  // Get shape function values at quadrature point
  const std::vector<std::vector<typename FEGenericBase<OutputShape>::OutputGradient> > &dphi = fe->get_dphi();

  // Loop over all the q_points in this finite element
  for (unsigned int qp=0; qp != du_vals.size(); qp++)
    {
      OutputType &du = du_vals[qp];

      // Compute the gradient at this q_point
      du = 0;

      for (unsigned int l=0; l != n_dofs; l++)
        du.add_scaled(dphi[l][qp], coef(l));
    }

  return;
}

#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
Tensor FEMContext::interior_hessian(unsigned int var, unsigned int qp) const
{
  Tensor d2u;

  this->interior_hessian( var, qp, d2u );

  return d2u;
}

template<typename OutputType>
void FEMContext::interior_hessian(unsigned int var, unsigned int qp,
                                  OutputType& d2u) const
{
  this->some_hessian
    <OutputType,
     &FEMContext::get_element_fe
       <typename TensorTools::MakeReal
         <typename TensorTools::DecrementRank
           <typename TensorTools::DecrementRank
             <OutputType>::type>::type>::type>,
     &DiffContext::get_elem_solution>
    (var, qp, d2u);
}


template<typename OutputType>
void FEMContext::interior_hessians
(unsigned int var,
 const NumericVector<Number> & _system_vector,
 std::vector<OutputType>& d2u_vals) const
{
  typedef typename TensorTools::MakeReal<
    typename TensorTools::DecrementRank<
    typename TensorTools::DecrementRank<
    OutputType>::type>::type>::type
    OutputShape;

  // Get local-to-global dof index lookup
  libmesh_assert_greater (this->get_dof_indices().size(), var);
  const unsigned int n_dofs = cast_int<unsigned int>
    (this->get_dof_indices(var).size());

  // Get current local coefficients
  const DenseSubVector<Number> &coef = get_localized_subvector(_system_vector, var);

  // Get finite element object
  FEGenericBase<OutputShape>* fe = NULL;
  this->get_element_fe<OutputShape>( var, fe );

  // Get shape function values at quadrature point
  const std::vector<std::vector<typename FEGenericBase<OutputShape>::OutputTensor> > &d2phi = fe->get_d2phi();

  // Loop over all the q_points in this finite element
  for (unsigned int qp=0; qp != d2u_vals.size(); qp++)
    {
      OutputType &d2u = d2u_vals[qp];

      // Compute the gradient at this q_point
      d2u = 0;

      for (unsigned int l=0; l != n_dofs; l++)
        d2u.add_scaled(d2phi[l][qp], coef(l));
    }

  return;
}


#endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES


template<typename OutputType>
void FEMContext::interior_curl(unsigned int var, unsigned int qp,
                               OutputType& curl_u) const
{
  typedef typename TensorTools::MakeReal<OutputType>::type OutputShape;

  // Get local-to-global dof index lookup
  libmesh_assert_greater (this->get_dof_indices().size(), var);
  const unsigned int n_dofs = cast_int<unsigned int>
    (this->get_dof_indices(var).size());

  // Get current local coefficients
  libmesh_assert_greater (this->_elem_subsolutions.size(), var);
  libmesh_assert(&(this->get_elem_solution(var)));
  const DenseSubVector<Number> &coef = this->get_elem_solution(var);

  // Get finite element object
  FEGenericBase<OutputShape>* fe = NULL;
  this->get_element_fe<OutputShape>( var, fe );

  // Get shape function values at quadrature point
  const std::vector<std::vector<typename FEGenericBase<OutputShape>::OutputShape> > &curl_phi = fe->get_curl_phi();

  // Accumulate solution curl
  curl_u = 0.;

  for (unsigned int l=0; l != n_dofs; l++)
    curl_u.add_scaled(curl_phi[l][qp], coef(l));

  return;
}


template<typename OutputType>
void FEMContext::interior_div(unsigned int var, unsigned int qp,
                              OutputType& div_u) const
{
  typedef typename
    TensorTools::IncrementRank
    <typename TensorTools::MakeReal<OutputType>::type>::type OutputShape;

  // Get local-to-global dof index lookup
  libmesh_assert_greater (this->get_dof_indices().size(), var);
  const unsigned int n_dofs = cast_int<unsigned int>
    (this->get_dof_indices(var).size());

  // Get current local coefficients
  libmesh_assert_greater (this->_elem_subsolutions.size(), var);
  libmesh_assert(&(this->get_elem_solution(var)));
  const DenseSubVector<Number> &coef = this->get_elem_solution(var);

  // Get finite element object
  FEGenericBase<OutputShape>* fe = NULL;
  this->get_element_fe<OutputShape>( var, fe );

  // Get shape function values at quadrature point
  const std::vector<std::vector<typename FEGenericBase<OutputShape>::OutputDivergence> > &div_phi = fe->get_div_phi();

  // Accumulate solution curl
  div_u = 0.;

  for (unsigned int l=0; l != n_dofs; l++)
    div_u += div_phi[l][qp] * coef(l);

  return;
}


Number FEMContext::side_value(unsigned int var, unsigned int qp) const
{
  Number u = 0.;

  this->side_value( var, qp, u );

  return u;
}


template<typename OutputType>
void FEMContext::side_value(unsigned int var, unsigned int qp,
                            OutputType& u) const
{
  this->some_value
    <OutputType,
     &FEMContext::get_side_fe
       <typename TensorTools::MakeReal<OutputType>::type>,
     &DiffContext::get_elem_solution>
    (var, qp, u);
}


template<typename OutputType>
void FEMContext::side_values
(unsigned int var,
 const NumericVector<Number> & _system_vector,
 std::vector<OutputType>& u_vals) const
{
  typedef typename TensorTools::MakeReal<OutputType>::type OutputShape;

  // Get local-to-global dof index lookup
  libmesh_assert_greater (this->get_dof_indices().size(), var);
  const unsigned int n_dofs = cast_int<unsigned int>
    (this->get_dof_indices(var).size());

  // Get current local coefficients
  const DenseSubVector<Number> &coef = get_localized_subvector(_system_vector, var);

  // Get the finite element object
  FEGenericBase<OutputShape>* the_side_fe = NULL;
  this->get_side_fe<OutputShape>( var, the_side_fe );

  // Get shape function values at quadrature point
  const std::vector<std::vector<OutputShape> > &phi = the_side_fe->get_phi();

  // Loop over all the q_points on this element
  for (unsigned int qp=0; qp != u_vals.size(); qp++)
    {
      OutputType &u = u_vals[qp];

      // Compute the value at this q_point
      u = 0.;

      for (unsigned int l=0; l != n_dofs; l++)
        u += phi[l][qp] * coef(l);
    }

  return;
}

Gradient FEMContext::side_gradient(unsigned int var, unsigned int qp) const
{
  Gradient du;

  this->side_gradient( var, qp, du );

  return du;
}


template<typename OutputType>
void FEMContext::side_gradient(unsigned int var, unsigned int qp,
                               OutputType& du) const
{
  typedef typename TensorTools::MakeReal
    <typename TensorTools::DecrementRank<OutputType>::type>::type
    OutputShape;

  // Get local-to-global dof index lookup
  libmesh_assert_greater (this->get_dof_indices().size(), var);
  const unsigned int n_dofs = cast_int<unsigned int>
    (this->get_dof_indices(var).size());

  // Get current local coefficients
  libmesh_assert_greater (this->_elem_subsolutions.size(), var);
  libmesh_assert(&(this->get_elem_solution(var)));
  const DenseSubVector<Number> &coef = this->get_elem_solution(var);

  // Get finite element object
  FEGenericBase<OutputShape>* the_side_fe = NULL;
  this->get_side_fe<OutputShape>( var, the_side_fe );

  // Get shape function values at quadrature point
  const std::vector<std::vector< typename FEGenericBase<OutputShape>::OutputGradient> > &dphi = the_side_fe->get_dphi();

  // Accumulate solution derivatives
  du = 0.;

  for (unsigned int l=0; l != n_dofs; l++)
    du.add_scaled(dphi[l][qp], coef(l));

  return;
}



template<typename OutputType>
void FEMContext::side_gradients
(unsigned int var,
 const NumericVector<Number> & _system_vector,
 std::vector<OutputType>& du_vals) const
{
  typedef typename TensorTools::MakeReal
    <typename TensorTools::DecrementRank<OutputType>::type>::type
    OutputShape;

  // Get local-to-global dof index lookup
  libmesh_assert_greater (this->get_dof_indices().size(), var);
  const unsigned int n_dofs = cast_int<unsigned int>
    (this->get_dof_indices(var).size());

  // Get current local coefficients
  const DenseSubVector<Number> &coef = get_localized_subvector(_system_vector, var);

  // Get finite element object
  FEGenericBase<OutputShape>* the_side_fe = NULL;
  this->get_side_fe<OutputShape>( var, the_side_fe );

  // Get shape function values at quadrature point
  const std::vector<std::vector<typename FEGenericBase<OutputShape>::OutputGradient> > &dphi = the_side_fe->get_dphi();

  // Loop over all the q_points in this finite element
  for (unsigned int qp=0; qp != du_vals.size(); qp++)
    {
      OutputType &du = du_vals[qp];

      du = 0;

      // Compute the gradient at this q_point
      for (unsigned int l=0; l != n_dofs; l++)
        du.add_scaled(dphi[l][qp], coef(l));
    }

  return;
}

#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
Tensor FEMContext::side_hessian(unsigned int var, unsigned int qp) const
{
  Tensor d2u;

  this->side_hessian( var, qp, d2u );

  return d2u;
}



template<typename OutputType>
void FEMContext::side_hessian(unsigned int var, unsigned int qp,
                              OutputType& d2u) const
{
  this->some_hessian
    <OutputType,
     &FEMContext::get_side_fe
       <typename TensorTools::MakeReal
         <typename TensorTools::DecrementRank
           <typename TensorTools::DecrementRank
             <OutputType>::type>::type>::type>,
     &DiffContext::get_elem_solution>
    (var, qp, d2u);
}



template<typename OutputType>
void FEMContext::side_hessians
(unsigned int var,
 const NumericVector<Number> & _system_vector,
 std::vector<OutputType>& d2u_vals) const
{
  typedef typename TensorTools::MakeReal<
    typename TensorTools::DecrementRank<
    typename TensorTools::DecrementRank<
    OutputType>::type>::type>::type
    OutputShape;

  // Get local-to-global dof index lookup
  libmesh_assert_greater (this->get_dof_indices().size(), var);
  const unsigned int n_dofs = cast_int<unsigned int>
    (this->get_dof_indices(var).size());

  // Get current local coefficients
  const DenseSubVector<Number> &coef = get_localized_subvector(_system_vector, var);

  // Get finite element object
  FEGenericBase<OutputShape>* the_side_fe = NULL;
  this->get_side_fe<OutputShape>( var, the_side_fe );

  // Get shape function values at quadrature point
  const std::vector<std::vector<typename FEGenericBase<OutputShape>::OutputTensor> > &d2phi = the_side_fe->get_d2phi();

  // Loop over all the q_points in this finite element
  for (unsigned int qp=0; qp != d2u_vals.size(); qp++)
    {
      OutputType &d2u = d2u_vals[qp];

      // Compute the gradient at this q_point
      d2u = 0;

      for (unsigned int l=0; l != n_dofs; l++)
        d2u.add_scaled(d2phi[l][qp], coef(l));
    }

  return;
}



#endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES



Number FEMContext::point_value(unsigned int var, const Point &p) const
{
  Number u = 0.;

  this->point_value( var, p, u );

  return u;
}

template<typename OutputType>
void FEMContext::point_value(unsigned int var, const Point &p,
                             OutputType& u) const
{
  typedef typename TensorTools::MakeReal<OutputType>::type OutputShape;

  // Get local-to-global dof index lookup
  libmesh_assert_greater (this->get_dof_indices().size(), var);
  const unsigned int n_dofs = cast_int<unsigned int>
    (this->get_dof_indices(var).size());

  // Get current local coefficients
  libmesh_assert_greater (this->_elem_subsolutions.size(), var);
  libmesh_assert(&(this->get_elem_solution(var)));
  const DenseSubVector<Number> &coef = this->get_elem_solution(var);

  // Get finite element object
  FEGenericBase<OutputShape>* fe = NULL;
  this->get_element_fe<OutputShape>( var, fe );

  // Build a FE for calculating u(p)
  UniquePtr<FEGenericBase<OutputShape> > fe_new = this->build_new_fe( fe, p );

  // Get the values of the shape function derivatives
  const std::vector<std::vector<OutputShape> >&  phi = fe_new->get_phi();

  u = 0.;

  for (unsigned int l=0; l != n_dofs; l++)
    u += phi[l][0] * coef(l);

  return;
}



Gradient FEMContext::point_gradient(unsigned int var, const Point &p) const
{
  Gradient grad_u;

  this->point_gradient( var, p, grad_u );

  return grad_u;
}



template<typename OutputType>
void FEMContext::point_gradient(unsigned int var, const Point &p,
                                OutputType& grad_u) const
{
  typedef typename TensorTools::MakeReal
    <typename TensorTools::DecrementRank<OutputType>::type>::type
    OutputShape;

  // Get local-to-global dof index lookup
  libmesh_assert_greater (this->get_dof_indices().size(), var);
  const unsigned int n_dofs = cast_int<unsigned int>
    (this->get_dof_indices(var).size());

  // Get current local coefficients
  libmesh_assert_greater (this->_elem_subsolutions.size(), var);
  libmesh_assert(&(this->get_elem_solution(var)));
  const DenseSubVector<Number> &coef = this->get_elem_solution(var);

  // Get finite element object
  FEGenericBase<OutputShape>* fe = NULL;
  this->get_element_fe<OutputShape>( var, fe );

  // Build a FE for calculating u(p)
  UniquePtr<FEGenericBase<OutputShape> > fe_new = this->build_new_fe( fe, p );

  // Get the values of the shape function derivatives
  const std::vector<std::vector<typename FEGenericBase<OutputShape>::OutputGradient> >&  dphi = fe_new->get_dphi();

  grad_u = 0.0;

  for (unsigned int l=0; l != n_dofs; l++)
    grad_u.add_scaled(dphi[l][0], coef(l));

  return;
}



#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES

Tensor FEMContext::point_hessian(unsigned int var, const Point &p) const
{
  Tensor hess_u;

  this->point_hessian( var, p, hess_u );

  return hess_u;
}


template<typename OutputType>
void FEMContext::point_hessian(unsigned int var, const Point &p,
                               OutputType& hess_u) const
{
  typedef typename TensorTools::MakeReal<
    typename TensorTools::DecrementRank<
    typename TensorTools::DecrementRank<
    OutputType>::type>::type>::type
    OutputShape;

  // Get local-to-global dof index lookup
  libmesh_assert_greater (this->get_dof_indices().size(), var);
  const unsigned int n_dofs = cast_int<unsigned int>
    (this->get_dof_indices(var).size());

  // Get current local coefficients
  libmesh_assert_greater (this->_elem_subsolutions.size(), var);
  libmesh_assert(&(this->get_elem_solution(var)));
  const DenseSubVector<Number> &coef = this->get_elem_solution(var);

  // Get finite element object
  FEGenericBase<OutputShape>* fe = NULL;
  this->get_element_fe<OutputShape>( var, fe );

  // Build a FE for calculating u(p)
  UniquePtr<FEGenericBase<OutputShape> > fe_new = this->build_new_fe( fe, p );

  // Get the values of the shape function derivatives
  const std::vector<std::vector<typename FEGenericBase<OutputShape>::OutputTensor> >&  d2phi = fe_new->get_d2phi();

  hess_u = 0.0;

  for (unsigned int l=0; l != n_dofs; l++)
    hess_u.add_scaled(d2phi[l][0], coef(l));

  return;
}

#endif // LIBMESH_ENABLE_SECOND_DERIVATIVES


template<typename OutputType>
void FEMContext::point_curl(unsigned int var, const Point &p,
                            OutputType& curl_u) const
{
  typedef typename TensorTools::MakeReal<OutputType>::type OutputShape;

  // Get local-to-global dof index lookup
  libmesh_assert_greater (this->get_dof_indices().size(), var);
  const unsigned int n_dofs = cast_int<unsigned int>
    (this->get_dof_indices(var).size());

  // Get current local coefficients
  libmesh_assert_greater (this->_elem_subsolutions.size(), var);
  libmesh_assert(&(this->get_elem_solution(var)));
  const DenseSubVector<Number> &coef = this->get_elem_solution(var);

  // Get finite element object
  FEGenericBase<OutputShape>* fe = NULL;
  this->get_element_fe<OutputShape>( var, fe );

  // Build a FE for calculating u(p)
  UniquePtr<FEGenericBase<OutputShape> > fe_new = this->build_new_fe( fe, p );

  // Get the values of the shape function derivatives
  const std::vector<std::vector<typename FEGenericBase<OutputShape>::OutputShape> >&  curl_phi = fe_new->get_curl_phi();

  curl_u = 0.0;

  for (unsigned int l=0; l != n_dofs; l++)
    curl_u.add_scaled(curl_phi[l][0], coef(l));

  return;
}



Number FEMContext::fixed_interior_value(unsigned int var, unsigned int qp) const
{
  Number u = 0.;

  this->fixed_interior_value( var, qp, u );

  return u;
}



template<typename OutputType>
void FEMContext::fixed_interior_value(unsigned int var, unsigned int qp,
                                      OutputType& u) const
{
  this->some_value
    <OutputType,
     &FEMContext::get_element_fe
       <typename TensorTools::MakeReal<OutputType>::type>,
     &DiffContext::get_elem_fixed_solution>
    (var, qp, u);
}



Gradient FEMContext::fixed_interior_gradient(unsigned int var, unsigned int qp) const
{
  Gradient du;

  this->fixed_interior_gradient( var, qp, du );

  return du;
}


template<typename OutputType>
void FEMContext::FEMContext::fixed_interior_gradient(unsigned int var, unsigned int qp,
                                                     OutputType& du) const
{
  this->some_gradient
    <OutputType,
     &FEMContext::get_element_fe
       <typename TensorTools::MakeReal
         <typename TensorTools::DecrementRank
           <OutputType>::type>::type>,
     &DiffContext::get_elem_fixed_solution>
    (var, qp, du);
}



#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
Tensor FEMContext::fixed_interior_hessian(unsigned int var, unsigned int qp) const
{
  Tensor d2u;

  this->fixed_interior_hessian( var, qp, d2u );

  return d2u;
}


template<typename OutputType>
void FEMContext::fixed_interior_hessian(unsigned int var, unsigned int qp,
                                        OutputType& d2u) const
{
  this->some_hessian
    <OutputType,
     &FEMContext::get_element_fe
       <typename TensorTools::MakeReal
         <typename TensorTools::DecrementRank
           <typename TensorTools::DecrementRank
             <OutputType>::type>::type>::type>,
     &DiffContext::get_elem_fixed_solution>
    (var, qp, d2u);
}
#endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES



Number FEMContext::fixed_side_value(unsigned int var, unsigned int qp) const
{
  Number u = 0.;

  this->fixed_side_value( var, qp, u );

  return u;
}


template<typename OutputType>
void FEMContext::fixed_side_value(unsigned int var, unsigned int qp,
                                  OutputType& u) const
{
  this->some_value
    <OutputType,
     &FEMContext::get_side_fe
       <typename TensorTools::MakeReal<OutputType>::type>,
     &DiffContext::get_elem_fixed_solution>
    (var, qp, u);
}



Gradient FEMContext::fixed_side_gradient(unsigned int var, unsigned int qp) const
{
  Gradient du;

  this->fixed_side_gradient( var, qp, du );

  return du;
}


template<typename OutputType>
void FEMContext::FEMContext::fixed_side_gradient(unsigned int var, unsigned int qp,
                                                 OutputType& du) const
{
  this->some_gradient
    <OutputType,
     &FEMContext::get_side_fe
       <typename TensorTools::MakeReal
         <typename TensorTools::DecrementRank
           <OutputType>::type>::type>,
     &DiffContext::get_elem_fixed_solution>
    (var, qp, du);
}



#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
Tensor FEMContext::fixed_side_hessian(unsigned int var, unsigned int qp) const
{
  Tensor d2u;

  this->fixed_side_hessian( var, qp, d2u );

  return d2u;
}

template<typename OutputType>
void FEMContext::fixed_side_hessian(unsigned int var, unsigned int qp,
                                    OutputType& d2u) const
{
  this->some_hessian
    <OutputType,
     &FEMContext::get_side_fe
       <typename TensorTools::MakeReal
         <typename TensorTools::DecrementRank
           <typename TensorTools::DecrementRank
             <OutputType>::type>::type>::type>,
     &DiffContext::get_elem_fixed_solution>
    (var, qp, d2u);
}
#endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES



Number FEMContext::fixed_point_value(unsigned int var, const Point &p) const
{
  Number u = 0.;

  this->fixed_point_value( var, p, u );

  return u;
}

template<typename OutputType>
void FEMContext::fixed_point_value(unsigned int var, const Point &p,
                                   OutputType& u) const
{
  typedef typename TensorTools::MakeReal<OutputType>::type OutputShape;

  // Get local-to-global dof index lookup
  libmesh_assert_greater (this->get_dof_indices().size(), var);
  const unsigned int n_dofs = cast_int<unsigned int>
    (this->get_dof_indices(var).size());

  // Get current local coefficients
  libmesh_assert_greater (_elem_fixed_subsolutions.size(), var);
  libmesh_assert(&(this->get_elem_fixed_solution(var)));
  const DenseSubVector<Number> &coef = this->get_elem_fixed_solution(var);

  // Get finite element object
  FEGenericBase<OutputShape>* fe = NULL;
  this->get_element_fe<OutputShape>( var, fe );

  // Build a FE for calculating u(p)
  UniquePtr<FEGenericBase<OutputShape> > fe_new = this->build_new_fe( fe, p );

  // Get the values of the shape function derivatives
  const std::vector<std::vector<OutputShape> >&  phi = fe_new->get_phi();

  u = 0.;

  for (unsigned int l=0; l != n_dofs; l++)
    u += phi[l][0] * coef(l);

  return;
}



Gradient FEMContext::fixed_point_gradient(unsigned int var, const Point &p) const
{
  Gradient grad_u;

  this->fixed_point_gradient( var, p, grad_u );

  return grad_u;
}



template<typename OutputType>
void FEMContext::fixed_point_gradient(unsigned int var, const Point &p,
                                      OutputType& grad_u) const
{
  typedef typename TensorTools::MakeReal
    <typename TensorTools::DecrementRank<OutputType>::type>::type
    OutputShape;

  // Get local-to-global dof index lookup
  libmesh_assert_greater (this->get_dof_indices().size(), var);
  const unsigned int n_dofs = cast_int<unsigned int>
    (this->get_dof_indices(var).size());

  // Get current local coefficients
  libmesh_assert_greater (_elem_fixed_subsolutions.size(), var);
  libmesh_assert(&(this->get_elem_fixed_solution(var)));
  const DenseSubVector<Number> &coef = this->get_elem_fixed_solution(var);

  // Get finite element object
  FEGenericBase<OutputShape>* fe = NULL;
  this->get_element_fe<OutputShape>( var, fe );

  // Build a FE for calculating u(p)
  UniquePtr<FEGenericBase<OutputShape> > fe_new = this->build_new_fe( fe, p );

  // Get the values of the shape function derivatives
  const std::vector<std::vector<typename FEGenericBase<OutputShape>::OutputGradient> >&  dphi = fe_new->get_dphi();

  grad_u = 0.0;

  for (unsigned int l=0; l != n_dofs; l++)
    grad_u.add_scaled(dphi[l][0], coef(l));

  return;
}


#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES

Tensor FEMContext::fixed_point_hessian(unsigned int var, const Point &p) const
{
  Tensor hess_u;

  this->fixed_point_hessian( var, p, hess_u );

  return hess_u;
}



template<typename OutputType>
void FEMContext::fixed_point_hessian(unsigned int var, const Point &p,
                                     OutputType& hess_u) const
{
  typedef typename TensorTools::MakeReal<
    typename TensorTools::DecrementRank<
    typename TensorTools::DecrementRank<
    OutputType>::type>::type>::type
    OutputShape;

  // Get local-to-global dof index lookup
  libmesh_assert_greater (this->get_dof_indices().size(), var);
  const unsigned int n_dofs = cast_int<unsigned int>
    (this->get_dof_indices(var).size());

  // Get current local coefficients
  libmesh_assert_greater (_elem_fixed_subsolutions.size(), var);
  libmesh_assert(&(this->get_elem_fixed_solution(var)));
  const DenseSubVector<Number> &coef = this->get_elem_fixed_solution(var);<commit_msg>Silence clang 3.5.1 warnings.<commit_after>// The libMesh Finite Element Library.
// Copyright (C) 2002-2014 <NAME>, <NAME>, <NAME>

// This library is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.

// This library is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// Lesser General Public License for more details.

// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA



#include "libmesh/boundary_info.h"
#include "libmesh/dof_map.h"
#include "libmesh/elem.h"
#include "libmesh/fe_base.h"
#include "libmesh/fe_interface.h"
#include "libmesh/fem_context.h"
#include "libmesh/libmesh_logging.h"
#include "libmesh/mesh_base.h"
#include "libmesh/quadrature.h"
#include "libmesh/system.h"
#include "libmesh/time_solver.h"
#include "libmesh/unsteady_solver.h" // For euler_residual

namespace libMesh
{

FEMContext::FEMContext (const System &sys)
  : DiffContext(sys),
    _mesh_sys(NULL),
    _mesh_x_var(0),
    _mesh_y_var(0),
    _mesh_z_var(0),
    side(0), edge(0),
    _boundary_info(sys.get_mesh().get_boundary_info()),
    _elem(NULL),
    _dim(sys.get_mesh().mesh_dimension()),
    _elem_dim(0), /* This will be reset in set_elem(). */
    _elem_dims(sys.get_mesh().elem_dimensions()),
    _element_qrule(4,NULL),
    _side_qrule(4,NULL),
    _edge_qrule(NULL)
{
  // Reserve space for the FEAbstract and QBase objects for each
  // element dimension possiblity (0,1,2,3)
  _element_fe.resize(4);
  _side_fe.resize(4);
  _element_fe_var.resize(4);
  _side_fe_var.resize(4);

  // We need to know which of our variables has the hardest
  // shape functions to numerically integrate.

  unsigned int nv = sys.n_vars();

  libmesh_assert (nv);
  FEType hardest_fe_type = sys.variable_type(0);

  bool have_scalar = false;

  for (unsigned int i=0; i != nv; ++i)
    {
      FEType fe_type = sys.variable_type(i);

      // FIXME - we don't yet handle mixed finite elements from
      // different families which require different quadrature rules
      // libmesh_assert_equal_to (fe_type.family, hardest_fe_type.family);

      if (fe_type.order > hardest_fe_type.order)
        hardest_fe_type = fe_type;

      // We need to detect SCALAR's so we can prepare FE objects for them.
      if( fe_type.family == SCALAR )
        have_scalar = true;
    }

  if(have_scalar)
    // SCALAR FEs have dimension 0 by assumption
    _elem_dims.insert(0);

  for( std::set<unsigned char>::const_iterator dim_it = _elem_dims.begin();
       dim_it != _elem_dims.end(); ++dim_it )
    {
      const unsigned char dim = *dim_it;

      // Create an adequate quadrature rule
      _element_qrule[dim] = hardest_fe_type.default_quadrature_rule
        (dim, sys.extra_quadrature_order).release();
      _side_qrule[dim] = hardest_fe_type.default_quadrature_rule
        (dim-1, sys.extra_quadrature_order).release();
      if (dim == 3)
        _edge_qrule = hardest_fe_type.default_quadrature_rule
          (1, sys.extra_quadrature_order).release();

      // Next, create finite element objects
      _element_fe_var[dim].resize(nv);
      _side_fe_var[dim].resize(nv);
      if (dim == 3)
        _edge_fe_var.resize(nv);


      for (unsigned int i=0; i != nv; ++i)
        {
          FEType fe_type = sys.variable_type(i);

          if ( _element_fe[dim][fe_type] == NULL )
            {
              _element_fe[dim][fe_type] = FEAbstract::build(dim, fe_type).release();
              _element_fe[dim][fe_type]->attach_quadrature_rule(_element_qrule[dim]);
              _side_fe[dim][fe_type] = FEAbstract::build(dim, fe_type).release();
              _side_fe[dim][fe_type]->attach_quadrature_rule(_side_qrule[dim]);

              if (this->_dim == 3)
                {
                  _edge_fe[fe_type] = FEAbstract::build(dim, fe_type).release();
                  _edge_fe[fe_type]->attach_quadrature_rule(_edge_qrule);
                }
            }

          _element_fe_var[dim][i] = _element_fe[dim][fe_type];
          _side_fe_var[dim][i] = _side_fe[dim][fe_type];
          if ((dim) == 3)
            _edge_fe_var[i] = _edge_fe[fe_type];

        }
    }
}

FEMContext::~FEMContext()
{
  // We don't want to store UniquePtrs in STL containers, but we don't
  // want to leak memory either
  for (std::vector<std::map<FEType, FEAbstract *> >::iterator d = _element_fe.begin();
       d != _element_fe.end(); ++d)
    for (std::map<FEType, FEAbstract *>::iterator i = d->begin();
         i != d->end(); ++i)
      delete i->second;

  for (std::vector<std::map<FEType, FEAbstract *> >::iterator d = _side_fe.begin();
       d != _side_fe.end(); ++d)
    for (std::map<FEType, FEAbstract *>::iterator i = d->begin();
         i != d->end(); ++i)
      delete i->second;

  for (std::map<FEType, FEAbstract *>::iterator i = _edge_fe.begin();
       i != _edge_fe.end(); ++i)
    delete i->second;
  _edge_fe.clear();

  for (std::vector<QBase*>::iterator i = _element_qrule.begin();
       i != _element_qrule.end(); ++i)
    delete *i;
  _element_qrule.clear();

  for (std::vector<QBase*>::iterator i = _side_qrule.begin();
       i != _side_qrule.end(); ++i)
    delete *i;
  _side_qrule.clear();

  delete _edge_qrule;
  _edge_qrule = NULL;
}



bool FEMContext::has_side_boundary_id(boundary_id_type id) const
{
  return _boundary_info.has_boundary_id(&(this->get_elem()), side, id);
}


std::vector<boundary_id_type> FEMContext::side_boundary_ids() const
{
  return _boundary_info.boundary_ids(&(this->get_elem()), side);
}



template<typename OutputType,
         typename FEMContext::FENeeded<OutputType>::value_getter fe_getter,
         FEMContext::diff_subsolution_getter subsolution_getter>
void FEMContext::some_value(unsigned int var, unsigned int qp, OutputType& u) const
{
  // Get local-to-global dof index lookup
  libmesh_assert_greater (this->get_dof_indices().size(), var);
  const unsigned int n_dofs = cast_int<unsigned int>
    (this->get_dof_indices(var).size());

  // Get current local coefficients
  const DenseSubVector<Number> &coef = (this->*subsolution_getter)(var);

  // Get finite element object
  typename FENeeded<OutputType>::value_base* fe = NULL;
  (this->*fe_getter)( var, fe );

  // Get shape function values at quadrature point
  const std::vector<std::vector
    <typename FENeeded<OutputType>::value_shape> > &phi = fe->get_phi();

  // Accumulate solution value
  u = 0.;

  for (unsigned int l=0; l != n_dofs; l++)
    u += phi[l][qp] * coef(l);
}



template<typename OutputType,
         typename FEMContext::FENeeded<OutputType>::grad_getter fe_getter,
         FEMContext::diff_subsolution_getter subsolution_getter>
void FEMContext::some_gradient(unsigned int var, unsigned int qp, OutputType& du) const
{
  // Get local-to-global dof index lookup
  libmesh_assert_greater (this->get_dof_indices().size(), var);
  const unsigned int n_dofs = cast_int<unsigned int>
    (this->get_dof_indices(var).size());

  // Get current local coefficients
  const DenseSubVector<Number> &coef = (this->*subsolution_getter)(var);

  // Get finite element object
  typename FENeeded<OutputType>::grad_base* fe = NULL;
  (this->*fe_getter)( var, fe );

  // Get shape function values at quadrature point
  const std::vector<std::vector
    <typename FENeeded<OutputType>::grad_base::OutputGradient> >
    &dphi = fe->get_dphi();

  // Accumulate solution derivatives
  du = 0;

  for (unsigned int l=0; l != n_dofs; l++)
    du.add_scaled(dphi[l][qp], coef(l));

  return;
}



#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
template<typename OutputType,
         typename FEMContext::FENeeded<OutputType>::hess_getter fe_getter,
         FEMContext::diff_subsolution_getter subsolution_getter>
void FEMContext::some_hessian(unsigned int var, unsigned int qp, OutputType& d2u) const
{
  // Get local-to-global dof index lookup
  libmesh_assert_greater (this->get_dof_indices().size(), var);
  const unsigned int n_dofs = cast_int<unsigned int>
    (this->get_dof_indices(var).size());

  // Get current local coefficients
  const DenseSubVector<Number> &coef = (this->*subsolution_getter)(var);

  // Get finite element object
  typename FENeeded<OutputType>::hess_base* fe = NULL;
  (this->*fe_getter)( var, fe );

  // Get shape function values at quadrature point
  const std::vector<std::vector
    <typename FENeeded<OutputType>::hess_base::OutputTensor> >
    &d2phi = fe->get_d2phi();

  // Accumulate solution second derivatives
  d2u = 0.0;

  for (unsigned int l=0; l != n_dofs; l++)
    d2u.add_scaled(d2phi[l][qp], coef(l));

  return;
}
#endif



Number FEMContext::interior_value(unsigned int var, unsigned int qp) const
{
  Number u;

  this->interior_value( var, qp, u );

  return u;
}

template<typename OutputType>
void FEMContext::interior_value(unsigned int var, unsigned int qp,
                                OutputType& u) const
{
  this->some_value
    <OutputType,
     &FEMContext::get_element_fe
       <typename TensorTools::MakeReal<OutputType>::type>,
     &DiffContext::get_elem_solution>
    (var, qp, u);
}


template<typename OutputType>
void FEMContext::interior_values (unsigned int var,
                                  const NumericVector<Number> & _system_vector,
                                  std::vector<OutputType>& u_vals) const
{
  typedef typename TensorTools::MakeReal<OutputType>::type OutputShape;

  // Get local-to-global dof index lookup
  libmesh_assert_greater (this->get_dof_indices().size(), var);
  const unsigned int n_dofs = cast_int<unsigned int>
    (this->get_dof_indices(var).size());

  // Get current local coefficients
  const DenseSubVector<Number> &coef = get_localized_subvector(_system_vector, var);

  // Get the finite element object
  FEGenericBase<OutputShape>* fe = NULL;
  this->get_element_fe<OutputShape>( var, fe );

  // Get shape function values at quadrature point
  const std::vector<std::vector<OutputShape> > &phi = fe->get_phi();

  // Loop over all the q_points on this element
  for (unsigned int qp=0; qp != u_vals.size(); qp++)
    {
      OutputType &u = u_vals[qp];

      // Compute the value at this q_point
      u = 0.;

      for (unsigned int l=0; l != n_dofs; l++)
        u += phi[l][qp] * coef(l);
    }

  return;
}

Gradient FEMContext::interior_gradient(unsigned int var, unsigned int qp) const
{
  Gradient du;

  this->interior_gradient( var, qp, du );

  return du;
}



template<typename OutputType>
void FEMContext::interior_gradient(unsigned int var, unsigned int qp,
                                   OutputType& du) const
{
  this->some_gradient
    <OutputType,
     &FEMContext::get_element_fe
       <typename TensorTools::MakeReal
         <typename TensorTools::DecrementRank
           <OutputType>::type>::type>,
     &DiffContext::get_elem_solution>
    (var, qp, du);
}



template<typename OutputType>
void FEMContext::interior_gradients
(unsigned int var,
 const NumericVector<Number> & _system_vector,
 std::vector<OutputType>& du_vals) const
{
  typedef typename TensorTools::MakeReal
    <typename TensorTools::DecrementRank<OutputType>::type>::type
    OutputShape;

  // Get local-to-global dof index lookup
  libmesh_assert_greater (this->get_dof_indices().size(), var);
  const unsigned int n_dofs = cast_int<unsigned int>
    (this->get_dof_indices(var).size());

  // Get current local coefficients
  const DenseSubVector<Number> &coef = get_localized_subvector(_system_vector, var);

  // Get finite element object
  FEGenericBase<OutputShape>* fe = NULL;
  this->get_element_fe<OutputShape>( var, fe );

  // Get shape function values at quadrature point
  const std::vector<std::vector<typename FEGenericBase<OutputShape>::OutputGradient> > &dphi = fe->get_dphi();

  // Loop over all the q_points in this finite element
  for (unsigned int qp=0; qp != du_vals.size(); qp++)
    {
      OutputType &du = du_vals[qp];

      // Compute the gradient at this q_point
      du = 0;

      for (unsigned int l=0; l != n_dofs; l++)
        du.add_scaled(dphi[l][qp], coef(l));
    }

  return;
}

#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
Tensor FEMContext::interior_hessian(unsigned int var, unsigned int qp) const
{
  Tensor d2u;

  this->interior_hessian( var, qp, d2u );

  return d2u;
}

template<typename OutputType>
void FEMContext::interior_hessian(unsigned int var, unsigned int qp,
                                  OutputType& d2u) const
{
  this->some_hessian
    <OutputType,
     &FEMContext::get_element_fe
       <typename TensorTools::MakeReal
         <typename TensorTools::DecrementRank
           <typename TensorTools::DecrementRank
             <OutputType>::type>::type>::type>,
     &DiffContext::get_elem_solution>
    (var, qp, d2u);
}


template<typename OutputType>
void FEMContext::interior_hessians
(unsigned int var,
 const NumericVector<Number> & _system_vector,
 std::vector<OutputType>& d2u_vals) const
{
  typedef typename TensorTools::MakeReal<
    typename TensorTools::DecrementRank<
    typename TensorTools::DecrementRank<
    OutputType>::type>::type>::type
    OutputShape;

  // Get local-to-global dof index lookup
  libmesh_assert_greater (this->get_dof_indices().size(), var);
  const unsigned int n_dofs = cast_int<unsigned int>
    (this->get_dof_indices(var).size());

  // Get current local coefficients
  const DenseSubVector<Number> &coef = get_localized_subvector(_system_vector, var);

  // Get finite element object
  FEGenericBase<OutputShape>* fe = NULL;
  this->get_element_fe<OutputShape>( var, fe );

  // Get shape function values at quadrature point
  const std::vector<std::vector<typename FEGenericBase<OutputShape>::OutputTensor> > &d2phi = fe->get_d2phi();

  // Loop over all the q_points in this finite element
  for (unsigned int qp=0; qp != d2u_vals.size(); qp++)
    {
      OutputType &d2u = d2u_vals[qp];

      // Compute the gradient at this q_point
      d2u = 0;

      for (unsigned int l=0; l != n_dofs; l++)
        d2u.add_scaled(d2phi[l][qp], coef(l));
    }

  return;
}


#endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES


template<typename OutputType>
void FEMContext::interior_curl(unsigned int var, unsigned int qp,
                               OutputType& curl_u) const
{
  typedef typename TensorTools::MakeReal<OutputType>::type OutputShape;

  // Get local-to-global dof index lookup
  libmesh_assert_greater (this->get_dof_indices().size(), var);
  const unsigned int n_dofs = cast_int<unsigned int>
    (this->get_dof_indices(var).size());

  // Get current local coefficients
  libmesh_assert_greater (this->_elem_subsolutions.size(), var);
  const DenseSubVector<Number> &coef = this->get_elem_solution(var);

  // Get finite element object
  FEGenericBase<OutputShape>* fe = NULL;
  this->get_element_fe<OutputShape>( var, fe );

  // Get shape function values at quadrature point
  const std::vector<std::vector<typename FEGenericBase<OutputShape>::OutputShape> > &curl_phi = fe->get_curl_phi();

  // Accumulate solution curl
  curl_u = 0.;

  for (unsigned int l=0; l != n_dofs; l++)
    curl_u.add_scaled(curl_phi[l][qp], coef(l));

  return;
}


template<typename OutputType>
void FEMContext::interior_div(unsigned int var, unsigned int qp,
                              OutputType& div_u) const
{
  typedef typename
    TensorTools::IncrementRank
    <typename TensorTools::MakeReal<OutputType>::type>::type OutputShape;

  // Get local-to-global dof index lookup
  libmesh_assert_greater (this->get_dof_indices().size(), var);
  const unsigned int n_dofs = cast_int<unsigned int>
    (this->get_dof_indices(var).size());

  // Get current local coefficients
  libmesh_assert_greater (this->_elem_subsolutions.size(), var);
  const DenseSubVector<Number> &coef = this->get_elem_solution(var);

  // Get finite element object
  FEGenericBase<OutputShape>* fe = NULL;
  this->get_element_fe<OutputShape>( var, fe );

  // Get shape function values at quadrature point
  const std::vector<std::vector<typename FEGenericBase<OutputShape>::OutputDivergence> > &div_phi = fe->get_div_phi();

  // Accumulate solution curl
  div_u = 0.;

  for (unsigned int l=0; l != n_dofs; l++)
    div_u += div_phi[l][qp] * coef(l);

  return;
}


Number FEMContext::side_value(unsigned int var, unsigned int qp) const
{
  Number u = 0.;

  this->side_value( var, qp, u );

  return u;
}


template<typename OutputType>
void FEMContext::side_value(unsigned int var, unsigned int qp,
                            OutputType& u) const
{
  this->some_value
    <OutputType,
     &FEMContext::get_side_fe
       <typename TensorTools::MakeReal<OutputType>::type>,
     &DiffContext::get_elem_solution>
    (var, qp, u);
}


template<typename OutputType>
void FEMContext::side_values
(unsigned int var,
 const NumericVector<Number> & _system_vector,
 std::vector<OutputType>& u_vals) const
{
  typedef typename TensorTools::MakeReal<OutputType>::type OutputShape;

  // Get local-to-global dof index lookup
  libmesh_assert_greater (this->get_dof_indices().size(), var);
  const unsigned int n_dofs = cast_int<unsigned int>
    (this->get_dof_indices(var).size());

  // Get current local coefficients
  const DenseSubVector<Number> &coef = get_localized_subvector(_system_vector, var);

  // Get the finite element object
  FEGenericBase<OutputShape>* the_side_fe = NULL;
  this->get_side_fe<OutputShape>( var, the_side_fe );

  // Get shape function values at quadrature point
  const std::vector<std::vector<OutputShape> > &phi = the_side_fe->get_phi();

  // Loop over all the q_points on this element
  for (unsigned int qp=0; qp != u_vals.size(); qp++)
    {
      OutputType &u = u_vals[qp];

      // Compute the value at this q_point
      u = 0.;

      for (unsigned int l=0; l != n_dofs; l++)
        u += phi[l][qp] * coef(l);
    }

  return;
}

Gradient FEMContext::side_gradient(unsigned int var, unsigned int qp) const
{
  Gradient du;

  this->side_gradient( var, qp, du );

  return du;
}


template<typename OutputType>
void FEMContext::side_gradient(unsigned int var, unsigned int qp,
                               OutputType& du) const
{
  typedef typename TensorTools::MakeReal
    <typename TensorTools::DecrementRank<OutputType>::type>::type
    OutputShape;

  // Get local-to-global dof index lookup
  libmesh_assert_greater (this->get_dof_indices().size(), var);
  const unsigned int n_dofs = cast_int<unsigned int>
    (this->get_dof_indices(var).size());

  // Get current local coefficients
  libmesh_assert_greater (this->_elem_subsolutions.size(), var);
  const DenseSubVector<Number> &coef = this->get_elem_solution(var);

  // Get finite element object
  FEGenericBase<OutputShape>* the_side_fe = NULL;
  this->get_side_fe<OutputShape>( var, the_side_fe );

  // Get shape function values at quadrature point
  const std::vector<std::vector< typename FEGenericBase<OutputShape>::OutputGradient> > &dphi = the_side_fe->get_dphi();

  // Accumulate solution derivatives
  du = 0.;

  for (unsigned int l=0; l != n_dofs; l++)
    du.add_scaled(dphi[l][qp], coef(l));

  return;
}



template<typename OutputType>
void FEMContext::side_gradients
(unsigned int var,
 const NumericVector<Number> & _system_vector,
 std::vector<OutputType>& du_vals) const
{
  typedef typename TensorTools::MakeReal
    <typename TensorTools::DecrementRank<OutputType>::type>::type
    OutputShape;

  // Get local-to-global dof index lookup
  libmesh_assert_greater (this->get_dof_indices().size(), var);
  const unsigned int n_dofs = cast_int<unsigned int>
    (this->get_dof_indices(var).size());

  // Get current local coefficients
  const DenseSubVector<Number> &coef = get_localized_subvector(_system_vector, var);

  // Get finite element object
  FEGenericBase<OutputShape>* the_side_fe = NULL;
  this->get_side_fe<OutputShape>( var, the_side_fe );

  // Get shape function values at quadrature point
  const std::vector<std::vector<typename FEGenericBase<OutputShape>::OutputGradient> > &dphi = the_side_fe->get_dphi();

  // Loop over all the q_points in this finite element
  for (unsigned int qp=0; qp != du_vals.size(); qp++)
    {
      OutputType &du = du_vals[qp];

      du = 0;

      // Compute the gradient at this q_point
      for (unsigned int l=0; l != n_dofs; l++)
        du.add_scaled(dphi[l][qp], coef(l));
    }

  return;
}

#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
Tensor FEMContext::side_hessian(unsigned int var, unsigned int qp) const
{
  Tensor d2u;

  this->side_hessian( var, qp, d2u );

  return d2u;
}



template<typename OutputType>
void FEMContext::side_hessian(unsigned int var, unsigned int qp,
                              OutputType& d2u) const
{
  this->some_hessian
    <OutputType,
     &FEMContext::get_side_fe
       <typename TensorTools::MakeReal
         <typename TensorTools::DecrementRank
           <typename TensorTools::DecrementRank
             <OutputType>::type>::type>::type>,
     &DiffContext::get_elem_solution>
    (var, qp, d2u);
}



template<typename OutputType>
void FEMContext::side_hessians
(unsigned int var,
 const NumericVector<Number> & _system_vector,
 std::vector<OutputType>& d2u_vals) const
{
  typedef typename TensorTools::MakeReal<
    typename TensorTools::DecrementRank<
    typename TensorTools::DecrementRank<
    OutputType>::type>::type>::type
    OutputShape;

  // Get local-to-global dof index lookup
  libmesh_assert_greater (this->get_dof_indices().size(), var);
  const unsigned int n_dofs = cast_int<unsigned int>
    (this->get_dof_indices(var).size());

  // Get current local coefficients
  const DenseSubVector<Number> &coef = get_localized_subvector(_system_vector, var);

  // Get finite element object
  FEGenericBase<OutputShape>* the_side_fe = NULL;
  this->get_side_fe<OutputShape>( var, the_side_fe );

  // Get shape function values at quadrature point
  const std::vector<std::vector<typename FEGenericBase<OutputShape>::OutputTensor> > &d2phi = the_side_fe->get_d2phi();

  // Loop over all the q_points in this finite element
  for (unsigned int qp=0; qp != d2u_vals.size(); qp++)
    {
      OutputType &d2u = d2u_vals[qp];

      // Compute the gradient at this q_point
      d2u = 0;

      for (unsigned int l=0; l != n_dofs; l++)
        d2u.add_scaled(d2phi[l][qp], coef(l));
    }

  return;
}



#endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES



Number FEMContext::point_value(unsigned int var, const Point &p) const
{
  Number u = 0.;

  this->point_value( var, p, u );

  return u;
}

template<typename OutputType>
void FEMContext::point_value(unsigned int var, const Point &p,
                             OutputType& u) const
{
  typedef typename TensorTools::MakeReal<OutputType>::type OutputShape;

  // Get local-to-global dof index lookup
  libmesh_assert_greater (this->get_dof_indices().size(), var);
  const unsigned int n_dofs = cast_int<unsigned int>
    (this->get_dof_indices(var).size());

  // Get current local coefficients
  libmesh_assert_greater (this->_elem_subsolutions.size(), var);
  const DenseSubVector<Number> &coef = this->get_elem_solution(var);

  // Get finite element object
  FEGenericBase<OutputShape>* fe = NULL;
  this->get_element_fe<OutputShape>( var, fe );

  // Build a FE for calculating u(p)
  UniquePtr<FEGenericBase<OutputShape> > fe_new = this->build_new_fe( fe, p );

  // Get the values of the shape function derivatives
  const std::vector<std::vector<OutputShape> >&  phi = fe_new->get_phi();

  u = 0.;

  for (unsigned int l=0; l != n_dofs; l++)
    u += phi[l][0] * coef(l);

  return;
}



Gradient FEMContext::point_gradient(unsigned int var, const Point &p) const
{
  Gradient grad_u;

  this->point_gradient( var, p, grad_u );

  return grad_u;
}



template<typename OutputType>
void FEMContext::point_gradient(unsigned int var, const Point &p,
                                OutputType& grad_u) const
{
  typedef typename TensorTools::MakeReal
    <typename TensorTools::DecrementRank<OutputType>::type>::type
    OutputShape;

  // Get local-to-global dof index lookup
  libmesh_assert_greater (this->get_dof_indices().size(), var);
  const unsigned int n_dofs = cast_int<unsigned int>
    (this->get_dof_indices(var).size());

  // Get current local coefficients
  libmesh_assert_greater (this->_elem_subsolutions.size(), var);
  const DenseSubVector<Number> &coef = this->get_elem_solution(var);

  // Get finite element object
  FEGenericBase<OutputShape>* fe = NULL;
  this->get_element_fe<OutputShape>( var, fe );

  // Build a FE for calculating u(p)
  UniquePtr<FEGenericBase<OutputShape> > fe_new = this->build_new_fe( fe, p );

  // Get the values of the shape function derivatives
  const std::vector<std::vector<typename FEGenericBase<OutputShape>::OutputGradient> >&  dphi = fe_new->get_dphi();

  grad_u = 0.0;

  for (unsigned int l=0; l != n_dofs; l++)
    grad_u.add_scaled(dphi[l][0], coef(l));

  return;
}



#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES

Tensor FEMContext::point_hessian(unsigned int var, const Point &p) const
{
  Tensor hess_u;

  this->point_hessian( var, p, hess_u );

  return hess_u;
}


template<typename OutputType>
void FEMContext::point_hessian(unsigned int var, const Point &p,
                               OutputType& hess_u) const
{
  typedef typename TensorTools::MakeReal<
    typename TensorTools::DecrementRank<
    typename TensorTools::DecrementRank<
    OutputType>::type>::type>::type
    OutputShape;

  // Get local-to-global dof index lookup
  libmesh_assert_greater (this->get_dof_indices().size(), var);
  const unsigned int n_dofs = cast_int<unsigned int>
    (this->get_dof_indices(var).size());

  // Get current local coefficients
  libmesh_assert_greater (this->_elem_subsolutions.size(), var);
  const DenseSubVector<Number> &coef = this->get_elem_solution(var);

  // Get finite element object
  FEGenericBase<OutputShape>* fe = NULL;
  this->get_element_fe<OutputShape>( var, fe );

  // Build a FE for calculating u(p)
  UniquePtr<FEGenericBase<OutputShape> > fe_new = this->build_new_fe( fe, p );

  // Get the values of the shape function derivatives
  const std::vector<std::vector<typename FEGenericBase<OutputShape>::OutputTensor> >&  d2phi = fe_new->get_d2phi();

  hess_u = 0.0;

  for (unsigned int l=0; l != n_dofs; l++)
    hess_u.add_scaled(d2phi[l][0], coef(l));

  return;
}

#endif // LIBMESH_ENABLE_SECOND_DERIVATIVES


template<typename OutputType>
void FEMContext::point_curl(unsigned int var, const Point &p,
                            OutputType& curl_u) const
{
  typedef typename TensorTools::MakeReal<OutputType>::type OutputShape;

  // Get local-to-global dof index lookup
  libmesh_assert_greater (this->get_dof_indices().size(), var);
  const unsigned int n_dofs = cast_int<unsigned int>
    (this->get_dof_indices(var).size());

  // Get current local coefficients
  libmesh_assert_greater (this->_elem_subsolutions.size(), var);
  const DenseSubVector<Number> &coef = this->get_elem_solution(var);

  // Get finite element object
  FEGenericBase<OutputShape>* fe = NULL;
  this->get_element_fe<OutputShape>( var, fe );

  // Build a FE for calculating u(p)
  UniquePtr<FEGenericBase<OutputShape> > fe_new = this->build_new_fe( fe, p );

  // Get the values of the shape function derivatives
  const std::vector<std::vector<typename FEGenericBase<OutputShape>::OutputShape> >&  curl_phi = fe_new->get_curl_phi();

  curl_u = 0.0;

  for (unsigned int l=0; l != n_dofs; l++)
    curl_u.add_scaled(curl_phi[l][0], coef(l));

  return;
}



Number FEMContext::fixed_interior_value(unsigned int var, unsigned int qp) const
{
  Number u = 0.;

  this->fixed_interior_value( var, qp, u );

  return u;
}



template<typename OutputType>
void FEMContext::fixed_interior_value(unsigned int var, unsigned int qp,
                                      OutputType& u) const
{
  this->some_value
    <OutputType,
     &FEMContext::get_element_fe
       <typename TensorTools::MakeReal<OutputType>::type>,
     &DiffContext::get_elem_fixed_solution>
    (var, qp, u);
}



Gradient FEMContext::fixed_interior_gradient(unsigned int var, unsigned int qp) const
{
  Gradient du;

  this->fixed_interior_gradient( var, qp, du );

  return du;
}


template<typename OutputType>
void FEMContext::FEMContext::fixed_interior_gradient(unsigned int var, unsigned int qp,
                                                     OutputType& du) const
{
  this->some_gradient
    <OutputType,
     &FEMContext::get_element_fe
       <typename TensorTools::MakeReal
         <typename TensorTools::DecrementRank
           <OutputType>::type>::type>,
     &DiffContext::get_elem_fixed_solution>
    (var, qp, du);
}



#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
Tensor FEMContext::fixed_interior_hessian(unsigned int var, unsigned int qp) const
{
  Tensor d2u;

  this->fixed_interior_hessian( var, qp, d2u );

  return d2u;
}


template<typename OutputType>
void FEMContext::fixed_interior_hessian(unsigned int var, unsigned int qp,
                                        OutputType& d2u) const
{
  this->some_hessian
    <OutputType,
     &FEMContext::get_element_fe
       <typename TensorTools::MakeReal
         <typename TensorTools::DecrementRank
           <typename TensorTools::DecrementRank
             <OutputType>::type>::type>::type>,
     &DiffContext::get_elem_fixed_solution>
    (var, qp, d2u);
}
#endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES



Number FEMContext::fixed_side_value(unsigned int var, unsigned int qp) const
{
  Number u = 0.;

  this->fixed_side_value( var, qp, u );

  return u;
}


template<typename OutputType>
void FEMContext::fixed_side_value(unsigned int var, unsigned int qp,
                                  OutputType& u) const
{
  this->some_value
    <OutputType,
     &FEMContext::get_side_fe
       <typename TensorTools::MakeReal<OutputType>::type>,
     &DiffContext::get_elem_fixed_solution>
    (var, qp, u);
}



Gradient FEMContext::fixed_side_gradient(unsigned int var, unsigned int qp) const
{
  Gradient du;

  this->fixed_side_gradient( var, qp, du );

  return du;
}


template<typename OutputType>
void FEMContext::FEMContext::fixed_side_gradient(unsigned int var, unsigned int qp,
                                                 OutputType& du) const
{
  this->some_gradient
    <OutputType,
     &FEMContext::get_side_fe
       <typename TensorTools::MakeReal
         <typename TensorTools::DecrementRank
           <OutputType>::type>::type>,
     &DiffContext::get_elem_fixed_solution>
    (var, qp, du);
}



#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
Tensor FEMContext::fixed_side_hessian(unsigned int var, unsigned int qp) const
{
  Tensor d2u;

  this->fixed_side_hessian( var, qp, d2u );

  return d2u;
}

template<typename OutputType>
void FEMContext::fixed_side_hessian(unsigned int var, unsigned int qp,
                                    OutputType& d2u) const
{
  this->some_hessian
    <OutputType,
     &FEMContext::get_side_fe
       <typename TensorTools::MakeReal
         <typename TensorTools::DecrementRank
           <typename TensorTools::DecrementRank
             <OutputType>::type>::type>::type>,
     &DiffContext::get_elem_fixed_solution>
    (var, qp, d2u);
}
#endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES



Number FEMContext::fixed_point_value(unsigned int var, const Point &p) const
{
  Number u = 0.;

  this->fixed_point_value( var, p, u );

  return u;
}

template<typename OutputType>
void FEMContext::fixed_point_value(unsigned int var, const Point &p,
                                   OutputType& u) const
{
  typedef typename TensorTools::MakeReal<OutputType>::type OutputShape;

  // Get local-to-global dof index lookup
  libmesh_assert_greater (this->get_dof_indices().size(), var);
  const unsigned int n_dofs = cast_int<unsigned int>
    (this->get_dof_indices(var).size());

  // Get current local coefficients
  libmesh_assert_greater (_elem_fixed_subsolutions.size(), var);
  const DenseSubVector<Number> &coef = this->get_elem_fixed_solution(var);

  // Get finite element object
  FEGenericBase<OutputShape>* fe = NULL;
  this->get_element_fe<OutputShape>( var, fe );

  // Build a FE for calculating u(p)
  UniquePtr<FEGenericBase<OutputShape> > fe_new = this->build_new_fe( fe, p );

  // Get the values of the shape function derivatives
  const std::vector<std::vector<OutputShape> >&  phi = fe_new->get_phi();

  u = 0.;

  for (unsigned int l=0; l != n_dofs; l++)
    u += phi[l][0] * coef(l);

  return;
}



Gradient FEMContext::fixed_point_gradient(unsigned int var, const Point &p) const
{
  Gradient grad_u;

  this->fixed_point_gradient( var, p, grad_u );

  return grad_u;
}



template<typename OutputType>
void FEMContext::fixed_point_gradient(unsigned int var, const Point &p,
                                      OutputType& grad_u) const
{
  typedef typename TensorTools::MakeReal
    <typename TensorTools::DecrementRank<OutputType>::type>::type
    OutputShape;

  // Get local-to-global dof index lookup
  libmesh_assert_greater (this->get_dof_indices().size(), var);
  const unsigned int n_dofs = cast_int<unsigned int>
    (this->get_dof_indices(var).size());

  // Get current local coefficients
  libmesh_assert_greater (_elem_fixed_subsolutions.size(), var);
  const DenseSubVector<Number> &coef = this->get_elem_fixed_solution(var);

  // Get finite element object
  FEGenericBase<OutputShape>* fe = NULL;
  this->get_element_fe<OutputShape>( var, fe );

  // Build a FE for calculating u(p)
  UniquePtr<FEGenericBase<OutputShape> > fe_new = this->build_new_fe( fe, p );

  // Get the values of the shape function derivatives
  const std::vector<std::vector<typename FEGenericBase<OutputShape>::OutputGradient> >&  dphi = fe_new->get_dphi();

  grad_u = 0.0;

  for (unsigned int l=0; l != n_dofs; l++)
    grad_u.add_scaled(dphi[l][0], coef(l));

  return;
}


#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES

Tensor FEMContext::fixed_point_hessian(unsigned int var, const Point &p) const
{
  Tensor hess_u;

  this->fixed_point_hessian( var, p, hess_u );

  return hess_u;
}



template<typename OutputType>
void FEMContext::fixed_point_hessian(unsigned int var, const Point &p,
                                     OutputType& hess_u) const
{
  typedef typename TensorTools::MakeReal<
    typename TensorTools::DecrementRank<
    typename TensorTools::DecrementRank<
    OutputType>::type>::type>::type
    OutputShape;

  // Get local-to-global dof index lookup
  libmesh_assert_greater (this->get_dof_indices().size(), var);
  const unsigned int n_dofs = cast_int<unsigned int>
    (this->get_dof_indices(var).size());

  // Get current local coefficients
  libmesh_assert_greater (_elem_fixed_subsolutions.size(), var);
  const DenseSubVector<Number> &coef = this->get_elem_fixed_solution(var);
{'fraction_non_alphanumeric': 0.11982327866075713, 'fraction_numerical': 0.004446263726957062, 'mean_word_length': 3.7957083869220876, 'pattern_counts': {'":': 0, '<': 535, '<?xml version=': 0, '>': 831, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 30, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '2612498', 'n_tokens_mistral': 24095, 'n_tokens_neox': 22987, 'n_words': 5488}