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

void
MoosePreconditioner::copyVarValues(MeshBase & mesh,
                                   const unsigned int from_system,
                                   const unsigned int from_var,
                                   const NumericVector<Number> & from_vector,
                                   const unsigned int to_system,
                                   const unsigned int to_var,
                                   NumericVector<Number> & to_vector)
{
  {
    MeshBase::node_iterator it = mesh.local_nodes_begin();
    MeshBase::node_iterator it_end = mesh.local_nodes_end();

    for (; it != it_end; ++it)
    {
      Node * node = *it;

      unsigned int n_comp = node->n_comp(from_system, from_var);

      mooseAssert(node->n_comp(from_system, from_var) == node->n_comp(to_system, to_var),
                  "Number of components does not match in each system");

      for (unsigned int i = 0; i < n_comp; i++)
      {
        dof_id_type from_dof = node->dof_number(from_system, from_var, i);
        dof_id_type to_dof = node->dof_number(to_system, to_var, i);

        to_vector.set(to_dof, from_vector(from_dof));
      }
    }
  }
  {
    MeshBase::element_iterator it = mesh.local_elements_begin();
    MeshBase::element_iterator it_end = mesh.local_elements_end();

    for (; it != it_end; ++it)
    {
      Elem * elem = *it;

      unsigned int n_comp = elem->n_comp(from_system, from_var);

      mooseAssert(elem->n_comp(from_system, from_var) == elem->n_comp(to_system, to_var),
                  "Number of components does not match in each system");

      for (unsigned int i = 0; i < n_comp; i++)
      {
        dof_id_type from_dof = elem->dof_number(from_system, from_var, i);
        dof_id_type to_dof = elem->dof_number(to_system, to_var, i);

        to_vector.set(to_dof, from_vector(from_dof));
      }
    }
  }
}

void
MooseVariable::insert(NumericVector<Number> & residual)
{
  if (_has_nodal_value)
  {
    for (unsigned int i=0; i<_nodal_u.size(); i++)
      residual.set(_dof_indices[i], _nodal_u[i]);
  }

  if (_has_nodal_value_neighbor)
  {
    for (unsigned int i=0; i<_nodal_u_neighbor.size(); i++)
      residual.set(_dof_indices_neighbor[i], _nodal_u_neighbor[0]);
  }
}

void
NodeFaceConstraint::computeSlaveValue(NumericVector<Number> & current_solution)
{
  dof_id_type & dof_idx = _var.nodalDofIndex();
  _qp = 0;
  current_solution.set(dof_idx, computeQpSlaveValue());
}

void
PresetNodalBC::computeValue(NumericVector<Number> & current_solution)
{
  dof_id_type & dof_idx = _var.nodalDofIndex();
  _qp = 0;
  current_solution.set(dof_idx, computeQpValue());
}

bool
RichardsMultiphaseProblem::updateSolution(NumericVector<Number>& vec_solution, NumericVector<Number>& ghosted_solution)
{
  bool updatedSolution = false;  // this gets set to true if we needed to enforce the bound at any node

  unsigned int sys_num = getNonlinearSystem().number();

  // For parallel procs i believe that i have to use local_nodes_begin, rather than just nodes_begin
  // _mesh comes from SystemBase (_mesh = getNonlinearSystem().subproblem().mesh(), and subproblem is this object)
  MeshBase::node_iterator nit = _mesh.getMesh().local_nodes_begin();
  const MeshBase::node_iterator nend = _mesh.getMesh().local_nodes_end();

  for ( ; nit != nend; ++nit)
  {
    const Node & node = *(*nit);

    // dofs[0] is the dof number of the bounded variable at this node
    // dofs[1] is the dof number of the lower variable at this node
    std::vector<unsigned int> dofs(2);
    dofs[0] = node.dof_number(sys_num, _bounded_var_num, 0);
    dofs[1] = node.dof_number(sys_num, _lower_var_num, 0);

    // soln[0] is the value of the bounded variable at this node
    // soln[1] is the value of the lower variable at this node
    std::vector<Number> soln(2);
    vec_solution.get(dofs, soln);

    // do the bounding
    if (soln[0] < soln[1])
    {
      /*
      dof_id_type nd = node.id();
      Moose::out << "nd = " << nd << " dof_bounded = " << dofs[0] << " dof_lower = " << dofs[1] << "\n";
      Moose::out << " bounded_value = " << soln[0] << " lower_value = " << soln[1] << "\n";
      */
      vec_solution.set(dofs[0], soln[1]); // set the bounded variable equal to the lower value
      updatedSolution = true;
    }

  }

  // The above vec_solution.set calls potentially added "set" commands to a queue
  // The following actions the queue (doing MPI commands if necessary), so
  // vec_solution will actually be modified by this following command
  vec_solution.close();

  // if any proc updated the solution, all procs will know about it
  _communicator.max(updatedSolution);

  if (updatedSolution)
  {
    ghosted_solution = vec_solution;
    ghosted_solution.close();
  }

  return updatedSolution;

}

