Assembly.h

Go to the documentation of this file.

//* This file is part of the MOOSE framework
//* https://mooseframework.inl.gov
//*
//* All rights reserved, see COPYRIGHT for full restrictions
//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
//*
//* Licensed under LGPL 2.1, please see LICENSE for details
//* https://www.gnu.org/licenses/lgpl-2.1.html

#pragma once

#include "DenseMatrix.h"
#include "MooseArray.h"
#include "MooseTypes.h"
#include "MooseVariableFE.h"
#include "MoosePassKey.h"
#include "ArbitraryQuadrature.h"

#include "libmesh/dense_vector.h"
#include "libmesh/enum_quadrature_type.h"
#include "libmesh/fe_type.h"
#include "libmesh/point.h"
#include "libmesh/fe_base.h"
#include "libmesh/numeric_vector.h"
#include "libmesh/elem_side_builder.h"

#include <unordered_map>

// libMesh forward declarations
namespace libMesh
{
class DofMap;
class CouplingMatrix;
class Elem;
template <typename>
class VectorValue;
typedef VectorValue<Real> RealVectorValue;
template <typename T>
class FEGenericBase;
typedef FEGenericBase<Real> FEBase;
typedef FEGenericBase<VectorValue<Real>> FEVectorBase;
class Node;
template <typename T>
class NumericVector;
template <typename T>
class SparseMatrix;
class StaticCondensation;
}

// MOOSE Forward Declares
class FaceInfo;
class MooseMesh;
class ArbitraryQuadrature;
class SystemBase;
class MooseVariableFieldBase;
class MooseVariableBase;
template <typename>
class MooseVariableFE;
class MooseVariableScalar;
typedef MooseVariableFE<Real> MooseVariable;
typedef MooseVariableFE<RealVectorValue> VectorMooseVariable;
typedef MooseVariableFE<RealEigenVector> ArrayMooseVariable;
class XFEMInterface;
class SubProblem;
class NodeFaceConstraint;

#ifdef MOOSE_KOKKOS_ENABLED
namespace Moose::Kokkos
{
class Assembly;
}
#endif

// Assembly.h does not import Moose.h nor libMeshReducedNamespace.h
using libMesh::FEBase;
using libMesh::FEFamily;
using libMesh::FEType;
using libMesh::FEVectorBase;
using libMesh::LAGRANGE_VEC;
using libMesh::Order;
using libMesh::QuadratureType;

template <typename P, typename C>
void coordTransformFactor(const SubProblem & s,
                          SubdomainID sub_id,
                          const P & point,
                          C & factor,
                          SubdomainID neighbor_sub_id = libMesh::Elem::invalid_subdomain_id);

template <typename P, typename C>
void coordTransformFactor(const MooseMesh & mesh,
                          SubdomainID sub_id,
                          const P & point,
                          C & factor,
                          SubdomainID neighbor_sub_id = libMesh::Elem::invalid_subdomain_id);

class Assembly
{
public:
  Assembly(SystemBase & sys, THREAD_ID tid);
  virtual ~Assembly();

  template <typename T>
  static const T * const & constify_ref(T * const & inref)
  {
    const T * const * ptr = &inref;
    return *ptr;
  }

  const FEBase * const & getFE(FEType type, unsigned int dim) const
  {
    buildFE(type);
    return constify_ref(_fe[dim][type]);
  }

  const FEBase * const & getFENeighbor(FEType type, unsigned int dim) const
  {
    buildNeighborFE(type);
    return constify_ref(_fe_neighbor[dim][type]);
  }

  const FEBase * const & getFEFace(FEType type, unsigned int dim) const
  {
    buildFaceFE(type);
    return constify_ref(_fe_face[dim][type]);
  }

  const FEBase * const & getFEFaceNeighbor(FEType type, unsigned int dim) const
  {
    buildFaceNeighborFE(type);
    return constify_ref(_fe_face_neighbor[dim][type]);
  }

  const FEVectorBase * const & getVectorFE(FEType type, unsigned int dim) const
  {
    buildVectorFE(type);
    return constify_ref(_vector_fe[dim][type]);
  }

  const FEVectorBase * const & getVectorFENeighbor(FEType type, unsigned int dim) const
  {
    buildVectorNeighborFE(type);
    return constify_ref(_vector_fe_neighbor[dim][type]);
  }

  const FEVectorBase * const & getVectorFEFace(FEType type, unsigned int dim) const
  {
    buildVectorFaceFE(type);
    return constify_ref(_vector_fe_face[dim][type]);
  }

  const FEVectorBase * const & getVectorFEFaceNeighbor(FEType type, unsigned int dim) const
  {
    buildVectorFaceNeighborFE(type);
    return constify_ref(_vector_fe_face_neighbor[dim][type]);
  }

#ifdef MOOSE_KOKKOS_ENABLED

  using InternalDataKey = Moose::PassKey<Moose::Kokkos::Assembly>;
#endif

  const libMesh::QBase * const & qRule() const { return constify_ref(_current_qrule); }

  libMesh::QBase * const & writeableQRule() { return _current_qrule; }

#ifdef MOOSE_KOKKOS_ENABLED

  libMesh::QBase * writeableQRule(unsigned int dim, SubdomainID block, InternalDataKey)
  {
    return qrules(dim, block).vol.get();
  }
#endif

  const MooseArray<Point> & qPoints() const { return _current_q_points; }

  const std::vector<Point> & qPointsMortar() const { return _fe_msm->get_xyz(); }

  const MooseArray<Point> & physicalPoints() const { return _current_physical_points; }

  const MooseArray<Real> & JxW() const { return _current_JxW; }

  const MooseArray<ADReal> & adJxW() const { return _ad_JxW; }

  const MooseArray<ADReal> & adJxWFace() const { return _ad_JxW_face; }

  const MooseArray<ADReal> & adCurvatures() const;

  const MooseArray<Real> & coordTransformation() const { return _coord; }

  const MooseArray<Real> & mortarCoordTransformation() const { return _coord_msm; }

  const MooseArray<ADReal> & adCoordTransformation() const
  {
    // Coord values for non-cartesian coordinate systems are functions of the locations of the
    // quadrature points in physical space. We also have no way of knowing whether this was called
    // from a volumetric or face object so we should set both volumetric and face xyz to true
    _calculate_xyz = true;
    _calculate_face_xyz = true;

    _calculate_ad_coord = true;
    return _ad_coord;
  }

  const Moose::CoordinateSystemType & coordSystem() const { return _coord_type; }

  const libMesh::QBase * const & qRuleFace() const { return constify_ref(_current_qrule_face); }

  libMesh::QBase * const & writeableQRuleFace() { return _current_qrule_face; }

#ifdef MOOSE_KOKKOS_ENABLED

  libMesh::QBase * writeableQRuleFace(unsigned int dim, SubdomainID block, InternalDataKey)
  {
    return qrules(dim, block).face.get();
  }
#endif

  const MooseArray<Point> & qPointsFace() const { return _current_q_points_face; }

  const MooseArray<Real> & JxWFace() const { return _current_JxW_face; }

  const MooseArray<Point> & normals() const { return _current_normals; }

  /***
   * Returns the array of normals for quadrature points on a current side
   */
  const std::vector<Eigen::Map<RealDIMValue>> & mappedNormals() const { return _mapped_normals; }

  const MooseArray<std::vector<Point>> & tangents() const { return _current_tangents; }

  unsigned int numExtraElemIntegers() const { return _extra_elem_ids.size() - 1; }

  const dof_id_type & extraElemID(unsigned int id) const
  {
    mooseAssert(id < _extra_elem_ids.size(), "An invalid extra element integer id");
    return _extra_elem_ids[id];
  }

  const dof_id_type & extraElemIDNeighbor(unsigned int id) const
  {
    mooseAssert(id < _neighbor_extra_elem_ids.size(), "An invalid extra element integer id");
    return _neighbor_extra_elem_ids[id];
  }

  const MooseArray<ADPoint> & adNormals() const { return _ad_normals; }

  const MooseArray<ADPoint> & adQPoints() const
  {
    _calculate_xyz = true;
    return _ad_q_points;
  }

  const MooseArray<ADPoint> & adQPointsFace() const
  {
    _calculate_face_xyz = true;
    return _ad_q_points_face;
  }

  template <bool is_ad>
  const MooseArray<Moose::GenericType<Point, is_ad>> & genericQPoints() const;

  const Elem * const & elem() const { return _current_elem; }

  const SubdomainID & currentSubdomainID() const { return _current_subdomain_id; }

  void setCurrentSubdomainID(SubdomainID i) { _current_subdomain_id = i; }

  const BoundaryID & currentBoundaryID() const { return _current_boundary_id; }

  void setCurrentBoundaryID(BoundaryID i) { _current_boundary_id = i; }

  const Real & elemVolume() const { return _current_elem_volume; }

  const unsigned int & side() const { return _current_side; }

  const unsigned int & neighborSide() const { return _current_neighbor_side; }

  const Elem * const & sideElem() const { return _current_side_elem; }

  const Real & sideElemVolume() const { return _current_side_volume; }

  const Elem * const & neighbor() const { return _current_neighbor_elem; }

  const Elem * const & lowerDElem() const { return _current_lower_d_elem; }

  const Elem * const & neighborLowerDElem() const { return _current_neighbor_lower_d_elem; }

  /*
   * @return The current lower-dimensional element volume
   */
  const Real & lowerDElemVolume() const;

  /*
   * @return The current neighbor lower-dimensional element volume
   */
  const Real & neighborLowerDElemVolume() const;

  const SubdomainID & currentNeighborSubdomainID() const { return _current_neighbor_subdomain_id; }

  void setCurrentNeighborSubdomainID(SubdomainID i) { _current_neighbor_subdomain_id = i; }

  const Real & neighborVolume()
  {
    _need_neighbor_elem_volume = true;
    return _current_neighbor_volume;
  }

  const libMesh::QBase * const & qRuleNeighbor() const
  {
    return constify_ref(_current_qrule_neighbor);
  }

  libMesh::QBase * const & writeableQRuleNeighbor() { return _current_qrule_neighbor; }

  const MooseArray<Real> & JxWNeighbor() const;

  const MooseArray<Point> & qPointsFaceNeighbor() const { return _current_q_points_face_neighbor; }

  const Node * const & node() const { return _current_node; }

  const Node * const & nodeNeighbor() const { return _current_neighbor_node; }

