C++ DiffContext类代码示例

bool ElasticitySystem::mass_residual(bool request_jacobian,
                                     DiffContext & context)
{
  FEMContext & c = cast_ref<FEMContext &>(context);

  FEBase * u_elem_fe;
  c.get_element_fe(u_var, u_elem_fe);

  // The number of local degrees of freedom in each variable
  const unsigned int n_u_dofs = c.get_dof_indices(u_var).size();

  // Element Jacobian * quadrature weights for interior integration
  const std::vector<Real> & JxW = u_elem_fe->get_JxW();

  const std::vector<std::vector<Real> > & phi = u_elem_fe->get_phi();

  // Residuals that we're populating
  libMesh::DenseSubVector<libMesh::Number> & Fu = c.get_elem_residual(u_var);
  libMesh::DenseSubVector<libMesh::Number> & Fv = c.get_elem_residual(v_var);
  libMesh::DenseSubVector<libMesh::Number> & Fw = c.get_elem_residual(w_var);

  libMesh::DenseSubMatrix<libMesh::Number> & Kuu = c.get_elem_jacobian(u_var, u_var);
  libMesh::DenseSubMatrix<libMesh::Number> & Kvv = c.get_elem_jacobian(v_var, v_var);
  libMesh::DenseSubMatrix<libMesh::Number> & Kww = c.get_elem_jacobian(w_var, w_var);

  unsigned int n_qpoints = c.get_element_qrule().n_points();

  for (unsigned int qp=0; qp != n_qpoints; qp++)
    {
      libMesh::Number u_ddot, v_ddot, w_ddot;
      c.interior_accel(u_var, qp, u_ddot);
      c.interior_accel(v_var, qp, v_ddot);
      c.interior_accel(w_var, qp, w_ddot);

      for (unsigned int i=0; i != n_u_dofs; i++)
        {
          Fu(i) += _rho*u_ddot*phi[i][qp]*JxW[qp];
          Fv(i) += _rho*v_ddot*phi[i][qp]*JxW[qp];
          Fw(i) += _rho*w_ddot*phi[i][qp]*JxW[qp];

          if (request_jacobian)
            {
              for (unsigned int j=0; j != n_u_dofs; j++)
                {
                  libMesh::Real jac_term = _rho*phi[i][qp]*phi[j][qp]*JxW[qp];
                  jac_term *= context.get_elem_solution_accel_derivative();

                  Kuu(i,j) += jac_term;
                  Kvv(i,j) += jac_term;
                  Kww(i,j) += jac_term;
                }
            }
        }
    }

  // If the Jacobian was requested, we computed it. Otherwise, we didn't.
  return request_jacobian;
}

bool L2System::element_time_derivative (bool request_jacobian,
                                        DiffContext & context)
{
  FEMContext & c = cast_ref<FEMContext &>(context);

  // First we get some references to cell-specific data that
  // will be used to assemble the linear system.

  // Element Jacobian * quadrature weights for interior integration
  const std::vector<Real> & JxW = c.get_element_fe(0)->get_JxW();

  const std::vector<std::vector<Real>> & phi = c.get_element_fe(0)->get_phi();

  const std::vector<Point> & xyz = c.get_element_fe(0)->get_xyz();

  // The number of local degrees of freedom in each variable
  const unsigned int n_u_dofs = c.n_dof_indices(0);

  // The subvectors and submatrices we need to fill:
  DenseSubMatrix<Number> & K = c.get_elem_jacobian(0, 0);
  DenseSubVector<Number> & F = c.get_elem_residual(0);

  unsigned int n_qpoints = c.get_element_qrule().n_points();

  libmesh_assert (input_contexts.find(&c) != input_contexts.end());

  FEMContext & input_c = *input_contexts[&c];
  input_c.pre_fe_reinit(*input_system, &c.get_elem());
  input_c.elem_fe_reinit();

  for (unsigned int qp=0; qp != n_qpoints; qp++)
    {
      Number u = c.interior_value(0, qp);

      Number ufunc = (*goal_func)(input_c, xyz[qp]);

      for (unsigned int i=0; i != n_u_dofs; i++)
        F(i) += JxW[qp] * ((u - ufunc) * phi[i][qp]);
      if (request_jacobian)
        {
          const Number JxWxD = JxW[qp] *
            context.get_elem_solution_derivative();

          for (unsigned int i=0; i != n_u_dofs; i++)
            for (unsigned int j=0; j != n_u_dofs; ++j)
              K(i,j) += JxWxD * (phi[i][qp] * phi[j][qp]);
        }
    } // end of the quadrature point qp-loop

  return request_jacobian;
}

