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 umpc_mean_rot(x,u,f,a,jdof,n,force,iit,idiscon, 20 & jnode,ikmpc,nmpc,ikboun,nboun,label) 21! 22! updates the coefficients in a mean rotation mpc 23! 24! INPUT: 25! 26! x(3,1..n) Carthesian coordinates of the nodes in the 27! user mpc. 28! u(3,1..n) Actual displacements of the nodes in the 29! user mpc. 30! jdof(1..n) Actual degrees of freedom of the mpc terms 31! n number of terms in the user mpc 32! force Actual value of the mpc force 33! iit iteration number 34! 35! OUTPUT: 36! 37! f Actual value of the mpc. If the mpc is 38! exactly satisfied, this value is zero 39! a(n) coefficients of the linearized mpc 40! jdof(1..n) Corrected degrees of freedom of the mpc terms 41! idiscon 0: no discontinuity 42! 1: discontinuity 43! If a discontinuity arises the previous 44! results are not extrapolated at the start of 45! a new increment 46! 47 implicit none 48! 49 character*20 label 50! 51 integer jdof(*),n,nkn,i,j,k,imax,iit,idiscon,node,nodeold, 52 & nodemax,idirold,idirmax,jnode(*),idof,id,iold(3),inew(3), 53 & ikmpc(*),nmpc,ikboun(*),nboun,idir,nendnode 54! 55 real*8 x(3,*),u(3,*),f,a(*),aa(3),cgx(3),cgu(3),pi(3),b(3), 56 & xi(3),dd,al,a1,amax,c1,c2,c3,c4,c9,c10,force,xcoef, 57 & transcoef(3),stdev 58! 59 nkn=(n-1)/3 60 if(3*nkn.ne.n-1) then 61 write(*,*) 62 & '*ERROR in meanrotmpc: MPC has wrong number of terms' 63 call exit(201) 64 endif 65! 66! normal along the rotation axis 67! 68 dd=0.d0 69 do i=1,3 70 aa(i)=x(i,n) 71 dd=dd+aa(i)**2 72 enddo 73 dd=dsqrt(dd) 74 if(dd.lt.1.d-10) then 75 write(*,*) 76 & '*ERROR in meanrotmpc: rotation vector has zero length' 77 call exit(201) 78 endif 79 do i=1,3 80 aa(i)=aa(i)/dd 81 enddo 82c write(*,*) 'umpc_mean_rot aa',(aa(i),i=1,3) 83! 84! finding the center of gravity of the position and the 85! displacements of the nodes involved in the MPC 86! 87 do i=1,3 88 cgx(i)=0.d0 89 cgu(i)=0.d0 90 enddo 91! 92c write(*,*) 'umpc_mean_rot, nkn',nkn 93 do i=1,nkn 94 do j=1,3 95 cgx(j)=cgx(j)+x(j,3*i-2) 96 cgu(j)=cgu(j)+u(j,3*i-2) 97 enddo 98c write(*,*) 'umpc_mean_rot x',i,(x(j,3*i-2),j=1,3) 99c write(*,*) 'umpc_mean_rot u',i,(u(j,3*i-2),j=1,3) 100 enddo 101! 102 do i=1,3 103 cgx(i)=cgx(i)/nkn 104 cgu(i)=cgu(i)/nkn 105 enddo 106c write(*,*) 'umpc_mean_rot cgx',(cgx(i),i=1,3) 107! 108! calculating a standard deviation; this quantity will 109! serve as a limit for checking the closeness of individual 110! nodes to the center of gravity 111! 112 stdev=0.d0 113 do i=1,nkn 114 do j=1,3 115 stdev=stdev+(x(j,3*i-2)-cgx(j))**2 116 enddo 117 enddo 118 stdev=stdev/nkn 119c write(*,*) 'umpc_mean_rot stdev',stdev 120! 121! initializing a 122! 123 do i=1,n 124 a(i)=0.d0 125 enddo 126! 127! calculating the partial derivatives and storing them in a 128! 129 f=0.d0 130 do i=1,nkn 131! 132! relative positions 133! 134 do j=1,3 135 pi(j)=x(j,3*i-2)-cgx(j) 136 xi(j)=u(j,3*i-2)-cgu(j)+pi(j) 137 enddo 138c write(*,*) 'umpc_mean_rot pi',(pi(j),j=1,3) 139! 140! projection on a plane orthogonal to the rotation vector 141! 142 c1=pi(1)*aa(1)+pi(2)*aa(2)+pi(3)*aa(3) 143 do j=1,3 144 pi(j)=pi(j)-c1*aa(j) 145 enddo 146c write(*,*) 'umpc_mean_rot c1 pi',c1,(pi(j),j=1,3) 147! 148 c1=xi(1)*aa(1)+xi(2)*aa(2)+xi(3)*aa(3) 149 do j=1,3 150 xi(j)=xi(j)-c1*aa(j) 151 enddo 152! 153 c1=pi(1)*pi(1)+pi(2)*pi(2)+pi(3)*pi(3) 154c if(c1.lt.1.d-20) then 155 if(c1.lt.stdev*1.d-10) then 156 if(label(8:9).ne.'BS') then 157 write(*,*) '*WARNING in meanrotmpc: node ',jnode(3*i-2) 158 write(*,*) ' is very close to the ' 159 write(*,*) ' rotation axis through the' 160 write(*,*) ' center of gravity of' 161 write(*,*) ' the nodal cloud in a' 162 write(*,*) ' mean rotation MPC.' 163 write(*,*) ' This node is not taken' 164 write(*,*) ' into account in the MPC' 165 endif 166 cycle 167 endif 168 c3=xi(1)*xi(1)+xi(2)*xi(2)+xi(3)*xi(3) 169 c2=dsqrt(c1*c3) 170! 171 al=(aa(1)*pi(2)*xi(3)+aa(2)*pi(3)*xi(1)+aa(3)*pi(1)*xi(2) 172 & -aa(3)*pi(2)*xi(1)-aa(1)*pi(3)*xi(2)-aa(2)*pi(1)*xi(3)) 173 & /c2 174! 175 f=f+dasin(al) 176! 177 do j=1,3 178 if(j.eq.1) then 179 c4=aa(2)*pi(3)-aa(3)*pi(2) 180 elseif(j.eq.2) then 181 c4=aa(3)*pi(1)-aa(1)*pi(3) 182 else 183 c4=aa(1)*pi(2)-aa(2)*pi(1) 184 endif 185 c9=(c4/c2-al*xi(j)/c3)/dsqrt(1.d0-al*al) 186c write(*,*) 'umpc_mean_rot c4,c9',j,c4,c9 187! 188 do k=1,nkn 189 if(i.eq.k) then 190 c10=c9*(1.d0-1.d0/real(nkn)) 191 else 192 c10=-c9/real(nkn) 193 endif 194 a(k*3-3+j)=a(k*3-3+j)+c10 195c write(*,*) 'umpc_mean_rot c10',j,c10,a(k*3-3+j) 196 enddo 197 enddo 198 enddo 199c do j=1,n 200c write(*,*) 'umpc_mean_rot a ',a(j) 201c enddo 202 a(n)=-nkn 203 f=f-nkn*u(1,n) 204! 205! assigning the degrees of freedom 206! the first three dofs may not be in ascending order 207! 208 do i=1,3 209 b(i)=a(i) 210 enddo 211 do i=1,3 212 a(i)=b(jdof(i)) 213 enddo 214! 215 do i=4,nkn 216 jdof(i*3-2)=1 217 jdof(i*3-1)=2 218 jdof(i*3)=3 219 enddo 220! 221 jdof(n)=1 222! 223! looking for the maximum tangent to decide which DOF should be 224! taken to be the dependent one 225! 226 if(dabs(a(1)).lt.1.d-5) then 227! 228! special treatment of beam and shell mean rotation MPC's 229! cf. usermpc.f 230! 231 if(label(8:9).eq.'BS') then 232 do i=1,3 233 transcoef(i)=a(n-4+i) 234 enddo 235 endif 236! 237 imax=0 238 amax=1.d-5 239! 240 if(label(8:9).eq.'BS') then 241 nendnode=n-4 242 else 243 nendnode=n-1 244 endif 245 do i=2,nendnode 246 if(dabs(a(i)).gt.amax) then 247 idir=jdof(i) 248! 249 if(label(8:9).eq.'BS') then 250 if(a(i)*transcoef(idir).gt.0) cycle 251 endif 252! 253 node=jnode(i) 254 idof=8*(node-1)+idir 255! 256! check whether dof is not used in another MPC 257! 258 call nident(ikmpc,idof,nmpc,id) 259 if(id.gt.0) then 260 if(ikmpc(id).eq.idof) cycle 261 endif 262! 263! check whether dof is not used in a SPC 264! 265 call nident(ikboun,idof,nboun,id) 266 if(id.gt.0) then 267 if(ikboun(id).eq.idof) cycle 268 endif 269! 270 amax=dabs(a(i)) 271 imax=i 272 endif 273 enddo 274! 275 if(imax.eq.0) then 276 write(*,*) '*ERROR in umpc_mean_rot' 277 write(*,*) ' no mean rotation MPC can be' 278 write(*,*) ' generated for the MPC containing' 279 write(*,*) ' node ',jnode(1) 280 call exit(201) 281 endif 282! 283 nodeold=jnode(1) 284 idirold=jdof(1) 285 nodemax=jnode(imax) 286 idirmax=jdof(imax) 287! 288! exchange the node information, if needed 289! 290 if(nodemax.ne.nodeold) then 291 do i=1,3 292 iold(jdof(i))=i 293 enddo 294 do i=4,n-4 295 if(jnode(i).eq.nodemax) then 296 if(jdof(i).eq.1) then 297 inew(1)=i 298 elseif(jdof(i).eq.2) then 299 inew(2)=i 300 else 301 inew(3)=i 302 endif 303 endif 304 enddo 305! 306 do j=1,3 307 node=jnode(iold(j)) 308 idir=jdof(iold(j)) 309 xcoef=a(iold(j)) 310 jnode(iold(j))=jnode(inew(j)) 311 jdof(iold(j))=jdof(inew(j)) 312 a(iold(j))=a(inew(j)) 313 jnode(inew(j))=node 314 jdof(inew(j))=idir 315 a(inew(j))=xcoef 316 enddo 317 endif 318! 319! exchange the direction information, if needed 320! 321 iold(1)=1 322 do i=1,3 323 if(jdof(i).eq.idirmax) then 324 inew(1)=i 325 exit 326 endif 327 enddo 328! 329 if(iold(1).ne.inew(1)) then 330 do j=1,1 331 idir=jdof(iold(j)) 332 xcoef=a(iold(j)) 333 jdof(iold(j))=jdof(inew(j)) 334 a(iold(j))=a(inew(j)) 335 jdof(inew(j))=idir 336 a(inew(j))=xcoef 337 enddo 338 endif 339 endif 340! 341 return 342 end 343