SPEED
MAKE_SPX_GRID_LOC.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
58
59
60 subroutine make_spx_grid_loc(nn_loc, xx_loc, yy_loc, zz_loc, &
61 cs_nnz_loc, cs_loc, nm, tm, sd, ne_loc, loc_n_num, &
62 alfa11,alfa12,alfa13, &
63 alfa21,alfa22,alfa23, &
64 alfa31,alfa32,alfa33, &
65 beta11,beta12,beta13, &
66 beta21,beta22,beta23, &
67 beta31,beta32,beta33, &
68 gamma1,gamma2,gamma3, &
69 delta1,delta2,delta3 )
70
71
72 use speed_exit_codes
73
74 implicit none
75
76 integer*4 :: nn_loc, cs_nnz_loc,nm, ne_loc
77 integer*4 :: n1,n2,n3,n4,n5,n6,n7,n8
78 integer*4 :: ic1,ic2,ic3,ic4,ic5,ic6,ic7,ic8
79 integer*4 :: im,ie,i,j,k,nn,ip, ic
80
81 integer*4, dimension(nn_loc) :: loc_n_num
82 integer*4, dimension(0:cs_nnz_loc) :: cs_loc
83 integer*4, dimension(nm) :: tm
84 integer*4, dimension(nm) :: sd
85
86 real*8 :: xp,yp,zp,jac
87 real*8 :: x1,x2,x3,x4,x5,x6,x7,x8
88 real*8 :: y1,y2,y3,y4,y5,y6,y7,y8
89 real*8 :: z1,z2,z3,z4,z5,z6,z7,z8
90
91 real*8, dimension(:), allocatable :: ct,ww
92 real*8, dimension(ne_loc) :: alfa11,alfa12,alfa13
93 real*8, dimension(ne_loc) :: alfa21,alfa22,alfa23
94 real*8, dimension(ne_loc) :: alfa31,alfa32,alfa33
95 real*8, dimension(ne_loc) :: beta11,beta12,beta13
96 real*8, dimension(ne_loc) :: beta21,beta22,beta23
97 real*8, dimension(ne_loc) :: beta31,beta32,beta33
98 real*8, dimension(ne_loc) :: gamma1,gamma2,gamma3
99 real*8, dimension(ne_loc) :: delta1,delta2,delta3
100 real*8, dimension(nn_loc), intent(inout) :: xx_loc,yy_loc,zz_loc
101
102 real*8, dimension(:,:), allocatable :: dd
103
104 nn = 2
105 allocate(ct(nn),ww(nn),dd(nn,nn))
106 call make_lgl_nw(nn,ct,ww,dd)
107
108 do im = 1,nm
109 if ((sd(im) +1).ne.nn) then
110 deallocate(ct,ww,dd)
111 nn = sd(im) +1
112 allocate(ct(nn),ww(nn),dd(nn,nn))
113 call make_lgl_nw(nn,ct,ww,dd)
114 endif
115
116 do ie = 1,ne_loc
117 if (cs_loc(cs_loc(ie -1) +0).eq.tm(im)) then
118 n1 = nn*nn*(1 -1) +nn*(1 -1) +1
119 n2 = nn*nn*(1 -1) +nn*(1 -1) +nn
120 n3 = nn*nn*(1 -1) +nn*(nn -1) +nn
121 n4 = nn*nn*(1 -1) +nn*(nn -1) +1
122 n5 = nn*nn*(nn -1) +nn*(1 -1) +1
123 n6 = nn*nn*(nn -1) +nn*(1 -1) +nn
124 n7 = nn*nn*(nn -1) +nn*(nn -1) +nn
125 n8 = nn*nn*(nn -1) +nn*(nn -1) +1
126
127
128 ic1 = cs_loc(cs_loc(ie -1) +n1)
129 ic2 = cs_loc(cs_loc(ie -1) +n2)
130 ic3 = cs_loc(cs_loc(ie -1) +n3)
131 ic4 = cs_loc(cs_loc(ie -1) +n4)
132 ic5 = cs_loc(cs_loc(ie -1) +n5)
133 ic6 = cs_loc(cs_loc(ie -1) +n6)
134 ic7 = cs_loc(cs_loc(ie -1) +n7)
135 ic8 = cs_loc(cs_loc(ie -1) +n8)
136
137
138 x1 = xx_loc(ic1)
139 y1 = yy_loc(ic1)
140 z1 = zz_loc(ic1)
141
142 x2 = xx_loc(ic2)
143 y2 = yy_loc(ic2)
144 z2 = zz_loc(ic2)
145
146 x3 = xx_loc(ic3)
147 y3 = yy_loc(ic3)
148 z3 = zz_loc(ic3)
149
150 x4 = xx_loc(ic4)
151 y4 = yy_loc(ic4)
152 z4 = zz_loc(ic4)
153
154 x5 = xx_loc(ic5)
155 y5 = yy_loc(ic5)
156 z5 = zz_loc(ic5)
157
158 x6 = xx_loc(ic6)
159 y6 = yy_loc(ic6)
160 z6 = zz_loc(ic6)
161
162 x7 = xx_loc(ic7)
163 y7 = yy_loc(ic7)
164 z7 = zz_loc(ic7)
165
166 x8 = xx_loc(ic8)
167 y8 = yy_loc(ic8)
168 z8 = zz_loc(ic8)
169
170
171
172
173 alfa11(ie) = 0.125d0*(-x1 +x2 +x3 -x4 -x5 +x6 +x7 -x8) !1/8 * e1
174 alfa21(ie) = 0.125d0*(-y1 +y2 +y3 -y4 -y5 +y6 +y7 -y8) !1/8 * e2
175 alfa31(ie) = 0.125d0*(-z1 +z2 +z3 -z4 -z5 +z6 +z7 -z8) !1/8 * e3
176
177 alfa12(ie) = 0.125d0*(-x1 -x2 +x3 +x4 -x5 -x6 +x7 +x8) !1/8 * f1
178 alfa22(ie) = 0.125d0*(-y1 -y2 +y3 +y4 -y5 -y6 +y7 +y8) !1/8 * f2
179 alfa32(ie) = 0.125d0*(-z1 -z2 +z3 +z4 -z5 -z6 +z7 +z8) !1/8 * f3
180
181 alfa13(ie) = 0.125d0*(-x1 -x2 -x3 -x4 +x5 +x6 +x7 +x8) !1/8 * g1
182 alfa23(ie) = 0.125d0*(-y1 -y2 -y3 -y4 +y5 +y6 +y7 +y8) !1/8 * g2
183 alfa33(ie) = 0.125d0*(-z1 -z2 -z3 -z4 +z5 +z6 +z7 +z8) !1/8 * g3
184
185 beta11(ie) = 0.125d0*(+x1 +x2 -x3 -x4 -x5 -x6 +x7 +x8) !1/8 * c1
186 beta21(ie) = 0.125d0*(+y1 +y2 -y3 -y4 -y5 -y6 +y7 +y8) !1/8 * c2
187 beta31(ie) = 0.125d0*(+z1 +z2 -z3 -z4 -z5 -z6 +z7 +z8) !1/8 * c3
188
189 beta12(ie) = 0.125d0*(+x1 -x2 -x3 +x4 -x5 +x6 +x7 -x8) !1/8 * d1
190 beta22(ie) = 0.125d0*(+y1 -y2 -y3 +y4 -y5 +y6 +y7 -y8) !1/8 * d2
191 beta32(ie) = 0.125d0*(+z1 -z2 -z3 +z4 -z5 +z6 +z7 -z8) !1/8 * d3
192
193 beta13(ie) = 0.125d0*(+x1 -x2 +x3 -x4 +x5 -x6 +x7 -x8) !1/8 * b1
194 beta23(ie) = 0.125d0*(+y1 -y2 +y3 -y4 +y5 -y6 +y7 -y8) !1/8 * b2
195 beta33(ie) = 0.125d0*(+z1 -z2 +z3 -z4 +z5 -z6 +z7 -z8) !1/8 * b3
196
197 gamma1(ie) = 0.125d0*(-x1 +x2 -x3 +x4 +x5 -x6 +x7 -x8) !1/8 * a1
198 gamma2(ie) = 0.125d0*(-y1 +y2 -y3 +y4 +y5 -y6 +y7 -y8) !1/8 * a2
199 gamma3(ie) = 0.125d0*(-z1 +z2 -z3 +z4 +z5 -z6 +z7 -z8) !1/8 * a3
200
201 delta1(ie) = 0.125d0*(+x1 +x2 +x3 +x4 +x5 +x6 +x7 +x8) !1/8 * h1
202 delta2(ie) = 0.125d0*(+y1 +y2 +y3 +y4 +y5 +y6 +y7 +y8) !1/8 * h2
203 delta3(ie) = 0.125d0*(+z1 +z2 +z3 +z4 +z5 +z6 +z7 +z8) !1/8 * h3
204
205
206 jac = alfa31(ie)*(alfa12(ie)*alfa23(ie) -alfa13(ie)*alfa22(ie)) &
207 - alfa32(ie)*(alfa11(ie)*alfa23(ie) -alfa13(ie)*alfa21(ie)) &
208 + alfa33(ie)*(alfa11(ie)*alfa22(ie) -alfa12(ie)*alfa21(ie))
209
210
211
212
213 if (jac.le.0.d0) then
214 write(*,*)'Error ! Orientation non-conforming !'
215 write(*,*)'(element ',ie,' is clockwise oriented)'
216 call exit(exit_elem_orient)
217 endif
218
219 do k = 1,nn
220 do j = 1,nn
221 do i = 1,nn
222 xp = alfa11(ie)*ct(i) + alfa12(ie)*ct(j) &
223 + alfa13(ie)*ct(k) + beta11(ie)*ct(j)*ct(k) &
224 + beta12(ie)*ct(i)*ct(k) + beta13(ie)*ct(i)*ct(j) &
225 + gamma1(ie)*ct(i)*ct(j)*ct(k) + delta1(ie)
226
227 yp = alfa21(ie)*ct(i) + alfa22(ie)*ct(j) &
228 + alfa23(ie)*ct(k) + beta21(ie)*ct(j)*ct(k) &
229 + beta22(ie)*ct(i)*ct(k) + beta23(ie)*ct(i)*ct(j) &
230 + gamma2(ie)*ct(i)*ct(j)*ct(k) + delta2(ie)
231
232 zp = alfa31(ie)*ct(i) + alfa32(ie)*ct(j) &
233 + alfa33(ie)*ct(k) + beta31(ie)*ct(j)*ct(k) &
234 + beta32(ie)*ct(i)*ct(k) + beta33(ie)*ct(i)*ct(j) &
235 + gamma3(ie)*ct(i)*ct(j)*ct(k) + delta3(ie)
236
237 ip = nn*nn*(k -1) +nn*(j -1) +i
238
239 xx_loc(cs_loc(cs_loc(ie -1) + ip)) = xp
240 yy_loc(cs_loc(cs_loc(ie -1) + ip)) = yp
241 zz_loc(cs_loc(cs_loc(ie -1) + ip)) = zp
242
243 enddo
244 enddo
245 enddo
246
247 endif
248 enddo
249 enddo
250
251 return
252
253 end subroutine make_spx_grid_loc
254
make_lgl_nw
subroutine make_lgl_nw(nb_pnt, xq, wq, dd)
Makes Gauss-Legendre-Lobatto nodes, weigths and spectral derivatives.
Definition MAKE_LGL_NW.f90:29
make_spx_grid_loc
subroutine make_spx_grid_loc(nn_loc, xx_loc, yy_loc, zz_loc, cs_nnz_loc, cs_loc, nm, tm, sd, ne_loc, l
Makes local spectral grid.
Definition MAKE_SPX_GRID_LOC.f90:62
speed_exit_codes
SPEED exit codes.
Definition MODULES.f90:25
speed_exit_codes::exit_elem_orient
integer, parameter exit_elem_orient
Definition MODULES.f90:40