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 equationfs(inpc,textpart,ipompc,nodempc,coefmpc,
20     &  nmpc,nmpc_,mpcfree,co,trab,ntrans,ikmpc,ilmpc,
21     &  labmpc,istep,istat,n,iline,ipol,inl,ipoinp,inp,ipoinpc,
22     &  lakon,ne,nload,sideload,ipkon,kon,nelemload,ier)
23!
24!     reading the input deck: *EQUATIONF
25!
26      implicit none
27!
28      character*1 inpc(*)
29      character*8 lakon(*)
30      character*20 labmpc(*),label,sideload(*)
31      character*132 textpart(16)
32!
33      integer ipompc(*),nodempc(3,*),nmpc,nmpc_,mpcfree,istep,istat,
34     &  n,i,j,ii,key,nterm,number,ntrans,ndir,indexe,
35     &  mpcfreeold,ikmpc(*),ilmpc(*),id,idof,itr,iline,ipol,inl,
36     &  ipoinp(2,*),inp(3,*),ipoinpc(0:*),ier,
37     &  k,m,iface,ifacel,nelem,ifaceq(8,6),ifacet(6,4),ifacew(8,5),
38     &  loadid,ne,nope,nopes,ipkon(*),nload,nelemload(2,*),kon(*)
39!
40      real*8 coefmpc(*),co(3,*),trab(7,*),a(3,3),x,cg(3)
41!
42      data ifaceq /4,3,2,1,11,10,9,12,
43     &            5,6,7,8,13,14,15,16,
44     &            1,2,6,5,9,18,13,17,
45     &            2,3,7,6,10,19,14,18,
46     &            3,4,8,7,11,20,15,19,
47     &            4,1,5,8,12,17,16,20/
48      data ifacet /1,3,2,7,6,5,
49     &             1,2,4,5,9,8,
50     &             2,3,4,6,10,9,
51     &             1,4,3,8,10,7/
52      data ifacew /1,3,2,9,8,7,0,0,
53     &             4,5,6,10,11,12,0,0,
54     &             1,2,5,4,7,14,10,13,
55     &             2,3,6,5,8,15,11,14,
56     &             4,6,3,1,12,15,9,13/
57!
58      do m=2,n
59         write(*,*)
60     &        '*WARNING reading *EQUATIONF: parameter not recognized:'
61         write(*,*) '         ',
62     &        textpart(m)(1:index(textpart(m),' ')-1)
63         call inputwarning(inpc,ipoinpc,iline,
64     &"*EQUATIONF%")
65      enddo
66!
67      if(istep.gt.0) then
68         write(*,*)
69     &     '*ERROR reading *EQUATIONF: *EQUATIONF should be placed'
70         write(*,*) '  before all step definitions'
71         ier=1
72         return
73      endif
74!
75      do
76         call getnewline(inpc,textpart,istat,n,key,iline,ipol,inl,
77     &        ipoinp,inp,ipoinpc)
78         if((istat.lt.0).or.(key.eq.1)) return
79         read(textpart(1)(1:10),'(i10)',iostat=istat) nterm
80!
81         nmpc=nmpc+1
82         if(nmpc.gt.nmpc_) then
83            write(*,*) '*ERROR reading *EQUATIONF: increase nmpc_'
84            ier=1
85            return
86         endif
87!
88         labmpc(nmpc)='FLUID               '
89         ipompc(nmpc)=mpcfree
90         ii=0
91!
92         do
93            call getnewline(inpc,textpart,istat,n,key,iline,ipol,inl,
94     &           ipoinp,inp,ipoinpc)
95            if((istat.lt.0).or.(key.eq.1)) then
96               write(*,*) '*ERROR reading *EQUATIONF: mpc definition ',
97     &             nmpc
98               write(*,*) '  is not complete. '
99               call inputerror(inpc,ipoinpc,iline,
100     &              "*EQUATIONF%",ier)
101               return
102            endif
103!
104            do i=1,n/4
105!
106               read(textpart((i-1)*4+1)(1:10),'(i10)',iostat=istat)nelem
107               if(istat.gt.0) then
108                  call inputerror(inpc,ipoinpc,iline,
109     &                 "*EQUATIONF%",ier)
110                  return
111               endif
112               if((nelem.gt.ne).or.(nelem.le.0)) then
113                  write(*,*) '*ERROR reading *EQUATIONF:'
114                  write(*,*) '       element ',nelem,' is not defined'
115                  ier=1
116                  return
117               endif
118!
119               read(textpart((i-1)*4+2)(2:2),'(i1)',iostat=istat)
120     &             ifacel
121               if(istat.gt.0) then
122                  call inputerror(inpc,ipoinpc,iline,
123     &                 "*EQUATIONF%",ier)
124                  return
125               endif
126               iface=10*nelem+ifacel
127!
128               read(textpart((i-1)*4+3)(1:10),'(i10)',iostat=istat) ndir
129               if(istat.gt.0) then
130                  call inputerror(inpc,ipoinpc,iline,
131     &                 "*EQUATIONF%",ier)
132                  return
133               endif
134               if(ndir.le.3) then
135               elseif(ndir.eq.4) then
136                  write(*,*) '*ERROR reading *EQUATIONF: an equation'
137                  write(*,*) '       on DOF 4 is not allowed'
138                  ier=1
139                  return
140               elseif(ndir.eq.5) then
141                  write(*,*) '*ERROR reading *EQUATIONF: an equation'
142                  write(*,*) '       on DOF 5 is not allowed'
143                  ier=1
144                  return
145               elseif(ndir.eq.6) then
146                  write(*,*) '*ERROR reading *EQUATIONF: an equation'
147                  write(*,*) '       on DOF 6 is not allowed'
148                  ier=1
149                  return
150               elseif(ndir.eq.8) then
151                  ndir=4
152               elseif(ndir.eq.11) then
153                  ndir=0
154               else
155                  write(*,*) '*ERROR reading *EQUATIONF:'
156                  write(*,*) '       direction',ndir,' is not defined'
157                  ier=1
158                  return
159               endif
160!
161               read(textpart((i-1)*4+4)(1:20),'(f20.0)',iostat=istat) x
162               if(istat.gt.0) then
163                  call inputerror(inpc,ipoinpc,iline,
164     &                 "*EQUATIONF%",ier)
165                  return
166               endif
167!
168!              check whether the face is transformed
169!
170               if(ntrans.le.0) then
171                  itr=0
172               else
173                  label(1:20)='T                   '
174                  write(label(2:2),'(i1)') ifacel
175                  call identifytransform(nelem,label,nelemload,sideload,
176     &                 nload,loadid)
177                  if(loadid.eq.0) then
178                     itr=0
179                  else
180                     itr=nelemload(2,loadid)
181                  endif
182               endif
183!
184               if((itr.eq.0).or.(ndir.eq.0).or.(ndir.eq.4)) then
185                  nodempc(1,mpcfree)=iface
186                  nodempc(2,mpcfree)=ndir
187                  coefmpc(mpcfree)=x
188!
189!                 updating ikmpc and ilmpc
190!
191                  if(ii.eq.0) then
192                     idof=-(8*(iface-1)+ndir)
193                     call nident(ikmpc,idof,nmpc-1,id)
194                     if(id.gt.0) then
195                        if(ikmpc(id).eq.idof) then
196                           write(*,100)
197     &                   (ikmpc(id))/8+1,ikmpc(id)-8*((ikmpc(id))/8)
198                           ier=1
199                           return
200                        endif
201                     endif
202                     do j=nmpc,id+2,-1
203                        ikmpc(j)=ikmpc(j-1)
204                        ilmpc(j)=ilmpc(j-1)
205                     enddo
206                     ikmpc(id+1)=idof
207                     ilmpc(id+1)=nmpc
208                  endif
209!
210                  mpcfreeold=mpcfree
211                  mpcfree=nodempc(3,mpcfree)
212                  if(mpcfree.eq.0) then
213                     write(*,*)
214     &                 '*ERROR reading *EQUATIONF: increase memmpc_'
215                     ier=1
216                     return
217                  endif
218               else
219!
220!        transformation applies
221!
222!        number of nodes belonging to the face
223!
224                  if(lakon(nelem)(4:4).eq.'8') then
225                     nope=8
226                     nopes=4
227                  elseif(lakon(nelem)(4:4).eq.'4') then
228                     nope=4
229                     nopes=3
230                  elseif(lakon(nelem)(4:4).eq.'6') then
231                     nope=6
232                     if(ifacel.le.2) then
233                        nopes=3
234                     else
235                        nopes=4
236                     endif
237                  endif
238!
239!     determining the center of gravity
240!
241                  do j=1,3
242                     cg(j)=0.d0
243                  enddo
244!
245                  indexe=ipkon(nelem)
246                  if(nope.eq.8) then
247                     do k=1,nopes
248                        do j=1,3
249                           cg(j)=cg(j)+
250     &                           co(j,kon(indexe+ifaceq(k,ifacel)))
251                        enddo
252                     enddo
253                  elseif(nope.eq.4) then
254                     do k=1,nopes
255                        do j=1,3
256                           cg(j)=cg(j)+
257     &                           co(j,kon(indexe+ifacet(k,ifacel)))
258                        enddo
259                     enddo
260                  else
261                     do k=1,nopes
262                        do j=1,3
263                           cg(j)=cg(j)+
264     &                           co(j,kon(indexe+ifacew(k,ifacel)))
265                        enddo
266                     enddo
267                  endif
268                  do j=1,3
269                     cg(j)=cg(j)/nopes
270                  enddo
271!
272!     determining the transformation coefficients at the center
273!     of gravity
274!
275                  call transformatrix(trab(1,itr),
276     &                 cg,a)
277!
278                  number=ndir-1
279                  if(ii.eq.0) then
280!
281!                    determining which direction to use for the
282!                    dependent side: should not occur on the dependent
283!                    side in another MPC and should have a nonzero
284!                    coefficient
285!
286                     do j=1,3
287                        number=number+1
288                        if(number.gt.3) number=1
289                        idof=-(8*(iface-1)+number)
290                        call nident(ikmpc,idof,nmpc-1,id)
291                        if(id.gt.0) then
292                           if(ikmpc(id).eq.idof) then
293                              cycle
294                           endif
295                        endif
296                        if(dabs(a(number,ndir)).lt.1.d-5) cycle
297                        exit
298                     enddo
299                     if(j.gt.3) then
300                        write(*,*)
301     &                       '*ERROR reading *EQUATIONF: MPC on face'
302                        write(*,*) ifacel,' of element',nelem
303                        write(*,*) ' in transformed coordinates'
304                        write(*,*) ' cannot be converted in MPC: all'
305                        write(*,*) ' DOFs in the node are used as'
306                        write(*,*) ' dependent nodes in other MPCs'
307                        ier=1
308                        return
309                     endif
310                     number=number-1
311!
312!                    updating ikmpc and ilmpc
313!
314                     do j=nmpc,id+2,-1
315                        ikmpc(j)=ikmpc(j-1)
316                        ilmpc(j)=ilmpc(j-1)
317                     enddo
318                     ikmpc(id+1)=idof
319                     ilmpc(id+1)=nmpc
320                  endif
321!
322                  do j=1,3
323                     number=number+1
324                     if(number.gt.3) number=1
325                     if(dabs(a(number,ndir)).lt.1.d-5) cycle
326                     nodempc(1,mpcfree)=iface
327                     nodempc(2,mpcfree)=number
328                     coefmpc(mpcfree)=x*a(number,ndir)
329                     mpcfreeold=mpcfree
330                     mpcfree=nodempc(3,mpcfree)
331                     if(mpcfree.eq.0) then
332                        write(*,*)
333     &                    '*ERROR reading *EQUATIONF: increase memmpc_'
334                        ier=1
335                        return
336                     endif
337                  enddo
338               endif
339!
340               ii=ii+1
341            enddo
342!
343            if(ii.eq.nterm) then
344               nodempc(3,mpcfreeold)=0
345               exit
346            endif
347         enddo
348      enddo
349!
350 100  format(/,'*ERROR reading *EQUATIONF: the DOF corresponding to',
351     &           /,'iface ',i1,' of element',i10,' in direction',
352     &           i5,' is detected on',
353     &           /,'the dependent side of two different MPC''s')
354      return
355      end
356
357