本文整理汇总了C++中WeakForm类的典型用法代码示例。如果您正苦于以下问题:C++ WeakForm类的具体用法?C++ WeakForm怎么用?C++ WeakForm使用的例子?那么恭喜您, 这里精选的类代码示例或许可以为您提供帮助。
在下文中一共展示了WeakForm类的20个代码示例,这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞,您的评价将有助于我们的系统推荐出更棒的C++代码示例。
示例1: 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(1, 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);
// Previous time level solution (initialized by the initial condition).
Solution u_prev_time(&mesh, init_cond);
// Initialize the weak formulation.
WeakForm wf;
wf.add_matrix_form(callback(jac), HERMES_NONSYM, HERMES_ANY);
wf.add_vector_form(callback(res), HERMES_ANY, &u_prev_time);
// Project the initial condition on the FE space to obtain initial
// coefficient vector for the Newton's method.
info("Projecting initial condition to obtain initial vector for the Newton's method.");
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);
// Set up the solver, matrix, and rhs according to the solver selection.
SparseMatrix* matrix = create_matrix(matrix_solver);
Vector* rhs = create_vector(matrix_solver);
Solver* solver = create_linear_solver(matrix_solver, matrix, rhs);
// Initialize views.
ScalarView sview("Solution", new WinGeom(0, 0, 500, 400));
OrderView oview("Mesh", new WinGeom(510, 0, 460, 400));
oview.show(&space);
// Time stepping loop:
double current_time = 0.0; int ts = 1;
do
{
info("---- Time step %d, t = %g s.", ts, current_time); ts++;
// Perform Newton's iteration.
info("Solving on coarse mesh:");
bool verbose = true;
if (!solve_newton(coeff_vec, &dp, solver, matrix, rhs,
NEWTON_TOL, NEWTON_MAX_ITER, verbose)) error("Newton's iteration failed.");
// Update previous time level solution.
Solution::vector_to_solution(coeff_vec, &space, &u_prev_time);
// Update time.
current_time += TAU;
// Show the new time level solution.
char title[100];
sprintf(title, "Solution, t = %g", current_time);
sview.set_title(title);
sview.show(&u_prev_time);
oview.show(&space);
}
while (current_time < T_FINAL);
// Cleanup.
delete [] coeff_vec;
delete matrix;
delete rhs;
delete solver;
// Wait for all views to be closed.
View::wait();
return 0;
}
开发者ID:andreslsuave,项目名称:hermes,代码行数:89,代码来源:main.cpp
示例2: main
int main(int argc, char* argv[])
{
// Time measurement
TimePeriod cpu_time;
cpu_time.tick();
// Load the mesh.
Mesh mesh;
H2DReader mloader;
mloader.load("lshape3q.mesh", &mesh); // quadrilaterals
//mloader.load("lshape3t.mesh", &mesh); // triangles
// Perform initial mesh refinemets.
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::Tuple<int>(BDY_1, BDY_6));
bc_types.add_bc_newton(Hermes::Tuple<int>(BDY_2, BDY_3, BDY_4, BDY_5));
// Enter Dirichlet boundary values.
BCValues bc_values;
bc_values.add_zero(Hermes::Tuple<int>(BDY_1, BDY_6));
// Create an Hcurl space with default shapeset.
HcurlSpace space(&mesh, &bc_types, &bc_values, P_INIT);
// Initialize the weak formulation.
WeakForm wf;
wf.add_matrix_form(callback(bilinear_form), HERMES_SYM);
wf.add_matrix_form_surf(callback(bilinear_form_surf));
wf.add_vector_form_surf(linear_form_surf, linear_form_surf_ord);
// Initialize coarse and reference mesh solutions.
Solution sln, ref_sln;
// Initialize exact solution.
ExactSolution sln_exact(&mesh, exact);
// Initialize refinement selector.
HcurlProjBasedSelector selector(CAND_LIST, CONV_EXP, H2DRS_DEFAULT_ORDER);
// DOF and CPU convergence graphs.
SimpleGraph graph_dof_est, graph_cpu_est,
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);
// Assemble the reference problem.
info("Solving on reference mesh.");
bool is_linear = true;
DiscreteProblem* dp = new DiscreteProblem(&wf, ref_space, is_linear);
SparseMatrix* matrix = create_matrix(matrix_solver);
Vector* rhs = create_vector(matrix_solver);
Solver* solver = create_linear_solver(matrix_solver, matrix, rhs);
dp->assemble(matrix, rhs);
// Time measurement.
cpu_time.tick();
// Solve the linear system of the reference problem. If successful, obtain the solution.
if(solver->solve()) Solution::vector_to_solution(solver->get_solution(), ref_space, &ref_sln);
else error ("Matrix solver failed.\n");
// Time measurement.
cpu_time.tick();
// Project the fine mesh solution onto the coarse mesh.
info("Projecting reference solution on coarse mesh.");
OGProjection::project_global(&space, &ref_sln, &sln, matrix_solver);
// Calculate element errors and total error estimate.
info("Calculating error estimate and exact error.");
Adapt* adaptivity = new Adapt(&space);
double err_est_rel = adaptivity->calc_err_est(&sln, &ref_sln) * 100;
// Calculate exact error,
bool solutions_for_adapt = false;
double err_exact_rel = adaptivity->calc_err_exact(&sln, &sln_exact, 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_exact_rel: %g%%", err_est_rel, err_exact_rel);
// Time measurement.
cpu_time.tick();
// Add entry to DOF and CPU convergence graphs.
graph_dof_est.add_values(Space::get_num_dofs(&space), err_est_rel);
graph_dof_est.save("conv_dof_est.dat");
graph_cpu_est.add_values(cpu_time.accumulated(), err_est_rel);
//.........这里部分代码省略.........
开发者ID:sriharifez,项目名称:hermes,代码行数:101,代码来源:main.cpp
示例3: main
int main(int argc, char* argv[])
{
// This is a hack I used to run the code a dozen of times when plotting convergence graphs.
if (argc > 1) {
if (argv[1][0] == 'e') method = IE;
else if (argv[1][0] == 's') method = SDIRK;
else error("what are you doing?");
}
if (argc > 2) {
TAU = std::atof(argv[2]);
}
// This is important to make sure we compare solution at exact same point in time when studying convergence.
int N_STEP = std::ceil(T_FINAL / TAU);
if (fabs(T_FINAL - N_STEP * TAU) > 1e-10) {
error("bad choice of TAU");
}
info("t_final = %g, tau = %g, n = %i", T_FINAL, TAU, N_STEP);
// 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 boudnary values.
BCValues bc_values;
bc_values.add_zero(BDY_DIRICHLET);
// 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);
// Previous time level solution (initialized by the initial condition).
Solution u_prev_time(&mesh, exact_solution);
// Project the initial condition on the FE space to obtain initial
// coefficient vector for the Newton's method.
info("Projecting initial condition to obtain initial vector for the Newton's method.");
scalar* coeff_vec = new scalar[ndof];
OGProjection::project_global(&space, &u_prev_time, coeff_vec, matrix_solver);
// Set up the solver, matrix, and rhs according to the solver selection.
SparseMatrix* matrix = create_matrix(matrix_solver);
Vector* rhs = create_vector(matrix_solver);
Solver* solver = create_linear_solver(matrix_solver, matrix, rhs);
// Initialize views.
ScalarView sview("Solution", new WinGeom(0, 0, 500, 400));
OrderView oview("Mesh", new WinGeom(520, 0, 450, 400));
oview.show(&space);
if (method == IE) {
info("IMPLICIT EULER METHOD");
// Initialize the weak formulation.
WeakForm wf;
wf.add_matrix_form(callback(jac), HERMES_NONSYM, HERMES_ANY);
wf.add_vector_form(callback(res), HERMES_ANY, &u_prev_time);
// Initialize the FE problem.
bool is_linear = false;
DiscreteProblem dp(&wf, &space, is_linear);
// Time stepping loop:
int ts = 0;
do {
info("---- Time step %d, t = %g s.", ++ts, TIME);
info("We are computing solution at next time step TIME+TAU = %g s.", TIME+TAU);
// Perform Newton's iteration.
info("Solving nonlinear problem:");
bool verbose = true;
if (!solve_newton(coeff_vec, &dp, solver, matrix, rhs,
NEWTON_TOL, NEWTON_MAX_ITER, verbose)) error("Newton's iteration failed.");
// Update previous time level solution.
Solution::vector_to_solution(coeff_vec, &space, &u_prev_time);
// Update time.
TIME = TIME + TAU;
// Compute exact error.
Solution exact_sln(&mesh, exact_solution);
double exact_l2_error = calc_abs_error(&u_prev_time, &exact_sln, HERMES_L2_NORM);
info("TIME: %g s.", TIME);
info("Exact error in l2-norm: %g.", exact_l2_error);
// Show the new time level solution.
char title[100];
//.........这里部分代码省略.........
开发者ID:andreslsuave,项目名称:hermes,代码行数:101,代码来源:main.cpp
示例4: main
int main(int argc, char* argv[])
{
// Choose a Butcher's table or define your own.
ButcherTable* bt = new ButcherTable(butcher_table_type);
if (bt->is_explicit()) info("Using a %d-stage explicit R-K method.", bt->get_size());
if (bt->is_diagonally_implicit()) info("Using a %d-stage diagonally implicit R-K method.", bt->get_size());
if (bt->is_fully_implicit()) info("Using a %d-stage fully implicit R-K method.", bt->get_size());
// Turn off adaptive time stepping if R-K method is not embedded.
if (bt->is_embedded() == false && ADAPTIVE_TIME_STEP_ON == true) {
warn("R-K method not embedded, turning off adaptive time stepping.");
ADAPTIVE_TIME_STEP_ON = false;
}
// Load the mesh.
Mesh mesh, basemesh;
H2DReader mloader;
mloader.load("square.mesh", &basemesh);
// Perform initial mesh refinements.
for(int i = 0; i < INIT_REF_NUM; i++) basemesh.refine_all_elements();
mesh.copy(&basemesh);
// 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);
// Convert initial condition into a Solution.
Solution* sln_prev_time = new Solution(&mesh, init_cond);
// Initialize the weak formulation.
WeakForm wf;
wf.add_matrix_form(callback(stac_jacobian), HERMES_NONSYM, HERMES_ANY, sln_prev_time);
wf.add_vector_form(callback(stac_residual), HERMES_ANY, sln_prev_time);
// Initialize the discrete problem.
bool is_linear = false;
DiscreteProblem dp_coarse(&wf, &space, is_linear);
// Create a refinement selector.
H1ProjBasedSelector selector(CAND_LIST, CONV_EXP, H2DRS_DEFAULT_ORDER);
// Graph for time step history.
SimpleGraph time_step_graph;
if (ADAPTIVE_TIME_STEP_ON) info("Time step history will be saved to file time_step_history.dat.");
// Time stepping loop.
double current_time = time_step; int ts = 1;
do
{
info("Begin time step %d.", ts);
// Periodic global derefinement.
if (ts > 1 && ts % UNREF_FREQ == 0)
{
info("Global mesh derefinement.");
if (UNREF_LEVEL == 1) mesh.unrefine_all_elements();
else mesh.copy(&basemesh);
space.set_uniform_order(P_INIT);
ndof = Space::get_num_dofs(&space);
}
// Spatial adaptivity loop. Note: sln_prev_time must not be
// changed during spatial adaptivity.
Solution ref_sln;
Solution* time_error_fn;
if (bt->is_embedded() == true) time_error_fn = new Solution(&mesh);
else time_error_fn = NULL;
bool done = false; int as = 1;
double err_est;
do {
// Construct globally refined reference mesh and setup reference space.
Space* ref_space = construct_refined_space(&space);
// Initialize discrete problem on reference mesh.
DiscreteProblem* ref_dp = new DiscreteProblem(&wf, ref_space);
// Runge-Kutta step on the fine mesh.
info("Runge-Kutta time step on fine mesh (t = %g s, tau = %g s, stages: %d).",
current_time, time_step, bt->get_size());
bool verbose = true;
bool is_linear = false;
if (!rk_time_step(current_time, time_step, bt, sln_prev_time, &ref_sln, time_error_fn,
ref_dp, matrix_solver, verbose, is_linear, NEWTON_TOL_FINE, NEWTON_MAX_ITER)) {
error("Runge-Kutta time step failed, try to decrease time step size.");
}
/* If ADAPTIVE_TIME_STEP_ON == true, estimate temporal error.
If too large or too small, then adjust it and restart the time step. */
double rel_err_time;
if (bt->is_embedded() == true) {
info("Calculating temporal error estimate.");
//.........这里部分代码省略.........
开发者ID:alieed,项目名称:hermes,代码行数:101,代码来源:main.cpp
示例5: main
int main(int argc, char **argv) {
int res = ERR_SUCCESS;
#ifdef WITH_PETSC
PetscInitialize(&argc, &argv, (char *) PETSC_NULL, PETSC_NULL);
#endif
set_verbose(false);
if (argc < 3) error("Not enough parameters");
HcurlShapesetLobattoHex shapeset;
printf("* Loading mesh '%s'\n", argv[1]);
Mesh mesh;
Mesh3DReader mesh_loader;
if (!mesh_loader.load(argv[1], &mesh)) error("Loading mesh file '%s'\n", argv[1]);
printf("* Setting the space up\n");
HcurlSpace space(&mesh, &shapeset);
space.set_bc_types(bc_types);
int order;
sscanf(argv[2], "%d", &order);
int dir_x = order, dir_y = order, dir_z = order;
order3_t o(dir_x, dir_y, dir_z);
printf(" - Setting uniform order to (%d, %d, %d)\n", o.x, o.y ,o.z);
space.set_uniform_order(o);
int ndofs = space.assign_dofs();
printf(" - Number of DOFs: %d\n", ndofs);
printf("* Calculating a solution\n");
#if defined WITH_UMFPACK
UMFPackMatrix mat;
UMFPackVector rhs;
UMFPackLinearSolver solver(&mat, &rhs);
#elif defined WITH_PARDISO
PardisoMatrix mat;
PardisoVector rhs;
PardisoSolver 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
WeakForm wf;
wf.add_matrix_form(bilinear_form<double, scalar>, bilinear_form<ord_t, ord_t>, SYM);
wf.add_matrix_form_surf(bilinear_form_surf<double, scalar>, bilinear_form_surf<ord_t, ord_t>);
wf.add_vector_form(linear_form<double, scalar>, linear_form<ord_t, ord_t>);
wf.add_vector_form_surf(linear_form_surf<double, scalar>, linear_form_surf<ord_t, ord_t>);
LinearProblem lp(&wf, &space);
// assemble stiffness matrix
Timer assemble_timer("Assembling stiffness matrix");
assemble_timer.start();
lp.assemble(&mat, &rhs);
assemble_timer.stop();
// solve the stiffness matrix
Timer solve_timer("Solving stiffness matrix");
solve_timer.start();
bool solved = solver.solve();
solve_timer.stop();
//#ifdef OUTPUT_DIR
mat.dump(stdout, "a");
rhs.dump(stdout, "b");
//#endif
if (solved) {
scalar *s = solver.get_solution();
Solution sln(&mesh);
sln.set_coeff_vector(&space, s);
printf("* Solution:\n");
for (int i = 1; i <= ndofs; i++) {
printf(" x[% 3d] = " SCALAR_FMT "\n", i, SCALAR(s[i]));
}
// output the measured values
printf("%s: %s (%lf secs)\n", assemble_timer.get_name(), assemble_timer.get_human_time(), assemble_timer.get_seconds());
printf("%s: %s (%lf secs)\n", solve_timer.get_name(), solve_timer.get_human_time(), solve_timer.get_seconds());
// norm
ExactSolution ex_sln(&mesh, exact_solution);
double hcurl_sln_norm = hcurl_norm(&sln);
double hcurl_err_norm = hcurl_error(&sln, &ex_sln);
printf(" - Hcurl solution norm: % le\n", hcurl_sln_norm);
printf(" - Hcurl error norm: % le\n", hcurl_err_norm);
double l2_sln_norm = l2_norm_hcurl(&sln);
double l2_err_norm = l2_error_hcurl(&sln, &ex_sln);
//.........这里部分代码省略.........
开发者ID:regmi,项目名称:hermes,代码行数:101,代码来源:hex-hcurl-cplx.cpp
示例6: main
int main(int argc, char* argv[])
{
// Load the mesh.
Mesh mesh;
ExodusIIReader mloader;
if (!mloader.load("iron-water.e", &mesh)) error("ExodusII mesh load failed.");
// Perform initial uniform mesh refinement.
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>(WATER_2, IRON));
bc_types.add_bc_neumann(WATER_1);
// Enter Dirichlet boundary values.
BCValues bc_values;
bc_values.add_zero(Hermes::vector<int>(WATER_2, IRON));
// Create an H1 space with default shapeset.
H1Space space(&mesh, &bc_types, &bc_values, P_INIT);
// Initialize the weak formulation
WeakForm wf;
wf.add_matrix_form(bilinear_form_water, bilinear_form_ord, HERMES_SYM, WATER_1);
wf.add_matrix_form(bilinear_form_water, bilinear_form_ord, HERMES_SYM, WATER_2);
wf.add_matrix_form(bilinear_form_iron, bilinear_form_ord, HERMES_SYM, IRON);
wf.add_vector_form(linear_form_source, linear_form_ord, WATER_1);
// Initialize coarse and reference mesh solution.
Solution sln, ref_sln;
// Initialize refinement selector.
H1ProjBasedSelector selector(CAND_LIST, CONV_EXP, H2DRS_DEFAULT_ORDER);
// DOF and CPU convergence graphs initialization.
SimpleGraph graph_dof, graph_cpu;
// Time measurement.
TimePeriod cpu_time;
cpu_time.tick();
// Adaptivity loop:
int as = 1;
bool done = false;
do
{
info("---- Adaptivity step %d:", as);
// Construct globally refined reference mesh and setup reference space.
Space* ref_space = construct_refined_space(&space);
// Assemble the reference problem.
info("Solving on reference mesh.");
bool is_linear = true;
DiscreteProblem* dp = new DiscreteProblem(&wf, ref_space, is_linear);
SparseMatrix* matrix = create_matrix(matrix_solver);
Vector* rhs = create_vector(matrix_solver);
Solver* solver = create_linear_solver(matrix_solver, matrix, rhs);
dp->assemble(matrix, rhs);
// Time measurement.
cpu_time.tick();
// Solve the linear system of the reference problem.
// If successful, obtain the solution.
if(solver->solve()) Solution::vector_to_solution(solver->get_solution(), ref_space, &ref_sln);
else error ("Matrix solver failed.\n");
// Project the fine mesh solution onto the coarse mesh.
info("Projecting reference solution on coarse mesh.");
OGProjection::project_global(&space, &ref_sln, &sln, matrix_solver);
// Calculate element errors and total error estimate.
info("Calculating error estimate.");
Adapt* adaptivity = new Adapt(&space);
double err_est_rel = adaptivity->calc_err_est(&sln, &ref_sln) * 100;
// Report results.
info("ndof_coarse: %d, ndof_fine: %d, err_est_rel: %g%%",
Space::get_num_dofs(&space), Space::get_num_dofs(ref_space), err_est_rel);
// Time measurement.
cpu_time.tick();
// Add entry to DOF and CPU convergence graphs.
graph_dof.add_values(Space::get_num_dofs(&space), err_est_rel);
graph_dof.save("conv_dof_est.dat");
graph_cpu.add_values(cpu_time.accumulated(), err_est_rel);
graph_cpu.save("conv_cpu_est.dat");
// If err_est too large, adapt the mesh.
if (err_est_rel < ERR_STOP) done = true;
else
{
info("Adapting coarse mesh.");
done = adaptivity->adapt(&selector, THRESHOLD, STRATEGY, MESH_REGULARITY);
// Increase the counter of performed adaptivity steps.
if (done == false) as++;
//.........这里部分代码省略.........
开发者ID:alieed,项目名称:hermes,代码行数:101,代码来源:main.cpp
示例7: main
int main(int argc, char **args)
{
// Test variable.
int success_test = 1;
// Load the initial mesh.
Mesh mesh;
H3DReader mesh_loader;
mesh_loader.load("../hexahedron.mesh3d", &mesh);
// Perform initial mesh refinement.
for (int i=0; i < INIT_REF_NUM; i++) mesh.refine_all_elements(H3D_H3D_H3D_REFT_HEX_XYZ);
// Create H1 space with default shapeset.
H1Space space(&mesh, bc_types, essential_bc_values, Ord3(P_INIT_X, P_INIT_Y, P_INIT_Z));
// Construct initial solution and set it to zero.
Solution sln_prev(&mesh);
sln_prev.set_zero();
// Initialize weak formulation.
WeakForm wf;
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>, HERMES_ANY_INT, &sln_prev);
// Initialize discrete 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).
}
// Exact error for testing purposes.
double err_exact;
// Time stepping.
int nsteps = (int) (FINAL_TIME/TAU + 0.5);
for (int ts = 0; ts < nsteps; ts++)
{
info("---- Time step %d, time %3.5f.", ts, TIME);
// Assemble the linear problem.
info("Assembling the linear problem (ndof: %d).", Space::get_num_dofs(&space));
if (ts == 0) dp.assemble(matrix, rhs);
else dp.assemble(NULL, rhs);
// Solve the linear system. If successful, obtain the solution.
info("Solving the linear problem.");
Solution sln(space.get_mesh());
if(solver->solve()) Solution::vector_to_solution(solver->get_solution(), &space, &sln);
else error ("Matrix solver failed.\n");
// Output solution.
if (solution_output)
out_fn_vtk(&sln, "sln", ts);
// Calculate exact error.
ExactSolution esln(&mesh, fndd);
info("Calculating exact error.");
Adapt *adaptivity = new Adapt(&space, HERMES_H1_NORM);
bool solutions_for_adapt = false;
err_exact = adaptivity->calc_err_exact(&sln, &esln, solutions_for_adapt, HERMES_TOTAL_ERROR_ABS) * 100;
info("Err. exact: %g%%.", err_exact);
// Next time step.
sln_prev = sln;
TIME += TAU;
// Cleanup.
delete adaptivity;
}
if(err_exact > 3.00)
success_test = 0;
// Clean up.
delete matrix;
delete rhs;
delete solver;
if (success_test) {
info("Success!");
return ERR_SUCCESS;
}
else {
info("Failure!");
return ERR_FAILURE;
//.........这里部分代码省略.........
开发者ID:B-Rich,项目名称:hermes-legacy,代码行数:101,代码来源:main.cpp
示例8: main
int main(int argc, char **args)
{
// Test variable.
int success_test = 1;
// Check the number of command-line parameters.
if (argc < 2) {
info("Use x, y, z, xy, xz, yz, or xyz as a command-line parameter.");
error("Not enough command-line parameters.");
}
// Determine anisotropy type from the command-line parameter.
ANISO_TYPE = parse_aniso_type(args[1]);
// Load the mesh.
Mesh mesh;
H3DReader mloader;
mloader.load("hex-0-1.mesh3d", &mesh);
// Assign the lowest possible directional polynomial degrees so that the problem's NDOF >= 1.
assign_poly_degrees();
// Create an H1 space with default shapeset.
info("Setting directional polynomial degrees %d, %d, %d.", P_INIT_X, P_INIT_Y, P_INIT_Z);
H1Space space(&mesh, bc_types, essential_bc_values, Ord3(P_INIT_X, P_INIT_Y, P_INIT_Z));
// Initialize weak formulation.
WeakForm wf;
wf.add_matrix_form(bilinear_form<double, scalar>, bilinear_form<Ord, Ord>, HERMES_SYM, HERMES_ANY);
wf.add_vector_form(linear_form<double, scalar>, linear_form<Ord, Ord>, HERMES_ANY);
// Set exact solution.
ExactSolution exact(&mesh, fndd);
// 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);
// Initialize the FE problem.
bool is_linear = true;
DiscreteProblem dp(&wf, ref_space, is_linear);
// Set up the solver, matrix, and rhs according to the solver selection.
SparseMatrix* matrix = create_matrix(matrix_solver);
Vector* rhs = create_vector(matrix_solver);
Solver* solver = create_linear_solver(matrix_solver, matrix, rhs);
// Initialize the preconditioner in the case of SOLVER_AZTECOO.
if (matrix_solver == SOLVER_AZTECOO)
{
((AztecOOSolver*) solver)->set_solver(iterative_method);
((AztecOOSolver*) solver)->set_precond(preconditioner);
// Using default iteration parameters (see solver/aztecoo.h).
}
// Assemble the reference problem.
info("Assembling on reference mesh (ndof: %d).", Space::get_num_dofs(ref_space));
dp.assemble(matrix, rhs);
// Time measurement.
cpu_time.tick();
// Solve the linear system on reference mesh. If successful, obtain the solution.
info("Solving on reference mesh.");
Solution ref_sln(ref_space->get_mesh());
if(solver->solve()) Solution::vector_to_solution(solver->get_solution(), ref_space, &ref_sln);
else {
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();
// Output solution and mesh with polynomial orders.
if (solution_output)
{
out_fn_vtk(&sln, "sln", as);
out_orders_vtk(&space, "order", as);
//.........这里部分代码省略.........
开发者ID:Veix123,项目名称:hermes,代码行数:101,代码来源:main.cpp
示例9: main
int main(int argc, char **args)
{
// Test variable.
int success_test = 1;
if (argc < 3) error("Not enough parameters.");
// Load the mesh.
Mesh mesh;
H3DReader mloader;
if (!mloader.load(args[1], &mesh)) error("Loading mesh file '%s'.", args[1]);
// Initialize the space according to the
// command-line parameters passed.
int o;
sscanf(args[2], "%d", &o);
Ord3 order(o, o, o);
HcurlSpace space(&mesh, bc_types, NULL, order);
// Initialize the weak formulation.
WeakForm wf;
wf.add_matrix_form(callback(bilinear_form), HERMES_SYM);
wf.add_vector_form(callback(linear_form));
// Initialize the FE problem.
bool is_linear = true;
DiscreteProblem dp(&wf, &space, is_linear);
// Set up the solver, matrix, and rhs according to the solver selection.
SparseMatrix* matrix = create_matrix(matrix_solver);
Vector* rhs = create_vector(matrix_solver);
Solver* solver = create_linear_solver(matrix_solver, matrix, rhs);
// Initialize the preconditioner in the case of SOLVER_AZTECOO.
if (matrix_solver == SOLVER_AZTECOO)
{
((AztecOOSolver*) solver)->set_solver(iterative_method);
((AztecOOSolver*) solver)->set_precond(preconditioner);
// Using default iteration parameters (see solver/aztecoo.h).
}
// Assemble the linear problem.
info("Assembling (ndof: %d).", Space::get_num_dofs(&space));
dp.assemble(matrix, rhs);
// Solve the linear system. If successful, obtain the solution.
info("Solving.");
Solution sln(&mesh);
if(solver->solve()) Solution::vector_to_solution(solver->get_solution(), &space, &sln);
else error ("Matrix solver failed.\n");
ExactSolution ex_sln(&mesh, exact_solution);
// Calculate exact error.
info("Calculating exact error.");
Adapt *adaptivity = new Adapt(&space, HERMES_HCURL_NORM);
bool solutions_for_adapt = false;
double err_exact = adaptivity->calc_err_exact(&sln, &ex_sln, solutions_for_adapt, HERMES_TOTAL_ERROR_ABS);
if (err_exact > EPS)
// Calculated solution is not precise enough.
success_test = 0;
// Clean up.
delete matrix;
delete rhs;
delete solver;
delete adaptivity;
if (success_test) {
info("Success!");
return ERR_SUCCESS;
}
else {
info("Failure!");
return ERR_FAILURE;
}
}
开发者ID:Zhonghua,项目名称:hermes-dev,代码行数:78,代码来源:main.cpp
示例10: main
int main(int argc, char* argv[])
{
// Choose a Butcher's table or define your own.
ButcherTable* bt = new ButcherTable(butcher_table_type);
if (bt->is_explicit()) info("Using a %d-stage explicit R-K method.", bt->get_size());
if (bt->is_diagonally_implicit()) info("Using a %d-stage diagonally implicit R-K method.", bt->get_size());
if (bt->is_fully_implicit()) info("Using a %d-stage fully implicit R-K method.", bt->get_size());
// Turn off adaptive time stepping if R-K method is not embedded.
if (bt->is_embedded() == false && ADAPTIVE_TIME_STEP_ON == true) {
warn("R-K method not embedded, turning off adaptive time stepping.");
ADAPTIVE_TIME_STEP_ON = false;
}
// Load the mesh.
Mesh mesh, basemesh;
H2DReader mloader;
mloader.load("wall.mesh", &basemesh);
// Perform initial mesh refinements.
for(int i = 0; i < INIT_REF_NUM; i++) basemesh.refine_all_elements();
basemesh.refine_towards_boundary(BDY_BOTTOM, INIT_REF_NUM_BDY);
mesh.copy(&basemesh);
// Enter boundary markers.
BCTypes bc_types;
bc_types.add_bc_neumann(Hermes::vector<int>(BDY_RIGHT, BDY_LEFT));
bc_types.add_bc_newton(Hermes::vector<int>(BDY_BOTTOM, BDY_TOP));
// Initialize an H1 space with default shapeset.
H1Space space(&mesh, &bc_types, NULL, P_INIT);
int ndof = Space::get_num_dofs(&space);
info("ndof = %d.", ndof);
// Convert initial condition into a Solution.
Solution* sln_prev_time = new Solution(&mesh, TEMP_INIT);
// Initialize weak formulation.
WeakForm wf;
wf.add_matrix_form(stac_jacobian_vol, stac_jacobian_vol_ord, HERMES_NONSYM, HERMES_ANY, sln_prev_time);
wf.add_vector_form(stac_residual_vol, stac_residual_vol_ord, HERMES_ANY, sln_prev_time);
wf.add_matrix_form_surf(stac_jacobian_bottom, stac_jacobian_bottom_ord, BDY_BOTTOM, sln_prev_time);
wf.add_vector_form_surf(stac_residual_bottom, stac_residual_bottom_ord, BDY_BOTTOM, sln_prev_time);
wf.add_matrix_form_surf(stac_jacobian_top, stac_jacobian_top_ord, BDY_TOP, sln_prev_time);
wf.add_vector_form_surf(stac_residual_top, stac_residual_top_ord, BDY_TOP, sln_prev_time);
// Initialize the FE problem.
bool is_linear = true;
DiscreteProblem dp(&wf, &space, is_linear);
// Create a refinement selector.
H1ProjBasedSelector selector(CAND_LIST, CONV_EXP, H2DRS_DEFAULT_ORDER);
// Visualize initial condition.
char title[100];
ScalarView sln_view("Initial condition", new WinGeom(0, 0, 1500, 360));
OrderView ordview("Initial mesh", new WinGeom(0, 410, 1500, 360));
ScalarView time_error_view("Temporal error", new WinGeom(0, 800, 1500, 360));
time_error_view.fix_scale_width(40);
ScalarView space_error_view("Spatial error", new WinGeom(0, 1220, 1500, 360));
space_error_view.fix_scale_width(40);
sln_view.show(sln_prev_time, HERMES_EPS_VERYHIGH);
ordview.show(&space);
// Graph for time step history.
SimpleGraph time_step_graph;
if (ADAPTIVE_TIME_STEP_ON) info("Time step history will be saved to file time_step_history.dat.");
// Time stepping loop:
double current_time = 0; int ts = 1;
do
{
info("Begin time step %d.", ts);
// Periodic global derefinement.
if (ts > 1 && ts % UNREF_FREQ == 0)
{
info("Global mesh derefinement.");
if (UNREF_LEVEL == 1) mesh.unrefine_all_elements();
else mesh.copy(&basemesh);
space.set_uniform_order(P_INIT);
ndof = Space::get_num_dofs(&space);
}
// Spatial adaptivity loop. Note: sln_prev_time must not be
// changed during spatial adaptivity.
Solution ref_sln;
Solution* time_error_fn;
if (bt->is_embedded() == true) time_error_fn = new Solution(&mesh);
else time_error_fn = NULL;
bool done = false; int as = 1;
double err_est;
do {
// Construct globally refined reference mesh and setup reference space.
Space* ref_space = construct_refined_space(&space);
// Initialize discrete problem on reference mesh.
DiscreteProblem* ref_dp = new DiscreteProblem(&wf, ref_space);
// Runge-Kutta step on the fine mesh.
info("Runge-Kutta time step on fine mesh (t = %g s, tau = %g s, stages: %d).",
//.........这里部分代码省略.........
开发者ID:sajedi,项目名称:hermes,代码行数:101,代码来源:main.cpp
示例11: main
int main(int argc, char* argv[])
{
// Choose a Butcher's table or define your own.
ButcherTable bt(butcher_table_type);
if (bt.is_explicit()) info("Using a %d-stage explicit R-K method.", bt.get_size());
if (bt.is_diagonally_implicit()) info("Using a %d-stage diagonally implicit R-K method.", bt.get_size());
if (bt.is_fully_implicit()) info("Using a %d-stage fully implicit R-K method.", bt.get_size());
// Load the mesh.
Mesh mesh;
H2DReader mloader;
mloader.load("cathedral.mesh", &mesh);
// Perform initial mesh refinements.
for(int i = 0; i < INIT_REF_NUM; i++) mesh.refine_all_elements();
mesh.refine_towards_boundary(BDY_AIR, INIT_REF_NUM_BDY);
mesh.refine_towards_boundary(BDY_GROUND, INIT_REF_NUM_BDY);
// Enter boundary markers.
BCTypes bc_types;
bc_types.add_bc_dirichlet(Hermes::vector<std::string>(BDY_GROUND));
bc_types.add_bc_newton(BDY_AIR);
// Enter Dirichlet boundary values.
BCValues bc_values;
bc_values.add_const(BDY_GROUND, TEMP_INIT);
// Initialize an H1 space with default shapeset.
H1Space space(&mesh, &bc_types, &bc_values, P_INIT);
int ndof = Space::get_num_dofs(&space);
info("ndof = %d.", ndof);
// Previous time level solution (initialized by the external temperature).
Solution u_prev_time(&mesh, TEMP_INIT);
// Initialize weak formulation.
WeakForm wf;
wf.add_matrix_form(callback(stac_jacobian_vol));
wf.add_vector_form(callback(stac_residual_vol));
wf.add_matrix_form_surf(callback(stac_jacobian_surf), BDY_AIR);
wf.add_vector_form_surf(callback(stac_residual_surf), BDY_AIR);
// 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);
// Initialize views.
ScalarView Tview("Temperature", new WinGeom(0, 0, 450, 600));
Tview.set_min_max_range(0,20);
Tview.fix_scale_width(30);
// 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;
bool is_linear = true;
if (!rk_time_step(current_time, time_step, &bt, coeff_vec, &dp, matrix_solver,
verbose, is_linear)) {
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;
// Show the new time level solution.
char title[100];
sprintf(title, "Time %3.2f, exterior temperature %3.5f", current_time, temp_ext(current_time));
Tview.set_title(title);
Tview.show(&u_prev_time);
// Increase counter of time steps.
ts++;
}
while (current_time < T_FINAL);
// Cleanup.
delete [] coeff_vec;
// Wait for the view to be closed.
View::wait();
return 0;
}
开发者ID:schnepp,项目名称:hermes,代码行数:94,代码来源:main.cpp
示例12: main
// Main function.
int main(int argc, char* argv[])
{
// Either use exact constitutive relations (slow) (method 0) or precalculate
// their linear approximations (faster) (method 1) or
// precalculate their quintic polynomial approximations (method 2) -- managed by following loop "Initializing polynomial approximation".
if (CONSTITUTIVE_TABLE_METHOD == 1)
CONSTITUTIVE_TABLES_READY = get_constitutive_tables(ITERATIVE_METHOD);
// Points to be used for polynomial approximation of K(h).
double* points = new double[NUM_OF_INSIDE_PTS];
// The van Genuchten + Mualem K(h) function is approximated by polynomials close to zero in case of CONSTITUTIVE_TABLE_METHOD=1.
// In case of CONSTITUTIVE_TABLE_METHOD=2, all constitutive functions are approximated by polynomials.
info("Initializing polynomial approximations.");
for (int i=0; i < MATERIAL_COUNT; i++) {
info("Processing layer %d", i);
init_polynomials(6 + NUM_OF_INSIDE_PTS, LOW_LIMIT, points, NUM_OF_INSIDE_PTS, i);
}
POLYNOMIALS_READY = true;
if (CONSTITUTIVE_TABLE_METHOD == 2) {
CONSTITUTIVE_TABLES_READY = true ;
//Assign table limit to global definition.
TABLE_LIMIT = INTERVALS_4_APPROX[NUM_OF_INTERVALS-1];
}
// Time measurement.
TimePeriod cpu_time;
cpu_time.tick();
// Load the mesh.
Mesh mesh, basemesh;
H2DReader mloader;
mloader.load(mesh_file, &basemesh);
// Initial refinements.
//basemesh.refine_towards_boundary(1, 1);
//basemesh.refine_towards_boundary(3, 1);
//basemesh.refine_towards_boundary(4, 1);
// Perform initial mesh refinements.
mesh.copy(&basemesh);
for(int i = 0; i < INIT_REF_NUM; i++) mesh.refine_all_elements();
// Enter boundary markers.
BCTypes bc_types;
bc_types.add_bc_dirichlet(BDY_1);
bc_types.add_bc_neumann(Hermes::Tuple<int>(BDY_2, BDY_3, BDY_4));
// Enter Dirichlet boundary values.
BCValues bc_values(&TIME);
bc_values.add_timedep_function(BDY_1, 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);
// Create a selector which will select optimal candidate.
H1ProjBasedSelector selector(CAND_LIST, CONV_EXP, H2DRS_DEFAULT_ORDER);
// Solutions for the time stepping and the Newton's method.
Solution sln, ref_sln, sln_prev_time, sln_prev_iter;
// Assign the function f() to the fine mesh.
sln_prev_time.set_exact(&mesh, init_cond);
sln_prev_iter.set_exact(&mesh, init_cond);
// Initialize the weak formulation.
WeakForm wf;
if (ITERATIVE_METHOD == 1) {
if (TIME_INTEGRATION == 1) {
info("Registering forms for the Newton's method (implicit Euler in time).");
wf.add_matrix_form(jac_form_vol_euler, jac_form_vol_ord, HERMES_NONSYM, HERMES_ANY,
&sln_prev_time);
wf.add_vector_form(res_form_vol_euler, res_form_vol_ord, HERMES_ANY,
&sln_prev_time);
}
else {
info("Registering forms for the Newton's method (Crank-Nicolson in time).");
wf.add_matrix_form(jac_form_vol_cranic, jac_form_vol_ord, HERMES_NONSYM, HERMES_ANY,
&sln_prev_time);
wf.add_vector_form(res_form_vol_cranic, res_form_vol_ord, HERMES_ANY,
&sln_prev_time);
}
}
else {
if (TIME_INTEGRATION == 1) {
info("Registering forms for the Picard's method (implicit Euler in time).");
wf.add_matrix_form(bilinear_form_picard_euler, bilinear_form_picard_euler_ord, HERMES_NONSYM, HERMES_ANY,
&sln_prev_iter);
wf.add_vector_form(linear_form_picard_euler, linear_form_picard_euler_ord, HERMES_ANY,
Hermes::Tuple<MeshFunction*>(&sln_prev_iter, &sln_prev_time));
}
else {
info("Registering forms for the Picard's method (Crank-Nicolson in time).");
error("Not implemented yet.");
wf.add_matrix_form(bilinear_form_picard_euler, bilinear_form_picard_euler_ord, HERMES_NONSYM, HERMES_ANY,
&sln_prev_iter);
//.........这里部分代码省略.........
开发者ID:sriharifez,项目名称:hermes,代码行数:101,代码来源:main.cpp
示例13: main
int main()
{
// Create space, set Dirichlet BC, enumerate basis functions.
Space* space = new Space(A, B, NELEM, DIR_BC_LEFT, DIR_BC_RIGHT, P_INIT, NEQ);
int ndof = Space::get_num_dofs(space);
info("ndof: %d", ndof);
// Initialize the weak formulation.
WeakForm wf;
wf.add_matrix_form(jacobian_vol);
wf.add_vector_form(residual_vol);
wf.add_vector_form_surf(0, residual_surf_right, BOUNDARY_RIGHT);
// 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");
// Cleaning
delete dp;
delete rhs;
delete solver;
delete[] coeff_vec;
delete space;
delete bc;
delete matrix;
info("Done.");
return 0;
}
开发者ID:ppkk,项目名称:hermes-legacy,
|
请发表评论