module GeometryClass
   use, intrinsic :: iso_fortran_env
   use MathClass
   use ArrayClass
   implicit none

   type Point_
      real(real64), allocatable :: coord(:)
      character*30 :: name
   contains
      procedure :: Init => InitPoint
      procedure :: set => setPoint
      procedure :: show => showPoint
   end type

   type Line_
      real(real64), allocatable :: coord(:, :)
   contains
      procedure :: Init => InitLine
      procedure :: setNode => SetNodeLine
      procedure :: import => importLine
      procedure :: show => showLine
   end type

   type Circle_
      real(real64) :: radius
      real(real64), allocatable :: center(:)
   contains
      procedure :: Init => InitCircle
      procedure :: SetCenter => InitSetCenterCircle
      procedure :: SetRadius => InitSetRadiusCircle
      procedure :: getArea => getAreaCircle
      procedure :: show => showCircle
   end type

   type Sphere_
      real(real64) :: radius
      real(real64) :: center(3)
   contains
      procedure :: Init => InitSphere
      procedure :: SetCenter => InitSetCenterSphere
      procedure :: SetRadius => InitSetRadiusSphere
      procedure :: show => showSphere
   end type

   type Triangle_
      real(real64), allocatable :: NodCoord(:, :)
      real(real64), allocatable :: OuterNormal(:)
      real(real64), allocatable :: Center(:)
   contains
      procedure :: Init => InitTriangle
      procedure :: setNode => setNodeTriangle
      procedure :: import => importTriangle
      procedure :: getCircle => getCircleTriangle
      procedure :: getArea => getAreaTriangle
      procedure :: show => showTriangle
      procedure :: GetOuterNormal => GetOuterNormalTriangle
   end type

   type Rectangle_
      real(real64), allocatable :: NodCoord(:, :)
   contains
      procedure :: Init => InitRectangle
      procedure :: create => createRectangle
      procedure :: move => moveRectangle
      procedure :: setNode => setNodeRectangle
      procedure :: import => importRectangle
      procedure :: getCircle => getCircleRectangle
      procedure :: getArea => getAreaRectangle
      procedure :: show => showRectangle
      procedure :: contact => contactRectangle
   end type

   type Tetragon_
      real(real64), allocatable :: NodCoord(:, :)
   end type

   type Tetrahedron_
      real(real64) :: NodCoord(4, 3)
      real(real64) :: radius
      real(real64) :: center(3)
   contains
      procedure :: Init => InitTetrahedron
      procedure :: getCircle => getCircleTetrahedron
   end type

   type Octahedron_
      real(real64), allocatable :: NodCoord(:, :)
   end type
contains

!#########################################################
   subroutine InitPoint(obj, dim)
      class(Point_), intent(inout)::obj
      integer(int32), optional, intent(in)::dim

      ! default == 3
      if (allocated(obj%coord)) then
         deallocate (obj%coord)
      end if
      allocate (obj%coord(input(default=3, option=dim)))
      obj%coord(:) = 0.0d0

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

!#########################################################
   subroutine showPoint(obj, Name)
      class(Point_), intent(in)::obj
      character(*), optional, intent(in)::Name

      if (present(Name)) then
         open (10, file=Name)
         write (10, *) obj%coord(:)
         close (10)
      else
         print *, "Point> x(:) = ", obj%coord(:)
      end if
   end subroutine
!#########################################################

!#########################################################
   subroutine setPoint(obj, x, y, z, xvec)
      class(Point_), intent(inout)::obj
      real(real64), optional, intent(in)::x, y, z
      real(real64), optional, intent(in)::xvec(:)
      integer(int32) :: n

      n = size(obj%coord)
      if (present(xvec)) then
         if (n /= size(xvec)) then
            print *, "ERROR :: setPoint :: n/=size(xvec)"
            return
         end if
         obj%coord(:) = xvec(:)
         return
      end if

      if (present(x)) then
         obj%coord(1) = x

      end if

      if (present(y)) then
         obj%coord(2) = y

      end if

      if (size(obj%coord) <= 2) then
         return
      end if

      if (present(z)) then
         obj%coord(3) = z
      end if

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

!#########################################################
   subroutine InitLine(obj, dim)
      class(Line_), intent(inout)::obj
      integer(int32), optional, intent(inout)::dim

      ! default = 3D
      if (allocated(obj%coord)) then
         deallocate (obj%coord)
      end if
      allocate (obj%coord(2, input(default=3, option=dim)))
      obj%coord(:, :) = 0.0d0

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

