From: bangerth Date: Mon, 19 May 2008 20:21:24 +0000 (+0000) Subject: Make some things a bit more generic. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=583689a12eba5ebecda165fc7ede30ff4fe06469;p=dealii-svn.git Make some things a bit more generic. git-svn-id: https://svn.dealii.org/trunk@16125 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/examples/step-33/step-33.cc b/deal.II/examples/step-33/step-33.cc index 1ed175886b..66db71abf1 100644 --- a/deal.II/examples/step-33/step-33.cc +++ b/deal.II/examples/step-33/step-33.cc @@ -84,16 +84,11 @@ #include - // And this again is C++ as well as two - // include files from the BOOST library that - // provide us with counted pointers and - // arrays of fixed size: + // And this again is C++: #include #include #include - -#include -#include +#include // To end this section, introduce everythin // in the dealii library into the current @@ -122,7 +117,7 @@ using namespace dealii; template struct EulerEquations { - // First a few variables that + // First a few variables that // describe the various components of our // solution vector in a generic way. This // includes the number of components in the @@ -146,6 +141,50 @@ struct EulerEquations static const unsigned int density_component = dim; static const unsigned int energy_component = dim+1; + // When generating graphical + // output way down in this + // program, we need to specify + // the names of the solution + // variables as well as how the + // various components group into + // vector and scalar fields. We + // could describe this there, but + // in order to keep things that + // have to do with the Euler + // equation localized here and + // the rest of the program as + // generic as possible, we + // provide this sort of + // information in the following + // two functions: + static + std::vector + component_names () + { + std::vector names (dim, "momentum"); + names.push_back ("density"); + names.push_back ("energy_density"); + + return names; + } + + + static + std::vector + component_interpretation () + { + std::vector + data_component_interpretation + (dim, DataComponentInterpretation::component_is_part_of_vector); + data_component_interpretation + .push_back (DataComponentInterpretation::component_is_scalar); + data_component_interpretation + .push_back (DataComponentInterpretation::component_is_scalar); + + return data_component_interpretation; + } + + // Next, we define the gas // constant. We will set it to 1.4 // in its definition immediately @@ -222,9 +261,8 @@ struct EulerEquations (*(W.begin() + energy_component) - compute_kinetic_energy(W))); } - - + // We define the flux function // $F(W)$ as one large matrix. // Each row of this matrix @@ -554,11 +592,8 @@ EulerEquations::Postprocessor:: get_data_component_interpretation () const { std::vector - interpretation; - - for (unsigned int d=0; d solve (Vector &solution); - void refine_grid (); + void compute_refinement_indicators (Vector &indicator) const; + void refine_grid (const Vector &indicator); void output_results () const; - void compute_refinement_indicator (); // The first few member variables @@ -1406,11 +1441,59 @@ class ConservationLaw Vector right_hand_side; - Epetra_SerialComm communicator; - - Epetra_Map *Map; - Epetra_CrsMatrix *Matrix; - Vector indicator; + // This final set of member + // variables (except for the + // object holding all run-time + // parameters at the very bottom) + // deals with the interface we + // have in this program to the + // Trilinos library that provides + // us with linear solvers. + // + // Trilinos is designed to be a + // library that also runs in + // parallel on distributed memory + // systems, so matrices and + // vectors need two things: (i) a + // communicator object that + // facilitates sending messages + // to remote machines, and (ii) a + // description which elements of + // a vector or matrix reside + // locally on a machine and which + // are stored remotely. + // + // We do not actually run the + // current program in parallel, + // and so the objects we use here + // are pretty much dummy objects + // for this purpose: the + // communicator below represents + // a system that includes only a + // single machine, and the index + // map encodes that all elements + // are stored + // locally. Nevertheless, we need + // them. + // + // Furthermore, we need a matrix + // object for the system matrix + // to be used in each Newton + // step. Note that map and matrix + // need to be updated for their + // sizes whenever we refine the + // mesh. In Trilinos, this is + // easiest done by simply + // deleting the previous object + // and creating a new one. To + // minimize hassle and avoid + // memory leaks, we use a + // std::auto_ptr + // instead of a plain pointer for + // this. + Epetra_SerialComm communicator; + std::auto_ptr Map; + std::auto_ptr Matrix; Parameters::AllParameters parameters; }; @@ -1424,9 +1507,7 @@ ConservationLaw::ConservationLaw (const char *input_filename) fe (FE_Q(1), EulerEquations::n_components), dof_handler (triangulation), quadrature (2), - face_quadrature (2), - Map(NULL), - Matrix(NULL) + face_quadrature (2) { ParameterHandler prm; Parameters::AllParameters::declare_parameters (prm); @@ -2119,8 +2200,7 @@ void ConservationLaw::setup_system () // Rebuild the map. In serial this doesn't do much, // but is needed. In parallel, this would desribe // the parallel dof layout. - if (Map) delete Map; - Map = new Epetra_Map(dof_handler.n_dofs(), 0, communicator); + Map.reset (new Epetra_Map(dof_handler.n_dofs(), 0, communicator)); // Epetra can build a more efficient matrix if // one knows ahead of time the maximum number of @@ -2133,9 +2213,7 @@ void ConservationLaw::setup_system () // the constructor that optimizes // with the existing lengths per row // variable. - if (Matrix != 0) - delete Matrix; - Matrix = new Epetra_CrsMatrix(Copy, *Map, &row_lengths[0], true); + Matrix.reset (new Epetra_CrsMatrix(Copy, *Map, &row_lengths[0], true)); // We add the sparsity pattern to the matrix by // inserting zeros. @@ -2188,7 +2266,7 @@ ConservationLaw::solve (Vector &newton_update) // changing th string given to the // Create function. Epetra_LinearProblem prob; - prob.SetOperator(Matrix); + prob.SetOperator(Matrix.get()); Amesos_BaseSolver *solver = Amesos().Create ("Amesos_Klu", prob); Assert (solver != NULL, ExcInternalError()); @@ -2255,7 +2333,7 @@ ConservationLaw::solve (Vector &newton_update) Solver.SetAztecParam(AZ_ilut_fill, parameters.ilut_fill); Solver.SetAztecParam(AZ_athresh, parameters.ilut_atol); Solver.SetAztecParam(AZ_rthresh, parameters.ilut_rtol); - Solver.SetUserMatrix(Matrix); + Solver.SetUserMatrix(Matrix.get()); // Run the solver iteration. Collect the number // of iterations and the residual. @@ -2273,7 +2351,9 @@ ConservationLaw::solve (Vector &newton_update) // simply use the density squared, which selects // shocks with some success. template -void ConservationLaw::compute_refinement_indicator () +void +ConservationLaw:: +compute_refinement_indicators (Vector &refinement_indicators) const { const unsigned int dofs_per_cell = dof_handler.get_fe().dofs_per_cell; std::vector dofs (dofs_per_cell); @@ -2303,21 +2383,21 @@ void ConservationLaw::compute_refinement_indicator () fe_v.get_function_values(predictor, U); fe_v.get_function_grads(predictor, dU); - indicator(cell_no) = 0; + refinement_indicators(cell_no) = 0; for (unsigned int q = 0; q < n_q_points; q++) { double ng = 0; for (unsigned int d = 0; d < dim; d++) ng += dU[q][EulerEquations::density_component][d]*dU[q][EulerEquations::density_component][d]; - indicator(cell_no) += std::log(1+std::sqrt(ng)); + refinement_indicators(cell_no) += std::log(1+std::sqrt(ng)); } - indicator(cell_no) /= n_q_points; + refinement_indicators(cell_no) /= n_q_points; } } template -void ConservationLaw::refine_grid () +void ConservationLaw::refine_grid (const Vector &refinement_indicators) { SolutionTransfer soltrans(dof_handler); @@ -2335,11 +2415,11 @@ void ConservationLaw::refine_grid () cell->clear_coarsen_flag(); cell->clear_refine_flag(); if (cell->level() < parameters.shock_levels && - std::fabs(indicator(cell_no)) > parameters.shock_val ) { + std::fabs(refinement_indicators(cell_no)) > parameters.shock_val ) { cell->set_refine_flag(); } else { if (cell->level() > 0 && - std::fabs(indicator(cell_no)) < 0.75*parameters.shock_val) + std::fabs(refinement_indicators(cell_no)) < 0.75*parameters.shock_val) cell->set_coarsen_flag(); } } @@ -2382,9 +2462,6 @@ void ConservationLaw::refine_grid () current_solution.reinit(dof_handler.n_dofs()); current_solution = old_solution; right_hand_side.reinit (dof_handler.n_dofs()); - - indicator.reinit(triangulation.n_active_cells()); - } @@ -2397,26 +2474,14 @@ void ConservationLaw::output_results () const DataOut data_out; data_out.attach_dof_handler (dof_handler); - std::vector old_solution_names (dim, "momentum"); - old_solution_names.push_back ("density"); - old_solution_names.push_back ("energy_density"); - - std::vector - data_component_interpretation - (dim, DataComponentInterpretation::component_is_part_of_vector); - data_component_interpretation - .push_back (DataComponentInterpretation::component_is_scalar); - data_component_interpretation - .push_back (DataComponentInterpretation::component_is_scalar); - data_out.add_data_vector (old_solution, old_solution_names, + data_out.add_data_vector (old_solution, + EulerEquations::component_names (), DataOut::type_dof_data, - data_component_interpretation); + EulerEquations::component_interpretation ()); data_out.add_data_vector (old_solution, postprocessor); - data_out.add_data_vector (indicator, "error"); - data_out.build_patches (); static unsigned int output_file_number = 0; @@ -2461,7 +2526,6 @@ void ConservationLaw::run () current_solution.reinit (dof_handler.n_dofs()); predictor.reinit (dof_handler.n_dofs()); right_hand_side.reinit (dof_handler.n_dofs()); - indicator.reinit(triangulation.n_active_cells()); setup_system(); @@ -2476,8 +2540,9 @@ void ConservationLaw::run () if (parameters.do_refine == true) for (unsigned int i = 0; i < parameters.shock_levels; i++) { - compute_refinement_indicator(); - refine_grid(); + Vector refinement_indicators (triangulation.n_active_cells()); + compute_refinement_indicators(refinement_indicators); + refine_grid(refinement_indicators); setup_system(); VectorTools::interpolate(dof_handler, @@ -2577,7 +2642,8 @@ void ConservationLaw::run () old_solution = current_solution; - compute_refinement_indicator(); + Vector refinement_indicators (triangulation.n_active_cells()); + compute_refinement_indicators(refinement_indicators); time += parameters.time_step; @@ -2593,7 +2659,7 @@ void ConservationLaw::run () // Refine, if refinement is selected. if (parameters.do_refine == true) { - refine_grid(); + refine_grid(refinement_indicators); setup_system(); newton_update.reinit (dof_handler.n_dofs());