diff --git a/Apps/Common/EqualOrderOperators.C b/Apps/Common/EqualOrderOperators.C index d3b5c39c..b1be3dee 100644 --- a/Apps/Common/EqualOrderOperators.C +++ b/Apps/Common/EqualOrderOperators.C @@ -13,6 +13,7 @@ #include "EqualOrderOperators.h" #include "FiniteElement.h" #include "Vec3.h" +#include "Vec3Oper.h" //! \brief Helper for adding an element matrix to several components. //! \param[out] EM The element matrix to add to. @@ -211,6 +212,19 @@ void EqualOrderOperators::Weak::Source(Vector& EV, const FiniteElement& fe, } +void EqualOrderOperators::Residual::Advection(Vector& EV, const FiniteElement& fe, + const Vec3& AC, const Tensor& g, + double scale, int basis) +{ + size_t nsd = fe.grad(basis).cols(); + for (size_t k = 1; k <= nsd; ++k) { + double ag = g[k-1]*AC; + for (size_t i = 1; i <= fe.basis(basis).size(); ++i) + EV((i-1)*nsd+k) = ag*scale*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, diff --git a/Apps/Common/EqualOrderOperators.h b/Apps/Common/EqualOrderOperators.h index 4df1b06e..239e073b 100644 --- a/Apps/Common/EqualOrderOperators.h +++ b/Apps/Common/EqualOrderOperators.h @@ -142,6 +142,16 @@ public: //! \brief Common weak residual operators using equal-ordered discretizations. class Residual { public: + //! \brief Compute an advection term. + //! \param[out] EM The element vector 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(Vector& EV, const FiniteElement& fe, + const Vec3& AC, const Tensor& g, + double scale = 1.0, int basis=1); + //! \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