text: /*
   * The MIT License
   *
   * Copyright (c) 2015-2017 The University of Utah
   *
    * Permission is hereby granted, free of charge, to any person obtaining a copy
   * of this software and associated documentation files (the "Software"), to
   * deal in the Software without restriction, including without limitation the
   * rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
   * sell copies of the Software, and to permit persons to whom the Software is
   * furnished to do so, subject to the following conditions:
   *
   * The above copyright notice and this permission notice shall be included in
   * all copies or substantial portions of the Software.
   *
   * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
   * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
   * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
   * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
   * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
   * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS
   * IN THE SOFTWARE.
   */
  
  #ifndef BoundaryConditionExpression_Expr_h
  #define BoundaryConditionExpression_Expr_h
  
  #include <expression/Expression.h>
  #include <spatialops/structured/stencil/FVStaggeredBCOp.h>
  #include <spatialops/OperatorDatabase.h>
  #include <boost/shared_ptr.hpp>
  
  namespace Expr{
  
  /**
   *  \class BoundaryConditionExpression
   *  \tparam FieldT the type of field associated with this expression
   *  \tparam DirT the direction type for the boundary condition (SpatialOps::XDIR, SpatialOps::YDIR, SpatialOps::ZDIR).
   */
  template< typename FieldT, typename DirT >
  class BoundaryConditionExpression : public Expr::Expression<FieldT>
  {
  private:
    typedef SpatialOps::BCOpTypesFromDirection<FieldT,DirT> BCOpT;
  
    void err_check( const SpatialOps::BCType bcType, const Tag& tag )
    {
      if( bcType != SpatialOps::DIRICHLET && bcType != SpatialOps::NEUMANN ){
        std::ostringstream msg;
        msg << "ERROR: An unsupported boundary type was specified for expression '" << tag << "'\n"
            << "\tMust use either 'SpatialOps::DIRICHLET' or 'SpatialOps::NEUMANN'\n"
            << __FILE__ << " : " << __LINE__ << std::endl;
        throw std::invalid_argument( msg.str() );
      }
    }
  
  public:
  
    typedef typename BCOpT::Dirichlet  DirichletOpT;
    typedef typename BCOpT::Neumann    NeumannOpT;
  
    typedef typename DirichletOpT::OperatorType::DestFieldType BCValFieldT;
  
    typedef SpatialOps::SpatialMask<BCValFieldT>  MaskType;
    typedef boost::shared_ptr<MaskType>           MaskPtr;
  
  protected:
  
    const SpatialOps::BCType bcType_;
    const SpatialOps::BCSide side_;
    const MaskPtr mask_;
  
    const DirichletOpT* dirichletOp_;
    const NeumannOpT* neumannOp_;
  
  public:
  
    /**
     * \brief obtain the mask to apply this expression on.
     */
    const MaskType& get_mask() const{ return *mask_; }
  
    /**
     * \brief reset the mask on this expression
     */
    void reset_mask( MaskPtr mask ){ mask_ = mask; }
  
    BoundaryConditionExpression( MaskPtr mask,
                                 const SpatialOps::BCType bcType,
                                 const SpatialOps::BCSide side )
    : Expr::Expression<FieldT>(),
      bcType_( bcType ),
      side_( side ),
      mask_( mask )
    {
      err_check( bcType, this->get_tag() );
    }
  
    virtual ~BoundaryConditionExpression(){}
  
    void bind_operators( const SpatialOps::OperatorDatabase& opDB )
    {
      switch( bcType_ ){
        case SpatialOps::DIRICHLET: dirichletOp_ = opDB.retrieve_operator<DirichletOpT>(); break;
        case SpatialOps::NEUMANN  : neumannOp_   = opDB.retrieve_operator<NeumannOpT  >(); break;
        case SpatialOps::UNSUPPORTED_BCTYPE:
        default: break; // error trapping was already handled
      }
    }
  };
  
  //--------------------------------------------------------------------
  
  /**
   * \class ConstantBCOpExpression
   * \author James C. Sutherland
   * \date November 2014
   * \brief Provides support for operator inversion BCs on masked fields.
   */
  template<typename FieldT,typename DirT>
  class ConstantBCOpExpression : public BoundaryConditionExpression<FieldT,DirT>
  {
    const double value_;
  
    ConstantBCOpExpression( typename BoundaryConditionExpression<FieldT,DirT>::MaskPtr mask,
                            const SpatialOps::BCType bcType,
                            const SpatialOps::BCSide side,
                            const double value )
    : BoundaryConditionExpression<FieldT,DirT>(mask,bcType,side),
      value_(value)
    {
      this->set_gpu_runnable(true);
    }
  
  public:
  
    class Builder : public ExpressionBuilder
    {
    public:
      typedef typename BoundaryConditionExpression<FieldT,DirT>::MaskPtr  MaskPtr;
      typedef typename BoundaryConditionExpression<FieldT,DirT>::MaskType MaskType;
      /**
       *  @brief Build a ConstantBCExpression expression
       *  @param resultTag the tag for the value that this expression computes
       *  @param mask the point(s) to apply the BC on
       *  @param bcType DIRICHLET or NEUMANN
       *  @param side the side of the domain we are applying BCs MINUS_SIDE or PLUS_SIDE
       *  @param value the value to set for the BC
       */
      Builder( const Tag& resultTag,
               typename BoundaryConditionExpression<FieldT,DirT>::MaskPtr mask,
               const SpatialOps::BCType bcType,
               const SpatialOps::BCSide side,
               const double value )
      : ExpressionBuilder( resultTag ),
        mask_(mask), bcType_(bcType), side_(side), value_(value)
      {
        assert( !mask->points().empty() );  // why would we create a BC expression with no points?
      }
  
      Expr::ExpressionBase* build() const{ return new ConstantBCOpExpression<FieldT,DirT>(mask_,bcType_,side_,value_); }
  
    private:
      typename BoundaryConditionExpression<FieldT,DirT>::MaskPtr mask_;
      const SpatialOps::BCType bcType_;
      const SpatialOps::BCSide side_;
      const double value_;
    };
  
    ~ConstantBCOpExpression(){}
    void evaluate()
    {
      FieldT& result = this->value();
      switch( this->bcType_ ){
        case SpatialOps::DIRICHLET: (*this->dirichletOp_)( this->get_mask(), result, value_, this->side_ == SpatialOps::MINUS_SIDE ); break;
        case SpatialOps::NEUMANN  : (*this->neumannOp_  )( this->get_mask(), result, value_, this->side_ == SpatialOps::MINUS_SIDE ); break;
        case SpatialOps::UNSUPPORTED_BCTYPE: assert( false ); // shouldn't get here.
      }
    }
  
  };
  
  //--------------------------------------------------------------------
  
  } // namespace Expr
  
  #endif // BoundaryConditionExpression_Expr_h
  