D:/demona/2.0/manual/doc/2-31set/set.f90

Go to the documentation of this file.

subroutine set

use global


implicit none

! model related pointers

integer,pointer :: IM(:,:,:,:)

! setSet related pointers
integer, pointer :: bcVec(:)
integer, pointer :: maxBC


integer,pointer :: minI,minJ,maxI,maxJ,maxT

integer,pointer :: myMaxI,myMaxJ,myMaxT

integer,pointer :: allBC

integer i,j,t,bc

integer dummy

integer myI,myJ,myT,myBC

integer bcVal

integer,allocatable :: jVecTest(:,:),iVecTest(:),tVecTest(:) ! to use first sweep more efficiently.

integer,allocatable :: bcVecTest(:) ! to detect number of bc used in the model - consider all unk

integer,pointer :: mysetVec(:)

integer,allocatable :: setDofVec(:)


!print*,''
!print*,'-----------------'
!print*,'set started'
!print*,'-----------------'


IM=>myModel%physics(currentPhysics)%domain(currentDomain)%dofMatrix

minI=>myModel%physics(currentPhysics)%domain(currentDomain)%minI
minJ=>myModel%physics(currentPhysics)%domain(currentDomain)%minJ
maxI=>myModel%physics(currentPhysics)%domain(currentDomain)%maxI
maxJ=>myModel%physics(currentPhysics)%domain(currentDomain)%maxJ

!print*,currentPhysics

!print*,minI,minJ,maxI,maxJ

dof=>myModel%physics(currentPhysics)%dof



allocate(setDofVec(dof))

setDofVec=0

!print*,dof

!do j=5,1,-1
!       print*,IM(:,j,2,boundary)
!end do

!stop


!-------------------------------------------------------------------------
!< detect the bc's that are in the set.
!< allocate an array of allBc length (index starts with 0 since default is interior and interior is 1)
!< if there is a bc in the set for any unk, 
allBC=>myModel%physics(currentPhysics)%domain(currentDomain)%allBC
allocate(bcVecTest(allBc))
bcVecTest=0


do j=minJ,maxJ
        
        do i=minI,maxI

                do t=1,size(IM,dim=3)

                        !if (IM(i,j,t,flag).eq.1) then
                        if ((IM(i,j,t,flag).eq.1).and.(IM(i,j,t,absolute).ne.0)) then   ! or one check for bc value and skip directly if zero
                                                                                                                                                        ! this line is added to avoid out-of domain cells in FV

                                !bcVecTest(IM(i,j,t,boundary))=1
                                bcVecTest(IM(i,j,1,boundary))=1  ! currently only first unk is used as reference.

                        end if

                end do


        end do


end do
!< bc's are detected
!-------------------------------------------------------------------------
!< detect the size of bcVec
maxBC=>myModel%physics(currentPhysics)%domain(currentDomain)%set%maxBC
maxBC=0
do i=1,allBc

        if (bcVecTest(i).eq.1) then
        
                maxBc=maxBc+1

        end if

end do
!-------------------------------------------------------------------------
!< allocate bcVec
allocate(myModel%physics(currentPhysics)%domain(currentDomain)%set%bcVec(maxBC))

bcVec=>myModel%physics(currentPhysics)%domain(currentDomain)%set%bcVec
bcVec=0

myModel%physics(currentPhysics)%domain(currentDomain)%set%bcVec=0

dummy=1
!< fill bcVec with the "bc index" like 0,2,5..
do i=1,allBc

        if (bcVecTest(i).eq.1) then
        
                bcVec(dummy)=i
                dummy=dummy+1

        end if

end do
!print*,'print bc types'
!print*,myModel%set%bcVec
!print*,'-----------------'
!stop
!-------------------------------------------------------------------------
!< allocate bc
!allocate(myModel%set%bc(0:maxBC-1))
!allocate(myModel%set%bc(maxBC))
allocate(myModel%physics(currentPhysics)%domain(currentDomain)%set%bc(allBc))
!< cannot be initialized directly because this is a derived type
!< all components are initialized separately.
!print*,'type-1'
!print*,'sweep the domain and calculate maxJ of each myModel%set%bc'
!do myT=1,maxBc
!       bcVal=bcVec(myT)
!       myModel%set%bc(bcVal)%maxJ=0
!end do
!< initialize maxJ of each bc
!do myT=1,maxBc
!       myModel%set%bc(myT)%maxJ=0
!end do
do myT=1,allBc
        myModel%physics(currentPhysics)%domain(currentDomain)%set%bc(myT)%maxJ=0
