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 modifympc(inodestet,nnodestet,co,doubleglob,
20     &     integerglob,ipompc,nodempc,coefmpc,nmpc,nmpc_,labmpc,mpcfree,
21     &     ikmpc,ilmpc,jq,irow,icol,loc,irowt,jqt,itemp,au,ixcol,
22     &     ikboun,nboun,nodeboun,mpcrfna,mpcrfnb,nodempcref,
23     &     coefmpcref,memmpcref_,mpcfreeref,maxlenmpcref,memmpc_,
24     &     maxlenmpc,istep)
25!
26!     generating MPC's connecting the nodes in the old tet mesh in
27!     which SPC's or point forces were defined with the nodes in the
28!     refined tet mesh
29!
30      implicit none
31!
32      logical spc
33!
34      character*20 labmpc(*)
35!
36      integer inodestet(*),nnodestet,integerglob(*),nktet,netet,ne,nkon,
37     &     nfaces,nfield,nselect,imastset,iselect(1),nterms,idir1,
38     &     nelem,ialset(1),mpcfreeold,jq(*),irow(*),icol(*),
39     &     iendset(1),istartset(1),konl(20),loopa,loc(*),irowt(*),
40     &     node,i,j,k,m,ipompc(*),nodempc(3,*),nmpc,nmpc_,mpcfree,
41     &     ikmpc(*),ilmpc(*),idof,id,jqt(*),nzs,nosol,nodeboun(*),
42     &     kflag,itemp(*),maxrow,ixcol(*),isoltot,newnode,iy,ic,idummy,
43     &     nlength,kopt,ixcolnew,idnew,i1,ikboun(*),nboun,mpcrfna,
44     &     mpcrfnb,index,index1,node1,mpc,idof1,idof2,id1,id2,
45     &     nodempcref(3,*),memmpcref_,mpcfreeref,maxlenmpcref,memmpc_,
46     &     maxlenmpc,istep
47!
48      real*8 co(3,*),doubleglob(*),coords(3),ratio(20),value,au(*),
49     &     coefmpc(*),depcoef,dist,xopt,coef1,coefmpcref(*)
50!
51      if(istep.gt.1) then
52!
53!     restoring the original equations
54!
55        memmpc_=memmpcref_
56        mpcfree=mpcfreeref
57        maxlenmpc=maxlenmpcref
58!
59!     mpcfreeref=-1 indicates that the MPC's have changed and that
60!     a new copy is to be taken before the next call to cascade.c
61!
62        mpcfreeref=-1
63!
64        do i=1,memmpc_
65          do j=1,3
66            nodempc(j,i)=nodempcref(j,i)
67          enddo
68          coefmpc(i)=coefmpcref(i)
69        enddo
70      endif
71!
72!     restore the original form of the MPC's
73!     an MPC with number in between mpcrfna and mpcrfnb (i.e. MPC's
74!     corresponding to the connection of the old mesh with the new
75!     mesh) is in its orignal form if the coefficient of the first
76!     term is exactly 1.d0
77!
78      do i=mpcrfna,mpcrfnb
79        index=ipompc(i)
80        if(coefmpc(index).eq.1.d0) cycle
81!
82!     MPC has been reordered and has to be brought in its original
83!     form
84!
85        node1=nodempc(1,index)
86        idir1=nodempc(2,index)
87        coef1=coefmpc(index)
88        idof1=8*(node1-1)+idir1
89        index1=index
90        call nident(ikmpc,idof1,nmpc,id1)
91!
92        index=nodempc(3,index)
93!
94        do
95          if(index.eq.0) then
96            write(*,*) '*ERROR in modifympc: equation ',i
97            write(*,*) '       cannot be modified'
98            call exit(201)
99          endif
100          if(coefmpc(index).eq.1.d0) then
101!
102!     restore the leading term of the equation
103!
104            nodempc(1,index1)=nodempc(1,index)
105            coefmpc(index1)=coefmpc(index)
106!
107!     restore the independent term of the equation
108!     id1>id2 since the refined mesh has higher node numbers
109!     than the unrefined mesh
110!
111            idof2=8*(nodempc(1,index)-1)+nodempc(2,index)
112            call nident(ikmpc,idof2,nmpc,id2)
113            do j=id1,id2+2,-1
114              ikmpc(j)=ikmpc(j-1)
115              ilmpc(j)=ilmpc(j-1)
116            enddo
117            ikmpc(id2+1)=idof2
118            ilmpc(id2+1)=i
119!
120            nodempc(1,index)=node1
121            coefmpc(index)=coef1
122            exit
123          else
124            index=nodempc(3,index)
125          endif
126        enddo
127      enddo
128!
129!     compile all dependent nodes in jq, irow..
130!     taking the equations for the first dof (x-direction) only
131!
132      nzs=0
133      nnodestet=0
134      maxrow=0
135!
136      do i=mpcrfna,mpcrfnb
137        index=ipompc(i)
138!
139!       x-direction
140!
141        if(nodempc(2,index).ne.1) cycle
142!
143        nnodestet=nnodestet+1
144        inodestet(nnodestet)=nodempc(1,index)
145        jq(nnodestet)=nzs+1
146        index=nodempc(3,index)
147        do
148          if(index.eq.0) exit
149          nzs=nzs+1
150          irow(nzs)=nodempc(1,index)
151          maxrow=max(maxrow,irow(nzs))
152          au(nzs)=coefmpc(index)
153          irowt(nzs)=nnodestet
154          loc(nzs)=nzs
155          index=nodempc(3,index)
156        enddo
157      enddo
158      jq(nnodestet+1)=nzs+1
159!
160!     determine the nodes for which the dependent term in the MPC's
161!     has to be modified: these are MPC's for which the dependent term
162!     is subject to other SPC's or MPC's
163!
164!     SPC's
165!
166      do i=1,nboun
167        node=nodeboun(i)
168        call nident(inodestet,node,nnodestet,id)
169        if(id.gt.0) then
170          if(inodestet(id).eq.node) then
171            icol(id)=1
172          endif
173        endif
174      enddo
175!
176!     MPC's
177!
178      do i=1,mpcrfna-1
179        index=ipompc(i)
180        node=nodempc(1,index)
181        call nident(inodestet,node,nnodestet,id)
182        if(id.gt.0) then
183          if(inodestet(id).eq.node) then
184            icol(id)=1
185          endif
186        endif
187      enddo
188!
189!     for the nodes i for which the equations have to be modified
190!     icol(i) is nonzero
191!
192      isoltot=0
193      do i=1,nnodestet
194        if(icol(i).eq.1) then
195          icol(i)=jq(i+1)-jq(i)
196        else
197          isoltot=isoltot+1
198        endif
199      enddo
200!
201!     copying irow into itemp, sorting itemp and field loc and
202!     irowt along
203!
204      do j=1,nzs
205        itemp(j)=irow(j)
206      enddo
207!
208      kflag=2
209      call isortiii(itemp,loc,irowt,nzs,kflag)
210!
211!     determining jqt (columns per row); jqt and irowt is the
212!     equivalent (for row per row treatment) of jq and irow (for
213!     column per column treatment); one can consider jqt and irowt
214!     to be the jq and irow fields for the transpose of the matrix
215!     (therefore the appended letter "t")
216!
217      j=1
218      do m=1,nzs
219        if(itemp(m).ge.j) then
220          do k=j,itemp(m)
221            jqt(k)=m
222          enddo
223          j=itemp(m)+1
224        endif
225      enddo
226      jqt(j)=nzs+1
227!
228!     creating a unique number consisting of the column number and the
229!     number of rows per column; if sorted, the result is the same as
230!     for sorting icol.
231!
232      do j=1,nnodestet
233        if(icol(j).gt.0) then
234          ixcol(j)=(nnodestet+1)*icol(j)+j
235        endif
236      enddo
237!
238!     sorting ixcol
239!
240      kflag=1
241      call isortii(ixcol,idummy,nnodestet,kflag)
242!
243!     reordering the equations: the first term is to be a
244!     new node instead of an old node; starting with columns
245!     with only 1 nonzero value
246!
247      nosol=0
248      do
249        if(isoltot.eq.nnodestet) exit
250        isoltot=isoltot+1
251!
252!       local dependent node number
253!
254        i=ixcol(isoltot)
255     &       -int(ixcol(isoltot)/(nnodestet+1))*(nnodestet+1)
256        ixcol(isoltot)=0
257        icol(i)=0
258!
259!       global dependent node number
260!
261        node=inodestet(i)
262!
263!       determining the optimal row for the exchange
264!
265        kopt=0
266        xopt=0.d0
267        do k=jq(i),jq(i+1)-1
268          if(irow(k).gt.0) then
269            if(dabs(au(k)).gt.xopt) then
270              kopt=k
271              xopt=dabs(au(k))
272            endif
273          endif
274        enddo
275        k=kopt
276!
277        if(k.gt.0) then
278!
279!     row dof
280!
281          newnode=irow(k)
282          depcoef=au(k)
283!
284!         generating MPC's
285!
286          do j=1,3
287!
288!     modify ikmpc to take the switch from id2 to id1 as
289!     dependent node into account
290!
291            idof1=8*(node-1)+j
292            call nident(ikmpc,idof1,nmpc,id1)
293            do
294              if(labmpc(ilmpc(id1)).eq.'RM                  ') exit
295              id=id-1
296            enddo
297            mpc=ilmpc(id1)
298            index=ipompc(mpc)
299            idof2=8*(newnode-1)+j
300            call nident(ikmpc,idof2,nmpc,id2)
301            do m=id1,id2-1
302              ikmpc(m)=ikmpc(m+1)
303              ilmpc(m)=ilmpc(m+1)
304            enddo
305            ikmpc(id2)=idof2
306            ilmpc(id2)=mpc
307!
308!     inserting the new dependent term
309!
310            nodempc(1,index)=newnode
311            coefmpc(index)=depcoef
312!
313!     inserting the new independent term
314!
315            do
316              index=nodempc(3,index)
317              if(index.eq.0) then
318                write(*,*) '*ERROR in modifympc: equation ',mpc
319                write(*,*) '       cannot be modified'
320                call exit(201)
321              endif
322              if(nodempc(1,index).eq.newnode) then
323                nodempc(1,index)=node
324                coefmpc(index)=1.d0
325                exit
326              endif
327            enddo
328          enddo
329!
330!     remove the indpendent node from the data base:
331!     tagging the row dof in all other columns
332!
333          do m=jqt(newnode),jqt(newnode+1)-1
334            if(irowt(m).eq.i) cycle
335            ic=irowt(m)
336!
337!           no SPC/MPC or already treated
338!
339            if(icol(ic).eq.0) cycle
340!
341!     setting negative sign in irow to show that the row
342!     is not available any more (already used as dependent dof)
343!
344            irow(loc(m))=-irow(loc(m))
345            iy=(nnodestet+1)*icol(ic)+ic
346            nlength=nnodestet-isoltot
347            call nident(ixcol(isoltot+1),iy,nlength,id)
348            if(ixcol(isoltot+id).ne.iy) then
349              write(*,*) '*ERROR in modifympc:'
350              write(*,*) '       data base corrupt'
351              call exit(201)
352            endif
353            icol(ic)=icol(ic)-1
354            ixcolnew=ixcol(isoltot+id)-(nnodestet+1)
355!
356!           not treated equation has no independent terms any more
357!
358            if(ixcolnew.lt.(nnodestet+1)) then
359              nosol=nosol+1
360              if(nosol.eq.1) then
361                write(*,*) '*WARNING in modifympc; failed to connect'
362                write(*,*) '         node ',inodestet(ic),
363     &               ' of unrefined mesh'
364                write(*,*)
365     &'         to the refined mesh; loads and boundary conditions'
366                write(*,*)
367     &'         in this node are not taken into account; a list of'
368                write(*,*)
369     &'         not connected nodes is stored in'
370                write(*,*)
371     &'         WarnNodeMissRefineConnection.nam'
372              else
373                write(*,*) '*WARNING in modifympc; failed to connect'
374                write(*,*) '         node ',inodestet(ic),
375     &               ' of unrefined mesh'
376              endif
377              write(23,*)
378     &             '*NSET,NSET=WarnNodeMissRefineConnection'
379              write(23,*) inodestet(ic)
380            endif
381            call nident(ixcol(isoltot+1),ixcolnew,id-1,idnew)
382            do i1=isoltot+id,isoltot+idnew+2,-1
383              ixcol(i1)=ixcol(i1-1)
384            enddo
385            ixcol(isoltot+idnew+1)=ixcolnew
386          enddo
387        endif
388      enddo
389      write(*,*)
390      write(*,*) '*INFO in modifympc: no solution for ',nosol,' nodes.'
391      write(*,*)
392      close(23)
393!
394      return
395      end
396