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