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 geocfd(nef,ipkonf,konf,lakonf,co,coel,cofa,nface, 20 & ielfa,area,ipnei,neiel,xxn,xxi,xle,xlen,xlet,xrlfa,cosa, 21 & volume,neifa,xxj,cosb,hmin,ifatie,cs,tieset,icyclic,c,neij, 22 & physcon,isolidsurf,nsolidsurf,dy,xxni,xxnj,xxicn,nflnei, 23 & iturbulent,rf,yy,vel,velo,veloo,xxna,ale,alet,h) 24! 25! calculating geometric variables of the cells and their faces 26! 27 implicit none 28! 29 character*8 lakonf(*) 30 character*81 tieset(3,*) 31! 32 integer nef,ipkonf(*),konf(*),nface,ielfa(4,*),ipnei(*),neiel(*), 33 & ifaceq(8,6),i,j,k,indexe,kflag,index1,index2,j1,nope, 34 & nodes(4),iel1,iel2,iel3,iface,indexf,neifa(*),nf(5),ifacet(7,4), 35 & ifacew(8,5),ied4(2,6),ied6(2,9),ied8(2,12),ifatie(*), 36 & ics,itie,neighface,ifirst_occurrence,icyclic,neij(*),jface, 37 & isolidsurf(*),nsolidsurf,iopp8(4,6),iopp6(3,5),iopp4(1,4), 38 & node,nelem,nflnei,iturbulent,nx(nsolidsurf),ny(nsolidsurf), 39 & nz(nsolidsurf),kneigh,neigh(1),nfa(nsolidsurf) 40! 41 real*8 co(3,*),coel(3,*),cofa(3,*),area(*),xxn(3,*),xxi(3,*), 42 & xle(*),xlen(*),xlet(*),xrlfa(3,*),cosa(*),xsj2(3),xi,et, 43 & shp2(7,4),xs2(3,7),xl2(3,8),xl13,volume(*),dxsj2,xl(3,8), 44 & xxj(3,*),cosb(*),hmin,cs(17,*),xn(3),theta,pi,dc,ds,dd, 45 & c(3,3),diff(3),p(3),q(3),a(3),physcon(*),dy(*),xs(3,3), 46 & aa,bb,cc,dist,xxni(3,*),xxnj(3,*),xxicn(3,*),rf(3,*),x13(3), 47 & yy(*),x(nsolidsurf),y(nsolidsurf),z(nsolidsurf),xo(nsolidsurf), 48 & yo(nsolidsurf),zo(nsolidsurf),xp,yp,zp,vel(nef,0:7),h(*), 49 & velo(nef,0:7),veloo(nef,0:7),areal,xxna(3,*),ale(*),alet(*) 50! 51! nodes belonging to the cell faces 52! 53 data ifaceq /4,3,2,1,11,10,9,12, 54 & 5,6,7,8,13,14,15,16, 55 & 1,2,6,5,9,18,13,17, 56 & 2,3,7,6,10,19,14,18, 57 & 3,4,8,7,11,20,15,19, 58 & 4,1,5,8,12,17,16,20/ 59 data ifacet /1,3,2,7,6,5,11, 60 & 1,2,4,5,9,8,12, 61 & 2,3,4,6,10,9,13, 62 & 1,4,3,8,10,7,14/ 63 data ifacew /1,3,2,9,8,7,0,0, 64 & 4,5,6,10,11,12,0,0, 65 & 1,2,5,4,7,14,10,13, 66 & 2,3,6,5,8,15,11,14, 67 & 4,6,3,1,12,15,9,13/ 68 data nf /3,3,4,4,4/ 69 data ied4 /1,2,2,3,3,1,1,4,2,4,3,4/ 70 data ied6 /1,2,2,3,3,1,4,5,5,6,6,4,1,4,2,5,3,6/ 71 data ied8 /1,2,2,3,3,4,4,1,5,6,6,7,7,8,8,1,1,5,2,6,3,7,4,8/ 72 data iopp4 /4,3,1,2/ 73 data iopp6 /4,5,6,1,3,2,3,6,0,1,4,0,2,5,0/ 74 data iopp8 /5,6,7,8,4,3,2,1,3,4,8,7,4,1,5,8,1,2,6,5,2,3,7,6/ 75! 76 ifirst_occurrence=1 77 icyclic=0 78! 79! coordinates of the center of the cells 80! 81 do i=1,nef 82 if(ipkonf(i).lt.0) cycle 83 if(lakonf(i)(1:1).ne.'F') cycle 84 indexe=ipkonf(i) 85 if(lakonf(i)(4:4).eq.'8') then 86 nope=8 87 else if(lakonf(i)(4:4).eq.'6') then 88 nope=6 89 else 90 nope=4 91 endif 92 do j=1,3 93 do k=1,nope 94 coel(j,i)=coel(j,i)+co(j,konf(indexe+k)) 95 enddo 96 coel(j,i)=coel(j,i)/nope 97 enddo 98 enddo 99! 100 kflag=2 101! 102! loop over all faces 103! 104 do i=1,nface 105! 106! check for cyclic symmetry 107! 108 if(ifatie(i).ne.0) then 109 ics=abs(ifatie(i)) 110 itie=int(cs(17,ics)) 111 if(tieset(1,itie)(81:81).eq.'P') then 112 if(ifirst_occurrence.eq.1) then 113 do k=1,3 114 diff(k)=-cs(5+k,ics) 115 enddo 116 ifirst_occurrence=0 117 endif 118 elseif(tieset(1,itie)(81:81).eq.'Z') then 119 if(ifirst_occurrence.eq.1) then 120 icyclic=1 121 pi=4.d0*datan(1.d0) 122! 123! normal along the cyclic symmetry axis such that 124! the slave surface is rotated clockwise through the 125! body into the master surface while looking in 126! the direction of xn 127! 128 do k=1,3 129 a(k)=cs(5+k,ics) 130 xn(k)=cs(8+k,ics)-a(k) 131 enddo 132 dd=dsqrt(xn(1)*xn(1)+xn(2)*xn(2)+xn(3)*xn(3)) 133 do k=1,3 134 xn(k)=xn(k)/dd 135 enddo 136! 137! angle from the master to the slave surface 138! 139 theta=-2.d0*pi/cs(1,ics) 140! 141! rotation matrix rotating a vector in the master 142! surface into a vector in the slave surface 143! 144 dc=dcos(theta) 145 ds=dsin(theta) 146! 147! C-matrix from Guido Dhondt, The Finite Element 148! Method for Three-Dimensional Thermomechanical 149! Applications p 158 150! 151 c(1,1)=dc+(1.d0-dc)*xn(1)*xn(1) 152 c(1,2)= (1.d0-dc)*xn(1)*xn(2)-ds*xn(3) 153 c(1,3)= (1.d0-dc)*xn(1)*xn(3)+ds*xn(2) 154 c(2,1)= (1.d0-dc)*xn(2)*xn(1)+ds*xn(3) 155 c(2,2)=dc+(1.d0-dc)*xn(2)*xn(2) 156 c(2,3)= (1.d0-dc)*xn(2)*xn(3)-ds*xn(1) 157 c(3,1)= (1.d0-dc)*xn(3)*xn(1)-ds*xn(2) 158 c(3,2)= (1.d0-dc)*xn(3)*xn(2)+ds*xn(1) 159 c(3,3)=dc+(1.d0-dc)*xn(3)*xn(3) 160 ifirst_occurrence=0 161 endif 162 else 163 write(*,*) '*ERROR in geocfd' 164 write(*,*) ' kind of cyclic symmetry' 165 write(*,*) ' not known' 166 stop 167 endif 168 endif 169! 170 iel1=ielfa(1,i) 171 indexe=ipkonf(iel1) 172 j1=ielfa(4,i) 173 if(lakonf(iel1)(4:4).eq.'8') then 174! 175! hexahedral element 176! 177! coordinates of the face centers 178! 179 do j=1,4 180 nodes(j)=konf(indexe+ifaceq(j,j1)) 181 do k=1,3 182 xl2(k,j)=co(k,nodes(j)) 183 cofa(k,i)=cofa(k,i)+xl2(k,j) 184 enddo 185 enddo 186 do k=1,3 187 cofa(k,i)=cofa(k,i)/4.d0 188 enddo 189! 190 xi=0.d0 191 et=0.d0 192 call shape4q(xi,et,xl2,xsj2,xs2,shp2,kflag) 193! 194! area of the face 195! 196 dxsj2=dsqrt(xsj2(1)*xsj2(1)+xsj2(2)*xsj2(2)+ 197 & xsj2(3)*xsj2(3)) 198 area(i)=4.d0*dxsj2 199! 200 neighface=6 201! 202 else if(lakonf(iel1)(4:4).eq.'6') then 203! 204! wedge element 205! 206! coordinates of the face centers 207! 208 do j=1,nf(j1) 209 nodes(j)=konf(indexe+ifacew(j,j1)) 210 do k=1,3 211 xl2(k,j)=co(k,nodes(j)) 212 cofa(k,i)=cofa(k,i)+xl2(k,j) 213 enddo 214 enddo 215 do k=1,3 216 cofa(k,i)=cofa(k,i)/nf(j1) 217 enddo 218! 219 xi=0.d0 220 et=0.d0 221 if(nf(j1).eq.3) then 222 call shape3tri(xi,et,xl2,xsj2,xs2,shp2,kflag) 223 else 224 call shape4q(xi,et,xl2,xsj2,xs2,shp2,kflag) 225 endif 226! 227! area of the face 228! 229 dxsj2=dsqrt(xsj2(1)*xsj2(1)+xsj2(2)*xsj2(2)+ 230 & xsj2(3)*xsj2(3)) 231 if(nf(j1).eq.3) then 232 area(i)=dxsj2/2.d0 233 else 234 area(i)=4.d0*dxsj2 235 endif 236! 237 neighface=5 238! 239 else 240! 241! tetrahedral element 242! 243! coordinates of the face centers 244! 245 do j=1,3 246 nodes(j)=konf(indexe+ifacet(j,j1)) 247 do k=1,3 248 xl2(k,j)=co(k,nodes(j)) 249 cofa(k,i)=cofa(k,i)+xl2(k,j) 250 enddo 251 enddo 252 do k=1,3 253 cofa(k,i)=cofa(k,i)/3.d0 254 enddo 255! 256 xi=0.d0 257 et=0.d0 258 call shape3tri(xi,et,xl2,xsj2,xs2,shp2,kflag) 259! 260! area of the face 261! 262 dxsj2=dsqrt(xsj2(1)*xsj2(1)+xsj2(2)*xsj2(2)+ 263 & xsj2(3)*xsj2(3)) 264 area(i)=dxsj2/2.d0 265! 266 neighface=4 267! 268 endif 269! 270! normal and xi-vector on face viewed from cell 1 271! 272 index1=ipnei(iel1)+j1 273 do k=1,3 274 xxn(k,index1)=xsj2(k)/dxsj2 275 xxi(k,index1)=cofa(k,i)-coel(k,iel1) 276 rf(k,i)=xxi(k,index1) 277 enddo 278! 279! distance from face center to the center of cell 1 280! 281 xle(index1)=dsqrt(xxi(1,index1)**2+xxi(2,index1)**2+ 282 & xxi(3,index1)**2) 283 do k=1,3 284 xxi(k,index1)=xxi(k,index1)/xle(index1) 285 enddo 286! 287! angle between the normal and the xi-vector 288! 289 cosa(index1)=xxn(1,index1)*xxi(1,index1)+ 290 & xxn(2,index1)*xxi(2,index1)+ 291 & xxn(3,index1)*xxi(3,index1) 292! 293 iel2=ielfa(2,i) 294! 295! check whether there is an adjacent cell 296! 297 if(iel2.ne.0) then 298 index2=ipnei(iel2)+neij(index1) 299! 300! normal and xi-vector on face viewed from cell 2 301! 302 if(ifatie(i).eq.0) then 303! 304! genuine neighbor 305! 306 do k=1,3 307 xxi(k,index2)=cofa(k,i)-coel(k,iel2) 308 xxn(k,index2)=-xxn(k,index1) 309 enddo 310! 311 xle(index2)=dsqrt(xxi(1,index2)**2+xxi(2,index2)**2+ 312 & xxi(3,index2)**2) 313 do k=1,3 314 xxi(k,index2)=xxi(k,index2)/xle(index2) 315 enddo 316! 317! angle between the normal and the xi-vector: xxn.xxi 318! 319 cosa(index2)=xxn(1,index2)*xxi(1,index2)+ 320 & xxn(2,index2)*xxi(2,index2)+ 321 & xxn(3,index2)*xxi(3,index2) 322! 323! distance from the face center to the center of the 324! adjacent cell 325! 326 xlen(index1)=xle(index2) 327 xlen(index2)=xle(index1) 328! 329 do k=1,3 330 xxj(k,index2)=coel(k,iel1)-coel(k,iel2) 331 enddo 332! 333! distance between the cell center and the center of the 334! adjacent cell 335! 336 xlet(index1)=dsqrt(xxj(1,index2)**2+xxj(2,index2)**2 337 & +xxj(3,index2)**2) 338 xlet(index2)=xlet(index1) 339! 340! xxj is the unit vector connecting neighboring cell centers 341! 342 do k=1,3 343 xxj(k,index2)=xxj(k,index2)/xlet(index2) 344 xxj(k,index1)=-xxj(k,index2) 345 enddo 346! 347! xxn.xxj 348! 349 cosb(index2)=xxn(1,index1)*xxj(1,index1)+ 350 & xxn(2,index1)*xxj(2,index1)+ 351 & xxn(3,index1)*xxj(3,index1) 352 cosb(index1)=cosb(index2) 353! 354c xrlfa(1,i)=xle(index2)/(xle(index1)+xle(index2)) 355c xrlfa(2,i)=xle(index1)/(xle(index1)+xle(index2)) 356 else 357! 358! cyclic symmetry face: some quantities are 359! calculated on the cyclic symmetric face 360! 361 xlen(index2)=xle(index1) 362! 363! rotational cyclic symmetry 364! 365! vector from the axis to the center of the 366! cyclic symmetric cell and orthogonal to the 367! axis 368! 369 if(tieset(1,itie)(81:81).eq.'Z') then 370 do k=1,3 371 p(k)=coel(k,iel2)-a(k) 372 enddo 373 dd=p(1)*xn(1)+p(2)*xn(2)+p(3)*xn(3) 374 do k=1,3 375 p(k)=p(k)-dd*xn(k) 376 enddo 377! 378 if(ifatie(i).gt.0) then 379! 380! vector rotated in the direction of the slave surface 381! (iel2 is adjacent to the master surface) 382! 383 do k=1,3 384 q(k)=c(k,1)*p(1)+c(k,2)*p(2)+c(k,3)*p(3) 385 enddo 386 else 387! 388! vector rotated in the direction of the master surface 389! (iel2 is adjacent to the slave surface) 390! 391 do k=1,3 392 q(k)=c(1,k)*p(1)+c(2,k)*p(2)+c(3,k)*p(3) 393 enddo 394 endif 395! 396! vector connecting the center of the cyclic 397! symmetry cell with the center of its rotated 398! ghost cell 399! 400 do k=1,3 401 diff(k)=q(k)-p(k) 402 enddo 403 endif 404! 405 do k=1,3 406 xxj(k,index1)=coel(k,iel2)-coel(k,iel1) 407 & +diff(k) 408 enddo 409! 410! distance between the cell center and the center of the 411! adjacent cell 412! 413 xlet(index1)=dsqrt(xxj(1,index1)**2+xxj(2,index1)**2 414 & +xxj(3,index1)**2) 415! 416! xxj is the unit vector connecting neighboring cell centers 417! 418 do k=1,3 419 xxj(k,index1)=xxj(k,index1)/xlet(index1) 420 enddo 421! 422! xxn.xxj 423! 424 cosb(index1)=xxn(1,index1)*xxj(1,index1)+ 425 & xxn(2,index1)*xxj(2,index1)+ 426 & xxn(3,index1)*xxj(3,index1) 427 endif 428! 429! calculating rf = the shortest vector between the 430! center of the face and the line connecting 431! the centers of the adjacent elements. The 432! corresponding point on this line is called p 433! calculating xrlfa = the ratio of the segments created by 434! the point p to the total distance between 435! the centers of the adjacent elements of 436! the face 437! 438 dd=rf(1,i)*xxj(1,index1)+ 439 & rf(2,i)*xxj(2,index1)+ 440 & rf(3,i)*xxj(3,index1) 441! 442 do k=1,3 443 rf(k,i)=rf(k,i)-dd*xxj(k,index1) 444 enddo 445! 446 xrlfa(2,i)=dd/xlet(index1) 447 xrlfa(1,i)=1.d0-xrlfa(2,i) 448 else 449! 450! xxi and xxj coincide 451! 452 do k=1,3 453 xxj(k,index1)=xxi(k,index1) 454 enddo 455 cosb(index1)=cosa(index1) 456! 457! external face: determining the cell next to the 458! adjacent cell 459! 460 iel3=ielfa(3,i) 461 if(iel3.eq.0) cycle 462c xl13=dsqrt((coel(1,iel1)-coel(1,iel3))**2+ 463c & (coel(2,iel1)-coel(2,iel3))**2+ 464c & (coel(3,iel1)-coel(3,iel3))**2) 465c xrlfa(1,i)=(xl13+xle(index1))/xl13 466c xrlfa(3,i)=1.d0-xrlfa(1,i) 467! 468! unit vector pointing from the center in element iel3 469! to the center in element iel1 470! 471 do k=1,3 472 x13(k)=coel(k,iel1)-coel(k,iel3) 473 enddo 474! 475 xl13=dsqrt(x13(1)*x13(1)+x13(2)*x13(2)+x13(3)*x13(3)) 476! 477 do k=1,3 478 x13(k)=x13(k)/xl13 479 enddo 480! 481 dd=rf(1,i)*x13(1)+rf(2,i)*x13(2)+rf(3,i)*x13(3) 482! 483 do k=1,3 484 rf(k,i)=rf(k,i)-dd*x13(k) 485 enddo 486! 487 xrlfa(3,i)=-dd/xl13 488 xrlfa(1,i)=1.d0-xrlfa(3,i) 489! 490 endif 491 enddo 492! 493! for cyclic symmetric faces xrlfa has not been filled yet 494! 495c if(ifirst_occurrence.eq.0) then 496c do i=1,nface 497c if(ifatie(i).ne.0) then 498c index1=ipnei(ielfa(1,i))+ielfa(4,i) 499c xrlfa(1,i)=xlen(index1)/(xle(index1)+xlen(index1)) 500c xrlfa(2,i)=1.d0-xrlfa(1,i) 501c endif 502c enddo 503c endif 504! 505! calculation of the volume of the elements 506! 507 do i=1,nef 508 if(ipkonf(i).lt.0) cycle 509 if(lakonf(i)(1:1).ne.'F') cycle 510 indexf=ipnei(i) 511 volume(i)=0.d0 512 do j=1,ipnei(i+1)-ipnei(i) 513 iface=neifa(indexf+j) 514 volume(i)=volume(i)+ 515 & area(iface)*cofa(1,iface)*xxn(1,indexf+j) 516 enddo 517 enddo 518! 519! calculation of the minimum length within the cells 520! 521c hmin=1.d30 522c do i=1,nef 523c indexe=ipkonf(i) 524c read(lakonf(i)(4:4),'(i1)') nope 525c do j=1,nope 526c do k=1,3 527c xl(k,j)=co(k,konf(indexe+j)) 528c enddo 529c enddo 530c if(nope.eq.4) then 531c do j=1,6 532c hmin=min(hmin,(xl(1,ied4(1,j))-xl(1,ied4(2,j)))**2+ 533c & (xl(2,ied4(1,j))-xl(2,ied4(2,j)))**2+ 534c & (xl(3,ied4(1,j))-xl(3,ied4(2,j)))**2) 535c enddo 536c elseif(nope.eq.6) then 537c do j=1,9 538c hmin=min(hmin,(xl(1,ied6(1,j))-xl(1,ied6(2,j)))**2+ 539c & (xl(2,ied6(1,j))-xl(2,ied6(2,j)))**2+ 540c & (xl(3,ied6(1,j))-xl(3,ied6(2,j)))**2) 541c enddo 542c else 543c do j=1,12 544c hmin=min(hmin,(xl(1,ied8(1,j))-xl(1,ied8(2,j)))**2+ 545c & (xl(2,ied8(1,j))-xl(2,ied8(2,j)))**2+ 546c & (xl(3,ied8(1,j))-xl(3,ied8(2,j)))**2) 547c enddo 548c endif 549c enddo 550c hmin=dsqrt(hmin) 551c write(*,*) 'hmin first ',hmin 552! 553! calculation of the minimum height within the cells 554! 555 hmin=1.d30 556 do i=1,nef 557 if(ipkonf(i).lt.0) cycle 558 if(lakonf(i)(1:1).ne.'F') cycle 559 if(lakonf(i)(4:4).eq.'8') then 560 nope=8 561 else if(lakonf(i)(4:4).eq.'6') then 562 nope=6 563 else 564 nope=4 565 endif 566 indexf=ipnei(i) 567c write(*,*) 'geocfd ',i,volume(i) 568 do j=1,ipnei(i+1)-ipnei(i) 569 indexf=indexf+1 570 iface=neifa(indexf) 571 if(nope.eq.4) then 572 h(indexf)=volume(i)*3.d0/area(iface) 573 elseif(nope.eq.6) then 574 if(j.le.2) then 575 h(indexf)=volume(i)/area(iface) 576 else 577 h(indexf)=volume(i)*2.d0/area(iface) 578 endif 579 else 580 h(indexf)=volume(i)/area(iface) 581c write(*,*) 'geocfd ',indexf,area(iface),h(indexf) 582 endif 583 hmin=min(h(indexf),hmin) 584 enddo 585 enddo 586c write(*,*) 'hmin second ',hmin 587! 588! calculate the distance to the nearest node for solid surface 589! faces 590! 591 if(iturbulent.gt.0) then 592! 593 if(dabs(physcon(5)).le.0.d0) then 594 write(*,*) '*ERROR in geocfd: velocity at infinity' 595 write(*,*) ' is nonpositive;' 596 write(*,*) ' wrong *VALUES AT INFINITY' 597 call exit(201) 598 endif 599! 600 if(dabs(physcon(7)).le.0.d0) then 601 write(*,*) '*ERROR in geocfd: density at infinity' 602 write(*,*) ' is nonpositive;' 603 write(*,*) ' wrong *VALUES AT INFINITY' 604 call exit(201) 605 endif 606! 607 if(dabs(physcon(8)).le.0.d0) then 608 write(*,*) '*ERROR in geocfd: length of the ' 609 write(*,*) ' computational domain is nonpositive;' 610 write(*,*) ' wrong *VALUES AT INFINITY' 611 call exit(201) 612 endif 613! 614 do i=1,nsolidsurf 615 iface=isolidsurf(i) 616 nelem=int(iface/10) 617 indexe=ipkonf(nelem) 618 jface=iface-nelem*10 619! 620! xl contains the coordinates of the nodes belonging 621! to the face 622! 623 if(lakonf(nelem)(4:4).eq.'8') then 624 do j=1,4 625 node=konf(indexe+ifaceq(j,jface)) 626 do k=1,3 627 xl(k,j)=co(k,node) 628 enddo 629 enddo 630 elseif(lakonf(nelem)(4:4).eq.'6') then 631 if(jface.le.2) then 632 do j=1,3 633 node=konf(indexe+ifacew(j,jface)) 634 do k=1,3 635 xl(k,j)=co(k,node) 636 enddo 637 enddo 638 else 639 do j=1,4 640 node=konf(indexe+ifacew(j,jface)) 641 do k=1,3 642 xl(k,j)=co(k,node) 643 enddo 644 enddo 645 endif 646 else 647 do j=1,3 648 node=konf(indexe+ifacet(j,jface)) 649 do k=1,3 650 xl(k,j)=co(k,node) 651 enddo 652 enddo 653 endif 654! 655! determine the face number in field ielfa 656! 657 nfa(i)=neifa(ipnei(nelem)+jface) 658! 659! determine the plane through the face (exact for 660! 3-node face, approximate for a 4-node face) 661! 662 if((lakonf(nelem)(4:4).eq.'8').or. 663 & ((lakonf(nelem)(4:4).eq.'6').and.(jface.gt.2))) then 664! 665! computation of the local derivative of the global coordinates 666! (xs) 667! 668 do j=1,3 669 xs(j,1)=-xl(j,1)+xl(j,2)+xl(j,3)-xl(j,4) 670 xs(j,2)=-xl(j,1)-xl(j,2)+xl(j,3)+xl(j,4) 671 enddo 672! 673! computation of the jacobian vector for xi,et=0 674! 675 aa=xs(2,1)*xs(3,2)-xs(3,1)*xs(2,2) 676 bb=xs(1,2)*xs(3,1)-xs(3,2)*xs(1,1) 677 cc=xs(1,1)*xs(2,2)-xs(2,1)*xs(1,2) 678 dd=dsqrt(aa*aa+bb*bb+cc*cc) 679 aa=aa/dd 680 bb=bb/dd 681 cc=cc/dd 682 dd=-(aa*(xl(1,1)+xl(1,2)+xl(1,3)+xl(1,4)) 683 & +bb*(xl(2,1)+xl(2,2)+xl(2,3)+xl(2,4)) 684 & +cc*(xl(3,1)+xl(3,2)+xl(3,3)+xl(3,4)))/4.d0 685 else 686! 687! computation of the local derivative of the global coordinates 688! (xs) 689! 690 do j=1,3 691 xs(j,1)=-xl(j,1)+xl(j,2) 692 xs(j,2)=-xl(j,1)+xl(j,3) 693 enddo 694! 695! computation of the jacobian vector (unique for triangle) 696! 697 aa=xs(2,1)*xs(3,2)-xs(3,1)*xs(2,2) 698 bb=xs(1,2)*xs(3,1)-xs(3,2)*xs(1,1) 699 cc=xs(1,1)*xs(2,2)-xs(2,1)*xs(1,2) 700 dd=dsqrt(aa*aa+bb*bb+cc*cc) 701 aa=aa/dd 702 bb=bb/dd 703 cc=cc/dd 704 dd=-(aa*xl(1,1)+bb*xl(2,1)+cc*xl(3,1)) 705 endif 706! 707! determine the shortest distance within the element 708! from the solid surface face 709! 710 dist=1.d30 711 if(lakonf(nelem)(4:4).eq.'8') then 712 do j=1,4 713 node=konf(indexe+iopp8(j,jface)) 714 dist=min(dist,-(aa*co(1,node)+bb*co(2,node) 715 & +cc*co(3,node)+dd)) 716 enddo 717 elseif(lakonf(nelem)(4:4).eq.'6') then 718 if(jface.le.2) then 719 do j=1,3 720 node=konf(indexe+iopp6(j,jface)) 721 dist=min(dist,-(aa*co(1,node)+bb*co(2,node) 722 & +cc*co(3,node)+dd)) 723 enddo 724 else 725 do j=1,2 726 node=konf(indexe+iopp6(j,jface)) 727 dist=min(dist,-(aa*co(1,node)+bb*co(2,node) 728 & +cc*co(3,node)+dd)) 729 enddo 730 endif 731 else 732 node=konf(indexe+iopp4(1,jface)) 733 dist=min(dist,-(aa*co(1,node)+bb*co(2,node) 734 & +cc*co(3,node)+dd)) 735 endif 736! 737! 60.d0/(0.075d0*delta(y)**2) 738! 739 dy(i)=800.d0/(dist*dist) 740 enddo 741 endif 742! 743! calculate for each element center the shortest distance to a solid 744! surface (only for the BSL and SST turbulence model) 745! 746 if(iturbulent.gt.2) then 747! 748 do i=1,nsolidsurf 749 iface=nfa(i) 750 x(i)=cofa(1,iface) 751 y(i)=cofa(2,iface) 752 z(i)=cofa(3,iface) 753 xo(i)=x(i) 754 yo(i)=y(i) 755 zo(i)=z(i) 756 nx(i)=i 757 ny(i)=i 758 nz(i)=i 759 enddo 760! 761 kflag=2 762 call dsort(x,nx,nsolidsurf,kflag) 763 call dsort(y,ny,nsolidsurf,kflag) 764 call dsort(z,nz,nsolidsurf,kflag) 765! 766 kneigh=1 767 do i=1,nef 768 xp=coel(1,i) 769 yp=coel(2,i) 770 zp=coel(3,i) 771 call near3d(xo,yo,zo,x,y,z,nx,ny,nz,xp,yp,zp,nsolidsurf, 772 & neigh,kneigh) 773 yy(i)=dsqrt((xp-xo(neigh(1)))**2+ 774 & (yp-yo(neigh(1)))**2+ 775 & (zp-zo(neigh(1)))**2) 776 enddo 777 endif 778! 779! auxiliary fields 780! 781 do i=1,nflnei 782 areal=area(neifa(i)) 783 do k=1,3 784 xxni(k,i)=xxn(k,i)-xxi(k,i) 785 xxnj(k,i)=(xxn(k,i)-xxj(k,i))*areal 786 xxicn(k,i)=xxn(k,i)-xxj(k,i)/cosb(i) 787 xxna(k,i)=xxn(k,i)*areal 788 enddo 789 ale(i)=areal/xle(i) 790 alet(i)=areal/xlet(i) 791 cosa(i)=ale(i)/cosa(i) 792 enddo 793c 794c initial conditions 795c 796c do i=1,nef 797c vel(i,0)=1.d0+coel(2,i)*(1.d0-coel(2,i))/2.d0 798c vel(i,1)=coel(2,i) 799c vel(i,2)=0.d0 800c vel(i,3)=0.d0 801c vel(i,4)=1.d0 802cc do j=0,4 803cc velo(i,j)=vel(i,j) 804cc veloo(i,j)=vel(i,j) 805cc enddo 806c enddo 807c write(*,*) 'geocfd neifa,neiel' 808c do i=1,6*nef 809c write(*,*) (i-1)/6+1,i-6*((i-1)/6),neifa(i),neiel(i) 810c enddo 811c write(*,*) 'geocfd xle,xlen,xlet' 812c do i=1,6*nef 813c write(*,*) (i-1)/6+1,i-6*((i-1)/6),xle(i),xlen(i),xlet(i) 814c enddo 815c write(*,*) 'geocfd xxn' 816c do i=1,6*nef 817c write(*,*) (i-1)/6+1,i-6*((i-1)/6),(xxn(j,i),j=1,3) 818c enddo 819c write(*,*) 'geocfd xxi' 820c do i=1,6*nef 821c write(*,*) (i-1)/6+1,i-6*((i-1)/6),(xxi(j,i),j=1,3) 822c enddo 823c write(*,*) 'geocfd xxj' 824c do i=1,6*nef 825c write(*,*) (i-1)/6+1,i-6*((i-1)/6),(xxj(j,i),j=1,3) 826c enddo 827c write(*,*) 'geocfd cosa,cosb' 828c do i=1,6*nef 829c write(*,*) (i-1)/6+1,i-6*((i-1)/6),cosa(i),cosb(i) 830c enddo 831c do i=1,nef 832c write(*,*) 'coef ',i,(coel(j,i),j=1,3) 833c enddo 834c do i=1,nface 835c write(*,*) 'cofa ',i,(cofa(j,i),j=1,3) 836c enddo 837c do i=1,nface 838c write(*,*) 'rf ',i,(rf(j,i),j=1,3) 839c enddo 840! 841 return 842 end 843