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 equationfs(inpc,textpart,ipompc,nodempc,coefmpc, 20 & nmpc,nmpc_,mpcfree,co,trab,ntrans,ikmpc,ilmpc, 21 & labmpc,istep,istat,n,iline,ipol,inl,ipoinp,inp,ipoinpc, 22 & lakon,ne,nload,sideload,ipkon,kon,nelemload,ier) 23! 24! reading the input deck: *EQUATIONF 25! 26 implicit none 27! 28 character*1 inpc(*) 29 character*8 lakon(*) 30 character*20 labmpc(*),label,sideload(*) 31 character*132 textpart(16) 32! 33 integer ipompc(*),nodempc(3,*),nmpc,nmpc_,mpcfree,istep,istat, 34 & n,i,j,ii,key,nterm,number,ntrans,ndir,indexe, 35 & mpcfreeold,ikmpc(*),ilmpc(*),id,idof,itr,iline,ipol,inl, 36 & ipoinp(2,*),inp(3,*),ipoinpc(0:*),ier, 37 & k,m,iface,ifacel,nelem,ifaceq(8,6),ifacet(6,4),ifacew(8,5), 38 & loadid,ne,nope,nopes,ipkon(*),nload,nelemload(2,*),kon(*) 39! 40 real*8 coefmpc(*),co(3,*),trab(7,*),a(3,3),x,cg(3) 41! 42 data ifaceq /4,3,2,1,11,10,9,12, 43 & 5,6,7,8,13,14,15,16, 44 & 1,2,6,5,9,18,13,17, 45 & 2,3,7,6,10,19,14,18, 46 & 3,4,8,7,11,20,15,19, 47 & 4,1,5,8,12,17,16,20/ 48 data ifacet /1,3,2,7,6,5, 49 & 1,2,4,5,9,8, 50 & 2,3,4,6,10,9, 51 & 1,4,3,8,10,7/ 52 data ifacew /1,3,2,9,8,7,0,0, 53 & 4,5,6,10,11,12,0,0, 54 & 1,2,5,4,7,14,10,13, 55 & 2,3,6,5,8,15,11,14, 56 & 4,6,3,1,12,15,9,13/ 57! 58 do m=2,n 59 write(*,*) 60 & '*WARNING reading *EQUATIONF: parameter not recognized:' 61 write(*,*) ' ', 62 & textpart(m)(1:index(textpart(m),' ')-1) 63 call inputwarning(inpc,ipoinpc,iline, 64 &"*EQUATIONF%") 65 enddo 66! 67 if(istep.gt.0) then 68 write(*,*) 69 & '*ERROR reading *EQUATIONF: *EQUATIONF should be placed' 70 write(*,*) ' before all step definitions' 71 ier=1 72 return 73 endif 74! 75 do 76 call getnewline(inpc,textpart,istat,n,key,iline,ipol,inl, 77 & ipoinp,inp,ipoinpc) 78 if((istat.lt.0).or.(key.eq.1)) return 79 read(textpart(1)(1:10),'(i10)',iostat=istat) nterm 80! 81 nmpc=nmpc+1 82 if(nmpc.gt.nmpc_) then 83 write(*,*) '*ERROR reading *EQUATIONF: increase nmpc_' 84 ier=1 85 return 86 endif 87! 88 labmpc(nmpc)='FLUID ' 89 ipompc(nmpc)=mpcfree 90 ii=0 91! 92 do 93 call getnewline(inpc,textpart,istat,n,key,iline,ipol,inl, 94 & ipoinp,inp,ipoinpc) 95 if((istat.lt.0).or.(key.eq.1)) then 96 write(*,*) '*ERROR reading *EQUATIONF: mpc definition ', 97 & nmpc 98 write(*,*) ' is not complete. ' 99 call inputerror(inpc,ipoinpc,iline, 100 & "*EQUATIONF%",ier) 101 return 102 endif 103! 104 do i=1,n/4 105! 106 read(textpart((i-1)*4+1)(1:10),'(i10)',iostat=istat)nelem 107 if(istat.gt.0) then 108 call inputerror(inpc,ipoinpc,iline, 109 & "*EQUATIONF%",ier) 110 return 111 endif 112 if((nelem.gt.ne).or.(nelem.le.0)) then 113 write(*,*) '*ERROR reading *EQUATIONF:' 114 write(*,*) ' element ',nelem,' is not defined' 115 ier=1 116 return 117 endif 118! 119 read(textpart((i-1)*4+2)(2:2),'(i1)',iostat=istat) 120 & ifacel 121 if(istat.gt.0) then 122 call inputerror(inpc,ipoinpc,iline, 123 & "*EQUATIONF%",ier) 124 return 125 endif 126 iface=10*nelem+ifacel 127! 128 read(textpart((i-1)*4+3)(1:10),'(i10)',iostat=istat) ndir 129 if(istat.gt.0) then 130 call inputerror(inpc,ipoinpc,iline, 131 & "*EQUATIONF%",ier) 132 return 133 endif 134 if(ndir.le.3) then 135 elseif(ndir.eq.4) then 136 write(*,*) '*ERROR reading *EQUATIONF: an equation' 137 write(*,*) ' on DOF 4 is not allowed' 138 ier=1 139 return 140 elseif(ndir.eq.5) then 141 write(*,*) '*ERROR reading *EQUATIONF: an equation' 142 write(*,*) ' on DOF 5 is not allowed' 143 ier=1 144 return 145 elseif(ndir.eq.6) then 146 write(*,*) '*ERROR reading *EQUATIONF: an equation' 147 write(*,*) ' on DOF 6 is not allowed' 148 ier=1 149 return 150 elseif(ndir.eq.8) then 151 ndir=4 152 elseif(ndir.eq.11) then 153 ndir=0 154 else 155 write(*,*) '*ERROR reading *EQUATIONF:' 156 write(*,*) ' direction',ndir,' is not defined' 157 ier=1 158 return 159 endif 160! 161 read(textpart((i-1)*4+4)(1:20),'(f20.0)',iostat=istat) x 162 if(istat.gt.0) then 163 call inputerror(inpc,ipoinpc,iline, 164 & "*EQUATIONF%",ier) 165 return 166 endif 167! 168! check whether the face is transformed 169! 170 if(ntrans.le.0) then 171 itr=0 172 else 173 label(1:20)='T ' 174 write(label(2:2),'(i1)') ifacel 175 call identifytransform(nelem,label,nelemload,sideload, 176 & nload,loadid) 177 if(loadid.eq.0) then 178 itr=0 179 else 180 itr=nelemload(2,loadid) 181 endif 182 endif 183! 184 if((itr.eq.0).or.(ndir.eq.0).or.(ndir.eq.4)) then 185 nodempc(1,mpcfree)=iface 186 nodempc(2,mpcfree)=ndir 187 coefmpc(mpcfree)=x 188! 189! updating ikmpc and ilmpc 190! 191 if(ii.eq.0) then 192 idof=-(8*(iface-1)+ndir) 193 call nident(ikmpc,idof,nmpc-1,id) 194 if(id.gt.0) then 195 if(ikmpc(id).eq.idof) then 196 write(*,100) 197 & (ikmpc(id))/8+1,ikmpc(id)-8*((ikmpc(id))/8) 198 ier=1 199 return 200 endif 201 endif 202 do j=nmpc,id+2,-1 203 ikmpc(j)=ikmpc(j-1) 204 ilmpc(j)=ilmpc(j-1) 205 enddo 206 ikmpc(id+1)=idof 207 ilmpc(id+1)=nmpc 208 endif 209! 210 mpcfreeold=mpcfree 211 mpcfree=nodempc(3,mpcfree) 212 if(mpcfree.eq.0) then 213 write(*,*) 214 & '*ERROR reading *EQUATIONF: increase memmpc_' 215 ier=1 216 return 217 endif 218 else 219! 220! transformation applies 221! 222! number of nodes belonging to the face 223! 224 if(lakon(nelem)(4:4).eq.'8') then 225 nope=8 226 nopes=4 227 elseif(lakon(nelem)(4:4).eq.'4') then 228 nope=4 229 nopes=3 230 elseif(lakon(nelem)(4:4).eq.'6') then 231 nope=6 232 if(ifacel.le.2) then 233 nopes=3 234 else 235 nopes=4 236 endif 237 endif 238! 239! determining the center of gravity 240! 241 do j=1,3 242 cg(j)=0.d0 243 enddo 244! 245 indexe=ipkon(nelem) 246 if(nope.eq.8) then 247 do k=1,nopes 248 do j=1,3 249 cg(j)=cg(j)+ 250 & co(j,kon(indexe+ifaceq(k,ifacel))) 251 enddo 252 enddo 253 elseif(nope.eq.4) then 254 do k=1,nopes 255 do j=1,3 256 cg(j)=cg(j)+ 257 & co(j,kon(indexe+ifacet(k,ifacel))) 258 enddo 259 enddo 260 else 261 do k=1,nopes 262 do j=1,3 263 cg(j)=cg(j)+ 264 & co(j,kon(indexe+ifacew(k,ifacel))) 265 enddo 266 enddo 267 endif 268 do j=1,3 269 cg(j)=cg(j)/nopes 270 enddo 271! 272! determining the transformation coefficients at the center 273! of gravity 274! 275 call transformatrix(trab(1,itr), 276 & cg,a) 277! 278 number=ndir-1 279 if(ii.eq.0) then 280! 281! determining which direction to use for the 282! dependent side: should not occur on the dependent 283! side in another MPC and should have a nonzero 284! coefficient 285! 286 do j=1,3 287 number=number+1 288 if(number.gt.3) number=1 289 idof=-(8*(iface-1)+number) 290 call nident(ikmpc,idof,nmpc-1,id) 291 if(id.gt.0) then 292 if(ikmpc(id).eq.idof) then 293 cycle 294 endif 295 endif 296 if(dabs(a(number,ndir)).lt.1.d-5) cycle 297 exit 298 enddo 299 if(j.gt.3) then 300 write(*,*) 301 & '*ERROR reading *EQUATIONF: MPC on face' 302 write(*,*) ifacel,' of element',nelem 303 write(*,*) ' in transformed coordinates' 304 write(*,*) ' cannot be converted in MPC: all' 305 write(*,*) ' DOFs in the node are used as' 306 write(*,*) ' dependent nodes in other MPCs' 307 ier=1 308 return 309 endif 310 number=number-1 311! 312! updating ikmpc and ilmpc 313! 314 do j=nmpc,id+2,-1 315 ikmpc(j)=ikmpc(j-1) 316 ilmpc(j)=ilmpc(j-1) 317 enddo 318 ikmpc(id+1)=idof 319 ilmpc(id+1)=nmpc 320 endif 321! 322 do j=1,3 323 number=number+1 324 if(number.gt.3) number=1 325 if(dabs(a(number,ndir)).lt.1.d-5) cycle 326 nodempc(1,mpcfree)=iface 327 nodempc(2,mpcfree)=number 328 coefmpc(mpcfree)=x*a(number,ndir) 329 mpcfreeold=mpcfree 330 mpcfree=nodempc(3,mpcfree) 331 if(mpcfree.eq.0) then 332 write(*,*) 333 & '*ERROR reading *EQUATIONF: increase memmpc_' 334 ier=1 335 return 336 endif 337 enddo 338 endif 339! 340 ii=ii+1 341 enddo 342! 343 if(ii.eq.nterm) then 344 nodempc(3,mpcfreeold)=0 345 exit 346 endif 347 enddo 348 enddo 349! 350 100 format(/,'*ERROR reading *EQUATIONF: the DOF corresponding to', 351 & /,'iface ',i1,' of element',i10,' in direction', 352 & i5,' is detected on', 353 & /,'the dependent side of two different MPC''s') 354 return 355 end 356 357