C++ MultiLevelMesh::PrintInfo方法代码示例

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

  // init Petsc-MPI communicator
  FemusInit mpinit(argc, args, MPI_COMM_WORLD);

  // define multilevel mesh
  MultiLevelMesh mlMsh;
  double scalingFactor = 1.;
  //read coarse level mesh and generate finers level meshes
  //mlMsh.ReadCoarseMesh("./input/box.neu", "seventh", scalingFactor);
  //mlMsh.ReadCoarseMesh("./input/square.neu", "seventh", scalingFactor);
  //mlMsh.ReadCoarseMesh("./input/square_tri.neu", "seventh", scalingFactor);
  mlMsh.ReadCoarseMesh("./input/cube_mixed.neu", "seventh", scalingFactor);
  /* "seventh" is the order of accuracy that is used in the gauss integration scheme
      probably in the furure it is not going to be an argument of this function   */
  unsigned numberOfUniformLevels = 1;
  unsigned numberOfSelectiveLevels = 3;
  mlMsh.RefineMesh(numberOfUniformLevels + numberOfSelectiveLevels, numberOfUniformLevels , SetRefinementFlag);
  mlMsh.PrintInfo();

  // define the multilevel solution and attach the mlMsh object to it
  MultiLevelSolution mlSol(&mlMsh);

  // add variables to mlSol
  //mlSol.AddSolution("U", LAGRANGE, FIRST);
  mlSol.AddSolution("U", LAGRANGE, SECOND);
//   mlSol.AddSolution("V", LAGRANGE, SERENDIPITY);
//   mlSol.AddSolution("W", LAGRANGE, SECOND);
  mlSol.AddSolution("P", DISCONTINOUS_POLYNOMIAL, ZERO);
  mlSol.AddSolution("T", DISCONTINOUS_POLYNOMIAL, FIRST);
//
  mlSol.Initialize("All");    // initialize all varaibles to zero
//
  mlSol.Initialize("U", InitalValueU);
  mlSol.Initialize("P", InitalValueP);
  mlSol.Initialize("T", InitalValueT);    // note that this initialization is the same as piecewise constant element
//
//   // print solutions
  
  mlSol.AttachSetBoundaryConditionFunction(SetBoundaryCondition);
  mlSol.GenerateBdc("All");
  
  std::vector < std::string > variablesToBePrinted;
  variablesToBePrinted.push_back("U");
  variablesToBePrinted.push_back("P");
  variablesToBePrinted.push_back("T");

  VTKWriter vtkIO(&mlSol);
  vtkIO.write(DEFAULT_OUTPUTDIR, "biquadratic", variablesToBePrinted);

  GMVWriter gmvIO(&mlSol);
  variablesToBePrinted.push_back("all");
  gmvIO.SetDebugOutput(true);
  gmvIO.write(DEFAULT_OUTPUTDIR, "biquadratic", variablesToBePrinted);

  return 0;
}

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

  // init Petsc-MPI communicator
  FemusInit mpinit(argc, args, MPI_COMM_WORLD);


  // define multilevel mesh
  MultiLevelMesh mlMsh;
  // read coarse level mesh and generate finers level meshes
  double scalingFactor = 1.;
  mlMsh.ReadCoarseMesh("./input/square_quad.neu", "seventh", scalingFactor);
  /* "seventh" is the order of accuracy that is used in the gauss integration scheme
    probably in the future it is not going to be an argument of this function   */

  unsigned maxNumberOfMeshes = 5;

  vector < vector < double > > l2Norm;
  l2Norm.resize(maxNumberOfMeshes);

  vector < vector < double > > semiNorm;
  semiNorm.resize(maxNumberOfMeshes);

  for (unsigned i = 0; i < maxNumberOfMeshes; i++) {   // loop on the mesh level

    unsigned numberOfUniformLevels = i + 1;
    unsigned numberOfSelectiveLevels = 0;
    mlMsh.RefineMesh( numberOfUniformLevels + numberOfSelectiveLevels, numberOfUniformLevels , NULL);

    // erase all the coarse mesh levels
    mlMsh.EraseCoarseLevels(numberOfUniformLevels - 1);

    // print mesh info
    mlMsh.PrintInfo();

    FEOrder feOrder[3] = {FIRST, SERENDIPITY, SECOND};
    l2Norm[i].resize(3);
    semiNorm[i].resize(3);

    for (unsigned j = 0; j < 3; j++) {   // loop on the FE Order

      std::cout << "level = " << i << " FEM = " << j << std::endl;

      // define the multilevel solution and attach the mlMsh object to it
      MultiLevelSolution mlSol(&mlMsh);

      // add variables to mlSol
      mlSol.AddSolution("u", LAGRANGE, feOrder[j]);
      mlSol.AddSolution("v", LAGRANGE, feOrder[j]);
      mlSol.Initialize("All");

      // attach the boundary condition function and generate boundary data
      mlSol.AttachSetBoundaryConditionFunction(SetBoundaryCondition);
      mlSol.GenerateBdc("u");
      mlSol.GenerateBdc("v");

      // define the multilevel problem attach the mlSol object to it
      MultiLevelProblem mlProb(&mlSol);

      // add system Poisson in mlProb as a Linear Implicit System
      NonLinearImplicitSystem& system = mlProb.add_system < NonLinearImplicitSystem > ("Poisson");

      // add solution "u" to system
      system.AddSolutionToSystemPDE("u");
      system.AddSolutionToSystemPDE("v");

      // attach the assembling function to system
      system.SetAssembleFunction(AssembleBilaplaceProblem_AD);

      // initilaize and solve the system
      system.init();
      system.MGsolve();

      std::pair< double , double > norm = GetErrorNorm(&mlSol);
      l2Norm[i][j]  = norm.first;
      semiNorm[i][j] = norm.second;
      // print solutions
      std::vector < std::string > variablesToBePrinted;
      variablesToBePrinted.push_back("All");

      VTKWriter vtkIO(&mlSol);
      vtkIO.write(DEFAULT_OUTPUTDIR, "biquadratic", variablesToBePrinted, i);

    }
  }

  // print the seminorm of the error and the order of convergence between different levels
  std::cout << std::endl;
  std::cout << std::endl;
  std::cout << "l2 ERROR and ORDER OF CONVERGENCE:\n\n";
  std::cout << "LEVEL\tFIRST\t\t\tSERENDIPITY\t\tSECOND\n";

  for (unsigned i = 0; i < maxNumberOfMeshes; i++) {
    std::cout << i + 1 << "\t";
    std::cout.precision(14);

    for (unsigned j = 0; j < 3; j++) {
      std::cout << l2Norm[i][j] << "\t";
    }

    std::cout << std::endl;
//.........这里部分代码省略.........

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



  // init Petsc-MPI communicator
  FemusInit mpinit(argc, args, MPI_COMM_WORLD);

  // define multilevel mesh
  MultiLevelMesh mlMsh;
  // read coarse level mesh and generate finers level meshes
  double scalingFactor = 1.;
  mlMsh.ReadCoarseMesh("./input/cube_hex.neu", "seventh", scalingFactor);
  //mlMsh.ReadCoarseMesh ( "./input/square_quad.neu", "seventh", scalingFactor );
  /* "seventh" is the order of accuracy that is used in the gauss integration scheme
     probably in the furure it is not going to be an argument of this function   */
  unsigned dim = mlMsh.GetDimension();

  unsigned numberOfUniformLevels = 3;
  unsigned numberOfSelectiveLevels = 0;
  mlMsh.RefineMesh(numberOfUniformLevels , numberOfUniformLevels + numberOfSelectiveLevels, NULL);

  // erase all the coarse mesh levels
  mlMsh.EraseCoarseLevels(numberOfUniformLevels - 1);

  // print mesh info
  mlMsh.PrintInfo();

  MultiLevelSolution mlSol(&mlMsh);

  // add variables to mlSol
  mlSol.AddSolution("U", LAGRANGE, SECOND);
  mlSol.AddSolution("V", LAGRANGE, SECOND);

  if (dim == 3) mlSol.AddSolution("W", LAGRANGE, SECOND);

  mlSol.AddSolution("P", LAGRANGE, FIRST);
  mlSol.Initialize("All");

  // attach the boundary condition function and generate boundary data
  mlSol.AttachSetBoundaryConditionFunction(SetBoundaryCondition);
  mlSol.GenerateBdc("All");

  // define the multilevel problem attach the mlSol object to it
  MultiLevelProblem mlProb(&mlSol);

  // add system Poisson in mlProb as a Linear Implicit System
  NonLinearImplicitSystem& system = mlProb.add_system < NonLinearImplicitSystem > ("NS");

  // add solution "u" to system
  system.AddSolutionToSystemPDE("U");
  system.AddSolutionToSystemPDE("V");

  if (dim == 3) system.AddSolutionToSystemPDE("W");

  system.AddSolutionToSystemPDE("P");

  // attach the assembling function to system
  system.SetAssembleFunction(AssembleBoussinesqAppoximation_AD);

  // initilaize and solve the system
  system.init();
  system.solve();

  // print solutions
  std::vector < std::string > variablesToBePrinted;
  variablesToBePrinted.push_back("All");

  VTKWriter vtkIO(&mlSol);
  vtkIO.write(DEFAULT_OUTPUTDIR, "biquadratic", variablesToBePrinted);

  return 0;
}

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

  // init Petsc-MPI communicator
  FemusInit mpinit (argc, args, MPI_COMM_WORLD);

  // define multilevel mesh
  MultiLevelMesh mlMsh;
  // read coarse level mesh and generate finers level meshes
  double scalingFactor = 1.;
  //mlMsh.ReadCoarseMesh("./input/cube_hex.neu","seventh",scalingFactor);
  //mlMsh.ReadCoarseMesh("./input/square_quad.neu", "seventh", scalingFactor);
  mlMsh.ReadCoarseMesh ("./input/quadAMR.neu", "seventh", scalingFactor);
  /* "seventh" is the order of accuracy that is used in the gauss integration scheme
     probably in the furure it is not going to be an argument of this function   */
  unsigned dim = mlMsh.GetDimension();

