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 updategeodata(nktet,netet_,h,d,dmin,ipoed,iedg,cotet, 20 & planfa,bc,cg,kontet,ifac,ipofa,doubleglob,integerglob,ipoeln) 21! 22! calculating the desired size of h in all nodes of the actual 23! mesh 24! 25 implicit none 26! 27 integer nktet,netet_,ipoed(*),iedg(3,*),kontet(4,*),ifac(4,*), 28 & ipofa(*),integerglob(*),ipoeln(*), 29 & nodef(3),nodes(4),iselect(1),nselect,nterms,nktri,nkon, 30 & nfield,nfaces,netet,nelem,ne,n1,n2,i,j,ialset(1),ielement, 31 & iendset(1),iface,imastset,istartset(1),index,konl(20), 32 & loopa 33! 34 real*8 h(*),d(*),dmin,cotet(3,*),planfa(4,*),bc(4,*),cg(3,*), 35 & doubleglob(*),x,y,z,r,ratio(20),p1x,p1y,p1z,p2x,p2y,p2z,p3x, 36 & p3y,p3z,p4x,p4y,p4z,a11,a12,a13,a21,a22,a23,a31,a32,a33,d1, 37 & d2,d3,det,coords(3),dist 38! 39! determine the size of all edges 40! 41 dmin=1.d30 42! 43 loop: do i=1,nktet 44 index=ipoed(i) 45 do 46 if(index.eq.0) cycle loop 47! 48 n1=iedg(1,index) 49 n2=iedg(2,index) 50! 51 d(index)=dsqrt((cotet(1,n1)-cotet(1,n2))**2+ 52 & (cotet(2,n1)-cotet(2,n2))**2+ 53 & (cotet(3,n1)-cotet(3,n2))**2) 54! 55 if(d(index).lt.dmin) dmin=d(index) 56! 57 index=iedg(3,index) 58 enddo 59 enddo loop 60! 61! calculating the desired edge size through interpolation 62! 63! initializing fields 64! 65 nktri=integerglob(1) 66 netet=integerglob(2) 67 ne=integerglob(3) 68 nkon=integerglob(4) 69 nfaces=integerglob(5) 70 nfield=1 71 nselect=1 72 iselect(1)=1 73 imastset=0 74! 75 do i=1,nktet 76 if(ipoeln(i).eq.0) cycle 77! 78! perform the interpolation for the internal node 79! 80 do j=1,3 81 coords(j)=cotet(j,i) 82 enddo 83 loopa=8 84 call basis(doubleglob(1),doubleglob(netet+1), 85 & doubleglob(2*netet+1), 86 & doubleglob(3*netet+1),doubleglob(4*netet+1), 87 & doubleglob(5*netet+1),integerglob(6), 88 & integerglob(netet+6), 89 & integerglob(2*netet+6),doubleglob(6*netet+1), 90 & integerglob(3*netet+6),nktri,netet, 91 & doubleglob(4*nfaces+6*netet+1),nfield, 92 & doubleglob(nktri+4*nfaces+6*netet+1), 93 & integerglob(7*netet+6),integerglob(ne+7*netet+6), 94 & integerglob(2*ne+7*netet+6), 95 & integerglob(nkon+2*ne+7*netet+6), 96 & coords(1),coords(2),coords(3),h(i),ratio,iselect, 97 & nselect,istartset,iendset,ialset,imastset, 98 & integerglob(nkon+2*ne+8*netet+6),nterms,konl,nelem, 99 & loopa,dist) 100 enddo 101! 102! updating the cg and bc fields 103! 104 do i=1,netet_ 105 ielement=i 106! 107 if(kontet(1,ielement).eq.0) cycle 108! 109 do j=1,4 110 nodes(j)=kontet(j,ielement) 111 enddo 112! 113! finding the center and radius of the circumscribed sphere 114! 115 p1x=cotet(1,nodes(1)) 116 p1y=cotet(2,nodes(1)) 117 p1z=cotet(3,nodes(1)) 118! 119 p2x=cotet(1,nodes(2)) 120 p2y=cotet(2,nodes(2)) 121 p2z=cotet(3,nodes(2)) 122! 123 p3x=cotet(1,nodes(3)) 124 p3y=cotet(2,nodes(3)) 125 p3z=cotet(3,nodes(3)) 126! 127 p4x=cotet(1,nodes(4)) 128 p4y=cotet(2,nodes(4)) 129 p4z=cotet(3,nodes(4)) 130! 131 a11=p1x-p2x 132 a12=p1y-p2y 133 a13=p1z-p2z 134 d1=((p1x+p2x)*a11+(p1y+p2y)*a12+(p1z+p2z)*a13)/2.d0 135! 136 a21=p1x-p3x 137 a22=p1y-p3y 138 a23=p1z-p3z 139 d2=((p1x+p3x)*a21+(p1y+p3y)*a22+(p1z+p3z)*a23)/2.d0 140! 141 a31=p1x-p4x 142 a32=p1y-p4y 143 a33=p1z-p4z 144 d3=((p1x+p4x)*a31+(p1y+p4y)*a32+(p1z+p4z)*a33)/2.d0 145! 146 det=a11*(a22*a33-a32*a23)-a12*(a21*a33-a31*a23)+ 147 & a13*(a21*a32-a31*a22) 148 x=(d1*(a22*a33-a23*a32)-d2*(a12*a33-a32*a13)+ 149 & d3*(a12*a23-a22*a13))/det 150 y=(-d1*(a21*a33-a31*a23)+d2*(a11*a33-a31*a13)- 151 & d3*(a11*a23-a21*a13))/det 152 z=(d1*(a21*a32-a31*a22)-d2*(a11*a32-a31*a12)+ 153 & d3*(a11*a22-a21*a12))/det 154! 155 r=dsqrt((x-p1x)**2+(y-p1y)**2+(z-p1z)**2) 156! 157 bc(1,ielement)=x 158 bc(2,ielement)=y 159 bc(3,ielement)=z 160 bc(4,ielement)=r 161! 162! determining the center of gravity of the element 163! 164 cg(1,ielement)=(p1x+p2x+p3x+p4x)/4.d0 165 cg(2,ielement)=(p1y+p2y+p3y+p4y)/4.d0 166 cg(3,ielement)=(p1z+p2z+p3z+p4z)/4.d0 167! 168 enddo 169! 170! update planfa field 171! 172 loop2:do i=1,nktet 173 iface=ipofa(i) 174 do 175 if(iface.eq.0) cycle loop2 176 nodef(1)=ifac(1,iface) 177 nodef(2)=ifac(2,iface) 178 nodef(3)=ifac(3,iface) 179! 180 call planeeq(cotet,nodef,planfa(1,iface)) 181! 182 iface=ifac(4,iface) 183 enddo 184 enddo loop2 185! 186 return 187 end 188 189