MAKE_DAMPING_MATRIX.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/>.
 
 
 
      subroutine make_damping_matrix(nn,ct,ww,dd,rho,gamma,&                                        
                           A11,A12,A13,A21,A22,A23,&        
                           A31,A32,A33,B11,B12,B13,&        
                           B21,B22,B23,B31,B32,B33,&        
                           GG1,GG2,GG3,DD1,DD2,DD3,&        
                           mc_el,mck_el)
      
 
      implicit none 
      
      integer*4 :: nn 
      integer*4 :: i,j,k 
 
      real*8 :: a11,a12,a13,a21,a22,a23,a31,a32,a33 
      real*8 :: b11,b12,b13,b21,b22,b23,b31,b32,b33 
      real*8 :: gg1,gg2,gg3,dd1,dd2,dd3 
      real*8 :: dxdx,dxdy,dxdz,dydx,dydy,dydz,dzdx,dzdy,dzdz,det_j 
 
      real*8, dimension(nn) :: ct,ww 
 
      real*8, dimension(nn,nn) :: dd 
 
      real*8, dimension(nn,nn,nn) :: rho,gamma
      real*8, dimension(nn,nn,nn) :: mc_el,mck_el 
 
     
      
      do k = 1,nn 
         do j = 1,nn 
            do i = 1,nn 
               dxdx = a11 +b12*ct(k) +b13*ct(j) & 
                    + gg1*ct(j)*ct(k) 
               dydx = a21 +b22*ct(k) +b23*ct(j) & 
                    + gg2*ct(j)*ct(k) 
               dzdx = a31 +b32*ct(k) +b33*ct(j) & 
                    + gg3*ct(j)*ct(k) 
               
               dxdy = a12 +b11*ct(k) +b13*ct(i) & 
                    + gg1*ct(k)*ct(i) 
               dydy = a22 +b21*ct(k) +b23*ct(i) & 
                    + gg2*ct(k)*ct(i) 
               dzdy = a32 +b31*ct(k) +b33*ct(i) & 
                    + gg3*ct(k)*ct(i) 
               
               dxdz = a13 +b11*ct(j) +b12*ct(i) & 
                    + gg1*ct(i)*ct(j) 
               dydz = a23 +b21*ct(j) +b22*ct(i) & 
                    + gg2*ct(i)*ct(j) 
               dzdz = a33 +b31*ct(j) +b32*ct(i) & 
                    + gg3*ct(i)*ct(j) 
               
               det_j = dxdz * (dydx*dzdy - dzdx*dydy) & 
                     - dydz * (dxdx*dzdy - dzdx*dxdy) & 
                     + dzdz * (dxdx*dydy - dydx*dxdy) 
 
 
             
                    mc_el(i,j,k) = 2 * gamma(i,j,k) * rho(i,j,k) &                                                
                                 * det_j * ww(i) * ww(j) * ww(k) 
                                   
                   mck_el(i,j,k) = (gamma(i,j,k)**2) * rho(i,j,k) & 
                                * det_j * ww(i) * ww(j) * ww(k) 
 
 
            enddo 
         enddo 
      enddo 
      
      return 
      
      end subroutine make_damping_matrix