//   unsigned numberOfUniformLevels = 3;
//   unsigned numberOfSelectiveLevels = 0;
//   mlMsh.RefineMesh(numberOfUniformLevels , numberOfUniformLevels + numberOfSelectiveLevels, NULL);

  unsigned numberOfUniformLevels = 4;
  unsigned numberOfSelectiveLevels = 3;
  mlMsh.RefineMesh (numberOfUniformLevels + numberOfSelectiveLevels, numberOfUniformLevels , SetRefinementFlag);

  // erase all the coarse mesh levels
  //mlMsh.EraseCoarseLevels(numberOfUniformLevels - 3);

  // print mesh info
  mlMsh.PrintInfo();

  MultiLevelSolution mlSol (&mlMsh);

  // add variables to mlSol


  mlSol.AddSolution ("U", LAGRANGE, SERENDIPITY);
  mlSol.AddSolution ("V", LAGRANGE, SECOND);


  mlSol.Initialize ("All");

  // attach the boundary condition function and generate boundary data
  mlSol.AttachSetBoundaryConditionFunction (SetBoundaryCondition);

  mlSol.GenerateBdc ("All");

  // define the multilevel problem attach the mlSol object to it
  MultiLevelProblem mlProb (&mlSol);

  // add system Poisson in mlProb as a Linear Implicit System
  NonLinearImplicitSystem& system = mlProb.add_system < NonLinearImplicitSystem > ("Poisson");

  // add solution "u" to system
  system.AddSolutionToSystemPDE ("U");
  system.AddSolutionToSystemPDE ("V");

  //system.SetLinearEquationSolverType(FEMuS_DEFAULT);
  system.SetLinearEquationSolverType (FEMuS_ASM); // Additive Swartz Method
  // attach the assembling function to system
  system.SetAssembleFunction (AssemblePoisson_AD);


  system.SetMaxNumberOfNonLinearIterations(10);
  system.SetMaxNumberOfLinearIterations(3);
  system.SetAbsoluteLinearConvergenceTolerance(1.e-12);
  system.SetNonLinearConvergenceTolerance(1.e-8);
  system.SetMgType(F_CYCLE); // Q1 What's F cycle

  system.SetNumberPreSmoothingStep (0);
  system.SetNumberPostSmoothingStep (2);
  // initilaize and solve the system
  system.init();

  system.SetSolverFineGrids (GMRES);
  system.SetPreconditionerFineGrids (ILU_PRECOND);
  system.SetTolerances (1.e-3, 1.e-20, 1.e+50, 5);

  system.SetNumberOfSchurVariables (1);
  system.SetElementBlockNumber (4);
  //system.SetDirichletBCsHandling(ELIMINATION);
  //system.solve();
  system.MGsolve();

  // print solutions
  std::vector < std::string > variablesToBePrinted;
  variablesToBePrinted.push_back ("All");

  VTKWriter vtkIO (&mlSol);
  vtkIO.Write (DEFAULT_OUTPUTDIR, "biquadratic", variablesToBePrinted);

  GMVWriter gmvIO (&mlSol);
  variablesToBePrinted.push_back ("all");
  gmvIO.SetDebugOutput (true);
  gmvIO.Write (DEFAULT_OUTPUTDIR, "biquadratic", variablesToBePrinted);


  return 0;
}

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

  SlepcInitialize(&argc, &args, PETSC_NULL, PETSC_NULL);

  // init Petsc-MPI communicator
  FemusInit mpinit(argc, args, MPI_COMM_WORLD);

  // define multilevel mesh
  MultiLevelMesh mlMsh;
  double scalingFactor = 1.;



  unsigned numberOfUniformLevels = 1;
  unsigned numberOfSelectiveLevels = 0;

  unsigned nx = static_cast<unsigned>(floor(pow(2.,11) + 0.5));
  
  double length = 2. * 1465700.;

  mlMsh.GenerateCoarseBoxMesh(nx, 0, 0, -length / 2, length / 2, 0., 0., 0., 0., EDGE3, "seventh");

  //mlMsh.RefineMesh(numberOfUniformLevels + numberOfSelectiveLevels, numberOfUniformLevels , NULL);
  mlMsh.PrintInfo();

  // define the multilevel solution and attach the mlMsh object to it
  MultiLevelSolution mlSol(&mlMsh);

  for(unsigned i = 0; i < NumberOfLayers; i++) {
    char name[10];
    sprintf(name, "h%d", i);
    mlSol.AddSolution(name, DISCONTINUOUS_POLYNOMIAL, ZERO);
    sprintf(name, "v%d", i);
    mlSol.AddSolution(name, LAGRANGE, FIRST);
    sprintf(name, "T%d", i);
    mlSol.AddSolution(name, DISCONTINUOUS_POLYNOMIAL, ZERO);
    sprintf(name, "HT%d", i);
    mlSol.AddSolution(name, DISCONTINUOUS_POLYNOMIAL, ZERO);
  }
  
  mlSol.AddSolution("b", DISCONTINUOUS_POLYNOMIAL, ZERO, 1, false);

  mlSol.AddSolution("eta", DISCONTINUOUS_POLYNOMIAL, ZERO, 1, false);

  mlSol.Initialize("All");
  
  
  mlSol.Initialize("h0",InitalValueH0);
  mlSol.Initialize("T0",InitalValueT0);
  if(NumberOfLayers > 1){
    mlSol.Initialize("h1",InitalValueH1);
    mlSol.Initialize("T1",InitalValueT1);
    if(NumberOfLayers > 2){
      mlSol.Initialize("h2",InitalValueH2);
      mlSol.Initialize("T2",InitalValueT2);
    }
  }
  
  for(unsigned i = 0; i < NumberOfLayers; i++) {
    char name[10];
    sprintf(name, "v%d", i);
    mlSol.Initialize(name, InitalValueV);
  }

  mlSol.Initialize("b", InitalValueB);
  
  mlSol.AttachSetBoundaryConditionFunction(SetBoundaryCondition);
  mlSol.GenerateBdc("All");

  MultiLevelProblem ml_prob(&mlSol);

  // ******* Add FEM system to the MultiLevel problem *******
  LinearImplicitSystem& system = ml_prob.add_system < LinearImplicitSystem > ("SW");
  for(unsigned i = 0; i < NumberOfLayers; i++) {
    char name[10];
    sprintf(name, "h%d", i);
    system.AddSolutionToSystemPDE(name);
    sprintf(name, "v%d", i);
    system.AddSolutionToSystemPDE(name);
    sprintf(name, "HT%d", i);
    system.AddSolutionToSystemPDE(name);
  }
  system.init();

  mlSol.SetWriter(VTK);
  std::vector<std::string> print_vars;
  print_vars.push_back("All");
  //mlSol.GetWriter()->SetDebugOutput(true);
  mlSol.GetWriter()->Write(DEFAULT_OUTPUTDIR, "linear", print_vars, 0);

  unsigned numberOfTimeSteps = 2000;
  for(unsigned i = 0; i < numberOfTimeSteps; i++) {
    ETD(ml_prob);
    mlSol.GetWriter()->Write(DEFAULT_OUTPUTDIR, "linear", print_vars, (i + 1)/1);
  }
  return 0;
}

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

  // init Petsc-MPI communicator
  FemusInit mpinit (argc, args, MPI_COMM_WORLD);

  // define multilevel mesh
  MultiLevelMesh mlMsh;
  // read coarse level mesh and generate finers level meshes
  double scalingFactor = 1.;
  //mlMsh.ReadCoarseMesh("./input/cube_hex.neu","seventh",scalingFactor);
  //mlMsh.ReadCoarseMesh("./input/square_quad.neu", "seventh", scalingFactor);
  mlMsh.ReadCoarseMesh ("./input/quadAMR.neu", "seventh", scalingFactor);
  /* "seventh" is the order of accuracy that is used in the gauss integration scheme
     probably in the furure it is not going to be an argument of this function   */
  unsigned dim = mlMsh.GetDimension();

  unsigned maxNumberOfMeshes = 5;


  vector < vector < double > > l2Norm;
  l2Norm.resize (maxNumberOfMeshes);

  vector < vector < double > > semiNorm;
  semiNorm.resize (maxNumberOfMeshes);

