1 2!------------------------------------------------------------------------------ 3!> Map results from mesh to mesh. The from-Mesh is stored in an octree from 4!> which it is relatively fast to find the to-nodes. When the node is found 5!> interpolation is performed. Optionally there may be an existing projector 6!> that speeds up the interpolation. 7!------------------------------------------------------------------------------ 8 SUBROUTINE InterpolateMeshToMesh( OldMesh, NewMesh, OldVariables, & 9 NewVariables, UseQuadrantTree, Projector, MaskName, UnfoundNodes ) 10!------------------------------------------------------------------------------ 11 USE Lists 12 USE SParIterComm 13 USE Interpolation 14 USE CoordinateSystems 15 USE MeshUtils, ONLY: ReleaseMesh 16!------------------------------------------------------------------------------- 17 TYPE(Mesh_t), TARGET :: OldMesh, NewMesh 18 TYPE(Variable_t), POINTER, OPTIONAL :: OldVariables, NewVariables 19 LOGICAL, OPTIONAL :: UseQuadrantTree 20 TYPE(Projector_t), POINTER, OPTIONAL :: Projector 21 CHARACTER(LEN=*),OPTIONAL :: MaskName 22 LOGICAL, POINTER, OPTIONAL :: UnfoundNodes(:) 23!------------------------------------------------------------------------------- 24 INTEGER, ALLOCATABLE :: perm(:), vperm(:) 25 INTEGER, POINTER :: nperm(:) 26 LOGICAL, ALLOCATABLE :: FoundNodes(:), FoundNodesPar(:) 27 TYPE(Mesh_t), POINTER :: nMesh 28 TYPE(VAriable_t), POINTER :: Var, nVar 29 INTEGER :: i,j,k,l,nfound,maxrecv,n,ierr,nvars,npart,proc,status(MPI_STATUS_SIZE) 30 REAL(KIND=dp) :: myBB(6), eps2, dn 31 REAL(KIND=dp), POINTER :: store(:) 32 REAL(KIND=dp), ALLOCATABLE, TARGET :: astore(:),vstore(:,:), BB(:,:), & 33 nodes_x(:),nodes_y(:),nodes_z(:), xpart(:), ypart(:), zpart(:) 34 LOGICAL :: al, Stat 35 INTEGER :: PassiveCoordinate 36 37 TYPE ProcRecv_t 38 INTEGER :: n = 0 39 REAL(KIND=dp), ALLOCATABLE :: nodes_x(:),nodes_y(:),nodes_z(:) 40 END TYPE ProcRecv_t 41 TYPE(ProcRecv_t), ALLOCATABLE, TARGET :: ProcRecv(:) 42 43 TYPE ProcSend_t 44 INTEGER :: n = 0 45 INTEGER, ALLOCATABLE :: perm(:) 46 END TYPE ProcSend_t 47 TYPE(ProcSend_t), ALLOCATABLE :: ProcSend(:) 48 49!------------------------------------------------------------------------------- 50 INTERFACE 51 SUBROUTINE InterpolateMeshToMeshQ( OldMesh, NewMesh, OldVariables, NewVariables, & 52 UseQuadrantTree, Projector, MaskName, FoundNodes, NewMaskPerm, KeepUnfoundNodes ) 53 USE Types 54 TYPE(Variable_t), POINTER, OPTIONAL :: OldVariables, NewVariables 55 TYPE(Mesh_t), TARGET :: OldMesh, NewMesh 56 LOGICAL, OPTIONAL :: UseQuadrantTree,FoundNodes(:) 57 CHARACTER(LEN=*),OPTIONAL :: MaskName 58 TYPE(Projector_t), POINTER, OPTIONAL :: Projector 59 INTEGER, OPTIONAL, POINTER :: NewMaskPerm(:) !< Mask the new variable set by the given MaskName when trying to define the interpolation. 60 LOGICAL, OPTIONAL :: KeepUnfoundNodes !< Do not disregard unfound nodes from projector 61 END SUBROUTINE InterpolateMeshToMeshQ 62 END INTERFACE 63!------------------------------------------------------------------------------- 64 65 ALLOCATE( FoundNodes(NewMesh % NumberOfNodes) ); FoundNodes=.FALSE. 66 67 IF(PRESENT(UnfoundNodes)) THEN 68 IF(ASSOCIATED(UnfoundNodes)) DEALLOCATE(UnfoundNodes) 69 ALLOCATE(UnfoundNodes(NewMesh % NumberOfNodes)) 70 END IF 71 72 ! In serial interpolation is simple 73 !------------------------------------ 74 IF ( ParEnv % PEs<=1 ) THEN 75 CALL InterpolateMeshToMeshQ( OldMesh, NewMesh, OldVariables, & 76 NewVariables, UseQuadrantTree, Projector, MaskName, FoundNodes ) 77 78 IF( InfoActive(20) ) THEN 79 n = COUNT(.NOT. FoundNodes ) 80 CALL Info('InterpolateMeshToMesh','Number of unfound nodes in serial: '//TRIM(I2S(n))) 81 END IF 82 83 IF(PRESENT(UnfoundNodes)) UnfoundNodes = .NOT. FoundNodes 84 RETURN 85 END IF 86 87 ! Passive coordinate is needed also here in order not to use that direction 88 ! for the bounding box checks. 89 !--------------------------------------------------------------------------- 90 PassiveCoordinate = ListGetInteger( CurrentModel % Simulation, & 91 'Interpolation Passive Coordinate', Stat ) 92 IF (.NOT. Stat .AND. ASSOCIATED(CurrentModel % Solver)) THEN 93 PassiveCoordinate = ListGetInteger( CurrentModel % Solver % Values, & 94 'Interpolation Passive Coordinate', Stat ) 95 END IF 96 97 ! Interpolate within our own partition, flag the points we found: 98 ! --------------------------------------------------------------- 99 CALL InterpolateMeshToMeshQ( OldMesh, NewMesh, OldVariables, & 100 NewVariables, UseQuadrantTree, MaskName=MaskName, FoundNodes=FoundNodes ) 101 102 IF(PRESENT(UnfoundNodes)) UnfoundNodes = .NOT. FoundNodes 103 104 ! special case "all found": 105 !-------------------------- 106 n = COUNT(.NOT.FoundNodes); dn = n 107 108 IF( InfoActive(20) ) THEN 109 CALL Info('InterpolateMeshToMesh','Number of unfound nodes in own partition: '//TRIM(I2S(n))) 110 END IF 111 112 AL = .FALSE. 113 IF (.NOT.ASSOCIATED(ParEnv % Active) ) THEN 114 ALLOCATE(Parenv % Active(PArEnv % PEs)) 115 AL = .TRUE. 116 ParEnv % Active = .TRUE. 117 END IF 118 119 CALL SParActiveSUM(dn,2) 120 IF ( dn==0 ) RETURN 121 122 ! No use to continue even in parallel, since the OldMeshes are all the same! 123 IF( OldMesh % SingleMesh ) THEN 124 CALL Warn('InterpolateMeshToMesh','Could not find all dofs in single mesh: '//TRIM(I2S(NINT(dn)))) 125 RETURN 126 END IF 127 128 ! Exchange partition bounding boxes: 129 ! This is needed to eliminate the amount of data to send among partitions. 130 ! ------------------------------------------------------------------------ 131 myBB = HUGE(mybb(1)) 132 IF(OldMesh % NumberOfNodes /= 0) THEN 133 myBB(1) = MINVAL(OldMesh % Nodes % x) 134 myBB(2) = MINVAL(OldMesh % Nodes % y) 135 myBB(3) = MINVAL(OldMesh % Nodes % z) 136 myBB(4) = MAXVAL(OldMesh % Nodes % x) 137 myBB(5) = MAXVAL(OldMesh % Nodes % y) 138 myBB(6) = MAXVAL(OldMesh % Nodes % z) 139 140 eps2 = 0.1_dp * MAXVAL(myBB(4:6)-myBB(1:3)) 141 myBB(1:3) = myBB(1:3) - eps2 142 myBB(4:6) = myBB(4:6) + eps2 143 END IF 144 145 ALLOCATE(BB(6,ParEnv % PEs)) 146 DO i=1,ParEnv % PEs 147 IF ( Parenv % mype == i-1 .OR. .NOT. ParEnv % Active(i) ) CYCLE 148 proc = i-1 149 CALL MPI_BSEND( myBB, 6, MPI_DOUBLE_PRECISION, proc, & 150 999, ELMER_COMM_WORLD, ierr ) 151 END DO 152 DO i=1,COUNT(ParEnv % Active)-1 153 CALL MPI_RECV( myBB, 6, MPI_DOUBLE_PRECISION, MPI_ANY_SOURCE, & 154 999, ELMER_COMM_WORLD, status, ierr ) 155 proc = status(MPI_SOURCE) 156 BB(:,proc+1) = myBB 157 END DO 158 159 CALL CheckBuffer((n*(3 * 2)) + 2) !3 x double precision coord, 2 x count 160 161 IF ( n==0 ) THEN 162 ! We have found all nodes, nothing to do except sent the info to others! 163 !---------------------------------------------------------------------- 164 DEALLOCATE(BB) 165 DO i=1,ParEnv % PEs 166 IF ( Parenv % mype == i-1 .OR. .NOT. ParEnv % Active(i) ) CYCLE 167 proc = i-1 168 CALL MPI_BSEND( n, 1, MPI_INTEGER, proc, & 169 1001, ELMER_COMM_WORLD, ierr ) 170 END DO 171 ELSE 172 ! Extract nodes that we didn't find from our own partition... 173 ! ------------------------------------------------------------ 174 ALLOCATE( Perm(n), nodes_x(n), nodes_y(n),nodes_z(n) ); Perm=0 175 j = 0 176 DO i=1,NewMesh % NumberOfNodes 177 IF ( FoundNodes(i) ) CYCLE 178 j = j + 1 179 perm(j) = i 180 nodes_x(j) = NewMesh % Nodes % x(i) 181 nodes_y(j) = NewMesh % Nodes % y(i) 182 nodes_z(j) = NewMesh % Nodes % z(i) 183 END DO 184 185 ! ...and ask those from others 186 ! ------------------------------- 187 ALLOCATE(ProcSend(ParEnv % PEs)) 188 DO i=1,ParEnv % PEs 189 IF ( Parenv % mype == i-1 .OR. .NOT. ParEnv % Active(i) ) CYCLE 190 proc = i-1 191 192 ! extract those of the missing nodes that are within the other 193 ! partions bounding box: 194 ! ------------------------------------------------------------ 195 myBB = BB(:,i) 196 npart = 0 197 DO j=1,n 198 IF ( ( nodes_x(j)<myBB(1) .OR. nodes_x(j)>myBB(4) ) .AND. PassiveCoordinate /= 1 ) CYCLE 199 IF ( ( nodes_y(j)<myBB(2) .OR. nodes_y(j)>myBB(5) ) .AND. PassiveCoordinate /= 2 ) CYCLE 200 IF ( ( nodes_z(j)<myBB(3) .OR. nodes_z(j)>myBB(6) ) .AND. PassiveCoordinate /= 3 ) CYCLE 201 npart = npart+1 202 END DO 203 ProcSend(proc+1) % n = npart 204 IF ( npart>0 ) THEN 205 ALLOCATE( xpart(npart),ypart(npart),zpart(npart),ProcSend(proc+1) % perm(npart) ) 206 npart = 0 207 DO j=1,n 208 IF ( ( nodes_x(j)<myBB(1) .OR. nodes_x(j)>myBB(4) ) .AND. PassiveCoordinate /= 1 ) CYCLE 209 IF ( ( nodes_y(j)<myBB(2) .OR. nodes_y(j)>myBB(5) ) .AND. PassiveCoordinate /= 2 ) CYCLE 210 IF ( ( nodes_z(j)<myBB(3) .OR. nodes_z(j)>myBB(6) ) .AND. PassiveCoordinate /= 3 ) CYCLE 211 npart=npart+1 212 ProcSend(proc+1) % perm(npart)=j 213 xpart(npart) = Nodes_x(j) 214 ypart(npart) = Nodes_y(j) 215 zpart(npart) = Nodes_z(j) 216 END DO 217 END IF 218 219 ! send count... 220 ! ------------- 221 CALL MPI_BSEND( npart, 1, MPI_INTEGER, proc, & 222 1001, ELMER_COMM_WORLD, ierr ) 223 224 IF ( npart==0 ) CYCLE 225 226 ! ...and points 227 ! ------------- 228 CALL MPI_BSEND( xpart, npart, MPI_DOUBLE_PRECISION, proc, & 229 1002, ELMER_COMM_WORLD, ierr ) 230 CALL MPI_BSEND( ypart, npart, MPI_DOUBLE_PRECISION, proc, & 231 1003, ELMER_COMM_WORLD, ierr ) 232 CALL MPI_BSEND( zpart, npart, MPI_DOUBLE_PRECISION, proc, & 233 1004, ELMER_COMM_WORLD, ierr ) 234 235 DEALLOCATE(xpart,ypart,zpart) 236 END DO 237 DEALLOCATE(nodes_x,nodes_y,nodes_z,BB) 238 END IF 239 240 241 ! receive points from others: 242 ! ---------------------------- 243 ALLOCATE(ProcRecv(Parenv % Pes)) 244 DO i=1,COUNT(ParEnv % Active)-1 245 CALL MPI_RECV( n, 1, MPI_INTEGER, MPI_ANY_SOURCE, & 246 1001, ELMER_COMM_WORLD, status, ierr ) 247 248 proc = status(MPI_SOURCE) 249 ProcRecv(proc+1) % n = n 250 251 IF ( n<=0 ) CYCLE 252 253 ALLOCATE(ProcRecv(proc+1) % Nodes_x(n), & 254 ProcRecv(proc+1) % Nodes_y(n),ProcRecv(proc+1) % Nodes_z(n)) 255 256 CALL MPI_RECV( ProcRecv(proc+1) % nodes_x, n, MPI_DOUBLE_PRECISION, proc, & 257 1002, ELMER_COMM_WORLD, status, ierr ) 258 CALL MPI_RECV( ProcRecv(proc+1) % nodes_y, n, MPI_DOUBLE_PRECISION, proc, & 259 1003, ELMER_COMM_WORLD, status, ierr ) 260 CALL MPI_RECV( ProcRecv(proc+1) % nodes_z, n, MPI_DOUBLE_PRECISION, proc, & 261 1004, ELMER_COMM_WORLD, status, ierr ) 262 END DO 263 264 ! Count variables and received nodes, and check MPI buffer is 265 ! sufficiently large: 266 ! ----------------------------------------------------------- 267 Var => OldVariables 268 nvars = 0 269 DO WHILE(ASSOCIATED(Var)) 270 nvars = nvars + 1 271 Var => Var % Next 272 END DO 273 274 maxrecv = 0 275 DO i=1,SIZE(ProcRecv) 276 maxrecv = MAX(maxrecv, ProcRecv(i) % n) 277 END DO 278 279 !For each node, we send a single integer perm and 280 !a real(dp) per variable. Also sending two counts 281 CALL CheckBuffer(SIZE(ProcRecv) * maxrecv * ((2 * nvars) + 1) + 2) 282 283 ! Check the received points and extract values for the to-be-interpolated- 284 ! variables, if we have the points within our domain: 285 ! ------------------------------------------------------------------------ 286 DO i=1,ParEnv % PEs 287 IF ( Parenv % mype == i-1 .OR. .NOT. ParEnv % Active(i) ) CYCLE 288 289 proc = i-1 290 n = ProcRecv(i) % n 291 292 IF ( n==0 ) THEN 293 CALL MPI_BSEND( n, 1, MPI_INTEGER, proc, & 294 2001, ELMER_COMM_WORLD, ierr ) 295 CYCLE 296 END IF 297 298 ! Construct temporary mesh structure for the received points: 299 ! ----------------------------------------------------------- 300 Nmesh => AllocateMesh() 301 Nmesh % Nodes % x => ProcRecv(i) % nodes_x 302 Nmesh % Nodes % y => ProcRecv(i) % nodes_y 303 Nmesh % Nodes % z => ProcRecv(i) % nodes_z 304 Nmesh % NumberOfNodes = n 305 306 ALLOCATE(nperm(n)) 307 DO j=1,n 308 nPerm(j)=j 309 END DO 310 311 Var => OldVariables 312 nvars = 0 313 DO WHILE(ASSOCIATED(Var)) 314 IF ( Var % DOFs==1 .AND. ASSOCIATED(Var % Perm).AND..NOT.Var % Secondary ) THEN 315 ALLOCATE(store(n)); store=0 316 nvars = nvars+1 317 CALL VariableAdd(nMesh % Variables,nMesh,Var % Solver, & 318 Var % Name,1,store,nperm ) 319 IF ( ASSOCIATED(Var % PrevValues) ) THEN 320 j = SIZE(Var % PrevValues,2) 321 nvars = nvars+j 322 Nvar => VariableGet( Nmesh % Variables,Var % Name,ThisOnly=.TRUE.) 323 ALLOCATE(Nvar % PrevValues(n,j)) 324 END IF 325 END IF 326 Var => Var % Next 327 END DO 328 329 ! try interpolating values for the points: 330 ! ---------------------------------------- 331 ALLOCATE( FoundNodesPar(n) ); FoundNodesPar=.FALSE. 332 333 CALL InterpolateMeshToMeshQ( OldMesh, nMesh, OldVariables, & 334 nMesh % Variables, UseQuadrantTree, MaskName=MaskName, FoundNodes=FoundNodesPar ) 335 336 nfound = COUNT(FoundNodesPar) 337 338 CALL MPI_BSEND( nfound, 1, MPI_INTEGER, proc, & 339 2001, ELMER_COMM_WORLD, ierr ) 340 341 ! send interpolated values back to the owner: 342 ! ------------------------------------------- 343 IF ( nfound>0 ) THEN 344 ALLOCATE(vstore(nfound,nvars), vperm(nfound)); vstore=0 345 k = 0 346 DO j=1,n 347 IF ( .NOT.FoundNodesPar(j)) CYCLE 348 k = k + 1 349 vperm(k) = j 350 Var => OldVariables 351 nvars = 0 352 DO WHILE(ASSOCIATED(Var)) 353 IF ( Var % DOFs==1 .AND. ASSOCIATED(Var % Perm).AND..NOT.Var % Secondary) THEN 354 Nvar => VariableGet( Nmesh % Variables,Var % Name,ThisOnly=.TRUE.) 355 nvars=nvars+1 356 vstore(k,nvars)=Nvar % Values(j) 357 IF ( ASSOCIATED(Var % PrevValues) ) THEN 358 DO l=1,SIZE(Var % PrevValues,2) 359 nvars = nvars+1 360 vstore(k,nvars)=Nvar % PrevValues(j,l) 361 END DO 362 END IF 363 END IF 364 Var => Var % Next 365 END DO 366 END DO 367 368 CALL MPI_BSEND( vperm, nfound, MPI_INTEGER, proc, & 369 2002, ELMER_COMM_WORLD, status, ierr ) 370 371 DO j=1,nvars 372 CALL MPI_BSEND( vstore(:,j), nfound,MPI_DOUBLE_PRECISION, proc, & 373 2002+j, ELMER_COMM_WORLD,ierr ) 374 END DO 375 376 DEALLOCATE(vstore, vperm) 377 END IF 378 379 !Pointers to ProcRev, deallocated automatically below 380 NULLIFY(Nmesh % Nodes % x,& 381 Nmesh % Nodes % y,& 382 Nmesh % Nodes % z) 383 384 CALL ReleaseMesh(Nmesh) 385 DEALLOCATE(FoundNodesPar, Nmesh) 386 END DO 387 DEALLOCATE(ProcRecv) 388 389 ! Receive interpolated values: 390 ! ---------------------------- 391 DO i=1,COUNT(ParEnv % Active)-1 392 393 ! recv count: 394 ! ----------- 395 CALL MPI_RECV( n, 1, MPI_INTEGER, MPI_ANY_SOURCE, & 396 2001, ELMER_COMM_WORLD, status, ierr ) 397 398 proc = status(MPI_SOURCE) 399 IF ( n<=0 ) THEN 400 IF ( ALLOCATED(ProcSend) ) THEN 401 IF ( ALLOCATED(ProcSend(proc+1) % Perm)) & 402 DEALLOCATE(ProcSend(proc+1) % Perm) 403 END IF 404 CYCLE 405 END IF 406 407 ALLOCATE(astore(n),vperm(n)) 408 409 ! recv permutation (where in the original array the 410 ! points the partition found are): 411 ! -------------------------------------------------- 412 CALL MPI_RECV( vperm, n, MPI_INTEGER, proc, & 413 2002, ELMER_COMM_WORLD, status, ierr ) 414 415 !Mark nodes as found 416 DO j=1,n 417 k=perm(ProcSend(proc+1) % Perm(vperm(j))) 418 FoundNodes(k) = .TRUE. 419 END DO 420 421 ! recv values and store: 422 ! ---------------------- 423 Var => OldVariables 424 nvars=0 425 DO WHILE(ASSOCIATED(Var)) 426 IF ( Var % DOFs==1 .AND. ASSOCIATED(Var % Perm) .AND..NOT.Var % Secondary ) THEN 427 428 nvars=nvars+1 429 CALL MPI_RECV( astore, n, MPI_DOUBLE_PRECISION, proc, & 430 2002+nvars, ELMER_COMM_WORLD, status, ierr ) 431 432 Nvar => VariableGet( NewMesh % Variables,Var % Name,ThisOnly=.TRUE.) 433 434 IF ( ASSOCIATED(Nvar) ) THEN 435 DO j=1,n 436 k=perm(ProcSend(proc+1) % Perm(vperm(j))) 437 IF ( Nvar % perm(k)>0 ) & 438 Nvar % Values(Nvar % Perm(k)) = astore(j) 439 END DO 440 END IF 441 442 IF ( ASSOCIATED(Var % PrevValues) ) THEN 443 DO l=1,SIZE(Var % PrevValues,2) 444 nvars=nvars+1 445 CALL MPI_RECV( astore, n, MPI_DOUBLE_PRECISION, proc, & 446 2002+nvars, ELMER_COMM_WORLD, status, ierr ) 447 448 IF ( ASSOCIATED(Nvar) ) THEN 449 DO j=1,n 450 k=perm(ProcSend(proc+1) % Perm(vperm(j))) 451 IF ( Nvar % perm(k)>0 ) & 452 Nvar % PrevValues(Nvar % Perm(k),l) = astore(j) 453 END DO 454 END IF 455 END DO 456 END IF 457 END IF 458 Var => Var % Next 459 END DO 460 DEALLOCATE(astore,vperm,ProcSend(proc+1) % perm) 461 END DO 462 IF ( ALLOCATED(Perm) ) DEALLOCATE(Perm,ProcSend) 463 464 CALL MPI_BARRIER(ParEnv % ActiveComm,ierr) 465 466 IF(AL) THEN 467 DEALLOCATE(Parenv % Active) 468 ParEnv % Active => NULL() 469 END IF 470 471 n = COUNT(.NOT. FoundNodes ) 472 CALL Info('InterpolateMeshToMesh',& 473 'Number of unfound nodes in all partitions: '//TRIM(I2S(n)),Level=6) 474 475 IF(PRESENT(UnfoundNodes)) UnfoundNodes = .NOT. FoundNodes 476 DEALLOCATE( FoundNodes ) 477 478 479CONTAINS 480 481!------------------------------------------------------------------------------ 482 FUNCTION AllocateMesh() RESULT(Mesh) 483!------------------------------------------------------------------------------ 484 TYPE(Mesh_t), POINTER :: Mesh 485!------------------------------------------------------------------------------ 486 INTEGER :: istat 487 488 ALLOCATE( Mesh, STAT=istat ) 489 IF ( istat /= 0 ) & 490 CALL Fatal( 'AllocateMesh', 'Unable to allocate a few bytes of memory?' ) 491 492! Nothing computed on this mesh yet! 493! ---------------------------------- 494 Mesh % SavesDone = 0 495 Mesh % OutputActive = .FALSE. 496 497 Mesh % AdaptiveDepth = 0 498 Mesh % Changed = .FALSE. ! TODO: Change this sometime 499 500 Mesh % Stabilize = .FALSE. 501 502 Mesh % Variables => NULL() 503 Mesh % Parent => NULL() 504 Mesh % Child => NULL() 505 Mesh % Next => NULL() 506 Mesh % RootQuadrant => NULL() 507 Mesh % Elements => NULL() 508 Mesh % Edges => NULL() 509 Mesh % Faces => NULL() 510 Mesh % Projector => NULL() 511 Mesh % NumberOfEdges = 0 512 Mesh % NumberOfFaces = 0 513 Mesh % NumberOfNodes = 0 514 Mesh % NumberOfBulkElements = 0 515 Mesh % NumberOfBoundaryElements = 0 516 517 Mesh % MaxFaceDOFs = 0 518 Mesh % MaxEdgeDOFs = 0 519 Mesh % MaxBDOFs = 0 520 Mesh % MaxElementDOFs = 0 521 Mesh % MaxElementNodes = 0 522 523 Mesh % ViewFactors => NULL() 524 525 ALLOCATE( Mesh % Nodes, STAT=istat ) 526 IF ( istat /= 0 ) & 527 CALL Fatal( 'AllocateMesh', 'Unable to allocate a few bytes of memory?' ) 528 NULLIFY( Mesh % Nodes % x ) 529 NULLIFY( Mesh % Nodes % y ) 530 NULLIFY( Mesh % Nodes % z ) 531 Mesh % Nodes % NumberOfNodes = 0 532 533 Mesh % ParallelInfo % NumberOfIfDOFs = 0 534 NULLIFY( Mesh % ParallelInfo % GlobalDOFs ) 535 NULLIFY( Mesh % ParallelInfo % INTERFACE ) 536 NULLIFY( Mesh % ParallelInfo % NeighbourList ) 537 538 END FUNCTION AllocateMesh 539!------------------------------------------------------------------------------- 540END SUBROUTINE InterpolateMeshToMesh 541!------------------------------------------------------------------------------- 542 543 544!------------------------------------------------------------------------------ 545!> Interpolates values of all variables from a mesh associated with 546!> the old model to the mesh of the new model. 547!------------------------------------------------------------------------------ 548 SUBROUTINE InterpolateMeshToMeshQ( OldMesh, NewMesh, OldVariables, NewVariables, & 549 UseQuadrantTree, Projector, MaskName, FoundNodes, NewMaskPerm, KeepUnfoundNodes ) 550!------------------------------------------------------------------------------ 551 USE DefUtils 552!------------------------------------------------------------------------------- 553 TYPE(Mesh_t), TARGET :: OldMesh !< Old mesh structure 554 TYPE(Mesh_t), TARGET :: NewMesh !< New mesh structure 555 TYPE(Variable_t), POINTER, OPTIONAL :: OldVariables !< Old model variable structure 556 TYPE(Variable_t), POINTER, OPTIONAL :: NewVariables !< New model variable structure. NewVariables defines the variables to be interpolated 557 LOGICAL, OPTIONAL :: UseQuadrantTree !< If true use the RootQuadrant of the old mesh in interpolation. 558 TYPE(Projector_t), POINTER, OPTIONAL :: Projector !< Use projector between meshes for interpolation, if available 559 CHARACTER(LEN=*),OPTIONAL :: MaskName !< Mask the old variable set by the given MaskName when trying to define the interpolation. 560 LOGICAL, OPTIONAL :: FoundNodes(:) !< List of nodes where the interpolation was a success 561 INTEGER, OPTIONAL, POINTER :: NewMaskPerm(:) !< Mask the new variable set by the given MaskName when trying to define the interpolation. 562 LOGICAL, OPTIONAL :: KeepUnfoundNodes !< Do not disregard unfound nodes from projector 563!------------------------------------------------------------------------------ 564 INTEGER :: dim 565 TYPE(Nodes_t) :: ElementNodes 566 INTEGER :: nBulk, i, j, k, l, n, np, bf_id, QTreeFails, TotFails, FoundCnt 567 REAL(KIND=dp), DIMENSION(3) :: Point 568 INTEGER, POINTER :: Indexes(:) 569 REAL(KIND=dp), DIMENSION(3) :: LocalCoordinates 570 TYPE(Variable_t), POINTER :: OldSol, NewSol, Var 571 INTEGER, POINTER :: OldPerm(:) 572 REAL(KIND=dp), POINTER :: OldValue(:), NewValue(:), ElementValues(:) 573 TYPE(Quadrant_t), POINTER :: LeafQuadrant 574 TYPE(Element_t),POINTER :: Element, Parent 575 576 REAL(KIND=dp), ALLOCATABLE :: Basis(:),Vals(:),dVals(:,:), & 577 RotWBasis(:,:), WBasis(:,:) 578 REAL(KIND=dp) :: BoundingBox(6), detJ, u,v,w,s,val,rowsum, F(3,3), G(3,3) 579 580 LOGICAL :: UseQTree, TryQTree, Stat, UseProjector, EdgeBasis, PiolaT, Parallel, & 581 TryLinear, KeepUnfoundNodesL 582 TYPE(Quadrant_t), POINTER :: RootQuadrant 583 584 INTEGER, POINTER CONTIG :: Rows(:), Cols(:) 585 INTEGER, POINTER :: Diag(:) 586 587 TYPE Epntr_t 588 TYPE(Element_t), POINTER :: Element 589 END TYPE Epntr_t 590 591 TYPE(Epntr_t), ALLOCATABLE :: ElemPtrs(:) 592 593 INTEGER, ALLOCATABLE:: RInd(:) 594 LOGICAL :: Found, EpsAbsGiven,EpsRelGiven, MaskExists, CylProject, ProjectorAllocated 595 INTEGER :: eps_tries, nrow, PassiveCoordinate 596 REAL(KIND=dp) :: eps1 = 0.1, eps2, eps_global, eps_local, eps_basis,eps_numeric 597 REAL(KIND=dp), POINTER CONTIG :: Values(:) 598 REAL(KIND=dp), POINTER :: LocalU(:), LocalV(:), LocalW(:) 599 600 TYPE(Nodes_t), SAVE :: Nodes 601 602 !$OMP THREADPRIVATE(eps1,Nodes) 603 604!------------------------------------------------------------------------------ 605 606 Parallel = (ParEnv % PEs > 1) 607 608 FoundCnt = 0 609 IF ( OldMesh % NumberOfNodes == 0 ) RETURN 610! 611! If projector argument given, search for existing 612! projector matrix, or generate new projector, if 613! not already there: 614! ------------------------------------------------ 615 IF ( PRESENT(Projector) ) THEN 616 Projector => NewMesh % Projector 617 618 DO WHILE( ASSOCIATED( Projector ) ) 619 IF ( ASSOCIATED(Projector % Mesh, OldMesh) ) THEN 620 CALL Info('InterpolateMesh2Mesh','Applying exiting projector in interpolation',Level=12) 621 IF ( PRESENT(OldVariables) ) CALL ApplyProjector() 622 RETURN 623 END IF 624 Projector => Projector % Next 625 END DO 626 627 n = NewMesh % NumberOfNodes 628 ALLOCATE( LocalU(n), LocalV(n), LocalW(n), ElemPtrs(n) ) 629 DO i=1,n 630 NULLIFY( ElemPtrs(i) % Element ) 631 END DO 632 END IF 633! 634! Check if using the spatial division hierarchy for the search: 635! ------------------------------------------------------------- 636 637 RootQuadrant => OldMesh % RootQuadrant 638 dim = CoordinateSystemDimension() 639 640 IF ( .NOT. PRESENT( UseQuadrantTree ) ) THEN 641 UseQTree = .TRUE. 642 ELSE 643 UseQTree = UseQuadrantTree 644 ENDIF 645 646 IF ( UseQTree ) THEN 647 IF ( .NOT.ASSOCIATED( RootQuadrant ) ) THEN 648 BoundingBox(1) = MINVAL(OldMesh % Nodes % x) 649 BoundingBox(2) = MINVAL(OldMesh % Nodes % y) 650 BoundingBox(3) = MINVAL(OldMesh % Nodes % z) 651 BoundingBox(4) = MAXVAL(OldMesh % Nodes % x) 652 BoundingBox(5) = MAXVAL(OldMesh % Nodes % y) 653 BoundingBox(6) = MAXVAL(OldMesh % Nodes % z) 654 655 eps2 = 0.1_dp * MAXVAL(BoundingBox(4:6)-BoundingBox(1:3)) 656 BoundingBox(1:3) = BoundingBox(1:3) - eps2 657 BoundingBox(4:6) = BoundingBox(4:6) + eps2 658 659 CALL BuildQuadrantTree( OldMesh,BoundingBox,OldMesh % RootQuadrant) 660 RootQuadrant => OldMesh % RootQuadrant 661 END IF 662 END IF 663 664! Use mask or not 665!--------------------------------------- 666 MaskExists = PRESENT( MaskName ) 667 668!------------------------------------------------------------------------------ 669 670 n = OldMesh % MaxElementNodes 671 ALLOCATE( ElementNodes % x(n), ElementNodes % y(n), & 672 ElementNodes % z(n), ElementValues(n) ) 673 674 eps_global = ListGetConstReal( CurrentModel % Simulation, & 675 'Interpolation Global Epsilon', Stat) 676 IF(.NOT. Stat) eps_global = 2.0d-10 677 678 eps_local = ListGetConstReal( CurrentModel % Simulation, & 679 'Interpolation Local Epsilon', Stat ) 680 IF(.NOT. Stat) eps_local = 1.0d-10 681 682 eps_tries = ListGetInteger( CurrentModel % Simulation, & 683 'Interpolation Max Iterations', Stat ) 684 IF(.NOT. Stat) eps_tries = 12 685 686 eps_numeric = ListGetConstReal( CurrentModel % Simulation, & 687 'Interpolation Numeric Epsilon', Stat) 688 IF(.NOT. Stat) eps_numeric = 1.0e-10 689 690 PassiveCoordinate = ListGetInteger( CurrentModel % Simulation, & 691 'Interpolation Passive Coordinate', Stat ) 692 IF (.NOT. Stat .AND. ASSOCIATED(CurrentModel % Solver)) THEN 693 PassiveCoordinate = ListGetInteger( CurrentModel % Solver % Values, & 694 'Interpolation Passive Coordinate', Stat ) 695 END IF 696 697 CylProject = ListGetLogical( CurrentModel % Simulation, & 698 'Interpolation Cylindric', Stat ) 699 IF (.NOT. Stat .AND. ASSOCIATED(CurrentModel % Solver)) THEN 700 CylProject = ListGetLogical( CurrentModel % Solver % Values, & 701 'Interpolation Cylindric', Stat ) 702 END IF 703 704 QTreeFails = 0 705 TotFails = 0 706 707 EdgeBasis = .FALSE. 708 IF (ASSOCIATED(CurrentModel % Solver)) & 709 EdgeBasis = ListGetLogical(CurrentModel % Solver % Values,'Edge Basis',Found) 710 711 PiolaT = .FALSE. 712 IF (EdgeBasis) THEN 713 PiolaT = ListGetLogical(CurrentModel % Solver % Values,'Use Piola Transform',Found) 714 END IF 715 716 TryLinear = ListGetLogical( CurrentModel % Simulation, 'Try Linear Search If Qtree Fails', Found) 717 IF(.NOT.Found) TryLinear = .TRUE. 718 719 IF ( PRESENT(KeepUnfoundNodes) ) THEN 720 KeepUnfoundNodesL = KeepUnfoundNodes 721 ELSE 722 KeepUnfoundNodesL = .TRUE. 723 END IF 724 725 FoundCnt = 0 726!------------------------------------------------------------------------------ 727! Loop over all nodes in the new mesh 728!------------------------------------------------------------------------------ 729 DO i=1,NewMesh % NumberOfNodes 730!------------------------------------------------------------------------------ 731 732 ! Only get the variable for the requested nodes 733 IF( PRESENT( NewMaskPerm ) ) THEN 734 IF( NewMaskPerm(i) == 0 ) CYCLE 735 END IF 736 737 Point(1) = NewMesh % Nodes % x(i) 738 Point(2) = NewMesh % Nodes % y(i) 739 Point(3) = NewMesh % Nodes % z(i) 740 741 IF( PassiveCoordinate /= 0 ) THEN 742 Point(PassiveCoordinate) = 0.0_dp 743 END IF 744 745 IF( CylProject ) THEN 746 Point(1) = SQRT( Point(1)**2 + Point(2)**2 ) 747 Point(2) = Point(3) 748 Point(3) = 0.0_dp 749 END IF 750 751!------------------------------------------------------------------------------ 752! Find in which old mesh bulk element the point belongs to 753!------------------------------------------------------------------------------ 754 Found = .FALSE. 755 TryQTree = ASSOCIATED(RootQuadrant) .AND. UseQTree 756 757 IF( TryQTree ) THEN 758!------------------------------------------------------------------------------ 759! Find the last existing quadrant that the point belongs to 760!------------------------------------------------------------------------------ 761 Element => NULL() 762 CALL FindLeafElements(Point, dim, RootQuadrant, LeafQuadrant) 763 764 IF ( ASSOCIATED(LeafQuadrant) ) THEN 765 ! Go through the bulk elements in the last ChildQuadrant 766 ! only. Try to find matching element with progressively 767 ! sloppier tests. Allow at most 100 % of slack: 768 ! ------------------------------------------------------- 769 Eps1 = eps_global 770 Eps2 = eps_local 771 772 DO j=1,eps_tries 773 DO k=1, LeafQuadrant % NElemsInQuadrant 774 Element => OldMesh % Elements(LeafQuadrant % Elements(k)) 775 776 IF( MaskExists ) THEN 777 bf_id = ListGetInteger( CurrentModel % Bodies(Element % BodyId) % Values, & 778 'Body Force', Found ) 779 IF( .NOT. Found ) CYCLE 780 IF(.NOT. ListCheckPresent( & 781 CurrentModel % BodyForces(bf_id) % Values,MaskName) ) CYCLE 782 END IF 783 784 Indexes => Element % NodeIndexes 785 n = Element % TYPE % NumberOfNodes 786 787 ElementNodes % x(1:n) = OldMesh % Nodes % x(Indexes) 788 ElementNodes % y(1:n) = OldMesh % Nodes % y(Indexes) 789 ElementNodes % z(1:n) = OldMesh % Nodes % z(Indexes) 790 791 Found = PointInElement( Element, ElementNodes, & 792 Point, LocalCoordinates, Eps1, Eps2, NumericEps=eps_numeric,EdgeBasis=PiolaT) 793 IF ( Found ) EXIT 794 END DO 795 IF ( Found ) EXIT 796 797 Eps1 = 10 * Eps1 798 Eps2 = 10 * Eps2 799 IF( Eps1 > 1.0_dp ) EXIT 800 END DO 801 END IF 802 END IF 803 804 IF( .NOT. TryQTree .OR. (.NOT. Found .AND. .NOT. Parallel .AND. TryLinear ) ) THEN 805 !------------------------------------------------------------------------------ 806 ! Go through all old mesh bulk elements 807 !------------------------------------------------------------------------------ 808 DO k=1,OldMesh % NumberOfBulkElements 809 Element => OldMesh % Elements(k) 810 811 n = Element % TYPE % NumberOfNodes 812 Indexes => Element % NodeIndexes 813 814 ElementNodes % x(1:n) = OldMesh % Nodes % x(Indexes) 815 ElementNodes % y(1:n) = OldMesh % Nodes % y(Indexes) 816 ElementNodes % z(1:n) = OldMesh % Nodes % z(Indexes) 817 818 Found = PointInElement( Element, ElementNodes, & 819 Point, LocalCoordinates ) 820 IF( Found ) THEN 821 IF( TryQTree ) QTreeFails = QtreeFails + 1 822 EXIT 823 END IF 824 END DO 825 END IF 826 827 IF (.NOT.Found) THEN 828 Element => NULL() 829 IF (.NOT. Parallel ) THEN 830 WRITE( Message,'(A,I0,A)' ) 'Point ',i,' was not found in any of the elements!' 831 CALL Info( 'InterpolateMeshToMesh', Message, Level=20 ) 832 TotFails = TotFails + 1 833 END IF 834 CYCLE 835 END IF 836 IF ( PRESENT(FoundNodes) ) FoundNodes(i) = .TRUE. 837 838!------------------------------------------------------------------------------ 839! 840! Found Element in OldModel: 841! --------------------------------- 842 843 IF ( PRESENT(Projector) ) THEN 844 FoundCnt = FoundCnt + 1 845 IF ( KeepUnfoundNodesL ) THEN 846 ElemPtrs(i) % Element => Element 847 LocalU(i) = LocalCoordinates(1) 848 LocalV(i) = LocalCoordinates(2) 849 LocalW(i) = LocalCoordinates(3) 850 ELSE 851 ElemPtrs(FoundCnt) % Element => Element 852 LocalU(FoundCnt) = LocalCoordinates(1) 853 LocalV(FoundCnt) = LocalCoordinates(2) 854 LocalW(FoundCnt) = LocalCoordinates(3) 855 END IF 856 END IF 857 858 IF ( .NOT.PRESENT(OldVariables) .OR. PRESENT(Projector) ) CYCLE 859!------------------------------------------------------------------------------ 860! 861! Go through all variables to be interpolated: 862! -------------------------------------------- 863 Var => OldVariables 864 DO WHILE( ASSOCIATED( Var ) ) 865 866 IF( SIZE( Var % Values ) == Var % DOFs ) THEN 867 Var => Var % Next 868 CYCLE 869 END IF 870 871 IF( Var % Secondary ) THEN 872 Var => Var % Next 873 CYCLE 874 END IF 875 876 IF ( Var % DOFs == 1 .AND. & 877 Var % Name(1:10) /= 'coordinate') THEN 878 879!------------------------------------------------------------------------------ 880! 881! Interpolate variable at Point in Element: 882! ------------------------------------------------ 883 884 NewSol => VariableGet( NewMesh % Variables, Var % Name, .TRUE. ) 885 IF ( .NOT. ASSOCIATED( NewSol ) ) THEN 886 Var => Var % Next 887 CYCLE 888 END IF 889 OldSol => VariableGet( OldMesh % Variables, Var % Name, .TRUE. ) 890 891 892 ! Check that the node was found in the old mesh: 893 ! ---------------------------------------------- 894 IF ( ASSOCIATED (Element) ) THEN 895 !------------------------------------------------------------------------------ 896! 897! Check for rounding errors: 898! -------------------------- 899 IF( OldSol % TYPE == Variable_on_nodes_on_elements ) THEN 900 Indexes => Element % DGIndexes 901 ELSE 902 Indexes => Element % NodeIndexes 903 END IF 904 905 906 IF ( ALL(OldSol % Perm(Indexes) > 0) ) THEN 907 IF ( NewSol % Perm(i) /= 0 ) THEN 908 ElementValues(1:n) = & 909 OldSol % Values(OldSol % Perm(Indexes)) 910 911 val = InterpolateInElement( Element, ElementValues, & 912 LocalCoordinates(1), LocalCoordinates(2), LocalCoordinates(3) ) 913 914 NewSol % Values(NewSol % Perm(i)) = val 915 916 IF ( ASSOCIATED( OldSol % PrevValues ) ) THEN 917 DO j=1,SIZE(OldSol % PrevValues,2) 918 ElementValues(1:n) = & 919 OldSol % PrevValues(OldSol % Perm(Indexes),j) 920 921 val = InterpolateInElement( Element, ElementValues, & 922 LocalCoordinates(1), LocalCoordinates(2), LocalCoordinates(3) ) 923 924 NewSol % PrevValues(NewSol % Perm(i),j) = val 925 END DO 926 END IF 927 END IF 928 END IF 929 ELSE 930 IF ( NewSol % Perm(i)/=0 ) NewValue(NewSol % Perm(i))=0.0_dp 931 END IF 932 933!------------------------------------------------------------------------------ 934 END IF 935 Var => Var % Next 936 END DO 937!------------------------------------------------------------------------------ 938 END DO 939 940 IF( .NOT. Parallel ) THEN 941 IF( QtreeFails > 0 ) THEN 942 WRITE( Message,'(A,I0)' ) 'Number of points not found in quadtree: ',QtreeFails 943 CALL Info( 'InterpolateMeshToMesh', Message ) 944 IF( TotFails == 0 ) THEN 945 CALL Info( 'InterpolateMeshToMesh','All nodes still found by N^2 dummy search!' ) 946 END IF 947 END IF 948 IF( TotFails == 0 ) THEN 949 CALL Info('InterpolateMeshToMesh','Found all nodes in the target mesh',Level=6) 950 ELSE 951 WRITE( Message,'(A,I0,A,I0,A)') 'Points not found: ',TotFails,' (found ',& 952 NewMesh % NumberOfNodes - TotFails,')' 953 CALL Warn( 'InterpolateMeshToMesh', Message ) 954 END IF 955 END IF 956 957!------------------------------------------------------------------------------ 958! 959! Construct mesh projector, if requested. Next time around 960! will use the existing projector to interpolate values: 961! --------------------------------------------------------- 962 IF ( PRESENT(Projector) ) THEN 963 964 IF ( KeepUnfoundNodesL ) THEN 965 n = NewMesh % NumberOfNodes 966 ELSE 967 n = FoundCnt 968 END IF 969 ALLOCATE( Basis(100),Vals(100), Indexes(100)) 970 971 ! The critical value of basis function that is accepted to the 972 ! projector. Note that the sum of weights is one, so this 973 ! we know the scale for this one. 974 eps_basis = ListGetConstReal( CurrentModel % Simulation, & 975 'Interpolation Basis Epsilon', Stat ) 976 IF(.NOT. Stat) eps_basis = 0.0d-12 977 978 979 ALLOCATE( Rows(n+1) ) 980 Rows(1) = 1 981 ProjectorAllocated = .FALSE. 982 983100 nrow = 1 984 985 DO i=1,n 986 987 IF(EdgeBasis.AND.ASSOCIATED(OldMesh % Parent)) THEN 988 Element => OldMesh % Parent % Faces(ElemPtrs(i) % Element % ElementIndex) 989 IF(ASSOCIATED(Element % BoundaryInfo)) THEN 990 Parent => Element % BoundaryInfo% Left 991 IF (ASSOCIATED(Parent)) THEN 992 k = Element % TYPE % NumberOfNodes 993 np = Parent % TYPE % NumberOfNodes 994 END IF 995 END IF 996 ELSE 997 Element => ElemPtrs(i) % Element 998 END IF 999 Found = ASSOCIATED( Element ) 1000 1001 IF( .NOT. Found ) THEN 1002 ! It seems unnecessary to make a matrix entry in case no target element is found! 1003 IF(.FALSE.) THEN 1004 IF( ProjectorAllocated ) THEN 1005 Cols(nrow) = 1 1006 Values(nrow) = 0._dp 1007 END IF 1008 nrow = nrow + 1 1009 END IF 1010 ELSE 1011 1012 u = LocalU(i) 1013 v = LocalV(i) 1014 w = LocalW(i) 1015 1016 IF(EdgeBasis) THEN 1017 CALL GetElementNodes(Nodes,Element) 1018 ELSE 1019 CALL GetElementNodes(Nodes,Element,UMesh=OldMesh) 1020 END IF 1021 1022 k = GetElementDOFs( Indexes, Element, NotDG=ASSOCIATED(CurrentModel % Solver)) 1023 1024 np = GetElementNOFNodes(Element) 1025 IF (ANY(Indexes(1:np)>Element % NodeIndexes)) np=0 1026 1027 IF( EdgeBasis) THEN 1028 IF(.NOT.ALLOCATED(dVals)) THEN 1029 ALLOCATE(dVals(k,3),WBasis(k,3),RotWBasis(k,3)) 1030 ELSE IF(SIZE(dVals,1)<k) THEN 1031 DEALLOCATE(dVals,WBasis,RotWBasis) 1032 ALLOCATE(dVals(k,3),WBasis(k,3),RotWBasis(k,3)) 1033 END IF 1034 1035 IF(PiolaT) THEN 1036 stat = ElementInfo(Element,Nodes,u,v,w,detJ,Vals,EdgeBasis=WBasis ) 1037 ELSE 1038 stat = ElementInfo(Element,Nodes,u,v,w,detJ,Vals,dVals) 1039 CALL GetEdgeBasis(Element,WBasis,RotWBasis,Vals,dVals) 1040 END IF 1041 ELSE 1042 stat = ElementInfo(Element,Nodes,u,v,w,detJ,Vals) 1043 END IF 1044 1045 1046 rowsum = 0.0_dp 1047 DO j=1,k 1048! IF( ABS( vals(j) ) < eps_basis ) CYCLE 1049 IF(j<=np) rowsum = rowsum + vals(j) 1050 IF (.NOT. ProjectorAllocated) THEN 1051 IF (.NOT.EdgeBasis.OR.(EdgeBasis.AND.j<=np)) THEN 1052 nrow = nrow+1 1053 ELSE 1054 nrow = nrow+3 1055 END IF 1056 END IF 1057 END DO 1058 1059 1060 IF( ProjectorAllocated ) THEN 1061 DO j=1,k 1062! IF( ABS( vals(j) ) < eps_basis ) CYCLE 1063 IF(.NOT.EdgeBasis) Rind(Indexes(j)) = Rind(Indexes(j)) + 1 1064 1065 ! Always normalize the weights to one! 1066 IF (.NOT.EdgeBasis.OR.(EdgeBasis.AND.j<=np)) THEN 1067 Cols(nrow) = Indexes(j) 1068 Values(nrow) = vals(j) / rowsum 1069 nrow = nrow + 1 1070 ELSE 1071 Cols(nrow) = -Indexes(j) 1072 Values(nrow) = WBasis(j-np,1) 1073 nrow = nrow + 1 1074 Cols(nrow) = -Indexes(j) 1075 Values(nrow) = WBasis(j-np,2) 1076 nrow = nrow + 1 1077 Cols(nrow) = -Indexes(j) 1078 Values(nrow) = WBasis(j-np,3) 1079 nrow = nrow + 1 1080 END IF 1081 END DO 1082 END IF 1083 END IF 1084 1085 Rows(i+1) = nrow 1086 END DO 1087 1088 IF( .NOT. ProjectorAllocated ) THEN 1089 ALLOCATE( Cols(Rows(n+1)-1), Values(Rows(n+1)-1) ) 1090 Cols = 0 1091 Values = 0 1092 1093 ALLOCATE( Projector ) 1094 Projector % Matrix => AllocateMatrix() 1095 Projector % Matrix % NumberOfRows = n 1096 Projector % Matrix % Rows => Rows 1097 Projector % Matrix % Cols => Cols 1098 Projector % Matrix % Values => Values 1099 1100 Projector % Next => NewMesh % Projector 1101 NewMesh % Projector => Projector 1102 NewMesh % Projector % Mesh => OldMesh 1103 1104 IF( .NOT.EdgeBasis) THEN 1105 ALLOCATE(Rind(OldMesh % NumberOfNodes)); Rind = 0 1106 END IF 1107 1108 ProjectorAllocated = .TRUE. 1109 1110 GOTO 100 1111 END IF 1112 1113 DEALLOCATE( Basis, Vals, ElemPtrs, LocalU, LocalV, LocalW, Indexes ) 1114 1115! Store also the transpose of the projector: 1116! ------------------------------------------ 1117 Projector % TMatrix => NULL() 1118 IF(.NOT.EdgeBasis) THEN 1119 IF ( Found ) THEN 1120 n = OldMesh % NumberOfNodes 1121 ! Needed for some matrices 1122 n = MAX( n, MAXVAL( Projector % Matrix % Cols ) ) 1123 1124 ALLOCATE( Rows(n+1) ) 1125 Rows(1) = 1 1126 DO i=2,n+1 1127 Rows(i) = Rows(i-1)+Rind(i-1) 1128 END DO 1129 1130 ALLOCATE( Cols(Rows(n+1)-1), Values(Rows(n+1)-1) ) 1131 Projector % TMatrix => AllocateMatrix() 1132 Projector % TMatrix % NumberOfRows = n 1133 Projector % TMatrix % Rows => Rows 1134 Projector % TMatrix % Cols => Cols 1135 Projector % TMatrix % Values => Values 1136 1137 RInd = 0 1138 DO i=1,Projector % Matrix % NumberOfRows 1139 DO j=Projector % Matrix % Rows(i), Projector % Matrix % Rows(i+1)-1 1140 k = Projector % Matrix % Cols(j) 1141 l = Rows(k) + Rind(k) 1142 Rind(k) = Rind(k) + 1 1143 Cols(l) = i 1144 Values(l) = Projector % Matrix % Values(j) 1145 END DO 1146 END DO 1147 END IF 1148 1149 DEALLOCATE(Rind) 1150 END IF 1151 1152 IF ( PRESENT(OldVariables) ) CALL ApplyProjector 1153 END IF 1154 1155 DEALLOCATE( ElementNodes % x, ElementNodes % y, & 1156 ElementNodes % z, ElementValues ) 1157 1158CONTAINS 1159 1160 1161!------------------------------------------------------------------------------ 1162 SUBROUTINE ApplyProjector 1163!------------------------------------------------------------------------------ 1164 INTEGER :: i 1165 TYPE(Variable_t), POINTER :: Var 1166!------------------------------------------------------------------------------ 1167 Var => OldVariables 1168 DO WHILE( ASSOCIATED(Var) ) 1169 IF( SIZE( Var % Values ) == Var % DOFs ) THEN 1170 Var => Var % Next 1171 CYCLE 1172 END IF 1173 1174 IF( Var % Secondary ) THEN 1175 Var => Var % Next 1176 CYCLE 1177 END IF 1178 1179 IF ( Var % DOFs == 1 .AND. & 1180 Var % Name(1:10) /= 'coordinate') THEN 1181 1182 OldSol => VariableGet( OldMesh % Variables, Var % Name, .TRUE. ) 1183 NewSol => VariableGet( NewMesh % Variables, Var % Name, .TRUE. ) 1184 IF ( .NOT. (ASSOCIATED (NewSol) ) ) THEN 1185 Var => Var % Next 1186 CYCLE 1187 END IF 1188 1189 CALL CRS_ApplyProjector( Projector % Matrix, & 1190 OldSol % Values, OldSol % Perm, & 1191 NewSol % Values, NewSol % Perm ) 1192 1193 IF ( ASSOCIATED( OldSol % PrevValues ) ) THEN 1194 DO i=1,SIZE(OldSol % PrevValues,2) 1195 CALL CRS_ApplyProjector( Projector % Matrix, & 1196 OldSol % PrevValues(:,i), OldSol % Perm, & 1197 NewSol % PrevValues(:,i), NewSol % Perm ) 1198 END DO 1199 END IF 1200 END IF 1201 Var => Var % Next 1202 END DO 1203!------------------------------------------------------------------------------ 1204 END SUBROUTINE ApplyProjector 1205!------------------------------------------------------------------------------ 1206 END SUBROUTINE InterpolateMeshToMeshQ 1207!------------------------------------------------------------------------------ 1208 1209 1210 1211 !--------------------------------------------------------------------------- 1212 !> Create a projector for mapping between interfaces using the Galerkin method 1213 !> A temporal mesh structure with a node for each Gaussian integration point is 1214 !> created. The this projector matrix is transferred to a projector on the nodal 1215 !> coordinates. 1216 !--------------------------------------------------------------------------- 1217 FUNCTION WeightedProjector(BMesh2, BMesh1, InvPerm2, InvPerm1, & 1218 UseQuadrantTree, Repeating, AntiRepeating, PeriodicScale, & 1219 NodalJump ) & 1220 RESULT ( Projector ) 1221 !--------------------------------------------------------------------------- 1222 USE DefUtils 1223 IMPLICIT NONE 1224 1225 TYPE(Mesh_t), POINTER :: BMesh1, BMesh2 1226 REAL(KIND=dp) :: PeriodicScale 1227 INTEGER, POINTER :: InvPerm1(:), InvPerm2(:) 1228 LOGICAL :: UseQuadrantTree, Repeating, AntiRepeating 1229 TYPE(Matrix_t), POINTER :: Projector 1230 LOGICAL :: NodalJump 1231 !-------------------------------------------------------------------------- 1232 LOGICAL, ALLOCATABLE :: MirrorNode(:) 1233 TYPE(Matrix_t), POINTER :: GaussProjector 1234 TYPE(Nodes_t), POINTER :: GaussNodes, RealNodes, ElementNodes 1235 TYPE(Element_t), POINTER :: Element 1236 INTEGER :: i,j,k,l,n,p,q,NoNodes, NoGaussPoints,Indexes(100),& 1237 nodesize, totsize, eqindsize, RelOrder 1238 INTEGER, POINTER :: NodeIndexes(:),Rows(:),Cols(:) 1239 REAL(KIND=dp), POINTER :: Basis(:), Values(:) 1240 REAL(KIND=dp) :: u,v,w,val,detJ,weight,x 1241 TYPE(GaussIntegrationPoints_t) :: IntegStuff 1242 LOGICAL :: Stat, EdgeBasis,Found,AxisSym, PiolaT 1243 TYPE(Nodes_t), SAVE :: Nodes 1244 1245 REAL(KIND=dp) :: vq(3),wq(3),f(3,3),g(3,3) 1246 REAL(KIND=dp), ALLOCATABLE ::WBasis(:,:),RotWbasis(:,:),dBasisdx(:,:) 1247 1248 INTEGER :: ind,np,qq,pp 1249 INTEGER, ALLOCATABLE :: Eqind(:), xPerm(:) 1250 1251 1252 CALL Info('WeightedProjector','Creating Galerkin projector between two boundaries',Level=7) 1253 1254 ALLOCATE( GaussNodes, ElementNodes ) 1255 RealNodes => Bmesh1 % Nodes 1256 NoNodes = Bmesh1 % NumberOfNodes 1257 1258 EdgeBasis = .FALSE. 1259 IF (ASSOCIATED(CurrentModel % Solver)) THEN 1260 EdgeBasis = ListGetLogical(CurrentModel % Solver % Values,'Edge Basis',Found) 1261 END IF 1262 1263 PiolaT = .FALSE. 1264 IF( EdgeBasis ) THEN 1265 PiolaT = ListGetLogical( CurrentModel % Solver % Values, 'Use Piola Transform', Found) 1266 CALL Info('weightedProjector','Accounting for edge elements in projector.',Level=7) 1267 END IF 1268 1269 RelOrder = ListGetInteger( CurrentModel % Solver % Values, & 1270 'Projector Relative Integration Order', Found, minv=-1,maxv=1) 1271 1272 ! Calculate the total number of Gaussian integration points 1273 ! and allocate space for the node structures. 1274 !---------------------------------------------------------- 1275 NoGaussPoints = 0 1276 DO i=1, BMesh1 % NumberOfBulkElements 1277 Element => BMesh1 % Elements(i) 1278 IntegStuff = GaussPoints( Element, RelOrder=RelOrder, EdgeBasis=PiolaT ) 1279 NoGaussPoints = NoGaussPoints + IntegStuff % n 1280 END DO 1281 1282 WRITE( Message,'(A,I0,A,I0)') 'Number of nodes and gauss points:'& 1283 ,NoNodes,' and ',NoGaussPoints 1284 CALL Info('WeightedProjector',Message,Level=10) 1285 1286 ALLOCATE( GaussNodes % x(NoGaussPoints), GaussNodes % y(NoGaussPoints), GaussNodes % z(NoGaussPoints)) 1287 1288 ! Change the local coordinates of the BMesh2 to match to corresponding faces: 1289 ! --------------------------------------------------------------------------- 1290 IF(EdgeBasis) THEN 1291 ALLOCATE(xPerm(Bmesh2 % Parent % NumberofNodes)); xPerm=0 1292 DO i=1,SIZE(InvPerm2) 1293 xPerm(InvPerm2(i)) = i 1294 END DO 1295 1296 DO i=1, BMesh2 % NumberOfBulkElements 1297 Element => BMesh2 % Parent % Faces(BMesh2 % Elements(i) % ElementIndex) 1298 BMesh2 % Elements(i) % NodeIndexes = xPerm(Element % NodeIndexes) 1299 END DO 1300 END IF 1301 1302 AxisSym = .FALSE. 1303 IF ( CurrentCoordinateSystem() == AxisSymmetric .OR. & 1304 CurrentCoordinateSystem() == CylindricSymmetric ) THEN 1305 IF( NodalJump ) THEN 1306 AxisSym = .TRUE. 1307 ELSE IF (ASSOCIATED(CurrentModel % Solver)) THEN 1308 AxisSym = ListGetLogical(CurrentModel % Solver % Values,'Projector Metrics',Found) 1309 END IF 1310 IF( AxisSym ) CALL Info('weightedProjector','Projector will be weighted for axi symmetry',Level=7) 1311 END IF 1312 1313 1314 nodesize = BMesh1 % Parent % NumberOfNodes 1315 totsize = BMesh1 % Parent % NumberOfNodes + BMesh1 % Parent % NumberOfEdges 1316 IF(ASSOCIATED(CurrentModel % Solver)) THEN 1317 totsize = CurrentModel % Solver % Matrix % NumberOfRows / & 1318 CurrentModel % Solver % Variable % Dofs 1319 END IF 1320 1321 IF(EdgeBasis) THEN 1322 n = BMesh1 % Parent % MaxElementDOFs 1323 ALLOCATE( ElementNodes % x(n), ElementNodes % y(n), ElementNodes % z(n), & 1324 Basis(n), dBasisdx(n,3), WBasis(n,3), RotWBasis(n,3) ) 1325 eqindsize = totsize 1326 ELSE 1327 n = BMesh1 % MaxElementNodes 1328 ALLOCATE( ElementNodes % x(n), ElementNodes % y(n), ElementNodes % z(n), Basis(n) ) 1329 eqindsize = BMesh1 % NumberOfNodes 1330 END IF 1331 1332 eqindsize = 0 1333 DO i=1, BMesh1 % NumberOfBulkElements 1334 IF(EdgeBasis) THEN 1335 Element => BMesh1 % Parent % Faces(BMesh1 % Elements(i) % ElementIndex) 1336 n = GetElementDOFs(Indexes,Element) 1337 np = GetElementNOFNodes(Element) 1338 ELSE 1339 Element => BMesh1 % Elements(i) 1340 n = Element % TYPE % NumberOfNodes 1341 np = n 1342 Indexes(1:n) = Element % NodeIndexes 1343 END IF 1344 eqindsize = MAX( eqindsize, MAXVAL(Indexes(1:n)) ) 1345 END DO 1346 1347 1348 ! Create the nodal coordinates for all Gaussian integration points 1349 !----------------------------------------------------------------- 1350 NoGaussPoints = 0 1351 DO i=1, BMesh1 % NumberOfBulkElements 1352 Element => BMesh1 % Elements(i) 1353 n = Element % TYPE % NumberOfNodes 1354 NodeIndexes => Element % NodeIndexes 1355 ElementNodes % x(1:n) = RealNodes % x(NodeIndexes(1:n)) 1356 ElementNodes % y(1:n) = RealNodes % y(NodeIndexes(1:n)) 1357 ElementNodes % z(1:n) = RealNodes % z(NodeIndexes(1:n)) 1358 1359 IntegStuff = GaussPoints( Element, RelOrder=RelOrder, EdgeBasis=PiolaT ) 1360 DO j=1,IntegStuff % n 1361 NoGaussPoints = NoGaussPoints + 1 1362 u = IntegStuff % u(j) 1363 v = IntegStuff % v(j) 1364 w = IntegStuff % w(j) 1365 IF(PiolaT) THEN 1366 stat = ElementInfo( Element, ElementNodes, u,v,w, & 1367 detJ, Basis, EdgeBasis=WBasis ) 1368 ELSE 1369 Stat = ElementInfo( Element, ElementNodes, u, v, w, detJ, Basis ) 1370 END IF 1371 GaussNodes % x(NoGaussPoints) = SUM( Basis(1:n) * ElementNodes % x(1:n) ) 1372 GaussNodes % y(NoGaussPoints) = SUM( Basis(1:n) * ElementNodes % y(1:n) ) 1373 GaussNodes % z(NoGaussPoints) = SUM( Basis(1:n) * ElementNodes % z(1:n) ) 1374 END DO 1375 END DO 1376 1377 BMesh1 % Nodes => GaussNodes 1378 BMesh1 % NumberOfNodes = NoGaussPoints 1379 1380 ! Create the mirror node flag and map the nodes of Mesh1 to be 1381 ! in the interval of Mesh2. 1382 !----------------------------------------------------------------- 1383 IF( Repeating ) THEN 1384 IF( AntiRepeating ) THEN 1385 ALLOCATE( MirrorNode( BMesh1 % NumberOfNodes ) ) 1386 MirrorNode = .FALSE. 1387 END IF 1388 CALL PreRotationalProjector(BMesh1, BMesh2, MirrorNode ) 1389 END IF 1390 1391 ! Create the projector for Gaussian integration points 1392 !----------------------------------------------------------------- 1393 GaussProjector => MeshProjector( BMesh2, BMesh1, UseQuadrantTree ) 1394 Rows => GaussProjector % Rows 1395 Cols => GaussProjector % Cols 1396 Values => GaussProjector % Values 1397 1398 ! If there are mirror nodes change the sign 1399 !----------------------------------------------------------------------------- 1400 IF( AntiRepeating ) THEN 1401 CALL PostRotationalProjector( GaussProjector, MirrorNode ) 1402 IF( ALLOCATED( MirrorNode) ) DEALLOCATE( MirrorNode ) 1403 END IF 1404 1405 ! Transfer the projector on the Gaussian points to that on 1406 ! nodal points utilizing the flexibility of the list matrix. 1407 !----------------------------------------------------------- 1408 Projector => AllocateMatrix() 1409 Projector % FORMAT = MATRIX_LIST 1410 Projector % ProjectorType = PROJECTOR_TYPE_GALERKIN 1411 1412 ALLOCATE(Eqind(eqindsize)) 1413 EqInd = 0 1414 1415 ALLOCATE(Projector % InvPerm(eqindsize)) 1416 Projector % InvPerm = 0 1417 1418 Ind = 0 1419 1420 NoGaussPoints = 0 1421 DO i=1, BMesh1 % NumberOfBulkElements 1422 IF(EdgeBasis) THEN 1423 Element => BMesh1 % Parent % Faces(BMesh1 % Elements(i) % ElementIndex) 1424 n = GetElementDOFs(Indexes,Element) 1425 np = GetElementNOFNodes(Element) 1426 IF(ANY(Indexes(1:np)>Element % NodeIndexes)) np=0 1427 CALL GetElementNodes(Nodes,Element) 1428 ELSE 1429 Element => BMesh1 % Elements(i) 1430 n = Element % TYPE % NumberOfNodes 1431 np = n 1432 Indexes(1:n) = Element % NodeIndexes 1433 ElementNodes % x(1:n) = RealNodes % x(Indexes(1:n)) 1434 ElementNodes % y(1:n) = RealNodes % y(Indexes(1:n)) 1435 ElementNodes % z(1:n) = RealNodes % z(Indexes(1:n)) 1436 END IF 1437 1438 IntegStuff = GaussPoints(Element, RelOrder=RelOrder, EdgeBasis=PiolaT ) 1439 DO j=1,IntegStuff % n 1440 NoGaussPoints = NoGaussPoints + 1 1441 u = IntegStuff % u(j) 1442 v = IntegStuff % v(j) 1443 w = IntegStuff % w(j) 1444 1445 IF(EdgeBasis) THEN 1446 IF(PiolaT) THEN 1447 stat = ElementInfo( Element,Nodes,u,v,w,detJ,Basis,EdgeBasis=WBasis) 1448 ELSE 1449 Stat = ElementInfo(Element, Nodes, u, v, w, detJ, Basis,dBasisdx) 1450 CALL GetEdgeBasis(Element,WBasis,RotWBasis,Basis,dBasisdx) 1451 END IF 1452 ELSE 1453 Stat = ElementInfo(Element, ElementNodes, u, v, w, detJ, Basis) 1454 END IF 1455 1456 ! Modify weight so that the projector is consistent with the coordinate system. 1457 weight = detJ * IntegStuff % s(j) 1458 IF( AxisSym ) THEN 1459 IF( EdgeBasis ) THEN 1460 x = SUM( Basis(1:np) * Nodes % x(1:np) ) 1461 ELSE 1462 x = SUM( Basis(1:np) * ElementNodes % x(1:np) ) 1463 END IF 1464 weight = weight * x 1465 END IF 1466 1467 1468 ! Do the numbering of new dofs 1469 ! This needs to be done here because the nodal jump 1470 ! needs the index related to (p,q) pair. 1471 DO p=1,np 1472 IF (EQind(Indexes(p))==0) THEN 1473 Ind = Ind+1 1474 EQind(Indexes(p)) = Ind 1475 IF( EdgeBasis ) THEN 1476 Projector % InvPerm(Ind) = Indexes(p) 1477 ELSE 1478 Projector % InvPerm(Ind) = InvPerm1(Indexes(p)) 1479 END IF 1480 END IF 1481 END DO 1482 1483 DO p=1,np 1484 val = weight * Basis(p) 1485 1486 DO q=1,np 1487 qq = Indexes(q) 1488 IF(.NOT.EdgeBasis) qq=InvPerm1(qq) 1489 CALL List_AddToMatrixElement(Projector % ListMatrix, EQind(Indexes(p)), qq, Basis(q) * val ) 1490 1491 ! Add a diagonal entry to the future constrained system. 1492 ! This will enable a jump to the discontinuous boundary. 1493 ! So far no value is added just the sparse matrix entry. 1494 !IF( NodalJump ) THEN 1495 ! IF( Indexes(p) <= nodesize .AND. Indexes(q) <= nodesize ) THEN 1496 ! CALL List_AddToMatrixElement(Projector % ListMatrix, EQind(Indexes(p)),& 1497 ! totsize + EQInd(Indexes(q)), 0.0_dp ) 1498 ! END IF 1499 !END IF 1500 END DO 1501 1502 DO q = Rows(NoGaussPoints), Rows(NoGaussPoints+1)-1 1503 qq = Cols(q) 1504 IF (qq<=0) EXIT 1505 IF(.NOT.EdgeBasis) qq=InvPerm2(qq) 1506 CALL List_AddToMatrixElement(Projector % ListMatrix, & 1507 EQind(Indexes(p)), qq, -PeriodicScale * Values(q) * val ) 1508 END DO 1509 END DO 1510 1511 IF(EdgeBasis)THEN 1512 DO p=np+1,n 1513 pp=p-np 1514 1515 wq = WBasis(pp,:) 1516 1517 IF (EQind(Indexes(p))==0) THEN 1518 Ind = Ind+1 1519 EQind(Indexes(p)) = Ind 1520 Projector % InvPerm(Ind) = Indexes(p) 1521 END IF 1522 1523 DO q=np+1,n 1524 qq = q-np 1525 1526 vq = Wbasis(qq,:) 1527 CALL List_AddToMatrixElement(Projector % ListMatrix, EQind(Indexes(p)),& 1528 Indexes(q), weight * SUM(vq*wq)) 1529 END DO 1530 1531 1532 DO q = Rows(NoGaussPoints)+np, Rows(NoGaussPoints+1)-1,3 1533 IF(Cols(q)>=0) STOP 'q' 1534 1535 vq(1) = Values(q) 1536 vq(2) = Values(q+1) 1537 vq(3) = Values(q+2) 1538 1539 CALL List_AddToMatrixElement( Projector % ListMatrix,& 1540 EQind(Indexes(p)), -Cols(q), -PeriodicScale * weight * SUM(vq*wq)) 1541 END DO 1542 END DO 1543 ENDIF 1544 END DO 1545 END DO 1546 1547 1548 CALL List_ToCRSMatrix(Projector) 1549 1550 BMesh1 % Nodes => RealNodes 1551 BMesh1 % NumberOfNodes = NoNodes 1552 1553 ! Finally, deallocate permanent storage 1554 !---------------------------------------------------------------- 1555 DEALLOCATE( ElementNodes % x, ElementNodes % y, ElementNodes % z ) 1556 DEALLOCATE( ElementNodes ) 1557 1558 DEALLOCATE( GaussNodes % x, GaussNodes % y, GaussNodes % z) 1559 DEALLOCATE( GaussNodes ) 1560 1561 DEALLOCATE( Basis ) 1562 IF(EdgeBasis) DEALLOCATE( dBasisdx, WBasis, RotWBasis ) 1563 1564!------------------------------------------------------------------------------ 1565 END FUNCTION WeightedProjector 1566!------------------------------------------------------------------------------ 1567