C++ ErrorVector类代码示例

void MeshRefinement::flag_elements_by_mean_stddev (const ErrorVector & error_per_cell,
                                                   const Real refine_frac,
                                                   const Real coarsen_frac,
                                                   const unsigned int max_l)
{
  // The function arguments are currently just there for
  // backwards_compatibility
  if (!_use_member_parameters)
    {
      // If the user used non-default parameters, lets warn
      // that they're deprecated
      if (refine_frac != 0.3 ||
          coarsen_frac != 0.0 ||
          max_l != libMesh::invalid_uint)
        libmesh_deprecated();

      _refine_fraction = refine_frac;
      _coarsen_fraction = coarsen_frac;
      _max_h_level = max_l;
    }

  // Get the mean value from the error vector
  const Real mean = error_per_cell.mean();

  // Get the standard deviation.  This equals the
  // square-root of the variance
  const Real stddev = std::sqrt (error_per_cell.variance());

  // Check for valid fractions
  libmesh_assert_greater_equal (_refine_fraction, 0);
  libmesh_assert_less_equal (_refine_fraction, 1);
  libmesh_assert_greater_equal (_coarsen_fraction, 0);
  libmesh_assert_less_equal (_coarsen_fraction, 1);

  // The refine and coarsen cutoff
  const Real refine_cutoff  =  mean + _refine_fraction  * stddev;
  const Real coarsen_cutoff =  std::max(mean - _coarsen_fraction * stddev, 0.);

  // Loop over the elements and flag them for coarsening or
  // refinement based on the element error
  for (auto & elem : _mesh.active_element_ptr_range())
    {
      const dof_id_type id  = elem->id();

      libmesh_assert_less (id, error_per_cell.size());

      const ErrorVectorReal elem_error = error_per_cell[id];

      // Possibly flag the element for coarsening ...
      if (elem_error <= coarsen_cutoff)
        elem->set_refinement_flag(Elem::COARSEN);

      // ... or refinement
      if ((elem_error >= refine_cutoff) && (elem->level() < _max_h_level))
        elem->set_refinement_flag(Elem::REFINE);
    }
}

realv Layer::errorWeighting(ErrorVector _deltas, Mat _weights) {
	if (_deltas.getLength() != ((uint) _weights.rows)) {
		throw length_error("Layer : Uncorrect length between deltas and weights");
	}
	realv sum = 0;
	for (uint i = 0; i < _deltas.getLength(); i++) {
		sum += _deltas[i] * _weights.at<realv>(i, 0);
	}
	return sum;
}

void ErrorEstimator::estimate_errors(const EquationSystems & equation_systems,
                                     ErrorVector & error_per_cell,
                                     const std::map<const System *, SystemNorm> & error_norms,
                                     const std::map<const System *, const NumericVector<Number> *> * solution_vectors,
                                     bool estimate_parent_error)
{
  SystemNorm old_error_norm = this->error_norm;

  // Sum the error values from each system
  for (unsigned int s = 0; s != equation_systems.n_systems(); ++s)
    {
      ErrorVector system_error_per_cell;
      const System & sys = equation_systems.get_system(s);
      if (error_norms.find(&sys) == error_norms.end())
        this->error_norm = old_error_norm;
      else
        this->error_norm = error_norms.find(&sys)->second;

      const NumericVector<Number> * solution_vector = nullptr;
      if (solution_vectors &&
          solution_vectors->find(&sys) != solution_vectors->end())
        solution_vector = solution_vectors->find(&sys)->second;

      this->estimate_error(sys, system_error_per_cell,
                           solution_vector, estimate_parent_error);

      if (s)
        {
          libmesh_assert_equal_to (error_per_cell.size(), system_error_per_cell.size());
          for (std::size_t i=0; i != error_per_cell.size(); ++i)
            error_per_cell[i] += system_error_per_cell[i];
        }
      else
        error_per_cell = system_error_per_cell;
    }

  // Restore our old state before returning
  this->error_norm = old_error_norm;
}

