]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Reindent everything.
authorbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Sat, 10 May 2008 23:50:22 +0000 (23:50 +0000)
committerbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Sat, 10 May 2008 23:50:22 +0000 (23:50 +0000)
git-svn-id: https://svn.dealii.org/trunk@16071 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/examples/step-33/step-33.cc

index 7bd3a0e930941b226fd20f41244f6c597f11aa7b..234d66ea6467ced9a10046fb83a3dfa92a1eebd9 100644 (file)
 /*    to the file deal.II/doc/license.html for the  text  and     */
 /*    further information on this license.                        */
 
-                                  // This program solves the Euler equations
-                                  // of gas dynamics for a given configuration
-                                  // file.  It uses a standard Galerkin approach
-                                  // with weakly applied boundary conditions.
-                                  // @sect3{Include files}
+                                // @sect3{Include files}
 
                                  // Here we have the necessary TRILINOS includes.
                                  //
 #include <fstream>
 #include <vector>
 
-                                // Introduce the dealii library into the current namespace.
+                                // Introduce the dealii library into the current namespace.
 using namespace dealii;
 
 #define DIMENSION 2
 
-                                // We define a shorter name for the automatic differentiation
-                                // type.
+                                // We define a shorter name for the automatic differentiation
+                                // type.
 typedef Sacado::Fad::DFad<double> fad_double;
 typedef unsigned int UInt;
-                                // The Epetra library requires a 'communicator', which describes
-                                // the layout of a parallel (or serial) set of processors.
+                                // The Epetra library requires a 'communicator', which describes
+                                // the layout of a parallel (or serial) set of processors.
 Epetra_SerialComm *Comm;
 
-                                // @sect3{Flux function definition}
-                                // Here we define the flux function for this system of conservation
-                                // laws.  Note: it would be terribly difficult to use this example
-                                // to solve some other system of conservation laws.
-                                //
-                                // We define the number of components in the system.  Euler's has
-                                // one entry for momenta in each spatial direction, plus the energy
-                                // and density components.
+                                // @sect3{Flux function definition}
+                                // Here we define the flux function for this system of conservation
+                                // laws.  Note: it would be terribly difficult to use this example
+                                // to solve some other system of conservation laws.
+                                //
+                                // We define the number of components in the system.  Euler's has
+                                // one entry for momenta in each spatial direction, plus the energy
+                                // and density components.
 #define N_COMP (2 + DIMENSION)
-                                // Define a handle to the density and energy indices.  We have arrange
-                                // the momenta to be first, then density, and, lastly, energy.
+                                // Define a handle to the density and energy indices.  We have arrange
+                                // the momenta to be first, then density, and, lastly, energy.
 #define DENS_IDX DIMENSION
 #define ENERGY_IDX (DIMENSION+1)
 
-                                // The gas constant.  This value is representative of air.
+                                // The gas constant.  This value is representative of air.
 const double GAMMA = 1.4;
-                                // We define the flux functions as one large matrix.  Each row of this
-                                // matrix represents a scalar conservation law for the component in
-                                // that row.  We template the numerical type of the flux function
-                                // so that we may use the automatic differentiation type here.
-                                // The flux functions are defined in terms of the
-                                // conserved variables $\rho w_0, \dots, \rho w_{d-1}, \rho, E$,
-                                // so they do not look exactly like the Euler equations one is
-                                // used to seeing.  We evaluate the flux at a single quadrature
-                                // point.
+                                // We define the flux functions as one large matrix.  Each row of this
+                                // matrix represents a scalar conservation law for the component in
+                                // that row.  We template the numerical type of the flux function
+                                // so that we may use the automatic differentiation type here.
+                                // The flux functions are defined in terms of the
+                                // conserved variables $\rho w_0, \dots, \rho w_{d-1}, \rho, E$,
+                                // so they do not look exactly like the Euler equations one is
+                                // used to seeing.  We evaluate the flux at a single quadrature
+                                // point.
 template <typename number, int dim>
 void Flux(std::vector<std::vector<number> >  &flux, 
           const Point<dim> &/*point*/,
           const std::vector<number> &W)
 {
 
-                               // Pressure is a dependent variable: $p = 
-                               // (\gamma - 1)(E-\frac{1}{2} \rho |v|^2)$.
-    number rho_normVsqr;
-    for (unsigned int d0 = 0; d0 < dim; d0++)
-      rho_normVsqr += W[d0]*W[d0];
-                               // Since W are $\rho v$, we get a $\rho^2$ in the
-                               // numerator, so dividing a $\rho$ out gives the desired $ \rho |v|^2$.
-    rho_normVsqr /= W[DENS_IDX];
-
-    number pressure = (GAMMA-1.0)*(W[ENERGY_IDX] - number(0.5)*(rho_normVsqr));
-
-                               // We compute the momentum terms.  We divide by the
-                               // density here to get $v_i \rho v_j$
-    for (unsigned int d = 0; d < dim; d++) {
-      for (unsigned int d1 = 0; d1 < dim; d1++) {
-        flux[d][d1] = W[d]*W[d1]/W[DENS_IDX];
-      }
-                              // The pressure contribution, along the diagonal:
-      flux[d][d] += pressure;
-                              // Advection/conservation of density:
-      flux[DENS_IDX][d] = W[d]; 
-                              // And, lastly, conservation of energy.
-      flux[ENERGY_IDX][d] = W[d]/W[DENS_IDX]*
-                               (W[ENERGY_IDX] + pressure); // energy
+                                  // Pressure is a dependent variable: $p = 
+                                  // (\gamma - 1)(E-\frac{1}{2} \rho |v|^2)$.
+  number rho_normVsqr;
+  for (unsigned int d0 = 0; d0 < dim; d0++)
+    rho_normVsqr += W[d0]*W[d0];
+                                  // Since W are $\rho v$, we get a $\rho^2$ in the
+                                  // numerator, so dividing a $\rho$ out gives the desired $ \rho |v|^2$.
+  rho_normVsqr /= W[DENS_IDX];
+
+  number pressure = (GAMMA-1.0)*(W[ENERGY_IDX] - number(0.5)*(rho_normVsqr));
+
+                                  // We compute the momentum terms.  We divide by the
+                                  // density here to get $v_i \rho v_j$
+  for (unsigned int d = 0; d < dim; d++) {
+    for (unsigned int d1 = 0; d1 < dim; d1++) {
+      flux[d][d1] = W[d]*W[d1]/W[DENS_IDX];
     }
+                                    // The pressure contribution, along the diagonal:
+    flux[d][d] += pressure;
+                                    // Advection/conservation of density:
+    flux[DENS_IDX][d] = W[d]; 
+                                    // And, lastly, conservation of energy.
+    flux[ENERGY_IDX][d] = W[d]/W[DENS_IDX]*
+                         (W[ENERGY_IDX] + pressure); // energy
+  }
 }
 
-                              // On the boundaries of the domain and across `hanging nodes` we use
-                              // a numerical flux function to enforce boundary conditions.  This routine
-                              // is the basic Lax-Friedrich's flux with a stabilization parameter
-                              // $\alpha$.
+                                // On the boundaries of the domain and across `hanging nodes` we use
+                                // a numerical flux function to enforce boundary conditions.  This routine
+                                // is the basic Lax-Friedrich's flux with a stabilization parameter
+                                // $\alpha$.
 template <typename number, int dim>
 void LFNumFlux(
-            std::vector<std::vector<fad_double> > &nflux,
-            const std::vector<Point<dim> > &points, 
-            const std::vector<Point<dim> > &normals,
-            const std::vector<std::vector<number> > &Wplus,
-            const std::vector<std::vector<number> > &Wminus,
-            double alpha)
+  std::vector<std::vector<fad_double> > &nflux,
+  const std::vector<Point<dim> > &points, 
+  const std::vector<Point<dim> > &normals,
+  const std::vector<std::vector<number> > &Wplus,
+  const std::vector<std::vector<number> > &Wminus,
+  double alpha)
 {
   const unsigned int n_q_points = points.size();
 
-                             // We evaluate the flux at each of the quadrature points.
+                                  // We evaluate the flux at each of the quadrature points.
   for (unsigned int q = 0; q < n_q_points; q++) {
     std::vector<std::vector<fad_double> > iflux(N_COMP,
-                                std::vector<fad_double>(dim, 0));
+                                               std::vector<fad_double>(dim, 0));
     std::vector<std::vector<fad_double> > oflux(N_COMP,
-                                std::vector<fad_double>(dim, 0));
+                                               std::vector<fad_double>(dim, 0));
 
     Flux<number, dim>(iflux, points[q], Wplus[q]);
     Flux<number, dim>(oflux, points[q], Wminus[q]);
@@ -181,28 +176,28 @@ void LFNumFlux(
       for (unsigned int d = 0; d < dim; d++) {
         nflux[q][di] += 0.5*(iflux[di][d] + oflux[di][d])*normals[q](d);
       }
-        nflux[q][di] += 0.5*alpha*(Wplus[q][di] - Wminus[q][di]);
+      nflux[q][di] += 0.5*alpha*(Wplus[q][di] - Wminus[q][di]);
     }
   }
 
 }
 
-                            // @sect3{Initial and side condition parsing}
-                            // For the initial condition we use the expression parser function
-                            // object.
+                                // @sect3{Initial and side condition parsing}
+                                // For the initial condition we use the expression parser function
+                                // object.
 template <int dim>
 class InitialCondition :  public FunctionParser<dim> 
 {
   public:
     InitialCondition ();
     
-                            // This function should be called after parsing, but before using
-                            // the object.  It formalizes the expressions and initializes the
-                            // function parser with the appropriate expressions.
+                                    // This function should be called after parsing, but before using
+                                    // the object.  It formalizes the expressions and initializes the
+                                    // function parser with the appropriate expressions.
     void Init();
 
-                            // During parsing we call this function as the initial condition
-                            // for one of the $\mathbf{w}$ variables is encountered.
+                                    // During parsing we call this function as the initial condition
+                                    // for one of the $\mathbf{w}$ variables is encountered.
 
     void set_ic(int _row, std::string &expr) {
       expressions[_row] = expr;
@@ -210,8 +205,8 @@ class InitialCondition :  public FunctionParser<dim>
 
     virtual void vector_value_list (const std::vector<Point<dim> > &points,
                                    std::vector<Vector<double> >   &value_list) const;
- private:
- std::vector<std::string> expressions;
 private:
   std::vector<std::string> expressions;
 };
 
 template <int dim>
@@ -220,21 +215,21 @@ InitialCondition<dim>::InitialCondition () :
                 expressions(N_COMP, "0.0")
 {}
 
-                            // Here we set up x,y,z as the variables that one should use in the input
-                            // deck to describe their initial condition.
+                                // Here we set up x,y,z as the variables that one should use in the input
+                                // deck to describe their initial condition.
 template<int dim>
 void InitialCondition<dim>::Init() {
- std::map<std::string, double> constants;
- constants["M_PI"] =  M_PI;
- std::string variables = (dim == 2 ? "x,y" : "x,y,z");
 std::map<std::string, double> constants;
 constants["M_PI"] =  M_PI;
 std::string variables = (dim == 2 ? "x,y" : "x,y,z");
 
- FunctionParser<dim>::initialize(variables, expressions, constants);
 FunctionParser<dim>::initialize(variables, expressions, constants);
 
 }
 
 template <int dim>
 void InitialCondition<dim>::vector_value_list (const std::vector<Point<dim> > &points,
-                                           std::vector<Vector<double> >   &value_list) const 
+                                              std::vector<Vector<double> >   &value_list) const 
 {
   const unsigned int n_points = points.size();
 
@@ -243,10 +238,10 @@ void InitialCondition<dim>::vector_value_list (const std::vector<Point<dim> > &p
 
   for (unsigned int p=0; p<n_points; ++p)
     InitialCondition<dim>::vector_value (points[p],
-                                     value_list[p]);
+                                        value_list[p]);
 }
 
