static void process_region(img** volume, img** orig, ushort thresh,
	vector<ushort>& rar, vector<ushort>& car, vector<ushort>& zar)
{
	static num dmap[6][512][512];	
	
	// if the region is too small...
	if (rar.size() < 600) {
		cout << "--> deleted an unlikely region (size: " << rar.size() << ")\n";
		return;
	}
	
	// if the region spans too many slices vertically...
	int region_zmin = *min_element(zar.begin(), zar.end());
	if (*max_element(zar.begin(), zar.end()) - region_zmin > 5) {
		dp("--> deleted an unlikely region (z-interval > 5)");
		return;
	}
	
	// find the boundary of the region
	vector<ushort> edge_rar, edge_car, edge_zar;
	
	for (ushort i=0; i < zar.size(); ++i) {
		for (ushort q=zar[i]-1; q < zar[i]+2; ++q) {
			for (ushort k=rar[i]-1; k < rar[i]+2; ++k) {
				for (ushort j=car[i]-1; j < car[i]+2; ++j) {
					if (volume[q]->pix[k][j] == white) {
						edge_rar.push_back(rar[i]);
						edge_car.push_back(car[i]);
						edge_zar.push_back(zar[i]);
						volume[zar[i]]->pix[rar[i]][car[i]] = 0;
						break;
					}
				}
				
				if (volume[zar[i]]->pix[rar[i]][car[i]] == 0) {
					break;
				}
			}
			
			if (volume[zar[i]]->pix[rar[i]][car[i]] == 0) {
				break;
			}
		}
	}
	
	ar("--> region-boundary z-vec", edge_zar, edge_zar.size()); 
	
	// create a distance map
	uint blen = edge_rar.size();
	
	num *erar, *ecar, *ezar, *dvec, *tmp;
	try {
		erar = new num[blen];
		ecar = new num[blen];
		ezar = new num[blen];
		dvec = new num[blen];
		tmp = new num[blen];
	} catch (bad_alloc& err) {
		dp("Fatal memory allocation error; results may be unreliable!");
		return;
	}
	
	for (uint i=0; i < blen; ++i) {
		// copy the boundary voxels into num-arrays to make things simple
		erar[i] = (num) edge_rar[i];
		ecar[i] = (num) edge_car[i];
		ezar[i] = (num) edge_zar[i];
	}
	
	// find the minimum distance from the boundary for each point
	for (uint i=0; i < zar.size(); ++i) {
		mat_op(erar, (num) rar[i], blen, sub, tmp);
		mat_op(tmp, dx, blen, mul, tmp);
		m_apply(tmp, square, blen, dvec);

		mat_op(ecar, (num) car[i], blen, sub, tmp);
		mat_op(tmp, dy, blen, mul, tmp);
		m_apply(tmp, square, blen, tmp);
		mat_op(tmp, dvec, blen, add, dvec);
		
		mat_op(ezar, (num) zar[i], blen, sub, tmp);
		mat_op(tmp, dz, blen, mul, tmp);
		m_apply(tmp, square, blen, tmp);
		mat_op(tmp, dvec, blen, add, dvec);
		
		m_apply(dvec, sqrt, blen, dvec);
		
		// store the distance into the distance map
		dmap[zar[i] - region_zmin][rar[i]][car[i]] = *min_element(dvec, dvec+blen);
	}
	
	// find the 'centroid' of the region
	num cen_max = 0;
	ushort cenz=0, cenx=0, ceny=0;
	
	for (ushort i=0; i < zar.size(); ++i) {
		num cur_val = dmap[zar[i] - region_zmin][rar[i]][car[i]];
		if (cur_val > cen_max) {
			cen_max = cur_val;
			cenz = zar[i] - region_zmin;
//.........这里部分代码省略.........

int main(int argc, char* argv[])
{
  // Load the mesh.
  MeshSharedPtr mesh(new Mesh);
  MeshReaderH2D mloader;
  mloader.load("square.mesh", mesh);

  // Initial mesh refinements.
  for (int i = 0; i < INIT_GLOB_REF_NUM; i++) mesh->refine_all_elements();
  mesh->refine_towards_boundary("Top", INIT_REF_NUM_BDY);

  // Initialize boundary conditions.
  CustomEssentialBCNonConst bc_essential({ "Bottom", "Right", "Top", "Left" });
  EssentialBCs<double> bcs(&bc_essential);

  // Create an H1 space with default shapeset.
  SpaceSharedPtr<double> space(new H1Space<double>(mesh, &bcs, P_INIT));
  int ndof = space->get_num_dofs();
  Hermes::Mixins::Loggable::Static::info("ndof = %d.", ndof);

  // Zero initial solutions. This is why we use H_OFFSET.
  MeshFunctionSharedPtr<double> h_time_prev(new ZeroSolution<double>(mesh));

  // Initialize views.
  ScalarView view("Initial condition", new WinGeom(0, 0, 600, 500));
  view.fix_scale_width(80);

  // Visualize the initial condition.
  view.show(h_time_prev);

  // Initialize the constitutive relations.
  ConstitutiveRelations* constitutive_relations;
  if (constitutive_relations_type == CONSTITUTIVE_GENUCHTEN)
    constitutive_relations = new ConstitutiveRelationsGenuchten(ALPHA, M, N, THETA_S, THETA_R, K_S, STORATIVITY);
  else
    constitutive_relations = new ConstitutiveRelationsGardner(ALPHA, THETA_S, THETA_R, K_S);

  // Initialize the weak formulation.
  double current_time = 0;
  WeakFormSharedPtr<double> wf(new CustomWeakFormRichardsIE(time_step, h_time_prev, constitutive_relations));

  // Initialize the FE problem.
  DiscreteProblem<double> dp(wf, space);

  // Initialize Newton solver.
  NewtonSolver<double> newton(&dp);
  newton.set_verbose_output(true);

  // Time stepping:
  int ts = 1;
  do
  {
    Hermes::Mixins::Loggable::Static::info("---- Time step %d, time %3.5f s", ts, current_time);

    // Perform Newton's iteration.
    try
    {
      newton.set_max_allowed_iterations(NEWTON_MAX_ITER);
      newton.solve();
    }
    catch (Hermes::Exceptions::Exception e)
    {
      e.print_msg();
      throw Hermes::Exceptions::Exception("Newton's iteration failed.");
    };

    // Translate the resulting coefficient vector into the Solution<double> sln->
    Solution<double>::vector_to_solution(newton.get_sln_vector(), space, h_time_prev);

    // Visualize the solution.
    char title[100];
    sprintf(title, "Time %g s", current_time);
    view.set_title(title);
    view.show(h_time_prev);

    // Increase current time and time step counter.
    current_time += time_step;
    ts++;
  } while (current_time < T_FINAL);

  // Wait for the view to be closed.
  View::wait();
  return 0;
}

int main(int argc, char **args) 
{
  // Test variable.
  int success_test = 1;

	if (argc < 3) error("Not enough parameters.");

  // Load the mesh.
	Mesh mesh;
  H3DReader mloader;
  if (!mloader.load(args[1], &mesh)) error("Loading mesh file '%s'.", args[1]);

  // Initialize the space according to the
  // command-line parameters passed.
	int o;
	sscanf(args[2], "%d", &o);
	Ord3 order(o, o, o);
	HcurlSpace space(&mesh, bc_types, NULL, order);

  // Initialize the weak formulation.
	WeakForm wf;
	wf.add_matrix_form(callback(bilinear_form), HERMES_SYM);
	wf.add_vector_form(callback(linear_form));

  // Initialize the FE problem.
  bool is_linear = true;
  DiscreteProblem dp(&wf, &space, is_linear);

  // Set up the solver, matrix, and rhs according to the solver selection.
  SparseMatrix* matrix = create_matrix(matrix_solver);
  Vector* rhs = create_vector(matrix_solver);
  Solver* solver = create_linear_solver(matrix_solver, matrix, rhs);

  // Initialize the preconditioner in the case of SOLVER_AZTECOO.
  if (matrix_solver == SOLVER_AZTECOO) 
  {
    ((AztecOOSolver*) solver)->set_solver(iterative_method);
    ((AztecOOSolver*) solver)->set_precond(preconditioner);
    // Using default iteration parameters (see solver/aztecoo.h).
  }

  // Assemble the linear problem.
  info("Assembling (ndof: %d).", Space::get_num_dofs(&space));
  dp.assemble(matrix, rhs);

  // Solve the linear system. If successful, obtain the solution.
  info("Solving.");
		Solution sln(&mesh);
  if(solver->solve()) Solution::vector_to_solution(solver->get_solution(), &space, &sln);
  else error ("Matrix solver failed.\n");

  ExactSolution ex_sln(&mesh, exact_solution);

  // Calculate exact error.
  info("Calculating exact error.");
  Adapt *adaptivity = new Adapt(&space, HERMES_HCURL_NORM);
  bool solutions_for_adapt = false;
  double err_exact = adaptivity->calc_err_exact(&sln, &ex_sln, solutions_for_adapt, HERMES_TOTAL_ERROR_ABS);

  if (err_exact > EPS)
		// Calculated solution is not precise enough.
		success_test = 0;

  // Clean up.
  delete matrix;
  delete rhs;
  delete solver;
  delete adaptivity;
  
  if (success_test) {
    info("Success!");
    return ERR_SUCCESS;
	}
	else {
    info("Failure!");
    return ERR_FAILURE;
	}
}

int main(int argc, char **args)
{
  // Check the number of command-line parameters.
  if (argc < 2) {
    info("Use x, y, z, xy, xz, yz, or xyz as a command-line parameter.");
    error("Not enough command-line parameters.");
  }

  // Determine anisotropy type from the command-line parameter.
  ANISO_TYPE = parse_aniso_type(args[1]);

  // Load the mesh.
  Mesh mesh;
  H3DReader mesh_loader;
  mesh_loader.load("hex-0-1.mesh3d", &mesh);

  // Assign the lowest possible directional polynomial degrees so that the problem's NDOF >= 1.
  assign_poly_degrees();

  // Create an H1 space with default shapeset.
  info("Setting directional polynomial degrees %d, %d, %d.", P_INIT_X, P_INIT_Y, P_INIT_Z);
  H1Space space(&mesh, bc_types, essential_bc_values, Ord3(P_INIT_X, P_INIT_Y, P_INIT_Z));

  // Initialize weak formulation.
  WeakForm wf;
  wf.add_matrix_form(bilinear_form<double, scalar>, bilinear_form<Ord, Ord>, HERMES_SYM, HERMES_ANY);
  wf.add_vector_form(linear_form<double, scalar>, linear_form<Ord, Ord>, HERMES_ANY);

  // Set exact solution.
  ExactSolution exact(&mesh, fndd);

  // DOF and CPU convergence graphs.
  SimpleGraph graph_dof_est, graph_cpu_est, graph_dof_exact, graph_cpu_exact;

  // Time measurement.
  TimePeriod cpu_time;
  cpu_time.tick();

  // Adaptivity loop. 
  int as = 1; 
  bool done = false;
  do 
  {
    info("---- Adaptivity step %d:", as);

    // Construct globally refined reference mesh and setup reference space.
    Space* ref_space = construct_refined_space(&space,1 , H3D_H3D_H3D_REFT_HEX_XYZ);

    // Initialize discrete problem.
    bool is_linear = true;
    DiscreteProblem dp(&wf, ref_space, is_linear);

    // Set up the solver, matrix, and rhs according to the solver selection.
    SparseMatrix* matrix = create_matrix(matrix_solver);
    Vector* rhs = create_vector(matrix_solver);
    Solver* solver = create_linear_solver(matrix_solver, matrix, rhs);
    
    // Initialize the preconditioner in the case of SOLVER_AZTECOO.
    if (matrix_solver == SOLVER_AZTECOO) 
    {
      ((AztecOOSolver*) solver)->set_solver(iterative_method);
      ((AztecOOSolver*) solver)->set_precond(preconditioner);
      // Using default iteration parameters (see solver/aztecoo.h).
    }
  
    // Assemble the reference problem.
    info("Assembling on reference mesh (ndof: %d).", Space::get_num_dofs(ref_space));
    dp.assemble(matrix, rhs);

    // Time measurement.
    cpu_time.tick();

    // Solve the linear system on reference mesh. If successful, obtain the solution.
    info("Solving on reference mesh.");
    Solution ref_sln(ref_space->get_mesh());
    if(solver->solve()) Solution::vector_to_solution(solver->get_solution(), ref_space, &ref_sln);
    else error ("Matrix solver failed.\n");

    // Time measurement.
    cpu_time.tick();

    // Project the reference solution on the coarse mesh.
    Solution sln(space.get_mesh());
    info("Projecting reference solution on coarse mesh.");
    OGProjection::project_global(&space, &ref_sln, &sln, matrix_solver);

    // Time measurement.
    cpu_time.tick();

    // Output solution and mesh with polynomial orders.
    if (solution_output) 
    {
      out_fn_vtk(&sln, "sln", as);
      out_orders_vtk(&space, "order", as);
    }

    // Skip the visualization time.
    cpu_time.tick(HERMES_SKIP);

    // Calculate element errors and total error estimate.
//.........这里部分代码省略.........

int main(int argc, char* argv[])
{
  // Load the mesh.
  Mesh mesh;
  H2DReader mloader;
  mloader.load("cathedral.mesh", &mesh);

  // Perform initial mesh refinements.
  for(int i = 0; i < INIT_REF_NUM; i++) mesh.refine_all_elements();
  mesh.refine_towards_boundary(BDY_AIR, INIT_REF_NUM_BDY);

  // Enter boundary markers.
  BCTypes bc_types;
  bc_types.add_bc_dirichlet(BDY_GROUND);
  bc_types.add_bc_newton(BDY_AIR);

  // Enter Dirichlet boundary values.
  BCValues bc_values;
  bc_values.add_const(BDY_GROUND, T_INIT);

  // Initialize an H1 space with default shapeset.
  H1Space space(&mesh, &bc_types, &bc_values, P_INIT);
  int ndof = Space::get_num_dofs(&space);
  info("ndof = %d.", ndof);
 
  // Initialize the solution.
  Solution tsln;

  // Set the initial condition.
  tsln.set_const(&mesh, T_INIT);

  // Initialize weak formulation.
  WeakForm wf;
  wf.add_matrix_form(bilinear_form<double, double>, bilinear_form<Ord, Ord>);
  wf.add_matrix_form_surf(bilinear_form_surf<double, double>, bilinear_form_surf<Ord, Ord>, BDY_AIR);
  wf.add_vector_form(linear_form<double, double>, linear_form<Ord, Ord>, HERMES_ANY, &tsln);
  wf.add_vector_form_surf(linear_form_surf<double, double>, linear_form_surf<Ord, Ord>, BDY_AIR);

  // Initialize the FE problem.
  bool is_linear = true;
  DiscreteProblem dp(&wf, &space, is_linear);

  // Set up the solver, matrix, and rhs according to the solver selection.
  SparseMatrix* matrix = create_matrix(matrix_solver);
  Vector* rhs = create_vector(matrix_solver);
  Solver* solver = create_linear_solver(matrix_solver, matrix, rhs);

  // Initialize views.
  ScalarView Tview("Temperature", new WinGeom(0, 0, 450, 600));
  char title[100];
  sprintf(title, "Time %3.5f, exterior temperature %3.5f", TIME, temp_ext(TIME));
  Tview.set_min_max_range(0,20);
  Tview.set_title(title);
  Tview.fix_scale_width(3);

  // Time stepping:
  int nsteps = (int)(FINAL_TIME/TAU + 0.5);
  bool rhs_only = false;
  for(int ts = 1; ts <= nsteps; ts++)
  {
    info("---- Time step %d, time %3.5f, ext_temp %g", ts, TIME, temp_ext(TIME));

    // First time assemble both the stiffness matrix and right-hand side vector,
    // then just the right-hand side vector.
    if (rhs_only == false) info("Assembling the stiffness matrix and right-hand side vector.");
    else info("Assembling the right-hand side vector (only).");
    dp.assemble(matrix, rhs, rhs_only);
    rhs_only = true;

    // Solve the linear system and if successful, obtain the solution.
    info("Solving the matrix problem.");
    if(solver->solve())
      Solution::vector_to_solution(solver->get_solution(), &space, &tsln);
    else 
      error ("Matrix solver failed.\n");

    // Update the time variable.
    TIME += TAU;

    // Visualize the solution.
    sprintf(title, "Time %3.2f, exterior temperature %3.5f", TIME, temp_ext(TIME));
    Tview.set_title(title);
    Tview.show(&tsln);
  }

  // Wait for the view to be closed.
  View::wait();

  // Clean up.
  delete solver;
  delete matrix;
  delete rhs;

  return 0;
}

int main(){
	int numTests;
	cin >> numTests;
	while(numTests--){
		int N; //  assume knowing the length of the array as N
		cin >> N; 
		std::vector<int> v;
		for (int i = 0; i < N; ++i)
		{
			int t;
			cin >> t;
			v.push_back(t);
		}

		if(v.size() <= 1){
			// trivial case
		}
		if(v.size() == 2){
			cout << max(v[0],v[1]) << endl;
			cout << v[0] < v[1] ? 1 : 0 << endl;
			cout << v[0] < v[1] ? v[1] : v[0] << endl;
		}

		// dp[i] is the maximum ending at ith house
		std::vector<int> dp(N,0);
		std::vector<std::vector<int> > record(N, std::vector<int>());
		dp[0] = v[0];
		dp[1] = v[1];
		record[0].push_back(0);
		record[1].push_back(1);

		int final_max = 0, final_idx;
		for (int i = 2; i < N; ++i)
		{
			// dp[i] = max(dp[i-1], dp[i-2]+v[i]);
			
			if (dp[i-1] > dp[i-2] + v[i])
			{

				dp[i] = dp[i-1];
				record[i] = record[i-1];
			}
			else{

				// if not considering the last one, it is exact the same as rob I
				// considering last if only the first and the last is possible exist together

				if (i == N - 1  && record[i-2][0] == 0) // the last one and the first one is already chosen
				{
					if (v[i] > v[0]) // if last one is bigger than first one, simple replace it
					{
						dp[i] = dp[i-2] + v[i] - v[0];
						record[i] = record[i-2];
						record[i].erase(record[i].begin());
						record[i].push_back(i);
					}
					// otherwise do nothing
				}
				else{
					dp[i] = dp[i-2] + v[i];
					record[i] = record[i-2];
					record[i].push_back(i);
				}
			}
			if (final_max < dp[i])
			{
				final_max = dp[i];
				final_idx = i;
			}
		}
		// output the results
		cout << final_max << endl;
		for (int i = 0; i < record[final_idx].size(); ++i)
		{
			cout << record[final_idx][i] << " ";
		}
		cout << endl
		for (int i = 0; i < record[final_idx].size(); ++i)
		{
			cout << v[record[final_idx][i]] << " ";
		}
		cout << endl;

	}
}

int main(int argc, char* argv[])
{
  // Time measurement.
  Hermes::Mixins::TimeMeasurable cpu_time;
  cpu_time.tick();

  // Load the mesh.
  Mesh mesh;
  if (USE_XML_FORMAT == true)
  {
    MeshReaderH2DXML mloader;  
    Hermes::Mixins::Loggable::Static::info("Reading mesh in XML format.");
    mloader.load("domain.xml", &mesh);
  }
  else 
  {
    MeshReaderH2D mloader;
    Hermes::Mixins::Loggable::Static::info("Reading mesh in original format.");
    mloader.load("domain.mesh", &mesh);
  }

  // Perform initial mesh refinements.
  for (int i = 0; i < INIT_REF_NUM; i++) mesh.refine_all_elements();
  
  // Initialize boundary conditions
  CustomEssentialBCNonConst bc_essential("Horizontal");
  EssentialBCs<double> bcs(&bc_essential);

  // Create an H1 space with default shapeset.
  H1Space<double> space(&mesh, &bcs, P_INIT);
  int ndof = space.get_num_dofs();
  Hermes::Mixins::Loggable::Static::info("ndof = %d", ndof);

  // Initialize the weak formulation.
  CustomWeakFormGeneral wf("Horizontal");

  // Initialize the FE problem.
  DiscreteProblem<double> dp(&wf, &space);

  // Initialize Newton solver.
  NewtonSolver<double> newton(&dp);

  // Perform Newton's iteration.
  try
  {
    newton.solve();
  }
  catch(std::exception& e)
  {
    std::cout << e.what();
    
  }

  // Translate the resulting coefficient vector into a Solution.
  Solution<double> sln;
  Solution<double>::vector_to_solution(newton.get_sln_vector(), &space, &sln);

  // Time measurement.
  cpu_time.tick();

  // View the solution and mesh.
  ScalarView sview("Solution", new WinGeom(0, 0, 440, 350));
  sview.show(&sln);
  OrderView oview("Polynomial orders", new WinGeom(450, 0, 405, 350));
  oview.show(&space);

  // Print timing information.
  Hermes::Mixins::Loggable::Static::info("Total running time: %g s", cpu_time.accumulated());

  // Wait for all views to be closed.
  View::wait();

  return 0;
}

int main(int argc, char **args) 
{
  // Time measurement.
  TimePeriod cpu_time;
  cpu_time.tick();

  // Load the mesh. 
  Mesh mesh;
  ExodusIIReader mloader;
  mloader.load("brick_with_holes_tet.e", &mesh);

  // Create H1 space with default shapeset for x-displacement component. 
  H1Space xdisp(&mesh, bc_types_x, essential_bc_values, Ord3(P_INIT));
  
  // Create H1 space with default shapeset for y-displacement component. 
  H1Space ydisp(&mesh, bc_types_y, essential_bc_values, Ord3(P_INIT));
  
  // Create H1 space with default shapeset for z-displacement component. 
  H1Space zdisp(&mesh, bc_types_z, essential_bc_values, Ord3(P_INIT));
  
  // Initialize weak formulation.
  WeakForm wf(3);
  wf.add_matrix_form(0, 0, callback(bilinear_form_0_0), HERMES_SYM);
  wf.add_matrix_form(0, 1, callback(bilinear_form_0_1), HERMES_SYM);
  wf.add_matrix_form(0, 2, callback(bilinear_form_0_2), HERMES_SYM);
  wf.add_vector_form_surf(0, callback(surf_linear_form_x), bdy_force);

  wf.add_matrix_form(1, 1, callback(bilinear_form_1_1), HERMES_SYM);
  wf.add_matrix_form(1, 2, callback(bilinear_form_1_2), HERMES_SYM);
  wf.add_vector_form_surf(1, callback(surf_linear_form_y), bdy_force);

  wf.add_matrix_form(2, 2, callback(bilinear_form_2_2), HERMES_SYM);
  wf.add_vector_form_surf(2, callback(surf_linear_form_z), bdy_force);

  // Initialize discrete problem.
  bool is_linear = true;
  DiscreteProblem dp(&wf, Tuple<Space *>(&xdisp, &ydisp, &zdisp), is_linear);

  // Initialize the solver in the case of SOLVER_PETSC or SOLVER_MUMPS.
  initialize_solution_environment(matrix_solver, argc, args);

  // Set up the solver, matrix, and rhs according to the solver selection.
  SparseMatrix* matrix = create_matrix(matrix_solver);
  Vector* rhs = create_vector(matrix_solver);
  Solver* solver = create_linear_solver(matrix_solver, matrix, rhs);

  // Initialize the preconditioner in the case of SOLVER_AZTECOO.
  if (matrix_solver == SOLVER_AZTECOO) 
  {
    ((AztecOOSolver*) solver)->set_solver(iterative_method);
    ((AztecOOSolver*) solver)->set_precond(preconditioner);
    // Using default iteration parameters (see solver/aztecoo.h).
  }

  // Assemble stiffness matrix and load vector.
  info("Assembling the linear problem (ndof: %d).", Space::get_num_dofs(Tuple<Space *>(&xdisp, &ydisp, &zdisp)));
  dp.assemble(matrix, rhs);

  // Solve the linear system. If successful, obtain the solution.
  info("Solving the linear problem.");
  Solution xsln(xdisp.get_mesh());
  Solution ysln(ydisp.get_mesh());
  Solution zsln(zdisp.get_mesh());
  if(solver->solve()) Solution::vector_to_solutions(solver->get_solution(), 
                      Tuple<Space *>(&xdisp, &ydisp, &zdisp), Tuple<Solution *>(&xsln, &ysln, &zsln));
  else error ("Matrix solver failed.\n");

  // Output all components of the solution.
  if (solution_output) out_fn_vtk(&xsln, &ysln, &zsln, "sln");
  
  // Time measurement.
  cpu_time.tick();

  // Print timing information.
  info("Solutions saved. Total running time: %g s.", cpu_time.accumulated());

  // Clean up.
  delete matrix;
  delete rhs;
  delete solver;
  
  // Properly terminate the solver in the case of SOLVER_PETSC or SOLVER_MUMPS.
  finalize_solution_environment(matrix_solver);

  return 0;
}

int main(int argc, char* argv[])
{
  // Instantiate a class with global functions.
  Hermes2D hermes2d;

  // Load the mesh.
  Mesh mesh, mesh1;
  H2DReader mloader;
  mloader.load("domain.mesh", &mesh);

  // Perform uniform mesh refinement.
  mesh.refine_all_elements();

  // Initialize boundary conditions.
  DefaultEssentialBCConst zero_disp("Bottom", 0.0);
  EssentialBCs bcs(&zero_disp);

  // Create x- and y- displacement space using the default H1 shapeset.
  H1Space u1_space(&mesh, &bcs, P_INIT);
  H1Space u2_space(&mesh, &bcs, P_INIT);
  int ndof = Space::get_num_dofs(Hermes::vector<Space *>(&u1_space, &u2_space));
  info("ndof = %d", ndof);

  // Initialize the weak formulation.
  CustomWeakFormLinearElasticity wf(E, nu, rho*g1, "Top", f0, f1);

  // Initialize the FE problem.
  DiscreteProblem dp(&wf, Hermes::vector<Space *>(&u1_space, &u2_space));

  // Set up the solver, matrix, and rhs according to the solver selection.
  SparseMatrix* matrix = create_matrix(matrix_solver);
  Vector* rhs = create_vector(matrix_solver);
  Solver* solver = create_linear_solver(matrix_solver, matrix, rhs);

  // Initial coefficient vector for the Newton's method.  
  scalar* coeff_vec = new scalar[ndof];
  memset(coeff_vec, 0, ndof*sizeof(scalar));

  // Perform Newton's iteration.
  bool verbose = true;
  bool jacobian_changed = true;
  if (!hermes2d.solve_newton(coeff_vec, &dp, solver, matrix, rhs, jacobian_changed,
      NEWTON_TOL, NEWTON_MAX_ITER, verbose)) error("Newton's iteration failed.");

  // Translate the resulting coefficient vector into the Solution sln.
  Solution u1_sln, u2_sln;
  Solution::vector_to_solutions(coeff_vec, Hermes::vector<Space *>(&u1_space, &u2_space), 
                                Hermes::vector<Solution *>(&u1_sln, &u2_sln));
  
  // Visualize the solution.
  ScalarView view("Von Mises stress [Pa]", new WinGeom(0, 0, 800, 400));
  double lambda = (E * nu) / ((1 + nu) * (1 - 2*nu));  // First Lame constant.
  double mu = E / (2*(1 + nu));                        // Second Lame constant.
  VonMisesFilter stress(Hermes::vector<MeshFunction *>(&u1_sln, &u2_sln), lambda, mu);
  view.show_mesh(false);
  view.show(&stress, HERMES_EPS_HIGH, H2D_FN_VAL_0, &u1_sln, &u2_sln, 1.5e5);

  // Wait for the view to be closed.
  View::wait();

  // Clean up.
  delete [] coeff_vec;
  delete solver;
  delete matrix;
  delete rhs;

  return 0;
}

int main(int argc, char* argv[])
{
  // Choose a Butcher's table or define your own.
  ButcherTable bt(butcher_table_type);
  if (bt.is_explicit()) info("Using a %d-stage explicit R-K method.", bt.get_size());
  if (bt.is_diagonally_implicit()) info("Using a %d-stage diagonally implicit R-K method.", bt.get_size());
  if (bt.is_fully_implicit()) info("Using a %d-stage fully implicit R-K method.", bt.get_size());

  // Load the mesh.
  Mesh mesh, basemesh;
  MeshReaderH2D mloader;
  mloader.load("square.mesh", &basemesh);
  mesh.copy(&basemesh);

  // Initial mesh refinements.
  for(int i = 0; i < INIT_GLOB_REF_NUM; i++) mesh.refine_all_elements();
  mesh.refine_towards_boundary("Top", INIT_REF_NUM_BDY);

  // Initialize boundary conditions.
  CustomEssentialBCNonConst bc_essential(Hermes::vector<std::string>("Bottom", "Right", "Top", "Left"));
  EssentialBCs<double> bcs(&bc_essential);

  // Create an H1 space with default shapeset.
  H1Space<double> space(&mesh, &bcs, P_INIT);
  int ndof_coarse = Space<double>::get_num_dofs(&space);
  info("ndof_coarse = %d.", ndof_coarse);

  // Zero initial solution. This is why we use H_OFFSET.
  ZeroSolution h_time_prev(&mesh), h_time_new(&mesh);

  // Initialize the constitutive relations.
  ConstitutiveRelations* constitutive_relations;
  if(constitutive_relations_type == CONSTITUTIVE_GENUCHTEN)
    constitutive_relations = new ConstitutiveRelationsGenuchten(ALPHA, M, N, THETA_S, THETA_R, K_S, STORATIVITY);
  else
    constitutive_relations = new ConstitutiveRelationsGardner(ALPHA, THETA_S, THETA_R, K_S);

  // Initialize the weak formulation.
  CustomWeakFormRichardsRK wf(constitutive_relations);

  // Initialize the FE problem.
  DiscreteProblem<double> dp(&wf, &space);

  // Create a refinement selector.
  H1ProjBasedSelector<double> selector(CAND_LIST, CONV_EXP, H2DRS_DEFAULT_ORDER);

  // Visualize initial condition.
  char title[100];
  ScalarView view("Initial condition", new WinGeom(0, 0, 440, 350));
  OrderView ordview("Initial mesh", new WinGeom(445, 0, 440, 350));
  view.show(&h_time_prev);
  ordview.show(&space);

  // DOF and CPU convergence graphs initialization.
  SimpleGraph graph_dof, graph_cpu;
  
  // Time measurement.
  TimePeriod cpu_time;
  cpu_time.tick();
  
  // Time stepping loop.
  double current_time = 0; int ts = 1;
  do 
  {
    // Periodic global derefinement.
    if (ts > 1 && ts % UNREF_FREQ == 0) 
    {
      info("Global mesh derefinement.");
      switch (UNREF_METHOD) {
        case 1: mesh.copy(&basemesh);
                space.set_uniform_order(P_INIT);
                break;
        case 2: mesh.unrefine_all_elements();
                space.set_uniform_order(P_INIT);
                break;
        case 3: space.unrefine_all_mesh_elements();
                space.adjust_element_order(-1, -1, P_INIT, P_INIT);
                break;
        default: error("Wrong global derefinement method.");
      }

      ndof_coarse = Space<double>::get_num_dofs(&space);
    }

    // Spatial adaptivity loop. Note: h_time_prev must not be changed 
    // during spatial adaptivity. 
    bool done = false; int as = 1;
    double err_est;
    do {
      info("Time step %d, adaptivity step %d:", ts, as);

      // Construct globally refined reference mesh and setup reference space.
      Space<double>* ref_space = Space<double>::construct_refined_space(&space);
      int ndof_ref = Space<double>::get_num_dofs(ref_space);

      // Time measurement.
      cpu_time.tick();

      // Initialize Runge-Kutta time stepping.
      RungeKutta<double> runge_kutta(&wf, ref_space, &bt, matrix_solver);
//.........这里部分代码省略.........

dgFloat32 dgCollisionChamferCylinder::RayCast(const dgVector& q0, const dgVector& q1, dgFloat32 maxT, dgContactPoint& contactOut, const dgBody* const body, void* const userData, OnRayPrecastAction preFilter) const
{
	if (q0.m_x > m_height) {
		if (q1.m_x < m_height) {
			dgFloat32 t1 = (m_height - q0.m_x) / (q1.m_x - q0.m_x);
			dgFloat32 y = q0.m_y + (q1.m_y - q0.m_y) * t1;
			dgFloat32 z = q0.m_z + (q1.m_z - q0.m_z) * t1;
			if ((y * y + z * z) < m_radius * m_radius) {
				contactOut.m_normal = dgVector(dgFloat32(1.0f), dgFloat32(0.0f), dgFloat32(0.0f), dgFloat32(0.0f));
				return t1;
			}
		}
	}

	if (q0.m_x < -m_height) {
		if (q1.m_x > -m_height) {
			dgFloat32 t1 = (-m_height - q0.m_x) / (q1.m_x - q0.m_x);
			dgFloat32 y = q0.m_y + (q1.m_y - q0.m_y) * t1;
			dgFloat32 z = q0.m_z + (q1.m_z - q0.m_z) * t1;
			if ((y * y + z * z) < m_radius * m_radius) {
				contactOut.m_normal = dgVector(dgFloat32(-1.0f), dgFloat32(0.0f), dgFloat32(0.0f), dgFloat32(0.0f));
				return t1;
			}
		}
	}

	dgVector dq((q1 - q0) & dgVector::m_triplexMask);

	// avoid NaN as a result of a division by zero
	if ((dq % dq) <= 0.0f) {
		return dgFloat32(1.2f);
	}

	dgVector dir(dq.CompProduct4(dq.InvMagSqrt()));
	if (dgAbsf(dir.m_x) > 0.9999f) {
		return dgCollisionConvex::RayCast(q0, q1, maxT, contactOut, body, NULL, NULL);
	}

	dgVector p0(q0 & dgVector::m_triplexMask);
	dgVector p1(q1 & dgVector::m_triplexMask);

	p0.m_x = dgFloat32 (0.0f);
	p1.m_x = dgFloat32 (0.0f);

	dgVector dp (p1 - p0);
	dgFloat32 a = dp.DotProduct4(dp).GetScalar();
	dgFloat32 b = dgFloat32 (2.0f) * dp.DotProduct4(p0).GetScalar();
	dgFloat32 c = p0.DotProduct4(p0).GetScalar() - m_radius * m_radius;

	dgFloat32 disc = b * b - dgFloat32 (4.0f) * a * c;
	if (disc >= dgFloat32 (0.0f)) {
		disc = dgSqrt (disc);
		dgVector origin0(p0 + dp.Scale4 ((-b + disc) / (dgFloat32 (2.0f) * a)));
		dgVector origin1(p0 + dp.Scale4 ((-b - disc) / (dgFloat32 (2.0f) * a)));
		dgFloat32 t0 = dgRayCastSphere(q0, q1, origin0, m_height);
		dgFloat32 t1 = dgRayCastSphere(q0, q1, origin1, m_height);
		if(t1 < t0) {
			t0 = t1;
			origin0 = origin1;
		}

		if ((t0 >= 0.0f) && (t0 <= 1.0f)) {
			contactOut.m_normal = q0 + dq.Scale4(t0) - origin0;
			dgAssert(contactOut.m_normal.m_w == dgFloat32(0.0f));
			contactOut.m_normal = contactOut.m_normal.CompProduct4(contactOut.m_normal.DotProduct4(contactOut.m_normal).InvSqrt());
			return t0;
		}
	} else {
		dgVector origin0 (dgPointToRayDistance (dgVector::m_zero, p0, p1)); 
		origin0 = origin0.Scale4(m_radius / dgSqrt(origin0.DotProduct4(origin0).GetScalar()));
		dgFloat32 t0 = dgRayCastSphere(q0, q1, origin0, m_height);
		if ((t0 >= 0.0f) && (t0 <= 1.0f)) {
			contactOut.m_normal = q0 + dq.Scale4(t0) - origin0;
			dgAssert(contactOut.m_normal.m_w == dgFloat32(0.0f));
			contactOut.m_normal = contactOut.m_normal.CompProduct4(contactOut.m_normal.DotProduct4(contactOut.m_normal).InvSqrt());
			return t0;
		}
	}

	return dgFloat32(1.2f);
}

int main(int argc, char* argv[])
{
  // Time measurement.
  TimePeriod cpu_time;
  cpu_time.tick();

  // Load the mesh.
  Mesh mesh;
  H2DReader mloader;
  mloader.load("domain.mesh", &mesh);

  // Perform initial mesh refinements.
  for (int i=0; i < INIT_REF_NUM; i++) mesh.refine_all_elements();

  // Enter boundary markers.
  BCTypes bc_types;
  bc_types.add_bc_dirichlet(BDY_HORIZONTAL);
  bc_types.add_bc_neumann(BDY_VERTICAL);

  BCValues bc_values;
  bc_values.add_function(BDY_HORIZONTAL, essential_bc_values);

  // Create an H1 space with default shapeset.
  H1Space space(&mesh, &bc_types, &bc_values, P_INIT);
  int ndof = Space::get_num_dofs(&space);
  info("ndof = %d", ndof);

  // Initialize the weak formulation.
  bool matrix_free = false;
  WeakForm wf;
  wf.add_matrix_form(bilinear_form, bilinear_form_ord, HERMES_SYM);
  wf.add_vector_form(linear_form, linear_form_ord);
  wf.add_vector_form_surf(linear_form_surf, linear_form_surf_ord, BDY_VERTICAL);

  // Initialize the FE problem.
  bool is_linear = true;
  DiscreteProblem dp(&wf, &space, is_linear);
   
  SparseMatrix* matrix = create_matrix(matrix_solver);
  Vector* rhs = create_vector(matrix_solver);
  Solver* solver = create_linear_solver(matrix_solver, matrix, rhs);

  if (matrix_solver == SOLVER_AZTECOO) 
  {
    ((AztecOOSolver*) solver)->set_solver(iterative_method);
    ((AztecOOSolver*) solver)->set_precond(preconditioner);
    // Using default iteration parameters (see solver/aztecoo.h).
  }
    
  // Initialize the solution.
  Solution sln;
  
  // Assemble the stiffness matrix and right-hand side vector.
  info("Assembling the stiffness matrix and right-hand side vector.");
  dp.assemble(matrix, rhs);

  // Solve the linear system and if successful, obtain the solution.
  info("Solving the matrix problem.");
  if(solver->solve())
    Solution::vector_to_solution(solver->get_solution(), &space, &sln);
  else
    error ("Matrix solver failed.\n");
  
  // Time measurement.
  cpu_time.tick();
  
  // Clean up.
  delete solver;
  delete matrix;
  delete rhs;

  // View the solution and mesh.
  ScalarView sview("Solution", new WinGeom(0, 0, 440, 350));
  sview.show(&sln);
  OrderView  oview("Polynomial orders", new WinGeom(450, 0, 400, 350));
  oview.show(&space);

  // Skip visualization time.
  cpu_time.tick(HERMES_SKIP);

  // Print timing information.
  verbose("Total running time: %g s", cpu_time.accumulated());

  // Wait for all views to be closed.
  View::wait();
  
  return 0;
}

void StrandBlockSolver::lhsDissipation()
{
  int jj=nPstr+2,nn=nFaces+nBedges;
  Array4D<double> aa(nq,nq,jj,nn);
  for (int n=0; n<nFaces+nBedges; n++)
    for (int j=0; j<nPstr+2; j++)
      for (int k=0; k<nq; k++)
	for (int l=0; l<nq; l++) aa(l,k,j,n) = 0.;


  // unstructured faces
  int c1,c2,n1,n2,jp,fc,m,mm,m1l,m1u,m2l,m2u,npts=1;
  double a[nq*nq];
  for (int n=0; n<nEdges; n++){
    c1            = edge(0,n);
    c2            = edge(1,n);
    m1l           = ncsc(c1  );
    m1u           = ncsc(c1+1);
    m2l           = ncsc(c2  );
    m2u           = ncsc(c2+1);
    fc            = fClip(c1);
    if (fClip(c2) > fc) fc = fClip(c2);
  for (int j=1; j<fc+1; j++){
    sys->lhsDisFluxJacobian(npts,&facs(0,j,n),&xvs(j,n),
			    &q(0,j,c1),&q(0,j,c2),&a[0]);

    for (int k=0; k<nq*nq; k++) a[k] *= .5;

    // diagonal contributions
    for (int k=0; k<nq; k++){
      m             = k*nq;
    for (int l=0; l<nq; l++){
      dd(l,k,j,c1) += a[m+l];
      dd(l,k,j,c2) += a[m+l];
    }}

    // off-diagonal contributions
    for (int k=0; k<nq; k++){
      m             = k*nq;
    for (int l=0; l<nq; l++){
      aa(l,k,j,c1)  = a[m+l];
      aa(l,k,j,c2)  = a[m+l];
    }}

    for (int mm=m1l; mm<m1u; mm++){
      nn            = csc(mm);
      for (int k=0; k<nq; k++){
	m           = k*nq;
      for (int l=0; l<nq; l++)
	bu(l,k,j,mm) -= aa(l,k,j,nn);
      }}

    for (int mm=m2l; mm<m2u; mm++){
      nn            = csc(mm);
      for (int k=0; k<nq; k++){
	m           = k*nq;
      for (int l=0; l<nq; l++)
	bu(l,k,j,mm) -= aa(l,k,j,nn);
      }}

    // reset working array to zero
    for (int k=0; k<nq; k++){
      m             = k*nq;
    for (int l=0; l<nq; l++){
      aa(l,k,j,c1)  = 0.;
      aa(l,k,j,c2)  = 0.;
    }}
  }}

  aa.deallocate();


  // structured faces
  double Ax,Ay,ds,eps=1.e-14;
  for (int n=0; n<nFaces-nGfaces; n++){
    n1            = face(0,n);
    n2            = face(1,n);
  for (int j=0; j<fClip(n)+1; j++){
    jp            = j+1;
    Ax            = facu(0,j,n);
    Ay            = facu(1,j,n);
    ds            = Ax*Ax+Ay*Ay;
    if (fabs(ds) < eps){
      for (int k=0; k<nq; k++){
	m         = k*nq;
      for (int l=0; l<nq; l++)
	a[m+l]    = 0.;
      }}
    else
      sys->lhsDisFluxJacobian(npts,&facu(0,j,n),&xvu(j,n),
			      &q(0,j,n),&q(0,jp,n),&a[0]);
    for (int k=0; k<nq; k++){
      m             = k*nq;
    for (int l=0; l<nq; l++){
      dd(l,k,j ,n) += (a[m+l]*.5);
      dp(l,k,j ,n) -= (a[m+l]*.5);
      dd(l,k,jp,n) += (a[m+l]*.5);
      dm(l,k,jp,n) -= (a[m+l]*.5);
    }}}}
}

int main(int argc, char **args)
{
  // Test variable.
  int success_test = 1;

	if (argc < 2) error("Not enough parameters");

  // Load the mesh.
	Mesh mesh;
  H3DReader mloader;
  if (!mloader.load(args[1], &mesh)) error("Loading mesh file '%s'\n", args[1]);

  // Initialize the space 1.
	Ord3 o1(2, 2, 2);
	H1Space space1(&mesh, bc_types, essential_bc_values, o1);

	// Initialize the space 2.
	Ord3 o2(4, 4, 4);
	H1Space space2(&mesh, bc_types, essential_bc_values, o2);

  // Initialize the weak formulation.
  WeakForm wf(2);
	wf.add_matrix_form(0, 0, bilinear_form_1<double, scalar>, bilinear_form_1<Ord, Ord>, HERMES_SYM);
	wf.add_vector_form(0, linear_form_1<double, scalar>, linear_form_1<Ord, Ord>);
	wf.add_matrix_form(1, 1, bilinear_form_2<double, scalar>, bilinear_form_2<Ord, Ord>, HERMES_SYM);
	wf.add_vector_form(1, linear_form_2<double, scalar>, linear_form_2<Ord, Ord>);

  // Initialize the FE problem.
  bool is_linear = true;
  DiscreteProblem dp(&wf, Tuple<Space *>(&space1, &space2), is_linear);

  // Initialize the solver in the case of SOLVER_PETSC or SOLVER_MUMPS.
  initialize_solution_environment(matrix_solver, argc, args);

  // Set up the solver, matrix, and rhs according to the solver selection.
  SparseMatrix* matrix = create_matrix(matrix_solver);
  Vector* rhs = create_vector(matrix_solver);
  Solver* solver = create_linear_solver(matrix_solver, matrix, rhs);

  // Initialize the preconditioner in the case of SOLVER_AZTECOO.
  if (matrix_solver == SOLVER_AZTECOO) 
  {
    ((AztecOOSolver*) solver)->set_solver(iterative_method);
    ((AztecOOSolver*) solver)->set_precond(preconditioner);
    // Using default iteration parameters (see solver/aztecoo.h).
  }

  // Assemble the linear problem.
  info("Assembling (ndof: %d).", Space::get_num_dofs(Tuple<Space *>(&space1, &space2)));
  dp.assemble(matrix, rhs);
    
  // Solve the linear system. If successful, obtain the solution.
  info("Solving.");
  Solution sln1(&mesh);
  Solution sln2(&mesh);
  if(solver->solve()) Solution::vector_to_solutions(solver->get_solution(), Tuple<Space *>(&space1, &space2), Tuple<Solution *>(&sln1, &sln2));
  else error ("Matrix solver failed.\n");
    
  ExactSolution ex_sln1(&mesh, exact_sln_fn_1);
  ExactSolution ex_sln2(&mesh, exact_sln_fn_2);

  // Calculate exact error.
  info("Calculating exact error.");
  Adapt *adaptivity = new Adapt(Tuple<Space *>(&space1, &space2), Tuple<ProjNormType>(HERMES_H1_NORM, HERMES_H1_NORM));
  bool solutions_for_adapt = false;
  double err_exact = adaptivity->calc_err_exact(Tuple<Solution *>(&sln1, &sln2), Tuple<Solution *>(&ex_sln1, &ex_sln2), solutions_for_adapt, HERMES_TOTAL_ERROR_ABS);

  if (err_exact > EPS)
		// Calculated solution is not precise enough.
		success_test = 0;

  // Clean up.
  delete matrix;
  delete rhs;
  delete solver;
  delete adaptivity;

  // Properly terminate the solver in the case of SOLVER_PETSC or SOLVER_MUMPS.
  finalize_solution_environment(matrix_solver);
  
  if (success_test) {
    info("Success!");
    return ERR_SUCCESS;
  }
  else {
    info("Failure!");
    return ERR_FAILURE;
  }
}

void writeFrames(const Kinect::FrameSource::IntrinsicParameters& ip,const Kinect::FrameBuffer& color,const Kinect::MeshBuffer& mesh,const char* lwoFileName)
	{
	/* Create the texture file name: */
	std::string textureFileName(lwoFileName,Misc::getExtension(lwoFileName));
	textureFileName.append("-color.png");
	
	/* Write the color frame as a texture image: */
	{
	Images::RGBImage texImage(color.getSize(0),color.getSize(1));
	Images::RGBImage::Color* tiPtr=texImage.modifyPixels();
	const unsigned char* cfPtr=color.getTypedBuffer<unsigned char>();
	for(int y=0;y<color.getSize(1);++y)
		for(int x=0;x<color.getSize(0);++x,++tiPtr,cfPtr+=3)
			*tiPtr=Images::RGBImage::Color(cfPtr);
	
	Images::writeImageFile(texImage,textureFileName.c_str());
	}
	
	/* Open the LWO file: */
	IO::FilePtr lwoFile=IO::openFile(lwoFileName,IO::File::WriteOnly);
	lwoFile->setEndianness(Misc::BigEndian);
	
	/* Create the LWO file structure via the FORM chunk: */
	{
	IFFChunkWriter form(lwoFile,"FORM");
	form.write<char>("LWO2",4);
	
	/* Create the TAGS chunk: */
	{
	IFFChunkWriter tags(&form,"TAGS");
	tags.writeString("ColorImage");
	tags.writeChunk();
	}
	
	/* Create the LAYR chunk: */
	{
	IFFChunkWriter layr(&form,"LAYR");
	layr.write<Misc::UInt16>(0U);
	layr.write<Misc::UInt16>(0x0U);
	for(int i=0;i<3;++i)
		layr.write<Misc::Float32>(0.0f);
	layr.writeString("DepthImage");
	layr.writeChunk();
	}
	
	/* Create an index map for all vertices to omit unused vertices: */
	unsigned int* indices=new unsigned int[mesh.numVertices];
	for(unsigned int i=0;i<mesh.numVertices;++i)
		indices[i]=~0x0U;
	unsigned int numUsedVertices=0;
	
	/* Create the PNTS, BBOX and VMAP chunks in one go: */
	{
	typedef Kinect::FrameSource::IntrinsicParameters::PTransform PTransform;
	typedef PTransform::Point Point;
	typedef Geometry::Box<Point::Scalar,3> Box;
	
	IFFChunkWriter bbox(&form,"BBOX");
	IFFChunkWriter pnts(&form,"PNTS");
	IFFChunkWriter vmap(&form,"VMAP");
	
	/* Write the VMAP header: */
	vmap.write<char>("TXUV",4);
	vmap.write<Misc::UInt16>(2U);
	vmap.writeString("ColorImageUV");
	
	/* Process all triangle vertices: */
	Box pBox=Box::empty;
	const Kinect::MeshBuffer::Vertex* vertices=mesh.getVertices();
	const Kinect::MeshBuffer::Index* tiPtr=mesh.getTriangleIndices();
	for(unsigned int i=0;i<mesh.numTriangles*3;++i,++tiPtr)
		{