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

  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 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_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;
          }
      }
  }