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

Go to the source code of this file.

Functions/Subroutines

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

Function/Subroutine Documentation

◆ extrapol_vel()

subroutine extrapol_vel ( 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, dimension(*), intent(in)  ipnei,
integer, intent(in)  nef,
integer, intent(in)  icyclic,
real*8, dimension(3,3), intent(in)  c,
integer, dimension(*), intent(in)  ifatie,
real*8, dimension(3,*), intent(in)  xxn,
real*8, dimension(3,3,*), intent(inout)  gradvel,
real*8, dimension(3,3,*), intent(inout)  gradvfa,
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,
real*8, dimension(3,*), intent(in)  xxj,
real*8, dimension(*), intent(in)  xlet,
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 ! inter/extrapolation of v at the center of the elements
25 ! to the center of the faces
26 !
27  implicit none
28 !
29  character*20 labmpc(*)
30 !
31  integer nface,ielfa(4,*),ifabou(*),iel1,iel2,iel3,i,j,ipointer,
32  & indexf,ipnei(*),nef,icyclic,ifatie(*),k,l,m,ifa,neifa(*),
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  & c(3,3),xxn(3,*),dd,vfap(0:7,nface),gradvel(3,3,*),
37  & gradvfa(3,3,*),
38  & rf(3),area(*),volume(*),xle(*),xxi(3,*),gradnor,xxj(3,*),
39  & xlet(*),coefmpc(*)
40 !
41  intent(in) nface,ielfa,xrlfa,vel,
42  & ifabou,xbounact,ipnei,nef,icyclic,c,ifatie,xxn,
43  & neifa,rf,area,volume,xle,xxi,xxj,xlet,
44  & coefmpc,nmpc,labmpc,ipompc,nodempc,ifaext,nfaext,nactdoh
45 !
46  intent(inout) vfa,gradvel,gradvfa
47 !
48 ! initialization of the facial velocities
49 !
50 !c$omp parallel default(none)
51 !c$omp& shared(nface,ielfa,xrlfa,vfap,vel,ipnei,ifabou,xbounact,
52 !c$omp& icyclic,c,ifatie,xxn)
53 !c$omp& private(i,iel1,xl1,iel2,xl2,j,iel3,ipointer,indexf,dd)
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  xl2=xrlfa(2,i)
64  if((icyclic.eq.0).or.(ifatie(i).eq.0)) then
65  do j=1,3
66  vfap(j,i)=xl1*vel(iel1,j)+xl2*vel(iel2,j)
67  enddo
68  elseif(ifatie(i).gt.0) then
69  do j=1,3
70  vfap(j,i)=xl1*vel(iel1,j)+xl2*
71  & (c(j,1)*vel(iel2,1)
72  & +c(j,2)*vel(iel2,2)
73  & +c(j,3)*vel(iel2,3))
74  enddo
75  else
76  do j=1,3
77  vfap(j,i)=xl1*vel(iel1,j)+xl2*
78  & (c(1,j)*vel(iel2,1)
79  & +c(2,j)*vel(iel2,2)
80  & +c(3,j)*vel(iel2,3))
81  enddo
82  endif
83  elseif(ielfa(3,i).gt.0) then
84 !
85 ! boundary face; no zero gradient
86 !
87  iel3=ielfa(3,i)
88  ipointer=-iel2
89 !
90 ! global x-direction
91 !
92  if(ifabou(ipointer+1).gt.0) then
93 !
94 ! v_1 given
95 !
96  vfap(1,i)=xbounact(ifabou(ipointer+1))
97  else
98 !
99 ! extrapolation
100 !
101  vfap(1,i)=xl1*vel(iel1,1)+xrlfa(3,i)*vel(iel3,1)
102  endif
103 !
104 ! global y-direction
105 !
106  if(ifabou(ipointer+2).gt.0) then
107 !
108 ! v_2 given
109 !
110  vfap(2,i)=xbounact(ifabou(ipointer+2))
111  else
112 !
113 ! extrapolation
114 !
115  vfap(2,i)=xl1*vel(iel1,2)+xrlfa(3,i)*vel(iel3,2)
116  endif
117 !
118 ! global z-direction
119 !
120  if(ifabou(ipointer+3).gt.0) then
121 !
122 ! v_3 given
123 !
124  vfap(3,i)=xbounact(ifabou(ipointer+3))
125  else
126 !
127 ! extrapolation
128 !
129  vfap(3,i)=xl1*vel(iel1,3)+xrlfa(3,i)*vel(iel3,3)
130  endif
131 !
132 ! correction for sliding boundary conditions
133 !
134  if(ifabou(ipointer+5).lt.0) then
135  indexf=ipnei(iel1)+ielfa(4,i)
136  dd=vfap(1,i)*xxn(1,indexf)+
137  & vfap(2,i)*xxn(2,indexf)+
138  & vfap(3,i)*xxn(3,indexf)
139  do j=1,3
140  vfap(j,i)=vfap(j,i)-dd*xxn(j,indexf)
141  enddo
142  endif
143 !
144  else
145 !
146 ! boundary face; zero gradient
147 !
148  do j=1,3
149  vfap(j,i)=vel(iel1,j)
150  enddo
151  endif
152  enddo
153 !c$omp end do
154 !c$omp end parallel
155 !
156 ! Multiple point constraints
157 !
158  if(nmpc.gt.0) then
159  is=1
160  ie=3
161  call applympc(nface,ielfa,is,ie,ifabou,ipompc,vfap,coefmpc,
162  & nodempc,ipnei,neifa,labmpc,xbounact,nactdoh,
163  & ifaext,nfaext)
164  endif
165 !
166 ! calculate the gradient of the velocities at the center of
167 ! the elements
168 !
169 !c$omp parallel default(none)
170 !c$omp& shared(nef,ipnei,neifa,gradvel,vfap,area,xxn,volume)
171 !c$omp& private(i,indexf,ifa,k,l)
172 !c$omp do
173  do i=1,nef
174 !
175 ! initialization
176 !
177  do k=1,3
178  do l=1,3
179  gradvel(k,l,i)=0.d0
180  enddo
181  enddo
182 !
183  do indexf=ipnei(i)+1,ipnei(i+1)
184  ifa=neifa(indexf)
185  do k=1,3
186  do l=1,3
187  gradvel(k,l,i)=gradvel(k,l,i)+
188  & vfap(k,ifa)*area(ifa)*xxn(l,indexf)
189  enddo
190  enddo
191  enddo
192 !
193 ! dividing by the volume of the element
194 !
195  do k=1,3
196  do l=1,3
197  gradvel(k,l,i)=gradvel(k,l,i)/volume(i)
198  enddo
199  enddo
200  enddo
201 !c$omp end do
202 !c$omp end parallel
203 !
204 ! interpolate/extrapolate the velocity gradient from the
205 ! center of the elements to the center of the faces
206 !
207 !c$omp parallel default(none)
208 !c$omp& shared(nface,ielfa,xrlfa,gradvfa,gradvel,icyclic,c,ifatie,
209 !c$omp& indexf,ipnei,xxn)
210 !c$omp& private(i,iel1,xl1,iel2,k,l,xl2,gradnor)
211 !c$omp do
212  do i=1,nface
213  iel1=ielfa(1,i)
214  xl1=xrlfa(1,i)
215  iel2=ielfa(2,i)
216  if(iel2.gt.0) then
217 !
218 ! face between two elements
219 !
220  xl2=xrlfa(2,i)
221  if((icyclic.eq.0).or.(ifatie(i).eq.0)) then
222  do k=1,3
223  do l=1,3
224  gradvfa(k,l,i)=xl1*gradvel(k,l,iel1)+
225  & xl2*gradvel(k,l,iel2)
226  enddo
227  enddo
228  elseif(ifatie(i).gt.0) then
229  do k=1,3
230  do l=1,3
231  gradvfa(k,l,i)=xl1*gradvel(k,l,iel1)+xl2*
232  & (c(k,1)*gradvel(1,1,iel2)*c(l,1)+
233  & c(k,1)*gradvel(1,2,iel2)*c(l,2)+
234  & c(k,1)*gradvel(1,3,iel2)*c(l,3)+
235  & c(k,2)*gradvel(2,1,iel2)*c(l,1)+
236  & c(k,2)*gradvel(2,2,iel2)*c(l,2)+
237  & c(k,2)*gradvel(2,3,iel2)*c(l,3)+
238  & c(k,3)*gradvel(3,1,iel2)*c(l,1)+
239  & c(k,3)*gradvel(3,2,iel2)*c(l,2)+
240  & c(k,3)*gradvel(3,3,iel2)*c(l,3))
241  enddo
242  enddo
243  else
244  do k=1,3
245  do l=1,3
246  gradvfa(k,l,i)=xl1*gradvel(k,l,iel1)+xl2*
247  & (c(1,k)*gradvel(1,1,iel2)*c(1,l)+
248  & c(1,k)*gradvel(1,2,iel2)*c(2,l)+
249  & c(1,k)*gradvel(1,3,iel2)*c(3,l)+
250  & c(2,k)*gradvel(2,1,iel2)*c(1,l)+
251  & c(2,k)*gradvel(2,2,iel2)*c(2,l)+
252  & c(2,k)*gradvel(2,3,iel2)*c(3,l)+
253  & c(3,k)*gradvel(3,1,iel2)*c(1,l)+
254  & c(3,k)*gradvel(3,2,iel2)*c(2,l)+
255  & c(3,k)*gradvel(3,3,iel2)*c(3,l))
256  enddo
257  enddo
258  endif
259  elseif(ielfa(3,i).gt.0) then
260 !
261 ! boundary face; no zero gradient
262 !
263  do k=1,3
264  do l=1,3
265  gradvfa(k,l,i)=xl1*gradvel(k,l,iel1)+
266  & xrlfa(3,i)*gradvel(k,l,ielfa(3,i))
267  enddo
268  enddo
269  else
270 !
271 ! boundary face; zero gradient in i-direction
272 !
273  indexf=ipnei(iel1)+ielfa(4,i)
274  do k=1,3
275  gradnor=gradvel(k,1,iel1)*xxi(1,indexf)
276  & +gradvel(k,2,iel1)*xxi(2,indexf)
277  & +gradvel(k,3,iel1)*xxi(3,indexf)
278  do l=1,3
279  gradvfa(k,l,i)=gradvel(k,l,iel1)
280  & -gradnor*xxi(l,indexf)
281  enddo
282  enddo
283  endif
284  enddo
285 !c$omp end do
286 !c$omp end parallel
287 !
288 ! correction loops
289 !
290  do m=1,2
291 !
292 ! Moukalled et al. p 279
293 !
294 !c$omp parallel default(none)
295 !c$omp& shared(nface,ielfa,vfa,vfap,gradvfa,rf,ifabou,xxn,ipnei,xxi,
296 !c$omp& xle,vel)
297 !c$omp& private(i,iel1,iel2,ipointer,dd,j,indexf)
298 !c$omp do
299  do i=1,nface
300  iel1=ielfa(1,i)
301  iel2=ielfa(2,i)
302  if(iel2.gt.0) then
303 !
304 ! face between two elements
305 !
306  do j=1,3
307  vfa(j,i)=vfap(j,i)+gradvfa(j,1,i)*rf(1)+
308  & gradvfa(j,2,i)*rf(2)+
309  & gradvfa(j,3,i)*rf(3)
310  enddo
311  elseif(ielfa(3,i).gt.0) then
312 !
313 ! boundary face; no zero gradient
314 !
315  ipointer=-iel2
316 !
317 ! x-direction
318 !
319  if(ifabou(ipointer+1).gt.0) then
320  vfa(1,i)=vfap(1,i)
321  else
322  vfa(1,i)=vfap(1,i)+gradvfa(1,1,i)*rf(1)+
323  & gradvfa(1,2,i)*rf(2)+
324  & gradvfa(1,3,i)*rf(3)
325  endif
326 !
327 ! y-direction
328 !
329  if(ifabou(ipointer+2).gt.0) then
330  vfa(2,i)=vfap(2,i)
331  else
332  vfa(2,i)=vfap(2,i)+gradvfa(2,1,i)*rf(1)+
333  & gradvfa(2,2,i)*rf(2)+
334  & gradvfa(2,3,i)*rf(3)
335  endif
336 !
337 ! z-direction
338 !
339  if(ifabou(ipointer+3).gt.0) then
340  vfa(3,i)=vfap(3,i)
341  else
342  vfa(3,i)=vfap(3,i)+gradvfa(3,1,i)*rf(1)+
343  & gradvfa(3,2,i)*rf(2)+
344  & gradvfa(3,3,i)*rf(3)
345  endif
346 !
347 ! correction for sliding boundary conditions
348 !
349  if(ifabou(ipointer+5).lt.0) then
350  indexf=ipnei(iel1)+ielfa(4,i)
351  dd=vfa(1,i)*xxn(1,indexf)+
352  & vfa(2,i)*xxn(2,indexf)+
353  & vfa(3,i)*xxn(3,indexf)
354  do j=1,3
355  vfa(j,i)=vfa(j,i)-dd*xxn(j,indexf)
356  enddo
357  endif
358  else
359 !
360 ! boundary face; zero gradient
361 !
362  indexf=ipnei(iel1)+ielfa(4,i)
363  do j=1,3
364  vfa(j,i)=vel(iel1,j)+
365  & (gradvfa(j,1,i)*xxi(1,indexf)+
366  & gradvfa(j,2,i)*xxi(2,indexf)+
367  & gradvfa(j,3,i)*xxi(3,indexf))*xle(indexf)
368  enddo
369  endif
370  enddo
371 !c$omp end do
372 !c$omp end parallel
373 !
374 ! Multiple point constraints
375 !
376  if(nmpc.gt.0) then
377  is=1
378  ie=3
379  call applympc(nface,ielfa,is,ie,ifabou,ipompc,vfa,coefmpc,
380  & nodempc,ipnei,neifa,labmpc,xbounact,nactdoh,
381  & ifaext,nfaext)
382  endif
383 !
384 ! calculate the gradient of the velocities at the center of
385 ! the elements
386 !
387 !c$omp parallel default(none)
388 !c$omp& shared(nef,ipnei,neifa,gradvel,vfa,area,xxn,volume)
389 !c$omp& private(i,indexf,ifa,k,l)
390 !c$omp do
391  do i=1,nef
392 !
393 ! initialization
394 !
395  do k=1,3
396  do l=1,3
397  gradvel(k,l,i)=0.d0
398  enddo
399  enddo
400 !
401  do indexf=ipnei(i)+1,ipnei(i+1)
402  ifa=neifa(indexf)
403  do k=1,3
404  do l=1,3
405  gradvel(k,l,i)=gradvel(k,l,i)+
406  & vfa(k,ifa)*area(ifa)*xxn(l,indexf)
407  enddo
408  enddo
409  enddo
410 !
411 ! dividing by the volume of the element
412 !
413  do k=1,3
414  do l=1,3
415  gradvel(k,l,i)=gradvel(k,l,i)/volume(i)
416  enddo
417  enddo
418  enddo
419 !c$omp end do
420 !c$omp end parallel
421 !
422 ! interpolate/extrapolate the velocity gradient from the
423 ! center of the elements to the center of the faces
424 !
425 !c$omp parallel default(none)
426 !c$omp& shared(nface,ielfa,xrlfa,gradvfa,gradvel,icyclic,c,ifatie,
427 !c$omp& indexf,ipnei,xxn)
428 !c$omp& private(i,iel1,xl1,iel2,k,l,xl2)
429 !c$omp do
430  do i=1,nface
431  iel1=ielfa(1,i)
432  xl1=xrlfa(1,i)
433  iel2=ielfa(2,i)
434  if(iel2.gt.0) then
435 !
436 ! face in between two elements
437 !
438  xl2=xrlfa(2,i)
439  if((icyclic.eq.0).or.(ifatie(i).eq.0)) then
440  do k=1,3
441  do l=1,3
442  gradvfa(k,l,i)=xl1*gradvel(k,l,iel1)+
443  & xl2*gradvel(k,l,iel2)
444  enddo
445  enddo
446  elseif(ifatie(i).gt.0) then
447  do k=1,3
448  do l=1,3
449  gradvfa(k,l,i)=xl1*gradvel(k,l,iel1)+xl2*
450  & (c(k,1)*gradvel(1,1,iel2)*c(l,1)+
451  & c(k,1)*gradvel(1,2,iel2)*c(l,2)+
452  & c(k,1)*gradvel(1,3,iel2)*c(l,3)+
453  & c(k,2)*gradvel(2,1,iel2)*c(l,1)+
454  & c(k,2)*gradvel(2,2,iel2)*c(l,2)+
455  & c(k,2)*gradvel(2,3,iel2)*c(l,3)+
456  & c(k,3)*gradvel(3,1,iel2)*c(l,1)+
457  & c(k,3)*gradvel(3,2,iel2)*c(l,2)+
458  & c(k,3)*gradvel(3,3,iel2)*c(l,3))
459  enddo
460  enddo
461  else
462  do k=1,3
463  do l=1,3
464  gradvfa(k,l,i)=xl1*gradvel(k,l,iel1)+xl2*
465  & (c(1,k)*gradvel(1,1,iel2)*c(1,l)+
466  & c(1,k)*gradvel(1,2,iel2)*c(2,l)+
467  & c(1,k)*gradvel(1,3,iel2)*c(3,l)+
468  & c(2,k)*gradvel(2,1,iel2)*c(1,l)+
469  & c(2,k)*gradvel(2,2,iel2)*c(2,l)+
470  & c(2,k)*gradvel(2,3,iel2)*c(3,l)+
471  & c(3,k)*gradvel(3,1,iel2)*c(1,l)+
472  & c(3,k)*gradvel(3,2,iel2)*c(2,l)+
473  & c(3,k)*gradvel(3,3,iel2)*c(3,l))
474  enddo
475  enddo
476  endif
477  elseif(ielfa(3,i).gt.0) then
478 !
479 ! boundary face; no zero gradient
480 !
481  do k=1,3
482  do l=1,3
483  gradvfa(k,l,i)=xl1*gradvel(k,l,iel1)+
484  & xrlfa(3,i)*gradvel(k,l,ielfa(3,i))
485  enddo
486  enddo
487  else
488 !
489 ! boundary face; zero gradient in i-direction
490 !
491  indexf=ipnei(iel1)+ielfa(4,i)
492  do k=1,3
493  gradnor=gradvel(k,1,iel1)*xxi(1,indexf)
494  & +gradvel(k,2,iel1)*xxi(2,indexf)
495  & +gradvel(k,3,iel1)*xxi(3,indexf)
496  do l=1,3
497  gradvfa(k,l,i)=gradvel(k,l,iel1)
498  & -gradnor*xxi(l,indexf)
499  enddo
500  enddo
501  endif
502  enddo
503 !c$omp end do
504 !c$omp end parallel
505  enddo
506 !
507 ! correct the facial velocity gradients:
508 ! Moukalled et al. p 289
509 !
510 !c$omp parallel default(none)
511 !c$omp& shared(nface,ielfa,ipnei,vel,xlet,gradvfa,xxj)
512 !c$omp& private(i,iel2,iel1,indexf,dd,k,l)
513 !c$omp do
514  do i=1,nface
515  iel2=ielfa(2,i)
516  if(iel2.gt.0) then
517  iel1=ielfa(1,i)
518  indexf=ipnei(iel1)+ielfa(4,i)
519  do l=1,3
520  dd=(vel(iel2,l)-vel(iel1,l))/xlet(indexf)
521  & -gradvfa(l,1,i)*xxj(1,indexf)
522  & -gradvfa(l,2,i)*xxj(2,indexf)
523  & -gradvfa(l,3,i)*xxj(3,indexf)
524  do k=1,3
525  gradvfa(l,k,i)=gradvfa(l,k,i)+dd*xxj(k,indexf)
526  enddo
527  enddo
528  endif
529  enddo
530 !c$omp end do
531 !c$omp end parallel
532 !
533  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)