-                             // As above, we use the expression function parser for boundary conditions.
+                                // As above, we use the expression function parser for boundary conditions.
 template <int dim>
 class SideCondition :  public FunctionParser<dim> 
 {
@@ -254,15 +249,15 @@ class SideCondition :  public FunctionParser<dim>
     SideCondition (int ncomp);
     ~SideCondition ();
 
-                            // As above.
+                                    // As above.
     void Init();
-                           // As above.
+                                    // As above.
     void set_coeff_row(int _row_n, std::string &expr);
     
     virtual void vector_value_list (const std::vector<Point<dim> > &points,
                                    std::vector<Vector<double> >   &value_list) const;
   private:
-  std::vector<std::string> expressions;
+    std::vector<std::string> expressions;
 };
 
 template <int dim>
@@ -274,16 +269,16 @@ SideCondition<dim>::SideCondition (int ncomp) :
 template <int dim>
 void SideCondition<dim>::set_coeff_row (int _row_n, std::string &expr) 
 {
- expressions[_row_n] = expr;
 expressions[_row_n] = expr;
 }
 
 template <int dim>
 void SideCondition<dim>::Init() {
- std::map<std::string, double> constants;
- constants["M_PI"] =  M_PI;
- std::string variables = (dim == 2 ? "x,y" : "x,y,z");
 std::map<std::string, double> constants;
 constants["M_PI"] =  M_PI;
 std::string variables = (dim == 2 ? "x,y" : "x,y,z");
 
- FunctionParser<dim>::initialize(variables, expressions, constants);
 FunctionParser<dim>::initialize(variables, expressions, constants);
 
 }
 
@@ -306,9 +301,9 @@ void SideCondition<dim>::vector_value_list (const std::vector<Point<dim> > &poin
                                      value_list[p]);
 }
                                 // @sect3{Conservation Law class}
-                           // Here we define a Conservation Law class that helps group
-                           // operations and data for our Euler equations into a manageable
-                           // entity.  Functions will be described as their definitions appear.
+                                // Here we define a Conservation Law class that helps group
+                                // operations and data for our Euler equations into a manageable
+                                // entity.  Functions will be described as their definitions appear.
 template <int dim>
 class ConsLaw
 {
@@ -364,138 +359,138 @@ class ConsLaw
     void assemble_cell_term(const FEValues<dim>& fe_v,
                             std::vector<unsigned int> &dofs,
                             unsigned int cell_no
-                            );
+    );
     
     void assemble_face_term(
-                            int face_no,
-                            const FEFaceValuesBase<dim>& fe_v,
-                            const FEFaceValuesBase<dim>& fe_v_neighbor,
-                             std::vector<unsigned int> &dofs,
-                             std::vector<unsigned int> &dofs_neighbor,
-                             int boundary = -1
-                             );
+      int face_no,
+      const FEFaceValuesBase<dim>& fe_v,
+      const FEFaceValuesBase<dim>& fe_v_neighbor,
+      std::vector<unsigned int> &dofs,
+      std::vector<unsigned int> &dofs_neighbor,
+      int boundary = -1
+    );
 
     unsigned int get_n_components() const { return N_COMP;}
 
   private:
-                                    // T = current time, dT = time step, TF = final time.
+                                    // T = current time, dT = time step, TF = final time.
     double T, dT, TF;
     double face_diameter;
     double cell_diameter;
-                                    // An object to handle parsing the input deck.
+                                    // An object to handle parsing the input deck.
     ParameterHandler prm;
-                                    // Name of the mesh to read in.
+                                    // Name of the mesh to read in.
     string mesh;
     InitialCondition<dim> ic;
 
-                                    // Enums for the various supported boundary conditions.
+                                    // Enums for the various supported boundary conditions.
     typedef enum {INFLOW_BC = 1, OUTFLOW_BC=2, NO_PENETRATION_BC=3, PRESSURE_BC=4} bc_type;
 
-                                    // For each boundary we store a map from boundary # to the type
-                                    // of boundary condition.  If the boundary condition is prescribed,
-                                    // we store a pointer to a function object that will hold the expression
-                                    // for that boundary condition.
+                                    // For each boundary we store a map from boundary # to the type
+                                    // of boundary condition.  If the boundary condition is prescribed,
+                                    // we store a pointer to a function object that will hold the expression
+                                    // for that boundary condition.
     typedef typename std::map<unsigned int, std::pair<std::vector<bc_type>, Function<dim>*> > bdry_map_type;
     bdry_map_type bdry_map;
 
     void add_boundary(unsigned int bd, std::vector<bc_type>& flags, Function<dim> *bf);
 
-                                    // An object to store parameter information about the Aztec solver.
+                                    // An object to store parameter information about the Aztec solver.
     typedef struct {
-      int LIN_OUTPUT;
-      typedef enum { GMRES = 0, DIRECT = 1} solver_type;
-      solver_type SOLVER;
-      typedef enum { QUIET = 0, VERBOSE = 1 } output_type;
-      output_type OUTPUT;
-                                    // Linear residual tolerance.
-      double RES;
-      int MAX_ITERS;
-                                    // We use the ILUT preconditioner.  This is similar
-                                    // to the ILU.  FILL is the number of extra entries
-                                    // to add when forming the ILU decomposition.
-      double ILUT_FILL;
-                                    // When forming the preconditioner, for certain problems
-                                    // bad conditioning (or just bad luck) can cause the
-                                    // preconditioner to be very poorly conditioned.  Hence
-                                    // it can help to add diagonal perturbations to the
-                                    // original matrix and form the preconditioner for this
-                                    // slightly better matrix.  ATOL is an absolute perturbation
-                                    // that is added to the diagonal before forming the
-                                    // prec, and RTOL is a scaling factor $rtol >= 1$.
-      double ILUT_ATOL;
-      double ILUT_RTOL;
-                                    // The ILUT will drop any values that have magnitude less
-                                    // than this value.  This is a way to
-                                    // manage the amount of memory used by this preconditioner.
-      double ILUT_DROP;
+       int LIN_OUTPUT;
+       typedef enum { GMRES = 0, DIRECT = 1} solver_type;
+       solver_type SOLVER;
+       typedef enum { QUIET = 0, VERBOSE = 1 } output_type;
+       output_type OUTPUT;
+                                        // Linear residual tolerance.
+       double RES;
+       int MAX_ITERS;
+                                        // We use the ILUT preconditioner.  This is similar
+                                        // to the ILU.  FILL is the number of extra entries
+                                        // to add when forming the ILU decomposition.
+       double ILUT_FILL;
+                                        // When forming the preconditioner, for certain problems
+                                        // bad conditioning (or just bad luck) can cause the
+                                        // preconditioner to be very poorly conditioned.  Hence
+                                        // it can help to add diagonal perturbations to the
+                                        // original matrix and form the preconditioner for this
+                                        // slightly better matrix.  ATOL is an absolute perturbation
+                                        // that is added to the diagonal before forming the
+                                        // prec, and RTOL is a scaling factor $rtol >= 1$.
+       double ILUT_ATOL;
+       double ILUT_RTOL;
+                                        // The ILUT will drop any values that have magnitude less
+                                        // than this value.  This is a way to
+                                        // manage the amount of memory used by this preconditioner.
+       double ILUT_DROP;
     } solver_params_type;
 
     solver_params_type solver_params;
 
-                                    // Some refinement parameters.
+                                    // Some refinement parameters.
     typedef struct {
-      typedef enum { NONE = 0, FIXED_NUMBER = 1, SHOCK = 2} refine_type;
-      double high_frac;
-      double low_frac;
-      refine_type refine;
-      double high_frac_sav;
-      double max_cells;
-      double shock_val;
-      double shock_levels;
+       typedef enum { NONE = 0, FIXED_NUMBER = 1, SHOCK = 2} refine_type;
+       double high_frac;
+       double low_frac;
+       refine_type refine;
+       double high_frac_sav;
+       double max_cells;
+       double shock_val;
+       double shock_levels;
     } refinement_params_type;
 
     refinement_params_type refinement_params;
 
-                                    // The user can set the stabilization parameter $\alpha$ 
-                                    // in the Lax-Friedrich's flux.  
+                                    // The user can set the stabilization parameter $\alpha$ 
+                                    // in the Lax-Friedrich's flux.  
     typedef struct {
-      typedef enum {CONSTANT=1,MESH=2} LF_stab_type;
-      LF_stab_type LF_stab;
-      double LF_stab_value;
+       typedef enum {CONSTANT=1,MESH=2} LF_stab_type;
+       LF_stab_type LF_stab;
+       double LF_stab_value;
     } flux_params_type;
 
     flux_params_type flux_params;
 
     bool is_stationary;
 
-                                    // Power for the mesh stabilization term.
+                                    // Power for the mesh stabilization term.
     double diffusion_power;
     double gravity;
-                                    // If true, we output the squared gradient of the
-                                    // density instead of density.  Using this one can
-                                    // create shock plots.
+                                    // If true, we output the squared gradient of the
+                                    // density instead of density.  Using this one can
+                                    // create shock plots.
     bool schlieren_plot;
-                                    // How often to create an output file.
+                                    // How often to create an output file.
     double output_step;
 
     Epetra_Map         *Map;
     Epetra_CrsMatrix   *Matrix;
     Vector<double>      indicator;
  
-                                   // Crank-Nicolson value
+                                    // Crank-Nicolson value
     const double        theta; 
 
 };
 
 
-                                    // Asign a row of the conservation law a specified
-                                    // boundary type and (possibly) function.
+                                // Asign a row of the conservation law a specified
+                                // boundary type and (possibly) function.
 template <int dim>
 void ConsLaw<dim>::add_boundary(unsigned int bd,
-        std::vector<bc_type> &flags, Function<dim> *bf) {
+                               std::vector<bc_type> &flags, Function<dim> *bf) {
 
   std::pair<std::vector<bc_type>, Function<dim> *> entry(flags, bf);
   bdry_map[bd] = entry;
 }
 
 
-                                    // Apply the initialial condition.  Simultaneously
-                                    // initialize the non-linear solution.
+                                // Apply the initialial condition.  Simultaneously
+                                // initialize the non-linear solution.
 template <int dim>
 void ConsLaw<dim>::initialize() {
- VectorTools::interpolate(dof_handler,
 VectorTools::interpolate(dof_handler,
                            ic, solution);
- VectorTools::interpolate(dof_handler,
 VectorTools::interpolate(dof_handler,
                            ic, nlsolution);
 
 }
@@ -511,63 +506,63 @@ void ConsLaw<dim>::assemble_cell_term(
   const FEValues<dim> &fe_v,
   std::vector<unsigned int> &dofs,
   unsigned int /*cell_no*/
-  
+) 
 {
-                                 // The residual for each row (i) will be accumulating 
-                                 // into this fad variable.  At the end of the assembly
-                                 // for this row, we will query for the sensitivities
-                                 // to this variable and add them into the Jacobian.
+                                  // The residual for each row (i) will be accumulating 
+                                  // into this fad variable.  At the end of the assembly
+                                  // for this row, we will query for the sensitivities
+                                  // to this variable and add them into the Jacobian.
   fad_double F_i;
   unsigned int dofs_per_cell = fe_v.dofs_per_cell;
   unsigned int n_q_points = fe_v.n_quadrature_points;
 
-                                 // We will define the dofs on this cell in these fad variables.
+                                  // We will define the dofs on this cell in these fad variables.
   std::vector<fad_double> DOF(dofs_per_cell);
 
-                                 // Values of the conservative variables at the quadrature points.
+                                  // Values of the conservative variables at the quadrature points.
   std::vector<std::vector<fad_double > > W (n_q_points,
-                                               std::vector<fad_double >(get_n_components()));
+                                           std::vector<fad_double >(get_n_components()));
 
-                                 // Values at the last time step of the conservative variables.
-                                 // Note that these do not use fad variables, since they do
-                                 // not depend on the 'variables to be sought'=DOFS.
+                                  // Values at the last time step of the conservative variables.
+                                  // Note that these do not use fad variables, since they do
+                                  // not depend on the 'variables to be sought'=DOFS.
   std::vector<std::vector<double > > Wl (n_q_points,
-                                               std::vector<double >(get_n_components()));
+                                        std::vector<double >(get_n_components()));
 
-                                 // Here we will hold the averaged values of the conservative
-                                 // variables that we will linearize around (cn=Crank Nicholson).
+                                  // Here we will hold the averaged values of the conservative
+                                  // variables that we will linearize around (cn=Crank Nicholson).
   std::vector<std::vector<fad_double > > Wcn (n_q_points,
-                                               std::vector<fad_double >(get_n_components()));
+                                             std::vector<fad_double >(get_n_components()));
 
