Makes local spectral grid.
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
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))
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))
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)
174 alfa21(ie) = 0.125d0*(-y1 +y2 +y3 -y4 -y5 +y6 +y7 -y8)
175 alfa31(ie) = 0.125d0*(-z1 +z2 +z3 -z4 -z5 +z6 +z7 -z8)
176
177 alfa12(ie) = 0.125d0*(-x1 -x2 +x3 +x4 -x5 -x6 +x7 +x8)
178 alfa22(ie) = 0.125d0*(-y1 -y2 +y3 +y4 -y5 -y6 +y7 +y8)
179 alfa32(ie) = 0.125d0*(-z1 -z2 +z3 +z4 -z5 -z6 +z7 +z8)
180
181 alfa13(ie) = 0.125d0*(-x1 -x2 -x3 -x4 +x5 +x6 +x7 +x8)
182 alfa23(ie) = 0.125d0*(-y1 -y2 -y3 -y4 +y5 +y6 +y7 +y8)
183 alfa33(ie) = 0.125d0*(-z1 -z2 -z3 -z4 +z5 +z6 +z7 +z8)
184
185 beta11(ie) = 0.125d0*(+x1 +x2 -x3 -x4 -x5 -x6 +x7 +x8)
186 beta21(ie) = 0.125d0*(+y1 +y2 -y3 -y4 -y5 -y6 +y7 +y8)
187 beta31(ie) = 0.125d0*(+z1 +z2 -z3 -z4 -z5 -z6 +z7 +z8)
188
189 beta12(ie) = 0.125d0*(+x1 -x2 -x3 +x4 -x5 +x6 +x7 -x8)
190 beta22(ie) = 0.125d0*(+y1 -y2 -y3 +y4 -y5 +y6 +y7 -y8)
191 beta32(ie) = 0.125d0*(+z1 -z2 -z3 +z4 -z5 +z6 +z7 -z8)
192
193 beta13(ie) = 0.125d0*(+x1 -x2 +x3 -x4 +x5 -x6 +x7 -x8)
194 beta23(ie) = 0.125d0*(+y1 -y2 +y3 -y4 +y5 -y6 +y7 -y8)
195 beta33(ie) = 0.125d0*(+z1 -z2 +z3 -z4 +z5 -z6 +z7 -z8)
196
197 gamma1(ie) = 0.125d0*(-x1 +x2 -x3 +x4 +x5 -x6 +x7 -x8)
198 gamma2(ie) = 0.125d0*(-y1 +y2 -y3 +y4 +y5 -y6 +y7 -y8)
199 gamma3(ie) = 0.125d0*(-z1 +z2 -z3 +z4 +z5 -z6 +z7 -z8)
200
201 delta1(ie) = 0.125d0*(+x1 +x2 +x3 +x4 +x5 +x6 +x7 +x8)
202 delta2(ie) = 0.125d0*(+y1 +y2 +y3 +y4 +y5 +y6 +y7 +y8)
203 delta3(ie) = 0.125d0*(+z1 +z2 +z3 +z4 +z5 +z6 +z7 +z8)
204
205
206 jac = alfa31(ie)*(alfa12(ie)*alfa23(ie) -alfa13(ie)*alfa22
207 - alfa32(ie)*(alfa11(ie)*alfa23(ie) -alfa13(ie)*alfa21
208 + alfa33(ie)*(alfa11(ie)*alfa22(ie) -alfa12(ie)*alfa21
209
210
211
212
213 if (jac.le.0.d0) then
214 write(*,*)'Error ! Orientation non-conforming !'
215 write(*,*)'(element ',ie,' is clockwise oriented)'
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
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
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
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
subroutine make_lgl_nw(nb_pnt, xq, wq, dd)
Makes Gauss-Legendre-Lobatto nodes, weigths and spectral derivatives.
integer, parameter exit_elem_orient