//   unsigned numberOfUniformLevels = 3;
//   unsigned numberOfSelectiveLevels = 0;
//   mlMsh.RefineMesh(numberOfUniformLevels , numberOfUniformLevels + numberOfSelectiveLevels, NULL);

  for (unsigned i = 1; i < maxNumberOfMeshes; i++) {

    unsigned numberOfUniformLevels = i + 3;
    unsigned numberOfSelectiveLevels = 0;
    //mlMsh.RefineMesh (numberOfUniformLevels + numberOfSelectiveLevels, numberOfUniformLevels , SetRefinementFlag);
    mlMsh.RefineMesh (numberOfUniformLevels + numberOfSelectiveLevels, numberOfUniformLevels , NULL);
    // erase all the coarse mesh levels
    //mlMsh.EraseCoarseLevels(numberOfUniformLevels - 3);

    // print mesh info
    mlMsh.PrintInfo();

    FEOrder feOrder[3] = {FIRST, SERENDIPITY, SECOND};
    l2Norm[i].resize (3);
    semiNorm[i].resize (3);

    for (unsigned j = 0; j < 3; j++) {

      MultiLevelSolution mlSol (&mlMsh);

      // add variables to mlSol

      mlSol.AddSolution("Flag",  DISCONTINOUS_POLYNOMIAL, ZERO);

      mlSol.AddSolution ("U", LAGRANGE, feOrder[j]);

      mlSol.Initialize ("All");

      // attach the boundary condition function and generate boundary data
      mlSol.AttachSetBoundaryConditionFunction (SetBoundaryCondition);

      mlSol.GenerateBdc ("All");

      // define the multilevel problem attach the mlSol object to it
      MultiLevelProblem mlProb (&mlSol);

      // add system Poisson in mlProb as a Linear Implicit System
      NonLinearImplicitSystem& system = mlProb.add_system < NonLinearImplicitSystem > ("Poisson");

      // add solution "u" to system
      system.AddSolutionToSystemPDE ("U");

      //system.SetMgSmoother(GMRES_SMOOTHER);
      system.SetMgSmoother (ASM_SMOOTHER); // Additive Swartz Method
      // attach the assembling function to system
      system.SetAssembleFunction (AssemblePoisson_AD);

      system.SetMaxNumberOfNonLinearIterations (10);
      system.SetMaxNumberOfLinearIterations (3);
      system.SetAbsoluteLinearConvergenceTolerance (1.e-12);
      system.SetNonLinearConvergenceTolerance (1.e-8);
      system.SetMgType (F_CYCLE);

      system.SetNumberPreSmoothingStep (0);
      system.SetNumberPostSmoothingStep (2);
      // initilaize and solve the system
      system.init();

      system.SetSolverFineGrids (GMRES);
      system.SetPreconditionerFineGrids (ILU_PRECOND);
      system.SetTolerances (1.e-3, 1.e-20, 1.e+50, 5);

      system.SetNumberOfSchurVariables (1);
      system.SetElementBlockNumber (4);
      //system.SetDirichletBCsHandling(ELIMINATION);
      //system.solve();
      system.MGsolve();

      std::pair< double , double > norm = GetErrorNorm (&mlSol);
      l2Norm[i][j]  = norm.first;
      semiNorm[i][j] = norm.second;
//.........这里部分代码省略.........

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

  unsigned precType = 0;

  if(argc >= 2) {
    if(!strcmp("FS_VT", args[1])) precType = FS_VTp;
    else if(!strcmp("FS_TV", args[1])) precType = FS_TVp;
    else if(!strcmp("ASM_VT", args[1])) precType = ASM_VTp;
    else if(!strcmp("ASM_TV", args[1])) precType = ASM_TVp;
    else if(!strcmp("ILU_VT", args[1])) precType = ILU_VTp;

    if(!strcmp("ILU_TV", args[1])) precType = ILU_TVp;

    if(precType == 0) {
      std::cout << "wrong input arguments!" << std::endl;
      abort();
    }
  }
  else {
    std::cout << "No input argument set default preconditioner = NS+T" << std::endl;
    precType = FS_VTp;
  }
  
  if(argc >= 3) {
    Prandtl = strtod(args[2], NULL);
    std::cout << Prandtl<<std::endl;
  }
  
  
  if(argc >= 4) {
    Rayleigh = strtod(args[3], NULL);
    std::cout << Rayleigh <<std::endl;
  }
    
  // init Petsc-MPI communicator
  FemusInit mpinit(argc, args, MPI_COMM_WORLD);

  // define multilevel mesh
  MultiLevelMesh mlMsh;
  // read coarse level mesh and generate finers level meshes
  double scalingFactor = 1.;
  //mlMsh.ReadCoarseMesh("./input/cube_hex.neu","seventh",scalingFactor);
  mlMsh.ReadCoarseMesh("./input/square_quad.neu", "seventh", scalingFactor);
  /* "seventh" is the order of accuracy that is used in the gauss integration scheme
     probably in the furure it is not going to be an argument of this function   */
  unsigned dim = mlMsh.GetDimension();

  unsigned numberOfUniformLevels = 8;
  unsigned numberOfSelectiveLevels = 0;
  mlMsh.RefineMesh(numberOfUniformLevels , numberOfUniformLevels + numberOfSelectiveLevels, NULL);
  // erase all the coarse mesh levels
  //mlMsh.EraseCoarseLevels(2);

  // print mesh info
  mlMsh.PrintInfo();

  MultiLevelSolution mlSol(&mlMsh);

  // add variables to mlSol
  mlSol.AddSolution("T", LAGRANGE, SECOND);
  mlSol.AddSolution("U", LAGRANGE, SECOND);
  mlSol.AddSolution("V", LAGRANGE, SECOND);

  if(dim == 3) mlSol.AddSolution("W", LAGRANGE, SECOND);

  mlSol.AddSolution("P",  DISCONTINOUS_POLYNOMIAL, FIRST);

  mlSol.AssociatePropertyToSolution("P", "Pressure");
  mlSol.Initialize("All");
  mlSol.Initialize("T",InitalValueT);

  // attach the boundary condition function and generate boundary data
  mlSol.AttachSetBoundaryConditionFunction(SetBoundaryCondition);
  mlSol.FixSolutionAtOnePoint("P");
  mlSol.GenerateBdc("All");

  // define the multilevel problem attach the mlSol object to it
  MultiLevelProblem mlProb(&mlSol);

  // add system Poisson in mlProb as a Linear Implicit System
  NonLinearImplicitSystem& system = mlProb.add_system < NonLinearImplicitSystem > ("NS");

  if(precType == FS_TVp || precType == ASM_TVp || precType == ILU_TVp)
    system.AddSolutionToSystemPDE("T");

  // add solution "u" to system
  system.AddSolutionToSystemPDE("U");
  system.AddSolutionToSystemPDE("V");

  if(precType == ASM_VTp)
    system.AddSolutionToSystemPDE("T");

  if(dim == 3) system.AddSolutionToSystemPDE("W");

  system.AddSolutionToSystemPDE("P");

  if(precType == FS_VTp || precType == ILU_VTp) system.AddSolutionToSystemPDE("T");


  //BEGIN buid fieldSplitTree (only for FieldSplitPreconditioner)
//.........这里部分代码省略.........

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

  unsigned precType = 0;

  if(argc >= 2) {
    Prandtl = strtod(args[1], NULL);
    std::cout << Prandtl << std::endl;
  }

  if(argc >= 3) {
    Rayleigh = strtod(args[2], NULL);
    std::cout << Rayleigh << std::endl;
  }

  std::cout << Prandtl<<" "<< Rayleigh << std::endl; 
  
  // init Petsc-MPI communicator
  FemusInit mpinit(argc, args, MPI_COMM_WORLD);
  // define multilevel mesh
  MultiLevelMesh mlMsh;
  // read coarse level mesh and generate finers level meshes
  double scalingFactor = 1.;
  //mlMsh.ReadCoarseMesh("./input/cube_hex.neu","seventh",scalingFactor);
  mlMsh.ReadCoarseMesh("./input/hex_cube.neu", "seventh", scalingFactor);
  /* "seventh" is the order of accuracy that is used in the gauss integration scheme
     probably in the furure it is not going to be an argument of this function   */
  unsigned dim = mlMsh.GetDimension();

  unsigned numberOfUniformLevels = 4;
  unsigned numberOfSelectiveLevels = 0;
 // mlMsh.RefineMesh(numberOfUniformLevels , numberOfUniformLevels + numberOfSelectiveLevels, NULL);
  mlMsh.RefineMesh(numberOfUniformLevels + numberOfSelectiveLevels, numberOfUniformLevels, SetRefinementFlag);
  // erase all the coarse mesh levels
  //mlMsh.EraseCoarseLevels(0);

  // print mesh info
  mlMsh.PrintInfo();
  MultiLevelSolution mlSol(&mlMsh);

  // add variables to mlSol
  mlSol.AddSolution("T", LAGRANGE, FIRST);
  mlSol.AddSolution("U", LAGRANGE, SECOND);
  mlSol.AddSolution("V", LAGRANGE, SECOND);

  if(dim == 3) mlSol.AddSolution("W", LAGRANGE, SECOND);

  mlSol.AddSolution("P",  DISCONTINUOUS_POLYNOMIAL, FIRST);
  mlSol.AssociatePropertyToSolution("P", "Pressure");
  mlSol.Initialize("All");
  mlSol.Initialize("T", InitalValueT);

  // attach the boundary condition function and generate boundary data
  mlSol.AttachSetBoundaryConditionFunction(SetBoundaryCondition);
  mlSol.FixSolutionAtOnePoint("P");
  mlSol.GenerateBdc("All");

  // define the multilevel problem attach the mlSol object to it
  MultiLevelProblem mlProb(&mlSol);

  // add system Poisson in mlProb as a Linear Implicit System
  NonLinearImplicitSystem& system = mlProb.add_system < NonLinearImplicitSystem > ("NS");
  // add solution "u" to system
  system.AddSolutionToSystemPDE("U");
  system.AddSolutionToSystemPDE("V");
  if(dim == 3) system.AddSolutionToSystemPDE("W");
  system.AddSolutionToSystemPDE("P");
  system.AddSolutionToSystemPDE("T");
  system.SetLinearEquationSolverType(FEMuS_DEFAULT);

  // attach the assembling function to system
  system.SetAssembleFunction(AssembleBoussinesqAppoximation);

  system.SetMaxNumberOfNonLinearIterations(20);
  system.SetNonLinearConvergenceTolerance(1.e-8);
  //system.SetMaxNumberOfResidualUpdatesForNonlinearIteration(10);
  //system.SetResidualUpdateConvergenceTolerance(1.e-15);

  //system.SetMaxNumberOfLinearIterations(10);
  //system.SetAbsoluteLinearConvergenceTolerance(1.e-15);

  system.SetMaxNumberOfLinearIterations(1);
  system.SetAbsoluteLinearConvergenceTolerance(1.e-15);

  system.SetMgType(V_CYCLE);
  system.SetNumberPreSmoothingStep(1);
  system.SetNumberPostSmoothingStep(1);
  // initilaize and solve the system
  system.init();

  system.SetSolverFineGrids(RICHARDSON);
  system.SetRichardsonScaleFactor(.6);
  system.SetPreconditionerFineGrids(ILU_PRECOND);

  system.SetTolerances(1.e-8, 1.e-15, 1.e+50, 30, 30); //GMRES tolerances
  system.ClearVariablesToBeSolved();
  system.AddVariableToBeSolved("All");
  system.SetNumberOfSchurVariables(1);
  system.SetElementBlockNumber("All");

  system.MGsolve();
//.........这里部分代码省略.........

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

  // init Petsc-MPI communicator
  FemusInit mpinit(argc, args, MPI_COMM_WORLD);

  // define multilevel mesh
  MultiLevelMesh mlMsh;
  // read coarse level mesh and generate finers level meshes
  double scalingFactor = 1.;
  //mlMsh.ReadCoarseMesh("./input/cube_hex.neu","seventh",scalingFactor);
  mlMsh.ReadCoarseMesh("./input/square_quad.neu", "seventh", scalingFactor);
  /* "seventh" is the order of accuracy that is used in the gauss integration scheme
     probably in the furure it is not going to be an argument of this function   */
  unsigned dim = mlMsh.GetDimension();

  unsigned numberOfUniformLevels = 7;
  unsigned numberOfSelectiveLevels = 0;
  mlMsh.RefineMesh(numberOfUniformLevels , numberOfUniformLevels + numberOfSelectiveLevels, NULL);

  // erase all the coarse mesh levels
  //mlMsh.EraseCoarseLevels(numberOfUniformLevels - 3);

  // print mesh info
  mlMsh.PrintInfo();

  MultiLevelSolution mlSol(&mlMsh);

  // add variables to mlSol
  mlSol.AddSolution("T", LAGRANGE, FIRST);
  mlSol.AddSolution("U", LAGRANGE, SECOND);
  mlSol.AddSolution("V", LAGRANGE, SECOND);

  if (dim == 3) mlSol.AddSolution("W", LAGRANGE, SECOND);

  //mlSol.AddSolution("P", LAGRANGE, FIRST);
  mlSol.AddSolution("P",  DISCONTINOUS_POLYNOMIAL, FIRST);

  mlSol.AssociatePropertyToSolution("P", "Pressure");
  mlSol.Initialize("All");

  // attach the boundary condition function and generate boundary data
  mlSol.AttachSetBoundaryConditionFunction(SetBoundaryCondition);
  mlSol.FixSolutionAtOnePoint("P");
  mlSol.GenerateBdc("All");

  // define the multilevel problem attach the mlSol object to it
  MultiLevelProblem mlProb(&mlSol);

  // add system Poisson in mlProb as a Linear Implicit System
  NonLinearImplicitSystem& system = mlProb.add_system < NonLinearImplicitSystem > ("NS");

  // add solution "u" to system
  system.AddSolutionToSystemPDE("U");
  system.AddSolutionToSystemPDE("V");
  system.AddSolutionToSystemPDE("P");

  if (dim == 3) system.AddSolutionToSystemPDE("W");

  system.AddSolutionToSystemPDE("T");

  //system.SetMgSmoother(GMRES_SMOOTHER);
  system.SetMgSmoother(FIELDSPLIT_SMOOTHER); // Additive Swartz Method
  //system.SetMgSmoother(ASM_SMOOTHER); // Additive Swartz Method
  // attach the assembling function to system
  system.SetAssembleFunction(AssembleBoussinesqAppoximation_AD);

  system.SetMaxNumberOfNonLinearIterations(20);
  system.SetMaxNumberOfLinearIterations(3);
  system.SetLinearConvergenceTolerance(1.e-12);
  system.SetNonLinearConvergenceTolerance(1.e-8);
  system.SetMgType(F_CYCLE);

  system.SetNumberPreSmoothingStep(0);
  system.SetNumberPostSmoothingStep(2);
  // initilaize and solve the system
  system.init();

  system.SetSolverFineGrids(GMRES);
  system.SetPreconditionerFineGrids(ILU_PRECOND);
  //system.SetTolerances(1.e-20, 1.e-20, 1.e+50, 40);
  system.SetTolerances(1.e-3, 1.e-20, 1.e+50, 5);


  system.ClearVariablesToBeSolved();
  system.AddVariableToBeSolved("All");
  system.SetNumberOfSchurVariables(1);
  system.SetElementBlockNumber(4);
  //system.SetDirichletBCsHandling(ELIMINATION);
  //system.solve();
  system.MGsolve();

  // print solutions
  std::vector < std::string > variablesToBePrinted;
  variablesToBePrinted.push_back("All");

  VTKWriter vtkIO(&mlSol);
  vtkIO.write(DEFAULT_OUTPUTDIR, "biquadratic", variablesToBePrinted);

  return 0;
}

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

    // init Petsc-MPI communicator
    FemusInit mpinit(argc, args, MPI_COMM_WORLD);


    // define multilevel mesh


    unsigned maxNumberOfMeshes;
    maxNumberOfMeshes = 1;

    vector < vector < double > > l2Norm;
    l2Norm.resize(maxNumberOfMeshes);

    vector < vector < double > > semiNorm;
    semiNorm.resize(maxNumberOfMeshes);

    for (unsigned i = 0; i < maxNumberOfMeshes; i++) {   // loop on the mesh level

        std::ostringstream filename;

        filename << "./input/sphere.neu";

        MultiLevelMesh mlMsh;
        // read coarse level mesh and generate finers level meshes
        double scalingFactor = 2.;
        //mlMsh.ReadCoarseMesh("./input/circle_quad.neu","seventh", scalingFactor);
        mlMsh.ReadCoarseMesh(filename.str().c_str(), "seventh", scalingFactor);
        /* "seventh" is the order of accuracy that is used in the gauss integration scheme
        probably in the furure it is not going to be an argument of this function   */
        unsigned dim = mlMsh.GetDimension();

        unsigned numberOfUniformLevels = 4;
        unsigned numberOfSelectiveLevels = 0;
        mlMsh.RefineMesh(numberOfUniformLevels , numberOfUniformLevels + numberOfSelectiveLevels, NULL);

        // erase all the coarse mesh levels
        mlMsh.EraseCoarseLevels(numberOfUniformLevels - 1);

        // print mesh info
        mlMsh.PrintInfo();

        FEOrder feOrder = SECOND;
        l2Norm[i].resize(1);
        semiNorm[i].resize(1);

        // define the multilevel solution and attach the mlMsh object to it
        MultiLevelSolution mlSol(&mlMsh);

        // add variables to mlSol
        mlSol.AddSolution("X", LAGRANGE, feOrder, 2);
        mlSol.AddSolution("Y", LAGRANGE, feOrder, 2);
        mlSol.AddSolution("Z", LAGRANGE, feOrder, 2);
        mlSol.AddSolution("H", LAGRANGE, feOrder, 2);

        mlSol.Initialize("X", InitalValueXTorus);
        mlSol.Initialize("Y", InitalValueYTorus);
        mlSol.Initialize("Z", InitalValueZTorus);
        mlSol.Initialize("H", InitalValueHTorus);
        // attach the boundary condition function and generate boundary data
        mlSol.AttachSetBoundaryConditionFunction(SetBoundaryConditionTorus);

//     mlSol.Initialize("X", InitalValueXSphere);
//     mlSol.Initialize("Y", InitalValueYSphere);
//     mlSol.Initialize("Z", InitalValueZSphere);
//     mlSol.Initialize("H", InitalValueHSphere);
//         // attach the boundary condition function and generate boundary data
//     mlSol.AttachSetBoundaryConditionFunction(SetBoundaryConditionSphere);


        mlSol.GenerateBdc("X", "Steady");
        mlSol.GenerateBdc("Y", "Steady");
        mlSol.GenerateBdc("Z", "Steady");
        mlSol.GenerateBdc("H", "Steady");

        // define the multilevel problem attach the mlSol object to it
        MultiLevelProblem mlProb(&mlSol);

        // add system Wilmore in mlProb as a Linear Implicit System
        TransientNonlinearImplicitSystem& system = mlProb.add_system < TransientNonlinearImplicitSystem > ("Willmore");

        // add solution "X", "Y", "Z" and "H" to the system
        system.AddSolutionToSystemPDE("X");
        system.AddSolutionToSystemPDE("Y");
        system.AddSolutionToSystemPDE("Z");
        system.AddSolutionToSystemPDE("H");

        system.SetMaxNumberOfNonLinearIterations(6);

        // attach the assembling function to system
        system.SetAssembleFunction(AssembleWillmoreFlow_AD);

        // initilaize and solve the system
        system.init();

        std::vector < std::string > variablesToBePrinted;
        variablesToBePrinted.push_back("All");

        std::vector < std::string > surfaceVariables;
//.........这里部分代码省略.........

//.........这里部分代码省略.........
    system.init();

    system.SetSolverFineGrids(GMRES);
    system.SetPreconditionerFineGrids(ILU_PRECOND);
    system.SetTolerances(1.e-3, 1.e-20, 1.e+50, 5);

    system.SetNumberOfSchurVariables(1);
    system.SetElementBlockNumber(4);
   
    system.MGsolve();
    
    GetError(&mlSol);
    
    std::pair< double , double > norm = GetError (&mlSol);
      H1normE[i][0]  = norm.first;
      H1norm[i][0] = norm.second;
    
    // print solutions
    std::vector < std::string > variablesToBePrinted;
    variablesToBePrinted.push_back("All");
    VTKWriter vtkIO(&mlSol);
    vtkIO.SetDebugOutput(true);
    vtkIO.Write(DEFAULT_OUTPUTDIR, "biquadratic", variablesToBePrinted);
    
    //refine the mesh
    MeshRefinement meshcoarser(*mlMsh.GetLevel(numberOfUniformLevels-1));
    bool elementsHaveBeenRefined = meshcoarser.FlagElementsToBeRefined(0.005, mlSol.GetSolutionLevel(numberOfUniformLevels-1)->GetSolutionName("Error"));  //non-uniform
    
    //bool elementsHaveBeenRefined = true; //uniform
    //meshcoarser.FlagAllElementsToBeRefined();//uniform
    if( !elementsHaveBeenRefined ){
      std::cout << " the solution has converged\n"; 
      maxNumberOfMeshes = i + 1;
      break;
    }
    mlMsh.AddAMRMeshLevel();
    mlSol.AddSolutionLevel();
    mlSol.RefineSolution(numberOfUniformLevels);
    //}
    numberOfUniformLevels += 1;
    
    

  }
  
  
  // print the seminorm of the error and the order of convergence between different levels
  std::cout << std::endl;
  std::cout << std::endl;
  std::cout << "H1 ERROR and ORDER OF CONVERGENCE:\n\n";
  std::cout << "LEVEL \n";

  for (unsigned i = 0; i < maxNumberOfMeshes; i++) {
    std::cout << i + 1 << "\t";
    std::cout.precision (14);

    //for (unsigned j = 0; j < 3; j++) {
      std::cout << H1normE[i][0] << "\t";
    //}

    std::cout << std::endl;

    if (i < maxNumberOfMeshes - 1) {
      std::cout.precision(3);
      std::cout << "\t\t";

      //for (unsigned j = 0; j < 3; j++) {
        std::cout << log(H1normE[i][0] / H1normE[i + 1][0]) / log(2.) << "\t\t\t";
      //}

      std::cout << std::endl;
    }     


  }

  std::cout << std::endl;
  std::cout << std::endl;
  std::cout << "H1 Relative ERROR and ORDER OF CONVERGENCE:\n\n";
  std::cout << "LEVEL \n";

  for (unsigned i = 0; i < maxNumberOfMeshes; i++) {
    std::cout << i + 1 << "\t";
    std::cout.precision (14);

    //for (unsigned j = 0; j < 3; j++) {
      std::cout << H1normE[i][0] / H1norm[i][0] << "\t";
    //}

    std::cout << std::endl;



  }


   mlMsh.PrintInfo();
  
  return 0;
}

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

  // init Petsc-MPI communicator
  FemusInit mpinit(argc, args, MPI_COMM_WORLD);
  // define multilevel mesh
  MultiLevelMesh mlMsh;

  // read coarse level mesh and generate finers level meshes
  double scalingFactor = 1.;
  //mlMsh.ReadCoarseMesh("./input/cube_hex.neu","seventh",scalingFactor);
  mlMsh.ReadCoarseMesh("./input/rectangle_w1_h8.neu", "seventh", scalingFactor);
  /* "seventh" is the order of accuracy that is used in the gauss integration scheme
     probably in the furure it is not going to be an argument of this function   */
  unsigned dim = mlMsh.GetDimension();

  unsigned numberOfUniformLevels = 4;
  unsigned numberOfSelectiveLevels = 0;
  mlMsh.RefineMesh(numberOfUniformLevels , numberOfUniformLevels + numberOfSelectiveLevels, NULL);

  // erase all the coarse mesh levels
  //mlMsh.EraseCoarseLevels(numberOfUniformLevels - 3);

  // print mesh info
  mlMsh.PrintInfo();
  MultiLevelSolution mlSol(&mlMsh);

  // add variables to mlSol
  mlSol.AddSolution("T", LAGRANGE, SERENDIPITY, 2);
  mlSol.AddSolution("U", LAGRANGE, SECOND, 2);
  mlSol.AddSolution("V", LAGRANGE, SECOND, 2);

  if(dim == 3) mlSol.AddSolution("W", LAGRANGE, SECOND, 2);

  //mlSol.AddSolution("P", LAGRANGE, FIRST);
  mlSol.AddSolution("P",  DISCONTINOUS_POLYNOMIAL, FIRST, 2);

  mlSol.AssociatePropertyToSolution("P", "Pressure");
  mlSol.Initialize("All");

  mlSol.Initialize("T", InitalValueT);

  // attach the boundary condition function and generate boundary data
  mlSol.AttachSetBoundaryConditionFunction(SetBoundaryCondition);
  mlSol.FixSolutionAtOnePoint("P");
  mlSol.GenerateBdc("U");
  mlSol.GenerateBdc("V");
  mlSol.GenerateBdc("P");
  mlSol.GenerateBdc("T", "Time_dependent");

  // define the multilevel problem attach the mlSol object to it
  MultiLevelProblem mlProb(&mlSol);

  // add system Poisson in mlProb as a Linear Implicit System
  TransientNonlinearImplicitSystem& system = mlProb.add_system < TransientNonlinearImplicitSystem > ("NS");

  // add solution "u" to system
  system.AddSolutionToSystemPDE("U");
  system.AddSolutionToSystemPDE("V");
  system.AddSolutionToSystemPDE("P");

  if(dim == 3) system.AddSolutionToSystemPDE("W");

  system.AddSolutionToSystemPDE("T");


  std::vector < unsigned > fieldUVP(3);
  fieldUVP[0] = system.GetSolPdeIndex("U");
  fieldUVP[1] = system.GetSolPdeIndex("V");
  fieldUVP[2] = system.GetSolPdeIndex("P");

  std::vector < unsigned > solutionTypeUVP(3);
  solutionTypeUVP[0] = mlSol.GetSolutionType("U");
  solutionTypeUVP[1] = mlSol.GetSolutionType("V");
  solutionTypeUVP[2] = mlSol.GetSolutionType("P");

  FieldSplitTree FS_NS(PREONLY, ASM_PRECOND, fieldUVP, solutionTypeUVP, "Navier-Stokes");
  FS_NS.SetAsmBlockSize(4);
  FS_NS.SetAsmNumeberOfSchurVariables(1);