-                                 // Gradients of the current variables.  It is a
-                                 // bit of a shame that we have to compute these; we almost don't.
-                                 // The nice thing about a simple conservation law is that the
-                                 // the flux doesn't generally involve any gradients.  We do
-                                 // need these, however, for the diffusion stabilization. 
-   std::vector<std::vector<std::vector<fad_double> > > Wgrads (n_q_points,
-            std::vector<std::vector<fad_double> >(get_n_components(),
-                            std::vector<fad_double>(dim)));
+                                  // Gradients of the current variables.  It is a
+                                  // bit of a shame that we have to compute these; we almost don't.
+                                  // The nice thing about a simple conservation law is that the
+                                  // the flux doesn't generally involve any gradients.  We do
+                                  // need these, however, for the diffusion stabilization. 
+  std::vector<std::vector<std::vector<fad_double> > > Wgrads (n_q_points,
+                                                             std::vector<std::vector<fad_double> >(get_n_components(),
+                                                                                                   std::vector<fad_double>(dim)));
 
 
   const std::vector<double> &JxW = fe_v.get_JxW_values ();
 
   
-                                  // Here is the magical point where we declare a subset
-                                  // of the fad variables as degrees of freedom.  All 
-                                  // calculations that reference these variables (either
-                                  // directly or indirectly) will accumulate sensitivies
-                                  // with respect to these dofs.
+                                  // Here is the magical point where we declare a subset
+                                  // of the fad variables as degrees of freedom.  All 
+                                  // calculations that reference these variables (either
+                                  // directly or indirectly) will accumulate sensitivies
+                                  // with respect to these dofs.
   for (unsigned int in = 0; in < dofs_per_cell; in++) {
-      DOF[in] = nlsolution(dofs[in]);
-      DOF[in].diff(in, dofs_per_cell);
+    DOF[in] = nlsolution(dofs[in]);
+    DOF[in].diff(in, dofs_per_cell);
   }
 
