C++ vector::clear方法代码示例

bool calc_errors(Hermes::vector<Solution* > left, Hermes::vector<Solution *> right, Hermes::vector<double> & err_abs, Hermes::vector<double> & norm_vals,
                 double & err_abs_total, double & norm_total, double & err_rel_total, Hermes::vector<ProjNormType> norms)
{
  bool default_norms = false;
  // Checks.
  if(left.size() != right.size())
    return false;
  if (norms != Hermes::vector<ProjNormType>())
  {
    if(left.size() != norms.size())
      return false;
  }
  else
    default_norms = true;

  // Zero the resulting Tuples.
  err_abs.clear();
  norm_vals.clear();

  // Zero the sums.
  err_abs_total = 0;
  norm_total = 0;
  err_rel_total = 0;

  // Calculation.
  for(unsigned int i = 0; i < left.size(); i++)
  {
    err_abs.push_back(calc_abs_error(left[i], right[i], default_norms ? HERMES_H1_NORM : norms[i]));
    norm_vals.push_back(calc_norm(right[i], default_norms ? HERMES_H1_NORM : norms[i]));
    err_abs_total += err_abs[i] * err_abs[i];
    norm_total += norm_vals[i] * norm_vals[i];
  }

  err_abs_total = sqrt(err_abs_total);
  norm_total = sqrt(norm_total);
  err_rel_total = err_abs_total / norm_total * 100.;

  // Everything went well, return appropriate flag.
  return true;
}

