C++ Mesh类代码示例

 DLL_HEADER int Ng_GetNSeg_2D (Ng_Mesh * mesh)
 {
    Mesh * m = (Mesh*)mesh;
    return m->GetNSeg();
 }

Mesh* FloorBuildingF::toMesh()
{
	//ROOF6_LEFT_UL_X;

	Mesh* m = new Mesh();
	m->addVertex(ul, ROOF6_LEFT_UL_X/divide, ROOF6_LEFT_UL_Y/divide); //0
	m->addVertex(ur, ROOF6_RIGHT_UR_X/divide, ROOF6_RIGHT_UR_Y/divide); //1
	m->addVertex(dr, ROOF6_RIGHT_DR_X/divide, ROOF6_RIGHT_DR_Y/divide); //2
	m->addVertex(dl, ROOF6_LEFT_DL_X/divide, ROOF6_LEFT_DL_Y/divide); //3

	Vector border_ul, border_ur, border_dr, border_dl, norm;
	norm = ur - ul;
	norm = Normalized(norm);
	border_ul = ul + norm*BORDER;
	//border_ul[2] += BORDER;
	border_dl = border_ul;
	border_dl[1] = dl[1];
	border_ur = ur;
	//border_ur[2] -= BORDER;
	border_ur = ur - norm*BORDER;
	border_dr = border_ur;
	border_dr[1] = dr[1];

	m->addVertex(border_ul, ROOF6_LEFT_UR_X/divide, ROOF6_LEFT_UR_Y/divide); //4
	m->addVertex(border_ur, ROOF6_RIGHT_UL_X/divide, ROOF6_RIGHT_UL_Y/divide); //5
	m->addVertex(border_dr, ROOF6_RIGHT_DL_X/divide, ROOF6_RIGHT_DL_Y/divide); //6
	m->addVertex(border_dl, ROOF6_LEFT_DR_X/divide, ROOF6_LEFT_DR_Y/divide); //7

	// m->addVertex(border_ul, ROOF6_MID_UL_X/divide, ROOF6_MID_UL_Y/divide); //8
	// m->addVertex(border_ur, ROOF6_MID_UR_X/divide, ROOF6_MID_UR_Y/divide); //9
	// m->addVertex(border_dr, ROOF6_MID_DR_X/divide, ROOF6_MID_DR_Y/divide); //10
	// m->addVertex(border_dl, ROOF6_MID_DL_X/divide, ROOF6_MID_DL_Y/divide); //11


	// m->addFace(0, 1, 3);
	// m->addFace(1, 2, 3);
/*

	m->addVertex(ul, ROOF6_LEFT_UL_X/divide, ROOF6_LEFT_UL_Y/divide); //0
	m->addVertex(ur, ROOF6_RIGHT_UR_X/divide, ROOF6_RIGHT_UR_Y/divide); //1
	m->addVertex(dr, ROOF6_RIGHT_DR_X/divide, ROOF6_RIGHT_DR_Y/divide); //2
	m->addVertex(dl, ROOF6_LEFT_DL_X/divide, ROOF6_LEFT_DL_Y/divide); //3

	m->addVertex(border_ul, ROOF6_LEFT_UR_X/divide, ROOF6_LEFT_UR_Y/divide); //4
	m->addVertex(border_ur, ROOF6_RIGHT_UL_X/divide, ROOF6_RIGHT_UL_Y/divide); //5
	m->addVertex(border_dr, ROOF6_RIGHT_DL_X/divide, ROOF6_RIGHT_DL_Y/divide); //6
	m->addVertex(border_dl, ROOF6_LEFT_DR_X/divide, ROOF6_LEFT_DR_Y/divide); //7

	m->addVertex(border_ul, ROOF6_MID_UL_X/divide, ROOF6_MID_UL_Y/divide); //8
	m->addVertex(border_ur, ROOF6_MID_UR_X/divide, ROOF6_MID_UR_Y/divide); //9
	m->addVertex(border_dr, ROOF6_MID_DR_X/divide, ROOF6_MID_DR_Y/divide); //10
	m->addVertex(border_dl, ROOF6_MID_DL_X/divide, ROOF6_MID_DL_Y/divide); //11

*/
	m->addFace(0, 4, 3);//, ROOF6_LEFT_UL_X/divide, ROOF6_LEFT_UL_Y/divide, ROOF6_LEFT_UR_X/divide, ROOF6_LEFT_UR_Y/divide, ROOF6_LEFT_DL_X/divide, ROOF6_LEFT_DL_Y/divide);
	m->addFace(3, 4, 7);//, ROOF6_LEFT_DL_X/divide, ROOF6_LEFT_DL_Y/divide, border_ul, ROOF6_LEFT_UR_X/divide, ROOF6_LEFT_UR_Y/divide, ROOF6_LEFT_DR_X/divide, ROOF6_LEFT_DR_Y/divide);
	m->addFace(5, 1, 6);//, ROOF6_RIGHT_UL_X/divide, ROOF6_RIGHT_UL_Y/divide, ur, ROOF6_RIGHT_UR_X/divide, ROOF6_RIGHT_UR_Y/divide, ROOF6_RIGHT_DL_X/divide, ROOF6_RIGHT_DL_Y/divide);
	m->addFace(6, 1, 2);//, ROOF6_RIGHT_DL_X/divide, ROOF6_RIGHT_DL_Y/divide, ROOF6_RIGHT_UR_X/divide, ROOF6_RIGHT_UR_Y/divide, ROOF6_RIGHT_DR_X/divide, ROOF6_RIGHT_DR_Y/divide);
	m->addFace(4, 5, 7);//, ROOF6_MID_UL_X/divide, ROOF6_MID_UL_Y/divide, ROOF6_MID_UR_X/divide, ROOF6_MID_UR_Y/divide, ROOF6_MID_DL_X/divide, ROOF6_MID_DL_Y/divide);
	m->addFace(7, 5, 6);//, ROOF6_MID_DL_X/divide, ROOF6_MID_DL_Y/divide, ROOF6_MID_UR_X/divide, ROOF6_MID_UR_Y/divide, ROOF6_MID_DR_X/divide, ROOF6_MID_DR_Y/divide);

	return m;

	//m->addVertex(0.0, height/2, 0.0);
	//m->addFace(0, 2+i, 2+i+1);
}

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

  // Enter boundary markers.
  BCTypes bc_types;
  bc_types.add_bc_dirichlet(Hermes::vector<int>(BDY_BOTTOM, BDY_OUTER, BDY_LEFT, BDY_INNER));

  // Enter Dirichlet boudnary values.
  BCValues bc_values;
  bc_values.add_function(Hermes::vector<int>(BDY_BOTTOM, BDY_OUTER, BDY_LEFT, BDY_INNER), 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.
  WeakForm wf;
  wf.add_matrix_form(callback(bilinear_form));
  wf.add_vector_form(callback(linear_form));

  // Testing n_dof and correctness of solution vector
  // for p_init = 1, 2, ..., 10
  int success = 1;
  Solution sln;
  for (int p_init = 1; p_init <= 10; p_init++) {

    printf("********* p_init = %d *********\n", p_init);
    space.set_uniform_order(p_init);

    // 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 solution.
    Solution sln;

    // Assemble the stiffness matrix and right-hand side vector.
    info("Assembling the stiffness matrix and right-hand side vector.");
    bool rhsonly = false;
    dp.assemble(matrix, rhs, rhsonly);

    // 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");

    int ndof = Space::get_num_dofs(&space);
    printf("ndof = %d\n", ndof);
    double sum = 0;
    for (int i=0; i < ndof; i++) sum += solver->get_solution()[i];
    printf("coefficient sum = %g\n", sum);

    // Actual test. The values of 'sum' depend on the
    // current shapeset. If you change the shapeset,
    // you need to correct these numbers.
    if (p_init == 1 && fabs(sum - 1.7251) > 1e-3) success = 0;
    if (p_init == 2 && fabs(sum - 3.79195) > 1e-3) success = 0;
    if (p_init == 3 && fabs(sum - 3.80206) > 1e-3) success = 0;
    if (p_init == 4 && fabs(sum - 3.80156) > 1e-3) success = 0;
    if (p_init == 5 && fabs(sum - 3.80155) > 1e-3) success = 0;
    if (p_init == 6 && fabs(sum - 3.80154) > 1e-3) success = 0;
    if (p_init == 7 && fabs(sum - 3.80154) > 1e-3) success = 0;
    if (p_init == 8 && fabs(sum - 3.80153) > 1e-3) success = 0;
    if (p_init == 9 && fabs(sum - 3.80152) > 1e-3) success = 0;
    if (p_init == 10 && fabs(sum - 3.80152) > 1e-3) success = 0;
  }

  if (success == 1) {
    printf("Success!\n");
    return ERR_SUCCESS;
  }
  else {
    printf("Failure!\n");
    return ERR_FAILURE;
  }
}

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

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

  // Boundary condition types;
  BCTypes bc_types;

  // Initialize boundary condition types and spaces with default shapesets.
  bc_types.add_bc_neumann(Hermes::vector<int>(BDY_SOLID_WALL, BDY_INLET_OUTLET));
  L2Space space_rho(&mesh, &bc_types, P_INIT);
  L2Space space_rho_v_x(&mesh, &bc_types, P_INIT);
  L2Space space_rho_v_y(&mesh, &bc_types, P_INIT);
  L2Space space_e(&mesh, &bc_types, P_INIT);

  // Initialize solutions, set initial conditions.
  Solution prev_rho, prev_rho_v_x, prev_rho_v_y, prev_e;
  prev_rho.set_exact(&mesh, ic_density);
  prev_rho_v_x.set_exact(&mesh, ic_density_vel_x);
  prev_rho_v_y.set_exact(&mesh, ic_density_vel_y);
  prev_e.set_exact(&mesh, ic_energy);

  // Solutions for the time derivative estimate.
  Solution sln_temp_rho, sln_temp_rho_v_x, sln_temp_rho_v_y, sln_temp_e;

  // Initialize weak formulation.
  bool is_matrix_free = true;
  WeakForm wf(4, is_matrix_free);
  
  // Volumetric linear forms.
  wf.add_vector_form(0, callback(linear_form_0_time));
  wf.add_vector_form(1, callback(linear_form_1_time));
  wf.add_vector_form(2, callback(linear_form_2_time));
  wf.add_vector_form(3, callback(linear_form_3_time));

  // Volumetric linear forms.
  // Linear forms coming from the linearization by taking the Eulerian fluxes' Jacobian matrices 
  // from the previous time step.
  // Unnecessary for FVM.
  if(P_INIT.order_h > 0 || P_INIT.order_v > 0) {
    // First flux.
    wf.add_vector_form(0, callback(linear_form_0_1), HERMES_ANY);
    
    wf.add_vector_form(1, callback(linear_form_1_0_first_flux), HERMES_ANY);
    wf.add_vector_form(1, callback(linear_form_1_1_first_flux), HERMES_ANY);
    wf.add_vector_form(1, callback(linear_form_1_2_first_flux), HERMES_ANY);
    wf.add_vector_form(1, callback(linear_form_1_3_first_flux), HERMES_ANY);
    wf.add_vector_form(2, callback(linear_form_2_0_first_flux), HERMES_ANY);
    wf.add_vector_form(2, callback(linear_form_2_1_first_flux), HERMES_ANY);
    wf.add_vector_form(2, callback(linear_form_2_2_first_flux), HERMES_ANY);
    wf.add_vector_form(2, callback(linear_form_2_3_first_flux), HERMES_ANY);
    wf.add_vector_form(3, callback(linear_form_3_0_first_flux), HERMES_ANY);
    wf.add_vector_form(3, callback(linear_form_3_1_first_flux), HERMES_ANY);
    wf.add_vector_form(3, callback(linear_form_3_2_first_flux), HERMES_ANY);
    wf.add_vector_form(3, callback(linear_form_3_3_first_flux), HERMES_ANY);

    // Second flux.
    wf.add_vector_form(0, callback(linear_form_0_2), HERMES_ANY);
    wf.add_vector_form(1, callback(linear_form_1_0_second_flux), HERMES_ANY);
    wf.add_vector_form(1, callback(linear_form_1_1_second_flux), HERMES_ANY);
    wf.add_vector_form(1, callback(linear_form_1_2_second_flux), HERMES_ANY);
    wf.add_vector_form(1, callback(linear_form_1_3_second_flux), HERMES_ANY);
    wf.add_vector_form(2, callback(linear_form_2_0_second_flux), HERMES_ANY);
    wf.add_vector_form(2, callback(linear_form_2_1_second_flux), HERMES_ANY);
    wf.add_vector_form(2, callback(linear_form_2_2_second_flux), HERMES_ANY);
    wf.add_vector_form(2, callback(linear_form_2_3_second_flux), HERMES_ANY);
    wf.add_vector_form(3, callback(linear_form_3_0_second_flux), HERMES_ANY);
    wf.add_vector_form(3, callback(linear_form_3_1_second_flux), HERMES_ANY);
    wf.add_vector_form(3, callback(linear_form_3_2_second_flux), HERMES_ANY);
    wf.add_vector_form(3, callback(linear_form_3_3_second_flux), HERMES_ANY);
  }

  // Volumetric linear forms coming from the time discretization.
  wf.add_vector_form(0, linear_form_time, linear_form_order, HERMES_ANY, &prev_rho);
  wf.add_vector_form(1, linear_form_time, linear_form_order, HERMES_ANY, &prev_rho_v_x);
  wf.add_vector_form(2, linear_form_time, linear_form_order, HERMES_ANY, &prev_rho_v_y);
  wf.add_vector_form(3, linear_form_time, linear_form_order, HERMES_ANY, &prev_e);

  // Surface linear forms - inner edges coming from the DG formulation.
  wf.add_vector_form_surf(0, linear_form_interface_0, linear_form_order, H2D_DG_INNER_EDGE);
  wf.add_vector_form_surf(1, linear_form_interface_1, linear_form_order, H2D_DG_INNER_EDGE);
  wf.add_vector_form_surf(2, linear_form_interface_2, linear_form_order, H2D_DG_INNER_EDGE);
  wf.add_vector_form_surf(3, linear_form_interface_3, linear_form_order, H2D_DG_INNER_EDGE);

  // Surface linear forms - inlet / outlet edges.
  wf.add_vector_form_surf(0, bdy_flux_inlet_outlet_comp_0, linear_form_order, BDY_INLET_OUTLET);
  wf.add_vector_form_surf(1, bdy_flux_inlet_outlet_comp_1, linear_form_order, BDY_INLET_OUTLET);
  wf.add_vector_form_surf(2, bdy_flux_inlet_outlet_comp_2, linear_form_order, BDY_INLET_OUTLET);
  wf.add_vector_form_surf(3, bdy_flux_inlet_outlet_comp_3, linear_form_order, BDY_INLET_OUTLET);

  // Surface linear forms - Solid wall edges.
  wf.add_vector_form_surf(0, bdy_flux_solid_wall_comp_0, linear_form_order, BDY_SOLID_WALL);
  wf.add_vector_form_surf(1, bdy_flux_solid_wall_comp_1, linear_form_order, BDY_SOLID_WALL);
  wf.add_vector_form_surf(2, bdy_flux_solid_wall_comp_2, linear_form_order, BDY_SOLID_WALL);
  wf.add_vector_form_surf(3, bdy_flux_solid_wall_comp_3, linear_form_order, BDY_SOLID_WALL);
