GET_DIME_SISM.f90

Go to the documentation of this file.

!    Copyright (C) 2012 The SPEED FOUNDATION
!    Author: Ilario Mazzieri
!
!    This file is part of SPEED.
!
!    SPEED is free software; you can redistribute it and/or modify it
!    under the terms of the GNU Affero General Public License as
!    published by the Free Software Foundation, either version 3 of the
!    License, or (at your option) any later version.
!
!    SPEED 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
!    Affero General Public License for more details.
!
!    You should have received a copy of the GNU Affero General Public License
!    along with SPEED.  If not, see <http://www.gnu.org/licenses/>.
 
! i.e When a node is at interface of two partitions, it is conunted only in one partition.
 
 
      subroutine get_dime_sism(Xipo,Yipo,Zipo, &
                                 X1,Y1,Z1, &       
                                 X2,Y2,Z2, &       
                                 X3,Y3,Z3, &       
                                 nnod,xs,ys,zs, &       
                                 node_sism, nn_loc, mpi_id, ind, loc_n_num)                                    
         
      implicit none
      
      integer*4 :: isn,i,nn_loc,mpi_id,ind
      integer*4 :: nnod,node_sism
      
      integer*4, dimension(nn_loc) :: loc_n_num
 
      real*8 :: xipo,yipo,zipo,x1,y1,z1,x2,y2,z2,x3,y3,z3,x4,y4,z4,ux,uy,uz,vx,vy,vz,tol
      real*8 :: p1x,p1y,p1z,p2x,p2y,p2z
      real*8 :: epsilon, twopi, anglesum, costheta,m1,m2
 
      real*8, dimension(4) :: x,y,z
      real*8, dimension(nn_loc) :: xs,ys,zs
 
      node_sism = 0
 
      tol = 6.28
      epsilon = 1.d-8
      twopi = 6.283185307179586476925287
      
      
      x(1) = x1
      y(1) = y1
      z(1) = z1
 
      x(2) = x2
      y(2) = y2
      z(2) = z2
 
      x(3) = x3
      y(3) = y3
      z(3) = z3
 
      x(4) = x1
      y(4) = y1
      z(4) = z1
      
      do isn = 1,nnod
         anglesum = 0.0d0
         do i = 1,3
            p1x = x(i) - xs(isn)
            p1y = y(i) - ys(isn)
            p1z = z(i) - zs(isn) 
            p2x = x(i+1) - xs(isn) 
            p2y = y(i+1) - ys(isn)
            p2z = z(i+1) - zs(isn)
 
            m1 = dsqrt(p1x*p1x + p1y*p1y + p1z*p1z)
            m2 = dsqrt(p2x*p2x + p2y*p2y + p2z*p2z)
      
            if ((m1*m2).le.epsilon) then
              anglesum = twopi
            else
              costheta = (p1x*p2x + p1y*p2y + p1z*p2z)/(m1*m2)
              anglesum = anglesum + dacos(costheta)
            endif
 
                        
        enddo
 
 
        if (anglesum.ge.tol) then
           node_sism = node_sism + 1
        endif
 
      enddo
 
 
      return
      end subroutine get_dime_sism