From: bangerth Date: Sun, 18 May 2008 00:21:14 +0000 (+0000) Subject: Hook the postprocessor into the program. Other minor changes. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=d31fcb24c4eb2a530fe2ee08afb894cbfa870803;p=dealii-svn.git Hook the postprocessor into the program. Other minor changes. git-svn-id: https://svn.dealii.org/trunk@16114 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 99be279ecd..43c5633c22 100644 --- a/deal.II/examples/step-33/step-33.cc +++ b/deal.II/examples/step-33/step-33.cc @@ -95,7 +95,7 @@ using namespace dealii; - // @sect3{Flux function definition} + // @sect3{Euler equation specifics} // Here we define the flux function for this // particular system of conservation laws, @@ -591,10 +591,6 @@ n_output_variables () const } -template class EulerEquations<2>::Postprocessor; - - - // @sect3{Run time parameter handling} // Our next job is to define a few @@ -987,7 +983,6 @@ class ConsLaw void output_results (const unsigned int cycle) const; void initialize(); void estimate(); - void postprocess(); void compute_predictor(); static const unsigned int max_n_boundaries = 10; @@ -1011,8 +1006,6 @@ class ConsLaw // An estimate of the next time value; used for adaptivity and as a // guess for the next Newton iteration. Vector predictor; - // Values after post-processing (used to output the physical variables). - Vector ppsolution; // The solution to the linear problem during the Newton iteration Vector dsolution; Vector right_hand_side; @@ -1762,7 +1755,6 @@ void ConsLaw::initialize_system () solution.reinit (dof_handler.n_dofs()); nlsolution.reinit (dof_handler.n_dofs()); predictor.reinit (dof_handler.n_dofs()); - ppsolution.reinit (dof_handler.n_dofs()); dsolution.reinit (dof_handler.n_dofs()); right_hand_side.reinit (dof_handler.n_dofs()); indicator.reinit(triangulation.n_active_cells()); @@ -1942,98 +1934,6 @@ void ConsLaw::solve (Vector &dsolution, int &niter, double &lin_res niter = Solver.NumIters(); lin_residual = Solver.TrueResidual(); } -} - - // @sect3{Postprocessing and Output} Recover - // the physical variables from the - // conservative variables so that output will - // be (perhaps) more meaningfull. -template -void ConsLaw::postprocess() { - const unsigned int dofs_per_cell = dof_handler.get_fe().dofs_per_cell; - std::vector dofs (dofs_per_cell); - UpdateFlags update_flags = update_values - | update_gradients - | update_q_points - | update_JxW_values; - UpdateFlags update_flags1 = update_values - | update_gradients - | update_q_points - | update_JxW_values; - - QGauss quadrature_formula(4); - - const std::vector > &us = fe.base_element(0).get_unit_support_points(); - - - Quadrature unit_support(us); - - int n_q_points = quadrature_formula.n_quadrature_points; - int n_uq_points = unit_support.n_quadrature_points; - - FEValues fe_v ( - mapping, fe, quadrature_formula, update_flags); - - FEValues fe_v_unit ( - mapping, fe, unit_support, update_flags1); - - std::vector > U(n_uq_points, - Vector(EulerEquations::n_components)); - std::vector > UU(n_q_points, - Vector(EulerEquations::n_components)); - std::vector > > dU(n_uq_points, - std::vector >(EulerEquations::n_components)); - - typename DoFHandler::active_cell_iterator - cell = dof_handler.begin_active(), - endc = dof_handler.end(); - - // Loop the cells - for (unsigned int cell_no=0; cell!=endc; ++cell, ++cell_no) { - cell->get_dof_indices (dofs); - fe_v_unit.reinit(cell); - fe_v.reinit(cell); - - fe_v_unit.get_function_values(solution, U); - fe_v_unit.get_function_grads(solution, dU); - fe_v.get_function_values(solution, UU); - - for (unsigned int q = 0; q < fe_v.get_fe().base_element(0).n_dofs_per_cell(); q++) - { - unsigned int didx - = fe_v.get_fe().component_to_system_index(EulerEquations::density_component, q); - unsigned int eidx - = fe_v.get_fe().component_to_system_index(EulerEquations::energy_component, q); - double rho_normVsqr = 0; - for (unsigned int d = 0; d < dim; d++) - { - unsigned int vidx = fe_v.get_fe().component_to_system_index(d, q); - ppsolution(dofs[vidx]) = solution(dofs[vidx])/solution(dofs[didx]); - rho_normVsqr += solution(dofs[vidx])*solution(dofs[vidx]); - } - rho_normVsqr /= solution(dofs[didx]); - // Pressure - ppsolution(dofs[eidx]) - = (EulerEquations::gas_gamma-1.0)*(solution(dofs[eidx]) - 0.5*rho_normVsqr); - - // Either output density or gradient - // squared of density, depending on - // what the user wants. -//TODO: if schlieren plot then simply use a postprocessor - if (output_params.schlieren_plot == false) - ppsolution(dofs[didx]) = solution(dofs[didx]); - else - { - double ng = 0; - for (unsigned int i = 0; i < dim; i++) - ng += dU[q][EulerEquations::density_component][i]*dU[q][EulerEquations::density_component][i]; - ng = std::sqrt(ng); - ppsolution(dofs[didx]) = ng; - } - } - - } // cell - } // Loop and assign a value for refinement. We @@ -2147,7 +2047,6 @@ void ConsLaw::refine_grid () // resize these vectors for the new grid. nlsolution.reinit(dof_handler.n_dofs()); - ppsolution.reinit(dof_handler.n_dofs()); nlsolution = solution; dsolution.reinit (dof_handler.n_dofs()); right_hand_side.reinit (dof_handler.n_dofs()); @@ -2159,15 +2058,16 @@ void ConsLaw::refine_grid () template void ConsLaw::output_results (const unsigned int cycle) const { - char filename[512]; - std::sprintf(filename, "solution-%03d.vtk", cycle); - std::ofstream output (filename); + std::string filename = "solution-" + + Utilities::int_to_string (cycle, 3) + + ".vtk"; + std::ofstream output (filename.c_str()); DataOut data_out; data_out.attach_dof_handler (dof_handler); - std::vector solution_names (dim, "velocity"); + std::vector solution_names (dim, "momentum"); solution_names.push_back ("density"); - solution_names.push_back ("pressure"); + solution_names.push_back ("energy_density"); std::vector data_component_interpretation @@ -2177,15 +2077,18 @@ void ConsLaw::output_results (const unsigned int cycle) const data_component_interpretation .push_back (DataComponentInterpretation::component_is_scalar); - data_out.add_data_vector (ppsolution, solution_names, + data_out.add_data_vector (solution, solution_names, DataOut::type_dof_data, data_component_interpretation); + typename EulerEquations::Postprocessor + postprocessor (output_params.schlieren_plot); + data_out.add_data_vector (solution, postprocessor); + data_out.add_data_vector (indicator, "error"); + data_out.build_patches (); data_out.write_vtk (output); - - output.close(); } // @sect3{Parsing the Input Deck} @@ -2430,7 +2333,7 @@ void ConsLaw::run () initialize(); predictor = solution; } - postprocess(); + output_results (nstep); // Determine when we will output next. @@ -2509,11 +2412,8 @@ void ConsLaw::run () solution = nlsolution; - estimate(); - postprocess(); - T += dT; // Output if it is time.