//   std::vector < unsigned > fieldUV(2);
//   fieldUV[0] = system.GetSolPdeIndex("U");
//   fieldUV[1] = system.GetSolPdeIndex("V");
//
//   FieldSplitTree FS_UV( PREONLY, ILU_PRECOND, fieldUV, "Velocity");
//
//   std::vector < unsigned > fieldP(1);
//   fieldP[0] = system.GetSolPdeIndex("P");
//
//   FieldSplitTree FS_P( PREONLY, ILU_PRECOND, fieldP, "pressure");
//
//    std::vector < FieldSplitTree *> FS1;
//
//   FS1.reserve(2);
//   FS1.push_back(&FS_UV);
//   FS1.push_back(&FS_P);
//   FieldSplitTree FS_NS( GMRES, FS_SCHUR_PRECOND, FS1, "Navier-Stokes");



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

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

  // init Petsc-MPI communicator
  FemusInit mpinit(argc, args, MPI_COMM_WORLD);

  // define multilevel mesh
  MultiLevelMesh mlMsh;
  // read coarse level mesh and generate finers level meshes
  double scalingFactor = 1.;
  //mlMsh.ReadCoarseMesh("./input/cube_hex.neu","seventh",scalingFactor);
  mlMsh.ReadCoarseMesh("./input/square_quad.neu", "seventh", scalingFactor);
  /* "seventh" is the order of accuracy that is used in the gauss integration scheme
     probably in the furure it is not going to be an argument of this function   */
  unsigned dim = mlMsh.GetDimension();

  unsigned numberOfUniformLevels = 7;
  unsigned numberOfSelectiveLevels = 0;
  mlMsh.RefineMesh(numberOfUniformLevels , numberOfUniformLevels + numberOfSelectiveLevels, NULL);

  // erase all the coarse mesh levels
  //mlMsh.EraseCoarseLevels(numberOfUniformLevels - 1);

  // print mesh info
  mlMsh.PrintInfo();

  MultiLevelSolution mlSol(&mlMsh);

  // add variables to mlSol
  mlSol.AddSolution("U", LAGRANGE, SECOND);
  mlSol.AddSolution("V", LAGRANGE, SECOND);

  if (dim == 3) mlSol.AddSolution("W", LAGRANGE, SECOND);

  //mlSol.AddSolution("P", LAGRANGE, FIRST);
  mlSol.AddSolution("P",  DISCONTINOUS_POLYNOMIAL, FIRST);

  mlSol.AssociatePropertyToSolution("P", "Pressure");
  mlSol.Initialize("All");

  // attach the boundary condition function and generate boundary data
  mlSol.AttachSetBoundaryConditionFunction(SetBoundaryCondition);
  mlSol.FixSolutionAtOnePoint("P");
  mlSol.GenerateBdc("All");

  // define the multilevel problem attach the mlSol object to it
  MultiLevelProblem mlProb(&mlSol);

  // add system Poisson in mlProb as a Linear Implicit System
  NonLinearImplicitSystem& system = mlProb.add_system < NonLinearImplicitSystem > ("NS");

  // add solution "u" to system
  system.AddSolutionToSystemPDE("U");
  system.AddSolutionToSystemPDE("V");
  system.AddSolutionToSystemPDE("P");

  if (dim == 3) system.AddSolutionToSystemPDE("W");
 
  std::vector < unsigned > fieldUV(2);
  fieldUV[0] = system.GetSolPdeIndex("U");
  fieldUV[1] = system.GetSolPdeIndex("V");
  FieldSplitTree FS_UV( PREONLY, ILU_PRECOND, fieldUV , "Velocity");
  FS_UV.SetupKSPTolerances(1.e-3,1.e-20,1.e+50, 1); // changed by Guoyi Ke 

  std::vector < unsigned > fieldP(1);
  fieldP[0] = system.GetSolPdeIndex("P");

  //FS_P.SetFieldSplitSchurFactType{PC_FIELDSPLIT_SCHUR_FACT_LOWER}; //changed by Guoyi Ke
  FieldSplitTree FS_P(PREONLY, ILU_PRECOND, fieldP, "Pressure");
  FS_P.SetupKSPTolerances(1.e-3,1.e-20,1.e+50, 1); //changed by Guoyi Ke
  
  std::vector < FieldSplitTree *> FS1;
  FS1.reserve(2);
  FS1.push_back(&FS_UV);
  FS1.push_back(&FS_P);
  
  FieldSplitTree FS_NS(GMRES, FS_SCHUR_PRECOND, FS1, "Navier-Stokes");
  FS_NS.SetupSchurFactorizationType(SCHUR_FACT_UPPER); // SCHUR_FACT_UPPER, SCHUR_FACT_LOWER,SCHUR_FACT_FULL; how to use if FS_SCHUR_PRECOND? Guoyike
  FS_NS.SetupSchurPreType(SCHUR_PRE_SELFP);// SCHUR_PRE_SELF, SCHUR_PRE_SELFP, SCHUR_PRE_USER, SCHUR_PRE_A11,SCHUR_PRE_FULL;

  //system.SetMgSmoother(GMRES_SMOOTHER);
  system.SetMgSmoother(FIELDSPLIT_SMOOTHER); // Additive Swartz Method

  //system.SetMgSmoother(ASM_SMOOTHER); // Additive Swartz Method
  // attach the assembling function to system
  system.SetAssembleFunction(AssembleBoussinesqAppoximation_AD);

  system.SetMaxNumberOfNonLinearIterations(20);
  system.SetMaxNumberOfLinearIterations(3);
  system.SetAbsoluteLinearConvergenceTolerance(1.e-12);
  system.SetNonLinearConvergenceTolerance(1.e-8);
  system.SetMgType(F_CYCLE);

  system.SetNumberPreSmoothingStep(2);
  system.SetNumberPostSmoothingStep(2);

  // initilaize and solve the system
  system.init();

  system.SetSolverFineGrids(GMRES);
  system.SetPreconditionerFineGrids(ILU_PRECOND);
