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