C++ LibMeshInit类代码示例

int main(int argc, char** argv)
{
  LibMeshInit init (argc, argv);

  Mesh mesh(init.comm());
  EquationSystems es(mesh);

  libMesh::out << "Usage: " << argv[0]
               << " mesh solution" << std::endl;

  libMesh::out << "Loading..." << std::endl;

  mesh.read(argv[1]);
  libMesh::out << "Loaded mesh " << argv[1] << std::endl;
  mesh.print_info();

  es.read(argv[2], EquationSystems::READ_HEADER |
          EquationSystems::READ_DATA |
          EquationSystems::READ_ADDITIONAL_DATA |
          EquationSystems::READ_BASIC_ONLY);
  libMesh::out << "Loaded solution " << argv[2] << std::endl;
  es.print_info();

  libMesh::out.precision(16);

  for (unsigned int i = 0; i != es.n_systems(); ++i)
    {
      System &sys = es.get_system(i);

      output_norms(sys, *sys.solution, std::string("solution"));
      for (unsigned int j = 0; j != sys.n_vectors(); ++j)
        output_norms(sys, sys.get_vector(j), sys.vector_name(j));
    }
}

int main (int argc, char ** argv)
{
  LibMeshInit init (argc, argv);

  if (libMesh::on_command_line("--help") || argc < 3)
    {
      libMesh::out << "Example: " << argv[0] << " --mesh=filename.e --n-procs='4 8 16' "
                                                "[--num-ghost-layers <n>] [--dry-run] [--ascii]\n\n"
                   << "--mesh             Full name of the mesh file to read in. \n"
                   << "--n-procs          Vector of number of processors.\n"
                   << "--num-ghost-layers Number of layers to ghost when partitioning (Default: 1).\n"
                   << "--dry-run          Only test the partitioning, don't write any files.\n"
                   << "--ascii            Write ASCII cpa files rather than binary cpr files.\n"
                   << std::endl;

      return 0;
    }

  std::string filename = libMesh::command_line_value("--mesh", std::string());

  std::vector<processor_id_type> all_n_procs;
  libMesh::command_line_vector("--n-procs", all_n_procs);

  unsigned int num_ghost_layers = libMesh::command_line_value("--num-ghost-layers", 1);

  ReplicatedMesh mesh(init.comm());

  // If the user has requested additional ghosted layers, we need to add a ghosting functor.
  DefaultCoupling default_coupling;
  if (num_ghost_layers > 1)
    {
      default_coupling.set_n_levels(num_ghost_layers);
      mesh.add_ghosting_functor(default_coupling);
    }

  libMesh::out << "Reading " << filename << std::endl;

  mesh.read(filename);

  for (const auto & n_procs : all_n_procs)
    {
      libMesh::out << "splitting " << n_procs << " ways..." << std::endl;

      auto cpr = split_mesh(mesh, n_procs);

      if (!libMesh::on_command_line("--dry-run"))
        {
          libMesh::out << "    * writing " << cpr->current_processor_ids().size() << " files per process..." << std::endl;

          const bool binary = !libMesh::on_command_line("--ascii");

          cpr->binary() = binary;
          std::ostringstream outputname;
          outputname << remove_extension(filename) << (binary ? ".cpr" : ".cpa");
          cpr->write(outputname.str());
        }
    }

  return 0;
}

int main(int argc, char** argv)
{
  // Initialize libMesh.
  LibMeshInit init (argc, argv);
  if (on_command_line("--help"))
    print_help(argc, argv);
  else
    {
#if !defined(LIBMESH_ENABLE_SECOND_DERIVATIVES)
      libmesh_example_requires(false, "--enable-second");
#elif !defined(LIBMESH_ENABLE_PERIODIC)
      libmesh_example_requires(false, "--enable-periodic");
#endif

      // This is a PETSc-specific solver
      libmesh_example_requires(libMesh::default_solver_package() == PETSC_SOLVERS, "--enable-petsc");

      const int dim = command_line_value("dim",1);

      // Skip higher-dimensional examples on a lower-dimensional libMesh build
      libmesh_example_requires(dim <= LIBMESH_DIM, "2D/3D support");

      Biharmonic* biharmonic;
      Biharmonic::Create(&biharmonic,init.comm());
      biharmonic->viewParameters();
      biharmonic->init();
      biharmonic->run();
      Biharmonic::Destroy(&biharmonic);
    }
  return 0;
}