end do


!allocate(jVecTest(maxJ-minJ+1,maxBc))
allocate(jVecTest(minJ:maxJ,maxBc))
jVecTest=0


!< this loop is performed only for BC's in the set
do j=minJ,maxJ

        
        if (maxval(IM(:,j,:,flag)).eq.0) cycle !fast cycle  useful in smoothers - since they are based on single or dof (if block) number of unknowns


        ! detect maxJ and jVec entries using jVecTest
        do myBC=1,maxBc  ! myT is not 


                do i=minI,maxI




                
                        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
                        dummy=IM(i,j,absolute,boundary) !< we assume same bc but can be defined seperately too
                        ! if different bc's will be defined for unks then unk%bc%j%iVec should be used i.e type3 [that does not exist]
                        ! since bc should be a sub of unk
                        ! so keep third index as 1 and use
                        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

                        !< note that, for domain with holes, holes should be default zeros
                        !< we need to check the absolute level whether its unk index is greater than 1
                        if (maxVal(IM(i,j,:,absolute)).eq.0) cycle !< if not in the domain, cycle

                        if (maxVal(IM(i,j,:,flag)).eq.0) cycle !< if not in the set, cycle


                        bc=bcVec(myBC)

                        myMaxJ=>myModel%physics(currentPhysics)%domain(currentDomain)%set%bc(bc)%maxJ
                                
                        if (dummy.eq.bcVec(myBC)) then

                                !myModel%set%bc(bc)%maxJ=myModel%set%bc(bc)%maxJ+1
                                myMaxJ=myMaxJ+1

                                !print*,j,myBC

                                jVecTest(j,myBC)=1  ! will this change to for full bc

                                exit  ! should exit i loop

                        end if

        
                end do

        
        end do


end do
!< maxJ and jVec entries are found.



!< allocate jVec's - fill the entries from jVecTest
do myBC=1,maxBc 

        bc=bcVec(myBC)

        myMaxJ=>myModel%physics(currentPhysics)%domain(currentDomain)%set%bc(bc)%maxJ

        !print*,myMaxJ
        !stop

        !< allocate jVec
        allocate(myModel%physics(currentPhysics)%domain(currentDomain)%set%bc(bc)%jVec(myMaxJ))
        myModel%physics(currentPhysics)%domain(currentDomain)%set%bc(bc)%jVec(myMaxJ)=0

        !<fill the jVec

        dummy=1

        do myJ=minJ,maxJ

                if (jVecTest(myJ,myBC).eq.1)  then

                        !print*,myJ

                        myModel%physics(currentPhysics)%domain(currentDomain)%set%bc(bc)%jVec(dummy)=myJ

                        dummy=dummy+1

                end if

        end do

        !< allocate j
        allocate(myModel%physics(currentPhysics)%domain(currentDomain)%set%bc(bc)%j(myMaxJ))
        !< myModel%set%bc(t)%j's cannot be initialized directly since its a type
        !< its elements will be initialized

        ! sweep over j's of the set only
        do myJ=1,myMaxJ
                
                j=myModel%physics(currentPhysics)%domain(currentDomain)%set%bc(bc)%jVec(myJ)


                !print*,myModel%set%bc(bc)%jVec
                !stop
                

                myMaxI=>myModel%physics(currentPhysics)%domain(currentDomain)%set%bc(bc)%j(myJ)%maxI
                myMaxI=0

                !allocate(iVecTest(maxI-minI+1))
                allocate(iVecTest(minI:maxI))
                iVecTest=0

                !< sweep i nodes
                do i=minI,maxI



                        if (maxVal(IM(i,j,:,absolute)).eq.0) cycle !< if not in the domain, cycle

                        if (maxVal(IM(i,j,:,flag)).eq.0) cycle !< if not in the set, cycle

                        if (IM(i,j,absolute,boundary).ne.bc) cycle

                        !print*,i,j

                        myMaxI=myMaxI+1

                        iVecTest(i)=1  ! will this change to for full bc

                end do
                !< maxI and iVec entries are found

                !print*,myMaxJ,j
                !stop
                
                allocate(myModel%physics(currentPhysics)%domain(currentDomain)%set%bc(bc)%j(myJ)%i(myMaxI))
                allocate(myModel%physics(currentPhysics)%domain(currentDomain)%set%bc(bc)%j(myJ)%iVec(myMaxI))

                myModel%physics(currentPhysics)%domain(currentDomain)%set%bc(bc)%j(myJ)%iVec=0

                
                !<fill the iVec's
                
                dummy=1

                do myI=minI,maxI

                        if (iVecTest(myI).eq.1)  then

                                myModel%physics(currentPhysics)%domain(currentDomain)%set%bc(bc)%j(myJ)%iVec(dummy)=myI

                                dummy=dummy+1

                        end if

                end do

                deallocate(iVecTest)

        end do

