SPEED
Functions/Subroutines
MAKE_INTERNAL_FORCE.f90 File Reference

Go to the source code of this file.

Functions/Subroutines

subroutine make_internal_force (nn, xq, wq, dd, aa11, aa12, aa13, aa21, aa22, aa23, aa31, aa32, aa33, bb11, bb12, bb13, bb21, bb22, bb23, bb31, bb32, bb33, gg1, gg2, gg3, dd1, dd2, dd3, sxx, syy, szz, syz, szx, sxy, fx, fy, fz)
 Makes internal forces.
 

Function/Subroutine Documentation

◆ make_internal_force()

subroutine make_internal_force ( integer*4  nn,
real*8, dimension(nn)  xq,
real*8, dimension(nn)  wq,
real*8, dimension(nn,nn)  dd,
real*8  aa11,
real*8  aa12,
real*8  aa13,
real*8  aa21,
real*8  aa22,
real*8  aa23,
real*8  aa31,
real*8  aa32,
real*8  aa33,
real*8  bb11,
real*8  bb12,
real*8  bb13,
real*8  bb21,
real*8  bb22,
real*8  bb23,
real*8  bb31,
real*8  bb32,
real*8  bb33,
real*8  gg1,
real*8  gg2,
real*8  gg3,
real*8  dd1,
real*8  dd2,
real*8  dd3,
real*8, dimension(nn,nn,nn)  sxx,
real*8, dimension(nn,nn,nn)  syy,
real*8, dimension(nn,nn,nn)  szz,
real*8, dimension(nn,nn,nn)  syz,
real*8, dimension(nn,nn,nn)  szx,
real*8, dimension(nn,nn,nn)  sxy,
real*8, dimension(nn,nn,nn)  fx,
real*8, dimension(nn,nn,nn)  fy,
real*8, dimension(nn,nn,nn)  fz 
)

Makes internal forces.

Author
Ilario Mazzieri
Date
September, 2013
Version
1.0
Parameters
[in]nnnumber of 1-D GLL nodes
[in]xqGLL nodes
[in]wqGLL weights
[in]ddspectral derivative matrix
[in]A11costant values for the bilinear map
[in]A12costant values for the bilinear map
[in]A13costant values for the bilinear map
[in]A21costant values for the bilinear map
[in]A22costant values for the bilinear map
[in]A23costant values for the bilinear map
[in]A31costant values for the bilinear map
[in]A32costant values for the bilinear map
[in]A33costant values for the bilinear map
[in]B11costant values for the bilinear map
[in]B12costant values for the bilinear map
[in]B13costant values for the bilinear map
[in]B21costant values for the bilinear map
[in]B22costant values for the bilinear map
[in]B23costant values for the bilinear map
[in]B31costant values for the bilinear map
[in]B32costant values for the bilinear map
[in]B33costant values for the bilinear map
[in]GG1costant values for the bilinear map
[in]GG2costant values for the bilinear map
[in]GG3costant values for the bilinear map
[in]DD1costant values for the bilinear map
[in]DD2costant values for the bilinear map
[in]DD3costant values for the bilinear map
[in]sxxnodal values for the stress tensor
[in]syynodal values for the stress tensor
[in]szznodal values for the stress tensor
[in]syznodal values for the stress tensor
[in]szxnodal values for the stress tensor
[in]sxynodal values for the stress tensor
[out]fxx-componnent for internal forces
[out]fyy-componnent for internal forces
[out]fzz-componnent for internal forces

Definition at line 61 of file MAKE_INTERNAL_FORCE.f90.

