]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Remove gravity from input files. Move handling boundary conditions to a function...
authorbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 20 May 2008 20:49:10 +0000 (20:49 +0000)
committerbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 20 May 2008 20:49:10 +0000 (20:49 +0000)
git-svn-id: https://svn.dealii.org/trunk@16146 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/examples/step-33/doc/intro.dox
deal.II/examples/step-33/doc/results.dox
deal.II/examples/step-33/input.prm
deal.II/examples/step-33/step-33.cc

index ef696f4c252bb644c1c079fbd7e73f16b039b908..f4f6999dbb401ea12855b92246db37051def6a6e 100644 (file)
@@ -8,6 +8,11 @@ The code solves the basic Euler equations of gas dynamics, by using a
 fully implicit Newton iteration (inspired by Sandia's Aria code).  The
 code may be configured by an input deck to run different simulations
 on different meshes, with differing boundary conditions.
+<br>
+The original code and documentation was later slightly modified by Wolfgang
+Bangerth to make it more modular and allow replacing the parts that are
+specific to the Euler equations by other hyperbolic conservation laws without
+too much trouble.
 </i>
 
 <b>Note:</b>The program uses the <a
index b214cdedb78f00c246f08ff62cdc79d52c807d7e..458b755113c5a159e17443d6e301e9e346f3c5dc 100644 (file)
@@ -14,9 +14,6 @@ set mesh = slide.inp
 # Stabilization parameter
 set diffusion power = 2.0
 
-# Scaled value for gravity.  Positive means gravity points down.
-set gravity = 1.0
-
 # --------------------------------------------------
 # Boundary conditions
 # We may specify boundary conditions for up to MAX_BD boundaries.
index eb7711041bdd26d8f0f7e9baed27fda1ecfc8485..c0ae906b650cb0b1bce2f28293f053bc07726fc5 100644 (file)
@@ -7,9 +7,6 @@ set mesh = slide.inp
 # Stabilization parameter
 set diffusion power = 2.0
 
-# Scaled value for gravity.  Positive means gravity points down.
-set gravity = 1.0
-
 # --------------------------------------------------
 # Boundary conditions
 # We may specify boundary conditions for up to MAX_BD boundaries.
index 2409380ce1a3ebcc49d86a4484a127e59b9fa15d..3257bb4ed305d778cec9cfc0d7234e6aed89d8be 100644 (file)
@@ -120,6 +120,8 @@ using namespace dealii;
 template <int dim>
 struct EulerEquations
 {
+                                    // @sect4{Component description}
+    
                                     // First a few variables that
                                     // describe the various components of our
                                     // solution vector in a generic way. This
@@ -188,6 +190,8 @@ struct EulerEquations
       }
     
     
+                                    // @sect4{Transformations between variables}
+    
                                     // Next, we define the gas
                                     // constant. We will set it to 1.4
                                     // in its definition immediately
@@ -206,7 +210,7 @@ struct EulerEquations
                                     // and $O_2$.
     static const double gas_gamma;
 
-    
+
                                     // In the following, we will need to
                                     // compute the kinetic energy and the
                                     // pressure from a vector of conserved
@@ -266,6 +270,8 @@ struct EulerEquations
       }        
 
 
+                                    // @sect4{EulerEquations::compute_flux_matrix}    
+    
                                     // We define the flux function
                                     // $F(W)$ as one large matrix.
                                     // Each row of this matrix
@@ -335,6 +341,8 @@ struct EulerEquations
       }
 
 
+                                    // @sect4{EulerEquations::compute_normal_flux}
+
                                     // On the boundaries of the
                                     // domain and across hanging
                                     // nodes we use a numerical flux
@@ -345,12 +353,12 @@ struct EulerEquations
                                     // $\alpha$. It's form has also
                                     // been given already in the
                                     // introduction:
