//#####################################################################
// Copyright 2017, Ounan Ding, Lin Huang, Craig Schroeder.
// This file is part of PhysBAM whose distribution is governed by the license contained in the accompanying file PHYSBAM_COPYRIGHT.txt.
//#####################################################################
// Class FACE_ITERATOR_THREADED
//##################################################################### 
#include <Core/Grids/FACE_ITERATOR_THREADED.h>
#ifdef USE_OPENMP
#include <omp.h>
#endif
using namespace PhysBAM;
template<class TV> FACE_ITERATOR_THREADED<TV>::
FACE_ITERATOR_THREADED(const GRID<TV>& grid_input,const RANGE<IV>& domain,bool whole_region)
    :grid(grid_input)
{
#ifdef USE_OPENMP
    int threads=omp_get_num_threads();
    int tid=omp_get_thread_num();
#else
    int threads=1;
    int tid=0;
#endif
    
    Reset_Regions();
    int x_cells=domain.max_corner.x-domain.min_corner.x;

    int first=TV::m*x_cells*tid/threads;
    int last=TV::m*x_cells*(tid+1)/threads;
    for(int a=0;a<TV::m;a++){
        if(first<x_cells){
            RANGE<IV> axis_domain(domain);
            axis_domain.min_corner(0)=domain.min_corner(0)+first;
            if(whole_region) axis_domain.max_corner(a)++;
            else if(!first) axis_domain.min_corner(a)++;
            if(last<x_cells) axis_domain.max_corner(0)=domain.min_corner(0)+last;
            axes[number_of_regions]=a;
#pragma omp critical
            if(axis_domain.min_corner(0)<axis_domain.max_corner(0))
                Add_Region(axis_domain);}
        last-=x_cells;
        if(last<0) break;
        first-=x_cells;
        if(first<0) first=0;}
    axis=axes[0];
    Reset();
}
//#####################################################################
// Constructor
//#####################################################################
// axis_input==0 means iterate through faces in all dimensions
template<class TV> FACE_ITERATOR_THREADED<TV>::
FACE_ITERATOR_THREADED(const GRID<TV>& grid_input,int number_of_ghost_cells,T_REGION region_type)
    :FACE_ITERATOR_THREADED(grid_input,grid_input.Cell_Indices(number_of_ghost_cells),region_type==WHOLE_REGION)
{
    PHYSBAM_ASSERT(region_type==WHOLE_REGION || region_type==INTERIOR_REGION);
}
//#####################################################################
namespace PhysBAM{
template class FACE_ITERATOR_THREADED<VECTOR<float,1> >;
template class FACE_ITERATOR_THREADED<VECTOR<float,2> >;
template class FACE_ITERATOR_THREADED<VECTOR<float,3> >;
template class FACE_ITERATOR_THREADED<VECTOR<double,1> >;
template class FACE_ITERATOR_THREADED<VECTOR<double,2> >;
template class FACE_ITERATOR_THREADED<VECTOR<double,3> >;
}
