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

Go to the source code of this file.

Functions/Subroutines

subroutine extrapolatecontact (yi, yn, ipkon, inum, kon, lakon, nfield, nk, ne, mi, ndim, co, cflag, vold, force, pslavsurf, islavact, islavnode, nslavnode, ntie, islavsurf, ielprop, prop, ielmat, ne0)
 

Function/Subroutine Documentation

◆ extrapolatecontact()

subroutine extrapolatecontact ( real*8, dimension(ndim,mi(1),*)  yi,
real*8, dimension(nfield,*)  yn,
integer, dimension(*)  ipkon,
integer, dimension(*)  inum,
integer, dimension(*)  kon,
character*8, dimension(*)  lakon,
integer  nfield,
integer  nk,
integer  ne,
integer, dimension(*)  mi,
integer  ndim,
real*8, dimension(3,*)  co,
character*1  cflag,
real*8, dimension(0:mi(2),*)  vold,
logical  force,
real*8, dimension(3,*)  pslavsurf,
integer, dimension(*)  islavact,
integer, dimension(*)  islavnode,
integer, dimension(*)  nslavnode,
integer  ntie,
integer, dimension(2,*)  islavsurf,
integer, dimension(*)  ielprop,
real*8, dimension(*)  prop,
integer, dimension(mi(3),*)  ielmat,
integer  ne0 
)
22 !
23 ! extrapolates contact values at the integration points to the
24 ! nodes (for face-to-face penalty contact)
25 !
26  implicit none
27 !
28  logical force
29 !
30  character*1 cflag
31  character*8 lakon(*),lakonl
32 !
33  integer ipkon(*),inum(*),kon(*),mi(*),ne,indexe,nope,
34  & nonei20(3,12),nfield,nonei10(3,6),nk,i,j,k,l,ndim,
35  & nonei15(3,9),m,iflag,jj,indexc,islavsurf(2,*),ll,
36  & indexcj,ifacew1(4,5),islavnode(*),nslavnode(*),ntie,
37  & nlayer,nopeexp,ifacew2(8,5),ifaceq(8,6),ifacet(6,4),
38  & nopespring,ifaces,nopespringj,ifacej,jfaces,n,nelems,
39  & idgn,idgnr,idglda,idgip(4),idgldb,idginfo,igauss,islavact(*),
40  & konl(26),node,nopes,ielprop(*),ielmat(mi(3),*),ne0,
41  & nodes(8)
42 !
43  real*8 yi(ndim,mi(1),*),yn(nfield,*),field(nfield,20*mi(3)),
44  & co(3,*),xi,et,vold(0:mi(2),*),xs2(3,7),xsj2(3),shp2(7,8),
45  & aa(4,4),bb(4,6),pl(3,8),pslavsurf(3,*),xslavnor(3,nk),
46  & cc(3,4),dd(3,6),c_limit(2,nfield),nodepos(4,2),xn(3),
47  & t1(3),t2(3),trac(3),xquad(2,9),xtri(2,7),xl2s(3,9),
48  & stn(6,nk),dt1,dl,a2(6,2),a4(4,4),a27(20,27),a9(6,9),
49  & a8(8,8),prop(*)
50 !
51  data nonei10 /5,1,2,6,2,3,7,3,1,8,1,4,9,2,4,10,3,4/
52 !
53  data nonei15 /7,1,2,8,2,3,9,3,1,10,4,5,11,5,6,12,6,4,
54  & 13,1,4,14,2,5,15,3,6/
55 !
56  data nonei20 /9,1,2,10,2,3,11,3,4,12,4,1,
57  & 13,5,6,14,6,7,15,7,8,16,8,5,
58  & 17,1,5,18,2,6,19,3,7,20,4,8/
59 !
60 ! nodes per face for hex elements
61 !
62  data ifaceq /4,3,2,1,11,10,9,12,
63  & 5,6,7,8,13,14,15,16,
64  & 1,2,6,5,9,18,13,17,
65  & 2,3,7,6,10,19,14,18,
66  & 3,4,8,7,11,20,15,19,
67  & 4,1,5,8,12,17,16,20/
68 !
69 ! nodes per face for tet elements
70 !
71  data ifacet /1,3,2,7,6,5,
72  & 1,2,4,5,9,8,
73  & 2,3,4,6,10,9,
74  & 1,4,3,8,10,7/
75 !
76 ! nodes per face for linear wedge elements
77 !
78  data ifacew1 /1,3,2,0,
79  & 4,5,6,0,
80  & 1,2,5,4,
81  & 2,3,6,5,
82  & 3,1,4,6/
83 !
84 ! nodes per face for quadratic wedge elements
85 !
86  data ifacew2 /1,3,2,9,8,7,0,0,
87  & 4,5,6,10,11,12,0,0,
88  & 1,2,5,4,7,14,10,13,
89  & 2,3,6,5,8,15,11,14,
90  & 3,1,4,6,9,13,12,15/
91 !
92  data iflag /2/
93 !
94  if(nfield.eq.0) return
95 !
96  do i=1,nk
97  inum(i)=0
98  enddo
99 !
100  do i=1,nk
101  do j=1,nfield
102  yn(j,i)=0.d0
103  enddo
104  enddo
105 !
106  i=0
107  do
108  i=i+1
109  if(i.gt.ne) exit
110 !
111  if(ipkon(i).lt.0) cycle
112  indexc=ipkon(i)
113 !
114  if((lakon(i)(1:1).ne.'E').or.(lakon(i)(7:7).ne.'C')) cycle
115  nopespring=kon(indexc)
116  ifaces=islavsurf(1,kon(indexc+nopespring+2))
117 !
118  nelems=int(ifaces/10.d0)
119  lakonl=lakon(nelems)
120 !
121 ! determining all contact elements belonging to slave face
122 ! "ifaces"
123 !
124  n=i
125 !
126  do
127  n=n+1
128  if(n.gt.ne) exit
129  indexcj=ipkon(n)
130  nopespringj=kon(indexcj)
131  ifacej=islavsurf(1,kon(indexcj+nopespringj+2))
132  if(ifaces.ne.ifacej) exit
133  enddo
134  n=n-1
135  jfaces=ifaces-10*int(ifaces/10.d0)
136  indexe=ipkon(nelems)
137 !
138  if(lakonl(4:4).eq.'2') then
139  nope=20
140  elseif(lakonl(4:4).eq.'8') then
141  nope=8
142  elseif(lakonl(4:5).eq.'10') then
143  nope=10
144  elseif(lakonl(4:4).eq.'4') then
145  nope=4
146  elseif(lakonl(4:5).eq.'15') then
147  nope=15
148  elseif(lakonl(4:4).eq.'6') then
149  nope=6
150  elseif((lakon(i)(1:1).eq.'E').and.(lakon(i)(7:7).eq.'A')) then
151  inum(kon(indexe+1))=inum(kon(indexe+1))+1
152  inum(kon(indexe+2))=inum(kon(indexe+2))+1
153  cycle
154  else
155  cycle
156  endif
157 !
158  if((lakonl(4:5).eq.'20').or.(lakonl(4:4).eq.'8').or.
159  & (((lakonl(4:5).eq.'15').or.(lakonl(4:4).eq.'6')).and.
160  & (jfaces.gt.2))) then
161  if(lakonl(7:8).ne.'LC') then
162  field(1:nfield,1:20)=0.d0
163  do k=1,4
164  do l=1,4
165  aa(k,l)=0.d0
166  enddo
167  do l=1,nfield
168  bb(k,l)=0.d0
169  enddo
170  enddo
171  nodepos(1,1)=-1.d0
172  nodepos(1,2)=-1.d0
173  nodepos(2,1)=1.d0
174  nodepos(2,2)=-1.d0
175  nodepos(3,1)=1.d0
176  nodepos(3,2)=1.d0
177  nodepos(4,1)=-1.d0
178  nodepos(4,2)=1.d0
179  do j=i,n
180  do k=1,nfield
181  if(j.eq.i) then
182  c_limit(1,k)=yi(k,1,j)
183  c_limit(2,k)=yi(k,1,j)
184  endif
185  if(c_limit(1,k).lt.yi(k,1,j)) then
186  c_limit(1,k)=yi(k,1,j)
187  endif
188  if(c_limit(2,k).gt.yi(k,1,j)) then
189  c_limit(2,k)=yi(k,1,j)
190  endif
191  enddo
192  indexcj=ipkon(j)
193  nopespringj=kon(indexcj)
194  igauss=kon(indexcj+nopespringj+1)
195  xi=pslavsurf(1,igauss)
196  et=pslavsurf(2,igauss)
197  if((n-i).lt.2) exit
198  if((n-i).eq.2) then
199  aa(j+1-i,1)=xi
200  aa(j+1-i,2)=et
201  do k=1,nfield
202  dd(j+1-i,k)=yi(k,1,j)
203  enddo
204  elseif((n-i).eq.3) then
205  call shape4q(xi,et,pl,xsj2,xs2,shp2,iflag)
206  do k=1,4
207  aa(j+1-i,k)=shp2(4,k)
208  enddo
209  do k=1,nfield
210  bb(j+1-i,k)=yi(k,1,j)
211  enddo
212  else
213  call shape4q(xi,et,pl,xsj2,xs2,shp2,iflag)
214  do k=1,4
215  do l=k,4
216  aa(k,l)=aa(k,l)+shp2(4,k)*shp2(4,l)
217  enddo
218  enddo
219  do k=1,4
220  do l=1,nfield
221  bb(k,l)=bb(k,l)+shp2(4,k)*yi(l,1,j)
222  enddo
223  enddo
224  endif
225  enddo
226  if(n.eq.i) then
227  do k=1,4
228  do l=1,nfield
229  bb(k,l)=bb(k,l)+yi(l,1,n)
230  enddo
231  enddo
232  elseif((n-i).eq.1) then
233  do j=i,n
234  do k=1,4
235  do l=1,nfield
236  bb(k,l)=bb(k,l)+yi(l,1,j)*0.5d0
237  enddo
238  enddo
239  enddo
240  elseif((n-i).eq.2) then
241  do k=1,4
242  do l=1,nfield
243  call plane_eq(aa(1,1),aa(1,2),dd(1,l),
244  & aa(2,1),aa(2,2),dd(2,l),aa(3,1),
245  & aa(3,2),dd(3,l),nodepos(k,1),
246  & nodepos(k,2),bb(k,l))
247  enddo
248  enddo
249  else
250  if((n-i).ne.3) then
251  do k=1,4
252  do l=1,k-1
253  aa(k,l)=aa(l,k)
254  enddo
255  enddo
256  endif
257  idgn=4
258  idgnr=nfield
259  idglda=4
260  idgldb=4
261  call dgesv(idgn,idgnr,aa,idglda,idgip,bb,idgldb,
262  & idginfo)
263  endif
264  if((lakonl(4:4).eq.'6').or.(lakonl(4:5).eq.'15')) then
265  do j=1,4
266  do k=1,nfield
267  if((c_limit(1,k).gt.bb(j,k)).and.
268  & (c_limit(2,k).lt.bb(j,k))) then
269  field(k,ifacew1(j,jfaces))=bb(j,k)
270  endif
271  if(c_limit(1,k).lt.bb(j,k)) then
272  field(k,ifacew1(j,jfaces))=c_limit(1,k)
273  endif
274  if(c_limit(2,k).gt.bb(j,k)) then
275  field(k,ifacew1(j,jfaces))=c_limit(2,k)
276  endif
277  enddo
278  enddo
279  else
280  do j=1,4
281  do k=1,nfield
282  if((c_limit(1,k).gt.bb(j,k)).and.
283  & (c_limit(2,k).lt.bb(j,k))) then
284  field(k,ifaceq(j,jfaces))=bb(j,k)
285  endif
286  if(c_limit(1,k).lt.bb(j,k)) then
287  field(k,ifaceq(j,jfaces))=c_limit(1,k)
288  endif
289  if(c_limit(2,k).gt.bb(j,k)) then
290  field(k,ifaceq(j,jfaces))=c_limit(2,k)
291  endif
292  enddo
293  enddo
294  endif
295  endif
296  elseif((lakonl(4:5).eq.'10').or.(lakonl(4:4).eq.'4').or.
297  & (((lakonl(4:5).eq.'15').or.(lakonl(4:4).eq.'6')).and.
298  & (jfaces.le.2))) then
299  field(1:nfield,1:15)=0.d0
300  if(lakonl(7:8).ne.'LC') then
301  do k=1,3
302  do l=1,3
303  cc(k,l)=0.d0
304  enddo
305  do l=1,nfield
306  dd(k,l)=0.d0
307  enddo
308  enddo
309  do j=i,n
310  do k=1,nfield
311  if(j.eq.i) then
312  c_limit(1,k)=yi(k,1,j)
313  c_limit(2,k)=yi(k,1,j)
314  endif
315  if(c_limit(1,k).lt.yi(k,1,j)) then
316  c_limit(1,k)=yi(k,1,j)
317  endif
318  if(c_limit(2,k).gt.yi(k,1,j)) then
319  c_limit(2,k)=yi(k,1,j)
320  endif
321  enddo
322  indexcj=ipkon(j)
323  nopespringj=kon(indexcj)
324  igauss=kon(indexcj+nopespringj+1)
325  xi=pslavsurf(1,igauss)
326  et=pslavsurf(2,igauss)
327  if((n-i).lt.2) exit
328  if((n-i).eq.2) then
329  call shape3tri(xi,et,pl,xsj2,xs2,shp2,iflag)
330  do k=1,3
331  cc(j+1-i,k)=shp2(4,k)
332  enddo
333  do k=1,nfield
334  dd(j+1-i,k)=yi(k,1,j)
335  enddo
336  else
337  call shape3tri(xi,et,pl,xsj2,xs2,shp2,iflag)
338  do k=1,3
339  do l=1,3
340  cc(k,l)=cc(k,l)+shp2(4,k)*shp2(4,l)
341  enddo
342  enddo
343  do k=1,3
344  do l=1,nfield
345  dd(k,l)=dd(k,l)+shp2(4,k)*yi(l,1,j)
346  enddo
347  enddo
348  endif
349  enddo
350  if(n.eq.i) then
351  do k=1,3
352  do l=1,nfield
353  dd(k,l)=dd(k,l)+yi(l,1,j)
354  enddo
355  enddo
356  elseif((n-i).eq.1) then
357  do j=i,n
358  do k=1,3
359  do l=1,nfield
360  dd(k,l)=dd(k,l)+yi(l,1,j)*0.5d0
361  enddo
362  enddo
363  enddo
364  else
365  idgn=3
366  idgnr=nfield
367  idglda=3
368  idgldb=3
369  call dgesv(idgn,idgnr,cc,idglda,idgip,dd,idgldb,
370  & idginfo)
371  endif
372  if((lakonl(4:4).eq.'6').or.(lakonl(4:5).eq.'15')) then
373  do j=1,3
374  do k=1,nfield
375  if((c_limit(1,k).gt.dd(j,k)).and.
376  & (c_limit(2,k).lt.dd(j,k))) then
377  field(k,ifacew1(j,jfaces))=dd(j,k)
378  endif
379  if(c_limit(1,k).lt.dd(j,k)) then
380  field(k,ifacew1(j,jfaces))=c_limit(1,k)
381  endif
382  if(c_limit(2,k).gt.dd(j,k)) then
383  field(k,ifacew1(j,jfaces))=c_limit(2,k)
384  endif
385  enddo
386  enddo
387  else
388  do j=1,3
389  do k=1,nfield
390  if((c_limit(1,k).gt.dd(j,k)).and.
391  & (c_limit(2,k).lt.dd(j,k))) then
392  field(k,ifacet(j,jfaces))=dd(j,k)
393  endif
394  if(c_limit(1,k).lt.dd(j,k)) then
395  field(k,ifacet(j,jfaces))=c_limit(1,k)
396  endif
397  if(c_limit(2,k).gt.dd(j,k)) then
398  field(k,ifacet(j,jfaces))=c_limit(2,k)
399  endif
400  enddo
401  enddo
402  endif
403  endif
404  endif
405 !
406 ! determining the field values in the midside nodes of the
407 ! slave face
408 !
409  if(lakonl(4:5).eq.'20') then
410  if(lakonl(7:8).ne.'LC') then
411  do j=5,8
412  do k=1,nfield
413  field(k,ifaceq(j,jfaces))=
414  & (field(k,nonei20(2,ifaceq(j,jfaces)-8))+
415  & field(k,nonei20(3,ifaceq(j,jfaces)-8)))/2.d0
416  enddo
417  enddo
418  else
419  nlayer=0
420  do j=1,mi(3)
421  if(ielmat(j,i).gt.0) then
422  nlayer=nlayer+1
423  else
424  exit
425  endif
426  enddo
427  do m=1,nlayer
428  jj=20*(m-1)
429  do j=9,20
430  do k=1,nfield
431  field(k,jj+j)=(field(k,jj+nonei20(2,j-8))
432  & +field(k,jj+nonei20(3,j-8)))/2.d0
433  enddo
434  enddo
435  enddo
436  endif
437  elseif(lakonl(4:5).eq.'10') then
438  do j=4,6
439  do k=1,nfield
440  field(k,ifacet(j,jfaces))=
441  & (field(k,nonei10(2,ifacet(j,jfaces)-4))+
442  & field(k,nonei10(3,ifacet(j,jfaces)-4)))/2.d0
443  enddo
444  enddo
445  elseif(lakonl(4:5).eq.'15') then
446  if(lakonl(7:8).ne.'LC') then
447  if(jfaces.le.2) then
448  do j=4,6
449  do k=1,nfield
450  field(k,ifacew2(j,jfaces))=
451  & (field(k,nonei15(2,ifacew2(j,jfaces)-6))+
452  & field(k,nonei15(3,ifacew2(j,jfaces)-6)))
453  & /2.d0
454  enddo
455  enddo
456  else
457  do j=5,8
458  do k=1,nfield
459  field(k,ifacew2(j,jfaces))=
460  & (field(k,nonei15(2,ifacew2(j,jfaces)-6))+
461  & field(k,nonei15(3,ifacew2(j,jfaces)-6)))
462  & /2.d0
463  enddo
464  enddo
465  endif
466  else
467  nlayer=0
468  do j=1,mi(3)
469  if(ielmat(j,i).gt.0) then
470  nlayer=nlayer+1
471  else
472  exit
473  endif
474  enddo
475  do m=1,nlayer
476  jj=15*(m-1)
477  do j=7,15
478  do k=1,nfield
479  field(k,jj+j)=(field(k,jj+nonei15(2,j-6))
480  & +field(k,jj+nonei15(3,j-6)))/2.d0
481  enddo
482  enddo
483  enddo
484  endif
485  endif
486 !
487 ! transferring the field values into yn
488 !
489  if(lakonl(7:8).ne.'LC') then
490  do j=1,nope
491  do k=1,nfield-2
492  yn(k,kon(indexe+j))=yn(k,kon(indexe+j))+
493  & field(k,j)
494  enddo
495 !
496 ! interchanging positions of two last rows of yn
497 !
498  yn(nfield-1,kon(indexe+j))=yn(nfield-1,kon(indexe+j))+
499  & field(nfield,j)
500  yn(nfield,kon(indexe+j))=yn(nfield,kon(indexe+j))+
501  & field(nfield-1,j)
502  inum(kon(indexe+j))=inum(kon(indexe+j))+1
503  enddo
504  else
505  do j=1,nope*nlayer
506  do k=1,nfield
507  yn(k,kon(indexe+nopeexp+j))=
508  & yn(k,kon(indexe+nopeexp+j))+field(k,j)
509  enddo
510  inum(kon(indexe+nopeexp+j))=inum(kon(indexe+nopeexp+j))+1
511  enddo
512  endif
513 !
514 c Bernhardi start
515 c incompatible modes elements
516  if(lakonl(1:5).eq.'C3D8I') then
517  do j=1,3
518  do k=1,nfield
519  yn(k,kon(indexe+nope+j))=0.0d0
520  enddo
521 c inum(kon(indexe+nope+j))=inum(kon(indexe+nope+j))+1
522  enddo
523  endif
524 c Bernhardi end
525 !
526  i=n
527 !
528  enddo
529 !
530 ! taking the mean of nodal contributions coming from different
531 ! elements having the node in common
532 !
533  do i=1,nk
534  if(inum(i).gt.0) then
535  do j=1,nfield
536  yn(j,i)=yn(j,i)/inum(i)
537  enddo
538  endif
539  enddo
540 !
541 ! zeroing nonactive nodes
542 !
543  do i=1,nslavnode(ntie+1)
544  if(islavact(i).ne.1) then
545  do j=1,nfield
546  yn(j,islavnode(i))=0.d0
547  enddo
548  endif
549  enddo
550 !
551 ! for 1d and 2d elements only:
552 ! finding the solution in the original nodes
553 !
554  if((cflag.ne.' ').and.(cflag.ne.'E')) then
555  call map3dto1d2d(yn,ipkon,inum,kon,lakon,nfield,nk,ne,cflag,co,
556  & vold,force,mi)
557  endif
558 !
559  return
plane_eq
subroutine plane_eq(x1, y1, z1, x2, y2, z2, x3, y3, z3, x0, y0, output)
Definition: plane_eq.f:33
shape3tri
subroutine shape3tri(xi, et, xl, xsj, xs, shp, iflag)
Definition: shape3tri.f:20
dgesv
subroutine dgesv(N, NRHS, A, LDA, IPIV, B, LDB, INFO)
Definition: dgesv.f:58
shape4q
subroutine shape4q(xi, et, xl, xsj, xs, shp, iflag)
Definition: shape4q.f:20
map3dto1d2d
subroutine map3dto1d2d(yn, ipkon, inum, kon, lakon, nfield, nk, ne, cflag, co, vold, force, mi)
Definition: map3dto1d2d.f:21
nodes
subroutine nodes(inpc, textpart, co, nk, nk_, set, istartset, iendset, ialset, nset, nset_, nalset, nalset_, istep, istat, n, iline, ipol, inl, ipoinp, inp, ipoinpc)
Definition: nodes.f:22
Hosted by OpenAircraft.com, (Michigan UAV, LLC)