C++ Files::RedirectCout方法代码示例

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

    // ======= Init ========================
  // init Petsc-MPI communicator
  FemusInit mpinit(argc, args, MPI_COMM_WORLD);
  
    // ======= Files ========================
  Files files; 
        files.CheckIODirectories();
        files.RedirectCout();

    // ======= Quad Rule ========================
  std::string fe_quad_rule("seventh");
 /* "seventh" is the order of accuracy that is used in the gauss integration scheme
    In the future it is not going to be an argument of the mesh function   */
  
    // ======= Mesh ========================
  MultiLevelMesh ml_mesh;
  ml_mesh.GenerateCoarseBoxMesh(NSUB_X,NSUB_Y,0,0.,1.,0.,1.,0.,0.,QUAD9,fe_quad_rule.c_str());
  unsigned numberOfUniformLevels = 1;
  unsigned numberOfSelectiveLevels = 0;
  ml_mesh.RefineMesh(numberOfUniformLevels , numberOfUniformLevels + numberOfSelectiveLevels, NULL);
  ml_mesh.PrintInfo();

    // ======= Solution ========================
  MultiLevelSolution ml_sol(&ml_mesh);  // define the multilevel solution and attach the ml_mesh object to it

  // add variables to ml_sol
  ml_sol.AddSolution("state",   LAGRANGE, FIRST);
  ml_sol.AddSolution("control", LAGRANGE, FIRST);
  ml_sol.AddSolution("adjoint", LAGRANGE, FIRST);
  ml_sol.AddSolution("mu",      LAGRANGE, FIRST);  
  ml_sol.AddSolution("TargReg", DISCONTINUOUS_POLYNOMIAL, ZERO); //this variable is not solution of any eqn, it's just a given field
  ml_sol.AddSolution("ContReg", DISCONTINUOUS_POLYNOMIAL, ZERO); //this variable is not solution of any eqn, it's just a given field

  const unsigned int fake_time_dep_flag = 2;  //this is needed to be able to use _SolOld
  const std::string act_set_flag_name = "act_flag";
  ml_sol.AddSolution(act_set_flag_name.c_str(), LAGRANGE, FIRST,fake_time_dep_flag);               

    // ======= Problem ========================
  MultiLevelProblem ml_prob(&ml_sol);  // define the multilevel problem attach the ml_sol object to it

  ml_prob.SetQuadratureRuleAllGeomElems(fe_quad_rule);
  ml_prob.SetFilesHandler(&files);
  
    // ======= Initial values ========================
  ml_sol.Initialize("All");    // initialize all variables to zero

//   ml_sol.Initialize("All", SetInitialCondition, &ml_prob); //unfortunately if I do this it sets all to zero //I would like to do an attach function similar to the BC
  ml_sol.Initialize("state",   SetInitialCondition, &ml_prob);
  ml_sol.Initialize("control", SetInitialCondition, &ml_prob);
  ml_sol.Initialize("adjoint", SetInitialCondition, &ml_prob);
  ml_sol.Initialize("mu",      SetInitialCondition, &ml_prob);
  ml_sol.Initialize("TargReg", SetInitialCondition, &ml_prob);
  ml_sol.Initialize("ContReg", SetInitialCondition, &ml_prob);
  ml_sol.Initialize(act_set_flag_name.c_str(),  SetInitialCondition, &ml_prob);

    // ======= Boundary Conditions ========================
  ml_sol.AttachSetBoundaryConditionFunction(SetBoundaryCondition);  // attach the boundary condition function and generate boundary data

