NSPressureNeumannBC::NSPressureNeumannBC(const InputParameters & parameters) :
    NSIntegratedBC(parameters),
    _pressure(coupledValue(NS::pressure)),
    _component(getParam<unsigned>("component")),
    _pressure_derivs(*this)
{
}

GBEvolutionBase::GBEvolutionBase(const InputParameters & parameters)
  : DerivativeMaterialInterface<Material>(parameters),
    _f0s(getParam<Real>("f0s")),
    _wGB(getParam<Real>("wGB")),
    _length_scale(getParam<Real>("length_scale")),
    _time_scale(getParam<Real>("time_scale")),
    _GBmob0(getParam<Real>("GBmob0")),
    _Q(getParam<Real>("Q")),
    _GBMobility(getParam<Real>("GBMobility")),
    _molar_vol(getParam<Real>("molar_volume")),
    _T(coupledValue("T")),
    _sigma(declareProperty<Real>("sigma")),
    _M_GB(declareProperty<Real>("M_GB")),
    _kappa(declareProperty<Real>("kappa_op")),
    _gamma(declareProperty<Real>("gamma_asymm")),
    _L(declareProperty<Real>("L")),
    _dLdT(parameters.hasDefaultCoupledValue("T")
              ? nullptr
              : &declarePropertyDerivative<Real>("L", getVar("T", 0)->name())),
    _l_GB(declareProperty<Real>("l_GB")),
    _mu(declareProperty<Real>("mu")),
    _entropy_diff(declareProperty<Real>("entropy_diff")),
    _molar_volume(declareProperty<Real>("molar_volume")),
    _act_wGB(declareProperty<Real>("act_wGB")),
    _kb(8.617343e-5),     // Boltzmann constant in eV/K
    _JtoeV(6.24150974e18) // Joule to eV conversion
{
  if (_GBMobility == -1 && _GBmob0 == 0)
    mooseError("Either a value for GBMobility or for GBmob0 and Q must be provided");
}

TensorMechanicsMaterial::TensorMechanicsMaterial(const std::string & name,
                                             InputParameters parameters)
    : Material(name, parameters),
      _grad_disp_x(coupledGradient("disp_x")),
      _grad_disp_y(coupledGradient("disp_y")),
      _grad_disp_z(_mesh.dimension() == 3 ? coupledGradient("disp_z") : _grad_zero),
      _grad_disp_x_old(_fe_problem.isTransient() ? coupledGradientOld("disp_x") : _grad_zero),
      _grad_disp_y_old(_fe_problem.isTransient() ? coupledGradientOld("disp_y") : _grad_zero),
      _grad_disp_z_old(_fe_problem.isTransient() && _mesh.dimension() == 3 ? coupledGradientOld("disp_z") : _grad_zero),
      _stress(declareProperty<RankTwoTensor>("stress")),
      _elastic_strain(declareProperty<RankTwoTensor>("elastic_strain")),
      _elasticity_tensor(declareProperty<ElasticityTensorR4>("elasticity_tensor")),
      _Jacobian_mult(declareProperty<ElasticityTensorR4>("Jacobian_mult")),
      //_d_stress_dT(declareProperty<RankTwoTensor>("d_stress_dT")),
      _euler_angle_1(getParam<Real>("euler_angle_1")),
      _euler_angle_2(getParam<Real>("euler_angle_2")),
      _euler_angle_3(getParam<Real>("euler_angle_3")),
      _Cijkl_vector(getParam<std::vector<Real> >("C_ijkl")),
      _all_21(getParam<bool>("all_21")),
      _Cijkl(),
      _Euler_angles(_euler_angle_1, _euler_angle_2, _euler_angle_3),
      _has_T(isCoupled("temperature")),
      _T(_has_T ? &coupledValue("temperature") : NULL)
{
  // fill in the local tensors from the input vector information

  _Cijkl.fillFromInputVector(_Cijkl_vector, _all_21);
}

