NieVelocityImposition.h

// Copyright (C) 2015 Technische Universitaet Muenchen
// This file is part of the Mamico project. For conditions of distribution
// and use, please see the copyright notice in Mamico's main folder, or at
// www5.in.tum.de/mamico
#ifndef _MOLECULARDYNAMICS_COUPLING_NIEVELOCITYIMPOSITION_H_
#define _MOLECULARDYNAMICS_COUPLING_NIEVELOCITYIMPOSITION_H_
 
#include "coupling/MomentumInsertion.h"
#include "coupling/cell-mappings/ComputeAvgForceAndVelocity.h"
#include "coupling/cell-mappings/NieVelocityImpositionMapping.h"
#include "coupling/interface/MDSolverInterface.h"
 
namespace coupling {
template <class LinkedCell, unsigned int dim> class NieVelocityImposition;
}

template <class LinkedCell, unsigned int dim> class coupling::NieVelocityImposition : public coupling::MomentumInsertion<LinkedCell, dim> {
public:
  NieVelocityImposition(coupling::interface::MDSolverInterface<LinkedCell, dim>* const mdSolverInterface, const unsigned int& outermostLayer,
                        const unsigned int& innermostLayer, const tarch::la::Vector<2 * dim, bool> impositionEnabled)
      : coupling::MomentumInsertion<LinkedCell, dim>(mdSolverInterface), _outermostLayer(outermostLayer), _innermostLayer(innermostLayer),
        _impositionEnabled(impositionEnabled), _enableInnerImposition(false) {}

  virtual ~NieVelocityImposition() {}

  unsigned int getTimeIntervalPerMomentumInsertion() const override { return 1; }

  void insertMomentum(coupling::datastructures::CouplingCellWithLinkedCells<LinkedCell, dim>& cell, I02 idx) const override {
    // nop if this is not an imposition cell
    if (!isInsideImpositionLayer(idx)) {
      return;
    }
    // set continuum velocity
    tarch::la::Vector<dim, double> continuumVelocity(cell.getMicroscopicMomentum());
 
    coupling::cellmappings::ComputeAvgForceAndVelocity<LinkedCell, dim> computeForceAndVelocity(
        coupling::MomentumInsertion<LinkedCell, dim>::_mdSolverInterface);
    cell.iterateConstCells(computeForceAndVelocity);
    const tarch::la::Vector<dim, double> avgVel(computeForceAndVelocity.getAvgVelocity());
    const tarch::la::Vector<dim, double> avgF(computeForceAndVelocity.getAvgForce());
 
    coupling::cellmappings::NieVelocityImpositionMapping<LinkedCell, dim> velocityImposition(continuumVelocity, avgVel, avgF,
                                                                                             coupling::MomentumInsertion<LinkedCell, dim>::_mdSolverInterface);
    cell.iterateCells(velocityImposition);
  }
 
  void setInnerImposition(bool enable) override { _enableInnerImposition = enable; }
 
private:
  bool isInsideImpositionLayer(I01 globalCellIndex) const {
    bool inner = true;
    bool outer = false;
    tarch::la::Vector<dim, unsigned int> globalIndexUnsigned{globalCellIndex.get()};
    for (unsigned int d = 0; d < dim; d++) {
      if (_impositionEnabled[d * 2]) { // negative direction
        inner = inner && (globalIndexUnsigned[d] > _innermostLayer);
        outer = outer || (globalIndexUnsigned[d] < _outermostLayer);
      }
      if (_impositionEnabled[d * 2 + 1]) { // positive direction
        inner = inner && (globalIndexUnsigned[d] < 1 + I09::numberCellsInDomain[d] - _innermostLayer);
        outer = outer || (globalIndexUnsigned[d] > 1 + I09::numberCellsInDomain[d] - _outermostLayer);
      }
    }
    return (_enableInnerImposition || !inner) && !outer;
  }

  const unsigned int _outermostLayer;
  const unsigned int _innermostLayer;
  const tarch::la::Vector<2 * dim, bool> _impositionEnabled;
  bool _enableInnerImposition;
};
 
#endif // _MOLECULARDYNAMICS_COUPLING_VELOCITYGRADIENTRELAXATION_H_