From 180bc079a4d8473d6b837677d02914810921d32e Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Mon, 27 Aug 2007 20:16:18 +0000 Subject: [PATCH] add git-svn-id: https://svn.dealii.org/trunk@15076 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-29/step-29.cc | 455 +++++++++++++++++++++++++++ deal.II/examples/step-29/step-29.prm | 109 +++++++ 2 files changed, 564 insertions(+) create mode 100644 deal.II/examples/step-29/step-29.cc create mode 100644 deal.II/examples/step-29/step-29.prm diff --git a/deal.II/examples/step-29/step-29.cc b/deal.II/examples/step-29/step-29.cc new file mode 100644 index 0000000000..cafd632328 --- /dev/null +++ b/deal.II/examples/step-29/step-29.cc @@ -0,0 +1,455 @@ +/* $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. */ + +#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 + + +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; +}; + + +template +inline +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 +{ + Assert (value_list.size() == points.size(), + ExcDimensionMismatch (value_list.size(), points.size())); + + for (unsigned int p=0; p::vector_value (points[p], value_list[p]); +} + + +class ParameterReader : public Subscriptor +{ + public: + ParameterReader(ParameterHandler &); + void read_parameters(); + + private: + void declare_parameters(); + ParameterHandler *prm; +}; + + +ParameterReader::ParameterReader(ParameterHandler ¶mhandler) + : + prm(¶mhandler) +{} + + +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->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( "c", "1.5e5", + Patterns::Double(), + "Wave speed"); + + prm->declare_entry( "omega", "1.5e7", + Patterns::Double(), + "Frequency"); + + prm->declare_entry( "Output file", "solution", + Patterns::Anything(), + "Name of the output file (without extension)"); + + DataOutInterface<1>::declare_parameters (*prm); +} + + +void ParameterReader::read_parameters() +{ + declare_parameters(); + + const std::string parameter_file = "step-29.prm"; + prm->read_input (parameter_file); +} + + +template +class UltrasoundProblem +{ + public: + UltrasoundProblem (ParameterHandler &); + ~UltrasoundProblem (); + void run (); + + private: + void make_grid (); + void setup_system (); + void assemble_system (); + void solve (); + void postprocess (); + void output_results () const; + + ParameterHandler &prm; + + Triangulation triangulation; + DoFHandler dof_handler; + FESystem fe; + + SparsityPattern sparsity_pattern; + SparseMatrix system_matrix; + + Vector solution, system_rhs, intensity; + + 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")) +{} + + +template +UltrasoundProblem::~UltrasoundProblem () +{ + dof_handler.clear(); +} + + +template +void UltrasoundProblem::make_grid () +{ + 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")); + + 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); + + const HyperBallBoundary boundary(focal_point, radius); + triangulation.set_boundary(1, boundary); + + triangulation.refine_global (prm.get_integer("Number of refinements")); + + deallog << " Number of active cells: " + << triangulation.n_active_cells() + << std::endl; + + triangulation.set_boundary(1); +} + + +template +void UltrasoundProblem::setup_system () +{ + dof_handler.distribute_dofs (fe); + + deallog << " Number of degrees of freedom: " + << dof_handler.n_dofs() + << std::endl; + + 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; + + const unsigned int dofs_per_cell = fe.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); + + FEFaceValues fe_face_values (fe, face_quadrature_formula, + update_values | update_JxW_values); + + 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; i boundary_values; + VectorTools::interpolate_boundary_values (dof_handler, + 1, + DirichletBoundaryValues(), + boundary_values); + MatrixTools::apply_boundary_values (boundary_values, + system_matrix, + solution, + system_rhs); +} + + +template +void UltrasoundProblem::solve () +{ + SparseDirectUMFPACK A_direct; + + A_direct.initialize(system_matrix); + A_direct.vmult(solution,system_rhs); +} + + +template +void UltrasoundProblem::postprocess () +{ + QMidpoint midpoint_rule; + FEValues fe_values(fe, midpoint_rule, update_values); + + std::vector > values; + values.resize(1); + values[0].reinit(2); + + 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); + + intensity(i) = sqrt( values[0](0)*values[0](0) + + values[0](1)*values[0](1)); + } +} + + +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); + + const std::string filename = output_file + + DataOutBase::default_suffix(format); + + std::ofstream output (filename.c_str()); + + DataOut data_out; + data_out.parse_parameters(prm); + data_out.attach_dof_handler (dof_handler); + + 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"); + + data_out.build_patches (); + + data_out.write (output, format); +} + + +template +void UltrasoundProblem::run () +{ + make_grid (); + setup_system (); + assemble_system (); + solve (); + postprocess (); + output_results (); +} + + +int main () +{ + try + { + 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; + } + catch (...) + { + 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 new file mode 100644 index 0000000000..9d50174c58 --- /dev/null +++ b/deal.II/examples/step-29/step-29.prm @@ -0,0 +1,109 @@ +# 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 + +# 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 +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 + + # 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 +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 + + # 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 +end + + -- 2.39.5