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

void
LStableDirk4::postStep(NumericVector<Number> & residual)
{
    // Error if _stage got messed up somehow.
    if (_stage > _n_stages)
        mooseError("LStableDirk4::postStep(): Member variable _stage can only have values 1-" << _n_stages << ".");

    // In the standard RK notation, the residual of stage 1 of s is given by:
    //
    // R := M*(Y_i - y_n)/dt - \sum_{j=1}^s a_{ij} * f(t_n + c_j*dt, Y_j) = 0
    //
    // where:
    // .) M is the mass matrix
    // .) Y_i is the stage solution
    // .) dt is the timestep, and is accounted for in the _Re_time residual.
    // .) f are the "non-time" residuals evaluated for a given stage solution.
    // .) The minus signs are already "baked in" to the residuals and so do not appear below.

    // Store this stage's non-time residual.  We are calling operator=
    // here, and that calls close().
    *_stage_residuals[_stage-1] = _Re_non_time;

    // Build up the residual for this stage.
    residual.add(1., _Re_time);
    for (unsigned int j = 0; j < _stage; ++j)
        residual.add(_a[_stage-1][j], *_stage_residuals[j]);
    residual.close();
}

void
NodalConstraint::computeResidual(NumericVector<Number> & residual)
{
  Real scaling_factor = _var.scalingFactor();
  _qp = 0;
  // master node
  dof_id_type & dof_idx = _var.nodalDofIndex();
  residual.add(dof_idx, scaling_factor * computeQpResidual(Moose::Master));
  // slave node
  dof_id_type & dof_idx_neighbor = _var.nodalDofIndexNeighbor();
  residual.add(dof_idx_neighbor, scaling_factor * computeQpResidual(Moose::Slave));
}

void
AStableDirk4::postResidual(NumericVector<Number> & residual)
{
  if (_t_step == 1 && _safe_start)
    _bootstrap_method->postResidual(residual);

  else
  {
    // Error if _stage got messed up somehow.
    if (_stage > 4)
      mooseError("AStableDirk4::postResidual(): Member variable _stage can only have values 1-4.");

    if (_stage < 4)
    {
      // In the standard RK notation, the residual of stage 1 of s is given by:
      //
      // R := M*(Y_i - y_n)/dt - \sum_{j=1}^s a_{ij} * f(t_n + c_j*dt, Y_j) = 0
      //
      // where:
      // .) M is the mass matrix
      // .) Y_i is the stage solution
      // .) dt is the timestep, and is accounted for in the _Re_time residual.
      // .) f are the "non-time" residuals evaluated for a given stage solution.
      // .) The minus signs are already "baked in" to the residuals and so do not appear below.

      // Store this stage's non-time residual.  We are calling operator=
      // here, and that calls close().
      *_stage_residuals[_stage - 1] = _Re_non_time;

      // Build up the residual for this stage.
      residual.add(1., _Re_time);
      for (unsigned int j = 0; j < _stage; ++j)
        residual.add(_a[_stage - 1][j], *_stage_residuals[j]);
      residual.close();
    }
    else
    {
      // The update step is a final solve:
      //
      // R := M*(y_{n+1} - y_n)/dt - \sum_{j=1}^s b_j * f(t_n + c_j*dt, Y_j) = 0
      //
      // We could potentially fold _b up into an extra row of _a and
      // just do one more stage, but I think handling it separately like
      // this is easier to understand and doesn't create too much code
      // repitition.
      residual.add(1., _Re_time);
      for (unsigned int j = 0; j < 3; ++j)
        residual.add(_b[j], *_stage_residuals[j]);
      residual.close();
    }
  }
}