!#########################################################
   subroutine SetNodeLine(obj, point, position)
      class(Line_), intent(inout)::obj
      class(Point_), intent(in)::Point
      integer(int32), intent(in)::position

      if (position <= 0 .or. position >= 3) then
         print *, "ERROR :: Line%SetNode => (position <=0 .or. position >= 3)", position
         return
      end if

      if (size(obj%coord, 2) /= size(point%coord)) then
         print *, "ERROR :: Line%SetNode => (size(obj%coord,2)/=size(point%coord) )"
         return
      else
         obj%coord(position, :) = point%coord(:)
      end if

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

!#########################################################
   subroutine importLine(obj, NodCoord)
      class(Line_), intent(inout)::obj
      integer(int32), intent(in)::NodCoord(:, :)

      if (size(obj%coord, 2) /= size(NodCoord, 2)) then
         print *, "ERROR :: importLine :: size(obj%NodCoord,2)/=size(NodCoord,2)"
         return
      end if
      obj%coord(1:2, :) = NodCoord(1:2, :)

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

!#########################################################
   subroutine showLine(obj, Name)
      class(Line_), intent(in)::obj
      character(*), optional, intent(in)::Name

      if (present(Name)) then
         open (10, file=Name)
         write (10, *) obj%coord(1, :)
         write (10, *) obj%coord(2, :)
         close (10)
      else
         print *, "Line> x1(:) = ", obj%coord(1, :)
         print *, "Line> x2(:) = ", obj%coord(2, :)
      end if
   end subroutine
!#########################################################

!#########################################################
   subroutine InitCircle(obj, dim)
      class(Circle_), intent(inout)::obj
      integer(int32), optional, intent(inout)::dim

      ! default = 3D
      if (allocated(obj%center)) then
         deallocate (obj%center)
      end if
      allocate (obj%center(input(default=3, option=dim)))
      obj%center(:) = 0.0d0
      obj%radius = 1.0d0

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

!#########################################################
   subroutine InitSetCenterCircle(obj, point)
      class(Circle_), intent(inout)::obj
      class(Point_), intent(in)::point

      if (size(obj%center) /= size(point%coord)) then
         print *, "ERROR ::InitSetCenterCircle ::  ( size(obj%center)/=size(point%coord) )"
         return
      end if

      obj%center(:) = point%coord(:)

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

!#########################################################
   subroutine InitSetRadiusCircle(obj, radius)
      class(Circle_), intent(inout)::obj
      real(real64), intent(in)::radius

      obj%radius = radius

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

!#########################################################
   function getAreaCircle(obj) result(area)
      class(Circle_), intent(inout)::obj
      real(real64) :: area, pi
      pi = 3.14159260d0

      area = obj%radius*obj%radius*pi

   end function
!#########################################################

!#########################################################
   subroutine showCircle(obj, Name)
      class(Circle_), intent(in)::obj
      character(*), optional, intent(in)::Name
      real(real64) :: angle, dtheta, pi
      integer(int32) :: i

      pi = 3.14159260d0

      if (present(Name)) then
         open (10, file=Name)
         angle = 0.0d0
         dtheta = 2.0d0*pi/360.0
         do i = 1, 360
            write (10, *) obj%center(1) + obj%radius*cos(angle), obj%center(2) + obj%radius*sin(angle)
            angle = angle + dtheta
         end do
         close (10)
      else
         print *, "center> O(:) = ", obj%center(:)
         print *, "radius = ", obj%radius
      end if

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

!#########################################################
   subroutine InitSphere(obj, dim)
      class(Sphere_), intent(inout)::obj
      integer(int32), optional, intent(inout)::dim

      obj%center(:) = 0.0d0
      obj%radius = 1.0d0

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

!#########################################################
   subroutine InitSetCenterSphere(obj, point)
      class(Sphere_), intent(inout)::obj
      class(Point_), intent(in)::point

      if (size(obj%center) /= size(point%coord)) then
         print *, "ERROR ::InitSetCenterSphere ::  ( size(obj%center)/=size(point%coord) )"
         return
      end if

      obj%center(:) = point%coord(:)

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

!#########################################################
   subroutine InitSetRadiusSphere(obj, radius)
      class(Sphere_), intent(inout)::obj
      real(real64), intent(in)::radius

      obj%radius = radius

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

