SPEED
MAKE_STRAIN_ROT_STRESS_TENSOR.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_strain_rot_stress_tensor(nn,ct,ww,dd,&
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 ux,uy,uz,div,&
79 lambda,mu,&
80 rotx,roty,rotz,&
81 strainxx,strainyy,strainzz,&
82 strainxy,strainyz,strainzx,&
83 stressxx,stressyy,stresszz,&
84 stressxy,stressyz,stresszx)
85
86
87
88 implicit none
89
90 integer*4 :: nn
91 integer*4 :: i,j,k,p,q,r,L
92
93 real*8 :: alfa11,alfa12,alfa13,alfa21,alfa22,alfa23,alfa31,alfa32,alfa33
94 real*8 :: beta11,beta12,beta13,beta21,beta22,beta23,beta31,beta32,beta33
95 real*8 :: gamma1,gamma2,gamma3,delta1,delta2,delta3,dphi,phi
96 real*8 :: dxdx,dxdy,dxdz,dydx,dydy,dydz,dzdx,dzdy,dzdz,det_j
97 real*8 :: duxdx,duxdy,duxdz,duydx,duydy,duydz,duzdx,duzdy,duzdz
98 real*8 :: t1ux,t2ux,t3ux,t1uy,t2uy,t3uy,t1uz,t2uz,t3uz, dpdcsi, dpdeta, dpdzeta
99
100 real*8, dimension(nn) :: ct,ww
101
102 real*8, dimension(nn,nn) :: dd
103
104 real*8, dimension(nn,nn,nn) :: ux,uy,uz
105 real*8, dimension(nn,nn,nn) :: div,rotx,roty,rotz
106 real*8, dimension(nn,nn,nn) :: strainxx,strainyy,strainzz
107 real*8, dimension(nn,nn,nn) :: strainxy,strainyz,strainzx
108 real*8, dimension(nn,nn,nn) :: stressxx,stressyy,stresszz
109 real*8, dimension(nn,nn,nn) :: stressxy,stressyz,stresszx
110 real*8, dimension(nn,nn,nn) :: lambda,mu
111
112 real*8, dimension(nn**3) :: derx, dery, derz, ux_el, uy_el, uz_el
113
114! STRESS CALCULATION
115
116
117
118 do r = 1,nn
119 do q = 1,nn
120 do p = 1,nn
121 !M = p + (q-1)*nn + (r-1)*nn*nn
122
123 do k = 1, nn
124 do j = 1, nn
125 do i = 1, nn
126 l = i + (j-1)*nn + (k-1)*nn*nn
127
128 ux_el(l) = ux(i,j,k)
129 uy_el(l) = uy(i,j,k)
130 uz_el(l) = uz(i,j,k)
131
132
133 dxdx = alfa11 + beta12*ct(r) + beta13*ct(q) &
134 + gamma1*ct(q)*ct(r)
135 dydx = alfa21 + beta22*ct(r) + beta23*ct(q) &
136 + gamma2*ct(q)*ct(r)
137 dzdx = alfa31 + beta32*ct(r) + beta33*ct(q) &
138 + gamma3*ct(q)*ct(r)
139
140 dxdy = alfa12 + beta11*ct(r) + beta13*ct(p) &
141 + gamma1*ct(r)*ct(p)
142 dydy = alfa22 + beta21*ct(r) + beta23*ct(p) &
143 + gamma2*ct(r)*ct(p)
144 dzdy = alfa32 + beta31*ct(r) + beta33*ct(p) &
145 + gamma3*ct(r)*ct(p)
146
147 dxdz = alfa13 + beta11*ct(q) + beta12*ct(p) &
148 + gamma1*ct(p)*ct(q)
149 dydz = alfa23 + beta21*ct(q) + beta22*ct(p) &
150 + gamma2*ct(p)*ct(q)
151 dzdz = alfa33 + beta31*ct(q) + beta32*ct(p) &
152 + gamma3*ct(p)*ct(q)
153
154 det_j = dxdz * (dydx*dzdy - dzdx*dydy) &
155 - dydz * (dxdx*dzdy - dzdx*dxdy) &
156 + dzdz * (dxdx*dydy - dydx*dxdy)
157
158
159 dpdcsi = dphi(nn-1, ct(p), ct(i)) * phi(nn-1, ct(q), ct(j)) * phi(nn-1, ct(r), ct(k))
160 dpdeta = phi(nn-1, ct(p), ct(i)) * dphi(nn-1, ct(q), ct(j)) * phi(nn-1, ct(r), ct(k))
161 dpdzeta = phi(nn-1, ct(p), ct(i)) * phi(nn-1, ct(q), ct(j)) * dphi(nn-1, ct(r), ct(k))
162
163
164 derx(l) = (dydx*(dzdy*dpdzeta - dzdz*dpdeta) - dydy*(dzdx*dpdzeta - dzdz*dpdcsi) + &
165 dydz*(dzdx*dpdeta - dzdy*dpdcsi))/det_j
166
167 dery(l) = (dzdx*(dxdy*dpdzeta - dxdz*dpdeta) - dzdy*(dxdx*dpdzeta - dxdz*dpdcsi) + &
168 dzdz*(dxdx*dpdeta - dxdy*dpdcsi))/det_j
169
170 derz(l) = (dxdx*(dydy*dpdzeta - dydz*dpdeta) - dxdy*(dydx*dpdzeta - dydz*dpdcsi) + &
171 dxdz*(dydx*dpdeta - dydy*dpdcsi))/det_j
172
173
174 enddo
175 enddo
176 enddo
177
178
179
180 strainxx(p,q,r) = dot_product(derx,ux_el)
181 strainyy(p,q,r) = dot_product(dery,uy_el)
182 strainyy(p,q,r) = dot_product(derz,uz_el)
183 strainxy(p,q,r) = 0.5d0*(dot_product(derx,uy_el) + dot_product(dery,ux_el))
184 strainyz(p,q,r) = 0.5d0*(dot_product(derz,uy_el) + dot_product(dery,uz_el))
185 strainzx(p,q,r) = 0.5d0*(dot_product(derx,uz_el) + dot_product(derz,ux_el))
186
187 rotx(p,q,r) = 0.5d0*(dot_product(dery,uz_el) - dot_product(derz,uy_el))
188 roty(p,q,r) = 0.5d0*(dot_product(derz,ux_el) - dot_product(derx,uz_el))
189 rotz(p,q,r) = 0.5d0*(dot_product(derx,uy_el) - dot_product(dery,ux_el))
190
191 stressxx(p,q,r) = lambda(p,q,r) * (dot_product(derx,ux_el) + dot_product(dery,uy_el) &
192 + dot_product(derz,uz_el)) + 2.0d0*mu(p,q,r) * dot_product(derx,ux_el)
193 stressyy(p,q,r) = lambda(p,q,r) * (dot_product(derx,ux_el) + dot_product(dery,uy_el) &
194 + dot_product(derz,uz_el)) + 2.0d0*mu(p,q,r) * dot_product(dery,uy_el)
195 stresszz(p,q,r) = lambda(p,q,r) * (dot_product(derx,ux_el) + dot_product(dery,uy_el) &
196 + dot_product(derz,uz_el)) + 2.0d0*mu(p,q,r) * dot_product(derz,uz_el)
197
198 stressyz(p,q,r) = mu(p,q,r) * (dot_product(derz,uy_el) + dot_product(dery,uz_el))
199 stresszx(p,q,r) = mu(p,q,r) * (dot_product(derx,uz_el) + dot_product(derz,ux_el))
200 stressxy(p,q,r) = mu(p,q,r) * (dot_product(derx,uy_el) + dot_product(dery,ux_el))
201
202
203 enddo
204 enddo
205 enddo
206
207
208 return
209
210 end subroutine make_strain_rot_stress_tensor
make_strain_rot_stress_tensor
subroutine make_strain_rot_stress_tensor(nn, ct, ww, dd, alfa11, alfa12, alfa13, alfa21, alfa
Computes nodal values for ratiotional, strain and stress tensor.
Definition MAKE_STRAIN_ROT_STRESS_TENSOR.f90:75
phi
real *8 function phi(ng, csii, csij)
Evaluates Lengendre basis function on a specific point.
Definition PHI.f90:31