1!/*****************************************************************************/ 2! * 3! * Elmer, A Finite Element Software for Multiphysical Problems 4! * 5! * Copyright 1st April 1995 - , CSC - IT Center for Science Ltd., Finland 6! * 7! * This library is free software; you can redistribute it and/or 8! * modify it under the terms of the GNU Lesser General Public 9! * License as published by the Free Software Foundation; either 10! * version 2.1 of the License, or (at your option) any later version. 11! * 12! * This library is distributed in the hope that it will be useful, 13! * but WITHOUT ANY WARRANTY; without even the implied warranty of 14! * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 15! * Lesser General Public License for more details. 16! * 17! * You should have received a copy of the GNU Lesser General Public 18! * License along with this library (in file ../LGPL-2.1); if not, write 19! * to the Free Software Foundation, Inc., 51 Franklin Street, 20! * Fifth Floor, Boston, MA 02110-1301 USA 21! * 22! *****************************************************************************/ 23! 24!/****************************************************************************** 25! * 26! * FreeSurface utilities 27! ! TODO: Not quite finished 28! ! TODO1: outdated should be removed... (12.12.2003, Juha) 29! * 30! ****************************************************************************** 31! * 32! * Authors: Juha Ruokolainen 33! * Email: Juha.Ruokolainen@csc.fi 34! * Web: http://www.csc.fi/elmer 35! * Address: CSC - IT Center for Science Ltd. 36! * Keilaranta 14 37! * 02101 Espoo, Finland 38! * 39! * Original Date: 02 Jun 1997 40! * 41! ****************************************************************************/ 42 43!> Internal free surface utilities. 44!> \deprecated 45!> \ingroup ElmerLib 46!> \{ 47 48MODULE FreeSurface 49 50 USE DirectSolve 51 USE IterSolve 52 USE ElementUtils 53 54 IMPLICIT NONE 55 56CONTAINS 57!------------------------------------------------------------------------------- 58! 59! 60 SUBROUTINE MeanCurvature( Model ) 61!------------------------------------------------------------------------------- 62 TYPE(Model_t) :: Model 63!------------------------------------------------------------------------------- 64 INTEGER :: i,j,k,n,t,BC,NodeIndexes(16) 65 INTEGER, POINTER :: Reorder(:),Visited(:) 66 LOGICAL :: L 67 68 REAL(KIND=dp) :: ddxddu,ddyddu,ddzddu,ddxdudv,ddydudv,ddzdudv, & 69 ddxddv,ddyddv,ddzddv,Auu,Auv,Avv,Buu,Buv,Bvv,detA,u,v,x,y,z 70 71 REAL(KIND=dp), ALLOCATABLE :: dxdu(:),dydu(:),dzdu(:) 72 REAL(KIND=dp), ALLOCATABLE :: dxdv(:),dydv(:),dzdv(:) 73 74 REAL(KIND=dp), TARGET :: nx(16),ny(16),nz(16),Nrm(3) 75 REAL(KIND=dp), POINTER :: Curvature(:) 76 77 REAL(KIND=dp) :: Basis(16),dBasisdx(16,3),ddBasisddx(16,3,3) 78 TYPE(Nodes_t) :: Nodes 79 TYPE(Element_t), POINTER :: Boundary 80 81 LOGICAL :: stat 82 83 REAL(KIND=dp) :: Metric(3,3),SqrtMetric,Symbols(3,3,3),dSymbols(3,3,3,3) 84!------------------------------------------------------------------------------- 85 ALLOCATE( Visited(Model % NumberOfNodes) ) 86 87 IF ( .NOT.ASSOCIATED( Model % FreeSurfaceNodes) ) THEN 88 ALLOCATE( Model % FreeSurfaceNodes( Model % NumberOfNodes ) ) 89 Model % FreeSurfaceNodes = 0 90 END IF 91 Reorder => Model % FreeSurfaceNodes 92 93 ! Count nodes on free boundaries, and make a permutation table for nodes, 94 ! should perhaps be saved instead of recomputed every time: 95 !------------------------------------------------------------------------ 96 Visited = 0 97 Reorder = 0 98 k = 0 99 DO BC = 1,Model % NumberOfBCs 100 IF (.NOT. ListGetLogical(Model % BCs(BC) % Values, 'Free Surface',L)) CYCLE 101 102 DO t=Model % NumberOfBulkElements+1,Model % NumberOfBulkElements + & 103 Model % NumberOfBoundaryElements 104 105 Boundary => Model % Elements(t) 106 IF ( Boundary % BoundaryInfo % Constraint /= BC ) CYCLE 107 108 n = Boundary % Type % NumberOfNodes 109 DO i=1,n 110 j = Boundary % NodeIndexes(i) 111 IF ( Visited(j) == 0 ) THEN 112 k = k + 1 113 Reorder(j) = k 114 END IF 115 Visited(j) = Visited(j) + 1 116 END DO 117 END DO 118 END DO 119!------------------------------------------------------------------------------- 120! If no free surfaces, return 121!------------------------------------------------------------------------------- 122 IF ( k == 0 ) THEN 123 DEALLOCATE( Visited ) 124 RETURN 125 END IF 126!------------------------------------------------------------------------------- 127! Allocate some memories... 128!------------------------------------------------------------------------------- 129 IF ( .NOT.ASSOCIATED( Model % BoundaryCurvatures ) ) THEN 130 ALLOCATE( Model % BoundaryCurvatures(k) ) 131 END IF 132 Curvature => Model % BoundaryCurvatures 133 Curvature = 0.0D0 134 135 ALLOCATE( dxdu(k),dydu(k),dzdu(k),dxdv(k),dydv(k),dzdv(k) ) 136 dxdu = 0.0D0 137 dydu = 0.0D0 138 dzdu = 0.0D0 139 dxdv = 0.0D0 140 dydv = 0.0D0 141 dzdv = 0.0D0 142!------------------------------------------------------------------------------- 143! Compute sum of first derivatives on nodes of the boundaries 144!------------------------------------------------------------------------------- 145 DO BC = 1,Model % NumberOfBCs 146 147 IF (.NOT. ListGetLogical(Model % BCs(BC) % Values, 'Free Surface',L)) CYCLE 148 149 DO t=Model % NumberOfBulkElements+1,Model % NumberOfBulkElements + & 150 Model % NumberOfBoundaryElements 151 152 Boundary => Model % Elements(t) 153 IF ( Boundary % BoundaryInfo % Constraint /= BC ) CYCLE 154 155 n = Boundary % Type % NumberOfNodes 156 nx(1:n) = Model % Nodes % x(Boundary % NodeIndexes) 157 ny(1:n) = Model % Nodes % y(Boundary % NodeIndexes) 158 nz(1:n) = Model % Nodes % z(Boundary % NodeIndexes) 159#if 1 160nodes % x => nx(1:n) 161nodes % y => ny(1:n) 162nodes % z => nz(1:n) 163#endif 164 165 DO i=1,n 166 j = Reorder(Boundary % NodeIndexes(i)) 167 168 IF ( Boundary % Type % Dimension == 1 ) THEN 169 u = Boundary % Type % NodeU(i) 170 171#if 0 172 dxdu(j) = dxdu(j) + FirstDerivative1D( Boundary,nx(1:n),u ) 173 dydu(j) = dydu(j) + FirstDerivative1D( Boundary,ny(1:n),u ) 174#else 175stat = ElementInfo( Boundary, Nodes, u, 0.0d0, 0.0d0, detA, & 176 Basis, dBasisdx ) 177dydu(j) = dydu(j) + SUM( dBasisdx(1:n,1) * ny(1:n) ) 178#endif 179 ELSE 180 u = Boundary % Type % NodeU(i) 181 v = Boundary % Type % NodeV(i) 182 183 dxdu(j) = dxdu(j) + FirstDerivativeInU2D( Boundary,nx(1:n),u,v ) 184 dydu(j) = dydu(j) + FirstDerivativeInU2D( Boundary,ny(1:n),u,v ) 185 dzdu(j) = dzdu(j) + FirstDerivativeInU2D( Boundary,nz(1:n),u,v ) 186 187 dxdv(j) = dxdv(j) + FirstDerivativeInV2D( Boundary,nx(1:n),u,v ) 188 dydv(j) = dydv(j) + FirstDerivativeInV2D( Boundary,ny(1:n),u,v ) 189 dzdv(j) = dzdv(j) + FirstDerivativeInV2D( Boundary,nz(1:n),u,v ) 190 END IF 191 END DO 192 END DO 193 END DO 194 195!------------------------------------------------------------------------------- 196! Average of the derivatives at nodes... 197!------------------------------------------------------------------------------- 198 DO BC=1,Model % NumberOfBCs 199 IF (.NOT. ListGetLogical(Model % BCs(BC) % Values, 'Free Surface',L)) CYCLE 200 201 DO t=Model % NumberOfBulkElements+1,Model % NumberOfBulkElements + & 202 Model % NumberOfBoundaryElements 203 204 Boundary => Model % Elements(t) 205 IF ( Boundary % BoundaryInfo % Constraint /= BC ) CYCLE 206 207 n = Boundary % Type % NumberOfNodes 208 DO i=1,n 209 j = Boundary % NodeIndexes(i) 210 IF ( Visited(j) > 0 ) THEN 211 dxdu(Reorder(j)) = dxdu(Reorder(j)) / Visited(j) 212 dydu(Reorder(j)) = dydu(Reorder(j)) / Visited(j) 213 dzdu(Reorder(j)) = dzdu(Reorder(j)) / Visited(j) 214 dxdv(Reorder(j)) = dxdu(Reorder(j)) / Visited(j) 215 dydv(Reorder(j)) = dydv(Reorder(j)) / Visited(j) 216 dzdv(Reorder(j)) = dzdv(Reorder(j)) / Visited(j) 217 Visited(j) = -Visited(j) 218 END IF 219 END DO 220 END DO 221 END DO 222 223 Visited = ABS( Visited ) 224 225!------------------------------------------------------------------------------- 226! The curvature computation begins 227!------------------------------------------------------------------------------- 228 DO BC=1,Model % NumberOfBCs 229 230 IF (.NOT. ListGetLogical(Model % BCs(BC) % Values, 'Free Surface',L)) CYCLE 231 232 DO t=Model % NumberOfBulkElements+1,Model % NumberOfBulkElements + & 233 Model % NumberOfBoundaryElements 234 235 Boundary => Model % Elements(t) 236 237 IF ( Boundary % BoundaryInfo % Constraint /= BC ) CYCLE 238 239 n = Boundary % Type % NumberOfNodes 240 NodeIndexes(1:n) = Boundary % NodeIndexes 241!------------------------------------------------------------------------------- 242! Go through element nodal points 243!------------------------------------------------------------------------------- 244#if 1 245nx(1:n) = Model % Nodes % x(NodeIndexes(1:n)) 246ny(1:n) = Model % Nodes % y(NodeIndexes(1:n)) 247nz(1:n) = Model % Nodes % z(NodeIndexes(1:n)) 248nodes % x => nx(1:n) 249nodes % y => ny(1:n) 250nodes % z => nz(1:n) 251#endif 252 DO i=1,n 253 j = Reorder( NodeIndexes(i) ) 254 255 u = Boundary % Type % NodeU(i) 256 v = 0.0D0 257 IF ( Boundary % Type % Dimension > 1 ) v = Boundary % Type % NodeV(i) 258 259 SqrtMetric = 1.0D0 260 x = Model % Nodes % x(NodeIndexes(i)) 261 y = Model % Nodes % y(NodeIndexes(i)) 262 z = Model % Nodes % z(NodeIndexes(i)) 263 264 IF (CurrentCoordinateSystem() /= Cartesian ) THEN 265 CALL CoordinateSystemInfo( Metric,SqrtMetric,Symbols,dSymbols,X,Y,Z ) 266 END IF 267!------------------------------------------------------------------------------- 268! 2D case, compute the curvature 269!------------------------------------------------------------------------------- 270 IF ( Boundary % Type % Dimension == 1 ) THEN 271#if 0 272!------------------------------------------------------------------------------- 273! Second partial derivatives of the space coordinates with respect to 274! curve coordinate 275!------------------------------------------------------------------------------- 276 ddxddu = FirstDerivative1D( Boundary,dxdu(Reorder(NodeIndexes(1:n))),u ) 277 ddyddu = FirstDerivative1D( Boundary,dydu(Reorder(NodeIndexes(1:n))),u ) 278!------------------------------------------------------------------------------- 279! curve 'metric' 280!------------------------------------------------------------------------------- 281 Auu = dxdu(j)*dxdu(j) + dydu(j)*dydu(j) 282 detA = 1.0d0 / SQRT(Auu) 283 284 Nrm(1) = -dydu(j) 285 Nrm(2) = dxdu(j) 286 Nrm(3) = 0.0D0 287 Nrm = Nrm / SQRT( SUM( Nrm**2 ) ) 288 CALL CheckNormalDirection( Boundary,Nrm,x,y,z ) 289!------------------------------------------------------------------------------- 290! and finally the curvature 291!------------------------------------------------------------------------------- 292 Curvature(j) = Curvature(j) + & 293 0.5d0 * Auu * ( ddxddu*Nrm(1) + ddyddu*Nrm(2) ) 294#else 295stat = ElementInfo( Boundary, Nodes, u, 0.0d0, 0.0d0, detA, & 296 Basis, dBasisdx, ddBasisddx, .TRUE. ) 297ddyddu = SUM( dBasisdx(1:n,1) * dydu(Reorder(NodeIndexes(1:n))) ) 298Curvature(j) = Curvature(j) + ( y * ddyddu - (1+dydu(j)**2) ) / ( y * ( 1+dydu(j)**2 )**(3.0d0/2.0d0) ) 299#endif 300 301 302 ELSE 303!------------------------------------------------------------------------------- 304! 3D case, compute the curvature 305!------------------------------------------------------------------------------- 306 u = Boundary % Type % NodeU(i) 307 v = Boundary % Type % NodeV(i) 308!------------------------------------------------------------------------------- 309! Second partial derivatives of the space coordinates with respect to 310! surface coordinates 311!------------------------------------------------------------------------------- 312 ddxddu = FirstDerivativeInU2D( Boundary, dxdu(NodeIndexes(1:n)),u,v ) 313 ddyddu = FirstDerivativeInU2D( Boundary, dydu(NodeIndexes(1:n)),u,v ) 314 ddzddu = FirstDerivativeInU2D( Boundary, dzdu(NodeIndexes(1:n)),u,v ) 315 316 ddxdudv = FirstDerivativeInU2D( Boundary, dxdv(NodeIndexes(1:n)),u,v ) 317 ddydudv = FirstDerivativeInU2D( Boundary, dydv(NodeIndexes(1:n)),u,v ) 318 ddzdudv = FirstDerivativeInU2D( Boundary, dzdv(NodeIndexes(1:n)),u,v ) 319 320 ddxddv = FirstDerivativeInV2D( Boundary, dxdv(NodeIndexes(1:n)),u,v ) 321 ddyddv = FirstDerivativeInV2D( Boundary, dydv(NodeIndexes(1:n)),u,v ) 322 ddzddv = FirstDerivativeInV2D( Boundary, dzdv(NodeIndexes(1:n)),u,v ) 323!------------------------------------------------------------------------------- 324! Surface metric 325!------------------------------------------------------------------------------- 326 Auu = dxdu(k)*dxdu(k) + dydu(k)*dydu(k) + dzdu(k)*dzdu(k) 327 Auv = dxdu(k)*dxdv(k) + dydu(k)*dydv(k) + dzdu(k)*dzdv(k) 328 Avv = dxdv(k)*dxdv(k) + dydv(k)*dydv(k) + dzdv(k)*dzdv(k) 329 330 detA = 1.0D0 / SQRT(Auu*Avv - Auv*Auv) 331!------------------------------------------------------------------------------- 332! Change metric to contravariant form 333!------------------------------------------------------------------------------- 334 u = Auu 335 Auu = Avv * detA 336 Auv = -Auv * detA 337 Avv = u * detA 338!------------------------------------------------------------------------------- 339! Normal vector to surface 340!------------------------------------------------------------------------------- 341 Nrm(1) = (dydu(k) * dzdv(k) - dydv(k) * dzdu(k)) * detA 342 Nrm(2) = (dxdv(k) * dzdu(k) - dxdu(k) * dzdv(k)) * detA 343 Nrm(3) = (dxdu(k) * dydv(k) - dxdv(k) * dydu(k)) * detA 344 Nrm = Nrm / SQRT( SUM( Nrm**2 ) ) 345 346 CALL CheckNormalDirection( Boundary,Nrm,x,y,z ) 347!------------------------------------------------------------------------------- 348! The second fundamental form of the surface 349!------------------------------------------------------------------------------- 350 Buu = ddxddu * Nrm(1) + ddyddu * Nrm(2) + ddzddu * Nrm(3) 351 Buv = ddxdudv * Nrm(1) + ddydudv * Nrm(2) + ddzdudv * Nrm(3) 352 Bvv = ddxddv * Nrm(1) + ddyddv * Nrm(2) + ddzddv * Nrm(3) 353!------------------------------------------------------------------------------- 354! And finally, the curvature 355!------------------------------------------------------------------------------- 356 Curvature(j) = Curvature(j) + (Auu*Buu + 2*Auv*Buv + Avv*Bvv) 357 END IF 358 END DO 359 END DO 360 END DO 361visited=abs(visited) 362!------------------------------------------------------------------------------- 363! Average to nodes 364!------------------------------------------------------------------------------- 365 DO BC=1,Model % NumberOfBCs 366 IF (.NOT. ListGetLogical(Model % BCs(BC) % Values, 'Free Surface',L)) CYCLE 367 368 DO t=Model % NumberOfBulkElements+1,Model % NumberOfBulkElements + & 369 Model % NumberOfBoundaryElements 370 371 Boundary => Model % Elements(t) 372 IF ( Boundary % BoundaryInfo % Constraint /= BC ) CYCLE 373 374 n = Boundary % Type % NumberOfNodes 375 DO i=1,n 376 j = Boundary % NodeIndexes(i) 377 IF ( Visited(j) > 0 ) THEN 378 Curvature(Reorder(j)) = Curvature(Reorder(j)) / Visited(j) 379write(*,*) j,curvature(Reorder(j)) 380 Visited(j) = -Visited(j) 381 END IF 382 END DO 383 END DO 384 END DO 385!------------------------------------------------------------------------------- 386! We are done here. 387!------------------------------------------------------------------------------- 388 DEALLOCATE( Visited,dxdu,dydu,dzdu,dxdv,dydv,dzdv ) 389print*,'----------------------' 390!------------------------------------------------------------------------------- 391 END SUBROUTINE MeanCurvature 392!------------------------------------------------------------------------------- 393 394 395 396 397!------------------------------------------------------------------------------- 398 SUBROUTINE MoveBoundary( Model,Relax ) 399!------------------------------------------------------------------------------- 400 USE DefUtils 401 TYPE(Model_t) :: Model 402 REAL(KIND=dp) :: Relax 403!------------------------------------------------------------------------------- 404 INTEGER :: i,ii,j,k,m,n,t,Which, FIter 405 LOGICAL, ALLOCATABLE :: Visited(:),Turned(:) 406 LOGICAL :: L 407 408 REAL(KIND=dp) :: Auu,Auv,Avv,Buu,Buv,Bvv,detA,u,v,x,y,z,S,R,Ux,Uy,Uz 409 410 REAL(KIND=dp) :: dxdu,dydu,dzdu, FEPS 411 REAL(KIND=dp) :: dxdv,dydv,dzdv,x1,x2,y1,y2,z1,z2 412 413 REAL(KIND=dp), TARGET :: N1,N2,N3,Nrm(3) 414 REAL(KIND=dp), POINTER :: Curvature(:) 415 416 TYPE(Element_t), POINTER :: Boundary, Element 417 TYPE(Nodes_t) :: BoundaryNodes, ElementNodes 418 419 REAL(KIND=dp), ALLOCATABLE :: XCoord(:),YCoord(:),ZCoord(:), AvarageNormal(:,:) 420 LOGICAL :: XMoved, YMoved, ZMoved 421 422 TYPE( Variable_t), POINTER :: Velocity1,Velocity2,Velocity3 423 424 TYPE(Mesh_t), POINTER :: Mesh 425 TYPE(ValueList_t), POINTER :: BC 426 427 REAL(KIND=dp) :: SqrtMetric, dLBasisdx(16,2),NodalBasis(16) 428 429 SAVE BoundaryNodes, ElementNodes 430!------------------------------------------------------------------------------- 431 432 FEPS = ListGetConstReal( Model % Solver % Values, & 433 'Free Surface Convergence Tolerance', L ) 434 IF ( .NOT. L ) FEPS = 1.0d-12 435 436 FIter = ListGetInteger( Model % Solver % Values, & 437 'Free Surface Max Iterations', L ) 438 IF ( .NOT. L ) FIter = 100 439 440 Velocity1 => VariableGet( Model % Variables, 'Velocity 1' ) 441 Velocity2 => VariableGet( Model % Variables, 'Velocity 2' ) 442 Velocity3 => VariableGet( Model % Variables, 'Velocity 3' ) 443!------------------------------------------------------------------------------- 444 ALLOCATE( Visited(Model % NumberOfNodes),Turned(Model % NumberOfNodes) ) 445 Visited = .FALSE. 446 Turned = .FALSE. 447!------------------------------------------------------------------------------- 448 449 ALLOCATE( XCoord(Model % NumberOfNodes) ) 450 XMoved = .FALSE. 451 XCoord = Model % Nodes % x 452 453 ALLOCATE( YCoord(Model % NumberOfNodes) ) 454 YMoved = .FALSE. 455 YCoord = Model % Nodes % y 456 457 ALLOCATE( ZCoord(Model % NumberOfNodes) ) 458 ZMoved = .FALSE. 459 ZCoord = Model % Nodes % z 460 461 ALLOCATE( AvarageNormal(Model % NumberOfBCs,3) ) 462 AvarageNormal = 0 463 464 Ux = 0; Uy = 0; Uz = 0; 465!------------------------------------------------------------------------------- 466! Check normal direction first in a separate loop, cause within the node 467! moving loop, the parent elements might be a little funny...! 468!------------------------------------------------------------------------------- 469 Mesh => Model % Solver % Mesh 470 DO t = 1,Mesh % NumberOfBoundaryElements 471 Boundary => GetBoundaryElement(t) 472 IF ( .NOT. ActiveBoundaryElement() ) CYCLE 473 474 BC => GetBC() 475 IF ( .NOT. ASSOCIATED(BC) ) CYCLE 476 IF ( .NOT.GetLogical(BC, 'Free Surface',L) ) CYCLE 477 478 n = GetElementNOFNodes() 479 CALL GetElementNodes( BoundaryNodes ) 480!------------------------------------------------------------------------------- 481! Go through boundary element nodes 482!------------------------------------------------------------------------------- 483 DO i=1,n 484 k = Boundary % NodeIndexes(i) 485 486 IF ( Boundary % Type % Dimension == 1 ) THEN 487 IF ( i==1 ) THEN 488 j=2 489 ELSE 490 j=1 491 END IF 492 j = Boundary % NodeIndexes(j) 493 x1 = XCoord(j) - XCoord(k) 494 y1 = YCoord(j) - YCoord(k) 495 496 Ux = Velocity1 % Values( Velocity1 % Perm(k) ) 497 Uy = Velocity2 % Values( Velocity2 % Perm(k) ) 498 499 IF ( Ux*x1 + Uy*y1 > 0 ) CYCLE 500 END IF 501 502!------------------------------------------------------------------------------- 503! Should not move the same node twice,so check if already done 504!------------------------------------------------------------------------------- 505 IF ( .NOT.Visited(k) ) THEN 506 Visited(k) = .TRUE. 507!------------------------------------------------------------------------------- 508! 2D case, compute normal 509!------------------------------------------------------------------------------- 510 IF ( Boundary % Type % Dimension == 1 ) THEN 511 u = Boundary % Type % NodeU(i) 512!------------------------------------------------------------------------------ 513! Basis function derivatives with respect to local coordinates 514!------------------------------------------------------------------------------ 515 NodalBasis(1:n) = 0.0d0 516 DO m=1,n 517 NodalBasis(m) = 1.0d0 518 dLBasisdx(m,1) = FirstDerivative1D( Boundary,NodalBasis,u ) 519 NodalBasis(m) = 0.0d0 520 END DO 521 522 Nrm(1) = -SUM( BoundaryNodes % y(1:n)*dLBasisdx(1:n,1) ) 523 Nrm(2) = SUM( BoundaryNodes % x(1:n)*dLBasisdx(1:n,1) ) 524 Nrm(3) = 0.0D0 525 ELSE 526!------------------------------------------------------------------------------- 527! 3D case, compute normal 528!------------------------------------------------------------------------------- 529 u = Boundary % Type % NodeU(i) 530 v = Boundary % Type % NodeV(i) 531 532 NodalBasis(1:n) = 0.0D0 533 DO m=1,N 534 NodalBasis(m) = 1.0D0 535 dLBasisdx(m,1) = FirstDerivativeInU2D( Boundary,NodalBasis,u,v ) 536 dLBasisdx(m,2) = FirstDerivativeInV2D( Boundary,NodalBasis,u,v ) 537 NodalBasis(m) = 0.0D0 538 END DO 539 540 dxdu = SUM( BoundaryNodes % x(1:n) * dLBasisdx(1:n,1) ) 541 dydu = SUM( BoundaryNodes % y(1:n) * dLBasisdx(1:n,1) ) 542 dzdu = SUM( BoundaryNodes % z(1:n) * dLBasisdx(1:n,1) ) 543 544 dxdv = SUM( BoundaryNodes % x(1:n) * dLBasisdx(1:n,2) ) 545 dydv = SUM( BoundaryNodes % y(1:n) * dLBasisdx(1:n,2) ) 546 dzdv = SUM( BoundaryNodes % z(1:n) * dLBasisdx(1:n,2) ) 547 548 Nrm(1) = dydu * dzdv - dydv * dzdu 549 Nrm(2) = dxdv * dzdu - dxdu * dzdv 550 Nrm(3) = dxdu * dydv - dxdv * dydu 551 END IF 552!------------------------------------------------------------------------------- 553! Turn the normal to point outwards, or towards less dense material 554!------------------------------------------------------------------------------- 555 x = Model % Nodes % x(k) 556 y = Model % Nodes % y(k) 557 z = Model % Nodes % z(k) 558 CALL CheckNormalDirection( Boundary,Nrm,x,y,z,Turned(k) ) 559 560 m = Boundary % BoundaryInfo % Constraint 561 R = SQRT(SUM(Nrm**2)) 562 AvarageNormal(m,1:3) = AvarageNormal(m,1:3) + ABS(Nrm)/R 563 END IF 564 END DO 565 END DO 566 567!------------------------------------------------------------------------------- 568! Iterate until convergence 569!------------------------------------------------------------------------------- 570 DO ii=1,FIter 571 S = 0.0d0 572 Visited = .FALSE. 573 574 DO t=1,Mesh % NumberOfBoundaryElements 575 Boundary => GetBoundaryElement(t) 576 IF ( .NOT. ActiveBoundaryElement() ) CYCLE 577 578 BC => GetBC() 579 IF (.NOT.GetLogical(BC, 'Free Surface',L)) CYCLE 580 581 i = CoordinateSystemDimension() 582 Which = ListGetInteger( BC, 'Free Coordinate', L,minv=1, maxv=i ) 583 IF ( .NOT. L ) THEN 584 i = Boundary % BoundaryInfo % Constraint 585 Nrm = AvarageNormal(i,:) 586 Which = 1 587 DO i=2,3 588 IF ( Nrm(i) > Nrm(Which) ) Which=i 589 END DO 590 END IF 591 592 n = GetElementNOFNOdes() 593 CALL GetElementNodes( BoundaryNodes ) 594 595 DO i=1,n 596 k = Boundary % NodeIndexes(i) 597 598 IF ( Boundary % Type % Dimension == 1 ) THEN 599 IF ( i==1 ) THEN 600 j=2 601 ELSE 602 j=1 603 END IF 604 j = Boundary % NodeIndexes(j) 605 x1 = XCoord(j) - XCoord(k) 606 y1 = YCoord(j) - YCoord(k) 607 608 Ux = Velocity1 % Values( Velocity1 % Perm(k) ) 609 Uy = Velocity2 % Values( Velocity2 % Perm(k) ) 610 611 IF ( Ux*x1 + Uy*y1 > 0 ) CYCLE 612 END IF 613 614!------------------------------------------------------------------------------- 615! Should not move the same node twice,so check if already done 616!------------------------------------------------------------------------------- 617 IF ( .NOT.Visited(k) ) THEN 618 Visited(k) = .TRUE. 619!------------------------------------------------------------------------------- 620! 2D case, compute normal 621!------------------------------------------------------------------------------- 622 IF ( Boundary % Type % Dimension == 1 ) THEN 623 u = Boundary % Type % NodeU(i) 624!------------------------------------------------------------------------------ 625! Basis function derivatives with respect to local coordinates 626!------------------------------------------------------------------------------ 627 NodalBasis(1:n) = 0.0D0 628 DO m=1,n 629 NodalBasis(m) = 1.0D0 630 dLBasisdx(m,1) = FirstDerivative1D( Boundary,NodalBasis,u ) 631 NodalBasis(m) = 0.0D0 632 END DO 633 634 Nrm(1) = -SUM( BoundaryNodes % y(1:n)*dLBasisdx(1:n,1) ) 635 Nrm(2) = SUM( BoundaryNodes % x(1:n)*dLBasisdx(1:n,1) ) 636 Nrm(3) = 0.0D0 637 ELSE 638!---------------------------------------------------------------------------- 639! 3D case, compute normal 640!---------------------------------------------------------------------------- 641 u = Boundary % Type % NodeU(i) 642 v = Boundary % Type % NodeV(i) 643 644 NodalBasis(1:n) = 0.0D0 645 DO j=1,n 646 NodalBasis(j) = 1.0D0 647 dLBasisdx(j,1) = FirstDerivativeInU2D( Boundary,NodalBasis,u,v ) 648 dLBasisdx(j,2) = FirstDerivativeInV2D( Boundary,NodalBasis,u,v ) 649 NodalBasis(j) = 0.0D0 650 END DO 651 652 dxdu = SUM( BoundaryNodes % x(1:n) * dLBasisdx(1:n,1) ) 653 dydu = SUM( BoundaryNodes % y(1:n) * dLBasisdx(1:n,1) ) 654 dzdu = SUM( BoundaryNodes % z(1:n) * dLBasisdx(1:n,1) ) 655 656 dxdv = SUM( BoundaryNodes % x(1:n) * dLBasisdx(1:n,2) ) 657 dydv = SUM( BoundaryNodes % y(1:n) * dLBasisdx(1:n,2) ) 658 dzdv = SUM( BoundaryNodes % z(1:n) * dLBasisdx(1:n,2) ) 659 660 Nrm(1) = dydu * dzdv - dydv * dzdu 661 Nrm(2) = dxdv * dzdu - dxdu * dzdv 662 Nrm(3) = dxdu * dydv - dxdv * dydu 663 END IF 664!------------------------------------------------------------------------------- 665! Turn the normal to point outwards, or towards less dense material 666!------------------------------------------------------------------------------- 667 IF ( Turned(k) ) Nrm = -Nrm 668!------------------------------------------------------------------------------- 669! Now then, lets move the node so that u.n will be reduced. 670!------------------------------------------------------------------------------- 671 IF ( Boundary % Type % Dimension == 1 ) THEN 672!------------------------------------------------------------------------------- 673! TODO :: This won't handle the three node line 674! 2D case, move the nodes.. 675!------------------------------------------------------------------------------- 676 R = Ux * Nrm(1) + Uy * Nrm(2) 677 IF ( Which == 2 ) THEN 678 IF ( ABS(Ux) > AEPS ) THEN 679 IF ( Turned(k) ) THEN 680 Model % Nodes % y(k) = Model % Nodes % y(k) - & 681 R / ( Ux*dLBasisdx(i,1) ) 682 ELSE 683 Model % Nodes % y(k) = Model % Nodes % y(k) + & 684 R / ( Ux*dLBasisdx(i,1) ) 685 END IF 686 YMoved = .TRUE. 687 END IF 688 ELSE 689 IF ( ABS(Uy) > AEPS ) THEN 690 IF ( Turned(k) ) THEN 691 Model % Nodes % x(k) = Model % Nodes % x(k) + & 692 R / ( Uy*dLBasisdx(i,1) ) 693 ELSE 694 Model % Nodes % x(k) = Model % Nodes % x(k) - & 695 R / ( Uy*dLBasisdx(i,1) ) 696 END IF 697 XMoved = .TRUE. 698 END IF 699 END IF 700 ELSE 701!------------------------------------------------------------------------------- 702! 3D case, move the nodes.. 703! TODO :: This is just guesswork, no testing done... 704!---------------------------------------------------------------------------- 705 Ux = Velocity1 % Values( Velocity1 % Perm(k) ) 706 Uy = Velocity2 % Values( Velocity2 % Perm(k) ) 707 Uz = Velocity3 % Values( Velocity3 % Perm(k) ) 708 709 R = Ux * Nrm(1) + Uy * Nrm(2) + Uz * Nrm(3) 710 IF ( Which == 1 ) THEN 711 IF ( ABS(Uy) > AEPS .OR. ABS(Uz) > AEPS ) THEN 712 Model % Nodes % x(k) = Model % Nodes % x(k) + R / & 713 ( (dzdu*dLBasisdx(i,2) - dzdv*dLBasisdx(i,1))*Uy + & 714 (dydv*dLBasisdx(i,1) - dydu*dLBasisdx(i,2))*Uz ) 715 XMoved = .TRUE. 716 END IF 717 ELSE IF ( Which == 2 ) THEN 718 IF ( ABS(Ux) > AEPS .OR. ABS(Uz) > AEPS ) THEN 719 Model % Nodes % y(k) = Model % Nodes % y(k) + R / & 720 ( (dzdv*dLBasisdx(i,1) - dzdu*dLBasisdx(i,2))*Ux + & 721 (dxdu*dLBasisdx(i,2) - dxdv*dLBasisdx(i,1))*Uz ) 722 YMoved = .TRUE. 723 END IF 724 ELSE 725 IF ( ABS(Ux) > AEPS .OR. ABS(Uy) > AEPS ) THEN 726 Model % Nodes % z(k) = Model % Nodes % z(k) + R / & 727 ( (dydu*dLBasisdx(i,2) - dydv*dLBasisdx(i,1))*Ux + & 728 (dxdv*dLBasisdx(i,1) - dxdu*dLBasisdx(i,2))*Uy ) 729 ZMoved = .TRUE. 730 END IF 731 END IF 732 END IF 733!------------------------------------------------------------------------------- 734 S = S + ( (Ux*Nrm(1)+Uy*Nrm(2)+Uz*Nrm(3))/SQRT(SUM(Nrm**2)) )**2 735!------------------------------------------------------------------------------- 736 END IF 737 END DO 738 END DO 739!------------------------------------------------------------------------------- 740 S = SQRT(S) 741 WRITE( Message, '(a,i4,a,e21.12)' ) 'Iter: ', ii, ' Free surface Residual: ',S 742 CALL Info( 'FreeSurface', Message, Level=4 ) 743 IF ( S < FEPS ) EXIT 744 END DO 745!------------------------------------------------------------------------------- 746 IF ( XMoved ) THEN 747 CALL PoissonSolve( Model,XCoord,YCoord,ZCoord,Model % Nodes % x,1 ) 748 IF ( Relax /= 1.0d0 ) & 749 Model % Nodes % x = Relax * Model % Nodes % x + (1-Relax)*XCoord 750 END IF 751 752 IF ( YMoved ) THEN 753 CALL PoissonSolve( Model,XCoord,YCoord,ZCoord,Model % Nodes % y,2 ) 754 IF ( Relax /= 1.0d0 ) & 755 Model % Nodes % y = Relax * Model % Nodes % y + (1-Relax)*YCoord 756 END IF 757 758 IF ( ZMoved ) THEN 759 CALL PoissonSolve( Model,XCoord,YCoord,ZCoord,Model % Nodes % z,3 ) 760 IF ( Relax /= 1.0d0 ) & 761 Model % Nodes % z = Relax * Model % Nodes % z + (1-Relax)*ZCoord 762 END IF 763 764 DEALLOCATE( Visited,Turned,XCoord,YCoord,ZCoord,AvarageNormal ) 765 766 ALLOCATE( ElementNodes % x(Model % MaxElementNodes) ) 767 ALLOCATE( ElementNodes % y(Model % MaxElementNodes) ) 768 ALLOCATE( ElementNodes % z(Model % MaxElementNodes) ) 769 770 DO i=1,Model % NumberOfBulkElements 771 Element => Model % Elements(i) 772 n = Element % Type % NumberOfNodes 773 ElementNodes % x(1:n) = Model % Nodes % x(Element % NodeIndexes) 774 ElementNodes % y(1:n) = Model % Nodes % y(Element % NodeIndexes) 775 ElementNodes % z(1:n) = Model % Nodes % z(Element % NodeIndexes) 776 CALL StabParam( Element, ElementNodes, n, & 777 Element % StabilizationMK, Element % hK ) 778 END DO 779 780 DEALLOCATE( ElementNodes % x, ElementNodes % y, ElementNodes % z) 781!------------------------------------------------------------------------------- 782 END SUBROUTINE MoveBoundary 783!------------------------------------------------------------------------------- 784 785 786 787!------------------------------------------------------------------------------- 788 SUBROUTINE PoissonSolve( Model,NX,NY,NZ,Solution,Moved ) 789!------------------------------------------------------------------------------- 790 TYPE(Model_t) :: Model 791 INTEGER :: Moved 792 REAL(KIND=dp) :: NX(:),NY(:),NZ(:),Solution(:) 793!------------------------------------------------------------------------------- 794 795 TYPE(Matrix_t), POINTER :: CMatrix 796 INTEGER, POINTER :: CPerm(:) 797 798 LOGICAL :: FirstTime = .TRUE. 799 REAL(KIND=dp), ALLOCATABLE :: ForceVector(:) 800 801 REAL(KIND=dp) :: Basis(16),dBasisdx(16,3),ddBasisddx(16,3,3) 802 REAL(KIND=dp) :: SqrtElementMetric,u,v,w,s,A 803 LOGICAL :: Stat 804 805 INTEGER :: i,j,k,n,q,p,t,dim 806 807 TYPE(Solver_t), POINTER :: Solver 808 TYPE(Nodes_t) :: Nodes 809 TYPE(Element_t), POINTER :: Element 810 811 INTEGER :: N_Integ 812 813 REAL(KIND=dp), POINTER :: U_Integ(:),V_Integ(:),W_Integ(:),S_Integ(:), & 814 LocalMatrix(:,:) 815 816 TYPE(GaussIntegrationPoints_t), TARGET :: IntegStuff 817 818 SAVE CMatrix,CPerm,FirstTime,ForceVector,Solver,Nodes,LocalMatrix 819!------------------------------------------------------------------------------- 820 821 IF ( FirstTime ) THEN 822 FirstTime = .FALSE. 823 824 ALLOCATE( Solver ) 825 Solver % Values => ListAllocate() 826 Solver % Mesh => CurrentModel % Mesh 827 828 CALL ListAddString( Solver % Values, & 829 'Linear System Iterative Method', 'CGS' ) 830 CALL ListAddInteger( Solver % Values, & 831 'Linear System Max Iterations', 500 ) 832 CALL ListAddConstReal( Solver % Values, & 833 'Linear System Convergence Tolerance', 1.0D-9 ) 834 CALL ListAddString( Solver % Values, & 835 'Linear System Preconditioning', 'ILU0' ) 836 CALL ListAddInteger( Solver % Values, & 837 'Linear System Residual Output', 1 ) 838 839 ALLOCATE( ForceVector( Model % NumberOfNodes ), & 840 CPerm( Model % NumberOfNodes) ) 841 842 CMatrix => CreateMatrix( CurrentModel,Solver,Solver % Mesh,CPerm,1,MATRIX_CRS,.FALSE. ) 843 844 ALLOCATE( LocalMatrix( Model % MaxElementNodes,Model % MaxElementNodes ) ) 845 ALLOCATE( Nodes % x(Model % MaxElementNodes) ) 846 ALLOCATE( Nodes % y(Model % MaxElementNodes) ) 847 ALLOCATE( Nodes % z(Model % MaxElementNodes) ) 848 849 Solver % TimeOrder = 0 850 END IF 851 852!------------------------------------------------------------------------------ 853 854 CALL CRS_ZeroMatrix( CMatrix ) 855 ForceVector = 0.0D0 856 857!------------------------------------------------------------------------------ 858 DO i=1,Model % NumberOfBulkElements 859!------------------------------------------------------------------------------ 860 Element => Model % Elements(i) 861 862 dim = Element % Type % Dimension 863 n = Element % Type % NumberOfNodes 864 865 Nodes % x(1:n) = nx( Element % NodeIndexes ) 866 Nodes % y(1:n) = ny( Element % NodeIndexes ) 867 Nodes % z(1:n) = nz( Element % NodeIndexes ) 868 869 IntegStuff = GaussPoints( Element ) 870 U_Integ => IntegStuff % u 871 V_Integ => IntegStuff % v 872 W_Integ => IntegStuff % w 873 S_Integ => IntegStuff % s 874 N_Integ = IntegStuff % n 875!------------------------------------------------------------------------------ 876 LocalMatrix = 0.0D0 877 DO t=1,N_Integ 878 u = U_Integ(t) 879 v = V_Integ(t) 880 w = W_Integ(t) 881!------------------------------------------------------------------------------ 882! Basis function values & derivatives at the integration point 883!------------------------------------------------------------------------------ 884 stat = ElementInfo( Element,Nodes,u,v,w,SqrtElementMetric, & 885 Basis,dBasisdx ) 886 887 s = SqrtElementMetric * S_Integ(t) 888 DO p=1,N 889 DO q=1,N 890 A = dBasisdx(p,Moved) * dBasisdx(q,Moved) 891 LocalMatrix(p,q) = LocalMatrix(p,q) + s*A 892 END DO 893 END DO 894 END DO 895 896 CALL CRS_GlueLocalMatrix( CMatrix,n,1,Element % NodeIndexes,LocalMatrix ) 897 END DO 898 899!------------------------------------------------------------------------------- 900 DO t = Model % NumberOfBulkElements + 1, Model % NumberOfBulkElements + & 901 Model % NumberOfBoundaryElements 902!------------------------------------------------------------------------------- 903 Element => Model % Elements(t) 904 n = Element % Type % NumberOfNodes 905 906 DO i=1,Model % NumberOfBCs 907 IF ( Element % BoundaryInfo % Constraint == Model % BCs(i) % Tag ) THEN 908 IF ( .NOT.ListGetLogical(Model % BCs(i) % Values,'Free Moving',stat) ) THEN 909 DO j=1,n 910 k = Element % NodeIndexes(j) 911 912 ForceVector(k) = Solution(k) 913 CALL CRS_ZeroRow( CMatrix,k ) 914 CALL CRS_SetMatrixElement( CMatrix,k,k,1.0D0 ) 915 END DO 916 END IF 917 END IF 918 END DO 919!------------------------------------------------------------------------------- 920 END DO 921!------------------------------------------------------------------------------- 922 923 CALL IterSolver( CMatrix,Solution,ForceVector,Solver ) 924 925!------------------------------------------------------------------------------- 926 END SUBROUTINE PoissonSolve 927!------------------------------------------------------------------------------- 928 929END MODULE FreeSurface 930 931!> \} 932