617 lines
19 KiB
Fortran
617 lines
19 KiB
Fortran
!
|
|
! Copyright (C) 1996-2016 The SIESTA group
|
|
! This file is distributed under the terms of the
|
|
! GNU General Public License: see COPYING in the top directory
|
|
! or http://www.gnu.org/copyleft/gpl.txt.
|
|
! See Docs/Contributors.txt for a list of contributors.
|
|
!
|
|
! This code segment has been fully created by:
|
|
! Nick Papior Andersen, 2013, nickpapior@gmail.com
|
|
! Please conctact the author, prior to re-using this code.
|
|
!
|
|
|
|
! This module contains the geometrical object:
|
|
! SQUARE
|
|
|
|
module m_geom_square
|
|
|
|
use m_geom_aux
|
|
use m_geom_plane, only : geo_plane_delta
|
|
use m_geom_plane, only : geo_plane_gauss
|
|
use m_geom_plane, only : geo_plane_exp
|
|
use m_geom_plane, only : voxel_in_plane
|
|
use m_geom_plane, only : voxel_val_plane
|
|
use m_geom_plane, only : correct_plane
|
|
|
|
implicit none
|
|
|
|
private
|
|
|
|
! Types
|
|
public :: geo_square_delta
|
|
public :: geo_square_gauss
|
|
public :: geo_square_exp
|
|
! Interfaces
|
|
public :: fgeo_read_square
|
|
public :: correct_square
|
|
public :: voxel_in_square
|
|
public :: voxel_val_square
|
|
|
|
! Explicit routines
|
|
public :: fgeo_read_square_delta
|
|
public :: fgeo_read_square_gauss
|
|
public :: fgeo_read_square_exp
|
|
public :: correct_square_delta
|
|
public :: correct_square_gauss
|
|
public :: correct_square_exp
|
|
public :: voxel_in_square_delta
|
|
public :: voxel_in_square_gauss
|
|
public :: voxel_in_square_exp
|
|
public :: voxel_val_square_delta
|
|
public :: voxel_val_square_gauss
|
|
public :: voxel_val_square_exp
|
|
|
|
type geo_square_delta
|
|
! Describing a bounded plane which can be described by the infinite plane
|
|
! One can retrieve the normal vector of the plane spanned by: v1 and v2 by
|
|
! n = (v1 - o) \cross (v2 - o)
|
|
! 'o' describes the origo of the
|
|
! plane spanned by the vectors, v1 and v2.
|
|
! With n we describe the infinite plane as in the previous case
|
|
! and subsequently bound it by:
|
|
! --------------------- --------
|
|
! | v1_1 | v2_1 | n_1 | | x1_1 |
|
|
! | v1_2 | v2_2 | n_2 | x = | x1_2 |
|
|
! | v1_3 | v2_3 | n_3 | | x1_3 |
|
|
! --------------------- --------
|
|
! and x = (a,b,c)^T.
|
|
! Now to be in the bounded plane, a and b must both satisfy:
|
|
! \{a,b\} \in [0;1], we check whether c is small enough via
|
|
! the inifinite plane method.
|
|
sequence
|
|
! The coordinate from which the plane is spanned
|
|
! by the two vectors v1 and v2 is assigned to p%c
|
|
! The vectors which span the plane
|
|
real(dp) :: v1(3), v2(3)
|
|
! Infinite plane
|
|
type(geo_plane_delta) :: p
|
|
end type geo_square_delta
|
|
|
|
type geo_square_gauss
|
|
sequence
|
|
real(dp) :: v1(3), v2(3)
|
|
type(geo_plane_gauss) :: p
|
|
end type geo_square_gauss
|
|
|
|
type geo_square_exp
|
|
sequence
|
|
real(dp) :: v1(3), v2(3)
|
|
type(geo_plane_exp) :: p
|
|
end type geo_square_exp
|
|
|
|
interface fgeo_read_square
|
|
module procedure fgeo_read_square_delta
|
|
module procedure fgeo_read_square_gauss
|
|
module procedure fgeo_read_square_exp
|
|
end interface
|
|
|
|
interface correct_square
|
|
module procedure correct_square_delta
|
|
module procedure correct_square_gauss
|
|
module procedure correct_square_exp
|
|
end interface
|
|
|
|
interface voxel_in_square
|
|
module procedure voxel_in_square_delta
|
|
module procedure voxel_in_square_gauss
|
|
module procedure voxel_in_square_exp
|
|
end interface
|
|
|
|
interface voxel_val_square
|
|
module procedure voxel_val_square_delta
|
|
module procedure voxel_val_square_gauss
|
|
module procedure voxel_val_square_exp
|
|
end interface
|
|
|
|
contains
|
|
|
|
! This subroutine returns .true. in has if the box spanned in
|
|
! lr,lr+dx,lr+dx+dy,..., has any points inside the plane.
|
|
! - lr is the "lower-left" of the box.
|
|
! - d is the extension of the box
|
|
function voxel_in_square_delta(plane,ll,d) result(has)
|
|
type(geo_square_delta), intent(in) :: plane
|
|
real(dp), intent(in) :: ll(3), d(3)
|
|
logical :: has
|
|
|
|
! We retrieve whether the box passes through the
|
|
! infinite plane spanned by the local plane
|
|
has = voxel_in_plane(plane%p,ll,d)
|
|
|
|
! If it does not intersect, we can safely return
|
|
if ( .not. has ) return
|
|
|
|
has = P_voxel_square(ll,d,plane%v1,plane%v2,plane%p%n,plane%p%c)
|
|
|
|
end function voxel_in_square_delta
|
|
|
|
function voxel_in_square_gauss(plane,ll,d) result(has)
|
|
type(geo_square_gauss), intent(in) :: plane
|
|
real(dp), intent(in) :: ll(3), d(3)
|
|
logical :: has
|
|
|
|
has = voxel_in_plane(plane%p,ll,d)
|
|
|
|
if ( .not. has ) return
|
|
|
|
has = P_voxel_square(ll,d,plane%v1,plane%v2,plane%p%n,plane%p%c)
|
|
|
|
end function voxel_in_square_gauss
|
|
|
|
function voxel_in_square_exp(plane,ll,d) result(has)
|
|
type(geo_square_exp), intent(in) :: plane
|
|
real(dp), intent(in) :: ll(3), d(3)
|
|
logical :: has
|
|
|
|
has = voxel_in_plane(plane%p,ll,d)
|
|
|
|
if ( .not. has ) return
|
|
|
|
has = P_voxel_square(ll,d,plane%v1,plane%v2,plane%p%n,plane%p%c)
|
|
|
|
end function voxel_in_square_exp
|
|
|
|
|
|
function P_voxel_square(ll,d,v1,v2,n,o) result(has)
|
|
real(dp), intent(in) :: ll(3), d(3), v1(3), v2(3), n(3), o(3)
|
|
logical :: has
|
|
integer :: ipiv(4)
|
|
real(dp) :: sys(3,4)
|
|
|
|
! We need to check whether the bounds of the voxel
|
|
! lies within the plane...
|
|
! Thus we solve the linear equation:
|
|
! Ax=B
|
|
! where A = [v1,v2,n] and B = ll+.5*d
|
|
sys(:,1) = v1
|
|
sys(:,2) = v2
|
|
! We do need the normal vector, otherwise
|
|
! we would not be able to retrieve the
|
|
! solution to the system. (i.e. we do not
|
|
! know if the center cuts the plane)
|
|
sys(:,3) = n
|
|
! Move the point down to origo
|
|
sys(:,4) = ll + .5_dp * d - o
|
|
call dgesv(3,1,sys(1,1),3,ipiv,sys(1,4),3,ipiv(4))
|
|
if ( ipiv(4) /= 0 ) then
|
|
write(*,*) ' LAPACK-error: Could not solve the linear &
|
|
&equation for the voxel'
|
|
! This should never happen as any point in space
|
|
! should be reacheable by the vectors
|
|
! (we however, need that v1 and v2 are linearly
|
|
! independent)
|
|
end if
|
|
|
|
! We already know that it lies in the plane
|
|
! so now we need to ensure that it lies within
|
|
! the two spanning vectors plane...
|
|
has= 0._dp <= sys(1,4) .and. sys(1,4) <= 1.0_dp .and. &
|
|
0._dp <= sys(2,4) .and. sys(2,4) <= 1.0_dp
|
|
|
|
end function P_voxel_square
|
|
|
|
function voxel_val_square_delta(plane,ll,d) result(vol)
|
|
type(geo_square_delta), intent(in) :: plane
|
|
real(dp), intent(in) :: ll(3), d(3)
|
|
real(dp) :: vol
|
|
vol = 1._dp
|
|
end function voxel_val_square_delta
|
|
|
|
! This subroutine returns the value of the charge
|
|
! that is attributed from the gaussian function about the
|
|
! finite plane
|
|
! - ll is the "lower-left" of the box.
|
|
! - d is the extension of the box
|
|
function voxel_val_square_gauss(plane,ll,d) result(vol)
|
|
type(geo_square_gauss), intent(in) :: plane
|
|
real(dp), intent(in) :: ll(3), d(3)
|
|
real(dp) :: vol
|
|
|
|
! It has already been checked that it exists
|
|
vol = voxel_val_plane(plane%p,ll,d)
|
|
|
|
end function voxel_val_square_gauss
|
|
|
|
function voxel_val_square_exp(plane,ll,d) result(vol)
|
|
type(geo_square_exp), intent(in) :: plane
|
|
real(dp), intent(in) :: ll(3), d(3)
|
|
real(dp) :: vol
|
|
|
|
vol = voxel_val_plane(plane%p,ll,d)
|
|
|
|
end function voxel_val_square_exp
|
|
|
|
|
|
! Subroutine for initializing the planes by transforming them
|
|
! to the corresponding cell
|
|
subroutine correct_square_delta(cell,plane)
|
|
real(dp), intent(in) :: cell(3,3)
|
|
type(geo_square_delta), intent(inout) :: plane
|
|
real(dp), parameter :: EPS = 1e-5
|
|
real(dp) :: fac
|
|
logical :: inde
|
|
integer :: i
|
|
|
|
! The simple test is that if any-one coordinate
|
|
! is zero, and the other is not, then
|
|
! they are linearly independent
|
|
inde = .false.
|
|
do i = 1 , 3
|
|
if ( abs(plane%v1(i)) < EPS .and. &
|
|
abs(plane%v2(i)) > EPS ) then
|
|
inde = .true.
|
|
else if ( abs(plane%v2(i)) < EPS .and. &
|
|
abs(plane%v1(i)) > EPS ) then
|
|
inde = .true.
|
|
else if ( abs(plane%v1(i)) > EPS .and. &
|
|
abs(plane%v2(i)) > EPS ) then
|
|
! We create the factor here (as we must not divide by 0)
|
|
fac = plane%v1(i) / plane%v2(i)
|
|
end if
|
|
end do
|
|
|
|
if ( .not. inde ) then
|
|
inde = .not. ( &
|
|
count( abs(fac * plane%v2 - plane%v1) <= EPS ) == 3 )
|
|
end if
|
|
|
|
|
|
! Kill if the vectors are linearly dependent...
|
|
if ( .not. inde ) then
|
|
call die('The geometric finite plane has linearly &
|
|
&dependent spanning vectors. This is not allowed.')
|
|
end if
|
|
|
|
! Correct the infinite plane that the thing lies in
|
|
call correct_plane(cell,plane%p)
|
|
|
|
end subroutine correct_square_delta
|
|
|
|
! Subroutine for initializing the planes by transforming them
|
|
! to the corresponding cell
|
|
subroutine correct_square_gauss(cell,plane)
|
|
real(dp), intent(in) :: cell(3,3)
|
|
type(geo_square_gauss), intent(inout) :: plane
|
|
real(dp), parameter :: EPS = 1e-5
|
|
type(geo_square_delta) :: t_plane
|
|
type(geo_plane_delta) :: ti_plane
|
|
|
|
t_plane%v1 = plane%v1
|
|
t_plane%v2 = plane%v2
|
|
ti_plane%n = plane%p%n
|
|
ti_plane%c = plane%p%c
|
|
t_plane%p = ti_plane
|
|
call correct_square_delta(cell,t_plane)
|
|
plane%p%n = t_plane%p%n
|
|
plane%p%d_A = t_plane%p%d + plane%p%cut_off
|
|
plane%p%d_B = t_plane%p%d - plane%p%cut_off
|
|
|
|
end subroutine correct_square_gauss
|
|
|
|
! Subroutine for initializing the planes by transforming them
|
|
! to the corresponding cell
|
|
subroutine correct_square_exp(cell,plane)
|
|
real(dp), intent(in) :: cell(3,3)
|
|
type(geo_square_exp), intent(inout) :: plane
|
|
real(dp), parameter :: EPS = 1e-5
|
|
type(geo_square_delta) :: t_plane
|
|
type(geo_plane_delta) :: ti_plane
|
|
|
|
t_plane%v1 = plane%v1
|
|
t_plane%v2 = plane%v2
|
|
ti_plane%n = plane%p%n
|
|
ti_plane%c = plane%p%c
|
|
t_plane%p = ti_plane
|
|
call correct_square_delta(cell,t_plane)
|
|
plane%p%n = t_plane%p%n
|
|
plane%p%d_A = t_plane%p%d + plane%p%cut_off
|
|
plane%p%d_B = t_plane%p%d - plane%p%cut_off
|
|
|
|
end subroutine correct_square_exp
|
|
|
|
|
|
! Reading in information from block about the square-delta object
|
|
subroutine fgeo_read_square_delta(bName, ngeom,geom,params,par_unit)
|
|
|
|
use intrinsic_missing, only : VNORM
|
|
use fdf
|
|
|
|
! ********************
|
|
! * INPUT variables *
|
|
! ********************
|
|
character(len=*), intent(in) :: bName
|
|
integer, intent(in) :: ngeom
|
|
character(len=*), intent(in), optional :: par_unit
|
|
|
|
! ********************
|
|
! * OUTPUT variables *
|
|
! ********************
|
|
type(geo_square_delta), intent(out) :: geom(ngeom)
|
|
real(dp), intent(out) :: params(ngeom)
|
|
|
|
! ********************
|
|
! * LOCAL variables *
|
|
! ********************
|
|
type(block_fdf) :: bfdf
|
|
type(parsed_line), pointer :: pline
|
|
integer :: t, ip
|
|
real(dp) :: v
|
|
|
|
! if none found, simply return
|
|
if ( ngeom <= 0 ) return
|
|
|
|
! Count the number of planes
|
|
call fgeo_count(bName,GEOM_SQUARE_DELTA,t)
|
|
|
|
if ( t /= ngeom ) &
|
|
call die('Could not find any delta squares')
|
|
|
|
! Do the processing
|
|
if ( .not. fdf_block(bName,bfdf) ) &
|
|
call die('Could not find the block again...?')
|
|
|
|
ip = 0
|
|
count_geom: do while ( ip < ngeom )
|
|
|
|
! First loop until we have the geometry
|
|
call fgeo_next(bfdf,t,v,unit=par_unit)
|
|
|
|
! If it is not the correct geometry, cycle...
|
|
if ( t /= GEOM_SQUARE_DELTA ) cycle
|
|
|
|
! Counter for geometry
|
|
ip = ip + 1
|
|
|
|
! The parameter for the geometry
|
|
params(ip) = v
|
|
|
|
! Step delta
|
|
if ( .not. fdf_bline(bfdf,pline) ) &
|
|
call die('Could not step the delta mark')
|
|
|
|
! A plane is defined by the lower-left corner
|
|
if ( .not. fdf_bline(bfdf,pline) ) &
|
|
call die('Could not read the lower-left corner &
|
|
&of a plane geometry object')
|
|
call fgeo_read_vals(pline,geom(ip)%p%c,units=.true.)
|
|
|
|
! A plane is defined by two spanning vectors
|
|
if ( .not. fdf_bline(bfdf,pline) ) &
|
|
call die('Could not read the first vector &
|
|
&of a plane geometry object')
|
|
call fgeo_read_vals(pline,geom(ip)%v1,units=.true.)
|
|
|
|
if ( .not. fdf_bline(bfdf,pline) ) &
|
|
call die('Could not read the second vector &
|
|
&of a plane geometry object')
|
|
call fgeo_read_vals(pline,geom(ip)%v2,units=.true.)
|
|
|
|
! We need to calculate the equivalent
|
|
! infinite plane for this finite plane
|
|
geom(ip)%p%n(1) = &
|
|
geom(ip)%v1(2)*geom(ip)%v2(3) - &
|
|
geom(ip)%v1(3)*geom(ip)%v2(2)
|
|
geom(ip)%p%n(2) = &
|
|
geom(ip)%v1(3)*geom(ip)%v2(1) - &
|
|
geom(ip)%v1(1)*geom(ip)%v2(3)
|
|
geom(ip)%p%n(3) = &
|
|
geom(ip)%v1(1)*geom(ip)%v2(2) - &
|
|
geom(ip)%v1(2)*geom(ip)%v2(1)
|
|
! Normalize the normal vector
|
|
geom(ip)%p%n = geom(ip)%p%n / VNORM(geom(ip)%p%n)
|
|
|
|
end do count_geom
|
|
|
|
end subroutine fgeo_read_square_delta
|
|
|
|
|
|
subroutine fgeo_read_square_gauss(bName, ngeom,geom,params,par_unit)
|
|
|
|
use intrinsic_missing, only : VNORM
|
|
use fdf
|
|
|
|
! ********************
|
|
! * INPUT variables *
|
|
! ********************
|
|
character(len=*), intent(in) :: bName
|
|
integer, intent(in) :: ngeom
|
|
character(len=*), intent(in), optional :: par_unit
|
|
|
|
! ********************
|
|
! * OUTPUT variables *
|
|
! ********************
|
|
type(geo_square_gauss), intent(out) :: geom(ngeom)
|
|
real(dp), intent(out) :: params(ngeom)
|
|
|
|
! ********************
|
|
! * LOCAL variables *
|
|
! ********************
|
|
type(block_fdf) :: bfdf
|
|
type(parsed_line), pointer :: pline
|
|
integer :: t, ip
|
|
real(dp) :: v
|
|
|
|
! if none found, simply return
|
|
if ( ngeom <= 0 ) return
|
|
|
|
! Count the number of planes
|
|
call fgeo_count(bName,GEOM_SQUARE_GAUSS,t)
|
|
|
|
if ( t /= ngeom ) &
|
|
call die('Could not find any Gaussian squares')
|
|
|
|
! Do the processing
|
|
if ( .not. fdf_block(bName,bfdf) ) &
|
|
call die('Could not find the block again...?')
|
|
|
|
ip = 0
|
|
count_geom: do while ( ip < ngeom )
|
|
|
|
! First loop until we have the geometry
|
|
call fgeo_next(bfdf,t,v,unit=par_unit)
|
|
|
|
! If it is not the correct geometry, cycle...
|
|
if ( t /= GEOM_SQUARE_GAUSS ) cycle
|
|
|
|
! Counter for geometry
|
|
ip = ip + 1
|
|
|
|
! The parameter for the geometry
|
|
params(ip) = v
|
|
|
|
! Step to gauss
|
|
if ( .not. fdf_bline(bfdf,pline) ) &
|
|
call die('Could not step the gauss mark')
|
|
|
|
! We need to read in the \sigma and cut_off
|
|
call fgeo_read_vals(pline,geom(ip)%p%c(1:2),units=.true.,unit_offset=1)
|
|
v = geom(ip)%p%c(1)
|
|
geom(ip)%p%inv_2var = 1._dp/(2._dp*v**2)
|
|
geom(ip)%p%cut_off = geom(ip)%p%c(2)
|
|
|
|
! A plane is defined by the lower-left corner
|
|
if ( .not. fdf_bline(bfdf,pline) ) &
|
|
call die('Could not read the lower-left corner &
|
|
&of a plane geometry object')
|
|
call fgeo_read_vals(pline,geom(ip)%p%c,units=.true.)
|
|
|
|
! A plane is defined by two spanning vectors
|
|
if ( .not. fdf_bline(bfdf,pline) ) &
|
|
call die('Could not read the first vector &
|
|
&of a plane geometry object')
|
|
call fgeo_read_vals(pline,geom(ip)%v1,units=.true.)
|
|
|
|
if ( .not. fdf_bline(bfdf,pline) ) &
|
|
call die('Could not read the second vector &
|
|
&of a plane geometry object')
|
|
call fgeo_read_vals(pline,geom(ip)%v2,units=.true.)
|
|
|
|
! We need to calculate the equivalent
|
|
! infinite plane for this finite plane
|
|
geom(ip)%p%n(1) = &
|
|
geom(ip)%v1(2)*geom(ip)%v2(3) - &
|
|
geom(ip)%v1(3)*geom(ip)%v2(2)
|
|
geom(ip)%p%n(2) = &
|
|
geom(ip)%v1(3)*geom(ip)%v2(1) - &
|
|
geom(ip)%v1(1)*geom(ip)%v2(3)
|
|
geom(ip)%p%n(3) = &
|
|
geom(ip)%v1(1)*geom(ip)%v2(2) - &
|
|
geom(ip)%v1(2)*geom(ip)%v2(1)
|
|
! Normalize the normal vector
|
|
geom(ip)%p%n = geom(ip)%p%n / VNORM(geom(ip)%p%n)
|
|
|
|
end do count_geom
|
|
|
|
end subroutine fgeo_read_square_gauss
|
|
|
|
subroutine fgeo_read_square_exp(bName, ngeom,geom,params,par_unit)
|
|
|
|
use intrinsic_missing, only : VNORM
|
|
use fdf
|
|
|
|
! ********************
|
|
! * INPUT variables *
|
|
! ********************
|
|
character(len=*), intent(in) :: bName
|
|
integer, intent(in) :: ngeom
|
|
character(len=*), intent(in), optional :: par_unit
|
|
|
|
! ********************
|
|
! * OUTPUT variables *
|
|
! ********************
|
|
type(geo_square_exp), intent(out) :: geom(ngeom)
|
|
real(dp), intent(out) :: params(ngeom)
|
|
|
|
! ********************
|
|
! * LOCAL variables *
|
|
! ********************
|
|
type(block_fdf) :: bfdf
|
|
type(parsed_line), pointer :: pline
|
|
integer :: t, ip
|
|
real(dp) :: v
|
|
|
|
! if none found, simply return
|
|
if ( ngeom <= 0 ) return
|
|
|
|
! Count the number of planes
|
|
call fgeo_count(bName,GEOM_SQUARE_EXP,t)
|
|
|
|
if ( t /= ngeom ) &
|
|
call die('Could not find any Gaussian squares')
|
|
|
|
! Do the processing
|
|
if ( .not. fdf_block(bName,bfdf) ) &
|
|
call die('Could not find the block again...?')
|
|
|
|
ip = 0
|
|
count_geom: do while ( ip < ngeom )
|
|
|
|
! First loop until we have the geometry
|
|
call fgeo_next(bfdf,t,v,unit=par_unit)
|
|
|
|
! If it is not the correct geometry, cycle...
|
|
if ( t /= GEOM_SQUARE_EXP ) cycle
|
|
|
|
! Counter for geometry
|
|
ip = ip + 1
|
|
|
|
! The parameter for the geometry
|
|
params(ip) = v
|
|
|
|
! Step to gauss
|
|
if ( .not. fdf_bline(bfdf,pline) ) &
|
|
call die('Could not step the gauss mark')
|
|
|
|
! We need to read in the half-length and cut_off
|
|
call fgeo_read_vals(pline,geom(ip)%p%c(1:2),units=.true.,unit_offset=1)
|
|
v = geom(ip)%p%c(1)
|
|
geom(ip)%p%h = log(2._dp) / v
|
|
geom(ip)%p%cut_off = geom(ip)%p%c(2)
|
|
|
|
! A plane is defined by the lower-left corner
|
|
if ( .not. fdf_bline(bfdf,pline) ) &
|
|
call die('Could not read the lower-left corner &
|
|
&of a plane geometry object')
|
|
call fgeo_read_vals(pline,geom(ip)%p%c,units=.true.)
|
|
|
|
! A plane is defined by two spanning vectors
|
|
if ( .not. fdf_bline(bfdf,pline) ) &
|
|
call die('Could not read the first vector &
|
|
&of a plane geometry object')
|
|
call fgeo_read_vals(pline,geom(ip)%v1,units=.true.)
|
|
|
|
if ( .not. fdf_bline(bfdf,pline) ) &
|
|
call die('Could not read the second vector &
|
|
&of a plane geometry object')
|
|
call fgeo_read_vals(pline,geom(ip)%v2,units=.true.)
|
|
|
|
! We need to calculate the equivalent
|
|
! infinite plane for this finite plane
|
|
geom(ip)%p%n(1) = &
|
|
geom(ip)%v1(2)*geom(ip)%v2(3) - &
|
|
geom(ip)%v1(3)*geom(ip)%v2(2)
|
|
geom(ip)%p%n(2) = &
|
|
geom(ip)%v1(3)*geom(ip)%v2(1) - &
|
|
geom(ip)%v1(1)*geom(ip)%v2(3)
|
|
geom(ip)%p%n(3) = &
|
|
geom(ip)%v1(1)*geom(ip)%v2(2) - &
|
|
geom(ip)%v1(2)*geom(ip)%v2(1)
|
|
! Normalize the normal vector
|
|
geom(ip)%p%n = geom(ip)%p%n / VNORM(geom(ip)%p%n)
|
|
|
|
end do count_geom
|
|
|
|
end subroutine fgeo_read_square_exp
|
|
|
|
end module m_geom_square
|