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 solidsections(inpc,textpart,set,istartset,iendset,
20     &  ialset,nset,ielmat,matname,nmat,ielorien,orname,norien,
21     &  lakon,thicke,kon,ipkon,irstrt,istep,istat,n,iline,ipol,inl,
22     &  ipoinp,inp,cs,mcs,iaxial,ipoinpc,mi,co,ixfree,xnor,iponor,
23     &  ier,orab)
24!
25!     reading the input deck: *SOLID SECTION
26!
27      implicit none
28!
29      logical nodalthickness
30!
31      character*1 inpc(*)
32      character*8 lakon(*)
33      character*80 matname(*),orname(*),material,orientation
34      character*81 set(*),elset
35      character*132 textpart(16)
36!
37      integer mi(*),istartset(*),iendset(*),ialset(*),ielmat(mi(3),*),
38     &  ielorien(mi(3),*),kon(*),ipkon(*),indexe,irstrt(*),nset,nmat,
39     &  norien,ielem,node1,node2,m,indexx,ixfree,iponor(2,*),id,
40     &  istep,istat,n,key,i,j,k,l,imaterial,iorientation,ipos,
41     &  iline,ipol,inl,ipoinp(2,*),inp(3,*),mcs,iaxial,ipoinpc(0:*),
42     &  ier,numnod
43!
44      real*8 thicke(mi(3),*),thickness,pi,cs(17,*),xn(3),co(3,*),p(3),
45     &     dd,xnor(*),orab(7,*)
46!
47      if((istep.gt.0).and.(irstrt(1).ge.0)) then
48         write(*,*) '*ERROR reading *SOLID SECTION: *SOLID SECTION'
49         write(*,*)'       should be placed before all step definitions'
50         ier=1
51         return
52      endif
53!
54      nodalthickness=.false.
55      pi=4.d0*datan(1.d0)
56!
57      orientation='
58     &                           '
59      elset='
60     &                      '
61      ipos=0
62!
63      do i=2,n
64         if(textpart(i)(1:9).eq.'MATERIAL=') then
65            material=textpart(i)(10:89)
66         elseif(textpart(i)(1:12).eq.'ORIENTATION=') then
67            orientation=textpart(i)(13:92)
68         elseif(textpart(i)(1:6).eq.'ELSET=') then
69            elset=textpart(i)(7:86)
70            elset(81:81)=' '
71            ipos=index(elset,' ')
72            elset(ipos:ipos)='E'
73         elseif(textpart(i)(1:14).eq.'NODALTHICKNESS') then
74            nodalthickness=.true.
75         else
76            write(*,*)
77     &      '*WARNING reading *SOLID SECTION: parameter not recognized:'
78            write(*,*) '         ',
79     &                 textpart(i)(1:index(textpart(i),' ')-1)
80            call inputwarning(inpc,ipoinpc,iline,
81     &"*SOLID SECTION%")
82         endif
83      enddo
84!
85!     check for the existence of the material
86!
87      do i=1,nmat
88         if(matname(i).eq.material) exit
89      enddo
90      if(i.gt.nmat) then
91         do i=1,nmat
92            if(matname(i)(1:11).eq.'ANISO_CREEP') then
93               if(matname(i)(12:20).eq.material(1:9)) exit
94            elseif(matname(i)(1:10).eq.'ANISO_PLAS') then
95               if(matname(i)(11:20).eq.material(1:10)) exit
96            endif
97         enddo
98      endif
99      if(i.gt.nmat) then
100         write(*,*)
101     &      '*ERROR reading *SOLID SECTION: nonexistent material'
102         write(*,*) '  '
103         call inputerror(inpc,ipoinpc,iline,
104     &        "*SOLID SECTION%",ier)
105         return
106      endif
107      imaterial=i
108!
109!     check for the existence of the orientation
110!
111      if(orientation.eq.'                    ') then
112         iorientation=0
113      else
114         do i=1,norien
115            if(orname(i).eq.orientation) exit
116         enddo
117         if(i.gt.norien) then
118            write(*,*)
119     &       '*ERROR reading *SOLID SECTION: nonexistent orientation'
120            write(*,*) '  '
121            call inputerror(inpc,ipoinpc,iline,
122     &           "*SOLID SECTION%",ier)
123            return
124         endif
125         iorientation=i
126      endif
127!
128!     check for the existence of the set
129!
130      if(ipos.eq.0) then
131         write(*,*) '*ERROR reading *SOLID SECTION: no element set ',
132     &        elset
133         write(*,*) '       was been defined. '
134         call inputerror(inpc,ipoinpc,iline,
135     &        "*SOLID SECTION%",ier)
136         return
137      endif
138c      do i=1,nset
139c         if(set(i).eq.elset) exit
140c      enddo
141      call cident81(set,elset,nset,id)
142      i=nset+1
143      if(id.gt.0) then
144        if(elset.eq.set(id)) then
145          i=id
146        endif
147      endif
148      if(i.gt.nset) then
149         elset(ipos:ipos)=' '
150         write(*,*) '*ERROR reading *SOLID SECTION: element set ',elset
151         write(*,*) '  has not yet been defined. '
152         call inputerror(inpc,ipoinpc,iline,
153     &        "*SOLID SECTION%",ier)
154         return
155      endif
156!
157!     assigning the elements of the set the appropriate material
158!     and orientation number
159!
160      do j=istartset(i),iendset(i)
161         if(ialset(j).gt.0) then
162            k=ialset(j)
163            if((lakon(k)(1:1).eq.'B').or.
164     &         (lakon(k)(1:1).eq.'S')) then
165               write(*,*)
166     &          '*ERROR reading *SOLID SECTION: *SOLID SECTION can'
167               write(*,*) '       not be used for beam or shell elements
168     &'
169               write(*,*) '       Faulty element: ',k
170               ier=1
171               return
172            endif
173            ielmat(1,k)=imaterial
174            if(ielorien(1,k).lt.0) then
175!
176!              an orientation based on a distribution has the
177!              same name as the distribution. However, to a
178!              distribution which is not used in any orientation
179!              definition no local coordinate system has been
180!              assigned (i.e. orab(7,..)=0.d0)
181!
182               if((orname(-ielorien(1,k)).eq.orientation).and.
183     &              (orab(7,-ielorien(1,k)).ne.0.d0)) then
184                  ielorien(1,k)=-ielorien(1,k)
185               else
186                  ielorien(1,k)=iorientation
187               endif
188            else
189               ielorien(1,k)=iorientation
190            endif
191         else
192            k=ialset(j-2)
193            do
194               k=k-ialset(j)
195               if(k.ge.ialset(j-1)) exit
196               if((lakon(k)(1:1).eq.'B').or.
197     &              (lakon(k)(1:1).eq.'S')) then
198                  write(*,*) '*ERROR reading *SOLID SECTION: *SOLID SECT
199     &ION can'
200                  write(*,*) '       not be used for beam or shell eleme
201     &nts'
202                  write(*,*) '       Faulty element: ',k
203                  ier=1
204                  return
205               endif
206               ielmat(1,k)=imaterial
207               if(ielorien(1,k).lt.0) then
208!
209!                 an orientation based on a distribution has the
210!                 same name as the distribution. However, to a
211!                 distribution which is not used in any orientation
212!                 definition no local coordinate system has been
213!                 assigned (i.e. orab(7,..)=0.d0)
214!
215                  if((orname(-ielorien(1,k)).eq.orientation).and.
216     &                 (orab(7,-ielorien(1,k)).ne.0.d0)) then
217                     ielorien(1,k)=-ielorien(1,k)
218                  else
219                     ielorien(1,k)=iorientation
220                  endif
221               else
222                  ielorien(1,k)=iorientation
223               endif
224            enddo
225         endif
226      enddo
227!
228      call getnewline(inpc,textpart,istat,n,key,iline,ipol,inl,
229     &     ipoinp,inp,ipoinpc)
230      if(istat.lt.0) return
231!
232!     assigning a thickness to plane stress/strain elements and an angle to
233!     axisymmetric elements
234!
235      if(key.eq.0) then
236!
237!        second line with thickness given
238!
239         read(textpart(1)(1:20),'(f20.0)',iostat=istat) thickness
240         if(istat.gt.0) then
241            call inputerror(inpc,ipoinpc,iline,
242     &           "*SOLID SECTION%",ier)
243            return
244         endif
245!
246!        for axial symmetric structures:
247!           thickness for axial symmetric elements: 2 degrees
248!           thickness for plane stress elements: reduced by 180
249!           thickness for plane strain elements: reduced by 180
250!
251         if(.not.nodalthickness) then
252            if(iaxial.eq.180) then
253               if(lakon(ialset(istartset(i)))(1:2).eq.'CA') then
254                  thickness=datan(1.d0)*8.d0/iaxial
255               elseif(lakon(ialset(istartset(i)))(1:3).eq.'CPS') then
256                  thickness=thickness/iaxial
257               elseif(lakon(ialset(istartset(i)))(1:3).eq.'CPE') then
258                  thickness=thickness/iaxial
259               endif
260            endif
261         else
262!
263!           for those elements for which nodal thickness is activated
264!           the thickness is set to -1.d0
265!
266            thickness=-1.d0
267         endif
268!
269!        assigning the thickness to each node of the corresponding
270!        elements (thickness specified)
271!
272         do j=istartset(i),iendset(i)
273            if(ialset(j).gt.0) then
274               if((lakon(ialset(j))(1:2).eq.'CP').or.
275     &              (lakon(ialset(j))(1:2).eq.'CA')) then
276!
277!                 plane stress/strain or axisymmetric elements
278!
279                  indexe=ipkon(ialset(j))
280                  read(lakon(ialset(j))(4:4),'(i1)') numnod
281                  do l=1,numnod
282                     thicke(1,indexe+l)=thickness
283                  enddo
284               elseif(lakon(ialset(j))(1:1).eq.'T') then
285                  ielem=ialset(j)
286!
287!                 default cross section for trusses is the
288!                 rectangular cross section
289!
290                  lakon(ielem)(8:8)='R'
291                  indexe=ipkon(ielem)
292                  node1=kon(indexe+1)
293                  if(lakon(ielem)(4:4).eq.'2') then
294                     node2=kon(indexe+2)
295                  else
296                     node2=kon(indexe+3)
297                  endif
298!
299!                 determining a vector orthogonal to the truss
300!                 element
301!
302                  do l=1,3
303                     xn(l)=co(l,node2)-co(l,node1)
304                  enddo
305                  if(dabs(xn(1)).gt.0.d0) then
306                     p(1)=-xn(3)
307                     p(2)=0.d0
308                     p(3)=xn(1)
309                  elseif(dabs(xn(2)).gt.0.d0) then
310                     p(1)=xn(2)
311                     p(2)=-xn(1)
312                     p(3)=0.d0
313                  else
314                     p(1)=0.d0
315                     p(2)=xn(3)
316                     p(3)=-xn(2)
317                  endif
318                  dd=dsqrt(p(1)*p(1)+p(2)*p(2)+p(3)*p(3))
319                  if(dd.lt.1.d-10) then
320                     write(*,*)
321     &                    '*ERROR reading *SOLID SECTION: normal'
322                     write(*,*) '       in direction 1 has zero size'
323                     ier=1
324                     return
325                  endif
326                  do l=1,3
327                     p(l)=p(l)/dd
328                  enddo
329                  do l=1,3
330                     thicke(1,indexe+l)=dsqrt(thickness)
331                     thicke(2,indexe+l)=dsqrt(thickness)
332                     indexx=ixfree
333                     do m=1,3
334                        xnor(indexx+m)=p(m)
335                     enddo
336                     ixfree=ixfree+6
337                     iponor(1,indexe+l)=indexx
338                  enddo
339               endif
340            else
341               k=ialset(j-2)
342               do
343                  k=k-ialset(j)
344                  if(k.ge.ialset(j-1)) exit
345                  if((lakon(k)(1:2).eq.'CP').or.
346     &                 (lakon(k)(1:2).eq.'CA')) then
347                     indexe=ipkon(k)
348                     read(lakon(k)(4:4),'(i1)') numnod
349                     do l=1,numnod
350                        thicke(1,indexe+l)=thickness
351                     enddo
352                  elseif(lakon(k)(1:1).eq.'T') then
353!
354!                    default cross section for trusses is the
355!                    rectangular cross section
356!
357                     lakon(k)(8:8)='R'
358                     indexe=ipkon(k)
359                     node1=kon(indexe+1)
360                     if(lakon(k)(4:4).eq.'2') then
361                        node2=kon(indexe+2)
362                     else
363                        node2=kon(indexe+3)
364                     endif
365!
366!                    determining a vector orthogonal to the truss
367!                    element
368!
369                     do l=1,3
370                        xn(l)=co(l,node2)-co(l,node1)
371                     enddo
372                     if(dabs(xn(1)).gt.0.d0) then
373                        p(1)=-xn(3)
374                        p(2)=0.d0
375                        p(3)=xn(1)
376                     elseif(dabs(xn(2)).gt.0.d0) then
377                        p(1)=xn(2)
378                        p(2)=-xn(1)
379                        p(3)=0.d0
380                     else
381                        p(1)=0.d0
382                        p(2)=xn(3)
383                        p(3)=-xn(2)
384                     endif
385                     dd=dsqrt(p(1)*p(1)+p(2)*p(2)+p(3)*p(3))
386                     if(dd.lt.1.d-10) then
387                        write(*,*)
388     &                       '*ERROR reading *SOLID SECTION: normal'
389                        write(*,*) '       in direction 1 has zero size'
390                        ier=1
391                        return
392                     endif
393                     do l=1,3
394                        p(l)=p(l)/dd
395                     enddo
396                     do l=1,3
397                        thicke(1,indexe+l)=dsqrt(thickness)
398                        thicke(2,indexe+l)=dsqrt(thickness)
399                        indexx=ixfree
400                        do m=1,3
401                           xnor(indexx+m)=p(m)
402                        enddo
403                        ixfree=ixfree+6
404                        iponor(1,indexe+l)=indexx
405                     enddo
406                  endif
407               enddo
408            endif
409         enddo
410      else
411!
412!        no second line (no thickness given)
413!
414!        assigning the thickness to each node of the corresponding
415!        elements (thickness not specified: only axisymmetric elements
416!        or plane stress/strain elements with nodal thickness)
417!
418         thickness=datan(1.d0)*8.d0/iaxial
419         do j=istartset(i),iendset(i)
420            if(ialset(j).gt.0) then
421               if(lakon(ialset(j))(1:2).eq.'CA') then
422!
423!                 axisymmetric elements
424!
425                  indexe=ipkon(ialset(j))
426                  read(lakon(ialset(j))(4:4),'(i1)') numnod
427                  do l=1,numnod
428                     thicke(1,indexe+l)=thickness
429                  enddo
430               elseif((lakon(ialset(j))(1:3).eq.'CPS').or.
431     &                 (lakon(ialset(j))(1:3).eq.'CPE')) then
432                  if(nodalthickness) then
433                     indexe=ipkon(ialset(j))
434                     read(lakon(ialset(j))(4:4),'(i1)') numnod
435                     do l=1,numnod
436                        thicke(1,indexe+l)=-1.d0
437                     enddo
438                  endif
439               endif
440            else
441               k=ialset(j-2)
442               do
443                  k=k-ialset(j)
444                  if(k.ge.ialset(j-1)) exit
445                  if(lakon(k)(1:2).eq.'CA') then
446!
447!                 axisymmetric elements
448!
449                     indexe=ipkon(k)
450                     read(lakon(k)(4:4),'(i1)') numnod
451                     do l=1,numnod
452                        thicke(1,indexe+l)=thickness
453                     enddo
454                  elseif((lakon(k)(1:3).eq.'CPS').or.
455     &                    (lakon(k)(1:3).eq.'CPE')) then
456                     if(nodalthickness) then
457                        indexe=ipkon(k)
458                        read(lakon(k)(4:4),'(i1)') numnod
459                        do l=1,numnod
460                           thicke(1,indexe+l)=-1.d0
461                        enddo
462                     endif
463                  endif
464               enddo
465            endif
466         enddo
467      endif
468!
469!        defining cyclic symmetric conditions for axisymmetric
470!        elements (needed for cavity radiation)
471!
472      do j=istartset(i),iendset(i)
473         if(ialset(j).gt.0) then
474            if(lakon(ialset(j))(1:2).eq.'CA') then
475               if(mcs.gt.1) then
476                  write(*,*) '*ERROR reading *SOLID SECTION: '
477                  write(*,*) '       axisymmetric elements cannot be
478     &combined with cyclic symmetry'
479                  ier=1
480                  return
481               elseif(mcs.eq.1) then
482                  if(int(cs(1,1)).ne.int(2.d0*pi/thickness+0.5d0))
483     &                 then
484                     write(*,*) '*ERROR reading *SOLID SECTION: '
485                     write(*,*) '       it is not allowed to define t
486     &wo different'
487                     write(*,*) '       angles for an axisymmetric st
488     &ructure'
489                     ier=1
490                     return
491                  else
492                     exit
493                  endif
494               endif
495               mcs=1
496               cs(1,1)=2.d0*pi/thickness+0.5d0
497               cs(2,1)=-0.5d0
498               cs(3,1)=-0.5d0
499               cs(5,1)=1.5d0
500               do k=6,9
501                  cs(k,1)=0.d0
502               enddo
503               cs(10,1)=1.d0
504               cs(11,1)=0.d0
505               cs(12,1)=-1.d0
506               cs(14,1)=0.5d0
507               cs(15,1)=dcos(thickness)
508               cs(16,1)=dsin(thickness)
509               exit
510            endif
511         else
512            k=ialset(j-2)
513            do
514               k=k-ialset(j)
515               if(k.ge.ialset(j-1)) exit
516               if(lakon(k)(1:2).eq.'CA') then
517                  if(mcs.gt.1) then
518                     write(*,*) '*ERROR reading *SOLID SECTION: '
519                     write(*,*) '       axisymmetric elements cannot
520     &be combined with cyclic symmetry'
521                     ier=1
522                     return
523                  elseif(mcs.eq.1) then
524                     if(int(cs(1,1)).ne.int(2.d0*pi/thickness+0.5d0))
525     &                    then
526                        write(*,*) '*ERROR reading *SOLID SECTION: '
527                        write(*,*) '       it is not allowed to defin
528     &e two different'
529                        write(*,*) '       angles for an axisymmetric
530     &structure'
531                        ier=1
532                        return
533                     else
534                        exit
535                     endif
536                  endif
537                  mcs=1
538                  cs(1,1)=2.d0*pi/thickness+0.5d0
539                  cs(2,1)=-0.5d0
540                  cs(3,1)=-0.5d0
541                  cs(5,1)=1.5d0
542                  do k=6,9
543                     cs(k,1)=0.d0
544                  enddo
545                  cs(10,1)=1.d0
546                  cs(11,1)=0.d0
547                  cs(12,1)=-1.d0
548                  cs(14,1)=0.5d0
549                  cs(15,1)=dcos(thickness)
550                  cs(16,1)=dsin(thickness)
551                  exit
552               endif
553            enddo
554         endif
555      enddo
556!
557      if(key.eq.0) then
558         call getnewline(inpc,textpart,istat,n,key,iline,ipol,inl,
559     &        ipoinp,inp,ipoinpc)
560      endif
561!
562      return
563      end
564
565