CalculiX  2.13
A Free Software Three-Dimensional Structural Finite Element Program
Functions/Subroutines
extrapol_oel.f File Reference

Go to the source code of this file.

Functions/Subroutines

subroutine extrapol_oel (nface, ielfa, xrlfa, vel, vfa, ifabou, xbounact, nef, gradoel, gradofa, neifa, rf, area, volume, xle, xxi, icyclic, xxn, ipnei, ifatie, xlet, xxj, coefmpc, nmpc, labmpc, ipompc, nodempc, ifaext, nfaext, nactdoh, umfa, physcon, dy)
 

Function/Subroutine Documentation

◆ extrapol_oel()

subroutine extrapol_oel ( integer, intent(in)  nface,
integer, dimension(4,*), intent(in)  ielfa,
real*8, dimension(3,*), intent(in)  xrlfa,
real*8, dimension(nef,0:7), intent(in)  vel,
real*8, dimension(0:7,*), intent(inout)  vfa,
integer, dimension(*), intent(in)  ifabou,
real*8, dimension(*), intent(in)  xbounact,
integer, intent(in)  nef,
real*8, dimension(3,*), intent(inout)  gradoel,
real*8, dimension(3,*), intent(inout)  gradofa,
integer, dimension(*), intent(in)  neifa,
real*8, dimension(3), intent(in)  rf,
real*8, dimension(*), intent(in)  area,
real*8, dimension(*), intent(in)  volume,
real*8, dimension(*), intent(in)  xle,
real*8, dimension(3,*), intent(in)  xxi,
integer, intent(in)  icyclic,
real*8, dimension(3,*), intent(in)  xxn,
integer, dimension(*), intent(in)  ipnei,
integer, dimension(*), intent(in)  ifatie,
real*8, dimension(*), intent(in)  xlet,
real*8, dimension(3,*), intent(in)  xxj,
real*8, dimension(*), intent(in)  coefmpc,
integer, intent(in)  nmpc,
character*20, dimension(*), intent(in)  labmpc,
integer, dimension(*), intent(in)  ipompc,
integer, dimension(3,*), intent(in)  nodempc,
integer, dimension(*), intent(in)  ifaext,
integer, intent(in)  nfaext,
integer, dimension(*), intent(in)  nactdoh,
real*8, dimension(*), intent(in)  umfa,
real*8, dimension(*), intent(in)  physcon,
real*8, dimension(*)  dy 
)
24 !
25 ! extrapolation of temperature values to the faces
26 !
27  implicit none
28 !
29  character*20 labmpc(*)
30 !
31  integer nface,ielfa(4,*),ifabou(*),i,iel1,iel2,nef,
32  & neifa(*),icyclic,ifa,indexf,k,l,m,ipnei(*),ifatie(*),ipointer,
33  & is,ie,nmpc,ipompc(*),nodempc(3,*),ifaext(*),nfaext,nactdoh(*)
34 !
35  real*8 xrlfa(3,*),vel(nef,0:7),vfa(0:7,*),xbounact(*),xl1,xl2,
36  & vfap(0:7,nface),gradoel(3,*),gradofa(3,*),rf(3),area(*),
37  & volume(*),xle(*),xxi(3,*),c(3,3),gradnor,xxn(3,*),umfa(*),
38  & xxj(3,*),dd,xlet(*),coefmpc(*),cinf,physcon(*),dy(*)
39 !
40  intent(in) nface,ielfa,xrlfa,vel,umfa,physcon,
41  & ifabou,xbounact,nef,neifa,rf,area,volume,
42  & xle,xxi,icyclic,xxn,ipnei,ifatie,xlet,xxj,
43  & coefmpc,nmpc,labmpc,ipompc,nodempc,ifaext,nfaext,nactdoh
44 !
45  intent(inout) vfa,gradoel,gradofa
46 !
47 ! 5.5 x 10**(-3.5) = 1.7393e-3
48 !
49  cinf=5.5d0*physcon(5)/physcon(8)
50 !
51 !c$omp parallel default(none)
52 !c$omp& shared(nface,ielfa,xrlfa,vfap,vel,ifabou,xbounact,cinf,
53 !c$omp& umfa)
54 !c$omp& private(i,iel1,xl1,iel2,ibou)
55 !c$omp do
56  do i=1,nface
57  iel1=ielfa(1,i)
58  xl1=xrlfa(1,i)
59  iel2=ielfa(2,i)
60  if(iel2.gt.0) then
61 !
62 ! face between two elements: interpolation
63 !
64  vfap(7,i)=xl1*vel(iel1,7)+xrlfa(2,i)*vel(iel2,7)
65  elseif(ielfa(3,i).gt.0) then
66 !
67 ! boundary face; no zero gradient
68 !
69 ! iel2=0 is not possible: if iel2=0, there are no
70 ! boundary conditions on the face, hence it is an
71 ! exit, which means zero gradient and
72 ! ielfa(3,i) <= 0
73 !
74  ipointer=-iel2
75 !
76  if(ifabou(ipointer+5).gt.0) then
77 !
78 ! wall: turbulent dissipation rate known
79 !
80  vfap(7,i)=dy(ifabou(ipointer+5))*umfa(i)/vfa(5,i)
81  elseif(((ifabou(ipointer+1).gt.0).and.
82  & (ifabou(ipointer+2).gt.0).and.
83  & (ifabou(ipointer+3).gt.0)).or.
84  & (ifabou(ipointer+5).lt.0)) then
85 !
86 ! inlet or sliding conditions: turbulent dissipation rate known
87 !
88  vfap(7,i)=cinf*umfa(i)
89  else
90 !
91 ! extrapolation
92 !
93  vfap(7,i)=xl1*vel(iel1,7)+xrlfa(3,i)*vel(ielfa(3,i),7)
94  endif
95  else
96 !
97 ! boundary face; zero gradient
98 !
99  vfap(7,i)=vel(iel1,7)
100  endif
101  enddo
102 !c$omp end do
103 !c$omp end parallel
104 !
105 ! Multiple point constraints
106 !
107  if(nmpc.gt.0) then
108  is=6
109  ie=6
110  call applympc(nface,ielfa,is,ie,ifabou,ipompc,vfap,coefmpc,
111  & nodempc,ipnei,neifa,labmpc,xbounact,nactdoh,
112  & ifaext,nfaext)
113  endif
114 !
115 ! calculate the gradient of the temperature at the center of
116 ! the elements
117 !
118 !c$omp parallel default(none)
119 !c$omp& shared(nef,ipnei,neifa,gradoel,vfap,area,xxn,volume)
120 !c$omp& private(i,indexf,ifa)
121 !c$omp do
122  do i=1,nef
123 !
124 ! initialization
125 !
126  do l=1,3
127  gradoel(l,i)=0.d0
128  enddo
129 !
130  do indexf=ipnei(i)+1,ipnei(i+1)
131  ifa=neifa(indexf)
132  do l=1,3
133  gradoel(l,i)=gradoel(l,i)+
134  & vfap(7,ifa)*area(ifa)*xxn(l,indexf)
135  enddo
136  enddo
137 !
138 ! dividing by the volume of the element
139 !
140  do l=1,3
141  gradoel(l,i)=gradoel(l,i)/volume(i)
142  enddo
143  enddo
144 !c$omp end do
145 !c$omp end parallel
146 !
147 ! interpolate/extrapolate the temperature gradient from the
148 ! center of the elements to the center of the faces
149 !
150 !c$omp parallel default(none)
151 !c$omp& shared(nface,ielfa,xrlfa,gradofa,gradoel,icyclic,c,ifatie,
152 !c$omp& ipnei,xxn,ifabou)
153 !c$omp& private(i,iel1,xl1,iel2,l,xl2,indexf,gradnor)
154 !c$omp do
155  do i=1,nface
156  iel1=ielfa(1,i)
157  xl1=xrlfa(1,i)
158  iel2=ielfa(2,i)
159  if(iel2.gt.0) then
160 !
161 ! face between two elements
162 !
163  xl2=xrlfa(2,i)
164  if((icyclic.eq.0).or.(ifatie(i).eq.0)) then
165  do l=1,3
166  gradofa(l,i)=xl1*gradoel(l,iel1)+
167  & xl2*gradoel(l,iel2)
168  enddo
169  elseif(ifatie(i).gt.0) then
170  do l=1,3
171  gradofa(l,i)=xl1*gradoel(l,iel1)+xl2*
172  & (gradoel(1,iel2)*c(l,1)+
173  & gradoel(2,iel2)*c(l,2)+
174  & gradoel(3,iel2)*c(l,3))
175  enddo
176  else
177  do l=1,3
178  gradofa(l,i)=xl1*gradoel(l,iel1)+xl2*
179  & (gradoel(1,iel2)*c(1,l)+
180  & gradoel(2,iel2)*c(2,l)+
181  & gradoel(3,iel2)*c(3,l))
182  enddo
183  endif
184  elseif(ielfa(3,i).gt.0) then
185 !
186 ! the above condition implies iel2!=0
187 ! if iel2 were zero, no b.c. would apply
188 ! which means exit and hence zero gradient:
189 ! ielfa(3,i) would be zero or negative
190 !
191 ! boundary face; no zero gradient
192 !
193  do l=1,3
194  gradofa(l,i)=xl1*gradoel(l,iel1)+
195  & xrlfa(3,i)*gradoel(l,abs(ielfa(3,i)))
196  enddo
197  else
198 !
199 ! boundary face; zero gradient in i-direction
200 !
201  indexf=ipnei(iel1)+ielfa(4,i)
202  gradnor=gradoel(1,iel1)*xxi(1,indexf)+
203  & gradoel(2,iel1)*xxi(2,indexf)+
204  & gradoel(3,iel1)*xxi(3,indexf)
205  do l=1,3
206  gradofa(l,i)=gradoel(l,iel1)
207  & -gradnor*xxi(l,indexf)
208  enddo
209  endif
210  enddo
211 !c$omp end do
212 !c$omp end parallel
213 !
214 ! correction loops
215 !
216  do m=1,2
217 !
218 ! Moukalled et al. p 279
219 !
220 !c$omp parallel default(none)
221 !c$omp& shared(nface,ielfa,xrlfa,vfa,vfap,gradofa,rf,ifabou,xxi,xle,
222 !c$omp& vel,ipnei)
223 !c$omp& private(i,iel1,iel2,xl1,indexf)
224 !c$omp do
225  do i=1,nface
226  iel1=ielfa(1,i)
227  iel2=ielfa(2,i)
228  xl1=xrlfa(1,i)
229  if(iel2.gt.0) then
230 !
231 ! interpolation
232 !
233  vfa(7,i)=vfap(7,i)+gradofa(1,i)*rf(1)
234  & +gradofa(2,i)*rf(2)
235  & +gradofa(3,i)*rf(3)
236  elseif(ielfa(3,i).gt.0) then
237 !
238 ! no implicit zero gradient
239 !
240  ipointer=-iel2
241 !
242  if(ifabou(ipointer+5).gt.0) then
243 !
244 ! wall: turbulent dissipation rate known
245 !
246  vfa(7,i)=vfap(7,i)
247  elseif(((ifabou(ipointer+1).gt.0).and.
248  & (ifabou(ipointer+2).gt.0).and.
249  & (ifabou(ipointer+3).gt.0)).or.
250  & (ifabou(ipointer+5).lt.0)) then
251 !
252 ! inlet or sliding conditions: turbulent dissipation rate known
253 !
254  vfa(7,i)=vfap(7,i)
255  else
256 !
257 ! turbulent kinetic energy is not given
258 !
259  vfa(7,i)=vfap(7,i)+gradofa(1,i)*rf(1)+
260  & gradofa(2,i)*rf(2)+
261  & gradofa(3,i)*rf(3)
262  endif
263  else
264 !
265 ! zero gradient in i-direction
266 !
267  indexf=ipnei(iel1)+ielfa(4,i)
268  vfa(7,i)=vel(iel1,7)
269  & +(gradofa(1,i)*xxi(1,indexf)+
270  & gradofa(2,i)*xxi(2,indexf)+
271  & gradofa(3,i)*xxi(3,indexf))*xle(indexf)
272  endif
273  enddo
274 !c$omp end do
275 !c$omp end parallel
276 !
277 ! Multiple point constraints
278 !
279  if(nmpc.gt.0) then
280  is=6
281  ie=6
282  call applympc(nface,ielfa,is,ie,ifabou,ipompc,vfa,coefmpc,
283  & nodempc,ipnei,neifa,labmpc,xbounact,nactdoh,
284  & ifaext,nfaext)
285  endif
286 !
287 ! calculate the gradient of the temperature at the center of
288 ! the elements
289 !
290 !c$omp parallel default(none)
291 !c$omp& shared(nef,ipnei,neifa,gradoel,vfa,area,xxn,volume)
292 !c$omp& private(i,indexf,ifa)
293 !c$omp do
294  do i=1,nef
295 !
296 ! initialization
297 !
298  do l=1,3
299  gradoel(l,i)=0.d0
300  enddo
301 !
302  do indexf=ipnei(i)+1,ipnei(i+1)
303  ifa=neifa(indexf)
304  do l=1,3
305  gradoel(l,i)=gradoel(l,i)+
306  & vfa(7,ifa)*area(ifa)*xxn(l,indexf)
307  enddo
308  enddo
309 !
310 ! dividing by the volume of the element
311 !
312  do l=1,3
313  gradoel(l,i)=gradoel(l,i)/volume(i)
314  enddo
315  enddo
316 !c$omp end do
317 !c$omp end parallel
318 !
319 ! interpolate/extrapolate the temperature gradient from the
320 ! center of the elements to the center of the faces
321 !
322 !c$omp parallel default(none)
323 !c$omp& shared(nface,ielfa,xrlfa,gradofa,gradoel,icyclic,c,ifatie,
324 !c$omp& ipnei,xxn,ifabou)
325 !c$omp& private(i,iel1,xl1,iel2,l,xl2,indexf,gradnor)
326 !c$omp do
327  do i=1,nface
328  iel1=ielfa(1,i)
329  xl1=xrlfa(1,i)
330  iel2=ielfa(2,i)
331  if(iel2.gt.0) then
332 !
333 ! face in between two elements
334 !
335  xl2=xrlfa(2,i)
336  if((icyclic.eq.0).or.(ifatie(i).eq.0)) then
337  do l=1,3
338  gradofa(l,i)=xl1*gradoel(l,iel1)+
339  & xl2*gradoel(l,iel2)
340  enddo
341  elseif(ifatie(i).gt.0) then
342  do l=1,3
343  gradofa(l,i)=xl1*gradoel(l,iel1)+xl2*
344  & (gradoel(1,iel2)*c(l,1)+
345  & gradoel(2,iel2)*c(l,2)+
346  & gradoel(3,iel2)*c(l,3))
347  enddo
348  else
349  do l=1,3
350  gradofa(l,i)=xl1*gradoel(l,iel1)+xl2*
351  & (gradoel(1,iel2)*c(1,l)+
352  & gradoel(2,iel2)*c(2,l)+
353  & gradoel(3,iel2)*c(3,l))
354  enddo
355  endif
356  elseif(ielfa(3,i).gt.0) then
357 !
358 ! boundary face; no zero gradient
359 !
360  do l=1,3
361  gradofa(l,i)=xl1*gradoel(l,iel1)+
362  & xrlfa(3,i)*gradoel(l,abs(ielfa(3,i)))
363  enddo
364  else
365 !
366 ! boundary face; zero gradient in i-direction
367 !
368  indexf=ipnei(iel1)+ielfa(4,i)
369  gradnor=gradoel(1,iel1)*xxi(1,indexf)+
370  & gradoel(2,iel1)*xxi(2,indexf)+
371  & gradoel(3,iel1)*xxi(3,indexf)
372  do l=1,3
373  gradofa(l,i)=gradoel(l,iel1)
374  & -gradnor*xxi(l,indexf)
375  enddo
376  endif
377  enddo
378 !c$omp end do
379 !c$omp end parallel
380 !
381  enddo
382 !
383 ! correct the facial temperature gradients:
384 ! Moukalled et al. p 289
385 !
386 !c$omp parallel default(none)
387 !c$omp& shared(nface,ielfa,ipnei,vel,xlet,gradofa,xxj)
388 !c$omp& private(i,iel2,iel1,indexf,dd,k)
389 !c$omp do
390  do i=1,nface
391  iel2=ielfa(2,i)
392  if(iel2.gt.0) then
393  iel1=ielfa(1,i)
394  indexf=ipnei(iel1)+ielfa(4,i)
395  dd=(vel(iel2,7)-vel(iel1,7))/xlet(indexf)
396  & -gradofa(1,i)*xxj(1,indexf)
397  & -gradofa(2,i)*xxj(2,indexf)
398  & -gradofa(3,i)*xxj(3,indexf)
399  do k=1,3
400  gradofa(k,i)=gradofa(k,i)+dd*xxj(k,indexf)
401  enddo
402  endif
403  enddo
404 !c$omp end do
405 !c$omp end parallel
406 !
407  return
applympc
subroutine applympc(nface, ielfa, is, ie, ifabou, ipompc, vfa, coefmpc, nodempc, ipnei, neifa, labmpc, xbounact, nactdoh, ifaext, nfaext)
Definition: applympc.f:21
Hosted by OpenAircraft.com, (Michigan UAV, LLC)