MODULE Icosahedron !Tools for creation of global finite element grids by subdivision of the !facets of the icosahedron into quasi-equilateral spherical triangles. ! !by Peter Bird, UCLA; written in BASIC ~1982; translated to Fortran 90 2001.11 !Modified by Zhen Liu, UCLA; on Make_Global_Grid !---------------------------------------------------------------------------- ! ! USER ROUTINES: ! ! SUBROUTINE Make_Global_Grid ! SUBROUTINE Write_Global_Grid ! ! UTILITY ROUTINES (called by the user routines): ! ! RECURSIVE SUBROUTINE Divide ! SUBROUTINE Element_Out ! SUBROUTINE Unitize !---------------------------------------------------------------------------- CONTAINS SUBROUTINE Make_Global_Grid (n_slice) ! only input(!) !--------------------------------------------------------------- !Generate a finite element grid of spherical triangles by subdivision of the icosahedron: !Level 0 (generated here) has 12 vertices, 30 edges, and 20 triangles. !Then, subdivide each face n_slice times, and output triangular elements, ! by using SUBR Divide, which calls itself recursively! !--------------------------------------------------------------- IMPLICIT NONE INTEGER, INTENT(IN) :: n_slice ! subdivision level; 0 or higher !--------------------------------------------------------------- INTEGER :: m_slice ! copy of n_slice, allowing it to be changed (counted-down) INTEGER :: facets_done, i, j, k INTEGER, DIMENSION(3) :: node_number REAL, DIMENSION(3) :: rx, ry, rz DOUBLE PRECISION, PARAMETER :: s = 1.107148719D0 DOUBLE PRECISION :: dot1, dot2, dot3, x1, x2, x3, y1, y2, y3, z1, z2, z3 DOUBLE PRECISION, DIMENSION(12) :: lat, lon DOUBLE PRECISION, DIMENSION(3) :: v1, v2, v3 DOUBLE PRECISION, DIMENSION(3, 12) :: abg !Cartesian (alpha, beta, gamma) coordinates of these vertices. SAVE ! implicit in Digital Fortran; must be requested in Intel Fortran. !generate basic form with a vertex (a 5-fold axis) up; highest symmetry axis lat(1) = 1.570796327D0 lon(1) = 0.0D0 DO i = 2, 6 lat(i) = lat(1) - s lon(i) = (i - 2.0D0) * 1.256637061D0 END DO DO i = 7, 11 lat(i) = -lat(1) + s lon(i) = (i - 7.0D0) * 1.256637061D0 + .628318531D0 END DO lat(12) = -lat(1) lon(12) = 0.0D0 DO i = 1, 12 abg(1, i) = COS(lat(i)) * COS(lon(i)) abg(2, i) = COS(lat(i)) * SIN(lon(i)) abg(3, i) = SIN(lat(i)) END DO !------------------------------------------------------- !create output file for dumping results as they are found: OPEN (UNIT = 729, FILE = "Icosahedron.tmp", FORM = "UNFORMATTED") ! unconditional; overwrites any existing file !------------------------------------------------------- !find all 20 faces and subdivide each into four spherical triangles; ! WRITE (*, "(' Creating global grid by level-',I1,' subdivision of icosahedron facets:')") n_slice ! WRITE (*, *) ! advance, because next WRITE will not m_slice = n_slice facets_done = 0 DO i = 1, 10 DO j = (i + 1), 11 DO k = (j + 1), 12 dot1 = abg(1, i) * abg(1, j) + abg(2, i) * abg(2, j) + abg(3, i) * abg(3, j) dot2 = abg(1, j) * abg(1, k) + abg(2, j) * abg(2, k) + abg(3, j) * abg(3, k) dot3 = abg(1, k) * abg(1, i) + abg(2, k) * abg(2, i) + abg(3, k) * abg(3, i) IF ((dot1 > 0.3D0) .AND. (dot2 > 0.3D0) .AND. (dot3 > 0.3D0)) THEN x1 = abg(1, i) x2 = abg(1, j) x3 = abg(1, k) y1 = abg(2, i) y2 = abg(2, j) y3 = abg(2, k) z1 = abg(3, i) z2 = abg(3, j) z3 = abg(3, k) !Note: Divide will call itself, recursively. !Therefore, all inputs are simple numbers (not vectors), !to go on stack as values, not addresses. !Also, note that output is sent to a temporary file, to avoid multiple copies on stack! CALL Divide(x1, y1, z1, x2, y2, z2, x3, y3, z3, m_slice) ! using copy m_slice = n_slice facets_done = facets_done + 1 ! WRITE (*, "('+',I8,' facets out of 20 divided into elements')") facets_done END IF END DO END DO END DO !----------------------------------------------------------------------- !read binary file, extracting groups of 3 uvecs, and assigning to lists: CLOSE (UNIT = 729) END SUBROUTINE Make_Global_Grid SUBROUTINE Write_Global_Grid (path, & ! [Drive:][\Path\] for output (or null) & n_slice, & & numnod, node_uvec, & & numel, nodes) ! all of these are INTENT(IN); get values ! by calling Make_Global_Grid first IMPLICIT NONE CHARACTER*(*), INTENT(IN) :: path ! [Drive:][\Path\] for output (or null) INTEGER, INTENT(IN) :: n_slice ! subdivision level; 0 or higher INTEGER, INTENT(IN) :: numnod ! number of nodes REAL, DIMENSION(:,:), INTENT(IN) :: node_uvec ! (3, numnod) INTEGER, INTENT(IN) :: numel ! number of elements INTEGER, DIMENSION(:,:), INTENT(IN) :: nodes ! (3, numel) !--------------------------------------------------------------- CHARACTER*1 :: c1 CHARACTER*11 :: grid_file CHARACTER*80 :: grid_pathfile INTEGER :: i REAL, PARAMETER :: Pi = 3.141592654 REAL, PARAMETER :: degrees_per_radian = 57.2957795 REAL :: equat, equat2, lat, lon, phi, theta REAL, DIMENSION(3) :: uvec SAVE ! implicit in Digital Fortran; must be requested in Intel Fortran. !--------------------------------------------------------------- WRITE (c1, "(I1)") n_slice grid_file = "global" // c1 // ".feg" grid_pathfile = TRIM(path) // grid_file WRITE (*, "(' Writing ',A,'...')") TRIM(grid_file) OPEN (UNIT = 730, FILE = grid_pathfile) ! unconditional; overwrites any older file WRITE (730, "(A)") grid_file WRITE (730, "(I6,I6,' 0 1000000 F')") numnod, numnod DO i = 1, numnod uvec(1:3) = node_uvec(1:3, i) !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - !Note: NOT calling Uvec_2_LonLat, because that is in Sphere or Map_Projections/Adobe_Illustrator, ! and in many projects would cause confusing duplication of code. However, steps exactly the same. equat2 = uvec(1)*uvec(1) + uvec(2)*uvec(2) IF (equat2 == 0.) THEN phi = 0. ! actually undefined; default 0. IF (uvec(3) > 0.) THEN theta = 0. ! N pole ELSE theta = Pi ! S pole END IF ELSE equat = SQRT(equat2) theta = ATAN2(equat, uvec(3)) phi = ATAN2(uvec(2), uvec(1)) END IF lat = 90. - degrees_per_radian * ABS(theta) lat = MAX (lat, -90.) lon = degrees_per_radian * phi IF (lon > 180.) lon = lon - 360. IF (lon <= -180.) lon = lon + 360. !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - WRITE (730, "(I6,F9.3,F8.3,' 0.0 0.0 0.0 0.0')") i, lon, lat END DO WRITE (730, *) numel DO i = 1, numel WRITE (730, "(4I6)") i, nodes(1, i), nodes(2, i), nodes(3, i) END DO WRITE (730, "(' 0')") CLOSE (UNIT = 730) WRITE (*, "('+Writing ',A,'...DONE')") grid_file END SUBROUTINE Write_Global_Grid RECURSIVE SUBROUTINE Divide (x1, y1, z1, x2, y2, z2, x3, y3, z3, m_slice) !CALLed only by SUBR Make_Global_Grid (but called multiple times) and by itself (ditto). !Accepts 3 unit vectors of triangle corners as input, subdivides the !triangle m_slice times into 4 triangles, and outputs results to a file. !NOTE that Divide calls itself recursively, so all the arguments should !be simple numbers passed by value, not by address. !To avoid complications, output is dumped into an open binary file (UNIT = 729). !------------------------------------------------------- IMPLICIT NONE DOUBLE PRECISION, INTENT(IN) :: x1, y1, z1, x2, y2, z2, x3, y3, z3 INTEGER, INTENT(IN) :: m_slice !------------------------------------------------------- INTEGER :: l_slice DOUBLE PRECISION :: xe1, xe2, xe3, ye1, ye2, ye3, ze1, ze2, ze3 !------------------------------------------------------- xe1 = 0.5D0 * (x2 + x3) xe2 = 0.5D0 * (x3 + x1) xe3 = 0.5D0 * (x1 + x2) ye1 = 0.5D0 * (y2 + y3) ye2 = 0.5D0 * (y3 + y1) ye3 = 0.5D0 * (y1 + y2) ze1 = 0.5D0 * (z2 + z3) ze2 = 0.5D0 * (z3 + z1) ze3 = 0.5D0 * (z1 + z2) CALL Unitize(xe1, ye1, ze1) CALL Unitize(xe2, ye2, ze2) CALL Unitize(xe3, ye3, ze3) IF (m_slice > 1) THEN l_slice = m_slice - 1 CALL Divide( x1, y1, z1, xe2, ye2, ze2, xe3, ye3, ze3, l_slice) CALL Divide( x2, y2, z2, xe1, ye1, ze1, xe3, ye3, ze3, l_slice) CALL Divide( x3, y3, z3, xe1, ye1, ze1, xe2, ye2, ze2, l_slice) CALL Divide(xe1, ye1, ze1, xe2, ye2, ze2, xe3, ye3, ze3, l_slice) ELSE IF (m_slice == 1) THEN CALL ElementOut( x1, y1, z1, xe2, ye2, ze2, xe3, ye3, ze3) CALL ElementOut( x2, y2, z2, xe1, ye1, ze1, xe3, ye3, ze3) CALL ElementOut( x3, y3, z3, xe1, ye1, ze1, xe2, ye2, ze2) CALL ElementOut(xe1, ye1, ze1, xe2, ye2, ze2, xe3, ye3, ze3) ELSE ! m_slice == 0 CALL ElementOut(x1, y1, z1, x2, y2, z2, x3, y3, z3) END IF END SUBROUTINE Divide SUBROUTINE ElementOut (x1, y1, z1, x2, y2, z2, x3, y3, z3) !CALLed only by SUBR Divide (but called multiple times). !Outputs an element defined by 3 unit radius vectors to its vertices. !The output (binary) file must already be open, as device 729. IMPLICIT NONE DOUBLE PRECISION, INTENT(IN) :: x1, y1, z1, x2, y2, z2, x3, y3, z3 REAL :: rx1, rx2, rx3, ry1, ry2, ry3, rz1, rz2, rz3 SAVE ! implicit in Digital Fortran; must be requested in Intel Fortran. !Convert to single-precision: rx1 = x1; rx2 = x2; rx3 = x3 ry1 = y1; ry2 = y2; ry3 = y3 rz1 = z1; rz2 = z2; rz3 = z3 !Write to binary file: WRITE (729) rx1, ry1, rz1, rx2, ry2, rz2, rx3, ry3, rz3 END SUBROUTINE ElementOut ! ! Unitize already contained in ORBLIB.f90 ! ! ! SUBROUTINE Unitize (x, y, z) ! !Converts any vector in 3-component double-precision form to a unit vector. ! IMPLICIT NONE ! DOUBLE PRECISION, INTENT(INOUT) :: x, y, z ! DOUBLE PRECISION :: r ! r = SQRT(x * x + y * y + z * z) ! IF (r > 0.0D0) THEN ! x = x / r ! y = y / r ! z = z / r ! ELSE ! x = 1.0D0 ! y = 0.0D0 ! z = 0.0D0 ! END IF ! END SUBROUTINE Unitize END MODULE Icosahedron