module InsectClass use fem use LeafClass implicit none type :: Insect_ type(Leaf_),pointer :: Leaf real(real64) :: x(3) ! position: x, y, z real(real64) :: volume ! m^3 real(real64) :: eatSpeed ! m^3/day contains procedure :: create => createInsect procedure :: set => setInsect procedure :: eat => eatInsect end type contains ! ################################################# subroutine createInsect(obj,x,y,z,volume,eatSpeed) class(Insect_),intent(inout) :: obj real(real64),optional,intent(in) :: x,y,z,volume,eatSpeed obj%x(1) = input(default=0.0d0, option=x) obj%x(2) = input(default=0.0d0, option=y) obj%x(3) = input(default=0.0d0, option=z) obj%volume = input(default=0.0d0, option=volume) obj%eatSpeed = input(default=1.0d0, option=eatSpeed) end subroutine ! ################################################# ! ################################################# subroutine setInsect(obj,leaf) class(Insect_),intent(inout) :: obj type(Leaf_),target,intent(in) :: leaf ! set insect on domains if(associated(obj%leaf) )then nullify(obj%leaf) endif obj%leaf => leaf end subroutine ! ################################################# ! ################################################# subroutine eatInsect(obj, dt,debug) class(Insect_),intent(inout) :: obj real(real64),intent(in) :: dt type(Random_) :: random logical,optional,intent(in) :: debug integer(int32) :: i,j,n,elemid,elem_num,nodeid,itr,nodeid_tr,k,old_elemid integer(int32),allocatable :: eatElemList(:),elemnod(:,:),neiElemList(:) real(real64),allocatable :: x(:) real(real64) :: relem,dv,dist,dist_tr if(.not. associated(obj%leaf) )then print *, "Caution :: InsectClass :: No leaf to eat!!" return endif ! Get a random element elem_num = obj%leaf%femdomain%ne() relem = dble(elem_num-1) * random%random() + 1.0d0 elemid = int(relem) obj%x(:) = obj%leaf%femdomain%mesh%getCenterCoordinate(elemid=elemid) allocate(eatElemList(elem_num)) eatElemList(:) = 0 ! start from elemid itr=0 do itr = itr + 1 ! try to eat untill full dv = obj%leaf%femdomain%getVolume(elemid) if(obj%volume+dv < dt*obj%eatSpeed)then print *,obj%volume,dv endif if(obj%volume+dv > dt*obj%eatSpeed)then exit else eatElemList(elemid) = 1 obj%volume = obj%volume + dv ! move to next element neiElemList = obj%leaf%femdomain%mesh%getNeighboringElement(elemid) old_elemid = elemid do i=1,size(neiElemList) ! search neighbor elements if( eatElemList( neiElemList(i) ) ==0 )then elemid = neiElemList(i) exit endif enddo if(old_elemid == elemid)then ! search over all elements do i=1,size(eatElemList) if(eatElemList(i) ==0)then elemid = eatElemList(i) exit endif enddo endif if(old_elemid == elemid)then ! all elements are eaten exit endif endif enddo do i=size(eatElemList),1,-1 if(eatElemList(i)==1 )then call removeArray(& mat=obj%leaf%femdomain%mesh%elemnod,& remove1stColumn=.true.,& NextOf=i-1 & ) endif enddo call obj%leaf%femdomain%mesh%clean() return ! regacy ! find the nearest element => not implemented now. ! start from a random element ! elem_num = obj%leaf%femdomain%ne() ! relem = dble(elem_num-1) * random%random() + 1.0d0 ! ! elemid = int(relem) ! ! nodeid = obj%leaf%femdomain%mesh%elemnod(elemid,1) ! ! obj%x(:) = obj%leaf%femdomain%mesh%nodcoord(nodeid,:) ! ! allocate(eatElemList(elem_num)) ! eatElemList(:) = 0 ! itr = 0 ! ! do ! itr = itr + 1 ! dv = 0.0d0 ! !nodeid = 0 ! if(minval(eatElemList)==1 )then ! exit ! endif ! k=0 ! ! greedy method ! ! find nearest element ! do i=1,elem_num ! k=k+1 ! nodeid_tr = obj%leaf%femdomain%mesh%elemnod(i,1) ! if(nodeid_tr==nodeid)then ! cycle ! endif ! ! if(k==1)then ! nodeid_tr = obj%leaf%femdomain%mesh%elemnod(i,1) ! ! x(:) = obj%leaf%femdomain%mesh%nodcoord(nodeid_tr,:) ! dist = sqrt(dot_product(x(1:3)-obj%x(1:3),x-obj%x(1:3) ) ) ! elemid = i ! cycle ! endif ! ! if(eatElemList(i)==1 )then ! cycle ! else ! nodeid_tr = obj%leaf%femdomain%mesh%elemnod(i,1) ! x(:) = obj%leaf%femdomain%mesh%nodcoord(nodeid_tr,:) ! dist_tr = sqrt(dot_product(x(1:3)-obj%x(1:3),x-obj%x(1:3) ) ) ! if(dist_tr < dist)then ! dist = dist_tr ! elemid = i ! endif ! endif ! enddo ! ! ! ! try to eat untill full ! dv = obj%leaf%femdomain%getVolume(elemid) ! eatElemList(elemid) = 1 ! ! !print *, itr, dv,nodeid,dist_tr,obj%x(:) ! !stop ! !print *, nodeid,"-",nodeid_tr ! ! !print *, elemid,sum(eatElemList) ! !if(itr==3)then ! ! stop ! !endif ! ! ! if(sum(eatElemList) == size(eatElemList) )then ! print *, "all of the leaf was eaten" ! exit ! endif ! ! if(dv+obj%volume > dt*obj%eatSpeed)then ! eatElemList(elemid) = 0 ! exit ! else ! ! insect eats ! obj%volume = obj%volume + dv ! ! if(present(debug) )then ! if(debug .eqv. .true.)then ! if(mod(itr,100)==0 )then ! print *, obj%volume,"/",dt*obj%eatSpeed ! endif ! endif ! endif ! eatElemList(elemid) = 1 ! ! !print *, elemid,sum(eatElemList) ! ! insect moves ! nodeid = obj%leaf%femdomain%mesh%elemnod(elemid,1) ! obj%x(:) = obj%leaf%femdomain%mesh%nodcoord(nodeid,:) ! endif ! ! ! enddo ! ! ! !print *, itr, obj%volume ! !print *, sum(eatElemList) ! ! ! do i=1,size(eatElemList) ! if(eatElemList(i)==1 )then ! call removeArray(& ! mat=obj%leaf%femdomain%mesh%elemnod,& ! remove1stColumn=.true.,& ! NextOf=i-1 & ! ) ! endif ! enddo ! ! call obj%leaf%femdomain%mesh%clean() end subroutine eatInsect ! ################################################# end module