  void createQRules(QuadratureType type,
                    Order order,
                    Order volume_order,
                    Order face_order,
                    SubdomainID block,
                    bool allow_negative_qweights = true);

  void bumpVolumeQRuleOrder(Order volume_order, SubdomainID block);

  void bumpAllQRuleOrder(Order order, SubdomainID block);

  void setVolumeQRule(libMesh::QBase * qrule, unsigned int dim);

  void setFaceQRule(libMesh::QBase * qrule, unsigned int dim);

  void setMortarQRule(Order order);

  void activateDual() { _need_dual = true; }

  bool needDual() const { return _need_dual; }

  void clearCachedQRules();

private:
  void setLowerQRule(libMesh::QBase * qrule, unsigned int dim);

public:
  void setNeighborQRule(libMesh::QBase * qrule, unsigned int dim);

  void reinit(const Elem * elem);

  void setVolumeQRule(const Elem * elem);

  void reinitElemFaceRef(const Elem * elem,
                         unsigned int elem_side,
                         Real tolerance,
                         const std::vector<Point> * const pts = nullptr,
                         const std::vector<Real> * const weights = nullptr);

  void reinitNeighborFaceRef(const Elem * neighbor_elem,
                             unsigned int neighbor_side,
                             Real tolerance,
                             const std::vector<Point> * const pts,
                             const std::vector<Real> * const weights = nullptr);

  void reinitDual(const Elem * elem, const std::vector<Point> & pts, const std::vector<Real> & JxW);

  void reinitLowerDElem(const Elem * elem,
                        const std::vector<Point> * const pts = nullptr,
                        const std::vector<Real> * const weights = nullptr);

  void reinitNeighborLowerDElem(const Elem * elem);

  void reinitMortarElem(const Elem * elem);

  const std::vector<Real> & jxWMortar() const { return *_JxW_msm; }

  const libMesh::QBase * const & qRuleMortar() const { return constify_ref(_qrule_msm); }

private:
  void computeADFace(const Elem & elem, const unsigned int side);

public:
  void reinitAtPhysical(const Elem * elem, const std::vector<Point> & physical_points);

  void reinit(const Elem * elem, const std::vector<Point> & reference_points);

  void setFaceQRule(const Elem * const elem, const unsigned int side);

  void reinit(const Elem * elem, unsigned int side);

  void reinit(const Elem * elem, unsigned int side, const std::vector<Point> & reference_points);

  void reinitFVFace(const FaceInfo & fi);

  void reinitElemAndNeighbor(const Elem * elem,
                             unsigned int side,
                             const Elem * neighbor,
                             unsigned int neighbor_side,
                             const std::vector<Point> * neighbor_reference_points = nullptr);

  void reinitNeighborAtPhysical(const Elem * neighbor,
                                unsigned int neighbor_side,
                                const std::vector<Point> & physical_points);

  void reinitNeighborAtPhysical(const Elem * neighbor, const std::vector<Point> & physical_points);

  void reinitNeighbor(const Elem * neighbor, const std::vector<Point> & reference_points);

  void reinit(const Node * node);

  void init(const libMesh::CouplingMatrix * cm);

  void initNonlocalCoupling();

  void prepareJacobianBlock();

  void prepareResidual();

  void prepare();
  void prepareNonlocal();

  void prepareVariable(MooseVariableFieldBase * var);
  void prepareVariableNonlocal(MooseVariableFieldBase * var);
  void prepareNeighbor();

  void prepareLowerD();

  void prepareBlock(unsigned int ivar, unsigned jvar, const std::vector<dof_id_type> & dof_indices);
  void prepareBlockNonlocal(unsigned int ivar,
                            unsigned jvar,
                            const std::vector<dof_id_type> & idof_indices,
                            const std::vector<dof_id_type> & jdof_indices);
  void prepareScalar();
  void prepareOffDiagScalar();

  template <typename T>
  void copyShapes(MooseVariableField<T> & v);
  void copyShapes(unsigned int var);

  template <typename T>
  void copyFaceShapes(MooseVariableField<T> & v);
  void copyFaceShapes(unsigned int var);

  template <typename T>
  void copyNeighborShapes(MooseVariableField<T> & v);
  void copyNeighborShapes(unsigned int var);

  class GlobalDataKey
  {
    // Blessed classes
    friend class Assembly;
    friend class SubProblem;
    friend class FEProblemBase;
    friend class DisplacedProblem;
    friend class ComputeMortarFunctor;
    friend class NonlinearSystemBase;
    GlobalDataKey() {}
    GlobalDataKey(const GlobalDataKey &) {}
  };

  class LocalDataKey
  {
    // Blessed classes
    friend class Assembly;
    friend class TaggingInterface;
    LocalDataKey() {}
    LocalDataKey(const LocalDataKey &) {}
  };

  void addResidual(GlobalDataKey, const std::vector<VectorTag> & vector_tags);
  void addResidualNeighbor(GlobalDataKey, const std::vector<VectorTag> & vector_tags);
  void addResidualLower(GlobalDataKey, const std::vector<VectorTag> & vector_tags);

  void addResidualScalar(GlobalDataKey, const std::vector<VectorTag> & vector_tags);

  void cacheResidual(GlobalDataKey, const std::vector<VectorTag> & tags);

  void cacheResidualNeighbor(GlobalDataKey, const std::vector<VectorTag> & tags);

  void cacheResidualLower(GlobalDataKey, const std::vector<VectorTag> & tags);

  void addCachedResiduals(GlobalDataKey, const std::vector<VectorTag> & tags);

  void clearCachedResiduals(GlobalDataKey);

  void addCachedResidualDirectly(NumericVector<Number> & residual,
                                 GlobalDataKey,
                                 const VectorTag & vector_tag);

  void setResidual(NumericVector<Number> & residual, GlobalDataKey, const VectorTag & vector_tag);

  void setResidualNeighbor(NumericVector<Number> & residual,
                           GlobalDataKey,
                           const VectorTag & vector_tag);

  void addJacobian(GlobalDataKey);

  void addJacobianNonlocal(GlobalDataKey);

  void addJacobianNeighbor(GlobalDataKey);

  void addJacobianScalar(GlobalDataKey);

  void addJacobianOffDiagScalar(unsigned int ivar, GlobalDataKey);

  void addJacobianBlock(libMesh::SparseMatrix<Number> & jacobian,
                        unsigned int ivar,
                        unsigned int jvar,
                        const libMesh::DofMap & dof_map,
                        std::vector<dof_id_type> & dof_indices,
                        GlobalDataKey,
                        TagID tag);

  void addJacobianBlockTags(libMesh::SparseMatrix<Number> & jacobian,
                            unsigned int ivar,
                            unsigned int jvar,
                            const libMesh::DofMap & dof_map,
                            std::vector<dof_id_type> & dof_indices,
                            GlobalDataKey,
                            const std::set<TagID> & tags);

  void addJacobianBlockNonlocal(libMesh::SparseMatrix<Number> & jacobian,
                                unsigned int ivar,
                                unsigned int jvar,
                                const libMesh::DofMap & dof_map,
                                const std::vector<dof_id_type> & idof_indices,
                                const std::vector<dof_id_type> & jdof_indices,
                                GlobalDataKey,
                                TagID tag);

  void addJacobianBlockNonlocalTags(libMesh::SparseMatrix<Number> & jacobian,
                                    unsigned int ivar,
                                    unsigned int jvar,
                                    const libMesh::DofMap & dof_map,
                                    const std::vector<dof_id_type> & idof_indices,
                                    const std::vector<dof_id_type> & jdof_indices,
                                    GlobalDataKey,
                                    const std::set<TagID> & tags);

  void addJacobianNeighborLowerD(GlobalDataKey);

  void addJacobianLowerD(GlobalDataKey);

  void cacheJacobianMortar(GlobalDataKey);

  void addJacobianNeighbor(libMesh::SparseMatrix<Number> & jacobian,
                           unsigned int ivar,
                           unsigned int jvar,
                           const libMesh::DofMap & dof_map,
                           std::vector<dof_id_type> & dof_indices,
                           std::vector<dof_id_type> & neighbor_dof_indices,
                           GlobalDataKey,
                           TagID tag);

  void addJacobianNeighborTags(libMesh::SparseMatrix<Number> & jacobian,
                               unsigned int ivar,
                               unsigned int jvar,
                               const libMesh::DofMap & dof_map,
                               std::vector<dof_id_type> & dof_indices,
                               std::vector<dof_id_type> & neighbor_dof_indices,
                               GlobalDataKey,
                               const std::set<TagID> & tags);

  void cacheJacobian(GlobalDataKey);

  void cacheJacobianNonlocal(GlobalDataKey);

  void cacheJacobianNeighbor(GlobalDataKey);

  void addCachedJacobian(GlobalDataKey);

  void setCachedJacobian(GlobalDataKey);

  void zeroCachedJacobian(GlobalDataKey);

  DenseVector<Number> & residualBlock(unsigned int var_num, LocalDataKey, TagID tag_id)
  {
    return _sub_Re[tag_id][var_num];
  }

  DenseVector<Number> & residualBlockNeighbor(unsigned int var_num, LocalDataKey, TagID tag_id)
  {
    return _sub_Rn[tag_id][var_num];
  }

  DenseVector<Number> & residualBlockLower(unsigned int var_num, LocalDataKey, TagID tag_id)
  {
    return _sub_Rl[tag_id][var_num];
  }

  DenseMatrix<Number> & jacobianBlock(unsigned int ivar, unsigned int jvar, LocalDataKey, TagID tag)
  {
    jacobianBlockUsed(tag, ivar, jvar, true);
    return _sub_Kee[tag][ivar][_block_diagonal_matrix ? 0 : jvar];
  }

  DenseMatrix<Number> &
  jacobianBlockNonlocal(unsigned int ivar, unsigned int jvar, LocalDataKey, TagID tag)
  {
    jacobianBlockNonlocalUsed(tag, ivar, jvar, true);
    return _sub_Keg[tag][ivar][_block_diagonal_matrix ? 0 : jvar];
  }

  DenseMatrix<Number> & jacobianBlockNeighbor(
      Moose::DGJacobianType type, unsigned int ivar, unsigned int jvar, LocalDataKey, TagID tag);

  DenseMatrix<Number> & jacobianBlockMortar(Moose::ConstraintJacobianType type,
                                            unsigned int ivar,
                                            unsigned int jvar,
                                            LocalDataKey,
                                            TagID tag);