//   ml_sol.GenerateBdc("All");  //this would do it also for the non-equation-related variables
  ml_sol.GenerateBdc("state");
  ml_sol.GenerateBdc("control");
  ml_sol.GenerateBdc("adjoint");
  ml_sol.GenerateBdc("mu");  //we need this for all Pde variables to make the matrix iterations work... but this should be related to the matrix and not to the sol... The same for the initial condition

    // ======= System ========================
  NonLinearImplicitSystemWithPrimalDualActiveSetMethod& system = ml_prob.add_system < NonLinearImplicitSystemWithPrimalDualActiveSetMethod > ("OptSys");

  system.SetActiveSetFlagName(act_set_flag_name);

  //here we decide the order in the matrix!
  const std::vector < std::string > sol_matrix_pos = {"state","control","adjoint","mu"};
  for (unsigned k = 0; k < sol_matrix_pos.size(); k++)  system.AddSolutionToSystemPDE(sol_matrix_pos[k].c_str());  
  
  // attach the assembling function to system
  system.SetAssembleFunction(AssembleProblem);
  
  ml_sol.SetWriter(VTK);   //need to move this here for the DebugNonlinear function
  ml_sol.GetWriter()->SetDebugOutput(true);
  
  system.SetDebugNonlinear(true);
  system.SetDebugFunction(ComputeIntegral);
  //   system.SetMaxNumberOfNonLinearIterations(2);

    // ======= Solve ========================
  system.init();    // initialize and solve the system
  system.MGsolve();
  
//   ComputeIntegral(ml_prob);
 
    // ======= Final Print ========================
  std::vector < std::string > variablesToBePrinted;
  variablesToBePrinted.push_back("all");
  ml_sol.GetWriter()->Write(files.GetOutputPath()/*DEFAULT_OUTPUTDIR*/, "biquadratic", variablesToBePrinted);    // print solutions


  return 0;
}