-    template <typename number>
+    template <typename InputVector>
     static
     void numerical_normal_flux(const Point<dim>          &normal,
-                              const std::vector<number> &Wplus,
-                              const std::vector<number> &Wminus,
-                              const double alpha,
+                              const InputVector         &Wplus,
+                              const InputVector         &Wminus,
+                              const double               alpha,
                               Sacado::Fad::DFad<double> (&normal_flux)[n_components])
       {
        Sacado::Fad::DFad<double> iflux[n_components][dim];
@@ -363,12 +371,13 @@ struct EulerEquations
          {
            normal_flux[di] = 0;
            for (unsigned int d=0; d<dim; ++d) 
-             normal_flux[di] += 0.5*(iflux[di][d] + oflux[di][d]) * normal(d);
+             normal_flux[di] += 0.5*(iflux[di][d] + oflux[di][d]) * normal[d];
              
            normal_flux[di] += 0.5*alpha*(Wplus[di] - Wminus[di]);
          }
       }
 
+                                    // @sect4{EulerEquations::compute_forcing_vector}
 
                                     // In the same way as describing the flux
                                     // function $\mathbf F(\mathbf w)$, we
@@ -377,10 +386,10 @@ struct EulerEquations
                                     // mentioned in the introduction, we
                                     // consider only gravity here, which
                                     // leads to the specific form $\mathbf
-                                    // G(\mathbf w) = \left( \begin{array}{c}
-                                    // g_1\rho \\ g_2\rho \\ g_3\rho \\ 0
-                                    // \\ \rho \mathbf g \cdot \mathbf v
-                                    // \end{array} \right)$, shown here for
+                                    // G(\mathbf w) = \left(
+                                    // g_1\rho, g_2\rho, g_3\rho, 0,
+                                    // \rho \mathbf g \cdot \mathbf v
+                                    // \right)^T$, shown here for
                                     // the 3d case. More specifically, we
                                     // will consider only $\mathbf
                                     // g=(0,0,-1)^T$ in 3d, or $\mathbf
@@ -408,8 +417,181 @@ struct EulerEquations
                    forcing[c] = 0;
            }
       }