void
ExplicitTVDRK2::postResidual(NumericVector<Number> & residual)
{
  if (_stage == 1)
  {
    // If postResidual() is called before solve(), _stage==1 and we don't
    // need to do anything.
  }
  else if (_stage == 2)
  {
    // In the standard RK notation, the stage 2 residual is given by:
    //
    // R := M*(Y_2 - y_n)/dt - f(t_n, Y_1) = 0
    //
    // where:
    // .) M is the mass matrix.
    // .) f(t_n, Y_1) is the residual we are currently computing,
    //    since this method is intended to be used with "implicit=false"
    //    kernels.
    // .) M*(Y_2 - y_n)/dt corresponds to the residual of the time kernels.
    // .) The minus signs are "baked in" to the non-time residuals, so
    //    they do not appear here.
    // .) The current non-time residual is saved for the next stage.
    _residual_old = _Re_non_time;
    _residual_old.close();

    residual.add(1.0, _Re_time);
    residual.add(1.0, _residual_old);
    residual.close();
  }
  else if (_stage == 3)
  {
    // In the standard RK notation, the update step residual is given by:
    //
    // R := M*(2*y_{n+1} - Y_2 - y_n)/(2*dt) - (1/2)*f(t_n+dt/2, Y_2) = 0
    //
    // where:
    // .) M is the mass matrix.
    // .) f(t_n+dt/2, Y_2) is the residual from stage 2.
    // .) The minus sign is already baked in to the non-time
    //    residuals, so it does not appear here.
    // .) Although this is an update step, we have to do a "solve"
    //    using the mass matrix.
    residual.add(1.0, _Re_time);
    residual.add(0.5, _Re_non_time);
    residual.close();
  }
  else
    mooseError(
        "ExplicitTVDRK2::postResidual(): _stage = ", _stage, ", only _stage = 1-3 is allowed.");
}

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

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

void
LStableDirk2::postStep(NumericVector<Number> & residual)
{
  if (_stage == 1)
  {
    // In the standard RK notation, the stage 1 residual is given by:
    //
    // R := (Y_1 - y_n)/dt - alpha*f(t_n + alpha*dt, Y_1) = 0
    //
    // where:
    // .) f(t_n + alpha*dt, Y_1) corresponds to the residuals of the
    //    non-time kernels.  We'll save this as "_residual_stage1" to use
    //    later.
    // .) (Y_1 - y_n)/dt corresponds to the residual of the time kernels.
    // .) The minus sign in front of alpha is already "baked in" to
    //    the non-time residuals, so it does not appear here.
    _residual_stage1 = _Re_non_time;
    _residual_stage1.close();

    residual.add(1., _Re_time);
    residual.add(_alpha, _residual_stage1);
    residual.close();
  }
  else if (_stage == 2)
  {
    // In the standard RK notation, the stage 2 residual is given by:
    //
    // R := (Y_2 - y_n)/dt - (1-alpha)*f(t_n + alpha*dt, Y_1) - alpha*f(t_n + dt, Y_2) = 0
    //
    // where:
    // .) f(t_n + alpha*dt, Y_1) has already been saved as _residual_stage1.
    // .) f(t_n + dt, Y_2) will now be saved as "_residual_stage2".
    // .) (Y_2 - y_n)/dt corresponds to the residual of the time kernels.
    // .) The minus signs are once again "baked in" to the non-time
    //    residuals, so they do not appear here.
    //
    // The solution at the end of stage 2, i.e. Y_2, is also the final solution.
    _residual_stage2 = _Re_non_time;
    _residual_stage2.close();

    residual.add(1., _Re_time);
    residual.add(1.-_alpha, _residual_stage1);
    residual.add(_alpha, _residual_stage2);
    residual.close();
  }
  else
    mooseError("LStableDirk2::postStep(): Member variable _stage can only have values 1 or 2.");
}

// Jacobian assembly function for the Laplace-Young system
void LaplaceYoung::postcheck (const NumericVector<Number> & old_soln,
                              NumericVector<Number> & search_direction,
                              NumericVector<Number> & new_soln,
                              bool & /*changed_search_direction*/,
                              bool & changed_new_soln,
                              NonlinearImplicitSystem & /*S*/)
{
  // Back up along the search direction by some amount.  Since Newton
  // already works well for this problem, the only affect of this is
  // to degrade the rate of convergence.
  //
  // The minus sign is due to the sign of the "search_direction"
  // vector which comes in from the Newton solve. The RHS of the
  // nonlinear system, i.e. the residual, is generally not multiplied
  // by -1, so the solution vector, i.e. the search_direction, has a
  // factor -1.
  if (_gamma != 1.0)
    {
      new_soln = old_soln;
      new_soln.add(-_gamma, search_direction);
      changed_new_soln = true;
    }
  else
    changed_new_soln = false;
}