!#########################################################
   subroutine showSphere(obj, Name)
      class(Sphere_), intent(in)::obj
      character(*), optional, intent(in)::Name
      real(real64) :: angle1, angle2, dtheta, pi
      integer(int32) :: i, j

      pi = 3.14159260d0

      if (present(Name)) then
         open (10, file=Name)
         angle1 = 0.0d0
         angle2 = 0.0d0
         dtheta = 2.0d0*pi/360.0
         do i = 1, 360
            do j = 1, 360
               write (10, *) obj%center(1) + obj%radius*sin(angle1)*cos(angle2), &
                  obj%center(2) + obj%radius*sin(angle1)*sin(angle2), &
                  obj%center(3) + obj%radius*cos(angle1)
            end do
         end do
         close (10)
      else
         print *, "center> O(:) = ", obj%center(:)
         print *, "radius = ", obj%radius
      end if

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

!#########################################################
   subroutine InitTriangle(obj, dim)
      class(Triangle_), intent(inout)::obj
      integer(int32), optional, intent(in)::dim

      if (allocated(obj%NodCoord)) then
         deallocate (obj%NodCoord)
      end if
      allocate (obj%NodCoord(3, input(default=3, option=dim)))
      obj%NodCoord(:, :) = 0.0d0

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

!#########################################################
   subroutine setNodeTriangle(obj, point, order)
      class(Triangle_), intent(inout)::obj
      class(Point_), intent(in)::point
      integer(int32), intent(in)::order

      if (size(obj%NodCoord, 2) /= size(point%coord)) then
         print *, "ERROR ::InitSetNodeTriangle ::  ( size(obj%NodCoord,2)/=size(point%coord)  )"
         return
      end if
      obj%NodCoord(order, :) = point%coord(:)
   end subroutine
!#########################################################

!#########################################################
   subroutine importTriangle(obj, NodCoord, FileName)
      class(Triangle_), intent(inout)::obj
      integer(int32), optional, intent(in)::NodCoord(:, :)

      integer(int32) :: n, i
      character(*), optional, intent(in)::FileName

      if (present(FileName)) then
         open (30, file=FileName)
         read (30, *) n
         call obj%init(dim=n)
         do i = 1, size(obj%NodCoord, 1)
            read (30, *) obj%NodCoord(i, :)
         end do
         close (30)
         print *, "Imported triangle from ", FileName
         return
      end if

      if (present(NodCoord)) then

         if (size(obj%NodCoord, 2) /= size(NodCoord, 2)) then
            print *, "ERROR :: importTriangle :: size(obj%NodCoord,2)/=size(NodCoord,2)"
            return
         end if
         obj%NodCoord(1:3, :) = NodCoord(1:3, :)
      end if

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

