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

Go to the source code of this file.

Functions/Subroutines

subroutine pretensionsections (inpc, textpart, ipompc, nodempc, coefmpc, nmpc, nmpc_, mpcfree, nk, ikmpc, ilmpc, labmpc, istep, istat, n, iline, ipol, inl, ipoinp, inp, ipoinpc, lakon, kon, ipkon, set, nset, istartset, iendset, ialset, co, ics, dcs, t0, ithermal, ne)
 

Function/Subroutine Documentation

◆ pretensionsections()

subroutine pretensionsections ( character*1, dimension(*)  inpc,
character*132, dimension(16)  textpart,
integer, dimension(*)  ipompc,
integer, dimension(3,*)  nodempc,
real*8, dimension(*)  coefmpc,
integer  nmpc,
integer  nmpc_,
integer  mpcfree,
integer  nk,
integer, dimension(*)  ikmpc,
integer, dimension(*)  ilmpc,
character*20, dimension(*)  labmpc,
integer  istep,
integer  istat,
integer  n,
integer  iline,
integer  ipol,
integer  inl,
integer, dimension(2,*)  ipoinp,
integer, dimension(3,*)  inp,
integer, dimension(0:*)  ipoinpc,
character*8, dimension(*)  lakon,
integer, dimension(*)  kon,
integer, dimension(*)  ipkon,
character*81, dimension(*)  set,
integer  nset,
integer, dimension(*)  istartset,
integer, dimension(*)  iendset,
integer, dimension(*)  ialset,
real*8, dimension(3,*)  co,
integer, dimension(2,*)  ics,
real*8, dimension(*)  dcs,
real*8, dimension(*)  t0,
integer, dimension(2)  ithermal,
integer  ne 
)
24 !
25 ! reading the input deck: *PRE-TENSION SECTION
26 !
27  implicit none
28 !
29  logical twod
30 !
31  character*1 inpc(*)
32  character*8 lakon(*)
33  character*20 labmpc(*)
34  character*81 surface,set(*)
35  character*132 textpart(16)
36 !
37  integer ipompc(*),nodempc(3,*),nmpc,nmpc_,mpcfree,istep,istat,
38  & n,i,j,key,nk,node,ifacequad(3,4),ifacetria(3,3),npt,
39  & mpcfreeold,ikmpc(*),ilmpc(*),id,idof,iline,ipol,inl,
40  & ipoinp(2,*),inp(3,*),ipoinpc(0:*),irefnode,lathyp(3,6),inum,
41  & jn,jt,jd,iside,nelem,jface,nopes,nface,nodef(8),nodel(8),
42  & ifaceq(8,6),ifacet(6,4),ifacew1(4,5),ifacew2(8,5),indexpret,
43  & k,ipos,nkold,nope,m,kon(*),ipkon(*),indexe,iset,nset,idir,
44  & istartset(*),iendset(*),ialset(*),index1,ics(2,*),mpcpret,
45  & mint,iflag,ithermal(2),ielem,three,in(3),node1,node2,isign,
46  & ndep,nind,kflag,ne,nkref,noderef
47 !
48  real*8 coefmpc(*),xn(3),xt(3),xd(3),dd,co(3,*),dcs(*),area,
49  & areanodal(8),xl2(3,8),xi,et,weight,shp2(7,8),t0(*),
50  & xs2(3,2),xsj2(3),xsj,yn(3),r
51 !
52  include "gauss.f"
53 !
54 ! latin hypercube positions in a 3 x 3 matrix
55 !
56  data lathyp /1,2,3,1,3,2,2,1,3,2,3,1,3,1,2,3,2,1/
57 !
58 ! nodes per face for hex elements
59 !
60  data ifaceq /4,3,2,1,11,10,9,12,
61  & 5,6,7,8,13,14,15,16,
62  & 1,2,6,5,9,18,13,17,
63  & 2,3,7,6,10,19,14,18,
64  & 3,4,8,7,11,20,15,19,
65  & 4,1,5,8,12,17,16,20/
66 !
67 ! nodes per face for tet elements
68 !
69  data ifacet /1,3,2,7,6,5,
70  & 1,2,4,5,9,8,
71  & 2,3,4,6,10,9,
72  & 1,4,3,8,10,7/
73 !
74 ! nodes per face for linear wedge elements
75 !
76  data ifacew1 /1,3,2,0,
77  & 4,5,6,0,
78  & 1,2,5,4,
79  & 2,3,6,5,
80  & 3,1,4,6/
81 !
82 ! nodes per face for quadratic wedge elements
83 !
84  data ifacew2 /1,3,2,9,8,7,0,0,
85  & 4,5,6,10,11,12,0,0,
86  & 1,2,5,4,7,14,10,13,
87  & 2,3,6,5,8,15,11,14,
88  & 3,1,4,6,9,13,12,15/
89 !
90 ! nodes per face for quad elements
91 !
92  data ifacequad /1,5,2,
93  & 2,6,3,
94  & 3,7,4,
95  & 4,8,1/
96 !
97 ! nodes per face for tria elements
98 !
99  data ifacetria /1,4,2,
100  & 2,5,3,
101  & 3,6,1/
102 !
103 ! flag for shape functions
104 !
105  data iflag /2/
106 !
107  if(istep.gt.0) then
108  write(*,*)
109  & '*ERROR reading *PRE-TENSION SECTION: *EQUATION should'
110  write(*,*) ' be placed before all step definitions'
111  call exit(201)
112  endif
113 !
114  ielem=0
115  do i=1,81
116  surface(i:i)=' '
117  enddo
118 !
119  do i=2,n
120  if(textpart(i)(1:8).eq.'SURFACE=') then
121  if(ielem.ne.0) then
122  write(*,*) '*ERROR reading PRE-TENSION SECTION:'
123  write(*,*) ' ELEMENT and SURFACE are'
124  write(*,*) ' mutually exclusive'
125  call inputerror(inpc,ipoinpc,iline,
126  &"*PRE-TENSION SECTION%")
127  endif
128  surface=textpart(i)(9:88)
129  ipos=index(surface,' ')
130  surface(ipos:ipos)='T'
131  elseif(textpart(i)(1:5).eq.'NODE=') then
132  read(textpart(i)(6:15),'(i10)',iostat=istat) irefnode
133  if(istat.gt.0) call inputerror(inpc,ipoinpc,iline,
134  &"*PRE-TENSION SECTION%")
135  if((irefnode.gt.nk).or.(irefnode.le.0)) then
136  write(*,*) '*ERROR reading *PRE-TENSION SECTION:'
137  write(*,*) ' node ',irefnode,' is not defined'
138  call exit(201)
139  endif
140  elseif(textpart(i)(1:8).eq.'ELEMENT=') then
141  if(surface(1:1).ne.' ') then
142  write(*,*) '*ERROR reading PRE-TENSION SECTION:'
143  write(*,*) ' ELEMENT and SURFACE are'
144  write(*,*) ' mutually exclusive'
145  call inputerror(inpc,ipoinpc,iline,
146  &"*PRE-TENSION SECTION%")
147  endif
148  read(textpart(i)(9:18),'(i10)',iostat=istat) ielem
149  if(istat.gt.0) then
150  write(*,*) '*ERROR reading PRE-TENSION SECTION:'
151  write(*,*) ' cannot read element number'
152  call inputerror(inpc,ipoinpc,iline,
153  &"*PRE-TENSION SECTION%")
154  endif
155  if((ielem.gt.ne).or.(ielem.le.0)) then
156  write(*,*) '*ERROR reading PRE-TENSION SECTION:'
157  write(*,*) ' element',ielem,'is not defined'
158  call inputerror(inpc,ipoinpc,iline,
159  &"*PRE-TENSION SECTION%")
160  endif
161  else
162  write(*,*)
163  & '*WARNING reading *PRE-TENSION SECTION: parameter not recog
164  &nized:'
165  write(*,*) ' ',
166  & textpart(i)(1:index(textpart(i),' ')-1)
167  call inputwarning(inpc,ipoinpc,iline,
168  &"*PRE-TENSION SECTION%")
169  endif
170  enddo
171 !
172 ! checking whether the surface exists and is an element face
173 ! surface
174 !
175  if(surface(1:1).ne.' ') then
176  iset=0
177  do i=1,nset
178  if(set(i).eq.surface) then
179  iset=i
180  exit
181  endif
182  enddo
183  if(iset.eq.0) then
184  write(*,*)
185  & '*ERROR reading *PRE-TENSION SECTION: nonexistent surface'
186  write(*,*) ' or surface consists of nodes'
187  call inputerror(inpc,ipoinpc,iline,
188  &"*PRE-TENSION SECTION%")
189  endif
190  elseif(ielem.gt.0) then
191  if(lakon(ielem)(1:3).ne.'B31') then
192  write(*,*) '*ERROR reading PRE-TENSION SECTION:'
193  write(*,*) ' element',ielem,' is not a linear'
194  write(*,*) ' beam element'
195  call exit(201)
196  endif
197  else
198  write(*,*) '*ERROR reading PRE-TENSION SECTION:'
199  write(*,*) ' either the parameter SURFACE or the'
200  write(*,*) ' parameter ELEMENT must be used'
201  call inputerror(inpc,ipoinpc,iline,
202  &"*PRE-TENSION SECTION%")
203  endif
204 !
205 ! reading the normal vector and normalizing
206 !
207  call getnewline(inpc,textpart,istat,n,key,iline,ipol,inl,
208  & ipoinp,inp,ipoinpc)
209 !
210  do i=1,3
211  read(textpart(i)(1:20),'(f20.0)',iostat=istat) xn(i)
212  if(istat.gt.0) call inputerror(inpc,ipoinpc,iline,
213  &"*PRE-TENSION SECTION%")
214  enddo
215  dd=dsqrt(xn(1)*xn(1)+xn(2)*xn(2)+xn(3)*xn(3))
216  do i=1,3
217  xn(i)=xn(i)/dd
218  enddo
219 !
220 ! beam pretension
221 !
222  if(ielem.gt.0) then
223  indexe=ipkon(ielem)
224 cc
225 !
226 ! removing the beam element (is not needed any more; if not removed
227 ! the beam should have a small cross section compared to the
228 ! pre-tensioned structure in order to get the same results as in
229 ! ABAQUS
230 !
231  ipkon(ielem)=-1
232 cc
233  node1=kon(indexe+1)
234  node2=kon(indexe+2)
235  do i=1,3
236  yn(i)=dabs(xn(i))
237  in(i)=i
238  enddo
239 !
240 ! sorting yn in descending order
241 !
242  three=3
243  kflag=-2
244  call dsort(yn,in,three,kflag)
245 !
246 ! equation: u_2*n-u_1*n+d=0
247 !
248 ! looking for the largest coefficient
249 !
250  do i=1,3
251  if(yn(i).lt.1.d-10) cycle
252 !
253 ! default dependent and independent nodes
254 !
255  ndep=node2
256  nind=node1
257  isign=1
258 !
259  idof=8*(node2-1)+in(i)
260  call nident(ikmpc,idof,nmpc,id)
261  if(id.ne.0) then
262  if(ikmpc(id).eq.idof) then
263  idof=8*(node1-1)+in(i)
264  call nident(ikmpc,idof,nmpc,id)
265  if(id.ne.0) then
266  if(ikmpc(id).eq.idof) then
267  cycle
268  endif
269  endif
270  ndep=node1
271  nind=node2
272  isign=-1
273  endif
274  endif
275 !
276  nmpc=nmpc+1
277  if(nmpc.gt.nmpc_) then
278  write(*,*) '*ERROR in equations: increase nmpc_'
279  call exit(201)
280  endif
281  labmpc(nmpc)='PRETENSION '
282  ipompc(nmpc)=mpcfree
283 !
284 ! updating ikmpc and ilmpc
285 !
286  do j=nmpc,id+2,-1
287  ikmpc(j)=ikmpc(j-1)
288  ilmpc(j)=ilmpc(j-1)
289  enddo
290  ikmpc(id+1)=idof
291  ilmpc(id+1)=nmpc
292 !
293 c jn=in(i)
294  idir=in(i)
295 !
296  if(dabs(xn(idir)).gt.1.d-10) then
297  nodempc(1,mpcfree)=ndep
298  nodempc(2,mpcfree)=idir
299  coefmpc(mpcfree)=isign*xn(idir)
300  mpcfree=nodempc(3,mpcfree)
301 !
302  nodempc(1,mpcfree)=nind
303  nodempc(2,mpcfree)=idir
304  coefmpc(mpcfree)=-isign*xn(idir)
305  mpcfree=nodempc(3,mpcfree)
306  endif
307 !
308  idir=idir+1
309  if(idir.eq.4) idir=1
310  if(dabs(xn(idir)).gt.1.d-10) then
311  nodempc(1,mpcfree)=ndep
312  nodempc(2,mpcfree)=idir
313  coefmpc(mpcfree)=isign*xn(idir)
314  mpcfree=nodempc(3,mpcfree)
315 !
316  nodempc(1,mpcfree)=nind
317  nodempc(2,mpcfree)=idir
318  coefmpc(mpcfree)=-isign*xn(idir)
319  mpcfree=nodempc(3,mpcfree)
320  endif
321 !
322  idir=idir+1
323  if(idir.eq.4) idir=1
324  if(dabs(xn(idir)).gt.1.d-10) then
325  nodempc(1,mpcfree)=ndep
326  nodempc(2,mpcfree)=idir
327  coefmpc(mpcfree)=isign*xn(idir)
328  mpcfree=nodempc(3,mpcfree)
329 !
330  nodempc(1,mpcfree)=nind
331  nodempc(2,mpcfree)=idir
332  coefmpc(mpcfree)=-isign*xn(idir)
333  mpcfree=nodempc(3,mpcfree)
334  endif
335 !
336 ! inhomogeneous term
337 !
338  nodempc(1,mpcfree)=irefnode
339  nodempc(2,mpcfree)=1
340  coefmpc(mpcfree)=1.d0
341  mpcfreeold=mpcfree
342  mpcfree=nodempc(3,mpcfree)
343  nodempc(3,mpcfreeold)=0
344 !
345  call getnewline(inpc,textpart,istat,n,key,iline,ipol,inl,
346  & ipoinp,inp,ipoinpc)
347 !
348  return
349  enddo
350  write(*,*) '*ERROR reading *PRE-TENSION SECTION'
351  write(*,*) ' all DOFS of the beam elements'
352  write(*,*) ' have been used previously'
353  call exit(201)
354  endif
355 !
356 ! finding a unit vector xt perpendicular to the normal vector
357 ! using a unit vector in x or in y
358 !
359  if(dabs(xn(1)).lt.0.95d0) then
360  xt(1)=1.d0-xn(1)*xn(1)
361  xt(2)=-xn(1)*xn(2)
362  xt(3)=-xn(1)*xn(3)
363  else
364  xt(1)=-xn(2)*xn(1)
365  xt(2)=1.d0-xn(2)*xn(2)
366  xt(3)=-xn(2)*xn(3)
367  endif
368  dd=dsqrt(xt(1)*xt(1)+xt(2)*xt(2)+xt(3)*xt(3))
369  do i=1,3
370  xt(i)=xt(i)/dd
371  enddo
372 !
373 ! xd=xn x xt
374 !
375  xd(1)=xn(2)*xt(3)-xn(3)*xt(2)
376  xd(2)=xn(3)*xt(1)-xn(1)*xt(3)
377  xd(3)=xn(1)*xt(2)-xn(2)*xt(1)
378 !
379 ! generating a Latin hypercube
380 ! checking which DOF's of xn, xt and xd are nonzero
381 !
382  do inum=1,6
383  if((dabs(xn(lathyp(1,inum))).gt.1.d-3).and.
384  & (dabs(xt(lathyp(2,inum))).gt.1.d-3).and.
385  & (dabs(xd(lathyp(3,inum))).gt.1.d-3)) exit
386  enddo
387  jn=lathyp(1,inum)
388  jt=lathyp(2,inum)
389  jd=lathyp(3,inum)
390 !
391 ! generating the MPCs
392 !
393  indexpret=0
394  nkold=nk
395  m=iendset(iset)-istartset(iset)+1
396 !
397 ! number of distinct pre-strain nodes for the present keyword
398 !
399  npt=0
400  area=0.d0
401 !
402 ! loop over all element faces belonging to the surface
403 !
404  do k=1,m
405  twod=.false.
406  iside=ialset(istartset(iset)+k-1)
407  nelem=int(iside/10.d0)
408  indexe=ipkon(nelem)
409  jface=iside-10*nelem
410 !
411 ! nodes: #nodes in the face
412 ! the nodes are stored in nodef(*)
413 !
414  if(lakon(nelem)(4:4).eq.'2') then
415  nopes=8
416  nface=6
417  elseif(lakon(nelem)(3:4).eq.'D8') then
418  nopes=4
419  nface=6
420  elseif(lakon(nelem)(4:5).eq.'10') then
421  nopes=6
422  nface=4
423  nope=10
424  elseif(lakon(nelem)(4:4).eq.'4') then
425  nopes=3
426  nface=4
427  nope=4
428  elseif(lakon(nelem)(4:5).eq.'15') then
429  if(jface.le.2) then
430  nopes=6
431  else
432  nopes=8
433  endif
434  nface=5
435  nope=15
436  elseif(lakon(nelem)(3:4).eq.'D6') then
437  if(jface.le.2) then
438  nopes=3
439  else
440  nopes=4
441  endif
442  nface=5
443  nope=6
444  elseif((lakon(nelem)(2:2).eq.'8').or.
445  & (lakon(nelem)(4:4).eq.'8')) then
446 !
447 ! 8-node 2-D elements
448 !
449  nopes=3
450  nface=4
451  nope=8
452  if(lakon(nelem)(4:4).eq.'8') then
453  twod=.true.
454  jface=jface-2
455  endif
456  elseif((lakon(nelem)(2:2).eq.'6').or.
457  & (lakon(nelem)(4:4).eq.'6')) then
458 !
459 ! 6-node 2-D elements
460 !
461  nopes=3
462  nface=3
463  if(lakon(nelem)(4:4).eq.'6') then
464  twod=.true.
465  jface=jface-2
466  endif
467  else
468  cycle
469  endif
470 !
471 ! determining the nodes of the face
472 !
473  if(nface.eq.3) then
474  do i=1,nopes
475  nodef(i)=kon(indexe+ifacetria(i,jface))
476  nodel(i)=ifacetria(i,jface)
477  enddo
478  elseif(nface.eq.4) then
479  if(nope.eq.8) then
480  do i=1,nopes
481  nodef(i)=kon(indexe+ifacequad(i,jface))
482  nodel(i)=ifacequad(i,jface)
483  enddo
484  else
485  do i=1,nopes
486  nodef(i)=kon(indexe+ifacet(i,jface))
487  nodel(i)=ifacet(i,jface)
488  enddo
489  endif
490  elseif(nface.eq.5) then
491  if(nope.eq.6) then
492  do i=1,nopes
493  nodef(i)=kon(indexe+ifacew1(i,jface))
494  nodel(i)=ifacew1(i,jface)
495  enddo
496  elseif(nope.eq.15) then
497  do i=1,nopes
498  nodef(i)=kon(indexe+ifacew2(i,jface))
499  nodel(i)=ifacew2(i,jface)
500  enddo
501  endif
502  elseif(nface.eq.6) then
503  do i=1,nopes
504  nodef(i)=kon(indexe+ifaceq(i,jface))
505  nodel(i)=ifaceq(i,jface)
506  enddo
507  endif
508 !
509 ! loop over the nodes belonging to the face
510 ! ics(1,*): pretension node
511 ! ics(2,*): corresponding partner node
512 ! dcs(*): area corresponding to pretension node
513 !
514  do i=1,nopes
515  node=nodef(i)
516  call nident2(ics,node,npt,id)
517  if(id.gt.0) then
518  if(ics(1,id).eq.node) then
519 !
520 ! node was already treated: replacing the node
521 ! by the partner node
522 !
523  kon(indexe+nodel(i))=ics(2,id)
524  cycle
525  endif
526  endif
527 !
528 ! generating a partner node
529 !
530  nk=nk+1
531 !
532 ! coordinates for the new node
533 !
534  do j=1,3
535  co(j,nk)=co(j,node)
536  enddo
537 !
538 ! updating the topology
539 !
540  kon(indexe+nodel(i))=nk
541 !
542 ! updating ics
543 !
544  npt=npt+1
545  do j=npt,id+2,-1
546  ics(1,j)=ics(1,j-1)
547  ics(2,j)=ics(2,j-1)
548  dcs(j)=dcs(j-1)
549  enddo
550  ics(1,id+1)=node
551  ics(2,id+1)=nk
552  dcs(id+1)=0.d0
553 !
554 ! first MPC perpendicular to the normal direction
555 !
556  idof=8*(nk-1)+jt
557  call nident(ikmpc,idof,nmpc,id)
558 !
559  nmpc=nmpc+1
560  if(nmpc.gt.nmpc_) then
561  write(*,*) '*ERROR in equations: increase nmpc_'
562  call exit(201)
563  endif
564  ipompc(nmpc)=mpcfree
565  labmpc(nmpc)=' '
566 !
567 ! updating ikmpc and ilmpc
568 !
569  do j=nmpc,id+2,-1
570  ikmpc(j)=ikmpc(j-1)
571  ilmpc(j)=ilmpc(j-1)
572  enddo
573  ikmpc(id+1)=idof
574  ilmpc(id+1)=nmpc
575 !
576  idir=jt
577  if(dabs(xt(idir)).gt.1.d-10) then
578  nodempc(1,mpcfree)=nk
579  nodempc(2,mpcfree)=idir
580  coefmpc(mpcfree)=-xt(idir)
581  mpcfreeold=mpcfree
582  mpcfree=nodempc(3,mpcfree)
583  endif
584 !
585  idir=idir+1
586  if(idir.eq.4) idir=1
587  if(dabs(xt(idir)).gt.1.d-10) then
588  nodempc(1,mpcfree)=nk
589  nodempc(2,mpcfree)=idir
590  coefmpc(mpcfree)=-xt(idir)
591  mpcfreeold=mpcfree
592  mpcfree=nodempc(3,mpcfree)
593  endif
594 !
595  idir=idir+1
596  if(idir.eq.4) idir=1
597  if(dabs(xt(idir)).gt.1.d-10) then
598  nodempc(1,mpcfree)=nk
599  nodempc(2,mpcfree)=idir
600  coefmpc(mpcfree)=-xt(idir)
601  mpcfreeold=mpcfree
602  mpcfree=nodempc(3,mpcfree)
603  endif
604 !
605  idir=jt
606  if(dabs(xt(idir)).gt.1.d-10) then
607  nodempc(1,mpcfree)=node
608  nodempc(2,mpcfree)=jt
609  coefmpc(mpcfree)=xt(idir)
610  mpcfreeold=mpcfree
611  mpcfree=nodempc(3,mpcfree)
612  endif
613 !
614  idir=idir+1
615  if(idir.eq.4) idir=1
616  if(dabs(xt(idir)).gt.1.d-10) then
617  nodempc(1,mpcfree)=node
618  nodempc(2,mpcfree)=idir
619  coefmpc(mpcfree)=xt(idir)
620  mpcfreeold=mpcfree
621  mpcfree=nodempc(3,mpcfree)
622  endif
623 !
624  idir=idir+1
625  if(idir.eq.4) idir=1
626  if(dabs(xt(idir)).gt.1.d-10) then
627  nodempc(1,mpcfree)=node
628  nodempc(2,mpcfree)=idir
629  coefmpc(mpcfree)=xt(idir)
630  mpcfreeold=mpcfree
631  mpcfree=nodempc(3,mpcfree)
632  endif
633  nodempc(3,mpcfreeold)=0
634 !
635 ! second MPC perpendicular to the normal direction
636 !
637  if(.not.twod) then
638  idof=8*(nk-1)+jd
639  call nident(ikmpc,idof,nmpc,id)
640 !
641  nmpc=nmpc+1
642  if(nmpc.gt.nmpc_) then
643  write(*,*) '*ERROR in equations: increase nmpc_'
644  call exit(201)
645  endif
646  labmpc(nmpc)=' '
647  ipompc(nmpc)=mpcfree
648 !
649 ! updating ikmpc and ilmpc
650 !
651  do j=nmpc,id+2,-1
652  ikmpc(j)=ikmpc(j-1)
653  ilmpc(j)=ilmpc(j-1)
654  enddo
655  ikmpc(id+1)=idof
656  ilmpc(id+1)=nmpc
657 !
658  idir=jd
659  if(dabs(xd(idir)).gt.1.d-10) then
660  nodempc(1,mpcfree)=nk
661  nodempc(2,mpcfree)=idir
662  coefmpc(mpcfree)=-xd(idir)
663  mpcfreeold=mpcfree
664  mpcfree=nodempc(3,mpcfree)
665  endif
666 !
667  idir=idir+1
668  if(idir.eq.4) idir=1
669  if(dabs(xd(idir)).gt.1.d-10) then
670  nodempc(1,mpcfree)=nk
671  nodempc(2,mpcfree)=idir
672  coefmpc(mpcfree)=-xd(idir)
673  mpcfreeold=mpcfree
674  mpcfree=nodempc(3,mpcfree)
675  endif
676 !
677  idir=idir+1
678  if(idir.eq.4) idir=1
679  if(dabs(xd(idir)).gt.1.d-10) then
680  nodempc(1,mpcfree)=nk
681  nodempc(2,mpcfree)=idir
682  coefmpc(mpcfree)=-xd(idir)
683  mpcfreeold=mpcfree
684  mpcfree=nodempc(3,mpcfree)
685  endif
686 !
687  idir=jd
688  if(dabs(xd(idir)).gt.1.d-10) then
689  nodempc(1,mpcfree)=node
690  nodempc(2,mpcfree)=idir
691  coefmpc(mpcfree)=xd(idir)
692  mpcfreeold=mpcfree
693  mpcfree=nodempc(3,mpcfree)
694  endif
695 !
696  idir=idir+1
697  if(idir.eq.4) idir=1
698  if(dabs(xd(idir)).gt.1.d-10) then
699  nodempc(1,mpcfree)=node
700  nodempc(2,mpcfree)=idir
701  coefmpc(mpcfree)=xd(idir)
702  mpcfreeold=mpcfree
703  mpcfree=nodempc(3,mpcfree)
704  endif
705 !
706  idir=idir+1
707  if(idir.eq.4) idir=1
708  if(dabs(xd(idir)).gt.1.d-10) then
709  nodempc(1,mpcfree)=node
710  nodempc(2,mpcfree)=idir
711  coefmpc(mpcfree)=xd(idir)
712  mpcfreeold=mpcfree
713  mpcfree=nodempc(3,mpcfree)
714  endif
715  nodempc(3,mpcfreeold)=0
716  endif
717 !
718 ! MPC in normal direction
719 !
720 ! check whether initialized
721 !
722  if(indexpret.eq.0) then
723  idof=8*(nk-1)+jn
724  call nident(ikmpc,idof,nmpc,id)
725 !
726  nmpc=nmpc+1
727  if(nmpc.gt.nmpc_) then
728  write(*,*) '*ERROR in equations: increase nmpc_'
729  call exit(201)
730  endif
731  labmpc(nmpc)='PRETENSION '
732  ipompc(nmpc)=mpcfree
733  mpcpret=nmpc
734 !
735 ! taking the first node as reference node for the condition:
736 ! "the displacement in pre-tension direction in all nodes
737 ! should be the same as for the reference node" (pre-tension
738 ! surfaces stay parallel to each other)
739 !
740  nkref=nk
741  noderef=node
742 !
743 ! updating ikmpc and ilmpc
744 !
745  do j=nmpc,id+2,-1
746  ikmpc(j)=ikmpc(j-1)
747  ilmpc(j)=ilmpc(j-1)
748  enddo
749  ikmpc(id+1)=idof
750  ilmpc(id+1)=nmpc
751 c else
752 c nodempc(3,indexpret)=mpcfree
753  endif
754 !
755 ! MPC's specifying that the pre-tension surfaces should stay
756 ! parallel
757 !
758  if(indexpret.ne.0) then
759  idof=8*(nk-1)+jn
760  call nident(ikmpc,idof,nmpc,id)
761 !
762  nmpc=nmpc+1
763  if(nmpc.gt.nmpc_) then
764  write(*,*) '*ERROR in equations: increase nmpc_'
765  call exit(201)
766  endif
767  ipompc(nmpc)=mpcfree
768  labmpc(nmpc)=' '
769 !
770 ! updating ikmpc and ilmpc
771 !
772  do j=nmpc,id+2,-1
773  ikmpc(j)=ikmpc(j-1)
774  ilmpc(j)=ilmpc(j-1)
775  enddo
776  ikmpc(id+1)=idof
777  ilmpc(id+1)=nmpc
778 !
779  idir=jn
780  if(dabs(xn(idir)).gt.1.d-10) then
781  nodempc(1,mpcfree)=nk
782  nodempc(2,mpcfree)=idir
783  coefmpc(mpcfree)=-xn(idir)
784  mpcfreeold=mpcfree
785  mpcfree=nodempc(3,mpcfree)
786  endif
787 !
788  idir=idir+1
789  if(idir.eq.4) idir=1
790  if(dabs(xn(idir)).gt.1.d-10) then
791  nodempc(1,mpcfree)=nk
792  nodempc(2,mpcfree)=idir
793  coefmpc(mpcfree)=-xn(idir)
794  mpcfreeold=mpcfree
795  mpcfree=nodempc(3,mpcfree)
796  endif
797 !
798  idir=idir+1
799  if(idir.eq.4) idir=1
800  if(dabs(xn(idir)).gt.1.d-10) then
801  nodempc(1,mpcfree)=nk
802  nodempc(2,mpcfree)=idir
803  coefmpc(mpcfree)=-xn(idir)
804  mpcfreeold=mpcfree
805  mpcfree=nodempc(3,mpcfree)
806  endif
807 !
808  idir=jn
809  if(dabs(xn(idir)).gt.1.d-10) then
810  nodempc(1,mpcfree)=node
811  nodempc(2,mpcfree)=idir
812  coefmpc(mpcfree)=xn(idir)
813  mpcfreeold=mpcfree
814  mpcfree=nodempc(3,mpcfree)
815  endif
816 !
817  idir=idir+1
818  if(idir.eq.4) idir=1
819  if(dabs(xn(idir)).gt.1.d-10) then
820  nodempc(1,mpcfree)=node
821  nodempc(2,mpcfree)=idir
822  coefmpc(mpcfree)=xn(idir)
823  mpcfreeold=mpcfree
824  mpcfree=nodempc(3,mpcfree)
825  endif
826 !
827  idir=idir+1
828  if(idir.eq.4) idir=1
829  if(dabs(xn(idir)).gt.1.d-10) then
830  nodempc(1,mpcfree)=node
831  nodempc(2,mpcfree)=idir
832  coefmpc(mpcfree)=xn(idir)
833  mpcfreeold=mpcfree
834  mpcfree=nodempc(3,mpcfree)
835  endif
836 !
837 ! corresponding terms in the reference node
838 !
839  idir=jn
840  if(dabs(xn(idir)).gt.1.d-10) then
841  nodempc(1,mpcfree)=nkref
842  nodempc(2,mpcfree)=idir
843  coefmpc(mpcfree)=xn(idir)
844  mpcfreeold=mpcfree
845  mpcfree=nodempc(3,mpcfree)
846  endif
847 !
848  idir=idir+1
849  if(idir.eq.4) idir=1
850  if(dabs(xn(idir)).gt.1.d-10) then
851  nodempc(1,mpcfree)=nkref
852  nodempc(2,mpcfree)=idir
853  coefmpc(mpcfree)=xn(idir)
854  mpcfreeold=mpcfree
855  mpcfree=nodempc(3,mpcfree)
856  endif
857 !
858  idir=idir+1
859  if(idir.eq.4) idir=1
860  if(dabs(xn(idir)).gt.1.d-10) then
861  nodempc(1,mpcfree)=nkref
862  nodempc(2,mpcfree)=idir
863  coefmpc(mpcfree)=xn(idir)
864  mpcfreeold=mpcfree
865  mpcfree=nodempc(3,mpcfree)
866  endif
867 !
868  idir=jn
869  if(dabs(xn(idir)).gt.1.d-10) then
870  nodempc(1,mpcfree)=noderef
871  nodempc(2,mpcfree)=idir
872  coefmpc(mpcfree)=-xn(idir)
873  mpcfreeold=mpcfree
874  mpcfree=nodempc(3,mpcfree)
875  endif
876 !
877  idir=idir+1
878  if(idir.eq.4) idir=1
879  if(dabs(xn(idir)).gt.1.d-10) then
880  nodempc(1,mpcfree)=noderef
881  nodempc(2,mpcfree)=idir
882  coefmpc(mpcfree)=-xn(idir)
883  mpcfreeold=mpcfree
884  mpcfree=nodempc(3,mpcfree)
885  endif
886 !
887  idir=idir+1
888  if(idir.eq.4) idir=1
889  if(dabs(xn(idir)).gt.1.d-10) then
890  nodempc(1,mpcfree)=noderef
891  nodempc(2,mpcfree)=idir
892  coefmpc(mpcfree)=-xn(idir)
893  mpcfreeold=mpcfree
894  mpcfree=nodempc(3,mpcfree)
895  endif
896  nodempc(3,mpcfreeold)=0
897 !
898  nodempc(3,indexpret)=mpcfree
899 c endif
900  else
901 !
902 ! governing pre-tension equation
903 !
904  idir=jn
905  if(dabs(xn(idir)).gt.1.d-10) then
906  nodempc(1,mpcfree)=nk
907  nodempc(2,mpcfree)=idir
908  coefmpc(mpcfree)=-xn(idir)
909  indexpret=mpcfree
910  mpcfree=nodempc(3,mpcfree)
911 !
912  nodempc(1,mpcfree)=node
913  nodempc(2,mpcfree)=idir
914  coefmpc(mpcfree)=xn(idir)
915  indexpret=mpcfree
916  mpcfree=nodempc(3,mpcfree)
917  endif
918 !
919  idir=idir+1
920  if(idir.eq.4) idir=1
921  if(dabs(xn(idir)).gt.1.d-10) then
922  nodempc(1,mpcfree)=nk
923  nodempc(2,mpcfree)=idir
924  coefmpc(mpcfree)=-xn(idir)
925  indexpret=mpcfree
926  mpcfree=nodempc(3,mpcfree)
927 !
928  nodempc(1,mpcfree)=node
929  nodempc(2,mpcfree)=idir
930  coefmpc(mpcfree)=xn(idir)
931  indexpret=mpcfree
932  mpcfree=nodempc(3,mpcfree)
933  endif
934 !
935  idir=idir+1
936  if(idir.eq.4) idir=1
937  if(dabs(xn(idir)).gt.1.d-10) then
938  nodempc(1,mpcfree)=nk
939  nodempc(2,mpcfree)=idir
940  coefmpc(mpcfree)=-xn(idir)
941  indexpret=mpcfree
942  mpcfree=nodempc(3,mpcfree)
943 !
944  nodempc(1,mpcfree)=node
945  nodempc(2,mpcfree)=idir
946  coefmpc(mpcfree)=xn(idir)
947  indexpret=mpcfree
948  mpcfree=nodempc(3,mpcfree)
949  endif
950  endif
951 !
952 ! thermal MPC
953 !
954  if(ithermal(2).gt.0) then
955  idof=8*(nk-1)
956  call nident(ikmpc,idof,nmpc,id)
957 !
958  nmpc=nmpc+1
959  if(nmpc.gt.nmpc_) then
960  write(*,*) '*ERROR in equations: increase nmpc_'
961  call exit(201)
962  endif
963 !
964 ! the label for purely mechanical calculations is
965 ! necessary in order to copy t1 at the newly created
966 ! nodes in file calinput.f
967 !
968  if(ithermal(2).gt.1) then
969  labmpc(nmpc)=' '
970  else
971  labmpc(nmpc)='THERMALPRET '
972  endif
973 !
974  ipompc(nmpc)=mpcfree
975 !
976 ! updating ikmpc and ilmpc
977 !
978  do j=nmpc,id+2,-1
979  ikmpc(j)=ikmpc(j-1)
980  ilmpc(j)=ilmpc(j-1)
981  enddo
982  ikmpc(id+1)=idof
983  ilmpc(id+1)=nmpc
984 !
985  nodempc(1,mpcfree)=nk
986  nodempc(2,mpcfree)=0
987  coefmpc(mpcfree)=1.d0
988  mpcfree=nodempc(3,mpcfree)
989 !
990  nodempc(1,mpcfree)=node
991  nodempc(2,mpcfree)=0
992  coefmpc(mpcfree)=-1.d0
993 !
994  mpcfreeold=mpcfree
995  mpcfree=nodempc(3,mpcfree)
996  nodempc(3,mpcfreeold)=0
997  endif
998 !
999 ! initial value for new node
1000 !
1001  if(ithermal(1).gt.0) then
1002  t0(nk)=t0(node)
1003  endif
1004 !
1005  enddo
1006 !
1007 ! calculating the area of the face and its contributions
1008 ! to the facial nodes
1009 !
1010 ! number of integration points
1011 !
1012  if(lakon(nelem)(3:5).eq.'D8R') then
1013  mint=1
1014  elseif(lakon(nelem)(3:4).eq.'D8') then
1015  mint=4
1016  elseif(lakon(nelem)(4:6).eq.'20R') then
1017  mint=4
1018  elseif(lakon(nelem)(4:4).eq.'2') then
1019  mint=9
1020  elseif(lakon(nelem)(4:5).eq.'10') then
1021  mint=3
1022  elseif(lakon(nelem)(4:4).eq.'4') then
1023  mint=1
1024  elseif(lakon(nelem)(3:4).eq.'D6') then
1025  mint=1
1026  elseif(lakon(nelem)(4:5).eq.'15') then
1027  if(jface.le.2) then
1028  mint=3
1029  else
1030  mint=4
1031  endif
1032 !
1033 ! faces of 2-D elements
1034 !
1035  elseif((lakon(nelem)(3:3).eq.'R').or.
1036  & (lakon(nelem)(5:5).eq.'R')) then
1037  mint=2
1038  else
1039  mint=3
1040  endif
1041 !
1042  do i=1,nopes
1043  areanodal(i)=0.d0
1044  do j=1,3
1045  xl2(j,i)=co(j,nodef(i))
1046  enddo
1047  enddo
1048 !
1049  do m=1,mint
1050  if((lakon(nelem)(3:5).eq.'D8R').or.
1051  & ((lakon(nelem)(3:4).eq.'D6').and.(nopes.eq.4))) then
1052  xi=gauss2d1(1,m)
1053  et=gauss2d1(2,m)
1054  weight=weight2d1(m)
1055  elseif((lakon(nelem)(3:4).eq.'D8').or.
1056  & (lakon(nelem)(4:6).eq.'20R').or.
1057  & ((lakon(nelem)(4:5).eq.'15').and.
1058  & (nopes.eq.8))) then
1059  xi=gauss2d2(1,m)
1060  et=gauss2d2(2,m)
1061  weight=weight2d2(m)
1062  elseif(lakon(nelem)(4:4).eq.'2') then
1063  xi=gauss2d3(1,m)
1064  et=gauss2d3(2,m)
1065  weight=weight2d3(m)
1066  elseif((lakon(nelem)(4:5).eq.'10').or.
1067  & ((lakon(nelem)(4:5).eq.'15').and.
1068  & (nopes.eq.6))) then
1069  xi=gauss2d5(1,m)
1070  et=gauss2d5(2,m)
1071  weight=weight2d5(m)
1072  elseif((lakon(nelem)(4:4).eq.'4').or.
1073  & ((lakon(nelem)(3:4).eq.'D6').and.
1074  & (nopes.eq.3))) then
1075  xi=gauss2d4(1,m)
1076  et=gauss2d4(2,m)
1077  weight=weight2d4(m)
1078 !
1079 ! faces of 2-D elements
1080 !
1081  elseif((lakon(nelem)(3:3).eq.'R').or.
1082  & (lakon(nelem)(5:5).eq.'R')) then
1083  xi=gauss1d2(1,m)
1084  weight=weight1d2(m)
1085  else
1086  xi=gauss1d3(1,m)
1087  weight=weight1d3(m)
1088  endif
1089 !
1090  if(nopes.eq.8) then
1091  call shape8q(xi,et,xl2,xsj2,xs2,shp2,iflag)
1092  elseif(nopes.eq.4) then
1093  call shape4q(xi,et,xl2,xsj2,xs2,shp2,iflag)
1094  elseif(nopes.eq.6) then
1095  call shape6tri(xi,et,xl2,xsj2,xs2,shp2,iflag)
1096  elseif((nopes.eq.3).and.(.not.twod)) then
1097  call shape3tri(xi,et,xl2,xsj2,xs2,shp2,iflag)
1098  else
1099 !
1100 ! 3-node line
1101 !
1102  call shape3l(xi,xl2,xsj2,xs2,shp2,iflag)
1103  endif
1104 !
1105 ! calculating the total area and nodal area
1106 !
1107  if(.not.twod) then
1108  xsj=weight*dsqrt(xsj2(1)**2+xsj2(2)**2+xsj2(3)**2)
1109  elseif(lakon(nelem)(1:3).eq.'CAX') then
1110 !
1111 ! radial distance is taken into account for the area
1112 !
1113  r=0.d0
1114  do i=1,nopes
1115  r=r+shp2(4,i)*xl2(1,i)
1116  enddo
1117  xsj=weight*xsj2(1)*r
1118  else
1119  xsj=weight*xsj2(1)
1120  endif
1121  area=area+xsj
1122  do i=1,nopes
1123  areanodal(i)=areanodal(i)+xsj*shp2(4,i)
1124  enddo
1125 !
1126  enddo
1127 !
1128 ! inserting the nodal area into field dcs
1129 !
1130  do i=1,nopes
1131  node=nodef(i)
1132  call nident2(ics,node,npt,id)
1133  dcs(id)=dcs(id)+areanodal(i)
1134  enddo
1135 !
1136  enddo
1137 !
1138  nodempc(3,indexpret)=mpcfree
1139  nodempc(1,mpcfree)=irefnode
1140  nodempc(2,mpcfree)=1
1141  coefmpc(mpcfree)=area
1142  mpcfreeold=mpcfree
1143  mpcfree=nodempc(3,mpcfree)
1144  nodempc(3,mpcfreeold)=0
1145 !
1146 ! changing the coefficients of the pretension MPC
1147 !
1148  index1=ipompc(mpcpret)
1149  do
1150  node=nodempc(1,nodempc(3,index1))
1151  call nident2(ics,node,npt,id)
1152  do j=1,2
1153 c coefmpc(index1)=coefmpc(index1)*dcs(id)
1154  coefmpc(index1)=coefmpc(index1)*area
1155  index1=nodempc(3,index1)
1156  enddo
1157  if(nodempc(1,index1).eq.irefnode) exit
1158  enddo
1159 !
1160  call getnewline(inpc,textpart,istat,n,key,iline,ipol,inl,
1161  & ipoinp,inp,ipoinpc)
1162 !
1163 c do i=1,nmpc
1164 c call writempc(ipompc,nodempc,coefmpc,labmpc,i)
1165 c enddo
1166 c do i=1,nmpc
1167 c write(*,*) i,ikmpc(i),ilmpc(i)
1168 c enddo
1169 !
1170  return
shape8q
subroutine shape8q(xi, et, xl, xsj, xs, shp, iflag)
Definition: shape8q.f:20
shape3l
subroutine shape3l(xi, xl, xsj, xs, shp, iflag)
Definition: shape3l.f:20
shape3tri
subroutine shape3tri(xi, et, xl, xsj, xs, shp, iflag)
Definition: shape3tri.f:20
inputwarning
subroutine inputwarning(inpc, ipoinpc, iline, text)
Definition: inputwarning.f:20
nident2
subroutine nident2(x, px, n, id)
Definition: nident2.f:27
shape4q
subroutine shape4q(xi, et, xl, xsj, xs, shp, iflag)
Definition: shape4q.f:20
getnewline
subroutine getnewline(inpc, textpart, istat, n, key, iline, ipol, inl, ipoinp, inp, ipoinpc)
Definition: getnewline.f:21
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
inputerror
subroutine inputerror(inpc, ipoinpc, iline, text)
Definition: inputerror.f:20
Hosted by OpenAircraft.com, (Michigan UAV, LLC)