SPEED
MAKE_DAMPING_MATRIX_NLE.f90
Go to the documentation of this file.
1! Copyright (C) 2012 The SPEED FOUNDATION
2! Author: Ilario Mazzieri
3!
4! This file is part of SPEED.
5!
6! SPEED is free software; you can redistribute it and/or modify it
7! under the terms of the GNU Affero General Public License as
8! published by the Free Software Foundation, either version 3 of the
9! License, or (at your option) any later version.
10!
11! SPEED is distributed in the hope that it will be useful, but
12! WITHOUT ANY WARRANTY; without even the implied warranty of
13! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14! Affero General Public License for more details.
15!
16! You should have received a copy of the GNU Affero General Public License
17! along with SPEED. If not, see <http://www.gnu.org/licenses/>.
18
72
73 subroutine make_damping_matrix_nle(nn,ct,ww,dd,rho,gamma0,&
74 alfa11,alfa12,alfa13,alfa21,alfa22,alfa23,&
75 alfa31,alfa32,alfa33,beta11,beta12,beta13,&
76 beta21,beta22,beta23,beta31,beta32,beta33,&
77 gamma1,gamma2,gamma3,delta1,delta2,delta3,&
78 mc_el,mck_el,&
79 cs_nnz,cs,ielem,&
80 vcase,R_el,fpeak, &
81 func_type,func_indx,func_data,nfdata,nf,t_stress,tag_func,yon,gamma_new,&
82 tcase, nnod_loc, vs_tria)
83
84
85 implicit none
86
87 integer*4 :: nn
88 integer*4 :: i,j,k
89 integer*4 :: nnt,tagmat
90 integer*4 :: cs_nnz
91 integer*4 :: vcase,nfdata
92 integer*4 :: nf,fn,ielem,tcase,nnod_loc,in,is
93
94 integer*4, dimension(nf) :: func_type
95 integer*4, dimension(nf +1) :: func_indx
96 integer*4, dimension(0:cs_nnz) :: cs
97 integer*4, dimension(nf) :: tag_func
98
99 integer*4, dimension(nn,nn,nn) :: yon
100
101 real*8 :: alfa11,alfa12,alfa13,alfa21,alfa22,alfa23,alfa31,alfa32,alfa33
102 real*8 :: beta11,beta12,beta13,beta21,beta22,beta23,beta31,beta32,beta33
103 real*8 :: gamma1,gamma2,gamma3,delta1,delta2,delta3
104 real*8 :: dxdx,dxdy,dxdz,dydx,dydy,dydz,dzdx,dzdy,dzdz,det_j
105 real*8 :: depth
106 real*8 :: lambda, mu
107 real*8 :: t_stress
108 real*8 :: get_func_value
109 real*8 :: csi,fpeak,q,pi,gamma
110
111 real*8, dimension(nfdata) :: func_data
112 real*8, dimension(nn) :: ct,ww
113 real*8, dimension(nnod_loc) :: vs_tria
114
115 real*8, dimension(nn,nn) :: dd
116
117 real*8, dimension(nn,nn,nn) :: rho,gamma0,gamma_new
118 real*8, dimension(nn,nn,nn) :: mc_el,mck_el
119 real*8, dimension(nn,nn,nn) :: r_el
120
121
122 gamma_new = gamma0;
123
124! ELEMENT NODAL MASS CALCULATION
125
126 do k = 1,nn
127 do j = 1,nn
128 do i = 1,nn
129
130 if (yon(i,j,k).eq.1) then
131
132 dxdx = alfa11 +beta12*ct(k) +beta13*ct(j) &
133 + gamma1*ct(j)*ct(k)
134 dydx = alfa21 +beta22*ct(k) +beta23*ct(j) &
135 + gamma2*ct(j)*ct(k)
136 dzdx = alfa31 +beta32*ct(k) +beta33*ct(j) &
137 + gamma3*ct(j)*ct(k)
138 dxdy = alfa12 +beta11*ct(k) +beta13*ct(i) &
139 + gamma1*ct(k)*ct(i)
140 dydy = alfa22 +beta21*ct(k) +beta23*ct(i) &
141 + gamma2*ct(k)*ct(i)
142 dzdy = alfa32 +beta31*ct(k) +beta33*ct(i) &
143 + gamma3*ct(k)*ct(i)
144 dxdz = alfa13 +beta11*ct(j) +beta12*ct(i) &
145 + gamma1*ct(i)*ct(j)
146 dydz = alfa23 +beta21*ct(j) +beta22*ct(i) &
147 + gamma2*ct(i)*ct(j)
148 dzdz = alfa33 +beta31*ct(j) +beta32*ct(i) &
149 + gamma3*ct(i)*ct(j)
150 det_j = dxdz * (dydx*dzdy - dzdx*dydy) &
151 - dydz * (dxdx*dzdy - dzdx*dxdy) &
152 + dzdz * (dxdx*dydy - dydx*dxdy)
153
154 is = nn*nn*(k -1) +nn*(j -1) +i
155 in = cs(cs(ielem -1) +is)
156 !-------------------------------------------------------------------
157 !
158 if(tcase .ne. 16) then
159 do fn = 1,nf
160 if (vcase.eq.tag_func(fn)) then
161 if (func_type(fn) .eq. 61) then
162 csi = get_func_value(nf,func_type,func_indx,func_data, nfdata,&
163 fn,r_el(i,j,k))
164
165
166 q = 1.d0/(2.d0*csi)
167 pi = 4.0d0 * datan(1.0d0)
168 gamma = pi * fpeak / q
169
170
171 if (gamma0(i,j,k) .ge. gamma) gamma = gamma0(i,j,k)
172 gamma_new(i,j,k) = gamma
173
174 mc_el(i,j,k) = 2.d0 * gamma * rho(i,j,k) &
175 * det_j * ww(i) * ww(j) * ww(k)
176 mck_el(i,j,k) = (gamma**2) * rho(i,j,k) &
177 * det_j * ww(i) * ww(j) * ww(k)
178
179 endif
180 endif
181 enddo
182 elseif (tcase .eq. 16) then
183
184 do fn = 1,nf
185 if (vcase .eq. tag_func(fn)) then
186 if (func_type(fn) .eq. 61 .and. vs_tria(in) .lt. 325.d0) then
187
188 csi = get_func_value(nf,func_type,func_indx,func_data, nfdata,&
189 fn,r_el(i,j,k))
190
191
192 q = 1.d0/(2.d0*csi)
193 pi = 4.0d0 * datan(1.0d0)
194 gamma = pi * fpeak / q
195
196
197 if (gamma0(i,j,k) .ge. gamma) gamma = gamma0(i,j,k)
198 gamma_new(i,j,k) = gamma
199
200 mc_el(i,j,k) = 2.d0 * gamma * rho(i,j,k) &
201 * det_j * ww(i) * ww(j) * ww(k)
202 mck_el(i,j,k) = (gamma**2) * rho(i,j,k) &
203 * det_j * ww(i) * ww(j) * ww(k)
204
205
206 elseif (func_type(fn) .eq. 63 .and. vs_tria(in) .lt. 450.d0) then
207 csi = get_func_value(nf,func_type,func_indx,func_data, nfdata,&
208 fn,r_el(i,j,k))
209
210
211 q = 1.d0/(2.d0*csi)
212 pi = 4.0d0 * datan(1.0d0)
213 gamma = pi * fpeak / q
214
215
216 if (gamma0(i,j,k) .ge. gamma) gamma = gamma0(i,j,k)
217 gamma_new(i,j,k) = gamma
218
219 mc_el(i,j,k) = 2.d0 * gamma * rho(i,j,k) &
220 * det_j * ww(i) * ww(j) * ww(k)
221 mck_el(i,j,k) = (gamma**2) * rho(i,j,k) &
222 * det_j * ww(i) * ww(j) * ww(k)
223
224
225 endif
226 endif
227 enddo
228
229 endif!tcase = 16
230
231 !
232 !-------------------------------------------------------------------
233
234 endif !if (yon(i,j,k).eq.1) then
235
236 enddo
237 enddo
238 enddo
239
240
241 return
242
243 end subroutine make_damping_matrix_nle
244
get_func_value
real *8 function get_func_value(nb_fnc, type_fnc, ind_fnc, data_fnc, nb_data_fnc, id_fnc, time
Computes time evolution function.
Definition GET_FUNC_VALUE.f90:38
make_damping_matrix_nle
subroutine make_damping_matrix_nle(nn, ct, ww, dd, rho, gamma0, alfa11, alfa12, alfa13, alfa21, alfa22, alfa23, alfa31, alfa32, alfa33, beta11, beta12, beta13, beta21, beta22, beta23, beta31, beta32, beta33, gamma1, gamma2, gamma3, delta1, delta2, delta3, mc_el, mck_el, cs_nnz, cs, ielem, vcase, r_el, fpeak, func_type, func_indx, func_data, nfdata, nf, t_str
Make damping matrices for non-linear elastic materials.
Definition MAKE_DAMPING_MATRIX_NLE.f90:82