  void cacheResidualNodes(const DenseVector<Number> & res,
                          const std::vector<dof_id_type> & dof_index,
                          LocalDataKey,
                          TagID tag);

  void
  cacheJacobian(numeric_index_type i, numeric_index_type j, Real value, LocalDataKey, TagID tag);

  void cacheJacobian(numeric_index_type i,
                     numeric_index_type j,
                     Real value,
                     LocalDataKey,
                     const std::set<TagID> & tags);

  void cacheJacobianBlock(DenseMatrix<Number> & jac_block,
                          const std::vector<dof_id_type> & idof_indices,
                          const std::vector<dof_id_type> & jdof_indices,
                          Real scaling_factor,
                          LocalDataKey,
                          TagID tag);

  template <typename Residuals, typename Indices>
  void cacheResiduals(const Residuals & residuals,
                      const Indices & row_indices,
                      Real scaling_factor,
                      LocalDataKey,
                      const std::set<TagID> & vector_tags);

  template <typename Residuals, typename Indices>
  void cacheJacobian(const Residuals & residuals,
                     const Indices & row_indices,
                     Real scaling_factor,
                     LocalDataKey,
                     const std::set<TagID> & matrix_tags);

  template <typename Residuals, typename Indices>
  void cacheResidualsWithoutConstraints(const Residuals & residuals,
                                        const Indices & row_indices,
                                        Real scaling_factor,
                                        LocalDataKey,
                                        const std::set<TagID> & vector_tags);

  template <typename Residuals, typename Indices>
  void cacheJacobianWithoutConstraints(const Residuals & residuals,
                                       const Indices & row_indices,
                                       Real scaling_factor,
                                       LocalDataKey,
                                       const std::set<TagID> & matrix_tags);

  std::vector<std::pair<MooseVariableFieldBase *, MooseVariableFieldBase *>> & couplingEntries()
  {
    return _cm_ff_entry;
  }
  const std::vector<std::pair<MooseVariableFieldBase *, MooseVariableFieldBase *>> &
  couplingEntries() const
  {
    return _cm_ff_entry;
  }
  std::vector<std::pair<MooseVariableFieldBase *, MooseVariableFieldBase *>> &
  nonlocalCouplingEntries()
  {
    return _cm_nonlocal_entry;
  }
  const std::vector<std::pair<MooseVariableFieldBase *, MooseVariableScalar *>> &
  fieldScalarCouplingEntries() const
  {
    return _cm_fs_entry;
  }
  const std::vector<std::pair<MooseVariableScalar *, MooseVariableFieldBase *>> &
  scalarFieldCouplingEntries() const
  {
    return _cm_sf_entry;
  }

  // Read-only references
  const VariablePhiValue & phi() const { return _phi; }
  template <typename T>
  const ADTemplateVariablePhiGradient<T> & adGradPhi(const MooseVariableFE<T> & v) const
  {
    return _ad_grad_phi_data.at(v.feType());
  }
  const VariablePhiValue & phi(const MooseVariableField<Real> &) const { return _phi; }
  const VariablePhiGradient & gradPhi() const { return _grad_phi; }
  const VariablePhiGradient & gradPhi(const MooseVariableField<Real> &) const { return _grad_phi; }
  const VariablePhiSecond & secondPhi() const { return _second_phi; }
  const VariablePhiSecond & secondPhi(const MooseVariableField<Real> &) const
  {
    return _second_phi;
  }

  const VariablePhiValue & phiFace() const { return _phi_face; }
  const VariablePhiValue & phiFace(const MooseVariableField<Real> &) const { return _phi_face; }
  const VariablePhiGradient & gradPhiFace() const { return _grad_phi_face; }
  const VariablePhiGradient & gradPhiFace(const MooseVariableField<Real> &) const
  {
    return _grad_phi_face;
  }
  const VariablePhiSecond & secondPhiFace(const MooseVariableField<Real> &) const
  {
    return _second_phi_face;
  }

  const VariablePhiValue & phiNeighbor(const MooseVariableField<Real> &) const
  {
    return _phi_neighbor;
  }
  const VariablePhiGradient & gradPhiNeighbor(const MooseVariableField<Real> &) const
  {
    return _grad_phi_neighbor;
  }
  const VariablePhiSecond & secondPhiNeighbor(const MooseVariableField<Real> &) const
  {
    return _second_phi_neighbor;
  }

  const VariablePhiValue & phiFaceNeighbor(const MooseVariableField<Real> &) const
  {
    return _phi_face_neighbor;
  }
  const VariablePhiGradient & gradPhiFaceNeighbor(const MooseVariableField<Real> &) const
  {
    return _grad_phi_face_neighbor;
  }
  const VariablePhiSecond & secondPhiFaceNeighbor(const MooseVariableField<Real> &) const
  {
    return _second_phi_face_neighbor;
  }

  const VectorVariablePhiValue & phi(const MooseVariableField<RealVectorValue> &) const
  {
    return _vector_phi;
  }
  const VectorVariablePhiGradient & gradPhi(const MooseVariableField<RealVectorValue> &) const
  {
    return _vector_grad_phi;
  }
  const VectorVariablePhiSecond & secondPhi(const MooseVariableField<RealVectorValue> &) const
  {
    return _vector_second_phi;
  }
  const VectorVariablePhiCurl & curlPhi(const MooseVariableField<RealVectorValue> &) const
  {
    return _vector_curl_phi;
  }
  const VectorVariablePhiDivergence & divPhi(const MooseVariableField<RealVectorValue> &) const
  {
    return _vector_div_phi;
  }

  const VectorVariablePhiValue & phiFace(const MooseVariableField<RealVectorValue> &) const
  {
    return _vector_phi_face;
  }
  const VectorVariablePhiGradient & gradPhiFace(const MooseVariableField<RealVectorValue> &) const
  {
    return _vector_grad_phi_face;
  }
  const VectorVariablePhiSecond & secondPhiFace(const MooseVariableField<RealVectorValue> &) const
  {
    return _vector_second_phi_face;
  }
  const VectorVariablePhiCurl & curlPhiFace(const MooseVariableField<RealVectorValue> &) const
  {
    return _vector_curl_phi_face;
  }
  const VectorVariablePhiDivergence & divPhiFace(const MooseVariableField<RealVectorValue> &) const
  {
    return _vector_div_phi_face;
  }

  const VectorVariablePhiValue & phiNeighbor(const MooseVariableField<RealVectorValue> &) const
  {
    return _vector_phi_neighbor;
  }
  const VectorVariablePhiGradient &
  gradPhiNeighbor(const MooseVariableField<RealVectorValue> &) const
  {
    return _vector_grad_phi_neighbor;
  }
  const VectorVariablePhiSecond &
  secondPhiNeighbor(const MooseVariableField<RealVectorValue> &) const
  {
    return _vector_second_phi_neighbor;
  }
  const VectorVariablePhiCurl & curlPhiNeighbor(const MooseVariableField<RealVectorValue> &) const
  {
    return _vector_curl_phi_neighbor;
  }
  const VectorVariablePhiDivergence &
  divPhiNeighbor(const MooseVariableField<RealVectorValue> &) const
  {
    return _vector_div_phi_neighbor;
  }

  const VectorVariablePhiValue & phiFaceNeighbor(const MooseVariableField<RealVectorValue> &) const
  {
    return _vector_phi_face_neighbor;
  }
  const VectorVariablePhiGradient &
  gradPhiFaceNeighbor(const MooseVariableField<RealVectorValue> &) const
  {
    return _vector_grad_phi_face_neighbor;
  }
  const VectorVariablePhiSecond &
  secondPhiFaceNeighbor(const MooseVariableField<RealVectorValue> &) const
  {
    return _vector_second_phi_face_neighbor;
  }
  const VectorVariablePhiCurl &
  curlPhiFaceNeighbor(const MooseVariableField<RealVectorValue> &) const
  {
    return _vector_curl_phi_face_neighbor;
  }
  const VectorVariablePhiDivergence &
  divPhiFaceNeighbor(const MooseVariableField<RealVectorValue> &) const
  {
    return _vector_div_phi_face_neighbor;
  }

  // Writeable references
  VariablePhiValue & phi(const MooseVariableField<Real> &) { return _phi; }
  VariablePhiGradient & gradPhi(const MooseVariableField<Real> &) { return _grad_phi; }
  VariablePhiSecond & secondPhi(const MooseVariableField<Real> &) { return _second_phi; }

  VariablePhiValue & phiFace(const MooseVariableField<Real> &) { return _phi_face; }
  VariablePhiGradient & gradPhiFace(const MooseVariableField<Real> &) { return _grad_phi_face; }
  VariablePhiSecond & secondPhiFace(const MooseVariableField<Real> &) { return _second_phi_face; }

  VariablePhiValue & phiNeighbor(const MooseVariableField<Real> &) { return _phi_neighbor; }
  VariablePhiGradient & gradPhiNeighbor(const MooseVariableField<Real> &)
  {
    return _grad_phi_neighbor;
  }
  VariablePhiSecond & secondPhiNeighbor(const MooseVariableField<Real> &)
  {
    return _second_phi_neighbor;
  }

  VariablePhiValue & phiFaceNeighbor(const MooseVariableField<Real> &)
  {
    return _phi_face_neighbor;
  }
  VariablePhiGradient & gradPhiFaceNeighbor(const MooseVariableField<Real> &)
  {
    return _grad_phi_face_neighbor;
  }
  VariablePhiSecond & secondPhiFaceNeighbor(const MooseVariableField<Real> &)
  {
    return _second_phi_face_neighbor;
  }

  // Writeable references with vector variable
  VectorVariablePhiValue & phi(const MooseVariableField<RealVectorValue> &) { return _vector_phi; }
  VectorVariablePhiGradient & gradPhi(const MooseVariableField<RealVectorValue> &)
  {
    return _vector_grad_phi;
  }
  VectorVariablePhiSecond & secondPhi(const MooseVariableField<RealVectorValue> &)
  {
    return _vector_second_phi;
  }
  VectorVariablePhiCurl & curlPhi(const MooseVariableField<RealVectorValue> &)
  {
    return _vector_curl_phi;
  }
  VectorVariablePhiDivergence & divPhi(const MooseVariableField<RealVectorValue> &)
  {
    return _vector_div_phi;
  }

