1!/*****************************************************************************/ 2! * 3! * Elmer, A Finite Element Software for Multiphysical Problems 4! * 5! * Copyright 1st April 1995 - , CSC - IT Center for Science Ltd., Finland 6! * 7! * This library is free software; you can redistribute it and/or 8! * modify it under the terms of the GNU Lesser General Public 9! * License as published by the Free Software Foundation; either 10! * version 2.1 of the License, or (at your option) any later version. 11! * 12! * This library is distributed in the hope that it will be useful, 13! * but WITHOUT ANY WARRANTY; without even the implied warranty of 14! * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 15! * Lesser General Public License for more details. 16! * 17! * You should have received a copy of the GNU Lesser General Public 18! * License along with this library (in file ../LGPL-2.1); if not, write 19! * to the Free Software Foundation, Inc., 51 Franklin Street, 20! * Fifth Floor, Boston, MA 02110-1301 USA 21! * 22! *****************************************************************************/ 23! 24!/****************************************************************************** 25! * 26! * Authors: Juha Ruokolainen 27! * Email: Juha.Ruokolainen@csc.fi 28! * Web: http://www.csc.fi/elmer 29! * Address: CSC - IT Center for Science Ltd. 30! * Keilaranta 14 31! * 02101 Espoo, Finland 32! * 33! * Original Date: 12 Jun 2003 34! * 35! ******************************************************************************/ 36 37!> \defgroup ElmerLib Elmer library 38!> \{ 39!> \defgroup DefUtils Default API 40!> \{ 41 42!-------------------------------------------------------------------------------- 43!> Module containing utility subroutines with default values for various 44!> system subroutine arguments. For user defined solvers and subroutines 45!> this module should provide all the needed functionality for typical finite 46!> element procedures. 47!-------------------------------------------------------------------------------- 48MODULE DefUtils 49 50#include "../config.h" 51 52 USE Adaptive 53 USE SolverUtils 54 55 IMPLICIT NONE 56 57 INTERFACE DefaultUpdateEquations 58 MODULE PROCEDURE DefaultUpdateEquationsR, DefaultUpdateEquationsC 59 END INTERFACE 60 61 INTERFACE DefaultUpdatePrec 62 MODULE PROCEDURE DefaultUpdatePrecR, DefaultUpdatePrecC 63 END INTERFACE 64 65 INTERFACE DefaultUpdateMass 66 MODULE PROCEDURE DefaultUpdateMassR, DefaultUpdateMassC 67 END INTERFACE 68 69 INTERFACE DefaultUpdateBulk 70 MODULE PROCEDURE DefaultUpdateBulkR, DefaultUpdateBulkC 71 END INTERFACE 72 73 INTERFACE DefaultUpdateDamp 74 MODULE PROCEDURE DefaultUpdateDampR, DefaultUpdateDampC 75 END INTERFACE 76 77 INTERFACE DefaultUpdateForce 78 MODULE PROCEDURE DefaultUpdateForceR, DefaultUpdateForceC 79 END INTERFACE 80 81 INTERFACE DefaultUpdateTimeForce 82 MODULE PROCEDURE DefaultUpdateTimeForceR, DefaultUpdateTimeForceC 83 END INTERFACE 84 85 INTERFACE Default1stOrderTime 86 MODULE PROCEDURE Default1stOrderTimeR, Default1stOrderTimeC 87 END INTERFACE 88 89 INTERFACE Default2ndOrderTime 90 MODULE PROCEDURE Default2ndOrderTimeR, Default2ndOrderTimeC 91 END INTERFACE 92 93 INTERFACE GetLocalSolution 94 MODULE PROCEDURE GetScalarLocalSolution, GetVectorLocalSolution 95 END INTERFACE 96 97 INTERFACE GetLocalEigenmode 98 MODULE PROCEDURE GetScalarLocalEigenmode, GetVectorLocalEigenmode 99 END INTERFACE 100 101 INTEGER, ALLOCATABLE, TARGET, PRIVATE :: IndexStore(:), VecIndexStore(:) 102 REAL(KIND=dp), ALLOCATABLE, TARGET, PRIVATE :: ValueStore(:) 103 !$OMP THREADPRIVATE(IndexStore, VecIndexStore, ValueStore) 104 105 TYPE(Element_t), POINTER :: CurrentElementThread => NULL() 106 !$OMP THREADPRIVATE(CurrentElementThread) 107 108 ! TODO: Get actual values for these from mesh 109 INTEGER, PARAMETER, PRIVATE :: ISTORE_MAX_SIZE = 1024 110 INTEGER, PARAMETER, PRIVATE :: VSTORE_MAX_SIZE = 1024 111 PRIVATE :: GetIndexStore, GetVecIndexStore, GetValueStore 112CONTAINS 113 114 115 FUNCTION GetVersion() RESULT(ch) 116 CHARACTER(LEN=:), ALLOCATABLE :: ch 117 ch = VERSION 118 END FUNCTION GetVersion 119 120 FUNCTION GetSifName(Found) RESULT(ch) 121 CHARACTER(LEN=:), ALLOCATABLE :: ch 122 LOGICAL, OPTIONAL :: Found 123 ch = ListGetString( CurrentModel % Simulation,'Solver Input File', Found ) 124 END FUNCTION GetSifName 125 126 FUNCTION GetRevision(Found) RESULT(ch) 127 CHARACTER(LEN=:), ALLOCATABLE :: ch 128 LOGICAL, OPTIONAL :: Found 129#ifdef REVISION 130 ch = REVISION 131 IF(PRESENT(Found)) Found = .TRUE. 132#else 133 ch = "unknown" 134 IF(PRESENT(Found)) Found = .FALSE. 135#endif 136 END FUNCTION GetRevision 137 138 FUNCTION GetCompilationDate(Found) RESULT(ch) 139 CHARACTER(LEN=:), ALLOCATABLE :: ch 140 LOGICAL, OPTIONAL :: Found 141#ifdef COMPILATIONDATE 142 ch = COMPILATIONDATE 143 IF(PRESENT(Found)) Found = .TRUE. 144#else 145 ch = "unknown" 146 IF(PRESENT(Found)) Found = .FALSE. 147#endif 148 END FUNCTION GetCompilationDate 149 150 FUNCTION GetIndexStore() RESULT(ind) 151 IMPLICIT NONE 152 INTEGER, POINTER CONTIG :: ind(:) 153 INTEGER :: istat 154 155 IF ( .NOT. ALLOCATED(IndexStore) ) THEN 156 ALLOCATE( IndexStore(ISTORE_MAX_SIZE), STAT=istat ) 157 IndexStore = 0 158 IF ( Istat /= 0 ) CALL Fatal( 'GetIndexStore', & 159 'Memory allocation error.' ) 160 END IF 161 ind => IndexStore 162 END FUNCTION GetIndexStore 163 164 FUNCTION GetVecIndexStore() RESULT(ind) 165 IMPLICIT NONE 166 INTEGER, POINTER CONTIG :: ind(:) 167 INTEGER :: istat 168 169 IF ( .NOT. ALLOCATED(VecIndexStore) ) THEN 170 ALLOCATE( VecIndexStore(ISTORE_MAX_SIZE), STAT=istat ) 171 VecIndexStore = 0 172 IF ( istat /= 0 ) CALL Fatal( 'GetVecIndexStore', & 173 'Memory allocation error.' ) 174 END IF 175 ind => VecIndexStore 176 END FUNCTION GetVecIndexStore 177 178 FUNCTION GetValueStore(n) RESULT(val) 179 IMPLICIT NONE 180 REAL(KIND=dp), POINTER CONTIG :: val(:) 181 INTEGER :: n, istat 182 183 IF ( .NOT.ALLOCATED(ValueStore) ) THEN 184 ALLOCATE( ValueStore(VSTORE_MAX_SIZE), STAT=istat ) 185 ValueStore = REAL(0, dp) 186 IF ( Istat /= 0 ) CALL Fatal( 'GetValueStore', & 187 'Memory allocation error.' ) 188 END IF 189 IF (n > VSTORE_MAX_SIZE) THEN 190 CALL Fatal( 'GetValueStore', 'Not enough memory allocated for store.' ) 191 END IF 192 val => ValueStore(1:n) 193 END FUNCTION GetValueStore 194 195!> Returns handle to the active solver 196 FUNCTION GetSolver() RESULT( Solver ) 197 TYPE(Solver_t), POINTER :: Solver 198 Solver => CurrentModel % Solver 199 END FUNCTION GetSolver 200 201!> Returns handle to the active matrix 202 FUNCTION GetMatrix( USolver ) RESULT( Matrix ) 203 TYPE(Matrix_t), POINTER :: Matrix 204 TYPE(Solver_t), OPTIONAL, TARGET :: USolver 205 206 IF ( PRESENT( USolver ) ) THEN 207 Matrix => USolver % Matrix 208 ELSE 209 Matrix => CurrentModel % Solver % Matrix 210 END IF 211 END FUNCTION GetMatrix 212 213!> Returns handle to the active mesh 214 FUNCTION GetMesh( USolver ) RESULT( Mesh ) 215 TYPE(Mesh_t), POINTER :: Mesh 216 TYPE(Solver_t), OPTIONAL, TARGET :: USolver 217 218 IF ( PRESENT( USolver ) ) THEN 219 Mesh => USolver % MEsh 220 ELSE 221 Mesh => CurrentModel % Solver % Mesh 222 END IF 223 END FUNCTION GetMesh 224 225!> Returns handle to the active element 226 FUNCTION GetCurrentElement(Element) RESULT(Ret_Element) 227 IMPLICIT NONE 228 TYPE(Element_t), OPTIONAL, TARGET :: Element 229 TYPE(Element_t), POINTER :: Ret_Element 230 231 IF (PRESENT(Element)) THEN 232 Ret_Element=>Element 233 ELSE 234#ifdef _OPENMP 235 IF (omp_in_parallel()) THEN 236 Ret_Element=>CurrentElementThread 237 ELSE 238 Ret_Element=>CurrentModel % CurrentElement 239 END IF 240#else 241 Ret_Element => CurrentModel % CurrentElement 242#endif 243 END IF 244 END FUNCTION GetCurrentElement 245 246!> Sets handle to the active element of the current thread. 247!> Old handle is given as a return value as what would be returned 248!> by a call to GetCurrentElement 249 FUNCTION SetCurrentElement(Element) RESULT(OldElement) 250 IMPLICIT NONE 251 TYPE(Element_t), TARGET :: Element 252 TYPE(Element_t), POINTER :: OldElement 253 254#ifdef _OPENMP 255 IF (omp_in_parallel()) THEN 256 OldElement => CurrentElementThread 257 CurrentElementThread => Element 258 ELSE 259 OldElement => CurrentModel % CurrentElement 260 CurrentModel % CurrentElement => Element 261 END IF 262#else 263 OldElement => CurrentModel % CurrentElement 264 CurrentModel % CurrentElement => Element 265#endif 266 END FUNCTION SetCurrentElement 267 268!> Returns handle to the index of the current element 269 FUNCTION GetElementIndex(Element) RESULT(Indx) 270 TYPE(Element_t), OPTIONAL :: Element 271 INTEGER :: Indx 272 TYPE(Element_t), POINTER :: CurrElement 273 274 CurrElement => GetCurrentElement(Element) 275 Indx = CurrElement % ElementIndex 276 END FUNCTION GetElementIndex 277 278 SUBROUTINE GetElementNodeIndex(i, Element, n, FOUND) 279 IMPLICIT None 280 281 ! variables in function header 282 INTEGER :: i, n 283 TYPE(Element_t), POINTER :: Element 284 Logical :: FOUND 285 286 DO i=1, SIZE(Element%NodeIndexes) 287 IF (n == Element%NodeIndexes(i)) THEN 288 FOUND=.TRUE. 289 EXIT 290 END IF 291 END DO 292 END SUBROUTINE GetElementNodeIndex 293 294 FUNCTION GetIPIndex( LocalIp, USolver, Element, IpVar ) RESULT ( GlobalIp ) 295 INTEGER :: LocalIp, GlobalIp 296 297 TYPE(Solver_t), OPTIONAL, TARGET :: USolver 298 TYPE(Element_t), OPTIONAL :: Element 299 TYPE(Variable_t), POINTER, OPTIONAL :: IpVar 300 301 TYPE(Solver_t), POINTER :: Solver 302 TYPE(Element_t), POINTER :: CurrElement 303 INTEGER :: n, m 304 INTEGER, POINTER :: IpPerm(:) 305 306 IF ( PRESENT( USolver ) ) THEN 307 Solver => USolver 308 ELSE 309 Solver => CurrentModel % Solver 310 END IF 311 312 CurrElement => GetCurrentElement(Element) 313 n = CurrElement % ElementIndex 314 GlobalIp = 0 315 316 IF( PRESENT( IpVar ) ) THEN 317 IF( IpVar % TYPE /= Variable_on_gauss_points ) THEN 318 CALL Fatal('GetIpIndex','Variable is not of type gauss points!') 319 END IF 320 321 IpPerm => IpVar % Perm 322 m = IpPerm(n+1) - IpPerm(n) 323 324 ! This is a sign that the variable is not active at the element 325 IF( m == 0 ) RETURN 326 ELSE 327 IF( .NOT. ASSOCIATED( Solver % IpTable ) ) THEN 328 CALL Fatal('GetIpIndex','Cannot access index of gaussian point!') 329 END IF 330 331 IpPerm => Solver % IpTable % IpOffset 332 m = IpPerm(n+1) - IpPerm(n) 333 END IF 334 335 ! There are not sufficient number of gauss points in the permutation table to have a 336 ! local index this big. 337 IF( m < LocalIp ) THEN 338 CALL Warn('GetIpIndex','Inconsistent number of IP points!') 339 RETURN 340 END IF 341 342 GlobalIp = IpPerm(n) + LocalIp 343 344 END FUNCTION GetIPIndex 345 346 347 348 FUNCTION GetIPCount( USolver, IpVar ) RESULT ( IpCount ) 349 INTEGER :: IpCount 350 TYPE(Solver_t), OPTIONAL, TARGET :: USolver 351 TYPE(Variable_t), OPTIONAL, POINTER :: IpVar 352 353 TYPE(Solver_t), POINTER :: Solver 354 INTEGER, POINTER :: IpPerm 355 356 IF ( PRESENT( USolver ) ) THEN 357 Solver => USolver 358 ELSE 359 Solver => CurrentModel % Solver 360 END IF 361 362 IF( PRESENT( IpVar ) ) THEN 363 IF( IpVar % TYPE /= Variable_on_gauss_points ) THEN 364 CALL Fatal('GetIpCount','Variable is not of type gauss points!') 365 END IF 366 IpCount = SIZE( IpVar % Values ) / IpVar % Dofs 367 ELSE 368 IF( .NOT. ASSOCIATED( Solver % IpTable ) ) THEN 369 CALL Fatal('GetIpCount','Gauss point table not initialized') 370 END IF 371 IpCount = Solver % IpTable % IpCount 372 END IF 373 374 END FUNCTION GetIPCount 375 376 377!> Returns the number of active elements for the current solver 378 FUNCTION GetNOFActive( USolver ) RESULT(n) 379 INTEGER :: n 380 TYPE(Solver_t), OPTIONAL, TARGET :: USolver 381 TYPE(Solver_t), POINTER :: Solver 382 383 IF ( PRESENT( USolver ) ) THEN 384 Solver => USolver 385 ELSE 386 Solver => CurrentModel % Solver 387 END IF 388 389 IF( ASSOCIATED( Solver % ColourIndexList ) ) THEN 390 Solver % CurrentColour = Solver % CurrentColour + 1 391 n = Solver % ColourIndexList % ptr(Solver % CurrentColour+1) & 392 - Solver % ColourIndexList % ptr(Solver % CurrentColour) 393 CALL Info('GetNOFActive','Number of active elements: '& 394 //TRIM(I2S(n))//' in colour '//TRIM(I2S(Solver % CurrentColour)),Level=20) 395 ELSE 396 n = Solver % NumberOfActiveElements 397 CALL Info('GetNOFActive','Number of active elements: '& 398 //TRIM(I2S(n)),Level=20) 399 END IF 400 401 END FUNCTION GetNOFActive 402 403 !> Return number of boundary elements of the current boundary colour 404 !> and increments the colour counter 405 FUNCTION GetNOFBoundaryActive( USolver ) RESULT(n) 406 INTEGER :: n 407 TYPE(Solver_t), OPTIONAL, TARGET :: USolver 408 TYPE(Solver_t), POINTER :: Solver 409 410 IF ( PRESENT( USolver ) ) THEN 411 Solver => USolver 412 ELSE 413 Solver => CurrentModel % Solver 414 END IF 415 416 IF( ASSOCIATED( Solver % BoundaryColourIndexList ) ) THEN 417 Solver % CurrentBoundaryColour = Solver % CurrentBoundaryColour + 1 418 n = Solver % BoundaryColourIndexList % ptr(Solver % CurrentBoundaryColour+1) & 419 - Solver % BoundaryColourIndexList % ptr(Solver % CurrentBoundaryColour) 420 CALL Info('GetNOFBoundaryActive','Number of boundary elements: '& 421 //TRIM(I2S(n))//' in colour '//TRIM(I2S(Solver % CurrentBoundaryColour)),Level=20) 422 ELSE 423 n = Solver % Mesh % NumberOfBoundaryElements 424 CALL Info('GetNOFBoundaryActive','Number of active elements: '& 425 //TRIM(I2S(n)),Level=20) 426 END IF 427 428 END FUNCTION GetNOFBoundaryActive 429 430!> Returns the current time 431 FUNCTION GetTime() RESULT(st) 432 REAL(KIND=dp) :: st 433 TYPE(Variable_t), POINTER :: v 434 435 v => CurrentModel % Solver % Mesh % Variables 436 v => VariableGet( v, 'time' ) 437 st = v % Values(1) 438 END FUNCTION GetTime 439 440!> Returns the current periodic time 441 FUNCTION GetPeriodicTime() RESULT(st) 442 REAL(KIND=dp) :: st 443 TYPE(Variable_t), POINTER :: v 444 445 v => CurrentModel % Solver % Mesh % Variables 446 v => VariableGet( v, 'periodic time' ) 447 st = v % Values(1) 448 END FUNCTION GetPeriodicTime 449 450!> Returns the current timestep 451 FUNCTION GetTimeStep() RESULT(st) 452 INTEGER :: st 453 TYPE(Variable_t), POINTER :: v 454 455 v => CurrentModel % Solver % Mesh % Variables 456 v => VariableGet( v, 'timestep' ) 457 st = NINT(v % Values(1)) 458 END FUNCTION GetTimestep 459 460!> Returns the current timestep interval 461 FUNCTION GetTimeStepInterval() RESULT(st) 462 INTEGER :: st 463 TYPE(Variable_t), POINTER :: v 464 465 v => CurrentModel % Solver % Mesh % Variables 466 v => VariableGet( v, 'timestep interval') 467 st = NINT(v % Values(1)) 468 END FUNCTION GetTimestepInterval 469 470!> Returns the current timestep size 471 FUNCTION GetTimestepSize() RESULT(st) 472 REAL(KIND=dp) :: st 473 TYPE(Variable_t), POINTER :: v 474 475 v => CurrentModel % Solver % Mesh % Variables 476 v => VariableGet( v, 'timestep size') 477 st = v % Values(1) 478 END FUNCTION GetTimestepSize 479 480!> Returns the angular frequency 481 FUNCTION GetAngularFrequency(ValueList,Found, UElement) RESULT(w) 482 REAL(KIND=dp) :: w 483 TYPE(ValueList_t), POINTER, OPTIONAL :: ValueList 484 LOGICAL, OPTIONAL :: Found 485 TYPE(Element_t), POINTER, OPTIONAL :: UElement 486 487 w = ListGetAngularFrequency( ValueList, Found, UElement ) 488 END FUNCTION GetAngularFrequency 489 490!> Returns the current coupled system iteration loop count 491 FUNCTION GetCoupledIter() RESULT(st) 492 INTEGER :: st 493 TYPE(Variable_t), POINTER :: v 494 495 v => CurrentModel % Solver % Mesh % Variables 496 v => VariableGet( v, 'coupled iter') 497 st = NINT(v % Values(1)) 498 END FUNCTION GetCoupledIter 499 500!> Returns the current nonlinear system iteration loop count 501 FUNCTION GetNonlinIter() RESULT(st) 502 INTEGER :: st 503 TYPE(Variable_t), POINTER :: v 504 505 v => CurrentModel % Solver % Mesh % Variables 506 v => VariableGet( v, 'nonlin iter') 507 st = NINT(v % Values(1)) 508 END FUNCTION GetNonlinIter 509 510!> Returns the number of boundary elements 511 FUNCTION GetNOFBoundaryElements( UMesh ) RESULT(n) 512 INTEGER :: n 513 TYPE(Mesh_t), OPTIONAL :: UMesh 514 515 IF ( PRESENT( UMesh ) ) THEN 516 n = UMesh % NumberOfBoundaryElements 517 ELSE 518 n = CurrentModel % Mesh % NumberOfBoundaryElements 519 END IF 520 END FUNCTION GetNOFBoundaryElements 521 522!> Returns a scalar field in the nodes of the element 523 SUBROUTINE GetScalarLocalSolution( x,name,UElement,USolver,tStep, UVariable ) 524 REAL(KIND=dp) :: x(:) 525 CHARACTER(LEN=*), OPTIONAL :: name 526 TYPE(Solver_t) , OPTIONAL, TARGET :: USolver 527 TYPE(Element_t), OPTIONAL, TARGET :: UElement 528 TYPE(Variable_t), OPTIONAL, TARGET :: UVariable 529 INTEGER, OPTIONAL :: tStep 530 531 REAL(KIND=dp), POINTER :: Values(:) 532 TYPE(Variable_t), POINTER :: Variable 533 TYPE(Solver_t) , POINTER :: Solver 534 TYPE(Element_t), POINTER :: Element 535 536 INTEGER :: i, j, k, n 537 INTEGER, POINTER :: Indexes(:) 538 539 Solver => CurrentModel % Solver 540 IF ( PRESENT(USolver) ) Solver => USolver 541 542 x = 0.0d0 543 544 IF(.NOT. PRESENT(UVariable)) THEN 545 Variable => Solver % Variable 546 ELSE 547 Variable => UVariable 548 END IF 549 550 IF ( PRESENT(name) ) THEN 551 Variable => VariableGet( Solver % Mesh % Variables, name ) 552 END IF 553 IF ( .NOT. ASSOCIATED( Variable ) ) RETURN 554 555 Element => GetCurrentElement(UElement) 556 557 Values => Variable % Values 558 IF ( PRESENT(tStep) ) THEN 559 IF ( tStep<0 ) THEN 560 IF ( ASSOCIATED(Variable % PrevValues) ) THEN 561 IF( -tStep<=SIZE(Variable % PrevValues,2)) & 562 Values => Variable % PrevValues(:,-tStep) 563 END IF 564 END IF 565 END IF 566 567 ! If variable is defined on gauss points return that instead 568 IF( Variable % TYPE == Variable_on_gauss_points ) THEN 569 j = Element % ElementIndex 570 n = Variable % Perm(j+1) - Variable % Perm(j) 571 DO i=1,n 572 x(i) = Values(Variable % Perm(j) + i) 573 END DO 574 RETURN 575 END IF 576 577 578 Indexes => GetIndexStore() 579 IF ( ASSOCIATED(Variable % Solver) ) THEN 580 n = GetElementDOFs( Indexes, Element, Variable % Solver ) 581 ELSE 582 n = GetElementDOFs( Indexes, Element, Solver ) 583 END IF 584 n = MIN( n, SIZE(x) ) 585 586 587 IF ( ASSOCIATED( Variable % Perm ) ) THEN 588 IF( Variable % PeriodicFlipActive ) THEN 589 DO i=1,n 590 j = Indexes(i) 591 IF ( j>0 .AND. j<=SIZE(Variable % Perm) ) THEN 592 k = Variable % Perm(j) 593 IF ( k>0 ) THEN 594 x(i) = Values(k) 595 IF( CurrentModel % Mesh % PeriodicFlip(j) ) x(i) = -x(i) 596 END IF 597 END IF 598 END DO 599 ELSE 600 DO i=1,n 601 j = Indexes(i) 602 IF ( j>0 .AND. j<=SIZE(Variable % Perm) ) THEN 603 j = Variable % Perm(j) 604 IF ( j>0 ) x(i) = Values(j) 605 END IF 606 END DO 607 END IF 608 ELSE 609 DO i=1,n 610 j = Indexes(i) 611 IF ( j>0 .AND. j<=SIZE(Variable % Values) ) & 612 x(i) = Values(Indexes(i)) 613 END DO 614 END IF 615 END SUBROUTINE GetScalarLocalSolution 616 617 618 619!> Returns a vector field in the nodes of the element 620 SUBROUTINE GetVectorLocalSolution( x,name,UElement,USolver,tStep, UVariable ) 621 REAL(KIND=dp) :: x(:,:) 622 CHARACTER(LEN=*), OPTIONAL :: name 623 TYPE(Solver_t), OPTIONAL, TARGET :: USolver 624 TYPE(Element_t), OPTIONAL, TARGET :: UElement 625 TYPE(Variable_t), OPTIONAL, TARGET :: UVariable 626 INTEGER, OPTIONAL :: tStep 627 628 TYPE(Variable_t), POINTER :: Variable 629 TYPE(Solver_t) , POINTER :: Solver 630 TYPE(Element_t), POINTER :: Element 631 632 INTEGER :: i, j, k, l, n 633 INTEGER, POINTER :: Indexes(:) 634 REAL(KIND=dp), POINTER :: Values(:) 635 636 Solver => CurrentModel % Solver 637 IF ( PRESENT(USolver) ) Solver => USolver 638 639 x = 0.0d0 640 641 IF(.NOT. PRESENT(UVariable)) THEN 642 Variable => Solver % Variable 643 ELSE 644 Variable => UVariable 645 END IF 646 647 IF ( PRESENT(name) ) THEN 648 Variable => VariableGet( Solver % Mesh % Variables, name ) 649 END IF 650 IF ( .NOT. ASSOCIATED( Variable ) ) RETURN 651 652 Element => GetCurrentElement(UElement) 653 654 Indexes => GetIndexStore() 655 IF ( ASSOCIATED(Variable % Solver ) ) THEN 656 n = GetElementDOFs( Indexes, Element, Variable % Solver ) 657 ELSE 658 n = GetElementDOFs( Indexes, Element, Solver ) 659 END IF 660 n = MIN( n, SIZE(x,2) ) 661 662 Values => Variable % Values 663 IF ( PRESENT(tStep) ) THEN 664 IF ( tStep<0 ) THEN 665 IF ( ASSOCIATED(Variable % PrevValues) ) THEN 666 IF ( -tStep<=SIZE(Variable % PrevValues,2)) & 667 Values => Variable % PrevValues(:,-tStep) 668 END IF 669 END IF 670 END IF 671 672 673 ! If variable is defined on gauss points return that instead 674 IF( Variable % TYPE == Variable_on_gauss_points ) THEN 675 ASSOCIATE(dofs => variable % dofs) 676 j = Element % ElementIndex 677 n = Variable % Perm(j+1) - Variable % Perm(j) 678 IF (size(x,1) < dofs .or. size(x,2) < n) THEN 679 write (message,*) 'Attempting to get IP solution to a too small array of size', & 680 shape(x), '. Required size:', dofs, n 681 CALL Fatal('GetVectorLocalSolution', message) 682 END IF 683 DO i=1,n 684 ASSOCIATE(p => variable % perm(j) + i) 685 DO k=1,dofs 686 x(k, i) = Values((p-1)*dofs + k) 687 END DO 688 END ASSOCIATE 689 END DO 690 RETURN 691 END ASSOCIATE 692 END IF 693 694 DO i=1,Variable % DOFs 695 IF ( ASSOCIATED( Variable % Perm ) ) THEN 696 IF( Variable % PeriodicFlipActive ) THEN 697 DO j=1,n 698 k = Indexes(j) 699 IF ( k>0 .AND. k<=SIZE(Variable % Perm) ) THEN 700 l = Variable % Perm(k) 701 IF( l>0 ) THEN 702 x(i,j) = Values(Variable % DOFs*(l-1)+i) 703 IF( CurrentModel % Mesh % PeriodicFlip(k) ) x(i,j) = -x(i,j) 704 END IF 705 END IF 706 END DO 707 ELSE 708 DO j=1,n 709 k = Indexes(j) 710 IF ( k>0 .AND. k<=SIZE(Variable % Perm) ) THEN 711 l = Variable % Perm(k) 712 IF (l>0) x(i,j) = Values(Variable % DOFs*(l-1)+i) 713 END IF 714 END DO 715 END IF 716 ELSE 717 DO j=1,n 718 IF ( Variable % DOFs*(Indexes(j)-1)+i <= & 719 SIZE( Variable % Values ) ) THEN 720 x(i,j) = Values(Variable % DOFs*(Indexes(j)-1)+i) 721 END IF 722 END DO 723 END IF 724 END DO 725 END SUBROUTINE GetVectorLocalSolution 726 727 728! Eigenmodes may be used as a basis of sensitivity analysis, model reduction etc. 729! then these subroutines may be used to obtain the local eigenmodes 730!------------------------------------------------------------------------------- 731 732!> Returns the number of eigenmodes 733 FUNCTION GetNofEigenModes( name,USolver) RESULT (NofEigenModes) 734 735 CHARACTER(LEN=*), OPTIONAL :: name 736 TYPE(Solver_t) , OPTIONAL, TARGET :: USolver 737 INTEGER :: NofEigenModes 738 739 REAL(KIND=dp), POINTER :: Values(:) 740 TYPE(Variable_t), POINTER :: Variable 741 TYPE(Solver_t) , POINTER :: Solver 742 743 NofEigenModes = 0 744 745 Solver => CurrentModel % Solver 746 IF ( PRESENT(USolver) ) Solver => USolver 747 748 Variable => Solver % Variable 749 IF ( PRESENT(name) ) THEN 750 Variable => VariableGet( Solver % Mesh % Variables, name ) 751 END IF 752 753 IF ( .NOT. ASSOCIATED( Variable ) ) RETURN 754 IF ( .NOT. ASSOCIATED( Variable % EigenValues ) ) RETURN 755 756 NofEigenModes = SIZE( Variable % EigenValues, 1) 757 END FUNCTION GetNofEigenModes 758 759 760!> Returns the desired eigenmode as a scalar field in an element 761 SUBROUTINE GetScalarLocalEigenmode( x,name,UElement,USolver,NoEigen,ComplexPart ) 762 REAL(KIND=dp) :: x(:) 763 CHARACTER(LEN=*), OPTIONAL :: name 764 TYPE(Solver_t) , OPTIONAL, TARGET :: USolver 765 TYPE(Element_t), OPTIONAL, TARGET :: UElement 766 INTEGER, OPTIONAL :: NoEigen 767 LOGICAL, OPTIONAL :: ComplexPart 768 769 COMPLEX(KIND=dp), POINTER :: Values(:) 770 TYPE(Variable_t), POINTER :: Variable 771 TYPE(Solver_t) , POINTER :: Solver 772 TYPE(Element_t), POINTER :: Element 773 LOGICAL :: IsComplex 774 775 INTEGER :: i, j, n 776 INTEGER, POINTER :: Indexes(:) 777 778 Solver => CurrentModel % Solver 779 IF ( PRESENT(USolver) ) Solver => USolver 780 781 x = 0.0d0 782 783 Variable => Solver % Variable 784 IF ( PRESENT(name) ) THEN 785 Variable => VariableGet( Solver % Mesh % Variables, name ) 786 END IF 787 IF ( .NOT. ASSOCIATED( Variable ) ) RETURN 788 IF ( .NOT. ASSOCIATED( Variable % EigenVectors ) ) RETURN 789 790 IsComplex = .FALSE. 791 IF( PRESENT( ComplexPart) ) IsComplex = ComplexPart 792 793 Element => GetCurrentElement(UElement) 794 795 IF ( ASSOCIATED( Variable ) ) THEN 796 Indexes => GetIndexStore() 797 IF ( ASSOCIATED(Variable % Solver ) ) THEN 798 n = GetElementDOFs( Indexes, Element, Variable % Solver ) 799 ELSE 800 n = GetElementDOFs( Indexes, Element, Solver ) 801 END IF 802 n = MIN( n, SIZE(x) ) 803 804 Values => Variable % EigenVectors( :, NoEigen ) 805 806 IF ( ASSOCIATED( Variable % Perm ) ) THEN 807 DO i=1,n 808 j = Indexes(i) 809 IF ( j>0 .AND. j<= SIZE(Variable % Perm)) THEN 810 j = Variable % Perm(j) 811 IF ( j>0 ) THEN 812 IF ( IsComplex ) THEN 813 x(i) = AIMAG(Values(j)) 814 ELSE 815 x(i) = REAL(Values(j)) 816 END IF 817 END IF 818 END IF 819 END DO 820 ELSE 821 x(1:n) = Values(Indexes(1:n)) 822 END IF 823 END IF 824 END SUBROUTINE GetScalarLocalEigenmode 825 826 827 828!> Returns the desired eigenmode as a vector field in an element 829 SUBROUTINE GetVectorLocalEigenmode( x,name,UElement,USolver,NoEigen,ComplexPart ) 830 REAL(KIND=dp) :: x(:,:) 831 CHARACTER(LEN=*), OPTIONAL :: name 832 TYPE(Solver_t), OPTIONAL, TARGET :: USolver 833 TYPE(Element_t), OPTIONAL, TARGET :: UElement 834 INTEGER, OPTIONAL :: NoEigen 835 LOGICAL, OPTIONAL :: ComplexPart 836 837 TYPE(Variable_t), POINTER :: Variable 838 TYPE(Solver_t) , POINTER :: Solver 839 TYPE(Element_t), POINTER :: Element 840 LOGICAL :: IsComplex 841 842 INTEGER :: i, j, k, n 843 INTEGER, POINTER :: Indexes(:) 844 COMPLEX(KIND=dp), POINTER :: Values(:) 845 846 Solver => CurrentModel % Solver 847 IF ( PRESENT(USolver) ) Solver => USolver 848 849 IsComplex = .FALSE. 850 IF( PRESENT( ComplexPart) ) IsComplex = ComplexPart 851 852 x = 0.0d0 853 854 Variable => Solver % Variable 855 IF ( PRESENT(name) ) THEN 856 Variable => VariableGet( Solver % Mesh % Variables, name ) 857 END IF 858 IF ( .NOT. ASSOCIATED( Variable ) ) RETURN 859 IF ( .NOT. ASSOCIATED( Variable % EigenVectors ) ) RETURN 860 861 Element => GetCurrentElement(UElement) 862 863 IF ( ASSOCIATED( Variable ) ) THEN 864 Indexes => GetIndexStore() 865 IF ( ASSOCIATED(Variable % Solver ) ) THEN 866 n = GetElementDOFs( Indexes, Element, Variable % Solver ) 867 ELSE 868 n = GetElementDOFs( Indexes, Element, Solver ) 869 END IF 870 n = MIN( n, SIZE(x) ) 871 872 Values => Variable % EigenVectors( :, NoEigen ) 873 874 DO i=1,Variable % DOFs 875 IF ( ASSOCIATED( Variable % Perm ) ) THEN 876 DO j=1,n 877 k = Indexes(j) 878 IF ( k>0 .AND. k<= SIZE(Variable % Perm)) THEN 879 k = Variable % Perm(k) 880 IF ( k>0 ) THEN 881 IF ( IsComplex ) THEN 882 x(i,j) = AIMAG(Values(Variable % DOFs*(k-1)+i)) 883 ELSE 884 x(i,j) = REAL(Values(Variable % DOFs*(k-1)+i)) 885 END IF 886 END IF 887 END IF 888 END DO 889 ELSE 890 DO j=1,n 891 IF( IsComplex ) THEN 892 x(i,j) = AIMAG( Values(Variable % DOFs*(Indexes(j)-1)+i) ) 893 ELSE 894 x(i,j) = REAL( Values(Variable % DOFs*(Indexes(j)-1)+i) ) 895 END IF 896 END DO 897 END IF 898 END DO 899 END IF 900 END SUBROUTINE GetVectorLocalEigenmode 901 902 903 FUNCTION DefaultVariableGet( Name, ThisOnly, USolver ) RESULT ( Var ) 904 905 CHARACTER(LEN=*) :: Name 906 LOGICAL, OPTIONAL :: ThisOnly 907 TYPE(Solver_t), POINTER, OPTIONAL :: USolver 908 TYPE(Variable_t), POINTER :: Var 909!------------------------------------------------------------------------------ 910 TYPE(Variable_t), POINTER :: Variables 911 912 IF( PRESENT( USolver ) ) THEN 913 Variables => USolver % Mesh % Variables 914 ELSE 915 Variables => CurrentModel % Solver % Mesh % Variables 916 END IF 917 918 Var => VariableGet( Variables, Name, ThisOnly ) 919 920 END FUNCTION DefaultVariableGet 921 922 923!------------------------------------------------------------------------------ 924!> Add variable to the default variable list. 925!------------------------------------------------------------------------------ 926 SUBROUTINE DefaultVariableAdd( Name, DOFs, Perm, Values,& 927 Output,Secondary,VariableType,Global,InitValue,USolver,Var ) 928 929 CHARACTER(LEN=*) :: Name 930 INTEGER, OPTIONAL :: DOFs 931 REAL(KIND=dp), OPTIONAL, POINTER :: Values(:) 932 LOGICAL, OPTIONAL :: Output 933 INTEGER, OPTIONAL, POINTER :: Perm(:) 934 LOGICAL, OPTIONAL :: Secondary 935 INTEGER, OPTIONAL :: VariableType 936 LOGICAL, OPTIONAL :: Global 937 REAL(KIND=dp), OPTIONAL :: InitValue 938 TYPE(Solver_t), OPTIONAL, TARGET :: USolver 939 TYPE(Variable_t), OPTIONAL, POINTER :: Var 940 !------------------------------------------------------------------------------ 941 TYPE(Variable_t), POINTER :: Variables 942 TYPE(Mesh_t), POINTER :: Mesh 943 TYPE(Solver_t), POINTER :: Solver 944 945 IF( PRESENT( USolver ) ) THEN 946 Solver => USolver 947 ELSE 948 Solver => CurrentModel % Solver 949 END IF 950 Mesh => Solver % Mesh 951 Variables => Mesh % Variables 952 953 CALL VariableAddVector( Variables,Mesh,Solver,Name,DOFs,Values,& 954 Perm,Output,Secondary,VariableType,Global,InitValue ) 955 956 IF( PRESENT( Var ) ) THEN 957 Var => VariableGet( Variables, Name ) 958 END IF 959 960 END SUBROUTINE DefaultVariableAdd 961!------------------------------------------------------------------------------ 962 963 964!> Returns a string by its name if found in the list structure 965 FUNCTION GetString( List, Name, Found ) RESULT(str) 966 TYPE(ValueList_t), POINTER :: List 967 CHARACTER(LEN=*) :: Name 968 LOGICAL, OPTIONAL :: Found 969 CHARACTER(LEN=MAX_NAME_LEN) :: str 970 971 str = ListGetString( List, Name, Found ) 972 END FUNCTION GetString 973 974 975!> Returns an integer by its name if found in the list structure 976 FUNCTION GetInteger( List, Name, Found ) RESULT(i) 977 TYPE(ValueList_t), POINTER :: List 978 CHARACTER(LEN=*) :: Name 979 LOGICAL, OPTIONAL :: Found 980 981 INTEGER :: i 982 983 i = ListGetInteger( List, Name, Found ) 984 END FUNCTION GetInteger 985 986 987!> Returns a logical flag by its name if found in the list structure, otherwise false 988 FUNCTION GetLogical( List, Name, Found ) RESULT(l) 989 TYPE(ValueList_t), POINTER :: List 990 CHARACTER(LEN=*) :: Name 991 LOGICAL, OPTIONAL :: Found 992 993 LOGICAL :: l 994 995 l = ListGetLogical( List, Name, Found ) 996 END FUNCTION GetLogical 997 998 999!> Returns a constant real by its name if found in the list structure 1000 RECURSIVE FUNCTION GetConstReal( List, Name, Found,x,y,z ) RESULT(r) 1001 TYPE(ValueList_t), POINTER :: List 1002 CHARACTER(LEN=*) :: Name 1003 LOGICAL, OPTIONAL :: Found 1004 REAL(KIND=dp), OPTIONAL :: x,y,z 1005 1006 REAL(KIND=dp) :: r,xx,yy,zz 1007 1008 xx = 0 1009 yy = 0 1010 zz = 0 1011 IF ( PRESENT( x ) ) xx = x 1012 IF ( PRESENT( y ) ) yy = y 1013 IF ( PRESENT( z ) ) zz = z 1014 1015 r = ListGetConstReal( List, Name, Found,xx,yy,zz ) 1016 END FUNCTION GetConstReal 1017 1018 1019!> Returns a real that may depend on global variables such as time, or timestep size, 1020!! by its name if found in the list structure 1021 RECURSIVE FUNCTION GetCReal( List, Name, Found ) RESULT(s) 1022 TYPE(ValueList_t), POINTER :: List 1023 CHARACTER(LEN=*) :: Name 1024 LOGICAL, OPTIONAL :: Found 1025 INTEGER, TARGET :: Dnodes(1) 1026 INTEGER, POINTER :: NodeIndexes(:) 1027 1028 REAL(KIND=dp) :: s 1029 REAL(KIND=dp), POINTER CONTIG :: x(:) 1030 TYPE(Element_t), POINTER :: Element 1031 1032 INTEGER :: n, nthreads, thread, istat 1033 1034 IF ( PRESENT( Found ) ) Found = .FALSE. 1035 1036 NodeIndexes => Dnodes 1037 n = 1 1038 NodeIndexes(n) = 1 1039 1040 x => GetValueStore(n) 1041 x(1:n) = REAL(0, dp) 1042 IF( ASSOCIATED(List) ) THEN 1043 IF ( ASSOCIATED(List % Head) ) THEN 1044 x(1:n) = ListGetReal( List, Name, n, NodeIndexes, Found ) 1045 END IF 1046 END IF 1047 s = x(1) 1048 END FUNCTION GetCReal 1049 1050 1051!> Returns a real by its name if found in the list structure, and in the active element. 1052 RECURSIVE FUNCTION GetReal( List, Name, Found, UElement ) RESULT(x) 1053 IMPLICIT NONE 1054 TYPE(ValueList_t), POINTER :: List 1055 CHARACTER(LEN=*) :: Name 1056 LOGICAL, OPTIONAL :: Found 1057 INTEGER, TARGET :: Dnodes(1) 1058 INTEGER, POINTER :: NodeIndexes(:) 1059 TYPE(Element_t), OPTIONAL, TARGET :: UElement 1060 1061 REAL(KIND=dp), POINTER CONTIG :: x(:) 1062 TYPE(Element_t), POINTER :: Element 1063 1064 INTEGER :: n, istat 1065 1066 IF ( PRESENT( Found ) ) Found = .FALSE. 1067 1068 Element => GetCurrentElement(UElement) 1069 1070 IF ( ASSOCIATED(Element) ) THEN 1071 n = GetElementNOFNodes(Element) 1072 NodeIndexes => Element % NodeIndexes 1073 ELSE 1074 n = 1 1075 NodeIndexes => Dnodes 1076 NodeIndexes(1) = 1 1077 END IF 1078 1079 x => GetValueStore(n) 1080 x(1:n) = REAL(0, dp) 1081 IF( ASSOCIATED(List) ) THEN 1082 IF ( ASSOCIATED(List % Head) ) THEN 1083 x(1:n) = ListGetReal( List, Name, n, NodeIndexes, Found ) 1084 END IF 1085 END IF 1086 END FUNCTION GetReal 1087 1088 RECURSIVE SUBROUTINE GetRealValues( List, Name, Values, Found, UElement ) 1089 IMPLICIT NONE 1090 TYPE(ValueList_t), POINTER :: List 1091 CHARACTER(LEN=*) :: Name 1092 REAL(KIND=dp) CONTIG :: Values(:) 1093 LOGICAL, OPTIONAL :: Found 1094 TYPE(Element_t), OPTIONAL, TARGET :: UElement 1095 1096 ! Variables 1097 INTEGER, TARGET :: Dnodes(1) 1098 INTEGER, POINTER CONTIG :: NodeIndexes(:) 1099 TYPE(Element_t), POINTER :: Element 1100 INTEGER :: n, istat 1101 1102 IF ( PRESENT( Found ) ) Found = .FALSE. 1103 1104 Element => GetCurrentElement(UElement) 1105 1106 IF ( ASSOCIATED(Element) ) THEN 1107 n = GetElementNOFNodes(Element) 1108 NodeIndexes => Element % NodeIndexes 1109 ELSE 1110 n = 1 1111 NodeIndexes => Dnodes 1112 NodeIndexes(1) = 1 1113 END IF 1114 1115 IF( ASSOCIATED(List) ) THEN 1116 IF ( ASSOCIATED(List % Head) ) THEN 1117 Values(1:n) = ListGetReal( List, Name, n, NodeIndexes, Found ) 1118 END IF 1119 END IF 1120 END SUBROUTINE GetRealValues 1121 1122 1123!> Returns a material property from either of the parents of the current boundary element 1124 RECURSIVE FUNCTION GetParentMatProp( Name, UElement, Found, UParent ) RESULT(x) 1125 CHARACTER(LEN=*) :: Name 1126 TYPE(Element_t), OPTIONAL, TARGET :: UElement 1127 LOGICAL, OPTIONAL :: Found 1128 TYPE(Element_t), OPTIONAL, POINTER :: UParent 1129 1130 REAL(KIND=dp), POINTER CONTIG :: x(:) 1131 INTEGER, POINTER :: Indexes(:) 1132 LOGICAL :: GotIt, GotMat 1133 INTEGER :: n, leftright, mat_id 1134 TYPE(ValueList_t), POINTER :: Material 1135 TYPE(Element_t), POINTER :: Element, Parent 1136 1137 Element => GetCurrentElement(Uelement) 1138 1139 IF( .NOT. ASSOCIATED( Element ) ) THEN 1140 CALL Warn('GetParentMatProp','Element not associated!') 1141 END IF 1142 1143 IF( PRESENT(UParent) ) NULLIFY( UParent ) 1144 1145 n = GetElementNOFNodes(Element) 1146 Indexes => Element % NodeIndexes 1147 1148 x => GetValueStore(n) 1149 x(1:n) = REAL(0, dp) 1150 1151 IF(.NOT. ASSOCIATED( Element % BoundaryInfo ) ) THEN 1152 CALL Warn('GetParentMatProp','Boundary element needs parent information!') 1153 RETURN 1154 END IF 1155 1156 1157 Gotit = .FALSE. 1158 DO leftright = 1, 2 1159 1160 IF( leftright == 1) THEN 1161 Parent => Element % BoundaryInfo % Left 1162 ELSE 1163 Parent => Element % BoundaryInfo % Right 1164 END IF 1165 1166 IF( ASSOCIATED(Parent) ) THEN 1167 1168 GotMat = .FALSE. 1169 IF( Parent % BodyId > 0 .AND. Parent % BodyId <= CurrentModel % NumberOfBodies ) THEN 1170 mat_id = ListGetInteger( CurrentModel % Bodies(Parent % BodyId) % Values,'Material',GotMat) 1171 ELSE 1172 CALL Warn('GetParentMatProp','Invalid parent BodyId '//TRIM(I2S(Parent % BodyId))//& 1173 ' for element '//TRIM(I2S(Parent % ElementIndex))) 1174 CYCLE 1175 END IF 1176 1177 IF(.NOT. GotMat) THEN 1178 CALL Warn('GetParentMatProp','Parent body '//TRIM(I2S(Parent % BodyId))//' does not have material associated!') 1179 END IF 1180 1181 IF( mat_id > 0 .AND. mat_id <= CurrentModel % NumberOfMaterials ) THEN 1182 Material => CurrentModel % Materials(mat_id) % Values 1183 ELSE 1184 CALL Warn('GetParentMatProp','Material index '//TRIM(I2S(mat_id))//' not associated to material list') 1185 CYCLE 1186 END IF 1187 1188 IF( .NOT. ASSOCIATED( Material ) ) CYCLE 1189 1190 IF ( ListCheckPresent( Material,Name) ) THEN 1191 x(1:n) = ListGetReal(Material, Name, n, Indexes) 1192 IF( PRESENT( UParent ) ) UParent => Parent 1193 Gotit = .TRUE. 1194 EXIT 1195 END IF 1196 END IF 1197 END DO 1198 1199 IF( PRESENT( Found ) ) THEN 1200 Found = GotIt 1201 ELSE IF(.NOT. GotIt) THEN 1202 CALL Warn('GetParentMatProp','Property '//TRIM(Name)//' not in either parents!') 1203 END IF 1204 1205 END FUNCTION GetParentMatProp 1206 1207 1208!> Returns a constant real array by its name if found in the list structure. 1209 RECURSIVE SUBROUTINE GetConstRealArray( List, x, Name, Found, UElement ) 1210 TYPE(ValueList_t), POINTER :: List 1211 REAL(KIND=dp), POINTER :: x(:,:) 1212 CHARACTER(LEN=*) :: Name 1213 LOGICAL, OPTIONAL :: Found 1214 TYPE(Element_t), OPTIONAL, TARGET :: UElement 1215 1216 IF ( PRESENT( Found ) ) Found = .FALSE. 1217 IF(ASSOCIATED(List)) THEN 1218 IF ( ASSOCIATED(List % Head) ) THEN 1219 x => ListGetConstRealArray( List, Name, Found ) 1220 END IF 1221 END IF 1222 END SUBROUTINE GetConstRealArray 1223 1224!> Returns a real array by its name if found in the list structure, and in the active element. 1225 RECURSIVE SUBROUTINE GetRealArray( List, x, Name, Found, UElement ) 1226 REAL(KIND=dp), POINTER :: x(:,:,:) 1227 TYPE(ValueList_t), POINTER :: List 1228 CHARACTER(LEN=*) :: Name 1229 LOGICAL, OPTIONAL :: Found 1230 TYPE(Element_t), OPTIONAL, TARGET :: UElement 1231 1232 TYPE(Element_t), POINTER :: Element 1233 1234 INTEGER :: n 1235 1236 IF ( PRESENT( Found ) ) Found = .FALSE. 1237 1238 Element => GetCurrentElement(UElement) 1239 1240 n = GetElementNOFNodes( Element ) 1241 IF ( ASSOCIATED(List) ) THEN 1242 IF ( ASSOCIATED(List % Head) ) THEN 1243 CALL ListGetRealArray( List, Name, x, n, Element % NodeIndexes, Found ) 1244 END IF 1245 END IF 1246 END SUBROUTINE GetRealArray 1247 1248!> Returns a real vector by its name if found in the list structure, and in the active element. 1249 RECURSIVE SUBROUTINE GetRealVector( List, x, Name, Found, UElement ) 1250 REAL(KIND=dp) :: x(:,:) 1251 TYPE(ValueList_t), POINTER :: List 1252 CHARACTER(LEN=*) :: Name 1253 LOGICAL, OPTIONAL :: Found 1254 TYPE(Element_t), OPTIONAL, TARGET :: UElement 1255 1256 TYPE(Element_t), POINTER :: Element 1257 1258 INTEGER :: n 1259 1260 x = 0._dp 1261 IF ( PRESENT( Found ) ) Found = .FALSE. 1262 1263 Element => GetCurrentElement(UElement) 1264 1265 n = GetElementNOFNodes( Element ) 1266 IF ( ASSOCIATED(List) ) THEN 1267 IF ( ASSOCIATED(List % Head) ) THEN 1268 CALL ListGetRealvector( List, Name, x, n, Element % NodeIndexes, Found ) 1269 END IF 1270 END IF 1271 END SUBROUTINE GetRealVector 1272 1273!> Returns a complex vector by its name if found in the list structure, and in the active element. 1274 RECURSIVE SUBROUTINE GetComplexVector( List, x, Name, Found, UElement ) 1275 COMPLEX(KIND=dp) :: x(:,:) 1276 TYPE(ValueList_t), POINTER :: List 1277 CHARACTER(LEN=*) :: Name 1278 LOGICAL, OPTIONAL :: Found 1279 TYPE(Element_t), OPTIONAL, TARGET :: UElement 1280 1281 TYPE(Element_t), POINTER :: Element 1282 LOGICAL :: lFound 1283 INTEGER :: n 1284 REAL(KIND=dp), ALLOCATABLE :: xr(:,:) 1285 1286 x = 0._dp 1287 IF ( PRESENT( Found ) ) Found = .FALSE. 1288 1289 Element => GetCurrentElement(UElement) 1290 1291 n = GetElementNOFNodes( Element ) 1292 IF ( ASSOCIATED(List) ) THEN 1293 IF ( ASSOCIATED(List % Head) ) THEN 1294 ALLOCATE(xr(SIZE(x,1),SIZE(x,2))) 1295 CALL ListGetRealvector( List, Name, xr, n, & 1296 Element % NodeIndexes, lFound ) 1297 IF(PRESENT(Found)) Found=lFound 1298 x = xr 1299 CALL ListGetRealvector( List, TRIM(Name)//' im', & 1300 xr, n, Element % NodeIndexes, lFound ) 1301 IF(PRESENT(Found)) Found=Found.OR.lFound 1302 x = CMPLX(REAL(x), xr) 1303 END IF 1304 END IF 1305 END SUBROUTINE GetComplexVector 1306 1307 1308!> Set a named elementwise property (real-valued) to the active element or 1309!> given element 1310 SUBROUTINE SetElementProperty( Name, Values, UElement ) 1311 CHARACTER(LEN=*) :: Name 1312 REAL(KIND=dp) :: Values(:) 1313 TYPE(Element_t), POINTER, OPTIONAL :: UElement 1314 1315 TYPE(ElementData_t), POINTER :: p 1316 1317 TYPE(Element_t), POINTER :: Element 1318 1319 Element => GetCurrentElement(UElement) 1320 1321 p => Element % PropertyData 1322 DO WHILE( ASSOCIATED(p) ) 1323 IF ( Name==p % Name ) EXIT 1324 p => p % Next 1325 END DO 1326 1327 IF ( ASSOCIATED(p) ) THEN 1328 IF ( SIZE(P % Values) == SIZE(Values) ) THEN 1329 P % Values = Values 1330 ELSE 1331 DEALLOCATE( P % Values ) 1332 ALLOCATE( P % Values(SIZE(Values)) ) 1333 P % Values = Values 1334 END IF 1335 ELSE 1336 ALLOCATE(p) 1337 ALLOCATE( P % Values(SIZE(Values)) ) 1338 p % Values = Values 1339 p % Name = Name 1340 p % Next => Element % PropertyData 1341 Element % PropertyData => p 1342 END IF 1343 END SUBROUTINE SetElementProperty 1344 1345!> Get a named elementwise property (real-valued) from the active element or 1346!> from a given element 1347 FUNCTION GetElementProperty( Name, UElement ) RESULT(Values) 1348 CHARACTER(LEN=*) :: Name 1349 REAL(KIND=dp), POINTER :: Values(:) 1350 TYPE(Element_t), POINTER, OPTIONAL :: UElement 1351 1352 TYPE(ElementData_t), POINTER :: p 1353 1354 TYPE(Element_t), POINTER :: Element 1355 1356 Element => GetCurrentElement(UElement) 1357 1358 Values => NULL() 1359 p=> Element % PropertyData 1360 1361 DO WHILE( ASSOCIATED(p) ) 1362 IF ( Name==p % Name ) THEN 1363 Values => p % Values 1364 RETURN 1365 END IF 1366 p => p % Next 1367 END DO 1368 END FUNCTION GetElementProperty 1369 1370 1371!> Get a handle to the active element from the list of all active elements 1372 FUNCTION GetActiveElement(t,USolver) RESULT(Element) 1373 INTEGER :: t 1374 TYPE(Element_t), POINTER :: Element 1375 TYPE( Solver_t ), OPTIONAL, TARGET :: USolver 1376 1377 TYPE( Solver_t ), POINTER :: Solver 1378 INTEGER :: ind 1379 1380 Solver => CurrentModel % Solver 1381 IF ( PRESENT( USolver ) ) Solver => USolver 1382 1383 IF ( t > 0 .AND. t <= Solver % NumberOfActiveElements ) THEN 1384 ! Check if colouring is really used by the solver 1385 IF( Solver % CurrentColour > 0 .AND. & 1386 ASSOCIATED( Solver % ColourIndexList ) ) THEN 1387 ind = Solver % ActiveElements( & 1388 Solver % ColourIndexList % ind(& 1389 Solver % ColourIndexList % ptr(Solver % CurrentColour)+(t-1) ) ) 1390 ELSE 1391 ind = Solver % ActiveElements(t) 1392 END IF 1393 1394 Element => Solver % Mesh % Elements( ind ) 1395 1396#ifdef _OPENMP 1397 IF (omp_in_parallel()) THEN 1398 CurrentElementThread => Element 1399 ELSE 1400 ! May be used by user functions, not thread safe 1401 CurrentModel % CurrentElement => Element 1402 END IF 1403#else 1404 ! May be used by user functions, not thread safe 1405 CurrentModel % CurrentElement => Element 1406#endif 1407 ELSE 1408 WRITE( Message, * ) 'Invalid element number requested: ', t 1409 CALL Fatal( 'GetActiveElement', Message ) 1410 END IF 1411 END FUNCTION GetActiveElement 1412 1413 1414!> Get a handle to a boundary element from the list of all boundary elements 1415 FUNCTION GetBoundaryElement(t,USolver) RESULT(Element) 1416 INTEGER :: t 1417 TYPE(Element_t), POINTER :: Element 1418 TYPE( Solver_t ), OPTIONAL, TARGET :: USolver 1419 TYPE( Solver_t ), POINTER :: Solver 1420 INTEGER :: ind 1421 1422 Solver => CurrentModel % Solver 1423 IF ( PRESENT( USolver ) ) Solver => USolver 1424 1425 IF ( t > 0 .AND. t <= Solver % Mesh % NumberOfBoundaryElements ) THEN 1426 ! Check if colouring is really used by the solver 1427 IF( Solver % CurrentBoundaryColour > 0 .AND. & 1428 ASSOCIATED( Solver % BoundaryColourIndexList ) ) THEN 1429 ind = Solver % BoundaryColourIndexList % ind( & 1430 Solver % BoundaryColourIndexList % ptr(Solver % CurrentBoundaryColour)+(t-1)) 1431 ELSE 1432 ind = t 1433 END IF 1434 1435 ! Element => Solver % Mesh % Elements( Solver % Mesh % NumberOfBulkElements+t ) 1436 Element => Solver % Mesh % Elements( Solver % Mesh % NumberOfBulkElements + ind ) 1437#ifdef _OPENMP 1438 IF (omp_in_parallel()) THEN 1439 ! May be used by user functions, thread safe 1440 CurrentElementThread => Element 1441 ELSE 1442 CurrentModel % CurrentElement => Element 1443 END IF 1444#else 1445 CurrentModel % CurrentElement => Element 1446#endif 1447 ELSE 1448 WRITE( Message, * ) 'Invalid element number requested: ', t 1449 CALL Fatal( 'GetBoundaryElement', Message ) 1450 END IF 1451 END FUNCTION GetBoundaryElement 1452 1453 1454!> Check if the boundary element is active in the current solve 1455 FUNCTION ActiveBoundaryElement(UElement,USolver,DGBoundary) RESULT(l) 1456 TYPE(Element_t), OPTIONAL, TARGET :: UElement 1457 TYPE(Solver_t), OPTIONAL, TARGET :: USolver 1458 LOGICAL, OPTIONAL :: DGBoundary 1459 1460 LOGICAL :: l, DGb 1461 INTEGER :: n, n2 1462 INTEGER, POINTER :: Indexes(:) 1463 1464 TYPE( Solver_t ), POINTER :: Solver 1465 TYPE(Element_t), POINTER :: Element, P1, P2 1466 1467 Solver => CurrentModel % Solver 1468 IF ( PRESENT( USolver ) ) Solver => USolver 1469 1470 Element => GetCurrentElement(UElement) 1471 1472 Indexes => GetIndexStore() 1473 n = GetElementDOFs( Indexes, Element, Solver ) 1474 1475 DGb = Solver % DG .AND. PRESENT(DGboundary) 1476 IF(DGb) DGb = DGboundary 1477 1478 IF (DGb) THEN 1479 P1 => Element % BoundaryInfo % Left 1480 P2 => Element % BoundaryInfo % Right 1481 IF ( ASSOCIATED(P1).AND.ASSOCIATED(P2) ) THEN 1482 n = P1 % Type % NumberOfNodes 1483 l = ALL(Solver % Variable % Perm(Indexes(1:n)) > 0) 1484 IF (.NOT.l) THEN 1485 n2 = P2 % Type % NumberOfNodes 1486 l = ALL(Solver % Variable % Perm(Indexes(n+1:n+n2)) > 0) 1487 END IF 1488 ELSE 1489 l = ALL(Solver % Variable % Perm(Indexes(1:n)) > 0) 1490 END IF 1491 ELSE 1492 IF (isActivePElement(Element)) n=GetElementNOFNOdes(Element) 1493 l = ALL(Solver % Variable % Perm(Indexes(1:n)) > 0) 1494 END IF 1495 END FUNCTION ActiveBoundaryElement 1496 1497 1498!> Return the element code in Elmer convention of the active element 1499 FUNCTION GetElementCode( Element ) RESULT(etype) 1500 INTEGER :: etype 1501 TYPE(Element_t), OPTIONAL :: Element 1502 TYPE(Element_t), POINTER :: CurrElement 1503 1504 CurrElement => GetCurrentElement(Element) 1505 etype = CurrElement % TYPE % ElementCode 1506 END FUNCTION GetElementCode 1507 1508!> Return the element dimension in Elmer convention of the active element 1509 FUNCTION GetElementDim( Element ) RESULT(edim) 1510 INTEGER :: edim 1511 TYPE(Element_t), OPTIONAL :: Element 1512 TYPE(Element_t), POINTER :: CurrElement 1513 INTEGER :: etype 1514 1515 CurrElement => GetCurrentElement(Element) 1516 etype = CurrElement % TYPE % ElementCode 1517 IF( etype >= 500 ) THEN 1518 edim = 3 1519 ELSE IF( etype >= 300 ) THEN 1520 edim = 2 1521 ELSE IF( etype >= 200 ) THEN 1522 edim = 1 1523 ELSE 1524 edim = 0 1525 END IF 1526 END FUNCTION GetElementDim 1527 1528 1529!> Return the element family in Elmer convention of the active element 1530 FUNCTION GetElementFamily( Element ) RESULT(family) 1531 INTEGER :: family 1532 TYPE(Element_t), OPTIONAL :: Element 1533 TYPE(Element_t), POINTER :: CurrElement 1534 1535 CurrElement => GetCurrentElement(Element) 1536 family = CurrElement % TYPE % ElementCode / 100 1537 END FUNCTION GetElementFamily 1538 1539 1540!> Return the number of corners nodes i.e. the number of dofs for the lowest order element 1541 FUNCTION GetElementCorners( Element ) RESULT(corners) 1542 INTEGER :: corners 1543 TYPE(Element_t), OPTIONAL :: Element 1544 TYPE(Element_t), POINTER :: CurrElement 1545 1546 CurrElement => GetCurrentElement(Element) 1547 corners = CurrElement % TYPE % ElementCode / 100 1548 IF( corners >= 5 .AND. corners <= 7 ) THEN 1549 corners = corners - 1 1550 END IF 1551 END FUNCTION GetElementCorners 1552 1553!> Return true if the element is a possible flux element 1554!> Needed to skip nodal elements in 2D and 3D boundary condition setting. 1555 FUNCTION PossibleFluxElement( Element, Mesh ) RESULT(possible) 1556 LOGICAL :: possible 1557 TYPE(Element_t), OPTIONAL :: Element 1558 TYPE(Mesh_t), OPTIONAL :: Mesh 1559 INTEGER :: MeshDim, family 1560 1561 ! Orphan elements are not currently present in the mesh so any 1562 ! boundary condition that exists is a possible flux element also. 1563 ! Thus this routine is more or less obsolete. 1564 possible = .TRUE. 1565 1566 RETURN 1567 1568 1569 IF( PRESENT( Mesh ) ) THEN 1570 MeshDim = Mesh % MeshDim 1571 ELSE 1572 MeshDim = CurrentModel % Solver % Mesh % MeshDim 1573 END IF 1574 1575 family = GetElementFamily( Element ) 1576 1577 ! This is not a generic rule but happens to be true for all combinations 1578 ! 3D: families 3 and 4 1579 ! 2D: family 2 1580 ! 1D: family 1 1581 possible = ( MeshDim <= family ) 1582 1583 END FUNCTION PossibleFluxElement 1584 1585 1586!> Return the number of nodes in the active element 1587 FUNCTION GetElementNOFNodes( Element ) RESULT(n) 1588 INTEGER :: n 1589 TYPE(Element_t), OPTIONAL :: Element 1590 TYPE(Element_t), POINTER :: CurrElement 1591 1592 CurrElement => GetCurrentElement(Element) 1593 n = CurrElement % TYPE % NumberOfNodes 1594 END FUNCTION GetELementNOFNodes 1595 1596 1597!> Return the number of element degrees of freedom 1598 FUNCTION GetElementNOFDOFs( UElement,USolver ) RESULT(n) 1599 INTEGER :: n 1600 TYPE(Solver_t), OPTIONAL, TARGET :: USolver 1601 TYPE(Element_t), OPTIONAL, TARGET :: UElement 1602 1603 TYPE(Element_t), POINTER :: Element 1604 TYPE(Solver_t), POINTER :: Solver 1605 1606 INTEGER :: i,j, id, ElemFamily 1607 LOGICAL :: Found, GB, NeedEdges 1608 1609 Element => GetCurrentElement( UElement ) 1610 ElemFamily = GetElementFamily(Element) 1611 1612 IF ( PRESENT( USolver ) ) THEN 1613 Solver => USolver 1614 ELSE 1615 Solver => CurrentModel % Solver 1616 END IF 1617 1618 n = 0 1619 IF( Solver % DG ) THEN 1620 n = Element % DGDOFs 1621 IF ( n>0 ) RETURN 1622 END IF 1623 1624 1625 id =Element % BodyId 1626 IF ( Id==0 .AND. ASSOCIATED(Element % BoundaryInfo) ) THEN 1627 IF ( ASSOCIATED(Element % BoundaryInfo % Left) ) & 1628 id = Element % BoundaryInfo % Left % BodyId 1629 1630 IF ( ASSOCIATED(Element % BoundaryInfo % Right) ) & 1631 id = Element % BoundaryInfo % Right % BodyId 1632 END IF 1633 IF ( Id==0 ) id=1 1634 1635 IF ( Solver % Def_Dofs(ElemFamily,id,1)>0 ) n = Element % NDOFs 1636 1637 NeedEdges = .FALSE. 1638 DO i=2,SIZE(Solver % Def_Dofs,3) 1639 IF (Solver % Def_Dofs(ElemFamily, id, i)>=0) THEN 1640 NeedEdges = .TRUE. 1641 EXIT 1642 END IF 1643 END DO 1644 IF ( .NOT. NeedEdges ) RETURN 1645 1646 IF ( ASSOCIATED( Element % EdgeIndexes ) ) THEN 1647 IF ( Solver % Mesh % MaxEdgeDOFs == Solver % Mesh % MinEdgeDOFs ) THEN 1648 n = n + Element % Type % NumberOfEdges * Solver % Mesh % MaxEdgeDOFs 1649 ELSE 1650!DIR$ IVDEP 1651 DO j=1,Element % Type % NumberOFEdges 1652 n = n + Solver % Mesh % Edges(Element % EdgeIndexes(j)) % BDOFs 1653 END DO 1654 END IF 1655 END IF 1656 1657 IF ( ASSOCIATED( Element % FaceIndexes ) ) THEN 1658 IF ( Solver % Mesh % MaxFaceDOFs == Solver % Mesh % MinFaceDOFs ) THEN 1659 n = n + Element % Type % NumberOfFaces * Solver % Mesh % MaxFaceDOFs 1660 ELSE 1661!DIR$ IVDEP 1662 DO j=1,Element % Type % NumberOFFaces 1663 n = n + Solver % Mesh % Faces( Element % FaceIndexes(j) ) % BDOFs 1664 END DO 1665 END IF 1666 END IF 1667 1668 !GB = ListGetLogical( Solver % Values, 'Bubbles in Global System', Found ) 1669 !IF (.NOT.Found) GB = .TRUE. 1670 GB = Solver % GlobalBubbles 1671 1672 IF ( GB .OR. ASSOCIATED(Element % BoundaryInfo) ) n=n+MAX(0,Element % BDOFs) 1673 END FUNCTION GetElementNOFDOFs 1674 1675 1676!> In addition to returning the number of degrees of freedom associated with 1677!> the element, the indexes of the degrees of freedom are also returned. 1678 FUNCTION GetElementDOFs( Indexes, UElement, USolver,NotDG ) RESULT(NB) 1679 TYPE(Element_t), OPTIONAL, TARGET :: UElement 1680 TYPE(Solver_t), OPTIONAL, TARGET :: USolver 1681 INTEGER :: Indexes(:) 1682 LOGICAL, OPTIONAL :: NotDG 1683 1684 TYPE(Solver_t), POINTER :: Solver 1685 TYPE(Element_t), POINTER :: Element, Parent, Edge, Face 1686 1687 LOGICAL :: Found, GB, DGdisable, NeedEdges 1688 INTEGER :: nb,i,j,k,id,EDOFs, FDOFs, BDOFs,FaceDOFs, EdgeDOFs, BubbleDOFs, Ind, ElemFamily 1689 1690 Element => GetCurrentElement(UElement) 1691 ElemFamily = GetElementFamily(Element) 1692 1693 IF ( PRESENT( USolver ) ) THEN 1694 Solver => USolver 1695 ELSE 1696 Solver => CurrentModel % Solver 1697 END IF 1698 1699 NB = 0 1700 1701 DGDisable=.FALSE. 1702 IF (PRESENT(NotDG)) DGDisable=NotDG 1703 1704 IF ( .NOT. DGDisable .AND. Solver % DG ) THEN 1705 DO i=1,Element % DGDOFs 1706 NB = NB + 1 1707 Indexes(NB) = Element % DGIndexes(i) 1708 END DO 1709 1710 IF ( ASSOCIATED( Element % BoundaryInfo ) ) THEN 1711 IF ( ASSOCIATED( Element % BoundaryInfo % Left ) ) THEN 1712 DO i=1,Element % BoundaryInfo % Left % DGDOFs 1713 NB = NB + 1 1714 Indexes(NB) = Element % BoundaryInfo % Left % DGIndexes(i) 1715 END DO 1716 END IF 1717 IF ( ASSOCIATED( Element % BoundaryInfo % Right ) ) THEN 1718 DO i=1,Element % BoundaryInfo % Right % DGDOFs 1719 NB = NB + 1 1720 Indexes(NB) = Element % BoundaryInfo % Right % DGIndexes(i) 1721 END DO 1722 END IF 1723 END IF 1724 1725 IF ( NB > 0 ) RETURN 1726 END IF 1727 1728 id =Element % BodyId 1729 IF ( Id==0 .AND. ASSOCIATED(Element % BoundaryInfo) ) THEN 1730 IF ( ASSOCIATED(Element % BoundaryInfo % Left) ) & 1731 id = Element % BoundaryInfo % Left % BodyId 1732 1733 IF ( ASSOCIATED(Element % BoundaryInfo % Right) ) & 1734 id = Element % BoundaryInfo % Right % BodyId 1735 END IF 1736 IF ( id == 0 ) id=1 1737 1738 IF ( Solver % Def_Dofs(ElemFamily,id,1)>0 ) THEN 1739 DO i=1,Element % NDOFs 1740 NB = NB + 1 1741 Indexes(NB) = Element % NodeIndexes(i) 1742 END DO 1743 END IF 1744 1745 ! default for nodal elements, if no solver active: 1746 ! ------------------------------------------------ 1747 IF(.NOT.ASSOCIATED(Solver)) RETURN 1748 IF(.NOT.ASSOCIATED(Solver % Mesh)) RETURN 1749 1750 NeedEdges = .FALSE. 1751 DO i=2,SIZE(Solver % Def_Dofs,3) 1752 IF (Solver % Def_Dofs(ElemFamily, id, i)>=0) THEN 1753 NeedEdges = .TRUE. 1754 EXIT 1755 END IF 1756 END DO 1757 IF ( .NOT. NeedEdges ) RETURN 1758 1759 FaceDOFs = Solver % Mesh % MaxFaceDOFs 1760 EdgeDOFs = Solver % Mesh % MaxEdgeDOFs 1761 BubbleDOFs = Solver % Mesh % MaxBDOFs 1762 1763 IF ( ASSOCIATED(Element % EdgeIndexes) ) THEN 1764 DO j=1,Element % TYPE % NumberOFEdges 1765 EDOFs = Solver % Mesh % Edges(Element % EdgeIndexes(j)) % BDOFs 1766 DO i=1,EDOFs 1767 NB = NB + 1 1768 Indexes(NB) = EdgeDOFs*(Element % EdgeIndexes(j)-1) + & 1769 i + Solver % Mesh % NumberOfNodes 1770 END DO 1771 END DO 1772 END IF 1773 1774 IF ( ASSOCIATED( Element % FaceIndexes ) ) THEN 1775 DO j=1,Element % TYPE % NumberOFFaces 1776 FDOFs = Solver % Mesh % Faces( Element % FaceIndexes(j) ) % BDOFs 1777 DO i=1,FDOFs 1778 NB = NB + 1 1779 Indexes(NB) = FaceDOFs*(Element % FaceIndexes(j)-1) + i + & 1780 Solver % Mesh % NumberOfNodes + EdgeDOFs*Solver % Mesh % NumberOfEdges 1781 END DO 1782 END DO 1783 END IF 1784 1785 GB = Solver % GlobalBubbles 1786 1787 IF ( ASSOCIATED(Element % BoundaryInfo) ) THEN 1788 Parent => Element % BoundaryInfo % Left 1789 IF (.NOT.ASSOCIATED(Parent) ) & 1790 Parent => Element % BoundaryInfo % Right 1791 IF (.NOT.ASSOCIATED(Parent) ) RETURN 1792 1793 SELECT CASE(GetElementFamily(Element)) 1794 CASE(2) 1795 IF ( ASSOCIATED(Parent % EdgeIndexes) ) THEN 1796 IF ( isActivePElement(Element) ) THEN 1797 Ind=Element % PDefs % LocalNumber 1798 ELSE 1799 DO Ind=1,Parent % TYPE % NumberOfEdges 1800 Edge => Solver % Mesh % Edges(Parent % EdgeIndexes(ind)) 1801 k = 0 1802 DO i=1,Edge % TYPE % NumberOfNodes 1803 DO j=1,Element % TYPE % NumberOfNodes 1804 IF ( Edge % NodeIndexes(i)==Element % NodeIndexes(j) ) k=k+1 1805 END DO 1806 END DO 1807 IF ( k==Element % TYPE % NumberOfNodes) EXIT 1808 END DO 1809 END IF 1810 1811 EDOFs = Element % BDOFs 1812 DO i=1,EDOFs 1813 NB = NB + 1 1814 Indexes(NB) = EdgeDOFs*(Parent % EdgeIndexes(Ind)-1) + & 1815 i + Solver % Mesh % NumberOfNodes 1816 END DO 1817 END IF 1818 1819 CASE(3,4) 1820 IF ( ASSOCIATED( Parent % FaceIndexes ) ) THEN 1821 IF ( isActivePElement(Element) ) THEN 1822 Ind=Element % PDefs % LocalNumber 1823 ELSE 1824 DO Ind=1,Parent % TYPE % NumberOfFaces 1825 Face => Solver % Mesh % Faces(Parent % FaceIndexes(ind)) 1826 k = 0 1827 DO i=1,Face % TYPE % NumberOfNodes 1828 DO j=1,Element % TYPE % NumberOfNodes 1829 IF ( Face % NodeIndexes(i)==Element % NodeIndexes(j)) k=k+1 1830 END DO 1831 END DO 1832 IF ( k==Face % TYPE % NumberOfNodes) EXIT 1833 END DO 1834 END IF 1835 1836 FDOFs = Element % BDOFs 1837 DO i=1,FDOFs 1838 NB = NB + 1 1839 Indexes(NB) = FaceDOFs*(Parent % FaceIndexes(Ind)-1) + i + & 1840 Solver % Mesh % NumberOfNodes + EdgeDOFs*Solver % Mesh % NumberOfEdges 1841 END DO 1842 END IF 1843 END SELECT 1844 ELSE IF ( GB ) THEN 1845 IF ( ASSOCIATED(Element % BubbleIndexes) ) THEN 1846 DO i=1,Element % BDOFs 1847 NB = NB + 1 1848 Indexes(NB) = FaceDOFs*Solver % Mesh % NumberOfFaces + & 1849 Solver % Mesh % NumberOfNodes + EdgeDOFs*Solver % Mesh % NumberOfEdges + & 1850 Element % BubbleIndexes(i) 1851 END DO 1852 END IF 1853 END IF 1854 END FUNCTION GetElementDOFs 1855 1856 1857! ----------------------------------------------------------------------------- 1858!> Returns the number of bubble degrees of freedom in the active element. 1859!> If the sif file contains more than one solver section 1860!> with each of them having their own specification of the "Element" 1861!> keyword, the returned value may not be the number of bubbles that 1862!> should be assigned to the solver. With the optional argument 1863!> Update = .TRUE., the correct solver-wise bubble count can be returned and 1864!> the bubble count assigned to the Element argument is updated. 1865! ----------------------------------------------------------------------------- 1866 FUNCTION GetElementNOFBDOFs( Element, USolver, Update ) RESULT(n) 1867! ----------------------------------------------------------------------------- 1868 INTEGER :: n 1869 TYPE(Element_t), OPTIONAL :: Element 1870 TYPE(Solver_t), OPTIONAL, POINTER :: USolver 1871 LOGICAL, OPTIONAL :: Update 1872 1873 TYPE(Element_t), POINTER :: CurrElement 1874 TYPE(Solver_t), POINTER :: Solver 1875 LOGICAL :: Found, GB, UpdateRequested 1876 INTEGER :: k 1877 1878 IF ( PRESENT( USolver ) ) THEN 1879 Solver => USolver 1880 ELSE 1881 Solver => CurrentModel % Solver 1882 END IF 1883 1884 UpdateRequested = .FALSE. 1885 IF ( PRESENT(Update) ) UpdateRequested = Update 1886 1887 !GB = ListGetLogical( Solver % Values, 'Bubbles in Global System', Found ) 1888 !IF (.NOT.Found) GB = .TRUE. 1889 GB = Solver % GlobalBubbles 1890 1891 n = 0 1892 IF ( .NOT. GB ) THEN 1893 CurrElement => GetCurrentElement(Element) 1894 IF (UpdateRequested) THEN 1895 n = Solver % Def_Dofs(GetElementFamily(CurrElement), & 1896 CurrElement % Bodyid, 5) 1897 IF ( n>=0 ) THEN 1898 CurrElement % BDOFs = n 1899 ELSE 1900 n = CurrElement % BDOFs 1901 END IF 1902 ELSE 1903 n = CurrElement % BDOFs 1904 END IF 1905 ELSE 1906 ! Rectify the bubble count assigned to the Element argument in case 1907 ! some other solver has tampered it: 1908 IF (UpdateRequested) THEN 1909 CurrElement => GetCurrentElement(Element) 1910 k = Solver % Def_Dofs(GetElementFamily(CurrElement), & 1911 CurrElement % Bodyid, 5) 1912 IF ( k>=0 ) CurrElement % BDOFs = k 1913 END IF 1914 END IF 1915 END FUNCTION GetElementNOFBDOFs 1916 1917 1918!> Returns the nodal coordinate values in the active element 1919 SUBROUTINE GetElementNodes( ElementNodes, UElement, USolver, UMesh ) 1920 TYPE(Nodes_t) :: ElementNodes 1921 TYPE(Solver_t), OPTIONAL, TARGET :: USolver 1922 TYPE(Mesh_t), OPTIONAL, TARGET :: UMesh 1923 TYPE(Element_t), OPTIONAL, TARGET :: UElement 1924 1925 INTEGER :: i,n,nd,sz,sz1 1926 INTEGER, POINTER :: Indexes(:) 1927 1928 TYPE(Solver_t), POINTER :: Solver 1929 TYPE(Mesh_t), POINTER :: Mesh 1930 TYPE(Element_t), POINTER :: Element 1931 1932 Solver => CurrentModel % Solver 1933 IF ( PRESENT( USolver ) ) Solver => USolver 1934 1935 Element => GetCurrentElement(UElement) 1936 1937 Mesh => Solver % Mesh 1938 IF ( PRESENT( UMesh ) ) Mesh => UMesh 1939 n = MAX(Mesh % MaxElementNodes,Mesh % MaxElementDOFs) 1940 1941 IF ( .NOT. ASSOCIATED( ElementNodes % x ) ) THEN 1942 ALLOCATE( ElementNodes % x(n), ElementNodes % y(n),ElementNodes % z(n) ) 1943 ELSE IF ( SIZE(ElementNodes % x)<n ) THEN 1944 DEALLOCATE(ElementNodes % x, ElementNodes % y, ElementNodes % z) 1945 ALLOCATE( ElementNodes % x(n), ElementNodes % y(n),ElementNodes % z(n) ) 1946 END IF 1947 1948 n = Element % TYPE % NumberOfNodes 1949 ElementNodes % x(1:n) = Mesh % Nodes % x(Element % NodeIndexes) 1950 ElementNodes % y(1:n) = Mesh % Nodes % y(Element % NodeIndexes) 1951 ElementNodes % z(1:n) = Mesh % Nodes % z(Element % NodeIndexes) 1952 1953 sz = SIZE(ElementNodes % x) 1954 IF ( sz > n ) THEN 1955 ElementNodes % x(n+1:sz) = 0.0_dp 1956 ElementNodes % y(n+1:sz) = 0.0_dp 1957 ElementNodes % z(n+1:sz) = 0.0_dp 1958 END IF 1959 1960 sz1 = SIZE(Mesh % Nodes % x) 1961 IF (sz1 > Mesh % NumberOfNodes) THEN 1962 Indexes => GetIndexStore() 1963 nd = GetElementDOFs(Indexes,Element,NotDG=.TRUE.) 1964 DO i=n+1,nd 1965 IF ( Indexes(i)>0 .AND. Indexes(i)<=sz1 ) THEN 1966 ElementNodes % x(i) = Mesh % Nodes % x(Indexes(i)) 1967 ElementNodes % y(i) = Mesh % Nodes % y(Indexes(i)) 1968 ElementNodes % z(i) = Mesh % Nodes % z(Indexes(i)) 1969 END IF 1970 END DO 1971 END IF 1972 END SUBROUTINE GetElementNodes 1973 1974!> Returns the nodal coordinate values in the active element 1975 SUBROUTINE GetElementNodesVec( ElementNodes, UElement, USolver, UMesh ) 1976 TYPE(Nodes_t), TARGET :: ElementNodes 1977 TYPE(Solver_t), OPTIONAL, TARGET :: USolver 1978 TYPE(Mesh_t), OPTIONAL, TARGET :: UMesh 1979 TYPE(Element_t), OPTIONAL, TARGET :: UElement 1980 1981 INTEGER :: padn, dum 1982 1983 INTEGER :: i,n,nd,sz,sz1 1984 INTEGER, POINTER CONTIG :: Indexes(:) 1985 1986 TYPE(Solver_t), POINTER :: Solver 1987 TYPE(Mesh_t), POINTER :: Mesh 1988 TYPE(Element_t), POINTER :: Element 1989 1990 Solver => CurrentModel % Solver 1991 IF ( PRESENT( USolver ) ) Solver => USolver 1992 1993 Element => GetCurrentElement(UElement) 1994 1995 IF ( PRESENT( UMesh ) ) THEN 1996 Mesh => UMesh 1997 ELSE 1998 Mesh => Solver % Mesh 1999 END IF 2000 2001 n = MAX(Mesh % MaxElementNodes,Mesh % MaxElementDOFs) 2002 padn = n 2003 2004 ! Here we could pad beginning of columns of xyz to 64-byte 2005 ! boundaries if needed as follows 2006 ! padn=NBytePad(n,STORAGE_SIZE(REAL(1,dp))/8,64) 2007 2008 IF (.NOT. ALLOCATED( ElementNodes % xyz)) THEN 2009 IF (ASSOCIATED(ElementNodes % x)) DEALLOCATE(ElementNodes % x) 2010 IF (ASSOCIATED(ElementNodes % y)) DEALLOCATE(ElementNodes % y) 2011 IF (ASSOCIATED(ElementNodes % z)) DEALLOCATE(ElementNodes % z) 2012 2013 ALLOCATE(ElementNodes % xyz(padn,3)) 2014 ElementNodes % xyz = REAL(0,dp) 2015 ElementNodes % x => ElementNodes % xyz(1:n,1) 2016 ElementNodes % y => ElementNodes % xyz(1:n,2) 2017 ElementNodes % z => ElementNodes % xyz(1:n,3) 2018 ELSE IF (SIZE(ElementNodes % xyz, 1)<padn) THEN 2019 DEALLOCATE(ElementNodes % xyz) 2020 ALLOCATE(ElementNodes % xyz(padn,3)) 2021 ElementNodes % xyz = REAL(0,dp) 2022 ElementNodes % x => ElementNodes % xyz(1:n,1) 2023 ElementNodes % y => ElementNodes % xyz(1:n,2) 2024 ElementNodes % z => ElementNodes % xyz(1:n,3) 2025 ELSE 2026 ElementNodes % x => ElementNodes % xyz(1:n,1) 2027 ElementNodes % y => ElementNodes % xyz(1:n,2) 2028 ElementNodes % z => ElementNodes % xyz(1:n,3) 2029 END IF 2030 2031 n = Element % TYPE % NumberOfNodes 2032!DIR$ IVDEP 2033 DO i=1,n 2034 ElementNodes % x(i) = Mesh % Nodes % x(Element % NodeIndexes(i)) 2035 ElementNodes % y(i) = Mesh % Nodes % y(Element % NodeIndexes(i)) 2036 ElementNodes % z(i) = Mesh % Nodes % z(Element % NodeIndexes(i)) 2037 END DO 2038 2039 sz = SIZE(ElementNodes % xyz,1) 2040 IF ( sz > n ) THEN 2041 ElementNodes % xyz(n+1:sz,1) = 0.0d0 2042 ElementNodes % xyz(n+1:sz,2) = 0.0d0 2043 ElementNodes % xyz(n+1:sz,3) = 0.0d0 2044 END IF 2045 2046 sz1 = SIZE(Mesh % Nodes % x) 2047 IF (sz1 > Mesh % NumberOfNodes) THEN 2048 Indexes => GetIndexStore() 2049 nd = GetElementDOFs(Indexes,Element,NotDG=.TRUE.) 2050!DIR$ IVDEP 2051 DO i=n+1,nd 2052 IF ( Indexes(i)>0 .AND. Indexes(i)<=sz1 ) THEN 2053 ElementNodes % x(i) = Mesh % Nodes % x(Indexes(i)) 2054 ElementNodes % y(i) = Mesh % Nodes % y(Indexes(i)) 2055 ElementNodes % z(i) = Mesh % Nodes % z(Indexes(i)) 2056 END IF 2057 END DO 2058 END IF 2059 END SUBROUTINE GetElementNodesVec 2060 2061!> Get element body id 2062!------------------------------------------------------------------------------ 2063 FUNCTION GetBody( Element ) RESULT(body_id) 2064!------------------------------------------------------------------------------ 2065 INTEGER::Body_id 2066 TYPE(Element_t), OPTIONAL :: Element 2067!------------------------------------------------------------------------------ 2068 TYPE(Element_t), POINTER :: el 2069!------------------------------------------------------------------------------ 2070 el => GetCurrentElement(Element) 2071 body_id= el % BodyId 2072!------------------------------------------------------------------------------ 2073 END FUNCTION GetBody 2074!------------------------------------------------------------------------------ 2075 2076 2077!> Get element body parameters 2078!------------------------------------------------------------------------------ 2079 FUNCTION GetBodyParams(Element) RESULT(lst) 2080!------------------------------------------------------------------------------ 2081 TYPE(ValueList_t), POINTER :: Lst 2082 TYPE(Element_t), OPTIONAL :: Element 2083!------------------------------------------------------------------------------ 2084 TYPE(Element_t), POINTER :: el 2085!------------------------------------------------------------------------------ 2086 lst => CurrentModel % Bodies(GetBody(Element)) % Values 2087!------------------------------------------------------------------------------ 2088 END FUNCTION GetBodyParams 2089!------------------------------------------------------------------------------ 2090 2091 2092!> Get the body force index of the active element 2093!------------------------------------------------------------------------------ 2094 FUNCTION GetBodyForceId( Element, Found ) RESULT(bf_id) 2095!------------------------------------------------------------------------------ 2096 LOGICAL, OPTIONAL :: Found 2097 TYPE(Element_t), OPTIONAL :: Element 2098 TYPE(Element_t), POINTER :: CurrElement 2099 2100 INTEGER :: bf_id, body_id 2101 2102 CurrElement => GetCurrentElement(Element) 2103 body_id = CurrElement % BodyId 2104 2105 IF ( PRESENT( Found ) ) THEN 2106 bf_id = ListGetInteger( CurrentModel % Bodies(body_id) % Values, & 2107 'Body Force', Found, minv=1,maxv=CurrentModel % NumberOfBodyForces ) 2108 ELSE 2109 bf_id = ListGetInteger( CurrentModel % Bodies(body_id) % Values, & 2110 'Body Force', minv=1,maxv=CurrentModel % NumberOfBodyForces ) 2111 END IF 2112!------------------------------------------------------------------------------ 2113 END FUNCTION GetBodyForceId 2114!------------------------------------------------------------------------------ 2115 2116 2117 2118!------------------------------------------------------------------------------ 2119!> Returns the material index of the active element 2120 FUNCTION GetMaterialId( Element, Found ) RESULT(mat_id) 2121!------------------------------------------------------------------------------ 2122 LOGICAL, OPTIONAL :: Found 2123 TYPE(Element_t), OPTIONAL :: Element 2124 TYPE(Element_t), POINTER :: CurrElement 2125 2126 INTEGER :: mat_id, body_id 2127 2128 CurrElement => GetCurrentElement(Element) 2129 body_id = CurrElement % BodyId 2130 2131 IF( body_id <= 0 ) THEN 2132 mat_id = 0 2133 IF( PRESENT( Found ) ) Found = .FALSE. 2134 RETURN 2135 END IF 2136 2137 IF ( PRESENT( Found ) ) THEN 2138 mat_id = ListGetInteger( CurrentModel % Bodies(body_id) % Values, & 2139 'Material', Found, minv=1,maxv=CurrentModel % NumberOfMaterials ) 2140 ELSE 2141 mat_id = ListGetInteger( CurrentModel % Bodies(body_id) % Values, & 2142 'Material', minv=1,maxv=CurrentModel % NumberOfMaterials ) 2143 END IF 2144!------------------------------------------------------------------------------ 2145 END FUNCTION GetMaterialId 2146!------------------------------------------------------------------------------ 2147 2148 2149!------------------------------------------------------------------------------ 2150!> Get component list given component id 2151 FUNCTION GetComponent(i) RESULT(list) 2152!------------------------------------------------------------------------------ 2153 INTEGER :: i 2154 TYPE(ValueList_t), POINTER :: list 2155 2156 List => Null() 2157 IF(i>=0 .AND. i<=SIZE(CurrentModel % Components)) list=> & 2158 CurrentModel % Components(i) % Values 2159!------------------------------------------------------------------------------ 2160 END FUNCTION GetComponent 2161!------------------------------------------------------------------------------ 2162 2163 2164!------------------------------------------------------------------------------ 2165!> Returns the equation index of the active element 2166 FUNCTION GetEquationId( Element, Found ) RESULT(eq_id) 2167!------------------------------------------------------------------------------ 2168 LOGICAL, OPTIONAL :: Found 2169 TYPE(Element_t), OPTIONAL :: Element 2170 TYPE(Element_t), POINTER :: CurrElement 2171 2172 INTEGER :: eq_id, body_id 2173 2174 CurrElement => GetCurrentElement(Element) 2175 body_id = CurrElement % BodyId 2176 2177 IF ( PRESENT( Found ) ) THEN 2178 eq_id = ListGetInteger( CurrentModel % Bodies(body_id) % Values, & 2179 'Equation', Found, minv=1,maxv=CurrentModel % NumberOfEquations ) 2180 ELSE 2181 eq_id = ListGetInteger( CurrentModel % Bodies(body_id) % Values, & 2182 'Equation', minv=1,maxv=CurrentModel % NumberOfEquations ) 2183 END IF 2184!------------------------------------------------------------------------------ 2185 END FUNCTION GetEquationId 2186!------------------------------------------------------------------------------ 2187 2188 2189 2190!------------------------------------------------------------------------------ 2191!> Returns handle to the Simulation value list 2192 FUNCTION GetSimulation() RESULT(Simulation) 2193!------------------------------------------------------------------------------ 2194 TYPE(ValueList_t), POINTER :: Simulation 2195 Simulation => CurrentModel % Simulation 2196!------------------------------------------------------------------------------ 2197 END FUNCTION GetSimulation 2198!------------------------------------------------------------------------------ 2199 2200 2201 2202!------------------------------------------------------------------------------ 2203!> Returns handle to the Constants value list 2204 FUNCTION GetConstants() RESULT(Constants) 2205!------------------------------------------------------------------------------ 2206 TYPE(ValueList_t), POINTER :: Constants 2207 Constants => CurrentModel % Constants 2208!------------------------------------------------------------------------------ 2209 END FUNCTION GetConstants 2210!------------------------------------------------------------------------------ 2211 2212 2213 2214!------------------------------------------------------------------------------ 2215!> Returns handle to the Solver value list of the active solver 2216 FUNCTION GetSolverParams(Solver) RESULT(SolverParam) 2217!------------------------------------------------------------------------------ 2218 TYPE(ValueList_t), POINTER :: SolverParam 2219 TYPE(Solver_t), OPTIONAL :: Solver 2220 2221 SolverParam => ListGetSolverParams(Solver) 2222!------------------------------------------------------------------------------ 2223 END FUNCTION GetSolverParams 2224!------------------------------------------------------------------------------ 2225 2226 2227 2228!------------------------------------------------------------------------------ 2229!> Returns handle to Material value list of the active element 2230 FUNCTION GetMaterial( Element, Found ) RESULT(Material) 2231!------------------------------------------------------------------------------ 2232 TYPE(Element_t), OPTIONAL :: Element 2233 LOGICAL, OPTIONAL :: Found 2234 2235 TYPE(ValueList_t), POINTER :: Material 2236 2237 LOGICAL :: L 2238 INTEGER :: mat_id 2239 2240 IF ( PRESENT( Element ) ) THEN 2241 mat_id = GetMaterialId( Element, L ) 2242 ELSE 2243 mat_id = GetMaterialId( Found=L ) 2244 END IF 2245 2246 Material => Null() 2247 IF ( L ) Material => CurrentModel % Materials(mat_id) % Values 2248 IF ( PRESENT( Found ) ) Found = L 2249!------------------------------------------------------------------------------ 2250 END FUNCTION GetMaterial 2251!------------------------------------------------------------------------------ 2252 2253 2254!------------------------------------------------------------------------------ 2255!> Returns handle to Parent element of a boundary element with a larger body id. 2256!------------------------------------------------------------------------------ 2257 FUNCTION GetBulkElementAtBoundary( Element, Found ) RESULT(BulkElement) 2258!------------------------------------------------------------------------------ 2259 TYPE(Element_t), OPTIONAL :: Element 2260 LOGICAL, OPTIONAL :: Found 2261 TYPE(element_t), POINTER :: BulkElement 2262!------------------------------------------------------------------------------ 2263 TYPE(element_t), POINTER :: BulkElementL, BulkElementR, BoundaryElement 2264 LOGICAL :: L 2265 INTEGER :: mat_id, BodyIdL, BodyIdR 2266 2267 BulkElement => NULL() 2268 2269 BoundaryElement => GetCurrentElement(Element) 2270 2271 IF ( .NOT. ASSOCIATED(BoundaryElement % boundaryinfo)) RETURN 2272 BulkElementR => BoundaryElement % boundaryinfo % right 2273 BulkElementL => BoundaryElement % boundaryinfo % left 2274 BodyIdR = 0; BodyIdL = 0 2275 2276 IF (ASSOCIATED(BulkElementR)) BodyIdR = BulkElementR % BodyId 2277 IF (ASSOCIATED(BulkElementL)) BodyIdL = BulkElementL % BodyId 2278 2279 IF (BodyIdR == 0 .AND. BodyIdL == 0) THEN 2280 RETURN 2281 ELSE IF (BodyIdR > BodyIdL) THEN 2282 BulkElement => BulkElementR 2283 ELSE IF (bodyIdL >= BodyIdR) THEN 2284 BulkElement => BulkElementL 2285 END IF 2286 2287 IF( PRESENT( Found ) ) Found = ASSOCIATED( BulkElement ) 2288 2289!------------------------------------------------------------------------------ 2290 END FUNCTION GetBulkElementAtBoundary 2291!------------------------------------------------------------------------------ 2292 2293 2294!------------------------------------------------------------------------------ 2295!> Returns handle to Material value list of the bulk material meeting 2296!> element with larger body id. Typically Element is a boundary element. 2297 FUNCTION GetBulkMaterialAtBoundary( Element, Found ) RESULT(Material) 2298!------------------------------------------------------------------------------ 2299 TYPE(Element_t), OPTIONAL :: Element 2300 LOGICAL, OPTIONAL :: Found 2301 TYPE(ValueList_t), POINTER :: Material 2302!------------------------------------------------------------------------------ 2303 TYPE(element_t), POINTER :: BulkElement 2304 LOGICAL :: L 2305 INTEGER :: mat_id 2306 2307 Material => NULL() 2308 2309 BulkElement => GetBulkElementAtBoundary(Element, Found) 2310 2311 IF( ASSOCIATED( BulkElement ) ) THEN 2312 mat_id = GetMaterialId( BulkElement, L ) 2313 IF ( L ) Material => CurrentModel % Materials(mat_id) % Values 2314 ELSE 2315 L = .FALSE. 2316 END IF 2317 2318 IF ( PRESENT( Found ) ) Found = L 2319!------------------------------------------------------------------------------ 2320 END FUNCTION GetBulkMaterialAtBoundary 2321!------------------------------------------------------------------------------ 2322 2323!------------------------------------------------------------------------------ 2324!> Return handle to the Body Force value list of the active element 2325 FUNCTION GetBodyForce( Element, Found ) RESULT(BodyForce) 2326!------------------------------------------------------------------------------ 2327 TYPE(Element_t), OPTIONAL :: Element 2328 LOGICAL, OPTIONAL :: Found 2329 2330 TYPE(ValueList_t), POINTER :: BodyForce 2331 2332 LOGICAL :: l 2333 INTEGER :: bf_id 2334 2335 IF ( PRESENT( Element ) ) THEN 2336 bf_id = GetBodyForceId( Element, L ) 2337 ELSE 2338 bf_id = GetBodyForceId( Found=L ) 2339 END IF 2340 2341 BodyForce => Null() 2342 IF ( L ) BodyForce => CurrentModel % BodyForces(bf_id) % Values 2343 IF ( PRESENT( Found ) ) Found = L 2344!------------------------------------------------------------------------------ 2345 END FUNCTION GetBodyForce 2346!------------------------------------------------------------------------------ 2347 2348 2349!> Is the actice solver solved in the frequency space 2350!------------------------------------------------------------------------------ 2351 FUNCTION EigenOrHarmonicAnalysis(Usolver) RESULT(L) 2352 LOGICAL :: L 2353 TYPE(Solver_t), OPTIONAL,TARGET :: USolver 2354!------------------------------------------------------------------------------ 2355 TYPE(Solver_t), POINTER :: Solver 2356 2357 IF (PRESENT(USolver)) THEN 2358 Solver => USolver 2359 ELSE 2360 Solver => CurrentModel % Solver 2361 END IF 2362 L = Solver % NOFEigenValues > 0 2363!------------------------------------------------------------------------------ 2364 END FUNCTION EigenOrHarmonicAnalysis 2365!------------------------------------------------------------------------------ 2366 2367 2368!> Returns the handle to the equation where the active element belongs to 2369!------------------------------------------------------------------------------ 2370 FUNCTION GetEquation( Element, Found ) RESULT(Equation) 2371!------------------------------------------------------------------------------ 2372 TYPE(Element_t), OPTIONAL :: Element 2373 LOGICAL, OPTIONAL :: Found 2374 2375 TYPE(ValueList_t), POINTER :: Equation 2376 2377 LOGICAL :: L 2378 INTEGER :: eq_id 2379 2380 2381 IF ( PRESENT( Element ) ) THEN 2382 eq_id = GetEquationId( Element, L ) 2383 ELSE 2384 eq_id = GetEquationId( Found=L ) 2385 END IF 2386 2387 NULLIFY( Equation ) 2388 IF ( L ) Equation => CurrentModel % Equations(eq_id) % Values 2389 IF ( PRESENT( Found ) ) Found = L 2390!------------------------------------------------------------------------------ 2391 END FUNCTION GetEquation 2392!------------------------------------------------------------------------------ 2393 2394 2395 2396!> Returns the Boundary Condition index of the active element 2397!------------------------------------------------------------------------------ 2398 FUNCTION GetBCId( UElement ) RESULT(bc_id) 2399!------------------------------------------------------------------------------ 2400 TYPE(Element_t), OPTIONAL, TARGET :: UElement 2401 2402 INTEGER :: bc_id 2403 2404 TYPE(Element_t), POINTER :: Element 2405 2406 Element => GetCurrentElement( UElement ) 2407 2408 IF(.NOT. ASSOCIATED( Element % BoundaryInfo ) ) THEN 2409 bc_id = 0 2410 RETURN 2411 END IF 2412 2413 DO bc_id=1,CurrentModel % NumberOfBCs 2414 IF ( Element % BoundaryInfo % Constraint == CurrentModel % BCs(bc_id) % Tag ) EXIT 2415 END DO 2416 IF ( bc_id > CurrentModel % NumberOfBCs ) bc_id=0 2417!------------------------------------------------------------------------------ 2418 END FUNCTION GetBCId 2419!------------------------------------------------------------------------------ 2420 2421 2422!> Returns handle to the value list of the Boundary Condition where the active element belongs to 2423!------------------------------------------------------------------------------ 2424 FUNCTION GetBC( UElement ) RESULT(bc) 2425!------------------------------------------------------------------------------ 2426 TYPE(Element_t), OPTIONAL, TARGET :: UElement 2427 TYPE(ValueList_t), POINTER :: BC 2428 2429 INTEGER :: bc_id 2430 2431 TYPE(Element_t), POINTER :: Element 2432 2433 Element => GetCurrentElement( UElement ) 2434 2435 BC => Null() 2436 bc_id = GetBCId( Element ) 2437 IF ( bc_id > 0 ) BC => CurrentModel % BCs(bc_id) % Values 2438!------------------------------------------------------------------------------ 2439 END FUNCTION GetBC 2440!------------------------------------------------------------------------------ 2441 2442 2443!> Returns the index of the Initial Condition of the active element 2444!------------------------------------------------------------------------------ 2445 FUNCTION GetICId( Element, Found ) RESULT(ic_id) 2446!------------------------------------------------------------------------------ 2447 LOGICAL, OPTIONAL :: Found 2448 TYPE(Element_t), OPTIONAL :: Element 2449 2450 TYPE(Element_t), POINTER :: CElement 2451 INTEGER :: ic_id, body_id 2452 2453 CElement => GetCurrentElement( Element ) 2454 body_id = CElement % BodyId 2455 2456 IF ( PRESENT( Found ) ) THEN 2457 ic_id = ListGetInteger( CurrentModel % Bodies(body_id) % Values, & 2458 'Initial Condition', Found, minv=1,maxv=CurrentModel % NumberOfICs ) 2459 ELSE 2460 ic_id = ListGetInteger( CurrentModel % Bodies(body_id) % Values, & 2461 'Initial Condition', minv=1,maxv=CurrentModel % NumberOfICs ) 2462 END IF 2463!------------------------------------------------------------------------------ 2464 END FUNCTION GetIcId 2465!------------------------------------------------------------------------------ 2466 2467!> Returns handle to the value list of the Initial Condition where the active element belongs to 2468!------------------------------------------------------------------------------ 2469 FUNCTION GetIC( Element, Found ) RESULT(IC) 2470!------------------------------------------------------------------------------ 2471 TYPE(Element_t), OPTIONAL :: Element 2472 LOGICAL, OPTIONAL :: Found 2473 2474 TYPE(ValueList_t), POINTER :: IC 2475 2476 LOGICAL :: L 2477 INTEGER :: ic_id 2478 2479 IF ( PRESENT( Element ) ) THEN 2480 ic_id = GetICId( Element, L ) 2481 ELSE 2482 ic_id = GetICId( Found=L ) 2483 END IF 2484 2485 IC => Null() 2486 IF ( L ) IC => CurrentModel % ICs(ic_id) % Values 2487 IF ( PRESENT( Found ) ) Found = L 2488!------------------------------------------------------------------------------ 2489 END FUNCTION GetIC 2490!------------------------------------------------------------------------------ 2491 2492!> Add the local matrix entries to for real valued equations that are of first order in time 2493!------------------------------------------------------------------------------ 2494 SUBROUTINE Default1stOrderTimeR( M, A, F, UElement, USolver ) 2495!------------------------------------------------------------------------------ 2496 REAL(KIND=dp) :: M(:,:),A(:,:), F(:) 2497 TYPE(Solver_t), OPTIONAL, TARGET :: USolver 2498 TYPE(Element_t), OPTIONAL, TARGET :: UElement 2499 2500 LOGICAL :: Found 2501 TYPE(ValueList_t), POINTER :: Params 2502 2503 TYPE(Solver_t), POINTER :: Solver 2504 TYPE(Variable_t), POINTER :: x 2505 TYPE(Element_t), POINTER :: Element 2506 2507 INTEGER :: n 2508 REAL(KIND=dp) :: dt 2509 INTEGER, POINTER :: Indexes(:) 2510 2511 IF ( PRESENT(USolver) ) THEN 2512 Solver => USolver 2513 ELSE 2514 Solver => CurrentModel % Solver 2515 END IF 2516 2517 Params => GetSolverParams(Solver) 2518 2519 ! Antiperiodic elimination and FCT always use this 2520 IF (GetLogical(Params,'Use Global Mass Matrix',Found)) THEN 2521 CALL DefaultUpdateMass(M,UElement,USolver) 2522 RETURN 2523 END IF 2524 2525 Element => GetCurrentElement( UElement ) 2526 2527 x => Solver % Variable 2528 2529 dt = Solver % dt 2530 Indexes => GetIndexStore() 2531 n = GetElementDOFs( Indexes,Element,Solver ) 2532 2533 CALL Add1stOrderTime( M, A, F, dt, n, x % DOFs, & 2534 x % Perm(Indexes(1:n)), Solver, UElement=Element ) 2535 2536!------------------------------------------------------------------------------ 2537 END SUBROUTINE Default1stOrderTimeR 2538!------------------------------------------------------------------------------ 2539 2540!> Add the local matrix entries to for complex valued equations that are of first order in time 2541!------------------------------------------------------------------------------ 2542 SUBROUTINE Default1stOrderTimeC( MC, AC, FC, UElement, USolver ) 2543!------------------------------------------------------------------------------ 2544 COMPLEX(KIND=dp) :: MC(:,:),AC(:,:), FC(:) 2545 TYPE(Solver_t), OPTIONAL, TARGET :: USolver 2546 TYPE(Element_t), OPTIONAL, TARGET :: UElement 2547 2548 TYPE(Solver_t), POINTER :: Solver 2549 TYPE(Variable_t), POINTER :: x 2550 TYPE(Element_t), POINTER :: Element 2551 2552 REAL(KIND=dp), ALLOCATABLE :: M(:,:),A(:,:), F(:) 2553 2554 LOGICAL :: Found 2555 TYPE(ValueList_t), POINTER :: Params 2556 2557 INTEGER :: i,j,n,DOFs 2558 REAL(KIND=dp) :: dt 2559 INTEGER, POINTER :: Indexes(:) 2560 2561 IF ( PRESENT(USolver) ) THEN 2562 Solver => USolver 2563 ELSE 2564 Solver => CurrentModel % Solver 2565 END IF 2566 2567 Params=>GetSolverParams(Solver) 2568 2569 IF (GetLogical(Params,'Use Global Mass Matrix',Found)) THEN 2570 CALL DefaultUpdateMass(M,UElement,USolver) 2571 RETURN 2572 END IF 2573 2574 Element => GetCurrentElement( UElement ) 2575 2576 x => Solver % Variable 2577 2578 dt = Solver % dt 2579 DOFs = x % DOFs 2580 Indexes => GetIndexStore() 2581 n = GetElementDOFs( Indexes,Element,Solver ) 2582 2583 ALLOCATE( M(DOFs*n,DOFs*n), A(DOFs*n,DOFs*n), F(DOFs*n) ) 2584 DO i=1,n*DOFs/2 2585 F( 2*(i-1)+1 ) = REAL( FC(i) ) 2586 F( 2*(i-1)+2 ) = AIMAG( FC(i) ) 2587 2588 DO j=1,n*DOFs/2 2589 M( 2*(i-1)+1, 2*(j-1)+1 ) = REAL( MC(i,j) ) 2590 M( 2*(i-1)+1, 2*(j-1)+2 ) = -AIMAG( MC(i,j) ) 2591 M( 2*(i-1)+2, 2*(j-1)+1 ) = AIMAG( MC(i,j) ) 2592 M( 2*(i-1)+2, 2*(j-1)+2 ) = REAL( MC(i,j) ) 2593 A( 2*(i-1)+1, 2*(j-1)+1 ) = REAL( AC(i,j) ) 2594 A( 2*(i-1)+1, 2*(j-1)+2 ) = -AIMAG( AC(i,j) ) 2595 A( 2*(i-1)+2, 2*(j-1)+1 ) = AIMAG( AC(i,j) ) 2596 A( 2*(i-1)+2, 2*(j-1)+2 ) = REAL( AC(i,j) ) 2597 END DO 2598 END DO 2599 2600 CALL Add1stOrderTime( M, A, F, dt, n, x % DOFs, & 2601 x % Perm(Indexes(1:n)), Solver, UElement=Element ) 2602 2603 DO i=1,n*DOFs/2 2604 FC(i) = CMPLX( F(2*(i-1)+1), F(2*(i-1)+2),KIND=dp ) 2605 DO j=1,n*DOFs/2 2606 MC(i,j) = CMPLX(M(2*(i-1)+1,2*(j-1)+1), -M(2*(i-1)+1,2*(j-1)+2), KIND=dp) 2607 AC(i,j) = CMPLX(A(2*(i-1)+1,2*(j-1)+1), -A(2*(i-1)+1,2*(j-1)+2), KIND=dp) 2608 END DO 2609 END DO 2610 2611 DEALLOCATE( M, A, F ) 2612!------------------------------------------------------------------------------ 2613 END SUBROUTINE Default1stOrderTimeC 2614!------------------------------------------------------------------------------ 2615 2616 2617!------------------------------------------------------------------------------ 2618 SUBROUTINE Default1stOrderTimeGlobal(USolver) 2619!------------------------------------------------------------------------------ 2620 TYPE(Solver_t), OPTIONAL, TARGET :: USolver 2621!------------------------------------------------------------------------------ 2622 CHARACTER(LEN=MAX_NAME_LEN) :: Method 2623 TYPE(Solver_t), POINTER :: Solver 2624 INTEGER :: i,j,k,l,n,Order 2625 REAL(KIND=dp), POINTER :: SaveValues(:) => NULL() 2626 REAL(KIND=dp) :: FORCE(1), Dts(16) 2627 LOGICAL :: ConstantDt, Found, HasMass, HasFCT 2628 TYPE(Variable_t), POINTER :: DtVar 2629 SAVE STIFF, MASS, X 2630 REAL(KIND=dp), ALLOCATABLE :: STIFF(:,:),MASS(:,:), X(:,:) 2631 2632 !$OMP THREADPRIVATE(SaveValues) 2633 2634 IF ( PRESENT(USolver) ) THEN 2635 Solver => Usolver 2636 ELSE 2637 Solver => CurrentModel % solver 2638 END IF 2639 2640 Order = MAX( MIN( Solver % DoneTime, Solver % Order ), 1) 2641 HasMass = ASSOCIATED( Solver % Matrix % MassValues ) 2642 2643 HasFCT = ListGetLogical( Solver % Values,'Linear System FCT', Found ) 2644 2645 IF( HasFCT ) THEN 2646 IF( .NOT. HasMass ) THEN 2647 CALL Fatal('Default1stOrderTimeGlobal','FCT only makes sense if there is a mass matrix!') 2648 ELSE 2649 IF(.NOT. ASSOCIATED( Solver % Matrix % MassValuesLumped ) ) THEN 2650 CALL Fatal('Default1stOrderTimeGlobal','FCT requires a lumped mass matrix!') 2651 END IF 2652 HasMass = .FALSE. 2653 END IF 2654 END IF 2655 2656 ! This is now the default global time integration routine but the old hack may still be called 2657 !--------------------------------------------------------------------------------------------- 2658 IF( .NOT. ListGetLogical( Solver % Values,'Old Global Time Integration',Found ) ) THEN 2659 CALL Add1stOrderTime_CRS( Solver % Matrix, Solver % Matrix % rhs, & 2660 Solver % dt, Solver ) 2661 RETURN 2662 END IF 2663 2664 2665 ! The rest of the code in this subroutine is obsolete 2666 IF ( .NOT.ASSOCIATED(Solver % Variable % Values, SaveValues) ) THEN 2667 IF ( ALLOCATED(STIFF) ) DEALLOCATE( STIFF,MASS,X ) 2668 n = 0 2669 DO i=1,Solver % Matrix % NumberOfRows 2670 n = MAX( n,Solver % Matrix % Rows(i+1)-Solver % Matrix % Rows(i) ) 2671 END DO 2672 k = SIZE(Solver % Variable % PrevValues,2) 2673 ALLOCATE( STIFF(1,n),MASS(1,n),X(n,k) ) 2674 SaveValues => Solver % Variable % Values 2675 END IF 2676 2677 STIFF = 0.0_dp 2678 MASS = 0.0_dp 2679 X = 0.0_dp 2680 2681 Method = ListGetString( Solver % Values, 'Timestepping Method', Found ) 2682 IF ( Method == 'bdf' ) THEN 2683 Dts(1) = Solver % Dt 2684 ConstantDt = .TRUE. 2685 IF(Order > 1) THEN 2686 DtVar => VariableGet( Solver % Mesh % Variables, 'Timestep size' ) 2687 DO i=2,Order 2688 Dts(i) = DtVar % PrevValues(1,i-1) 2689 IF(ABS(Dts(i)-Dts(1)) > 1.0d-6 * Dts(1)) ConstantDt = .FALSE. 2690 END DO 2691 END IF 2692 END IF 2693 2694 DO i=1,Solver % Matrix % NumberOFRows 2695 n = 0 2696 k = 0 2697 2698 DO j=Solver % Matrix % Rows(i),Solver % Matrix % Rows(i+1)-1 2699 n = n+1 2700 STIFF(1,n) = Solver % Matrix % Values(j) 2701 IF( HasMass ) THEN 2702 MASS(1,n) = Solver % Matrix % MassValues(j) 2703 ELSE IF( HasFCT ) THEN 2704 IF( j == Solver % Matrix % Diag(i) ) k = n 2705 END IF 2706 X(n,:) = Solver % Variable % PrevValues(Solver % Matrix % Cols(j),:) 2707 END DO 2708 2709 ! Use lumped mass in lower order fct 2710 IF( HasFCT ) THEN 2711 IF( k == 0 ) THEN 2712 CALL Fatal('Default1stOrderTimeGlobal','Could not find diagonal entry for fct') 2713 ELSE 2714 MASS(1,k) = Solver % Matrix % MassValuesLumped(i) 2715 END IF 2716 END IF 2717 2718 FORCE(1) = Solver % Matrix % RHS(i) 2719 Solver % Matrix % Force(i,1) = FORCE(1) 2720 2721 SELECT CASE( Method ) 2722 CASE( 'fs' ) 2723 CALL FractionalStep( n, Solver % dt, MASS, STIFF, FORCE, & 2724 X(:,1), Solver % Beta, Solver ) 2725 2726 CASE('bdf') 2727 IF(ConstantDt) THEN 2728 CALL BDFLocal( n, Solver % dt, MASS, STIFF, FORCE, X, Order ) 2729 ELSE 2730 CALL VBDFLocal(n, Dts, MASS, STIFF, FORCE, X, Order ) 2731 END IF 2732 2733 CASE DEFAULT 2734 CALL NewmarkBeta( n, Solver % dt, MASS, STIFF, FORCE, & 2735 X(:,1), Solver % Beta ) 2736 END SELECT 2737 2738 IF( HasFCT ) MASS(1,k) = 0.0_dp 2739 2740 n = 0 2741 DO j=Solver % Matrix % Rows(i),Solver % Matrix % Rows(i+1)-1 2742 n=n+1 2743 Solver % Matrix % Values(j) = STIFF(1,n) 2744 END DO 2745 Solver % Matrix % RHS(i) = FORCE(1) 2746 END DO 2747 2748!---------------------------------------------------------------------------- 2749 END SUBROUTINE Default1stOrderTimeGlobal 2750!---------------------------------------------------------------------------- 2751 2752!------------------------------------------------------------------------------ 2753 SUBROUTINE Default2ndOrderTimeGlobal(USolver) 2754!------------------------------------------------------------------------------ 2755 TYPE(Solver_t), OPTIONAL, TARGET :: USolver 2756!------------------------------------------------------------------------------ 2757 TYPE(Solver_t), POINTER :: Solver 2758 INTEGER :: i,j,k,l,n 2759 REAL(KIND=dp), POINTER :: SaveValues(:) => NULL() 2760 REAL(KIND=dp) :: FORCE(1) 2761 LOGICAL :: Found, HasDamping, HasMass 2762 REAL(KIND=dp), ALLOCATABLE, SAVE :: STIFF(:,:),MASS(:,:), DAMP(:,:), X(:,:) 2763 !OMP THREADPRIVATE(SaveValues) 2764 2765 IF ( PRESENT(USolver) ) THEN 2766 Solver => Usolver 2767 ELSE 2768 Solver => CurrentModel % solver 2769 END IF 2770 2771 IF ( .NOT.ASSOCIATED(Solver % Variable % Values, SaveValues) ) THEN 2772 IF ( ALLOCATED(STIFF) ) DEALLOCATE( STIFF,MASS,DAMP,X ) 2773 n = 0 2774 DO i=1,Solver % Matrix % NumberOfRows 2775 n = MAX( n,Solver % Matrix % Rows(i+1)-Solver % Matrix % Rows(i) ) 2776 END DO 2777 k = SIZE(Solver % Variable % PrevValues,2) 2778 ALLOCATE( STIFF(1,n),MASS(1,n),DAMP(1,n),X(n,k) ) 2779 SaveValues => Solver % Variable % Values 2780 2781 STIFF = 0.0_dp 2782 MASS = 0.0_dp 2783 DAMP = 0.0_dp 2784 X = 0.0_dp 2785 END IF 2786 2787 HasDamping = ASSOCIATED(Solver % Matrix % DampValues ) 2788 HasMass = ASSOCIATED(Solver % Matrix % MassValues ) 2789 2790 DO i=1,Solver % Matrix % NumberOFRows 2791 n = 0 2792 DO j=Solver % Matrix % Rows(i),Solver % Matrix % Rows(i+1)-1 2793 n=n+1 2794 IF( HasMass ) MASS(1,n) = Solver % Matrix % MassValues(j) 2795 IF( HasDamping ) DAMP(1,n) = Solver % Matrix % DampValues(j) 2796 STIFF(1,n) = Solver % Matrix % Values(j) 2797 X(n,:) = Solver % Variable % PrevValues(Solver % Matrix % Cols(j),:) 2798 END DO 2799 FORCE(1) = Solver % Matrix % RHS(i) 2800 Solver % Matrix % Force(i,1) = FORCE(1) 2801 2802 CALL Bossak2ndOrder( n, Solver % dt, MASS, DAMP, STIFF, & 2803 FORCE, X(1:n,3), X(1:n,4), X(1:n,5), Solver % Alpha ) 2804 2805 n = 0 2806 DO j=Solver % Matrix % Rows(i),Solver % Matrix % Rows(i+1)-1 2807 n=n+1 2808 Solver % Matrix % Values(j) = STIFF(1,n) 2809 END DO 2810 Solver % Matrix % RHS(i) = FORCE(1) 2811 END DO 2812!---------------------------------------------------------------------------- 2813 END SUBROUTINE Default2ndOrderTimeGlobal 2814!---------------------------------------------------------------------------- 2815 2816 2817!------------------------------------------------------------------------------ 2818 SUBROUTINE Default2ndOrderTimeR( M, B, A, F, UElement, USolver ) 2819!------------------------------------------------------------------------------ 2820 REAL(KIND=dp) :: M(:,:), B(:,:), A(:,:), F(:) 2821 TYPE(Solver_t), OPTIONAL, TARGET :: USolver 2822 TYPE(Element_t), OPTIONAL, TARGET :: UElement 2823 2824 TYPE(Solver_t), POINTER :: Solver 2825 TYPE(Variable_t), POINTER :: x 2826 TYPE(Element_t), POINTER :: Element 2827 2828 LOGICAL :: Found 2829 TYPE(ValueList_t), POINTER :: Params 2830 2831 INTEGER :: n 2832 REAL(KIND=dp) :: dt 2833 INTEGER, POINTER :: Indexes(:) 2834 2835 Solver => CurrentModel % Solver 2836 IF ( PRESENT(USolver) ) Solver => USolver 2837 2838 Params=>GetSolverParams(Solver) 2839 2840 IF (GetLogical(Params,'Use Global Mass Matrix',Found)) THEN 2841 CALL DefaultUpdateMass(M,UElement,USolver) 2842 CALL DefaultUpdateDamp(B,UElement,USolver) 2843 RETURN 2844 END IF 2845 2846 Element => GetCurrentElement( UElement ) 2847 2848 x => Solver % Variable 2849 2850 dt = Solver % dt 2851 Indexes => GetIndexStore() 2852 n = GetElementDOFs( Indexes, Element, Solver ) 2853 2854 CALL Add2ndOrderTime( M, B, A, F, dt, n, x % DOFs, & 2855 x % Perm(Indexes(1:n)), Solver ) 2856!------------------------------------------------------------------------------ 2857 END SUBROUTINE Default2ndOrderTimeR 2858!------------------------------------------------------------------------------ 2859 2860 2861 2862!------------------------------------------------------------------------------ 2863 SUBROUTINE Default2ndOrderTimeC( MC, BC, AC, FC, UElement, USolver ) 2864!------------------------------------------------------------------------------ 2865 COMPLEX(KIND=dp) :: MC(:,:), BC(:,:), AC(:,:), FC(:) 2866 TYPE(Solver_t), OPTIONAL, TARGET :: USolver 2867 TYPE(Element_t), OPTIONAL, TARGET :: UElement 2868 2869 TYPE(Solver_t), POINTER :: Solver 2870 TYPE(Variable_t), POINTER :: x 2871 TYPE(Element_t), POINTER :: Element 2872 REAL(KIND=dp), ALLOCATABLE :: M(:,:), B(:,:), A(:,:), F(:) 2873 2874 LOGICAL :: Found 2875 TYPE(ValueList_t), POINTER :: Params 2876 2877 INTEGER :: i,j,n,DOFs 2878 REAL(KIND=dp) :: dt 2879 INTEGER, POINTER :: Indexes(:) 2880 2881 Solver => CurrentModel % Solver 2882 IF ( PRESENT(USolver) ) Solver => USolver 2883 2884 Params=>GetSolverParams(Solver) 2885 2886 IF (GetLogical(Params,'Use Global Mass Matrix',Found)) THEN 2887 CALL DefaultUpdateMass(M,UElement,USolver) 2888 CALL DefaultUpdateDamp(B,UElement,USolver) 2889 RETURN 2890 END IF 2891 2892 Element => GetCurrentElement( UElement ) 2893 2894 x => Solver % Variable 2895 2896 dt = Solver % dt 2897 DOFs = x % DOFs 2898 Indexes => GetIndexStore() 2899 n = GetElementDOFs( Indexes, Element, Solver ) 2900 2901 ALLOCATE( M(DOFs*n,DOFs*n), A(DOFs*n,DOFs*n), B(DOFs*n,DOFs*n), F(DOFs*n) ) 2902 DO i=1,n*DOFs/2 2903 F( 2*(i-1)+1 ) = REAL( FC(i) ) 2904 F( 2*(i-1)+2 ) = AIMAG( FC(i) ) 2905 2906 DO j=1,n*DOFs/2 2907 M(2*(i-1)+1, 2*(j-1)+1) = REAL( MC(i,j) ) 2908 M(2*(i-1)+1, 2*(j-1)+2) = -AIMAG( MC(i,j) ) 2909 M(2*(i-1)+2, 2*(j-1)+1) = AIMAG( MC(i,j) ) 2910 M(2*(i-1)+2, 2*(j-1)+2) = REAL( MC(i,j) ) 2911 B(2*(i-1)+1, 2*(j-1)+1) = REAL( BC(i,j) ) 2912 B(2*(i-1)+1, 2*(j-1)+2) = -AIMAG( BC(i,j) ) 2913 B(2*(i-1)+2, 2*(j-1)+1) = AIMAG( BC(i,j) ) 2914 B(2*(i-1)+2, 2*(j-1)+2) = REAL( BC(i,j) ) 2915 A(2*(i-1)+1, 2*(j-1)+1) = REAL( AC(i,j) ) 2916 A(2*(i-1)+1, 2*(j-1)+2) = -AIMAG( AC(i,j) ) 2917 A(2*(i-1)+2, 2*(j-1)+1) = AIMAG( AC(i,j) ) 2918 A(2*(i-1)+2, 2*(j-1)+2) = REAL( AC(i,j) ) 2919 END DO 2920 END DO 2921 2922 CALL Add2ndOrderTime( M, B, A, F, dt, n, x % DOFs, & 2923 x % Perm(Indexes(1:n)), Solver ) 2924 2925 DO i=1,n*DOFs/2 2926 FC(i) = CMPLX( F(2*(i-1)+1), F(2*(i-1)+2), KIND=dp ) 2927 DO j=1,n*DOFs/2 2928 MC(i,j) = CMPLX( M(2*(i-1)+1, 2*(j-1)+1), -M(2*(i-1)+1, 2*(j-1)+2), KIND=dp ) 2929 BC(i,j) = CMPLX( B(2*(i-1)+1, 2*(j-1)+1), -B(2*(i-1)+1, 2*(j-1)+2), KIND=dp ) 2930 AC(i,j) = CMPLX( A(2*(i-1)+1, 2*(j-1)+1), -A(2*(i-1)+1, 2*(j-1)+2), KIND=dp ) 2931 END DO 2932 END DO 2933 2934 DEALLOCATE( M, B, A, F ) 2935!------------------------------------------------------------------------------ 2936 END SUBROUTINE Default2ndOrderTimeC 2937!------------------------------------------------------------------------------ 2938 2939 2940!-------------------------------------------------------------------------------- 2941!> One can enforce weak coupling by calling a dependent solver a.k.a. slave solver 2942!> at different stages of the master solver: e.g. before and after the solver. 2943!> The strategy can be particularly efficient for nonlinear problems when the 2944!> slave solver is cheap and a stepsize control is applied 2945!> Also one can easily make postprocessing steps just at the correct slot. 2946!----------------------------------------------------------------------------- 2947 RECURSIVE SUBROUTINE DefaultSlaveSolvers( Solver, SlaveSolverStr) 2948!------------------------------------------------------------------------------ 2949 TYPE(Solver_t), POINTER :: Solver 2950 CHARACTER(LEN=*) :: SlaveSolverStr 2951 2952 TYPE(Solver_t), POINTER :: SlaveSolver 2953 TYPE(ValueList_t), POINTER :: Params 2954 TYPE(Variable_t), POINTER :: iterV 2955 INTEGER, POINTER :: SlaveSolverIndexes(:) 2956 INTEGER :: j,k,iter 2957 REAL(KIND=dp) :: dt 2958 LOGICAL :: Transient, Found, alloc_parenv 2959 2960 TYPE(ParEnv_t) :: SParEnv 2961 2962 INTERFACE 2963 SUBROUTINE SolverActivate_x(Model,Solver,dt,Transient) 2964 USE Types 2965 TYPE(Model_t)::Model 2966 TYPE(Solver_t),POINTER::Solver 2967 REAL(KIND=dp) :: dt 2968 LOGICAL :: Transient 2969 END SUBROUTINE SolverActivate_x 2970 END INTERFACE 2971 2972 SlaveSolverIndexes => ListGetIntegerArray( Solver % Values,& 2973 SlaveSolverStr,Found ) 2974 IF(.NOT. Found ) RETURN 2975 2976 CALL Info('DefaultSlaveSolvers','Executing slave solvers: '// & 2977 TRIM(SlaveSolverStr),Level=5) 2978 2979 dt = GetTimeStep() 2980 Transient = GetString(CurrentModel % Simulation,'Simulation type',Found)=='transient' 2981 2982 ! store the nonlinear iteration at the outer loop 2983 iterV => VariableGet( Solver % Mesh % Variables, 'nonlin iter' ) 2984 iter = NINT(iterV % Values(1)) 2985 2986 DO j=1,SIZE(SlaveSolverIndexes) 2987 k = SlaveSolverIndexes(j) 2988 SlaveSolver => CurrentModel % Solvers(k) 2989 2990 CALL Info('DefaultSlaveSolvers','Calling slave solver: '//TRIM(I2S(k)),Level=8) 2991 2992 IF(ParEnv % PEs>1) THEN 2993 SParEnv = ParEnv 2994 2995 IF(ASSOCIATED(SlaveSolver % Matrix)) THEN 2996 IF(ASSOCIATED(SlaveSolver % Matrix % ParMatrix) ) THEN 2997 ParEnv = SlaveSolver % Matrix % ParMatrix % ParEnv 2998 ELSE 2999 ParEnv % ActiveComm = SlaveSolver % Matrix % Comm 3000 END IF 3001 ELSE 3002 CALL ListAddLogical( SlaveSolver % Values, 'Slave not parallel', .TRUE.) 3003 END IF 3004 END IF 3005 3006 CurrentModel % Solver => SlaveSolver 3007 CALL SolverActivate_x( CurrentModel,SlaveSolver,dt,Transient) 3008 3009 IF(ParEnv % PEs>1) THEN 3010 ParEnv = SParEnv 3011 END IF 3012 END DO 3013 iterV % Values = iter 3014 CurrentModel % Solver => Solver 3015 3016 END SUBROUTINE DefaultSlaveSolvers 3017!------------------------------------------------------------------------------ 3018 3019 3020 3021!> Performs initialization for matrix equation related to the active solver 3022!------------------------------------------------------------------------------ 3023 RECURSIVE SUBROUTINE DefaultInitialize( USolver, UseConstantBulk ) 3024!------------------------------------------------------------------------------ 3025 TYPE(Solver_t), OPTIONAL, TARGET, INTENT(IN) :: USolver 3026 LOGICAL, OPTIONAL :: UseConstantBulk 3027!------------------------------------------------------------------------------ 3028 TYPE(Solver_t), POINTER :: Solver 3029 INTEGER :: i,n 3030 LOGICAL :: Found 3031 3032 IF ( PRESENT( USolver ) ) THEN 3033 Solver => USolver 3034 ELSE 3035 Solver => CurrentModel % Solver 3036 END IF 3037 3038 IF( PRESENT( UseConstantBulk ) ) THEN 3039 IF ( UseConstantBulk ) THEN 3040 CALL Info('DefaultInitialize','Using constant bulk matrix',Level=8) 3041 IF( .NOT. ASSOCIATED( Solver % Matrix % BulkValues ) ) THEN 3042 CALL Warn('DefaultInitialize','Constant bulk system requested but not associated!') 3043 RETURN 3044 END IF 3045 3046 n = SIZE(Solver % Matrix % Values) 3047 DO i=1,n 3048 Solver % Matrix % Values(i) = Solver % Matrix % BulkValues(i) 3049 END DO 3050 3051 IF( ASSOCIATED( Solver % Matrix % BulkMassValues ) ) THEN 3052 DO i=1,n 3053 Solver % Matrix % MassValues(i) = Solver % Matrix % BulkMassValues(i) 3054 END DO 3055 END IF 3056 IF( ASSOCIATED( Solver % Matrix % BulkDampValues ) ) THEN 3057 DO i=1,n 3058 Solver % Matrix % DampValues(i) = Solver % Matrix % BulkDampValues(i) 3059 END DO 3060 END IF 3061 IF( ASSOCIATED( Solver % Matrix % BulkRhs ) ) THEN 3062 n = SIZE(Solver % Matrix % RHS) 3063 DO i=1,n 3064 Solver % Matrix % rhs(i) = Solver % Matrix % BulkRhs(i) 3065 END DO 3066 END IF 3067 RETURN 3068 END IF 3069 END IF 3070 3071 3072 CALL DefaultSlaveSolvers(Solver,'Slave Solvers') ! this is the initial name of the slot 3073 CALL DefaultSlaveSolvers(Solver,'Nonlinear Pre Solvers') 3074 3075 3076 ! If we changed the system last time to harmonic one then revert back the real system 3077 IF( ListGetLogical( Solver % Values,'Harmonic Mode',Found ) ) THEN 3078 CALL ChangeToHarmonicSystem( Solver, .TRUE. ) 3079 END IF 3080 3081 3082 IF(.NOT. ASSOCIATED( Solver % Matrix ) ) THEN 3083 CALL Fatal('DefaultInitialize','No matrix exists, cannot initialize!') 3084 END IF 3085 3086 CALL InitializeToZero( Solver % Matrix, Solver % Matrix % RHS ) 3087 3088 IF( ALLOCATED(Solver % Matrix % ConstrainedDOF) ) THEN 3089 Solver % Matrix % ConstrainedDOF = .FALSE. 3090 END IF 3091 3092 IF( ALLOCATED(Solver % Matrix % Dvalues) ) THEN 3093 Solver % Matrix % Dvalues = 0._dp 3094 END IF 3095 3096 IF( ListGetLogical( Solver % Values,'Bulk Assembly Timing',Found ) ) THEN 3097 CALL ResetTimer('BulkAssembly'//GetVarName(Solver % Variable) ) 3098 END IF 3099 3100 ! This is a slot for calling solver that contribute to the assembly 3101 CALL DefaultSlaveSolvers(Solver,'Assembly Solvers') 3102 3103!------------------------------------------------------------------------------ 3104 END SUBROUTINE DefaultInitialize 3105!------------------------------------------------------------------------------ 3106 3107 3108 3109!> Performs pre-steps related to the the active solver 3110!------------------------------------------------------------------------------ 3111 RECURSIVE SUBROUTINE DefaultStart( USolver ) 3112!------------------------------------------------------------------------------ 3113 TYPE(Solver_t), OPTIONAL, TARGET, INTENT(IN) :: USolver 3114 3115 TYPE(Solver_t), POINTER :: Solver 3116 LOGICAL :: Found 3117 TYPE(ValueList_t), POINTER :: Params 3118 3119 IF ( PRESENT( USolver ) ) THEN 3120 Solver => USolver 3121 ELSE 3122 Solver => CurrentModel % Solver 3123 END IF 3124 3125 Params => Solver % Values 3126 3127 CALL Info('DefaultStart','Starting solver: '//& 3128 TRIM(ListGetString(Params,'Equation')),Level=10) 3129 3130 ! When Newton linearization is used we may reset it after previously visiting the solver 3131 IF( Solver % NewtonActive ) THEN 3132 IF( ListGetLogical( Params,'Nonlinear System Reset Newton', Found) ) Solver % NewtonActive = .FALSE. 3133 END IF 3134 3135 ! If we changed the system last time to harmonic one then revert back the real system 3136 IF( ListGetLogical( Params,'Harmonic Mode',Found ) ) THEN 3137 CALL ChangeToHarmonicSystem( Solver, .TRUE. ) 3138 END IF 3139 3140 ! One can run preprocessing solver in this slot. 3141 !----------------------------------------------------------------------------- 3142 CALL DefaultSlaveSolvers(Solver,'Pre Solvers') 3143 3144!------------------------------------------------------------------------------ 3145 END SUBROUTINE DefaultStart 3146!------------------------------------------------------------------------------ 3147 3148 3149 3150!> Performs finalizing steps related to the the active solver 3151!------------------------------------------------------------------------------ 3152 RECURSIVE SUBROUTINE DefaultFinish( USolver ) 3153!------------------------------------------------------------------------------ 3154 TYPE(Solver_t), OPTIONAL, TARGET, INTENT(IN) :: USolver 3155 3156 TYPE(Solver_t), POINTER :: Solver 3157 3158 IF ( PRESENT( USolver ) ) THEN 3159 Solver => USolver 3160 ELSE 3161 Solver => CurrentModel % Solver 3162 END IF 3163 3164 ! One can run postprocessing solver in this slot. 3165 !----------------------------------------------------------------------------- 3166 CALL DefaultSlaveSolvers(Solver,'Post Solvers') 3167 3168 CALL Info('DefaultFinish','Finished solver: '//& 3169 TRIM(ListGetString(Solver % Values,'Equation')),Level=8) 3170 3171!------------------------------------------------------------------------------ 3172 END SUBROUTINE DefaultFinish 3173!------------------------------------------------------------------------------ 3174 3175 3176!> Solver the matrix equation related to the active solver 3177!------------------------------------------------------------------------------ 3178 RECURSIVE FUNCTION DefaultSolve( USolver, BackRotNT ) RESULT(Norm) 3179!------------------------------------------------------------------------------ 3180 TYPE(Solver_t), OPTIONAL, TARGET, INTENT(in) :: USolver 3181 REAL(KIND=dp) :: Norm 3182 LOGICAL, OPTIONAL, INTENT(in) :: BackRotNT 3183 3184 TYPE(Matrix_t), POINTER :: A 3185 TYPE(Variable_t), POINTER :: x 3186 REAL(KIND=dp), POINTER CONTIG :: b(:) 3187 REAL(KIND=dp), POINTER CONTIG :: SOL(:) 3188! REAL(KIND=dp), POINTER :: SOL(:) 3189 3190 LOGICAL :: Found, BackRot 3191 3192 TYPE(ValueList_t), POINTER :: Params 3193 TYPE(Solver_t), POINTER :: Solver 3194 TYPE(Matrix_t), POINTER :: Ctmp 3195 CHARACTER(LEN=MAX_NAME_LEN) :: linsolver, precond, dumpfile, saveslot 3196 INTEGER :: NameSpaceI, Count, MaxCount, i 3197 LOGICAL :: LinearSystemTrialing, SourceControl 3198 3199 CALL Info('DefaultSolve','Solving linear system with default routines',Level=10) 3200 3201 Solver => CurrentModel % Solver 3202 Norm = REAL(0, dp) 3203 IF ( PRESENT( USolver ) ) Solver => USolver 3204 3205 Params => GetSolverParams(Solver) 3206 3207 NameSpaceI = NINT( ListGetCReal( Params,'Linear System Namespace Number', Found ) ) 3208 LinearSystemTrialing = ListGetLogical( Params,'Linear System Trialing', Found ) 3209 IF( LinearSystemTrialing ) NameSpaceI = MAX( 1, NameSpaceI ) 3210 3211 IF( NameSpaceI > 0 ) THEN 3212 CALL Info('DefaultSolve','Linear system namespace number: '//TRIM(I2S(NameSpaceI)),Level=7) 3213 CALL ListPushNamespace('linsys'//TRIM(I2S(NameSpaceI))//':') 3214 END IF 3215 3216 IF ( ListGetLogical( Params,'Linear System Save',Found )) THEN 3217 saveslot = GetString( Params,'Linear System Save Slot', Found ) 3218 IF(.NOT. Found .OR. TRIM( saveslot ) == 'solve') THEN 3219 CALL SaveLinearSystem( Solver ) 3220 END IF 3221 END IF 3222 3223 IF (PRESENT(BackRotNT)) THEN 3224 BackRot=GetLogical(Params,'Back Rotate N-T Solution',Found) 3225 IF(.NOT.Found) BackRot=.TRUE. 3226 3227 IF (BackRot.NEQV.BackRotNT) & 3228 CALL ListAddLogical(Params,'Back Rotate N-T Solution',BackRotNT) 3229 END IF 3230 3231 3232 IF( ListGetLogical( Params,'Harmonic Mode',Found ) ) THEN 3233 CALL ChangeToHarmonicSystem( Solver ) 3234 END IF 3235 3236 ! Combine the individual projectors into one massive projector 3237 CALL GenerateConstraintMatrix( CurrentModel, Solver ) 3238 3239 IF( GetLogical(Params,'Linear System Solver Disabled',Found) ) THEN 3240 CALL Info('DefaultSolve','Solver disabled, exiting early!',Level=10) 3241 RETURN 3242 END IF 3243 3244 SourceControl = ListGetLogical( Params,'Apply Source Control',Found ) 3245 IF(SourceControl) CALL ControlLinearSystem( Solver,PreSolve=.TRUE. ) 3246 3247 CALL Info('DefaultSolve','Calling SolveSystem for linear solution',Level=20) 3248 3249 A => Solver % Matrix 3250 x => Solver % Variable 3251 b => A % RHS 3252 SOL => x % Values 3253 3254 ! Debugging stuff activated only when "Max Output Level" >= 20 3255 IF( InfoActive( 20 ) ) THEN 3256 PRINT *,'range b'//TRIM(I2S(ParEnv % MyPe))//':', & 3257 MINVAL( b ), MAXVAL( b ), SUM( b ), SUM( ABS( b ) ) 3258 PRINT *,'range A'//TRIM(I2S(ParEnv % MyPe))//':', & 3259 MINVAL( A % Values ), MAXVAL( A % Values ), SUM( A % Values ), SUM( ABS(A % Values) ) 3260 END IF 3261 3262 326310 CONTINUE 3264 3265 CALL SolveSystem(A,ParMatrix,b,SOL,x % Norm,x % DOFs,Solver) 3266 3267 IF( InfoActive( 20 ) ) THEN 3268 PRINT *,'range x'//TRIM(I2S(ParEnv % MyPe))//':', & 3269 MINVAL( SOL ), MAXVAL( SOL ), SUM( SOL ), SUM( ABS( SOL ) ) 3270 END IF 3271 3272 3273 IF( LinearSystemTrialing ) THEN 3274 IF( x % LinConverged > 0 ) THEN 3275 IF( ListGetLogical( Params,'Linear System Trialing Conserve',Found ) ) THEN 3276 MaxCount = ListGetInteger( Params,'Linear System Trialing Conserve Rounds',Found ) 3277 IF( Found ) THEN 3278 i = NINT( ListGetConstReal( Params,'Linear System Namespace Number',Found ) ) 3279 IF( i == NameSpaceI ) THEN 3280 Count = 1 + ListGetInteger( Params,'Linear System Namespace Conserve Count',Found ) 3281 ELSE 3282 Count = 0 3283 END IF 3284 IF( Count > MaxCount ) THEN 3285 NameSpaceI = 0 3286 Count = 0 3287 END IF 3288 CALL ListAddInteger( Params,'Linear System Namespace Conserve Count',Count ) 3289 END IF 3290 CALL ListAddConstReal( Params,'Linear System Namespace Number', 1.0_dp *NameSpaceI ) 3291 END IF 3292 ELSE 3293 NameSpaceI = NameSpaceI + 1 3294 IF( .NOT. ListCheckPrefix( Params,'linsys'//TRIM(I2S(NameSpaceI)) ) ) THEN 3295 CALL Fatal('DefaultSolve','Exhausted all linear system strategies!') 3296 END IF 3297 CALL ListPopNamespace() 3298 CALL Info('DefaultSolve','Linear system namespace number: '//TRIM(I2S(NameSpaceI)),Level=7) 3299 CALL ListPushNamespace('linsys'//TRIM(I2S(NameSpaceI))//':') 3300 GOTO 10 3301 END IF 3302 END IF 3303 3304 IF(SourceControl) CALL ControlLinearSystem( Solver,PreSolve=.FALSE. ) 3305 3306 IF ( ListGetLogical( Params,'Linear System Save',Found )) THEN 3307 saveslot = GetString( Params,'Linear System Save Slot', Found ) 3308 IF( Found .AND. TRIM( saveslot ) == 'after') THEN 3309 CALL SaveLinearSystem( Solver ) 3310 END IF 3311 END IF 3312 3313 3314 ! If flux corrected transport is used then apply the corrector to the system 3315 IF( GetLogical( Params,'Linear System FCT',Found ) ) THEN 3316 CALL FCT_Correction( Solver ) 3317 END IF 3318 3319 ! Backchange the linear system 3320 IF( ListGetLogical( Params,'Harmonic Mode',Found ) ) THEN 3321 CALL ChangeToHarmonicSystem( Solver, .TRUE. ) 3322 END IF 3323 3324 IF (PRESENT(BackRotNT)) THEN 3325 IF (BackRot.NEQV.BackRotNT) & 3326 CALL ListAddLogical(Params,'Back Rotate N-T Solution',BackRot) 3327 END IF 3328 3329 Norm = x % Norm 3330 3331 IF( NameSpaceI > 0 ) CALL ListPopNamespace() 3332 3333 ! One can run postprocessing solver in this slot in every nonlinear iteration. 3334 !----------------------------------------------------------------------------- 3335 CALL DefaultSlaveSolvers(Solver,'Nonlinear Post Solvers') 3336 3337 3338 ! This could be somewhere else too. Now it is here for debugging. 3339 CALL SaveParallelInfo( Solver ) 3340 3341!------------------------------------------------------------------------------ 3342 END FUNCTION DefaultSolve 3343!------------------------------------------------------------------------------ 3344 3345 3346!------------------------------------------------------------------------------ 3347!> Is the system converged. Wrapper to hide the dirty test. 3348!------------------------------------------------------------------------------ 3349 FUNCTION DefaultConverged( USolver ) RESULT( Converged ) 3350!------------------------------------------------------------------------------ 3351 TYPE(Solver_t), OPTIONAL, TARGET, INTENT(in) :: USolver 3352 TYPE(Solver_t), POINTER :: Solver 3353 LOGICAL :: Converged 3354 3355 Solver => CurrentModel % Solver 3356 IF ( PRESENT( USolver ) ) Solver => USolver 3357 3358 Converged = ( Solver % Variable % NonlinConverged > 0 ) 3359 3360 END FUNCTION DefaultConverged 3361!------------------------------------------------------------------------------ 3362 3363 3364!------------------------------------------------------------------------------ 3365 FUNCTION DefaultLinesearch( Converged, USolver, FirstIter, nsize, values, values0 ) RESULT( ReduceStep ) 3366!------------------------------------------------------------------------------ 3367 LOGICAL, OPTIONAL :: Converged 3368 TYPE(Solver_t), OPTIONAL, TARGET :: USolver 3369 LOGICAL, OPTIONAL :: FirstIter 3370 INTEGER, OPTIONAL :: nsize 3371 REAL(KIND=dp), OPTIONAL, TARGET :: values(:), values0(:) 3372 LOGICAL :: ReduceStep 3373 3374 LOGICAL :: stat, First, Last, DoLinesearch 3375 TYPE(Solver_t), POINTER :: Solver 3376 TYPE(Variable_t), POINTER :: iterV 3377 INTEGER :: iter, previter, MaxIter 3378 REAL(KIND=dp) :: LinesearchCond 3379 3380 SAVE :: previter 3381 3382 IF ( PRESENT( USolver ) ) THEN 3383 Solver => USolver 3384 ELSE 3385 Solver => CurrentModel % Solver 3386 END IF 3387 3388 DoLinesearch = .FALSE. 3389 IF( ListCheckPrefix( Solver % Values,'Nonlinear System Linesearch') ) THEN 3390 LineSearchCond = ListGetCReal( Solver % Values,& 3391 'Nonlinear System Linesearch Condition', Stat ) 3392 IF( Stat ) THEN 3393 DoLinesearch = ( LineSearchCond > 0.0_dp ) 3394 CALL ListAddLogical( Solver % Values,'Nonlinear System Linesearch', DoLinesearch ) 3395 ELSE 3396 DoLinesearch = ListGetLogical( Solver % Values,'Nonlinear System Linesearch',Stat) 3397 END IF 3398 END IF 3399 3400 ! This routine might be called for convenience also without checking 3401 ! first whether it is needed. 3402 IF(.NOT. DoLinesearch ) THEN 3403 ReduceStep = .FALSE. 3404 IF( PRESENT( Converged ) ) Converged = .FALSE. 3405 RETURN 3406 END IF 3407 3408 IF( PRESENT( FirstIter ) ) THEN 3409 First = FirstIter 3410 Last = .FALSE. 3411 ELSE 3412 ! This is the first trial if we are the first nonlinear iteration 3413 ! for the first time. 3414 iterV => VariableGet( Solver % Mesh % Variables, 'nonlin iter' ) 3415 iter = NINT(iterV % Values(1)) 3416 MaxIter = ListGetInteger( Solver % Values,'Nonlinear System Max Iterations',Stat) 3417 First = (iter == 1 ) .AND. (iter /= previter) 3418 Last = (iter == MaxIter ) 3419 previter = iter 3420 END IF 3421 3422 ReduceStep = CheckStepSize(Solver,First,nsize,values,values0) 3423 3424 IF( Last .AND. .NOT. ReduceStep ) THEN 3425 CALL Info('DefaultLinesearch',& 3426 'Maximum number of nonlinear iterations reached, giving up after linesearch',Level=6) 3427 END IF 3428 3429 IF( PRESENT( Converged ) ) THEN 3430 Converged = ( Solver % Variable % NonlinConverged == 1 ) .OR. Last 3431 END IF 3432 3433 END FUNCTION DefaultLinesearch 3434 3435 3436 3437!------------------------------------------------------------------------------ 3438 SUBROUTINE DefaultUpdateEquationsR( G, F, UElement, USolver, BulkUpdate, VecAssembly ) 3439!------------------------------------------------------------------------------ 3440 TYPE(Solver_t), OPTIONAL, TARGET :: USolver 3441 TYPE(Element_t), OPTIONAL, TARGET :: UElement 3442 REAL(KIND=dp) :: G(:,:), f(:) 3443 LOGICAL, OPTIONAL :: BulkUpdate, VecAssembly 3444 3445 TYPE(Solver_t), POINTER :: Solver 3446 TYPE(Matrix_t), POINTER :: A 3447 TYPE(Variable_t), POINTER :: x 3448 TYPE(Element_t), POINTER :: Element, P1, P2 3449 REAL(KIND=dp), POINTER CONTIG :: b(:), svalues(:) 3450 3451 CHARACTER(LEN=MAX_NAME_LEN) :: str 3452 3453 LOGICAL :: Found, BUpd, VecAsm, MCAsm 3454 3455 INTEGER :: i, j, n, nd 3456 INTEGER(KIND=AddrInt) :: Proc 3457 INTEGER, POINTER CONTIG :: Indexes(:), PermIndexes(:) 3458 3459 IF ( PRESENT( USolver ) ) THEN 3460 Solver => USolver 3461 ELSE 3462 Solver => CurrentModel % Solver 3463 END IF 3464 A => Solver % Matrix 3465 x => Solver % Variable 3466 b => A % RHS 3467 3468 Element => GetCurrentElement( UElement ) 3469 3470 VecAsm = .FALSE. 3471 IF ( PRESENT( VecAssembly )) THEN 3472 VecAsm = VecAssembly 3473 END IF 3474 3475 IF ( ASSOCIATED(Element % BoundaryInfo) ) THEN 3476 Proc = Solver % BoundaryElementProcedure 3477 ELSE 3478 Proc = Solver % BulkElementProcedure 3479 END IF 3480 IF ( Proc /= 0 ) THEN 3481 n = GetElementNOFNodes( Element ) 3482 nd = GetElementNOFDOFs( Element, Solver ) 3483 CALL ExecLocalProc( Proc, CurrentModel, Solver, & 3484 G, F, Element, n, nd ) 3485 END IF 3486 3487 IF ( ParEnv % PEs > 1 ) THEN 3488 3489 IF ( ASSOCIATED(Element % BoundaryInfo) ) THEN 3490 P1 => Element % BoundaryInfo % Left 3491 P2 => Element % BoundaryInfo % Right 3492 IF ( ASSOCIATED(P1) .AND. ASSOCIATED(P2) ) THEN 3493 IF ( P1 % PartIndex /= ParEnv % myPE .AND. & 3494 P2 % PartIndex /= ParEnv % myPE )RETURN 3495 3496 IF ( P1 % PartIndex /= ParEnv % myPE .OR. & 3497 P2 % PartIndex /= ParEnv % myPE ) THEN 3498 G=G/2; F=F/2; 3499 END IF 3500 ELSE IF ( ASSOCIATED(P1) ) THEN 3501 IF ( P1 % PartIndex /= ParEnv % myPE ) RETURN 3502 ELSE IF ( ASSOCIATED(P2) ) THEN 3503 IF ( P2 % PartIndex /= ParEnv % myPE ) RETURN 3504 END IF 3505 ELSE IF ( Element % PartIndex/=ParEnv % myPE ) THEN 3506 IF(GetLogical(Solver % Values,'Linear System FCT',Found)) THEN 3507 Indexes => GetIndexStore() 3508 n = GetElementDOFs( Indexes, Element, Solver ) 3509 IF(.NOT.ASSOCIATED(A % HaloValues)) THEN 3510 ALLOCATE(A % HaloValues(SIZE(A % Values))); A % HaloValues=0._dp 3511 END IF 3512 CALL UpdateGlobalEquations( A,G,b,0._dp*f,n,x % DOFs, & 3513 x % Perm(Indexes(1:n)),UElement=Element,GlobalValues=A % HaloValues ) 3514 END IF 3515 RETURN 3516 END IF 3517 END IF 3518 3519 ! Vectorized version of the glueing process requested 3520 IF (VecAsm) THEN 3521#ifdef _OPENMP 3522 IF (OMP_GET_NUM_THREADS() == 1) THEN 3523 MCAsm = .TRUE. 3524 ELSE 3525 ! Check if multicoloured assembly is in use 3526 MCAsm = (Solver % CurrentColour > 0) .AND. & 3527 ASSOCIATED(Solver % ColourIndexList) 3528 END IF 3529#else 3530 MCAsm = .TRUE. 3531#endif 3532 Indexes => GetIndexStore() 3533 n = GetElementDOFs( Indexes, Element, Solver ) 3534 3535 PermIndexes => GetVecIndexStore() 3536 ! Get permuted indices 3537!DIR$ IVDEP 3538 DO j=1,n 3539 PermIndexes(j) = Solver % Variable % Perm(Indexes(j)) 3540 END DO 3541 3542 CALL UpdateGlobalEquationsVec( A, G, b, f, n, & 3543 x % DOFs, PermIndexes, & 3544 UElement=Element, MCAssembly=MCAsm ) 3545 ELSE 3546 Indexes => GetIndexStore() 3547 n = GetElementDOFs( Indexes, Element, Solver ) 3548 3549 ! If we have any antiperiodic entries we need to check them all! 3550 IF( Solver % PeriodicFlipActive ) THEN 3551 CALL FlipPeriodicLocalMatrix( Solver, n, Indexes, x % dofs, G ) 3552 CALL FlipPeriodicLocalForce( Solver, n, Indexes, x % dofs, f ) 3553 END IF 3554 3555 IF(Solver % DirectMethod == DIRECT_PERMON) THEN 3556 CALL UpdateGlobalEquations( A,G,b,f,n,x % DOFs, & 3557 x % Perm(Indexes(1:n)), UElement=Element ) 3558 CALL UpdatePermonMatrix( A, G, n, x % DOFs, x % Perm(Indexes(1:n)) ) 3559 ELSE 3560 CALL UpdateGlobalEquations( A,G,b,f,n,x % DOFs, & 3561 x % Perm(Indexes(1:n)), UElement=Element ) 3562 END IF 3563 3564 ! backflip, in case G is needed again 3565 IF( Solver % PeriodicFlipActive ) THEN 3566 CALL FlipPeriodicLocalMatrix( Solver, n, Indexes, x % dofs, G ) 3567 CALL FlipPeriodicLocalForce( Solver, n, Indexes, x % dofs, f ) 3568 END IF 3569 3570 END IF 3571!------------------------------------------------------------------------------ 3572 END SUBROUTINE DefaultUpdateEquationsR 3573!------------------------------------------------------------------------------ 3574 3575 3576!------------------------------------------------------------------------------ 3577 SUBROUTINE UpdatePermonMatrix(A,G,n,dofs,nind) 3578!------------------------------------------------------------------------------ 3579#ifdef HAVE_FETI4I 3580 use feti4i 3581#endif 3582 3583 TYPE(Matrix_t) :: A 3584 INTEGER :: n, dofs, nInd(:) 3585 REAL(KIND=dp) :: G(:,:) 3586!------------------------------------------------------------------------------ 3587 REAL(KIND=C_DOUBLE), ALLOCATABLE :: vals(:) 3588 INTEGER, POINTER :: ptr 3589 INTEGER :: i,j,k,l,k1,k2 3590 INTEGER :: matrixType, eType 3591 INTEGER(C_INT), ALLOCATABLE :: ind(:) 3592 3593 TYPE(Element_t), POINTER :: CElement 3594 3595#ifdef HAVE_FETI4I 3596!!$ INTERFACE 3597!!$ FUNCTION Permon_InitMatrix(n) RESULT(handle) BIND(C,Name="permon_init") 3598!!$ USE, INTRINSIC :: ISO_C_BINDING 3599!!$ TYPE(C_PTR) :: Handle 3600!!$ INTEGER(C_INT), VALUE :: n 3601!!$ END FUNCTION Permon_InitMatrix 3602!!$ 3603!!$ SUBROUTINE Permon_UpdateMatrix(handle,n,inds,vals) BIND(C,Name="permon_update") 3604!!$ USE, INTRINSIC :: ISO_C_BINDING 3605!!$ TYPE(C_PTR), VALUE :: Handle 3606!!$ INTEGER(C_INT), VALUE :: n 3607!!$ INTEGER(C_INT) :: inds(*) 3608!!$ REAL(C_DOUBLE) :: vals(*) 3609!!$ END SUBROUTINE Permon_UpdateMatrix 3610!!$ END INTERFACE 3611 3612 IF(.NOT. C_ASSOCIATED(A % PermonMatrix)) THEN 3613 A % NoDirichlet = .TRUE. 3614 !! A % PermonMatrix = Permon_InitMatrix(A % NumberOFRows) 3615 !! TODO: get correct matrix type 3616 matrixType = 0 !! symmetric positive definite (for other types see feti4i.h) 3617 CALL FETI4ICreateStiffnessMatrix(A % PermonMatrix, matrixType, 1) !TODO add number of rows A % NumberOFRows 3618 END IF 3619 3620 3621 ALLOCATE(vals(n*n*dofs*dofs), ind(n*dofs)) 3622 DO i=1,n 3623 DO j=1,dofs 3624 k1 = (i-1)*dofs + j 3625 DO k=1,n 3626 DO l=1,dofs 3627 k2 = (k-1)*dofs + l 3628 vals(dofs*n*(k1-1)+k2) = G(k1,k2) 3629 END DO 3630 END DO 3631 ind(k1) = dofs*(nInd(i)-1)+j 3632 END DO 3633 END DO 3634 3635 !CALL Permon_UpdateMatrix( A % PermonMatrix, n*dofs, ind, vals ) 3636 3637 CElement => GetCurrentElement() 3638 eType = ElementDim( CElement ) 3639 ! type of the element is the same as its dimension 3640 3641 CALL FETI4IAddElement(A % PermonMatrix, eType, n, nInd, n*dofs, ind, vals) 3642 3643#endif 3644 3645!------------------------------------------------------------------------------ 3646 END SUBROUTINE UpdatePermonMatrix 3647!------------------------------------------------------------------------------ 3648 3649 3650!------------------------------------------------------------------------------ 3651 SUBROUTINE DefaultUpdateEquationsC( GC, FC, UElement, USolver, BulkUpdate, MCAssembly ) 3652!------------------------------------------------------------------------------ 3653 TYPE(Solver_t), OPTIONAL, TARGET :: USolver 3654 TYPE(Element_t), OPTIONAL, TARGET :: UElement 3655 COMPLEX(KIND=dp) :: GC(:,:), FC(:) 3656 LOGICAL, OPTIONAL :: BulkUpdate, MCAssembly 3657 3658 TYPE(Solver_t), POINTER :: Solver 3659 TYPE(Matrix_t), POINTER :: A 3660 TYPE(Variable_t), POINTER :: x 3661 TYPE(Element_t), POINTER :: Element, P1, P2 3662 REAL(KIND=dp), POINTER CONTIG :: b(:) 3663 3664 REAL(KIND=dp), POINTER :: G(:,:), F(:) 3665 3666 LOGICAL :: Found, BUpd 3667 3668 INTEGER :: i,j,n,DOFs 3669 INTEGER, POINTER :: Indexes(:) 3670 3671 IF ( PRESENT( USolver ) ) THEN 3672 Solver => USolver 3673 ELSE 3674 Solver => CurrentModel % Solver 3675 END IF 3676 A => Solver % Matrix 3677 x => Solver % Variable 3678 b => A % RHS 3679 3680 Element => GetCurrentElement( UElement ) 3681 3682 DOFs = x % DOFs 3683 Indexes => GetIndexStore() 3684 n = GetElementDOFs( Indexes, Element, Solver ) 3685 3686 IF ( ParEnv % PEs > 1 ) THEN 3687 IF ( ASSOCIATED(Element % BoundaryInfo) ) THEN 3688 P1 => Element % BoundaryInfo % Left 3689 P2 => Element % BoundaryInfo % Right 3690 IF ( ASSOCIATED(P1) .AND. ASSOCIATED(P2) ) THEN 3691 IF ( P1 % PartIndex/=ParEnv % myPE .AND. & 3692 P2 % PartIndex/=ParEnv % myPE )RETURN 3693 3694 IF ( P1 % PartIndex/=ParEnv % myPE .OR. & 3695 P2 % PartIndex/=ParEnv % myPE ) THEN 3696 GC=GC/2; FC=FC/2; 3697 END IF 3698 ELSE IF ( ASSOCIATED(P1) ) THEN 3699 IF ( P1 % PartIndex/=ParEnv % myPE ) RETURN 3700 ELSE IF ( ASSOCIATED(P2) ) THEN 3701 IF ( P2 % PartIndex/=ParEnv % myPE ) RETURN 3702 END IF 3703 ELSE IF ( Element % PartIndex/=ParEnv % myPE ) THEN 3704 RETURN 3705 END IF 3706 END IF 3707 3708 ALLOCATE( G(DOFs*n,DOFs*n), F(DOFs*n) ) 3709 DO i=1,n*DOFs/2 3710 F( 2*(i-1)+1 ) = REAL( FC(i) ) 3711 F( 2*(i-1)+2 ) = AIMAG( FC(i) ) 3712 3713 DO j=1,n*DOFs/2 3714 G( 2*(i-1)+1, 2*(j-1)+1 ) = REAL( GC(i,j) ) 3715 G( 2*(i-1)+1, 2*(j-1)+2 ) = -AIMAG( GC(i,j) ) 3716 G( 2*(i-1)+2, 2*(j-1)+1 ) = AIMAG( GC(i,j) ) 3717 G( 2*(i-1)+2, 2*(j-1)+2 ) = REAL( GC(i,j) ) 3718 END DO 3719 END DO 3720 3721 ! If we have any antiperiodic entries we need to check them all! 3722 IF( Solver % PeriodicFlipActive ) THEN 3723 CALL FlipPeriodicLocalMatrix( Solver, n, Indexes, x % dofs, G ) 3724 CALL FlipPeriodicLocalForce( Solver, n, Indexes, x % dofs, f ) 3725 END IF 3726 3727 CALL UpdateGlobalEquations( A,G,b,f,n,x % DOFs,x % Perm(Indexes(1:n)) ) 3728 3729 DEALLOCATE( G, F) 3730!------------------------------------------------------------------------------ 3731 END SUBROUTINE DefaultUpdateEquationsC 3732!------------------------------------------------------------------------------ 3733 3734!> Adds the elementwise contribution the right-hand-side of the real valued matrix equation 3735!------------------------------------------------------------------------------ 3736 SUBROUTINE DefaultUpdateForceR( F, UElement, USolver, BulkUpdate ) 3737!------------------------------------------------------------------------------ 3738 REAL(KIND=dp) :: F(:) 3739 TYPE(Solver_t), OPTIONAL, TARGET :: USolver 3740 TYPE(Element_t), OPTIONAL, TARGET :: UElement 3741 LOGICAL, OPTIONAL :: BulkUpdate 3742 3743 TYPE(Solver_t), POINTER :: Solver 3744 TYPE(Variable_t), POINTER :: x 3745 TYPE(Element_t), POINTER :: Element, P1, P2 3746 3747 LOGICAL :: Found, BUpd 3748 3749 INTEGER :: n 3750 INTEGER, POINTER :: Indexes(:) 3751 3752 Solver => CurrentModel % Solver 3753 IF ( PRESENT(USolver) ) Solver => USolver 3754 3755 Element => GetCurrentElement( UElement ) 3756 3757 x => Solver % Variable 3758 Indexes => GetIndexStore() 3759 n = GetElementDOFs( Indexes, Element, Solver ) 3760 3761 IF ( ParEnv % PEs > 1 ) THEN 3762 IF ( ASSOCIATED(Element % BoundaryInfo) ) THEN 3763 P1 => Element % BoundaryInfo % Left 3764 P2 => Element % BoundaryInfo % Right 3765 IF ( ASSOCIATED(P1) .AND. ASSOCIATED(P2) ) THEN 3766 IF ( P1 % PartIndex/=ParEnv % myPE .AND. & 3767 P2 % PartIndex/=ParEnv % myPE )RETURN 3768 3769 IF ( P1 % PartIndex/=ParEnv % myPE .OR. & 3770 P2 % PartIndex/=ParEnv % myPE ) F=F/2; 3771 ELSE IF ( ASSOCIATED(P1) ) THEN 3772 IF ( P1 % PartIndex /= ParEnv % myPE ) RETURN 3773 ELSE IF ( ASSOCIATED(P2) ) THEN 3774 IF ( P2 % PartIndex /= ParEnv % myPE ) RETURN 3775 END IF 3776 ELSE IF ( Element % PartIndex/=ParEnv % myPE ) THEN 3777 RETURN 3778 END IF 3779 END IF 3780 3781 ! If we have any antiperiodic entries we need to check them all! 3782 IF( Solver % PeriodicFlipActive ) THEN 3783 CALL FlipPeriodicLocalForce( Solver, n, Indexes, x % dofs, f ) 3784 END IF 3785 3786 CALL UpdateGlobalForce( Solver % Matrix % RHS, & 3787 F, n, x % DOFs, x % Perm(Indexes(1:n)), UElement=Element) 3788 3789 IF( Solver % PeriodicFlipActive ) THEN 3790 CALL FlipPeriodicLocalForce( Solver, n, Indexes, x % dofs, f ) 3791 END IF 3792 3793 3794!------------------------------------------------------------------------------ 3795 END SUBROUTINE DefaultUpdateForceR 3796!------------------------------------------------------------------------------ 3797 3798 3799!> Adds the elementwise contribution the right-hand-side of the complex valued matrix equation 3800!------------------------------------------------------------------------------ 3801 SUBROUTINE DefaultUpdateForceC( FC, UElement, USolver ) 3802!------------------------------------------------------------------------------ 3803 COMPLEX(KIND=dp) :: FC(:) 3804 TYPE(Solver_t), OPTIONAL, TARGET :: USolver 3805 TYPE(Element_t), OPTIONAL, TARGET :: UElement 3806 3807 TYPE(Solver_t), POINTER :: Solver 3808 TYPE(Variable_t), POINTER :: x 3809 TYPE(Element_t), POINTER :: Element, P1, P2 3810 3811 REAL(KIND=dp), ALLOCATABLE :: F(:) 3812 3813 INTEGER :: i,n,DOFs 3814 INTEGER, POINTER :: Indexes(:) 3815 3816 Solver => CurrentModel % Solver 3817 IF ( PRESENT(USolver) ) Solver => USolver 3818 3819 Element => GetCurrentElement( UElement ) 3820 3821 x => Solver % Variable 3822 DOFs = x % DOFs 3823 Indexes => GetIndexStore() 3824 n = GetElementDOFs( Indexes, Element, Solver ) 3825 3826 IF ( ParEnv % PEs > 1 ) THEN 3827 IF ( ASSOCIATED(Element % BoundaryInfo) ) THEN 3828 P1 => Element % BoundaryInfo % Left 3829 P2 => Element % BoundaryInfo % Right 3830 IF ( ASSOCIATED(P1) .AND. ASSOCIATED(P2) ) THEN 3831 IF ( P1 % PartIndex/=ParEnv % myPE .AND. & 3832 P2 % PartIndex/=ParEnv % myPE )RETURN 3833 3834 IF ( P1 % PartIndex/=ParEnv % myPE .OR. & 3835 P2 % PartIndex/=ParEnv % myPE ) FC=FC/2; 3836 ELSE IF ( ASSOCIATED(P1) ) THEN 3837 IF ( P1 % PartIndex/=ParEnv % myPE ) RETURN 3838 ELSE IF ( ASSOCIATED(P2) ) THEN 3839 IF ( P2 % PartIndex/=ParEnv % myPE ) RETURN 3840 END IF 3841 ELSE IF ( Element % PartIndex/=ParEnv % myPE ) THEN 3842 RETURN 3843 END IF 3844 END IF 3845 3846 3847 ALLOCATE( F(DOFs*n) ) 3848 DO i=1,n*DOFs/2 3849 F( 2*(i-1) + 1 ) = REAL(FC(i)) 3850 F( 2*(i-1) + 2 ) = AIMAG(FC(i)) 3851 END DO 3852 3853 IF( Solver % PeriodicFlipActive ) THEN 3854 CALL FlipPeriodicLocalForce( Solver, n, Indexes, x % dofs, f ) 3855 END IF 3856 3857 CALL UpdateGlobalForce( Solver % Matrix % RHS, & 3858 F, n, x % DOFs, x % Perm(Indexes(1:n)) ) 3859 3860 DEALLOCATE( F ) 3861!------------------------------------------------------------------------------ 3862 END SUBROUTINE DefaultUpdateForceC 3863!------------------------------------------------------------------------------ 3864 3865 3866!------------------------------------------------------------------------------ 3867 SUBROUTINE DefaultUpdateTimeForceR( F, UElement, USolver ) 3868!------------------------------------------------------------------------------ 3869 REAL(KIND=dp) :: F(:) 3870 TYPE(Solver_t), OPTIONAL, TARGET :: USolver 3871 TYPE(Element_t), OPTIONAL, TARGET :: UElement 3872 3873 TYPE(Solver_t), POINTER :: Solver 3874 TYPE(Variable_t), POINTER :: x 3875 TYPE(Element_t), POINTER :: Element, P1, P2 3876 3877 INTEGER :: n 3878 INTEGER, POINTER :: Indexes(:) 3879 3880 Solver => CurrentModel % Solver 3881 IF ( PRESENT(USolver) ) Solver => USolver 3882 3883 Element => GetCurrentElement( UElement ) 3884 3885 x => Solver % Variable 3886 Indexes => GetIndexStore() 3887 n = GetElementDOFs( Indexes, Element, Solver ) 3888 3889 IF ( ParEnv % PEs > 1 ) THEN 3890 IF ( ASSOCIATED(Element % BoundaryInfo) ) THEN 3891 P1 => Element % BoundaryInfo % Left 3892 P2 => Element % BoundaryInfo % Right 3893 IF ( ASSOCIATED(P1) .AND. ASSOCIATED(P2) ) THEN 3894 IF ( P1 % PartIndex/=ParEnv % myPE .AND. & 3895 P2 % PartIndex/=ParEnv % myPE )RETURN 3896 3897 IF ( P1 % PartIndex/=ParEnv % myPE .OR. & 3898 P2 % PartIndex/=ParEnv % myPE ) F=F/2; 3899 ELSE IF ( ASSOCIATED(P1) ) THEN 3900 IF ( P1 % PartIndex /= ParEnv % myPE ) RETURN 3901 ELSE IF ( ASSOCIATED(P2) ) THEN 3902 IF ( P2 % PartIndex /= ParEnv % myPE ) RETURN 3903 END IF 3904 ELSE IF ( Element % PartIndex/=ParEnv % myPE ) THEN 3905 RETURN 3906 END IF 3907 END IF 3908 3909 IF( Solver % PeriodicFlipActive ) THEN 3910 CALL FlipPeriodicLocalForce( Solver, n, Indexes, x % dofs, f ) 3911 END IF 3912 3913 CALL UpdateTimeForce( Solver % Matrix,Solver % Matrix % RHS, & 3914 F, n, x % DOFs, x % Perm(Indexes(1:n)) ) 3915 3916 IF( Solver % PeriodicFlipActive ) THEN 3917 CALL FlipPeriodicLocalForce( Solver, n, Indexes, x % dofs, f ) 3918 END IF 3919 3920!------------------------------------------------------------------------------ 3921 END SUBROUTINE DefaultUpdateTimeForceR 3922!------------------------------------------------------------------------------ 3923 3924 3925 3926!------------------------------------------------------------------------------ 3927 SUBROUTINE DefaultUpdateTimeForceC( FC, UElement, USolver ) 3928!------------------------------------------------------------------------------ 3929 COMPLEX(KIND=dp) :: FC(:) 3930 TYPE(Solver_t), OPTIONAL, TARGET :: USolver 3931 TYPE(Element_t), OPTIONAL, TARGET :: UElement 3932 3933 TYPE(Solver_t), POINTER :: Solver 3934 TYPE(Variable_t), POINTER :: x 3935 TYPE(Element_t), POINTER :: Element, P1, P2 3936 3937 REAL(KIND=dp), ALLOCATABLE :: F(:) 3938 3939 INTEGER :: i,n,DOFs 3940 INTEGER, POINTER :: Indexes(:) 3941 3942 Solver => CurrentModel % Solver 3943 IF ( PRESENT(USolver) ) Solver => USolver 3944 3945 Element => GetCurrentElement( UElement ) 3946 3947 x => Solver % Variable 3948 DOFs = x % DOFs 3949 Indexes => GetIndexStore() 3950 n = GetElementDOFs( Indexes, Element, Solver ) 3951 3952 IF ( ParEnv % PEs > 1 ) THEN 3953 IF ( ASSOCIATED(Element % BoundaryInfo) ) THEN 3954 P1 => Element % BoundaryInfo % Left 3955 P2 => Element % BoundaryInfo % Right 3956 IF ( ASSOCIATED(P1) .AND. ASSOCIATED(P2) ) THEN 3957 IF ( P1 % PartIndex/=ParEnv % myPE .AND. & 3958 P2 % PartIndex/=ParEnv % myPE )RETURN 3959 3960 IF ( P1 % PartIndex/=ParEnv % myPE .OR. & 3961 P2 % PartIndex/=ParEnv % myPE ) FC=FC/2; 3962 ELSE IF ( ASSOCIATED(P1) ) THEN 3963 IF ( P1 % PartIndex/=ParEnv % myPE ) RETURN 3964 ELSE IF ( ASSOCIATED(P2) ) THEN 3965 IF ( P2 % PartIndex/=ParEnv % myPE ) RETURN 3966 END IF 3967 ELSE IF ( Element % PartIndex/=ParEnv % myPE ) THEN 3968 RETURN 3969 END IF 3970 END IF 3971 3972 ALLOCATE( F(DOFs*n) ) 3973 DO i=1,n*DOFs/2 3974 F( 2*(i-1) + 1 ) = REAL(FC(i)) 3975 F( 2*(i-1) + 2 ) = -AIMAG(FC(i)) 3976 END DO 3977 3978 IF( Solver % PeriodicFlipActive ) THEN 3979 CALL FlipPeriodicLocalForce( Solver, n, Indexes, x % dofs, f ) 3980 END IF 3981 3982 CALL UpdateTimeForce( Solver % Matrix,Solver % Matrix % RHS, & 3983 F, n, x % DOFs, x % Perm(Indexes(1:n)) ) 3984 3985 DEALLOCATE( F ) 3986!------------------------------------------------------------------------------ 3987 END SUBROUTINE DefaultUpdateTimeForceC 3988!------------------------------------------------------------------------------ 3989 3990 3991 3992!------------------------------------------------------------------------------ 3993 SUBROUTINE DefaultUpdatePrecR( M, UElement, USolver ) 3994!------------------------------------------------------------------------------ 3995 TYPE(Solver_t), OPTIONAL,TARGET :: USolver 3996 TYPE(Element_t), OPTIONAL, TARGET :: UElement 3997 REAL(KIND=dp) :: M(:,:) 3998 3999 TYPE(Solver_t), POINTER :: Solver 4000 TYPE(Matrix_t), POINTER :: A 4001 TYPE(Variable_t), POINTER :: x 4002 TYPE(Element_t), POINTER :: Element, P1, P2 4003 4004 INTEGER :: i,j,n 4005 INTEGER, POINTER :: Indexes(:) 4006 4007 IF ( PRESENT( USolver ) ) THEN 4008 Solver => USolver 4009 A => Solver % Matrix 4010 x => Solver % Variable 4011 ELSE 4012 Solver => CurrentModel % Solver 4013 A => Solver % Matrix 4014 x => Solver % Variable 4015 END IF 4016 4017 Element => GetCurrentElement( UElement ) 4018 4019 Indexes => GetIndexStore() 4020 n = GetElementDOFs( Indexes, Element, Solver ) 4021 4022 IF ( ParEnv % PEs > 1 ) THEN 4023 IF ( ASSOCIATED(Element % BoundaryInfo) ) THEN 4024 P1 => Element % BoundaryInfo % Left 4025 P2 => Element % BoundaryInfo % Right 4026 IF ( ASSOCIATED(P1) .AND. ASSOCIATED(P2) ) THEN 4027 IF ( P1 % PartIndex/=ParEnv % myPE .AND. & 4028 P2 % PartIndex/=ParEnv % myPE )RETURN 4029 4030 IF ( P1 % PartIndex/=ParEnv % myPE .OR. & 4031 P2 % PartIndex/=ParEnv % myPE ) M=M/2 4032 ELSE IF ( ASSOCIATED(P1) ) THEN 4033 IF ( P1 % PartIndex /= ParEnv % myPE ) RETURN 4034 ELSE IF ( ASSOCIATED(P2) ) THEN 4035 IF ( P2 % PartIndex /= ParEnv % myPE ) RETURN 4036 END IF 4037 ELSE IF ( Element % PartIndex/=ParEnv % myPE ) THEN 4038 RETURN 4039 END IF 4040 END IF 4041 4042!$OMP CRITICAL 4043 IF ( .NOT. ASSOCIATED( A % PrecValues ) ) THEN 4044 ALLOCATE( A % PrecValues(SIZE(A % Values)) ) 4045 A % PrecValues = 0.0d0 4046 END IF 4047!$OMP END CRITICAL 4048 4049 ! flip mass matrix for periodic elimination 4050 IF( Solver % PeriodicFlipActive ) THEN 4051 CALL FlipPeriodicLocalMatrix( Solver, n, Indexes, x % dofs, M ) 4052 END IF 4053 4054 CALL UpdateMassMatrix( A, M, n, x % DOFs, x % Perm(Indexes(1:n)), & 4055 A % PrecValues ) 4056 4057 IF( Solver % PeriodicFlipActive ) THEN 4058 CALL FlipPeriodicLocalMatrix( Solver, n, Indexes, x % dofs, M ) 4059 END IF 4060 4061!------------------------------------------------------------------------------ 4062 END SUBROUTINE DefaultUpdatePrecR 4063!------------------------------------------------------------------------------ 4064 4065!------------------------------------------------------------------------------ 4066 SUBROUTINE DefaultUpdatePrecC( MC, UElement, USolver ) 4067!------------------------------------------------------------------------------ 4068 TYPE(Solver_t), OPTIONAL,TARGET :: USolver 4069 TYPE(Element_t), OPTIONAL, TARGET :: UElement 4070 COMPLEX(KIND=dp) :: MC(:,:) 4071 4072 TYPE(Solver_t), POINTER :: Solver 4073 TYPE(Matrix_t), POINTER :: A 4074 TYPE(Variable_t), POINTER :: x 4075 TYPE(Element_t), POINTER :: Element, P1, P2 4076 4077 REAL(KIND=dp), ALLOCATABLE :: M(:,:) 4078 4079 INTEGER :: i,j,n,DOFs 4080 INTEGER, POINTER :: Indexes(:) 4081 4082 IF ( PRESENT( USolver ) ) THEN 4083 Solver => USolver 4084 A => Solver % Matrix 4085 x => Solver % Variable 4086 ELSE 4087 Solver => CurrentModel % Solver 4088 A => Solver % Matrix 4089 x => Solver % Variable 4090 END IF 4091 4092 Element => GetCurrentElement( UElement ) 4093 4094 DOFs = x % DOFs 4095 Indexes => GetIndexStore() 4096 n = GetElementDOFs( Indexes, Element, Solver ) 4097 4098 IF ( ParEnv % PEs > 1 ) THEN 4099 IF ( ASSOCIATED(Element % BoundaryInfo) ) THEN 4100 P1 => Element % BoundaryInfo % Left 4101 P2 => Element % BoundaryInfo % Right 4102 IF ( ASSOCIATED(P1) .AND. ASSOCIATED(P2) ) THEN 4103 IF ( P1 % PartIndex/=ParEnv % myPE .AND. & 4104 P2 % PartIndex/=ParEnv % myPE )RETURN 4105 4106 IF ( P1 % PartIndex/=ParEnv % myPE .OR. & 4107 P2 % PartIndex/=ParEnv % myPE ) MC=MC/2 4108 ELSE IF ( ASSOCIATED(P1) ) THEN 4109 IF ( P1 % PartIndex/=ParEnv % myPE ) RETURN 4110 ELSE IF ( ASSOCIATED(P2) ) THEN 4111 IF ( P2 % PartIndex/=ParEnv % myPE ) RETURN 4112 END IF 4113 ELSE IF ( Element % PartIndex/=ParEnv % myPE ) THEN 4114 RETURN 4115 END IF 4116 END IF 4117 4118!$OMP CRITICAL 4119 IF ( .NOT. ASSOCIATED( A % PrecValues ) ) THEN 4120 ALLOCATE( A % PrecValues(SIZE(A % Values)) ) 4121 A % PrecValues = 0.0d0 4122 END IF 4123!$OMP END CRITICAL 4124 4125 ALLOCATE( M(DOFs*n,DOFs*n) ) 4126 DO i=1,n*DOFs/2 4127 DO j=1,n*DOFs/2 4128 M(2*(i-1)+1, 2*(j-1)+1) = REAL( MC(i,j) ) 4129 M(2*(i-1)+1, 2*(j-1)+2) = -AIMAG( MC(i,j) ) 4130 M(2*(i-1)+2, 2*(j-1)+1) = AIMAG( MC(i,j) ) 4131 M(2*(i-1)+2, 2*(j-1)+2) = REAL( MC(i,j) ) 4132 END DO 4133 END DO 4134 4135 ! flip preconditioning matrix for periodic elimination 4136 IF( Solver % PeriodicFlipActive ) THEN 4137 CALL FlipPeriodicLocalMatrix( Solver, n, Indexes, x % dofs, M ) 4138 END IF 4139 4140 CALL UpdateMassMatrix( A, M, n, x % DOFs, x % Perm(Indexes(1:n)), & 4141 A % PrecValues ) 4142 DEALLOCATE( M ) 4143!------------------------------------------------------------------------------ 4144 END SUBROUTINE DefaultUpdatePrecC 4145!------------------------------------------------------------------------------ 4146 4147!------------------------------------------------------------------------------ 4148 SUBROUTINE DefaultUpdateMassR( M, UElement, USolver ) 4149!------------------------------------------------------------------------------ 4150 TYPE(Solver_t), OPTIONAL,TARGET :: USolver 4151 TYPE(Element_t), OPTIONAL, TARGET :: UElement 4152 REAL(KIND=dp) :: M(:,:) 4153 4154 TYPE(Solver_t), POINTER :: Solver 4155 TYPE(Matrix_t), POINTER :: A 4156 TYPE(Variable_t), POINTER :: x 4157 TYPE(Element_t), POINTER :: Element, P1, P2 4158 4159 LOGICAL :: Found 4160 4161 INTEGER :: i,j,n 4162 INTEGER, POINTER :: Indexes(:) 4163 4164 IF ( PRESENT( USolver ) ) THEN 4165 Solver => USolver 4166 A => Solver % Matrix 4167 x => Solver % Variable 4168 ELSE 4169 Solver => CurrentModel % Solver 4170 A => Solver % Matrix 4171 x => Solver % Variable 4172 END IF 4173 4174 Element => GetCurrentElement( UElement ) 4175 4176 Indexes => GetIndexStore() 4177 n = GetElementDOFs( Indexes, Element, Solver ) 4178 4179 IF ( ParEnv % PEs > 1 ) THEN 4180 IF ( ASSOCIATED(Element % BoundaryInfo) ) THEN 4181 P1 => Element % BoundaryInfo % Left 4182 P2 => Element % BoundaryInfo % Right 4183 IF ( ASSOCIATED(P1) .AND. ASSOCIATED(P2) ) THEN 4184 IF ( P1 % PartIndex/=ParEnv % myPE .AND. & 4185 P2 % PartIndex/=ParEnv % myPE )RETURN 4186 4187 IF ( P1 % PartIndex/=ParEnv % myPE .OR. & 4188 P2 % PartIndex/=ParEnv % myPE ) M=M/2 4189 ELSE IF ( ASSOCIATED(P1) ) THEN 4190 IF ( P1 % PartIndex /= ParEnv % myPE ) RETURN 4191 ELSE IF ( ASSOCIATED(P2) ) THEN 4192 IF ( P2 % PartIndex /= ParEnv % myPE ) RETURN 4193 END IF 4194 ELSE IF ( Element % PartIndex/=ParEnv % myPE ) THEN 4195 IF (ListGetLogical(Solver % Values, 'Linear System FCT', Found)) THEN 4196 Indexes => GetIndexStore() 4197 n = GetElementDOFs( Indexes, Element, Solver ) 4198 IF(.NOT.ASSOCIATED(A % HaloMassValues)) THEN 4199 ALLOCATE(A % HaloMassValues(SIZE(A % Values))); A % HaloMassValues=0._dp 4200 END IF 4201 CALL UpdateMassMatrix( A, M, n, x % DOFs, x % Perm(Indexes(1:n)), & 4202 A % HaloMassValues ) 4203 END IF 4204 RETURN 4205 END IF 4206 END IF 4207 4208!$OMP CRITICAL 4209 IF ( .NOT. ASSOCIATED( A % MassValues ) ) THEN 4210 ALLOCATE( A % MassValues(SIZE(A % Values)) ) 4211 A % MassValues = 0.0_dp 4212 END IF 4213!$OMP END CRITICAL 4214 4215 4216 ! flip mass matrix for periodic elimination 4217 IF( Solver % PeriodicFlipActive ) THEN 4218 CALL FlipPeriodicLocalMatrix( Solver, n, Indexes, x % dofs, M ) 4219 END IF 4220 4221 CALL UpdateMassMatrix( A, M, n, x % DOFs, x % Perm(Indexes(1:n)), & 4222 A % MassValues ) 4223 4224 ! backflip to be on the safe side 4225 IF( Solver % PeriodicFlipActive ) THEN 4226 CALL FlipPeriodicLocalMatrix( Solver, n, Indexes, x % Dofs, M ) 4227 END IF 4228 4229 4230!------------------------------------------------------------------------------ 4231 END SUBROUTINE DefaultUpdateMassR 4232!------------------------------------------------------------------------------ 4233 4234!------------------------------------------------------------------------------ 4235 SUBROUTINE DefaultUpdateMassC( MC, UElement, USolver ) 4236!------------------------------------------------------------------------------ 4237 TYPE(Solver_t), OPTIONAL,TARGET :: USolver 4238 TYPE(Element_t), OPTIONAL, TARGET :: UElement 4239 COMPLEX(KIND=dp) :: MC(:,:) 4240 4241 TYPE(Solver_t), POINTER :: Solver 4242 TYPE(Matrix_t), POINTER :: A 4243 TYPE(Variable_t), POINTER :: x 4244 TYPE(Element_t), POINTER :: Element, P1, P2 4245 4246 REAL(KIND=dp), ALLOCATABLE :: M(:,:) 4247 4248 INTEGER :: i,j,n,DOFs 4249 INTEGER, POINTER :: Indexes(:) 4250 4251 IF ( PRESENT( USolver ) ) THEN 4252 Solver => USolver 4253 A => Solver % Matrix 4254 x => Solver % Variable 4255 ELSE 4256 Solver => CurrentModel % Solver 4257 A => Solver % Matrix 4258 x => Solver % Variable 4259 END IF 4260 4261 Element => GetCurrentElement( UElement ) 4262 4263 DOFs = x % DOFs 4264 Indexes => GetIndexStore() 4265 n = GetElementDOFs( Indexes, Element, Solver ) 4266 4267 IF ( ParEnv % PEs > 1 ) THEN 4268 IF ( ASSOCIATED(Element % BoundaryInfo) ) THEN 4269 P1 => Element % BoundaryInfo % Left 4270 P2 => Element % BoundaryInfo % Right 4271 IF ( ASSOCIATED(P1) .AND. ASSOCIATED(P2) ) THEN 4272 IF ( P1 % PartIndex/=ParEnv % myPE .AND. & 4273 P2 % PartIndex/=ParEnv % myPE )RETURN 4274 4275 IF ( P1 % PartIndex/=ParEnv % myPE .OR. & 4276 P2 % PartIndex/=ParEnv % myPE ) MC=MC/2 4277 ELSE IF ( ASSOCIATED(P1) ) THEN 4278 IF ( P1 % PartIndex/=ParEnv % myPE ) RETURN 4279 ELSE IF ( ASSOCIATED(P2) ) THEN 4280 IF ( P2 % PartIndex/=ParEnv % myPE ) RETURN 4281 END IF 4282 ELSE IF ( Element % PartIndex/=ParEnv % myPE ) THEN 4283 RETURN 4284 END IF 4285 END IF 4286 4287!$OMP CRITICAL 4288 IF ( .NOT. ASSOCIATED( A % MassValues ) ) THEN 4289 ALLOCATE( A % MassValues(SIZE(A % Values)) ) 4290 A % MassValues = 0.0d0 4291 END IF 4292!$OMP END CRITICAL 4293 4294 ALLOCATE( M(DOFs*n,DOFs*n) ) 4295 DO i=1,n*DOFs/2 4296 DO j=1,n*DOFs/2 4297 M(2*(i-1)+1, 2*(j-1)+1) = REAL( MC(i,j) ) 4298 M(2*(i-1)+1, 2*(j-1)+2) = -AIMAG( MC(i,j) ) 4299 M(2*(i-1)+2, 2*(j-1)+1) = AIMAG( MC(i,j) ) 4300 M(2*(i-1)+2, 2*(j-1)+2) = REAL( MC(i,j) ) 4301 END DO 4302 END DO 4303 4304 ! flip mass matrix for periodic elimination 4305 IF( Solver % PeriodicFlipActive ) THEN 4306 CALL FlipPeriodicLocalMatrix( Solver, n, Indexes, x % dofs, M ) 4307 END IF 4308 4309 CALL UpdateMassMatrix( A, M, n, x % DOFs, x % Perm(Indexes(1:n)), & 4310 A % MassValues ) 4311 DEALLOCATE( M ) 4312!------------------------------------------------------------------------------ 4313 END SUBROUTINE DefaultUpdateMassC 4314!------------------------------------------------------------------------------ 4315 4316 4317 4318!------------------------------------------------------------------------------ 4319 RECURSIVE SUBROUTINE DefaultUpdateDampR( B, UElement, USolver ) 4320!------------------------------------------------------------------------------ 4321 TYPE(Solver_t), OPTIONAL, TARGET :: USolver 4322 TYPE(Element_t), OPTIONAL, TARGET :: UElement 4323 REAL(KIND=dp) :: B(:,:) 4324 4325 TYPE(Solver_t), POINTER :: Solver 4326 TYPE(Matrix_t), POINTER :: A 4327 TYPE(Variable_t), POINTER :: x 4328 TYPE(Element_t), POINTER :: Element, P1, P2 4329 4330 INTEGER :: i,j,n 4331 INTEGER, POINTER :: Indexes(:) 4332 4333 IF ( PRESENT( USolver ) ) THEN 4334 Solver => USolver 4335 ELSE 4336 Solver => CurrentModel % Solver 4337 END IF 4338 4339 A => Solver % Matrix 4340 x => Solver % Variable 4341 4342 Element => GetCurrentElement( UElement ) 4343 4344 Indexes => GetIndexStore() 4345 n = GetElementDOFs( Indexes, Element, Solver ) 4346 4347 IF ( ParEnv % PEs > 1 ) THEN 4348 IF ( ASSOCIATED(Element % BoundaryInfo) ) THEN 4349 P1 => Element % BoundaryInfo % Left 4350 P2 => Element % BoundaryInfo % Right 4351 IF ( ASSOCIATED(P1) .AND. ASSOCIATED(P2) ) THEN 4352 IF ( P1 % PartIndex/=ParEnv % myPE .AND. & 4353 P2 % PartIndex/=ParEnv % myPE )RETURN 4354 4355 IF ( P1 % PartIndex/=ParEnv % myPE .OR. & 4356 P2 % PartIndex/=ParEnv % myPE ) B=B/2; 4357 ELSE IF ( ASSOCIATED(P1) ) THEN 4358 IF ( P1 % PartIndex /= ParEnv % myPE ) RETURN 4359 ELSE IF ( ASSOCIATED(P2) ) THEN 4360 IF ( P2 % PartIndex /= ParEnv % myPE ) RETURN 4361 END IF 4362 ELSE IF ( Element % PartIndex/=ParEnv % myPE ) THEN 4363 RETURN 4364 END IF 4365 END IF 4366 4367!$OMP CRITICAL 4368 IF ( .NOT. ASSOCIATED( A % DampValues ) ) THEN 4369 ALLOCATE( A % DampValues(SIZE(A % Values)) ) 4370 A % DampValues = 0.0d0 4371 END IF 4372!$OMP END CRITICAL 4373 4374 ! flip damp matrix for periodic elimination 4375 IF( Solver % PeriodicFlipActive ) THEN 4376 CALL FlipPeriodicLocalMatrix( Solver, n, Indexes, x % dofs, B ) 4377 END IF 4378 4379 CALL UpdateMassMatrix( A, B, n, x % DOFs, x % Perm(Indexes(1:n)), & 4380 A % DampValues ) 4381 4382 IF( Solver % PeriodicFlipActive ) THEN 4383 CALL FlipPeriodicLocalMatrix( Solver, n, Indexes, x % dofs, B ) 4384 END IF 4385!------------------------------------------------------------------------------ 4386 END SUBROUTINE DefaultUpdateDampR 4387!------------------------------------------------------------------------------ 4388 4389 4390 4391!------------------------------------------------------------------------------ 4392 SUBROUTINE DefaultUpdateDampC( BC, UElement, USolver ) 4393!------------------------------------------------------------------------------ 4394 TYPE(Solver_t), OPTIONAL, TARGET :: USolver 4395 TYPE(Element_t), OPTIONAL, TARGET :: UElement 4396 COMPLEX(KIND=dp) :: BC(:,:) 4397 4398 TYPE(Solver_t), POINTER :: Solver 4399 TYPE(Matrix_t), POINTER :: A 4400 TYPE(Variable_t), POINTER :: x 4401 TYPE(Element_t), POINTER :: Element, P1, P2 4402 4403 REAL(KIND=dp), ALLOCATABLE :: B(:,:) 4404 4405 INTEGER :: i,j,n,DOFs 4406 INTEGER, POINTER :: Indexes(:) 4407 4408 IF ( PRESENT( USolver ) ) THEN 4409 Solver => USolver 4410 ELSE 4411 Solver => CurrentModel % Solver 4412 END IF 4413 4414 A => Solver % Matrix 4415 x => Solver % Variable 4416 4417 Element => GetCurrentElement( UElement ) 4418 4419 DOFs = x % DOFs 4420 Indexes => GetIndexStore() 4421 n = GetElementDOFs( Indexes, Element, Solver ) 4422 4423 IF ( ParEnv % PEs > 1 ) THEN 4424 IF ( ASSOCIATED(Element % BoundaryInfo) ) THEN 4425 P1 => Element % BoundaryInfo % Left 4426 P2 => Element % BoundaryInfo % Right 4427 IF ( ASSOCIATED(P1) .AND. ASSOCIATED(P2) ) THEN 4428 IF ( P1 % PartIndex/=ParEnv % myPE .AND. & 4429 P2 % PartIndex/=ParEnv % myPE )RETURN 4430 4431 IF ( P1 % PartIndex/=ParEnv % myPE .OR. & 4432 P2 % PartIndex/=ParEnv % myPE ) BC=BC/2 4433 ELSE IF ( ASSOCIATED(P1) ) THEN 4434 IF ( P1 % PartIndex/=ParEnv % myPE ) RETURN 4435 ELSE IF ( ASSOCIATED(P2) ) THEN 4436 IF ( P2 % PartIndex/=ParEnv % myPE ) RETURN 4437 END IF 4438 ELSE IF ( Element % PartIndex/=ParEnv % myPE ) THEN 4439 RETURN 4440 END IF 4441 END IF 4442 4443!$OMP CRITICAL 4444 IF ( .NOT. ASSOCIATED( A % DampValues ) ) THEN 4445 ALLOCATE( A % DampValues(SIZE(A % Values)) ) 4446 A % DampValues = 0.0d0 4447 END IF 4448!$OMP END CRITICAL 4449 4450 ALLOCATE( B(DOFs*n, DOFs*n) ) 4451 DO i=1,n*DOFs/2 4452 DO j=1,n*DOFs/2 4453 B(2*(i-1)+1, 2*(j-1)+1) = REAL( BC(i,j) ) 4454 B(2*(i-1)+1, 2*(j-1)+2) = -AIMAG( BC(i,j) ) 4455 B(2*(i-1)+2, 2*(j-1)+1) = AIMAG( BC(i,j) ) 4456 B(2*(i-1)+2, 2*(j-1)+2) = REAL( BC(i,j) ) 4457 END DO 4458 END DO 4459 4460 IF( Solver % PeriodicFlipActive ) THEN 4461 CALL FlipPeriodicLocalMatrix( Solver, n, Indexes, x % dofs, B ) 4462 END IF 4463 4464 CALL UpdateMassMatrix( A, B, n, x % DOFs, x % Perm(Indexes(1:n)), & 4465 A % DampValues ) 4466 DEALLOCATE( B ) 4467!------------------------------------------------------------------------------ 4468 END SUBROUTINE DefaultUpdateDampC 4469!------------------------------------------------------------------------------ 4470 4471 4472!------------------------------------------------------------------------------ 4473 SUBROUTINE DefaultUpdateBulkR( B, F, UElement, USolver ) 4474!------------------------------------------------------------------------------ 4475 TYPE(Solver_t), OPTIONAL, TARGET :: USolver 4476 TYPE(Element_t), OPTIONAL, TARGET :: UElement 4477 REAL(KIND=dp) :: B(:,:), F(:) 4478 4479 TYPE(Solver_t), POINTER :: Solver 4480 TYPE(Matrix_t), POINTER :: A 4481 TYPE(Variable_t), POINTER :: x 4482 TYPE(Element_t), POINTER :: Element, P1, P2 4483 4484 INTEGER :: i,j,n 4485 INTEGER, POINTER :: Indexes(:) 4486 4487 IF ( PRESENT( USolver ) ) THEN 4488 Solver => USolver 4489 ELSE 4490 Solver => CurrentModel % Solver 4491 END IF 4492 4493 A => Solver % Matrix 4494 x => Solver % Variable 4495 4496 Element => GetCurrentElement( UElement ) 4497 4498 Indexes => GetIndexStore() 4499 n = GetElementDOFs( Indexes, Element, Solver ) 4500 4501 IF ( ParEnv % PEs > 1 ) THEN 4502 IF ( ASSOCIATED(Element % BoundaryInfo) ) THEN 4503 P1 => Element % BoundaryInfo % Left 4504 P2 => Element % BoundaryInfo % Right 4505 IF ( ASSOCIATED(P1) .AND. ASSOCIATED(P2) ) THEN 4506 IF ( P1 % PartIndex/=ParEnv % myPE .AND. & 4507 P2 % PartIndex/=ParEnv % myPE )RETURN 4508 4509 IF ( P1 % PartIndex/=ParEnv % myPE .OR. & 4510 P2 % PartIndex/=ParEnv % myPE ) THEN 4511 B=B/2; F=F/2; 4512 END IF 4513 ELSE IF ( ASSOCIATED(P1) ) THEN 4514 IF ( P1 % PartIndex /= ParEnv % myPE ) RETURN 4515 ELSE IF ( ASSOCIATED(P2) ) THEN 4516 IF ( P2 % PartIndex /= ParEnv % myPE ) RETURN 4517 END IF 4518 ELSE IF ( Element % PartIndex/=ParEnv % myPE ) THEN 4519 RETURN 4520 END IF 4521 END IF 4522 4523!$OMP CRITICAL 4524 IF ( .NOT. ASSOCIATED( A % BulkValues ) ) THEN 4525 ALLOCATE( A % BulkValues(SIZE(A % Values)) ) 4526 A % BulkValues = 0.0_dp 4527 END IF 4528!$OMP END CRITICAL 4529 4530!$OMP CRITICAL 4531 IF ( .NOT. ASSOCIATED( A % BulkRHS ) ) THEN 4532 ALLOCATE( A % BulkRHS(SIZE(A % RHS)) ) 4533 A % BulkRHS = 0.0_dp 4534 END IF 4535!$OMP END CRITICAL 4536 4537 4538 ! If we have any antiperiodic entries we need to check them all! 4539 IF( Solver % PeriodicFlipActive ) THEN 4540 CALL FlipPeriodicLocalMatrix( Solver, n, Indexes, x % dofs, B ) 4541 CALL FlipPeriodicLocalForce( Solver, n, Indexes, x % dofs, f ) 4542 END IF 4543 4544 CALL UpdateGlobalEquations( A,B,A % BulkRHS,f,n,x % DOFs,x % Perm(Indexes(1:n)), & 4545 GlobalValues=A % BulkValues ) 4546 4547 IF( Solver % PeriodicFlipActive ) THEN 4548 CALL FlipPeriodicLocalMatrix( Solver, n, Indexes, x % dofs, B ) 4549 CALL FlipPeriodicLocalForce( Solver, n, Indexes, x % dofs, f ) 4550 END IF 4551 4552!------------------------------------------------------------------------------ 4553 END SUBROUTINE DefaultUpdateBulkR 4554!------------------------------------------------------------------------------ 4555 4556 4557 4558!------------------------------------------------------------------------------ 4559 SUBROUTINE DefaultUpdateBulkC( BC, FC, UElement, USolver ) 4560!------------------------------------------------------------------------------ 4561 TYPE(Solver_t), OPTIONAL, TARGET :: USolver 4562 TYPE(Element_t), OPTIONAL, TARGET :: UElement 4563 COMPLEX(KIND=dp) :: BC(:,:), FC(:) 4564 4565 TYPE(Solver_t), POINTER :: Solver 4566 TYPE(Matrix_t), POINTER :: A 4567 TYPE(Variable_t), POINTER :: x 4568 TYPE(Element_t), POINTER :: Element, P1, P2 4569 4570 REAL(KIND=dp), ALLOCATABLE :: B(:,:),F(:) 4571 4572 INTEGER :: i,j,n,DOFs 4573 INTEGER, POINTER :: Indexes(:) 4574 4575 IF ( PRESENT( USolver ) ) THEN 4576 Solver => USolver 4577 ELSE 4578 Solver => CurrentModel % Solver 4579 END IF 4580 4581 A => Solver % Matrix 4582 x => Solver % Variable 4583 4584 Element => GetCurrentElement( UElement ) 4585 4586 DOFs = x % DOFs 4587 Indexes => GetIndexStore() 4588 n = GetElementDOFs( Indexes, Element, Solver ) 4589 4590 IF ( ParEnv % PEs > 1 ) THEN 4591 IF ( ASSOCIATED(Element % BoundaryInfo) ) THEN 4592 P1 => Element % BoundaryInfo % Left 4593 P2 => Element % BoundaryInfo % Right 4594 IF ( ASSOCIATED(P1) .AND. ASSOCIATED(P2) ) THEN 4595 IF ( P1 % PartIndex/=ParEnv % myPE .AND. & 4596 P2 % PartIndex/=ParEnv % myPE )RETURN 4597 4598 IF ( P1 % PartIndex/=ParEnv % myPE .OR. & 4599 P2 % PartIndex/=ParEnv % myPE ) THEN 4600 BC=BC/2; FC=FC/2; 4601 END IF 4602 ELSE IF ( ASSOCIATED(P1) ) THEN 4603 IF ( P1 % PartIndex/=ParEnv % myPE ) RETURN 4604 ELSE IF ( ASSOCIATED(P2) ) THEN 4605 IF ( P2 % PartIndex/=ParEnv % myPE ) RETURN 4606 END IF 4607 ELSE IF ( Element % PartIndex/=ParEnv % myPE ) THEN 4608 RETURN 4609 END IF 4610 END IF 4611 4612 4613!$OMP CRITICAL 4614 IF ( .NOT. ASSOCIATED( A % BulkValues ) ) THEN 4615 ALLOCATE( A % BulkValues(SIZE(A % Values)) ) 4616 A % BulkValues = 0.0_dp 4617 END IF 4618!$OMP END CRITICAL 4619 4620!$OMP CRITICAL 4621 IF ( .NOT. ASSOCIATED( A % BulkRHS ) ) THEN 4622 ALLOCATE( A % BulkRHS(SIZE(A % RHS)) ) 4623 A % BulkRHS = 0.0_dp 4624 END IF 4625!$OMP END CRITICAL 4626 4627 ALLOCATE( B(DOFs*n, DOFs*n), F(DOFs*n) ) 4628 DO i=1,n*DOFs/2 4629 DO j=1,n*DOFs/2 4630 F( 2*(i-1)+1 ) = REAL( FC(i) ) 4631 F( 2*(i-1)+2 ) = AIMAG( FC(i) ) 4632 4633 B(2*(i-1)+1, 2*(j-1)+1) = REAL( BC(i,j) ) 4634 B(2*(i-1)+1, 2*(j-1)+2) = -AIMAG( BC(i,j) ) 4635 B(2*(i-1)+2, 2*(j-1)+1) = AIMAG( BC(i,j) ) 4636 B(2*(i-1)+2, 2*(j-1)+2) = REAL( BC(i,j) ) 4637 END DO 4638 END DO 4639 4640 IF( Solver % PeriodicFlipActive ) THEN 4641 CALL FlipPeriodicLocalMatrix( Solver, n, Indexes, x % dofs, B ) 4642 CALL FlipPeriodicLocalForce( Solver, n, Indexes, x % dofs, f ) 4643 END IF 4644 4645 CALL UpdateGlobalEquations( A,B,A % BulkRHS,f,n,x % DOFs,x % Perm(Indexes(1:n)), & 4646 GlobalValues=A % BulkValues ) 4647 4648 DEALLOCATE( B, F ) 4649!------------------------------------------------------------------------------ 4650 END SUBROUTINE DefaultUpdateBulkC 4651!------------------------------------------------------------------------------ 4652 4653 4654 SUBROUTINE DefaultUpdateDirichlet( Dvals, UElement, USolver, UIndexes, Dset ) 4655!------------------------------------------------------------------------------ 4656 REAL(KIND=dp) :: Dvals(:) 4657 TYPE(Solver_t), OPTIONAL,TARGET :: USolver 4658 TYPE(Element_t), OPTIONAL, TARGET :: UElement 4659 INTEGER, OPTIONAL, TARGET :: UIndexes(:) 4660 LOGICAL, OPTIONAL :: Dset(:) 4661 4662 TYPE(Solver_t), POINTER :: Solver 4663 TYPE(Matrix_t), POINTER :: A 4664 TYPE(Variable_t), POINTER :: x 4665 TYPE(Element_t), POINTER :: Element 4666 4667 INTEGER :: i,j,n 4668 INTEGER, POINTER :: Indexes(:) 4669 4670 4671 IF ( PRESENT( USolver ) ) THEN 4672 Solver => USolver 4673 ELSE 4674 Solver => CurrentModel % Solver 4675 END IF 4676 4677 A => Solver % Matrix 4678 x => Solver % Variable 4679 4680 Element => GetCurrentElement( UElement ) 4681 4682 IF( PRESENT( UIndexes ) ) THEN 4683 Indexes => Uindexes 4684 n = SIZE( Uindexes ) 4685 ELSE 4686 Indexes => GetIndexStore() 4687 n = GetElementDOFs( Indexes, Element, Solver ) 4688 END IF 4689 4690 DO i=1,n 4691 IF( PRESENT( Dset ) ) THEN 4692 IF( .NOT. Dset(i) ) CYCLE 4693 END IF 4694 CALL UpdateDirichletDof( A, Indexes(i), DVals(i) ) 4695 END DO 4696 4697 END SUBROUTINE DefaultUpdateDirichlet 4698 4699 4700 4701 4702!> Sets the Dirichlet conditions related to the variables of the active solver. 4703!------------------------------------------------------------------------------------------ 4704 SUBROUTINE DefaultDirichletBCs( USolver,Ux,UOffset,OffDiagonalMatrix) 4705!------------------------------------------------------------------------------------------ 4706 USE ElementDescription, ONLY: FaceElementOrientation 4707 IMPLICIT NONE 4708 4709 INTEGER, OPTIONAL :: UOffset 4710 LOGICAL, OPTIONAL :: OffDiagonalMatrix 4711 TYPE(Variable_t), OPTIONAL, TARGET :: Ux 4712 TYPE(Solver_t), OPTIONAL, TARGET :: USolver 4713!-------------------------------------------------------------------------------------------- 4714 TYPE(Matrix_t), POINTER :: A 4715 TYPE(Variable_t), POINTER :: x 4716 TYPE(Solver_t), POINTER :: Solver 4717 TYPE(ValueListEntry_t), POINTER :: ptr 4718 TYPE(ValueList_t), POINTER :: BC, Params 4719 TYPE(Element_t), POINTER :: Element, Parent, Edge, Face, SaveElement 4720 4721 REAL(KIND=dp), ALLOCATABLE :: Work(:), STIFF(:,:) 4722 REAL(KIND=dp), POINTER :: b(:) 4723 REAL(KIND=dp), POINTER :: DiagScaling(:) 4724 REAL(KIND=dp) :: xx, s, dval 4725 REAL(KIND=dp) :: DefaultDOFs(3) 4726 4727 INTEGER, ALLOCATABLE :: lInd(:), gInd(:) 4728 INTEGER :: FDofMap(4,3) 4729 INTEGER :: i, j, k, kk, l, m, n, nd, nb, np, mb, nn, ni, nj, i0 4730 INTEGER :: EDOFs, FDOFs, DOF, local, numEdgeDofs, istat, n_start, Offset 4731 INTEGER :: ActiveFaceId 4732 4733 LOGICAL :: RevertSign(4) ! Size for tetrahedra, update for other element shapes 4734 LOGICAL :: Flag,Found, ConstantValue, ScaleSystem, DirichletComm 4735 LOGICAL :: BUpd, PiolaTransform, QuadraticApproximation, SecondKindBasis 4736 4737 CHARACTER(LEN=MAX_NAME_LEN) :: name 4738 4739 SAVE gInd, lInd, STIFF, Work 4740!-------------------------------------------------------------------------------------------- 4741 4742 IF ( PRESENT( USolver ) ) THEN 4743 Solver => USolver 4744 ELSE 4745 Solver => CurrentModel % Solver 4746 END IF 4747 4748 Params => GetSolverParams(Solver) 4749 4750 IF ( GetString(Params,'Linear System Solver',Found)=='feti') THEN 4751 IF ( GetLogical(Params,'Total FETI', Found)) RETURN 4752 END IF 4753 4754 A => Solver % Matrix 4755 b => A % RHS 4756 IF ( PRESENT(Ux) ) THEN 4757 x => Ux 4758 ELSE 4759 x => Solver % Variable 4760 END IF 4761 4762 ! Create soft limiters to be later applied by the Dirichlet conditions 4763 ! This is done only once for each solver, hence the complex logic. 4764 !--------------------------------------------------------------------- 4765 IF( ListGetLogical( Solver % Values,'Apply Limiter',Found) ) THEN 4766 CALL DetermineSoftLimiter( Solver ) 4767 END IF 4768 4769 IF(.NOT.ALLOCATED(A % ConstrainedDOF)) THEN 4770 ALLOCATE(A % ConstrainedDOF(A % NumberOfRows)) 4771 A % ConstrainedDOF = .FALSE. 4772 END IF 4773 4774 IF(.NOT.ALLOCATED(A % Dvalues)) THEN 4775 ALLOCATE(A % Dvalues(A % NumberOfRows)) 4776 A % Dvalues = 0._dp 4777 END IF 4778 4779 Offset = 0 4780 IF(PRESENT(UOffset)) Offset=UOffset 4781 4782 n = Solver % Mesh % MaxElementDOFs 4783 IF ( .NOT. ALLOCATED( gInd ) ) THEN 4784 ALLOCATE( gInd(n), lInd(n), STIFF(n,n), Work(n), stat=istat ) 4785 IF ( istat /= 0 ) & 4786 CALL Fatal('DefUtils::DefaultDirichletBCs','Memory allocation failed.' ) 4787 ELSE IF ( SIZE(gInd) < n ) THEN 4788 DEALLOCATE( gInd, lInd, STIFF, Work ) 4789 ALLOCATE( gInd(n), lInd(n), STIFF(n,n), Work(n), stat=istat ) 4790 IF ( istat /= 0 ) & 4791 CALL Fatal('DefUtils::DefaultDirichletBCs','Memory allocation failed.' ) 4792 END IF 4793 4794 4795 IF ( x % DOFs > 1 ) THEN 4796 CALL SetDirichletBoundaries( CurrentModel,A, b, GetVarName(x),-1,x % DOFs,x % Perm ) 4797 END IF 4798 4799 CALL Info('DefUtils::DefaultDirichletBCs', & 4800 'Setting Dirichlet boundary conditions', Level=6) 4801 4802 ! ---------------------------------------------------------------------- 4803 ! Perform some preparations if BCs for p-approximation will be handled: 4804 ! ---------------------------------------------------------------------- 4805 ConstantValue = .FALSE. 4806 DO DOF=1,x % DOFs 4807 name = x % name 4808 IF ( x % DOFs > 1 ) name = ComponentName(name,DOF) 4809 4810 IF( .NOT. ListCheckPresentAnyBC( CurrentModel, name ) ) CYCLE 4811 4812 CALL Info('DefUtils::DefaultDirichletBCs', & 4813 'p-element preparations: '//TRIM(name), Level=15) 4814 4815 ! Clearing for p-approximation dofs associated with faces & edges: 4816 SaveElement => GetCurrentElement() 4817 DO i=1,Solver % Mesh % NumberOfBoundaryElements 4818 Element => GetBoundaryElement(i) 4819 IF ( .NOT. ActiveBoundaryElement() ) CYCLE 4820 4821 ! Get parent element: 4822 ! ------------------- 4823 Parent => Element % BoundaryInfo % Left 4824 IF ( .NOT. ASSOCIATED( Parent ) ) & 4825 Parent => Element % BoundaryInfo % Right 4826 4827 IF ( .NOT. ASSOCIATED(Parent) ) CYCLE 4828 4829 BC => GetBC() 4830 IF ( .NOT.ASSOCIATED(BC) ) CYCLE 4831 4832 ptr => ListFind(BC, Name,Found ) 4833 IF ( .NOT. ASSOCIATED(ptr) ) CYCLE 4834 4835 Constantvalue = ( ptr % type /= LIST_TYPE_CONSTANT_SCALAR_PROC ) 4836 4837 4838 IF ( isActivePElement(Parent)) THEN 4839 n = GetElementNOFNodes() 4840 ! Get indexes of boundary dofs: 4841 CALL getBoundaryIndexes( Solver % Mesh, Element, Parent, gInd, numEdgeDofs ) 4842 4843 DO k=n+1,numEdgeDofs 4844 nb = x % Perm( gInd(k) ) 4845 IF ( nb <= 0 ) CYCLE 4846 nb = Offset + x % DOFs * (nb-1) + DOF 4847 IF ( ConstantValue ) THEN 4848 A % ConstrainedDOF(nb) = .TRUE. 4849 A % Dvalues(nb) = 0._dp 4850 ELSE 4851 CALL ZeroRow( A, nb ) 4852 A % RHS(nb) = 0._dp 4853 END IF 4854 END DO 4855 ELSE 4856 ! To do: Check whether BCs for edge/face elements must be set via L2 projection. 4857 CYCLE 4858 END IF 4859 END DO 4860 4861 SaveElement => SetCurrentElement(SaveElement) 4862 END DO 4863 4864 4865 ! ------------------------------------------------------------------------------------- 4866 ! Set BCs for fields which are approximated using H1-conforming basis functions 4867 ! (either Lagrange basis or hierarchic p-basis): 4868 ! ------------------------------------------------------------------------------------- 4869 DO DOF=1,x % DOFs 4870 name = x % name 4871 IF (x % DOFs>1) name=ComponentName(name,DOF) 4872 4873 CALL SetNodalLoads( CurrentModel,A, b, & 4874 Name,DOF,x % DOFs,x % Perm ) ! , Offset ) not yet ? 4875 4876 CALL SetDirichletBoundaries( CurrentModel, A, b, & 4877 Name, DOF, x % DOFs, x % Perm, Offset, OffDiagonalMatrix ) 4878 4879 ! ---------------------------------------------------------------------------- 4880 ! Set Dirichlet BCs for edge and face dofs which come from approximating with 4881 ! p-elements: 4882 ! ---------------------------------------------------------------------------- 4883 IF( .NOT. ListCheckPresentAnyBC( CurrentModel, name ) ) CYCLE 4884 4885 CALL Info('DefUtils::DefaultDirichletBCs', & 4886 'p-element condition setup: '//TRIM(name), Level=15) 4887 4888 SaveElement => GetCurrentElement() 4889 DO i=1,Solver % Mesh % NumberOfBoundaryElements 4890 Element => GetBoundaryElement(i) 4891 IF ( .NOT. ActiveBoundaryElement() ) CYCLE 4892 4893 BC => GetBC() 4894 IF ( .NOT.ASSOCIATED(BC) ) CYCLE 4895 IF ( .NOT. ListCheckPresent(BC, Name) ) CYCLE 4896 4897 ! Get parent element: 4898 ! ------------------- 4899 Parent => Element % BoundaryInfo % Left 4900 IF ( .NOT. ASSOCIATED( Parent ) ) THEN 4901 Parent => Element % BoundaryInfo % Right 4902 END IF 4903 IF ( .NOT. ASSOCIATED( Parent ) ) CYCLE 4904 4905 ! Here set constraints for p-approximation only: 4906 ! ----------------------------------------------------- 4907 IF (.NOT.isActivePElement(Parent)) CYCLE 4908 4909 ptr => ListFind(BC, Name,Found ) 4910 Constantvalue = Ptr % Type /= LIST_TYPE_CONSTANT_SCALAR_PROC 4911 4912 IF ( ConstantValue ) CYCLE 4913 4914 4915 SELECT CASE(Parent % TYPE % DIMENSION) 4916 4917 CASE(2) 4918 ! If no edges do not try to set boundary conditions 4919 ! @todo This should changed to EXIT 4920 IF ( .NOT. ASSOCIATED( Solver % Mesh % Edges ) ) CYCLE 4921 4922 ! If boundary edge has no dofs move on to next edge 4923 IF (Element % BDOFs <= 0) CYCLE 4924 4925 ! Number of nodes for this element 4926 n = Element % TYPE % NumberOfNodes 4927 4928 ! Get indexes for boundary and values for dofs associated to them 4929 CALL getBoundaryIndexes( Solver % Mesh, Element, Parent, gInd, numEdgeDofs ) 4930 CALL LocalBcBDOFs( BC, Element, numEdgeDofs, Name, STIFF, Work ) 4931 4932 IF ( Solver % Matrix % Symmetric ) THEN 4933 4934 DO l=1,n 4935 nb = x % Perm( gInd(l) ) 4936 IF ( nb <= 0 ) CYCLE 4937 nb = Offset + x % DOFs * (nb-1) + DOF 4938 4939 s = A % Dvalues(nb) 4940 DO k=n+1,numEdgeDOFs 4941 Work(k) = Work(k) - s*STIFF(k,l) 4942 END DO 4943 END DO 4944 4945 DO k=n+1,numEdgeDOFs 4946 DO l=n+1,numEdgeDOFs 4947 STIFF(k-n,l-n) = STIFF(k,l) 4948 END DO 4949 Work(k-n) = Work(k) 4950 END DO 4951 l = numEdgeDOFs-n 4952 IF ( l==1 ) THEN 4953 Work(1) = Work(1)/STIFF(1,1) 4954 ELSE 4955 CALL SolveLinSys(STIFF(1:l,1:l),Work(1:l),l) 4956 END IF 4957 4958 DO k=n+1,numEdgeDOFs 4959 nb = x % Perm( gInd(k) ) 4960 IF ( nb <= 0 ) CYCLE 4961 nb = Offset + x % DOFs * (nb-1) + DOF 4962 4963 A % ConstrainedDOF(nb) = .TRUE. 4964 A % Dvalues(nb) = Work(k-n) 4965 END DO 4966 ELSE 4967 4968 ! Contribute this boundary to global system 4969 ! (i.e solve global boundary problem) 4970 DO k=n+1,numEdgeDofs 4971 nb = x % Perm( gInd(k) ) 4972 IF ( nb <= 0 ) CYCLE 4973 nb = Offset + x % DOFs * (nb-1) + DOF 4974 A % RHS(nb) = A % RHS(nb) + Work(k) 4975 DO l=1,numEdgeDofs 4976 mb = x % Perm( gInd(l) ) 4977 IF ( mb <= 0 ) CYCLE 4978 mb = Offset + x % DOFs * (mb-1) + DOF 4979 DO kk=A % Rows(nb)+DOF-1,A % Rows(nb+1)-1,x % DOFs 4980 IF ( A % Cols(kk) == mb ) THEN 4981 A % Values(kk) = A % Values(kk) + STIFF(k,l) 4982 EXIT 4983 END IF 4984 END DO 4985 END DO 4986 END DO 4987 END IF 4988 4989 CASE(3) 4990 ! If no faces present do not try to set boundary conditions 4991 ! @todo This should be changed to EXIT 4992 IF ( .NOT. ASSOCIATED(Solver % Mesh % Faces) ) CYCLE 4993 4994 ! Parameters of element 4995 n = Element % TYPE % NumberOfNodes 4996 4997 ! Get global boundary indexes and solve dofs associated to them 4998 CALL getBoundaryIndexes( Solver % Mesh, Element, & 4999 Parent, gInd, numEdgeDofs ) 5000 5001 ! If boundary face has no dofs skip to next boundary element 5002 IF (numEdgeDOFs == n) CYCLE 5003 5004 ! Get local solution 5005 CALL LocalBcBDofs( BC, Element, numEdgeDofs, Name, STIFF, Work ) 5006 5007 n_start = 1 5008 IF ( Solver % Matrix % Symmetric ) THEN 5009 DO l=1,n 5010 nb = x % Perm( gInd(l) ) 5011 IF ( nb <= 0 ) CYCLE 5012 nb = Offset + x % DOFs * (nb-1) + DOF 5013 5014 s = A % Dvalues(nb) 5015 DO k=n+1,numEdgeDOFs 5016 Work(k) = Work(k) - s*STIFF(k,l) 5017 END DO 5018 END DO 5019 n_start=n+1 5020 END IF 5021 5022 ! Contribute this entry to global boundary problem 5023 DO k=n+1,numEdgeDOFs 5024 nb = x % Perm( gInd(k) ) 5025 IF ( nb <= 0 ) CYCLE 5026 nb = Offset + x % DOFs * (nb-1) + DOF 5027 A % RHS(nb) = A % RHS(nb) + Work(k) 5028 DO l=n_start,numEdgeDOFs 5029 mb = x % Perm( gInd(l) ) 5030 IF ( mb <= 0 ) CYCLE 5031 mb = Offset + x % DOFs * (mb-1) + DOF 5032 DO kk=A % Rows(nb)+DOF-1,A % Rows(nb+1)-1,x % DOFs 5033 IF ( A % Cols(kk) == mb ) THEN 5034 A % Values(kk) = A % Values(kk) + STIFF(k,l) 5035 EXIT 5036 END IF 5037 END DO 5038 END DO 5039 END DO 5040 END SELECT 5041 END DO 5042 5043 SaveElement => SetCurrentElement(SaveElement) 5044 END DO 5045 5046 ! 5047 ! Apply special couple loads for 3-D models of solids: 5048 ! 5049 CALL SetCoupleLoads(CurrentModel, x % Perm, A, b, x % DOFs ) 5050 5051 ! ---------------------------------------------------------------------------- 5052 ! Set Dirichlet BCs for edge and face dofs which arise from approximating with 5053 ! edge (curl-conforming) or face (div-conforming) elements: 5054 ! ---------------------------------------------------------------------------- 5055 QuadraticApproximation = ListGetLogical(Solver % Values, 'Quadratic Approximation', Found) 5056 SecondKindBasis = ListGetLogical(Solver % Values, 'Second Kind Basis', Found) 5057 DO DOF=1,x % DOFs 5058 name = x % name 5059 IF (x % DOFs>1) name=ComponentName(name,DOF) 5060 5061 IF ( .NOT. ListCheckPrefixAnyBC(CurrentModel, TRIM(Name)//' {e}') .AND. & 5062 .NOT. ListCheckPrefixAnyBC(CurrentModel, TRIM(Name)//' {f}') ) CYCLE 5063 5064 CALL Info('SetDefaultDirichlet','Setting edge and face dofs',Level=15) 5065 5066 SaveElement => GetCurrentElement() 5067 DO i=1,Solver % Mesh % NumberOfBoundaryElements 5068 Element => GetBoundaryElement(i) 5069 5070 BC => GetBC() 5071 IF ( .NOT.ASSOCIATED(BC) ) CYCLE 5072 IF ( .NOT. ListCheckPrefix(BC, TRIM(Name)//' {e}') .AND. & 5073 .NOT. ListCheckPrefix(BC, TRIM(Name)//' {f}') ) CYCLE 5074 5075 ! Get parent element: 5076 ! ------------------- 5077 Parent => Element % BoundaryInfo % Left 5078 IF ( .NOT. ASSOCIATED( Parent ) ) THEN 5079 Parent => Element % BoundaryInfo % Right 5080 END IF 5081 IF ( .NOT. ASSOCIATED( Parent ) ) CYCLE 5082 np = Parent % TYPE % NumberOfNodes 5083 5084 IF ( ListCheckPrefix(BC, TRIM(Name)//' {e}') ) THEN 5085 !-------------------------------------------------------------------------------- 5086 ! We now devote this branch for handling edge (curl-conforming) finite elements 5087 ! which, in addition to edge DOFs, may also have DOFs associated with faces. 5088 !-------------------------------------------------------------------------------- 5089 IF ( ASSOCIATED( Solver % Mesh % Edges ) ) THEN 5090 SELECT CASE(GetElementFamily()) 5091 CASE(2) 5092 DO j=1,Parent % TYPE % NumberOfEdges 5093 Edge => Solver % Mesh % Edges(Parent % EdgeIndexes(j)) 5094 n = 0 5095 DO k=1,Element % TYPE % NumberOfNodes 5096 DO l=1,Edge % TYPE % NumberOfNodes 5097 IF ( Element % NodeIndexes(k)==Edge % NodeIndexes(l)) n=n+1 5098 END DO 5099 END DO 5100 IF ( n==Element % TYPE % NumberOfNodes ) EXIT 5101 END DO 5102 5103 IF ( .NOT. ActiveBoundaryElement(Edge) ) CYCLE 5104 5105 EDOFs = Edge % BDOFs ! The number of DOFs associated with edges 5106 n = Edge % TYPE % NumberOfNodes 5107 CALL VectorElementEdgeDOFs(BC,Edge,n,Parent,np,TRIM(Name)//' {e}',Work, & 5108 EDOFs, SecondKindBasis) 5109 5110 n=GetElementDOFs(gInd,Edge) 5111 5112 n_start = Solver % Def_Dofs(2,Parent % BodyId,1)*Edge % NDOFs 5113 DO j=1,EDOFs 5114 k = n_start + j 5115 nb = x % Perm(gInd(k)) 5116 IF ( nb <= 0 ) CYCLE 5117 nb = Offset + x % DOFs*(nb-1) + DOF 5118 5119 A % ConstrainedDOF(nb) = .TRUE. 5120 A % Dvalues(nb) = Work(j) 5121 END DO 5122 5123 CASE(3,4) 5124 !Check that solver % mesh % faces exists? 5125!-- Better: IF ( ASSOCIATED(Parent % FaceIndexes) ) THEN 5126 DO j=1,Parent % TYPE % NumberOfFaces 5127 Face => Solver % Mesh % Faces(Parent % FaceIndexes(j)) 5128 IF ( GetElementFamily(Element)==GetElementFamily(Face) ) THEN 5129 n = 0 5130 DO k=1,Element % TYPE % NumberOfNodes 5131 DO l=1,Face % TYPE % NumberOfNodes 5132 IF ( Element % NodeIndexes(k)==Face % NodeIndexes(l)) n=n+1 5133 END DO 5134 END DO 5135 IF ( n==Face % TYPE % NumberOfNodes ) EXIT 5136 END IF 5137 END DO 5138 5139 IF ( .NOT. ActiveBoundaryElement(Face) ) CYCLE 5140 5141 ! --------------------------------------------------------------------- 5142 ! Set first constraints for DOFs associated with edges. Save the values 5143 ! of DOFs in the array Work(:), so that the possible remaining DOFs 5144 ! associated with the face can be computed after this. 5145 ! --------------------------------------------------------------------- 5146 i0 = 0 5147 DO l=1,Face % TYPE % NumberOfEdges 5148 Edge => Solver % Mesh % Edges(Face % EdgeIndexes(l)) 5149 EDOFs = Edge % BDOFs 5150 n = Edge % TYPE % NumberOfNodes 5151 5152 CALL VectorElementEdgeDOFs(BC, Edge, n, Parent, np, TRIM(Name)//' {e}', & 5153 Work(i0+1:i0+EDOFs), EDOFs, SecondKindBasis) 5154 5155 n = GetElementDOFs(gInd,Edge) 5156 5157 n_start = Solver % Def_Dofs(2,Parent % BodyId,1)*Edge % NDOFs 5158 5159 DO j=1,EDOFs 5160 5161 k = n_start + j 5162 nb = x % Perm(gInd(k)) 5163 IF ( nb <= 0 ) CYCLE 5164 nb = Offset + x % DOFs*(nb-1) + DOF 5165 5166 A % ConstrainedDOF(nb) = .TRUE. 5167 A % Dvalues(nb) = Work(i0+j) 5168 END DO 5169 i0 = i0 + EDOFs 5170 END DO 5171 5172 ! --------------------------------------------------------------------- 5173 ! Set constraints for face DOFs via seeking the best approximation in L2. 5174 ! We use the variational equation (u x n,v') = (g x n - u0 x n,v) where 5175 ! u0 denotes the part of the interpolating function u+u0 which is already 5176 ! known and v is a test function for the Galerkin method. 5177 ! --------------------------------------------------------------------- 5178 IF (Face % BDOFs > 0) THEN 5179 EDOFs = i0 ! The count of edge DOFs set so far 5180 n = Face % TYPE % NumberOfNodes 5181 5182 CALL SolveLocalFaceDOFs(BC, Face, n, TRIM(Name)//' {e}', Work, EDOFs, & 5183 Face % BDOFs, QuadraticApproximation) 5184 5185 n = GetElementDOFs(GInd,Face) 5186 DO j=1,Face % BDOFs 5187 nb = x % Perm(GInd(n-Face % BDOFs+j)) ! The last entries should be face-DOF indices 5188 IF ( nb <= 0 ) CYCLE 5189 nb = Offset + x % DOFs*(nb-1) + DOF 5190 5191 A % ConstrainedDOF(nb) = .TRUE. 5192 A % Dvalues(nb) = Work(EDOFs+j) 5193 END DO 5194 END IF 5195 5196 END SELECT 5197 END IF 5198 ELSE IF ( ListCheckPrefix(BC, TRIM(Name)//' {f}') ) THEN 5199 !-------------------------------------------------------------------------- 5200 ! This branch should be able to handle BCs for face (div-conforming) 5201 ! elements. Now this works only for RT(0), ABF(0) and BMD(1) in 2D and 5202 ! for the Nedelec tetrahedron of the first and second kind in 3D. 5203 !-------------------------------------------------------------------------- 5204 SELECT CASE(GetElementFamily()) 5205 CASE(2) 5206 IF ( ASSOCIATED(Parent % EdgeIndexes) ) THEN 5207 DO j=1,Parent % TYPE % NumberOfEdges 5208 Edge => Solver % Mesh % Edges(Parent % EdgeIndexes(j)) 5209 n = 0 5210 DO k=1,Element % TYPE % NumberOfNodes 5211 DO l=1,Edge % TYPE % NumberOfNodes 5212 IF ( Element % NodeIndexes(k)==Edge % NodeIndexes(l)) n=n+1 5213 END DO 5214 END DO 5215 IF ( n==Element % TYPE % NumberOfNodes ) EXIT 5216 END DO 5217 5218 IF (n /= Element % TYPE % NumberOfNodes) CYCLE 5219 IF ( .NOT. ActiveBoundaryElement(Edge) ) CYCLE 5220 5221 EDOFs = Edge % BDOFs ! The number of DOFs associated with edges 5222 n = Edge % TYPE % NumberOfNodes 5223 CALL VectorElementEdgeDOFs(BC,Edge,n,Parent,np,TRIM(Name)//' {f}',Work, & 5224 EDOFs, SecondKindBasis, FaceElement=.TRUE.) 5225 5226 n=GetElementDOFs(gInd,Edge) 5227 5228 n_start = Solver % Def_Dofs(2,Parent % BodyId,1)*Edge % NDOFs 5229 DO j=1,EDOFs 5230 k = n_start + j 5231 nb = x % Perm(gInd(k)) 5232 IF ( nb <= 0 ) CYCLE 5233 nb = Offset + x % DOFs*(nb-1) + DOF 5234 5235 A % ConstrainedDOF(nb) = .TRUE. 5236 A % Dvalues(nb) = Work(j) 5237 END DO 5238 5239 END IF 5240 5241 CASE(3) 5242 IF ( ASSOCIATED( Parent % FaceIndexes ) ) THEN 5243 DO ActiveFaceId=1,Parent % TYPE % NumberOfFaces 5244 Face => Solver % Mesh % Faces(Parent % FaceIndexes(ActiveFaceId)) 5245 IF ( GetElementFamily(Element)==GetElementFamily(Face) ) THEN 5246 n = 0 5247 DO k=1,Element % TYPE % NumberOfNodes 5248 DO l=1,Face % TYPE % NumberOfNodes 5249 IF ( Element % NodeIndexes(k)==Face % NodeIndexes(l)) n=n+1 5250 END DO 5251 END DO 5252 IF ( n==Face % TYPE % NumberOfNodes ) EXIT 5253 END IF 5254 END DO 5255 5256 IF (n /= Element % TYPE % NumberOfNodes) CYCLE 5257 IF ( .NOT. ActiveBoundaryElement(Face) ) CYCLE 5258 5259 FDOFs = Face % BDOFs 5260 5261 IF (FDOFs > 0) THEN 5262 CALL FaceElementOrientation(Parent, RevertSign, ActiveFaceId) 5263 IF (SecondKindBasis) & 5264 CALL FaceElementBasisOrdering(Parent, FDofMap, ActiveFaceId) 5265 n = Face % TYPE % NumberOfNodes 5266 5267 CALL FaceElementDOFs(BC, Face, n, Parent, ActiveFaceId, & 5268 TRIM(Name)//' {f}', Work, FDOFs, SecondKindBasis) 5269 5270 IF (SecondKindBasis) THEN 5271 ! 5272 ! Conform to the orientation and ordering used in the 5273 ! assembly of the global equations 5274 ! 5275 DefaultDOFs(1:FDOFs) = Work(1:FDOFs) 5276 IF (RevertSign(ActiveFaceId)) THEN 5277 S = -1.0d0 5278 ELSE 5279 S = 1.0d0 5280 END IF 5281 5282 DO j=1,FDOFs 5283 k = FDofMap(ActiveFaceId,j) 5284 Work(j) = S * DefaultDOFs(k) 5285 END DO 5286 ELSE 5287 IF (RevertSign(ActiveFaceId)) Work(1:FDOFs) = -1.0d0*Work(1:FDOFs) 5288 END IF 5289 5290 n = GetElementDOFs(GInd,Face) 5291 ! 5292 ! Make an offset by the count of nodal DOFs. This provides 5293 ! the right starting point if edge DOFs are not present. 5294 ! 5295 n_start = Solver % Def_Dofs(3,Parent % BodyId,1) * Face % NDOFs 5296 ! 5297 ! Check if we need to increase the offset by the count of 5298 ! edge DOFs: 5299 ! 5300 IF ( ASSOCIATED(Face % EdgeIndexes) ) THEN 5301 EDOFs = 0 5302 DO l=1,Face % TYPE % NumberOfEdges 5303 Edge => Solver % Mesh % Edges(Face % EdgeIndexes(l)) 5304 EDOFs = EDOFs + Edge % BDOFs 5305 END DO 5306 n_start = n_start + Solver % Def_Dofs(3,Parent % BodyId,2) * EDOFs 5307 END IF 5308 5309 DO j=1,FDOFs 5310 k = n_start + j 5311 nb = x % Perm(gInd(k)) 5312 IF ( nb <= 0 ) CYCLE 5313 nb = Offset + x % DOFs*(nb-1) + DOF 5314 5315 A % ConstrainedDOF(nb) = .TRUE. 5316 A % Dvalues(nb) = Work(j) 5317 END DO 5318 END IF 5319 END IF 5320 5321 CASE DEFAULT 5322 CALL Warn('DefaultDirichletBCs', 'Cannot set face element DOFs for this element shape') 5323 END SELECT 5324 END IF 5325 END DO 5326 SaveElement => SetCurrentElement(SaveElement) 5327 END DO 5328 5329 5330 ! Add the possible constraint modes structures 5331 !---------------------------------------------------------- 5332 IF ( GetLogical(Solver % Values,'Constraint Modes Analysis',Found) ) THEN 5333 CALL SetConstraintModesBoundaries( CurrentModel, A, b, x % Name, x % DOFs, x % Perm ) 5334 END IF 5335 5336 5337#ifdef HAVE_FETI4I 5338 IF(C_ASSOCIATED(A % PermonMatrix)) THEN 5339 CALL Info('DefUtils::DefaultDirichletBCs','Permon matrix, Dirichlet conditions registered but not set!', Level=5) 5340 RETURN 5341 END IF 5342#endif 5343 5344 ! This is set outside so that it can be called more flexibilly 5345 CALL EnforceDirichletConditions( Solver, A, b ) 5346 5347 5348 CALL Info('DefUtils::DefaultDirichletBCs','Dirichlet boundary conditions set', Level=10) 5349!------------------------------------------------------------------------------ 5350 END SUBROUTINE DefaultDirichletBCs 5351!------------------------------------------------------------------------------ 5352 5353 5354! Solves a small dense linear system using Lapack routines 5355!------------------------------------------------------------------------------ 5356 SUBROUTINE SolveLinSys( A, x, n ) 5357!------------------------------------------------------------------------------ 5358 INTEGER :: n 5359 REAL(KIND=dp) :: A(n,n), x(n), b(n) 5360 5361 INTERFACE 5362 SUBROUTINE SolveLapack( N,A,x ) 5363 INTEGER N 5364 DOUBLE PRECISION A(n*n),x(n) 5365 END SUBROUTINE 5366 END INTERFACE 5367 5368!------------------------------------------------------------------------------ 5369 SELECT CASE(n) 5370 CASE(1) 5371 x(1) = x(1) / A(1,1) 5372 CASE(2) 5373 b = x 5374 CALL SolveLinSys2x2(A,x,b) 5375 CASE(3) 5376 b = x 5377 CALL SolveLinSys3x3(A,x,b) 5378 CASE DEFAULT 5379 CALL SolveLapack(n,A,x) 5380 END SELECT 5381!------------------------------------------------------------------------------ 5382 END SUBROUTINE SolveLinSys 5383!------------------------------------------------------------------------------ 5384 5385 5386!------------------------------------------------------------------------------ 5387!> This subroutine computes the values of DOFs that are associated with 5388!> mesh edges in the case of vector-valued (edge or face) finite elements, so that 5389!> the vector-valued interpolant of the BC data can be constructed. 5390!> The values of the DOFs are defined as D = S*(g.e,v)_E where the unit vector e 5391!> can be either tangential or normal to the edge, g is vector-valued data, 5392!> v is a polynomial on the edge E, and S reverts sign if necessary. 5393!------------------------------------------------------------------------------ 5394 SUBROUTINE VectorElementEdgeDOFs(BC, Element, n, Parent, np, Name, Integral, EDOFs, & 5395 SecondFamily, FaceElement) 5396!------------------------------------------------------------------------------ 5397 USE ElementDescription, ONLY: GetEdgeMap 5398 IMPLICIT NONE 5399 5400 TYPE(ValueList_t), POINTER :: BC !< The list of boundary condition values 5401 TYPE(Element_t), POINTER :: Element !< The boundary element handled 5402 INTEGER :: n !< The number of boundary element nodes 5403 TYPE(Element_t) :: Parent !< The parent element of the boundary element 5404 INTEGER :: np !< The number of parent element nodes 5405 CHARACTER(LEN=*) :: Name !< The name of boundary condition 5406 REAL(KIND=dp) :: Integral(:) !< The values of DOFs 5407 INTEGER, OPTIONAL :: EDOFs !< The number of DOFs 5408 LOGICAL, OPTIONAL :: SecondFamily !< To select the element family 5409 LOGICAL, OPTIONAL :: FaceElement !< If .TRUE., e is normal to the edge 5410!------------------------------------------------------------------------------ 5411 TYPE(Nodes_t), SAVE :: Nodes, Pnodes 5412 TYPE(ElementType_t), POINTER :: SavedType 5413 TYPE(GaussIntegrationPoints_t) :: IP 5414 5415 LOGICAL :: Lstat, RevertSign, SecondKindBasis, DivConforming 5416 INTEGER, POINTER :: Edgemap(:,:) 5417 INTEGER :: i,j,k,p,DOFs 5418 INTEGER :: i1,i2,i3 5419 5420 REAL(KIND=dp) :: Basis(n),Load(n),Vload(3,1:n),VL(3),e(3),d(3) 5421 REAL(KIND=dp) :: E21(3),E32(3) 5422 REAL(KIND=dp) :: u,v,L,s,DetJ 5423!------------------------------------------------------------------------------ 5424 DOFs = 1 5425 IF (PRESENT(EDOFs)) THEN 5426 IF (EDOFs > 2) THEN 5427 CALL Fatal('VectorElementEdgeDOFs','Cannot handle more than 2 DOFs per edge') 5428 ELSE 5429 DOFs = EDOFs 5430 END IF 5431 END IF 5432 5433 IF (PRESENT(SecondFamily)) THEN 5434 SecondKindBasis = SecondFamily 5435 IF (SecondKindBasis .AND. (DOFs /= 2) ) & 5436 CALL Fatal('VectorElementEdgeDOFs','2 DOFs per edge expected') 5437 ELSE 5438 SecondKindBasis = .FALSE. 5439 END IF 5440 5441 IF (PRESENT(FaceElement)) THEN 5442 DivConforming = FaceElement 5443 ELSE 5444 DivConforming = .FALSE. 5445 END IF 5446 5447 ! Get the nodes of the boundary and parent elements: 5448 CALL GetElementNodes(Nodes, Element) 5449 CALL GetElementNodes(PNodes, Parent) 5450 5451 RevertSign = .FALSE. 5452 EdgeMap => GetEdgeMap(GetElementFamily(Parent)) 5453 DO i=1,SIZE(EdgeMap,1) 5454 j=EdgeMap(i,1) 5455 k=EdgeMap(i,2) 5456 IF ( Parent % NodeIndexes(j)==Element % NodeIndexes(1) .AND. & 5457 Parent % NodeIndexes(k)==Element % NodeIndexes(2) ) THEN 5458 EXIT 5459 ELSE IF (Parent % NodeIndexes(j)==Element % NodeIndexes(2) .AND. & 5460 Parent % NodeIndexes(k)==Element % NodeIndexes(1) ) THEN 5461 ! This is the right edge but has opposite orientation as compared 5462 ! with the listing of the parent element edges 5463 RevertSign = .TRUE. 5464 EXIT 5465 END IF 5466 END DO 5467 5468 Load(1:n) = GetReal( BC, Name, Lstat, Element ) 5469 5470 i = LEN_TRIM(Name) 5471 VLoad(1,1:n)=GetReal(BC,Name(1:i)//' 1',Lstat,element) 5472 VLoad(2,1:n)=GetReal(BC,Name(1:i)//' 2',Lstat,element) 5473 VLoad(3,1:n)=GetReal(BC,Name(1:i)//' 3',Lstat,element) 5474 5475 e(1) = PNodes % x(k) - PNodes % x(j) 5476 e(2) = PNodes % y(k) - PNodes % y(j) 5477 e(3) = PNodes % z(k) - PNodes % z(j) 5478 e = e/SQRT(SUM(e**2)) 5479 IF (DivConforming) THEN 5480 ! The boundary normal is needed instead of the tangent vector. 5481 ! First, find the element director d that makes the parent 5482 ! element an oriented surface. 5483 i1 = EdgeMap(1,1) 5484 i2 = EdgeMap(1,2) 5485 i3 = EdgeMap(2,2) 5486 E21(1) = PNodes % x(i2) - PNodes % x(i1) 5487 E21(2) = PNodes % y(i2) - PNodes % y(i1) 5488 E21(3) = PNodes % z(i2) - PNodes % z(i1) 5489 E32(1) = PNodes % x(i3) - PNodes % x(i2) 5490 E32(2) = PNodes % y(i3) - PNodes % y(i2) 5491 E32(3) = PNodes % z(i3) - PNodes % z(i2) 5492 d = CrossProduct(E21, E32) 5493 d = d/SQRT(SUM(d**2)) 5494 ! Set e to be the outward normal to the parent element: 5495 e = CrossProduct(e, d) 5496 END IF 5497 5498 ! Is this element type stuff needed and for what? 5499 SavedType => Element % TYPE 5500 IF ( GetElementFamily()==1 ) Element % TYPE=>GetElementType(202) 5501 5502 Integral(1:DOFs) = 0._dp 5503 IP = GaussPoints(Element) 5504 DO p=1,IP % n 5505 Lstat = ElementInfo( Element, Nodes, IP % u(p), & 5506 IP % v(p), IP % w(p), DetJ, Basis ) 5507 s = IP % s(p) * DetJ 5508 5509 L = SUM(Load(1:n)*Basis(1:n)) 5510 VL = MATMUL(Vload(:,1:n),Basis(1:n)) 5511 5512 IF (SecondKindBasis) THEN 5513 u = IP % u(p) 5514 v = 0.5d0*(1.0d0-sqrt(3.0d0)*u) 5515 Integral(1)=Integral(1)+s*(L+SUM(VL*e))*v 5516 v = 0.5d0*(1.0d0+sqrt(3.0d0)*u) 5517 Integral(2)=Integral(2)+s*(L+SUM(VL*e))*v 5518 ELSE 5519 Integral(1)=Integral(1)+s*(L+SUM(VL*e)) 5520 5521 IF (.NOT. DivConforming) THEN 5522 ! This branch is concerned with the second-order curl-conforming elements 5523 IF (DOFs>1) THEN 5524 v = Basis(2)-Basis(1) 5525 ! The parent element must define the default for the positive tangent associated 5526 ! with the edge. Thus, if the boundary element handled has an opposite orientation, 5527 ! the sign must be reverted to get the positive coordinate associated with the 5528 ! parent element edge. 5529 IF (RevertSign) v = -1.0d0*v 5530 Integral(2)=Integral(2)+s*(L+SUM(VL*e))*v 5531 END IF 5532 END IF 5533 END IF 5534 END DO 5535 Element % TYPE => SavedType 5536 5537 j = Parent % NodeIndexes(j) 5538 IF ( ParEnv % PEs>1 ) & 5539 j=CurrentModel % Mesh % ParallelInfo % GlobalDOFs(j) 5540 5541 k = Parent % NodeIndexes(k) 5542 IF ( ParEnv % PEs>1 ) & 5543 k=CurrentModel % Mesh % ParallelInfo % GlobalDOFs(k) 5544 5545 IF (k < j) THEN 5546 IF (SecondKindBasis) THEN 5547 Integral(1)=-Integral(1) 5548 Integral(2)=-Integral(2) 5549 ELSE 5550 Integral(1)=-Integral(1) 5551 END IF 5552 END IF 5553!------------------------------------------------------------------------------ 5554 END SUBROUTINE VectorElementEdgeDOFs 5555!------------------------------------------------------------------------------ 5556 5557 5558!------------------------------------------------------------------------------ 5559!> This subroutine computes the values of DOFs that are associated with 5560!> mesh faces in the case of curl-conforming (edge) finite elements, so that 5561!> the edge finite element interpolant of the BC data can be constructed. 5562!> The values of the DOFs are obtained as the best approximation in L2 when 5563!> the values of the DOFs associated with edges are given. 5564!------------------------------------------------------------------------------ 5565 SUBROUTINE SolveLocalFaceDOFs(BC, Element, n, Name, DOFValues, & 5566 EDOFs, FDOFs, QuadraticApproximation) 5567!------------------------------------------------------------------------------ 5568 IMPLICIT NONE 5569 5570 TYPE(ValueList_t), POINTER :: BC !< The list of boundary condition values 5571 TYPE(Element_t), POINTER :: Element !< The boundary element handled 5572 INTEGER :: n !< The number of boundary element nodes 5573 CHARACTER(LEN=*) :: Name !< The name of boundary condition 5574 REAL(KIND=dp) :: DOFValues(:) !< The values of DOFs 5575 INTEGER :: EDOFs !< The number of edge DOFs 5576 INTEGER :: FDOFs !< The number of face DOFs 5577 LOGICAL :: QuadraticApproximation !< Use second-order edge element basis 5578!------------------------------------------------------------------------------ 5579 TYPE(Nodes_t), SAVE :: Nodes 5580 TYPE(GaussIntegrationPoints_t) :: IP 5581 5582 LOGICAL :: Lstat 5583 5584 INTEGER :: i,j,p,DOFs,BasisDegree 5585 5586 REAL(KIND=dp) :: Basis(n),Vload(3,n),VL(3),Normal(3) 5587 REAL(KIND=dp) :: EdgeBasis(EDOFs+FDOFs,3) 5588 REAL(KIND=dp) :: Mass(FDOFs,FDOFs), Force(FDOFs) 5589 REAL(KIND=dp) :: v,s,DetJ 5590!------------------------------------------------------------------------------ 5591 IF (QuadraticApproximation) THEN 5592 BasisDegree = 2 5593 ELSE 5594 BasisDegree = 1 5595 END IF 5596 5597 Mass = 0.0d0 5598 Force = 0.0d0 5599 5600 CALL GetElementNodes(Nodes, Element) 5601 5602 i = LEN_TRIM(Name) 5603 VLoad(1,1:n)=GetReal(BC,Name(1:i)//' 1',Lstat,element) 5604 VLoad(2,1:n)=GetReal(BC,Name(1:i)//' 2',Lstat,element) 5605 VLoad(3,1:n)=GetReal(BC,Name(1:i)//' 3',Lstat,element) 5606 5607 IP = GaussPoints(Element) 5608 DO p=1,IP % n 5609 5610 Lstat = EdgeElementInfo( Element, Nodes, IP % u(p), IP % v(p), IP % w(p), & 5611 DetF=DetJ, Basis=Basis, EdgeBasis=EdgeBasis, BasisDegree=BasisDegree, & 5612 ApplyPiolaTransform=.TRUE., TangentialTrMapping=.TRUE.) 5613 5614 Normal = NormalVector(Element, Nodes, IP % u(p), IP % v(p), .FALSE.) 5615 5616 VL = MATMUL(Vload(:,1:n),Basis(1:n)) 5617 5618 s = IP % s(p) * DetJ 5619 5620 DO i=1,FDOFs 5621 DO j=1,FDOFs 5622 Mass(i,j) = Mass(i,j) + SUM(EdgeBasis(EDOFs+i,:) * EdgeBasis(EDOFs+j,:)) * s 5623 END DO 5624 Force(i) = Force(i) + SUM(CrossProduct(VL,Normal) * EdgeBasis(EDOFs+i,:)) * s 5625 DO j=1,EDOFs 5626 Force(i) = Force(i) - DOFValues(j) * SUM(EdgeBasis(j,:) * EdgeBasis(EDOFs+i,:)) * s 5627 END DO 5628 END DO 5629 END DO 5630 5631 CALL LUSolve(FDOFs, Mass(1:FDOFs,1:FDOFs), Force(1:FDOFs)) 5632 DOFValues(EDOFs+1:EDOFs+FDOFs) = Force(1:FDOFs) 5633!------------------------------------------------------------------------------ 5634 END SUBROUTINE SolveLocalFaceDOFs 5635!------------------------------------------------------------------------------ 5636 5637!------------------------------------------------------------------------------ 5638!> This subroutine computes the values of DOFs that are associated with 5639!> mesh faces in the case of face (vector-valued) finite elements, so that 5640!> the vector-valued interpolant of the BC data can be constructed. 5641!> The values of the DOFs are defined as D = S*(g.n,v)_F where the unit vector n 5642!> is normal to the face, g is vector-valued data, v is a polynomial on the face F, 5643!> and S reverts sign if necessary. This subroutine performs neither sign 5644!> reversions nor the permutations of DOFs, i.e. the DOFs are returned in 5645!> the default form. 5646! TO DO: This can now handle only tetrahedral faces and needs an update when 5647! new 3-D element types are added. 5648!------------------------------------------------------------------------------ 5649 SUBROUTINE FaceElementDOFs(BC, Element, n, Parent, FaceId, Name, Integral, & 5650 FDOFs, SecondFamily) 5651!------------------------------------------------------------------------------ 5652 IMPLICIT NONE 5653 5654 TYPE(ValueList_t), POINTER, INTENT(IN) :: BC !< The list of boundary condition values 5655 TYPE(Element_t), POINTER, INTENT(IN) :: Element !< The boundary element handled 5656 INTEGER, INTENT(IN) :: n !< The number of boundary element nodes 5657 TYPE(Element_t), POINTER, INTENT(IN) :: Parent !< The parent element of the boundary element 5658 INTEGER, INTENT(IN) :: FaceId !< The parent element face corresponding to Element 5659 CHARACTER(LEN=*), INTENT(IN) :: Name !< The name of boundary condition 5660 REAL(KIND=dp), INTENT(OUT) :: Integral(:) !< The values of DOFs 5661 INTEGER, OPTIONAL, INTENT(IN) :: FDOFs !< The number of DOFs 5662 LOGICAL, OPTIONAL, INTENT(IN) :: SecondFamily !< To select the element family 5663!------------------------------------------------------------------------------ 5664 TYPE(Nodes_t), SAVE :: Nodes 5665 TYPE(Element_t), POINTER :: ElementCopy 5666 TYPE(GaussIntegrationPoints_t) :: IP 5667 LOGICAL :: SecondKindBasis, stat, ElementCopyCreated 5668 INTEGER :: TetraFaceMap(4,3), ActiveFaceMap(3) 5669 INTEGER :: DOFs, i, p 5670 REAL(KIND=dp) :: VLoad(3,n), VL(3), Normal(3), Basis(n), DetJ, s 5671 REAL(KIND=dp) :: f(3), u, v 5672!------------------------------------------------------------------------------ 5673 IF (GetElementFamily(Parent) /= 5) & 5674 CALL Fatal('FaceElementDOFs','A tetrahedral parent element supposed') 5675 5676 DOFs = 1 5677 IF (PRESENT(FDOFs)) THEN 5678 IF (FDOFs > 3) THEN 5679 CALL Fatal('FaceElementDOFs','Cannot yet handle more than 3 DOFs per face') 5680 ELSE 5681 DOFs = FDOFs 5682 END IF 5683 END IF 5684 5685 IF (PRESENT(SecondFamily)) THEN 5686 SecondKindBasis = SecondFamily 5687 IF (SecondKindBasis .AND. (DOFs /= 3) ) & 5688 CALL Fatal('FaceElementDOFs','3 DOFs per face expected') 5689 ELSE 5690 SecondKindBasis = .FALSE. 5691 END IF 5692 IF (.NOT. SecondKindBasis .AND. DOFs > 1) & 5693 CALL Fatal('FaceElementDOFs','An unexpected DOFs count per face') 5694 5695 ElementCopyCreated= .FALSE. 5696 IF (SecondKindBasis) THEN 5697 TetraFaceMap(1,:) = [ 2, 1, 3 ] 5698 TetraFaceMap(2,:) = [ 1, 2, 4 ] 5699 TetraFaceMap(3,:) = [ 2, 3, 4 ] 5700 TetraFaceMap(4,:) = [ 3, 1, 4 ] 5701 5702 ActiveFaceMap(1:3) = TetraFaceMap(FaceId,1:3) 5703 5704 IF (ANY(Element % NodeIndexes(1:3) /= Parent % NodeIndexes(ActiveFaceMap(1:3)))) THEN 5705 ! 5706 ! The parent element face is indexed differently than the boundary element. 5707 ! Create a copy of the boundary element which is indexed as the parent element 5708 ! face so that we can return the values of DOFs in the default order. 5709 ! Reordering is supposed to be done outside this subroutine. 5710 ! 5711 ElementCopyCreated = .TRUE. 5712 ElementCopy => AllocateElement() 5713 ElementCopy % Type => Element % Type 5714 ALLOCATE(ElementCopy % NodeIndexes(3)) 5715 ElementCopy % NodeIndexes(1:3) = Parent % NodeIndexes(ActiveFaceMap(1:3)) 5716 ElementCopy % BodyId = Element % BodyId 5717 ElementCopy % BoundaryInfo => Element % BoundaryInfo 5718 ELSE 5719 ElementCopy => Element 5720 END IF 5721 ELSE 5722 ElementCopy => Element 5723 END IF 5724 CALL GetElementNodes(Nodes, ElementCopy) 5725 5726 i = LEN_TRIM(Name) 5727 VLoad(1,1:n) = GetReal(BC, Name(1:i)//' 1', stat, ElementCopy) 5728 VLoad(2,1:n) = GetReal(BC, Name(1:i)//' 2', stat, ElementCopy) 5729 VLoad(3,1:n) = GetReal(BC, Name(1:i)//' 3', stat, ElementCopy) 5730 5731 IP = GaussPoints(ElementCopy, 3) ! Feasible for a triangular face 5732 Integral(:) = 0.0d0 5733 DO p=1,IP % n 5734 stat = ElementInfo(ElementCopy, Nodes, IP % u(p), & 5735 IP % v(p), IP % w(p), DetJ, Basis) 5736 ! 5737 ! We need a normal that points outwards from the parent element. 5738 ! The following function call should be consistent with this goal 5739 ! in the case of a volume-vacuum interface provided a target body 5740 ! for the normal has not been given to blur the situation. 5741 ! TO DO: Modify to allow other scenarios 5742 ! 5743 Normal = NormalVector(ElementCopy, Nodes, IP % u(p), IP % v(p), .TRUE.) 5744 5745 VL = MATMUL(VLoad(:,1:n), Basis(1:n)) 5746 s = IP % s(p) * DetJ 5747 5748 IF (SecondKindBasis) THEN 5749 ! Standard coordinates mapped to the p-element coordinates: 5750 u = -1.0d0 + 2.0d0*IP % u(p) + IP % v(p) 5751 v = sqrt(3.0d0)*IP % v(p) 5752 ! 5753 ! The weight functions for the evaluation of DOFs: 5754 f(1) = sqrt(3.0d0) * 0.5d0 * (1.0d0 - 2.0d0*u + 1.0d0/3.0d0 - 2.0d0/sqrt(3.0d0)*v) 5755 f(2) = sqrt(3.0d0) * 0.5d0 * (1.0d0 + 2.0d0*u + 1.0d0/3.0d0 - 2.0d0/sqrt(3.0d0)*v) 5756 f(3) = sqrt(3.0d0) * (-1.0d0/3.0d0 + 2.0d0/sqrt(3.0d0)*v) 5757 5758 DO i=1,DOFs 5759 Integral(i) = Integral(i) + SUM(VL(1:3) * Normal(1:3)) * f(i) * s 5760 END DO 5761 ELSE 5762 Integral(1) = Integral(1) + SUM(VL(1:3) * Normal(1:3)) * s 5763 END IF 5764 END DO 5765 IF (ElementCopyCreated) DEALLOCATE(ElementCopy % NodeIndexes) 5766!------------------------------------------------------------------------------ 5767 END SUBROUTINE FaceElementDOFs 5768!------------------------------------------------------------------------------ 5769 5770 5771!> In the case of p-approximation, compute the element stiffness matrix and 5772!> force vector in order to assemble a system of equations for approximating 5773!> a given Dirichlet condition 5774!------------------------------------------------------------------------------ 5775 SUBROUTINE LocalBcBDOFs(BC, Element, nd, Name, STIFF, Force ) 5776!------------------------------------------------------------------------------ 5777 5778 IMPLICIT NONE 5779 5780 TYPE(ValueList_t), POINTER :: BC !< The list of boundary condition values 5781 TYPE(Element_t), POINTER :: Element !< The boundary element handled 5782 INTEGER :: nd !< The number of DOFs in the boundary element 5783 CHARACTER(LEN=MAX_NAME_LEN) :: Name !< The name of boundary condition 5784 REAL(KIND=dp) :: STIFF(:,:) !< The element stiffness matrix 5785 REAL(KIND=dp) :: Force(:) !< The element force vector 5786!------------------------------------------------------------------------------ 5787 TYPE(GaussIntegrationPoints_t) :: IP 5788 INTEGER :: p,q,t 5789 REAL(KIND=dp) :: Basis(nd) 5790 REAL(KIND=dp) :: xip,yip,zip,s,DetJ,Load 5791 LOGICAL :: stat 5792 TYPE(Nodes_t) :: Nodes 5793 SAVE Nodes 5794!------------------------------------------------------------------------------ 5795 5796 ! Get nodes of boundary elements parent and gauss points for boundary 5797 CALL GetElementNodes( Nodes, Element ) 5798 IP = GaussPoints( Element ) 5799 5800 FORCE(1:nd) = 0.0d0 5801 STIFF(1:nd,1:nd) = 0.0d0 5802 5803 DO t=1,IP % n 5804 stat = ElementInfo( Element, Nodes, IP % u(t), & 5805 IP % v(t), IP % w(t), DetJ, Basis ) 5806 5807 s = IP % s(t) * DetJ 5808 5809 ! Get value of boundary condition 5810 xip = SUM( Basis(1:nd) * Nodes % x(1:nd) ) 5811 yip = SUM( Basis(1:nd) * Nodes % y(1:nd) ) 5812 zip = SUM( Basis(1:nd) * Nodes % z(1:nd) ) 5813 Load = ListGetConstReal( BC, Name, x=xip,y=yip,z=zip ) 5814 5815 ! Build local stiffness matrix and force vector 5816 DO p=1,nd 5817 DO q=1,nd 5818 STIFF(p,q) = STIFF(p,q) + s * Basis(p)*Basis(q) 5819 END DO 5820 FORCE(p) = FORCE(p) + s * Load * Basis(p) 5821 END DO 5822 END DO 5823!------------------------------------------------------------------------------ 5824 END SUBROUTINE LocalBcBDOFs 5825!------------------------------------------------------------------------------ 5826 5827 5828!------------------------------------------------------------------------------ 5829!> Finished the bulk assembly of the matrix equation. 5830!> Optionally save the matrix for later use. 5831!------------------------------------------------------------------------------ 5832 SUBROUTINE DefaultFinishBulkAssembly( Solver, BulkUpdate ) 5833!------------------------------------------------------------------------------ 5834 TYPE(Solver_t), OPTIONAL, TARGET :: Solver 5835 LOGICAL, OPTIONAL :: BulkUpdate 5836 5837 TYPE(Solver_t), POINTER :: PSolver 5838 TYPE(ValueList_t), POINTER :: Params 5839 LOGICAL :: Bupd, Found 5840 INTEGER :: n 5841 CHARACTER(LEN=MAX_NAME_LEN) :: str 5842 LOGICAL :: Transient 5843 REAL(KIND=dp) :: SScond 5844 INTEGER :: Order 5845 5846 IF( PRESENT( Solver ) ) THEN 5847 PSolver => Solver 5848 ELSE 5849 PSolver => CurrentModel % Solver 5850 END IF 5851 5852 Params => GetSolverParams( PSolver ) 5853 5854 IF( ListGetLogical( Params,'Bulk Assembly Timing',Found ) ) THEN 5855 CALL CheckTimer('BulkAssembly'//GetVarName(PSolver % Variable), Level=5, Delete=.TRUE. ) 5856 END IF 5857 5858 ! Reset colouring 5859 PSolver % CurrentColour = 0 5860 5861 BUpd = .FALSE. 5862 IF ( PRESENT(BulkUpdate) ) THEN 5863 BUpd = BulkUpdate 5864 ELSE 5865 BUpd = GetLogical( Params,'Calculate Loads', Found ) 5866 IF( BUpd ) THEN 5867 str = GetString( Params,'Calculate Loads Slot', Found ) 5868 IF(Found) THEN 5869 BUpd = ( TRIM( str ) == 'bulk assembly') 5870 END IF 5871 END IF 5872 BUpd = BUpd .OR. GetLogical( Params,'Constant Bulk System', Found ) 5873 BUpd = BUpd .OR. GetLogical( Params,'Save Bulk System', Found ) 5874 BUpd = BUpd .OR. GetLogical( Params,'Constant Bulk Matrix', Found ) 5875 BUpd = BUpd .OR. GetLogical( Params,'Constraint Modes Analysis',Found) 5876 BUpd = BUpd .OR. GetLogical( Params,'Control Use Loads',Found ) 5877 END IF 5878 5879 IF( BUpd ) THEN 5880 str = GetString( Params,'Equation',Found) 5881 CALL Info('DefaultFinishBulkAssembly','Saving bulk values for: '//TRIM(str), Level=6 ) 5882 IF( GetLogical( Params,'Constraint Modes Mass Lumping',Found) ) THEN 5883 CALL CopyBulkMatrix( PSolver % Matrix, BulkMass = .TRUE. ) 5884 ELSE 5885 CALL CopyBulkMatrix( PSolver % Matrix, BulkMass = ASSOCIATED(PSolver % Matrix % MassValues), & 5886 BulkDamp = ASSOCIATED(PSolver % Matrix % DampValues) ) 5887 END IF 5888 END IF 5889 5890 IF( GetLogical( Params,'Bulk System Multiply',Found ) ) THEN 5891 CALL Info('DefaultFinishAssembly','Multiplying matrix equation',Level=10) 5892 CALL LinearSystemMultiply( PSolver ) 5893 END IF 5894 5895 IF ( ListGetLogical( Params,'Linear System Save',Found )) THEN 5896 str = GetString( Params,'Linear System Save Slot', Found ) 5897 IF(Found .AND. TRIM( str ) == 'bulk assembly') THEN 5898 CALL SaveLinearSystem( PSolver ) 5899 END IF 5900 END IF 5901 5902 IF( ListGetLogical( Params,'Linear System Remove Zeros',Found ) ) THEN 5903 CALL CRS_RemoveZeros( PSolver % Matrix ) 5904 END IF 5905 5906 IF( ListGetLogical( PSolver % Values,'Boundary Assembly Timing',Found ) ) THEN 5907 CALL ResetTimer('BoundaryAssembly'//GetVarName(PSolver % Variable) ) 5908 END IF 5909 5910 END SUBROUTINE DefaultFinishBulkAssembly 5911 5912 5913!------------------------------------------------------------------------------ 5914!> Finished the boundary assembly of the matrix equation. 5915!> Optionally save the matrix for later use. 5916!------------------------------------------------------------------------------ 5917 SUBROUTINE DefaultFinishBoundaryAssembly( Solver, BulkUpdate ) 5918!------------------------------------------------------------------------------ 5919 TYPE(Solver_t), OPTIONAL, TARGET :: Solver 5920 LOGICAL, OPTIONAL :: BulkUpdate 5921 TYPE(Solver_t), POINTER :: PSolver 5922 TYPE(ValueList_t), POINTER :: Params 5923 LOGICAL :: Bupd, Found 5924 INTEGER :: n 5925 CHARACTER(LEN=MAX_NAME_LEN) :: str 5926 5927 IF( PRESENT( Solver ) ) THEN 5928 PSolver => Solver 5929 ELSE 5930 PSolver => CurrentModel % Solver 5931 END IF 5932 5933 Params => GetSolverParams(PSolver) 5934 5935 IF( ListGetLogical( Params,'Boundary Assembly Timing',Found ) ) THEN 5936 CALL CheckTimer('BoundaryAssembly'//GetVarName(PSolver % Variable), Level=5, Delete=.TRUE. ) 5937 END IF 5938 5939 ! Reset colouring 5940 PSolver % CurrentBoundaryColour = 0 5941 5942 BUpd = .FALSE. 5943 IF ( PRESENT(BulkUpdate) ) THEN 5944 BUpd = BulkUpdate 5945 IF ( .NOT. BUpd ) RETURN 5946 5947 ELSE 5948 BUpd = GetLogical( Params,'Calculate Loads', Found ) 5949 IF( BUpd ) THEN 5950 str = GetString( Params,'Calculate Loads Slot', Found ) 5951 IF(Found) THEN 5952 BUpd = ( TRIM( str ) == 'boundary assembly') 5953 ELSE 5954 BUpd = .FALSE. 5955 END IF 5956 BUpd = BUpd .OR. GetLogical( Params,'Constant System', Found ) 5957 END IF 5958 END IF 5959 5960 IF( BUpd ) THEN 5961 CALL Info('DefaultFinishBoundaryAssembly','Saving system values for Solver: '& 5962 //TRIM(PSolver % Variable % Name), Level=8) 5963 CALL CopyBulkMatrix( PSolver % Matrix ) 5964 END IF 5965 5966 IF ( ListGetLogical( Params,'Linear System Save',Found )) THEN 5967 str = GetString( Params,'Linear System Save Slot', Found ) 5968 IF(Found .AND. TRIM( str ) == 'boundary assembly') THEN 5969 CALL SaveLinearSystem( PSolver ) 5970 END IF 5971 END IF 5972 5973 ! Create contact BCs using mortar conditions. 5974 !--------------------------------------------------------------------- 5975 IF( ListGetLogical( PSolver % Values,'Apply Contact BCs',Found) ) THEN 5976 CALL DetermineContact( PSolver ) 5977 END IF 5978 5979 5980 END SUBROUTINE DefaultFinishBoundaryAssembly 5981 5982 5983 5984!------------------------------------------------------------------------------ 5985!> Finished the assembly of the matrix equation, mainly effects in transient simulation 5986!> Also may be used to set implicit relaxation on the linear system before 5987!> applying Dirichlet conditions. If flux corrected transport is applied then 5988!> make the initial linear system to be of low order. 5989!------------------------------------------------------------------------------ 5990 SUBROUTINE DefaultFinishAssembly( Solver ) 5991!------------------------------------------------------------------------------ 5992 TYPE(Solver_t), OPTIONAL, TARGET :: Solver 5993 5994 INTEGER :: order, n 5995 LOGICAL :: Found, Transient 5996 TYPE(ValueList_t), POINTER :: Params 5997 TYPE(Solver_t), POINTER :: PSolver 5998 TYPE(Matrix_t), POINTER :: A 5999 CHARACTER(LEN=MAX_NAME_LEN) :: str 6000 REAL(KIND=dp) :: sscond 6001 6002 IF( PRESENT( Solver ) ) THEN 6003 PSolver => Solver 6004 ELSE 6005 PSolver => CurrentModel % Solver 6006 END IF 6007 6008 Params => GetSolverParams(PSolver) 6009 6010 ! Nonlinear timestepping needs a copy of the linear system from previous 6011 ! timestep. Hence the saving of the linear system is enforced. 6012 IF( ListGetLogical( Params,'Nonlinear Timestepping', Found ) ) THEN 6013 CALL Info('DefaultFinishAssembly','Saving system values for Solver: '& 6014 //TRIM(PSolver % Variable % Name), Level=8) 6015 CALL CopyBulkMatrix( PSolver % Matrix ) 6016 END IF 6017 6018 ! Makes a low order matrix of the initial one saving original values 6019 ! to BulkValues. Also created a lumped mass matrix. 6020 IF( ListGetLogical( Params,'Linear System FCT',Found ) ) THEN 6021 IF( PSolver % Variable % Dofs == 1 ) THEN 6022 CALL CRS_FCTLowOrder( PSolver % Matrix ) 6023 ELSE 6024 CALL Fatal('DefaultFinishAssembly','FCT scheme implemented only for one dof') 6025 END IF 6026 END IF 6027 6028 IF(GetLogical(Params,'Use Global Mass Matrix',Found)) THEN 6029 6030 Transient = ( ListGetString( CurrentModel % Simulation, 'Simulation Type' ) == 'transient') 6031 IF( Transient ) THEN 6032 SSCond = ListGetCReal( PSolver % Values,'Steady State Condition',Found ) 6033 IF( Found .AND. SSCond > 0.0_dp ) Transient = .FALSE. 6034 END IF 6035 6036 IF( Transient ) THEN 6037 order = GetInteger(Params,'Time Derivative Order',Found) 6038 IF(.NOT.Found) Order = PSolver % TimeOrder 6039 6040 SELECT CASE(order) 6041 6042 CASE(1) 6043 CALL Default1stOrderTimeGlobal(PSolver) 6044 6045 CASE(2) 6046 CALL Default2ndOrderTimeGlobal(PSolver) 6047 END SELECT 6048 END IF 6049 END IF 6050 6051 CALL FinishAssembly( PSolver, PSolver % Matrix % RHS ) 6052 6053 IF( GetLogical( Params,'Linear System Multiply',Found ) ) THEN 6054 CALL Info('DefaultFinishAssembly','Multiplying matrix equation',Level=10) 6055 CALL LinearSystemMultiply( PSolver ) 6056 END IF 6057 6058 IF( ListCheckPrefix( Params,'Linear System Diagonal Min') ) THEN 6059 CALL LinearSystemMinDiagonal( PSolver ) 6060 END IF 6061 6062 6063 IF ( ListGetLogical( Params,'Linear System Save',Found )) THEN 6064 str = GetString( Params,'Linear System Save Slot', Found ) 6065 IF(Found .AND. TRIM( str ) == 'assembly') THEN 6066 CALL SaveLinearSystem( PSolver ) 6067 END IF 6068 END IF 6069 6070!------------------------------------------------------------------------------ 6071 END SUBROUTINE DefaultFinishAssembly 6072!------------------------------------------------------------------------------ 6073 6074 6075 6076!> Returns integration points for edge or face of p element 6077!------------------------------------------------------------------------------ 6078 FUNCTION GaussPointsBoundary(Element, boundary, np) RESULT(gaussP) 6079!------------------------------------------------------------------------------ 6080 USE PElementMaps, ONLY : getElementBoundaryMap 6081 USE Integration 6082 IMPLICIT NONE 6083 6084 ! Parameters 6085 TYPE(Element_t) :: Element 6086 INTEGER, INTENT(IN) :: boundary, np 6087 6088 TYPE( GaussIntegrationPoints_t ) :: gaussP 6089 TYPE(Nodes_t) :: bNodes 6090 TYPE(Element_t) :: mapElement 6091 TYPE(Element_t), POINTER :: RefElement 6092 INTEGER :: i, n, eCode, bMap(4) 6093 REAL(KIND=dp), TARGET :: x(4), y(4), z(4) 6094 REAL(KIND=dp), POINTER CONTIG :: xP(:), yP(:), zP(:) 6095 6096 SELECT CASE(Element % TYPE % ElementCode / 100) 6097 ! Triangle and Quadrilateral 6098 CASE (3,4) 6099 n = 2 6100 eCode = 202 6101 ! Tetrahedron 6102 CASE (5) 6103 n = 3 6104 eCode = 303 6105 ! Pyramid 6106 CASE (6) 6107 ! Select edge element by boundary 6108 IF (boundary == 1) THEN 6109 n = 4 6110 eCode = 404 6111 ELSE 6112 n = 3 6113 eCode = 303 6114 END IF 6115 ! Wedge 6116 CASE (7) 6117 ! Select edge element by boundary 6118 SELECT CASE (boundary) 6119 CASE (1,2) 6120 n = 3 6121 eCode = 303 6122 CASE (3,4,5) 6123 n = 4 6124 eCode = 404 6125 END SELECT 6126 ! Brick 6127 CASE (8) 6128 n = 4 6129 eCode = 404 6130 CASE DEFAULT 6131 WRITE (*,*) 'DefUtils::GaussPointsBoundary: Unsupported element type' 6132 END SELECT 6133 6134 ! Get element boundary map 6135 bMap(1:4) = getElementBoundaryMap(Element, boundary) 6136 ! Get ref nodes for element 6137 xP => x 6138 yP => y 6139 zP => z 6140 CALL GetRefPElementNodes( Element % Type,xP,yP,zP ) 6141 ALLOCATE(bNodes % x(n), bNodes % y(n), bNodes % z(n)) 6142 6143 ! Set coordinate points of destination 6144 DO i=1,n 6145 IF (bMap(i) == 0) CYCLE 6146 bNodes % x(i) = x(bMap(i)) 6147 bNodes % y(i) = y(bMap(i)) 6148 bNodes % z(i) = z(bMap(i)) 6149 END DO 6150 6151 ! Get element to map from 6152 mapElement % TYPE => GetElementType(eCode) 6153 CALL AllocateVector(mapElement % NodeIndexes, mapElement % TYPE % NumberOfNodes) 6154 6155 ! Get gauss points and map them to given element 6156 gaussP = GaussPoints( mapElement, np ) 6157 6158 CALL MapGaussPoints( mapElement, mapElement % TYPE % NumberOfNodes, gaussP, bNodes ) 6159 6160 ! Deallocate memory 6161 DEALLOCATE(bNodes % x, bNodes % y, bNodes % z, MapElement % NodeIndexes) 6162!------------------------------------------------------------------------------ 6163 END FUNCTION GaussPointsBoundary 6164!------------------------------------------------------------------------------ 6165 6166 6167!------------------------------------------------------------------------------ 6168 SUBROUTINE MapGaussPoints( Element, n, gaussP, Nodes ) 6169!------------------------------------------------------------------------------ 6170 IMPLICIT NONE 6171 6172 TYPE(Element_t) :: Element 6173 TYPE(GaussIntegrationPoints_t) :: gaussP 6174 TYPE(Nodes_t) :: Nodes 6175 INTEGER :: n 6176 6177 INTEGER :: i 6178 REAL(KIND=dp) :: xh,yh,zh,sh, DetJ 6179 REAL(KIND=dp) :: Basis(n) 6180 LOGICAL :: stat 6181 6182 ! Map each gauss point from reference element to given nodes 6183 DO i=1,gaussP % n 6184 stat = ElementInfo( Element, Nodes, gaussP % u(i), gaussP % v(i), gaussP % w(i), & 6185 DetJ, Basis ) 6186 6187 IF (.NOT. stat) THEN 6188 CALL Fatal( 'DefUtils::MapGaussPoints', 'Element to map degenerate') 6189 END IF 6190 6191 ! Get mapped points 6192 sh = gaussP % s(i) * DetJ 6193 xh = SUM( Basis(1:n) * Nodes % x(1:n) ) 6194 yh = SUM( Basis(1:n) * Nodes % y(1:n) ) 6195 zh = SUM( Basis(1:n) * Nodes % z(1:n) ) 6196 ! Set mapped points 6197 gaussP % u(i) = xh 6198 gaussP % v(i) = yh 6199 gaussP % w(i) = zh 6200 gaussP % s(i) = sh 6201 END DO 6202!------------------------------------------------------------------------------ 6203 END SUBROUTINE MapGaussPoints 6204!------------------------------------------------------------------------------ 6205 6206!> Calculate global indexes of boundary dofs for given p-element lying on 6207!> a boundary. 6208!------------------------------------------------------------------------------ 6209 SUBROUTINE getBoundaryIndexes( Mesh, Element, Parent, Indexes, indSize ) 6210!------------------------------------------------------------------------------ 6211! 6212! Type(Mesh_t) :: Mesh 6213! INPUT: Finite element mesh containing edges and faces of elements 6214! 6215! Type(Element_t) :: Element 6216! INPUT: Boundary element to get indexes for 6217! 6218! Type(Element_t) :: Parent 6219! INPUT: Parent of boundary element to get indexes for 6220! 6221! INTEGER :: Indexes(:) 6222! OUTPUT: Calculated indexes of boundary element in global system 6223! 6224! INTEGER :: indSize 6225! OUTPUT: Size of created index vector, i.e. how many indexes were created 6226! starting from index 1 6227!------------------------------------------------------------------------------ 6228 IMPLICIT NONE 6229 6230 ! Parameters 6231 TYPE(Mesh_t) :: Mesh 6232 TYPE(Element_t) :: Parent 6233 TYPE(Element_t), POINTER :: Element 6234 INTEGER :: indSize, Indexes(:) 6235 6236 ! Variables 6237 TYPE(Element_t), POINTER :: Edge, Face 6238 INTEGER :: i,j,n 6239 6240 ! Clear indexes 6241 Indexes = 0 6242 n = Element % TYPE % NumberOfNodes 6243 6244 ! Nodal indexes 6245 Indexes(1:n) = Element % NodeIndexes(1:n) 6246 6247 ! Assign rest of indexes if necessary 6248 SELECT CASE(Parent % TYPE % DIMENSION) 6249 CASE (1) 6250 indSize = n 6251 CASE (2) 6252 ! Add index for each bubble dof in edge 6253 DO i=1,Element % BDOFs 6254 n = n+1 6255 6256 IF (SIZE(Indexes) < n) THEN 6257 CALL Warn('DefUtils::getBoundaryIndexes','Not enough space reserved for indexes') 6258 RETURN 6259 END IF 6260 6261 Indexes(n) = Mesh % NumberOfNodes + & 6262 (Parent % EdgeIndexes(Element % PDefs % localNumber)-1) * Mesh % MaxEdgeDOFs + i 6263 END DO 6264 6265 indSize = n 6266 CASE (3) 6267 ! Get boundary face 6268 Face => Mesh % Faces( Parent % FaceIndexes(Element % PDefs % localNumber) ) 6269 6270 ! Add indexes of faces edges 6271 DO i=1, Face % TYPE % NumberOfEdges 6272 Edge => Mesh % Edges( Face % EdgeIndexes(i) ) 6273 6274 ! If edge has no dofs jump to next edge 6275 IF (Edge % BDOFs <= 0) CYCLE 6276 6277 DO j=1,Edge % BDOFs 6278 n = n + 1 6279 6280 IF (SIZE(Indexes) < n) THEN 6281 CALL Warn('DefUtils::getBoundaryIndexes','Not enough space reserved for indexes') 6282 RETURN 6283 END IF 6284 6285 Indexes(n) = Mesh % NumberOfNodes +& 6286 ( Face % EdgeIndexes(i)-1)*Mesh % MaxEdgeDOFs + j 6287 END DO 6288 END DO 6289 6290 ! Add indexes of faces bubbles 6291 DO i=1,Face % BDOFs 6292 n = n + 1 6293 6294 IF (SIZE(Indexes) < n) THEN 6295 CALL Warn('DefUtils::getBoundaryIndexes','Not enough space reserved for indexes') 6296 RETURN 6297 END IF 6298 6299 Indexes(n) = Mesh % NumberOfNodes + & 6300 Mesh % NumberOfEdges * Mesh % MaxEdgeDOFs + & 6301 (Parent % FaceIndexes( Element % PDefs % localNumber )-1) * Mesh % MaxFaceDOFs + i 6302 END DO 6303 6304 indSize = n 6305 CASE DEFAULT 6306 CALL Fatal('DefUtils::getBoundaryIndexes','Unsupported dimension') 6307 END SELECT 6308!------------------------------------------------------------------------------ 6309 END SUBROUTINE getBoundaryIndexes 6310!------------------------------------------------------------------------------ 6311 6312 6313!> Calculate global AND local indexes of boundary dofs for given p-element 6314!> lying on a boundary. 6315!------------------------------------------------------------------------------ 6316 SUBROUTINE getBoundaryIndexesGL( Mesh, Element, BElement, lIndexes, gIndexes, indSize ) 6317!------------------------------------------------------------------------------ 6318! 6319! ARGUMENTS: 6320! 6321! Type(Mesh_t) :: Mesh 6322! INPUT: Finite element mesh containing edges and faces of elements 6323! 6324! Type(Element_t) :: Element 6325! INPUT: Parent of boundary element to get indexes for 6326! 6327! Type(Element_t) :: BElement 6328! INPUT: Boundary element to get indexes for 6329! 6330! INTEGER :: lIndexes(:), gIndexes(:) 6331! OUTPUT: Calculated indexes of boundary element in local and 6332! global system 6333! 6334! INTEGER :: indSize 6335! OUTPUT: Size of created index vector, i.e. how many indexes were created 6336! starting from index 1 6337! 6338!------------------------------------------------------------------------------ 6339 IMPLICIT NONE 6340 6341 ! Parameters 6342 TYPE(Mesh_t) :: Mesh 6343 TYPE(Element_t) :: Element 6344 TYPE(Element_t), POINTER :: BElement 6345 INTEGER :: indSize, lIndexes(:), gIndexes(:) 6346 ! Variables 6347 TYPE(Element_t), POINTER :: Edge, Face 6348 INTEGER :: i,j,k,n,edgeDofSum, faceOffSet, edgeOffSet(12), localBoundary, nNodes, bMap(4), & 6349 faceEdgeMap(4) 6350 LOGICAL :: stat 6351 6352 ! Clear indexes 6353 lIndexes = 0 6354 gIndexes = 0 6355 6356 ! Get boundary map and number of nodes on boundary 6357 localBoundary = BElement % PDefs % localNumber 6358 nNodes = BElement % TYPE % NumberOfNodes 6359 bMap(1:4) = getElementBoundaryMap(Element, localBoundary) 6360 n = nNodes + 1 6361 6362 ! Assign local and global node indexes 6363 lIndexes(1:nNodes) = bMap(1:nNodes) 6364 gIndexes(1:nNodes) = Element % NodeIndexes(lIndexes(1:nNodes)) 6365 6366 ! Assign rest of indexes 6367 SELECT CASE(Element % TYPE % DIMENSION) 6368 CASE (2) 6369 edgeDofSum = Element % TYPE % NumberOfNodes 6370 6371 IF (SIZE(lIndexes) < nNodes + Mesh % MaxEdgeDOFs) THEN 6372 WRITE (*,*) 'DefUtils::getBoundaryIndexes: Not enough space reserved for edge indexes' 6373 RETURN 6374 END IF 6375 6376 DO i=1,Element % TYPE % NumberOfEdges 6377 Edge => Mesh % Edges( Element % EdgeIndexes(i) ) 6378 6379 ! For boundary edge add local and global indexes 6380 IF (localBoundary == i) THEN 6381 DO j=1,Edge % BDOFs 6382 lIndexes(n) = edgeDofSum + j 6383 gIndexes(n) = Mesh % NumberOfNodes + & 6384 (Element % EdgeIndexes(localBoundary)-1) * Mesh % MaxEdgeDOFs + j 6385 n = n+1 6386 END DO 6387 EXIT 6388 END IF 6389 6390 edgeDofSum = edgeDofSum + Edge % BDOFs 6391 END DO 6392 6393 indSize = n - 1 6394 CASE (3) 6395 IF (SIZE(lIndexes) < nNodes + (Mesh % MaxEdgeDOFs * BElement % TYPE % NumberOfEdges) +& 6396 Mesh % MaxFaceDofs) THEN 6397 WRITE (*,*) 'DefUtils::getBoundaryIndexes: Not enough space reserved for edge indexes' 6398 RETURN 6399 END IF 6400 6401 ! Get offsets for each edge 6402 edgeOffSet = 0 6403 faceOffSet = 0 6404 edgeDofSum = 0 6405 DO i=1,Element % TYPE % NumberOfEdges 6406 Edge => Mesh % Edges( Element % EdgeIndexes(i) ) 6407 edgeOffSet(i) = edgeDofSum 6408 edgeDofSum = edgeDofSum + Edge % BDOFs 6409 END DO 6410 6411 ! Get offset for faces 6412 faceOffSet = edgeDofSum 6413 6414 ! Add element edges to local indexes 6415 faceEdgeMap(1:4) = getFaceEdgeMap(Element, localBoundary) 6416 Face => Mesh % Faces( Element % FaceIndexes(localBoundary) ) 6417 DO i=1,Face % TYPE % NumberOfEdges 6418 Edge => Mesh % Edges( Face % EdgeIndexes(i) ) 6419 6420 IF (Edge % BDOFs <= 0) CYCLE 6421 6422 DO j=1,Edge % BDOFs 6423 lIndexes(n) = Element % TYPE % NumberOfNodes + edgeOffSet(faceEdgeMap(i)) + j 6424 gIndexes(n) = Mesh % NumberOfNodes +& 6425 ( Face % EdgeIndexes(i)-1)*Mesh % MaxEdgeDOFs + j 6426 n=n+1 6427 END DO 6428 END DO 6429 6430 DO i=1,Element % TYPE % NumberOfFaces 6431 Face => Mesh % Faces( Element % FaceIndexes(i) ) 6432 6433 IF (Face % BDOFs <= 0) CYCLE 6434 6435 ! For boundary face add local and global indexes 6436 IF (localBoundary == i) THEN 6437 DO j=1,Face % BDOFs 6438 lIndexes(n) = Element % TYPE % NumberOfNodes + faceOffSet + j 6439 gIndexes(n) = Mesh % NumberOfNodes + & 6440 Mesh % NumberOfEdges * Mesh % MaxEdgeDOFs + & 6441 (Element % FaceIndexes(localBoundary)-1) * Mesh % MaxFaceDOFs + j 6442 n=n+1 6443 END DO 6444 EXIT 6445 END IF 6446 6447 faceOffSet = faceOffSet + Face % BDOFs 6448 END DO 6449 6450 indSize = n - 1 6451 END SELECT 6452 END SUBROUTINE getBoundaryIndexesGL 6453 6454 6455 6456!------------------------------------------------------------------------------ 6457 SUBROUTINE GetParentUVW( Element,n,Parent,np,U,V,W,Basis ) 6458!------------------------------------------------------------------------------ 6459 TYPE(Element_t) :: Element, Parent 6460 INTEGER :: n, np 6461 REAL(KIND=dp) :: U,V,W,Basis(:) 6462!------------------------------------------------------------------------------ 6463 INTEGER :: i,j 6464 REAL(KIND=dp), POINTER :: LU(:), LV(:), LW(:) 6465 6466 LU => Parent % TYPE % NodeU 6467 LV => Parent % TYPE % NodeV 6468 LW => Parent % TYPE % NodeW 6469 6470 U = 0.0_dp 6471 V = 0.0_dp 6472 W = 0.0_dp 6473 DO i = 1,n 6474 DO j = 1,np 6475 IF ( Element % NodeIndexes(i) == Parent % NodeIndexes(j) ) THEN 6476 U = U + Basis(i) * LU(j) 6477 V = V + Basis(i) * LV(j) 6478 W = W + Basis(i) * LW(j) 6479 EXIT 6480 END IF 6481 END DO 6482 END DO 6483!------------------------------------------------------------------------------ 6484 END SUBROUTINE GetParentUVW 6485!------------------------------------------------------------------------------ 6486 6487 6488!> Returns flag telling whether Newton linearization is active 6489!------------------------------------------------------------------------------ 6490 FUNCTION GetNewtonActive( USolver ) RESULT( NewtonActive ) 6491 LOGICAL :: NewtonActive 6492 TYPE(Solver_t), OPTIONAL, TARGET :: USolver 6493 6494 IF ( PRESENT( USolver ) ) THEN 6495 NewtonActive = USolver % NewtonActive 6496 ELSE 6497 NewtonActive = CurrentModel % Solver % NewtonActive 6498 END IF 6499 END FUNCTION GetNewtonActive 6500 6501 6502!----------------------------------------------------------------------- 6503!> This routine may be used to terminate the program in the case of an error. 6504!----------------------------------------------------------------------- 6505 SUBROUTINE Assert(Condition, Caller, ErrorMessage) 6506!----------------------------------------------------------------------- 6507 CHARACTER(LEN=*), OPTIONAL :: Caller, ErrorMessage 6508 LOGICAL :: Condition 6509!----------------------------------------------------------------------- 6510 IF ( .NOT. OutputLevelMask(0) ) STOP EXIT_ERROR 6511 6512 IF(Condition) RETURN !Assertion passed 6513 6514 WRITE( Message, '(A)') 'ASSERTION ERROR' 6515 6516 IF(PRESENT(Caller)) THEN 6517 WRITE( Message, '(A,A,A)') TRIM(Message),': ',TRIM(Caller) 6518 END IF 6519 6520 IF(PRESENT(ErrorMessage)) THEN 6521 WRITE( Message, '(A,A,A)') TRIM(Message),': ',TRIM(ErrorMessage) 6522 END IF 6523 6524 WRITE( *, '(A)', ADVANCE='YES' ) Message 6525 6526 !Provide a stack trace if no caller info provided 6527#ifdef __GFORTRAN__ 6528 IF(.NOT.PRESENT(Caller)) CALL BACKTRACE 6529#endif 6530 6531 STOP EXIT_ERROR 6532!----------------------------------------------------------------------- 6533 END SUBROUTINE Assert 6534!----------------------------------------------------------------------- 6535 6536 FUNCTION GetNOFColours(USolver) RESULT( ncolours ) 6537 IMPLICIT NONE 6538 TYPE(Solver_t), TARGET, OPTIONAL :: USolver 6539 INTEGER :: ncolours 6540 6541 ncolours = 1 6542 IF ( PRESENT( USolver ) ) THEN 6543 IF( ASSOCIATED( USolver % ColourIndexList ) ) THEN 6544 ncolours = USolver % ColourIndexList % n 6545 USolver % CurrentColour = 0 6546 END IF 6547 ELSE 6548 IF( ASSOCIATED( CurrentModel % Solver % ColourIndexList ) ) THEN 6549 ncolours = CurrentModel % Solver % ColourIndexList % n 6550 CurrentModel % Solver % CurrentColour = 0 6551 END IF 6552 END IF 6553 6554 CALL Info('GetNOFColours','Number of colours: '//TRIM(I2S(ncolours)),Level=12) 6555 END FUNCTION GetNOFColours 6556 6557 FUNCTION GetNOFBoundaryColours(USolver) RESULT( ncolours ) 6558 IMPLICIT NONE 6559 TYPE(Solver_t), TARGET, OPTIONAL :: USolver 6560 INTEGER :: ncolours 6561 6562 ncolours = 1 6563 IF ( PRESENT( USolver ) ) THEN 6564 IF( ASSOCIATED( USolver % BoundaryColourIndexList ) ) THEN 6565 ncolours = USolver % BoundaryColourIndexList % n 6566 USolver % CurrentBoundaryColour = 0 6567 END IF 6568 ELSE 6569 IF( ASSOCIATED( CurrentModel % Solver % BoundaryColourIndexList ) ) THEN 6570 ncolours = CurrentModel % Solver % BoundaryColourIndexList % n 6571 CurrentModel % Solver % CurrentBoundaryColour = 0 6572 END IF 6573 END IF 6574 6575 CALL Info('GetNOFBoundaryColours','Number of colours: '//TRIM(I2S(ncolours)),Level=12) 6576 END FUNCTION GetNOFBoundaryColours 6577 6578 ! Check given colourings are valid and see if they are free of race conditions. 6579 SUBROUTINE CheckColourings(Solver) 6580 IMPLICIT NONE 6581 TYPE(Solver_t) :: Solver 6582 6583 TYPE(Mesh_t), POINTER :: Mesh 6584 TYPE(Graph_t), POINTER :: Colours 6585 TYPE(Graph_t), POINTER :: BoundaryColours 6586 6587 TYPE(Element_t), POINTER :: Element 6588 6589 INTEGER, ALLOCATABLE :: Indexes(:), DOFIndexes(:) 6590 INTEGER :: col, elem, belem, NDOF, dof 6591 LOGICAL :: errors 6592 6593 errors = .FALSE. 6594 6595 Mesh => Solver % Mesh 6596 Colours => Solver % ColourIndexList 6597 BoundaryColours => Solver % BoundaryColourIndexList 6598 6599 ! Allocate workspace and initialize it 6600 ALLOCATE(Indexes(MAX(Mesh % NumberOfNodes,& 6601 Mesh % NumberOfBulkElements*Mesh % MaxElementDOFs,& 6602 Mesh%NumberOfBoundaryElements*Mesh % MaxElementDOFs)), & 6603 DOFIndexes(Mesh % MaxElementDOFs)) 6604 Indexes = 0 6605 6606 ! Check that every element has a colour 6607 IF (ASSOCIATED(Colours)) THEN 6608 DO col=1,Colours % N 6609 DO elem=Colours % Ptr(col), Colours%Ptr(col+1)-1 6610 Indexes(Colours % Ind(elem))=Indexes(Colours % Ind(elem))+1 6611 END DO 6612 END DO 6613 DO elem=1,Mesh % NumberOfBulkElements 6614 IF (Indexes(elem) < 1 .OR. Indexes(elem) > 1) THEN 6615 CALL Warn('CheckColourings','Element not colored correctly: '//i2s(elem)) 6616 errors = .TRUE. 6617 END IF 6618 END DO 6619 6620 Indexes = 0 6621 ! Check that colouring is free of race conditions 6622 DO col=1,Colours % N 6623 DO elem=Colours % Ptr(col), Colours%Ptr(col+1)-1 6624 Element => Mesh % Elements(Colours % Ind(elem)) 6625 NDOF = GetElementDOFs( DOFIndexes, Element, Solver ) 6626 DO dof=1,NDOF 6627 Indexes(DOFIndexes(dof))=Indexes(DOFIndexes(dof))+1 6628 END DO 6629 END DO 6630 ! Check colouring 6631 DO dof=1,Mesh % NumberOfBulkElements*Mesh % MaxElementDOFs 6632 IF (Indexes(dof)>1) THEN 6633 CALL Warn('CheckColourings','DOF not colored correctly: '//i2s(dof)) 6634 errors = .TRUE. 6635 END IF 6636 Indexes(dof)=0 6637 END DO 6638 END DO 6639 END IF 6640 6641 ! Check that every boundary element has a colour 6642 IF (ASSOCIATED(BoundaryColours)) THEN 6643 6644 DO col=1,BoundaryColours % N 6645 DO elem=BoundaryColours % Ptr(col), BoundaryColours%Ptr(col+1)-1 6646 Indexes(BoundaryColours % Ind(elem))=Indexes(BoundaryColours % Ind(elem))+1 6647 END DO 6648 END DO 6649 DO elem=Mesh % NumberOfBulkElements+1,Mesh % NumberOfBulkElements + Mesh % NumberOfBoundaryElements 6650 belem = elem - Mesh % NumberOfBulkElements 6651 IF (Indexes(belem) < 1 .OR. Indexes(belem) > 1) THEN 6652 CALL Warn('CheckColourings','Boundary element not colored correctly: '//i2s(belem)) 6653 errors = .TRUE. 6654 END IF 6655 END DO 6656 6657 Indexes = 0 6658 ! Check that colouring is free of race conditions 6659 DO col=1,BoundaryColours % N 6660 DO elem=BoundaryColours % Ptr(col), BoundaryColours%Ptr(col+1)-1 6661 Element => Mesh % Elements(Mesh % NumberOfBulkElements + BoundaryColours % Ind(elem)) 6662 NDOF = GetElementDOFs( DOFIndexes, Element, Solver, NotDG=.TRUE. ) 6663 ! WRITE (*,'(2(A,I0))') 'BELEM=', elem, ', CMAP=', BoundaryColours % Ind(elem) 6664 ! WRITE (*,*) DOFIndexes(1:NDOF) 6665 DO dof=1,NDOF 6666 Indexes(DOFIndexes(dof))=Indexes(DOFIndexes(dof))+1 6667 ! WRITE (*,'(4(A,I0))') 'EID=', Element % ElementIndex,', dof=', dof, & 6668 ! ', ind=', DOFIndexes(dof), ', colour=', col 6669 END DO 6670 END DO 6671 ! Check colouring 6672 DO dof=1,Mesh % NumberOfBulkElements*Mesh % MaxElementDOFs 6673 IF (Indexes(dof)>1) THEN 6674 CALL Warn('CheckColourings','Boundary DOF not colored correctly: '//i2s(dof)) 6675 errors = .TRUE. 6676 END IF 6677 Indexes(dof)=0 6678 END DO 6679 END DO 6680 END IF 6681 6682 IF (errors) THEN 6683 CALL Warn('CheckColourings','Mesh colouring contained errors') 6684 END IF 6685 6686 DEALLOCATE(Indexes, DOFIndexes) 6687 END SUBROUTINE CheckColourings 6688 6689 6690 6691!------------------------------------------------------------------------------ 6692!> Assemble coupling matrices related to shell-solid interaction by 6693!> utilizing the director data of the shell model. 6694!> A possible scenario is that the diagonal blocks are the matrices of the 6695!> solvers listed using the keyword "Block Solvers". The (1,1)-block is then 6696!> tied up with the value of the first entry in the "Block Solvers" array. 6697!> Here it is assumed that the (2,2)-block is a shell stiffness matrix. 6698!> NOTE: This is still under construction and doesn't couple forces from 6699!> the shell model to the solid model. 6700!------------------------------------------------------------------------------ 6701 SUBROUTINE StructureCouplingAssembly_defutils( Solver, FVar, SVar, A_f, A_s, A_fs, A_sf, & 6702 IsSolid, IsPlate, IsShell, IsBeam ) 6703!------------------------------------------------------------------------------ 6704 TYPE(Solver_t) :: Solver !< The leading solver defining block structure 6705 TYPE(Variable_t), POINTER :: FVar !< Slave structure variable 6706 TYPE(Variable_t), POINTER :: SVar !< Master structure variable 6707 TYPE(Matrix_t), POINTER :: A_f !< (2,2)-block for the "slave" variable 6708 TYPE(Matrix_t), POINTER :: A_s !< (1,1)-block for the "master" variable 6709 TYPE(Matrix_t), POINTER :: A_fs !< (2,1)-block for interaction 6710 TYPE(Matrix_t), POINTER :: A_sf !< (1,2)-block for interaction 6711 LOGICAL :: IsSolid, IsPlate, IsShell, IsBeam !< The type of the slave variable 6712 !------------------------------------------------------------------------------ 6713 LOGICAL, POINTER :: ConstrainedF(:), ConstrainedS(:) 6714 INTEGER, POINTER :: FPerm(:), SPerm(:) 6715 INTEGER :: FDofs, SDofs 6716 TYPE(Mesh_t), POINTER :: Mesh 6717 INTEGER :: i,j,k,jf,js,kf,ks,nf,ns,dim,ncount 6718 REAL(KIND=dp) :: vdiag 6719 !------------------------------------------------------------------------------ 6720 6721 CALL Info('StructureCouplingAssembly','Creating coupling matrix for structures',Level=6) 6722 6723 Mesh => Solver % Mesh 6724 dim = Mesh % MeshDim 6725 6726 ! S refers to the first and F to the second block (was fluid): 6727 FPerm => FVar % Perm 6728 SPerm => SVar % Perm 6729 6730 fdofs = FVar % Dofs 6731 sdofs = SVar % Dofs 6732 6733 IF( IsSolid ) CALL Info('StructureCouplingAssembly','Assuming coupling with solid solver',Level=8) 6734 IF( IsBeam ) CALL Info('StructureCouplingAssembly','Assuming coupling with beam solver',Level=8) 6735 IF( IsPlate ) CALL Info('StructureCouplingAssembly','Assuming coupling with plate solver',Level=8) 6736 IF( IsShell ) CALL Info('StructureCouplingAssembly','Assuming coupling with shell solver',Level=8) 6737 6738 ConstrainedF => A_f % ConstrainedDof 6739 ConstrainedS => A_s % ConstrainedDof 6740 6741 nf = SIZE( FVar % Values ) 6742 ns = SIZE( SVar % Values ) 6743 6744 CALL Info('StructureCouplingAssembly','Slave structure dofs '//TRIM(I2S(nf))//& 6745 ' with '//TRIM(I2S(fdofs))//' components',Level=10) 6746 CALL Info('StructureCouplingAssembly','Master structure dofs '//TRIM(I2S(ns))//& 6747 ' with '//TRIM(I2S(sdofs))//' components',Level=10) 6748 CALL Info('StructureCouplingAssembly','Assuming '//TRIM(I2S(dim))//& 6749 ' active spatial dimensions',Level=10) 6750 6751 IF( A_fs % FORMAT == MATRIX_LIST ) THEN 6752 ! Add the largest entry that allocates the whole list matrix structure 6753 CALL AddToMatrixElement(A_fs,nf,ns,0.0_dp) 6754 CALL AddToMatrixElement(A_sf,ns,nf,0.0_dp) 6755 ELSE 6756 ! If we are revisiting then initialize the CRS matrices to zero 6757 A_fs % Values = 0.0_dp 6758 A_sf % Values = 0.0_dp 6759 END IF 6760 6761 6762 IF (IsShell) THEN 6763 ! 6764 ! The (2,2)-block is a shell matrix. The (2,1)-block should define 6765 ! Dirichlet constraints and the (1,2)-block should apply forces for the 6766 ! (1,1)-block. 6767 ! 6768 ! The following block tries to mimic the functionality of the subroutine 6769 ! SetSolidCouplingBCs in ShellSolver.F90 6770 ! 6771 SetSolidCouplingBCs: BLOCK 6772 TYPE(Variable_t), POINTER :: Displacement 6773 TYPE(Matrix_t), POINTER :: ShellMatrix, A 6774 TYPE(ValueList_t), POINTER :: ValueList 6775 TYPE(Element_t), POINTER :: Element 6776 6777 INTEGER, ALLOCATABLE, TARGET :: BoundaryNodes(:) 6778 INTEGER, ALLOCATABLE :: NearNodes(:) 6779 INTEGER, POINTER :: Perm(:), NodeIndices(:) 6780 INTEGER, POINTER :: Cols(:), Rows(:), Diag(:) 6781 INTEGER :: TargetCount, TargetNode, TargetInd, Row, ShellDOFs, DOFs 6782 INTEGER :: i, j, k, l, n, p, jz, lz, np, i0 6783 INTEGER :: ju, jv, ku, kv 6784 6785 REAL(KIND=dp), ALLOCATABLE :: NearCoordinates(:,:), AllDirectors(:,:) 6786 REAL(KIND=dp), POINTER :: DirectorValues(:) 6787 REAL(KIND=dp) :: res_z, maxres_z, minres_z, h_eff 6788 REAL(KIND=dp) :: d(3), e3(3), d_h(3), v(3) 6789 6790 IF (.NOT. SVar % DOFs == 3) CALL Fatal('StructureCouplingAssembly', & 6791 'Solid-shell coupling possible in 3D only') 6792 DOFs = 3 6793 6794 ShellMatrix => A_f 6795 ShellDOFs = FVar % DOFs 6796 IF (.NOT. ALLOCATED(ShellMatrix % ConstrainedDOF)) & 6797 ALLOCATE(ShellMatrix % ConstrainedDOF(ShellMatrix % NumberOfRows)) 6798 6799 A => A_s 6800 Diag => A % Diag 6801 Rows => A % Rows 6802 Cols => A % Cols 6803 Perm => SVar % Perm 6804 6805 IF (.NOT. ASSOCIATED(A % InvPerm)) THEN 6806 ALLOCATE(A % InvPerm(A % NumberOfRows)) 6807 DO i = 1,SIZE(Perm) 6808 IF (Perm(i) > 0) THEN 6809 A % InvPerm(Perm(i)) = i 6810 END IF 6811 END DO 6812 END IF 6813 6814 ! --------------------------------------------------------- 6815 ! Count nodes where the coupling will be activated and 6816 ! allocate arrays for saving the directors at these nodes: 6817 ! --------------------------------------------------------- 6818 p = 0 6819 DO i=1,Mesh % NumberOfNodes 6820 jf = FPerm(i) 6821 js = SPerm(i) 6822 IF (jf > 0 .AND. js > 0 ) p = p + 1 6823 END DO 6824 6825 ALLOCATE(BoundaryNodes(p)) 6826 ALLOCATE(AllDirectors(3,p)) 6827 BoundaryNodes = 0 6828 6829 ! ----------------------------------------------------------- 6830 ! Try to figure out the director from the shell element data: 6831 ! ----------------------------------------------------------- 6832 l = 0 6833 DO K=1,Mesh % NumberOfBulkElements 6834 Element => Mesh % Elements(K) 6835 NodeIndices => Element % NodeIndexes 6836 n = Element % TYPE % NumberOfNodes 6837 ! 6838 ! Proceed with shell elements only 6839 ! 6840 IF (ANY(FPerm(NodeIndices(1:n)) == 0)) CYCLE 6841 6842 DirectorValues => NULL() 6843 DirectorValues => GetElementalDirectorInt(Mesh,Element) 6844 6845 IF (.NOT. ASSOCIATED(DirectorValues)) THEN 6846 CALL Fatal('StructureCouplingAssembly', & 6847 'Director cannot be found from shell elements') 6848 ! ELSE 6849 ! PRINT *, 'ok, DIRECTOR DATA FOUND FOR ELEMENT = ', K 6850 END IF 6851 6852 !print *, 'Nodes are = ', NodeIndices(1:n) 6853 DO i=1,n 6854 IF (SPerm(NodeIndices(i)) > 0) THEN 6855 ! This is a common node. If the node hasn't yet been listed, 6856 ! do it now: 6857 IF (ANY(BoundaryNodes(:) == NodeIndices(i))) THEN 6858 ! PRINT *, 'Skipping already listed node' 6859 CYCLE 6860 END IF 6861 l = l + 1 6862 BoundaryNodes(l) = NodeIndices(i) 6863 i0 = 3*(i-1) 6864 AllDirectors(1:3,l) = DirectorValues(i0+1:i0+3) 6865 !ELSE 6866 ! PRINT *, 'This node is not shared with solid, index = ', NodeIndices(i) 6867 END IF 6868 END DO 6869 6870 END DO 6871 NodeIndices => BoundaryNodes(:) 6872 TargetCount = l 6873 6874 IF (TargetCount /= SIZE(BoundaryNodes)) CALL Fatal('StructureCouplingAssembly', & 6875 'Error in retrieving director on solid-shell interface') 6876 6877 ncount = 0 6878 6879 Generate_Dirichlet_Block: DO p=1,TargetCount 6880 TargetNode = NodeIndices(p) 6881 TargetInd = Perm(TargetNode) 6882 ! IF (TargetInd == 0) CYCLE 6883 !------------------------------------------------------------------------------ 6884 ! Find nodes which can potentially be used to calculate the normal derivative 6885 ! of the 3-D solution: 6886 !------------------------------------------------------------------------------ 6887 Row = TargetInd * DOFs 6888 n = (Rows(Row+1)-1 - Rows(Row)-Dofs+1)/DOFs + 1 6889 ALLOCATE(NearNodes(n), NearCoordinates(3,n)) 6890 6891 k = 0 6892 DO i = Rows(Row)+Dofs-1, Rows(Row+1)-1, Dofs 6893 j = Cols(i)/Dofs 6894 k = k + 1 6895 NearNodes(k) = A % InvPerm(j) 6896 END DO 6897 ! PRINT *, 'POTENTIAL NODE CONNECTIONS:' 6898 ! print *, 'Nodes near target=', NearNodes(1:k) 6899 6900 ! 6901 ! The position vectors for the potential nodes: 6902 ! 6903 NearCoordinates(1,1:n) = Mesh % Nodes % x(NearNodes(1:n)) - Mesh % Nodes % x(TargetNode) 6904 NearCoordinates(2,1:n) = Mesh % Nodes % y(NearNodes(1:n)) - Mesh % Nodes % y(TargetNode) 6905 NearCoordinates(3,1:n) = Mesh % Nodes % z(NearNodes(1:n)) - Mesh % Nodes % z(TargetNode) 6906 6907 d = AllDirectors(:,p) 6908 e3 = d/SQRT(DOT_PRODUCT(d,d)) 6909 !------------------------------------------------------------------------------ 6910 ! Seek for nodes which are closest to be parallel to d and have a non-negligible 6911 ! component with respect to d 6912 !------------------------------------------------------------------------------ 6913 maxres_z = 0.0d0 6914 minres_z = 0.0d0 6915 jz = 0 6916 lz = 0 6917 DO i=1,n 6918 IF (NearNodes(i) == TargetNode) CYCLE 6919 6920 res_z = DOT_PRODUCT(e3(:), NearCoordinates(:,i)) / & 6921 SQRT(DOT_PRODUCT(NearCoordinates(:,i), NearCoordinates(:,i))) 6922 ! 6923 ! Skip nearly orthogonal couplings: 6924 ! 6925 IF (ABS(res_z) < 2.0d-2) CYCLE 6926 6927 IF (res_z > 0.0d0) THEN 6928 ! 6929 ! A near node is on +d side 6930 ! 6931 IF (res_z > maxres_z) THEN 6932 jz = NearNodes(i) 6933 maxres_z = res_z 6934 END IF 6935 ELSE 6936 ! 6937 ! A near node is on -d side 6938 ! 6939 IF (res_z < minres_z) THEN 6940 lz = NearNodes(i) 6941 minres_z = res_z 6942 END IF 6943 END IF 6944 END DO 6945 6946 IF (jz == 0) jz = TargetNode 6947 IF (lz == 0) lz = TargetNode 6948 IF (jz == lz) CALL Fatal('StructureCouplingAssembly', & 6949 'No solid nodes to span the director') 6950 6951 6952 !PRINT *, 'HANDLING NODE = ', TargetNode 6953 !PRINT *, 'UPPER NODE = ', JZ 6954 !PRINT *, 'LOWER NODE = ', LZ 6955 6956 ! Now, evaluate the directional derivative DNU(:) in the normal direction: 6957! i = Perm(lz) 6958! j = Perm(jz) 6959! k = Perm(TargetNode) 6960! U_lower(1:3) = SVar % Values(i*DOFs-2:i*DOFs) 6961! U_upper(1:3) = SVar % Values(j*DOFs-2:j*DOFs) 6962! U_mid(1:3) = SVar % Values(k*DOFs-2:k*DOFs) 6963 6964 v(1:3) = [Mesh % Nodes % x(jz) - Mesh % Nodes % x(lz), & 6965 Mesh % Nodes % y(jz) - Mesh % Nodes % y(lz), & 6966 Mesh % Nodes % z(jz) - Mesh % Nodes % z(lz)] 6967 h_eff = SQRT(DOT_PRODUCT(v,v)) 6968! DNU(:) = -1.0d0/h_eff * (U_upper(:) - U_lower(:)) 6969 6970 d_h = v/SQRT(DOT_PRODUCT(v,v)) 6971 IF (ABS(DOT_PRODUCT(d_h,e3)) < 0.98d0) THEN 6972 CALL Warn('StructureCouplingAssembly', & 6973 'A coupling omitted: Solid-model nodes does not span the director') 6974 CYCLE 6975 END IF 6976 6977 ! 6978 ! Finally, constrain the shell to follow the deformation of the solid: 6979 ! 6980 jv = FPerm(TargetNode) 6981 DO j=1,DOFs 6982 ju = SPerm(TargetNode) 6983 ku = sdofs*(ju-1)+j 6984 kv = fdofs*(jv-1)+j 6985 6986 DO k = A_f % Rows(kv),A_f % Rows(kv+1)-1 6987 IF (.NOT. ConstrainedF(ku)) THEN 6988 ! 6989 ! TO DO: Add shell forces to solid 6990 ! 6991 END IF 6992 ! 6993 ! Erase matrix values as a preparation for setting Dirichlet constraints: 6994 ! 6995 A_f % Values(k) = 0.0_dp 6996 END DO 6997 ! 6998 ! Create a Dirichlet constraint to make the shell translation to 6999 ! follow the translation of the solid: 7000 ! 7001 A_f % rhs(kv) = 0.0_dp 7002 A_f % Values(A_f % Diag(kv)) = 1.0_dp 7003 CALL AddToMatrixElement(A_fs, kv, ku, -1.0_dp) 7004 7005 !--------------------- 7006 ! The rotational DOFs 7007 !--------------------- 7008 kv = fdofs*(jv-1)+3+j 7009 DO k = A_f % Rows(kv),A_f % Rows(kv+1)-1 7010 ! 7011 ! TO DO: Add shell moments to solid 7012 ! 7013 7014 ! 7015 ! Erase values as a preparation for setting Dirichlet constraints: 7016 ! 7017 A_f % Values(k) = 0.0_dp 7018 END DO 7019 ! 7020 ! Create a Dirichlet constraint to make the shell rotation to 7021 ! follow the deformation of the solid: 7022 ! 7023 A_f % rhs(kv) = 0.0_dp 7024 A_f % Values(A_f % Diag(kv)) = 1.0_dp 7025 ju = SPerm(lz) 7026 ku = sdofs*(ju-1)+j 7027 CALL AddToMatrixElement(A_fs, kv, ku, -1.0d0/h_eff) 7028 ju = SPerm(jz) 7029 ku = sdofs*(ju-1)+j 7030 CALL AddToMatrixElement(A_fs, kv, ku, 1.0d0/h_eff) 7031 END DO 7032 7033 DEALLOCATE(NearNodes, NearCoordinates) 7034 ncount = ncount + 1 7035 END DO Generate_Dirichlet_Block 7036 7037 IF (TargetCount /= ncount) CALL Fatal('StructureCouplingAssembly', & 7038 'Constraint setting fails for some nodes') 7039 7040 IF (ALLOCATED(AllDirectors)) DEALLOCATE(AllDirectors) 7041 IF (ALLOCATED(BoundaryNodes)) DEALLOCATE(BoundaryNodes) 7042 7043 END BLOCK SetSolidCouplingBCs 7044 7045 ELSE 7046 CALL Fatal('StructureCouplingAssembly','Coupling type not implemented yet!') 7047 END IF 7048 7049 IF( A_fs % FORMAT == MATRIX_LIST ) THEN 7050 CALL List_toCRSMatrix(A_fs) 7051 CALL List_toCRSMatrix(A_sf) 7052 END IF 7053 7054 !PRINT *,'interface fs sum:',SUM(A_fs % Values), SUM( ABS( A_fs % Values ) ) 7055 !PRINT *,'interface sf sum:',SUM(A_sf % Values), SUM( ABS( A_sf % Values ) ) 7056 7057 CALL Info('StructureCouplingAssembly','Number of nodes on interface: '& 7058 //TRIM(I2S(ncount)),Level=10) 7059 CALL Info('StructureCouplingAssembly','Number of entries in slave-master coupling matrix: '& 7060 //TRIM(I2S(SIZE(A_fs % Values))),Level=10) 7061 CALL Info('StructureCouplingAssembly','Number of entries in master-slave coupling matrix: '& 7062 //TRIM(I2S(SIZE(A_sf % Values))),Level=10) 7063 7064 CALL Info('StructureCouplingAssembly','All done',Level=20) 7065 7066 7067!------------------------------------------------------------------------------- 7068 END SUBROUTINE StructureCouplingAssembly_defutils 7069!------------------------------------------------------------------------------- 7070 7071END MODULE DefUtils 7072 7073!> \} // end of subgroup 7074!> \} // end of group 7075 7076