!#########################################################
   subroutine getCircleTriangle(obj, type_of_circle, circle)
      class(Triangle_), intent(in)::obj
      type(Circle_), intent(inout) :: circle
      character(*), intent(in) :: type_of_circle
      real(real64) :: x1, x2, x3, y1, y2, y3
      real(real64), allocatable :: a(:), b(:), c(:)
      integer(int32) :: i, n

      n = size(obj%NodCoord, 2)
      call initCircle(circle, dim=n)

      if (type_of_circle == "center_of_gravity") then

         do i = 1, 3
            circle%center(:) = circle%center(:) + 1.0d0/3.0d0*obj%NodCoord(i, :)
         end do
         circle%radius = 1.0d0

      elseif (type_of_circle == "circumcenter") then

         ! 外心を計算する
         allocate (a(size(obj%NodCoord, 2)))
         a(:) = obj%NodCoord(1, :)
         x1 = obj%NodCoord(1, 1)
         y1 = obj%NodCoord(1, 2)
         x2 = obj%NodCoord(2, 1)
         y2 = obj%NodCoord(2, 2)
         x3 = obj%NodCoord(3, 1)
         y3 = obj%NodCoord(3, 2)
         circle%center(1) = 0.50d0*( &
                            (x1*x1 + y1*y1)*(y2 - y3) + &
                            (x2*x2 + y2*y2)*(y3 - y1) + &
                            (x3*x3 + y3*y3)*(y1 - y2) &
                            )
         circle%center(1) = circle%center(1)/( &
                            x1*(y2 - y3) + &
                            x2*(y3 - y1) + &
                            x3*(y1 - y2) &
                            )
         circle%center(2) = 0.50d0*( &
                            (x1*x1 + y1*y1)*(x2 - x3) + &
                            (x2*x2 + y2*y2)*(x3 - x1) + &
                            (x3*x3 + y3*y3)*(x1 - x2) &
                            )
         circle%center(2) = -circle%center(2)/( &
                            x1*(y2 - y3) + &
                            x2*(y3 - y1) + &
                            x3*(y1 - y2) &
                            )
         circle%radius = distance(circle%center, a)

      elseif (type_of_circle == "innter_center") then

         ! 外心を計算する
         allocate (a(size(obj%NodCoord, 2)))
         allocate (b(size(obj%NodCoord, 2)))
         allocate (c(size(obj%NodCoord, 2)))
         a(:) = obj%NodCoord(1, :)
         b(:) = obj%NodCoord(2, :)
         c(:) = obj%NodCoord(3, :)
         x1 = obj%NodCoord(1, 1)
         y1 = obj%NodCoord(1, 2)
         x2 = obj%NodCoord(2, 1)
         y2 = obj%NodCoord(2, 2)
         x3 = obj%NodCoord(3, 1)
         y3 = obj%NodCoord(3, 2)
         circle%center(1) = sqrt((x3 - x2)*(x3 - x2) + (y3 - y2)*(y3 - y2))*x1 + &
                            sqrt((x1 - x3)*(x1 - x3) + (y1 - y3)*(y1 - y3))*x2 + &
                            sqrt((x2 - x1)*(x2 - x1) + (y2 - y1)*(y2 - y1))*x3
         circle%center(1) = circle%center(1)/( &
                            sqrt((x2 - x1)*(x2 - x1) + (y2 - y1)*(y2 - y1)) + &
                            sqrt((x3 - x2)*(x3 - x2) + (y3 - y2)*(y3 - y2)) + &
                            sqrt((x1 - x3)*(x1 - x3) + (y1 - y3)*(y1 - y3)) &
                            )
         circle%center(2) = sqrt((x3 - x2)*(x3 - x2) + (y3 - y2)*(y3 - y2))*y1 + &
                            sqrt((x1 - x3)*(x1 - x3) + (y1 - y3)*(y1 - y3))*y2 + &
                            sqrt((x2 - x1)*(x2 - x1) + (y2 - y1)*(y2 - y1))*y3
         circle%center(2) = circle%center(2)/( &
                            sqrt((x2 - x1)*(x2 - x1) + (y2 - y1)*(y2 - y1)) + &
                            sqrt((x3 - x2)*(x3 - x2) + (y3 - y2)*(y3 - y2)) + &
                            sqrt((x1 - x3)*(x1 - x3) + (y1 - y3)*(y1 - y3)) &
                            )
         circle%radius = 2.0d0*obj%getArea()/(distance(a, b) + distance(b, c) + distance(c, a))

      elseif (type_of_circle == "excenter") then
         print *, "ERROR :: getCircleTriangle :: type_of_circle/excenter is now being implemented"
      elseif (type_of_circle == "orthocenter") then
         print *, "ERROR :: getCircleTriangle :: type_of_circle/orthocenter is now being implemented"
      end if

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

!#########################################################
   function getAreaTriangle(obj) result(area)
      class(Triangle_), intent(in)::obj
      real(real64) :: area
      real(real64) :: x1, x2, x3, y1, y2, y3
      real(real64) :: a(3), b(3), c(3)

      if (size(obj%NodCoord, 2) == 2) then

         x1 = obj%NodCoord(1, 1)
         y1 = obj%NodCoord(1, 2)
         x2 = obj%NodCoord(2, 1)
         y2 = obj%NodCoord(2, 2)
         x3 = obj%NodCoord(3, 1)
         y3 = obj%NodCoord(3, 2)
         area = abs(0.50d0*(x1*y2 + x2*y3 + x3*y1 - y1*x2 - y2*x3 - y3*x1))

      elseif (size(obj%NodCoord, 2) == 3) then
         a(:) = obj%NodCoord(1, :)
         b(:) = obj%NodCoord(2, :)
         c(:) = obj%NodCoord(3, :)
         a(:) = a(:) - c(:)
         b(:) = b(:) - c(:)
         area = (a(1)*a(1) + a(2)*a(2) + a(3)*a(3))*(b(1)*b(1) + b(2)*b(2) + b(3)*b(3))
         area = area - (a(1)*b(1) + a(2)*b(2) + a(3)*b(3))*(a(1)*b(1) + a(2)*b(2) + a(3)*b(3))
         area = 0.50d0*area
         area = sqrt(area)
      else
         print *, "ERROR :: getAreaTriangle, size(obj%NodCoord,2)==2 "
         return
      end if

   end function