//.........这里部分代码省略.........

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

  unsigned precType = 0;

  if (argc >= 2) {
    Miu = strtod(args[1], NULL);
    std::cout << Miu << std::endl;
  }

  // init Petsc-MPI communicator
  FemusInit mpinit(argc, args, MPI_COMM_WORLD);

  // define multilevel mesh
  MultiLevelMesh mlMsh;
  // read coarse level mesh and generate finers level meshes
  double scalingFactor = 1.;
  //mlMsh.ReadCoarseMesh("./input/cube_hex.neu","seventh",scalingFactor);
  mlMsh.ReadCoarseMesh("./input/quad_square.neu", "seventh", scalingFactor);
  /* "seventh" is the order of accuracy that is used in the gauss integration scheme
     probably in the furure it is not going to be an argument of this function   */
  unsigned dim = mlMsh.GetDimension();

  unsigned numberOfUniformLevels = 8;
  unsigned numberOfSelectiveLevels = 0;
  mlMsh.RefineMesh(numberOfUniformLevels + numberOfSelectiveLevels, numberOfUniformLevels, SetRefinementFlag);
 // mlMsh.RefineMesh(numberOfUniformLevels , numberOfUniformLevels + numberOfSelectiveLevels, NULL);
  // erase all the coarse mesh levels
  mlMsh.EraseCoarseLevels(3);

  // print mesh info
  mlMsh.PrintInfo();
  MultiLevelSolution mlSol(&mlMsh);

  // add variables to mlSol
  mlSol.AddSolution("U", LAGRANGE, SECOND);
  mlSol.AddSolution("V", LAGRANGE, SECOND);
  if (dim == 3) mlSol.AddSolution("W", LAGRANGE, SECOND);
  mlSol.AddSolution("P",  DISCONTINUOUS_POLYNOMIAL, FIRST);
  mlSol.AssociatePropertyToSolution("P", "Pressure");
  mlSol.Initialize("All");
