1!/*****************************************************************************/ 2! * 3! * Elmer/Ice, a glaciological add-on to Elmer 4! * http://elmerice.elmerfem.org 5! * 6! * 7! * This program is free software; you can redistribute it and/or 8! * modify it under the terms of the GNU General Public License 9! * as published by the Free Software Foundation; either version 2 10! * of the License, or (at your option) any later version. 11! * 12! * This program 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 15! * GNU General Public License for more details. 16! * 17! * You should have received a copy of the GNU General Public License 18! * along with this program (in file fem/GPL-2); if not, write to the 19! * Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, 20! * Boston, MA 02110-1301, USA. 21! * 22! *****************************************************************************/ 23! ****************************************************************************** 24! * 25! * Authors: Joe Todd 26! * Email: 27! * Web: http://elmerice.elmerfem.org 28! * 29! * 30! ***************************************************************************** 31 32!This moduled, loosely named 'CalvingGeometry' is for basically any 33!reusable routines for the 3D calving model. 34 35MODULE CalvingGeometry 36 37 USE Types 38 USE SParIterComm 39 USE MainUtils 40 41 IMPLICIT NONE 42 43 INTERFACE DoubleIntVectorSize 44 MODULE PROCEDURE DoubleIntVectorSizeP, DoubleIntVectorSizeA 45 END INTERFACE 46 47!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 48 ! Derived type for 3D crevasse group info 49 ! 50 ! Using the => Next, => Prev format like 51 ! variables_t, because there's no way of 52 ! knowing, a priori, how many we need. 53 ! 54 ! Actually the only use of this is borrowed by BasalMelt3D.F90, so its misnamed... 55!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 56 TYPE CrevasseGroup3D_t 57 INTEGER :: NumberOfNodes = 0, ID = 0 58 INTEGER, POINTER :: NodeNumbers(:) => NULL() 59 INTEGER, POINTER :: BoundaryNodes(:) => NULL(), FrontNodes(:) => NULL() !allocatable too? 60 REAL(KIND=dp) :: BoundingBox(4) !min_x, max_x, min_y, max_y 61 62 LOGICAL :: FrontConnected !Does the group touch the terminus? 63 TYPE(CrevasseGroup3D_t), POINTER :: Next => NULL(), Prev => NULL() 64 END TYPE CrevasseGroup3D_t 65 66!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 67 ! Derived type for a calving path defined by 68 ! the IsoSurface/Line solver. 69 ! (doubly linked list) 70!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 71 TYPE CrevassePath_t 72 INTEGER :: NumberOfNodes = 0, NumberOfElements = 0, ID = 0 73 INTEGER, POINTER :: NodeNumbers(:) => NULL(), ElementNumbers(:)=>NULL() 74! INTEGER :: Ends(2) 75 REAL(KIND=dp) :: Left, Right, Extent 76 TYPE(CrevassePath_t), POINTER :: Next => NULL(), Prev => NULL() 77 LOGICAL :: Valid = .TRUE. 78 END TYPE CrevassePath_t 79 80CONTAINS 81 82 83 !Returns the neighbours of a specified node using the matrix 84 !provided. 85 !Note the current definition of neighbours: 86 !Two nodes are neighbours if they are in the same bulk element 87 !NOT ONLY if they are joined by a bar... 88 !User may provide an inverse perm (InvPerm_in), or else this will recomputed 89 !each time (which would be pretty inefficient) 90 FUNCTION FindNodeNeighbours(NodeNumber, Matrix, Perm, DOFs, InvPerm_in) RESULT (Neighbours) 91 INTEGER :: NodeNumber, NoNeighbours, i, count, DOFs !<---!!! 92 TYPE(Matrix_t), POINTER :: Matrix 93 INTEGER, POINTER :: Perm(:), Neighbours(:), InvPerm(:) 94 INTEGER, POINTER, OPTIONAL, INTENT(IN) :: InvPerm_in(:) 95 LOGICAL :: Debug 96 Debug = .FALSE. 97 98 IF(PRESENT(InvPerm_in)) THEN 99 InvPerm => InvPerm_in 100 ELSE 101 IF(Debug) PRINT *, 'Debug FindNodeNeighbours, creating InvPerm' 102 InvPerm => CreateInvPerm(Perm) 103 END IF 104 105 NoNeighbours = Matrix % Rows((Perm(NodeNumber)*DOFs)+1) & 106 - Matrix % Rows(Perm(NodeNumber)*DOFs) 107 108 IF(MOD(NoNeighbours, DOFs).NE. 0) & 109 CALL FATAL("Geometry","This shouldn't have happened...") 110 111 !Each neighbour appears once per DOF, and there's also the current node thus: (x/DOFS) - 1... 112 NoNeighbours = (NoNeighbours / DOFs) - 1 113 114 ALLOCATE(Neighbours(NoNeighbours)) 115 Neighbours = 0 116 117 count = 0 118 119 DO i=Matrix % Rows(Perm(NodeNumber)*DOFs),& 120 (Matrix % Rows((Perm(NodeNumber)*DOFs)+1)-1) !move along the row 121 IF(MOD(i,DOFs) /= 0) CYCLE !Stored DOF1, DOF2, DOF3, only need every DOFth 122 IF(MOD(Matrix % Cols(i), DOFs) /= 0) CALL Fatal("Geometry:FindNodeNeighbours", & 123 "This is a programming error, Matrix structure is not what was expected.") 124 125 IF(InvPerm(Matrix % Cols(i)/DOFs) == NodeNumber) CYCLE !Not our own neighbour 126 count = count + 1 127 Neighbours(count) = & 128 InvPerm(Matrix % Cols(i)/DOFs) 129 END DO 130 131 IF(.NOT. PRESENT(InvPerm_in)) DEALLOCATE(InvPerm) 132 133 END FUNCTION FindNodeNeighbours 134 135 136 !----------------------------------------------------------------------------- 137 !Returns the 2D (x,y) Cartesian distance between two nodes 138 !NOTE: This isn't well programmed, should probably pass nodes... 139 FUNCTION NodeDist2D(Nodes, NodeNum1, NodeNum2 ) RESULT (dist) 140 TYPE(Nodes_t) :: Nodes 141 INTEGER :: NodeNum1, NodeNum2 142 REAL(KIND=dp) :: xdist,ydist,dist 143 !Pythagoras in 2D 144 xdist = Nodes % x(NodeNum1)& 145 - Nodes % x(NodeNum2) 146 ydist = Nodes % y(NodeNum1)& 147 - Nodes % y(NodeNum2) 148 !TODO: Can this be simplified? See Interpolation.f90 149 dist = ((xdist**2) + (ydist**2))**0.5 150 END FUNCTION NodeDist2D 151 152 !----------------------------------------------------------------------------- 153 !Returns the 3D Cartesian distance between two nodes 154 !NOTE: This isn't well programmed, should probably pass nodes... 155 FUNCTION NodeDist3D( Nodes, Node1, Node2 ) RESULT (dist) 156 TYPE(Nodes_t) :: Nodes 157 INTEGER :: Node1, Node2 158 REAL(KIND=dp) :: xdist,ydist,zdist,xydist,dist 159 !Pythagoras in 3D 160 xdist = Nodes % x(Node1)& 161 - Nodes % x(Node2) 162 ydist = Nodes % y(Node1)& 163 - Nodes % y(Node2) 164 zdist = Nodes % z(Node1)& 165 - Nodes % z(Node2) 166 !TODO: Can this be simplified? See Interpolation.f90 167 xydist = ((xdist**2) + (ydist**2))**0.5 168 dist = ((xydist**2) + (zdist**2))**0.5 169 END FUNCTION NodeDist3D 170 171 !----------------------------------------------------------------------------- 172 !Returns the inverse permutation table for a given perm and DOFs 173 !NOTE, differs from the definition of InvPerm currently used in 174 !Calving.F90 175 FUNCTION CreateInvPerm(Perm) RESULT(InvPerm) 176 INTEGER, POINTER :: Perm(:), InvPerm(:) 177 INTEGER :: i, j 178 179 ALLOCATE(InvPerm(MAXVAL(Perm))) 180 181 j = 0 182 DO i=1,SIZE(Perm) 183 IF(Perm(i) == 0) CYCLE 184 j = j + 1 185 InvPerm( Perm(i) ) = j 186 END DO 187 188 END FUNCTION CreateInvPerm 189 190 !----------------------------------------------------------------------------- 191 !Returns dx/dy for two given nodes 192 FUNCTION NodesGradXY(Nodes, Node1, Node2)RESULT(dxdy) 193 INTEGER :: Node1, Node2 194 TYPE(Nodes_t) :: Nodes 195 REAL(KIND=dp) :: dx,dy,dxdy 196 197 dx = Nodes % x(Node1) - Nodes % x(Node2) 198 dy = Nodes % y(Node1) - Nodes % y(Node2) 199 dxdy = dx/dy 200 END FUNCTION NodesGradXY 201 202 !----------------------------------------------------------------------------- 203 !Returns the number of decimal places of a real number 204 !which has been read from a text file (.e.g mesh.nodes) 205 !this differs from intrinsic PRECISION() because these 206 !numbers often have trailing 000s or 999s 207 FUNCTION RealAeps(in)RESULT(myaeps) 208 REAL(KIND=dp) :: in, toler, x, myaeps 209 INTEGER :: sigs, mag, decs 210 211 !Find how many decimal places 212 mag = FLOOR(LOG10(ABS(in))) + 1 !Order of magnitude of number 213 decs = PRECISION(in) - mag !total digits - magnitude = decimal places 214 215 toler = 10.0_dp**(-decs) 216 sigs = 0 217 x = in 218 219 DO WHILE (.TRUE.) 220 IF(ABS(x - NINT(x)) < toler) THEN !found the precision limit 221 EXIT 222 ELSE 223 sigs = sigs + 1 224 x = x * 10 !move the decimal point along 225 x = x - FLOOR(x) !convert number to O(1) so FLOOR doesn't reach integer limit 226 toler = toler * 10.0_dp !1 fewer remaining decimal places 227 END IF 228 END DO 229 myaeps = 10.0**(-sigs) 230 END FUNCTION RealAeps 231 232 !----------------------------------------------------------------------------- 233 ! Constructs paths of connected isoline (202) elements which intersect the 234 ! front. Each path will begin and end with a node where OnFront=.TRUE. 235 !----------------------------------------------------------------------------- 236 SUBROUTINE FindCrevassePaths(IsoMesh, OnFront, CrevassePaths, PathCount) 237 IMPLICIT NONE 238 TYPE(Mesh_t), POINTER :: IsoMesh 239 LOGICAL, ALLOCATABLE :: OnFront(:) 240 TYPE(CrevassePath_t), POINTER :: CrevassePaths 241 INTEGER :: PathCount 242 !---------------------------------------------- 243 TYPE(CrevassePath_t), POINTER :: CurrentPath 244 LOGICAL :: Found, Debug 245 INTEGER :: i,j,NodeCount,ElemCount, NextElem 246 INTEGER, ALLOCATABLE :: WorkElems(:), WorkNodes(:) 247 248 Debug = .FALSE. 249 PathCount = 1 250 251 !TODO assert all 202 elements 252 253 ALLOCATE(CrevassePaths) 254 CurrentPath => CrevassePaths 255 256 ALLOCATE(WorkElems(100), WorkNodes(100)) 257 WorkElems = 0; WorkNodes = 0 258 259 DO i=1, IsoMesh % NumberOfBulkElements 260 261 IF(ANY(OnFront(Isomesh % Elements(i) % NodeIndexes))) THEN 262 !Found an element with one node on calving front 263 264 IF(ElementPathID(CrevassePaths, i) /= 0) CYCLE !already in a path 265 266 !Starting a new group... 267 CurrentPath % ID = PathCount 268 IF(Debug) PRINT *, 'Potential calving isomesh element: ',i 269 270 ElemCount = 1 271 NextElem = i 272 273 !Identify which of the two nodes are on the front... 274 DO j=1,2 275 IF(OnFront(IsoMesh % Elements(i) % NodeIndexes(j))) EXIT 276 END DO 277 IF(j==3) CALL Fatal("FindCrevassePaths", "Couldn't find node on boundary") 278 279 !... and put it first in the list 280 WorkNodes(1) = IsoMesh % Elements(i) % NodeIndexes(j) 281 NodeCount = 2 282 283 !Follow the chain 284 DO WHILE(.TRUE.) 285 286 WorkElems(ElemCount) = NextElem 287 ElemCount = ElemCount + 1 288 !Put the other node into the list 289 DO j=1,2 290 IF(ANY(WorkNodes == IsoMesh % Elements(NextElem) % NodeIndexes(j))) CYCLE 291 WorkNodes(NodeCount) = IsoMesh % Elements(NextElem) % NodeIndexes(j) 292 NodeCount = NodeCount + 1 293 EXIT 294 END DO 295 296 !Look for element which contains previous element's node 297 Found = .FALSE. 298 DO j=1,IsoMesh % NumberOfBulkElements 299 IF(ANY(IsoMesh % Elements(j) % NodeIndexes == WorkNodes(NodeCount-1))) THEN 300 301 !already in another group (is this possible?) 302 IF(ElementPathID(CrevassePaths, j ) /= 0) CYCLE 303 !Already in current group 304 IF(ANY(WorkElems == j)) CYCLE 305 306 NextElem = j 307 Found = .TRUE. 308 EXIT 309 END IF 310 END DO 311 312 IF(.NOT. Found) EXIT 313 314 IF(ElemCount > SIZE(WorkElems)) THEN 315 IF(Debug) PRINT *, 'FindCrevassePaths, doubling size of element array.' 316 CALL DoubleIntVectorSize(WorkElems) 317 END IF 318 IF(NodeCount > SIZE(WorkNodes)) THEN 319 IF(Debug) PRINT *, 'FindCrevassePaths, doubling size of node array.' 320 CALL DoubleIntVectorSize(WorkNodes) 321 END IF 322 END DO 323 324 ElemCount = ElemCount - 1 325 NodeCount = NodeCount - 1 326 327 CurrentPath % NumberOfNodes = NodeCount 328 CurrentPath % NumberOfElements = ElemCount 329 330 ALLOCATE(CurrentPath % ElementNumbers(ElemCount), & 331 CurrentPath % NodeNumbers(NodeCount)) 332 333 CurrentPath % NodeNumbers = WorkNodes(1:NodeCount) 334 CurrentPath % ElementNumbers = WorkElems(1:ElemCount) 335 336 WorkNodes = 0 337 WorkElems = 0 338 339 ALLOCATE(CurrentPath % Next) 340 CurrentPath % Next % Prev => CurrentPath 341 CurrentPath => CurrentPath % Next 342 PathCount = PathCount + 1 343 END IF 344 END DO 345 346 !We always overshoot by one 347 PathCount = PathCount - 1 348 349 IF(PathCount > 0) THEN 350 PRINT *,'Number of crevasse paths: ', PathCount 351 CurrentPath % Prev % Next => NULL() 352 DEALLOCATE(CurrentPath) 353 ELSE 354 PRINT *,'No crevasse paths' 355 DEALLOCATE(CrevassePaths) 356 END IF 357 358 DEALLOCATE(WorkNodes, WorkElems) 359 360 END SUBROUTINE FindCrevassePaths 361 362 !Removes a CrevassePath from a linked list of CrevassePaths 363 SUBROUTINE RemoveCrevassePath(Path) 364 IMPLICIT NONE 365 TYPE(CrevassePath_t), POINTER :: Path 366 !------------------------------------------------ 367 IF(ASSOCIATED(Path % Prev)) Path % Prev % Next => Path % Next 368 IF(ASSOCIATED(Path % Next)) Path % Next % Prev => Path % Prev 369 370 IF(ASSOCIATED(Path % NodeNumbers)) DEALLOCATE(Path % NodeNumbers) 371 IF(ASSOCIATED(Path % ElementNumbers)) DEALLOCATE(Path % ElementNumbers) 372 DEALLOCATE(Path) 373 374 END SUBROUTINE RemoveCrevassePath 375 376 !-------------------------------------------------------------------- 377 ! 'Tidies up' isomesh and the CrevassePaths found by FindCrevassePaths 378 !-------------------------------------------------------------------- 379 ! This involves removing duplicate nodes, taking care to replace node 380 ! indexes in affected elements. This then allows easy removal of 381 ! 202 elements with zero length. 382 ! 383 ! Closed loops are removed from crevasse paths 384 !-------------------------------------------------------------------- 385 386 SUBROUTINE CheckCrevasseNodes(Mesh, CrevassePaths) 387 IMPLICIT NONE 388 TYPE(Mesh_t), POINTER :: Mesh 389 TYPE(CrevassePath_t), POINTER :: CrevassePaths 390 !------------------------------------------------- 391 TYPE(CrevassePath_t), POINTER :: CurrentPath,WorkPath 392 INTEGER :: i,j,ElNo,counter, ElementNumbers(2) 393 INTEGER, ALLOCATABLE :: ReplaceWithNode(:),WorkInt(:) 394 LOGICAL :: Debug 395 LOGICAL, ALLOCATABLE :: RemoveElement(:), RemoveNode(:), PathRemoveElement(:) 396 397 Debug = .FALSE. 398 399 ALLOCATE(RemoveNode(Mesh % NumberOfNodes),& 400 ReplaceWithNode(Mesh % NumberOfNodes),& 401 RemoveElement(Mesh % NumberOfBulkElements)) 402 RemoveNode = .FALSE. 403 RemoveElement = .FALSE. 404 ReplaceWithNode = 0 405 406 !Cycle mesh NODES, looking for duplicates and marking them for deletion 407 DO i=1,Mesh % NumberOfNodes 408 IF(RemoveNode(i)) CYCLE 409 DO j=1,Mesh % NumberOfNodes 410 IF(i==j .OR. RemoveNode(j)) CYCLE 411 IF(Mesh % Nodes % x(i) == Mesh % Nodes % x(j) .AND.& 412 Mesh % Nodes % y(i) == Mesh % Nodes % y(j) .AND.& 413 Mesh % Nodes % z(i) == Mesh % Nodes % z(j)) THEN 414 RemoveNode(j) = .TRUE. 415 ReplaceWithNode(j) = i 416 END IF 417 END DO 418 END DO 419 420 !Replace element nodeindexes where nodes are removed 421 DO i=1,Mesh % NumberOfBulkElements 422 DO j=1,SIZE(Mesh % Elements(i) % NodeIndexes) 423 IF(RemoveNode(Mesh % Elements(i) % NodeIndexes(j))) & 424 Mesh % Elements(i) % NodeIndexes(j) = & 425 ReplaceWithNode(Mesh % Elements(i) % NodeIndexes(j)) 426 END DO 427 END DO 428 429 !Mark elements with zero length (duplicate node indexes) for removal 430 DO i=1,Mesh % NumberOfBulkElements 431 IF(Mesh % Elements(i) % NodeIndexes(1) == Mesh % Elements(i) % NodeIndexes(2)) THEN 432 RemoveElement(i) = .TRUE. 433 IF(Debug) PRINT *,'debug, removing element: ',i,' with identical nodes: ',& 434 Mesh % Elements(i) % NodeIndexes(1) 435 END IF 436 END DO 437 438 IF(Debug) PRINT *,'Debug, removing ',COUNT(RemoveElement),' of ',SIZE(RemoveElement),' elements' 439 440 !Cycle paths, looking for nodes which are identical and removing them, joining up elements etc 441 CurrentPath => CrevassePaths 442 DO WHILE(ASSOCIATED(CurrentPath)) 443 444 IF(Debug) PRINT *,'Debug, Path: ',CurrentPath % ID,'initial no elems: ',& 445 CurrentPath % NumberOfElements,& 446 ' no nodes: ', CurrentPath % NumberOfNodes 447 448 ALLOCATE(WorkInt(CurrentPath % NumberOfElements)) 449 WorkInt = 0 450 counter = 0 451 452 !Mark pairs of duplicate elements in path for removal 453 ALLOCATE(PathRemoveElement(CurrentPath % NumberOfElements)) 454 PathRemoveElement = .FALSE. 455 456 IF(CurrentPath % NumberOfElements == 1) THEN 457 !Only has one element, remove 458 PathRemoveElement = .TRUE. 459 ELSE 460 DO i=1,CurrentPath % NumberOfElements-1 461 462 IF(PathRemoveElement(i)) CYCLE 463 ElementNumbers(1) = CurrentPath % ElementNumbers(i) 464 IF(RemoveElement(ElementNumbers(1))) CYCLE 465 466 j = i+1 467 IF(PathRemoveElement(j)) CYCLE 468 ElementNumbers(2) = CurrentPath % ElementNumbers(j) 469 IF(RemoveElement(ElementNumbers(2))) CYCLE 470 471 IF( ANY(Mesh % Elements(ElementNumbers(1)) % NodeIndexes == & 472 Mesh % Elements(ElementNumbers(2)) % NodeIndexes(1)) .AND. & 473 ANY(Mesh % Elements(ElementNumbers(1)) % NodeIndexes == & 474 Mesh % Elements(ElementNumbers(2)) % NodeIndexes(2)) ) THEN 475 PathRemoveElement(j) = .TRUE. 476 PathRemoveElement(i) = .TRUE. 477 IF(Debug) PRINT *,'Path: ',CurrentPath % ID,' removing identical elements: ',i,' ',j 478 END IF 479 480 END DO 481 482 483 !Check if entire crevasse group is a closed loop 484 ElementNumbers(1) = CurrentPath % ElementNumbers(1) 485 ElementNumbers(2) = CurrentPath % ElementNumbers(CurrentPath % NumberOfElements) 486 DO i=1,2 487 IF(.NOT. ANY(Mesh % Elements(CurrentPath % ElementNumbers(2)) % NodeIndexes == & 488 Mesh % Elements(ElementNumbers(1)) % NodeIndexes(i))) EXIT 489 END DO 490 IF(i==3) CALL Fatal("CheckCrevassePaths","Programming error: unable to determine first node") 491 IF(ANY(Mesh % Elements(ElementNumbers(2)) % NodeIndexes == & 492 Mesh % Elements(ElementNumbers(1)) % NodeIndexes(i))) THEN 493 PathRemoveElement = .TRUE. 494 IF(Debug) PRINT *,'Debug, removing path ',CurrentPath % ID,' because its entirely closed.' 495 END IF 496 497 !For each element 'i' in turn, cycle backwards through element list looking 498 !for element(i)'s nodes. If found, this indicates a closed loop which should 499 !be removed. 500 DO i=1,CurrentPath % NumberOfElements 501 IF(PathRemoveElement(i)) CYCLE 502 IF(RemoveElement(CurrentPath % ElementNumbers(i))) CYCLE 503 ElementNumbers(1) = CurrentPath % ElementNumbers(i) 504 505 DO j=CurrentPath % NumberOfElements,i+1,-1 !cycle backwards from end to i+1 506 IF(PathRemoveElement(j)) CYCLE 507 IF(RemoveElement(CurrentPath % ElementNumbers(j))) CYCLE 508 ElementNumbers(2) = CurrentPath % ElementNumbers(j) 509 510 IF( ANY(Mesh % Elements(ElementNumbers(1)) % NodeIndexes == & 511 Mesh % Elements(ElementNumbers(2)) % NodeIndexes(1)) .OR. & 512 ANY(Mesh % Elements(ElementNumbers(1)) % NodeIndexes == & 513 Mesh % Elements(ElementNumbers(2)) % NodeIndexes(2)) ) THEN 514 PathRemoveElement(i+1:j-1) = .TRUE. 515 IF(Debug) PRINT *,'CheckCrevasseNodes, & 516 &Removing a closed loop from ',i+1,' to ',j-1 517 END IF 518 519 END DO 520 521 END DO 522 END IF 523 524 !Replace CrevassePath % ElementNumbers based on previous removals 525 DO i=1,CurrentPath % NumberOfElements 526 IF(.NOT.RemoveElement(CurrentPath % ElementNumbers(i)) .AND. & 527 .NOT.PathRemoveElement(i)) THEN 528 counter = counter + 1 529 WorkInt(counter) = CurrentPath % ElementNumbers(i) 530 IF(Debug) THEN 531 PRINT *,'Debug, keeping element: ',i,' from path: ',CurrentPath % ID 532 PRINT *,'Debug, element global: ',CurrentPath % ElementNumbers(i),' and nodes :',& 533 Mesh % Elements(CurrentPath % ElementNumbers(i)) % NodeIndexes 534 END IF 535 ELSE 536 IF(Debug) THEN 537 PRINT *,'Debug, removing element: ',i,' from path: ',CurrentPath % ID 538 PRINT *,'Debug, element global: ',CurrentPath % ElementNumbers(i),' and nodes :',& 539 Mesh % Elements(CurrentPath % ElementNumbers(i)) % NodeIndexes 540 END IF 541 END IF 542 END DO 543 IF(counter < CurrentPath % NumberOfElements) THEN 544 IF(Debug) PRINT *,'debug, path loses ',CurrentPath % NumberOfElements - counter,& 545 ' of ',CurrentPath % NumberOfElements,' elements.' 546 547 CurrentPath % NumberOfElements = counter 548 DEALLOCATE(CurrentPath % ElementNumbers) 549 ALLOCATE(CurrentPath % ElementNumbers(counter)) 550 551 CurrentPath % ElementNumbers = WorkInt(1:counter) 552 END IF 553 DEALLOCATE(WorkInt,PathRemoveElement) 554 555 IF (CurrentPath % NumberOfElements <= 0) THEN 556 WorkPath => CurrentPath % Next 557 558 IF(ASSOCIATED(CurrentPath,CrevassePaths)) CrevassePaths => WorkPath 559 CALL RemoveCrevassePath(CurrentPath) 560 IF(Debug) CALL Info("CheckCrevasseNodes",& 561 "Removing a crevasse path with no elements") 562 CurrentPath => WorkPath 563 CYCLE 564 END IF 565 566 !Now reconstruct node list for path: 567 DEALLOCATE(CurrentPath % NodeNumbers) 568 CurrentPath % NumberOfNodes = CurrentPath % NumberOfElements + 1 569 ALLOCATE(CurrentPath % NodeNumbers(CurrentPath % NumberOfNodes)) 570 CurrentPath % NodeNumbers = 0 571 572 !First node 573 IF(CurrentPath % NumberOfElements >= 2) THEN 574 DO i=1,2 575 IF( ANY(Mesh % Elements(CurrentPath % ElementNumbers(2)) % NodeIndexes == & 576 Mesh % Elements(CurrentPath % ElementNumbers(1)) % NodeIndexes(i))) CYCLE 577 CurrentPath % NodeNumbers(1) = & 578 Mesh % Elements(CurrentPath % ElementNumbers(1)) % NodeIndexes(i) 579 580 IF(i==2) THEN !Reorder so that nodeindexes(1) and (2) are in chain order 581 Mesh % Elements(CurrentPath % ElementNumbers(1)) % NodeIndexes(2) = & 582 Mesh % Elements(CurrentPath % ElementNumbers(1)) % NodeIndexes(1) 583 Mesh % Elements(CurrentPath % ElementNumbers(1)) % NodeIndexes(1) = & 584 CurrentPath % NodeNumbers(1) 585 END IF 586 EXIT 587 END DO 588 ELSE !Rare, single element path, choice of first node is arbitrary... 589 CurrentPath % NodeNumbers(1) = & 590 Mesh % Elements(CurrentPath % ElementNumbers(1)) % NodeIndexes(1) 591 END IF 592 593 IF(Debug) PRINT *,'Path ',CurrentPath % ID,' has first node: ',CurrentPath % NodeNumbers(1) 594 595 !Follow the chain... 596 DO i=1,CurrentPath % NumberOfElements 597 ElNo = CurrentPath % ElementNumbers(i) 598 DO j=1,2 599 IF(ANY(CurrentPath % NodeNumbers == Mesh % Elements(ElNo) % NodeIndexes(j))) CYCLE 600 CurrentPath % NodeNumbers(i+1) = Mesh % Elements(ElNo) % NodeIndexes(j) 601 602 IF(j==1) THEN !Reorder so that nodeindexes(1) and (2) are in chain order 603 Mesh % Elements(CurrentPath % ElementNumbers(i)) % NodeIndexes(1) = & 604 Mesh % Elements(CurrentPath % ElementNumbers(i)) % NodeIndexes(2) 605 Mesh % Elements(CurrentPath % ElementNumbers(i)) % NodeIndexes(2) = & 606 CurrentPath % NodeNumbers(i+1) 607 END IF 608 609 EXIT 610 END DO 611 END DO 612 613 IF(Debug) PRINT *,'Debug, path ',CurrentPath % ID,' has nodes: ',CurrentPath % NodeNumbers 614 IF(ANY(CurrentPath % NodeNumbers == 0)) CALL Fatal("CheckCrevasseNodes","Failed to fill node indexes") 615 CurrentPath => CurrentPath % Next 616 END DO 617 618 END SUBROUTINE CheckCrevasseNodes 619 620 !---------------------------------------------------- 621 ! Checks paths for projectability and overlap 622 ! In case of overlap, smaller enclosed path is deleted 623 ! In case of unprojectability, nodes are moved laterally 624 ! to restore projectability. 625 !---------------------------------------------------- 626 ! NOTE: if this breaks, it could be due to two paths 627 ! sharing a node. Thinking about it, I see no reason 628 ! this should be an issue, but we'll see... 629 SUBROUTINE ValidateCrevassePaths(Mesh, CrevassePaths, FrontOrientation, PathCount) 630 IMPLICIT NONE 631 TYPE(Mesh_t), POINTER :: Mesh 632 TYPE(CrevassePath_t), POINTER :: CrevassePaths 633 REAL(KIND=dp) :: FrontOrientation(3) 634 INTEGER :: PathCount, First, Last, LeftIdx, RightIdx 635 !--------------------------------------------------- 636 REAL(KIND=dp) :: RotationMatrix(3,3), UnRotationMatrix(3,3), FrontDist, MaxDist, & 637 ShiftTo, Dir1(2), Dir2(2),a1(2),a2(2),b1(2),b2(2),intersect(2) 638 REAL(KIND=dp), ALLOCATABLE :: ConstrictDirection(:,:) 639 TYPE(CrevassePath_t), POINTER :: CurrentPath, OtherPath, WorkPath, LeftPath, RightPath 640 INTEGER :: i,j,k,n,ElNo,ShiftToMe, NodeNums(2),A,B,FirstIndex, LastIndex,Start 641 INTEGER, ALLOCATABLE :: WorkInt(:) 642 LOGICAL :: Debug, Shifted, CCW, ToLeft, Snakey, OtherRight, ShiftRightPath,headland 643 LOGICAL, ALLOCATABLE :: PathMoveNode(:), DeleteElement(:), BreakElement(:), & 644 FarNode(:), DeleteNode(:), Constriction(:) 645 646 Debug = .FALSE. 647 Snakey = .TRUE. 648 649 RotationMatrix = ComputeRotationMatrix(FrontOrientation) 650 UnRotationMatrix = TRANSPOSE(RotationMatrix) 651 652 ! Temporarily rotate the mesh 653 CALL RotateMesh(Mesh, RotationMatrix) 654 655 ! Find path %left, %right, %extent (width) 656 CALL ComputePathExtent(CrevassePaths, Mesh % Nodes, .TRUE.) 657 658 IF(Snakey) THEN 659 !----------------------------------------------------- 660 ! Paths should not 'snake' inwards in a narrow slit... 661 !----------------------------------------------------- 662 663 !it's insufficient to require that no nodes be 664 !further away than the two edge nodes. 665 !Instead, must ensure that no nodes are further away than any 666 !surrounding nodes. 667 668 !First need to determine path orientation 669 !with respect to front.... 670 671 CurrentPath => CrevassePaths 672 DO WHILE(ASSOCIATED(CurrentPath)) 673 674 !First and last node on path 675 First = CurrentPath % NodeNumbers(1) 676 Last = CurrentPath % NodeNumbers(CurrentPath % NumberOfNodes) 677 678 !if ToLeft, the crevasse path goes from right to left, from the 679 !perspective of someone sitting in the fjord, looking at the front 680 ToLeft = Mesh % Nodes % y(Last) > Mesh % Nodes % y(First) 681 682 IF(Debug) THEN 683 FrontDist = NodeDist3D(Mesh % Nodes,First, Last) 684 PRINT *,'PATH: ', CurrentPath % ID, ' FrontDist: ',FrontDist 685 PRINT *,'PATH: ', CurrentPath % ID, & 686 ' nonodes: ',CurrentPath % NumberOfNodes,& 687 ' noelems: ',CurrentPath % NumberOfElements 688 END IF 689 690 !Cycle path nodes, finding those which are too far away 691 ALLOCATE(FarNode(CurrentPath % NumberOfNodes), & 692 Constriction(CurrentPath % NumberOfNodes),& 693 ConstrictDirection(CurrentPath % NumberOfNodes,2)) 694 FarNode = .FALSE. 695 Constriction = .FALSE. 696 ConstrictDirection = 0.0_dp 697 698 !Determine which nodes have the potential to be constriction (based on angle) 699 !and compute constriction direction (i.e. which way the 'pointy bit' points...') 700 DO i=2,CurrentPath % NumberOfNodes-1 701 First = CurrentPath % NodeNumbers(i-1) 702 Last = CurrentPath % NodeNumbers(i+1) 703 n = CurrentPath % NodeNumbers(i) 704 705 CCW = ((Mesh % Nodes % y(n) - Mesh % Nodes % y(First)) * & 706 (Mesh % Nodes % z(Last) - Mesh % Nodes % z(First))) > & 707 ((Mesh % Nodes % z(n) - Mesh % Nodes % z(First)) * & 708 (Mesh % Nodes % y(Last) - Mesh % Nodes % y(First))) 709 710 IF(CCW .NEQV. ToLeft) THEN 711 Constriction(i) = .TRUE. 712 !Calculate constriction direction 713 714 Dir1(1) = Mesh % Nodes % y(n) - Mesh % Nodes % y(First) 715 Dir1(2) = Mesh % Nodes % z(n) - Mesh % Nodes % z(First) 716 Dir1 = Dir1 / ((Dir1(1)**2.0 + Dir1(2)**2.0) ** 0.5) 717 718 Dir2(1) = Mesh % Nodes % y(n) - Mesh % Nodes % y(Last) 719 Dir2(2) = Mesh % Nodes % z(n) - Mesh % Nodes % z(Last) 720 Dir2 = Dir2 / ((Dir2(1)**2.0 + Dir2(2)**2.0) ** 0.5) 721 722 ConstrictDirection(i,1) = Dir1(1) + Dir2(1) 723 ConstrictDirection(i,2) = Dir1(2) + Dir2(2) 724 ConstrictDirection(i,:) = ConstrictDirection(i,:) / & 725 ((ConstrictDirection(i,1)**2.0 + ConstrictDirection(i,2)**2.0) ** 0.5) 726 727 IF(Debug) PRINT *, 'Debug, node ',i,' dir1,2: ',Dir1, Dir2 728 IF(Debug) PRINT *, 'Debug, node ',i,' constriction direction: ',ConstrictDirection(i,:) 729 IF(Debug) PRINT *, 'Debug, node ',i,' xyz: ',Mesh % Nodes % x(n),Mesh % Nodes % y(n),Mesh % Nodes % z(n) 730 END IF 731 END DO 732 733 !First and last can always be constriction 734 Constriction(1) = .TRUE. 735 Constriction(SIZE(Constriction)) = .TRUE. 736 737 !Compute constriction direction for first and last 738 !We don't have info about the third node, so take orthogonal to 2 739 Last = CurrentPath % NodeNumbers(2) 740 n = CurrentPath % NodeNumbers(1) 741 Dir1(1) = Mesh % Nodes % y(n) - Mesh % Nodes % y(Last) 742 Dir1(2) = Mesh % Nodes % z(n) - Mesh % Nodes % z(Last) 743 Dir1 = Dir1 / ((Dir1(1)**2.0 + Dir1(2)**2.0) ** 0.5) 744 745 !Depending on which end of the chain we are, 746 !we take either the right or left orthogonal vector 747 IF(ToLeft) THEN 748 ConstrictDirection(1,1) = Dir1(2) 749 ConstrictDirection(1,2) = -1.0 * Dir1(1) 750 ELSE 751 ConstrictDirection(1,1) = -1.0 * Dir1(2) 752 ConstrictDirection(1,2) = Dir1(1) 753 END IF 754 IF(Debug) PRINT *, 'Debug, node 1 constriction direction: ',ConstrictDirection(1,:) 755 756 Last = CurrentPath % NodeNumbers(CurrentPath % NumberOfNodes - 1) 757 n = CurrentPath % NodeNumbers(CurrentPath % NumberOfNodes) 758 759 Dir1(1) = Mesh % Nodes % y(n) - Mesh % Nodes % y(Last) 760 Dir1(2) = Mesh % Nodes % z(n) - Mesh % Nodes % z(Last) 761 Dir1 = Dir1 / ((Dir1(1)**2.0 + Dir1(2)**2.0) ** 0.5) 762 763 IF(.NOT. ToLeft) THEN 764 ConstrictDirection(CurrentPath % NumberOfNodes,1) = Dir1(2) 765 ConstrictDirection(CurrentPath % NumberOfNodes,2) = -1.0 * Dir1(1) 766 ELSE 767 ConstrictDirection(CurrentPath % NumberOfNodes,1) = -1.0 * Dir1(2) 768 ConstrictDirection(CurrentPath % NumberOfNodes,2) = Dir1(1) 769 END IF 770 IF(Debug) PRINT *, 'Debug, node last constriction direction: ',& 771 ConstrictDirection(CurrentPath % NumberOfNodes,:) 772 773 !--------------------------------------- 774 ! Now that we have constrictions marked and directions computed, cycle nodes 775 776 DO i=1,CurrentPath % NumberOfNodes 777 IF(.NOT. Constriction(i)) CYCLE 778 779 DO j=CurrentPath % NumberOfNodes,i+1,-1 780 IF(.NOT. Constriction(j)) CYCLE 781 782 783 First = CurrentPath % NodeNumbers(i) 784 Last = CurrentPath % NodeNumbers(j) 785 786 !Check that these constrictions 'face' each other via dot product 787 Dir1(1) = Mesh % Nodes % y(Last) - Mesh % Nodes % y(First) 788 Dir1(2) = Mesh % Nodes % z(Last) - Mesh % Nodes % z(First) 789 Dir2(1) = -Dir1(1) 790 Dir2(2) = -Dir1(2) 791 792 !If the two constrictions aren't roughly facing each other: 793 ! < > rather than > < 794 ! then skip this combo 795 IF(SUM(ConstrictDirection(i,:)*Dir1) < 0) THEN 796 IF(Debug) PRINT *,'Constrictions ',i,j,' do not face each other 1: ',& 797 SUM(ConstrictDirection(i,:)*Dir1) 798 CYCLE 799 END IF 800 801 IF(SUM(ConstrictDirection(j,:)*Dir2) < 0) THEN 802 IF(Debug) PRINT *,'Constrictions ',j,i,' do not face each other 2: ',& 803 SUM(ConstrictDirection(j,:)*Dir2) 804 CYCLE 805 END IF 806 807 IF(Debug) PRINT *,'Constrictions ',i,j,' do face each other: ',& 808 SUM(ConstrictDirection(i,:)*Dir1) 809 810 !test that the line drawn between the constriction doesn't intersect 811 !any intermediate elements as this indicates 812 !crossing a headland (difficult to draw - but it's bad news) 813 ! 814 ! - --- ---- - 815 ! \/ \ / \/ 816 ! ---- 817 ! 818 819 a1(1) = Mesh % Nodes % y(First) 820 a1(2) = Mesh % Nodes % z(First) 821 a2(1) = Mesh % Nodes % y(Last) 822 a2(2) = Mesh % Nodes % z(Last) 823 headland = .FALSE. 824 DO k=i+1,j-2 825 b1(1) = Mesh % Nodes % y(k) 826 b1(2) = Mesh % Nodes % z(k) 827 b2(1) = Mesh % Nodes % y(k+1) 828 b2(2) = Mesh % Nodes % z(k+1) 829 830 CALL LineSegmentsIntersect(a1,a2,b1,b2,intersect,headland) 831 IF(headland) EXIT 832 END DO 833 IF(headland) CYCLE 834 835 MaxDist = NodeDist3D(Mesh % Nodes,First, Last) 836 837 DO k=i+1,j-1 838 IF(FarNode(k)) CYCLE 839 840 n = CurrentPath % NodeNumbers(k) 841 842 IF((NodeDist3D(Mesh % Nodes, First, n) <= MaxDist) .AND. & 843 (NodeDist3D(Mesh % Nodes, Last, n) <= MaxDist)) CYCLE !within range 844 845 FarNode(k) = .TRUE. 846 IF(Debug) PRINT *,'Far node: ',k,' xyz: ',Mesh % Nodes % x(n),& 847 Mesh % Nodes % y(n),Mesh % Nodes % z(n) 848 849 END DO 850 END DO 851 END DO 852 853 !Cycle elements, marking those which need to be adjusted 854 ALLOCATE(BreakElement(CurrentPath % NumberOfElements),& 855 DeleteElement(CurrentPath % NumberOfElements)) 856 BreakElement = .FALSE. 857 DeleteElement = .FALSE. 858 859 DO i=1,CurrentPath % NumberOfElements 860 IF(ANY(FarNode(i:i+1))) THEN 861 IF(ALL(FarNode(i:i+1))) THEN 862 DeleteElement(i) = .TRUE. 863 IF(Debug) PRINT *,'PATH: ', CurrentPath % ID, ' element: ',i,' is deleted.' 864 ELSE 865 BreakElement(i) = .TRUE. 866 IF(Debug) PRINT *,'PATH: ', CurrentPath % ID, ' element: ',i,' is broken.' 867 END IF 868 END IF 869 END DO 870 871 DO i=1,CurrentPath % NumberOfElements 872 IF((.NOT. BreakElement(i)) .OR. DeleteElement(i)) CYCLE 873 !This element needs to be adjusted 874 DO j=i+1,CurrentPath % NumberOfElements 875 IF(.NOT. (BreakElement(j) .OR. DeleteElement(j))) & 876 CALL Fatal("ValidateCrevasseGroups","Programming error in maintaining aspect ratio") 877 IF(DeleteElement(j)) CYCLE 878 !This is the next 'break element' after i 879 !Determine which nodes we keep 880 881 IF((CurrentPath % NodeNumbers(j) /= & 882 Mesh % Elements(CurrentPath % ElementNumbers(j)) % NodeIndexes(1)) .OR. & 883 (CurrentPath % NodeNumbers(j+1) /= & 884 Mesh % Elements(CurrentPath % ElementNumbers(j)) % NodeIndexes(2))) THEN 885 886 CALL Fatal("ValidateCrevassePaths", "Chain building error") 887 END IF 888 889 Mesh % Elements(CurrentPath % ElementNumbers(i)) % NodeIndexes(2) = & 890 Mesh % Elements(CurrentPath % ElementNumbers(j)) % NodeIndexes(2) 891 892 !We now want to delete it, because we only keep one from each broken pair 893 DeleteElement(j) = .TRUE. 894 EXIT !we paired this one, move on 895 END DO 896 END DO 897 898 !Delete the elements and nodes 899 IF(COUNT(DeleteElement) > 0) THEN 900 !elements 901 ALLOCATE(WorkInt(COUNT(.NOT. DeleteElement))) 902 WorkInt = PACK(CurrentPath % ElementNumbers,.NOT.DeleteElement) 903 904 DEALLOCATE(CurrentPath % ElementNumbers) 905 ALLOCATE(CurrentPath % ElementNumbers(SIZE(WorkInt))) 906 907 CurrentPath % ElementNumbers = WorkInt 908 CurrentPath % NumberOfElements = SIZE(WorkInt) 909 DEALLOCATE(WorkInt) 910 911 !nodes 912 ALLOCATE(WorkInt(COUNT(.NOT. FarNode))) 913 WorkInt = PACK(CurrentPath % NodeNumbers, .NOT.FarNode) 914 915 DEALLOCATE(CurrentPath % NodeNumbers) 916 ALLOCATE(CurrentPath % NodeNumbers(SIZE(WorkInt))) 917 918 CurrentPath % NodeNumbers = WorkInt 919 CurrentPath % NumberOfNodes = SIZE(WorkInt) 920 DEALLOCATE(WorkInt) 921 END IF 922 923 DEALLOCATE(FarNode, Constriction, ConstrictDirection, BreakElement, DeleteElement) 924 CurrentPath => CurrentPath % Next 925 END DO 926 END IF !Snakey 927 928 !Update Left, Right & Extent 929 CALL ComputePathExtent(CrevassePaths, Mesh % Nodes, .TRUE.) 930 931 !----------------------------------------------------- 932 ! Move nodes from crevassepaths which aren't projectable 933 !----------------------------------------------------- 934 ! 1) Path elements are ordered as a chain 935 ! 2) Path % Element(i) has nodes i, i+1 936 ! 937 ! Go through CrevassePath nodes, marking those 938 ! which are 'shadowed' by further away elements. 939 !----------------------------------------------------- 940 941 CurrentPath => CrevassePaths 942 DO WHILE(ASSOCIATED(CurrentPath)) 943 944 ALLOCATE(PathMoveNode(CurrentPath % NumberOfNodes)) 945 PathMoveNode = .FALSE. 946 947 DO i=1,CurrentPath % NumberOfNodes 948 n = CurrentPath % NodeNumbers(i) 949 DO j=1,CurrentPath % NumberOfElements 950 ElNo = CurrentPath % ElementNumbers(j) 951 NodeNums = Mesh % Elements(ElNo) % NodeIndexes 952 IF(ANY(NodeNums == n)) CYCLE !Node is in element, skip 953 !Check if node lies between element nodes 954 IF( (Mesh % Nodes % y(NodeNums(1)) > Mesh % Nodes % y(n)) .NEQV. & 955 (Mesh % Nodes % y(NodeNums(2)) > Mesh % Nodes % y(n)) ) THEN 956 !Check the node is in front of the element 957 958 A = MINLOC(Mesh % Nodes % z(NodeNums),1) 959 B = MAXLOC(Mesh % Nodes % z(NodeNums),1) 960 CCW = ((Mesh % Nodes % y(n) - Mesh % Nodes % y(NodeNums(A))) * & 961 (Mesh % Nodes % z(NodeNums(B)) - Mesh % Nodes % z(NodeNums(A)))) > & 962 ((Mesh % Nodes % z(n) - Mesh % Nodes % z(NodeNums(A))) * & 963 (Mesh % Nodes % y(NodeNums(B)) - Mesh % Nodes % y(NodeNums(A)))) 964 965 ToLeft = Mesh % Nodes % y(NodeNums(A)) > Mesh % Nodes % y(NodeNums(B)) 966 967 IF(CCW .EQV. ToLeft) THEN 968 !Node should be removed 969 PathMoveNode(i) = .TRUE. 970 EXIT 971 END IF 972 973 END IF 974 END DO 975 END DO 976 977 IF(Debug) THEN 978 PRINT *,'Path ',CurrentPath % ID,' has ',& 979 COUNT(PathMoveNode),' nodes which need to be shifted.' 980 981 DO i=1,CurrentPath % NumberOfNodes 982 IF(.NOT. PathMoveNode(i)) CYCLE 983 PRINT *,'Need to move node: ',i,' y: ',& 984 Mesh % Nodes % y(CurrentPath % NodeNumbers(i)),& 985 ' z: ',Mesh % Nodes % z(CurrentPath % NodeNumbers(i)) 986 987 END DO 988 END IF 989 990 !Now that nodes have been marked as shadowed, identify chains 991 !and the location of the node to which these groups of nodes should be moved. 992 Shifted = .TRUE. 993 Start = 1 994 DO WHILE(Shifted) 995 Shifted = .FALSE. 996 997 DO i=Start,CurrentPath % NumberOfNodes 998 IF(PathMoveNode(i)) THEN 999 IF(.NOT. Shifted) THEN 1000 Shifted = .TRUE. 1001 FirstIndex = i 1002 END IF 1003 LastIndex = i 1004 ELSE 1005 IF(Shifted) EXIT 1006 END IF 1007 END DO 1008 IF(.NOT. Shifted) EXIT 1009 1010 !We have identified a chain from FirstIndex to LastIndex which need to be moved. 1011 !They should be moved to either FirstIndex-1 or LastIndex+1 1012 !(Whichever is further back) 1013 !Note special case at start and end of path 1014 IF(FirstIndex == 1) THEN 1015 ShiftToMe = CurrentPath % NodeNumbers(LastIndex+1) 1016 ELSE IF(LastIndex == CurrentPath % NumberOfNodes) THEN 1017 ShiftToMe = CurrentPath % NodeNumbers(FirstIndex-1) 1018 ELSE IF(Mesh % Nodes % z(CurrentPath % NodeNumbers(FirstIndex-1)) <& 1019 Mesh % Nodes % z(CurrentPath % NodeNumbers(LastIndex+1))) THEN 1020 ShiftToMe = CurrentPath % NodeNumbers(FirstIndex-1) 1021 ELSE 1022 ShiftToMe = CurrentPath % NodeNumbers(LastIndex+1) 1023 END IF 1024 1025 Mesh % Nodes % y(CurrentPath % NodeNumbers(FirstIndex:LastIndex)) = & 1026 Mesh % Nodes % y(ShiftToMe) 1027 1028 IF(Debug) PRINT *,'Shifting nodes ',FirstIndex,' to ',LastIndex,& 1029 ' to point: ',Mesh % Nodes % y(ShiftToMe) 1030 Start = LastIndex + 1 1031 END DO 1032 1033 DEALLOCATE(PathMoveNode) 1034 CurrentPath => CurrentPath % Next 1035 END DO 1036 1037 !NOTE: probably not really necessary here, Shifted nodes don't extend 1038 !the extent 1039 !Update Left, Right & Extent 1040 CALL ComputePathExtent(CrevassePaths, Mesh % Nodes, .TRUE.) 1041 1042 !-------------------------------------------------------- 1043 ! Remove crevassepaths which are contained within others. 1044 !-------------------------------------------------------- 1045 ! 1) All crevasse paths start and end on the calving front. 1046 ! 2) Crevasse paths can't cross each other. 1047 ! 1048 ! Thus, iff a crevasse path is surrounded laterally by 1049 ! another single crevasse path, we remove it, because 1050 ! it must be contained by the larger one. 1051 !-------------------------------------------------------- 1052 1053 CurrentPath => CrevassePaths 1054 DO WHILE(ASSOCIATED(CurrentPath)) 1055 1056 OtherPath => CrevassePaths 1057 DO WHILE(ASSOCIATED(OtherPath)) 1058 IF(ASSOCIATED(OtherPath, CurrentPath)) THEN 1059 OtherPath => OtherPath % Next 1060 CYCLE 1061 END IF 1062 1063 IF((CurrentPath % Left >= OtherPath % Left) .AND. & 1064 (CurrentPath % Right <= OtherPath % Right)) THEN!contained within 1065 CurrentPath % Valid = .FALSE. 1066 IF(Debug) PRINT *,'Debug, marked path ',CurrentPath % ID,' for deletion & 1067 &because its contained within path ',OtherPath % ID 1068 END IF 1069 OtherPath => OtherPath % Next 1070 END DO 1071 1072 CurrentPath => CurrentPath % Next 1073 END DO 1074 1075 !Actually remove previous marked 1076 CurrentPath => CrevassePaths 1077 DO WHILE(ASSOCIATED(CurrentPath)) 1078 WorkPath => CurrentPath % Next 1079 1080 IF(.NOT. CurrentPath % Valid) THEN 1081 IF(ASSOCIATED(CurrentPath,CrevassePaths)) CrevassePaths => WorkPath 1082 CALL RemoveCrevassePath(CurrentPath) 1083 IF(Debug) CALL Info("ValidateCrevassePaths","Removing a crevasse path") 1084 END IF 1085 CurrentPath => WorkPath 1086 END DO 1087 1088 !------------------------------------------------- 1089 ! Check for paths partly obscuring each other 1090 ! (fully obscured are dealt with above) 1091 !------------------------------------------------- 1092 ! If paths partially overlap, the overlapping nodes 1093 ! of whichever path is seaward are moved. 1094 ! i.e. the larger calving event takes precedent 1095 !------------------------------------------------- 1096 1097 CurrentPath => CrevassePaths 1098 DO WHILE(ASSOCIATED(CurrentPath)) 1099 1100 OtherPath => CrevassePaths 1101 DO WHILE(ASSOCIATED(OtherPath)) 1102 IF(ASSOCIATED(OtherPath, CurrentPath)) THEN 1103 OtherPath => OtherPath % Next 1104 CYCLE 1105 END IF 1106 1107 IF((CurrentPath % Left < OtherPath % Right) .EQV. & 1108 (OtherPath % Left < CurrentPath % Right)) THEN !overlap 1109 1110 IF(Debug) PRINT *,'Debug, paths: ',CurrentPath % ID, OtherPath % ID,' partially overlap' 1111 1112 !Is the other path to the right or left? 1113 OtherRight = CurrentPath % Right < OtherPath % Right 1114 1115 !Check not fully contained - should have been dealt with above 1116 IF((CurrentPath % Right > OtherPath % Right) .NEQV. & 1117 (CurrentPath % Left > OtherPath % Left)) THEN 1118 CALL Warn("ValidateCrevassePaths","Encountered full overlap which & 1119 &should already have been taken care of! OK if this is rare, & 1120 &otherwise maybe programming error") 1121 END IF 1122 1123 IF(OtherRight) THEN 1124 RightPath => OtherPath 1125 LeftPath => CurrentPath 1126 ELSE 1127 RightPath => CurrentPath 1128 LeftPath => OtherPath 1129 END IF 1130 1131 !Find the left and rightmost nodes of the two paths 1132 DO i=1,LeftPath % NumberOfNodes 1133 IF(Debug) PRINT *,'Debug, node ',i,' of leftpath: ',& 1134 Mesh % Nodes % y(LeftPath % NodeNumbers(i)), LeftPath % Right 1135 1136 IF(Mesh % Nodes % y(LeftPath % NodeNumbers(i)) >= LeftPath % Right) LeftIdx = i 1137 END DO 1138 1139 DO i=1,RightPath % NumberOfNodes 1140 IF(Debug) PRINT *,'Debug, node ',i,' of rightpath: ',& 1141 Mesh % Nodes % y(RightPath % NodeNumbers(i)), RightPath % Left 1142 1143 IF(Mesh % Nodes % y(RightPath % NodeNumbers(i)) <= RightPath % Left) RightIdx = i 1144 END DO 1145 1146 !See which is further forward. 1147 ShiftRightPath = Mesh % Nodes % z(LeftPath % NodeNumbers(LeftIdx)) < & 1148 Mesh % Nodes % z(RightPath % NodeNumbers(RightIdx)) 1149 1150 IF(ShiftRightPath) THEN 1151 ShiftTo = Mesh % Nodes % y(LeftPath % NodeNumbers(LeftIdx)) 1152 DO i=1,RightPath % NumberOfNodes 1153 IF(Mesh % Nodes % y(RightPath % NodeNumbers(i)) < ShiftTo) THEN 1154 IF(Debug) PRINT *,'Debug, overlap shifting right node ',i,' path '& 1155 ,RightPath % ID,' from ', Mesh % Nodes % y(RightPath % NodeNumbers(i)),& 1156 ' to ',ShiftTo 1157 Mesh % Nodes % y(RightPath % NodeNumbers(i)) = ShiftTo 1158 END IF 1159 END DO 1160 CALL ComputePathExtent(RightPath, Mesh % Nodes, .FALSE.) 1161 1162 ELSE 1163 ShiftTo = Mesh % Nodes % y(RightPath % NodeNumbers(RightIdx)) 1164 DO i=1,LeftPath % NumberOfNodes 1165 IF(Mesh % Nodes % y(LeftPath % NodeNumbers(i)) > ShiftTo) THEN 1166 IF(Debug) PRINT *,'Debug, overlap shifting left node ',i,' path ',& 1167 LeftPath % ID,' from ',Mesh % Nodes % y(LeftPath % NodeNumbers(i)),& 1168 ' to ',ShiftTo 1169 Mesh % Nodes % y(LeftPath % NodeNumbers(i)) = ShiftTo 1170 END IF 1171 END DO 1172 CALL ComputePathExtent(LeftPath, Mesh % Nodes, .FALSE.) 1173 1174 END IF 1175 END IF 1176 1177 OtherPath => OtherPath % Next 1178 END DO 1179 1180 CurrentPath => CurrentPath % Next 1181 END DO 1182 1183 !----------------------------------------------------------------------- 1184 ! Remove elements whose nodes are in a vertical line 1185 ! (to prevent potential issues in interp) 1186 !----------------------------------------------------------------------- 1187 ! This occurs due to the shifting which occurs above. 1188 ! NOTE: This breaks the assumption that element(i) has nodes (i) & (i+1) 1189 ! It also breaks the chain! Currently OK but don't rely on this below this 1190 ! point, or in Calving3D.F90 1191 !----------------------------------------------------------------------- 1192 1193 CurrentPath => CrevassePaths 1194 DO WHILE(ASSOCIATED(CurrentPath)) 1195 1196 ALLOCATE(DeleteElement(CurrentPath % NumberOfElements),& 1197 DeleteNode(CurrentPath % NumberOfNodes)) 1198 DeleteElement = .FALSE. 1199 DeleteNode = .FALSE. 1200 1201 DO i=1,CurrentPath % NumberOfElements 1202 !Element i is composed of nodes i,i+1 1203 IF(Mesh % Nodes % y(CurrentPath % NodeNumbers(i)) == & 1204 Mesh % Nodes % y(CurrentPath % NodeNumbers(i+1))) THEN 1205 DeleteElement(i) = .TRUE. 1206 IF(Debug) PRINT *,'Debug, deleting element: ',i,' from path: ',& 1207 CurrentPath % ID,' because its a straight line' 1208 END IF 1209 END DO 1210 1211 IF(DeleteElement(1)) DeleteNode(1) = .TRUE. 1212 IF(DeleteElement(SIZE(DeleteElement))) DeleteNode(SIZE(DeleteNode)) = .TRUE. 1213 1214 DO i=2,CurrentPath % NumberOfNodes-1 1215 IF(DeleteElement(i-1) .AND. DeleteElement(i)) DeleteNode(i) = .TRUE. 1216 END DO 1217 1218 !Delete them 1219 IF(COUNT(DeleteElement) > 0) THEN 1220 !elements 1221 ALLOCATE(WorkInt(COUNT(.NOT. DeleteElement))) 1222 WorkInt = PACK(CurrentPath % ElementNumbers,.NOT.DeleteElement) 1223 1224 DEALLOCATE(CurrentPath % ElementNumbers) 1225 ALLOCATE(CurrentPath % ElementNumbers(SIZE(WorkInt))) 1226 1227 CurrentPath % ElementNumbers = WorkInt 1228 CurrentPath % NumberOfElements = SIZE(WorkInt) 1229 DEALLOCATE(WorkInt) 1230 1231 !nodes 1232 ALLOCATE(WorkInt(COUNT(.NOT. DeleteNode))) 1233 WorkInt = PACK(CurrentPath % NodeNumbers, .NOT.DeleteNode) 1234 1235 DEALLOCATE(CurrentPath % NodeNumbers) 1236 ALLOCATE(CurrentPath % NodeNumbers(SIZE(WorkInt))) 1237 1238 CurrentPath % NodeNumbers = WorkInt 1239 CurrentPath % NumberOfNodes = SIZE(WorkInt) 1240 DEALLOCATE(WorkInt) 1241 END IF 1242 1243 DEALLOCATE(DeleteElement, DeleteNode) 1244 CurrentPath => CurrentPath % Next 1245 END DO 1246 1247 !-------------------------------------------------------- 1248 ! Put the mesh back 1249 !-------------------------------------------------------- 1250 CALL RotateMesh(Mesh, UnRotationMatrix) 1251 1252 END SUBROUTINE ValidateCrevassePaths 1253 1254 !Calculates the left and rightmost extent, and the difference (width) of 1255 !Path, given the node locations in Nodes. 1256 SUBROUTINE ComputePathExtent(CrevassePaths, Nodes, DoAll) 1257 TYPE(CrevassePath_t), POINTER :: CrevassePaths 1258 TYPE(Nodes_t), POINTER :: Nodes 1259 LOGICAL :: DoAll 1260 !----------------------------------------------- 1261 TYPE(CrevassePath_t), POINTER :: CurrentPath 1262 INTEGER :: n 1263 1264 CurrentPath => CrevassePaths 1265 DO WHILE(ASSOCIATED(CurrentPath)) 1266 CurrentPath % Left = HUGE(1.0_dp) 1267 CurrentPath % Right = -1.0*HUGE(1.0_dp) 1268 1269 n = CurrentPath % NumberOfNodes 1270 1271 CurrentPath % Left = MINVAL(Nodes % y(CurrentPath % NodeNumbers)) 1272 CurrentPath % Right = MAXVAL(Nodes % y(CurrentPath % NodeNumbers)) 1273 1274 CurrentPath % Extent = CurrentPath % Right - CurrentPath % Left 1275 1276 CurrentPath => CurrentPath % Next 1277 1278 IF(.NOT. DoAll) EXIT 1279 END DO 1280 1281 END SUBROUTINE ComputePathExtent 1282 1283 !----------------------------------------------------------------------------- 1284 ! Returns the Path ID of the CrevassePath_t which contains the given element 1285 ! 0 if not found 1286 !----------------------------------------------------------------------------- 1287 FUNCTION ElementPathID(CrevassePaths, ElementNo) RESULT(ID) 1288 TYPE(CrevassePath_t), POINTER :: CrevassePaths 1289 INTEGER :: ElementNo, ID 1290 !---------------------------------------------- 1291 TYPE(CrevassePath_t), POINTER :: CurrentPath 1292 1293 ID = 0 1294 1295 CurrentPath => CrevassePaths 1296 DO WHILE(ASSOCIATED(CurrentPath)) 1297 IF(ASSOCIATED(CurrentPath % ElementNumbers)) THEN 1298 IF(ANY(CurrentPath % ElementNumbers == ElementNo)) THEN 1299 ID = CurrentPath % ID 1300 EXIT 1301 END IF 1302 END IF 1303 CurrentPath => CurrentPath % Next 1304 END DO 1305 1306 END FUNCTION ElementPathID 1307 1308 !----------------------------------------------------------------------------- 1309 ! Constructs groups of nodes which fall below a given threshold for some variable 1310 ! Used with the result of ProjectCalving, it groups nodes which have crevasse 1311 ! penetration beyond the threshold. 1312 !----------------------------------------------------------------------------- 1313 SUBROUTINE FindCrevasseGroups(Mesh, Variable, Neighbours, Threshold, Groups) 1314 IMPLICIT NONE 1315 1316 TYPE(Mesh_t), POINTER :: Mesh 1317 TYPE(Variable_t), POINTER :: Variable 1318 INTEGER, POINTER :: Neighbours(:,:) 1319 TYPE(CrevasseGroup3D_t), POINTER :: Groups, CurrentGroup 1320 REAL(KIND=dp) :: Threshold 1321 !--------------------------------------- 1322 INTEGER :: i, ID 1323 REAL(KIND=dp), POINTER :: Values(:) 1324 INTEGER, POINTER :: VPerm(:) 1325 INTEGER, ALLOCATABLE :: WorkInt(:) 1326 LOGICAL, ALLOCATABLE :: Condition(:) 1327 LOGICAL :: First, Debug 1328 1329 Debug = .FALSE. 1330 1331 Values => Variable % Values 1332 VPerm => Variable % Perm 1333 1334 ALLOCATE(Condition(Mesh % NumberOfNodes)) 1335 DO i=1, Mesh % NumberOfNodes 1336 1337 IF(VPerm(i) <= 0) THEN 1338 Condition(i) = .FALSE. 1339 ELSE IF(Values(VPerm(i)) < Threshold) THEN 1340 Condition(i) = .TRUE. 1341 ELSE 1342 Condition(i) = .FALSE. 1343 END IF 1344 1345 END DO 1346 1347 First = .TRUE. 1348 ID = 1 1349 DO i=1,Mesh % NumberOfNodes 1350 IF(.NOT. Condition(i)) CYCLE 1351 1352 IF(Debug) PRINT *,'PE:', ParEnv % MyPE,' debug, new group' 1353 1354 IF(First) THEN 1355 ALLOCATE(CurrentGroup) 1356 Groups => CurrentGroup 1357 First = .FALSE. 1358 ELSE 1359 ALLOCATE(CurrentGroup % Next) 1360 CurrentGroup % Next % Prev => CurrentGroup 1361 CurrentGroup => CurrentGroup % Next 1362 END IF 1363 1364 CurrentGroup % ID = ID 1365 ID = ID + 1 1366 1367 ALLOCATE(CurrentGroup % NodeNumbers(500)) 1368 CurrentGroup % NumberOfNodes = 1 1369 1370 !Add node to group and switch it off 1371 CurrentGroup % NodeNumbers(CurrentGroup % NumberOfNodes) = i 1372 Condition(i) = .FALSE. 1373 1374 !Search neighbours 1375 CALL SearchNeighbours(i, Neighbours, CurrentGroup, Condition) 1376 1377 ALLOCATE(WorkInt(CurrentGroup % NumberOfNodes)) 1378 WorkInt = CurrentGroup % NodeNumbers(1:CurrentGroup % NumberOfNodes) 1379 DEALLOCATE(CurrentGroup % NodeNumbers) 1380 ALLOCATE(CurrentGroup % NodeNumbers(CurrentGroup % NumberOfNodes)) 1381 CurrentGroup % NodeNumbers = WorkInt 1382 DEALLOCATE(WorkInt) 1383 1384 CALL UpdateCGrpBB(CurrentGroup, Mesh) 1385 END DO 1386 1387 IF(Debug) THEN 1388 CurrentGroup => Groups 1389 i=1 1390 DO WHILE(ASSOCIATED(CurrentGroup)) 1391 PRINT *,'group: ',i,' has ', CurrentGroup % NumberOfNodes,' nodes.' 1392 i = i + 1 1393 CurrentGroup => CurrentGroup % Next 1394 END DO 1395 END IF 1396 END SUBROUTINE FindCrevasseGroups 1397 1398 SUBROUTINE DeallocateCrevasseGroup(CGrp) 1399 TYPE(CrevasseGroup3D_t), POINTER :: CGrp 1400 1401 IF(ASSOCIATED(CGrp % Next)) CGrp % Next % Prev => CGrp % Prev 1402 IF(ASSOCIATED(CGrp % Prev)) CGrp % Prev % Next => CGrp % Next 1403 1404 IF(ASSOCIATED(CGrp % NodeNumbers)) DEALLOCATE(CGrp % NodeNumbers) 1405 IF(ASSOCIATED(CGrp % FrontNodes)) DEALLOCATE(CGrp % FrontNodes) 1406 IF(ASSOCIATED(CGrp % BoundaryNodes)) DEALLOCATE(CGrp % BoundaryNodes) 1407 1408 DEALLOCATE(CGrp) 1409 1410 END SUBROUTINE DeallocateCrevasseGroup 1411 1412 !Update the Bounding Box of a CrevasseGroup 1413 SUBROUTINE UpdateCGrpBB(CGrp, Mesh) 1414 TYPE(CrevasseGroup3D_t), POINTER :: CGrp 1415 TYPE(Mesh_t), POINTER :: Mesh 1416 1417 CGrp % BoundingBox(1) = MINVAL(Mesh % Nodes % x(CGrp % NodeNumbers)) 1418 CGrp % BoundingBox(2) = MAXVAL(Mesh % Nodes % x(CGrp % NodeNumbers)) 1419 CGrp % BoundingBox(3) = MINVAL(Mesh % Nodes % y(CGrp % NodeNumbers)) 1420 CGrp % BoundingBox(4) = MAXVAL(Mesh % Nodes % y(CGrp % NodeNumbers)) 1421 1422 END SUBROUTINE UpdateCGrpBB 1423 1424 !Add a list of points to a CrevasseGroup3D object 1425 !Don't need to pass the mesh because we're just adding 1426 !point indices 1427 SUBROUTINE AddNodesToGroup(Group, Points, PointCount) 1428 TYPE(CrevasseGroup3D_t), POINTER :: Group 1429 INTEGER :: Points(:) 1430 INTEGER, POINTER :: NewNodeNumbers(:) 1431 INTEGER :: PointCount, NewNumberOfNodes 1432 1433 NewNumberOfNodes = Group % NumberOfNodes + PointCount 1434 ALLOCATE(NewNodeNumbers(NewNumberOfNodes)) 1435 1436 NewNodeNumbers(1:Group % NumberOfNodes) = Group % NodeNumbers 1437 NewNodeNumbers(Group % NumberOfNodes+1:NewNumberOfNodes) = Points(1:PointCount) 1438 1439 !Update count 1440 Group % NumberOfNodes = NewNumberOfNodes 1441 1442 !Point Group to new node list 1443 DEALLOCATE(Group % NodeNumbers) 1444 Group % NodeNumbers => NewNodeNumbers 1445 NULLIFY(NewNodeNumbers) 1446 END SUBROUTINE AddNodesToGroup 1447 1448 !------------------------------------------------------------ 1449 ! Routine to recursively search neighbours and put them 1450 ! in the current group 1451 ! Adapted from 2D Calving 1452 !------------------------------------------------------------ 1453 RECURSIVE SUBROUTINE SearchNeighbours(nodenum, Neighbours, Group, Condition) 1454 INTEGER :: nodenum 1455 INTEGER, POINTER :: Neighbours(:,:) 1456 TYPE(CrevasseGroup3D_t), POINTER :: Group 1457 LOGICAL, ALLOCATABLE :: Condition(:) 1458 !------------------------------------------------ 1459 INTEGER :: i, neighbourindex, NoNeighbours 1460 1461 NoNeighbours = COUNT(Neighbours(nodenum,:) > 0) 1462 DO i = 1,NoNeighbours 1463 neighbourindex = Neighbours(nodenum,i) 1464 IF(.NOT. Condition(neighbourindex)) CYCLE 1465 1466 Group % NumberOfNodes = Group % NumberOfNodes + 1 1467 1468 !check space 1469 IF(Group % NumberOfNodes > SIZE(Group % NodeNumbers)) THEN 1470 PRINT *, 'Debug, need more space, allocating: ', 2*SIZE(Group % NodeNumbers) 1471 CALL DoubleIntVectorSize(Group % NodeNumbers) 1472 PRINT *, 'Debug, new size: ', SIZE(Group % NodeNumbers) 1473 END IF 1474 1475 Group % NodeNumbers(Group % NumberOfNodes) = neighbourindex 1476 1477 !Switch it off so it doesn't get readded 1478 Condition(neighbourindex) = .FALSE. 1479 1480 CALL SearchNeighbours(neighbourindex, Neighbours, Group, Condition) 1481 END DO 1482 1483 END SUBROUTINE SearchNeighbours 1484 1485 !Marks recursive neighbours with same int 1486 RECURSIVE SUBROUTINE MarkNeighbours(nodenum, Neighbours, Array, Mark) 1487 INTEGER :: nodenum 1488 INTEGER, POINTER :: Array(:) 1489 LOGICAL, POINTER :: Neighbours(:,:) 1490 !------------------------------------------------ 1491 INTEGER :: i, Mark 1492 1493 DO i = 1,SIZE(Neighbours,1) 1494 IF(.NOT. Neighbours(nodenum,i)) CYCLE 1495 IF(Array(i)==Mark) CYCLE !already got 1496 1497 Array(i) = Mark 1498 CALL MarkNeighbours(i, Neighbours, Array, Mark) 1499 END DO 1500 1501 END SUBROUTINE MarkNeighbours 1502 1503 !------------------------------------------------------------- 1504 ! Given a CrevasseGroup3D object, finds and stores boundary nodes 1505 ! BoundaryMask is a logical array TRUE where node sits on a 1506 ! mesh (not group) boundary 1507 ! Note: Not used 1508 !------------------------------------------------------------- 1509 SUBROUTINE GetGroupBoundaryNodes(Group, Neighbours, BoundaryMask) 1510 TYPE(CrevasseGroup3D_t), POINTER :: Group 1511 INTEGER, POINTER :: Neighbours(:,:) 1512 LOGICAL :: BoundaryMask(:) 1513 !----------------------------------------- 1514 INTEGER :: i, j, node, BNodes, NoNeighbours, neighbour 1515 INTEGER, ALLOCATABLE :: WorkInt(:) 1516 LOGICAL :: IsBoundaryNode 1517 1518 IF(ASSOCIATED(Group % BoundaryNodes)) & 1519 DEALLOCATE(Group % BoundaryNodes) 1520 1521 ALLOCATE(Group % BoundaryNodes(100)) 1522 Group % BoundaryNodes = 0 1523 BNodes = 0 1524 1525 DO i=1, Group % NumberOfNodes 1526 IsBoundaryNode = .FALSE. 1527 node = Group % NodeNumbers(i) 1528 1529 IF(BoundaryMask(node)) THEN 1530 IsBoundaryNode = .TRUE. 1531 ELSE 1532 NoNeighbours = COUNT(Neighbours(node, :) > 0) 1533 DO j=1,NoNeighbours 1534 neighbour = Neighbours(node, j) 1535 IF(ANY(Group % NodeNumbers == neighbour)) CYCLE 1536 1537 !Only get here if there's a node NOT in the group 1538 IsBoundaryNode = .TRUE. 1539 EXIT 1540 END DO 1541 END IF 1542 1543 IF(IsBoundaryNode) THEN 1544 BNodes = BNodes + 1 1545 IF(BNodes > SIZE(Group % BoundaryNodes)) & 1546 CALL DoubleIntVectorSize(Group % BoundaryNodes) 1547 Group % BoundaryNodes(BNodes) = node 1548 END IF 1549 END DO 1550 1551 ALLOCATE(WorkInt(BNodes)) 1552 WorkInt = Group % BoundaryNodes(1:BNodes) 1553 DEALLOCATE(Group % BoundaryNodes) 1554 ALLOCATE(Group % BoundaryNodes(BNodes)) 1555 Group % BoundaryNodes = WorkInt 1556 DEALLOCATE(WorkInt) 1557 1558 !TODO: Order boundary nodes (clockwise?) 1559 END SUBROUTINE GetGroupBoundaryNodes 1560 1561!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1562 ! Function to detect if a given node lies within 1563 ! a 3D crevasse group (physically, not 'graph'ically 1564 ! Note: not used... 1565!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1566 FUNCTION NodeInCrevasseGroup(NodeNumber, Nodes, CrevasseGroup) RESULT(InGroup) 1567 INTEGER :: NodeNumber 1568 TYPE(Nodes_t) :: Nodes 1569 TYPE(CrevasseGroup3D_t) :: CrevasseGroup 1570!-------------------------------------------------- 1571 LOGICAL :: InGroup 1572 REAL(KIND=dp) :: node_x,node_y, BB(4) 1573 1574 IF(ANY(CrevasseGroup % NodeNumbers == NodeNumber)) & 1575 CALL Fatal("NodeInCrevasseGroup", "Scanning for node which& 1576 &belongs to CrevasseGroup. This is not intended usage!") 1577 1578 node_x = Nodes % x(NodeNumber) 1579 node_y = Nodes % y(NodeNumber) 1580 1581 BB = CrevasseGroup % BoundingBox 1582 1583 IF(node_x < BB(1) .OR. node_x > BB(2) .OR. & 1584 node_y < BB(3) .OR. node_y > BB(4)) THEN 1585 1586 InGroup = .FALSE. 1587 RETURN 1588 1589 END IF 1590 1591 CALL Fatal("NodeInCrevasseGroup", "Haven't finished implementing this yet!") 1592 1593 !Recursively look at node neighbours, stopping when we reach a group member, 1594 !until we reach freedom (a mesh boundary node) or give up (node is contained 1595 !within crevassegroup, and this tells us about the topology of the group) 1596 1597 !RETURN should not just be a logical, this should be repurposed to inform 1598 !about which boundary it reached. 1599 1600 END FUNCTION NodeInCrevasseGroup 1601 1602 !Doubles the size of a pointer integer array 1603 !This version takes a Pointer argument, should 1604 !be used with care... 1605 SUBROUTINE DoubleIntVectorSizeP(Vec, fill) 1606 INTEGER, POINTER :: Vec(:) 1607 INTEGER, OPTIONAL :: fill 1608 !---------------------------------------- 1609 INTEGER, ALLOCATABLE :: WorkVec(:) 1610 1611 ALLOCATE(WorkVec(SIZE(Vec))) 1612 WorkVec = Vec 1613 1614 DEALLOCATE(Vec) 1615 ALLOCATE(Vec(2*SIZE(WorkVec))) 1616 1617 IF(PRESENT(fill)) THEN 1618 Vec = fill 1619 ELSE 1620 Vec = 0 1621 END IF 1622 1623 Vec(1:SIZE(WorkVec)) = WorkVec 1624 1625 END SUBROUTINE DoubleIntVectorSizeP 1626 1627 !Doubles the size of a pointer integer array 1628 !Allocatable array version 1629 SUBROUTINE DoubleIntVectorSizeA(Vec, fill) 1630 INTEGER, ALLOCATABLE :: Vec(:) 1631 INTEGER, OPTIONAL :: fill 1632 !---------------------------------------- 1633 INTEGER, ALLOCATABLE :: WorkVec(:) 1634 1635 ALLOCATE(WorkVec(SIZE(Vec))) 1636 WorkVec = Vec 1637 1638 DEALLOCATE(Vec) 1639 ALLOCATE(Vec(2*SIZE(WorkVec))) 1640 1641 IF(PRESENT(fill)) THEN 1642 Vec = fill 1643 ELSE 1644 Vec = 0 1645 END IF 1646 1647 Vec(1:SIZE(WorkVec)) = WorkVec 1648 1649 END SUBROUTINE DoubleIntVectorSizeA 1650 1651 1652 !----------------------------------------------------------------------------- 1653 !Given a Nodes_t object, removes the nodes specified by RemoveLogical array 1654 !Optionally, user may provide a list of node numbers (NodeNums), from which 1655 !relevant nodes will also be removed 1656 SUBROUTINE RemoveNodes(InNodes, RemoveLogical, NodeNums) 1657 TYPE(Nodes_t) :: InNodes, WorkNodes 1658 LOGICAL, ALLOCATABLE :: RemoveLogical(:) 1659 INTEGER :: i,counter 1660 INTEGER, POINTER, OPTIONAL :: NodeNums(:) 1661 INTEGER, ALLOCATABLE :: WorkNodeNums(:) 1662 1663 WorkNodes % NumberOfNodes = SIZE(InNodes % x) - COUNT(RemoveLogical) 1664 1665 ALLOCATE(WorkNodes % x(WorkNodes % NumberOfNodes),& 1666 WorkNodes % y(WorkNodes % NumberOfNodes),& 1667 WorkNodes % z(WorkNodes % NumberOfNodes)) 1668 IF(PRESENT(NodeNums)) ALLOCATE(WorkNodeNums(WorkNodes % NumberOfNodes)) 1669 1670 counter = 1 1671 DO i=1,InNodes % NumberOfNodes 1672 IF(.NOT. RemoveLogical(i)) THEN 1673 WorkNodes % x(counter) = InNodes % x(i) 1674 WorkNodes % y(counter) = InNodes % y(i) 1675 WorkNodes % z(counter) = InNodes % z(i) 1676 WorkNodeNums(counter) = NodeNums(i) 1677 1678 counter = counter + 1 1679 END IF 1680 END DO 1681 1682 DEALLOCATE(InNodes % x, InNodes % y, InNodes % z ) 1683 ALLOCATE(InNodes % x(WorkNodes % NumberOfNodes), & 1684 InNodes % y(WorkNodes % NumberOfNodes), & 1685 InNodes % z(WorkNodes % NumberOfNodes)) 1686 1687 IF(PRESENT(NodeNums)) THEN 1688 DEALLOCATE(NodeNums) 1689 ALLOCATE(NodeNums(WorkNodes % NumberOfNodes)) 1690 END IF 1691 1692 InNodes % NumberOfNodes = WorkNodes % NumberOfNodes 1693 InNodes % x = WorkNodes % x 1694 InNodes % y = WorkNodes % y 1695 InNodes % z = WorkNodes % z 1696 NodeNums = WorkNodeNums 1697 1698 DEALLOCATE(WorkNodes % x, WorkNodes % y, WorkNodes % z) 1699 IF(PRESENT(NodeNums)) DEALLOCATE(WorkNodeNums) 1700 END SUBROUTINE RemoveNodes 1701 1702 !------------------------------------------------------------------------------ 1703 !> Sort an index array, and change the order of an real array accordingly. 1704 !> Stolen from GeneralUtils, modified so as to leave the initial index array in order 1705 !------------------------------------------------------------------------------ 1706 SUBROUTINE MySortF( n,c,b ) 1707 !------------------------------------------------------------------------------ 1708 INTEGER :: n,c(:) 1709 INTEGER, ALLOCATABLE :: a(:) 1710 REAL(KIND=dp) :: b(:) 1711 !------------------------------------------------------------------------------ 1712 1713 INTEGER :: i,j,l,ir,ra 1714 REAL(KIND=dp) :: rb 1715 !------------------------------------------------------------------------------ 1716 1717 ALLOCATE(a(SIZE(c))) 1718 a = c 1719 1720 IF ( n <= 1 ) RETURN 1721 1722 l = n / 2 + 1 1723 ir = n 1724 DO WHILE( .TRUE. ) 1725 1726 IF ( l > 1 ) THEN 1727 l = l - 1 1728 ra = a(l) 1729 rb = b(l) 1730 ELSE 1731 ra = a(ir) 1732 rb = b(ir) 1733 a(ir) = a(1) 1734 b(ir) = b(1) 1735 ir = ir - 1 1736 IF ( ir == 1 ) THEN 1737 a(1) = ra 1738 b(1) = rb 1739 RETURN 1740 END IF 1741 END IF 1742 i = l 1743 j = l + l 1744 DO WHILE( j <= ir ) 1745 IF ( j<ir ) THEN 1746 IF ( a(j)<a(j+1) ) j = j+1 1747 END IF 1748 IF ( ra<a(j) ) THEN 1749 a(i) = a(j) 1750 b(i) = b(j) 1751 i = j 1752 j = j + i 1753 ELSE 1754 j = ir + 1 1755 END IF 1756 a(i) = ra 1757 b(i) = rb 1758 END DO 1759 END DO 1760 1761 DEALLOCATE(a) 1762 1763 !------------------------------------------------------------------------------ 1764 END SUBROUTINE MySortF 1765 !------------------------------------------------------------------------------ 1766 1767 1768 !If EdgeMaskName is not provided, returns the ring of nodes which define the extent 1769 !of the upper surface of the mesh, arbitrarily beginning with the nodes from the lowest 1770 !partition (PE). 1771 !If EdgeMaskName is provided, this specifies a lateral margin. Then this returns an 1772 !ordered list of nodenumbers which specify an edge of a domain, 1773 !where the edge is determined by the overlap between the two provided permutations 1774 !NOTE: Returned domain edge is valid only on boss partition (PE=0) 1775 SUBROUTINE GetDomainEdge(Model, Mesh, TopPerm, OrderedNodes, OrderedNodeNums, Parallel, & 1776 EdgeMaskName, Simplify, MinDist) 1777 1778 IMPLICIT NONE 1779 1780 TYPE(Model_t) :: Model 1781 TYPE(Mesh_t), POINTER :: Mesh 1782 INTEGER, POINTER :: TopPerm(:) 1783 TYPE(Nodes_t) :: OrderedNodes, UnorderedNodes 1784 LOGICAL :: Parallel 1785 CHARACTER(MAX_NAME_LEN), OPTIONAL :: EdgeMaskName 1786 LOGICAL, OPTIONAL :: Simplify 1787 REAL(KIND=dp), OPTIONAL :: MinDist 1788 !---------------------------------------------------------------- 1789 TYPE(Element_t), POINTER :: Element 1790 TYPE(NeighbourList_T), ALLOCATABLE :: PartNeighbourList(:) 1791 INTEGER :: i,j,k,m,n,prev,next,part_start,find_start,find_fin,find_stride,put_start,& 1792 put_fin, counter,NoNodes, NoNodesOnEdge, NoNeighbours, neigh, Segments, TotSegSplits, & 1793 direction, index, segnum, soff, foff, target_nodenum, next_nodenum, EdgeBCtag, GlobalNN 1794 INTEGER :: comm, ierr !MPI stuff 1795 INTEGER, POINTER :: UnorderedNodeNums(:)=>NULL(), OrderedNodeNums(:), & 1796 UOGlobalNodeNums(:)=>NULL(), OrderedGlobalNodeNums(:)=>NULL() 1797 INTEGER, ALLOCATABLE :: NeighbourPartsList(:), PartNodesOnEdge(:), & 1798 disps(:), nodenum_disps(:), PartOrder(:,:), MyCornerNodes(:), MyNeighbourParts(:), & 1799 NewSegStart(:), PartSegments(:), SegStarts_Gather(:), WorkInt(:), NodeNeighbours(:,:), & 1800 GlobalCorners(:), CornerParts(:), PCornerCounts(:) 1801 LOGICAL :: Debug, ActivePart, Boss, Simpl, NotThis, Found, ThisBC, FullBoundary 1802 LOGICAL, ALLOCATABLE :: OnEdge(:), ActivePartList(:), RemoveNode(:), IsCornerNode(:) 1803 REAL(KIND=dp) :: prec 1804 REAL(KIND=dp), ALLOCATABLE :: WorkReal(:,:) 1805 CHARACTER(MAX_NAME_LEN) :: FuncName 1806 1807 TYPE AllocIntList_t 1808 INTEGER, DIMENSION(:), POINTER :: Indices 1809 END TYPE AllocIntList_t 1810 TYPE(AllocIntList_t), ALLOCATABLE :: PartSegStarts(:) 1811 1812 !------------------------------------------------ 1813 ! Change in strategy: 1814 ! 1815 ! Previously, used stiffness matrix to determine connectivity, but 1816 ! this causes problems when multiple nodes on the boundary reside 1817 ! in the same top surface tri element: 1818 ! 1819 ! *===*===*---*===*===* 1820 ! from this one---^\ / 1821 ! * 1822 ! ^-- we want this one 1823 ! 1824 ! Various versions of this issue can occur... 1825 ! 1826 ! SO, instead of using the stiffness matrix, we should 1827 ! check all the boundary elements on the relevant SIDE 1828 ! boundary (e.g. calving front, right sidewall...), 1829 ! looking for elements containing nodes for which the 1830 ! top mask is true. 1831 ! 1832 ! Because of the extruded structure of the mesh, nodes 1833 ! within the same boundary quad will always be neighbours, 1834 ! and each node shall have no more than 2 neighbours. 1835 !---------------------------------------------------- 1836 1837 FuncName = "GetDomainEdge" 1838 Debug = .FALSE. 1839 ActivePart = .TRUE. 1840 1841 NoNodes = SIZE(TopPerm) !total number of nodes in domain/partition 1842 1843 IF(Parallel) THEN 1844 comm = ELMER_COMM_WORLD 1845 Boss = (ParEnv % MyPE == 0) 1846 ELSE 1847 Boss = .TRUE. !only one part in serial, so it's in charge of computation 1848 END IF 1849 1850 IF(Boss .AND. Debug) THEN 1851 PRINT *, '=================================================' 1852 PRINT *, ' Locating domain edge for ',TRIM(EdgeMaskName) 1853 PRINT *, '=================================================' 1854 END IF 1855 1856 IF(PRESENT(Simplify)) THEN 1857 Simpl = Simplify 1858 ELSE 1859 Simpl = .FALSE. 1860 END IF 1861 1862 ALLOCATE(OnEdge(NoNodes), NodeNeighbours(NoNodes,2)) 1863 OnEdge = .FALSE. 1864 NodeNeighbours = -1 1865 1866 FullBoundary = .NOT.(PRESENT(EdgeMaskName)) 1867 IF(.NOT. FullBoundary) THEN 1868 !Find correct BC from logical 1869 DO i=1,Model % NumberOfBCs 1870 ThisBC = ListGetLogical(Model % BCs(i) % Values,EdgeMaskName,Found) 1871 IF((.NOT. Found) .OR. (.NOT. ThisBC)) CYCLE 1872 EdgeBCtag = Model % BCs(i) % Tag 1873 EXIT 1874 END DO 1875 END IF 1876 1877 !Cycle boundary elements, marking nodes on edge and finding neighbours 1878 DO i=Mesh % NumberOfBulkElements+1, & 1879 Mesh % NumberOfBulkElements+Mesh % NumberOfBoundaryElements 1880 Element => Mesh % Elements(i) 1881 1882 IF((.NOT. FullBoundary) .AND. Element % BoundaryInfo % Constraint /= EdgeBCtag) & 1883 CYCLE !elem not on lateral boundary 1884 1885 IF(ALL(TopPerm(Element % NodeIndexes) > 0)) CYCLE !not a lateral element 1886 IF(.NOT. ANY(TopPerm(Element % NodeIndexes) > 0)) CYCLE !elem contains no nodes on top 1887 !Logic gates above should leave only lateral elements with some nodes on top. 1888 1889 IF(GetElementFamily(Element) == 1) & 1890 CALL Fatal(FuncName, "101 Elements are supposed to be a thing of the past!") 1891 1892 !Cycle nodes in element 1893 DO j=1,Element % TYPE % NumberOfNodes 1894 IF(.NOT. TopPerm(Element % NodeIndexes(j)) > 0) CYCLE 1895 OnEdge(Element % NodeIndexes(j)) = .TRUE. 1896 1897 !Cycle nodes in element 1898 DO k=1,Element % TYPE % NumberOfNodes 1899 IF(j==k) CYCLE 1900 IF(.NOT. TopPerm(Element % NodeIndexes(k))>0) CYCLE 1901 DO m=1,2 !fill NodeNeighbours 1902 IF(NodeNeighbours(Element % NodeIndexes(j),m) /= -1) CYCLE 1903 NodeNeighbours(Element % NodeIndexes(j),m) = Element % NodeIndexes(k) 1904 EXIT 1905 END DO 1906 IF(.NOT. ANY(NodeNeighbours(Element % NodeIndexes(j),:) == Element % NodeIndexes(k))) & 1907 CALL Fatal(FuncName,'Identified more than two neighbours') 1908 END DO 1909 END DO 1910 1911 END DO 1912 1913 NoNodesOnEdge = COUNT(OnEdge) 1914 IF(NoNodesOnEdge == 1) THEN 1915 CALL Fatal(FuncName, "A single node identified on boundary, should not be possible. & 1916 &Someone is messing around with 101 elements.") 1917 END IF 1918 1919 ALLOCATE(UnorderedNodeNums(NoNodesOnEdge),& 1920 OrderedNodeNums(NoNodesOnEdge)) 1921 OrderedNodeNums = -1 !initialize to invalid value 1922 1923 j = 0 1924 DO i=1,NoNodes 1925 IF(.NOT. OnEdge(i)) CYCLE 1926 j = j + 1 1927 UnorderedNodeNums(j) = i 1928 END DO 1929 1930 !Cycle nodes on edge, looking for one with only one neighbour (a corner) 1931 !Edge case = serial fullboundary run, no corner exists, choose arbitrarily 1932 !Rare case (not dealt with!! TODO) = parallel fullboundary, no corners 1933 ! (whole mesh edge in one partition) 1934 IF(NoNodesOnEdge > 1) THEN 1935 1936 ALLOCATE(IsCornerNode(NoNodesOnEdge)) 1937 IsCornerNode = .FALSE. 1938 1939 DO i=1,NoNodesOnEdge 1940 IsCornerNode(i) = COUNT(NodeNeighbours(UnOrderedNodeNums(i),:) == -1) == 1 1941 IF(COUNT(NodeNeighbours(UnOrderedNodeNums(i),:) == -1) == 2) & 1942 CALL Fatal(FuncName, "Found an isolated node on edge") 1943 END DO 1944 1945 IF(MOD(COUNT(IsCornerNode),2) /= 0) THEN 1946 WRITE(Message,'(A,i0)') "Found an odd number of& 1947 & corner nodes in partition: ",ParEnv % MyPE 1948 CALL Fatal(FuncName, Message) 1949 END IF 1950 1951 IF(FullBoundary .AND. .NOT. Parallel) THEN 1952 1953 !If serial FullBoundary request, no corner exists so just choose the first 1954 !unordered node in the list and loop from there 1955 Segments = 1 1956 ALLOCATE(MyCornerNodes(2)) 1957 MyCornerNodes(1) = 1 1958 1959 ELSE 1960 1961 Segments = COUNT(IsCornerNode) / 2 1962 IF(Debug .AND. Segments > 1) PRINT *, & 1963 'Partition ',ParEnv % MyPE, ' has ',Segments,' line segments on boundary.' 1964 1965 ALLOCATE(NewSegStart(Segments-1)) 1966 ALLOCATE(MyCornerNodes(COUNT(IsCornerNode))) 1967 1968 counter = 1 1969 DO i=1,NoNodesOnEdge 1970 IF(IsCornerNode(i)) THEN 1971 MyCornerNodes(counter) = i 1972 counter = counter + 1 1973 END IF 1974 END DO 1975 1976 END IF 1977 1978 counter = 1 1979 DO k=1,Segments 1980 1981 IF(k==1) THEN 1982 OrderedNodeNums(counter) = UnorderedNodeNums(MyCornerNodes(1)) 1983 ELSE 1984 DO i=2, SIZE(MyCornerNodes) 1985 IF(ANY(OrderedNodeNums == UnorderedNodeNums(MyCornerNodes(i)))) THEN 1986 CYCLE 1987 ELSE 1988 OrderedNodeNums(counter) = UnorderedNodeNums(MyCornerNodes(i)) 1989 EXIT 1990 END IF 1991 END DO 1992 END IF 1993 counter = counter + 1 1994 1995 !---------------------------------------------------- 1996 ! Move along from corner, filling in order 1997 !---------------------------------------------------- 1998 DO i=counter,NoNodesOnEdge 1999 Found = .FALSE. 2000 IF(OrderedNodeNums(i-1) == -1) CALL Abort() 2001 2002 DO j=1,2 2003 IF(NodeNeighbours(OrderedNodeNums(i-1),j) == -1) CYCLE !First and last nodes, corner 2004 IF(ANY(OrderedNodeNums(1:i-1) == NodeNeighbours(OrderedNodeNums(i-1),j))) & 2005 CYCLE !already in list 2006 2007 OrderedNodeNums(i) = NodeNeighbours(OrderedNodeNums(i-1),j) 2008 Found = .TRUE. 2009 END DO 2010 2011 IF(.NOT. Found) EXIT 2012 END DO 2013 2014 counter = i 2015 2016 IF(counter >= NoNodesOnEdge) EXIT !this should be redundant... 2017 NewSegStart(k) = counter 2018 END DO 2019 2020 ELSE !Either 1 or 0 nodes found, not an active boundary partition 2021 !0 node case, obvious 2022 !1 node case, if THIS partition only has one node on the boundary, 2023 !this same node must be caught by two other partitions, so we aren't needed. 2024 ALLOCATE(NewSegStart(0), MyCornerNodes(0)) 2025 ActivePart = .FALSE. 2026 Segments = 0 2027 NoNodesOnEdge = 0 !simplifies mpi comms later 2028 IF(.NOT.Parallel) CALL Fatal(FuncName,& 2029 "Found either 1 or 0 nodes in a serial run, this isn't a valid boundary edge!") 2030 END IF 2031 2032 2033 !Remember that, in parallel, we're using local rather than global node numbers 2034 IF(Parallel) THEN 2035 2036 !gather corner count - replaces 101 element detection 2037 ALLOCATE(PCornerCounts(ParEnv % PEs),disps(ParEnv % PEs)) 2038 2039 CALL MPI_AllGather(SIZE(MyCornerNodes), 1, MPI_INTEGER, PCornerCounts, & 2040 1, MPI_INTEGER, ELMER_COMM_WORLD, ierr) 2041 2042 disps(1) = 0 2043 DO i=2, ParEnv % PEs 2044 disps(i) = disps(i-1) + PCornerCounts(i-1) 2045 END DO 2046 2047 ALLOCATE(GlobalCorners(SUM(PCornerCounts)),& 2048 CornerParts(SUM(PCornerCounts))) 2049 2050 !gather corner nodenums 2051 CALL MPI_AllGatherV(Mesh % ParallelInfo % GlobalDOFs(UnorderedNodeNums(MyCornerNodes)), & 2052 SIZE(MyCornerNodes), MPI_INTEGER, GlobalCorners, PCornerCounts, disps, & 2053 MPI_INTEGER, ELMER_COMM_WORLD, ierr) 2054 2055 !note which partition sent each corner node 2056 counter = 1 2057 DO i=1,ParEnv % PEs 2058 IF(PCornerCounts(i) == 0) CYCLE 2059 CornerParts(counter:counter+PCornerCounts(i)-1) = i-1 2060 counter = counter + PCornerCounts(i) 2061 END DO 2062 2063 !Quick check: 2064 DO i=1,SIZE(GlobalCorners) 2065 counter = COUNT(GlobalCorners == GlobalCorners(i)) 2066 IF(counter > 2) CALL Fatal(FuncName,"Programming error in partition & 2067 &segment detection, node found too many times!") 2068 END DO 2069 !Now GlobalCorners and CornerParts tell us which partitions found corner nodes 2070 !(i.e. nodes which will join other segments) 2071 2072 IF(ActivePart) THEN 2073 ALLOCATE(MyNeighbourParts(Segments*2)) 2074 2075 DO i=1,Segments*2 !Find neighbour partition numbers 2076 2077 IF(i==1) THEN 2078 n = OrderedNodeNums(1) 2079 ELSE IF(i==Segments*2) THEN 2080 n = OrderedNodeNums(NoNodesOnEdge) 2081 ELSE IF(MOD(i,2)==0) THEN 2082 n = OrderedNodeNums(NewSegStart(i/2)-1) 2083 ELSE 2084 n = OrderedNodeNums(NewSegStart(i/2)) 2085 END IF 2086 2087 MyNeighbourParts(i) = -1 !default if not caught in loop below 2088 GlobalNN = Mesh % ParallelInfo % GlobalDOFs(n) 2089 DO j=1,SIZE(GlobalCorners) 2090 IF(GlobalCorners(j) /= GlobalNN) CYCLE 2091 IF(CornerParts(j) == ParEnv % MyPE) CYCLE 2092 MyNeighbourParts(i) = CornerParts(j) 2093 IF( .NOT. (ANY(Mesh % ParallelInfo % NeighbourList(n) % Neighbours & 2094 == MyNeighbourParts(i)))) CALL Fatal(FuncName, & 2095 "Failed sanity check on neighbour partition detection.") 2096 END DO 2097 END DO 2098 ELSE 2099 ALLOCATE(MyNeighbourParts(0)) 2100 END IF 2101 2102 IF(Boss) ALLOCATE(PartSegments(ParEnv % PEs)) 2103 2104 CALL MPI_GATHER(Segments, 1, MPI_INTEGER, PartSegments, & 2105 1, MPI_INTEGER, 0, comm, ierr) 2106 2107 IF(Boss) THEN 2108 2109 TotSegSplits = 0 2110 DO i=1,SIZE(PartSegments) 2111 TotSegSplits = TotSegSplits + MAX(PartSegments(i)-1,0) 2112 END DO 2113 2114 ALLOCATE(nodenum_disps(ParEnv % PEs), & 2115 PartNodesOnEdge(ParEnv % PEs), & 2116 NeighbourPartsList(SUM(PartSegments)*2), & 2117 PartNeighbourList(ParEnv % PEs), & 2118 SegStarts_Gather(TotSegSplits)) 2119 2120 DO i=1,ParEnv % PEs 2121 ALLOCATE(PartNeighbourList(i) % Neighbours(PartSegments(i)*2)) 2122 END DO 2123 2124 disps(1) = 0 2125 DO i=2, ParEnv % PEs 2126 disps(i) = disps(i-1) + MAX(PartSegments(i-1)-1,0) 2127 END DO 2128 2129 END IF 2130 2131 !Get found count from each part to boss 2132 CALL MPI_GATHER(NoNodesOnEdge, 1, MPI_INTEGER, PartNodesOnEdge, & 2133 1, MPI_INTEGER, 0, comm ,ierr) 2134 IF(ierr /= MPI_SUCCESS) CALL Fatal(FuncName,"MPI Error!") 2135 CALL MPI_BARRIER(ELMER_COMM_WORLD, ierr) 2136 2137 IF(Debug .AND. Boss) THEN 2138 PRINT *, 'boss size(SegStarts_Gather): ', SIZE(SegStarts_Gather) 2139 PRINT *, 'boss PartSegments: ', PartSegments 2140 PRINT *, 'boss disps:', disps 2141 DO i=1,ParEnv % PEs 2142 IF(PartNodesOnEdge(i) == 0) CYCLE 2143 PRINT *, 'partition ',i-1,' NoNodesOnEdge: ',PartNodesOnEdge(i) 2144 END DO 2145 END IF 2146 2147 IF(Boss) THEN 2148 ALLOCATE(WorkInt(ParEnv % PEs)) 2149 WorkInt = MAX(PartSegments-1,0) 2150 END IF 2151 2152 CALL MPI_GATHERV(NewSegStart, MAX(Segments-1,0), MPI_INTEGER, SegStarts_Gather, & 2153 WorkInt, disps, MPI_INTEGER, 0, comm, ierr) 2154 IF(ierr /= MPI_SUCCESS) CALL Fatal(FuncName,"MPI Error!") 2155 CALL MPI_BARRIER(ELMER_COMM_WORLD, ierr) 2156 2157 IF(Boss) THEN 2158 ALLOCATE(PartSegStarts(ParEnv % PEs)) 2159 DO i=1,ParEnv % PEs 2160 j = PartSegments(i) 2161 ALLOCATE( PartSegStarts(i) % Indices(MAX((j - 1),0))) 2162 IF(j > 1) THEN 2163 IF(Debug) PRINT *, 'debug disps(i),j', disps(i),j 2164 PartSegStarts(i) % Indices = SegStarts_Gather(1+disps(i) : (1+disps(i) + (j-1)-1) ) 2165 END IF 2166 IF(Debug) PRINT *, i,' partsegstarts: ', PartSegStarts(i) % Indices 2167 END DO 2168 2169 disps(1) = 0 2170 DO i=2, ParEnv % PEs 2171 disps(i) = disps(i-1) + PartSegments(i-1)*2 2172 END DO 2173 2174 WorkInt = PartSegments*2 2175 END IF 2176 2177 !Get neighbour part numbers from each part to boss 2178 CALL MPI_GATHERV(MyNeighbourParts, Segments*2, MPI_INTEGER, NeighbourPartsList, & 2179 WorkInt, disps, MPI_INTEGER, 0, comm ,ierr) 2180 IF(ierr /= MPI_SUCCESS) CALL Fatal(FuncName,"MPI Error!") 2181 CALL MPI_BARRIER(ELMER_COMM_WORLD, ierr) 2182 2183 IF(Debug .AND. Boss) PRINT *, 'DEBUG, NewSegStart: ', NewSegStart 2184 2185 IF(Boss) THEN 2186 ActivePartList = (PartNodesOnEdge > 0) 2187 2188 !Here we account for shared nodes on partition boundaries 2189 OrderedNodes % NumberOfNodes = SUM(PartNodesOnEdge) - (SIZE(NeighbourPartsList)/2 - 1) 2190 !but they are still present when gathered... 2191 UnorderedNodes % NumberOfNodes = SUM(PartNodesOnEdge) 2192 2193 ALLOCATE(PartOrder(SIZE(NeighbourPartsList)/2,2),& 2194 OrderedNodes % x(OrderedNodes % NumberOfNodes),& 2195 OrderedNodes % y(OrderedNodes % NumberOfNodes),& 2196 OrderedNodes % z(OrderedNodes % NumberOfNodes),& 2197 UnorderedNodes % x(UnorderedNodes % NumberOfNodes),& 2198 UnorderedNodes % y(UnorderedNodes % NumberOfNodes),& 2199 UnorderedNodes % z(UnorderedNodes % NumberOfNodes),& 2200 UOGlobalNodeNums(UnorderedNodes % NumberOfNodes),& 2201 OrderedGlobalNodeNums(OrderedNodes % NumberOfNodes)) 2202 2203 nodenum_disps(1) = 0 2204 DO i=2, ParEnv % PEs 2205 nodenum_disps(i) = nodenum_disps(i-1) + PartNodesOnEdge(i-1) 2206 END DO 2207 2208 IF(Debug) THEN 2209 PRINT *, 'debug disps: ', disps 2210 PRINT *, 'debug nodenum_disps: ', nodenum_disps 2211 PRINT *, 'debug neighbourpartslist: ',NeighbourPartsList 2212 PRINT *, 'Partition Segments: ',PartSegments 2213 END IF 2214 END IF 2215 2216 !----------------------------------------------------------- 2217 ! Gather node coords from all partitions 2218 ! Note, they're going into 'UnorderedNodes': though they are ordered 2219 ! within their partition, the partitions aren't ordered... 2220 !----------------------------------------------------------- 2221 2222 !Global Node Numbers 2223 CALL MPI_GATHERV(Mesh % ParallelInfo % GlobalDOFs(OrderedNodeNums),& 2224 NoNodesOnEdge,MPI_INTEGER,& 2225 UOGlobalNodeNums,PartNodesOnEdge,& 2226 nodenum_disps,MPI_INTEGER,0,comm, ierr) 2227 IF(ierr /= MPI_SUCCESS) CALL Fatal(FuncName,"MPI Error!") 2228 CALL MPI_BARRIER(ELMER_COMM_WORLD, ierr) 2229 2230 !X coords 2231 CALL MPI_GATHERV(Mesh % Nodes % x(OrderedNodeNums),& 2232 NoNodesOnEdge,MPI_DOUBLE_PRECISION,& 2233 UnorderedNodes % x,PartNodesOnEdge,& 2234 nodenum_disps,MPI_DOUBLE_PRECISION,0,comm, ierr) 2235 IF(ierr /= MPI_SUCCESS) CALL Fatal(FuncName,"MPI Error!") 2236 CALL MPI_BARRIER(ELMER_COMM_WORLD, ierr) 2237 2238 !Y coords 2239 CALL MPI_GATHERV(Mesh % Nodes % y(OrderedNodeNums),& 2240 NoNodesOnEdge,MPI_DOUBLE_PRECISION,& 2241 UnorderedNodes % y,PartNodesOnEdge,& 2242 nodenum_disps,MPI_DOUBLE_PRECISION,0,comm, ierr) 2243 IF(ierr /= MPI_SUCCESS) CALL Fatal(FuncName,"MPI Error!") 2244 CALL MPI_BARRIER(ELMER_COMM_WORLD, ierr) 2245 2246 !Z coords 2247 CALL MPI_GATHERV(Mesh % Nodes % z(OrderedNodeNums),& 2248 NoNodesOnEdge,MPI_DOUBLE_PRECISION,& 2249 UnorderedNodes % z,PartNodesOnEdge,& 2250 nodenum_disps,MPI_DOUBLE_PRECISION,0,comm, ierr) 2251 IF(ierr /= MPI_SUCCESS) CALL Fatal(FuncName,"MPI Error!") 2252 CALL MPI_BARRIER(ELMER_COMM_WORLD, ierr) 2253 2254 !----------------------------------------------------------- 2255 ! Determine order of partitions by linking neighbours and 2256 ! checking globalnodenumbers where appropriate 2257 !----------------------------------------------------------- 2258 2259 IF(Boss) THEN 2260 !Notes: NeighbourPartsList is zero indexed, like PEs 2261 !PartOrder is 1 indexed 2262 !disps is 1 indexed. So disps(NeighbourPartsList+1) 2263 2264 PartOrder = 0 !init 2265 direction = 0 2266 prev = -1 2267 next = 0 2268 2269 !First fill in PartNeighbourList % Neighbours 2270 DO i=1,ParEnv % PEs 2271 IF(PartSegments(i)==0) CYCLE 2272 PartNeighbourList(i) % Neighbours = & 2273 NeighbourPartsList( (1+disps(i)) : (1+disps(i) + (PartSegments(i)*2) - 1) ) 2274 !There is the possibility of missing an end (-1) due to partition 2275 !landing right on corner 2276 DO j=1,SIZE(PartNeighbourList(i) % Neighbours) 2277 IF(PartNeighbourList(i) % Neighbours(j) == -1) CYCLE 2278 IF(PartSegments(PartNeighbourList(i) % Neighbours(j)+1) < 1) THEN 2279 IF(Debug) PRINT *, 'Neighbour ',PartNeighbourList(i) % Neighbours(j)+1,& 2280 "isn't really on boundary, so changing to -1" 2281 PartNeighbourList(i) % Neighbours(j) = -1 2282 END IF 2283 END DO 2284 2285 IF(Debug) PRINT *, i-1, ': Neighbours: ', PartNeighbourList(i) % Neighbours 2286 !find a corner partition 2287 IF(ANY(PartNeighbourList(i) % Neighbours == prev)) next = i 2288 END DO 2289 2290 !No partition had corner (-1) 2291 IF(next==0) THEN 2292 IF(FullBoundary) THEN !this is expected, a closed loop so no -1 2293 DO i=1,ParEnv % PEs 2294 IF(PartSegments(i)>0) THEN 2295 next = i 2296 prev = PartNeighbourList(i) % Neighbours(1) 2297 EXIT 2298 END IF 2299 END DO 2300 ELSE 2301 CALL Fatal(FuncName,"Error finding corner of requested boundary in partitions.") 2302 END IF 2303 ELSE IF(FullBoundary) THEN 2304 CALL Fatal(FuncName,"Error - found corner but requested FullBoundary& 2305 &- programming mistake.") 2306 END IF 2307 2308 IF(Debug) THEN 2309 PRINT *, 'Debug GetDomainEdge, globalno, unorderednodes % x: ' 2310 DO i=1,SIZE(UOGlobalNodeNums) 2311 PRINT *, i, UOGlobalNodeNums(i), UnorderedNodes % x(i) 2312 END DO 2313 2314 PRINT *, 'debug nodenum_disps: ' 2315 DO i=1, SIZE(nodenum_disps) 2316 PRINT *, i,' ',nodenum_disps(i) 2317 END DO 2318 END IF 2319 2320 2321 counter = 1 2322 2323 DO WHILE(.TRUE.) 2324 IF(Debug) PRINT *,'Next Partition is: ',next 2325 IF((COUNT(PartNeighbourList(next) % Neighbours == prev) == 1) .OR. & 2326 (prev == -1)) THEN 2327 DO j=1,SIZE(PartNeighbourList(next) % Neighbours) 2328 IF(PartNeighbourList(next) % Neighbours(j) == prev) THEN 2329 index = j 2330 EXIT 2331 END IF 2332 END DO 2333 ELSE !Neighbours on both sides, so need to inspect globalnodenumbers 2334 IF(Debug) PRINT *, 'debug, two matches' 2335 DO j=1,SIZE(PartNeighbourList(next) % Neighbours) 2336 IF(PartNeighbourList(next) % Neighbours(j) == prev) THEN 2337 2338 segnum = ((j-1)/2) + 1 2339 direction = (2 * MOD(j, 2)) - 1 2340 2341 IF(segnum == 1) THEN 2342 soff = 0 2343 ELSE 2344 soff = PartSegStarts(next) % Indices(segnum - 1) - 1 2345 END IF 2346 IF(segnum == PartSegments(next)) THEN 2347 foff = 0 2348 ELSE 2349 foff = -1 * (PartNodesOnEdge(next) - PartSegStarts(next) % Indices(segnum) + 1) 2350 END IF 2351 2352 IF(direction > 0) THEN 2353 next_nodenum = UOGlobalNodeNums(1 + nodenum_disps(next) + soff) 2354 ELSE 2355 !one node before (-1) the next partition's (+1) nodes 2356 IF(next == ParEnv % PEs) THEN 2357 k = SIZE(UOGlobalNodeNums) 2358 ELSE 2359 k = 1 + nodenum_disps(next+1) - 1 2360 END IF 2361 next_nodenum = UOGlobalNodeNums(k + foff) 2362 END IF 2363 IF(Debug) THEN 2364 PRINT *, 'debug, next_nodenum: ', next_nodenum 2365 PRINT *, 'debug, target_nodenum: ', target_nodenum 2366 END IF 2367 IF(next_nodenum == target_nodenum) THEN 2368 index = j 2369 EXIT 2370 END IF 2371 END IF 2372 END DO 2373 END IF 2374 2375 segnum = ((index-1)/2) + 1 !1,2 -> 1, 3,4 -> 2 2376 direction = (2 * MOD(index, 2)) - 1 2377 PartOrder(counter,1) = next - 1 2378 PartOrder(counter,2) = direction * segnum 2379 counter = counter + 1 2380 2381 IF(Debug) THEN 2382 PRINT *, 'index: ', index 2383 PRINT *, 'segnum: ', segnum 2384 PRINT *, 'direction: ',direction 2385 PRINT *, 'next: ', next 2386 PRINT *, 'prev: ', prev 2387 END IF 2388 2389 prev = next - 1 2390 j = next 2391 next = PartNeighbourList(next) % Neighbours(index + direction) 2392 2393 !In case of two matches, need a target node to find 2394 IF(segnum == 1) THEN 2395 soff = 0 2396 ELSE 2397 soff = PartSegStarts(j) % Indices(segnum - 1) - 1 2398 END IF 2399 IF(segnum == PartSegments(j)) THEN 2400 foff = 0 2401 ELSE 2402 foff = -1 * (PartNodesOnEdge(j) - PartSegStarts(j) % Indices(segnum) + 1) 2403 END IF 2404 2405 IF(direction < 0) THEN 2406 target_nodenum = UOGlobalNodeNums(1 + nodenum_disps(prev+1) + soff) 2407 ELSE 2408 IF(prev + 1 == ParEnv % PEs) THEN 2409 k = SIZE(UOGlobalNodeNums) 2410 ELSE 2411 k = 1 + nodenum_disps(prev+1+1) - 1 2412 END IF 2413 !one node before (-1) the next partition's (+1) nodes 2414 target_nodenum = UOGlobalNodeNums(k + foff) 2415 END IF 2416 2417 !wipe them out so we don't accidentally come back this way 2418 PartNeighbourList(j) % Neighbours(index:index+direction:direction) = -2 2419 2420 IF(FullBoundary) THEN 2421 IF(Debug) THEN 2422 PRINT *, 'new index: ', index 2423 PRINT *, 'new segnum: ', segnum 2424 PRINT *, 'new direction: ',direction 2425 PRINT *, 'new next: ', next 2426 PRINT *, 'new prev: ', prev 2427 PRINT *, 'new neighbours: ', PartNeighbourList(next+1) % Neighbours 2428 END IF 2429 2430 IF(ALL(PartNeighbourList(next+1) % Neighbours == -2)) THEN 2431 IF(Debug) PRINT *,'Finished cycling neighbours in FullBoundary' 2432 EXIT 2433 END IF 2434 ELSE IF(next == -1) THEN 2435 EXIT 2436 END IF 2437 2438 next = next + 1 2439 END DO 2440 2441 IF(Debug) PRINT *, 'Debug GetDomainEdge, part order:', PartOrder 2442 2443 END IF 2444 2445 !----------------------------------------------------------- 2446 ! Put nodes collected from partitions into order 2447 !----------------------------------------------------------- 2448 2449 IF(Boss) THEN 2450 put_start = 1 2451 2452 DO i=1,SIZE(PartOrder,1) 2453 j = PartOrder(i,1) + 1 2454 segnum = PartOrder(i,2) 2455 2456 IF(j==0) CALL Abort() 2457 2458 foff = 0 2459 soff = 0 2460 IF(PartSegments(j) > 1) THEN 2461 IF(Debug) THEN 2462 PRINT *, 'Debug GetDomainEdge, extracting nodes from segmented partition' 2463 PRINT *, 'Debug GetDomainEdge, segnum: ', segnum 2464 PRINT *, 'Debug GetDomainEdge, partnodes: ', PartNodesOnEdge(j) 2465 PRINT *, 'Debug GetDomainEdge, PartSegStarts(j) % Indices: ',& 2466 PartSegStarts(j) % Indices 2467 PRINT *, 'Debug GetDomainEdge, nodenum_disps(j): ',nodenum_disps(j) 2468 END IF 2469 2470 IF(ABS(segnum) == 1) THEN 2471 2472 soff = 0 2473 ELSE 2474 soff = PartSegStarts(j) % Indices(ABS(segnum) - 1) - 1 2475 END IF 2476 IF(ABS(segnum) == PartSegments(j)) THEN 2477 foff = 0 2478 ELSE 2479 foff = -1 * (PartNodesOnEdge(j) - PartSegStarts(j) % Indices(ABS(segnum)) + 1) 2480 END IF 2481 END IF 2482 2483 part_start = 1 + nodenum_disps(j) !where are this partitions nodes? 2484 IF(segnum > 0) THEN 2485 find_start = part_start + soff 2486 find_fin = part_start + PartNodesOnEdge(j) - 1 + foff 2487 find_stride = 1 2488 ELSE 2489 find_fin = part_start + soff 2490 find_start = part_start + PartNodesOnEdge(j) - 1 + foff 2491 find_stride = -1 2492 END IF 2493 2494 put_fin = put_start + ABS(find_start - find_fin) 2495 IF(Debug) THEN 2496 PRINT *, 'Debug, find start, end: ',find_start, find_fin, find_stride 2497 PRINT *, 'Debug, put start, end: ',put_start, put_fin 2498 PRINT *, 'Total slots: ',SIZE(OrderedNodes % x) 2499 END IF 2500 2501 OrderedNodes % x(put_start:put_fin) = & 2502 UnorderedNodes % x(find_start:find_fin:find_stride) 2503 OrderedNodes % y(put_start:put_fin) = & 2504 UnorderedNodes % y(find_start:find_fin:find_stride) 2505 OrderedNodes % z(put_start:put_fin) = & 2506 UnorderedNodes % z(find_start:find_fin:find_stride) 2507 OrderedGlobalNodeNums(put_start:put_fin) = & 2508 UOGlobalNodeNums(find_start:find_fin:find_stride) 2509 2510 put_start = put_fin !1 node overlap 2511 END DO 2512 2513 IF(FullBoundary) THEN 2514 !In the full boundary case, we've inadvertently saved the first node twice 2515 ! (once at the end too) - this sorts that out 2516 n = OrderedNodes % NumberOfNodes - 1 2517 OrderedNodes % NumberOfNodes = n 2518 2519 ALLOCATE(WorkReal(n,3)) 2520 WorkReal(:,1) = OrderedNodes % x(1:n) 2521 WorkReal(:,2) = OrderedNodes % y(1:n) 2522 WorkReal(:,3) = OrderedNodes % z(1:n) 2523 DEALLOCATE(OrderedNodes % x, OrderedNodes % y, OrderedNodes % z) 2524 ALLOCATE(OrderedNodes % x(n), OrderedNodes % y(n), OrderedNodes % z(n)) 2525 OrderedNodes % x(1:n) = WorkReal(:,1) 2526 OrderedNodes % y(1:n) = WorkReal(:,2) 2527 OrderedNodes % z(1:n) = WorkReal(:,3) 2528 DEALLOCATE(WorkReal) 2529 END IF 2530 2531 DEALLOCATE(OrderedNodeNums) 2532 ALLOCATE(OrderedNodeNums(OrderedNodes % NumberOfNodes)) 2533 OrderedNodeNums = OrderedGlobalNodeNums(1:OrderedNodes % NumberOfNodes) 2534 2535 IF(Debug) THEN 2536 PRINT *, 'Debug GetDomainEdge, globalno, orderednodes % x: ' 2537 DO i=1,SIZE(OrderedNodes % x) 2538 PRINT *, OrderedNodeNums(i), OrderedNodes % x(i) 2539 END DO 2540 END IF 2541 END IF 2542 2543 ELSE !serial 2544 OrderedNodes % NumberOfNodes = NoNodesOnEdge 2545 ALLOCATE(OrderedNodes % x(OrderedNodes % NumberOfNodes),& 2546 OrderedNodes % y(OrderedNodes % NumberOfNodes),& 2547 OrderedNodes % z(OrderedNodes % NumberOfNodes)) 2548 2549 OrderedNodes % x = Mesh % Nodes % x(OrderedNodeNums) 2550 OrderedNodes % y = Mesh % Nodes % y(OrderedNodeNums) 2551 OrderedNodes % z = Mesh % Nodes % z(OrderedNodeNums) 2552 2553 !No action required on OrderedNodeNums... 2554 END IF 2555 2556 !------------------------------------------------------------- 2557 ! Simplify geometry by removing interior nodes on any straight 2558 ! lines if requested 2559 !------------------------------------------------------------- 2560 IF(Simpl .AND. Boss) THEN 2561 ALLOCATE(RemoveNode(OrderedNodes % NumberOfNodes)) 2562 RemoveNode = .FALSE. 2563 2564 DO i=2,OrderedNodes % NumberOfNodes-1 !Test all interior nodes 2565 IF(Debug) THEN 2566 PRINT *, (NodesGradXY(OrderedNodes,i,i-1)) 2567 PRINT *, (NodesGradXY(OrderedNodes,i+1,i)) 2568 PRINT *, 'DIFF: ',ABS(NodesGradXY(OrderedNodes,i,i-1) -& 2569 NodesGradXY(OrderedNodes,i+1,i)) 2570 PRINT *, '' 2571 END IF 2572 2573 !Need to determine numerical precision of input datapoints 2574 !i.e. after how many decimal places are values constant 2575 !e.g. 0.23000000... or 99999... 2576 prec = MAX(RealAeps(OrderedNodes % x(i)),RealAeps(OrderedNodes % y(i))) 2577 2578 IF(ABS(NodesGradXY(OrderedNodes,i,i-1) - NodesGradXY(OrderedNodes,i+1,i)) < prec) THEN 2579 RemoveNode(i) = .TRUE. 2580 END IF 2581 END DO 2582 2583 IF(COUNT(RemoveNode) > 0) THEN 2584 2585 CALL RemoveNodes(OrderedNodes, RemoveNode, OrderedNodeNums) 2586 2587 IF(Debug) THEN 2588 PRINT *, 'Debug GetDomainEdge, Simplify removing: ', COUNT(RemoveNode), ' nodes' 2589 DO i=1,OrderedNodes % NumberOfNodes 2590 PRINT *, 'Debug GetDomainEdge, node: ',i 2591 PRINT *, 'x: ',OrderedNodes % x(i),'y: ',OrderedNodes % y(i) 2592 END DO 2593 END IF !debug 2594 2595 END IF !removing any nodes 2596 DEALLOCATE(RemoveNode) 2597 END IF !simplify 2598 2599 !------------------------------------------------------------- 2600 ! Remove any nodes which are closer together than MinDist, if 2601 ! this is specified. 2602 !------------------------------------------------------------- 2603 IF(PRESENT(MinDist) .AND. Boss) THEN 2604 !Cycle all nodes, remove any too close together 2605 !This won't guarantee that the new domain edge is *within* the old one 2606 !but could be adapted to do so 2607 ALLOCATE(RemoveNode(OrderedNodes % NumberOfNodes)) 2608 RemoveNode = .FALSE. 2609 DO i=2,OrderedNodes % NumberOfNodes-1 !Test all interior nodes 2610 j = i - 1 2611 DO WHILE(RemoveNode(j)) 2612 j = j-1 2613 END DO 2614 2615 IF(NodeDist2D(OrderedNodes, i, j) < MinDist) THEN 2616 RemoveNode(i) = .TRUE. 2617 IF(Debug) THEN 2618 PRINT *, 'Debug GetDomainEdge, MinDist, removing node ',i,' too close to: ', j 2619 PRINT *, 'Debug GetDomainEdge, MinDist, dist: ',NodeDist2D(OrderedNodes, i, j) 2620 END IF 2621 END IF 2622 END DO 2623 2624 IF(COUNT(RemoveNode) > 0) THEN 2625 2626 CALL RemoveNodes(OrderedNodes, RemoveNode, OrderedNodeNums) 2627 2628 IF(Debug) THEN 2629 PRINT *, 'Debug GetDomainEdge, MinDist removing: ', COUNT(RemoveNode), ' nodes' 2630 DO i=1,OrderedNodes % NumberOfNodes 2631 PRINT *, 'Debug GetDomainEdge, node: ',i 2632 PRINT *, 'x: ',OrderedNodes % x(i),'y: ',OrderedNodes % y(i) 2633 END DO 2634 END IF !debug 2635 2636 END IF !removing any nodes 2637 DEALLOCATE(RemoveNode) 2638 END IF !MinDist 2639 2640 !------------ DEALLOCATIONS ------------------ 2641 2642 DEALLOCATE(OnEdge, UnorderedNodeNums, GlobalCorners, CornerParts, PCornerCounts) 2643 2644 IF(Boss .AND. Parallel) THEN !Deallocations 2645 DEALLOCATE(UnorderedNodes % x, & 2646 UnorderedNodes % y, & 2647 UnorderedNodes % z, & 2648 PartNodesOnEdge, & 2649 disps, nodenum_disps, & 2650 PartOrder, & 2651 UOGlobalNodeNums, & 2652 OrderedGlobalNodeNums) 2653 END IF 2654 2655 END SUBROUTINE GetDomainEdge 2656 2657 ! Copies over time variables and creates coordinate vars. Basically pinched 2658 ! from AddMeshCoordinatesAndTime() and Multigrid 2659 SUBROUTINE CopyIntrinsicVars(OldMesh, NewMesh) 2660 IMPLICIT NONE 2661 2662 TYPE(Mesh_t), POINTER :: OldMesh, NewMesh 2663 TYPE(Solver_t), POINTER :: Solver 2664 TYPE(Variable_t), POINTER :: WorkVar 2665 !---------------------------------------------------------- 2666 NULLIFY( Solver ) 2667 2668 CALL VariableAdd( NewMesh % Variables, NewMesh,Solver, & 2669 'Coordinate 1',1,NewMesh % Nodes % x ) 2670 2671 CALL VariableAdd(NewMesh % Variables,NewMesh,Solver, & 2672 'Coordinate 2',1,NewMesh % Nodes % y ) 2673 2674 CALL VariableAdd(NewMesh % Variables,NewMesh,Solver, & 2675 'Coordinate 3',1,NewMesh % Nodes % z ) 2676 2677 WorkVar => VariableGet( OldMesh % Variables, 'Time', ThisOnly=.TRUE.) 2678 CALL VariableAdd( NewMesh % Variables, NewMesh, Solver, 'Time', 1, WorkVar % Values ) 2679 2680 WorkVar => VariableGet( OldMesh % Variables, 'Periodic Time', ThisOnly=.TRUE.) 2681 CALL VariableAdd( NewMesh % Variables, NewMesh, Solver, 'Periodic Time', 1, WorkVar % Values ) 2682 2683 WorkVar => VariableGet( OldMesh % Variables, 'Timestep', ThisOnly=.TRUE.) 2684 CALL VariableAdd( NewMesh % Variables, NewMesh, Solver, 'Timestep', 1, WorkVar % Values ) 2685 2686 WorkVar => VariableGet( OldMesh % Variables, 'Timestep size', ThisOnly=.TRUE.) 2687 CALL VariableAdd( NewMesh % Variables, NewMesh, Solver, 'Timestep size', 1, WorkVar % Values ) 2688 2689 WorkVar => VariableGet( OldMesh % Variables, 'Timestep interval', ThisOnly=.TRUE.) 2690 CALL VariableAdd( NewMesh % Variables, NewMesh, Solver, 'Timestep interval', 1, WorkVar % Values ) 2691 2692 WorkVar => VariableGet( OldMesh % Variables, 'Coupled iter', ThisOnly=.TRUE.) 2693 CALL VariableAdd( NewMesh % Variables, NewMesh, Solver, 'Coupled iter', 1, WorkVar % Values ) 2694 2695 WorkVar => VariableGet( OldMesh % Variables, 'Nonlin iter', ThisOnly=.TRUE.) 2696 CALL VariableAdd( NewMesh % Variables, NewMesh, Solver, 'Nonlin iter', 1, WorkVar % Values ) 2697 2698 END SUBROUTINE CopyIntrinsicVars 2699 2700!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 2701 !Function to rotate a mesh by rotationmatrix 2702!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 2703 SUBROUTINE RotateMesh(Mesh, RotationMatrix) 2704 2705 IMPLICIT NONE 2706 2707 TYPE(Mesh_t) :: Mesh 2708 REAL(KIND=dp) :: RotationMatrix(3,3), NodeHolder(3) 2709 INTEGER :: i 2710 2711 DO i=1,Mesh % NumberOfNodes 2712 NodeHolder(1) = Mesh % Nodes % x(i) 2713 NodeHolder(2) = Mesh % Nodes % y(i) 2714 NodeHolder(3) = Mesh % Nodes % z(i) 2715 2716 NodeHolder = MATMUL(RotationMatrix,NodeHolder) 2717 2718 Mesh % Nodes % x(i) = NodeHolder(1) 2719 Mesh % Nodes % y(i) = NodeHolder(2) 2720 Mesh % Nodes % z(i) = NodeHolder(3) 2721 END DO 2722 2723 END SUBROUTINE RotateMesh 2724 2725 SUBROUTINE DeallocateElement(Element) 2726 2727 IMPLICIT NONE 2728 TYPE(Element_t) :: Element 2729 2730 IF ( ASSOCIATED( Element % NodeIndexes ) ) & 2731 DEALLOCATE( Element % NodeIndexes ) 2732 Element % NodeIndexes => NULL() 2733 2734 IF ( ASSOCIATED( Element % EdgeIndexes ) ) & 2735 DEALLOCATE( Element % EdgeIndexes ) 2736 Element % EdgeIndexes => NULL() 2737 2738 IF ( ASSOCIATED( Element % FaceIndexes ) ) & 2739 DEALLOCATE( Element % FaceIndexes ) 2740 Element % FaceIndexes => NULL() 2741 2742 IF ( ASSOCIATED( Element % DGIndexes ) ) & 2743 DEALLOCATE( Element % DGIndexes ) 2744 Element % DGIndexes => NULL() 2745 2746 IF ( ASSOCIATED( Element % BubbleIndexes ) ) & 2747 DEALLOCATE( Element % BubbleIndexes ) 2748 Element % BubbleIndexes => NULL() 2749 2750 IF ( ASSOCIATED( Element % PDefs ) ) & 2751 DEALLOCATE( Element % PDefs ) 2752 Element % PDefs => NULL() 2753 2754 END SUBROUTINE DeallocateElement 2755 2756 !Identify front elements connected to the bed, which are sufficiently horizontal 2757 !to warrant reclassification as basal elements. 2758 !Note, only does elements currently connected to the bed. i.e. one row per dt 2759 !Returns: 2760 ! NewBasalNode(:), LOGICAL true where frontal node becomes basal 2761 ! ExFrontalNode(:), LOGICAL true where a frontal node no longer 2762 ! belongs to its front column (though it may still be on the front...) 2763 ! 2764 ! NOTE, if an error in this subroutine, could be element 2765 ! which sits between 2 NewBasalElems 2766 2767 SUBROUTINE ConvertFrontalToBasal(Model, Mesh, FrontMaskName, BotMaskName, & 2768 ZThresh, NewBasalNode, FoundSome) 2769 2770 TYPE(Model_t) :: Model 2771 TYPE(Mesh_t), POINTER :: Mesh 2772 REAL(KIND=dp) :: ZThresh 2773 LOGICAL :: FoundSome 2774 LOGICAL, POINTER :: NewBasalNode(:), ExFrontalNode(:), NewBasalElem(:) 2775 CHARACTER(MAX_NAME_LEN) :: FrontMaskName, BotMaskName 2776 !------------------------------------------------------- 2777 TYPE(Nodes_t) :: Nodes 2778 TYPE(Solver_t), POINTER :: NullSolver => NULL() 2779 TYPE(Element_t), POINTER :: Element, New303Elements(:,:), WorkElements(:) 2780 INTEGER :: i,j,k,n,dummyint, ierr, FrontBCtag, BasalBCtag, count303, & 2781 CountSharedExFrontal, CountSharedNewBasal, SharedExGlobal(2), & 2782 SharedNewGlobal(2), OldElemCount, NewElemCount 2783 INTEGER, POINTER :: NodeIndexes(:), AllSharedExGlobal(:)=>NULL(), & 2784 AllSharedNewGlobal(:)=>NULL(), FrontPerm(:), BotPerm(:) 2785 REAL(KIND=dp) :: Normal(3) 2786 LOGICAL :: ThisBC, Found, Debug 2787 CHARACTER(MAX_NAME_LEN) :: FuncName 2788 2789 FoundSome = .FALSE. 2790 FuncName = "ConvertFrontalToBasal" 2791 Debug = .FALSE. 2792 2793 n = Mesh % NumberOfNodes 2794 ALLOCATE(NewBasalNode(n),& 2795 ExFrontalNode(n),& 2796 FrontPerm(n),& 2797 BotPerm(n),& 2798 NewBasalElem(Mesh % NumberOfBulkElements+1: & 2799 Mesh % NumberOfBulkElements+Mesh % NumberOfBoundaryElements)) 2800 2801 NewBasalNode = .FALSE. 2802 ExFrontalNode = .FALSE. 2803 NewBasalElem = .FALSE. 2804 2805 CALL MakePermUsingMask( Model, NullSolver, Mesh, BotMaskName, & 2806 .FALSE., BotPerm, dummyint) 2807 CALL MakePermUsingMask( Model, NullSolver, Mesh, FrontMaskName, & 2808 .FALSE., FrontPerm, dummyint) 2809 2810 !Find frontal BC from logical 2811 DO i=1,Model % NumberOfBCs 2812 ThisBC = ListGetLogical(Model % BCs(i) % Values,FrontMaskName,Found) 2813 IF((.NOT. Found) .OR. (.NOT. ThisBC)) CYCLE 2814 FrontBCtag = Model % BCs(i) % Tag 2815 EXIT 2816 END DO 2817 2818 !Find basal BC from logical 2819 DO i=1,Model % NumberOfBCs 2820 ThisBC = ListGetLogical(Model % BCs(i) % Values,BotMaskName,Found) 2821 IF((.NOT. Found) .OR. (.NOT. ThisBC)) CYCLE 2822 BasalBCtag = Model % BCs(i) % Tag 2823 EXIT 2824 END DO 2825 2826 CountSharedExFrontal = 0 2827 CountSharedNewBasal = 0 2828 SharedExGlobal = 0 2829 SharedNewGlobal = 0 2830 2831 !--------------------------------------------------- 2832 ! Find elements for conversion, and set node switches 2833 !--------------------------------------------------- 2834 DO i=Mesh % NumberOfBulkElements + 1, & 2835 Mesh % NumberOfBulkElements + Mesh % NumberOfBoundaryElements 2836 2837 Element => Mesh % Elements(i) 2838 IF(Element % BoundaryInfo % Constraint /= FrontBCtag) CYCLE !not on front 2839 IF(Element % TYPE % ElementCode == 101) CYCLE 2840 2841 NodeIndexes => Element % NodeIndexes 2842 2843 IF(.NOT. (ANY(BotPerm(NodeIndexes) > 0) )) CYCLE !not connected to bed 2844 2845 n = Element % TYPE % NumberOfNodes 2846 2847 ALLOCATE(Nodes % x(n), Nodes % y(n), Nodes % z(n)) 2848 2849 Nodes % x = Mesh % Nodes % x(NodeIndexes) 2850 Nodes % y = Mesh % Nodes % y(NodeIndexes) 2851 Nodes % z = Mesh % Nodes % z(NodeIndexes) 2852 2853 Normal = NormalVector(Element, Nodes) 2854 2855 !compare element normal to threshold 2856 IF(Normal(3) < ZThresh) THEN 2857 FoundSome = .TRUE. 2858 2859 !Nodes currently on bed become 'ex frontal nodes' 2860 !Nodes not currently on bed become 'new basal nodes' 2861 DO j=1,SIZE(NodeIndexes) 2862 2863 IF(BotPerm(NodeIndexes(j)) > 0) THEN 2864 IF(.NOT. ExFrontalNode(NodeIndexes(j))) THEN !maybe already got in another elem 2865 2866 ExFrontalNode(NodeIndexes(j)) = .TRUE. 2867 2868 !If node is in another partition, need to pass this info 2869 IF(SIZE(Mesh % ParallelInfo % NeighbourList(NodeIndexes(j)) % Neighbours)>1) THEN 2870 CountSharedExFrontal = CountSharedExFrontal + 1 2871 IF(CountSharedExFrontal > 2) CALL Fatal(FuncName, & 2872 "Found more than 2 ExFrontalNodes on partition boundary...") 2873 2874 SharedExGlobal(CountSharedExFrontal) = Mesh % ParallelInfo % GlobalDofs(NodeIndexes(j)) 2875 END IF 2876 END IF 2877 ELSE 2878 IF(.NOT. NewBasalNode(NodeIndexes(j))) THEN !maybe already got in another elem 2879 2880 NewBasalNode(NodeIndexes(j)) = .TRUE. 2881 2882 !If node is in another partition, need to pass this info 2883 IF(SIZE(Mesh % ParallelInfo % NeighbourList(NodeIndexes(j)) % Neighbours)>1) THEN 2884 CountSharedNewBasal = CountSharedNewBasal + 1 2885 IF(CountSharedNewBasal > 2) CALL Fatal(FuncName, & 2886 "Found more than 2 NewBasalNodes on partition boundary...") 2887 2888 SharedNewGlobal(CountSharedNewBasal) = & 2889 Mesh % ParallelInfo % GlobalDofs(NodeIndexes(j)) 2890 END IF 2891 END IF 2892 2893 2894 END IF 2895 2896 END DO 2897 2898 NewBasalElem(i) = .TRUE. 2899 IF(Debug) PRINT *, ParEnv % MyPE, 'Debug, converting element: ',i,& 2900 ' with nodes: ', NodeIndexes 2901 END IF 2902 2903 DEALLOCATE(Nodes % x, Nodes % y, Nodes % z) 2904 END DO 2905 2906 !Distribute information about shared frontal nodes 2907 !which are no longer on the front. 2908 !NOTE: we may also need to pass NewBasalNodes... 2909 IF(Debug) PRINT *,ParEnv % MyPE, ' Debug, shared ex frontal nodes: ',SharedExGlobal 2910 IF(Debug) PRINT *,ParEnv % MyPE, ' Debug, shared new basal nodes: ',SharedNewGlobal 2911 2912 ALLOCATE(AllSharedExGlobal(2*ParEnv % PEs),& 2913 AllSharedNewGlobal(2*ParEnv % PEs)) 2914 2915 CALL MPI_ALLGATHER(SharedExGlobal,2,MPI_INTEGER,& 2916 AllSharedExGlobal,2,MPI_INTEGER, ELMER_COMM_WORLD, ierr) 2917 CALL MPI_ALLGATHER(SharedNewGlobal,2,MPI_INTEGER,& 2918 AllSharedNewGlobal,2,MPI_INTEGER, ELMER_COMM_WORLD, ierr) 2919 2920 DO i=1,Mesh % NumberOfNodes 2921 IF(FrontPerm(i) <= 0) CYCLE 2922 IF(ANY(AllSharedExGlobal == Mesh % ParallelInfo % GlobalDOFs(i))) THEN 2923 ExFrontalNode(i) = .TRUE. 2924 FoundSome = .TRUE. 2925 IF(Debug) PRINT *, ParEnv % MyPE, ' Debug, received shared exfrontalnode: ',i 2926 END IF 2927 IF(ANY(AllSharedNewGlobal == Mesh % ParallelInfo % GlobalDOFs(i))) THEN 2928 NewBasalNode(i) = .TRUE. 2929 FoundSome = .TRUE. 2930 IF(Debug) PRINT *, ParEnv % MyPE, ' Debug, received shared newbasalnode: ',i 2931 END IF 2932 END DO 2933 2934 !------------------------------------------------------------------------------ 2935 ! Cycle front elements, looking for those to convert 404 -> 303 for front interp 2936 ! And, also, a rare case where one element is sandwiched between shared ExFrontalNodes 2937 ! In this case, cycle 2938 !------------------------------------------------------------------------------ 2939 DO j=1,2 2940 2941 count303 = 0 2942 DO i=Mesh % NumberOfBulkElements + 1, & 2943 Mesh % NumberOfBulkElements + Mesh % NumberOfBoundaryElements 2944 2945 Element => Mesh % Elements(i) 2946 IF(Element % BoundaryInfo % Constraint /= FrontBCtag) CYCLE !not on front 2947 IF(Element % TYPE % ElementCode == 101) CYCLE 2948 IF(NewBasalElem(i)) CYCLE !element disappears from front entirely 2949 2950 NodeIndexes => Element % NodeIndexes 2951 2952 IF(.NOT. (ANY(BotPerm(NodeIndexes) > 0) )) CYCLE 2953 IF(.NOT. ANY(ExFrontalNode(NodeIndexes))) CYCLE !Not affected 2954 2955 IF(j==2 .AND. Debug) PRINT *, ParEnv % MyPE, ' Debug, switching element: ',& 2956 i,' with nodeindexes ', NodeIndexes 2957 2958 IF(COUNT(ExFrontalNode(NodeIndexes)) /= 1) CYCLE 2959 2960 !iff only change one row of elements at at time, we only get here 2961 !through elements to the side which become 303 2962 count303 = count303 + 1 2963 2964 !First time we just count and allocate... 2965 IF(j==2) THEN 2966 DO k=1,2 2967 New303Elements(count303,k) % TYPE => GetElementType( 303, .FALSE. ) 2968 New303Elements(count303,k) % NDOFs = 3 2969 New303Elements(count303,k) % ElementIndex = i 2970 New303Elements(count303,k) % BodyID = Element % BodyID 2971 2972 ALLOCATE(New303Elements(count303,k) % NodeIndexes(3)) 2973 END DO 2974 2975 !The temporary frontal element 2976 New303Elements(count303,1) % NodeIndexes = & 2977 PACK(NodeIndexes, (.NOT. ExFrontalNode(NodeIndexes))) 2978 2979 !The temporary basal element 2980 New303Elements(count303,2) % NodeIndexes = & 2981 PACK(NodeIndexes, ( (BotPerm(NodeIndexes)>0) .OR. NewBasalNode(NodeIndexes) ) ) 2982 2983 DO k=1,2 2984 ALLOCATE(New303Elements(count303,k) % BoundaryInfo) 2985 New303Elements(count303,k) % BoundaryInfo % Left => Element % BoundaryInfo % Left 2986 New303Elements(count303,k) % BoundaryInfo % Right => Element % BoundaryInfo % Right 2987 2988 IF(k==1) THEN 2989 n = FrontBCtag 2990 ELSE 2991 n = BasalBCtag 2992 END IF 2993 2994 New303Elements(count303,k) % BoundaryInfo % Constraint = n 2995 END DO 2996 2997 IF(Debug) PRINT *, ParEnv % MyPE, ' debug, new frontal element ',i,' has nodes: ', & 2998 New303Elements(count303,1) % NodeIndexes 2999 3000 IF(Debug) PRINT *, ParEnv % MyPE, ' debug, new basal element ',i,' has nodes: ', & 3001 New303Elements(count303,2) % NodeIndexes 3002 END IF 3003 END DO 3004 3005 IF(j==1) THEN 3006 ALLOCATE(New303Elements(count303,2)) 3007 END IF 3008 3009 END DO 3010 3011 !------------------------------------------------------- 3012 ! Now modify mesh % elements accordingly 3013 !------------------------------------------------------- 3014 IF(FoundSome) THEN 3015 3016 OldElemCount = Mesh % NumberOfBulkElements + & 3017 Mesh % NumberOfBoundaryElements 3018 NewElemCount = OldElemCount + count303 3019 3020 ALLOCATE(WorkElements(NewElemCount)) 3021 WorkElements(1:OldElemCount) = Mesh % Elements(1:OldElemCount) 3022 3023 DO i=1,count303 3024 n = New303Elements(i,1) % ElementIndex 3025 3026 Element => WorkElements(n) 3027 3028 CALL FreeElementStuff(Element) 3029 3030 Element = New303Elements(i,1) 3031 Element => WorkElements(OldElemCount + i) 3032 3033 Element = New303Elements(i,2) 3034 Element % ElementIndex = OldElemCount + i 3035 END DO 3036 3037 ! Change constraint on NewBasalElem 3038 DO i=LBOUND(NewBasalElem,1),UBOUND(NewBasalElem,1) 3039 IF(.NOT. NewBasalElem(i)) CYCLE 3040 WorkElements(i) % BoundaryInfo % Constraint = BasalBCtag 3041 END DO 3042 3043 DEALLOCATE(Mesh % Elements) 3044 Mesh % NumberOfBoundaryElements = Mesh % NumberOfBoundaryElements + count303 3045 Mesh % Elements => WorkElements 3046 END IF 3047 3048 CALL SParIterAllReduceOR(FoundSome) 3049 3050 NULLIFY(WorkElements) 3051 3052 !TODO: Free New303Elements 3053 DEALLOCATE(AllSharedExGlobal, AllSharedNewGlobal, & 3054 NewBasalElem, FrontPerm, BotPerm, ExFrontalNode) 3055 3056 END SUBROUTINE ConvertFrontalToBasal 3057 3058 SUBROUTINE FreeElementStuff(Element) 3059 TYPE(Element_t), POINTER :: Element 3060 IF(ASSOCIATED(Element % NodeIndexes)) DEALLOCATE(Element % NodeIndexes) 3061 IF(ASSOCIATED(Element % EdgeIndexes)) DEALLOCATE(Element % EdgeIndexes) 3062 IF(ASSOCIATED(Element % FaceIndexes)) DEALLOCATE(Element % FaceIndexes) 3063 IF(ASSOCIATED(Element % BubbleIndexes)) DEALLOCATE(Element % BubbleIndexes) 3064 IF(ASSOCIATED(Element % DGIndexes)) DEALLOCATE(Element % DGIndexes) 3065 IF(ASSOCIATED(Element % PDefs)) DEALLOCATE(Element % PDefs) 3066 END SUBROUTINE FreeElementStuff 3067 3068 3069 !Turns off (or back on) a specified solver, and adds a string "Save Exec When" 3070 ! to solver % values to allow it to be switched back on to the correct setting. 3071 SUBROUTINE SwitchSolverExec(Solver, Off) 3072 3073 IMPLICIT NONE 3074 3075 TYPE(Solver_t) :: Solver 3076 LOGICAL :: Off 3077 !----------------------------------------- 3078 CHARACTER(MAX_NAME_LEN) :: SaveExecWhen 3079 LOGICAL :: Found 3080 3081 SaveExecWhen = ListGetString(Solver % Values, "Save Exec When", Found) 3082 IF(.NOT. Found) THEN 3083 SaveExecWhen = ListGetString(Solver % Values, 'Exec Solver', Found) 3084 IF(.NOT. Found) SaveExecWhen = 'always' 3085 CALL ListAddString(Solver % Values, 'Save Exec When', SaveExecWhen) 3086 END IF 3087 3088 IF(Off) THEN 3089 3090 !Turning the solver off 3091 Solver % SolverExecWhen = SOLVER_EXEC_NEVER 3092 CALL ListAddString(Solver % Values, 'Exec Solver', 'Never') 3093 3094 ELSE 3095 3096 CALL ListAddString(Solver % Values, 'Exec Solver', SaveExecWhen) 3097 3098 SELECT CASE( SaveExecWhen ) 3099 CASE( 'never' ) 3100 Solver % SolverExecWhen = SOLVER_EXEC_NEVER 3101 CASE( 'always' ) 3102 Solver % SolverExecWhen = SOLVER_EXEC_ALWAYS 3103 CASE( 'after simulation', 'after all' ) 3104 Solver % SolverExecWhen = SOLVER_EXEC_AFTER_ALL 3105 CASE( 'before simulation', 'before all' ) 3106 Solver % SolverExecWhen = SOLVER_EXEC_AHEAD_ALL 3107 CASE( 'before timestep' ) 3108 Solver % SolverExecWhen = SOLVER_EXEC_AHEAD_TIME 3109 CASE( 'after timestep' ) 3110 Solver % SolverExecWhen = SOLVER_EXEC_AFTER_TIME 3111 CASE( 'before saving' ) 3112 Solver % SolverExecWhen = SOLVER_EXEC_AHEAD_SAVE 3113 CASE( 'after saving' ) 3114 Solver % SolverExecWhen = SOLVER_EXEC_AFTER_SAVE 3115 CASE DEFAULT 3116 CALL Fatal("SwitchSolverExec","Programming error here...") 3117 END SELECT 3118 3119 END IF 3120 3121 END SUBROUTINE SwitchSolverExec 3122 3123 SUBROUTINE PlanePointIntersection ( pp, pnorm, p1, p2, p_intersect, found_intersection ) 3124 !Get the intersection point between a line and plane in 3D 3125 ! Plane defined by point "pp" and norm "pnorm", line defined by points "p1" and "p2" 3126 ! Intersection returned in p_intersect 3127 !found_intersection = .FALSE. if they happen to be parallel 3128 3129 REAL(KIND=dp) :: pp(3), pnorm(3), p1(3), p2(3), p_intersect(3) 3130 LOGICAL :: found_intersection 3131 !---------------------------- 3132 REAL(KIND=dp) :: pl(3), dist 3133 3134 pl = p2 - p1 3135 3136 IF(ABS(DOT_PRODUCT(pl,pnorm)) < EPSILON(1.0_dp)) THEN 3137 !Line and plane are parallel... 3138 found_intersection = .FALSE. 3139 RETURN 3140 END IF 3141 3142 dist = DOT_PRODUCT((pp - p1), pnorm) / DOT_PRODUCT(pl,pnorm) 3143 3144 p_intersect = p1 + dist*pl 3145 found_intersection = .TRUE. 3146 3147 END SUBROUTINE PlanePointIntersection 3148 3149 SUBROUTINE LineSegmentsIntersect ( a1, a2, b1, b2, intersect_point, does_intersect ) 3150 ! Find if two 2D line segments intersect 3151 ! Line segment 'a' runs from point a1 => a2, same for b 3152 3153 IMPLICIT NONE 3154 3155 REAL(KIND=dp) :: a1(2), a2(2), b1(2), b2(2), intersect_point(2) 3156 LOGICAL :: does_intersect 3157 !----------------------- 3158 REAL(KIND=dp) :: r(2), s(2), rxs, bma(2), t, u 3159 3160 3161 does_intersect = .FALSE. 3162 intersect_point = 0.0_dp 3163 3164 r = a2 - a1 3165 s = b2 - b1 3166 3167 rxs = VecCross2D(r,s) 3168 3169 IF(rxs == 0.0_dp) RETURN 3170 3171 bma = b1 - a1 3172 3173 t = VecCross2D(bma,s) / rxs 3174 u = VecCross2D(bma,r) / rxs 3175 3176 IF(t < 0.0_dp .OR. t > 1.0_dp .OR. u < 0.0_dp .OR. u > 1.0_dp) RETURN 3177 3178 intersect_point = a1 + (t * r) 3179 does_intersect = .TRUE. 3180 3181 END SUBROUTINE LineSegmentsIntersect 3182 3183 SUBROUTINE LinesIntersect ( a1, a2, b1, b2, intersect_point, does_intersect ) 3184 ! Find where two 2D lines intersect 3185 ! Line 'a' explicitly defined by points a1, a2 which lie on line, same for b 3186 ! based on LineSegmentsIntersect above 3187 3188 IMPLICIT NONE 3189 3190 REAL(KIND=dp) :: a1(2), a2(2), b1(2), b2(2), intersect_point(2) 3191 LOGICAL :: does_intersect 3192 !----------------------- 3193 REAL(KIND=dp) :: r(2), s(2), rxs, bma(2), t, u 3194 3195 3196 does_intersect = .TRUE. 3197 3198 intersect_point = 0.0_dp 3199 3200 r = a2 - a1 3201 s = b2 - b1 3202 3203 rxs = VecCross2D(r,s) 3204 3205 IF(rxs == 0.0_dp) THEN 3206 does_intersect = .FALSE. 3207 RETURN 3208 ENDIF 3209 3210 bma = b1 - a1 3211 3212 t = VecCross2D(bma,s) / rxs 3213 u = VecCross2D(bma,r) / rxs 3214 3215 intersect_point = a1 + (t * r) 3216 3217 END SUBROUTINE LinesIntersect 3218 3219 FUNCTION VecCross2D(a, b) RESULT (c) 3220 REAL(KIND=dp) :: a(2), b(2), c 3221 3222 c = a(1)*b(2) - a(2)*b(1) 3223 3224 END FUNCTION VecCross2D 3225 3226 !This subroutine should identify discrete calving events for the 3227 !purposes of local remeshing. For now it returns 1 3228 SUBROUTINE CountCalvingEvents(Model, Mesh,CCount) 3229 TYPE(Model_t) :: Model 3230 TYPE(Mesh_t),POINTER :: Mesh 3231 INTEGER :: CCount 3232 3233 Ccount = 1 3234 END SUBROUTINE CountCalvingEvents 3235 3236 ! shortest distance of c to segment ab, a b and c are in 2D 3237 FUNCTION PointLineSegmDist2D(a, b, c) RESULT (pdis) 3238 REAL(KIND=dp) :: a(2), b(2), c(2), n(2), v(2), dd, t, pdis 3239 n=b-a ! Vector ab 3240 dd = (n(1)**2.+n(2)**2.) ! Length of ab squared 3241 dd = DOT_PRODUCT(n,n) ! alternative writing 3242 t = DOT_PRODUCT(c-a,b-a)/dd 3243 dd = MAXVAL( (/0.0_dp, MINVAL( (/1.0_dp,t/) ) /) ) 3244 v = c - a - dd * n 3245 pdis=sqrt(v(1)**2.+v(2)**2.) 3246 END FUNCTION PointLineSegmDist2D 3247 3248 ! Takes two meshes which are assumed to represent the same domain 3249 ! and interpolates variables between them. Uses full dimension 3250 ! interpolation (InterpolateMeshToMesh) for all nodes, then picks 3251 ! up missing boundary nodes using reduced dim 3252 ! (InterpolateVarToVarReduced) 3253 SUBROUTINE SwitchMesh(Model, Solver, OldMesh, NewMesh) 3254 3255 IMPLICIT NONE 3256 3257 TYPE(Model_t) :: Model 3258 TYPE(Solver_t) :: Solver 3259 TYPE(Mesh_t), POINTER :: OldMesh, NewMesh 3260 !------------------------------------------------- 3261 TYPE(Solver_t), POINTER :: WorkSolver 3262 TYPE(Variable_t), POINTER :: Var=>NULL(), NewVar=>NULL(), WorkVar=>NULL() 3263 TYPE(Valuelist_t), POINTER :: Params 3264 TYPE(Matrix_t), POINTER :: WorkMatrix=>NULL() 3265 LOGICAL :: Found, Global, GlobalBubbles, Debug, DoPrevValues, & 3266 NoMatrix, DoOptimizeBandwidth, PrimaryVar, HasValuesInPartition, & 3267 PrimarySolver 3268 LOGICAL, POINTER :: UnfoundNodes(:)=>NULL() 3269 INTEGER :: i,j,k,DOFs, nrows,n 3270 INTEGER, POINTER :: WorkPerm(:)=>NULL(), SolversToIgnore(:)=>NULL() 3271 REAL(KIND=dp), POINTER :: WorkReal(:)=>NULL(), WorkReal2(:)=>NULL(), PArray(:,:) => NULL() 3272 REAL(KIND=dp) :: FrontOrientation(3), RotationMatrix(3,3), UnRotationMatrix(3,3), & 3273 globaleps, localeps 3274 CHARACTER(LEN=MAX_NAME_LEN) :: SolverName, WorkName 3275 3276 INTERFACE 3277 SUBROUTINE InterpolateMeshToMesh( OldMesh, NewMesh, OldVariables, & 3278 NewVariables, UseQuadrantTree, Projector, MaskName, UnfoundNodes ) 3279 !------------------------------------------------------------------------------ 3280 USE Lists 3281 USE SParIterComm 3282 USE Interpolation 3283 USE CoordinateSystems 3284 !------------------------------------------------------------------------------- 3285 TYPE(Mesh_t), TARGET :: OldMesh, NewMesh 3286 TYPE(Variable_t), POINTER, OPTIONAL :: OldVariables, NewVariables 3287 LOGICAL, OPTIONAL :: UseQuadrantTree 3288 LOGICAL, POINTER, OPTIONAL :: UnfoundNodes(:) 3289 TYPE(Projector_t), POINTER, OPTIONAL :: Projector 3290 CHARACTER(LEN=*),OPTIONAL :: MaskName 3291 END SUBROUTINE InterpolateMeshToMesh 3292 END INTERFACE 3293 3294 SolverName = "SwitchMesh" 3295 Debug = .FALSE. 3296 Params => Solver % Values 3297 CALL Info( 'Remesher', ' ',Level=4 ) 3298 CALL Info( 'Remesher', '-------------------------------------',Level=4 ) 3299 CALL Info( 'Remesher', ' Switching from old to new mesh...',Level=4 ) 3300 CALL Info( 'Remesher', '-------------------------------------',Level=4 ) 3301 CALL Info( 'Remesher', ' ',Level=4 ) 3302 3303 IF(ASSOCIATED(NewMesh % Variables)) CALL Fatal(SolverName,& 3304 "New mesh already has variables associated!") 3305 3306 !interpolation epsilons 3307 globaleps = 1.0E-2_dp 3308 localeps = 1.0E-2_dp 3309 3310 !---------------------------------------------- 3311 ! Get the orientation of the calving front 3312 ! & compute rotation matrix 3313 !---------------------------------------------- 3314 PArray => ListGetConstRealArray( Model % Constants,'Front Orientation', & 3315 Found, UnfoundFatal=.TRUE.) 3316 DO i=1,3 3317 FrontOrientation(i) = PArray(i,1) 3318 END DO 3319 RotationMatrix = ComputeRotationMatrix(FrontOrientation) 3320 UnRotationMatrix = TRANSPOSE(RotationMatrix) 3321 3322 !---------------------------------------------- 3323 ! Action 3324 !---------------------------------------------- 3325 3326 CALL CopyIntrinsicVars(OldMesh, NewMesh) 3327 3328 !---------------------------------------------- 3329 ! Add Variables to NewMesh 3330 !---------------------------------------------- 3331 3332 Var => OldMesh % Variables 3333 DO WHILE( ASSOCIATED(Var) ) 3334 3335 DoPrevValues = ASSOCIATED(Var % PrevValues) 3336 WorkSolver => Var % Solver 3337 HasValuesInPartition = .TRUE. 3338 3339 !Do nothing if it already exists 3340 !e.g. it's a DOF component added previously 3341 NewVar => VariableGet( NewMesh % Variables, Var % Name, ThisOnly = .TRUE.) 3342 IF(ASSOCIATED(NewVar)) THEN 3343 NULLIFY(NewVar) 3344 Var => Var % Next 3345 CYCLE 3346 END IF 3347 3348 DOFs = Var % DOFs 3349 Global = (SIZE(Var % Values) .EQ. DOFs) 3350 3351 !Allocate storage for values and perm 3352 IF(Global) THEN 3353 ALLOCATE(WorkReal(DOFs)) 3354 WorkReal = Var % Values 3355 3356 CALL VariableAdd( NewMesh % Variables, NewMesh, & 3357 Var % Solver, TRIM(Var % Name), & 3358 Var % DOFs, WorkReal) 3359 3360 ELSE !Regular field variable 3361 ALLOCATE(WorkPerm(NewMesh % NumberOfNodes)) 3362 3363 IF(.NOT. ASSOCIATED(WorkSolver)) THEN 3364 WRITE(Message, '(a,a,a)') "Variable ",Var % Name," has no solver, unexpected." 3365 CALL Fatal(SolverName, Message) 3366 END IF 3367 3368 PrimaryVar = ASSOCIATED(WorkSolver % Variable, Var) 3369 3370 IF(PrimaryVar) THEN !Take care of the matrix 3371 NoMatrix = ListGetLogical( WorkSolver % Values, 'No matrix',Found) 3372 !Issue here, this will recreate matrix for every variable associated w/ solver. 3373 3374 IF(.NOT. NoMatrix) THEN 3375 IF(ParEnv % MyPE == 0) PRINT *, 'Computing matrix for variable: ',TRIM(Var % Name) 3376 3377 DoOptimizeBandwidth = ListGetLogical( WorkSolver % Values, & 3378 'Optimize Bandwidth', Found ) 3379 IF ( .NOT. Found ) DoOptimizeBandwidth = .TRUE. 3380 3381 GlobalBubbles = ListGetLogical( WorkSolver % Values, & 3382 'Bubbles in Global System', Found ) 3383 IF ( .NOT. Found ) GlobalBubbles = .TRUE. 3384 3385 WorkMatrix => CreateMatrix(Model, WorkSolver, & 3386 NewMesh, WorkPerm, DOFs, MATRIX_CRS, DoOptimizeBandwidth, & 3387 ListGetString( WorkSolver % Values, 'Equation' ), & 3388 GlobalBubbles = GlobalBubbles ) 3389 3390 IF(ASSOCIATED(WorkMatrix)) THEN 3391 WorkMatrix % Comm = ELMER_COMM_WORLD 3392 3393 WorkMatrix % Symmetric = ListGetLogical( WorkSolver % Values, & 3394 'Linear System Symmetric', Found ) 3395 3396 WorkMatrix % Lumped = ListGetLogical( WorkSolver % Values, & 3397 'Lumped Mass Matrix', Found ) 3398 3399 CALL AllocateVector( WorkMatrix % RHS, WorkMatrix % NumberOfRows ) 3400 WorkMatrix % RHS = 0.0_dp 3401 WorkMatrix % RHS_im => NULL() 3402 3403 ALLOCATE(WorkMatrix % Force(WorkMatrix % NumberOfRows, WorkSolver % TimeOrder+1)) 3404 WorkMatrix % Force = 0.0_dp 3405 ELSE 3406 !No nodes in this partition now 3407 NoMatrix = .TRUE. 3408 END IF 3409 END IF 3410 3411 IF ( ASSOCIATED(Var % EigenValues) ) THEN 3412 n = SIZE(Var % EigenValues) 3413 3414 IF ( n > 0 ) THEN 3415 WorkSolver % NOFEigenValues = n 3416 CALL AllocateVector( NewVar % EigenValues,n ) 3417 CALL AllocateArray( NewVar % EigenVectors, n, & 3418 SIZE(NewVar % Values) ) 3419 3420 NewVar % EigenValues = 0.0d0 3421 NewVar % EigenVectors = 0.0d0 3422 IF(.NOT.NoMatrix) THEN 3423 CALL AllocateVector( WorkMatrix % MassValues, SIZE(WorkMatrix % Values) ) 3424 WorkMatrix % MassValues = 0.0d0 3425 END IF 3426 END IF 3427 END IF 3428 3429 !Check for duplicate solvers with same var 3430 !Nullify/deallocate and repoint the matrix 3431 !Note: previously this DO loop was after the FreeMatrix 3432 !and pointing below, but this caused double free errors 3433 DO j=1,Model % NumberOfSolvers 3434 IF(ASSOCIATED(WorkSolver, Model % Solvers(j))) CYCLE 3435 IF(.NOT. ASSOCIATED(Model % Solvers(j) % Variable)) CYCLE 3436 IF( TRIM(Model % Solvers(j) % Variable % Name) /= TRIM(Var % Name)) CYCLE 3437 3438 !If the other solver's matrix is the same as WorkSolver matrix, we just 3439 !nullify, otherwise we deallocate. After the first timestep, solvers 3440 !with the same variable will have the same matrix 3441 IF(ASSOCIATED(Model % Solvers(j) % Matrix, WorkSolver % Matrix)) THEN 3442 Model % Solvers(j) % Matrix => NULL() 3443 ELSE 3444 CALL FreeMatrix(Model % Solvers(j) % Matrix) 3445 END IF 3446 !Point this other solver % matrix to the matrix we just created 3447 Model % Solvers(j) % Matrix => WorkMatrix 3448 END DO 3449 3450 !Deallocate the old matrix & repoint 3451 IF(ASSOCIATED(WorkSolver % Matrix)) CALL FreeMatrix(WorkSolver % Matrix) 3452 WorkSolver % Matrix => WorkMatrix 3453 3454 NULLIFY(WorkMatrix) 3455 3456 !NOTE: We don't switch Solver % Variable here, because 3457 !Var % Solver % Var doesn't necessarily point to self 3458 !if solver has more than one variable. We do this below. 3459 ELSE 3460 k = InitialPermutation(WorkPerm, Model, WorkSolver, & 3461 NewMesh, ListGetString(WorkSolver % Values,'Equation')) 3462 END IF !Primary var 3463 3464 HasValuesInPartition = COUNT(WorkPerm>0) > 0 3465 IF(HasValuesInPartition) THEN 3466 ALLOCATE(WorkReal(COUNT(WorkPerm>0)*DOFs)) 3467 ELSE 3468 !this is silly but it matches AddEquationBasics 3469 ALLOCATE(WorkReal(NewMesh % NumberOfNodes * DOFs)) 3470 END IF 3471 3472 WorkReal = 0.0_dp 3473 CALL VariableAdd( NewMesh % Variables, NewMesh, & 3474 Var % Solver, TRIM(Var % Name), & 3475 Var % DOFs, WorkReal, WorkPerm, & 3476 Var % Output, Var % Secondary, Var % TYPE ) 3477 3478 END IF !Not global 3479 3480 NewVar => VariableGet( NewMesh % Variables, Var % Name, ThisOnly = .TRUE. ) 3481 IF(.NOT.ASSOCIATED(NewVar)) CALL Fatal(SolverName,& 3482 "Problem creating variable on new mesh.") 3483 3484 IF(DoPrevValues) THEN 3485 ALLOCATE(NewVar % PrevValues( SIZE(NewVar % Values), SIZE(Var % PrevValues,2) )) 3486 END IF 3487 3488 !Add the components of variables with more than one DOF 3489 !NOTE, this implementation assumes the vector variable 3490 !comes before the scalar components in the list. 3491 !e.g., we add Mesh Update and so here we add MU 1,2,3 3492 !SO: next time round, new variable (MU 1) already exists 3493 !and so it's CYCLE'd 3494 IF((DOFs > 1) .AND. (.NOT.Global)) THEN 3495 nrows = SIZE(WorkReal) 3496 DO i=1,DOFs 3497 3498 WorkReal2 => WorkReal( i:nrows-DOFs+i:DOFs ) 3499 WorkName = ComponentName(TRIM(Var % Name),i) 3500 CALL VariableAdd( NewMesh % Variables, NewMesh, & 3501 Var % Solver, WorkName, & 3502 1, WorkReal2, WorkPerm, & 3503 Var % Output, Var % Secondary, Var % TYPE ) 3504 3505 IF(DoPrevValues) THEN 3506 WorkVar => VariableGet( NewMesh % Variables, WorkName, .TRUE. ) 3507 IF(.NOT. ASSOCIATED(WorkVar)) CALL Fatal(SolverName, & 3508 "Error allocating Remesh Update PrevValues.") 3509 3510 NULLIFY(WorkVar % PrevValues) 3511 WorkVar % PrevValues => NewVar % PrevValues(i:nrows-DOFs+i:DOFs,:) 3512 END IF 3513 3514 NULLIFY(WorkReal2) 3515 END DO 3516 END IF 3517 3518 NULLIFY(WorkReal, WorkPerm) 3519 Var => Var % Next 3520 END DO 3521 3522 !Go back through and set non-primary variables to have same % perm as the primary var. 3523 !Bit of a hack - would be nice to somehow do this in one loop... 3524 !Set perms equal if: variable has solver, solver has variable, both variables have perm 3525 Var => NewMesh % Variables 3526 DO WHILE (ASSOCIATED(Var)) 3527 3528 WorkSolver => Var % Solver 3529 IF(ASSOCIATED(WorkSolver)) THEN 3530 IF(ASSOCIATED(WorkSolver % Variable % Perm)) THEN 3531 WorkVar => VariableGet(NewMesh % Variables, & 3532 WorkSolver % Variable % Name, .TRUE., UnfoundFatal=.TRUE.) 3533 PrimaryVar = ASSOCIATED(WorkSolver % Variable, Var) 3534 IF(ASSOCIATED(WorkVar) .AND. .NOT. PrimaryVar) THEN 3535 IF(ASSOCIATED(WorkVar % Perm) .AND. ASSOCIATED(Var % Perm)) THEN 3536 Var % Perm = WorkVar % Perm 3537 END IF 3538 END IF 3539 END IF 3540 END IF 3541 3542 Var => Var % Next 3543 END DO 3544 3545 !set partitions to active, so variable can be -global -nooutput 3546 CALL ParallelActive(.TRUE.) 3547 !MPI_BSend buffer issue in this call to InterpolateMeshToMesh 3548 !Free quadrant tree to ensure its rebuilt in InterpolateMeshToMesh (bug fix) 3549 CALL FreeQuadrantTree(OldMesh % RootQuadrant) 3550 CALL InterpolateMeshToMesh( OldMesh, NewMesh, OldMesh % Variables, UnfoundNodes=UnfoundNodes) 3551 IF(ANY(UnfoundNodes)) THEN 3552 PRINT *, ParEnv % MyPE, ' missing ', COUNT(UnfoundNodes),' out of ',SIZE(UnfoundNodes),& 3553 ' nodes in SwitchMesh.' 3554 END IF 3555 3556 !--------------------------------------------------------- 3557 ! For top, bottom and calving front BC, do reduced dim 3558 ! interpolation to avoid epsilon problems 3559 !--------------------------------------------------------- 3560 3561 CALL InterpMaskedBCReduced(Model, Solver, OldMesh, NewMesh, OldMesh % Variables, & 3562 "Top Surface Mask",globaleps=globaleps,localeps=localeps) 3563 CALL InterpMaskedBCReduced(Model, Solver, OldMesh, NewMesh, OldMesh % Variables, & 3564 "Bottom Surface Mask",globaleps=globaleps,localeps=localeps) 3565 3566 CALL RotateMesh(OldMesh, RotationMatrix) 3567 CALL RotateMesh(NewMesh, RotationMatrix) 3568 3569 !CHANGE - need to delete UnfoundNOtes from this statement, or front 3570 !variables not copied across. If you get some odd interpolation artefact, 3571 !suspect this 3572 CALL InterpMaskedBCReduced(Model, Solver, OldMesh, NewMesh, OldMesh % Variables, & 3573 "Calving Front Mask",globaleps=globaleps,localeps=localeps) 3574 3575 !NOTE: InterpMaskedBCReduced on the calving front will most likely fail to 3576 ! find a few points, due to vertical adjustment to account for GroundedSolver. 3577 ! Briefly, the 'DoGL' sections of CalvingRemesh adjust the Z coordinate of 3578 ! basal nodes which are grounded, to ensure they match the bed dataset. 3579 ! Thus, it's not impossible for points on the new mesh to sit slightly outside 3580 ! the old. 3581 ! However, these points should sit behind or on the old calving front, so 3582 ! InterpMaskedBC... on the bed should get them. Thus the only thing that may 3583 ! be missed would be variables defined solely on the front. Currently, none 3584 ! of these are important for the next timestep, so this should be fine. 3585 3586 CALL RotateMesh(NewMesh, UnrotationMatrix) 3587 CALL RotateMesh(OldMesh, UnrotationMatrix) 3588 3589 !----------------------------------------------- 3590 ! Point solvers at the correct mesh and variable 3591 !----------------------------------------------- 3592 3593 !CHANGE 3594 !Needs to be told to ignore certain solvers if using multiple meshes 3595 SolversToIgnore => ListGetIntegerArray(Params, 'Solvers To Ignore') 3596 3597 DO i=1,Model % NumberOfSolvers 3598 WorkSolver => Model % Solvers(i) 3599 3600 !CHANGE - see above 3601 IF (ASSOCIATED(SolversToIgnore)) THEN 3602 IF(ANY(SolversToIgnore(1:SIZE(SolversToIgnore))==i)) CYCLE 3603 END IF 3604 3605 WorkSolver % Mesh => NewMesh !note, assumption here that there's only one active mesh 3606 3607 !hack to get SingleSolver to recompute 3608 !should be taken care of by Mesh % Changed, but 3609 !this is reset by CoupledSolver for some reason 3610 WorkSolver % NumberOfActiveElements = -1 3611 3612 IF(.NOT. ASSOCIATED(WorkSolver % Variable)) CYCLE 3613 IF(WorkSolver % Variable % NameLen == 0) CYCLE !dummy !invalid read 3614 3615 !Check for multiple solvers with same var: 3616 !If one of the duplicate solvers is only executed before the simulation (or never), 3617 !then we don't point the variable at this solver. (e.g. initial groundedmask). 3618 !If both solvers are executed during each timestep, we have a problem. 3619 !If neither are, it doesn't matter, and so the the later occurring solver will have 3620 !the variable pointed at it (arbitrary). 3621 PrimarySolver = .TRUE. 3622 DO j=1,Model % NumberOfSolvers 3623 IF(j==i) CYCLE 3624 IF(.NOT. ASSOCIATED(Model % Solvers(j) % Variable)) CYCLE 3625 IF(TRIM(Model % Solvers(j) % Variable % Name) == WorkSolver % Variable % Name) THEN 3626 3627 IF( (WorkSolver % SolverExecWhen == SOLVER_EXEC_NEVER) .OR. & 3628 (WorkSolver % SolverExecWhen == SOLVER_EXEC_AHEAD_ALL) ) THEN 3629 IF((Model % Solvers(j) % SolverExecWhen == SOLVER_EXEC_NEVER) .OR. & 3630 (Model % Solvers(j) % SolverExecWhen == SOLVER_EXEC_AHEAD_ALL) ) THEN 3631 PrimarySolver = .TRUE. 3632 ELSE 3633 PrimarySolver = .FALSE. 3634 WorkSolver % Matrix => NULL() 3635 EXIT 3636 END IF 3637 ELSE 3638 IF( (Model % Solvers(j) % SolverExecWhen == SOLVER_EXEC_NEVER) .OR. & 3639 (Model % Solvers(j) % SolverExecWhen == SOLVER_EXEC_AHEAD_ALL) ) THEN 3640 PrimarySolver = .TRUE. 3641 EXIT 3642 ELSE 3643 WRITE(Message, '(A,A)') "Unable to determine main solver for variable: ", & 3644 TRIM(WorkSolver % Variable % Name) 3645 CALL Fatal(SolverName, Message) 3646 END IF 3647 END IF 3648 3649 END IF 3650 END DO 3651 3652 WorkVar => VariableGet(NewMesh % Variables, & 3653 WorkSolver % Variable % Name, .TRUE.) !invalid read 3654 3655 IF(ASSOCIATED(WorkVar)) THEN 3656 WorkSolver % Variable => WorkVar 3657 IF(PrimarySolver) WorkVar % Solver => WorkSolver 3658 ELSE 3659 WRITE(Message, '(a,a,a)') "Variable ",WorkSolver % Variable % Name," wasn't & 3660 &correctly switched to the new mesh." !invalid read 3661 PRINT *, i,' debug, solver equation: ', ListGetString(WorkSolver % Values, "Equation") 3662 CALL Fatal(SolverName, Message) 3663 END IF 3664 3665 END DO 3666 3667 3668 NewMesh % Next => OldMesh % Next 3669 Model % Meshes => NewMesh 3670 Model % Mesh => NewMesh 3671 Model % Variables => NewMesh % Variables 3672 3673 !Free old mesh and associated variables 3674 CALL ReleaseMesh(OldMesh) 3675 DEALLOCATE(OldMesh) 3676 DEALLOCATE(UnfoundNodes) 3677 3678 OldMesh => Model % Meshes 3679 3680 END SUBROUTINE SwitchMesh 3681 3682 SUBROUTINE InterpMaskedBCReduced(Model, Solver, OldMesh, NewMesh, Variables, MaskName, & 3683 SeekNodes, globaleps, localeps) 3684 3685 USE InterpVarToVar 3686 3687 IMPLICIT NONE 3688 3689 TYPE(Model_t) :: Model 3690 TYPE(Solver_t) :: Solver 3691 TYPE(Mesh_t), POINTER :: OldMesh, NewMesh 3692 TYPE(Variable_t), POINTER :: Variables 3693 INTEGER, POINTER :: OldMaskPerm(:)=>NULL(), NewMaskPerm(:)=>NULL() 3694 INTEGER, POINTER :: InterpDim(:) 3695 INTEGER :: i,dummyint 3696 REAL(KIND=dp), OPTIONAL :: globaleps,localeps 3697 REAL(KIND=dp) :: geps,leps 3698 LOGICAL :: Debug 3699 LOGICAL, POINTER :: OldMaskLogical(:), NewMaskLogical(:), UnfoundNodes(:)=>NULL() 3700 LOGICAL, POINTER, OPTIONAL :: SeekNodes(:) 3701 CHARACTER(LEN=*) :: MaskName 3702 3703 CALL MakePermUsingMask( Model, Solver, NewMesh, MaskName, & 3704 .FALSE., NewMaskPerm, dummyint) 3705 3706 CALL MakePermUsingMask( Model, Solver, OldMesh, MaskName, & 3707 .FALSE., OldMaskPerm, dummyint) 3708 3709 ALLOCATE(OldMaskLogical(SIZE(OldMaskPerm)),& 3710 NewMaskLogical(SIZE(NewMaskPerm))) 3711 3712 OldMaskLogical = (OldMaskPerm <= 0) 3713 NewMaskLogical = (NewMaskPerm <= 0) 3714 IF(PRESENT(SeekNodes)) NewMaskLogical = & 3715 NewMaskLogical .OR. .NOT. SeekNodes 3716 3717 IF(PRESENT(globaleps)) THEN 3718 geps = globaleps 3719 ELSE 3720 geps = 1.0E-4 3721 END IF 3722 3723 IF(PRESENT(localeps)) THEN 3724 leps = localeps 3725 ELSE 3726 leps = 1.0E-4 3727 END IF 3728 3729 IF(Debug) PRINT *, ParEnv % MyPE,'Debug, on boundary: ',TRIM(MaskName),' seeking ',& 3730 COUNT(.NOT. NewMaskLogical),' of ',SIZE(NewMaskLogical),' nodes.' 3731 3732 ALLOCATE(InterpDim(1)) 3733 InterpDim(1) = 3 3734 3735 CALL ParallelActive(.TRUE.) 3736 CALL InterpolateVarToVarReduced(OldMesh, NewMesh, "remesh update 1", InterpDim, & 3737 UnfoundNodes, OldMaskLogical, NewMaskLogical, Variables=OldMesh % Variables, & 3738 GlobalEps=geps, LocalEps=leps) 3739 3740 IF(ANY(UnfoundNodes)) THEN 3741 !NewMaskLogical changes purpose, now it masks supporting nodes 3742 NewMaskLogical = (NewMaskPerm <= 0) 3743 3744 DO i=1, SIZE(UnfoundNodes) 3745 IF(UnfoundNodes(i)) THEN 3746 PRINT *,ParEnv % MyPE,'Didnt find point: ', i, & 3747 ' x:', NewMesh % Nodes % x(i),& 3748 ' y:', NewMesh % Nodes % y(i),& 3749 ' z:', NewMesh % Nodes % z(i) 3750 3751 CALL InterpolateUnfoundPoint( i, NewMesh, "remesh update 1", InterpDim, & 3752 NodeMask=NewMaskLogical, Variables=NewMesh % Variables ) 3753 END IF 3754 END DO 3755 3756 WRITE(Message, '(i0,a,a,a,i0,a,i0,a)') ParEnv % MyPE,& 3757 ' Failed to find all points on face: ',MaskName, ', ',& 3758 COUNT(UnfoundNodes),' of ',COUNT(.NOT. NewMaskLogical),' missing points.' 3759 CALL Warn("InterpMaskedBCReduced", Message) 3760 END IF 3761 3762 DEALLOCATE(OldMaskLogical, & 3763 NewMaskLogical, NewMaskPerm, & 3764 OldMaskPerm, UnfoundNodes) 3765 3766 END SUBROUTINE InterpMaskedBCReduced 3767 3768 !Function to return the orientation of a calving front 3769 !If specified in SIF, returns this, otherwise computes it 3770 FUNCTION GetFrontOrientation(Model) RESULT (Orientation) 3771 TYPE(Model_t) :: Model 3772 !-------------------------- 3773 INTEGER :: i 3774 REAL(KIND=dp) :: Orientation(3),OrientSaved(3) 3775 REAL(KIND=dp), POINTER :: PArray(:,:) => NULL() 3776 3777 LOGICAL :: FirstTime=.TRUE.,Constant 3778 3779 SAVE :: FirstTime,Constant,PArray,OrientSaved 3780 3781 IF(FirstTime) THEN 3782 FirstTime = .FALSE. 3783 !TODO - this will need to be defined on individual boundary conditions 3784 !if we want to handle multiple calving fronts in same simulation. 3785 PArray => ListGetConstRealArray( Model % Constants,'Front Orientation', & 3786 Constant) 3787 DO i=1,3 3788 OrientSaved(i) = PArray(i,1) 3789 END DO 3790 IF(Constant) THEN 3791 CALL Info("GetFrontOrientation","Using predefined Front Orientation from SIF.", Level=6) 3792 ELSE 3793 CALL Info("GetFrontOrientation","No predefined Front Orientation, computing instead.", Level=6) 3794 END IF 3795 END IF 3796 3797 IF(Constant) THEN 3798 Orientation = OrientSaved 3799 RETURN 3800 ELSE 3801 !Not implemented yet 3802 END IF 3803 3804 END FUNCTION GetFrontOrientation 3805 3806END MODULE CalvingGeometry 3807 3808