+
+
+                                    // @sect4{Dealing with boundary conditions}
+
+                                    // Another thing we have to deal with is
+                                    // boundary conditions. To this end, let
+                                    // us first define the kinds of boundary
+                                    // conditions we currently know how to
+                                    // deal with:
+    enum BoundaryKind
+    {
+         inflow_boundary,
+         outflow_boundary,
+         no_penetration_boundary,
+         pressure_boundary
+    };
+
+
+                                    // The next part is to actually decide
+                                    // what to do at each kind of
+                                    // boundary. To this end, remember from
+                                    // the introduction that boundary
+                                    // conditions are specified by choosing a
+                                    // value $\mathbf w^-$ on the outside of
+                                    // a boundary given an inhomogeneity
+                                    // $\mathbf j$ and possibly the
+                                    // solution's value $\mathbf w^+$ on the
+                                    // inside. Both are then passed to the
+                                    // numerical flux $\mathbf
+                                    // H(\mathbf{w}^+, \mathbf{w}^-,
+                                    // \mathbf{n})$ to define boundary
+                                    // contributions to the bilinear form.
+                                    //
+                                    // Boundary conditions can in some cases
+                                    // be specified for each component of the
+                                    // solution vector independently. For
+                                    // example, if component $c$ is marked
+                                    // for inflow, then $w^-_c = j_c$. If it
+                                    // is an outflow, then $w^-_c =
+                                    // w^+_c$. These two simple cases are
+                                    // handled first in the function below.
+                                    //
+                                    // There is a little snag that makes this
+                                    // function unpleasant from a C++
+                                    // language viewpoint: The output vector
+                                    // <code>Wminus</code> will of course be
+                                    // modified, so it shouldn't be a
+                                    // <code>const</code> argument. Yet it is
+                                    // in the implementation below, and needs
+                                    // to be in order to allow the code to
+                                    // compile. The reason is that we call
+                                    // this function at a place where
+                                    // <code>Wminus</code> is of type
+                                    // <code>Table@<2,Sacado::Fad::DFad@<double@>
+                                    // @></code>, this being 2d table with
+                                    // indices representing the quadrature
+                                    // point and the vector component,
+                                    // respectively. We call this function
+                                    // with <code>Wminus[q]</code> as last
+                                    // argument; subscripting a 2d table
+                                    // yields a temporary accessor object
+                                    // representing a 1d vector, just what we
+                                    // want here. The problem is that a
+                                    // temporary accessor object can't be
+                                    // bound to a non-const reference
+                                    // argument of a function, as we would
+                                    // like here, according to the C++ 1998
+                                    // and 2003 standards (something that
+                                    // will be fixed with the next standard
+                                    // in the form of rvalue references).  We
+                                    // get away with making the output
+                                    // argument here a constant because it is
+                                    // the <i>accessor</i> object that's
+                                    // constant, not the table it points to:
+                                    // that one can still be written to. The
+                                    // hack is unpleasant nevertheless
+                                    // because it restricts the kind of data
+                                    // types that may be used as template
+                                    // argument to this function: a regular
+                                    // vector isn't going to do because that
+                                    // one can not be written to when marked
+                                    // <code>const</code>. With no good
+                                    // solution around at the moment, we'll
+                                    // go with the pragmatic, even if not
+                                    // pretty, solution shown here:
+    template <typename DataVector>
+    static
+    void
+    compute_Wminus (const BoundaryKind  (&boundary_kind)[n_components],
+                   const Point<dim>     &normal_vector,
+                   const DataVector     &Wplus,
+                   const Vector<double> &boundary_values,
+                   const DataVector     &Wminus)
+      {
+       for (unsigned int c = 0; c < n_components; c++)
+         switch (boundary_kind[c])
+           {
+             case inflow_boundary:
+             {
+               Wminus[c] = boundary_values(c);
+               break;
+             }
+
+             case outflow_boundary:
+             {
+               Wminus[c] = Wplus[c];
+               break;
+             }     
+               
+                                              // Prescribed pressure boundary
+                                              // conditions are a bit more
+                                              // complicated by the fact that
+                                              // even though the pressure is
+                                              // prescribed, we really are
+                                              // setting the energy component
+                                              // 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):
+             case pressure_boundary:
+             {
+               const typename DataVector::value_type
+                 density = (boundary_kind[density_component] ==
+                            inflow_boundary
+                            ?
+                            boundary_values(density_component)
+                            :
+                            Wplus[density_component]);
+
+               typename DataVector::value_type kinetic_energy = 0;
+               for (unsigned int d=0; d<dim; ++d)
+                 if (boundary_kind[d] == inflow_boundary)
+                   kinetic_energy += boundary_values(d)*boundary_values(d);
+                 else
+                   kinetic_energy += Wplus[d]*Wplus[d];
+               kinetic_energy *= 1./2./density;
+
+               Wminus[c] = boundary_values(c) / (gas_gamma-1.0) +
+                           kinetic_energy;
+                 
+               break;
+             }
+
+             case 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.
+               Sacado::Fad::DFad<double> vdotn = 0;
+               for (unsigned int d = 0; d < dim; d++) {
+                 vdotn += Wplus[d]*normal_vector[d];
+               }
+
+               Wminus[c] = Wplus[c] - 2.0*vdotn*normal_vector[c];
+               break;
+             }
+
+             default:
+                   Assert (false, ExcNotImplemented());
+           }    
+      }
+
+
     
     
+                                    // @sect4{EulerEquations::Postprocessor}
     
                                     // Finally, we declare a class that
                                     // implements a postprocessing of data