//   mlSol.Initialize("U", InitalValueU);

  // attach the boundary condition function and generate boundary data
  mlSol.AttachSetBoundaryConditionFunction(SetBoundaryCondition);
  mlSol.FixSolutionAtOnePoint("P");
  mlSol.GenerateBdc("All");

  // define the multilevel problem attach the mlSol object to it
  MultiLevelProblem mlProb(&mlSol);

  // add system Poisson in mlProb as a Linear Implicit System
  NonLinearImplicitSystem& system = mlProb.add_system < NonLinearImplicitSystem > ("NS");

  // add solution "u" to system
  system.AddSolutionToSystemPDE("U");
  system.AddSolutionToSystemPDE("V");
  if (dim == 3) system.AddSolutionToSystemPDE("W");
  system.AddSolutionToSystemPDE("P");
  //system.SetLinearEquationSolverType(FEMuS_DEFAULT);
  system.SetLinearEquationSolverType(FEMuS_ASM);
  
  // attach the assembling function to system
  system.SetAssembleFunction(AssembleBoussinesqAppoximation);

  system.SetMaxNumberOfNonLinearIterations(20);
  system.SetNonLinearConvergenceTolerance(1.e-8);
  //system.SetMaxNumberOfResidualUpdatesForNonlinearIteration(10);
  //system.SetResidualUpdateConvergenceTolerance(1.e-15);

  system.SetMaxNumberOfLinearIterations(1);
  system.SetAbsoluteLinearConvergenceTolerance(1.e-15);


  system.SetMgType(V_CYCLE);
  system.SetNumberPreSmoothingStep(1);
  system.SetNumberPostSmoothingStep(1);
  // initilaize and solve the system
  system.init();

  system.SetSolverFineGrids(RICHARDSON);
  system.SetPreconditionerFineGrids(MLU_PRECOND);

  system.SetTolerances(1.e-5, 1.e-8, 1.e+50, 30, 30); //GMRES tolerances

  system.ClearVariablesToBeSolved();
  system.AddVariableToBeSolved("All");
  system.SetNumberOfSchurVariables(1);
  system.SetElementBlockNumber("All");

  system.MGsolve();
  // print solutions
  std::vector < std::string > variablesToBePrinted;
  variablesToBePrinted.push_back("All");

  VTKWriter vtkIO(&mlSol);
  vtkIO.Write(DEFAULT_OUTPUTDIR, "biquadratic", variablesToBePrinted);

  mlMsh.PrintInfo();
  