PorousFlowAqueousPreDisMineral::PorousFlowAqueousPreDisMineral(const InputParameters & parameters)
  : PorousFlowMaterialVectorBase(parameters),
    _num_reactions(_dictator.numAqueousKinetic()),
    _sec_conc(_nodal_material
                  ? declareProperty<std::vector<Real>>("PorousFlow_mineral_concentration_nodal")
                  : declareProperty<std::vector<Real>>("PorousFlow_mineral_concentration_qp")),

    _porosity_old(_nodal_material ? getMaterialPropertyOld<Real>("PorousFlow_porosity_nodal")
                                  : getMaterialPropertyOld<Real>("PorousFlow_porosity_qp")),
    _sec_conc_old(
        _nodal_material
            ? getMaterialPropertyOld<std::vector<Real>>("PorousFlow_mineral_concentration_nodal")
            : getMaterialPropertyOld<std::vector<Real>>("PorousFlow_mineral_concentration_qp")),
    _reaction_rate(
        _nodal_material
            ? getMaterialProperty<std::vector<Real>>("PorousFlow_mineral_reaction_rate_nodal")
            : getMaterialProperty<std::vector<Real>>("PorousFlow_mineral_reaction_rate_qp")),

    _initial_conc_supplied(isParamValid("initial_concentrations")),
    _num_initial_conc(_initial_conc_supplied ? coupledComponents("initial_concentrations")
                                             : _num_reactions)
{
  if (_num_initial_conc != _dictator.numAqueousKinetic())
    mooseError("PorousFlowAqueousPreDisMineral: The number of initial concentrations is ",
               _num_initial_conc,
               " but the Dictator knows that the number of aqueous kinetic "
               "(precipitation-dissolution) reactions is ",
               _dictator.numAqueousKinetic());

  _initial_conc.resize(_num_initial_conc);
  if (_initial_conc_supplied)
    for (unsigned r = 0; r < _num_reactions; ++r)
      _initial_conc[r] = (_nodal_material ? &coupledNodalValue("initial_concentrations", r)
                                          : &coupledValue("initial_concentrations", r));
}

TensorMechanicsMaterial::TensorMechanicsMaterial(const std::string & name,
                                                 InputParameters parameters) :
    Material(name, parameters),
    _grad_disp_x(coupledGradient("disp_x")),
    _grad_disp_y(coupledGradient("disp_y")),
    _grad_disp_z(_mesh.dimension() == 3 ? coupledGradient("disp_z") : _grad_zero),
    _grad_disp_x_old(_fe_problem.isTransient() ? coupledGradientOld("disp_x") : _grad_zero),
    _grad_disp_y_old(_fe_problem.isTransient() ? coupledGradientOld("disp_y") : _grad_zero),
    _grad_disp_z_old(_fe_problem.isTransient() && _mesh.dimension() == 3 ? coupledGradientOld("disp_z") : _grad_zero),
    _stress(declareProperty<RankTwoTensor>("stress")),
    _total_strain(declareProperty<RankTwoTensor>("total_strain")),
    _elastic_strain(declareProperty<RankTwoTensor>("elastic_strain")),
    _elasticity_tensor(declareProperty<ElasticityTensorR4>("elasticity_tensor")),
    _Jacobian_mult(declareProperty<ElasticityTensorR4>("Jacobian_mult")),
    // _d_stress_dT(declareProperty<RankTwoTensor>("d_stress_dT")),
    _euler_angle_1(getParam<Real>("euler_angle_1")),
    _euler_angle_2(getParam<Real>("euler_angle_2")),
    _euler_angle_3(getParam<Real>("euler_angle_3")),
    _Cijkl_vector(getParam<std::vector<Real> >("C_ijkl")),
    _Cijkl(),
    _Euler_angles(_euler_angle_1, _euler_angle_2, _euler_angle_3),
    _has_T(isCoupled("temperature")),
    _T(_has_T ? &coupledValue("temperature") : NULL),
    _fill_method((RankFourTensor::FillMethod)(int)getParam<MooseEnum>("fill_method"))
{
  _Cijkl.fillFromInputVector(_Cijkl_vector, _fill_method);

  const std::vector<FunctionName> & fcn_names( getParam<std::vector<FunctionName> >("initial_stress") );
  const unsigned num = fcn_names.size();
  if (!(num == 0 || num == 3*3))
    mooseError("Either zero or " << 3*3 << " initial stress functions must be provided to TensorMechanicsMaterial.  You supplied " << num << "\n");
  _initial_stress.resize(num);
  for (unsigned i = 0 ; i < num ; ++i)
    _initial_stress[i] = &getFunctionByName(fcn_names[i]);
}

NSGravityPower::NSGravityPower(const InputParameters & parameters) :
    Kernel(parameters),
    _momentum_var(coupled("momentum")),
    _momentum(coupledValue("momentum")),
    _acceleration(getParam<Real>("acceleration"))
{
}

