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

  void BoundaryConditions::pin_value( AssemblyContext& context,
                                      const CachedValues& /*cache*/,
                                      const bool request_jacobian,
                                      const VariableIndex var, 
                                      const double pin_value,
                                      const libMesh::Point& pin_location, 
                                      const double penalty )
  {
    if (context.get_elem().contains_point(pin_location))
      {
        libMesh::FEGenericBase<libMesh::Real>* elem_fe = NULL; 
        context.get_element_fe( var, elem_fe );

	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

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

	libMesh::Number var_value = context.point_value(var, pin_location);

	libMesh::FEType fe_type = elem_fe->get_fe_type();
      
	libMesh::Point point_loc_in_masterelem = 
	  libMesh::FEInterface::inverse_map(context.get_dim(), fe_type, &context.get_elem(), pin_location);

	std::vector<libMesh::Real> phi(n_var_dofs);

	for (unsigned int i=0; i != n_var_dofs; i++)
          {
            phi[i] = libMesh::FEInterface::shape( context.get_dim(), fe_type, &context.get_elem(), i, 
                                                  point_loc_in_masterelem );
          }
      
	for (unsigned int i=0; i != n_var_dofs; i++)
	  {
	    F_var(i) += penalty*(var_value - pin_value)*phi[i];
	  
	    /** \todo What the hell is the context.get_elem_solution_derivative() all about? */
	    if (request_jacobian && context.get_elem_solution_derivative())
	      {
		libmesh_assert (context.get_elem_solution_derivative() == 1.0);
	      
		for (unsigned int j=0; j != n_var_dofs; j++)
		  K_var(i,j) += penalty*phi[i]*phi[j];

	      } // End if request_jacobian
	  } // End i loop
      } // End if pin_location

    return;
  }

 void MultiphysicsSystem::compute_postprocessed_quantity( unsigned int quantity_index,
                                                          const AssemblyContext& context,
                                                          const libMesh::Point& point,
                                                          libMesh::Real& value )
 {
   for( PhysicsListIter physics_iter = _physics_list.begin();
        physics_iter != _physics_list.end();
        physics_iter++ )
     {
       // Only compute if physics is active on current subdomain or globally
       if( (physics_iter->second)->enabled_on_elem( &context.get_elem() ) )
         {
           (physics_iter->second)->compute_postprocessed_quantity( quantity_index, context, point, value );
         }
     }
   return;
 }

  void SpalartAllmarasSPGSMStabilization<Mu>::element_time_derivative
  ( bool compute_jacobian,
    AssemblyContext & context )
  {
    // Get a pointer to the current element, we need this for computing the distance to wall for the
    // quadrature points
    libMesh::Elem &elem_pointer = context.get_elem();

    // The number of local degrees of freedom in each variable.
    const unsigned int n_nu_dofs = context.get_dof_indices(this->_turbulence_vars.nu()).size();

    // Element Jacobian * quadrature weights for interior integration.
    const std::vector<libMesh::Real> &JxW =
      context.get_element_fe(this->_turbulence_vars.nu())->get_JxW();

    // The viscosity shape function gradients (in global coords.)
    // at interior quadrature points.
    const std::vector<std::vector<libMesh::RealGradient> >& nu_gradphi =
      context.get_element_fe(this->_turbulence_vars.nu())->get_dphi();

    // Quadrature point locations
    //const std::vector<libMesh::Point>& nu_qpoint =
    //context.get_element_fe(this->_turbulence_vars.nu())->get_xyz();

    //libMesh::DenseSubMatrix<libMesh::Number> &Knunu = context.get_elem_jacobian(this->_turbulence_vars.nu(), this->_turbulence_vars.nu()); // R_{nu},{nu}

    libMesh::DenseSubVector<libMesh::Number> &Fnu = context.get_elem_residual(this->_turbulence_vars.nu()); // R_{nu}

    libMesh::FEBase* fe = context.get_element_fe(this->_turbulence_vars.nu());

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

    // Auto pointer to distance fcn evaluated at quad points
    std::unique_ptr< libMesh::DenseVector<libMesh::Real> > distance_qp;

    // Fill the vector of distances to quadrature points
    distance_qp = this->distance_function->interpolate(&elem_pointer, context.get_element_qrule().get_points());

    for (unsigned int qp=0; qp != n_qpoints; qp++)
      {
        libMesh::Gradient grad_nu;
        grad_nu = context.interior_gradient(this->_turbulence_vars.nu(), qp);

        libMesh::Real jac = JxW[qp];

        // The physical viscosity
        libMesh::Real _mu_qp = this->_mu(context, qp);

        // To be fixed
        // For the channel flow we will just set the distance function analytically
        //(*distance_qp)(qp) = std::min(fabs(y),fabs(1 - y));

        // The flow velocity
        libMesh::Number u,v;
        u = context.interior_value(this->_flow_vars.u(), qp);
        v = context.interior_value(this->_flow_vars.v(), qp);

        libMesh::NumberVectorValue U(u,v);
        if (this->_flow_vars.dim() == 3)
          U(2) = context.interior_value(this->_flow_vars.w(), qp);

        // Stabilization terms

        libMesh::RealGradient g = this->_stab_helper.compute_g( fe, context, qp );
        libMesh::RealTensor G = this->_stab_helper.compute_G( fe, context, qp );

        libMesh::Real tau_spalart = this->_stab_helper.compute_tau_spalart( context, qp, g, G, this->_rho, U, _mu_qp, this->_is_steady );

        libMesh::Number RM_spalart = this->_stab_helper.compute_res_spalart_steady( context, qp, this->_rho, _mu_qp, (*distance_qp)(qp), this->_infinite_distance );

        for (unsigned int i=0; i != n_nu_dofs; i++)
          {
            Fnu(i) += jac*( -tau_spalart*RM_spalart*this->_rho*(U*nu_gradphi[i][qp]) );
          }

        if( compute_jacobian )
          {
            libmesh_not_implemented();
          }

      }
  }

	void PracticeCDRinv::element_time_derivative( bool compute_jacobian,
						AssemblyContext& context,
						CachedValues& /*cache*/ ){
	
		// The number of local degrees of freedom in each variable.
    const unsigned int n_c_dofs = context.get_dof_indices(_c_var).size();

    // 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(_c_var)->get_JxW();

    // The temperature shape function gradients (in global coords.)
    // at interior quadrature points.
    const std::vector<std::vector<libMesh::RealGradient> >& dphi =
      context.get_element_fe(_c_var)->get_dphi();
    const std::vector<std::vector<libMesh::Real> >& phi = context.get_element_fe(_c_var)->get_phi();

    const std::vector<libMesh::Point>& q_points = 
      context.get_element_fe(_c_var)->get_xyz();
    
  	libMesh::DenseSubMatrix<libMesh::Number> &J_c_zc = context.get_elem_jacobian(_c_var, _zc_var);
		libMesh::DenseSubMatrix<libMesh::Number> &J_c_c = context.get_elem_jacobian(_c_var, _c_var);
	
		libMesh::DenseSubMatrix<libMesh::Number> &J_zc_c = context.get_elem_jacobian(_zc_var, _c_var);
		libMesh::DenseSubMatrix<libMesh::Number> &J_zc_fc = context.get_elem_jacobian(_zc_var, _fc_var);
	
		libMesh::DenseSubMatrix<libMesh::Number> &J_fc_zc = context.get_elem_jacobian(_fc_var, _zc_var);
		libMesh::DenseSubMatrix<libMesh::Number> &J_fc_fc = context.get_elem_jacobian(_fc_var, _fc_var);
		
		libMesh::DenseSubVector<libMesh::Number> &Rc = context.get_elem_residual( _c_var );;
		libMesh::DenseSubVector<libMesh::Number> &Rzc = context.get_elem_residual( _zc_var );
		libMesh::DenseSubVector<libMesh::Number> &Rfc = context.get_elem_residual( _fc_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 = context.get_element_qrule().n_points();

    for (unsigned int qp=0; qp != n_qpoints; qp++){

			libMesh::Number 
	      c = context.interior_value(_c_var, qp),
	      zc = context.interior_value(_zc_var, qp),
	      fc = context.interior_value(_fc_var, qp);
	    libMesh::Gradient 
	      grad_c = context.interior_gradient(_c_var, qp),
	      grad_zc = context.interior_gradient(_zc_var, qp),
	      grad_fc = context.interior_gradient(_fc_var, qp);
			
	  	//location of quadrature point
	  	const libMesh::Real ptx = q_points[qp](0);
	  	const libMesh::Real pty = q_points[qp](1);
			
   		int xind, yind;
   		libMesh::Real xdist = 1.e10; libMesh::Real ydist = 1.e10;
   		for(int ii=0; ii<x_pts.size(); ii++){
   			libMesh::Real tmp = std::abs(ptx - x_pts[ii]);
   			if(xdist > tmp){
   				xdist = tmp;
   				xind = ii;
   			}
   			else
   				break;
   		} 
   		for(int jj=0; jj<y_pts[xind].size(); jj++){
   			libMesh::Real tmp = std::abs(pty - y_pts[xind][jj]);
   			if(ydist > tmp){
   				ydist = tmp;
   				yind = jj;
   			}
   			else
   				break;
   		}
   		libMesh::Real u = vel_field[xind][yind](0);
   		libMesh::Real v = vel_field[xind][yind](1);

	    libMesh::NumberVectorValue U     (u,     v);

	
			// First, an i-loop over the  degrees of freedom.
			for (unsigned int i=0; i != n_c_dofs; i++){
				
				Rc(i) += JxW[qp]*(-_k*grad_zc*dphi[i][qp] + U*grad_zc*phi[i][qp] + 2*_R*zc*c*phi[i][qp]);
	      Rzc(i) += JxW[qp]*(-_k*grad_c*dphi[i][qp] - U*grad_c*phi[i][qp] + _R*c*c*phi[i][qp] + fc*phi[i][qp]);
     		Rfc(i) += JxW[qp]*(_beta*grad_fc*dphi[i][qp] + zc*phi[i][qp]);

				if (compute_jacobian){
					for (unsigned int j=0; j != n_c_dofs; j++){
						J_c_zc(i,j) += JxW[qp]*(-_k*dphi[j][qp]*dphi[i][qp] + U*dphi[j][qp]*phi[i][qp] 
															+ 2*_R*phi[j][qp]*c*phi[i][qp]);
						J_c_c(i,j) += JxW[qp]*(2*_R*zc*phi[j][qp]*phi[i][qp]);

						J_zc_c(i,j) += JxW[qp]*(-_k*dphi[j][qp]*dphi[i][qp] - U*dphi[j][qp]*phi[i][qp] 
																+ 2*_R*c*phi[j][qp]*phi[i][qp]);
//.........这里部分代码省略.........