changed: make WeakOperators a class with static functions

allows for templating over operator classes

Apps/Common/EqualOrderOperators.C

@@ -0,0 +1,238 @@
+//==============================================================================
+//!
+//! \file EqualOrderOperators.C
+//!
+//! \date Jul 22 2015
+//!
+//! \author Arne Morten Kvarving / SINTEF
+//!
+//! \brief Various discrete equal-ordered operators.
+//!
+//==============================================================================
+
+#include "EqualOrderOperators.h"
+#include "FiniteElement.h"
+#include "Vec3.h"
+
+//! \brief Helper for adding an element matrix to several components.
+//! \param[out] EM The element matrix to add to.
+//! \param[in] A The scalar element matrix to add.
+//! \param[in] cmp Number of components to add matrix to
+//! \param[in] nf Number of components in total matrix.
+//! \param[in] scmp Index of first component to add matrix to.
+static void addComponents(Matrix& EM, const Matrix& A,
+                          size_t cmp, size_t nf, size_t scmp)
+{
+  if (cmp == 1 && nf == 1)
+    EM += A;
+  else
+    for (size_t i = 1; i <= A.rows(); ++i)
+      for (size_t j = 1; j <= A.cols(); ++j)
+        for (size_t k = 1; k <= cmp; ++k)
+          EM(nf*(i-1)+k+scmp,nf*(j-1)+k+scmp) += A(i, j);
+}
+
+
+//! \brief Helper applying a divergence (1) or a gradient (2) operation
+template<int Operation>
+static void DivGrad(Matrix& EM, const FiniteElement& fe,
+            double scale, int basis, int tbasis)
+{
+  size_t nsd = fe.grad(basis).cols();
+  for (size_t i = 1; i <= fe.basis(tbasis).size();++i)
+    for (size_t j = 1; j <= fe.basis(basis).size();++j)
+      for (size_t k = 1; k <= nsd; ++k) {
+        double div = fe.basis(basis)(j)*fe.grad(tbasis)(i,k)*fe.detJxW;
+        if (Operation == 2)
+          EM((i-1)*nsd+k,j) += -scale*div;
+        if (Operation == 1)
+          EM(j, (i-1)*nsd+k) += scale*div;
+      }
+}
+
+
+void EqualOrderOperators::Weak::Advection(Matrix& EM, const FiniteElement& fe,
+                                          const Vec3& AC,
+                                          double scale, int basis)
+{
+  Matrix C(fe.basis(basis).size(), fe.basis(basis).size());
+  size_t ncmp = EM.rows() / C.rows();
+  for (size_t i = 1; i <= fe.basis(basis).size(); ++i) {
+    for (size_t j = 1; j <= fe.basis(basis).size(); ++j) {
+      // Sum convection for each direction
+      for (size_t k = 1; k <= fe.grad(basis).cols(); ++k)
+        C(i,j) += AC[k-1]*fe.grad(basis)(j,k);
+
+      C(i,j) *= scale*fe.basis(basis)(i)*fe.detJxW;
+    }
+  }
+  addComponents(EM, C, ncmp, ncmp, 0);
+}
+
+
+void EqualOrderOperators::Weak::Convection(Matrix& EM, const FiniteElement& fe,
+                                           const Vec3& U, const Tensor& dUdX,
+                                           double scale, bool conservative, int basis)
+{
+  size_t cmp = EM.rows() / fe.basis(basis).size();
+  if (conservative) {
+    Advection(EM, fe, U, -scale, basis);
+    for (size_t i = 1;i <= fe.basis(basis).size();i++)
+      for (size_t j = 1;j <= fe.basis(basis).size();j++) {
+        for (size_t k = 1;k <= cmp;k++) {
+          for (size_t l = 1;l <= cmp;l++)
+            EM((j-1)*cmp+l,(i-1)*cmp+k) -= scale*U[l-1]*fe.basis(basis)(i)*fe.grad(basis)(j,k)*fe.detJxW;
+        }
+      }
+  }
+  else {
+    Advection(EM, fe, U, scale, basis);
+    for (size_t i = 1;i <= fe.basis(basis).size();i++)
+      for (size_t j = 1;j <= fe.basis(basis).size();j++) {
+        for (size_t k = 1;k <= cmp;k++) {
+          for (size_t l = 1;l <= cmp;l++)
+            EM((j-1)*cmp+l,(i-1)*cmp+k) += scale*dUdX(l,k)*fe.basis(basis)(j)*fe.basis(basis)(i)*fe.detJxW;
+        }
+      }
+  }
+}
+
+
+void EqualOrderOperators::Weak::Divergence(Matrix& EM, const FiniteElement& fe,
+                                           double scale, int basis, int tbasis)
+{
+  DivGrad<1>(EM,fe,scale,basis,tbasis);
+}
+
+
+void EqualOrderOperators::Weak::Gradient(Matrix& EM, const FiniteElement& fe,
+                                         double scale, int basis, int tbasis)
+{
+  DivGrad<2>(EM,fe,scale,basis,tbasis);
+}
+
+
+void EqualOrderOperators::Weak::Divergence(Vector& EV, const FiniteElement& fe,
+                                           const Vec3& D, double scale, int basis)
+{
+  for (size_t i = 1; i <= fe.basis(basis).size(); ++i) {
+    double div=0.0;
+    for (size_t k = 1; k <= fe.grad(basis).cols(); ++k)
+      div += D[k-1]*fe.grad(basis)(i,k);
+    EV(i) += scale*div*fe.detJxW;
+  }
+}
+
+
+void EqualOrderOperators::Weak::Gradient(Vector& EV, const FiniteElement& fe,
+                                         double scale, int basis)
+{
+  size_t nsd = fe.grad(basis).cols();
+  for (size_t i = 1; i <= fe.basis(basis).size(); ++i)
+    for (size_t k = 1; k <= nsd; ++k)
+      EV((i-1)*nsd+k) += scale*fe.grad(basis)(i,k)*fe.detJxW;
+}
+
+
+void EqualOrderOperators::Weak::Laplacian(Matrix& EM, const FiniteElement& fe,
+                                          double scale, bool stress, int basis)
+{
+  size_t cmp = EM.rows() / fe.basis(basis).size();
+  Matrix A;
+  A.multiply(fe.grad(basis),fe.grad(basis),false,true);
+  A *= scale*fe.detJxW;
+  addComponents(EM, A, cmp, cmp, 0);
+  if (stress)
+    for (size_t i = 1; i <= fe.basis(basis).size(); i++)
+      for (size_t j = 1; j <= fe.basis(basis).size(); j++)
+        for (size_t k = 1; k <= cmp; k++)
+          for (size_t l = 1; l <= cmp; l++)
+            EM(cmp*(j-1)+k,cmp*(i-1)+l) += scale*fe.grad(basis)(i,k)*fe.grad(basis)(j,l)*fe.detJxW;
+}
+
+
+void EqualOrderOperators::Weak::LaplacianCoeff(Matrix& EM, const Matrix& K,
+                                               const FiniteElement& fe,
+                                               double scale, int basis)
+{
+  Matrix KB;
+  KB.multiply(K,fe.grad(basis),false,true).multiply(scale*fe.detJxW);
+  EM.multiply(fe.grad(basis),KB,false,false,true);
+}
+
+
+void EqualOrderOperators::Weak::Mass(Matrix& EM, const FiniteElement& fe,
+                                     double scale, int basis)
+{
+  size_t ncmp = EM.rows()/fe.basis(basis).size();
+  Matrix A;
+  A.outer_product(fe.basis(basis),fe.basis(basis),false);
+  A *= scale*fe.detJxW;
+  addComponents(EM, A, ncmp, ncmp, 0);
+}
+
+
+void EqualOrderOperators::Weak::Source(Vector& EV, const FiniteElement& fe,
+                                       double scale, int cmp, int basis)
+{
+  size_t ncmp = EV.size() / fe.basis(basis).size();
+  if (cmp == 1 && ncmp == 1)
+    EV.add(fe.basis(basis), scale*fe.detJxW);
+  else {
+    for (size_t i = 1; i <= fe.basis(basis).size(); ++i)
+      for (size_t k  = (cmp == 0 ? 1: cmp);
+                  k <= (cmp == 0 ? ncmp : cmp); ++k)
+        EV(ncmp*(i-1)+k) += scale*fe.basis(basis)(i)*fe.detJxW;
+  }
+}
+
+
+void EqualOrderOperators::Weak::Source(Vector& EV, const FiniteElement& fe,
+                                       const Vec3& f, double scale, int basis)
+{
+  size_t cmp = EV.size() / fe.basis(basis).size();
+  for (size_t i = 1; i <= fe.basis(basis).size(); ++i)
+    for (size_t k = 1; k <= cmp; ++k)
+      EV(cmp*(i-1)+k) += scale*f[k-1]*fe.basis(basis)(i)*fe.detJxW;
+}
+
+
+void EqualOrderOperators::Residual::Convection(Vector& EV, const FiniteElement& fe,
+                                                const Vec3& U, const Tensor& dUdX,
+                                               const Vec3& UC, double scale,
+                                               bool conservative, int basis)
+{
+  size_t cmp = EV.size() / fe.basis(basis).size();
+  if (conservative) {
+    for (size_t i = 1;i <= fe.basis(basis).size();i++)
+      for (size_t k = 1;k <= cmp;k++)
+        for (size_t l = 1;l <= cmp;l++)
+          // Convection
+          EV((i-1)*cmp + k) += scale*U[k-1]*UC[l-1]*fe.grad(basis)(i,l)*fe.detJxW;
+  }
+  else {
+    for (size_t k = 1;k <= cmp;k++) {
+      // Convection
+      double conv = 0.0;
+      for (size_t l = 1;l <= cmp;l++)
+        conv += UC[l-1]*dUdX(k,l);
+      conv *= scale*fe.detJxW;
+
+      for (size_t i = 1;i <= fe.basis(basis).size();i++)
+        EV((i-1)*cmp + k) -= conv*fe.basis(basis)(i);
+    }
+  }
+}
+
+
+void EqualOrderOperators::Residual::Divergence(Vector& EV, const FiniteElement& fe,
+                                                const Tensor& dUdX,
+                                                double scale, size_t basis)
+{
+  for (size_t i = 1; i <= fe.basis(basis).size(); ++i) {
+    double div=0.0;
+    for (size_t k = 1; k <= fe.grad(basis).cols(); ++k)
+      div += dUdX(k,k);
+    EV(i) += scale*div*fe.basis(basis)(i)*fe.detJxW;
+  }
+}

