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 equations(inpc,textpart,ipompc,nodempc,coefmpc, 20 & nmpc,nmpc_,mpcfree,nk,co,trab,inotr,ntrans,ikmpc,ilmpc, 21 & labmpc,istep,istat,n,iline,ipol,inl,ipoinp,inp,ipoinpc, 22 & set,istartset,iendset,ialset,nset,nodempcref,coefmpcref, 23 & ikmpcref,memmpcref_,mpcfreeref,maxlenmpcref,memmpc_, 24 & maxlenmpc,ier) 25! 26! reading the input deck: *EQUATION 27! 28 implicit none 29! 30 character*1 inpc(*) 31 character*20 labmpc(*) 32 character*81 set(*),noset 33 character*132 textpart(16) 34! 35 integer ipompc(*),nodempc(3,*),nmpc,nmpc_,mpcfree,istep,istat, 36 & n,i,j,ii,key,nterm,number,nk,inotr(2,*),ntrans,node,ndir, 37 & mpcfreeold,ikmpc(*),ilmpc(*),id,idof,itr,iline,ipol,inl, 38 & ipoinp(2,*),inp(3,*),ipoinpc(0:*),impcstart,impcend,i1, 39 & istartset(*),iendset(*),ialset(*),nset,k,l,m,index1,ipos, 40 & impc,nodempcref(3,*),ikmpcref(*),memmpcref_,mpcfreeref, 41 & maxlenmpcref,memmpc_,maxlenmpc,ier 42! 43 real*8 coefmpc(*),co(3,*),trab(7,*),a(3,3),x,coefmpcref(*) 44! 45 do m=2,n 46 if(textpart(m)(1:9).eq.'REMOVEALL') then 47! 48 if(istep.eq.1) then 49 write(*,*) '*ERROR reading *EQUATION' 50 write(*,*) ' removing equations is not allowed' 51 write(*,*) ' in the first step' 52 ier=1 53 return 54 endif 55! 56 do j=1,nmpc 57 index1=ipompc(j) 58 if(index1.eq.0) cycle 59 do 60 if(nodempc(3,index1).eq.0) then 61 nodempc(3,index1)=mpcfree 62 mpcfree=ipompc(j) 63 exit 64 endif 65 index1=nodempc(3,index1) 66 enddo 67 ipompc(j)=0 68 ikmpc(j)=0 69 ilmpc(j)=0 70 enddo 71 nmpc=0 72 call getnewline(inpc,textpart,istat,n,key,iline,ipol,inl, 73 & ipoinp,inp,ipoinpc) 74 return 75 elseif(textpart(m)(1:6).eq.'REMOVE') then 76! 77 if(istep.eq.1) then 78 write(*,*) '*ERROR reading *EQUATION' 79 write(*,*) ' removing equations is not allowed' 80 write(*,*) ' in the first step' 81 ier=1 82 return 83 endif 84! 85 do i=1,nmpc 86 if(ikmpcref(i).ne.ikmpc(i)) then 87 write(*,*) '*ERROR reading *EQUATION' 88 write(*,*) ' The dependent terms in some' 89 write(*,*) ' of the nonlinear equations have' 90 write(*,*) ' changed since the start of the' 91 write(*,*) ' calculation. Removing equations' 92 write(*,*) ' does not work' 93 ier=1 94 return 95 endif 96 enddo 97! 98! restoring the original equations (before the first call to 99! cascade) 100! 101 memmpc_=memmpcref_ 102 mpcfree=mpcfreeref 103 maxlenmpc=maxlenmpcref 104! 105! mpcfreeref=-1 indicates that the MPC's have changed and that 106! a new copy is to be taken before the next call to cascade.c 107! 108 mpcfreeref=-1 109! 110 do i=1,memmpc_ 111 do j=1,3 112 nodempc(j,i)=nodempcref(j,i) 113 enddo 114 coefmpc(i)=coefmpcref(i) 115 enddo 116! 117 do 118 call getnewline(inpc,textpart,istat,n,key,iline,ipol,inl, 119 & ipoinp,inp,ipoinpc) 120 if((istat.lt.0).or.(key.eq.1)) return 121! 122 read(textpart(2)(1:10),'(i10)',iostat=istat) impcstart 123 if(istat.gt.0) then 124 call inputerror(inpc,ipoinpc,iline, 125 & "*EQUATION%",ier) 126 return 127 endif 128! 129 if(textpart(3)(1:1).eq.' ') then 130 impcend=impcstart 131 else 132 read(textpart(3)(1:10),'(i10)',iostat=istat) impcend 133 if(istat.gt.0) then 134 call inputerror(inpc,ipoinpc,iline, 135 & "*EQUATION%",ier) 136 return 137 endif 138 endif 139! 140 read(textpart(1)(1:10),'(i10)',iostat=istat) l 141 if(istat.eq.0) then 142 if((l.gt.nk).or.(l.le.0)) then 143 write(*,*) '*ERROR reading *BOUNDARY:' 144 write(*,*) ' node ',l,' is not defined' 145 ier=1 146 return 147 endif 148 do i1=impcstart,impcend 149 idof=8*(l-1)+i1 150 call nident(ikmpc,idof,nmpc,id) 151 if(id.gt.0) then 152 if(ikmpc(id).eq.idof) then 153 impc=ilmpc(id) 154 call mpcrem(impc,mpcfree,nodempc,nmpc, 155 & ikmpc,ilmpc,labmpc,coefmpc,ipompc) 156 cycle 157 endif 158 endif 159 write(*,*) 160 & '*WARNING reading *EQUATION: MPC to remove' 161 write(*,*) ' is not defined; node:',l 162 write(*,*) ' degree of freedom:',i1 163 enddo 164 else 165 read(textpart(1)(1:80),'(a80)',iostat=istat) noset 166 noset(81:81)=' ' 167 ipos=index(noset,' ') 168 noset(ipos:ipos)='N' 169c do i=1,nset 170c if(set(i).eq.noset) exit 171c enddo 172 call cident81(set,noset,nset,id) 173 i=nset+1 174 if(id.gt.0) then 175 if(noset.eq.set(id)) then 176 i=id 177 endif 178 endif 179 if(i.gt.nset) then 180 noset(ipos:ipos)=' ' 181 write(*,*) '*ERROR reading *BOUNDARY: node set ', 182 & noset 183 write(*,*) ' has not yet been defined. ' 184 call inputerror(inpc,ipoinpc,iline, 185 & "*EQUATION%",ier) 186 return 187 endif 188 do j=istartset(i),iendset(i) 189 if(ialset(j).gt.0) then 190 k=ialset(j) 191 do i1=impcstart,impcend 192 idof=8*(k-1)+i1 193 call nident(ikmpc,idof,nmpc,id) 194 if(id.gt.0) then 195 if(ikmpc(id).eq.idof) then 196 impc=ilmpc(id) 197 call mpcrem(impc,mpcfree,nodempc, 198 & nmpc,ikmpc,ilmpc,labmpc,coefmpc, 199 & ipompc) 200 cycle 201 endif 202 endif 203 write(*,*) 204 & '*WARNING reading *EQUATION: MPC to remove' 205 write(*,*) ' is not defined; node:',k 206 write(*,*) ' degree of freedom:',i1 207 enddo 208 else 209 k=ialset(j-2) 210 do 211 k=k-ialset(j) 212 if(k.ge.ialset(j-1)) exit 213 do i1=impcstart,impcend 214 idof=8*(k-1)+i1 215 call nident(ikmpc,idof,nmpc,id) 216 if(id.gt.0) then 217 if(ikmpc(id).eq.idof) then 218 impc=ilmpc(id) 219 call mpcrem(impc,mpcfree, 220 & nodempc,nmpc,ikmpc,ilmpc,labmpc, 221 & coefmpc,ipompc) 222 cycle 223 endif 224 endif 225 write(*,*) 226 & '*WARNING reading *EQUATION: MPC to remove' 227 write(*,*) 228 & ' is not defined; node:',k 229 write(*,*)' degree of freedom:',i1 230 enddo 231 enddo 232 endif 233 enddo 234 endif 235 enddo 236 return 237 else 238 write(*,*) 239 & '*WARNING reading *EQUATION: parameter not recognized:' 240 write(*,*) ' ', 241 & textpart(m)(1:index(textpart(m),' ')-1) 242 call inputwarning(inpc,ipoinpc,iline, 243 & "*EQUATION%") 244 endif 245 enddo 246! 247 if(istep.gt.0) then 248 write(*,*) 249 & '*ERROR reading *EQUATION: *EQUATION should be placed' 250 write(*,*) ' before all step definitions' 251 ier=1 252 return 253 endif 254! 255 do 256 call getnewline(inpc,textpart,istat,n,key,iline,ipol,inl, 257 & ipoinp,inp,ipoinpc) 258 if((istat.lt.0).or.(key.eq.1)) return 259 read(textpart(1)(1:10),'(i10)',iostat=istat) nterm 260! 261 nmpc=nmpc+1 262 if(nmpc.gt.nmpc_) then 263 write(*,*) '*ERROR reading *EQUATION: increase nmpc_' 264 ier=1 265 return 266 endif 267! 268 labmpc(nmpc)=' ' 269 ipompc(nmpc)=mpcfree 270 ii=0 271! 272 do 273 call getnewline(inpc,textpart,istat,n,key,iline,ipol,inl, 274 & ipoinp,inp,ipoinpc) 275 if((istat.lt.0).or.(key.eq.1)) then 276 write(*,*) '*ERROR reading *EQUATION: mpc definition ', 277 & nmpc 278 write(*,*) ' is not complete. ' 279 call inputerror(inpc,ipoinpc,iline, 280 & "*EQUATION%",ier) 281 return 282 endif 283! 284 do i=1,n/3 285! 286 read(textpart((i-1)*3+1)(1:10),'(i10)',iostat=istat) node 287 if(istat.gt.0) then 288 call inputerror(inpc,ipoinpc,iline, 289 & "*EQUATION%",ier) 290 return 291 endif 292 if((node.gt.nk).or.(node.le.0)) then 293 write(*,*) '*ERROR reading *EQUATION:' 294 write(*,*) ' node ',node,' is not defined' 295 ier=1 296 return 297 endif 298! 299 read(textpart((i-1)*3+2)(1:10),'(i10)',iostat=istat) ndir 300 if(istat.gt.0) then 301 call inputerror(inpc,ipoinpc,iline, 302 & "*EQUATION%",ier) 303 return 304 endif 305 if(ndir.le.6) then 306c elseif(ndir.eq.4) then 307c ndir=5 308c elseif(ndir.eq.5) then 309c ndir=6 310c elseif(ndir.eq.6) then 311c ndir=7 312 elseif(ndir.eq.8) then 313 ndir=4 314 elseif(ndir.eq.11) then 315 ndir=0 316 else 317 write(*,*) '*ERROR reading *EQUATION:' 318 write(*,*) ' direction',ndir,' is not defined' 319 ier=1 320 return 321 endif 322! 323 read(textpart((i-1)*3+3)(1:20),'(f20.0)',iostat=istat) x 324 if(istat.gt.0) then 325 call inputerror(inpc,ipoinpc,iline, 326 & "*EQUATION%",ier) 327 return 328 endif 329! 330! check whether the node is transformed 331! 332 if(ntrans.le.0) then 333 itr=0 334 elseif(inotr(1,node).eq.0) then 335 itr=0 336 else 337 itr=inotr(1,node) 338 endif 339! 340 if((itr.eq.0).or.(ndir.eq.0).or.(ndir.eq.4)) then 341 nodempc(1,mpcfree)=node 342 nodempc(2,mpcfree)=ndir 343 coefmpc(mpcfree)=x 344! 345! updating ikmpc and ilmpc 346! 347 if(ii.eq.0) then 348 idof=8*(node-1)+ndir 349 call nident(ikmpc,idof,nmpc-1,id) 350 if(id.gt.0) then 351 if(ikmpc(id).eq.idof) then 352 write(*,100) 353 & (ikmpc(id))/8+1,ikmpc(id)-8*((ikmpc(id))/8) 354 ier=1 355 return 356 endif 357 endif 358 do j=nmpc,id+2,-1 359 ikmpc(j)=ikmpc(j-1) 360 ilmpc(j)=ilmpc(j-1) 361 enddo 362 ikmpc(id+1)=idof 363 ilmpc(id+1)=nmpc 364 endif 365! 366 mpcfreeold=mpcfree 367 mpcfree=nodempc(3,mpcfree) 368 if(mpcfree.eq.0) then 369 write(*,*) 370 & '*ERROR reading *EQUATION: increase memmpc_' 371 ier=1 372 return 373 endif 374 else 375 call transformatrix(trab(1,inotr(1,node)), 376 & co(1,node),a) 377! 378 number=ndir-1 379 if(ii.eq.0) then 380! 381! determining which direction to use for the 382! dependent side: should not occur on the dependent 383! side in another MPC and should have a nonzero 384! coefficient 385! 386 do j=1,3 387 number=number+1 388 if(number.gt.3) number=1 389 idof=8*(node-1)+number 390 call nident(ikmpc,idof,nmpc-1,id) 391 if(id.gt.0) then 392 if(ikmpc(id).eq.idof) then 393 cycle 394 endif 395 endif 396 if(dabs(a(number,ndir)).lt.1.d-5) cycle 397 exit 398 enddo 399 if(j.gt.3) then 400 write(*,*) 401 & '*ERROR reading *EQUATION: SPC in node' 402 write(*,*) node,' in transformed coordinates' 403 write(*,*) ' cannot be converted in MPC: all' 404 write(*,*) ' DOFs in the node are used as' 405 write(*,*) ' dependent nodes in other MPCs' 406 ier=1 407 return 408 endif 409 number=number-1 410! 411! updating ikmpc and ilmpc 412! 413 do j=nmpc,id+2,-1 414 ikmpc(j)=ikmpc(j-1) 415 ilmpc(j)=ilmpc(j-1) 416 enddo 417 ikmpc(id+1)=idof 418 ilmpc(id+1)=nmpc 419 endif 420! 421 do j=1,3 422 number=number+1 423 if(number.gt.3) number=1 424 if(dabs(a(number,ndir)).lt.1.d-5) cycle 425 nodempc(1,mpcfree)=node 426 nodempc(2,mpcfree)=number 427 coefmpc(mpcfree)=x*a(number,ndir) 428 mpcfreeold=mpcfree 429 mpcfree=nodempc(3,mpcfree) 430 if(mpcfree.eq.0) then 431 write(*,*) 432 & '*ERROR reading *EQUATION: increase memmpc_' 433 ier=1 434 return 435 endif 436 enddo 437 endif 438! 439 ii=ii+1 440 enddo 441! 442 if(ii.eq.nterm) then 443 nodempc(3,mpcfreeold)=0 444 exit 445 endif 446 enddo 447 enddo 448! 449 100 format(/,'*ERROR reading *EQUATION: the DOF corresponding to', 450 & /,'node ',i10,' in direction',i1,' is detected on', 451 & /,'the dependent side of two different MPC''s') 452 return 453 end 454 455