  VectorVariablePhiValue & phiFace(const MooseVariableField<RealVectorValue> &)
  {
    return _vector_phi_face;
  }
  VectorVariablePhiGradient & gradPhiFace(const MooseVariableField<RealVectorValue> &)
  {
    return _vector_grad_phi_face;
  }
  VectorVariablePhiSecond & secondPhiFace(const MooseVariableField<RealVectorValue> &)
  {
    return _vector_second_phi_face;
  }
  VectorVariablePhiCurl & curlPhiFace(const MooseVariableField<RealVectorValue> &)
  {
    return _vector_curl_phi_face;
  }
  VectorVariablePhiDivergence & divPhiFace(const MooseVariableField<RealVectorValue> &)
  {
    return _vector_div_phi_face;
  }

  VectorVariablePhiValue & phiNeighbor(const MooseVariableField<RealVectorValue> &)
  {
    return _vector_phi_neighbor;
  }
  VectorVariablePhiGradient & gradPhiNeighbor(const MooseVariableField<RealVectorValue> &)
  {
    return _vector_grad_phi_neighbor;
  }
  VectorVariablePhiSecond & secondPhiNeighbor(const MooseVariableField<RealVectorValue> &)
  {
    return _vector_second_phi_neighbor;
  }
  VectorVariablePhiCurl & curlPhiNeighbor(const MooseVariableField<RealVectorValue> &)
  {
    return _vector_curl_phi_neighbor;
  }
  VectorVariablePhiDivergence & divPhiNeighbor(const MooseVariableField<RealVectorValue> &)
  {
    return _vector_div_phi_neighbor;
  }
  VectorVariablePhiValue & phiFaceNeighbor(const MooseVariableField<RealVectorValue> &)
  {
    return _vector_phi_face_neighbor;
  }
  VectorVariablePhiGradient & gradPhiFaceNeighbor(const MooseVariableField<RealVectorValue> &)
  {
    return _vector_grad_phi_face_neighbor;
  }
  VectorVariablePhiSecond & secondPhiFaceNeighbor(const MooseVariableField<RealVectorValue> &)
  {
    return _vector_second_phi_face_neighbor;
  }
  VectorVariablePhiCurl & curlPhiFaceNeighbor(const MooseVariableField<RealVectorValue> &)
  {
    return _vector_curl_phi_face_neighbor;
  }
  VectorVariablePhiDivergence & divPhiFaceNeighbor(const MooseVariableField<RealVectorValue> &)
  {
    return _vector_div_phi_face_neighbor;
  }

  // Writeable references with array variable
  VariablePhiValue & phi(const MooseVariableField<RealEigenVector> &) { return _phi; }
  VariablePhiGradient & gradPhi(const MooseVariableField<RealEigenVector> &) { return _grad_phi; }
  VariablePhiSecond & secondPhi(const MooseVariableField<RealEigenVector> &) { return _second_phi; }

  VariablePhiValue & phiFace(const MooseVariableField<RealEigenVector> &) { return _phi_face; }
  VariablePhiGradient & gradPhiFace(const MooseVariableField<RealEigenVector> &)
  {
    return _grad_phi_face;
  }
  VariablePhiSecond & secondPhiFace(const MooseVariableField<RealEigenVector> &)
  {
    return _second_phi_face;
  }

  VariablePhiValue & phiNeighbor(const MooseVariableField<RealEigenVector> &)
  {
    return _phi_neighbor;
  }
  VariablePhiGradient & gradPhiNeighbor(const MooseVariableField<RealEigenVector> &)
  {
    return _grad_phi_neighbor;
  }
  VariablePhiSecond & secondPhiNeighbor(const MooseVariableField<RealEigenVector> &)
  {
    return _second_phi_neighbor;
  }

  VariablePhiValue & phiFaceNeighbor(const MooseVariableField<RealEigenVector> &)
  {
    return _phi_face_neighbor;
  }
  VariablePhiGradient & gradPhiFaceNeighbor(const MooseVariableField<RealEigenVector> &)
  {
    return _grad_phi_face_neighbor;
  }
  VariablePhiSecond & secondPhiFaceNeighbor(const MooseVariableField<RealEigenVector> &)
  {
    return _second_phi_face_neighbor;
  }

  template <typename OutputType>
  const typename OutputTools<OutputType>::VariablePhiValue & fePhi(FEType type) const
  {
    buildFE(type);
    return _fe_shape_data[type]->_phi;
  }

  template <typename OutputType>
  const typename OutputTools<OutputType>::VariablePhiGradient & feGradPhi(FEType type) const
  {
    buildFE(type);
    return _fe_shape_data[type]->_grad_phi;
  }

  template <typename OutputType>
  const ADTemplateVariablePhiGradient<OutputType> & feADGradPhi(FEType type) const
  {
    return _ad_grad_phi_data[type];
  }

  template <typename OutputType>
  const typename OutputTools<OutputType>::VariablePhiSecond & feSecondPhi(FEType type) const
  {
    _need_second_derivative.insert(type);
    buildFE(type);
    return _fe_shape_data[type]->_second_phi;
  }

  template <typename OutputType>
  const typename OutputTools<OutputType>::VariablePhiValue & fePhiLower(FEType type) const;

  template <typename OutputType>
  const typename OutputTools<OutputType>::VariablePhiValue & feDualPhiLower(FEType type) const;

  template <typename OutputType>
  const typename OutputTools<OutputType>::VariablePhiGradient & feGradPhiLower(FEType type) const;

  template <typename OutputType>
  const typename OutputTools<OutputType>::VariablePhiGradient &
  feGradDualPhiLower(FEType type) const;

  template <typename OutputType>
  const typename OutputTools<OutputType>::VariablePhiValue & fePhiFace(FEType type) const
  {
    buildFaceFE(type);
    return _fe_shape_data_face[type]->_phi;
  }

  template <typename OutputType>
  const typename OutputTools<OutputType>::VariablePhiGradient & feGradPhiFace(FEType type) const
  {
    buildFaceFE(type);
    return _fe_shape_data_face[type]->_grad_phi;
  }

  template <typename OutputType>
  const ADTemplateVariablePhiGradient<OutputType> & feADGradPhiFace(FEType type) const
  {
    return _ad_grad_phi_data_face[type];
  }

  template <typename OutputType>
  const typename OutputTools<OutputType>::VariablePhiSecond & feSecondPhiFace(FEType type) const
  {
    _need_second_derivative.insert(type);
    buildFaceFE(type);
    return _fe_shape_data_face[type]->_second_phi;
  }

  template <typename OutputType>
  const typename OutputTools<OutputType>::VariablePhiValue & fePhiNeighbor(FEType type) const
  {
    buildNeighborFE(type);
    return _fe_shape_data_neighbor[type]->_phi;
  }

  template <typename OutputType>
  const typename OutputTools<OutputType>::VariablePhiGradient & feGradPhiNeighbor(FEType type) const
  {
    buildNeighborFE(type);
    return _fe_shape_data_neighbor[type]->_grad_phi;
  }

  template <typename OutputType>
  const typename OutputTools<OutputType>::VariablePhiSecond & feSecondPhiNeighbor(FEType type) const
  {
    _need_second_derivative_neighbor.insert(type);
    buildNeighborFE(type);
    return _fe_shape_data_neighbor[type]->_second_phi;
  }

  template <typename OutputType>
  const typename OutputTools<OutputType>::VariablePhiValue & fePhiFaceNeighbor(FEType type) const
  {
    buildFaceNeighborFE(type);
    return _fe_shape_data_face_neighbor[type]->_phi;
  }

  template <typename OutputType>
  const typename OutputTools<OutputType>::VariablePhiGradient &
  feGradPhiFaceNeighbor(FEType type) const
  {
    buildFaceNeighborFE(type);
    return _fe_shape_data_face_neighbor[type]->_grad_phi;
  }

  template <typename OutputType>
  const typename OutputTools<OutputType>::VariablePhiSecond &
  feSecondPhiFaceNeighbor(FEType type) const
  {
    _need_second_derivative_neighbor.insert(type);
    buildFaceNeighborFE(type);
    return _fe_shape_data_face_neighbor[type]->_second_phi;
  }

  template <typename OutputType>
  const typename OutputTools<OutputType>::VariablePhiCurl & feCurlPhi(FEType type) const
  {
    _need_curl.insert(type);
    buildFE(type);
    return _fe_shape_data[type]->_curl_phi;
  }

  template <typename OutputType>
  const typename OutputTools<OutputType>::VariablePhiCurl & feCurlPhiFace(FEType type) const
  {
    _need_curl.insert(type);
    buildFaceFE(type);
    return _fe_shape_data_face[type]->_curl_phi;
  }

  template <typename OutputType>
  const typename OutputTools<OutputType>::VariablePhiCurl & feCurlPhiNeighbor(FEType type) const
  {
    _need_curl.insert(type);
    buildNeighborFE(type);
    return _fe_shape_data_neighbor[type]->_curl_phi;
  }

  template <typename OutputType>
  const typename OutputTools<OutputType>::VariablePhiCurl & feCurlPhiFaceNeighbor(FEType type) const
  {
    _need_curl.insert(type);
    buildFaceNeighborFE(type);
    return _fe_shape_data_face_neighbor[type]->_curl_phi;
  }

  template <typename OutputType>
  const typename OutputTools<OutputType>::VariablePhiDivergence & feDivPhi(FEType type) const
  {
    buildFE(type);
    return _fe_shape_data[type]->_div_phi;
  }

  template <typename OutputType>
  const typename OutputTools<OutputType>::VariablePhiDivergence & feDivPhiFace(FEType type) const
  {
    buildFaceFE(type);
    return _fe_shape_data_face[type]->_div_phi;
  }

  template <typename OutputType>
  const typename OutputTools<OutputType>::VariablePhiDivergence &
  feDivPhiNeighbor(FEType type) const
  {
    buildNeighborFE(type);
    return _fe_shape_data_neighbor[type]->_div_phi;
  }

  template <typename OutputType>
  const typename OutputTools<OutputType>::VariablePhiDivergence &
  feDivPhiFaceNeighbor(FEType type) const
  {
    buildFaceNeighborFE(type);
    return _fe_shape_data_face_neighbor[type]->_div_phi;
  }

  Real elementVolume(const Elem * elem) const;

  void setXFEM(std::shared_ptr<XFEMInterface> xfem) { _xfem = xfem; }

  void assignDisplacements(
      std::vector<std::pair<unsigned int, unsigned short>> && disp_numbers_and_directions);