!#########################################################

!#########################################################
   subroutine showTriangle(obj, Name, option)
      class(Triangle_), intent(in)::obj
      character(*), optional, intent(in)::Name, option

      if (present(Name)) then
         if (present(option)) then
            open (10, file=Name)
            write (10, *) obj%NodCoord(1, :)
            write (10, *) obj%NodCoord(2, :)
            write (10, *) obj%NodCoord(3, :)
            write (10, *) obj%NodCoord(1, :)
            close (10)
         else
            open (10, file=Name)
            write (10, *) obj%NodCoord(1, :), obj%NodCoord(2, :) - obj%NodCoord(1, :)
            write (10, *) obj%NodCoord(2, :), obj%NodCoord(3, :) - obj%NodCoord(2, :)
            write (10, *) obj%NodCoord(3, :), obj%NodCoord(1, :) - obj%NodCoord(3, :)
            close (10)
         end if
      else
         print *, "Triangle> x1(:) = ", obj%NodCoord(1, :)
         print *, "Triangle> x2(:) = ", obj%NodCoord(2, :)
         print *, "Triangle> x3(:) = ", obj%NodCoord(3, :)
      end if
   end subroutine
!#########################################################

!#########################################################
   subroutine GetOuterNormalTriangle(obj)
      class(Triangle_), intent(inout) :: obj
      real(real64) :: x1(3), x2(3), x3(3), on(3)
      integer(int32) :: n

      if (.not. allocated(obj%OuterNormal)) then
         n = size(obj%NodCoord, 2)
         allocate (obj%OuterNormal(n))
         obj%OuterNormal(:) = 0.0d0
      end if

      x1(:) = 0.0d0
      x2(:) = 0.0d0
      x1(:) = obj%NodCoord(2, 1:n) - obj%NodCoord(1, 1:n)
      x2(:) = obj%NodCoord(3, 1:n) - obj%NodCoord(1, 1:n)
      on(:) = 1.0d0/norm(cross_product(x1, x2))*cross_product(x1, x2)
      obj%OuterNormal(1:n) = on(1:n)

      x1(1:n) = obj%NodCoord(1, 1:n)
      x2(1:n) = obj%NodCoord(2, 1:n)
      x3(1:n) = obj%NodCoord(3, 1:n)
      if (allocated(obj%center)) then
         deallocate (obj%center)
      end if
      allocate (obj%center(n))
      obj%center(:) = 0.0d0
      obj%center(1:n) = obj%center(1:n) + x1(1:n)
      obj%center(1:n) = obj%center(1:n) + x2(1:n)
      obj%center(1:n) = obj%center(1:n) + x3(1:n)
      obj%center(1:n) = 1.0d0/3.0d0*obj%center(1:n)

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

!#########################################################
   subroutine InitRectangle(obj, dim)
      class(Rectangle_), intent(inout)::obj
      integer(int32), optional, intent(in)::dim

      if (allocated(obj%NodCoord)) then
         deallocate (obj%NodCoord)
      end if
      allocate (obj%NodCoord(4, input(default=3, option=dim)))
      obj%NodCoord(:, :) = 0.0d0

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

!#########################################################
   subroutine createRectangle(obj)
      class(Rectangle_), intent(inout) :: obj

      ! create unit one
      allocate (obj%NodCoord(8, 3))
      obj%NodCoord(1, 1) = -1.0d0; obj%NodCoord(1, 2) = -1.0d0; obj%NodCoord(1, 3) = -1.0d0; 
      obj%NodCoord(2, 1) = 1.0d0; obj%NodCoord(2, 2) = -1.0d0; obj%NodCoord(2, 3) = -1.0d0; 
      obj%NodCoord(3, 1) = 1.0d0; obj%NodCoord(3, 2) = 1.0d0; obj%NodCoord(3, 3) = -1.0d0; 
      obj%NodCoord(4, 1) = -1.0d0; obj%NodCoord(4, 2) = 1.0d0; obj%NodCoord(4, 3) = -1.0d0; 
      obj%NodCoord(5, 1) = -1.0d0; obj%NodCoord(5, 2) = -1.0d0; obj%NodCoord(5, 3) = 1.0d0; 
      obj%NodCoord(6, 1) = 1.0d0; obj%NodCoord(6, 2) = -1.0d0; obj%NodCoord(6, 3) = 1.0d0; 
      obj%NodCoord(7, 1) = 1.0d0; obj%NodCoord(7, 2) = 1.0d0; obj%NodCoord(7, 3) = 1.0d0; 
      obj%NodCoord(8, 1) = -1.0d0; obj%NodCoord(8, 2) = 1.0d0; obj%NodCoord(8, 3) = 1.0d0; 
   end subroutine createRectangle
