1 files changed, 427 insertions, 0 deletions
diff --git a/src/boundary/momentaOnBoundaries.hh b/src/boundary/momentaOnBoundaries.hh
new file mode 100644
index 0000000..3ec9f94
--- /dev/null
+++ b/src/boundary/momentaOnBoundaries.hh
@@ -0,0 +1,427 @@
+/*  This file is part of the OpenLB library
+ *
+ *  Copyright (C) 2006, 2007 Jonas Latt
+ *  E-mail contact: info@openlb.net
+ *  The most recent release of OpenLB can be downloaded at
+ *  <http://www.openlb.net/>
+ *
+ *  This program is free software; you can redistribute it and/or
+ *  modify it under the terms of the GNU General Public License
+ *  as published by the Free Software Foundation; either version 2
+ *  of the License, or (at your option) any later version.
+ *
+ *  This program is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public
+ *  License along with this program; if not, write to the Free
+ *  Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
+ *  Boston, MA  02110-1301, USA.
+*/
+
+/** \file
+ * Implementation of boundary cell dynamics -- generic implementation.
+ */
+#ifndef MOMENTA_ON_BOUNDARIES_HH
+#define MOMENTA_ON_BOUNDARIES_HH
+
+#include "momentaOnBoundaries.h"
+#include "core/util.h"
+#include "dynamics/lbHelpers.h"
+
+namespace olb {
+
+////////////////////// Class EquilibriumBM //////////////////////
+
+template<typename T, typename DESCRIPTOR>
+EquilibriumBM<T,DESCRIPTOR>::EquilibriumBM()
+{
+  _rho = (T)1;
+  for (int iD=0; iD<DESCRIPTOR::d; ++iD) {
+    _u[iD] = T();
+  }
+}
+
+template<typename T, typename DESCRIPTOR>
+EquilibriumBM<T,DESCRIPTOR>::EquilibriumBM(T rho, const T u[DESCRIPTOR::d])
+{
+  _rho = rho;
+  for (int iD=0; iD<DESCRIPTOR::d; ++iD) {
+    _u[iD] = u[iD];
+  }
+}
+
+template<typename T, typename DESCRIPTOR>
+T EquilibriumBM<T,DESCRIPTOR>::computeRho( Cell<T,DESCRIPTOR> const& cell ) const
+{
+  return _rho;
+}
+
+template<typename T, typename DESCRIPTOR>
+void EquilibriumBM<T,DESCRIPTOR>::computeU(
+  Cell<T,DESCRIPTOR> const& cell, T u[DESCRIPTOR::d]) const
+{
+  for (int iD=0; iD<DESCRIPTOR::d; ++iD) {
+    u[iD] = _u[iD];
+  }
+}
+
+template<typename T, typename DESCRIPTOR>
+void EquilibriumBM<T,DESCRIPTOR>::computeJ (
+  Cell<T,DESCRIPTOR> const& cell, T j[DESCRIPTOR::d]) const
+{
+  for (int iD=0; iD<DESCRIPTOR::d; ++iD) {
+    j[iD] = _u[iD]*_rho;
+  }
+}
+
+template<typename T, typename DESCRIPTOR>
+void EquilibriumBM<T,DESCRIPTOR>::computeStress( Cell<T,DESCRIPTOR> const& cell,
+    T rho, const T u[DESCRIPTOR::d], T pi[util::TensorVal<DESCRIPTOR >::n] ) const
+{
+  for (int iPi=0; iPi<util::TensorVal<DESCRIPTOR >::n; ++iPi) {
+    pi[iPi] = T();
+  }
+}
+
+
+template<typename T, typename DESCRIPTOR>
+void EquilibriumBM<T,DESCRIPTOR>::defineRho( Cell<T,DESCRIPTOR>& cell, T rho )
+{
+  _rho = rho;
+}
+
+template<typename T, typename DESCRIPTOR>
+void EquilibriumBM<T,DESCRIPTOR>::defineU(
+  Cell<T,DESCRIPTOR>& cell, const T u[DESCRIPTOR::d])
+{
+  for (int iD=0; iD<DESCRIPTOR::d; ++iD) {
+    _u[iD] = u[iD];
+  }
+}
+
+template<typename T, typename DESCRIPTOR>
+void EquilibriumBM<T,DESCRIPTOR>::defineAllMomenta( Cell<T,DESCRIPTOR>& cell,
+    T rho, const T u[DESCRIPTOR::d], const T pi[util::TensorVal<DESCRIPTOR >::n] )
+{
+  defineRho(cell, rho);
+  defineU(cell, u);
+}
+
+
+
+////////////////////// Class VelocityBM //////////////////////
+
+template<typename T, typename DESCRIPTOR, int direction, int orientation>
+VelocityBM<T,DESCRIPTOR,direction,orientation>::VelocityBM()
+{
+  for (int iD=0; iD<DESCRIPTOR::d; ++iD) {
+    _u[iD] = T();
+  }
+}
+
+/** It takes as argument the value of the velocity to be imposed on the
+ * boundary.
+ */
+template<typename T, typename DESCRIPTOR, int direction, int orientation>
+VelocityBM<T,DESCRIPTOR,direction,orientation>::VelocityBM(const T u[DESCRIPTOR::d])
+{
+  for (int iD=0; iD<DESCRIPTOR::d; ++iD) {
+    _u[iD] = u[iD];
+  }
+}
+
+template<typename T, typename DESCRIPTOR, int direction, int orientation>
+T VelocityBM<T,DESCRIPTOR,direction,orientation>::computeRho( Cell<T,DESCRIPTOR> const& cell ) const
+{
+  return velocityBMRho<T, DESCRIPTOR, direction, orientation>(cell, _u);
+}
+
+template<typename T, typename DESCRIPTOR, int direction, int orientation>
+void VelocityBM<T,DESCRIPTOR,direction,orientation>::computeU (
+  Cell<T,DESCRIPTOR> const& cell, T u[DESCRIPTOR::d]) const
+{
+  computeU(u);
+}
+
+template<typename T, typename DESCRIPTOR, int direction, int orientation>
+void VelocityBM<T,DESCRIPTOR,direction,orientation>::computeJ (
+  Cell<T,DESCRIPTOR> const& cell, T j[DESCRIPTOR::d]) const
+{
+  T rho = computeRho(cell);
+  for (int iD=0; iD<DESCRIPTOR::d; ++iD) {
+    j[iD] = _u[iD]*rho;
+  }
+}
+
+template<typename T, typename DESCRIPTOR, int direction, int orientation>
+void VelocityBM<T,DESCRIPTOR,direction,orientation>::computeU( T u[DESCRIPTOR::d] ) const
+{
+  for (int iD=0; iD<DESCRIPTOR::d; ++iD) {
+    u[iD] = _u[iD];
+  }
+}
+
+template<typename T, typename DESCRIPTOR, int direction, int orientation>
+void VelocityBM<T,DESCRIPTOR,direction,orientation>::defineRho(Cell<T,DESCRIPTOR>& cell, T rho )
+{
+}
+
+template<typename T, typename DESCRIPTOR, int direction, int orientation>
+void VelocityBM<T,DESCRIPTOR,direction,orientation>::defineU (
+  Cell<T,DESCRIPTOR>& cell, const T u[DESCRIPTOR::d])
+{
+  defineU(u);
+}
+
+template<typename T, typename DESCRIPTOR, int direction, int orientation>
+void VelocityBM<T,DESCRIPTOR,direction,orientation>::defineU(const T u[DESCRIPTOR::d])
+{
+  for (int iD=0; iD<DESCRIPTOR::d; ++iD) {
+    _u[iD] = u[iD];
+  }
+}
+
+template<typename T, typename DESCRIPTOR, int direction, int orientation>
+void VelocityBM<T,DESCRIPTOR,direction,orientation>::defineAllMomenta (
+  Cell<T,DESCRIPTOR>& cell, T rho, const T u[DESCRIPTOR::d], const T pi[util::TensorVal<DESCRIPTOR >::n] )
+{
+  this->defineRhoU(cell, rho, u);
+}
+
+
+////////////////////// Class PressureBM //////////////////////
+
+template<typename T, typename DESCRIPTOR, int direction, int orientation>
+PressureBM<T,DESCRIPTOR,direction,orientation>::PressureBM()
+{
+  for (int iD=0; iD<DESCRIPTOR::d; ++iD) {
+    _values[iD] = T();
+  }
+  _values[direction] = 1.;
+}
+
+/** It takes as argument the value of the tangential velocity
+ * components, and the pressure, to be imposed on the boundary.
+ */
+template<typename T, typename DESCRIPTOR, int direction, int orientation>
+PressureBM<T,DESCRIPTOR,direction,orientation>::PressureBM(const T values[DESCRIPTOR::d])
+{
+  for (int iD=0; iD<DESCRIPTOR::d; ++iD) {
+    _values[iD] = values[iD];
+  }
+}
+
+template<typename T, typename DESCRIPTOR, int direction, int orientation>
+T PressureBM<T,DESCRIPTOR,direction,orientation>::computeRho( Cell<T,DESCRIPTOR> const& cell ) const
+{
+  return computeRho();
+}
+
+template<typename T, typename DESCRIPTOR, int direction, int orientation>
+T PressureBM<T,DESCRIPTOR,direction,orientation>::computeRho() const
+{
+  return _values[direction];
+}
+
+template<typename T, typename DESCRIPTOR, int direction, int orientation>
+void PressureBM<T,DESCRIPTOR,direction,orientation>::computeU (
+  Cell<T,DESCRIPTOR> const& cell, T u[DESCRIPTOR::d]) const
+{
+  for (int iD=0; iD<DESCRIPTOR::d; ++iD) {
+    u[iD] = _values[iD];
+  }
+  T rho = _values[direction];
+
+  std::vector<int> const& onWallIndices
+    = util::subIndex<DESCRIPTOR, direction, 0>();
+
+  std::vector<int> const& normalIndices
+    = util::subIndex<DESCRIPTOR, direction, orientation>();
+
+  T rhoOnWall = T();
+  for (auto & e : onWallIndices) {
+    rhoOnWall += cell[e];
+  }
+
+  T rhoNormal = T();
+  for (auto & e : normalIndices) {
+    rhoNormal += cell[e];
+  }
+
+  u[direction] = (T)orientation*( ((T)2*rhoNormal+rhoOnWall+(T)1 ) / rho-(T)1 );
+}
+
+template<typename T, typename DESCRIPTOR, int direction, int orientation>
+void PressureBM<T,DESCRIPTOR,direction,orientation>::computeJ (
+  Cell<T,DESCRIPTOR> const& cell, T j[DESCRIPTOR::d]) const
+{
+  computeU(cell, j);
+  T rho = computeRho(cell);
+  for (int iD=0; iD<DESCRIPTOR::d; ++iD) {
+    j[iD] *= rho;
+  }
+}
+
+template<typename T, typename DESCRIPTOR, int direction, int orientation>
+void PressureBM<T,DESCRIPTOR,direction,orientation>::defineRho(Cell<T,DESCRIPTOR>& cell, T rho)
+{
+  defineRho(rho);
+}
+
+template<typename T, typename DESCRIPTOR, int direction, int orientation>
+void PressureBM<T,DESCRIPTOR,direction,orientation>::defineRho(T rho )
+{
+  _values[direction] = rho;
+}
+
+template<typename T, typename DESCRIPTOR, int direction, int orientation>
+void PressureBM<T,DESCRIPTOR,direction,orientation>::defineU(Cell<T,DESCRIPTOR>& cell, const T u[DESCRIPTOR::d])
+{
+  for (int iD=0; iD<DESCRIPTOR::d; ++iD) {
+    if (iD != direction) {
+      _values[iD] = u[iD];
+    }
+  }
+}
+
+template<typename T, typename DESCRIPTOR, int direction, int orientation>
+void PressureBM<T,DESCRIPTOR,direction,orientation>::defineAllMomenta(
+  Cell<T,DESCRIPTOR>& cell, T rho, const T u[DESCRIPTOR::d], const T pi[util::TensorVal<DESCRIPTOR >::n] )
+{
+  this->defineRhoU(cell, rho, u);
+}
+
+
+////////  FreeStressBM //////////////////////////////////////////////
+
+template<typename T, typename DESCRIPTOR>
+void FreeStressBM<T,DESCRIPTOR>::computeStress (
+  Cell<T,DESCRIPTOR> const& cell, T rho, const T u[DESCRIPTOR::d], T pi[util::TensorVal<DESCRIPTOR >::n] ) const
+{
+  lbHelpers<T,DESCRIPTOR>::computeStress(cell, rho, u, pi);
+}
+
+
+
+////////////////////// Class RegularizedBM //////////////////////
+
+template<typename T, typename DESCRIPTOR, int direction, int orientation>
+void RegularizedBM<T,DESCRIPTOR,direction,orientation>::computeStress (
+  Cell<T,DESCRIPTOR> const& cell, T rho, const T u[DESCRIPTOR::d],
+  T pi[util::TensorVal<DESCRIPTOR >::n] ) const
+{
+  BoundaryHelpers<T,DESCRIPTOR,direction,orientation>::computeStress(cell, rho, u, pi);
+}
+
+////////////////////// Class FixedVelocityBM //////////////////////////
+
+template<typename T, typename DESCRIPTOR>
+T FixedVelocityBM<T,DESCRIPTOR>::computeRho(Cell<T,DESCRIPTOR> const& cell) const
+{
+  return _basicMomenta.computeRho(cell);
+}
+
+template<typename T, typename DESCRIPTOR>
+void FixedVelocityBM<T,DESCRIPTOR>::computeU(Cell<T,DESCRIPTOR> const& cell, T u[DESCRIPTOR::d]) const
+{
+  for (int iD=0; iD<DESCRIPTOR::d; ++iD) {
+    u[iD] = _fixU[iD];
+  }
+}
+
+template<typename T, typename DESCRIPTOR>
+void FixedVelocityBM<T,DESCRIPTOR>::computeJ(Cell<T,DESCRIPTOR> const& cell, T j[DESCRIPTOR::d]) const
+{
+  T rho = computeRho(cell);
+  for (int iD=0; iD<DESCRIPTOR::d; ++iD) {
+    j[iD] = _fixU[iD]*rho;
+  }
+}
+
+template<typename T, typename DESCRIPTOR>
+void FixedVelocityBM<T,DESCRIPTOR>::computeStress (
+  Cell<T,DESCRIPTOR> const& cell, T rho, const T u[DESCRIPTOR::d], T pi[util::TensorVal<DESCRIPTOR >::n] ) const
+{
+  _basicMomenta.computeStress(cell, rho, u, pi);
+}
+
+template<typename T, typename DESCRIPTOR>
+void FixedVelocityBM<T,DESCRIPTOR>::computeRhoU (
+  Cell<T,DESCRIPTOR> const& cell, T& rho, T u[DESCRIPTOR::d] ) const
+{
+  rho = computeRho(cell);
+  computeU(cell,u);
+}
+
+template<typename T, typename DESCRIPTOR>
+void FixedVelocityBM<T,DESCRIPTOR>::computeAllMomenta (
+  Cell<T,DESCRIPTOR> const& cell, T& rho, T u[DESCRIPTOR::d], T pi[util::TensorVal<DESCRIPTOR >::n] ) const
+{
+  _basicMomenta.computeAllMomenta(cell, rho, u, pi);
+  computeU(cell, u);
+}
+
+template<typename T, typename DESCRIPTOR>
+void FixedVelocityBM<T,DESCRIPTOR>::defineRho(Cell<T,DESCRIPTOR>& cell, T rho)
+{
+  _basicMomenta.defineRho(cell, rho);
+}
+
+template<typename T, typename DESCRIPTOR>
+void FixedVelocityBM<T,DESCRIPTOR>::defineU(Cell<T,DESCRIPTOR>& cell, const T u[DESCRIPTOR::d])
+{
+  for (int iD=0; iD<DESCRIPTOR::d; ++iD) {
+    _fixU[iD] = u[iD];
+  }
+}
+
+template<typename T, typename DESCRIPTOR>
+void FixedVelocityBM<T,DESCRIPTOR>::defineRhoU(Cell<T,DESCRIPTOR>& cell, T rho, const T u[DESCRIPTOR::d])
+{
+  defineRho(cell,rho);
+  defineU(cell,u);
+}
+
+template<typename T, typename DESCRIPTOR>
+void FixedVelocityBM<T,DESCRIPTOR>::defineAllMomenta( Cell<T,DESCRIPTOR>& cell, T rho, const T u[DESCRIPTOR::d],
+    const T pi[util::TensorVal<DESCRIPTOR >::n] )
+{
+  _basicMomenta.defineAllMomenta(cell, rho, u, pi);
+  defineU(cell,u);
+}
+
+////////////////////// Extracted helper functions //////////////////////////
+
+template<typename T, typename DESCRIPTOR, int direction, int orientation>
+T velocityBMRho( Cell<T,DESCRIPTOR> const& cell, const T* u )
+{
+  std::vector<int> const& onWallIndices
+    = util::subIndex<DESCRIPTOR, direction, 0>();
+
+  std::vector<int> const& normalIndices
+    = util::subIndex<DESCRIPTOR, direction, orientation>();
+
+  T rhoOnWall = T();
+  for (auto & e : onWallIndices) {
+    rhoOnWall += cell[e];
+  }
+
+  T rhoNormal = T();
+  for (auto & e : normalIndices) {
+    rhoNormal += cell[e];
+  }
+
+  T rho =((T)2*rhoNormal+rhoOnWall+(T)1) /
+         ((T)1+(T)orientation*u[direction]);
+
+  return rho;
+}
+
+
+
+}  // namespace olb
+
+#endif