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! 20! Subroutine x_interpolate.f 21! 22! Triangulates the face and interpolates xstate variables according to 23! the plane equations of these planes 24! 25! by: Jaro Hokkanen 26! 27! 28 subroutine interpolateinface(kk,xstate,xstateini,numpts,nstate_, 29 & mi,islavsurf,pslavsurf, 30 & ne0,islavsurfold,pslavsurfold) 31! 32 implicit none 33! 34 integer numpts,i_int,n_int,i,j,k,kk,l,ll,nn, 35 & koncont(3,2*numpts+1),itri,kflag,neigh(1),kneigh, 36 & imastop(3,2*numpts+1),indexcj,nopespringj,list(numpts), 37 & igauss,mi(*),nstate_,itriangle(100),itriold, 38 & ifaceq(8,6),ip(numpts),ne0,itrinew,ntriangle, 39 & ifacet(6,4),ifacew1(4,5),ifacew2(8,5),n,islavsurf(2,*), 40 & ibin(numpts),ivert1,ntriangle_,nterms,m,islavsurfold(2,*), 41 & nx(2*numpts+1),ny(2*numpts+1),isol,id 42! 43 real*8 xstate(nstate_,mi(1),*),p(3),pslavsurfold(3,*), 44 & xstateini(nstate_,mi(1),*),coi(2,numpts+3),pneigh(3,3), 45 & cg(2,2*numpts+1),pslavsurf(3,*),xil,etl,ratio(3), 46 & z(9,3),x(2*numpts+1),xo(2*numpts+1),dis(3), 47 & y(2*numpts+1),yo(2*numpts+1),straight(9,2*numpts+1),dist 48! 49! nodes per face for hex elements 50! 51 data ifaceq /4,3,2,1,11,10,9,12, 52 & 5,6,7,8,13,14,15,16, 53 & 1,2,6,5,9,18,13,17, 54 & 2,3,7,6,10,19,14,18, 55 & 3,4,8,7,11,20,15,19, 56 & 4,1,5,8,12,17,16,20/ 57! 58! nodes per face for tet elements 59! 60 data ifacet /1,3,2,7,6,5, 61 & 1,2,4,5,9,8, 62 & 2,3,4,6,10,9, 63 & 1,4,3,8,10,7/ 64! 65! nodes per face for linear wedge elements 66! 67 data ifacew1 /1,3,2,0, 68 & 4,5,6,0, 69 & 1,2,5,4, 70 & 2,3,6,5, 71 & 3,1,4,6/ 72! 73! nodes per face for quadratic wedge elements 74! 75 data ifacew2 /1,3,2,9,8,7,0,0, 76 & 4,5,6,10,11,12,0,0, 77 & 1,2,5,4,7,14,10,13, 78 & 2,3,6,5,8,15,11,14, 79 & 3,1,4,6,9,13,12,15/ 80! 81 kneigh=1 82! 83 do i=1,numpts 84 list(i)=i 85 enddo 86! 87! Loop over the old integration points within the face 88! 89 ll=0 90 do l=islavsurfold(2,kk)+1,islavsurfold(2,kk+1) 91 ll=ll+1 92 coi(1,ll)=pslavsurfold(1,l) 93 coi(2,ll)=pslavsurfold(2,l) 94 ip(ll)=ne0+l 95 enddo 96! 97! Calling deltri for triangulation 98! 99 kflag=2 100! 101 call deltri(numpts,numpts,coi(1,1:numpts+3),coi(2,1:numpts+3), 102 & list,ibin,koncont,imastop,n) 103! 104! rearranging imastop field: imastop(1,i) is the triangle opposite 105! of local node 1 in triangle i etc.. 106! 107 nn=0 108 do i=1,n 109 ivert1=imastop(1,i) 110 imastop(1,i)=imastop(2,i) 111 imastop(2,i)=imastop(3,i) 112 imastop(3,i)=ivert1 113 enddo 114! 115! determining the center of gravity and the bounding planes 116! of the triangles 117! 118 call updatecont2d(koncont,n,coi,cg,straight) 119! 120! sorting the centers of gravity 121! 122 do l=1,n 123 xo(l)=cg(1,l) 124 x(l)=xo(l) 125 nx(l)=l 126 yo(l)=cg(2,l) 127 y(l)=yo(l) 128 ny(l)=l 129 enddo 130 kflag=2 131 call dsort(x,nx,n,kflag) 132 call dsort(y,ny,n,kflag) 133! 134! Loop over the new integration points 135! 136 do igauss=islavsurf(2,kk)+1,islavsurf(2,kk+1) 137! 138! coordinates of the new integration point 139! 140 p(1)=pslavsurf(1,igauss) 141 p(2)=pslavsurf(2,igauss) 142! 143! closest triangle center of gravity 144! 145 call near2d(xo,yo,x,y,nx,ny,p(1),p(2), 146 & n,neigh,kneigh) 147! 148 isol=0 149! 150 itriold=0 151 itri=neigh(1) 152 ntriangle=0 153 ntriangle_=100 154! 155 loop1: do 156 do l=1,3 157 ll=3*l-2 158 dist=straight(ll,itri)*p(1)+ 159 & straight(ll+1,itri)*p(2)+ 160 & straight(ll+2,itri) 161! 162 if(dist.gt.0.d0) then 163 itrinew=imastop(l,itri) 164 if(itrinew.eq.0) then 165c write(*,*) '**border reached' 166 exit loop1 167 elseif(itrinew.eq.itriold) then 168c write(*,*) '**solution in between triangles' 169 isol=itri 170 exit loop1 171 else 172 call nident(itriangle,itrinew,ntriangle,id) 173 if(id.gt.0) then 174 if(itriangle(id).eq.itrinew) then 175c write(*,*) '**circular path;no solution' 176 exit loop1 177 endif 178 endif 179 ntriangle=ntriangle+1 180 if(ntriangle.gt.ntriangle_) then 181c write(*,*) '**too many iterations' 182 exit loop1 183 endif 184 do k=ntriangle,id+2,-1 185 itriangle(k)=itriangle(k-1) 186 enddo 187 itriangle(id+1)=itrinew 188 itriold=itri 189 itri=itrinew 190 cycle loop1 191 endif 192 elseif(l.eq.3) then 193c write(*,*) '**regular solution' 194 isol=itri 195 exit loop1 196 endif 197 enddo 198 enddo loop1 199! 200 if(isol.eq.0) then 201! 202! mapping the new integration point onto the nearest 203! triangle 204! 205! calculating the distance from each edge of the triangle 206! dis(i) is the distance of the new integration point fromt 207! the edge opposite to local node i of the triangle itri 208! 209 do l=1,3 210 ll=3*l-2 211 dis(l)=straight(ll,itri)*p(1)+ 212 & straight(ll+1,itri)*p(2)+ 213 & straight(ll+2,itri) 214 enddo 215! 216 if(dis(1).gt.0.d0) then 217 if(dis(2).gt.0.d0) then 218! 219! projection is node 3 of the triangle 220! 221 p(1)=coi(1,koncont(3,itri)) 222 p(2)=coi(2,koncont(3,itri)) 223 elseif(dis(3).gt.0.d0) then 224! 225! projection is node 2 of the triangle 226! 227 p(1)=coi(1,koncont(2,itri)) 228 p(2)=coi(2,koncont(2,itri)) 229 else 230! 231! projection onto local edge 1 (opposite to local node 1) 232! 233 p(1)=p(1)-dis(1)*straight(1,itri) 234 p(2)=p(2)-dis(1)*straight(2,itri) 235 endif 236 elseif(dis(2).gt.0.d0) then 237 if(dis(3).gt.0.d0) then 238! 239! projection is node 1 of the triangle 240! 241 p(1)=coi(1,koncont(1,itri)) 242 p(2)=coi(2,koncont(1,itri)) 243 else 244! 245! projection onto local edge 2 (opposite to local node 2) 246! 247 p(1)=p(1)-dis(2)*straight(4,itri) 248 p(2)=p(2)-dis(2)*straight(5,itri) 249 endif 250 else 251! 252! projection onto local edge 3 (opposite to local node 3) 253! 254 p(1)=p(1)-dis(3)*straight(7,itri) 255 p(2)=p(2)-dis(3)*straight(8,itri) 256 endif 257c do k=1,3 258c do m=1,2 259c pneigh(m,k)=coi(m,koncont(k,itri)) 260c enddo 261c pneigh(3,k)=0.d0 262c enddo 263c p(3)=0.d0 264c nterms=3 265c ! 266c call attach_2d(pneigh,p,nterms,ratio,dist,xil,etl) 267c ! 268c do m=1,2 269c p(m)=0.d0 270c do k=1,3 271c p(m)=p(m)+ratio(k)*pneigh(m,k) 272c enddo 273c enddo 274c 275 endif 276! 277! Assigning xstate values from the vertices of the triangle 278! to the field z (integration point values from xstateini) 279! 280 do k=1,3 281 do i=1,9 282 z(i,k)=xstateini(i,1,ip(koncont(k,itri))) 283 enddo 284 enddo 285! 286! Calling plane_eq to interpolate values using plane equation 287! 288 do i=1,9 289 call plane_eq( 290 & coi(1,koncont(1,itri)),coi(2,koncont(1,itri)),z(i,1), 291 & coi(1,koncont(2,itri)),coi(2,koncont(2,itri)),z(i,2), 292 & coi(1,koncont(3,itri)),coi(2,koncont(3,itri)),z(i,3), 293 & p(1),p(2),xstate(i,1,ne0+igauss)) 294 enddo 295! 296 enddo 297! 298 return 299 end 300 301