@@ -1149,17 +1331,11 @@ namespace Parameters
   {
       static const unsigned int max_n_boundaries = 10;
 
-      enum BoundaryKind
-      {
-           inflow_boundary,
-           outflow_boundary,
-           no_penetration_boundary,
-           pressure_boundary
-      };
-
       struct BoundaryConditions
       {
-         BoundaryKind        kind[EulerEquations<dim>::n_components];
+         typename EulerEquations<dim>::BoundaryKind
+         kind[EulerEquations<dim>::n_components];
+         
          FunctionParser<dim> values;
 
          BoundaryConditions ();
@@ -1169,7 +1345,6 @@ namespace Parameters
       AllParameters ();
       
       double diffusion_power;
-      double gravity;
 
       double time_step, final_time;
       double theta;
@@ -1212,10 +1387,6 @@ namespace Parameters
                      Patterns::Double(),
                      "power of mesh size for diffusion");
 
-    prm.declare_entry("gravity", "0.0",
-                     Patterns::Double(),
-                     "gravity forcing");
-
     prm.enter_subsection("time stepping");
     {
       prm.declare_entry("time step", "0.1",
@@ -1282,7 +1453,6 @@ namespace Parameters
   {
     mesh_filename = prm.get("mesh");
     diffusion_power = prm.get_double("diffusion power");
-    gravity = prm.get_double("gravity");
 
     prm.enter_subsection("time stepping");
     {
@@ -1318,13 +1488,17 @@ namespace Parameters
                = prm.get("w_" + Utilities::int_to_string(di));
 
              if ((di < dim) && (no_penetration == true))
-               boundary_conditions[boundary_id].kind[di] = no_penetration_boundary;
+               boundary_conditions[boundary_id].kind[di]
+                 = EulerEquations<dim>::no_penetration_boundary;
              else if (boundary_type == "inflow")
-               boundary_conditions[boundary_id].kind[di] = inflow_boundary;
+               boundary_conditions[boundary_id].kind[di]
+                 = EulerEquations<dim>::inflow_boundary;
              else if (boundary_type == "pressure")
-               boundary_conditions[boundary_id].kind[di] = pressure_boundary;
+               boundary_conditions[boundary_id].kind[di]
+                 = EulerEquations<dim>::pressure_boundary;
              else if (boundary_type == "outflow")
-               boundary_conditions[boundary_id].kind[di] = outflow_boundary;
+               boundary_conditions[boundary_id].kind[di]
+                 = EulerEquations<dim>::outflow_boundary;
              else
                AssertThrow (false, ExcNotImplemented());
 
@@ -1397,14 +1571,14 @@ class ConservationLaw
     void assemble_system ();
     void assemble_cell_term (const FEValues<dim>             &fe_v,
                             const std::vector<unsigned int> &dofs);
-    void assemble_face_term(const unsigned int           face_no,
-                           const FEFaceValuesBase<dim> &fe_v,
-                           const FEFaceValuesBase<dim> &fe_v_neighbor,
-                           const std::vector<unsigned int>   &dofs,
-                           const std::vector<unsigned int>   &dofs_neighbor,
-                           const bool                   external_face,
-                           const unsigned int           boundary_id,
-                           const double                 face_diameter);
+    void assemble_face_term (const unsigned int           face_no,
+                            const FEFaceValuesBase<dim> &fe_v,
+                            const FEFaceValuesBase<dim> &fe_v_neighbor,
+                            const std::vector<unsigned int>   &dofs,
+                            const std::vector<unsigned int>   &dofs_neighbor,
+                            const bool                   external_face,
+                            const unsigned int           boundary_id,
+                            const double                 face_diameter);
 
     std::pair<unsigned int, double> solve (Vector<double> &solution);
 
@@ -1836,13 +2010,13 @@ void ConservationLaw<dim>::assemble_system ()
 
                    neighbor_child->get_dof_indices (dof_indices_neighbor);
 
-                   assemble_face_term(face_no, fe_v_subface,
-                                      fe_v_face_neighbor,
-                                      dof_indices,
-                                      dof_indices_neighbor,
-                                      false,
-                                      numbers::invalid_unsigned_int,
-                                      neighbor_child->diameter());                   
+                   assemble_face_term (face_no, fe_v_subface,
+                                       fe_v_face_neighbor,
+                                       dof_indices,
+                                       dof_indices_neighbor,
+                                       false,
+                                       numbers::invalid_unsigned_int,
+                                       neighbor_child->diameter());                  
                  }
              }
 
@@ -1882,13 +2056,13 @@ void ConservationLaw<dim>::assemble_system ()
                                              neighbor_face_no,
                                              neighbor_subface_no);
                      
-               assemble_face_term(face_no, fe_v_face,
-                                  fe_v_subface_neighbor,
-                                  dof_indices,
-                                  dof_indices_neighbor,
-                                  false,
-                                  numbers::invalid_unsigned_int,
-                                  cell->face(face_no)->diameter());
+               assemble_face_term (face_no, fe_v_face,
+                                   fe_v_subface_neighbor,
+                                   dof_indices,
+                                   dof_indices_neighbor,
+                                   false,
+                                   numbers::invalid_unsigned_int,
+                                   cell->face(face_no)->diameter());
              }
          }
     } 
