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

Go to the source code of this file.

Functions/Subroutines

subroutine gencontelem_f2f (tieset, ntie, itietri, ne, ipkon, kon, lakon, cg, straight, ifree, koncont, co, vold, xo, yo, zo, x, y, z, nx, ny, nz, ielmat, elcon, istep, iinc, iit, ncmat_, ntmat_, mi, imastop, islavsurf, itiefac, springarea, tietol, reltime, filab, nasym, pslavsurf, pmastsurf, clearini, theta, xstateini, xstate, nstate_, ne0, icutb, ialeatoric, nmethod, jobnamef)
 

Function/Subroutine Documentation

◆ gencontelem_f2f()

subroutine gencontelem_f2f ( character*81, dimension(3,*)  tieset,
integer  ntie,
integer, dimension(2,ntie)  itietri,
integer  ne,
integer, dimension(*)  ipkon,
integer, dimension(*)  kon,
character*8, dimension(*)  lakon,
real*8, dimension(3,*)  cg,
real*8, dimension(16,*)  straight,
integer  ifree,
integer, dimension(4,*)  koncont,
real*8, dimension(3,*)  co,
real*8, dimension(0:mi(2),*)  vold,
real*8, dimension(*)  xo,
real*8, dimension(*)  yo,
real*8, dimension(*)  zo,
real*8, dimension(*)  x,
real*8, dimension(*)  y,
real*8, dimension(*)  z,
integer, dimension(*)  nx,
integer, dimension(*)  ny,
integer, dimension(*)  nz,
integer, dimension(mi(3),*)  ielmat,
real*8, dimension(0:ncmat_,ntmat_,*)  elcon,
integer  istep,
integer  iinc,
integer  iit,
integer  ncmat_,
integer  ntmat_,
integer, dimension(*)  mi,
integer, dimension(3,*)  imastop,
integer, dimension(2,*)  islavsurf,
integer, dimension(2,*)  itiefac,
real*8, dimension(2,*)  springarea,
real*8, dimension(3,*)  tietol,
real*8  reltime,
character*87, dimension(*)  filab,
integer  nasym,
real*8, dimension(3,*)  pslavsurf,
real*8, dimension(6,*)  pmastsurf,
real*8, dimension(3,9,*)  clearini,
real*8  theta,
real*8, dimension(nstate_,mi(1),*)  xstateini,
real*8, dimension(nstate_,mi(1),*)  xstate,
integer  nstate_,
integer  ne0,
integer  icutb,
integer  ialeatoric,
integer  nmethod,
character*132, dimension(*)  jobnamef 
)
25 !
26 ! generate contact elements for the slave contact nodes
27 !
28  implicit none
29 !
30  character*8 lakon(*)
31  character*33 cfile
32  character*81 tieset(3,*),setname
33  character*87 filab(*)
34  character*132 jobnamef(*)
35 !
36  integer ntie,ifree,nasym,nstate_,ne0,
37  & itietri(2,ntie),ipkon(*),kon(*),koncont(4,*),ne,
38  & neigh(1),iflag,kneigh,i,j,k,l,jj,nn,isol,icutb,
39  & itri,ll,kflag,n,nx(*),ny(*),istep,iinc,mi(*),
40  & nz(*),nstart,ielmat(mi(3),*),imat,ifaceq(8,6),ifacet(6,4),
41  & ifacew1(4,5),ifacew2(8,5),nelemm,jfacem,indexe,iit,
42  & nface,nope,nodefm(9),ncmat_,ntmat_,number(4),lenset,
43  & iteller,ifaces,jfaces,ifacem,indexel,iloop,iprev,iact,
44  & imastop(3,*), itriangle(100),ntriangle,ntriangle_,itriold,
45  & itrinew,id,islavsurf(2,*),itiefac(2,*),nelems,m,mint2d,nopes,
46  & iloc,nopem,nodefs(9),indexf,ialeatoric,nmethod
47 !
48  real*8 cg(3,*),straight(16,*),co(3,*),vold(0:mi(2),*),p(3),
49  & dist,xo(*),yo(*),zo(*),x(*),y(*),z(*),clearini(3,9,*),
50  & elcon(0:ncmat_,ntmat_,*),weight,theta,harvest,
51  & springarea(2,*),xl2(3,9),area,xi,et,shp2(7,9),
52  & xs2(3,2),xsj2(3),tietol(3,*),reltime,xstate(nstate_,mi(1),*),
53  & clear,ratio(9),pl(3,9),xstateini(nstate_,mi(1),*),
54  & pproj(3),al(3),xn(3),xm(3),dm,pslavsurf(3,*),pmastsurf(6,*)
55 !
56 ! nodes per face for hex elements
57 !
58  data ifaceq /4,3,2,1,11,10,9,12,
59  & 5,6,7,8,13,14,15,16,
60  & 1,2,6,5,9,18,13,17,
61  & 2,3,7,6,10,19,14,18,
62  & 3,4,8,7,11,20,15,19,
63  & 4,1,5,8,12,17,16,20/
64 !
65 ! nodes per face for tet elements
66 !
67  data ifacet /1,3,2,7,6,5,
68  & 1,2,4,5,9,8,
69  & 2,3,4,6,10,9,
70  & 1,4,3,8,10,7/
71 !
72 ! nodes per face for linear wedge elements
73 !
74  data ifacew1 /1,3,2,0,
75  & 4,5,6,0,
76  & 1,2,5,4,
77  & 2,3,6,5,
78  & 3,1,4,6/
79 !
80 ! nodes per face for quadratic wedge elements
81 !
82  data ifacew2 /1,3,2,9,8,7,0,0,
83  & 4,5,6,10,11,12,0,0,
84  & 1,2,5,4,7,14,10,13,
85  & 2,3,6,5,8,15,11,14,
86  & 3,1,4,6,9,13,12,15/
87 !
88 ! flag for shape functions
89 !
90  data iflag /2/
91  data indexel /0/
92 !
93  save indexel
94 !
95  include "gauss.f"
96 !
97  if((iinc.eq.1).and.(iit.le.0)) call random_seed()
98 !
99 ! opening a file to store the contact spring elements
100 !
101  if(filab(1)(3:3).eq.'C') then
102  iteller=iteller+1
103 !
104  do i=1,132
105  if(jobnamef(1)(i:i).eq.' ') exit
106  enddo
107  i=i-1
108  cfile=jobnamef(1)(1:i)//'.cel'
109  open(27,file=cfile,status='unknown',position='append')
110 !
111  setname(1:15)='contactelements'
112  lenset=15
113  number(1)=istep
114  number(2)=iinc
115  number(3)=icutb+1
116  number(4)=iit
117  if(number(4).le.0) number(4)=1
118  do i=1,4
119  setname(lenset+1:lenset+1)='_'
120  if(i.eq.1) then
121  setname(lenset+2:lenset+3)='st'
122  elseif(i.eq.2) then
123  setname(lenset+2:lenset+3)='in'
124  elseif(i.eq.3) then
125  setname(lenset+2:lenset+3)='at'
126  else
127  setname(lenset+2:lenset+3)='it'
128  endif
129  if(number(i).lt.10) then
130  write(setname(lenset+4:lenset+4),'(i1)') number(i)
131  lenset=lenset+4
132  elseif(number(i).lt.100) then
133  write(setname(lenset+4:lenset+5),'(i2)') number(i)
134  lenset=lenset+5
135  elseif(number(i).lt.1000) then
136  write(setname(lenset+4:lenset+6),'(i3)') number(i)
137  lenset=lenset+6
138  else
139  write(*,*) '*ERROR in gencontelem_f2f: no more than 1000'
140  write(*,*) ' steps/increments/cutbacks/iterations'
141  write(*,*) ' allowed (for output in '
142  write(*,*) ' contactelements.inp)'
143  stop
144  endif
145  enddo
146  endif
147 !
148  iloc=0
149 !
150 ! loop over all active contact ties
151 !
152  do i=1,ntie
153  if(tieset(1,i)(81:81).ne.'C') cycle
154  imat=int(tietol(2,i))
155 !
156 ! sorting the centers of gravity of the triangulation
157 ! of the master surface in order to find the master
158 ! location corresponding to each slave integration point
159 ! (only at the start of an increment)
160 !
161  if(iit.le.0) then
162  kneigh=1
163  nstart=itietri(1,i)-1
164  n=itietri(2,i)-nstart
165  if(n.lt.kneigh) kneigh=n
166  do j=1,n
167  xo(j)=cg(1,nstart+j)
168  x(j)=xo(j)
169  nx(j)=j
170  yo(j)=cg(2,nstart+j)
171  y(j)=yo(j)
172  ny(j)=j
173  zo(j)=cg(3,nstart+j)
174  z(j)=zo(j)
175  nz(j)=j
176  enddo
177  kflag=2
178  call dsort(x,nx,n,kflag)
179  call dsort(y,ny,n,kflag)
180  call dsort(z,nz,n,kflag)
181  endif
182 !
183 ! loop over all slave faces
184 !
185 ! iprev is the number of contact elements at the end
186 ! of the previous increment at the start of a
187 ! new increment, else zero
188 ! iact is the number of actual contact elements
189 !
190  do iloop=1,2
191  iprev=0
192  iact=0
193  do jj=itiefac(1,i), itiefac(2,i)
194  ifaces=islavsurf(1,jj)
195  nelems=int(ifaces/10)
196  jfaces=ifaces-nelems*10
197 !
198  if(lakon(nelems)(4:5).eq.'8R') then
199  mint2d=1
200  nopes=4
201  nope=8
202  elseif(lakon(nelems)(4:4).eq.'8') then
203  mint2d=4
204  nopes=4
205  nope=8
206  elseif(lakon(nelems)(4:6).eq.'20R') then
207  mint2d=4
208  nopes=8
209  nope=20
210  elseif(lakon(nelems)(4:5).eq.'20') then
211  mint2d=9
212  nopes=8
213  nope=20
214  elseif(lakon(nelems)(4:5).eq.'10') then
215  mint2d=3
216  nopes=6
217  nope=10
218  elseif(lakon(nelems)(4:4).eq.'4') then
219  mint2d=1
220  nopes=3
221  nope=4
222 !
223 ! treatment of wedge faces
224 !
225  elseif(lakon(nelems)(4:4).eq.'6') then
226  mint2d=1
227  nope=6
228  if(jfaces.le.2) then
229  nopes=3
230  else
231  nopes=4
232  endif
233  elseif(lakon(nelems)(4:5).eq.'15') then
234  nope=15
235  if(jfaces.le.2) then
236  mint2d=3
237  nopes=6
238  else
239  mint2d=4
240  nopes=8
241  endif
242  endif
243 !
244 ! actual position of the nodes belonging to the
245 ! slave surface
246 !
247  if((nope.eq.20).or.(nope.eq.8)) then
248  do m=1,nopes
249  nodefs(m)=kon(ipkon(nelems)+ifaceq(m,jfaces))
250  do j=1,3
251  xl2(j,m)=co(j,nodefs(m))+clearini(j,m,jj)
252  & *reltime+vold(j,nodefs(m))
253  enddo
254  enddo
255  elseif((nope.eq.10).or.(nope.eq.4)) then
256  do m=1,nopes
257  nodefs(m)=kon(ipkon(nelems)+ifacet(m,jfaces))
258  do j=1,3
259  xl2(j,m)=co(j,nodefs(m))+clearini(j,m,jj)
260  & *reltime+vold(j,nodefs(m))
261  enddo
262  enddo
263  elseif(nope.eq.15) then
264  do m=1,nopes
265  nodefs(m)=kon(ipkon(nelems)+ifacew2(m,jfaces))
266  do j=1,3
267  xl2(j,m)=co(j,nodefs(m))+clearini(j,m,jj)
268  & *reltime+vold(j,nodefs(m))
269  enddo
270  enddo
271  else
272  do m=1,nopes
273  nodefs(m)=kon(ipkon(nelems)+ifacew1(m,jfaces))
274  do j=1,3
275  xl2(j,m)=co(j,nodefs(m))+clearini(j,m,jj)
276  & *reltime+vold(j,nodefs(m))
277  enddo
278  enddo
279  endif
280 !
281 ! loop over all integration points in the slave surface
282 !
283 ! actions which are only done in the first iteration of an
284 ! increment (for each slave face integration point)
285 !
286 ! - determine an opposite master face
287 ! - determine the local coordinates of the opposite master
288 ! location
289 ! - determine the local normal on the opposite master
290 ! location
291 ! - determine the representative slave area for the slave
292 ! integration point
293 !
294  area=0.d0
295  mint2d=islavsurf(2,jj+1)-islavsurf(2,jj)
296  if(mint2d.eq.0) cycle
297  indexf=islavsurf(2,jj)
298 !
299  do m=1,mint2d
300  iloc=indexf+m
301  xi=pslavsurf(1,indexf+m)
302  et=pslavsurf(2,indexf+m)
303  weight=pslavsurf(3,indexf+m)
304 !
305  if(nopes.eq.9) then
306  call shape9q(xi,et,xl2,xsj2,xs2,shp2,iflag)
307  elseif(nopes.eq.8) then
308  call shape8q(xi,et,xl2,xsj2,xs2,shp2,iflag)
309  elseif(nopes.eq.4) then
310  call shape4q(xi,et,xl2,xsj2,xs2,shp2,iflag)
311  elseif(nopes.eq.6) then
312  call shape6tri(xi,et,xl2,xsj2,xs2,shp2,iflag)
313  elseif(nopes.eq.7) then
314  call shape7tri(xi,et,xl2,xsj2,xs2,shp2,iflag)
315  else
316  call shape3tri(xi,et,xl2,xsj2,xs2,shp2,iflag)
317  endif
318  if(iit.le.0) then
319  area=dsqrt(xsj2(1)**2+xsj2(2)**2+xsj2(3)**2)*weight
320  endif
321 !
322 ! search a master face for each gauss point and generate a contact
323 ! spring element if successful
324 !
325  do k=1,3
326  p(k)=0.d0
327  do j=1,nopes
328  p(k)=p(k)+shp2(4,j)*xl2(k,j)
329  enddo
330  enddo
331  if(iit.gt.0) then
332  if(int(pmastsurf(3,iloc)).eq.0) then
333  isol=0
334  else
335  isol=1
336  endif
337  else
338 !
339 ! determining the kneigh neighboring master contact
340 ! triangle centers of gravity
341 !
342  call near3d(xo,yo,zo,x,y,z,nx,ny,nz,p(1),p(2),p(3),
343  & n,neigh,kneigh)
344 !
345  isol=0
346 !
347  itriold=0
348  itri=neigh(1)+itietri(1,i)-1
349  ntriangle=0
350  ntriangle_=100
351 !
352  loop1: do
353  do l=1,3
354  ll=4*l-3
355  dist=straight(ll,itri)*p(1)+
356  & straight(ll+1,itri)*p(2)+
357  & straight(ll+2,itri)*p(3)+
358  & straight(ll+3,itri)
359 !
360 ! 1.d-6 was increased to 1.d-3 on 19/04/2012
361 ! this is important for 2d-calculations or
362 ! calculations for which structures fit exactly
363 ! at their boundaries
364 !
365  if(dist.gt.1.d-3*dsqrt(area)) then
366  itrinew=imastop(l,itri)
367  if(itrinew.eq.0) then
368 c write(*,*) '**border reached'
369  exit loop1
370  elseif((itrinew.lt.itietri(1,i)).or.
371  & (itrinew.gt.itietri(2,i))) then
372 c write(*,*) '**border reached'
373  exit loop1
374  elseif(itrinew.eq.itriold) then
375 c write(*,*) '**solution in between triangles'
376  isol=itri
377  exit loop1
378  else
379  call nident(itriangle,itrinew,
380  & ntriangle,id)
381  if(id.gt.0) then
382  if(itriangle(id).eq.itrinew) then
383 c write(*,*) '**circular path;no solution'
384  exit loop1
385  endif
386  endif
387  ntriangle=ntriangle+1
388  if(ntriangle.gt.ntriangle_) then
389 c write(*,*) '**too many iterations'
390  exit loop1
391  endif
392  do k=ntriangle,id+2,-1
393  itriangle(k)=itriangle(k-1)
394  enddo
395  itriangle(id+1)=itrinew
396  itriold=itri
397  itri=itrinew
398  cycle loop1
399  endif
400  elseif(l.eq.3) then
401 c write(*,*) '**regular solution'
402  isol=itri
403  exit loop1
404  endif
405  enddo
406  enddo loop1
407  endif
408 !
409 ! integration point is catalogued if opposite master
410 ! face was detected
411 !
412  if(isol.eq.0) then
413  if(iit.le.0) then
414  pmastsurf(3,iloc)=0.5d0
415  endif
416  else
417 !
418 ! determining the clearance
419 !
420 ! identifying the element face to which the
421 ! triangle belongs
422 !
423  if(iit.le.0) then
424  springarea(1,iloc)=area
425  springarea(2,iloc)=0.d0
426  ifacem=koncont(4,itri)
427  pmastsurf(3,iloc)=ifacem+0.5d0
428  else
429  ifacem=int(pmastsurf(3,iloc))
430  endif
431  nelemm=int(ifacem/10.d0)
432  jfacem=ifacem-10*nelemm
433 !
434  indexe=ipkon(nelemm)
435  if(lakon(nelemm)(4:5).eq.'20') then
436  nopem=8
437  nface=6
438  elseif(lakon(nelemm)(4:4).eq.'8') then
439  nopem=4
440  nface=6
441  elseif(lakon(nelemm)(4:5).eq.'10') then
442  nopem=6
443  nface=4
444  elseif(lakon(nelemm)(4:4).eq.'4') then
445  nopem=3
446  nface=4
447  elseif(lakon(nelemm)(4:5).eq.'15') then
448  if(jfacem.le.2) then
449  nopem=6
450  else
451  nopem=8
452  endif
453  nface=5
454  nope=15
455  elseif(lakon(nelemm)(4:4).eq.'6') then
456  if(jfacem.le.2) then
457  nopem=3
458  else
459  nopem=4
460  endif
461  nface=5
462  nope=6
463  else
464  cycle
465  endif
466 !
467 ! determining the nodes of the master face
468 !
469  if(nface.eq.4) then
470  do k=1,nopem
471  nodefm(k)=kon(indexe+ifacet(k,jfacem))
472  enddo
473  elseif(nface.eq.5) then
474  if(nope.eq.6) then
475  do k=1,nopem
476  nodefm(k)=kon(indexe+ifacew1(k,jfacem))
477  enddo
478  elseif(nope.eq.15) then
479  do k=1,nopem
480  nodefm(k)=kon(indexe+ifacew2(k,jfacem))
481  enddo
482  endif
483  elseif(nface.eq.6) then
484  do k=1,nopem
485  nodefm(k)=kon(indexe+ifaceq(k,jfacem))
486  enddo
487  endif
488 !
489 ! total number of nodes belonging to the contact
490 ! element
491 !
492  nope=nopem+nopes
493 !
494 ! orthogonal projection on the master element face
495 !
496  do k=1,nopem
497  do nn=1,3
498  pl(nn,k)=co(nn,nodefm(k))+vold(nn,nodefm(k))
499  enddo
500  enddo
501 !
502  if(iit.le.0) then
503  do nn=1,3
504  pproj(nn)=p(nn)
505  enddo
506  call attach(pl,pproj,nopem,ratio,dist,xi,et)
507  pmastsurf(1,indexf+m)=xi
508  pmastsurf(2,indexf+m)=et
509  else
510  xi=pmastsurf(1,indexf+m)
511  et=pmastsurf(2,indexf+m)
512  endif
513 !
514  if(nopem.eq.9) then
515  call shape9q(xi,et,pl,xm,xs2,shp2,iflag)
516  elseif(nopem.eq.8) then
517  call shape8q(xi,et,pl,xm,xs2,shp2,iflag)
518  elseif(nopem.eq.4) then
519  call shape4q(xi,et,pl,xm,xs2,shp2,iflag)
520  elseif(nopem.eq.6) then
521  call shape6tri(xi,et,pl,xm,xs2,shp2,iflag)
522  elseif(nopem.eq.7) then
523  call shape7tri(xi,et,pl,xm,xs2,shp2,iflag)
524  else
525  call shape3tri(xi,et,pl,xm,xs2,shp2,iflag)
526  endif
527 !
528  if(iit.gt.0) then
529  do nn=1,3
530  pproj(nn)=0.d0
531  do k=1,nopem
532  pproj(nn)=pproj(nn)+shp2(4,k)*pl(nn,k)
533  enddo
534  enddo
535  endif
536 !
537  do nn=1,3
538  al(nn)=p(nn)-pproj(nn)
539  enddo
540 !
541 ! normal on the surface
542 !
543  if(iit.le.0) then
544  dm=dsqrt(xm(1)*xm(1)+xm(2)*xm(2)+xm(3)*xm(3))
545  do nn=1,3
546  xn(nn)=xm(nn)/dm
547  enddo
548  pmastsurf(4,iloc)=xn(1)
549  pmastsurf(5,iloc)=xn(2)
550  pmastsurf(6,iloc)=xn(3)
551  else
552  xn(1)=pmastsurf(4,iloc)
553  xn(2)=pmastsurf(5,iloc)
554  xn(3)=pmastsurf(6,iloc)
555  endif
556 !
557 ! distance from surface along normal (= clearance)
558 !
559  clear=al(1)*xn(1)+al(2)*xn(2)+al(3)*xn(3)
560 !
561  if((nmethod.eq.4).or.(nmethod.eq.2)) then
562 c if(nmethod.eq.4) then
563 !
564 ! dynamic calculation
565 !
566  if((clear.gt.0.d0).and.
567  & (int(elcon(3,1,imat)).ne.4)) then
568  isol=0
569  endif
570  else
571 !
572 ! static calculation: more sophisticated
573 ! algorithm to detect contact
574 !
575  if((istep.eq.1).and.(iit.le.0.d0)) then
576  if(clear.lt.0.d0) then
577  springarea(2,iloc)=clear/(1.d0-theta)
578  elseif(clear.lt.1.d0/elcon(2,1,imat)) then
579  clear=0.d0
580  endif
581  endif
582  clear=clear-springarea(2,iloc)*(1.d0-reltime)
583 !
584 ! if iloop=1 AND a new increment was started the
585 ! number of contact elements at the end of the
586 ! previous increment is counted (using xstateini)
587 !
588 ! furthermore, the regular penetration criterion is
589 ! used to decide on the creation of a contact element
590 !
591  if(iloop.eq.1) then
592  if((isol.ne.0).and.(int(elcon(3,1,imat)).ne.4))
593  & then
594  if(((istep.gt.1).or.(iinc.gt.1)).and.
595  & (iit.le.0).and.(ncmat_.ge.7).and.
596  & (elcon(6,1,imat).gt.0.d0)) then
597  if(dsqrt(xstateini(4,1,ne0+iloc)**2+
598  & xstateini(5,1,ne0+iloc)**2+
599  & xstateini(6,1,ne0+iloc)**2)
600  & .ge.1.d-30) iprev=iprev+1
601  endif
602 !
603 ! if a local minimum was detected in checkconvergence:
604 ! remove in an aleatoric way 10 % of the contact elements
605 !
606  if(ialeatoric.eq.1) then
607  call random_number(harvest)
608  if(harvest.gt.0.9) isol=0
609  endif
610 !
611  if(isol.ne.0) then
612  if((icutb.eq.0).or.(ncmat_.lt.7).or.
613  & (elcon(6,1,imat).le.0.)) then
614 !
615 ! this is no cut-back: only use a negative
616 ! clearance as contact criterion
617 ! (this also applies if no friction is
618 ! defined)
619 !
620  if(clear.gt.0.d0) then
621  isol=0
622  else
623  iact=iact+1
624  endif
625  else
626 !
627 ! this is the first iteration in a cut-back:
628 ! all contact elements from the end of the
629 ! previous increment are included as well
630 !
631  if((dsqrt(xstateini(4,1,ne0+iloc)**2+
632  & xstateini(5,1,ne0+iloc)**2+
633  & xstateini(6,1,ne0+iloc)**2)
634  & .lt.1.d-30).and.(clear.gt.0.d0))
635  & then
636  isol=0
637  else
638  iact=iact+1
639  endif
640  endif
641  endif
642  endif
643  else
644 !
645 ! if iloop=2 it means that at the start of a new
646 ! increment no contact elements were generated
647 ! based on penetration, although there were contact
648 ! elements at the end of the last increment. In that
649 ! case the contact elements from the end of the last
650 ! increment are taken.
651 !
652  if((isol.ne.0).and.(int(elcon(3,1,imat)).ne.4))
653  & then
654  if(dsqrt(xstateini(4,1,ne0+iloc)**2+
655  & xstateini(5,1,ne0+iloc)**2+
656  & xstateini(6,1,ne0+iloc)**2)
657  & .lt.1.d-30) isol=0
658  endif
659  endif
660 !
661  endif
662 !
663  endif
664 !
665  if(isol.ne.0) then
666 !
667 ! generation of a contact spring element
668 !
669  ne=ne+1
670  ipkon(ne)=ifree+1
671  lakon(ne)='ESPRNGC '
672  ielmat(1,ne)=imat
673 !
674 ! nasym indicates whether at least one contact
675 ! spring elements exhibits friction in the present
676 ! step. If so, nasym=1, else nasym=0; nasym=1
677 ! triggers the asymmetric equation solver
678 !
679  if(ncmat_.ge.7) then
680  if((elcon(6,1,imat).gt.0).and.
681  & (int(elcon(3,1,imat)).ne.4)) then
682  nasym=1
683  endif
684  endif
685  if(ncmat_.ge.8) then
686  if(elcon(8,1,imat).gt.0) then
687  nasym=1
688  endif
689  endif
690 !
691  kon(ifree+1)=nopes+nopem
692  ifree=ifree+1
693 !
694  do k=1,nopem
695  kon(ifree+k)=nodefm(k)
696  enddo
697  ifree=ifree+nopem
698  do k=1,nopes
699  kon(ifree+k)=nodefs(k)
700  enddo
701  ifree=ifree+nopes
702  kon(ifree+1)=iloc
703  kon(ifree+2)=jj
704  kon(ifree+3)=indexf+m
705 !
706  ifree=ifree+3
707 !
708  write(lakon(ne)(8:8),'(i1)') nopem
709 !
710  indexel=indexel+1
711 !
712  if((nopem.eq.4).or.(nopem.eq.8)) then
713  if((nopes.eq.4).or.(nopes.eq.8)) then
714  if(filab(1)(3:3).eq.'C') then
715  write(27,100) setname(1:lenset)
716  100 format('*ELEMENT,TYPE=C3D8,ELSET=',a)
717  write(27,*) ne0+indexel,',',nodefm(1),',',
718  & nodefm(2),',',nodefm(3),',',nodefm(4),
719  & ',',nodefs(2),',',nodefs(1),',',
720  & nodefs(4),',',nodefs(3)
721  endif
722  endif
723  if((nopes.eq.3).or.(nopes.eq.6)) then
724  if(filab(1)(3:3).eq.'C') then
725  write(27,101) setname(1:lenset)
726  101 format('*ELEMENT,TYPE=C3D8,ELSET=',a)
727  write(27,*) ne0+indexel,',',nodefm(1),',',
728  & nodefm(2),',',nodefm(3),',',nodefm(4),
729  & ',',nodefs(2),',',nodefs(1),',',
730  & nodefs(3),',',nodefs(3)
731  endif
732  endif
733  endif
734  if((nopem.eq.3).or.(nopem.eq.6)) then
735  if((nopes.eq.4).or.(nopes.eq.8)) then
736  if(filab(1)(3:3).eq.'C') then
737  write(27,102) setname(1:lenset)
738  102 format('*ELEMENT,TYPE=C3D8,ELSET=',a)
739  write(27,*) ne0+indexel,',',nodefm(1),',',
740  & nodefm(2),',',nodefm(3),',',nodefm(3),
741  & ',',nodefs(2),',',nodefs(1),',',
742  & nodefs(4),',',nodefs(3)
743  endif
744  endif
745  if((nopes.eq.3).or.(nopes.eq.6)) then
746  if(filab(1)(3:3).eq.'C') then
747  write(27,103) setname(1:lenset)
748  103 format('*ELEMENT,TYPE=C3D6,ELSET=',a)
749  write(27,*) ne0+indexel,',',nodefm(1),',',
750  & nodefm(2),',',nodefm(3),',',nodefs(2),
751  & ',',nodefs(1),',',nodefs(3)
752  endif
753  endif
754  endif
755  else
756 !
757 ! no contact: set internal variables at end of increment
758 ! to zero
759 !
760 ! next line is inserted because xstate is not yet
761 ! allocated for iit <= 0
762 !
763  if(iit.gt.0) then
764  do k=1,nstate_
765  xstate(k,1,ne0+iloc)=0.d0
766  enddo
767  endif
768 !
769  endif
770 !
771  enddo
772  enddo
773  if((iact.ne.0).or.(iprev.eq.0).or.(nmethod.eq.4)) exit
774  write(*,*)'*INFO in gencontelem_f2f: contact lost at the'
775  write(*,*)' start of a new increment; contact'
776  write(*,*)' elements from end of previous increment'
777  write(*,*)' are kept.'
778  write(*,*)' number of previous contact elements:',iprev
779  write(*,*)' number of actual contact elements:',iact
780  enddo
781  enddo
782 !
783 ! closing the file containing the contact elements
784 !
785  if(filab(1)(3:3).eq.'C') then
786  close(27)
787  endif
788 !
789  return
near3d
subroutine near3d(xo, yo, zo, x, y, z, nx, ny, nz, xp, yp, zp, n, neighbor, k)
Definition: near3d.f:20
shape9q
subroutine shape9q(xi, et, xl, xsj, xs, shp, iflag)
Definition: shape9q.f:20
shape8q
subroutine shape8q(xi, et, xl, xsj, xs, shp, iflag)
Definition: shape8q.f:20
shape3tri
subroutine shape3tri(xi, et, xl, xsj, xs, shp, iflag)
Definition: shape3tri.f:20
shape7tri
subroutine shape7tri(xi, et, xl, xsj, xs, shp, iflag)
Definition: shape7tri.f:20
stop
subroutine stop()
Definition: stop.f:20
dist
static double * dist
Definition: radflowload.c:42
shape4q
subroutine shape4q(xi, et, xl, xsj, xs, shp, iflag)
Definition: shape4q.f:20
nident
subroutine nident(x, px, n, id)
Definition: nident.f:26
dsort
subroutine dsort(dx, iy, n, kflag)
Definition: dsort.f:6
shape6tri
subroutine shape6tri(xi, et, xl, xsj, xs, shp, iflag)
Definition: shape6tri.f:20
attach
subroutine attach(pneigh, pnode, nterms, ratio, dist, xil, etl)
Definition: attach.f:20
Hosted by OpenAircraft.com, (Michigan UAV, LLC)