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