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