-                                  // Here we compute the shape function values and gradients
-                                  // at the quadrature points.  Ideally, we could call into 
-                                  // something like get_function_values, get_function_grads,
-                                  // but since we don't want to make the entire solution vector
-                                  // fad types, only the local cell variables, we explicitly
-                                  // code this loop;
+                                  // Here we compute the shape function values and gradients
+                                  // at the quadrature points.  Ideally, we could call into 
+                                  // something like get_function_values, get_function_grads,
+                                  // but since we don't want to make the entire solution vector
+                                  // fad types, only the local cell variables, we explicitly
+                                  // code this loop;
   for (unsigned int q = 0; q < n_q_points; q++) {
     for (unsigned int di = 0; di < get_n_components(); di++) {
       W[q][di] = 0;
@@ -578,18 +573,18 @@ void ConsLaw<dim>::assemble_cell_term(
       }
     }
     for (unsigned int sf = 0; sf < dofs_per_cell; sf++) {
-     int di = fe_v.get_fe().system_to_component_index(sf).first;
-     W[q][di] +=
-                DOF[sf]*fe_v.shape_value_component(sf, q, di);
-     Wl[q][di] +=
-                solution(dofs[sf])*fe_v.shape_value_component(sf, q, di);
-     Wcn[q][di] +=
-                (theta*DOF[sf]+(1-theta)*solution(dofs[sf]))*fe_v.shape_value_component(sf, q, di);
-
-     for (unsigned int d = 0; d < dim; d++) {
-       Wgrads[q][di][d] += DOF[sf]*
-                 fe_v.shape_grad_component(sf, q, di)[d];
-     } // for d
+      int di = fe_v.get_fe().system_to_component_index(sf).first;
+      W[q][di] +=
+       DOF[sf]*fe_v.shape_value_component(sf, q, di);
+      Wl[q][di] +=
+       solution(dofs[sf])*fe_v.shape_value_component(sf, q, di);
+      Wcn[q][di] +=
+       (theta*DOF[sf]+(1-theta)*solution(dofs[sf]))*fe_v.shape_value_component(sf, q, di);
+
+      for (unsigned int d = 0; d < dim; d++) {
+       Wgrads[q][di][d] += DOF[sf]*
+                           fe_v.shape_grad_component(sf, q, di)[d];
+      } // for d
 
     }
 
@@ -601,90 +596,90 @@ void ConsLaw<dim>::assemble_cell_term(
                                    // this could be done in a better way, since this
                                    // could be a rather large object, but for now it 
                                    // seems to work just fine.
-   std::vector<std::vector<std::vector<fad_double> > > flux(n_q_points, 
-                                std::vector<std::vector<fad_double> >(get_n_components(),
-                                std::vector<fad_double>(dim, 0)));
+  std::vector<std::vector<std::vector<fad_double> > > flux(n_q_points, 
+                                                          std::vector<std::vector<fad_double> >(get_n_components(),
+                                                                                                std::vector<fad_double>(dim, 0)));
 
-    for (unsigned int q=0; q < n_q_points; ++q) {
-      Flux<fad_double, dim>(flux[q], fe_v.get_quadrature_points()[q], Wcn[q]);
-    }
+  for (unsigned int q=0; q < n_q_points; ++q) {
+    Flux<fad_double, dim>(flux[q], fe_v.get_quadrature_points()[q], Wcn[q]);
+  }
 
-                                    // We now have all of the function values/grads/fluxes,
-                                    // so perform the assembly.  We have an outer loop
-                                    // through the components of the system, and an
-                                    // inner loop over the quadrature points, where we
-                                    // accumulate contributions to the ith residual.
-                                    //
-                                    // We initialy sum all contributions of the residual
-                                    // in the positive sense, so that we don't need to
-                                    // negative the Jacobian entries.  Then, when we sum
-                                    // into the <code> right_hand_side </code> vector,
-                                    // we negate this residual.
-    for (unsigned int i=0; i<fe_v.dofs_per_cell; ++i) 
-      {
-                                    // Find which component this dof contributes to.
-        const unsigned int
-          component_i = fe_v.get_fe().system_to_component_index(i).first;
+                                  // We now have all of the function values/grads/fluxes,
+                                  // so perform the assembly.  We have an outer loop
+                                  // through the components of the system, and an
+                                  // inner loop over the quadrature points, where we
+                                  // accumulate contributions to the ith residual.
+                                  //
+                                  // We initialy sum all contributions of the residual
+                                  // in the positive sense, so that we don't need to
+                                  // negative the Jacobian entries.  Then, when we sum
+                                  // into the <code> right_hand_side </code> vector,
+                                  // we negate this residual.
+  for (unsigned int i=0; i<fe_v.dofs_per_cell; ++i) 
+    {
+                                      // Find which component this dof contributes to.
+      const unsigned int
+       component_i = fe_v.get_fe().system_to_component_index(i).first;
 
-                                    // Initialize the fad residual to zero (removes
-                                    // any previous sensitivities.
-        F_i = 0;
+                                      // Initialize the fad residual to zero (removes
+                                      // any previous sensitivities.
+      F_i = 0;
 
-                                    // Loop quadrature points.
-        for (unsigned int point=0; point<fe_v.n_quadrature_points; ++point) {
+                                      // Loop quadrature points.
+      for (unsigned int point=0; point<fe_v.n_quadrature_points; ++point) {
 
-          fad_double fdotgv = 0;
+       fad_double fdotgv = 0;
 
-                                    // Integrate the flux times gradient of the test function
-          for (unsigned int d = 0; d < dim; d++) 
-            fdotgv += flux[point][component_i][d]*fe_v.shape_grad_component(i, point, component_i)[d];
+                                        // Integrate the flux times gradient of the test function
+       for (unsigned int d = 0; d < dim; d++) 
+         fdotgv += flux[point][component_i][d]*fe_v.shape_grad_component(i, point, component_i)[d];
            
-          F_i -= fdotgv*JxW[point];
+       F_i -= fdotgv*JxW[point];
 
-                                    // The mass term (if the simulation is non-stationary).
-          fad_double delta_t= 1.0/dT*(W[point][component_i] - Wl[point][component_i]);
-          if (!is_stationary) F_i += delta_t*
-                 fe_v.shape_value_component(i, point, component_i)*JxW[point];
+                                        // The mass term (if the simulation is non-stationary).
+       fad_double delta_t= 1.0/dT*(W[point][component_i] - Wl[point][component_i]);
+       if (!is_stationary) F_i += delta_t*
+                                  fe_v.shape_value_component(i, point, component_i)*JxW[point];
 
-                                   // Stabilization (cell wise diffusion)
-          fad_double guv = 0;
-          for (unsigned int d = 0; d < dim; d++) {
-            guv += fe_v.shape_grad_component(i, point, component_i)[d]*
-                      Wgrads[point][component_i][d];
-          }
+                                        // Stabilization (cell wise diffusion)
+       fad_double guv = 0;
+       for (unsigned int d = 0; d < dim; d++) {
+         guv += fe_v.shape_grad_component(i, point, component_i)[d]*
+                Wgrads[point][component_i][d];
+       }
 
-            F_i += 1.0*std::pow(cell_diameter, diffusion_power)*guv*JxW[point];
+       F_i += 1.0*std::pow(cell_diameter, diffusion_power)*guv*JxW[point];
           
-                                   // The gravity component only enters into the energy 
-                                   // equation and into the vertical component of the 
-                                   // velocity.
-          if (component_i == dim - 1) {
-            F_i += gravity*Wcn[point][DENS_IDX]*fe_v.shape_value_component(i,point, component_i)*JxW[point];
-          } else if (component_i == ENERGY_IDX) {
-            F_i += gravity*Wcn[point][DENS_IDX]*Wcn[point][dim-1]*
-                   fe_v.shape_value_component(i,point, component_i)*JxW[point];
-          }
-        } // for q
-
-                                   // Here we gain access to the array of sensitivities
-                                   // of the residual.  We then sum these into the
-                                   // Epetra matrix.
-        double *values = &(F_i.fastAccessDx(0));
-        Matrix->SumIntoGlobalValues(dofs[i],
-            dofs_per_cell, &values[0], reinterpret_cast<int*>(&dofs[0]));
+                                        // The gravity component only enters into the energy 
+                                        // equation and into the vertical component of the 
+                                        // velocity.
+       if (component_i == dim - 1) {
+         F_i += gravity*Wcn[point][DENS_IDX]*fe_v.shape_value_component(i,point, component_i)*JxW[point];
+       } else if (component_i == ENERGY_IDX) {
+         F_i += gravity*Wcn[point][DENS_IDX]*Wcn[point][dim-1]*
+                fe_v.shape_value_component(i,point, component_i)*JxW[point];
+       }
+      } // for q
+
+                                      // Here we gain access to the array of sensitivities
+                                      // of the residual.  We then sum these into the
+                                      // Epetra matrix.
+      double *values = &(F_i.fastAccessDx(0));
+      Matrix->SumIntoGlobalValues(dofs[i],
+                                 dofs_per_cell, &values[0], reinterpret_cast<int*>(&dofs[0]));
  
-                                   // Add minus the residual to the right hand side.
-        right_hand_side(dofs[i]) -= F_i.val();
+                                      // Add minus the residual to the right hand side.
+      right_hand_side(dofs[i]) -= F_i.val();
 
-      } // for i
+    } // for i
 }
-                                   // @sect4{%Function: assemble_face_term}
-                                   // These are either
-                                   // boundary terms or terms across differing 
-                                   // levels of refinement.  In the first case,
-                                   // fe_v==fe_v_neighbor and dofs==dofs_neighbor.
-                                   // The int boundary < 0 if not at a boundary,
-                                   // otherwise it is the boundary indicator.
+                                // @sect4{%Function: assemble_face_term}
+                                // These are either
+                                // boundary terms or terms across differing 
+                                // levels of refinement.  In the first case,
+                                // fe_v==fe_v_neighbor and dofs==dofs_neighbor.
+                                // The int boundary < 0 if not at a boundary,
+                                // otherwise it is the boundary indicator.
 template <int dim>
 void ConsLaw<dim>::assemble_face_term(
   int face_no,
@@ -693,7 +688,7 @@ void ConsLaw<dim>::assemble_face_term(
   std::vector<unsigned int> &dofs,
   std::vector<unsigned int> &dofs_neighbor,
   int boundary
-  
+) 
 {
   fad_double F_i;
   const unsigned int n_q_points = fe_v.n_quadrature_points;
@@ -702,101 +697,101 @@ void ConsLaw<dim>::assemble_face_term(
   Assert(dofs_per_cell == ndofs_per_cell,
         ExcDimensionMismatch(dofs_per_cell, ndofs_per_cell));
 
-                                  // As above, the fad degrees of freedom
+                                  // As above, the fad degrees of freedom
   std::vector<fad_double> DOF(dofs_per_cell+ndofs_per_cell);
 
-                                  // The conservative variables for this cell,
-                                  // and for 
+                                  // The conservative variables for this cell,
+                                  // and for 
   std::vector<std::vector<fad_double > > Wplus (n_q_points,
-                                               std::vector<fad_double >(get_n_components()));
+                                               std::vector<fad_double >(get_n_components()));
   std::vector<std::vector<fad_double > > Wminus (n_q_points,
-                                               std::vector<fad_double >(get_n_components()));
+                                                std::vector<fad_double >(get_n_components()));
 
 
   const std::vector<double> &JxW = fe_v.get_JxW_values ();
   const std::vector<Point<dim> > &normals = fe_v.get_normal_vectors ();
 
 
-                                  // If we are at a boundary, then dofs_neighbor are
-                                  // the same as dofs, so we do not want to duplicate them.
-                                  // If there is a neighbor cell, then we want to include 
-                                  // them.
+                                  // If we are at a boundary, then dofs_neighbor are
+                                  // the same as dofs, so we do not want to duplicate them.
+                                  // If there is a neighbor cell, then we want to include 
+                                  // them.
   int ndofs = (boundary < 0 ? dofs_per_cell + ndofs_per_cell : dofs_per_cell);
-                                  // Set the local DOFS.
+                                  // Set the local DOFS.
   for (unsigned int in = 0; in < dofs_per_cell; in++) {
-      DOF[in] = nlsolution(dofs[in]);
-      DOF[in].diff(in, ndofs);
+    DOF[in] = nlsolution(dofs[in]);
+    DOF[in].diff(in, ndofs);
   }
-                                  // If present, set the neighbor dofs.
+                                  // If present, set the neighbor dofs.
   if (boundary < 0)
-  for (unsigned int in = 0; in < ndofs_per_cell; in++) {
+    for (unsigned int in = 0; in < ndofs_per_cell; in++) {
       DOF[in+dofs_per_cell] = nlsolution(dofs_neighbor[in]);
       DOF[in+dofs_per_cell].diff(in+dofs_per_cell, ndofs);
-  }
+    }
 
-                                  // Set the values of the local conservative variables.
-                                  // Initialize all variables to zero.
+                                  // Set the values of the local conservative variables.
+                                  // Initialize all variables to zero.
   for (unsigned int q = 0; q < n_q_points; q++) {
     for (unsigned int di = 0; di < get_n_components(); di++) {
-           Wplus[q][di] = 0;
-           Wminus[q][di] = 0;
+      Wplus[q][di] = 0;
+      Wminus[q][di] = 0;
     }
     for (unsigned int sf = 0; sf < dofs_per_cell; sf++) {
-     int di = fe_v.get_fe().system_to_component_index(sf).first;
-     Wplus[q][di] +=
-                (theta*DOF[sf]+(1.0-theta)*solution(dofs[sf]))*fe_v.shape_value_component(sf, q, di);
+      int di = fe_v.get_fe().system_to_component_index(sf).first;
+      Wplus[q][di] +=
+       (theta*DOF[sf]+(1.0-theta)*solution(dofs[sf]))*fe_v.shape_value_component(sf, q, di);
     }
 
 
-                                 // If there is a cell across, then initialize
-                                 // the exterior trace as a function of the other
-                                 // cell degrees of freedom.
+                                    // If there is a cell across, then initialize
+                                    // the exterior trace as a function of the other
+                                    // cell degrees of freedom.
     if (boundary < 0) {
       for (unsigned int sf = 0; sf < ndofs_per_cell; sf++) {
-       int di = fe_v_neighbor.get_fe().system_to_component_index(sf).first;
-       Wminus[q][di] +=
-                (theta*DOF[sf+dofs_per_cell]+(1.0-theta)*solution(dofs_neighbor[sf]))*
-                fe_v_neighbor.shape_value_component(sf, q, di);
+       int di = fe_v_neighbor.get_fe().system_to_component_index(sf).first;
+       Wminus[q][di] +=
+         (theta*DOF[sf+dofs_per_cell]+(1.0-theta)*solution(dofs_neighbor[sf]))*
+         fe_v_neighbor.shape_value_component(sf, q, di);
       }
     } 
-   } // for q
-
-                               // If this is a boundary, then the values of $W^-$ will
-                               // be either functions of $W^+$, or they will be prescribed.
-                               // This switch sets them appropriately.  Since we are
-                               // using fad variables here, sensitivities will be updated 
-                               // appropriately.  These sensitivities would be tremendously
-                               // difficult to manage without fad!!!
-   if (boundary >= 0) {
-                               // Get the boundary descriptor.
-     typename bdry_map_type::iterator bme = bdry_map.find(boundary);
-     assert(bme != bdry_map.end());
-
-                             // Evaluate the function object.  This is a bit
-                             // tricky; a given boundary might have both prescribed
-                             // and implicit values.  If a particular component is not
-                             // prescribed, the values evaluate to zero and are
-                             // ignored, below.
-     std::vector<Vector<double> > bvals(n_q_points, Vector<double>(N_COMP));
-     bme->second.second->vector_value_list(fe_v.get_quadrature_points(), bvals);
-
-                             // We loop the quadrature points, and we treat each
-                             // component individualy.
-     for (unsigned int q = 0; q < n_q_points; q++) {
+  } // for q
+
+                                  // If this is a boundary, then the values of $W^-$ will
+                                  // be either functions of $W^+$, or they will be prescribed.
+                                  // This switch sets them appropriately.  Since we are
+                                  // using fad variables here, sensitivities will be updated 
+                                  // appropriately.  These sensitivities would be tremendously
+                                  // difficult to manage without fad!!!
+  if (boundary >= 0) {
+                                    // Get the boundary descriptor.
+    typename bdry_map_type::iterator bme = bdry_map.find(boundary);
+    assert(bme != bdry_map.end());
+
+                                    // Evaluate the function object.  This is a bit
+                                    // tricky; a given boundary might have both prescribed
+                                    // and implicit values.  If a particular component is not
+                                    // prescribed, the values evaluate to zero and are
+                                    // ignored, below.
+    std::vector<Vector<double> > bvals(n_q_points, Vector<double>(N_COMP));
+    bme->second.second->vector_value_list(fe_v.get_quadrature_points(), bvals);
+
+                                    // We loop the quadrature points, and we treat each
+                                    // component individualy.
+    for (unsigned int q = 0; q < n_q_points; q++) {
       for (unsigned int di = 0; di < get_n_components(); di++) {
 
-                             // An inflow/dirichlet type of boundary condition
+                                        // An inflow/dirichlet type of boundary condition
         if (bme->second.first[di] == INFLOW_BC) {
           Wminus[q][di] = bvals[q](di);
         } else if (bme->second.first[di] == PRESSURE_BC) {
-                             // A prescribed pressure boundary condition.  This boundary
-                             // condition is complicated by the fact that even though
-                             // the pressure is prescribed, we really are setting
-                             // the energy index here, which will depend on velocity
-                             // and pressure. So even though this seems like a dirichlet
-                             // type boundary condition, we get sensitivities of
-                             // energy to velocity and density (unless these
-                             // are also prescribed.
+                                          // A prescribed pressure boundary condition.  This boundary
+                                          // condition is complicated by the fact that even though
+                                          // the pressure is prescribed, we really are setting
+                                          // the energy index here, which will depend on velocity
+                                          // and pressure. So even though this seems like a dirichlet
+                                          // type boundary condition, we get sensitivities of
+                                          // energy to velocity and density (unless these
+                                          // are also prescribed.
           fad_double rho_vel_sqr = 0;
           fad_double dens;
           
@@ -810,22 +805,22 @@ void ConsLaw<dim>::assemble_face_term(
               rho_vel_sqr += Wplus[q][d]*Wplus[q][d];
           }
           rho_vel_sqr /= dens;
-                             // Finally set the energy value as determined by the
-                             // prescribed pressure and the other variables.
+                                          // Finally set the energy value as determined by the
+                                          // prescribed pressure and the other variables.
           Wminus[q][di] = bvals[q](di)/(GAMMA-1.0) +
-                             0.5*rho_vel_sqr;
+                         0.5*rho_vel_sqr;
 
         } else if (bme->second.first[di] == OUTFLOW_BC) {
-                            // A free/outflow boundary, very simple.
+                                          // A free/outflow boundary, very simple.
           Wminus[q][di] = Wplus[q][di];
 
         } else { 
-                            // We must be at a no-penetration boundary.  We
-                            // prescribe the velocity (we are dealing with a
-                            // particular component here so that the average
-                            // of the velocities is orthogonal to the surface
-                            // normal.  This creates sensitivies of across
-                            // the velocity components.
+                                          // We must be at a no-penetration boundary.  We
+                                          // prescribe the velocity (we are dealing with a
+                                          // particular component here so that the average
+                                          // of the velocities is orthogonal to the surface
+                                          // normal.  This creates sensitivies of across
+                                          // the velocity components.
           fad_double vdotn = 0;
           for (unsigned int d = 0; d < dim; d++) {
             vdotn += Wplus[q][d]*normals[q](d);
@@ -834,60 +829,60 @@ void ConsLaw<dim>::assemble_face_term(
           Wminus[q][di] = Wplus[q][di] - 2.0*vdotn*normals[q](di);
         }
       }
-     } // for q
-   } // b>= 0
+    } // for q
+  } // b>= 0
    
-                           // Determine the Lax-Friedrich's stability parameter,
-                           // and evaluate the numerical flux function at the quadrature points
-   std::vector<std::vector<fad_double> > nflux(n_q_points, std::vector<fad_double>(get_n_components(), 0));
-     double alpha = 1;
-
-     switch(flux_params.LF_stab) {
-       case flux_params_type::CONSTANT:
-         alpha = flux_params.LF_stab_value;
-       break;
-       case flux_params_type::MESH:
-         alpha = face_diameter/(2.0*dT);
-       break;
-     }
-
-     LFNumFlux<fad_double, dim>(nflux, fe_v.get_quadrature_points(), normals, Wplus, Wminus,
-                     alpha);
-
-                          // Now assemble the face term
-     for (unsigned int i=0; i<fe_v.dofs_per_cell; ++i) {
-        if (!fe_v.get_fe().has_support_on_face(i, face_no)) continue;
-        F_i = 0;
-        for (unsigned int point=0; point<n_q_points; ++point)
-        {
-          const unsigned int
-            component_i = fe_v.get_fe().system_to_component_index(i).first;
-
-          F_i += nflux[point][component_i]*fe_v.shape_value_component(i, point, component_i)*JxW[point];
-
-         } 
-
-                          // Retrieve a pointer to the jacobian.
-        double *values = &(F_i.fastAccessDx(0));
-
-                          // Honestly, I forget why this can happen, but 
-                          // for some reason it can!!
-        if (!values) continue;
-
-                          // Update the matrix.  Depending on whether there
-                          // is/isn't a neighboring cell, we add more/less
-                          // entries.
-        Matrix->SumIntoGlobalValues(dofs[i],
-          dofs_per_cell, &values[0], reinterpret_cast<int*>(&dofs[0]));
-        if (boundary < 0) {
-          Matrix->SumIntoGlobalValues(dofs[i],
-            dofs_per_cell, &values[dofs_per_cell], reinterpret_cast<int*>(&dofs_neighbor[0]));
-        }
+                                  // Determine the Lax-Friedrich's stability parameter,
+                                  // and evaluate the numerical flux function at the quadrature points
+  std::vector<std::vector<fad_double> > nflux(n_q_points, std::vector<fad_double>(get_n_components(), 0));
+  double alpha = 1;
+
+  switch(flux_params.LF_stab) {
+    case flux_params_type::CONSTANT:
+         alpha = flux_params.LF_stab_value;
+         break;
+    case flux_params_type::MESH:
+         alpha = face_diameter/(2.0*dT);
+         break;
+  }
+
+  LFNumFlux<fad_double, dim>(nflux, fe_v.get_quadrature_points(), normals, Wplus, Wminus,
+                            alpha);
+
+                                  // Now assemble the face term
+  for (unsigned int i=0; i<fe_v.dofs_per_cell; ++i) {
+    if (!fe_v.get_fe().has_support_on_face(i, face_no)) continue;
+    F_i = 0;
+    for (unsigned int point=0; point<n_q_points; ++point)
+      {
+       const unsigned int
+         component_i = fe_v.get_fe().system_to_component_index(i).first;
+
+       F_i += nflux[point][component_i]*fe_v.shape_value_component(i, point, component_i)*JxW[point];
 
-                          // And add into the residual
-        right_hand_side(dofs[i]) -= F_i.val();
       } 
 
+                                    // Retrieve a pointer to the jacobian.
+    double *values = &(F_i.fastAccessDx(0));
+
+                                    // Honestly, I forget why this can happen, but 
+                                    // for some reason it can!!
+    if (!values) continue;
+
+                                    // Update the matrix.  Depending on whether there
+                                    // is/isn't a neighboring cell, we add more/less
+                                    // entries.
+    Matrix->SumIntoGlobalValues(dofs[i],
+                               dofs_per_cell, &values[0], reinterpret_cast<int*>(&dofs[0]));
+    if (boundary < 0) {
+      Matrix->SumIntoGlobalValues(dofs[i],
+                                 dofs_per_cell, &values[dofs_per_cell], reinterpret_cast<int*>(&dofs_neighbor[0]));
+    }
+
+                                    // And add into the residual
+    right_hand_side(dofs[i]) -= F_i.val();
+  } 
+
 }
                                  // @sect4{Assembling the whole system}
                                  // Now we put all of the assembly pieces together
@@ -900,8 +895,8 @@ void ConsLaw<dim>::assemble_system (double &res_norm)
   FESystem<dim> &fe = *fe_ptr;
   const unsigned int dofs_per_cell = dof_handler.get_fe().dofs_per_cell;
 
-                                  // We track the dofs on this cell and (if necessary)
-                                  // the adjacent cell.
+                                  // We track the dofs on this cell and (if necessary)
+                                  // the adjacent cell.
   std::vector<unsigned int> dofs (dofs_per_cell);
   std::vector<unsigned int> dofs_neighbor (dofs_per_cell);
 
@@ -1020,11 +1015,11 @@ void ConsLaw<dim>::assemble_system (double &res_norm)
                                                // fe_v and dofs as described
                                                // in the assembly routine.
              assemble_face_term(
-                                   face_no, fe_v_face,
-                                   fe_v_face,
-                                   dofs,
-                                   dofs,
-                                   face->boundary_indicator());
+               face_no, fe_v_face,
+               fe_v_face,
+               dofs,
+               dofs,
+               face->boundary_indicator());
            }
          else
            {
@@ -1038,7 +1033,7 @@ void ConsLaw<dim>::assemble_system (double &res_norm)
 
              if (face->has_children())
                {
-                  // case I: This cell refined compared to neighbor
+                                                  // case I: This cell refined compared to neighbor
 
                  const unsigned int neighbor2=
                    cell->neighbor_of_neighbor(face_no);
@@ -1064,13 +1059,13 @@ void ConsLaw<dim>::assemble_system (double &res_norm)
                      fe_v_face_neighbor.reinit (neighbor_child, neighbor2);
                      neighbor_child->get_dof_indices (dofs_neighbor);
 
-                                               // Assemble as if we are working with
-                                               // a DG element.
+                                                      // Assemble as if we are working with
+                                                      // a DG element.
                      assemble_face_term(
-                                             face_no, fe_v_subface,
-                                            fe_v_face_neighbor,
-                                             dofs,
-                                             dofs_neighbor);
+                       face_no, fe_v_subface,
+                       fe_v_face_neighbor,
+                       dofs,
+                       dofs_neighbor);
                      
                    }
                                                   // End of ``if
@@ -1078,13 +1073,13 @@ void ConsLaw<dim>::assemble_system (double &res_norm)
                }
              else
                {
-                                                  // We have no children, but 
-                                                  // the neighbor cell may be refine
-                                                  // compared to use
+                                                  // We have no children, but 
+                                                  // the neighbor cell may be refine
+                                                  // compared to use
                  neighbor->get_dof_indices (dofs_neighbor);
                  if (neighbor->level() != cell->level()) 
                    {
-                      // case II: This is refined compared to neighbor
+                                                      // case II: This is refined compared to neighbor
                      Assert(neighbor->level() < cell->level(), ExcInternalError());
                      const std::pair<unsigned int, unsigned int> faceno_subfaceno=
                        cell->neighbor_of_coarser_neighbor(face_no);
@@ -1107,31 +1102,31 @@ void ConsLaw<dim>::assemble_system (double &res_norm)
                                                    neighbor_subface_no);
                      
                      assemble_face_term(
-                                             face_no, fe_v_face,
-                                            fe_v_subface_neighbor,
-                                             dofs,
-                                             dofs_neighbor);
+                       face_no, fe_v_face,
+                       fe_v_subface_neighbor,
+                       dofs,
+                       dofs_neighbor);
 
                    }
 
                } 
