int main(int argc, char* argv[])
{
  // Load the mesh.
  Mesh mesh_whole_domain, mesh_with_hole;
  Hermes::vector<Mesh*> meshes (&mesh_whole_domain, &mesh_with_hole);
  MeshReaderH2DXML mloader;
  mloader.load("domain.xml", meshes);

  // Temperature mesh: Initial uniform mesh refinements in graphite.
  meshes[0]->refine_by_criterion(element_in_graphite, INIT_REF_NUM_TEMPERATURE_GRAPHITE);

  // Temperature mesh: Initial uniform mesh refinements in fluid.
  meshes[0]->refine_by_criterion(element_in_fluid, INIT_REF_NUM_TEMPERATURE_FLUID);

  // Fluid mesh: Initial uniform mesh refinements.
  for(int i = 0; i < INIT_REF_NUM_FLUID; i++)
    meshes[1]->refine_all_elements();

  // Initial refinements towards boundary of graphite.
  for(unsigned int meshes_i = 0; meshes_i < meshes.size(); meshes_i++)
    meshes[meshes_i]->refine_towards_boundary("Inner Wall", INIT_REF_NUM_BDY_GRAPHITE);

  // Initial refinements towards the top and bottom edges.
  for(unsigned int meshes_i = 0; meshes_i < meshes.size(); meshes_i++)
    meshes[meshes_i]->refine_towards_boundary("Outer Wall", INIT_REF_NUM_BDY_WALL);

  /* View both meshes. */
  MeshView m1("Mesh for temperature"), m2("Mesh for fluid");
  m1.show(&mesh_whole_domain);
  m2.show(&mesh_with_hole);

  // Initialize boundary conditions.
  EssentialBCNonConst bc_inlet_vel_x("Inlet", VEL_INLET, H, STARTUP_TIME);
  DefaultEssentialBCConst<double> bc_other_vel_x(Hermes::vector<std::string>("Outer Wall", "Inner Wall"), 0.0);
  EssentialBCs<double> bcs_vel_x(Hermes::vector<EssentialBoundaryCondition<double> *>(&bc_inlet_vel_x, &bc_other_vel_x));
  DefaultEssentialBCConst<double> bc_vel_y(Hermes::vector<std::string>("Inlet", "Outer Wall", "Inner Wall"), 0.0);
  EssentialBCs<double> bcs_vel_y(&bc_vel_y);
  EssentialBCs<double> bcs_pressure;
  DefaultEssentialBCConst<double> bc_temperature(Hermes::vector<std::string>("Outer Wall", "Inlet"), 20.0);
  EssentialBCs<double> bcs_temperature(&bc_temperature);

  // Spaces for velocity components, pressure and temperature.
  H1Space<double> xvel_space(&mesh_with_hole, &bcs_vel_x, P_INIT_VEL);
  H1Space<double> yvel_space(&mesh_with_hole, &bcs_vel_y, P_INIT_VEL);
#ifdef PRESSURE_IN_L2
  L2Space<double> p_space(&mesh_with_hole, P_INIT_PRESSURE);
#else
  H1Space<double> p_space(&mesh_with_hole, &bcs_pressure, P_INIT_PRESSURE);
#endif
  H1Space<double> temperature_space(&mesh_whole_domain, &bcs_temperature, P_INIT_TEMPERATURE);
  Hermes::vector<Space<double> *> all_spaces(&xvel_space, 
      &yvel_space, &p_space, &temperature_space);
  Hermes::vector<const Space<double> *> all_spaces_const(&xvel_space, 
      &yvel_space, &p_space, &temperature_space);

  // Calculate and report the number of degrees of freedom.
  int ndof = Space<double>::get_num_dofs(Hermes::vector<const Space<double> *>(&xvel_space, 
      &yvel_space, &p_space, &temperature_space));
  info("ndof = %d.", ndof);

  // Define projection norms.
  ProjNormType vel_proj_norm = HERMES_H1_NORM;
#ifdef PRESSURE_IN_L2
  ProjNormType p_proj_norm = HERMES_L2_NORM;
#else
  ProjNormType p_proj_norm = HERMES_H1_NORM;
#endif
  ProjNormType temperature_proj_norm = HERMES_H1_NORM;
  Hermes::vector<ProjNormType> all_proj_norms = Hermes::vector<ProjNormType>(vel_proj_norm, 
      vel_proj_norm, p_proj_norm, temperature_proj_norm);

  // Initial conditions and such.
  info("Setting initial conditions.");
  ZeroSolution xvel_prev_time(&mesh_with_hole), yvel_prev_time(&mesh_with_hole), p_prev_time(&mesh_with_hole);
  CustomInitialConditionTemperature temperature_init_cond(&mesh_whole_domain, HOLE_MID_X, HOLE_MID_Y, 
      0.5*OBSTACLE_DIAMETER, TEMPERATURE_INIT_FLUID, TEMPERATURE_INIT_GRAPHITE); 
  Solution<double> temperature_prev_time;
  Hermes::vector<Solution<double> *> all_solutions = Hermes::vector<Solution<double> *>(&xvel_prev_time, 
      &yvel_prev_time, &p_prev_time, &temperature_prev_time);
  Hermes::vector<MeshFunction<double> *> all_meshfns = Hermes::vector<MeshFunction<double> *>(&xvel_prev_time, 
      &yvel_prev_time, &p_prev_time, &temperature_init_cond);

  // Project all initial conditions on their FE spaces to obtain aninitial
  // coefficient vector for the Newton's method. We use local projection
  // to avoid oscillations in temperature on the graphite-fluid interface
  // FIXME - currently the LocalProjection only does the lowest-order part (linear
  // interpolation) at the moment. Higher-order part needs to be added.
  double* coeff_vec = new double[ndof];
  info("Projecting initial condition to obtain initial vector for the Newton's method.");
  //OGProjection<double>::project_global(all_spaces, all_meshfns, coeff_vec, matrix_solver, all_proj_norms);
  LocalProjection<double>::project_local(all_spaces_const, all_meshfns, coeff_vec, matrix_solver, all_proj_norms);

  // Translate the solution vector back to Solutions. This is needed to replace
  // the discontinuous initial condition for temperature_prev_time with its projection.
  Solution<double>::vector_to_solutions(coeff_vec, all_spaces_const, all_solutions);

  // Calculate Reynolds number.
  double reynolds_number = VEL_INLET * OBSTACLE_DIAMETER / KINEMATIC_VISCOSITY_FLUID;
  info("RE = %g", reynolds_number);
  if (reynolds_number < 1e-8) error("Re == 0 will not work - the equations use 1/Re.");
//.........这里部分代码省略.........

    bool Adapt<Scalar>::adapt(Hermes::vector<RefinementSelectors::Selector<Scalar> *> refinement_selectors, double thr, int strat,
      int regularize, double to_be_processed)
    {
      error_if(!have_errors, "element errors have to be calculated first, call Adapt<Scalar>::calc_err_est().");
      error_if(refinement_selectors == Hermes::vector<RefinementSelectors::Selector<Scalar> *>(), "selector not provided");
      if (spaces.size() != refinement_selectors.size()) error("Wrong number of refinement selectors.");
      Hermes::TimePeriod cpu_time;

      //get meshes
      int max_id = -1;
      Mesh* meshes[H2D_MAX_COMPONENTS];
      for (int j = 0; j < this->num; j++) 
      {
        meshes[j] = this->spaces[j]->get_mesh();
        if (rsln[j] != NULL) 
        {
          rsln[j]->set_quad_2d(&g_quad_2d_std);
          rsln[j]->enable_transform(false);
        }
        if (meshes[j]->get_max_element_id() > max_id)
          max_id = meshes[j]->get_max_element_id();
      }

      //reset element refinement info
      int** idx = new int*[max_id];
      for(int i = 0; i < max_id; i++)
        idx[i] = new int[num];

      for(int j = 0; j < max_id; j++)
        for(int l = 0; l < this->num; l++)
          idx[j][l] = -1; // element not refined

      double err0_squared = 1000.0;
      double processed_error_squared = 0.0;

      std::vector<ElementToRefine> elem_inx_to_proc; //list of indices of elements that are going to be processed
      elem_inx_to_proc.reserve(num_act_elems);

      //adaptivity loop
      double error_squared_threshod = -1; //an error threshold that breaks the adaptivity loop in a case of strategy 1
      int num_exam_elem = 0; //a number of examined elements
      int num_ignored_elem = 0; //a number of ignored elements
      int num_not_changed = 0; //a number of element that were not changed
      int num_priority_elem = 0; //a number of elements that were processed using priority queue

      bool first_regular_element = true; //true if first regular element was not processed yet
      int inx_regular_element = 0;
      while (inx_regular_element < num_act_elems || !priority_queue.empty())
      {
        int id, comp, inx_element;

        //get element identification
        if (priority_queue.empty()) 
        {
          id = regular_queue[inx_regular_element].id;
          comp = regular_queue[inx_regular_element].comp;
          inx_element = inx_regular_element;
          inx_regular_element++;
        }
        else 
        {
          id = priority_queue.front().id;
          comp = priority_queue.front().comp;
          inx_element = -1;
          priority_queue.pop();
          num_priority_elem++;
        }
        num_exam_elem++;

        //get info linked with the element
        double err_squared = errors[comp][id];
        Mesh* mesh = meshes[comp];
        Element* e = mesh->get_element(id);

        if (!should_ignore_element(inx_element, mesh, e)) 
        {
          //check if adaptivity loop should end
          if (inx_element >= 0) 
          {
            //prepare error threshold for strategy 1
            if (first_regular_element) 
            {
              error_squared_threshod = thr * err_squared;
              first_regular_element = false;
            }

            // first refinement strategy:
            // refine elements until prescribed amount of error is processed
            // if more elements have similar error refine all to keep the mesh symmetric
            if ((strat == 0) && (processed_error_squared > sqrt(thr) * errors_squared_sum)
              && fabs((err_squared - err0_squared)/err0_squared) > 1e-3) break;

            // second refinement strategy:
            // refine all elements whose error is bigger than some portion of maximal error
            if ((strat == 1) && (err_squared < error_squared_threshod)) break;

            if ((strat == 2) && (err_squared < thr)) break;

            if ((strat == 3) &&
              ( (err_squared < error_squared_threshod) ||
//.........这里部分代码省略.........

    void NeighborSearch<Scalar>::clear_initial_sub_idx()
    {
      if(neighborhood_type != H2D_DG_GO_DOWN)
        return;
      // Obtain the transformations sequence.
      Hermes::vector<unsigned int> transformations = get_transforms(original_central_el_transform);

      Hermes::vector<unsigned int> updated_transformations;
      for(int i = 0; i < transformations.size(); i++)
      {
        if(! ((active_edge == 0 && transformations[i] == 4) || (active_edge == 1 && transformations[i] == 7) || (active_edge == 2 && transformations[i] == 5) || (active_edge == 3 && transformations[i] == 6)) )
        {
          if(active_edge == 0 && transformations[i] == 6)
            updated_transformations.push_back(0);
          else if(active_edge == 0 && transformations[i] == 7)
            updated_transformations.push_back(1);
          else if(active_edge == 1 && transformations[i] == 4)
            updated_transformations.push_back(1);
          else if(active_edge == 1 && transformations[i] == 5)
            updated_transformations.push_back(2);
          else if(active_edge == 2 && transformations[i] == 6)
            updated_transformations.push_back(3);
          else if(active_edge == 2 && transformations[i] == 7)
            updated_transformations.push_back(2);
          else if(active_edge == 3 && transformations[i] == 4)
            updated_transformations.push_back(0);
          else if(active_edge == 3 && transformations[i] == 5)
            updated_transformations.push_back(3);
          else
            updated_transformations.push_back(transformations[i]);
        }
      }

      // Test for active element.
      if(updated_transformations.size() == 0)
        return;

      for(unsigned int i = 0; i < n_neighbors; i++)
      {
        // Find the index where the additional subelement mapping (on top of the initial one from assembling) starts.
        unsigned int j = 0;
        // Note that we do not have to test if central_transformations is empty or how long it is, because it has to be
        // longer than transformations (and that is tested).
        // Also the function compatible_transformations() does not have to be used, as now the array central_transformations
        // has been adjusted so that it contains the array transformations.
        while(central_transformations.get(i)->transf[j] == updated_transformations[j])
          if(++j > updated_transformations.size() - 1)
            break;
        if(j > central_transformations.get(i)->num_levels)
          j = central_transformations.get(i)->num_levels;

        for(unsigned int level = central_transformations.get(i)->num_levels; level < updated_transformations.size(); level++)
        {
          if(!neighbor_transformations.present(i))
            neighbor_transformations.add(new Transformations, i);

          Transformations* neighbor_transforms = neighbor_transformations.get(i);

          // Triangles.
          if(central_el->get_mode() == HERMES_MODE_TRIANGLE)
            if((active_edge == 0 && updated_transformations[level] == 0) ||
              (active_edge == 1 && updated_transformations[level] == 1) ||
              (active_edge == 2 && updated_transformations[level] == 2))
              neighbor_transforms->transf[neighbor_transforms->num_levels++] = (!neighbor_edge.orientation ? neighbor_edge.local_num_of_edge : (neighbor_edge.local_num_of_edge + 1) % 3);
            else
              neighbor_transforms->transf[neighbor_transforms->num_levels++] = (neighbor_edge.orientation ? neighbor_edge.local_num_of_edge : (neighbor_edge.local_num_of_edge + 1) % 3);
          // Quads.
          else
            if((active_edge == 0 && (updated_transformations[level] == 0 || updated_transformations[level] == 6)) ||
              (active_edge == 1 && (updated_transformations[level] == 1 || updated_transformations[level] == 4)) ||
              (active_edge == 2 && (updated_transformations[level] == 2 || updated_transformations[level] == 7)) ||
              (active_edge == 3 && (updated_transformations[level] == 3 || updated_transformations[level] == 5)))
              neighbor_transforms->transf[neighbor_transforms->num_levels++] = (!neighbor_edge.orientation ? neighbor_edge.local_num_of_edge : (neighbor_edge.local_num_of_edge + 1) % 4);
            else if((active_edge == 0 && (updated_transformations[level] == 1 || updated_transformations[level] == 7)) ||
              (active_edge == 1 && (updated_transformations[level] == 2 || updated_transformations[level] == 5)) ||
              (active_edge == 2 && (updated_transformations[level] == 3 || updated_transformations[level] == 6)) ||
              (active_edge == 3 && (updated_transformations[level] == 0 || updated_transformations[level] == 4)))
              neighbor_transforms->transf[neighbor_transforms->num_levels++] = (neighbor_edge.orientation ? neighbor_edge.local_num_of_edge : (neighbor_edge.local_num_of_edge + 1) % 4);
        }

        central_transformations.get(i)->strip_initial_transformations(j);
      }
    }

  bool select_refinement(Element* element, int order, MeshFunction<complex>* rsln, ElementToRefine& refinement)
  {
    switch(strategy)
    {
    case(noSelectionH):
      {
        refinement.split = H2D_REFINEMENT_H;
        refinement.best_refinement_polynomial_order_type[H2D_REFINEMENT_H][0] = 
          refinement.best_refinement_polynomial_order_type[H2D_REFINEMENT_H][1] = 
          refinement.best_refinement_polynomial_order_type[H2D_REFINEMENT_H][2] = 
          refinement.best_refinement_polynomial_order_type[H2D_REFINEMENT_H][3] = 
          order;
        ElementToRefine::copy_orders(refinement.refinement_polynomial_order, refinement.best_refinement_polynomial_order_type[H2D_REFINEMENT_H]);
        return true;
      }
      break;
    case(noSelectionHP):
      {
        int max_allowed_order = this->max_order;
        if(this->max_order == H2DRS_DEFAULT_ORDER)
          max_allowed_order = H2DRS_MAX_ORDER;
        int order_h = H2D_GET_H_ORDER(order), order_v = H2D_GET_V_ORDER(order);
        int increased_order_h = std::min(max_allowed_order, order_h + 1), increased_order_v = std::min(max_allowed_order, order_v + 1);
        int increased_order;
        if(element->is_triangle())
          increased_order = refinement.best_refinement_polynomial_order_type[H2D_REFINEMENT_H][0] = H2D_MAKE_QUAD_ORDER(increased_order_h, increased_order_h);
        else
          increased_order = refinement.best_refinement_polynomial_order_type[H2D_REFINEMENT_H][0] = H2D_MAKE_QUAD_ORDER(increased_order_h, increased_order_v);

        refinement.split = H2D_REFINEMENT_H;
        refinement.best_refinement_polynomial_order_type[H2D_REFINEMENT_H][0] = 
          refinement.best_refinement_polynomial_order_type[H2D_REFINEMENT_H][1] = 
          refinement.best_refinement_polynomial_order_type[H2D_REFINEMENT_H][2] = 
          refinement.best_refinement_polynomial_order_type[H2D_REFINEMENT_H][3] = 
          increased_order;
        ElementToRefine::copy_orders(refinement.refinement_polynomial_order, refinement.best_refinement_polynomial_order_type[H2D_REFINEMENT_H]);
        return true;
      }
      case(hXORpSelectionBasedOnError):
      {
        //make an uniform order in a case of a triangle
        int order_h = H2D_GET_H_ORDER(order), order_v = H2D_GET_V_ORDER(order);

        int current_min_order, current_max_order;
        this->get_current_order_range(element, current_min_order, current_max_order);

        if(current_max_order < std::max(order_h, order_v))
          current_max_order = std::max(order_h, order_v);

        int last_order_h = std::min(current_max_order, order_h + 1), last_order_v = std::min(current_max_order, order_v + 1);
        int last_order = H2D_MAKE_QUAD_ORDER(last_order_h, last_order_v);

        //build candidates.
        Hermes::vector<Cand> candidates;
        candidates.push_back(Cand(H2D_REFINEMENT_P, last_order));
        candidates.push_back(Cand(H2D_REFINEMENT_H, order, order, order, order));
        
        this->evaluate_cands_error(candidates, element, rsln);
        
        Cand* best_candidate = (candidates[0].error < candidates[1].error) ? &candidates[0] : &candidates[1];
        Cand* best_candidates_specific_type[4];
        best_candidates_specific_type[H2D_REFINEMENT_P] = &candidates[0];
        best_candidates_specific_type[H2D_REFINEMENT_H] = &candidates[1];
        best_candidates_specific_type[2] = NULL;
        best_candidates_specific_type[3] = NULL;

        //copy result to output
        refinement.split = best_candidate->split;
        ElementToRefine::copy_orders(refinement.refinement_polynomial_order, best_candidate->p);
        for(int i = 0; i < 4; i++)
          if(best_candidates_specific_type[i] != NULL)
            ElementToRefine::copy_orders(refinement.best_refinement_polynomial_order_type[i], best_candidates_specific_type[i]->p);

        ElementToRefine::copy_errors(refinement.errors, best_candidate->errors);

        //modify orders in a case of a triangle such that order_v is zero
        if(element->is_triangle())
          for(int i = 0; i < H2D_MAX_ELEMENT_SONS; i++)
            refinement.refinement_polynomial_order[i] = H2D_MAKE_QUAD_ORDER(H2D_GET_H_ORDER(refinement.refinement_polynomial_order[i]), 0);

        return true;
      }
    default:
      H1ProjBasedSelector<complex>::select_refinement(element, order, rsln, refinement);
      return true;
      break;
    }
  }

double KellyTypeAdapt::eval_interface_estimator(KellyTypeAdapt::ErrorEstimatorForm* err_est_form,
                                                RefMap *rm, SurfPos* surf_pos,
                                                LightArray<NeighborSearch*>& neighbor_searches, int neighbor_index)
{
  NeighborSearch* nbs = neighbor_searches.get(neighbor_index);
  Hermes::vector<MeshFunction*> slns;
  for (int i = 0; i < num; i++)
    slns.push_back(this->sln[i]);
  
  // Determine integration order.
  ExtData<Ord>* fake_ui = dp.init_ext_fns_ord(slns, neighbor_searches);
  
  // Order of additional external functions.
  // ExtData<Ord>* fake_ext = dp.init_ext_fns_ord(err_est_form->ext, nbs);

  // Order of geometric attributes (eg. for multiplication of a solution with coordinates, normals, etc.).
  Geom<Ord>* fake_e = new InterfaceGeom<Ord>(init_geom_ord(), nbs->neighb_el->marker, nbs->neighb_el->id, nbs->neighb_el->get_diameter());
  double fake_wt = 1.0;
  Ord o = err_est_form->ord(1, &fake_wt, fake_ui->fn, fake_ui->fn[err_est_form->i], fake_e, NULL);

  int order = rm->get_inv_ref_order();
  order += o.get_order();

  limit_order(order);

  // Clean up.
  if (fake_ui != NULL)
  {
    for (int i = 0; i < num; i++)
      delete fake_ui->fn[i];
    fake_ui->free_ord();
    delete fake_ui;
  }
  
  delete fake_e;
  
  //delete fake_ext;
  
  Quad2D* quad = this->sln[err_est_form->i]->get_quad_2d();
  int eo = quad->get_edge_points(surf_pos->surf_num, order);
  int np = quad->get_num_points(eo);
  double3* pt = quad->get_points(eo);
  
  // Init geometry and jacobian*weights (do not use the NeighborSearch caching mechanism).
  double3* tan = rm->get_tangent(surf_pos->surf_num, eo);
  double* jwt = new double[np];
  for(int i = 0; i < np; i++)
    jwt[i] = pt[i][2] * tan[i][2];
  
  Geom<double>* e = new InterfaceGeom<double>(init_geom_surf(rm, surf_pos, eo), 
                                              nbs->neighb_el->marker, 
                                              nbs->neighb_el->id, 
                                              nbs->neighb_el->get_diameter());
    
  // function values
  ExtData<scalar>* ui = dp.init_ext_fns(slns, neighbor_searches, order);
  //ExtData<scalar>* ext = dp.init_ext_fns(err_est_form->ext, nbs);

  scalar res = interface_scaling_const *
                err_est_form->value(np, jwt, ui->fn, ui->fn[err_est_form->i], e, NULL);

  if (ui != NULL) { ui->free(); delete ui; }
  //if (ext != NULL) { ext->free(); delete ext; }
  e->free(); delete e;
  delete [] jwt;

  return std::abs(0.5*res);   // Edges are parameterized from 0 to 1 while integration weights
                              // are defined in (-1, 1). Thus multiplying with 0.5 to correct
                              // the weights.
}

    void NewtonSolver<Scalar>::solve(Scalar* coeff_vec, double newton_tol, int newton_max_iter, bool residual_as_function)
    {
      _F_
      // Obtain the number of degrees of freedom.
      int ndof = this->dp->get_num_dofs();

      // Delete the old solution vector, if there is any.
      if(this->sln_vector != NULL)
        delete [] this->sln_vector;

      this->sln_vector = new Scalar[ndof];
      if(coeff_vec == NULL) 
        memset(this->sln_vector, 0, ndof*sizeof(Scalar));
      else
        for (int i = 0; i < ndof; i++)
          this->sln_vector[i] = coeff_vec[i];
          
      // The Newton's loop.
      double residual_norm;
      int it = 1;

      bool delete_timer = false;
      if (this->timer == NULL)
      {
        this->timer = new TimePeriod;
        delete_timer = true;
      }

      this->timer->tick();
      setup_time += this->timer->last();

      while (true)
      {
        // Assemble just the residual vector.
        this->dp->assemble(this->sln_vector, residual);

        this->timer->tick();
        assemble_time += this->timer->last();

        // Measure the residual norm.
        if (residual_as_function)
        {
          // Prepare solutions for measuring residual norm.
          Hermes::vector<Solution<Scalar>*> solutions;
          Hermes::vector<bool> dir_lift_false;
          for (unsigned int i = 0; i < static_cast<DiscreteProblem<Scalar>*>(this->dp)->get_spaces().size(); i++) {
            solutions.push_back(new Solution<Scalar>());
            dir_lift_false.push_back(false);
          }

          Solution<Scalar>::vector_to_solutions(residual, 
              static_cast<DiscreteProblem<Scalar>*>(this->dp)->get_spaces(), solutions, dir_lift_false);

          // Calculate the norm.
          residual_norm = Global<Scalar>::calc_norms(solutions);

          // Clean up.
          for (unsigned int i = 0; i < static_cast<DiscreteProblem<Scalar>*>(this->dp)->get_spaces().size(); i++)
            delete solutions[i];
        }
        else
        {
          // Calculate the l2-norm of residual vector, this is the traditional way.
          residual_norm = Global<Scalar>::get_l2_norm(residual);
        }

        // Info for the user.
        if(it == 1) {
          if(this->verbose_output)
            info("---- Newton initial residual norm: %g", residual_norm);
        }
        else
          if(this->verbose_output)
            info("---- Newton iter %d, residual norm: %g", it - 1, residual_norm);

        // If maximum allowed residual norm is exceeded, fail.
        if (residual_norm > max_allowed_residual_norm)
        {
          throw Exceptions::ValueException("residual norm", residual_norm, max_allowed_residual_norm);
        }

        // If residual norm is within tolerance, return 'true'.
        // This is the only correct way of ending.
        if (residual_norm < newton_tol && it > 1) 
        {
          this->timer->tick();
          solve_time += this->timer->last();

          if (delete_timer)
          {
            delete this->timer;
            this->timer = NULL;
          }

          return;
        }

        this->timer->tick();
        solve_time += this->timer->last();

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

// Filter for entropy which uses the constants defined above.
static void calc_entropy_estimate_func(int n, Hermes::vector<scalar*> scalars, scalar* result)
{
  for (int i = 0; i < n; i++)
    result[i] = std::log((calc_pressure(scalars.at(0)[i], scalars.at(1)[i], scalars.at(2)[i], scalars.at(3)[i]) / P_EXT)
    / pow((scalars.at(0)[i] / RHO_EXT), KAPPA));
};

      void L2ProjBasedSelector<Scalar>::precalc_ortho_shapes(const double3* gip_points, const int num_gip_points, const Trf* trfs, const int num_noni_trfs, const Hermes::vector<typename OptimumSelector<Scalar>::ShapeInx>& shapes, const int max_shape_inx, typename ProjBasedSelector<Scalar>::TrfShape& svals) 
      {
        //calculate values
        precalc_shapes(gip_points, num_gip_points, trfs, num_noni_trfs, shapes, max_shape_inx, svals);

        //calculate orthonormal basis
        const int num_shapes = (int)shapes.size();
        for(int i = 0; i < num_shapes; i++) 
        {
          const int inx_shape_i = shapes[i].inx;

          //orthogonalize
          for(int j = 0; j < i; j++) 
          {
            const int inx_shape_j = shapes[j].inx;

            //calculate product of non-transformed functions
            double product = 0.0;
            for(int k = 0; k < num_gip_points; k++) 
            {
              double sum = 0.0;
              sum += svals[H2D_TRF_IDENTITY][inx_shape_i][H2D_L2FE_VALUE][k] * svals[H2D_TRF_IDENTITY][inx_shape_j][H2D_L2FE_VALUE][k];
              product += gip_points[k][H2D_GIP2D_W] * sum;
            }

            //for all transformations
            int inx_trf = 0;
            bool done = false;
            while (!done && inx_trf < H2D_TRF_NUM) 
            {
              //for all integration points
              for(int k = 0; k < num_gip_points; k++) 
              {
                svals[inx_trf][inx_shape_i][H2D_L2FE_VALUE][k] -= product * svals[inx_trf][inx_shape_j][H2D_L2FE_VALUE][k];
              }

              //move to the next transformation
              if (inx_trf == H2D_TRF_IDENTITY)
                done = true;
              else 
              {
                inx_trf++;
                if (inx_trf >= num_noni_trfs) //if all transformations were processed, move to the identity transformation
                  inx_trf = H2D_TRF_IDENTITY;
              }
            }
            error_if(!done, "All transformation processed but identity transformation not found."); //identity transformation has to be the last transformation
          }

          //normalize
          //calculate norm
          double norm_squared = 0.0;
          for(int k = 0; k < num_gip_points; k++) 
          {
            double sum = 0.0;
            sum += sqr(svals[H2D_TRF_IDENTITY][inx_shape_i][H2D_L2FE_VALUE][k]);
            norm_squared += gip_points[k][H2D_GIP2D_W] * sum;
          }
          double norm = sqrt(norm_squared);
          assert_msg(finite(1/norm), "Norm (%g) is almost zero.", norm);

          //for all transformations: normalize
          int inx_trf = 0;
          bool done = false;
          while (!done && inx_trf < H2D_TRF_NUM) 
          {
            //for all integration points
            for(int k = 0; k < num_gip_points; k++) 
            {
              svals[inx_trf][inx_shape_i][H2D_L2FE_VALUE][k] /= norm;
            }

            //move to the next transformation
            if (inx_trf == H2D_TRF_IDENTITY)
              done = true;
            else 
            {
              inx_trf++;
              if (inx_trf >= num_noni_trfs) //if all transformations were processed, move to the identity transformation
                inx_trf = H2D_TRF_IDENTITY;
            }
          }
          error_if(!done, "All transformation processed but identity transformation not found."); //identity transformation has to be the last transformation
        }
      }

QList<SolutionArray *> SolutionAgros::solveSolutioArray(Hermes::vector<EssentialBCs> bcs)
{
    QTime time;

    // solution agros array
    QList<SolutionArray *> solutionArrayList;

    // load the mesh file
    mesh = readMeshFromFile(tempProblemFileName() + ".mesh");
    refineMesh(mesh, true, true);

    // create an H1 space
    Hermes::vector<Space *> space;
    // create hermes solution array
    Hermes::vector<Solution *> solution;
    // create reference solution
    Hermes::vector<Solution *> solutionReference;

    // projection norms
    Hermes::vector<ProjNormType> projNormType;

    // prepare selector
    Hermes::vector<RefinementSelectors::Selector *> selector;

    // error marker
    bool isError = false;

    RefinementSelectors::Selector *select = NULL;
    switch (adaptivityType)
    {
    case AdaptivityType_H:
        select = new RefinementSelectors::HOnlySelector();
        break;
    case AdaptivityType_P:
        select = new RefinementSelectors::H1ProjBasedSelector(RefinementSelectors::H2D_P_ANISO,
                                                              Util::config()->convExp,
                                                              H2DRS_DEFAULT_ORDER);
        break;
    case AdaptivityType_HP:
        select = new RefinementSelectors::H1ProjBasedSelector(RefinementSelectors::H2D_HP_ANISO,
                                                              Util::config()->convExp,
                                                              H2DRS_DEFAULT_ORDER);
        break;
    }

    for (int i = 0; i < numberOfSolution; i++)
    {
        space.push_back(new H1Space(mesh, &bcs[i], polynomialOrder));

        // set order by element
        for (int j = 0; j < Util::scene()->labels.count(); j++)
            if (Util::scene()->labels[j]->material != Util::scene()->materials[0])
                space.at(i)->set_uniform_order(Util::scene()->labels[j]->polynomialOrder > 0 ? Util::scene()->labels[j]->polynomialOrder : polynomialOrder,
                                               QString::number(j).toStdString());

        // solution agros array
        solution.push_back(new Solution());

        if (adaptivityType != AdaptivityType_None)
        {
            // add norm
            projNormType.push_back(Util::config()->projNormType);
            // add refinement selector
            selector.push_back(select);
            // reference solution
            solutionReference.push_back(new Solution());
        }
    }

    // check for DOFs
    if (Space::get_num_dofs(space) == 0)
    {
        m_progressItemSolve->emitMessage(QObject::tr("DOF is zero"), true);
    }
    else
    {
        for (int i = 0; i < numberOfSolution; i++)
        {
            // transient
            if (analysisType == AnalysisType_Transient)
            {
                // constant initial solution
                solution.at(i)->set_const(mesh, initialCondition);
                solutionArrayList.append(solutionArray(solution.at(i)));
            }

            // nonlinear
            if ((linearityType != LinearityType_Linear) && (analysisType != AnalysisType_Transient))
            {
                solution.at(i)->set_const(mesh, 0.0);
            }
        }

        actualTime = 0.0;

        // update time function
        Util::scene()->problemInfo()->hermes()->updateTimeFunctions(actualTime);

        m_wf->set_current_time(actualTime);
        m_wf->solution = solution;
//.........这里部分代码省略.........

bool rk_time_step(double current_time, double time_step, ButcherTable* const bt,
                  Solution* sln_time_prev, Solution* sln_time_new, Solution* error_fn, 
                  DiscreteProblem* dp, MatrixSolverType matrix_solver,
                  bool verbose, bool is_linear, double newton_tol, int newton_max_iter,
                  double newton_damping_coeff, double newton_max_allowed_residual_norm)
{
  // Check for not implemented features.
  if (matrix_solver != SOLVER_UMFPACK)
    error("Sorry, rk_time_step() still only works with UMFpack.");
  if (dp->get_weak_formulation()->get_neq() > 1)
    error("Sorry, rk_time_step() does not work with systems yet.");

  // Get number of stages from the Butcher's table.
  int num_stages = bt->get_size();

  // Check whether the user provided a nonzero B2-row if he wants temporal error estimation.
  if(error_fn != NULL) if (bt->is_embedded() == false) {
    error("rk_time_step(): R-K method must be embedded if temporal error estimate is requested.");
  }

  // Matrix for the time derivative part of the equation (left-hand side).
  UMFPackMatrix* matrix_left = new UMFPackMatrix();

  // Matrix and vector for the rest (right-hand side).
  UMFPackMatrix* matrix_right = new UMFPackMatrix();
  UMFPackVector* vector_right = new UMFPackVector();

  // Create matrix solver.
  Solver* solver = create_linear_solver(matrix_solver, matrix_right, vector_right);

  // Get space, mesh, and ndof for the stage solutions in the R-K method (K_i vectors).
  Space* K_space = dp->get_space(0);
  Mesh* K_mesh = K_space->get_mesh();
  int ndof = K_space->get_num_dofs();

  // Create spaces for stage solutions K_i. This is necessary
  // to define a num_stages x num_stages block weak formulation.
  Hermes::vector<Space*> stage_spaces;
  stage_spaces.push_back(K_space);
  for (int i = 1; i < num_stages; i++) {
    stage_spaces.push_back(K_space->dup(K_mesh));
  }
  Space::assign_dofs(stage_spaces);

  // Create a multistage weak formulation.
  WeakForm stage_wf_left;                   // For the matrix M (size ndof times ndof).
  WeakForm stage_wf_right(num_stages);      // For the rest of equation (written on the right),
                                            // size num_stages*ndof times num_stages*ndof.

  Solution** stage_time_sol = new Solution*[num_stages];
                                            // This array will be filled by artificially created
                                            // solutions to represent stage times.
  create_stage_wf(current_time, time_step, bt, dp, &stage_wf_left, &stage_wf_right, stage_time_sol); 

  // Initialize discrete problems for the assembling of the
  // matrix M and the stage Jacobian matrix and residual.
  DiscreteProblem stage_dp_left(&stage_wf_left, K_space);
  DiscreteProblem stage_dp_right(&stage_wf_right, stage_spaces);

  // Vector K_vector of length num_stages * ndof. will represent
  // the 'K_i' vectors in the usual R-K notation.
  scalar* K_vector = new scalar[num_stages*ndof];
  memset(K_vector, 0, num_stages * ndof * sizeof(scalar));

  // Vector u_ext_vec will represent h \sum_{j=1}^s a_{ij} K_i.
  scalar* u_ext_vec = new scalar[num_stages*ndof];

  // Vector for the left part of the residual.
  scalar* vector_left = new scalar[num_stages*ndof];

  // Prepare residuals of stage solutions.
  Hermes::vector<Solution*> residuals;
  Hermes::vector<bool> add_dir_lift;
  for (int i = 0; i < num_stages; i++) {
    residuals.push_back(new Solution(K_mesh));
    add_dir_lift.push_back(false);
  }

  // Assemble the block-diagonal mass matrix M of size ndof times ndof.
  // The corresponding part of the global residual vector is obtained 
  // just by multiplication.
  stage_dp_left.assemble(matrix_left);

  // The Newton's loop.
  double residual_norm;
  int it = 1;
  while (true)
  {
    // Prepare vector h\sum_{j=1}^s a_{ij} K_j.
    for (int i = 0; i < num_stages; i++) {                // block row
      for (int idx = 0; idx < ndof; idx++) {
        scalar increment = 0;
        for (int j = 0; j < num_stages; j++) {
          increment += bt->get_A(i, j) * K_vector[j*ndof + idx];
        }
        u_ext_vec[i*ndof + idx] = time_step * increment;
      }
    }

    multiply_as_diagonal_block_matrix(matrix_left, num_stages, K_vector, vector_left);
//.........这里部分代码省略.........

bool rk_time_step(double current_time, double time_step, ButcherTable* const bt,
                  scalar* coeff_vec, scalar* err_vec, DiscreteProblem* dp, MatrixSolverType matrix_solver,
                  bool verbose, bool is_linear, double newton_tol, int newton_max_iter,
                  double newton_damping_coeff, double newton_max_allowed_residual_norm)
{
  // Check for not implemented features.
  if (matrix_solver != SOLVER_UMFPACK)
    error("Sorry, rk_time_step() still only works with UMFpack.");
  if (dp->get_weak_formulation()->get_neq() > 1)
    error("Sorry, rk_time_step() does not work with systems yet.");

  // Get number of stages from the Butcher's table.
  int num_stages = bt->get_size();

  // Check whether the user provided a second B-row if he wants 
  // err_vec.
  if(err_vec != NULL) {
    double b2_coeff_sum = 0;
    for (int i=0; i < num_stages; i++) b2_coeff_sum += fabs(bt->get_B2(i)); 
    if (b2_coeff_sum < 1e-10) 
      error("err_vec != NULL but the B2 row in the Butcher's table is zero in rk_time_step().");
  }

  // Matrix for the time derivative part of the equation (left-hand side).
  UMFPackMatrix* matrix_left = new UMFPackMatrix();

  // Matrix and vector for the rest (right-hand side).
  UMFPackMatrix* matrix_right = new UMFPackMatrix();
  UMFPackVector* vector_right = new UMFPackVector();

  // Create matrix solver.
  Solver* solver = create_linear_solver(matrix_solver, matrix_right, vector_right);

  // Get original space, mesh, and ndof.
  dp->get_space(0);
  Mesh* mesh = dp->get_space(0)->get_mesh();
  int ndof = dp->get_space(0)->get_num_dofs();

  // Create spaces for stage solutions. This is necessary
  // to define a num_stages x num_stages block weak formulation.
  Hermes::vector<Space*> stage_spaces;
  stage_spaces.push_back(dp->get_space(0));
  for (int i = 1; i < num_stages; i++) {
    stage_spaces.push_back(dp->get_space(0)->dup(mesh));
  }
  Space::assign_dofs(stage_spaces);

  // Create a multistage weak formulation.
  WeakForm stage_wf_left;                   // For the matrix M (size ndof times ndof).
  WeakForm stage_wf_right(num_stages);      // For the rest of equation (written on the right),
                                            // size num_stages*ndof times num_stages*ndof.
  create_stage_wf(current_time, time_step, bt, dp, &stage_wf_left, &stage_wf_right); 

  // Initialize discrete problems for the assembling of the
  // matrix M and the stage Jacobian matrix and residual.
  DiscreteProblem stage_dp_left(&stage_wf_left, dp->get_space(0));
  DiscreteProblem stage_dp_right(&stage_wf_right, stage_spaces);

  // Vector K_vector of length num_stages * ndof. will represent
  // the 'k_i' vectors in the usual R-K notation.
  scalar* K_vector = new scalar[num_stages*ndof];
  memset(K_vector, 0, num_stages * ndof * sizeof(scalar));

  // Vector u_prev_vec will represent y_n + h \sum_{j=1}^s a_{ij}k_i
  // in the usual R-K notation.
  scalar* u_prev_vec = new scalar[num_stages*ndof];

  // Vector for the left part of the residual.
  scalar* vector_left = new scalar[num_stages*ndof];

  // Prepare residuals of stage solutions.
  Hermes::vector<Solution*> residuals;
  Hermes::vector<bool> add_dir_lift;
  for (int i = 0; i < num_stages; i++) {
    residuals.push_back(new Solution(mesh));
    add_dir_lift.push_back(false);
  }

  // Assemble the block-diagonal mass matrix M of size ndof times ndof.
  // The corresponding part of the global residual vector is obtained 
  // just by multiplication.
  stage_dp_left.assemble(matrix_left);

  // The Newton's loop.
  double residual_norm;
  int it = 1;
  while (true)
  {
    // Prepare vector Y_n + h\sum_{j=1}^s a_{ij} K_j.
    for (int i = 0; i < num_stages; i++) {                // block row
      for (int idx = 0; idx < ndof; idx++) {
        scalar increment = 0;
        for (int j = 0; j < num_stages; j++) {
          increment += bt->get_A(i, j) * K_vector[j*ndof + idx];
        }
        u_prev_vec[i*ndof + idx] = coeff_vec[idx] + time_step * increment;
      }
    }

    multiply_as_diagonal_block_matrix(matrix_left, num_stages, 
//.........这里部分代码省略.........

int main() 
{
  // Time measurement.
  TimePeriod cpu_time;
  cpu_time.tick();

  // Create space, set Dirichlet BC, enumerate basis functions.
  Space* space = new Space(A, B, NELEM, DIR_BC_LEFT, DIR_BC_RIGHT, 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;
  bool success = false;
  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(!(success = 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++;
  }
  info("Total running time: %g s", cpu_time.accumulated());

  // Test variable.
  info("ndof = %d.", Space::get_num_dofs(space));

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

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

    bool NewtonSolver<Scalar>::sol