C++ NumericVector::insert方法代码示例

//' LULU smoother.
//'
//' Performs LULU smoothing of the provided vector.
//'
//' @param x A real-valued vector.
//'
//' @return The LULU-smoothed version of \code{x}.
//'
//' @examples
//' x <- rnorm(10)
//' lulu(x)
//' @export
//[[Rcpp::export]]
NumericVector lulu(NumericVector x)
{
int n = x.size();
NumericVector x_out(n);

x.insert(0,0);
x.insert(n+1,0);

for (int i=1; i<=n; i++)
{
if ((x[i] < x[i-1]) & (x[i] < x[i+1]))
{
x_out[i-1] = std::min(x[i-1], x[i+1]);
}
else if ((x[i] > x[i-1]) & (x[i] > x[i+1]))
{
x_out[i-1] = std::max(x[i-1], x[i+1]);
}
else
{
x_out[i-1] = x[i];
}
}

return x_out;
}

void
MooseVariable::insert(NumericVector<Number> & residual)
{
  if (_has_nodal_value)
    residual.insert(&_nodal_u[0], _dof_indices);

  if (_has_nodal_value_neighbor)
    residual.insert(&_nodal_u_neighbor[0], _dof_indices_neighbor);
}

void
MooseVariableScalar::insert(NumericVector<Number> & soln)
{
  // We may have redundantly computed this value on many different
  // processors, but only the processor which actually owns it should
  // be saving it to the solution vector, to avoid O(N_scalar_vars)
  // unnecessary communication.

  const dof_id_type first_dof = _dof_map.first_dof();
  const dof_id_type end_dof = _dof_map.end_dof();
  if (_dof_indices.size() > 0 && first_dof <= _dof_indices[0] && _dof_indices[0] < end_dof)
    soln.insert(&_u[0], _dof_indices);
}