//.........这里部分代码省略.........
                    // emit signal
                    m_progressItemSolve->emitMessage(QObject::tr("Adaptivity rel. error (step: %2/%3, DOFs: %4/%5): %1%").
                                                     arg(error, 0, 'f', 3).
                                                     arg(i + 1).
                                                     arg(maxAdaptivitySteps).
                                                     arg(Space::get_num_dofs(space)).
                                                     arg(Space::get_num_dofs(spaceReference)), false, 1);
                    // add error to the list
                    m_progressItemSolve->addAdaptivityError(error, Space::get_num_dofs(space));

                    if (error < adaptivityTolerance || Space::get_num_dofs(space) >= adaptivityMaxDOFs)
                    {
                        break;
                    }
                    if (i != maxAdaptivitySteps-1) adaptivity.adapt(selector,
                                                                    Util::config()->threshold,
                                                                    Util::config()->strategy,
                                                                    Util::config()->meshRegularity);
                    actualAdaptivitySteps = i+1;
                }

                if (m_progressItemSolve->isCanceled())
                {
                    isError = true;
                    break;
                }

                // delete reference space
                for (int i = 0; i < spaceReference.size(); i++)
                {
                    delete spaceReference.at(i)->get_mesh();
                    delete spaceReference.at(i);
                }
                spaceReference.clear();
            }

            // clean up.
            delete solver;
            delete matrix;
            delete rhs;
        }

        // delete reference solution
        for (int i = 0; i < solutionReference.size(); i++)
            delete solutionReference.at(i);
        solutionReference.clear();

        // delete selector
        if (select) delete select;
        selector.clear();

        // timesteps
        if (!isError)
        {
            SparseMatrix *matrix = NULL;
            Vector *rhs = NULL;
            Solver *solver = NULL;

            // allocate dp for transient solution
            DiscreteProblem *dpTran = NULL;
            if (analysisType == AnalysisType_Transient)
            {
                // set up the solver, matrix, and rhs according to the solver selection.
                matrix = create_matrix(matrixSolver);
                rhs = create_vector(matrixSolver);
                solver = create_linear_solver(matrixSolver, matrix, rhs);

int main() 
{
  // Create space, set Dirichlet BC, enumerate basis functions.
  Space* space = new Space(A, B, NELEM, DIR_BC_LEFT, Hermes::vector<BCSpec *>(), P_INIT, NEQ, NEQ);
 
  // Enumerate basis functions, info for user.
  int ndof = Space::get_num_dofs(space);
  info("ndof: %d", ndof);

  // Initialize the weak formulation.
  WeakForm wf(2);
  wf.add_matrix_form(0, 0, jacobian_0_0);
  wf.add_matrix_form(0, 1, jacobian_0_1);
  wf.add_matrix_form(1, 0, jacobian_1_0);
  wf.add_matrix_form(1, 1, jacobian_1_1);
  wf.add_vector_form(0, residual_0);
  wf.add_vector_form(1, residual_1);

  // Initialize the FE problem.
  bool is_linear = false;
  DiscreteProblem *dp = new DiscreteProblem(&wf, space, is_linear);

  // Newton's loop.
  // Fill vector coeff_vec using dof and coeffs arrays in elements.
  double *coeff_vec = new double[Space::get_num_dofs(space)];
  get_coeff_vector(space, coeff_vec);

  // Set up the solver, matrix, and rhs according to the solver selection.
  SparseMatrix* matrix = create_matrix(matrix_solver);
  Vector* rhs = create_vector(matrix_solver);
  Solver* solver = create_linear_solver(matrix_solver, matrix, rhs);

  int it = 1;
  while (1) 
  {
    // Obtain the number of degrees of freedom.
    int ndof = Space::get_num_dofs(space);

    // Assemble the Jacobian matrix and residual vector.
    dp->assemble(coeff_vec, matrix, rhs);

    // Calculate the l2-norm of residual vector.
    double res_l2_norm = get_l2_norm(rhs);

    // Info for user.
    info("---- Newton iter %d, ndof %d, res. l2 norm %g", it, Space::get_num_dofs(space), res_l2_norm);

    // If l2 norm of the residual vector is within tolerance, then quit.
    // NOTE: at least one full iteration forced
    //       here because sometimes the initial
    //       residual on fine mesh is too small.
    if(res_l2_norm < NEWTON_TOL && it > 1) break;

    // Multiply the residual vector with -1 since the matrix 
    // equation reads J(Y^n) \deltaY^{n+1} = -F(Y^n).
    for(int i=0; i<ndof; i++) rhs->set(i, -rhs->get(i));

    // Solve the linear system.
    if(!solver->solve())
      error ("Matrix solver failed.\n");

    // Add \deltaY^{n+1} to Y^n.
    for (int i = 0; i < ndof; i++) coeff_vec[i] += solver->get_solution()[i];

    // If the maximum number of iteration has been reached, then quit.
    if (it >= NEWTON_MAX_ITER) error ("Newton method did not converge.");
    
    // Copy coefficients from vector y to elements.
    set_coeff_vector(coeff_vec, space);

    it++;
  }
  
  // Plot the solution.
  Linearizer l(space);
  l.plot_solution("solution.gp");

  // Plot the resulting space.
  space->plot("space.gp");

  // Cleanup.
  for(unsigned i = 0; i < DIR_BC_LEFT.size(); i++)
      delete DIR_BC_LEFT[i];
  DIR_BC_LEFT.clear();

  for(unsigned i = 0; i < DIR_BC_RIGHT.size(); i++)
      delete DIR_BC_RIGHT[i];
  DIR_BC_RIGHT.clear();

  delete matrix;
  delete rhs;
  delete solver;
  delete[] coeff_vec;
  delete dp;
  delete space;

  info("Done.");
  return 0;
}

    void WeakForm<Scalar>::get_stages(Hermes::vector<const Space<Scalar> *> spaces, Hermes::vector<Solution<Scalar> *>& u_ext,
      Hermes::vector<Stage<Scalar> >& stages, bool want_matrix, bool want_vector, bool one_stage) const
    {
      _F_;

      if (!want_matrix && !want_vector) return;

      unsigned int i;
      stages.clear();

      if (want_matrix || want_vector)
      {
        // This is because of linear problems where
        // matrix terms with the Dirichlet lift go to rhs.
        // Process volume matrix forms.
        for (i = 0; i < mfvol.size(); i++)
        {
          unsigned int ii = mfvol[i]->i, jj = mfvol[i]->j;
          Mesh* m1 = spaces[ii]->get_mesh();
          Mesh* m2 = spaces[jj]->get_mesh();
          Stage<Scalar>* s = find_stage(stages, ii, jj, m1, m2, mfvol[i]->ext, u_ext, one_stage);
          s->mfvol.push_back(mfvol[i]);
        }

        // Process surface matrix forms.
        for (i = 0; i < mfsurf.size(); i++)
        {
          unsigned int ii = mfsurf[i]->i, jj = mfsurf[i]->j;
          Mesh* m1 = spaces[ii]->get_mesh();
          Mesh* m2 = spaces[jj]->get_mesh();
          Stage<Scalar>* s = find_stage(stages, ii, jj, m1, m2, mfsurf[i]->ext, u_ext, one_stage);
          s->mfsurf.push_back(mfsurf[i]);
        }

        // Multi component forms.
        for (unsigned i = 0; i < mfvol_mc.size(); i++)
        {
          Mesh* the_one_mesh = spaces[mfvol_mc.at(i)->coordinates.at(0).first]->get_mesh();
          for(unsigned int form_i = 0; form_i < mfvol_mc.at(i)->coordinates.size(); form_i++)
          {
            if(spaces[mfvol_mc.at(i)->coordinates.at(form_i).first]->get_mesh()->get_seq() != the_one_mesh->get_seq())
              error("When using multi-component forms, the Meshes have to be identical.");
            if(spaces[mfvol_mc.at(i)->coordinates.at(form_i).second]->get_mesh()->get_seq() != the_one_mesh->get_seq())
              error("When using multi-component forms, the Meshes have to be identical.");
          }

          Stage<Scalar>* s = find_stage(stages, mfvol_mc.at(i)->coordinates, the_one_mesh, the_one_mesh, mfvol_mc[i]->ext, u_ext, one_stage);
          s->mfvol_mc.push_back(mfvol_mc[i]);
        }
        for (unsigned i = 0; i < mfsurf_mc.size(); i++)
        {
          Mesh* the_one_mesh = spaces[mfsurf_mc.at(i)->coordinates.at(0).first]->get_mesh();
          for(unsigned int form_i = 0; form_i < mfsurf_mc.at(i)->coordinates.size(); form_i++)
          {
            if(spaces[mfsurf_mc.at(i)->coordinates.at(form_i).first]->get_mesh()->get_seq() != the_one_mesh->get_seq())
              error("When using multi-component forms, the Meshes have to be identical.");
            if(spaces[mfsurf_mc.at(i)->coordinates.at(form_i).second]->get_mesh()->get_seq() != the_one_mesh->get_seq())
              error("When using multi-component forms, the Meshes have to be identical.");
          }

          Stage<Scalar>* s = find_stage(stages, mfsurf_mc.at(i)->coordinates, the_one_mesh, the_one_mesh, mfsurf_mc[i]->ext, u_ext, one_stage);
          s->mfsurf_mc.push_back(mfsurf_mc[i]);
        }
      }
      if (want_vector)
      {
        // Process volume vector forms.
        for (unsigned i = 0; i < vfvol.size(); i++)
        {
          unsigned int ii = vfvol[i]->i;
          Mesh *m = spaces[ii]->get_mesh();
          Stage<Scalar>*s = find_stage(stages, ii, ii, m, m, vfvol[i]->ext, u_ext, one_stage);
          s->vfvol.push_back(vfvol[i]);
        }

        // Process surface vector forms.
        for (unsigned i = 0; i < vfsurf.size(); i++)
        {
          unsigned int ii = vfsurf[i]->i;
          Mesh *m = spaces[ii]->get_mesh();
          Stage<Scalar>*s = find_stage(stages, ii, ii, m, m, vfsurf[i]->ext, u_ext, one_stage);
          s->vfsurf.push_back(vfsurf[i]);
        }

        // Multi component forms.
        for (unsigned i = 0; i < vfvol_mc.size(); i++)
        {
          Mesh* the_one_mesh = spaces[vfvol_mc.at(i)->coordinates.at(0)]->get_mesh();
          for(unsigned int form_i = 0; form_i < vfvol_mc.at(i)->coordinates.size(); form_i++)
            if(spaces[vfvol_mc.at(i)->coordinates.at(form_i)]->get_mesh()->get_seq() != the_one_mesh->get_seq())
              error("When using multi-component forms, the Meshes have to be identical.");

          Stage<Scalar>*s = find_stage(stages, vfvol_mc.at(i)->coordinates, the_one_mesh, the_one_mesh, vfvol_mc[i]->ext, u_ext, one_stage);
          s->vfvol_mc.push_back(vfvol_mc[i]);
        }

        for (unsigned i = 0; i < vfsurf_mc.size(); i++)
        {
          Mesh* the_one_mesh = spaces[vfsurf_mc.at(i)->coordinates.at(0)]->get_mesh();
          for(unsigned int form_i = 0; form_i < vfsurf_mc.at(i)->coordinates.size(); form_i++)
//.........这里部分代码省略.........