From fbb6e29445f52d5fc5fcc7fd87914790aa9fa456 Mon Sep 17 00:00:00 2001
From: hartmann <hartmann@0785d39b-7218-0410-832d-ea1e28bc413d>
Date: Fri, 9 Nov 2001 13:14:19 +0000
Subject: [PATCH] Linear transport in non-conservative form.

git-svn-id: https://svn.dealii.org/trunk@5187 0785d39b-7218-0410-832d-ea1e28bc413d
---
 deal.II/examples/step-12/step-12.cc | 716 ++++++++++++++++++++++++++++
 1 file changed, 716 insertions(+)
 create mode 100644 deal.II/examples/step-12/step-12.cc

diff --git a/deal.II/examples/step-12/step-12.cc b/deal.II/examples/step-12/step-12.cc
new file mode 100644
index 0000000000..ff3e5eeaf1
--- /dev/null
+++ b/deal.II/examples/step-12/step-12.cc
@@ -0,0 +1,716 @@
+/* $Id$ */
+/* Author: Ralf Hartmann, University of Heidelberg, 2000 */
+
+				 // The first few files have already
+				 // been covered in previous examples
+				 // and will thus not be further
+				 // commented on.
+#include <base/quadrature_lib.h>
+#include <base/function.h>
+#include <lac/vector.h>
+#include <lac/sparse_matrix.h>
+#include <lac/vector_memory.h>
+#include <grid/tria.h>
+#include <dofs/dof_handler.h>
+#include <grid/grid_generator.h>
+#include <grid/tria_accessor.h>
+#include <grid/tria_iterator.h>
+#include <dofs/dof_accessor.h>
+#include <dofs/dof_tools.h>
+#include <fe/fe_values.h>
+#include <numerics/data_out.h>
+#include <grid/grid_out.h>
+#include <grid/grid_refinement.h>
+#include <numerics/error_estimator.h>
+
+#include <fe/mapping_q1.h>
+#include <fe/fe_dgq.h>
+#include <lac/precondition_block.h>
+#include <lac/solver_richardson.h>
+
+
+#include <fstream>
+
+
+template <int dim>
+class Beta
+{
+  public:
+    Beta () {};
+
+    void value_list (const std::vector<Point<dim> > &points,
+		     std::vector<Point<dim> > &values) const;
+};
+
+
+template <int dim>
+class RHS:  public Function<dim>
+{
+  public:
+    RHS() {};
+    
+    virtual void value_list (const std::vector<Point<dim> > &points,
+			     std::vector<double> &values,
+			     const unsigned int component=0) const;
+};
+
+
+template <int dim>
+class BoundaryFunction:  public Function<dim>
+{
+  public:
+    BoundaryFunction() {};
+    
+    virtual void value_list (const std::vector<Point<dim> > &points,
+			     std::vector<double> &values,
+			     const unsigned int component=0) const;
+};
+
+template <int dim>
+class DGAssembler
+{
+  public:
+    DGAssembler() {};
+
+    void assemble_cell_term(const FEValuesBase<dim>& fe_v,
+			    FullMatrix<double> &cell_matrix,
+			    Vector<double> &cell_vector);
+    
+    void assemble_face_term(const FEFaceValuesBase<dim>& fe_v,
+			    const FEFaceValuesBase<dim>& fe_v_neighbor,
+			    FullMatrix<double> &cell_matrix,
+			    FullMatrix<double> &cell_inflow_matrix,
+			    Vector<double> &cell_vector);
+    
+  private:
+    Beta<dim> beta_function;
+    RHS<dim> rhs_function;
+    BoundaryFunction<dim> boundary_function;
+};
+
+				 // The main class is again almost
+				 // unchanged. Two additions, however,
+				 // are made: we have added the
+				 // ``refine'' function, which is used
+				 // to adaptively refine the grid
+				 // (instead of the global refinement
+				 // in the previous examples), and a
+				 // variable which will hold the
+				 // constraints associated to the
+				 // hanging nodes.
+template <int dim>
+class TransportProblem
+{
+  public:
+    TransportProblem ();
+    ~TransportProblem ();
+
+    void run ();
+    
+  private:
+    void setup_system ();
+    void assemble_system ();
+    void solve ();
+    void refine_grid ();
+    void output_results (const unsigned int cycle) const;
+
+    Triangulation<dim>   triangulation;
+    MappingQ1<dim>       mapping;
+    
+				     // We need a finite element
+				     // again. This time, we will want
+				     // to use quadratic polynomials
+				     // (but this is only specified in
+				     // the constructor):
+    FE_DGQ<dim>          fe;
+    DoFHandler<dim>      dof_handler;
+
+    SparsityPattern      sparsity_pattern;
+    SparseMatrix<double> system_matrix;
+    
+    Vector<double>       solution;
+    Vector<double>       right_hand_side;
+
+    DGAssembler<dim>     dg_assembler;
+};
+
+
+template <>
+void Beta<2>::value_list(const std::vector<Point<2> > &points,
+			 std::vector<Point<2> > &values) const
+{
+  Assert(values.size()==points.size(), ExcDimensionMismatch(values.size(),points.size()));
+  for (unsigned int i=0; i<points.size(); ++i)
+    {
+      const Point<2> &p=points[i];
+      Point<2> &beta=values[i];
+
+      beta(0)=-p(1);
+      beta(1)=p(0);
+      beta/=sqrt(beta*beta);
+    }
+}
+
+
+
+
+template <int dim>
+void RHS<dim>::value_list(const std::vector<Point<dim> > &,
+			  std::vector<double> &values,
+			  const unsigned int) const
+{
+  for (unsigned int i=0; i<values.size(); ++i)
+    values[i]=0;
+}
+
+
+
+
+template <int dim>
+void BoundaryFunction<dim>::value_list(const std::vector<Point<dim> > &points,
+				       std::vector<double> &values,
+				       const unsigned int) const
+{
+  Assert(values.size()==points.size(), ExcDimensionMismatch(values.size(),points.size()));
+  for (unsigned int i=0; i<values.size(); ++i)
+    {
+      if (points[i](0)<0.5)
+	values[i]=1.;
+      else
+	values[i]=0.;
+    }
+}
+
+
+
+
+template <int dim>
+void DGAssembler<dim>::assemble_cell_term(const FEValuesBase<dim>& fe_v,
+					  FullMatrix<double> &cell_matrix,
+					  Vector<double> &cell_vector)
+{
+  const vector<vector<Tensor<1,2> > > &grad_v = fe_v.get_shape_grads ();
+  const FullMatrix<double> &v = fe_v.get_shape_values ();
+  const vector<double> &JxW = fe_v.get_JxW_values ();
+
+  vector<Point<dim> > beta (fe_v.n_quadrature_points);
+  vector<double> rhs (fe_v.n_quadrature_points);
+  
+  beta_function.value_list (fe_v.get_quadrature_points(), beta);
+  rhs_function.value_list (fe_v.get_quadrature_points(), rhs);
+  
+  for (unsigned int point=0; point<fe_v.n_quadrature_points; ++point)
+    for (unsigned int i=0; i<fe_v.dofs_per_cell; ++i) 
+      {
+	for (unsigned int j=0; j<fe_v.dofs_per_cell; ++j)
+	  cell_matrix(i,j) += beta[point]*grad_v[j][point]*
+			      v(i,point) *
+			      JxW[point];
+	
+	cell_vector(i) += rhs[point] *v(i,point) * JxW[point];
+      }
+}
+
+
+template <int dim>
+void DGAssembler<dim>::assemble_face_term(const FEFaceValuesBase<dim>& fe_v,
+					  const FEFaceValuesBase<dim>& fe_v_neighbor,      
+					  FullMatrix<double> &cell_matrix,
+					  FullMatrix<double> &cell_inflow_matrix,
+					  Vector<double> &cell_vector)
+{
+  DoFHandler<dim>::face_iterator face=fe_v.get_face();
+  
+  const FullMatrix<double> &v = fe_v.get_shape_values ();
+  const FullMatrix<double> &v_neighbor = fe_v_neighbor.get_shape_values ();  
+  const vector<double> &JxW = fe_v.get_JxW_values ();
+  const vector<Point<dim> > &normals = fe_v.get_normal_vectors ();
+
+  vector<Point<dim> > beta (fe_v.n_quadrature_points);
+  vector<double> g(fe_v.n_quadrature_points);
+  
+  beta_function.value_list (fe_v.get_quadrature_points(), beta);
+
+  if (face->at_boundary())
+    boundary_function.value_list (fe_v.get_quadrature_points(), g);
+
+  for (unsigned int point=0; point<fe_v.n_quadrature_points; ++point)
+    {
+      double beta_n=beta[point] * normals[point];
+      if (beta_n<0)
+	{
+	  for (unsigned int i=0; i<fe_v.dofs_per_cell; ++i)
+	    {
+	      for (unsigned int j=0; j<fe_v.dofs_per_cell; ++j)
+		cell_matrix(i,j) -= beta_n *
+				    v(j,point) *
+				    v(i,point) *
+				    JxW[point];
+
+	      if (!face->at_boundary())
+		for (unsigned int k=0; k<fe_v_neighbor.dofs_per_cell; ++k)
+		  cell_inflow_matrix(i,k) += beta_n *
+					     v_neighbor(k,point) *
+					     v(i,point) *
+					     JxW[point];
+	      else
+		cell_vector(i) -= beta_n *
+				  g[point] *
+				  v(i,point) *
+				  JxW[point];
+	    }
+	}
+    }
+}
+
+  
+template <int dim>
+TransportProblem<dim>::TransportProblem () :
+                fe (1),
+		dof_handler (triangulation)
+{}
+
+
+template <int dim>
+TransportProblem<dim>::~TransportProblem () 
+{
+  dof_handler.clear ();
+};
+
+
+
+template <int dim>
+void TransportProblem<dim>::setup_system ()
+{
+				   // To distribute degrees of
+				   // freedom, the ``dof_handler''
+				   // variable takes only the finite
+				   // element object. In this case, it
+				   // will distribute four degrees of
+				   // freedom per cell.
+  dof_handler.distribute_dofs (fe);
+
+  sparsity_pattern.reinit (dof_handler.n_dofs(),
+			   dof_handler.n_dofs(),
+			   dof_handler.max_couplings_between_dofs());
+  DoFTools::make_flux_sparsity_pattern (dof_handler, sparsity_pattern);
+  sparsity_pattern.compress();
+  
+  system_matrix.reinit (sparsity_pattern);
+
+  solution.reinit (dof_handler.n_dofs());
+  right_hand_side.reinit (dof_handler.n_dofs());
+};
+
+
+
+template <int dim>
+void TransportProblem<dim>::assemble_system () 
+{
+				   // See Cockburn paper for the proper quadrature.
+  QGauss2<dim>  quadrature;
+  QGauss2<dim-1>  face_quadrature;
+  
+  const unsigned int dofs_per_cell = dof_handler.get_fe().dofs_per_cell;
+  vector<unsigned int> dofs (dofs_per_cell);
+  vector<unsigned int> dofs_neighbor (dofs_per_cell);
+
+  UpdateFlags update_flags = UpdateFlags(update_values
+					 | update_gradients
+					 | update_q_points
+					 | update_JxW_values);
+  
+  UpdateFlags face_update_flags = UpdateFlags(update_values
+					      | update_q_points
+					      | update_JxW_values
+					      | update_normal_vectors);
+   
+
+  FEValues<dim> fe_v (
+    mapping, fe, quadrature, update_flags);
+  FEFaceValues<dim> fe_v_face (
+    mapping, fe, face_quadrature, face_update_flags);
+  FESubfaceValues<dim> fe_v_subface (
+    mapping, fe, face_quadrature, face_update_flags);
+  FEFaceValues<dim> fe_v_face_neighbor (
+    mapping, fe, face_quadrature, UpdateFlags(update_values | update_default));
+  FESubfaceValues<dim> fe_v_subface_neighbor (
+    mapping, fe, face_quadrature, UpdateFlags(update_values | update_default));
+
+				   // includes the u and v terms
+  FullMatrix<double> cell_matrix (dofs_per_cell, dofs_per_cell);
+				   // includes u_hat and v terms
+  FullMatrix<double> cell_inflow_matrix (dofs_per_cell, dofs_per_cell);
+
+  Vector<double>  cell_vector (dofs_per_cell);
+
+  DoFHandler<dim>::active_cell_iterator cell = dof_handler.begin_active(),
+					endc = dof_handler.end();
+  DoFHandler<dim>::face_iterator face;
+  DoFHandler<dim>::cell_iterator neighbor;
+  DoFHandler<dim>::cell_iterator neighbor_child;
+
+  for (;cell!=endc; ++cell) 
+    {
+				       // re-init fe values for this cell
+      fe_v.reinit (cell);
+
+      cell_matrix.clear ();
+      cell_vector.clear ();
+
+      dg_assembler.assemble_cell_term(fe_v,
+				      cell_matrix,
+				      cell_vector);
+      
+      cell->get_dof_indices (dofs);
+
+      for (unsigned int face_no=0; face_no<GeometryInfo<dim>::faces_per_cell; ++face_no)
+	{
+	  face = cell->face(face_no);
+	  
+	  cell_inflow_matrix.clear();
+
+	  if (face->at_boundary())
+	    {
+	      fe_v_face.reinit (cell, face_no);
+
+	      dg_assembler.assemble_face_term(fe_v_face,
+					      fe_v_face,
+					      cell_matrix,
+					      cell_inflow_matrix,
+					      cell_vector);
+	    }
+	  else // if (!face->at_boundary())
+	    {
+	      Assert (cell->neighbor(face_no).state() == valid, ExcInternalError());
+	      neighbor = cell->neighbor(face_no);
+	      
+	      if (face->has_children())  // i.e. neighbor is one level more refined than cell
+		{
+						   // store which number #cell# has in the
+						   // list of neighbors of #neighbor#
+		  const unsigned int neighbor2=cell->neighbor_of_neighbor(face_no);
+		  
+		  
+						   // loop over all subfaces
+		  for (unsigned int subface_no=0; subface_no<GeometryInfo<dim>::subfaces_per_face;
+		       ++subface_no)
+		    {
+						       // get an iterator pointing to the
+						       // cell behind the present subface
+		      neighbor_child = neighbor->child(GeometryInfo<dim>::
+						       child_cell_on_face(neighbor2,subface_no));
+		      Assert (neighbor_child->face(neighbor2) == face->child(subface_no),
+			      ExcInternalError());
+		      Assert (!neighbor_child->has_children(), ExcInternalError());
+
+		      fe_v_subface.reinit (cell, face_no, subface_no);
+		      fe_v_face_neighbor.reinit (neighbor_child, neighbor2);
+	
+		      cell_inflow_matrix.clear();
+	    
+		      dg_assembler.assemble_face_term(fe_v_subface,
+						      fe_v_face_neighbor,
+						      cell_matrix,
+						      cell_inflow_matrix,
+						      cell_vector);
+		  
+						       // get indices of dofs of neighbor_child cell
+		      neighbor_child->get_dof_indices (dofs_neighbor);
+		      						
+						       // distribute cell matrix
+		      for (unsigned int i=0; i<dofs_per_cell; ++i)
+			for (unsigned int k=0; k<dofs_per_cell; ++k)
+			  system_matrix.add(dofs[i], dofs_neighbor[k],
+					    cell_inflow_matrix(i,k));
+		    }
+		}
+	      else // if (!face->has_children())
+		{
+		  if (neighbor->level() == cell->level()) 
+		    {
+						       // store which number #cell# has in the
+						       // list of neighbors of #neighbor#
+		      const unsigned int neighbor2=cell->neighbor_of_neighbor(face_no);
+
+		      fe_v_face.reinit (cell, face_no);
+		      fe_v_face_neighbor.reinit (neighbor, neighbor2);
+		      
+		      dg_assembler.assemble_face_term(fe_v_face,
+						      fe_v_face_neighbor,
+						      cell_matrix,
+						      cell_inflow_matrix,
+						      cell_vector);
+		    }
+		  else // if (neighbor->level() < cell->level()) i.e. neighbor is one level coarser than cell
+		    {
+		      Assert(neighbor->level() < cell->level(), ExcInternalError());
+
+		      const std::pair<unsigned int, unsigned int> faceno_subfaceno=
+			cell->neighbor_of_coarser_neighbor(face_no);
+		      const unsigned int neighbor_face_no=faceno_subfaceno.first,
+				      neighbor_subface_no=faceno_subfaceno.second;
+
+		      Assert (neighbor->neighbor(neighbor_face_no)
+			      ->child(GeometryInfo<dim>::child_cell_on_face(
+				face_no,neighbor_subface_no)) == cell, ExcInternalError());
+			
+						       // now 'neighbor_face_no' stores the number
+						       // of a face in the list of faces of 'neighbor'.
+						       // This face has got a subface that is 
+						       // between 'cell' and 'neighbor'.
+						       // 'neighbor_subface_no' stores the number
+						       // of this subface in the list of subfaces of this
+						       // face 'neighbor->face(neighbor_face_no)'
+						       // that is between 'cell' and 'neighbor'
+		      fe_v_face.reinit (cell, face_no);
+		      fe_v_subface_neighbor.reinit (neighbor, neighbor_face_no,
+						    neighbor_subface_no);
+		      
+		      dg_assembler.assemble_face_term(fe_v_face,
+						      fe_v_subface_neighbor,
+						      cell_matrix,
+						      cell_inflow_matrix,
+						      cell_vector);
+		    } // else // if (neighbor->level() < cell->level())
+
+						   // get indices of dofs of neighbor_child cell
+		  neighbor->get_dof_indices (dofs_neighbor);
+		 						
+						   // distribute cell_inflow_matrix
+		  for (unsigned int i=0; i<dofs_per_cell; ++i)
+		    for (unsigned int k=0; k<dofs_per_cell; ++k)
+		      system_matrix.add(dofs[i], dofs_neighbor[k],
+					cell_inflow_matrix(i,k));
+		} // else // if (!face->has_children())
+	    }  // else // if (!face->at_boundary())
+	} //for (face_no...)
+      
+				       // distribute cell matrix
+      for (unsigned int i=0; i<dofs_per_cell; ++i)
+	for (unsigned int j=0; j<dofs_per_cell; ++j)
+	  system_matrix.add(dofs[i], dofs[j], cell_matrix(i,j));
+      
+				       // distribute cell vector
+      for (unsigned int i=0; i<dofs_per_cell; ++i)
+	right_hand_side(dofs[i]) += cell_vector(i);
+    }  // for (cell...)
+};
+
+
+
+template <int dim>
+void TransportProblem<dim>::solve () 
+{  
+  SolverControl           solver_control (1000, 1e-12);
+  PrimitiveVectorMemory<> vector_memory;
+  SolverRichardson<>      solver (solver_control, vector_memory);
+
+  PreconditionBlockSSOR<double> preconditioner;
+  preconditioner.initialize(system_matrix, fe.dofs_per_cell);
+  preconditioner.invert_diagblocks();
+  
+  solver.solve (system_matrix, solution, right_hand_side,
+		preconditioner);
+};
+
+
+template <int dim>
+void TransportProblem<dim>::refine_grid ()
+{
+  Vector<float> estimated_error_per_cell (triangulation.n_active_cells());
+
+  FunctionMap<dim>::type neumann_boundary;
+
+  KellyErrorEstimator<dim>::estimate (dof_handler,
+				      QGauss3<dim-1>(),
+				      neumann_boundary,
+				      solution,
+				      estimated_error_per_cell);
+
+  GridRefinement::refine_and_coarsen_fixed_number (triangulation,
+						   estimated_error_per_cell,
+						   0.3, 0.03);
+
+  triangulation.execute_coarsening_and_refinement ();
+};
+
+
+
+template <int dim>
+void TransportProblem<dim>::output_results (const unsigned int cycle) const
+{
+				   // We want to write the grid in
+				   // each cycle. Here is another way
+				   // to quickly produce a filename
+				   // based on the cycle number. It
+				   // assumes that the numbers `0'
+				   // through `9' are represented
+				   // consecutively in the character
+				   // set (which is the case in all
+				   // known character sets). However,
+				   // this will only work if the cycle
+				   // number is less than ten, which
+				   // we check by an assertion.
+  std::string filename = "grid-";
+  filename += ('0' + cycle);
+  Assert (cycle < 10, ExcInternalError());
+  
+  filename += ".eps";
+  std::ofstream eps_output (filename.c_str());
+
+				   // Using this filename, we write
+				   // each grid as a postscript file.
+  GridOut grid_out;
+  grid_out.write_eps (triangulation, eps_output);
+
+				   // output of the solution
+  filename = "sol-";
+  filename += ('0' + cycle);
+  Assert (cycle < 10, ExcInternalError());
+  
+  filename += ".gnuplot";
+  std::ofstream gnuplot_output (filename.c_str());
+  
+  DataOut<dim> data_out;
+  data_out.attach_dof_handler (dof_handler);
+  data_out.add_data_vector (solution, "u");
+
+  data_out.build_patches ();
+  
+  data_out.write_gnuplot(gnuplot_output);
+};
+
+
+
+template <int dim>
+void TransportProblem<dim>::run () 
+{
+  for (unsigned int cycle=0; cycle<3; ++cycle)
+    {
+      std::cout << "Cycle " << cycle << ':' << std::endl;
+
+      if (cycle == 0)
+	{
+	  GridGenerator::hyper_cube (triangulation);
+
+	  triangulation.refine_global (3);
+	}
+      else
+					 // In case this is not the
+					 // first cycle, we want to
+					 // refine the grid. Unlike
+					 // the global refinement
+					 // employed in the last
+					 // example, we now use the
+					 // adaptive procedure
+					 // described in the function
+					 // which we now call:
+	{
+	  refine_grid ();
+	};
+      
+
+      std::cout << "   Number of active cells:       "
+		<< triangulation.n_active_cells()
+		<< std::endl;
+
+      setup_system ();
+
+      std::cout << "   Number of degrees of freedom: "
+		<< dof_handler.n_dofs()
+		<< std::endl;
+      
+      assemble_system ();
+      solve ();
+      output_results (cycle);
+    }
+}
+
+int main () 
+{
+
+				   // The general idea behind the
+				   // layout of this function is as
+				   // follows: let's try to run the
+				   // program as we did before...
+  try
+    {
+      TransportProblem<2> Transport_problem_2d;
+      Transport_problem_2d.run ();
+    }
+				   // ...and if this should fail, try
+				   // to gather as much information as
+				   // possible. Specifically, if the
+				   // exception that was thrown is an
+				   // object of a class that is
+				   // derived from the C++ standard
+				   // class ``exception'', then we can
+				   // use the ``what'' member function
+				   // to get a string which describes
+				   // the reason why the exception was
+				   // thrown. 
+				   //
+				   // The deal.II exception classes
+				   // are all derived from the
+				   // standard class, and in
+				   // particular, the ``exc.what()''
+				   // function will return
+				   // approximately the same string as
+				   // would be generated if the
+				   // exception was thrown using the
+				   // ``Assert'' macro. You have seen
+				   // the output of such an exception
+				   // in the previous example, and you
+				   // then know that it contains the
+				   // file and line number of where
+				   // the exception occured, and some
+				   // other information. This is also
+				   // what would be printed in the
+				   // following.
+  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;
+				       // We can't do much more than
+				       // printing as much information
+				       // as we can get to, so abort
+				       // with error:
+      return 1;
+    }
+				   // If the exception that was thrown
+				   // somewhere was not an object of a
+				   // class derived from the standard
+				   // ``exception'' class, then we
+				   // can't do anything at all. We
+				   // then simply print an error
+				   // message and exit.
+  catch (...) 
+    {
+      std::cerr << std::endl << std::endl
+		<< "----------------------------------------------------"
+		<< std::endl;
+      std::cerr << "Unknown exception!" << std::endl
+		<< "Aborting!" << std::endl
+		<< "----------------------------------------------------"
+		<< std::endl;
+      return 1;
+    };
+
+				   // If we got to this point, there
+				   // was no exception which
+				   // propagated up to the main
+				   // function (maybe there were some,
+				   // but they were caught somewhere
+				   // in the program or the
+				   // library). Therefore, the program
+				   // performed as was expected and we
+				   // can return without error.
+  return 0;
+};
-- 
2.39.5