@@ -2264,12 +2438,15 @@ assemble_cell_term (const FEValues<dim>             &fe_v,
 
                                 // @sect4{ConservationLaw::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.
+                                // Here, we do essentially the same as in the
+                                // previous function. t the top, we introduce
+                                // the independent variables. Because the
+                                // current function is also used if we are
+                                // working on an internal face between two
+                                // cells, the independent variables are not
+                                // only the degrees of freedom on the current
+                                // cell but in the case of an interior face
+                                // also the ones on the neighbor.
 template <int dim>
 void
 ConservationLaw<dim>::assemble_face_term(const unsigned int           face_no,
@@ -2282,160 +2459,126 @@ ConservationLaw<dim>::assemble_face_term(const unsigned int           face_no,
                                         const double                 face_diameter) 
 {
   const unsigned int n_q_points = fe_v.n_quadrature_points;
-  const unsigned int dofs_per_cell = fe_v.get_fe().dofs_per_cell;
-  const unsigned int ndofs_per_cell = fe_v_neighbor.get_fe().dofs_per_cell;
-  Assert(dofs_per_cell == ndofs_per_cell,
-        ExcDimensionMismatch(dofs_per_cell, ndofs_per_cell));
-
-                                  // As above, the fad degrees of freedom
-  std::vector<Sacado::Fad::DFad<double> > independent_local_dof_values(dofs_per_cell+ndofs_per_cell);
-
-                                  // The conservative variables for this cell,
-                                  // and for 
-  std::vector<std::vector<Sacado::Fad::DFad<double> > > Wplus (n_q_points,
-                                                              std::vector<Sacado::Fad::DFad<double> >(EulerEquations<dim>::n_components));
-  std::vector<std::vector<Sacado::Fad::DFad<double> > > Wminus (n_q_points,
-                                                               std::vector<Sacado::Fad::DFad<double> >(EulerEquations<dim>::n_components));
-
-
-  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.
-  int ndofs = (external_face == false ? dofs_per_cell + ndofs_per_cell : dofs_per_cell);
-                                  // Set the local independent_local_dof_valuesS.
-  for (unsigned int in = 0; in < dofs_per_cell; in++) {
-    independent_local_dof_values[in] = current_solution(dof_indices[in]);
-    independent_local_dof_values[in].diff(in, ndofs);
-  }
-                                  // If present, set the neighbor dofs.
-  if (external_face == false)
-    for (unsigned int in = 0; in < ndofs_per_cell; in++) {
-      independent_local_dof_values[in+dofs_per_cell] = current_solution(dof_indices_neighbor[in]);
-      independent_local_dof_values[in+dofs_per_cell].diff(in+dofs_per_cell, ndofs);
-    }
+  const unsigned int dofs_per_cell = fe_v.dofs_per_cell;
 
-                                  // 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 < EulerEquations<dim>::n_components; di++) {
-      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] +=
-       (parameters.theta*independent_local_dof_values[sf]+(1.0-parameters.theta)*old_solution(dof_indices[sf]))*fe_v.shape_value_component(sf, q, di);
+  std::vector<Sacado::Fad::DFad<double> >
+    independent_local_dof_values (dofs_per_cell),
+    independent_neighbor_dof_values (external_face == false ?
+                                    dofs_per_cell :
+                                    0);
+
+  const unsigned int n_independent_variables = (external_face == false ?
+                                               2 * dofs_per_cell :
+                                               dofs_per_cell);
+  
+  for (unsigned int i = 0; i < dofs_per_cell; i++)
+    {
+      independent_local_dof_values[i] = current_solution(dof_indices[i]);
+      independent_local_dof_values[i].diff(i, n_independent_variables);
     }
 
+  if (external_face == false)
+    for (unsigned int i = 0; i < dofs_per_cell; i++)
+      {
+       independent_neighbor_dof_values[i]
+         = current_solution(dof_indices_neighbor[i]);
+       independent_neighbor_dof_values[i]
+         .diff(i+dofs_per_cell, n_independent_variables);
+      }
+
+
+                                  // Next, we need to define the values of
+                                  // the conservative variables $\tilde
+                                  // {\mathbf W}$ on this side of the face
+                                  // ($\tilde {\mathbf W}^+$) and on the
+                                  // opposite side ($\tilde {\mathbf
+                                  // W}^-$). The former can be computed in
+                                  // exactly the same way as in the previous
+                                  // function, but note that the
+                                  // <code>fe_v</code> variable now is of
+                                  // type FEFaceValues or FESubfaceValues:
+  Table<2,Sacado::Fad::DFad<double> >
+    Wplus (n_q_points, EulerEquations<dim>::n_components),
+    Wminus (n_q_points, EulerEquations<dim>::n_components);
 
-                                    // If there is a cell across, then initialize
-                                    // the exterior trace as a function of the other
-                                    // cell degrees of freedom.
-    if (external_face == false) {
-      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] +=
-         (parameters.theta*independent_local_dof_values[sf+dofs_per_cell]+(1.0-parameters.theta)*old_solution(dof_indices_neighbor[sf]))*
-         fe_v_neighbor.shape_value_component(sf, q, di);
+  for (unsigned int q=0; q<n_q_points; ++q)
+    for (unsigned int i=0; i<dofs_per_cell; ++i)
+      {
+       const unsigned int component_i = fe_v.get_fe().system_to_component_index(i).first;
+       Wplus[q][component_i] += (parameters.theta *
+                                 independent_local_dof_values[i]
+                                 +
+                                 (1.0-parameters.theta) *
+                                 old_solution(dof_indices[i])) *
+                                fe_v.shape_value_component(i, q, component_i);
       }
-    } 
-  } // for q
 