-                                    // End of ``face not at boundary'':
+                                              // End of ``face not at boundary'':
            }
-                                    // End of loop over all faces:
+                                          // End of loop over all faces:
        } 
       
-                                     // End iteration through cells.
+                                      // End iteration through cells.
     } 
 
-                                     // Notify Epetra that the matrix is done.
-    Matrix->FillComplete();
+                                  // Notify Epetra that the matrix is done.
+  Matrix->FillComplete();
 
-                                    // Compute the nonlinear residual.
-    res_norm = right_hand_side.l2_norm();
+                                  // Compute the nonlinear residual.
+  res_norm = right_hand_side.l2_norm();
     
 }
 
-                                    // Create a conservation law with some defaults.
+                                // Create a conservation law with some defaults.
 template <int dim>
 ConsLaw<dim>::ConsLaw ()
                :
@@ -1149,14 +1144,14 @@ ConsLaw<dim>::ConsLaw ()
                 theta(0.5)
 {}
 
-                        // At one time this example could work for both DG and
-                        // continuous finite elements.  The choice was made here.
+                                // At one time this example could work for both DG and
+                                // continuous finite elements.  The choice was made here.
 template <int dim>
 void ConsLaw<dim>::build_fe() {
   fe_ptr = new FESystem<dim>(FE_Q<dim>(1), N_COMP);
 }
 
-                        // Bye bye Conservation law.
+                                // Bye bye Conservation law.
 template <int dim>
 ConsLaw<dim>::~ConsLaw () 
 {
@@ -1164,13 +1159,13 @@ ConsLaw<dim>::~ConsLaw ()
   delete fe_ptr;
 }
 
-                        // @sect3{Initialize System}
-                        // Sizes all of the vectors and sets up the
-                        // sparsity patter.  This function is called at
-                        // the very beginning of a simulation.  The function
-                        // <code> setup_system </code> repeats some of these
-                        // chores and is called after adaptivity in leiu
-                        // of this function.
+                                // @sect3{Initialize System}
+                                // Sizes all of the vectors and sets up the
+                                // sparsity patter.  This function is called at
+                                // the very beginning of a simulation.  The function
+                                // <code> setup_system </code> repeats some of these
+                                // chores and is called after adaptivity in leiu
+                                // of this function.
 template <int dim>
 void ConsLaw<dim>::initialize_system ()
 {
@@ -1189,9 +1184,9 @@ void ConsLaw<dim>::initialize_system ()
   indicator.reinit(triangulation.n_active_cells());
 }
 
-                                  // @sect3{Setup System}
-                                  // We call this function to build the sparsity
-                                  // and the matrix.
+                                // @sect3{Setup System}
+                                // We call this function to build the sparsity
+                                // and the matrix.
 template <int dim>
 void ConsLaw<dim>::setup_system ()
 {
@@ -1230,15 +1225,15 @@ void ConsLaw<dim>::setup_system ()
   for (unsigned int i=0; i<dof_handler.n_dofs(); ++i)
     row_lengths[i] = sparsity_pattern.row_length (i);
 
-                                  // Now we build the matrix, using
-                                  // the constructor that optimizes
-                                  // with the existing lengths per row
-                                  // variable.
+                                  // Now we build the matrix, using
+                                  // the constructor that optimizes
+                                  // with the existing lengths per row
+                                  // variable.
   if (Matrix) delete Matrix;
   Matrix = new Epetra_CrsMatrix(Copy, *Map, &row_lengths[0], true);
 
-                                 // We add the sparsity pattern to the matrix by
-                                 // inserting zeros.
+                                  // We add the sparsity pattern to the matrix by
+                                  // inserting zeros.
   const unsigned int max_nonzero_entries = *std::max_element (row_lengths.begin(),
                                                              row_lengths.end());
   std::vector<double> vals(max_nonzero_entries, 0);
@@ -1247,20 +1242,20 @@ void ConsLaw<dim>::setup_system ()
   unsigned int cur_row = 0;
   unsigned int cur_col = 0;
   for (SparsityPattern::iterator s_i = sparsity_pattern.begin(); 
-     s_i != sparsity_pattern.end(); s_i++) {
+       s_i != sparsity_pattern.end(); s_i++) {
     if (s_i->row() != cur_row) {
       Matrix->InsertGlobalValues(cur_row, cur_col, &vals[0], &row_indices[0]);
       cur_col = 0;
       cur_row = s_i->row();
     }
-  row_indices[cur_col++] = s_i->column();
+    row_indices[cur_col++] = s_i->column();
   }
-                                 // The last row.
+                                  // The last row.
   Matrix->InsertGlobalValues(cur_row, cur_col, &vals[0], &row_indices[0]);
 
-                                 // Epetra requires this function after building or
-                                 // filling a matrix.  It typically does some parallel
-                                 // bookeeping; perhaps more.
+                                  // Epetra requires this function after building or
+                                  // filling a matrix.  It typically does some parallel
+                                  // bookeeping; perhaps more.
   Matrix->FillComplete();
 
 }
