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 distributingcouplings(inpc,textpart,ipompc,nodempc,
20     &  coefmpc,nmpc,nmpc_,mpcfree,nk,ikmpc,ilmpc,
21     &  labmpc,istep,istat,n,iline,ipol,inl,ipoinp,inp,ipoinpc,lakon,
22     &  kon,ipkon,set,nset,istartset,iendset,ialset,co,ier)
23!
24!     reading the input deck: *DISTRIBUTING COUPLING
25!
26      implicit none
27!
28      character*1 inpc(*)
29      character*8 lakon(*)
30      character*20 labmpc(*)
31      character*81 set(*),elset,noset
32      character*132 textpart(16)
33!
34      integer ipompc(*),nodempc(3,*),nmpc,nmpc_,mpcfree,istep,istat,
35     &  n,i,j,key,nk,node,ier,
36     &  mpcfreeold,ikmpc(*),ilmpc(*),id,idof,iline,ipol,inl,
37     &  ipoinp(2,*),inp(3,*),ipoinpc(0:*),irefnode,
38     &  k,ipos,kon(*),ipkon(*),nset,idir,newmpc,
39     &  istartset(*),iendset(*),ialset(*),indexm,ielem
40!
41      real*8 coefmpc(*),co(3,*),weight,totweight
42!
43      elset(1:1)=' '
44      do i=2,n
45         if(textpart(i)(1:6).eq.'ELSET=') then
46            elset=textpart(i)(7:86)
47            elset(81:81)=' '
48            ipos=index(elset,' ')
49            elset(ipos:ipos)='E'
50         else
51            write(*,*) '*WARNING reading *DISTRIBUTING COUPLING:'
52            write(*,*) '         parameter not recognized:'
53            write(*,*) '         ',
54     &           textpart(i)(1:index(textpart(i),' ')-1)
55            call inputwarning(inpc,ipoinpc,iline,
56     &"*DISTRIBUTING COUPLING%")
57         endif
58      enddo
59!
60      if(elset(1:1).eq.' ') then
61         write(*,*) '*ERROR reading *DISTRIBUTING COUPLING:'
62         write(*,*) '       no element set given'
63         call inputerror(inpc,ipoinpc,iline,
64     &        "*DISTRIBUTING COUPLING%",ier)
65         return
66      endif
67!
68!     check whether the element set exists
69!
70c      do i=1,nset
71c         if(set(i).eq.elset) exit
72c      enddo
73      call cident81(set,elset,nset,id)
74      i=nset+1
75      if(id.gt.0) then
76        if(elset.eq.set(id)) then
77          i=id
78        endif
79      endif
80      if(i.gt.nset) then
81         write(*,*) '*ERROR reading *DISTRIBUTING COUPLING:'
82         write(*,*) '       element set ',elset(1:ipos-1),
83     &         ' is not defined'
84         ier=1
85         return
86      endif
87!
88!     check whether only one element belongs
89!     to the set
90!
91      if(istartset(i).ne.iendset(i)) then
92         write(*,*) '*ERROR reading *DISTRIBUTING COUPLING:'
93         write(*,*) '       element set ',elset(1:ipos-1),
94     &        ' contains more than one element'
95         ier=1
96         return
97      endif
98!
99!     check whether the element is a DCOUP3D element
100!
101      ielem=ialset(istartset(i))
102      if(lakon(ielem)(1:7).ne.'DCOUP3D') then
103         write(*,*) '*ERROR reading *DISTRIBUTING COUPLING:'
104         write(*,*) '       element ',ielem,' is not a'
105         write(*,*) '       DCOUP3D element'
106         ier=1
107         return
108      endif
109!
110!     the reference node belongs to the DCOUP3D element
111!
112      irefnode=kon(ipkon(ielem)+1)
113      newmpc=0
114      totweight=0.d0
115!
116!     generate a MPC for dof 1
117!
118      do
119         call getnewline(inpc,textpart,istat,n,key,iline,ipol,inl,
120     &        ipoinp,inp,ipoinpc)
121         if((istat.lt.0).or.(key.eq.1)) exit
122!
123         read(textpart(1)(1:10),'(i10)',iostat=istat) node
124         if(istat.eq.0) then
125            if(node.gt.nk) then
126               write(*,*) '*ERROR reading *DISTRIBUTING COUPLING:'
127               write(*,*) '       node ',node,' is not defined'
128               ier=1
129               return
130            endif
131!
132!           if first node : new MPC
133!
134            if(newmpc.eq.0) then
135               idof=8*(node-1)+1
136               call nident(ikmpc,idof,nmpc,id)
137               if(id.gt.0) then
138                  if(ikmpc(id).eq.idof) then
139                     write(*,*) '*ERROR reading *DISTRIBUTING COUPLING:'
140                     write(*,*) '       dof 1 of node ',node,
141     &                    ' is already used'
142                     ier=1
143                     return
144                  endif
145               endif
146!
147               nmpc=nmpc+1
148               if(nmpc.gt.nmpc_) then
149                  write(*,*) '*ERROR reading *DISTRIBUTING COUPLING:'
150                  write(*,*) '       increase nmpc_'
151                  ier=1
152                  return
153               endif
154               ipompc(nmpc)=mpcfree
155               labmpc(nmpc)='                    '
156               ipompc(nmpc)=mpcfree
157!
158!              updating ikmpc and ilmpc
159!
160               do j=nmpc,id+2,-1
161                  ikmpc(j)=ikmpc(j-1)
162                  ilmpc(j)=ilmpc(j-1)
163               enddo
164               ikmpc(id+1)=idof
165               ilmpc(id+1)=nmpc
166!
167               newmpc=1
168            endif
169!
170!           reading the weight
171!
172            read(textpart(2)(1:20),'(f20.0)',iostat=istat) weight
173            if(istat.gt.0) then
174               call inputerror(inpc,ipoinpc,iline,
175     &              "*DISTRIBUTING COUPLING%",ier)
176               return
177            endif
178            totweight=totweight+weight
179!
180!           new term in MPC
181!
182            nodempc(1,mpcfree)=node
183            nodempc(2,mpcfree)=1
184            coefmpc(mpcfree)=weight
185            mpcfree=nodempc(3,mpcfree)
186!
187         else
188!
189!           node set
190!
191            read(textpart(1)(1:80),'(a80)',iostat=istat) noset
192            read(textpart(2)(1:20),'(f20.0)',iostat=istat) weight
193            if(istat.gt.0) then
194               call inputerror(inpc,ipoinpc,iline,
195     &              "*DISTRIBUTING COUPLING%",ier)
196               return
197            endif
198            noset(81:81)=' '
199            ipos=index(noset,' ')
200            noset(ipos:ipos)='N'
201c            do i=1,nset
202c               if(set(i).eq.noset) exit
203c            enddo
204            call cident81(set,noset,nset,id)
205            i=nset+1
206            if(id.gt.0) then
207              if(noset.eq.set(id)) then
208                i=id
209              endif
210            endif
211            if(i.gt.nset) then
212               noset(ipos:ipos)=' '
213               write(*,*) '*ERROR reading *DISTRIBUTING COUPLING:'
214               write(*,*) '       node set ',noset
215               write(*,*) '       has not yet been defined. '
216               call inputerror(inpc,ipoinpc,iline,
217     &              "*DISTRIBUTING COUPLING%",ier)
218               return
219            endif
220            do j=istartset(i),iendset(i)
221               if(ialset(j).gt.0) then
222                  node=ialset(j)
223                  totweight=totweight+weight
224!
225                  if(newmpc.eq.0) then
226                     idof=8*(node-1)+1
227                     call nident(ikmpc,idof,nmpc,id)
228                     if(id.gt.0) then
229                        if(ikmpc(id).eq.idof) then
230                           write(*,*)
231     &                       '*ERROR reading *DISTRIBUTING COUPLING:'
232                           write(*,*) '       dof 1 of node ',node,
233     &                          ' is already used'
234                           ier=1
235                           return
236                        endif
237                     endif
238!
239                     nmpc=nmpc+1
240                     if(nmpc.gt.nmpc_) then
241                        write(*,*)
242     &                    '*ERROR reading *DISTRIBUTING COUPLING:'
243                        write(*,*) '       increase nmpc_'
244                        ier=1
245                        return
246                     endif
247                     ipompc(nmpc)=mpcfree
248                     labmpc(nmpc)='                    '
249                     ipompc(nmpc)=mpcfree
250!
251!              updating ikmpc and ilmpc
252!
253                     do k=nmpc,id+2,-1
254                        ikmpc(k)=ikmpc(k-1)
255                        ilmpc(k)=ilmpc(k-1)
256                     enddo
257                     ikmpc(id+1)=idof
258                     ilmpc(id+1)=nmpc
259!
260                     newmpc=1
261                  endif
262!
263!           new term in MPC
264!
265                  nodempc(1,mpcfree)=node
266                  nodempc(2,mpcfree)=1
267                  coefmpc(mpcfree)=weight
268                  mpcfree=nodempc(3,mpcfree)
269!
270               else
271                  node=ialset(j-2)
272                  do
273                     node=node-ialset(j)
274                     if(node.ge.ialset(j-1)) exit
275                     totweight=totweight+weight
276!
277!                    new term in MPC
278!
279                     nodempc(1,mpcfree)=node
280                     nodempc(2,mpcfree)=1
281                     coefmpc(mpcfree)=weight
282                     mpcfree=nodempc(3,mpcfree)
283                  enddo
284               endif
285            enddo
286         endif
287      enddo
288!
289!     reference node
290!
291      nodempc(1,mpcfree)=irefnode
292      nodempc(2,mpcfree)=1
293      coefmpc(mpcfree)=-totweight
294      mpcfreeold=mpcfree
295      mpcfree=nodempc(3,mpcfree)
296      nodempc(3,mpcfreeold)=0
297!
298!     dofs 2 and 3
299!
300      do idir=2,3
301!
302         indexm=ipompc(nmpc)
303         node=nodempc(1,indexm)
304!
305         idof=8*(node-1)+idir
306         call nident(ikmpc,idof,nmpc,id)
307         if(id.gt.0) then
308            if(ikmpc(id).eq.idof) then
309               write(*,*) '*ERROR reading *DISTRIBUTING COUPLING:'
310               write(*,*) '       dof',idir,' of node ',node,
311     &              ' is already used'
312               ier=1
313               return
314            endif
315         endif
316!
317         nmpc=nmpc+1
318         if(nmpc.gt.nmpc_) then
319            write(*,*) '*ERROR reading *DISTRIBUTING COUPLING:'
320            write(*,*) '       increase nmpc_'
321            ier=1
322            return
323         endif
324         ipompc(nmpc)=mpcfree
325         labmpc(nmpc)='                    '
326         ipompc(nmpc)=mpcfree
327!
328!     updating ikmpc and ilmpc
329!
330         do j=nmpc,id+2,-1
331            ikmpc(j)=ikmpc(j-1)
332            ilmpc(j)=ilmpc(j-1)
333         enddo
334         ikmpc(id+1)=idof
335         ilmpc(id+1)=nmpc
336!
337         do
338            nodempc(1,mpcfree)=nodempc(1,indexm)
339            nodempc(2,mpcfree)=idir
340            coefmpc(mpcfree)=coefmpc(indexm)
341            if(nodempc(3,indexm).eq.0) then
342               mpcfreeold=mpcfree
343               mpcfree=nodempc(3,mpcfree)
344               nodempc(3,mpcfreeold)=0
345               exit
346            else
347               mpcfree=nodempc(3,mpcfree)
348               indexm=nodempc(3,indexm)
349            endif
350         enddo
351      enddo
352!
353      return
354      end
355
356