-                                  // If this is a boundary, then the values
+                                  // Computing $\tilde {\mathbf W}^-$ is a
+                                  // bit more complicated. If this is an
+                                  // internal face, we can compute it as
+                                  // above by simply using the independent
+                                  // variables from the neighbor:
+  if (external_face == false)
+    {
+      for (unsigned int q=0; q<n_q_points; ++q)
+       for (unsigned int i=0; i<dofs_per_cell; ++i)
+         {
+           const unsigned int component_i = fe_v_neighbor.get_fe().
+                                            system_to_component_index(i).first;
+           Wminus[q][component_i] += (parameters.theta *
+                                      independent_neighbor_dof_values[i]
+                                      +
+                                      (1.0-parameters.theta) *
+                                      old_solution(dof_indices_neighbor[i]))*
+                                     fe_v_neighbor.shape_value_component(i, q, component_i);
+         }
+    }
+                                  // On the other hand, if this is an
+                                  // external boundary face, 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,
+                                  // $W^+$, or they will be prescribed,
+                                  // depending on the kind of boundary
+                                  // condition imposed here. 
+                                  //
+                                  // To start the evaluation, let us ensure
+                                  // that the boundary id specified for this
+                                  // boundary is one for which we actually
+                                  // have data in the parameters
+                                  // object. Next, we evaluate the function
+                                  // object for the inhomogeneity.  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.
+                                  //
+                                  // The rest is done by a function that
+                                  // actually knows the specifics of Euler
+                                  // equation boundary conditions. Note that
+                                  // since we are using fad variables here,
                                   // sensitivities will be updated
