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! * Generic utilities related to components 27! * 28! ****************************************************************************** 29! * 30! * Authors: Juha Ruokolainen 31! * Email: Juha.Ruokolainen@csc.fi 32! * Web: http://www.csc.fi/elmer 33! * Address: CSC - IT Center for Science Ltd. 34! * Keilaranta 14 35! * 02101 Espoo, Finland 36! * 37! * Original Date: 28 Sep 1998 38! * 39! *****************************************************************************/ 40 41!> Generic utilities related to components 42!------------------------------------------------------------------------------ 43 44!> \ingroup ElmerLib 45!> \{ 46 47 48MODULE ComponentUtils 49 50 USE ElementUtils 51 USE ModelDescription 52 IMPLICIT NONE 53 54 CONTAINS 55 56 57!------------------------------------------------------------------------------ 58!> Compute reduction operators for a given component with given nodal vector field. 59!> Force is simple sum of nodal forces 60!> Moment is moment of nodal forces about a given center point 61!> Torque is moment of nodal forces about a given rotational axis 62!> If the given field is a elemental (DG) field it may be reduced by 63!> optional SetPerm reordering for minimal discontinuous set. 64!------------------------------------------------------------------------------ 65 SUBROUTINE ComponentNodalForceReduction(Model, Mesh, CompParams, NF, & 66 Force, Moment, Torque, SetPerm ) 67!------------------------------------------------------------------------------ 68 TYPE(Model_t) :: Model 69 TYPE(Mesh_t), POINTER :: Mesh 70 TYPE(ValueList_t), POINTER :: CompParams 71 TYPE(Variable_t), POINTER :: NF 72 REAL(KIND=dp), OPTIONAL :: Moment(3), Force(3), Torque 73 INTEGER, POINTER, OPTIONAL :: SetPerm(:) 74!------------------------------------------------------------------------------ 75! Local variables 76!------------------------------------------------------------------------------ 77 TYPE(Element_t), POINTER :: Element 78 LOGICAL, ALLOCATABLE :: VisitedNode(:) 79 REAL(KIND=dp) :: Origin(3), Axis(3), P(3), F(3), v1(3), v2(3) 80 REAL(KIND=dp), POINTER :: Pwrk(:,:) 81 INTEGER :: t, i, j, k, n, dofs, globalnode 82 LOGICAL :: ElementalVar, Found, NeedLocation 83 INTEGER, POINTER :: MasterEntities(:),NodeIndexes(:),DofIndexes(:) 84 LOGICAL :: VisitNodeOnlyOnce 85 INTEGER :: FirstElem, LastElem 86 LOGICAL :: BcMode 87 88 89 CALL Info('ComponentNodalForceReduction','Performing reduction for component: '& 90 //TRIM(ListGetString(CompParams,'Name')),Level=10) 91 92 IF(.NOT. (PRESENT(Torque) .OR. PRESENT(Moment) .OR. PRESENT(Force) ) ) THEN 93 CALL Warn('ComponentNodalForceReduction','Nothing to compute!') 94 RETURN 95 END IF 96 97 IF( PRESENT(Torque)) Torque = 0.0_dp 98 IF( PRESENT(Moment)) Moment = 0.0_dp 99 IF( PRESENT(Force)) Force = 0.0_dp 100 101 BcMode = .FALSE. 102 MasterEntities => ListGetIntegerArray( CompParams,'Master Bodies',Found ) 103 IF( .NOT. Found ) THEN 104 MasterEntities => ListGetIntegerArray( CompParams,'Master Boundaries',Found ) 105 BcMode = .TRUE. 106 END IF 107 108 IF(.NOT. Found ) THEN 109 CALL Warn('ComponentNodalForceReduction',& 110 '> Master Bodies < or > Master Boundaries < not given') 111 RETURN 112 END IF 113 114 NeedLocation = PRESENT( Moment ) .OR. PRESENT( Torque ) 115 116 ! User may specific origin and axis for torque computation 117 ! By default (0,0,0) is the origin, and (0,0,1) the axis. 118 Pwrk => ListGetConstRealArray( CompParams,'Torque Origin',Found ) 119 IF( Found ) THEN 120 IF( SIZE(Pwrk,1) /= 3 .OR. SIZE(Pwrk,2) /= 1 ) THEN 121 CALL Fatal('ComponentNodalForceReduction','Size of > Torque Origin < should be 3!') 122 END IF 123 Origin = Pwrk(1:3,1) 124 ELSE 125 Origin = 0.0_dp 126 END IF 127 Pwrk => ListGetConstRealArray( CompParams,'Torque Axis',Found ) 128 IF( Found ) THEN 129 IF( SIZE(Pwrk,1) /= 3 .OR. SIZE(Pwrk,2) /= 1 ) THEN 130 CALL Fatal('ComponentNodalForceReduction','Size of > Torque Axis < should be 3!') 131 END IF 132 Axis = Pwrk(1:3,1) 133 ! Normalize axis is it should just be used for the direction 134 Axis = Axis / SQRT( SUM( Axis*Axis ) ) 135 ELSE 136 Axis = 0.0_dp 137 Axis(3) = 1.0_dp 138 END IF 139 140 ElementalVar = ( NF % TYPE == Variable_on_nodes_on_elements ) 141 IF( PRESENT( SetPerm ) .AND. .NOT. ElementalVar ) THEN 142 CALL Fatal('ComponentNodalForceReduction','SetPerm is usable only for elemental fields') 143 END IF 144 145 dofs = NF % Dofs 146 IF( dofs == 2 ) F(3) = 0.0_dp 147 148 ! For nodal field compute only once each node 149 ! For DG field each node is visited only once by construction 150 VisitNodeOnlyOnce = .NOT. ElementalVar .OR. PRESENT(SetPerm) 151 IF( VisitNodeOnlyOnce ) THEN 152 IF( PRESENT( SetPerm ) ) THEN 153 n = MAXVAL( SetPerm ) 154 ELSE 155 n = Mesh % NumberOfNodes 156 END IF 157 ALLOCATE(VisitedNode( n ) ) 158 VisitedNode = .FALSE. 159 END IF 160 161 IF( BcMode ) THEN 162 FirstElem = Mesh % NumberOfBulkElements + 1 163 LastElem = Mesh % NumberOfBulkElements + Mesh % NumberOfBoundaryElements 164 ELSE 165 FirstElem = 1 166 LastElem = Mesh % NumberOfBulkElements 167 END IF 168 169 170 DO t=FirstElem,LastElem 171 Element => Mesh % Elements(t) 172 173 IF( BcMode ) THEN 174 IF( ALL( MasterEntities /= Element % BoundaryInfo % Constraint ) ) CYCLE 175 ELSE 176 IF( ALL( MasterEntities /= Element % BodyId ) ) CYCLE 177 END IF 178 179 n = Element % TYPE % NumberOfNodes 180 NodeIndexes => Element % NodeIndexes 181 IF( ElementalVar ) THEN 182 DofIndexes => Element % DGIndexes 183 ELSE 184 DofIndexes => NodeIndexes 185 END IF 186 187 DO i=1,n 188 j = DofIndexes(i) 189 k = NF % Perm(j) 190 IF( k == 0 ) CYCLE 191 192 IF( VisitNodeOnlyOnce ) THEN 193 IF( PRESENT( SetPerm ) ) j = SetPerm(j) 194 IF( VisitedNode(j) ) CYCLE 195 VisitedNode(j) = .TRUE. 196 END IF 197 198 globalnode = NodeIndexes(i) 199 200 ! Only compute the parallel reduction once 201 IF( ParEnv % PEs > 1 ) THEN 202 IF( Mesh % ParallelInfo % NeighbourList(globalnode) % Neighbours(1) /= ParEnv % MyPE ) CYCLE 203 END IF 204 205 F(1) = NF % Values( dofs*(k-1) + 1) 206 F(2) = NF % Values( dofs*(k-1) + 2) 207 IF( dofs == 3 ) THEN 208 F(3) = NF % Values( dofs*(k-1) + 3) 209 END IF 210 211 IF( PRESENT( Force ) ) THEN 212 ! calculate simple sum 213 Force = Force + F 214 END IF 215 216 IF( NeedLocation ) THEN 217 P(1) = Mesh % Nodes % x(globalnode) 218 P(2) = Mesh % Nodes % y(globalnode) 219 P(3) = Mesh % Nodes % z(globalnode) 220 221 v1 = P - Origin 222 223 ! Calculate moment 224 IF( PRESENT( Moment ) ) THEN 225 Moment = Moment + CrossProduct(v1,F) 226 END IF 227 228 ! Calculate torque around an axis 229 IF( PRESENT( Torque ) ) THEN 230 v1 = (1.0_dp - SUM(Axis*v1) ) * v1 231 v2 = CrossProduct(v1,F) 232 Torque = Torque + SUM(Axis*v2) 233 END IF 234 END IF 235 236 END DO 237 END DO 238 239 IF( PRESENT( Force ) ) THEN 240 DO i=1,3 241 Force(i) = ParallelReduction(Force(i)) 242 END DO 243 END IF 244 245 IF( PRESENT( Moment ) ) THEN 246 DO i=1,3 247 Moment(i) = ParallelReduction(Moment(i)) 248 END DO 249 END IF 250 251 IF( PRESENT( Torque ) ) THEN 252 Torque = ParallelReduction(Torque) 253 END IF 254 255!------------------------------------------------------------------------------ 256 END SUBROUTINE ComponentNodalForceReduction 257!------------------------------------------------------------------------------ 258 259 260 261 262!------------------------------------------------------------------------------ 263!> Perform reduction from distributed data to components. 264!> These are similar operations as the stastistical operations in SaveScalars. 265!------------------------------------------------------------------------------ 266 FUNCTION ComponentNodalReduction(Model, Mesh, CompParams, Var, OperName ) RESULT ( OperX ) 267!------------------------------------------------------------------------------ 268 TYPE(Model_t) :: Model 269 TYPE(Mesh_t), POINTER :: Mesh 270 TYPE(ValueList_t), POINTER :: CompParams 271 TYPE(Variable_t), POINTER :: Var 272 CHARACTER(LEN=MAX_NAME_LEN) :: OperName 273 REAL(KIND=dp) :: OperX 274!------------------------------------------------------------------------------ 275! Local variables 276!------------------------------------------------------------------------------ 277 TYPE(Element_t), POINTER :: Element 278 LOGICAL, ALLOCATABLE :: VisitedNode(:) 279 INTEGER :: t, i, j, k, l, n, NoDofs, globalnode, sumi 280 REAL(KIND=dp) :: X, Minimum, Maximum, AbsMinimum, AbsMaximum, SumX, SumXX, SumAbsX 281 LOGICAL :: ElementalVar, Found 282 INTEGER, POINTER :: MasterEntities(:),NodeIndexes(:),DofIndexes(:) 283 LOGICAL :: VisitNodeOnlyOnce, Initialized 284 INTEGER :: FirstElem, LastElem 285 LOGICAL :: BcMode 286 287 288 CALL Info('ComponentNodalReduction','Performing reduction for component: '& 289 //TRIM(ListGetString(CompParams,'Name')),Level=10) 290 291 OperX = 0.0_dp 292 293 BcMode = .FALSE. 294 MasterEntities => ListGetIntegerArray( CompParams,'Master Bodies',Found ) 295 IF( .NOT. Found ) THEN 296 MasterEntities => ListGetIntegerArray( CompParams,'Master Boundaries',Found ) 297 BcMode = .TRUE. 298 END IF 299 300 IF(.NOT. Found ) THEN 301 CALL Warn('ComponentNodalReduction',& 302 '> Master Bodies < or > Master Boundaries < not given') 303 RETURN 304 END IF 305 306 IF( BcMode ) THEN 307 FirstElem = Mesh % NumberOfBulkElements + 1 308 LastElem = Mesh % NumberOfBulkElements + Mesh % NumberOfBoundaryElements 309 ELSE 310 FirstElem = 1 311 LastElem = Mesh % NumberOfBulkElements 312 END IF 313 314 ElementalVar = ( Var % TYPE == Variable_on_nodes_on_elements ) 315 NoDofs = Var % Dofs 316 Initialized = .FALSE. 317 318 ! For nodal field compute only once each node 319 ! For DG field each node is visited only once by construction 320 VisitNodeOnlyOnce = .NOT. ElementalVar 321 IF( VisitNodeOnlyOnce ) THEN 322 n = Mesh % NumberOfNodes 323 ALLOCATE(VisitedNode( n ) ) 324 VisitedNode = .FALSE. 325 END IF 326 327 sumi = 0 328 sumx = 0.0_dp 329 sumxx = 0.0_dp 330 sumabsx = 0.0_dp 331 Maximum = 0.0_dp 332 Minimum = 0.0_dp 333 AbsMaximum = 0.0_dp 334 AbsMinimum = 0.0_dp 335 336 DO t=FirstElem,LastElem 337 Element => Mesh % Elements(t) 338 339 IF( BcMode ) THEN 340 IF( ALL( MasterEntities /= Element % BoundaryInfo % Constraint ) ) CYCLE 341 ELSE 342 IF( ALL( MasterEntities /= Element % BodyId ) ) CYCLE 343 END IF 344 345 n = Element % TYPE % NumberOfNodes 346 NodeIndexes => Element % NodeIndexes 347 IF( ElementalVar ) THEN 348 DofIndexes => Element % DGIndexes 349 ELSE 350 DofIndexes => NodeIndexes 351 END IF 352 353 DO i=1,n 354 j = DofIndexes(i) 355 IF( ASSOCIATED( Var % Perm ) ) THEN 356 k = Var % Perm(j) 357 IF( k == 0 ) CYCLE 358 ELSE 359 k = j 360 END IF 361 362 IF( VisitNodeOnlyOnce ) THEN 363 IF( VisitedNode(j) ) CYCLE 364 VisitedNode(j) = .TRUE. 365 END IF 366 367 globalnode = NodeIndexes(i) 368 369 ! Only compute the parallel reduction once 370 IF( ParEnv % PEs > 1 ) THEN 371 IF( Mesh % ParallelInfo % NeighbourList(globalnode) % Neighbours(1) /= ParEnv % MyPE ) CYCLE 372 END IF 373 374 IF(NoDofs <= 1) THEN 375 x = Var % Values(k) 376 ELSE 377 x = 0.0d0 378 DO l=1,NoDofs 379 x = x + Var % Values(NoDofs*(k-1)+l) ** 2 380 END DO 381 x = SQRT(x) 382 END IF 383 384 IF(.NOT. Initialized) THEN 385 Initialized = .TRUE. 386 Maximum = x 387 Minimum = x 388 AbsMaximum = x 389 AbsMinimum = x 390 END IF 391 392 sumi = sumi + 1 393 sumx = sumx + x 394 sumxx = sumxx + x*x 395 sumabsx = sumabsx + ABS( x ) 396 Maximum = MAX(x,Maximum) 397 Minimum = MIN(x,Minimum) 398 IF(ABS(x) > ABS(AbsMaximum) ) AbsMaximum = x 399 IF(ABS(x) < ABS(AbsMinimum) ) AbsMinimum = x 400 END DO 401 END DO 402 403 404 sumi = NINT( ParallelReduction(1.0_dp * sumi) ) 405 IF( sumi == 0 ) THEN 406 CALL Warn('ComponentNodalReduction','No active nodes to reduced!') 407 RETURN 408 END IF 409 410 411 SELECT CASE(OperName) 412 413 CASE ('sum') 414 sumx = ParallelReduction(sumx) 415 operx = sumx 416 417 CASE ('sum abs') 418 sumx = ParallelReduction(sumabsx) 419 operx = sumabsx 420 421 CASE ('min') 422 minimum = ParallelReduction(minimum,1) 423 operx = Minimum 424 425 CASE ('max') 426 maximum = ParallelReduction(maximum,2) 427 operx = Maximum 428 429 CASE ('min abs') 430 Absminimum = ParallelReduction(AbsMinimum,1) 431 operx = AbsMinimum 432 433 CASE ('max abs') 434 Absmaximum = ParallelReduction(AbsMaximum,2) 435 operx = AbsMaximum 436 437 CASE ('range') 438 minimum = ParallelReduction(minimum,1) 439 maximum = ParallelReduction(maximum,2) 440 operx = Maximum - Minimum 441 442 CASE ('mean') 443 sumx = ParallelReduction(sumx) 444 operx = sumx / sumi 445 446 CASE ('mean abs') 447 sumx = ParallelReduction(sumabsx) 448 operx = sumabsx / sumi 449 450 CASE ('variance') 451 sumx = ParallelReduction(sumx) 452 sumxx = ParallelReduction(sumxx) 453 Operx = SQRT( sumxx/sumi-(sumx*sumx)/(sumi*sumi) ) 454 455 CASE DEFAULT 456 CALL Warn('ComponentNodalReduction','Unknown statistical operator!') 457 458 END SELECT 459 460 CALL Info('ComponentNodalReduction','Reduction operator finished',Level=12) 461 462!------------------------------------------------------------------------------ 463 END FUNCTION ComponentNodalReduction 464!------------------------------------------------------------------------------ 465 466 467!------------------------------------------------------------------------------ 468!> Perform reduction from distributed data to components. 469!> These are similar operations as the stastistical operations in SaveScalars. 470!------------------------------------------------------------------------------ 471 FUNCTION ComponentIntegralReduction(Model, Mesh, CompParams, Var, & 472 OperName, CoeffName, GotCoeff ) RESULT ( OperX ) 473!------------------------------------------------------------------------------ 474 TYPE(Model_t) :: Model 475 TYPE(Mesh_t), POINTER :: Mesh 476 TYPE(ValueList_t), POINTER :: CompParams 477 TYPE(Variable_t), POINTER :: Var 478 CHARACTER(LEN=MAX_NAME_LEN) :: OperName, CoeffName 479 LOGICAL :: GotCoeff 480 REAL(KIND=dp) :: OperX 481!------------------------------------------------------------------------------ 482! Local variables 483!------------------------------------------------------------------------------ 484 TYPE(Element_t), POINTER :: Element 485 INTEGER :: t, i, j, k, n, NoDofs 486 INTEGER, POINTER :: NodeIndexes(:), DofIndexes(:) 487 REAL(KIND=dp) :: SqrtElementMetric,U,V,W,S,x,Grad(3) 488 REAL(KIND=dp) :: func, CoeffAtIp, integral, vol 489 LOGICAL :: ElementalVar, Found, Stat 490 INTEGER :: PermIndexes(Model % MaxElementNodes) 491 REAL(KIND=dp) :: Basis(Model % MaxElementNodes), dBasisdx(Model % MaxElementNodes,3) 492 REAL(KIND=dp) :: Coeff(Model % MaxElementNodes) 493 TYPE(GaussIntegrationPoints_t) :: IntegStuff 494 INTEGER, POINTER :: MasterEntities(:) 495 TYPE(Nodes_t), SAVE :: ElementNodes 496 LOGICAL, SAVE :: AllocationsDone = .FALSE. 497 INTEGER :: FirstElem, LastElem 498 LOGICAL :: BcMode 499 500 501 CALL Info('ComponentIntegralReduction','Performing reduction for component: '& 502 //TRIM(ListGetString(CompParams,'Name')),Level=10) 503 504 OperX = 0.0_dp 505 vol = 0.0_dp 506 integral = 0.0_dp 507 508 BcMode = .FALSE. 509 MasterEntities => ListGetIntegerArray( CompParams,'Master Bodies',Found ) 510 IF( .NOT. Found ) THEN 511 MasterEntities => ListGetIntegerArray( CompParams,'Master Boundaries',Found ) 512 BcMode = .TRUE. 513 END IF 514 515 IF(.NOT. Found ) THEN 516 CALL Warn('ComponentIntegralReduction',& 517 '> Master Bodies < or > Master Boundaries < not given') 518 RETURN 519 END IF 520 521 IF( BcMode ) THEN 522 FirstElem = Mesh % NumberOfBulkElements + 1 523 LastElem = Mesh % NumberOfBulkElements + Mesh % NumberOfBoundaryElements 524 ELSE 525 FirstElem = 1 526 LastElem = Mesh % NumberOfBulkElements 527 END IF 528 529 IF(.NOT. AllocationsDone ) THEN 530 n = Model % MaxElementNodes 531 ALLOCATE( ElementNodes % x(n), ElementNodes % y(n), ElementNodes % z(n) ) 532 AllocationsDone = .TRUE. 533 END IF 534 535 ElementalVar = ( Var % TYPE == Variable_on_nodes_on_elements ) 536 NoDofs = Var % Dofs 537 CoeffAtIp = 1.0_dp 538 539 DO t=FirstElem,LastElem 540 Element => Mesh % Elements(t) 541 IF( BcMode ) THEN 542 IF( ALL( MasterEntities /= Element % BoundaryInfo % Constraint ) ) CYCLE 543 ELSE 544 IF( ALL( MasterEntities /= Element % BodyId ) ) CYCLE 545 END IF 546 547 n = Element % TYPE % NumberOfNodes 548 NodeIndexes => Element % NodeIndexes 549 IF( ElementalVar ) THEN 550 DofIndexes => Element % DGIndexes 551 ELSE 552 DofIndexes => NodeIndexes 553 END IF 554 555 IF( ASSOCIATED( Var % Perm ) ) THEN 556 PermIndexes = Var % Perm( DofIndexes ) 557 ELSE 558 PermIndexes = DofIndexes 559 END IF 560 561 IF ( ANY(PermIndexes == 0 ) ) CYCLE 562 563 ElementNodes % x(1:n) = Mesh % Nodes % x(NodeIndexes(1:n)) 564 ElementNodes % y(1:n) = Mesh % Nodes % y(NodeIndexes(1:n)) 565 ElementNodes % z(1:n) = Mesh % Nodes % z(NodeIndexes(1:n)) 566 567 IF(GotCoeff) THEN 568 k = ListGetInteger( Model % Bodies( Element % BodyId ) % Values, & 569 'Material', Found ) 570 Coeff(1:n) = ListGetReal( Model % Materials(k) % Values, & 571 CoeffName, n, NodeIndexes(1:n) ) 572 END IF 573 574!------------------------------------------------------------------------------ 575! Numerical integration 576!------------------------------------------------------------------------------ 577 IntegStuff = GaussPoints( Element ) 578 579 DO i=1,IntegStuff % n 580 U = IntegStuff % u(i) 581 V = IntegStuff % v(i) 582 W = IntegStuff % w(i) 583!------------------------------------------------------------------------------ 584! Basis function values & derivatives at the integration point 585!------------------------------------------------------------------------------ 586 stat = ElementInfo( Element,ElementNodes,U,V,W,SqrtElementMetric, & 587 Basis,dBasisdx ) 588!------------------------------------------------------------------------------ 589! Coordinatesystem dependent info 590!------------------------------------------------------------------------------ 591 s = SqrtElementMetric * IntegStuff % s(i) 592 IF ( CurrentCoordinateSystem() /= Cartesian ) THEN 593 x = 2 * PI * SUM( ElementNodes % x(1:n)*Basis(1:n) ) 594 END IF 595 596 IF( GotCoeff ) THEN 597 CoeffAtIp = SUM( Coeff(1:n) * Basis(1:n) ) 598 END IF 599 vol = vol + S 600 601 SELECT CASE(OperName) 602 603 CASE ('volume') 604 integral = integral + coeffAtIp * S 605 606 CASE ('int','int mean') 607 func = SUM( Var % Values(PermIndexes(1:n)) * Basis(1:n) ) 608 integral = integral + S * coeffAtIp * func 609 610 CASE ('int abs','int abs mean') 611 func = ABS( SUM( Var % Values(PermIndexes(1:n)) * Basis(1:n) ) ) 612 integral = integral + S * coeffAtIp * func 613 614 CASE ('diffusive energy') 615 DO j = 1, 3 616 Grad(j) = SUM( dBasisdx(1:n,j) * Var % Values(PermIndexes(1:n) ) ) 617 END DO 618 integral = integral + s * CoeffAtIp * SUM( Grad * Grad ) 619 620 CASE ('convective energy') 621 func = SUM( Var % Values(PermIndexes(1:n)) * Basis(1:n) ) 622 623 IF(NoDofs == 1) THEN 624 func = SUM( Var % Values(PermIndexes(1:n)) * Basis(1:n) ) 625 integral = integral + s * coeffAtIp * func**2 626 ELSE 627 func = 0.0d0 628 DO j=1,NoDofs 629 func = SUM( Var % Values(NoDofs*(PermIndexes(1:n)-1)+j) * Basis(1:n) ) 630 integral = integral + s * coeffAtIp * func**2 631 END DO 632 END IF 633 634 CASE ('potential energy') 635 func = SUM( Var % Values(PermIndexes(1:n)) * Basis(1:n) ) 636 integral = integral + s * coeffAtIp * func 637 638 CASE DEFAULT 639 CALL Warn('ComponentIntegralReduction','Unknown operator') 640 641 END SELECT 642 643 END DO 644 645 END DO 646 647 integral = ParallelReduction( integral ) 648 649 SELECT CASE(OperName) 650 651 CASE ('volume') 652 operx = integral 653 654 CASE ('int') 655 operx = integral 656 657 CASE ('int abs') 658 operx = integral 659 660 CASE ('int mean') 661 vol = ParallelReduction( vol ) 662 operx = integral / vol 663 664 CASE ('int abs mean') 665 vol = ParallelReduction( vol ) 666 operx = integral / vol 667 668 CASE ('diffusive energy') 669 operx = 0.5d0 * integral 670 671 CASE ('convective energy') 672 operx = 0.5d0 * integral 673 674 CASE ('potential energy') 675 operx = integral 676 677 END SELECT 678 679 CALL Info('ComponentIntegralReduction','Reduction operator finished',Level=12) 680 681!------------------------------------------------------------------------------ 682 END FUNCTION ComponentIntegralReduction 683!------------------------------------------------------------------------------ 684 685 686!------------------------------------------------------------------------------ 687!> Each solver may include a list of dependent components that are updated 688!> after the solver (or the nonlinear iteration related to it) has been executed. 689!------------------------------------------------------------------------------ 690 SUBROUTINE UpdateDependentComponents( ComponentList ) 691 INTEGER, POINTER :: ComponentList(:) 692 693 INTEGER :: i,j,NoVar 694 CHARACTER(LEN=MAX_NAME_LEN) :: OperName, VarName, CoeffName, TmpOper 695 LOGICAL :: GotVar, GotOper, GotCoeff, VectorResult 696 TYPE(ValueList_t), POINTER :: CompParams 697 REAL(KIND=dp) :: ScalarVal, VectorVal(3), Power, Voltage 698 TYPE(Variable_t), POINTER :: Var 699 700 CALL Info('UpdateDepedentComponents','Updating Components to reflect new solution',Level=6) 701 702 DO j=1,CurrentModel % NumberOfComponents 703 704 IF( ALL( ComponentList /= j ) ) CYCLE 705 706 CALL Info('UpdateDepedentComponents','Updating component: '//TRIM(I2S(j))) 707 CompParams => CurrentModel % Components(j) % Values 708 709 710 NoVar = 0 711 DO WHILE( .TRUE. ) 712 NoVar = NoVar + 1 713 OperName = ListGetString( CompParams,'Operator '//TRIM(I2S(NoVar)), GotOper) 714 VarName = ListGetString( CompParams,'Variable '//TRIM(I2S(NoVar)), GotVar) 715 CoeffName = ListGetString( CompParams,'Coeffcient '//TRIM(I2S(NoVar)), GotCoeff) 716 717 IF(.NOT. GotVar .AND. GotOper .AND. OperName == 'electric resistance') THEN 718 VarName = 'Potential' 719 GotVar = .TRUE. 720 CALL Info('UpdateDependentComponents',& 721 'Defaulting field to > Potential < for operator: '//TRIM(OperName),Level=8) 722 END IF 723 724 IF(.NOT. (GotVar .AND. GotOper ) ) EXIT 725 726 Var => VariableGet( CurrentModel % Mesh % Variables, VarName ) 727 IF( .NOT. ASSOCIATED( Var ) ) THEN 728 CALL Info('UpdateDependentComponents','Variable not available: '//TRIM(VarName)) 729 CYCLE 730 END IF 731 VectorResult = .FALSE. 732 733 SELECT CASE( OperName ) 734 735 CASE('electric resistance') 736 IF(.NOT. GotCoeff ) THEN 737 CoeffName = 'electric conductivity' 738 GotCoeff = .TRUE. 739 END IF 740 TmpOper = 'diffusive energy' 741 Power = ComponentIntegralReduction(CurrentModel, CurrentModel % Mesh, CompParams, Var, & 742 TmpOper, CoeffName, GotCoeff ) 743 TmpOper = 'range' 744 Voltage = ComponentNodalReduction(CurrentModel, CurrentModel % Mesh, CompParams, Var, & 745 TmpOper ) 746 ScalarVal = Voltage**2 / Power 747 CALL ListAddConstReal( CompParams,'res: '//TRIM(OperName),ScalarVal ) 748 749 CASE ('sum','sum abs','min','max','min abs','max abs','range','mean','mean abs','variance') 750 ScalarVal = ComponentNodalReduction(CurrentModel, CurrentModel % Mesh, CompParams, Var, & 751 OperName ) 752 CALL ListAddConstReal( CompParams,'res: '//TRIM(OperName)//' '//TRIM(VarName),ScalarVal ) 753 754 CASE ('volume','int','int abs','int mean','int abs mean','diffusive energy',& 755 'convective energy','potential energy') 756 ScalarVal = ComponentIntegralReduction(CurrentModel, CurrentModel % Mesh, CompParams, Var, & 757 OperName, CoeffName, GotCoeff ) 758 CALL ListAddConstReal( CompParams,'res: '//TRIM(OperName)//' '//TRIM(VarName),ScalarVal ) 759 760 CASE('force') 761 CALL ComponentNodalForceReduction(CurrentModel, CurrentModel % Mesh, CompParams, Var, & 762 Force = VectorVal ) 763 VectorResult = .TRUE. 764 765 CASE('moment') 766 CALL ComponentNodalForceReduction(CurrentModel, CurrentModel % Mesh, CompParams, Var, & 767 Moment = VectorVal ) 768 VectorResult = .TRUE. 769 770 CASE('torque') 771 CALL ComponentNodalForceReduction(CurrentModel, CurrentModel % Mesh, CompParams, Var, & 772 Torque = ScalarVal ) 773 774 CASE DEFAULT 775 CALL Fatal('UpdateDependentComponents','Uknown operator: '//TRIM(OperName)) 776 777 END SELECT 778 779 IF( VectorResult ) THEN 780 DO i=1,3 781 WRITE( Message,'(A,ES15.6)') TRIM(OperName)//': '//TRIM(VarName)//' '& 782 //TRIM(I2S(i))//': ',ScalarVal 783 CALL Info('UpdateDependentComponents',Message,Level=5) 784 CALL ListAddConstReal( CompParams,'res: '//TRIM(OperName)//': '& 785 //TRIM(VarName)//' '//TRIM(I2S(i)),VectorVal(i) ) 786 END DO 787 ELSE 788 WRITE( Message,'(A,ES15.6)') & 789 'comp '//TRIM(I2S(j))//': '//TRIM(OperName)//': '//TRIM(VarName)//': ',ScalarVal 790 CALL Info('UpdateDependentComponents',Message,Level=5) 791 CALL ListAddConstReal( CurrentModel % Simulation, & 792 'res: comp '//TRIM(I2S(j))//': '//TRIM(OperName)//' '//TRIM(VarName),ScalarVal ) 793 END IF 794 795 END DO 796 END DO 797 798!------------------------------------------------------------------------------ 799 END SUBROUTINE UpdateDependentComponents 800!------------------------------------------------------------------------------ 801 802 803 END MODULE ComponentUtils 804!------------------------------------------------------------------------------ 805 806 807!> \} 808 809 810