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 equations(inpc,textpart,ipompc,nodempc,coefmpc,
20     &     nmpc,nmpc_,mpcfree,nk,co,trab,inotr,ntrans,ikmpc,ilmpc,
21     &     labmpc,istep,istat,n,iline,ipol,inl,ipoinp,inp,ipoinpc,
22     &     set,istartset,iendset,ialset,nset,nodempcref,coefmpcref,
23     &     ikmpcref,memmpcref_,mpcfreeref,maxlenmpcref,memmpc_,
24     &     maxlenmpc,ier)
25!
26!     reading the input deck: *EQUATION
27!
28      implicit none
29!
30      character*1 inpc(*)
31      character*20 labmpc(*)
32      character*81 set(*),noset
33      character*132 textpart(16)
34!
35      integer ipompc(*),nodempc(3,*),nmpc,nmpc_,mpcfree,istep,istat,
36     &     n,i,j,ii,key,nterm,number,nk,inotr(2,*),ntrans,node,ndir,
37     &     mpcfreeold,ikmpc(*),ilmpc(*),id,idof,itr,iline,ipol,inl,
38     &     ipoinp(2,*),inp(3,*),ipoinpc(0:*),impcstart,impcend,i1,
39     &     istartset(*),iendset(*),ialset(*),nset,k,l,m,index1,ipos,
40     &     impc,nodempcref(3,*),ikmpcref(*),memmpcref_,mpcfreeref,
41     &     maxlenmpcref,memmpc_,maxlenmpc,ier
42!
43      real*8 coefmpc(*),co(3,*),trab(7,*),a(3,3),x,coefmpcref(*)
44!
45      do m=2,n
46        if(textpart(m)(1:9).eq.'REMOVEALL') then
47!
48          if(istep.eq.1) then
49            write(*,*) '*ERROR reading *EQUATION'
50            write(*,*) '       removing equations is not allowed'
51            write(*,*) '       in the first step'
52            ier=1
53            return
54          endif
55!
56          do j=1,nmpc
57            index1=ipompc(j)
58            if(index1.eq.0) cycle
59            do
60              if(nodempc(3,index1).eq.0) then
61                nodempc(3,index1)=mpcfree
62                mpcfree=ipompc(j)
63                exit
64              endif
65              index1=nodempc(3,index1)
66            enddo
67            ipompc(j)=0
68            ikmpc(j)=0
69            ilmpc(j)=0
70          enddo
71          nmpc=0
72          call getnewline(inpc,textpart,istat,n,key,iline,ipol,inl,
73     &         ipoinp,inp,ipoinpc)
74          return
75        elseif(textpart(m)(1:6).eq.'REMOVE') then
76!
77          if(istep.eq.1) then
78            write(*,*) '*ERROR reading *EQUATION'
79            write(*,*) '       removing equations is not allowed'
80            write(*,*) '       in the first step'
81            ier=1
82            return
83          endif
84!
85          do i=1,nmpc
86            if(ikmpcref(i).ne.ikmpc(i)) then
87              write(*,*) '*ERROR reading *EQUATION'
88              write(*,*) '       The dependent terms in some'
89              write(*,*) '       of the nonlinear equations have'
90              write(*,*) '       changed since the start of the'
91              write(*,*) '       calculation. Removing equations'
92              write(*,*) '       does not work'
93              ier=1
94              return
95            endif
96          enddo
97!
98!     restoring the original equations (before the first call to
99!     cascade)
100!
101          memmpc_=memmpcref_
102          mpcfree=mpcfreeref
103          maxlenmpc=maxlenmpcref
104!
105!     mpcfreeref=-1 indicates that the MPC's have changed and that
106!     a new copy is to be taken before the next call to cascade.c
107!
108          mpcfreeref=-1
109!
110          do i=1,memmpc_
111            do j=1,3
112              nodempc(j,i)=nodempcref(j,i)
113            enddo
114            coefmpc(i)=coefmpcref(i)
115          enddo
116!
117          do
118            call getnewline(inpc,textpart,istat,n,key,iline,ipol,inl,
119     &           ipoinp,inp,ipoinpc)
120            if((istat.lt.0).or.(key.eq.1)) return
121!
122            read(textpart(2)(1:10),'(i10)',iostat=istat) impcstart
123            if(istat.gt.0) then
124              call inputerror(inpc,ipoinpc,iline,
125     &             "*EQUATION%",ier)
126              return
127            endif
128!
129            if(textpart(3)(1:1).eq.' ') then
130              impcend=impcstart
131            else
132              read(textpart(3)(1:10),'(i10)',iostat=istat) impcend
133              if(istat.gt.0) then
134                call inputerror(inpc,ipoinpc,iline,
135     &               "*EQUATION%",ier)
136                return
137              endif
138            endif
139!
140            read(textpart(1)(1:10),'(i10)',iostat=istat) l
141            if(istat.eq.0) then
142              if((l.gt.nk).or.(l.le.0)) then
143                write(*,*) '*ERROR reading *BOUNDARY:'
144                write(*,*) '       node ',l,' is not defined'
145                ier=1
146                return
147              endif
148              do i1=impcstart,impcend
149                idof=8*(l-1)+i1
150                call nident(ikmpc,idof,nmpc,id)
151                if(id.gt.0) then
152                  if(ikmpc(id).eq.idof) then
153                    impc=ilmpc(id)
154                    call mpcrem(impc,mpcfree,nodempc,nmpc,
155     &                   ikmpc,ilmpc,labmpc,coefmpc,ipompc)
156                    cycle
157                  endif
158                endif
159                write(*,*)
160     &               '*WARNING reading *EQUATION: MPC to remove'
161                write(*,*) '         is not defined; node:',l
162                write(*,*) '         degree of freedom:',i1
163              enddo
164            else
165              read(textpart(1)(1:80),'(a80)',iostat=istat) noset
166              noset(81:81)=' '
167              ipos=index(noset,' ')
168              noset(ipos:ipos)='N'
169c              do i=1,nset
170c                if(set(i).eq.noset) exit
171c              enddo
172              call cident81(set,noset,nset,id)
173              i=nset+1
174              if(id.gt.0) then
175                if(noset.eq.set(id)) then
176                  i=id
177                endif
178              endif
179              if(i.gt.nset) then
180                noset(ipos:ipos)=' '
181                write(*,*) '*ERROR reading *BOUNDARY: node set ',
182     &               noset
183                write(*,*) '  has not yet been defined. '
184                call inputerror(inpc,ipoinpc,iline,
185     &               "*EQUATION%",ier)
186                return
187              endif
188              do j=istartset(i),iendset(i)
189                if(ialset(j).gt.0) then
190                  k=ialset(j)
191                  do i1=impcstart,impcend
192                    idof=8*(k-1)+i1
193                    call nident(ikmpc,idof,nmpc,id)
194                    if(id.gt.0) then
195                      if(ikmpc(id).eq.idof) then
196                        impc=ilmpc(id)
197                        call mpcrem(impc,mpcfree,nodempc,
198     &                       nmpc,ikmpc,ilmpc,labmpc,coefmpc,
199     &                       ipompc)
200                        cycle
201                      endif
202                    endif
203                    write(*,*)
204     &                   '*WARNING reading *EQUATION: MPC to remove'
205                    write(*,*) '         is not defined; node:',k
206                    write(*,*) '         degree of freedom:',i1
207                  enddo
208                else
209                  k=ialset(j-2)
210                  do
211                    k=k-ialset(j)
212                    if(k.ge.ialset(j-1)) exit
213                    do i1=impcstart,impcend
214                      idof=8*(k-1)+i1
215                      call nident(ikmpc,idof,nmpc,id)
216                      if(id.gt.0) then
217                        if(ikmpc(id).eq.idof) then
218                          impc=ilmpc(id)
219                          call mpcrem(impc,mpcfree,
220     &                         nodempc,nmpc,ikmpc,ilmpc,labmpc,
221     &                         coefmpc,ipompc)
222                          cycle
223                        endif
224                      endif
225                      write(*,*)
226     &                     '*WARNING reading *EQUATION: MPC to remove'
227                      write(*,*)
228     &                     '         is not defined; node:',k
229                      write(*,*)'         degree of freedom:',i1
230                    enddo
231                  enddo
232                endif
233              enddo
234            endif
235          enddo
236          return
237        else
238          write(*,*)
239     &         '*WARNING reading *EQUATION: parameter not recognized:'
240          write(*,*) '         ',
241     &         textpart(m)(1:index(textpart(m),' ')-1)
242          call inputwarning(inpc,ipoinpc,iline,
243     &         "*EQUATION%")
244        endif
245      enddo
246!
247      if(istep.gt.0) then
248        write(*,*)
249     &       '*ERROR reading *EQUATION: *EQUATION should be placed'
250        write(*,*) '  before all step definitions'
251        ier=1
252        return
253      endif
254!
255      do
256        call getnewline(inpc,textpart,istat,n,key,iline,ipol,inl,
257     &       ipoinp,inp,ipoinpc)
258        if((istat.lt.0).or.(key.eq.1)) return
259        read(textpart(1)(1:10),'(i10)',iostat=istat) nterm
260!
261        nmpc=nmpc+1
262        if(nmpc.gt.nmpc_) then
263          write(*,*) '*ERROR reading *EQUATION: increase nmpc_'
264          ier=1
265          return
266        endif
267!
268        labmpc(nmpc)='                    '
269        ipompc(nmpc)=mpcfree
270        ii=0
271!
272        do
273          call getnewline(inpc,textpart,istat,n,key,iline,ipol,inl,
274     &         ipoinp,inp,ipoinpc)
275          if((istat.lt.0).or.(key.eq.1)) then
276            write(*,*) '*ERROR reading *EQUATION: mpc definition ',
277     &           nmpc
278            write(*,*) '  is not complete. '
279            call inputerror(inpc,ipoinpc,iline,
280     &           "*EQUATION%",ier)
281            return
282          endif
283!
284          do i=1,n/3
285!
286            read(textpart((i-1)*3+1)(1:10),'(i10)',iostat=istat) node
287            if(istat.gt.0) then
288              call inputerror(inpc,ipoinpc,iline,
289     &             "*EQUATION%",ier)
290              return
291            endif
292            if((node.gt.nk).or.(node.le.0)) then
293              write(*,*) '*ERROR reading *EQUATION:'
294              write(*,*) '       node ',node,' is not defined'
295              ier=1
296              return
297            endif
298!
299            read(textpart((i-1)*3+2)(1:10),'(i10)',iostat=istat) ndir
300            if(istat.gt.0) then
301              call inputerror(inpc,ipoinpc,iline,
302     &             "*EQUATION%",ier)
303              return
304            endif
305            if(ndir.le.6) then
306c     elseif(ndir.eq.4) then
307c     ndir=5
308c     elseif(ndir.eq.5) then
309c     ndir=6
310c     elseif(ndir.eq.6) then
311c     ndir=7
312            elseif(ndir.eq.8) then
313              ndir=4
314            elseif(ndir.eq.11) then
315              ndir=0
316            else
317              write(*,*) '*ERROR reading *EQUATION:'
318              write(*,*) '       direction',ndir,' is not defined'
319              ier=1
320              return
321            endif
322!
323            read(textpart((i-1)*3+3)(1:20),'(f20.0)',iostat=istat) x
324            if(istat.gt.0) then
325              call inputerror(inpc,ipoinpc,iline,
326     &             "*EQUATION%",ier)
327              return
328            endif
329!
330!     check whether the node is transformed
331!
332            if(ntrans.le.0) then
333              itr=0
334            elseif(inotr(1,node).eq.0) then
335              itr=0
336            else
337              itr=inotr(1,node)
338            endif
339!
340            if((itr.eq.0).or.(ndir.eq.0).or.(ndir.eq.4)) then
341              nodempc(1,mpcfree)=node
342              nodempc(2,mpcfree)=ndir
343              coefmpc(mpcfree)=x
344!
345!     updating ikmpc and ilmpc
346!
347              if(ii.eq.0) then
348                idof=8*(node-1)+ndir
349                call nident(ikmpc,idof,nmpc-1,id)
350                if(id.gt.0) then
351                  if(ikmpc(id).eq.idof) then
352                    write(*,100)
353     &                   (ikmpc(id))/8+1,ikmpc(id)-8*((ikmpc(id))/8)
354                    ier=1
355                    return
356                  endif
357                endif
358                do j=nmpc,id+2,-1
359                  ikmpc(j)=ikmpc(j-1)
360                  ilmpc(j)=ilmpc(j-1)
361                enddo
362                ikmpc(id+1)=idof
363                ilmpc(id+1)=nmpc
364              endif
365!
366              mpcfreeold=mpcfree
367              mpcfree=nodempc(3,mpcfree)
368              if(mpcfree.eq.0) then
369                write(*,*)
370     &               '*ERROR reading *EQUATION: increase memmpc_'
371                ier=1
372                return
373              endif
374            else
375              call transformatrix(trab(1,inotr(1,node)),
376     &             co(1,node),a)
377!
378              number=ndir-1
379              if(ii.eq.0) then
380!
381!     determining which direction to use for the
382!     dependent side: should not occur on the dependent
383!     side in another MPC and should have a nonzero
384!     coefficient
385!
386                do j=1,3
387                  number=number+1
388                  if(number.gt.3) number=1
389                  idof=8*(node-1)+number
390                  call nident(ikmpc,idof,nmpc-1,id)
391                  if(id.gt.0) then
392                    if(ikmpc(id).eq.idof) then
393                      cycle
394                    endif
395                  endif
396                  if(dabs(a(number,ndir)).lt.1.d-5) cycle
397                  exit
398                enddo
399                if(j.gt.3) then
400                  write(*,*)
401     &                 '*ERROR reading *EQUATION: SPC in node'
402                  write(*,*) node,' in transformed coordinates'
403                  write(*,*) ' cannot be converted in MPC: all'
404                  write(*,*) ' DOFs in the node are used as'
405                  write(*,*) ' dependent nodes in other MPCs'
406                  ier=1
407                  return
408                endif
409                number=number-1
410!
411!     updating ikmpc and ilmpc
412!
413                do j=nmpc,id+2,-1
414                  ikmpc(j)=ikmpc(j-1)
415                  ilmpc(j)=ilmpc(j-1)
416                enddo
417                ikmpc(id+1)=idof
418                ilmpc(id+1)=nmpc
419              endif
420!
421              do j=1,3
422                number=number+1
423                if(number.gt.3) number=1
424                if(dabs(a(number,ndir)).lt.1.d-5) cycle
425                nodempc(1,mpcfree)=node
426                nodempc(2,mpcfree)=number
427                coefmpc(mpcfree)=x*a(number,ndir)
428                mpcfreeold=mpcfree
429                mpcfree=nodempc(3,mpcfree)
430                if(mpcfree.eq.0) then
431                  write(*,*)
432     &                 '*ERROR reading *EQUATION: increase memmpc_'
433                  ier=1
434                  return
435                endif
436              enddo
437            endif
438!
439            ii=ii+1
440          enddo
441!
442          if(ii.eq.nterm) then
443            nodempc(3,mpcfreeold)=0
444            exit
445          endif
446        enddo
447      enddo
448!
449 100  format(/,'*ERROR reading *EQUATION: the DOF corresponding to',
450     &     /,'node ',i10,' in direction',i1,' is detected on',
451     &     /,'the dependent side of two different MPC''s')
452      return
453      end
454
455