//.........这里部分代码省略.........

int main(int argc, char** args) {
  
  unsigned precType = 0;

  if(argc >= 2) {
    if(!strcmp("FS_VT", args[1])) precType = FS_VTp;
    else if(!strcmp("FS_TV", args[1])) precType = FS_TVp;
    else if(!strcmp("ASM_VT", args[1])) precType = ASM_VTp;
    else if(!strcmp("ASM_TV", args[1])) precType = ASM_TVp;
    else if(!strcmp("ILU_VT", args[1])) precType = ILU_VTp;

    if(!strcmp("ILU_TV", args[1])) precType = ILU_TVp;

    if(precType == 0) {
      std::cout << "wrong input arguments!" << std::endl;
      abort();
    }
  }
  else {
    std::cout << "No input argument set default preconditioner = NS+T" << std::endl;
    precType = FS_VTp;
  }

  if(argc >= 3) {
    Prandtl = strtod(args[2], NULL);
    std::cout << Prandtl << std::endl;
  }


  if(argc >= 4) {
    Rayleigh = strtod(args[3], NULL);
    std::cout << Rayleigh << std::endl;
  }
  
  
  // init Petsc-MPI communicator
  FemusInit mpinit(argc, args, MPI_COMM_WORLD);
  // define multilevel mesh
  MultiLevelMesh mlMsh;

  // read coarse level mesh and generate finers level meshes
  double scalingFactor = 1.;
  //mlMsh.ReadCoarseMesh("./input/cube_hex.neu","seventh",scalingFactor);
  mlMsh.ReadCoarseMesh("./input/rectangle_w4_h1.neu", "seventh", scalingFactor);
  /* "seventh" is the order of accuracy that is used in the gauss integration scheme
     probably in the furure it is not going to be an argument of this function   */
  unsigned dim = mlMsh.GetDimension();

  unsigned numberOfUniformLevels = 7;
  unsigned numberOfSelectiveLevels = 0;
  mlMsh.RefineMesh(numberOfUniformLevels , numberOfUniformLevels + numberOfSelectiveLevels, NULL);

  // erase all the coarse mesh levels
  //mlMsh.EraseCoarseLevels(numberOfUniformLevels - 3);
  mlMsh.EraseCoarseLevels(3);

  // print mesh info
  mlMsh.PrintInfo();
  MultiLevelSolution mlSol(&mlMsh);

  // add variables to mlSol
  mlSol.AddSolution("T", LAGRANGE, SERENDIPITY, 2);
  mlSol.AddSolution("U", LAGRANGE, SECOND, 2);
  mlSol.AddSolution("V", LAGRANGE, SECOND, 2);

  if(dim == 3) mlSol.AddSolution("W", LAGRANGE, SECOND, 2);

  //mlSol.AddSolution("P", LAGRANGE, FIRST);
  mlSol.AddSolution("P",  DISCONTINUOUS_POLYNOMIAL, FIRST, 2);

  mlSol.AssociatePropertyToSolution("P", "Pressure");
  mlSol.Initialize("All");
  mlSol.Initialize("T", InitalValueT);

  // attach the boundary condition function and generate boundary data
  mlSol.AttachSetBoundaryConditionFunction(SetBoundaryCondition);
  mlSol.FixSolutionAtOnePoint("P");
  mlSol.GenerateBdc("U");
  mlSol.GenerateBdc("V");
  mlSol.GenerateBdc("P");
  mlSol.GenerateBdc("T", "Time_dependent");

  // define the multilevel problem attach the mlSol object to it
  MultiLevelProblem mlProb(&mlSol);

  // add system Poisson in mlProb as a Linear Implicit System
  TransientNonlinearImplicitSystem& system = mlProb.add_system < TransientNonlinearImplicitSystem > ("NS");

  
  if(precType == FS_TVp || precType == ASM_TVp || precType == ILU_TVp) system.AddSolutionToSystemPDE("T");
  // add solution "u" to system
  system.AddSolutionToSystemPDE("U");
  system.AddSolutionToSystemPDE("V");
  
  if(precType == ASM_VTp) system.AddSolutionToSystemPDE("T");
  
  if(dim == 3) system.AddSolutionToSystemPDE("W");
  
  system.AddSolutionToSystemPDE("P");
  
//.........这里部分代码省略.........