Next: I.1.4 Finalize_Multi_Mesh Procedure Up: I.1 Multi_Mesh Class Code Previous: I.1.2 Initialize_Uniform_Multi_Mesh Procedure
The main documentation of the Initialize_Orthogonal_Multi_Mesh Procedure contains additional explanation of this code listing.
subroutine Initialize_Orthogonal_MMesh (Mesh, NDimensions, &
Coordinates_Nodes_X, Coordinates_Nodes_Y, Coordinates_Nodes_Z, &
Mesh_Name, status)
! Input variables.
type(integer), intent(in) :: NDimensions ! Number of Dimensions.
type(real,1) :: Coordinates_Nodes_X ! X-coordinates of the nodes.
type(real,1), optional :: Coordinates_Nodes_Y ! Y-coordinates of the nodes.
type(real,1), optional :: Coordinates_Nodes_Z ! Z-coordinates of the nodes.
type(character,*), intent(in), optional :: Mesh_Name ! Mesh name.
! Output variables.
! Multi_Mesh to be initialized.
type(Multi_Mesh_type), intent(inout) :: Mesh
type(Status_type), intent(out), optional :: status ! Exit status.
! Internal variables.
! Total number of cells in the X-, Y-, and Z-directions.
type(integer) :: NCells_X_total, NCells_Y_total, NCells_Z_total
type(character,name_length) :: Shape ! Cell shape.
type(character,name_length) :: Geometry ! Cell geometry (Cartesian).
type(character,name_length) :: Uniformity ! Set to "Uniform".
type(character,name_length) :: Orthogonality ! Set to "Orthogonal".
type(character,name_length) :: Structure ! Set to "Structured".
type(logical) :: AMR ! Set to false.
type(Status_type), dimension(20) :: allocate_status ! Allocation Status.
type(Status_type) :: consolidated_status ! Consolidated Status.
type(integer) :: NPEs_X ! Number of PEs in X.
type(integer) :: NPEs_Y ! Number of PEs in Y.
type(integer) :: NPEs_Z ! Number of PEs in Z.
type(real,1) :: Lengths ! Physical extent of the
! domain in each direction.
! Structure length vectors, which give numbers for all PEs.
type(integer,1) :: NCells_Vector ! Number of cells.
type(integer,1) :: NCells_X_Vector ! Number of cells in the X direction.
type(integer,1) :: NCells_Y_Vector ! Number of cells in the Y direction.
type(integer,1) :: NCells_Z_Vector ! Number of cells in the Z direction.
type(integer,1) :: NFaces_Vector ! Number of faces.
type(integer,1) :: NNodes_Vector ! Number of nodes.
type(integer,1) :: NNodes_X_Vector ! Number of nodes in the X direction.
type(integer,1) :: NNodes_Y_Vector ! Number of nodes in the Y direction.
type(integer,1) :: NNodes_Z_Vector ! Number of nodes in the Z direction.
! Location of this_PE in the PE-mesh.
type(integer) :: this_PE_X, this_PE_Y, this_PE_Z
type(integer) :: node, pe_x, pe_y, pe_z ! Loop parameters.
! Offsets - starting points for this PE.
type(integer) :: Offset_PE_X, Offset_PE_Y, Offset_PE_Z ! 1st node indices.
! Mesh coordinates and indices.
! The coordinates of the nodes on this PE.
type(real,2) :: Coordinates_Nodes_PE
! The nodes for the cells on this PE.
type(integer,2) :: Nodes_of_Cells_PE
! The cells for the cells (that is, across each face) on this PE.
type(integer,2) :: Cells_of_Cells_PE
! Face Flags for Structured Meshes.
type(integer,2) :: Flag_Faces_of_Cells
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! Verify requirements.
! Mesh type info is valid.
VERIFY(NDimensions .InInterval. (/1, 3/),5)
! Allocations and initializations.
call Initialize (allocate_status)
call Initialize (consolidated_status)
! Set NCells_D_total and physical dimensions of the mesh - Lengths.
call Initialize (Lengths, NDimensions, allocate_status(1))
select case (NDimensions)
case (1)
NCells_X_total = SIZE(Coordinates_Nodes_X) - 1
NCells_Y_total = 1
NCells_Z_total = 1
Lengths(1) = Coordinates_Nodes_X(SIZE(Coordinates_Nodes_X)) &
- Coordinates_Nodes_X(1)
case (2)
NCells_X_total = SIZE(Coordinates_Nodes_X) - 1
NCells_Y_total = SIZE(Coordinates_Nodes_Y) - 1
NCells_Z_total = 1
Lengths(1) = Coordinates_Nodes_X(SIZE(Coordinates_Nodes_X)) &
- Coordinates_Nodes_X(1)
Lengths(2) = Coordinates_Nodes_Y(SIZE(Coordinates_Nodes_Y)) &
- Coordinates_Nodes_Y(1)
case (3)
NCells_X_total = SIZE(Coordinates_Nodes_X) - 1
NCells_Y_total = SIZE(Coordinates_Nodes_Y) - 1
NCells_Z_total = SIZE(Coordinates_Nodes_Z) - 1
Lengths(1) = Coordinates_Nodes_X(SIZE(Coordinates_Nodes_X)) &
- Coordinates_Nodes_X(1)
Lengths(2) = Coordinates_Nodes_Y(SIZE(Coordinates_Nodes_Y)) &
- Coordinates_Nodes_Y(1)
Lengths(3) = Coordinates_Nodes_Z(SIZE(Coordinates_Nodes_Z)) &
- Coordinates_Nodes_Z(1)
end select
! Generate the connectivity.
call Gen_StructureMesh_Connectivity (NDimensions, Lengths, &
NCells_X_total, NCells_Y_total, NCells_Z_total, Shape, Structure, AMR, &
NPEs_X, NPEs_Y, NPEs_Z, this_PE_X, this_PE_Y, this_PE_Z, &
NNodes_Vector, NCells_Vector, NFaces_Vector, &
NCells_X_Vector, NCells_Y_Vector, NCells_Z_Vector, &
NNodes_X_Vector, NNodes_Y_Vector, NNodes_Z_Vector, &
Nodes_of_Cells_PE, Cells_of_Cells_PE, Flag_Faces_of_Cells, &
allocate_status(2))
! Set mesh type info for an Orthogonal Mesh.
Geometry = 'Cartesian'
Uniformity = 'Nonuniform'
Orthogonality = 'Orthogonal'
! Set Offsets for coordinates on this PE. Note that this is
! the integer offset of the indices, not the physical distance
! offset as is used in the Uniform mesh initialization.
if (this_PE_X == 1) then
Offset_PE_X = 0
else
Offset_PE_X = SUM(NCells_X_Vector(1:this_PE_X-1))
end if
if (NDimensions >= 2) then
if (this_PE_Y == 1) then
Offset_PE_Y = 0
else
Offset_PE_Y = SUM(NCells_Y_Vector(1:this_PE_Y-1))
end if
end if
if (NDimensions == 3) then
if (this_PE_Z == 1) then
Offset_PE_Z = 0
else
Offset_PE_Z = SUM(NCells_Z_Vector(1:this_PE_Z-1))
end if
end if
! Set Node Coordinates on this PE.
call Initialize (Coordinates_Nodes_PE, NDimensions, &
NNodes_Vector(this_PE), allocate_status(3))
node = 1
do pe_z = 1, NNodes_Z_Vector(this_PE_Z)
do pe_y = 1, NNodes_Y_Vector(this_PE_Y)
do pe_x = 1, NNodes_X_Vector(this_PE_X)
Coordinates_Nodes_PE(1,node) = &
Coordinates_Nodes_X(Offset_PE_X + pe_x)
if (NDimensions >= 2) then
Coordinates_Nodes_PE(2,node) = &
Coordinates_Nodes_Y(Offset_PE_Y + pe_y)
end if
if (NDimensions == 3) then
Coordinates_Nodes_PE(3,node) = &
Coordinates_Nodes_Z(Offset_PE_Z + pe_z)
end if
! Increment node.
node = node + 1
end do
end do
end do
VERIFY((node-1)==NNodes_Vector(this_PE),5)
! Initialize the Multi-Mesh object.
call Initialize_Base_Multi_Mesh (Mesh, NDimensions, Geometry, &
Uniformity, Orthogonality, Structure, AMR, Shape, &
NNodes_Vector, NCells_Vector, NFaces_Vector, &
Coordinates_Nodes_PE, Nodes_of_Cells_PE, Mesh_Name, &
allocate_status(4))
! Set physical dimensions of the mesh (Lengths).
call Initialize (Mesh%Lengths, NDimensions, allocate_status(5))
Mesh%Lengths = Lengths
! Set Mesh%Cells_of_Cells_Index and Mesh%Flag_Faces_of_Cells.
call Initialize (Mesh%Cells_of_Cells_Index, Mesh%Cell_Structure, &
Mesh%Cell_Structure, &
Many_of_One_Array=Cells_of_Cells_PE, &
status=allocate_status(6))
call Initialize (Mesh%Flag_Faces_of_Cells, NCells_Vector(this_PE), &
NDimensions*2, allocate_status(7))
Mesh%Flag_Faces_of_Cells = Flag_Faces_of_Cells
! Finalize temporary variables.
call Finalize (Cells_of_Cells_PE, allocate_status(8))
call Finalize (Coordinates_Nodes_PE, allocate_status(9))
call Finalize (Flag_Faces_of_Cells, allocate_status(10))
call Finalize (NCells_Vector, allocate_status(11))
call Finalize (NCells_X_Vector, allocate_status(12))
call Finalize (NCells_Y_Vector, allocate_status(13))
call Finalize (NCells_Z_Vector, allocate_status(14))
call Finalize (NFaces_Vector, allocate_status(15))
call Finalize (NNodes_Vector, allocate_status(16))
call Finalize (NNodes_X_Vector, allocate_status(17))
call Finalize (NNodes_Y_Vector, allocate_status(18))
call Finalize (NNodes_Z_Vector, allocate_status(19))
call Finalize (Nodes_of_Cells_PE, allocate_status(20))
! Consolidate and handle status.
consolidated_status = allocate_status
if (PRESENT(status)) then
WARN_IF(Error(consolidated_status), 5)
status = consolidated_status
else
VERIFY(Normal(consolidated_status), 5)
end if
call Finalize (consolidated_status)
call Finalize (allocate_status)
end subroutine Initialize_Orthogonal_MMesh
subroutine Gen_StructureMesh_Connectivity (NDimensions, Lengths, &
NCells_X_total, NCells_Y_total, NCells_Z_total, Shape, Structure, AMR, &
NPEs_X, NPEs_Y, NPEs_Z, this_PE_X, this_PE_Y, this_PE_Z, &
NNodes_Vector, NCells_Vector, NFaces_Vector, &
NCells_X_Vector, NCells_Y_Vector, NCells_Z_Vector, &
NNodes_X_Vector, NNodes_Y_Vector, NNodes_Z_Vector, &
Nodes_of_Cells_PE, Cells_of_Cells_PE, Flag_Faces_of_Cells, status)
! Use associations.
use Caesar_Mathematics_Module
!~~~~~~~~~~~~~~~~~
! Input variables.
!~~~~~~~~~~~~~~~~~
! Mesh type information.
type(integer), intent(in) :: NDimensions ! Number of Dimensions.
type(real,1,np), intent(in) :: Lengths ! Physical extent of the
! domain in each direction.
! Total number of cells in the X-, Y-, and Z-directions.
type(integer), intent(in) :: NCells_X_total
type(integer), intent(inout), optional :: NCells_Y_total
type(integer), intent(inout), optional :: NCells_Z_total
!~~~~~~~~~~~~~~~~~~
! Output variables.
!~~~~~~~~~~~~~~~~~~
type(Status_type), intent(out), optional :: status ! Exit status.
type(character,name_length), intent(out) :: Shape ! Cell shape.
type(character,name_length), intent(out) :: Structure ! Set to Structured.
type(logical), intent(out) :: AMR ! Set to false.
! Structure length vectors, which give numbers for all PEs.
type(integer,1) :: NCells_Vector ! Number of cells.
type(integer,1) :: NFaces_Vector ! Number of faces.
type(integer,1) :: NNodes_Vector ! Number of nodes.
type(integer,1) :: NCells_X_Vector ! Number of cells in the X direction.
type(integer,1) :: NCells_Y_Vector ! Number of cells in the Y direction.
type(integer,1) :: NCells_Z_Vector ! Number of cells in the Z direction.
type(integer,1) :: NNodes_X_Vector ! Number of nodes in the X direction.
type(integer,1) :: NNodes_Y_Vector ! Number of nodes in the Y direction.
type(integer,1) :: NNodes_Z_Vector ! Number of nodes in the Z direction.
type(integer), intent(out) :: NPEs_X ! Number of PEs in X.
type(integer), intent(out) :: NPEs_Y ! Number of PEs in Y.
type(integer), intent(out) :: NPEs_Z ! Number of PEs in Z.
! Location of this_PE in the PE-mesh.
type(integer), intent(out) :: this_PE_X, this_PE_Y, this_PE_Z
! The nodes for the cells on this PE.
type(integer,2) :: Nodes_of_Cells_PE
! The cells for the cells (that is, across each face) on this PE.
type(integer,2) :: Cells_of_Cells_PE
! Face Flags for Structured Meshes. This is set to: 0 for internal,
! 1 for left (-x), 2 for right (+x), 3 for front (-y), 4 for back (+y),
! 5 for bottom (-z), 6 for top (+z).
type(integer,2) :: Flag_Faces_of_Cells
!~~~~~~~~~~~~~~~~~~~~
! Internal variables.
!~~~~~~~~~~~~~~~~~~~~
type(Status_type), dimension(18) :: allocate_status ! Allocation Status.
type(Status_type) :: consolidated_status ! Consolidated Status.
type(integer), dimension(32) :: NPE_Factors ! Prime factors of NPEs.
type(integer) :: NNPE_Factors ! Number of NPE Factors.
! Structure length vectors, which give numbers for all PEs.
type(integer) :: NCells_total ! Total number of cells on all PEs.
type(integer) :: NFaces_total ! Total number of faces on all PEs.
type(integer) :: NNodes_total ! Total number of nodes on all PEs.
type(integer) :: i, j, k ! Loop parameters.
type(integer) :: cell, cell_x, cell_y, cell_z ! Loop parameters.
type(integer) :: pe, pe_x, pe_y, pe_z ! Loop parameters.
! Nodes_of_Cells evaluation point.
type(integer) :: node_cell_x, node_cell_y, node_cell_z, node_PE
type(integer) :: node_PE_X, node_PE_Y, node_PE_Z
! Cells_of_Cells evaluation point.
type(integer) :: other_cell_x, other_cell_y, other_cell_z, other_cell_PE
type(integer) :: other_cell_PE_X, other_cell_PE_Y, other_cell_PE_Z
! Offsets - starting points for this PE.
type(integer,1) :: Offset_Cells_Vector ! Global number of first cell.
type(integer,1) :: Offset_Faces_Vector ! Global number of first face.
type(integer,1) :: Offset_Nodes_Vector ! Global number of first node.
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! Verify requirements.
! Mesh type info is valid.
VERIFY(NDimensions .InInterval. (/1, 3/),5)
VERIFY(SIZE(Lengths) == NDimensions,5) ! Lengths is correct size.
! Allocations and initializations.
call Initialize (allocate_status)
call Initialize (consolidated_status)
! Initialize vectors.
call Initialize (NCells_Vector, NPEs, allocate_status(1))
call Initialize (NFaces_Vector, NPEs, allocate_status(2))
call Initialize (NNodes_Vector, NPEs, allocate_status(3))
call Initialize (Offset_Cells_Vector, NPEs, allocate_status(4))
call Initialize (Offset_Faces_Vector, NPEs, allocate_status(5))
call Initialize (Offset_Nodes_Vector, NPEs, allocate_status(6))
! Set Shape, NCells_total based on NDimensions.
select case (NDimensions)
case (1)
Shape = "Segmented "
NCells_total = NCells_X_total
NCells_Y_total = 1
NCells_Z_total = 1
NFaces_total = NCells_total + 1
NNodes_total = NCells_total + 1
case (2)
Shape = "Quadrilateral"
NCells_total = NCells_X_total * NCells_Y_total
NCells_Z_total = 1
NFaces_total = 2*NCells_total + NCells_X_total + NCells_Y_total
NNodes_total = (NCells_X_total + 1) * (NCells_Y_total + 1)
case (3)
Shape = "Hexahedral "
NCells_total = NCells_X_total * NCells_Y_total * NCells_Z_total
NFaces_total = 3*NCells_total + NCells_X_total * NCells_Y_total &
+ NCells_Y_total * NCells_Z_total &
+ NCells_Z_total * NCells_X_total
NNodes_total = (NCells_X_total + 1) * &
(NCells_Y_total + 1) * &
(NCells_Z_total + 1)
end select
! Set mesh type info for a Structured Mesh.
AMR = .false.
Structure = 'Structured'
! Determine parallel distributions: the number of PEs in each
! direction. The method used is to divide the problem between PEs so
! that the PE blocks look like a structured mesh with dimensions
! NPEs_X x NPEs_Y x NPEs_Z. Each block will have roughly the same number
! of cells, so load-balancing will be good (in 3D, each block will have
! roughly NCells_PE = (NCells_X_total/NPEs_X) * (NCells_Y_total/NPEs_Y) *
! (NCells_Z_total/NPEs_Z) ).
!
! To balance potential communication costs, the number of PEs in each
! direction are chosen to keep the physical dimensions of the blocks as
! equal as possible. The algorithm chosen is to prime factorize the
! total NPEs, then to iteratively add the largest remaining factor to the
! direction with the largest value of Length/NPEs until all the factors
! are used. This should make the blocks more equal (cubical), but it may
! not be the absolute optimal method. An additional requirement is that
! the NPEs in any direction must remain less than or equal to the number
! of cells in that direction.
!
! Better partitioning could be done if the overlying physics were known,
! however the mesh has no knowledge of the physics, and could indeed be
! used with different physics at different times.
if (NDimensions /= 1) then
call Prime_Factors (NPEs, NNPE_Factors, NPE_Factors)
end if
select case (NDimensions)
case (1)
NPEs_X = NPEs
NPEs_Y = 1
NPEs_Z = 1
VERIFY(NPEs_X <= NCells_X_total,5)
case (2)
NPEs_X = 1
NPEs_Y = 1
NPEs_Z = 1
do while (NNPE_Factors > 0)
! Increase NPEs of direction with largest Length/NPEs, unless
! that would make NPEs > NCells in that direction.
if (Lengths(1)/NPEs_X > Lengths(2)/NPEs_Y .and. &
NPEs_X*NPE_Factors(NNPE_Factors) <= NCells_X_total) then
NPEs_X = NPEs_X * NPE_Factors(NNPE_Factors)
NNPE_Factors = NNPE_Factors - 1
else if (NPEs_Y*NPE_Factors(NNPE_Factors) <= NCells_Y_total) then
NPEs_Y = NPEs_Y * NPE_Factors(NNPE_Factors)
NNPE_Factors = NNPE_Factors - 1
! If unsuccessful, then try to increase NPEs in other directions,
! retaining condition that NPEs must be less than or equal to
! NCells in that direction.
else if (NPEs_X*NPE_Factors(NNPE_Factors) <= NCells_X_total) then
NPEs_X = NPEs_X * NPE_Factors(NNPE_Factors)
NNPE_Factors = NNPE_Factors - 1
! This condition is only reached if this mesh cannot be split up
! on the required number of PEs without leaving some PEs with
! no cells.
else
VERIFY(.false.,1)
end if
end do
case (3)
NPEs_X = 1
NPEs_Y = 1
NPEs_Z = 1
do while (NNPE_Factors > 0)
! Increase NPEs of direction with largest Length/NPEs, unless
! that would make NPEs > NCells in that direction.
if (Lengths(1)/NPEs_X > Lengths(2)/NPEs_Y .and. &
Lengths(1)/NPEs_X > Lengths(3)/NPEs_Z .and. &
NPEs_X*NPE_Factors(NNPE_Factors) <= NCells_X_total) then
NPEs_X = NPEs_X * NPE_Factors(NNPE_Factors)
NNPE_Factors = NNPE_Factors - 1
else if (Lengths(2)/NPEs_Y > Lengths(1)/NPEs_X .and. &
Lengths(2)/NPEs_Y > Lengths(3)/NPEs_Z .and. &
NPEs_Y*NPE_Factors(NNPE_Factors) <= NCells_Y_total) then
NPEs_Y = NPEs_Y * NPE_Factors(NNPE_Factors)
NNPE_Factors = NNPE_Factors - 1
else if (NPEs_Z*NPE_Factors(NNPE_Factors) <= NCells_Z_total) then
NPEs_Z = NPEs_Z * NPE_Factors(NNPE_Factors)
NNPE_Factors = NNPE_Factors - 1
! If unsuccessful, then try to increase NPEs in other directions,
! retaining condition that NPEs must be less than or equal to
! NCells in that direction.
else if (NPEs_X*NPE_Factors(NNPE_Factors) <= NCells_X_total) then
NPEs_X = NPEs_X * NPE_Factors(NNPE_Factors)
NNPE_Factors = NNPE_Factors - 1
else if (NPEs_Y*NPE_Factors(NNPE_Factors) <= NCells_Y_total) then
NPEs_Y = NPEs_Y * NPE_Factors(NNPE_Factors)
NNPE_Factors = NNPE_Factors - 1
! This condition is only reached if this mesh cannot be split up
! on the required number of PEs without leaving some PEs with
! no cells.
else
VERIFY(.false.,1)
end if
end do
end select
! Initialize NCells_D_Vectors, NNodes_D_Vectors.
call Initialize (NCells_X_Vector, NPEs_X, allocate_status(7))
call Initialize (NCells_Y_Vector, NPEs_Y, allocate_status(8))
call Initialize (NCells_Z_Vector, NPEs_Z, allocate_status(9))
call Initialize (NNodes_X_Vector, NPEs_X, allocate_status(10))
call Initialize (NNodes_Y_Vector, NPEs_Y, allocate_status(11))
call Initialize (NNodes_Z_Vector, NPEs_Z, allocate_status(12))
! Set NCells_D_Vectors (divide evenly among the PEs).
call Generate_Even_Distribution (NCells_X_Vector, NCells_X_total)
call Generate_Even_Distribution (NCells_Y_Vector, NCells_Y_total)
call Generate_Even_Distribution (NCells_Z_Vector, NCells_Z_total)
! Set NNodes_D_Vectors -- extra node on last PE in each direction.
NNodes_X_Vector = NCells_X_Vector
NNodes_Y_Vector = NCells_Y_Vector
NNodes_Z_Vector = NCells_Z_Vector
NNodes_X_Vector(NPEs_X) = NNodes_X_Vector(NPEs_X) + 1
if (NDimensions > 1) then
NNodes_Y_Vector(NPEs_Y) = NNodes_Y_Vector(NPEs_Y) + 1
if (NDimensions > 2) then
NNodes_Z_Vector(NPEs_Z) = NNodes_Z_Vector(NPEs_Z) + 1
end if
end if
! Set NCells_Vector, NFaces_Vector, and NNodes_Vector.
pe = 1
do pe_z = 1, NPEs_Z
do pe_y = 1, NPEs_Y
do pe_x = 1, NPEs_X
! Set NCells_Vector.
NCells_Vector(pe) = NCells_X_Vector(pe_x) * &
NCells_Y_Vector(pe_y) * &
NCells_Z_Vector(pe_z)
! Set NFaces_Vector.
NFaces_Vector(pe) = NDimensions * NCells_Vector(pe)
! Add extra YZ faces for pe_x = NPEs_X plane.
if (pe_x == NPEs_X) then
NFaces_Vector(pe) = NFaces_Vector(pe) + &
NCells_Y_Vector(pe_y) * &
NCells_Z_Vector(pe_z)
end if
! Add extra XZ faces for pe_y = NPEs_Y plane in 2-D and 3-D.
if (pe_y == NPEs_Y .and. NDimensions /= 1) then
NFaces_Vector(pe) = NFaces_Vector(pe) + &
NCells_X_Vector(pe_x) * &
NCells_Z_Vector(pe_z)
end if
! Add extra XY faces for pe_z = NPEs_Z plane in 3-D.
if (pe_z == NPEs_Z .and. NDimensions == 3) then
NFaces_Vector(pe) = NFaces_Vector(pe) + &
NCells_X_Vector(pe_x) * &
NCells_Y_Vector(pe_y)
end if
! Set NNodes_Vector.
NNodes_Vector(pe) = NNodes_X_Vector(pe_x) * &
NNodes_Y_Vector(pe_y) * &
NNodes_Z_Vector(pe_z)
! Capture PE info.
if (pe == this_PE) then
this_PE_X = pe_x
this_PE_Y = pe_y
this_PE_Z = pe_z
end if
! Increment pe.
pe = pe + 1
end do
end do
end do
VERIFY((pe-1)==NPEs,8)
VERIFY(SUM(NCells_Vector)==NCells_total,8)
VERIFY(SUM(NFaces_Vector)==NFaces_total,8)
VERIFY(SUM(NNodes_Vector)==NNodes_total,8)
! Set Offset Vectors.
Offset_Cells_Vector(1) = 1
Offset_Nodes_Vector(1) = 1
Offset_Faces_Vector(1) = 1
if (NPEs > 1) then
pe = 1
do pe_z = 1, NPEs_Z
do pe_y = 1, NPEs_Y
do pe_x = 1, NPEs_X
! Skip first pe -- it's already done.
if (pe/=1) then
! Set Offset_Cells_Vector.
Offset_Cells_Vector(pe) = Offset_Cells_Vector(pe-1) + &
NCells_Vector(pe-1)
! Set Offset_Nodes_Vector.
Offset_Nodes_Vector(pe) = Offset_Nodes_Vector(pe-1) + &
NNodes_Vector(pe-1)
! Set Offset_Faces_Vector.
Offset_Faces_Vector(pe) = Offset_Faces_Vector(pe-1) + &
NFaces_Vector(pe-1)
end if
! Increment pe.
pe = pe + 1
end do
end do
end do
end if
! Set Nodes_of_Cells_PE.
call Initialize (Nodes_of_Cells_PE, NCells_Vector(this_PE), &
2**NDimensions, allocate_status(13))
cell = 1
do cell_z = 1, NCells_Z_Vector(this_PE_Z)
do cell_y = 1, NCells_Y_Vector(this_PE_Y)
do cell_x = 1, NCells_X_Vector(this_PE_X)
! These next three loops loop over the nodes in a cell.
do k = 0, 1
do j = 0, 1
do i = 0, 1
! Determine which PE the node is on, and set
! node_cell_D and node_PE to reflect this.
node_cell_x = cell_x + i
node_cell_y = cell_y + j
node_cell_z = cell_z + k
node_PE = this_PE
node_PE_X = this_PE_X
node_PE_Y = this_PE_Y
node_PE_Z = this_PE_Z
if (node_cell_x > NNodes_X_Vector(this_PE_X)) then
node_cell_x = 1
node_PE = node_PE + 1
node_PE_X = node_PE_X + 1
end if
if (node_cell_y > NNodes_Y_Vector(this_PE_Y)) then
node_cell_y = 1
node_PE = node_PE + NPEs_X
node_PE_Y = node_PE_Y + 1
end if
if (node_cell_z > NNodes_Z_Vector(this_PE_Z)) then
node_cell_z = 1
node_PE = node_PE + NPEs_X * NPEs_Y
node_PE_Z = node_PE_Z + 1
end if
! Use formula to set node number for this node within the
! cell, now that the PE that the node is on is known. Note
! that the generalized multi-dimensional node numbering
! (within a cell) is used.
Nodes_of_Cells_PE(cell,i + 2*j + 4*k + 1) = &
Offset_Nodes_Vector(node_PE) + &
(node_cell_x - 1) + &
(node_cell_y - 1) * NNodes_X_Vector(node_PE_X) + &
(node_cell_z - 1) * NNodes_X_Vector(node_PE_X) * &
NNodes_Y_Vector(node_PE_Y)
end do
if (NDimensions < 2) exit
end do
if (NDimensions < 3) exit
end do
! Increment cell.
cell = cell + 1
end do
end do
end do
! Set Cells_of_Cells_PE and Flag_Faces_of_Cells.
call Initialize (Cells_of_Cells_PE, NCells_Vector(this_PE), &
NDimensions*2, allocate_status(14))
call Initialize (Flag_Faces_of_Cells, NCells_Vector(this_PE), &
NDimensions*2, allocate_status(15))
cell = 0
do cell_z = 1, NCells_Z_Vector(this_PE_Z)
do cell_y = 1, NCells_Y_Vector(this_PE_Y)
do cell_x = 1, NCells_X_Vector(this_PE_X)
! Increment cell.
cell = cell + 1
! Loop over X Faces.
do i = -1, 1, 2
! Determine which PE the "other" cell is on, and set
! other_cell_D and other_cell_PE to reflect this. Notice
! the loop-around which is done to allow periodic boundary
! conditions.
other_cell_x = cell_x + i
other_cell_y = cell_y
other_cell_z = cell_z
other_cell_PE = this_PE
other_cell_PE_X = this_PE_X
other_cell_PE_Y = this_PE_Y
other_cell_PE_Z = this_PE_Z
if (other_cell_x > NCells_X_Vector(this_PE_X)) then
other_cell_PE = other_cell_PE + 1
other_cell_PE_X = other_cell_PE_X + 1
if (other_cell_PE_X > NPEs_X) then
other_cell_PE_X = other_cell_PE_X - NPEs_X ! Should = 1.
other_cell_PE = other_cell_PE - NPEs_X
Flag_Faces_of_Cells(cell,2) = 2
end if
other_cell_x = 1
end if
if (other_cell_x < 1) then
other_cell_PE = other_cell_PE - 1
other_cell_PE_X = other_cell_PE_X - 1
if (other_cell_PE_X < 1) then
other_cell_PE_X = other_cell_PE_X + NPEs_X ! Should = NPEs_X.
other_cell_PE = other_cell_PE + NPEs_X
Flag_Faces_of_Cells(cell,1) = 1
end if
other_cell_x = NCells_X_Vector(other_cell_PE_X)
end if
! Use formula to set cell number for the "other" cell, now
! that the PE that the other cell is on is known. Note
! that the generalized multi-dimensional node numbering
! (within a cell) is used. This sets values for faces 1 and 2.
Cells_of_Cells_PE(cell,(i + 3)/2) = &
Offset_Cells_Vector(other_cell_PE) + &
(other_cell_x - 1) + &
(other_cell_y - 1) * NCells_X_Vector(other_cell_PE_X) + &
(other_cell_z - 1) * NCells_X_Vector(other_cell_PE_X) * &
NCells_Y_Vector(other_cell_PE_Y)
end do
if (NDimensions < 2) cycle
! Loop over Y Faces.
do j = -1, 1, 2
! Determine which PE the "other" cell is on, and set
! other_cell_D and other_cell_PE to reflect this. Notice
! the loop-around which is done to allow periodic boundary
! conditions.
other_cell_x = cell_x
other_cell_y = cell_y + j
other_cell_z = cell_z
other_cell_PE = this_PE
other_cell_PE_X = this_PE_X
other_cell_PE_Y = this_PE_Y
other_cell_PE_Z = this_PE_Z
if (other_cell_y > NCells_Y_Vector(this_PE_Y)) then
other_cell_PE = other_cell_PE + NPEs_X
other_cell_PE_Y = other_cell_PE_Y + 1
if (other_cell_PE_Y > NPEs_Y) then
other_cell_PE_Y = other_cell_PE_Y - NPEs_Y ! Should = 1.
other_cell_PE = other_cell_PE - NPEs_Y * NPEs_X
Flag_Faces_of_Cells(cell,4) = 4
end if
other_cell_y = 1
end if
if (other_cell_y < 1) then
other_cell_PE = other_cell_PE - NPEs_X
other_cell_PE_Y = other_cell_PE_Y - 1
if (other_cell_PE_Y < 1) then
other_cell_PE_Y = other_cell_PE_Y + NPEs_Y ! Should = NPEs_Y.
other_cell_PE = other_cell_PE + NPEs_Y * NPEs_X
Flag_Faces_of_Cells(cell,3) = 3
end if
other_cell_y = NCells_Y_Vector(other_cell_PE_Y)
end if
! Use formula to set cell number for the "other" cell, now
! that the PE that the other cell is on is known. Note
! that the generalized multi-dimensional node numbering
! (within a cell) is used. This sets values for faces 3 and 4.
Cells_of_Cells_PE(cell,(j + 7)/2) = &
Offset_Cells_Vector(other_cell_PE) + &
(other_cell_x - 1) + &
(other_cell_y - 1) * NCells_X_Vector(other_cell_PE_X) + &
(other_cell_z - 1) * NCells_X_Vector(other_cell_PE_X) * &
NCells_Y_Vector(other_cell_PE_Y)
end do
if (NDimensions < 3) cycle
! Loop over Z Faces.
do k = -1, 1, 2
! Determine which PE the "other" cell is on, and set
! other_cell_D and other_cell_PE to reflect this. Notice
! the loop-around which is done to allow periodic boundary
! conditions.
other_cell_x = cell_x
other_cell_y = cell_y
other_cell_z = cell_z + k
other_cell_PE = this_PE
other_cell_PE_X = this_PE_X
other_cell_PE_Y = this_PE_Y
other_cell_PE_Z = this_PE_Z
if (other_cell_z > NCells_Z_Vector(this_PE_Z)) then
other_cell_PE = other_cell_PE + NPEs_X * NPEs_Y
other_cell_PE_Z = other_cell_PE_Z + 1
if (other_cell_PE_Z > NPEs_Z) then
other_cell_PE_Z = other_cell_PE_Z - NPEs_Z ! Should = 1.
other_cell_PE = other_cell_PE - NPEs_Z * NPEs_X * NPEs_Y
Flag_Faces_of_Cells(cell,6) = 6
end if
other_cell_z = 1
end if
if (other_cell_z < 1) then
other_cell_PE = other_cell_PE - NPEs_X * NPEs_Y
other_cell_PE_Z = other_cell_PE_Z - 1
if (other_cell_PE_Z < 1) then
other_cell_PE_Z = other_cell_PE_Z + NPEs_Z ! Should = NPEs_Z.
other_cell_PE = other_cell_PE + NPEs_Z * NPEs_X * NPEs_Y
Flag_Faces_of_Cells(cell,5) = 5
end if
other_cell_z = NCells_Z_Vector(other_cell_PE_Z)
end if
! Use formula to set cell number for the "other" cell, now
! that the PE that the other cell is on is known. Note
! that the generalized multi-dimensional node numbering
! (within a cell) is used. This sets values for faces 5 and 6.
Cells_of_Cells_PE(cell,(k + 11)/2) = &
Offset_Cells_Vector(other_cell_PE) + &
(other_cell_x - 1) + &
(other_cell_y - 1) * NCells_X_Vector(other_cell_PE_X) + &
(other_cell_z - 1) * NCells_X_Vector(other_cell_PE_X) * &
NCells_Y_Vector(other_cell_PE_Y)
end do
end do
end do
end do
! Finalize temporary variables.
call Finalize (Offset_Cells_Vector, allocate_status(16))
call Finalize (Offset_Faces_Vector, allocate_status(17))
call Finalize (Offset_Nodes_Vector, allocate_status(18))
! Consolidate and handle status.
consolidated_status = allocate_status
if (PRESENT(status)) then
WARN_IF(Error(consolidated_status), 5)
status = consolidated_status
else
VERIFY(Normal(consolidated_status), 5)
end if
call Finalize (consolidated_status)
call Finalize (allocate_status)
return
end subroutine Gen_StructureMesh_Connectivity