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