  void saveLocalArrayResidual(DenseVector<Number> & re,
                              unsigned int i,
                              unsigned int ntest,
                              const RealEigenVector & v) const
  {
    for (unsigned int j = 0; j < v.size(); ++j, i += ntest)
      re(i) += v(j);
  }

  void saveLocalADArray(std::vector<ADReal> & re,
                        unsigned int i,
                        unsigned int ntest,
                        const ADRealEigenVector & v) const;

  void saveDiagLocalArrayJacobian(DenseMatrix<Number> & ke,
                                  unsigned int i,
                                  unsigned int ntest,
                                  unsigned int j,
                                  unsigned int nphi,
                                  unsigned int ivar,
                                  const RealEigenVector & v) const
  {
    unsigned int pace = (_component_block_diagonal[ivar] ? 0 : nphi);
    for (unsigned int k = 0; k < v.size(); ++k, i += ntest, j += pace)
      ke(i, j) += v(k);
  }

  void saveFullLocalArrayJacobian(DenseMatrix<Number> & ke,
                                  unsigned int i,
                                  unsigned int ntest,
                                  unsigned int j,
                                  unsigned int nphi,
                                  unsigned int ivar,
                                  unsigned int jvar,
                                  const RealEigenMatrix & v) const
  {
    if (ivar == jvar && _component_block_diagonal[ivar])
    {
      for (unsigned int k = 0; k < v.rows(); ++k, i += ntest)
        ke(i, j) += v(k, k);
    }
    else
    {
      const unsigned int saved_j = j;
      for (unsigned int k = 0; k < v.rows(); ++k, i += ntest)
      {
        j = saved_j;
        for (unsigned int l = 0; l < v.cols(); ++l, j += nphi)
          ke(i, j) += v(k, l);
      }
    }
  }

  DenseVector<Real> getJacobianDiagonal(DenseMatrix<Number> & ke)
  {
    unsigned int rows = ke.m();
    unsigned int cols = ke.n();
    DenseVector<Real> diag(rows);
    for (unsigned int i = 0; i < rows; i++)
      // % operation is needed to account for cases of no component coupling of array variables
      diag(i) = ke(i, i % cols);
    return diag;
  }

  inline const libMesh::QBase * attachQRuleElem(unsigned int dim, FEBase & fe)
  {
    auto qrule = qrules(dim).vol.get();
    fe.attach_quadrature_rule(qrule);
    return qrule;
  }

  inline const libMesh::QBase * attachQRuleFace(unsigned int dim, FEBase & fe)
  {
    auto qrule = qrules(dim).face.get();
    fe.attach_quadrature_rule(qrule);
    return qrule;
  }

  void hasScalingVector();

  void modifyArbitraryWeights(const std::vector<Real> & weights);

  bool computingResidual() const { return _computing_residual; }

  bool computingJacobian() const { return _computing_jacobian; }

  bool computingResidualAndJacobian() const { return _computing_residual_and_jacobian; }

  const Elem * const & msmElem() const { return _msm_elem; }

  void havePRefinement(const std::unordered_set<FEFamily> & disable_p_refinement_for_families);

  void setCurrentLowerDElem(const Elem * const lower_d_elem);

private:
  void reinitFE(const Elem * elem);

  void reinitFEFace(const Elem * elem, unsigned int side);

  void computeFaceMap(const Elem & elem, const unsigned int side, const std::vector<Real> & qw);

  void reinitFEFaceNeighbor(const Elem * neighbor, const std::vector<Point> & reference_points);

  void reinitFENeighbor(const Elem * neighbor, const std::vector<Point> & reference_points);

  template <typename Points, typename Coords>
  void setCoordinateTransformation(const libMesh::QBase * qrule,
                                   const Points & q_points,
                                   Coords & coord,
                                   SubdomainID sub_id);

  void computeCurrentElemVolume();

  void computeCurrentFaceVolume();

  void computeCurrentNeighborVolume();

  void modifyWeightsDueToXFEM(const Elem * elem);

  void modifyFaceWeightsDueToXFEM(const Elem * elem, unsigned int side = 0);

  template <typename OutputType>
  void computeGradPhiAD(const Elem * elem,
                        unsigned int n_qp,
                        ADTemplateVariablePhiGradient<OutputType> & grad_phi,
                        libMesh::FEGenericBase<OutputType> * fe);

  void resizeADMappingObjects(unsigned int n_qp, unsigned int dim);

  void
  computeSinglePointMapAD(const Elem * elem, const std::vector<Real> & qw, unsigned p, FEBase * fe);

  void addResidual(const VectorTag & vector_tag);
  void addResidualNeighbor(const VectorTag & vector_tag);
  void addResidualLower(const VectorTag & vector_tag);
  void addResidualScalar(const VectorTag & vector_tag);

  void clearCachedResiduals(const VectorTag & vector_tag);

  void cacheResidual(dof_id_type dof, Real value, TagID tag_id);

  void cacheResidual(dof_id_type dof, Real value, const std::set<TagID> & tags);

  void processLocalResidual(DenseVector<Number> & res_block,
                            std::vector<dof_id_type> & dof_indices,
                            const std::vector<Real> & scaling_factor);

  void addResidualBlock(NumericVector<Number> & residual,
                        DenseVector<Number> & res_block,
                        const std::vector<dof_id_type> & dof_indices,
                        const std::vector<Real> & scaling_factor);

  void cacheResidualBlock(std::vector<Real> & cached_residual_values,
                          std::vector<dof_id_type> & cached_residual_rows,
                          DenseVector<Number> & res_block,
                          const std::vector<dof_id_type> & dof_indices,
                          const std::vector<Real> & scaling_factor);

  void setResidualBlock(NumericVector<Number> & residual,
                        DenseVector<Number> & res_block,
                        const std::vector<dof_id_type> & dof_indices,
                        const std::vector<Real> & scaling_factor);

  void addJacobianBlock(libMesh::SparseMatrix<Number> & jacobian,
                        DenseMatrix<Number> & jac_block,
                        const MooseVariableBase & ivar,
                        const MooseVariableBase & jvar,
                        const std::vector<dof_id_type> & idof_indices,
                        const std::vector<dof_id_type> & jdof_indices);

  void cacheJacobianBlock(DenseMatrix<Number> & jac_block,
                          const MooseVariableBase & ivar,
                          const MooseVariableBase & jvar,
                          const std::vector<dof_id_type> & idof_indices,
                          const std::vector<dof_id_type> & jdof_indices,
                          TagID tag);

  void cacheJacobianBlockNonzero(DenseMatrix<Number> & jac_block,
                                 const MooseVariableBase & ivar,
                                 const MooseVariableBase & jvar,
                                 const std::vector<dof_id_type> & idof_indices,
                                 const std::vector<dof_id_type> & jdof_indices,
                                 TagID tag);

  void addJacobianCoupledVarPair(const MooseVariableBase & ivar, const MooseVariableBase & jvar);

  void cacheJacobianCoupledVarPair(const MooseVariableBase & ivar, const MooseVariableBase & jvar);

  void clearCachedJacobian();

  void buildFE(FEType type) const;

  void buildFaceFE(FEType type) const;

  void buildNeighborFE(FEType type) const;

  void buildFaceNeighborFE(FEType type) const;

  void buildLowerDFE(FEType type) const;

  void buildLowerDDualFE(FEType type) const;

  void buildVectorFE(FEType type) const;

  void buildVectorFaceFE(FEType type) const;

  void buildVectorNeighborFE(FEType type) const;

  void buildVectorFaceNeighborFE(FEType type) const;

  void buildVectorLowerDFE(FEType type) const;
  void buildVectorDualLowerDFE(FEType type) const;

  void jacobianBlockUsed(TagID tag, unsigned int ivar, unsigned int jvar, bool used)
  {
    _jacobian_block_used[tag][ivar][_block_diagonal_matrix ? 0 : jvar] = used;
  }

  char jacobianBlockUsed(TagID tag, unsigned int ivar, unsigned int jvar) const
  {
    return _jacobian_block_used[tag][ivar][_block_diagonal_matrix ? 0 : jvar];
  }

  void jacobianBlockNeighborUsed(TagID tag, unsigned int ivar, unsigned int jvar, bool used)
  {
    _jacobian_block_neighbor_used[tag][ivar][_block_diagonal_matrix ? 0 : jvar] = used;
  }

  char jacobianBlockNeighborUsed(TagID tag, unsigned int ivar, unsigned int jvar) const
  {
    return _jacobian_block_neighbor_used[tag][ivar][_block_diagonal_matrix ? 0 : jvar];
  }

  void jacobianBlockLowerUsed(TagID tag, unsigned int ivar, unsigned int jvar, bool used)
  {
    _jacobian_block_lower_used[tag][ivar][_block_diagonal_matrix ? 0 : jvar] = used;
  }

  char jacobianBlockLowerUsed(TagID tag, unsigned int ivar, unsigned int jvar) const
  {
    return _jacobian_block_lower_used[tag][ivar][_block_diagonal_matrix ? 0 : jvar];
  }

  void jacobianBlockNonlocalUsed(TagID tag, unsigned int ivar, unsigned int jvar, bool used)
  {
    _jacobian_block_nonlocal_used[tag][ivar][_block_diagonal_matrix ? 0 : jvar] = used;
  }

  char jacobianBlockNonlocalUsed(TagID tag, unsigned int ivar, unsigned int jvar) const
  {
    return _jacobian_block_nonlocal_used[tag][ivar][_block_diagonal_matrix ? 0 : jvar];
  }

  void helpersRequestData();

  SystemBase & _sys;
  SubProblem & _subproblem;

  const bool _displaced;

  const libMesh::CouplingMatrix * _cm;
  const libMesh::CouplingMatrix & _nonlocal_cm;

  const bool & _computing_residual;

  const bool & _computing_jacobian;

  const bool & _computing_residual_and_jacobian;

  std::vector<std::pair<MooseVariableFieldBase *, MooseVariableFieldBase *>> _cm_ff_entry;
  std::vector<std::pair<MooseVariableFieldBase *, MooseVariableScalar *>> _cm_fs_entry;
  std::vector<std::pair<MooseVariableScalar *, MooseVariableFieldBase *>> _cm_sf_entry;
  std::vector<std::pair<MooseVariableScalar *, MooseVariableScalar *>> _cm_ss_entry;
  std::vector<std::pair<MooseVariableFieldBase *, MooseVariableFieldBase *>> _cm_nonlocal_entry;
  std::vector<std::vector<std::vector<unsigned char>>> _jacobian_block_used;
  std::vector<std::vector<std::vector<unsigned char>>> _jacobian_block_nonlocal_used;
  std::vector<std::vector<std::vector<unsigned char>>> _jacobian_block_neighbor_used;
  std::vector<std::vector<std::vector<unsigned char>>> _jacobian_block_lower_used;
  const libMesh::DofMap & _dof_map;
  THREAD_ID _tid;