int main (int argc, char** argv){

	LibMeshInit init (argc, argv); //initialize libmesh library

	std::cout << "Running " << argv[0];
	for (int i=1; i<argc; i++)
		std::cout << " " << argv[i];
	std::cout << std::endl << std::endl;

	Mesh mesh(init.comm());
	
	mesh.read("channel_long.exo");
	
	std::string stash_assign = "divvy.txt";
	std::ofstream output(stash_assign.c_str());

	MeshBase::element_iterator       elem_it  = mesh.elements_begin();
  const MeshBase::element_iterator elem_end = mesh.elements_end();
  for (; elem_it != elem_end; ++elem_it){
    Elem* elem = *elem_it;
    		
    if(output.is_open()){
    	output << elem->id() << " " << elem->subdomain_id() << "\n";
    }
  }
  output.close();
	
	return 0;
}

// Begin the main program.
int main (int argc, char** argv)
{
    // Initialize libMesh and any dependent libaries, like in example 2.
    LibMeshInit init (argc, argv);

    // This example requires Adaptive Mesh Refinement support - although
    // it only refines uniformly, the refinement code used is the same
    // underneath.  It also requires libmesh support for Triangle and
    // Tetgen, which means that libmesh must be configured with
    // --disable-strict-lgpl for this example to run.
#if !defined(LIBMESH_HAVE_TRIANGLE) || !defined(LIBMESH_HAVE_TETGEN) || !defined(LIBMESH_ENABLE_AMR)
    libmesh_example_requires(false, "--disable-strict-lgpl --enable-amr");
#else

    // Skip this 2D example if libMesh was compiled as 1D-only.
    libmesh_example_requires(2 <= LIBMESH_DIM, "2D support");

    // Read the mesh from file.  This is the coarse mesh that will be used
    // in example 10 to demonstrate adaptive mesh refinement.  Here we will
    // simply read it in and uniformly refine it 5 times before we compute
    // with it.
    Mesh mesh(init.comm());

    {
        unsigned int dim=2;

        if (argc == 3 && std::atoi(argv[2]) == 3)
        {
            libmesh_here();
            dim=3;
        }

        mesh.read ((dim==2) ? "mesh.xda" : "hybrid_3d.xda");
    }

    // Create a MeshRefinement object to handle refinement of our mesh.
    // This class handles all the details of mesh refinement and coarsening.
    MeshRefinement mesh_refinement (mesh);

    // Uniformly refine the mesh 4 times.  This is the
    // first time we use the mesh refinement capabilities
    // of the library.
    mesh_refinement.uniformly_refine (4);

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

    // integrate the desired function
    integrate_function (mesh);

    // All done.
    return 0;
#endif
}

int main (int argc, char** argv)
{
  // Initialize the library.  This is necessary because the library
  // may depend on a number of other libraries (i.e. MPI and PETSc)
  // that require initialization before use.  When the LibMeshInit
  // object goes out of scope, other libraries and resources are
  // finalized.
  LibMeshInit init (argc, argv);

  // Check for proper usage. The program is designed to be run
  // as follows:
  // ./ex1 -d DIM input_mesh_name [-o output_mesh_name]
  // where [output_mesh_name] is an optional parameter giving
  // a filename to write the mesh into.
  if (argc < 4)
    {
      if (libMesh::processor_id() == 0)
        std::cerr << "Usage: " << argv[0] << " -d 2 in.mesh [-o out.mesh]"
                  << std::endl;

      // This handy function will print the file name, line number,
      // and then abort.  Currently the library does not use C++
      // exception handling.
      libmesh_error();
    }

  // Get the dimensionality of the mesh from argv[2]
  const unsigned int dim = std::atoi(argv[2]);

  // Skip higher-dimensional examples on a lower-dimensional libMesh build
  libmesh_example_assert(dim <= LIBMESH_DIM, "2D/3D support");

  // Create a mesh, with dimension to be overridden later, on the
  // default MPI communicator.
  Mesh mesh(init.comm());

  // Read the input mesh.
  mesh.read (argv[3]);

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

  // Write the output mesh if the user specified an
  // output file name.
  if (argc >= 6 && std::string("-o") == argv[4])
    mesh.write (argv[5]);

  // All done.  libMesh objects are destroyed here.  Because the
  // LibMeshInit object was created first, its destruction occurs
  // last, and it's destructor finalizes any external libraries and
  // checks for leaked memory.
  return 0;
}

