SPEED
Functions/Subroutines
MAKE_DAMPING_MATRIX_NLE.f90 File Reference

Go to the source code of this file.

Functions/Subroutines

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.
 

Function/Subroutine Documentation

◆ make_damping_matrix_nle()

subroutine make_damping_matrix_nle ( 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, dimension(nn,nn,nn)  gamma0,
real*8  alfa11,
real*8  alfa12,
real*8  alfa13,
real*8  alfa21,
real*8  alfa22,
real*8  alfa23,
real*8  alfa31,
real*8  alfa32,
  alfa33,
real*8  beta11,
real*8  beta12,
real*8  beta13,
real*8  beta21,
real*8  beta22,
real*8  beta23,
real*8  beta31,
real*8  beta32,
  beta33,
real*8  gamma1,
real*8  gamma2,
real*8  gamma3,
real*8  delta1,
real*8  delta2,
real*8  delta3,
real*8, dimension(nn,nn,nn)  mc_el,
real*8, dimension(nn,nn,nn)  mck_el,
integer*4  cs_nnz,
integer*4, dimension(0:cs_nnz)  cs,
integer*4  ielem,
integer*4  vcase,
real*8, dimension(nn,nn,nn)  r_el,
real*8  fpeak,
integer*4, dimension(nf)  func_type,
integer*4, dimension(nf +1)  func_indx,
real*8, dimension(nfdata)  func_data,
integer*4  nfdata,
integer*4  nf,
  t_str 
)

Make damping matrices for non-linear elastic materials.

Author
Ilario Mazzieri
Date
November, 2014
Version
1.2
Parameters
[in]nnnuber of 1D GLL nodes
[in]ct1D GLL nodes
[in]ww1D GLL weights
[in]ddspectral derivative matrix
[in]rhomaterial density
[in]gamma0damping coefficient
[in]alfa11costant values for the bilinear map
[in]alfa12costant values for the bilinear map
[in]alfa13costant values for the bilinear map
[in]alfa21costant values for the bilinear map
[in]alfa22costant values for the bilinear map
[in]alfa23costant values for the bilinear map
[in]alfa31costant values for the bilinear map
[in]alfa32costant values for the bilinear map
[in]alfa33costant values for the bilinear map
[in]beta11costant values for the bilinear map
[in]beta12costant values for the bilinear map
[in]beta13costant values for the bilinear map
[in]beta21costant values for the bilinear map
[in]beta22costant values for the bilinear map
[in]beta23costant values for the bilinear map
[in]beta31costant values for the bilinear map
[in]beta32costant values for the bilinear map
[in]beta33costant values for the bilinear map
[in]gamma1costant values for the bilinear map
[in]gamma2costant values for the bilinear map
[in]gamma3costant values for the bilinear map
[in]delta1costant values for the bilinear map
[in]delta2costant values for the bilinear map
[in]delta3costant values for the bilinear map
[in]cs_nnzlength of cs
[in]csconnectivity vector
[in]vcasedepth for non linear block
[in]R_elreduction factor associated to strain invariant
[in]fpeakpeak frequency
[in]nfnumber of functions
[in]func_typefunction type
[in]func_indxindices for the data
[in]nfdatanumber of data for each function
[in]func_datadata for the calculation (depending on type)
[in]t_stresstime
[in]tag_funcfunction label
[in]yon1 if the node is non-linear elastic, 0 otherwise
[in]tcasenot-honoring test case
[in]nnod_locnumber of local nodes
[in]vs_triavs values for each node
[out]mc_elnon-linear damping matrix
[out]mck_elnon-linear damping matrix
[out]gamma_newgamma value for debug mode

Definition at line 73 of file MAKE_DAMPING_MATRIX_NLE.f90.

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
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

References get_func_value().

Here is the call graph for this function: