PreProcessingClass.f90 Source File


Contents


Source Code

module PreprocessingClass
    use, intrinsic :: iso_fortran_env
    use std
    use FEMDomainClass
    use PostProcessingClass

    
    implicit none    
    
    type :: PreProcessing_
        type(FEMDomain_) :: FEMDomain
        type(FEMDomain_),pointer :: pFEMDomain
        character*200   :: PictureName
        character*200   :: RGBDataName,PixcelSizeDataName
        integer(int32)         :: PixcelSize(2),num_of_pixcel
        integer(int32)         :: ColorRGB(3)
    contains
        !procedure :: set                => setPreprocessing
        procedure :: getScfFromImage    => getScfFromImagePreProcessing
        procedure :: Init               => InitializePrePro
        procedure :: finalize           => finalizePrePro
        procedure :: ImportPictureName  => ImportPictureName
        procedure :: importPixcelAsNode => importPixcelAsNodePreProcessing
        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 :: modifySuefaceNode => modifySuefaceNodePrepro
        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 :: ConvertGeo2VTK     => ConvertGeo2VTK
        procedure :: Export             => ExportPreProcessing
        procedure :: ExportAsLodgingSim => ExportAsLodgingSimProcessing
        procedure :: import             => importPreProcessing
        procedure :: Reverse            => ReversePreProcessing
        procedure :: SetDataType        => SetDataTypeFEMDomain
        procedure :: SetSolver          => SetSolverPreProcessing 
        procedure :: SetUp              => SetUpPreprocessing
        procedure :: SetScale           => SetScalePreProcessing
        procedure :: SetBC              => SetBoundaryConditionPrePro
        procedure :: removeBC           => removeBoundaryConditionPrePro
        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 
        procedure :: meshing            => meshingPreProcessing
        
    end type
    
contains