//.........这里部分代码省略.........

bool FileSystem::parseObj(const char *path, Model &m)
{
	std::string line;
	std::ifstream fileStream(path);
	unsigned int ind = 0; //amout of indices
	Mesh mesh;

	std::vector<glm::vec3> v, n;
	std::vector<glm::vec2> t;

	if (fileStream.is_open())
	{
		while (fileStream.good())
		{
			getline(fileStream, line);

			std::vector<std::string> tokens = split(line, ' ');

			if (line.substr(0, 2) == "v ")
			{
				float x = (float)atof(tokens[1].c_str());
				float y = (float)atof(tokens[2].c_str());
				float z = (float)atof(tokens[3].c_str());
				v.push_back(glm::vec3(x, y, z));
			}

			else if (line.substr(0, 2) == "vt")
			{
				float x = (float)atof(tokens[1].c_str());
				float y = (float)atof(tokens[2].c_str());
				t.push_back(glm::vec2(x, y));
			}

			else if (line.substr(0, 2) == "vn")
			{
				float x = (float)atof(tokens[1].c_str());
				float y = (float)atof(tokens[2].c_str());
				float z = (float)atof(tokens[3].c_str());
				n.push_back(glm::vec3(x, y, z));
			}

			else if (line.substr(0, 2) == "f ")
			{
				glm::vec3 vertex, normal;
				glm::vec2 texCoord;

				std::vector<std::string> v1 = split(tokens[1], '/');
				int i = atoi(v1[0].c_str()) - 1;
				vertex.x = v[i].x; vertex.y = v[i].y; vertex.z = v[i].z;
				i = atoi(v1[1].c_str()) - 1;
				texCoord.x = t[i].x; texCoord.y = t[i].y;
				i = atoi(v1[2].c_str()) - 1;
				normal.x = n[i].x; normal.y = n[i].y; normal.z = n[i].z;
				mesh.addVertex(vertex, normal, texCoord);

				std::vector<std::string> v2 = split(tokens[2], '/');
				i = atoi(v2[0].c_str()) - 1;
				vertex.x = v[i].x; vertex.y = v[i].y; vertex.z = v[i].z;
				i = atoi(v2[1].c_str()) - 1;
				texCoord.x = t[i].x; texCoord.y = t[i].y;
				i = atoi(v2[2].c_str()) - 1;
				normal.x = n[i].x; normal.y = n[i].y; normal.z = n[i].z;
				mesh.addVertex(vertex, normal, texCoord);

				std::vector<std::string> v3 = split(tokens[3], '/');
				i = atoi(v3[0].c_str()) - 1;
				vertex.x = v[i].x; vertex.y = v[i].y; vertex.z = v[i].z;
				i = atoi(v3[1].c_str()) - 1;
				texCoord.x = t[i].x; texCoord.y = t[i].y;
				i = atoi(v3[2].c_str()) - 1;
				normal.x = n[i].x; normal.y = n[i].y; normal.z = n[i].z;
				mesh.addVertex(vertex, normal, texCoord);
			}
		}

		std::shared_ptr<Material> mat(new Material);
		Texture tex;
		if (!loadTexture("resource/white_D.png", tex))
			return false;
		mat->setDifTex(tex);
		mesh.addMaterial(mat);

		m.addMesh(std::make_shared<Mesh>(mesh));

		fileStream.close();
		return true;
	}

	std::cout << "Unable to open file " << path << std::endl;
	return false;
}

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 refinemets.
  for (int i=0; i < INIT_REF_NUM; i++)  mesh.refine_all_elements();

  // Initialize boundary conditions.
  BCTypes bc_types;
  bc_types.add_bc_dirichlet(BDY_LEFT);
  bc_types.add_bc_neumann(BDY_NEUMANN);
  bc_types.add_bc_newton(BDY_COOLED);

  // Enter Dirichlet boundary values.
  BCValues bc_values_t;
  bc_values_t.add_const(BDY_LEFT, 1.0);

  BCValues bc_values_c;
  bc_values_c.add_zero(BDY_LEFT);

  // Create H1 spaces with default shapesets.
  H1Space* t_space = new H1Space(&mesh, &bc_types, &bc_values_t, P_INIT);
  H1Space* c_space = new H1Space(&mesh, &bc_types, &bc_values_c, P_INIT);
  int ndof = Space::get_num_dofs(Hermes::vector<Space *>(t_space, c_space));
  info("ndof = %d.", ndof);

  // Define initial conditions.
  Solution t_prev_time_1, c_prev_time_1, t_prev_time_2, 
           c_prev_time_2, t_iter, c_iter, t_prev_newton, c_prev_newton;
  t_prev_time_1.set_exact(&mesh, temp_ic);  
  c_prev_time_1.set_exact(&mesh, conc_ic);
  t_prev_time_2.set_exact(&mesh, temp_ic);  
  c_prev_time_2.set_exact(&mesh, conc_ic);
  t_iter.set_exact(&mesh, temp_ic);   
  c_iter.set_exact(&mesh, conc_ic);

  // Filters for the reaction rate omega and its derivatives.
  DXDYFilter omega(omega_fn, Hermes::vector<MeshFunction*>(&t_prev_time_1, &c_prev_time_1));
  DXDYFilter omega_dt(omega_dt_fn, Hermes::vector<MeshFunction*>(&t_prev_time_1, &c_prev_time_1));
  DXDYFilter omega_dc(omega_dc_fn, Hermes::vector<MeshFunction*>(&t_prev_time_1, &c_prev_time_1));

  // Initialize weak formulation.
  WeakForm wf(2, JFNK ? true : false);
  if (!JFNK || (JFNK && PRECOND == 1))
  {
    wf.add_matrix_form(0, 0, callback(newton_bilinear_form_0_0), HERMES_NONSYM, HERMES_ANY, &omega_dt);
    wf.add_matrix_form_surf(0, 0, callback(newton_bilinear_form_0_0_surf), 3);
    wf.add_matrix_form(1, 1, callback(newton_bilinear_form_1_1), HERMES_NONSYM, HERMES_ANY, &omega_dc);
    wf.add_matrix_form(0, 1, callback(newton_bilinear_form_0_1), HERMES_NONSYM, HERMES_ANY, &omega_dc);
    wf.add_matrix_form(1, 0, callback(newton_bilinear_form_1_0), HERMES_NONSYM, HERMES_ANY, &omega_dt);
  }
  else if (PRECOND == 2)
  {
    wf.add_matrix_form(0, 0, callback(precond_0_0));
    wf.add_matrix_form(1, 1, callback(precond_1_1));
  }
  wf.add_vector_form(0, callback(newton_linear_form_0), HERMES_ANY, 
                     Hermes::vector<MeshFunction*>(&t_prev_time_1, &t_prev_time_2, &omega));
  wf.add_vector_form_surf(0, callback(newton_linear_form_0_surf), 3);
  wf.add_vector_form(1, callback(newton_linear_form_1), HERMES_ANY, 
                     Hermes::vector<MeshFunction*>(&c_prev_time_1, &c_prev_time_2, &omega));

  // Project the functions "t_iter" and "c_iter" on the FE space 
  // in order to obtain initial vector for NOX. 
  info("Projecting initial solutions on the FE meshes.");
  scalar* coeff_vec = new scalar[ndof];
  OGProjection::project_global(Hermes::vector<Space *>(t_space, c_space), Hermes::vector<MeshFunction*>(&t_prev_time_1, &c_prev_time_1),
    coeff_vec);

  // Measure the projection time.
  double proj_time = cpu_time.tick().last();

  // Initialize finite element problem.
  DiscreteProblem dp(&wf, Hermes::vector<Space*>(t_space, c_space));

  // Initialize NOX solver and preconditioner.
  NoxSolver solver(&dp);
  RCP<Precond> pc = rcp(new MlPrecond("sa"));
  if (PRECOND)
  {
    if (JFNK) solver.set_precond(pc);
    else solver.set_precond("Ifpack");
  }
  if (TRILINOS_OUTPUT)
    solver.set_output_flags(NOX::Utils::Error | NOX::Utils::OuterIteration |
                            NOX::Utils::OuterIterationStatusTest |
                            NOX::Utils::LinearSolverDetails);

  // Time stepping loop:
  double total_time = 0.0;
  cpu_time.tick_reset();
  for (int ts = 1; total_time <= T_FINAL; ts++)
  {
//.........这里部分代码省略.........

bool Adapt::adapt(Tuple<RefinementSelectors::Selector *> refinement_selectors, double thr, int strat, 
            int regularize, double to_be_processed)
{
  error_if(!have_errors, "element errors have to be calculated first, call calc_elem_errors().");
  error_if(refinement_selectors == Tuple<RefinementSelectors::Selector *>(), "selector not provided");
  if (spaces.size() != refinement_selectors.size()) error("Wrong number of refinement selectors.");
  TimePeriod cpu_time;

  //get meshes
  int max_id = -1;
  Mesh* meshes[H2D_MAX_COMPONENTS];
  for (int j = 0; j < this->neq; j++) {
    meshes[j] = this->spaces[j]->get_mesh();
    rsln[j]->set_quad_2d(&g_quad_2d_std);
    rsln[j]->enable_transform(false);
    if (meshes[j]->get_max_element_id() > max_id)
      max_id = meshes[j]->get_max_element_id();
  }

  //reset element refinement info
  AUTOLA2_OR(int, idx, max_id + 1, this->neq + 1);
  for(int j = 0; j < max_id; j++)
    for(int l = 0; l < this->neq; l++)
      idx[j][l] = -1; // element not refined

  double err0_squared = 1000.0;
  double processed_error_squared = 0.0;

  vector<ElementToRefine> elem_inx_to_proc; //list of indices of elements that are going to be processed
  elem_inx_to_proc.reserve(num_act_elems);

  //adaptivity loop
  double error_squared_threshod = -1; //an error threshold that breaks the adaptivity loop in a case of strategy 1
  int num_exam_elem = 0; //a number of examined elements
  int num_ignored_elem = 0; //a number of ignored elements
  int num_not_changed = 0; //a number of element that were not changed
  int num_priority_elem = 0; //a number of elements that were processed using priority queue

  bool first_regular_element = true; //true if first regular element was not processed yet
  int inx_regular_element = 0;
  while (inx_regular_element < num_act_elems || !priority_queue.empty())
  {
    int id, comp, inx_element;

    //get element identification
    if (priority_queue.empty()) {
      id = regular_queue[inx_regular_element].id;
      comp = regular_queue[inx_regular_element].comp;
      inx_element = inx_regular_element;
      inx_regular_element++;
    }
    else {
      id = priority_queue.front().id;
      comp = priority_queue.front().comp;
      inx_element = -1;
      priority_queue.pop();
      num_priority_elem++;
    }
    num_exam_elem++;

    //get info linked with the element
    double err_squared = errors_squared[comp][id];
    Mesh* mesh = meshes[comp];
    Element* e = mesh->get_element(id);

    if (!should_ignore_element(inx_element, mesh, e)) {
      //check if adaptivity loop should end
      if (inx_element >= 0) {
        //prepare error threshold for strategy 1
        if (first_regular_element) {
          error_squared_threshod = thr * err_squared;
          first_regular_element = false;
        }

        // first refinement strategy:
        // refine elements until prescribed amount of error is processed
        // if more elements have similar error refine all to keep the mesh symmetric
        if ((strat == 0) && (processed_error_squared > sqrt(thr) * errors_squared_sum) 
                         && fabs((err_squared - err0_squared)/err0_squared) > 1e-3) break;

        // second refinement strategy:
        // refine all elements whose error is bigger than some portion of maximal error
        if ((strat == 1) && (err_squared < error_squared_threshod)) break;

        if ((strat == 2) && (err_squared < thr)) break;

        if ((strat == 3) &&
          ( (err_squared < error_squared_threshod) ||
          ( processed_error_squared > 1.5 * to_be_processed )) ) break;
      }

      // get refinement suggestion
      ElementToRefine elem_ref(id, comp);
      int current = this->spaces[comp]->get_element_order(id);
      bool refined = refinement_selectors[comp]->select_refinement(e, current, rsln[comp], elem_ref);

      //add to a list of elements that are going to be refined
      if (can_refine_element(mesh, e, refined, elem_ref) ) {
        idx[id][comp] = (int)elem_inx_to_proc.size();
        elem_inx_to_proc.push_back(elem_ref);
//.........这里部分代码省略.........

void SubDomain::apply_markers(std::map<std::size_t, std::size_t>& sub_domains,
                              std::size_t dim,
                              T sub_domain,
                              const Mesh& mesh,
                              bool check_midpoint) const
{
  // FIXME: This function can probably be folded into the above
  //        function operator[] in std::map and MeshFunction.

  log(TRACE, "Computing sub domain markers for sub domain %d.", sub_domain);

  // Compute facet - cell connectivity if necessary
  const std::size_t D = mesh.topology().dim();
  if (dim == D - 1)
  {
    mesh.init(D - 1);
    mesh.init(D - 1, D);
  }

  // Set geometric dimension (needed for SWIG interface)
  _geometric_dimension = mesh.geometry().dim();

  // Speed up the computation by only checking each vertex once (or
  // twice if it is on the boundary for some but not all facets).
  RangedIndexSet boundary_visited(mesh.num_vertices());
  RangedIndexSet interior_visited(mesh.num_vertices());
  std::vector<bool> boundary_inside(mesh.num_vertices());
  std::vector<bool> interior_inside(mesh.num_vertices());

  // Always false when not marking facets
  bool on_boundary = false;

  // Compute sub domain markers
  Progress p("Computing sub domain markers", mesh.num_entities(dim));
  for (MeshEntityIterator entity(mesh, dim); !entity.end(); ++entity)
  {
    // Check if entity is on the boundary if entity is a facet
    if (dim == D - 1)
      on_boundary = (entity->num_global_entities(D) == 1);

    // Select the visited-cache to use for this facet (or entity)
    RangedIndexSet&    is_visited = (on_boundary ? boundary_visited : interior_visited);
    std::vector<bool>& is_inside  = (on_boundary ? boundary_inside  : interior_inside);

    // Start by assuming all points are inside
    bool all_points_inside = true;

    // Check all incident vertices if dimension is > 0 (not a vertex)
    if (entity->dim() > 0)
    {
      for (VertexIterator vertex(*entity); !vertex.end(); ++vertex)
      {
        if (is_visited.insert(vertex->index()))
        {
          Array<double> x(_geometric_dimension, const_cast<double*>(vertex->x()));
          is_inside[vertex->index()] = inside(x, on_boundary);
        }

        if (!is_inside[vertex->index()])
        {
          all_points_inside = false;
          break;
        }
      }
    }

    // Check midpoint (works also in the case when we have a single vertex)
    if (all_points_inside && check_midpoint)
    {
      Array<double> x(_geometric_dimension,
                      const_cast<double*>(entity->midpoint().coordinates()));
      if (!inside(x, on_boundary))
        all_points_inside = false;
    }

    // Mark entity with all vertices inside
    if (all_points_inside)
      sub_domains[entity->index()] =  sub_domain;

    p++;
  }
}

//-----------------------------------------------------------------------------
void SubDomain::mark_facets(Mesh& mesh,
                            std::size_t sub_domain,
                            bool check_midpoint) const
{
  mark(mesh, mesh.topology().dim() - 1, sub_domain, check_midpoint);
}

// Main function.
int main(int argc, char* argv[])
{
  // Check number of command-line parameters.
  if (argc < 2) 
    error("Not enough parameters: Provide a Butcher's table type.");

  int b_type = atoi(argv[1]);
  info ("%d", b_type);

  switch (b_type)
  {
    case 1: butcher_table_type = Explicit_RK_1; break;
    case 2: butcher_table_type = Implicit_RK_1; break;

    case 3: butcher_table_type = Explicit_RK_2; break;
    case 4: butcher_table_type = Implicit_Crank_Nicolson_2_2; break;
    case 5: butcher_table_type = Implicit_SDIRK_2_2; break;
    case 6: butcher_table_type = Implicit_Lobatto_IIIA_2_2; break;
    case 7: butcher_table_type = Implicit_Lobatto_IIIB_2_2; break;
    case 8: butcher_table_type = Implicit_Lobatto_IIIC_2_2; break;

    case 9: butcher_table_type = Explicit_RK_3; break;

    case 10: butcher_table_type = Explicit_RK_4; break;
    case 11: butcher_table_type = Implicit_Lobatto_IIIA_3_4; break;
    case 12: butcher_table_type = Implicit_Lobatto_IIIB_3_4; break;
    case 13: butcher_table_type = Implicit_Lobatto_IIIC_3_4; break;

    case 14: butcher_table_type = Implicit_Radau_IIA_3_5; break;

    default: error("Admissible command-line options are from 1 to 14.");
  }

  // Choose a Butcher's table or define your own.
  ButcherTable bt(butcher_table_type);

  // Load the mesh.
  Mesh mesh;
  H2DReader 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(BDY_DIRICHLET, INIT_BDY_REF_NUM);

  // Enter boundary markers.
  BCTypes bc_types;
  bc_types.add_bc_dirichlet(BDY_DIRICHLET);

  // Enter Dirichlet boundary values.
  BCValues bc_values;
  bc_values.add_function(BDY_DIRICHLET, essential_bc_values);   

  // Create an H1 space with default shapeset.
  H1Space* space = new H1Space(&mesh, &bc_types, &bc_values, P_INIT);

  int ndof = Space::get_num_dofs(space);
  info("ndof = %d.", ndof);

  // Previous time level solution (initialized by the initial condition).
  Solution* u_prev_time = new Solution(&mesh, init_cond);

  // Initialize the weak formulation.
  WeakForm wf;
  wf.add_matrix_form(callback(stac_jacobian));
  wf.add_vector_form(callback(stac_residual));

  // Project the initial condition on the FE space to obtain initial solution coefficient vector.
  info("Projecting initial condition to translate initial condition into a vector.");
  scalar* coeff_vec = new scalar[ndof];
  OGProjection::project_global(space, u_prev_time, coeff_vec, matrix_solver);

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

  // Time stepping loop:
  double current_time = 0.0; int ts = 1;
  do 
  {
    // Perform one Runge-Kutta time step according to the selected Butcher's table.
    info("Runge-Kutta time step (t = %g, tau = %g, stages: %d).", 
         current_time, time_step, bt.get_size());
    bool verbose = true;
    if (!rk_time_step(current_time, time_step, &bt, coeff_vec, &dp, matrix_solver,
		      verbose, NEWTON_TOL, NEWTON_MAX_ITER)) {
      error("Runge-Kutta time step failed, try to decrease time step size.");
    }

    // Convert coeff_vec into a new time level solution.
    Solution::vector_to_solution(coeff_vec, space, u_prev_time);

    // Update time.
    current_time += time_step;

    // Increase counter of time steps.
    ts++;
  } 
  while (current_time < T_FINAL);
//.........这里部分代码省略.........

   /*! GMSH v2.xx mesh format export function
   *
   *  This function extends the export capabilities of
   *  Netgen to include the GMSH v2.xx File Format.
   *
   *  Current features of this function include:
   *
   *  1. Exports Triangles, Quadrangles and Tetrahedra \n
   *  2. Supports upto second order elements of each type
   *
   */
   void WriteGmsh2Format (const Mesh & mesh,
      const CSGeometry & geom,
      const string & filename)
   {
      ofstream outfile (filename.c_str());
      outfile.precision(6);
      outfile.setf (ios::fixed, ios::floatfield);
      outfile.setf (ios::showpoint);

      int np = mesh.GetNP();  /// number of points in mesh
      int ne = mesh.GetNE();  /// number of 3D elements in mesh
      int nse = mesh.GetNSE();  /// number of surface elements (BC)
      int i, j, k, l;


      /*
      * 3D section : Volume elements (currently only tetrahedra)
      */

      if ((ne > 0)
         && (mesh.VolumeElement(1).GetNP() <= 10)
         && (mesh.SurfaceElement(1).GetNP() <= 6))
      {
         cout << "Write GMSH v2.xx Format \n";
         cout << "The GMSH v2.xx export is currently available for elements upto 2nd Order\n" << endl;

         int inverttets = mparam.inverttets;
         int invertsurf = mparam.inverttrigs;

         /// Prepare GMSH 2.0 file (See GMSH 2.0 Documentation)
         outfile << "$MeshFormat\n";
         outfile << (float)2.0 << " "
            << (int)0 << " "
            << (int)sizeof(double) << "\n";
         outfile << "$EndMeshFormat\n";

         /// Write nodes
         outfile << "$Nodes\n";
         outfile << np << "\n";

         for (i = 1; i <= np; i++)
         {
            const Point3d & p = mesh.Point(i);
            outfile << i << " "; /// node number
            outfile << p.X() << " ";
            outfile << p.Y() << " ";
            outfile << p.Z() << "\n";
         }

         outfile << "$EndNodes\n";

         /// write elements (both, surface elements and volume elements)
         outfile << "$Elements\n";
         outfile << ne + nse << "\n";  ////  number of elements + number of surfaces BC

         for (i = 1; i <= nse; i++)
         {
            int elType = 0;

            Element2d el = mesh.SurfaceElement(i);
            if(invertsurf) el.Invert();

            if(el.GetNP() == 3) elType = GMSH_TRIG;	//// GMSH Type for a 3 node triangle
            if(el.GetNP() == 6) elType = GMSH_TRIG6;  //// GMSH Type for a 6 node triangle
            if(elType == 0)
            {
               cout << " Invalid surface element type for Gmsh 2.0 3D-Mesh Export Format !\n";
               return;
            }

            outfile << i;
            outfile << " ";
            outfile << elType;
            outfile << " ";
            outfile << "2";                  //// Number of tags (2 => Physical and elementary entities)
            outfile << " ";
            outfile << mesh.GetFaceDescriptor (el.GetIndex()).BCProperty() << " ";
            /// that means that physical entity = elementary entity (arbitrary approach)
            outfile << mesh.GetFaceDescriptor (el.GetIndex()).BCProperty() << " ";
            for (j = 1; j <= el.GetNP(); j++)
            {
               outfile << " ";
               outfile << el.PNum(triGmsh[j]);
            }
            outfile << "\n";
         }


         for (i = 1; i <= ne; i++)
//.........这里部分代码省略.........

 // Create a new netgen mesh object
 DLL_HEADER Ng_Mesh * Ng_NewMesh ()
 {
    Mesh * mesh = new Mesh;  
    mesh->AddFaceDescriptor (FaceDescriptor (1, 1, 0, 1));
    return (Ng_Mesh*)(void*)mesh;
 }

// *****************************************************************
void	VertexPaint::fillSelectionColor()
{
	int		mci;

	// Put Data in Undo/Redo List.
	if(!theHold.Holding())
		theHold.Begin();
	
	ModContextList	modContexts;
	INodeTab		nodeTab;
	
	GetCOREInterface()->GetModContexts(modContexts, nodeTab);
	
	for (mci=0; mci<modContexts.Count(); mci++) 
	{
		ModContext *mc = modContexts[mci];
		if(mc && mc->localData)
			theHold.Put(new VertexPaintRestore((VertexPaintData*)mc->localData, this));
	}

	theHold.Accept(GetString(IDS_RESTORE_FILL));


	// Which Component to change??
	VertexPaintData::TComponent	whichComponent;
	switch(getEditionType())
	{
	case 0: whichComponent= VertexPaintData::Red; break;
	case 1: whichComponent= VertexPaintData::Green; break;
	case 2: whichComponent= VertexPaintData::Blue; break;
	}


	// Modify all meshes.
	for (mci=0; mci<modContexts.Count(); mci++) 
	{
		ModContext *mc = modContexts[mci];
		if(mc && mc->localData)
		{
			VertexPaintData* d = (VertexPaintData*)mc->localData;
			Mesh*		mesh = 	d->GetMesh();
			if (mesh && mesh->vertCol) 
			{
				// For all faces of the mesh
				for(int fi=0; fi<mesh->getNumFaces(); fi++)
				{
					Face* f = &mesh->faces[fi];

					for (int i=0; i<3; i++) 
					{
						// Skip painting because not selected??
						if(mesh->selLevel == MESH_VERTEX && !mesh->VertSel()[f->v[i]] )
							continue;
						if(mesh->selLevel == MESH_FACE && !mesh->FaceSel()[fi])
							continue;
						// Also skip if face is hidden.
						if(f->Hidden())
							continue;

						// Apply painting
						d->SetColor(f->v[i],  1, GetActiveColor(), whichComponent);
					}
				}
			}
		}
	}

	// refresh
	NotifyDependents(FOREVER, PART_VERTCOLOR, REFMSG_CHANGE);
	ip->RedrawViews(ip->GetTime());
}

int main(int argc, char* argv[])
{
  // Load the mesh
  Mesh mesh;
  H2DReader mloader;
  mloader.load("domain.mesh", &mesh);
  // initial uniform mesh refinement
  for (int i=0; i < INIT_REF_NUM; i++) mesh.refine_all_elements();

  // Initialize the shapeset and the cache
  H1Shapeset shapeset;
  PrecalcShapeset pss(&shapeset);

  // Create finite element space
  H1Space space(&mesh, &shapeset);
  space.set_bc_types(bc_types);
  space.set_bc_values(bc_values);
  space.set_uniform_order(P_INIT);

  // Enumerate basis functions
  space.assign_dofs();

  // initialize the weak formulation
  WeakForm wf(1);
  wf.add_biform(0, 0, bilinear_form_1, bilinear_form_ord, SYM, 1);
  wf.add_biform(0, 0, bilinear_form_2, bilinear_form_ord, SYM, 2);
  wf.add_biform(0, 0, bilinear_form_3, bilinear_form_ord, SYM, 3);
  wf.add_biform(0, 0, bilinear_form_4, bilinear_form_ord, SYM, 4);
  wf.add_biform(0, 0, bilinear_form_5, bilinear_form_ord, SYM, 5);
  wf.add_liform(0, linear_form_1, linear_form_ord, 1);
  wf.add_liform(0, linear_form_3, linear_form_ord, 3);

  // Visualize solution and mesh
  ScalarView sview("Coarse solution", 0, 100, 798, 700);
  OrderView  oview("Polynomial orders", 800, 100, 798, 700);

  // Matrix solver
  UmfpackSolver solver;

  // DOF and CPU convergence graphs
  SimpleGraph graph_dof, graph_cpu;

  // prepare selector
  RefinementSelectors::H1NonUniformHP selector(ISO_ONLY, ADAPT_TYPE, 1.0, H2DRS_DEFAULT_ORDER, &shapeset);

  // Adaptivity loop
  int it = 1, ndofs;
  bool done = false;
  double cpu = 0.0;
  Solution sln_coarse, sln_fine;
  do
    {
    info("\n---- Adaptivity step %d ---------------------------------------------\n", it++);

    // Time measurement
    begin_time();

    // Solve the coarse mesh problem
    LinSystem ls(&wf, &solver);
    ls.set_spaces(1, &space);
    ls.set_pss(1, &pss);
    ls.assemble();
    ls.solve(1, &sln_coarse);

    // Time measurement
    cpu += end_time();

    // View the solution and mesh
    sview.show(&sln_coarse);
    oview.show(&space);

    // Time measurement
    begin_time();

    // Solve the fine mesh problem
    RefSystem rs(&ls);
    rs.assemble();
    rs.solve(1, &sln_fine);

    // Calculate error estimate wrt. fine mesh solution
    H1AdaptHP hp(1, &space);
    double err_est = hp.calc_error(&sln_coarse, &sln_fine) * 100;
    info("Estimate of error: %g%%", err_est);

    // add entry to DOF convergence graph
    graph_dof.add_values(space.get_num_dofs(), err_est);
    graph_dof.save("conv_dof.dat");

    // add entry to CPU convergence graph
    graph_cpu.add_values(cpu, err_est);
    graph_cpu.save("conv_cpu.dat");

    // If err_est too large, adapt the mesh
    if (err_est < ERR_STOP) done = true;
    else {
      hp.adapt(THRESHOLD, STRATEGY, &selector, MESH_REGULARITY);
      ndofs = space.assign_dofs();
      if (ndofs >= NDOF_STOP) done = true;
    }

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

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

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

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

  // Convert initial condition into a Solution<std::complex<double> >.
  CustomInitialCondition psi_time_prev(&mesh);
  Solution<std::complex<double> > psi_time_new(&mesh);

  // Initialize the weak formulation.
  double current_time = 0;

  CustomWeakFormGPRK wf(h, m, g, omega);
  
  // Initialize boundary conditions.
  DefaultEssentialBCConst<std::complex<double> > bc_essential("Bdy", 0.0);
  EssentialBCs<std::complex<double> > bcs(&bc_essential);

  // Create an H1 space with default shapeset.
  H1Space<std::complex<double> > space(&mesh, &bcs, P_INIT);
  int ndof = space.get_num_dofs();
  Hermes::Mixins::Loggable::Static::info("ndof = %d", ndof);
 
  // Initialize the FE problem.
  DiscreteProblem<std::complex<double> > dp(&wf, &space);

  // Initialize views.
  ScalarView sview_real("Solution - real part", new WinGeom(0, 0, 600, 500));
  ScalarView sview_imag("Solution - imaginary part", new WinGeom(610, 0, 600, 500));
  sview_real.fix_scale_width(80);
  sview_imag.fix_scale_width(80);

  // Initialize Runge-Kutta time stepping.
  RungeKutta<std::complex<double> > runge_kutta(&wf, &space, &bt);
  
  // Time stepping:
  int ts = 1;
  int nstep = (int)(T_FINAL/time_step + 0.5);
  for(int ts = 1; ts <= nstep; ts++)
  {
    // Perform one Runge-Kutta time step according to the selected Butcher's table.
    Hermes::Mixins::Loggable::Static::info("Runge-Kutta time step (t = %g s, time step = %g s, stages: %d).", 
         current_time, time_step, bt.get_size());
    
    try
    {
      runge_kutta.setTime(current_time);
      runge_kutta.setTimeStep(time_step);
      runge_kutta.rk_time_step_newton(&psi_time_prev, &psi_time_new);
    }
    catch(Exceptions::Exception& e)
    {
      e.printMsg();
      throw Hermes::Exceptions::Exception("Runge-Kutta time step failed");
    }

    // Show the new time level solution.
    char title[100];
    sprintf(title, "Solution - real part, Time %3.2f s", current_time);
    sview_real.set_title(title);
    sprintf(title, "Solution - imaginary part, Time %3.2f s", current_time);
    sview_imag.set_title(title);
    RealFilter real(&psi_time_new);
    ImagFilter imag(&psi_time_new);
    sview_real.show(&real);
    sview_imag.show(&imag);

    // Copy solution for the new time step.
    psi_time_prev.copy(&psi_time_new);

    // Increase current time and time step counter.
    current_time += time_step;
    ts++;
  }

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