-                                  // appropriately.  These sensitivities
-                                  // would be tremendously difficult to
-                                  // manage without fad!!!
-  if (external_face == true)
+                                  // appropriately, a process that would
+                                  // otherwise be tremendously complicated.
+  else
     {
       Assert (boundary_id < Parameters::AllParameters<dim>::max_n_boundaries,
              ExcIndexRange (boundary_id, 0,
                             Parameters::AllParameters<dim>::max_n_boundaries));
 
-                                      // 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>(EulerEquations<dim>::n_components));
-      parameters.boundary_conditions[boundary_id].values.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 < EulerEquations<dim>::n_components; di++) {
-
-                                          // An inflow/dirichlet type of boundary condition
-         if (parameters.boundary_conditions[boundary_id].kind[di] == Parameters::AllParameters<dim>::inflow_boundary) {
-           Wminus[q][di] = bvals[q](di);
-         } else if (parameters.boundary_conditions[boundary_id].kind[di] == Parameters::AllParameters<dim>::pressure_boundary) {
-                                            // 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.
-           Sacado::Fad::DFad<double> rho_vel_sqr = 0;
-           Sacado::Fad::DFad<double> dens;
-          
-           dens = parameters.boundary_conditions[boundary_id].kind[EulerEquations<dim>::density_component] == Parameters::AllParameters<dim>::inflow_boundary ? bvals[q](EulerEquations<dim>::density_component) :
-                  Wplus[q][EulerEquations<dim>::density_component];
-
-           for (unsigned int d=0; d < dim; d++) {
-             if (parameters.boundary_conditions[boundary_id].kind[d] == Parameters::AllParameters<dim>::inflow_boundary)
-               rho_vel_sqr += bvals[q](d)*bvals[q](d);
-             else
-               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.
-           Wminus[q][di] = bvals[q](di)/(EulerEquations<dim>::gas_gamma-1.0) +
-                           0.5*rho_vel_sqr;
-
-         } else if (parameters.boundary_conditions[boundary_id].kind[di] == Parameters::AllParameters<dim>::outflow_boundary) {
-                                            // 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.
-           Sacado::Fad::DFad<double> vdotn = 0;
-           for (unsigned int d = 0; d < dim; d++) {
-             vdotn += Wplus[q][d]*normals[q](d);
-           }
-
-           Wminus[q][di] = Wplus[q][di] - 2.0*vdotn*normals[q](di);
-         }
-       }
-      } // for q
-    } // b>= 0
-   
+      std::vector<Vector<double> >
+       boundary_values(n_q_points, Vector<double>(EulerEquations<dim>::n_components));
+      parameters.boundary_conditions[boundary_id]
+       .values.vector_value_list(fe_v.get_quadrature_points(),
+                                 boundary_values);
+
+      for (unsigned int q = 0; q < n_q_points; q++)
+       EulerEquations<dim>::compute_Wminus (parameters.boundary_conditions[boundary_id].kind,
+                                            fe_v.normal_vector(q),
+                                            Wplus[q],
+                                            boundary_values[q],
+                                            Wminus[q]);
+    }
+  
                                   // Determine the Lax-Friedrich's stability parameter,
                                   // and evaluate the numerical flux function at the quadrature points
   typedef Sacado::Fad::DFad<double> NormalFlux[EulerEquations<dim>::n_components];
@@ -2457,7 +2600,8 @@ ConservationLaw<dim>::assemble_face_term(const unsigned int           face_no,
     }
 
   for (unsigned int q=0; q<n_q_points; ++q)
-    EulerEquations<dim>::numerical_normal_flux(normals[q], Wplus[q], Wminus[q], alpha,
+    EulerEquations<dim>::numerical_normal_flux(fe_v.normal_vector(q),
+                                              Wplus[q], Wminus[q], alpha,
                                               normal_fluxes[q]);
 
                                   // Now assemble the face term

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.