bool DifferentiablePhysics::nonlocal_mass_residual(bool request_jacobian,
                                                   DiffContext &c)
{
  FEMContext &context = cast_ref<FEMContext&>(c);

  for (unsigned int var = 0; var != context.n_vars(); ++var)
    {
      if (!this->is_time_evolving(var))
        continue;

      if (c.get_system().variable(var).type().family != SCALAR)
        continue;

      const std::vector<dof_id_type>& dof_indices =
        context.get_dof_indices(var);

      const unsigned int n_dofs = cast_int<unsigned int>
        (dof_indices.size());

      DenseSubVector<Number> &Fs = context.get_elem_residual(var);
      DenseSubMatrix<Number> &Kss = context.get_elem_jacobian( var, var );

      const libMesh::DenseSubVector<libMesh::Number> &Us =
        context.get_elem_solution(var);

      for (unsigned int i=0; i != n_dofs; ++i)
        {
          Fs(i) -= Us(i);

          if (request_jacobian)
            Kss(i,i) -= context.elem_solution_rate_derivative;
        }
    }

  return request_jacobian;
}

bool Euler2Solver::_general_residual (bool request_jacobian,
                                      DiffContext & context,
                                      ResFuncType mass,
                                      ResFuncType damping,
                                      ResFuncType time_deriv,
                                      ResFuncType constraint,
                                      ReinitFuncType reinit_func,
                                      bool compute_second_order_eqns)
{
  unsigned int n_dofs = context.get_elem_solution().size();

  // Local nonlinear solution at old timestep
  DenseVector<Number> old_elem_solution(n_dofs);
  for (unsigned int i=0; i != n_dofs; ++i)
    old_elem_solution(i) =
      old_nonlinear_solution(context.get_dof_indices()[i]);

  // Local time derivative of solution
  context.get_elem_solution_rate() = context.get_elem_solution();
  context.get_elem_solution_rate() -= old_elem_solution;
  context.elem_solution_rate_derivative = 1 / _system.deltat;
  context.get_elem_solution_rate() *=
    context.elem_solution_rate_derivative;

  // Our first evaluations are at the final elem_solution
  context.elem_solution_derivative = 1.0;

  // If a fixed solution is requested, we'll use the elem_solution
  // at the new timestep
  // FIXME - should this be the theta solution instead?
  if (_system.use_fixed_solution)
    context.get_elem_fixed_solution() = context.get_elem_solution();

  context.fixed_solution_derivative = 1.0;

  // We need to save the old jacobian and old residual since we'll be
  // multiplying some of the new contributions by theta or 1-theta
  DenseMatrix<Number> old_elem_jacobian(n_dofs, n_dofs);
  DenseVector<Number> old_elem_residual(n_dofs);
  old_elem_residual.swap(context.get_elem_residual());
  if (request_jacobian)
    old_elem_jacobian.swap(context.get_elem_jacobian());

  // Local time derivative of solution
  context.get_elem_solution_rate() = context.get_elem_solution();
  context.get_elem_solution_rate() -= old_elem_solution;
  context.elem_solution_rate_derivative = 1 / _system.deltat;
  context.get_elem_solution_rate() *=
    context.elem_solution_rate_derivative;

  // If we are asked to compute residuals for second order variables,
  // we also populate the acceleration part so the user can use that.
  if (compute_second_order_eqns)
    this->prepare_accel(context);

  // Move the mesh into place first if necessary, set t = t_{n+1}
  (context.*reinit_func)(1.);

  // First, evaluate time derivative at the new timestep.
  // The element should already be in the proper place
  // even for a moving mesh problem.
  bool jacobian_computed =
    (_system.get_physics()->*time_deriv)(request_jacobian, context);

  // Next, evaluate the mass residual at the new timestep

  jacobian_computed = (_system.get_physics()->*mass)(jacobian_computed, context) &&
    jacobian_computed;

  // If we have second-order variables, we need to get damping terms
  // and the velocity equations
  if (compute_second_order_eqns)
    {
      jacobian_computed = (_system.get_physics()->*damping)(jacobian_computed, context) &&
        jacobian_computed;

      jacobian_computed = this->compute_second_order_eqns(jacobian_computed, context) &&
        jacobian_computed;
    }

  // Add the constraint term
  jacobian_computed = (_system.get_physics()->*constraint)(jacobian_computed, context) &&
    jacobian_computed;

  // The new solution's contribution is scaled by theta
  context.get_elem_residual() *= theta;
  context.get_elem_jacobian() *= theta;

  // Save the new solution's term
  DenseMatrix<Number> elem_jacobian_newterm(n_dofs, n_dofs);
  DenseVector<Number> elem_residual_newterm(n_dofs);
  elem_residual_newterm.swap(context.get_elem_residual());
  if (request_jacobian)
    elem_jacobian_newterm.swap(context.get_elem_jacobian());

  // Add the time-dependent term for the old solution

  // Make sure elem_solution is set up for elem_reinit to use
  // Move elem_->old_, old_->elem_
  context.get_elem_solution().swap(old_elem_solution);
//.........这里部分代码省略.........

//for non-Dirichlet boundary conditions and the bit from diffusion term
bool ConvDiff_AuxSadjSys::side_time_derivative(bool request_jacobian, DiffContext & context)
{
    const unsigned int dim = this->get_mesh().mesh_dimension();

    FEMContext &ctxt = cast_ref<FEMContext&>(context);

    // First we get some references to cell-specific data that
    // will be used to assemble the linear system.
    FEBase* side_fe = NULL;
    ctxt.get_side_fe(aux_c_var, side_fe );

    // Element Jacobian * quadrature weights for interior integration
    const std::vector<Real> &JxW = side_fe->get_JxW();

    // Side basis functions
    const std::vector<std::vector<Real> > &phi = side_fe->get_phi();

    // Side Quadrature points
    const std::vector<Point > &qside_point = side_fe->get_xyz();

    //normal vector
    const std::vector<Point> &face_normals = side_fe->get_normals();

    // The number of local degrees of freedom in each variable
    const unsigned int n_c_dofs = ctxt.get_dof_indices(aux_c_var).size();

    // The subvectors and submatrices we need to fill:
    DenseSubMatrix<Number> &J_c_auxz = ctxt.get_elem_jacobian(aux_c_var, aux_zc_var);

    DenseSubVector<Number> &Rc = ctxt.get_elem_residual( aux_c_var );
    //Rf gets no contribution from sides

    unsigned int n_qpoints = ctxt.get_side_qrule().n_points();

    bool isEast = false;
    if (dim == 3) {
        isEast = ctxt.has_side_boundary_id(2);
    }
    else if (dim == 2) {
        isEast = ctxt.has_side_boundary_id(1);
    }

    //set (in)flux boundary condition on west side
    //homogeneous neumann (Danckwerts) outflow boundary condition on east side
    //no-flux (equivalently, homoegenous neumann) boundary conditions on north, south, top, bottom sides
    //"strong" enforcement of boundary conditions
    for (unsigned int qp=0; qp != n_qpoints; qp++)
    {
        Number aux_z = ctxt.side_value(aux_zc_var, qp);

        //velocity vector
        NumberVectorValue U(porosity*vx, 0., 0.);

        for (unsigned int i=0; i != n_c_dofs; i++)
        {
            if(isEast)
                Rc(i) += JxW[qp]*(-U*face_normals[qp]*aux_z)*phi[i][qp];

            if(request_jacobian && context.get_elem_solution_derivative())
            {
                for (unsigned int j=0; j != n_c_dofs; j++)
                {
                    if(isEast)
                        J_c_auxz(i,j) += JxW[qp]*(-U*face_normals[qp]*phi[j][qp])*phi[i][qp];
                }
            } // end - if (request_jacobian && context.get_elem_solution_derivative())
        } //end of outer dof (i) loop
    }

    return request_jacobian;
}

bool HeatSystem::element_time_derivative (bool request_jacobian,
                                          DiffContext & context)
{
  FEMContext & c = libmesh_cast_ref<FEMContext &>(context);

  const unsigned int mesh_dim =
    c.get_system().get_mesh().mesh_dimension();

  // First we get some references to cell-specific data that
  // will be used to assemble the linear system.
  const unsigned int dim = c.get_elem().dim();
  FEBase * fe = libmesh_nullptr;
  c.get_element_fe(T_var, fe, dim);

  // Element Jacobian * quadrature weights for interior integration
  const std::vector<Real> & JxW = fe->get_JxW();

  const std::vector<Point> & xyz = fe->get_xyz();

  const std::vector<std::vector<Real> > & phi = fe->get_phi();

  const std::vector<std::vector<RealGradient> > & dphi = fe->get_dphi();

  // The number of local degrees of freedom in each variable
  const unsigned int n_T_dofs = c.get_dof_indices(T_var).size();

  // The subvectors and submatrices we need to fill:
  DenseSubMatrix<Number> & K = c.get_elem_jacobian(T_var, T_var);
  DenseSubVector<Number> & F = c.get_elem_residual(T_var);

  // Now we will build the element Jacobian and residual.
  // Constructing the residual requires the solution and its
  // gradient from the previous timestep.  This must be
  // calculated at each quadrature point by summing the
  // solution degree-of-freedom values by the appropriate
  // weight functions.
  unsigned int n_qpoints = c.get_element_qrule().n_points();

  for (unsigned int qp=0; qp != n_qpoints; qp++)
    {
      // Compute the solution gradient at the Newton iterate
      Gradient grad_T = c.interior_gradient(T_var, qp);

      const Number k = _k[dim];

      const Point & p = xyz[qp];

      // solution + laplacian depend on problem dimension
      const Number u_exact = (mesh_dim == 2) ?
        std::sin(libMesh::pi*p(0)) * std::sin(libMesh::pi*p(1)) :
        std::sin(libMesh::pi*p(0)) * std::sin(libMesh::pi*p(1)) *
        std::sin(libMesh::pi*p(2));

      // Only apply forcing to interior elements
      const Number forcing = (dim == mesh_dim) ?
        -k * u_exact * (dim * libMesh::pi * libMesh::pi) : 0;

      const Number JxWxNK = JxW[qp] * -k;

      for (unsigned int i=0; i != n_T_dofs; i++)
        F(i) += JxWxNK * (grad_T * dphi[i][qp] + forcing * phi[i][qp]);
      if (request_jacobian)
        {
          const Number JxWxNKxD = JxWxNK *
            context.get_elem_solution_derivative();

          for (unsigned int i=0; i != n_T_dofs; i++)
            for (unsigned int j=0; j != n_T_dofs; ++j)
              K(i,j) += JxWxNKxD * (dphi[i][qp] * dphi[j][qp]);
        }
    } // end of the quadrature point qp-loop

  return request_jacobian;
}

bool NewmarkSolver::_general_residual (bool request_jacobian,
                                       DiffContext & context,
                                       ResFuncType mass,
                                       ResFuncType damping,
                                       ResFuncType time_deriv,
                                       ResFuncType constraint)
{
  unsigned int n_dofs = context.get_elem_solution().size();

  // We might need to save the old jacobian in case one of our physics
  // terms later is unable to update it analytically.
  DenseMatrix<Number> old_elem_jacobian(n_dofs, n_dofs);

  // Local velocity at old time step
  DenseVector<Number> old_elem_solution_rate(n_dofs);
  for (unsigned int i=0; i != n_dofs; ++i)
    old_elem_solution_rate(i) =
      old_solution_rate(context.get_dof_indices()[i]);

  // The user is computing the initial acceleration
  // So upstream we've swapped _system.solution and _old_local_solution_accel
  // So we need to give the context the correct entries since we're solving for
  // acceleration here.
  if( _is_accel_solve )
    {
      // System._solution is actually the acceleration right now so we need
      // to reset the elem_solution to the right thing, which in this case
      // is the initial guess for displacement, which got swapped into
      // _old_solution_accel vector
      DenseVector<Number> old_elem_solution(n_dofs);
      for (unsigned int i=0; i != n_dofs; ++i)
        old_elem_solution(i) =
          old_solution_accel(context.get_dof_indices()[i]);

      context.elem_solution_derivative = 0.0;
      context.elem_solution_rate_derivative = 0.0;
      context.elem_solution_accel_derivative = 1.0;

      // Acceleration is currently the unknown so it's already sitting
      // in elem_solution() thanks to FEMContext::pre_fe_reinit
      context.get_elem_solution_accel() = context.get_elem_solution();

      // Now reset elem_solution() to what the user is expecting
      context.get_elem_solution() = old_elem_solution;

      context.get_elem_solution_rate() = old_elem_solution_rate;

      // The user's Jacobians will be targeting derivatives w.r.t. u_{n+1}.
      // Although the vast majority of cases will have the correct analytic
      // Jacobians in this iteration, since we reset elem_solution_derivative*,
      // if there are coupled/overlapping problems, there could be
      // mismatches in the Jacobian. So we force finite differencing for
      // the first iteration.
      request_jacobian = false;
    }
  // Otherwise, the unknowns are the displacements and everything is straight
  // foward and is what you think it is
  else
    {
      if (request_jacobian)
        old_elem_jacobian.swap(context.get_elem_jacobian());

      // Local displacement at old timestep
      DenseVector<Number> old_elem_solution(n_dofs);
      for (unsigned int i=0; i != n_dofs; ++i)
        old_elem_solution(i) =
          old_nonlinear_solution(context.get_dof_indices()[i]);

      // Local acceleration at old time step
      DenseVector<Number> old_elem_solution_accel(n_dofs);
      for (unsigned int i=0; i != n_dofs; ++i)
        old_elem_solution_accel(i) =
          old_solution_accel(context.get_dof_indices()[i]);

      // Convenience
      libMesh::Real dt = _system.deltat;

      context.elem_solution_derivative = 1.0;

      // Local velocity at current time step
      // v_{n+1} = gamma/(beta*Delta t)*(x_{n+1}-x_n)
      //         + (1-(gamma/beta))*v_n
      //         + (1-gamma/(2*beta))*(Delta t)*a_n
      context.elem_solution_rate_derivative = (_gamma/(_beta*dt));

      context.get_elem_solution_rate()  = context.get_elem_solution();
      context.get_elem_solution_rate() -= old_elem_solution;
      context.get_elem_solution_rate() *= context.elem_solution_rate_derivative;
      context.get_elem_solution_rate().add( (1.0-_gamma/_beta), old_elem_solution_rate);
      context.get_elem_solution_rate().add( (1.0-_gamma/(2.0*_beta))*dt, old_elem_solution_accel);



      // Local acceleration at current time step
      // a_{n+1} = (1/(beta*(Delta t)^2))*(x_{n+1}-x_n)
      //         - 1/(beta*Delta t)*v_n
      //         - (1/(2*beta)-1)*a_n
      context.elem_solution_accel_derivative = 1.0/(_beta*dt*dt);

      context.get_elem_solution_accel()  = context.get_elem_solution();
//.........这里部分代码省略.........

UniquePtr<DiffContext> DifferentiableSystem::build_context ()
{
  DiffContext* context = new DiffContext(*this);
  context->set_deltat_pointer( &this->deltat );
  return UniquePtr<DiffContext>(context);
}

bool ElasticitySystem::mass_residual(bool request_jacobian,
                                     DiffContext & context)
{
  FEMContext & c = cast_ref<FEMContext &>(context);

  // We need to extract the corresponding velocity variable.
  // This allows us to use either a FirstOrderUnsteadySolver
  // or a SecondOrderUnsteadySolver. That is, we get back the velocity variable
  // index for FirstOrderUnsteadySolvers or, if it's a SecondOrderUnsteadySolver,
  // this is actually just giving us back the same variable index.

  // If we only wanted to use a SecondOrderUnsteadySolver, then this
  // step would be unnecessary and we would just
  // populate the _u_var, etc. blocks of the residual and Jacobian.
  unsigned int u_dot_var = this->get_second_order_dot_var(_u_var);
  unsigned int v_dot_var = this->get_second_order_dot_var(_v_var);
  unsigned int w_dot_var = this->get_second_order_dot_var(_w_var);

  FEBase * u_elem_fe;
  c.get_element_fe(u_dot_var, u_elem_fe);

  // The number of local degrees of freedom in each variable
  const unsigned int n_u_dofs = c.get_dof_indices(u_dot_var).size();

  // Element Jacobian * quadrature weights for interior integration
  const std::vector<Real> & JxW = u_elem_fe->get_JxW();

  const std::vector<std::vector<Real>> & phi = u_elem_fe->get_phi();

  // Residuals that we're populating
  libMesh::DenseSubVector<libMesh::Number> & Fu = c.get_elem_residual(u_dot_var);
  libMesh::DenseSubVector<libMesh::Number> & Fv = c.get_elem_residual(v_dot_var);
  libMesh::DenseSubVector<libMesh::Number> & Fw = c.get_elem_residual(w_dot_var);

  libMesh::DenseSubMatrix<libMesh::Number> & Kuu = c.get_elem_jacobian(u_dot_var, u_dot_var);
  libMesh::DenseSubMatrix<libMesh::Number> & Kvv = c.get_elem_jacobian(v_dot_var, v_dot_var);
  libMesh::DenseSubMatrix<libMesh::Number> & Kww = c.get_elem_jacobian(w_dot_var, w_dot_var);

  unsigned int n_qpoints = c.get_element_qrule().n_points();

  for (unsigned int qp=0; qp != n_qpoints; qp++)
    {
      // If we only cared about using FirstOrderUnsteadySolvers for time-stepping,
      // then we could actually just use interior rate, but using interior_accel
      // allows this assembly to function for both FirstOrderUnsteadySolvers
      // and SecondOrderUnsteadySolvers
      libMesh::Number u_ddot, v_ddot, w_ddot;
      c.interior_accel(u_dot_var, qp, u_ddot);
      c.interior_accel(v_dot_var, qp, v_ddot);
      c.interior_accel(w_dot_var, qp, w_ddot);

      for (unsigned int i=0; i != n_u_dofs; i++)
        {
          Fu(i) += _rho*u_ddot*phi[i][qp]*JxW[qp];
          Fv(i) += _rho*v_ddot*phi[i][qp]*JxW[qp];
          Fw(i) += _rho*w_ddot*phi[i][qp]*JxW[qp];

          if (request_jacobian)
            {
              for (unsigned int j=0; j != n_u_dofs; j++)
                {
                  libMesh::Real jac_term = _rho*phi[i][qp]*phi[j][qp]*JxW[qp];
                  jac_term *= context.get_elem_solution_accel_derivative();

                  Kuu(i,j) += jac_term;
                  Kvv(i,j) += jac_term;
                  Kww(i,j) += jac_term;
                }
            }
        }
    }

  // If the Jacobian was requested, we computed it. Otherwise, we didn't.
  return request_jacobian;
}

bool EulerSolver::element_residual (bool request_jacobian,
                                    DiffContext &context)
{
    unsigned int n_dofs = context.get_elem_solution().size();

    // Local nonlinear solution at old timestep
    DenseVector<Number> old_elem_solution(n_dofs);
    for (unsigned int i=0; i != n_dofs; ++i)
        old_elem_solution(i) =
            old_nonlinear_solution(context.get_dof_indices()[i]);

    // Local nonlinear solution at time t_theta
    DenseVector<Number> theta_solution(context.get_elem_solution());
    theta_solution *= theta;
    theta_solution.add(1. - theta, old_elem_solution);

    // Technically the elem_solution_derivative is either theta
    // or -1.0 in this implementation, but we scale the former part
    // ourselves
    context.elem_solution_derivative = 1.0;

// Technically the fixed_solution_derivative is always theta,
// but we're scaling a whole jacobian by theta after these first
// evaluations
    context.fixed_solution_derivative = 1.0;

    // If a fixed solution is requested, we'll use theta_solution
    if (_system.use_fixed_solution)
        context.get_elem_fixed_solution() = theta_solution;

    // Move theta_->elem_, elem_->theta_
    context.get_elem_solution().swap(theta_solution);

    // Move the mesh into place first if necessary
    context.elem_reinit(theta);

    // We're going to compute just the change in elem_residual
    // (and possibly elem_jacobian), then add back the old values
    DenseVector<Number> old_elem_residual(context.get_elem_residual());
    DenseMatrix<Number> old_elem_jacobian;
    if (request_jacobian)
    {
        old_elem_jacobian = context.get_elem_jacobian();
        context.get_elem_jacobian().zero();
    }
    context.get_elem_residual().zero();

    // Get the time derivative at t_theta
    bool jacobian_computed =
        _system.element_time_derivative(request_jacobian, context);

    // For a moving mesh problem we may need the pseudoconvection term too
    jacobian_computed =
        _system.eulerian_residual(jacobian_computed, context) && jacobian_computed;

    // Scale the time-dependent residual and jacobian correctly
    context.get_elem_residual() *= _system.deltat;
    if (jacobian_computed)
        context.get_elem_jacobian() *= (theta * _system.deltat);

    // The fixed_solution_derivative is always theta,
    // and now we're done scaling jacobians
    context.fixed_solution_derivative = theta;

    // We evaluate mass_residual with the change in solution
    // to get the mass matrix, reusing old_elem_solution to hold that
    // delta_solution.  We're solving dt*F(u) - du = 0, so our
    // delta_solution is old_solution - new_solution.
    // We're still keeping elem_solution in theta_solution for now
    old_elem_solution -= theta_solution;

    // Move old_->elem_, theta_->old_
    context.get_elem_solution().swap(old_elem_solution);

    // We do a trick here to avoid using a non-1
    // elem_solution_derivative:
    context.get_elem_jacobian() *= -1.0;
    context.fixed_solution_derivative *= -1.0;
    jacobian_computed = _system.mass_residual(jacobian_computed, context) &&
                        jacobian_computed;
    context.get_elem_jacobian() *= -1.0;
    context.fixed_solution_derivative *= -1.0;

    // Move elem_->elem_, old_->theta_
    context.get_elem_solution().swap(theta_solution);

    // Restore the elem position if necessary
    context.elem_reinit(1.);

    // Add the constraint term
    jacobian_computed = _system.element_constraint(jacobian_computed, context) &&
                        jacobian_computed;

    // Add back the old residual and jacobian
    context.get_elem_residual() += old_elem_residual;
    if (request_jacobian)
    {
        if (jacobian_computed)
            context.get_elem_jacobian() += old_elem_jacobian;
        else
//.........这里部分代码省略.........

bool Euler2Solver::element_residual (bool request_jacobian,
                                     DiffContext &context)
{
  unsigned int n_dofs = context.elem_solution.size();

  // Local nonlinear solution at old timestep
  DenseVector<Number> old_elem_solution(n_dofs);
  for (unsigned int i=0; i != n_dofs; ++i)
    old_elem_solution(i) =
      old_nonlinear_solution(context.dof_indices[i]);

  // We evaluate mass_residual with the change in solution
  // to get the mass matrix, reusing old_elem_solution to hold that
  // delta_solution.
  DenseVector<Number> delta_elem_solution(context.elem_solution);
  delta_elem_solution -= old_elem_solution;

  // Our first evaluations are at the true elem_solution
  context.elem_solution_derivative = 1.0;

  // If a fixed solution is requested, we'll use the elem_solution
  // at the new timestep
  if (_system.use_fixed_solution)
    context.elem_fixed_solution = context.elem_solution;

  context.fixed_solution_derivative = 1.0;

  // We're going to compute just the change in elem_residual
  // (and possibly elem_jacobian), then add back the old values
  DenseVector<Number> total_elem_residual(context.elem_residual);
  DenseMatrix<Number> old_elem_jacobian, total_elem_jacobian;
  if (request_jacobian)
    {
      old_elem_jacobian = context.elem_jacobian;
      total_elem_jacobian = context.elem_jacobian;
      context.elem_jacobian.zero();
    }
  context.elem_residual.zero();

  // First, evaluate time derivative at the new timestep.
  // The element should already be in the proper place
  // even for a moving mesh problem.
  bool jacobian_computed =
    _system.element_time_derivative(request_jacobian, context);

  // For a moving mesh problem we may need the pseudoconvection term too
  jacobian_computed =
    _system.eulerian_residual(jacobian_computed, context) && jacobian_computed;

  // Scale the new time-dependent residual and jacobian correctly
  context.elem_residual *= (theta * _system.deltat);
  total_elem_residual += context.elem_residual;
  context.elem_residual.zero();

  if (jacobian_computed)
    {
      context.elem_jacobian *= (theta * _system.deltat);
      total_elem_jacobian += context.elem_jacobian;
      context.elem_jacobian.zero();
    }

  // Next, evaluate the mass residual at the new timestep,
  // with the delta_solution.
  // Evaluating the mass residual at both old and new timesteps will be
  // redundant in most problems but may be necessary for time accuracy
  // or stability in moving mesh problems or problems with user-overridden
  // mass_residual functions

  // Move elem_->delta_, delta_->elem_
  context.elem_solution.swap(delta_elem_solution);

  jacobian_computed = _system.mass_residual(jacobian_computed, context) &&
    jacobian_computed;

  context.elem_residual *= -theta;
  total_elem_residual += context.elem_residual;
  context.elem_residual.zero();

  if (jacobian_computed)
    {
      // The minus sign trick here is to avoid using a non-1
      // elem_solution_derivative:
      context.elem_jacobian *= -theta;
      total_elem_jacobian += context.elem_jacobian;
      context.elem_jacobian.zero();
    }

  // Add the time-dependent term for the old solution

  // Make sure elem_solution is set up for elem_reinit to use
  // Move delta_->old_, old_->elem_
  context.elem_solution.swap(old_elem_solution);

  // Move the mesh into place first if necessary
  context.elem_reinit(0.);

  if (_system.use_fixed_solution)
    {
      context.elem_solution_derivative = 0.0;
      jacobian_computed =
//.........这里部分代码省略.........

bool EulerSolver::_general_residual (bool request_jacobian,
                                     DiffContext &context,
                                     ResFuncType mass,
                                     ResFuncType time_deriv,
                                     ResFuncType constraint,
                                     ReinitFuncType reinit_func)
{
  unsigned int n_dofs = context.get_elem_solution().size();

  // We might need to save the old jacobian in case one of our physics
  // terms later is unable to update it analytically.
  DenseMatrix<Number> old_elem_jacobian(n_dofs, n_dofs);
  if (request_jacobian)
    old_elem_jacobian.swap(context.get_elem_jacobian());

  // Local nonlinear solution at old timestep
  DenseVector<Number> old_elem_solution(n_dofs);
  for (unsigned int i=0; i != n_dofs; ++i)
    old_elem_solution(i) =
      old_nonlinear_solution(context.get_dof_indices()[i]);

  // Local time derivative of solution
  context.get_elem_solution_rate() = context.get_elem_solution();
  context.get_elem_solution_rate() -= old_elem_solution;
  context.elem_solution_rate_derivative = 1 / _system.deltat;
  context.get_elem_solution_rate() *=
    context.elem_solution_rate_derivative;

  // Local nonlinear solution at time t_theta
  DenseVector<Number> theta_solution(context.get_elem_solution());
  theta_solution *= theta;
  theta_solution.add(1. - theta, old_elem_solution);

  context.elem_solution_derivative = theta;
  context.fixed_solution_derivative = theta;

  // If a fixed solution is requested, we'll use theta_solution
  if (_system.use_fixed_solution)
    context.get_elem_fixed_solution() = theta_solution;

  // Move theta_->elem_, elem_->theta_
  context.get_elem_solution().swap(theta_solution);

  // Move the mesh into place first if necessary
  (context.*reinit_func)(theta);

  // Get the time derivative at t_theta
  bool jacobian_computed =
    (_system.*time_deriv)(request_jacobian, context);

  jacobian_computed = (_system.*mass)(jacobian_computed, context) &&
    jacobian_computed;

  // Restore the elem position if necessary
  (context.*reinit_func)(1);

  // Move elem_->elem_, theta_->theta_
  context.get_elem_solution().swap(theta_solution);
  context.elem_solution_derivative = 1;

  // Add the constraint term
  jacobian_computed = (_system.*constraint)(jacobian_computed, context) &&
    jacobian_computed;

  // Add back (or restore) the old jacobian
  if (request_jacobian)
    {
      if (jacobian_computed)
        context.get_elem_jacobian() += old_elem_jacobian;
      else
        context.get_elem_jacobian().swap(old_elem_jacobian);
    }

  return jacobian_computed;
}