@@ -1272,93 +1267,93 @@ template <int dim>
 void ConsLaw<dim>::solve (Vector<double> &dsolution, int &niter, double &lin_residual) 
 {
 
-                                 // We must hand the solvers Epetra vectors.
-                                 // Luckily, they support the concept of a 
-                                 // 'view', so we just send in a pointer to our
-                                 // dealii vectors.
-    Epetra_Vector x(View, *Map, dsolution.begin());
-    Epetra_Vector b(View, *Map, right_hand_side.begin());
+                                  // We must hand the solvers Epetra vectors.
+                                  // Luckily, they support the concept of a 
+                                  // 'view', so we just send in a pointer to our
+                                  // dealii vectors.
+  Epetra_Vector x(View, *Map, dsolution.begin());
+  Epetra_Vector b(View, *Map, right_hand_side.begin());
 
-                                 // The Direct option selects the Amesos solver.
+                                  // The Direct option selects the Amesos solver.
   if (solver_params.SOLVER == solver_params_type::DIRECT) {
    
-                                 // Setup for solving with
-                                 // Amesos. Other solvers are
-                                 // available and may be selected by
-                                 // changing th string given to the
-                                 // <code>Create</code> function.
-     Epetra_LinearProblem prob;
-     prob.SetOperator(Matrix);
-     Amesos_BaseSolver *solver = Amesos().Create ("Amesos_Klu", prob);
-
-     Assert (solver != NULL, ExcInternalError());
-
-                                 // There are two parts to the direct solve.
-                                 // As I understand, the symbolic part figures
-                                 // out the sparsity patterns, and then the
-                                 // numerical part actually performs Gaussian
-                                 // elimination or whatever the approach is.
-     if (solver_params.OUTPUT == solver_params_type::VERBOSE)
-       std::cout << "Starting Symbolic fact\n" << std::flush;
+                                    // Setup for solving with
+                                    // Amesos. Other solvers are
+                                    // available and may be selected by
+                                    // changing th string given to the
+                                    // <code>Create</code> function.
+    Epetra_LinearProblem prob;
+    prob.SetOperator(Matrix);
+    Amesos_BaseSolver *solver = Amesos().Create ("Amesos_Klu", prob);
+
+    Assert (solver != NULL, ExcInternalError());
+
+                                    // There are two parts to the direct solve.
+                                    // As I understand, the symbolic part figures
+                                    // out the sparsity patterns, and then the
+                                    // numerical part actually performs Gaussian
+                                    // elimination or whatever the approach is.
+    if (solver_params.OUTPUT == solver_params_type::VERBOSE)
+      std::cout << "Starting Symbolic fact\n" << std::flush;
 
-     solver->SymbolicFactorization();
+    solver->SymbolicFactorization();
 
-     if (solver_params.OUTPUT == solver_params_type::VERBOSE)
-         std::cout << "Starting Numeric fact\n" << std::flush;
+    if (solver_params.OUTPUT == solver_params_type::VERBOSE)
+      std::cout << "Starting Numeric fact\n" << std::flush;
 
-     solver->NumericFactorization();
+    solver->NumericFactorization();
 
     
-                                 // Define the linear problem by setting the
-                                 // right hand and left hand sides.
-     prob.SetRHS(&b);
-     prob.SetLHS(&x);
-                                 // And finally solve the problem.
-     if (solver_params.OUTPUT == solver_params_type::VERBOSE)
-       std::cout << "Starting solve\n" << std::flush;
-     solver->Solve();
-     niter = 0;
-     lin_residual = 0;
-
-                                 // We must free the solver that was created
-                                 // for us.
-     delete solver;
+                                    // Define the linear problem by setting the
+                                    // right hand and left hand sides.
+    prob.SetRHS(&b);
+    prob.SetLHS(&x);
+                                    // And finally solve the problem.
+    if (solver_params.OUTPUT == solver_params_type::VERBOSE)
+      std::cout << "Starting solve\n" << std::flush;
+    solver->Solve();
+    niter = 0;
+    lin_residual = 0;
+
+                                    // We must free the solver that was created
+                                    // for us.
+    delete solver;
 
   } else if (solver_params.SOLVER == solver_params_type::GMRES) {
 
-                                 // For the iterative solvers, we use Aztec.
+                                    // For the iterative solvers, we use Aztec.
     AztecOO Solver;
 
-                                 // Select the appropriate level of verbosity.
+                                    // Select the appropriate level of verbosity.
     if (solver_params.OUTPUT == solver_params_type::QUIET)
       Solver.SetAztecOption(AZ_output, AZ_none);
 
     if (solver_params.OUTPUT == solver_params_type::VERBOSE)
       Solver.SetAztecOption(AZ_output, AZ_all);
 
-                                 // Select gmres.  Other solvers are available.
+                                    // Select gmres.  Other solvers are available.
     Solver.SetAztecOption(AZ_solver, AZ_gmres);
     Solver.SetRHS(&b);
     Solver.SetLHS(&x);
 
-                                 // Set up the ILUT preconditioner.  I do not know
-                                 // why, but we must pretend like we are in parallel
-                                 // using domain decomposition or the preconditioner
-                                 // refuses to activate.
+                                    // Set up the ILUT preconditioner.  I do not know
+                                    // why, but we must pretend like we are in parallel
+                                    // using domain decomposition or the preconditioner
+                                    // refuses to activate.
     Solver.SetAztecOption(AZ_precond, AZ_dom_decomp);
     Solver.SetAztecOption(AZ_subdomain_solve, AZ_ilut);
     Solver.SetAztecOption(AZ_overlap, 0);
     Solver.SetAztecOption(AZ_reorder, 0);
 
-                                 // ILUT parameters as described above.
+                                    // ILUT parameters as described above.
     Solver.SetAztecParam(AZ_drop, solver_params.ILUT_DROP);
     Solver.SetAztecParam(AZ_ilut_fill, solver_params.ILUT_FILL);
     Solver.SetAztecParam(AZ_athresh, solver_params.ILUT_ATOL);
     Solver.SetAztecParam(AZ_rthresh, solver_params.ILUT_RTOL);
     Solver.SetUserMatrix(Matrix);
 
-                                 // Run the solver iteration.  Collect the number
-                                 // of iterations and the residual.
+                                    // Run the solver iteration.  Collect the number
+                                    // of iterations and the residual.
     Solver.Iterate(solver_params.MAX_ITERS, solver_params.RES);
     niter = Solver.NumIters();
     lin_residual = Solver.TrueResidual();
@@ -1378,19 +1373,19 @@ void ConsLaw<dim>::postprocess() {
                             | update_q_points
                             | update_JxW_values;
   UpdateFlags update_flags1 = update_values
-                            | update_gradients
-                            | update_q_points
-                            | update_JxW_values;
+                             | update_gradients
+                             | update_q_points
+                             | update_JxW_values;
 
- QGauss<dim>  quadrature_formula(4);
 QGauss<dim>  quadrature_formula(4);
 
- const std::vector<Point<dim> > &us = fe_ptr->base_element(0).get_unit_support_points();
 const std::vector<Point<dim> > &us = fe_ptr->base_element(0).get_unit_support_points();
 
 
- Quadrature<dim>  unit_support(us);
 Quadrature<dim>  unit_support(us);
 
- int n_q_points = quadrature_formula.n_quadrature_points;
- int n_uq_points = unit_support.n_quadrature_points;
 int n_q_points = quadrature_formula.n_quadrature_points;
 int n_uq_points = unit_support.n_quadrature_points;
 
   FEValues<dim> fe_v (
     mapping, *fe_ptr, quadrature_formula, update_flags);
@@ -1401,15 +1396,15 @@ void ConsLaw<dim>::postprocess() {
   std::vector<Vector<double> > U(n_uq_points,
                                  Vector<double>(get_n_components()));
   std::vector<Vector<double> > UU(n_q_points,
-                                 Vector<double>(get_n_components()));
+                                 Vector<double>(get_n_components()));
   std::vector<std::vector<Tensor<1,dim> > > dU(n_uq_points,
-                                            std::vector<Tensor<1,dim> >(get_n_components()));
+                                              std::vector<Tensor<1,dim> >(get_n_components()));
   
   typename DoFHandler<dim>::active_cell_iterator
     cell = dof_handler.begin_active(),
     endc = dof_handler.end();
 
-                                    // Loop the cells
+                                  // Loop the cells
   for (unsigned int cell_no=0; cell!=endc; ++cell, ++cell_no) {
     cell->get_dof_indices (dofs);
     fe_v_unit.reinit(cell);
@@ -1429,11 +1424,11 @@ void ConsLaw<dim>::postprocess() {
         rho_normVsqr += solution(dofs[vidx])*solution(dofs[vidx]);
       }
       rho_normVsqr /= solution(dofs[didx]);
-                                 // Pressure
+                                      // Pressure
       ppsolution(dofs[eidx]) = (GAMMA-1.0)*(solution(dofs[eidx]) - 0.5*rho_normVsqr);
 
-                                 // Either output density or gradient squared of density,
-                                 // depending on what the user wants.
+                                      // Either output density or gradient squared of density,
+                                      // depending on what the user wants.
       if (!schlieren_plot) {
         ppsolution(dofs[didx]) = solution(dofs[didx]);
       } else {
@@ -1448,9 +1443,9 @@ void ConsLaw<dim>::postprocess() {
 
 }
 
-                            // Loop and assign a value for refinement.  We
-                            // simply use the density squared, which selects
-                            // shocks with some success.
+                                // Loop and assign a value for refinement.  We
+                                // simply use the density squared, which selects
+                                // shocks with some success.
 template <int dim>
 void ConsLaw<dim>::estimate() {
   
@@ -1471,7 +1466,7 @@ void ConsLaw<dim>::estimate() {
   std::vector<Vector<double> > U(n_q_points,
                                  Vector<double>(get_n_components()));
   std::vector<std::vector<Tensor<1,dim> > > dU(n_q_points,
-                                            std::vector<Tensor<1,dim> >(get_n_components()));
+                                              std::vector<Tensor<1,dim> >(get_n_components()));
   
   typename DoFHandler<dim>::active_cell_iterator
     cell = dof_handler.begin_active(),
@@ -1505,11 +1500,11 @@ void ConsLaw<dim>::refine_grid ()
     cell = dof_handler.begin_active(),
     endc = dof_handler.end();
 
-                                  // Loop cells.  If the indicator
-                                  // for the cell matches the refinement criterion,
-                                  // refine, else unrefine.  The unrefinement has
-                                  // a slight hysterisis to avoid 'flashing' from refined
-                                  // to unrefined.
+                                  // Loop cells.  If the indicator
+                                  // for the cell matches the refinement criterion,
+                                  // refine, else unrefine.  The unrefinement has
+                                  // a slight hysterisis to avoid 'flashing' from refined
+                                  // to unrefined.
   for (unsigned int cell_no=0; cell!=endc; ++cell, ++cell_no) {
     cell->clear_coarsen_flag();
     cell->clear_refine_flag();
@@ -1518,13 +1513,13 @@ void ConsLaw<dim>::refine_grid ()
       cell->set_refine_flag();
     } else {
       if (cell->level() > 0 &&
-         std::fabs(indicator(cell_no)) < 0.75*refinement_params.shock_val)
-           cell->set_coarsen_flag();
+         std::fabs(indicator(cell_no)) < 0.75*refinement_params.shock_val)
+       cell->set_coarsen_flag();
     }
   }
 
-                                  // The following code prolongs the solution
-                                  // to the new grid and carries out the refinement.
+                                  // The following code prolongs the solution
+                                  // to the new grid and carries out the refinement.
   std::vector<Vector<double> > interp_in;
   std::vector<Vector<double> > interp_out;
 
@@ -1540,24 +1535,24 @@ void ConsLaw<dim>::refine_grid ()
   dof_handler.distribute_dofs (*fe_ptr);
 
   {
-  Vector<double> new_solution(1);
-  Vector<double> new_predictor(1);
+    Vector<double> new_solution(1);
+    Vector<double> new_predictor(1);
 
-  interp_out.push_back(new_solution);
-  interp_out.push_back(new_predictor);
-  interp_out[0].reinit(dof_handler.n_dofs());
-  interp_out[1].reinit(dof_handler.n_dofs());
+    interp_out.push_back(new_solution);
+    interp_out.push_back(new_predictor);
+    interp_out[0].reinit(dof_handler.n_dofs());
+    interp_out[1].reinit(dof_handler.n_dofs());
   }
 
   soltrans.interpolate(interp_in, interp_out);
   
-                               // Let the vector delete a very small vector
+                                  // Let the vector delete a very small vector
   solution.reinit(1);
   predictor.reinit(1);
   solution.swap(interp_out[0]);
   predictor.swap(interp_out[1]);
 
-                               // resize these vectors for the new grid.
+                                  // resize these vectors for the new grid.
   nlsolution.reinit(dof_handler.n_dofs());
   ppsolution.reinit(dof_handler.n_dofs());
   nlsolution = solution;
@@ -1600,15 +1595,15 @@ void ConsLaw<dim>::output_results (const unsigned int cycle) const
   output.close();
 }
 
-                                   // @sect3{Parsing the Input Deck}
-                                   // Declare the parameters for the
-                                   // input deck.  We assume a certain
-                                   // maximum number of boundaries and process
-                                   // any boundary the user supplies up to
-                                   // that maximum number.  We
-                                   // leave a detailed explanation of these 
-                                   // parameters to our description of the input
-                                   // sample file.
+                                // @sect3{Parsing the Input Deck}
+                                // Declare the parameters for the
+                                // input deck.  We assume a certain
+                                // maximum number of boundaries and process
+                                // any boundary the user supplies up to
+                                // that maximum number.  We
+                                // leave a detailed explanation of these 
+                                // parameters to our description of the input
+                                // sample file.
 const UInt MAX_BD = 10;
 template <int dim>
 void ConsLaw<dim>::declare_parameters() {
@@ -1619,160 +1614,160 @@ void ConsLaw<dim>::declare_parameters() {
                     "intput file");
 
   prm.declare_entry("diffusion power", "2.0",
-                     Patterns::Double(),
-                     "power of mesh size for diffusion");
+                   Patterns::Double(),
+                   "power of mesh size for diffusion");
 
   prm.declare_entry("gravity", "0.0",
-                     Patterns::Double(),
-                     "gravity forcing");
+                   Patterns::Double(),
+                   "gravity forcing");
 
                                    // Time stepping block
   prm.enter_subsection("time stepping");
-    prm.declare_entry("time step", "0.1",
-                     Patterns::Double(),
-                     "simulation time step");
-    prm.declare_entry("final time", "10.0",
-                     Patterns::Double(),
-                     "simulation end time");
+  prm.declare_entry("time step", "0.1",
+                   Patterns::Double(),
+                   "simulation time step");
+  prm.declare_entry("final time", "10.0",
+                   Patterns::Double(),
+                   "simulation end time");
   prm.leave_subsection();
 
 
-                                  // Declare the boundary parameters
+                                  // Declare the boundary parameters
   for (unsigned int b = 0; b < MAX_BD; b++) {
     char bd[512];
     std::sprintf(bd, "boundary_%d", b);
     prm.enter_subsection(bd);
     prm.declare_entry("no penetration", "false",
-                       Patterns::Selection("true|false"),
-                       "<true|false>");
-                                  // declare a slot for each of the conservative
-                                  // variables.
+                     Patterns::Selection("true|false"),
+                     "<true|false>");
+                                    // declare a slot for each of the conservative
+                                    // variables.
     for (unsigned int di = 0; di < N_COMP; di++) {
       char var[512];
       std::sprintf(var, "w_%d", di);
       prm.declare_entry(var, "outflow",
-                     Patterns::Selection(
-        "inflow|outflow|pressure"),
-        "<inflow|outflow|pressure>");
+                       Patterns::Selection(
+                         "inflow|outflow|pressure"),
+                       "<inflow|outflow|pressure>");
       
-                                   // for dirichlet, a function in x,y,z
+                                      // for dirichlet, a function in x,y,z
       std::sprintf(var, "w_%d value", di);
       prm.declare_entry(var, "0.0",
-                     Patterns::Anything(),
-                 "expression in x,y,z");
+                       Patterns::Anything(),
+                       "expression in x,y,z");
     }
 
     prm.leave_subsection();
   }
 
