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 attach_3d(pneigh,pnode,nterms,ratio,dist,xil,etl,zel, 20 & loopa) 21! 22! ataches node with coordinates in "pnode" to the face containing 23! "nterms" nodes with coordinates in field "pneigh" (nterms < 20). 24! cave: the coordinates are stored in pneigh(1..3,*) 25! 26 implicit none 27! 28 integer nterms,i,j,k,m,imin,jmin,kmin,im,jm,km,n,loopa 29! 30 real*8 ratio(20),pneigh(3,20),pnode(3),a,xi(-1:1,-1:1,-1:1), 31 & et(-1:1,-1:1,-1:1),ze(-1:1,-1:1,-1:1),p(3),distmin,d1,dist, 32 & xil,etl,zel,dx(3),al,dummy 33! 34 n=3 35! 36 d1=1.d0 37! 38 xi(0,0,0)=0.d0 39 et(0,0,0)=0.d0 40 ze(0,0,0)=0.d0 41 call distattach_3d(xi(0,0,0),et(0,0,0),ze(0,0,0),pneigh,pnode,a,p, 42 & ratio,nterms) 43 distmin=a 44 imin=0 45 jmin=0 46 kmin=0 47! 48 do m=1,loopa 49! 50! initialisation 51! 52 d1=d1/10.d0 53! 54 do i=-1,1 55 do j=-1,1 56 do k=-1,1 57 if((i.eq.0).and.(j.eq.0).and.(k.eq.0)) cycle 58! 59 xi(i,j,k)=xi(0,0,0)+i*d1 60 et(i,j,k)=et(0,0,0)+j*d1 61 ze(i,j,k)=ze(0,0,0)+k*d1 62! 63! check whether inside the (-1,1)x(-1,1)x(-1,1) domain 64! 65 if((xi(i,j,k).le.1.d0).and. 66 & (xi(i,j,k).ge.-1.d0).and. 67 & (et(i,j,k).le.1.d0).and. 68 & (et(i,j,k).ge.-1.d0).and. 69 & (ze(i,j,k).le.1.d0).and. 70 & (ze(i,j,k).ge.-1.d0)) then 71 call distattach_3d(xi(i,j,k),et(i,j,k),ze(i,j,k), 72 & pneigh,pnode,a,p,ratio,nterms) 73! 74! checking for smallest initial distance 75! 76 if(a.lt.distmin) then 77 distmin=a 78 imin=i 79 jmin=j 80 kmin=k 81 endif 82 endif 83! 84 enddo 85 enddo 86 enddo 87! 88! minimizing the distance from the face to the node 89! 90 do 91! 92! exit if minimum found 93! 94c write(*,*) 'attach_3d',m,imin,jmin,kmin 95 if((imin.eq.0).and.(jmin.eq.0).and.(kmin.eq.0)) exit 96! 97! new center of 3x3x3 matrix 98! 99 xi(0,0,0)=xi(imin,jmin,kmin) 100 et(0,0,0)=et(imin,jmin,kmin) 101 ze(0,0,0)=ze(imin,jmin,kmin) 102! 103 im=imin 104 jm=jmin 105 km=kmin 106! 107 imin=0 108 jmin=0 109 kmin=0 110! 111 do i=-1,1 112 do j=-1,1 113 do k=-1,1 114 if((i+im.lt.-1).or.(i+im.gt.1).or. 115 & (j+jm.lt.-1).or.(j+jm.gt.1).or. 116 & (k+km.lt.-1).or.(k+km.gt.1)) then 117! 118 xi(i,j,k)=xi(0,0,0)+i*d1 119 et(i,j,k)=et(0,0,0)+j*d1 120 ze(i,j,k)=ze(0,0,0)+k*d1 121! 122! check whether inside the (-1,1)x(-1,1)x(-1,1) domain 123! 124 if((xi(i,j,k).le.1.d0).and. 125 & (xi(i,j,k).ge.-1.d0).and. 126 & (et(i,j,k).le.1.d0).and. 127 & (et(i,j,k).ge.-1.d0).and. 128 & (ze(i,j,k).le.1.d0).and. 129 & (ze(i,j,k).ge.-1.d0)) then 130 call distattach_3d(xi(i,j,k),et(i,j,k), 131 & ze(i,j,k),pneigh,pnode,a,p,ratio,nterms) 132! 133! check for new minimum 134! 135 if(a.lt.distmin) then 136 distmin=a 137 imin=i 138 jmin=j 139 kmin=k 140 endif 141 endif 142! 143 endif 144 enddo 145 enddo 146 enddo 147 enddo 148 enddo 149! 150 call distattach_3d(xi(0,0,0),et(0,0,0),ze(0,0,0),pneigh,pnode,a,p, 151 & ratio,nterms) 152! 153 do i=1,3 154 pnode(i)=p(i) 155 enddo 156! 157 dist=dsqrt(a) 158! 159 if((nterms.eq.20).or.(nterms.eq.8)) then 160 xil=xi(0,0,0) 161 etl=et(0,0,0) 162 zel=ze(0,0,0) 163 elseif((nterms.eq.4).or.(nterms.eq.10)) then 164 xil=(xi(0,0,0)+1.d0)/2.d0 165 etl=(et(0,0,0)+1.d0)/2.d0 166 zel=(ze(0,0,0)+1.d0)/2.d0 167 dx(1)=xil 168 dx(2)=etl 169 dx(3)=zel 170 call insertsortd(dx,n) 171c call dsort(dx,iy,n,kflag) 172 if(dx(3).gt.1.d-30) then 173 al=dx(3)/(xil+etl+zel) 174 xil=al*xil 175 etl=al*etl 176 zel=al*zel 177 endif 178 elseif((nterms.eq.6).or.(nterms.eq.15)) then 179 xil=(xi(0,0,0)+1.d0)/2.d0 180 etl=(et(0,0,0)+1.d0)/2.d0 181 if(xil+etl.gt.1.d0) then 182 dummy=xil 183 xil=1.d0-etl 184 etl=1.d0-dummy 185 endif 186 zel=ze(0,0,0) 187 endif 188! 189 return 190 end 191 192