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 rigidbodys(inpc,textpart,set,istartset,iendset,
20     &     ialset,nset,nset_,nalset,nalset_,ipompc,nodempc,coefmpc,
21     &     labmpc,nmpc,nmpc_,mpcfree,ikmpc,ilmpc,lakon,ipkon,kon,nk,nk_,
22     &     nodeboun,ndirboun,ikboun,ilboun,nboun,nboun_,iperturb,ne_,
23     &     ctrl,typeboun,istep,istat,n,iline,ipol,inl,ipoinp,inp,co,
24     &     ipoinpc,ier)
25!
26!     reading the input deck: *RIGID BODY
27!
28      implicit none
29!
30      character*1 typeboun(*),inpc(*)
31      character*8 lakon(*)
32      character*20 labmpc(*)
33      character*81 set(*),elset,noset
34      character*132 textpart(16)
35!
36      integer istartset(*),iendset(*),ialset(*),ipompc(*),
37     &     nodempc(3,*),
38     &     nset,nset_,nalset,nalset_,nmpc,nmpc_,mpcfree,nk,nk_,ikmpc(*),
39     &     ilmpc(*),ipkon(*),kon(*),inoset,ielset,i,node,ielement,id,
40     &     indexe,nope,istep,istat,n,irefnode,irotnode,ne_,
41     &     j,idof,k,nodeboun(*),ndirboun(*),ikboun(*),ilboun(*),
42     &     nboun,nboun_,key,iperturb(*),ipos,iline,ipol,inl,ipoinp(2,*),
43     &     inp(3,*),ipoinpc(0:*),jmin,jmax,ier
44!
45      real*8 coefmpc(3,*),ctrl(*),co(3,*)
46!
47      jmin=1
48      jmax=3
49!
50      if(istep.gt.0) then
51        write(*,*)
52     &       '*ERROR reading *RIGID BODY: *RIGID BODY should be placed'
53        write(*,*) '  before all step definitions'
54        ier=1
55        return
56      endif
57!
58!     the *RIGID BODY option implies a nonlinear geometric
59!     calculation
60!
61      if(iperturb(1).eq.1) then
62        write(*,*) '*ERROR reading *RIGID BODY: the *RIGID BODY option'
63        write(*,*) '       cannot be used in a perturbation step'
64        ier=1
65        return
66      endif
67!
68      elset='
69     &'
70      noset='
71     &'
72      irefnode=0
73      irotnode=0
74!
75      do i=2,n
76        if(textpart(i)(1:6).eq.'ELSET=') then
77          if(noset(1:1).eq.' ') then
78            elset(1:80)=textpart(i)(7:86)
79            ipos=index(elset,' ')
80            elset(ipos:ipos)='E'
81          else
82            write(*,*) '*ERROR reading *RIGID BODY: either NSET or'
83            write(*,*) '       ELSET can be specified, not both'
84            ier=1
85            return
86          endif
87        elseif(textpart(i)(1:8).eq.'PINNSET=') then
88          if(elset(1:1).eq.' ') then
89            noset(1:80)=textpart(i)(9:88)
90            ipos=index(noset,' ')
91            noset(ipos:ipos)='N'
92          else
93            write(*,*) '*ERROR reading *RIGID BODY: either NSET or'
94            write(*,*) '       ELSET can be specified, not both'
95            ier=1
96            return
97          endif
98        elseif(textpart(i)(1:5).eq.'NSET=') then
99          if(elset(1:1).eq.' ') then
100            noset(1:80)=textpart(i)(6:85)
101            ipos=index(noset,' ')
102            noset(ipos:ipos)='N'
103          else
104            write(*,*) '*ERROR reading *RIGID BODY: either NSET or'
105            write(*,*) '       ELSET can be specified, not both'
106            ier=1
107            return
108          endif
109        elseif(textpart(i)(1:8).eq.'REFNODE=') then
110          read(textpart(i)(9:18),'(i10)',iostat=istat) irefnode
111          if(istat.gt.0) then
112            call inputerror(inpc,ipoinpc,iline,
113     &           "*RIGID BODY%",ier)
114            return
115          endif
116          if(irefnode.gt.nk) then
117            write(*,*) '*ERROR reading *RIGID BODY: ref node',
118     &           irefnode
119            write(*,*) '       has not been defined'
120            ier=1
121            return
122          endif
123        elseif(textpart(i)(1:8).eq.'ROTNODE=') then
124          read(textpart(i)(9:18),'(i10)',iostat=istat) irotnode
125          if(istat.gt.0) then
126            call inputerror(inpc,ipoinpc,iline,
127     &           "*RIGID BODY%",ier)
128            return
129          endif
130          if(irotnode.gt.nk) then
131            write(*,*) '*ERROR reading *RIGID BODY: rot node',
132     &           irotnode
133            write(*,*) '       has not been defined'
134            ier=1
135            return
136          endif
137        else
138          write(*,*)
139     &         '*WARNING reading *RIGID BODY: parameter not recognized:'
140          write(*,*) '         ',
141     &         textpart(i)(1:index(textpart(i),' ')-1)
142          call inputwarning(inpc,ipoinpc,iline,
143     &         "*RIGID BODY%")
144        endif
145      enddo
146!
147!     check whether a set was defined
148!
149      if((elset(1:1).eq.' ').and.
150     &     (noset(1:1).eq.' ')) then
151        write(*,*) '*WARNING reading *RIGID BODY: no set defined'
152        call getnewline(inpc,textpart,istat,n,key,iline,ipol,inl,
153     &       ipoinp,inp,ipoinpc)
154        return
155      endif
156!
157      inoset=0
158      ielset=0
159!
160!     checking whether the set exists
161!
162      if(noset(1:1).ne.' ') then
163c     do i=1,nset
164c     if(set(i).eq.noset) then
165c     inoset=i
166c     exit
167c     endif
168c     enddo
169        call cident81(set,noset,nset,id)
170        if(id.gt.0) then
171          if(noset.eq.set(id)) then
172            inoset=id
173          endif
174        endif
175        if(inoset.eq.0) then
176          write(*,*) '*WARNING reading *RIGID BODY: node set ',noset
177          write(*,*) '         does not exist'
178          call getnewline(inpc,textpart,istat,n,key,iline,ipol,inl,
179     &         ipoinp,inp,ipoinpc)
180          return
181        endif
182      endif
183!
184      if(elset(1:1).ne.' ') then
185c     do i=1,nset
186c     if(set(i).eq.elset) then
187c     ielset=i
188c     exit
189c     endif
190c     enddo
191        call cident81(set,elset,nset,id)
192        if(id.gt.0) then
193          if(elset.eq.set(id)) then
194            ielset=id
195          endif
196        endif
197        if(ielset.eq.0) then
198          write(*,*) '*WARNING reading *RIGID BODY: element set ',
199     &         elset
200          write(*,*) '         does not exist'
201          call getnewline(inpc,textpart,istat,n,key,iline,ipol,inl,
202     &         ipoinp,inp,ipoinpc)
203          return
204        endif
205      endif
206!
207!     check for the existence of irefnode and irotnode; if none were
208!     defined, new nodes are generated
209!
210      if(irefnode.eq.0) then
211        nk=nk+1
212        if(nk.gt.nk_) then
213          write(*,*) '*ERROR reading *RIGID BODY: increase nk_'
214          ier=1
215          return
216        endif
217        irefnode=nk
218!
219!     default position of the reference node is the origin
220!
221        co(1,nk)=0.d0
222        co(2,nk)=0.d0
223        co(3,nk)=0.d0
224      endif
225!
226      if(irotnode.eq.0) then
227        nk=nk+1
228        if(nk.gt.nk_) then
229          write(*,*) '*ERROR reading *RIGID BODY: increase nk_'
230          ier=1
231          return
232        endif
233        irotnode=nk
234      endif
235!
236!     check whether other equations apply to the dependent nodes
237!
238      if(inoset.ne.0) then
239        do i=istartset(inoset),iendset(inoset)
240          node=ialset(i)
241          if(node.gt.nk_) then
242            write(*,*) '*ERROR reading *RIGID BODY: node ',node
243            write(*,*) '       belonging to set ',noset
244            write(*,*) '       has not been defined'
245            ier=1
246            return
247          endif
248          if((node.eq.irefnode).or.(node.eq.irotnode)) cycle
249          do j=1,3
250            idof=8*(node-1)+j
251            call nident(ikmpc,idof,nmpc,id)
252            if(id.gt.0) then
253              if(ikmpc(id).eq.idof) then
254                write(*,*) '*WARNING reading *RIGID BODY: dof ',j
255                write(*,*) '         of node ',node,' belonging'
256                write(*,*) '         to a rigid body is detected'
257                write(*,*) '         on the dependent side of '
258                write(*,*) '         another equation; no rigid'
259                write(*,*) '         body constrained applied'
260              endif
261            endif
262          enddo
263        enddo
264      endif
265!
266      if(ielset.ne.0) then
267        do i=istartset(ielset),iendset(ielset)
268          ielement=ialset(i)
269          if(ielement.gt.ne_) then
270            write(*,*) '*ERROR reading *RIGID BODY: element ',
271     &           ielement
272            write(*,*) '       belonging to set ',elset
273            write(*,*) '       has not been defined'
274            ier=1
275            return
276          endif
277          if(ipkon(ielement).lt.0) cycle
278          indexe=ipkon(ielement)
279          if(lakon(ielement)(4:4).eq.'2') then
280            nope=20
281          elseif(lakon(ielement)(4:4).eq.'8') then
282            nope=8
283          elseif(lakon(ielement)(4:5).eq.'10') then
284            nope=10
285          elseif(lakon(ielement)(4:4).eq.'4') then
286            nope=4
287          elseif(lakon(ielement)(4:5).eq.'15') then
288            nope=15
289          else
290            nope=6
291          endif
292          do k=indexe+1,indexe+nope
293            node=kon(k)
294            if((node.eq.irefnode).or.(node.eq.irotnode)) cycle
295            do j=1,3
296              idof=8*(node-1)+j
297              call nident(ikmpc,idof,nmpc,id)
298              if(id.gt.0) then
299                if(ikmpc(id).eq.idof) then
300                  write(*,*)'*WARNING reading *RIGID BODY: dof ',
301     &                 j,'of node ',node,' belonging to a'
302                  write(*,*)'         rigid body is detected on th
303     &e dependent side of another'
304                  write(*,*)'         equation; no rigid body cons
305     &trained applied'
306                endif
307              endif
308            enddo
309          enddo
310        enddo
311      endif
312!
313!     generating the equations in basis form
314!
315!     node set
316!
317      if(inoset.ne.0) then
318        do i=istartset(inoset),iendset(inoset)
319          node=ialset(i)
320          if(node.gt.0) then
321            if((node.eq.irefnode).or.(node.eq.irotnode)) cycle
322            call rigidmpc(ipompc,nodempc,coefmpc,irefnode,irotnode,
323     &           labmpc,nmpc,nmpc_,mpcfree,ikmpc,ilmpc,nk,nk_,
324     &           nodeboun,ndirboun,ikboun,ilboun,nboun,nboun_,node,
325     &           typeboun,co,jmin,jmax)
326          else
327            node=ialset(i-2)
328            do
329              node=node-ialset(i)
330              if(node.ge.ialset(i-1)) exit
331              if((node.eq.irefnode).or.(node.eq.irotnode)) cycle
332              call rigidmpc(ipompc,nodempc,coefmpc,irefnode,
333     &             irotnode,labmpc,nmpc,nmpc_,mpcfree,ikmpc,ilmpc,
334     &             nk,nk_,nodeboun,ndirboun,ikboun,ilboun,nboun,
335     &             nboun_,node,typeboun,co,jmin,jmax)
336            enddo
337          endif
338        enddo
339      endif
340!
341!     element set
342!
343      if(ielset.ne.0) then
344        do i=istartset(ielset),iendset(ielset)
345          ielement=ialset(i)
346          if(ielement.gt.0) then
347            if(ipkon(ielement).lt.0) cycle
348            indexe=ipkon(ielement)
349            if(lakon(ielement)(4:4).eq.'2') then
350              nope=20
351            elseif(lakon(ielement)(4:4).eq.'8') then
352              nope=8
353            elseif(lakon(ielement)(4:5).eq.'10') then
354              nope=10
355            elseif(lakon(ielement)(4:4).eq.'4') then
356              nope=4
357            elseif(lakon(ielement)(4:5).eq.'15') then
358              nope=15
359            else
360              nope=6
361            endif
362            do k=indexe+1,indexe+nope
363              node=kon(k)
364              if((node.eq.irefnode).or.(node.eq.irotnode)) cycle
365              call rigidmpc(ipompc,nodempc,coefmpc,irefnode,
366     &             irotnode,labmpc,nmpc,nmpc_,mpcfree,ikmpc,ilmpc,
367     &             nk,nk_,nodeboun,ndirboun,ikboun,ilboun,nboun,
368     &             nboun_,node,typeboun,co,jmin,jmax)
369            enddo
370          else
371            ielement=ialset(i-2)
372            do
373              ielement=ielement-ialset(i)
374              if(ielement.ge.ialset(i-1)) exit
375              if(ipkon(ielement).lt.0) cycle
376              indexe=ipkon(ielement)
377              if(lakon(ielement)(4:4).eq.'2') then
378                nope=20
379              elseif(lakon(ielement)(4:4).eq.'8') then
380                nope=8
381              elseif(lakon(ielement)(4:5).eq.'10') then
382                nope=10
383              elseif(lakon(ielement)(4:4).eq.'4') then
384                nope=4
385              elseif(lakon(ielement)(4:5).eq.'15') then
386                nope=15
387              else
388                nope=6
389              endif
390              do k=indexe+1,indexe+nope
391                node=kon(k)
392                if((node.eq.irefnode).or.(node.eq.irotnode)) cycle
393                call rigidmpc(ipompc,nodempc,coefmpc,irefnode,
394     &               irotnode,labmpc,nmpc,nmpc_,mpcfree,ikmpc,
395     &               ilmpc,nk,nk_,nodeboun,ndirboun,ikboun,ilboun,
396     &               nboun,nboun_,node,typeboun,co,jmin,jmax)
397              enddo
398            enddo
399          endif
400        enddo
401      endif
402!
403      call getnewline(inpc,textpart,istat,n,key,iline,ipol,inl,
404     &     ipoinp,inp,ipoinpc)
405!
406      return
407      end
408
409