SPEED
MAKE_DIV_ROT.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
19
59
60 subroutine make_div_rot(nn,ct,ww,dd,&
61 alfa11,alfa12,alfa13,alfa21,alfa22,alfa23,&
62 alfa31,alfa32,alfa33,beta11,beta12,beta13,&
63 beta21,beta22,beta23,beta31,beta32,beta33,&
64 gamma1,gamma2,gamma3,delta1,delta2,delta3,&
65 ux,uy,uz,div,rotx,roty,rotz)
66
67 implicit none
68
69 integer*4 :: nn
70 integer*4 :: i,j,k,p,q,r
71
72 real*8 :: alfa11,alfa12,alfa13,alfa21,alfa22,alfa23,alfa31,alfa32,alfa33
73 real*8 :: beta11,beta12,beta13,beta21,beta22,beta23,beta31,beta32,beta33
74 real*8 :: gamma1,gamma2,gamma3,delta1,delta2,delta3
75 real*8 :: dxdx,dxdy,dxdz,dydx,dydy,dydz,dzdx,dzdy,dzdz,det_j
76 real*8 :: duxdx,duxdy,duxdz,duydx,duydy,duydz,duzdx,duzdy,duzdz
77 real*8 :: t1ux,t2ux,t3ux,t1uy,t2uy,t3uy,t1uz,t2uz,t3uz
78
79 real*8, dimension(nn) :: ct,ww
80
81 real*8, dimension(nn,nn) :: dd
82
83 real*8, dimension(nn,nn,nn) :: ux,uy,uz
84 real*8, dimension(nn,nn,nn) :: div,rotx,roty,rotz
85
86! STRESS CALCULATION
87
88 do r = 1,nn
89 do q = 1,nn
90 do p = 1,nn
91 t1ux = 0.d0
92 t1uy = 0.d0
93 t1uz = 0.d0
94 t2ux = 0.d0
95 t2uy = 0.d0
96 t2uz = 0.d0
97 t3ux = 0.d0
98 t3uy = 0.d0
99 t3uz = 0.d0
100
101 dxdx = alfa11 + beta12*ct(r) + beta13*ct(q) &
102 + gamma1*ct(q)*ct(r)
103 dydx = alfa21 + beta22*ct(r) + beta23*ct(q) &
104 + gamma2*ct(q)*ct(r)
105 dzdx = alfa31 + beta32*ct(r) + beta33*ct(q) &
106 + gamma3*ct(q)*ct(r)
107
108 dxdy = alfa12 + beta11*ct(r) + beta13*ct(p) &
109 + gamma1*ct(r)*ct(p)
110 dydy = alfa22 + beta21*ct(r) + beta23*ct(p) &
111 + gamma2*ct(r)*ct(p)
112 dzdy = alfa32 + beta31*ct(r) + beta33*ct(p) &
113 + gamma3*ct(r)*ct(p)
114
115 dxdz = alfa13 + beta11*ct(q) + beta12*ct(p) &
116 + gamma1*ct(p)*ct(q)
117 dydz = alfa23 + beta21*ct(q) + beta22*ct(p) &
118 + gamma2*ct(p)*ct(q)
119 dzdz = alfa33 + beta31*ct(q) + beta32*ct(p) &
120 + gamma3*ct(p)*ct(q)
121
122 det_j = dxdz * (dydx*dzdy - dzdx*dydy) &
123 - dydz * (dxdx*dzdy - dzdx*dxdy) &
124 + dzdz * (dxdx*dydy - dydx*dxdy)
125
126
127 do i = 1,nn
128 t1ux = t1ux + ux(i,q,r) * dd(p,i)
129 t1uy = t1uy + uy(i,q,r) * dd(p,i)
130 t1uz = t1uz + uz(i,q,r) * dd(p,i)
131 enddo
132
133 do j = 1,nn
134 t2ux = t2ux + ux(p,j,r) * dd(q,j)
135 t2uy = t2uy + uy(p,j,r) * dd(q,j)
136 t2uz = t2uz + uz(p,j,r) * dd(q,j)
137 enddo
138
139 do k = 1,nn
140 t3ux = t3ux + ux(p,q,k) * dd(r,k)
141 t3uy = t3uy + uy(p,q,k) * dd(r,k)
142 t3uz = t3uz + uz(p,q,k) * dd(r,k)
143 enddo
144
145
146 duxdx = 1.0d0 / det_j *(&
147 ((dydy*dzdz - dydz*dzdy) * t1ux) + &
148 ((dydz*dzdx - dydx*dzdz) * t2ux) + &
149 ((dydx*dzdy - dydy*dzdx) * t3ux))
150
151 duydx = 1.0d0 / det_j *(&
152 ((dydy*dzdz - dydz*dzdy) * t1uy) + &
153 ((dydz*dzdx - dydx*dzdz) * t2uy) + &
154 ((dydx*dzdy - dydy*dzdx) * t3uy))
155
156 duzdx = 1.0d0 / det_j *(&
157 ((dydy*dzdz - dydz*dzdy) * t1uz) + &
158 ((dydz*dzdx - dydx*dzdz) * t2uz) + &
159 ((dydx*dzdy - dydy*dzdx) * t3uz))
160
161 duxdy = 1.0d0 / det_j *(&
162 ((dzdy*dxdz - dzdz*dxdy) * t1ux) + &
163 ((dzdz*dxdx - dzdx*dxdz) * t2ux) + &
164 ((dzdx*dxdy - dzdy*dxdx) * t3ux))
165
166 duydy = 1.0d0 / det_j *(&
167 ((dzdy*dxdz - dzdz*dxdy) * t1uy) + &
168 ((dzdz*dxdx - dzdx*dxdz) * t2uy) + &
169 ((dzdx*dxdy - dzdy*dxdx) * t3uy))
170
171 duzdy = 1.0d0 / det_j *(&
172 ((dzdy*dxdz - dzdz*dxdy) * t1uz) + &
173 ((dzdz*dxdx - dzdx*dxdz) * t2uz) + &
174 ((dzdx*dxdy - dzdy*dxdx) * t3uz))
175
176 duxdz = 1.0d0 / det_j *(&
177 ((dxdy*dydz - dxdz*dydy) * t1ux) + &
178 ((dxdz*dydx - dxdx*dydz) * t2ux) + &
179 ((dxdx*dydy - dxdy*dydx) * t3ux))
180
181 duydz = 1.0d0 / det_j *(&
182 ((dxdy*dydz - dxdz*dydy) * t1uy) + &
183 ((dxdz*dydx - dxdx*dydz) * t2uy) + &
184 ((dxdx*dydy - dxdy*dydx) * t3uy))
185
186 duzdz = 1.0d0 / det_j *(&
187 ((dxdy*dydz - dxdz*dydy) * t1uz) + &
188 ((dxdz*dydx - dxdx*dydz) * t2uz) + &
189 ((dxdx*dydy - dxdy*dydx) * t3uz))
190
191
192 div(p,q,r) = duxdx +duydy +duzdz
193 rotx(p,q,r) = duzdy -duydz
194 roty(p,q,r) = duxdz -duzdx
195 rotz(p,q,r) = duydx -duxdy
196
197 enddo
198 enddo
199 enddo
200
201 return
202
203 end subroutine make_div_rot
204
make_div_rot
subroutine make_div_rot(nn, ct, ww, dd, alfa11, alfa12, alfa13, alfa21, alfa22, alfa23, alfa31, alfa32, alfa33, beta11, beta12, beta13, beta21, beta22, beta23, beta31, beta32, beta33, gamma1, gamma2, gamma3, delta1, delta2, delta3, ux, uy, uz, div, rotx, roty, rotz)
Computes div(u) and rot(u).
Definition MAKE_DIV_ROT.f90:66