1!/*****************************************************************************/ 2! * 3! * Elmer/Ice, a glaciological add-on to Elmer 4! * http://elmerice.elmerfem.org 5! * 6! * 7! * This program is free software; you can redistribute it and/or 8! * modify it under the terms of the GNU General Public License 9! * as published by the Free Software Foundation; either version 2 10! * of the License, or (at your option) any later version. 11! * 12! * This program is distributed in the hope that it will be useful, 13! * but WITHOUT ANY WARRANTY; without even the implied warranty of 14! * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 15! * GNU General Public License for more details. 16! * 17! * You should have received a copy of the GNU General Public License 18! * along with this program (in file fem/GPL-2); if not, write to the 19! * Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, 20! * Boston, MA 02110-1301, USA. 21! * 22! *****************************************************************************/ 23! ****************************************************************************** 24! * 25! * Authors: Juha Ruokolainen, Olivier Gagliardini, Fabien Gillet-Chaulet 26! * Email: Juha.Ruokolainen@csc.fi 27! * Web: http://elmerice.elmerfem.org 28! * Address: CSC - IT Center for Science Ltd. 29! * Keilaranta 14 30! * 02101 Espoo, Finland 31! * 32! * Date of modification: 03/05 33! * 34! *****************************************************************************/ 35!> Solver for fabric parameter equations 36!------------------------------------------------------------------------------ 37 RECURSIVE SUBROUTINE FabricSolver( Model,Solver,dt,TransientSimulation ) 38!------------------------------------------------------------------------------ 39 40 USE DefUtils 41 42 IMPLICIT NONE 43!------------------------------------------------------------------------------ 44!****************************************************************************** 45! 46! Solve stress equations for one timestep 47! 48! ARGUMENTS: 49! 50! TYPE(Model_t) :: Model, 51! INPUT: All model information (mesh,materials,BCs,etc...) 52! 53! TYPE(Solver_t) :: Solver 54! INPUT: Linear equation solver options 55! 56! REAL(KIND=dp) :: dt, 57! INPUT: Timestep size for time dependent simulations (NOTE: Not used 58! currently) 59! 60!****************************************************************************** 61 62 TYPE(Model_t) :: Model 63 TYPE(Solver_t), TARGET :: Solver 64 65 LOGICAL :: TransientSimulation 66 REAL(KIND=dp) :: dt 67!------------------------------------------------------------------------------ 68! Local variables 69!------------------------------------------------------------------------------ 70 TYPE(Solver_t), POINTER :: PSolver 71 72 TYPE(Matrix_t),POINTER :: StiffMatrix 73 74 INTEGER :: dim,n1,n2,i,j,k,l,n,t,iter,NDeg,STDOFs,LocalNodes,istat 75 76 TYPE(ValueList_t),POINTER :: Material, BC 77 TYPE(Nodes_t) :: ElementNodes 78 TYPE(Element_t),POINTER :: CurrentElement, Element, & 79 ParentElement, LeftParent, RightParent, Edge 80 81 REAL(KIND=dp) :: RelativeChange,UNorm,PrevUNorm,Gravity(3), & 82 Tdiff,Normal(3),NewtonTol,NonlinearTol,s,Wn(7) 83 84 85 INTEGER :: NewtonIter,NonlinearIter 86 87 TYPE(Variable_t), POINTER :: FabricSol, TempSol, FabricVariable, FlowVariable, & 88 EigenFabricVariable, MeshVeloVariable 89 90 REAL(KIND=dp), POINTER :: Temperature(:),Fabric(:), & 91 FabricValues(:), FlowValues(:), EigenFabricValues(:), & 92 MeshVeloValues(:), Solution(:), Ref(:) 93 94 INTEGER, POINTER :: TempPerm(:),FabricPerm(:),NodeIndexes(:), & 95 FlowPerm(:), MeshVeloPerm(:), EigenFabricPerm(:) 96 97 REAL(KIND=dp) :: rho,lambda !Interaction parameter,diffusion parameter 98 REAL(KIND=dp) :: A1plusA2 99 Real(KIND=dp), parameter :: Rad2deg=180._dp/Pi 100 REAL(KIND=dp) :: a2(6) 101 REAL(KIND=dp) :: ai(3), Angle(3) 102 103 LOGICAL :: GotForceBC,GotIt,NewtonLinearization = .FALSE.,UnFoundFatal=.TRUE. 104 105 INTEGER :: body_id,bf_id,eq_id, comp, Indexes(128) 106! 107 INTEGER :: old_body = -1 108 109 REAL(KIND=dp) :: FabricGrid(4879) 110 111 LOGICAL :: AllocationsDone = .FALSE., FreeSurface 112 113 TYPE(Variable_t), POINTER :: TimeVar 114 115 REAL(KIND=dp), ALLOCATABLE:: MASS(:,:), STIFF(:,:), & 116 LocalFluidity(:), LOAD(:,:),Force(:), LocalTemperature(:), & 117 Alpha(:,:),Beta(:), K1(:), K2(:), E1(:), E2(:), E3(:), & 118 Velocity(:,:), MeshVelocity(:,:) 119 120 SAVE MASS, STIFF, LOAD, Force,ElementNodes,Alpha,Beta, & 121 LocalTemperature, LocalFluidity, AllocationsDone, K1, K2, & 122 E1, E2, E3, Wn, FabricGrid, rho, lambda, Velocity, & 123 MeshVelocity, old_body, dim, comp 124!------------------------------------------------------------------------------ 125 CHARACTER(LEN=MAX_NAME_LEN) :: viscosityFile 126 127 REAL(KIND=dp) :: Bu,Bv,Bw,RM(3,3), SaveTime = -1 128 REAL(KIND=dp), POINTER :: PrevFabric(:),CurrFabric(:),TempFabVal(:) 129 130 SAVE ViscosityFile, PrevFabric, CurrFabric,TempFabVal 131#ifdef USE_ISO_C_BINDINGS 132 REAL(KIND=dp) :: at, at0 133#else 134 REAL(KIND=dp) :: at, at0, CPUTime, RealTime 135#endif 136!------------------------------------------------------------------------------ 137 INTERFACE 138 Subroutine R2Ro(a2,dim,ai,angle) 139 USE Types 140 REAL(KIND=dp),intent(in) :: a2(6) 141 Integer :: dim 142 REAL(KIND=dp),intent(out) :: ai(3), Angle(3) 143 End Subroutine R2Ro 144 End Interface 145!------------------------------------------------------------------------------ 146! Read constants from constants section of SIF file 147!------------------------------------------------------------------------------ 148 149 Wn(7) = ListGetConstReal( Model % Constants, 'Gas Constant', GotIt,UnFoundFatal=UnFoundFatal ) 150 !Previous default value: Wn(7) = 8.314 151 WRITE(Message,'(A,F10.4)')'Gas Constant =',Wn(7) 152 CALL INFO('FabricSolve',Message,Level=4) 153!------------------------------------------------------------------------------ 154! Get variables needed for solution 155!------------------------------------------------------------------------------ 156 IF ( .NOT. ASSOCIATED( Solver % Matrix ) ) RETURN 157 158 Solution => Solver % Variable % Values 159 STDOFs = Solver % Variable % DOFs 160 161 FabricSol => VariableGet( Solver % Mesh % Variables, 'Fabric' ) 162 FabricPerm => FabricSol % Perm 163 FabricValues => FabricSol % Values 164 165 TempSol => VariableGet( Solver % Mesh % Variables, 'Temperature' ) 166 IF ( ASSOCIATED( TempSol) ) THEN 167 TempPerm => TempSol % Perm 168 Temperature => TempSol % Values 169 END IF 170 171 FlowVariable => VariableGet( Solver % Mesh % Variables, 'AIFlow' ) 172 IF ( ASSOCIATED( FlowVariable ) ) THEN 173 FlowPerm => FlowVariable % Perm 174 FlowValues => FlowVariable % Values 175 END IF 176 177!!!!! Mesh Velo 178 MeshVeloVariable => VariableGet( Solver % Mesh % Variables, & 179 'Mesh Velocity' ) 180 181 IF ( ASSOCIATED( MeshVeloVariable ) ) THEN 182 MeshVeloPerm => MeshVeloVariable % Perm 183 MeshVeloValues => MeshVeloVariable % Values 184 END IF 185 186 187 StiffMatrix => Solver % Matrix 188 ! UNorm = Solver % Variable % Norm 189 Unorm = SQRT( SUM( FabricValues**2 ) / SIZE(FabricValues) ) 190!------------------------------------------------------------------------------ 191 192!------------------------------------------------------------------------------ 193! Allocate some permanent storage, this is done first time only 194!------------------------------------------------------------------------------ 195 IF ( .NOT. AllocationsDone .OR. Solver % Mesh % Changed) THEN 196 N = Model % MaxElementNodes 197 198 dim = CoordinateSystemDimension() 199 200 IF ( AllocationsDone ) THEN 201 DEALLOCATE( LocalTemperature, & 202 K1,K2,E1,E2,E3, & 203 Force, LocalFluidity, & 204 Velocity,MeshVelocity, & 205 MASS,STIFF, & 206 LOAD, Alpha, Beta, & 207 CurrFabric, TempFabVal ) 208 END IF 209 210 ALLOCATE( LocalTemperature( N ), LocalFluidity( N ), & 211 K1( N ), K2( N ), E1( N ), E2( N ), E3( N ), & 212 Force( 2*STDOFs*N ), & 213 Velocity(4, N ),MeshVelocity(3,N), & 214 MASS( 2*STDOFs*N,2*STDOFs*N ), & 215 STIFF( 2*STDOFs*N,2*STDOFs*N ), & 216 LOAD( 4,N ), Alpha( 3,N ), Beta( N ), & 217 CurrFabric( 5*SIZE(Solver % Variable % Values)), & 218 TempFabVal( SIZE(FabricValues)), & 219 STAT=istat ) 220 221 222 IF ( istat /= 0 ) THEN 223 CALL Fatal( 'FabricSolve', 'Memory allocation error.' ) 224 END IF 225 226 CurrFabric = 0. 227 TempFabVal = 0. 228 IF ( TransientSimulation ) THEN 229 IF (AllocationsDone ) DEALLOCATE (PrevFabric) 230 ALLOCATE( PrevFabric( 5*SIZE(Solver % Variable % Values)) ) 231 PrevFabric = 0. 232 END IF 233 234 DO i=1,Solver % NumberOFActiveElements 235 CurrentElement => GetActiveElement(i) 236 n = GetElementDOFs( Indexes ) 237 n = GetElementNOFNodes() 238 NodeIndexes => CurrentElement % NodeIndexes 239 Indexes(1:n) = Solver % Variable % Perm( Indexes(1:n) ) 240 DO COMP=1,5 241 IF ( TransientSimulation ) THEN 242 PrevFabric(5*(Indexes(1:n)-1)+COMP) = & 243 FabricValues(5*(FabricPerm(NodeIndexes(1:n))-1)+COMP) 244 END IF 245 CurrFabric(5*(Indexes(1:n)-1)+COMP) = & 246 FabricValues(5*(FabricPerm(NodeIndexes(1:n))-1)+COMP) 247 END DO 248 END DO 249 250 AllocationsDone = .TRUE. 251 END IF 252 253 IF( TransientSimulation ) THEN 254 TimeVar => VariableGet( Solver % Mesh % Variables, 'Time' ) 255 IF ( SaveTime /= TimeVar % Values(1) ) THEN 256 SaveTime = TimeVar % Values(1) 257 PrevFabric = CurrFabric 258 END IF 259 END IF 260 261!------------------------------------------------------------------------------ 262! Do some additional initialization, and go for it 263!------------------------------------------------------------------------------ 264 265!------------------------------------------------------------------------------ 266 NonlinearTol = ListGetConstReal( Solver % Values, & 267 'Nonlinear System Convergence Tolerance' ) 268 269 NewtonTol = ListGetConstReal( Solver % Values, & 270 'Nonlinear System Newton After Tolerance' ) 271 272 NewtonIter = ListGetInteger( Solver % Values, & 273 'Nonlinear System Newton After Iterations' ) 274 275 NonlinearIter = ListGetInteger( Solver % Values, & 276 'Nonlinear System Max Iterations',GotIt ) 277 278 IF ( .NOT.GotIt ) NonlinearIter = 1 279!------------------------------------------------------------------------------ 280 281!------------------------------------------------------------------------------ 282 283!------------------------------------------------------------------------------ 284 DO iter=1,NonlinearIter 285!------------------------------------------------------------------------------ 286 at = CPUTime() 287 at0 = RealTime() 288 289 290 CALL Info( 'FabricSolve', ' ', Level=4 ) 291 CALL Info( 'FabricSolve', ' ', Level=4 ) 292 CALL Info( 'FabricSolve', & 293 '-------------------------------------',Level=4 ) 294 WRITE( Message, * ) 'Fabric solver iteration', iter 295 CALL Info( 'FabricSolve', Message,Level=4 ) 296 CALL Info( 'FabricSolve', & 297 '-------------------------------------',Level=4 ) 298 CALL Info( 'FabricSolve', ' ', Level=4 ) 299 CALL Info( 'FabricSolve', 'Starting assembly...',Level=4 ) 300 301!------------------------------------------------------------------------------ 302 303 PrevUNorm = UNorm 304 305 DO COMP=1,2*dim-1 306 307 Solver % Variable % Values = CurrFabric( COMP::5 ) 308 IF ( TransientSimulation ) THEN 309 Solver % Variable % PrevValues(:,1) = PrevFabric( COMP::5 ) 310 END IF 311 312 CALL DefaultInitialize() 313!------------------------------------------------------------------------------ 314 DO t=1,Solver % NumberOFActiveElements 315!------------------------------------------------------------------------------ 316 317 IF ( RealTime() - at0 > 1.0 ) THEN 318 WRITE(Message,'(a,i3,a)' ) ' Assembly: ', INT(100.0 - 100.0 * & 319 (Solver % NumberOfActiveElements-t) / & 320 (1.0*Solver % NumberOfActiveElements)), ' % done' 321 322 CALL Info( 'FabricSolve', Message, Level=5 ) 323 at0 = RealTime() 324 END IF 325 326 CurrentElement => GetActiveElement(t) 327 CALL GetElementNodes( ElementNodes ) 328 n = GetElementDOFs( Indexes ) 329 n = GetElementNOFNodes() 330 NodeIndexes => CurrentElement % NodeIndexes 331 332 333 Material => GetMaterial() 334 body_id = CurrentElement % BodyId 335!------------------------------------------------------------------------------ 336! Read in material constants from Material section 337!------------------------------------------------------------------------------ 338 IF (body_id /= old_body) Then 339 old_body = body_id 340 CALL GetMaterialDefs() 341 END IF 342 343 LocalFluidity(1:n) = ListGetReal( Material, & 344 'Fluidity Parameter', n, NodeIndexes, GotIt,& 345 UnFoundFatal=UnFoundFatal) 346 !Previous default value: LocalFluidity(1:n) = 1.0 347!------------------------------------------------------------------------------ 348! Get element local stiffness & mass matrices 349!------------------------------------------------------------------------------ 350 LocalTemperature = 0.0D0 351 IF ( ASSOCIATED(TempSol) ) THEN 352 DO i=1,n 353 k = TempPerm(NodeIndexes(i)) 354 LocalTemperature(i) = Temperature(k) 355 END DO 356 ELSE 357 LocalTemperature(1:n) = 0.0d0 358 END IF 359 360 K1(1:n) = CurrFabric( 5*(Solver % Variable % Perm(Indexes(1:n))-1)+1 ) 361 K2(1:n) = CurrFabric( 5*(Solver % Variable % Perm(Indexes(1:n))-1)+2 ) 362 E1(1:n) = CurrFabric( 5*(Solver % Variable % Perm(Indexes(1:n))-1)+3 ) 363 E2(1:n) = CurrFabric( 5*(Solver % Variable % Perm(Indexes(1:n))-1)+4 ) 364 E3(1:n) = CurrFabric( 5*(Solver % Variable % Perm(Indexes(1:n))-1)+5 ) 365 366 k = FlowVariable % DOFs 367 Velocity = 0.0d0 368 DO i=1,k 369 Velocity(i,1:n) = FlowValues(k*(FlowPerm(NodeIndexes)-1)+i) 370 END DO 371 372!------------meshvelocity 373 MeshVelocity=0._dp 374 IF (ASSOCIATED(MeshVeloVariable)) Then 375 k = MeshVeloVariable % DOFs 376 DO i=1,k 377 MeshVelocity(i,1:n) = MeshVeloValues(k*(MeshVeloPerm(NodeIndexes)-1)+i) 378 END DO 379 EndIF 380!---------------------------------- 381 382 CALL LocalMatrix( COMP, MASS, STIFF, FORCE, LOAD, K1, K2, E1, & 383 E2, E3, LocalTemperature, LocalFluidity, Velocity, & 384 MeshVelocity, CurrentElement, n, ElementNodes, Wn, rho, lambda ) 385 386!------------------------------------------------------------------------------ 387! Update global matrices from local matrices 388!------------------------------------------------------------------------------ 389 IF ( TransientSimulation ) CALL Default1stOrderTime(MASS,STIFF,FORCE) 390 CALL DefaultUpdateEquations( STIFF, FORCE ) 391!------------------------------------------------------------------------------ 392 END DO 393 CALL Info( 'FabricSolve', 'Assembly done', Level=4 ) 394!------------------------------------------------------------------------------ 395 396!------------------------------------------------------------------------------ 397! Assembly of the edge terms 398!------------------------------------------------------------------------------ 399! 400!3D => Edges => Faces 401 If (dim.eq.3) then 402 DO t=1,Solver % Mesh % NumberOfFaces 403 Edge => Solver % Mesh % Faces(t) 404 IF ( .NOT. ActiveBoundaryElement(Edge) ) CYCLE 405 406 LeftParent => Edge % BoundaryInfo % Left 407 RightParent => Edge % BoundaryInfo % Right 408 IF ( ASSOCIATED(RightParent) .AND. ASSOCIATED(LeftParent) ) THEN 409 n = GetElementNOFNodes( Edge ) 410 n1 = GetElementNOFNodes( LeftParent ) 411 n2 = GetElementNOFNodes( RightParent ) 412 413 k = FlowVariable % DOFs 414 Velocity = 0.0d0 415 DO i=1,k 416 Velocity(i,1:n) = FlowValues(k*(FlowPerm(Edge % NodeIndexes)-1)+i) 417 END DO 418 419 !-------------------mesh velo 420 MeshVelocity=0._dp 421 IF ( ASSOCIATED( MeshVeloVariable ) ) THEN 422 k = MeshVeloVariable % DOFs 423 DO i=1,k 424 MeshVelocity(i,1:n) = MeshVeloValues(k*(MeshVeloPerm(Edge % NodeIndexes)-1)+i) 425 END DO 426 END IF 427 !-------------------------- 428 429 FORCE = 0.0d0 430 MASS = 0.0d0 431 CALL LocalJumps( STIFF,Edge,n,LeftParent,n1,RightParent,n2,Velocity,MeshVelocity ) 432 IF ( TransientSimulation ) CALL Default1stOrderTime(MASS, STIFF, FORCE) 433 CALL DefaultUpdateEquations( STIFF, FORCE, Edge ) 434 END IF 435 END DO 436! 437! 2D 438 Else 439 440 DO t=1,Solver % Mesh % NumberOfEdges 441 Edge => Solver % Mesh % Edges(t) 442 IF ( .NOT. ActiveBoundaryElement(Edge) ) CYCLE 443 444 LeftParent => Edge % BoundaryInfo % Left 445 RightParent => Edge % BoundaryInfo % Right 446 IF ( ASSOCIATED(RightParent) .AND. ASSOCIATED(LeftParent) ) THEN 447 n = GetElementNOFNodes( Edge ) 448 n1 = GetElementNOFNodes( LeftParent ) 449 n2 = GetElementNOFNodes( RightParent ) 450 451 k = FlowVariable % DOFs 452 Velocity = 0.0d0 453 DO i=1,k 454 Velocity(i,1:n) = FlowValues(k*(FlowPerm(Edge % NodeIndexes(1:n))-1)+i) 455 END DO 456 457 !-------------------mesh velo 458 MeshVelocity=0._dp 459 IF ( ASSOCIATED( MeshVeloVariable ) ) THEN 460 k = MeshVeloVariable % DOFs 461 DO i=1,k 462 MeshVelocity(i,1:n) = MeshVeloValues(k*(MeshVeloPerm(Edge % NodeIndexes)-1)+i) 463 END DO 464 END IF 465 !-------------------------- 466 467 FORCE = 0.0d0 468 MASS = 0.0d0 469 CALL LocalJumps( STIFF,Edge,n,LeftParent,n1,RightParent,n2,Velocity,MeshVelocity ) 470 IF ( TransientSimulation ) CALL Default1stOrderTime(MASS, STIFF, FORCE) 471 CALL DefaultUpdateEquations( STIFF, FORCE, Edge ) 472 END IF 473 END DO 474 475 END IF 476 477 CALL DefaultFinishBulkAssembly() 478!------------------------------------------------------------------------------ 479! Loop over the boundary elements 480!------------------------------------------------------------------------------ 481 DO t = 1, Solver % Mesh % NumberOfBoundaryElements 482!------------------------------------------------------------------------------ 483 484 Element => GetBoundaryElement(t) 485 IF( .NOT. ActiveBoundaryElement() ) CYCLE 486 IF( GetElementFamily(Element) == 1 ) CYCLE 487 488 ParentElement => Element % BoundaryInfo % Left 489 IF ( .NOT. ASSOCIATED( ParentElement ) ) & 490 ParentElement => Element % BoundaryInfo % Right 491 492 n = GetElementNOFNodes( Element ) 493 n1 = GetElementNOFnodes( ParentElement ) 494 495 k = FlowVariable % DOFs 496 Velocity = 0.0d0 497 DO i=1,k 498 Velocity(i,1:n) = FlowValues(k*(FlowPerm(Element % NodeIndexes(1:n))-1)+i) 499 END DO 500 501!-------------------mesh velo 502 MeshVelocity=0._dp 503 IF ( ASSOCIATED( MeshVeloVariable ) ) THEN 504 k = MeshVeloVariable % DOFs 505 DO i=1,k 506 MeshVelocity(i,1:n) = MeshVeloValues(k*(MeshVeloPerm(Element % NodeIndexes)-1)+i) 507 End do 508 END IF 509!-------------------------- 510 511 512 BC => GetBC() 513 LOAD = 0.0d0 514 GotIt = .FALSE. 515 IF ( ASSOCIATED(BC) ) THEN 516 LOAD(1,1:n) = GetReal( BC, ComponentName('Fabric', Comp) , GotIt ) 517 END IF 518 519 MASS = 0.0d0 520 CALL LocalMatrixBoundary( STIFF, FORCE, LOAD(1,1:n), & 521 Element, n, ParentElement, n1, Velocity,MeshVelocity, GotIt ) 522 523 IF ( TransientSimulation ) CALL Default1stOrderTime(MASS, STIFF, FORCE) 524 CALL DefaultUpdateEquations( STIFF, FORCE ) 525 END DO 526 527 CALL DefaultFinishAssembly() 528!------------------------------------------------------------------------------ 529 CALL Info( 'FabricSolve', 'Set boundaries done', Level=4 ) 530 531!------------------------------------------------------------------------------ 532! Solve the system and check for convergence 533!------------------------------------------------------------------------------ 534 Unorm = DefaultSolve() 535! CurrFabric( COMP::5 ) = Solver % Variable % Values 536 WRITE(Message,*) 'solve done', minval( solver % variable % values), maxval( Solver % variable % values) 537 CALL Info( 'FabricSolve', Message, Level=4 ) 538 539 n1 = Solver % Mesh % NumberOfNodes 540 ALLOCATE( Ref(n1) ) 541 Ref = 0 542 ! 543 ! FabricValues( COMP::5 ) = 0 !fab sinon remet toutes les 544 ! fabriques � 0 dans le cas ou on a 2 domaines dont l'un a la 545 ! fabrique fixe 546 TempFabVal(COMP::5 ) = 0. !fab 547 548 DO t=1,Solver % NumberOfActiveElements 549 Element => GetActiveElement(t) 550 n = GetElementDOFs( Indexes ) 551 n = GetElementNOFNodes() 552 553 DO i=1,n 554 k = Element % NodeIndexes(i) 555 TempFabVal( 5*(FabricPerm(k)-1) + COMP ) = & 556 TempFabVal( 5*(FabricPerm(k)-1) + COMP ) + & 557 Solver % Variable % Values( Solver % Variable % Perm(Indexes(i)) ) 558 FabricValues( 5*(FabricPerm(k)-1) + COMP ) = & 559 TempFabVal(5*(FabricPerm(k)-1) + COMP ) 560 Ref(k) = Ref(k) + 1 561 END DO 562 END DO 563 564 DO i=1,n1 565 j=FabricPerm(i) 566 IF (j < 1) CYCLE 567 IF ( Ref(i) > 0 ) THEN 568 FabricValues( 5*(j-1)+COMP ) = & 569 FabricValues( 5*(j-1)+COMP ) / Ref(i) 570 END IF 571 END DO 572 573 DEALLOCATE( Ref ) 574 575 SELECT CASE( Comp ) 576 CASE(1) 577 FabricValues( COMP:SIZE(FabricValues):5 ) = & 578 MIN(MAX( FabricValues( COMP:SIZE(FabricValues):5 ) , 0._dp),1._dp) 579 580 CASE(2) 581 FabricValues( COMP:SIZE(FabricValues):5 ) = & 582 MIN(MAX( FabricValues( COMP:SIZE(FabricValues):5 ) , 0._dp),1._dp) 583 584 ! !DO i=1,SIZE(FabricValues),5 585 ! ! IF((FabricValues(i)+FabricValues(i+1)).GT.1._dp) THEN 586 ! ! A1plusA2=FabricValues(i)+FabricValues(i+1) 587 ! ! FabricValues(i)= & 588 ! ! FabricValues(i)/A1plusA2 589 ! ! FabricValues(i+1)= & 590 ! ! FabricValues(i+1)/A1plusA2 591 ! ! END IF 592 ! ! END DO 593 ! 594 595! CASE(3:5) 596! DO i=COMP,SIZE(FabricValues),5 597! IF(FabricValues(i).GT.0._dp) THEN 598! FabricValues(i) = MIN( FabricValues(i) , 0.5_dp) 599! ELSE 600! FabricValues(i) = MAX( FabricValues(i) , -0.5_dp) 601! END IF 602! END DO 603 END SELECT 604 605 END DO ! End DO Comp 606 607 DO i=1,Solver % NumberOFActiveElements 608 CurrentElement => GetActiveElement(i) 609 n = GetElementDOFs( Indexes ) 610 n = GetElementNOFNodes() 611 NodeIndexes => CurrentElement % NodeIndexes 612 Indexes(1:n) = Solver % Variable % Perm( Indexes(1:n) ) 613 DO COMP=1,5 614 CurrFabric(5*(Indexes(1:n)-1)+COMP) = & 615 FabricValues(5*(FabricPerm(NodeIndexes(1:n))-1)+COMP) 616 END DO 617 END DO 618 619 620 Unorm = SQRT( SUM( FabricValues**2 ) / SIZE(FabricValues) ) 621 Solver % Variable % Norm = Unorm 622 623 IF ( PrevUNorm + UNorm /= 0.0d0 ) THEN 624 RelativeChange = 2.0d0 * ABS( PrevUNorm - UNorm) / ( PrevUnorm + UNorm) 625 ELSE 626 RelativeChange = 0.0d0 627 END IF 628 629 WRITE( Message, * ) 'Result Norm : ',UNorm 630 CALL Info( 'FabricSolve', Message, Level=4 ) 631 WRITE( Message, * ) 'Relative Change : ',RelativeChange 632 CALL Info( 'FabricSolve', Message, Level=4 ) 633 634 635!------------------------------------------------------------------------------ 636 IF ( RelativeChange < NewtonTol .OR. & 637 iter > NewtonIter ) NewtonLinearization = .TRUE. 638 639 IF ( RelativeChange < NonLinearTol ) EXIT 640 641!------------------------------------------------------------------------------ 642 EigenFabricVariable => & 643 VariableGet( Solver % Mesh % Variables, 'EigenV' ) 644 IF ( ASSOCIATED( EigenFabricVariable ) ) THEN 645 EigenFabricPerm => EigenFabricVariable % Perm 646 EigenFabricValues => EigenFabricVariable % Values 647 648 DO t=1,Solver % NumberOFActiveElements 649 650 CurrentElement => GetActiveElement(t) 651 n = GetElementNOFNodes() 652 NodeIndexes => CurrentElement % NodeIndexes 653 654 K1(1:n) = FabricValues( 5*(FabricPerm(NodeIndexes(1:n))-1)+1 ) 655 K2(1:n) = FabricValues( 5*(FabricPerm(NodeIndexes(1:n))-1)+2 ) 656 E1(1:n) = FabricValues( 5*(FabricPerm(NodeIndexes(1:n))-1)+3 ) 657 E2(1:n) = FabricValues( 5*(FabricPerm(NodeIndexes(1:n))-1)+4 ) 658 E3(1:n) = FabricValues( 5*(FabricPerm(NodeIndexes(1:n))-1)+5 ) 659 660 Do i=1,n 661 a2(1)=K1(i) 662 a2(2)=K2(i) 663 a2(3)=1._dp-a2(1)-a2(2) 664 a2(4)=E1(i) 665 a2(5)=E2(i) 666 a2(6)=E3(i) 667 668 call R2Ro(a2,dim,ai,angle) 669 670 angle(:)=angle(:)*rad2deg 671 If (angle(1).gt.90._dp) angle(1)=angle(1)-180._dp 672 If (angle(1).lt.-90._dp) angle(1)=angle(1)+180._dp 673 674 EigenFabricValues( 6 * (EigenFabricPerm(NodeIndexes(i))-1) + 1)=ai(1) 675 EigenFabricValues( 6 * (EigenFabricPerm(NodeIndexes(i))-1) + 2)=ai(2) 676 EigenFabricValues( 6 * (EigenFabricPerm(NodeIndexes(i))-1) + 3 )=ai(3) 677 EigenFabricValues( 6 * (EigenFabricPerm(NodeIndexes(i))-1) + 4 )=angle(1) 678 EigenFabricValues( 6 * (EigenFabricPerm(NodeIndexes(i))-1) + 5 )=angle(2) 679 EigenFabricValues( 6 * (EigenFabricPerm(NodeIndexes(i))-1) + 6 )=angle(3) 680 End do 681 END DO 682 683 END IF 684!------------------------------------------------------------------------------ 685 END DO ! of nonlinear iter 686!------------------------------------------------------------------------------ 687CONTAINS 688 689 SUBROUTINE GetMaterialDefs() 690 691 viscosityFile = ListGetString( Material ,'Viscosity File',GotIt, UnFoundFatal) 692 OPEN( 1, File = viscosityFile) 693 DO i=1,813 694 READ( 1, '(6(e14.8))' ) FabricGrid( 6*(i-1)+1:6*(i-1)+6 ) 695 END DO 696 READ(1 , '(e14.8)' ) FabricGrid(4879) 697 CLOSE(1) 698 699 rho = ListGetConstReal(Material, 'Interaction Parameter', GotIt ) 700 IF (.NOT.GotIt) THEN 701 WRITE(Message,'(A)') 'Interaction Parameter notfound. & 702 &Setting to the value in ViscosityFile' 703 CALL INFO('AIFlowSolve', Message, Level = 20) 704 rho = FabricGrid(4879) 705 ELSE 706 WRITE(Message,'(A,F10.4)') 'Interaction Parameter = ', rho 707 CALL INFO('AIFlowSolve', Message, Level = 20) 708 END IF 709 710 lambda = ListGetConstReal( Material, 'Diffusion Parameter', GotIt,UnFoundFatal=UnFoundFatal) 711 !Previous default value: lambda = 0.0_dp 712 WRITE(Message,'(A,F10.4)') 'Diffusion Parameter = ', lambda 713 CALL INFO('AIFlowSolve', Message, Level = 20) 714 715 Wn(2) = ListGetConstReal( Material , 'Powerlaw Exponent', GotIt,UnFoundFatal=UnFoundFatal) 716 !Previous default value: Wn(2) = 1.0 717 WRITE(Message,'(A,F10.4)') 'Powerlaw Exponent = ', Wn(2) 718 CALL INFO('AIFlowSolve', Message, Level = 20) 719 720 Wn(3) = ListGetConstReal( Material, 'Activation Energy 1', GotIt,UnFoundFatal=UnFoundFatal) 721 !Previous default value: Wn(3) = 1.0 722 WRITE(Message,'(A,F10.4)') 'Activation Energy 1 = ', Wn(3) 723 CALL INFO('AIFlowSolve', Message, Level = 20) 724 725 Wn(4) = ListGetConstReal( Material, 'Activation Energy 2', GotIt,UnFoundFatal=UnFoundFatal) 726 !Previous default value: Wn(4) = 1.0 727 WRITE(Message,'(A,F10.4)') 'Activation Energy 2 = ', Wn(4) 728 CALL INFO('AIFlowSolve', Message, Level = 20) 729 730 Wn(5) = ListGetConstReal(Material, 'Reference Temperature', GotIt,UnFoundFatal=UnFoundFatal) 731 !Previous default value: Wn(5) = -10.0 732 WRITE(Message,'(A,F10.4)') 'Reference Temperature = ', Wn(5) 733 CALL INFO('AIFlowSolve', Message, Level = 20) 734 735 Wn(6) = ListGetConstReal( Material, 'Limit Temperature', GotIt,UnFoundFatal=UnFoundFatal) 736 !Previous default value: Wn(6) = -10.0 737 WRITE(Message,'(A,F10.4)') 'Limit Temperature = ', Wn(6) 738 CALL INFO('AIFlowSolve', Message, Level = 20) 739!------------------------------------------------------------------------------ 740 END SUBROUTINE GetMaterialDefs 741!------------------------------------------------------------------------------ 742 743 744!------------------------------------------------------------------------------ 745 SUBROUTINE LocalMatrix( Comp, MASS, STIFF, FORCE, LOAD, & 746 NodalK1, NodalK2, NodalEuler1, NodalEuler2, NodalEuler3, & 747 NodalTemperature, NodalFluidity, NodalVelo, NodMeshVel, & 748 Element, n, Nodes, Wn, rho,lambda ) 749!------------------------------------------------------------------------------ 750 751 REAL(KIND=dp) :: STIFF(:,:),MASS(:,:) 752 REAL(KIND=dp) :: LOAD(:,:), NodalVelo(:,:),NodMeshVel(:,:) 753 REAL(KIND=dp), DIMENSION(:) :: FORCE, NodalK1, NodalK2, NodalEuler1, & 754 NodalEuler2, NodalEuler3, NodalTemperature, NodalFluidity 755 756 TYPE(Nodes_t) :: Nodes 757 TYPE(Element_t) :: Element 758 INTEGER :: n, Comp 759!------------------------------------------------------------------------------ 760! 761 REAL(KIND=dp) :: Basis(2*n),ddBasisddx(1,1,1) 762 REAL(KIND=dp) :: dBasisdx(2*n,3),SqrtElementMetric 763 764 REAL(KIND=dp) :: A1, A2, A3, E1, E2, E3, Theta 765 766 REAL(KIND=dp) :: A,M, hK,tau,pe1,pe2,unorm,C0, SU(n), SW(n) 767 REAL(KIND=dp) :: LoadAtIp, Temperature 768 REAL(KIND=dp) :: rho,lambda,Deq, ai(6),a4(9),hmax 769 770 INTEGER :: i,j,k,p,q,t,dim,NBasis,ind(3), DOFs = 1 771 772 REAL(KIND=dp) :: s,u,v,w, Radius, B(6,3), G(3,6), C44,C55,C66 773 REAL(KIND=dp) :: Wn(:),Velo(3),DStress(6),StrainR(6),Spin(3),SD(6) 774 775 REAL(KIND=dp) :: LGrad(3,3),StrainRate(3,3),D(6),angle(3),epsi 776 REAL(KIND=dp) :: ap(3),C(6,6),Spin1(3,3),Stress(3,3) 777 Integer :: INDi(6),INDj(6) 778 LOGICAL :: CSymmetry 779 780 LOGICAL :: Fond 781 782! 783 INTEGER :: N_Integ 784 REAL(KIND=dp), DIMENSION(:), POINTER :: U_Integ,V_Integ,W_Integ,S_Integ 785 786 LOGICAL :: stat 787 TYPE(GaussIntegrationPoints_t), TARGET :: IntegStuff 788 789 INTERFACE 790 FUNCTION BGlenT( Tc, W) 791 USE Types 792 REAL(KIND=dp) :: BGlenT,Tc,W(7) 793 END FUNCTION 794 795 SUBROUTINE IBOF(ai,a4) 796 USE Types 797 REAL(KIND=dp),intent(in) :: ai(6) 798 REAL(KIND=dp),intent(out) :: a4(9) 799 END SUBROUTINE 800 801 Subroutine R2Ro(ai,dim,a2,angle) 802 USE Types 803 REAL(KIND=dp),intent(in) :: ai(6) 804 Integer :: dim 805 REAL(KIND=dp),intent(out) :: a2(3), Angle(3) 806 End Subroutine R2Ro 807 808 Subroutine OPILGGE_ai_nl(a2,Angle,etaI,eta36) 809 USE Types 810 REAL(kind=dp), INTENT(in), DIMENSION(3) :: a2 811 REAL(kind=dp), INTENT(in), DIMENSION(3) :: Angle 812 REAL(kind=dp), INTENT(in), DIMENSION(:) :: etaI 813 REAL(kind=dp), INTENT(out), DIMENSION(6,6) :: eta36 814 END SUBROUTINE OPILGGE_ai_nl 815 816 END INTERFACE 817!------------------------------------------------------------------------------ 818 Fond=.False. 819 820 hmax = maxval (Nodes % y(1:n)) 821 822 dim = CoordinateSystemDimension() 823 824 FORCE = 0.0D0 825 MASS = 0.0D0 826 STIFF = 0.0D0 827! 828! Integration stuff: 829! ------------------ 830 NBasis = n 831 IntegStuff = GaussPoints( Element ) 832 833 U_Integ => IntegStuff % u 834 V_Integ => IntegStuff % v 835 W_Integ => IntegStuff % w 836 S_Integ => IntegStuff % s 837 N_Integ = IntegStuff % n 838 839 hk = ElementDiameter( Element, Nodes ) 840! 841! Now we start integrating: 842! ------------------------- 843 DO t=1,N_Integ 844 845 u = U_Integ(t) 846 v = V_Integ(t) 847 w = W_Integ(t) 848 849!------------------------------------------------------------------------------ 850! Basis function values & derivatives at the integration point 851!------------------------------------------------------------------------------ 852 stat = ElementInfo( Element,Nodes,u,v,w,SqrtElementMetric, & 853 Basis, dBasisdx, ddBasisddx, .FALSE. ) 854 855 s = SqrtElementMetric * S_Integ(t) 856!------------------------------------------------------------------------------ 857! 858! Orientation parameters at the integration point: 859! ------------------------------------------------ 860 A1 = SUM( NodalK1(1:n) * Basis(1:n) ) 861 A2 = SUM( NodalK2(1:n) * Basis(1:n) ) 862 A3 = 1._dp - A1 - A2 863 864 E1 = SUM( NodalEuler1(1:n) * Basis(1:n) ) 865 E2 = SUM( NodalEuler2(1:n) * Basis(1:n) ) 866 E3 = SUM( NodalEuler3(1:n) * Basis(1:n) ) 867! 868! Fluidity at the integration point: 869!--------------------------------------------- 870! Wn(1)=SUM( NodalFluidity(1:n)*Basis(1:n) ) 871! 872! Temperature at the integration point: 873! ------------------------------------- 874! Temperature = SUM( NodalTemperature(1:n)*Basis(1:n) ) 875! 876! Get theta parameter: (the (fluidity in the basal plane)/2, 877! function of the Temperature ) 878! ----------------------------------------------------- 879 Theta = 1._dp / ( FabricGrid(5) + FabricGrid(6) ) 880 !Theta = Theta 881 882! Strain-Rate, Stresses and Spin 883 884 CSymmetry = CurrentCoordinateSystem() == AxisSymmetric 885 886 Stress = 0.0 887 StrainRate = 0.0 888 Spin1 = 0.0 889! 890! Material parameters at that point 891! --------------------------------- 892! 893 ai(1)=A1 894 ai(2)=A2 895 ai(3)=A3 896 ai(4)=E1 897 ai(5)=E2 898 ai(6)=E3 899 900! fourth order orientation tensor 901 call IBOF(ai,a4) 902 903! A2 expressed in the orthotropic frame 904! 905 call R2Ro(ai,dim,ap,angle) 906 907! Get viscosity 908 909 CALL OPILGGE_ai_nl(ap, Angle, FabricGrid, C) 910 911! Compute strainRate and Spin : 912! ----------------------------- 913 914 LGrad = MATMUL( NodalVelo(1:3,1:n), dBasisdx(1:n,1:3) ) 915 916 StrainRate = 0.5 * ( LGrad + TRANSPOSE(LGrad) ) 917 918 Spin1 = 0.5 * ( LGrad - TRANSPOSE(LGrad) ) 919 920 IF ( CSymmetry ) THEN 921 StrainRate(1,3) = 0.0 922 StrainRate(2,3) = 0.0 923 StrainRate(3,1) = 0.0 924 StrainRate(3,2) = 0.0 925 StrainRate(3,3) = 0.0 926 927 Radius = SUM( Nodes % x(1:n) * Basis(1:n) ) 928 929 IF ( Radius > 10*AEPS ) THEN 930 StrainRate(3,3) = SUM( Nodalvelo(1,1:n) * Basis(1:n) ) / Radius 931 END IF 932 933 epsi = StrainRate(1,1)+StrainRate(2,2)+StrainRate(3,3) 934 DO i=1,3 935 StrainRate(i,i) = StrainRate(i,i) - epsi/3.0 936 END DO 937 938 ELSE 939 epsi = StrainRate(1,1)+StrainRate(2,2)+StrainRate(3,3) 940 DO i=1,dim 941 StrainRate(i,i) = StrainRate(i,i) - epsi/dim 942 END DO 943 944 END IF 945 946! 947! Compute deviatoric stresses: 948! ---------------------------- 949 D(1) = StrainRate(1,1) 950 D(2) = StrainRate(2,2) 951 D(3) = StrainRate(3,3) 952 D(4) = 2. * StrainRate(1,2) 953 D(5) = 2. * StrainRate(2,3) 954 D(6) = 2. * StrainRate(3,1) 955 956 INDi(1:6) = (/ 1, 2, 3, 1, 2, 3 /) 957 INDj(1:6) = (/ 1, 2, 3, 2, 3, 1 /) 958 DO k = 1, 2*dim 959 DO j = 1, 2*dim 960 Stress( INDi(k),INDj(k) ) = & 961 Stress( INDi(k),INDj(k) ) + C(k,j) * D(j) 962 END DO 963 IF (k > 3) Stress( INDj(k),INDi(k) ) = Stress( INDi(k),INDj(k) ) 964 END DO 965 966 967! SD=(1-r)D + r psi/2 S : 968! ----------------------- 969 SD=0._dp 970 DO i=1,2*dim 971 SD(i)= (1._dp - rho)*StrainRate(INDi(i),INDj(i)) + rho *& 972 Theta * Stress(INDi(i),INDj(i)) 973 END DO 974 Do i=1,2*dim-3 975 Spin(i)=Spin1(INDi(i+3),INDj(i+3)) 976 End do 977 978 Deq=sqrt(2._dp*(SD(1)*SD(1)+SD(2)*SD(2)+SD(3)*SD(3)+2._dp* & 979 (SD(4)*SD(4)+SD(5)*SD(5)+SD(6)*SD(6)))/3._dp) 980! 981! Velocity : 982! ---------- 983 Velo = 0.0d0 984 DO i=1,dim 985 Velo(i) = SUM( Basis(1:n) * (NodalVelo(i,1:n) - NodMeshVel(i,1:n)) ) 986 END DO 987 Unorm = SQRT( SUM( Velo**2._dp ) ) 988 989! 990! Reaction coefficient: 991! --------------------- 992 SELECT CASE(comp) 993 CASE(1) 994 !C0 = -2._dp*(SD(1)-SD(3)) 995 C0=-2._dp*SD(1)-3._dp*lambda*Deq 996 CASE(2) 997 !C0 = -2._dp*(SD(2)-SD(3)) 998 C0 = -2._dp*SD(2)-3._dp*lambda*Deq 999 1000 CASE(3) 1001 !C0 = -(SD(1)+SD(2)-2._dp*SD(3)) 1002 C0 = -(SD(1)+SD(2))-3._dp*lambda*Deq 1003 1004 CASE(4) 1005 C0 = -(SD(2)-SD(3))-3._dp*lambda*Deq 1006 1007 CASE(5) 1008 C0 = -(SD(1)-SD(3))-3._dp*lambda*Deq 1009 1010 1011 END SELECT 1012 1013 If (Fond) C0=0._dp 1014 1015! Loop over basis functions (of both unknowns and weights): 1016! --------------------------------------------------------- 1017 DO p=1,NBasis 1018 DO q=1,NBasis 1019 A = 0.0d0 1020 M = Basis(p) * Basis(q) 1021! 1022! Reaction terms: 1023! --------------- 1024 A = A - C0 * Basis(q) * Basis(p) 1025 1026 ! 1027 ! Advection terms: 1028 ! ---------------- 1029 DO j=1,dim 1030 A = A - Velo(j) * Basis(q) * dBasisdx(p,j) 1031 END DO 1032 1033! Add nodal matrix to element matrix: 1034! ----------------------------------- 1035 MASS( p,q ) = MASS( p,q ) + s * M 1036 STIFF( p,q ) = STIFF( p,q ) + s * A 1037 END DO 1038 1039 1040! 1041! The righthand side...: 1042! ---------------------- 1043 SELECT CASE(comp) 1044 1045 CASE(1) 1046 1047 LoadAtIp = 2._dp*( Spin(1)*E1 + SD(4)*(2._dp*a4(7)-E1) + & 1048 SD(1)*a4(1) + SD(2)*a4(3) - SD(3)*(-A1+a4(1)+a4(3)) ) + & 1049 lambda*Deq 1050 1051 IF(dim == 3) THEN 1052 LoadAtIp = LoadAtIp + 2._dp*( -Spin(3)*E3 + & 1053 SD(6)*(2._dp*a4(6)-E3) + 2._dp*SD(5)*a4(4) ) 1054 END IF 1055 1056 1057 CASE(2) 1058 LoadAtIp = 2._dp*( -Spin(1)*E1 + SD(4)*(2._dp*a4(9)-E1) + & 1059 SD(2)*a4(2) + SD(1)*a4(3) - SD(3)*(-A2+a4(2)+a4(3)) ) + & 1060 lambda*Deq 1061 1062 IF(dim == 3) THEN 1063 LoadAtIp = LoadAtIp + 2._dp*( Spin(2)*E2 + & 1064 SD(5)*(2._dp*a4(8)-E2) + 2._dp*SD(6)*a4(5) ) 1065 END IF 1066 1067 1068 CASE(3) 1069 LoadAtIp = Spin(1)*(A2-A1) + SD(4)*(4._dp*a4(3)-A1-A2) + & 1070 2._dp* ( SD(1)*a4(7) + SD(2)*a4(9) - SD(3)*(-E1+a4(7)+a4(9)) ) 1071 1072 IF(dim == 3) THEN 1073 LoadAtIp = LoadAtIp - Spin(3)*E2 + Spin(2)*E3 & 1074 + SD(6)*(4._dp*a4(4)-E2) + SD(5)*(4._dp*a4(5)-E3) 1075 END IF 1076 1077 CASE(4) 1078 LoadAtIp = Spin(2)*(A3-A2) + Spin(3)*E1 - Spin(1)*E3 +& 1079 SD(4)*(4._dp*a4(5)-E3) + SD(5)*(3._dp*A2-A3-4._dp*(a4(2)+a4(3))) & 1080 + SD(6)*(3._dp*E1-4._dp*(a4(7)+a4(9))) + & 1081 2._dp*( SD(1)*a4(4) + SD(2)*a4(8) - SD(3)*(a4(4)+a4(8)) ) 1082 1083 CASE(5) 1084 LoadAtIp = Spin(3)*(A1-A3) + Spin(1)*E2 - Spin(2)*E1 +& 1085 SD(4)*(4._dp*a4(4)-E2) + & 1086 SD(6)*(3._dp*A1-A3-4.*(a4(1)+a4(3))) + SD(5)*(3._dp*E1-4._dp*(a4(7)+a4(9))) + & 1087 2._dp*( SD(1)*a4(6) + SD(2)*a4(5) - SD(3)*(a4(6)+a4(5)) ) 1088 1089 1090 1091 END SELECT 1092 1093 If (Fond) LoadAtIp=0._dp 1094 LoadAtIp= LoadAtIp * Basis(p) 1095 FORCE(p) = FORCE(p) + s*LoadAtIp 1096 END DO 1097 1098 1099 END DO 1100 1101 1000 FORMAT((a),x,i2,x,6(e13.5,x)) 1102 1001 FORMAT(6(e13.5,x)) 1103!------------------------------------------------------------------------------ 1104 END SUBROUTINE LocalMatrix 1105!------------------------------------------------------------------------------ 1106 1107!------------------------------------------------------------------------------ 1108 SUBROUTINE FindParentUVW( Edge, nEdge, Parent, nParent, U, V, W, Basis ) 1109!------------------------------------------------------------------------------ 1110 IMPLICIT NONE 1111 TYPE(Element_t), POINTER :: Edge, Parent 1112 INTEGER :: nEdge, nParent 1113 REAL( KIND=dp ) :: U, V, W, Basis(:) 1114!------------------------------------------------------------------------------ 1115 INTEGER :: i, j,l 1116 REAL(KIND=dp) :: NodalParentU(nEdge),NodalParentV(nEdge),NodalParentW(nEdge) 1117!------------------------------------------------------------------------------ 1118 DO i = 1,nEdge 1119 DO j = 1,nParent 1120 IF ( Edge % NodeIndexes(i) == Parent % NodeIndexes(j) ) THEN 1121 NodalParentU(i) = Parent % Type % NodeU(j) 1122 NodalParentV(i) = Parent % Type % NodeV(j) 1123 NodalParentW(i) = Parent % Type % NodeW(j) 1124 EXIT 1125 END IF 1126 END DO 1127 END DO 1128 U = SUM( Basis(1:nEdge) * NodalParentU(1:nEdge) ) 1129 V = SUM( Basis(1:nEdge) * NodalParentV(1:nEdge) ) 1130 W = SUM( Basis(1:nEdge) * NodalParentW(1:nEdge) ) 1131!------------------------------------------------------------------------------ 1132 END SUBROUTINE FindParentUVW 1133!------------------------------------------------------------------------------ 1134 1135 1136!------------------------------------------------------------------------------ 1137 SUBROUTINE LocalJumps( STIFF,Edge,n,LeftParent,n1,RightParent,n2,Velo,MeshVelo ) 1138!------------------------------------------------------------------------------ 1139 REAL(KIND=dp) :: STIFF(:,:), Velo(:,:),MeshVelo(:,:) 1140 INTEGER :: n,n1,n2 1141 TYPE(Element_t), POINTER :: Edge, LeftParent, RightParent 1142!------------------------------------------------------------------------------ 1143 REAL(KIND=dp) :: EdgeBasis(n), EdgedBasisdx(n,3), EdgeddBasisddx(n,3,3) 1144 REAL(KIND=dp) :: LeftBasis(n1), LeftdBasisdx(n1,3), LeftddBasisddx(n1,3,3) 1145 REAL(KIND=dp) :: RightBasis(n2), RightdBasisdx(n2,3), RightddBasisddx(n2,3,3) 1146 REAL(KIND=dp) :: Jump(n1+n2), Average(n1+n2) 1147 REAL(KIND=dp) :: detJ, U, V, W, S, Udotn, xx, yy 1148 LOGICAL :: Stat 1149 INTEGER :: i, j, p, q, dim, t, nEdge, nParent 1150 TYPE(GaussIntegrationPoints_t) :: IntegStuff 1151 REAL(KIND=dp) :: hE, Normal(3), cu(3), LeftOut(3) 1152 1153 TYPE(Nodes_t) :: EdgeNodes, LeftParentNodes, RightParentNodes 1154 1155 Save EdgeNodes, LeftParentNodes, RightParentNodes 1156!------------------------------------------------------------------------------ 1157 dim = CoordinateSystemDimension() 1158 STIFF = 0.0d0 1159 1160 CALL GetElementNodes( EdgeNodes, Edge ) 1161 CALL GetElementNodes( LeftParentNodes, LeftParent ) 1162 CALL GetElementNodes( RightParentNodes, RightParent ) 1163!------------------------------------------------------------------------------ 1164! Numerical integration over the edge 1165!------------------------------------------------------------------------------ 1166 IntegStuff = GaussPoints( Edge ) 1167 1168 LeftOut(1) = SUM( LeftParentNodes % x(1:n1) ) / n1 1169 LeftOut(2) = SUM( LeftParentNodes % y(1:n1) ) / n1 1170 LeftOut(3) = SUM( LeftParentNodes % z(1:n1) ) / n1 1171 LeftOut(1) = SUM( EdgeNodes % x(1:n) ) / n - LeftOut(1) 1172 LeftOut(2) = SUM( EdgeNodes % y(1:n) ) / n - LeftOut(2) 1173 LeftOut(3) = SUM( EdgeNodes % z(1:n) ) / n - LeftOut(3) 1174 1175 DO t=1,IntegStuff % n 1176 U = IntegStuff % u(t) 1177 V = IntegStuff % v(t) 1178 W = IntegStuff % w(t) 1179 S = IntegStuff % s(t) 1180 1181 ! Basis function values & derivatives at the integration point: 1182 !-------------------------------------------------------------- 1183 stat = ElementInfo( Edge, EdgeNodes, U, V, W, detJ, & 1184 EdgeBasis, EdgedBasisdx, EdgeddBasisddx, .FALSE. ) 1185 1186 S = S * detJ 1187 1188 Normal = NormalVector( Edge, EdgeNodes, U, V, .FALSE. ) 1189 IF ( SUM( LeftOut*Normal ) < 0 ) Normal = -Normal 1190 1191 ! Find basis functions for the parent elements: 1192 ! --------------------------------------------- 1193 CALL FindParentUVW( Edge,n,LeftParent,n1,U,V,W,EdgeBasis ) 1194 stat = ElementInfo( LeftParent, LeftParentNodes, U, V, W, detJ, & 1195 LeftBasis, LeftdBasisdx, LeftddBasisddx, .FALSE. ) 1196 1197 CALL FindParentUVW( Edge,n,RightParent,n2,U,V,W,EdgeBasis ) 1198 stat = ElementInfo( RightParent, RightParentNodes, U, V, W, detJ, & 1199 RightBasis, RightdBasisdx, RightddBasisddx, .FALSE. ) 1200 1201 ! Integrate jump terms: 1202 ! --------------------- 1203 Jump(1:n1) = LeftBasis(1:n1) 1204 Jump(n1+1:n1+n2) = -RightBasis(1:n2) 1205 1206 Average(1:n1) = LeftBasis(1:n1) / 2 1207 Average(n1+1:n1+n2) = RightBasis(1:n2) / 2 1208 1209 cu = 0.0d0 1210 DO i=1,dim 1211 cu(i) = SUM( (Velo(i,1:n)-MeshVelo(i,1:n)) * EdgeBasis(1:n) ) 1212 END DO 1213 Udotn = SUM( Normal * cu ) 1214 1215 DO p=1,n1+n2 1216 DO q=1,n1+n2 1217 STIFF(p,q) = STIFF(p,q) + s * Udotn * Average(q) * Jump(p) 1218 STIFF(p,q) = STIFF(p,q) + s * ABS(Udotn)/2 * Jump(q) * Jump(p) 1219 END DO 1220 END DO 1221 END DO 1222!------------------------------------------------------------------------------ 1223 END SUBROUTINE LocalJumps 1224!------------------------------------------------------------------------------ 1225 1226 1227 1228!------------------------------------------------------------------------------ 1229 SUBROUTINE LocalMatrixBoundary( STIFF, FORCE, LOAD, & 1230 Element, n, ParentElement, np, Velo,MeshVelo, InFlowBC ) 1231!------------------------------------------------------------------------------ 1232 REAL(KIND=dp) :: STIFF(:,:), FORCE(:), LOAD(:), Velo(:,:),MeshVelo(:,:) 1233 INTEGER :: n, np 1234 LOGICAL :: InFlowBC 1235 TYPE(Element_t), POINTER :: Element, ParentElement 1236!------------------------------------------------------------------------------ 1237 REAL(KIND=dp) :: Basis(n), dBasisdx(n,3), ddBasisddx(n,3,3) 1238 REAL(KIND=dp) :: ParentBasis(np), ParentdBasisdx(np,3), ParentddBasisddx(np,3,3) 1239 INTEGER :: i,j,p,q,t,dim 1240 1241 REAL(KIND=dp) :: Normal(3), g, L, Udotn, UdotnA, cu(3), detJ,U,V,W,S 1242 LOGICAL :: Stat 1243 TYPE(GaussIntegrationPoints_t) :: IntegStuff 1244 1245 TYPE(Nodes_t) :: Nodes, ParentNodes 1246 SAVE Nodes, ParentNodes 1247!------------------------------------------------------------------------------ 1248 dim = CoordinateSystemDimension() 1249 FORCE = 0.0d0 1250 STIFF = 0.0d0 1251 1252 CALL GetElementNodes( Nodes, Element ) 1253 CALL GetElementNodes( ParentNodes, ParentElement ) 1254 1255 ! Numerical integration: 1256 !----------------------- 1257 IntegStuff = GaussPoints( Element ) 1258! 1259! Compute the average velocity.dot.Normal 1260! 1261 UdotnA = 0.0 1262 DO t=1,IntegStuff % n 1263 U = IntegStuff % u(t) 1264 V = IntegStuff % v(t) 1265 W = IntegStuff % w(t) 1266 S = IntegStuff % s(t) 1267 1268 Normal = NormalVector( Element, Nodes, U, V, .TRUE. ) 1269 1270 ! Basis function values & derivatives at the integration point: 1271 ! ------------------------------------------------------------- 1272 stat = ElementInfo( Element, Nodes, U, V, W, detJ, & 1273 Basis, dBasisdx, ddBasisddx, .FALSE. ) 1274 S = S * detJ 1275 cu = 0.0d0 1276 DO i=1,dim 1277 cu(i) = SUM( (Velo(i,1:n)-MeshVelo(i,1:n)) * Basis(1:n) ) 1278 END DO 1279 UdotnA = UdotnA + s*SUM( Normal * cu ) 1280 1281 END DO 1282 1283 DO t=1,IntegStuff % n 1284 U = IntegStuff % u(t) 1285 V = IntegStuff % v(t) 1286 W = IntegStuff % w(t) 1287 S = IntegStuff % s(t) 1288 1289 Normal = NormalVector( Element, Nodes, U, V, .TRUE. ) 1290 1291 ! Basis function values & derivatives at the integration point: 1292 ! ------------------------------------------------------------- 1293 stat = ElementInfo( Element, Nodes, U, V, W, detJ, & 1294 Basis, dBasisdx, ddBasisddx, .FALSE. ) 1295 S = S * detJ 1296 1297 CALL FindParentUVW( Element, n, ParentElement, np, U, V, W, Basis ) 1298 stat = ElementInfo( ParentElement, ParentNodes, U, V, W, & 1299 detJ, ParentBasis, ParentdBasisdx, ParentddBasisddx, .FALSE. ) 1300 1301 L = SUM( LOAD(1:n) * Basis(1:n) ) 1302 cu = 0.0d0 1303 DO i=1,dim 1304 cu(i) = SUM( (Velo(i,1:n)-MeshVelo(i,1:n)) * Basis(1:n) ) 1305 END DO 1306 Udotn = SUM( Normal * cu ) 1307 1308 DO p = 1,np 1309 IF (InFlowBC .And. (UdotnA < 0.) ) THEN 1310 FORCE(p) = FORCE(p) - s * Udotn*L*ParentBasis(p) 1311 ELSE 1312 DO q=1,np 1313 STIFF(p,q) = STIFF(p,q) + s*Udotn*ParentBasis(q)*ParentBasis(p) 1314 END DO 1315 END IF 1316 END DO 1317 END DO 1318!------------------------------------------------------------------------------ 1319 END SUBROUTINE LocalMatrixBoundary 1320!------------------------------------------------------------------------------ 1321 1322 1323 1324!------------------------------------------------------------------------------ 1325 END SUBROUTINE FabricSolver 1326!------------------------------------------------------------------------------ 1327 1328