Apps/Common/EqualOrderOperators.h

@@ -0,0 +1,191 @@
+//==============================================================================
+//!
+//! \file EqualOrderOperators.h
+//!
+//! \date Jul 22 2015
+//!
+//! \author Arne Morten Kvarving / SINTEF
+//!
+//! \brief Various equal-ordered discrete operators.
+//!
+//==============================================================================
+
+#ifndef EQUAL_ORDER_OPERATORS_H
+#define EQUAL_ORDER_OPERATORS_H
+
+class Vec3;
+
+#include "FiniteElement.h"
+#include "MatVec.h"
+
+/*! \brief Common discrete operators using equal-ordered discretizations.
+ */
+
+class EqualOrderOperators
+{
+public:
+  //! \brief Common weak operators using equal-ordered discretizations.
+  class Weak {
+  public:
+    //! \brief Compute an advection term.
+    //! \param[out] EM The element matrix to add contribution to
+    //! \param[in] fe The finite element to evaluate for
+    //! \param[in] AC Advecting field
+    //! \param[in] scale Scaling factor for contribution
+    //! \param[in] basis Basis to use
+    static void Advection(Matrix& EM, const FiniteElement& fe,
+                          const Vec3& AC, double scale=1.0, int basis=1);
+
+    //! \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(Matrix& EM, const FiniteElement& fe,
+                           const Vec3& U, const Tensor& dUdX, double scale,
+                           bool conservative=false, int basis=1);
+
+    //! \brief Compute a divergence term.
+    //! \param[out] EM The element matrix to add contribution to
+    //! \param[in] fe The finite element to evaluate for
+    //! \param[in] scale Scaling factor for contribution
+    //! \param[in] basis Basis for field
+    //! \param[in] tbasis Test function basis
+    static void Divergence(Matrix& EM, const FiniteElement& fe,
+                           double scale=1.0, int basis=1, int tbasis=1);
+
+    //! \brief Compute a divergence term.
+    //! \param[out] EV The element vector to add contribution to
+    //! \param[in] fe The finite element to evaluate for
+    //! \param[in] D Divergence of field
+    //! \param[in] scale Scaling factor for contribution
+    //! \param[in] basis Test function basis
+    static void Divergence(Vector& EV, const FiniteElement& fe,
+                           const Vec3& D, double scale=1.0, int basis=1);
+
+    //! \brief Compute a gradient term for a (potentially mixed) vector/scalar field.
+    //! \param[out] EM The element matrix to add contribution to
+    //! \param[in] fe The finite element to evaluate for
+    //! \param[in] scale Scaling factor for contribution
+    //! \param[in] basis Basis for field
+    //! \param[in] tbasis Test function basis
+    static void Gradient(Matrix& EM, const FiniteElement& fe,
+                         double scale=1.0, int basis=1, int tbasis=1);
+
+    //! \brief Compute a gradient term.
+    //! \param[out] EV The element vector to add contribution to
+    //! \param[in] fe The finite element to evaluate for
+    //! \param[in] scale Scaling factor for contribution
+    //! \param[in] tbasis Test function basis
+    static void Gradient(Vector& EV, const FiniteElement& fe,
+                         double scale=1.0, int basis=1);
+
+    //! \brief Compute a laplacian.
+    //! \param[out] EM The element matrix to add contribution to
+    //! \param[in] fe The finite element to evaluate for
+    //! \param[in] scale Scaling factor for contribution
+    //! \param[in] stress Whether to add extra stress formulation terms
+    //! \param[in] basis Basis to use
+    static void Laplacian(Matrix& EM, const FiniteElement& fe,
+                          double scale=1.0, bool stress=false, int basis=1);
+
+    //! \brief Compute a heteregenous coefficient laplacian.
+    //! \param[out] EM The element matrix to add contribution to
+    //! \param[out] K The coefficient matrix
+    //! \param[in] fe The finite element to evaluate for
+    //! \param[in] scale Scaling factor for contribution
+    static void LaplacianCoeff(Matrix& EM, const Matrix& K, const FiniteElement& fe,
+                               double scale=1.0, int basis=1);
+
+    //! \brief Compute a mass term.
+    //! \param[out] EM The element matrix to add contribution to
+    //! \param[in] fe The finite element to evaluate for
+    //! \param[in] scale Scaling factor for contribution
+    //! \param[in] basis Basis to use
+    static void Mass(Matrix& EM, const FiniteElement& fe,
+                     double scale=1.0, int basis=1);
+
+    //! \brief Compute a source term.
+    //! \param[out] EV The element vector to add contribution to
+    //! \param[in] fe The finite element to evaluate for
+    //! \param[in] scale Scaling factor for contribution
+    //! \param[in] basis Basis to use
+    //! \param[in] cmp Component to add (0 for all)
+    static void Source(Vector& EV, const FiniteElement& fe,
+                       double scale=1.0, int cmp=1, int basis=1);
+
+    //! \brief Compute a vector-source term.
+    //! \param[out] EV The element vector to add contribution to
+    //! \param[in] fe The finite element to evaluate for
+    //! \param[in] scale Vector with contributions
+    //! \param[in] scale Scaling factor for contribution
+    //! \param[in] basis Basis to use
+    static void Source(Vector& EV, const FiniteElement& fe,
+                       const Vec3& f, double scale=1.0, int basis=1);
+  };
+
+  //! \brief Common weak residual operators using equal-ordered discretizations.
+  class Residual {
+  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(Vector& EV, const FiniteElement& fe,
+                           const Vec3& U, const Tensor& dUdX, const Vec3& UC,
+                           double scale, bool conservative=false, int basis=1);
+
+    //! \brief Compute a divergence 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] dUdX Gradient of field
+    //! \param[in] scale Scaling factor for contribution
+    //! \param[in] basis Basis to use
+    static void Divergence(Vector& EV, const FiniteElement& fe,
+                           const Tensor& dUdX, double scale=1.0,
+                           size_t basis=1);
+
+    //! \brief Compute a laplacian 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] dUdX Current solution gradient
+    //! \param[in] scale Scaling factor for contribution
+    //! \param[in] stress Whether to add extra stress formulation terms
+    //! \param[in] basis Basis to use
+    template<class T>
+    static void Laplacian(Vector& EV, const FiniteElement& fe,
+                          const T& dUdX, double scale=1.0,
+                          bool stress=false, int basis=1)
+    {
+      size_t cmp = EV.size() / fe.basis(basis).size();
+      for (size_t i = 1;i <= fe.basis(basis).size();i++) {
+        for (size_t k = 1;k <= cmp;k++) {
+          double diff = 0.0;
+          for (size_t l = 1;l <= cmp;l++)
+            diff += dUdX(k,l)*fe.grad(basis)(i,l);
+          diff *= scale*fe.detJxW;
+
+          // Add negative residual to rhs of momentum equation
+          EV((i-1)*cmp + k) -= diff;
+        }
+      }
+
+      // Use stress formulation
+      if (stress) {
+        for (size_t i = 1;i <= fe.basis(basis).size();i++)
+          for (size_t k = 1;k <= cmp;k++)
+            for (size_t l = 1;l <= cmp;l++)
+              // Diffusion
+              EV((i-1)*cmp + k) -= scale*dUdX(l,k)*fe.grad(basis)(i,l)*fe.detJxW;
+      }
+    }
+  };
+};
+
+#endif