void
Assembly::setResidualBlock(NumericVector<Number> & residual, DenseVector<Number> & res_block, std::vector<dof_id_type> & dof_indices, Real scaling_factor)
{
  if (dof_indices.size() > 0)
  {
    std::vector<dof_id_type> di(dof_indices);
    _dof_map.constrain_element_vector(res_block, di, false);

    if (scaling_factor != 1.0)
    {
      _tmp_Re = res_block;
      _tmp_Re *= scaling_factor;
      for(unsigned int i=0; i<di.size(); i++)
        residual.set(di[i], _tmp_Re(i));
    }
    else
      for(unsigned int i=0; i<di.size(); i++)
        residual.set(di[i], res_block(i));
  }
}

void
dataLoad(std::istream & stream, NumericVector<Real> & v, void * /*context*/)
{
  numeric_index_type size = v.local_size();

  for (numeric_index_type i = v.first_local_index(); i < v.first_local_index() + size; i++)
  {
    Real r = 0;
    stream.read((char *) &r, sizeof(r));
    v.set(i, r);
  }

  v.close();
}

void AssembleOptimization::inequality_constraints (const NumericVector<Number> & X,
                                                   NumericVector<Number> & C_ineq,
                                                   OptimizationSystem & /*sys*/)
{
  C_ineq.zero();

  std::unique_ptr<NumericVector<Number>> X_localized =
    NumericVector<Number>::build(X.comm());
  X_localized->init(X.size(), false, SERIAL);
  X.localize(*X_localized);

  std::vector<Number> constraint_values(1);
  constraint_values[0] = (*X_localized)(200)*(*X_localized)(200) + (*X_localized)(201) - 5.;

  for (std::size_t i=0; i<constraint_values.size(); i++)
    if ((C_ineq.first_local_index() <= i) && (i < C_ineq.last_local_index()))
      C_ineq.set(i, constraint_values[i]);
}

void
NodalBC::computeResidual(NumericVector<Number> & residual)
{
    if (_var.isNodalDefined())
    {
        dof_id_type & dof_idx = _var.nodalDofIndex();
        _qp = 0;
        Real res = computeQpResidual();
        residual.set(dof_idx, res);

        if (_has_save_in)
        {
            Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
            for (unsigned int i=0; i<_save_in.size(); i++)
                _save_in[i]->sys().solution().set(_save_in[i]->nodalDofIndex(), res);
        }
    }
}

void AssembleOptimization::equality_constraints (const NumericVector<Number> & X,
                                                 NumericVector<Number> & C_eq,
                                                 OptimizationSystem & /*sys*/)
{
  C_eq.zero();

  UniquePtr<NumericVector<Number> > X_localized =
    NumericVector<Number>::build(X.comm());
  X_localized->init(X.size(), false, SERIAL);
  X.localize(*X_localized);

  std::vector<Number> constraint_values(3);
  constraint_values[0] = (*X_localized)(17);
  constraint_values[1] = (*X_localized)(23);
  constraint_values[2] = (*X_localized)(98) + (*X_localized)(185);

  for (unsigned int i=0; i<constraint_values.size(); i++)
    if ((C_eq.first_local_index() <= i) &&
        (i < C_eq.last_local_index()))
      C_eq.set(i, constraint_values[i]);
}

void
GrainTracker::swapSolutionValuesHelper(Node * curr_node, unsigned int curr_var_idx, unsigned int new_var_idx, NumericVector<Real> & solution, NumericVector<Real> & solution_old, NumericVector<Real> & solution_older)
{
  if (curr_node && curr_node->processor_id() == processor_id())
  {
    // Reinit the node so we can get and set values of the solution here
    _subproblem.reinitNode(curr_node, 0);

    // Swap the values from one variable to the other
    {
      VariableValue & value = _vars[curr_var_idx]->nodalSln();
      VariableValue & value_old = _vars[curr_var_idx]->nodalSlnOld();
      VariableValue & value_older = _vars[curr_var_idx]->nodalSlnOlder();

      // Copy Value from intersecting variable to new variable
      dof_id_type & dof_index = _vars[new_var_idx]->nodalDofIndex();

      // Set the only DOF for this variable on this node
      solution.set(dof_index, value[0]);
      solution_old.set(dof_index, value_old[0]);
      solution_older.set(dof_index, value_older[0]);
    }
    {
      VariableValue & value = _vars[new_var_idx]->nodalSln();
      VariableValue & value_old = _vars[new_var_idx]->nodalSlnOld();
      VariableValue & value_older = _vars[new_var_idx]->nodalSlnOlder();

      // Copy Value from variable to the intersecting variable
      dof_id_type & dof_index = _vars[curr_var_idx]->nodalDofIndex();

      // Set the only DOF for this variable on this node
      solution.set(dof_index, value[0]);
      solution_old.set(dof_index, value_old[0]);
      solution_older.set(dof_index, value_older[0]);
    }
  }
}

