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 applympc_dpel(nface,ielfa,xrlfa,vel,vfa,
20     &  ifabou,xbounact,nef,gradpcel,gradpcfa,neifa,rf,area,volume,
21     &  xle,xxi,icyclic,xxn,ipnei,ifatie,
22     &  coefmpc,nmpc,labmpc,ipompc,nodempc,ifaext,nfaext,nactdoh,
23     &  iflag,xxj,xlet)
24!
25!     applying the MPC's in the calculation of the pressure
26!     correction
27!
28      implicit none
29!
30      character*20 labmpc(*)
31!
32      integer nface,ielfa(4,*),ifabou(*),i,iel1,iel2,nef,ibou,
33     &  neifa(*),icyclic,ifa,indexf,l,m,ipnei(*),ifatie(*),
34     &  is,ie,nmpc,ipompc(*),nodempc(3,*),ifaext(*),nfaext,nactdoh(*),
35     &  iel,index,mpc,ipointer,k,ielorig,iface,iflag
36!
37      real*8 xrlfa(3,*),vel(nef,0:7),vfa(0:7,*),xbounact(*),xl1,xl2,
38     &   vfap(0:7,nface),gradpcel(3,*),gradpcfa(3,*),rf(3,*),area(*),
39     &   volume(*),xle(*),xxi(3,*),c(3,3),gradnor,xxn(3,*),coefmpc(*),
40     &   coefnorm,sum,xxj(3,*),xlet(*),dd
41!
42!
43!
44!     Multiple point constraints
45!
46      if(nmpc.gt.0) then
47         do k=1,nfaext
48            i=ifaext(k)
49            if(ielfa(2,i).ge.0) cycle
50            ipointer=-ielfa(2,i)
51            if(ifabou(ipointer+4).ge.0) cycle
52            mpc=-ifabou(ipointer+4)
53!
54            index=ipompc(mpc)
55            sum=0.d0
56            coefnorm=0.d0
57            do
58               if(index.eq.0) exit
59               if(nodempc(1,index).lt.0) then
60!
61!                 a negative number refers to a boundary
62!                 condition (fields nodeboun, ndirboun..)
63!                 resulting from a SPC in local coordinates
64!
65                  sum=sum+0.d0
66               else
67!
68!                 face term
69!
70                  ielorig=int(nodempc(1,index)/10.d0)
71                  iel=nactdoh(ielorig)
72                  iface=nodempc(1,index)-10*ielorig
73                  sum=sum+coefmpc(index)
74     &                 *vfa(nodempc(2,index),neifa(ipnei(iel)+iface))
75                  coefnorm=coefnorm+coefmpc(index)**2
76               endif
77               index=nodempc(3,index)
78            enddo
79!
80!           distribute the sum across all terms which are not
81!           fixed by boundary conditions
82!
83            index=ipompc(mpc)
84            do
85               if(index.eq.0) exit
86               if(nodempc(1,index).gt.0) then
87                  ielorig=int(nodempc(1,index)/10.d0)
88                  iel=nactdoh(ielorig)
89                  iface=nodempc(1,index)-10*ielorig
90                  vfa(nodempc(2,index),neifa(ipnei(iel)+iface))=
91     &                 vfa(nodempc(2,index),neifa(ipnei(iel)+iface))-
92     &                 sum*coefmpc(index)/coefnorm
93               endif
94               index=nodempc(3,index)
95            enddo
96         enddo
97      endif
98!
99      return
100      end
101