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