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 genmpc(inodestet,nnodestet,co,doubleglob,integerglob,
20     &     ipompc,nodempc,coefmpc,nmpc,nmpc_,labmpc,mpcfree,ikmpc,
21     &     ilmpc)
22!
23!     generating MPC's connecting the nodes in the old tet mesh in
24!     which SPC's or point forces were defined with the nodes in the
25!     refined tet mesh
26!
27      implicit none
28!
29      character*20 labmpc(*)
30!
31      integer inodestet(*),nnodestet,integerglob(*),nktet,netet,ne,nkon,
32     &     nfaces,nfield,nselect,imastset,iselect(1),nterms,
33     &     nelem,ialset(1),mpcfreeold,
34     &     iendset(1),istartset(1),konl(20),loopa,
35     &     node,i,j,k,m,ipompc(*),nodempc(3,*),nmpc,nmpc_,mpcfree,
36     &     ikmpc(*),ilmpc(*),idof,id
37!
38      real*8 co(3,*),doubleglob(*),coords(3),ratio(20),value,
39     &     coefmpc(*),dist
40!
41      nktet=integerglob(1)
42      netet=integerglob(2)
43      ne=integerglob(3)
44      nkon=integerglob(4)
45      nfaces=integerglob(5)
46      nfield=0
47      nselect=0
48      imastset=0
49      loopa=8
50!
51      loop1:do i=1,nnodestet
52        node =inodestet(i)
53!
54        do j=1,3
55          coords(j)=co(j,node)
56        enddo
57!
58        call basis(doubleglob(1),doubleglob(netet+1),
59     &       doubleglob(2*netet+1),
60     &       doubleglob(3*netet+1),doubleglob(4*netet+1),
61     &       doubleglob(5*netet+1),integerglob(6),integerglob(netet+6),
62     &       integerglob(2*netet+6),doubleglob(6*netet+1),
63     &       integerglob(3*netet+6),nktet,netet,
64     &       doubleglob(4*nfaces+6*netet+1),nfield,
65     &       doubleglob(13*nktet+4*nfaces+6*netet+1),
66     &       integerglob(7*netet+6),integerglob(ne+7*netet+6),
67     &       integerglob(2*ne+7*netet+6),
68     &       integerglob(nkon+2*ne+7*netet+6),
69     &       coords(1),coords(2),coords(3),value,ratio,iselect,nselect,
70     &       istartset,iendset,ialset,imastset,
71     &       integerglob(nkon+2*ne+8*netet+6),nterms,konl,nelem,loopa,
72     &       dist)
73!
74!     if the old node number was kept in the new mesh it means that
75!     this node was not moved and no MPC's are needed
76!
77        do j=1,nterms
78          if(konl(j).eq.node) then
79c            write(*,*) '*INFO in genmpc: no MPC generated for node',node
80            cycle loop1
81          endif
82        enddo
83!
84!     check whether SPC was applied in the node; if not, the
85!     old node is used as dependent node
86!
87        do j=1,3
88          idof=8*(node-1)+j
89          call nident(ikmpc,idof,nmpc,id)
90!
91!     temporarily more than one MPC with the same idof is allowed
92!
93c          if(id.gt.0) then
94c            if(ikmpc(id).eq.idof) then
95c              write(*,*) '*ERROR in genmpc: a MPC was defined'
96c              write(*,*) '       in the dependent node'
97c              call exit(201)
98c            endif
99c          endif
100!
101          nmpc=nmpc+1
102!
103          if(nmpc.gt.nmpc_) then
104            write(*,*) '*ERROR reading *EQUATION: increase nmpc_'
105            return
106          endif
107!
108!     MPC is characterized by the label RM (refine mesh)
109!
110          labmpc(nmpc)='RM                  '
111          ipompc(nmpc)=mpcfree
112!
113!     updating ikmpc and ilmpc
114!
115          do m=nmpc,id+2,-1
116            ikmpc(m)=ikmpc(m-1)
117            ilmpc(m)=ilmpc(m-1)
118          enddo
119          ikmpc(id+1)=idof
120          ilmpc(id+1)=nmpc
121!
122          nodempc(1,mpcfree)=node
123          nodempc(2,mpcfree)=j
124          coefmpc(mpcfree)=1.d0
125          mpcfree=nodempc(3,mpcfree)
126!
127          do k=1,nterms
128            if(dabs(ratio(k)).gt.1.d-10) then
129              nodempc(1,mpcfree)=konl(k)
130              nodempc(2,mpcfree)=j
131              coefmpc(mpcfree)=-ratio(k)
132              mpcfreeold=mpcfree
133              mpcfree=nodempc(3,mpcfree)
134!
135              if(mpcfree.eq.0) then
136                write(*,*)
137     &               '*ERROR reading *EQUATION: increase memmpc_'
138                return
139              endif
140            endif
141          enddo
142          nodempc(3,mpcfreeold)=0
143        enddo
144        cycle
145      enddo loop1
146!
147      return
148      end
149