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 bounadd(node,is,ie,val,nodeboun,ndirboun,xboun,
20     &  nboun,nboun_,iamboun,iamplitude,nam,ipompc,nodempc,
21     &  coefmpc,nmpc,nmpc_,mpcfree,inotr,trab,
22     &  ntrans,ikboun,ilboun,ikmpc,ilmpc,co,nk,nk_,labmpc,type,
23     &  typeboun,nmethod,iperturb,fixed,vold,nodetrue,mi,label)
24!
25!     adds a boundary condition to the data base
26!
27      implicit none
28!
29      logical fixed,rottracoupling
30!
31      character*1 type,typeboun(*)
32      character*20 labmpc(*),label
33!
34      integer nodeboun(*),ndirboun(*),node,is,ie,nboun,nboun_,i,j,
35     &  iamboun(*),iamplitude,nam,ipompc(*),nodempc(3,*),nmpc,nmpc_,
36     &  mpcfree,inotr(2,*),ntrans,ikboun(*),ilboun(*),ikmpc(*),
37     &  ilmpc(*),itr,idof,newnode,number,id,idofnew,idnew,nk,nk_,
38     &  mpcfreenew,nmethod,iperturb(*),ii,nodetrue,mi(*),three,kflag,
39     &  iy(3),inumber,irotnode(11),irotdof(11)
40!
41      real*8 xboun(*),val,coefmpc(*),trab(7,*),a(3,3),co(3,*),
42     &  vold(0:mi(2),*),dx(3)
43!
44      if(ntrans.le.0) then
45         itr=0
46      elseif(inotr(1,node).eq.0) then
47         itr=0
48      else
49         itr=inotr(1,node)
50      endif
51!
52!     checking for boundary conditions on rotational dofs of
53!     distributing couplings
54!
55      rottracoupling=.false.
56      if((ie.ge.4).and.(ie.le.6)) then
57!
58!        rotational dof
59!
60         do ii=is,ie
61            irotnode(ii)=node
62            if(ii.gt.3) then
63c               idof=8*(node-1)+ii+1
64               idof=8*(node-1)+ii
65               call nident(ikmpc,idof,nmpc,id)
66               if(id.gt.0) then
67                  if(ikmpc(id).eq.idof) then
68                     if(labmpc(ilmpc(id))(1:14).eq.'ROTTRACOUPLING')then
69                        rottracoupling=.true.
70                        irotnode(ii)=
71     &                     nodempc(1,nodempc(3,ipompc(ilmpc(id))))
72                        irotdof(ii)=
73     &                     nodempc(2,nodempc(3,ipompc(ilmpc(id))))
74                        itr=0
75                     endif
76                  endif
77               endif
78            endif
79         enddo
80      endif
81!
82      loop: do ii=is,ie
83!
84!     change: transformations on rotations are taken into account
85!     by the normal of the mean rotation MPC, not by expanding the
86!     MPC in Carthesian coordinates
87!
88         if((itr.eq.0).or.(ii.eq.0).or.(ii.gt.3)) then
89!
90!        no transformation applies: simple SPC
91!
92            if(rottracoupling) then
93               node=irotnode(ii)
94               if(ii.gt.3) then
95                  i=irotdof(ii)
96               else
97                  i=ii
98               endif
99            else
100c               if(ii.le.3) then
101               if(ii.le.6) then
102                  i=ii
103c               elseif(ii.eq.4) then
104c                  i=5
105c               elseif(ii.eq.5) then
106c                  i=6
107c               elseif(ii.eq.6) then
108c                  i=7
109               elseif(ii.eq.8) then
110                  i=4
111               elseif(ii.eq.11) then
112                  i=0
113               else
114                  write(*,*) '*ERROR in bounadd: unknown DOF: ',
115     &                 ii
116                  call exit(201)
117               endif
118            endif
119!
120            if((fixed).and.(i.lt.5)) then
121               val=vold(i,nodetrue)
122            elseif(fixed) then
123               write(*,*) '*ERROR in bounadd: parameter FIXED cannot'
124               write(*,*) '       be used for rotations'
125               call exit(201)
126            endif
127            idof=8*(node-1)+i
128            call nident(ikboun,idof,nboun,id)
129            if(id.gt.0) then
130               if(ikboun(id).eq.idof) then
131                  j=ilboun(id)
132                  if(typeboun(j).ne.type) cycle loop
133                  xboun(j)=val
134                  if(nam.gt.0) iamboun(j)=iamplitude
135                  cycle loop
136               endif
137            endif
138            nboun=nboun+1
139            if(nboun.gt.nboun_) then
140               write(*,*) '*ERROR in bounadd: increase nboun_'
141               call exit(201)
142            endif
143            if((nmethod.eq.4).and.(iperturb(1).le.1)) then
144               write(*,*) '*ERROR in bounadd: in a modal dynamic step'
145               write(*,*) '       new SPCs are not allowed'
146               call exit(201)
147            endif
148            nodeboun(nboun)=node
149            ndirboun(nboun)=i
150            xboun(nboun)=val
151            typeboun(nboun)=type
152            if(nam.gt.0) iamboun(nboun)=iamplitude
153!
154!           updating ikboun and ilboun
155!
156            do j=nboun,id+2,-1
157               ikboun(j)=ikboun(j-1)
158               ilboun(j)=ilboun(j-1)
159            enddo
160            ikboun(id+1)=idof
161            ilboun(id+1)=nboun
162         else
163!
164!        transformation applies: SPC is MPC in global carthesian
165!        coordinates
166!
167            call transformatrix(trab(1,itr),co(1,node),a)
168c            if(ii.le.3) then
169            if(ii.le.6) then
170               i=ii
171c            elseif(ii.eq.4) then
172c               i=5
173c            elseif(ii.eq.5) then
174c               i=6
175c            elseif(ii.eq.6) then
176c               i=7
177            elseif(ii.eq.8) then
178               i=4
179            elseif(ii.eq.11) then
180               i=0
181            else
182               write(*,*) '*ERROR in bounadd: unknown DOF: ',
183     &              ii
184               call exit(201)
185            endif
186            if((fixed).and.(i.lt.5)) then
187               val=vold(i,nodetrue)
188            elseif(fixed) then
189               write(*,*) '*ERROR in bounadd: parameter FIXED cannot'
190               write(*,*) '       be used for rotations'
191               call exit(201)
192            endif
193            if(inotr(2,node).ne.0) then
194               newnode=inotr(2,node)
195               idofnew=8*(newnode-1)+i
196               call nident(ikboun,idofnew,nboun,idnew)
197               if(idnew.gt.0) then
198                  if(ikboun(idnew).eq.idofnew) then
199                     j=ilboun(idnew)
200c
201                     if(typeboun(j).ne.type) cycle
202c
203                     xboun(j)=val
204c                     typeboun(j)=type
205                     if(nam.gt.0) iamboun(j)=iamplitude
206                     cycle
207                  endif
208               endif
209            else
210!
211!              new node is generated for the inhomogeneous MPC term
212!
213               if((nmethod.eq.4).and.(iperturb(1).le.1)) then
214                  write(*,*)'*ERROR in bounadd: in a modal dynamic step'
215                  write(*,*) '       new SPCs are not allowed'
216                  call exit(201)
217               endif
218               nk=nk+1
219               if(nk.gt.nk_) then
220                  write(*,*) '*ERROR in bounadd: increase nk_'
221                  call exit(201)
222               endif
223               newnode=nk
224               inotr(2,node)=newnode
225               idofnew=8*(newnode-1)+i
226               idnew=nboun
227!
228!              copying the initial conditions from node into newnode
229!
230               do j=0,mi(2)
231                  vold(j,newnode)=vold(j,node)
232               enddo
233            endif
234!
235!           new mpc
236!
237            iy(1)=1
238            iy(2)=2
239            iy(3)=3
240            dx(1)=dabs(a(1,i))
241            dx(2)=dabs(a(2,i))
242            dx(3)=dabs(a(3,i))
243            three=3
244            kflag=-2
245            call dsort(dx,iy,three,kflag)
246            do inumber=1,3
247               number=iy(inumber)
248               idof=8*(node-1)+number
249               call nident(ikmpc,idof,nmpc,id)
250               if(id.ne.0) then
251                  if(ikmpc(id).eq.idof) cycle
252               endif
253               if(dabs(a(number,i)).lt.1.d-5) cycle
254               nmpc=nmpc+1
255               if(nmpc.gt.nmpc_) then
256                  write(*,*) '*ERROR in bounadd: increase nmpc_'
257                  call exit(201)
258               endif
259               labmpc(nmpc)=label
260               ipompc(nmpc)=mpcfree
261               do j=nmpc,id+2,-1
262                  ikmpc(j)=ikmpc(j-1)
263                  ilmpc(j)=ilmpc(j-1)
264               enddo
265               ikmpc(id+1)=idof
266               ilmpc(id+1)=nmpc
267               exit
268            enddo
269!
270!           check whether a dependent term was found; if none was
271!           found this can be due to the fact that:
272!           - all dofs were used by other MPC's
273!           - the MPC coefficients were too small
274!           - or a combination of both
275!
276            if(inumber.gt.3) then
277               write(*,*) '*ERROR in bounadd'
278               write(*,*) '       SPC in node',node
279               write(*,*) '       and local direction',ii
280               write(*,*) '       cannot be applied: all'
281               write(*,*) '       degrees of freedom have'
282               write(*,*) '       been used by other MPCs'
283               write(*,*) '       or the coefficient is'
284               write(*,*) '       too small'
285               call exit(201)
286            endif
287!
288            inumber=inumber-1
289            do j=1,3
290               inumber=inumber+1
291               if(inumber.gt.3) inumber=1
292               number=iy(inumber)
293c               if(dabs(a(number,i)).lt.1.d-5) cycle
294               if(dabs(a(number,i)).lt.1.d-30) cycle
295               nodempc(1,mpcfree)=node
296               nodempc(2,mpcfree)=number
297               coefmpc(mpcfree)=a(number,i)
298               mpcfree=nodempc(3,mpcfree)
299               if(mpcfree.eq.0) then
300                  write(*,*) '*ERROR in bounadd: increase memmpc_'
301                  call exit(201)
302               endif
303            enddo
304            nodempc(1,mpcfree)=newnode
305            nodempc(2,mpcfree)=i
306            coefmpc(mpcfree)=-1.d0
307            mpcfreenew=nodempc(3,mpcfree)
308            if(mpcfreenew.eq.0) then
309               write(*,*) '*ERROR in bounadd: increase nmpc_'
310               call exit(201)
311            endif
312            nodempc(3,mpcfree)=0
313            mpcfree=mpcfreenew
314!
315!           nonhomogeneous term
316!
317            nboun=nboun+1
318            if(nboun.gt.nboun_) then
319               write(*,*) '*ERROR in bounadd: increase nboun_'
320               call exit(201)
321            endif
322            nodeboun(nboun)=newnode
323            ndirboun(nboun)=i
324            xboun(nboun)=val
325            typeboun(nboun)=type
326            if(nam.gt.0) iamboun(nboun)=iamplitude
327!
328!           updating ikboun and ilboun
329!
330            do j=nboun,idnew+2,-1
331               ikboun(j)=ikboun(j-1)
332               ilboun(j)=ilboun(j-1)
333            enddo
334            ikboun(idnew+1)=idofnew
335            ilboun(idnew+1)=nboun
336!
337c         enddo
338         endif
339      enddo loop
340!
341      return
342      end
343
344