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