  MooseMesh & _mesh;

  unsigned int _mesh_dimension;

  const FEType _helper_type;

  mutable bool _user_added_fe_of_helper_type;
  mutable bool _user_added_fe_face_of_helper_type;
  mutable bool _user_added_fe_face_neighbor_of_helper_type;
  mutable bool _user_added_fe_neighbor_of_helper_type;
  mutable bool _user_added_fe_lower_of_helper_type;

  std::vector<std::unique_ptr<FEBase>> _unique_fe_helper;
  std::vector<std::unique_ptr<FEBase>> _unique_fe_face_helper;
  std::vector<std::unique_ptr<FEBase>> _unique_fe_face_neighbor_helper;
  std::vector<std::unique_ptr<FEBase>> _unique_fe_neighbor_helper;
  std::vector<std::unique_ptr<FEBase>> _unique_fe_lower_helper;

  bool _building_helpers;

  std::shared_ptr<XFEMInterface> _xfem;

  std::map<FEType, FEBase *> _current_fe;
  std::map<FEType, FEBase *> _current_fe_face;
  std::map<FEType, FEBase *> _current_fe_neighbor;
  std::map<FEType, FEBase *> _current_fe_face_neighbor;

  std::map<FEType, FEVectorBase *> _current_vector_fe;
  std::map<FEType, FEVectorBase *> _current_vector_fe_face;
  std::map<FEType, FEVectorBase *> _current_vector_fe_neighbor;
  std::map<FEType, FEVectorBase *> _current_vector_fe_face_neighbor;

  /**** Volume Stuff ****/

  mutable std::map<unsigned int, std::map<FEType, FEBase *>> _fe;
  mutable std::map<unsigned int, std::map<FEType, FEVectorBase *>> _vector_fe;
  std::map<unsigned int, FEBase *> _holder_fe_helper;
  FEBase * _current_fe_helper;
  libMesh::QBase * _current_qrule;
  libMesh::QBase * _current_qrule_volume;
  ArbitraryQuadrature * _current_qrule_arbitrary;
  ArbitraryQuadrature * _current_qrule_arbitrary_face;
  MooseArray<Point> _current_q_points;
  MooseArray<Real> _current_JxW;
  Moose::CoordinateSystemType _coord_type;
  MooseArray<Real> _coord;
  MooseArray<ADReal> _ad_coord;

  struct QRules
  {
    QRules()
      : vol(nullptr),
        face(nullptr),
        arbitrary_vol(nullptr),
        arbitrary_face(nullptr),
        neighbor(nullptr)
    {
    }

    std::unique_ptr<libMesh::QBase> vol;
    std::unique_ptr<libMesh::QBase> face;
    std::unique_ptr<libMesh::QBase> fv_face;
    std::unique_ptr<ArbitraryQuadrature> arbitrary_vol;
    std::unique_ptr<ArbitraryQuadrature> arbitrary_face;
    std::unique_ptr<ArbitraryQuadrature> neighbor;
  };

  std::unordered_map<SubdomainID, std::vector<QRules>> _qrules;

  libMesh::QBase * qruleFace(const Elem * elem, unsigned int side);
  ArbitraryQuadrature * qruleArbitraryFace(const Elem * elem, unsigned int side);

  template <typename T>
  T * qruleFaceHelper(const Elem * elem, unsigned int side, std::function<T *(QRules &)> rule_fn)
  {
    auto dim = elem->dim();
    auto neighbor = elem->neighbor_ptr(side);
    auto q = rule_fn(qrules(dim, elem->subdomain_id()));
    if (!neighbor)
      return q;

    // find the maximum face quadrature order for all blocks the face is in
    auto neighbor_block = neighbor->subdomain_id();
    if (neighbor_block == elem->subdomain_id())
      return q;

    auto q_neighbor = rule_fn(qrules(dim, neighbor_block));
    if (q->get_order() > q_neighbor->get_order())
      return q;
    return q_neighbor;
  }

  inline QRules & qrules(unsigned int dim) { return qrules(dim, _current_subdomain_id); }

  inline QRules & qrules(unsigned int dim, SubdomainID block)
  {
    if (_qrules.find(block) == _qrules.end())
    {
      mooseAssert(_qrules.find(Moose::ANY_BLOCK_ID) != _qrules.end(),
                  "missing quadrature rules for specified block");
      mooseAssert(_qrules[Moose::ANY_BLOCK_ID].size() > dim,
                  "quadrature rules not sized property for dimension");
      return _qrules[Moose::ANY_BLOCK_ID][dim];
    }
    mooseAssert(_qrules.find(block) != _qrules.end(),
                "missing quadrature rules for specified block");
    mooseAssert(_qrules[block].size() > dim, "quadrature rules not sized property for dimension");
    return _qrules[block][dim];
  }

  /**** Face Stuff ****/

  mutable std::map<unsigned int, std::map<FEType, FEBase *>> _fe_face;
  mutable std::map<unsigned int, std::map<FEType, FEVectorBase *>> _vector_fe_face;
  std::map<unsigned int, FEBase *> _holder_fe_face_helper;
  FEBase * _current_fe_face_helper;
  libMesh::QBase * _current_qrule_face;
  ArbitraryQuadrature * _current_qface_arbitrary;
  MooseArray<Point> _current_q_points_face;
  MooseArray<Real> _current_JxW_face;
  MooseArray<Point> _current_normals;
  std::vector<Eigen::Map<RealDIMValue>> _mapped_normals;
  MooseArray<std::vector<Point>> _current_tangents;

  std::vector<dof_id_type> _extra_elem_ids;
  std::vector<dof_id_type> _neighbor_extra_elem_ids;
  std::map<unsigned int, const std::vector<Point> *> _holder_normals;

  /**** Neighbor Stuff ****/

  mutable std::map<unsigned int, std::map<FEType, FEBase *>> _fe_neighbor;
  mutable std::map<unsigned int, std::map<FEType, FEBase *>> _fe_face_neighbor;
  mutable std::map<unsigned int, std::map<FEType, FEVectorBase *>> _vector_fe_neighbor;
  mutable std::map<unsigned int, std::map<FEType, FEVectorBase *>> _vector_fe_face_neighbor;

  std::map<unsigned int, FEBase *> _holder_fe_neighbor_helper;
  std::map<unsigned int, FEBase *> _holder_fe_face_neighbor_helper;

  mutable std::map<unsigned int, std::map<FEType, FEBase *>> _fe_lower;
  mutable std::map<unsigned int, std::map<FEType, FEVectorBase *>> _vector_fe_lower;
  std::map<unsigned int, FEBase *> _holder_fe_lower_helper;

  libMesh::QBase * _current_qrule_neighbor;
  MooseArray<Point> _current_q_points_face_neighbor;
  mutable bool _need_JxW_neighbor;
  MooseArray<Real> _current_JxW_neighbor;
  MooseArray<Real> _coord_neighbor;
  MooseArray<Real> _coord_msm;

  /********** mortar stuff *************/

  const std::vector<Real> * _JxW_msm;
  std::unique_ptr<FEBase> _fe_msm;
  libMesh::QBase * _qrule_msm;
  bool _custom_mortar_qrule;

  libMesh::QBase * _current_qrule_lower;

protected:
  const Elem * _current_elem;
  SubdomainID _current_subdomain_id;
  BoundaryID _current_boundary_id;
  Real _current_elem_volume;
  unsigned int _current_side;
  const Elem * _current_side_elem;
  Real _current_side_volume;
  const Elem * _current_neighbor_elem;
  SubdomainID _current_neighbor_subdomain_id;
  unsigned int _current_neighbor_side;
  const Elem * _current_neighbor_side_elem;
  mutable bool _need_neighbor_elem_volume;
  Real _current_neighbor_volume;
  const Node * _current_node;
  const Node * _current_neighbor_node;
  bool _current_elem_volume_computed;
  bool _current_side_volume_computed;

  const Elem * _current_lower_d_elem;
  const Elem * _current_neighbor_lower_d_elem;
  mutable bool _need_lower_d_elem_volume;
  Real _current_lower_d_elem_volume;
  mutable bool _need_neighbor_lower_d_elem_volume;
  Real _current_neighbor_lower_d_elem_volume;
  bool _need_dual;

  MooseArray<Point> _current_physical_points;

  /*
   * Residual contributions <tag_index, ivar>
   *
   * tag_index is the index into _residual_vector_tags, that is, _sub_Re[0] corresponds to the tag
   * with TagID _residual_vector_tags[0]._id
   *
   * When ivar corresponds to an array variable, the dense vector is in size of ndof * count,
   * where count is the number of components of the array variable. The local residual is ordered
   * as (r_i,j, i = 1,...,ndof; j = 1,...,count).
   *
   * Dense vectors for variables (ivar+i, i = 1,...,count) are empty.
   */
  std::vector<std::vector<DenseVector<Number>>> _sub_Re;
  std::vector<std::vector<DenseVector<Number>>> _sub_Rn;
  std::vector<std::vector<DenseVector<Number>>> _sub_Rl;

  DenseVector<Number> _tmp_Re;

  /*
   * Jacobian contributions <Tag, ivar, jvar>
   * When ivar corresponds to an array variable, the number of rows of the dense matrix is in size
   * of indof * icount, where icount is the number of components of ivar. When jvar corresponds to
   * an array variable, the number of columns of the dense matrix is in size of jndof * jcount,
   * where jcount is the number of components of jvar. The local residual is ordered as
   * (K_(i,j,k,l), k=1,...,jndof; l = 1,...,jcout; i = 1,...,indof; j = 1,...,icount).
   *
   * Dense matrices for variables (ivar+i, i = 1,...,icount) or (jvar+j, j = 1,...,jcount) are
   * empty.
   */
  std::vector<std::vector<std::vector<DenseMatrix<Number>>>> _sub_Kee;
  std::vector<std::vector<std::vector<DenseMatrix<Number>>>> _sub_Keg;

