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 mafillsmforc(nforc,ndirforc,nodeforc,xforc,nactdof,
20     &  fext,ipompc,nodempc,coefmpc,mi,rhsi,fnext,
21     &  nmethod,ntrans,inotr,trab,co)
22!
23!     including point forces into the external force vector
24!
25      implicit none
26!
27      integer nforc,ndirforc(*),nodeforc(2,*),mi(*),nactdof(0:mi(2),*),
28     &  ipompc(*),nodempc(3,*),i,jdof,id,ist,
29     &  index,nmethod,rhsi,ntrans,inotr(2,*),node,idir,j,itr
30!
31      real*8 xforc(*),fext(*),coefmpc(*),fnext(0:mi(2),*),trab(7,*),
32     &  a(3,3),co(3,*)
33!
34      if(rhsi.eq.1) then
35!
36!        point forces
37!
38         do i=1,nforc
39            if(ndirforc(i).gt.mi(2)) cycle
40!
41            node=nodeforc(1,i)
42!
43!           check for transformation
44!
45            if(ntrans.eq.0) then
46               itr=0
47            else
48               itr=inotr(1,node)
49            endif
50!
51            if(itr.eq.0) then
52!
53!              no transformation
54!
55!              updating the external force vector for dynamic
56!              calculations
57!
58               if(nmethod.eq.4) fnext(ndirforc(i),nodeforc(1,i))=
59     &              fnext(ndirforc(i),nodeforc(1,i))+xforc(i)
60!
61               jdof=nactdof(ndirforc(i),nodeforc(1,i))
62               if(jdof.gt.0) then
63                  fext(jdof)=fext(jdof)+xforc(i)
64               else
65!
66!                 node is a dependent node of a MPC: distribute
67!                 the forces among the independent nodes
68!                 (proportional to their coefficients)
69!
70                  if(jdof.ne.2*(jdof/2)) then
71                     id=(-jdof+1)/2
72                     ist=ipompc(id)
73                     index=nodempc(3,ist)
74                     if(index.eq.0) cycle
75                     do
76                        jdof=nactdof(nodempc(2,index),nodempc(1,index))
77                        if(jdof.gt.0) then
78                           fext(jdof)=fext(jdof)-
79     &                          coefmpc(index)*xforc(i)/coefmpc(ist)
80                        endif
81                        index=nodempc(3,index)
82                        if(index.eq.0) exit
83                     enddo
84                  endif
85               endif
86            else
87!
88!     transformation
89!
90               call transformatrix(trab(1,itr),co(1,node),a)
91!
92               idir=ndirforc(i)
93!
94!     updating the external force vector for dynamic
95!     calculations
96!
97               if(nmethod.eq.4) then
98                  do j=1,3
99                     fnext(j,node)=fnext(j,node)+a(j,idir)*xforc(i)
100                  enddo
101               endif
102!
103               do j=1,3
104                  jdof=nactdof(j,node)
105                  if(jdof.gt.0) then
106                     fext(jdof)=fext(jdof)+a(j,idir)*xforc(i)
107                  else
108!
109!     node is a dependent node of a MPC: distribute
110!     the forces among the independent nodes
111!     (proportional to their coefficients)
112!
113                     if(jdof.ne.2*(jdof/2)) then
114                        id=(-jdof+1)/2
115                        ist=ipompc(id)
116                        index=nodempc(3,ist)
117                        if(index.eq.0) cycle
118                        do
119                           jdof=nactdof(nodempc(2,index),
120     &                                  nodempc(1,index))
121                           if(jdof.gt.0) then
122                              fext(jdof)=fext(jdof)-
123     &                             coefmpc(index)*a(j,idir)*xforc(i)
124     &                             /coefmpc(ist)
125                           endif
126                           index=nodempc(3,index)
127                           if(index.eq.0) exit
128                        enddo
129                     endif
130                  endif
131               enddo
132            endif
133         enddo
134!
135      endif
136!
137      return
138      end
139
140