end do


!print*,mymodel%set%bc(3)%j(2)%iVec

!stop




do myBC=1,maxBC
        
        bc=myModel%physics(currentPhysics)%domain(currentDomain)%set%bcVec(myBc)

        myMaxJ=>myModel%physics(currentPhysics)%domain(currentDomain)%set%bc(bc)%maxJ

        do myJ=1,myMaxJ

                j=myModel%physics(currentPhysics)%domain(currentDomain)%set%bc(bc)%jVec(myJ)

                myMaxI=>myModel%physics(currentPhysics)%domain(currentDomain)%set%bc(bc)%j(myJ)%maxI


                
                do myI=1,myMaxI

                        i=myModel%physics(currentPhysics)%domain(currentDomain)%set%bc(bc)%j(myJ)%iVec(myI)

                        !< the two if statements below are not needed since the set is already created for in-set variables
                        !if (maxVal(IM(i,j,:,absolute)).eq.0) cycle !< if not in the domain, cycle

                        !if (maxVal(IM(i,j,:,flag)).eq.0) cycle !< if not in the set, cycle


                        allocate(myModel%physics(currentPhysics)%domain(currentDomain)%set%bc(bc)%j(myJ)%i(myI)%endIndex(dof))
                        myModel%physics(currentPhysics)%domain(currentDomain)%set%bc(bc)%j(myJ)%i(myI)%endIndex=0
                        
                        myMaxT=>myModel%physics(currentPhysics)%domain(currentDomain)%set%bc(bc)%j(myJ)%i(myI)%maxT
                        myMaxT=0

        
                        allocate(tVecTest(dof))
                        tVecTest=0

                        do t=1,dof ! dof is a global variable for localdof

                        
                                if (IM(i,j,t,flag).eq.1) then

                                        myModel%physics(currentPhysics)%domain(currentDomain)%set%bc(bc)%j(myJ)%i(myI)%endIndex(t)=1

                                        myMaxT=myMaxT+1

                                        tVecTest(t)=1

                                        setDofVec(t)=1

                                end if


                        end do

                        allocate(myModel%physics(currentPhysics)%domain(currentDomain)%set%bc(bc)%j(myJ)%i(myI)%tVec(myMaxT))
                        myModel%physics(currentPhysics)%domain(currentDomain)%set%bc(bc)%j(myJ)%i(myI)%tVec=0

                        !<fill the iVec's
                        
                        dummy=1

                        do t=1,dof

                                if (tVecTest(t).eq.1)  then

                                        myModel%physics(currentPhysics)%domain(currentDomain)%set%bc(bc)%j(myJ)%i(myI)%tVec(dummy)=t

                                        dummy=dummy+1

                                end if

                        end do


                        deallocate(tVecTest)

                
                end do

        

        end do


end do


!print*,IM(:,3,1,boundary)
!print*,IM(:,2,1,boundary)
!print*,IM(:,1,1,boundary)

!stop

!print*,mymodel%set%bc(1)%j(1)%iVec
!stop

!print*,mymodel%set%bc(3)%j(2)%i(1)%endindex


setDof=>myModel%physics(currentPhysics)%domain(currentDomain)%setDof

setDof=0

do t=1,dof

        if (setDofVec(t).eq.1) setDof=setDof+1

end do


!if (allocated(myModel%physics(currentPhysics)%domain(currentDomain)%dofVec)) deallocate(myModel%physics(currentPhysics)%domain(currentDomain)%dofVec)

allocate(myModel%physics(currentPhysics)%domain(currentDomain)%dofVec(setDof))
myModel%physics(currentPhysics)%domain(currentDomain)%dofVec=0


dummy=1
do t=1,dof
        
        if (setDofVec(t).eq.1) then

                myModel%physics(currentPhysics)%domain(currentDomain)%dofVec(dummy)=t
                dummy=dummy+1

        end if

