SPEED
NEWTON_RAPSON.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
61
62 subroutine newton_rapson(xnod, ynod, znod, &
63 alfa11,alfa12,alfa13, &
64 alfa21,alfa22,alfa23, &
65 alfa31,alfa32,alfa33, &
66 beta11,beta12,beta13, &
67 beta21,beta22,beta23, &
68 beta31,beta32,beta33, &
69 gamma1,gamma2,gamma3, &
70 delta1,delta2,delta3, &
71 tt, csi_s, eta_s, zeta_s, nofi, id_mpi,&
72 ipiu, imeno, toll1, toll2 , flag)
73
74 implicit none
75
76 integer*4 :: nofiter
77 integer*4 :: per_csi, per_eta, per_zeta
78 integer*4 :: iic, il,ih
79 integer*4, intent(out) :: tt
80 integer*4, intent(in) :: nofi,id_mpi,ipiu,imeno, flag
81
82 real*8 :: ll,lm ,xvt1,yvt1,xvt2,yvt2,xvc1,yvc1,xvc2,yvc2
83 real*8 :: a1,a2,a3,b1,b2,b3,c1,c2,c3,d1,d2,d3
84 real*8 :: e1,e2,e3,f1,f2,f3,g1,g2,g3,h1,h2,h3
85 real*8 :: delta_csi, delta_eta, delta_zeta
86 real*8 :: a,b,c,d,e,f,g,h,p, det
87 real*8 :: a11, a12, a13, a21, a22, a23, a31, a32, a33
88 real*8 :: ff1, ff2, ff3
89 real*8 :: alfa11,alfa12,alfa13,alfa21,alfa22,alfa23,alfa31,alfa32,alfa33
90 real*8 :: beta11,beta12,beta13,beta21,beta22,beta23,beta31,beta32,beta33
91 real*8 :: gamma1,gamma2,gamma3,delta1,delta2,delta3
92 real*8, intent(in) :: xnod, ynod, znod, toll1, toll2
93 real*8, intent(out) :: csi_s, eta_s, zeta_s
94
95
96 a1 = 8.d0*gamma1
97 a2 = 8.d0*gamma2
98 a3 = 8.d0*gamma3
99 b1 = 8.d0*beta13
100 b2 = 8.d0*beta23
101 b3 = 8.d0*beta33
102 c1 = 8.d0*beta11
103 c2 = 8.d0*beta21
104 c3 = 8.d0*beta31
105 d1 = 8.d0*beta12
106 d2 = 8.d0*beta22
107 d3 = 8.d0*beta32
108 e1 = 8.d0*alfa11
109 e2 = 8.d0*alfa21
110 e3 = 8.d0*alfa31
111 f1 = 8.d0*alfa12
112 f2 = 8.d0*alfa22
113 f3 = 8.d0*alfa32
114 g1 = 8.d0*alfa13
115 g2 = 8.d0*alfa23
116 g3 = 8.d0*alfa33
117 h1 = 8.d0*delta1
118 h2 = 8.d0*delta2
119 h3 = 8.d0*delta3
120
121 csi_s = 0.d0
122 eta_s = 0.d0
123 zeta_s = 0.d0
124 delta_csi = 0.0d0
125 delta_eta = 0.0d0
126 delta_zeta = 0.0d0
127
128 tt = 0
129 nofiter = 2000
130
131 per_csi = 0
132 per_eta = 0
133 per_zeta = 0
134
135
136 do il = 1, nofiter
137
138 a = a1*eta_s*zeta_s + b1*eta_s + d1*zeta_s + e1
139 b = a1*csi_s*zeta_s + b1*csi_s + c1*zeta_s + f1
140 c = a1*csi_s*eta_s + c1*eta_s + d1*csi_s + g1
141
142 d = a2*eta_s*zeta_s + b2*eta_s + d2*zeta_s + e2
143 e = a2*csi_s*zeta_s + b2*csi_s + c2*zeta_s + f2
144 f = a2*csi_s*eta_s + c2*eta_s + d2*csi_s + g2
145
146 g = a3*eta_s*zeta_s + b3*eta_s + d3*zeta_s + e3
147 h = a3*csi_s*zeta_s + b3*csi_s + c3*zeta_s + f3
148 p = a3*csi_s*eta_s + c3*eta_s + d3*csi_s + g3
149
150 det = a*(e*p - f*h) - b*(d*p - f*g) + c*(d*h -e*g)
151
152 if(det .eq. 0.d0) then
153 return
154 endif
155
156 a11 = e*p - f*h
157 a12 = -d*p + f*g
158 a13 = d*h - e*g
159 a21 = -b*p + c*h
160 a22 = a*p - c*g
161 a23 = -a*h + b*g
162 a31 = b*f - c*e
163 a32 = -a*f + c*d
164 a33 = a*e - b*d
165
166 ff1 = a1*csi_s*eta_s*zeta_s + b1*csi_s*eta_s + c1*eta_s*zeta_s + d1*csi_s*zeta_s &
167 + e1*csi_s + f1*eta_s + g1*zeta_s + h1 -8.d0*xnod
168
169 ff2 = a2*csi_s*eta_s*zeta_s + b2*csi_s*eta_s + c2*eta_s*zeta_s + d2*csi_s*zeta_s &
170 + e2*csi_s + f2*eta_s + g2*zeta_s + h2 -8.d0*ynod
171
172 ff3 = a3*csi_s*eta_s*zeta_s + b3*csi_s*eta_s + c3*eta_s*zeta_s + d3*csi_s*zeta_s &
173 + e3*csi_s + f3*eta_s + g3*zeta_s + h3 -8.d0*znod
174
175
176 delta_csi = (-1.d0/det)*(a11*ff1 + a21*ff2 + a31*ff3)
177
178 delta_eta = (-1.d0/det)*(a12*ff1 + a22*ff2 + a32*ff3)
179
180 delta_zeta = (-1.d0/det)*(a13*ff1 + a23*ff2 + a33*ff3)
181
182
183
184 if (abs(delta_csi).le. toll1 .and. abs(delta_eta).le. toll1 .and. abs(delta_zeta).le. toll1 ) then
185 tt = 1
186
187 if(flag.eq.1) then
188 if (dabs(csi_s) .gt. toll2 .or. dabs(eta_s) .gt. toll2 .or. dabs(zeta_s).gt. toll2 ) then
189 tt = 0
190 return
191 endif
192 endif
193
194 return
195 endif
196
197
198 csi_s = csi_s + delta_csi
199 eta_s = eta_s + delta_eta
200 zeta_s = zeta_s + delta_zeta
201
202 if (csi_s .gt. 1.d0 .AND. csi_s .lt. 2.d0 .AND. per_csi .eq. 0) then
203 csi_s = 1.d0
204 per_csi = 1
205 endif
206
207 if (csi_s .lt. -1.d0 .AND. csi_s .gt. -2.d0 .AND. per_csi .eq. 0) then
208 csi_s = -1.d0
209 per_csi = 1
210 endif
211
212
213 if (eta_s .gt. 1.d0 .AND. eta_s .lt. 2.d0 .AND. per_eta .eq. 0) then
214 eta_s = 1.d0
215 per_eta = 1
216 endif
217
218 if (eta_s .lt. -1.d0 .AND. eta_s .gt. -2.d0 .AND. per_eta .eq. 0) then
219 eta_s = -1.d0
220 per_eta = 1
221 endif
222
223 if (zeta_s .gt. 1.d0 .AND. zeta_s .lt. 2.d0 .AND. per_zeta .eq. 0) then
224 zeta_s = 1.d0
225 per_zeta = 1
226 endif
227
228 if (zeta_s .lt. -1.d0 .AND. zeta_s .gt. -2.d0 .AND. per_zeta .eq. 0) then
229 zeta_s = -1.d0
230 per_zeta = 1
231 endif
232
233
234 enddo
235
236 if(il .ge. nofiter) then
237 tt = 0
238 return
239 endif
240
241 end subroutine newton_rapson
newton_rapson
subroutine newton_rapson(xnod, ynod, znod, alfa11, alfa12, alfa13, alfa21, alfa22, alfa23, alfa31, alfa32, alfa33, beta11, beta12, beta13, beta21, beta22, beta23, beta31, beta32, beta33, gamma1, gamma2, gamma3, delta1, delta2, delta3, tt, csi_s, eta_s, zeta_s, nofi, id_mpi, ipiu, imeno, toll1, toll2, flag)
Performs the Newton Rapson algorithm.
Definition NEWTON_RAPSON.f90:73