! #########################################################
subroutine getScfFromImagePreProcessing(obj,project,ElemType,MPIData,R,G,B,scalex,scaley,&
    Soilfile,sR,SG,sB,SolverName)
    class(PreProcessing_),intent(inout) :: obj
    class(MPI_),intent(inout)           :: MPIData

    type(Dictionary_)       :: InfileList,DBoundlist,NBoundlist,Materialist
    type(PreProcessing_)    :: leaf,soil
    character(*),intent(in) :: project,elemtype,SolverName
    character(*),optional,intent(in) :: Soilfile
    integer(int32),intent(in) :: R,G,B
    integer(int32),optional,intent(in) :: sR,SG,sB
    real(real64),intent(in) :: scalex,scaley
    real(real64) :: Dbound_val,Nbound_val,xratio,yratio
    character(200)         :: name,name1,name2,name3,name4,str_id,&
        sname,dirichlet,neumann,materials,parameters
    integer(int32) :: NumOfImages,i,id,num_d,num_n,DBoundRGB(3),Dbound_xyz,NBoundRGB(3),Nbound_xyz
    integer(int32) :: NumOfMaterial,NumOfparameter,matid,MaterialRGB(3)
    real(real64),allocatable :: matpara(:)

    

    if(trim(ElemType) /= "LinearRectangularGp4")then
        print *, "ERROR :: now only LinearRectangularGp4 is available."
        return
    endif


    ! get paths for Image lists
    open(50,file=trim(project)//"filenamelist.txt")
    read(50,*) NumOfImages
    call InfileList%Init(NumOfImages)
    do i=1,NumOfImages
        read(50,'(A)' ) name
        call InfileList%Input(i, trim(name) )
    enddo
    close(50)
    call MPIData%createStack(total=NumOfImages)

    ! get boundary information list
    open(60, file=trim(project)//"boundcondlist.txt")
    read(60,*) dirichlet
    read(60,*) num_d
    call DBoundlist%Init(num_d)
    do i=1,num_d
        read(60,'(A)' ) name
        read(60,*)  DBoundRGB(1:3)
        read(60,*)  Dbound_xyz, Dbound_val
        call DBoundlist%Input(i, content=trim(name) )
        call DBoundlist%Input(i, intlist=DBoundRGB )
        call DBoundlist%Input(i, IntValue=Dbound_xyz )
        call DBoundlist%Input(i, RealValue=Dbound_val )
    enddo
    read(60,*) neumann
    read(60,*) num_n
    do i=1,num_n
        read(60,'(A)' ) name
        read(60,*)  NBoundRGB(1:3)
        read(60,*)  Nbound_xyz, Nbound_val
        call NBoundlist%Input(i, content=trim(name) )
        call NBoundlist%Input(i, Intlist=NBoundRGB )
        call NBoundlist%Input(i, IntValue=Nbound_xyz )
        call NBoundlist%Input(i, RealValue=Nbound_val )
    enddo
    close(60)


    ! get paths for material information lists
    open(70,file=trim(project)//"materialist.txt")
    read(70, '(A)' ) materials
    read(70,*) NumOfMaterial
    read(70, '(A)' ) parameters
    read(70,*) NumOfparameter
    allocate(matpara(NumOfparameter))
    
    call Materialist%Init(NumOfMaterial)
    do i=1,NumOfMaterial
        read(70,'(A)' ) name
        read(70,*) MaterialRGB(1:3)
        read(70,*) matpara(1:NumOfparameter)
        call materialist%Input(i, content=trim(name) )
        call materialist%Input(i, Intlist=materialRGB )
        call materialist%Input(i, Realist=matpara )
    !    call Materialist%Input(i, trim(name) )
    enddo
    close(70)

    
    
    do i=1,size(MPIData%LocalStack)
        id=MPIData%LocalStack(i)
        str_id= trim(adjustl(fstring(id)))
        name=trim(InfileList%get( MPIData%LocalStack(i)) )
        print *, "MyRank",MPIData%MyRank,"|",trim(name)
        
        ! Get Pixcel
        call leaf%ImportPictureName(name)
        call leaf%GetPixcelSize(MPIData)
        call leaf%SetColor(R,G,B)
        call leaf%GetPixcelByRGB(MPIData,err=5,onlycoord=.true.)
        ! Get Outline
        call leaf%GetSurfaceNode(MPIData)
        call leaf%AssembleSurfaceElement(MPIData,dim=2,threshold=5,DelRange=5)
        
        ! Convert SurfaceNod to .geo
        call leaf%ExportGeoFile(MPIData,Name=trim(project)//"mesh"//trim(str_id)//".geo" )
        
        ! Run Gmsh to convert .geo to .msh
        call leaf%ConvertGeo2Msh(MPIData ,Name=trim(project)//"mesh"//trim(str_id)//".geo" )
        call leaf%ConvertGeo2Inp(MPIData ,Name=trim(project)//"mesh"//trim(str_id)//".geo" )
        call leaf%ConvertGeo2Mesh(MPIData,Name=trim(project)//"mesh"//trim(str_id)//".geo" )
        
        ! Convert .msh to .scf
        
        call leaf%ConvertMesh2Scf(MPIData,ElementType=ElemType,&
            Name=trim(project)//"mesh"//trim(str_id)//".mesh" )
        call leaf%FEMDomain%checkconnectivity(fix=.true.)

        !call leaf%Convert3Dto2D()


        call leaf%SetSolver(InSolverType=SolverName)
        call leaf%SetUp(NoFacetMode=.true.)

        call leaf%SetMatPara(materialist=materialist,simple=.true.,MaterialID=1)
        
        call leaf%setBC(MPIData=MPIData,dirichlet=.true.,Boundinfo=DBoundlist)
        call leaf%setBC(MPIData=MPIData,neumann=.true.,Boundinfo=NBoundlist)
        call leaf%SetControlPara(OptionalItrTol=100,OptionalTimestep=100,OptionalSimMode=1)
        
        
        !call leaf%Export(Name=trim(project)//"root"//trim(str_id)//".geo")

        ! get soil mesh
        if(present(Soilfile) )then
            sname=trim(Soilfile)
            call soil%ImportPictureName(sname)
            call soil%GetPixcelSize(MPIData)
            call soil%SetColor(sR,sG,sB)
            call soil%GetPixcelByRGB(MPIData,err=5,onlycoord=.true.)

            ! Get Outline (simple outline)
            ! see soil as a box
            call soil%GetSurfaceNode(MPIData,box=.true.)
            call soil%modifySuefaceNode(Mesh=leaf%FEMDomain%Mesh,boolean="diff")
            
            ! Convert SurfaceNod to .geo
            call soil%ExportGeoFile(MPIData,Name=trim(project)//"soil"//trim(str_id)//".geo" )
            
            ! Run Gmsh to convert .geo to .msh
            call soil%ConvertGeo2Msh(MPIData ,Name=trim(project)//"soil"//trim(str_id)//".geo" )
            call soil%ConvertGeo2Inp(MPIData ,Name=trim(project)//"soil"//trim(str_id)//".geo" )
            call soil%ConvertGeo2Mesh(MPIData,Name=trim(project)//"soil"//trim(str_id)//".geo" )
            ! Convert .msh to .scf
            call soil%ConvertMesh2Scf(MPIData,ElementType=ElemType,&
            Name=trim(project)//"soil"//trim(str_id)//".mesh")
            
            call soil%FEMDomain%checkconnectivity(fix=.true.)
            
            call soil%SetSolver(InSolverType=SolverName)
            call soil%SetUp(NoFacetMode=.true.)

            call soil%SetMatPara(materialist=materialist,simple=.true.,MaterialID=2)
            

            ! setup boundary conditions
            call soil%setBC(MPIData=MPIData,dirichlet=.true.,Boundinfo=DBoundlist)
            call soil%setBC(MPIData=MPIData,Neumann=.true.,Boundinfo=NBoundlist)

            call soil%SetControlPara(OptionalItrTol=100,OptionalTimestep=100,OptionalSimMode=1)
            
            !call soil%Export(Name=trim(project)//"soil"//trim(str_id)//".geo")

        endif
        xratio=scalex/leaf%PixcelSize(1)
        yratio=scaley/leaf%PixcelSize(2)
        call soil%SetScale(xratio=xratio,yratio=yratio)
        call leaf%SetScale(xratio=xratio,yratio=yratio)
        call soil%Reverse()
        call leaf%Reverse()
        call leaf%Export(with=soil,Name=trim(project)//"rootandsoil"//trim(str_id)//".scf",regacy=.true.)
        
        

        ! destructor
        call leaf%finalize()
        call soil%finalize()

        
    enddo




    
    
    
    ! Export Object
    !call leaf%FEMDomain%GmshPlotVector(Name="Tutorial/InputData/grass_leaf",step=0,&
    !    withMsh=.true.,FieldName="DispBound",NodeWize=.true.,onlyDirichlet=.true.)
    
    


end subroutine
! #########################################################

! #########################################################
subroutine InitializePrePro(obj,Default)
    class(PreProcessing_),intent(inout)::obj
    logical,optional,intent(in)::Default

    call obj%FEMDomain%Init(Default)
end subroutine
! #########################################################

! #########################################################
subroutine finalizePrePro(obj)
    class(PreProcessing_),intent(inout)::obj

    call obj%FEMDomain%delete()
    
    obj%PictureName = ""
    obj%RGBDataName = ""
    obj%PixcelSizeDataName = ""
    obj%PixcelSize(:)=0
    obj%num_of_pixcel=0
    obj%ColorRGB(:) =0 
end subroutine
! #########################################################


! #########################################################
subroutine ImportPictureName(obj,Name)
    class(PreProcessing_)::obj
    character(*),intent(in)::Name
    obj%PictureName=trim(Name)
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_),intent(inout)          :: MPIData
    character(*),optional,intent(in)   :: Name
    character *30       :: pid
    character *200      :: python_script
    character *200      :: python_buffer
    character *200      :: command
    integer(int32)              :: fh


    call MPIData%GetInfo()

    pid = trim(adjustl(fstring(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




    obj%PixcelSizeDataName=python_buffer
    !print *, trim(python_script)
    
    ! using python script
    ! python imaging library is to be installed.
    fh=MPIData%MyRank+100


    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)//'")'
    !print *, command
    write(fh,'(A)') adjustl(trim(command))
    command = 'python_buffer = open("'//trim(python_buffer)//'","w")'
    write(fh,'(A)') adjustl(trim(command))
    !print *, command
    ! get pixcel size
    command = "rgb_im = img_in.convert('RGB')"
    !print *, command
    write(fh,'(A)') adjustl(trim(command))
    command = "size = rgb_im.size"
    !print *, command
    write(fh,'(A)') adjustl(trim(command))
    command = "print( str(size[0]), ' ',str(size[1])  ) "
    !print *, command
    write(fh,'(A)') adjustl(trim(command))
    
    ! write size
    command = "python_buffer.write( str(size[0]))"
    !print *, command
    write(fh,'(A)') adjustl(trim(command))
    command = "python_buffer.write('\n')"
    !print *, command
    write(fh,'(A)') adjustl(trim(command))
    command = "python_buffer.write( str(size[1]))"
    !print *, command
    write(fh,'(A)') adjustl(trim(command))
    
    ! close
    command = "img_in.close()"
    !print *, command
    write(fh,'(A)') adjustl(trim(command))
    command = "python_buffer.close()"
    !print *, command
    write(fh,'(A)') adjustl(trim(command))
    close(fh)

    command = "python3 "//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(int32)              :: 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 = "python3 "//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(int32),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(int32)              :: 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"
    obj%RGBDataName=python_buffer
    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 = "python3 "//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")
    if(sizeofpc==0)then
        print *, "ERROR :: GetPixcelByRGB >> no such color"
        stop 
    endif
    do i=1,sizeofpc
        read(fh,*)obj%FEMDomain%Mesh%NodCoord(i,1:2)
    enddo
    obj%FEMDomain%Mesh%NodCoord(:,2)=-1.0d0*obj%FEMDomain%Mesh%NodCoord(:,2)
    close(fh)


end subroutine
! #########################################################


! #########################################################
subroutine SetColor(obj,Red,Green,Blue)
    class(PreProcessing_),intent(inout):: obj
    integer(int32),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,convex,division,box)
    class(PreProcessing_),intent(inout):: obj
    class(MPI_),intent(inout)          :: MPIData
    character(*),optional,intent(in)   :: Name
    integer(int32),optional,intent(in)        :: r,NumOfMaxNod,division
    logical,optional,intent(in)        :: convex,box

    character*200   :: python_buffer
    character*20    :: pid
    
    integer(int32),allocatable :: KilledPixcel(:)
    
    integer(int32) :: i,j,k,n,node_id,check_node_id,point_count,fh,MaxNod,dv,itr
    real(real64) :: x_real,y_real,z_real,xmin,xmax,ymin,ymax,dly,dlx,xwidth,ywidth,xm,ym
    real(real64) :: x_real_tr,y_real_tr,z_real_tr,diff_real,max_r,x_tr,y_tr,ymax_tr,ymin_tr
    real(real64),allocatable :: buffer(:,:),NodCoordDiv(:,:),xmaxval(:),ymaxval(:),xminmval(:),yminmval(:)
    ! in case of box
    if(present(box) )then
        if(box .eqv. .true.)then
            print *, "Notice :: obj%GetSurfaceNode as Box"
            xm=minval(obj%FEMDomain%Mesh%NodCoord(:,1))
            ym=minval(obj%FEMDomain%Mesh%NodCoord(:,2)) 
            dv = input(default=10,option=division)        
            xwidth = maxval(obj%FEMDomain%Mesh%NodCoord(:,1))-minval(obj%FEMDomain%Mesh%NodCoord(:,1))
            ywidth = maxval(obj%FEMDomain%Mesh%NodCoord(:,2))-minval(obj%FEMDomain%Mesh%NodCoord(:,2)) 
            dlx = xwidth/dble(dv)
            dly = ywidth/dble(dv)
            xmax=maxval(obj%FEMDomain%Mesh%NodCoord(:,1))
            ymax=maxval(obj%FEMDomain%Mesh%NodCoord(:,2))
            xmin=minval(obj%FEMDomain%Mesh%NodCoord(:,1))
            ymin=minval(obj%FEMDomain%Mesh%NodCoord(:,2))
            allocate(buffer(4*dv,2) )
            buffer(:,:)=0.0d0
            
            do i=1,dv
                buffer(i,1) = dlx*dble(i-1) + xm 
                buffer(i,2) = ymin
            enddo

            do i=1,dv
                buffer(i+dv,1) = xmax 
                buffer(i+dv,2) = dly*dble(i-1) + ym
            enddo

            do i=1,dv
                buffer(i+dv*2,1) = xmax - dlx*dble(i-1)  
                buffer(i+dv*2,2) = ymax
            enddo

            do i=1,dv
                buffer(i+dv*3,1) = xmin 
                buffer(i+dv*3,2) = ymax - dly*dble(i-1)
            enddo

            deallocate(obj%FEMDomain%Mesh%NodCoord)
            allocate(obj%FEMDomain%Mesh%NodCoord(size(buffer,1),size(buffer,2) )  )
            do i=1,size(buffer,1)
                obj%FEMDomain%Mesh%NodCoord(i,:)=buffer(size(buffer,1)-i+1,:)
            enddo
            return
        endif
    endif

    ! in case of convex
    ! in case of convex
    !if(present(convex) )then
    !    if(convex .eqv. .true.)then
    !        xm=minval(obj%FEMDomain%Mesh%NodCoord(:,1))
    !        ym=minval(obj%FEMDomain%Mesh%NodCoord(:,2)) 
    !        ! get discrete points with interval by division
    !        dv = input(default=10,option=division)
    !        xwidth = maxval(obj%FEMDomain%Mesh%NodCoord(:,1))-minval(obj%FEMDomain%Mesh%NodCoord(:,1))
    !        ywidth = maxval(obj%FEMDomain%Mesh%NodCoord(:,2))-minval(obj%FEMDomain%Mesh%NodCoord(:,2)) 
    !        dlx = xwidth/dble(dv)
    !        dly = ywidth/dble(dv)
    !        allocate(buffer(dv*4,2) )
    !        allocate(xmaxval(dv) )
    !        allocate(xmaxval(dv) )
    !        allocate(yminval(dv) )
    !        allocate(yminval(dv) )
!
    !        ! get x-max values
    !        do i=1,dv
    !            ymin=dble(i-1)*ywidth+ym
    !            ymax=dble(i)*ywidth+ym
    !            itr=1
    !            do j=1,size(obj%FEMDomain%Mesh%NodCoord,1)
    !                y_tr=obj%FEMDomain%Mesh%NodCoord(j,2)
    !                if(ymin <= y_tr .and. y_tr <= ymax  )then
    !                    if(itr==1)then
    !                        itr=itr+1
    !                        xmax_tr=obj%FEMDomain%Mesh%NodCoord(j,1)
    !                        xmin_tr=obj%FEMDomain%Mesh%NodCoord(j,1)
    !                    else
    !                        if(xmax_tr<=obj%FEMDomain%Mesh%NodCoord(j,1) )then
    !                            xmax_tr=obj%FEMDomain%Mesh%NodCoord(j,1) ! update y-max
    !                            xmaxval(i)=xmax_tr
    !                        endif
    !                        if(xmin_tr>=obj%FEMDomain%Mesh%NodCoord(j,1) )then
    !                            xmin_tr=obj%FEMDomain%Mesh%NodCoord(j,1) ! update y-max
    !                            xminval(i)=xmin_tr
    !                        endif
    !                    endif
    !                else
    !                    cycle
    !                endif
    !            enddo
    !        enddo
!
    !        do i=1,dv
    !            buffer(i,2)=ymaxval(i)
    !        enddo
!
    !        do i=1,dv
    !            buffer(i+dv*2,dv*3)=yminval(i)
    !        enddo
    !        
    !        do i=1,dv
    !            buffer(i,2)=ymaxval(i)
    !        enddo
!
    !        do i=1,dv
    !            buffer(i+dv*2,dv*3)=yminval(i)
    !        enddo
    !        
!
    !    endif
    !endif
    ! in case of convex
    ! in case of convex


    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
            if(size(obj%FEMDomain%Mesh%NodCoord,1 ) < i )then
                print *, "ERROR"
                exit
            endif
            if(size(buffer,1 ) < point_count )then
                print *, "ERROR"
                exit
            endif
            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(:,:)
    

    ! remove clossing 
    call unwindLine(obj%FEMDomain%Mesh%NodCoord)
    


end subroutine
! #########################################################



! #########################################################
subroutine AssembleSurfaceElement(obj,MPIData,dim,threshold,DelRange,Name)
    class(PreProcessing_),intent(inout):: obj
    class(MPI_),intent(inout)          :: MPIData
    integer(int32),optional,intent(in)        :: dim,threshold,DelRange
    character(*),optional,intent(in)   :: Name
    character*200   :: python_buffer
    character*20    :: pid
    real(real64),allocatable     :: buffer(:,:),r_value(:),line_segment_len(:)
    integer(int32),allocatable     :: checked(:),line_segment(:,:),kill(:)
    real(real64)                 :: x(3),x_tr(3),r_tr,r_ref,a1(2),a2(2),b1(2),b2(2),c1,c2,c3
    real(real64) :: bufvec(2),middle(2)
    integer(int32) :: i,j,k,n,trial_num,id_tr,fh,r_threshold,drange,nn,ys_solved,m


    
    if(present(threshold) )then
        r_threshold=threshold
    else
        r_threshold=5
    endif

    if(present(DelRange) )then
        drange=DelRange
    else
        drange=5
    endif

    if( present(dim) .and. dim/=2 )then
        call MPIData%End()
        print *, "AssembleSurfaceElement :: >> only 2-D is available."
        stop   
    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)



    ! modifier to remove invalid surface nodes
    ! remove solitary island
    ! for 2-D cases
    ! get median of line segments
    n=size(obj%FEMDomain%Mesh%NodCoord,1)
    allocate(line_segment(n,3),line_segment_len(n) ) 
    line_segment(:,:)=0
    do i=1,n-1
        line_segment(i,1) = i
        line_segment(i,2) = i+1
    enddo
    line_segment(n,1) = n
    line_segment(n,2) = 1

    do i=1,n
        a1(1:2)=obj%FEMDomain%Mesh%NodCoord( line_segment(i,1) ,1:2 )
        a2(1:2)=obj%FEMDomain%Mesh%NodCoord( line_segment(i,2) ,1:2 )
        line_segment_len(i)=dsqrt(dot_product(a2-a1,a2-a1)  )
    enddo

    ! write a operation to remove invalid nodes.


    ! remove crossing surface
    

    do i=1,n-2
        a1(1:2)=obj%FEMDomain%Mesh%NodCoord( line_segment(i,1) ,1:2 )
        a2(1:2)=obj%FEMDomain%Mesh%NodCoord( line_segment(i,2) ,1:2 )-a1(1:2)
        ys_solved = 0
        do j=i+2,n
            b1(1:2)=obj%FEMDomain%Mesh%NodCoord( line_segment(j,1) ,1:2 )-a1(1:2)
            b2(1:2)=obj%FEMDomain%Mesh%NodCoord( line_segment(j,2) ,1:2 )-a1(1:2)
            ! detect crossing by cross product
            c1=( a2(1)*b1(2) - a2(2)*b1(1))*( a2(1)*b2(2) - a2(2)*b2(1))
            if(c1>0.0d0)then
                cycle
            else
                b1(1:2)=b1(1:2)+a1(1:2)
                b2(1:2)=b2(1:2)+a1(1:2)
                a2(1:2)=a2(1:2)+a1(1:2)

                b2(1:2)=b2(1:2)-b1(1:2)
                a1(1:2)=a1(1:2)-b1(1:2)
                a2(1:2)=a2(1:2)-b1(1:2)
                c2=( b2(1)*a1(2)-b2(2)*a1(1) )*( b2(1)*a2(2)-b2(2)*a2(1) )
                if(c2>0.0d0)then
                    cycle
                else
                    ! crossed
                    
                    nn= line_segment(j,1)
                    k=i
                    
                    do 
                        if( line_segment(k,2) >= line_segment(nn,1) )then
                            ys_solved=1
                            exit
                        endif
                        if( abs(line_segment(k,2) - line_segment(nn,1))<=1 )then
                            ys_solved=1
                            exit
                        endif
                        
                        print *, "surface-line segments are Crossed >> modification ",line_segment(k,2) ," to ",line_segment(nn,1)
                        
                        bufvec(1:2)=obj%FEMDomain%Mesh%NodCoord(  line_segment(k,2) ,1:2)
                        obj%FEMDomain%Mesh%NodCoord( line_segment(k,2) ,1:2 ) = &
                            obj%FEMDomain%Mesh%NodCoord( line_segment(nn,1) ,1:2 )
                        obj%FEMDomain%Mesh%NodCoord( line_segment(nn,1) ,1:2 ) = bufvec(1:2)
                        
                        if(  line_segment(nn,1) - line_segment(k,2)  <= 1)then
                            ys_solved=1
                            exit
                        endif

                        nn=nn-1
                        k=k+1

                    enddo
                endif
            endif

            if(ys_solved == 1)then
                exit
            endif
        enddo
    enddo

    ! reverse
    nn=size(obj%FEMDomain%Mesh%NodCoord,1)
    do i=1,size(obj%FEMDomain%Mesh%NodCoord,1)
        if(nn <= i )then
            exit
        endif
        bufvec(1:2)=obj%FEMDomain%Mesh%NodCoord(i,1:2)
        obj%FEMDomain%Mesh%NodCoord(i,1:2)=obj%FEMDomain%Mesh%NodCoord(nn,1:2)
        obj%FEMDomain%Mesh%NodCoord(nn,1:2)=bufvec(1:2)
        nn=nn-1
    enddo

    ! kill invalid nodes
    nn=size(obj%FEMDomain%Mesh%NodCoord,1)
    allocate(kill(size(obj%FEMDomain%Mesh%NodCoord,1) ) )
    kill(:)=0
    do i=1,nn-1
        a1(1:2)=obj%FEMDomain%Mesh%NodCoord( line_segment(i,1) ,1:2 ) !11
        a2(1:2)=obj%FEMDomain%Mesh%NodCoord( line_segment(i,2) ,1:2 ) !9
        b1(1:2)=obj%FEMDomain%Mesh%NodCoord( line_segment(i+1,2) ,1:2 ) !10
        middle(1:2)=0.50d0*a2(1:2)+0.50d0*a1(1:2) !11-9
        if( dot_product(middle-b1,middle-b1) == 0.0d0 )then
            ! invalid node
            kill( line_segment(i+1,2) ) = 1
            print *, line_segment(i,2) , " will be killed."
        endif
    enddo

    if(allocated(buffer) )then
        deallocate(buffer)
    endif

    allocate(buffer( nn-sum(kill),2 ) )
    k=0
    do i=1,nn
        if(kill(i)==1 )then
            cycle
        else
            k=k+1
            buffer(k,1:2)=obj%FEMDomain%Mesh%NodCoord(i,1:2)
        endif
    enddo
    deallocate(obj%FEMDomain%Mesh%NodCoord)
    n=size(buffer,1)
    m=size(buffer,2)
    allocate(obj%FEMDomain%Mesh%NodCoord(n,m)  )
    do i=1,n
        obj%FEMDomain%Mesh%NodCoord(i,1:2)=buffer(i,1:2)
    enddo






end subroutine
! #########################################################



! #########################################################
subroutine ReduceSize(obj,MPIData,interval,Name,auto,curvetol)
    class(PreProcessing_),intent(inout):: obj
    class(MPI_),intent(inout)          :: MPIData
    character(*),optional,intent(in)    :: Name
    integer(int32),optional,intent(in) :: interval
    character*200   :: python_buffer
    character*20    :: pid
    real(real64),allocatable:: buffer(:,:)
    integer(int32),allocatable:: chosen(:),kill(:)
    logical,optional,intent(in) :: auto
    integer(int32) :: i,j,k,n,m,fh,itr,id1,id2,id3,killcount,killcount_b
    real(real64),optional,intent(in) :: curvetol
    real(real64) :: x1(2),x2(2),x3(2),x4(2),dp1,dp2,tol

    tol=input(default=0.010d0,option=curvetol)
    if(present(auto) )then
        if(auto .eqv. .true.)then
            
            n=size(obj%FEMDomain%Mesh%NodCoord,1)  
            m=size(obj%FEMDomain%Mesh%NodCoord,2)  
            allocate(kill(n) )
            kill(:)=0
            do i=1,n-2
                if(kill(i)==1 )then
                    cycle
                endif
                do 
                    x1(1:2)=obj%FEMDomain%Mesh%NodCoord(i  ,1:2)    
                    x2(1:2)=obj%FEMDomain%Mesh%NodCoord(i+1,1:2)  
                    x3(1:2)=obj%FEMDomain%Mesh%NodCoord(i+2,1:2)  
                    x4(1:2)=0.50d0*x1(1:2)+0.50d0*x2(1:2)
                    dp1=dot_product(x3-x4,x3-x4)
                    dp2=dot_product(x1-x4,x1-x4)
                    if(dp1/dp2 < tol )then
                        kill(i+1)=1
                        cycle
                    else
                        exit
                    endif
                enddo
            enddo

            ! kill nodes
            allocate( buffer( n-sum(kill) ,2) )
            itr=1
            do i=1,n
                if(kill(i)==1 )then
                    buffer(itr,1:2)=obj%FEMDomain%Mesh%NodCoord(i  ,1:2)
                    itr=itr+1
                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(:,:)
            return
        endif
    endif


    n=size(obj%FEMDomain%Mesh%NodCoord,1)  
    m=size(obj%FEMDomain%Mesh%NodCoord,2)   
    allocate(chosen(n) )
    chosen(:)=0
    k=0
    do i=1,n
        k=k+1
        if(k==interval)then
            chosen(i)=1
            k=0
        else
            cycle
        endif 
    enddo

    allocate(buffer( sum(chosen),3 ))

    k=0
    do i=1,n
        if(chosen(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(int32) :: i,j,k,n,fh,xsize
    real(real64) :: 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=trim(Name)
        
    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,clmin,clmax)
    class(PreProcessing_),intent(inout):: obj
    class(MPI_),intent(inout)          :: MPIData
    character(*),optional,intent(in)    :: Name
    real(real64),optional,intent(in) :: clmin,clmax
    character*200   :: python_buffer
    character*200   :: command
    character*20    :: pid
    integer(int32) :: i,j,k,n,fh
    real(real64) :: cmin,cmax

    write(pid,*) MPIData%MyRank

    cmin=input(default=100.0d0,option=clmin)
    cmax=input(default=100000.0d0,option=clmax)
    fh=MPIData%MyRank+10
    python_buffer=" "
    command  = " "
    python_buffer="GetSurface_pid_"//trim(adjustl(pid))//".geo"
    if(present(Name) )then
        python_buffer=trim(Name)
    endif
    command="gmsh "//trim(python_buffer)//" -2 -algo del2d -clmin "//trim(fstring(cmin))&
        //" -clmax "//trim(fstring(cmax))

    writE(*,'(A)') trim(command)
    
    call execute_command_line(command)

    !call execute_command_line("sh ./MakeMesh.sh")
    


end subroutine
! #########################################################


! #########################################################
subroutine ConvertGeo2VTK(obj,MPIData,Name,clmin,clmax)
    class(PreProcessing_),intent(inout):: obj
    class(MPI_),intent(inout)          :: MPIData
    character(*),optional,intent(in)    :: Name
    real(real64),optional,intent(in) :: clmin,clmax
    character*200   :: python_buffer
    character*200   :: command
    character*20    :: pid
    integer(int32) :: i,j,k,n,fh
    real(real64) :: cmin,cmax

    write(pid,*) MPIData%MyRank

    cmin=input(default=100.0d0,option=clmin)
    cmax=input(default=100000.0d0,option=clmax)
    fh=MPIData%MyRank+10
    python_buffer=" "
    command  = " "
    python_buffer="GetSurface_pid_"//trim(adjustl(pid))//".geo"
    if(present(Name) )then
        python_buffer=trim(Name)
    endif
    command="gmsh "//trim(python_buffer)//" -2 -algo del2d -clmin "//trim(fstring(cmin))&
        //" -clmax "//trim(fstring(cmax))//" -format vtk"

    writE(*,'(A)') trim(command)
    
    call execute_command_line(command)

    !call execute_command_line("sh ./MakeMesh.sh")
    


end subroutine
! #########################################################


! #########################################################
subroutine ConvertGeo2Inp(obj,MPIData,Name,clmin,clmax)
    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(int32) :: i,j,k,n,fh
    real(real64) :: cmin,cmax
    real(real64),optional,intent(in) :: clmin,clmax

    write(pid,*) MPIData%MyRank
    cmin=input(default=100.0d0,option=clmin)
    cmax=input(default=100000.0d0,option=clmax)

    fh=MPIData%MyRank+10
    python_buffer=" "
    command  = " "
    python_buffer="GetSurface_pid_"//trim(adjustl(pid))//".geo"
    if(present(Name) )then
        python_buffer=trim(Name)
    endif
    !command="gmsh "//trim(python_buffer)//" -2 -algo del2d -clmin 100 -clmax 100000 -format inp" 
    command="gmsh "//trim(python_buffer)//" -2 -algo del2d -clmin "//trim(fstring(cmin))//&
        " -clmax "//trim(fstring(cmax))//" -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,clmin,clmax)
    class(PreProcessing_),intent(inout):: obj
    class(MPI_),intent(inout)          :: MPIData
    integer(int32),optional,intent(in) :: SizePara
    character(*),optional,intent(in)    :: Name
    character*200   :: python_buffer
    character*200   :: command
    character*20    :: pid,a
    integer(int32) :: i,j,k,n,fh,sp
    real(real64) :: cmin,cmax
    real(real64),optional,intent(in) :: clmin,clmax
    
    cmin=input(default=100.0d0,option=clmin)
    cmax=input(default=100000.0d0,option=clmax)

    if(present(SizePara) )then
        sp=SizePara
    else
        sp=100
    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=trim(Name)
    endif
    !command="gmsh "//trim(python_buffer)//" -2 -algo del2d -clmin"//trim(a)//" -clmax 100000  -format mesh" 
    command="gmsh "//trim(python_buffer)//" -2 -algo del2d -clmin "//trim(fstring(cmin))&
        //" -clmax "//trim(fstring(cmax))//" -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(int32),allocatable :: elem1(:),surface_nod(:)
    integer(int32) :: i,j,k,n,n1,n2,fh,a,nm,mm,nod_num,nn,elem_num,surf_num
    integer(int32) :: elem_num_all,n3,n4,n5,n6,n7,elemnod_num,startfrom
    real(real64) :: 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_),optional,intent(inout)          :: MPIData
    type(Mesh_) :: tobj
    character(*),optional,intent(in) :: Name
    character*200   :: python_buffer
    character*200   :: command,infile,outfile
    character(*),optional,intent(in) :: ElementType
    character*20    :: pid
    character*200 MeshFormat
	character*14 EndMeshFormat
	character*6  Nodes
	character*200  EndNodes,Elements
    character*12  EndElements	
    integer(int32),allocatable :: elem1(:),surface_nod(:),triangle(:,:),devide_line(:,:),buffer(:,:)
    integer(int32) :: i,j,k,n,n1,n2,fh,a,nm,mm,nod_num,nn,elem_num,surf_num,l,numnum
    integer(int32) :: elem_num_all,n3,n4,n5,n6,n7,elemnod_num,startfrom,node1,node2,tr1,tr2
    real(real64) :: re1,re2
    

    print *,  "ConvertMesh2Scf >>>> only for 2D"
    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
    ! ======================================================
    


    ! ======================================================
    numnum=input(default=1,option=MPIData%MyRank)
    write(pid,*) numnum
    fh=input(default=1,option=MPIData%MyRank+10)
    infile="GetSurface_pid_"//trim(adjustl(pid))//".mesh"
    outfile = "GetSurface_pid_"//trim(adjustl(pid))//".scf"
    if(present(Name) )then
        infile  = Name
        outfile = Name//".scf"
        
    endif
	
    open(fh,file=infile,status="old")
    print *, "Opening ",trim(infile)
	! ======================================================
    
    ! ======================================================
	!read file to get nod and elem number
    read(fh,*)MeshFormat
    if( trim(adjustl(MeshFormat) )/= "MeshVersionFormatted")then
        print *, "ConvertMesh2Scf ERROR :: ",MeshFormat
    endif
    read(fh,*) MeshFormat
    read(fh,*) mm
    do 
        MeshFormat=" "
        read(fh,*) MeshFormat
        if( trim(adjustl(MeshFormat) ) == "Vertices" )then
            print *, "ConvertMesh2Scf reading ",trim(adjustl(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,2) )
            obj%FEMDomain%Mesh%NodCoord(:,:)=0.0d0
	        do i=1,nod_num
	        	read(fh,*) obj%FEMDomain%Mesh%NodCoord(i,1:2)
            enddo

            cycle
            ! ======================================================
        elseif( trim(adjustl(MeshFormat) ) == "Quadrilaterals" )then
            ! ======================================================
            print *, "ConvertMesh2Scf reading ",trim(adjustl(MeshFormat) )
            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
            exit
            ! ======================================================
        elseif( trim(adjustl(MeshFormat) ) == "Triangles" )then
            
            ! ======================================================
            print *, "ConvertMesh2Scf reading ",trim(adjustl(MeshFormat) )
            read(fh,*)mm
            allocate(tobj%ElemNod(mm,3) )

            do i=1,mm
                read(fh,*) tobj%ElemNod(i,1:3)
            enddo
            ! ======================================================
        else
            print *, "ConvertMesh2Scf Skipped",trim(adjustl(MeshFormat))
            if(trim(adjustl(MeshFormat)) == "End")then
                exit
            endif
            read(fh,*)mm
            do i=1,mm
                read(fh,*) n
            enddo
        endif
        
    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 
    if(.not. allocated(tobj%ElemNod) )then
        print *, "No triangles"
        return
    endif

    allocate( tobj%NodCoord(size(obj%FEMDomain%Mesh%NodCoord,1  ),2  ) )
    tobj%NodCoord(:,1:2)=obj%FEMDomain%Mesh%NodCoord(:,1:2)

    call tobj%convertMeshType(option="convertTriangleToRectangular")
        
    if(allocated(obj%FEMDomain%Mesh%ElemNod)  )then
        print *, "triangular and rectangurar => ignore triangular"
        return
    else
        print *, "triangular => converted."
        allocate( obj%FEMDomain%Mesh%ElemNod( sizE(tobj%ElemNod,1),1:4  )   )
        deallocate(obj%FEMDomain%Mesh%NodCoord )
        allocate(obj%FEMDomain%Mesh%NodCoord(size(tobj%NodCoord,1),size(tobj%NodCoord,2) ) )
        obj%FEMDomain%Mesh%NodCoord(:,:)=tobj%NodCoord(:,:)
        obj%FEMDomain%Mesh%ElemNod(:,:)=tobj%ElemNod(:,:)

    endif

    print *, "surface information is updated"
    call obj%FEMDomain%Mesh%GetSurface()

    return

    
    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
    


end subroutine
! #########################################################









! #########################################################
subroutine ExportPreProcessing(obj,MPIData,FileName,MeshDimension,Name,regacy,with)
    class(PreProcessing_),intent(inout):: obj
    class(PreProcessing_),optional,intent(inout):: with
    class(MPI_),optional,intent(inout) :: MPIData
    character*200,optional,intent(in)  :: FileName
    character(*),optional,intent(in)   :: Name
    integer(int32),optional,intent(in) :: MeshDimension
    logical,optional,intent(in)::regacy
    character*200   :: python_buffer
    character*200   :: command
    character*20    :: pid
    integer(int32) :: 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,Name=Name,regacy=regacy&
    ,with=with%FEMDomain)
    
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(real64),allocatable,optional,intent(inout)::MatPara(:,:)
    real(real64),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(:,:)=1.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,xratio,yratio)
    class(PreProcessing_),intent(inout)::obj
    real(real64),optional,intent(in)::scalex,scaley,scalez
    real(real64),optional,intent(in)::picscalex,picscaley,picscalez,xratio,yratio
    real(real64) :: lx,ly,lz
    integer(int32) :: i

    if(present(xratio) )then
        do i=1,size(obj%FEMDomain%Mesh%NodCoord,1)
            obj%FEMDomain%Mesh%NodCoord(i,1)=xratio*obj%FEMDomain%Mesh%NodCoord(i,1)
        enddo
    endif
    if(present(yratio) )then
        do i=1,size(obj%FEMDomain%Mesh%NodCoord,1)
            obj%FEMDomain%Mesh%NodCoord(i,2)=yratio*obj%FEMDomain%Mesh%NodCoord(i,2)
        enddo
    endif
    if(present(yratio) .or. present(xratio))then
        return
    endif

    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,BoundInfo,MPIData)
    class(PreProcessing_),intent(inout)::obj
    type(preprocessing_) :: DBC
    class(Dictionary_),optional,intent(in) :: BoundInfo
    class(MPI_),optional,intent(inout) :: MPIData
    real(real64),optional,intent(in)::xmin,xmax
    real(real64),optional,intent(in)::ymin,ymax
    real(real64) :: x_min,x_max
    real(real64) :: y_min,y_max
    real(real64),optional,intent(in)::zmin,zmax
    real(real64),optional,intent(in)::tmin,tmax
    logical,optional,intent(in)::Dirichlet,Neumann,Initial
    integer(int32),optional,intent(in)::NumOfValPerNod,val_id
    real(real64),optional,intent(in)::val
    integer(int32) :: i,j,n,k,l

    if(present(BoundInfo) )then
        if(.not.allocated(BoundInfo%Dictionary))then
            return
        endif

        do i=1,BoundInfo%sizeof()
            ! list of boundary conditions is given as figures
            print *, "now dbc"
            if(.not. present(MPIData) )then
                print *, "ERROR :: setBC :: MPIData should be imported."
                return
            endif
            print *,trim(BoundInfo%content(i) )

            call DBC%ImportPictureName( trim(BoundInfo%content(i) ) )
            call DBC%GetPixcelSize(MPIData, name="DBC")
            call DBC%SetColor(BoundInfo%IntList(i,1),&
                BoundInfo%IntList(i,2),BoundInfo%IntList(i,3))
            
            call DBC%GetPixcelByRGB(MPIData,err=5,onlycoord=.true.)
            
            x_min=minval(DBC%FEMDomain%Mesh%NodCoord(:,1) )
            x_max=maxval(DBC%FEMDomain%Mesh%NodCoord(:,1) )
            y_min=minval(DBC%FEMDomain%Mesh%NodCoord(:,2) )
            y_max=maxval(DBC%FEMDomain%Mesh%NodCoord(:,2) )

            ! debug
            !print *, minval(obj%FEMDomain%Mesh%NodCoord(:,1)),&
            !maxval(obj%FEMDomain%Mesh%NodCoord(:,1)),&
            !minval(obj%FEMDomain%Mesh%NodCoord(:,2)),&
            !maxval(obj%FEMDomain%Mesh%NodCoord(:,2))
            n=size(obj%FEMDomain%Mesh%NodCoord,2)
            call obj%setBC(Dirichlet=.true.,xmin=x_min,xmax=x_max,ymin=y_min,ymax=y_max,&
            val_id=BoundInfo%intvalue(i),val=BoundInfo%realvalue(i),NumOfValPerNod=n )
            
            print *, "boundary condition ID : ",i
            call DBC%finalize()
        enddo
        return
    endif

    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)
            print *, "boundary conditions are added."
            if(.not. allocated(obj%FEMDomain%Boundary%DBoundNum) )then
                n=size(obj%FEMDomain%Boundary%DBoundNodID,2)
                allocate(obj%FEMDomain%Boundary%DBoundNum(n) )
                print *, "caution .not. allocated(obj%FEMDomain%Boundary%DBoundNum "
            endif
            do i=1,size(obj%FEMDomain%Boundary%DBoundNum)
                k=countif(Array=obj%FEMDomain%Boundary%DBoundNodID(:,i),Value=-1,notEqual=.true.)
                obj%FEMDomain%Boundary%DBoundNum(i)=k
            enddo

        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(int32),optional,intent(in)::NumOfValue
    logical :: DB,NB,IC
    integer(int32) :: 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(int32) :: 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(real64),allocatable::buffer(:,:)
    integer(int32) :: 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)
    if(allocated(obj%FEMDomain%Mesh%ElemMat)  )then
        if(size(obj%FEMDomain%Mesh%ElemMat)/=size(obj%FEMDomain%Mesh%ElemNod,1) )then
            deallocate(obj%FEMDomain%Mesh%ElemMat)
        endif
    endif
    if(.not.allocated(obj%FEMDomain%Mesh%ElemMat) )then
        n=size(obj%FEMDomain%Mesh%ElemNod,1)
        allocate(obj%FEMDomain%Mesh%ElemMat(n) )
        obj%FEMDomain%Mesh%ElemMat(:)=1
    endif
end subroutine
!##################################################


!##################################################
subroutine Convert2Dto3D(obj,Thickness,division)
    class(PreProcessing_),intent(inout)::obj
    real(real64),allocatable::buffer(:,:)
    real(real64),optional,intent(in)::Thickness
    integer(int32),optional,intent(in)::division
    real(real64) :: Tn
    integer(int32) :: 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(real64),optional,intent(in)::OptionalTol
    integer(int32),optional,intent(in)::OptionalSimMode,OptionalItrTol,OptionalTimestep
    
    call SetControlPara(obj%FEMDomain%ControlPara,OptionalTol,OptionalItrTol,OptionalTimestep,OptionalSimMode)
end subroutine
!##################################################

!##################################################
subroutine SetMatParaPreProcessing(obj,MaterialID,ParameterID,Val,materialist,simple)
    class(PreProcessing_),intent(inout)::obj
    class(Dictionary_),optional,intent(inout) :: materialist
    integer(int32),optional,intent(in)::MaterialID,ParameterID
    real(real64),optional,intent(in)::Val
    logical,optional,intent(in) :: simple
    integer(int32) ::i,n,m,p(2),mm

    if(present(materialist) )then
        ! import material information from list
        print *, "total ",materialist%sizeof()," materials are imported."
        n=materialist%sizeof()
        m=size(materialist%Dictionary(1)%Realist)
        if(allocated(obj%FEMDomain%MaterialProp%MatPara) )then
            deallocate(obj%FEMDomain%MaterialProp%MatPara)
        endif
        allocate(obj%FEMDomain%MaterialProp%MatPara(n,m) )
        do i=1,materialist%sizeof()
            obj%FEMDomain%MaterialProp%MatPara(i,:)=materialist%Dictionary(i)%Realist(:)
        enddo
    endif

    if(present(simple) )then
        if(simple .eqv. .true.)then
            if(.not.allocated(obj%FEMDomain%Mesh%ElemMat) )then
                n=size(obj%FEMDomain%Mesh%ElemNod,1)
                allocate(obj%FEMDomain%Mesh%ElemMat(n) )
            else
                deallocate(obj%FEMDomain%Mesh%ElemMat)
                n=size(obj%FEMDomain%Mesh%ElemNod,1)
                allocate(obj%FEMDomain%Mesh%ElemMat(n) )
            endif
            obj%FEMDomain%Mesh%ElemMat(:)=input(default=1,option=MaterialID)
        endif
        return
    endif



    if(.not.allocated(obj%FEMDomain%MaterialProp%MatPara) )then
        allocate(obj%FEMDomain%MaterialProp%MatPara(MaterialID,ParameterID) )
        obj%FEMDomain%MaterialProp%MatPara(MaterialID,ParameterID)=val
        return
    else
        mm=size(obj%FEMDomain%MaterialProp%MatPara,1)
        if(size(obj%FEMDomain%MaterialProp%MatPara,1)<MaterialID)then
            do i=1,MaterialID-size(obj%FEMDomain%MaterialProp%MatPara,1)
                call extendArray(obj%FEMDomain%MaterialProp%MatPara,extend1stColumn=.true.)
                obj%FEMDomain%MaterialProp%MatPara(i+mm,:)=0.0d0
            enddo
        endif

        mm=size(obj%FEMDomain%MaterialProp%MatPara,2)
        if(size(obj%FEMDomain%MaterialProp%MatPara,2)<ParameterID)then
            do i=1,ParameterID-size(obj%FEMDomain%MaterialProp%MatPara,2)
                call extendArray(obj%FEMDomain%MaterialProp%MatPara,extend2ndColumn=.true.)
                obj%FEMDomain%MaterialProp%MatPara(:,i+mm)=0.0d0
            enddo
        endif
        obj%FEMDomain%MaterialProp%MatPara(MaterialID,ParameterID)=val
    endif

!    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(real64),optional,intent(in)::xmin,xmax
    real(real64),optional,intent(in)::ymin,ymax
    real(real64),optional,intent(in)::zmin,zmax
    real(real64),optional,intent(in)::tmin,tmax
    integer(int32),optional,intent(in)::MaterialID
    integer(int32) :: 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(real64),optional,intent(in) :: Radius,XSize,YSize,ZSize,Xloc,Yloc,Zloc
    integer(int32) :: 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(int32),optional,intent(in) :: XDiv,Ydic,Zdiv
    real(real64) :: ground_level
    real(real64),allocatable::RSInterface(:,:),xmin,xmax,ymin,ymax,NodCoord(:,:)
    integer(int32) :: ground_surface_id,n,m,itr,k,i,j,buf(2)
    integer(int32) :: NodeNum,DimNum,ElemNum,ElemNodNum,startnode,endnode,newnodenum
    integer(int32),allocatable::RSIElemID(:),RSINodeID(:),RSIElemNod(:,:),AvailFE(:)
    integer(int32),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


            
            if(allocated(obj%FEMDomain%Mesh%ElemMat)  )then
                if(size(obj%FEMDomain%Mesh%ElemMat)/=size(obj%FEMDomain%Mesh%ElemNod,1) )then
                    deallocate(obj%FEMDomain%Mesh%ElemMat)
                endif
            endif
            if(.not.allocated(obj%FEMDomain%Mesh%ElemMat) )then
                n=size(obj%FEMDomain%Mesh%ElemNod,1)
                allocate(obj%FEMDomain%Mesh%ElemMat(n) )
                obj%FEMDomain%Mesh%ElemMat(:)=1
            endif

        
            
            





            !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,withNeumannBC,withDirichletBC,withMaterial)
    class(PreProcessing_),intent(inout)::obj
    character(*),optional,intent(in):: Name
    integer(int32),optional,intent(in):: Step
	logical,optional,intent(in)::withNeumannBC,withDirichletBC,withMaterial

    integer(int32) :: stp


    if(present(Step) )then
        stp=step
    else
        stp=0
    endif

    call GmshPlotMesh(obj%FEMDomain,OptionalStep=stp,Name=trim(Name),withNeumannBC=withNeumannBC,&
        withDirichletBC=withDirichletBC,withMaterial=withMaterial)

end subroutine
!##################################################


!##################################################
subroutine meshingPreProcessing(obj)
    class(PreProcessing_),intent(inout)::obj

    call obj%FEMDomain%meshing()
end subroutine
!##################################################

!##################################################
subroutine importPixcelAsNodePreProcessing(obj,interval)
    class(PreProcessing_),intent(inout)::Obj
    integer(int32),optional,intent(in):: interval
    integer(int32) ::i,j, n,k,l,m
    real(real64),allocatable :: random(:)
    real(real64),allocatable :: NewNodCoord(:,:)
    !integer(int32) :: fh1,fh2,xsize,ysize,final_size,interval_,i,j,k,xpixcel,ypixcel

    m=size(obj%FEMDomain%Mesh%NodCoord,1)
    n=input(default=1,option=interval)
    k=size(obj%FEMDomain%Mesh%NodCoord,1)
    allocate(random(k/n) )
    call random_number(random)
    allocate(NewNodCoord(k/n,2) )

    random(:)=random(:)*dble(m)

    do i=1,size(NewNodCoord,1)
        NewNodCoord(i,1:2)=obj%FEMDomain%Mesh%NodCoord(int(random(i) ),1:2)
    enddo
    deallocate(obj%FEMDomain%Mesh%NodCoord)
    allocate(obj%FEMDomain%Mesh%NodCoord(size(NewNodCoord,1),size(NewNodCoord,2) ) )
    obj%FEMDomain%Mesh%NodCoord(:,:)=NewNodCoord(:,:)
    
    return
!    fh1=10
!    open(fh1,file=trim(obj%PixcelSizeDataName),status="old" )
!
!    read(fh1,*) xsize
!    read(fh1,*) ysize  
!    close(fh1)
!
!    fh2=10
!    open(fh2,file=trim(obj%RGBDataName),status="old" )
!    if(allocated(obj%FEMDomain%Mesh%NodCoord) )then
!        deallocate(obj%FEMDomain%Mesh%NodCoord)
!    endif
!    interval_=input(default=1,option=interval)
!    final_size=xsize*ysize/interval_-1
!    allocate(obj%FEMDomain%Mesh%NodCoord(final_size,2 ) )
!    k=1
!    j=1
!    do i=1,xsize*ysize
!        print *, i,"/",xsize*ysize
!        if(j>final_size)then
!            exit
!        endif
!        read(fh2,*) xpixcel,ypixcel
!        obj%FEMDomain%Mesh%NodCoord(j,1)=dble(xpixcel)
!        obj%FEMDomain%Mesh%NodCoord(j,2)=dble(ypixcel)
!        k=k+1
!        if(k==interval_)then
!            k=1
!            j=j+1
!        endif
!    enddo
!    close(fh2)


end subroutine
!##################################################


!##################################################
subroutine removeBoundaryConditionPrePro(obj,Dirichlet,Neumann,Initial)
    class(PreProcessing_),intent(inout)::obj
    logical,optional,intent(in)::Dirichlet,Neumann,Initial
    integer(int32) :: i,j,n


    if(present(Dirichlet) )then
        if(Dirichlet .eqv. .true.)then
            
            call removeDBoundCondition(obj%FEMDomain)
            return

        endif
    endif    
    
    if(present(Neumann) )then
        if(Neumann .eqv. .true.)then
            
            call removeNBoundCondition(obj%FEMDomain)
            
            return
        endif
    endif
    
    if(present(Initial) )then
        if(Initial .eqv. .true.)then
            
            call removeTBoundCondition(obj%FEMDomain)
            
            return
        endif
    endif    
end subroutine
!##################################################



!##################################################
subroutine importPreProcessing(obj,Name,FileHandle,Mesh)
    class(PreProcessing_),intent(inout)::obj
    type(Mesh_),optional,intent(in)::Mesh
    character(*),optional,intent(in)::Name
    integer(int32),optional,intent(in)::FileHandle

    if(present(Mesh) )then
        call obj%FEMDomain%import(Mesh=Mesh)
        return
    endif
    call obj%FEMDomain%import(OptionalProjectName=Name,FileHandle=FileHandle)

end subroutine
!##################################################


!##################################################
subroutine modifySuefaceNodePrepro(obj,Mesh,boolean)
    class(PreProcessing_),intent(inout) :: obj
    class(Mesh_),intent(inout) :: Mesh
    character(*),intent(in) :: boolean
    real(real64),allocatable :: surfacenod_m(:,:), surfacenod(:,:),buffer(:,:),surf_nod_buffer(:,:)
    integer(int32) :: i,j,n,itr,cross1,cross2,cross3,cross4,end1,end2,cases
    integer(int32),allocatable :: in_out(:), in_out_m(:)
    integer(int32) :: s(4),s_m(4),non
    real(real64) :: xmax,ymax,xmin,ymin,x_tr,y_tr,direct,direct_m
    real(real64) :: xmax_m,ymax_m,xmin_m,ymin_m,end1_m,end2_m,xe1,xe2,ye1,ye2

    if(boolean == "diff" .or.boolean == "Diff"  )then
        cross1=0
        ! only for 2D
        if(.not. allocated(Mesh%SurfaceLine2D) )then
            call Mesh%GetSurface()
        endif

        ! Only for box-shaped soil and roots:
        ! It should be boolean operation
        n=size(Mesh%SurfaceLine2D)
        xmin=minval(obj%FEMDomain%Mesh%NodCoord(:,1))
        ymin=minval(obj%FEMDomain%Mesh%NodCoord(:,2))
        xmax=maxval(obj%FEMDomain%Mesh%NodCoord(:,1))
        ymax=maxval(obj%FEMDomain%Mesh%NodCoord(:,2))
        xmin_m=minval(Mesh%NodCoord(:,1))
        ymin_m=minval(Mesh%NodCoord(:,2))
        xmax_m=maxval(Mesh%NodCoord(:,1))
        ymax_m=maxval(Mesh%NodCoord(:,2))
        
        allocate(in_out(size(Mesh%SurfaceLine2D,1) ) )
        in_out(:)=0
        itr=0
        end1=1
        end2=n
        do i=1,n
            x_tr=Mesh%NodCoord( Mesh%SurfaceLine2D(i) ,1)
            y_tr=Mesh%NodCoord( Mesh%SurfaceLine2D(i) ,2)
            if( xmin <= x_tr .and. x_tr<= xmax )then
                if( ymin <= y_tr .and. y_tr<= ymax )then
                    ! in 
                    itr=itr+1
                    in_out(i)=1
                else
                    ! out
                    in_out(i)=0
                    if( ymin > y_tr  )then
                        cross1=1
                    else
                        cross1=3
                    endif
                endif
            else
                !out
                in_out(i)=0
                cross1=2
                if( xmin > x_tr  )then
                    cross1=4
                else
                    cross1=2
                endif
            endif
        enddo
        allocate(surfacenod(itr,2 ))
        
        ! only for roots surrounded by soils
        ! detect two ends
        do i=1,n-1
            if(in_out(i)==0 .and. in_out(i+1)==1 )then
                end1=i
            endif
            if(in_out(i)==1 .and. in_out(i+1)==0 )then
                end2=i
            endif
        enddo    

        itr=0
        do i=1,n
            x_tr=Mesh%NodCoord( Mesh%SurfaceLine2D(i) ,1)
            y_tr=Mesh%NodCoord( Mesh%SurfaceLine2D(i) ,2)
            if( xmin <= x_tr .and. x_tr<= xmax )then
                if( ymin <= y_tr .and. y_tr<= ymax )then
                    itr=itr+1       
                    surfacenod(itr,1)=x_tr
                    surfacenod(itr,2)=y_tr
                else
                    cycle
                endif
            else
                cycle
            endif
        enddo

        ! remove overlapped soil surface
        !allocate( in_out_m(size(obj%FEMDomain%Mesh%NodCoord,1) ))
        !in_out_m(:)=0
        !itr=0
        !do i=1,size(in_out_m)
        !    x_tr=obj%FEMDomain%Mesh%NodCoord( i ,1)
        !    y_tr=obj%FEMDomain%Mesh%NodCoord( i ,2)
        !    if( xmin_m <= x_tr .and. x_tr<= xmax_m )then
        !        if( ymin_m <= y_tr .and. y_tr<= ymax_m )then
        !            ! in 
        !            itr=itr+1
        !            in_out_m(i)=0
        !        else
        !            ! out
        !            in_out_m(i)=1
        !            
        !        endif
        !    else
        !        !out
        !        in_out_m(i)=1
        !        
        !    endif
        !enddo
        !allocate(surfacenod_m(itr,2 ))
        ! only for roots surrounded by soils
        ! detect two ends
        !do i=1,n-1
        !    if(in_out_m(i)==0 .and. in_out_m(i+1)==1 )then
        !        end1_m=i
        !    endif
        !    if(in_out_m(i)==1 .and. in_out_m(i+1)==0 )then
        !        end2_m=i
        !    endif
        !enddo 
!
        !itr=0
        !do i=1,n
        !    x_tr=obj%FEMDomain%Mesh%NodCoord( i ,1)
        !    y_tr=obj%FEMDomain%Mesh%NodCoord( i ,2)
        !    if( xmin_m <= x_tr .and. x_tr<= xmax_m )then
        !        if( ymin_m <= y_tr .and. y_tr<= ymax_m )then
        !            itr=itr+1       
        !            surfacenod_m(itr,1)=x_tr
        !            surfacenod_m(itr,2)=y_tr
        !        else
        !            cycle
        !        endif
        !    else
        !        cycle
        !    endif
        !enddo

        ! add soil surface and root surface
        n=4+sum(in_out)
        allocate(buffer(size(obj%FEMDomain%Mesh%NodCoord,1),2 ) )
        buffer(:,1:2)=obj%FEMDomain%Mesh%NodCoord(:,1:2)
        deallocate( obj%FEMDomain%Mesh%NodCoord )
        allocate( obj%FEMDomain%Mesh%NodCoord(n,2 ) )

        ! get subdomain of root domain which is in soil domain 
        ! end1 to end2
        allocate( surf_nod_buffer ( sum(in_out),2 ) )
        itr=0
        i=end1-1

        do 
            itr=itr+1
            i=i+1
            if(i>size(Mesh%SurfaceLine2D,1) )then
                i=1
            endif

            if(in_out(i)==0 )then
                itr=itr-1
                cycle
            endif

            surf_nod_buffer(itr,1:2)=Mesh%NodCoord(Mesh%SurfaceLine2D(i),1:2)
            if(itr>=size(surf_nod_buffer,1))then
                exit
            endif
        enddo

        ! add soil surface (box-shaped)

        xe1=Mesh%NodCoord(Mesh%SurfaceLine2D(end1),1)
        xe2=Mesh%NodCoord(Mesh%SurfaceLine2D(end2),1)
        ye1=Mesh%NodCoord(Mesh%SurfaceLine2D(end1),2)
        ye2=Mesh%NodCoord(Mesh%SurfaceLine2D(end2),2)

        ! get anti-clockwize surface-line
        non=size(surf_nod_buffer,1)
        if(cross1==1 )then
            ! head is down-ward
            print *, "head is down-ward"
            
            if(xe1 < xe2  )then
                do i=1,non
                    obj%FEMDomain%Mesh%NodCoord( i ,1:2) = surf_nod_buffer(i,1:2)
                enddo
            else
                do i=1,non
                    obj%FEMDomain%Mesh%NodCoord( i ,1:2) = surf_nod_buffer(non - i+1,1:2 )
                enddo
            endif
            obj%FEMDomain%Mesh%NodCoord( non+1 ,1) = xmax
            obj%FEMDomain%Mesh%NodCoord( non+2 ,1) = xmax
            obj%FEMDomain%Mesh%NodCoord( non+3 ,1) = xmin
            obj%FEMDomain%Mesh%NodCoord( non+4 ,1) = xmin

            obj%FEMDomain%Mesh%NodCoord( non+1 ,2) = ymin
            obj%FEMDomain%Mesh%NodCoord( non+2 ,2) = ymax
            obj%FEMDomain%Mesh%NodCoord( non+3 ,2) = ymax
            obj%FEMDomain%Mesh%NodCoord( non+4 ,2) = ymin
        elseif(cross1==2)then
            ! head is right-side
            print *, "head is right-side"
            
            if(ye1 < ye2  )then
                do i=1,non
                    obj%FEMDomain%Mesh%NodCoord( i ,1:2) = surf_nod_buffer(i,1:2)
                enddo
            else
                do i=1,non
                    obj%FEMDomain%Mesh%NodCoord( i ,1:2) = surf_nod_buffer(non - i+1,1:2 )
                enddo
            endif
            obj%FEMDomain%Mesh%NodCoord( non+1 ,1) = xmax
            obj%FEMDomain%Mesh%NodCoord( non+2 ,1) = xmin
            obj%FEMDomain%Mesh%NodCoord( non+3 ,1) = xmin
            obj%FEMDomain%Mesh%NodCoord( non+4 ,1) = xmax
            obj%FEMDomain%Mesh%NodCoord( non+1 ,2) = ymax
            obj%FEMDomain%Mesh%NodCoord( non+2 ,2) = ymax
            obj%FEMDomain%Mesh%NodCoord( non+3 ,2) = ymin
            obj%FEMDomain%Mesh%NodCoord( non+4 ,2) = ymin
        elseif(cross1==3)then
            ! head is upper-side
            print *, "head is upper-side"
            
            if(xe1 > xe2  )then
                do i=1,non
                    obj%FEMDomain%Mesh%NodCoord( i ,1:2) = surf_nod_buffer(i,1:2)
                enddo
            else
                do i=1,non
                    obj%FEMDomain%Mesh%NodCoord( i ,1:2) = surf_nod_buffer(non - i+1,1:2 )
                enddo
            endif
            obj%FEMDomain%Mesh%NodCoord( non+1 ,1) = xmin
            obj%FEMDomain%Mesh%NodCoord( non+2 ,1) = xmin
            obj%FEMDomain%Mesh%NodCoord( non+3 ,1) = xmax
            obj%FEMDomain%Mesh%NodCoord( non+4 ,1) = xmax

            obj%FEMDomain%Mesh%NodCoord( non+1 ,2) = ymax
            obj%FEMDomain%Mesh%NodCoord( non+2 ,2) = ymin
            obj%FEMDomain%Mesh%NodCoord( non+3 ,2) = ymin
            obj%FEMDomain%Mesh%NodCoord( non+4 ,2) = ymax
        elseif(cross1==4)then
            ! head is upper-side
            print *, " head is upper-side"
            
            if(ye1 > ye2  )then
                do i=1,non
                    obj%FEMDomain%Mesh%NodCoord( i ,1:2) = surf_nod_buffer(i,1:2)
                enddo
            else
                do i=1,non
                    obj%FEMDomain%Mesh%NodCoord( i ,1:2) = surf_nod_buffer(non - i+1,1:2 )
                enddo
            endif
            obj%FEMDomain%Mesh%NodCoord( non+1 ,1) = xmin
            obj%FEMDomain%Mesh%NodCoord( non+2 ,1) = xmax
            obj%FEMDomain%Mesh%NodCoord( non+3 ,1) = xmax
            obj%FEMDomain%Mesh%NodCoord( non+4 ,1) = xmin

            obj%FEMDomain%Mesh%NodCoord( non+1 ,2) = ymin
            obj%FEMDomain%Mesh%NodCoord( non+2 ,2) = ymin
            obj%FEMDomain%Mesh%NodCoord( non+3 ,2) = ymax
            obj%FEMDomain%Mesh%NodCoord( non+4 ,2) = ymax
        else
            ! inside
            print *, "none of them"
            print *, cross1
            
            obj%FEMDomain%Mesh%NodCoord( 1:non ,1:2) = surf_nod_buffer(1:non,1:2)
            obj%FEMDomain%Mesh%NodCoord( non+1 ,1) = xmin
            obj%FEMDomain%Mesh%NodCoord( non+2 ,1) = xmax
            obj%FEMDomain%Mesh%NodCoord( non+3 ,1) = xmax
            obj%FEMDomain%Mesh%NodCoord( non+4 ,1) = xmin

            obj%FEMDomain%Mesh%NodCoord( non+1 ,2) = ymin
            obj%FEMDomain%Mesh%NodCoord( non+2 ,2) = ymin
            obj%FEMDomain%Mesh%NodCoord( non+3 ,2) = ymax
            obj%FEMDomain%Mesh%NodCoord( non+4 ,2) = ymax
        endif



        return

!        do i=1,size(Mesh%SurfaceLine2D,1)
!            if(Mesh%NodCoord( Mesh%SurfaceLine2D(i),1)==maxval( Mesh%NodCoord(:,1) ) )then
!                s(1)=i
!            elseif(Mesh%NodCoord( Mesh%SurfaceLine2D(i) ,2)==maxval( Mesh%NodCoord(:,2) ) )then
!                s(2)=i
!            elseif(Mesh%NodCoord( Mesh%SurfaceLine2D(i) ,1)==minval( Mesh%NodCoord(:,1) ) )then
!                s(3)=i
!            elseif(Mesh%NodCoord( Mesh%SurfaceLine2D(i) ,2)==minval( Mesh%NodCoord(:,2) ) )then    
!                s(4)=i
!            else
!                cycle
!            endif
!        enddo
!        
!        do i=1,size(obj%FEMDomain%Mesh%NodCoord,1)
!            if(obj%FEMDomain%Mesh%NodCoord(i,1)==maxval( obj%FEMDomain%Mesh%NodCoord(:,1) ) )then
!                s_m(1)=i
!            elseif(obj%FEMDomain%Mesh%NodCoord(i,2)==maxval( obj%FEMDomain%Mesh%NodCoord(:,2) ) )then
!                s_m(2)=i
!            elseif(obj%FEMDomain%Mesh%NodCoord(i,1)==minval( obj%FEMDomain%Mesh%NodCoord(:,1) ) )then
!                s_m(3)=i
!            elseif(obj%FEMDomain%Mesh%NodCoord(i,2)==minval( obj%FEMDomain%Mesh%NodCoord(:,2) ) )then    
!                s_m(4)=i
!            else
!                cycle
!            endif
!        enddo
!
!        direct=dble(s(4)-s(3))/abs(s(4)-s(3) )+dble(s(3)-s(2))/abs(s(3)-s(2) )&
!            +dble(s(2)-s(1))/abs(s(2)-s(1) )+dble(s(1)-s(4))/abs(s(1)-s(4) )
!        direct_m=dble(s_m(4)-s_m(3))/abs(s_m(4)-s_m(3) )+dble(s_m(3)-s_m(2))/abs(s_m(3)-s_m(2) )&
!            +dble(s_m(2)-s_m(1))/abs(s_m(2)-s_m(1) )+dble(s_m(1)-s_m(4))/abs(s_m(1)-s_m(4) )
!
!
!
!        if(direct * direct_m <= 0.0d0)then
!            print *, "opposite direction"
!            do i=1,size(buffer,1)
!                buffer(i,1:2)=obj%FEMDomain%Mesh%NodCoord( size(buffer,1)-i+1    ,1:2)
!            enddo
!        else
!            print *, "same direction"
!            buffer(:,1:2)=obj%FEMDomain%Mesh%NodCoord(:,1:2)
!        endif
!        



        !scall showarray(obj%FEMDomain%Mesh%NodCoord)
        !call showarray(Mesh%NodCoord)
!        i=end1_m
!        itr=0
!        do 
!            itr=itr+1
!            i=i+1
!            if(i>size(buffer,1) )then
!                i=1
!            endif
!
!            if(in_out_m(i)==0 )then
!                itr=itr-1
!                cycle
!            endif
!
!            obj%FEMDomain%Mesh%NodCoord(itr,1:2)=buffer(i,1:2)
!            if(itr>=sum(in_out_m) )then
!                exit
!            endif
!        enddo
!        itr=itr+1
!        obj%FEMDomain%Mesh%NodCoord(itr,1:2)=Mesh%NodCoord(Mesh%SurfaceLine2D(end1),1:2)
!        i=end1



    endif
    
end subroutine
!##################################################


! #########################################################
subroutine ExportAsLodgingSimProcessing(obj,soil,Name,penalypara,displacement)
    class(PreProcessing_),intent(inout):: obj
    type(PreProcessing_),intent(inout):: soil
    real(real64),optional,intent(in)::penalypara,displacement
    character(*),intent(in) :: Name
    real(real64) :: x_max,x_min, y_max, y_min
    integer(int32) :: number_of_node,fh,i
    integer(int32),allocatable :: top_root_list(:)
    integer(int32),allocatable :: rightside_soil_list(:)
    integer(int32),allocatable :: leftside_soil_list(:)
    integer(int32),allocatable :: bottom_soil_list(:)
    integer(int32) :: total_x_cond, total_y_cond
    ! Export As Lodging Simulator 25 file.

    ! bug exists

    ! give ux=0, uy=0 at toproot.
    y_max=maxval(obj%FEMDomain%Mesh%NodCoord(:,2))
    number_of_node=countif(Array=obj%FEMDomain%Mesh%NodCoord(:,2), Equal=.true., Value=y_max)
    allocate(top_root_list(number_of_node))
    top_root_list(:)=getif(Array=obj%FEMDomain%Mesh%NodCoord(:,2),Value=y_max)
    
    ! give ux=0, uy=0 at right-side of soil.
    x_max=maxval(Soil%FEMDomain%Mesh%NodCoord(:,1))
    number_of_node=countif(Array=Soil%FEMDomain%Mesh%NodCoord(:,1), Equal=.true., Value=x_max)
    allocate(rightside_soil_list(number_of_node))
    rightside_soil_list(:)=0
    rightside_soil_list(:)=getif(Array=Soil%FEMDomain%Mesh%NodCoord(:,1),Value=x_max)
    rightside_soil_list(:)=rightside_soil_list(:)+size(obj%FEMDomain%Mesh%NodCoord,1)

    ! give ux=0, uy=0 at left-side of soil.
    x_min=minval(Soil%FEMDomain%Mesh%NodCoord(:,1))
    number_of_node=countif(Array=Soil%FEMDomain%Mesh%NodCoord(:,1), Equal=.true., Value=x_min)
    allocate(leftside_soil_list(number_of_node))
    leftside_soil_list(:)=0
    leftside_soil_list(:)=getif(Array=Soil%FEMDomain%Mesh%NodCoord(:,1),Value=x_min)
    leftside_soil_list(:)=leftside_soil_list(:)+size(obj%FEMDomain%Mesh%NodCoord,1)


    ! give ux=0, uy=0 at the bottom of soil.
    y_min=minval(Soil%FEMDomain%Mesh%NodCoord(:,2))
    number_of_node=countif(Array=Soil%FEMDomain%Mesh%NodCoord(:,2), Equal=.true., Value=y_min)
    allocate(bottom_soil_list(number_of_node))
    bottom_soil_list(:)=0
    bottom_soil_list(:)=getif(Array=Soil%FEMDomain%Mesh%NodCoord(:,2),Value=y_min)
    bottom_soil_list(:)=bottom_soil_list(:)+size(obj%FEMDomain%Mesh%NodCoord,1)

    total_x_cond = size(top_root_list)+size(rightside_soil_list)+size(leftside_soil_list)+size(bottom_soil_list)
    total_y_cond = size(top_root_list)+size(rightside_soil_list)+size(leftside_soil_list)+size(bottom_soil_list)

    ! export info
    fh=22
    open(fh,file=Name)
    write(fh,*) 2
    write(fh,*) 1, size(obj%FEMDomain%Mesh%NodCoord,1)
    write(fh,*) size(obj%FEMDomain%Mesh%NodCoord,1)+1, &
        size(obj%FEMDomain%Mesh%NodCoord,1)+size(soil%FEMDomain%Mesh%NodCoord,1)
    write(fh,*) size(obj%FEMDomain%Mesh%ElemNod,1)
    write(fh,*) size(soil%FEMDomain%Mesh%ElemNod,1)
    write(fh,*) size(obj%FEMDomain%Mesh%NodCoord,1)+size(soil%FEMDomain%Mesh%NodCoord,1)
    call showArray(Mat=obj%FEMDomain%Mesh%NodCoord,FileHandle=fh)
    call showArray(Mat=Soil%FEMDomain%Mesh%NodCoord,FileHandle=fh)
    write(fh,*) size(obj%FEMDomain%Mesh%ElemNod,1)+size(soil%FEMDomain%Mesh%ElemNod,1),&
        size(soil%FEMDomain%Mesh%ElemNod,2) 
    call showArray(Mat=obj%FEMDomain%Mesh%ElemNod,FileHandle=fh)
    call showArray(Mat=Soil%FEMDomain%Mesh%ElemNod,FileHandle=fh,Add=size(obj%FEMDomain%Mesh%NodCoord,1))
    call showArray(Mat=obj%FEMDomain%Mesh%ElemMat,FileHandle=fh)
    Soil%FEMDomain%Mesh%ElemMat(:)=4
    call showArray(Mat=Soil%FEMDomain%Mesh%ElemMat,FileHandle=fh)
    write(fh,*) 4
    write(fh,*) "6038.93994      0.349999994       0.00000000       1.00000002E+20  0.777999997       0.00000000"    
    write(fh,*) "6038.93994      0.349999994       0.00000000       1.00000002E+20  0.777999997       0.00000000"    
    write(fh,*) "60000.0000      0.349999994       0.00000000       1.00000002E+20  0.777999997       0.00000000"    
    write(fh,*) "60000.0000      0.349999994       0.00000000       1.00000002E+20  0.777999997       0.00000000"    
    write(fh,*) total_x_cond, total_y_cond
    call showArray(Mat=top_root_list      ,FileHandle=fh)
    call showArray(Mat=rightside_soil_list,FileHandle=fh)
    call showArray(Mat=leftside_soil_list ,FileHandle=fh)
    call showArray(Mat=bottom_soil_list   ,FileHandle=fh)
    do i=1,size(top_root_list)
        write(fh,*) 0.0d0
    enddo
    do i=1,size(rightside_soil_list)
        write(fh,*) 0.0d0
    enddo
    do i=1,size(leftside_soil_list)
        write(fh,*) 0.0d0
    enddo
    do i=1,size(bottom_soil_list)
        write(fh,*) 0.0d0
    enddo
    call showArray(Mat=top_root_list      ,FileHandle=fh)
    call showArray(Mat=rightside_soil_list,FileHandle=fh)
    call showArray(Mat=leftside_soil_list ,FileHandle=fh)
    call showArray(Mat=bottom_soil_list   ,FileHandle=fh)
    do i=1,size(top_root_list)
        write(fh,*) -1.0d0
    enddo
    do i=1,size(rightside_soil_list)
        write(fh,*) 0.0d0
    enddo
    do i=1,size(leftside_soil_list)
        write(fh,*) 0.0d0
    enddo
    do i=1,size(bottom_soil_list)
        write(fh,*) 0.0d0
    enddo
    write(fh,*) 0
    call obj%FEMDomain%Mesh%getSurface()
    call Soil%FEMDomain%Mesh%getSurface()
    write(fh,*) size(obj%FEMDomain%Mesh%SurfaceLine2D)+size(Soil%FEMDomain%Mesh%SurfaceLine2D)
    call showArray(Mat= obj%FEMDomain%Mesh%SurfaceLine2D,FileHandle=fh)
    call showArray(Mat=Soil%FEMDomain%Mesh%SurfaceLine2D,FileHandle=fh,Add=size(obj%FEMDomain%Mesh%NodCoord,1))
    write(fh,*) 1, size(obj%FEMDomain%Mesh%SurfaceLine2D)
    write(fh,*) size(obj%FEMDomain%Mesh%SurfaceLine2D)+1, &
        size(obj%FEMDomain%Mesh%SurfaceLine2D)+size(Soil%FEMDomain%Mesh%SurfaceLine2D)
    write(fh,*) "0.1000000000000E-01   0.1000000000000E-01"
    write(fh,*) "1  1"
    write(fh,*) 1, size(obj%FEMDomain%Mesh%SurfaceLine2D)+size(Soil%FEMDomain%Mesh%SurfaceLine2D),1 
    write(fh,*) 1
    write(fh,*) "0.5000000000000E+05   0.5000000000000E+05   0.2402100000000E+01   0.5404000000000E+00"
    write(fh,*) "1  800  1"
    close(fh) 

end subroutine ExportAsLodgingSimProcessing
! #########################################################


end module