Apps/Common/ResidualOperators.h

@@ -1,114 +0,0 @@
-//==============================================================================
-//!
-//! \file ResidualOperators.h
-//!
-//! \date Aug 19 2015
-//!
-//! \author Arne Morten Kvarving / SINTEF
-//!
-//! \brief Various weak, discrete residual operators
-//!
-//==============================================================================
-
-#ifndef RESIDUALOPERATORS_H_
-#define RESIDUALOPERATORS_H_
-
-class Vec3;
-
-#include "FiniteElement.h"
-#include "MatVec.h"
-
-namespace ResidualOperators
-{
-  //! \brief Compute a convection term in a residual vector.
-  //! \param[out] 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
-  template<class T>
-  void Convection(Vector& EV, const FiniteElement& fe,
-                  const Vec3& U, const T& dUdX, const Vec3& UC,
-                  double scale, bool conservative=false, int basis=1)
-  {
-    size_t cmp = EV.size() / fe.basis(basis).size();
-    if (conservative) {
-      for (size_t i = 1;i <= fe.basis(basis).size();i++)
-        for (size_t k = 1;k <= cmp;k++)
-          for (size_t l = 1;l <= cmp;l++)
-            // Convection
-            EV((i-1)*cmp + k) += scale*U[k-1]*UC[l-1]*fe.grad(basis)(i,l)*fe.detJxW;
-    }
-    else {
-      for (size_t k = 1;k <= cmp;k++) {
-        // Convection
-        double conv = 0.0;
-        for (size_t l = 1;l <= cmp;l++)
-          conv += UC[l-1]*dUdX(k,l);
-        conv *= scale*fe.detJxW;
-
-        for (size_t i = 1;i <= fe.basis(basis).size();i++)
-          EV((i-1)*cmp + k) -= conv*fe.basis(basis)(i);
-      }
-    }
-  }
-
-  //! \brief Compute a divergence term in a residual vector.
-  //! \param[out] EV The element vector to add contribution to
-  //! \param[in] fe The finite element to evaluate for
-  //! \param[in] dUdX Gradient of field
-  //! \param[in] scale Scaling factor for contribution
-  //! \param[in] basis Basis to use
-  template<class T>
-  void Divergence(Vector& EV, const FiniteElement& fe,
-                  const T& dUdX, double scale=1.0,
-                  size_t basis=1)
-  {
-    for (size_t i = 1; i <= fe.basis(basis).size(); ++i) {
-      double div=0.0;
-      for (size_t k = 1; k <= fe.grad(basis).cols(); ++k)
-        div += dUdX(k,k);
-      EV(i) += scale*div*fe.basis(basis)(i)*fe.detJxW;
-    }
-  }
-
-  //! \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
-  //! \param[in] dUdX Current solution gradient
-  //! \param[in] scale Scaling factor for contribution
-  //! \param[in] stress Whether to add extra stress formulation terms
-  //! \param[in] basis Basis to use
-  template<class T>
-  void Laplacian(Vector& EV, const FiniteElement& fe,
-                 const T& dUdX, double scale=1.0,
-                 bool stress=false, int basis=1)
-  {
-    size_t cmp = EV.size() / fe.basis(basis).size();
-    for (size_t i = 1;i <= fe.basis(basis).size();i++) {
-      for (size_t k = 1;k <= cmp;k++) {
-        double diff = 0.0;
-        for (size_t l = 1;l <= cmp;l++)
-          diff += dUdX(k,l)*fe.grad(basis)(i,l);
-        diff *= scale*fe.detJxW;
-
-        // Add negative residual to rhs of momentum equation
-        EV((i-1)*cmp + k) -= diff;
-      }
-    }
-
-    // Use stress formulation
-    if (stress) {
-      for (size_t i = 1;i <= fe.basis(basis).size();i++)
-	for (size_t k = 1;k <= cmp;k++)
-	  for (size_t l = 1;l <= cmp;l++)
-	    // Diffusion
-	    EV((i-1)*cmp + k) -= scale*dUdX(l,k)*fe.grad(basis)(i,l)*fe.detJxW;
-    }
-  }
-}
-
-#endif

