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 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! * 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: 08 Jun 1997 34! * 35! *****************************************************************************/ 36 37!> \ingroup Solvers 38 39!------------------------------------------------------------------------------ 40SUBROUTINE MagnetoDynamicsCalcFields_Init0(Model,Solver,dt,Transient) 41!------------------------------------------------------------------------------ 42 USE MagnetoDynamicsUtils 43 44 IMPLICIT NONE 45!------------------------------------------------------------------------------ 46 TYPE(Solver_t), TARGET :: Solver 47 TYPE(Model_t) :: Model 48 REAL(KIND=dp) :: dt 49 LOGICAL :: Transient 50!------------------------------------------------------------------------------ 51 CHARACTER(LEN=MAX_NAME_LEN) :: sname,pname 52 LOGICAL :: Found, ElementalFields, RealField, FoundVar, Hcurl 53 INTEGER, POINTER :: Active(:) 54 INTEGER :: mysolver,i,j,k,l,n,m,vDOFs, soln,pIndex 55 TYPE(ValueList_t), POINTER :: SolverParams, DGSolverParams 56 TYPE(Solver_t), POINTER :: Solvers(:), PSolver 57 58 ! This is really using DG so we don't need to make any dirty tricks to create DG fields 59 ! as is done in this initialization. 60 SolverParams => GetSolverParams() 61 62 ! The only purpose of this parsing of the variable name is to identify 63 ! whether the field is real or complex. As the variable has not been 64 ! created at this stage we have to do some dirty parsing. 65 pname = GetString(SolverParams, 'Potential variable', Found) 66 vdofs = 0 67 pIndex = 0 68 FoundVar = .FALSE. 69 70 IF( Found ) THEN 71 DO i=1,Model % NumberOfSolvers 72 sname = GetString(Model % Solvers(i) % Values, 'Variable', Found) 73 74 J=INDEX(sname,'[')-1 75 IF ( j<=0 ) j=LEN_TRIM(sname) 76 IF ( sname(1:j) == pname(1:LEN_TRIM(pname)) )THEN 77 k = 0 78 vDofs = 0 79 j=INDEX(sname,':') 80 DO WHILE(j>0) 81 Vdofs=Vdofs+ICHAR(sname(j+k+1:j+k+1))-ICHAR('0') 82 k = k+j 83 IF(k<LEN(sname)) j=INDEX(sname(k+1:),':') 84 END DO 85 pIndex = i 86 FoundVar = .TRUE. 87 EXIT 88 END IF 89 END DO 90 91 IF(.NOT. FoundVar ) THEN 92 CALL Fatal('MagnetoDynamicsCalcFields_Init0','Could not find solver for variable: '//TRIM(sname)) 93 END IF 94 END IF 95 96 97 ! When we created the case for GUI where "av" is not given in sif then it is impossible to 98 ! determine from the variable declaration what kind of solver we have. 99 IF( .NOT. FoundVar ) THEN 100 Hcurl = .FALSE. 101 DO i=Model % NumberOfSolvers,1,-1 102 sname = GetString(Model % Solvers(i) % Values, 'Procedure', Found) 103 104 j = INDEX( sname,'WhitneyAVHarmonicSolver') 105 IF( j > 0 ) THEN 106 Hcurl = .TRUE. 107 vDofs = 2 108 EXIT 109 END IF 110 111 j = INDEX( sname,'MagnetoDynamics2DHarmonic') 112 IF( j > 0 ) THEN 113 Vdofs = 2 114 EXIT 115 END IF 116 117 j = INDEX( sname,'WhitneyAVSolver') 118 IF( j > 0 ) THEN 119 Hcurl = .TRUE. 120 vDofs = 1 121 EXIT 122 END IF 123 124 j = INDEX( sname,'MagnetoDynamics2D') 125 IF( j > 0 ) THEN 126 Vdofs = 1 127 EXIT 128 END IF 129 END DO 130 131 IF( Vdofs == 0 ) THEN 132 CALL Fatal('MagnetoDynamicsCalcFields_Init0','Could not determine target variable type (real or complex)') 133 END IF 134 pIndex = i 135 END IF 136 137 RealField = ( Vdofs /= 2 ) 138 IF( RealField ) THEN 139 CALL Info('MagnetoDynamicsCalcFields_Init0','The target solver seems to be real valued',Level=12) 140 ELSE 141 CALL Info('MagnetoDynamicsCalcFields_Init0','The target solver seems to be complex valued',Level=12) 142 END IF 143 144 CALL ListAddNewLogical( SolverParams, 'Target Variable Real Field', RealField ) 145 CALL Info('MagnetoDynamicsCalcFields_Init0','Target Variable Solver Index: '& 146 //TRIM(I2S(pIndex)),Level=12) 147 CALL ListAddNewInteger( SolverParams, 'Target Variable Solver Index', pIndex ) 148 149!------------------------------------------------------------------------------ 150END SUBROUTINE MagnetoDynamicsCalcFields_Init0 151!------------------------------------------------------------------------------ 152 153 154!------------------------------------------------------------------------------ 155!> \ingroup Solvers 156!------------------------------------------------------------------------------ 157SUBROUTINE MagnetoDynamicsCalcFields_Init(Model,Solver,dt,Transient) 158!------------------------------------------------------------------------------ 159 USE MagnetoDynamicsUtils 160 161 IMPLICIT NONE 162!------------------------------------------------------------------------------ 163 TYPE(Solver_t) :: Solver 164 TYPE(Model_t) :: Model 165 REAL(KIND=dp) :: dt 166 LOGICAL :: Transient 167!------------------------------------------------------------------------------ 168 169 INTEGER :: i 170 LOGICAL :: Found, FluxFound, NodalFields, ElementalFields, & 171 RealField, ComplexField, LorentzConductivity 172 TYPE(ValueList_t), POINTER :: EQ, SolverParams 173 174 LorentzConductivity = ListCheckPrefixAnyBodyForce(Model, "Angular Velocity") .or. & 175 ListCheckPrefixAnyBodyForce(Model, "Lorentz Velocity") 176 177 IF(.NOT.ASSOCIATED(Solver % Values)) Solver % Values=>ListAllocate() 178 SolverParams => GetSolverParams() 179 180 ! Inherit this from the _init0 solver. Hence we know it must exist! 181 RealField = ListGetLogical( SolverParams,'Target Variable Real Field') 182 ComplexField = .NOT. RealField 183 184 CALL ListAddString( SolverParams, 'Variable', '-nooutput hr_dummy' ) 185 186 CALL ListAddLogical( SolverParams, 'Linear System refactorize', .FALSE.) 187 188 ! add these in the beginning, so that SaveData sees these existing, even 189 ! if executed before the actual computations... 190 ! ----------------------------------------------------------------------- 191 CALL ListAddConstReal(Model % Simulation,'res: Eddy current power',0._dp) 192 193 IF( ListGetLogical( SolverParams,'Separate Magnetic Energy',Found ) ) THEN 194 CALL ListAddConstReal(Model % Simulation,'res: Electric Field Energy',0._dp) 195 CALL ListAddConstReal(Model % Simulation,'res: Magnetic Field Energy',0._dp) 196 ELSE 197 CALL ListAddConstReal(Model % Simulation,'res: ElectroMagnetic Field Energy',0._dp) 198 END IF 199 200 IF (GetLogical(SolverParams,'Show Angular Frequency',Found)) & 201 CALL ListAddConstReal(Model % Simulation,'res: Angular Frequency',0._dp) 202 203 ! add these in the beginning only if the Magnetix Flux Average is computed 204 ! ------------------------------------------------------------------------- 205 IF (ListGetLogicalAnyBC( Model,'Magnetic Flux Average')) THEN 206 CALL ListAddConstReal( Model % Simulation,'res: Magnetic Flux Average', 0._dp) 207 CALL ListAddConstReal(Model % Simulation, & 208 'res: Magnetic Flux Density Average',0._dp) 209 210 IF (.NOT. RealField ) THEN 211 CALL ListAddConstReal(Model % Simulation,'res: Magnetic Flux im Average',0._dp) 212 CALL ListAddConstReal(Model % Simulation, & 213 'res: Magnetic Flux Density im Average', 0._dp ) 214 END IF 215 216 CALL ListAddConstReal(Model % Simulation,'res: Magnetic Flux Area',0._dp) 217 END IF 218 219 NodalFields = .NOT. GetLogical( SolverParams, 'Skip Nodal Fields', Found) 220 IF(.NOT. Found ) NodalFields = GetLogical( SolverParams, 'Calculate Nodal Fields', Found) 221 IF(.NOT. Found ) NodalFields = .TRUE. 222 223 i=1 224 DO WHILE(.TRUE.) 225 IF ( .NOT. ListCheckPresent(SolverParams,"Exported Variable "//TRIM(i2s(i))) ) EXIT 226 i = i + 1 227 END DO 228 i = i - 1 229 230 IF( NodalFields ) THEN 231 i = i + 1 232 IF ( RealField ) THEN 233 CALL ListAddString( SolverParams, "Exported Variable "//TRIM(i2s(i)), & 234 "Magnetic Flux Density[Magnetic Flux Density:3]" ) 235 ELSE 236 CALL ListAddString( SolverParams, "Exported Variable "//TRIM(i2s(i)), & 237 "Magnetic Flux Density[Magnetic Flux Density re:3 Magnetic Flux Density im:3]" ) 238 END IF 239 240 IF (GetLogical(SolverParams,'Calculate Magnetic Vector Potential',Found)) THEN 241 i = i + 1 242 IF ( RealField ) THEN 243 CALL ListAddString( SolverParams, "Exported Variable "//TRIM(i2s(i)), & 244 "Magnetic Vector Potential[Magnetic Vector Potential:3]" ) 245 ELSE 246 CALL ListAddString( SolverParams, "Exported Variable "//TRIM(i2s(i)), & 247 "Magnetic Vector Potential[Magnetic Vector Potential re:3 Magnetic Vector Potential im:3]") 248 END IF 249 END IF 250 251 IF (GetLogical(SolverParams,'Calculate Magnetic Field Strength',Found)) THEN 252 i = i + 1 253 IF ( RealField ) THEN 254 CALL ListAddString( SolverParams, "Exported Variable "//TRIM(i2s(i)), & 255 "Magnetic Field Strength[Magnetic Field Strength:3]" ) 256 ELSE 257 CALL ListAddString( SolverParams, "Exported Variable "//TRIM(i2s(i)), & 258 "Magnetic Field Strength[Magnetic Field Strength re:3 Magnetic Field Strength im:3]") 259 END IF 260 END IF 261 262 IF (GetLogical(SolverParams,'Calculate JxB',Found)) THEN 263 i = i + 1 264 IF ( RealField ) THEN 265 CALL ListAddString( SolverParams, "Exported Variable "//TRIM(i2s(i)), & 266 "JxB[JxB:3]" ) 267 ELSE 268 CALL ListAddString( SolverParams, "Exported Variable "//TRIM(i2s(i)), & 269 "JxB[JxB re:3 JxB im:3]") 270 END IF 271 END IF 272 273 IF ( GetLogical( SolverParams, 'Calculate Maxwell Stress', Found ) ) THEN 274 i = i + 1 275 IF ( RealField ) THEN 276 CALL ListAddString( SolverParams, "Exported Variable "//TRIM(i2s(i)), & 277 "Maxwell Stress[Maxwell Stress:6]" ) 278 ELSE 279 CALL ListAddString( SolverParams, "Exported Variable "//TRIM(i2s(i)), & 280 "Maxwell Stress[Maxwell Stress re:6 Maxwell Stress im:6]" ) 281 END IF 282 END IF 283 284 IF ( GetLogical( SolverParams, 'Calculate Current Density', Found ) ) THEN 285 i = i + 1 286 IF ( RealField ) THEN 287 CALL ListAddString( SolverParams, "Exported Variable "//TRIM(i2s(i)), & 288 "Current Density[Current Density:3]" ) 289 ELSE 290 CALL ListAddString( SolverParams, "Exported Variable "//TRIM(i2s(i)), & 291 "Current Density[Current Density re:3 Current Density im:3]" ) 292 END IF 293 END IF 294 295 IF ( GetLogical( SolverParams, 'Calculate Joule Heating', Found ) ) THEN 296 i = i + 1 297 CALL ListAddString( SolverParams, "Exported Variable "//TRIM(i2s(i)), & 298 "Joule Heating" ) 299 END IF 300 301 IF ( GetLogical( SolverParams, 'Calculate Harmonic Loss', Found ) ) THEN 302 IF( RealField ) THEN 303 CALL Warn('MagnetcDynamicsCalcFields',& 304 'Harmonic loss computation only available for complex systems!') 305 ELSE 306 i = i + 1 307 CALL ListAddString( SolverParams, "Exported Variable "//TRIM(i2s(i)), & 308 "Harmonic Loss Linear" ) 309 i = i + 1 310 CALL ListAddString( SolverParams, "Exported Variable "//TRIM(i2s(i)), & 311 "Harmonic Loss Quadratic" ) 312 END IF 313 END IF 314 315 IF ( Transient .OR. .NOT. RealField .OR. LorentzConductivity) THEN 316 IF ( GetLogical( SolverParams, 'Calculate Electric Field', Found ) ) THEN 317 i = i + 1 318 IF ( RealField ) THEN 319 CALL ListAddString( SolverParams, "Exported Variable "//TRIM(i2s(i)), & 320 "Electric Field[Electric Field:3]" ) 321 ELSE 322 CALL ListAddString( SolverParams, "Exported Variable "//TRIM(i2s(i)), & 323 "Electric Field[Electric Field re:3 Electric Field im:3]" ) 324 END IF 325 END IF 326 327 IF ( GetLogical( SolverParams, 'Calculate Winding Voltage', Found ) ) THEN 328 i = i + 1 329 IF ( RealField ) THEN 330 CALL ListAddString( SolverParams, "Exported Variable "//TRIM(i2s(i)), & 331 "Winding Voltage" ) 332 ELSE 333 CALL ListAddString( SolverParams, "Exported Variable "//TRIM(i2s(i)), & 334 "Winding Voltage[Winding Voltage re:1 Winding Voltage im:1]" ) 335 END IF 336 END IF 337 END IF 338 339 IF ( GetLogical( SolverParams, 'Calculate Nodal Heating', Found ) ) THEN 340 i = i + 1 341 CALL ListAddString( SolverParams, "Exported Variable "//TRIM(i2s(i)), & 342 "Nodal Joule Heating" ) 343 END IF 344 345 IF (GetLogical(SolverParams, 'Calculate Nodal Forces', Found) ) THEN 346 IF( RealField ) THEN 347 i = i + 1 348 CALL ListAddString( SolverParams, "Exported Variable "//TRIM(i2s(i)), & 349 "Nodal Force[Nodal Force:3]" ) 350 ELSE 351 i = i + 1 352 CALL ListAddString( SolverParams, "Exported Variable "//TRIM(i2s(i)), & 353 "Nodal Force[Nodal Force:3]" ) 354 CALL Warn('MagnetcDynamicsCalcFields',& 355 'Calculating experimental average nodal forces. Use at own risk.') 356 END IF 357 END IF 358 END IF 359 360 ! If we have DG for the standard fields they are already elemental... 361 IF (GetLogical(SolverParams,'Discontinuous Galerkin',Found)) RETURN 362 363 ! Choose elemental if not otherwise specified. 364 ElementalFields = .NOT. GetLogical( SolverParams, 'Skip Elemental Fields', Found) 365 IF(.NOT. Found ) ElementalFields = GetLogical( SolverParams, 'Calculate Elemental Fields', Found) 366 IF(.NOT. Found ) ElementalFields = .TRUE. 367 368 IF( ElementalFields ) THEN 369 i = i + 1 370 IF ( RealField ) THEN 371 CALL ListAddString( SolverParams, "Exported Variable "//TRIM(i2s(i)), & 372 "-dg Magnetic Flux Density E[Magnetic Flux Density E:3]" ) 373 ELSE 374 CALL ListAddString( SolverParams, "Exported Variable "//TRIM(i2s(i)), & 375 "-dg Magnetic Flux Density E[Magnetic Flux Density re E:3 Magnetic Flux Density im E:3]" ) 376 END IF 377 378 IF (GetLogical(SolverParams,'Calculate Magnetic Vector Potential',Found)) THEN 379 i = i + 1 380 IF ( RealField ) THEN 381 CALL ListAddString( SolverParams, "Exported Variable "//TRIM(i2s(i)), & 382 "-dg Magnetic Vector Potential E[Magnetic Vector Potential E:3]" ) 383 ELSE 384 CALL ListAddString( SolverParams, "Exported Variable "//TRIM(i2s(i)), & 385 "-dg Magnetic Vector Potential E[Magnetic Vector Potential re E:3 Magnetic Vector Potential im E:3]" ) 386 END IF 387 END IF 388 389 IF (GetLogical(SolverParams,'Calculate Magnetic Field Strength',Found)) THEN 390 i = i + 1 391 IF ( RealField ) THEN 392 CALL ListAddString( SolverParams, "Exported Variable "//TRIM(i2s(i)), & 393 "-dg Magnetic Field Strength E[Magnetic Field Strength E:3]" ) 394 ELSE 395 CALL ListAddString( SolverParams, "Exported Variable "//TRIM(i2s(i)), & 396 "-dg Magnetic Field Strength E[Magnetic Field Strength re E:3 Magnetic Field Strength im E:3]" ) 397 END IF 398 END IF 399 400 IF (GetLogical(SolverParams,'Calculate JxB',Found)) THEN 401 i = i + 1 402 IF ( RealField ) THEN 403 CALL ListAddString( SolverParams, "Exported Variable "//TRIM(i2s(i)), & 404 "-dg JxB E[JxB E:3]" ) 405 ELSE 406 CALL ListAddString( SolverParams, "Exported Variable "//TRIM(i2s(i)), & 407 "-dg JxB E[JxB re E:3 JxB im E:3]" ) 408 END IF 409 END IF 410 411 IF ( GetLogical( SolverParams, 'Calculate Maxwell Stress', Found ) ) THEN 412 i = i + 1 413 IF ( RealField ) THEN 414 CALL ListAddString( SolverParams, "Exported Variable "//TRIM(i2s(i)), & 415 "-dg Maxwell Stress E[Maxwell Stress E:6]" ) 416 ELSE 417 CALL ListAddString( SolverParams, "Exported Variable "//TRIM(i2s(i)), & 418 "-dg Maxwell Stress E[Maxwell Stress re E:6 Maxwell Stress im E:6]" ) 419 END IF 420 END IF 421 422 IF ( GetLogical( SolverParams, 'Calculate Current Density', Found ) ) THEN 423 i = i + 1 424 IF ( RealField ) THEN 425 CALL ListAddString( SolverParams, "Exported Variable "//TRIM(i2s(i)), & 426 "-dg Current Density E[Current Density E:3]" ) 427 ELSE 428 CALL ListAddString( SolverParams, "Exported Variable "//TRIM(i2s(i)), & 429 "-dg Current Density E[Current Density re E:3 Current Density im E:3]" ) 430 END IF 431 END IF 432 433 IF ( GetLogical( SolverParams, 'Calculate Joule Heating', Found ) ) THEN 434 i = i + 1 435 CALL ListAddString( SolverParams, "Exported Variable "//TRIM(i2s(i)), & 436 "-dg Joule Heating E" ) 437 END IF 438 439 IF ( GetLogical( SolverParams, 'Calculate Harmonic Loss', Found ) ) THEN 440 IF( RealField ) THEN 441 CALL Warn('MagnetoDynamicsCalcFields',& 442 'Harmonic loss computation only available for complex systems!') 443 ELSE 444 i = i + 1 445 CALL ListAddString( SolverParams, "Exported Variable "//TRIM(i2s(i)), & 446 "-dg Harmonic Loss Linear E" ) 447 i = i + 1 448 CALL ListAddString( SolverParams, "Exported Variable "//TRIM(i2s(i)), & 449 "-dg Harmonic Loss Quadratic E" ) 450 END IF 451 END IF 452 453 IF ( Transient .OR. ComplexField .OR. LorentzConductivity ) THEN 454 IF ( GetLogical( SolverParams, 'Calculate Electric Field', Found ) ) THEN 455 i = i + 1 456 IF ( RealField ) THEN 457 CALL ListAddString( SolverParams, "Exported Variable "//TRIM(i2s(i)), & 458 "-dg Electric Field E[Electric Field E:3]" ) 459 ELSE 460 CALL ListAddString( SolverParams, "Exported Variable "//TRIM(i2s(i)), & 461 "-dg Electric Field E[Electric Field re E:3 Electric Field im E:3]" ) 462 END IF 463 END IF 464 465 IF ( GetLogical( SolverParams, 'Calculate Winding Voltage', Found ) ) THEN 466 i = i + 1 467 IF ( RealField ) THEN 468 CALL ListAddString( SolverParams, "Exported Variable "//TRIM(i2s(i)), & 469 "-dg Winding Voltage E" ) 470 ELSE 471 CALL ListAddString( SolverParams, "Exported Variable "//TRIM(i2s(i)), & 472 "-dg Winding Voltage E[Winding Voltage re E:1 Winding Voltage im E:1]" ) 473 END IF 474 END IF 475 END IF 476 477 IF (GetLogical(SolverParams, 'Calculate Nodal Forces', Found) ) THEN 478 IF( RealField ) THEN 479 i = i + 1 480 CALL ListAddString( SolverParams, "Exported Variable "//TRIM(i2s(i)), & 481 "-dg Nodal Force E[Nodal Force E:3]" ) 482 ELSE 483 i = i + 1 484 CALL ListAddString( SolverParams, "Exported Variable "//TRIM(i2s(i)), & 485 "-dg Nodal Force E[Nodal Force E:3]" ) 486 CALL Warn('MagnetcDynamicsCalcFields',& 487 'Calculating experimental average nodal forces. Use at own risk.') 488 END IF 489 END IF 490 END IF 491 492!------------------------------------------------------------------------------ 493END SUBROUTINE MagnetoDynamicsCalcFields_Init 494!------------------------------------------------------------------------------ 495 496!------------------------------------------------------------------------------ 497!> Calculate fields resulting from the edge element formulation of the magnetic 498!> field equations. 499!> \ingroup Solvers 500!------------------------------------------------------------------------------ 501 SUBROUTINE MagnetoDynamicsCalcFields(Model,Solver,dt,Transient) 502!------------------------------------------------------------------------------ 503 USE MagnetoDynamicsUtils 504 USE CircuitUtils 505 USE Zirka 506 use zirkautils 507 508 IMPLICIT NONE 509!------------------------------------------------------------------------------ 510 TYPE(Solver_t) :: Solver 511 TYPE(Model_t) :: Model 512 REAL(KIND=dp) :: dt 513 LOGICAL :: Transient 514!------------------------------------------------------------------------------ 515! The following arrays have hard-coded sizes which may need to be altered if 516! new finite elements are added. Current defaults are 54 edge finite element 517! DOFs and 27 nodal DOFs at maximum (obtained for the second-order brick over 518! a background element of type 827): 519!------------------------------------------------------------------------------ 520 REAL(KIND=dp) :: WBasis(54,3), RotWBasis(54,3), Basis(27), dBasisdx(27,3) 521 REAL(KIND=dp) :: SOL(2,81), PSOL(81), ElPotSol(1,27), C(27) 522 REAL(KIND=dp) :: Wbase(27), alpha(27), NF_ip(27,3) 523 REAL(KIND=dp) :: PR(27), omega_velo(3,27), lorentz_velo(3,27) 524 COMPLEX(KIND=dp) :: Magnetization(3,27), BodyForceCurrDens(3,27) 525 COMPLEX(KIND=dp) :: R_Z(27) 526!------------------------------------------------------------------------------ 527 REAL(KIND=dp) :: s,u,v,w, Norm 528 REAL(KIND=dp) :: B(2,3), E(2,3), JatIP(2,3), VP_ip(2,3), JXBatIP(2,3), CC_J(2,3), HdotB 529 REAL(KIND=dp) :: detJ, C_ip, PR_ip, ST(3,3), Omega, ThinLinePower, Power, Energy(3), w_dens 530 REAL(KIND=dp) :: Freq, FreqPower, FieldPower, LossCoeff, ValAtIP 531 REAL(KIND=dp) :: Freq2, FreqPower2, FieldPower2, LossCoeff2 532 REAL(KIND=dp) :: ComponentLoss(2,2), rot_velo(3), angular_velo(3) 533 REAL(KIND=dp) :: Coeff, Coeff2, TotalLoss(3), LumpedForce(3), localAlpha, localV(2), nofturns, coilthickness 534 REAL(KIND=dp) :: Flux(2), AverageFluxDensity(2), Area, N_j, wvec(3), PosCoord(3), TorqueDeprecated(3) 535 REAL(KIND=dp) :: R_ip, mu_r 536 537 COMPLEX(KIND=dp) :: MG_ip(3), BodyForceCurrDens_ip(3) 538 COMPLEX(KIND=dp) :: CST(3,3) 539 COMPLEX(KIND=dp) :: CMat_ip(3,3) 540 COMPLEX(KIND=dp) :: imag_value, Zs 541 COMPLEX(KIND=dp), ALLOCATABLE :: Tcoef(:,:,:) 542 COMPLEX(KIND=dp), POINTER, SAVE :: Reluct_Z(:,:,:) => NULL() 543 COMPLEX(KIND=dp) :: R_ip_Z, Nu(3,3) 544 545 INTEGER, PARAMETER :: ind1(6) = [1,2,3,1,2,1] 546 INTEGER, PARAMETER :: ind2(6) = [1,2,3,2,3,3] 547 548 TYPE(Variable_t), POINTER :: Var, MFD, MFS, CD, EF, MST, & 549 JH, NJH, VP, FWP, JXB, ML, ML2, LagrangeVar, NF 550 TYPE(Variable_t), POINTER :: EL_MFD, EL_MFS, EL_CD, EL_EF, & 551 EL_MST, EL_JH, EL_VP, EL_FWP, EL_JXB, EL_ML, EL_ML2, & 552 EL_NF 553 554 INTEGER :: Active,i,j,k,l,m,n,nd,np,p,q,DOFs,vDOFs,dim,BodyId,& 555 VvarDofs,VvarId,IvarId,Reindex,Imindex,EdgeBasisDegree 556 557 TYPE(Solver_t), POINTER :: pSolver, ElPotSolver 558 CHARACTER(LEN=MAX_NAME_LEN) :: Pname, CoilType, ElectricPotName, LossFile, CurrPathPotName 559 560 TYPE(ValueList_t), POINTER :: Material, BC, BodyForce, BodyParams, SolverParams 561 562 LOGICAL :: Found, FoundMagnetization, stat, Cubic, LossEstimation, & 563 CalcFluxLogical, CoilBody, PreComputedElectricPot, ImposeCircuitCurrent, & 564 ItoJCoeffFound, ImposeBodyForceCurrent, HasVelocity, HasAngularVelocity, & 565 HasLorenzVelocity, HaveAirGap, UseElementalNF, HasTensorReluctivity, & 566 ImposeBodyForcePotential, JouleHeatingFromCurrent, HasZirka 567 LOGICAL :: PiolaVersion, ElementalFields, NodalFields, RealField, SecondOrder 568 LOGICAL :: CSymmetry, HBCurve, LorentzConductivity, HasThinLines=.FALSE. 569 570 TYPE(GaussIntegrationPoints_t) :: IP 571 TYPE(Nodes_t), SAVE :: Nodes 572 TYPE(Element_t), POINTER :: Element 573 574 INTEGER, ALLOCATABLE :: Pivot(:), TorqueGroups(:) 575 INTEGER, POINTER :: MasterBodies(:) 576 577 REAL(KIND=dp), POINTER CONTIG :: Fsave(:) 578 REAL(KIND=dp), POINTER :: HB(:,:)=>NULL(), CubicCoeff(:)=>NULL(), & 579 HBBVal(:), HBCval(:), HBHval(:) 580 REAL(KIND=dp) :: Babs 581 TYPE(Mesh_t), POINTER :: Mesh 582 REAL(KIND=dp), ALLOCATABLE, TARGET :: Gforce(:,:), MASS(:,:), FORCE(:,:) 583 REAL(KIND=dp), ALLOCATABLE :: BodyLoss(:,:), RotM(:,:,:), Torque(:) 584 585 REAL(KIND=dp), ALLOCATABLE :: ThinLineCrossect(:),ThinLineCond(:) 586 587 REAL(KIND=DP), POINTER :: Cwrk(:,:,:)=>NULL(), Cwrk_im(:,:,:)=>NULL() 588 589 REAL(KIND=dp) :: ItoJCoeff, CircuitCurrent, CircEqVoltageFactor 590 TYPE(ValueList_t), POINTER :: CompParams 591 REAL(KIND=dp) :: DetF, F(3,3), G(3,3), GT(3,3) 592 REAL(KIND=dp), ALLOCATABLE :: EBasis(:,:), CurlEBasis(:,:) 593 594 REAL(KIND=dp) :: xcoord, grads_coeff, val 595 TYPE(ValueListEntry_t), POINTER :: HBLst 596 REAL(KIND=dp) :: HarmPowerCoeff 597 REAL(KIND=dp) :: line_tangent(3) 598 INTEGER :: IOUnit, pIndex 599 REAL(KIND=dp) :: SaveNorm 600 INTEGER :: NormIndex 601 602 INTEGER, POINTER, SAVE :: SetPerm(:) => NULL() 603!------------------------------------------------------------------------------------------- 604 IF ( .NOT. ASSOCIATED( Solver % Matrix ) ) RETURN 605 606 CALL Info('MagnetoDynamicsCalcFields','------------------------------',Level=6) 607 CALL Info('MagnetoDynamicsCalcFields','Computing postprocessed fields',Level=5) 608 609 dim = CoordinateSystemDimension() 610 SolverParams => GetSolverParams() 611 612 ! This is a hack to be able to control the norm that is tested for 613 NormIndex = GetInteger(SolverParams,'Show Norm Index',Found ) 614 SaveNorm = 0.0_dp 615 616 IF (GetLogical(SolverParams, 'Calculate harmonic peak power', Found)) THEN 617 HarmPowerCoeff = 1.0_dp 618 ELSE 619 HarmPowerCoeff = 0.5_dp 620 END IF 621 622 pIndex = ListGetInteger( SolverParams,'Target Variable Solver Index',UnfoundFatal=.TRUE.) 623 pSolver => Model % Solvers(pIndex) 624 pname = getVarName(pSolver % Variable) 625 626 CALL Info('MagnetoDynamicsCalcFields','Using potential variable: '//TRIM(Pname),Level=7) 627 628 ! Inherit the solution basis from the primary solver 629 vDOFs = pSolver % Variable % DOFs 630 SecondOrder = GetLogical( pSolver % Values, 'Quadratic Approximation', Found ) 631 IF (SecondOrder) THEN 632 EdgeBasisDegree = 2 633 ELSE 634 EdgeBasisDegree = 1 635 END IF 636 637 IF( SecondOrder ) THEN 638 PiolaVersion = .TRUE. 639 ELSE 640 PiolaVersion = GetLogical( pSolver % Values,'Use Piola Transform', Found ) 641 END IF 642 643 IF (PiolaVersion) & 644 CALL Info('MagnetoDynamicsCalcFields', & 645 'Using Piola transformed finite elements',Level=5) 646 647 ElectricPotName = GetString(SolverParams, 'Precomputed Electric Potential', PrecomputedElectricPot) 648 IF (PrecomputedElectricPot) THEN 649 DO i=1, Model % NumberOfSolvers 650 ElPotSolver => Model % Solvers(i) 651 IF (ElectricPotName==getVarName(ElPotSolver % Variable)) EXIT 652 END DO 653 END IF 654 655 ! Do we have a real or complex valued primary field? 656 RealField = ( vDofs == 1 ) 657 658 LorentzConductivity = ListCheckPrefixAnyBodyForce(Model, "Angular Velocity") .or. & 659 ListCheckPrefixAnyBodyForce(Model, "Lorentz Velocity") 660 661 Mesh => GetMesh() 662 LagrangeVar => VariableGet( Solver % Mesh % Variables,'LagrangeMultiplier', ThisOnly=.TRUE.) 663 664 MFD => VariableGet( Mesh % Variables, 'Magnetic Flux Density' ) 665 EL_MFD => VariableGet( Mesh % Variables, 'Magnetic Flux Density E' ) 666 667 MFS => VariableGet( Mesh % Variables, 'Magnetic Field Strength') 668 EL_MFS => VariableGet( Mesh % Variables, 'Magnetic Field Strength E') 669 670 VP => VariableGet( Mesh % Variables, 'Magnetic Vector Potential') 671 EL_VP => VariableGet( Mesh % Variables, 'Magnetic Vector Potential E') 672 673 IF( .NOT. PreComputedElectricPot ) THEN 674 ImposeBodyForcePotential = GetLogical(SolverParams, 'Impose Body Force Potential', Found) 675 IF (.NOT. Found) ImposeBodyForcePotential = & 676 ListCheckPresentAnyBodyForce( Model,'Electric Potential') 677 ELSE 678 ImposeBodyForcePotential = .FALSE. 679 END IF 680 681 ImposeBodyForceCurrent = GetLogical(SolverParams, 'Impose Body Force Current', Found) 682 IF (.NOT. Found) ImposeBodyForceCurrent = ListCheckPrefixAnyBodyForce( Model,'Current Density') 683 684 ImposeCircuitCurrent = GetLogical(SolverParams, 'Impose Circuit Current', Found) 685 CurrPathPotName = GetString(SolverParams, 'Circuit Current Path Potential Name', Found) 686 IF (.NOT. Found) CurrPathPotName = 'W' 687 688 EF => NULL(); EL_EF => NULL(); 689 CD => NULL(); EL_CD => NULL(); 690 JH => NULL(); EL_JH => NULL(); 691 FWP => NULL(); EL_FWP => NULL(); 692 JXB => NULL(); EL_JXB => NULL(); 693 ML => NULL(); EL_ML => NULL(); 694 ML2 => NULL(); EL_ML2 => NULL(); 695 NF => NULL(); EL_NF => NULL(); 696 NJH => NULL() 697 698 IF ( Transient .OR. .NOT. RealField .OR. LorentzConductivity ) THEN 699 EF => VariableGet( Mesh % Variables, 'Electric Field' ) 700 FWP => VariableGet( Mesh % Variables, 'Winding Voltage' ) 701 702 EL_EF => VariableGet( Mesh % Variables, 'Electric Field E' ) 703 EL_FWP => VariableGet( Mesh % Variables, 'Winding Voltage E' ) 704 END IF 705 706 NF => VariableGet( Mesh % Variables, 'Nodal Force') 707 EL_NF => VariableGet( Mesh % Variables, 'Nodal Force E') 708 709 CD => VariableGet( Mesh % Variables, 'Current Density' ) 710 EL_CD => VariableGet( Mesh % Variables, 'Current Density E' ) 711 712 JH => VariableGet( Mesh % Variables, 'Joule Heating' ) 713 EL_JH => VariableGet( Mesh % Variables, 'Joule Heating E' ) 714 715 NJH => VariableGet( Mesh % Variables, 'Nodal Joule Heating' ) 716 717 IF(.NOT. RealField ) THEN 718 ML => VariableGet( Mesh % Variables, 'Harmonic Loss Linear') 719 EL_ML => VariableGet( Mesh % Variables, 'Harmonic Loss Linear E') 720 ML2 => VariableGet( Mesh % Variables, 'Harmonic Loss Quadratic') 721 EL_ML2 => VariableGet( Mesh % Variables, 'Harmonic Loss Quadratic E') 722 END IF 723 724 JXB => VariableGet( Mesh % Variables, 'JxB') 725 EL_JXB => VariableGet( Mesh % Variables, 'JxB E') 726 727 MST => variableGet( Mesh % Variables, 'Maxwell stress' ) 728 EL_MST => variableGet( Mesh % Variables, 'Maxwell stress E' ) 729 730 DOFs = 0 731 IF ( ASSOCIATED(MFD) ) DOFs=DOFs+3 732 IF ( ASSOCIATED(MFS) ) DOFs=DOFs+3 733 IF ( ASSOCIATED(VP) ) DOFs=DOFs+3 734 IF ( ASSOCIATED(CD) ) DOFs=DOFs+3 735 IF ( ASSOCIATED(FWP) ) DOFs=DOFs+1 736 IF ( ASSOCIATED(EF) ) DOFs=DOFs+3 737 IF ( ASSOCIATED(JXB) ) DOFs=DOFs+3 738 IF ( ASSOCIATED(MST) ) DOFs=DOFs+6 739 IF ( ASSOCIATED(NF) ) DOFs=DOFs+3 740 DOFs = DOFs*vDOFs 741 IF ( ASSOCIATED(JH) .OR. ASSOCIATED(NJH)) DOFs=DOFs+1 742 IF ( ASSOCIATED(ML) ) DOFs=DOFs+1 743 IF ( ASSOCIATED(ML2) ) DOFs=DOFs+1 744 NodalFields = DOFs > 0 745 746 IF(NodalFields) THEN 747 ALLOCATE(GForce(SIZE(Solver % Matrix % RHS),DOFs)); Gforce=0._dp 748 ELSE 749 DOFs = 0 750 IF ( ASSOCIATED(EL_MFD) ) DOFs=DOFs+3 751 IF ( ASSOCIATED(EL_MFS) ) DOFs=DOFs+3 752 IF ( ASSOCIATED(EL_VP) ) DOFs=DOFs+3 753 IF ( ASSOCIATED(EL_CD) ) DOFs=DOFs+3 754 IF ( ASSOCIATED(EL_FWP) ) DOFs=DOFs+1 755 IF ( ASSOCIATED(EL_EF) ) DOFs=DOFs+3 756 IF ( ASSOCIATED(EL_JXB) ) DOFs=DOFs+3 757 IF ( ASSOCIATED(EL_MST) ) DOFs=DOFs+6 758 IF ( ASSOCIATED(EL_NF) ) DOFs=DOFs+3 759 DOFs = DOFs*vDOFs 760 IF ( ASSOCIATED(EL_NF) ) DOFs=DOFs+3 761 IF ( ASSOCIATED(EL_JH) ) DOFs=DOFs+1 762 IF ( ASSOCIATED(EL_ML) ) DOFs=DOFs+1 763 IF ( ASSOCIATED(EL_ML2) ) DOFs=DOFs+1 764 END IF 765 766 ElementalFields = .FALSE. 767 IF ( ASSOCIATED(EL_MFD) ) ElementalFields=.TRUE. 768 IF ( ASSOCIATED(EL_MFS) ) ElementalFields=.TRUE. 769 IF ( ASSOCIATED(EL_VP) ) ElementalFields=.TRUE. 770 IF ( ASSOCIATED(EL_CD) ) ElementalFields=.TRUE. 771 IF ( ASSOCIATED(EL_FWP) ) ElementalFields=.TRUE. 772 IF ( ASSOCIATED(EL_EF) ) ElementalFields=.TRUE. 773 IF ( ASSOCIATED(EL_JXB) ) ElementalFields=.TRUE. 774 IF ( ASSOCIATED(EL_MST) ) ElementalFields=.TRUE. 775 IF ( ASSOCIATED(EL_NF) ) ElementalFields=.TRUE. 776 IF ( ASSOCIATED(EL_JH) ) ElementalFields=.TRUE. 777 IF ( ASSOCIATED(EL_ML) ) ElementalFields=.TRUE. 778 IF ( ASSOCIATED(EL_ML2) ) ElementalFields=.TRUE. 779 780 n = Mesh % MaxElementDOFs 781 ALLOCATE( MASS(n,n), FORCE(n,DOFs), Tcoef(3,3,n), RotM(3,3,n), Pivot(n)) 782 783 SOL = 0._dp; PSOL=0._dp 784 785 LossEstimation = GetLogical(SolverParams,'Loss Estimation',Found) & 786 .OR. ASSOCIATED( ML ) .OR. ASSOCIATED( EL_ML ) & 787 .OR. ASSOCIATED( ML2 ) .OR. ASSOCIATED( EL_ML2 ) 788 789 IF (LossEstimation) THEN 790 FreqPower = GetCReal( SolverParams,'Harmonic Loss Linear Frequency Exponent',Found ) 791 IF( .NOT. Found ) FreqPower = 1.0_dp 792 793 FreqPower2 = GetCReal( SolverParams,'Harmonic Loss Quadratic Frequency Exponent',Found ) 794 IF( .NOT. Found ) FreqPower2 = 2.0_dp 795 796 FieldPower = GetCReal( SolverParams,'Harmonic Loss Linear Exponent',Found ) 797 IF( .NOT. Found ) FieldPower = 2.0_dp 798 FieldPower = FieldPower / 2.0_dp 799 800 FieldPower2 = GetCReal( SolverParams,'Harmonic Loss Quadratic Exponent',Found ) 801 IF( .NOT. Found ) FieldPower2 = 2.0_dp 802 FieldPower2 = FieldPower2 / 2.0_dp 803 804 IF(.NOT. ListCheckPresentAnyMaterial( Model,'Harmonic Loss Linear Coefficient') ) THEN 805 CALL Warn('MagnetoDynamicsCalcFields',& 806 'Harmonic loss requires > Harmonic Loss Linear Coefficient < in material section!') 807 END IF 808 809 IF(.NOT. ListCheckPresentAnyMaterial( Model,'Harmonic Loss Quadratic Coefficient') ) THEN 810 CALL Warn('MagnetoDynamicsCalcFields',& 811 'Harmonic loss requires > Harmonic Loss Quadratic Coefficient < in material section!') 812 END IF 813 814 ComponentLoss = 0.0_dp 815 ALLOCATE( BodyLoss(3,Model % NumberOfBodies) ) 816 BodyLoss = 0.0_dp 817 END IF 818 819 820 C = 0._dp; PR=0._dp 821 Magnetization = 0._dp 822 823 Power = 0._dp; Energy = 0._dp 824 CALL DefaultInitialize() 825 826 827 DO i = 1, GetNOFActive() 828 Element => GetActiveElement(i) 829 n = GetElementNOFNodes() 830 np = n*pSolver % Def_Dofs(GetElementFamily(Element),Element % BodyId,1) 831 nd = GetElementNOFDOFs(uSolver=pSolver) 832 833 IF (SIZE(Tcoef,3) /= n) THEN 834 DEALLOCATE(Tcoef) 835 ALLOCATE(Tcoef(3,3,n)) 836 END IF 837 838 CALL GetElementNodes( Nodes ) 839 840 ! If potential is not available we have to use given current directly to estimate Joule losses 841 JouleHeatingFromCurrent = ( np == 0 .AND. & 842 .NOT. ( PreComputedElectricPot .OR. ImposeBodyForcePotential ) ) 843 844 BodyId = GetBody() 845 Material => GetMaterial() 846 BodyForce => GetBodyForce() 847 848 ItoJCoeffFound = .FALSE. 849 IF(ImposeCircuitCurrent) THEN 850 ItoJCoeff = ListGetConstReal(GetBodyParams(Element), & 851 'Current to density coefficient', ItoJCoeffFound) 852 IF(ItoJCoeffFound) THEN 853 CALL GetLocalSolution(Wbase,CurrPathPotName) 854 IvarId = GetInteger(GetBodyParams(Element), 'Circuit Current Variable Id', Found) 855 IF (.NOT. Found) CALL Fatal ('MagnetoDynamicsCalcFields', 'Circuit Current Variable Id not found!') 856 CircuitCurrent = LagrangeVar % Values(IvarId) 857 END IF 858 END IF 859 860 861 CALL GetVectorLocalSolution(SOL,Pname,uSolver=pSolver) 862 IF (PrecomputedElectricPot) THEN 863 CALL GetScalarLocalSolution(ElPotSol(1,:),ElectricPotName,uSolver=ElPotSolver) 864 END IF 865 866 IF( ImposeBodyForcePotential ) THEN 867 ElPotSol(1,:) = GetReal(BodyForce,'Electric Potential',Found) 868 END IF 869 870 IF ( Transient ) THEN 871 CALL GetScalarLocalSolution(PSOL,Pname,uSolver=pSolver,Tstep=-1) 872 PSOL(1:nd)=(SOL(1,1:nd)-PSOL(1:nd))/dt 873 END IF 874 875 Omega = GetAngularFrequency(pSOlver % Values,Found,Element) 876 IF( .NOT. ( RealField .OR. Found ) ) THEN 877 CALL Fatal('MagnetoDynamicsCalcFields',& 878 '(Angular) Frequency must be given for complex fields!') 879 END IF 880 Freq = Omega / (2*PI) 881 882 IF ( ASSOCIATED(MFS) ) THEN 883 FoundMagnetization = .FALSE. 884 IF(ASSOCIATED(BodyForce)) THEN 885 CALL GetComplexVector( BodyForce,Magnetization(1:3,1:n),'Magnetization',FoundMagnetization) 886 END IF 887 888 IF(.NOT.FoundMagnetization) THEN 889 CALL GetComplexVector( BodyForce,Magnetization(1:3,1:n),'Magnetization',FoundMagnetization) 890 END IF 891 END IF 892 893 IF (ImposeBodyForceCurrent ) THEN 894 BodyForceCurrDens = 0._dp 895 IF ( ASSOCIATED(BodyForce) ) THEN 896 SELECT CASE(DIM) 897 CASE(3) 898 CALL GetComplexVector( BodyForce, BodyForceCurrDens(1:3,1:n), 'Current Density', Found ) 899 CASE(2) 900 IF (Vdofs == 2) THEN 901 BodyForceCurrDens(3,1:n) = CMPLX( ListGetReal(BodyForce, 'Current Density', n, Element % NodeIndexes, Found), & 902 ListGetReal(BodyForce, 'Current Density im', n, Element % NodeIndexes, Found), KIND=dp) 903 ELSE 904 BodyForceCurrDens(3,1:n) = CMPLX( ListGetReal(BodyForce, 'Current Density', & 905 n, Element % NodeIndexes, Found), 0, KIND=dp) 906 END IF 907 END SELECT 908 909 END IF 910 END IF 911 912 CALL GetPermittivity(Material,PR,n) 913 914 CoilBody = .FALSE. 915 CompParams => GetComponentParams( Element ) 916 CoilType = '' 917 RotM = 0._dp 918 IF (ASSOCIATED(CompParams)) THEN 919 CoilType = GetString(CompParams, 'Coil Type', Found) 920 IF (Found) CoilBody = .TRUE. 921 CircEqVoltageFactor = GetConstReal(CompParams, 'Circuit Equation Voltage Factor', Found) 922 IF (.NOT. Found) CircEqVoltageFactor = 1._dp 923 END IF 924 925 !------------------------------------------------------------------------------ 926 ! Read conductivity values (might be a tensor) 927 !------------------------------------------------------------------------------ 928 !C(1:n) = GetReal(Material,'Electric Conductivity',Found) 929 Tcoef = GetCMPLXElectricConductivityTensor(Element, n, CoilBody, CoilType) 930 931 dim = CoordinateSystemDimension() 932 CSymmetry = ( CurrentCoordinateSystem() == AxisSymmetric .OR. & 933 CurrentCoordinateSystem() == CylindricSymmetric ) 934 935 IF (CoilBody) THEN 936 937 Call GetWPotential(Wbase) 938 939 SELECT CASE (CoilType) 940 CASE ('stranded') 941 IvarId = GetInteger (CompParams, 'Circuit Current Variable Id', Found) 942 IF (.NOT. Found) CALL Fatal ('MagnetoDynamicsCalcFields', 'Circuit Current Variable Id not found!') 943 944 N_j = GetConstReal (CompParams, 'Stranded Coil N_j', Found) 945 IF (.NOT. Found) CALL Fatal ('MagnetoDynamicsCalcFields', 'Stranded Coil N_j not found!') 946 947 nofturns = GetConstReal(CompParams, 'Number of Turns', Found) 948 IF (.NOT. Found) CALL Fatal('MagnetoDynamicsCalcFields','Stranded Coil: Number of Turns not found!') 949 CASE ('massive') 950 VvarId = GetInteger (CompParams, 'Circuit Voltage Variable Id', Found) 951 IF (.NOT. Found) CALL Fatal ('MagnetoDynamicsCalcFields', 'Circuit Voltage Variable Id not found!') 952 953 CASE ('foil winding') 954 CALL GetLocalSolution(alpha,'Alpha') 955 956 IF (dim == 3) CALL GetElementRotM(Element, RotM, n) 957 958 VvarId = GetInteger (CompParams, 'Circuit Voltage Variable Id', Found) 959 IF (.NOT. Found) CALL Fatal ('MagnetoDynamicsCalcFields', 'Circuit Voltage Variable Id not found!') 960 961 coilthickness = GetConstReal(CompParams, 'Coil Thickness', Found) 962 IF (.NOT. Found) CALL Fatal('MagnetoDynamicsCalcFields','Foil Winding: Coil Thickness not found!') 963 964 nofturns = GetConstReal(CompParams, 'Number of Turns', Found) 965 IF (.NOT. Found) CALL Fatal('MagnetoDynamicsCalcFields','Foil Winding: Number of Turns not found!') 966 967 VvarDofs = GetInteger (CompParams, 'Circuit Voltage Variable dofs', Found) 968 IF (.NOT. Found) CALL Fatal ('MagnetoDynamicsCalcFields', 'Circuit Voltage Variable dofs not found!') 969 ! in case of a foil winding, transform the conductivity tensor: 970 ! ------------------------------------------------------------- 971 972 IF (dim == 3) THEN 973 DO k = 1,n 974 Tcoef(1:3,1:3,k) = MATMUL(MATMUL(RotM(1:3,1:3,k), Tcoef(1:3,1:3,k)), TRANSPOSE(RotM(1:3,1:3,k))) 975 END DO 976 END IF 977 CASE DEFAULT 978 CALL Fatal ('MagnetoDynamicsCalcFields', 'Non existent Coil Type Chosen!') 979 END SELECT 980 END IF 981 982 983 !--------------------------------------------------------------------------------------------- 984 R_Z = CMPLX(0.0_dp, 0.0_dp, kind=dp) 985 HasTensorReluctivity = .FALSE. 986 CALL GetConstRealArray( Material, HB, 'H-B curve', Found ) 987 IF ( ASSOCIATED(HB) ) THEN 988 Cubic = GetLogical( Material, 'Cubic spline for H-B curve', Found) 989 l = SIZE(HB,1) 990 HBBval => HB(:,1) 991 HBHval => HB(:,2) 992 IF(l>1) THEN 993 IF (Cubic.AND..NOT.ASSOCIATED(CubicCoeff) ) THEN 994 ALLOCATE(CubicCoeff(l)) 995 CALL CubicSpline(l,HB(:,1),HB(:,2),CubicCoeff) 996 END IF 997 HBCval => CubicCoeff 998 END IF 999 1000 IF(l<=1) THEN 1001 HBLst => ListFind(Material,'H-B Curve',HBcurve) 1002 IF(HBcurve) THEN 1003 HBCval => HBLst % CubicCoeff 1004 HBBval => HBLst % TValues 1005 HBHval => HBLst % FValues(1,1,:) 1006 END IF 1007 END IF 1008 ELSE 1009 ! 1010 ! Seek reluctivity as complex-valued: A given reluctivity can be a tensor 1011 ! 1012 CALL GetReluctivity(Material,Reluct_Z,n,HasTensorReluctivity) 1013 IF (HasTensorReluctivity) THEN 1014 IF (SIZE(Reluct_Z,1)==1 .AND. SIZE(Reluct_Z,2)==1) THEN 1015 l = MIN(SIZE(R_Z), SIZE(Reluct_Z,3)) 1016 R_Z(1:l) = Reluct_Z(1,1,1:l) 1017 HasTensorReluctivity = .FALSE. 1018 ELSE 1019 R_Z = CMPLX(0.0_dp, 0.0_dp, kind=dp) 1020 END IF 1021 ELSE 1022 ! Seek via a given permeability: In this case the reluctivity will be 1023 ! a complex scalar: 1024 CALL GetReluctivity(Material,R_Z,n) 1025 END IF 1026 END IF 1027 1028 HasVelocity = .FALSE. 1029 HasLorenzVelocity = .FALSE. 1030 HasAngularVelocity = .FALSE. 1031 IF(ASSOCIATED(BodyForce)) THEN 1032 CALL GetRealVector( BodyForce, omega_velo, 'Angular velocity', HasAngularVelocity) 1033 CALL GetRealVector( BodyForce, lorentz_velo, 'Lorentz velocity', HasLorenzVelocity) 1034 HasVelocity = HasAngularVelocity .OR. HasLorenzVelocity 1035 END IF 1036 1037 1038 ! Calculate nodal fields: 1039 ! ----------------------- 1040 IF (SecondOrder) THEN 1041 IP = GaussPoints(Element, EdgeBasis=dim==3, PReferenceElement=PiolaVersion, EdgeBasisDegree=EdgeBasisDegree) 1042 ELSE 1043 IP = GaussPoints(Element, EdgeBasis=dim==3, PReferenceElement=PiolaVersion) 1044 END IF 1045 1046 MASS = 0._dp 1047 FORCE = 0._dp 1048 E = 0._dp; B=0._dp 1049 1050 haszirka = .false. 1051 if(ASSOCIATED(MFS) .OR. ASSOCIATED(el_MFS)) THEN 1052 CALL GetHystereticMFS(Element, force(:,4:6), pSolver, HasZirka, CSymmetry=CSymmetry) 1053 end if 1054 1055 DO j = 1,IP % n 1056 u = IP % U(j) 1057 v = IP % V(j) 1058 w = IP % W(j) 1059 1060 IF (PiolaVersion) THEN 1061 stat = EdgeElementInfo( Element, Nodes, u, v, w, DetF=DetJ, Basis=Basis, & 1062 EdgeBasis=WBasis, RotBasis=RotWBasis, dBasisdx=dBasisdx, & 1063 BasisDegree = EdgeBasisDegree, ApplyPiolaTransform = .TRUE.) 1064 ELSE 1065 stat=ElementInfo(Element,Nodes,u,v,w,detJ,Basis,dBasisdx) 1066 IF( dim == 3 ) THEN 1067 CALL GetEdgeBasis(Element,WBasis,RotWBasis,Basis,dBasisdx) 1068 END IF 1069 END IF 1070 1071 s = IP % s(j) * detJ 1072 1073 grads_coeff = -1._dp/GetCircuitModelDepth() 1074 IF( CSymmetry ) THEN 1075 xcoord = SUM( Basis(1:n) * Nodes % x(1:n) ) 1076 grads_coeff = grads_coeff/xcoord 1077 s = s * xcoord 1078 END IF 1079 1080 1081 DO k=1,vDOFs 1082 SELECT CASE(dim) 1083 CASE(2) 1084 ! This has been done with the same sign convention as in MagnetoDynamics2D: 1085 ! ------------------------------------------------------------------------- 1086 IF ( CSymmetry ) THEN 1087 B(k,1) = -SUM( SOL(k,1:nd) * dBasisdx(1:nd,2) ) 1088 B(k,2) = SUM( SOL(k,1:nd) * dBasisdx(1:nd,1) ) & 1089 + SUM( SOL(k,1:nd) * Basis(1:nd) ) / xcoord 1090 B(k,3) = 0._dp 1091 ELSE 1092 B(k,1) = SUM( SOL(k,1:nd) * dBasisdx(1:nd,2) ) 1093 B(k,2) = -SUM( SOL(k,1:nd) * dBasisdx(1:nd,1) ) 1094 B(k,3) = 0._dp 1095 END IF 1096 CASE(3) 1097 B(k,:) = MATMUL( SOL(k,np+1:nd), RotWBasis(1:nd-np,:) ) 1098 END SELECT 1099 END DO 1100 IF(ImposeCircuitCurrent .and. ItoJCoeffFound) THEN 1101 wvec = -MATMUL(Wbase(1:n), dBasisdx(1:n,:)) 1102 IF(SUM(wvec**2._dp) .GE. AEPS) THEN 1103 wvec = wvec/SQRT(SUM(wvec**2._dp)) 1104 ELSE 1105 wvec = [0.0_dp, 0.0_dp, 1.0_dp] 1106 END IF 1107 END IF 1108 1109 ! Compute the velocity field in the form v + w x r: 1110 ! ------------------------------------------------- 1111 IF(HasVelocity) THEN 1112 rot_velo = 0.0_dp 1113 angular_velo = 0.0_dp 1114 IF( HasAngularVelocity ) THEN 1115 angular_velo(1) = SUM(basis(1:n)*omega_velo(1,1:n)) 1116 angular_velo(2) = SUM(basis(1:n)*omega_velo(2,1:n)) 1117 angular_velo(3) = SUM(basis(1:n)*omega_velo(3,1:n)) 1118 DO k=1,n 1119 rot_velo(1:3) = rot_velo(1:3) + CrossProduct(omega_velo(1:3,k), [ & 1120 basis(k) * Nodes % x(k), & 1121 basis(k) * Nodes % y(k), & 1122 basis(k) * Nodes % z(k)]) 1123 END DO 1124 END IF 1125 IF( HasLorenzVelocity ) THEN 1126 rot_velo(1:3) = rot_velo(1:3) + [ & 1127 SUM(basis(1:n)*lorentz_velo(1,1:n)), & 1128 SUM(basis(1:n)*lorentz_velo(2,1:n)), & 1129 SUM(basis(1:n)*lorentz_velo(3,1:n))] 1130 END IF 1131 END IF 1132 !------------------------------- 1133 ! The conductivity as a tensor 1134 ! ------------------------------- 1135 !C_ip = SUM( Basis(1:n)*C(1:n) ) 1136 DO k=1,3 1137 DO l=1,3 1138 CMat_ip(k,l) = SUM( Tcoef(k,l,1:n) * Basis(1:n) ) 1139 END DO 1140 END DO 1141 BodyForceCurrDens_ip(1:3) = 0._dp 1142 IF(ImposeBodyForceCurrent) THEN 1143 DO l=1,3 1144 BodyForceCurrDens_ip(l) = SUM(BodyForceCurrDens(l,1:n)*Basis(1:n)) 1145 END DO 1146 END IF 1147 1148 IF (vDOFs > 1) THEN ! Complex case 1149 IF (CoilType /= 'stranded') THEN 1150 ! -j * Omega A 1151 SELECT CASE(dim) 1152 CASE(2) 1153 E(1,:) = 0._dp 1154 E(2,:) = 0._dp 1155 E(1,3) = Omega*SUM(SOL(2,1:nd) * Basis(1:nd)) 1156 E(2,3) = -Omega*SUM(SOL(1,1:nd) * Basis(1:nd)) 1157 CASE(3) 1158 E(1,:) = Omega*MATMUL(SOL(2,np+1:nd),WBasis(1:nd-np,:)) 1159 E(2,:) = -Omega*MATMUL(SOL(1,np+1:nd),WBasis(1:nd-np,:)) 1160 END SELECT 1161 ELSE 1162 E(1,:) = 0._dp 1163 E(2,:) = 0._dp 1164 END IF 1165 1166 localV=0._dp 1167 SELECT CASE (CoilType) 1168 1169 CASE ('stranded') 1170 SELECT CASE(dim) 1171 CASE(2) 1172 wvec = [0._dp, 0._dp, 1._dp] 1173 CASE(3) 1174 wvec = -MATMUL(Wbase(1:np), dBasisdx(1:np,:)) 1175 wvec = wvec/SQRT(SUM(wvec**2._dp)) 1176 END SELECT 1177 IF(CMat_ip(3,3) /= 0._dp ) THEN 1178 imag_value = LagrangeVar % Values(IvarId) + im * LagrangeVar % Values(IvarId+1) 1179 E(1,:) = E(1,:)+REAL(imag_value * N_j * wvec / CMat_ip(3,3)) 1180 E(2,:) = E(2,:)+AIMAG(imag_value * N_j * wvec / CMat_ip(3,3)) 1181 END IF 1182 1183 CASE ('massive') 1184 localV(1) = localV(1) + LagrangeVar % Values(VvarId) * CircEqVoltageFactor 1185 localV(2) = localV(2) + LagrangeVar % Values(VvarId+1) * CircEqVoltageFactor 1186 SELECT CASE(dim) 1187 CASE(2) 1188 E(1,3) = E(1,3)-localV(1) * grads_coeff 1189 E(2,3) = E(2,3)-localV(2) * grads_coeff 1190 CASE(3) 1191 E(1,:) = E(1,:)-localV(1) * MATMUL(Wbase(1:np), dBasisdx(1:np,:)) 1192 E(2,:) = E(2,:)-localV(2) * MATMUL(Wbase(1:np), dBasisdx(1:np,:)) 1193 END SELECT 1194 1195 CASE ('foil winding') 1196 localAlpha = coilthickness *SUM(alpha(1:np) * Basis(1:np)) 1197 DO k = 1, VvarDofs-1 1198 Reindex = 2*k 1199 Imindex = Reindex+1 1200 localV(1) = localV(1) + LagrangeVar % Values(VvarId+Reindex) * localAlpha**(k-1) * CircEqVoltageFactor 1201 localV(2) = localV(2) + LagrangeVar % Values(VvarId+Imindex) * localAlpha**(k-1) * CircEqVoltageFactor 1202 END DO 1203 SELECT CASE(dim) 1204 CASE(2) 1205 E(1,3) = E(1,3)-localV(1) * grads_coeff 1206 E(2,3) = E(2,3)-localV(2) * grads_coeff 1207 CASE(3) 1208 E(1,:) = E(1,:)-localV(1) * MATMUL(Wbase(1:np), dBasisdx(1:np,:)) 1209 E(2,:) = E(2,:)-localV(2) * MATMUL(Wbase(1:np), dBasisdx(1:np,:)) 1210 END SELECT 1211 1212 CASE DEFAULT 1213 SELECT CASE(dim) 1214 CASE(2) 1215 IF (HasLorenzVelocity) THEN 1216 ! Add v x curl A: 1217 IF (CSymmetry) THEN 1218 E(1,3) = E(1,3) - rot_velo(1) * B(1,2) + rot_velo(2) * B(1,1) 1219 E(2,3) = E(2,3) - rot_velo(1) * B(2,2) + rot_velo(2) * B(2,1) 1220 ELSE 1221 E(1,3) = E(1,3) + rot_velo(1) * B(1,2) - rot_velo(2) * B(1,1) 1222 E(2,3) = E(2,3) + rot_velo(1) * B(2,2) - rot_velo(2) * B(2,1) 1223 END IF 1224 END IF 1225 ! 1226 ! To make this perfect, the electric field corresponding to the source 1227 ! should be returned on the source region 1228 ! 1229 CASE(3) 1230 ! -Grad(V) 1231 E(1,:) = E(1,:)-MATMUL(SOL(1,1:np), dBasisdx(1:np,:)) 1232 E(2,:) = E(2,:)-MATMUL(SOL(2,1:np), dBasisdx(1:np,:)) 1233 1234 IF (HasVelocity) THEN 1235 ! 1236 ! Add v x curl A so as to handle the steady amplitude solution of 1237 ! the time harmonic equations. Multiplication with the electric 1238 ! conductivity will give the current density with respect to 1239 ! the fixed frame: 1240 ! 1241 E(1,:) = E(1,:) + CrossProduct(rot_velo, & 1242 MATMUL(SOL(1,np+1:nd), RotWBasis(1:nd-np,:))) 1243 E(2,:) = E(2,:) + CrossProduct(rot_velo, & 1244 MATMUL(SOL(2,np+1:nd), RotWBasis(1:nd-np,:))) 1245 END IF 1246 1247 IF( ImposeBodyForcePotential ) THEN 1248 E(1,:) = E(1,:) - MATMUL(ElPotSol(1,1:n), dBasisdx(1:n,:)) 1249 END IF 1250 END SELECT 1251 1252 END SELECT 1253 1254 ELSE ! Real case 1255 IF (CoilType /= 'stranded') THEN 1256 SELECT CASE(dim) 1257 CASE(2) 1258 E(1,1) = 0._dp 1259 E(1,2) = 0._dp 1260 E(1,3) = -SUM(PSOL(1:nd) * Basis(1:nd)) 1261 CASE(3) 1262 E(1,:) = -MATMUL(PSOL(np+1:nd), Wbasis(1:nd-np,:)) 1263 END SELECT 1264 ELSE 1265 E(1,:) = 0._dp 1266 END IF 1267 localV=0._dp 1268 1269 SELECT CASE (CoilType) 1270 1271 CASE ('stranded') 1272 SELECT CASE(dim) 1273 CASE(2) 1274 wvec = [0._dp, 0._dp, 1._dp] 1275 IF(CMat_ip(1,1) /= 0._dp ) & 1276 E(1,:) = E(1,:)+ LagrangeVar % Values(IvarId) * N_j * wvec / CMat_ip(1,1) 1277 CASE(3) 1278 wvec = -MATMUL(Wbase(1:np), dBasisdx(1:np,:)) 1279 wvec = wvec/SQRT(SUM(wvec**2._dp)) 1280 IF(CMat_ip(3,3) /= 0._dp ) & 1281 E(1,:) = E(1,:)+ LagrangeVar % Values(IvarId) * N_j * wvec / CMat_ip(3,3) 1282 END SELECT 1283 1284 CASE ('massive') 1285 localV(1) = localV(1) + LagrangeVar % Values(VvarId) * CircEqVoltageFactor 1286 SELECT CASE(dim) 1287 CASE(2) 1288 E(1,3) = E(1,3)-localV(1) * grads_coeff 1289 CASE(3) 1290 E(1,:) = E(1,:)-localV(1) * MATMUL(Wbase(1:np), dBasisdx(1:np,:)) 1291 END SELECT 1292 1293 CASE ('foil winding') 1294 localAlpha = coilthickness *SUM(alpha(1:np) * Basis(1:np)) 1295 DO k = 1, VvarDofs-1 1296 localV(1) = localV(1) + LagrangeVar % Values(VvarId+k) * localAlpha**(k-1) * CircEqVoltageFactor 1297 END DO 1298 SELECT CASE(dim) 1299 CASE(2) 1300 E(1,3) = E(1,3)-localV(1) * grads_coeff 1301 CASE(3) 1302 E(1,:) = E(1,:)-localV(1) * MATMUL(Wbase(1:np), dBasisdx(1:np,:)) 1303 END SELECT 1304 1305 CASE DEFAULT 1306 SELECT CASE(dim) 1307 CASE(2) 1308 IF (HasLorenzVelocity) THEN 1309 ! Add v x curl A 1310 IF (CSymmetry) THEN 1311 E(1,3) = E(1,3) - rot_velo(1) * B(1,2) + rot_velo(2) * B(1,1) 1312 ELSE 1313 E(1,3) = E(1,3) + rot_velo(1) * B(1,2) - rot_velo(2) * B(1,1) 1314 END IF 1315 END IF 1316 ! 1317 ! To make this perfect, the electric field corresponding to the source 1318 ! should be returned on the source region 1319 ! 1320 CASE(3) 1321 IF (Transient) THEN 1322 E(1,:) = E(1,:)-MATMUL(SOL(1,1:np), dBasisdx(1:np,:)) 1323 END IF 1324 1325 IF (np > 0 .AND. .NOT. Transient) THEN 1326 E(1,:) = -MATMUL(SOL(1,1:np), dBasisdx(1:np,:)) 1327 1328 IF (HasVelocity) THEN 1329 ! 1330 ! Add v x curl A so as to handle the steady state solution of 1331 ! the evolutionary equations. Multiplication with the electric 1332 ! conductivity will give the current density with respect to 1333 ! the fixed frame: 1334 ! 1335 E(1,:) = E(1,:) + CrossProduct(rot_velo, & 1336 MATMUL(SOL(1,np+1:nd), RotWBasis(1:nd-np,:))) 1337 END IF 1338 ELSE IF ( PrecomputedElectricPot ) THEN 1339 E(1,:) = -MATMUL(ElPotSol(1,1:n), dBasisdx(1:n,:)) 1340 END IF 1341 1342 IF ( ImposeBodyForcePotential ) THEN 1343 E(1,:) = E(1,:) - MATMUL(ElPotSol(1,1:n), dBasisdx(1:n,:)) 1344 END IF 1345 END SELECT 1346 1347 1348 END SELECT 1349 END IF 1350 1351 Nu = CMPLX(0.0d0, 0.0d0, kind=dp) 1352 IF ( ASSOCIATED(HB) ) THEN 1353 IF (RealField) THEN 1354 Babs=SQRT(SUM(B(1,:)**2)) 1355 ELSE 1356 Babs = SQRT(SUM(B(1,:)**2 + B(2,:)**2)) 1357 END IF 1358 Babs = MAX(Babs, 1.d-8) 1359 R_ip = InterpolateCurve(HBBval,HBHval,Babs,HBCval)/Babs 1360 w_dens = IntegrateCurve(HBBval,HBHval,HBCval,0._dp,Babs) 1361 DO k=1,3 1362 Nu(k,k) = CMPLX(R_ip, 0.0d0, kind=dp) 1363 END DO 1364 ELSE 1365 IF (HasTensorReluctivity) THEN 1366 IF (SIZE(Reluct_Z,2) == 1) THEN 1367 DO k = 1, MIN(3, SIZE(Reluct_Z,1)) 1368 Nu(k,k) = SUM(Basis(1:n)*Reluct_Z(k,1,1:n)) 1369 END DO 1370 ELSE 1371 DO k = 1, MIN(3, SIZE(Reluct_Z,1)) 1372 DO l = 1, MIN(3, SIZE(Reluct_Z,2)) 1373 Nu(k,l) = sum(Basis(1:n)*Reluct_Z(k,l,1:n)) 1374 END DO 1375 END DO 1376 END IF 1377 R_ip = 0.0d0 1378 ELSE 1379 R_ip_Z = SUM(Basis(1:n)*R_Z(1:n)) 1380 DO k=1,3 1381 Nu(k,k) = R_ip_Z 1382 END DO 1383 ! 1384 ! The calculation of the Maxwell stress tensor doesn't yet support 1385 ! a tensor-form reluctivity. Create the scalar reluctivity parameter 1386 ! so that the Maxwell stress tensor may be calculated. The complex 1387 ! part will be ignored. 1388 ! 1389 R_ip = REAL(R_ip_Z) 1390 END IF 1391 IF (RealField) THEN 1392 w_dens = 0.5*SUM(B(1,:)*MATMUL(REAL(Nu), B(1,:))) 1393 ELSE 1394 ! This yields twice the time average: 1395 w_dens = 0.5*( SUM(MATMUL(REAL(Nu), B(1,:)) * B(1,:)) + & 1396 SUM(MATMUL(REAL(Nu), B(2,:)) * B(2,:)) ) 1397 END IF 1398 END IF 1399 PR_ip = SUM( Basis(1:n)*PR(1:n) ) 1400 1401 IF ( ASSOCIATED(MFS).OR.ASSOCIATED(EL_MFS) ) THEN 1402 DO l=1,3 1403 MG_ip(l) = SUM( Magnetization(l,1:n)*Basis(1:n) ) 1404 END DO 1405 END IF 1406 1407 IF( ASSOCIATED(VP).OR.ASSOCIATED(EL_VP) ) THEN 1408 DO l=1,vDOFs 1409 SELECT CASE(dim) 1410 CASE(2) 1411 VP_ip(l,1) = 0._dp 1412 VP_ip(l,2) = 0._dp 1413 VP_ip(l,3) = SUM(SOL(l,1:nd) * Basis(1:nd)) 1414 CASE(3) 1415 VP_ip(l,:)=MATMUL(SOL(l,np+1:nd),WBasis(1:nd-np,:)) 1416 END SELECT 1417 END DO 1418 END IF 1419 1420 IF (RealField) THEN 1421 HdotB = SUM(MATMUL(REAL(Nu), B(1,:)) * B(1,:)) 1422 ELSE 1423 HdotB = SUM(MATMUL(REAL(Nu), B(1,:)) * B(1,:)) + & 1424 SUM(MATMUL(REAL(Nu), B(2,:)) * B(2,:)) 1425 END IF 1426 1427 IF (ASSOCIATED(NF).OR.ASSOCIATED(EL_NF)) THEN 1428 NF_ip = 0._dp 1429 DO k=1,n 1430 DO l=1,3 1431 val = SUM(dBasisdx(k,1:3)*B(1,1:3)) 1432 NF_ip(k,l) = NF_ip(k,l) - SUM(REAL(Nu(l,1:3)) * B(1,1:3)) * val + & 1433 (HdotB-w_dens)*dBasisdx(k,l) 1434 END DO 1435 END DO 1436 1437 IF (.NOT. RealField) THEN 1438 DO k=1,n 1439 DO l=1,3 1440 val = SUM(dBasisdx(k,1:3)*B(2,1:3)) 1441 NF_ip(k,l) = NF_ip(k,l) - SUM(REAL(Nu(l,1:3)) * B(2,1:3)) * val 1442 END DO 1443 END DO 1444 END IF 1445 END IF 1446 1447 Energy(1) = Energy(1) + s*0.5*PR_ip*SUM(E**2) 1448 Energy(2) = Energy(2) + s*w_dens 1449 Energy(3) = Energy(3) + (HdotB - w_dens) * s 1450 1451 DO p=1,n 1452 DO q=1,n 1453 MASS(p,q)=MASS(p,q)+s*Basis(p)*Basis(q) 1454 END DO 1455 k = 0 1456 DO l=1,vDOFs 1457 FORCE(p,k+1:k+3) = FORCE(p,k+1:k+3)+s*B(l,:)*Basis(p) 1458 k = k+3 1459 END DO 1460 1461 IF ( (ASSOCIATED(MFS).OR.ASSOCIATED(EL_MFS)) .and. .not. HasZirka) THEN 1462 IF(.NOT. HasZirka) then 1463 IF (RealField) THEN 1464 FORCE(p,k+1:k+3) = FORCE(p,k+1:k+3) + & 1465 s * (MATMUL(REAL(Nu), B(1,:)) - REAL(MG_ip)) * Basis(p) 1466 k = k+3 1467 ELSE 1468 FORCE(p,k+1:k+3) = FORCE(p,k+1:k+3) + s * & 1469 (MATMUL(REAL(Nu), B(1,:)) - MATMUL(AIMAG(Nu), B(2,:)) - REAL(MG_ip)) * Basis(p) 1470 k = k+3 1471 FORCE(p,k+1:k+3) = FORCE(p,k+1:k+3) + s * & 1472 (MATMUL(AIMAG(Nu), B(1,:)) + MATMUL(REAL(Nu), B(2,:)) - AIMAG(MG_ip)) * Basis(p) 1473 k = k+3 1474 END IF 1475 ELSE 1476 ! Never here? 1477 FORCE(p,k+1:k+3) = FORCE(p,k+1:k+3)-s*(REAL(MG_ip))*Basis(p) 1478 END IF 1479 END IF 1480 IF ( ASSOCIATED(VP).OR.ASSOCIATED(EL_VP)) THEN 1481 DO l=1,vDOFs 1482 FORCE(p,k+1:k+3) = FORCE(p,k+1:k+3)+s*VP_ip(l,:)*Basis(p) 1483 k = k+3 1484 END DO 1485 END IF 1486 IF ( ASSOCIATED(EF).OR.ASSOCIATED(EL_EF)) THEN 1487 DO l=1,vDOFs 1488 FORCE(p,k+1:k+3) = FORCE(p,k+1:k+3)+s*E(l,:)*Basis(p) 1489 k = k+3 1490 END DO 1491 END IF 1492 1493 IF ( ASSOCIATED(CD).OR.ASSOCIATED(EL_CD)) THEN 1494 IF (ItoJCoeffFound) THEN 1495 IF (Vdofs == 1) THEN 1496 DO l=1,3 1497 CC_J(1,l) = ItoJCoeff*wvec(l)*CircuitCurrent 1498 END DO 1499 ELSE 1500 CALL Fatal('MagnetoDynamicsCalcFields','Complex circuit current imposing is not implemented') 1501 END IF 1502 ELSE 1503 CC_J(1,:) = 0.0_dp 1504 END IF 1505 1506 IF (Vdofs == 1) THEN 1507 DO l=1,3 1508 JatIP(1,l) = SUM( REAL(CMat_ip(l,1:3)) * E(1,1:3) ) + CC_J(1,l) + REAL(BodyForceCurrDens_ip(l)) 1509 ! 1510 ! No need for a "HasVelocity" check: the effect of v x B is already inbuilt into 1511 ! the definition of the E-field 1512 ! IF( HasVelocity ) THEN 1513 ! JatIP(1,l) = JatIP(1,l) + SUM( REAL(CMat_ip(l,1:3)) * CrossProduct(rot_velo, B(1,1:3))) 1514 ! END IF 1515 ! 1516 FORCE(p,k+l) = FORCE(p,k+l)+s*JatIp(1,l)*Basis(p) 1517 END DO 1518 k = k+3 1519 ELSE 1520 DO l=1,3 1521 JatIp(1,l) = SUM( REAL(CMat_ip(l,1:3)) * E(1,1:3) ) - & 1522 SUM( AIMAG(CMat_ip(l,1:3)) * E(2,1:3) ) + REAL(BodyForceCurrDens_ip(l)) 1523 ! 1524 ! No need for a "HasVelocity" check: the effect of v x B is already inbuilt into 1525 ! the definition of the E-field 1526 ! IF( HasVelocity ) THEN 1527 ! JatIp(1,l) = JatIp(1,l) + SUM( REAL(CMat_ip(l,1:3)) * CrossProduct(rot_velo, B(1,1:3)) ) - & 1528 ! SUM( AIMAG(CMat_ip(l,1:3)) * CrossProduct(rot_velo, B(2,1:3)) ) 1529 ! END IF 1530 ! 1531 FORCE(p,k+l) = FORCE(p,k+l)+s*JatIp(1,l)*Basis(p) 1532 END DO 1533 k = k+3 1534 DO l=1,3 1535 JatIp(2,l) = SUM( AIMAG(CMat_ip(l,1:3)) * E(1,1:3) ) + & 1536 SUM( REAL(CMat_ip(l,1:3)) * E(2,1:3) ) + AIMAG(BodyForceCurrDens_ip(l)) 1537 ! 1538 ! No need for a "HasVelocity" check: the effect of v x B is already inbuilt into 1539 ! the definition of the E-field 1540 ! IF( HasVelocity ) THEN 1541 ! JatIp(2,l) = JatIp(2,l) + SUM( AIMAG(CMat_ip(l,1:3)) * CrossProduct(rot_velo, B(1,1:3)) ) + & 1542 ! SUM( REAL(CMat_ip(l,1:3)) * CrossProduct(rot_velo, B(2,1:3)) ) 1543 ! END IF 1544 ! 1545 FORCE(p,k+l) = FORCE(p,k+l)+s*JatIp(2,l)*Basis(p) 1546 END DO 1547 k = k+3 1548 END IF 1549 END IF 1550 1551 IF ( ASSOCIATED(JXB).OR.ASSOCIATED(EL_JXB)) THEN 1552 IF (.NOT. ASSOCIATED(CD) .AND. .NOT. ASSOCIATED(EL_CD)) THEN 1553 CALL Warn('MagnetoDynamicsCalcFields', 'Cannot Calculate JxB since Current Density is not calculated!') 1554 ELSE 1555 IF (Vdofs == 1) THEN 1556 JXBatIP(1,:) = crossproduct(JatIP(1,:),B(1,:)) 1557 DO l=1,dim 1558 FORCE(p,k+l) = FORCE(p,k+l)+s*JXBatIP(1,l)*Basis(p) 1559 END DO 1560 k = k+3 1561 ELSE 1562 IF( CSymmetry ) THEN 1563 ! TODO: Have to figure out why cylindrical coords have opposite sign 1564 JXBatIP(1,1) = JatIP(2,3)*B(2,2) + JatIP(1,3)*B(1,2) 1565 JXBatIP(1,2) = - JatIP(2,3)*B(2,1) - JatIP(1,3)*B(1,1) 1566 JXBatIP(1,3) = 0.0_dp 1567 1568 JXBatIP(2,1) = JatIP(2,3)*B(1,2) - JatIP(1,3)*B(2,2) 1569 JXBatIP(2,2) = - JatIP(2,3)*B(1,1) + JatIP(1,3)*B(2,1) 1570 JXBatIP(2,3) = 0.0_dp 1571 ELSE 1572 JXBatIP(1,1) = JatIP(2,2)*B(2,3) - JatIP(2,3)*B(2,2) + JatIP(1,2)*B(1,3) - JatIP(1,3)*B(1,2) 1573 JXBatIP(1,2) = JatIP(2,3)*B(2,1) - JatIP(2,1)*B(2,3) + JatIP(1,3)*B(1,1) - JatIP(1,1)*B(1,3) 1574 JXBatIP(1,3) = JatIP(2,1)*B(2,2) - JatIP(2,2)*B(2,1) + JatIP(1,1)*B(1,2) - JatIP(1,2)*B(1,1) 1575 1576 JXBatIP(2,1) = JatIP(2,2)*B(1,3) - JatIP(2,3)*B(1,2) - JatIP(1,2)*B(2,3) + JatIP(1,3)*B(2,2) 1577 JXBatIP(2,2) = JatIP(2,3)*B(1,1) - JatIP(2,1)*B(1,3) - JatIP(1,3)*B(2,1) + JatIP(1,1)*B(2,3) 1578 JXBatIP(2,3) = JatIP(2,1)*B(1,2) - JatIP(2,2)*B(1,1) - JatIP(1,1)*B(2,2) + JatIP(1,2)*B(2,1) 1579 END IF 1580 1581 JXBatIP = 0.5_dp*JXBatIP 1582 1583 DO l=1,dim 1584 FORCE(p,k+l) = FORCE(p,k+l)+s*JXBatIP(1,l)*Basis(p) 1585 END DO 1586 k = k+3 1587 DO l=1,dim 1588 FORCE(p,k+l) = FORCE(p,k+l)+s*JXBatIP(2,l)*Basis(p) 1589 END DO 1590 k = k+3 1591 END IF 1592 END IF 1593 END IF 1594 1595 IF ( ASSOCIATED(FWP).OR.ASSOCIATED(EL_FWP)) THEN 1596 IF (Vdofs == 1) THEN 1597 FORCE(p,k+1) = FORCE(p,k+1)+s*LocalV(1)*Basis(p) 1598 k = k+1 1599 ELSE 1600 FORCE(p,k+1) = FORCE(p,k+1)+s*LocalV(1) * Basis(p) 1601 k = k+1 1602 FORCE(p,k+1) = FORCE(p,k+1)+s*LocalV(2) * Basis(p) 1603 k = k+1 1604 END IF 1605 END IF 1606 1607 IF (vDOFS == 1) THEN 1608 IF ( JouleHeatingFromCurrent .AND. (ASSOCIATED(CD).OR.ASSOCIATED(EL_CD)) ) THEN 1609 ! The Joule heating power per unit volume: J.E = J.J/sigma 1610 Coeff = 0.0_dp 1611 DO l=1,3 1612 IF( REAL( CMat_ip(l,l) ) > EPSILON( Coeff ) ) THEN 1613 Coeff = Coeff + JatIP(1,l) * JatIP(1,l) / REAL( CMat_ip(l,l) ) * & 1614 Basis(p) * s 1615 END IF 1616 END DO 1617 ELSE 1618 ! The Joule heating power per unit volume: J.E = (sigma * E).E 1619 Coeff = SUM( MATMUL( REAL(CMat_ip(1:3,1:3)), TRANSPOSE(E(1:1,1:3)) ) * & 1620 TRANSPOSE(E(1:1,1:3)) ) * Basis(p) * s 1621 END IF 1622 ! 1623 ! 1624 ! No need for a "HasVelocity" check: the effect of v x B is already inbuilt into 1625 ! the definition of the J-field via J's dependence on the E-field 1626 ! 1627 ! IF (HasVelocity) THEN 1628 ! Coeff = Coeff + SUM(MATMUL(REAL(CMat_ip), CrossProduct(rot_velo, B(1,:))) * & 1629 ! CrossProduct(rot_velo,B(1,:)))*Basis(p)*s 1630 ! END IF 1631 ! 1632 ELSE 1633 ! Now Power = J.conjugate(E), with the possible imaginary component neglected. 1634 Coeff = HarmPowerCoeff * (SUM( MATMUL( REAL(CMat_ip(1:3,1:3)), TRANSPOSE(E(1:1,1:3)) ) * & 1635 TRANSPOSE(E(1:1,1:3)) ) * Basis(p) * s - & 1636 SUM( MATMUL( AIMAG(CMat_ip(1:3,1:3)), TRANSPOSE(E(2:2,1:3)) ) * & 1637 TRANSPOSE(E(1:1,1:3)) ) * Basis(p) * s + & 1638 SUM( MATMUL( AIMAG(CMat_ip(1:3,1:3)), TRANSPOSE(E(1:1,1:3)) ) * & 1639 TRANSPOSE(E(2:2,1:3)) ) * Basis(p) * s + & 1640 SUM( MATMUL( REAL(CMat_ip(1:3,1:3)), TRANSPOSE(E(2:2,1:3)) ) * & 1641 TRANSPOSE(E(2:2,1:3)) ) * Basis(p) * s) 1642 ! 1643 ! Again, no need for a "HasVelocity" check 1644 ! IF (HasVelocity) THEN 1645 ! Coeff = Coeff + HarmPowerCoeff * (SUM( MATMUL( REAL(CMat_ip(1:3,1:3)), CrossProduct(rot_velo, B(1,:)) ) * & 1646 ! CrossProduct(rot_velo, B(1,:)) ) * Basis(p) * s - & 1647 ! SUM( MATMUL( AIMAG(CMat_ip(1:3,1:3)), CrossProduct(rot_velo, B(2,:)) ) * & 1648 ! CrossProduct(rot_velo, B(1,:)) ) * Basis(p) * s + & 1649 ! SUM( MATMUL( AIMAG(CMat_ip(1:3,1:3)), CrossProduct(rot_velo, B(1,:)) ) * & 1650 ! CrossProduct(rot_velo, B(2,:)) ) * Basis(p) * s + & 1651 ! SUM( MATMUL( REAL(CMat_ip(1:3,1:3)), CrossProduct(rot_velo, B(2,:)) ) * & 1652 ! CrossProduct(rot_velo, B(2,:)) ) * Basis(p) * s) 1653 ! END IF 1654 ! 1655 END IF 1656 1657 IF(ALLOCATED(BodyLoss)) BodyLoss(3,BodyId) = BodyLoss(3,BodyId) + Coeff 1658 Power = Power + Coeff 1659 IF ( ASSOCIATED(JH) .OR. ASSOCIATED(EL_JH) .OR. ASSOCIATED(NJH) ) THEN 1660 FORCE(p,k+1) = FORCE(p,k+1) + Coeff 1661 k = k+1 1662 END IF 1663 1664 1665 1666 !------------------------------------------------- 1667 ! Compute a loss estimate for cos and sin modes: 1668 !------------------------------------------------- 1669 IF (LossEstimation) THEN 1670 LossCoeff = ListGetFun( Material,'Harmonic Loss Linear Coefficient',Freq,Found ) 1671 LossCoeff2 = ListGetFun( Material,'Harmonic Loss Quadratic Coefficient',Freq,Found ) 1672 ! No losses to add if loss coefficient is not given 1673 IF( Found ) THEN 1674 DO l=1,2 1675 ValAtIP = SUM( B(l,1:3) ** 2 ) 1676 Coeff = s * Basis(p) * LossCoeff * ( Freq ** FreqPower ) * ( ValAtIp ** FieldPower ) 1677 Coeff2 = s * Basis(p) * LossCoeff2 * ( Freq ** FreqPower2 ) * ( ValAtIp ** FieldPower2 ) 1678 ComponentLoss(1,l) = ComponentLoss(1,l) + Coeff 1679 BodyLoss(1,BodyId) = BodyLoss(1,BodyId) + Coeff 1680 ComponentLoss(2,l) = ComponentLoss(2,l) + Coeff2 1681 BodyLoss(2,BodyId) = BodyLoss(2,BodyId) + Coeff2 1682 END DO 1683 ELSE 1684 Coeff = 0.0_dp 1685 Coeff2 = 0.0_dp 1686 END IF 1687 1688 IF ( ASSOCIATED(ML) .OR. ASSOCIATED(EL_ML) ) THEN 1689 FORCE(p,k+1) = FORCE(p,k+1) + Coeff 1690 k = k + 1 1691 END IF 1692 IF ( ASSOCIATED(ML2) .OR. ASSOCIATED(EL_ML2) ) THEN 1693 FORCE(p,k+1) = FORCE(p,k+1) + Coeff2 1694 k = k + 1 1695 END IF 1696 END IF 1697 1698 IF ( ASSOCIATED(MST).OR.ASSOCIATED(EL_MST)) THEN 1699 IF ( Vdofs==1 ) THEN 1700 DO l=1,3 1701 DO m=l,3 1702 ST(l,m)=PR_ip*E(1,l)*E(1,m)+R_ip*B(1,l)*B(1,m) 1703 END DO 1704 ST(l,l)=ST(l,l)-(PR_ip*SUM(E(1,:)**2)+R_ip*SUM(B(1,:)**2))/2 1705 END DO 1706 DO l=1,6 1707 FORCE(p,k+l)=FORCE(p,k+l) + s*ST(ind1(l),ind2(l))*Basis(p) 1708 END DO 1709 k = k + 6 1710 ELSE 1711 DO l=1,3 1712 DO m=l,3 1713 CST(l,m) = PR_ip*CMPLX(E(1,l),E(2,l),KIND=dp) * & 1714 CMPLX(E(1,m),E(2,m),KIND=dp) 1715 CST(l,m) = CST(l,m) + & 1716 R_ip*CMPLX(B(1,l),B(2,l),KIND=dp) * & 1717 CMPLX(B(1,m),B(2,m),KIND=dp) 1718 END DO 1719 CST(l,l) = CST(l,l) - & 1720 (PR_ip*SUM(ABS(CMPLX(E(1,:),E(2,:)))**2)+ & 1721 R_ip*SUM(ABS(CMPLX(B(1,:),B(2,:)))**2))/2 1722 END DO 1723 DO l=1,6 1724 FORCE(p,k+l)=FORCE(p,k+l) + s*REAL(CST(ind1(l),ind2(l)))*Basis(p) 1725 END DO 1726 k = k + 6 1727 DO l=1,6 1728 FORCE(p,k+l)=FORCE(p,k+l) + s*AIMAG(CST(ind1(l),ind2(l)))*Basis(p) 1729 END DO 1730 k = k + 6 1731 END IF 1732 END IF 1733 IF (ASSOCIATED(NF).OR.ASSOCIATED(EL_NF)) THEN 1734 IF(RealField) THEN 1735 FORCE(p,k+1:k+3) = FORCE(p,k+1:k+3) + s*NF_ip(p,1:3) 1736 ELSE 1737 FORCE(p,k+1:k+3) = FORCE(p,k+1:k+3) + 0.5*s*NF_ip(p,1:3) 1738 END IF 1739 k = k + 3 1740 END IF 1741 END DO ! p 1742 END DO ! j 1743 1744 1745 IF(NodalFields) THEN 1746 CALL DefaultUpdateEquations( MASS,Force(:,1)) 1747 Fsave => Solver % Matrix % RHS 1748 DO l=1,k 1749 Solver % Matrix % RHS => GForce(:,l) 1750 CALL DefaultUpdateForce(Force(:,l)) 1751 END DO 1752 Solver % Matrix % RHS => Fsave 1753 END IF 1754 1755 IF(ElementalFields) THEN 1756 dofs = 0 1757 CALL LUdecomp(MASS,n,pivot) 1758 CALL LocalSol(EL_MFD, 3*vdofs, n, MASS, FORCE, pivot, Dofs) 1759 CALL LocalSol(EL_MFS, 3*vdofs, n, MASS, FORCE, pivot, Dofs) 1760 CALL LocalSol(EL_VP, 3*vdofs, n, MASS, FORCE, pivot, Dofs) 1761 CALL LocalSol(EL_EF, 3*vdofs, n, MASS, FORCE, pivot, Dofs) 1762 CALL LocalSol(EL_CD, 3*vdofs, n, MASS, FORCE, pivot, Dofs) 1763 CALL LocalSol(EL_JXB, 3*vdofs, n, MASS, FORCE, pivot, Dofs) 1764 CALL LocalSol(EL_FWP, 1*vdofs, n, MASS, FORCE, pivot, Dofs) 1765 CALL LocalSol(EL_JH, 1, n, MASS, FORCE, pivot, Dofs) 1766 CALL LocalSol(EL_ML, 1, n, MASS, FORCE, pivot, Dofs) 1767 CALL LocalSol(EL_ML2, 1, n, MASS, FORCE, pivot, Dofs) 1768 CALL LocalSol(EL_MST, 6*vdofs, n, MASS, FORCE, pivot, Dofs) 1769 1770 ! This is a nodal quantity 1771 CALL LocalCopy(EL_NF, 3, n, FORCE, Dofs) 1772 END IF 1773 END DO 1774 1775 1776 ! Assembly of the face terms: 1777 !---------------------------- 1778 IF (GetLogical(SolverParams,'Discontinuous Galerkin',Found)) THEN 1779 IF (GetLogical(SolverParams,'Average Within Materials',Found)) THEN 1780 FORCE = 0.0_dp 1781 CALL AddLocalFaceTerms( MASS, FORCE(:,1) ) 1782 END IF 1783 END IF 1784 1785 1786 IF(NodalFields) THEN 1787 Fsave => Solver % Matrix % RHS 1788 DOFs = 0 1789 CALL GlobalSol(MFD, 3*vdofs, Gforce, Dofs) 1790 CALL GlobalSol(MFS, 3*vdofs, Gforce, Dofs) 1791 CALL GlobalSol(VP , 3*vdofs, Gforce, Dofs) 1792 CALL GlobalSol(EF, 3*vdofs, Gforce, Dofs) 1793 CALL GlobalSol(CD, 3*vdofs, Gforce, Dofs) 1794 CALL GlobalSol(JXB, 3*vdofs, Gforce, Dofs) 1795 CALL GlobalSol(FWP, 1*vdofs, Gforce, Dofs) 1796 1797 ! Nodal heating directly uses the loads 1798 IF (ASSOCIATED(NJH)) THEN 1799 NJH % Values = Gforce(:,dofs+1) 1800 ! Update the dofs only if it not used as the r.h.s. for the following field 1801 IF(.NOT. ASSOCIATED(JH) ) dofs = dofs + 1 1802 END IF 1803 CALL GlobalSol(JH , 1 , Gforce, Dofs) 1804 1805 CALL GlobalSol(ML , 1 , Gforce, Dofs) 1806 CALL GlobalSol(ML2, 1 , Gforce, Dofs) 1807 CALL GlobalSol(MST, 6*vdofs, Gforce, Dofs) 1808 IF (ASSOCIATED(NF)) THEN 1809 DO i=1,3 1810 dofs = dofs + 1 1811 NF % Values(i::3) = Gforce(:,dofs) 1812 END DO 1813 END IF 1814 Solver % Matrix % RHS => Fsave 1815 END IF 1816 1817 1818 ! Lump componentwise forces and torques. 1819 ! Prefer DG nodal force variable if air gap is present 1820 1821 ! Warn if user has air gaps and no "nodal force e" 1822 HaveAirGap = ListCheckPresentAnyBC( Model, 'Air Gap Length' ) 1823 UseElementalNF = ASSOCIATED( EL_NF ) .AND. ( .NOT. ASSOCIATED( NF ) .OR. HaveAirGap ) 1824 1825 1826 IF( UseElementalNF ) THEN 1827 1828 ! Collect nodal forces from airgaps 1829 CALL CalcBoundaryModels() 1830 1831 ! Create a minimal discontinuous set such that discontinuity is only created 1832 ! when body has an air gap boundary condition. Only do the reduction for the 1st time. 1833 IF( .NOT. ASSOCIATED( SetPerm ) ) THEN 1834 CALL Info('MagnetoDynamicsCalcFields','Creating minimal elemental set',Level=10) 1835 SetPerm => MinimalElementalSet( Mesh,'db', Solver % Variable % Perm, & 1836 BcFlag = 'Air Gap Length', & 1837 NonGreedy = ListGetLogical( Solver % Values,'Nongreedy Jump',Found) ) 1838 END IF 1839 1840 ! Sum up (no averaging) the elemental fields such that each elemental nodes has also 1841 ! the contributions of the related nodes in other elements 1842 CALL ReduceElementalVar( Mesh, EL_NF, SetPerm, TakeAverage = .FALSE.) 1843 1844 DO j=1,Model % NumberOfComponents 1845 CompParams => Model % Components(j) % Values 1846 1847 IF ( ListGetLogical( CompParams,'Calculate Magnetic Force', Found ) ) THEN 1848 1849 CALL ComponentNodalForceReduction(Model, Mesh, CompParams, EL_NF, & 1850 Force = LumpedForce, SetPerm = SetPerm ) 1851 1852 WRITE( Message,'(A,3ES15.6)') 'Magnetic force reduced: > '& 1853 //TRIM(ListGetString(CompParams,'Name'))//' < :', LumpedForce 1854 CALL Info('MagnetoDynamicsCalcFields',Message,Level=6) 1855 1856 DO i=1,3 1857 CALL ListAddConstReal( CompParams,'res: magnetic force '//TRIM(I2S(i)), LumpedForce(i) ) 1858 END DO 1859 1860 END IF 1861 1862 IF( ListGetLogical( CompParams,'Calculate Magnetic Torque', Found ) ) THEN 1863 1864 CALL ComponentNodalForceReduction(Model, Mesh, CompParams, EL_NF, & 1865 Torque = val, SetPerm = SetPerm ) 1866 1867 WRITE( Message,'(A,ES15.6)') 'Magnetic torque reduced: > '& 1868 //TRIM(ListGetString(CompParams,'Name'))//' < :', val 1869 CALL Info('MagnetoDynamicsCalcFields',Message,Level=6) 1870 1871 CALL ListAddConstReal( CompParams,'res: magnetic torque', val ) 1872 END IF 1873 END DO 1874 ELSE 1875 DO j=1,Model % NumberOfComponents 1876 CompParams => Model % Components(j) % Values 1877 1878 IF( ListGetLogical( CompParams,'Calculate Magnetic Force', Found ) ) THEN 1879 1880 ! fail if there is no nodal force variable available 1881 IF (.NOT. ASSOCIATED(NF)) THEN 1882 CALL Warn('MagnetoDynamicsCalcFields','Unable to calculated lumped & 1883 &forces because nodal forces are not present. Use keyword & 1884 &"Calculate Nodal Forces = true" in MagnetoDynamicsCalcFields solver.') 1885 EXIT 1886 END IF 1887 1888 CALL ComponentNodalForceReduction(Model, Mesh, CompParams, NF, & 1889 Force = LumpedForce ) 1890 1891 WRITE( Message,'(A,3ES15.6)') 'Magnetic force reduced: > '& 1892 //TRIM(ListGetString(CompParams,'Name'))//' < :', LumpedForce 1893 CALL Info('MagnetoDynamicsCalcFields',Message,Level=6) 1894 1895 DO i=1,3 1896 CALL ListAddConstReal( CompParams,'res: magnetic force '//TRIM(I2S(i)), LumpedForce(i) ) 1897 END DO 1898 1899 END IF 1900 1901 IF( ListGetLogical( CompParams,'Calculate Magnetic Torque', Found ) ) THEN 1902 1903 ! fail if there is no nodal force variable available 1904 IF (.NOT. ASSOCIATED(NF)) THEN 1905 CALL Warn('MagnetoDynamicsCalcFields','Unable to calculated lumped & 1906 &forces because nodal forces are not present. Use keyword & 1907 &"Calculate Nodal Forces = true" in MagnetoDynamicsCalcFields solver.') 1908 EXIT 1909 END IF 1910 1911 ! Warn if user has air gaps and no "nodal force e" is available 1912 IF ( HaveAirGap ) THEN 1913 CALL Warn('MagnetoDynamicsCalcFields', 'Cannot calculate air gap & 1914 &forces correctly because elemental field "Nodal Force e" is not & 1915 &present.') 1916 END IF 1917 1918 CALL ComponentNodalForceReduction(Model, Mesh, CompParams, NF, & 1919 Torque = val ) 1920 1921 WRITE( Message,'(A,ES15.6)') 'Magnetic torque reduced: > '& 1922 //TRIM(ListGetString(CompParams,'Name'))//' < :', val 1923 CALL Info('MagnetoDynamicsCalcFields',Message,Level=6) 1924 1925 CALL ListAddConstReal( CompParams,'res: magnetic torque', val ) 1926 END IF 1927 END DO 1928 END IF 1929 1930 1931 ! Perform parallel reductions 1932 Power = ParallelReduction(Power) 1933 Energy(1) = ParallelReduction(Energy(1)) 1934 Energy(2) = ParallelReduction(Energy(2)) 1935 Energy(3) = ParallelReduction(Energy(3)) 1936 1937 IF (LossEstimation) THEN 1938 DO j=1,2 1939 DO i=1,2 1940 ComponentLoss(j,i) = ParallelReduction(ComponentLoss(j,i)) 1941 END DO 1942 END DO 1943 1944 DO j=1,3 1945 DO i=1,Model % NumberOfBodies 1946 BodyLoss(j,i) = ParallelReduction(BodyLoss(j,i)) 1947 END DO 1948 TotalLoss(j) = SUM( BodyLoss(j,:) ) 1949 END DO 1950 END IF 1951 1952 1953 WRITE(Message,*) 'Eddy current power: ', Power 1954 CALL Info( 'MagnetoDynamicsCalcFields', Message ) 1955 CALL ListAddConstReal( Model % Simulation, 'res: Eddy current power', Power ) 1956 1957 IF ( ListGetLogical( SolverParams,'Separate Magnetic Energy',Found ) ) THEN 1958 WRITE(Message,'(A,ES12.3)') 'Electric Field Energy: ', Energy(1) 1959 CALL Info( 'MagnetoDynamicsCalcFields', Message ) 1960 WRITE(Message,'(A,ES12.3)') 'Magnetic Field Energy: ', Energy(2) 1961 CALL Info( 'MagnetoDynamicsCalcFields', Message ) 1962 WRITE(Message,'(A,ES12.3)') 'Magnetic Coenergy: ', Energy(3) 1963 CALL Info( 'MagnetoDynamicsCalcFields', Message ) 1964 CALL ListAddConstReal(Model % Simulation,'res: Electric Field Energy',Energy(1)) 1965 CALL ListAddConstReal(Model % Simulation,'res: Magnetic Field Energy',Energy(2)) 1966 CALL ListAddConstReal(Model % Simulation,'res: Magnetic Coenergy',Energy(3)) 1967 ELSE 1968 WRITE(Message,'(A,ES12.3)') 'ElectroMagnetic Field Energy: ',SUM(Energy(1:2)) 1969 CALL Info( 'MagnetoDynamicsCalcFields', Message ) 1970 CALL ListAddConstReal(Model % Simulation,'res: ElectroMagnetic Field Energy',SUM(Energy(1:2))) 1971 END IF 1972 1973 IF(ALLOCATED(Gforce)) DEALLOCATE(Gforce) 1974 DEALLOCATE( MASS,FORCE,Tcoef,RotM ) 1975 1976 IF (LossEstimation) THEN 1977 CALL ListAddConstReal( Model % Simulation,'res: harmonic loss linear',TotalLoss(1) ) 1978 CALL ListAddConstReal( Model % Simulation,'res: harmonic loss quadratic',TotalLoss(2) ) 1979 CALL ListAddConstReal( Model % Simulation,'res: joule loss',TotalLoss(3) ) 1980 1981 DO k=1,2 1982 IF( k == 1 ) THEN 1983 CALL Info('MagnetoDynamicsCalcFields','Harmonic Loss Linear by components',Level=6) 1984 ELSE 1985 CALL Info('MagnetoDynamicsCalcFields','Harmonic Loss Quadratic by components',Level=6) 1986 END IF 1987 WRITE( Message,'(A,ES12.3)') 'Loss for cos mode: ', ComponentLoss(k,1) 1988 CALL Info('MagnetoDynamicsCalcFields', Message, Level=6 ) 1989 WRITE( Message,'(A,ES12.3)') 'Loss for sin mode: ', ComponentLoss(k,2) 1990 CALL Info('MagnetoDynamicsCalcFields', Message, Level=6 ) 1991 WRITE( Message,'(A,ES12.3)') 'Total loss: ',TotalLoss(k) 1992 CALL Info('MagnetoDynamicsCalcFields',Message, Level=5 ) 1993 END DO 1994 1995 DO k=1,3 1996 IF( TotalLoss(k) < TINY( TotalLoss(k) ) ) CYCLE 1997 IF( k == 1 ) THEN 1998 CALL Info('MagnetoDynamicsCalcFields','Harmonic Loss Linear by bodies',Level=6) 1999 ELSE IF( k == 2 ) THEN 2000 CALL Info('MagnetoDynamicsCalcFields','Harmonic Loss Quadratic by bodies',Level=6) 2001 ELSE 2002 CALL Info('MagnetoDynamicsCalcFields','Joule Loss by bodies',Level=6) 2003 END IF 2004 2005 DO j=1,Model % NumberOfBodies 2006 IF( BodyLoss(k,j) < TINY( TotalLoss(k) ) ) CYCLE 2007 WRITE( Message,'(A,I0,A,ES12.3)') 'Body ',j,' : ',BodyLoss(k,j) 2008 CALL Info('MagnetoDynamicsCalcFields', Message, Level=6 ) 2009 END DO 2010 2011 ! Save losses to components if requested. 2012 !--------------------------------------------------------------------------- 2013 DO j=1,Model % NumberOfComponents 2014 CompParams => Model % Components(j) % Values 2015 IF( ListGetLogical( CompParams,'Calculate Magnetic Losses', Found ) ) THEN 2016 MasterBodies => ListGetIntegerArray( CompParams,'Master Bodies',Found ) 2017 IF(.NOT. Found ) CYCLE 2018 val = SUM( BodyLoss(k,MasterBodies) ) 2019 IF( k == 1 ) THEN 2020 CALL ListAddConstReal( CompParams,'res: harmonic loss linear',val ) 2021 ELSE IF( k == 2 ) THEN 2022 CALL ListAddConstReal( CompParams,'res: harmonic loss quadratic',val ) 2023 ELSE 2024 CALL ListAddConstReal( CompParams,'res: joule loss',val ) 2025 END IF 2026 END IF 2027 END DO 2028 END DO 2029 2030 IF( ParEnv % MyPe == 0 ) THEN 2031 LossFile = ListGetString(SolverParams,'Harmonic Loss Filename',Found ) 2032 IF( Found ) THEN 2033 OPEN(NEWUNIT=IOUnit, FILE=LossFile) 2034 WRITE(IOUnit,'(A)') '!body_id harmonic(1) harmonic(2) joule' 2035 DO j=1,Model % NumberOfBodies 2036 IF( SUM(BodyLoss(1:3,j)) < TINY( TotalLoss(1) ) ) CYCLE 2037 WRITE(IOUnit,'(I0,T10,3ES17.9)') j, BodyLoss(1:3,j) 2038 END DO 2039 CALL Info('MagnetoDynamicsCalsFields', & 2040 'Harmonic loss for bodies was saved to file: '//TRIM(LossFile),Level=6 ) 2041 CLOSE(IOUnit) 2042 END IF 2043 END IF 2044 2045 DEALLOCATE( BodyLoss ) 2046 END IF 2047 2048 IF (GetLogical(SolverParams,'Show Angular Frequency',Found)) THEN 2049 WRITE(Message,*) 'Angular Frequency: ', Omega 2050 CALL Info( 'MagnetoDynamicsCalcFields', Message ) 2051 CALL ListAddConstReal(Model % Simulation,'res: Angular Frequency', Omega) 2052 END IF 2053 2054 IF(ASSOCIATED(NF)) THEN 2055 CALL NodalTorque(Torque, TorqueGroups) 2056 DO i=1,SIZE(TorqueGroups) 2057 j = TorqueGroups(i) 2058 WRITE (Message,'(A)') 'res: Group '//TRIM(I2S(j))//' torque' 2059 CALL ListAddConstReal(Model % Simulation, TRIM(Message), Torque(i)) 2060 WRITE (Message,'(A,F0.8)') 'Torque Group '//TRIM(I2S(j))//' torque:', Torque(i) 2061 CALL Info( 'MagnetoDynamicsCalcFields', Message) 2062 END DO 2063 2064 CALL NodalTorqueDeprecated(TorqueDeprecated, Found) 2065 IF (Found) THEN 2066 WRITE(Message,*) 'Torque over defined bodies', TorqueDeprecated 2067 CALL Info( 'MagnetoDynamicsCalcFields', Message ) 2068 CALL Warn( 'MagnetoDynamicsCalcFields', 'Keyword "Calculate Torque over body" is deprecated, use Torque Groups instead') 2069 CALL ListAddConstReal(Model % Simulation, 'res: x-axis torque over defined bodies', TorqueDeprecated(1)) 2070 CALL ListAddConstReal(Model % Simulation, 'res: y-axis torque over defined bodies', TorqueDeprecated(2)) 2071 CALL ListAddConstReal(Model % Simulation, 'res: z-axis torque over defined bodies', TorqueDeprecated(3)) 2072 END IF 2073 END IF 2074 2075 ! Flux On Boundary: 2076 !------------------ 2077 2078 CalcFluxLogical = .FALSE. 2079 Flux = 0._dp 2080 Area = 0._dp 2081 AverageFluxDensity = 0._dp 2082 2083 IF (ListGetLogicalAnyBC( Model,'Magnetic Flux Average')) THEN 2084 IF (PiolaVersion) THEN 2085 CALL Warn('MagnetoDynamicsCalcFields', & 2086 'Magnetic Flux Average: The feature is not yet available for Piola transformed basis functions') 2087 ELSE 2088 DO i=1,GetNOFBoundaryElements() 2089 Element => GetBoundaryElement(i) 2090 BC=>GetBC() 2091 IF (.NOT. ASSOCIATED(BC) ) CYCLE 2092 2093 SELECT CASE(GetElementFamily()) 2094 CASE(1) 2095 CYCLE 2096 CASE(2) 2097 k = GetBoundaryEdgeIndex(Element,1); Element => Mesh % Edges(k) 2098 CASE(3,4) 2099 k = GetBoundaryFaceIndex(Element) ; Element => Mesh % Faces(k) 2100 END SELECT 2101 IF (.NOT. ActiveBoundaryElement(Element)) CYCLE 2102 2103 IF (ASSOCIATED(Element % BoundaryInfo % Right)) THEN 2104 BodyId = Element % BoundaryInfo % Right % BodyID 2105 ELSE IF (ASSOCIATED(Element % BoundaryInfo % Left)) THEN 2106 BodyId = Element % BoundaryInfo % Left % BodyID 2107 ELSE 2108 CALL Fatal ('MagnetoDynamicsCalcFields', 'Magnetic Flux Average: Boundary Element has not got a parent element.') 2109 END IF 2110 2111 n = GetElementNOFNodes() 2112 np = n*pSolver % Def_Dofs(GetElementFamily(Element),BodyId,1) 2113 nd = GetElementNOFDOFs(uElement=Element, uSolver=pSolver) 2114 CALL GetVectorLocalSolution(SOL,Pname,uElement=Element,uSolver=pSolver) 2115 2116 CalcFluxLogical = GetLogical( BC, 'Magnetic Flux Average', Found) 2117 IF (Found .AND. CalcFluxLogical) CALL calcAverageFlux(Flux, Area, Element, n, nd, np, SOL, vDOFs) 2118 END DO 2119 Flux(1) = ParallelReduction(Flux(1)) 2120 Flux(2) = ParallelReduction(Flux(2)) 2121 Area = ParallelReduction(Area) 2122 2123 IF( Area < EPSILON( Area ) ) THEN 2124 CALL WARN('MagnetoDynamicsCalcFields', 'Magnetic Flux Average Computation: Area < Epsilon(Area)') 2125 RETURN 2126 END IF 2127 2128 AverageFluxDensity = Flux / Area 2129 2130 WRITE(Message,*) 'Magnetic Flux Average: ', Flux(1) 2131 CALL Info( 'MagnetoDynamicsCalcFields', Message ) 2132 CALL ListAddConstReal( Model % Simulation, 'res: Magnetic Flux Average', Flux(1) ) 2133 2134 IF (vDOFs == 2) THEN 2135 WRITE(Message,*) 'Magnetic Flux im Average: ', Flux(2) 2136 CALL Info( 'MagnetoDynamicsCalcFields', Message ) 2137 CALL ListAddConstReal( Model % Simulation, 'res: Magnetic Flux im Average', Flux(2) ) 2138 END IF 2139 2140 WRITE(Message,*) 'Magnetic Flux Density Average: ', AverageFluxDensity(1) 2141 CALL Info( 'MagnetoDynamicsCalcFields', Message ) 2142 CALL ListAddConstReal( Model % Simulation,'res: Magnetic Flux Density Average', & 2143 AverageFluxDensity(1)) 2144 2145 IF (vDOFs == 2) THEN 2146 WRITE(Message,*) 'Magnetic Flux Density im Average: ', AverageFluxDensity(2) 2147 CALL Info( 'MagnetoDynamicsCalcFields', Message ) 2148 CALL ListAddConstReal( Model % Simulation,'res: Magnetic Flux Density im Average', & 2149 AverageFluxDensity(2)) 2150 END IF 2151 2152 WRITE(Message,*) 'Magnetic Flux Area: ', Area 2153 CALL Info( 'MagnetoDynamicsCalcFields', Message ) 2154 CALL ListAddConstReal( Model % Simulation,'res: Magnetic Flux Area', Area ) 2155 END IF 2156 END IF 2157 2158 ! 2159 ! Some postprocessing of surface currents generated by skin BCs (in time-harmonic cases): 2160 ! 2161 IF (ListCheckPresentAnyBC(Model, 'Layer Electric Conductivity') .AND. vdofs==2) THEN 2162 Power = 0.0_dp 2163 DO i=1,GetNOFBoundaryElements() 2164 Element => GetBoundaryElement(i) 2165 BC => GetBC() 2166 IF (.NOT. ASSOCIATED(BC)) CYCLE 2167 2168 SELECT CASE(GetElementFamily()) 2169 CASE(1) 2170 CYCLE 2171 CASE(2) 2172 k = GetBoundaryEdgeIndex(Element,1) 2173 Element => Mesh % Edges(k) 2174 CASE(3,4) 2175 k = GetBoundaryFaceIndex(Element) 2176 Element => Mesh % Faces(k) 2177 END SELECT 2178 IF (.NOT. ActiveBoundaryElement(Element)) CYCLE 2179 2180 C_ip = GetConstReal(BC, 'Layer Electric Conductivity', Found) 2181 IF (ABS(C_ip) > AEPS) THEN 2182 mu_r = GetConstReal(BC, 'Layer Relative Permeability', Found) 2183 IF (.NOT. Found) mu_r = 1.0_dp 2184 ELSE 2185 CYCLE 2186 END IF 2187 2188 n = GetElementNOFNodes(Element) 2189 nd = GetElementNOFDOFs(uElement=Element, uSolver=pSolver) 2190 np = n ! Note: the scalar potential should be present by default in the time-harmonic case 2191 2192 CALL GetVectorLocalSolution(SOL, Pname, uElement=Element, uSolver=pSolver) 2193 CALL GetElementNodes(Nodes, Element) 2194 2195 IP = GaussPoints(Element, EdgeBasis=.TRUE., PReferenceElement=PiolaVersion, & 2196 EdgeBasisDegree=EdgeBasisDegree) 2197 2198 DO j=1,IP % n 2199 IF ( PiolaVersion ) THEN 2200 stat = EdgeElementInfo(Element, Nodes, IP % U(j), IP % V(j), IP % W(j), & 2201 DetF = DetJ, Basis = Basis, EdgeBasis = WBasis, dBasisdx = dBasisdx, & 2202 BasisDegree = EdgeBasisDegree, ApplyPiolaTransform = .TRUE.) 2203 ELSE 2204 stat = ElementInfo(Element, Nodes, IP % U(j), IP % V(j), IP % W(j), & 2205 detJ, Basis, dBasisdx) 2206 2207 CALL GetEdgeBasis(Element, WBasis, RotWBasis, Basis, dBasisdx) 2208 END IF 2209 2210 val = SQRT(2.0_dp/(C_ip * Omega * 4.0d0 * PI * 1d-7 * mu_r)) ! The layer thickness 2211 Zs = CMPLX(1.0_dp, 1.0_dp, KIND=dp) / (C_ip*val) 2212 2213 E(1,:) = Omega * MATMUL(SOL(2,np+1:nd), WBasis(1:nd-np,:)) - MATMUL(SOL(1,1:np), dBasisdx(1:np,:)) 2214 E(2,:) = -Omega * MATMUL(SOL(1,np+1:nd), WBasis(1:nd-np,:)) - MATMUL(SOL(2,1:np), dBasisdx(1:np,:)) 2215 2216 s = IP % s(j) * detJ 2217 IF( CSymmetry ) THEN 2218 xcoord = SUM( Basis(1:n) * Nodes % x(1:n) ) 2219 s = s * xcoord 2220 END IF 2221 2222 ! Compute the (real) power to maintain the surface current j_S in terms of 2223 ! the surface impedance from the power density P_S = 1/2 Real(1/Zs) E.conjugate(E) 2224 Power = Power + HarmPowerCoeff * REAL(1.0_dp/Zs) * & 2225 (SUM(E(1,:)**2) + SUM(E(2,:)**2)) * s 2226 2227 ! The total power required to maintain the current in the layer when the current density is 2228 ! assumed to be constant through the layer thickness: 2229 !Power = Power + HarmPowerCoeff * C_ip * (SUM(E(1,:)**2) + SUM(E(2,:)**2)) * val * detJ * IP % s(j) 2230 2231 END DO 2232 END DO 2233 Power = ParallelReduction(Power) 2234 2235 WRITE(Message,*) 'Surface current power (the Joule effect): ', Power 2236 CALL Info( 'MagnetoDynamicsCalcFields', Message ) 2237 CALL ListAddConstReal(Model % Simulation, 'res: Surface current power', Power) 2238 2239 END IF 2240 2241 HasThinLines = ListCheckPresentAnyBC(Model,'Thin Line Crossection Area') 2242 IF(HasThinLines) THEN 2243 ThinLinePower = 0._dp 2244 ALLOCATE(ThinLineCrossect(n), ThinLineCond(n)) 2245 Active = GetNOFBoundaryElements() 2246 DO i=1,Active 2247 Element => GetBoundaryElement(i) 2248 BC=>GetBC() 2249 2250 ThinLineCrossect = GetReal( BC, 'Thin Line Crossection Area', Found) 2251 2252 IF (Found) THEN 2253 CALL Info("CalcFields", "Found a Thin Line Element", level=10) 2254 ThinLineCond = GetReal(BC, 'Thin Line Conductivity', Found) 2255 IF (.NOT. Found) CALL Fatal('CalcFields','Thin Line Conductivity not found!') 2256 HasThinLines = .TRUE. 2257 ELSE 2258 CYCLE 2259 END IF 2260 2261 IF (.NOT. ASSOCIATED(BC) ) CYCLE 2262 SELECT CASE(GetElementFamily()) 2263 CASE(1) 2264 CYCLE 2265 CASE(2) 2266 k = GetBoundaryEdgeIndex(Element,1); Element => Mesh % Edges(k) 2267 CASE(3,4) 2268 CYCLE 2269 END SELECT 2270 IF (.NOT. ActiveBoundaryElement(Element)) CYCLE 2271 2272 2273 Model % CurrentElement => Element 2274 nd = GetElementNOFDOFs(Element) 2275 n = GetElementNOFNodes(Element) 2276 CALL GetElementNodes(Nodes, Element) 2277 ! line_tangent = 0._dp 2278 ! line_tangent(1) = Nodes % x(2) - Nodes % x(1) 2279 ! line_tangent(2) = Nodes % y(2) - Nodes % y(1) 2280 ! line_tangent(3) = Nodes % z(2) - Nodes % z(1) 2281 ! line_tangent(:) = line_tangent(:) / SQRT(SUM(line_tangent(:)**2._dp)) 2282 2283 CALL GetVectorLocalSolution(SOL, Pname, uElement=Element, uSolver=pSolver) 2284 IF (Transient) THEN 2285 CALL GetScalarLocalSolution(PSOL,Pname,uSolver=pSolver,Tstep=-1) 2286 PSOL(1:nd)=(SOL(1,1:nd)-PSOL(1:nd))/dt 2287 END IF 2288 2289 ! Numerical integration: 2290 !----------------------- 2291 IP = GaussPoints(Element) 2292 2293 np = n*MAXVAL(Solver % Def_Dofs(GetElementFamily(Element),:,1)) 2294 2295 DO j=1,IP % n 2296 stat = EdgeElementInfo( Element, Nodes, IP % U(j), IP % V(j), & 2297 IP % W(j), DetF = DetJ, Basis = Basis, EdgeBasis = WBasis, & 2298 dBasisdx = dBasisdx, BasisDegree = 1, & 2299 ApplyPiolaTransform = .TRUE.) 2300 2301 C_ip = SUM(Basis(1:n) * ThinLineCond(1:n)) 2302 Area = SUM(Basis(1:n) * ThinLineCrossect(1:n)) 2303 s = detJ*IP % s(j) 2304 2305 IF (vDOFS == 1) THEN 2306 !da/dt part 2307 IF (Transient) THEN 2308 E(1,:) = -MATMUL(PSOL(np+1:nd), Wbasis(1:nd-np,:)) 2309 END IF 2310 2311 !grad V part 2312 E(1,:) = E(1,:)-MATMUL(SOL(1,1:np), dBasisdx(1:np,:)) 2313 2314 ! The Joule heating power per unit volume: J.E = (sigma * E).E 2315 ! Coeff = Area * C_ip * SUM(line_tangent(:) * E(1,:)) ** 2._dp * s 2316 Coeff = Area * C_ip * SUM(E(1,:) ** 2._dp) * s 2317 ELSE 2318 !da/dt part 2319 E(1,:) = Omega*MATMUL(SOL(2,np+1:nd),WBasis(1:nd-np,:)) 2320 !grad V part 2321 E(2,:) = -Omega*MATMUL(SOL(1,np+1:nd),WBasis(1:nd-np,:)) 2322 2323 E(1,:) = E(1,:)-MATMUL(SOL(1,1:np), dBasisdx(1:np,:)) 2324 E(2,:) = E(2,:)-MATMUL(SOL(2,1:np), dBasisdx(1:np,:)) 2325 CALL Warn('CalcFields', 'Power loss not implemented for harmonic case') 2326 Coeff = 0._dp 2327 ! Now Power = J.conjugate(E), with the possible imaginary component neglected. 2328 !Coeff = HarmPowerCoeff * (SUM( MATMUL( REAL(CMat_ip(1:3,1:3)), TRANSPOSE(E(1:1,1:3)) ) * & 2329 ! TRANSPOSE(E(1:1,1:3)) ) * Basis(p) * s - & 2330 ! SUM( MATMUL( AIMAG(CMat_ip(1:3,1:3)), TRANSPOSE(E(2:2,1:3)) ) * & 2331 ! TRANSPOSE(E(1:1,1:3)) ) * Basis(p) * s + & 2332 ! SUM( MATMUL( AIMAG(CMat_ip(1:3,1:3)), TRANSPOSE(E(1:1,1:3)) ) * & 2333 ! TRANSPOSE(E(2:2,1:3)) ) * Basis(p) * s + & 2334 ! SUM( MATMUL( REAL(CMat_ip(1:3,1:3)), TRANSPOSE(E(2:2,1:3)) ) * & 2335 ! TRANSPOSE(E(2:2,1:3)) ) * Basis(p) * s) 2336 END IF 2337 2338 ThinLinePower = ThinLinePower + Coeff 2339 2340 END DO 2341 END DO 2342 2343 ThinLinePower = ParallelReduction(ThinLinePower) 2344 WRITE(Message,*) 'Total thin line power (the Joule effect): ', ThinLinePower 2345 CALL Info( 'MagnetoDynamicsCalcFields', Message ) 2346 CALL ListAddConstReal(Model % Simulation, 'res: thin line power', ThinLinePower) 2347 2348 DEALLOCATE(ThinLineCrossect, ThinLineCond) 2349 END IF 2350 2351 IF( NormIndex > 0 ) THEN 2352 WRITE(Message,*) 'Reverting norm to: ', SaveNorm 2353 CALL Info( 'MagnetoDynamicsCalcFields', Message ) 2354 Solver % Variable % Norm = SaveNorm 2355 END IF 2356 2357CONTAINS 2358 2359!------------------------------------------------------------------- 2360 SUBROUTINE SumElementalVariable(Var, Values, BodyId, uAdditive) 2361!------------------------------------------------------------------- 2362 IMPLICIT NONE 2363 TYPE(Variable_t), POINTER :: Var 2364 REAL(KIND=dp), OPTIONAL, TARGET :: Values(:) 2365 INTEGER, OPTIONAL :: BodyId 2366 LOGICAL, OPTIONAL :: uAdditive 2367 2368 TYPE(Element_t), POINTER :: Element 2369 REAL(KIND=dp), ALLOCATABLE :: NodeSum(:) 2370 INTEGER :: n, j, k, l, nodeind, dgind, bias 2371 LOGICAL, ALLOCATABLE :: AirGapNode(:) 2372 LOGICAL :: Additive 2373 REAL(KIND=dp), POINTER :: ValuesSource(:) 2374 2375 2376 Additive = .FALSE. 2377 IF(PRESENT(uAdditive)) Additive = uAdditive 2378 2379 IF(PRESENT(Values)) THEN 2380 ValuesSource => Values 2381 ELSE 2382 IF( .NOT. ASSOCIATED( Var ) ) RETURN 2383 ValuesSource => Var % Values 2384 END IF 2385 2386 n = Mesh % NumberOFNodes 2387 ALLOCATE(NodeSum(n), AirGapNode(n)) 2388 AirGapNode = .FALSE. 2389 2390 ! Collect nodal sum of DG elements 2391 DO k=1,Var % Dofs 2392 NodeSum = 0.0_dp 2393 2394 ! Collect DG data to nodal vector 2395 DO j=1, Mesh % NumberOfBulkElements 2396 Element => Mesh % Elements(j) 2397 IF(PRESENT(BodyID)) THEN 2398 IF(Element % BodyID /= BodyID) CYCLE 2399 END IF 2400 DO l = 1, Element % TYPE % NumberOfNodes 2401 nodeind = Element % NodeIndexes(l) 2402 dgind = Var % Perm(Element % DGIndexes(l)) 2403 IF( dgind > 0 ) THEN 2404 NodeSum( nodeind ) = NodeSum( nodeind ) + & 2405 ValuesSource(Var % DOFs*( dgind-1)+k ) 2406 END IF 2407 END DO 2408 END DO 2409 2410 ! Sum nodal data to elements 2411 DO j=1, Mesh % NumberOfBulkElements 2412 Element => Mesh % Elements(j) 2413 IF(PRESENT(BodyID)) THEN 2414 IF(Element % BodyID /= BodyID) CYCLE 2415 END IF 2416 DO l=1,Element%TYPE%NumberofNodes 2417 nodeind = Element % NodeIndexes(l) 2418 dgind = Var % Perm(Element % DGIndexes(l)) 2419 IF( dgind > 0 ) THEN 2420 IF( Additive) THEN 2421 Var % Values( var % DOFs*(dgind-1)+k) = NodeSum(nodeind) + & 2422 Var % Values( var % DOFs*(dgind-1)+k) 2423 ELSE 2424 Var % Values( var % DOFs*(dgind-1)+k) = NodeSum(nodeind) 2425 END IF 2426 END IF 2427 END DO 2428 END DO 2429 2430 END DO 2431!------------------------------------------------------------------- 2432 END SUBROUTINE SumElementalVariable 2433!------------------------------------------------------------------- 2434 2435 2436!------------------------------------------------------------------- 2437 SUBROUTINE CalcBoundaryModels( ) 2438!------------------------------------------------------------------- 2439 IMPLICIT NONE 2440!------------------------------------------------------------------- 2441 REAL(KIND=dp) :: GapLength(27), AirGapMu(27) 2442 2443!------------------------------------------------------------------- 2444 LOGICAL :: FirstTime = .TRUE. 2445 REAL(KIND=dp) :: B2, GapLength_ip, LeftCenter(3), & 2446 RightCenter(3), BndCenter(3), LeftNormal(3), RightNormal(3), & 2447 NF_ip_l(27,3), NF_ip_r(27,3), xcoord 2448 TYPE(Element_t), POINTER :: LeftParent, RightParent, BElement 2449 TYPE(Nodes_t), SAVE :: LPNodes, RPNodes 2450 REAL(KIND=dp) :: F(3,3) 2451 INTEGER :: n_lp, n_rp, LeftBodyID, RightBodyID 2452 REAL(KIND=dp), ALLOCATABLE :: LeftFORCE(:,:), RightFORCE(:,:), & 2453 AirGapForce(:,:), ForceValues(:) 2454 INTEGER, ALLOCATABLE :: RightMap(:), LeftMap(:) 2455 REAL(KIND=dp) :: ParentNodalU(n), parentNodalV(n), ParentNodalW(n) 2456 REAL(KIND=dp) :: Normal(3) 2457 REAL(KIND=dp), SAVE :: mu0 = 1.2566370614359173e-6_dp 2458 LOGICAL, ALLOCATABLE :: BodyMask(:) 2459 LOGICAL :: HasLeft, HasRight 2460 2461 n = Mesh % MaxElementDOFs 2462 2463 ALLOCATE(LeftFORCE(n,3), RightForce(n,3), RightMap(n), LeftMap(n), & 2464 AirGapForce(3,Mesh % NumberOfNodes) ) 2465 2466 IF ( FirstTime ) THEN 2467 mu0 = GetConstReal(CurrentModel % Constants, & 2468 'Permeability of Vacuum', Found) 2469 IF(.NOT. Found) mu0 = 1.2566370614359173e-6 2470 END IF 2471 2472 LeftBodyID = -1 2473 RightBodyID = -1 2474 2475 DO i = 1, GetNOFBoundaryElements() 2476 BElement => GetBoundaryElement(i, uSolver=pSolver) 2477 BC => GetBC(BElement) 2478 IF (.NOT. ASSOCIATED(BC) ) CYCLE 2479 2480 n = GetElementNOFNodes(BElement) 2481 GapLength(1:n) = GetReal(BC, 'Air Gap Length', Found) 2482 IF(.NOT. Found) CYCLE 2483 2484 IF( .NOT. ASSOCIATED( BElement % BoundaryInfo ) ) CYCLE 2485 2486 HasLeft = ASSOCIATED(BElement % BoundaryInfo % Left) 2487 HasRight = ASSOCIATED(BElement % BoundaryInfo % Right) 2488 IF( .NOT. (HasLeft .OR. HasRight)) THEN 2489 CALL Warn('MagnetoDynamicsCalcFields', 'Airgap Length given on orphan boundary') 2490 CYCLE 2491 END IF 2492 2493 IF(.NOT. (HasLeft .AND. HasRight)) & 2494 CALL Warn('MagnetoDynamicsCalcFields', 'Onesided airgap force calculation is untested.') 2495 2496 BElement => Mesh % Faces(GetBoundaryFaceIndex(BElement)) 2497 IF(.NOT. ActiveBoundaryElement(BElement, uSolver=pSolver)) CYCLE 2498 2499 LeftBodyID = BElement % BoundaryInfo % Left % BodyID 2500 RightBodyID = BElement % BoundaryInfo % Right % BodyID 2501 IF(LeftBodyID == RightBodyID) THEN 2502 CALL Warn('MagnetoDynamicsCalcFields', 'Airgap in the middle of single body Id') 2503 CYCLE 2504 END IF 2505 2506 IF(HasLeft) LeftParent => BElement % BoundaryInfo % Left 2507 IF(HasRight) RightParent => BElement % BoundaryInfo % Right 2508 2509 IF(HasLeft) n_lp = GetElementNOFNodes(LeftParent) 2510 if(HasRight) n_rp = GetElementNOFNodes(RightParent) 2511 2512 CALL GetElementNodes(Nodes, BElement) 2513 IF(HasLeft) CALL GetElementNodes(LPNodes, LeftParent) 2514 IF(HasRight) CALL GetElementNodes(RPNodes, RightParent) 2515 2516 CALL GetVectorLocalSolution(SOL,Pname,uElement=BElement, uSolver=pSolver) 2517 2518 IF(HasLeft) LeftCenter(1:3) = & 2519 [ SUM(LPNodes % x(1:n_lp)), SUM(LPNodes % y(1:n_lp)), SUM(LPNodes % z(1:n_lp)) ] / n_lp 2520 IF(HasRight) RightCenter(1:3) = & 2521 [ SUM(RPNodes % x(1:n_rp)), SUM(RPNodes % y(1:n_rp)), SUM(RPNodes % z(1:n_rp)) ] / n_rp 2522 BndCenter(1:3) = [ sum(Nodes % x(1:n)), sum(Nodes % y(1:n)), sum(Nodes % z(1:n)) ] / n 2523 2524 np = n*MAXVAL(pSolver % Def_Dofs(GetElementFamily(BElement),:,1)) 2525 nd = GetElementNOFDOFs(uElement=BElement, uSolver=pSolver) 2526 2527 2528 DO k = 1,n 2529 IF(HasLeft) THEN 2530 DO l = 1,n_lp 2531 IF(LeftParent % NodeIndexes(l) == BElement % NodeIndexes(k)) LeftMap(k) = l 2532 END DO 2533 END IF 2534 IF(HasRight) THEN 2535 DO l = 1,n_rp 2536 IF(RightParent % NodeIndexes(l) == BElement % NodeIndexes(k)) RightMap(k) = l 2537 END DO 2538 END IF 2539 END DO 2540 2541 AirGapMu(1:n) = GetReal(BC, 'Air Gap Relative Permeability', Found) 2542 IF(.NOT. Found) AirGapMu(1:n) = 1.0_dp 2543 2544 LeftFORCE = 0.0_dp 2545 RightFORCE = 0.0_dp 2546 2547 IF (SecondOrder) THEN 2548 IP = GaussPoints(BElement, EdgeBasis=dim==3, PReferenceElement=PiolaVersion, EdgeBasisDegree=EdgeBasisDegree) 2549 ELSE 2550 IP = GaussPoints(BElement, EdgeBasis=dim==3, PReferenceElement=PiolaVersion) 2551 END IF 2552 2553 2554 DO j = 1,IP % n 2555 IF ( PiolaVersion ) THEN 2556 stat = EdgeElementInfo( BElement, Nodes, IP % U(j), IP % V(j), IP % W(j), & 2557 F = F, DetF = DetJ, Basis = Basis, EdgeBasis = WBasis, RotBasis = RotWBasis, & 2558 dBasisdx=dBasisdx, BasisDegree = EdgeBasisDegree, ApplyPiolaTransform = .TRUE.) 2559 ELSE 2560 stat = ElementInfo( BElement, Nodes, IP % U(j), IP % V(j), & 2561 IP % W(j), detJ, Basis, dBasisdx ) 2562 2563 CALL GetEdgeBasis(BElement, WBasis, RotWBasis, Basis, dBasisdx) 2564 END IF 2565 2566 R_ip = SUM( Basis(1:n)/(mu0*AirGapMu(1:n)) ) 2567 GapLength_ip = SUM( Basis(1:n)*GapLength(1:n) ) 2568 2569 s = detJ * IP % s(j) 2570 IF ( CSymmetry ) THEN 2571 xcoord = SUM( Basis(1:n) * Nodes % x(1:n) ) 2572 s = s * xcoord 2573 END IF 2574 2575 Normal = NormalVector(BElement, Nodes, IP% U(j), IP % V(j)) 2576 IF(HasLeft) THEN 2577 IF( SUM(normal*(LeftCenter - bndcenter)) >= 0 ) THEN 2578 LeftNormal = -Normal 2579 ELSE 2580 LeftNormal = Normal 2581 END IF 2582 END IF 2583 2584 IF(HasRight) THEN 2585 IF( SUM(normal*(RightCenter - bndcenter)) >= 0 ) THEN 2586 RightNormal = -Normal 2587 ELSE 2588 RightNormal = Normal 2589 END IF 2590 END IF 2591 2592 2593 DO k=1,vDOFs 2594 SELECT CASE(dim) 2595 CASE(2) 2596 ! This has been done with the same sign convention as in MagnetoDynamics2D: 2597 ! ------------------------------------------------------------------------- 2598 IF ( CSymmetry ) THEN 2599 B(k,1) = -SUM( SOL(k,1:nd) * dBasisdx(1:nd,2) ) 2600 B(k,2) = SUM( SOL(k,1:nd) * dBasisdx(1:nd,1) ) & 2601 + SUM( SOL(k,1:nd) * Basis(1:nd) ) / xcoord 2602 B(k,3) = 0._dp 2603 ELSE 2604 B(k,1) = SUM( SOL(k,1:nd) * dBasisdx(1:nd,2) ) 2605 B(k,2) = -SUM( SOL(k,1:nd) * dBasisdx(1:nd,1) ) 2606 B(k,3) = 0._dp 2607 END IF 2608 CASE(3) 2609 B(k,:) = normal*sum( SOL(k,np+1:nd)* RotWBasis(1:nd-np,3) ) 2610 END SELECT 2611 END DO 2612 2613 2614 2615 B2 = SUM(B(1,:)*B(1,:) + B(2,:)*B(2,:)) 2616 IF (ASSOCIATED(NF).OR.ASSOCIATED(EL_NF)) THEN 2617 NF_ip_r = 0._dp 2618 NF_ip_l = 0._dp 2619 DO k=1,n 2620 DO l=1,3 2621 DO m=1,3 2622 IF(HasLeft) NF_ip_l(k,l) = NF_ip_l(k,l) + R_ip*B(1,l)*B(1,m)*(LeftNormal(m)*Basis(k)) 2623 IF(HasRight) NF_ip_r(k,l) = NF_ip_r(k,l) + R_ip*B(1,l)*B(1,m)*(RightNormal(m)*Basis(k)) 2624 END DO 2625 IF(HasLeft) NF_ip_l(k,l) = NF_ip_l(k,l) - 0.5*R_ip*B2*(LeftNormal(l)*Basis(k)) 2626 IF(HasRight) NF_ip_r(k,l) = NF_ip_r(k,l) - 0.5*R_ip*B2*(RightNormal(l)*Basis(k)) 2627 END DO 2628 END DO 2629 END IF 2630 2631 Energy(2) = Energy(2) + GapLength_ip*s*0.5*R_ip*B2 2632 2633 DO p=1,n 2634 IF(HasLeft) LeftFORCE(LeftMap(p), 1:3) = LeftFORCE(LeftMap(p), 1:3) + s*NF_ip_l(p,1:3) 2635 IF(HasRight) RightFORCE(RightMap(p), 1:3) = RightFORCE(RightMap(p), 1:3) + s*NF_ip_r(p,1:3) 2636 END DO 2637 END DO ! Integration points 2638 2639 IF(ElementalFields) THEN 2640 IF(HasLeft) CALL LocalCopy(EL_NF, 3, n_lp, LeftFORCE, 0, UElement=LeftParent, uAdditive=.TRUE.) 2641 IF(HasRight) CALL LocalCopy(EL_NF, 3, n_rp, RightFORCE, 0, UElement=RightParent, uAdditive=.TRUE.) 2642 END IF 2643 END DO ! Boundary elements 2644 2645 DEALLOCATE( LeftFORCE, RightFORCE, RightMap, LeftMap, AirGapForce ) 2646!------------------------------------------------------------------- 2647 END SUBROUTINE CalcBoundaryModels 2648!------------------------------------------------------------------- 2649 2650 2651!------------------------------------------------------------------------------ 2652 SUBROUTINE NodalTorqueDeprecated(T, FoundOne) 2653!------------------------------------------------------------------------------ 2654 IMPLICIT NONE 2655 REAL(KIND=dp), INTENT(OUT) :: T(3) 2656 LOGICAL, INTENT(OUT) :: FoundOne 2657!------------------------------------------------------------------------------ 2658 REAL(KIND=dp) :: P(3), F(3) 2659 TYPE(Element_t), POINTER :: Element 2660 TYPE(Variable_t), POINTER :: CoordVar 2661 LOGICAL :: VisitedNode(Mesh % NumberOfNodes) 2662 INTEGER :: pnodal, nnt, ElemNodeDofs(27), ndofs, globalnode, m, n 2663 LOGICAL :: ONCE=.TRUE., Found 2664 2665 VisitedNode = .FALSE. 2666 FoundOne = .FALSE. 2667 2668 DO n=1,size(Model % bodies) 2669 IF(GetLogical(Model % bodies(n) % Values, 'Calculate Torque over body', FoundOne)) EXIT 2670 END DO 2671 IF(.not. FoundOne) RETURN 2672 T = 0._dp 2673 P = 0._dp 2674 2675 DO pnodal=1,GetNOFActive() 2676 Element => GetActiveElement(pnodal) 2677 IF(GetLogical(GetBodyParams(Element), 'Calculate Torque over body', Found)) THEN 2678 ndofs = GetElementDOFs(ElemNodeDofs) 2679 DO nnt=1,ndofs 2680 globalnode = ElemNodeDofs(nnt) 2681 IF (.NOT. VisitedNode(globalnode)) THEN 2682 F(1) = NF % Values( 3*(NF % Perm((globalnode))-1) + 1) 2683 F(2) = NF % Values( 3*(NF % Perm((globalnode))-1) + 2) 2684 F(3) = NF % Values( 3*(NF % Perm((globalnode))-1) + 3) 2685 P(1) = Mesh % Nodes % x(globalnode) 2686 P(2) = Mesh % Nodes % y(globalnode) 2687 P(3) = Mesh % Nodes % z(globalnode) 2688 T(1) = T(1) + P(2)*F(3)-P(3)*F(2) 2689 T(2) = T(2) + P(3)*F(1)-P(1)*F(3) 2690 T(3) = T(3) + P(1)*F(2)-P(2)*F(1) 2691 VisitedNode(globalnode) = .TRUE. 2692 END IF 2693 END DO ! nnt 2694 END IF 2695 END DO ! pnodal 2696 T(1) = ParallelReduction(T(1)) 2697 T(2) = ParallelReduction(T(2)) 2698 T(3) = ParallelReduction(T(3)) 2699!------------------------------------------------------------------------------ 2700 END SUBROUTINE NodalTorqueDeprecated 2701!------------------------------------------------------------------------------ 2702 2703!------------------------------------------------------------------------------ 2704 SUBROUTINE NodalTorque(T, TorqueGroups) 2705!------------------------------------------------------------------------------ 2706 IMPLICIT NONE 2707 INTEGER, ALLOCATABLE, INTENT(OUT) :: TorqueGroups(:) 2708 REAL(KIND=dp), ALLOCATABLE, INTENT(OUT) :: T(:) 2709!------------------------------------------------------------------------------ 2710! Local variables 2711!------------------------------------------------------------------------------ 2712 2713 REAL(KIND=dp), POINTER :: origins(:,:), omegas(:,:) 2714 TYPE(Element_t), POINTER :: Element 2715 TYPE(Variable_t), POINTER :: CoordVar 2716 TYPE(ValueList_t), POINTER :: BodyParams, SolverParams 2717 TYPE(BodyArray_t), POINTER :: bodies(:) 2718 INTEGER, POINTER :: LocalGroups(:) 2719 REAL(KIND=dp), ALLOCATABLE :: axes(:,:) 2720 2721 LOGICAL, ALLOCATABLE :: VisitedNode(:,:) 2722 INTEGER, ALLOCATABLE :: AllGroups(:) 2723 2724 REAL(KIND=dp) :: origin(3), axisvector(3), P(3), F(3), v1(3), v2(3), nrm 2725 INTEGER :: pnodal, nnt, ElemNodeDofs(27), ndofs, globalnode, ng, ngroups, & 2726 n, maxngroups, m, k, num_origins, num_axes, pivot 2727 LOGICAL :: Found 2728 2729 ! make union of body-wise declared torque groups. \TODO: make this abstract and move to generalutils 2730 bodies => Model % bodies 2731 maxngroups = 0 2732 DO n=1,size(bodies) 2733 LocalGroups => ListGetIntegerArray(bodies(n) % Values, "Torque Groups", Found) 2734 IF (Found) THEN 2735 maxngroups = maxngroups + size(LocalGroups) 2736 END IF 2737 END DO 2738 IF(maxngroups .eq. 0) THEN 2739 ALLOCATE(TorqueGroups(0), T(0)) 2740 RETURN 2741 END IF 2742 2743 ALLOCATE(AllGroups(maxngroups)) 2744 AllGroups = -1 2745 ngroups = 0 2746 DO n=1,size(bodies) 2747 LocalGroups => ListGetIntegerArray(bodies(n) % Values, "Torque Groups", Found) 2748 IF (Found) THEN 2749 AllGroups((ngroups+1):(ngroups+size(LocalGroups))) = LocalGroups(1:size(LocalGroups)) 2750 ngroups = ngroups + size(LocalGroups) 2751 END IF 2752 END DO 2753 call SORT(size(AllGroups), AllGroups) 2754 pivot = AllGroups(1) 2755 k = 1 2756 m = 1 2757 do while(pivot .ne. -1) 2758 AllGroups(k) = pivot 2759 DO n=m,size(AllGroups) 2760 IF (AllGroups(k) .ne. AllGroups(n)) then 2761 pivot = AllGroups(n) 2762 k = k + 1 2763 m = n 2764 exit 2765 end if 2766 pivot = -1 2767 END DO 2768 END DO 2769 ALLOCATE(TorqueGroups(k)) 2770 IF(k .eq. 0) RETURN 2771 2772 TorqueGroups = AllGroups(1:k) 2773 ! done making union 2774 2775 SolverParams => GetSolverParams() 2776 origins => ListGetConstRealArray(SolverParams, "Torque Group Origins", Found) 2777 IF (.NOT. Found) THEN 2778 num_origins = 0 2779 ELSE 2780 num_origins = SIZE(origins,1) 2781 END IF 2782 2783 omegas => ListGetConstRealArray(SolverParams, "Torque Group Axes", Found) 2784 IF (.NOT. Found) THEN 2785 num_axes = 0 2786 ELSE 2787 num_axes = SIZE(omegas,1) 2788 ALLOCATE(axes(num_axes, size(omegas, 2))) 2789 axes = omegas 2790 DO k = 1, num_axes 2791 nrm = sqrt(sum(axes(k,:)*axes(k,:))) 2792 IF (nrm .EQ. 0._dp) THEN 2793 CALL Warn('MagnetoDynamicsCalcFields',& 2794 'Axis for the torque group '//TRIM(I2S(k))//' is a zero vector') 2795 CYCLE 2796 END IF 2797 axes(k,:) = axes(k,:) / nrm 2798 END DO 2799 END IF 2800 2801 ng = size(TorqueGroups,1) 2802 ALLOCATE(T(ng*3)) 2803 ALLOCATE(VisitedNode(Mesh % NumberOfNodes, ng)) 2804 VisitedNode = .FALSE. 2805 T = 0._dp 2806 2807 2808 DO pnodal=1,GetNOFActive() 2809 Element => GetActiveElement(pnodal) 2810 BodyParams => GetBodyParams(Element) 2811 LocalGroups => ListGetIntegerArray(BodyParams, "Torque Groups", Found) 2812 IF(.not. Found) CYCLE 2813 ndofs = GetElementDOFs(ElemNodeDofs) 2814 DO nnt=1,ndofs 2815 globalnode = ElemNodeDofs(nnt) 2816 F(1) = NF % Values( 3*(NF % Perm((globalnode))-1) + 1) 2817 F(2) = NF % Values( 3*(NF % Perm((globalnode))-1) + 2) 2818 F(3) = NF % Values( 3*(NF % Perm((globalnode))-1) + 3) 2819 P(1) = Mesh % Nodes % x(globalnode) 2820 P(2) = Mesh % Nodes % y(globalnode) 2821 P(3) = Mesh % Nodes % z(globalnode) 2822 DO ng=1,size(LocalGroups) 2823 IF (.NOT. VisitedNode(globalnode, LocalGroups(ng))) THEN 2824 VisitedNode(globalnode, LocalGroups(ng)) = .TRUE. 2825 IF (LocalGroups(ng) .gt. num_origins) THEN 2826 origin = 0._dp 2827 ELSE 2828 origin = origins(LocalGroups(ng),1:3) 2829 END IF 2830 IF (LocalGroups(ng) .gt. num_axes) THEN 2831 axisvector = 0._dp 2832 axisvector(3) = 1._dp 2833 ELSE 2834 axisvector = axes(LocalGroups(ng), 1:3) 2835 END IF 2836 v1 = P - origin 2837 v1 = (1 - sum(axisvector*v1))*v1 2838 v2 = CrossProduct(v1,F) 2839 T(LocalGroups(ng)) = T(LocalGroups(ng)) + sum(axisvector*v2) 2840 END IF 2841 END DO 2842 END DO 2843 2844 END DO 2845 DO ng=1,size(TorqueGroups) 2846 T(ng) = ParallelReduction(T(ng)) 2847 END DO 2848 2849!------------------------------------------------------------------------------ 2850 END SUBROUTINE NodalTorque 2851!------------------------------------------------------------------------------ 2852 2853!------------------------------------------------------------------------------ 2854 SUBROUTINE GlobalSol(Var, m, b, dofs ) 2855!------------------------------------------------------------------------------ 2856 IMPLICIT NONE 2857 REAL(KIND=dp), TARGET CONTIG :: b(:,:) 2858 INTEGER :: m, dofs 2859 TYPE(Variable_t), POINTER :: Var 2860!------------------------------------------------------------------------------ 2861 INTEGER :: i 2862!------------------------------------------------------------------------------ 2863 IF(.NOT. ASSOCIATED(var)) RETURN 2864 2865 CALL Info('MagnetoDynamicsCalcFields','Solving for field: '//TRIM(Var % Name),Level=6) 2866 2867 DO i=1,m 2868 dofs = dofs+1 2869 Solver % Matrix % RHS => b(:,dofs) 2870 Solver % Variable % Values=0 2871 Norm = DefaultSolve() 2872 var % Values(i::m) = Solver % Variable % Values 2873 2874 IF( NormIndex == dofs ) SaveNorm = Norm 2875 END DO 2876!------------------------------------------------------------------------------ 2877 END SUBROUTINE GlobalSol 2878!------------------------------------------------------------------------------ 2879 2880 2881!------------------------------------------------------------------------------ 2882 SUBROUTINE LocalSol(Var, m, n, A, b, pivot, dofs ) 2883!------------------------------------------------------------------------------ 2884 IMPLICIT NONE 2885 TYPE(Variable_t), POINTER :: Var 2886 INTEGER :: pivot(:), m,n,dofs 2887 REAL(KIND=dp) :: b(:,:), A(:,:) 2888!------------------------------------------------------------------------------ 2889 INTEGER :: ind(n), i 2890 REAL(KIND=dp) :: x(n) 2891!------------------------------------------------------------------------------ 2892 IF(.NOT. ASSOCIATED(var)) RETURN 2893 2894 ind(1:n) = Var % Perm(Element % DGIndexes(1:n)) 2895 2896 IF( ANY( ind(1:n) <= 0 ) ) RETURN 2897 2898 ind(1:n) = Var % DOFs * (ind(1:n)-1) 2899 2900 DO i=1,m 2901 dofs = dofs+1 2902 x = b(1:n,dofs) 2903 CALL LUSolve(n,MASS,x,pivot) 2904 Var % Values(ind(1:n)+i) = x(1:n) 2905 END DO 2906!------------------------------------------------------------------------------ 2907 END SUBROUTINE LocalSol 2908!------------------------------------------------------------------------------ 2909 2910!------------------------------------------------------------------------------ 2911 SUBROUTINE LocalCopy(Var, m, n, b, bias, UElement, Values, uAdditive) 2912!------------------------------------------------------------------------------ 2913 IMPLICIT NONE 2914 TYPE(Variable_t), POINTER :: Var 2915 INTEGER, INTENT(IN) :: m,n,bias 2916 INTEGER :: dofs 2917 REAL(KIND=dp) :: b(:,:) 2918 TYPE(Element_t), POINTER, OPTIONAL :: UElement 2919 REAL(KIND=dp), OPTIONAL :: Values(:) 2920 LOGICAL, OPTIONAL :: uAdditive 2921!------------------------------------------------------------------------------ 2922 INTEGER :: ind(n), i 2923 LOGICAL :: Additive 2924!------------------------------------------------------------------------------ 2925 IF(.NOT. ASSOCIATED(var)) RETURN 2926 2927 IF(PRESENT(UElement)) THEN 2928 ind(1:n) = Var % Perm(UElement % DGIndexes(1:n)) 2929 ELSE 2930 ind(1:n) = Var % Perm(Element % DGIndexes(1:n)) 2931 END IF 2932 2933 IF( ANY(ind(1:n) == 0 ) ) RETURN 2934 2935 ind(1:n) = Var % Dofs * ( ind(1:n) - 1) 2936 2937 IF(PRESENT(uAdditive)) THEN 2938 Additive = uAdditive 2939 ELSE 2940 Additive = .FALSE. 2941 END IF 2942 2943 dofs = bias 2944 IF(PRESENT(Values)) THEN 2945 DO i=1,m 2946 dofs = dofs+1 2947 IF(Additive) THEN 2948 Values(ind(1:n)+i) = Values(ind(1:n)+i) + b(1:n,dofs) 2949 ELSE 2950 Values(ind(1:n)+i) = b(1:n,dofs) 2951 END IF 2952 END DO 2953 ELSE 2954 DO i=1,m 2955 dofs = dofs+1 2956 IF(Additive) THEN 2957 Var % Values(ind(1:n)+i) = Var % Values(ind(1:n)+i) + b(1:n,dofs) 2958 ELSE 2959 Var % Values(ind(1:n)+i) = b(1:n,dofs) 2960 END IF 2961 END DO 2962 END IF 2963!------------------------------------------------------------------------------ 2964 END SUBROUTINE LocalCopy 2965!------------------------------------------------------------------------------ 2966 2967!------------------------------------------------------------------------------ 2968 SUBROUTINE AddLocalFaceTerms(STIFF,FORCE) 2969!------------------------------------------------------------------------------ 2970 IMPLICIT NONE 2971 REAL(KIND=dp) :: STIFF(:,:), FORCE(:) 2972 2973 TYPE(Element_t),POINTER :: P1,P2,Face,Faces(:) 2974 INTEGER ::t,n,n1,n2,NumberOfFaces,dim 2975 2976 dim = CoordinateSystemDimension() 2977 2978 IF (dim==2) THEN 2979 Faces => Solver % Mesh % Edges 2980 NumberOfFaces = Solver % Mesh % NumberOfEdges 2981 ELSE 2982 Faces => Solver % Mesh % Faces 2983 NumberOfFaces = Solver % Mesh % NumberOfFaces 2984 END IF 2985 2986 DO t=1,NumberOfFaces 2987 Face => Faces(t) 2988 IF ( .NOT. ActiveBoundaryElement(Face) ) CYCLE 2989 2990 P1 => Face % BoundaryInfo % Left 2991 P2 => Face % BoundaryInfo % Right 2992 IF ( ASSOCIATED(P2) .AND. ASSOCIATED(P1) ) THEN 2993 IF(.NOT.ASSOCIATED(GetMaterial(P1),GetMaterial(P2))) CYCLE 2994 2995 n = GetElementNOFNodes(Face) 2996 n1 = GetElementNOFNodes(P1) 2997 n2 = GetElementNOFNodes(P2) 2998 2999 CALL LocalJumps( STIFF,Face,n,P1,n1,P2,n2) 3000 CALL DefaultUpdateEquations( STIFF, FORCE, Face ) 3001 END IF 3002 END DO 3003!------------------------------------------------------------------------------ 3004 END SUBROUTINE AddLocalFaceTerms 3005!------------------------------------------------------------------------------ 3006 3007 3008!------------------------------------------------------------------------------ 3009 SUBROUTINE LocalJumps( STIFF,Face,n,P1,n1,P2,n2) 3010!------------------------------------------------------------------------------ 3011 IMPLICIT NONE 3012 REAL(KIND=dp) :: STIFF(:,:) 3013 INTEGER :: n,n1,n2 3014 TYPE(Element_t), POINTER :: Face, P1, P2 3015!------------------------------------------------------------------------------ 3016 REAL(KIND=dp) :: FaceBasis(n), P1Basis(n1), P2Basis(n2) 3017 REAL(KIND=dp) :: Jump(n1+n2), detJ, U, V, W, S 3018 LOGICAL :: Stat 3019 INTEGER :: i, j, p, q, t, nFace, nParent 3020 TYPE(GaussIntegrationPoints_t) :: IntegStuff 3021 3022 TYPE(Nodes_t) :: FaceNodes, P1Nodes, P2Nodes 3023 SAVE FaceNodes, P1Nodes, P2Nodes 3024!------------------------------------------------------------------------------ 3025 STIFF = 0._dp 3026 3027 CALL GetElementNodes(FaceNodes, Face) 3028 CALL GetElementNodes(P1Nodes, P1) 3029 CALL GetElementNodes(P2Nodes, P2) 3030!------------------------------------------------------------------------------ 3031! Numerical integration over the edge 3032!------------------------------------------------------------------------------ 3033 IntegStuff = GaussPoints( Face ) 3034 3035 DO t=1,IntegStuff % n 3036 U = IntegStuff % u(t) 3037 V = IntegStuff % v(t) 3038 W = IntegStuff % w(t) 3039 S = IntegStuff % s(t) 3040 3041 ! Basis function values & derivatives at the integration point: 3042 !-------------------------------------------------------------- 3043 stat = ElementInfo(Face, FaceNodes, U, V, W, detJ, FaceBasis) 3044 3045 S = S * detJ 3046 3047 ! Find basis functions for the parent elements: 3048 ! --------------------------------------------- 3049 CALL GetParentUVW(Face, n, P1, n1, U, V, W, FaceBasis) 3050 stat = ElementInfo(P1, P1Nodes, U, V, W, detJ, P1Basis) 3051 3052 CALL GetParentUVW(Face, n, P2, n2, U, V, W, FaceBasis) 3053 stat = ElementInfo(P2, P2Nodes, U, V, W, detJ, P2Basis) 3054 3055 ! Integrate jump terms: 3056 ! --------------------- 3057 Jump(1:n1) = P1Basis(1:n1) 3058 Jump(n1+1:n1+n2) = -P2Basis(1:n2) 3059 3060 DO p=1,n1+n2 3061 DO q=1,n1+n2 3062 STIFF(p,q) = STIFF(p,q) + s * Jump(q)*Jump(p) 3063 END DO 3064 END DO 3065 END DO 3066!------------------------------------------------------------------------------ 3067 END SUBROUTINE LocalJumps 3068!------------------------------------------------------------------------------ 3069 3070!------------------------------------------------------------------------------ 3071 SUBROUTINE calcAverageFlux (Flux, Area, Element, n, nd, np, SOL, vDOFs) 3072!------------------------------------------------------------------------------ 3073 IMPLICIT NONE 3074 INTEGER :: n, nd 3075 TYPE(Element_t), POINTER :: Element 3076!------------------------------------------------------------------------------ 3077 REAL(KIND=dp) :: Basis(n),dBasisdx(n,3),DetJ,L(3),Normal(3) 3078 REAL(KIND=dp) :: WBasis(nd,3), RotWBasis(nd,3), B(2, 3), & 3079 SOL(2,32), Flux(2), Area 3080 LOGICAL :: Stat 3081 TYPE(GaussIntegrationPoints_t) :: IP 3082 INTEGER :: j, k, np, vDOFs 3083 3084 TYPE(Nodes_t), SAVE :: Nodes 3085!------------------------------------------------------------------------------ 3086 CALL GetElementNodes( Nodes ) 3087 IP = GaussPoints(Element) 3088 3089 IF( dim == 2 ) THEN 3090 CALL Warn('CalcAverageFlux','Not implemented for 2D problems yet!') 3091 END IF 3092 3093 B=0._dp 3094 3095 DO j=1,IP % n 3096 stat = ElementInfo( Element, Nodes, IP % U(j), IP % V(j), & 3097 IP % W(j), detJ, Basis, dBasisdx ) 3098 CALL GetEdgeBasis(Element, WBasis, RotWBasis, Basis, dBasisdx) 3099 Normal = NormalVector( Element, Nodes, IP % U(j), IP % V(j), .TRUE. ) 3100 3101 s = IP % s(j) * detJ 3102 3103 3104 DO k=1, vDOFs 3105 B(k,:) = MATMUL( SOL(k, np+1:nd), RotWBasis(1:nd-np,:) ) 3106 Flux(k) = Flux(k) + s * SUM(Normal * B(k,:)) 3107 END DO 3108 3109 Area = Area + s 3110 3111 END DO 3112!------------------------------------------------------------------------------ 3113 END SUBROUTINE calcAverageFlux 3114!------------------------------------------------------------------------------ 3115 3116!------------------------------------------------------------------------ 3117END SUBROUTINE MagnetoDynamicsCalcFields 3118!------------------------------------------------------------------------ 3119 3120