-                                // Initial condition block.
+                                  // Initial condition block.
   prm.enter_subsection("initial condition");
-    for (unsigned int di = 0; di < N_COMP; di++) {
-      char var[512];
-      std::sprintf(var, "w_%d", di);
+  for (unsigned int di = 0; di < N_COMP; di++) {
+    char var[512];
+    std::sprintf(var, "w_%d", di);
       
-      // for dirichlet, a function in x,y,z
-      std::sprintf(var, "w_%d value", di);
-      prm.declare_entry(var, "0.0",
-                     Patterns::Anything(),
-                 "expression in x,y,z");
-    }
+                                    // for dirichlet, a function in x,y,z
+    std::sprintf(var, "w_%d value", di);
+    prm.declare_entry(var, "0.0",
+                     Patterns::Anything(),
+                     "expression in x,y,z");
+  }
   prm.leave_subsection();
 
-                              // The linear solver block.
+                                  // The linear solver block.
   prm.enter_subsection("linear solver");
-    prm.declare_entry("output", "quiet",
-                     Patterns::Selection(
-                     "quiet|verbose"),
-                      "<quiet|verbose>");
-    prm.declare_entry("method", "gmres",
-                     Patterns::Selection(
-                     "gmres|direct"),
-                      "<gmres|direct>");
-    prm.declare_entry("residual", "1e-10",
-                     Patterns::Double(),
-                     "linear solver residual");
-    prm.declare_entry("max iters", "300",
-                     Patterns::Double(),
-                     "maximum solver iterations");
-    prm.declare_entry("ilut fill", "2",
-                     Patterns::Double(),
-                     "ilut preconditioner fill");
-    prm.declare_entry("ilut absolute tolerance", "1e-9",
-                     Patterns::Double(),
-                     "ilut preconditioner tolerance");
-    prm.declare_entry("ilut relative tolerance", "1.1",
-                     Patterns::Double(),
-                     "rel tol");
-    prm.declare_entry("ilut drop tolerance", "1e-10",
-                     Patterns::Double(),
-                     "ilut drop tol");
+  prm.declare_entry("output", "quiet",
+                   Patterns::Selection(
+                     "quiet|verbose"),
+                   "<quiet|verbose>");
+  prm.declare_entry("method", "gmres",
+                   Patterns::Selection(
+                     "gmres|direct"),
+                   "<gmres|direct>");
+  prm.declare_entry("residual", "1e-10",
+                   Patterns::Double(),
+                   "linear solver residual");
+  prm.declare_entry("max iters", "300",
+                   Patterns::Double(),
+                   "maximum solver iterations");
+  prm.declare_entry("ilut fill", "2",
+                   Patterns::Double(),
+                   "ilut preconditioner fill");
+  prm.declare_entry("ilut absolute tolerance", "1e-9",
+                   Patterns::Double(),
+                   "ilut preconditioner tolerance");
+  prm.declare_entry("ilut relative tolerance", "1.1",
+                   Patterns::Double(),
+                   "rel tol");
+  prm.declare_entry("ilut drop tolerance", "1e-10",
+                   Patterns::Double(),
+                   "ilut drop tol");
   prm.leave_subsection();
 
 
-                           // A refinement controller block.
+                                  // A refinement controller block.
   prm.enter_subsection("refinement");
-    prm.declare_entry("refinement", "none",
-                     Patterns::Selection(
-                     "none|fixed number|shock"),
-                      "<on|off>");
-    prm.declare_entry("refinement fraction", "0.1",
-                     Patterns::Double(),
-                     "Fraction of high refinement");
-    prm.declare_entry("unrefinement fraction", "0.1",
-                     Patterns::Double(),
-                     "Fraction of low unrefinement");
-    prm.declare_entry("max elements", "1000000",
-                     Patterns::Double(),
-                     "maximum number of elements");
-    prm.declare_entry("shock value", "4.0",
-                     Patterns::Double(),
-                     "value for shock indicator");
-    prm.declare_entry("shock levels", "3.0",
-                     Patterns::Double(),
-                     "number of shock refinement levels");
+  prm.declare_entry("refinement", "none",
+                   Patterns::Selection(
+                     "none|fixed number|shock"),
+                   "<on|off>");
+  prm.declare_entry("refinement fraction", "0.1",
+                   Patterns::Double(),
+                   "Fraction of high refinement");
+  prm.declare_entry("unrefinement fraction", "0.1",
+                   Patterns::Double(),
+                   "Fraction of low unrefinement");
+  prm.declare_entry("max elements", "1000000",
+                   Patterns::Double(),
+                   "maximum number of elements");
+  prm.declare_entry("shock value", "4.0",
+                   Patterns::Double(),
+                   "value for shock indicator");
+  prm.declare_entry("shock levels", "3.0",
+                   Patterns::Double(),
+                   "number of shock refinement levels");
   prm.leave_subsection();
 
-                      // Output control.
+                                  // Output control.
   prm.enter_subsection("output");
-    prm.declare_entry("density", "standard",
-                     Patterns::Selection(
-                     "standard|schlieren"),
-                      "<standard|schlieren>");
-    prm.declare_entry("step", "-1",
-                     Patterns::Double(),
-                     "output once per this period");
+  prm.declare_entry("density", "standard",
+                   Patterns::Selection(
+                     "standard|schlieren"),
+                   "<standard|schlieren>");
+  prm.declare_entry("step", "-1",
+                   Patterns::Double(),
+                   "output once per this period");
   prm.leave_subsection();
 
-                      // Flux control
+                                  // Flux control
   prm.enter_subsection("flux");
-    prm.declare_entry("stab", "alpha",
-                     Patterns::Selection(
-                     "alpha|constant|mesh"),
-                      "<alpha|constant|mesh>");
-    prm.declare_entry("stab value", "1",
-                     Patterns::Double(),
-                     "alpha stabilization");
+  prm.declare_entry("stab", "alpha",
+                   Patterns::Selection(
+                     "alpha|constant|mesh"),
+                   "<alpha|constant|mesh>");
+  prm.declare_entry("stab value", "1",
+                   Patterns::Double(),
+                   "alpha stabilization");
   prm.leave_subsection();
 
 
 }
 
