C++ AssemblyContext::interior_rate方法代码示例

  void ReactingLowMachNavierStokesStabilizationBase<Mixture,Evaluator>::compute_res_transient( AssemblyContext& context,
                                                                                               unsigned int qp,
                                                                                               libMesh::Real& RP_t,
                                                                                               libMesh::RealGradient& RM_t,
                                                                                               libMesh::Real& RE_t,
                                                                                               std::vector<libMesh::Real>& Rs_t )
  {
    libMesh::Real T = context.interior_value( this->_temp_vars.T(), qp );

    std::vector<libMesh::Real> ws(this->n_species());
    for(unsigned int s=0; s < this->_n_species; s++ )
      {
        ws[s] = context.interior_value(this->_species_vars.species(s), qp);
      }

    Evaluator gas_evaluator( this->_gas_mixture );
    const libMesh::Real R_mix = gas_evaluator.R_mix(ws);
    const libMesh::Real p0 = this->get_p0_transient(context,qp);
    const libMesh::Real rho = this->rho(T, p0, R_mix);
    const libMesh::Real cp = gas_evaluator.cp(T,p0,ws);
    const libMesh::Real M = gas_evaluator.M_mix( ws );

    // M_dot = -M^2 \sum_s w_dot[s]/Ms
    libMesh::Real M_dot = 0.0;
    std::vector<libMesh::Real> ws_dot(this->n_species());
    for(unsigned int s=0; s < this->n_species(); s++)
      {
        context.interior_rate(this->_species_vars.species(s), qp, ws_dot[s]);

        // Start accumulating M_dot
        M_dot += ws_dot[s]/this->_gas_mixture.M(s);
      }
    libMesh::Real M_dot_over_M = M_dot*(-M);

    libMesh::RealGradient u_dot;
    context.interior_rate(this->_flow_vars.u(), qp, u_dot(0));
    context.interior_rate(this->_flow_vars.v(), qp, u_dot(1));
    if(this->mesh_dim(context) == 3)
      context.interior_rate(this->_flow_vars.w(), qp, u_dot(2));

    libMesh::Real T_dot;
    context.interior_rate(this->_temp_vars.T(), qp, T_dot);

    RP_t = -T_dot/T + M_dot_over_M;
    RM_t = rho*u_dot;
    RE_t = rho*cp*T_dot;
    for(unsigned int s=0; s < this->n_species(); s++)
      {
        Rs_t[s] = rho*ws_dot[s];
      }

    return;
  }

  void HeatConduction<K>::mass_residual( bool compute_jacobian,
				      AssemblyContext& context,
				      CachedValues& /*cache*/ )
  {
    // 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<libMesh::Real> &JxW = 
      context.get_element_fe(_temp_vars.T_var())->get_JxW();

    // The shape functions at interior quadrature points.
    const std::vector<std::vector<libMesh::Real> >& phi = 
      context.get_element_fe(_temp_vars.T_var())->get_phi();

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

    // The subvectors and submatrices we need to fill:
    libMesh::DenseSubVector<libMesh::Real> &F =
      context.get_elem_residual(_temp_vars.T_var());

    libMesh::DenseSubMatrix<libMesh::Real> &M =
      context.get_elem_jacobian(_temp_vars.T_var(), _temp_vars.T_var());

    unsigned int n_qpoints = context.get_element_qrule().n_points();
    
    for (unsigned int qp = 0; qp != n_qpoints; ++qp)
      {
	// For the mass residual, we need to be a little careful.
	// The time integrator is handling the time-discretization
	// for us so we need to supply M(u_fixed)*u' for the residual.
	// u_fixed will be given by the fixed_interior_value function
	// while u' will be given by the interior_rate function.
        libMesh::Real T_dot;
        context.interior_rate(_temp_vars.T_var(), qp, T_dot);

	for (unsigned int i = 0; i != n_T_dofs; ++i)
	  {
	    F(i) -= JxW[qp]*(_rho*_Cp*T_dot*phi[i][qp] );

	    if( compute_jacobian )
              {
                for (unsigned int j=0; j != n_T_dofs; j++)
                  {
		    // We're assuming rho, cp are constant w.r.t. T here.
                    M(i,j) -=
                      context.get_elem_solution_rate_derivative()
                        * JxW[qp]*_rho*_Cp*phi[j][qp]*phi[i][qp] ;
                  }
              }// End of check on Jacobian

	  } // End of element dof loop

      } // End of the quadrature point loop

    return;
  }

  libMesh::Real HeatTransferStabilizationHelper::compute_res_energy_transient( AssemblyContext& context,
                                                                               unsigned int qp,
                                                                               const libMesh::Real rho,
                                                                               const libMesh::Real Cp ) const
  {
    libMesh::Real T_dot;
    context.interior_rate(this->_temp_vars.T(), qp, T_dot);

    return rho*Cp*T_dot;
  }

  void HeatTransferStabilizationHelper::compute_res_energy_transient_and_derivs
  ( AssemblyContext& context,
    unsigned int qp,
    const libMesh::Real rho,
    const libMesh::Real Cp,
    libMesh::Real &res,
    libMesh::Real &d_res_dTdot
    ) const
  {
    libMesh::Real T_dot;
    context.interior_rate(this->_temp_vars.T(), qp, T_dot);

    res = rho*Cp*T_dot;
    d_res_dTdot = rho*Cp;
  }