void AssembleOptimization::gradient (const NumericVector<Number> & soln,
                                     NumericVector<Number> & grad_f,
                                     OptimizationSystem & /*sys*/)
{
  grad_f.zero();

  A_matrix->vector_mult(grad_f, soln);
  grad_f.add(-1, *F_vector);
}

void
ImplicitMidpoint::postStep(NumericVector<Number> & residual)
{
  if (_stage == 1)
  {
    // In the standard RK notation, the stage 1 residual is given by:
    //
    // R := M*(Y_1 - y_n)/dt - (1/2)*f(t_n + dt/2, Y_1) = 0
    //
    // where:
    // .) M is the mass matrix
    // .) f(t_n + dt/2, Y_1) is saved in _residual_stage1
    // .) The minus sign is baked in to the non-time residuals, so it does not appear here.
    _residual_stage1 = _Re_non_time;
    _residual_stage1.close();

    residual.add(1., _Re_time);
    residual.add(0.5, _Re_non_time);
    residual.close();
  }
  else if (_stage == 2)
  {
    // The update step.  In the standard RK notation, the update
    // residual is given by:
    //
    // R := M*(y_{n+1} - y_n)/dt - f(t_n + dt/2, Y_1) = 0
    //
    // where:
    // .) M is the mass matrix.
    // .) f(t_n + dt/2, Y_1) is the residual from stage 1, it has already
    //    been saved as _residual_stage1.
    // .) The minus signs are "baked in" to the non-time residuals, so
    //    they do not appear here.
    residual.add(1., _Re_time);
    residual.add(1., _residual_stage1);
    residual.close();
  }
  else
    mooseError("ImplicitMidpoint::postStep(): _stage = " << _stage << ", only _stage = 1, 2 is allowed.");
}

// Jacobian assembly function for the Laplace-Young system
void LaplaceYoung::postcheck (const NumericVector<Number> & /*old_soln*/,
                              NumericVector<Number> & search_direction,
                              NumericVector<Number> & new_soln,
                              bool & /*changed_search_direction*/,
                              bool & changed_new_soln,
                              NonlinearImplicitSystem & /*S*/)
{
  // Back up along the search direction by some amount.  Since Newton
  // already works well for this problem, the only affect of this is
  // to degrade the rate of convergence.
  new_soln.add(_gamma - 1., search_direction);
  changed_new_soln = true;
}

