Module class_SpaceDimension
!!=====================================================
!!=====================================================
!! A SpaceDimension is simply a CoordinateMap 
!! combined with a spatial index 1,2 or 3.
!! Recall...
!! Coordinate maps describe the sampling distribution 
!! associated with a coordinate.
!!
!! Notation: coordinate x is parameterized via s, 
!! i.e. x=x(s). The discrete arrays [x] [s] and [dxds] 
!! define the mapping. Particular maps are specified  
!! via their "label". Some maps require parameters, 
!! which are specified via the 2 vector params.
!!
!! Currently implemented maps: 
!!        unit, identity, cheby, periodic
!!=====================================================
!!=====================================================
 use class_CoordinateMap
 implicit none
 public  :: Map2SpaceDimension
 public  :: DeleteSpaceDimension
 public  :: DescribeSpaceDimension
 !public  :: MergeDimensions
 
 type SpaceDimension  
  integer                    :: DirIndex
  type( CoordinateMap )      :: map
 end type SpaceDimension
 
 interface assignment (=)
  module procedure EqualSpaceDimensions
 end interface
 
Contains      
  
  subroutine Map2SpaceDimension(map,id,spacedim)   
   use class_CoordinateMap
   implicit none  
   type( CoordinateMap ),intent(in)       :: map
   type( SpaceDimension ),intent(inout)   :: spacedim
   integer                                :: id
   spacedim%DirIndex=id !! i.e. 1,2 or 3   
   spacedim%map = map   
  end subroutine Map2SpaceDimension
  
  
  
  
  subroutine DeleteSpaceDimension(dim)
   use class_CoordinateMap
   type( SpaceDimension )  :: dim
   call DeleteCoordinateMap( dim%map )  !!DeleteCoordinateMap checks allocation first
  end subroutine DeleteSpaceDimension
   
   
   
   
  subroutine EqualSpaceDimensions(new,old)
   use class_CoordinateMap
   type( SpaceDimension ),intent(inout) :: new
   type( SpaceDimension ),intent(in)    :: old 
   new%map = NewCoordinateMap(old%map%npts,       &
                              old%map%label,      &
                              old%map%x,          &
                              old%map%s,          &
                              old%map%dxds,       &
                              old%map%d2xds2 )
   new%DirIndex = old%DirIndex
  end subroutine EqualSpaceDimensions


  

  
!  function MergeDimensions(xdim,ydim,zdim) result(AllSpaceDims)
!   type( SpaceDimension ),intent(in)   :: xdim
!   type( SpaceDimension ),intent(in)   :: ydim
!   type( SpaceDimension ),intent(in)   :: zdim
!   type( SpaceDimension )   :: AllSpaceDims(3)
!    !! type SpaceDimension has an overloaded assignment operator
!     AllSpaceDims(1) = xdim
!     AllSpaceDims(2) = ydim
!     AllSpaceDims(3) = zdim    
!  end function MergeDimensions
    
  
  subroutine DescribeSpaceDimension(dim,docfile)
   use class_CoordinateMap
   implicit none
   type( SpaceDimension ), intent(in):: dim
   character (len=80 )               :: coordlabel
   character (len=80 )               :: coordname
   character (len=80 )               :: docfile
   integer                           :: myid
   integer                           :: ierr
   include 'mpif.h'
   
   call mpi_comm_rank(MPI_COMM_WORLD,myid,ierr) 
   
   select case( dim%DirIndex )
    case(1)
     coordlabel='x  (or x1 or xsi), parameterized by s -----> x(s)'
     coordname=' x '
    case(2)
     coordlabel='y  (or y1 or eta), parameterized by s -----> y(s)'
     coordname=' y '
    case(3)
     coordlabel='z  (or z1 or zeta), parameterized by s -----> z(s)'
     coordname=' z '
   end select
   
   if(myid==0) then
    open(1,file=docfile,position='append')
    write(1,*) ' '
    write(1,*) '============================================================'
    write(1,*) 'SpaceDimension  ',dim%DirIndex
    write(1,*)  trim(coordlabel)   
    call DescribeCoordinateMap(dim%map,coordname,docfile,myid)
    write(1,*) '  '
    close(1)
   endif
  end subroutine DescribeSpaceDimension
    
end Module class_SpaceDimension !============================================