initialize Subroutine

   subroutine initialize(self, refinement_levels, facet, bmin, bmax)
   !< Initialize AABB tree.
   class(aabb_tree_object), intent(inout)        :: self              !< AABB tree.
   integer(I4P),            intent(in)           :: refinement_levels !< Total number of refinement levels used.
   type(facet_object),      intent(in), optional :: facet(:)          !< Facets list.
   type(vector_R8P),        intent(in), optional :: bmin              !< Minimum point of AABB.
   type(vector_R8P),        intent(in), optional :: bmax              !< Maximum point of AABB.
   integer(I4P)                                  :: level             !< Counter.
   integer(I4P)                                  :: b, bb, bbb, bbbb  !< Counter.
   integer(I4P)                                  :: parent            !< Parent node index.
   type(aabb_object)                             :: octant(8)         !< AABB octants.
   type(facet_object), allocatable               :: facet_(:)         !< Facets list, local variable.

   call self%destroy
   self%refinement_levels = refinement_levels
   self%nodes_number = nodes_number(refinement_levels=self%refinement_levels)
   allocate(self%node(0:self%nodes_number-1))
   associate(node=>self%node)
      ! inizialize all tree nodes with only the bounding box
      call node(0)%initialize(facet=facet, bmin=bmin, bmax=bmax)
      do level=1, self%refinement_levels                                                             ! loop over refinement levels
         b = first_node(level=level)                                                                 ! first node at level
         do bb=1, nodes_number_at_level(level=level), TREE_RATIO                                     ! loop over nodes at level
            bbb = b + bb - 1                                                                         ! node numeration in tree
            parent = parent_node(node=bbb)                                                           ! parent of the current node
            if (node(parent)%is_allocated()) then                                                    ! create children nodes
               call node(parent)%compute_octants(octant=octant)                                      ! compute parent AABB octants
               do bbbb=0, TREE_RATIO-1                                                               ! loop over children
                  call node(bbb+bbbb)%initialize(bmin=octant(bbbb+1)%bmin, bmax=octant(bbbb+1)%bmax) ! initialize node
               enddo
            endif
         enddo
      enddo

      ! fill all tree nodes with facets
      if (present(facet)) then
         allocate(facet_, source=facet)

         ! add facets to nodes
         do level=self%refinement_levels, 0, -1           ! loop over refinement levels
            b = first_node(level=level)                   ! first node at level
            do bb=1, nodes_number_at_level(level=level)   ! loop over nodes at level
               bbb = b + bb - 1                           ! node numeration in tree
               if (allocated(facet_)) then                ! check if facets list still has facets
                  call node(bbb)%add_facets(facet=facet_) ! add facets to node and prune added facets from list
               endif
            enddo
         enddo

         ! destroy void nodes (except root node)
         do level=self%refinement_levels, 1, -1                        ! loop over refinement levels
            b = first_node(level=level)                                ! first node at level
            do bb=1, nodes_number_at_level(level=level)                ! loop over nodes at level
               bbb = b + bb - 1                                        ! node numeration in tree
               if (.not.node(bbb)%has_facets()) call node(bbb)%destroy ! destroy void node
            enddo
         enddo

         ! update AABB extents
         do level=self%refinement_levels, 1, -1         ! loop over refinement levels
            b = first_node(level=level)                 ! first node at level
            do bb=1, nodes_number_at_level(level=level) ! loop over nodes at level
               bbb = b + bb - 1                         ! node numeration in tree
               call node(bbb)%update_extents            ! update extents
            enddo
         enddo
      endif
   endassociate
   self%is_initialized = .true.
   endsubroutine initialize