]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
add
authorbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 27 Aug 2007 20:16:18 +0000 (20:16 +0000)
committerbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 27 Aug 2007 20:16:18 +0000 (20:16 +0000)
git-svn-id: https://svn.dealii.org/trunk@15076 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/examples/step-29/step-29.cc [new file with mode: 0644]
deal.II/examples/step-29/step-29.prm [new file with mode: 0644]

diff --git a/deal.II/examples/step-29/step-29.cc b/deal.II/examples/step-29/step-29.cc
new file mode 100644 (file)
index 0000000..cafd632
--- /dev/null
@@ -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 <base/quadrature_lib.h>
+#include <base/function.h>
+#include <base/logstream.h>
+#include <base/parameter_handler.h>
+#include <base/subscriptor.h>
+#include <lac/vector.h>
+#include <lac/full_matrix.h>
+#include <lac/sparse_matrix.h>
+#include <lac/sparse_direct.h>
+#include <grid/tria.h>
+#include <grid/grid_generator.h>
+#include <grid/tria_accessor.h>
+#include <grid/tria_iterator.h>
+#include <grid/tria_boundary_lib.h>
+#include <dofs/dof_handler.h>
+#include <dofs/dof_accessor.h>
+#include <dofs/dof_tools.h>
+#include <fe/fe_system.h>
+#include <fe/fe_q.h>
+#include <fe/fe_values.h>
+#include <numerics/matrices.h>
+#include <numerics/data_out.h>
+#include <numerics/vectors.h>
+#include <fstream>
+
+#define DIM 2
+
+
+using namespace dealii;
+
+template <int dim>
+class DirichletBoundaryValues : public Function<dim>
+{
+  public:
+    DirichletBoundaryValues() : Function<dim> (2) {};
+    
+    virtual void vector_value (        const Point<dim> &p,
+                                Vector<double>   &values) const;
+    virtual void vector_value_list (const std::vector<Point<dim> > &   points,
+                                    std::vector<Vector<double> > &             value_list) const;
+};
+
+
+template <int dim>
+inline
+void DirichletBoundaryValues<dim>::vector_value (      const Point<dim> &      /*p*/,
+                                                       Vector<double> &                values) const 
+{
+  Assert (values.size() == 2, ExcDimensionMismatch (values.size(), 2));
+       
+  values(0) = 1;
+  values(1) = 0;
+}
+
+
+template <int dim>
+void DirichletBoundaryValues<dim>::vector_value_list (const std::vector<Point<dim> > & points,
+                                                     std::vector<Vector<double> > &            value_list) const 
+{
+  Assert (value_list.size() == points.size(), 
+          ExcDimensionMismatch (value_list.size(), points.size()));
+
+  for (unsigned int p=0; p<points.size(); ++p)
+    DirichletBoundaryValues<dim>::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 &paramhandler)
+               :
+               prm(&paramhandler)
+{}
+
+
+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 <int dim>
+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<dim>                 triangulation;
+    DoFHandler<dim>         dof_handler;
+    FESystem<dim>                                              fe;
+    
+    SparsityPattern         sparsity_pattern;
+    SparseMatrix<double>    system_matrix;
+               
+    Vector<double>          solution, system_rhs, intensity;
+
+    const double                                               c, omega;
+};
+
+
+template <int dim>
+UltrasoundProblem<dim>::UltrasoundProblem (ParameterHandler&   param) 
+               :
+               prm(param),
+               dof_handler(triangulation),
+               fe(FE_Q<dim>(1), 2),
+               c(prm.get_double("c")), 
+               omega(prm.get_double("omega"))
+{}
+
+
+template <int dim>
+UltrasoundProblem<dim>::~UltrasoundProblem () 
+{
+  dof_handler.clear();
+}
+
+
+template <int dim>
+void UltrasoundProblem<dim>::make_grid ()
+{
+  GridGenerator::subdivided_hyper_cube (triangulation, 5, 0, 1);
+       
+  const Point<dim>     transducer = (dim == 2) ? 
+                                    Point<dim> (0.5, 0.0) :
+                                    Point<dim> (0.5, 0.5, 0.0), 
+                      focal_point = (dim == 2) ?
+                                    Point<dim> (0.5, prm.get_double("Focal distance")) :
+                                    Point<dim> (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<dim>::cell_iterator
+    cell = triangulation.begin (),
+    endc = triangulation.end();
+         
+  for (; cell!=endc; ++cell)
+    for (unsigned int face=0; face<GeometryInfo<dim>::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<dim> 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 <int dim>
+void UltrasoundProblem<dim>::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 <int dim>
+void UltrasoundProblem<dim>::assemble_system () 
+{ 
+  const double om2 = omega * omega;
+  const double c2 = c * c;
+       
+  QGauss<dim>          quadrature_formula(3);
+  QGauss<dim-1> 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<double>   cell_matrix (dofs_per_cell, dofs_per_cell);
+
+  std::vector<unsigned int> local_dof_indices (dofs_per_cell);
+  
+  FEValues<dim>  fe_values (fe, quadrature_formula, 
+                           update_values | update_gradients |
+                           update_JxW_values);
+
+  FEFaceValues<dim> fe_face_values (fe, face_quadrature_formula, 
+                                   update_values | update_JxW_values);
+
+  typename DoFHandler<dim>::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; i<dofs_per_cell; ++i)
+       {
+         for (unsigned int j=0; j<dofs_per_cell; ++j)
+           {
+             if (fe.system_to_component_index(i).first == 
+                 fe.system_to_component_index(j).first)
+               {
+                 for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
+                   cell_matrix(i,j) += (((fe_values.shape_value(i,q_point) *
+                                          fe_values.shape_value(j,q_point)) *
+                                         (- om2)
+                                         +
+                                         (fe_values.shape_grad(i,q_point) *
+                                          fe_values.shape_grad(j,q_point)) *
+                                         c2) *
+                                        fe_values.JxW(q_point));
+               }
+           }
+       }
+
+      if (cell->at_boundary())
+       for (unsigned int face=0; face<GeometryInfo<dim>::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; i<dofs_per_cell; ++i)
+               for (unsigned int j=0; j<dofs_per_cell; ++j)
+                 if ((fe.system_to_component_index(i).first != 
+                      fe.system_to_component_index(j).first) &&
+                     fe.has_support_on_face(i, face) &&
+                     fe.has_support_on_face(j, face))
+
+                   for (unsigned int q_point=0; q_point<n_face_q_points; ++q_point)
+                     cell_matrix(i,j) += ((fe.system_to_component_index(i).first) ? 1 : (-1)) * 
+                                         fe_face_values.shape_value(i,q_point) *
+                                         fe_face_values.shape_value(j,q_point) *
+                                         c * omega *
+                                         fe_face_values.JxW(q_point);
+           }
+
+      cell->get_dof_indices (local_dof_indices);
+
+      for (unsigned int i=0; i<dofs_per_cell; ++i)
+       for (unsigned int j=0; j<dofs_per_cell; ++j)
+         system_matrix.add (   local_dof_indices[i],
+                               local_dof_indices[j],
+                               cell_matrix(i,j));
+    }
+
+  std::map<unsigned int,double> boundary_values;
+  VectorTools::interpolate_boundary_values (dof_handler,
+                                           1,
+                                           DirichletBoundaryValues<dim>(),
+                                           boundary_values);
+  MatrixTools::apply_boundary_values (boundary_values,
+                                     system_matrix,
+                                     solution,
+                                     system_rhs);
+}
+
+
+template <int dim>
+void UltrasoundProblem<dim>::solve () 
+{
+  SparseDirectUMFPACK                  A_direct;
+
+  A_direct.initialize(system_matrix);
+  A_direct.vmult(solution,system_rhs);
+}
+
+
+template <int dim>
+void UltrasoundProblem<dim>::postprocess () 
+{
+  QMidpoint<dim>       midpoint_rule;
+  FEValues<dim> fe_values(fe, midpoint_rule, update_values);
+
+  std::vector<Vector<double> >         values;
+  values.resize(1);
+  values[0].reinit(2);
+
+  typename DoFHandler<dim>::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 <int dim>
+void UltrasoundProblem<dim>::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<dim> data_out;
+  data_out.parse_parameters(prm);
+  data_out.attach_dof_handler (dof_handler);
+
+  std::vector<std::string> 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 <int dim>
+void UltrasoundProblem<dim>::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<DIM>   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 (file)
index 0000000..9d50174
--- /dev/null
@@ -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
+
+

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.