GeometryClass.f90 Source File


Contents

Source Code


Source Code

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)
    endif
    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(:)
    endif
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
        endif
        obj%coord(:)=xvec(:)
        return
    endif

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


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

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

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

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)
    endif
    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
    endif

    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(:)
    endif

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
    endif
    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,:)
    endif
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)
    endif
    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
    endif

    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
        enddo
        close(10)
    else
        print *, "center> O(:) = ", obj%center(:)
        print *, "radius = ", obj%radius
    endif

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
    endif

    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)
            enddo
        enddo
        close(10)
    else
        print *, "center> O(:) = ",obj%center(:)
        print *, "radius = ",obj%radius
    endif
    
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)
    endif
    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
    endif
    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=trim(FileName))
        read(30,*) n 
        call obj%init(dim=n)
        do i=1,size(obj%NodCoord,1)
            read(30,*) obj%NodCoord(i,:)
        enddo
        close(30)
        print *, "Imported triangle from ",trim(FileName)
        return
    endif

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

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,:)
        enddo
        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"
    endif

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
    endif

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)
        endif
    else
        print *, "Triangle> x1(:) = ",obj%NodCoord(1,:)
        print *, "Triangle> x2(:) = ",obj%NodCoord(2,:)
        print *, "Triangle> x3(:) = ",obj%NodCoord(3,:)
    endif
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
    endif

    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)
    endif
    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)
    endif
    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
    endif
    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=trim(FileName))
        read(30,*) n 
        call obj%init(dim=n)
        do i=1,size(obj%NodCoord,1)
            read(30,*) obj%NodCoord(i,:)
        enddo
        close(30)
        print *, "Imported rectangle from ",trim(FileName)
        return
    endif


    if(present(NodCoord) )then
        if(size(obj%NodCoord,2)/=size(NodCoord,2) )then
            print *, "ERROR :: importRectangle :: size(obj%NodCoord,2)/=size(NodCoord,2)"
            return
        endif
        obj%NodCoord(1:4,:)=NodCoord(1:4,:)
    endif
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,:)
        enddo
        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"
    endif

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
    endif

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)
        endif
    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,:)
    endif
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))
        enddo
        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
            endif
        enddo
        if(inside_score >= th)then
            contact=.true.
        else
            contact=.false.
        endif
    else
        print *, "Please implement contactRectangle for",dim,"D"
        stop 
    endif
    

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