end do




IM(:,:,:,relative)=0

dummy=1

myModel%physics(currentPhysics)%domain(currentDomain)%firstIndex=dummy+globalCorrection

select case(relativeModel)


        case(relativeLex)

                minI=>myModel%physics(currentPhysics)%domain(currentDomain)%minI
                minJ=>myModel%physics(currentPhysics)%domain(currentDomain)%minJ
                maxI=>myModel%physics(currentPhysics)%domain(currentDomain)%maxI
                maxJ=>myModel%physics(currentPhysics)%domain(currentDomain)%maxJ


                do j=minJ,maxJ


                        do i=minI,maxI

                                do t=1,dof

                                        if (IM(i,j,t,absolute).eq.0) cycle ! to be sure that the point is in the domain
                
                                        if (IM(i,j,t,flag).eq.0) cycle
                                
                                        IM(i,j,t,relative)=dummy+globalCorrection

                                        dummy=dummy+1

                                end do

                        end do

                end do


        case(relativeBC)

                do myBc=1,maxBc
                        
                        bc=myModel%physics(currentPhysics)%domain(currentDomain)%set%bcVec(myBc)
                        maxJ=>myModel%physics(currentPhysics)%domain(currentDomain)%set%bc(bc)%maxJ  ! maxJ can be replaced now ! could also been replaced in the previous loop

                        !print*,'-------'
                        !print*,myBC
                        
                        do myJ=1,maxJ

                                j=myModel%physics(currentPhysics)%domain(currentDomain)%set%bc(bc)%jVec(myJ)
                                maxI=>myModel%physics(currentPhysics)%domain(currentDomain)%set%bc(bc)%j(myJ)%maxI

                                !print*,maxI

                                do myI=1,maxI

                                        i=myModel%physics(currentPhysics)%domain(currentDomain)%set%bc(bc)%j(myJ)%iVec(myI)
                                        maxT=>myModel%physics(currentPhysics)%domain(currentDomain)%set%bc(bc)%j(myJ)%i(myI)%maxT

                                        do myT=1,maxT
                                                
                                                t=myModel%physics(currentPhysics)%domain(currentDomain)%set%bc(bc)%j(myJ)%i(myI)%tVec(myT)

                                                IM(i,j,t,relative)=dummy+globalCorrection

                                                dummy=dummy+1

                                        end do

                                end do

                        end do

                        !pause


                end do

        case default

                print*,'invalid relative indexing option'
                print*,'demona-set stop!'
                stop


end select


currentSetLength=dummy-1

myModel%physics(currentPhysics)%domain(currentDomain)%lastIndex=currentSetLength


!print*,'-------'
!print*,'size of relative problem'
!print*,currentSetLength
!print*,'-------'
!print*,' '



allocate(myModel%physics(currentPhysics)%domain(currentDomain)%setVec(currentSetLength))

mysetVec=>myModel%physics(currentPhysics)%domain(currentDomain)%setVec

!< fill setVec with absolute indices of set values
do myBc=1,maxBc
                        
        bc=myModel%physics(currentPhysics)%domain(currentDomain)%set%bcVec(myBc)
        maxJ=>myModel%physics(currentPhysics)%domain(currentDomain)%set%bc(bc)%maxJ  ! maxJ can be replaced now ! could also been replaced in the previous loop

        !print*,'-------'
        !print*,myBC
                        
        do myJ=1,maxJ

                j=myModel%physics(currentPhysics)%domain(currentDomain)%set%bc(bc)%jVec(myJ)
                maxI=>myModel%physics(currentPhysics)%domain(currentDomain)%set%bc(bc)%j(myJ)%maxI

                !print*,maxI

                do myI=1,maxI

                        i=myModel%physics(currentPhysics)%domain(currentDomain)%set%bc(bc)%j(myJ)%iVec(myI)
                        maxT=>myModel%physics(currentPhysics)%domain(currentDomain)%set%bc(bc)%j(myJ)%i(myI)%maxT

                        do myT=1,maxT
                                                
                                t=myModel%physics(currentPhysics)%domain(currentDomain)%set%bc(bc)%j(myJ)%i(myI)%tVec(myT)

                                mysetVec(IM(i,j,t,relative)-globalCorrection)=IM(i,j,t,absolute)
                                                                                
                        end do

                end do

        end do

        !pause

end do

globalCorrection=currentSetLength+globalCorrection

end subroutine set

---