box_shape.F90

Go to the documentation of this file.

!! Copyright (C) 2021 M. Oliveira, K. Lively, A. Obzhirov, I. Albar
!!
!! This program is free software; you can redistribute it and/or modify
!! it under the terms of the GNU General Public License as published by
!! the Free Software Foundation; either version 2, or (at your option)
!! any later version.
!!
!! This program is distributed in the hope that it will be useful,
!! but WITHOUT ANY WARRANTY; without even the implied warranty of
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!! GNU General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with this program; if not, write to the Free Software
!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
!! 02110-1301, USA.
!!
 
#include "global.h"
 
module box_shape_oct_m
  use box_oct_m
  use basis_vectors_oct_m
  use debug_oct_m
  use global_oct_m
  use math_oct_m
  use messages_oct_m
  use namespace_oct_m
  use profiling_oct_m
 
  implicit none
 
  private
  public ::         &
    box_shape_t,    &
    box_shape_init, &
    box_shape_end
 
  type, abstract, extends(box_t) :: box_shape_t
    private
    real(real64), allocatable,    public :: center(:)
    type(basis_vectors_t), public :: axes
    real(real64), allocatable            :: bounding_box(:,:)
  contains
    procedure :: contains_points => box_shape_contains_points
    procedure :: bounds => box_shape_bounds
    procedure(box_shape_shape_contains_points), deferred :: shape_contains_points
  end type box_shape_t
 
  abstract interface
    function box_shape_shape_contains_points(this, nn, xx) result(contained)
      import :: box_shape_t
      import :: real64
      class(box_shape_t), intent(in) :: this
      integer,            intent(in) :: nn
      real(real64), contiguous,  intent(in) :: xx(:,:)
      logical :: contained(1:nn)
    end function box_shape_shape_contains_points
  end interface
 
contains
 
  !--------------------------------------------------------------
  recursive function box_shape_contains_points(this, nn, xx) result(contained)
    class(box_shape_t), intent(in)  :: this
    integer,            intent(in)  :: nn
    real(real64), contiguous,  intent(in)  :: xx(:,:)
    logical :: contained(1:nn)
 
    ! no push_sub because this function is called very frequently
 
    contained = this%shape_contains_points(nn, this%axes%from_cartesian(nn, xx))
 
  end function box_shape_contains_points
 
  !--------------------------------------------------------------
  subroutine box_shape_init(this, namespace, dim, center, bounding_box_min, bounding_box_max, axes)
    class(box_shape_t), intent(inout) :: this
    type(namespace_t),  intent(in)    :: namespace
    integer,            intent(in)    :: dim
    real(real64),       intent(in)    :: center(dim)
    real(real64),       intent(in)    :: bounding_box_min(dim)
    real(real64),       intent(in)    :: bounding_box_max(dim)
    real(real64), optional,    intent(in)    :: axes(dim, dim)
 
    integer :: idir
    real(real64) :: normalized_axes(dim, dim)
 
    this%dim = dim
 
    safe_allocate(this%center(1:dim))
    safe_allocate(this%bounding_box_l(1:dim))
    this%bounding_box_l = m_zero
 
    this%center(1:dim) = center(1:dim)
 
    if (present(axes)) then
      do idir = 1, dim
        normalized_axes(:, idir) = axes(:, idir)/norm2(axes(:, idir))
      end do
    else
      normalized_axes = diagonal_matrix(dim, m_one)
    end if
    this%axes = basis_vectors_t(namespace, dim, normalized_axes)
 
    if (allocated(this%bounding_box)) then
      safe_deallocate_a(this%bounding_box)
    end if
    safe_allocate(this%bounding_box(1:2, 1:dim))
    this%bounding_box(1,:) = bounding_box_min + this%center
    this%bounding_box(2,:) = bounding_box_max + this%center
 
  end subroutine box_shape_init
 
  !--------------------------------------------------------------
  function box_shape_bounds(this, axes) result(bounds)
    class(box_shape_t),               intent(in)  :: this
    class(basis_vectors_t), optional, intent(in)  :: axes
    real(real64) :: bounds(2, this%dim)
 
    integer :: idir, ib
    real(real64) :: vertex_red(this%dim), vertex_cart(this%dim)
 
    push_sub(box_shape_bounds)
 
    if (this%is_inside_out()) then
      ! The box is unbound.
      bounds(1, :) = -m_huge
      bounds(2, :) =  m_huge
    else
      ! The bounds must be located at the axis-aligned bounding box vertices, so
      ! we loop over the vertices, convert them to reduced coordinates, and find
      ! the maximum and minimum coordinates. This algorithm might not be
      ! optimal, but it will do for now.
      bounds(1, :) =  m_huge
      bounds(2, :) = -m_huge
      do idir = 1, this%dim
        do ib = 1, 2
          ! Vertex in box reduced coordinates
          vertex_red = m_zero
          vertex_red(idir) = this%bounding_box(ib, idir)
 
          ! Transform vertex to Cartesian coordinates
          vertex_cart = this%axes%to_cartesian(vertex_red)
 
          ! Convert to reduced coordinates of the input axes
          if (present(axes)) then
            vertex_red = axes%from_cartesian(vertex_cart)
          else
            vertex_red = vertex_cart
          end if
 
          where (vertex_red < bounds(1,:))
            bounds(1, :) = vertex_red
          elsewhere (vertex_red > bounds(2,:))
            bounds(2, :) = vertex_red
          end where
        end do
      end do
    end if
 
    pop_sub(box_shape_bounds)
  end function box_shape_bounds
 
  !--------------------------------------------------------------
  subroutine box_shape_end(this)
    class(box_shape_t), intent(inout) :: this
 
    safe_deallocate_a(this%center)
    safe_deallocate_a(this%bounding_box)
    safe_deallocate_a(this%bounding_box_l)
 
  end subroutine box_shape_end
 
end module box_shape_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: