From faa0e9adf04555d3d147216c3f652c6fd6d88069 Mon Sep 17 00:00:00 2001
From: Arne Morten Kvarving <arne.morten.kvarving@sintef.no>
Date: Thu, 1 Dec 2016 12:26:25 +0100
Subject: [PATCH] add CompatibleOperators::(Residual|Weak)::Convection

---
 Apps/Common/CompatibleOperators.C | 42 +++++++++++++++++++++++++++++++
 Apps/Common/CompatibleOperators.h | 27 ++++++++++++++++++--
 2 files changed, 67 insertions(+), 2 deletions(-)

diff --git a/Apps/Common/CompatibleOperators.C b/Apps/Common/CompatibleOperators.C
index df94934c..4a945b9d 100644
--- a/Apps/Common/CompatibleOperators.C
+++ b/Apps/Common/CompatibleOperators.C
@@ -28,6 +28,32 @@ void CompatibleOperators::Weak::Advection(std::vector<Matrix>& EM,
           EM[n](i,j) += scale*AC[k]*fe.grad(n)(j,k)*fe.basis(n)(i)*fe.detJxW;
 }
 
+
+void CompatibleOperators::Weak::Convection(std::vector<Matrix>& EM,
+                                           const FiniteElement& fe,
+                                           const Vec3& U,
+                                           const Tensor& dUdX,
+                                           double scale,
+                                           bool conservative)
+{
+  // Convection
+  static const double vidx[3][3] = {{1, 6, 7},
+                                    {10, 2, 11},
+                                    {14, 15, 3}};
+  size_t nsd = fe.grad(1).cols();
+  for (size_t m=1; m <= nsd; ++m)
+    for (size_t i=1; i <= fe.basis(m).size(); ++i)
+      for (size_t n = 1; n <= nsd; ++n)
+        for (size_t j=1; j <= fe.basis(n).size(); ++j) {
+          double conv = fe.basis(n)(j) * dUdX(m,n);
+          if (m == n)
+            for (size_t k = 1; k <= nsd; ++k)
+              conv += U[k-1] * fe.grad(n)(j,k);
+          EM[vidx[m-1][n-1]](i,j) += scale * conv * fe.basis(m)(i) * fe.detJxW;
+        }
+}
+
+
 void CompatibleOperators::Weak::Gradient(std::vector<Matrix>& EM,
                                         const FiniteElement& fe,
                                         double scale)
@@ -92,6 +118,22 @@ void CompatibleOperators::Weak::Source(Vectors& EV,
 }
 
 
+void CompatibleOperators::Residual::Convection(Vectors& EV, const FiniteElement& fe,
+                                               const Vec3& U, const Tensor& dUdX,
+                                               const Vec3& UC, double scale,
+                                               bool conservative)
+{
+  size_t nsd = fe.grad(1).cols();
+  for (size_t n=1; n <= nsd; ++n) {
+    double conv = 0.0;
+    for (size_t k = 1; k<= nsd; ++k)
+      conv += U[k-1]*dUdX(n, k);
+    for (size_t i=1; i <= fe.basis(n).size(); ++i)
+      EV[n](i) -= scale * conv * fe.basis(n)(i) * fe.detJxW;
+  }
+}
+
+
 void CompatibleOperators::Residual::Laplacian(Vectors& EV,
                                               const FiniteElement& fe,
                                               const Tensor& dUdX,
diff --git a/Apps/Common/CompatibleOperators.h b/Apps/Common/CompatibleOperators.h
index 6ca0926a..6ff3d423 100644
--- a/Apps/Common/CompatibleOperators.h
+++ b/Apps/Common/CompatibleOperators.h
@@ -28,7 +28,7 @@ class CompatibleOperators
 public:
   //! \brief Common weak operators using div-compatible discretizations.
   class Weak {
-    public:
+  public:
     //! \brief Compute an advection term.
     //! \param[out] EM The element matrices to add contribution to
     //! \param[in] fe The finite element to evaluate for
@@ -37,6 +37,16 @@ public:
     static void Advection(std::vector<Matrix>& EM, const FiniteElement& fe,
                           const Vec3& AC, double scale=1.0);
 
+    //! \brief Compute a (nonlinear) convection term.
+    //! \param[out] EM The element matrix to add contribution to
+    //! \param[in] fe The finite element to evaluate for
+    //! \param[in] U  Advecting field
+    //! \param[in] conservative True to use the conservative formulation
+    //! \param[in] basis Basis to use
+    static void Convection(std::vector<Matrix>& EM, const FiniteElement& fe,
+                           const Vec3& U, const Tensor& dUdX, double scale,
+                           bool conservative=false);
+
     //! \brief Compute a gradient term.
     //! \param[out] EM The element matrix to add contribution to
     //! \param[in] fe The finite element to evaluate for
@@ -86,7 +96,20 @@ public:
   //! \brief Common weak residual operators using div-compatible discretizations.
   class Residual
   {
-    public:
+  public:
+    //! \brief Compute a convection term in a residual vector.
+    //! \param EV The element vector to add contribution to
+    //! \param[in] fe The finite element to evaluate for
+    //! \param[in] U Advected field
+    //! \param[in] dUdX Advected field gradient
+    //! \param[in] UC Advecting field
+    //! \param[in] scale Scaling factor for contribution
+    //! \param[in] conservative True to use the conservative formulation
+    //! \param[in] basis Basis to use
+    static void Convection(Vectors& EV, const FiniteElement& fe,
+                           const Vec3& U, const Tensor& dUdX, const Vec3& UC,
+                           double scale, bool conservative=false);
+
     //! \brief Compute a laplacian term in a residual vector.
     //! \param[out] EV The element vector to add contribution to
     //! \param[in] fe The finite element to evaluate for