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 orientations(inpc,textpart,orname,orab,norien,
20     &  norien_,istep,istat,n,iline,ipol,inl,ipoinp,inp,ipoinpc,ier)
21!
22!     reading the input deck: *ORIENTATION
23!
24      implicit none
25!
26      logical distribution
27!
28      character*1 inpc(*)
29      character*80 orname(*),distname,oriename
30      character*132 textpart(16)
31!
32      integer norien,norien_,istep,istat,n,key,i,iline,ipol,inl,
33     &  ipoinp(2,*),inp(3,*),ipoinpc(0:*),iaxis,j,ier,idis,iorien
34!
35      real*8 orab(7,*),a(3,3),c(3,3),angle,p(3),dc,ds,pi,system,ab(6)
36!
37      if(istep.gt.0) then
38         write(*,*)
39     &       '*ERROR reading *ORIENTATION: *ORIENTATION should be'
40         write(*,*) '  placed before all step definitions'
41         ier=1
42         return
43      endif
44!
45      distribution=.false.
46!
47!     rectangular coordinate system: orab(7,norien)=1
48!     cylindrical coordinate system: orab(7,norien)=-1
49!     default is rectangular
50!
51      system=1.d0
52      iaxis=0
53!
54      do i=2,n
55         if(textpart(i)(1:5).eq.'NAME=') then
56            oriename=textpart(i)(6:85)
57            if(textpart(i)(86:86).ne.' ') then
58               write(*,*) '*ERROR reading *ORIENTATION: name too long'
59               write(*,*) '       (more than 80 characters)'
60               write(*,*) '       orientation name:',textpart(i)(1:132)
61               ier=1
62               return
63            endif
64         elseif(textpart(i)(1:7).eq.'SYSTEM=') then
65            if(textpart(i)(8:8).eq.'C') then
66               system=-1.d0
67            endif
68         else
69            write(*,*)
70     &        '*WARNING reading *ORIENTATION: parameter not recognized:'
71            write(*,*) '         ',
72     &                 textpart(i)(1:index(textpart(i),' ')-1)
73            call inputwarning(inpc,ipoinpc,iline,
74     &"*ORIENTATION%")
75         endif
76      enddo
77!
78      call getnewline(inpc,textpart,istat,n,key,iline,ipol,inl,
79     &     ipoinp,inp,ipoinpc)
80      if((istat.lt.0).or.(key.eq.1)) then
81         write(*,*)
82     &      '*ERROR reading *ORIENTATION: definition of the following'
83         write(*,*) '  orientation is not complete: ',orname(norien)
84         ier=1
85         return
86      endif
87!
88      do i=1,6
89         read(textpart(i)(1:20),'(f20.0)',iostat=istat) ab(i)
90         if(istat.gt.0) then
91            distribution=.true.
92            read(textpart(1)(1:80),'(a80)',iostat=istat) distname
93            exit
94         endif
95      enddo
96!
97      call getnewline(inpc,textpart,istat,n,key,iline,ipol,inl,
98     &     ipoinp,inp,ipoinpc)
99!
100      if((istat.ge.0).and.(key.ne.1)) then
101         read(textpart(1)(1:10),'(i10)',iostat=istat) iaxis
102         if(istat.gt.0) then
103            write(*,*) '*ERROR reading *ORIENTATION: expected',
104     &        'axis of rotation'
105            call inputerror(inpc,ipoinpc,iline,
106     &           "*ORIENTATION%",ier)
107            return
108         endif
109         read(textpart(2)(1:20),'(f20.0)',iostat=istat) angle
110         if(istat.gt.0) then
111            write(*,*) '*ERROR reading *ORIENTATION: expected',
112     &        'angle of rotation'
113            call inputerror(inpc,ipoinpc,iline,
114     &           "*ORIENTATION%",ier)
115            return
116         endif
117         call getnewline(inpc,textpart,istat,n,key,iline,ipol,inl,
118     &        ipoinp,inp,ipoinpc)
119      endif
120!
121      if(.not.distribution) then
122!
123         norien=norien+1
124         if(norien.gt.norien_) then
125            write(*,*) '*ERROR reading *ORIENTATION: increase norien_'
126            ier=1
127            return
128         endif
129         do i=1,6
130            orab(i,norien)=ab(i)
131         enddo
132         orab(7,norien)=system
133         orname(norien)=oriename
134         if(iaxis.eq.0) return
135!
136!        additional rotation about an angle only for rectangular
137!        coordinate systems
138!
139         if(orab(7,norien).lt.0.d0) then
140            write(*,*) '*ERROR reading *ORIENTATION'
141            write(*,*) '       additional rotation about an angle'
142            write(*,*) '       is only allowed for rectangular systems'
143            ier=1
144            return
145         endif
146!
147         call transformatrix(orab(1,norien),p,a)
148!
149!        vector on the rotation axis
150!
151         do i=1,3
152            p(i)=a(i,iaxis)
153         enddo
154!
155!        rotation matrix
156!
157         pi=4.d0*datan(1.d0)
158         angle=angle*pi/180.d0
159         dc=dcos(angle)
160         ds=dsin(angle)
161         c(1,1)=dc+(1.d0-dc)*p(1)*p(1)
162         c(1,2)=-ds*p(3)+(1.d0-dc)*p(1)*p(2)
163         c(1,3)=ds*p(2)+(1.d0-dc)*p(1)*p(3)
164         c(2,1)=ds*p(3)+(1.d0-dc)*p(2)*p(1)
165         c(2,2)=dc+(1.d0-dc)*p(2)*p(2)
166         c(2,3)=-ds*p(1)+(1.d0-dc)*p(2)*p(3)
167         c(3,1)=-ds*p(2)+(1.d0-dc)*p(3)*p(1)
168         c(3,2)=ds*p(1)+(1.d0-dc)*p(3)*p(2)
169         c(3,3)=dc+(1.d0-dc)*p(3)*p(3)
170!
171!        rotate vector along the local x-axis and store
172!        as first point in orab
173!
174         do i=1,3
175            orab(i,norien)=0.d0
176            do j=1,3
177               orab(i,norien)=orab(i,norien)+c(i,j)*a(j,1)
178            enddo
179         enddo
180!
181!        rotate vector along the local y-axis and store as
182!        second point in orab
183!
184         do i=1,3
185            orab(i+3,norien)=0.d0
186            do j=1,3
187               orab(i+3,norien)=orab(i+3,norien)+c(i,j)*a(j,2)
188            enddo
189         enddo
190      else
191!
192!        distribution
193!
194         idis=0
195         do iorien=1,norien
196            if(orname(iorien).eq.distname) then
197               idis=idis+1
198               orab(7,iorien)=system
199               orname(iorien)=oriename
200               if(iaxis.ne.0) then
201!
202!     additional rotation about an angle only for rectangular
203!     coordinate systems
204!
205                  if(orab(7,iorien).lt.0.d0) then
206                     write(*,*) '*ERROR reading *ORIENTATION'
207                     write(*,*)
208     &                '       additional rotation about an angle'
209                     write(*,*)
210     &                '       is only allowed for rectangular systems'
211                     ier=1
212                     return
213                  endif
214!
215                  call transformatrix(orab(1,iorien),p,a)
216!
217!     vector on the rotation axis
218!
219                  do i=1,3
220                     p(i)=a(i,iaxis)
221                  enddo
222!
223!     rotation matrix
224!
225                  pi=4.d0*datan(1.d0)
226                  angle=angle*pi/180.d0
227                  dc=dcos(angle)
228                  ds=dsin(angle)
229                  c(1,1)=dc+(1.d0-dc)*p(1)*p(1)
230                  c(1,2)=-ds*p(3)+(1.d0-dc)*p(1)*p(2)
231                  c(1,3)=ds*p(2)+(1.d0-dc)*p(1)*p(3)
232                  c(2,1)=ds*p(3)+(1.d0-dc)*p(2)*p(1)
233                  c(2,2)=dc+(1.d0-dc)*p(2)*p(2)
234                  c(2,3)=-ds*p(1)+(1.d0-dc)*p(2)*p(3)
235                  c(3,1)=-ds*p(2)+(1.d0-dc)*p(3)*p(1)
236                  c(3,2)=ds*p(1)+(1.d0-dc)*p(3)*p(2)
237                  c(3,3)=dc+(1.d0-dc)*p(3)*p(3)
238!
239!     rotate vector along the local x-axis and store
240!     as first point in orab
241!
242                  do i=1,3
243                     orab(i,iorien)=0.d0
244                     do j=1,3
245                        orab(i,iorien)=orab(i,iorien)+c(i,j)*a(j,1)
246                     enddo
247                  enddo
248!
249!     rotate vector along the local y-axis and store as
250!     second point in orab
251!
252                  do i=1,3
253                     orab(i+3,iorien)=0.d0
254                     do j=1,3
255                        orab(i+3,iorien)=orab(i+3,iorien)+c(i,j)*a(j,2)
256                     enddo
257                  enddo
258!
259               endif
260            endif
261         enddo
262!
263!        if no corresponding distribution is found: error
264!
265         if(idis.eq.0) then
266            write(*,*) '*ERROR reading *ORIENTATION: distribution',
267     &          distname,' has not been defined'
268            call inputerror(inpc,ipoinpc,iline,
269     &           "*ORIENTATION%",ier)
270            return
271         endif
272       endif
273c       write(*,*) (orab(i,norien),i=1,7)
274!
275      return
276      end
277
278
279