text: #include <vector>
  #include <iostream>
  #include <sstream>
  #include <iomanip>
  #include <limits>
  #include <cmath>
  
  using std::cout;
  using std::endl;
  
  #include <spatialops/structured/FVStaggered.h>
  
  #include <spatialops/structured/Grid.h>
  #include <test/TestHelper.h>
  
  #include <spatialops/structured/FieldHelper.h>
  using namespace SpatialOps;
  
  template<typename OpT, typename DirT>
  bool test_bc_conversion_to_mask
                         ( const OperatorDatabase& opDB,
                           const IntVec&dim,
                           const IntVec& bcFlag,
                           const IntVec ijk,
                           const BCSide side,
                           const double bcVal )
  {
    using namespace SpatialOps;
    
    BCSide conv_side = NO_SIDE;
    if(side == MINUS_SIDE) conv_side = PLUS_SIDE;
    else conv_side = MINUS_SIDE;
  
    typedef typename OpT::SrcFieldType  SrcFieldT;
    typedef typename OpT::DestFieldType DestFieldT;
    typedef NeboBoundaryConditionBuilder<OpT>  BCOpT;
  
    const   OpT&   op = *opDB.retrieve_operator<  OpT>();
    const BCOpT& bcOp = *opDB.retrieve_operator<BCOpT>();
  
    const GhostData  fghost(1);
    const GhostData dfghost(1);
    const BoundaryCellInfo  fbc = BoundaryCellInfo::build< SrcFieldT>( bcFlag,bcFlag );
    const BoundaryCellInfo dfbc = BoundaryCellInfo::build<DestFieldT>( bcFlag,bcFlag );
  
    SpatFldPtr<SrcFieldT > f  = SpatialFieldStore::get_from_window<SrcFieldT >( get_window_with_ghost(dim, fghost, fbc),  fbc,  fghost );
    SpatFldPtr<DestFieldT> df = SpatialFieldStore::get_from_window<DestFieldT>( get_window_with_ghost(dim,dfghost,dfbc), dfbc, dfghost );
  
    int icnt=0;
    for( typename SrcFieldT::iterator ifld=f->begin(); ifld!=f->end(); ++ifld,++icnt ) *ifld = icnt;
    *df <<= 0.0;
    IntVec shift;
    if(conv_side == PLUS_SIDE || conv_side == NO_SIDE) 
        shift = BCOpT::Shift::MinusPoint::int_vec();
    else 
        shift = BCOpT::Shift::PlusPoint::int_vec();
  
    IntVec mask_point(ijk[0] - shift[0],
                      ijk[1] - shift[1],
                      ijk[2] - shift[2]);
    std::vector<IntVec> maskPoints(1, mask_point);
    const SpatialMask<SrcFieldT> mask( *f, maskPoints );
  
    // apply BC
    bcOp(convert<DestFieldT, DirT>(mask, conv_side).to_mask(*df), *f, bcVal, conv_side != MINUS_SIDE);
  
    // calculate the dest field
    op.apply_to_field( *f, *df );
  
    // verify that the BC was set properly.
    const int ix = df->window_without_ghost().flat_index(ijk);
  
    const double abserr = std::abs( (*df)[ix] - bcVal );
    const double relerr = abserr/std::abs(bcVal);
  
    const bool isOkay = (abserr<2.0e-12 && relerr<2.0e-12);
  
  
    if( !isOkay )
      cout << "expected: " << bcVal << endl
           << "found: " << (*df)[ix] << endl
           << "abserr: " << abserr << ", relerr: " << relerr << endl
           << endl;
    
    // apply BC
    bcOp(convert<DestFieldT, DirT>(mask, conv_side).to_mask(*df), *f, 1.0 + bcVal, side == MINUS_SIDE);
    
    // calculate the dest field
    op.apply_to_field( *f, *df );
  
    const bool isOkay2 = (std::abs(bcVal + 1.0 - (*df)[ix]) < 1e-10);
  
    if( !isOkay2 )
      cout << "failed to use converted mask from other type" <<endl
             << "expected:" << 1.0 + bcVal <<endl
             << "found: " << (*df)[ix] <<endl;
    
    *df <<= cond(convert<DestFieldT, DirT>(mask, conv_side).to_mask(*df), -bcVal)
                (*df);
  
  
    const bool isCondOkay = (std::abs((*df)[ix] + bcVal) < 1e-10);
    if(!isCondOkay)
        cout << "failed to apply converted mask in cond!" <<endl
             << "expected:" << -bcVal <<endl
             << "found: " << (*df)[ix] <<endl;
  
  
    *df <<= cond(convert<DestFieldT, DirT>(mask, PLUS_SIDE).to_mask(*df), -567.00)
                (0);
    
    *df <<= cond(convert<DestFieldT, DirT>(mask, MINUS_SIDE).to_mask(*df), -789.00)
                (*df);
  
    *df <<= cond(convert<DestFieldT, DirT>(mask, MINUS_SIDE, PLUS_SIDE).to_mask(*df), 0.00)
                (*df);
  
    const bool isMutipleFacesOkay = nebo_norm(*df) < 1e-10;
    if(!isMutipleFacesOkay)
        cout << "failed to convert mask with multiple faces!" <<endl
             << "expected:" << 0.00 <<endl
             << "found: " << nebo_norm(*df) <<endl;
    
    masked_assign(convert<DestFieldT, DirT>(mask, conv_side), *df, bcVal);
    const bool isMaskedAssignOkay = (std::abs((*df)[ix] - bcVal) < 1e-10);
    if(!isMaskedAssignOkay)
        cout << "failed to use converted mask in masked_assign!" <<endl
             << "expected:" << bcVal <<endl
             << "found: " << (*df)[ix] <<endl;
    
    return isOkay && isCondOkay && isMutipleFacesOkay && isOkay2;
  }
  //--------------------------------------------------------------------
  template<typename OpT, typename DirT>
  bool test_bc_conversion( const OperatorDatabase& opDB,
                           const IntVec&dim,
                           const IntVec& bcFlag,
                           const IntVec ijk,
                           const BCSide side,
                           const double bcVal )
  {
    using namespace SpatialOps;
    
    BCSide conv_side = NO_SIDE;
    if(side == MINUS_SIDE) conv_side = PLUS_SIDE;
    else conv_side = MINUS_SIDE;
  
    typedef typename OpT::SrcFieldType  SrcFieldT;
    typedef typename OpT::DestFieldType DestFieldT;
    typedef NeboBoundaryConditionBuilder<OpT>  BCOpT;
  
    const   OpT&   op = *opDB.retrieve_operator<  OpT>();
    const BCOpT& bcOp = *opDB.retrieve_operator<BCOpT>();
  
    const GhostData  fghost(1);
    const GhostData dfghost(1);
    const BoundaryCellInfo  fbc = BoundaryCellInfo::build< SrcFieldT>( bcFlag, bcFlag );
    const BoundaryCellInfo dfbc = BoundaryCellInfo::build<DestFieldT>( bcFlag, bcFlag );
  
    SpatFldPtr<SrcFieldT > f  = SpatialFieldStore::get_from_window<SrcFieldT >( get_window_with_ghost(dim, fghost, fbc),  fbc,  fghost );
    SpatFldPtr<DestFieldT> df = SpatialFieldStore::get_from_window<DestFieldT>( get_window_with_ghost(dim,dfghost,dfbc), dfbc, dfghost );
  
    int icnt=0;
    for( typename SrcFieldT::iterator ifld=f->begin(); ifld!=f->end(); ++ifld,++icnt ) *ifld = icnt;
    *df <<= 0.0;
    IntVec shift;
    if(conv_side == PLUS_SIDE || conv_side == NO_SIDE) 
        shift = BCOpT::Shift::MinusPoint::int_vec();
    else 
        shift = BCOpT::Shift::PlusPoint::int_vec();
  
    IntVec mask_point(ijk[0] - shift[0],
                      ijk[1] - shift[1],
                      ijk[2] - shift[2]);
    std::vector<IntVec> maskPoints(1, mask_point);
    const SpatialMask<SrcFieldT> mask( *f, maskPoints );
  
    // apply BC
    bcOp(convert<DestFieldT, DirT>(mask, conv_side), *f, bcVal);
  
    // calculate the dest field
    op.apply_to_field( *f, *df );
  
    // verify that the BC was set properly.
    const int ix = df->window_without_ghost().flat_index(ijk);
  
    const double abserr = std::abs( (*df)[ix] - bcVal );
    const double relerr = abserr/std::abs(bcVal);
  
    const bool isOkay = (abserr<2.0e-12 && relerr<2.0e-12);
  
  
    if( !isOkay )
      cout << "expected: " << bcVal << endl
           << "found: " << (*df)[ix] << endl
           << "abserr: " << abserr << ", relerr: " << relerr << endl
           << endl;
    
    // apply BC
    bcOp(convert<DestFieldT, DirT>(mask, conv_side), *f, 1.0 + bcVal, side == MINUS_SIDE);
    
    // calculate the dest field
    op.apply_to_field( *f, *df );
  
    const bool isOkay2 = std::abs(bcVal + 1.0 - (*df)[ix]) < 1e-10;
  
    if( !isOkay2 )
      cout << "failed to use converted mask from other type" <<endl
             << "expected:" << 1.0 + bcVal <<endl
             << "found: " << (*df)[ix] <<endl;
    
    *df <<= cond(convert<DestFieldT, DirT>(mask, conv_side).mask(), -bcVal)
                (*df);
  
  
    const bool isCondOkay = std::abs((*df)[ix] + bcVal) < 1e-10;
    if(!isCondOkay)
        cout << "failed to apply converted mask in cond!" <<endl
             << "expected:" << -bcVal <<endl
             << "found: " << (*df)[ix] <<endl;
  
  
    *df <<= cond(convert<DestFieldT, DirT>(mask, PLUS_SIDE).mask(), -567.00)
                (0);
    
    *df <<= cond(convert<DestFieldT, DirT>(mask, MINUS_SIDE).mask(), -789.00)
                (*df);
  
    *df <<= cond(convert<DestFieldT, DirT>(mask, MINUS_SIDE, PLUS_SIDE).mask(), 0.00)
                (*df);
  
    const bool isMutipleFacesOkay = nebo_norm(*df) < 1e-10;
    if(!isMutipleFacesOkay)
        cout << "failed to convert mask with multiple faces!" <<endl
             << "expected:" << 0.00 <<endl
             << "found: " << nebo_norm(*df) <<endl;
    
    masked_assign(convert<DestFieldT, DirT>(mask, conv_side), *df, bcVal);
    const bool isMaskedAssignOkay = std::abs((*df)[ix] - bcVal) < 1e-10;
    if(!isMaskedAssignOkay)
        cout << "failed to use converted mask in masked_assign!" <<endl
             << "expected:" << bcVal <<endl
             << "found: " << (*df)[ix] <<endl;
    
    return isOkay && isCondOkay && isMutipleFacesOkay && isOkay2;
  }
  template<typename OpT>
  bool test_bc_helper( const OperatorDatabase& opDB,
                       const IntVec&dim,
                       const IntVec& bcFlag,
                       const IntVec ijk,
                       const BCSide side,
                       const double bcVal )
  {
    using namespace SpatialOps;
  
    typedef typename OpT::SrcFieldType  SrcFieldT;
    typedef typename OpT::DestFieldType DestFieldT;
    typedef NeboBoundaryConditionBuilder<OpT>  BCOpT;
  
    const   OpT&   op = *opDB.retrieve_operator<  OpT>();
    const BCOpT& bcOp = *opDB.retrieve_operator<BCOpT>();
  
    const GhostData  fghost(1);
    const GhostData dfghost(1);
    const BoundaryCellInfo  fbc = BoundaryCellInfo::build< SrcFieldT>( bcFlag,bcFlag );
    const BoundaryCellInfo dfbc = BoundaryCellInfo::build<DestFieldT>( bcFlag,bcFlag );
  
    SpatFldPtr<SrcFieldT > f  = SpatialFieldStore::get_from_window<SrcFieldT >( get_window_with_ghost(dim, fghost, fbc),  fbc,  fghost );
    SpatFldPtr<DestFieldT> df = SpatialFieldStore::get_from_window<DestFieldT>( get_window_with_ghost(dim,dfghost,dfbc), dfbc, dfghost );
  
    int icnt=0;
    for( typename SrcFieldT::iterator ifld=f->begin(); ifld!=f->end(); ++ifld,++icnt ) *ifld = icnt;
    *df <<= 0.0;
  
    std::vector<IntVec> maskPoints(1,ijk);
    const SpatialMask<DestFieldT> mask( *df, maskPoints );
  
    // apply BC
    bcOp( mask, *f, bcVal, side==MINUS_SIDE );
  
    // calculate the dest field
    op.apply_to_field( *f, *df );
  
    // verify that the BC was set properly.
    const int ix = df->window_without_ghost().flat_index(ijk);
  
    const double abserr = std::abs( (*df)[ix] - bcVal );
    const double relerr = abserr/std::abs(bcVal);
  
    const bool isOkay = (abserr<2.0e-12 && relerr<2.0e-12);
    if( !isOkay )
      cout << "expected: " << bcVal << endl
           << "found: " << (*df)[ix] << endl
           << "abserr: " << abserr << ", relerr: " << relerr << endl
           << endl;
    return isOkay;
  }
  
  //--------------------------------------------------------------------
  
  template<typename OpT>
  bool
  test_bc_loop( const int dir,
                const OperatorDatabase& opDB,
                const std::string opName,
                const IntVec& dim,
                const IntVec& bcFlag,
                const int bcFaceIndex,
                const BCSide side,
                const double bcVal )
  {
    TestHelper status(false);
  
    typedef typename GetOperatorDir<OpT>::result DirT;
  
    if( dir == 0 ){ //X direction
      TestHelper tmp(false);
      const int i=bcFaceIndex;
      for( int j=0; j<dim[1]; ++j ){
        for( int k=0; k<dim[2]; ++k ){
          tmp( test_bc_helper<OpT>( opDB, dim, bcFlag, IntVec(i,j,k), side, bcVal ) );
          tmp( test_bc_conversion<OpT, DirT>( opDB, dim, bcFlag, IntVec(i,j,k), side, bcVal ) );
          tmp( test_bc_conversion_to_mask<OpT, DirT>( opDB, dim, bcFlag, IntVec(i,j,k), side, bcVal ) );
        }
      }
      status( tmp.ok(), "X BC " + opName );
    }
    else if( dir == 1 ){ //Y direction
      TestHelper tmp(false);
      const int j=bcFaceIndex;
      for( int i=0; i<dim[0]; ++i ){
        for( int k=0; k<dim[2]; ++k ){
          tmp( test_bc_helper<OpT>( opDB, dim, bcFlag, IntVec(i,j,k), side, bcVal ) );
          tmp( test_bc_conversion<OpT, DirT>( opDB, dim, bcFlag, IntVec(i,j,k), side, bcVal ) );
          tmp( test_bc_conversion_to_mask<OpT, DirT>( opDB, dim, bcFlag, IntVec(i,j,k), side, bcVal ) );
        }
      }
      status( tmp.ok(), "Y BC " + opName );
    }
    else if( dir == 2 ){ //Z direction
      TestHelper tmp(false);
      const int k=bcFaceIndex;
      for( int i=0; i<dim[0]; ++i ){
        for( int j=0; j<dim[1]; ++j ){
          tmp( test_bc_helper<OpT>( opDB, dim, bcFlag, IntVec(i,j,k), side, bcVal ) );
          tmp( test_bc_conversion<OpT, DirT>( opDB, dim, bcFlag, IntVec(i,j,k), side, bcVal ) );
          tmp( test_bc_conversion_to_mask<OpT, DirT>( opDB, dim, bcFlag, IntVec(i,j,k), side, bcVal ) );
        }
      }
      status( tmp.ok(), "Z BC " + opName );
    }
    else{
      std::cout << "ERROR - invalid direction detected on operator!" << std::endl;
      status(false);
    }
  
    return status.ok();
  }
  
  //--------------------------------------------------------------------
  
  template< typename VolT >
  bool test_bc( OperatorDatabase& opDB,
                const Grid& g,
                const IntVec& bcFlag )
  {
    using namespace SpatialOps;
  
    typedef BasicOpTypes<VolT> Ops;
  
    const IntVec& dim = g.extent();
  
    TestHelper status(true);
  
    if( dim[0]>1 ){
      // X BCs - Left side
      int i=0;
      status( test_bc_loop<typename Ops::GradX     >( 0, opDB, "Grad   Vol->SurfX", dim, bcFlag, i,   MINUS_SIDE, 1.2345 ), "-X Grad   Vol->SurfX" );
      status( test_bc_loop<typename Ops::DivX      >( 0, opDB, "Div    SurfX->Vol", dim, bcFlag, i,   MINUS_SIDE, 6.2345 ), "-X Div    SurfX->Vol" );
      status( test_bc_loop<typename Ops::InterpC2FX>( 0, opDB, "Interp Vol->SurfX", dim, bcFlag, i,   MINUS_SIDE, 123.45 ), "-X Interp Vol->SurfX" );
  
      status( test_bc_loop<typename OperatorTypeBuilder<GradientX,   VolT,VolT>::type >( 0, opDB,"Grad   Vol->Vol", dim, bcFlag, i, NO_SIDE, 1.2345 ), "-X Grad   Vol->Vol" );
      status( test_bc_loop<typename OperatorTypeBuilder<InterpolantX,VolT,VolT>::type >( 0, opDB,"Interp Vol->Vol", dim, bcFlag, i, NO_SIDE, 1.2345 ), "-X Interp Vol->Vol" );
      cout << endl;
  
      // X BCs - Right side
      i=dim[0];
      status( test_bc_loop<typename Ops::GradX     >( 0, opDB, "Grad   Vol->SurfX", dim, bcFlag, i,   PLUS_SIDE, 1.2345 ), "+X Grad   Vol->SurfX" );
      status( test_bc_loop<typename Ops::DivX      >( 0, opDB, "Div    SurfX->Vol", dim, bcFlag, i,   PLUS_SIDE, 6.789 ),  "+X Div    SurfX->Vol" );
      status( test_bc_loop<typename Ops::InterpC2FX>( 0, opDB, "Interp Vol->SurfX", dim, bcFlag, i,   PLUS_SIDE, 123.45 ), "+X Interp Vol->SurfX" );
  
      status( test_bc_loop<typename OperatorTypeBuilder<GradientX,   VolT,VolT>::type >( 0, opDB,"Grad   Vol->Vol", dim, bcFlag, i-1, NO_SIDE, 1.2345 ), "+X Grad   Vol->Vol" );
      status( test_bc_loop<typename OperatorTypeBuilder<InterpolantX,VolT,VolT>::type >( 0, opDB,"Interp Vol->Vol", dim, bcFlag, i-1, NO_SIDE, 1.2345 ), "+X Interp Vol->Vol" );
      cout << endl;
    }
  
    if( dim[1]>1 ){
      // Y BCs - Left side
      int j=0;
      status( test_bc_loop<typename Ops::GradY     >( 1, opDB, "Grad   Vol->SurfY", dim, bcFlag, j,   MINUS_SIDE, 1.23456 ), "-Y Grad   Vol->SurfY" );
      status( test_bc_loop<typename Ops::DivY      >( 1, opDB, "Div    SurfY->Vol", dim, bcFlag, j,   MINUS_SIDE, 1.23456 ), "-Y Div    SurfY->Vol" );
      status( test_bc_loop<typename Ops::InterpC2FY>( 1, opDB, "Interp Vol->SurfY", dim, bcFlag, j,   MINUS_SIDE, 123.456 ), "-Y Interp Vol->SurfY" );
  
      status( test_bc_loop<typename OperatorTypeBuilder<GradientY,   VolT,VolT>::type >( 1, opDB,"Grad   Vol->Vol", dim, bcFlag, j, NO_SIDE, 1.2345 ), "-Y Grad   Vol->Vol" );
      status( test_bc_loop<typename OperatorTypeBuilder<InterpolantY,VolT,VolT>::type >( 1, opDB,"Interp Vol->Vol", dim, bcFlag, j, NO_SIDE, 1.2345 ), "-Y Interp Vol->Vol" );
      cout << endl;
  
      // Y BCs - Right side
      j=dim[1];
      status( test_bc_loop<typename Ops::GradY     >( 1, opDB, "Grad   Vol->SurfY", dim, bcFlag, j,   PLUS_SIDE, 6.54321 ), "+Y Grad   Vol->SurfY" );
      status( test_bc_loop<typename Ops::DivY      >( 1, opDB, "Div    SurfY->Vol", dim, bcFlag, j,   PLUS_SIDE, 1.23456 ), "+Y Div    SurfY->Vol" );
      status( test_bc_loop<typename Ops::InterpC2FY>( 1, opDB, "Interp Vol->SurfY", dim, bcFlag, j,   PLUS_SIDE, 123.456 ), "+Y Interp Vol->SurfY" );
  
      status( test_bc_loop<typename OperatorTypeBuilder<GradientY,   VolT,VolT>::type >( 1, opDB,"Grad   Vol->Vol", dim, bcFlag, j-1, NO_SIDE, 1.2345 ), "+Y Grad   Vol->Vol" );
      status( test_bc_loop<typename OperatorTypeBuilder<InterpolantY,VolT,VolT>::type >( 1, opDB,"Interp Vol->Vol", dim, bcFlag, j-1, NO_SIDE, 1.2345 ), "+Y Interp Vol->Vol" );
      cout << endl;
    }
  
    if( dim[2]>1 ){
      // Z BCs - Left side
      int k=0;
      status( test_bc_loop<typename Ops::GradZ     >( 2, opDB, "Grad   Vol->SurfZ", dim, bcFlag, k,   MINUS_SIDE, 1.2345  ), "-Z Grad   Vol->SurfZ" );
      status( test_bc_loop<typename Ops::DivZ      >( 2, opDB, "Div    SurfZ->Vol", dim, bcFlag, k,   MINUS_SIDE, 2.34567 ), "-Z Div    SurfZ->Vol" );
      status( test_bc_loop<typename Ops::InterpC2FZ>( 2, opDB, "Interp Vol->SurfZ", dim, bcFlag, k,   MINUS_SIDE, 123.456 ), "-Z Interp Vol->SurfZ" );
  
      status( test_bc_loop<typename OperatorTypeBuilder<GradientZ,   VolT,VolT>::type >( 2, opDB,"Grad   Vol->Vol", dim, bcFlag, k, NO_SIDE, 1.2345 ), "-Z Grad   Vol->Vol" );
      status( test_bc_loop<typename OperatorTypeBuilder<InterpolantZ,VolT,VolT>::type >( 2, opDB,"Interp Vol->Vol", dim, bcFlag, k, NO_SIDE, 1.2345 ), "-Z Interp Vol->Vol" );
      cout << endl;
  
      // Z BCs - Right side
      k=dim[2];
      status( test_bc_loop<typename Ops::GradZ     >( 2, opDB, "Grad   Vol->SurfZ", dim, bcFlag, k,   PLUS_SIDE, 1.2345  ), "+Z Grad   Vol->SurfZ" );
      status( test_bc_loop<typename Ops::DivZ      >( 2, opDB, "Div    SurfZ->Vol", dim, bcFlag, k,   PLUS_SIDE, 7.89012 ), "+Z Div    SurfZ->Vol" );
      status( test_bc_loop<typename Ops::InterpC2FZ>( 2, opDB, "Interp Vol->SurfZ", dim, bcFlag, k,   PLUS_SIDE, 123.456 ), "+Z Interp Vol->SurfZ" );
  
      status( test_bc_loop<typename OperatorTypeBuilder<GradientZ,   VolT,VolT>::type >( 2, opDB,"Grad   Vol->Vol", dim, bcFlag, k-1, NO_SIDE, 1.2345 ), "+Z Grad   Vol->Vol" );
      status( test_bc_loop<typename OperatorTypeBuilder<InterpolantZ,VolT,VolT>::type >( 2, opDB,"Interp Vol->Vol", dim, bcFlag, k-1, NO_SIDE, 12.345 ), "+Z Interp Vol->Vol" );
      cout << endl;
    }
  
    return status.ok();
  }
  
  //--------------------------------------------------------------------
  
  bool test_driver( const IntVec& dim )
  {
    const DoubleVec length(1,1,1);
    const IntVec bcFlag(true,true,true);
  
    OperatorDatabase opDB;
    build_stencils( dim[0], dim[1], dim[2], length[0], length[1], length[2], opDB );
  
    const Grid grid( dim, length );
  
    TestHelper status(true);
  
    status( test_bc<SVolField>( opDB, grid, bcFlag ), "SVolField Ops" );
  
    if( dim[0]>1 ) status( test_bc<XVolField>( opDB, grid, bcFlag ), "XVolField Ops" );
    if( dim[1]>1 ) status( test_bc<YVolField>( opDB, grid, bcFlag ), "YVolField Ops" );
    if( dim[2]>1 ) status( test_bc<ZVolField>( opDB, grid, bcFlag ), "ZVolField Ops" );
  
    return status.ok();
  }
  
  //--------------------------------------------------------------------
  
  int main()
  {
    TestHelper status(true);
  
    try{
      status( test_driver( IntVec(10,1 ,1 ) ), "Mesh: (10,1,1)\n" );
      status( test_driver( IntVec(1 ,10,1 ) ), "Mesh: (1,10,1)\n" );
      status( test_driver( IntVec(1 ,1 ,10) ), "Mesh: (1,1,10)\n" );
  
      status( test_driver( IntVec(10,10,1 ) ), "Mesh: (10,10,1)\n" );
      status( test_driver( IntVec(10,1 ,10) ), "Mesh: (10,1,10)\n" );
      status( test_driver( IntVec(1 ,10,10) ), "Mesh: (1,10,10)\n" );
  
      status( test_driver( IntVec(10,10,10) ), "Mesh: (10,10,10)\n" );
  
      cout << endl << "----------" << endl
          << "BC Op Test: ";
      if( status.ok() ){
        cout << "PASS" << endl << "----------" << endl;
        return 0;
      }
    }
    catch( std::exception& err ){
      cout << err.what() << std::endl;
    }
    cout << "FAIL" << endl << "----------" << endl;
    return -1;
  }
  