int main (int argc, char** argv){

	LibMeshInit init (argc, argv); //initialize libmesh library

	std::cout << "Running " << argv[0];
	for (int i=1; i<argc; i++)
		std::cout << " " << argv[i];
	std::cout << std::endl << std::endl;

	Mesh mesh(init.comm());
	
	GetPot infile("inputs.in");
	std::string find_mesh_here = infile("mesh","psiLF_mesh.xda");
    mesh.read(find_mesh_here);
    
    // Print information about the mesh to the screen.
    mesh.print_info();

    // Create an equation systems object.
    EquationSystems equation_systems (mesh);

    // Name system; need to have same name is system being read in
    //ConvDiff_MprimeSys & system = 
    //    equation_systems.add_system<ConvDiff_MprimeSys>("Diff_ConvDiff_MprimeSys"); //for diff-convdiff
    ConvDiff_MprimeSys & system = 
        equation_systems.add_system<ConvDiff_MprimeSys>("ConvDiff_MprimeSys"); //for scalar-field

    // Steady-state problem	
    system.time_solver = AutoPtr<TimeSolver>(new SteadySolver(system));

    // Read in all the equation systems data from the LF solve (system, solutions, rhs, etc)
    std::string find_psiLF_here = infile("psiLF_file","psiLF.xda");
    std::cout << "Looking for psiLF at: " << find_psiLF_here << "\n\n";

    equation_systems.read(find_psiLF_here, READ,
	    EquationSystems::READ_HEADER |
	    EquationSystems::READ_DATA |
	    EquationSystems::READ_ADDITIONAL_DATA);
	    
    if(equation_systems.n_systems() > 1){
        std::cout << "\nName of system being read in does not match..." << std::endl;
    }
    
#ifdef LIBMESH_HAVE_GMV
    GMVIO(equation_systems.get_mesh()).write_equation_systems(std::string("psi.gmv"), equation_systems);
#endif
#ifdef LIBMESH_HAVE_EXODUS_API
    ExodusII_IO (mesh).write_equation_systems("psi.exo",equation_systems);
#endif // #ifdef LIBMESH_HAVE_EXODUS_API

    return 0;
}

// Begin the main program.
int main (int argc, char** argv)
{
  // Initialize libMesh and any dependent libaries, like in example 2.
  LibMeshInit init (argc, argv);

  libmesh_example_assert(2 <= LIBMESH_DIM, "2D support");

  std::cout << "Triangulating an L-shaped domain with holes" << std::endl;

  // 1.) 2D triangulation of L-shaped domain with three holes of different shape
  triangulate_domain(init.comm());

  libmesh_example_assert(3 <= LIBMESH_DIM, "3D support");

  std::cout << "Tetrahedralizing a prismatic domain with a hole" << std::endl;

  // 2.) 3D tetrahedralization of rectangular domain with hole.
  tetrahedralize_domain(init.comm());

  return 0;
}

int main(int argc, char* argv[])
{
  LibMeshInit init (argc, argv);

#ifdef LIBMESH_HAVE_DTK

  Mesh from_mesh(init.comm());
  MeshTools::Generation::build_cube(from_mesh, 4, 4, 4, 0, 1, 0, 1, 0, 1, HEX8);
  from_mesh.print_info();
  EquationSystems from_es(from_mesh);
  System & from_sys = from_es.add_system<ExplicitSystem>("From");
  unsigned int from_var = from_sys.add_variable("from");
  from_sys.attach_init_function(initialize);
  from_es.init();

  ExodusII_IO(from_mesh).write_equation_systems("from.e", from_es);

  Mesh to_mesh(init.comm());
  MeshTools::Generation::build_cube(to_mesh, 5, 5, 5, 0, 1, 0, 1, 0, 1, TET4);
  to_mesh.print_info();
  EquationSystems to_es(to_mesh);
  System & to_sys = to_es.add_system<ExplicitSystem>("To");
  unsigned int to_var = to_sys.add_variable("to");
  to_es.init();

  DTKSolutionTransfer dtk_transfer(init.comm());

  dtk_transfer.transfer(from_sys.variable(from_var), to_sys.variable(to_var));

  to_es.update();
  ExodusII_IO(to_mesh).write_equation_systems("to.e", to_es);

#endif

  return 0;
}

