module PhysicalFieldClass
   use ArrayClass
   use IOClass
   use LinearSolverClass
   implicit none

   type :: PhysicalField_
      character(len=:), allocatable, public:: name
      real(real64), allocatable, public:: scalar(:)
      real(real64), allocatable, public:: vector(:, :)
      real(real64), allocatable, public:: tensor(:, :, :)
      integer(int32) :: attribute = 0 ! node=1, element=2, gausspoint=3
      integer(int32) :: datastyle = 0 !scalar=1, vector=2, tensor=3
   contains

      procedure, pass :: importPhysicalFieldScalar
      procedure, pass :: importPhysicalFieldVector
      procedure, pass :: importPhysicalFieldTensor
      generic, public :: import => importPhysicalFieldScalar, &
         importPhysicalFieldVector, &
         importPhysicalFieldTensor
      procedure :: clear => clearPhysicalField
      procedure :: init => clearPhysicalField
      procedure :: remove => clearPhysicalField
      procedure :: msh => mshPhysicalField

   end type
contains

   subroutine importPhysicalFieldScalar(obj, scalar, name)
      class(PhysicalField_), intent(inout) :: obj
      real(real64), intent(in) :: scalar(:)
      character(*), intent(in) :: name

      obj%name = "untitled"
      obj%scalar = scalar
      obj%name = name

   end subroutine

   subroutine importPhysicalFieldVector(obj, vector, name)
      class(PhysicalField_), intent(inout) :: obj
      real(real64), intent(in) :: vector(:, :)
      character(*), intent(in) :: name

      obj%name = "untitled"
      obj%vector = vector
      obj%name = name

   end subroutine

   subroutine importPhysicalFieldtensor(obj, tensor, name)
      class(PhysicalField_), intent(inout) :: obj
      real(real64), intent(in) :: tensor(:, :, :)
      character(*), intent(in) :: name

      obj%name = "untitled"
      obj%tensor = tensor
      obj%name = name

   end subroutine

   subroutine clearPhysicalField(obj)
      class(PhysicalField_), intent(inout) :: obj

      if (allocated(obj%scalar)) then
         deallocate (obj%scalar)
      end if
      if (allocated(obj%vector)) then
         deallocate (obj%vector)
      end if
      if (allocated(obj%tensor)) then
         deallocate (obj%tensor)
      end if
      obj%name = "untitled"
   end subroutine
! ########################################################

