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

Go to the source code of this file.

Functions/Subroutines

subroutine couplings (inpc, textpart, set, istartset, iendset, ialset, nset, nboun, nk, ipompc, nodempc, coefmpc, nmpc, nmpc_, mpcfree, ikboun, ikmpc, ilmpc, co, labmpc, istat, n, iline, ipol, inl, ipoinp, inp, ipoinpc, norien, orname, orab, irstrt, ipkon, kon, lakon, istep, ics, dcs, nk_, nboun_, nodeboun, ndirboun, typeboun, ilboun, xboun)
 

Function/Subroutine Documentation

◆ couplings()

subroutine couplings ( character*1, dimension(*)  inpc,
character*132, dimension(16)  textpart,
character*81, dimension(*)  set,
integer, dimension(*)  istartset,
integer, dimension(*)  iendset,
integer, dimension(*)  ialset,
integer  nset,
integer  nboun,
integer  nk,
integer, dimension(*)  ipompc,
integer, dimension(3,*)  nodempc,
real*8, dimension(*)  coefmpc,
integer  nmpc,
integer  nmpc_,
integer  mpcfree,
integer, dimension(*)  ikboun,
integer, dimension(*)  ikmpc,
integer, dimension(*)  ilmpc,
real*8, dimension(3,*)  co,
character*20, dimension(*)  labmpc,
integer  istat,
integer  n,
integer  iline,
integer  ipol,
integer  inl,
integer, dimension(2,*)  ipoinp,
integer, dimension(3,*)  inp,
integer, dimension(0:*)  ipoinpc,
integer  norien,
character*80, dimension(*)  orname,
real*8, dimension(7,*)  orab,
integer  irstrt,
integer, dimension(*)  ipkon,
integer, dimension(*)  kon,
character*8, dimension(*)  lakon,
integer  istep,
integer, dimension(2,*)  ics,
real*8, dimension(*)  dcs,
integer  nk_,
integer  nboun_,
integer, dimension(*)  nodeboun,
integer, dimension(*)  ndirboun,
character*1, dimension(*)  typeboun,
integer, dimension(*)  ilboun,
real*8, dimension(*)  xboun 
)
25 !
26 ! reading the input deck: *COUPLING in combination with
27 ! *KINEMATIC or *DISTRIBUTING
28 !
29  implicit none
30 !
31  character*1 inpc(*),surfkind,typeboun(*)
32  character*8 lakon(*)
33  character*20 labmpc(*),label
34  character*80 orname(*),orientation
35  character*81 set(*),surfset
36  character*132 textpart(16),name
37 !
38  integer istartset(*),iendset(*),ialset(*),norien,irstrt,nface,
39  & iorientation,iface,jface,nset,nboun,istat,n,i,j,k,ibounstart,
40  & ibounend,key,nk,ipompc(*),nodempc(3,*),nmpc,nmpc_,mpcfree,
41  & ikboun(*),ikmpc(*),ilmpc(*),ipos,m,node,iline,ipol,inl,nope,
42  & ipoinp(2,*),inp(3,*),ipoinpc(0:*),jsurf,irefnode,indexe,nopes,
43  & ifaceq(8,6),ifacet(6,4),ifacew1(4,5),ifacew2(8,5),ipkon(*),
44  & kon(*),istep,nelem,ics(2,*),nodef(8),nboun_,nodeboun(*),
45  & npt,mint,index1,mpcfreeold,id,idof,iflag,inhomnode,nk_,kk,
46  & irotnode(3),l,index2,indexold(3),indexnew(3),idirold,idirmax,
47  & idir,indexmax,nodeold,node1,nodemax,ndirboun(*),ilboun(*),
48  & irotnode_kin,idupnode
49 !
50  real*8 coefmpc(*),co(3,*),orab(7,*),dcs(*),areanodal(8),xl2(3,8),
51  & shp2(7,8),xsj2(3),xsj,xi,et,weight,xs2(3,2),area,a(3,3),cgx(3),
52  & aa(3),pi(3),c1,c4,c9,c10,amax,xcoef,coef(3),xboun(*)
53 !
54  include "gauss.f"
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 !
92  if((istep.gt.0).and.(irstrt.ge.0)) then
93  write(*,*) '*ERROR reading *COUPLING: *COUPLING'
94  write(*,*)' should be placed before all step definitions'
95  call exit(201)
96  endif
97 !
98  label=' '
99  orientation='
100  & '
101  do i=1,81
102  surfset(i:i)=' '
103  enddo
104 !
105  name(1:1)=' '
106  do i=2,n
107  if(textpart(i)(1:8).eq.'REFNODE=') then
108  read(textpart(i)(9:18),'(i10)',iostat=istat) irefnode
109  if(istat.gt.0) call inputerror(inpc,ipoinpc,iline,
110  &"*COUPLING%")
111  if(irefnode.gt.nk) then
112  write(*,*) '*ERROR reading *COUPLING: ref node',irefnode
113  write(*,*) ' has not been defined'
114  call exit(201)
115  endif
116  else if(textpart(i)(1:8).eq.'SURFACE=') then
117  surfset(1:80)=textpart(i)(9:88)
118 !
119  ipos=index(surfset,' ')
120  surfkind='S'
121  surfset(ipos:ipos)=surfkind
122  do j=1,nset
123  if(set(j).eq.surfset) exit
124  enddo
125  if(j.gt.nset) then
126  surfkind='T'
127  surfset(ipos:ipos)=surfkind
128  do j=1,nset
129  if(set(j).eq.surfset) exit
130  enddo
131  if(j.gt.nset) then
132  write(*,*) '*ERROR reading *COUPLING:'
133  write(*,*) ' surface ',surfset
134  write(*,*) ' has not yet been defined.'
135  call exit(201)
136  endif
137  endif
138  jsurf=j
139  elseif(textpart(i)(1:12).eq.'ORIENTATION=') then
140  orientation=textpart(i)(13:92)
141  elseif(textpart(i)(1:15).eq.'CONSTRAINTNAME=') then
142  name(1:117)=textpart(i)(16:132)
143  else
144  write(*,*)
145  & '*WARNING reading *COUPLING: parameter not recognized:'
146  write(*,*) ' ',
147  & textpart(i)(1:index(textpart(i),' ')-1)
148  call inputwarning(inpc,ipoinpc,iline,
149  &"*COUPLING%")
150  endif
151  enddo
152 !
153  if(name(1:1).eq.' ') then
154  write(*,*)
155  & '*ERROR reading *COUPLING: no CONTRAINT NAME given'
156  write(*,*) ' '
157  call inputerror(inpc,ipoinpc,iline,
158  &"*COUPLING%")
159  call exit(201)
160  endif
161 !
162  if(orientation.eq.' ') then
163  iorientation=0
164  else
165  do i=1,norien
166  if(orname(i).eq.orientation) exit
167  enddo
168  if(i.gt.norien) then
169  write(*,*)
170  & '*ERROR reading *COUPLING: nonexistent orientation'
171  write(*,*) ' '
172  call inputerror(inpc,ipoinpc,iline,
173  &"*COUPLING%")
174  call exit(201)
175  endif
176  iorientation=i
177  endif
178 !
179 ! next keyword should be *KINEMATIC or *DISTRIBUTING
180 !
181  call getnewline(inpc,textpart,istat,n,key,iline,ipol,inl,
182  & ipoinp,inp,ipoinpc)
183  if(textpart(1)(2:10).eq.'KINEMATIC') then
184 !
185 ! catalogueing the nodes
186 !
187  npt=0
188 !
189  do j=istartset(jsurf),iendset(jsurf)
190  if(ialset(j).gt.0) then
191  if(surfkind.eq.'T') then
192 !
193 ! facial surface
194 !
195  iface=ialset(j)
196  nelem=int(iface/10)
197  jface=iface-nelem*10
198  indexe=ipkon(nelem)
199 !
200  if(lakon(nelem)(4:5).eq.'20') then
201  nopes=8
202  elseif(lakon(nelem)(4:4).eq.'2') then
203  nopes=9
204  elseif(lakon(nelem)(4:4).eq.'8') then
205  nopes=4
206  elseif(lakon(nelem)(4:5).eq.'10') then
207  nopes=6
208  elseif(lakon(nelem)(4:4).eq.'4') then
209  nopes=3
210  endif
211 !
212  if(lakon(nelem)(4:4).eq.'6') then
213  if(jface.le.2) then
214  nopes=3
215  else
216  nopes=4
217  endif
218  endif
219  if(lakon(nelem)(4:5).eq.'15') then
220  if(jface.le.2) then
221  nopes=6
222  else
223  nopes=8
224  endif
225  endif
226  else
227 !
228 ! nodal surface
229 !
230  nopes=1
231  endif
232 !
233  do m=1,nopes
234  if(surfkind.eq.'T') then
235  if((lakon(nelem)(4:4).eq.'2').or.
236  & (lakon(nelem)(4:4).eq.'8')) then
237  node=kon(indexe+ifaceq(m,jface))
238  elseif((lakon(nelem)(4:4).eq.'4').or.
239  & (lakon(nelem)(4:5).eq.'10')) then
240  node=kon(indexe+ifacet(m,jface))
241  elseif(lakon(nelem)(4:4).eq.'6') then
242  node=kon(indexe+ifacew1(m,jface))
243  elseif(lakon(nelem)(4:5).eq.'15') then
244  node=kon(indexe+ifacew2(m,jface))
245  endif
246  else
247  node =ialset(j)
248  endif
249 !
250  call nident2(ics,node,npt,id)
251  if(id.gt.0) then
252  if(ics(1,id).eq.node) then
253  cycle
254  endif
255  endif
256 !
257 ! updating ics
258 !
259  npt=npt+1
260  do l=npt,id+2,-1
261  ics(1,l)=ics(1,l-1)
262  enddo
263  ics(1,id+1)=node
264  enddo
265  else
266 !
267 ! if a negative value occurs the surface has to be
268 ! nodal
269 !
270  k=ialset(j-2)
271  do
272  k=k-ialset(j)
273  if(k.ge.ialset(j-1)) exit
274  node=k
275  call nident2(ics,node,npt,id)
276  if(id.gt.0) then
277  if(ics(1,id).eq.node) then
278  cycle
279  endif
280  endif
281 !
282 ! updating ics
283 !
284  npt=npt+1
285  do l=npt,id+2,-1
286  ics(1,l)=ics(1,l-1)
287  enddo
288  ics(1,id+1)=node
289  enddo
290  endif
291  enddo
292 !
293 ! generating a rotational node and connecting the
294 ! rotational dofs of the reference node with the
295 ! translational dofs of the rotational node
296 !
297 ! generating a rotational reference node
298 !
299  nk=nk+1
300  if(nk.gt.nk_) then
301  write(*,*)
302  & '*ERROR reading *KINEMATIC: increase nk_'
303  call exit(201)
304  endif
305  irotnode_kin=nk
306  do l=1,3
307  co(l,nk)=co(l,irefnode)
308  enddo
309 !
310 ! generating connecting MPCs between the rotational
311 ! dofs of irefnode and the translational dofs of
312 ! irotnode_kin
313 !
314  do k=1,3
315 !
316  nmpc=nmpc+1
317  if(nmpc.gt.nmpc_) then
318  write(*,*)
319  & '*ERROR reading *KINEMATIC: increase nmpc_'
320  call exit(201)
321  endif
322 !
323 ! the internal dofs for rotation are 4, 5 and 6
324 !
325  ipompc(nmpc)=mpcfree
326  labmpc(nmpc)='ROTTRACOUPLING '
327  idof=8*(irefnode-1)+k+3
328  call nident(ikmpc,idof,nmpc-1,id)
329  do l=nmpc,id+2,-1
330  ikmpc(l)=ikmpc(l-1)
331  ilmpc(l)=ilmpc(l-1)
332  enddo
333  ikmpc(id+1)=idof
334  ilmpc(id+1)=nmpc
335 !
336  nodempc(1,mpcfree)=irefnode
337  nodempc(2,mpcfree)=k+4
338  coefmpc(mpcfree)=1.d0
339  mpcfree=nodempc(3,mpcfree)
340 !
341  nodempc(1,mpcfree)=irotnode_kin
342  nodempc(2,mpcfree)=k
343  coefmpc(mpcfree)=-1.d0
344  mpcfreeold=mpcfree
345  mpcfree=nodempc(3,mpcfree)
346  nodempc(3,mpcfreeold)=0
347  enddo
348 !
349  if(iorientation.gt.0) then
350 !
351 ! duplicating the nodes
352 ! generating rigid body MPC's for all dofs in the new nodes
353 !
354  do m=1,npt
355  nk=nk+1
356  if(nk.gt.nk_) then
357  write(*,*)
358  & '*ERROR reading *KINEMATIC: increase nk_'
359  call exit(201)
360  endif
361  ics(2,m)=nk
362  do k=1,3
363  co(k,nk)=co(k,ics(1,m))
364  enddo
365  node=nk
366  ibounstart=1
367  ibounend=3
368  call rigidmpc(ipompc,nodempc,coefmpc,irefnode,
369  & irotnode_kin,
370  & labmpc,nmpc,nmpc_,mpcfree,ikmpc,ilmpc,nk,nk_,
371  & nodeboun,ndirboun,ikboun,ilboun,nboun,nboun_,node,
372  & typeboun,co,ibounstart,ibounend)
373  enddo
374  endif
375 !
376 ! reading the degrees of freedom
377 !
378  do
379  call getnewline(inpc,textpart,istat,n,key,iline,ipol,inl,
380  & ipoinp,inp,ipoinpc)
381  if((istat.lt.0).or.(key.eq.1)) return
382 !
383  read(textpart(1)(1:10),'(i10)',iostat=istat) ibounstart
384  if(istat.gt.0) call inputerror(inpc,ipoinpc,iline,
385  &"*KINEMATIC%")
386  if(ibounstart.lt.1) then
387  write(*,*) '*ERROR reading *KINEMATIC'
388  write(*,*) ' starting degree of freedom cannot'
389  write(*,*) ' be less than 1'
390  write(*,*) ' '
391  call inputerror(inpc,ipoinpc,iline,
392  &"*KINEMATIC%")
393  call exit(201)
394  endif
395 !
396  if(textpart(2)(1:1).eq.' ') then
397  ibounend=ibounstart
398  else
399  read(textpart(2)(1:10),'(i10)',iostat=istat) ibounend
400  if(istat.gt.0) call inputerror(inpc,ipoinpc,iline,
401  &"*BOUNDARY%")
402  endif
403  if(ibounend.gt.3) then
404  write(*,*) '*ERROR reading *KINEMATIC'
405  write(*,*) ' final degree of freedom cannot'
406  write(*,*) ' exceed 3'
407  write(*,*) ' '
408  call inputerror(inpc,ipoinpc,iline,
409  &"*KINEMATIC%")
410  call exit(201)
411  elseif(ibounend.lt.ibounstart) then
412  write(*,*) '*ERROR reading *KINEMATIC'
413  write(*,*) ' initial degree of freedom cannot'
414  write(*,*) ' exceed final degree of freedom'
415  write(*,*) ' '
416  call inputerror(inpc,ipoinpc,iline,
417  &"*KINEMATIC%")
418  call exit(201)
419  endif
420 !
421 ! generating the MPCs
422 !
423  if(iorientation.eq.0) then
424 !
425 ! generating rigid body MPC's for the appropriate dofs
426 !
427  do j=1,npt
428  node=ics(1,j)
429  call rigidmpc(ipompc,nodempc,coefmpc,irefnode,
430  & irotnode_kin,
431  & labmpc,nmpc,nmpc_,mpcfree,ikmpc,ilmpc,nk,nk_,
432  & nodeboun,ndirboun,ikboun,ilboun,nboun,nboun_,
433  & node,typeboun,co,ibounstart,ibounend)
434  enddo
435  else
436 !
437 ! connecting the original nodes with the duplicated nodes
438 ! for the appropriate dofs
439 !
440  do j=1,npt
441  node=ics(1,j)
442  idupnode=ics(2,j)
443  call mpcadd(node,ibounstart,ibounend,nboun,ipompc,
444  & nodempc,coefmpc,nmpc,nmpc_,mpcfree,orab,ikboun,
445  & ikmpc,ilmpc,co,labmpc,label,idupnode,
446  & iorientation)
447  enddo
448  endif
449  enddo
450  elseif(textpart(1)(2:13).eq.'DISTRIBUTING') then
451  if(surfkind.eq.'S') then
452  write(*,*) '*ERROR reading *DISTRIBUTING'
453  write(*,*) ' a nodal surface is not allowed'
454  write(*,*) ' please use a facial surface on'
455  write(*,*) ' the *COUPLING card'
456  call exit(201)
457  endif
458 !
459  npt=0
460  area=0.d0
461 !
462 ! catalogueing the nodes belonging to the surface (ics(1,*))
463 ! catalogueing the area (dcs(*))
464 !
465  do k=istartset(jsurf),iendset(jsurf)
466 !
467 ! facial surface
468 !
469  iface=ialset(k)
470  nelem=int(iface/10)
471  jface=iface-nelem*10
472  indexe=ipkon(nelem)
473 !
474 ! nodes: #nodes in the face
475 ! the nodes are stored in nodef(*)
476 !
477  if(lakon(nelem)(4:4).eq.'2') then
478  nopes=8
479  nface=6
480  elseif(lakon(nelem)(3:4).eq.'D8') then
481  nopes=4
482  nface=6
483  elseif(lakon(nelem)(4:5).eq.'10') then
484  nopes=6
485  nface=4
486  nope=10
487  elseif(lakon(nelem)(4:4).eq.'4') then
488  nopes=3
489  nface=4
490  nope=4
491  elseif(lakon(nelem)(4:5).eq.'15') then
492  if(jface.le.2) then
493  nopes=6
494  else
495  nopes=8
496  endif
497  nface=5
498  nope=15
499  elseif(lakon(nelem)(3:4).eq.'D6') then
500  if(jface.le.2) then
501  nopes=3
502  else
503  nopes=4
504  endif
505  nface=5
506  nope=6
507  else
508  cycle
509  endif
510 !
511 ! determining the nodes of the face
512 !
513  if(nface.eq.4) then
514  do i=1,nopes
515  nodef(i)=kon(indexe+ifacet(i,jface))
516  enddo
517  elseif(nface.eq.5) then
518  if(nope.eq.6) then
519  do i=1,nopes
520  nodef(i)=kon(indexe+ifacew1(i,jface))
521  enddo
522  elseif(nope.eq.15) then
523  do i=1,nopes
524  nodef(i)=kon(indexe+ifacew2(i,jface))
525  enddo
526  endif
527  elseif(nface.eq.6) then
528  do i=1,nopes
529  nodef(i)=kon(indexe+ifaceq(i,jface))
530  enddo
531  endif
532 !
533 ! loop over the nodes belonging to the face
534 ! ics(1,*): pretension node
535 ! dcs(*): area corresponding to pretension node
536 !
537  do i=1,nopes
538  node=nodef(i)
539  call nident2(ics,node,npt,id)
540  if(id.gt.0) then
541  if(ics(1,id).eq.node) then
542  cycle
543  endif
544  endif
545 !
546 ! updating ics
547 !
548  npt=npt+1
549  do j=npt,id+2,-1
550  ics(1,j)=ics(1,j-1)
551  dcs(j)=dcs(j-1)
552  enddo
553  ics(1,id+1)=node
554  dcs(id+1)=0.d0
555  enddo
556 !
557 ! calculating the area of the face and its contributions
558 ! to the facial nodes
559 !
560 ! number of integration points
561 !
562  if(lakon(nelem)(3:5).eq.'D8R') then
563  mint=1
564  elseif(lakon(nelem)(3:4).eq.'D8') then
565  mint=4
566  elseif(lakon(nelem)(4:6).eq.'20R') then
567  mint=4
568  elseif(lakon(nelem)(4:4).eq.'2') then
569  mint=9
570  elseif(lakon(nelem)(4:5).eq.'10') then
571  mint=3
572  elseif(lakon(nelem)(4:4).eq.'4') then
573  mint=1
574  elseif(lakon(nelem)(3:4).eq.'D6') then
575  mint=1
576  elseif(lakon(nelem)(4:5).eq.'15') then
577  if(jface.le.2) then
578  mint=3
579  else
580  mint=4
581  endif
582  endif
583 !
584  do i=1,nopes
585  areanodal(i)=0.d0
586  do j=1,3
587  xl2(j,i)=co(j,nodef(i))
588  enddo
589  enddo
590 !
591  do m=1,mint
592  if((lakon(nelem)(3:5).eq.'D8R').or.
593  & ((lakon(nelem)(3:4).eq.'D6').and.(nopes.eq.4))) then
594  xi=gauss2d1(1,m)
595  et=gauss2d1(2,m)
596  weight=weight2d1(m)
597  elseif((lakon(nelem)(3:4).eq.'D8').or.
598  & (lakon(nelem)(4:6).eq.'20R').or.
599  & ((lakon(nelem)(4:5).eq.'15').and.
600  & (nopes.eq.8))) then
601  xi=gauss2d2(1,m)
602  et=gauss2d2(2,m)
603  weight=weight2d2(m)
604  elseif(lakon(nelem)(4:4).eq.'2') then
605  xi=gauss2d3(1,m)
606  et=gauss2d3(2,m)
607  weight=weight2d3(m)
608  elseif((lakon(nelem)(4:5).eq.'10').or.
609  & ((lakon(nelem)(4:5).eq.'15').and.
610  & (nopes.eq.6))) then
611  xi=gauss2d5(1,m)
612  et=gauss2d5(2,m)
613  weight=weight2d5(m)
614  elseif((lakon(nelem)(4:4).eq.'4').or.
615  & ((lakon(nelem)(3:4).eq.'D6').and.
616  & (nopes.eq.3))) then
617  xi=gauss2d4(1,m)
618  et=gauss2d4(2,m)
619  weight=weight2d4(m)
620  endif
621 !
622  if(nopes.eq.8) then
623  call shape8q(xi,et,xl2,xsj2,xs2,shp2,iflag)
624  elseif(nopes.eq.4) then
625  call shape4q(xi,et,xl2,xsj2,xs2,shp2,iflag)
626  elseif(nopes.eq.6) then
627  call shape6tri(xi,et,xl2,xsj2,xs2,shp2,iflag)
628  elseif(nopes.eq.3) then
629  call shape3tri(xi,et,xl2,xsj2,xs2,shp2,iflag)
630  endif
631 !
632 ! calculating the total area and nodal area
633 !
634  xsj=weight*dsqrt(xsj2(1)**2+xsj2(2)**2+xsj2(3)**2)
635  area=area+xsj
636  do i=1,nopes
637  areanodal(i)=areanodal(i)+xsj*shp2(4,i)
638  enddo
639 !
640  enddo
641 !
642 ! inserting the nodal area into field dcs
643 !
644  do i=1,nopes
645  node=nodef(i)
646  call nident2(ics,node,npt,id)
647  dcs(id)=dcs(id)+areanodal(i)
648  enddo
649 !
650  enddo
651 !
652 ! create for each node in the surface a new node.
653 ! the ratio of the displacements between both nodes
654 ! is governed by the area for which each node is
655 ! representative
656 !
657 ! initializing the location of the center of gravity
658 !
659  do k=1,3
660  cgx(k)=0.d0
661  enddo
662 !
663  do m=1,npt
664  nk=nk+1
665  if(nk.gt.nk_) then
666  write(*,*)
667  & '*ERROR reading *DISTRIBUTING: increase nk_'
668  call exit(201)
669  endif
670  node=ics(1,m)
671  ics(1,m)=nk
672  do k=1,3
673  co(k,nk)=co(k,node)
674  enddo
675 !
676 ! generate the connecting MPC's
677 !
678  do k=1,3
679 !
680 ! contribution to the location of the center of gravity
681 !
682  cgx(k)=cgx(k)+dcs(m)*co(k,node)
683 !
684  nmpc=nmpc+1
685  if(nmpc.gt.nmpc_) then
686  write(*,*)
687  & '*ERROR reading *DISTRIBUTING: increase nmpc_'
688  call exit(201)
689  endif
690 !
691 ! MPC: u(old node)=
692 ! (mean area)/(area_m) * u(new node) (in all directions)
693 !
694 ! check whether MPC was already used
695 !
696  ipompc(nmpc)=mpcfree
697  labmpc(nmpc)=' '
698  idof=8*(node-1)+k
699  call nident(ikmpc,idof,nmpc-1,id)
700  if(id.gt.0) then
701  if(ikmpc(id).eq.idof) then
702  write(*,*) '*ERROR in couplings: dof',k
703  write(*,*) ' in node ',node
704  write(*,*) ' was already used'
705  call exit(201)
706  endif
707  endif
708  do l=nmpc,id+2,-1
709  ikmpc(l)=ikmpc(l-1)
710  ilmpc(l)=ilmpc(l-1)
711  enddo
712  ikmpc(id+1)=idof
713  ilmpc(id+1)=nmpc
714 !
715 ! generating the terms of the MPC
716 !
717  nodempc(1,mpcfree)=node
718  nodempc(2,mpcfree)=k
719  coefmpc(mpcfree)=1.d0
720  mpcfree=nodempc(3,mpcfree)
721 !
722  nodempc(1,mpcfree)=nk
723  nodempc(2,mpcfree)=k
724  coefmpc(mpcfree)=-area/(npt*dcs(m))
725  mpcfreeold=mpcfree
726  mpcfree=nodempc(3,mpcfree)
727  nodempc(3,mpcfreeold)=0
728  enddo
729  enddo
730 !
731  do k=1,3
732  cgx(k)=cgx(k)/area
733  enddo
734 !
735 ! generating the translational MPC's (default)
736 !
737  do m=1,3
738  node=ics(1,1)
739  idof=8*(node-1)+m
740  call nident(ikmpc,idof,nmpc,id)
741 !
742  nmpc=nmpc+1
743  if(nmpc.gt.nmpc_) then
744  write(*,*) '*ERROR in equations: increase nmpc_'
745  call exit(201)
746  endif
747  labmpc(nmpc)=' '
748  ipompc(nmpc)=mpcfree
749 !
750 ! updating ikmpc and ilmpc
751 !
752  do j=nmpc,id+2,-1
753  ikmpc(j)=ikmpc(j-1)
754  ilmpc(j)=ilmpc(j-1)
755  enddo
756  ikmpc(id+1)=idof
757  ilmpc(id+1)=nmpc
758 !
759  do j=1,npt
760  node=ics(1,j)
761  nodempc(1,mpcfree)=node
762  nodempc(2,mpcfree)=m
763  coefmpc(mpcfree)=1.d0
764  mpcfree=nodempc(3,mpcfree)
765  enddo
766  nodempc(1,mpcfree)=irefnode
767  nodempc(2,mpcfree)=m
768  coefmpc(mpcfree)=-npt
769  mpcfreeold=mpcfree
770  mpcfree=nodempc(3,mpcfree)
771  nodempc(3,mpcfreeold)=0
772  enddo
773 !
774 ! generating the rotational MPC's
775 !
776  irotnode(1)=0
777 !
778 ! reading the degrees of freedom
779 !
780  do
781  call getnewline(inpc,textpart,istat,n,key,iline,ipol,inl,
782  & ipoinp,inp,ipoinpc)
783  if((istat.lt.0).or.(key.eq.1)) return
784 !
785  read(textpart(1)(1:10),'(i10)',iostat=istat) ibounstart
786  if(istat.gt.0) call inputerror(inpc,ipoinpc,iline,
787  &"*KINEMATIC%")
788  if(ibounstart.lt.1) then
789  write(*,*) '*ERROR reading *KINEMATIC'
790  write(*,*) ' starting degree of freedom cannot'
791  write(*,*) ' be less than 1'
792  write(*,*) ' '
793  call inputerror(inpc,ipoinpc,iline,
794  &"*KINEMATIC%")
795  call exit(201)
796  endif
797 !
798  if(textpart(2)(1:1).eq.' ') then
799  ibounend=ibounstart
800  else
801  read(textpart(2)(1:10),'(i10)',iostat=istat) ibounend
802  if(istat.gt.0) call inputerror(inpc,ipoinpc,iline,
803  &"*BOUNDARY%")
804  endif
805 !
806  ibounstart=max(4,ibounstart)
807  if(ibounend.gt.6) then
808  write(*,*) '*ERROR reading *DISTRIBUTING'
809  write(*,*) ' final degree of freedom cannot'
810  write(*,*) ' exceed 6'
811  write(*,*) ' '
812  call inputerror(inpc,ipoinpc,iline,
813  &"*DISTRIBUTING%")
814  call exit(201)
815  elseif(ibounend.lt.ibounstart) then
816  cycle
817  endif
818 !
819  if((ibounend.gt.3).and.(irotnode(1).eq.0)) then
820 !
821 ! check the orientation definition: the local reference
822 ! system is not allowed to be cylindrical
823 !
824  if(iorientation.ne.0) then
825  if(orab(7,iorientation).lt.0.d0) then
826  write(*,*) '*ERROR reading *DISTRIBUTING'
827  write(*,*) ' a cylindrical local coordinate'
828  write(*,*) ' system is not allowed'
829  call exit(201)
830  endif
831 !
832  call transformatrix(orab(1,iorientation),
833  & co(1,irefnode),a)
834  endif
835 !
836 ! generating a rotational reference node
837 !
838  do k=1,3
839  nk=nk+1
840  if(nk.gt.nk_) then
841  write(*,*)
842  & '*ERROR reading *DISTRIBUTING: increase nk_'
843  call exit(201)
844  endif
845  irotnode(k)=nk
846  do l=1,3
847  co(l,nk)=co(l,irefnode)
848  enddo
849  enddo
850 !
851 ! generating connecting MPCs between the rotational
852 ! dofs of irefnode and the translational dofs of
853 ! irotnode
854 !
855  do k=1,3
856 !
857  nmpc=nmpc+1
858  if(nmpc.gt.nmpc_) then
859  write(*,*)
860  & '*ERROR reading *DISTRIBUTING: increase nmpc_'
861  call exit(201)
862  endif
863 !
864 ! the internal dofs for rotation are 4, 5 and 7
865 !
866  ipompc(nmpc)=mpcfree
867  labmpc(nmpc)='ROTTRACOUPLING '
868  idof=8*(irefnode-1)+k+3
869  call nident(ikmpc,idof,nmpc-1,id)
870  do l=nmpc,id+2,-1
871  ikmpc(l)=ikmpc(l-1)
872  ilmpc(l)=ilmpc(l-1)
873  enddo
874  ikmpc(id+1)=idof
875  ilmpc(id+1)=nmpc
876 !
877  nodempc(1,mpcfree)=irefnode
878  nodempc(2,mpcfree)=k+4
879  coefmpc(mpcfree)=1.d0
880  mpcfree=nodempc(3,mpcfree)
881 !
882  nodempc(1,mpcfree)=irotnode(k)
883  nodempc(2,mpcfree)=1
884  coefmpc(mpcfree)=-1.d0
885  mpcfreeold=mpcfree
886  mpcfree=nodempc(3,mpcfree)
887  nodempc(3,mpcfreeold)=0
888  enddo
889 !
890 ! generating a node for the inhomogeneous term
891 !
892  nk=nk+1
893  if(nk.gt.nk_) then
894  write(*,*)
895  & '*ERROR reading *DISTRIBUTING: increase nk_'
896  call exit(201)
897  endif
898  inhomnode=nk
899 !
900 ! fictituous center of gravity (for compatibility
901 ! reasons with usermpc.f)
902 !
903 c do k=1,3
904 c cgx(k)=co(k,irefnode)
905 c enddo
906  endif
907 !
908 ! generating rotational MPC's
909 !
910  do kk=ibounstart,ibounend
911  k=kk-3
912 !
913 ! axis of rotation
914 !
915  if(iorientation.eq.0) then
916  do j=1,3
917  aa(j)=0.d0
918  enddo
919  aa(k)=1.d0
920  else
921  do j=1,3
922  aa(j)=a(j,k)
923  enddo
924  endif
925 !
926 ! storing the axis of rotation as coordinates of the
927 ! rotation nodes
928 !
929  do j=1,3
930  co(j,irotnode(k))=aa(j)
931  enddo
932 !
933  nmpc=nmpc+1
934  if(nmpc.gt.nmpc_) then
935  write(*,*)
936  & '*ERROR reading *DISTRIBUTING: increase nmpc_'
937  call exit(201)
938  endif
939 !
940  ipompc(nmpc)=mpcfree
941  labmpc(nmpc)='MEANROTDISTRIB '
942 !
943 ! defining the terms of the MPC
944 !
945  do m=1,npt
946  do j=1,3
947  nodempc(1,mpcfree)=ics(1,m)
948  nodempc(2,mpcfree)=0
949  coefmpc(mpcfree)=0.d0
950  mpcfree=nodempc(3,mpcfree)
951  enddo
952  enddo
953  nodempc(1,mpcfree)=irotnode(k)
954  nodempc(2,mpcfree)=1
955  coefmpc(mpcfree)=0.d0
956  mpcfree=nodempc(3,mpcfree)
957  nodempc(1,mpcfree)=inhomnode
958  nodempc(2,mpcfree)=k
959  coefmpc(mpcfree)=0.d0
960  mpcfreeold=mpcfree
961  mpcfree=nodempc(3,mpcfree)
962  nodempc(3,mpcfreeold)=0
963 !
964 ! storing the inhomogeneous value as a boundary
965 ! condition in the inhomogeneous node
966 !
967  idof=8*(inhomnode-1)+k
968  call nident(ikboun,idof,nboun,id)
969  nboun=nboun+1
970  if(nboun.gt.nboun_) then
971  write(*,*) '*ERROR in usermpc: increase nboun_'
972  call exit(201)
973  endif
974  nodeboun(nboun)=inhomnode
975  ndirboun(nboun)=k
976  xboun(nboun)=0.d0
977  typeboun(nboun)='U'
978  do l=nboun,id+2,-1
979  ikboun(l)=ikboun(l-1)
980  ilboun(l)=ilboun(l-1)
981  enddo
982  ikboun(id+1)=idof
983  ilboun(id+1)=nboun
984 !
985 ! calculating the coefficients of the rotational MPC
986 !
987 ! loop over all nodes belonging to the mean rotation MPC
988 !
989  index2=ipompc(nmpc)
990  do
991  node=nodempc(1,index2)
992  if(node.eq.irotnode(k)) exit
993 !
994 ! relative positions
995 !
996  do j=1,3
997  pi(j)=co(j,node)-cgx(j)
998  enddo
999 !
1000 ! projection on a plane orthogonal to the rotation vector
1001 !
1002  c1=pi(1)*aa(1)+pi(2)*aa(2)+pi(3)*aa(3)
1003  do j=1,3
1004  pi(j)=pi(j)-c1*aa(j)
1005  enddo
1006 !
1007  c1=pi(1)*pi(1)+pi(2)*pi(2)+pi(3)*pi(3)
1008  if(c1.lt.1.d-20) then
1009  write(*,*)
1010  & '*WARNING reading *DISTRIBUTING: node ',node
1011  write(*,*) ' is very close to the '
1012  write(*,*) ' rotation axis through the'
1013  write(*,*) ' center of gravity of'
1014  write(*,*) ' the nodal cloud in a'
1015  write(*,*) ' mean rotation MPC.'
1016  write(*,*) ' This node is not taken'
1017  write(*,*) ' into account in the MPC'
1018  index2=nodempc(3,nodempc(3,nodempc(3,index2)))
1019  cycle
1020  endif
1021 !
1022  do j=1,3
1023  if(j.eq.1) then
1024  c4=aa(2)*pi(3)-aa(3)*pi(2)
1025  elseif(j.eq.2) then
1026  c4=aa(3)*pi(1)-aa(1)*pi(3)
1027  else
1028  c4=aa(1)*pi(2)-aa(2)*pi(1)
1029  endif
1030  c9=c4/c1
1031 !
1032  index1=ipompc(nmpc)
1033  do
1034  node1=nodempc(1,index1)
1035  if(node1.eq.irotnode(k)) exit
1036  if(node1.eq.node) then
1037  c10=c9*(1.d0-1.d0/real(npt))
1038  else
1039  c10=-c9/real(npt)
1040  endif
1041  if(j.eq.1) then
1042  coefmpc(index1)=coefmpc(index1)+c10
1043  elseif(j.eq.2) then
1044  coefmpc(nodempc(3,index1))=
1045  & coefmpc(nodempc(3,index1))+c10
1046  else
1047  coefmpc(nodempc(3,nodempc(3,index1)))=
1048  & coefmpc(nodempc(3,nodempc(3,index1)))+c10
1049  endif
1050  index1=
1051  & nodempc(3,nodempc(3,nodempc(3,index1)))
1052  enddo
1053  enddo
1054  index2=nodempc(3,nodempc(3,nodempc(3,index2)))
1055  enddo
1056  coefmpc(index2)=-npt
1057 !
1058 ! assigning the degrees of freedom
1059 !
1060  j=0
1061  index2=ipompc(nmpc)
1062  do
1063  j=j+1
1064  if(j.gt.3) j=1
1065  nodempc(2,index2)=j
1066  index2=nodempc(3,index2)
1067  if(nodempc(1,index2).eq.irotnode(k)) exit
1068  enddo
1069 !
1070 ! look for biggest coefficient of all but the last
1071 ! regular node. The general form of the MPC is:
1072 ! regular nodes, rotationl dof and inhomogeneous dof
1073 !
1074  indexmax=0
1075  index2=ipompc(nmpc)
1076  amax=1.d-5
1077 !
1078 ! coefficients of the first terms in the translational
1079 ! dofs
1080 !
1081  coef(1)=coefmpc(index2)
1082  coef(2)=coefmpc(nodempc(3,index2))
1083  coef(3)=coefmpc(nodempc(3,nodempc(3,index2)))
1084 !
1085 ! loop over all regular nodes
1086 !
1087  do i=1,3*npt
1088  if(dabs(coefmpc(index2)).gt.amax) then
1089  idir=nodempc(2,index2)
1090 !
1091  if(dabs(coefmpc(index2)-coef(idir)).lt.1.d-10) then
1092  index2=nodempc(3,index2)
1093  cycle
1094  endif
1095 !
1096  node=nodempc(1,index2)
1097  idof=8*(node-1)+idir
1098 !
1099 ! check whether the node was used in another MPC
1100 !
1101  call nident(ikmpc,idof,nmpc-1,id)
1102  if(id.gt.0) then
1103  if(ikmpc(id).eq.idof) then
1104  index2=nodempc(3,index2)
1105  cycle
1106  endif
1107  endif
1108 !
1109 ! check whether the node was used in a SPC
1110 !
1111  call nident(ikboun,idof,nboun,id)
1112  if(id.gt.0) then
1113  if(ikboun(id).eq.idof) then
1114  index2=nodempc(3,index2)
1115  cycle
1116  endif
1117  endif
1118 !
1119  amax=dabs(coefmpc(index2))
1120  indexmax=index2
1121  endif
1122 !
1123 ! after each node has been treated, a check is performed
1124 ! whether indexmax is already nonzero. Taking too many
1125 ! nodes into account may lead to dependent equations if
1126 ! all rotational dofs are constrained.
1127 !
1128  if((i/3)*3.eq.i) then
1129  if(indexmax.ne.0) exit
1130  endif
1131 !
1132  index2=nodempc(3,index2)
1133  enddo
1134 !
1135  if(indexmax.eq.0) then
1136  write(*,*)
1137  & '*WARNING reading *DISTRIBUTING: no MPC is '
1138  write(*,*) ' generated for the mean'
1139  write(*,*) ' rotation in node ',irefnode
1140  write(*,*) ' and direction ',kk
1141 !
1142 ! removing the MPC
1143 !
1144  mpcfreeold=mpcfree
1145  index2=ipompc(nmpc)
1146  ipompc(nmpc)=0
1147  mpcfree=index2
1148 !
1149 ! removing the MPC from fields nodempc and coefmpc
1150 !
1151  do
1152  nodempc(1,index2)=0
1153  nodempc(2,index2)=0
1154  coefmpc(index2)=0
1155  if(nodempc(3,index2).ne.0) then
1156  index2=nodempc(3,index2)
1157  else
1158  nodempc(3,index2)=mpcfreeold
1159  exit
1160  endif
1161  enddo
1162 !
1163 ! decrementing nmpc
1164 !
1165  nmpc=nmpc-1
1166 !
1167  cycle
1168  endif
1169 !
1170  nodemax=nodempc(1,indexmax)
1171  idirmax=nodempc(2,indexmax)
1172  index2=ipompc(nmpc)
1173 !
1174  nodeold=nodempc(1,index2)
1175  idirold=nodempc(2,index2)
1176 !
1177 ! exchange the node information in the MPC
1178 !
1179  if(nodemax.ne.nodeold) then
1180 !
1181  indexold(1)=index2
1182  indexold(2)=nodempc(3,index2)
1183  indexold(3)=nodempc(3,indexold(2))
1184 !
1185  index2=nodempc(3,indexold(3))
1186  do
1187  if(nodempc(1,index2).eq.irotnode(k)) exit
1188  if(nodempc(1,index2).eq.nodemax) then
1189  if(nodempc(2,index2).eq.1) then
1190  indexnew(1)=index2
1191  elseif(nodempc(2,index2).eq.2) then
1192  indexnew(2)=index2
1193  else
1194  indexnew(3)=index2
1195  endif
1196  endif
1197  index2=nodempc(3,index2)
1198  enddo
1199 !
1200  do j=1,3
1201  node=nodempc(1,indexold(j))
1202  idir=nodempc(2,indexold(j))
1203  xcoef=coefmpc(indexold(j))
1204  nodempc(1,indexold(j))=nodempc(1,indexnew(j))
1205  nodempc(2,indexold(j))=nodempc(2,indexnew(j))
1206  coefmpc(indexold(j))=coefmpc(indexnew(j))
1207  nodempc(1,indexnew(j))=node
1208  nodempc(2,indexnew(j))=idir
1209  coefmpc(indexnew(j))=xcoef
1210  enddo
1211  endif
1212 !
1213 ! exchange the direction information in the MPC
1214 !
1215  index2=ipompc(nmpc)
1216  if(idirmax.ne.1) then
1217  indexold(1)=index2
1218  if(idirmax.eq.2) then
1219  indexnew(1)=nodempc(3,index2)
1220  else
1221  indexnew(1)=nodempc(3,nodempc(3,index2))
1222  endif
1223 !
1224  do j=1,1
1225  idir=nodempc(2,indexold(j))
1226  xcoef=coefmpc(indexold(j))
1227  nodempc(2,indexold(j))=nodempc(2,indexnew(j))
1228  coefmpc(indexold(j))=coefmpc(indexnew(j))
1229  nodempc(2,indexnew(j))=idir
1230  coefmpc(indexnew(j))=xcoef
1231  enddo
1232  endif
1233 !
1234 ! determining the dependent dof of the MPC
1235 !
1236  idof=8*(nodemax-1)+idirmax
1237  call nident(ikmpc,idof,nmpc-1,id)
1238  do l=nmpc,id+2,-1
1239  ikmpc(l)=ikmpc(l-1)
1240  ilmpc(l)=ilmpc(l-1)
1241  enddo
1242  ikmpc(id+1)=idof
1243  ilmpc(id+1)=nmpc
1244 !
1245  enddo
1246  enddo
1247  else
1248  write(*,*)
1249  & '*ERROR reading *COUPLING: the line following'
1250  write(*,*) ' *COUPLING must contain the'
1251  write(*,*) ' *KINEMATIC keyword or the'
1252  write(*,*) ' *DISTRIBUTING keyword'
1253  write(*,*) ' '
1254  call inputerror(inpc,ipoinpc,iline,
1255  &"*COUPLING%")
1256  call exit(201)
1257  endif
1258 !
1259  return
shape8q
subroutine shape8q(xi, et, xl, xsj, xs, shp, iflag)
Definition: shape8q.f:20
max
#define max(a, b)
Definition: cascade.c:32
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
mpcadd
subroutine mpcadd(nodedep, is, ie, nboun, ipompc, nodempc, coefmpc, nmpc, nmpc_, mpcfree, orab, ikboun, ikmpc, ilmpc, co, labmpc, label, nodeind, iorientation)
Definition: mpcadd.f:22
c1
static double * c1
Definition: mafillvcompmain.c:30
transformatrix
subroutine transformatrix(xab, p, a)
Definition: transformatrix.f:20
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
rigidmpc
subroutine rigidmpc(ipompc, nodempc, coefmpc, irefnode, irotnode, labmpc, nmpc, nmpc_, mpcfree, ikmpc, ilmpc, nk, nk_, nodeboun, ndirboun, ikboun, ilboun, nboun, nboun_, node, typeboun, co, jmin, jmax)
Definition: rigidmpc.f:22
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)