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 dloads(inpc,textpart,set,istartset,iendset,
20     &     ialset,nset,nelemload,sideload,xload,nload,nload_,
21     &     ielmat,iamload,amname,nam,lakon,ne,dload_flag,istep,
22     &     istat,n,iline,ipol,inl,ipoinp,inp,cbody,ibody,xbody,nbody,
23     &     nbody_,xbodyold,iperturb,physcon,nam_,namtot_,namta,amta,
24     &     nmethod,ipoinpc,maxsectors,mi,idefload,idefbody,ipkon,
25     &     thicke,iamplitudedefault,namtot,ier)
26!
27!     reading the input deck: *DLOAD
28!
29      implicit none
30!
31      logical dload_flag,submodel,edgeload,surface
32!
33      character*1 inpc(*)
34      character*8 lakon(*)
35      character*20 sideload(*),label
36      character*80 amname(*),amplitude
37      character*81 set(*),elset,cbody(*)
38      character*132 textpart(16)
39!
40      integer istartset(*),iendset(*),ialset(*),nelemload(2,*),mi(*),
41     &     ielmat(mi(3),*),nset,nload,nload_,istep,istat,n,i,j,l,key,
42     &     iamload(2,*),nam,iamplitude,ipos,ne,iline,ipol,iperturb(*),
43     &     inl,ipoinp(2,*),inp(3,*),ibody(3,*),nbody,nbody_,nam_,namtot,
44     &     namtot_,namta(3,*),idelay,nmethod,lc,isector,node,id,
45     &     ipoinpc(0:*),maxsectors,jsector,iglobstep,idefload(*),
46     %     idefbody(*),ipkon(*),k,indexe,iamplitudedefault,ier
47!
48      real*8 xload(2,*),xbody(7,*),xmagnitude,dd,p1(3),p2(3),bodyf(3),
49     &     xbodyold(7,*),physcon(*),amta(2,*),xxmagnitude,
50     &     thicke(mi(3),*),thickness
51!
52      iamplitude=iamplitudedefault
53      idelay=0
54      lc=1
55      isector=0
56      submodel=.false.
57      iglobstep=0
58      edgeload=.false.
59      surface=.false.
60!
61      if(istep.lt.1) then
62        write(*,*) '*ERROR reading *DLOAD: *DLOAD should only be used'
63        write(*,*) '  within a STEP'
64        ier=1
65        return
66      endif
67!
68      do i=2,n
69        if((textpart(i)(1:6).eq.'OP=NEW').and.(.not.dload_flag)) then
70          do j=1,nload
71            if(sideload(j)(1:1).eq.'P') then
72              sideload(j)(3:4)='  '
73              xload(1,j)=0.d0
74            endif
75          enddo
76          do j=1,nbody
77            xbody(1,j)=0.d0
78          enddo
79        elseif(textpart(i)(1:10).eq.'AMPLITUDE=') then
80          read(textpart(i)(11:90),'(a80)') amplitude
81          do j=1,nam
82            if(amname(j).eq.amplitude) then
83              iamplitude=j
84              exit
85            endif
86          enddo
87          if(j.gt.nam) then
88            write(*,*)'*ERROR reading *DLOAD: nonexistent amplitude'
89            write(*,*) '  '
90            call inputerror(inpc,ipoinpc,iline,
91     &           "*DLOAD%",ier)
92            return
93          endif
94          iamplitude=j
95        elseif(textpart(i)(1:10).eq.'TIMEDELAY=') THEN
96          if(idelay.ne.0) then
97            write(*,*)
98     &           '*ERROR reading *DLOAD: the parameter TIME DELAY'
99            write(*,*) '       is used twice in the same keyword'
100            write(*,*) '       '
101            call inputerror(inpc,ipoinpc,iline,
102     &           "*DLOAD%",ier)
103            return
104          else
105            idelay=1
106          endif
107          nam=nam+1
108          if(nam.gt.nam_) then
109            write(*,*) '*ERROR reading *DLOAD: increase nam_'
110            ier=1
111            return
112          endif
113          amname(nam)='
114     &'
115          if(iamplitude.eq.0) then
116            write(*,*) '*ERROR reading *DLOAD: time delay must be'
117            write(*,*) '       preceded by the amplitude parameter'
118            ier=1
119            return
120          endif
121          namta(3,nam)=sign(iamplitude,namta(3,iamplitude))
122          iamplitude=nam
123c     if(nam.eq.1) then
124c     namtot=0
125c     else
126c     namtot=namta(2,nam-1)
127c     endif
128          namtot=namtot+1
129          if(namtot.gt.namtot_) then
130            write(*,*) '*ERROR dloads: increase namtot_'
131            ier=1
132            return
133          endif
134          namta(1,nam)=namtot
135          namta(2,nam)=namtot
136c     call reorderampl(amname,namta,nam)
137          read(textpart(i)(11:30),'(f20.0)',iostat=istat)
138     &         amta(1,namtot)
139          if(istat.gt.0) then
140            call inputerror(inpc,ipoinpc,iline,
141     &           "*DLOAD%",ier)
142            return
143          endif
144        elseif(textpart(i)(1:9).eq.'LOADCASE=') then
145          read(textpart(i)(10:19),'(i10)',iostat=istat) lc
146          if(istat.gt.0) then
147            call inputerror(inpc,ipoinpc,iline,
148     &           "*DLOAD%",ier)
149            return
150          endif
151          if(nmethod.ne.5) then
152            write(*,*)
153     &           '*ERROR reading *DLOAD: the parameter LOAD CASE'
154            write(*,*) '       is only allowed in STEADY STATE'
155            write(*,*) '       DYNAMICS calculations'
156            ier=1
157            return
158          endif
159        elseif(textpart(i)(1:7).eq.'SECTOR=') then
160          read(textpart(i)(8:17),'(i10)',iostat=istat) isector
161          if(istat.gt.0) then
162            call inputerror(inpc,ipoinpc,iline,
163     &           "*DLOAD%",ier)
164            return
165          endif
166          if((nmethod.le.3).or.(iperturb(1).gt.1)) then
167            write(*,*) '*ERROR reading *DLOAD: the parameter SECTOR'
168            write(*,*) '       is only allowed in MODAL DYNAMICS or'
169            write(*,*) '       STEADY STATE DYNAMICS calculations'
170            ier=1
171            return
172          endif
173          if(isector.gt.maxsectors) then
174            write(*,*) '*ERROR reading *DLOAD: sector ',isector
175            write(*,*) '       exceeds number of sectors'
176            ier=1
177            return
178          endif
179          isector=isector-1
180        elseif(textpart(i)(1:8).eq.'SUBMODEL') then
181          submodel=.true.
182        elseif(textpart(i)(1:5).eq.'STEP=') then
183          read(textpart(i)(6:15),'(i10)',iostat=istat) iglobstep
184          if(istat.gt.0) then
185            call inputerror(inpc,ipoinpc,iline,
186     &           "*DLOAD%",ier)
187            return
188          endif
189        elseif(textpart(i)(1:8).eq.'DATASET=') then
190          read(textpart(i)(9:18),'(i10)',iostat=istat) iglobstep
191          if(istat.gt.0) then
192            call inputerror(inpc,ipoinpc,iline,
193     &           "*DLOAD%",ier)
194            return
195          endif
196!
197!     the mode number for submodels
198!     is stored as a negative global step
199!
200          iglobstep=-iglobstep
201        else
202          write(*,*)
203     &         '*WARNING reading *DLOAD: parameter not recognized:'
204          write(*,*) '         ',
205     &         textpart(i)(1:index(textpart(i),' ')-1)
206          call inputwarning(inpc,ipoinpc,iline,
207     &         "*DLOAD%")
208        endif
209      enddo
210!
211!     check for the presence of an amplitude in submodel cases
212!
213      if(submodel) then
214        if(iamplitude.ne.0) then
215          write(*,*) '*WARNING reading *DSLOAD:'
216          write(*,*) '         no amplitude definition is allowed'
217          write(*,*) '         in combination with a submodel'
218        endif
219      endif
220!
221!     check whether global step was specified for submodel
222!
223      if((submodel).and.(iglobstep.eq.0)) then
224        write(*,*) '*ERROR reading *DLOAD: no global step'
225        write(*,*) '       step specified for the submodel'
226        call inputerror(inpc,ipoinpc,iline,
227     &       "*DLOAD%",ier)
228        return
229      endif
230!
231      do
232        call getnewline(inpc,textpart,istat,n,key,iline,ipol,inl,
233     &       ipoinp,inp,ipoinpc)
234        if((istat.lt.0).or.(key.eq.1)) return
235!
236        read(textpart(2)(1:20),'(a20)',iostat=istat) label
237!
238!     for submodels the load label is modified and the global
239!     step is stored in iamload(1,*)
240!
241        if(submodel) then
242          label(3:4)='SM'
243          iamplitude=iglobstep
244        endif
245!
246        if(label(3:4).ne.'NP') then
247          read(textpart(3)(1:20),'(f20.0)',iostat=istat) xmagnitude
248        else
249          read(textpart(3)(1:10),'(i10)',iostat=istat) node
250        endif
251        if(istat.gt.0) then
252          call inputerror(inpc,ipoinpc,iline,
253     &         "*DLOAD%",ier)
254          return
255        endif
256        if(label(1:7).eq.'CENTRIF') then
257          do i=1,3
258            read(textpart(i+3)(1:20),'(f20.0)',iostat=istat) p1(i)
259            if(istat.gt.0) then
260              call inputerror(inpc,ipoinpc,iline,
261     &             "*DLOAD%",ier)
262              return
263            endif
264          enddo
265          do i=1,3
266            read(textpart(i+6)(1:20),'(f20.0)',iostat=istat) p2(i)
267            if(istat.gt.0) then
268              call inputerror(inpc,ipoinpc,iline,
269     &             "*DLOAD%",ier)
270              return
271            endif
272          enddo
273          dd=dsqrt(p2(1)**2+p2(2)**2+p2(3)**2)
274          do i=1,3
275            p2(i)=p2(i)/dd
276          enddo
277        elseif(label(1:4).eq.'GRAV') then
278          do i=1,3
279            read(textpart(i+3)(1:20),'(f20.0)',iostat=istat) bodyf(i)
280            if(istat.gt.0) then
281              call inputerror(inpc,ipoinpc,iline,
282     &             "*DLOAD%",ier)
283              return
284            endif
285          enddo
286        elseif(label(1:6).eq.'NEWTON') then
287          if(iperturb(1).le.1) then
288            write(*,*) '*ERROR reading *DLOAD: NEWTON gravity force'
289            write(*,*) '       can only be used in a nonlinear'
290            write(*,*) '       procedure'
291            ier=1
292            return
293          endif
294          if(physcon(3).le.0.d0) then
295            write(*,*) '*ERROR reading *DLOAD: NEWTON gravity force'
296            write(*,*) '       requires the definition of a'
297            write(*,*) '       positive gravity constant with'
298            write(*,*) '       a *PHYSICAL CONSTANTS card'
299            ier=1
300            return
301          endif
302        elseif(((label(1:2).ne.'P1').and.(label(1:2).ne.'P2').and.
303     &         (label(1:2).ne.'P3').and.(label(1:2).ne.'P4').and.
304     &         (label(1:2).ne.'P5').and.(label(1:2).ne.'P6').and.
305     &         (label(1:2).ne.'P ').and.(label(1:2).ne.'BX').and.
306     &         (label(1:2).ne.'BY').and.(label(1:2).ne.'BZ').and.
307c     BernhardiStart
308     &         (label(1:2).ne.'ED')).or.
309     &         ((label(3:6).ne.'NOR1').and.(label(3:6).ne.'NOR2').and.
310     &         (label(3:6).ne.'NOR3').and.(label(3:6).ne.'NOR4')).and.
311c     BernhardiEnd
312     &         ((label(3:4).ne.'  ').and.(label(3:4).ne.'NU').and.
313     &         (label(3:4).ne.'NP').and.(label(3:4).ne.'SM'))) then
314          call inputerror(inpc,ipoinpc,iline,
315     &         "*DLOAD%",ier)
316          return
317        endif
318!
319        read(textpart(1)(1:10),'(i10)',iostat=istat) l
320        if(istat.eq.0) then
321          if(l.gt.ne) then
322            write(*,*) '*ERROR reading *DLOAD: element ',l
323            write(*,*) '       is not defined'
324            ier=1
325            return
326          endif
327          if((label(1:7).eq.'CENTRIF').or.(label(1:4).eq.'GRAV').or.
328     &         (label(1:6).eq.'NEWTON')) then
329            elset(1:80)=textpart(1)(1:80)
330            elset(81:81)=' '
331            call bodyadd(cbody,ibody,xbody,nbody,nbody_,elset,label,
332     &           iamplitude,xmagnitude,p1,p2,bodyf,xbodyold,lc,idefbody)
333          else
334            xxmagnitude=xmagnitude
335            if((lakon(l)(1:2).eq.'CP').or.
336     &           (lakon(l)(2:2).eq.'A').or.
337     &           (lakon(l)(7:7).eq.'E').or.
338     &           (lakon(l)(7:7).eq.'S').or.
339     &           (lakon(l)(7:7).eq.'A')) then
340              if(label(1:2).eq.'P1') then
341                label(1:2)='P3'
342              elseif(label(1:2).eq.'P2') then
343                label(1:2)='P4'
344              elseif(label(1:2).eq.'P3') then
345                label(1:2)='P5'
346              elseif(label(1:2).eq.'P4') then
347                label(1:2)='P6'
348              endif
349            elseif((lakon(l)(1:1).eq.'B').or.
350     &             (lakon(l)(7:7).eq.'B')) then
351              if(label(1:2).eq.'P2') label(1:2)='P5'
352            elseif((lakon(l)(1:1).eq.'S').or.
353     &             (lakon(l)(7:7).eq.'L')) then
354c     BernhardiStart
355              if(label(1:6).eq.'EDNOR1') then
356                label(1:2)='P3'
357                edgeload=.true.
358              elseif(label(1:6).eq.'EDNOR2') then
359                label(1:2)='P4'
360                edgeload=.true.
361              elseif(label(1:6).eq.'EDNOR3') then
362                label(1:2)='P5'
363                edgeload=.true.
364              elseif(label(1:6).eq.'EDNOR4') then
365                label(1:2)='P6'
366                edgeload=.true.
367              else
368                label(1:2)='P1'
369              endif
370!
371!     EDNOR is an edge load
372!
373              if(edgeload) then
374                indexe=ipkon(l)
375                thickness=0.d0
376                do k=1,mi(3)
377                  if(ielmat(k,l).ne.0) then
378                    thickness=thickness+thicke(k,indexe+1)
379                  else
380                    exit
381                  endif
382                enddo
383                xxmagnitude=xmagnitude/thickness
384              endif
385c     BernhardiEnd
386            endif
387            if(lc.ne.1) then
388              jsector=isector+maxsectors
389            else
390              jsector=isector
391            endif
392            if(label(3:4).ne.'NP') then
393              call loadadd(l,label,xxmagnitude,nelemload,sideload,
394     &             xload,nload,nload_,iamload,iamplitude,
395     &             nam,jsector,idefload)
396            else
397              call loadaddp(l,label,nelemload,sideload,
398     &             xload,nload,nload_,iamload,iamplitude,
399     &             nam,node)
400            endif
401          endif
402        else
403          read(textpart(1)(1:80),'(a80)',iostat=istat) elset
404          elset(81:81)=' '
405          ipos=index(elset,' ')
406!
407!     check for element set
408!
409          elset(ipos:ipos)='E'
410c          do i=1,nset
411c            if(set(i).eq.elset) exit
412c          enddo
413          call cident81(set,elset,nset,id)
414          i=nset+1
415          if(id.gt.0) then
416            if(elset.eq.set(id)) then
417              i=id
418            endif
419          endif
420          if(i.gt.nset) then
421!
422!     check for facial surface
423!
424            surface=.true.
425            elset(ipos:ipos)='T'
426c            do i=1,nset
427c              if(set(i).eq.elset) exit
428c            enddo
429            call cident81(set,elset,nset,id)
430            i=nset+1
431            if(id.gt.0) then
432              if(elset.eq.set(id)) then
433                i=id
434              endif
435            endif
436            if(i.gt.nset) then
437              elset(ipos:ipos)=' '
438              write(*,*) '*ERROR reading *DLOAD: element set '
439              write(*,*) '       or facial surface ',elset
440              write(*,*) '       has not yet been defined. '
441              call inputerror(inpc,ipoinpc,iline,
442     &             "*DLOAD%",ier)
443              return
444            endif
445          endif
446!
447          if((label(1:7).eq.'CENTRIF').or.(label(1:4).eq.'GRAV').or.
448     &         (label(1:6).eq.'NEWTON')) then
449            call bodyadd(cbody,ibody,xbody,nbody,nbody_,elset,label,
450     &           iamplitude,xmagnitude,p1,p2,bodyf,xbodyold,lc,idefbody)
451          else
452            l=ialset(istartset(i))
453            if(surface) then
454              write(label(2:2),'(i1)') l-10*(l/10)
455              l=l/10
456            endif
457            if((lakon(l)(1:2).eq.'CP').or.
458     &           (lakon(l)(2:2).eq.'A').or.
459     &           (lakon(l)(7:7).eq.'E').or.
460     &           (lakon(l)(7:7).eq.'S').or.
461     &           (lakon(l)(7:7).eq.'A')) then
462              if(label(1:2).eq.'P1') then
463                label(1:2)='P3'
464              elseif(label(1:2).eq.'P2') then
465                label(1:2)='P4'
466              elseif(label(1:2).eq.'P3') then
467                label(1:2)='P5'
468              elseif(label(1:2).eq.'P4') then
469                label(1:2)='P6'
470              endif
471            elseif((lakon(l)(1:1).eq.'B').or.
472     &             (lakon(l)(7:7).eq.'B')) then
473              if(label(1:2).eq.'P2') label(1:2)='P5'
474            elseif((lakon(l)(1:1).eq.'S').or.
475     &             (lakon(l)(7:7).eq.'L')) then
476c     BernhardiStart
477              if(label(1:6).eq.'EDNOR1') then
478                label(1:2)='P3'
479                edgeload=.true.
480              elseif(label(1:6).eq.'EDNOR2') then
481                label(1:2)='P4'
482                edgeload=.true.
483              elseif(label(1:6).eq.'EDNOR3') then
484                label(1:2)='P5'
485                edgeload=.true.
486              elseif(label(1:6).eq.'EDNOR4') then
487                label(1:2)='P6'
488                edgeload=.true.
489              else
490                label(1:2)='P1'
491              endif
492c     BernhardiEnd
493            endif
494!
495            do j=istartset(i),iendset(i)
496              if(ialset(j).gt.0) then
497                l=ialset(j)
498                if(surface) then
499                  write(label(2:2),'(i1)') l-10*(l/10)
500                  l=l/10
501                endif
502                xxmagnitude=xmagnitude
503!
504!     EDNOR is an edge load
505!
506                if(edgeload) then
507                  indexe=ipkon(l)
508                  thickness=0.d0
509                  do k=1,mi(3)
510                    if(ielmat(k,l).ne.0) then
511                      thickness=thickness+thicke(k,indexe+1)
512                    else
513                      exit
514                    endif
515                  enddo
516                  xxmagnitude=xmagnitude/thickness
517                endif
518!
519                if(lc.ne.1) then
520                  jsector=isector+maxsectors
521                else
522                  jsector=isector
523                endif
524                if(label(3:4).ne.'NP') then
525                  call loadadd(l,label,xxmagnitude,nelemload,
526     &                 sideload,xload,nload,nload_,iamload,
527     &                 iamplitude,nam,jsector,idefload)
528                else
529                  call loadaddp(l,label,nelemload,
530     &                 sideload,xload,nload,nload_,iamload,
531     &                 iamplitude,nam,node)
532                endif
533              else
534                l=ialset(j-2)
535                do
536                  l=l-ialset(j)
537                  if(l.ge.ialset(j-1)) exit
538                  xxmagnitude=xmagnitude
539!
540!     EDNOR is an edge load
541!
542                  if(edgeload) then
543                    indexe=ipkon(l)
544                    thickness=0.d0
545                    do k=1,mi(3)
546                      if(ielmat(k,l).ne.0) then
547                        thickness=thickness+thicke(k,indexe+1)
548                      else
549                        exit
550                      endif
551                    enddo
552                    xxmagnitude=xmagnitude/thickness
553                  endif
554!
555                  if(lc.ne.1) then
556                    jsector=isector+maxsectors
557                  else
558                    jsector=isector
559                  endif
560                  if(label(3:4).ne.'NP') then
561                    call loadadd(l,label,xxmagnitude,nelemload,
562     &                   sideload,xload,nload,nload_,
563     &                   iamload,iamplitude,nam,jsector,idefload)
564                  else
565                    call loadaddp(l,label,nelemload,
566     &                   sideload,xload,nload,nload_,
567     &                   iamload,iamplitude,nam,node)
568                  endif
569                enddo
570              endif
571            enddo
572          endif
573        endif
574      enddo
575!
576      return
577      end
578