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 knotmpc(ipompc,nodempc,coefmpc,irefnode,irotnode, 20 & iexpnode, 21 & labmpc,nmpc,nmpc_,mpcfree,ikmpc,ilmpc,nk,nk_,nodeboun,ndirboun, 22 & ikboun,ilboun,nboun,nboun_,idepnodes,typeboun,co,xboun,istep, 23 & k,ndepnodes,idim,e1,e2,t1) 24! 25! generates three knot MPC's for node "node" about reference 26! (translational) node irefnode and rotational node irotnode 27! 28 implicit none 29! 30 character*1 typeboun(*) 31 character*20 labmpc(*) 32! 33 integer ipompc(*),nodempc(3,*),nmpc,nmpc_,mpcfree,nk,nk_, 34 & ikmpc(*),idepnodes(*),k,ndepnodes,idim,n,matz,ier,i, 35 & ilmpc(*),node,id,mpcfreeold,j,idof,l,nodeboun(*), 36 & ndirboun(*),ikboun(*),ilboun(*),nboun,nboun_,irefnode, 37 & irotnode,iexpnode,istep,ispcnode,inode 38! 39 real*8 coefmpc(*),co(3,*),xboun(*),e(3,3,3),dc(3,3,3),s(3,3), 40 & w(3),z(3,3),fv1(3),fv2(3),e1(3),e2(3),t1(3),sx,sy,sz,sxx, 41 & sxy,sxz,syy,syz,szz,u2(3,3),u3(3,3) 42! 43! e_ijk symbol 44! 45 data e /0.d0,0.d0,0.d0,0.d0,0.d0,-1.d0,0.d0,1.d0,0.d0, 46 & 0.d0,0.d0,1.d0,0.d0,0.d0,0.d0,-1.d0,0.d0,0.d0, 47 & 0.d0,-1.d0,0.d0,1.d0,0.d0,0.d0,0.d0,0.d0,0.d0/ 48! 49! dc_ijk=e_ikj 50! 51 data dc /0.d0,0.d0,0.d0,0.d0,0.d0,1.d0,0.d0,-1.d0,0.d0, 52 & 0.d0,0.d0,-1.d0,0.d0,0.d0,0.d0,1.d0,0.d0,0.d0, 53 & 0.d0,1.d0,0.d0,-1.d0,0.d0,0.d0,0.d0,0.d0,0.d0/ 54! 55c idim=1 56! 57! on entry: if idim=1: only 2-d elements 58! if idim=3: at least one 1-d element 59! 60 if((istep.gt.1).and.(k.eq.1)) then 61! 62! determine the dimensionality of the knot 63! (only for step 2 or higher since in step 1 the geometry 64! is not known yet) 65! 66! determine the area moments of inertia 67! 68 sx=0.d0 69 sy=0.d0 70 sz=0.d0 71 sxx=0.d0 72 sxy=0.d0 73 sxz=0.d0 74 syy=0.d0 75 syz=0.d0 76 szz=0.d0 77! 78 do i=1,ndepnodes 79 node=idepnodes(i) 80 sx=sx+co(1,node) 81 sy=sy+co(2,node) 82 sz=sz+co(3,node) 83 sxx=sxx+co(1,node)*co(1,node) 84 sxy=sxy+co(1,node)*co(2,node) 85 sxz=sxz+co(1,node)*co(3,node) 86 syy=syy+co(2,node)*co(2,node) 87 syz=syz+co(2,node)*co(3,node) 88 szz=szz+co(3,node)*co(3,node) 89 enddo 90! 91 sxx=sxx-sx*sx/ndepnodes 92 sxy=sxy-sx*sy/ndepnodes 93 sxz=sxz-sx*sz/ndepnodes 94 syy=syy-sy*sy/ndepnodes 95 syz=syz-sy*sz/ndepnodes 96 szz=szz-sz*sz/ndepnodes 97! 98c write(*,*) 'sxx...',sxx,sxy,sxz,syy,syz,szz 99! 100 s(1,1)=sxx 101 s(1,2)=sxy 102 s(1,3)=sxz 103 s(2,1)=sxy 104 s(2,2)=syy 105 s(2,3)=syz 106 s(3,1)=sxz 107 s(3,2)=syz 108 s(3,3)=szz 109! 110! determining the eigenvalues 111! 112 n=3 113 matz=1 114 ier=0 115 call rs(n,n,s,w,matz,z,fv1,fv2,ier) 116 if(ier.ne.0) then 117 write(*,*) '*ERROR in knotmpc while calculating the' 118 write(*,*) ' eigenvalues/eigenvectors' 119 call exit(201) 120 endif 121! 122! the eigenvalues are the moments of inertia w.r.t. the 123! plane orthogonal to the eigenvector 124! 125! dimension=1 if the two lowest eigenvalues are zero 126! dimension=2 if only the lowest eigenvalue is zero 127! else dimension=3 128! 129c write(*,*) 'eigenvalues ',w(1),w(2),w(3) 130 if((w(1).lt.1.d-10).and.(w(2).lt.1.d-10)) then 131 idim=min(idim,1) 132c idim=1 133 elseif(w(1).lt.1.d-10) then 134 idim=min(idim,2) 135c idim=2 136 else 137 idim=min(idim,1) 138c idim=3 139 endif 140c write(*,*) 'knotmpc irefnode= ',irefnode,' idim= ',idim 141! 142! defining a local coordinate system for idim=2 143! 144 if(idim.eq.2) then 145 do i=1,3 146 t1(i)=z(i,1) 147 e2(i)=z(i,2) 148 e1(i)=z(i,3) 149 enddo 150! 151! check whether e1-e2-t1 is a rhs system 152! 153 if(t1(1)*(e1(2)*e2(3)-e1(3)*e2(2))- 154 & t1(2)*(e1(1)*e2(3)-e1(3)*e2(1))+ 155 & t1(3)*(e1(1)*e2(2)-e1(2)*e2(1)).lt.0.d0) then 156 do i=1,3 157 t1(i)=-t1(i) 158 enddo 159 endif 160c write(*,*) 't1 ',t1(1),t1(2),t1(3) 161c write(*,*) 'e1 ',e1(1),e1(2),e1(3) 162c write(*,*) 'e2 ',e2(1),e2(2),e2(3) 163! 164! storing t1 and e1 as coordinates of irotnode and 165! iexpnode, respectively 166! 167 do i=1,3 168 co(i,irotnode)=t1(i) 169 co(i,iexpnode)=e1(i) 170 enddo 171 endif 172 endif 173! 174 inode=idepnodes(k) 175! 176 nk=nk+1 177 if(nk.gt.nk_) then 178 write(*,*) '*ERROR in knotmpc: increase nk_' 179 call exit(201) 180 endif 181! 182 ispcnode=nk 183! 184 if((idim.eq.1).or.(idim.eq.3)) then 185! 186! knot nodes lie on a line 187! 188 do j=1,3 189 idof=8*(inode-1)+j 190 call nident(ikmpc,idof,nmpc,id) 191 if(id.gt.0) then 192 if(ikmpc(id).eq.idof) then 193 cycle 194 endif 195 endif 196 nmpc=nmpc+1 197 if(nmpc.gt.nmpc_) then 198 write(*,*) '*ERROR in knotmpc: increase nmpc_' 199 call exit(201) 200 endif 201! 202 ipompc(nmpc)=mpcfree 203 labmpc(nmpc)='KNOT ' 204 write(labmpc(nmpc)(5:5),'(i1)') idim 205! 206 do l=nmpc,id+2,-1 207 ikmpc(l)=ikmpc(l-1) 208 ilmpc(l)=ilmpc(l-1) 209 enddo 210 ikmpc(id+1)=idof 211 ilmpc(id+1)=nmpc 212! 213 nodempc(1,mpcfree)=inode 214 nodempc(2,mpcfree)=j 215 coefmpc(mpcfree)=1.d0 216 mpcfree=nodempc(3,mpcfree) 217! 218! translation term 219! 220 nodempc(1,mpcfree)=irefnode 221 nodempc(2,mpcfree)=j 222 coefmpc(mpcfree)=-1.d0 223 mpcfree=nodempc(3,mpcfree) 224! 225! expansion term 226! 227 nodempc(1,mpcfree)=iexpnode 228 nodempc(2,mpcfree)=1 229 if(istep.gt.1) then 230 coefmpc(mpcfree)=co(j,irefnode)-co(j,inode) 231 endif 232 mpcfree=nodempc(3,mpcfree) 233! 234! rotation terms 235! 236 nodempc(1,mpcfree)=irotnode 237 nodempc(2,mpcfree)=1 238 if(istep.gt.1) then 239 coefmpc(mpcfree)=dc(j,1,1)*(co(1,irefnode)-co(1,inode))+ 240 & dc(j,2,1)*(co(2,irefnode)-co(2,inode))+ 241 & dc(j,3,1)*(co(3,irefnode)-co(3,inode)) 242 endif 243 mpcfree=nodempc(3,mpcfree) 244 nodempc(1,mpcfree)=irotnode 245 nodempc(2,mpcfree)=2 246 if(istep.gt.1) then 247 coefmpc(mpcfree)=dc(j,1,2)*(co(1,irefnode)-co(1,inode))+ 248 & dc(j,2,2)*(co(2,irefnode)-co(2,inode))+ 249 & dc(j,3,2)*(co(3,irefnode)-co(3,inode)) 250 endif 251 mpcfree=nodempc(3,mpcfree) 252 nodempc(1,mpcfree)=irotnode 253 nodempc(2,mpcfree)=3 254 if(istep.gt.1) then 255 coefmpc(mpcfree)=dc(j,1,3)*(co(1,irefnode)-co(1,inode))+ 256 & dc(j,2,3)*(co(2,irefnode)-co(2,inode))+ 257 & dc(j,3,3)*(co(3,irefnode)-co(3,inode)) 258 endif 259 mpcfree=nodempc(3,mpcfree) 260 nodempc(1,mpcfree)=ispcnode 261 nodempc(2,mpcfree)=j 262 coefmpc(mpcfree)=1.d0 263 mpcfreeold=mpcfree 264 mpcfree=nodempc(3,mpcfree) 265 nodempc(3,mpcfreeold)=0 266 idof=8*(ispcnode-1)+j 267 call nident(ikboun,idof,nboun,id) 268 nboun=nboun+1 269 if(nboun.gt.nboun_) then 270 write(*,*) '*ERROR in knotmpc: increase nboun_' 271 call exit(201) 272 endif 273 nodeboun(nboun)=ispcnode 274 ndirboun(nboun)=j 275 typeboun(nboun)='R' 276 if(istep.gt.1) then 277 xboun(nboun)=0.d0 278 endif 279 do l=nboun,id+2,-1 280 ikboun(l)=ikboun(l-1) 281 ilboun(l)=ilboun(l-1) 282 enddo 283 ikboun(id+1)=idof 284 ilboun(id+1)=nboun 285 enddo 286 elseif(idim.eq.2) then 287! 288! knot nodes lie in a plane 289! 290 do i=1,3 291 do j=1,3 292 u2(i,j)=2.d0*e1(i)*e1(j) 293 u3(i,j)=2.d0*e2(i)*e2(j) 294 enddo 295 enddo 296! 297 do j=1,3 298 idof=8*(inode-1)+j 299 call nident(ikmpc,idof,nmpc,id) 300 if(id.gt.0) then 301 if(ikmpc(id).eq.idof) then 302 cycle 303 endif 304 endif 305 nmpc=nmpc+1 306 if(nmpc.gt.nmpc_) then 307 write(*,*) '*ERROR in knotmpc: increase nmpc_' 308 call exit(201) 309 endif 310! 311 ipompc(nmpc)=mpcfree 312 labmpc(nmpc)='KNOT2 ' 313! 314 do l=nmpc,id+2,-1 315 ikmpc(l)=ikmpc(l-1) 316 ilmpc(l)=ilmpc(l-1) 317 enddo 318 ikmpc(id+1)=idof 319 ilmpc(id+1)=nmpc 320! 321 nodempc(1,mpcfree)=inode 322 nodempc(2,mpcfree)=j 323 coefmpc(mpcfree)=1.d0 324 mpcfree=nodempc(3,mpcfree) 325! 326! translation term 327! 328 nodempc(1,mpcfree)=irefnode 329 nodempc(2,mpcfree)=j 330 coefmpc(mpcfree)=-1.d0 331 mpcfree=nodempc(3,mpcfree) 332! 333! expansion terms 334! 335! first term is "removed" (= amalgated with the second term) 336! since u1 is the null matrix 337! 338 nodempc(1,mpcfree)=iexpnode 339 nodempc(2,mpcfree)=2 340 coefmpc(mpcfree)=0.d0 341 mpcfree=nodempc(3,mpcfree) 342! 343 nodempc(1,mpcfree)=iexpnode 344 nodempc(2,mpcfree)=2 345 coefmpc(mpcfree)=u2(j,1)*(co(1,irefnode)-co(1,inode))+ 346 & u2(j,2)*(co(2,irefnode)-co(2,inode))+ 347 & u2(j,3)*(co(3,irefnode)-co(3,inode)) 348 mpcfree=nodempc(3,mpcfree) 349! 350 nodempc(1,mpcfree)=iexpnode 351 nodempc(2,mpcfree)=3 352 coefmpc(mpcfree)=u3(j,1)*(co(1,irefnode)-co(1,inode))+ 353 & u3(j,2)*(co(2,irefnode)-co(2,inode))+ 354 & u3(j,3)*(co(3,irefnode)-co(3,inode)) 355 mpcfree=nodempc(3,mpcfree) 356! 357! rotation terms 358! 359 nodempc(1,mpcfree)=irotnode 360 nodempc(2,mpcfree)=1 361 if(istep.gt.1) then 362 coefmpc(mpcfree)=dc(j,1,1)*(co(1,irefnode)-co(1,inode))+ 363 & dc(j,2,1)*(co(2,irefnode)-co(2,inode))+ 364 & dc(j,3,1)*(co(3,irefnode)-co(3,inode)) 365 endif 366 mpcfree=nodempc(3,mpcfree) 367 nodempc(1,mpcfree)=irotnode 368 nodempc(2,mpcfree)=2 369 if(istep.gt.1) then 370 coefmpc(mpcfree)=dc(j,1,2)*(co(1,irefnode)-co(1,inode))+ 371 & dc(j,2,2)*(co(2,irefnode)-co(2,inode))+ 372 & dc(j,3,2)*(co(3,irefnode)-co(3,inode)) 373 endif 374 mpcfree=nodempc(3,mpcfree) 375 nodempc(1,mpcfree)=irotnode 376 nodempc(2,mpcfree)=3 377 if(istep.gt.1) then 378 coefmpc(mpcfree)=dc(j,1,3)*(co(1,irefnode)-co(1,inode))+ 379 & dc(j,2,3)*(co(2,irefnode)-co(2,inode))+ 380 & dc(j,3,3)*(co(3,irefnode)-co(3,inode)) 381 endif 382 mpcfree=nodempc(3,mpcfree) 383 nodempc(1,mpcfree)=ispcnode 384 nodempc(2,mpcfree)=j 385 coefmpc(mpcfree)=1.d0 386 mpcfreeold=mpcfree 387 mpcfree=nodempc(3,mpcfree) 388 nodempc(3,mpcfreeold)=0 389 idof=8*(ispcnode-1)+j 390 call nident(ikboun,idof,nboun,id) 391 nboun=nboun+1 392 if(nboun.gt.nboun_) then 393 write(*,*) '*ERROR in knotmpc: increase nboun_' 394 call exit(201) 395 endif 396 nodeboun(nboun)=ispcnode 397 ndirboun(nboun)=j 398 typeboun(nboun)='R' 399 if(istep.gt.1) then 400 xboun(nboun)=0.d0 401 endif 402 do l=nboun,id+2,-1 403 ikboun(l)=ikboun(l-1) 404 ilboun(l)=ilboun(l-1) 405 enddo 406 ikboun(id+1)=idof 407 ilboun(id+1)=nboun 408 enddo 409 else 410! 411! to do 412! 413 endif 414! 415 return 416 end 417 418 419