GrainAdvectionVelocity::GrainAdvectionVelocity(const InputParameters & parameters) :
   Material(parameters),
   _grain_data(getUserObject<ComputeGrainCenterUserObject>("grain_data")),
   _grain_volumes(_grain_data.getGrainVolumes()),
   _grain_centers(_grain_data.getGrainCenters()),
   _grain_force_torque(getUserObject<GrainForceAndTorqueInterface>("grain_force")),
   _grain_forces(_grain_force_torque.getForceValues()),
   _grain_torques(_grain_force_torque.getTorqueValues()),
   _grain_force_derivatives(_grain_force_torque.getForceDerivatives()),
   _grain_torque_derivatives(_grain_force_torque.getTorqueDerivatives()),
   _mt(getParam<Real>("translation_constant")),
   _mr(getParam<Real>("rotation_constant")),
   _ncrys(_grain_forces.size()),
   _vals(_ncrys),
   _grad_vals(_ncrys),
   _velocity_advection(declareProperty<std::vector<RealGradient> >("advection_velocity")),
   _div_velocity_advection(declareProperty<std::vector<Real> >("advection_velocity_divergence")),
   _velocity_advection_derivative_c(declareProperty<std::vector<RealGradient> >("advection_velocity_derivative_c")),
   _div_velocity_advection_derivative_c(declareProperty<std::vector<Real> >("advection_velocity_divergence_derivative_c")),
   _velocity_advection_derivative_eta(declareProperty<std::vector<RealGradient> >("advection_velocity_derivative_eta"))
{
  //Loop through grains and load coupled variables into the arrays
  for (unsigned int i = 0; i < _ncrys; ++i)
  {
    _vals[i] = &coupledValue("etas", i);
    _grad_vals[i] = &coupledGradient("etas",i);
  }
}

TotalFreeEnergyBase::TotalFreeEnergyBase(const InputParameters & parameters)
  : AuxKernel(parameters),
    _nvars(coupledComponents("interfacial_vars")),
    _vars(_nvars),
    _grad_vars(_nvars),
    _kappa_names(getParam<std::vector<MaterialPropertyName>>("kappa_names")),
    _nkappas(_kappa_names.size()),
    _additional_free_energy(coupledValue("additional_free_energy"))
{
  // Fetch coupled variables and their gradients
  for (unsigned int i = 0; i < _nvars; ++i)
  {
    _vars[i] = &coupledValue("interfacial_vars", i);
    _grad_vars[i] = &coupledGradient("interfacial_vars", i);
  }
}

MathFreeEnergy::MathFreeEnergy(const std::string & name,
                       InputParameters parameters) :
    DerivativeFunctionMaterialBase(name, parameters),
    _c(coupledValue("c")),
    _c_var(coupled("c"))
{
}

PorousFlowTemperature::PorousFlowTemperature(const InputParameters & parameters)
  : DerivativeMaterialInterface<PorousFlowMaterial>(parameters),

    _num_pf_vars(_dictator.numVariables()),
    _temperature_var(_nodal_material ? coupledNodalValue("temperature")
                                     : coupledValue("temperature")),
    _grad_temperature_var(_nodal_material ? nullptr : &coupledGradient("temperature")),
    _temperature_is_PF(_dictator.isPorousFlowVariable(coupled("temperature"))),
    _t_var_num(_temperature_is_PF ? _dictator.porousFlowVariableNum(coupled("temperature")) : 0),

    _temperature(_nodal_material ? declareProperty<Real>("PorousFlow_temperature_nodal")
                                 : declareProperty<Real>("PorousFlow_temperature_qp")),
    _dtemperature_dvar(
        _nodal_material ? declareProperty<std::vector<Real>>("dPorousFlow_temperature_nodal_dvar")
                        : declareProperty<std::vector<Real>>("dPorousFlow_temperature_qp_dvar")),
    _grad_temperature(_nodal_material
                          ? nullptr
                          : &declareProperty<RealGradient>("PorousFlow_grad_temperature_qp")),
    _dgrad_temperature_dgradv(_nodal_material ? nullptr
                                              : &declareProperty<std::vector<Real>>(
                                                    "dPorousFlow_grad_temperature_qp_dgradvar")),
    _dgrad_temperature_dv(_nodal_material ? nullptr
                                          : &declareProperty<std::vector<RealGradient>>(
                                                "dPorousFlow_grad_temperature_qp_dvar"))
{
}