//.........这里部分代码省略.........
      unsigned ndofP = msh->GetElementDofNumber(iel, solPType);
      unsigned ndofD = msh->GetElementDofNumber(iel, solDType);
      // resize

      phiV.resize(ndofV);
      gradphiV.resize(ndofV * dim);
      nablaphiV.resize(ndofV * (3 * (dim - 1) + !(dim - 1)));

      solP.resize(ndofP);
      for(int d = 0; d < dim; d++) {
        solV[d].resize(ndofV);
        x0[d].resize(ndofD);
        x[d].resize(ndofD);
      }
      // get local to global mappings
      for(unsigned i = 0; i < ndofD; i++) {
        unsigned idof = msh->GetSolutionDof(i, iel, solDType);
        for(unsigned d = 0; d < dim; d++) {
          x0[d][i] = (*msh->_topology->_Sol[d])(idof);

          x[d][i] = (*msh->_topology->_Sol[d])(idof) +
                    (*solution->_Sol[solDIndex[d]])(idof);
        }
      }

      for(unsigned i = 0; i < ndofV; i++) {
        unsigned idof = msh->GetSolutionDof(i, iel, solVType);    // global to global mapping between solution node and solution dof
        for(unsigned  d = 0; d < dim; d++) {
          solV[d][i] = (*solution->_Sol[solVIndex[d]])(idof);      // global extraction and local storage for the solution
        }
      }



      for(unsigned i = 0; i < ndofP; i++) {
        unsigned idof = msh->GetSolutionDof(i, iel, solPType);
        solP[i] = (*solution->_Sol[solPIndex])(idof);
      }


      for(unsigned ig = 0; ig < mlSol._mlMesh->_finiteElement[ielt][solVType]->GetGaussPointNumber(); ig++) {
        // *** get Jacobian and test function and test function derivatives ***
        msh->_finiteElement[ielt][solVType]->Jacobian(x0, ig, weight0, phiV, gradphiV, nablaphiV);
        msh->_finiteElement[ielt][solVType]->Jacobian(x, ig, weight, phiV, gradphiV, nablaphiV);
        phiP = msh->_finiteElement[ielt][solPType]->GetPhi(ig);

        vol0->add(iproc, weight0);
        vol->add(iproc, weight);

        std::vector < double> SolV2(dim, 0.);
        for(unsigned i = 0; i < ndofV; i++) {
          for(unsigned d = 0; d < dim; d++) {
            SolV2[d] += solV[d][i] * phiV[i];
          }
        }

        double V2 = 0.;
        for(unsigned d = 0; d < dim; d++) {
          V2 += SolV2[d] * SolV2[d];
        }
        v2->add(iproc, V2 * weight);

        double P2 = 0;
        for(unsigned i = 0; i < ndofP; i++) {
          P2 += solP[i] * phiP[i];
        }
        P2 *= P2;
        p2->add(iproc, P2 * weight);
      }
    }
  }

  p2->close();
  v2->close();
  vol0->close();
  vol->close();

  double p2_l2 = p2->l1_norm();
  double v2_l2 = v2->l1_norm();
  double VOL0 = vol0->l1_norm();
  double VOL = vol->l1_norm();

  std::cout.precision(14);
  std::scientific;
  std::cout << " vol0 = " << VOL0 << std::endl;
  std::cout << " vol = " << VOL << std::endl;
  std::cout << " (vol-vol0)/vol0 = " << (VOL - VOL0) / VOL0 << std::endl;
  std::cout << " p_l2 norm / vol = " << sqrt(p2_l2 / VOL)  << std::endl;
  std::cout << " v_l2 norm / vol = " << sqrt(v2_l2 / VOL)  << std::endl;

  data[1] = (VOL - VOL0) / VOL0;
  data[2] = VOL;
  data[3] = sqrt(p2_l2 / VOL);
  data[4] = sqrt(v2_l2 / VOL);

  delete p2;
  delete v2;
  delete vol;

}

void
FrictionalContactProblem::applySlip(NumericVector<Number>& vec_solution,
                                    NumericVector<Number>& ghosted_solution,
                                    std::vector<SlipData> & iterative_slip)
{
    NonlinearSystem & nonlinear_sys = getNonlinearSystem();
    unsigned int dim = nonlinear_sys.subproblem().mesh().dimension();
    AuxiliarySystem & aux_sys = getAuxiliarySystem();
    NumericVector<Number> & aux_solution = aux_sys.solution();

    MooseVariable * disp_x_var = &getVariable(0,_disp_x);
    MooseVariable * disp_y_var = &getVariable(0,_disp_y);
    MooseVariable * disp_z_var = NULL;
    MooseVariable * inc_slip_x_var = &getVariable(0,_inc_slip_x);
    MooseVariable * inc_slip_y_var = &getVariable(0,_inc_slip_y);
    MooseVariable * inc_slip_z_var = NULL;
    if (dim == 3)
    {
        disp_z_var = &getVariable(0,_disp_z);
        inc_slip_z_var = &getVariable(0,_inc_slip_z);
    }

    MooseVariable * disp_var;
    MooseVariable * inc_slip_var;

    for (unsigned int iislip=0; iislip<iterative_slip.size(); ++iislip)
    {
        const Node * node = iterative_slip[iislip]._node;
        const unsigned int dof = iterative_slip[iislip]._dof;
        const Real slip = iterative_slip[iislip]._slip;

        if (dof == 0)
        {
            disp_var = disp_x_var;
            inc_slip_var = inc_slip_x_var;
        }
        else if (dof == 1)
        {
            disp_var = disp_y_var;
            inc_slip_var = inc_slip_y_var;
        }
        else
        {
            disp_var = disp_z_var;
            inc_slip_var = inc_slip_z_var;
        }

        dof_id_type
        solution_dof = node->dof_number(nonlinear_sys.number(), disp_var->number(), 0),
        inc_slip_dof = node->dof_number(aux_sys.number(), inc_slip_var->number(), 0);

        vec_solution.add(solution_dof, slip);
        aux_solution.add(inc_slip_dof, slip);
    }

    aux_solution.close();
    vec_solution.close();
    unsigned int num_slipping_nodes = iterative_slip.size();
    _communicator.sum(num_slipping_nodes);
    if (num_slipping_nodes > 0)
    {
        ghosted_solution = vec_solution;
        ghosted_solution.close();

        updateContactPoints(ghosted_solution,false);

        enforceRateConstraint(vec_solution, ghosted_solution);
    }
}

