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 modifympc(inodestet,nnodestet,co,doubleglob, 20 & integerglob,ipompc,nodempc,coefmpc,nmpc,nmpc_,labmpc,mpcfree, 21 & ikmpc,ilmpc,jq,irow,icol,loc,irowt,jqt,itemp,au,ixcol, 22 & ikboun,nboun,nodeboun,mpcrfna,mpcrfnb,nodempcref, 23 & coefmpcref,memmpcref_,mpcfreeref,maxlenmpcref,memmpc_, 24 & maxlenmpc,istep) 25! 26! generating MPC's connecting the nodes in the old tet mesh in 27! which SPC's or point forces were defined with the nodes in the 28! refined tet mesh 29! 30 implicit none 31! 32 logical spc 33! 34 character*20 labmpc(*) 35! 36 integer inodestet(*),nnodestet,integerglob(*),nktet,netet,ne,nkon, 37 & nfaces,nfield,nselect,imastset,iselect(1),nterms,idir1, 38 & nelem,ialset(1),mpcfreeold,jq(*),irow(*),icol(*), 39 & iendset(1),istartset(1),konl(20),loopa,loc(*),irowt(*), 40 & node,i,j,k,m,ipompc(*),nodempc(3,*),nmpc,nmpc_,mpcfree, 41 & ikmpc(*),ilmpc(*),idof,id,jqt(*),nzs,nosol,nodeboun(*), 42 & kflag,itemp(*),maxrow,ixcol(*),isoltot,newnode,iy,ic,idummy, 43 & nlength,kopt,ixcolnew,idnew,i1,ikboun(*),nboun,mpcrfna, 44 & mpcrfnb,index,index1,node1,mpc,idof1,idof2,id1,id2, 45 & nodempcref(3,*),memmpcref_,mpcfreeref,maxlenmpcref,memmpc_, 46 & maxlenmpc,istep 47! 48 real*8 co(3,*),doubleglob(*),coords(3),ratio(20),value,au(*), 49 & coefmpc(*),depcoef,dist,xopt,coef1,coefmpcref(*) 50! 51 if(istep.gt.1) then 52! 53! restoring the original equations 54! 55 memmpc_=memmpcref_ 56 mpcfree=mpcfreeref 57 maxlenmpc=maxlenmpcref 58! 59! mpcfreeref=-1 indicates that the MPC's have changed and that 60! a new copy is to be taken before the next call to cascade.c 61! 62 mpcfreeref=-1 63! 64 do i=1,memmpc_ 65 do j=1,3 66 nodempc(j,i)=nodempcref(j,i) 67 enddo 68 coefmpc(i)=coefmpcref(i) 69 enddo 70 endif 71! 72! restore the original form of the MPC's 73! an MPC with number in between mpcrfna and mpcrfnb (i.e. MPC's 74! corresponding to the connection of the old mesh with the new 75! mesh) is in its orignal form if the coefficient of the first 76! term is exactly 1.d0 77! 78 do i=mpcrfna,mpcrfnb 79 index=ipompc(i) 80 if(coefmpc(index).eq.1.d0) cycle 81! 82! MPC has been reordered and has to be brought in its original 83! form 84! 85 node1=nodempc(1,index) 86 idir1=nodempc(2,index) 87 coef1=coefmpc(index) 88 idof1=8*(node1-1)+idir1 89 index1=index 90 call nident(ikmpc,idof1,nmpc,id1) 91! 92 index=nodempc(3,index) 93! 94 do 95 if(index.eq.0) then 96 write(*,*) '*ERROR in modifympc: equation ',i 97 write(*,*) ' cannot be modified' 98 call exit(201) 99 endif 100 if(coefmpc(index).eq.1.d0) then 101! 102! restore the leading term of the equation 103! 104 nodempc(1,index1)=nodempc(1,index) 105 coefmpc(index1)=coefmpc(index) 106! 107! restore the independent term of the equation 108! id1>id2 since the refined mesh has higher node numbers 109! than the unrefined mesh 110! 111 idof2=8*(nodempc(1,index)-1)+nodempc(2,index) 112 call nident(ikmpc,idof2,nmpc,id2) 113 do j=id1,id2+2,-1 114 ikmpc(j)=ikmpc(j-1) 115 ilmpc(j)=ilmpc(j-1) 116 enddo 117 ikmpc(id2+1)=idof2 118 ilmpc(id2+1)=i 119! 120 nodempc(1,index)=node1 121 coefmpc(index)=coef1 122 exit 123 else 124 index=nodempc(3,index) 125 endif 126 enddo 127 enddo 128! 129! compile all dependent nodes in jq, irow.. 130! taking the equations for the first dof (x-direction) only 131! 132 nzs=0 133 nnodestet=0 134 maxrow=0 135! 136 do i=mpcrfna,mpcrfnb 137 index=ipompc(i) 138! 139! x-direction 140! 141 if(nodempc(2,index).ne.1) cycle 142! 143 nnodestet=nnodestet+1 144 inodestet(nnodestet)=nodempc(1,index) 145 jq(nnodestet)=nzs+1 146 index=nodempc(3,index) 147 do 148 if(index.eq.0) exit 149 nzs=nzs+1 150 irow(nzs)=nodempc(1,index) 151 maxrow=max(maxrow,irow(nzs)) 152 au(nzs)=coefmpc(index) 153 irowt(nzs)=nnodestet 154 loc(nzs)=nzs 155 index=nodempc(3,index) 156 enddo 157 enddo 158 jq(nnodestet+1)=nzs+1 159! 160! determine the nodes for which the dependent term in the MPC's 161! has to be modified: these are MPC's for which the dependent term 162! is subject to other SPC's or MPC's 163! 164! SPC's 165! 166 do i=1,nboun 167 node=nodeboun(i) 168 call nident(inodestet,node,nnodestet,id) 169 if(id.gt.0) then 170 if(inodestet(id).eq.node) then 171 icol(id)=1 172 endif 173 endif 174 enddo 175! 176! MPC's 177! 178 do i=1,mpcrfna-1 179 index=ipompc(i) 180 node=nodempc(1,index) 181 call nident(inodestet,node,nnodestet,id) 182 if(id.gt.0) then 183 if(inodestet(id).eq.node) then 184 icol(id)=1 185 endif 186 endif 187 enddo 188! 189! for the nodes i for which the equations have to be modified 190! icol(i) is nonzero 191! 192 isoltot=0 193 do i=1,nnodestet 194 if(icol(i).eq.1) then 195 icol(i)=jq(i+1)-jq(i) 196 else 197 isoltot=isoltot+1 198 endif 199 enddo 200! 201! copying irow into itemp, sorting itemp and field loc and 202! irowt along 203! 204 do j=1,nzs 205 itemp(j)=irow(j) 206 enddo 207! 208 kflag=2 209 call isortiii(itemp,loc,irowt,nzs,kflag) 210! 211! determining jqt (columns per row); jqt and irowt is the 212! equivalent (for row per row treatment) of jq and irow (for 213! column per column treatment); one can consider jqt and irowt 214! to be the jq and irow fields for the transpose of the matrix 215! (therefore the appended letter "t") 216! 217 j=1 218 do m=1,nzs 219 if(itemp(m).ge.j) then 220 do k=j,itemp(m) 221 jqt(k)=m 222 enddo 223 j=itemp(m)+1 224 endif 225 enddo 226 jqt(j)=nzs+1 227! 228! creating a unique number consisting of the column number and the 229! number of rows per column; if sorted, the result is the same as 230! for sorting icol. 231! 232 do j=1,nnodestet 233 if(icol(j).gt.0) then 234 ixcol(j)=(nnodestet+1)*icol(j)+j 235 endif 236 enddo 237! 238! sorting ixcol 239! 240 kflag=1 241 call isortii(ixcol,idummy,nnodestet,kflag) 242! 243! reordering the equations: the first term is to be a 244! new node instead of an old node; starting with columns 245! with only 1 nonzero value 246! 247 nosol=0 248 do 249 if(isoltot.eq.nnodestet) exit 250 isoltot=isoltot+1 251! 252! local dependent node number 253! 254 i=ixcol(isoltot) 255 & -int(ixcol(isoltot)/(nnodestet+1))*(nnodestet+1) 256 ixcol(isoltot)=0 257 icol(i)=0 258! 259! global dependent node number 260! 261 node=inodestet(i) 262! 263! determining the optimal row for the exchange 264! 265 kopt=0 266 xopt=0.d0 267 do k=jq(i),jq(i+1)-1 268 if(irow(k).gt.0) then 269 if(dabs(au(k)).gt.xopt) then 270 kopt=k 271 xopt=dabs(au(k)) 272 endif 273 endif 274 enddo 275 k=kopt 276! 277 if(k.gt.0) then 278! 279! row dof 280! 281 newnode=irow(k) 282 depcoef=au(k) 283! 284! generating MPC's 285! 286 do j=1,3 287! 288! modify ikmpc to take the switch from id2 to id1 as 289! dependent node into account 290! 291 idof1=8*(node-1)+j 292 call nident(ikmpc,idof1,nmpc,id1) 293 do 294 if(labmpc(ilmpc(id1)).eq.'RM ') exit 295 id=id-1 296 enddo 297 mpc=ilmpc(id1) 298 index=ipompc(mpc) 299 idof2=8*(newnode-1)+j 300 call nident(ikmpc,idof2,nmpc,id2) 301 do m=id1,id2-1 302 ikmpc(m)=ikmpc(m+1) 303 ilmpc(m)=ilmpc(m+1) 304 enddo 305 ikmpc(id2)=idof2 306 ilmpc(id2)=mpc 307! 308! inserting the new dependent term 309! 310 nodempc(1,index)=newnode 311 coefmpc(index)=depcoef 312! 313! inserting the new independent term 314! 315 do 316 index=nodempc(3,index) 317 if(index.eq.0) then 318 write(*,*) '*ERROR in modifympc: equation ',mpc 319 write(*,*) ' cannot be modified' 320 call exit(201) 321 endif 322 if(nodempc(1,index).eq.newnode) then 323 nodempc(1,index)=node 324 coefmpc(index)=1.d0 325 exit 326 endif 327 enddo 328 enddo 329! 330! remove the indpendent node from the data base: 331! tagging the row dof in all other columns 332! 333 do m=jqt(newnode),jqt(newnode+1)-1 334 if(irowt(m).eq.i) cycle 335 ic=irowt(m) 336! 337! no SPC/MPC or already treated 338! 339 if(icol(ic).eq.0) cycle 340! 341! setting negative sign in irow to show that the row 342! is not available any more (already used as dependent dof) 343! 344 irow(loc(m))=-irow(loc(m)) 345 iy=(nnodestet+1)*icol(ic)+ic 346 nlength=nnodestet-isoltot 347 call nident(ixcol(isoltot+1),iy,nlength,id) 348 if(ixcol(isoltot+id).ne.iy) then 349 write(*,*) '*ERROR in modifympc:' 350 write(*,*) ' data base corrupt' 351 call exit(201) 352 endif 353 icol(ic)=icol(ic)-1 354 ixcolnew=ixcol(isoltot+id)-(nnodestet+1) 355! 356! not treated equation has no independent terms any more 357! 358 if(ixcolnew.lt.(nnodestet+1)) then 359 nosol=nosol+1 360 if(nosol.eq.1) then 361 write(*,*) '*WARNING in modifympc; failed to connect' 362 write(*,*) ' node ',inodestet(ic), 363 & ' of unrefined mesh' 364 write(*,*) 365 &' to the refined mesh; loads and boundary conditions' 366 write(*,*) 367 &' in this node are not taken into account; a list of' 368 write(*,*) 369 &' not connected nodes is stored in' 370 write(*,*) 371 &' WarnNodeMissRefineConnection.nam' 372 else 373 write(*,*) '*WARNING in modifympc; failed to connect' 374 write(*,*) ' node ',inodestet(ic), 375 & ' of unrefined mesh' 376 endif 377 write(23,*) 378 & '*NSET,NSET=WarnNodeMissRefineConnection' 379 write(23,*) inodestet(ic) 380 endif 381 call nident(ixcol(isoltot+1),ixcolnew,id-1,idnew) 382 do i1=isoltot+id,isoltot+idnew+2,-1 383 ixcol(i1)=ixcol(i1-1) 384 enddo 385 ixcol(isoltot+idnew+1)=ixcolnew 386 enddo 387 endif 388 enddo 389 write(*,*) 390 write(*,*) '*INFO in modifympc: no solution for ',nosol,' nodes.' 391 write(*,*) 392 close(23) 393! 394 return 395 end 396