本文整理汇总了C++中create_linear_solver函数的典型用法代码示例。如果您正苦于以下问题:C++ create_linear_solver函数的具体用法?C++ create_linear_solver怎么用?C++ create_linear_solver使用的例子?那么恭喜您, 这里精选的函数代码示例或许可以为您提供帮助。
在下文中一共展示了create_linear_solver函数的20个代码示例,这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞,您的评价将有助于我们的系统推荐出更棒的C++代码示例。
示例1: main
int main(int argc, char* argv[])
{
Hermes2D hermes_2D;
// Load the mesh.
Mesh mesh;
H2DReader mloader;
mloader.load("domain-excentric.mesh", &mesh);
//mloader.load("domain-concentric.mesh", &mesh);
// Initial mesh refinements.
for (int i=0; i < INIT_REF_NUM; i++) mesh.refine_all_elements();
mesh.refine_towards_boundary("Inner", INIT_BDY_REF_NUM_INNER, false); // true for anisotropic refinements
mesh.refine_towards_boundary("Outer", INIT_BDY_REF_NUM_OUTER, false); // false for isotropic refinements
// Initialize boundary conditions.
EssentialBCNonConstX bc_inner_vel_x(std::string("Inner"), VEL, STARTUP_TIME);
EssentialBCNonConstY bc_inner_vel_y(std::string("Inner"), VEL, STARTUP_TIME);
DefaultEssentialBCConst bc_outer_vel(std::string("Outer"), 0.0);
EssentialBCs bcs_vel_x(Hermes::vector<EssentialBoundaryCondition *>(&bc_inner_vel_x, &bc_outer_vel));
EssentialBCs bcs_vel_y(Hermes::vector<EssentialBoundaryCondition *>(&bc_inner_vel_y, &bc_outer_vel));
EssentialBCs bcs_pressure;
// Spaces for velocity components and pressure.
H1Space xvel_space(&mesh, &bcs_vel_x, P_INIT_VEL);
H1Space yvel_space(&mesh, &bcs_vel_y, P_INIT_VEL);
#ifdef PRESSURE_IN_L2
L2Space p_space(&mesh, &bcs_pressure, P_INIT_PRESSURE);
#else
H1Space p_space(&mesh, &bcs_pressure, P_INIT_PRESSURE);
#endif
Hermes::vector<Space *> spaces = Hermes::vector<Space *>(&xvel_space, &yvel_space, &p_space);
// Calculate and report the number of degrees of freedom.
int ndof = Space::get_num_dofs(spaces);
info("ndof = %d.", ndof);
// Define projection norms.
ProjNormType vel_proj_norm = HERMES_H1_NORM;
#ifdef PRESSURE_IN_L2
ProjNormType p_proj_norm = HERMES_L2_NORM;
#else
ProjNormType p_proj_norm = HERMES_H1_NORM;
#endif
// Solutions for the Newton's iteration and time stepping.
info("Setting initial conditions.");
Solution xvel_prev_time, yvel_prev_time, p_prev_time;
xvel_prev_time.set_zero(&mesh);
yvel_prev_time.set_zero(&mesh);
p_prev_time.set_zero(&mesh);
Hermes::vector<Solution*> slns = Hermes::vector<Solution*>(&xvel_prev_time, &yvel_prev_time,
&p_prev_time);
// Initialize weak formulation.
WeakForm* wf = new WeakFormNSNewton(STOKES, RE, TAU, &xvel_prev_time, &yvel_prev_time);
// Initialize the FE problem.
DiscreteProblem dp(wf, spaces);
// 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.
VectorView vview("velocity [m/s]", new WinGeom(0, 0, 600, 500));
ScalarView pview("pressure [Pa]", new WinGeom(610, 0, 600, 500));
//vview.set_min_max_range(0, 1.6);
vview.fix_scale_width(80);
//pview.set_min_max_range(-0.9, 1.0);
pview.fix_scale_width(80);
pview.show_mesh(true);
// Project the initial condition on the FE space to obtain initial
// coefficient vector for the Newton's method.
scalar* coeff_vec = new scalar[Space::get_num_dofs(spaces)];
// Newton's vector is set to zero (no OG projection needed).
memset(coeff_vec, 0, ndof * sizeof(double));
/*
// This can be used for more complicated initial conditions.
info("Projecting initial condition to obtain initial vector for the Newton's method.");
OGProjection::project_global(spaces, slns, coeff_vec, matrix_solver,
Hermes::vector<ProjNormType>(vel_proj_norm, vel_proj_norm, p_proj_norm));
*/
// Time-stepping loop:
char title[100];
int num_time_steps = T_FINAL / TAU;
for (int ts = 1; ts <= num_time_steps; ts++)
{
current_time += TAU;
info("---- Time step %d, time = %g:", ts, current_time);
// Update time-dependent essential BCs.
info("Updating time-dependent essential BC.");
Space::update_essential_bc_values(Hermes::vector<Space *>(&xvel_space, &yvel_space, &p_space), current_time);
// Perform Newton's iteration.
info("Solving nonlinear problem:");
//.........这里部分代码省略.........
开发者ID:Zhonghua,项目名称:hermes-legacy,代码行数:101,代码来源:main.cpp
示例2: main
int main(int argc, char* argv[])
{
// Instantiate a class with global functions.
Hermes2D hermes2d;
// 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();
// Initialize boundary conditions
CustomEssentialBCNonConst bc_essential("Boundary horizontal");
EssentialBCs bcs(&bc_essential);
// Create an H1 space with default shapeset.
H1Space space(&mesh, &bcs, P_INIT);
// Initialize the weak formulation.
CustomWeakFormGeneral wf("Boundary vertical");
// Testing n_dof and correctness of solution vector
// for p_init = 1, 2, ..., 10
int success = 1;
for (int p_init = 1; p_init <= 10; p_init++) {
printf("********* p_init = %d *********\n", p_init);
space.set_uniform_order(p_init);
int ndof = space.get_num_dofs();
info("ndof = %d", ndof);
// Initialize the FE problem.
DiscreteProblem dp(&wf, &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.
if (!hermes2d.solve_newton(coeff_vec, &dp, solver, matrix, rhs)) error("Newton's iteration failed.");
// Translate the resulting coefficient vector into the Solution sln.
Solution sln;
Solution::vector_to_solution(coeff_vec, &space, &sln);
double sum = 0;
for (int i=0; i < ndof; i++) sum += coeff_vec[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.72173) > 1e-2) success = 0;
if (p_init == 2 && fabs(sum - 0.639908) > 1e-2) success = 0;
if (p_init == 3 && fabs(sum - 0.826367) > 1e-2) success = 0;
if (p_init == 4 && fabs(sum - 0.629395) > 1e-2) success = 0;
if (p_init == 5 && fabs(sum - 0.574235) > 1e-2) success = 0;
if (p_init == 6 && fabs(sum - 0.62792) > 1e-2) success = 0;
if (p_init == 7 && fabs(sum - 0.701982) > 1e-2) success = 0;
if (p_init == 8 && fabs(sum - 0.7982) > 1e-2) success = 0;
if (p_init == 9 && fabs(sum - 0.895069) > 1e-2) success = 0;
if (p_init == 10 && fabs(sum - 1.03031) > 1e-2) success = 0;
delete [] coeff_vec;
}
if (success == 1) {
printf("Success!\n");
return ERR_SUCCESS;
}
else {
printf("Failure!\n");
return ERR_FAILURE;
}
}
开发者ID:Zhonghua,项目名称:hermes-legacy,代码行数:82,代码来源:main.cpp
示例3: main
int main()
{
// Time measurement.
TimePeriod cpu_time;
cpu_time.tick();
// Create coarse mesh, set Dirichlet BC, enumerate basis functions.
Space* space = new Space(A, B, NELEM, DIR_BC_LEFT, DIR_BC_RIGHT, P_INIT, NEQ, NEQ);
// Enumerate basis functions, info for user.
int ndof = Space::get_num_dofs(space);
info("ndof: %d", ndof);
// Initialize the weak formulation.
WeakForm wf(2);
wf.add_matrix_form(0, 0, jacobian_0_0);
wf.add_matrix_form(0, 1, jacobian_0_1);
wf.add_matrix_form(1, 0, jacobian_1_0);
wf.add_matrix_form(1, 1, jacobian_1_1);
wf.add_vector_form(0, residual_0);
wf.add_vector_form(1, residual_1);
// Initialize the FE problem.
bool is_linear = false;
DiscreteProblem *dp_coarse = new DiscreteProblem(&wf, space, is_linear);
// Newton's loop on coarse mesh.
// Fill vector coeff_vec using dof and coeffs arrays in elements.
double *coeff_vec_coarse = new double[Space::get_num_dofs(space)];
get_coeff_vector(space, coeff_vec_coarse);
// Set up the solver, matrix, and rhs according to the solver selection.
SparseMatrix* matrix_coarse = create_matrix(matrix_solver);
Vector* rhs_coarse = create_vector(matrix_solver);
Solver* solver_coarse = create_linear_solver(matrix_solver, matrix_coarse, rhs_coarse);
int it = 1;
while (1)
{
// Obtain the number of degrees of freedom.
int ndof_coarse = Space::get_num_dofs(space);
// Assemble the Jacobian matrix and residual vector.
dp_coarse->assemble(coeff_vec_coarse, matrix_coarse, rhs_coarse);
// Calculate the l2-norm of residual vector.
double res_l2_norm = get_l2_norm(rhs_coarse);
// Info for user.
info("---- Newton iter %d, ndof %d, res. l2 norm %g", it, Space::get_num_dofs(space), res_l2_norm);
// If l2 norm of the residual vector is within tolerance, then quit.
// NOTE: at least one full iteration forced
// here because sometimes the initial
// residual on fine mesh is too small.
if(res_l2_norm < NEWTON_TOL_COARSE && it > 1) break;
// Multiply the residual vector with -1 since the matrix
// equation reads J(Y^n) \deltaY^{n+1} = -F(Y^n).
for(int i=0; i<ndof_coarse; i++) rhs_coarse->set(i, -rhs_coarse->get(i));
// Solve the linear system.
if(!solver_coarse->solve())
error ("Matrix solver failed.\n");
// Add \deltaY^{n+1} to Y^n.
for (int i = 0; i < ndof_coarse; i++) coeff_vec_coarse[i] += solver_coarse->get_solution()[i];
// If the maximum number of iteration has been reached, then quit.
if (it >= NEWTON_MAX_ITER) error ("Newton method did not converge.");
// Copy coefficients from vector y to elements.
set_coeff_vector(coeff_vec_coarse, space);
it++;
}
// Cleanup.
delete matrix_coarse;
delete rhs_coarse;
delete solver_coarse;
delete [] coeff_vec_coarse;
delete dp_coarse;
// DOF and CPU convergence graphs.
SimpleGraph graph_dof_est, graph_cpu_est;
SimpleGraph graph_dof_exact, graph_cpu_exact;
// 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);
// Initialize the FE problem.
bool is_linear = false;
//.........这里部分代码省略.........
开发者ID:andreslsuave,项目名称:hermes,代码行数:101,代码来源:main.cpp
示例4: main
int main()
{
// Time measurement.
TimePeriod cpu_time;
cpu_time.tick();
// Create space, set Dirichlet BC, enumerate basis functions.
Space* space = new Space(A, B, NELEM, DIR_BC_LEFT, DIR_BC_RIGHT, P_INIT, NEQ, NEQ);
// Enumerate basis functions, info for user.
int ndof = Space::get_num_dofs(space);
info("ndof: %d", ndof);
// Initialize the weak formulation.
WeakForm wf(2);
wf.add_matrix_form(0, 0, jacobian_0_0);
wf.add_matrix_form(0, 1, jacobian_0_1);
wf.add_matrix_form(1, 0, jacobian_1_0);
wf.add_matrix_form(1, 1, jacobian_1_1);
wf.add_vector_form(0, residual_0);
wf.add_vector_form(1, residual_1);
// Initialize the FE problem.
bool is_linear = false;
DiscreteProblem *dp = new DiscreteProblem(&wf, space, is_linear);
// Newton's loop.
// Fill vector coeff_vec using dof and coeffs arrays in elements.
double *coeff_vec = new double[Space::get_num_dofs(space)];
get_coeff_vector(space, coeff_vec);
// 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);
int it = 1;
bool success = false;
while (1)
{
// Obtain the number of degrees of freedom.
int ndof = Space::get_num_dofs(space);
// Assemble the Jacobian matrix and residual vector.
dp->assemble(coeff_vec, matrix, rhs);
// Calculate the l2-norm of residual vector.
double res_l2_norm = get_l2_norm(rhs);
// Info for user.
info("---- Newton iter %d, ndof %d, res. l2 norm %g", it, Space::get_num_dofs(space), res_l2_norm);
// If l2 norm of the residual vector is within tolerance, then quit.
// NOTE: at least one full iteration forced
// here because sometimes the initial
// residual on fine mesh is too small.
if(res_l2_norm < NEWTON_TOL && it > 1) break;
// Multiply the residual vector with -1 since the matrix
// equation reads J(Y^n) \deltaY^{n+1} = -F(Y^n).
for(int i=0; i<ndof; i++) rhs->set(i, -rhs->get(i));
// Solve the linear system.
if(!(success = solver->solve()))
error ("Matrix solver failed.\n");
// Add \deltaY^{n+1} to Y^n.
for (int i = 0; i < ndof; i++) coeff_vec[i] += solver->get_solution()[i];
// If the maximum number of iteration has been reached, then quit.
if (it >= NEWTON_MAX_ITER) error ("Newton method did not converge.");
// Copy coefficients from vector y to elements.
set_coeff_vector(coeff_vec, space);
it++;
}
info("Total running time: %g s", cpu_time.accumulated());
// Test variable.
info("ndof = %d.", Space::get_num_dofs(space));
if (success)
{
info("Success!");
return ERROR_SUCCESS;
}
else
{
info("Failure!");
return ERROR_FAILURE;
}
}
开发者ID:alieed,项目名称:hermes,代码行数:92,代码来源:main.cpp
示例5: main
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_1, o1);
// Initialize the space 2.
Ord3 o2(4, 4, 4);
H1Space space2(&mesh, bc_types, essential_bc_values_2, 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;
}
}
开发者ID:navi-vonavi,项目名称:hermes,代码行数:89,代码来源:main.cpp
示例6: main
int main(int argc, char **args)
{
// Test variable.
int success_test = 1;
if (argc < 5) 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.
sscanf(args[2], "%d", &P_INIT_X);
sscanf(args[3], "%d", &P_INIT_Y);
sscanf(args[4], "%d", &P_INIT_Z);
Ord3 order(P_INIT_X, P_INIT_Y, P_INIT_Z);
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);
wf.add_vector_form(linear_form<double, scalar>, linear_form<Ord, Ord>, HERMES_ANY);
// Time measurement.
TimePeriod cpu_time;
cpu_time.tick();
// Initialize the solver in the case of SOLVER_PETSC or SOLVER_MUMPS.
initialize_solution_environment(matrix_solver, argc, args);
// 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);
out_orders_vtk(ref_space, "space", as);
// Initialize the FE problem.
bool is_linear = true;
DiscreteProblem lp(&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));
lp.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 {
printf("Matrix solver failed.\n");
success_test = 0;
}
// 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();
// Calculate element errors and total error estimate.
info("Calculating error estimate and exact error.");
Adapt *adaptivity = new Adapt(&space, HERMES_H1_NORM);
bool solutions_for_adapt = true;
double err_est_rel = adaptivity->calc_err_est(&sln, &ref_sln, solutions_for_adapt) * 100;
// Report results.
info("ndof_coarse: %d, ndof_fine: %d.", Space::get_num_dofs(&space), Space::get_num_dofs(ref_space));
info("err_est_rel: %g%%.", err_est_rel);
// If err_est_rel is too large, adapt the mesh.
if (err_est_rel < ERR_STOP) {
done = true;
//.........这里部分代码省略.........
开发者ID:Veix123,项目名称:hermes,代码行数:101,代码来源:main.cpp
示例7: main
int main()
{
// Create space, set Dirichlet BC, enumerate basis functions.
Space* space = new Space(A, B, NELEM, DIR_BC_LEFT, Hermes::vector<BCSpec *>(), P_INIT, NEQ, NEQ);
// Enumerate basis functions, info for user.
int ndof = Space::get_num_dofs(space);
info("ndof: %d", ndof);
// Initialize the weak formulation.
WeakForm wf(2);
wf.add_matrix_form(0, 0, jacobian_0_0);
wf.add_matrix_form(0, 1, jacobian_0_1);
wf.add_matrix_form(1, 0, jacobian_1_0);
wf.add_matrix_form(1, 1, jacobian_1_1);
wf.add_vector_form(0, residual_0);
wf.add_vector_form(1, residual_1);
// Initialize the FE problem.
bool is_linear = false;
DiscreteProblem *dp = new DiscreteProblem(&wf, space, is_linear);
// Newton's loop.
// Fill vector coeff_vec using dof and coeffs arrays in elements.
double *coeff_vec = new double[Space::get_num_dofs(space)];
get_coeff_vector(space, coeff_vec);
// 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);
int it = 1;
while (1)
{
// Obtain the number of degrees of freedom.
int ndof = Space::get_num_dofs(space);
// Assemble the Jacobian matrix and residual vector.
dp->assemble(coeff_vec, matrix, rhs);
// Calculate the l2-norm of residual vector.
double res_l2_norm = get_l2_norm(rhs);
// Info for user.
info("---- Newton iter %d, ndof %d, res. l2 norm %g", it, Space::get_num_dofs(space), res_l2_norm);
// If l2 norm of the residual vector is within tolerance, then quit.
// NOTE: at least one full iteration forced
// here because sometimes the initial
// residual on fine mesh is too small.
if(res_l2_norm < NEWTON_TOL && it > 1) break;
// Multiply the residual vector with -1 since the matrix
// equation reads J(Y^n) \deltaY^{n+1} = -F(Y^n).
for(int i=0; i<ndof; i++) rhs->set(i, -rhs->get(i));
// Solve the linear system.
if(!solver->solve())
error ("Matrix solver failed.\n");
// Add \deltaY^{n+1} to Y^n.
for (int i = 0; i < ndof; i++) coeff_vec[i] += solver->get_solution()[i];
// If the maximum number of iteration has been reached, then quit.
if (it >= NEWTON_MAX_ITER) error ("Newton method did not converge.");
// Copy coefficients from vector y to elements.
set_coeff_vector(coeff_vec, space);
it++;
}
// Plot the solution.
Linearizer l(space);
l.plot_solution("solution.gp");
// Plot the resulting space.
space->plot("space.gp");
// Cleanup.
for(unsigned i = 0; i < DIR_BC_LEFT.size(); i++)
delete DIR_BC_LEFT[i];
DIR_BC_LEFT.clear();
for(unsigned i = 0; i < DIR_BC_RIGHT.size(); i++)
delete DIR_BC_RIGHT[i];
DIR_BC_RIGHT.clear();
delete matrix;
delete rhs;
delete solver;
delete[] coeff_vec;
delete dp;
delete space;
info("Done.");
return 0;
}
开发者ID:B-Rich,项目名称:hermes-legacy,代码行数:99,代码来源:main.cpp
示例8: main
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)
//.........这里部分代码省略.........
开发者ID:Zhonghua,项目名称:hermes-dev,代码行数:101,代码来源:main.cpp
示例9: main
int main(int argc, char **args)
{
// Test variable.
int success_test = 1;
if (argc < 5) 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.
sscanf(args[2], "%d", &m);
sscanf(args[3], "%d", &n);
sscanf(args[4], "%d", &o);
int mx = maxn(4, m, n, o, 4);
Ord3 order(mx, mx, mx);
H1Space space(&mesh, bc_types, NULL, order);
// Initialize the weak formulation.
WeakForm wf;
wf.add_matrix_form(bilinear_form<double, scalar>, bilinear_form<Ord, Ord>, HERMES_SYM);
wf.add_matrix_form_surf(bilinear_form_surf<double, scalar>, bilinear_form_surf<Ord, Ord>);
wf.add_vector_form(linear_form<double, scalar>, linear_form<Ord, Ord>);
wf.add_vector_form_surf(linear_form_surf<double, scalar>, linear_form_surf<Ord, Ord>);
// 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_H1_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;
}
}
开发者ID:andreslsuave,项目名称:hermes,代码行数:82,代码来源:main.cpp
示例10: main
int main(int argc, char **args)
{
// Test variable.
int success_test = 1;
for (int i = 0; i < 48; i++) {
for (int j = 0; j < 48; j++) {
info("Config: %d, %d ", i, j);
Mesh mesh;
for (unsigned int k = 0; k < countof(vtcs); k++)
mesh.add_vertex(vtcs[k].x, vtcs[k].y, vtcs[k].z);
unsigned int h1[] = {
hexs[0][i][0] + 1, hexs[0][i][1] + 1, hexs[0][i][2] + 1, hexs[0][i][3] + 1,
hexs[0][i][4] + 1, hexs[0][i][5] + 1, hexs[0][i][6] + 1, hexs[0][i][7] + 1 };
mesh.add_hex(h1);
unsigned int h2[] = {
hexs[1][j][0] + 1, hexs[1][j][1] + 1, hexs[1][j][2] + 1, hexs[1][j][3] + 1,
hexs[1][j][4] + 1, hexs[1][j][5] + 1, hexs[1][j][6] + 1, hexs[1][j][7] + 1 };
mesh.add_hex(h2);
// bc
for (unsigned int k = 0; k < countof(bnd); k++) {
unsigned int facet_idxs[Quad::NUM_VERTICES] = { bnd[k][0] + 1, bnd[k][1] + 1, bnd[k][2] + 1, bnd[k][3] + 1 };
mesh.add_quad_boundary(facet_idxs, bnd[k][4]);
}
mesh.ugh();
// Initialize the space.
H1Space space(&mesh, bc_types, essential_bc_values);
#ifdef XM_YN_ZO
Ord3 ord(4, 4, 4);
#elif defined XM_YN_ZO_2
Ord3 ord(4, 4, 4);
#elif defined X2_Y2_Z2
Ord3 ord(2, 2, 2);
#endif
space.set_uniform_order(ord);
// Initialize the weak formulation.
WeakForm wf;
#ifdef DIRICHLET
wf.add_matrix_form(bilinear_form<double, scalar>, bilinear_form<Ord, Ord>, HERMES_SYM);
wf.add_vector_form(linear_form<double, scalar>, linear_form<Ord, Ord>);
#elif defined NEWTON
wf.add_matrix_form(bilinear_form<double, scalar>, bilinear_form<Ord, Ord>, HERMES_SYM);
wf.add_matrix_form_surf(bilinear_form_surf<double, scalar>, bilinear_form_surf<Ord, Ord>);
wf.add_vector_form(linear_form<double, scalar>, linear_form<Ord, Ord>);
wf.add_vector_form_surf(linear_form_surf<double, scalar>, linear_form_surf<Ord, Ord>);
#endif
// 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(space.get_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_H1_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;
info("failed, error:%g", err_exact);
}
else
info("passed");
// Clean up.
delete matrix;
//.........这里部分代码省略.........
开发者ID:B-Rich,项目名称:hermes-legacy,代码行数:101,代码来源:main.cpp
示例11: main
int main(int argc, char* argv[])
{
// Load the mesh.
Mesh mesh;
H2DReader mloader;
mloader.load("domain.mesh", &mesh);
// 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(Hermes::vector<int>(BDY_BOTTOM, BDY_RIGHT, BDY_TOP, BDY_LEFT));
// Enter Dirichlet boundary values.
BCValues bc_values;
bc_values.add_zero(Hermes::vector<int>(BDY_BOTTOM, BDY_RIGHT, BDY_TOP, BDY_LEFT));
// Create an H1 space.
H1Space* phi_space = new H1Space(&mesh, &bc_types, &bc_values, P_INIT);
H1Space* psi_space = new H1Space(&mesh, &bc_types, &bc_values, P_INIT);
int ndof = Space::get_num_dofs(Hermes::vector<Space *>(phi_space, psi_space));
info("ndof = %d.", ndof);
// Initialize previous time level solutions.
Solution phi_prev_time, psi_prev_time;
phi_prev_time.set_exact(&mesh, init_cond_phi);
psi_prev_time.set_exact(&mesh, init_cond_psi);
// Initialize the weak formulation.
WeakForm wf(2);
wf.add_matrix_form(0, 0, callback(biform_euler_0_0));
wf.add_matrix_form(0, 1, callback(biform_euler_0_1));
wf.add_matrix_form(1, 0, callback(biform_euler_1_0));
wf.add_matrix_form(1, 1, callback(biform_euler_1_1));
wf.add_vector_form(0, callback(liform_euler_0), HERMES_ANY, &phi_prev_time);
wf.add_vector_form(1, callback(liform_euler_1), HERMES_ANY, &psi_prev_time);
// Initialize views.
ScalarView view("Psi", new WinGeom(0, 0, 600, 500));
view.fix_scale_width(80);
// Time stepping loop:
int nstep = (int)(T_FINAL/TAU + 0.5);
for(int ts = 1; ts <= nstep; ts++)
{
info("Time step %d:", ts);
info("Solving linear system.");
// Initialize the FE problem.
bool is_linear = true;
DiscreteProblem dp(&wf, Hermes::vector<Space *>(phi_space, psi_space), is_linear);
SparseMatrix* matrix = create_matrix(matrix_solver);
Vector* rhs = create_vector(matrix_solver);
Solver* solver = create_linear_solver(matrix_solver, matrix, rhs);
// 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_solutions(solver->get_solution(), Hermes::vector<Space *>(phi_space, psi_space), Hermes::vector<Solution *>(&phi_prev_time, &psi_prev_time));
else
error ("Matrix solver failed.\n");
// Show the new time level solution.
char title[100];
sprintf(title, "Time step %d", ts);
view.set_title(title);
view.show(&psi_prev_time);
}
// Wait for all views to be closed.
View::wait();
return 0;
}
开发者ID:andreslsuave,项目名称:hermes,代码行数:80,代码来源:main.cpp
示例12: main
int main(int argc, char* argv[])
{
// 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(&mesh, &bc_types, &bc_values, P_INIT);
int ndof = Space::get_num_dofs(&space);
info("ndof = %d.", ndof);
// Solution for the previous time level.
Solution u_prev_time;
// Initialize the weak formulation.
WeakForm wf;
if (TIME_INTEGRATION == 1) {
wf.add_matrix_form(jac_euler, jac_ord, HERMES_NONSYM, HERMES_ANY, &u_prev_time);
wf.add_vector_form(res_euler, res_ord, HERMES_ANY, &u_prev_time);
}
else {
wf.add_matrix_form(jac_cranic, jac_ord, HERMES_NONSYM, HERMES_ANY, &u_prev_time);
wf.add_vector_form(res_cranic, res_ord, HERMES_ANY, &u_prev_time);
}
// Project the function init_cond() on the FE space
// to obtain initial coefficient vector for the Newton's method.
scalar* coeff_vec = new scalar[Space::get_num_dofs(&space)];
info("Projecting initial condition to obtain initial vector for the Newton's method.");
OGProjection::project_global(&space, init_cond, coeff_vec, matrix_solver);
Solution::vector_to_solution(coeff_vec, &space, &u_prev_time);
// Time stepping loop:
double current_time = 0.0;
int t_step = 1;
do {
info("---- Time step %d, t = %g s.", t_step, current_time); t_step++;
// Initialize the FE problem.
bool is_linear = false;
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);
// Perform Newton's iteration.
bool verbose = true;
if (!solve_newton(coeff_vec, &dp, solver, matrix, rhs,
NEWTON_TOL, NEWTON_MAX_ITER, verbose)) error("Newton's iteration failed.");
// Translate the resulting coefficient vector into the Solution sln.
Solution::vector_to_solution(coeff_vec, &space, &u_prev_time);
// Cleanup.
delete matrix;
delete rhs;
delete solver;
// Update time.
current_time += TAU;
} while (current_time < T_FINAL);
delete [] coeff_vec;
info("Coordinate ( 0, 0) value = %lf", u_prev_time.get_pt_value(0.0, 0.0));
info("Coordinate ( 25, 25) value = %lf", u_prev_time.get_pt_value(25.0, 25.0));
info("Coordinate ( 50, 50) value = %lf", u_prev_time.get_pt_value(50.0, 50.0));
info("Coordinate ( 75, 75) value = %lf", u_prev_time.get_pt_value(75.0, 75.0));
info("Coordinate (100, 100) value = %lf", u_prev_time.get_pt_value(100.0, 100.0));
double coor_x_y[5] = {0.0, 25.0, 50.0, 75.0, 100.0};
double value[5] = {-1000.000000, -969.316013, -836.504249, -651.433710, -1000.000000};
bool success = true;
for (int i = 0; i < 5; i++)
{
if (abs(value[i] - u_prev_time.get_pt_value(coor_x_y[i], coor_x_y[i])) > 1E-6)
success = false;
}
if (success) {
printf("Success!\n");
return ERR_SUCCESS;
}
else {
//.........这里部分代码省略.........
开发者ID:andreslsuave,项目名称:hermes,代码行数:101,代码来源:main.cpp
示例13: main
int main() {
// Three macroelements are defined above via the interfaces[] array.
// poly_orders[]... initial poly degrees of macroelements.
// material_markers[]... material markers of macroelements.
// subdivisions[]... equidistant subdivision of macroelements.
int poly_orders[N_MAT] = {P_init_inner, P_init_outer, P_init_reflector };
int material_markers[N_MAT] = {Marker_inner, Marker_outer, Marker_reflector };
int subdivisions[N_MAT] = {N_subdiv_inner, N_subdiv_outer, N_subdiv_reflector };
// Create space.
Space* space = new Space(N_MAT, interfaces, poly_orders, material_markers, subdivisions, N_GRP, N_SLN);
// Enumerate basis functions, info for user.
info("N_dof = %d", Space::get_num_dofs(space));
// Initial approximation: u = 1.
double K_EFF_old;
double init_val = 1.0;
set_vertex_dofs_constant(space, init_val, 0);
// Initialize the weak formulation.
WeakForm wf;
wf.add_matrix_form(jacobian_vol_inner, NULL, Marker_inner);
wf.add_matrix_form(jacobian_vol_outer, NULL, Marker_outer);
wf.add_matrix_form(jacobian_vol_reflector, NULL, Marker_reflector);
wf.add_vector_form(residual_vol_inner, NULL, Marker_inner);
wf.add_vector_form(residual_vol_outer, NULL, Marker_outer);
wf.add_vector_form(residual_vol_reflector, NULL, Marker_reflector);
wf.add_vector_form_surf(residual_surf_left, BOUNDARY_LEFT);
wf.add_matrix_form_surf(jacobian_surf_right, BOUNDARY_RIGHT);
wf.add_vector_form_surf(residual_surf_right, BOUNDARY_RIGHT);
// Initialize the FE problem.
bool is_linear = false;
DiscreteProblem *dp = new DiscreteProblem(&wf, space, is_linear);
// Source iteration (power method).
for (int i = 0; i < Max_SI; i++)
{
// Obtain fission source.
int current_solution = 0, previous_solution = 1;
copy_dofs(current_solution, previous_solution, space);
// Obtain the number of degrees of freedom.
int ndof = Space::get_num_dofs(space);
// Fill vector coeff_vec using dof and coeffs arrays in elements.
double *coeff_vec = new double[Space::get_num_dofs(space)];
solution_to_vector(space, coeff_vec);
// 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);
int it = 1;
while (1) {
// Obtain the number of degrees of freedom.
int ndof = Space::get_num_dofs(space);
// Assemble the Jacobian matrix and residual vector.
dp->assemble(matrix, rhs);
// Calculate the l2-norm of residual vector.
double res_norm = 0;
for(int i=0; i<ndof; i++) res_norm += rhs->get(i)*rhs->get(i);
res_norm = sqrt(res_norm);
// Info for user.
info("---- Newton iter %d, residual norm: %.15f", it, res_norm);
// If l2 norm of the residual vector is within tolerance, then quit.
// NOTE: at least one full iteration forced
// here because sometimes the initial
// residual on fine mesh is too small.
if(res_norm < NEWTON_TOL && it > 1) break;
// Multiply the residual vector with -1 since the matrix
// equation reads J(Y^n) \deltaY^{n+1} = -F(Y^n).
for(int i=0; i<ndof; i++) rhs->set(i, -rhs->get(i));
// Solve the linear system.
if(!solver->solve())
error ("Matrix solver failed.\n");
// Add \deltaY^{n+1} to Y^n.
for (int i = 0; i < ndof; i++) coeff_vec[i] += solver->get_solution()[i];
// If the maximum number of iteration has been reached, then quit.
if (it >= NEWTON_MAX_ITER) error ("Newton method did not converge.");
// Copy coefficients from vector y to elements.
vector_to_solution(coeff_vec, space);
it++;
}
// Cleanup.
delete matrix;
delete rhs;
delete solver;
//.........这里部分代码省略.........
开发者ID:kameari,项目名称:hermes,代码行数:101,代码来源:main.cpp
|
请发表评论