  std::vector<std::vector<std::vector<DenseMatrix<Number>>>> _sub_Ken;
  std::vector<std::vector<std::vector<DenseMatrix<Number>>>> _sub_Kne;
  std::vector<std::vector<std::vector<DenseMatrix<Number>>>> _sub_Knn;
  std::vector<std::vector<std::vector<DenseMatrix<Number>>>> _sub_Kll;
  std::vector<std::vector<std::vector<DenseMatrix<Number>>>> _sub_Kle;
  std::vector<std::vector<std::vector<DenseMatrix<Number>>>> _sub_Kln;
  std::vector<std::vector<std::vector<DenseMatrix<Number>>>> _sub_Kel;
  std::vector<std::vector<std::vector<DenseMatrix<Number>>>> _sub_Knl;

  DenseMatrix<Number> _tmp_Ke;

  // Shape function values, gradients. second derivatives
  VariablePhiValue _phi;
  VariablePhiGradient _grad_phi;
  VariablePhiSecond _second_phi;

  VariablePhiValue _phi_face;
  VariablePhiGradient _grad_phi_face;
  VariablePhiSecond _second_phi_face;

  VariablePhiValue _phi_neighbor;
  VariablePhiGradient _grad_phi_neighbor;
  VariablePhiSecond _second_phi_neighbor;

  VariablePhiValue _phi_face_neighbor;
  VariablePhiGradient _grad_phi_face_neighbor;
  VariablePhiSecond _second_phi_face_neighbor;

  // Shape function values, gradients, second derivatives
  VectorVariablePhiValue _vector_phi;
  VectorVariablePhiGradient _vector_grad_phi;
  VectorVariablePhiSecond _vector_second_phi;
  VectorVariablePhiCurl _vector_curl_phi;
  VectorVariablePhiDivergence _vector_div_phi;

  VectorVariablePhiValue _vector_phi_face;
  VectorVariablePhiGradient _vector_grad_phi_face;
  VectorVariablePhiSecond _vector_second_phi_face;
  VectorVariablePhiCurl _vector_curl_phi_face;
  VectorVariablePhiDivergence _vector_div_phi_face;

  VectorVariablePhiValue _vector_phi_neighbor;
  VectorVariablePhiGradient _vector_grad_phi_neighbor;
  VectorVariablePhiSecond _vector_second_phi_neighbor;
  VectorVariablePhiCurl _vector_curl_phi_neighbor;
  VectorVariablePhiDivergence _vector_div_phi_neighbor;

  VectorVariablePhiValue _vector_phi_face_neighbor;
  VectorVariablePhiGradient _vector_grad_phi_face_neighbor;
  VectorVariablePhiSecond _vector_second_phi_face_neighbor;
  VectorVariablePhiCurl _vector_curl_phi_face_neighbor;
  VectorVariablePhiDivergence _vector_div_phi_face_neighbor;

  class FEShapeData
  {
  public:
    VariablePhiValue _phi;
    VariablePhiGradient _grad_phi;
    VariablePhiSecond _second_phi;
    VariablePhiCurl _curl_phi;
    VariablePhiDivergence _div_phi;
  };

  class VectorFEShapeData
  {
  public:
    VectorVariablePhiValue _phi;
    VectorVariablePhiGradient _grad_phi;
    VectorVariablePhiSecond _second_phi;
    VectorVariablePhiCurl _curl_phi;
    VectorVariablePhiDivergence _div_phi;
  };

  mutable std::map<FEType, std::unique_ptr<FEShapeData>> _fe_shape_data;
  mutable std::map<FEType, std::unique_ptr<FEShapeData>> _fe_shape_data_face;
  mutable std::map<FEType, std::unique_ptr<FEShapeData>> _fe_shape_data_neighbor;
  mutable std::map<FEType, std::unique_ptr<FEShapeData>> _fe_shape_data_face_neighbor;
  mutable std::map<FEType, std::unique_ptr<FEShapeData>> _fe_shape_data_lower;
  mutable std::map<FEType, std::unique_ptr<FEShapeData>> _fe_shape_data_dual_lower;

  mutable std::map<FEType, std::unique_ptr<VectorFEShapeData>> _vector_fe_shape_data;
  mutable std::map<FEType, std::unique_ptr<VectorFEShapeData>> _vector_fe_shape_data_face;
  mutable std::map<FEType, std::unique_ptr<VectorFEShapeData>> _vector_fe_shape_data_neighbor;
  mutable std::map<FEType, std::unique_ptr<VectorFEShapeData>> _vector_fe_shape_data_face_neighbor;
  mutable std::map<FEType, std::unique_ptr<VectorFEShapeData>> _vector_fe_shape_data_lower;
  mutable std::map<FEType, std::unique_ptr<VectorFEShapeData>> _vector_fe_shape_data_dual_lower;

  mutable std::map<FEType, ADTemplateVariablePhiGradient<Real>> _ad_grad_phi_data;
  mutable std::map<FEType, ADTemplateVariablePhiGradient<RealVectorValue>> _ad_vector_grad_phi_data;
  mutable std::map<FEType, ADTemplateVariablePhiGradient<Real>> _ad_grad_phi_data_face;
  mutable std::map<FEType, ADTemplateVariablePhiGradient<RealVectorValue>>
      _ad_vector_grad_phi_data_face;

  const std::vector<VectorTag> & _residual_vector_tags;

  std::vector<std::vector<Real>> _cached_residual_values;

  std::vector<std::vector<dof_id_type>> _cached_residual_rows;

  unsigned int _max_cached_residuals;

  std::vector<std::vector<Real>> _cached_jacobian_values;
  std::vector<std::vector<dof_id_type>> _cached_jacobian_rows;
  std::vector<std::vector<dof_id_type>> _cached_jacobian_cols;

  unsigned int _max_cached_jacobians;

  bool _block_diagonal_matrix;
  std::vector<bool> _component_block_diagonal;

  std::vector<dof_id_type> _temp_dof_indices;

  std::vector<Point> _temp_reference_points;

  std::vector<VectorValue<ADReal>> _ad_dxyzdxi_map;
  std::vector<VectorValue<ADReal>> _ad_dxyzdeta_map;
  std::vector<VectorValue<ADReal>> _ad_dxyzdzeta_map;
  std::vector<VectorValue<ADReal>> _ad_d2xyzdxi2_map;
  std::vector<VectorValue<ADReal>> _ad_d2xyzdxideta_map;
  std::vector<VectorValue<ADReal>> _ad_d2xyzdeta2_map;
  std::vector<ADReal> _ad_jac;
  MooseArray<ADReal> _ad_JxW;
  MooseArray<VectorValue<ADReal>> _ad_q_points;
  std::vector<ADReal> _ad_dxidx_map;
  std::vector<ADReal> _ad_dxidy_map;
  std::vector<ADReal> _ad_dxidz_map;
  std::vector<ADReal> _ad_detadx_map;
  std::vector<ADReal> _ad_detady_map;
  std::vector<ADReal> _ad_detadz_map;
  std::vector<ADReal> _ad_dzetadx_map;
  std::vector<ADReal> _ad_dzetady_map;
  std::vector<ADReal> _ad_dzetadz_map;

  MooseArray<ADReal> _ad_JxW_face;
  MooseArray<VectorValue<ADReal>> _ad_normals;
  MooseArray<VectorValue<ADReal>> _ad_q_points_face;
  MooseArray<Real> _curvatures;
  MooseArray<ADReal> _ad_curvatures;

  std::vector<std::pair<unsigned int, unsigned short>> _disp_numbers_and_directions;

  mutable bool _calculate_xyz;
  mutable bool _calculate_face_xyz;
  mutable bool _calculate_curvatures;

  mutable bool _calculate_ad_coord;

  mutable std::set<FEType> _need_second_derivative;
  mutable std::set<FEType> _need_second_derivative_neighbor;
  mutable std::set<FEType> _need_curl;
  mutable std::set<FEType> _need_div;
  mutable std::set<FEType> _need_face_div;
  mutable std::set<FEType> _need_neighbor_div;
  mutable std::set<FEType> _need_face_neighbor_div;

  const NumericVector<Real> * _scaling_vector = nullptr;

  libMesh::ElemSideBuilder _current_side_elem_builder;
  libMesh::ElemSideBuilder _current_neighbor_side_elem_builder;
  libMesh::ElemSideBuilder _compute_face_map_side_elem_builder;

  const Elem * _msm_elem = nullptr;

  DenseVector<Number> _element_vector;

  DenseMatrix<Number> _element_matrix;

  std::vector<dof_id_type> _row_indices, _column_indices;

  bool _have_p_refinement;

  std::vector<Point> _current_neighbor_ref_points;
};

template <typename OutputType>
const typename OutputTools<OutputType>::VariablePhiValue &
Assembly::fePhiLower(FEType type) const
{
  buildLowerDFE(type);
  return _fe_shape_data_lower[type]->_phi;
}

template <typename OutputType>
const typename OutputTools<OutputType>::VariablePhiValue &
Assembly::feDualPhiLower(FEType type) const
{
  buildLowerDDualFE(type);
  return _fe_shape_data_dual_lower[type]->_phi;
}

template <typename OutputType>
const typename OutputTools<OutputType>::VariablePhiGradient &
Assembly::feGradPhiLower(FEType type) const
{
  buildLowerDFE(type);
  return _fe_shape_data_lower[type]->_grad_phi;
}

template <typename OutputType>
const typename OutputTools<OutputType>::VariablePhiGradient &
Assembly::feGradDualPhiLower(FEType type) const
{
  buildLowerDDualFE(type);
  return _fe_shape_data_dual_lower[type]->_grad_phi;
}

template <>
inline const ADTemplateVariablePhiGradient<RealVectorValue> &
Assembly::feADGradPhi<RealVectorValue>(FEType type) const
{
  return _ad_vector_grad_phi_data[type];
}

template <>
inline const ADTemplateVariablePhiGradient<RealVectorValue> &
Assembly::feADGradPhiFace<RealVectorValue>(FEType type) const
{
  return _ad_vector_grad_phi_data_face[type];
}

template <>
const typename OutputTools<VectorValue<Real>>::VariablePhiValue &
Assembly::fePhi<VectorValue<Real>>(FEType type) const;

template <>
const typename OutputTools<VectorValue<Real>>::VariablePhiGradient &
Assembly::feGradPhi<VectorValue<Real>>(FEType type) const;

template <>
const typename OutputTools<VectorValue<Real>>::VariablePhiSecond &
Assembly::feSecondPhi<VectorValue<Real>>(FEType type) const;

