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! *  Authors: Juha Ruokolainen
27! *  Email:   Juha.Ruokolainen@csc.fi
28! *  Web:     http://www.csc.fi/elmer
29! *  Address: CSC - IT Center for Science Ltd.
30! *           Keilaranta 14
31! *           02101 Espoo, Finland
32! *
33! *  Original Date: 01 Oct 1996
34! *
35! ****************************************************************************/
36
37!> \ingroup ElmerLib
38!> \{
39
40!-----------------------------------------------------------------------------
41!>  Module defining coordinate systems.
42!-----------------------------------------------------------------------------
43
44MODULE CoordinateSystems
45
46  USE Types
47
48  IMPLICIT NONE
49
50  INTEGER, PARAMETER :: Cartesian = 1
51  INTEGER, PARAMETER :: Cylindric = 2, CylindricSymmetric = 3, AxisSymmetric = 4
52  INTEGER, PARAMETER :: Polar = 5
53  INTEGER :: Coordinates = Cartesian
54
55!----------------------------------------------------------------------------
56CONTAINS
57!----------------------------------------------------------------------------
58
59
60!----------------------------------------------------------------------------
61  FUNCTION CylindricalMetric( r,z,t ) RESULT(metric)
62
63    REAL(KIND=dp) :: r,z,t
64    REAL(KIND=dp), DIMENSION(3,3) :: Metric
65
66    Metric = 0.0d0
67    Metric(1,1) = 1.0d0
68    Metric(2,2) = 1.0d0
69    Metric(3,3) = 1.0d0
70    IF ( r /= 0.0d0 ) Metric(3,3) = 1.0d0 / (r**2)
71  END FUNCTION CylindricalMetric
72!----------------------------------------------------------------------------
73
74
75!----------------------------------------------------------------------------
76  FUNCTION CylindricalSqrtMetric( r,z,t ) RESULT(s)
77
78    REAL(KIND=dp) :: r,z,t,s
79
80    s = r
81  END FUNCTION CylindricalSqrtMetric
82!----------------------------------------------------------------------------
83
84
85!----------------------------------------------------------------------------
86  FUNCTION CylindricalSymbols( r,z,t ) RESULT(symbols)
87    REAL(KIND=dp) :: r,z,t
88
89    REAL(KIND=dp), DIMENSION(3,3,3) :: symbols
90
91    Symbols = 0.0d0
92    Symbols(3,3,1) = -r
93
94    IF ( r /= 0.0d0 ) THEN
95       Symbols(1,3,3) = 1.0d0 / r
96       Symbols(3,1,3) = 1.0d0 / r
97    END IF
98  END FUNCTION CylindricalSymbols
99!----------------------------------------------------------------------------
100
101
102!----------------------------------------------------------------------------
103  FUNCTION CylindricalDerivSymbols( r,z,t ) RESULT(dsymbols)
104    REAL(KIND=dp) :: r,z,t
105
106    REAL(KIND=dp), DIMENSION(3,3,3,3) :: dsymbols
107
108    dSymbols = 0.0d0
109    dSymbols(3,3,1,1) = -1.0d0
110
111    IF ( r /= 0.0 ) THEN
112       dSymbols(1,3,3,1) = -1.0d0 / (r**2)
113       dSymbols(3,1,3,1) = -1.0d0 / (r**2)
114    END IF
115  END FUNCTION CylindricalDerivSymbols
116!----------------------------------------------------------------------------
117
118
119!----------------------------------------------------------------------------
120  FUNCTION PolarMetric( r,p,t ) RESULT(Metric)
121    REAL(KIND=dp) :: r,p,t
122    INTEGER :: i
123    REAL(KIND=dp), DIMENSION(3,3) :: Metric
124
125    Metric = 0.0d0
126    DO i=1,3
127       Metric(i,i) = 1.0d0
128    END DO
129
130    IF ( r /= 0.0d0 ) THEN
131       Metric(2,2) = 1.0d0 / (r**2 * COS(t)**2)
132       IF ( CoordinateSystemDimension() == 3 ) THEN
133          Metric(3,3) = 1.0d0 / r**2
134       END IF
135    END IF
136  END FUNCTION PolarMetric
137!----------------------------------------------------------------------------
138
139
140!----------------------------------------------------------------------------
141  FUNCTION PolarSqrtMetric( r,p,t ) RESULT(s)
142    REAL(KIND=dp) :: r,p,t,s
143
144    IF ( CoordinateSystemDimension() == 2 ) THEN
145       s = SQRT( r**2 * COS(t)**2 )
146    ELSE
147       s = SQRT( r**4 * COS(t)**2 )
148    END IF
149  END FUNCTION PolarSqrtMetric
150!----------------------------------------------------------------------------
151
152
153!----------------------------------------------------------------------------
154  FUNCTION PolarSymbols( r,p,t ) RESULT(symbols)
155    REAL(KIND=dp) :: r,p,t
156
157    REAL(KIND=dp), DIMENSION(3,3,3) :: symbols
158
159    Symbols = 0.0d0
160    Symbols(2,2,1) = -r * COS(t)**2
161    IF ( r /= 0.0d0 ) THEN
162       Symbols(1,2,2) = 1.0d0 / r
163       Symbols(2,1,2) = 1.0d0 / r
164    END IF
165
166    IF ( CoordinateSystemDimension() == 3 ) THEN
167       Symbols(3,3,1) = -r
168       Symbols(2,2,3) = SIN(t)*COS(t)
169
170       Symbols(2,3,2) = -TAN(t)
171       Symbols(3,2,2) = -TAN(t)
172
173       IF ( r /= 0.0d0 ) THEN
174          Symbols(3,1,3) = 1.0d0 / r
175          Symbols(1,3,3) = 1.0d0 / r
176       END IF
177    END IF
178  END FUNCTION PolarSymbols
179!----------------------------------------------------------------------------
180
181
182!----------------------------------------------------------------------------
183  FUNCTION PolarDerivSymbols( r,p,t ) RESULT(dsymbols)
184    REAL(KIND=dp) :: r,p,t
185
186    REAL(KIND=dp), DIMENSION(3,3,3,3) :: dSymbols
187
188    dSymbols = 0.0d0
189    dSymbols(2,2,1,1) = -COS(t)**2
190    IF ( r /= 0.0d0 ) THEN
191       dSymbols(1,2,2,1) = -1.0d0 / r**2
192       dSymbols(2,1,2,1) = -1.0d0 / r**2
193    END IF
194
195    IF ( CoordinateSystemDimension() == 3 ) THEN
196       dSymbols(2,2,1,3) = -2*r*SIN(t)*COS(t)
197       dSymbols(3,3,1,1) = -1
198       dSymbols(2,2,3,3) =  COS(t)**2 - SIN(t)**2
199
200       dSymbols(2,3,2,3) = -1.0d0 / COS(t)**2
201       dSymbols(3,2,2,3) = -1.0d0 / COS(t)**2
202
203       IF ( r /= 0.0d0 ) THEN
204          dSymbols(1,3,3,1) = -1.0d0 / r**2
205          dSymbols(3,1,3,1) = -1.0d0 / r**2
206       END IF
207    END IF
208  END FUNCTION PolarDerivSymbols
209!----------------------------------------------------------------------------
210
211
212
213!----------------------------------------------------------------------------
214  FUNCTION CoordinateSqrtMetric( X,Y,Z ) RESULT( SqrtMetric )
215    REAL(KIND=dp) :: X,Y,Z,SqrtMetric
216
217    IF ( Coordinates == Cartesian ) THEN
218       SqrtMetric = 1.0d0
219    ELSE IF ( Coordinates >= Cylindric .AND. &
220       Coordinates <= AxisSymmetric ) THEN
221       SqrtMetric = CylindricalSqrtMetric( X,Y,Z )
222    ELSE IF ( Coordinates == Polar ) THEN
223       SqrtMetric = PolarSqrtMetric( X,Y,Z )
224    END IF
225
226  END FUNCTION CoordinateSqrtMetric
227!----------------------------------------------------------------------------
228
229
230!----------------------------------------------------------------------------
231  FUNCTION CurrentCoordinateSystem() RESULT(Coords)
232    INTEGER :: Coords
233
234    Coords = Coordinates
235  END FUNCTION CurrentCoordinateSystem
236!----------------------------------------------------------------------------
237
238
239!----------------------------------------------------------------------------
240  SUBROUTINE CoordinateSystemInfo( Metric,SqrtMetric, &
241              Symbols,dSymbols,X,Y,Z )
242
243    REAL(KIND=dp) :: Metric(3,3),SqrtMetric, &
244        Symbols(3,3,3),dSymbols(3,3,3,3)
245
246    INTEGER :: i
247    REAL(KIND=dp) :: X,Y,Z
248
249    IF ( Coordinates == Cartesian ) THEN
250
251       Metric  = 0.0d0
252       DO i=1,3
253          Metric(i,i) = 1.0d0
254       END DO
255
256       SqrtMetric = 1.0d0
257       Symbols    = 0.0d0
258       dSymbols   = 0.0d0
259
260    ELSE IF ( Coordinates >= Cylindric .AND.  &
261         Coordinates <= AxisSymmetric ) THEN
262
263       SqrtMetric = CylindricalSqrtMetric( X,Y,Z )
264       Metric     = CylindricalMetric( X,Y,Z )
265       Symbols    = CylindricalSymbols( X,Y,Z )
266       dSymbols   = CylindricalDerivSymbols( X,Y,Z )
267
268    ELSE IF ( Coordinates == Polar ) THEN
269
270       SqrtMetric = PolarSqrtMetric( X,Y,Z )
271       Metric     = PolarMetric( X,Y,Z )
272       Symbols    = PolarSymbols( X,Y,Z )
273       dSymbols   = PolarDerivSymbols( X,Y,Z )
274
275    END IF
276
277  END SUBROUTINE CoordinateSystemInfo
278!----------------------------------------------------------------------------
279
280
281!----------------------------------------------------------------------------
282  FUNCTION CoordinateSystemDimension() RESULT( dim )
283    INTEGER :: dim ! ,csys
284
285!   csys = CurrentCoordinateSystem()
286    dim  = CurrentModel % Dimension
287  END FUNCTION CoordinateSystemDimension
288!----------------------------------------------------------------------------
289
290!----------------------------------------------------------------------------
291END MODULE CoordinateSystems
292!----------------------------------------------------------------------------
293
294!> \} ElmerLib