C++ Projection::get_error方法代码示例

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

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

  if (strcmp(args[1], "h1") != 0 && strcmp(args[1], "h1-ipol"))
    error("Unknown type of the projection.");

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

	// Refine the mesh.
  mesh.refine_all_elements(H3D_H3D_H3D_REFT_HEX_XYZ);

	// Initialize the space.
#if defined X2_Y2_Z2
  Ord3 order(2, 2, 2);
#elif defined X3_Y3_Z3
  Ord3 order(3, 3, 3);
#elif defined XN_YM_ZO
  Ord3 order(2, 3, 4);
#endif
  H1Space space(&mesh, bc_types, essential_bc_values, order);

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

  // 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.
  dp.assemble(matrix, rhs);

  // Solve the linear system. If successful, obtain the solution.
  info("Solving the linear problem.");
  Solution sln(&mesh);
  if(solver->solve()) Solution::vector_to_solution(solver->get_solution(), &space, &sln);
  else {
	  info("Matrix solver failed.");
	  success_test = 0;
  }

  unsigned int ne = mesh.get_num_base_elements();
  for(std::map<unsigned int, Element*>::iterator it = mesh.elements.begin(); it != mesh.elements.end(); it++) {
    // We are done with base elements.
    if(it->first > ne)
      break;
    Element *e = it->second;
    
    Ord3 order(4, 4, 4);
    double error;

    Projection *proj;
    if (strcmp(args[1], "h1") == 0) proj = new H1Projection(&sln, e, space.get_shapeset());
    else if (strcmp(args[1], "h1-ipol") == 0) proj = new H1ProjectionIpol(&sln, e, space.get_shapeset());
    else success_test = 0;

    error = 0.0;
    error += proj->get_error(H3D_REFT_HEX_NONE, -1, order);
    error = sqrt(error);
    
    if (error > EPS)
		    // Calculated solution is not precise enough.
		    success_test = 0;

    //
    error = 0.0;
    error += proj->get_error(H3D_REFT_HEX_X, 20, order);
    error += proj->get_error(H3D_REFT_HEX_X, 21, order);
    error = sqrt(error);
    if (error > EPS)
		    // Calculated solution is not precise enough.
		    success_test = 0;

    //
    error = 0.0;
    error += proj->get_error(H3D_REFT_HEX_Y, 22, order);
    error += proj->get_error(H3D_REFT_HEX_Y, 23, order);
    error = sqrt(error);
    if (error > EPS)
//.........这里部分代码省略.........

int main(int argc, char *argv[])
{
	_F_
	int ret = ERROR_SUCCESS;

	if (argc < 3) {
		fprintf(stderr, "ERROR: not enough parameters\n");
		return ERR_FAILURE;
	}

	if (strcmp(argv[1], "h1") != 0 && strcmp(argv[1], "h1-ipol")) {
		fprintf(stderr, "ERROR: unknown type of the projection\n");
		return ERR_FAILURE;
	}

#ifdef WITH_PETSC
	PetscInitialize(NULL, NULL, (char *) PETSC_NULL, PETSC_NULL);
#endif
	set_verbose(false);

	H1ShapesetLobattoHex shapeset;

	Mesh mesh;
	Mesh3DReader mloader;
	if (!mloader.load(argv[2], &mesh)) {
		fprintf(stderr, "ERROR: loading mesh file '%s'\n", argv[2]);
		return ERR_FAILURE;
	}

#if defined WITH_UMFPACK
	UMFPackMatrix mat;
	UMFPackVector rhs;
	UMFPackLinearSolver solver(&mat, &rhs);
#elif defined WITH_PARDISO
	PardisoMatrix mat;
	PardisoVector rhs;
	PardisoLinearSolver solver(&mat, &rhs);
#elif defined WITH_PETSC
	PetscMatrix mat;
	PetscVector rhs;
	PetscLinearSolver solver(&mat, &rhs);
#elif defined WITH_MUMPS
	MumpsMatrix mat;
	MumpsVector rhs;
	MumpsSolver solver(&mat, &rhs);
#endif

	Word_t ne = mesh.elements.count();
	// make the mesh for the ref. solution
	mesh.refine_all_elements(H3D_H3D_H3D_REFT_HEX_XYZ);

	H1Space space(&mesh, &shapeset);
	space.set_bc_types(bc_types);
	space.set_essential_bc_values(essential_bc_values);

#if defined X2_Y2_Z2
	order3_t o(2, 2, 2);
#elif defined X3_Y3_Z3
	order3_t o(3, 3, 3);
#elif defined XN_YM_ZO
	order3_t o(2, 3, 4);
#endif
	space.set_uniform_order(o);

	WeakForm wf;
	wf.add_matrix_form(bilinear_form<double, scalar>, bilinear_form<ord_t, ord_t>, SYM, ANY);
	wf.add_vector_form(linear_form<double, scalar>, linear_form<ord_t, ord_t>, ANY);

	LinearProblem lp(&wf, &space);

	space.assign_dofs();

	// assemble the stiffness matrix
	lp.assemble(&mat, &rhs);

	// solve the stiffness matrix
	solver.solve();

	Solution sln(&mesh);
	sln.set_coeff_vector(&space, solver.get_solution());

	for (Word_t idx = mesh.elements.first(); idx <= ne; idx = mesh.elements.next(idx)) {
		Element *e = mesh.elements[idx];

		order3_t order(4, 4, 4);
		double error;

		Projection *proj;
		if (strcmp(argv[1], "h1") == 0) proj = new H1Projection(&sln, e, &shapeset);
		else if (strcmp(argv[1], "h1-ipol") == 0) proj = new H1ProjectionIpol(&sln, e, &shapeset);
		else return ERR_FAILURE;

		//
		error = 0.0;
		error += proj->get_error(H3D_REFT_HEX_NONE, -1, order);
		error = sqrt(error);
		CHECK_ERROR;

		//
		error = 0.0;
//.........这里部分代码省略.........