-                     // Code to actually parse an input file.  This function
-                     // matches the declarations above.
+                                // Code to actually parse an input file.  This function
+                                // matches the declarations above.
 template <int dim>
 void ConsLaw<dim>::load_parameters(const char *infile){
 
   prm.read_input(infile);
 
-                     // The global parameters.
+                                  // The global parameters.
   mesh = prm.get("mesh");
 
   diffusion_power = prm.get_double("diffusion power");
 
   gravity = prm.get_double("gravity");
 
-                    // The time stepping.
+                                  // The time stepping.
   prm.enter_subsection("time stepping");
   dT = prm.get_double("time step");
   std::cout << "dT=" << dT << std::endl;
@@ -1786,12 +1781,12 @@ void ConsLaw<dim>::load_parameters(const char *infile){
   std::cout << "TF=" << TF << std::endl;
   prm.leave_subsection();
 
-                   // The boundary info
+                                  // The boundary info
   for (unsigned int b = 0; b < MAX_BD; b++) {
     std::vector<bc_type> flags(N_COMP, OUTFLOW_BC);
 
-                   // Define a parser for every boundary, though it may be
-                   // unused.
+                                    // Define a parser for every boundary, though it may be
+                                    // unused.
     SideCondition<dim> *sd = new SideCondition<dim>(N_COMP);
     char bd[512];
     std::sprintf(bd, "boundary_%d", b);
@@ -1799,7 +1794,7 @@ void ConsLaw<dim>::load_parameters(const char *infile){
 
     const std::string &nopen = prm.get("no penetration");
 
-                    // Determine how each component is handled.
+                                    // Determine how each component is handled.
     for (unsigned int di = 0; di < N_COMP; di++) {
       char var[512];
       std::sprintf(var, "w_%d", di);
@@ -1819,79 +1814,79 @@ void ConsLaw<dim>::load_parameters(const char *infile){
     } 
     prm.leave_subsection();
 
-                     // Add the boundary condition to the law.
+                                    // Add the boundary condition to the law.
     sd->Init();
     add_boundary(b, flags, sd);
-   }
+  }
 
-                     // Initial conditions.
-   prm.enter_subsection("initial condition");
-    for (unsigned int di = 0; di < N_COMP; di++) {
-      char var[512];
+                                  // Initial conditions.
+  prm.enter_subsection("initial condition");
+  for (unsigned int di = 0; di < N_COMP; di++) {
+    char var[512];
 
-      std::sprintf(var, "w_%d value", di);
-      std::string var_value = prm.get(var);
-      ic.set_ic(di, var_value);  
-    }
-    ic.Init();
-   prm.leave_subsection();
-
-                    // The linear solver.
- prm.enter_subsection("linear solver");
-    const std::string &op = prm.get("output");
-    if (op == "verbose") solver_params.OUTPUT = solver_params_type::VERBOSE;
-    if (op == "quiet") solver_params.OUTPUT = solver_params_type::QUIET;
-    const std::string &sv = prm.get("method");
-    if (sv == "direct") {
-      solver_params.SOLVER = solver_params_type::DIRECT;
-    } else if (sv == "gmres") {
-      solver_params.SOLVER = solver_params_type::GMRES;
-    } 
+    std::sprintf(var, "w_%d value", di);
+    std::string var_value = prm.get(var);
+    ic.set_ic(di, var_value);  
+  }
+  ic.Init();
+  prm.leave_subsection();
 
-    solver_params.RES = prm.get_double("residual");
-    solver_params.MAX_ITERS = (int) prm.get_double("max iters");
-    solver_params.ILUT_FILL = prm.get_double("ilut fill");
-    solver_params.ILUT_ATOL = prm.get_double("ilut absolute tolerance");
-    solver_params.ILUT_RTOL = prm.get_double("ilut relative tolerance");
-    solver_params.ILUT_DROP = prm.get_double("ilut drop tolerance");
-    solver_params.RES = prm.get_double("residual");
+                                  // The linear solver.
+  prm.enter_subsection("linear solver");
+  const std::string &op = prm.get("output");
+  if (op == "verbose") solver_params.OUTPUT = solver_params_type::VERBOSE;
+  if (op == "quiet") solver_params.OUTPUT = solver_params_type::QUIET;
+  const std::string &sv = prm.get("method");
+  if (sv == "direct") {
+    solver_params.SOLVER = solver_params_type::DIRECT;
+  } else if (sv == "gmres") {
+    solver_params.SOLVER = solver_params_type::GMRES;
+  } 
+
+  solver_params.RES = prm.get_double("residual");
+  solver_params.MAX_ITERS = (int) prm.get_double("max iters");
+  solver_params.ILUT_FILL = prm.get_double("ilut fill");
+  solver_params.ILUT_ATOL = prm.get_double("ilut absolute tolerance");
+  solver_params.ILUT_RTOL = prm.get_double("ilut relative tolerance");
+  solver_params.ILUT_DROP = prm.get_double("ilut drop tolerance");
+  solver_params.RES = prm.get_double("residual");
   prm.leave_subsection();
 
 
-                       // And refiement.
+                                  // And refiement.
   prm.enter_subsection("refinement");
-    const std::string &ref = prm.get("refinement");
-    if (ref == "none") {
-      refinement_params.refine = refinement_params_type::NONE;
-    } else if (ref == "fixed number") {
-      refinement_params.refine = refinement_params_type::FIXED_NUMBER;
-    } else if (ref == "shock") {
-      refinement_params.refine = refinement_params_type::SHOCK;
-    } else
+  const std::string &ref = prm.get("refinement");
+  if (ref == "none") {
+    refinement_params.refine = refinement_params_type::NONE;
+  } else if (ref == "fixed number") {
+    refinement_params.refine = refinement_params_type::FIXED_NUMBER;
+  } else if (ref == "shock") {
+    refinement_params.refine = refinement_params_type::SHOCK;
+  } else
     refinement_params.high_frac = prm.get_double("refinement fraction");
-    refinement_params.high_frac_sav = refinement_params.high_frac;
-    refinement_params.low_frac = prm.get_double("unrefinement fraction");
-    refinement_params.max_cells = prm.get_double("max elements");
-    refinement_params.shock_val = prm.get_double("shock value");
-    refinement_params.shock_levels = prm.get_double("shock levels");
+  refinement_params.high_frac_sav = refinement_params.high_frac;
+  refinement_params.low_frac = prm.get_double("unrefinement fraction");
+  refinement_params.max_cells = prm.get_double("max elements");
+  refinement_params.shock_val = prm.get_double("shock value");
+  refinement_params.shock_levels = prm.get_double("shock levels");
   prm.leave_subsection();
     
-                           // Output control.
+                                  // Output control.
   prm.enter_subsection("output");
-    const std::string &dens = prm.get("density");
-    schlieren_plot = dens == "schlieren" ? true : false;
-    output_step = prm.get_double("step");
+  const std::string &dens = prm.get("density");
+  schlieren_plot = dens == "schlieren" ? true : false;
+  output_step = prm.get_double("step");
   prm.leave_subsection();
 
-                           // Flux control.
+                                  // Flux control.
   prm.enter_subsection("flux");
-    const std::string &stab = prm.get("stab");
-    if (stab == "constant") {
-      flux_params.LF_stab = flux_params_type::CONSTANT;
-    } else if (stab == "mesh ") {
-      flux_params.LF_stab = flux_params_type::MESH;
-    }
-    flux_params.LF_stab_value = prm.get_double("stab value");
+  const std::string &stab = prm.get("stab");
+  if (stab == "constant") {
+    flux_params.LF_stab = flux_params_type::CONSTANT;
+  } else if (stab == "mesh ") {
+    flux_params.LF_stab = flux_params_type::MESH;
+  }
+  flux_params.LF_stab_value = prm.get_double("stab value");
   prm.leave_subsection();
 
 
@@ -1902,25 +1897,25 @@ void ConsLaw<dim>::zero_matrix() {
   Matrix->PutScalar(0); Matrix->FillComplete();
 }
 
-                          // We use a predictor to try and make adaptivity
-                          // work better.  The idea is to try and refine ahead
-                          // of a front, rather than stepping into a coarse
-                          // set of elements and smearing the solution.  This
-                          // simple time extrapolator does the job.
+                                // We use a predictor to try and make adaptivity
+                                // work better.  The idea is to try and refine ahead
+                                // of a front, rather than stepping into a coarse
+                                // set of elements and smearing the solution.  This
+                                // simple time extrapolator does the job.
 template<int dim>
 void ConsLaw<dim>::compute_predictor() {
   predictor = nlsolution;
   predictor.sadd(3/2.0, -1/2.0, solution);
 }
 
-                          // @sect3{Run the simulation}
-                         // Contains the initialization
-                          // the time loop, and the inner Newton iteration.
+                                // @sect3{Run the simulation}
+                                // Contains the initialization
+                                // the time loop, and the inner Newton iteration.
 template <int dim>
 void ConsLaw<dim>::run () 
 {
 
-                          // Open and load the mesh.
+                                  // Open and load the mesh.
   GridIn<dim> grid_in;
   grid_in.attach_triangulation(triangulation);
   std::cout << "Opening mesh <" << mesh << ">" << std::endl;
@@ -1935,30 +1930,30 @@ void ConsLaw<dim>::run ()
 
   unsigned int nstep = 0;
   
-                           // Initialize fields and matrices.
+                                  // Initialize fields and matrices.
   initialize_system (); 
   setup_system();
   initialize(); 
   predictor = solution;
 
-                          // Initial refinement.  We apply the ic,
-                          // estimate, refine, and repeat until
-                          // happy.
+                                  // Initial refinement.  We apply the ic,
+                                  // estimate, refine, and repeat until
+                                  // happy.
   if (refinement_params.refine != refinement_params_type::NONE)
-  for (unsigned int i = 0; i < refinement_params.shock_levels; i++) {
-    estimate();
-    refine_grid();
-    setup_system();
-    initialize(); 
-    predictor = solution;
-  }
+    for (unsigned int i = 0; i < refinement_params.shock_levels; i++) {
+      estimate();
+      refine_grid();
+      setup_system();
+      initialize(); 
+      predictor = solution;
+    }
   postprocess();
   output_results (nstep);
 
-                           // Determine when we will output next.
+                                  // Determine when we will output next.
   double next_output = T + output_step;
 
-                           // @sect4{Main time stepping loop}
+                                  // @sect4{Main time stepping loop}
   predictor = solution;
   while(T < TF)
     {
@@ -1978,8 +1973,8 @@ void ConsLaw<dim>::run ()
       double res_norm;
       int lin_iter;
 
-                              // Print some relevant information during the
-                              // Newton iteration.
+                                      // Print some relevant information during the
+                                      // Newton iteration.
       std::cout << "NonLin Res:       Lin Iter     Lin Res" << std::endl;
       std::cout << "______________________________________" << std::endl;
 
@@ -1987,35 +1982,35 @@ void ConsLaw<dim>::run ()
       unsigned int nonlin_iter = 0;
       double lin_res;
 
-                             // @sect5{Newton iteration}
+                                      // @sect5{Newton iteration}
       nlsolution = predictor;
       while (!nonlin_done) {
         lin_iter = 0;
         zero_matrix();
         right_hand_side = 0;
         assemble_system (res_norm);
-                            // Flash a star to the screen so one can
-                            // know when the assembly has stopped and the linear
-                            // solution is starting.
+                                        // Flash a star to the screen so one can
+                                        // know when the assembly has stopped and the linear
+                                        // solution is starting.
         std::cout << "* " << std::flush;
 
-                            // Test against a (hardcoded) nonlinear tolderance.
-                            // Do not solve the linear system at the last step 
-                            // (since it would be a waste).
+                                        // Test against a (hardcoded) nonlinear tolderance.
+                                        // Do not solve the linear system at the last step 
+                                        // (since it would be a waste).
                       
         if (fabs(res_norm) < 1e-10) {
           nonlin_done = true;
         } else {
-                            // Solve the linear system and update with the
-                            // delta.
-           dsolution = 0;
-           solve (dsolution, lin_iter, lin_res);
-           nlsolution.add(1.0, dsolution);
+                                          // Solve the linear system and update with the
+                                          // delta.
+         dsolution = 0;
+         solve (dsolution, lin_iter, lin_res);
+         nlsolution.add(1.0, dsolution);
         }
 
-                            // Print the residuals.
+                                        // Print the residuals.
         std::printf("%-16.3e %04d        %-5.2e\n",
-              res_norm, lin_iter, lin_res);
+                   res_norm, lin_iter, lin_res);
 
         ++nonlin_iter;
 
@@ -2023,7 +2018,7 @@ void ConsLaw<dim>::run ()
                     ExcMessage ("No convergence in nonlinear solver"));
       } 
 
-                           // Various post convergence tasks.
+                                      // Various post convergence tasks.
       compute_predictor();
 
       solution = nlsolution;
@@ -2035,7 +2030,7 @@ void ConsLaw<dim>::run ()
 
       T += dT;
 
-                          // Output if it is time.
+                                      // Output if it is time.
       if (output_step < 0) {
         output_results (++nstep);
       } else if (T >= next_output) {
@@ -2043,7 +2038,7 @@ void ConsLaw<dim>::run ()
         next_output += output_step;
       }
 
-                          // Refine, if refinement is selected.
+                                      // Refine, if refinement is selected.
       if (refinement_params.refine != refinement_params_type::NONE) {
         refine_grid();
         setup_system();

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.