//.........这里部分代码省略.........
      // Solve the system "Biharmonic", just like example 2.
      system.solve();

      std::cout << "Linear solver converged at step: "
                << system.n_linear_iterations()
                << ", final residual: "
                << system.final_linear_residual()
                << std::endl;

      // Compute the error.
      exact_sol.compute_error("Biharmonic", "u");
      // Compute the norm.
      zero_sol.compute_error("Biharmonic", "u");

      // Print out the error values
      std::cout << "L2-Norm is: "
                << zero_sol.l2_error("Biharmonic", "u")
                << std::endl;
      std::cout << "H1-Norm is: "
                << zero_sol.h1_error("Biharmonic", "u")
                << std::endl;
      std::cout << "H2-Norm is: "
                << zero_sol.h2_error("Biharmonic", "u")
                << std::endl
                << std::endl;
      std::cout << "L2-Error is: "
                << exact_sol.l2_error("Biharmonic", "u")
                << std::endl;
      std::cout << "H1-Error is: "
                << exact_sol.h1_error("Biharmonic", "u")
                << std::endl;
      std::cout << "H2-Error is: "
                << exact_sol.h2_error("Biharmonic", "u")
                << std::endl
                << std::endl;

      // Print to output file
      out << equation_systems.n_active_dofs() << " "
          << exact_sol.l2_error("Biharmonic", "u") << " "
          << exact_sol.h1_error("Biharmonic", "u") << " "
          << exact_sol.h2_error("Biharmonic", "u") << std::endl;

      // Possibly refine the mesh
      if (r_step+1 != max_r_steps)
        {
          std::cout << "  Refining the mesh..." << std::endl;

          if (uniform_refine == 0)
            {
              ErrorVector error;
              LaplacianErrorEstimator error_estimator;

              error_estimator.estimate_error(system, error);
              mesh_refinement.flag_elements_by_elem_fraction (error);

              std::cout << "Mean Error: " << error.mean() <<
                              std::endl;
              std::cout << "Error Variance: " << error.variance() <<
                              std::endl;

              mesh_refinement.refine_and_coarsen_elements();
            }
          else
            {
              mesh_refinement.uniformly_refine(1);
            }
              
          // This call reinitializes the \p EquationSystems object for
          // the newly refined mesh.  One of the steps in the
          // reinitialization is projecting the \p solution,
          // \p old_solution, etc... vectors from the old mesh to
          // the current one.
          equation_systems.reinit ();
        }
    }            
  
#ifdef LIBMESH_HAVE_EXODUS_API
  // Write out the solution
  // After solving the system write the solution
  // to a ExodusII-formatted plot file.
  ExodusII_IO (mesh).write_equation_systems (exd_file,
                                       equation_systems);
#endif // #ifdef LIBMESH_HAVE_EXODUS_API

  // Close up the output file.
  out << "];\n"
      << "hold on\n"
      << "loglog(e(:,1), e(:,2), 'bo-');\n"
      << "loglog(e(:,1), e(:,3), 'ro-');\n"
      << "loglog(e(:,1), e(:,4), 'go-');\n"
      << "xlabel('log(dofs)');\n"
      << "ylabel('log(error)');\n"
      << "title('C1 " << approx_type << " elements');\n"
      << "legend('L2-error', 'H1-error', 'H2-error');\n";
  
  // All done.  
  return 0;
#endif // #ifndef LIBMESH_ENABLE_SECOND_DERIVATIVES
#endif // #ifndef LIBMESH_ENABLE_AMR
}

