The main documentation of the Get_Area_Faces_of_Cells_Multi_Mesh Procedure contains additional explanation of this code listing.
subroutine Get_Area_Faces_of_Cells_MMesh (Area_Faces_of_Cells, Mesh)
! Input variable.
type(Multi_Mesh_type), intent(inout) :: Mesh ! Mesh object.
! Input/Output variable.
type(real,3) :: Area_Faces_of_Cells ! Area_Faces_of_Cells BNV.
! Internal variables.
type(integer) :: d ! Dimension loop parameter.
! BNV of coordinates of nodes of cells.
type(real,3) :: Coordinates_Nodes_of_Cells
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! Verify requirements.
VERIFY(Valid_State(Mesh),5) ! Mesh is valid.
VERIFY(Valid_State(Area_Faces_of_Cells),5) ! Area_Faces_of_Cells is valid.
! Area_Faces_of_Cells has correct dimensions.
VERIFY(SIZE(Area_Faces_of_Cells,1) == Mesh%NDimensions,5)
VERIFY(SIZE(Area_Faces_of_Cells,2) == Mesh%NCells_PE,5)
VERIFY(SIZE(Area_Faces_of_Cells,3) == Mesh%Faces_per_Cell,5)!Not for polys.
! Get Coordinates of the Nodes if needed.
if (Mesh%Uniformity == "Nonuniform") then
call Initialize (Coordinates_Nodes_of_Cells, Mesh%NDimensions, &
Mesh%NCells_PE, Mesh%Nodes_per_Cell)
call Get_Coordinates_Nodes_of_Cells (Coordinates_Nodes_of_Cells, Mesh)
end if
! Set the vector areas for each face.
Area_Faces_of_Cells = zero
if (Mesh%Uniformity == "Uniform") then
! For a uniform mesh, all cells have the same
! three areas, plus or minus.
do d = 1, Mesh%NDimensions
Area_Faces_of_Cells(d,:,2*d-1) = -Mesh%Area_All_Faces(d)
Area_Faces_of_Cells(d,:,2*d ) = Mesh%Area_All_Faces(d)
end do
else if (Mesh%Orthogonality == "Orthogonal") then
! For an orthogonal mesh, the face areas may be calculated from the
! product of 0, 1 or 2 orthogonal edge lengths. First, set positive
! (+X,+Y,+Z) face areas.
select case (Mesh%NDimensions)
case (1)
! Suppressed Y, Z.
Area_Faces_of_Cells(:,:,2) = one ! Right (+X) Face.
case (2)
! Suppressed Z.
Area_Faces_of_Cells(1,:,2) = & ! Right (+X) Face.
( Coordinates_Nodes_of_Cells(2,:,3) & ! Delta Y -
- Coordinates_Nodes_of_Cells(2,:,1) ) ! Nodes 1 & 3
Area_Faces_of_Cells(2,:,4) = & ! Back (+Y) Face.
( Coordinates_Nodes_of_Cells(1,:,2) & ! Delta X -
- Coordinates_Nodes_of_Cells(1,:,1) ) ! Nodes 1 & 2
case (3)
Area_Faces_of_Cells(1,:,2) = & ! Right (+X) Face.
( Coordinates_Nodes_of_Cells(2,:,3) & ! Delta Y -
- Coordinates_Nodes_of_Cells(2,:,1) ) * & ! Nodes 1 & 3
( Coordinates_Nodes_of_Cells(3,:,5) & ! Delta Z -
- Coordinates_Nodes_of_Cells(3,:,1) ) ! Nodes 1 & 5
Area_Faces_of_Cells(2,:,4) = & ! Back (+Y) Face.
( Coordinates_Nodes_of_Cells(1,:,2) & ! Delta X -
- Coordinates_Nodes_of_Cells(1,:,1) ) * & ! Nodes 1 & 2
( Coordinates_Nodes_of_Cells(3,:,5) & ! Delta Z -
- Coordinates_Nodes_of_Cells(3,:,1) ) ! Nodes 1 & 5
Area_Faces_of_Cells(3,:,6) = & ! Top (+Z) Face.
( Coordinates_Nodes_of_Cells(1,:,2) & ! Delta X -
- Coordinates_Nodes_of_Cells(1,:,1) ) * & ! Nodes 1 & 2
( Coordinates_Nodes_of_Cells(2,:,3) & ! Delta Y -
- Coordinates_Nodes_of_Cells(2,:,1) ) ! Nodes 1 & 3
end select
! Set opposing (-X,-Y,-Z) faces to negative of the positive faces.
do d = 1, Mesh%NDimensions
Area_Faces_of_Cells(d,:,2*d-1) = -Area_Faces_of_Cells(d,:,2*d)
end do
else
! Coding for other mesh types not implemented yet.
VERIFY(.false.,1)
end if
! Finalize temporaries.
if (Mesh%Uniformity == "Nonuniform") then
call Finalize (Coordinates_Nodes_of_Cells)
end if
! Verify guarantees.
VERIFY(Valid_State(Mesh),5) ! Mesh is valid.
VERIFY(Valid_State(Area_Faces_of_Cells),5) ! Area_Faces_of_Cells is valid.
return
end subroutine Get_Area_Faces_of_Cells_MMesh