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 distributingcouplings(inpc,textpart,ipompc,nodempc, 20 & coefmpc,nmpc,nmpc_,mpcfree,nk,ikmpc,ilmpc, 21 & labmpc,istep,istat,n,iline,ipol,inl,ipoinp,inp,ipoinpc,lakon, 22 & kon,ipkon,set,nset,istartset,iendset,ialset,co,ier) 23! 24! reading the input deck: *DISTRIBUTING COUPLING 25! 26 implicit none 27! 28 character*1 inpc(*) 29 character*8 lakon(*) 30 character*20 labmpc(*) 31 character*81 set(*),elset,noset 32 character*132 textpart(16) 33! 34 integer ipompc(*),nodempc(3,*),nmpc,nmpc_,mpcfree,istep,istat, 35 & n,i,j,key,nk,node,ier, 36 & mpcfreeold,ikmpc(*),ilmpc(*),id,idof,iline,ipol,inl, 37 & ipoinp(2,*),inp(3,*),ipoinpc(0:*),irefnode, 38 & k,ipos,kon(*),ipkon(*),nset,idir,newmpc, 39 & istartset(*),iendset(*),ialset(*),indexm,ielem 40! 41 real*8 coefmpc(*),co(3,*),weight,totweight 42! 43 elset(1:1)=' ' 44 do i=2,n 45 if(textpart(i)(1:6).eq.'ELSET=') then 46 elset=textpart(i)(7:86) 47 elset(81:81)=' ' 48 ipos=index(elset,' ') 49 elset(ipos:ipos)='E' 50 else 51 write(*,*) '*WARNING reading *DISTRIBUTING COUPLING:' 52 write(*,*) ' parameter not recognized:' 53 write(*,*) ' ', 54 & textpart(i)(1:index(textpart(i),' ')-1) 55 call inputwarning(inpc,ipoinpc,iline, 56 &"*DISTRIBUTING COUPLING%") 57 endif 58 enddo 59! 60 if(elset(1:1).eq.' ') then 61 write(*,*) '*ERROR reading *DISTRIBUTING COUPLING:' 62 write(*,*) ' no element set given' 63 call inputerror(inpc,ipoinpc,iline, 64 & "*DISTRIBUTING COUPLING%",ier) 65 return 66 endif 67! 68! check whether the element set exists 69! 70c do i=1,nset 71c if(set(i).eq.elset) exit 72c enddo 73 call cident81(set,elset,nset,id) 74 i=nset+1 75 if(id.gt.0) then 76 if(elset.eq.set(id)) then 77 i=id 78 endif 79 endif 80 if(i.gt.nset) then 81 write(*,*) '*ERROR reading *DISTRIBUTING COUPLING:' 82 write(*,*) ' element set ',elset(1:ipos-1), 83 & ' is not defined' 84 ier=1 85 return 86 endif 87! 88! check whether only one element belongs 89! to the set 90! 91 if(istartset(i).ne.iendset(i)) then 92 write(*,*) '*ERROR reading *DISTRIBUTING COUPLING:' 93 write(*,*) ' element set ',elset(1:ipos-1), 94 & ' contains more than one element' 95 ier=1 96 return 97 endif 98! 99! check whether the element is a DCOUP3D element 100! 101 ielem=ialset(istartset(i)) 102 if(lakon(ielem)(1:7).ne.'DCOUP3D') then 103 write(*,*) '*ERROR reading *DISTRIBUTING COUPLING:' 104 write(*,*) ' element ',ielem,' is not a' 105 write(*,*) ' DCOUP3D element' 106 ier=1 107 return 108 endif 109! 110! the reference node belongs to the DCOUP3D element 111! 112 irefnode=kon(ipkon(ielem)+1) 113 newmpc=0 114 totweight=0.d0 115! 116! generate a MPC for dof 1 117! 118 do 119 call getnewline(inpc,textpart,istat,n,key,iline,ipol,inl, 120 & ipoinp,inp,ipoinpc) 121 if((istat.lt.0).or.(key.eq.1)) exit 122! 123 read(textpart(1)(1:10),'(i10)',iostat=istat) node 124 if(istat.eq.0) then 125 if(node.gt.nk) then 126 write(*,*) '*ERROR reading *DISTRIBUTING COUPLING:' 127 write(*,*) ' node ',node,' is not defined' 128 ier=1 129 return 130 endif 131! 132! if first node : new MPC 133! 134 if(newmpc.eq.0) then 135 idof=8*(node-1)+1 136 call nident(ikmpc,idof,nmpc,id) 137 if(id.gt.0) then 138 if(ikmpc(id).eq.idof) then 139 write(*,*) '*ERROR reading *DISTRIBUTING COUPLING:' 140 write(*,*) ' dof 1 of node ',node, 141 & ' is already used' 142 ier=1 143 return 144 endif 145 endif 146! 147 nmpc=nmpc+1 148 if(nmpc.gt.nmpc_) then 149 write(*,*) '*ERROR reading *DISTRIBUTING COUPLING:' 150 write(*,*) ' increase nmpc_' 151 ier=1 152 return 153 endif 154 ipompc(nmpc)=mpcfree 155 labmpc(nmpc)=' ' 156 ipompc(nmpc)=mpcfree 157! 158! updating ikmpc and ilmpc 159! 160 do j=nmpc,id+2,-1 161 ikmpc(j)=ikmpc(j-1) 162 ilmpc(j)=ilmpc(j-1) 163 enddo 164 ikmpc(id+1)=idof 165 ilmpc(id+1)=nmpc 166! 167 newmpc=1 168 endif 169! 170! reading the weight 171! 172 read(textpart(2)(1:20),'(f20.0)',iostat=istat) weight 173 if(istat.gt.0) then 174 call inputerror(inpc,ipoinpc,iline, 175 & "*DISTRIBUTING COUPLING%",ier) 176 return 177 endif 178 totweight=totweight+weight 179! 180! new term in MPC 181! 182 nodempc(1,mpcfree)=node 183 nodempc(2,mpcfree)=1 184 coefmpc(mpcfree)=weight 185 mpcfree=nodempc(3,mpcfree) 186! 187 else 188! 189! node set 190! 191 read(textpart(1)(1:80),'(a80)',iostat=istat) noset 192 read(textpart(2)(1:20),'(f20.0)',iostat=istat) weight 193 if(istat.gt.0) then 194 call inputerror(inpc,ipoinpc,iline, 195 & "*DISTRIBUTING COUPLING%",ier) 196 return 197 endif 198 noset(81:81)=' ' 199 ipos=index(noset,' ') 200 noset(ipos:ipos)='N' 201c do i=1,nset 202c if(set(i).eq.noset) exit 203c enddo 204 call cident81(set,noset,nset,id) 205 i=nset+1 206 if(id.gt.0) then 207 if(noset.eq.set(id)) then 208 i=id 209 endif 210 endif 211 if(i.gt.nset) then 212 noset(ipos:ipos)=' ' 213 write(*,*) '*ERROR reading *DISTRIBUTING COUPLING:' 214 write(*,*) ' node set ',noset 215 write(*,*) ' has not yet been defined. ' 216 call inputerror(inpc,ipoinpc,iline, 217 & "*DISTRIBUTING COUPLING%",ier) 218 return 219 endif 220 do j=istartset(i),iendset(i) 221 if(ialset(j).gt.0) then 222 node=ialset(j) 223 totweight=totweight+weight 224! 225 if(newmpc.eq.0) then 226 idof=8*(node-1)+1 227 call nident(ikmpc,idof,nmpc,id) 228 if(id.gt.0) then 229 if(ikmpc(id).eq.idof) then 230 write(*,*) 231 & '*ERROR reading *DISTRIBUTING COUPLING:' 232 write(*,*) ' dof 1 of node ',node, 233 & ' is already used' 234 ier=1 235 return 236 endif 237 endif 238! 239 nmpc=nmpc+1 240 if(nmpc.gt.nmpc_) then 241 write(*,*) 242 & '*ERROR reading *DISTRIBUTING COUPLING:' 243 write(*,*) ' increase nmpc_' 244 ier=1 245 return 246 endif 247 ipompc(nmpc)=mpcfree 248 labmpc(nmpc)=' ' 249 ipompc(nmpc)=mpcfree 250! 251! updating ikmpc and ilmpc 252! 253 do k=nmpc,id+2,-1 254 ikmpc(k)=ikmpc(k-1) 255 ilmpc(k)=ilmpc(k-1) 256 enddo 257 ikmpc(id+1)=idof 258 ilmpc(id+1)=nmpc 259! 260 newmpc=1 261 endif 262! 263! new term in MPC 264! 265 nodempc(1,mpcfree)=node 266 nodempc(2,mpcfree)=1 267 coefmpc(mpcfree)=weight 268 mpcfree=nodempc(3,mpcfree) 269! 270 else 271 node=ialset(j-2) 272 do 273 node=node-ialset(j) 274 if(node.ge.ialset(j-1)) exit 275 totweight=totweight+weight 276! 277! new term in MPC 278! 279 nodempc(1,mpcfree)=node 280 nodempc(2,mpcfree)=1 281 coefmpc(mpcfree)=weight 282 mpcfree=nodempc(3,mpcfree) 283 enddo 284 endif 285 enddo 286 endif 287 enddo 288! 289! reference node 290! 291 nodempc(1,mpcfree)=irefnode 292 nodempc(2,mpcfree)=1 293 coefmpc(mpcfree)=-totweight 294 mpcfreeold=mpcfree 295 mpcfree=nodempc(3,mpcfree) 296 nodempc(3,mpcfreeold)=0 297! 298! dofs 2 and 3 299! 300 do idir=2,3 301! 302 indexm=ipompc(nmpc) 303 node=nodempc(1,indexm) 304! 305 idof=8*(node-1)+idir 306 call nident(ikmpc,idof,nmpc,id) 307 if(id.gt.0) then 308 if(ikmpc(id).eq.idof) then 309 write(*,*) '*ERROR reading *DISTRIBUTING COUPLING:' 310 write(*,*) ' dof',idir,' of node ',node, 311 & ' is already used' 312 ier=1 313 return 314 endif 315 endif 316! 317 nmpc=nmpc+1 318 if(nmpc.gt.nmpc_) then 319 write(*,*) '*ERROR reading *DISTRIBUTING COUPLING:' 320 write(*,*) ' increase nmpc_' 321 ier=1 322 return 323 endif 324 ipompc(nmpc)=mpcfree 325 labmpc(nmpc)=' ' 326 ipompc(nmpc)=mpcfree 327! 328! updating ikmpc and ilmpc 329! 330 do j=nmpc,id+2,-1 331 ikmpc(j)=ikmpc(j-1) 332 ilmpc(j)=ilmpc(j-1) 333 enddo 334 ikmpc(id+1)=idof 335 ilmpc(id+1)=nmpc 336! 337 do 338 nodempc(1,mpcfree)=nodempc(1,indexm) 339 nodempc(2,mpcfree)=idir 340 coefmpc(mpcfree)=coefmpc(indexm) 341 if(nodempc(3,indexm).eq.0) then 342 mpcfreeold=mpcfree 343 mpcfree=nodempc(3,mpcfree) 344 nodempc(3,mpcfreeold)=0 345 exit 346 else 347 mpcfree=nodempc(3,mpcfree) 348 indexm=nodempc(3,indexm) 349 endif 350 enddo 351 enddo 352! 353 return 354 end 355 356