Real NewtonSolver::line_search(Real tol,
                               Real last_residual,
                               Real &current_residual,
                               NumericVector<Number> &newton_iterate,
                               const NumericVector<Number> &linear_solution)
{
  // Take a full step if we got a residual reduction or if we
  // aren't substepping
  if ((current_residual < last_residual) ||
      (!require_residual_reduction &&
       (!require_finite_residual || !libmesh_isnan(current_residual))))
    return 1.;

  // The residual vector
  NumericVector<Number> &rhs = *(_system.rhs);

  Real ax = 0.;  // First abscissa, don't take negative steps
  Real cx = 1.;  // Second abscissa, don't extrapolate steps

  // Find bx, a step length that gives lower residual than ax or cx
  Real bx = 1.;

  while (libmesh_isnan(current_residual) ||
         (current_residual > last_residual &&
          require_residual_reduction))
    {
      // Reduce step size to 1/2, 1/4, etc.
      Real substepdivision;
      if (brent_line_search && !libmesh_isnan(current_residual))
        {
          substepdivision = std::min(0.5, last_residual/current_residual);
          substepdivision = std::max(substepdivision, tol*2.);
        }
      else
        substepdivision = 0.5;

      newton_iterate.add (bx * (1.-substepdivision),
                          linear_solution);
      newton_iterate.close();
      bx *= substepdivision;
      if (verbose)
        libMesh::out << "  Shrinking Newton step to "
                     << bx << std::endl;

      // Check residual with fractional Newton step
      _system.assembly (true, false);

      rhs.close();
      current_residual = rhs.l2_norm();
      if (verbose)
        libMesh::out << "  Current Residual: "
                     << current_residual << std::endl;

      if (bx/2. < minsteplength &&
          (libmesh_isnan(current_residual) ||
           (current_residual > last_residual)))
        {
          libMesh::out << "Inexact Newton step FAILED at step "
                       << _outer_iterations << std::endl;

          if (!continue_after_backtrack_failure)
            {
              libmesh_convergence_failure();
            }
          else
            {
              libMesh::out << "Continuing anyway ..." << std::endl;
              _solve_result = DiffSolver::DIVERGED_BACKTRACKING_FAILURE;
              return bx;
            }
        }
    } // end while (current_residual > last_residual)

  // Now return that reduced-residual step, or  use Brent's method to
  // find a more optimal step.

  if (!brent_line_search)
    return bx;

  // Brent's method adapted from Numerical Recipes in C, ch. 10.2
  Real e = 0.;

  Real x = bx, w = bx, v = bx;

  // Residuals at bx
  Real fx = current_residual,
    fw = current_residual,
    fv = current_residual;

  // Max iterations for Brent's method loop
  const unsigned int max_i = 20;

  // for golden ratio steps
  const Real golden_ratio = 1.-(std::sqrt(5.)-1.)/2.;

  for (unsigned int i=1; i <= max_i; i++)
    {
      Real xm = (ax+cx)*0.5;
      Real tol1 = tol * std::abs(x) + tol*tol;
      Real tol2 = 2.0 * tol1;
//.........这里部分代码省略.........

void TensorShellMatrix<T>::vector_mult_add (NumericVector<T>& dest,
                                            const NumericVector<T>& arg) const
{
  dest.add(_w.dot(arg),_v);
}