module FEMIfaceClass use, intrinsic :: iso_fortran_env use MathClass !use MPIClass use ArrayClass use ShapeFunctionClass use MeshClass use MaterialPropClass use BoundaryConditionClass use ControlParameterClass use FEMDomainClass implicit none type, extends(FEMDomain_) :: FEMDomainPointer_ type(FEMDomain_), pointer :: FEMDomainp end type type :: FEMIface_ type(ShapeFunction_) :: ShapeFunction1, ShapeFunction2 type(Mesh_) :: Mesh1, Mesh2 ! Mesh[12]%ElemNod is a LOCAL node pointer, not a DIRECT pointer for each domains type(MaterialProp_) :: MaterialProp type(ControlParameter_) :: ControlPara type(FEMDomainPointer_), allocatable :: FEMDomains(:) real(real64), allocatable :: NTN_NodCoord(:, :) real(real64), allocatable :: NTS_NodCoord(:, :) real(real64), allocatable :: STS_NodCoord(:, :) real(real64), allocatable :: NTN_Val(:, :) real(real64), allocatable :: NTS_Val(:, :) real(real64), allocatable :: STS_Val(:, :) integer(int32), allocatable :: NTN_ElemNod(:, :) integer(int32), allocatable :: NTS_ElemNod(:, :) integer(int32), allocatable :: STS_ElemNod(:, :) integer(int32), allocatable :: NTN_Active(:) integer(int32), allocatable :: NTS_Active(:) integer(int32), allocatable :: STS_Active(:) real(real64), allocatable :: NTN_Value(:, :) real(real64), allocatable :: NTS_Value(:, :) real(real64), allocatable :: STS_Value(:, :) integer(int32), allocatable :: NTS_SegmentID(:, :) integer(int32), allocatable :: GloNodPoint1(:, :) integer(int32), allocatable :: GloNodPoint2(:, :) integer(int32) :: DomainID1 integer(int32) :: DomainID2 integer(int32) :: DomainID3 integer(int32) :: TimeStep integer(int32) :: NumOfImportedDomain character*200 :: FilePathDomain1 character*200 :: FilePathDomain2 character*200 :: FilePath character*200 :: FileNameDomain1 character*200 :: FileNameDomain2 character*200 :: FileName character*9 :: Dtype character*20 :: SolverType contains procedure :: Init => InitializeFEMIface procedure :: setFEMDomain => setFEMDomainFEMIface procedure :: Delete => DeallocateFEMIface procedure :: Import => ImportFEMIface procedure :: GetFEMIface => GetFEMIfaceFromFEMDomains procedure :: Export => ExportFEMIface procedure :: GmshPlotMesh => GmshPlotMeshFEMIface procedure :: GmshPlotNTS => GmshPlotNTSFEMIface procedure :: GetNTNelement => GetNTNelement procedure :: GetNTSelement => GetNTSelement procedure :: GetGlobalNodePointer => GetGlobalNodePointerNTS procedure :: updateTimestep => updateTimestepIface end type contains ! ######################################################### subroutine InitializeFEMIface(obj, NumOfDomain) class(FEMIface_), intent(inout)::obj integer(int32), optional, intent(in)::NumOfDomain integer(int32) :: i call obj%Delete() if (allocated(obj%FEMDomains)) deallocate (obj%FEMDomains) if (allocated(obj%FEMDomains)) deallocate (obj%FEMDomains) if (present(NumOfDomain)) then allocate (obj%FEMDomains(NumOfDomain)) else allocate (obj%FEMDomains(2)) end if obj%NumOfImportedDomain = 0 obj%TimeStep = 0 end subroutine ! ######################################################### ! ######################################################### subroutine setFEMDomainFEMIface(obj, dobj, Name) class(FEMIface_), intent(inout)::obj class(FEMDomain_), target, intent(in)::dobj character(*), optional, intent(in) :: Name obj%NumOfImportedDomain = obj%NumOfImportedDomain + 1 if (size(obj%FEMDomains, 1) < obj%NumOfImportedDomain) then print *, "ERROR :: FEMIFace >> setFEMDomainFEMIface :: size(obj%FEMDomains,1) < obj%NumOfImportedDomain " else obj%FEMDomains(obj%NumOfImportedDomain)%FEMDomainp => dobj if (present(Name)) then obj%FEMDomains(obj%NumOfImportedDomain)%FEMDomainp%FileName = Name end if end if end subroutine ! ######################################################### ! ######################################################### subroutine GmshPlotMeshFEMIface(obj, Name, withNeumannBC, withDirichletBC) class(FEMIface_), intent(inout)::obj character(*), optional, intent(in) :: Name logical, optional, intent(in)::withNeumannBC, withDirichletBC integer(int32) :: i do i = 1, size(obj%FEMDomains, 1) if (present(Name)) then call obj%FEMDomains(i)%FEMDomainp%GmshPlotMesh(Name=Name//fstring(i), & withNeumannBC=withNeumannBC, withDirichletBC=withDirichletBC) else call obj%FEMDomains(i)%FEMDomainp%GmshPlotMesh(Name=fstring(i), & withNeumannBC=withNeumannBC, withDirichletBC=withDirichletBC) end if end do if (present(Name)) then call obj%GmshPlotNTS(Name=Name//"interface") else call obj%GmshPlotNTS(Name="interface") end if end subroutine ! ######################################################### ! ######################################################### subroutine DeallocateFEMIface(obj) class(FEMIface_), intent(inout)::obj if (allocated(obj%NTN_NodCoord)) then deallocate (obj%NTN_NodCoord) else end if if (allocated(obj%NTS_NodCoord)) then deallocate (obj%NTS_NodCoord) end if if (allocated(obj%STS_NodCoord)) then deallocate (obj%STS_NodCoord) end if if (allocated(obj%NTN_ElemNod)) then deallocate (obj%NTN_ElemNod) end if if (allocated(obj%NTS_ElemNod)) then deallocate (obj%NTS_ElemNod) end if if (allocated(obj%STS_ElemNod)) then deallocate (obj%STS_ElemNod) end if if (allocated(obj%FEMDomains)) then deallocate (obj%FEMDomains) end if end subroutine ! ######################################################### ! ######################################################### subroutine ImportFEMIface(obj, OptionalFileFormat, OptionalProjectName, FileHandle) class(FEMIface_), intent(inout)::obj character*4, optional, intent(in)::OptionalFileFormat character*70, optional, intent(in)::OptionalProjectName character*4::FileFormat character*70::ProjectName character*74 ::FileName character*9 :: DataType integer(int32), allocatable::IntMat(:, :) real(real64), allocatable::RealMat(:, :) integer(int32), optional, intent(in)::FileHandle integer(int32) :: fh, i, j, k, NumOfDomain, n, m, DimNum, GpNum, ierr character*70 Msg if (present(FileHandle)) then fh = FileHandle else fh = 104 end if if (present(OptionalFileFormat)) then FileFormat = OptionalFileFormat else FileFormat = ".scf" end if if (present(OptionalProjectName)) then ProjectName = OptionalProjectName else ProjectName = "untitled" end if obj%FileName = ProjectName FileName = ProjectName//FileFormat open (fh, file=FileName, status="old") if (FileFormat == ".scf") then read (fh, *) DataType if (DataType /= "interface" .and. DataType /= " interface") then return end if obj%Dtype = DataType read (fh, *) obj%FileNameDomain1 read (fh, *) obj%FileNameDomain2 read (fh, *) obj%SolverType end if close (fh) end subroutine ! ######################################################### !###################### Get FEM Interfaces ########################## subroutine GetFEMIfaceFromFEMDomains(obj, obj1, obj2, MasterID, SlaveID) class(FEMDomain_), optional, intent(inout)::obj1, obj2 class(FEMIface_), intent(inout)::obj integer(int32), optional, intent(in) ::MasterID, SlaveID integer(int32) :: i, j, n1, ierr, err if (.not. present(obj1) .and. .not. present(obj2)) then call GetFEMIfaceFromPointer(obj, MasterID, SlaveID) return end if print *, "object names #1 : ", obj1%FileName print *, "object names #2 : ", obj2%FileName call GetInterface(obj1%Mesh, obj2%Mesh, obj%Mesh1, obj%Mesh2, err) ! AABB algorithm if (err == 1) then print *, "no contact" return end if n1 = index(obj1%FileName, ".scf", back=.true.) obj%FilePathDomain1 = obj1%FilePath obj%FilePathDomain2 = obj2%FilePath obj%FilePath = obj1%FilePath obj%FileNameDomain1 = obj1%FileName obj%FileNameDomain2 = obj2%FileName obj%FileName = "Iface_"//obj1%FileName(1:n1 - 1)//"_"//obj2%FileName call obj%GetNTNelement() call obj%GetNTSelement() call ShowArray(obj%NTN_NodCoord, FileHandle=10) call ShowArray(obj%NTS_NodCoord, FileHandle=20) end subroutine !###################### Get FEM Interfaces ########################## !###################### Get FEM Interfaces ########################## subroutine GetFEMIfaceFromPointer(obj, MasterID, SlaveID) class(FEMIface_), intent(inout)::obj class(FEMDomain_), pointer::obj1, obj2 integer(int32), optional, intent(in) :: MasterID, SlaveID integer(int32) :: i, j, n1, ierr, err if (present(MasterID)) then i = MasterID else i = 1 end if if (present(SlaveID)) then j = SlaveID else j = 2 end if obj1 => obj%FEMDomains(i)%FEMDomainp obj2 => obj%FEMDomains(j)%FEMDomainp print *, "object names #1 : ", obj1%FileName print *, "object names #2 : ", obj2%FileName call GetInterface(obj1%Mesh, obj2%Mesh, obj%Mesh1, obj%Mesh2, err) ! AABB algorithm call obj%GetGlobalNodePointer() !if(err==1)then ! print *, "no contact" ! return !endif !n1 = index(obj1%FileName,".scf", back=.true. ) !obj%FilePathDomain1=obj1%FilePath !obj%FilePathDomain2=obj2%FilePath !obj%FilePath =obj1%FilePath !obj%FileNameDomain1=obj1%FileName !obj%FileNameDomain2=obj2%FileName !obj%FileName ="Iface_"//obj1%FileName(1:n1-1)//"_"//obj2%FileName call obj%GetNTNelement() call obj%GetNTSelement() call ShowArray(obj%NTN_NodCoord, FileHandle=10) call ShowArray(obj%NTS_NodCoord, FileHandle=20) end subroutine !###################### Get FEM Interfaces ########################## !###################### Get Node-To-Node Elements ########################## subroutine GetNTNelement(obj) class(FEMIface_), intent(inout)::obj real(real64), allocatable::x(:), xn(:) real(real64) :: dist integer(int32) :: node_num1, dim_num1, node_num2, dim_num2, i, j, n, node_num integer(int32) :: master, dim_num, id node_num1 = size(obj%Mesh1%NodCoord, 1) dim_num1 = size(obj%Mesh1%NodCoord, 2) node_num2 = size(obj%Mesh2%NodCoord, 1) dim_num2 = size(obj%Mesh2%NodCoord, 2) if (dim_num1 /= dim_num2) then stop "ERROR :: GetNTNelement dimension of domain1 and domain2 is not consistent" end if dim_num = dim_num1 allocate (x(dim_num1)) allocate (xn(dim_num1)) if (node_num1 >= node_num2) then node_num = node_num2 master = 2 else node_num = node_num1 master = 1 end if if (allocated(obj%NTN_NodCoord)) deallocate (obj%NTN_NodCoord) if (allocated(obj%NTN_ElemNod)) deallocate (obj%NTN_ElemNod) if (allocated(obj%NTN_Val)) deallocate (obj%NTN_Val) allocate (obj%NTN_NodCoord(node_num, dim_num*2)) !In terms of local IDs allocate (obj%NTN_ElemNod(node_num, 2)) !In terms of local IDs allocate (obj%NTN_Val(node_num, 2)) !In terms of local IDs do i = 1, node_num if (master == 1) then ! domain 1 is master >> search the pair from domain 2 x(:) = obj%Mesh1%NodCoord(i, :) id = SearchNearestCoord(obj%Mesh2%NodCoord, x) obj%NTN_ElemNod(i, 1) = i obj%NTN_ElemNod(i, 2) = id xn(:) = obj%Mesh2%NodCoord(id, :) obj%NTN_Val(i, 1) = dsqrt(dot_product(x - xn, x - xn)) obj%NTN_NodCoord(i, 1:dim_num) = x(:) obj%NTN_NodCoord(i, dim_num + 1:2*dim_num) = xn(:) else ! domain 1 is master >> search the pair from domain 2 x(:) = obj%Mesh2%NodCoord(i, :) id = SearchNearestCoord(obj%Mesh1%NodCoord, x) obj%NTN_ElemNod(i, 2) = i obj%NTN_ElemNod(i, 1) = id xn(:) = obj%Mesh1%NodCoord(id, :) obj%NTN_Val(i, 1) = dsqrt(dot_product(x - xn, x - xn)) obj%NTN_NodCoord(i, 1:dim_num) = xn(:) obj%NTN_NodCoord(i, dim_num + 1:2*dim_num) = x(:) end if end do end subroutine !###################### Get Node-To-Node Elements ########################## !###################### Get Node-To-Segment Elements ########################## subroutine GetNTSelement(obj) class(FEMIface_), intent(inout)::obj real(real64), allocatable::x(:), xn(:), ElemMidPointCoord(:, :) real(real64) :: dist integer(int32) :: node_num1, dim_num1, node_num2, dim_num2, i, j, n, node_num, elem_num2, elem_num integer(int32) :: slave, dim_num, id, elemnod_num2, elemnod_num node_num1 = size(obj%Mesh1%NodCoord, 1) dim_num1 = size(obj%Mesh1%NodCoord, 2) dim_num2 = size(obj%Mesh2%NodCoord, 2) elem_num2 = size(obj%Mesh2%ElemNod, 1) elemnod_num2 = size(obj%Mesh2%ElemNod, 2) if (dim_num1 /= dim_num2) then stop "ERROR :: GetNTSelement dimension of domain1 and domain2 is not consistent" end if dim_num = dim_num1 node_num = node_num1 ! slave node elem_num = elem_num2 ! master segment elemnod_num = elemnod_num2 ! master segment allocate (x(dim_num1)) allocate (xn(dim_num1)) allocate (ElemMidPointCoord(size(obj%Mesh2%ElemNod, 1), dim_num)) do i = 1, elem_num xn(:) = 0.0d0 do j = 1, elemnod_num xn(:) = xn(:) + 1.0d0/dble(elemnod_num)*obj%Mesh2%NodCoord(obj%Mesh2%ElemNod(i, j), :) end do ElemMidPointCoord(i, :) = xn(:) end do if (allocated(obj%NTS_NodCoord)) deallocate (obj%NTS_NodCoord) if (allocated(obj%NTS_ElemNod)) deallocate (obj%NTS_ElemNod) if (allocated(obj%NTS_Val)) deallocate (obj%NTS_Val) if (allocated(obj%NTS_SegmentID)) deallocate (obj%NTS_SegmentID) allocate (obj%NTS_NodCoord(node_num, dim_num*(1 + elemnod_num))) !In terms of local IDs allocate (obj%NTS_ElemNod(node_num, 1 + elemnod_num)) !In terms of local IDs allocate (obj%NTS_Val(node_num, 2)) !In terms of local IDs allocate (obj%NTS_SegmentID(node_num, 1)) obj%NTS_Val(:, :) = 0.0d0 do i = 1, node_num x(:) = obj%Mesh1%NodCoord(i, :) id = SearchNearestCoord(ElemMidPointCoord, x) obj%NTS_ElemNod(i, 1) = i obj%NTS_SegmentID(i, 1) = id obj%NTS_NodCoord(i, 1:dim_num) = x(:) do j = 1, elemnod_num obj%NTS_ElemNod(i, 1 + j) = obj%Mesh2%ElemNod(id, j) obj%NTS_NodCoord(i, dim_num*(j) + 1:dim_num*(j + 1)) = & obj%Mesh2%NodCoord(obj%Mesh2%ElemNod(id, j), :) end do obj%NTS_NodCoord(i, 1:dim_num) = x(:) end do end subroutine !###################### Get Node-To-Segment Elements ########################## !###################### Get Segment-To-Segment Elements ########################## subroutine GetSTSelement(obj) class(FEMIface_), intent(inout)::obj print *, "now, developping" end subroutine !###################### Get Segment-To-Segment Elements ########################## ! ######################################################### subroutine ExportFEMIface(obj, OptionalFileFormat, OptionalProjectName, FileHandle) class(FEMIface_), intent(inout)::obj character*4, optional, intent(in)::OptionalFileFormat character*70, optional, intent(in)::OptionalProjectName character*4::FileFormat character*70::ProjectName character*74 ::FileName character*9 :: DataType integer(int32), allocatable::IntMat(:, :) real(real64), allocatable::RealMat(:, :) integer(int32), optional, intent(in)::FileHandle integer(int32) :: fh, i, j, k, NumOfDomain, n, m, DimNum, GpNum character*70 Msg if (present(FileHandle)) then fh = FileHandle else fh = 104 end if if (present(OptionalFileFormat)) then FileFormat = OptionalFileFormat else FileFormat = ".scf" end if if (present(OptionalProjectName)) then ProjectName = OptionalProjectName else ProjectName = "untitled" end if FileName = ProjectName//FileFormat open (fh, file=FileName, status="replace") if (FileFormat == ".scf") then if (obj%Dtype /= "interface") then return end if write (fh, '(A)') obj%Dtype write (fh, '(A)') obj%FileNameDomain1 write (fh, '(A)') obj%FileNameDomain2 write (fh, '(A)') obj%SolverType end if close (fh) end subroutine ! ######################################################### ! ######################################################### subroutine GmshPlotNTSFEMIface(obj, Name) class(FEMIFace_), intent(in)::obj type(FEMDomain_) :: Fobj integer(int32) :: i, j, n, m, dim_num, k, elemnodnum, step character(*), optional, intent(in)::Name ! Visualize NTS elements step = obj%TimeStep dim_num = size(obj%FEMDomains(1)%FEMDomainp%Mesh%NodCoord, 2) print *, dim_num n = size(obj%NTS_ElemNod, 1) m = size(obj%NTS_ElemNod, 2) if (dim_num == 2) then ! 4 elemnodnum = 4 elseif (dim_num == 3) then ! 8 elemnodnum = 8 else print *, "ERROR dim_num of NTS_NodCoord is ", dim_num return end if allocate (Fobj%Mesh%NodCoord(n*m, dim_num)) allocate (Fobj%Mesh%ElemNod(n, elemnodnum)) allocate (Fobj%Mesh%ElemMat(n)) Fobj%FileName = obj%FileName Fobj%Mesh%ElemMat(:) = 1 Fobj%Mesh%NodCoord(:, :) = 0.0d0 Fobj%Mesh%ElemNod(:, :) = 0 ! create Fobj%Mesh%NodCoord k = 0 do i = 1, n ! number of NTS elements do j = 1, m ! over nodes in a NTS elements k = k + 1 ! ID if (elemnodnum == 4) then Fobj%Mesh%NodCoord(k, 1:2) = obj%NTS_NodCoord(i, dim_num*(j - 1) + 1:dim_num*(j)) ! need revision elseif (elemnodnum == 8) then if (j == 1) then Fobj%Mesh%ElemNod(i, j + elemnodnum/2) = k else Fobj%Mesh%ElemNod(i, j - 1) = k if (j - 1 /= 1) then Fobj%Mesh%ElemNod(i, j - 1 + elemnodnum/2) = k end if end if Fobj%Mesh%NodCoord(k, 1:3) = obj%NTS_NodCoord(i, dim_num*(j - 1) + 1:dim_num*(j)) else print *, "ERROR2 dim_num of NTS_NodCoord is ", dim_num return end if end do end do call Fobj%GmshPlotMesh(Name=Name, OptionalStep=step) end subroutine ! ######################################################### ! ######################################################### subroutine GetGlobalNodePointerNTS(obj) class(FEMIface_)::obj !type(MPI_)::mpidata integer(int32) :: i, j, k, n, m, NumElemIface1, NumElemIface2 real(real64), allocatable :: x(:), x_tr(:) ! get GlobalNodePointer of NTS element NumElemIface1 = size(obj%Mesh1%ElemNod, 1) NumElemIface2 = size(obj%Mesh2%ElemNod, 1) n = size(obj%Mesh1%ElemNod, 2) m = size(obj%Mesh2%ElemNod, 2) if (allocated(obj%GloNodPoint1)) then deallocate (obj%GloNodPoint1) end if if (allocated(obj%GloNodPoint2)) then deallocate (obj%GloNodPoint2) end if allocate (obj%GloNodPoint1(NumElemIface1, n), obj%GloNodPoint2(NumElemIface2, m)) ! for Domain1 do i = 1, NumElemIface1 do j = 1, n !print *, obj%Mesh1%GlobalNodID(obj%Mesh1%ElemNod(i,j) ),"/",size(obj%FEMDomains(1)%FEMDomainp%Mesh%NodCoord,1) obj%GloNodPoint1(i, j) = obj%Mesh1%GlobalNodID(obj%Mesh1%ElemNod(i, j)) end do end do ! ! !! for Domain2 do i = 1, NumElemIface2 do j = 1, m !print *, obj%Mesh2%GlobalNodID(obj%Mesh2%ElemNod(i,j) ),"/",size(obj%FEMDomains(2)%FEMDomainp%Mesh%NodCoord,1) obj%GloNodPoint2(i, j) = obj%Mesh2%GlobalNodID(obj%Mesh2%ElemNod(i, j)) end do end do end subroutine ! ######################################################### ! ######################################################### subroutine updateTimestepIface(obj, timestep) class(FEMIFace_), intent(inout)::obj integer(int32), optional, intent(in)::timestep integer(int32) :: dt dt = input(default=1, option=timestep) obj%TimeStep = obj%TimeStep + dt end subroutine ! ######################################################### end module