GBEvolution::GBEvolution(const std::string & name,
                 InputParameters parameters)
    :Material(name, parameters),
   _temp(getParam<Real>("temp")),
   //_cg(coupledValue("cg")),
   _f0s(getParam<Real>("f0s")),
   _wGB(getParam<Real>("wGB")),
   _length_scale(getParam<Real>("length_scale")),
   _time_scale(getParam<Real>("time_scale")),
   _GBmob0(getParam<Real>("GBmob0")),
   _Q(getParam<Real>("Q")),
   _GBenergy(getParam<Real>("GBenergy")),
   _has_T(isCoupled("T")),
   _GBMobility(getParam<Real>("GBMobility")),
   _molar_vol(getParam<Real>("molar_volume")),
   _T(_has_T ? &coupledValue("T") : NULL),
   _sigma(declareProperty<Real>("sigma")),
   _M_GB(declareProperty<Real>("M_GB")),
   _kappa(declareProperty<Real>("kappa_op")),
   _gamma(declareProperty<Real>("gamma_asymm")),
   _L(declareProperty<Real>("L")),
   _l_GB(declareProperty<Real>("l_GB")),
   _mu(declareProperty<Real>("mu")),
   _entropy_diff(declareProperty<Real>("entropy_diff")),
   _molar_volume(declareProperty<Real>("molar_volume")),
   _act_wGB(declareProperty<Real>("act_wGB")),
   _tgrad_corr_mult(declareProperty<Real>("tgrad_corr_mult")),
   _kb(8.617343e-5) //Boltzmann constant in eV/K
{
  //Moose::out << "GB mob = " << _GBMobility << ", GBmob0 = " << _GBmob0 << std::endl;

  if (_GBMobility == -1 && _GBmob0 == 0)
    mooseError("Either a value for GBMobility or for GBmob0 and Q must be provided");
}

PressureDarcy::PressureDarcy(const InputParameters & parameters)
  :Diffusion(parameters),
   _permeability(getMaterialProperty<Real>("permeability")),
   _mu_h2o(getMaterialProperty<Real>("WaterViscosity")),
   _porosity(coupledValue("porosity"))
//   _HBO2(coupledValue("HBO2"))
{}

RheaEnergy::RheaEnergy(const std::string & name,
                       InputParameters parameters) :
  Kernel(name, parameters),
    // Material values:
    _vel(coupledValue("velocity")),
    _temp(isCoupled("temperature") ? coupledValue("temperature") : _zero),
    _pressure(coupledValue("pressure")),
    // Radiation value:
    _epsilon(isCoupled("radiation") ? coupledValue("radiation") : _zero),
    _grad_eps(isCoupled("radiation") ?  coupledGradient("radiation") : _grad_zero),
    // Material property:
    _sigma_a(getMaterialProperty<Real>("sigma_a")),
    // Constant:
    _c(getParam<Real>("speed_of_light")),
    _a(getParam<Real>("a"))
{}

SahoteMaterial::SahoteMaterial(const std::string & name,
                           InputParameters parameters)
    :Material(name, parameters),

     // Get simple parameters from the input fil
     _temperature(coupledValue("temperature")),
     // Parameters specifying material properties used for both structural/sample materials and cooling fluids
     _density(declareProperty<Real>("Density")),
     _specific_heat_capacity(declareProperty<Real>("SpecificHeatCapacity")),
     _thermal_conductivity(declareProperty<Real>("ThermalConductivity")),


     // Paramaters used for cooling fluids only
     _boiling_point(declareProperty<Real>("BoilingPoint")),
     _melting_point(declareProperty<Real>("MeltingPoint")),
     _dynamic_viscosity(declareProperty<Real>("DynamicViscosity")),
     _surface_tension(declareProperty<Real>("SurfaceTension")),
     _latent_heat_of_vapourisation(declareProperty<Real>("LatentHeatOfVapuorisation")),
     _material_type(getParam<MooseEnum>("material_type"))
     




{
}

MuMag::MuMag(const InputParameters & parameters) :
  AuxKernel(parameters),
   _Bx(coupledValue("Bx")),
   _By(coupledValue("By")),
   _Bz(coupledValue("Bz")),
   _A(getParam<Real>("A")),
   _B(getParam<Real>("B")),
   _C(getParam<Real>("C")),
   _D(getParam<Real>("D")),
   _E(getParam<Real>("E")),
   _mu(getParam<Real>("mu"))
{
  std::cout<<"First, calculating the field from a coil located at loc_x, loc_y, loc_z"<<"\n";
  std::cout<<"Next, calculating the field from  "<<"\n";
  std::cout<<"Note, we are using MOOSE's API and the 'solve' is pointless!"<<"\n";
}