!#########################################################

!#########################################################
   subroutine moveRectangle(obj, x, y, z)
      class(Rectangle_), intent(inout) :: obj
      real(real64), optional, intent(in) :: x, y, z

      obj%NodCoord(:, 1) = obj%NodCoord(:, 1) + input(default=0.0d0, option=x)
      obj%NodCoord(:, 2) = obj%NodCoord(:, 2) + input(default=0.0d0, option=y)
      obj%NodCoord(:, 3) = obj%NodCoord(:, 3) + input(default=0.0d0, option=z)

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

!#########################################################
   subroutine setNodeRectangle(obj, point, order)
      class(Rectangle_), intent(inout)::obj
      class(Point_), intent(in)::point
      integer(int32), intent(in)::order

      if (size(obj%NodCoord, 2) /= size(point%coord)) then
         print *, "ERROR ::InitSetNodeRectangle ::  ( size(obj%NodCoord,2)/=size(point%coord)  )"
         return
      end if
      obj%NodCoord(order, :) = point%coord(:)
   end subroutine
!#########################################################

!#########################################################
   subroutine importRectangle(obj, NodCoord, FileName)
      class(Rectangle_), intent(inout)::obj
      integer(int32), optional, intent(in)::NodCoord(:, :)
      integer(int32) :: n, i
      character(*), optional, intent(in)::FileName

      if (present(FileName)) then
         open (30, file=FileName)
         read (30, *) n
         call obj%init(dim=n)
         do i = 1, size(obj%NodCoord, 1)
            read (30, *) obj%NodCoord(i, :)
         end do
         close (30)
         print *, "Imported rectangle from ", FileName
         return
      end if

      if (present(NodCoord)) then
         if (size(obj%NodCoord, 2) /= size(NodCoord, 2)) then
            print *, "ERROR :: importRectangle :: size(obj%NodCoord,2)/=size(NodCoord,2)"
            return
         end if
         obj%NodCoord(1:4, :) = NodCoord(1:4, :)
      end if
   end subroutine
!#########################################################

