From 485881faa632c9eef7d4bfea39170abd964e3a2e Mon Sep 17 00:00:00 2001 From: bangerth Date: Thu, 2 Dec 2010 19:22:33 +0000 Subject: [PATCH] Rather than clumsily computing derived quantities within the finite element spaces, use the DataPostprocessor facilities to do so. This also yields a nice implementation of the viscous heating terms. git-svn-id: https://svn.dealii.org/trunk@22907 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-32/step-32.cc | 277 +++++++++++++++++----------- 1 file changed, 165 insertions(+), 112 deletions(-) diff --git a/deal.II/examples/step-32/step-32.cc b/deal.II/examples/step-32/step-32.cc index 776417ae48..662a3879dc 100644 --- a/deal.II/examples/step-32/step-32.cc +++ b/deal.II/examples/step-32/step-32.cc @@ -998,6 +998,8 @@ class BoussinesqFlowProblem void copy_local_to_global_temperature_rhs (const Assembly::CopyData::TemperatureRHS &data); + + class Postprocessor; }; @@ -2881,9 +2883,133 @@ void BoussinesqFlowProblem::solve () } + // @sect4{BoussinesqFlowProblem::output_results} + +template +class BoussinesqFlowProblem::Postprocessor : public DataPostprocessor +{ + public: + Postprocessor (const unsigned int partition); + + virtual + void + compute_derived_quantities_vector (const std::vector > &uh, + const std::vector > > &duh, + const std::vector > > &dduh, + const std::vector > &normals, + const std::vector > &evaluation_points, + std::vector > &computed_quantities) const; + + virtual std::vector get_names () const; + + virtual unsigned int n_output_variables() const; + virtual + std::vector + get_data_component_interpretation () const; + + virtual UpdateFlags get_needed_update_flags () const; + + private: + const unsigned int partition; +}; + + +template +BoussinesqFlowProblem::Postprocessor::Postprocessor (const unsigned int partition) + : + partition (partition) +{} + + +template +std::vector +BoussinesqFlowProblem::Postprocessor::get_names() const +{ + std::vector solution_names (dim, "velocity"); + solution_names.push_back ("p"); + solution_names.push_back ("T"); + solution_names.push_back ("friction_heating"); + solution_names.push_back ("partition"); + return solution_names; +} + + +template +unsigned int +BoussinesqFlowProblem::Postprocessor::n_output_variables() const +{ + return dim + 1 + 1 + 1 + 1; +} + + +template +std::vector +BoussinesqFlowProblem::Postprocessor:: +get_data_component_interpretation () const +{ + std::vector + interpretation (dim, + DataComponentInterpretation::component_is_part_of_vector); + + interpretation.push_back (DataComponentInterpretation::component_is_scalar); + interpretation.push_back (DataComponentInterpretation::component_is_scalar); + interpretation.push_back (DataComponentInterpretation::component_is_scalar); + interpretation.push_back (DataComponentInterpretation::component_is_scalar); + + return interpretation; +} + + +template +UpdateFlags +BoussinesqFlowProblem::Postprocessor::get_needed_update_flags() const +{ + return update_values | update_gradients; +} + + +template +void +BoussinesqFlowProblem::Postprocessor:: +compute_derived_quantities_vector (const std::vector > &uh, + const std::vector > > &duh, + const std::vector > > &/*dduh*/, + const std::vector > &/*normals*/, + const std::vector > &/*evaluation_points*/, + std::vector > &computed_quantities) const +{ + const unsigned int n_quadrature_points = uh.size(); + Assert (duh.size() == n_quadrature_points, ExcInternalError()); + Assert (computed_quantities.size() == n_quadrature_points, ExcInternalError()); + Assert (uh[0].size() == dim+2, ExcInternalError()); + Assert (computed_quantities[0].size()==n_output_variables(),ExcInternalError()); + + for (unsigned int q=0; q grad_u; + for (unsigned int d=0; d strain = symmetrize (grad_u); + computed_quantities[q](dim+2) = EquationData::eta * 2 * strain * strain; + + computed_quantities[q](dim+3) = partition; + } +} - // @sect4{BoussinesqFlowProblem::output_results} // This function does mostly what the // corresponding one did in to @@ -2964,18 +3090,12 @@ void BoussinesqFlowProblem::output_results () estimated_error_per_cell); const FESystem joint_fe (stokes_fe, 1, - temperature_fe, 1, - FE_DGQ(0), 1, - FE_DGQ(0), 1, - temperature_fe, 1); + temperature_fe, 1); DoFHandler joint_dof_handler (triangulation); joint_dof_handler.distribute_dofs (joint_fe); Assert (joint_dof_handler.n_dofs() == - stokes_dof_handler.n_dofs() + - temperature_dof_handler.n_dofs() + - 2*triangulation.n_global_active_cells() + - temperature_dof_handler.n_dofs(), + stokes_dof_handler.n_dofs() + temperature_dof_handler.n_dofs(), ExcInternalError()); TrilinosWrappers::MPI::Vector joint_solution; @@ -2996,7 +3116,8 @@ void BoussinesqFlowProblem::output_results () joint_endc = joint_dof_handler.end(), stokes_cell = stokes_dof_handler.begin_active(), temperature_cell = temperature_dof_handler.begin_active(); - for (unsigned int cell_index=0; joint_cell!=joint_endc; ++joint_cell, ++stokes_cell, ++temperature_cell, ++cell_index) + for (; joint_cell!=joint_endc; + ++joint_cell, ++stokes_cell, ++temperature_cell) if (!joint_cell->is_artificial() && !joint_cell->is_ghost()) { joint_cell->get_dof_indices (local_joint_dof_indices); @@ -3004,114 +3125,46 @@ void BoussinesqFlowProblem::output_results () temperature_cell->get_dof_indices (local_temperature_dof_indices); for (unsigned int i=0; isubdomain_id(); - } - else if (joint_fe.system_to_base_index(i).first.first == 3) - { - Assert (joint_fe.system_to_base_index(i).first.first == 3, - ExcInternalError()); - Assert (joint_fe.system_to_base_index(i).second - == 0, - ExcInternalError()); - joint_solution(local_joint_dof_indices[i]) - = estimated_error_per_cell (cell_index); - } - else - { - Assert (joint_fe.system_to_base_index(i).first.first == 4, - ExcInternalError()); - - joint_solution(local_joint_dof_indices[i]) - = local_temperature_dof_indices - [joint_fe.system_to_base_index(i).second]; - - } - - } + if (joint_fe.system_to_base_index(i).first.first == 0) + { + Assert (joint_fe.system_to_base_index(i).second + < + local_stokes_dof_indices.size(), + ExcInternalError()); + + joint_solution(local_joint_dof_indices[i]) + = stokes_solution(local_stokes_dof_indices + [joint_fe.system_to_base_index(i).second]); + } + else + { + Assert (joint_fe.system_to_base_index(i).first.first == 1, + ExcInternalError()); + Assert (joint_fe.system_to_base_index(i).second + < + local_temperature_dof_indices.size(), + ExcInternalError()); + joint_solution(local_joint_dof_indices[i]) + = temperature_solution(local_temperature_dof_indices + [joint_fe.system_to_base_index(i).second]); + } } } - std::vector joint_solution_names (dim, "velocity"); - joint_solution_names.push_back ("p"); - joint_solution_names.push_back ("T"); - joint_solution_names.push_back ("partition"); - joint_solution_names.push_back ("error"); - joint_solution_names.push_back ("TDofIndex"); - DataOut data_out; - - data_out.attach_dof_handler (joint_dof_handler); - - std::vector - data_component_interpretation - (dim+5, DataComponentInterpretation::component_is_scalar); - for (unsigned int i=0; i::type_dof_data, - data_component_interpretation); - data_out.build_patches (1+std::min(stokes_degree, temperature_degree)); + + Postprocessor postprocessor (Utilities::System:: + get_this_mpi_process(MPI_COMM_WORLD)); + + DataOut data_out; + data_out.attach_dof_handler (joint_dof_handler); + data_out.add_data_vector (locally_relevant_joint_solution, postprocessor); + data_out.build_patches (); static int out_index=0; const std::string filename = ("solution-" + -- 2.39.5