int main (int argc, char ** argv)
{
  LibMeshInit init (argc, argv);

#ifndef LIBMESH_HAVE_PETSC_TAO

  libmesh_example_requires(false, "PETSc >= 3.5.0 with built-in TAO support");

#elif LIBMESH_USE_COMPLEX_NUMBERS

  // According to
  // http://www.mcs.anl.gov/research/projects/tao/documentation/installation.html
  // TAO & PETSc-complex are currently mutually exclusive
  libmesh_example_requires(false, "PETSc >= 3.5.0 with built-in TAO support & real-numbers only");

#endif

  // We use a 2D domain.
  libmesh_example_requires(LIBMESH_DIM > 1, "--disable-1D-only");

  // TAO is giving us problems in parallel?
  if (init.comm().size() != 1)
    {
      libMesh::out << "This example can currently only be run in serial." << std::endl;
      return 77;
    }

  GetPot infile("optimization_ex2.in");
  const std::string approx_order = infile("approx_order", "FIRST");
  const std::string fe_family = infile("fe_family", "LAGRANGE");
  const unsigned int n_elem = infile("n_elem", 10);

  Mesh mesh(init.comm());
  MeshTools::Generation::build_square (mesh,
                                       n_elem,
                                       n_elem,
                                       -1., 1.,
                                       -1., 1.,
                                       QUAD9);

  mesh.print_info();

  EquationSystems equation_systems (mesh);

  OptimizationSystem & system =
    equation_systems.add_system<OptimizationSystem> ("Optimization");

  // The default is to use PETSc/Tao solvers, but let the user change
  // the optimization solver package on the fly.
  {
    const std::string optimization_solver_type = infile("optimization_solver_type",
                                                        "PETSC_SOLVERS");
    SolverPackage sp = Utility::string_to_enum<SolverPackage>(optimization_solver_type);
    UniquePtr<OptimizationSolver<Number> > new_solver =
      OptimizationSolver<Number>::build(system, sp);
    system.optimization_solver.reset(new_solver.release());
  }

  // Set tolerances and maximum iteration counts directly on the optimization solver.
  system.optimization_solver->max_objective_function_evaluations = 128;
  system.optimization_solver->objective_function_relative_tolerance = 1.e-4;
  system.optimization_solver->verbose = true;

  AssembleOptimization assemble_opt(system);

  system.add_variable("u",
                      Utility::string_to_enum<Order>   (approx_order),
                      Utility::string_to_enum<FEFamily>(fe_family));

  system.optimization_solver->objective_object = &assemble_opt;
  system.optimization_solver->gradient_object = &assemble_opt;
  system.optimization_solver->hessian_object = &assemble_opt;
  system.optimization_solver->equality_constraints_object = &assemble_opt;
  system.optimization_solver->equality_constraints_jacobian_object = &assemble_opt;
  system.optimization_solver->inequality_constraints_object = &assemble_opt;
  system.optimization_solver->inequality_constraints_jacobian_object = &assemble_opt;
  system.optimization_solver->lower_and_upper_bounds_object = &assemble_opt;

  // system.matrix and system.rhs are used for the gradient and Hessian,
  // so in this case we add an extra matrix and vector to store A and F.
  // This makes it easy to write the code for evaluating the objective,
  // gradient, and hessian.
  system.add_matrix("A_matrix");
  system.add_vector("F_vector");
  assemble_opt.A_matrix = &system.get_matrix("A_matrix");
  assemble_opt.F_vector = &system.get_vector("F_vector");


  equation_systems.init();
  equation_systems.print_info();

  assemble_opt.assemble_A_and_F();

  {
    std::vector< std::set<numeric_index_type> > constraint_jac_sparsity;
    std::set<numeric_index_type> sparsity_row;
    sparsity_row.insert(17);
    constraint_jac_sparsity.push_back(sparsity_row);
    sparsity_row.clear();

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

// The main program.
int main(int argc, char ** argv)
{
  // Initialize libMesh.
  LibMeshInit init (argc, argv);

#if !defined(LIBMESH_HAVE_XDR)
  // We need XDR support to write out reduced bases
  libmesh_example_requires(false, "--enable-xdr");
#elif defined(LIBMESH_DEFAULT_SINGLE_PRECISION)
  // XDR binary support requires double precision
  libmesh_example_requires(false, "--disable-singleprecision");
#elif defined(LIBMESH_DEFAULT_TRIPLE_PRECISION)
  // I have no idea why long double isn't working here... [RHS]
  libmesh_example_requires(false, "double precision");
#endif

  // Eigen can take forever to solve the offline mode portion of this
  // example
  libmesh_example_requires(libMesh::default_solver_package() != EIGEN_SOLVERS, "--enable-petsc or --enable-laspack");

  // Trilinos reports "true residual is too large" in the offline mode
  // portion of this example
  libmesh_example_requires(libMesh::default_solver_package() != TRILINOS_SOLVERS, "--enable-petsc or --enable-laspack");

  // This example only works if libMesh was compiled for 3D
  const unsigned int dim = 3;
  libmesh_example_requires(dim == LIBMESH_DIM, "3D support");

  std::string parameters_filename = "reduced_basis_ex5.in";
  GetPot infile(parameters_filename);

  unsigned int n_elem_x = infile("n_elem_x", 0);
  unsigned int n_elem_y = infile("n_elem_y", 0);
  unsigned int n_elem_z = infile("n_elem_z", 0);
  Real x_size           = infile("x_size", 0.);
  Real y_size           = infile("y_size", 0.);
  Real z_size           = infile("z_size", 0.);

  bool store_basis_functions = infile("store_basis_functions", true);

  // Read the "online_mode" flag from the command line
  GetPot command_line(argc, argv);
  int online_mode = 0;
  if (command_line.search(1, "-online_mode"))
    online_mode = command_line.next(online_mode);


  Mesh mesh (init.comm(), dim);
  MeshTools::Generation::build_cube (mesh,
                                     n_elem_x,
                                     n_elem_y,
                                     n_elem_z,
                                     0., x_size,
                                     0., y_size,
                                     0., z_size,
                                     HEX8);

  // Let's add a node boundary condition so that we can impose a point
  // load.
  // Each processor should know about each boundary condition it can
  // see, so we loop over all elements, not just local elements.
  MeshBase::const_element_iterator       el     = mesh.elements_begin();
  const MeshBase::const_element_iterator end_el = mesh.elements_end();
  for ( ; el != end_el; ++el)
    {
      const Elem * elem = *el;

      unsigned int side_max_x = 0, side_max_y = 0, side_max_z = 0;
      bool found_side_max_x = false, found_side_max_y = false, found_side_max_z = false;
      for (unsigned int side=0; side<elem->n_sides(); side++)
        {
          if (mesh.get_boundary_info().has_boundary_id(elem, side, BOUNDARY_ID_MAX_X))
            {
              side_max_x = side;
              found_side_max_x = true;
            }

          if (mesh.get_boundary_info().has_boundary_id(elem, side, BOUNDARY_ID_MAX_Y))
            {
              side_max_y = side;
              found_side_max_y = true;
            }

          if (mesh.get_boundary_info().has_boundary_id(elem, side, BOUNDARY_ID_MAX_Z))
            {
              side_max_z = side;
              found_side_max_z = true;
            }
        }

      // If elem has sides on boundaries
      // BOUNDARY_ID_MAX_X, BOUNDARY_ID_MAX_Y, BOUNDARY_ID_MAX_Z
      // then let's set a node boundary condition
      if (found_side_max_x && found_side_max_y && found_side_max_z)
        {
          for (unsigned int n=0; n<elem->n_nodes(); n++)
            {
              if (elem->is_node_on_side(n, side_max_x) &&
                  elem->is_node_on_side(n, side_max_y) &&
//.........这里部分代码省略.........

int main (int argc, char** argv){

	LibMeshInit init (argc, argv); //initialize libmesh library

	std::cout << "Running " << argv[0];
	for (int i=1; i<argc; i++)
		std::cout << " " << argv[i];
	std::cout << std::endl << std::endl;

	Mesh mesh(init.comm());
	
	//T-channel
	//mesh.read("mesh.e");
	//MeshRefinement meshRefinement(mesh);
	//meshRefinement.uniformly_refine(0);
	
	//1D - for debugging
	//int n = 20;
	//MeshTools::Generation::build_line(mesh, n, 0.0, 1.0, EDGE2); //n linear elements from 0 to 1
	
	//nice geometry (straight channel) 
	MeshTools::Generation::build_square (mesh, 
																					250, 50,
                                         -0.0, 5.0,
                                         -0.0, 1.0,
                                         QUAD9);

	//read in subdomain assignments
	
	std::vector<double> prev_assign(mesh.n_elem(), 0.);
	std::string read_assign = "do_divvy.txt";
	if(FILE *fp=fopen(read_assign.c_str(),"r")){
		int flag = 1;
		int elemNum, assign;
		int ind = 0;
		while(flag != -1){
			flag = fscanf(fp, "%d %d",&elemNum,&assign);
			if(flag != -1){
				prev_assign[ind] = assign;
				ind += 1;
			}
		}
		fclose(fp);
	}
	

	//to stash subdomain assignments
	std::string stash_assign = "divvy.txt";
	std::ofstream output(stash_assign.c_str());

	MeshBase::element_iterator       elem_it  = mesh.elements_begin();
  const MeshBase::element_iterator elem_end = mesh.elements_end();
  for (; elem_it != elem_end; ++elem_it){
    Elem* elem = *elem_it;
    Point c = elem->centroid();
    //Point cshift1(c(0)-0.63, c(1)-0.5);

    elem->subdomain_id() = prev_assign[elem->id()];
    //TEST/DEBUG
    //if(fabs(c(0)-3.0) < 0.35 && fabs(c(1)-0.5) < 0.35)
    //  elem->subdomain_id() = 1;
    		
    if(output.is_open()){
    	output << elem->id() << " " << elem->subdomain_id() << "\n";
    }
  }
  output.close();


	EquationSystems equation_systems (mesh);

	#ifdef LIBMESH_HAVE_EXODUS_API
    ExodusII_IO (mesh).write_equation_systems("meep.exo",equation_systems);
  #endif // #ifdef LIBMESH_HAVE_EXODUS_API

	return 0;
} // end main

// The main program.
int main (int argc, char ** argv)
{
  // Initialize libMesh.
  LibMeshInit init (argc, argv);

  // Skip this 2D example if libMesh was compiled as 1D-only.
  libmesh_example_requires(2 <= LIBMESH_DIM, "2D support");

  // This example NaNs with the Eigen sparse linear solvers and
  // Trilinos solvers, but should work OK with either PETSc or
  // Laspack.
  libmesh_example_requires(libMesh::default_solver_package() != EIGEN_SOLVERS, "--enable-petsc or --enable-laspack");
  libmesh_example_requires(libMesh::default_solver_package() != TRILINOS_SOLVERS, "--enable-petsc or --enable-laspack");

  // Create a mesh, with dimension to be overridden later, distributed
  // across the default MPI communicator.
  Mesh mesh(init.comm());

  // Use the MeshTools::Generation mesh generator to create a uniform
  // 2D grid on the square [-1,1]^2.  We instruct the mesh generator
  // to build a mesh of 8x8 Quad9 elements in 2D.  Building these
  // higher-order elements allows us to use higher-order
  // approximation, as in example 3.
  MeshTools::Generation::build_square (mesh,
                                       20, 20,
                                       0., 1.,
                                       0., 1.,
                                       QUAD9);

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

  // Create an equation systems object.
  EquationSystems equation_systems (mesh);

  // Declare the system and its variables.
  // Creates a transient system named "Navier-Stokes"
  TransientLinearImplicitSystem & system =
    equation_systems.add_system<TransientLinearImplicitSystem> ("Navier-Stokes");

  // Add the variables "u" & "v" to "Navier-Stokes".  They
  // will be approximated using second-order approximation.
  system.add_variable ("u", SECOND);
  system.add_variable ("v", SECOND);

  // Add the variable "p" to "Navier-Stokes". This will
  // be approximated with a first-order basis,
  // providing an LBB-stable pressure-velocity pair.
  system.add_variable ("p", FIRST);

  // Add a scalar Lagrange multiplier to constrain the
  // pressure to have zero mean.
  system.add_variable ("alpha", FIRST, SCALAR);

  // Give the system a pointer to the matrix assembly
  // function.
  system.attach_assemble_function (assemble_stokes);

  // Initialize the data structures for the equation system.
  equation_systems.init ();

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

  // Create a performance-logging object for this example
  PerfLog perf_log("Systems Example 3");

  // Get a reference to the Stokes system to use later.
  TransientLinearImplicitSystem & navier_stokes_system =
    equation_systems.get_system<TransientLinearImplicitSystem>("Navier-Stokes");

  // Now we begin the timestep loop to compute the time-accurate
  // solution of the equations.
  const Real dt = 0.01;
  navier_stokes_system.time = 0.0;
  const unsigned int n_timesteps = 15;

  // The number of steps and the stopping criterion are also required
  // for the nonlinear iterations.
  const unsigned int n_nonlinear_steps = 15;
  const Real nonlinear_tolerance = 1.e-3;

  // We also set a standard linear solver flag in the EquationSystems object
  // which controls the maxiumum number of linear solver iterations allowed.
  equation_systems.parameters.set<unsigned int>("linear solver maximum iterations") = 250;

  // Tell the system of equations what the timestep is by using
  // the set_parameter function.  The matrix assembly routine can
  // then reference this parameter.
  equation_systems.parameters.set<Real> ("dt") = dt;

  // The first thing to do is to get a copy of the solution at
  // the current nonlinear iteration.  This value will be used to
  // determine if we can exit the nonlinear loop.
  UniquePtr<NumericVector<Number> >
    last_nonlinear_soln (navier_stokes_system.solution->clone());

  for (unsigned int t_step=0; t_step<n_timesteps; ++t_step)
    {
//.........这里部分代码省略.........

int main(int argc, char** argv)
{
  // Initialize the library.  This is necessary because the library
  // may depend on a number of other libraries (i.e. MPI and PETSc)
  // that require initialization before use.  When the LibMeshInit
  // object goes out of scope, other libraries and resources are
  // finalized.
  LibMeshInit init (argc, argv);

  // Skip adaptive examples on a non-adaptive libMesh build
#ifndef LIBMESH_ENABLE_AMR
  libmesh_example_requires(false, "--enable-amr");
#else

  // Create a mesh, with dimension to be overridden later, on the
  // default MPI communicator.
  Mesh mesh(init.comm());

  GetPot command_line (argc, argv);

  int n = 4;
  if ( command_line.search(1, "-n") )
    n = command_line.next(n);

  // Build a 1D mesh with 4 elements from x=0 to x=1, using
  // EDGE3 (i.e. quadratic) 1D elements. They are called EDGE3 elements
  // because a quadratic element contains 3 nodes.
  MeshTools::Generation::build_line(mesh,n,0.,1.,EDGE3);

  // Define the equation systems object and the system we are going
  // to solve. See Introduction Example 2 for more details.
  EquationSystems equation_systems(mesh);
  LinearImplicitSystem& system = equation_systems.add_system
    <LinearImplicitSystem>("1D");

  // Add a variable "u" to the system, using second-order approximation
  system.add_variable("u",SECOND);

  // Give the system a pointer to the matrix assembly function. This
  // will be called when needed by the library.
  system.attach_assemble_function(assemble_1D);

  // Define the mesh refinement object that takes care of adaptively
  // refining the mesh.
  MeshRefinement mesh_refinement(mesh);

  // These parameters determine the proportion of elements that will
  // be refined and coarsened. Any element within 30% of the maximum
  // error on any element will be refined, and any element within 30%
  // of the minimum error on any element might be coarsened
  mesh_refinement.refine_fraction()  = 0.7;
  mesh_refinement.coarsen_fraction() = 0.3;
  // We won't refine any element more than 5 times in total
  mesh_refinement.max_h_level()      = 5;

  // Initialize the data structures for the equation system.
  equation_systems.init();

  // Refinement parameters
  const unsigned int max_r_steps = 5; // Refine the mesh 5 times

  // Define the refinement loop
  for(unsigned int r_step=0; r_step<=max_r_steps; r_step++)
    {
      // Solve the equation system
      equation_systems.get_system("1D").solve();

      // We need to ensure that the mesh is not refined on the last iteration
      // of this loop, since we do not want to refine the mesh unless we are
      // going to solve the equation system for that refined mesh.
      if(r_step != max_r_steps)
        {
          // Error estimation objects, see Adaptivity Example 2 for details
          ErrorVector error;
          KellyErrorEstimator error_estimator;

          // Compute the error for each active element
          error_estimator.estimate_error(system, error);

          // Output error estimate magnitude
          libMesh::out << "Error estimate\nl2 norm = " << error.l2_norm() <<
            "\nmaximum = " << error.maximum() << std::endl;

          // Flag elements to be refined and coarsened
          mesh_refinement.flag_elements_by_error_fraction (error);

          // Perform refinement and coarsening
          mesh_refinement.refine_and_coarsen_elements();

          // Reinitialize the equation_systems object for the newly refined
          // mesh. One of the steps in this is project the solution onto the
          // new mesh
          equation_systems.reinit();
        }
    }

  // Construct gnuplot plotting object, pass in mesh, title of plot
  // and boolean to indicate use of grid in plot. The grid is used to
  // show the edges of each element in the mesh.
  GnuPlotIO plot(mesh,"Adaptivity Example 1", GnuPlotIO::GRID_ON);
//.........这里部分代码省略.........

// Begin the main program.  Note that this example only
// works correctly if complex numbers have been enabled
// in the library.  In order to link against the complex
// PETSc libraries, you must have built PETSc with the same
// C++ compiler that you used to build libMesh.  This is
// so that the name mangling will be the same for the
// routines in both libraries.
int main (int argc, char ** argv)
{
  // Initialize libraries, like in example 2.
  LibMeshInit init (argc, argv);

  // This example is designed for complex numbers.
#ifndef LIBMESH_USE_COMPLEX_NUMBERS
  libmesh_example_requires(false, "--enable-complex");
#else

  // Check for proper usage.
  if (argc < 3)
    libmesh_error_msg("Usage: " << argv[0] << " -f [frequency]");

  if (init.comm().size() > 1)
    {
      if (init.comm().rank() == 0)
        {
          libMesh::err << "ERROR: Skipping example 7.\n"
                       << "MeshData objects currently only work in serial."
                       << std::endl;
        }
      return 0;
    }

  // Tell the user what we are doing.
  else
    {
      libMesh::out << "Running " << argv[0];

      for (int i=1; i<argc; i++)
        libMesh::out << " " << argv[i];

      libMesh::out << std::endl << std::endl;
    }

  // Get the frequency from argv[2] as a <i>float</i>,
  // currently, solve for 1/3rd, 2/3rd and 1/1th of the given frequency
  const Real frequency_in = atof(argv[2]);
  const unsigned int n_frequencies = 3;

  // Skip this 2D example if libMesh was compiled as 1D-only.
  libmesh_example_requires(2 <= LIBMESH_DIM, "2D support");

  // Create a mesh, with dimension to be overridden later, distributed
  // across the default MPI communicator.
  Mesh mesh(init.comm());

  // Create a corresponding MeshData
  // and activate it. For more information on this object
  // cf. example 12.
  MeshData mesh_data(mesh);
  mesh_data.activate();

  // Read the mesh file. Here the file lshape.unv contains
  // an L--shaped domain in .unv format.
  mesh.read("lshape.unv", &mesh_data);

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

  // The load on the boundary of the domain is stored in
  // the .unv formated mesh data file lshape_data.unv.
  // At this, the data is given as complex valued normal
  // velocities.
  mesh_data.read("lshape_data.unv");

  // Print information about the mesh to the screen.
  mesh_data.print_info();

  // Create an equation systems object, which now handles
  // a frequency system, as opposed to previous examples.
  // Also pass a MeshData pointer so the data can be
  // accessed in the matrix and rhs assembly.
  EquationSystems equation_systems (mesh, &mesh_data);

  // Create a FrequencySystem named "Helmholtz" & store a
  // reference to it.
  FrequencySystem & f_system =
    equation_systems.add_system<FrequencySystem> ("Helmholtz");

  // Add the variable "p" to "Helmholtz".  "p"
  // will be approximated using second-order approximation.
  f_system.add_variable("p", SECOND);

  // Tell the frequency system about the two user-provided
  // functions.  In other circumstances, at least the
  // solve function has to be attached.
  f_system.attach_assemble_function (assemble_helmholtz);
  f_system.attach_solve_function    (add_M_C_K_helmholtz);

  // To enable the fast solution scheme, additional
  // <i>global</i> matrices and one global vector, all appropriately sized,
//.........这里部分代码省略.........