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 velinireltoabs(ibody,xbody,cbody,nbody,set,
20     &     istartset,iendset,ialset,nset,veold,mi,ipkon,kon,lakon,co,
21     &     itreated)
22!
23!     assigns the body forces to the elements by use of field ipobody
24!
25      implicit none
26!
27      character*8 lakon(*)
28      character*81 cbody(*),elset,set(*)
29!
30      integer ibody(3,*),i,j,k,l,m,istartset(*),nbody,kon(*),id,
31     &     iendset(*),ialset(*),nset,istat,indexe,mi(*),ipkon(*),
32     &     node,nope,itreated(*)
33!
34      real*8 xbody(7,*),veold(0:mi(2),*),co(3,*),om,a(3),xn(3),p(3),dd
35!
36      do k=1,nbody
37         elset=cbody(k)
38!
39!        check for centrifugal loads
40!
41         if(ibody(1,k).eq.1) then
42            ibody(1,k)=0
43         else
44            cycle
45         endif
46!
47!        determine
48!          om: rotational speed
49!          a: point on axis
50!          xn: unit vector on axis
51!
52         om=dsqrt(xbody(1,k))
53         a(1)=xbody(2,k)
54         a(2)=xbody(3,k)
55         a(3)=xbody(4,k)
56         xn(1)=xbody(5,k)
57         xn(2)=xbody(6,k)
58         xn(3)=xbody(7,k)
59!
60!     check whether element number or set name
61!
62         read(elset,'(i21)',iostat=istat) l
63         if(istat.eq.0) then
64!
65!           treat element l
66!
67            indexe=ipkon(l)
68!
69!           nope is the number of nodes belonging to the element
70!
71            if(lakon(l)(4:5).eq.'20') then
72               nope=20
73            elseif(lakon(l)(4:4).eq.'8') then
74               nope=8
75            elseif(lakon(l)(4:5).eq.'10') then
76               nope=10
77            elseif(lakon(l)(4:4).eq.'4') then
78               nope=4
79            elseif(lakon(l)(4:5).eq.'15') then
80               nope=15
81            elseif(lakon(l)(4:4).eq.'6') then
82               nope=6
83            elseif(lakon(l)(1:2).eq.'ES') then
84!
85!              contact elements need not be treated
86!
87               if(lakon(l)(7:7).ne.'C') then
88                  nope=ichar(lakon(l)(8:8))-47
89               endif
90            elseif(lakon(l)(1:4).eq.'MASS') then
91               nope=1
92            endif
93!
94            do i=1,nope
95               node=kon(indexe+i)
96               if(itreated(node).eq.1) cycle
97               do m=1,3
98                  p(m)=co(m,node)-a(m)
99               enddo
100               dd=p(1)*xn(1)+p(2)*xn(2)+p(3)*xn(3)
101               do m=1,3
102                  p(m)=(p(m)-dd*xn(m))*om
103               enddo
104               veold(1,node)=veold(1,node)+
105     &                       xn(2)*p(3)-xn(3)*p(2)
106               veold(2,node)=veold(2,node)+
107     &                       xn(3)*p(1)-xn(1)*p(3)
108               veold(3,node)=veold(3,node)+
109     &              xn(1)*p(2)-xn(2)*p(1)
110               itreated(node)=1
111            enddo
112         else
113!
114!     set name
115!
116c            do i=1,nset
117c               if(set(i).eq.elset) exit
118c            enddo
119           call cident81(set,elset,nset,id)
120           i=nset+1
121           if(id.gt.0) then
122             if(elset.eq.set(id)) then
123               i=id
124             endif
125           endif
126!
127            do j=istartset(i),iendset(i)
128               if(ialset(j).gt.0) then
129                  l=ialset(j)
130!
131!                 treat element l
132!
133                  indexe=ipkon(l)
134!
135!                 nope is the number of nodes belonging to the element
136!
137                  if(lakon(l)(4:5).eq.'20') then
138                     nope=20
139                  elseif(lakon(l)(4:4).eq.'8') then
140                     nope=8
141                  elseif(lakon(l)(4:5).eq.'10') then
142                     nope=10
143                  elseif(lakon(l)(4:4).eq.'4') then
144                     nope=4
145                  elseif(lakon(l)(4:5).eq.'15') then
146                     nope=15
147                  elseif(lakon(l)(4:4).eq.'6') then
148                     nope=6
149                  elseif(lakon(l)(1:2).eq.'ES') then
150!
151!                    contact elements need not be treated
152!
153                     if(lakon(l)(7:7).ne.'C') then
154                        nope=ichar(lakon(l)(8:8))-47
155                     endif
156                  elseif(lakon(l)(1:4).eq.'MASS') then
157                     nope=1
158                  endif
159!
160                  do i=1,nope
161                     node=kon(indexe+i)
162                     if(itreated(node).eq.1) cycle
163                     do m=1,3
164                        p(m)=co(m,node)-a(m)
165                     enddo
166                     dd=p(1)*xn(1)+p(2)*xn(2)+p(3)*xn(3)
167                     do m=1,3
168                        p(m)=(p(m)-dd*xn(m))*om
169                     enddo
170                     veold(1,node)=veold(1,node)+
171     &                    xn(2)*p(3)-xn(3)*p(2)
172                     veold(2,node)=veold(2,node)+
173     &                    xn(3)*p(1)-xn(1)*p(3)
174                     veold(3,node)=veold(3,node)+
175     &                    xn(1)*p(2)-xn(2)*p(1)
176                     itreated(node)=1
177                  enddo
178               else
179                  l=ialset(j-2)
180                  do
181                     l=l-ialset(j)
182                     if(l.ge.ialset(j-1)) exit
183!
184!                    treat element l
185!
186                     indexe=ipkon(l)
187!
188!                    nope is the number of nodes belonging to the element
189!
190                     if(lakon(l)(4:5).eq.'20') then
191                        nope=20
192                     elseif(lakon(l)(4:4).eq.'8') then
193                        nope=8
194                     elseif(lakon(l)(4:5).eq.'10') then
195                        nope=10
196                     elseif(lakon(l)(4:4).eq.'4') then
197                        nope=4
198                     elseif(lakon(l)(4:5).eq.'15') then
199                        nope=15
200                     elseif(lakon(l)(4:4).eq.'6') then
201                        nope=6
202                     elseif(lakon(l)(1:2).eq.'ES') then
203!
204!                       contact elements need not be treated
205!
206                        if(lakon(l)(7:7).ne.'C') then
207                           nope=ichar(lakon(l)(8:8))-47
208                        endif
209                     elseif(lakon(l)(1:4).eq.'MASS') then
210                        nope=1
211                     endif
212!
213                     do i=1,nope
214                        node=kon(indexe+i)
215                        if(itreated(node).eq.1) cycle
216                        do m=1,3
217                           p(m)=co(m,node)-a(m)
218                        enddo
219                        dd=p(1)*xn(1)+p(2)*xn(2)+p(3)*xn(3)
220                        do m=1,3
221                           p(m)=(p(m)-dd*xn(m))*om
222                        enddo
223                        veold(1,node)=veold(1,node)+
224     &                       xn(2)*p(3)-xn(3)*p(2)
225                        veold(2,node)=veold(2,node)+
226     &                       xn(3)*p(1)-xn(1)*p(3)
227                        veold(3,node)=veold(3,node)+
228     &                       xn(1)*p(2)-xn(2)*p(1)
229                        itreated(node)=1
230                     enddo
231                  enddo
232               endif
233            enddo
234         endif
235      enddo
236!
237      return
238      end
239
240