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 beamgeneralsections(inpc,textpart,set,istartset,
20     &     iendset,ialset,nset,ielmat,matname,nmat,ielorien,orname,
21     &     norien,thicke,ipkon,iponor,xnor,ixfree,offset,lakon,irstrt,
22     &     istep,istat,n,iline,ipol,inl,ipoinp,inp,ipoinpc,mi,ielprop,
23     &     nprop,nprop_,prop,nelcon,ier)
24!
25!     reading the input deck: *BEAM SECTION
26!     for sections of type PIPE, BOX and GENERAL
27!
28!     this routine is used for sections which cannot be described
29!     by two thicknesses alone -> prop array must be used
30!
31      implicit none
32!
33      character*1 inpc(*)
34      character*4 section
35      character*8 lakon(*)
36      character*80 matname(*),orname(*),material,orientation
37      character*81 set(*),elset
38      character*132 textpart(16)
39!
40      integer istartset(*),iendset(*),ialset(*),mi(*),ielmat(mi(3),*),
41     &     ipoinpc(0:*),ielorien(mi(3),*),ipkon(*),iline,ipol,inl,lprop,
42     &     ipoinp(2,*),inp(3,*),nset,nmat,norien,istep,istat,n,key,i,j,
43     &     imaterial,iorientation,ipos,m,iponor(2,*),ixfree,indexx,
44     &     irstrt(*),ielprop(*),nprop_,nprop,npropstart,ndprop,
45     &     nelcon(2,*),ier,id,k,l,indexe,ndpropread
46!
47      real*8 thicke(mi(3),*),thickness1,thickness2,p(3),xnor(*),
48     &     offset(2,*),offset1,offset2,dd,prop(*)
49!
50      if((istep.gt.0).and.(irstrt(1).ge.0)) then
51        write(*,*) '*ERROR reading *BEAM SECTION:'
52        write(*,*) '       *BEAM SECTION should'
53        write(*,*) '       be placed before all step definitions'
54        ier=1
55        return
56      endif
57!
58      offset1=0.d0
59      offset2=0.d0
60      orientation='
61     &'
62      section='    '
63      ipos=1
64      npropstart=nprop
65!
66      do i=2,n
67        if(textpart(i)(1:9).eq.'MATERIAL=') then
68          material=textpart(i)(10:89)
69        elseif(textpart(i)(1:12).eq.'ORIENTATION=') then
70          orientation=textpart(i)(13:92)
71        elseif(textpart(i)(1:6).eq.'ELSET=') then
72          elset=textpart(i)(7:86)
73          elset(21:21)=' '
74          ipos=index(elset,' ')
75          elset(ipos:ipos)='E'
76        elseif(textpart(i)(1:8).eq.'SECTION=') then
77          if(textpart(i)(9:12).eq.'PIPE') then
78            section='PIPE'
79            ndprop=2
80          elseif(textpart(i)(9:11).eq.'BOX') then
81            section='BOX'
82            ndprop=6
83          elseif(textpart(i)(9:15).eq.'GENERAL') then
84            section='GENE'
85            ndprop=5
86          else
87            write(*,*)
88     &           '*ERROR reading *BEAM SECTION: unknown section'
89            ier=1
90            return
91          endif
92        elseif(textpart(i)(1:8).eq.'OFFSET1=') then
93          read(textpart(i)(9:28),'(f20.0)',iostat=istat) offset1
94          if(istat.gt.0) then
95            call inputerror(inpc,ipoinpc,iline,
96     &           "*BEAM SECTION%",ier)
97            return
98          endif
99        elseif(textpart(i)(1:8).eq.'OFFSET2=') then
100          read(textpart(i)(9:28),'(f20.0)',iostat=istat) offset2
101          if(istat.gt.0) then
102            call inputerror(inpc,ipoinpc,iline,
103     &           "*BEAM SECTION%",ier)
104            return
105          endif
106        else
107          write(*,*) '*WARNING reading *BEAM SECTION:'
108          write(*,*) '         parameter not recognized:'
109          write(*,*) '         ',
110     &         textpart(i)(1:index(textpart(i),' ')-1)
111          call inputwarning(inpc,ipoinpc,iline,
112     &         "*BEAM SECTION%")
113        endif
114      enddo
115!
116!     check whether a sections was defined
117!
118      if(section.eq.'    ') then
119        write(*,*) '*ERROR reading *BEAM SECTION:'
120        write(*,*) '       no section defined'
121        ier=1
122        return
123      endif
124!
125!     check for the existence of the set,the material and orientation
126!
127      do i=1,nmat
128        if(matname(i).eq.material) exit
129      enddo
130      if(i.gt.nmat) then
131        write(*,*) '*ERROR reading *BEAM SECTION:'
132        write(*,*) '       nonexistent material'
133        write(*,*) '  '
134        call inputerror(inpc,ipoinpc,iline,
135     &       "*BEAM SECTION%",ier)
136        return
137      endif
138      imaterial=i
139!
140      if(orientation.eq.'
141     &') then
142        iorientation=0
143c     elseif(nelcon(1,i).eq.2) then
144c     write(*,*) '*INFO reading *BEAM SECTION: an orientation'
145c     write(*,*) '      is for isotropic materials irrelevant'
146c     call inputinfo(inpc,ipoinpc,iline,
147c     &"*BEAM SECTION%")
148c     iorientation=0
149      else
150        do i=1,norien
151          if(orname(i).eq.orientation) exit
152        enddo
153        if(i.gt.norien) then
154          write(*,*)
155     &         '*ERROR reading *BEAM SECTION: nonexistent orientation'
156          write(*,*) '  '
157          call inputerror(inpc,ipoinpc,iline,
158     &         "*BEAM SECTION%",ier)
159          return
160        endif
161        iorientation=i
162      endif
163!
164c     do i=1,nset
165c     if(set(i).eq.elset) exit
166c     enddo
167      call cident81(set,elset,nset,id)
168      i=nset+1
169      if(id.gt.0) then
170        if(elset.eq.set(id)) then
171          i=id
172        endif
173      endif
174      if(i.gt.nset) then
175        elset(ipos:ipos)=' '
176        write(*,*)'*ERROR reading *BEAM SECTION: element set ',
177     &       elset(1:ipos)
178        write(*,*) '  has not yet been defined. '
179        call inputerror(inpc,ipoinpc,iline,
180     &       "*BEAM SECTION%",ier)
181        return
182      endif
183!
184!     assigning the elements of the set the appropriate material,
185!     orientation number, section and offset(s)
186!
187      if(section.ne.'GENE') then
188        do j=istartset(i),iendset(i)
189          if(ialset(j).gt.0) then
190            if(lakon(ialset(j))(1:4).ne.'B32R') then
191              write(*,*) '*ERROR reading *BEAM SECTION:'
192              write(*,*) '       *BEAM SECTION of type ',
193     &             section,' can'
194              write(*,*) '       only be used for B32R elements.'
195              write(*,*) '       Element ',ialset(j),' is not a B32R
196     &element.'
197              ier=1
198              return
199            endif
200            ielmat(1,ialset(j))=imaterial
201            ielorien(1,ialset(j))=iorientation
202            offset(1,ialset(j))=offset1
203            offset(2,ialset(j))=offset2
204            if(section.eq.'PIPE') then
205              lakon(ialset(j))(8:8)='P'
206            elseif(section.eq.'BOX') then
207              lakon(ialset(j))(8:8)='B'
208            endif
209          else
210            k=ialset(j-2)
211            do
212              k=k-ialset(j)
213              if(k.ge.ialset(j-1)) exit
214              if(lakon(k)(1:1).ne.'B') then
215                write(*,*) '*ERROR reading *BEAM SECTION:'
216                write(*,*) '       *BEAM SECTION of type ',
217     &               section,' can'
218                write(*,*) '       only be used for beam elements.'
219                write(*,*) '       Element ',k,' is not a beam elem
220     &ent.'
221                ier=1
222                return
223              endif
224              ielmat(1,k)=imaterial
225              ielorien(1,k)=iorientation
226              offset(1,k)=offset1
227              offset(2,k)=offset2
228              if(section.eq.'PIPE') then
229                lakon(k)(8:8)='P'
230              elseif(section.eq.'BOX') then
231                lakon(k)(8:8)='B'
232              endif
233            enddo
234          endif
235        enddo
236      else
237!
238!     general section
239!
240        do j=istartset(i),iendset(i)
241          if(ialset(j).gt.0) then
242            if(lakon(ialset(j))(1:5).ne.'U1   ') then
243              write(*,*) '*ERROR reading *BEAM SECTION:'
244              write(*,*) '       *BEAM SECTION of type GENERAL can'
245              write(*,*) '       only be used for U1 elements.'
246              write(*,*) '       Element ',ialset(j),' is not a U1
247     &element.'
248              ier=1
249              return
250            endif
251            ielmat(1,ialset(j))=imaterial
252            ielorien(1,ialset(j))=iorientation
253          else
254            k=ialset(j-2)
255            do
256              k=k-ialset(j)
257              if(k.ge.ialset(j-1)) exit
258              if(lakon(k)(1:5).ne.'U1   ') then
259                write(*,*) '*ERROR reading *BEAM SECTION:'
260                write(*,*) '       *BEAM SECTION of type GENERAL'
261                write(*,*) '       can only be used for beam'
262                write(*,*) '       elements.'
263                write(*,*) '       Element ',k,' is not a beam elem
264     &ent.'
265                ier=1
266                return
267              endif
268              ielmat(1,k)=imaterial
269              ielorien(1,k)=iorientation
270            enddo
271          endif
272        enddo
273      endif
274!
275!     reading the properties
276!
277      lprop=0
278      ndpropread=ndprop
279      do j=1,(ndpropread-1)/8+1
280        call getnewline(inpc,textpart,istat,n,key,iline,ipol,inl,
281     &       ipoinp,inp,ipoinpc)
282        do k=1,8
283          lprop=lprop+1
284          if(lprop.gt.ndpropread) exit
285          read(textpart(k),'(f40.0)',iostat=istat)
286     &         prop(nprop+lprop)
287          if(istat.gt.0) then
288            call inputerror(inpc,ipoinpc,iline,
289     &           "*BEAM SECTION%",ier)
290            return
291          endif
292        enddo
293      enddo
294      nprop=nprop+ndprop
295!
296      if(section.eq.'GENE') then
297!
298        call getnewline(inpc,textpart,istat,n,key,iline,ipol,inl,
299     &       ipoinp,inp,ipoinpc)
300        if((istat.lt.0).or.(key.eq.1)) then
301!
302!     default 1-direction
303!
304          prop(nprop+1)=0.d0
305          prop(nprop+2)=0.d0
306          prop(nprop+3)=-1.d0
307        else
308!
309!     1-direction specified by the user
310!
311          read(textpart(1)(1:20),'(f20.0)',iostat=istat) p(1)
312          if(istat.gt.0) then
313            call inputerror(inpc,ipoinpc,iline,
314     &           "*BEAM SECTION%",ier)
315            return
316          endif
317          read(textpart(2)(1:20),'(f20.0)',iostat=istat) p(2)
318          if(istat.gt.0) then
319            call inputerror(inpc,ipoinpc,iline,
320     &           "*BEAM SECTION%",ier)
321            return
322          endif
323          read(textpart(3)(1:20),'(f20.0)',iostat=istat) p(3)
324          if(istat.gt.0) then
325            call inputerror(inpc,ipoinpc,iline,
326     &           "*BEAM SECTION%",ier)
327            return
328          endif
329          dd=dsqrt(p(1)*p(1)+p(2)*p(2)+p(3)*p(3))
330          if(dd.lt.1.d-10) then
331            write(*,*)
332     &           '*ERROR reading *BEAM SECTION: normal in direction 1'
333            write(*,*) '       has zero size'
334            ier=1
335            return
336          endif
337          do j=1,3
338            prop(nprop+j)=p(j)/dd
339          enddo
340        endif
341        nprop=nprop+3
342!
343        prop(nprop+1)=offset1
344        prop(nprop+2)=offset2
345        nprop=nprop+2
346      endif
347!
348      if(nprop.gt.nprop_) then
349        write(*,*)
350     &       '*ERROR reading *BEAM SECTION: increase nprop_'
351        ier=1
352        return
353      endif
354!
355!     calculating the dimensions of the rectangular parent beam
356!
357      if(section.eq.'PIPE') then
358        thickness1=2.d0*prop(npropstart+1)
359        thickness2=thickness1
360      elseif(section.eq.'BOX') then
361        thickness1=prop(npropstart+1)
362        thickness2=prop(npropstart+2)
363      endif
364!
365!     assigning the thickness and the properties to the elements
366!
367      if(section.ne.'GENE') then
368        do j=istartset(i),iendset(i)
369          if(ialset(j).gt.0) then
370            indexe=ipkon(ialset(j))
371            do l=1,8
372              thicke(1,indexe+l)=thickness1
373              thicke(2,indexe+l)=thickness2
374            enddo
375            ielprop(ialset(j))=npropstart
376          else
377            k=ialset(j-2)
378            do
379              k=k-ialset(j)
380              if(k.ge.ialset(j-1)) exit
381              indexe=ipkon(k)
382              do l=1,8
383                thicke(1,indexe+l)=thickness1
384                thicke(2,indexe+l)=thickness2
385              enddo
386              ielprop(k)=npropstart
387            enddo
388          endif
389        enddo
390      else
391        do j=istartset(i),iendset(i)
392          if(ialset(j).gt.0) then
393            ielprop(ialset(j))=npropstart
394          else
395            k=ialset(j-2)
396            do
397              k=k-ialset(j)
398              if(k.ge.ialset(j-1)) exit
399              ielprop(k)=npropstart
400            enddo
401          endif
402        enddo
403      endif
404!
405      call getnewline(inpc,textpart,istat,n,key,iline,ipol,inl,
406     &     ipoinp,inp,ipoinpc)
407      if((istat.lt.0).or.(key.eq.1)) return
408!
409!     assigning normal direction 1 for the beam
410!
411      indexx=-1
412      read(textpart(1)(1:20),'(f20.0)',iostat=istat) p(1)
413      if(istat.gt.0) then
414        call inputerror(inpc,ipoinpc,iline,
415     &       "*BEAM SECTION%",ier)
416        return
417      endif
418      read(textpart(2)(1:20),'(f20.0)',iostat=istat) p(2)
419      if(istat.gt.0) then
420        call inputerror(inpc,ipoinpc,iline,
421     &       "*BEAM SECTION%",ier)
422        return
423      endif
424      read(textpart(3)(1:20),'(f20.0)',iostat=istat) p(3)
425      if(istat.gt.0) then
426        call inputerror(inpc,ipoinpc,iline,
427     &       "*BEAM SECTION%",ier)
428        return
429      endif
430      dd=dsqrt(p(1)*p(1)+p(2)*p(2)+p(3)*p(3))
431      if(dd.lt.1.d-10) then
432        write(*,*)
433     &       '*ERROR reading *BEAM SECTION: normal in direction 1'
434        write(*,*) '       has zero size'
435        ier=1
436        return
437      endif
438      do j=1,3
439        p(j)=p(j)/dd
440      enddo
441      do j=istartset(i),iendset(i)
442        if(ialset(j).gt.0) then
443          indexe=ipkon(ialset(j))
444          do l=1,8
445            if(indexx.eq.-1) then
446              indexx=ixfree
447              do m=1,3
448                xnor(indexx+m)=p(m)
449              enddo
450              ixfree=ixfree+6
451            endif
452            iponor(1,indexe+l)=indexx
453          enddo
454        else
455          k=ialset(j-2)
456          do
457            k=k-ialset(j)
458            if(k.ge.ialset(j-1)) exit
459            indexe=ipkon(k)
460            do l=1,8
461              if(indexx.eq.-1) then
462                indexx=ixfree
463                do m=1,3
464                  xnor(indexx+m)=p(m)
465                enddo
466                ixfree=ixfree+6
467              endif
468              iponor(1,indexe+l)=indexx
469            enddo
470          enddo
471        endif
472      enddo
473!
474      call getnewline(inpc,textpart,istat,n,key,iline,ipol,inl,
475     &     ipoinp,inp,ipoinpc)
476!
477      return
478      end
479