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: 01 Oct 1996 34! * 35! ******************************************************************************/ 36 37 38!--------------------------------------------------------------------------------- 39!> Diffuse-convective local matrix computing in general euclidean coordinates. 40!--------------------------------------------------------------------------------- 41!> \ingroup ElmerLib 42!> \{ 43 44MODULE DiffuseConvectiveGeneral 45 46 USE Integration 47 USE MaterialModels 48 USE Differentials 49 50 IMPLICIT NONE 51 52 CONTAINS 53 54!------------------------------------------------------------------------------ 55!> Return element local matrices and RSH vector for diffusion-convection 56!> equation (genaral euclidean coordinate system): 57!------------------------------------------------------------------------------ 58 SUBROUTINE DiffuseConvectiveGenCompose( MassMatrix,StiffMatrix,ForceVector, & 59 LoadVector,NodalCT,NodalC0,NodalC1,NodalC2,PhaseChange,Temperature,Enthalpy,& 60 Ux,Uy,Uz,MUx,MUy, MUz, NodalViscosity,NodalDensity,NodalPressure, & 61 NodaldPressureDt,NodalPressureCoeff,Compressible, Stabilize,Element,n,Nodes ) 62!------------------------------------------------------------------------------ 63! 64! REAL(KIND=dp) :: MassMatrix(:,:) 65! OUTPUT: time derivative coefficient matrix 66! 67! REAL(KIND=dp) :: StiffMatrix(:,:) 68! OUTPUT: rest of the equation coefficients 69! 70! REAL(KIND=dp) :: ForceVector(:) 71! OUTPUT: RHS vector 72! 73! REAL(KIND=dp) :: LoadVector(:) 74! INPUT: 75! 76! REAL(KIND=dp) :: NodalCT,NodalC0,NodalC1 77! INPUT: Coefficient of the time derivative term, 0 degree term, and the 78! convection term respectively 79! 80! REAL(KIND=dp) :: NodalC2(:,:,:) 81! INPUT: Nodal values of the diffusion term coefficient tensor 82! 83! LOGICAL :: PhaseChange 84! INPUT: Do we model phase change here... 85! 86! REAL(KIND=dp) :: Temperature 87! INPUT: Temperature from previous iteration, needed if we model 88! phase change 89! 90! REAL(KIND=dp) :: Enthalpy 91! INPUT: Enthalpy from previous iteration, needed if we model 92! phase change 93! 94! REAL(KIND=dp) :: Ux(:),Uy(:),Uz(:) 95! INPUT: Nodal values of velocity components from previous iteration 96! used only if coefficient of the convection term (C1) is nonzero 97! 98! REAL(KIND=dp) :: NodalViscosity(:) 99! INPUT: Nodal values of the viscosity 100! 101! LOGICAL :: Stabilize 102! INPUT: Should stabilzation be used ? Used only if coefficient of the 103! convection term (C1) is nonzero 104! 105! TYPE(Element_t) :: Element 106! INPUT: Structure describing the element (dimension,nof nodes, 107! interpolation degree, etc...) 108! 109! TYPE(Nodes_t) :: Nodes 110! INPUT: Element node coordinates 111! 112!------------------------------------------------------------------------------ 113 114 REAL(KIND=dp), DIMENSION(:) :: ForceVector,Ux,Uy,Uz,MUx,MUy,MUz,LoadVector 115 REAL(KIND=dp), DIMENSION(:,:) :: MassMatrix,StiffMatrix 116 REAL(KIND=dp) :: NodalC0(:),NodalC1(:),NodalCT(:),NodalC2(:,:,:) 117 REAL(KIND=dp) :: Temperature(:),Enthalpy(:),NodalViscosity(:), & 118 NodaldPressuredt(:),NodalPressureCoeff(:),NodalPressure(:),dT, NodalDensity(:) 119 120 LOGICAL :: Stabilize,PhaseChange,Compressible 121 122 INTEGER :: n 123 124 TYPE(Nodes_t) :: Nodes 125 TYPE(Element_t), POINTER :: Element 126 127 128!------------------------------------------------------------------------------ 129! Local variables 130!------------------------------------------------------------------------------ 131! 132 REAL(KIND=dp) :: ddBasisddx(n,3,3),dNodalBasisdx(n,n,3) 133 REAL(KIND=dp) :: Basis(2*n) 134 REAL(KIND=dp) :: dBasisdx(2*n,3),detJ 135 136 REAL(KIND=dp) :: Velo(3),Force 137 138 REAL(KIND=dp) :: A,M 139 REAL(KIND=dp) :: Load 140 141 REAL(KIND=dp) :: VNorm,hK,mK 142 REAL(KIND=dp) :: Lambda=1.0,Pe,Pe1,Pe2,C00,Tau,Delta,x,y,z 143 144 INTEGER :: i,j,k,c,p,q,t,dim,N_Integ,NBasis 145 146 REAL(KIND=dp) :: s,u,v,w,dEnth,dTemp,Viscosity,Pressure,pCoeff,DivVelo,dVelodx(3,3) 147 148 REAL(KIND=dp) :: SqrtMetric,Metric(3,3),Symb(3,3,3),dSymb(3,3,3,3) 149 150 REAL(KIND=dp), DIMENSION(:), POINTER :: U_Integ,V_Integ,W_Integ,S_Integ 151 152 REAL(KIND=dp) :: C0,CT,CL,C1,C2(3,3),dC2dx(3,3,3),SU(n),SW(n),Density 153 154 TYPE(GaussIntegrationPoints_t), TARGET :: IntegStuff 155 156 LOGICAL :: stat,CylindricSymmetry,Convection,ConvectAndStabilize,Bubbles, & 157 FrictionHeat, Found 158 TYPE(ValueList_t), POINTER :: BodyForce, Material 159 160 LOGICAL :: GotCondModel 161 162!------------------------------------------------------------------------------ 163 164 CylindricSymmetry = (CurrentCoordinateSystem() == CylindricSymmetric .OR. & 165 CurrentCoordinateSystem() == AxisSymmetric) 166 167 168 IF ( CylindricSymmetry ) THEN 169 dim = 3 170 ELSE 171 dim = CoordinateSystemDimension() 172 END IF 173 n = element % Type % NumberOfNodes 174 175 ForceVector = 0.0D0 176 StiffMatrix = 0.0D0 177 MassMatrix = 0.0D0 178 Load = 0.0D0 179 180 Convection = ANY( NodalC1 /= 0.0d0 ) 181 NBasis = n 182 Bubbles = .FALSE. 183 IF ( Convection .AND. .NOT. Stabilize ) THEN 184 NBasis = 2*n 185 Bubbles = .TRUE. 186 END IF 187 188 Material => GetMaterial() 189 GotCondModel = ListCheckPresent( Material,'Heat Conductivity Model') 190 191!------------------------------------------------------------------------------ 192! Integration stuff 193!------------------------------------------------------------------------------ 194 IF ( Bubbles ) THEN 195 IntegStuff = GaussPoints( element, Element % Type % GaussPoints2 ) 196 ELSE 197 IntegStuff = GaussPoints( element ) 198 END IF 199 U_Integ => IntegStuff % u 200 V_Integ => IntegStuff % v 201 W_Integ => IntegStuff % w 202 S_Integ => IntegStuff % s 203 N_Integ = IntegStuff % n 204 205!------------------------------------------------------------------------------ 206! Stabilization parameters: hK, mK (take a look at Franca et.al.) 207! If there is no convection term we don't need stabilization. 208!------------------------------------------------------------------------------ 209 ConvectAndStabilize = .FALSE. 210 IF ( Stabilize .AND. ANY(NodalC1 /= 0.0D0) ) THEN 211 ConvectAndStabilize = .TRUE. 212 hK = element % hK 213 mK = element % StabilizationMK 214 dNodalBasisdx = 0._dp 215 DO p=1,n 216 u = Element % Type % NodeU(p) 217 v = Element % Type % NodeV(p) 218 w = Element % Type % NodeW(p) 219 stat = ElementInfo( Element, Nodes, u,v,w, detJ, Basis, dBasisdx ) 220 dNodalBasisdx(1:n,p,:) = dBasisdx(1:n,:) 221 END DO 222 END IF 223 224!------------------------------------------------------------------------------ 225 BodyForce => GetBodyForce() 226 FrictionHeat = .FALSE. 227 IF (ASSOCIATED(BodyForce)) & 228 FrictionHeat = GetLogical( BodyForce, 'Friction Heat', Found ) 229!------------------------------------------------------------------------------ 230! Now we start integrating 231!------------------------------------------------------------------------------ 232 DO t=1,N_Integ 233 234 u = U_Integ(t) 235 v = V_Integ(t) 236 w = W_Integ(t) 237!------------------------------------------------------------------------------ 238! Basis function values & derivatives at the integration point 239!------------------------------------------------------------------------------ 240 stat = ElementInfo( Element,Nodes,u,v,w,detJ, & 241 Basis,dBasisdx, Bubbles=Bubbles ) 242 243!------------------------------------------------------------------------------ 244! Coordinatesystem dependent info 245!------------------------------------------------------------------------------ 246 IF ( CurrentCoordinateSystem()/= Cartesian ) THEN 247 x = SUM( nodes % x(1:n)*Basis(1:n) ) 248 y = SUM( nodes % y(1:n)*Basis(1:n) ) 249 z = SUM( nodes % z(1:n)*Basis(1:n) ) 250 END IF 251 252 CALL CoordinateSystemInfo( Metric,SqrtMetric,Symb,dSymb,x,y,z ) 253 254 s = SqrtMetric * detJ * S_Integ(t) 255!------------------------------------------------------------------------------ 256! Coefficient of the convection and time derivative terms at the 257! integration point 258!------------------------------------------------------------------------------ 259 C0 = SUM( NodalC0(1:n)*Basis(1:n) ) 260 CT = SUM( NodalCT(1:n)*Basis(1:n) ) 261 C1 = SUM( NodalC1(1:n)*Basis(1:n) ) 262!------------------------------------------------------------------------------ 263! Compute effective heatcapacity, if modelling phase change, 264! at the integration point. 265! NOTE: This is for heat equation only, not generally for diff.conv. equ. 266!------------------------------------------------------------------------------ 267 IF ( PhaseChange ) THEN 268 dEnth = 0.0D0 269 dTemp = 0.0D0 270 DO i=1,3 271 dEnth = dEnth + SUM( Enthalpy(1:n) * dBasisdx(1:n,i) )**2 272 dTemp = dTemp + SUM( Temperature(1:n) * dBasisdx(1:n,i) )**2 273 END DO 274 275 CL = SQRT( dEnth / dTemp ) 276 277 CT = CT + CL 278 END IF 279!------------------------------------------------------------------------------ 280! Coefficient of the diffusion term & its derivatives at the 281! integration point 282!------------------------------------------------------------------------------ 283 Density = SUM( NodalDensity(1:n) * Basis(1:n) ) 284 DO i=1,dim 285 DO j=1,dim 286 C2(i,j) = SQRT(Metric(i,i)) * SQRT(Metric(j,j)) * & 287 SUM( NodalC2(i,j,1:n) * Basis(1:n) ) 288 END DO 289 END DO 290 291 IF( GotCondModel ) THEN 292 DO i=1,dim 293 C2(i,i) = EffectiveConductivity( C2(i,i), Density, Element, & 294 Temperature, UX,UY,UZ, Nodes, n, n, u, v, w ) 295 END DO 296 END IF 297!------------------------------------------------------------------------------ 298! If there's no convection term we don't need the velocities, and 299! also no need for stabilization 300!------------------------------------------------------------------------------ 301 Convection = .FALSE. 302 IF ( C1 /= 0.0D0 ) THEN 303 Convection = .TRUE. 304 IF ( PhaseChange ) C1 = C1 + CL 305!------------------------------------------------------------------------------ 306! Velocity and pressure (deviation) from previous iteration 307! at the integration point 308!------------------------------------------------------------------------------ 309 Velo = 0.0D0 310 Velo(1) = SUM( (Ux(1:n)-MUx(1:n))*Basis(1:n) ) 311 Velo(2) = SUM( (Uy(1:n)-MUy(1:n))*Basis(1:n) ) 312 IF ( dim > 2 .AND. CurrentCoordinateSystem()/= AxisSymmetric ) THEN 313 Velo(3) = SUM( (Uz(1:n)-MUz(1:n))*Basis(1:n) ) 314 END IF 315 316 IF ( Compressible ) THEN 317 Pressure = SUM( NodalPressure(1:n)*Basis(1:n) ) 318 319 dVelodx = 0.0D0 320 DO i=1,3 321 dVelodx(1,i) = SUM( Ux(1:n)*dBasisdx(1:n,i) ) 322 dVelodx(2,i) = SUM( Uy(1:n)*dBasisdx(1:n,i) ) 323 IF ( dim > 2 .AND. CurrentCoordinateSystem()/= AxisSymmetric ) & 324 dVelodx(3,i) = SUM( Uz(1:n)*dBasisdx(1:n,i) ) 325 END DO 326 327 DivVelo = 0.0D0 328 DO i=1,dim 329 DivVelo = DivVelo + dVelodx(i,i) 330 END DO 331 IF ( CurrentCoordinateSystem()>= Cylindric .AND. & 332 CurrentCoordinateSystem()<= AxisSymmetric ) THEN 333! Cylindrical coordinates 334 DivVelo = DivVelo + Velo(1)/x 335 ELSE 336! General coordinate system 337 DO i=1,dim 338 DO j=i,dim 339 DivVelo = DivVelo + Velo(j)*Symb(i,j,i) 340 END DO 341 END DO 342 END IF 343 END IF 344 345!------------------------------------------------------------------------------ 346! Stabilization parameters... 347!------------------------------------------------------------------------------ 348 IF ( Stabilize ) THEN 349! VNorm = SQRT( SUM(Velo(1:dim)**2) ) 350 351 Vnorm = 0.0D0 352 DO i=1,dim 353 Vnorm = Vnorm + Velo(i)*Velo(i) / Metric(i,i) 354 END DO 355 Vnorm = SQRT( Vnorm ) 356 357#if 1 358 Pe = MIN(1.0D0,mK*hK*C1*VNorm/(2*ABS(C2(1,1)))) 359 360 Tau = 0.0D0 361 IF ( VNorm /= 0.0D0 ) THEN 362 Tau = hK * Pe / (2 * C1 * VNorm) 363 END IF 364#else 365 C00 = C0 366 IF ( dT > 0 ) C00 = C0 + CT 367 368 Pe1 = 0.0d0 369 IF ( C00 > 0 ) THEN 370 Pe1 = 2 * ABS(C2(1,1)) / ( mK * C00 * hK**2 ) 371 Pe1 = C00 * hK**2 * MAX( 1.0d0, Pe1 ) 372 ELSE 373 Pe1 = 2 * ABS(C2(1,1)) / mK 374 END IF 375 376 Pe2 = 0.0d0 377 IF ( C2(1,1) /= 0.0d0 ) THEN 378 Pe2 = ( mK * C1 * VNorm * hK ) / ABS(C2(1,1)) 379 Pe2 = 2*ABS(C2(1,1)) * MAX( 1.0d0, Pe2 ) / mK 380 ELSE 381 Pe2 = 2 * hK * C1 * VNorm 382 END IF 383 384 Tau = hk**2 / ( Pe1 + Pe2 ) 385#endif 386 387!------------------------------------------------------------------------------ 388 DO i=1,dim 389 DO j=1,dim 390 DO k=1,3 391 dC2dx(i,j,k) = SQRT(Metric(i,i))*SQRT(Metric(j,j))* & 392 SUM(NodalC2(i,j,1:n)*dBasisdx(1:n,k)) 393 END DO 394 END DO 395 END DO 396!------------------------------------------------------------------------------ 397! Compute residual & stabilization weight vectors 398!------------------------------------------------------------------------------ 399 DO p=1,n 400 SU(p) = C0 * Basis(p) 401 DO i = 1,dim 402 SU(p) = SU(p) + C1 * dBasisdx(p,i) * Velo(i) 403 IF ( Element % Type % BasisFunctionDegree <= 1 ) CYCLE 404 405 DO j=1,dim 406 SU(p) = SU(p) - dC2dx(i,j,j) * dBasisdx(p,i) 407 SU(p) = SU(p) - C2(i,j) * SUM(dNodalBasisdx(p,1:n,i)*dBasisdx(1:n,j)) 408 DO k=1,dim 409 SU(p) = SU(p) + C2(i,j) * Symb(i,j,k) * dBasisdx(p,k) 410 SU(p) = SU(p) - C2(i,k) * Symb(k,j,j) * dBasisdx(p,i) 411 SU(p) = SU(p) - C2(k,j) * Symb(k,j,i) * dBasisdx(p,i) 412 END DO 413 END DO 414 END DO 415 416 SW(p) = C0 * Basis(p) 417 418 DO i = 1,dim 419 SW(p) = SW(p) + C1 * dBasisdx(p,i) * Velo(i) 420 IF ( Element % Type % BasisFunctionDegree <= 1 ) CYCLE 421 422 DO j=1,dim 423 SW(p) = SW(p) - dC2dx(i,j,j) * dBasisdx(p,i) 424 SW(p) = SW(p) - C2(i,j) * SUM(dNodalBasisdx(p,1:n,i)*dBasisdx(1:n,j)) 425 DO k=1,dim 426 SW(p) = SW(p) + C2(i,j) * Symb(i,j,k) * dBasisdx(p,k) 427 SW(p) = SW(p) - C2(i,k) * Symb(k,j,j) * dBasisdx(p,i) 428 SW(p) = SW(p) - C2(k,j) * Symb(k,j,i) * dBasisdx(p,i) 429 END DO 430 END DO 431 END DO 432 END DO 433 END IF 434 END IF 435!------------------------------------------------------------------------------ 436! Loop over basis functions of both unknowns and weights 437!------------------------------------------------------------------------------ 438 DO p=1,NBasis 439 DO q=1,NBasis 440!------------------------------------------------------------------------------ 441! The diffusive-convective equation without stabilization 442!------------------------------------------------------------------------------ 443 M = CT * Basis(q) * Basis(p) 444 A = C0 * Basis(q) * Basis(p) 445 DO i=1,dim 446 DO j=1,dim 447 A = A + C2(i,j) * dBasisdx(q,i) * dBasisdx(p,j) 448 END DO 449 END DO 450 451 IF ( Convection ) THEN 452 DO i=1,dim 453 A = A + C1 * Velo(i) * dBasisdx(q,i) * Basis(p) 454 END DO 455 456!------------------------------------------------------------------------------ 457! Next we add the stabilization... 458!------------------------------------------------------------------------------ 459 IF ( Stabilize ) THEN 460 A = A + Tau * SU(q) * SW(p) 461 M = M + Tau * CT * Basis(q) * SW(p) 462 END IF 463 END IF 464 465 StiffMatrix(p,q) = StiffMatrix(p,q) + s * A 466 MassMatrix(p,q) = MassMatrix(p,q) + s * M 467 END DO 468 END DO 469 470!------------------------------------------------------------------------------ 471! Force at the integration point 472!------------------------------------------------------------------------------ 473 Force = SUM( LoadVector(1:n)*Basis(1:n) ) + & 474 JouleHeat( Element,Nodes,u,v,w,n ) 475 IF ( Convection ) THEN 476! IF ( Compressible ) Force = Force - Pressure * DivVelo 477 Pcoeff = SUM( NodalPressureCoeff(1:n) * Basis(1:n) ) 478 IF ( Pcoeff /= 0.0_dp ) THEN 479 Force = Force + Pcoeff*SUM( NodalDPressureDt(1:n) * Basis(1:n) ) 480 DO i=1,dim 481 Force = Force + Pcoeff*Velo(i)*SUM( NodalPressure(1:n)*dBasisdx(1:n,i) ) 482 END DO 483 END IF 484 485 IF ( FrictionHeat ) THEN 486 Viscosity = SUM( NodalViscosity(1:n) * Basis(1:n) ) 487 Viscosity = EffectiveViscosity( Viscosity, Density, Ux, Uy, Uz, & 488 Element, Nodes, n, n, u, v, w, LocalIP=t ) 489 IF ( Viscosity > 0.0d0 ) THEN 490 IF ( .NOT.Compressible ) THEN 491 dVelodx = 0.0D0 492 DO i=1,3 493 dVelodx(1,i) = SUM( Ux(1:n)*dBasisdx(1:n,i) ) 494 dVelodx(2,i) = SUM( Uy(1:n)*dBasisdx(1:n,i) ) 495 IF ( dim > 2 .AND. CurrentCoordinateSystem()/= AxisSymmetric ) & 496 dVelodx(3,i) = SUM( Uz(1:n)*dBasisdx(1:n,i) ) 497 END DO 498 END IF 499 Force = Force + 0.5d0 * Viscosity * & 500 SecondInvariant( Velo,dVelodx,Metric,Symb ) 501 END IF 502 END IF 503 END IF 504!------------------------------------------------------------------------------ 505! The righthand side... 506!------------------------------------------------------------------------------ 507 DO p=1,NBasis 508 Load = Basis(p) 509 510 IF ( ConvectAndStabilize ) THEN 511 Load = Load + Tau * SW(p) 512 END IF 513 514 ForceVector(p) = ForceVector(p) + s * Load * Force 515 END DO 516 517 END DO 518!------------------------------------------------------------------------------ 519 END SUBROUTINE DiffuseConvectiveGenCompose 520!------------------------------------------------------------------------------ 521 522 523!------------------------------------------------------------------------------ 524!> Return element local matrices and RHS vector for boundary conditions 525!> of diffusion convection equation. 526!------------------------------------------------------------------------------ 527 SUBROUTINE DiffuseConvectiveGenBoundary( BoundaryMatrix,BoundaryVector, & 528 LoadVector,NodalAlpha,Element,n,Nodes) 529!------------------------------------------------------------------------------ 530! 531! REAL(KIND=dp) :: BoundaryMatrix(:,:) 532! OUTPUT: coefficient matrix if equations 533! 534! REAL(KIND=dp) :: BoundaryVector(:) 535! OUTPUT: RHS vector 536! 537! REAL(KIND=dp) :: LoadVector(:) 538! INPUT: coefficient of the force term 539! 540! REAL(KIND=dp) :: NodalAlpha 541! INPUT: coefficient for temperature dependent term 542! 543! TYPE(Element_t) :: Element 544! INPUT: Structure describing the element (dimension,nof nodes, 545! interpolation degree, etc...) 546! 547! INTEGER :: n 548! INPUT: Number of element nodes 549! 550! TYPE(Nodes_t) :: Nodes 551! INPUT: Element node coordinates 552! 553!------------------------------------------------------------------------------ 554 555 REAL(KIND=dp) :: BoundaryMatrix(:,:),BoundaryVector(:) 556 REAL(KIND=dp) :: LoadVector(:),NodalAlpha(:) 557 TYPE(Nodes_t) :: Nodes 558 TYPE(Element_t),POINTER :: Element 559 560 INTEGER :: n 561!------------------------------------------------------------------------------ 562 563 REAL(KIND=dp) :: ddBasisddx(n,3,3) 564 REAL(KIND=dp) :: Basis(n) 565 REAL(KIND=dp) :: dBasisdx(n,3),detJ 566 567 REAL(KIND=dp) :: u,v,w,s,x,y,z 568 REAL(KIND=dp) :: Force,Alpha 569 REAL(KIND=dp), POINTER :: U_Integ(:),V_Integ(:),W_Integ(:),S_Integ(:) 570 571 REAL(KIND=dp) :: SqrtMetric,Metric(3,3),Symb(3,3,3),dSymb(3,3,3,3),normal(3) 572 573 INTEGER :: i,t,q,p,N_Integ 574 575 LOGICAL :: stat,CylindricSymmetry 576 577 TYPE(GaussIntegrationPoints_t), TARGET :: IntegStuff 578!------------------------------------------------------------------------------ 579 580 BoundaryVector = 0.0D0 581 BoundaryMatrix = 0.0D0 582 583!------------------------------------------------------------------------------ 584! Integration stuff 585!------------------------------------------------------------------------------ 586 IntegStuff = GaussPoints( element ) 587 U_Integ => IntegStuff % u 588 V_Integ => IntegStuff % v 589 W_Integ => IntegStuff % w 590 S_Integ => IntegStuff % s 591 N_Integ = IntegStuff % n 592 593!------------------------------------------------------------------------------ 594! Now we start integrating 595!------------------------------------------------------------------------------ 596 DO t=1,N_Integ 597 u = U_Integ(t) 598 v = V_Integ(t) 599 w = W_Integ(t) 600 601!------------------------------------------------------------------------------ 602! Basis function values & derivatives at the integration point 603!------------------------------------------------------------------------------ 604 stat = ElementInfo( Element,Nodes,u,v,w,detJ, & 605 Basis,dBasisdx ) 606 607 s = S_Integ(t) * detJ 608!------------------------------------------------------------------------------ 609! Coordinatesystem dependent info 610!------------------------------------------------------------------------------ 611 IF ( CurrentCoordinateSystem()/= Cartesian ) THEN 612 x = SUM( Nodes % x(1:n)*Basis ) 613 y = SUM( Nodes % y(1:n)*Basis ) 614 z = SUM( Nodes % z(1:n)*Basis ) 615 s = s * CoordinateSqrtMetric( x,y,z ) 616 END IF 617!------------------------------------------------------------------------------ 618! Basis function values at the integration point 619!------------------------------------------------------------------------------ 620 Alpha = SUM( NodalAlpha(1:n)*Basis ) 621 Force = SUM( LoadVector(1:n)*Basis ) 622 623 DO p=1,N 624 DO q=1,N 625 BoundaryMatrix(p,q) = BoundaryMatrix(p,q) + & 626 s * Alpha * Basis(q) * Basis(p) 627 END DO 628 END DO 629 630 DO q=1,N 631 BoundaryVector(q) = BoundaryVector(q) + s * Basis(q) * Force 632 END DO 633 END DO 634 END SUBROUTINE DiffuseConvectiveGenBoundary 635!------------------------------------------------------------------------------ 636 637END MODULE DiffuseConvectiveGeneral 638 639!> \} 640