C++ MyMesh::all_second_order方法代码示例

// Begin main function
int main (int argc, char** argv){
//PerfLog P("Example 2 Program");
//P.push("init","Program Initilization");

	// Gather command-line options
	UserOptions user = gather_command_line(argc, argv);

    // Initialize libraries
    LibMeshInit init (argc, argv);
        
	// Initilize the mesh and order variables
	MyMesh mesh;
	Order order;
	
	// Patch test (Bhatti Example 8-1)
	if(user.get_flag("patch")){
		GMVIO(mesh).read("../data/fem/examples/input/example2.gmv");
		mesh.all_first_order();
		order = FIRST;
	
	// 2-D box
	} else if (user.get_flag("2D")){
		MeshTools::Generation::build_square(mesh, 10, 10, 0., 0.04, 0., 0.04, TRI6);
		mesh.all_second_order();
		order = SECOND;
		
	// 3-D box	
	} else if (user.get_flag("3D")){
		MeshTools::Generation::build_cube(mesh, 10, 10, 10, 0., 0.04, 0., 0.04, 0., 0.04, TET10);
		mesh.all_second_order();
		order = SECOND;
	
	// Multi-element implementation of Bhatti Exmample 8-1
	} else {
		mesh.set_mesh_dimension(2);
		mesh.add_point(Point(0,0));
		mesh.add_point(Point(0.02,0));
		mesh.add_point(Point(0.02,0.04));
		mesh.add_point(Point(0,0.02));
	
		TriangleInterface t(mesh);
		t.desired_area() = 1e-4;
		t.triangulation_type() = TriangleInterface::PSLG;
		t.smooth_after_generating() = true;
		t.triangulate();
		mesh.all_second_order();
		order = SECOND;
	}
	
	// Create an equation system
 	EquationSystems eq_sys(mesh); 
 
	// Create a HeatEq class
	boost::shared_ptr<HeatEq> heateq(new HeatEq(eq_sys, order));
	 	
 	// Define the material constants
	heateq->system().set_constant<Real>("k", user.get<Real>("conductivity"));
    heateq->system().set_constant<Real>("rho", user.get<Real>("density"));
    heateq->system().set_constant<Real>("cp", user.get<Real>("specific-heat"));

	// Link to the initialization function
	heateq->system().add_initial_function(initial_function);

    // Add boundary IDs (this is some custom functionality that I added)
    mesh.find_neighbors();
    mesh.boundary_info->clear();
    mesh.add_boundary_id(0, "y", 0.0);  // bottom
    mesh.add_boundary_id(1, "x", 0.02); // right
    mesh.add_boundary_id(2, "x", 0.0);  // left
    mesh.add_boundary_id(3); // top

	// Convection boundary at bottom (user-specified)
	pConvection pC = heateq->system().add_boundary<HeatEqBoundaryConvection>(0);
	pC->h_constant = user.get<Real>("h-coefficient");
	pC->Tinf_constant = user.get<Real>("Tinf");	

	// Flux boundary at right-side (user-specified)
	pNeumann pN = heateq->system().add_boundary<HeatEqBoundaryNeumann>(1);
	pN->q_constant = user.get<Real>("flux");
	
	// Flux boundary at left-side (symetry; defaults to q = 0)
	heateq->system().add_boundary<HeatEqBoundaryNeumann>(2);
	
	// Top constant temperature boundary
	pDirichlet pD = heateq->system().add_boundary<HeatEqBoundaryDirichlet>(3);
	pD->fptr = dirichlet_function; // links the boundary function

	// Initialize system
	heateq->system().init(0.0);

	// Define a general filename
    FileParts outfile("../data/fem/examples/output/example2.ex2");
       
    // Export the initial mesh
	//ExodusII_IO(mesh).write_equation_systems(outfile.add_tstep(0,3,"_"), eq_sys);
	MyVTKIO vtk("../data/fem/examples/output/example2.vtu", eq_sys);
	vtk.write(0.0);

	// Define time stepping variables
//.........这里部分代码省略.........

// Begin main function
int main (int argc, char** argv){
  
 /* These fail in libMesh's debug mode
 
	// Add command-line options
	UserOptions opt("Program options");
	opt.add_flag("help","List the available options");
	opt.add_option<int>("nx",10, "Number of elements in x direction");
	opt.add_option<int>("ny",10, "Number of elements in y direction");
	opt.add_option<int>("num-steps,N", 100, "Number of time steps");
	opt.add_option<double>("t-step,t", 0.01, "Time step (sec.)");
	opt.add_option<double>("t-start", 0, "Start time (sec.)");
	opt.apply_options(argc, argv);

	// Collect a few options
	int nx = opt.get<int>("nx");
	int ny = opt.get<int>("ny");
	int N  = opt.get<int>("num-steps");
	double dt = opt.get<double>("t-step");
	double tstart = opt.get<double>("t-start");
 */
 	// Debug version (libmesh debug build doesn't like boost::program_options)
	int nx = 10;
	int ny = 10;
	int N  = 100;
	double dt = 0.01;
	double tstart = 0;

    // Initialize libraries
    LibMeshInit init (argc, argv);
        
    // Create a mesh 
    MyMesh mesh;
	MeshTools::Generation::build_square(mesh, nx, ny, 0., 1., 0., 1., QUAD8);
	mesh.all_second_order();

    // Create a HeatEq object
    EquationSystems eq_sys(mesh); 
    boost::shared_ptr<HeatEq> heateq(new HeatEq(eq_sys));

	// Assign an initilization function
	heateq->system().add_initial_function(initial_function);

	// Define a parameters for this problem
	const double pi = boost::math::constants::pi<double>();  
    heateq->system().set_constant<Real>("k", 1 / (2 * pi * pi));
   
    // Add boundary IDs 
    mesh.find_neighbors();
    mesh.boundary_info->clear();
    mesh.add_boundary_id(0); // all boundaries

    // Dirichlet Boundary uusing the standard function pointer
	boost::shared_ptr<HeatEqBoundaryDirichlet> pBC0_func = 
		heateq->system().add_boundary<HeatEqBoundaryDirichlet>(0);
	pBC0_func->fptr = boundary_function;

	// Dirichlet Boundary using boost::function
	//boost::shared_ptr<HeatEqBoundaryDirichlet> pBC0 = eqObj.add_boundary<HeatEqBoundaryDirichlet>(1);
	//pBC0->T_constant = 1;
	
	// Method of initializing with a boost::function (overrides above)
	boost::function<void (DenseVector<Number>& output, const Point&, const Real)> init_ptr;
	init_ptr = exact_solution_vec;
	heateq->system().add_initial_function(init_ptr);
	
	// Initialize the equation system
	double time = tstart; 
	heateq->system().init(time);
	
	// Define a general filename
    FileParts outfile("../data/fem/examples/output/example1.ex2");
  
    // Export the initial mesh
	ExodusII_IO(mesh).write_equation_systems(outfile.add_tstep(0,3,"_"), eq_sys);
	
	// Loop through time
	for (int t_step = 0; t_step < N; t_step++){
		
		// Incremenet the time counter, set the time and the
		// time step size as parameters in the EquationSystem.
		time += dt;
	
		// Display a progress message
		printf("time = %f; step %d of %d\n", time, t_step, N);

		// Update the old solution vector
		heateq->system().update_solution(time, dt);

		// Assemble and solve the linear system
		heateq->system().solve();
 
		// Output evey 10 timesteps to file.
		if ( (t_step+1)%10 == 0){
			std::string out = outfile.add_tstep(t_step+1,3,"_");
			ExodusII_IO(mesh).write_equation_systems(out, eq_sys);
        }
    }

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