SPEED
MAKE_ABC_FORCE.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
69
70
71 subroutine make_abc_force(nn,xq,wq,dd,rho,lambda,mu,&
72 AA11,AA12,AA13,AA21,AA22,AA23,&
73 AA31,AA32,AA33,BB11,BB12,BB13,&
74 BB21,BB22,BB23,BB31,BB32,BB33,&
75 GG1,GG2,GG3,DD1,DD2,DD3,&
76 ia,ib,ja,jb,ka,kb,ux,uy,uz,vx,vy,vz,fx,fy,fz)
77
78!
79!**************************************************************************************************
80
81 implicit none
82
83 integer*4 :: nn
84 integer*4 :: ia,ib,ja,jb,ka,kb
85 integer*4 :: i,j,k,p,q,r
86
87 real*8 :: aa11,aa12,aa13,aa21,aa22,aa23,aa31,aa32,aa33
88 real*8 :: bb11,bb12,bb13,bb21,bb22,bb23,bb31,bb32,bb33
89 real*8 :: gg1,gg2,gg3,dd1,dd2,dd3
90 real*8 :: dxdx,dxdy,dxdz,dydx,dydy,dydz,dzdx,dzdy,dzdz,det_j
91 real*8 :: duxdx,duxdy,duxdz,duydx,duydy,duydz,duzdx,duzdy,duzdz
92 real*8 :: dut1dt1,dut2dt2,dut3dt1,dut3dt2
93 real*8 :: dxdu,dxdv,dydu,dydv,dzdu,dzdv
94 real*8 :: s1,s2,s3,vt1,vt2,vt3,r8t
95 real*8 :: t1x,t1y,t1z,t2x,t2y,t2z,t3x,t3y,t3z
96 real*8 :: t1ux,t2ux,t3ux,t1uy,t2uy,t3uy,t1uz,t2uz,t3uz
97
98 real*8, dimension(nn) :: xq,wq
99
100 real*8, dimension(nn,nn) :: dd
101
102 real*8, dimension(nn,nn,nn) :: rho,lambda,mu
103 real*8, dimension(nn,nn,nn) :: ux,uy,uz
104 real*8, dimension(nn,nn,nn) :: vx,vy,vz
105 real*8, dimension(nn,nn,nn) :: fx,fy,fz
106 real*8, dimension(nn,nn,nn) :: alfa,beta
107
108 do r = 1,nn
109 do q = 1,nn
110 do p = 1,nn
111 alfa(p,q,r) = dsqrt((lambda(p,q,r) + 2.0d0*mu(p,q,r))/rho(p,q,r))
112 beta(p,q,r) = dsqrt(mu(p,q,r)/rho(p,q,r))
113 enddo
114 enddo
115 enddo
116
117 do r = 1,nn
118 do q = 1,nn
119 do p = 1,nn
120 fx(p,q,r) = 0.0d0; fy(p,q,r) = 0.0d0; fz(p,q,r) = 0.0d0
121 enddo
122 enddo
123 enddo
124
125 do r = ka,kb
126 do q = ja,jb
127 do p = ia,ib
128 t1ux = 0.d0; t1uy = 0.d0; t1uz = 0.d0
129 t2ux = 0.d0; t2uy = 0.d0; t2uz = 0.d0
130 t3ux = 0.d0; t3uy = 0.d0; t3uz = 0.d0
131
132 dxdx = aa11 + bb12*xq(r) + bb13*xq(q) + gg1*xq(q)*xq(r)
133 dydx = aa21 + bb22*xq(r) + bb23*xq(q) + gg2*xq(q)*xq(r)
134 dzdx = aa31 + bb32*xq(r) + bb33*xq(q) + gg3*xq(q)*xq(r)
135
136 dxdy = aa12 + bb11*xq(r) + bb13*xq(p) + gg1*xq(r)*xq(p)
137 dydy = aa22 + bb21*xq(r) + bb23*xq(p) + gg2*xq(r)*xq(p)
138 dzdy = aa32 + bb31*xq(r) + bb33*xq(p) + gg3*xq(r)*xq(p)
139
140 dxdz = aa13 + bb11*xq(q) + bb12*xq(p) + gg1*xq(p)*xq(q)
141 dydz = aa23 + bb21*xq(q) + bb22*xq(p) + gg2*xq(p)*xq(q)
142 dzdz = aa33 + bb31*xq(q) + bb32*xq(p) + gg3*xq(p)*xq(q)
143
144 det_j = dxdz * (dydx*dzdy - dzdx*dydy) &
145 - dydz * (dxdx*dzdy - dzdx*dxdy) &
146 + dzdz * (dxdx*dydy - dydx*dxdy)
147
148
149 do i = 1,nn
150 t1ux = t1ux + ux(i,q,r) * dd(p,i)
151 t1uy = t1uy + uy(i,q,r) * dd(p,i)
152 t1uz = t1uz + uz(i,q,r) * dd(p,i)
153 enddo
154
155 do j = 1,nn
156 t2ux = t2ux + ux(p,j,r) * dd(q,j)
157 t2uy = t2uy + uy(p,j,r) * dd(q,j)
158 t2uz = t2uz + uz(p,j,r) * dd(q,j)
159 enddo
160
161 do k = 1,nn
162 t3ux = t3ux + ux(p,q,k) * dd(r,k)
163 t3uy = t3uy + uy(p,q,k) * dd(r,k)
164 t3uz = t3uz + uz(p,q,k) * dd(r,k)
165 enddo
166
167
168 duxdx = 1.0d0 / det_j *(&
169 ((dydy*dzdz - dydz*dzdy) * t1ux) + &
170 ((dydz*dzdx - dydx*dzdz) * t2ux) + &
171 ((dydx*dzdy - dydy*dzdx) * t3ux))
172
173 duydx = 1.0d0 / det_j *(&
174 ((dydy*dzdz - dydz*dzdy) * t1uy) + &
175 ((dydz*dzdx - dydx*dzdz) * t2uy) + &
176 ((dydx*dzdy - dydy*dzdx) * t3uy))
177
178 duzdx = 1.0d0 / det_j *(&
179 ((dydy*dzdz - dydz*dzdy) * t1uz) + &
180 ((dydz*dzdx - dydx*dzdz) * t2uz) + &
181 ((dydx*dzdy - dydy*dzdx) * t3uz))
182
183 duxdy = 1.0d0 / det_j *(&
184 ((dzdy*dxdz - dzdz*dxdy) * t1ux) + &
185 ((dzdz*dxdx - dzdx*dxdz) * t2ux) + &
186 ((dzdx*dxdy - dzdy*dxdx) * t3ux))
187
188 duydy = 1.0d0 / det_j *(&
189 ((dzdy*dxdz - dzdz*dxdy) * t1uy) + &
190 ((dzdz*dxdx - dzdx*dxdz) * t2uy) + &
191 ((dzdx*dxdy - dzdy*dxdx) * t3uy))
192
193 duzdy = 1.0d0 / det_j *(&
194 ((dzdy*dxdz - dzdz*dxdy) * t1uz) + &
195 ((dzdz*dxdx - dzdx*dxdz) * t2uz) + &
196 ((dzdx*dxdy - dzdy*dxdx) * t3uz))
197
198 duxdz = 1.0d0 / det_j *(&
199 ((dxdy*dydz - dxdz*dydy) * t1ux) + &
200 ((dxdz*dydx - dxdx*dydz) * t2ux) + &
201 ((dxdx*dydy - dxdy*dydx) * t3ux))
202
203 duydz = 1.0d0 / det_j *(&
204 ((dxdy*dydz - dxdz*dydy) * t1uy) + &
205 ((dxdz*dydx - dxdx*dydz) * t2uy) + &
206 ((dxdx*dydy - dxdy*dydx) * t3uy))
207
208 duzdz = 1.0d0 / det_j *(&
209 ((dxdy*dydz - dxdz*dydy) * t1uz) + &
210 ((dxdz*dydx - dxdx*dydz) * t2uz) + &
211 ((dxdx*dydy - dxdy*dydx) * t3uz))
212
213 if ((ia.eq.ib).and.(ia.eq.1)) then
214 dxdu = dxdy
215 dydu = dydy
216 dzdu = dzdy
217
218 dxdv = -dxdz
219 dydv = -dydz
220 dzdv = -dzdz
221 endif
222
223 if ((ja.eq.jb).and.(ja.eq.1)) then
224 dxdu = dxdz
225 dydu = dydz
226 dzdu = dzdz
227
228 dxdv = -dxdx
229 dydv = -dydx
230 dzdv = -dzdx
231 endif
232
233 if ((ka.eq.kb).and.(ka.eq.1)) then
234 dxdu = dxdx
235 dydu = dydx
236 dzdu = dzdx
237
238 dxdv = -dxdy
239 dydv = -dydy
240 dzdv = -dzdy
241 endif
242
243 if ((ia.eq.ib).and.(ia.eq.nn)) then
244 dxdu = dxdy
245 dydu = dydy
246 dzdu = dzdy
247
248 dxdv = dxdz
249 dydv = dydz
250 dzdv = dzdz
251 endif
252
253 if ((ja.eq.jb).and.(ja.eq.nn)) then
254 dxdu = dxdz
255 dydu = dydz
256 dzdu = dzdz
257
258 dxdv = dxdx
259 dydv = dydx
260 dzdv = dzdx
261 endif
262
263 if ((ka.eq.kb).and.(ka.eq.nn)) then
264 dxdu = dxdx
265 dydu = dydx
266 dzdu = dzdx
267
268 dxdv = dxdy
269 dydv = dydy
270 dzdv = dzdy
271 endif
272
273 t1x = dxdu
274 t1y = dydu
275 t1z = dzdu
276
277 t2x = dxdv
278 t2y = dydv
279 t2z = dzdv
280
281 t3x = t1y*t2z - t1z*t2y
282 t3y = t1z*t2x - t1x*t2z
283 t3z = t1x*t2y - t1y*t2x
284
285 r8t = dsqrt(t1x*t1x +t1y*t1y +t1z*t1z)
286 t1x = t1x/r8t
287 t1y = t1y/r8t
288 t1z = t1z/r8t
289
290 r8t = dsqrt(t3x*t3x +t3y*t3y +t3z*t3z)
291 t3x = t3x/r8t
292 t3y = t3y/r8t
293 t3z = t3z/r8t
294
295 t2x = t3y*t1z - t3z*t1y
296 t2y = t3z*t1x - t3x*t1z
297 t2z = t3x*t1y - t3y*t1x
298
299
300 dut1dt1 = t1x*t1x*duxdx + t1x*t1y*duydx + t1x*t1z*duzdx &
301 + t1y*t1x*duxdy + t1y*t1y*duydy + t1y*t1z*duzdy &
302 + t1z*t1x*duxdz + t1z*t1y*duydz + t1z*t1z*duzdz
303
304 dut2dt2 = t2x*t2x*duxdx + t2x*t2y*duydx + t2x*t2z*duzdx &
305 + t2y*t2x*duxdy + t2y*t2y*duydy + t2y*t2z*duzdy &
306 + t2z*t2x*duxdz + t2z*t2y*duydz + t2z*t2z*duzdz
307
308 dut3dt1 = t1x*t3x*duxdx + t1x*t3y*duydx + t1x*t3z*duzdx &
309 + t1y*t3x*duxdy + t1y*t3y*duydy + t1y*t3z*duzdy &
310 + t1z*t3x*duxdz + t1z*t3y*duydz + t1z*t3z*duzdz
311
312 dut3dt2 = t2x*t3x*duxdx + t2x*t3y*duydx + t2x*t3z*duzdx &
313 + t2y*t3x*duxdy + t2y*t3y*duydy + t2y*t3z*duzdy &
314 + t2z*t3x*duxdz + t2z*t3y*duydz + t2z*t3z*duzdz
315
316 vt1 = t1x*vx(p,q,r) + t1y*vy(p,q,r) + t1z*vz(p,q,r)
317 vt2 = t2x*vx(p,q,r) + t2y*vy(p,q,r) + t2z*vz(p,q,r)
318 vt3 = t3x*vx(p,q,r) + t3y*vy(p,q,r) + t3z*vz(p,q,r)
319
320 s1 = (mu(p,q,r)*(2.0d0*beta(p,q,r) -alfa(p,q,r))) &
321 /beta(p,q,r) * dut3dt1 - mu(p,q,r)/beta(p,q,r) * vt1
322 s2 = (mu(p,q,r)*(2.0d0*beta(p,q,r) -alfa(p,q,r))) &
323 /beta(p,q,r) * dut3dt2 - mu(p,q,r)/beta(p,q,r) * vt2
324 s3 = (lambda(p,q,r)*beta(p,q,r) +2.0d0*mu(p,q,r)* &
325 (beta(p,q,r) -alfa(p,q,r)))/alfa(p,q,r) &
326 * (dut1dt1 +dut2dt2) - (lambda(p,q,r) +2.0d0*mu(p,q,r))/alfa(p,q,r) * vt3
327
328 det_j = dsqrt(((dxdu*dxdu +dydu*dydu +dzdu*dzdu) &
329 * (dxdv*dxdv +dydv*dydv +dzdv*dzdv)) &
330 -((dxdu*dxdv +dydu*dydv +dzdu*dzdv) &
331 * (dxdu*dxdv +dydu*dydv +dzdu*dzdv)))
332
333 fx(p,q,r) = fx(p,q,r) - wq(p)*wq(q)*wq(r)*det_j/wq(1) &
334 * (s1*t1x + s2*t2x + s3*t3x)
335 fy(p,q,r) = fy(p,q,r) - wq(p)*wq(q)*wq(r)*det_j/wq(1) &
336 * (s1*t1y + s2*t2y + s3*t3y)
337 fz(p,q,r) = fz(p,q,r) - wq(p)*wq(q)*wq(r)*det_j/wq(1) &
338 * (s1*t1z + s2*t2z + s3*t3z)
339
340 enddo
341 enddo
342 enddo
343
344 return
345
346 end subroutine make_abc_force
347
make_abc_force
subroutine make_abc_force(nn, xq, wq, dd, rho, lambda, mu, aa11, aa12, aa13, aa21, aa22, aa23, aa31, aa32, aa33, bb11, bb12, bb13, bb21, bb22, bb23, bb31, bb32, bb33, gg1, gg2, gg3, dd1, dd2, dd3, ia, ib, ja, jb, ka, kb, ux, uy, uz, vx, vy, vz, fx, f
Computes ABC's.
Definition MAKE_ABC_FORCE.f90:77