MAKE_RAYLEIGH_DAMPING_MATRIX.f90 File Reference

Go to the source code of this file.

Functions/Subroutines
subroutine	make_rayleigh_damping_matrix (nn, ct, ww, dd, rho, 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)
	Compute the Rayleigh damping matrix C = A0 M + A1 K, M = mass matrix, K= stiffness matrix.

Function/Subroutine Documentation

make_rayleigh_damping_matrix()

subroutine make_rayleigh_damping_matrix	(	integer*4	nn,
		real*8, dimension(nn)	ct,
		real*8, dimension(nn)	ww,
		real*8, dimension(nn,nn)	dd,
		real*8, dimension(nn,nn,nn)	rho,
		real*8	a11,
		real*8	a12,
		real*8	a13,
		real*8	a21,
		real*8	a22,
		real*8	a23,
		real*8	a31,
		real*8	a32,
		real*8	a33,
		real*8	b11,
		real*8	b12,
		real*8	b13,
		real*8	b21,
		real*8	b22,
		real*8	b23,
		real*8	b31,
		real*8	b32,
		real*8	b33,
		real*8	gg1,
		real*8	gg2,
		real*8	gg3,
		real*8	dd1,
		real*8	dd2,
		real*8	dd3,
		real*8, dimension(nn,nn,nn)	mc_el
	)

Compute the Rayleigh damping matrix C = A0 M + A1 K, M = mass matrix, K= stiffness matrix.

Author

Ilario Mazzieri

Date

November, 2014

Version

1.0

Parameters

[in]	nn	polynomial degree + 1
[in]	ct	LGL nodes
[in]	ww	LGL weights
[in]	dd	matrix of spectral derivates
[in]	rho	mass density
[in]	A11,...,DD3	coefficients for the bilinear map
[out]	mc_el	Rayleigh damping matrix

Definition at line 32 of file MAKE_RAYLEIGH_DAMPING_MATRIX.f90.

      
 
      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) = rho(i,j,k) &                                                
                                 * det_j * ww(i) * ww(j) * ww(k) 
                                   
                   
 
            enddo 
         enddo 
      enddo 
      
      return