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

Go to the source code of this file.

Functions/Subroutines

subroutine updatecontpen (koncont, ncont, co, vold, cg, straight, mi, imastnode, nmastnode, xmastnor, ntie, tieset, nset, set, istartset, iendset, ialset, ipkon, lakon, kon, cs, mcs, ics)
 

Function/Subroutine Documentation

◆ updatecontpen()

subroutine updatecontpen ( integer, dimension(4,*)  koncont,
integer  ncont,
real*8, dimension(3,*)  co,
real*8, dimension(0:mi(2),*)  vold,
real*8, dimension(3,*)  cg,
real*8, dimension(16,*)  straight,
integer, dimension(*)  mi,
integer, dimension(*)  imastnode,
integer, dimension(*)  nmastnode,
real*8, dimension(3,*)  xmastnor,
integer  ntie,
character*81, dimension(3,*)  tieset,
integer  nset,
character*81, dimension(*)  set,
integer, dimension(*)  istartset,
integer, dimension(*)  iendset,
integer, dimension(*)  ialset,
integer, dimension(*)  ipkon,
character*8, dimension(*)  lakon,
integer, dimension(*)  kon,
real*8, dimension(17,*)  cs,
integer  mcs,
integer, dimension(*)  ics 
)
23 !
24 ! update geometric date of the contact master surface triangulation
25 !
26  implicit none
27 !
28  character*1 kind1,kind2
29  character*8 lakon(*)
30  character*81 tieset(3,*),mastset,set(*)
31 !
32  integer koncont(4,*),ncont,i,j,k,node,mi(*),imastnode(*),
33  & nmastnode(*),ntie,imast,iplaneaxial,noeq,
34  & istartset(*),iendset(*),ialset(*),ifacem,nelemm,
35  & jfacem,indexe,ipkon(*),nopem,nope,konl(26),kon(*),m,
36  & ifaceq(8,6),ifacet(6,4),ifacew1(4,5),
37  & ifacew2(8,5),id,index1,indexnode(9),l,iflag,nset,itriact,
38  & ipos,ntrifaces,noeq4(2),noeq8(6),mcs,ics(*),icyc(3),
39  & istart,ilength,noeq9(8)
40 !
41  real*8 co(3,*),vold(0:mi(2),*),cg(3,*),straight(16,*),col(3,3),
42  & xmastnor(3,*),xl2m(3,9),xi,et,dd,xsj2(3),shp2(7,9),
43  & xs2(3,2),xnor(3,3),cs(17,*),xquad(2,9),xtri(2,7)
44 !
45 ! nodes per face for hex elements
46 !
47  data ifaceq /4,3,2,1,11,10,9,12,
48  & 5,6,7,8,13,14,15,16,
49  & 1,2,6,5,9,18,13,17,
50  & 2,3,7,6,10,19,14,18,
51  & 3,4,8,7,11,20,15,19,
52  & 4,1,5,8,12,17,16,20/
53 !
54 ! nodes per face for tet elements
55 !
56  data ifacet /1,3,2,7,6,5,
57  & 1,2,4,5,9,8,
58  & 2,3,4,6,10,9,
59  & 1,4,3,8,10,7/
60 !
61 ! nodes per face for linear wedge elements
62 !
63  data ifacew1 /1,3,2,0,
64  & 4,5,6,0,
65  & 1,2,5,4,
66  & 2,3,6,5,
67  & 3,1,4,6/
68 !
69 ! nodes per face for quadratic wedge elements
70 !
71  data ifacew2 /1,3,2,9,8,7,0,0,
72  & 4,5,6,10,11,12,0,0,
73  & 1,2,5,4,7,14,10,13,
74  & 2,3,6,5,8,15,11,14,
75  & 3,1,4,6,9,13,12,15/
76 !
77  data iflag /2/
78 !
79 ! new added data for the local coodinates for nodes
80 !
81  data xquad /-1.d0,-1.d0,
82  & 1.d0,-1.d0,
83  & 1.d0,1.d0,
84  & -1.d0,1.d0,
85  & 0.d0,-1.d0,
86  & 1.d0,0.d0,
87  & 0.d0,1.d0,
88  & -1.d0,0.d0,
89  & 0.d0,0.d0/
90 !
91  data xtri /0.d0,0.d0,
92  & 1.d0,0.d0,
93  & 0.d0,1.d0,
94  & .5d0,0.d0,
95  & .5d0,.5d0,
96  & 0.d0,.5d0,
97  & 0.333333333333333d0,0.333333333333333d0/
98 !
99 ! relative node position in the triangle topology the side
100 ! opposite of which does not need to be checked for the
101 ! slave node/master triangle matching (because of expansion
102 ! of plane stress/strain/axisymmetric elements)
103 !
104  data noeq4 /3,1/
105  data noeq8 /3,3,2,1,0,0/
106  data noeq9 /3,3,3,3,3,3,3,3/
107 !
108  kind1='C'
109  kind2='-'
110 !
111 ! calculation of the normals on the master surface
112 ! at the vertex nodes of the master surface -> xmastnor
113 ! Iteration over all tiesets
114 !
115  do i=1,ntie
116 !
117 ! check for contact conditions
118 !
119  if((tieset(1,i)(81:81).eq.'C').or.
120  & (tieset(1,i)(81:81).eq.'-')) then
121 
122  mastset=tieset(3,i)
123  do j=1,nset
124  if(set(j).eq.mastset) exit
125  enddo
126  if(j.gt.nset) then
127  write(*,*) '*ERROR in tiefaccont: master surface'
128  write(*,*) ' does not exist'
129  call exit(201)
130  endif
131  imast=j
132 !
133 ! calculate the mean normal in the master nodes
134 !
135  do j=istartset(imast),iendset(imast)
136  if(ialset(j).gt.0) then
137 !
138  ifacem=ialset(j)
139  nelemm=int(ifacem/10)
140  jfacem=ifacem - nelemm*10
141  indexe=ipkon(nelemm)
142 !
143 ! nopem: # of nodes in the master face
144 ! nope: # of nodes in the element
145 !
146  if(lakon(nelemm)(4:5).eq.'8R') then
147  nopem=4
148  nope=8
149  elseif(lakon(nelemm)(4:4).eq.'8') then
150  nopem=4
151  nope=8
152  elseif(lakon(nelemm)(4:5).eq.'20') then
153  nopem=8
154  nope=20
155  elseif(lakon(nelemm)(4:5).eq.'10') then
156  nopem=6
157  nope=10
158  elseif(lakon(nelemm)(4:4).eq.'4') then
159  nopem=3
160  nope=4
161 !
162 ! treatment of wedge faces
163 !
164  elseif(lakon(nelemm)(4:4).eq.'6') then
165  nope=6
166  if(jfacem.le.2) then
167  nopem=3
168  else
169  nopem=4
170  endif
171  elseif(lakon(nelemm)(4:5).eq.'15') then
172  nope=15
173  if(jfacem.le.2) then
174  nopem=6
175  else
176  nopem=8
177  endif
178  endif
179 !
180 ! actual position of the nodes belonging to the
181 ! master surface
182 !
183  do k=1,nope
184  konl(k)=kon(ipkon(nelemm)+k)
185  enddo
186 !
187  if((nope.eq.20).or.(nope.eq.8)) then
188  do m=1,nopem
189  do k=1,3
190  xl2m(k,m)=co(k,konl(ifaceq(m,jfacem)))+
191  & vold(k,konl(ifaceq(m,jfacem)))
192  enddo
193  enddo
194  elseif((nope.eq.10).or.(nope.eq.4))
195  & then
196  do m=1,nopem
197  do k=1,3
198  xl2m(k,m)=co(k,konl(ifacet(m,jfacem)))+
199  & vold(k,konl(ifacet(m,jfacem)))
200  enddo
201  enddo
202  elseif(nope.eq.15) then
203  do m=1,nopem
204  do k=1,3
205  xl2m(k,m)=co(k,konl(ifacew2(m,jfacem)))+
206  & vold(k,konl(ifacew2(m,jfacem)))
207  enddo
208  enddo
209  else
210  do m=1,nopem
211  do k=1,3
212  xl2m(k,m)=co(k,konl(ifacew1(m,jfacem)))+
213  & vold(k,konl(ifacew1(m,jfacem)))
214  enddo
215  enddo
216  endif
217 
218 ! calculate the normal vector in the nodes belonging to the master surface
219 !
220  if(nopem.eq.8) then
221  do m=1,nopem
222  xi=xquad(1,m)
223  et=xquad(2,m)
224  call shape8q(xi,et,xl2m,xsj2,xs2,shp2,iflag)
225  dd=dsqrt(xsj2(1)*xsj2(1) + xsj2(2)*xsj2(2)
226  & + xsj2(3)*xsj2(3))
227  xsj2(1)=xsj2(1)/dd
228  xsj2(2)=xsj2(2)/dd
229  xsj2(3)=xsj2(3)/dd
230 !
231  if(nope.eq.20) then
232  node=konl(ifaceq(m,jfacem))
233  elseif(nope.eq.15) then
234  node=konl(ifacew2(m,jfacem))
235  endif
236 !
237  call nident(imastnode(nmastnode(i)+1),node,
238  & nmastnode(i+1)-nmastnode(i),id)
239  index1=nmastnode(i)+id
240  indexnode(m)=index1
241  xmastnor(1,index1)=xmastnor(1,index1)
242  & +xsj2(1)
243  xmastnor(2,index1)=xmastnor(2,index1)
244  & +xsj2(2)
245  xmastnor(3,index1)=xmastnor(3,index1)
246  & +xsj2(3)
247  enddo
248  elseif(nopem.eq.4) then
249  do m=1,nopem
250  xi=xquad(1,m)
251  et=xquad(2,m)
252  call shape4q(xi,et,xl2m,xsj2,xs2,shp2,iflag)
253  dd=dsqrt(xsj2(1)*xsj2(1) + xsj2(2)*xsj2(2)
254  & + xsj2(3)*xsj2(3))
255  xsj2(1)=xsj2(1)/dd
256  xsj2(2)=xsj2(2)/dd
257  xsj2(3)=xsj2(3)/dd
258 !
259  if(nope.eq.8) then
260  node=konl(ifaceq(m,jfacem))
261  elseif(nope.eq.6) then
262  node=konl(ifacew1(m,jfacem))
263  endif
264 !
265  call nident(imastnode(nmastnode(i)+1),node,
266  & nmastnode(i+1)-nmastnode(i),id)
267 !
268  index1=nmastnode(i)+id
269  indexnode(m)=index1
270  xmastnor(1,index1)=xmastnor(1,index1)
271  & +xsj2(1)
272  xmastnor(2,index1)=xmastnor(2,index1)
273  & +xsj2(2)
274  xmastnor(3,index1)=xmastnor(3,index1)
275  & +xsj2(3)
276  enddo
277  elseif(nopem.eq.6) then
278  do m=1,nopem
279  xi=xtri(1,m)
280  et=xtri(2,m)
281  call shape6tri(xi,et,xl2m,xsj2,xs2,shp2,iflag)
282  dd=dsqrt(xsj2(1)*xsj2(1) + xsj2(2)*xsj2(2)
283  & + xsj2(3)*xsj2(3))
284  xsj2(1)=xsj2(1)/dd
285  xsj2(2)=xsj2(2)/dd
286  xsj2(3)=xsj2(3)/dd
287 !
288  if(nope.eq.10) then
289  node=konl(ifacet(m,jfacem))
290  elseif(nope.eq.15) then
291  node=konl(ifacew2(m,jfacem))
292  endif
293 !
294  call nident(imastnode(nmastnode(i)+1),node,
295  & nmastnode(i+1)-nmastnode(i),id)
296  index1=nmastnode(i)+id
297  indexnode(m)=index1
298  xmastnor(1,index1)=xmastnor(1,index1)
299  & +xsj2(1)
300  xmastnor(2,index1)=xmastnor(2,index1)
301  & +xsj2(2)
302  xmastnor(3,index1)=xmastnor(3,index1)
303  & +xsj2(3)
304  enddo
305  else
306  do m=1,nopem
307  xi=xtri(1,m)
308  et=xtri(2,m)
309  call shape3tri(xi,et,xl2m,xsj2,xs2,shp2,iflag)
310  dd=dsqrt(xsj2(1)*xsj2(1) + xsj2(2)*xsj2(2)
311  & + xsj2(3)*xsj2(3))
312  xsj2(1)=xsj2(1)/dd
313  xsj2(2)=xsj2(2)/dd
314  xsj2(3)=xsj2(3)/dd
315 !
316  if(nope.eq.6) then
317  node=konl(ifacew1(m,jfacem))
318  elseif(nope.eq.4) then
319  node=konl(ifacet(m,jfacem))
320  endif
321 !
322  call nident(imastnode(nmastnode(i)+1),node,
323  & nmastnode(i+1)-nmastnode(i),id)
324  index1=nmastnode(i)+id
325  indexnode(m)=index1
326  xmastnor(1,nmastnode(i)+id)=xmastnor(1,index1)
327  & +xsj2(1)
328  xmastnor(2,nmastnode(i)+id)=xmastnor(2,index1)
329  & +xsj2(2)
330  xmastnor(3,nmastnode(i)+id)=xmastnor(3,index1)
331  & +xsj2(3)
332  enddo
333  endif
334  endif
335  enddo
336 !
337 ! normalizing the normals
338 !
339  do l=nmastnode(i)+1,nmastnode(i+1)
340  dd=dsqrt(xmastnor(1,l)**2+xmastnor(2,l)**2+
341  & xmastnor(3,l)**2)
342  do m=1,3
343  xmastnor(m,l)=xmastnor(m,l)/dd
344  enddo
345  enddo
346  endif
347  enddo
348 !
349 ! Loop over all tiesets
350 !
351  itriact=0
352  do i=1,ntie
353 !
354 ! check for contact conditions
355 !
356  if((tieset(1,i)(81:81).eq.kind1).or.
357  & (tieset(1,i)(81:81).eq.kind2)) then
358  mastset=tieset(3,i)
359 !
360 ! determining the master surface
361 !
362  do j=1,nset
363  if(set(j).eq.mastset) exit
364  enddo
365  if(j.gt.nset) then
366  ipos=index(mastset,' ')
367  write(*,*) '*ERROR in updatecont: master surface',
368  & mastset(1:ipos-2)
369  write(*,*) ' does not exist or does not contain'
370  write(*,*) ' element faces'
371  call exit(201)
372  endif
373  imast=j
374 
375  do j=istartset(imast),iendset(imast)
376  nelemm=int(ialset(j)/10.d0)
377  jfacem=ialset(j)-10*nelemm
378 !
379  if(lakon(nelemm)(4:5).eq.'20') then
380  ntrifaces=6
381  elseif(lakon(nelemm)(4:4).eq.'8') then
382  ntrifaces=2
383  elseif(lakon(nelemm)(4:5).eq.'10') then
384  ntrifaces=4
385  elseif(lakon(nelemm)(4:4).eq.'4') then
386  ntrifaces=1
387  elseif(lakon(nelemm)(4:5).eq.'15') then
388  if(jfacem.le.2) then
389  ntrifaces=4
390  else
391  ntrifaces=6
392  endif
393  elseif(lakon(nelemm)(4:4).eq.'6') then
394  if(jfacem.le.2) then
395  ntrifaces=1
396  else
397  ntrifaces=2
398  endif
399  endif
400 !
401 ! check whether plane stress, plane strain or
402 ! axisymmetric element (linear or quadratic, not
403 ! relevant for tets)
404 !
405  iplaneaxial=0
406  if((lakon(nelemm)(7:7).eq.'S').or.
407  & (lakon(nelemm)(7:7).eq.'E').or.
408  & (lakon(nelemm)(7:7).eq.'A')) then
409  if(ntrifaces.eq.2) then
410  iplaneaxial=4
411  elseif(ntrifaces.eq.6) then
412  iplaneaxial=8
413  elseif(ntrifaces.eq.8) then
414  iplaneaxial=9
415  endif
416 c if((lakon(nelemm)(4:4).eq.'2').or.
417 c & (lakon(nelemm)(4:5).eq.'15')) then
418 c iplaneaxial=8
419 c else
420 c iplaneaxial=4
421 c endif
422  endif
423 !
424 ! loop over the master triangles
425 !
426  do k=1,ntrifaces
427 !
428  itriact=itriact+1
429 !
430 ! check whether one of the triangle sides does not
431 ! need an equation. This applies to plane stress,
432 ! plane strain and axisymmetric elements on the sides
433 ! with constant z-coordinate or constant angle
434 !
435  if(iplaneaxial.eq.0) then
436  noeq=0
437  elseif(iplaneaxial.eq.4) then
438  noeq=noeq4(k)
439  elseif(iplaneaxial.eq.8) then
440  noeq=noeq8(k)
441  elseif(iplaneaxial.eq.9) then
442  noeq=noeq9(k)
443  endif
444 !
445 ! check for cyclic symmetric structures; triangle
446 ! edges in a cyclic symmetry plane get a zero
447 ! equation
448 !
449  if(mcs.gt.0) then
450  do l=1,3
451  node=koncont(l,itriact)
452  icyc(l)=0
453  do m=1,mcs
454  istart=int(cs(14,m))+1
455  ilength=int(cs(4,m))
456  call nident(ics(istart),node,ilength,id)
457  if(id.gt.0) then
458  if(ics(istart+id-1).eq.node) then
459  icyc(l)=m
460  exit
461  endif
462  endif
463  enddo
464  enddo
465 !
466  if((icyc(1).eq.0).and.
467  & ((icyc(2).ne.0).and.(icyc(2).eq.icyc(3))))
468  & then
469  noeq=1
470  elseif((icyc(2).eq.0).and.
471  & ((icyc(3).ne.0).and.(icyc(3).eq.icyc(1))))
472  & then
473  noeq=2
474  elseif((icyc(3).eq.0).and.
475  & ((icyc(1).ne.0).and.(icyc(1).eq.icyc(2))))
476  & then
477  noeq=3
478  endif
479  endif
480 !
481 c itriact=itriact+1
482  do l=1,3
483  node=koncont(l,itriact)
484  do m=1,3
485  col(m,l)=co(m,node)+vold(m,node)
486  enddo
487  call nident(imastnode(nmastnode(i)+1),node,
488  & nmastnode(i+1)-nmastnode(i),id)
489  index1=nmastnode(i)+id
490  do m=1,3
491  xnor(m,l)=xmastnor(m,index1)
492  enddo
493  enddo
494 !
495 ! center of gravity of the triangles
496 !
497  do l=1,3
498  cg(l,itriact)=col(l,1)
499  enddo
500  do m=2,3
501  do l=1,3
502  cg(l,itriact)=cg(l,itriact)+col(l,m)
503  enddo
504  enddo
505  do l=1,3
506  cg(l,itriact)=cg(l,itriact)/3.d0
507  enddo
508 !
509 ! calculating the equation of the triangle plane and the planes
510 ! through the edges of a triangle and perpendicular to the
511 ! representative normal of each triangle edge
512 !
513  call straighteq3dpen(col,straight(1,itriact),xnor,
514  & noeq)
515  enddo
516  enddo
517  endif
518  enddo
519 !
520  return
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
straighteq3dpen
subroutine straighteq3dpen(col, straight, xnor, noeq)
Definition: straighteq3dpen.f:20
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
shape6tri
subroutine shape6tri(xi, et, xl, xsj, xs, shp, iflag)
Definition: shape6tri.f:20
Hosted by OpenAircraft.com, (Michigan UAV, LLC)