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

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

//.........这里部分代码省略.........

    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]);
						J_zc_fc(i,j) += JxW[qp]*(phi[j][qp]*phi[i][qp]);
	
	       		J_fc_zc(i,j) += JxW[qp]*(phi[j][qp]*phi[i][qp]);
       			J_fc_fc(i,j) += JxW[qp]*(_beta*dphi[j][qp]*dphi[i][qp]);
					} // end of the inner dof (j) loop
			  } // end - if (compute_jacobian && context.get_elem_solution_derivative())

			} // end of the outer dof (i) loop
    } // end of the quadrature point (qp) loop
    
	  for(unsigned int dnum=0; dnum<datavals.size(); dnum++){
	  	libMesh::Point data_point = datapts[dnum];
	  	if(context.get_elem().contains_point(data_point)){
	  		libMesh::Number cpred = context.point_value(_c_var, data_point);
	  		libMesh::Number cstar = datavals[dnum];
	  		
	  		unsigned int dim = context.get_system().get_mesh().mesh_dimension();
		    libMesh::FEType fe_type = context.get_element_fe(_c_var)->get_fe_type();
		    
		    //go between physical and reference element
		    libMesh::Point c_master = libMesh::FEInterface::inverse_map(dim, fe_type, &context.get_elem(), data_point); 	
		    
        std::vector<libMesh::Real> point_phi(n_c_dofs);
      	for (unsigned int i=0; i != n_c_dofs; i++){
      		//get value of basis function at mapped point in reference (master) element
          point_phi[i] = libMesh::FEInterface::shape(dim, fe_type, &context.get_elem(), i, c_master); 
        }
        
        for (unsigned int i=0; i != n_c_dofs; i++){
  	  		Rc(i) += (cpred - cstar)*point_phi[i];
	  
					if (compute_jacobian){
						for (unsigned int j=0; j != n_c_dofs; j++)
							J_c_c(i,j) += point_phi[j]*point_phi[i] ;
				  }
	  
  			}
	  	}
	  }

    return;
	
	}