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

  void BoundaryConditions::apply_neumann_normal( AssemblyContext& context,
                                                 const VariableIndex var,
                                                 const libMesh::Real sign,
                                                 const FEShape& value ) const
  {
    libMesh::FEGenericBase<FEShape>* side_fe = NULL; 
    context.get_side_fe( var, side_fe );

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

    // Element Jacobian * quadrature weight for side integration.
    const std::vector<libMesh::Real> &JxW_side = side_fe->get_JxW();

    // The var shape functions at side quadrature points.
    const std::vector<std::vector<FEShape> >& var_phi_side = side_fe->get_phi();

    libMesh::DenseSubVector<libMesh::Number> &F_var = context.get_elem_residual(var); // residual

    unsigned int n_qpoints = context.get_side_qrule().n_points();
    for (unsigned int qp=0; qp != n_qpoints; qp++)
      {
	for (unsigned int i=0; i != n_var_dofs; i++)
	  {
	    F_var(i) += sign*value*var_phi_side[i][qp]*JxW_side[qp];
	  }
      }

    return;
  }

  void AverageNusseltNumber::side_qoi_derivative( AssemblyContext& context,
                                                  const unsigned int qoi_index )
  {

    for( std::set<libMesh::boundary_id_type>::const_iterator id = _bc_ids.begin();
	 id != _bc_ids.end(); id++ )
      {
	if( context.has_side_boundary_id( (*id) ) )
	  {
	    libMesh::FEBase* T_side_fe;
	    context.get_side_fe<libMesh::Real>(this->_T_var, T_side_fe);

	    const std::vector<libMesh::Real> &JxW = T_side_fe->get_JxW();
	    
	    const std::vector<libMesh::Point>& normals = T_side_fe->get_normals();

	    unsigned int n_qpoints = context.get_side_qrule().n_points();
	    
            const unsigned int n_T_dofs = context.get_dof_indices(_T_var).size();

            const std::vector<std::vector<libMesh::Gradient> >& T_gradphi = T_side_fe->get_dphi();

	    libMesh::DenseSubVector<libMesh::Number>& dQ_dT =
              context.get_qoi_derivatives(qoi_index, _T_var);

	    // Loop over quadrature points  
	    for (unsigned int qp = 0; qp != n_qpoints; qp++)
	      {
		// Get the solution value at the quadrature point
		libMesh::Gradient grad_T = 0.0; 
		context.side_gradient(this->_T_var, qp, grad_T);
		
		// Update the elemental increment dR for each qp
		//qoi += (this->_scaling)*(this->_k)*(grad_T*normals[qp])*JxW[qp];

                for( unsigned int i = 0; i != n_T_dofs; i++ )
                  {
                    dQ_dT(i) += _scaling*_k*T_gradphi[i][qp]*normals[qp]*JxW[qp];
                  }

	      } // quadrature loop

	  } // end check on boundary id

      }

    return;
  }

  void LowMachNavierStokes<Mu,SH,TC>::assemble_thermo_press_side_time_deriv( bool /*compute_jacobian*/,
									     AssemblyContext& context )
  {
    // The number of local degrees of freedom in each variable.
    const unsigned int n_p0_dofs = context.get_dof_indices(this->_p0_var).size();

    // Element Jacobian * quadrature weight for side integration.
    //const std::vector<libMesh::Real> &JxW_side = context.get_side_fe(this->_T_var)->get_JxW();

    //const std::vector<Point> &normals = context.get_side_fe(this->_T_var)->get_normals();

    //libMesh::DenseSubVector<libMesh::Number> &F_p0 = context.get_elem_residual(this->_p0_var); // residual

    // Physical location of the quadrature points
    //const std::vector<libMesh::Point>& qpoint = context.get_side_fe(this->_T_var)->get_xyz();

    unsigned int n_qpoints = context.get_side_qrule().n_points();
    for (unsigned int qp=0; qp != n_qpoints; qp++)
      {
	/*
	libMesh::Number T = context.side_value( this->_T_var, qp );
	libMesh::Gradient U = ( context.side_value( this->_u_var, qp ),
				context.side_value( this->_v_var, qp ) );
	libMesh::Gradient grad_T = context.side_gradient( this->_T_var, qp );

	
	libMesh::Number p0 = context.side_value( this->_p0_var, qp );

	libMesh::Number k = this->_k(T);
	libMesh::Number cp = this->_cp(T);

	libMesh::Number cv = cp + this->_R;
	libMesh::Number gamma = cp/cv;
	libMesh::Number gamma_ratio = gamma/(gamma-1.0);
	*/

	//std::cout << "U = " << U << std::endl;

	//std::cout << "x = " << qpoint[qp] << ", grad_T = " << grad_T << std::endl;

	for (unsigned int i=0; i != n_p0_dofs; i++)
	  {
	    //F_p0(i) += (k*grad_T*normals[qp] - p0*gamma_ratio*U*normals[qp]  )*JxW_side[qp];
	  }
      }
  
    return;
  }

  void AverageNusseltNumber::side_qoi( AssemblyContext& context,
                                       const unsigned int qoi_index )
  {
    bool on_correct_side = false;

    for (std::set<libMesh::boundary_id_type>::const_iterator id =
         _bc_ids.begin(); id != _bc_ids.end(); id++ )
      if( context.has_side_boundary_id( (*id) ) )
        {
          on_correct_side = true;
          break;
        }

    if (!on_correct_side)
      return;

    libMesh::FEBase* side_fe;
    context.get_side_fe<libMesh::Real>(this->_T_var, side_fe);

    const std::vector<libMesh::Real> &JxW = side_fe->get_JxW();

    const std::vector<libMesh::Point>& normals = side_fe->get_normals();

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

    libMesh::Number& qoi = context.get_qois()[qoi_index];

    // Loop over quadrature points

    for (unsigned int qp = 0; qp != n_qpoints; qp++)
      {
	// Get the solution value at the quadrature point
	libMesh::Gradient grad_T = 0.0;
	context.side_gradient(this->_T_var, qp, grad_T);

	// Update the elemental increment dR for each qp
	qoi += (this->_scaling)*(this->_k)*(grad_T*normals[qp])*JxW[qp];

      } // quadrature loop
  }

  void ParsedBoundaryQoI::side_qoi( AssemblyContext& context,
                                    const unsigned int qoi_index )
  {
    libMesh::FEBase* side_fe;
    context.get_side_fe<libMesh::Real>(0, side_fe);
    const std::vector<libMesh::Real> &JxW = side_fe->get_JxW();

    const std::vector<libMesh::Point>& x_qp = side_fe->get_xyz();

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

    /*! \todo Need to generalize this to the multiple QoI case */
    libMesh::Number& qoi = context.get_qois()[qoi_index];

    for( unsigned int qp = 0; qp != n_qpoints; qp++ )
      {
        const libMesh::Number func_val =
          (*qoi_functional)(context, x_qp[qp], context.get_time());

        qoi += func_val * JxW[qp];
      }
  }

  bool GasRecombinationCatalyticWall<Chemistry>::eval_flux( bool compute_jacobian,
                                                            AssemblyContext& context,
                                                            libMesh::Real sign,
                                                            bool is_axisymmetric )
  {
    libMesh::FEGenericBase<libMesh::Real>* side_fe = NULL;
    context.get_side_fe( _reactant_var_idx, side_fe );

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

    libmesh_assert_equal_to( n_var_dofs, context.get_dof_indices(_product_var_idx).size() );

    // Element Jacobian * quadrature weight for side integration.
    const std::vector<libMesh::Real> &JxW_side = side_fe->get_JxW();

    // The var shape functions at side quadrature points.
    const std::vector<std::vector<libMesh::Real> >& var_phi_side = side_fe->get_phi();

    // Physical location of the quadrature points
    const std::vector<libMesh::Point>& var_qpoint = side_fe->get_xyz();

    // reactant residual
    libMesh::DenseSubVector<libMesh::Number> &F_r_var = context.get_elem_residual(_reactant_var_idx);

    // product residual
    libMesh::DenseSubVector<libMesh::Number> &F_p_var = context.get_elem_residual(_product_var_idx);

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

    for (unsigned int qp=0; qp != n_qpoints; qp++)
      {
        libMesh::Real jac = JxW_side[qp];

        if(is_axisymmetric)
          {
            const libMesh::Number r = var_qpoint[qp](0);
            jac *= r;
          }

        std::vector<libMesh::Real> mass_fractions(this->_chem_ptr->n_species());
        for( unsigned int s = 0; s < this->_chem_ptr->n_species(); s++ )
          mass_fractions[s] = context.side_value(this->_species_vars[s], qp);

        libMesh::Real Y_r = mass_fractions[this->_reactant_species_idx];
        libMesh::Real T =  context.side_value(this->_T_var, qp);
        libMesh::Real R_mix = this->_chem_ptr->R_mix(mass_fractions);
        libMesh::Real rho = this->rho( T, this->_p0, R_mix );

         const libMesh::Real r_value = this->compute_reactant_mass_flux(rho, Y_r, T);

         const libMesh::Real p_value = -r_value;

         for (unsigned int i=0; i != n_var_dofs; i++)
          {
            F_r_var(i) += sign*r_value*var_phi_side[i][qp]*jac;

            F_p_var(i) += sign*p_value*var_phi_side[i][qp]*jac;

            if( compute_jacobian )
              libmesh_not_implemented();
          }
      }

    // We're not computing the Jacobian yet
    return false;
  }

  void ParsedBoundaryQoI::side_qoi_derivative( AssemblyContext& context,
                                               const unsigned int qoi_index )
  {
    bool on_correct_side = false;

    for (std::set<libMesh::boundary_id_type>::const_iterator id =
         _bc_ids.begin(); id != _bc_ids.end(); id++ )
      if( context.has_side_boundary_id( (*id) ) )
        {
          on_correct_side = true;
          break;
        }

    if (!on_correct_side)
      return;

    libMesh::FEBase* side_fe;
    context.get_side_fe<libMesh::Real>(0, side_fe);
    const std::vector<libMesh::Real> &JxW = side_fe->get_JxW();

    const std::vector<libMesh::Point>& x_qp = side_fe->get_xyz();

    // Local DOF count and quadrature point count
    const unsigned int n_u_dofs = context.get_dof_indices().size();

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

    // Local solution vector - non-const version for finite
    // differenting purposes
    libMesh::DenseVector<libMesh::Number>& elem_solution =
      const_cast<libMesh::DenseVector<libMesh::Number>&>
        (context.get_elem_solution());

    /*! \todo Need to generalize this to the multiple QoI case */
    libMesh::DenseVector<libMesh::Number> &Qu =
      context.get_qoi_derivatives()[qoi_index];

    for( unsigned int qp = 0; qp != n_qpoints; qp++ )
      {
        // Central finite differencing to approximate derivatives.
        // FIXME - we should hook the FParserAD stuff into
        // ParsedFEMFunction

        for( unsigned int i = 0; i != n_u_dofs; ++i )
          {
            libMesh::Number &current_solution = elem_solution(i);
            const libMesh::Number original_solution = current_solution;

            current_solution = original_solution + libMesh::TOLERANCE;

            const libMesh::Number plus_val =
              (*qoi_functional)(context, x_qp[qp], context.get_time());

            current_solution = original_solution - libMesh::TOLERANCE;

            const libMesh::Number minus_val =
              (*qoi_functional)(context, x_qp[qp], context.get_time());

            Qu(i) += (plus_val - minus_val) *
                     (0.5 / libMesh::TOLERANCE) * JxW[qp];

            // Don't forget to restore the correct solution...
            current_solution = original_solution;
          }
      }
  }

  void BoundaryConditions::apply_neumann_axisymmetric( AssemblyContext& context,
                                                       const CachedValues& cache,
                                                       const bool request_jacobian,
                                                       const VariableIndex var,
                                                       const libMesh::Real sign,
                                                       SharedPtr<NeumannFuncObj> neumann_func ) const
  {
    libMesh::FEGenericBase<libMesh::Real>* side_fe = NULL; 
    context.get_side_fe( var, side_fe );

    // The number of local degrees of freedom
    const unsigned int n_var_dofs = context.get_dof_indices(var).size();
  
    // Element Jacobian * quadrature weight for side integration.
    const std::vector<libMesh::Real> &JxW_side = side_fe->get_JxW();

    // The var shape functions at side quadrature points.
    const std::vector<std::vector<libMesh::Real> >& var_phi_side =
      side_fe->get_phi();

    // Physical location of the quadrature points
    const std::vector<libMesh::Point>& var_qpoint =
      side_fe->get_xyz();

    const std::vector<libMesh::Point> &normals = side_fe->get_normals();

    libMesh::DenseSubVector<libMesh::Number> &F_var = context.get_elem_residual(var); // residual
    libMesh::DenseSubMatrix<libMesh::Number> &K_var = context.get_elem_jacobian(var,var); // jacobian

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

    for (unsigned int qp=0; qp != n_qpoints; qp++)
      {
	const libMesh::Point bc_value = neumann_func->value( context, cache, qp );
        libMesh::Point jac_value;
        if (request_jacobian)
          {
            jac_value = neumann_func->derivative( context, cache, qp );
          }

	const libMesh::Number r = var_qpoint[qp](0);

	for (unsigned int i=0; i != n_var_dofs; i++)
	  {
	    F_var(i) += sign*r*JxW_side[qp]*bc_value*normals[qp]*var_phi_side[i][qp];

	    if (request_jacobian)
	      {
		for (unsigned int j=0; j != n_var_dofs; j++)
		  {
		    K_var(i,j) += sign*r*JxW_side[qp]*jac_value*normals[qp]*
		      var_phi_side[i][qp]*var_phi_side[j][qp];
		  }
	      }
	  }
      } // End quadrature loop

    // Now must take care of the case that the boundary condition depends on variables
    // other than var.
    std::vector<VariableIndex> other_jac_vars = neumann_func->get_other_jac_vars();

    if( (other_jac_vars.size() > 0) && request_jacobian )
      {
	for( std::vector<VariableIndex>::const_iterator var2 = other_jac_vars.begin();
	     var2 != other_jac_vars.end();
	     var2++ )
	  {
            libMesh::FEGenericBase<libMesh::Real>* side_fe2 = NULL; 
            context.get_side_fe( *var2, side_fe2 );

            libMesh::DenseSubMatrix<libMesh::Number> &K_var2 = context.get_elem_jacobian(var,*var2); // jacobian

	    const unsigned int n_var2_dofs = context.get_dof_indices(*var2).size();
	    const std::vector<std::vector<libMesh::Real> >& var2_phi_side =
              side_fe2->get_phi();

	    for (unsigned int qp=0; qp != n_qpoints; qp++)
	      {
		const libMesh::Number r = var_qpoint[qp](0);

		const libMesh::Point jac_value = neumann_func->derivative( context, cache, qp, *var2 );

		for (unsigned int i=0; i != n_var_dofs; i++)
		  {
		    for (unsigned int j=0; j != n_var2_dofs; j++)
		      {
			K_var2(i,j) += sign*r*JxW_side[qp]*jac_value*normals[qp]*
			  var_phi_side[i][qp]*var2_phi_side[j][qp];
		      }
		  }
	      }
	  } // End loop over auxillary Jacobian variables
      }
    return;
  }