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

Go to the source code of this file.

Functions/Subroutines

subroutine extrapol_pel (nface, ielfa, xrlfa, vel, vfa, ifabou, xbounact, nef, gradpel, gradpfa, neifa, rf, area, volume, xle, xxi, icyclic, xxn, ipnei, ifatie, coefmpc, nmpc, labmpc, ipompc, nodempc, ifaext, nfaext, nactdoh)
 

Function/Subroutine Documentation

◆ extrapol_pel()

subroutine extrapol_pel ( 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)  gradpel,
real*8, dimension(3,*), intent(inout)  gradpfa,
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)  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 pressure element values to the faces
25 !
26  implicit none
27 !
28  character*20 labmpc(*)
29 !
30  integer nface,ielfa(4,*),ifabou(*),i,iel1,iel2,nef,ibou,
31  & neifa(*),icyclic,ifa,indexf,l,m,ipnei(*),ifatie(*),
32  & is,ie,nmpc,ipompc(*),nodempc(3,*),ifaext(*),nfaext,nactdoh(*)
33 !
34  real*8 xrlfa(3,*),vel(nef,0:7),vfa(0:7,*),xbounact(*),xl1,xl2,
35  & vfap(0:7,nface),gradpel(3,*),gradpfa(3,*),rf(3),area(*),
36  & volume(*),xle(*),xxi(3,*),c(3,3),gradnor,xxn(3,*),coefmpc(*)
37 !
38  intent(in) nface,ielfa,xrlfa,vel,
39  & ifabou,xbounact,nef,neifa,rf,area,volume,
40  & xle,xxi,icyclic,xxn,ipnei,ifatie,
41  & coefmpc,nmpc,labmpc,ipompc,nodempc,ifaext,nfaext,nactdoh
42 !
43  intent(inout) vfa,gradpel,gradpfa
44 !
45 !c$omp parallel default(none)
46 !c$omp& shared(nface,ielfa,xrlfa,vfap,vel,ifabou,xbounact)
47 !c$omp& private(i,iel1,xl1,iel2,ibou)
48 !c$omp do
49  do i=1,nface
50  iel1=ielfa(1,i)
51  xl1=xrlfa(1,i)
52  iel2=ielfa(2,i)
53  if(iel2.gt.0) then
54 !
55 ! face between two elements: interpolation
56 !
57  vfap(4,i)=xl1*vel(iel1,4)+xrlfa(2,i)*vel(iel2,4)
58  elseif(ielfa(3,i).ne.0) then
59 !
60 ! boundary face; more than one layer
61 !
62  ibou=0
63  if(iel2.lt.0) then
64  if(ifabou(-iel2+4).gt.0) then
65  ibou=ifabou(-iel2+4)
66  endif
67  endif
68 !
69  if(ibou.gt.0) then
70 !
71 ! pressure boundary condition
72 !
73  vfap(4,i)=xbounact(ibou)
74  else
75 !
76 ! extrapolation
77 !
78  vfap(4,i)=xl1*vel(iel1,4)
79  & +xrlfa(3,i)*vel(abs(ielfa(3,i)),4)
80  endif
81  else
82 !
83 ! boundary face; one layer
84 !
85  vfap(4,i)=vel(iel1,4)
86  endif
87  enddo
88 !c$omp end do
89 !c$omp end parallel
90 !
91 ! Multiple point constraints
92 !
93  if(nmpc.gt.0) then
94  is=4
95  ie=4
96  call applympc(nface,ielfa,is,ie,ifabou,ipompc,vfap,coefmpc,
97  & nodempc,ipnei,neifa,labmpc,xbounact,nactdoh,
98  & ifaext,nfaext)
99  endif
100 !
101 ! calculate the gradient of the pressure at the center of
102 ! the elements
103 !
104 !c$omp parallel default(none)
105 !c$omp& shared(nef,ipnei,neifa,gradpel,vfap,area,xxn,volume)
106 !c$omp& private(i,indexf,ifa)
107 !c$omp do
108  do i=1,nef
109 !
110 ! initialization
111 !
112  do l=1,3
113  gradpel(l,i)=0.d0
114  enddo
115 !
116  do indexf=ipnei(i)+1,ipnei(i+1)
117  ifa=neifa(indexf)
118  do l=1,3
119  gradpel(l,i)=gradpel(l,i)+
120  & vfap(4,ifa)*area(ifa)*xxn(l,indexf)
121  enddo
122  enddo
123 !
124 ! dividing by the volume of the element
125 !
126  do l=1,3
127  gradpel(l,i)=gradpel(l,i)/volume(i)
128  enddo
129  enddo
130 !c$omp end do
131 !c$omp end parallel
132 !
133 ! interpolate/extrapolate the pressure gradient from the
134 ! center of the elements to the center of the faces
135 !
136 !c$omp parallel default(none)
137 !c$omp& shared(nface,ielfa,xrlfa,gradpfa,gradpel,icyclic,c,ifatie,
138 !c$omp& ipnei,xxn)
139 !c$omp& private(i,iel1,xl1,iel2,l,xl2,indexf,gradnor)
140 !c$omp do
141  do i=1,nface
142  iel1=ielfa(1,i)
143  xl1=xrlfa(1,i)
144  iel2=ielfa(2,i)
145  if(iel2.gt.0) then
146 !
147 ! face in between two elements
148 !
149  xl2=xrlfa(2,i)
150  if((icyclic.eq.0).or.(ifatie(i).eq.0)) then
151  do l=1,3
152  gradpfa(l,i)=xl1*gradpel(l,iel1)+
153  & xl2*gradpel(l,iel2)
154  enddo
155  elseif(ifatie(i).gt.0) then
156  do l=1,3
157  gradpfa(l,i)=xl1*gradpel(l,iel1)+xl2*
158  & (gradpel(1,iel2)*c(l,1)+
159  & gradpel(2,iel2)*c(l,2)+
160  & gradpel(3,iel2)*c(l,3))
161  enddo
162  else
163  do l=1,3
164  gradpfa(l,i)=xl1*gradpel(l,iel1)+xl2*
165  & (gradpel(1,iel2)*c(1,l)+
166  & gradpel(2,iel2)*c(2,l)+
167  & gradpel(3,iel2)*c(3,l))
168  enddo
169  endif
170  elseif(ielfa(3,i).ne.0) then
171 !
172 ! boundary face; more than one layer; extrapolation
173 !
174  do l=1,3
175  gradpfa(l,i)=xl1*gradpel(l,iel1)+
176  & xrlfa(3,i)*gradpel(l,abs(ielfa(3,i)))
177  enddo
178  else
179 !
180 ! boundary face; one layer
181 !
182  indexf=ipnei(iel1)+ielfa(4,i)
183  gradnor=gradpel(1,iel1)*xxi(1,indexf)+
184  & gradpel(2,iel1)*xxi(2,indexf)+
185  & gradpel(3,iel1)*xxi(3,indexf)
186  do l=1,3
187  gradpfa(l,i)=gradpel(l,iel1)
188  & -gradnor*xxi(l,indexf)
189  enddo
190  endif
191  enddo
192 !c$omp end do
193 !c$omp end parallel
194 !
195 ! correction loops
196 !
197  do m=1,2
198 !
199 ! Moukalled et al. p 279
200 !
201 !c$omp parallel default(none)
202 !c$omp& shared(nface,ielfa,xrlfa,vfa,vfap,gradpfa,rf,ifabou,ipnei,
203 !c$omp& xxi,xle,vel)
204 !c$omp& private(i,iel1,iel2,xl1,ibou,indexf)
205 !c$omp do
206  do i=1,nface
207  iel1=ielfa(1,i)
208  iel2=ielfa(2,i)
209  xl1=xrlfa(1,i)
210  if(iel2.gt.0) then
211 !
212 ! face between two elements
213 !
214  vfa(4,i)=vfap(4,i)+gradpfa(1,i)*rf(1)
215  & +gradpfa(2,i)*rf(2)
216  & +gradpfa(3,i)*rf(3)
217  elseif(ielfa(3,i).ne.0) then
218 !
219 ! boundary face; more than one layer
220 !
221  ibou=0
222  if(iel2.lt.0) then
223  if(ifabou(-iel2+4).gt.0) then
224  ibou=ifabou(-iel2+4)
225  endif
226  endif
227 !
228  if(ibou.gt.0) then
229 !
230 ! pressure given
231 !
232  vfa(4,i)=vfap(4,i)
233  else
234 !
235 ! extrapolation
236 !
237  vfa(4,i)=vfap(4,i)+gradpfa(1,i)*rf(1)
238  & +gradpfa(2,i)*rf(2)
239  & +gradpfa(3,i)*rf(3)
240  endif
241  else
242 !
243 ! boundary face; one layer
244 !
245  indexf=ipnei(iel1)+ielfa(4,i)
246  vfa(4,i)=vel(iel1,4)
247  & +(gradpfa(1,i)*xxi(1,indexf)+
248  & gradpfa(2,i)*xxi(2,indexf)+
249  & gradpfa(3,i)*xxi(3,indexf))*xle(indexf)
250  endif
251  enddo
252 !c$omp end do
253 !c$omp end parallel
254 !
255 ! Multiple point constraints
256 !
257  if(nmpc.gt.0) then
258  is=4
259  ie=4
260  call applympc(nface,ielfa,is,ie,ifabou,ipompc,vfa,coefmpc,
261  & nodempc,ipnei,neifa,labmpc,xbounact,nactdoh,
262  & ifaext,nfaext)
263  endif
264 !
265 ! calculate the gradient of the pressure at the center of
266 ! the elements
267 !
268 !c$omp parallel default(none)
269 !c$omp& shared(nef,ipnei,neifa,gradpel,vfa,area,xxn,volume)
270 !c$omp& private(i,indexf,ifa)
271 !c$omp do
272  do i=1,nef
273 !
274 ! initialization
275 !
276  do l=1,3
277  gradpel(l,i)=0.d0
278  enddo
279 !
280  do indexf=ipnei(i)+1,ipnei(i+1)
281  ifa=neifa(indexf)
282  do l=1,3
283  gradpel(l,i)=gradpel(l,i)+
284  & vfa(4,ifa)*area(ifa)*xxn(l,indexf)
285  enddo
286  enddo
287 !
288 ! dividing by the volume of the element
289 !
290  do l=1,3
291  gradpel(l,i)=gradpel(l,i)/volume(i)
292  enddo
293  enddo
294 !c$omp end do
295 !c$omp end parallel
296 !
297 ! interpolate/extrapolate the pressure gradient from the
298 ! center of the elements to the center of the faces
299 !
300 !c$omp parallel default(none)
301 !c$omp& shared(nface,ielfa,xrlfa,gradpfa,gradpel,icyclic,c,ifatie,
302 !c$omp& ipnei,xxn)
303 !c$omp& private(i,iel1,xl1,iel2,l,xl2,indexf,gradnor)
304 !c$omp do
305  do i=1,nface
306  iel1=ielfa(1,i)
307  xl1=xrlfa(1,i)
308  iel2=ielfa(2,i)
309  if(iel2.gt.0) then
310 !
311 ! face in between two elements
312 !
313  xl2=xrlfa(2,i)
314  if((icyclic.eq.0).or.(ifatie(i).eq.0)) then
315  do l=1,3
316  gradpfa(l,i)=xl1*gradpel(l,iel1)+
317  & xl2*gradpel(l,iel2)
318  enddo
319  elseif(ifatie(i).gt.0) then
320  do l=1,3
321  gradpfa(l,i)=xl1*gradpel(l,iel1)+xl2*
322  & (gradpel(1,iel2)*c(l,1)+
323  & gradpel(2,iel2)*c(l,2)+
324  & gradpel(3,iel2)*c(l,3))
325  enddo
326  else
327  do l=1,3
328  gradpfa(l,i)=xl1*gradpel(l,iel1)+xl2*
329  & (gradpel(1,iel2)*c(1,l)+
330  & gradpel(2,iel2)*c(2,l)+
331  & gradpel(3,iel2)*c(3,l))
332  enddo
333  endif
334  elseif(ielfa(3,i).ne.0) then
335 !
336 ! boundary face; linear extrapolation
337 !
338  do l=1,3
339  gradpfa(l,i)=xl1*gradpel(l,iel1)+
340  & xrlfa(3,i)*gradpel(l,abs(ielfa(3,i)))
341  enddo
342  else
343 !
344 ! boundary face; constant extrapolation (one element layer)
345 !
346  indexf=ipnei(iel1)+ielfa(4,i)
347  gradnor=gradpel(1,iel1)*xxi(1,indexf)+
348  & gradpel(2,iel1)*xxi(2,indexf)+
349  & gradpel(3,iel1)*xxi(3,indexf)
350  do l=1,3
351  gradpfa(l,i)=gradpel(l,iel1)
352  & -gradnor*xxi(l,indexf)
353  enddo
354  endif
355  enddo
356 !c$omp end do
357 !c$omp end parallel
358 !
359  enddo
360 !
361  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)