From: shakirbsm <shakirbsm@gmail.com>
Date: Wed, 26 Aug 2015 17:13:50 +0000 (+0200)
Subject: Advection residuals for boundary and face fixed
X-Git-Tag: v8.4.0-rc2~524^2
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=7c1e7ad36228cf17de64253e5c7f4b2ef045564f;p=dealii.git

Advection residuals for boundary and face fixed
---

diff --git a/include/deal.II/integrators/advection.h b/include/deal.II/integrators/advection.h
index 1e3f4f71bd..236cfe00df 100644
--- a/include/deal.II/integrators/advection.h
+++ b/include/deal.II/integrators/advection.h
@@ -353,15 +353,18 @@ namespace LocalIntegrators
 
     /**
      * Scalar case: Residual for upwind flux at the boundary for weak advection operator. This is the
-     * value of the trial function at the outflow boundary and zero else:
+     * value of the trial function at the outflow boundary and the value of the incoming boundary
+     * condition on the inflow boundary:
      * @f[
      * a_{ij} = \int_{\partial\Omega}
-     * [\mathbf w\cdot\mathbf n]_-
-     * (u_i-g) v_j \, ds
+     * (\mathbf w\cdot\mathbf n)
+     * \widehat u v_j \, ds
      * @f]
      *
-     * Here, <i>u</i> is the finite element function whose values are given in the argument <tt>input</tt>.
-     * <i>g</i> is the inhomogenous boundary value in the argument <tt>data</tt>.
+     * Here, the numerical flux $\widehat u$ is the upwind value at the face,
+     * namely the finite element function whose values are given in the argument
+     * `input` on the outflow boundary.
+     * On the inflow boundary, it is the inhomogenous boundary value in the argument `data`.
      *
      * The <tt>velocity</tt> is provided as a VectorSlice, having <tt>dim</tt>
      * vectors, one for each velocity component. Each of the vectors must
@@ -402,15 +405,13 @@ namespace LocalIntegrators
           for (unsigned int d=0; d<dim; ++d)
             nv += fe.normal_vector(k)[d] * velocity[d][k * v_increment];
 
-          if (nv > 0)
+          // Always use the upwind value
+          const double val = (nv > 0.) ? input[k] : -data[k];
+
+          for (unsigned i=0; i<n_dofs; ++i)
             {
-              for (unsigned i=0; i<n_dofs; ++i)
-                {
-                  const double v= fe.shape_value(i,k);
-                  const double u= input[k];
-                  const double g= data[k];
-                  result(i) += dx * nv *(u-g)*v;
-                }
+              const double v= fe.shape_value(i,k);
+              result(i) += dx * nv * val *v;
             }
         }
     }
@@ -419,13 +420,19 @@ namespace LocalIntegrators
 
     /**
      * Vector-valued case: Residual for upwind flux at the boundary for weak advection operator. This is the
-     * value of the trial function at the outflow boundary and zero else:
+     * value of the trial function at the outflow boundary and the value of the incoming boundary
+     * condition on the inflow boundary:
      * @f[
      * a_{ij} = \int_{\partial\Omega}
-     * [\mathbf w\cdot\mathbf n]_-
-     * (u_i-g) v_j \, ds
+     * (\mathbf w\cdot\mathbf n)
+     * \widehat u v_j \, ds
      * @f]
      *
+     * Here, the numerical flux $\widehat u$ is the upwind value at the face,
+     * namely the finite element function whose values are given in the argument
+     * `input` on the outflow boundary.
+     * On the inflow boundary, it is the inhomogenous boundary value in the argument `data`.
+     *
      * The <tt>velocity</tt> is provided as a VectorSlice, having <tt>dim</tt>
      * vectors, one for each velocity component. Each of the vectors must
      * either have only a single entry, if the advection velocity is
@@ -466,16 +473,16 @@ namespace LocalIntegrators
           for (unsigned int d=0; d<dim; ++d)
             nv += fe.normal_vector(k)[d] * velocity[d][k * v_increment];
 
-          if (nv > 0)
+          std::vector<double> val(n_comp);
+
+          for (unsigned int d=0; d<n_comp; ++d)
             {
+              val[d] = (nv > 0.) ? input[d][k] : -data[d][k];
               for (unsigned i=0; i<n_dofs; ++i)
-                for (unsigned int d=0; d<n_comp; ++d)
-                  {
-                    const double v= fe.shape_value_component(i,k,d);
-                    const double u= input[d][k];
-                    const double g= data[d][k];
-                    result(i) += dx * nv *(u-g)*v;
-                  }
+                {
+                  const double v= fe.shape_value_component(i,k,d);
+                  result(i) += dx * nv * val[d] *v;
+                }
             }
         }
     }