67
68 implicit none
69
70 integer*4 :: nn
71 integer*4 :: i,j,k,p,q,r
72
73 real*8 :: aa11,aa12,aa13,aa21,aa22,aa23,aa31,aa32,aa33
74 real*8 :: bb11,bb12,bb13,bb21,bb22,bb23,bb31,bb32,bb33
75 real*8 :: gg1,gg2,gg3,dd1,dd2,dd3
76 real*8 :: dxdx,dxdy,dxdz,dydx,dydy,dydz,dzdx,dzdy,dzdz,det_j
77 real*8 :: duxdx,duxdy,duxdz,duydx,duydy,duydz,duzdx,duzdy,duzdz
78 real*8 :: t1ux,t2ux,t3ux,t1uy,t2uy,t3uy,t1uz,t2uz,t3uz
79
80 real*8, dimension(nn) :: xq,wq
81
82 real*8, dimension(nn,nn) :: dd
83
84 real*8, dimension(nn,nn,nn) :: sxx,syy,szz,syz,szx,sxy
85 real*8, dimension(nn,nn,nn) :: fx,fy,fz
86
87! FORCE CALCULATION
88
89 do r = 1,nn
90 do q = 1,nn
91 do p = 1,nn
92 t1ux = 0.d0; t1uy = 0.d0; t1uz = 0.d0
93 t2ux = 0.d0; t2uy = 0.d0; t2uz = 0.d0
94 t3ux = 0.d0; t3uy = 0.d0; t3uz = 0.d0
95
96 do i = 1,nn
97 dxdy = aa12 + bb11*xq(r) + bb13*xq(i) + gg1*xq(r)*xq(i)
98 dydy = aa22 + bb21*xq(r) + bb23*xq(i) + gg2*xq(r)*xq(i)
99 dzdy = aa32 + bb31*xq(r) + bb33*xq(i) + gg3*xq(r)*xq(i)
100
101 dxdz = aa13 + bb11*xq(q) + bb12*xq(i) + gg1*xq(i)*xq(q)
102 dydz = aa23 + bb21*xq(q) + bb22*xq(i) + gg2*xq(i)*xq(q)
103 dzdz = aa33 + bb31*xq(q) + bb32*xq(i) + gg3*xq(i)*xq(q)
104
105 t1ux = t1ux + wq(i)*wq(q)*wq(r) * dd(i,p) * &
106 ((dydy*dzdz - dydz*dzdy) * sxx(i,q,r) + &
107 (dzdy*dxdz - dzdz*dxdy) * sxy(i,q,r) + &
108 (dxdy*dydz - dxdz*dydy) * szx(i,q,r))
109 t1uy = t1uy + wq(i)*wq(q)*wq(r) * dd(i,p) * &
110 ((dydy*dzdz - dydz*dzdy) * sxy(i,q,r) + &
111 (dzdy*dxdz - dzdz*dxdy) * syy(i,q,r) + &
112 (dxdy*dydz - dxdz*dydy) * syz(i,q,r))
113 t1uz = t1uz + wq(i)*wq(q)*wq(r) * dd(i,p) * &
114 ((dydy*dzdz - dydz*dzdy) * szx(i,q,r) + &
115 (dzdy*dxdz - dzdz*dxdy) * syz(i,q,r) + &
116 (dxdy*dydz - dxdz*dydy) * szz(i,q,r))
117 enddo
118
119 do j = 1,nn
120 dxdx = aa11 + bb12*xq(r) + bb13*xq(j) + gg1*xq(j)*xq(r)
121 dydx = aa21 + bb22*xq(r) + bb23*xq(j) + gg2*xq(j)*xq(r)
122 dzdx = aa31 + bb32*xq(r) + bb33*xq(j) + gg3*xq(j)*xq(r)
123
124 dxdz = aa13 + bb11*xq(j) + bb12*xq(p) + gg1*xq(p)*xq(j)
125 dydz = aa23 + bb21*xq(j) + bb22*xq(p) + gg2*xq(p)*xq(j)
126 dzdz = aa33 + bb31*xq(j) + bb32*xq(p) + gg3*xq(p)*xq(j)
127
128 t2ux = t2ux + wq(p)*wq(j)*wq(r) * dd(j,q) * &
129 ((dydz*dzdx - dydx*dzdz) * sxx(p,j,r) + &
130 (dzdz*dxdx - dzdx*dxdz) * sxy(p,j,r) + &
131 (dxdz*dydx - dxdx*dydz) * szx(p,j,r))
132 t2uy = t2uy + wq(p)*wq(j)*wq(r) * dd(j,q) * &
133 ((dydz*dzdx - dydx*dzdz) * sxy(p,j,r) + &
134 (dzdz*dxdx - dzdx*dxdz) * syy(p,j,r) + &
135 (dxdz*dydx - dxdx*dydz) * syz(p,j,r))
136 t2uz = t2uz + wq(p)*wq(j)*wq(r) * dd(j,q) * &
137 ((dydz*dzdx - dydx*dzdz) * szx(p,j,r) + &
138 (dzdz*dxdx - dzdx*dxdz) * syz(p,j,r) + &
139 (dxdz*dydx - dxdx*dydz) * szz(p,j,r))
140 enddo
141
142 do k = 1,nn
143 dxdx = aa11 + bb12*xq(k) + bb13*xq(q) + gg1*xq(q)*xq(k)
144 dydx = aa21 + bb22*xq(k) + bb23*xq(q) + gg2*xq(q)*xq(k)
145 dzdx = aa31 + bb32*xq(k) + bb33*xq(q) + gg3*xq(q)*xq(k)
146
147 dxdy = aa12 + bb11*xq(k) + bb13*xq(p) + gg1*xq(k)*xq(p)
148 dydy = aa22 + bb21*xq(k) + bb23*xq(p) + gg2*xq(k)*xq(p)
149 dzdy = aa32 + bb31*xq(k) + bb33*xq(p) + gg3*xq(k)*xq(p)
150
151 t3ux = t3ux + wq(p)*wq(q)*wq(k) * dd(k,r) * &
152 ((dydx*dzdy - dydy*dzdx) * sxx(p,q,k) + &
153 (dzdx*dxdy - dzdy*dxdx) * sxy(p,q,k) + &
154 (dxdx*dydy - dxdy*dydx) * szx(p,q,k))
155 t3uy = t3uy + wq(p)*wq(q)*wq(k) * dd(k,r) * &
156 ((dydx*dzdy - dydy*dzdx) * sxy(p,q,k) + &
157 (dzdx*dxdy - dzdy*dxdx) * syy(p,q,k) + &
158 (dxdx*dydy - dxdy*dydx) * syz(p,q,k))
159 t3uz = t3uz + wq(p)*wq(q)*wq(k) * dd(k,r) * &
160 ((dydx*dzdy - dydy*dzdx) * szx(p,q,k) + &
161 (dzdx*dxdy - dzdy*dxdx) * syz(p,q,k) + &
162 (dxdx*dydy - dxdy*dydx) * szz(p,q,k))
163 enddo
164
165 fx(p,q,r) = t1ux + t2ux + t3ux
166 fy(p,q,r) = t1uy + t2uy + t3uy
167 fz(p,q,r) = t1uz + t2uz + t3uz
168 enddo
169 enddo
170 enddo
171
172
173 return
174