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

Go to the source code of this file.

Functions/Subroutines

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

Function/Subroutine Documentation

◆ extrapol_kel()

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