bool
FrictionalContactProblem::enforceRateConstraint(NumericVector<Number>& vec_solution, NumericVector<Number>& ghosted_solution)
{
    NonlinearSystem & nonlinear_sys = getNonlinearSystem();
    unsigned int dim = nonlinear_sys.subproblem().mesh().dimension();

    _displaced_problem->updateMesh(ghosted_solution, *_aux.currentSolution());

    MooseVariable * disp_x_var = &getVariable(0,_disp_x);
    MooseVariable * disp_y_var = &getVariable(0,_disp_y);
    MooseVariable * disp_z_var = NULL;
    if (dim == 3)
        disp_z_var = &getVariable(0,_disp_z);

    bool updatedSolution = false;

    if (getDisplacedProblem() && _interaction_params.size() > 0)
    {
        GeometricSearchData & displaced_geom_search_data = getDisplacedProblem()->geomSearchData();
        std::map<std::pair<unsigned int, unsigned int>, PenetrationLocator *> * penetration_locators = &displaced_geom_search_data._penetration_locators;

        for (std::map<std::pair<unsigned int, unsigned int>, PenetrationLocator *>::iterator plit = penetration_locators->begin();
                plit != penetration_locators->end();
                ++plit)
        {
            PenetrationLocator & pen_loc = *plit->second;
            std::set<dof_id_type> & has_penetrated = pen_loc._has_penetrated;

            bool frictional_contact_this_interaction = false;

            std::map<std::pair<int,int>,InteractionParams>::iterator ipit;
            std::pair<int,int> ms_pair(pen_loc._master_boundary,pen_loc._slave_boundary);
            ipit = _interaction_params.find(ms_pair);

            if (ipit != _interaction_params.end())
                frictional_contact_this_interaction = true;

            if (frictional_contact_this_interaction)
            {
                std::vector<dof_id_type> & slave_nodes = pen_loc._nearest_node._slave_nodes;

                for (unsigned int i=0; i<slave_nodes.size(); i++)
                {
                    dof_id_type slave_node_num = slave_nodes[i];

                    if (pen_loc._penetration_info[slave_node_num])
                    {
                        PenetrationInfo & info = *pen_loc._penetration_info[slave_node_num];

                        std::set<dof_id_type>::iterator hpit = has_penetrated.find(slave_node_num);

                        if (hpit != has_penetrated.end())
                        {
//              _console<<"Slave node: "<<slave_node_num<<std::endl;
                            const Node * node = info._node;

                            VectorValue<dof_id_type> solution_dofs(node->dof_number(nonlinear_sys.number(), disp_x_var->number(), 0),
                                                                   node->dof_number(nonlinear_sys.number(), disp_y_var->number(), 0),
                                                                   (disp_z_var ? node->dof_number(nonlinear_sys.number(), disp_z_var->number(), 0) : 0));

                            //Get old parametric coords from info
                            //get old element from info
                            //Apply current configuration to that element and find xyz coords of that point
                            //find xyz coords of current point from original coords and displacement vector
                            //calculate difference between those two point coordinates to get correction

//              std::vector<Point>points(1);
//              points[0] = info._starting_closest_point_ref;
//              Elem *side = (info._starting_elem->build_side(info._starting_side_num,false)).release();
//              FEBase &fe = *(pen_loc._fe[0]);
//              fe.reinit(side, &points);
//              RealVectorValue correction = starting_point[0] - current_coords;

//              RealVectorValue current_coords = *node;
//              RealVectorValue correction = info._closest_point - current_coords;

                            const Node & undisp_node = _mesh.node(node->id());
                            RealVectorValue solution = info._closest_point - undisp_node;

                            for (unsigned int i=0; i<dim; ++i)
                            {
//                  vec_solution.add(solution_dofs(i), correction(i));
                                vec_solution.set(solution_dofs(i), solution(i));
                            }
                            info._distance = 0.0;
                        }

                    }
                }
            }

            updatedSolution = true;
        }

        vec_solution.close();

        _communicator.max(updatedSolution);

        if (updatedSolution)
        {
//.........这里部分代码省略.........