C++ START_LOG函数代码示例

void UNVIO::count_elements (std::istream& in_file)
{
  START_LOG("count_elements()","UNVIO");

  if (this->_n_elements != 0)
    {
      libMesh::err << "Error: Trying to scan elements twice!"
		    << std::endl;
      libmesh_error();
    }


  // Simply read the element
  // dataset for the @e only
  // purpose to count nodes!

  std::string data;
  unsigned int fe_id;

  while (!in_file.eof())
    {
      // read element label
      in_file >> data;

      // end of dataset?
      if (data == "-1")
	break;

      // read fe_id
      in_file >> fe_id;

      // Skip related data,
      // and node number list
      in_file.ignore (256,'\n');
      in_file.ignore (256,'\n');

      // For some elements the node numbers
      // are given more than one record

      // TET10 or QUAD9
      if (fe_id == 118 || fe_id == 300)
	  in_file.ignore (256,'\n');

      // HEX20
      if (fe_id == 116)
	{
	  in_file.ignore (256,'\n');
	  in_file.ignore (256,'\n');
	}

      this->_n_elements++;
    }


  if (in_file.eof())
    {
      libMesh::err << "ERROR: File ended before end of element dataset!"
		    << std::endl;
      libmesh_error();
    }

  if (this->verbose())
    libMesh::out << "  Elements: " << this->_n_elements << std::endl;

  STOP_LOG("count_elements()","UNVIO");
}

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 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;
  if (scale_by_n_flux_faces)
    n_flux_faces.resize(error_per_cell.size(), 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();
    }

  fine_context.reset(new FEMContext(system));
  coarse_context.reset(new FEMContext(system));

  // Loop over all the variables we've been requested to find jumps in, to
  // pre-request
  for (var=0; var<n_vars; var++)
    {
      // Possibly skip this variable
      if (error_norm.weight(var) == 0.0) continue;

      // FIXME: Need to generalize this to vector-valued elements. [PB]
      FEBase * side_fe = libmesh_nullptr;

      const std::set<unsigned char> & elem_dims =
        fine_context->elem_dimensions();

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

Real RBEvaluation::rb_solve(unsigned int N)
{
  START_LOG("rb_solve()", "RBEvaluation");

  if(N > get_n_basis_functions())
    libmesh_error_msg("ERROR: N cannot be larger than the number of basis functions in rb_solve");

  const RBParameters& mu = get_parameters();

  // Resize (and clear) the solution vector
  RB_solution.resize(N);

  // Assemble the RB system
  DenseMatrix<Number> RB_system_matrix(N,N);
  RB_system_matrix.zero();

  DenseMatrix<Number> RB_Aq_a;
  for(unsigned int q_a=0; q_a<rb_theta_expansion->get_n_A_terms(); q_a++)
    {
      RB_Aq_vector[q_a].get_principal_submatrix(N, RB_Aq_a);

      RB_system_matrix.add(rb_theta_expansion->eval_A_theta(q_a, mu), RB_Aq_a);
    }

  // Assemble the RB rhs
  DenseVector<Number> RB_rhs(N);
  RB_rhs.zero();

  DenseVector<Number> RB_Fq_f;
  for(unsigned int q_f=0; q_f<rb_theta_expansion->get_n_F_terms(); q_f++)
    {
      RB_Fq_vector[q_f].get_principal_subvector(N, RB_Fq_f);

      RB_rhs.add(rb_theta_expansion->eval_F_theta(q_f, mu), RB_Fq_f);
    }

  // Solve the linear system
  if(N > 0)
    {
      RB_system_matrix.lu_solve(RB_rhs, RB_solution);
    }

  // Evaluate RB outputs
  DenseVector<Number> RB_output_vector_N;
  for(unsigned int n=0; n<rb_theta_expansion->get_n_outputs(); n++)
    {
      RB_outputs[n] = 0.;
      for(unsigned int q_l=0; q_l<rb_theta_expansion->get_n_output_terms(n); q_l++)
        {
          RB_output_vectors[n][q_l].get_principal_subvector(N, RB_output_vector_N);
          RB_outputs[n] += rb_theta_expansion->eval_output_theta(n,q_l,mu)*RB_output_vector_N.dot(RB_solution);
        }
    }

  if(evaluate_RB_error_bound) // Calculate the error bounds
    {
      // Evaluate the dual norm of the residual for RB_solution_vector
      Real epsilon_N = compute_residual_dual_norm(N);

      // Get lower bound for coercivity constant
      const Real alpha_LB = get_stability_lower_bound();
      // alpha_LB needs to be positive to get a valid error bound
      libmesh_assert_greater ( alpha_LB, 0. );

      // Evaluate the (absolute) error bound
      Real abs_error_bound = epsilon_N / residual_scaling_denom(alpha_LB);

      // Now compute the output error bounds
      for(unsigned int n=0; n<rb_theta_expansion->get_n_outputs(); n++)
        {
          RB_output_error_bounds[n] = abs_error_bound * eval_output_dual_norm(n, mu);
        }

      STOP_LOG("rb_solve()", "RBEvaluation");

      return abs_error_bound;
    }
  else // Don't calculate the error bounds
    {
      STOP_LOG("rb_solve()", "RBEvaluation");
      // Just return -1. if we did not compute the error bound
      return -1.;
    }
}

void NloptOptimizationSolver<T>::solve ()
{
  START_LOG("solve()", "NloptOptimizationSolver");

  this->init ();

  unsigned int nlopt_size = this->system().solution->size();

  // We have to have an objective function
  libmesh_assert( this->objective_object );

  // Set routine for objective and (optionally) gradient evaluation
  {
    nlopt_result ierr =
      nlopt_set_min_objective(_opt,
                              __libmesh_nlopt_objective,
                              this);
    if (ierr < 0)
      libmesh_error_msg("NLopt failed to set min objective: " << ierr);
  }

  if (this->lower_and_upper_bounds_object)
    {
      // Need to actually compute the bounds vectors first
      this->lower_and_upper_bounds_object->lower_and_upper_bounds(this->system());

      std::vector<Real> nlopt_lb(nlopt_size);
      std::vector<Real> nlopt_ub(nlopt_size);
      for(unsigned int i=0; i<nlopt_size; i++)
        {
          nlopt_lb[i] = this->system().get_vector("lower_bounds")(i);
          nlopt_ub[i] = this->system().get_vector("upper_bounds")(i);
        }

      nlopt_set_lower_bounds(_opt, &nlopt_lb[0]);
      nlopt_set_upper_bounds(_opt, &nlopt_ub[0]);
    }

  // If we have an equality constraints object, tell NLopt about it.
  if (this->equality_constraints_object)
    {
      // NLopt requires a vector to specify the tolerance for each constraint.
      // NLopt makes a copy of this vector internally, so it's safe for us to
      // let it go out of scope.
      std::vector<double> equality_constraints_tolerances(this->system().C_eq->size(),
                                                          _constraints_tolerance);

      // It would be nice to call the C interface directly, at least it should return an error
      // code we could parse... unfortunately, there does not seem to be a way to extract
      // the underlying nlopt_opt object from the nlopt::opt class!
      nlopt_result ierr =
        nlopt_add_equality_mconstraint(_opt,
                                       equality_constraints_tolerances.size(),
                                       __libmesh_nlopt_equality_constraints,
                                       this,
                                       &equality_constraints_tolerances[0]);

      if (ierr < 0)
        libmesh_error_msg("NLopt failed to add equality constraint: " << ierr);
    }

  // If we have an inequality constraints object, tell NLopt about it.
  if (this->inequality_constraints_object)
    {
      // NLopt requires a vector to specify the tolerance for each constraint
      std::vector<double> inequality_constraints_tolerances(this->system().C_ineq->size(),
                                                            _constraints_tolerance);

      nlopt_add_inequality_mconstraint(_opt,
                                       inequality_constraints_tolerances.size(),
                                       __libmesh_nlopt_inequality_constraints,
                                       this,
                                       &inequality_constraints_tolerances[0]);
    }

  // Set a relative tolerance on the optimization parameters
  nlopt_set_ftol_rel(_opt, this->objective_function_relative_tolerance);

  // Set the maximum number of allowed objective function evaluations
  nlopt_set_maxeval(_opt, this->max_objective_function_evaluations);

  // Reset internal iteration counter
  this->_iteration_count = 0;

  // Perform the optimization
  std::vector<Real> x(nlopt_size);
  Real min_val = 0.;
  _result = nlopt_optimize(_opt, &x[0], &min_val);

  if (_result < 0)
    libMesh::out << "NLopt failed!" << std::endl;
  else
    libMesh::out << "NLopt obtained optimal value: "
                 << min_val
                 << " in "
                 << this->get_iteration_count()
                 << " iterations."
                 << std::endl;

  STOP_LOG("solve()", "NloptOptimizationSolver");
//.........这里部分代码省略.........

Real RBEIMConstruction::truth_solve(int plot_solution)
{
  START_LOG("truth_solve()", "RBEIMConstruction");

  int training_parameters_found_index = -1;
  if( _parametrized_functions_in_training_set_initialized )
    {
      // Check if parameters are in the training set. If so, we can just load the
      // solution from _parametrized_functions_in_training_set

      for(unsigned int i=0; i<get_n_training_samples(); i++)
        {
          if(get_parameters() == get_params_from_training_set(i))
            {
              training_parameters_found_index = i;
              break;
            }
        }
    }

  // If the parameters are in the training set, just copy the solution vector
  if(training_parameters_found_index >= 0)
    {
      *solution = *_parametrized_functions_in_training_set[training_parameters_found_index];
      update(); // put the solution into current_local_solution as well
    }
  // Otherwise, we have to compute the projection
  else
    {
      RBEIMEvaluation& eim_eval = cast_ref<RBEIMEvaluation&>(get_rb_evaluation());
      eim_eval.set_parameters( get_parameters() );

      // Compute truth representation via projection
      const MeshBase& mesh = this->get_mesh();

      UniquePtr<DGFEMContext> c = this->build_context();
      DGFEMContext &context  = cast_ref<DGFEMContext&>(*c);

      this->init_context(context);

      rhs->zero();

      MeshBase::const_element_iterator       el     = mesh.active_local_elements_begin();
      const MeshBase::const_element_iterator end_el = mesh.active_local_elements_end();

      for ( ; el != end_el; ++el)
        {
          context.pre_fe_reinit(*this, *el);
          context.elem_fe_reinit();

          // All variables should have the same quadrature rule, hence
          // we can get JxW and xyz based on first_elem_fe.
          FEBase* first_elem_fe = NULL;
          context.get_element_fe( 0, first_elem_fe );
          unsigned int n_qpoints = context.get_element_qrule().n_points();
          const std::vector<Real> &JxW = first_elem_fe->get_JxW();
          const std::vector<Point> &xyz = first_elem_fe->get_xyz();

          // Loop over qp before var because parametrized functions often use
          // some caching based on qp.
          for (unsigned int qp=0; qp<n_qpoints; qp++)
            {
              for (unsigned int var=0; var<n_vars(); var++)
                {
                  FEBase* elem_fe = NULL;
                  context.get_element_fe( var, elem_fe );
                  const std::vector<std::vector<Real> >& phi = elem_fe->get_phi();

                  DenseSubVector<Number>& subresidual_var = context.get_elem_residual( var );

                  unsigned int n_var_dofs = cast_int<unsigned int>(context.get_dof_indices( var ).size());

                  Number eval_result = eim_eval.evaluate_parametrized_function(var, xyz[qp], *(*el));
                  for (unsigned int i=0; i != n_var_dofs; i++)
                    subresidual_var(i) += JxW[qp] * eval_result * phi[i][qp];
                }
            }

          // Apply constraints, e.g. periodic constraints
          this->get_dof_map().constrain_element_vector(context.get_elem_residual(), context.get_dof_indices() );

          // Add element vector to global vector
          rhs->add_vector(context.get_elem_residual(), context.get_dof_indices() );
        }

      // Solve to find the best fit, then solution stores the truth representation
      // of the function to be approximated
      solve_for_matrix_and_rhs(*inner_product_solver, *inner_product_matrix, *rhs);

      if (assert_convergence)
        check_convergence(*inner_product_solver);
    }

  if(plot_solution > 0)
    {
#ifdef LIBMESH_HAVE_EXODUS_API
      ExodusII_IO(get_mesh()).write_equation_systems ("truth.exo",
                                                      this->get_equation_systems());
#endif
    }
//.........这里部分代码省略.........

void RBEIMEvaluation::legacy_read_offline_data_from_files(const std::string& directory_name,
                                                          bool read_error_bound_data,
                                                          const bool read_binary_data)
{
  START_LOG("legacy_read_offline_data_from_files()", "RBEIMEvaluation");

  Parent::legacy_read_offline_data_from_files(directory_name, read_error_bound_data);

  // First, find out how many basis functions we had when Greedy terminated
  // This was set in RBSystem::read_offline_data_from_files
  unsigned int n_bfs = this->get_n_basis_functions();

  // The writing mode: DECODE for binary, READ for ASCII
  XdrMODE mode = read_binary_data ? DECODE : READ;

  // The suffix to use for all the files that are written out
  const std::string suffix = read_binary_data ? ".xdr" : ".dat";

  // Stream for creating file names
  std::ostringstream file_name;

  // Read in the interpolation matrix
  file_name.str("");
  file_name << directory_name << "/interpolation_matrix" << suffix;
  assert_file_exists(file_name.str());

  Xdr interpolation_matrix_in(file_name.str(), mode);

  for(unsigned int i=0; i<n_bfs; i++)
    {
      for(unsigned int j=0; j<=i; j++)
        {
          Number value;
          interpolation_matrix_in >> value;
          interpolation_matrix(i,j) = value;
        }
    }
  interpolation_matrix_in.close();

  // Also, read in the "extra" row
  file_name.str("");
  file_name << directory_name << "/extra_interpolation_matrix_row" << suffix;
  assert_file_exists(file_name.str());

  Xdr extra_interpolation_matrix_row_in(file_name.str(), mode);

  for(unsigned int j=0; j<n_bfs; j++)
    {
      Number value;
      extra_interpolation_matrix_row_in >> value;
      extra_interpolation_matrix_row(j) = value;
    }
  extra_interpolation_matrix_row_in.close();

  // Next read in interpolation_points
  file_name.str("");
  file_name << directory_name << "/interpolation_points" << suffix;
  assert_file_exists(file_name.str());

  Xdr interpolation_points_in(file_name.str(), mode);

  for(unsigned int i=0; i<n_bfs; i++)
    {
      Real x_val, y_val, z_val = 0.;
      interpolation_points_in >> x_val;

      if(LIBMESH_DIM >= 2)
        interpolation_points_in >> y_val;

      if(LIBMESH_DIM >= 3)
        interpolation_points_in >> z_val;

      Point p(x_val, y_val, z_val);
      interpolation_points.push_back(p);
    }
  interpolation_points_in.close();

  // Also, read in the extra interpolation point
  file_name.str("");
  file_name << directory_name << "/extra_interpolation_point" << suffix;
  assert_file_exists(file_name.str());

  Xdr extra_interpolation_point_in(file_name.str(), mode);

  {
    Real x_val, y_val, z_val = 0.;
    extra_interpolation_point_in >> x_val;

    if(LIBMESH_DIM >= 2)
      extra_interpolation_point_in >> y_val;

    if(LIBMESH_DIM >= 3)
      extra_interpolation_point_in >> z_val;

    Point p(x_val, y_val, z_val);
    extra_interpolation_point = p;
  }
  extra_interpolation_point_in.close();


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

void __libmesh_nlopt_inequality_constraints(unsigned m,
                                            double * result,
                                            unsigned n,
                                            const double * x,
                                            double * gradient,
                                            void * data)
{
  START_LOG("inequality_constraints()", "NloptOptimizationSolver");

  libmesh_assert(data);

  // data should be a pointer to the solver (it was passed in as void *)
  NloptOptimizationSolver<Number> * solver =
    static_cast<NloptOptimizationSolver<Number> *> (data);

  OptimizationSystem & sys = solver->system();

  // We'll use current_local_solution below, so let's ensure that it's consistent
  // with the vector x that was passed in.
  if (sys.solution->size() != n)
    libmesh_error_msg("Error: Input vector x has different length than sys.solution!");

  for(unsigned int i=sys.solution->first_local_index(); i<sys.solution->last_local_index(); i++)
    sys.solution->set(i, x[i]);
  sys.solution->close();

  // Impose constraints on the solution vector
  sys.get_dof_map().enforce_constraints_exactly(sys);

  // Update sys.current_local_solution based on the solution vector
  sys.update();

  // Call the user's inequality constraints function if there is one.
  OptimizationSystem::ComputeInequalityConstraints * ineco = solver->inequality_constraints_object;
  if (ineco)
    {
      ineco->inequality_constraints(*sys.current_local_solution,
                                    *sys.C_ineq,
                                    sys);

      sys.C_ineq->close();

      // Copy the values out of ineq_constraints into 'result'.
      // TODO: Even better would be if we could use 'result' directly
      // as the storage of ineq_constraints.  Perhaps a serial-only
      // NumericVector variant which supports this option?
      for (unsigned i=0; i<m; ++i)
        result[i] = (*sys.C_ineq)(i);

      // If gradient != NULL, then the Jacobian matrix of the equality
      // constraints has been requested.  The incoming 'gradient'
      // array is of length m*n and d(c_i)/d(x_j) = gradient[n*i+j].
      if (gradient)
        {
          OptimizationSystem::ComputeInequalityConstraintsJacobian * ineco_jac =
            solver->inequality_constraints_jacobian_object;

          if (ineco_jac)
            {
              ineco_jac->inequality_constraints_jacobian(*sys.current_local_solution,
                                                         *sys.C_ineq_jac,
                                                         sys);

              sys.C_ineq_jac->close();

              // copy the Jacobian data to the gradient array
              for(numeric_index_type i=0; i<m; i++)
                {
                  std::set<numeric_index_type>::iterator it = sys.ineq_constraint_jac_sparsity[i].begin();
                  std::set<numeric_index_type>::iterator it_end = sys.ineq_constraint_jac_sparsity[i].end();
                  for( ; it != it_end; ++it)
                    {
                      numeric_index_type dof_index = *it;
                      gradient[n*i+dof_index] = (*sys.C_ineq_jac)(i,dof_index);
                    }
                }
            }
          else
            libmesh_error_msg("Jacobian function not defined in __libmesh_nlopt_inequality_constraints");
        }

    }
  else
    libmesh_error_msg("Constraints function not defined in __libmesh_nlopt_inequality_constraints");

  STOP_LOG("inequality_constraints()", "NloptOptimizationSolver");
}

void InfFE<Dim,T_radial,T_map>::init_shape_functions(const Elem * inf_elem)
{
  libmesh_assert(inf_elem);


  // Start logging the radial shape function initialization
  START_LOG("init_shape_functions()", "InfFE");


  // -----------------------------------------------------------------
  // fast access to some const int's for the radial data
  const unsigned int n_radial_mapping_sf =
    cast_int<unsigned int>(radial_map.size());
  const unsigned int n_radial_approx_sf  =
    cast_int<unsigned int>(mode.size());
  const unsigned int n_radial_qp         =
    cast_int<unsigned int>(som.size());


  // -----------------------------------------------------------------
  // initialize most of the things related to mapping

  // The element type and order to use in the base map
  const Order    base_mapping_order     ( base_elem->default_order() );
  const ElemType base_mapping_elem_type ( base_elem->type()          );

  // the number of base shape functions used to construct the map
  // (Lagrange shape functions are used for mapping in the base)
  unsigned int n_base_mapping_shape_functions = Base::n_base_mapping_sf(base_mapping_elem_type,
                                                                        base_mapping_order);

  const unsigned int n_total_mapping_shape_functions =
    n_radial_mapping_sf * n_base_mapping_shape_functions;



  // -----------------------------------------------------------------
  // initialize most of the things related to physical approximation

  unsigned int n_base_approx_shape_functions;
  if (Dim > 1)
    n_base_approx_shape_functions = base_fe->n_shape_functions();
  else
    n_base_approx_shape_functions = 1;


  const unsigned int n_total_approx_shape_functions =
    n_radial_approx_sf * n_base_approx_shape_functions;

  // update class member field
  _n_total_approx_sf = n_total_approx_shape_functions;


  // The number of the base quadrature points.
  const unsigned int        n_base_qp =  base_qrule->n_points();

  // The total number of quadrature points.
  const unsigned int        n_total_qp =  n_radial_qp * n_base_qp;


  // update class member field
  _n_total_qp = n_total_qp;



  // -----------------------------------------------------------------
  // initialize the node and shape numbering maps
  {
    // these vectors work as follows: the i-th entry stores
    // the associated base/radial node number
    _radial_node_index.resize    (n_total_mapping_shape_functions);
    _base_node_index.resize      (n_total_mapping_shape_functions);

    // similar for the shapes: the i-th entry stores
    // the associated base/radial shape number
    _radial_shape_index.resize   (n_total_approx_shape_functions);
    _base_shape_index.resize     (n_total_approx_shape_functions);

    const ElemType inf_elem_type (inf_elem->type());

    // fill the node index map
    for (unsigned int n=0; n<n_total_mapping_shape_functions; n++)
      {
        compute_node_indices (inf_elem_type,
                              n,
                              _base_node_index[n],
                              _radial_node_index[n]);
        libmesh_assert_less (_base_node_index[n], n_base_mapping_shape_functions);
        libmesh_assert_less (_radial_node_index[n], n_radial_mapping_sf);
      }

    // fill the shape index map
    for (unsigned int n=0; n<n_total_approx_shape_functions; n++)
      {
        compute_shape_indices (this->fe_type,
                               inf_elem_type,
                               n,
                               _base_shape_index[n],
                               _radial_shape_index[n]);
        libmesh_assert_less (_base_shape_index[n], n_base_approx_shape_functions);
//.........这里部分代码省略.........

void FEMap::compute_edge_map(int dim,
			     const std::vector<Real>& qw,
			     const Elem* edge)
{
  libmesh_assert (edge != NULL);

  if (dim == 2)
    {
      // A 2D finite element living in either 2D or 3D space.
      // The edges here are the sides of the element, so the
      // (misnamed) compute_face_map function does what we want
      this->compute_face_map(dim, qw, edge);
      return;
    }

  libmesh_assert (dim == 3);  // 1D is unnecessary and currently unsupported

  START_LOG("compute_edge_map()", "FEMap");

  // The number of quadrature points.
  const unsigned int n_qp = qw.size();

  // Resize the vectors to hold data at the quadrature points
  this->xyz.resize(n_qp);
  this->dxyzdxi_map.resize(n_qp);
  this->dxyzdeta_map.resize(n_qp);
  this->d2xyzdxi2_map.resize(n_qp);
  this->d2xyzdxideta_map.resize(n_qp);
  this->d2xyzdeta2_map.resize(n_qp);
  this->tangents.resize(n_qp);
  this->normals.resize(n_qp);
  this->curvatures.resize(n_qp);

  this->JxW.resize(n_qp);

  // Clear the entities that will be summed
  for (unsigned int p=0; p<n_qp; p++)
    {
      this->tangents[p].resize(1);
      this->xyz[p].zero();
      this->dxyzdxi_map[p].zero();
      this->dxyzdeta_map[p].zero();
      this->d2xyzdxi2_map[p].zero();
      this->d2xyzdxideta_map[p].zero();
      this->d2xyzdeta2_map[p].zero();
    }

  // compute x, dxdxi at the quadrature points
  for (unsigned int i=0; i<this->psi_map.size(); i++) // sum over the nodes
    {
      const Point& edge_point = edge->point(i);

      for (unsigned int p=0; p<n_qp; p++) // for each quadrature point...
        {
	  this->xyz[p].add_scaled             (edge_point, this->psi_map[i][p]);
	  this->dxyzdxi_map[p].add_scaled     (edge_point, this->dpsidxi_map[i][p]);
	  this->d2xyzdxi2_map[p].add_scaled   (edge_point, this->d2psidxi2_map[i][p]);
        }
    }

  // Compute the tangents at the quadrature point
  // FIXME: normals (plural!) and curvatures are uncalculated
  for (unsigned int p=0; p<n_qp; p++)
    {
      const Point n  = this->dxyzdxi_map[p].cross(this->dxyzdeta_map[p]);
      this->tangents[p][0] = this->dxyzdxi_map[p].unit();

      // compute the jacobian at the quadrature points
      const Real jac = std::sqrt(this->dxdxi_map(p)*this->dxdxi_map(p) +
				 this->dydxi_map(p)*this->dydxi_map(p) +
				 this->dzdxi_map(p)*this->dzdxi_map(p));

      libmesh_assert (jac > 0.);

      this->JxW[p] = jac*qw[p];
    }

  STOP_LOG("compute_edge_map()", "FEMap");
}

void InfFE<Dim,T_radial,T_map>::init_radial_shape_functions(const Elem * libmesh_dbg_var(inf_elem))
{
  libmesh_assert(radial_qrule);
  libmesh_assert(inf_elem);


  /**
   * Start logging the radial shape function initialization
   */
  START_LOG("init_radial_shape_functions()", "InfFE");


  // -----------------------------------------------------------------
  // initialize most of the things related to mapping

  // The order to use in the radial map (currently independent of the element type)
  const Order        radial_mapping_order             (Radial::mapping_order());
  const unsigned int n_radial_mapping_shape_functions (Radial::n_dofs(radial_mapping_order));



  // -----------------------------------------------------------------
  // initialize most of the things related to physical approximation

  const Order        radial_approx_order             (fe_type.radial_order);
  const unsigned int n_radial_approx_shape_functions (Radial::n_dofs(radial_approx_order));

  const unsigned int        n_radial_qp = radial_qrule->n_points();
  const std::vector<Point> &   radial_qp = radial_qrule->get_points();



  // -----------------------------------------------------------------
  // resize the radial data fields

  mode.resize      (n_radial_approx_shape_functions);       // the radial polynomials (eval)
  dmodedv.resize   (n_radial_approx_shape_functions);

  som.resize       (n_radial_qp);                           // the (1-v)/2 weight
  dsomdv.resize    (n_radial_qp);

  radial_map.resize    (n_radial_mapping_shape_functions);  // the radial map
  dradialdv_map.resize (n_radial_mapping_shape_functions);


  for (unsigned int i=0; i<n_radial_mapping_shape_functions; i++)
    {
      radial_map[i].resize    (n_radial_qp);
      dradialdv_map[i].resize (n_radial_qp);
    }


  for (unsigned int i=0; i<n_radial_approx_shape_functions; i++)
    {
      mode[i].resize    (n_radial_qp);
      dmodedv[i].resize (n_radial_qp);
    }


  // compute scalar values at radial quadrature points
  for (unsigned int p=0; p<n_radial_qp; p++)
    {
      som[p]       = Radial::decay       (radial_qp[p](0));
      dsomdv[p]    = Radial::decay_deriv (radial_qp[p](0));
    }


  // evaluate the mode shapes in radial direction at radial quadrature points
  for (unsigned int i=0; i<n_radial_approx_shape_functions; i++)
    for (unsigned int p=0; p<n_radial_qp; p++)
      {
        mode[i][p]    = InfFE<Dim,T_radial,T_map>::eval       (radial_qp[p](0), radial_approx_order, i);
        dmodedv[i][p] = InfFE<Dim,T_radial,T_map>::eval_deriv (radial_qp[p](0), radial_approx_order, i);
      }


  // evaluate the mapping functions in radial direction at radial quadrature points
  for (unsigned int i=0; i<n_radial_mapping_shape_functions; i++)
    for (unsigned int p=0; p<n_radial_qp; p++)
      {
        radial_map[i][p]    = InfFE<Dim,INFINITE_MAP,T_map>::eval       (radial_qp[p](0), radial_mapping_order, i);
        dradialdv_map[i][p] = InfFE<Dim,INFINITE_MAP,T_map>::eval_deriv (radial_qp[p](0), radial_mapping_order, i);
      }

  /**
   * Stop logging the radial shape function initialization
   */
  STOP_LOG("init_radial_shape_functions()", "InfFE");

}

Real RBEIMEvaluation::rb_solve(unsigned int N)
{
  // Short-circuit if we are using the same parameters and value of N
  if( (_previous_parameters == get_parameters()) &&
      (_previous_N == N) )
    {
      return _previous_error_bound;
    }

  // Otherwise, update _previous parameters, _previous_N
  _previous_parameters = get_parameters();
  _previous_N = N;

  START_LOG("rb_solve()", "RBEIMEvaluation");

  if(N > get_n_basis_functions())
    libmesh_error_msg("ERROR: N cannot be larger than the number of basis functions in rb_solve");

  if(N==0)
    libmesh_error_msg("ERROR: N must be greater than 0 in rb_solve");

  // Get the rhs by sampling parametrized_function
  // at the first N interpolation_points
  DenseVector<Number> EIM_rhs(N);
  for(unsigned int i=0; i<N; i++)
    {
      EIM_rhs(i) = evaluate_parametrized_function(interpolation_points_var[i],
                                                  interpolation_points[i],
                                                  *interpolation_points_elem[i]);
    }



  DenseMatrix<Number> interpolation_matrix_N;
  interpolation_matrix.get_principal_submatrix(N, interpolation_matrix_N);

  interpolation_matrix_N.lu_solve(EIM_rhs, RB_solution);

  // Evaluate an a posteriori error bound
  if(evaluate_RB_error_bound)
    {
      // Compute the a posteriori error bound
      // First, sample the parametrized function at x_{N+1}
      Number g_at_next_x;
      if(N == get_n_basis_functions())
        g_at_next_x = evaluate_parametrized_function(extra_interpolation_point_var,
                                                     extra_interpolation_point,
                                                     *extra_interpolation_point_elem);
      else
        g_at_next_x = evaluate_parametrized_function(interpolation_points_var[N],
                                                     interpolation_points[N],
                                                     *interpolation_points_elem[N]);

      // Next, evaluate the EIM approximation at x_{N+1}
      Number EIM_approx_at_next_x = 0.;
      for(unsigned int j=0; j<N; j++)
        {
          if(N == get_n_basis_functions())
            {
              EIM_approx_at_next_x += RB_solution(j) * extra_interpolation_matrix_row(j);
            }
          else
            {
              EIM_approx_at_next_x += RB_solution(j) * interpolation_matrix(N,j);
            }
        }

      Real error_estimate = std::abs(g_at_next_x - EIM_approx_at_next_x);

      STOP_LOG("rb_solve()", "RBEIMEvaluation");

      _previous_error_bound = error_estimate;
      return error_estimate;
    }
  else // Don't evaluate an error bound
    {
      STOP_LOG("rb_solve()", "RBEIMEvaluation");
      _previous_error_bound = -1.;
      return -1.;
    }

}

void PointLocatorList::init ()
{
    libmesh_assert (!this->_list);

    if (this->_initialized)
    {
        libMesh::err << "ERROR: Already initialized!  Will ignore this call..."
                     << std::endl;
    }

    else

    {

        if (this->_master == NULL)
        {
            START_LOG("init(no master)", "PointLocatorList");

            // We are the master, so we have to build the list.
            // First create it, then get a handy reference, and
            // then try to speed up by reserving space...
            this->_list = new std::vector<std::pair<Point, const Elem *> >;
            std::vector<std::pair<Point, const Elem *> >& my_list = *(this->_list);

            my_list.clear();
            my_list.reserve(this->_mesh.n_active_elem());

            // fill our list with the centroids and element
            // pointers of the mesh.  For this use the handy
            // element iterators.
// 	  const_active_elem_iterator       el (this->_mesh.elements_begin());
// 	  const const_active_elem_iterator end(this->_mesh.elements_end());

            MeshBase::const_element_iterator       el  = _mesh.active_elements_begin();
            const MeshBase::const_element_iterator end = _mesh.active_elements_end();

            for (; el!=end; ++el)
                my_list.push_back(std::make_pair((*el)->centroid(), *el));

            STOP_LOG("init(no master)", "PointLocatorList");
        }

        else

        {
            // We are _not_ the master.  Let our _list point to
            // the master's list.  But for this we first transform
            // the master in a state for which we are friends
            // (this should also beware of a bad master pointer?).
            // And make sure the master @e has a list!
            const PointLocatorList* my_master =
                libmesh_cast_ptr<const PointLocatorList*>(this->_master);

            if (my_master->initialized())
                this->_list = my_master->_list;
            else
            {
                libMesh::err << "ERROR: Initialize master first, then servants!"
                             << std::endl;
                libmesh_error();
            }
        }

    }


    // ready for take-off
    this->_initialized = true;
}

void ParmetisPartitioner::_do_repartition (MeshBase& mesh,
					   const unsigned int n_sbdmns)
{
  libmesh_assert_greater (n_sbdmns, 0);

  // Check for an easy return
  if (n_sbdmns == 1)
    {
      this->single_partition(mesh);
      return;
    }

  // This function must be run on all processors at once
  parallel_only();

// What to do if the Parmetis library IS NOT present
#ifndef LIBMESH_HAVE_PARMETIS

  libmesh_here();
  libMesh::err << "ERROR: The library has been built without"  << std::endl
	        << "Parmetis support.  Using a Metis"           << std::endl
	        << "partitioner instead!"                       << std::endl;

  MetisPartitioner mp;

  mp.partition (mesh, n_sbdmns);

// What to do if the Parmetis library IS present
#else

  // Revert to METIS on one processor.
  if (libMesh::n_processors() == 1)
    {
      MetisPartitioner mp;
      mp.partition (mesh, n_sbdmns);
      return;
    }

  START_LOG("repartition()", "ParmetisPartitioner");

  // Initialize the data structures required by ParMETIS
  this->initialize (mesh, n_sbdmns);

  // Make sure all processors have enough active local elements.
  // Parmetis tends to crash when it's given only a couple elements
  // per partition.
  {
    bool all_have_enough_elements = true;
    for (unsigned int pid=0; pid<_n_active_elem_on_proc.size(); pid++)
      if (_n_active_elem_on_proc[pid] < MIN_ELEM_PER_PROC)
	all_have_enough_elements = false;

    // Parmetis will not work unless each processor has some
    // elements. Specifically, it will abort when passed a NULL
    // partition array on *any* of the processors.
    if (!all_have_enough_elements)
      {
	// FIXME: revert to METIS, although this requires a serial mesh
        MeshSerializer serialize(mesh);

	STOP_LOG ("repartition()", "ParmetisPartitioner");

	MetisPartitioner mp;
	mp.partition (mesh, n_sbdmns);

	return;
      }
  }

  // build the graph corresponding to the mesh
  this->build_graph (mesh);


  // Partition the graph
  std::vector<int> vsize(_vwgt.size(), 1);
  float itr = 1000000.0;
  MPI_Comm mpi_comm = libMesh::COMM_WORLD;

  // Call the ParMETIS adaptive repartitioning method.  This respects the
  // original partitioning when computing the new partitioning so as to
  // minimize the required data redistribution.
  Parmetis::ParMETIS_V3_AdaptiveRepart(_vtxdist.empty() ? NULL : &_vtxdist[0],
				       _xadj.empty()    ? NULL : &_xadj[0],
				       _adjncy.empty()  ? NULL : &_adjncy[0],
				       _vwgt.empty()    ? NULL : &_vwgt[0],
				       vsize.empty()    ? NULL : &vsize[0],
				       NULL,
				       &_wgtflag,
				       &_numflag,
				       &_ncon,
				       &_nparts,
				       _tpwgts.empty()  ? NULL : &_tpwgts[0],
				       _ubvec.empty()   ? NULL : &_ubvec[0],
				       &itr,
				       &_options[0],
				       &_edgecut,
				       _part.empty()    ? NULL : &_part[0],
				       &mpi_comm);

  // Assign the returned processor ids
//.........这里部分代码省略.........

void StatisticsVector<T>::histogram(std::vector<unsigned int>& bin_members,
				    unsigned int n_bins)
{
  // Must have at least 1 bin
  libmesh_assert (n_bins>0);

  const unsigned int n   = this->size();

  std::sort(this->begin(), this->end());

  // The StatisticsVector can hold both integer and float types.
  // We will define all the bins, etc. using Reals.
  Real min      = static_cast<Real>(this->minimum());
  Real max      = static_cast<Real>(this->maximum());
  Real bin_size = (max - min) / static_cast<Real>(n_bins);

  START_LOG ("histogram()", "StatisticsVector");

  std::vector<Real> bin_bounds(n_bins+1);
  for (unsigned int i=0; i<bin_bounds.size(); i++)
    bin_bounds[i] = min + i * bin_size;

  // Give the last bin boundary a little wiggle room: we don't want
  // it to be just barely less than the max, otherwise our bin test below
  // may fail.
  bin_bounds.back() += 1.e-6 * bin_size;

  // This vector will store the number of members each bin has.
  bin_members.resize(n_bins);

  unsigned int data_index = 0;
  for (unsigned int j=0; j<bin_members.size(); j++) // bin vector indexing
    {
      // libMesh::out << "(debug) Filling bin " << j << std::endl;

      for (unsigned int i=data_index; i<n; i++) // data vector indexing
	{
	  //	libMesh::out << "(debug) Processing index=" << i << std::endl;
	  Real current_val = static_cast<Real>( (*this)[i] );

	  // There may be entries in the vector smaller than the value
	  // reported by this->minimum().  (e.g. inactive elements in an
	  // ErrorVector.)  We just skip entries like that.
	  if ( current_val < min )
	    {
	      // 	    libMesh::out << "(debug) Skipping entry v[" << i << "]="
	      // 		      << (*this)[i]
	      // 		      << " which is less than the min value: min="
	      // 		      << min << std::endl;
	      continue;
	    }

	  if ( current_val > bin_bounds[j+1] ) // if outside the current bin (bin[j] is bounded
	                                       // by bin_bounds[j] and bin_bounds[j+1])
	    {
	      // libMesh::out.precision(16);
	      // 	    libMesh::out.setf(std::ios_base::fixed);
	      // 	    libMesh::out << "(debug) (*this)[i]= " << (*this)[i]
	      // 		      << " is greater than bin_bounds[j+1]="
	      //		      << bin_bounds[j+1]	 << std::endl;
	      data_index = i; // start searching here for next bin
	      break; // go to next bin
	    }

	  // Otherwise, increment current bin's count
	  bin_members[j]++;
	  // libMesh::out << "(debug) Binned index=" << i << std::endl;
	}
    }

#ifdef DEBUG
  // Check the number of binned entries
  const unsigned int n_binned = std::accumulate(bin_members.begin(),
						bin_members.end(),
						static_cast<unsigned int>(0),
						std::plus<unsigned int>());

  if (n != n_binned)
    {
      libMesh::out << "Warning: The number of binned entries, n_binned="
		<< n_binned
		<< ", did not match the total number of entries, n="
		<< n << "." << std::endl;
      //libmesh_error();
    }
#endif


  STOP_LOG ("histogram()", "StatisticsVector");
}

double __libmesh_nlopt_objective(unsigned n,
                                 const double * x,
                                 double * gradient,
                                 void * data)
{
  START_LOG("objective()", "NloptOptimizationSolver");

  // ctx should be a pointer to the solver (it was passed in as void *)
  NloptOptimizationSolver<Number> * solver =
    static_cast<NloptOptimizationSolver<Number> *> (data);

  OptimizationSystem & sys = solver->system();

  // We'll use current_local_solution below, so let's ensure that it's consistent
  // with the vector x that was passed in.
  for (unsigned int i=sys.solution->first_local_index();
       i<sys.solution->last_local_index(); i++)
    sys.solution->set(i, x[i]);

  // Make sure the solution vector is parallel-consistent
  sys.solution->close();

  // Impose constraints on X
  sys.get_dof_map().enforce_constraints_exactly(sys);

  // Update sys.current_local_solution based on X
  sys.update();

  Real objective;
  if (solver->objective_object != NULL)
    {
      objective =
        solver->objective_object->objective(
                                            *(sys.current_local_solution), sys);
    }
  else
    {
      libmesh_error_msg("Objective function not defined in __libmesh_nlopt_objective");
    }

  // If the gradient has been requested, fill it in
  if (gradient)
    {
      if (solver->gradient_object != NULL)
        {
          solver->gradient_object->gradient(
                                            *(sys.current_local_solution), *(sys.rhs), sys);

          // we've filled up sys.rhs with the gradient data, now copy it
          // to the nlopt data structure
          libmesh_assert(sys.rhs->size() == n);

          std::vector<double> grad;
          sys.rhs->localize_to_one(grad);
          for (unsigned i=0; i<n; ++i)
            gradient[i] = grad[i];
        }
      else
        libmesh_error_msg("Gradient function not defined in __libmesh_nlopt_objective");
    }

  STOP_LOG("objective()", "NloptOptimizationSolver");

  // Increment the iteration count.
  solver->get_iteration_count()++;

  // Possibly print the current value of the objective function
  if (solver->verbose)
    libMesh::out << objective << std::endl;

  return objective;
}