module PreprocessingClass use mpiclass use termclass use FEMDomainClass use PostProcessingClass implicit none type :: PreProcessing_ type(FEMDomain_) :: FEMDomain character*200 :: PictureName integer :: PixcelSize(2) integer :: ColorRGB(3) contains procedure :: Init => InitializePrePro procedure :: ImportPictureName => ImportPictureName procedure :: ShowName => ShowPictureName procedure :: ShowPixcelSize => ShowPixcelSize procedure :: GetPixcelSize => GetPixcelSize procedure :: GetAllPointCloud => GetAllPointCloud procedure :: SetColor => SetColor procedure :: ShowColor => ShowColor procedure :: GetPixcelByRGB => GetPixcelByRGB procedure :: GetSurfaceNode => GetPixcelSurfaceNode procedure :: AssembleSurfaceElement => AssembleSurfaceElement procedure :: ReduceSize => ReduceSize procedure :: ExportGeoFile => ExportGeoFile procedure :: ConvertGeo2Msh => ConvertGeo2Msh procedure :: ConvertGeo2Inp => ConvertGeo2Inp procedure :: ConvertGeo2Mesh => ConvertGeo2Mesh procedure :: ConvertMsh2Scf => ConvertMsh2Scf procedure :: ConvertMesh2Scf => ConvertMesh2Scf procedure :: Export => ExportPreProcessing procedure :: Reverse => ReversePreProcessing procedure :: SetDataType => SetDataTypeFEMDomain procedure :: SetSolver => SetSolverPreProcessing procedure :: SetUp => SetUpPreprocessing procedure :: SetScale => SetScalePreProcessing procedure :: SetBC => SetBoundaryConditionPrePro procedure :: SetSizeOfBC => SetSizeOfBCPrePrecessing procedure :: SetMatPara => SetMatParaPreProcessing procedure :: SetMatID => SetMatIDPreProcessing procedure :: ShowBC => ShowBCPrePrecessing procedure :: Convert3Dto2D => Convert3Dto2D procedure :: Convert2Dto3D => Convert2Dto3D procedure :: SetControlPara => SetControlParaPrePro procedure :: getSkelton => getSkeltonPreProcessing procedure :: Boolean => BooleanModifyerPreProcessing procedure :: setEntity => setEntityPreProcessing procedure :: showMesh => showMeshPreProcessing end type contains ! ######################################################### subroutine InitializePrePro(obj,Default) class(PreProcessing_),intent(inout)::obj logical,optional,intent(in)::Default call obj%FEMDomain%Init(Default) end subroutine ! ######################################################### ! ######################################################### subroutine ImportPictureName(obj,InPictureName) class(PreProcessing_)::obj character*200,intent(in)::InPictureName obj%PictureName=trim(InPictureName) end subroutine ! ######################################################### ! ######################################################### subroutine ShowPictureName(obj) class(PreProcessing_)::obj print *, trim(obj%PictureName) end subroutine ! ######################################################### ! ######################################################### subroutine ShowPixcelSize(obj) class(PreProcessing_)::obj character *20 :: pix_x character *20 :: pix_y write(pix_x,*) obj%PixcelSize(1) write(pix_y,*) obj%PixcelSize(2) print *, "Pixcel size is :: ",trim(adjustl(pix_x) )," x ",& trim(adjustl(pix_y) ) end subroutine ! ######################################################### ! ######################################################### subroutine GetPixcelSize(obj,MPIData,Name) class(PreProcessing_),intent(inout):: obj class(MPI_),optional,intent(inout) :: MPIData character(*),optional,intent(in) :: Name character *20 :: pid character *200 :: python_script character *200 :: python_buffer character *200 :: command integer :: fh call MPIData%GetInfo() write(pid,*) MPIData%MyRank python_script="GetPixcelSize_pid_"//trim(adjustl(pid))//".py" python_buffer="GetPixcelSize_pid_"//trim(adjustl(pid))//".txt" if(present(Name) )then python_script=Name//"GetPixcelSize_pid_"//trim(adjustl(pid))//".py" python_buffer=Name//"GetPixcelSize_pid_"//trim(adjustl(pid))//".txt" endif !print *, trim(python_script) ! using python script ! python imaging library is to be installed. fh=MPIData%MyRank+10 open(fh,file=trim(python_script),status="replace") command = "from PIL import Image" write(fh,'(A)') adjustl(trim(command)) command = "import sys" write(fh,'(A)') adjustl(trim(command)) command = "import os" write(fh,'(A)') adjustl(trim(command)) ! open file command = 'img_in = Image.open("'//trim(obj%PictureName)//'")' write(fh,'(A)') adjustl(trim(command)) command = 'python_buffer = open("'//trim(python_buffer)//'","w")' write(fh,'(A)') adjustl(trim(command)) ! get pixcel size command = "rgb_im = img_in.convert('RGB')" write(fh,'(A)') adjustl(trim(command)) command = "size = rgb_im.size" write(fh,'(A)') adjustl(trim(command)) command = "print( str(size[0]), ' ',str(size[1]) ) " write(fh,'(A)') adjustl(trim(command)) ! write size command = "python_buffer.write( str(size[0]))" write(fh,'(A)') adjustl(trim(command)) command = "python_buffer.write('\n')" write(fh,'(A)') adjustl(trim(command)) command = "python_buffer.write( str(size[1]))" write(fh,'(A)') adjustl(trim(command)) ! close command = "img_in.close()" write(fh,'(A)') adjustl(trim(command)) command = "python_buffer.close()" write(fh,'(A)') adjustl(trim(command)) close(fh) command = "python "//trim(python_script) print *, trim(command) call execute_command_line(trim(command)) ! get pixcel size open(fh,file=trim(python_buffer),status="old") read(fh,*) obj%PixcelSize(1) read(fh,*) obj%PixcelSize(2) close(fh) end subroutine ! ######################################################### ! ######################################################### subroutine GetAllPointCloud(obj,MPIData,Name) class(PreProcessing_),intent(inout):: obj class(MPI_),intent(inout) :: MPIData character(*),optional,intent(in) :: Name character *20 :: pid character *200 :: python_script character *200 :: python_buffer character *200 :: command integer :: fh call MPIData%GetInfo() write(pid,*) MPIData%MyRank python_script="GetAllPointCloud_pid_"//trim(adjustl(pid))//".py" python_buffer="GetAllPointCloud_pid_"//trim(adjustl(pid))//".txt" if(present(Name) )then python_script=Name//"GetAllPointCloud_pid_"//trim(adjustl(pid))//".py" python_buffer=Name//"GetAllPointCloud_pid_"//trim(adjustl(pid))//".txt" endif print *, trim(python_script) ! using python script ! python imaging library is to be installed. fh=MPIData%MyRank+10 open(fh,file=trim(python_script),status="replace") command = "from PIL import Image" write(fh,'(A)') adjustl(trim(command)) command = "import sys" write(fh,'(A)') adjustl(trim(command)) command = "import os" write(fh,'(A)') adjustl(trim(command)) ! open file command = 'img_in = Image.open("'//trim(obj%PictureName)//'")' write(fh,'(A)') adjustl(trim(command)) command = 'python_buffer = open("'//trim(python_buffer)//'","w")' write(fh,'(A)') adjustl(trim(command)) ! get pixcel size command = "rgb_im = img_in.convert('RGB')" write(fh,'(A)') adjustl(trim(command)) command = "size = rgb_im.size" write(fh,'(A)') adjustl(trim(command)) command = "print( str(size[0]), ' ',str(size[1]) ) " write(fh,'(A)') adjustl(trim(command)) ! get rgb pixcel coordinates command = "width,height =img_in.size" write(fh,'(A)') adjustl(trim(command)) command = "for i in range(width):" write(fh,'(A)') adjustl(trim(command)) command = "for j in range(height):" write(fh,'(A)') " "//adjustl(trim(command)) command = "R,G,B=rgb_im.getpixel((i,j))" write(fh,'(A)') " "//adjustl(trim(command)) command = "python_buffer.write(str(i)+'\t')" write(fh,'(A)') " "//adjustl(trim(command)) command = "python_buffer.write(str(j)+'\t')" write(fh,'(A)') " "//adjustl(trim(command)) command = "python_buffer.write(str(R)+'\t')" write(fh,'(A)') " "//adjustl(trim(command)) command = "python_buffer.write(str(G)+'\t')" write(fh,'(A)') " "//adjustl(trim(command)) command = "python_buffer.write(str(B)+'\n')" write(fh,'(A)') " "//adjustl(trim(command)) ! close command = "img_in.close()" write(fh,'(A)') adjustl(trim(command)) command = "python_buffer.close()" write(fh,'(A)') adjustl(trim(command)) close(fh) command = "python "//trim(python_script) print *, trim(command) call execute_command_line(trim(command)) end subroutine ! ######################################################### ! ######################################################### subroutine GetPixcelByRGB(obj,MPIData,err,onlycoord,Name) class(PreProcessing_),intent(inout):: obj class(MPI_),intent(inout) :: MPIData integer,optional,intent(in) :: err logical,optional,intent(in) :: onlycoord character(*),optional,intent(in) :: Name character *20 :: pid character *20 :: Red,Green,Blue character *20 :: er character *200 :: python_script character *200 :: python_buffer character *200 :: python_buffer_size character *200 :: command integer :: fh,error,sizeofpc,i if(present(err) )then error=err else error=0 endif call MPIData%GetInfo() write(pid,*) MPIData%MyRank write(Red,*) obj%ColorRGB(1) write(Green,*) obj%ColorRGB(2) write(Blue,*) obj%ColorRGB(3) write(er,*) error python_script="GetPixcelByRGB_pid_"//trim(adjustl(pid))//".py" python_buffer="GetPixcelByRGB_pid_"//trim(adjustl(pid))//".txt" python_buffer_size="GetPixcelByRGB_size_pid_"//trim(adjustl(pid))//".txt" if(present(Name) )then python_script =Name//"GetPixcelByRGB_pid_"//trim(adjustl(pid))//".py" python_buffer =Name//"GetPixcelByRGB_pid_"//trim(adjustl(pid))//".txt" python_buffer_size =Name//"GetPixcelByRGB_size_pid_"//trim(adjustl(pid))//".txt" endif print *, trim(python_script) ! using python script ! python imaging library is to be installed. fh=MPIData%MyRank+10 open(fh,file=trim(python_script),status="replace") command = "from PIL import Image" write(fh,'(A)') adjustl(trim(command)) command = "import sys" write(fh,'(A)') adjustl(trim(command)) command = "import os" write(fh,'(A)') adjustl(trim(command)) ! open file command = 'img_in = Image.open("'//trim(obj%PictureName)//'")' write(fh,'(A)') adjustl(trim(command)) command = 'python_buffer = open("'//trim(python_buffer)//'","w")' write(fh,'(A)') adjustl(trim(command)) command = 'python_buffer_size = open("'//trim(python_buffer_size)//'","w")' write(fh,'(A)') adjustl(trim(command)) ! get pixcel size command = "rgb_im = img_in.convert('RGB')" write(fh,'(A)') adjustl(trim(command)) command = "size = rgb_im.size" write(fh,'(A)') adjustl(trim(command)) command = "print( str(size[0]), ' ',str(size[1]) ) " write(fh,'(A)') adjustl(trim(command)) ! get rgb pixcel coordinates command = "itr = 0" write(fh,'(A)') adjustl(trim(command)) command = "width,height =img_in.size" write(fh,'(A)') adjustl(trim(command)) command = "for i in range(width):" write(fh,'(A)') adjustl(trim(command)) command = "for j in range(height):" write(fh,'(A)') " "//adjustl(trim(command)) command = "R,G,B=rgb_im.getpixel((i,j))" write(fh,'(A)') " "//adjustl(trim(command)) command = "er=abs(R-"//adjustl(trim(Red))//& ")+abs(G-"//adjustl(trim(Green))//& ")+abs(B-"//adjustl(trim(Blue))//")" write(fh,'(A)') " "//adjustl(trim(command)) command = "if er <= "//adjustl(trim(er))// " :" write(fh,'(A)') " "//adjustl(trim(command)) command = "python_buffer.write(str(i)+'\t')" write(fh,'(A)') " "//adjustl(trim(command)) command = "itr=itr+1" write(fh,'(A)') " "//adjustl(trim(command)) if( onlycoord .eqv. .true. )then command = "python_buffer.write(str(j)+'\n')" write(fh,'(A)') " "//adjustl(trim(command)) else command = "python_buffer.write(str(j)+'\t')" write(fh,'(A)') " "//adjustl(trim(command)) command = "python_buffer.write(str(R)+'\t')" write(fh,'(A)') " "//adjustl(trim(command)) command = "python_buffer.write(str(G)+'\t')" write(fh,'(A)') " "//adjustl(trim(command)) command = "python_buffer.write(str(B)+'\n')" write(fh,'(A)') " "//adjustl(trim(command)) endif ! close command = "python_buffer_size.write(str(itr)+'\n')" write(fh,'(A)') adjustl(trim(command)) command = "img_in.close()" write(fh,'(A)') adjustl(trim(command)) command = "python_buffer.close()" write(fh,'(A)') adjustl(trim(command)) command = "python_buffer_size.close()" write(fh,'(A)') adjustl(trim(command)) close(fh) command = "python "//trim(python_script) print *, trim(command) call execute_command_line(trim(command)) open(fh,file=python_buffer_size,status="old") read(fh,*) sizeofpc close(fh) allocate(obj%FEMDomain%Mesh%NodCoord(sizeofpc,3) ) obj%FEMDomain%Mesh%NodCoord(:,3)=0.0d0 open(fh,file=python_buffer,status="old") do i=1,sizeofpc read(fh,*)obj%FEMDomain%Mesh%NodCoord(i,1:2) enddo obj%FEMDomain%Mesh%NodCoord(i,2)=-1.0d0*obj%FEMDomain%Mesh%NodCoord(i,2) close(fh) end subroutine ! ######################################################### ! ######################################################### subroutine SetColor(obj,Red,Green,Blue) class(PreProcessing_),intent(inout):: obj integer,intent(in) :: Red,Green,Blue obj%ColorRGB(1)=Red obj%ColorRGB(2)=Green obj%ColorRGB(3)=Blue end subroutine ! ######################################################### ! ######################################################### subroutine ShowColor(obj) class(PreProcessing_),intent(inout):: obj print *, "Object Name is : ",trim(obj%PictureName) print *, "Red : ",obj%ColorRGB(1) print *, "Green : ",obj%ColorRGB(2) print *, "Blue : ",obj%ColorRGB(3) end subroutine ! ######################################################### ! ######################################################### subroutine GetPixcelSurfaceNode(obj,MPIData,r,NumOfMaxNod,Name) class(PreProcessing_),intent(inout):: obj class(MPI_),intent(inout) :: MPIData character(*),optional,intent(in) :: Name integer,optional,intent(in) :: r,NumOfMaxNod character*200 :: python_buffer character*20 :: pid integer,allocatable :: KilledPixcel(:) integer :: i,j,n,node_id,check_node_id,point_count,fh,MaxNod real(8) :: x_real,y_real,z_real real(8) :: x_real_tr,y_real_tr,z_real_tr,diff_real,max_r real(8),allocatable :: buffer(:,:) if(present(r) )then max_r=r else max_r=sqrt(2.10d0) endif if(present(NumOfMaxNod) )then MaxNod=NumOfMaxNod else MaxNod=7 endif n=size(obj%FEMDomain%Mesh%NodCoord,1) allocate(KilledPixcel(n) ) KilledPixcel(:)=0 ! remove isolated pixcel and surrounded nodes do i=1,n point_count=0 do j=1,n check_node_id=j x_real=obj%FEMDomain%Mesh%NodCoord(i,1) y_real=obj%FEMDomain%Mesh%NodCoord(i,2) z_real=obj%FEMDomain%Mesh%NodCoord(i,3) x_real_tr=obj%FEMDomain%Mesh%NodCoord(j,1) y_real_tr=obj%FEMDomain%Mesh%NodCoord(j,2) z_real_tr=obj%FEMDomain%Mesh%NodCoord(j,3) diff_real=(x_real-x_real_tr)*(x_real-x_real_tr)+& (y_real-y_real_tr)*(y_real-y_real_tr)+& (z_real-z_real_tr)*(z_real-z_real_tr) diff_real=dsqrt(diff_real) if(diff_real < max_r)then point_count=point_count+1 endif if(point_count > MaxNod )then KilledPixcel(i)=1 exit endif enddo if(point_count ==0 )then KilledPixcel(i)=1 endif enddo call MPIData%GetInfo() write(pid,*) MPIData%MyRank fh=MPIData%MyRank+10 python_buffer="GetSurface_pid_"//trim(adjustl(pid))//".txt" if(present(Name) )then python_buffer=Name//"GetSurface_pid_"//trim(adjustl(pid))//".txt" endif open(fh,file=python_buffer,status="replace") do i=1,n if(KilledPixcel(i)==0 )then write(fh,*) obj%FEMDomain%Mesh%NodCoord(i,1:2) else cycle endif enddo close(fh) do i=1,n point_count = point_count + KilledPixcel(i) enddo allocate(buffer(n-point_count,3) ) point_count=0 do i=1,n if(KilledPixcel(i)==0)then point_count = point_count + 1 buffer(point_count,:)=obj%FEMDomain%Mesh%NodCoord(i,:) else cycle endif enddo deallocate(obj%FEMDomain%Mesh%NodCoord) allocate(obj%FEMDomain%Mesh%NodCoord(size(buffer,1),size(buffer,2) ) ) obj%FEMDomain%Mesh%NodCoord(:,:)=buffer(:,:) end subroutine ! ######################################################### ! ######################################################### subroutine AssembleSurfaceElement(obj,MPIData,dim,threshold,DelRange,Name) class(PreProcessing_),intent(inout):: obj class(MPI_),intent(inout) :: MPIData integer,optional,intent(in) :: dim,threshold,DelRange character(*),optional,intent(in) :: Name character*200 :: python_buffer character*20 :: pid real(8),allocatable :: buffer(:,:),r_value(:) integer,allocatable :: checked(:) real(8) :: x(3),x_tr(3),r_tr,r_ref integer :: i,j,k,n,trial_num,id_tr,fh,r_threshold,drange if(present(threshold) )then r_threshold=threshold else r_threshold=10 endif if(present(DelRange) )then drange=DelRange else drange=10 endif if( present(dim) .and. dim/=2 )then call MPIData%End() stop "AssembleSurfaceElement :: >> only 2-D is available." endif n=size(obj%FEMDomain%Mesh%NodCoord,1 ) allocate(buffer(n,3)) buffer(:,:)=0.0d0 allocate(checked(n),r_value(n) ) checked(:)=0 buffer(1,:)=obj%FEMDomain%Mesh%NodCoord(1,:) checked(1)=1 do i=1,n-1 ! get buffer(i+1,:) trial_num=0 do j=1,n if(checked(j)==1 )then cycle endif trial_num=trial_num+1 x(:) =buffer(i,:) x_tr(:)=obj%FEMDomain%Mesh%NodCoord(j,:) r_tr=dsqrt(dot_product(x-x_tr,x-x_tr ) ) if(trial_num==1)then r_ref=r_tr id_tr=j else if(r_ref > r_tr)then id_tr=j r_ref=r_tr else cycle endif endif enddo buffer(i+1,:)=obj%FEMDomain%Mesh%NodCoord(id_tr,:) checked(id_tr)=1 enddo ! remove unnatural pixcel checked(:)=0 do i=1,size(buffer,1)-1 x(:) =buffer(i ,:) x_tr(:)=buffer(i+1,:) r_tr=dsqrt(dot_product(x-x_tr,x-x_tr ) ) if(r_tr > dble(r_threshold))then checked(i)=1 do k=1,drange if(i+k <= size(buffer,1))then checked(i+k)=1 endif if(i-k >= 1)then checked(i-k)=1 endif enddo print *, i,i+1 endif enddo if(maxval(checked)/=0 )then deallocate(obj%FEMDomain%Mesh%NodCoord) allocate(obj%FEMDomain%Mesh%NodCoord( n - sum(checked) ,3) ) k=0 do i=1,n if(checked(i)/=0 )then cycle else k=k+1 obj%FEMDomain%Mesh%NodCoord(k,1:2)=buffer(i,1:2) endif enddo else obj%FEMDomain%Mesh%NodCoord(:,:)=buffer(:,:) endif write(pid,*) MPIData%MyRank fh=MPIData%MyRank+10 python_buffer="GetSurface_pid_"//trim(adjustl(pid))//".txt" if(present(Name) )then python_buffer=Name//"GetSurface_pid_"//trim(adjustl(pid))//".txt" endif open(fh,file=python_buffer,status="replace") do i=1,size(obj%FEMDomain%Mesh%NodCoord,1) write(fh,*) obj%FEMDomain%Mesh%NodCoord(i,1:2) enddo close(fh) end subroutine ! ######################################################### ! ######################################################### subroutine ReduceSize(obj,MPIData,interval,Name) class(PreProcessing_),intent(inout):: obj class(MPI_),intent(inout) :: MPIData character(*),optional,intent(in) :: Name integer,intent(in) :: interval character*200 :: python_buffer character*20 :: pid real(8),allocatable:: buffer(:,:) integer,allocatable:: selected(:) integer :: i,j,k,n,m,fh n=size(obj%FEMDomain%Mesh%NodCoord,1) m=size(obj%FEMDomain%Mesh%NodCoord,2) allocate(selected(n) ) selected(:)=0 k=0 do i=1,n k=k+1 if(k==interval)then selected(i)=1 k=0 else cycle endif enddo allocate(buffer( sum(selected),3 )) k=0 do i=1,n if(selected(i)==1 )then k=k+1 buffer(k,:)= obj%FEMDomain%Mesh%NodCoord(i,:) endif enddo deallocate(obj%FEMDomain%Mesh%NodCoord) allocate(obj%FEMDomain%Mesh%NodCoord( size(buffer,1),size(buffer,2) ) ) obj%FEMDomain%Mesh%NodCoord(:,:)=buffer(:,:) write(pid,*) MPIData%MyRank fh=MPIData%MyRank+10 python_buffer="GetSurface_pid_"//trim(adjustl(pid))//".txt" if(present(Name) )then python_buffer=Name//"GetSurface_pid_"//trim(adjustl(pid))//".txt" endif open(fh,file=python_buffer,status="replace") do i=1,size(obj%FEMDomain%Mesh%NodCoord,1) write(fh,*) obj%FEMDomain%Mesh%NodCoord(i,1:2) enddo close(fh) end subroutine ! ######################################################### ! ######################################################### subroutine ExportGeoFile(obj,MPIData,Name) class(PreProcessing_),intent(inout):: obj character(*),optional,intent(in) :: Name class(MPI_),intent(inout) :: MPIData character*200 :: python_buffer character*20 :: pid integer :: i,j,k,n,fh,xsize real(8) :: x_p,y_p write(pid,*) MPIData%MyRank fh=MPIData%MyRank+10 python_buffer="GetSurface_pid_"//trim(adjustl(pid))//".geo" if(present(Name) )then python_buffer=Name//"GetSurface_pid_"//trim(adjustl(pid))//".geo" endif xsize = size(obj%FEMDomain%Mesh%NodCoord,1) open(fh,file=python_buffer) do i=1,xsize x_p=obj%FEMDomain%Mesh%NodCoord(i,1) y_p=obj%FEMDomain%Mesh%NodCoord(i,2) write(fh,*)"Point(",i,") = {",x_p,",",y_p,",","0, 1.0};" enddo write(fh,*)"Line(1) = {",xsize,",1};" do i=2,xsize write(fh,*)"Line(",i,") = {",i-1,",",i,"};" enddo write(fh,*)"Line Loop(",xsize+1,") = {" do i=1,xsize-1 write(fh,*)i,"," enddo write(fh,*)xsize,"};" write(fh,*)"Plane Surface(",xsize+2,") = {",xsize+1,"};" flush(fh) close(fh) !print *, "Done !!" !print *, "Generating mesh..." !call execute_command_line("gmsh.exe pm2.geo -2 -algo del2d -clmin 40") end subroutine ! ######################################################### ! ######################################################### subroutine ConvertGeo2Msh(obj,MPIData,Name) class(PreProcessing_),intent(inout):: obj class(MPI_),intent(inout) :: MPIData character(*),optional,intent(in) :: Name character*200 :: python_buffer character*200 :: command character*20 :: pid integer :: i,j,k,n,fh write(pid,*) MPIData%MyRank fh=MPIData%MyRank+10 python_buffer=" " command = " " python_buffer="GetSurface_pid_"//trim(adjustl(pid))//".geo" if(present(Name) )then python_buffer=Name//"GetSurface_pid_"//trim(adjustl(pid))//".geo" endif command="gmsh.exe "//trim(python_buffer)//" -2 -algo del2d -clmin 100" writE(*,'(A)') trim(command) call execute_command_line(command) !call execute_command_line("sh ./MakeMesh.sh") end subroutine ! ######################################################### ! ######################################################### subroutine ConvertGeo2Inp(obj,MPIData,Name) class(PreProcessing_),intent(inout):: obj class(MPI_),intent(inout) :: MPIData character(*),optional,intent(in) :: Name character*200 :: python_buffer character*200 :: command character*20 :: pid integer :: i,j,k,n,fh write(pid,*) MPIData%MyRank fh=MPIData%MyRank+10 python_buffer=" " command = " " python_buffer="GetSurface_pid_"//trim(adjustl(pid))//".geo" if(present(Name) )then python_buffer=Name//"GetSurface_pid_"//trim(adjustl(pid))//".geo" endif command="gmsh.exe "//trim(python_buffer)//" -2 -algo del2d -clmin 100 -format inp" writE(*,'(A)') trim(command) call execute_command_line(trim(command)) !call execute_command_line("sh ./MakeMesh.sh") end subroutine ! ######################################################### ! ######################################################### subroutine ConvertGeo2Mesh(obj,MPIData,SizePara,Name) class(PreProcessing_),intent(inout):: obj class(MPI_),intent(inout) :: MPIData integer,optional,intent(in) :: SizePara character(*),optional,intent(in) :: Name character*200 :: python_buffer character*200 :: command character*20 :: pid,a integer :: i,j,k,n,fh,sp if(present(SizePara) )then sp=SizePara else sp=10 endif write (a,*) sp write(pid,*) MPIData%MyRank fh=MPIData%MyRank+10 python_buffer=" " command = " " python_buffer="GetSurface_pid_"//trim(adjustl(pid))//".geo" if(present(Name) )then python_buffer=Name//"GetSurface_pid_"//trim(adjustl(pid))//".geo" endif command="gmsh.exe "//trim(python_buffer)//" -2 -algo del2d -clmin"//trim(a)//" -format mesh" writE(*,'(A)') trim(command) call execute_command_line(command) !call execute_command_line("./MakeMesh.sh") end subroutine ! ######################################################### ! ######################################################### subroutine ConvertMsh2Scf(obj,MPIData,ElementType,Name) class(PreProcessing_),intent(inout):: obj class(MPI_),intent(inout) :: MPIData character*200 :: python_buffer character*200 :: command,infile,outfile character*200,optional,intent(in) :: ElementType character(*),optional,intent(in) :: Name character*20 :: pid character*11 MeshFormat character*14 EndMeshFormat character*6 Nodes character*9 EndNodes,Elements character*12 EndElements integer,allocatable :: elem1(:),surface_nod(:) integer :: i,j,k,n,n1,n2,fh,a,nm,mm,nod_num,nn,elem_num,surf_num integer :: elem_num_all,n3,n4,n5,n6,n7,elemnod_num,startfrom real(8) :: re1,re2 ! ====================================================== ! deallocate all if(allocated(obj%FEMDomain%Mesh%NodCoord ))then deallocate(obj%FEMDomain%Mesh%NodCoord) endif if(allocated( obj%FEMDomain%Mesh%ElemNod ))then deallocate(obj%FEMDomain%Mesh%ElemNod) endif ! ====================================================== ! ====================================================== write(pid,*) MPIData%MyRank fh=MPIData%MyRank+10 infile="GetSurface_pid_"//trim(adjustl(pid))//".msh" outfile = "GetSurface_pid_"//trim(adjustl(pid))//".scf" if(present(Name) )then infile = Name//"GetSurface_pid_"//trim(adjustl(pid))//".msh" outfile = Name//"GetSurface_pid_"//trim(adjustl(pid))//".scf" endif open(fh,file=infile,status="old") print *, "Opening ",trim(infile) ! ====================================================== ! ====================================================== !read file to get nod and elem number read(fh,*)MeshFormat read(fh,*)re1,nm,mm read(fh,*)EndMeshFormat read(fh,*)Nodes if(nodes/="$Nodes")then if(nodes == "$Entit")then do read(fh,*)nodes if(nodes == "$Nodes")then print *, "Read entities" exit endif enddo else stop "ERROR: invalid location:$Nodes" endif endif ! ====================================================== ! Number of nodes read(fh,*)nod_num !print *,"nod_number",nod_num if(allocated(obj%FEMDomain%Mesh%NodCoord))then deallocate(obj%FEMDomain%Mesh%NodCoord) endif allocate(obj%FEMDomain%Mesh%NodCoord(nod_num,3) ) do i=1,nod_num read(fh,*) a if(i/=a)then stop "ERROR :: msh2scf" endif enddo read(fh,*)EndNodes if(EndNodes/="$EndNodes")then stop "ERROR: invalid location:$EndNodes" endif ! ====================================================== ! ====================================================== read(fh,*)Elements if(Elements/="$Elements")then stop "ERROR: invalid location: $Elements" endif read(fh,*)elem_num if(present(ElementType) )then if(trim(ElementType)=="LinearRectangularGp4")then elemnod_num= 4 elseif(trim(ElementType)=="LinearHexahedralGp8")then elemnod_num= 8 else print *, "PreProcessingClass.f90 >> Element : ",ElementType,"is not defined." return endif else ! default elemnod_num=4 endif startfrom=0 k=0 do i=1,elem_num read(fh,*)n1,n2,n3 if(n2==3 .and. n3==2)then if(startfrom==0)then startfrom=i endif k=k+1 endif enddo allocate(obj%FEMDomain%Mesh%ElemNod(k,elemnod_num) ) allocate( elem1(elemnod_num) ) read(fh,*)EndElements if(EndElements/="$EndElements")then stop "ERROR: invalid location: $EndElements" endif close(fh) ! ======================================================== ! ====================================================== write(pid,*) MPIData%MyRank fh=MPIData%MyRank+10 infile="GetSurface_pid_"//trim(adjustl(pid))//".msh" outfile = "GetSurface_pid_"//trim(adjustl(pid))//".scf" if(present(Name) )then infile = Name//"GetSurface_pid_"//trim(adjustl(pid))//".msh" outfile = Name//"GetSurface_pid_"//trim(adjustl(pid))//".scf" endif print *, "File Information is imported." open(fh,file=infile,status="old") ! ====================================================== ! ====================================================== !read file to get nod and elem number read(fh,*)MeshFormat read(fh,*)re1,nm,mm read(fh,*)EndMeshFormat read(fh,*)Nodes if(nodes/="$Nodes")then stop "ERROR: invalid location:$Nodes" endif ! ====================================================== ! Number of nodes read(fh,*)nod_num !print *,"nod_number",nod_num if(allocated(obj%FEMDomain%Mesh%NodCoord))then deallocate(obj%FEMDomain%Mesh%NodCoord) endif allocate(obj%FEMDomain%Mesh%NodCoord(nod_num,3) ) do i=1,nod_num read(fh,*) n2,obj%FEMDomain%Mesh%NodCoord(i,:) enddo read(fh,*)EndNodes if(EndNodes/="$EndNodes")then stop "ERROR: invalid location:$EndNodes" endif ! ====================================================== ! ====================================================== read(fh,*)Elements if(Elements/="$Elements")then stop "ERROR: invalid location: $Elements" endif read(fh,*)elem_num if(present(ElementType) )then if(trim(ElementType)=="LinearRectangularGp4")then elemnod_num= 4 elseif(trim(ElementType)=="LinearHexahedralGp8")then elemnod_num= 8 else print *, "PreProcessingClass.f90 >> Element : ",ElementType,"is not defined." return endif else ! default elemnod_num=4 endif k=0 do i=1,elem_num if(startfrom > i )then read(fh,*)n1,n2,n3 cycle endif if(i <= size(obj%FEMDomain%Mesh%ElemNod,1)+startfrom-1 )then k=k+1 read(fh,*)n1,n2,n3,n4,n5,obj%FEMDomain%Mesh%ElemNod(k,1:elemnod_num) cycle endif if(size(obj%FEMDomain%Mesh%ElemNod,1)+startfrom >= i)then read(fh,*) n1 endif enddo read(fh,*)EndElements if(EndElements/="$EndElements")then stop "ERROR: invalid location: $EndElements" endif close(fh) ! ======================================================== ! Setup FEMDomain end subroutine ! ######################################################### ! ######################################################### subroutine ConvertMesh2Scf(obj,MPIData,ElementType,Name) class(PreProcessing_),intent(inout):: obj class(MPI_),intent(inout) :: MPIData character(*),optional,intent(in) :: Name character*200 :: python_buffer character*200 :: command,infile,outfile character*200,optional,intent(in) :: ElementType character*20 :: pid character*11 MeshFormat character*14 EndMeshFormat character*6 Nodes character*9 EndNodes,Elements character*12 EndElements integer,allocatable :: elem1(:),surface_nod(:),triangle(:,:),devide_line(:,:) integer :: i,j,k,n,n1,n2,fh,a,nm,mm,nod_num,nn,elem_num,surf_num,l integer :: elem_num_all,n3,n4,n5,n6,n7,elemnod_num,startfrom,node1,node2,tr1,tr2 real(8) :: re1,re2 obj%FEMDomain%Mesh%ElemType = "LinearRectangularGp4" ! ====================================================== ! deallocate all if(allocated(obj%FEMDomain%Mesh%NodCoord ))then deallocate(obj%FEMDomain%Mesh%NodCoord) endif if(allocated( obj%FEMDomain%Mesh%ElemNod ))then deallocate(obj%FEMDomain%Mesh%ElemNod) endif ! ====================================================== ! ====================================================== write(pid,*) MPIData%MyRank fh=MPIData%MyRank+10 infile="GetSurface_pid_"//trim(adjustl(pid))//".mesh" outfile = "GetSurface_pid_"//trim(adjustl(pid))//".scf" if(present(Name) )then infile = Name//"GetSurface_pid_"//trim(adjustl(pid))//".mesh" outfile = Name//"GetSurface_pid_"//trim(adjustl(pid))//".scf" endif open(fh,file=infile,status="old") print *, "Opening ",trim(infile) ! ====================================================== ! ====================================================== !read file to get nod and elem number read(fh,*)MeshFormat read(fh,*)MeshFormat read(fh,*)mm read(fh,*)MeshFormat ! ====================================================== ! Number of nodes read(fh,*)nod_num !print *,"nod_number",nod_num if(allocated(obj%FEMDomain%Mesh%NodCoord))then deallocate(obj%FEMDomain%Mesh%NodCoord) endif allocate(obj%FEMDomain%Mesh%NodCoord(nod_num,3) ) obj%FEMDomain%Mesh%NodCoord(:,:)=0.0d0 do i=1,nod_num read(fh,*) obj%FEMDomain%Mesh%NodCoord(i,:) enddo read(fh,*)EndNodes ! ====================================================== read(fh,*)mm do i=1,mm read(fh,*) elem_num enddo read(fh,*)EndNodes ! ====================================================== ! ====================================================== read(fh,*)mm allocate(triangle(mm,4),devide_line(mm,3) ) devide_line(:,:)=-1 do i=1,mm read(fh,*) triangle(i,1:3) triangle(i,4)=-1 enddo read(fh,*)EndNodes ! ====================================================== ! ====================================================== read(fh,*)elem_num allocate(obj%FEMDomain%Mesh%ElemNod(elem_num,4)) obj%FEMDomain%Mesh%ElemNod(:,:)=-1 do i=1,elem_num read(fh,*) obj%FEMDomain%Mesh%ElemNod(i,1:4) enddo ! ====================================================== ! ====================================================== if(allocated(obj%FEMDomain%Mesh%ElemMat) )then deallocate(obj%FEMDomain%Mesh%ElemMat) endif allocate(obj%FEMDomain%Mesh%ElemMat(elem_num)) obj%FEMDomain%Mesh%ElemMat(:)=1 ! ====================================================== ! convert triangle do i=1,size(devide_line,1) do j=1,3 if(i==1)then node1=triangle(i,1) node2=triangle(i,2) elseif(i==2)then node1=triangle(i,2) node2=triangle(i,3) else node1=triangle(i,3) node2=triangle(i,1) endif do k=1,size(obj%FEMDomain%Mesh%ElemNod,1) do l=1,size(obj%FEMDomain%Mesh%ElemNod,2) if(l==1)then tr1=obj%FEMDomain%Mesh%ElemNod(k, 1) tr2=obj%FEMDomain%Mesh%ElemNod(k, 2) elseif(l==2)then tr1=obj%FEMDomain%Mesh%ElemNod(k, 2) tr2=obj%FEMDomain%Mesh%ElemNod(k, 3) elseif(l==3)then tr1=obj%FEMDomain%Mesh%ElemNod(k, 3) tr2=obj%FEMDomain%Mesh%ElemNod(k, 4) elseif(l==4)then tr1=obj%FEMDomain%Mesh%ElemNod(k, 4) tr2=obj%FEMDomain%Mesh%ElemNod(k, 1) else stop "ERROR :: ConvertMesh2Scf" endif if(node2==tr1 .and. node1 == tr2)then devide_line(i,j)=1 endif enddo enddo enddo enddo !call ShowArray(obj%FEMDomain%Mesh%NodCoord,triangle(:,1:3),20) !call ShowArray(obj%FEMDomain%Mesh%NodCoord,obj%FEMDomain%Mesh%ElemNod,30) !do i=1,size(devide_line,1) ! print *, devide_line(i,:) !enddo !stop "now debugging" end subroutine ! ######################################################### ! ######################################################### subroutine ExportPreProcessing(obj,MPIData,FileName,MeshDimension,Name) class(PreProcessing_),intent(inout):: obj class(MPI_),optional,intent(inout) :: MPIData character*200,optional,intent(in) :: FileName character(*),optional,intent(in) :: Name integer,optional,intent(in) :: MeshDimension character*200 :: python_buffer character*200 :: command character*20 :: pid integer :: i,j,k,n,fh fh=11 if(present(MPIData) )then write(pid,*) MPIData%MyRank fh=MPIData%MyRank+120 if(present(Name) )then python_buffer=Name//"GetSurface_pid_"//trim(adjustl(pid)) else python_buffer="GetSurface_pid_"//trim(adjustl(pid)) endif elseif(present(FileName) )then if(present(Name) )then python_buffer=Name//trim(FileName) else python_buffer=trim(FileName) endif else if(present(Name) )then python_buffer=Name else python_buffer="NoName" endif endif call obj%FEMDomain%Export(OptionalProjectName=python_buffer,FileHandle=fh) end subroutine ! ######################################################### !################################################## subroutine SetSolverPreProcessing(obj,inSolverType) class(PreProcessing_),intent(inout)::obj character*200,optional,intent(in) :: inSolverType character*200 ::sn if( present(inSolverType) )then sn = inSolverType else sn = "Default" endif Obj%FEMDomain%SolverType=sn end subroutine !################################################## !################################################## subroutine SetDataTypeFEMDomain(obj,inDType) class(PreProcessing_),intent(inout)::obj character*200,optional,intent(in) :: inDType character*200 :: sn sn="" if( .not.present(inDType) )then sn = "FEMDomain" else sn = inDType endif call obj%FEMDomain%SetDataType(trim(sn) ) end subroutine !################################################## !################################################## subroutine SetUpPreprocessing(obj,DataType,SolverType,NoFacetMode,MatPara) class(PreProcessing_),intent(inout)::obj character*200,optional,intent(in) :: DataType character*200,optional,intent(in) :: SolverType logical,optional,intent(in) :: NoFacetMode real(8),allocatable,optional,intent(inout)::MatPara(:,:) real(8),allocatable::MatParaDef(:,:) character*200 :: sn if(present(DataType) )then call obj%SetDataType(DataType) else call obj%SetDataType() endif call obj%FEMDomain%Mesh%Init(NoFacetMode=NoFacetMode) if(.not.present(MatPara) )then allocate(MatParaDef(1,1) ) MatParaDef(:,:)=100.0d0 call obj%FEMDomain%MaterialProp%Init(MaterialParameters=MatParaDef) else call obj%FEMDomain%MaterialProp%Init(MaterialParameters=MatPara) endif end subroutine !################################################## !################################################## subroutine SetScalePreProcessing(obj,scalex,scaley,scalez& ,picscalex,picscaley,picscalez) class(PreProcessing_),intent(inout)::obj real(8),optional,intent(in)::scalex,scaley,scalez real(8),optional,intent(in)::picscalex,picscaley,picscalez real(8) :: lx,ly,lz integer :: i if(present(scalex) )then if(scalex == 0.0d0)then stop "ERROR :: SetScalePreProcessing >> scalex==0.0d0" return else lx=maxval(obj%FEMDomain%Mesh%NodCoord(:,1) ) -& minval(obj%FEMDomain%Mesh%NodCoord(:,1) ) if(lx==0.0d0)then stop "ERROR :: SetScalePreProcessing >> lx==0.0d0" endif obj%FEMDomain%Mesh%NodCoord(:,1)=& obj%FEMDomain%Mesh%NodCoord(:,1)/lx*scalex endif endif if(present(scaley) )then if(scaley == 0.0d0)then stop "ERROR :: SetScalePreProcessing >> scaley==0.0d0" return else ly=maxval(obj%FEMDomain%Mesh%NodCoord(:,2) ) -& minval(obj%FEMDomain%Mesh%NodCoord(:,2) ) if(ly==0.0d0)then stop "ERROR :: SetScalePreProcessing >> ly==0.0d0" endif obj%FEMDomain%Mesh%NodCoord(:,2)=& obj%FEMDomain%Mesh%NodCoord(:,2)/ly*scaley endif endif if(present(scalez) )then if(scalez == 0.0d0)then stop "ERROR :: SetScalePreProcessing >> scalez==0.0d0" return else lz=maxval(obj%FEMDomain%Mesh%NodCoord(:,3) ) -& minval(obj%FEMDomain%Mesh%NodCoord(:,3) ) if(lz==0.0d0)then stop "ERROR :: SetScalePreProcessing >> lz==0.0d0" endif obj%FEMDomain%Mesh%NodCoord(:,3)=& obj%FEMDomain%Mesh%NodCoord(:,3)/lz*scalez endif endif do i=1,size(obj%FEMDomain%Mesh%NodCoord,1) write(123,*) obj%FEMDomain%Mesh%NodCoord(i,:) enddo end subroutine !################################################## !################################################## subroutine SetBoundaryConditionPrePro(obj,Dirichlet,Neumann,Initial,xmin,xmax,ymin,ymax,zmin,zmax,& tmin,tmax,val,val_id,NumOfValPerNod) class(PreProcessing_),intent(inout)::obj real(8),optional,intent(in)::xmin,xmax real(8),optional,intent(in)::ymin,ymax real(8),optional,intent(in)::zmin,zmax real(8),optional,intent(in)::tmin,tmax logical,optional,intent(in)::Dirichlet,Neumann,Initial integer,optional,intent(in)::NumOfValPerNod,val_id real(8),optional,intent(in)::val integer :: i,j,n if(present(Dirichlet) )then if(Dirichlet .eqv. .true.)then call AddDBoundCondition(obj%FEMDomain,xmin=xmin,xmax=xmax,ymin=ymin,ymax=ymax& ,zmin=zmin,zmax=zmax,tmin=tmin,tmax=tmax,val=val,val_id=val_id,NumOfValPerNod=NumOfValPerNod) return endif endif if(present(Neumann) )then if(Neumann .eqv. .true.)then call AddNBoundCondition(obj%FEMDomain,xmin=xmin,xmax=xmax,ymin=ymin,ymax=ymax& ,zmin=zmin,zmax=zmax,tmin=tmin,tmax=tmax,val=val,val_id=val_id,NumOfValPerNod=NumOfValPerNod) return endif endif if(present(Initial) )then if(Initial .eqv. .true.)then call AddTBoundCondition(obj%FEMDomain,xmin=xmin,xmax=xmax,ymin=ymin,ymax=ymax& ,zmin=zmin,zmax=zmax,tmin=tmin,tmax=tmax,val=val,val_id=val_id,NumOfValPerNod=NumOfValPerNod) return endif endif end subroutine !################################################## !################################################## subroutine SetSizeOfBCPrePrecessing(obj,Dirichlet,Neumann,Initial,NumOfValue) class(PreProcessing_),intent(inout)::obj logical,optional,intent(in)::Dirichlet,Neumann,Initial integer,optional,intent(in)::NumOfValue logical :: DB,NB,IC integer :: coord_dim if(.not. present(Dirichlet) )then DB = .false. else DB=Dirichlet endif if(.not. present(Neumann) )then NB = .false. else NB=Neumann endif if(.not. present(Initial) )then IC = .false. else IC=Initial endif coord_dim = size(obj%FEMDomain%Mesh%NodCoord,2) if(DB .eqv. .true.)then call obj%FEMDomain%InitDBC(NumOfValue) elseif(NB .eqv. .true.)then call obj%FEMDomain%InitNBC(NumOfValue) elseif(IC .eqv. .true.)then call obj%FEMDomain%InitTBC(NumOfValue) else return endif end subroutine !################################################## !################################################## subroutine ShowBCPrePrecessing(obj,Dirichlet,Neumann,Initial) class(PreProcessing_),intent(inout)::obj logical,optional,intent(in)::Dirichlet,Neumann,Initial integer :: n,m1,m2,i if(Dirichlet .eqv. .true.)then n=size(obj%FEMDomain%Boundary%DBoundNum) m1=size(obj%FEMDomain%Boundary%DBoundNodID,1) m2=size(obj%FEMDomain%Boundary%DBoundVal,1) print *, "Number of boundary conditions are :: ", n print *, "obj%FEMDomain%Boundary%DBoundNodID" do i=1,m1 print *, obj%FEMDomain%Boundary%DBoundNodID(i,:) enddo do i=1,m2 print *, obj%FEMDomain%Boundary%DBoundVal(i,:) enddo endif end subroutine !################################################## !################################################## subroutine Convert3Dto2D(obj) class(PreProcessing_),intent(inout)::obj real(8),allocatable::buffer(:,:) integer :: i,n,m n=size(obj%FEMDomain%Mesh%NodCoord,1) m=size(obj%FEMDomain%Mesh%NodCoord,2) allocate(buffer(n,2)) do i=1,n buffer(i,1:2)=obj%FEMDomain%Mesh%NodCoord(i,1:2) enddo deallocate(obj%FEMDomain%Mesh%NodCoord) allocate(obj%FEMDomain%Mesh%NodCoord(n,2) ) obj%FEMDomain%Mesh%NodCoord(:,:)=buffer(:,:) deallocate(buffer) end subroutine !################################################## !################################################## subroutine Convert2Dto3D(obj,Thickness,division) class(PreProcessing_),intent(inout)::obj real(8),allocatable::buffer(:,:) real(8),optional,intent(in)::Thickness integer,optional,intent(in)::division real(8) :: Tn integer :: i,j,n,m,NumOfLayer,numnod ! only for linear elements if(present(Thickness))then if(Thickness==0.0d0)then print *, "ERROR :: Convert2Dto3D >> Thickness = 0" return else Tn=Thickness endif else Tn=1.0d0 endif if(present(division))then if(division==0)then print *, "ERROR :: Convert2Dto3D >> division = 0" return endif NumOfLayer=division else NumOfLayer=1 endif numnod=size(obj%FEMDomain%Mesh%NodCoord,1) n=size(obj%FEMDomain%Mesh%NodCoord,1) m=size(obj%FEMDomain%Mesh%NodCoord,2) allocate(buffer(n*(NumOfLayer+1),3)) do j=1,NumOfLayer+1 do i=1,n buffer( n*(j-1) + i ,1:2) = obj%FEMDomain%Mesh%NodCoord(i,1:2) buffer( n*(j-1) + i ,3) = Tn / dble(NumOfLayer)*dble(j-1) enddo enddo deallocate(obj%FEMDomain%Mesh%NodCoord) allocate(obj%FEMDomain%Mesh%NodCoord( size(buffer,1) ,size(buffer,2) ) ) obj%FEMDomain%Mesh%NodCoord(:,:)=buffer(:,:) deallocate(buffer) ! ElemNod if(.not.allocated(obj%FEMDomain%Mesh%ElemNod) )then print *, "Caution :: Convert2Dto3D >> ElemNod is not allocated = 0" return endif n=size(obj%FEMDomain%Mesh%ElemNod,1) m=size(obj%FEMDomain%Mesh%ElemNod,2) allocate(buffer(n*NumOfLayer,m*2)) do j=1,NumOfLayer do i=1,n buffer( n*(j-1)+i, 1:m ) = obj%FEMDomain%Mesh%ElemNod(i,1:m)+numnod*(j-1) buffer( n*(j-1)+i, m+1:2*m ) = obj%FEMDomain%Mesh%ElemNod(i,1:m)+numnod*(j) enddo enddo deallocate(obj%FEMDomain%Mesh%ElemNod) allocate(obj%FEMDomain%Mesh%ElemNod( size(buffer,1) ,size(buffer,2) ) ) obj%FEMDomain%Mesh%ElemNod(:,:)=buffer(:,:) deallocate(buffer) ! ElemMat if(.not.allocated(obj%FEMDomain%Mesh%ElemMat) )then print *, "Caution :: Convert2Dto3D >> ElemMat is not allocated = 0" return endif allocate(buffer(n*NumOfLayer,1)) do j=1,NumOfLayer do i=1,n buffer( n*(j-1)+i, 1 ) = obj%FEMDomain%Mesh%ElemMat(i) enddo enddo deallocate(obj%FEMDomain%Mesh%ElemMat) allocate(obj%FEMDomain%Mesh%ElemMat( size(buffer,1) ) ) obj%FEMDomain%Mesh%ElemMat(:)=buffer(:,1) deallocate(buffer) end subroutine !################################################## !################################################## subroutine SetControlParaPrePro(obj,OptionalTol,OptionalItrTol,OptionalTimestep,OptionalSimMode) class(PreProcessing_),intent(inout)::obj real(8),optional,intent(in)::OptionalTol integer,optional,intent(in)::OptionalSimMode,OptionalItrTol,OptionalTimestep call SetControlPara(obj%FEMDomain%ControlPara,OptionalTol,OptionalItrTol,OptionalTimestep,OptionalSimMode) end subroutine !################################################## !################################################## subroutine SetMatParaPreProcessing(obj,NumOfMaterial,NumOfPara,MaterialID,ParameterID,Val,AllVal) class(PreProcessing_),intent(inout)::obj integer,optional,intent(in)::NumOfMaterial,NumOfPara,MaterialID,ParameterID real(8),optional,intent(in)::Val,AllVal(:) integer ::i,n,m,p(2),mm if(.not.allocated(obj%FEMDomain%MaterialProp%MatPara) )then allocate(obj%FEMDomain%MaterialProp%MatPara(1,1) ) obj%FEMDomain%MaterialProp%MatPara(1,1)=0.0d0 endif if(present(NumOfMaterial) )then if(NumOfMaterial > size(obj%FEMDomain%MaterialProp%MatPara,1) )then mm=size(obj%FEMDomain%MaterialProp%MatPara,1) do i=1,NumOfMaterial - mm call extendArray(obj%FEMDomain%MaterialProp%MatPara,extend1stColumn=.true.) if(present(val) )then obj%FEMDomain%MaterialProp%MatPara(mm+i,:)=val endif n=size(obj%FEMDomain%MaterialProp%MatPara,2) if(present(AllVal) )then m=size(AllVal) p(1)=n p(2)=m obj%FEMDomain%MaterialProp%MatPara(mm+i,1: minval(p) )=AllVal(1: minval(p) ) endif enddo endif endif if(present(NumOfPara) )then if(NumOfPara > size(obj%FEMDomain%MaterialProp%MatPara,2) )then mm=size(obj%FEMDomain%MaterialProp%MatPara,2) do i=1,NumOfPara - mm call extendArray(obj%FEMDomain%MaterialProp%MatPara,extend2ndColumn=.true.) if(present(val) )then obj%FEMDomain%MaterialProp%MatPara(:,mm+i)=val endif n=size(obj%FEMDomain%MaterialProp%MatPara,1) if(present(AllVal) )then m=size(AllVal) p(1)=n p(2)=m obj%FEMDomain%MaterialProp%MatPara(1: minval(p) , mm+i)=AllVal(1: minval(p) ) endif enddo endif endif end subroutine !################################################## !################################################## subroutine SetMatIDPreProcessing(obj,xmin,xmax,ymin,ymax,zmin,zmax,& tmin,tmax,MaterialID) class(PreProcessing_),intent(inout)::obj real(8),optional,intent(in)::xmin,xmax real(8),optional,intent(in)::ymin,ymax real(8),optional,intent(in)::zmin,zmax real(8),optional,intent(in)::tmin,tmax integer,optional,intent(in)::MaterialID integer :: i,j,n ! Now implementing call AddMaterialID(obj%FEMDomain,xmin=xmin,xmax=xmax,ymin=ymin,ymax=ymax& ,zmin=zmin,zmax=zmax,tmin=tmin,tmax=tmax,MaterialID=MaterialID) end subroutine !################################################## !################################################## subroutine getSkeltonPreProcessing(obj) class(PreProcessing_),intent(inout)::obj call obj%FEMDomain%MeltingSkelton() end subroutine !################################################## !################################################## subroutine setEntityPreProcessing(obj,Circle,Rectangle,Plane,Cylinder,Box,& Radius,XSize,YSize,ZSize,Xloc,Yloc,Zloc) class(PreProcessing_),intent(inout)::obj logical,optional,intent(in) :: Circle,Rectangle,Plane,Cylinder,Box real(8),optional,intent(in) :: Radius,XSize,YSize,ZSize,Xloc,Yloc,Zloc integer :: i if( present(Circle) )then if(Circle .eqv. .true.)then print *, "Under construction" return endif endif if( present(Rectangle) )then if(Rectangle .eqv. .true.)then if(.not. present(Xsize) )then print *, "Error :: setEntity >> Please import XSize" return endif if(.not. present(Ysize) )then print *, "Error :: setEntity >> Please import YSize" return endif if(allocated(obj%FEMDomain%Mesh%NodCoord) )then deallocate(obj%FEMDomain%Mesh%NodCoord) endif if(allocated(obj%FEMDomain%Mesh%ElemNod) )then deallocate(obj%FEMDomain%Mesh%ElemNod) endif if(allocated(obj%FEMDomain%Mesh%ElemMat) )then deallocate(obj%FEMDomain%Mesh%ElemMat) endif allocate(obj%FEMDomain%Mesh%NodCoord(4,2) ) obj%FEMDomain%Mesh%NodCoord(1,1)=0.0d0; obj%FEMDomain%Mesh%NodCoord(1,2)=0.0d0 obj%FEMDomain%Mesh%NodCoord(2,1)=Xsize; obj%FEMDomain%Mesh%NodCoord(2,2)=0.0d0 obj%FEMDomain%Mesh%NodCoord(3,1)=Xsize; obj%FEMDomain%Mesh%NodCoord(3,2)=ysize obj%FEMDomain%Mesh%NodCoord(4,1)=0.0d0; obj%FEMDomain%Mesh%NodCoord(4,2)=ysize allocate(obj%FEMDomain%Mesh%ElemNod(1,4) ) do i=1,4 obj%FEMDomain%Mesh%ElemNod(1,i)=i enddo allocate(obj%FEMDomain%Mesh%ElemMat(1) ) obj%FEMDomain%Mesh%ElemMat(1)=1 if( present(Xloc) )then obj%FEMDomain%Mesh%NodCoord(:,1)=obj%FEMDomain%Mesh%NodCoord(:,1)+Xloc endif if( present(Yloc) )then obj%FEMDomain%Mesh%NodCoord(:,2)=obj%FEMDomain%Mesh%NodCoord(:,2)+Yloc endif return endif endif if( present(Plane) )then if(Plane .eqv. .true.)then print *, "Under construction" return endif endif if( present(Cylinder) )then if(Cylinder .eqv. .true.)then print *, "Under construction" return endif endif if( present(Box) )then if(Box .eqv. .true.)then if(.not. present(Xsize) )then print *, "Error :: setEntity >> Please import XSize" return endif if(.not. present(Ysize) )then print *, "Error :: setEntity >> Please import YSize" return endif if(.not. present(Zsize) )then print *, "Error :: setEntity >> Please import YSize" return endif if(allocated(obj%FEMDomain%Mesh%NodCoord) )then deallocate(obj%FEMDomain%Mesh%NodCoord) endif if(allocated(obj%FEMDomain%Mesh%ElemNod) )then deallocate(obj%FEMDomain%Mesh%ElemNod) endif if(allocated(obj%FEMDomain%Mesh%ElemMat) )then deallocate(obj%FEMDomain%Mesh%ElemMat) endif allocate(obj%FEMDomain%Mesh%NodCoord(8,3) ) obj%FEMDomain%Mesh%NodCoord(1,1)=0.0d0; obj%FEMDomain%Mesh%NodCoord(1,2)=0.0d0; obj%FEMDomain%Mesh%NodCoord(1,3)=0.0d0; obj%FEMDomain%Mesh%NodCoord(2,1)=Xsize; obj%FEMDomain%Mesh%NodCoord(2,2)=0.0d0; obj%FEMDomain%Mesh%NodCoord(2,3)=0.0d0; obj%FEMDomain%Mesh%NodCoord(3,1)=Xsize; obj%FEMDomain%Mesh%NodCoord(3,2)=ysize; obj%FEMDomain%Mesh%NodCoord(3,3)=0.0d0; obj%FEMDomain%Mesh%NodCoord(4,1)=0.0d0; obj%FEMDomain%Mesh%NodCoord(4,2)=ysize; obj%FEMDomain%Mesh%NodCoord(4,3)=0.0d0; obj%FEMDomain%Mesh%NodCoord(5,1)=0.0d0; obj%FEMDomain%Mesh%NodCoord(5,2)=0.0d0; obj%FEMDomain%Mesh%NodCoord(5,3)=ZSize; obj%FEMDomain%Mesh%NodCoord(6,1)=Xsize; obj%FEMDomain%Mesh%NodCoord(6,2)=0.0d0; obj%FEMDomain%Mesh%NodCoord(6,3)=ZSize; obj%FEMDomain%Mesh%NodCoord(7,1)=Xsize; obj%FEMDomain%Mesh%NodCoord(7,2)=ysize; obj%FEMDomain%Mesh%NodCoord(7,3)=ZSize; obj%FEMDomain%Mesh%NodCoord(8,1)=0.0d0; obj%FEMDomain%Mesh%NodCoord(8,2)=ysize; obj%FEMDomain%Mesh%NodCoord(8,3)=ZSize; allocate(obj%FEMDomain%Mesh%ElemNod(1,8) ) obj%FEMDomain%Mesh%ElemNod(1,1)=1 obj%FEMDomain%Mesh%ElemNod(1,2)=2 obj%FEMDomain%Mesh%ElemNod(1,3)=3 obj%FEMDomain%Mesh%ElemNod(1,4)=4 obj%FEMDomain%Mesh%ElemNod(1,5)=5 obj%FEMDomain%Mesh%ElemNod(1,6)=6 obj%FEMDomain%Mesh%ElemNod(1,7)=7 obj%FEMDomain%Mesh%ElemNod(1,8)=8 allocate(obj%FEMDomain%Mesh%ElemMat(1) ) obj%FEMDomain%Mesh%ElemMat(1)=1 if( present(Xloc) )then obj%FEMDomain%Mesh%NodCoord(:,1)=obj%FEMDomain%Mesh%NodCoord(:,1)+Xloc endif if( present(Yloc) )then obj%FEMDomain%Mesh%NodCoord(:,2)=obj%FEMDomain%Mesh%NodCoord(:,2)+Yloc endif if( present(Zloc) )then obj%FEMDomain%Mesh%NodCoord(:,3)=obj%FEMDomain%Mesh%NodCoord(:,3)+Zloc endif return endif endif end subroutine !################################################## !################################################## subroutine BooleanModifyerPreProcessing(obj,ModObj,XDiv,Ydic,Zdiv) class(PreProcessing_),intent(inout)::obj class(PreProcessing_),intent(inout)::ModObj integer,optional,intent(in) :: XDiv,Ydic,Zdiv real(8) :: ground_level real(8),allocatable::RSInterface(:,:),xmin,xmax,ymin,ymax,NodCoord(:,:) integer :: ground_surface_id,n,m,itr,k,i,j,buf(2) integer :: NodeNum,DimNum,ElemNum,ElemNodNum,startnode,endnode,newnodenum integer,allocatable::RSIElemID(:),RSINodeID(:),RSIElemNod(:,:),AvailFE(:) integer,allocatable::OldNodID(:),OldtoNewNodID(:),countnode(:,:) ! ###### This is for 4-node elements or 8-node box element if(.not.allocated(obj%FEMDomain%Mesh%NodCoord) )then print *, "Boolean :: >> Error Please import obj%FEMDomain%Mesh%NodCoord" return endif if(.not.allocated(Modobj%FEMDomain%Mesh%NodCoord) )then print *, "Boolean :: >> Error Please import Modobj%FEMDomain%Mesh%NodCoord" return endif if(size(obj%FEMDomain%Mesh%NodCoord,2)==2)then ! if rectangle => ok if(size(obj%FEMDomain%Mesh%ElemNod,2) == 4)then ! ### only structural is supported ### print *, "Boolean :: ### only structural is supported ###" xmin=minval(obj%FEMDomain%Mesh%NodCoord(:,1) ) xmax=maxval(obj%FEMDomain%Mesh%NodCoord(:,1) ) ymin=minval(obj%FEMDomain%Mesh%NodCoord(:,2) ) ymax=maxval(obj%FEMDomain%Mesh%NodCoord(:,2) ) ground_level=maxval(obj%FEMDomain%Mesh%NodCoord(:,2) ) NodeNum=size(obj%FEMDomain%Mesh%NodCoord,1) DimNum=size(obj%FEMDomain%Mesh%NodCoord,2) ElemNum=size(obj%FEMDomain%Mesh%ElemNod,1) ElemNodNum=size(obj%FEMDomain%Mesh%ElemNod,2) allocate(RSIElemID(ElemNum),RSINodeID(NodeNum) ) allocate(OldNodID(NodeNum),OldtoNewNodID(NodeNum)) OldNodID(:)=0 OldtoNewNodID(:)=0 RSIElemID(:)=1 RSINodeID(:)=-1 if( .not.allocated(ModObj%FEMDomain%Mesh%FacetElemNod) )then call GetSurface(ModObj%FEMDomain%Mesh) endif call Obj%FEMDomain%Mesh%Copy(ModObj%FEMDomain%Mesh,Minimum=.true.) ! call showArray(-Obj%FEMDomain%Mesh%NodCoord,Obj%FEMDomain%Mesh%FacetElemNod,FileHandle=224) n=size(ModObj%FEMDomain%Mesh%FacetElemNod,1) m=size(ModObj%FEMDomain%Mesh%FacetElemNod,2) allocate(AvailFE(n) ) allocate(countnode( size(ModObj%FEMDomain%Mesh%NodCoord,1) ,2) ) countnode(:,:)=0 AvailFE(:)=0 itr=0 do i=1,n do j=1,m if(Modobj%FEMDomain%Mesh%NodCoord(ModObj%FEMDomain%Mesh%FacetElemNod(i,j) ,2)<=ground_level )then ! utilize countnode(ModObj%FEMDomain%Mesh%FacetElemNod(i,j),1)=countnode(ModObj%FEMDomain%Mesh%FacetElemNod(i,j),1)+1 countnode(ModObj%FEMDomain%Mesh%FacetElemNod(i,j),2)=i AvailFE(i)=1 itr=itr+1 exit endif enddo enddo k=0 do i=1,size(AvailFE) if(k==0 .and. AvailFE(i)==1 )then cycle endif if(k==0 .and. AvailFE(i)==0)then startnode=i k=1 endif if(k==1 .and. AvailFE(i)==0)then cycle endif if(k==1 .and. AvailFE(i)==1)then endnode=i exit endif enddo buf(1)=startnode buf(2)=endnode n=size(Obj%FEMDomain%Mesh%FacetElemNod,1) m=size(Obj%FEMDomain%Mesh%NodCoord,2) allocate(NodCoord(itr+4,m) ) NodCoord(:,:)=0.0d0 itr=0 do i=1,n if( i > minval(buf) )then exit endif if(AvailFE(i)==1 )then itr=itr+1 NodCoord(itr,:)=Obj%FEMDomain%Mesh%NodCoord(Obj%FEMDomain%Mesh%FacetElemNod(i,1) ,: ) endif enddo NodCoord(1+itr,1)=xmax ;NodCoord(1+itr,2)=ymax ; NodCoord(2+itr,1)=xmax ;NodCoord(2+itr,2)=ymin ; NodCoord(3+itr,1)=xmin ;NodCoord(3+itr,2)=ymin ; NodCoord(4+itr,1)=xmin ;NodCoord(4+itr,2)=ymax ; itr=itr+4 do i=1,n if( i < maxval(buf) )then cycle endif if(AvailFE(i)==1 )then itr=itr+1 NodCoord(itr,:)=Obj%FEMDomain%Mesh%NodCoord(Obj%FEMDomain%Mesh%FacetElemNod(i,1) ,: ) endif if(itr > n)then exit endif enddo deallocate(Obj%FEMDomain%Mesh%NodCoord) allocate(Obj%FEMDomain%Mesh%NodCoord(size(NodCoord,1),size(NodCoord,2) ) ) do i=1,size(NodCoord,1) Obj%FEMDomain%Mesh%NodCoord(i,:)=NodCoord(i,:) enddo !call showArray(NodCoord,FileHandle=226) !print *," " !print *, startnode,endnode return if(allocated(obj%FEMDomain%Mesh%ElemNod) )then deallocate(obj%FEMDomain%Mesh%ElemNod) endif allocate(obj%FEMDomain%Mesh%FacetElemNod(itr+4,m) ) ! FacetElements of under-gournd part of root domain is copyed to soil domain itr=0 do i=1,n if(AvailFE(i)==1 )then itr=itr+1 obj%FEMDomain%Mesh%FacetElemNod(itr,:)=Modobj%FEMDomain%Mesh%FacetElemNod(i,:) do j=1,m OldNodID( Modobj%FEMDomain%Mesh%FacetElemNod(i,j) )=1 enddo endif enddo ! list up old to new itr=0 do i=1,NodeNum if(OldNodID(i)==1 )then itr=itr+1 OldtoNewNodID(i)=itr endif enddo do i=1,size(obj%FEMDomain%Mesh%FacetElemNod,1) do j=1,size(obj%FEMDomain%Mesh%FacetElemNod,2) obj%FEMDomain%Mesh%FacetElemNod(i,j)=OldtoNewNodID( obj%FEMDomain%Mesh%FacetElemNod(i,j) ) if(obj%FEMDomain%Mesh%FacetElemNod(i,j) == 0 )then stop "BooleanModifyerPreProcessing :: ERROR :: OldtoNewNodID is wrong " endif enddo enddo NewNodeNum=maxval(OldtoNewNodID) if(allocated(obj%FEMDomain%Mesh%NodCoord) )then deallocate(obj%FEMDomain%Mesh%NodCoord) endif allocate( obj%FEMDomain%Mesh%NodCoord(NewNodeNum,DimNum) ) itr=0 do i=1,NodeNum if( OldNodID(i)==1 )then itr=itr+1 obj%FEMDomain%Mesh%NodCoord(itr,:)=Modobj%FEMDomain%Mesh%NodCoord(i,:) endif enddo allocate(countnode(NewNodeNum,2) ) countnode(:,:)=0 ! NodCoord, FacetElemNod are imported. ! Find edge of SurfaceNod ! Sort itr=0 do i=1,size(obj%FEMDomain%Mesh%FacetElemNod,1) do j=1,size(obj%FEMDomain%Mesh%FacetElemNod,2) countnode(obj%FEMDomain%Mesh%FacetElemNod(i,j),1 )=& countnode(obj%FEMDomain%Mesh%FacetElemNod(i,j),1 )+1 k=countnode(obj%FEMDomain%Mesh%FacetElemNod(i,j),2 ) if(k<j)then countnode(obj%FEMDomain%Mesh%FacetElemNod(i,j),2 )=j endif enddo enddo if(maxval(countnode(:,1) )==3 )then print *, "Boolean >> ERROR :: Same node appears 3 times " stop elseif(minval(countnode(:,1) )==2)then print *, "Boolean >> ERROR :: minval(countnode)==2 " stop else print *,"Boolean >> ERROR :: minval(countnode)<=0" stop endif itr=0 do i=1,size(countnode,1) if(countnode(i,1) == 1)then itr=itr+1 if(countnode(i,2)==1 )then startnode=i elseif(countnode(i,2)==2)then endnode=i else stop "Boolean >> ERROR :: countnode(i,2) /= (1 or 2) " endif endif enddo if(itr /= 2)then print *,"Boolean >> ERROR :: itr /= 2" stop endif ! あとはSurfaceNodを出力するだけ open(1233,file="test.txt") do i=1, size(obj%FEMDomain%Mesh%FacetElemNod,1) n=obj%FEMDomain%Mesh%FacetElemNod(i,1) write(1233,*) obj%FEMDomain%Mesh%NodCoord(i,:) do j=1,size(obj%FEMDomain%Mesh%FacetElemNod,1) enddo enddo !do i=1,NodeNum ! if(obj%FEMDomain%Mesh%NodCoord(i,2) <= ground_level )then ! RSINodeID(i)=1 ! endif !enddo ! !itr=0 !do i=1,ElemNum ! do j=1,ElemNodNum ! if(obj%FEMDomain%Mesh%ElemNod(i,j) == -1 )then ! RSIElemID(:)=-1 ! itr=itr+1 ! exit ! endif ! enddo !enddo !allocate(RSIElemNod(ElemNum-itr,ElemNodNum) ) !do i=1,ElemNum ! if(RSIElemNod(i) == 1 )then ! RSIElemNod(i,:)=obj%FEMDomain%Mesh%ElemNod(i,:) ! endif !enddo return else print *, "Boolean :: >> Error :: only rectangle is supported" return endif elseif(size(obj%FEMDomain%Mesh%NodCoord,2)==3)then ! if box => ok if(size(obj%FEMDomain%Mesh%ElemNod) == 8)then call Modobj%FEMDomain%Mesh%GetSurface() stop "Now implementing box" else print *, "Boolean :: >> Error :: only box is supported" return endif else print *, "Boolean :: >> Error size(obj%FEMDomain%Mesh%NodCoord,2) should be 2 or 3" return endif end subroutine !################################################## subroutine ReversePreProcessing(obj) class(PreProcessing_),intent(inout)::obj obj%FEMDomain%Mesh%NodCoord(:,:) = -obj%FEMDomain%Mesh%NodCoord(:,:) end subroutine !################################################## !################################################## subroutine showMeshPreProcessing(obj,Step,Name) class(PreProcessing_),intent(inout)::obj character(*),optional,intent(in):: Name integer,optional,intent(in):: Step integer :: stp if(present(Step) )then stp=step else stp=0 endif call GmshPlotMesh(obj%FEMDomain,OptionalStep=stp,Name=Name) end subroutine !################################################## end module