NucleationLocationUserObject::NucleationLocationUserObject(const InputParameters & parameters) :
    ElementUserObject(parameters),
    _mesh(_subproblem.mesh()),
    //_coupled_probability(coupledValue("coupled_aux")),
    _n_coupled_aux(getParam<int>("n_coupled_aux")),
    _dwell_time(getParam<Real>("dwell_time")),
    _num_orientations(getParam<int>("num_orientations")),
    _boundary_width(getParam<Real>("boundary_width")),
    _random_seed(getParam<int>("random_seed")),
    // _counter(0),
    //make restartable
    _counter(declareRestartableData<int>("counter", 0)),
    _phase_gen_index(std::numeric_limits<unsigned int>::max()),
    // _nuclei(0),
    //make restartable
    _nuclei(declareRestartableData<std::vector<Nucleus> >("nuclei")),
    //_old_nucleus_list_size(0),
    //make restartable
    _old_nucleus_list_size(declareRestartableData<unsigned int>("old_nucleus_list_size", 0)),
    _has_new_nucleus(false),
    _master_random(-1000),
    _slave_random(-100)
{
    if(_n_coupled_aux != coupledComponents("coupled_aux_vars"))
    mooseError("Please specify the correct # of coupled probabilities (NucleationLocationUserObject).");

  _coupled_probability.resize(_n_coupled_aux);

  for(unsigned int i=0; i<_n_coupled_aux; i++)
    _coupled_probability[i] = &coupledValue("coupled_aux_vars", i);
}

MaterialTensorOnLine :: MaterialTensorOnLine(const std::string & name, InputParameters parameters) :
  ElementUserObject(name, parameters),

  _tensor( getMaterialProperty<SymmTensor>( getParam<std::string>("tensor") ) ),
  _index( getParam<int>("index") ),
  _quantity_moose_enum( getParam<MooseEnum>("quantity") ),
  _p1( getParam<RealVectorValue>("point1") ),
  _p2( getParam<RealVectorValue>("point2") ),
  _line_id( getParam<int>("line_id") ),
  _file_name( getParam<std::string>("filename") ),

  _stream_open(false),

  _elem_line_id(coupledValue("element_line_id"))

{

  MaterialTensorAux::checkMaterialTensorParams(_quantity,_quantity_moose_enum,_index,_name);

  if (!_stream_open && libMesh::processor_id() == 0)
  {
    _output_file.open(_file_name.c_str(), std::ios::trunc | std::ios::out);
    _stream_open = true;
  }

}

KKSSplitCHCRes::KKSSplitCHCRes(const InputParameters & parameters) :
    DerivativeMaterialInterface<JvarMapInterface<SplitCHBase> >(parameters),
    // number of coupled variables (ca, args_a[])
    _nvar(_coupled_moose_vars.size()),
    _ca_var(coupled("ca")),
    _ca_name(getVar("ca", 0)->name()),
    _cb_var(coupled("cb")),
    _cb_name(getVar("cb", 0)->name()),
    _prop_h(getMaterialProperty<Real>("h_name")),
    _first_derivative_Fa(getMaterialPropertyDerivative<Real>("fa_name", _ca_name)),
    _second_derivative_Fa(getMaterialPropertyDerivative<Real>("fa_name", _ca_name, _ca_name)),
    _second_derivative_Fb(getMaterialPropertyDerivative<Real>("fb_name", _cb_name, _cb_name)),
    _w_var(coupled("w")),
    _w(coupledValue("w"))
{
  // reserve space for derivatives
  _d2Fadcadarg.resize(_nvar);

  // Iterate over all coupled variables
  for (unsigned int i = 0; i < _nvar; ++i)
  {
    MooseVariable *cvar = this->_coupled_moose_vars[i];

    // get the second derivative material property
    _d2Fadcadarg[i] = &getMaterialPropertyDerivative<Real>("fa_name", _ca_name, cvar->name());
  }
}

