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) end if 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 end if ! 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 end if 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 end if end do if (old_elemid == elemid) then ! search over all elements do i = 1, size(eatElemList) if (eatElemList(i) == 0) then elemid = eatElemList(i) exit end if end do end if if (old_elemid == elemid) then ! all elements are eaten exit end if end if end do do i = size(eatElemList), 1, -1 if (eatElemList(i) == 1) then call removeArray( & mat=obj%leaf%femdomain%mesh%elemnod, & remove1stColumn=.true., & NextOf=i - 1 & ) end if end do 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