template <>
const typename OutputTools<VectorValue<Real>>::VariablePhiValue &
Assembly::fePhiLower<VectorValue<Real>>(FEType type) const;

template <>
const typename OutputTools<VectorValue<Real>>::VariablePhiValue &
Assembly::feDualPhiLower<VectorValue<Real>>(FEType type) const;

template <>
const typename OutputTools<VectorValue<Real>>::VariablePhiGradient &
Assembly::feGradPhiLower<VectorValue<Real>>(FEType type) const;

template <>
const typename OutputTools<VectorValue<Real>>::VariablePhiGradient &
Assembly::feGradDualPhiLower<VectorValue<Real>>(FEType type) const;

template <>
const typename OutputTools<VectorValue<Real>>::VariablePhiValue &
Assembly::fePhiFace<VectorValue<Real>>(FEType type) const;

template <>
const typename OutputTools<VectorValue<Real>>::VariablePhiGradient &
Assembly::feGradPhiFace<VectorValue<Real>>(FEType type) const;

template <>
const typename OutputTools<VectorValue<Real>>::VariablePhiSecond &
Assembly::feSecondPhiFace<VectorValue<Real>>(FEType type) const;

template <>
const typename OutputTools<VectorValue<Real>>::VariablePhiValue &
Assembly::fePhiNeighbor<VectorValue<Real>>(FEType type) const;

template <>
const typename OutputTools<VectorValue<Real>>::VariablePhiGradient &
Assembly::feGradPhiNeighbor<VectorValue<Real>>(FEType type) const;

template <>
const typename OutputTools<VectorValue<Real>>::VariablePhiSecond &
Assembly::feSecondPhiNeighbor<VectorValue<Real>>(FEType type) const;

template <>
const typename OutputTools<VectorValue<Real>>::VariablePhiValue &
Assembly::fePhiFaceNeighbor<VectorValue<Real>>(FEType type) const;

template <>
const typename OutputTools<VectorValue<Real>>::VariablePhiGradient &
Assembly::feGradPhiFaceNeighbor<VectorValue<Real>>(FEType type) const;

template <>
const typename OutputTools<VectorValue<Real>>::VariablePhiSecond &
Assembly::feSecondPhiFaceNeighbor<VectorValue<Real>>(FEType type) const;

template <>
const typename OutputTools<VectorValue<Real>>::VariablePhiCurl &
Assembly::feCurlPhi<VectorValue<Real>>(FEType type) const;

template <>
const typename OutputTools<VectorValue<Real>>::VariablePhiCurl &
Assembly::feCurlPhiFace<VectorValue<Real>>(FEType type) const;

template <>
const typename OutputTools<VectorValue<Real>>::VariablePhiCurl &
Assembly::feCurlPhiNeighbor<VectorValue<Real>>(FEType type) const;

template <>
const typename OutputTools<VectorValue<Real>>::VariablePhiCurl &
Assembly::feCurlPhiFaceNeighbor<VectorValue<Real>>(FEType type) const;

template <>
const typename OutputTools<VectorValue<Real>>::VariablePhiDivergence &
Assembly::feDivPhi<VectorValue<Real>>(FEType type) const;

template <>
const typename OutputTools<VectorValue<Real>>::VariablePhiDivergence &
Assembly::feDivPhiFace<VectorValue<Real>>(FEType type) const;

template <>
const typename OutputTools<VectorValue<Real>>::VariablePhiDivergence &
Assembly::feDivPhiNeighbor<VectorValue<Real>>(FEType type) const;

template <>
const typename OutputTools<VectorValue<Real>>::VariablePhiDivergence &
Assembly::feDivPhiFaceNeighbor<VectorValue<Real>>(FEType type) const;

template <>
inline const ADTemplateVariablePhiGradient<RealVectorValue> &
Assembly::adGradPhi<RealVectorValue>(const MooseVariableFE<RealVectorValue> & v) const
{
  return _ad_vector_grad_phi_data.at(v.feType());
}

template <typename Residuals, typename Indices>
void
Assembly::cacheResiduals(const Residuals & residuals,
                         const Indices & input_row_indices,
                         const Real scaling_factor,
                         LocalDataKey,
                         const std::set<TagID> & vector_tags)
{
  mooseAssert(residuals.size() == input_row_indices.size(),
              "The number of residuals should match the number of dof indices");
  mooseAssert(residuals.size() >= 1, "Why you calling me with no residuals?");

  if (!computingResidual() || vector_tags.empty())
    return;

  if (residuals.size() == 1)
  {
    // No constraining is required. (This is likely a finite volume computation if we only have a
    // single dof)
    cacheResidualsWithoutConstraints(
        residuals, input_row_indices, scaling_factor, LocalDataKey{}, vector_tags);
    return;
  }

  // Need to make a copy because we might modify this in constrain_element_vector
  _row_indices.assign(input_row_indices.begin(), input_row_indices.end());

  _element_vector.resize(_row_indices.size());
  for (const auto i : index_range(_row_indices))
    _element_vector(i) = MetaPhysicL::raw_value(residuals[i]) * scaling_factor;

  // At time of writing, this method doesn't do anything with the asymmetric_constraint_rows
  // argument, but we set it to false to be consistent with processLocalResidual
  _dof_map.constrain_element_vector(
      _element_vector, _row_indices, /*asymmetric_constraint_rows=*/false);

  for (const auto i : index_range(_row_indices))
    cacheResidual(_row_indices[i], _element_vector(i), vector_tags);
}

template <typename Residuals, typename Indices>
void
Assembly::cacheResidualsWithoutConstraints(const Residuals & residuals,
                                           const Indices & row_indices,
                                           const Real scaling_factor,
                                           LocalDataKey,
                                           const std::set<TagID> & vector_tags)
{
  mooseAssert(residuals.size() == row_indices.size(),
              "The number of residuals should match the number of dof indices");
  mooseAssert(residuals.size() >= 1, "Why you calling me with no residuals?");

  if (computingResidual() && !vector_tags.empty())
    for (const auto i : index_range(row_indices))
      cacheResidual(
          row_indices[i], MetaPhysicL::raw_value(residuals[i]) * scaling_factor, vector_tags);
}

template <typename Residuals, typename Indices>
void
Assembly::cacheJacobian(const Residuals & residuals,
                        const Indices & input_row_indices,
                        const Real scaling_factor,
                        LocalDataKey,
                        const std::set<TagID> & matrix_tags)
{
  if (!computingJacobian() || matrix_tags.empty())
    return;

  if (residuals.size() == 1)
  {
    // No constraining is required. (This is likely a finite volume computation if we only have a
    // single dof)
    cacheJacobianWithoutConstraints(
        residuals, input_row_indices, scaling_factor, LocalDataKey{}, matrix_tags);
    return;
  }

  const auto & compare_dofs = residuals[0].derivatives().nude_indices();
#ifndef NDEBUG
  auto compare_dofs_set = std::set<dof_id_type>(compare_dofs.begin(), compare_dofs.end());

  for (const auto i : make_range(decltype(residuals.size())(1), residuals.size()))
  {
    const auto & residual = residuals[i];
    auto current_dofs_set = std::set<dof_id_type>(residual.derivatives().nude_indices().begin(),
                                                  residual.derivatives().nude_indices().end());
    mooseAssert(compare_dofs_set == current_dofs_set,
                "We're going to see whether the dof sets are the same. IIRC the degree of freedom "
                "dependence (as indicated by the dof index set held by the ADReal) has to be the "
                "same for every residual passed to this method otherwise constrain_element_matrix "
                "will not work.");
  }
#endif
  _column_indices.assign(compare_dofs.begin(), compare_dofs.end());

  // If there's no derivatives then there is nothing to do. Moreover, if we pass zero size column
  // indices to constrain_element_matrix then we will potentially get errors out of BLAS
  if (!_column_indices.size())
    return;

  // Need to make a copy because we might modify this in constrain_element_matrix
  _row_indices.assign(input_row_indices.begin(), input_row_indices.end());

  _element_matrix.resize(_row_indices.size(), _column_indices.size());
  for (const auto i : index_range(_row_indices))
  {
    const auto & sparse_derivatives = residuals[i].derivatives();

    for (const auto j : index_range(_column_indices))
      _element_matrix(i, j) = sparse_derivatives[_column_indices[j]] * scaling_factor;
  }

  _dof_map.constrain_element_matrix(_element_matrix, _row_indices, _column_indices);

  for (const auto i : index_range(_row_indices))
    for (const auto j : index_range(_column_indices))
      cacheJacobian(_row_indices[i], _column_indices[j], _element_matrix(i, j), {}, matrix_tags);
}

template <typename Residuals, typename Indices>
void
Assembly::cacheJacobianWithoutConstraints(const Residuals & residuals,
                                          const Indices & row_indices,
                                          const Real scaling_factor,
                                          LocalDataKey,
                                          const std::set<TagID> & matrix_tags)
{
  mooseAssert(residuals.size() == row_indices.size(),
              "The number of residuals should match the number of dof indices");
  mooseAssert(residuals.size() >= 1, "Why you calling me with no residuals?");

  if (!computingJacobian() || matrix_tags.empty())
    return;

  for (const auto i : index_range(row_indices))
  {
    const auto row_index = row_indices[i];

    const auto & sparse_derivatives = residuals[i].derivatives();
    const auto & column_indices = sparse_derivatives.nude_indices();
    const auto & raw_derivatives = sparse_derivatives.nude_data();

    for (std::size_t j = 0; j < column_indices.size(); ++j)
      cacheJacobian(
          row_index, column_indices[j], raw_derivatives[j] * scaling_factor, {}, matrix_tags);
  }
}

inline const Real &
Assembly::lowerDElemVolume() const
{
  _need_lower_d_elem_volume = true;
  return _current_lower_d_elem_volume;
}

inline const Real &
Assembly::neighborLowerDElemVolume() const
{
  _need_neighbor_lower_d_elem_volume = true;
  return _current_neighbor_lower_d_elem_volume;
}

inline void
Assembly::assignDisplacements(
    std::vector<std::pair<unsigned int, unsigned short>> && disp_numbers_and_directions)
{
  _disp_numbers_and_directions = std::move(disp_numbers_and_directions);
}

inline void
Assembly::setCurrentLowerDElem(const Elem * const lower_d_elem)
{
  _current_lower_d_elem = lower_d_elem;
}