!#########################################################
   subroutine getCircleRectangle(obj, type_of_circle, circle)
      class(Rectangle_), intent(in)::obj
      type(Circle_), intent(inout) :: circle
      character(*), intent(in) :: type_of_circle
      real(real64) :: x1, x2, x3, x4, y1, y2, y3, y4
      real(real64), allocatable :: a(:), b(:), c(:)
      integer(int32) :: i, n

      n = size(obj%NodCoord, 2)
      call initCircle(circle, dim=n)

      if (type_of_circle == "center_of_gravity") then

         do i = 1, 4
            circle%center(:) = circle%center(:) + 1.0d0/4.0d0*obj%NodCoord(i, :)
         end do
         circle%radius = 1.0d0

      elseif (type_of_circle == "circumcenter") then

         ! 外心を計算する
         allocate (a(size(obj%NodCoord, 2)))
         a(:) = obj%NodCoord(1, :)
         x1 = obj%NodCoord(1, 1)
         y1 = obj%NodCoord(1, 2)
         x2 = obj%NodCoord(2, 1)
         y2 = obj%NodCoord(2, 2)
         x3 = obj%NodCoord(3, 1)
         y3 = obj%NodCoord(3, 2)
         x4 = obj%NodCoord(4, 1)
         y4 = obj%NodCoord(4, 2)

         circle%center(1) = 0.50d0*( &
                            (x1*x1 + y1*y1)*(y2 - y3) + &
                            (x2*x2 + y2*y2)*(y3 - y1) + &
                            (x3*x3 + y3*y3)*(y1 - y2) &
                            )
         circle%center(1) = circle%center(1)/( &
                            x1*(y2 - y3) + &
                            x2*(y3 - y1) + &
                            x3*(y1 - y2) &
                            )
         circle%center(2) = 0.50d0*( &
                            (x1*x1 + y1*y1)*(x2 - x3) + &
                            (x2*x2 + y2*y2)*(x3 - x1) + &
                            (x3*x3 + y3*y3)*(x1 - x2) &
                            )
         circle%center(2) = -circle%center(2)/( &
                            x1*(y2 - y3) + &
                            x2*(y3 - y1) + &
                            x3*(y1 - y2) &
                            )
         circle%radius = distance(circle%center, a)

      elseif (type_of_circle == "innter_center") then

         ! 外心を計算する
         allocate (a(size(obj%NodCoord, 2)))
         allocate (b(size(obj%NodCoord, 2)))
         allocate (c(size(obj%NodCoord, 2)))
         a(:) = obj%NodCoord(1, :)
         b(:) = obj%NodCoord(2, :)
         c(:) = obj%NodCoord(3, :)
         x1 = obj%NodCoord(1, 1)
         y1 = obj%NodCoord(1, 2)
         x2 = obj%NodCoord(2, 1)
         y2 = obj%NodCoord(2, 2)
         x3 = obj%NodCoord(3, 1)
         y3 = obj%NodCoord(3, 2)
         circle%center(1) = sqrt((x3 - x2)*(x3 - x2) + (y3 - y2)*(y3 - y2))*x1 + &
                            sqrt((x1 - x3)*(x1 - x3) + (y1 - y3)*(y1 - y3))*x2 + &
                            sqrt((x2 - x1)*(x2 - x1) + (y2 - y1)*(y2 - y1))*x3
         circle%center(1) = circle%center(1)/( &
                            sqrt((x2 - x1)*(x2 - x1) + (y2 - y1)*(y2 - y1)) + &
                            sqrt((x3 - x2)*(x3 - x2) + (y3 - y2)*(y3 - y2)) + &
                            sqrt((x1 - x3)*(x1 - x3) + (y1 - y3)*(y1 - y3)) &
                            )
         circle%center(2) = sqrt((x3 - x2)*(x3 - x2) + (y3 - y2)*(y3 - y2))*y1 + &
                            sqrt((x1 - x3)*(x1 - x3) + (y1 - y3)*(y1 - y3))*y2 + &
                            sqrt((x2 - x1)*(x2 - x1) + (y2 - y1)*(y2 - y1))*y3
         circle%center(2) = circle%center(2)/( &
                            sqrt((x2 - x1)*(x2 - x1) + (y2 - y1)*(y2 - y1)) + &
                            sqrt((x3 - x2)*(x3 - x2) + (y3 - y2)*(y3 - y2)) + &
                            sqrt((x1 - x3)*(x1 - x3) + (y1 - y3)*(y1 - y3)) &
                            )
         circle%radius = 2.0d0*obj%getArea()/(distance(a, b) + distance(b, c) + distance(c, a))

      elseif (type_of_circle == "excenter") then
         print *, "ERROR :: getCircleRectangle :: type_of_circle/excenter is now being implemented"
      elseif (type_of_circle == "orthocenter") then
         print *, "ERROR :: getCircleRectangle :: type_of_circle/orthocenter is now being implemented"
      end if

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

!#########################################################
   function getAreaRectangle(obj) result(area)
      class(Rectangle_), intent(in)::obj
      real(real64) :: area
      real(real64) :: x1, x2, x3, x4, y1, y2, y3, y4

      if (size(obj%NodCoord, 2) == 2) then

         x1 = obj%NodCoord(1, 1)
         y1 = obj%NodCoord(1, 2)
         x2 = obj%NodCoord(2, 1)
         y2 = obj%NodCoord(2, 2)
         x3 = obj%NodCoord(3, 1)
         y3 = obj%NodCoord(3, 2)
         area = abs(0.50d0*(x1*y2 + x2*y3 + x3*y1 - y1*x2 - y2*x3 - y3*x1))

         x2 = obj%NodCoord(1, 1)
         y2 = obj%NodCoord(1, 2)
         x3 = obj%NodCoord(2, 1)
         y3 = obj%NodCoord(2, 2)
         x4 = obj%NodCoord(3, 1)
         y4 = obj%NodCoord(3, 2)
         area = area + abs(0.50d0*(x1*y2 + x2*y3 + x3*y1 - y1*x2 - y2*x3 - y3*x1))

      elseif (size(obj%NodCoord, 2) == 3) then
         print *, " getAreaRectangle >> not ready for 3D"
         stop
      else
         print *, "ERROR :: getAreaRectangle, size(obj%NodCoord,2)==2 "
         return
      end if

   end function