//.........这里部分代码省略.........
     std::fill(ctrl_upper.begin(), ctrl_upper.end(), 0.);
   
    for (unsigned i = 0; i < sol_actflag.size(); i++) {
        std::vector<double> node_coords_i(dim,0.);
        for (unsigned d = 0; d < dim; d++) node_coords_i[d] = coords_at_dofs[d][i];
        ctrl_lower[i] = InequalityConstraint(node_coords_i,false);
        ctrl_upper[i] = InequalityConstraint(node_coords_i,true);
        
    if      ( (sol_eldofs[pos_mu][i] + c_compl * (sol_eldofs[pos_ctrl][i] - ctrl_lower[i] )) < 0 )  sol_actflag[i] = 1;
    else if ( (sol_eldofs[pos_mu][i] + c_compl * (sol_eldofs[pos_ctrl][i] - ctrl_upper[i] )) > 0 )  sol_actflag[i] = 2;
    }
 
 //************** act flag **************************** 
    unsigned nDof_act_flag  = msh->GetElementDofNumber(iel, solFEType_act_flag);    // number of solution element dofs
    
    for (unsigned i = 0; i < nDof_act_flag; i++) {
      unsigned solDof_mu = msh->GetSolutionDof(i, iel, solFEType_act_flag); 
      (sol->_Sol[solIndex_act_flag])->set(solDof_mu,sol_actflag[i]);     
    }    

    //******************** ALL VARS ********************* 
    unsigned nDof_AllVars = 0;
    for (unsigned  k = 0; k < n_unknowns; k++) { nDof_AllVars += Sol_n_el_dofs[k]; }
 // TODO COMPUTE MAXIMUM maximum number of element dofs for one scalar variable
    int nDof_max    =  0;   
      for (unsigned  k = 0; k < n_unknowns; k++)     {
          if(Sol_n_el_dofs[k] > nDof_max)    nDof_max = Sol_n_el_dofs[k];
      }
  
    Res.resize(nDof_AllVars);                   std::fill(Res.begin(), Res.end(), 0.);
    Jac.resize(nDof_AllVars * nDof_AllVars);    std::fill(Jac.begin(), Jac.end(), 0.);
    
    L2G_dofmap_AllVars.resize(0);
      for (unsigned  k = 0; k < n_unknowns; k++)     L2G_dofmap_AllVars.insert(L2G_dofmap_AllVars.end(),L2G_dofmap[k].begin(),L2G_dofmap[k].end());
 //*************************************************** 


 //***** set control flag ****************************
  int control_el_flag = 0;
      control_el_flag = ControlDomainFlag_internal_restriction(elem_center);
  std::vector<int> control_node_flag(Sol_n_el_dofs[pos_ctrl],0);
  if (control_el_flag == 1) std::fill(control_node_flag.begin(), control_node_flag.end(), 1);
 //*************************************************** 
  



      // *** Gauss point loop ***
      for (unsigned ig = 0; ig < msh->_finiteElement[ielGeom][solFEType_max]->GetGaussPointNumber(); ig++) {
	
      // *** get gauss point weight, test function and test function partial derivatives ***
      for(unsigned int k = 0; k < n_unknowns; k++) {
         msh->_finiteElement[ielGeom][SolFEType[k]]->Jacobian(coords_at_dofs, ig, weight_qp, phi_fe_qp[k], phi_x_fe_qp[k], phi_xx_fe_qp[k]);
      }
      
   for (unsigned  k = 0; k < n_unknowns; k++) {
      for(unsigned int d = 0; d < dim; d++) {
        for(unsigned int i = 0; i < Sol_n_el_dofs[k]; i++) {
          phi_x_fe_qp_vec[k][d][i] =  phi_x_fe_qp[k][i * dim + d];
       }
     }
   }
      
      
   //HAVE TO RECALL IT TO HAVE BIQUADRATIC JACOBIAN
    for (unsigned int d = 0; d < dim; d++) {

void Biharmonic::JR::bounds(NumericVector<Number> & XL,
                            NumericVector<Number> & XU,
                            NonlinearImplicitSystem &)
{
  // sys is actually ignored, since it should be the same as *this.

  // Declare a performance log.  Give it a descriptive
  // string to identify what part of the code we are
  // logging, since there may be many PerfLogs in an
  // application.
  PerfLog perf_log ("Biharmonic bounds", false);

  // A reference to the DofMap object for this system.  The DofMap
  // object handles the index translation from node and element numbers
  // to degree of freedom numbers.  We will talk more about the DofMap
  // in future examples.
  const DofMap & dof_map = get_dof_map();

  // Get a constant reference to the Finite Element type
  // for the first (and only) variable in the system.
  FEType fe_type = dof_map.variable_type(0);

  // Build a Finite Element object of the specified type.  Since the
  // FEBase::build() member dynamically creates memory we will
  // store the object as a UniquePtr<FEBase>.  This can be thought
  // of as a pointer that will clean up after itself.
  UniquePtr<FEBase> fe (FEBase::build(_biharmonic._dim, fe_type));

  // Define data structures to contain the bound vectors contributions.
  DenseVector<Number> XLe, XUe;

  // These vector will hold the degree of freedom indices for
  // the element.  These define where in the global system
  // the element degrees of freedom get mapped.
  std::vector<dof_id_type> dof_indices;

  MeshBase::const_element_iterator       el     = _biharmonic._mesh.active_local_elements_begin();
  const MeshBase::const_element_iterator end_el = _biharmonic._mesh.active_local_elements_end();

  for ( ; el != end_el; ++el)
    {
      // Extract the shape function to be evaluated at the nodes
      const std::vector<std::vector<Real> > & phi = fe->get_phi();

      // Get the degree of freedom indices for the current element.
      // They are in 1-1 correspondence with shape functions phi
      // and define where in the global vector this element will.
      dof_map.dof_indices (*el, dof_indices);

      // Resize the local bounds vectors (zeroing them out in the process).
      XLe.resize(dof_indices.size());
      XUe.resize(dof_indices.size());

      // Extract the element node coordinates in the reference frame
      std::vector<Point> nodes;
      fe->get_refspace_nodes((*el)->type(), nodes);

      // Evaluate the shape functions at the nodes
      fe->reinit(*el, &nodes);

      // Construct the bounds based on the value of the i-th phi at the nodes.
      // Observe that this doesn't really work in general: we rely on the fact
      // that for Hermite elements each shape function is nonzero at most at a
      // single node.
      // More generally the bounds must be constructed by inspecting a "mass-like"
      // matrix (m_{ij}) of the shape functions (i) evaluated at their corresponding nodes (j).
      // The constraints imposed on the dofs (d_i) are then are -1 \leq \sum_i d_i m_{ij} \leq 1,
      // since \sum_i d_i m_{ij} is the value of the solution at the j-th node.
      // Auxiliary variables will need to be introduced to reduce this to a "box" constraint.
      // Additional complications will arise since m might be singular (as is the case for Hermite,
      // which, however, is easily handled by inspection).
      for (unsigned int i=0; i<phi.size(); ++i)
        {
          // FIXME: should be able to define INF and pass it to the solve
          Real infinity = 1.0e20;
          Real bound = infinity;
          for (unsigned int j = 0; j < nodes.size(); ++j)
            {
              if (phi[i][j])
                {
                  bound = 1.0/std::abs(phi[i][j]);
                  break;
                }
            }

          // The value of the solution at this node must be between 1.0 and -1.0.
          // Based on the value of phi(i)(i) the nodal coordinate must be between 1.0/phi(i)(i) and its negative.
          XLe(i) = -bound;
          XUe(i) = bound;
        }
      // The element bound vectors are now built for this element.
      // Insert them into the global vectors, potentially overwriting
      // the same dof contributions from other elements: no matter --
      // the bounds are always -1.0 and 1.0.
      XL.insert(XLe, dof_indices);
      XU.insert(XUe, dof_indices);
    }
}

void
MooseVariableScalar::insert(NumericVector<Number> & soln)
{
  if (_dof_indices.size() > 0)
    soln.insert(&_u[0], _dof_indices);
}