int main(int argc,char **args) {
  
  /// Init Petsc-MPI communicator
  FemusInit mpinit(argc,args,MPI_COMM_WORLD);
  
  Files files; 
        files.CheckIODirectories();
        files.RedirectCout();
	
  /// INIT MESH =================================  
  
  unsigned short nm,nr;
  nm=4;
  std::cout<<"MULTIGRID levels: "<< nm << endl;

  nr=0;
  std::cout<<"MAX_REFINEMENT levels: " << nr << endl<< endl;
  
  int tmp=nm;  nm+=nr;  nr=tmp;
  
  char *infile = new char [50];
 
  sprintf(infile,"./input/nsbenchreg.neu");
  
  //Adimensional quantity (Lref,Uref)
  double Lref = 1.;
  double Uref = 1.;
  
  MultiLevelMesh ml_msh(nm,nr,infile,"seventh",Lref,NULL);
   
  MultiLevelSolution ml_sol(&ml_msh);
   
  // generate solution vector
  ml_sol.AddSolution("U",LAGRANGE,SECOND,2);
  ml_sol.AddSolution("V",LAGRANGE,SECOND,2);
  // the pressure variable should be the last for the Schur decomposition
  ml_sol.AddSolution("P",DISCONTINOUS_POLYNOMIAL,FIRST,1);
  ml_sol.AssociatePropertyToSolution("P","Pressure");
  
  //Initialize (update Init(...) function)
  ml_sol.Initialize("All");
  
  //Set Boundary (update Dirichlet(...) function)
  ml_sol.AttachSetBoundaryConditionFunction(SetBoundaryCondition);
  ml_sol.GenerateBdc("U","Time_dependent");
  ml_sol.GenerateBdc("V");
  ml_sol.GenerateBdc("P");

  
  MultiLevelProblem ml_prob(&ml_sol);
  
  
  Parameter parameter(Lref,Uref);
  
  // Generate fluid Object (Adimensional quantities,viscosity,density,fluid-model)
  Fluid fluid(parameter,0.001,1.,"Newtonian");
  cout << "Fluid properties: " << endl;
  cout << fluid << endl;  

  // add fluid material
  ml_prob.parameters.set<Fluid>("Fluid") = fluid;
  
  
 
  //create systems
  // add the system Navier-Stokes to the MultiLevel problem
  TransientNonlinearImplicitSystem & system = ml_prob.add_system<TransientNonlinearImplicitSystem> ("Navier-Stokes");
  system.AddSolutionToSystemPDE("U");
  system.AddSolutionToSystemPDE("V");
  system.AddSolutionToSystemPDE("P");
  
  // init all the systems
  system.init();
 
  // System Navier-Stokes
  system.SetAssembleFunction(AssembleMatrixResNS);  
  system.SetMaxNumberOfLinearIterations(1);
  system.SetLinearConvergenceTolerance(1.e-8);  
  system.SetMgType(V_CYCLE);
  system.SetMaxNumberOfNonLinearIterations(15);

  // time loop parameter
  system.SetIntervalTime(0.1);
  const unsigned int n_timesteps = 20;
  const unsigned int write_interval = 1;
  
  for (unsigned time_step = 0; time_step < n_timesteps; time_step++) {
   
    // Solving Navier-Stokes system
    std::cout << std::endl;
    std::cout << " *********** Navier-Stokes ************  " << std::endl;
    ml_prob.get_system("Navier-Stokes").solve();
   
    //update Solution
    ml_prob.get_system<TransientNonlinearImplicitSystem>("Navier-Stokes").UpdateSolution();

    // print solution
    if ( !(time_step%write_interval) ) {
        
      //print solution 
//.........这里部分代码省略.........

int main(int argc,char **args) {
  
  /// Init Petsc-MPI communicator
  FemusInit mpinit(argc,args,MPI_COMM_WORLD);

  Files files; 
        files.CheckIODirectories();
        files.RedirectCout();
	
  unsigned short nm,nr;
  std::cout<<"#MULTIGRID levels? (>=1) \n";
  //std::cin>>nm;
  nm=4;

  std::cout<<"#MAX_REFINEMENT levels? (>=0) \n";
  //std::cin>>nr;
  nr=0;
  int tmp=nm;
  nm+=nr;
  nr=tmp;

  const std::string infile = "./input/fsifirst.neu";
  
  double Lref = 1.;
  double Uref = 1.;
  double rhof = 1000.;
  double muf = 1.;
  double rhos = 1000;
  double ni = 0.4;
  double E = 5600000;
  
  MultiLevelMesh ml_msh(nm,nr,infile.c_str(),"fifth",Lref,SetRefinementFlag);
  
  MultiLevelSolution ml_sol(&ml_msh);
  
  //Start System Variables
  ml_sol.AddSolution("DX",LAGRANGE,SECOND,2);
  ml_sol.AddSolution("DY",LAGRANGE,SECOND,2);
//   ml_sol.AssociatePropertyToSolution("DX","Displacement"); // Add this line
//   ml_sol.AssociatePropertyToSolution("DY","Displacement"); // Add this line 
  ml_sol.AddSolution("U",LAGRANGE,SECOND,2);
  ml_sol.AddSolution("V",LAGRANGE,SECOND,2);
  
  //   ml_sol.PairSolution("U","DX"); // Add this line
  //   ml_sol.PairSolution("V","DY"); // Add this line 
  
  ml_sol.AddSolution("AX",LAGRANGE,SECOND,1,0);
  ml_sol.AddSolution("AY",LAGRANGE,SECOND,1,0);
  // Since the Pressure is a Lagrange multiplier it is used as an implicit variable
  ml_sol.AddSolution("P",DISCONTINOUS_POLYNOMIAL,FIRST,1);
  ml_sol.AssociatePropertyToSolution("P","Pressure"); // Add this line

  //Initialize (update Init(...) function)
  ml_sol.Initialize("All");

  //Set Boundary (update Dirichlet(...) function)
  ml_sol.AttachSetBoundaryConditionFunction(SetBoundaryCondition);
  ml_sol.GenerateBdc("DX","Steady");
  ml_sol.GenerateBdc("DY","Steady");
  ml_sol.GenerateBdc("U","Time_dependent");
  ml_sol.GenerateBdc("V","Steady");
  ml_sol.GenerateBdc("AX","Steady");
  ml_sol.GenerateBdc("AY","Steady");
  ml_sol.GenerateBdc("P","Steady");

  
  MultiLevelProblem ml_prob(&ml_sol);
  

  Parameter par(Lref,Uref);
  
  // Generate Solid Object
  Solid solid(par,E,ni,rhos,"Neo-Hookean");
  cout << "Solid properties: " << endl;
  cout << solid << endl;
  
  // Generate Fluid Object
  Fluid fluid(par,muf,rhof,"Newtonian");
  cout << "Fluid properties: " << endl;
  cout << fluid << endl;

  // Add fluid object
  ml_prob.parameters.set<Fluid>("Fluid") = fluid;
  
  // Add Solid Object
  ml_prob.parameters.set<Solid>("Solid") = solid;

  ml_msh.MarkStructureNode();
   
  //create systems
  // add the system FSI to the MultiLevel problem
  TransientMonolithicFSINonlinearImplicitSystem & system = ml_prob.add_system<TransientMonolithicFSINonlinearImplicitSystem> ("Fluid-Structure-Interaction");
  system.AddSolutionToSystemPDE("DX");
  system.AddSolutionToSystemPDE("DY");
  system.AddSolutionToSystemPDE("U");
  system.AddSolutionToSystemPDE("V");
  system.AddSolutionToSystemPDE("P");
  
  // init all the systems
  system.init();
//.........这里部分代码省略.........

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

#ifdef HAVE_LIBMESH
   libMesh::LibMeshInit init(argc,argv);
#else   
   FemusInit init(argc,argv);
#endif
   
 // ======= Files ========================
  Files files; 
        files.ConfigureRestart();
        files.CheckIODirectories();
        files.CopyInputFiles();
        files.RedirectCout();

  // ======= Physics Input Parser ========================
  FemusInputParser<double> physics_map("Physics",files.GetOutputPath());

  const double rhof   = physics_map.get("rho0");
  const double Uref   = physics_map.get("Uref");
  const double Lref   = physics_map.get("Lref");
  const double  muf   = physics_map.get("mu0");

  const double  _pref = rhof*Uref*Uref;           physics_map.set("pref",_pref);
  const double   _Re  = (rhof*Uref*Lref)/muf;     physics_map.set("Re",_Re);
  const double   _Fr  = (Uref*Uref)/(9.81*Lref);  physics_map.set("Fr",_Fr);
  const double   _Pr  = muf/rhof;                 physics_map.set("Pr",_Pr);

  // ======= Mesh =====
  const unsigned NoLevels = 3;
  const unsigned dim = 2;
  const GeomElType geomel_type = QUAD;
  GenCase mesh(NoLevels,dim,geomel_type,"inclQ2D2x2.gam");
          mesh.SetLref(1.);  
	  
  // ======= MyDomainShape  (optional, implemented as child of Domain) ====================
  FemusInputParser<double> box_map("Box",files.GetOutputPath());
  Box mybox(mesh.get_dim(),box_map);
      mybox.InitAndNondimensionalize(mesh.get_Lref());

          mesh.SetDomain(&mybox);    
	  
          mesh.GenerateCase(files.GetOutputPath());

          mesh.SetLref(Lref);
      mybox.InitAndNondimensionalize(mesh.get_Lref());
	  
          XDMFWriter::ReadMeshAndNondimensionalizeBiquadraticHDF5(files.GetOutputPath(),mesh); 
	  XDMFWriter::PrintMeshXDMF(files.GetOutputPath(),mesh,BIQUADR_FE);
          XDMFWriter::PrintMeshLinear(files.GetOutputPath(),mesh);
	  
  //gencase is dimensionalized, meshtwo is nondimensionalized
  //since the meshtwo is nondimensionalized, all the BC and IC are gonna be implemented on a nondimensionalized mesh
  //now, a mesh may or may not have an associated domain
  //moreover, a mesh may or may not be read from file
  //the generation is dimensional, the nondimensionalization occurs later
  //Both the Mesh and the optional domain must be nondimensionalized
  //first, we have to say if the mesh has a shape or not
  //that depends on the application, it must be put at the main level
  //then, after you know the shape, you may or may not generate the mesh with that shape 
  //the two things are totally independent, and related to the application, not to the library

  // ===== QuantityMap : this is like the MultilevelSolution =========================================
  QuantityMap  qty_map;
  qty_map.SetMeshTwo(&mesh);
  qty_map.SetInputParser(&physics_map);

  Pressure       pressure("Qty_Pressure",qty_map,1,LL);             qty_map.AddQuantity(&pressure);
  VelocityX     velocityX("Qty_Velocity0",qty_map,1,QQ);            qty_map.AddQuantity(&velocityX);
  VelocityY     velocityY("Qty_Velocity1",qty_map,1,QQ);            qty_map.AddQuantity(&velocityY);
  Temperature temperature("Qty_Temperature",qty_map,1,QQ);          qty_map.AddQuantity(&temperature);
  TempLift       templift("Qty_TempLift",qty_map,1,QQ);             qty_map.AddQuantity(&templift);  
  TempAdj         tempadj("Qty_TempAdj",qty_map,1,QQ);              qty_map.AddQuantity(&tempadj);  
  // ===== end QuantityMap =========================================
  
  // ====== Start new main =================================
  
  MultiLevelMesh ml_msh;
  ml_msh.GenerateCoarseBoxMesh(8,8,0,0,1,0,2,0,0,QUAD9,"fifth"); //   ml_msh.GenerateCoarseBoxMesh(numelemx,numelemy,numelemz,xa,xb,ya,yb,za,zb,elemtype,"fifth");
  ml_msh.RefineMesh(NoLevels,NoLevels,NULL);
  ml_msh.PrintInfo();
  
  ml_msh.SetWriter(XDMF);
  //ml_msh.GetWriter()->write(files.GetOutputPath(),"biquadratic");
  
  ml_msh.SetDomain(&mybox);    
	  
  MultiLevelSolution ml_sol(&ml_msh);
  ml_sol.AddSolution("Qty_Temperature",LAGRANGE,SECOND,0);
  ml_sol.AddSolution("Qty_TempLift",LAGRANGE,SECOND,0);
  ml_sol.AddSolution("Qty_TempAdj",LAGRANGE,SECOND,0);
  ml_sol.AddSolutionVector(ml_msh.GetDimension(),"Qty_Velocity",LAGRANGE,SECOND,0);
  ml_sol.AddSolution("Qty_Pressure",LAGRANGE,FIRST,0);
  ml_sol.AddSolution("Qty_TempDes",LAGRANGE,SECOND,0,false); //this is not going to be an Unknown! //moreover, this is not going to need any BC (i think they are excluded with "false") // I would like to have a Solution that is NOT EVEN related to the mesh at all... just like a function "on-the-fly"

  // ******* Set problem *******
  MultiLevelProblem ml_prob(&ml_sol);
  ml_prob.SetMeshTwo(&mesh);
  ml_prob.SetQruleAndElemType("fifth");
  ml_prob.SetInputParser(&physics_map);
//.........这里部分代码省略.........

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

  bool Vanka=0, Gmres=0, Asm=0;
  if(argc >= 2) {
    if( !strcmp("vanka",args[1])) 	Vanka=1;
    else if( !strcmp("gmres",args[1])) 	Gmres=1;
    else if( !strcmp("asm",args[1])) 	Asm=1;
    
    if(Vanka+Gmres+Asm==0) {
      cout << "wrong input arguments!" << endl;
      exit(0);
    }
  }
  else {
    cout << "No input argument set default smoother = Gmres" << endl;
    Gmres=1;
  }
  
  /// Init Petsc-MPI communicator
  FemusInit mpinit(argc,args,MPI_COMM_WORLD);
  
  Files files; 
        files.CheckIODirectories();
        files.RedirectCout();

  /// INIT MESH =================================  
  
  unsigned short nm,nr;
  nm=2;
  std::cout<<"MULTIGRID levels: "<< nm << endl;

  nr=0;
  std::cout<<"MAX_REFINEMENT levels: " << nr << endl<< endl;
  
  int tmp=nm;  nm+=nr;  nr=tmp;
  
  char *infile = new char [50];
 
  sprintf(infile,"./input/nsbenc.neu");
  
  //Adimensional quantity (Lref,Uref)
  double Lref = 1.;
  double Uref = 1.;
  
  //Steadystate NonLinearMultiLevelProblem  
  //MultiLevelMesh ml_msh(nm,nr,infile,"seventh",Lref,SetRefinementFlag); 
  MultiLevelMesh ml_msh;
  ml_msh.ReadCoarseMesh(infile,"seventh",Lref);
  ml_msh.RefineMesh(nm,nr,NULL);
  
  // ml_msh.EraseCoarseLevels(2);
  
  MultiLevelSolution ml_sol(&ml_msh);
  
  // generate solution vector
  ml_sol.AddSolution("T",LAGRANGE,SECOND);
  ml_sol.AddSolution("U",LAGRANGE,SECOND);
  ml_sol.AddSolution("V",LAGRANGE,SECOND);
  // the pressure variable should be the last for the Schur decomposition
  ml_sol.AddSolution("P",DISCONTINOUS_POLYNOMIAL,FIRST);
  ml_sol.AssociatePropertyToSolution("P","Pressure");
 
  //Initialize (update Init(...) function)
  ml_sol.Initialize("U",InitVariableU);
  ml_sol.Initialize("V");
  ml_sol.Initialize("P");
  ml_sol.Initialize("T");
  
  //Set Boundary (update Dirichlet(...) function)
  ml_sol.AttachSetBoundaryConditionFunction(SetBoundaryCondition);
  ml_sol.GenerateBdc("U");
  ml_sol.GenerateBdc("V");
  ml_sol.GenerateBdc("P");
  ml_sol.GenerateBdc("T");
  
  MultiLevelProblem ml_prob(&ml_sol);
  
  // add fluid material
  Parameter parameter(Lref,Uref);
  
  // Generate fluid Object (Adimensional quantities,viscosity,density,fluid-model)
  Fluid fluid(parameter,0.001,1,"Newtonian",0.001,1.);
  cout << "Fluid properties: " << endl;
  cout << fluid << endl;
  
  ml_prob.parameters.set<Fluid>("Fluid") = fluid;
   
  
  //BEGIN Navier-Stokes Multilevel Problem
  std::cout << std::endl;
  std::cout << " *********** Navier-Stokes ************  " << std::endl;
    
  NonLinearImplicitSystem & system1 = ml_prob.add_system<NonLinearImplicitSystem> ("Navier-Stokes");
  system1.AddSolutionToSystemPDE("U");
  system1.AddSolutionToSystemPDE("V");
  system1.AddSolutionToSystemPDE("P");
  
  // Set MG Options
  system1.SetAssembleFunction(AssembleMatrixResNS);  
  system1.SetMaxNumberOfNonLinearIterations(3);
//.........这里部分代码省略.........

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

  bool Vanka=0, Gmres=0, Asm=0;
  if(argc >= 2) {
    if( !strcmp("vanka",args[1])) 	Vanka=1;
    else if( !strcmp("gmres",args[1])) 	Gmres=1;
    else if( !strcmp("asm",args[1])) 	Asm=1;
    
    if(Vanka+Gmres+Asm==0) {
      cout << "wrong input arguments!" << endl;
      exit(0);
    }
  }
  else {
    cout << "No input argument set default smoother = Gmres" << endl;
    Gmres=1;
  }
  
  /// Init Petsc-MPI communicator
  FemusInit mpinit(argc,args,MPI_COMM_WORLD);
  
  Files files; 
        files.CheckIODirectories();
        files.RedirectCout();

  /// INIT MESH =================================  
  
  unsigned short nm,nr;
  nm=6;
  std::cout<<"MULTIGRID levels: "<< nm << endl;

  nr=0;
  std::cout<<"MAX_REFINEMENT levels: " << nr << endl<< endl;
  
  int tmp=nm;  nm+=nr;  nr=tmp;
  
  char *infile = new char [50];
 
  sprintf(infile,"./input/nsbenc.neu");
  
  //Adimensional quantity (Lref,Uref)
  double Lref = 1.;
  double Uref = 1.;
  
  MultiLevelMesh ml_msh;
  ml_msh.ReadCoarseMesh(infile,"seventh",Lref);
  ml_msh.RefineMesh(nm,nr,NULL);
  
  ml_msh.PrintInfo();
  
  MultiLevelSolution ml_sol(&ml_msh);
  
  // generate solution vector
//   ml_sol.AddSolution("T",LAGRANGE,SECOND);
//   ml_sol.AddSolution("U",LAGRANGE,SECOND);
//   ml_sol.AddSolution("V",LAGRANGE,SECOND);
//   // the pressure variable should be the last for the Schur decomposition
//   ml_sol.AddSolution("P",DISCONTINOUS_POLYNOMIAL,FIRST);
  
  ml_sol.AddSolution("T",LAGRANGE,FIRST);
  ml_sol.AddSolution("U",LAGRANGE,FIRST);
  ml_sol.AddSolution("V",LAGRANGE,FIRST);
  // the pressure variable should be the last for the Schur decomposition
  ml_sol.AddSolution("P",LAGRANGE,FIRST);
  ml_sol.AssociatePropertyToSolution("P","Pressure");
 
  //Initialize (update Init(...) function)
  ml_sol.Initialize("U",InitVariableU);
  ml_sol.Initialize("V");
  ml_sol.Initialize("P");
  ml_sol.Initialize("T");
  
  //Set Boundary (update Dirichlet(...) function)
  ml_sol.AttachSetBoundaryConditionFunction(SetBoundaryCondition);
  ml_sol.GenerateBdc("U");
  ml_sol.GenerateBdc("V");
  ml_sol.GenerateBdc("P");
  ml_sol.GenerateBdc("T");
  
  MultiLevelProblem ml_prob(&ml_sol);
  
  // add fluid material
  Parameter parameter(Lref,Uref);
  
  // Generate fluid Object (Adimensional quantities,viscosity,density,fluid-model)
  Fluid fluid(parameter,0.001,1,"Newtonian",0.001,1.);
  cout << "Fluid properties: " << endl;
  cout << fluid << endl;
  
  ml_prob.parameters.set<Fluid>("Fluid") = fluid;
   
  
  //BEGIN Stokes Multilevel Problem
  std::cout << std::endl;
  std::cout << " *********** Stokes ************  " << std::endl;
    
  LinearImplicitSystem & system1 = ml_prob.add_system<LinearImplicitSystem> ("Stokes");
  system1.AddSolutionToSystemPDE("U");
  system1.AddSolutionToSystemPDE("V");
  system1.AddSolutionToSystemPDE("P");
//.........这里部分代码省略.........

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

  // ======= Initialize ========================
  // init Petsc-MPI communicator
  FemusInit mpinit(argc, args, MPI_COMM_WORLD);
  
  // ======= Files ========================
  Files files; 
        files.CheckIODirectories();
        files.RedirectCout();

  // ======= Quad Rule ===================
  std::string fe_quad_rule("seventh");

  // ======= Mesh ========================
  //Nondimensional quantity (Lref)
  double Lref = 1.;
  
// 2d
//   const unsigned int nsub_x = 16;
//   const unsigned int nsub_y = 16;
//   const unsigned int nsub_z = 0;
//   const std::vector<double> xyz_min = {0.,0.,0.};
//   const std::vector<double> xyz_max = {1.,1.,0.};
//   const ElemType geom_elem_type = QUAD9;

// 1d
  const unsigned int nsub_x = 16;
  const unsigned int nsub_y = 0;
  const unsigned int nsub_z = 0;
  const std::vector<double> xyz_min = {0.,0.,0.};
  const std::vector<double> xyz_max = {1.,0.,0.};
  const ElemType geom_elem_type = EDGE3;
  
//    std::string input_file = "Lshape_longer_y.med";
//    std::string input_file = "Lshape.med";
//    std::string input_file = "circle_tri6.med";
//    std::string input_file = "ellipse_tri6.med";
//    std::string input_file = "ellipse_with_hole_tri6.med";
   std::string input_file = "interval.med";
   std::ostringstream mystream; mystream << "./" << DEFAULT_INPUTDIR << "/" << input_file;
  const std::string infile = mystream.str();
  
  MultiLevelMesh ml_msh;
    ml_msh.GenerateCoarseBoxMesh(nsub_x,nsub_y,nsub_z,xyz_min[0],xyz_max[0],xyz_min[1],xyz_max[1],xyz_min[2],xyz_max[2],geom_elem_type,fe_quad_rule.c_str());
//   ml_msh.ReadCoarseMesh(infile.c_str(),fe_quad_rule.c_str(),Lref);

  unsigned numberOfUniformLevels = 1;
  unsigned numberOfSelectiveLevels = 0;
  ml_msh.RefineMesh(numberOfUniformLevels , numberOfUniformLevels + numberOfSelectiveLevels, NULL);
  ml_msh.PrintInfo();

  // ======= Solution ========================
  MultiLevelSolution ml_sol(&ml_msh);  // define the multilevel solution and attach the mlMsh object to it
  const unsigned int time_dep_flag = 2;
  
  // ======= Unknowns and Fields ========================
  std::string unknown = "u";
  ml_sol.AddSolution(unknown.c_str(), LAGRANGE, FIRST, time_dep_flag);
  ml_sol.AddSolution("time", DISCONTINUOUS_POLYNOMIAL, ZERO, time_dep_flag);

  // ======= Problem ========================
  MultiLevelProblem ml_prob(&ml_sol);  // define the multilevel problem attach the ml_sol object to it

  ml_prob.SetQuadratureRuleAllGeomElems(fe_quad_rule);
  ml_prob.SetFilesHandler(&files);
  
  // ======= Initial values ========================
  ml_sol.Initialize("All");    // initialize all variables to zero
  ml_sol.Initialize(unknown.c_str(),   SetInitialCondition, &ml_prob);

  // ======= Boundary Conditions ========================
  ml_sol.AttachSetBoundaryConditionFunction(SetBoundaryCondition);  // attach the boundary condition function and generate boundary data
  ml_sol.GenerateBdc(unknown.c_str()); //"Time_dependent");

  // ======= System ========================
  TransientNonlinearImplicitSystem & system = ml_prob.add_system<TransientNonlinearImplicitSystem> ("Timedep");
  system.AddSolutionToSystemPDE(unknown.c_str());

  system.init();  // does it have to stay here or later?

  system.SetAssembleFunction(AssembleMatrixRes_VC);

  system.SetOuterSolver(PREONLY);
  system.SetSolverFineGrids(PREONLY);
  system.SetPreconditionerFineGrids(MLU_PRECOND);
  
//   system.SetMaxNumberOfLinearIterations(1);
//   system.SetAbsoluteLinearConvergenceTolerance(1.e-8);
//   system.SetMgType(V_CYCLE);
  system.SetMaxNumberOfNonLinearIterations(30);
  system.SetNonLinearConvergenceTolerance(1.e-8);

  //**************
  ml_sol.SetWriter(VTK);   //need to move this here for the DebugNonlinear function
  ml_sol.GetWriter()->SetDebugOutput(true);
//   system.SetDebugNonlinear(true);
  //**************
  
  const unsigned fine_lev = ml_sol._mlMesh->GetNumberOfLevels() - 1;
//.........这里部分代码省略.........