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 fillknotmpc(co,ipompc,nodempc,coefmpc,labmpc, 20 & nmpc,nmpcold,mpcfree,idim,e1,e2,t1) 21! 22! updates the coefficients in nonlinear MPC's 23! 24 implicit none 25! 26 character*20 labmpc(*) 27! 28 integer ipompc(*),nodempc(3,*),irefnode,irotnode,idir,idim,n, 29 & nmpc,index,ii,inode,nmpcold,iexpnode,irefnodeprev,i,ndepnodes, 30 & matz,ier,j,indexnext,node,mpcfree,nodeprev 31! 32 real*8 co(3,*),coefmpc(*),e(3,3,3),dc(3,3,3) ,sx,sy,sz,sxx, 33 & sxy,sxz,syy,syz,szz,s(3,3),w(3),z(3,3),fv1(3),fv2(3),e1(3), 34 & e2(3),t1(3),u2(3,3),u3(3,3) 35! 36! e_ijk symbol 37! 38 data e /0.d0,0.d0,0.d0,0.d0,0.d0,-1.d0,0.d0,1.d0,0.d0, 39 & 0.d0,0.d0,1.d0,0.d0,0.d0,0.d0,-1.d0,0.d0,0.d0, 40 & 0.d0,-1.d0,0.d0,1.d0,0.d0,0.d0,0.d0,0.d0,0.d0/ 41! 42! dc_ijk=e_ikj 43! 44 data dc /0.d0,0.d0,0.d0,0.d0,0.d0,1.d0,0.d0,-1.d0,0.d0, 45 & 0.d0,0.d0,-1.d0,0.d0,0.d0,0.d0,1.d0,0.d0,0.d0, 46 & 0.d0,1.d0,0.d0,-1.d0,0.d0,0.d0,0.d0,0.d0,0.d0/ 47! 48 irefnodeprev=0 49! 50 do ii=nmpcold+1,nmpc 51 if(labmpc(ii)(1:4).eq.'KNOT') then 52! 53! from labmpc: if idim=1: only 2-d elements 54! if idim=3: at least one 1-d element 55! 56 irefnode=nodempc(1,nodempc(3,ipompc(ii))) 57! 58 if(irefnode.ne.irefnodeprev) then 59! 60! new knot 61! 62 irefnodeprev=irefnode 63 read(labmpc(ii)(5:5),'(i1)') idim 64! 65! determine the area moments of inertia 66! 67 sx=0.d0 68 sy=0.d0 69 sz=0.d0 70 sxx=0.d0 71 sxy=0.d0 72 sxz=0.d0 73 syy=0.d0 74 syz=0.d0 75 szz=0.d0 76! 77 ndepnodes=0 78! 79 nodeprev=0 80 do i=ii,nmpc 81 if(labmpc(i)(1:4).eq.'KNOT') then 82 if(nodempc(1,nodempc(3,ipompc(i))).eq.irefnode)then 83! 84! node belonging to the same knot 85! 86 node=nodempc(1,ipompc(i)) 87! 88 if(node.ne.nodeprev) then 89 nodeprev=node 90 ndepnodes=ndepnodes+1 91! 92 sx=sx+co(1,node) 93 sy=sy+co(2,node) 94 sz=sz+co(3,node) 95 sxx=sxx+co(1,node)*co(1,node) 96 sxy=sxy+co(1,node)*co(2,node) 97 sxz=sxz+co(1,node)*co(3,node) 98 syy=syy+co(2,node)*co(2,node) 99 syz=syz+co(2,node)*co(3,node) 100 szz=szz+co(3,node)*co(3,node) 101 endif 102 else 103 exit 104 endif 105 else 106 exit 107 endif 108 enddo 109! 110 sxx=sxx-sx*sx/ndepnodes 111 sxy=sxy-sx*sy/ndepnodes 112 sxz=sxz-sx*sz/ndepnodes 113 syy=syy-sy*sy/ndepnodes 114 syz=syz-sy*sz/ndepnodes 115 szz=szz-sz*sz/ndepnodes 116! 117 s(1,1)=sxx 118 s(1,2)=sxy 119 s(1,3)=sxz 120 s(2,1)=sxy 121 s(2,2)=syy 122 s(2,3)=syz 123 s(3,1)=sxz 124 s(3,2)=syz 125 s(3,3)=szz 126! 127! determining the eigenvalues 128! 129 n=3 130 matz=1 131 ier=0 132 call rs(n,n,s,w,matz,z,fv1,fv2,ier) 133 if(ier.ne.0) then 134 write(*,*) '*ERROR in knotmpc while calculating the' 135 write(*,*) ' eigenvalues/eigenvectors' 136 call exit(201) 137 endif 138! 139! the eigenvalues are the moments of inertia w.r.t. the 140! plane orthogonal to the eigenvector 141! 142! dimension=1 if the two lowest eigenvalues are zero 143! dimension=2 if only the lowest eigenvalue is zero 144! else dimension=3 145! 146! the dimension cannot exceed the maximum dimension of 147! the elements connected to the knot (2d-element nodes: 148! dimension 1, 1d-element nodes: dimension 3) 149! 150c write(*,*) 'fillknotmpc eigenvalues ',w(1),w(2),w(3) 151 if((w(1).lt.1.d-10).and.(w(2).lt.1.d-10)) then 152 idim=min(idim,1) 153c idim=1 154 elseif(w(1).lt.1.d-10) then 155 idim=min(idim,2) 156c idim=2 157 else 158 idim=min(idim,1) 159c idim=3 160 endif 161c write(*,*) 'fillknotmpc iref= ',irefnode,' idim= ',idim 162! 163! defining a local coordinate system for idim=2 164! 165 if(idim.eq.2) then 166 do i=1,3 167 t1(i)=z(i,1) 168 e2(i)=z(i,2) 169 e1(i)=z(i,3) 170 enddo 171! 172! check whether e1-e2-t1 is a rhs system 173! 174 if(t1(1)*(e1(2)*e2(3)-e1(3)*e2(2))- 175 & t1(2)*(e1(1)*e2(3)-e1(3)*e2(1))+ 176 & t1(3)*(e1(1)*e2(2)-e1(2)*e2(1)).lt.0.d0) then 177 do i=1,3 178 t1(i)=-t1(i) 179 enddo 180 endif 181c write(*,*) 't1 ',t1(1),t1(2),t1(3) 182c write(*,*) 'e1 ',e1(1),e1(2),e1(3) 183c write(*,*) 'e2 ',e2(1),e2(2),e2(3) 184! 185! storing t1 and e1 as coordinates of irotnode and 186! iexpnode, respectively 187! 188 iexpnode=nodempc(1,nodempc(3,nodempc(3,ipompc(ii)))) 189 irotnode= 190 & nodempc(1,nodempc(3,nodempc(3,nodempc(3,ipompc(ii))))) 191 do i=1,3 192 co(i,irotnode)=t1(i) 193 co(i,iexpnode)=e1(i) 194 enddo 195 endif 196 endif 197! 198 if((idim.eq.1).or.(idim.eq.3)) then 199! 200! knot on a line 201! 202 labmpc(ii)(5:5)='1' 203! 204! dependent node 205! 206 index=ipompc(ii) 207 inode=nodempc(1,index) 208 idir=nodempc(2,index) 209! 210! translation node 211! 212 index=nodempc(3,index) 213 irefnode=nodempc(1,index) 214! 215! expansion node 216! 217 index=nodempc(3,index) 218 iexpnode=nodempc(1,index) 219 coefmpc(index)=co(idir,irefnode)-co(idir,inode) 220! 221! rotation node 222! 223 index=nodempc(3,index) 224 irotnode=nodempc(1,index) 225! 226! determining the coefficients of the rotational degrees 227! of freedom 228! 229 coefmpc(index)=dc(idir,1,1)*(co(1,irefnode)-co(1,inode))+ 230 & dc(idir,2,1)*(co(2,irefnode)-co(2,inode))+ 231 & dc(idir,3,1)*(co(3,irefnode)-co(3,inode)) 232! 233 index=nodempc(3,index) 234 coefmpc(index)=dc(idir,1,2)*(co(1,irefnode)-co(1,inode))+ 235 & dc(idir,2,2)*(co(2,irefnode)-co(2,inode))+ 236 & dc(idir,3,2)*(co(3,irefnode)-co(3,inode)) 237! 238 index=nodempc(3,index) 239 coefmpc(index)=dc(idir,1,3)*(co(1,irefnode)-co(1,inode))+ 240 & dc(idir,2,3)*(co(2,irefnode)-co(2,inode))+ 241 & dc(idir,3,3)*(co(3,irefnode)-co(3,inode)) 242! 243 elseif(idim.eq.2) then 244! 245! nodes of knot lie in a plane 246! 247 labmpc(ii)(5:5)='2' 248! 249 do i=1,3 250 do j=1,3 251 u2(i,j)=2.d0*e1(i)*e1(j) 252 u3(i,j)=2.d0*e2(i)*e2(j) 253 enddo 254 enddo 255! 256! dependent node 257! 258 index=ipompc(ii) 259 inode=nodempc(1,index) 260 idir=nodempc(2,index) 261! 262! translation node 263! 264 index=nodempc(3,index) 265 irefnode=nodempc(1,index) 266! 267! expansion node (first term is amalgated with the second 268! term since the coefficient matrix is zero) 269! 270 index=nodempc(3,index) 271 iexpnode=nodempc(1,index) 272 nodempc(2,index)=2 273 coefmpc(index)=0.d0 274 indexnext=nodempc(3,index) 275 nodempc(3,index)=mpcfree 276! 277 nodempc(1,mpcfree)=iexpnode 278 nodempc(2,mpcfree)=2 279 coefmpc(mpcfree)=u2(idir,1)*(co(1,irefnode)-co(1,inode))+ 280 & u2(idir,2)*(co(2,irefnode)-co(2,inode))+ 281 & u2(idir,3)*(co(3,irefnode)-co(3,inode)) 282 mpcfree=nodempc(3,mpcfree) 283! 284 nodempc(1,mpcfree)=iexpnode 285 nodempc(2,mpcfree)=3 286 coefmpc(mpcfree)=u3(idir,1)*(co(1,irefnode)-co(1,inode))+ 287 & u3(idir,2)*(co(2,irefnode)-co(2,inode))+ 288 & u3(idir,3)*(co(3,irefnode)-co(3,inode)) 289 index=mpcfree 290 mpcfree=nodempc(3,mpcfree) 291 nodempc(3,index)=indexnext 292! 293! rotation node 294! 295 index=indexnext 296 irotnode=nodempc(1,index) 297! 298! determining the coefficients of the rotational degrees 299! of freedom 300! 301 coefmpc(index)=dc(idir,1,1)*(co(1,irefnode)-co(1,inode))+ 302 & dc(idir,2,1)*(co(2,irefnode)-co(2,inode))+ 303 & dc(idir,3,1)*(co(3,irefnode)-co(3,inode)) 304! 305 index=nodempc(3,index) 306 coefmpc(index)=dc(idir,1,2)*(co(1,irefnode)-co(1,inode))+ 307 & dc(idir,2,2)*(co(2,irefnode)-co(2,inode))+ 308 & dc(idir,3,2)*(co(3,irefnode)-co(3,inode)) 309! 310 index=nodempc(3,index) 311 coefmpc(index)=dc(idir,1,3)*(co(1,irefnode)-co(1,inode))+ 312 & dc(idir,2,3)*(co(2,irefnode)-co(2,inode))+ 313 & dc(idir,3,3)*(co(3,irefnode)-co(3,inode)) 314! 315 endif 316 endif 317 enddo 318! 319 return 320 end 321