Apps/Common/Test/TestEqualOrderOperators.C

@@ -1,16 +1,16 @@
 //==============================================================================
 //!
-//! \file TestWeakoperators.C
+//! \file TestEqualOrderOperators.C
 //!
 //! \date Feb 16 2015
 //!
 //! \author Arne Morten Kvarving / SINTEF
 //!
-//! \brief Tests for various mesh quality indicators
+//! \brief Tests for various discrete equal-order operators.
 //!
 //==============================================================================
 
-#include "WeakOperators.h"
+#include "EqualOrderOperators.h"
 #include "gtest/gtest.h"
 #include "FiniteElement.h"
 
@@ -36,7 +36,7 @@ static FiniteElement getFE()
   return fe;
 }
 
-TEST(TestWeakOperators, Advection)
+TEST(TestEqualOrderOperators, Advection)
 {
   FiniteElement fe = getFE();
 
@@ -44,14 +44,14 @@ TEST(TestWeakOperators, Advection)
   Vec3 U;
   U[0] = 1.0; U[1] = 2.0;
   // First component only
-  WeakOperators::Advection(EM_scalar, fe, U);
+  EqualOrderOperators::Weak::Advection(EM_scalar, fe, U);
   const DoubleVec EM_scalar_ref = {{ 7.0, 10},
                                    {14.0, 20.0}};
   check_matrix_equal(EM_scalar, EM_scalar_ref);
 
   Matrix EM_vec(2*2, 2*2);
   U[0] = 3.0; U[1] = 4.0;
-  WeakOperators::Advection(EM_vec, fe, U, 2.0);
+  EqualOrderOperators::Weak::Advection(EM_vec, fe, U, 2.0);
   const DoubleVec EM_vec_ref = {{30.0,  0.0, 44.0,  0.0},
                                 { 0.0, 30.0,  0.0, 44.0},
                                 {60.0,  0.0, 88.0,  0.0},
@@ -60,12 +60,12 @@ TEST(TestWeakOperators, Advection)
 }
 
 
-TEST(TestWeakOperators, Divergence)
+TEST(TestEqualOrderOperators, Divergence)
 {
   FiniteElement fe = getFE();
 
   Matrix EM_scalar(2, 2*2);
-  WeakOperators::Divergence(EM_scalar, fe);
+  EqualOrderOperators::Weak::Divergence(EM_scalar, fe);
   const DoubleVec EM_scalar_ref = {{1.0, 3.0, 2.0, 4.0},
                                    {2.0, 6.0, 4.0, 8.0}};
   check_matrix_equal(EM_scalar, EM_scalar_ref);
@@ -73,7 +73,7 @@ TEST(TestWeakOperators, Divergence)
   Vector EV_scalar(4);
   Vec3 D;
   D[0] = 1.0; D[1] = 2.0;
-  WeakOperators::Divergence(EV_scalar, fe, D, 1.0);
+  EqualOrderOperators::Weak::Divergence(EV_scalar, fe, D, 1.0);
   ASSERT_NEAR(EV_scalar(1),  7.0, 1e-13);
   ASSERT_NEAR(EV_scalar(2), 10.0, 1e-13);
   ASSERT_NEAR(EV_scalar(3),  0.0, 1e-13);
@@ -81,13 +81,13 @@ TEST(TestWeakOperators, Divergence)
 }
 
 
-TEST(TestWeakOperators, Gradient)
+TEST(TestEqualOrderOperators, Gradient)
 {
   FiniteElement fe = getFE();
 
   Matrix EM_scalar(2*2,2);
   // single component + pressure
-  WeakOperators::Gradient(EM_scalar, fe);
+  EqualOrderOperators::Weak::Gradient(EM_scalar, fe);
   const DoubleVec EM_scalar_ref = {{-1.0,-2.0},
                                    {-3.0,-6.0},
                                    {-2.0,-4.0},
@@ -95,7 +95,7 @@ TEST(TestWeakOperators, Gradient)
   check_matrix_equal(EM_scalar, EM_scalar_ref);
 
   Vector EV_scalar(4);
-  WeakOperators::Gradient(EV_scalar, fe);
+  EqualOrderOperators::Weak::Gradient(EV_scalar, fe);
   ASSERT_NEAR(EV_scalar(1),  fe.dNdX(1,1), 1e-13);
   ASSERT_NEAR(EV_scalar(2),  fe.dNdX(1,2), 1e-13);
   ASSERT_NEAR(EV_scalar(3),  fe.dNdX(2,1), 1e-13);
@@ -103,13 +103,13 @@ TEST(TestWeakOperators, Gradient)
 }
 
 
-TEST(TestWeakOperators, Laplacian)
+TEST(TestEqualOrderOperators, Laplacian)
 {
   FiniteElement fe = getFE();
 
   // single scalar block
   Matrix EM_scalar(2,2);
-  WeakOperators::Laplacian(EM_scalar, fe);
+  EqualOrderOperators::Weak::Laplacian(EM_scalar, fe);
 
   const DoubleVec EM_scalar_ref = {{10.0, 14.0},
                                    {14.0, 20.0}};
@@ -118,7 +118,7 @@ TEST(TestWeakOperators, Laplacian)
 
   // multiple (2) blocks in 3 component element matrix
   Matrix EM_multi(2*2,2*2);
-  WeakOperators::Laplacian(EM_multi, fe, 1.0);
+  EqualOrderOperators::Weak::Laplacian(EM_multi, fe, 1.0);
 
   const DoubleVec EM_multi_ref = {{10.0,  0.0, 14.0,  0.0},
                                   { 0.0, 10.0,  0.0, 14.0},
@@ -129,7 +129,7 @@ TEST(TestWeakOperators, Laplacian)
 
   // stress formulation
   Matrix EM_stress(2*2,2*2);
-  WeakOperators::Laplacian(EM_stress, fe, 1.0, true);
+  EqualOrderOperators::Weak::Laplacian(EM_stress, fe, 1.0, true);
   const DoubleVec EM_stress_ref = {{11.0,  3.0, 16.0,  6.0},
                                    { 3.0, 19.0,  4.0, 26.0},
                                    {16.0,  4.0, 24.0,  8.0},
@@ -138,25 +138,25 @@ TEST(TestWeakOperators, Laplacian)
   check_matrix_equal(EM_stress, EM_stress_ref);
 
   Matrix EM_coeff(2, 2);
-  WeakOperators::LaplacianCoeff(EM_coeff, fe.dNdX, fe, 2.0);
+  EqualOrderOperators::Weak::LaplacianCoeff(EM_coeff, fe.dNdX, fe, 2.0);
   const DoubleVec EM_coeff_ref = {{104.0, 148.0},
                                   {152.0, 216.0}};
   check_matrix_equal(EM_coeff, EM_coeff_ref);
 }
 
 
-TEST(TestWeakOperators, Mass)
+TEST(TestEqualOrderOperators, Mass)
 {
   FiniteElement fe = getFE();
 
   Matrix EM_scalar(2,2);
-  WeakOperators::Mass(EM_scalar, fe, 2.0);
+  EqualOrderOperators::Weak::Mass(EM_scalar, fe, 2.0);
   const DoubleVec EM_scalar_ref = {{2.0, 4.0},
                                    {4.0, 8.0}};
   check_matrix_equal(EM_scalar, EM_scalar_ref);
 
   Matrix EM_vec(2*2,2*2);
-  WeakOperators::Mass(EM_vec, fe);
+  EqualOrderOperators::Weak::Mass(EM_vec, fe);
 
   const DoubleVec EM_vec_ref = {{1.0, 0.0, 2.0, 0.0},
                                 {0.0, 1.0, 0.0, 2.0},
@@ -166,37 +166,37 @@ TEST(TestWeakOperators, Mass)
 }
 
 
-TEST(TestWeakOperators, Source)
+TEST(TestEqualOrderOperators, Source)
 {
   FiniteElement fe = getFE();
 
   Vector EV_scalar(2);
-  WeakOperators::Source(EV_scalar, fe, 2.0);
+  EqualOrderOperators::Weak::Source(EV_scalar, fe, 2.0);
   ASSERT_NEAR(EV_scalar(1),  2.0, 1e-13);
   ASSERT_NEAR(EV_scalar(2),  4.0, 1e-13);
 
   Vector EV_vec(4);
   Vec3 f;
   f[0] = 1.0; f[1] = 2.0;
-  WeakOperators::Source(EV_vec, fe, f, 1.0);
+  EqualOrderOperators::Weak::Source(EV_vec, fe, f, 1.0);
   ASSERT_NEAR(EV_vec(1),  1.0, 1e-13);
   ASSERT_NEAR(EV_vec(2),  2.0, 1e-13);
   ASSERT_NEAR(EV_vec(3),  2.0, 1e-13);
   ASSERT_NEAR(EV_vec(4),  4.0, 1e-13);
   EV_vec.fill(0.0);
-  WeakOperators::Source(EV_vec, fe, 2.0, 0);
+  EqualOrderOperators::Weak::Source(EV_vec, fe, 2.0, 0);
   ASSERT_NEAR(EV_vec(1),  2.0, 1e-13);
   ASSERT_NEAR(EV_vec(2),  2.0, 1e-13);
   ASSERT_NEAR(EV_vec(3),  4.0, 1e-13);
   ASSERT_NEAR(EV_vec(4),  4.0, 1e-13);
   EV_vec.fill(0.0);
-  WeakOperators::Source(EV_vec, fe, 2.0, 1);
+  EqualOrderOperators::Weak::Source(EV_vec, fe, 2.0, 1);
   ASSERT_NEAR(EV_vec(1),  2.0, 1e-13);
   ASSERT_NEAR(EV_vec(2),  0.0, 1e-13);
   ASSERT_NEAR(EV_vec(3),  4.0, 1e-13);
   ASSERT_NEAR(EV_vec(4),  0.0, 1e-13);
   EV_vec.fill(0.0);
-  WeakOperators::Source(EV_vec, fe, 2.0, 2);
+  EqualOrderOperators::Weak::Source(EV_vec, fe, 2.0, 2);
   ASSERT_NEAR(EV_vec(1),  0.0, 1e-13);
   ASSERT_NEAR(EV_vec(2),  2.0, 1e-13);
   ASSERT_NEAR(EV_vec(3),  0.0, 1e-13);

Apps/Common/WeakOperators.C

@@ -1,170 +0,0 @@
-//==============================================================================
-//!
-//! \file WeakOperators.C
-//!
-//! \date Jul 22 2015
-//!
-//! \author Arne Morten Kvarving / SINTEF
-//!
-//! \brief Various weak, discrete operators.
-//!
-//==============================================================================
-
-#include "WeakOperators.h"
-#include "FiniteElement.h"
-#include "Vec3.h"
-
-namespace WeakOperators {
-  //! \brief Helper for adding an element matrix to several components.
-  //! \param[out] EM The element matrix to add to.
-  //! \param[in] A The scalar element matrix to add.
-  //! \param[in] cmp Number of components to add matrix to
-  //! \param[in] nf Number of components in total matrix.
-  //! \param[in] scmp Index of first component to add matrix to.
-  static void addComponents(Matrix& EM, const Matrix& A,
-                            size_t cmp, size_t nf, size_t scmp)
-  {
-    if (cmp == 1 && nf == 1)
-      EM += A;
-    else
-      for (size_t i = 1; i <= A.rows(); ++i)
-        for (size_t j = 1; j <= A.cols(); ++j)
-          for (size_t k = 1; k <= cmp; ++k)
-            EM(nf*(i-1)+k+scmp,nf*(j-1)+k+scmp) += A(i, j);
-  }
-
-
-  //! \brief Helper applying a divergence (1) or a gradient (2) or both (0) operation
-  template<int Operation>
-  static void DivGrad(Matrix& EM, const FiniteElement& fe,
-              double scale, int basis, int tbasis)
-  {
-    size_t nsd = fe.grad(basis).cols();
-    for (size_t i = 1; i <= fe.basis(tbasis).size();++i)
-      for (size_t j = 1; j <= fe.basis(basis).size();++j)
-        for (size_t k = 1; k <= nsd; ++k) {
-          double div = fe.basis(basis)(j)*fe.grad(tbasis)(i,k)*fe.detJxW;
-          if (Operation == 2)
-            EM((i-1)*nsd+k,j) += -scale*div;
-          if (Operation == 1)
-            EM(j, (i-1)*nsd+k) += scale*div;
-        }
-  }
-
-
-  void Advection(Matrix& EM, const FiniteElement& fe,
-                 const Vec3& AC, double scale, int basis)
-  {
-    Matrix C(fe.basis(basis).size(), fe.basis(basis).size());
-    size_t ncmp = EM.rows() / C.rows();
-    for (size_t i = 1; i <= fe.basis(basis).size(); ++i) {
-      for (size_t j = 1; j <= fe.basis(basis).size(); ++j) {
-        // Sum convection for each direction
-        for (size_t k = 1; k <= fe.grad(basis).cols(); ++k)
-          C(i,j) += AC[k-1]*fe.dNdX(j,k);
-
-        C(i,j) *= scale*fe.N(i)*fe.detJxW;
-      }
-    }
-    addComponents(EM, C, ncmp, ncmp, 0);
-  }
-
-
-  void Divergence(Matrix& EM, const FiniteElement& fe, double scale,
-                  int basis, int tbasis)
-  {
-    DivGrad<1>(EM,fe,scale,basis,tbasis);
-  }
-
-
-  void Gradient(Matrix& EM, const FiniteElement& fe, double scale,
-                int basis, int tbasis)
-  {
-    DivGrad<2>(EM,fe,scale,basis,tbasis);
-  }
-
-
-  void Divergence(Vector& EV, const FiniteElement& fe,
-                  const Vec3& D, double scale, int basis)
-  {
-    for (size_t i = 1; i <= fe.basis(basis).size(); ++i) {
-      double div=0.0;
-      for (size_t k = 1; k <= fe.grad(basis).cols(); ++k)
-        div += D[k-1]*fe.grad(basis)(i,k);
-      EV(i) += scale*div*fe.detJxW;
-    }
-  }
-
-
-  void Gradient(Vector& EV, const FiniteElement& fe,
-                double scale, int basis)
-  {
-    size_t nsd = fe.grad(basis).cols();
-    for (size_t i = 1; i <= fe.basis(basis).size(); ++i)
-      for (size_t k = 1; k <= nsd; ++k)
-        EV((i-1)*nsd+k) += scale*fe.grad(basis)(i,k)*fe.detJxW;
-  }
-
-
-  void Laplacian(Matrix& EM, const FiniteElement& fe,
-                 double scale, bool stress, int basis)
-  {
-    size_t cmp = EM.rows() / fe.basis(basis).size();
-    Matrix A;
-    A.multiply(fe.grad(basis),fe.grad(basis),false,true);
-    A *= scale*fe.detJxW;
-    addComponents(EM, A, cmp, cmp, 0);
-    if (stress)
-      for (size_t i = 1; i <= fe.basis(basis).size(); i++)
-        for (size_t j = 1; j <= fe.basis(basis).size(); j++)
-          for (size_t k = 1; k <= cmp; k++)
-            for (size_t l = 1; l <= cmp; l++)
-              EM(cmp*(j-1)+k,cmp*(i-1)+l) += scale*fe.grad(basis)(i,k)*fe.grad(basis)(j,l)*fe.detJxW;
-  }
-
-
-  void LaplacianCoeff(Matrix& EM, const Matrix& K,
-                      const FiniteElement& fe,
-                      double scale, int basis)
-  {
-    Matrix KB;
-    KB.multiply(K,fe.grad(basis),false,true).multiply(scale*fe.detJxW);
-    EM.multiply(fe.grad(basis),KB,false,false,true);
-  }
-
-
-  void Mass(Matrix& EM, const FiniteElement& fe,
-            double scale, int basis)
-  {
-    size_t ncmp = EM.rows()/fe.basis(basis).size();
-    Matrix A;
-    A.outer_product(fe.basis(basis),fe.basis(basis),false);
-    A *= scale*fe.detJxW;
-    addComponents(EM, A, ncmp, ncmp, 0);
-  }
-
-
-  void Source(Vector& EV, const FiniteElement& fe,
-              double scale, int cmp, int basis)
-  {
-    size_t ncmp = EV.size() / fe.basis(basis).size();
-    if (cmp == 1 && ncmp == 1)
-      EV.add(fe.basis(basis), scale*fe.detJxW);
-    else {
-      for (size_t i = 1; i <= fe.basis(basis).size(); ++i)
-        for (size_t k  = (cmp == 0 ? 1: cmp);
-                    k <= (cmp == 0 ? ncmp : cmp); ++k)
-          EV(ncmp*(i-1)+k) += scale*fe.basis(basis)(i)*fe.detJxW;
-    }
-  }
-
-
-  void Source(Vector& EV, const FiniteElement& fe,
-              const Vec3& f, double scale, int basis)
-  {
-    size_t cmp = EV.size() / fe.basis(basis).size();
-    for (size_t i = 1; i <= fe.basis(basis).size(); ++i)
-      for (size_t k = 1; k <= cmp; ++k)
-        EV(cmp*(i-1)+k) += scale*f[k-1]*fe.basis(basis)(i)*fe.detJxW;
-  }
-}

Apps/Common/WeakOperators.h

@@ -1,145 +0,0 @@
-//==============================================================================
-//!
-//! \file WeakOperators.h
-//!
-//! \date Jul 22 2015
-//!
-//! \author Arne Morten Kvarving / SINTEF
-//!
-//! \brief Various weak, discrete operators.
-//!
-//==============================================================================
-
-#ifndef WEAKOPERATORS_H_
-#define WEAKOPERATORS_H_
-
-class Vec3;
-
-#include "FiniteElement.h"
-#include "MatVec.h"
-
-namespace WeakOperators
-{
-  //! \brief Compute an advection term.
-  //! \param[out] EM The element matrix to add contribution to
-  //! \param[in] fe The finite element to evaluate for
-  //! \param[in] AC Advecting field
-  //! \param[in] scale Scaling factor for contribution
-  //! \param[in] basis Basis to use
-  void Advection(Matrix& EM, const FiniteElement& fe,
-                 const Vec3& AC, double scale=1.0, int basis=1);
-
-  //! \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
-  template<class T>
-  void Convection(Matrix& EM, const FiniteElement& fe,
-                  const Vec3& U, const T& dUdX, double scale,
-                  bool conservative=false, int basis=1)
-  {
-    size_t cmp = EM.rows() / fe.basis(basis).size();
-    if (conservative) {
-      Advection(EM, fe, U, -scale, basis);
-      for (size_t i = 1;i <= fe.basis(basis).size();i++)
-        for (size_t j = 1;j <= fe.basis(basis).size();j++) {
-          for (size_t k = 1;k <= cmp;k++) {
-            for (size_t l = 1;l <= cmp;l++)
-              EM((j-1)*cmp+l,(i-1)*cmp+k) -= scale*U[l-1]*fe.basis(basis)(i)*fe.grad(basis)(j,k)*fe.detJxW;
-          }
-        }
-    }
-    else {
-      Advection(EM, fe, U, scale, basis);
-      for (size_t i = 1;i <= fe.basis(basis).size();i++)
-        for (size_t j = 1;j <= fe.basis(basis).size();j++) {
-          for (size_t k = 1;k <= cmp;k++) {
-            for (size_t l = 1;l <= cmp;l++)
-              EM((j-1)*cmp+l,(i-1)*cmp+k) += scale*dUdX(l,k)*fe.basis(basis)(j)*fe.basis(basis)(i)*fe.detJxW;
-          }
-        }
-    }
-  }
-
-  //! \brief Compute a divergence term.
-  //! \param[out] EM The element matrix to add contribution to
-  //! \param[in] fe The finite element to evaluate for
-  //! \param[in] scale Scaling factor for contribution
-  //! \param[in] basis Basis for field
-  //! \param[in] tbasis Test function basis
-  void Divergence(Matrix& EM, const FiniteElement& fe,
-                  double scale=1.0, int basis=1, int tbasis=1);
-
-  //! \brief Compute a divergence term.
-  //! \param[out] EV The element vector to add contribution to
-  //! \param[in] fe The finite element to evaluate for
-  //! \param[in] D Divergence of field
-  //! \param[in] scale Scaling factor for contribution
-  //! \param[in] basis Test function basis
-  void Divergence(Vector& EV, const FiniteElement& fe,
-                  const Vec3& D, double scale=1.0, int basis=1);
-
-  //! \brief Compute a gradient term for a (potentially mixed) vector/scalar field.
-  //! \param[out] EM The element matrix to add contribution to
-  //! \param[in] fe The finite element to evaluate for
-  //! \param[in] scale Scaling factor for contribution
-  //! \param[in] basis Basis for field
-  //! \param[in] tbasis Test function basis
-  void Gradient(Matrix& EM, const FiniteElement& fe,
-                double scale=1.0, int basis=1, int tbasis=1);
-
-  //! \brief Compute a gradient term.
-  //! \param[out] EV The element vector to add contribution to
-  //! \param[in] fe The finite element to evaluate for
-  //! \param[in] scale Scaling factor for contribution
-  //! \param[in] tbasis Test function basis
-  void Gradient(Vector& EV, const FiniteElement& fe,
-                double scale=1.0, int basis=1);
-
-  //! \brief Compute a laplacian.
-  //! \param[out] EM The element matrix to add contribution to
-  //! \param[in] fe The finite element to evaluate for
-  //! \param[in] scale Scaling factor for contribution
-  //! \param[in] stress Whether to add extra stress formulation terms
-  //! \param[in] basis Basis to use
-  void Laplacian(Matrix& EM, const FiniteElement& fe,
-                 double scale=1.0, bool stress=false, int basis=1);
-
-  //! \brief Compute a heteregenous coefficient laplacian.
-  //! \param[out] EM The element matrix to add contribution to
-  //! \param[out] K The coefficient matrix
-  //! \param[in] fe The finite element to evaluate for
-  //! \param[in] scale Scaling factor for contribution
-  void LaplacianCoeff(Matrix& EM, const Matrix& K, const FiniteElement& fe,
-                      double scale=1.0, int basis=1);
-
-  //! \brief Compute a mass term.
-  //! \param[out] EM The element matrix to add contribution to
-  //! \param[in] fe The finite element to evaluate for
-  //! \param[in] scale Scaling factor for contribution
-  //! \param[in] basis Basis to use
-  void Mass(Matrix& EM, const FiniteElement& fe,
-            double scale=1.0, int basis=1);
-
-  //! \brief Compute a source term.
-  //! \param[out] EV The element vector to add contribution to
-  //! \param[in] fe The finite element to evaluate for
-  //! \param[in] scale Scaling factor for contribution
-  //! \param[in] basis Basis to use
-  //! \param[in] cmp Component to add (0 for all)
-  void Source(Vector& EV, const FiniteElement& fe,
-              double scale=1.0, int cmp=1, int basis=1);
-
-  //! \brief Compute a vector-source term.
-  //! \param[out] EV The element vector to add contribution to
-  //! \param[in] fe The finite element to evaluate for
-  //! \param[in] scale Vector with contributions
-  //! \param[in] scale Scaling factor for contribution
-  //! \param[in] basis Basis to use
-  void Source(Vector& EV, const FiniteElement& fe,
-              const Vec3& f, double scale=1.0, int basis=1);
-}
-
-#endif