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, Peter Råback 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: 10 May 2000 34! * Edited: 9.5.2012 35! * 36! *****************************************************************************/ 37 38!------------------------------------------------------------------------------ 39!> An alternative fork for mesh update solver. The intended use for this 40!> is, for example, in geometry optimization when the normal MeshUpdate is 41!> not always suitable as it may be needed to extent real physical displacements. 42!> This is a dynamically loaded solver with a standard interface. 43!> \ingroup Solvers 44!------------------------------------------------------------------------------ 45 SUBROUTINE MeshSolver( Model,Solver,dt,TransientSimulation ) 46!------------------------------------------------------------------------------ 47 USE DefUtils 48 49 IMPLICIT NONE 50!------------------------------------------------------------------------------ 51 TYPE(Model_t) :: Model 52 TYPE(Solver_t), TARGET :: Solver 53 LOGICAL :: TransientSimulation 54 REAL(KIND=dp) :: dt 55!------------------------------------------------------------------------------ 56! Local variables 57!------------------------------------------------------------------------------ 58 INTEGER :: i,j,k,n,nd,nb,t,ind,STDOFs,LocalNodes,istat,dim,NoActive 59 INTEGER :: VisitedTimes = 0 60 TYPE(Element_t),POINTER :: Element 61 TYPE(ValueList_t),POINTER :: Params, Material, BC 62 REAL(KIND=dp) :: UNorm, maxu, TargetCoeff, val, Relax 63 TYPE(Variable_t), POINTER :: MeshSol, TargetSol, GradSol 64 TYPE(Mesh_t), POINTER :: Mesh 65 REAL(KIND=dp), POINTER :: MeshUpdate(:), PrevMeshUpdate(:) 66 INTEGER, POINTER :: MeshPerm(:), NodeIndexes(:) 67 CHARACTER(LEN=MAX_NAME_LEN) :: TargetFieldName 68 LOGICAL :: AllocationsDone = .FALSE., Isotropic = .TRUE., & 69 GotForceBC, Found, MovingMesh, Cumulative,GotTargetField, & 70 GotTargetSurface, GotGradSol 71 REAL(KIND=dp),ALLOCATABLE:: STIFF(:,:),& 72 LOAD(:,:),FORCE(:), ElasticModulus(:),PoissonRatio(:), & 73 Alpha(:,:), Beta(:), Gamma(:), RefSurface(:) 74 REAL(KIND=dp), POINTER CONTIG :: OrigX(:), OrigY(:), OrigZ(:), & 75 TrueX(:), TrueY(:), TrueZ(:) 76 REAL(KIND=dp) :: at,at0 77 CHARACTER(*), PARAMETER :: Caller = 'NonPhysicalMeshSolve' 78 79 SAVE STIFF, LOAD, FORCE, AllocationsDone, & 80 ElasticModulus, PoissonRatio, & 81 OrigX, OrigY, OrigZ, TrueX, TrueY, TrueZ, PrevMeshUpdate, & 82 VisitedTimes,dim,Alpha,Beta,Gamma,RefSurface 83 84!------------------------------------------------------------------------------ 85 86 CALL Info( Caller, '-------------------------------------', Level=4 ) 87 CALL Info( Caller, 'Nonphysical Mesh Solver:', Level=4 ) 88 CALL Info( Caller, '-------------------------------------', Level=4 ) 89 90 91!------------------------------------------------------------------------------ 92! Get variables needed for solution 93!------------------------------------------------------------------------------ 94 IF ( .NOT. ASSOCIATED( Solver % Matrix ) ) RETURN 95 96 VisitedTimes = VisitedTimes + 1 97 98 dim = CoordinateSystemDimension() 99 MeshSol => Solver % Variable 100 MeshPerm => MeshSol % Perm 101 STDOFs = MeshSol % DOFs 102 MeshUpdate => MeshSol % Values 103 Params => GetSolverParams() 104 Mesh => Solver % Mesh 105 106 LocalNodes = COUNT( MeshPerm > 0 ) 107 108 IF ( LocalNodes <= 0 ) RETURN 109 Cumulative = ListGetLogical( Params,'Cumulative Displacement',Found) 110 111!------------------------------------------------------------------------------ 112! Allocate some permanent storage, this is done first time only 113!------------------------------------------------------------------------------ 114 IF ( .NOT. AllocationsDone ) THEN 115 N = Mesh % MaxElementDOFs 116 117 IF ( AllocationsDone ) THEN 118 DEALLOCATE( ElasticModulus, PoissonRatio, & 119 FORCE, STIFF, Alpha, Beta, Gamma, RefSurface, LOAD, STAT=istat ) 120 END IF 121 122 ALLOCATE( & 123 Alpha(3,N), Beta(N), Gamma(N), RefSurface(N), & 124 ElasticModulus( N ), PoissonRatio( N ), & 125 FORCE( STDOFs*N ), STIFF( STDOFs*N,STDOFs*N ), & 126 LOAD( 4,N ),STAT=istat ) 127 128 IF(.NOT. Cumulative) THEN 129 n = SIZE( Mesh % Nodes % x ) 130 ALLOCATE( OrigX(n), OrigY(n), OrigZ(n) ) 131 OrigX = Mesh % Nodes % x 132 OrigY = Mesh % Nodes % y 133 OrigZ = Mesh % Nodes % z 134 END IF 135 136 IF ( istat /= 0 ) THEN 137 CALL Fatal( Caller, 'Memory allocation error.' ) 138 END IF 139 140 AllocationsDone = .TRUE. 141!------------------------------------------------------------------------------ 142 END IF 143!------------------------------------------------------------------------------ 144 145 TrueX => Mesh % Nodes % x 146 TrueY => Mesh % Nodes % y 147 TrueZ => Mesh % Nodes % z 148 149 IF( .NOT. Cumulative ) THEN 150 Mesh % Nodes % x => OrigX 151 Mesh % Nodes % y => OrigY 152 Mesh % Nodes % z => OrigZ 153 END IF 154 155! This refers to an another mesh deformation routine moving the mesh 156! such that these two must coexist. For generality this is set true 157! by default. 158!-------------------------------------------------------------------- 159 MovingMesh = ListGetLogical( Params,'Moving Mesh', Found ) 160 IF(.NOT. Found) MovingMesh = .TRUE. 161 MovingMesh = MovingMesh .AND. ( VisitedTimes > 1 ) 162 163! This variable is needed to enable that the mesh is changed externally independent of this solver. 164! Then the variations of this solver must always be relative to the previous update. 165!--------------------------------------------------------------------------------- 166 IF( MovingMesh ) THEN 167 IF( VisitedTimes == 2 ) THEN 168 n = SIZE( MeshUpdate ) 169 ALLOCATE( PrevMeshUpdate(n) ) 170 PrevMeshUpdate = 0.0_dp 171 END IF 172 PrevMeshUpdate = MeshUpdate 173 END IF 174 175!------------------------------------------------------------------------------ 176! Do some additional initialization, and go for it 177!------------------------------------------------------------------------------ 178 at = CPUTime() 179 at0 = RealTime() 180 181 CALL DefaultInitialize() 182 CALL StartAdvanceOutput( 'NonphysicalMeshSolve', 'Assembly: ' ) 183!------------------------------------------------------------------------------ 184 185 NoActive = GetNOFActive() 186 187 188 DO t=1,NoActive 189 190 CALL AdvanceOutput( t, NoActive ) 191 192 Element => GetActiveElement(t) 193 nd = GetElementNOFDOFs() 194 nb = GetElementNOFBDOFs() 195 n = GetElementNOFNodes() 196 197 Material => GetMaterial() 198 199 ElasticModulus(1:n) = GetReal( Material,'Mesh Elastic Modulus', Found ) 200 IF ( .NOT. Found ) THEN 201 ElasticModulus(1:n) = GetReal( Material,'Youngs Modulus', Found ) 202 END IF 203 IF ( .NOT. Found ) ElasticModulus(1:n) = 1.0_dp 204 205 PoissonRatio(1:n) = GetReal( Material,'Mesh Poisson Ratio', Found ) 206 IF ( .NOT. Found ) THEN 207 PoissonRatio(1:n) = GetReal( Material,'Poisson Ratio', Found ) 208 END IF 209 IF ( .NOT. Found ) PoissonRatio(1:n) = 0.25_dp 210 211!------------------------------------------------------------------------------ 212! Get element local stiffness & mass matrices 213!------------------------------------------------------------------------------ 214 CALL LocalMatrix( STIFF, FORCE, ElasticModulus, & 215 PoissonRatio, .FALSE., Isotropic, Element, n, nd, nb ) 216 217!------------------------------------------------------------------------------ 218! Update global matrices from local matrices 219!------------------------------------------------------------------------------ 220 CALL DefaultUpdateEquations( STIFF, FORCE ) 221 END DO 222 CALL DefaultFinishBulkAssembly() 223 224 225!------------------------------------------------------------------------------ 226! Nodal displacements may be set by penalty also. This is the target field 227! for the displacements. 228!------------------------------------------------------------------------------ 229 230 TargetFieldName = GetString(Params,'Target Field',GotTargetField ) 231 IF( GotTargetField ) THEN 232 TargetSol => VariableGet( Model % Variables, TargetFieldName ) 233 IF( .NOT. ASSOCIATED(TargetSol)) THEN 234 CALL Fatal('NonphysicalMeshSolve',& 235 'Given > Target Field < does not exist: '//TRIM( TargetFieldName ) ) 236 END IF 237 END IF 238 239 TargetFieldName = GetString(Params,'Target Surface',GotTargetSurface ) 240 GotGradSol = .FALSE. 241 IF( GotTargetSurface ) THEN 242 IF( GotTargetField ) THEN 243 CALL Fatal('NonphysicalMeshSolve','Cannot have > Target Field < and > Target Surface < at same time!') 244 END IF 245 TargetSol => VariableGet( Model % Variables, TargetFieldName ) 246 IF( .NOT. ASSOCIATED(TargetSol)) THEN 247 CALL Fatal('NonphysicalMeshSolve',& 248 'Given > Target Surface < does not exist: '//TRIM( TargetFieldName ) ) 249 END IF 250 251 TargetFieldName = GetString(Params,'Grad Surface',GotGradSol ) 252 IF( GotGradSol ) THEN 253 GradSol => VariableGet( Model % Variables, TargetFieldName ) 254 IF( .NOT. ASSOCIATED(TargetSol)) THEN 255 CALL Fatal('NonphysicalMeshSolve',& 256 'Given > Grad Surface < does not exist: '//TRIM( TargetFieldName ) ) 257 END IF 258 END IF 259 ELSE 260 RefSurface = 0.0_dp 261 END IF 262 263 264!------------------------------------------------------------------------------ 265! Neumann & Newton boundary conditions 266!------------------------------------------------------------------------------ 267 DO t = 1, Mesh % NumberOfBoundaryElements 268 269 Element => GetBoundaryElement(t) 270 IF ( .NOT.ActiveBoundaryElement() ) CYCLE 271 272 BC => GetBC() 273 IF ( .NOT. ASSOCIATED(BC) ) CYCLE 274 275!------------------------------------------------------------------------------ 276! Force in given direction BC: \tau\cdot n = F 277!------------------------------------------------------------------------------ 278 nd = GetElementNOFDOFs() 279 n = GetElementNOFNodes() 280 nb = GetElementNOFBDOFs() 281 282 NodeIndexes => Element % NodeIndexes 283 284 LOAD = 0.0_dp 285 Alpha = 0.0_dp 286 Beta = 0.0_dp 287 Gamma = 0.0_dp 288 RefSurface = 0.0_dp 289 290 Alpha(1,1:n) = GetReal( BC, 'Mesh Coefficient 1', Found ) 291 GotForceBC = Found 292 Alpha(2,1:n) = GetReal( BC, 'Mesh Coefficient 2', Found ) 293 GotForceBC = GotForceBC .OR. Found 294 Alpha(3,1:n) = GetReal( BC, 'Mesh Coefficient 3', Found ) 295 GotForceBC = GotForceBC .OR. Found 296 297 LOAD(1,1:n) = GetReal( BC, 'Mesh Force 1', Found ) 298 GotForceBC = GotForceBC .OR. Found 299 LOAD(2,1:n) = GetReal( BC, 'Mesh Force 2', Found ) 300 GotForceBC = GotForceBC .OR. Found 301 LOAD(3,1:n) = GetReal( BC, 'Mesh Force 3', Found ) 302 GotForceBC = GotForceBC .OR. Found 303 304 Beta(1:n) = GetReal( BC, 'Mesh Normal Force',Found ) 305 GotForceBC = GotForceBC .OR. Found 306 307 Gamma(1:n) = GetReal( BC, 'Mesh Penalty Factor', Found ) 308 GotForceBC = GotForceBC .OR. Found 309 310 IF( GotTargetSurface ) THEN 311 RefSurface(1:n) = GetReal( BC,'Reference Surface',Found) 312 END IF 313 314 IF ( .NOT. GotForceBC ) CYCLE 315 316 CALL MeshBoundary( STIFF,FORCE,LOAD,Alpha,Beta,Gamma,RefSurface,Element,n,nd,nb ) 317 318!------------------------------------------------------------------------------ 319 320 CALL DefaultUpdateEquations( STIFF, FORCE ) 321 END DO 322!------------------------------------------------------------------------------ 323 324 CALL DefaultFinishAssembly() 325 CALL Info( Caller, 'Assembly done', Level=4 ) 326 327!------------------------------------------------------------------------------ 328! Set the nodal displacement for the nodes by using weights, if requested 329!------------------------------------------------------------------------------ 330 CALL NodalDisplacementPenalty() 331 332!------------------------------------------------------------------------------ 333! Dirichlet boundary conditions 334!------------------------------------------------------------------------------ 335 CALL DefaultDirichletBCs() 336 337!------------------------------------------------------------------------------ 338 CALL Info( Caller, 'Set boundaries done', Level=4 ) 339!------------------------------------------------------------------------------ 340! Solve the system and check for convergence 341!------------------------------------------------------------------------------ 342 343 UNorm = DefaultSolve() 344 345 Relax = ListGetCReal( Params,'Nonlinear System Relaxation Factor',Found) 346 IF( Found ) MeshUpdate = Relax * MeshUpdate 347 348 n = SIZE( MeshPerm ) 349 Relax = ListGetCReal( Params,'Mesh Relaxation Factor',Found) 350 IF( .NOT. Found ) Relax = 1.0_dp 351 352 IF( MovingMesh ) THEN 353 DO i=1,n 354 j = MeshPerm(i) 355 IF( j == 0 ) CYCLE 356 IF( dim == 2 ) THEN 357 Truex(i) = Truex(i) + Relax * ( MeshUpdate(2*j-1) - PrevMeshUpdate(2*j-1) ) 358 Truey(i) = Truey(i) + Relax * ( MeshUpdate(2*j) - PrevMeshUpdate(2*j) ) 359 ELSE 360 Truex(i) = Truex(i) + Relax * ( MeshUpdate(3*j-2) - PrevMeshUpdate(3*j-2) ) 361 Truey(i) = Truey(i) + Relax * ( MeshUpdate(3*j-1) - PrevMeshUpdate(3*j-1) ) 362 Truez(i) = Truez(i) + Relax * ( MeshUpdate(3*j) - PrevMeshUpdate(3*j) ) 363 END IF 364 END DO 365 ELSE 366 DO i=1,n 367 j = MeshPerm(i) 368 IF( j == 0 ) CYCLE 369 IF( dim == 2 ) THEN 370 Truex(i) = Truex(i) + Relax * MeshUpdate(2*j-1) 371 Truey(i) = Truey(i) + Relax * MeshUpdate(2*j) 372 ELSE 373 Truex(i) = Truex(i) + Relax * MeshUpdate(3*j-2) 374 Truey(i) = Truey(i) + Relax * MeshUpdate(3*j-1) 375 Truez(i) = Truez(i) + Relax * MeshUpdate(3*j) 376 END IF 377 END DO 378 END IF 379 380 IF(.NOT. Cumulative) THEN 381 Mesh % Nodes % x => TrueX 382 Mesh % Nodes % y => TrueY 383 Mesh % Nodes % z => TrueZ 384 END IF 385 386 387 388 CONTAINS 389 390!------------------------------------------------------------------------------ 391 SUBROUTINE LocalMatrix( STIFF,FORCE,NodalYoung, NodalPoisson, & 392 PlaneStress, Isotropic, Element,n, nd, nb ) 393!------------------------------------------------------------------------------ 394 IMPLICIT NONE 395 396 REAL(KIND=dp) :: NodalPoisson(:), NodalYoung(:) 397 REAL(KIND=dp), TARGET :: STIFF(:,:), FORCE(:) 398 INTEGER :: n,nd,nb 399 TYPE(Element_t) :: Element 400 LOGICAL :: PlaneStress, Isotropic 401!------------------------------------------------------------------------------ 402! 403 REAL(KIND=dp) :: Basis(nd) 404 REAL(KIND=dp) :: dBasisdx(nd,3),detJ 405 406 REAL(KIND=dp) :: NodalLame1(n),NodalLame2(n),Lame1,Lame2, & 407 Poisson, Young, Coeff 408 409 REAL(KIND=dp), POINTER :: A(:,:) 410 REAL(KIND=dp) :: s,u,v,w 411 INTEGER :: i,j,k,p,q,t,dim 412 LOGICAL :: stat 413 TYPE(GaussIntegrationPoints_t), TARGET :: IntegStuff 414 415 TYPE(Nodes_t) :: Nodes 416 SAVE Nodes 417!------------------------------------------------------------------------------ 418 419 CALL GetElementNodes( Nodes ) 420 dim = CoordinateSystemDimension() 421 422 Coeff = ListGetConstReal( Params,'Mass Coefficient',Stat) 423 424 IF ( PlaneStress ) THEN 425 NodalLame1(1:n) = NodalYoung(1:n) * NodalPoisson(1:n) / & 426 ((1.0d0 - NodalPoisson(1:n)**2)) 427 ELSE 428 NodalLame1(1:n) = NodalYoung(1:n) * NodalPoisson(1:n) / & 429 ((1.0d0 + NodalPoisson(1:n)) * (1.0d0 - 2.0d0*NodalPoisson(1:n))) 430 END IF 431 432 NodalLame2(1:n) = NodalYoung(1:n) / (2* (1.0d0 + NodalPoisson(1:n))) 433 434 STIFF = 0.0d0 435 FORCE = 0.0d0 436 437 ! Integration stuff: 438 ! ------------------ 439 IntegStuff = GaussPoints( Element ) 440 441 ! Now we start integrating: 442 ! ------------------------- 443 DO t=1,IntegStuff % n 444 u = IntegStuff % u(t) 445 v = IntegStuff % v(t) 446 w = IntegStuff % w(t) 447 448 ! Basis function values & derivatives at the integration point: 449 !-------------------------------------------------------------- 450 stat = ElementInfo( Element,Nodes, u, v, w, detJ, & 451 Basis, dBasisdx ) 452 453 s = detJ * IntegStuff % s(t) 454 455 ! Lame parameters at the integration point: 456 ! ----------------------------------------- 457 Lame1 = SUM( NodalLame1(1:n)*Basis(1:n) ) 458 Lame2 = SUM( NodalLame2(1:n)*Basis(1:n) ) 459 460 461 ! Loop over basis functions (of both unknowns and weights): 462 ! --------------------------------------------------------- 463 DO p=1,nd 464 DO q=p,nd 465 A => STIFF( dim*(p-1)+1:dim*p,dim*(q-1)+1:dim*q ) 466 DO i=1,dim 467 DO j = 1,dim 468 A(i,j) = A(i,j) + s * Lame1 * dBasisdx(q,j) * dBasisdx(p,i) 469 A(i,i) = A(i,i) + s * Lame2 * dBasisdx(q,j) * dBasisdx(p,j) 470 A(i,j) = A(i,j) + s * Lame2 * dBasisdx(q,i) * dBasisdx(p,j) 471 END DO 472 A(i,i) = A(i,i) + s * Coeff * dBasisdx(q,i) * dBasisdx(p,i) 473 END DO 474 END DO 475 END DO 476 END DO 477 478 ! Assign the symmetric block: 479 ! --------------------------- 480 DO p=1,dim*nd 481 DO q=1,p-1 482 STIFF(p,q)=STIFF(q,p) 483 END DO 484 END DO 485 486 IF ( nb == 0 ) THEN 487 DO p=nd-Element % BDOFs+1,nd 488 DO i=1,dim 489 j = (p-1)*dim + i 490 STIFF( j,: ) = 0.0d0 491 STIFF( :,j ) = 0.0d0 492 STIFF( j,j ) = 1.0d0 493 FORCE( j ) = 0.0d0 494 END DO 495 END DO 496 END IF 497!------------------------------------------------------------------------------ 498 END SUBROUTINE LocalMatrix 499!------------------------------------------------------------------------------ 500 501!------------------------------------------------------------------------------ 502 SUBROUTINE MeshBoundary( STIFF,FORCE,LOAD,NodalAlpha,NodalBeta,NodalGamma,& 503 NodalRefSurface,Element,n,nd,nb ) 504!------------------------------------------------------------------------------ 505 REAL(KIND=dp) :: STIFF(:,:),FORCE(:) 506 REAL(KIND=dp) :: NodalAlpha(:,:),NodalBeta(:),NodalGamma(:),NodalRefSurface(:),LOAD(:,:) 507 TYPE(Element_t),POINTER :: Element 508 INTEGER :: n,nd,nb 509!------------------------------------------------------------------------------ 510 REAL(KIND=dp) :: Basis(nd), dBasisdx(nd,3),detJ 511 REAL(KIND=dp) :: u,v,w,s 512 REAL(KIND=dp) :: Alpha(3),Beta,Gamma,Normal(3),LoadAtIP(3),ExtDisp(3),RefSurface, & 513 GradSurface(3), AbsGradSurface2, Surface 514 REAL(KIND=dp), ALLOCATABLE :: ParentBasis(:), dParentBasisdx(:,:) 515 REAL(KIND=dp) :: ParentdetJ 516 517 INTEGER :: i,t,q,p,dim,np 518 LOGICAL :: stat, AllocationsDone 519 TYPE(Element_t),POINTER :: Parent 520 521 TYPE(GaussIntegrationPoints_t), TARGET :: IntegStuff 522 523 TYPE(Nodes_t) :: Nodes, ParentNodes 524 SAVE Nodes, ParentNodes, AllocationsDone, ParentBasis, dParentBasisdx 525!------------------------------------------------------------------------------ 526 527 dim = Element % TYPE % DIMENSION + 1 528 CALL GetElementNodes( Nodes ) 529 530 FORCE = 0.0D0 531 STIFF = 0.0D0 532 533 IntegStuff = GaussPoints( Element ) 534 535 IF( GotTargetSurface .AND. .NOT. GotGradSol ) THEN 536 IF( .NOT. AllocationsDone ) THEN 537 np = Mesh % MaxElementNodes 538 ALLOCATE( ParentBasis(np), dParentBasisdx(np,3)) 539 AllocationsDone = .TRUE. 540 END IF 541 542 Parent => Element % BoundaryInfo % Left 543 IF ( .NOT. ASSOCIATED(Parent) ) & 544 Parent => Element % BoundaryInfo % Right 545 IF(.NOT.ASSOCIATED(Parent)) RETURN 546 547 np = GetElementNOFDOFs(Parent) 548 CALL GetElementNodes( ParentNodes, Parent ) 549 END IF 550 551 552 DO t=1,IntegStuff % n 553 554 u = IntegStuff % u(t) 555 v = IntegStuff % v(t) 556 w = IntegStuff % w(t) 557 558!------------------------------------------------------------------------------ 559! Basis function values & derivatives at the integration point 560!------------------------------------------------------------------------------ 561 stat = ElementInfo( Element, Nodes, u, v, w, detJ, & 562 Basis, dBasisdx ) 563 564 s = detJ * IntegStuff % s(t) 565!------------------------------------------------------------------------------ 566 LoadAtIP = 0.0D0 567 DO i=1,dim 568 LoadAtIP(i) = SUM( LOAD(i,1:n)*Basis ) 569 Alpha(i) = SUM( NodalAlpha(i,1:n)*Basis ) 570 END DO 571 572 Normal = NormalVector( Element,Nodes,u,v,.TRUE. ) 573 LoadAtIP = LoadAtIP + SUM( NodalBeta(1:n) * Basis ) * Normal 574 575 IF( GotTargetSurface .OR. GotTargetField ) THEN 576 IF( GotTargetField ) THEN 577 DO i=1,dim 578 ExtDisp(i) = SUM( Basis(1:n) * TargetSol % Values( & 579 dim * ( TargetSol % Perm( Element % NodeIndexes(1:n) ) - 1) + i ) ) 580 END DO 581 ELSE 582 RefSurface = SUM( Basis(1:n) * NodalRefSurface(1:n)) 583 Surface = SUM( Basis(1:n) * TargetSol % Values( Element % NodeIndexes(1:n) ) ) 584 585 IF( GotGradSol ) THEN 586 DO i=1,dim 587 GradSurface(i) = & 588 SUM( Basis(1:n) * GradSol % Values( dim*(Element % NodeIndexes(1:n)-1)+i ) ) 589 END DO 590 ELSE 591 CALL GetParentUVW( Element,n,Parent,np,U,V,W,Basis ) 592 stat = ElementInfo( Parent,ParentNodes,U,V,W,ParentdetJ, & 593 ParentBasis,dParentBasisdx ) 594 DO i=1,dim 595 GradSurface(i) = & 596 SUM( dParentBasisdx(1:np,i) * TargetSol % Values( Parent % NodeIndexes(1:np) ) ) 597 END DO 598 END IF 599 ! A suggested displacement to the direction of the normal. 600 ! The absolute gradient is used both in the normalization and 601 ! when making the unit vector and hence the SQRT operator is not 602 ! performed on purpose. 603 AbsGradSurface2 = SUM( GradSurface(1:dim)**2) 604 DO i=1,dim 605 ExtDisp(i) = (RefSurface-Surface) * GradSurface(i) / AbsGradSurface2 606 END DO 607 END IF 608 Gamma = SUM( NodalGamma(1:n) * Basis ) 609 LoadAtIP = LoadAtIP + ExtDisp * Gamma 610 ELSE 611 Gamma = 0.0_dp 612 END IF 613 614 615 DO p=1,nd 616 DO q=1,nd 617 DO i=1,dim 618 STIFF((p-1)*dim+i,(q-1)*dim+i) = & 619 STIFF((p-1)*dim+i,(q-1)*dim+i) + & 620 s * ( Alpha(i) + Gamma ) * Basis(q) * Basis(p) 621 END DO 622 END DO 623 END DO 624 625 DO q=1,nd 626 DO i=1,dim 627 FORCE((q-1)*dim+i) = FORCE((q-1)*dim+i) + & 628 s * Basis(q) * LoadAtIP(i) 629 END DO 630 END DO 631 632 END DO 633 634 IF ( nb == 0 ) THEN 635 DO p=nd-Element % BDOFs+1,nd 636 DO i=1,dim 637 j = (p-1)*dim + i 638 STIFF( j,: ) = 0.0d0 639 STIFF( :,j ) = 0.0d0 640 STIFF( j,j ) = 1.0d0 641 FORCE( j ) = 0.0d0 642 END DO 643 END DO 644 END IF 645 646!------------------------------------------------------------------------------ 647 END SUBROUTINE MeshBoundary 648!------------------------------------------------------------------------------ 649 650!----------------------------------------------------------- 651! Set the nodal coordinates by penalty 652!------------------------------------------------------------ 653 SUBROUTINE NodalDisplacementPenalty() 654 655 INTEGER :: i,j,k,j2 656 REAL(KIND=dp) :: ExtDisp(3),GradSurface(3),AbsGradSurface2,Surface,& 657 RefSurface,TargetCoeff,Coeff 658 LOGICAL :: LocalPenalty, GotWeightSol,Visited=.FALSE. 659 INTEGER, POINTER :: PenaltyPerm(:),WeightPerm(:) 660 TYPE(Variable_t), POINTER :: WeightSol 661 662 SAVE Visited 663 664 IF(.NOT. (GotTargetSurface .OR. GotTargetField) ) RETURN 665 666 TargetCoeff = GetCReal(Params,'Nodal Penalty Factor',Found) 667 IF( ABS( TargetCoeff ) < TINY(TargetCoeff) ) RETURN 668 669 GotWeightSol = .FALSE. 670 IF( GetLogical( Params,'Use Boundary Weights',Found) ) THEN 671 IF( .NOT. Visited ) THEN 672 CALL CalculateNodalWeights( Solver,.TRUE.,VarName = 'Boundary Weights') 673 WeightSol => VariableGet( Model % Variables,'Boundary Weights') 674 WeightSol % Output = .FALSE. 675 ELSE 676 WeightSol => VariableGet( Model % Variables,'Boundary Weights') 677 END IF 678 GotWeightSol = ASSOCIATED(WeightSol) 679 IF( GotWeightSol ) THEN 680 CALL Info(Caller,'Using Boundary Weights for setting nodal penalty',Level=8) 681 ELSE 682 CALL Fatal(Caller,'Boundary Weights not present!') 683 END IF 684 END IF 685 686 LocalPenalty = ListCheckPresentAnyBC( Model,'Apply Nodal Penalty') 687 IF(.NOT. LocalPenalty ) THEN 688 LocalPenalty = ListCheckPresentAnyBodyForce( Model,'Apply Nodal Penalty') 689 END IF 690 691 IF( LocalPenalty ) THEN 692 CALL Info(Caller,'Setting permutation for local penalties') 693 ALLOCATE( PenaltyPerm( Mesh % NumberOfNodes ) ) 694 CALL MakePermUsingMask( Model, Solver, Mesh,'Apply Nodal Penalty', .FALSE., PenaltyPerm, i ) 695 ! PRINT *,'Number of penalty nodes',i 696 END IF 697 698 IF( .NOT. GotTargetField ) THEN 699 RefSurface = ListGetCReal( Params,'Reference Surface') 700 END IF 701 702 DO i=1,Mesh % NumberOfNodes 703 704 j = MeshPerm(i) 705 IF(j==0) CYCLE 706 707 j2 = TargetSol % Perm(i) 708 IF(j2==0) CYCLE 709 710 IF(LocalPenalty) THEN 711 IF( PenaltyPerm(i) == 0 ) CYCLE 712 END IF 713 714 IF( GotTargetField ) THEN 715 DO k=1,dim 716 ExtDisp(k) = TargetSol % Values( dim * (j2-1) + k ) 717 END DO 718 ELSE 719 Surface = TargetSol % Values( j2 ) 720 IF( GotGradSol ) THEN 721 j2 = GradSol % Perm(i) 722 DO k=1,dim 723 GradSurface(k) = GradSol % Values( dim * (j2-1) + k ) 724 END DO 725 ELSE 726 CALL Fatal('NonphysicalMeshSolve','You should provide GradSol!') 727 END IF 728 AbsGradSurface2 = SUM( GradSurface(1:dim)**2) 729 DO k=1,dim 730 ExtDisp(k) = (RefSurface-Surface) * GradSurface(k) / AbsGradSurface2 731 END DO 732 END IF 733 734 DO k=1,dim 735 ind = 3 * (j-1) + k 736 Coeff = TargetCoeff 737 IF( GotWeightSol ) Coeff = Coeff * WeightSol % Values( WeightSol % Perm(i) ) 738 Solver % Matrix % rhs(ind) = Solver % Matrix % rhs(ind) + TargetCoeff * ExtDisp(k) 739 CALL CRS_AddToMatrixElement( Solver % Matrix,ind,ind,TargetCoeff ) 740 END DO 741 END DO 742 743 IF( LocalPenalty ) DEALLOCATE( PenaltyPerm ) 744 745 Visited = .TRUE. 746 747 END SUBROUTINE NodalDisplacementPenalty 748 749 750 751!------------------------------------------------------------------------------ 752END SUBROUTINE MeshSolver 753!------------------------------------------------------------------------------ 754