module ShapeFunctionClass use, intrinsic :: iso_fortran_env use MathClass use ArrayClass use IOClass implicit none type::ShapeFunction_ real(real64), allocatable::Nmat(:) real(real64), allocatable::dNdgzi(:, :) real(real64), allocatable::dNdgzidgzi(:, :) real(real64), allocatable::gzi(:) real(real64), allocatable::GaussPoint(:, :) real(real64), allocatable::GaussIntegWei(:) real(real64), allocatable::Jmat(:, :) ! d x/d xi real(real64), allocatable::JmatInv(:, :)! d xi/d x real(real64), allocatable::ElemCoord(:, :) real(real64), allocatable::ElemCoord_n(:, :) real(real64), allocatable::du(:, :) real(real64) :: detJ integer(int32) :: NumOfNode integer(int32) :: NumOfOrder integer(int32) :: NumOfDim = 0 integer(int32) :: NumOfGp = 0 integer(int32) :: GpID integer(int32) :: ierr integer(int32) :: currentGpID integer(int32) :: ElementID logical :: ReducedIntegration = .false. logical :: Empty = .true. character*70::ElemType character(len=60):: ErrorMsg contains procedure :: init => initShapeFunction procedure :: update => updateShapeFunction procedure :: SetType => SetShapeFuncType procedure :: GetAll => GetAllShapeFunc procedure :: get => GetAllShapeFunc procedure :: getOnlyNvec => GetShapeFunction procedure :: Deallocate => DeallocateShapeFunction procedure :: getType => getShapeFuncType procedure :: GetGaussPoint => GetGaussPoint procedure :: export => exportShapeFunction procedure :: remove => removeShapeFunction procedure :: save => saveShapeFunction procedure :: open => openShapeFunction end type ShapeFunction_ contains !! ##################################################### subroutine openShapeFunction(obj, path, name) class(ShapeFunction_), intent(inout)::obj character(*), intent(in) :: path character(*), optional, intent(in) :: name type(IO_) :: f if (present(name)) then call execute_command_line("mkdir -p "//path//"/"//name) call f%open(path//"/"//name//"/", "ShapeFunction", "prop") call openArray(f%fh, obj%Nmat) call openArray(f%fh, obj%dNdgzi) call openArray(f%fh, obj%dNdgzidgzi) call openArray(f%fh, obj%gzi) call openArray(f%fh, obj%GaussPoint) call openArray(f%fh, obj%GaussIntegWei) call openArray(f%fh, obj%Jmat) call openArray(f%fh, obj%JmatInv) call openArray(f%fh, obj%ElemCoord) call openArray(f%fh, obj%ElemCoord_n) call openArray(f%fh, obj%du) read (f%fh, *) obj%detJ read (f%fh, *) obj%NumOfNode read (f%fh, *) obj%NumOfOrder read (f%fh, *) obj%NumOfDim read (f%fh, *) obj%NumOfGp read (f%fh, *) obj%GpID read (f%fh, *) obj%ierr read (f%fh, *) obj%currentGpID read (f%fh, *) obj%ReducedIntegration read (f%fh, *) obj%ElemType read (f%fh, *) obj%ErrorMsg call f%close() else call execute_command_line("mkdir -p "//path//"/ShapeFunction") call f%open(path//"/", "ShapeFunction/ShapeFunction", "prop") call openArray(f%fh, obj%Nmat) call openArray(f%fh, obj%dNdgzi) call openArray(f%fh, obj%dNdgzidgzi) call openArray(f%fh, obj%gzi) call openArray(f%fh, obj%GaussPoint) call openArray(f%fh, obj%GaussIntegWei) call openArray(f%fh, obj%Jmat) call openArray(f%fh, obj%JmatInv) call openArray(f%fh, obj%ElemCoord) call openArray(f%fh, obj%ElemCoord_n) call openArray(f%fh, obj%du) read (f%fh, *) obj%detJ read (f%fh, *) obj%NumOfNode read (f%fh, *) obj%NumOfOrder read (f%fh, *) obj%NumOfDim read (f%fh, *) obj%NumOfGp read (f%fh, *) obj%GpID read (f%fh, *) obj%ierr read (f%fh, *) obj%currentGpID read (f%fh, *) obj%ReducedIntegration read (f%fh, *) obj%ElemType read (f%fh, *) obj%ErrorMsg call f%close() end if end subroutine !! ##################################################### !! ##################################################### subroutine saveShapeFunction(obj, path, name) class(ShapeFunction_), intent(inout)::obj character(*), intent(in) :: path character(*), optional, intent(in) :: name type(IO_) :: f if (present(name)) then call execute_command_line("mkdir -p "//path//"/"//name) call f%open(path//"/"//name//"/", "ShapeFunction", "prop") call writeArray(f%fh, obj%Nmat) call writeArray(f%fh, obj%dNdgzi) call writeArray(f%fh, obj%dNdgzidgzi) call writeArray(f%fh, obj%gzi) call writeArray(f%fh, obj%GaussPoint) call writeArray(f%fh, obj%GaussIntegWei) call writeArray(f%fh, obj%Jmat) call writeArray(f%fh, obj%JmatInv) call writeArray(f%fh, obj%ElemCoord) call writeArray(f%fh, obj%ElemCoord_n) call writeArray(f%fh, obj%du) write (f%fh, *) obj%detJ write (f%fh, *) obj%NumOfNode write (f%fh, *) obj%NumOfOrder write (f%fh, *) obj%NumOfDim write (f%fh, *) obj%NumOfGp write (f%fh, *) obj%GpID write (f%fh, *) obj%ierr write (f%fh, *) obj%currentGpID write (f%fh, *) obj%ReducedIntegration write (f%fh, *) obj%ElemType write (f%fh, *) obj%ErrorMsg call f%close() else call execute_command_line("mkdir -p "//path//"/ShapeFunction") call f%open(path//"/", "ShapeFunction/ShapeFunction", "prop") call writeArray(f%fh, obj%Nmat) call writeArray(f%fh, obj%dNdgzi) call writeArray(f%fh, obj%dNdgzidgzi) call writeArray(f%fh, obj%gzi) call writeArray(f%fh, obj%GaussPoint) call writeArray(f%fh, obj%GaussIntegWei) call writeArray(f%fh, obj%Jmat) call writeArray(f%fh, obj%JmatInv) call writeArray(f%fh, obj%ElemCoord) call writeArray(f%fh, obj%ElemCoord_n) call writeArray(f%fh, obj%du) write (f%fh, *) obj%detJ write (f%fh, *) obj%NumOfNode write (f%fh, *) obj%NumOfOrder write (f%fh, *) obj%NumOfDim write (f%fh, *) obj%NumOfGp write (f%fh, *) obj%GpID write (f%fh, *) obj%ierr write (f%fh, *) obj%currentGpID write (f%fh, *) obj%ReducedIntegration write (f%fh, *) obj%ElemType write (f%fh, *) obj%ErrorMsg call f%close() end if end subroutine !! ##################################################### subroutine removeShapeFunction(obj) class(ShapeFunction_), intent(inout)::obj if (allocated(obj%Nmat)) then deallocate (obj%Nmat) end if if (allocated(obj%dNdgzi)) then deallocate (obj%dNdgzi) end if if (allocated(obj%dNdgzidgzi)) then deallocate (obj%dNdgzidgzi) end if if (allocated(obj%gzi)) then deallocate (obj%gzi) end if if (allocated(obj%GaussPoint)) then deallocate (obj%GaussPoint) end if if (allocated(obj%GaussIntegWei)) then deallocate (obj%GaussIntegWei) end if if (allocated(obj%Jmat)) then deallocate (obj%Jmat) end if if (allocated(obj%ElemCoord)) then deallocate (obj%ElemCoord) end if if (allocated(obj%ElemCoord_n)) then deallocate (obj%ElemCoord_n) end if if (allocated(obj%du)) then deallocate (obj%du) end if obj%detJ = 0.0d0 obj%NumOfNode = 0 obj%NumOfOrder = 0 obj%NumOfDim = 0 obj%NumOfGp = 0 obj%GpID = 0 obj%ierr = 0 obj%currentGpID = 0 obj%ReducedIntegration = .false. obj%ElemType = " " obj%ErrorMsg = " " end subroutine !! ################################################## subroutine initShapeFunction(obj, ElemType) class(ShapeFunction_), intent(inout) :: obj character(*), intent(in), optional :: ElemType obj%ElemType = ElemType call obj%SetType() end subroutine !! ################################################## !! ################################################## subroutine updateShapeFunction(obj, ElemType, NodCoord, ElemNod, ElemID, GpID) class(ShapeFunction_), intent(inout) :: obj character(*), optional, intent(in) :: ElemType integer(int32), intent(in) ::ElemNod(:, :), ElemID, GpID real(real64), intent(in) ::NodCoord(:, :) integer(int32) :: i, j if (present(ElemType)) then call obj%init(ElemType) end if call obj%GetAll(elem_id=i, nod_coord=NodCoord, elem_nod=ElemNod, OptionalGpID=j) end subroutine !! ################################################## !! ################################################## subroutine SetShapeFuncType(obj, NumOfDim, NumOfNodePerElem, ReducedIntegration, NumOfGp) class(ShapeFunction_), intent(inout)::obj logical, optional, intent(in) :: ReducedIntegration character(:), allocatable :: TrimedElemType integer(int32), optional, intent(in) :: NumOfDim, NumOfNodePerElem integer(int32), optional, intent(in) :: NumOfGp if (present(NumOfDim)) then if (present(NumOfNodePerElem)) then call obj%getType(NumOfDim, NumOfNodePerElem, NumOfGp) end if end if if (present(ReducedIntegration)) then obj%ReducedIntegration = ReducedIntegration end if TrimedElemType = obj%ElemType if (TrimedElemType == "LinearRectangularGp4") then obj%NumOfNode = 4 obj%NumOfOrder = 1 obj%NumOfDim = 2 obj%NumOfGp = 4 elseif (TrimedElemType == "LinearHexahedralGp8") then obj%NumOfNode = 8 obj%NumOfOrder = 1 obj%NumOfDim = 3 obj%NumOfGp = 8 elseif (TrimedElemType == "Triangular") then obj%NumOfNode = 3 obj%NumOfOrder = 1 obj%NumOfDim = 2 obj%NumOfGp = 1 elseif (TrimedElemType == "LinearTetrahedral") then obj%NumOfNode = 4 obj%NumOfOrder = 1 obj%NumOfDim = 3 obj%NumOfGp = 1 elseif (TrimedElemType == "QuadHexahedralGp8") then obj%NumOfNode = 20 obj%NumOfOrder = 2 obj%NumOfDim = 3 obj%NumOfGp = 8 elseif (TrimedElemType == "QuadHexahedralGp27") then obj%NumOfNode = 20 obj%NumOfOrder = 2 obj%NumOfDim = 3 obj%NumOfGp = 27 else obj%ErrorMsg = "ShapeFunctionClass.f90 :: SetShapeFuncType >> Element : "//TrimedElemType//"is not defined." print *, "ShapeFunctionClass.f90 :: SetShapeFuncType >> Element : ", TrimedElemType, "is not defined." return end if end subroutine SetShapeFuncType !! ################################################## !! ################################################## subroutine getShapeFuncType(obj, NumOfDim, NumOfNodePerElem, NumOfGp) class(ShapeFunction_), intent(inout)::obj character*70 :: TrimedElemType integer(int32), intent(in) :: NumOfDim, NumOfNodePerElem integer(int32), optional, intent(in) :: NumOfGp if (NumOfDim == 2) then if (NumOfNodePerElem == 4) then obj%ElemType = "LinearRectangularGp4" else print *, "For 2D, NumOfNodePerElem = ", NumOfNodePerElem, " is not set." end if elseif (NumOfDim == 3) then if (NumOfNodePerElem == 8) then obj%ElemType = "LinearHexahedralGp8" elseif (NumOfNodePerElem == 20) then if (present(NumOfGP)) then if (NumOfGp == 8) then obj%ElemType = "QuadHexahedralGp8" elseif (NumOfGp == 27) then obj%ElemType = "QuadHexahedralGp27" else print *, "For 3D, 20-noded element,Number of Gauss point ", NumOfGP, " is not ready." end if else obj%ElemType = "QuadHexahedralGp8" end if else print *, "For 3D, NumOfNodePerElem = ", NumOfNodePerElem, " is not set." end if else print *, "NumOfDim = ", NumOfDim, " is not set." end if end subroutine getShapeFuncType !! ################################################## !! ################################################## subroutine GetAllShapeFunc(obj, elem_id, nod_coord, nod_coord_n, elem_nod, OptionalNumOfNode, OptionalNumOfOrder, & OptionalNumOfDim, OptionalNumOfGp, OptionalGpID, ReducedIntegration, NumOfDim, NumOfNodePerElem) class(ShapeFunction_), intent(inout)::obj integer(int32), optional, intent(in)::elem_nod(:, :), elem_id integer(int32), optional, intent(in) :: NumOfDim, NumOfNodePerElem real(real64), optional, intent(in)::nod_coord(:, :), nod_coord_n(:, :) integer(int32), optional, intent(in)::OptionalNumOfNode, OptionalNumOfOrder, & OptionalNumOfDim, OptionalNumOfGp, OptionalGpID logical, optional, intent(in) :: ReducedIntegration integer(int32) :: nd, ne nd = input(default=0, option=NumOfDim) ne = input(default=0, option=NumOfNodePerElem) if (present(nod_coord)) then nd = size(nod_coord, 2) end if if (present(nod_coord)) then ne = size(elem_nod, 2) end if if (obj%NumOfGp == 0 .and. obj%NumOfDim == 0) then if (nd == 0 .or. ne == 0) then print *, "ERROR :: GetAllShapeFunc >> please import NumOfDim and NumOfNodePerElem" stop end if call obj%SetType(NumOfDim=nd, NumOfNodePerElem=ne, ReducedIntegration=ReducedIntegration) end if if (present(ReducedIntegration)) then obj%ReducedIntegration = ReducedIntegration end if if (present(OptionalNumOfNode)) then obj%NumOfNode = OptionalNumOfNode end if if (present(OptionalNumOfOrder)) then obj%NumOfOrder = OptionalNumOfOrder end if if (present(OptionalNumOfDim)) then obj%NumOfDim = OptionalNumOfDim end if if (present(OptionalNumOfGp)) then obj%NumOfGp = OptionalNumOfGp end if if (present(OptionalGpID)) then obj%GpID = OptionalGpID end if call GetGaussPoint(obj) call SetGaussPoint(obj) call GetShapeFunction(obj) call GetShapeFuncDer1(obj) if (.not. present(nod_coord) .or. .not. present(elem_nod)) then obj%ErrorMsg = "Mesh%NodCoord and Mesh%ElemNod is necessary to get Jmat" print *, obj%ErrorMsg return end if call GetElemCoord(obj, nod_coord, elem_nod, elem_id) if (present(nod_coord_n)) then call GetElemCoord_n(obj, nod_coord_n, elem_nod, elem_id) call Getdu(obj) end if call GetJmat(obj) obj%detJ = det_mat(obj%Jmat, size(obj%Jmat, 1)) end subroutine GetAllShapeFunc !! ################################################## !! ################################################## subroutine DeallocateShapeFunction(obj) class(ShapeFunction_), intent(inout)::obj if (allocated(obj%Nmat)) deallocate (obj%Nmat) if (allocated(obj%dNdgzi)) deallocate (obj%dNdgzi) if (allocated(obj%dNdgzidgzi)) deallocate (obj%dNdgzidgzi) if (allocated(obj%gzi)) deallocate (obj%gzi) if (allocated(obj%GaussPoint)) deallocate (obj%GaussPoint) if (allocated(obj%GaussIntegWei)) deallocate (obj%GaussIntegWei) if (allocated(obj%Jmat)) deallocate (obj%Jmat) obj%ErrorMsg = "All allocatable entities are deallocated" end subroutine DeallocateShapeFunction !! ################################################## !! ###################################### subroutine GetGaussPoint(obj) class(ShapeFunction_), intent(inout)::obj real(real64) :: pval(1:3), wval(1:3) if (.not. allocated(obj%GaussPoint)) then allocate (obj%GaussPoint(obj%NumOfDim, obj%NumOfGp)) else if (size(obj%GaussPoint, 1) /= obj%NumOfDim .or. & size(obj%GaussPoint, 2) /= obj%NumOfGp) then deallocate (obj%GaussPoint) allocate (obj%GaussPoint(obj%NumOfDim, obj%NumOfGp)) end if end if if (.not. allocated(obj%GaussIntegWei)) then allocate (obj%GaussIntegWei(obj%NumOfGp)) else if (size(obj%GaussIntegWei) /= obj%NumOfGp) then deallocate (obj%GaussIntegWei) allocate (obj%GaussIntegWei(obj%NumOfGp)) end if end if if (size(obj%GaussPoint, 1) == 2 .and. size(obj%GaussPoint, 2) == 4) then obj%GaussPoint(1, 1) = -0.57735026918962576d0 obj%GaussPoint(1, 2) = 0.57735026918962576d0 obj%GaussPoint(1, 3) = 0.57735026918962576d0 obj%GaussPoint(1, 4) = -0.57735026918962576d0 obj%GaussPoint(2, 1) = -0.57735026918962576d0 obj%GaussPoint(2, 2) = -0.57735026918962576d0 obj%GaussPoint(2, 3) = 0.57735026918962576d0 obj%GaussPoint(2, 4) = 0.57735026918962576d0 obj%GaussIntegWei(:) = 1.0d0 if (obj%ReducedIntegration .eqv. .true.) then obj%GaussPoint(:, :) = 0.0d0 obj%GaussIntegWei(:) = 0.250d0 end if elseif (size(obj%GaussPoint, 1) == 3 .and. size(obj%GaussPoint, 2) == 8) then obj%GaussPoint(1, 1) = -0.57735026918962576d0 obj%GaussPoint(1, 2) = 0.57735026918962576d0 obj%GaussPoint(1, 3) = 0.57735026918962576d0 obj%GaussPoint(1, 4) = -0.57735026918962576d0 obj%GaussPoint(1, 5) = -0.57735026918962576d0 obj%GaussPoint(1, 6) = 0.57735026918962576d0 obj%GaussPoint(1, 7) = 0.57735026918962576d0 obj%GaussPoint(1, 8) = -0.57735026918962576d0 obj%GaussPoint(2, 1) = -0.57735026918962576d0 obj%GaussPoint(2, 2) = -0.57735026918962576d0 obj%GaussPoint(2, 3) = 0.57735026918962576d0 obj%GaussPoint(2, 4) = 0.57735026918962576d0 obj%GaussPoint(2, 5) = -0.57735026918962576d0 obj%GaussPoint(2, 6) = -0.57735026918962576d0 obj%GaussPoint(2, 7) = 0.57735026918962576d0 obj%GaussPoint(2, 8) = 0.57735026918962576d0 obj%GaussPoint(3, 1) = -0.57735026918962576d0 obj%GaussPoint(3, 2) = -0.57735026918962576d0 obj%GaussPoint(3, 3) = -0.57735026918962576d0 obj%GaussPoint(3, 4) = -0.57735026918962576d0 obj%GaussPoint(3, 5) = 0.57735026918962576d0 obj%GaussPoint(3, 6) = 0.57735026918962576d0 obj%GaussPoint(3, 7) = 0.57735026918962576d0 obj%GaussPoint(3, 8) = 0.57735026918962576d0 obj%GaussIntegWei(:) = 1.0d0 if (obj%ReducedIntegration .eqv. .true.) then obj%GaussPoint(:, :) = 0.0d0 obj%GaussIntegWei(:) = 0.1250d0 end if elseif (size(obj%GaussPoint, 1) == 3 .and. size(obj%GaussPoint, 2) == 27) then pval = [-sqrt(3.0d0/5.0d0), 0.0d0, sqrt(3.0d0/5.0d0)] wval = [5.0d0/9.0d0, 8.0d0/9.0d0, 5.0d0/9.0d0] obj%GaussPoint(1:3, 1) = [pval(1), pval(1), pval(1)]; obj%GaussIntegWei(1) = (wval(1)*wval(1)*wval(1)) obj%GaussPoint(1:3, 2) = [pval(2), pval(1), pval(1)]; obj%GaussIntegWei(2) = (wval(2)*wval(1)*wval(1)) obj%GaussPoint(1:3, 3) = [pval(3), pval(1), pval(1)]; obj%GaussIntegWei(3) = (wval(3)*wval(1)*wval(1)) obj%GaussPoint(1:3, 4) = [pval(1), pval(2), pval(1)]; obj%GaussIntegWei(4) = (wval(1)*wval(2)*wval(1)) obj%GaussPoint(1:3, 5) = [pval(2), pval(2), pval(1)]; obj%GaussIntegWei(5) = (wval(2)*wval(2)*wval(1)) obj%GaussPoint(1:3, 6) = [pval(3), pval(2), pval(1)]; obj%GaussIntegWei(6) = (wval(3)*wval(2)*wval(1)) obj%GaussPoint(1:3, 7) = [pval(1), pval(3), pval(1)]; obj%GaussIntegWei(7) = (wval(1)*wval(3)*wval(1)) obj%GaussPoint(1:3, 8) = [pval(2), pval(3), pval(1)]; obj%GaussIntegWei(8) = (wval(2)*wval(3)*wval(1)) obj%GaussPoint(1:3, 9) = [pval(3), pval(3), pval(1)]; obj%GaussIntegWei(9) = (wval(3)*wval(3)*wval(1)) obj%GaussPoint(1:3, 10) = [pval(1), pval(1), pval(2)]; obj%GaussIntegWei(10) = (wval(1)*wval(1)*wval(2)) obj%GaussPoint(1:3, 11) = [pval(2), pval(1), pval(2)]; obj%GaussIntegWei(11) = (wval(2)*wval(1)*wval(2)) obj%GaussPoint(1:3, 12) = [pval(3), pval(1), pval(2)]; obj%GaussIntegWei(12) = (wval(3)*wval(1)*wval(2)) obj%GaussPoint(1:3, 13) = [pval(1), pval(2), pval(2)]; obj%GaussIntegWei(13) = (wval(1)*wval(2)*wval(2)) obj%GaussPoint(1:3, 14) = [pval(2), pval(2), pval(2)]; obj%GaussIntegWei(14) = (wval(2)*wval(2)*wval(2)) obj%GaussPoint(1:3, 15) = [pval(3), pval(2), pval(2)]; obj%GaussIntegWei(15) = (wval(3)*wval(2)*wval(2)) obj%GaussPoint(1:3, 16) = [pval(1), pval(3), pval(2)]; obj%GaussIntegWei(16) = (wval(1)*wval(3)*wval(2)) obj%GaussPoint(1:3, 17) = [pval(2), pval(3), pval(2)]; obj%GaussIntegWei(17) = (wval(2)*wval(3)*wval(2)) obj%GaussPoint(1:3, 18) = [pval(3), pval(3), pval(2)]; obj%GaussIntegWei(18) = (wval(3)*wval(3)*wval(2)) obj%GaussPoint(1:3, 19) = [pval(1), pval(1), pval(3)]; obj%GaussIntegWei(19) = (wval(1)*wval(1)*wval(3)) obj%GaussPoint(1:3, 20) = [pval(2), pval(1), pval(3)]; obj%GaussIntegWei(20) = (wval(2)*wval(1)*wval(3)) obj%GaussPoint(1:3, 21) = [pval(3), pval(1), pval(3)]; obj%GaussIntegWei(21) = (wval(3)*wval(1)*wval(3)) obj%GaussPoint(1:3, 22) = [pval(1), pval(2), pval(3)]; obj%GaussIntegWei(22) = (wval(1)*wval(2)*wval(3)) obj%GaussPoint(1:3, 23) = [pval(2), pval(2), pval(3)]; obj%GaussIntegWei(23) = (wval(2)*wval(2)*wval(3)) obj%GaussPoint(1:3, 24) = [pval(3), pval(2), pval(3)]; obj%GaussIntegWei(24) = (wval(3)*wval(2)*wval(3)) obj%GaussPoint(1:3, 25) = [pval(1), pval(3), pval(3)]; obj%GaussIntegWei(25) = (wval(1)*wval(3)*wval(3)) obj%GaussPoint(1:3, 26) = [pval(2), pval(3), pval(3)]; obj%GaussIntegWei(26) = (wval(2)*wval(3)*wval(3)) obj%GaussPoint(1:3, 27) = [pval(3), pval(3), pval(3)]; obj%GaussIntegWei(27) = (wval(3)*wval(3)*wval(3)) elseif (size(obj%GaussPoint, 2) == 1) then !! Triangular or Tetrahedral if (allocated(obj%GaussPoint)) then deallocate (obj%GaussPoint) end if if (allocated(obj%GaussIntegWei)) then deallocate (obj%GaussIntegWei) end if obj%GaussPoint(1, 1) = 1.0d0 obj%GaussPoint(2, 1) = 1.0d0 obj%GaussPoint(3, 1) = 1.0d0 obj%GaussIntegWei(:) = 1.0d0 else print *, "ERROR :: ShapeFunctionClass/GetGaussPoint is not defined" end if end subroutine GetGaussPoint !! ###################################### !! ###################################### subroutine SetGaussPoint(obj) class(ShapeFunction_), intent(inout)::obj if (allocated(obj%gzi)) then if (size(obj%gzi, 1) /= obj%NumOfDim) then deallocate (obj%gzi) allocate (obj%gzi(obj%NumOfDim)) end if else allocate (obj%gzi(obj%NumOfDim)) end if if (obj%NumOfDim /= size(obj%GaussPoint, 1)) then print *, "ERROR::SetGaussPoint", obj%NumOfDim, size(obj%GaussPoint, 1) obj%ErrorMsg = "ERROR::SetGaussPoint" obj%ierr = 1 else obj%gzi(:) = obj%GaussPoint(:, obj%GpID) obj%ErrorMsg = "Succeed::SetGaussPoint" obj%ierr = 0 end if end subroutine SetGaussPoint !! ###################################### !! ###################################### subroutine GetShapeFunction(obj) class(ShapeFunction_), intent(inout)::obj if (allocated(obj%Nmat)) then if (size(obj%Nmat, 1) /= obj%NumOfNode) then deallocate (obj%Nmat) allocate (obj%Nmat(obj%NumOfNode)) end if else allocate (obj%Nmat(obj%NumOfNode)) end if if (obj%NumOfNode == 1) then !! ######################################################################### !! ####### ####### !! ####### + (1) ####### !! ####### ####### !! ####### ####### !! ######################################################################### print *, "ERROR::GetShapeFunction" obj%ErrorMsg = "ERROR::GetShapeFunction NumOfNode==1 " obj%ierr = 1 elseif (obj%NumOfNode == 2) then if (obj%NumOfDim == 1) then !! ######################################################################### !! ####### ####### !! ####### +-----------------------+ ####### !! ####### (1) (2) ####### !! ####### ####### !! ######################################################################### obj%Nmat(1) = 0.50d0*(1.0d0 - obj%gzi(1)) obj%Nmat(2) = 0.50d0*(-1.0d0 + obj%gzi(1)) obj%ErrorMsg = "Succeed::GetShapeFunction " obj%ierr = 0 else print *, "ERROR::GetShapeFunction" obj%ErrorMsg = "ERROR::GetShapeFunction NumOfNode==2 NumOfOrder/=1 " obj%ierr = 1 end if elseif (obj%NumOfNode == 3) then if (obj%NumOfDim == 1) then !! ######################################################################### !! ####### ####### !! ####### +-----------+-------------+ ####### !! ####### (1) (2) (3) ####### !! ####### ####### !! ######################################################################### !! allocate(obj%Nmat(obj%NumOfNode,obj%NumOfDim) ) !! !! obj%Nmat(1,1)=0.50d0*( 1.0d0-gzi(1)) !! obj%Nmat(1,2)=0.50d0*(-1.0d0+gzi(1)) !! obj%Nmat(1,3)=0.50d0*( 1.0d0-gzi(1)) !! !! obj%ErrorMsg="Succeed::GetShapeFunction " !! obj%ierr=0 print *, "ERROR::GetShapeFunction" obj%ErrorMsg = "ERROR::GetShapeFunction " obj%ierr = 1 elseif (obj%NumOfDim == 2) then !! ######################################################################### !! ####### ####### !! ####### (1) +-------------------------+ ####### !! ####### \ / (3) ####### !! ####### \ / ####### !! ####### \ / ####### !! ####### \ / ####### !! ####### \ / ####### !! ####### \ / ####### !! ####### (2) + ####### !! ######################################################################### print *, "ERROR::GetShapeFunction" obj%ErrorMsg = "ERROR::GetShapeFunction " obj%ierr = 1 else print *, "ERROR::GetShapeFunction" obj%ErrorMsg = "ERROR::GetShapeFunction " obj%ierr = 1 end if elseif (obj%NumOfNode == 4) then if (obj%NumOfDim == 1) then !! ######################################################################### !! ####### ####### !! ####### +-------+---------+-------+ ####### !! ####### (1) (2) (3) (4) ####### !! ####### ####### !! ######################################################################### !! obj%ErrorMsg="Succeed::GetShapeFunction " !! obj%ierr=0 print *, "ERROR::GetShapeFunction" obj%ErrorMsg = "ERROR::GetShapeFunction " obj%ierr = 1 elseif (obj%NumOfDim == 2) then !! ######################################################################### !! ####### ####### !! ####### (1) +-------------------------+ (4) ####### !! ####### ! ! ####### !! ####### ! ! ####### !! ####### ! ! ####### !! ####### ! ! ####### !! ####### ! ! ####### !! ####### ! ! ####### !! ####### (2) +-------------------------+ (3) ####### !! ######################################################################### obj%Nmat(1) = 0.250d0*(1.0d0 - obj%gzi(1))*(1.0d0 - obj%gzi(2)) obj%Nmat(2) = 0.250d0*(1.0d0 + obj%gzi(1))*(1.0d0 - obj%gzi(2)) obj%Nmat(3) = 0.250d0*(1.0d0 + obj%gzi(1))*(1.0d0 + obj%gzi(2)) obj%Nmat(4) = 0.250d0*(1.0d0 - obj%gzi(1))*(1.0d0 + obj%gzi(2)) obj%ErrorMsg = "Succeed::GetShapeFunction " obj%ierr = 0 elseif (obj%NumOfDim == 3) then !! ######################################################################### !! ####### (4) + ####### !! ####### / \ ####### !! ####### / \ ####### !! ####### / \ ####### !! ####### / ___+ (3) ####### !! ####### / ___---- \ ####### !! ####### (1) + - \ ####### !! ####### -----____ \ ####### !! ####### ----- + (2) ####### !! ######################################################################### print *, "ERROR::GetShapeFunction" obj%ErrorMsg = "ERROR::GetShapeFunction " obj%ierr = 1 else print *, "ERROR::GetShapeFunction" obj%ErrorMsg = "ERROR::GetShapeFunction " obj%ierr = 1 end if elseif (obj%NumOfNode == 5) then print *, "ERROR::GetShapeFunction" obj%ErrorMsg = "ERROR::GetShapeFunction " obj%ierr = 1 elseif (obj%NumOfNode == 6) then print *, "ERROR::GetShapeFunction" obj%ErrorMsg = "ERROR::GetShapeFunction " obj%ierr = 1 elseif (obj%NumOfNode == 7) then print *, "ERROR::GetShapeFunction" obj%ErrorMsg = "ERROR::GetShapeFunction " obj%ierr = 1 elseif (obj%NumOfNode == 8) then if (obj%NumOfDim == 3) then !! ######################################################################### !! ####### (8) +-------------------------+ (7) ####### !! ####### /! /! ####### !! ####### / ! (6) / ! ####### !! ####### (5) +--!----------------------+ ! ####### !! ####### ! ! ! ! ####### !! ####### ! +----------------------!--+ (3) ####### !! ####### ! / (4) ! / ####### !! ####### !/ !/ ####### !! ####### (1) +-------------------------+ (2) ####### !! ######################################################################### obj%Nmat(1) = 0.1250d0*(1.0d0 - obj%gzi(1))*(1.0d0 - obj%gzi(2))*(1.0d0 - obj%gzi(3)) obj%Nmat(2) = 0.1250d0*(1.0d0 + obj%gzi(1))*(1.0d0 - obj%gzi(2))*(1.0d0 - obj%gzi(3)) obj%Nmat(3) = 0.1250d0*(1.0d0 + obj%gzi(1))*(1.0d0 + obj%gzi(2))*(1.0d0 - obj%gzi(3)) obj%Nmat(4) = 0.1250d0*(1.0d0 - obj%gzi(1))*(1.0d0 + obj%gzi(2))*(1.0d0 - obj%gzi(3)) obj%Nmat(5) = 0.1250d0*(1.0d0 - obj%gzi(1))*(1.0d0 - obj%gzi(2))*(1.0d0 + obj%gzi(3)) obj%Nmat(6) = 0.1250d0*(1.0d0 + obj%gzi(1))*(1.0d0 - obj%gzi(2))*(1.0d0 + obj%gzi(3)) obj%Nmat(7) = 0.1250d0*(1.0d0 + obj%gzi(1))*(1.0d0 + obj%gzi(2))*(1.0d0 + obj%gzi(3)) obj%Nmat(8) = 0.1250d0*(1.0d0 - obj%gzi(1))*(1.0d0 + obj%gzi(2))*(1.0d0 + obj%gzi(3)) obj%ErrorMsg = "Succeed::GetShapeFunction " obj%ierr = 0 else print *, "ERROR::GetShapeFunction" obj%ErrorMsg = "ERROR::GetShapeFunction " obj%ierr = 1 end if elseif (obj%NumOfNode == 9) then print *, "ERROR::GetShapeFunction" obj%ErrorMsg = "ERROR::GetShapeFunction " obj%ierr = 1 elseif (obj%NumOfNode == 10) then print *, "ERROR::GetShapeFunction" obj%ErrorMsg = "ERROR::GetShapeFunction " obj%ierr = 1 elseif (obj%NumOfNode == 11) then print *, "ERROR::GetShapeFunction" obj%ErrorMsg = "ERROR::GetShapeFunction " obj%ierr = 1 elseif (obj%NumOfNode == 12) then print *, "ERROR::GetShapeFunction" obj%ErrorMsg = "ERROR::GetShapeFunction " obj%ierr = 1 elseif (obj%NumOfNode == 13) then print *, "ERROR::GetShapeFunction" obj%ErrorMsg = "ERROR::GetShapeFunction " obj%ierr = 1 elseif (obj%NumOfNode == 14) then print *, "ERROR::GetShapeFunction" obj%ErrorMsg = "ERROR::GetShapeFunction " obj%ierr = 1 elseif (obj%NumOfNode == 15) then print *, "ERROR::GetShapeFunction" obj%ErrorMsg = "ERROR::GetShapeFunction " obj%ierr = 1 elseif (obj%NumOfNode == 16) then print *, "ERROR::GetShapeFunction" obj%ErrorMsg = "ERROR::GetShapeFunction " obj%ierr = 1 elseif (obj%NumOfNode == 17) then print *, "ERROR::GetShapeFunction" obj%ErrorMsg = "ERROR::GetShapeFunction " obj%ierr = 1 elseif (obj%NumOfNode == 18) then print *, "ERROR::GetShapeFunction" obj%ErrorMsg = "ERROR::GetShapeFunction " obj%ierr = 1 elseif (obj%NumOfNode == 19) then print *, "ERROR::GetShapeFunction" obj%ErrorMsg = "ERROR::GetShapeFunction " obj%ierr = 1 elseif (obj%NumOfNode == 20) then if (obj%NumOfDim == 3) then !! ######################################################################### !! ####### (8) +----------(15)------------+ (7) ####### !! ####### (16) /! (14) /! ####### !! ####### / ! (20) (6) / ! (19) ####### !! ####### (5) +--!-------(13)------------+ ! ####### !! ####### ! ! ! ! ####### !! ####### (17)! +----------(11)--------!--+ (3) ####### !! ####### (12)! / (4) (18) ! / (10) ####### !! ####### !/ !/ ####### !! ####### (1) +----------(9)------------+ (2) ####### !! ######################################################################### ! Koers, R.W.J., 1989. Use of modified standard 20-node isoparametric brick elements for representing stress/strain fields at a crack tip for elastic and perfectly plastic material. Int. J. Fract. 40, 79–110. https://doi.org/10.1007/BF00963969 obj%Nmat(1) = -0.1250d0*(1.0d0 - obj%gzi(1))*(1.0d0 - obj%gzi(2))*(1.0d0 - obj%gzi(3)) & *(2.0d0 + obj%gzi(1) + obj%gzi(2) + obj%gzi(3)) obj%Nmat(2) = -0.1250d0*(1.0d0 + obj%gzi(1))*(1.0d0 - obj%gzi(2))*(1.0d0 - obj%gzi(3)) & *(2.0d0 - obj%gzi(1) + obj%gzi(2) + obj%gzi(3)) obj%Nmat(3) = -0.1250d0*(1.0d0 + obj%gzi(1))*(1.0d0 + obj%gzi(2))*(1.0d0 - obj%gzi(3)) & *(2.0d0 - obj%gzi(1) - obj%gzi(2) + obj%gzi(3)) obj%Nmat(4) = -0.1250d0*(1.0d0 - obj%gzi(1))*(1.0d0 + obj%gzi(2))*(1.0d0 - obj%gzi(3)) & *(2.0d0 + obj%gzi(1) - obj%gzi(2) + obj%gzi(3)) obj%Nmat(5) = -0.1250d0*(1.0d0 - obj%gzi(1))*(1.0d0 - obj%gzi(2))*(1.0d0 + obj%gzi(3)) & *(2.0d0 + obj%gzi(1) + obj%gzi(2) - obj%gzi(3)) obj%Nmat(6) = -0.1250d0*(1.0d0 + obj%gzi(1))*(1.0d0 - obj%gzi(2))*(1.0d0 + obj%gzi(3)) & *(2.0d0 - obj%gzi(1) + obj%gzi(2) - obj%gzi(3)) obj%Nmat(7) = -0.1250d0*(1.0d0 + obj%gzi(1))*(1.0d0 + obj%gzi(2))*(1.0d0 + obj%gzi(3)) & *(2.0d0 - obj%gzi(1) - obj%gzi(2) - obj%gzi(3)) obj%Nmat(8) = -0.1250d0*(1.0d0 - obj%gzi(1))*(1.0d0 + obj%gzi(2))*(1.0d0 + obj%gzi(3)) & *(2.0d0 + obj%gzi(1) - obj%gzi(2) - obj%gzi(3)) obj%Nmat(9) = 0.250d0*(1.0d0 - obj%gzi(1)**2)*(1.0d0 - obj%gzi(2))*(1.0d0 - obj%gzi(3)) obj%Nmat(10) = 0.250d0*(1.0d0 + obj%gzi(1))*(1.0d0 - obj%gzi(2)**2)*(1.0d0 - obj%gzi(3)) obj%Nmat(11) = 0.250d0*(1.0d0 - obj%gzi(1)**2)*(1.0d0 + obj%gzi(2))*(1.0d0 - obj%gzi(3)) obj%Nmat(12) = 0.250d0*(1.0d0 - obj%gzi(1))*(1.0d0 - obj%gzi(2)**2)*(1.0d0 - obj%gzi(3)) obj%Nmat(13) = 0.250d0*(1.0d0 - obj%gzi(1)**2)*(1.0d0 - obj%gzi(2))*(1.0d0 + obj%gzi(3)) obj%Nmat(14) = 0.250d0*(1.0d0 + obj%gzi(1))*(1.0d0 - obj%gzi(2)**2)*(1.0d0 + obj%gzi(3)) obj%Nmat(15) = 0.250d0*(1.0d0 - obj%gzi(1)**2)*(1.0d0 + obj%gzi(2))*(1.0d0 + obj%gzi(3)) obj%Nmat(16) = 0.250d0*(1.0d0 - obj%gzi(1))*(1.0d0 - obj%gzi(2)**2)*(1.0d0 + obj%gzi(3)) obj%Nmat(17) = 0.250d0*(1.0d0 - obj%gzi(1))*(1.0d0 - obj%gzi(2))*(1.0d0 - obj%gzi(3)**2) obj%Nmat(18) = 0.250d0*(1.0d0 + obj%gzi(1))*(1.0d0 - obj%gzi(2))*(1.0d0 - obj%gzi(3)**2) obj%Nmat(19) = 0.250d0*(1.0d0 + obj%gzi(1))*(1.0d0 + obj%gzi(2))*(1.0d0 - obj%gzi(3)**2) obj%Nmat(20) = 0.250d0*(1.0d0 - obj%gzi(1))*(1.0d0 + obj%gzi(2))*(1.0d0 - obj%gzi(3)**2) obj%ErrorMsg = "Succeed::GetShapeFunction " obj%ierr = 0 else print *, "ERROR::GetShapeFunction" obj%ErrorMsg = "ERROR::GetShapeFunction " obj%ierr = 1 end if elseif (obj%NumOfNode == 21) then print *, "ERROR::GetShapeFunction" obj%ErrorMsg = "ERROR::GetShapeFunction " obj%ierr = 1 elseif (obj%NumOfNode == 22) then print *, "ERROR::GetShapeFunction" obj%ErrorMsg = "ERROR::GetShapeFunction " obj%ierr = 1 elseif (obj%NumOfNode == 23) then print *, "ERROR::GetShapeFunction" obj%ErrorMsg = "ERROR::GetShapeFunction " obj%ierr = 1 elseif (obj%NumOfNode == 24) then print *, "ERROR::GetShapeFunction" obj%ErrorMsg = "ERROR::GetShapeFunction " obj%ierr = 1 else print *, "ERROR::GetShapeFunction" obj%ErrorMsg = "ERROR::GetShapeFunction " obj%ierr = 1 end if end subroutine GetShapeFunction !! ###################################### !! ###################################### subroutine GetShapeFuncDer1(obj) class(ShapeFunction_), intent(inout)::obj if (allocated(obj%dNdgzi)) then if (size(obj%dNdgzi, 1) /= obj%NumOfDim .or. size(obj%dNdgzi, 2) /= obj%NumOfNode) then deallocate (obj%dNdgzi) allocate (obj%dNdgzi(obj%NumOfDim, obj%NumOfNode)) end if else allocate (obj%dNdgzi(obj%NumOfDim, obj%NumOfNode)) end if if (obj%NumOfNode == 1) then !! ######################################################################### !! ####### ####### !! ####### + (1) ####### !! ####### ####### !! ####### ####### !! ######################################################################### print *, "ERROR::GetShapeFunction" obj%ErrorMsg = "ERROR::GetShapeFunction NumOfNode==1 " obj%ierr = 1 elseif (obj%NumOfNode == 2) then if (obj%NumOfDim == 1) then !! ######################################################################### !! ####### ####### !! ####### +-----------------------+ ####### !! ####### (1) (2) ####### !! ####### ####### !! ######################################################################### obj%dNdgzi(1, 1) = -0.50d0 obj%dNdgzi(1, 2) = 0.50d0 obj%ErrorMsg = "Succeed::GetShapeFunction " obj%ierr = 0 else print *, "ERROR::GetShapeFunction" obj%ErrorMsg = "ERROR::GetShapeFunction NumOfNode==2 NumOfOrder/=1 " obj%ierr = 1 end if elseif (obj%NumOfNode == 3) then if (obj%NumOfDim == 1) then !! ######################################################################### !! ####### ####### !! ####### +-----------+-------------+ ####### !! ####### (1) (2) (3) ####### !! ####### ####### !! ######################################################################### !! allocate(obj%dNdgzi(obj%NumOfNode,obj%NumOfDim) ) !! !! obj%dNdgzi(1,1)=0.50d0*( 1.0d0-gzi(1)) !! obj%dNdgzi(1,2)=0.50d0*(-1.0d0+gzi(1)) !! obj%dNdgzi(1,3)=0.50d0*( 1.0d0-gzi(1)) !! !! obj%ErrorMsg="Succeed::GetShapeFunction " !! obj%ierr=0 print *, "ERROR::GetShapeFunction" obj%ErrorMsg = "ERROR::GetShapeFunction " obj%ierr = 1 elseif (obj%NumOfDim == 2) then !! ######################################################################### !! ####### ####### !! ####### (1) +-------------------------+ ####### !! ####### \ / (3) ####### !! ####### \ / ####### !! ####### \ / ####### !! ####### \ / ####### !! ####### \ / ####### !! ####### \ / ####### !! ####### (2) + ####### !! ######################################################################### print *, "ERROR::GetShapeFunction" obj%ErrorMsg = "ERROR::GetShapeFunction " obj%ierr = 1 else print *, "ERROR::GetShapeFunction" obj%ErrorMsg = "ERROR::GetShapeFunction " obj%ierr = 1 end if elseif (obj%NumOfNode == 4) then if (obj%NumOfDim == 1) then !! ######################################################################### !! ####### ####### !! ####### +-------+---------+-------+ ####### !! ####### (1) (2) (3) (4) ####### !! ####### ####### !! ######################################################################### !! obj%ErrorMsg="Succeed::GetShapeFunction " !! obj%ierr=0 print *, "ERROR::GetShapeFunction" obj%ErrorMsg = "ERROR::GetShapeFunction " obj%ierr = 1 elseif (obj%NumOfDim == 2) then !! ######################################################################### !! ####### ####### !! ####### (4) +-------------------------+ (3) ####### !! ####### ! ! ####### !! ####### ! ! ####### !! ####### ! ! ####### !! ####### ! ! ####### !! ####### ! ! ####### !! ####### ! ! ####### !! ####### (1) +-------------------------+ (2) ####### !! ######################################################################### obj%dNdgzi(1, 1) = -0.250d0*(1.0d0 - obj%gzi(2)) obj%dNdgzi(1, 2) = 0.250d0*(1.0d0 - obj%gzi(2)) obj%dNdgzi(1, 3) = 0.250d0*(1.0d0 + obj%gzi(2)) obj%dNdgzi(1, 4) = -0.250d0*(1.0d0 + obj%gzi(2)) obj%dNdgzi(2, 1) = -0.250d0*(1.0d0 - obj%gzi(1)) obj%dNdgzi(2, 2) = -0.250d0*(1.0d0 + obj%gzi(1)) obj%dNdgzi(2, 3) = 0.250d0*(1.0d0 + obj%gzi(1)) obj%dNdgzi(2, 4) = 0.250d0*(1.0d0 - obj%gzi(1)) obj%ErrorMsg = "Succeed::GetShapeFunction " obj%ierr = 0 elseif (obj%NumOfDim == 3) then !! ######################################################################### !! ####### (4) + ####### !! ####### / \ ####### !! ####### / \ ####### !! ####### / \ ####### !! ####### / ___+ (3) ####### !! ####### / ___---- \ ####### !! ####### (1) + - \ ####### !! ####### -----____ \ ####### !! ####### ----- + (2) ####### !! ######################################################################### print *, "ERROR::GetShapeFunction" obj%ErrorMsg = "ERROR::GetShapeFunction " obj%ierr = 1 else print *, "ERROR::GetShapeFunction" obj%ErrorMsg = "ERROR::GetShapeFunction " obj%ierr = 1 end if elseif (obj%NumOfNode == 5) then print *, "ERROR::GetShapeFunction" obj%ErrorMsg = "ERROR::GetShapeFunction " obj%ierr = 1 elseif (obj%NumOfNode == 6) then print *, "ERROR::GetShapeFunction" obj%ErrorMsg = "ERROR::GetShapeFunction " obj%ierr = 1 elseif (obj%NumOfNode == 7) then print *, "ERROR::GetShapeFunction" obj%ErrorMsg = "ERROR::GetShapeFunction " obj%ierr = 1 elseif (obj%NumOfNode == 8) then if (obj%NumOfDim == 3) then !! ######################################################################### !! ####### (8) +-------------------------+ (7) ####### !! ####### /! /! ####### !! ####### / ! (6) / ! ####### !! ####### (5) +--!----------------------+ ! ####### !! ####### ! ! ! ! ####### !! ####### ! +----------------------!--+ (3) ####### !! ####### ! / (4) ! / ####### !! ####### !/ !/ ####### !! ####### (1) +-------------------------+ (2) ####### !! ######################################################################### obj%dNdgzi(1, 1) = -0.1250d0*(1.0d0 - obj%gzi(2))*(1.0d0 - obj%gzi(3)) obj%dNdgzi(1, 2) = +0.1250d0*(1.0d0 - obj%gzi(2))*(1.0d0 - obj%gzi(3)) obj%dNdgzi(1, 3) = +0.1250d0*(1.0d0 + obj%gzi(2))*(1.0d0 - obj%gzi(3)) obj%dNdgzi(1, 4) = -0.1250d0*(1.0d0 + obj%gzi(2))*(1.0d0 - obj%gzi(3)) obj%dNdgzi(1, 5) = -0.1250d0*(1.0d0 - obj%gzi(2))*(1.0d0 + obj%gzi(3)) obj%dNdgzi(1, 6) = +0.1250d0*(1.0d0 - obj%gzi(2))*(1.0d0 + obj%gzi(3)) obj%dNdgzi(1, 7) = +0.1250d0*(1.0d0 + obj%gzi(2))*(1.0d0 + obj%gzi(3)) obj%dNdgzi(1, 8) = -0.1250d0*(1.0d0 + obj%gzi(2))*(1.0d0 + obj%gzi(3)) obj%dNdgzi(2, 1) = -0.1250d0*(1.0d0 - obj%gzi(1))*(1.0d0 - obj%gzi(3)) obj%dNdgzi(2, 2) = -0.1250d0*(1.0d0 + obj%gzi(1))*(1.0d0 - obj%gzi(3)) obj%dNdgzi(2, 3) = +0.1250d0*(1.0d0 + obj%gzi(1))*(1.0d0 - obj%gzi(3)) obj%dNdgzi(2, 4) = +0.1250d0*(1.0d0 - obj%gzi(1))*(1.0d0 - obj%gzi(3)) obj%dNdgzi(2, 5) = -0.1250d0*(1.0d0 - obj%gzi(1))*(1.0d0 + obj%gzi(3)) obj%dNdgzi(2, 6) = -0.1250d0*(1.0d0 + obj%gzi(1))*(1.0d0 + obj%gzi(3)) obj%dNdgzi(2, 7) = +0.1250d0*(1.0d0 + obj%gzi(1))*(1.0d0 + obj%gzi(3)) obj%dNdgzi(2, 8) = +0.1250d0*(1.0d0 - obj%gzi(1))*(1.0d0 + obj%gzi(3)) obj%dNdgzi(3, 1) = -0.1250d0*(1.0d0 - obj%gzi(1))*(1.0d0 - obj%gzi(2)) obj%dNdgzi(3, 2) = -0.1250d0*(1.0d0 + obj%gzi(1))*(1.0d0 - obj%gzi(2)) obj%dNdgzi(3, 3) = -0.1250d0*(1.0d0 + obj%gzi(1))*(1.0d0 + obj%gzi(2)) obj%dNdgzi(3, 4) = -0.1250d0*(1.0d0 - obj%gzi(1))*(1.0d0 + obj%gzi(2)) obj%dNdgzi(3, 5) = +0.1250d0*(1.0d0 - obj%gzi(1))*(1.0d0 - obj%gzi(2)) obj%dNdgzi(3, 6) = +0.1250d0*(1.0d0 + obj%gzi(1))*(1.0d0 - obj%gzi(2)) obj%dNdgzi(3, 7) = +0.1250d0*(1.0d0 + obj%gzi(1))*(1.0d0 + obj%gzi(2)) obj%dNdgzi(3, 8) = +0.1250d0*(1.0d0 - obj%gzi(1))*(1.0d0 + obj%gzi(2)) obj%ErrorMsg = "Succeed::GetShapeFunction " obj%ierr = 0 else print *, "ERROR::GetShapeFunction" obj%ErrorMsg = "ERROR::GetShapeFunction " obj%ierr = 1 end if elseif (obj%NumOfNode == 9) then print *, "ERROR::GetShapeFunction" obj%ErrorMsg = "ERROR::GetShapeFunction " obj%ierr = 1 elseif (obj%NumOfNode == 10) then print *, "ERROR::GetShapeFunction" obj%ErrorMsg = "ERROR::GetShapeFunction " obj%ierr = 1 elseif (obj%NumOfNode == 11) then print *, "ERROR::GetShapeFunction" obj%ErrorMsg = "ERROR::GetShapeFunction " obj%ierr = 1 elseif (obj%NumOfNode == 12) then print *, "ERROR::GetShapeFunction" obj%ErrorMsg = "ERROR::GetShapeFunction " obj%ierr = 1 elseif (obj%NumOfNode == 13) then print *, "ERROR::GetShapeFunction" obj%ErrorMsg = "ERROR::GetShapeFunction " obj%ierr = 1 elseif (obj%NumOfNode == 14) then print *, "ERROR::GetShapeFunction" obj%ErrorMsg = "ERROR::GetShapeFunction " obj%ierr = 1 elseif (obj%NumOfNode == 15) then print *, "ERROR::GetShapeFunction" obj%ErrorMsg = "ERROR::GetShapeFunction " obj%ierr = 1 elseif (obj%NumOfNode == 16) then print *, "ERROR::GetShapeFunction" obj%ErrorMsg = "ERROR::GetShapeFunction " obj%ierr = 1 elseif (obj%NumOfNode == 17) then print *, "ERROR::GetShapeFunction" obj%ErrorMsg = "ERROR::GetShapeFunction " obj%ierr = 1 elseif (obj%NumOfNode == 18) then print *, "ERROR::GetShapeFunction" obj%ErrorMsg = "ERROR::GetShapeFunction " obj%ierr = 1 elseif (obj%NumOfNode == 19) then print *, "ERROR::GetShapeFunction" obj%ErrorMsg = "ERROR::GetShapeFunction " obj%ierr = 1 elseif (obj%NumOfNode == 20) then if (obj%NumOfDim == 3) then !! ######################################################################### !! ####### (8) +----------(15)------------+ (7) ####### !! ####### (16) /! (14) /! ####### !! ####### / ! (20) (6) / ! (19) ####### !! ####### (5) +--!-------(13)------------+ ! ####### !! ####### ! ! ! ! ####### !! ####### (17)! +----------(11)--------!--+ (3) ####### !! ####### (12)! / (4) (18) ! / (10) ####### !! ####### !/ !/ ####### !! ####### (1) +----------(9)------------+ (2) ####### !! ######################################################################### ! Koers, R.W.J., 1989. Use of modified standard 20-node isoparametric brick elements for representing stress/strain fields at a crack tip for elastic and perfectly plastic material. Int. J. Fract. 40, 79–110. https://doi.org/10.1007/BF00963969 obj%dNdgzi(1, 1) = -0.1250d0*(-1.0d0)*(1.0d0 - obj%gzi(2))*(1.0d0 - obj%gzi(3)) & *(2.0d0 + obj%gzi(1) + obj%gzi(2) + obj%gzi(3)) & - 0.1250d0*(1.0d0 - obj%gzi(1))*(1.0d0 - obj%gzi(2))*(1.0d0 - obj%gzi(3)) & *1.0d0 obj%dNdgzi(1, 2) = -0.1250d0*(1.0d0)*(1.0d0 - obj%gzi(2))*(1.0d0 - obj%gzi(3)) & *(2.0d0 - obj%gzi(1) + obj%gzi(2) + obj%gzi(3)) & - 0.1250d0*(1.0d0 + obj%gzi(1))*(1.0d0 - obj%gzi(2))*(1.0d0 - obj%gzi(3)) & *(-1.0d0) obj%dNdgzi(1, 3) = -0.1250d0*(1.0d0)*(1.0d0 + obj%gzi(2))*(1.0d0 - obj%gzi(3)) & *(2.0d0 - obj%gzi(1) - obj%gzi(2) + obj%gzi(3)) & - 0.1250d0*(1.0d0 + obj%gzi(1))*(1.0d0 + obj%gzi(2))*(1.0d0 - obj%gzi(3)) & *(-1.0d0) obj%dNdgzi(1, 4) = -0.1250d0*(-1.0d0)*(1.0d0 + obj%gzi(2))*(1.0d0 - obj%gzi(3)) & *(2.0d0 + obj%gzi(1) - obj%gzi(2) + obj%gzi(3)) & - 0.1250d0*(1.0d0 - obj%gzi(1))*(1.0d0 + obj%gzi(2))*(1.0d0 - obj%gzi(3)) & *(1.0d0) obj%dNdgzi(1, 5) = -0.1250d0*(-1.0d0)*(1.0d0 - obj%gzi(2))*(1.0d0 + obj%gzi(3)) & *(2.0d0 + obj%gzi(1) + obj%gzi(2) - obj%gzi(3)) & - 0.1250d0*(1.0d0 - obj%gzi(1))*(1.0d0 - obj%gzi(2))*(1.0d0 + obj%gzi(3)) & *(1.0d0) obj%dNdgzi(1, 6) = -0.1250d0*(1.0d0)*(1.0d0 - obj%gzi(2))*(1.0d0 + obj%gzi(3)) & *(2.0d0 - obj%gzi(1) + obj%gzi(2) - obj%gzi(3)) & - 0.1250d0*(1.0d0 + obj%gzi(1))*(1.0d0 - obj%gzi(2))*(1.0d0 + obj%gzi(3)) & *(-1.0d0) obj%dNdgzi(1, 7) = -0.1250d0*(1.0d0)*(1.0d0 + obj%gzi(2))*(1.0d0 + obj%gzi(3)) & *(2.0d0 - obj%gzi(1) - obj%gzi(2) - obj%gzi(3)) & - 0.1250d0*(1.0d0 + obj%gzi(1))*(1.0d0 + obj%gzi(2))*(1.0d0 + obj%gzi(3)) & *(-1.0d0) obj%dNdgzi(1, 8) = -0.1250d0*(-1.0d0)*(1.0d0 + obj%gzi(2))*(1.0d0 + obj%gzi(3)) & *(2.0d0 + obj%gzi(1) - obj%gzi(2) - obj%gzi(3)) & - 0.1250d0*(1.0d0 - obj%gzi(1))*(1.0d0 + obj%gzi(2))*(1.0d0 + obj%gzi(3)) & *(1.0d0) !<<<ok>>> obj%dNdgzi(1, 9) = 0.250d0*(-2.0d0*obj%gzi(1))*(1.0d0 - obj%gzi(2))*(1.0d0 - obj%gzi(3)) obj%dNdgzi(1, 10) = 0.250d0*(1.0d0)*(1.0d0 - obj%gzi(2)**2)*(1.0d0 - obj%gzi(3)) obj%dNdgzi(1, 11) = 0.250d0*(-2.0d0*obj%gzi(1))*(1.0d0 + obj%gzi(2))*(1.0d0 - obj%gzi(3)) obj%dNdgzi(1, 12) = 0.250d0*(-1.0d0)*(1.0d0 - obj%gzi(2)**2)*(1.0d0 - obj%gzi(3)) obj%dNdgzi(1, 13) = 0.250d0*(-2.0d0*obj%gzi(1))*(1.0d0 - obj%gzi(2))*(1.0d0 + obj%gzi(3)) obj%dNdgzi(1, 14) = 0.250d0*(1.0d0)*(1.0d0 - obj%gzi(2)**2)*(1.0d0 + obj%gzi(3)) obj%dNdgzi(1, 15) = 0.250d0*(-2.0d0*obj%gzi(1))*(1.0d0 + obj%gzi(2))*(1.0d0 + obj%gzi(3)) obj%dNdgzi(1, 16) = 0.250d0*(-1.0d0)*(1.0d0 - obj%gzi(2)**2)*(1.0d0 + obj%gzi(3)) obj%dNdgzi(1, 17) = 0.250d0*(-1.0d0)*(1.0d0 - obj%gzi(2))*(1.0d0 - obj%gzi(3)**2) obj%dNdgzi(1, 18) = 0.250d0*(1.0d0)*(1.0d0 - obj%gzi(2))*(1.0d0 - obj%gzi(3)**2) obj%dNdgzi(1, 19) = 0.250d0*(1.0d0)*(1.0d0 + obj%gzi(2))*(1.0d0 - obj%gzi(3)**2) obj%dNdgzi(1, 20) = 0.250d0*(-1.0d0)*(1.0d0 + obj%gzi(2))*(1.0d0 - obj%gzi(3)**2) !<<<ok>>> obj%dNdgzi(2, 1) = -0.1250d0*(1.0d0 - obj%gzi(1))*(-1.0d0)*(1.0d0 - obj%gzi(3)) & *(2.0d0 + obj%gzi(1) + obj%gzi(2) + obj%gzi(3)) & - 0.1250d0*(1.0d0 - obj%gzi(1))*(1.0d0 - obj%gzi(2))*(1.0d0 - obj%gzi(3)) & *(1.0d0) obj%dNdgzi(2, 2) = -0.1250d0*(1.0d0 + obj%gzi(1))*(-1.0d0)*(1.0d0 - obj%gzi(3)) & *(2.0d0 - obj%gzi(1) + obj%gzi(2) + obj%gzi(3)) & - 0.1250d0*(1.0d0 + obj%gzi(1))*(1.0d0 - obj%gzi(2))*(1.0d0 - obj%gzi(3)) & *(1.0d0) obj%dNdgzi(2, 3) = -0.1250d0*(1.0d0 + obj%gzi(1))*(1.0d0)*(1.0d0 - obj%gzi(3)) & *(2.0d0 - obj%gzi(1) - obj%gzi(2) + obj%gzi(3)) & - 0.1250d0*(1.0d0 + obj%gzi(1))*(1.0d0 + obj%gzi(2))*(1.0d0 - obj%gzi(3)) & *(-1.0d0) obj%dNdgzi(2, 4) = -0.1250d0*(1.0d0 - obj%gzi(1))*(1.0d0)*(1.0d0 - obj%gzi(3)) & *(2.0d0 + obj%gzi(1) - obj%gzi(2) + obj%gzi(3)) & - 0.1250d0*(1.0d0 - obj%gzi(1))*(1.0d0 + obj%gzi(2))*(1.0d0 - obj%gzi(3)) & *(-1.0d0) obj%dNdgzi(2, 5) = -0.1250d0*(1.0d0 - obj%gzi(1))*(-1.0d0)*(1.0d0 + obj%gzi(3)) & *(2.0d0 + obj%gzi(1) + obj%gzi(2) - obj%gzi(3)) & - 0.1250d0*(1.0d0 - obj%gzi(1))*(1.0d0 - obj%gzi(2))*(1.0d0 + obj%gzi(3)) & *(1.0d0) obj%dNdgzi(2, 6) = -0.1250d0*(1.0d0 + obj%gzi(1))*(-1.0d0)*(1.0d0 + obj%gzi(3)) & *(2.0d0 - obj%gzi(1) + obj%gzi(2) - obj%gzi(3)) & - 0.1250d0*(1.0d0 + obj%gzi(1))*(1.0d0 - obj%gzi(2))*(1.0d0 + obj%gzi(3)) & *(1.0d0) obj%dNdgzi(2, 7) = -0.1250d0*(1.0d0 + obj%gzi(1))*(1.0d0)*(1.0d0 + obj%gzi(3)) & *(2.0d0 - obj%gzi(1) - obj%gzi(2) - obj%gzi(3)) & - 0.1250d0*(1.0d0 + obj%gzi(1))*(1.0d0 + obj%gzi(2))*(1.0d0 + obj%gzi(3)) & *(-1.0d0) obj%dNdgzi(2, 8) = -0.1250d0*(1.0d0 - obj%gzi(1))*(1.0d0)*(1.0d0 + obj%gzi(3)) & *(2.0d0 + obj%gzi(1) - obj%gzi(2) - obj%gzi(3)) & - 0.1250d0*(1.0d0 - obj%gzi(1))*(1.0d0 + obj%gzi(2))*(1.0d0 + obj%gzi(3)) & *(-1.0d0) !<<<ok>>> obj%dNdgzi(2, 9) = 0.250d0*(1.0d0 - obj%gzi(1)**2)*(-1.0d0)*(1.0d0 - obj%gzi(3)) obj%dNdgzi(2, 10) = 0.250d0*(1.0d0 + obj%gzi(1))*(-2.0d0*obj%gzi(2))*(1.0d0 - obj%gzi(3)) obj%dNdgzi(2, 11) = 0.250d0*(1.0d0 - obj%gzi(1)**2)*(1.0d0)*(1.0d0 - obj%gzi(3)) obj%dNdgzi(2, 12) = 0.250d0*(1.0d0 - obj%gzi(1))*(-2.0d0*obj%gzi(2))*(1.0d0 - obj%gzi(3)) obj%dNdgzi(2, 13) = 0.250d0*(1.0d0 - obj%gzi(1)**2)*(-1.0d0)*(1.0d0 + obj%gzi(3)) obj%dNdgzi(2, 14) = 0.250d0*(1.0d0 + obj%gzi(1))*(-2.0d0*obj%gzi(2))*(1.0d0 + obj%gzi(3)) obj%dNdgzi(2, 15) = 0.250d0*(1.0d0 - obj%gzi(1)**2)*(1.0d0)*(1.0d0 + obj%gzi(3)) obj%dNdgzi(2, 16) = 0.250d0*(1.0d0 - obj%gzi(1))*(-2.0d0*obj%gzi(2))*(1.0d0 + obj%gzi(3)) obj%dNdgzi(2, 17) = 0.250d0*(1.0d0 - obj%gzi(1))*(-1.0d0)*(1.0d0 - obj%gzi(3)**2) obj%dNdgzi(2, 18) = 0.250d0*(1.0d0 + obj%gzi(1))*(-1.0d0)*(1.0d0 - obj%gzi(3)**2) obj%dNdgzi(2, 19) = 0.250d0*(1.0d0 + obj%gzi(1))*(1.0d0)*(1.0d0 - obj%gzi(3)**2) obj%dNdgzi(2, 20) = 0.250d0*(1.0d0 - obj%gzi(1))*(1.0d0)*(1.0d0 - obj%gzi(3)**2) !<<<ok>>> obj%dNdgzi(3, 1) = -0.1250d0*(1.0d0 - obj%gzi(1))*(1.0d0 - obj%gzi(2))*(-1.0d0) & *(2.0d0 + obj%gzi(1) + obj%gzi(2) + obj%gzi(3)) & - 0.1250d0*(1.0d0 - obj%gzi(1))*(1.0d0 - obj%gzi(2))*(1.0d0 - obj%gzi(3)) & *(1.0d0) obj%dNdgzi(3, 2) = -0.1250d0*(1.0d0 + obj%gzi(1))*(1.0d0 - obj%gzi(2))*(-1.0d0) & *(2.0d0 - obj%gzi(1) + obj%gzi(2) + obj%gzi(3)) & - 0.1250d0*(1.0d0 + obj%gzi(1))*(1.0d0 - obj%gzi(2))*(1.0d0 - obj%gzi(3)) & *(1.0d0) obj%dNdgzi(3, 3) = -0.1250d0*(1.0d0 + obj%gzi(1))*(1.0d0 + obj%gzi(2))*(-1.0d0) & *(2.0d0 - obj%gzi(1) - obj%gzi(2) + obj%gzi(3)) & - 0.1250d0*(1.0d0 + obj%gzi(1))*(1.0d0 + obj%gzi(2))*(1.0d0 - obj%gzi(3)) & *(1.0d0) obj%dNdgzi(3, 4) = -0.1250d0*(1.0d0 - obj%gzi(1))*(1.0d0 + obj%gzi(2))*(-1.0d0) & *(2.0d0 + obj%gzi(1) - obj%gzi(2) + obj%gzi(3)) & - 0.1250d0*(1.0d0 - obj%gzi(1))*(1.0d0 + obj%gzi(2))*(1.0d0 - obj%gzi(3)) & *(1.0d0) obj%dNdgzi(3, 5) = -0.1250d0*(1.0d0 - obj%gzi(1))*(1.0d0 - obj%gzi(2))*(1.0d0) & *(2.0d0 + obj%gzi(1) + obj%gzi(2) - obj%gzi(3)) & - 0.1250d0*(1.0d0 - obj%gzi(1))*(1.0d0 - obj%gzi(2))*(1.0d0 + obj%gzi(3)) & *(-1.0d0) obj%dNdgzi(3, 6) = -0.1250d0*(1.0d0 + obj%gzi(1))*(1.0d0 - obj%gzi(2))*(1.0d0) & *(2.0d0 - obj%gzi(1) + obj%gzi(2) - obj%gzi(3)) & - 0.1250d0*(1.0d0 + obj%gzi(1))*(1.0d0 - obj%gzi(2))*(1.0d0 + obj%gzi(3)) & *(-1.0d0) obj%dNdgzi(3, 7) = -0.1250d0*(1.0d0 + obj%gzi(1))*(1.0d0 + obj%gzi(2))*(1.0d0) & *(2.0d0 - obj%gzi(1) - obj%gzi(2) - obj%gzi(3)) & - 0.1250d0*(1.0d0 + obj%gzi(1))*(1.0d0 + obj%gzi(2))*(1.0d0 + obj%gzi(3)) & *(-1.0d0) obj%dNdgzi(3, 8) = -0.1250d0*(1.0d0 - obj%gzi(1))*(1.0d0 + obj%gzi(2))*(1.0d0) & *(2.0d0 + obj%gzi(1) - obj%gzi(2) - obj%gzi(3)) & - 0.1250d0*(1.0d0 - obj%gzi(1))*(1.0d0 + obj%gzi(2))*(1.0d0 + obj%gzi(3)) & *(-1.0d0) !<<<ok>>> obj%dNdgzi(3, 9) = 0.250d0*(1.0d0 - obj%gzi(1)**2)*(1.0d0 - obj%gzi(2))*(-1.0d0) obj%dNdgzi(3, 10) = 0.250d0*(1.0d0 + obj%gzi(1))*(1.0d0 - obj%gzi(2)**2)*(-1.0d0) obj%dNdgzi(3, 11) = 0.250d0*(1.0d0 - obj%gzi(1)**2)*(1.0d0 + obj%gzi(2))*(-1.0d0) obj%dNdgzi(3, 12) = 0.250d0*(1.0d0 - obj%gzi(1))*(1.0d0 - obj%gzi(2)**2)*(-1.0d0) obj%dNdgzi(3, 13) = 0.250d0*(1.0d0 - obj%gzi(1)**2)*(1.0d0 - obj%gzi(2))*(1.0d0) obj%dNdgzi(3, 14) = 0.250d0*(1.0d0 + obj%gzi(1))*(1.0d0 - obj%gzi(2)**2)*(1.0d0) obj%dNdgzi(3, 15) = 0.250d0*(1.0d0 - obj%gzi(1)**2)*(1.0d0 + obj%gzi(2))*(1.0d0) obj%dNdgzi(3, 16) = 0.250d0*(1.0d0 - obj%gzi(1))*(1.0d0 - obj%gzi(2)**2)*(1.0d0) obj%dNdgzi(3, 17) = 0.250d0*(1.0d0 - obj%gzi(1))*(1.0d0 - obj%gzi(2))*(-2.0d0*obj%gzi(3)) obj%dNdgzi(3, 18) = 0.250d0*(1.0d0 + obj%gzi(1))*(1.0d0 - obj%gzi(2))*(-2.0d0*obj%gzi(3)) obj%dNdgzi(3, 19) = 0.250d0*(1.0d0 + obj%gzi(1))*(1.0d0 + obj%gzi(2))*(-2.0d0*obj%gzi(3)) obj%dNdgzi(3, 20) = 0.250d0*(1.0d0 - obj%gzi(1))*(1.0d0 + obj%gzi(2))*(-2.0d0*obj%gzi(3)) obj%ErrorMsg = "Succeed::GetShapeFunction " obj%ierr = 0 else print *, "ERROR::GetShapeFunction" obj%ErrorMsg = "ERROR::GetShapeFunction " obj%ierr = 1 end if elseif (obj%NumOfNode == 21) then print *, "ERROR::GetShapeFunction" obj%ErrorMsg = "ERROR::GetShapeFunction " obj%ierr = 1 elseif (obj%NumOfNode == 22) then print *, "ERROR::GetShapeFunction" obj%ErrorMsg = "ERROR::GetShapeFunction " obj%ierr = 1 elseif (obj%NumOfNode == 23) then print *, "ERROR::GetShapeFunction" obj%ErrorMsg = "ERROR::GetShapeFunction " obj%ierr = 1 elseif (obj%NumOfNode == 24) then print *, "ERROR::GetShapeFunction" obj%ErrorMsg = "ERROR::GetShapeFunction " obj%ierr = 1 else print *, "ERROR::GetShapeFunction" obj%ErrorMsg = "ERROR::GetShapeFunction " obj%ierr = 1 end if end subroutine GetShapeFuncDer1 !! ###################################### !! ###################################### subroutine GetShapeFuncDer2(obj) class(ShapeFunction_), intent(inout)::obj end subroutine GetShapeFuncDer2 !! ###################################### !! ###################################### subroutine GetElemCoord(obj, nod_coord, elem_nod, elem_id) class(ShapeFunction_), intent(inout)::obj integer(int32), intent(in)::elem_nod(:, :), elem_id real(real64), intent(in)::nod_coord(:, :) integer(int32)::i, j, k, n, m n = size(elem_nod, 2) m = size(nod_coord, 2) if (allocated(obj%ElemCoord)) then if (n /= size(obj%ElemCoord, 1) .or. m /= size(obj%ElemCoord, 2)) then deallocate (obj%ElemCoord) allocate (obj%ElemCoord(n, m)) end if else allocate (obj%ElemCoord(n, m)) end if do j = 1, n obj%ElemCoord(j, 1:m) = nod_coord(elem_nod(elem_id, j), 1:m) !! print *, obj%ElemCoord(j,1:m) end do return end subroutine GetElemCoord !! ###################################### !! ###################################### subroutine GetElemCoord_n(obj, nod_coord_n, elem_nod, elem_id) class(ShapeFunction_), intent(inout)::obj integer(int32), intent(in)::elem_nod(:, :), elem_id real(real64), intent(in)::nod_coord_n(:, :) integer(int32)::i, j, k, n, m n = size(elem_nod, 2) m = size(nod_coord_n, 2) if (allocated(obj%ElemCoord_n)) then if (n /= size(obj%ElemCoord_n, 1) .or. m /= size(obj%ElemCoord_n, 2)) then deallocate (obj%ElemCoord_n) allocate (obj%ElemCoord_n(n, m)) end if else allocate (obj%ElemCoord_n(n, m)) end if do j = 1, n obj%ElemCoord_n(j, 1:m) = nod_coord_n(elem_nod(elem_id, j), 1:m) !! print *, obj%ElemCoord_n(j,1:m) end do return end subroutine !! ###################################### !! ###################################### subroutine getdu(obj) class(ShapeFunction_), intent(inout)::obj integer(int32)::i, j, k, n, m n = size(obj%ElemCoord, 1) m = size(obj%ElemCoord, 2) if (allocated(obj%du)) then if (n /= size(obj%du, 1) .or. m /= size(obj%du, 2)) then deallocate (obj%du) allocate (obj%du(n, m)) end if else allocate (obj%du(n, m)) end if do j = 1, n obj%du(j, 1:m) = obj%ElemCoord(j, 1:m) - obj%ElemCoord_n(j, 1:m) end do return end subroutine !! ###################################### !! ###################################### subroutine GetJmat(obj) class(ShapeFunction_), intent(inout)::obj integer(int32) n n = size(obj%ElemCoord, 2) if (allocated(obj%Jmat)) then if (n /= size(obj%Jmat, 1) .or. n /= size(obj%Jmat, 2)) then deallocate (obj%Jmat) allocate (obj%Jmat(n, n)) allocate (obj%JmatInv(n, n)) end if else allocate (obj%Jmat(n, n)) allocate (obj%JmatInv(n, n)) end if obj%Jmat(:, :) = matmul(obj%dNdgzi, obj%ElemCoord) call inverse_rank_2(obj%Jmat, obj%JmatInv) end subroutine GetJmat !! ###################################### subroutine exportShapeFunction(obj, restart, path) class(ShapeFunction_), intent(inout) :: obj logical, optional, intent(in) :: restart character(*), intent(in) :: path type(IO_) :: f if (present(restart)) then call execute_command_line("mkdir -p "//path//"/ShapeFunction") call f%open(path//"/ShapeFunction/", "ShapeFunction", ".res") write (f%fh, *) obj%Nmat(:) write (f%fh, *) obj%dNdgzi(:, :) write (f%fh, *) obj%dNdgzidgzi(:, :) write (f%fh, *) obj%gzi(:) write (f%fh, *) obj%gaussPoint(:, :) write (f%fh, *) obj%gaussIntegWei(:) write (f%fh, *) obj%Jmat(:, :) write (f%fh, *) obj%JmatInv(:, :) write (f%fh, *) obj%ElemCoord(:, :) write (f%fh, *) obj%ElemCoord_n(:, :) write (f%fh, *) obj%du(:, :) write (f%fh, *) obj%detJ write (f%fh, *) obj%NumOfNode write (f%fh, *) obj%NumOfOrder write (f%fh, *) obj%NumOfDim write (f%fh, *) obj%NumOfGp write (f%fh, *) obj%GpID write (f%fh, *) obj%ierr write (f%fh, *) obj%currentGpID write (f%fh, *) obj%ReducedIntegration write (f%fh, '(A)') obj%ElemType write (f%fh, '(A)') obj%ErrorMsg call f%close() end if end subroutine end module ShapeFunctionClass