add_facets Subroutine

   subroutine add_facets(self, facet)
   !< Add facets to AABB.
   !<
   !< @note Previously stored facets are lost.
   !<
   !< @note Facets added to AABB are removed to facets list that is also returned.
   class(aabb_object),              intent(inout) :: self              !< AABB.
   type(facet_object), allocatable, intent(inout) :: facet(:)          !< Facets list.
   integer(I4P)                                   :: scratch_unit_add  !< Scratch unit file for added facets.
   integer(I4P)                                   :: scratch_unit_rem  !< Scratch unit file for remaining facets.
   integer(I4P)                                   :: rem_facets_number !< Remaining facets number.
   integer(I4P)                                   :: f                 !< Counter.

   self%facets_number = 0
   rem_facets_number = 0
   if (allocated(self%facet)) deallocate(self%facet)
   open(newunit=scratch_unit_add, status='scratch', access='stream', form='unformatted')
   open(newunit=scratch_unit_rem, status='scratch', access='stream', form='unformatted')
   do f=1, size(facet, dim=1)
      if (self%is_inside(point=facet(f)%vertex_1).and.&
          self%is_inside(point=facet(f)%vertex_2).and.&
          self%is_inside(point=facet(f)%vertex_3)) then
         self%facets_number = self%facets_number + 1
         call facet(f)%save_into_file_binary(file_unit=scratch_unit_add)
      else
         rem_facets_number = rem_facets_number + 1
         call facet(f)%save_into_file_binary(file_unit=scratch_unit_rem)
      endif
   enddo
   if (self%facets_number > 0) then
      allocate(self%facet(1:self%facets_number))
      rewind(unit=scratch_unit_add)
      do f=1, self%facets_number
         call self%facet(f)%load_from_file_binary(file_unit=scratch_unit_add)
         call self%facet(f)%compute_metrix
      enddo
   endif
   close(unit=scratch_unit_add)
   deallocate(facet)
   if (rem_facets_number > 0) then
      allocate(facet(1:rem_facets_number))
      rewind(unit=scratch_unit_rem)
      do f=1, rem_facets_number
         call facet(f)%load_from_file_binary(file_unit=scratch_unit_rem)
         call facet(f)%compute_metrix
      enddo
   endif
   close(unit=scratch_unit_rem)
   endsubroutine add_facets