void AdjointRefinementEstimator::estimate_error (const System & _system,
                                                 ErrorVector & error_per_cell,
                                                 const NumericVector<Number> * solution_vector,
                                                 bool /*estimate_parent_error*/)
{
  // We have to break the rules here, because we can't refine a const System
  System & system = const_cast<System &>(_system);

  // An EquationSystems reference will be convenient.
  EquationSystems & es = system.get_equation_systems();

  // The current mesh
  MeshBase & mesh = es.get_mesh();

  // Get coarse grid adjoint solutions.  This should be a relatively
  // quick (especially with preconditioner reuse) way to get a good
  // initial guess for the fine grid adjoint solutions.  More
  // importantly, subtracting off a coarse adjoint approximation gives
  // us better local error indication, and subtracting off *some* lift
  // function is necessary for correctness if we have heterogeneous
  // adjoint Dirichlet conditions.

  // Solve the adjoint problem(s) on the coarse FE space
  // Only if the user didn't already solve it for us
  if (!system.is_adjoint_already_solved())
    system.adjoint_solve(_qoi_set);

  // Loop over all the adjoint problems and, if any have heterogenous
  // Dirichlet conditions, get the corresponding coarse lift
  // function(s)
  for (unsigned int j=0; j != system.qoi.size(); j++)
    {
      // Skip this QoI if it is not in the QoI Set or if there are no
      // heterogeneous Dirichlet boundaries for it
      if (_qoi_set.has_index(j) &&
          system.get_dof_map().has_adjoint_dirichlet_boundaries(j))
        {
          std::ostringstream liftfunc_name;
          liftfunc_name << "adjoint_lift_function" << j;
          NumericVector<Number> & liftvec =
            system.add_vector(liftfunc_name.str());

          system.get_dof_map().enforce_constraints_exactly
            (system, &liftvec, true);
        }
    }

  // We'll want to back up all coarse grid vectors
  std::map<std::string, NumericVector<Number> *> coarse_vectors;
  for (System::vectors_iterator vec = system.vectors_begin(); vec !=
         system.vectors_end(); ++vec)
    {
      // The (string) name of this vector
      const std::string & var_name = vec->first;

      coarse_vectors[var_name] = vec->second->clone().release();
    }
  // Back up the coarse solution and coarse local solution
  NumericVector<Number> * coarse_solution =
    system.solution->clone().release();
  NumericVector<Number> * coarse_local_solution =
    system.current_local_solution->clone().release();

  // And we'll need to temporarily change solution projection settings
  bool old_projection_setting;
  old_projection_setting = system.project_solution_on_reinit();

  // Make sure the solution is projected when we refine the mesh
  system.project_solution_on_reinit() = true;

  // And it'll be best to avoid any repartitioning
  UniquePtr<Partitioner> old_partitioner(mesh.partitioner().release());

  // And we can't allow any renumbering
  const bool old_renumbering_setting = mesh.allow_renumbering();
  mesh.allow_renumbering(false);

  // Use a non-standard solution vector if necessary
  if (solution_vector && solution_vector != system.solution.get())
    {
      NumericVector<Number> * newsol =
        const_cast<NumericVector<Number> *> (solution_vector);
      newsol->swap(*system.solution);
      system.update();
    }

  // Resize the error_per_cell vector to be
  // the number of elements, initialized to 0.
  error_per_cell.clear();
  error_per_cell.resize (mesh.max_elem_id(), 0.);

#ifndef NDEBUG
  // n_coarse_elem is only used in an assertion later so
  // avoid declaring it unless asserts are active.
  const dof_id_type n_coarse_elem = mesh.n_elem();
#endif

  // Uniformly refine the mesh
  MeshRefinement mesh_refinement(mesh);

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

void ExactErrorEstimator::estimate_error (const System & system,
                                          ErrorVector & error_per_cell,
                                          const NumericVector<Number> * solution_vector,
                                          bool estimate_parent_error)
{
  // Ignore the fact that this variable is unused when !LIBMESH_ENABLE_AMR
  libmesh_ignore(estimate_parent_error);

  // The current mesh
  const MeshBase & mesh = system.get_mesh();

  // The dimensionality of the mesh
  const unsigned int dim = mesh.mesh_dimension();

  // The number of variables in the system
  const unsigned int n_vars = system.n_vars();

  // The DofMap for this system
  const DofMap & dof_map = system.get_dof_map();

  // Resize the error_per_cell vector to be
  // the number of elements, initialize it to 0.
  error_per_cell.resize (mesh.max_elem_id());
  std::fill (error_per_cell.begin(), error_per_cell.end(), 0.);

  // Prepare current_local_solution to localize a non-standard
  // solution vector if necessary
  if (solution_vector && solution_vector != system.solution.get())
    {
      NumericVector<Number> * newsol =
        const_cast<NumericVector<Number> *>(solution_vector);
      System & sys = const_cast<System &>(system);
      newsol->swap(*sys.solution);
      sys.update();
    }

  // Loop over all the variables in the system
  for (unsigned int var=0; var<n_vars; var++)
    {
      // Possibly skip this variable
      if (error_norm.weight(var) == 0.0) continue;

      // The (string) name of this variable
      const std::string & var_name = system.variable_name(var);

      // The type of finite element to use for this variable
      const FEType & fe_type = dof_map.variable_type (var);

      UniquePtr<FEBase> fe (FEBase::build (dim, fe_type));

      // Build an appropriate Gaussian quadrature rule
      UniquePtr<QBase> qrule =
        fe_type.default_quadrature_rule (dim,
                                         _extra_order);

      fe->attach_quadrature_rule (qrule.get());

      // Prepare a global solution and a MeshFunction of the fine system if we need one
      UniquePtr<MeshFunction> fine_values;
      UniquePtr<NumericVector<Number> > fine_soln = NumericVector<Number>::build(system.comm());
      if (_equation_systems_fine)
        {
          const System & fine_system = _equation_systems_fine->get_system(system.name());

          std::vector<Number> global_soln;
          // FIXME - we're assuming that the fine system solution gets
          // used even when a different vector is used for the coarse
          // system
          fine_system.update_global_solution(global_soln);
          fine_soln->init
            (cast_int<numeric_index_type>(global_soln.size()), true,
             SERIAL);
          (*fine_soln) = global_soln;

          fine_values = UniquePtr<MeshFunction>
            (new MeshFunction(*_equation_systems_fine,
                              *fine_soln,
                              fine_system.get_dof_map(),
                              fine_system.variable_number(var_name)));
          fine_values->init();
        } else {
        // Initialize functors if we're using them
        for (unsigned int i=0; i != _exact_values.size(); ++i)
          if (_exact_values[i])
            _exact_values[i]->init();

        for (unsigned int i=0; i != _exact_derivs.size(); ++i)
          if (_exact_derivs[i])
            _exact_derivs[i]->init();

        for (unsigned int i=0; i != _exact_hessians.size(); ++i)
          if (_exact_hessians[i])
            _exact_hessians[i]->init();
      }

      // Request the data we'll need to compute with
      fe->get_JxW();
      fe->get_phi();
      fe->get_dphi();
#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
//.........这里部分代码省略.........

//.........这里部分代码省略.........
  // Mesh Refinement object - to test effect of constant refined mesh (not refined at every timestep)
  MeshRefinement mesh_refinement(mesh);
  mesh_refinement.refine_fraction() = refine_percentage;
  mesh_refinement.coarsen_fraction() = coarsen_percentage;
  mesh_refinement.max_h_level() = max_r_level;

  // Print information about the system to the screen.
  equation_systems.print_info();
  
  ExodusII_IO exodusIO = ExodusII_IO(mesh); //for writing multiple timesteps to one file
  
  for (unsigned int r_step=0; r_step<max_r_steps; r_step++)
    {
      std::cout << "\nBeginning Solve " << r_step+1 << std::endl;
      
      for (unsigned int t_step=0; t_step != n_timesteps; ++t_step)
        {
      
          std::cout << "\n\nSolving time step " << t_step << ", time = "
                    << system.time << std::endl;
          
          system.solve();
          system.postprocess();
          
          // Advance to the next timestep in a transient problem
          system.time_solver->advance_timestep();
   
      } //end stepping through time loop
      
      std::cout << "\n  Refining the mesh..." << std::endl;
      
      // The \p ErrorVector is a particular \p StatisticsVector
      // for computing error information on a finite element mesh.
      ErrorVector error;
      
      if (indicator_type == "patch")
        {
          // The patch recovery estimator should give a
          // good estimate of the solution interpolation
          // error.
          PatchRecoveryErrorEstimator error_estimator;
          error_estimator.set_patch_reuse(false); //anisotropy trips up reuse
          error_estimator.estimate_error (system, error);
        }
      else if (indicator_type == "kelly")
        {

          // The Kelly error estimator is based on
          // an error bound for the Poisson problem
          // on linear elements, but is useful for
          // driving adaptive refinement in many problems
          KellyErrorEstimator error_estimator;

          error_estimator.estimate_error (system, error);
        }
      
      // Write out the error distribution
      if(write_error){
        std::ostringstream ss;
        ss << r_step;
#ifdef LIBMESH_HAVE_EXODUS_API
        std::string error_output = "error_"+ss.str()+".e";
#else
        std::string error_output = "error_"+ss.str()+".gmv";
#endif
        error.plot_error( error_output, mesh );

//.........这里部分代码省略.........
    solver->relative_step_tolerance =
        infile("relative_step_tolerance", 1.e-3);
    solver->relative_residual_tolerance =
        infile("relative_residual_tolerance", 0.0);
    solver->absolute_residual_tolerance =
        infile("absolute_residual_tolerance", 0.0);

    // And the linear solver options
    solver->max_linear_iterations =
        infile("max_linear_iterations", 50000);
    solver->initial_linear_tolerance =
        infile("initial_linear_tolerance", 1.e-3);

    // Print information about the system to the screen.
    equation_systems.print_info();

    // Now we begin the timestep loop to compute the time-accurate
    // solution of the equations...not that this is transient, but eh, why not...
    for (unsigned int t_step=0; t_step != n_timesteps; ++t_step)
    {
        // A pretty update message
        std::cout << "\n\nSolving time step " << t_step << ", time = "
                  << system.time << std::endl;

        // Adaptively solve the timestep
        unsigned int a_step = 0;
        for (; a_step != max_adaptivesteps; ++a_step)
        {   // VESTIGIAL for now ('vestigial' eh ? ;) )

            std::cout << "\n\n I should be skipped what are you doing here lalalalalalala *!**!*!*!*!*!* \n\n";

            system.solve();
            system.postprocess();
            ErrorVector error;
            AutoPtr<ErrorEstimator> error_estimator;

            // To solve to a tolerance in this problem we
            // need a better estimator than Kelly
            if (global_tolerance != 0.)
            {
                // We can't adapt to both a tolerance and a mesh
                // size at once
                libmesh_assert_equal_to (nelem_target, 0);

                UniformRefinementEstimator *u =
                    new UniformRefinementEstimator;

                // The lid-driven cavity problem isn't in H1, so
                // lets estimate L2 error
                u->error_norm = L2;

                error_estimator.reset(u);
            }
            else
            {
                // If we aren't adapting to a tolerance we need a
                // target mesh size
                libmesh_assert_greater (nelem_target, 0);

                // Kelly is a lousy estimator to use for a problem
                // not in H1 - if we were doing more than a few
                // timesteps we'd need to turn off or limit the
                // maximum level of our adaptivity eventually
                error_estimator.reset(new KellyErrorEstimator);
            }

void AdjointResidualErrorEstimator::estimate_error (const System & _system,
                                                    ErrorVector & error_per_cell,
                                                    const NumericVector<Number> * solution_vector,
                                                    bool estimate_parent_error)
{
  LOG_SCOPE("estimate_error()", "AdjointResidualErrorEstimator");

  // The current mesh
  const MeshBase & mesh = _system.get_mesh();

  // Resize the error_per_cell vector to be
  // the number of elements, initialize it to 0.
  error_per_cell.resize (mesh.max_elem_id());
  std::fill (error_per_cell.begin(), error_per_cell.end(), 0.);

  // Get the number of variables in the system
  unsigned int n_vars = _system.n_vars();

  // We need to make a map of the pointer to the solution vector
  std::map<const System *, const NumericVector<Number> *>solutionvecs;
  solutionvecs[&_system] = _system.solution.get();

  // Solve the dual problem if we have to
  if (!_system.is_adjoint_already_solved())
    {
      // FIXME - we'll need to change a lot of APIs to make this trick
      // work with a const System...
      System &  system = const_cast<System &>(_system);
      system.adjoint_solve(_qoi_set);
    }

  // Flag to check whether we have not been asked to weight the variable error contributions in any specific manner
  bool error_norm_is_identity = error_norm.is_identity();

  // Create an ErrorMap/ErrorVector to store the primal, dual and total_dual variable errors
  ErrorMap primal_errors_per_cell;
  ErrorMap dual_errors_per_cell;
  ErrorMap total_dual_errors_per_cell;
  // Allocate ErrorVectors to this map if we're going to use it
  if (!error_norm_is_identity)
    for(unsigned int v = 0; v < n_vars; v++)
      {
        primal_errors_per_cell[std::make_pair(&_system, v)] = new ErrorVector;
        dual_errors_per_cell[std::make_pair(&_system, v)] = new ErrorVector;
        total_dual_errors_per_cell[std::make_pair(&_system, v)] = new ErrorVector;
      }
  ErrorVector primal_error_per_cell;
  ErrorVector dual_error_per_cell;
  ErrorVector total_dual_error_per_cell;

  // Have we been asked to weight the variable error contributions in any specific manner
  if(!error_norm_is_identity) // If we do
    {
      // Estimate the primal problem error for each variable
      _primal_error_estimator->estimate_errors
        (_system.get_equation_systems(), primal_errors_per_cell, &solutionvecs, estimate_parent_error);
    }
  else // If not
    {
      // Just get the combined error estimate
      _primal_error_estimator->estimate_error
        (_system, primal_error_per_cell, solution_vector, estimate_parent_error);
    }

  // Sum and weight the dual error estimate based on our QoISet
  for (unsigned int i = 0; i != _system.qoi.size(); ++i)
    {
      if (_qoi_set.has_index(i))
        {
          // Get the weight for the current QoI
          Real error_weight = _qoi_set.weight(i);

          // We need to make a map of the pointer to the adjoint solution vector
          std::map<const System *, const NumericVector<Number> *>adjointsolutionvecs;
          adjointsolutionvecs[&_system] = &_system.get_adjoint_solution(i);

          // Have we been asked to weight the variable error contributions in any specific manner
          if(!error_norm_is_identity) // If we have
            {
              _dual_error_estimator->estimate_errors
                (_system.get_equation_systems(), dual_errors_per_cell, &adjointsolutionvecs,
                 estimate_parent_error);
            }
          else // If not
            {
              // Just get the combined error estimate
              _dual_error_estimator->estimate_error
                (_system, dual_error_per_cell, &(_system.get_adjoint_solution(i)), estimate_parent_error);
            }

          std::size_t error_size;

          // Get the size of the first ErrorMap vector; this will give us the number of elements
          if(!error_norm_is_identity) // If in non default weights case
            {
              error_size = dual_errors_per_cell[std::make_pair(&_system, 0)]->size();
            }
          else // If in the standard default weights case
            {
              error_size = dual_error_per_cell.size();
//.........这里部分代码省略.........

//-----------------------------------------------------------------
// Mesh refinement methods
void MeshRefinement::flag_elements_by_error_fraction (const ErrorVector& error_per_cell,
						      const Real refine_frac,
						      const Real coarsen_frac,
						      const unsigned int max_l)
{
  parallel_only();

  // The function arguments are currently just there for
  // backwards_compatibility
  if (!_use_member_parameters)
  {
    // If the user used non-default parameters, lets warn
    // that they're deprecated
    if (refine_frac != 0.3 ||
	coarsen_frac != 0.0 ||
	max_l != libMesh::invalid_uint)
      libmesh_deprecated();

    _refine_fraction = refine_frac;
    _coarsen_fraction = coarsen_frac;
    _max_h_level = max_l;
  }

  // Check for valid fractions..
  // The fraction values must be in [0,1]
  libmesh_assert_greater_equal (_refine_fraction, 0);
  libmesh_assert_less_equal (_refine_fraction, 1);
  libmesh_assert_greater_equal (_coarsen_fraction, 0);
  libmesh_assert_less_equal (_coarsen_fraction, 1);

  // Clean up the refinement flags.  These could be left
  // over from previous refinement steps.
  this->clean_refinement_flags();

  // We're getting the minimum and maximum error values
  // for the ACTIVE elements
  Real error_min = 1.e30;
  Real error_max = 0.;

  // And, if necessary, for their parents
  Real parent_error_min = 1.e30;
  Real parent_error_max = 0.;

  // Prepare another error vector if we need to sum parent errors
  ErrorVector error_per_parent;
  if (_coarsen_by_parents)
  {
    create_parent_error_vector(error_per_cell,
			       error_per_parent,
			       parent_error_min,
			       parent_error_max);
  }

  // We need to loop over all active elements to find the minimum
  MeshBase::element_iterator       el_it  =
    _mesh.active_local_elements_begin();
  const MeshBase::element_iterator el_end =
    _mesh.active_local_elements_end();

  for (; el_it != el_end; ++el_it)
  {
    const unsigned int id  = (*el_it)->id();
    libmesh_assert_less (id, error_per_cell.size());

    error_max = std::max (error_max, error_per_cell[id]);
    error_min = std::min (error_min, error_per_cell[id]);
  }
  CommWorld.max(error_max);
  CommWorld.min(error_min);

  // Compute the cutoff values for coarsening and refinement
  const Real error_delta = (error_max - error_min);
  const Real parent_error_delta = parent_error_max - parent_error_min;

  const Real refine_cutoff  = (1.- _refine_fraction)*error_max;
  const Real coarsen_cutoff = _coarsen_fraction*error_delta + error_min;
  const Real parent_cutoff = _coarsen_fraction*parent_error_delta + error_min;

//   // Print information about the error
//   libMesh::out << " Error Information:"                     << std::endl
// 	    << " ------------------"                     << std::endl
// 	    << "   min:              " << error_min      << std::endl
// 	    << "   max:              " << error_max      << std::endl
// 	    << "   delta:            " << error_delta    << std::endl
// 	    << "     refine_cutoff:  " << refine_cutoff  << std::endl
// 	    << "     coarsen_cutoff: " << coarsen_cutoff << std::endl;



  // Loop over the elements and flag them for coarsening or
  // refinement based on the element error

  MeshBase::element_iterator       e_it  =
    _mesh.active_elements_begin();
  const MeshBase::element_iterator e_end =
    _mesh.active_elements_end();
  for (; e_it != e_end; ++e_it)
  {
//.........这里部分代码省略.........

Real ErrorVector::median() const
{
  ErrorVector ev = (*this);

  return ev.median();
}

//.........这里部分代码省略.........
    infile("max_nonlinear_iterations", 15);
  solver.relative_step_tolerance =
    infile("relative_step_tolerance", 1.e-3);
  solver.relative_residual_tolerance =
    infile("relative_residual_tolerance", 0.0);
  solver.absolute_residual_tolerance =
    infile("absolute_residual_tolerance", 0.0);
    
  // And the linear solver options
  solver.max_linear_iterations =
    infile("max_linear_iterations", 50000);
  solver.initial_linear_tolerance =
    infile("initial_linear_tolerance", 1.e-3);

  // Print information about the system to the screen.
  equation_systems.print_info();

  // Now we begin the timestep loop to compute the time-accurate
  // solution of the equations.
  for (unsigned int t_step=0; t_step != n_timesteps; ++t_step)
    {
      // A pretty update message
      std::cout << "\n\nSolving time step " << t_step << ", time = "
                << system.time << std::endl;

      // Adaptively solve the timestep
      unsigned int a_step = 0;
      for (; a_step != max_adaptivesteps; ++a_step)
        {
          system.solve();

          system.postprocess();

          ErrorVector error;

          AutoPtr<ErrorEstimator> error_estimator;

          // To solve to a tolerance in this problem we
          // need a better estimator than Kelly
          if (global_tolerance != 0.)
            {
              // We can't adapt to both a tolerance and a mesh
              // size at once
              libmesh_assert_equal_to (nelem_target, 0);

              UniformRefinementEstimator *u =
                new UniformRefinementEstimator;

              // The lid-driven cavity problem isn't in H1, so
              // lets estimate L2 error
              u->error_norm = L2;

              error_estimator.reset(u);
            }
          else
            {
              // If we aren't adapting to a tolerance we need a
              // target mesh size
              libmesh_assert_greater (nelem_target, 0);

              // Kelly is a lousy estimator to use for a problem
              // not in H1 - if we were doing more than a few
              // timesteps we'd need to turn off or limit the
              // maximum level of our adaptivity eventually
              error_estimator.reset(new KellyErrorEstimator);
            }

//.........这里部分代码省略.........
  // Solve this as a steady system
  system.time_solver = libmesh_make_unique<SteadySolver>(system);

  // Initialize the system
  equation_systems.init ();

  // Set the time stepping options
  system.deltat = deltat;

  // And the nonlinear solver options
  DiffSolver & solver = *(system.time_solver->diff_solver().get());
  solver.quiet = infile("solver_quiet", true);
  solver.verbose = !solver.quiet;
  solver.max_nonlinear_iterations = infile("max_nonlinear_iterations", 15);
  solver.relative_step_tolerance = infile("relative_step_tolerance", 1.e-3);
  solver.relative_residual_tolerance = infile("relative_residual_tolerance", 0.0);
  solver.absolute_residual_tolerance = infile("absolute_residual_tolerance", 0.0);

  // And the linear solver options
  solver.max_linear_iterations = infile("max_linear_iterations", 50000);
  solver.initial_linear_tolerance = infile("initial_linear_tolerance", 1.e-3);

  // Print information about the system to the screen.
  equation_systems.print_info();

  // Adaptively solve the steady solution
  unsigned int a_step = 0;
  for (; a_step != max_adaptivesteps; ++a_step)
    {
      system.solve();

      system.postprocess();

      ErrorVector error;

      std::unique_ptr<ErrorEstimator> error_estimator;

      // To solve to a tolerance in this problem we
      // need a better estimator than Kelly
      if (global_tolerance != 0.)
        {
          // We can't adapt to both a tolerance and a mesh
          // size at once
          libmesh_assert_equal_to (nelem_target, 0);

          UniformRefinementEstimator * u =
            new UniformRefinementEstimator;

          // The lid-driven cavity problem isn't in H1, so
          // lets estimate L2 error
          u->error_norm = L2;

          error_estimator.reset(u);
        }
      else
        {
          // If we aren't adapting to a tolerance we need a
          // target mesh size
          libmesh_assert_greater (nelem_target, 0);

          // Kelly is a lousy estimator to use for a problem
          // not in H1 - if we were doing more than a few
          // timesteps we'd need to turn off or limit the
          // maximum level of our adaptivity eventually
          error_estimator.reset(new KellyErrorEstimator);
        }

//.........这里部分代码省略.........
  solver->quiet = infile("solver_quiet", true);
  solver->verbose = !solver->quiet;
  solver->max_nonlinear_iterations =
    infile("max_nonlinear_iterations", 15);
  solver->relative_step_tolerance =
    infile("relative_step_tolerance", 1.e-3);
  solver->relative_residual_tolerance =
    infile("relative_residual_tolerance", 0.0);
  solver->absolute_residual_tolerance =
    infile("absolute_residual_tolerance", 0.0);

  // And the linear solver options
  solver->max_linear_iterations =
    infile("max_linear_iterations", 50000);
  solver->initial_linear_tolerance =
    infile("initial_linear_tolerance", 1.e-3);

  // Print information about the system to the screen.
  equation_systems.print_info();
  
  // Now we begin the timestep loop to compute the time-accurate
  // solution of the equations...not that this is transient, but eh, why not...
	for (unsigned int t_step=0; t_step != n_timesteps; ++t_step){
    // A pretty update message
    std::cout << "\n\nSolving time step " << t_step << ", time = "
              << system.time << std::endl;

    // Adaptively solve the timestep
    unsigned int a_step = 0;
    for (; a_step != max_adaptivesteps; ++a_step)
      {
        system.solve();
        system.postprocess();
        ErrorVector error;
        AutoPtr<ErrorEstimator> error_estimator;

        // To solve to a tolerance in this problem we
        // need a better estimator than Kelly
        if (global_tolerance != 0.)
          {
            // We can't adapt to both a tolerance and a mesh
            // size at once
            libmesh_assert_equal_to (nelem_target, 0);

            UniformRefinementEstimator *u =
              new UniformRefinementEstimator;

            // The lid-driven cavity problem isn't in H1, so
            // lets estimate L2 error
            u->error_norm = L2;

            error_estimator.reset(u);
          }
        else
          {
            // If we aren't adapting to a tolerance we need a
            // target mesh size
            libmesh_assert_greater (nelem_target, 0);

            // Kelly is a lousy estimator to use for a problem
            // not in H1 - if we were doing more than a few
            // timesteps we'd need to turn off or limit the
            // maximum level of our adaptivity eventually
            error_estimator.reset(new KellyErrorEstimator);
          }

void JumpErrorEstimator::estimate_error (const System& system,
					 ErrorVector& error_per_cell,
					 const NumericVector<Number>* solution_vector,
					 bool estimate_parent_error)
{
  START_LOG("estimate_error()", "JumpErrorEstimator");
  /*

  Conventions for assigning the direction of the normal:

  - e & f are global element ids

  Case (1.) Elements are at the same level, e<f
            Compute the flux jump on the face and
	    add it as a contribution to error_per_cell[e]
	    and error_per_cell[f]

                   ----------------------
		  |           |          |
		  |           |    f     |
		  |           |          |
		  |    e      |---> n    |
		  |           |          |
		  |           |          |
                   ----------------------


   Case (2.) The neighbor is at a higher level.
             Compute the flux jump on e's face and
	     add it as a contribution to error_per_cell[e]
	     and error_per_cell[f]

                   ----------------------
		  |     |     |          |
		  |     |  e  |---> n    |
		  |     |     |          |
		  |-----------|    f     |
		  |     |     |          |
		  |     |     |          |
		  |     |     |          |
                   ----------------------
  */

  // The current mesh
  const MeshBase& mesh = system.get_mesh();

  // The dimensionality of the mesh
  const unsigned int dim = mesh.mesh_dimension();

  // The number of variables in the system
  const unsigned int n_vars = system.n_vars();

  // The DofMap for this system
  const DofMap& dof_map = system.get_dof_map();

  // Resize the error_per_cell vector to be
  // the number of elements, initialize it to 0.
  error_per_cell.resize (mesh.max_elem_id());
  std::fill (error_per_cell.begin(), error_per_cell.end(), 0.);

  // Declare a vector of floats which is as long as
  // error_per_cell above, and fill with zeros.  This vector will be
  // used to keep track of the number of edges (faces) on each active
  // element which are either:
  // 1) an internal edge
  // 2) an edge on a Neumann boundary for which a boundary condition
  //    function has been specified.
  // The error estimator can be scaled by the number of flux edges (faces)
  // which the element actually has to obtain a more uniform measure
  // of the error.  Use floats instead of ints since in case 2 (above)
  // f gets 1/2 of a flux face contribution from each of his
  // neighbors
  std::vector<float> n_flux_faces (error_per_cell.size());

  // Prepare current_local_solution to localize a non-standard
  // solution vector if necessary
  if (solution_vector && solution_vector != system.solution.get())
    {
      NumericVector<Number>* newsol =
        const_cast<NumericVector<Number>*>(solution_vector);
      System &sys = const_cast<System&>(system);
      newsol->swap(*sys.solution);
      sys.update();
    }

  // Loop over all the variables in the system
  for (var=0; var<n_vars; var++)
    {
      // Possibly skip this variable
      if (error_norm.weight(var) == 0.0) continue;

      // The type of finite element to use for this variable
      const FEType& fe_type = dof_map.variable_type (var);

      // Finite element objects for the same face from
      // different sides
      fe_fine = FEBase::build (dim, fe_type);
      fe_coarse = FEBase::build (dim, fe_type);

      // Build an appropriate Gaussian quadrature rule
//.........这里部分代码省略.........