PorousFlowVariableBase::PorousFlowVariableBase(const InputParameters & parameters) :
    DerivativeMaterialInterface<Material>(parameters),

    _dictator(getUserObject<PorousFlowDictator>("PorousFlowDictator")),
    _num_phases(_dictator.numPhases()),
    _num_components(_dictator.numComponents()),
    _num_pf_vars(_dictator.numVariables()),
    _temperature_nodal_var(coupledNodalValue("temperature")),
    _temperature_qp_var(coupledValue("temperature")),
    _temperature_varnum(coupled("temperature")),

    _porepressure_nodal(declareProperty<std::vector<Real> >("PorousFlow_porepressure_nodal")),
    _porepressure_nodal_old(declarePropertyOld<std::vector<Real> >("PorousFlow_porepressure_nodal")),
    _porepressure_qp(declareProperty<std::vector<Real> >("PorousFlow_porepressure_qp")),
    _gradp_qp(declareProperty<std::vector<RealGradient> >("PorousFlow_grad_porepressure_qp")),
    _dporepressure_nodal_dvar(declareProperty<std::vector<std::vector<Real> > >("dPorousFlow_porepressure_nodal_dvar")),
    _dporepressure_qp_dvar(declareProperty<std::vector<std::vector<Real> > >("dPorousFlow_porepressure_qp_dvar")),
    _dgradp_qp_dgradv(declareProperty<std::vector<std::vector<Real> > >("dPorousFlow_grad_porepressure_qp_dgradvar")),
    _dgradp_qp_dv(declareProperty<std::vector<std::vector<RealGradient> > >("dPorousFlow_grad_porepressure_qp_dvar")),

    _saturation_nodal(declareProperty<std::vector<Real> >("PorousFlow_saturation_nodal")),
    _saturation_nodal_old(declarePropertyOld<std::vector<Real> >("PorousFlow_saturation_nodal")),
    _saturation_qp(declareProperty<std::vector<Real> >("PorousFlow_saturation_qp")),
    _grads_qp(declareProperty<std::vector<RealGradient> >("PorousFlow_grad_saturation_qp")),
    _dsaturation_nodal_dvar(declareProperty<std::vector<std::vector<Real> > >("dPorousFlow_saturation_nodal_dvar")),
    _dsaturation_qp_dvar(declareProperty<std::vector<std::vector<Real> > >("dPorousFlow_saturation_qp_dvar")),
    _dgrads_qp_dgradv(declareProperty<std::vector<std::vector<Real> > >("dPorousFlow_grad_saturation_qp_dgradvar")),
    _dgrads_qp_dv(declareProperty<std::vector<std::vector<RealGradient> > >("dPorousFlow_grad_saturation_qp_dv")),

    _temperature_nodal(declareProperty<std::vector<Real> >("PorousFlow_temperature_nodal")),
    _temperature_qp(declareProperty<std::vector<Real> >("PorousFlow_temperature_qp")),
    _dtemperature_nodal_dvar(declareProperty<std::vector<std::vector<Real> > >("dPorousFlow_temperature_nodal_dvar")),
    _dtemperature_qp_dvar(declareProperty<std::vector<std::vector<Real> > >("dPorousFlow_temperature_qp_dvar"))
{
}

RichardsPorepressureNames::RichardsPorepressureNames(const std::string & name, InputParameters parameters) :
  GeneralUserObject(name, parameters),
  Coupleable(parameters, false),
  ZeroInterface(parameters),
  _num_p(coupledComponents("porepressure_vars")),
  _the_names(std::string())

{
  unsigned int max_moose_var_num_seen(0);

  _moose_var_num.resize(_num_p);
  _moose_var_value.resize(_num_p);
  _moose_var_value_old.resize(_num_p);
  _moose_grad_var.resize(_num_p);
  _moose_raw_var.resize(_num_p);
  for (unsigned int i=0 ; i<_num_p; ++i)
    {
      _moose_var_num[i] = coupled("porepressure_vars", i);
      max_moose_var_num_seen = (max_moose_var_num_seen > _moose_var_num[i] ? max_moose_var_num_seen : _moose_var_num[i]);
      _moose_var_value[i] = &coupledValue("porepressure_vars", i); // coupledValue returns a reference (an alias) to a VariableValue, and the & turns it into a pointer
      _moose_var_value_old[i] = (_is_transient ? &coupledValueOld("porepressure_vars", i) : &_zero);
      _moose_grad_var[i] = &coupledGradient("porepressure_vars", i);
      _moose_raw_var[i] = getVar("porepressure_vars", i);
      _the_names += getVar("porepressure_vars", i)->name() + " ";
    }
  _the_names.erase(_the_names.end() - 1, _the_names.end()); // remove trailing space

  _pressure_var_num.resize(max_moose_var_num_seen + 1);
  for (unsigned int i=0 ; i<max_moose_var_num_seen+1 ; ++i)
    _pressure_var_num[i] = _num_p; // NOTE: indicates that i is not a porepressure variable
  for (unsigned int i=0 ; i<_num_p; ++i)
    _pressure_var_num[_moose_var_num[i]] = i;
}