! ########################################################
! export physical field by gmsh (.msh)
   subroutine mshPhysicalField(obj, name, caption)
      class(PhysicalField_), intent(inout) :: obj
      character(*), intent(in) :: name
      character(*), optional, intent(in) ::caption
      character(200) :: cap
      integer(int32) :: nb_scalar_points, nb_vector_points, nb_tensor_points, i, j
      real(real64), allocatable :: tensor(:, :), eigenVector(:, :)
      type(IO_) :: f
      ! doc: http://gmsh.info/doc/texinfo/gmsh.html#MSH-file-format-version-2-_0028Legacy_0029
      if (present(caption)) then
         cap = caption
      else
         cap = "untitled"
      end if
      if (.not. allocated(obj%scalar)) then
         nb_scalar_points = 0
      else
         nb_scalar_points = size(obj%scalar, 1)
      end if
      if (.not. allocated(obj%vector)) then
         nb_vector_points = 0
      else
         nb_vector_points = size(obj%vector, 1)
      end if

      if (.not. allocated(obj%tensor)) then
         nb_tensor_points = 0
      else
         nb_tensor_points = size(obj%tensor, 1)
      end if

      ! まずはスカラー
      if (nb_scalar_points /= 0 .and. obj%attribute == 1) then
         call f%open(name//"_node_scalar.msh")
         write (f%fh, '(a)') "$MeshFormat"
         write (f%fh, '(a)') "2.2 0 8"
         write (f%fh, '(a)') "$EndMeshFormat"
         write (f%fh, '(a)') "$NodeData"
         write (f%fh, '(a)') "1" ! one string tag
         write (f%fh, '(a)') cap ! the name of the view
         write (f%fh, '(a)') "1" ! one real tag
         write (f%fh, '(a)') "0.0" ! time value
         write (f%fh, '(a)') "3"   ! three integer tag
         write (f%fh, '(a)') "0"   ! the timestep (starts from 0)
         write (f%fh, '(a)') "1"! 1-component (scalar) field
         write (f%fh, '(a)') str(size(obj%scalar, 1))
         do i = 1, size(obj%scalar, 1)
            write (f%fh, '(a)') str(i)//" "//str(obj%scalar(i))
         end do
         write (f%fh, '(a)') "$EndNodeData"
         call f%close()
      end if
      ! 次にベクトル
      if (nb_vector_points /= 0 .and. obj%attribute == 1) then
         call f%open(name//"_node_vector.msh")
         write (f%fh, '(a)') "$MeshFormat"
         write (f%fh, '(a)') "2.2 0 8"
         write (f%fh, '(a)') "$EndMeshFormat"
         write (f%fh, '(a)') "$NodeData"
         write (f%fh, '(a)') "1" ! one string tag
         write (f%fh, '(a)') cap ! the name of the view
         write (f%fh, '(a)') "1" ! one real tag
         write (f%fh, '(a)') "0.0" ! time value
         write (f%fh, '(a)') "3"   ! three integer tag
         write (f%fh, '(a)') "0"   ! the timestep (starts from 0)
         write (f%fh, '(a)') str(size(obj%vector, 2))! n-component (vector) field
         write (f%fh, '(a)') str(size(obj%vector, 1))
         do i = 1, size(obj%vector, 1)
            write (f%fh, '(a)', advance="no") str(i)//" "
            do j = 1, size(obj%vector, 2) - 1
               write (f%fh, '(a)', advance="no") str(obj%vector(i, j))//" "
            end do
            j = size(obj%vector, 2)
            write (f%fh, '(a)', advance="yes") str(obj%vector(i, j))
         end do
         write (f%fh, '(a)') "$EndNodeData"
         call f%close()
      end if

      ! for element-data

      if (nb_scalar_points /= 0 .and. obj%attribute == 2) then
         call f%open(name//"_elem_scalar.msh")
         write (f%fh, '(a)') "$MeshFormat"
         write (f%fh, '(a)') "2.2 0 8"
         write (f%fh, '(a)') "$EndMeshFormat"
         write (f%fh, '(a)') "$ElementData"
         write (f%fh, '(a)') "1" ! one string tag
         write (f%fh, '(a)') cap ! the name of the view
         write (f%fh, '(a)') "1" ! one real tag
         write (f%fh, '(a)') "0.0" ! time value
         write (f%fh, '(a)') "3"   ! three integer tag
         write (f%fh, '(a)') "0"   ! the timestep (starts from 0)
         write (f%fh, '(a)') "1"! 1-component (scalar) field
         write (f%fh, '(a)') str(size(obj%scalar, 1))
         do i = 1, size(obj%scalar, 1)
            write (f%fh, '(a)') str(i)//" "//str(obj%scalar(i))
         end do
         write (f%fh, '(a)') "$EndElementData"
         call f%close()
      end if
      ! 次にベクトル
      if (nb_vector_points /= 0 .and. obj%attribute == 2) then
         call f%open(name//"_elem_vector.msh")
         write (f%fh, '(a)') "$MeshFormat"
         write (f%fh, '(a)') "2.2 0 8"
         write (f%fh, '(a)') "$EndMeshFormat"
         write (f%fh, '(a)') "$ElementData"
         write (f%fh, '(a)') "1" ! one string tag
         write (f%fh, '(a)') cap ! the name of the view
         write (f%fh, '(a)') "1" ! one real tag
         write (f%fh, '(a)') "0.0" ! time value
         write (f%fh, '(a)') "3"   ! three integer tag
         write (f%fh, '(a)') "0"   ! the timestep (starts from 0)
         write (f%fh, '(a)') str(size(obj%vector, 2))! n-component (vector) field
         write (f%fh, '(a)') str(size(obj%vector, 1))
         do i = 1, size(obj%vector, 1)
            write (f%fh, '(a)', advance="no") str(i)//" "
            do j = 1, size(obj%vector, 2) - 1
               write (f%fh, '(a)', advance="no") str(obj%vector(i, j))//" "
            end do
            j = size(obj%vector, 2)
            write (f%fh, '(a)', advance="yes") str(obj%vector(i, j))
         end do
         write (f%fh, '(a)') "$EndElementData"
         call f%close()
      end if

      ! テンソルの場合は実装中；
      if (nb_tensor_points /= 0 .and. obj%attribute == 3) then
         ! テンソルかつガウス点に定義（応力テンソルなど）
         ! 方針 >> 平均化して要素ごとに1つの主応力ベクトル図をかく
         i = size(obj%tensor, 2)
         j = size(obj%tensor, 3)
         if (i /= j) then
            print *, "mshFEMDomain >> ERROR >> size(obj%tensor,3)/=size(obj%tensor,2)"
            return
         end if
         allocate (tensor(i, j))
         if (i == 2) then
            call eigen_2d(tensor, eigenVector)
         elseif (i == 3) then
            allocate (eigenVector(3, 3))
            eigenVector(:, :) = eigen_3d(tensor)
         else
            print *, "ERROR :: mshPhysicalField for more dimension than 4-D>> not implemented yet."
            stop
         end if
      end if

   end subroutine
! ########################################################

end module
