From: allmaras Date: Wed, 12 Sep 2007 23:39:11 +0000 (+0000) Subject: Updated step-29, now using new DataPostprocessor class to output data X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=5fb4cb902ea4a187f9709fa73b585de6cf7fba2c;p=dealii-svn.git Updated step-29, now using new DataPostprocessor class to output data git-svn-id: https://svn.dealii.org/trunk@15201 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/examples/step-29/step-29.cc b/deal.II/examples/step-29/step-29.cc index cafd632328..2e6542348b 100644 --- a/deal.II/examples/step-29/step-29.cc +++ b/deal.II/examples/step-29/step-29.cc @@ -1,75 +1,81 @@ -/* $Id: step-28.cc 14713 2007-05-27 04:05:41Z bangerth $ */ -/* Version: $Name: $ */ -/* */ -/* Copyright (C) 2007 by the deal.II authors and Moritz Allmaras */ -/* */ -/* This file is subject to QPL and may not be distributed */ -/* without copyright and license information. Please refer */ -/* to the file deal.II/doc/license.html for the text and */ -/* further information on this license. */ + // @sect3{Include files} + + // The following header files are unchanged + // from step-7 and have been discussed before: #include #include #include -#include -#include + #include #include #include -#include + #include #include #include #include #include + #include #include #include -#include + #include #include + #include #include #include + #include -#define DIM 2 + // This header file is needed for + // the ParameterHandler class that we will + // use to read parameters from a + // configuration file during runtime: +#include +#include +#include + +#define DIM 2 using namespace dealii; + template class DirichletBoundaryValues : public Function { public: DirichletBoundaryValues() : Function (2) {}; - - virtual void vector_value ( const Point &p, - Vector &values) const; - - virtual void vector_value_list (const std::vector > & points, - std::vector > & value_list) const; + + virtual void vector_value (const Point &p, + Vector &values) const; + + virtual void vector_value_list (const std::vector > &points, + std::vector > &value_list) const; }; template inline -void DirichletBoundaryValues::vector_value ( const Point & /*p*/, - Vector & values) const +void DirichletBoundaryValues::vector_value (const Point &/*p*/, + Vector &values) const { Assert (values.size() == 2, ExcDimensionMismatch (values.size(), 2)); - + values(0) = 1; values(1) = 0; } template -void DirichletBoundaryValues::vector_value_list (const std::vector > & points, - std::vector > & value_list) const +void DirichletBoundaryValues::vector_value_list (const std::vector > &points, + std::vector > &value_list) const { Assert (value_list.size() == points.size(), - ExcDimensionMismatch (value_list.size(), points.size())); + ExcDimensionMismatch (value_list.size(), points.size())); for (unsigned int p=0; p::vector_value (points[p], value_list[p]); @@ -81,53 +87,142 @@ class ParameterReader : public Subscriptor public: ParameterReader(ParameterHandler &); void read_parameters(); - + private: void declare_parameters(); - ParameterHandler *prm; + ParameterHandler &prm; }; ParameterReader::ParameterReader(ParameterHandler ¶mhandler) - : - prm(¶mhandler) + : + prm(paramhandler) {} +void ParameterReader::read_parameters() +{ + declare_parameters(); + + const std::string parameter_file = "step-29.prm"; + prm.read_input (parameter_file); +} + + void ParameterReader::declare_parameters() { - prm->declare_entry( "Number of refinements", "5", - Patterns::Integer(1,10), - "Number of global mesh refinement steps " - "applied to initial coarse grid"); + prm.enter_subsection ("Mesh & geometry parameters"); + + prm.declare_entry("Number of refinements", "6", + Patterns::Integer(1,10), + "Number of global mesh refinement steps " + "applied to initial coarse grid"); + + prm.declare_entry("Focal distance", "0.3", + Patterns::Double(0), + "Distance of the focal point of the lens " + "to the x-axis (or xy-plane in 3D)"); - prm->declare_entry( "Focal distance", "0.3", - Patterns::Double(0), - "Distance of the focal point of the lens " - "to the x-axis (or xy-plane in 3D)"); + prm.leave_subsection (); - prm->declare_entry( "c", "1.5e5", - Patterns::Double(), - "Wave speed"); + prm.enter_subsection ("Physical constants"); - prm->declare_entry( "omega", "1.5e7", - Patterns::Double(), - "Frequency"); + prm.declare_entry("c", "1.5e5", + Patterns::Double(), + "Wave speed"); - prm->declare_entry( "Output file", "solution", - Patterns::Anything(), - "Name of the output file (without extension)"); + prm.declare_entry("omega", "5.0e7", + Patterns::Double(), + "Frequency"); - DataOutInterface<1>::declare_parameters (*prm); + prm.leave_subsection (); + + prm.enter_subsection ("Output parameters"); + + prm.declare_entry("Output file", "solution", + Patterns::Anything(), + "Name of the output file (without extension)"); + + DataOutInterface<1>::declare_parameters (prm); + + prm.leave_subsection (); } -void ParameterReader::read_parameters() +template +class Postprocessor : public DataPostprocessor { - declare_parameters(); + public: - const std::string parameter_file = "step-29.prm"; - prm->read_input (parameter_file); + using DataPostprocessor::compute_derived_quantities; + + void compute_derived_quantities ( + std::vector< Vector< double > > &, + const std::vector< Vector< double > > &, + const std::vector< std::vector< Tensor< 1, dim > > > &, + const std::vector< std::vector< Tensor< 2, dim > > > &, + const std::vector< Point< dim > > & + ) const; + + std::vector get_names () const; + UpdateFlags get_needed_update_flags () const; + unsigned int n_output_variables () const; +}; + + +template +std::vector +Postprocessor::get_names() const +{ + std::vector field_names; + + field_names.push_back("Re_u"); + field_names.push_back("Im_u"); + field_names.push_back("Intensity"); + + return field_names; +} + + +template +UpdateFlags +Postprocessor::get_needed_update_flags () const +{ + return update_values; +} + + +template +unsigned int +Postprocessor::n_output_variables () const +{ + return 3; +} + + +template +void +Postprocessor::compute_derived_quantities ( + std::vector< Vector< double > > &computed_quantities, + const std::vector< Vector< double > > &uh, + const std::vector< std::vector< Tensor< 1, dim > > > &/*duh*/, + const std::vector< std::vector< Tensor< 2, dim > > > &/*dduh*/, + const std::vector< Point< dim > > &/*normals*/ + ) const +{ + Assert(computed_quantities.size() == uh.size(), + ExcDimensionMismatch (computed_quantities.size(), uh.size())); + + for (unsigned int i=0; i triangulation; + DoFHandler dof_handler; + FESystem fe; - Triangulation triangulation; - DoFHandler dof_handler; - FESystem fe; - - SparsityPattern sparsity_pattern; - SparseMatrix system_matrix; - - Vector solution, system_rhs, intensity; + SparsityPattern sparsity_pattern; + SparseMatrix system_matrix; + Vector solution, system_rhs; - const double c, omega; + const double c, omega; }; template -UltrasoundProblem::UltrasoundProblem (ParameterHandler& param) - : - prm(param), - dof_handler(triangulation), - fe(FE_Q(1), 2), - c(prm.get_double("c")), - omega(prm.get_double("omega")) +double +UltrasoundProblem::get_omega (ParameterHandler &prm) +{ + prm.enter_subsection ("Physical constants"); + double omega_tmp = prm.get_double("omega"); + prm.leave_subsection (); + + return omega_tmp; +} + + +template +double +UltrasoundProblem::get_c (ParameterHandler &prm) +{ + prm.enter_subsection ("Physical constants"); + double c_tmp = prm.get_double("c"); + prm.leave_subsection (); + + return c_tmp; +} + + +template +UltrasoundProblem::UltrasoundProblem (ParameterHandler& param) + : + prm(param), + dof_handler(triangulation), + fe(FE_Q(1), 2), + c(get_c(prm)), + omega(get_omega(prm)) {} @@ -183,33 +303,41 @@ UltrasoundProblem::~UltrasoundProblem () template void UltrasoundProblem::make_grid () { + prm.enter_subsection ("Mesh & geometry parameters"); + + const double focal_distance = prm.get_double("Focal distance"); + const unsigned int N_ref = prm.get_integer("Number of refinements"); + + prm.leave_subsection (); + GridGenerator::subdivided_hyper_cube (triangulation, 5, 0, 1); - + const Point transducer = (dim == 2) ? Point (0.5, 0.0) : Point (0.5, 0.5, 0.0), - focal_point = (dim == 2) ? - Point (0.5, prm.get_double("Focal distance")) : - Point (0.5, 0.5, prm.get_double("Focal distance")); + focal_point = (dim == 2) ? + Point (0.5, focal_distance) : + Point (0.5, 0.5, focal_distance); - double radius = sqrt( (focal_point.distance(transducer) * + double radius = sqrt( (focal_point.distance(transducer) * focal_point.distance(transducer)) + ((dim==2) ? 0.01 : 0.02)); typename Triangulation::cell_iterator cell = triangulation.begin (), endc = triangulation.end(); - + for (; cell!=endc; ++cell) for (unsigned int face=0; face::faces_per_cell; ++face) if ( cell->face(face)->at_boundary() && ((cell->face(face)->center() - transducer).square() < 0.01) ) - cell->face(face)->set_boundary_indicator (1); + + cell->face(face)->set_boundary_indicator (1); const HyperBallBoundary boundary(focal_point, radius); triangulation.set_boundary(1, boundary); - triangulation.refine_global (prm.get_integer("Number of refinements")); + triangulation.refine_global (N_ref); deallog << " Number of active cells: " << triangulation.n_active_cells() @@ -225,40 +353,39 @@ void UltrasoundProblem::setup_system () dof_handler.distribute_dofs (fe); deallog << " Number of degrees of freedom: " - << dof_handler.n_dofs() - << std::endl; + << dof_handler.n_dofs() + << std::endl; + + sparsity_pattern.reinit (dof_handler.n_dofs(), + dof_handler.n_dofs(), + dof_handler.max_couplings_between_dofs()); - sparsity_pattern.reinit ( dof_handler.n_dofs(), - dof_handler.n_dofs(), - dof_handler.max_couplings_between_dofs()); DoFTools::make_sparsity_pattern (dof_handler, sparsity_pattern); sparsity_pattern.compress(); system_matrix.reinit (sparsity_pattern); system_rhs.reinit (dof_handler.n_dofs()); solution.reinit (dof_handler.n_dofs()); - intensity.reinit(triangulation.n_active_cells()); } template void UltrasoundProblem::assemble_system () -{ +{ const double om2 = omega * omega; const double c2 = c * c; - - QGauss quadrature_formula(3); - QGauss face_quadrature_formula(3); - const unsigned int n_q_points = quadrature_formula.n_quadrature_points; - const unsigned int n_face_q_points = face_quadrature_formula.n_quadrature_points; + QGauss quadrature_formula(2); + QGauss face_quadrature_formula(2); - const unsigned int dofs_per_cell = fe.dofs_per_cell; + const unsigned int n_q_points = quadrature_formula.n_quadrature_points, + n_face_q_points = face_quadrature_formula.n_quadrature_points, + dofs_per_cell = fe.dofs_per_cell; - FullMatrix cell_matrix (dofs_per_cell, dofs_per_cell); + FullMatrix cell_matrix (dofs_per_cell, dofs_per_cell); std::vector local_dof_indices (dofs_per_cell); - + FEValues fe_values (fe, quadrature_formula, update_values | update_gradients | update_JxW_values); @@ -269,67 +396,70 @@ void UltrasoundProblem::assemble_system () typename DoFHandler::active_cell_iterator cell = dof_handler.begin_active(), endc = dof_handler.end(); + for (; cell!=endc; ++cell) + { + cell_matrix = 0; + fe_values.reinit (cell); + + for (unsigned int i=0; iat_boundary()) - for (unsigned int face=0; face::faces_per_cell; ++face) - if ( cell->face(face)->at_boundary() && - (cell->face(face)->boundary_indicator() == 0) ) - { - fe_face_values.reinit (cell, face); - - for (unsigned int i=0; iget_dof_indices (local_dof_indices); - - for (unsigned int i=0; iat_boundary()) + for (unsigned int face=0; face::faces_per_cell; ++face) + if (cell->face(face)->at_boundary() && + (cell->face(face)->boundary_indicator() == 0) ) + { + fe_face_values.reinit (cell, face); + + for (unsigned int i=0; iget_dof_indices (local_dof_indices); + + for (unsigned int i=0; i boundary_values; VectorTools::interpolate_boundary_values (dof_handler, 1, DirichletBoundaryValues(), boundary_values); + MatrixTools::apply_boundary_values (boundary_values, system_matrix, solution, @@ -338,9 +468,9 @@ void UltrasoundProblem::assemble_system () template -void UltrasoundProblem::solve () +void UltrasoundProblem::solve () { - SparseDirectUMFPACK A_direct; + SparseDirectUMFPACK A_direct; A_direct.initialize(system_matrix); A_direct.vmult(solution,system_rhs); @@ -348,56 +478,31 @@ void UltrasoundProblem::solve () template -void UltrasoundProblem::postprocess () +void UltrasoundProblem::output_results () const { - QMidpoint midpoint_rule; - FEValues fe_values(fe, midpoint_rule, update_values); - - std::vector > values; - values.resize(1); - values[0].reinit(2); + Postprocessor pproc; + DataOut data_out; - typename DoFHandler::active_cell_iterator - cell = dof_handler.begin_active(), - endc = dof_handler.end(); - for (unsigned int i=0; cell!=endc; ++cell, ++i) - { - fe_values.reinit(cell); - - fe_values.get_function_values (solution, values); + data_out.attach_dof_handler (dof_handler); - intensity(i) = sqrt( values[0](0)*values[0](0) - + values[0](1)*values[0](1)); - } -} + prm.enter_subsection("Output parameters"); + const std::string output_file = prm.get("Output file"), + output_format = prm.get("Output format"); -template -void UltrasoundProblem::output_results () const -{ - const std::string output_file(prm.get("Output file")), - output_format(prm.get("Output format")); - - DataOutBase::OutputFormat format = DataOutBase::parse_output_format(output_format); + data_out.parse_parameters(prm); - const std::string filename = output_file + - DataOutBase::default_suffix(format); + prm.leave_subsection (); - std::ofstream output (filename.c_str()); + DataOutBase::OutputFormat format = DataOutBase::parse_output_format(output_format); - DataOut data_out; - data_out.parse_parameters(prm); - data_out.attach_dof_handler (dof_handler); + const std::string filename = output_file + + DataOutBase::default_suffix(format); - std::vector solution_names; - solution_names.push_back ("Re_u"); - solution_names.push_back ("Im_u"); - data_out.add_data_vector (solution, solution_names); - - data_out.add_data_vector (intensity, "intensity"); + std::ofstream output (filename.c_str()); + data_out.add_data_vector (solution, pproc); data_out.build_patches (); - data_out.write (output, format); } @@ -409,7 +514,6 @@ void UltrasoundProblem::run () setup_system (); assemble_system (); solve (); - postprocess (); output_results (); } @@ -417,39 +521,40 @@ void UltrasoundProblem::run () int main () { try - { - ParameterHandler prm; - ParameterReader param(prm); - param.read_parameters(); - - Assert (DIM > 1, ExcNotImplemented()); - - UltrasoundProblem ultrasound_problem (prm); - ultrasound_problem.run (); - } + { + ParameterHandler prm; + + ParameterReader param(prm); + param.read_parameters(); + + Assert (DIM > 1, ExcNotImplemented()); + + UltrasoundProblem ultrasound_problem (prm); + ultrasound_problem.run (); + } catch (std::exception &exc) - { - std::cerr << std::endl << std::endl - << "----------------------------------------------------" - << std::endl; - std::cerr << "Exception on processing: " << std::endl - << exc.what() << std::endl - << "Aborting!" << std::endl - << "----------------------------------------------------" - << std::endl; - return 1; - } + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + } catch (...) - { - std::cerr << std::endl << std::endl - << "----------------------------------------------------" - << std::endl; - std::cerr << "Unknown exception!" << std::endl - << "Aborting!" << std::endl - << "----------------------------------------------------" - << std::endl; - return 1; - } + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + } return 0; } diff --git a/deal.II/examples/step-29/step-29.prm b/deal.II/examples/step-29/step-29.prm index 9d50174c58..9f48aef3ca 100644 --- a/deal.II/examples/step-29/step-29.prm +++ b/deal.II/examples/step-29/step-29.prm @@ -1,109 +1,28 @@ # Listing of Parameters # --------------------- -# Distance of the focal point of the lens to the x-axis (or xy-plane in 3D) -set Focal distance = 0.3 -# Number of global mesh refinement steps applied to initial coarse grid -set Number of refinements = 5 +subsection Mesh & geometry parameters + # Distance of the focal point of the lens to the x-axis (or xy-plane in 3D) + set Focal distance = 0.3 -# Name of the output file (without extension) -set Output file = solution - -# A name for the output format to be used -set Output format = vtk - -# Number of subdivisions of each mesh cell -set Subdivisions = 1 - -# Wave speed -set c = 1.5e5 - -# Frequency -set omega = 1.5e7 - - -subsection DX output parameters - # Output format of vertex coordinates, which is either a text representation - # (ascii) or binary floating point values of 32 or 64 bits length - set Coordinates format = ascii - - # Output format of data values, which is either a text representation - # (ascii) or binary floating point values of 32 or 64 bits length - set Data format = ascii - - # Output format of integer numbers, which is either a text representation - # (ascii) or binary integer values of 32 or 64 bits length - set Integer format = ascii - - # A boolean field indicating whether neighborship information between cells - # is to be written to the OpenDX output file - set Write neighbors = true + # Number of global mesh refinement steps applied to initial coarse grid + set Number of refinements = 6 end -subsection Eps output parameters - # Angle of the viewing position against the vertical axis - set Azimut angle = 60 - - # Name of a color function used to colorize mesh lines and/or cell - # interiors - set Color function = default - - # Whether the interior of cells shall be shaded - set Color shading of interior of cells = true - - # Whether the mesh lines, or only the surface should be drawn - set Draw mesh lines = true - - # Whether only the mesh lines, or also the interior of cells should be - # plotted. If this flag is false, then one can see through the mesh - set Fill interior of cells = true +subsection Physical constants + # Wave speed + set c = 1.5e5 - # Number of the input vector that is to be used to generate color - # information - set Index of vector for color = 0 - - # Number of the input vector that is to be used to generate height - # information - set Index of vector for height = 0 - - # The width in which the postscript renderer is to plot lines - set Line widths in eps units = 0.5 - - # Whether width or height should be scaled to match the given size - set Scale to width or height = width - - # Scaling for the z-direction relative to the scaling used in x- and - # y-directions - set Scaling for z-axis = 1 - - # The size (width or height) to which the eps output file is to be scaled - set Size (width or height) in eps units = 300 - - # Angle of the viewing direction against the y-axis - set Turn angle = 30 + # Frequency + set omega = 5.0e7 end -subsection Povray output parameters - # Whether camera and lightling information should be put into an external - # file "data.inc" or into the POVRAY input file - set Include external file = true +subsection Output parameters + # Name of the output file (without extension) + set Output file = solution - # Whether POVRAY should use bicubic patches - set Use bicubic patches = false - - # A flag indicating whether POVRAY should use smoothed triangles instead of - # the usual ones - set Use smooth triangles = false -end - - -subsection UCD output parameters - # A flag indicating whether a comment should be written to the beginning of - # the output file indicating date and time of creation as well as the - # creating program - set Write preamble = true + # A name for the output format to be used + set Output format = gmv end - -