From: bangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Date: Mon, 27 Aug 2007 20:16:18 +0000 (+0000)
Subject: add
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=78e811feb062213f3dae5abad9ef89f7db17f773;p=dealii-svn.git

add

git-svn-id: https://svn.dealii.org/trunk@15076 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
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 <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
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
+
+