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! *  Subroutine for transferring material parameters into field values.
27! *
28! ******************************************************************************
29! *
30! *  Authors: Thomas Zwinger, Peter Råback
31! *  Email:   Peter.Raback@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: 20 Nov 2001
38! *
39! *****************************************************************************/
40
41!> \ingroup Solvers
42!> \{
43
44
45!------------------------------------------------------------------------------
46!> Routine for saving material parameters as fields.
47! Written by Thomas Zwinger
48!------------------------------------------------------------------------------
49SUBROUTINE SaveMaterials( Model,Solver,dt,TransientSimulation )
50  USE DefUtils
51  USE Types
52  USE Lists
53  USE Integration
54  USE ElementDescription
55  USE SolverUtils
56
57  IMPLICIT NONE
58!------------------------------------------------------------------------------
59  TYPE(Solver_t), TARGET :: Solver
60  TYPE(Model_t) :: Model
61  REAL(KIND=dp) :: dt
62  LOGICAL :: TransientSimulation
63!------------------------------------------------------------------------------
64! Local variables
65!------------------------------------------------------------------------------
66  TYPE(Solver_t), POINTER :: PointerToSolver
67  TYPE(Mesh_t), POINTER :: Mesh
68  TYPE(Element_t),POINTER :: CurrentElement
69  TYPE(ValueList_t), POINTER :: ValueList, Params
70  TYPE(Variable_t), POINTER :: Var
71  INTEGER :: NoParams, DIM, ParamNo, istat, LocalNodes, &
72       n, j, i, elementNumber, FieldNodes, BodyForceParams
73  INTEGER, POINTER :: NodeIndexes(:), FieldPerm(:)=>NULL()
74  REAL(KIND=dp), POINTER :: Field(:)
75  REAL(KIND=dp), ALLOCATABLE :: LocalParam(:)
76  CHARACTER(LEN=MAX_NAME_LEN) ::  ParamName(99), Name
77  LOGICAL :: GotCoeff, GotIt, GotOper, GotVar
78  INTEGER :: PrevNodes = 0
79
80  SAVE PrevNodes
81
82
83  CALL Info('SaveMaterials','Creating selected material parameters as fields')
84
85  !-------------------------------------------
86  ! Set some pointers and other initialization
87  !-------------------------------------------
88  Params => GetSolverParams()
89  PointerToSolver => Solver
90  Mesh => Solver % Mesh
91
92  DIM = CoordinateSystemDimension()
93  LocalNodes = Model % NumberOfNodes
94
95  IF( PrevNodes > 0 ) THEN
96    IF( LocalNodes /= PrevNodes ) THEN
97      CALL Warn('SaveMaterials','Solver cannot deal with changed mesh: Exiting!')
98      RETURN
99    END IF
100  END IF
101  PrevNodes = LocalNodes
102
103
104  n = Mesh % MaxElementNodes
105  ALLOCATE( LocalParam(n), STAT=istat )
106  IF( istat /= 0 ) CALL Fatal('SaveMaterials','Memory allocation error 1')
107
108  ! Find out how many variables should we saved
109  NoParams = 0
110  GotVar = .TRUE.
111
112  DO WHILE(GotVar)
113    NoParams = NoParams + 1
114    WRITE (Name,'(A,I0)') 'Parameter ',NoParams
115    ParamName(NoParams) = ListGetString( Params, TRIM(Name), GotVar )
116  END DO
117  NoParams = NoParams-1
118
119  IF( NoParams == 0) THEN
120    CALL WARN( 'SaveMaterials', 'No parameters found: No fields will be created')
121    RETURN
122  END IF
123
124  BodyForceParams = ListGetInteger( Params,'Body Force Parameters',GotIt)
125
126  !------------------
127  ! Add new Variables
128  ! -----------------
129  DO ParamNo = 1, NoParams
130    Var => VariableGet( Model % Variables, TRIM(ParamName(ParamNo)), .TRUE.)
131    IF( ASSOCIATED( Var ) ) CYCLE
132
133    ALLOCATE(FieldPerm(LocalNodes),STAT=istat)
134    IF( istat /= 0 ) CALL Fatal('SaveMaterials','Memory allocation error 2')
135
136    FieldPerm = 0
137
138    DO elementNumber=1,Mesh % NumberOfBulkElements + Mesh % NumberOfBoundaryElements
139
140      CurrentElement => Mesh % Elements(elementNumber)
141      !------------------------------------------------------------------
142      ! do nothing, if we are dealing with a halo-element in parallel run
143      !------------------------------------------------------------------
144      !IF (CurrentElement % PartIndex /= Parenv % mype) CYCLE
145
146      n = GetElementNOFNodes(CurrentElement)
147      NodeIndexes => CurrentElement % NodeIndexes
148      Model % CurrentElement => CurrentElement
149
150      IF( ParamNo > BodyForceParams ) THEN
151        ValueList => GetMaterial(CurrentElement)
152      ELSE
153        ValueList => GetBodyForce(CurrentElement)
154      END IF
155      IF(.NOT. ASSOCIATED( ValueList ) ) CYCLE
156      IF( ListCheckPresent(ValueList, TRIM(ParamName(ParamNo))) ) THEN
157        FieldPerm( NodeIndexes ) = 1
158      END IF
159    END DO
160
161    j = 0
162    DO i=1,LocalNodes
163      IF( FieldPerm(i) > 0 ) THEN
164        j = j + 1
165        FieldPerm(i) = j
166      END IF
167    END DO
168
169    IF( j == 0 ) THEN
170      CALL Warn('SaveMaterials',&
171          'Parameter '//TRIM(ParamName(ParamNo))//' not present in any material')
172    END IF
173    WRITE( Message,'(A,I0,A)') 'Parameter > '//TRIM(ParamName(ParamNo))&
174        //' < defined with ',j,' dofs'
175    CALL Info('SaveMaterials',Message)
176
177    IF( j > 0 ) THEN
178      ALLOCATE(Field(j),STAT=istat)
179      IF( istat /= 0 ) CALL Fatal('SaveMaterials','Memory allocation error 3')
180      Field = 0.0_dp
181    END IF
182
183    IF( j > 0 .OR. ParEnv % PEs > 1 ) THEN
184      ! Add this always in parallel since the variable may be active in some partition
185      CALL VariableAdd( Mesh % Variables, Mesh, PointerToSolver, &
186          TRIM(ParamName(ParamNo)), 1, Field, FieldPerm )
187    END IF
188
189    NULLIFY( Field )
190    NULLIFY( FieldPerm )
191  END DO
192
193  !-----------------------------------
194  ! loop all parameters to be exported
195  !-------------------------------------
196  DO ParamNo=1,NoParams
197
198    Var => VariableGet( Model % Variables, TRIM(ParamName(ParamNo)), .TRUE.)
199    IF( .NOT. ASSOCIATED( Var ) ) CYCLE
200
201    Field => Var % Values
202    IF( .NOT. ASSOCIATED( Field ) ) CYCLE
203
204    FieldPerm => Var % Perm
205    IF( .NOT. ASSOCIATED( FieldPerm ) ) CYCLE
206
207    IF( MAXVAL( FieldPerm ) <= 0 ) CYCLE
208
209    DO elementNumber=1,Mesh % NumberOfBulkElements + Mesh % NumberOfBoundaryElements
210
211      CurrentElement => Mesh % Elements(elementNumber)
212      IF (CurrentElement % PartIndex /= Parenv % mype) CYCLE
213
214      n = GetElementNOFNodes(CurrentElement)
215      NodeIndexes => CurrentElement % NodeIndexes
216
217      IF( .NOT. ALL(FieldPerm(NodeIndexes) > 0) ) CYCLE
218
219      Model % CurrentElement => CurrentElement
220
221      IF( ParamNo > BodyForceParams ) THEN
222        ValueList => GetMaterial(CurrentElement)
223      ELSE
224        ValueList => GetBodyForce(CurrentElement)
225      END IF
226      IF(.NOT. ASSOCIATED( ValueList ) ) CYCLE
227      LocalParam(1:n) = ListGetReal(ValueList, TRIM(ParamName(ParamNo)), &
228          n, NodeIndexes, GotIt)
229      IF(.NOT. GotIt) CYCLE
230
231      Field(FieldPerm(NodeIndexes(1:n))) = LocalParam(1:n)
232    END DO
233  END DO
234
235END SUBROUTINE SaveMaterials
236!------------------------------------------------------------------------------
237
238!> \}
239