!#########################################################

!#########################################################
   subroutine showRectangle(obj, Name, option)
      class(Rectangle_), intent(in)::obj
      character(*), optional, intent(in)::Name, option

      if (present(Name)) then
         if (present(option)) then
            open (10, file=Name)
            write (10, *) obj%NodCoord(1, :)
            write (10, *) obj%NodCoord(2, :)
            write (10, *) obj%NodCoord(3, :)
            write (10, *) obj%NodCoord(4, :)
            write (10, *) obj%NodCoord(1, :)
            close (10)
         else
            open (10, file=Name)
            write (10, *) obj%NodCoord(1, :), obj%NodCoord(2, :) - obj%NodCoord(1, :)
            write (10, *) obj%NodCoord(2, :), obj%NodCoord(3, :) - obj%NodCoord(2, :)
            write (10, *) obj%NodCoord(3, :), obj%NodCoord(4, :) - obj%NodCoord(3, :)
            write (10, *) obj%NodCoord(4, :), obj%NodCoord(1, :) - obj%NodCoord(4, :)
            close (10)
         end if
      else
         print *, "Rectangle> x1(:) = ", obj%NodCoord(1, :)
         print *, "Rectangle> x2(:) = ", obj%NodCoord(2, :)
         print *, "Rectangle> x3(:) = ", obj%NodCoord(3, :)
         print *, "Rectangle> x3(:) = ", obj%NodCoord(4, :)
      end if
   end subroutine
!#########################################################

!#########################################################
   function contactRectangle(obj, Rectangle, threshold) result(contact)
      class(Rectangle_), intent(in) :: obj
      class(Rectangle_), intent(in) :: Rectangle
      integer(int32), optional, intent(in) :: threshold
      logical :: contact, inside
      integer(int32) :: i, j, k, n, m, dim, inside_score, th
      real(real64)::dist1, dist2, dist3, x2d(2), x3d(3), x_max(3), x_min(3)

      th = input(default=1, option=threshold)

      ! Contact Rectangle to Rectangle
      dim = size(obj%NodCoord, 2)
      if (dim == 3) then
         ! get range
         do i = 1, dim
            x_max(i) = maxval(Rectangle%NodCoord(:, i))
            x_min(i) = minval(Rectangle%NodCoord(:, i))
         end do
         inside_score = 0
         inside = .false.
         do i = 1, size(obj%NodCoord, 1)
            x3d(1:3) = obj%NodCoord(i, 1:3)
            inside = InOrOut(x=x3d, xmax=x_max, xmin=x_min, DimNum=dim)
            if (inside .eqv. .true.) then
               inside_score = inside_score + 1
            end if
         end do
         if (inside_score >= th) then
            contact = .true.
         else
            contact = .false.
         end if
      else
         print *, "Please implement contactRectangle for", dim, "D"
         stop
      end if

   end function contactRectangle
!#########################################################

!#########################################################
   subroutine InitTetrahedron(obj, NodCoord)
      class(Tetrahedron_), intent(inout) :: obj
      real(real64), intent(in) :: NodCoord(4, 3)

      obj%NodCoord(:, :) = NodCoord(:, :)

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

!#########################################################
   subroutine getCircleTetrahedron(obj)
      class(Tetrahedron_), intent(inout) :: obj
      real(real64) :: a(3), b(3), c(3), d(3), e(3), f(3), g(3), s, t, r, N(3)
      real(real64) :: a_(3), b_(3), c_(3), d_(3), aA, aB, aC, aD, V, k

      a(:) = obj%NodCoord(1, :)
      b(:) = obj%NodCoord(2, :)
      c(:) = obj%NodCoord(3, :)
      d(:) = obj%NodCoord(4, :)

      a_(:) = a(:) - d(:)
      b_(:) = b(:) - d(:)
      c_(:) = c(:) - d(:)
      d_(:) = 0.0d0

      aA = 0.50d0*norm(cross_product(a_, b_))
      aB = 0.50d0*norm(cross_product(b_, c_))
      aC = 0.50d0*norm(cross_product(c_, d_))
      aD = 0.50d0*norm(cross_product(d_, a_))

      V = 1.0d0/6.0d0*dot_product(cross_product(a_(:), b_(:)), c_(:))
      r = 3.0d0*V/(aA + aB + aC + aD)
      obj%radius = r

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

!#########################################################
!subroutine
!end subroutine
!#########################################################
end module GeometryClass
