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