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

Go to the source code of this file.

Functions/Subroutines

subroutine extrapol_tel (nface, ielfa, xrlfa, vel, vfa, ifabou, xbounact, nef, gradtel, gradtfa, neifa, rf, area, volume, xle, xxi, icyclic, xxn, ipnei, ifatie, xload, xlet, xxj, coefmpc, nmpc, labmpc, ipompc, nodempc, ifaext, nfaext, nactdoh)
 

Function/Subroutine Documentation

◆ extrapol_tel()

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