diff --git a/README.md b/README.md
index ca5e024..28c6952 100644
--- a/README.md
+++ b/README.md
@@ -33,6 +33,8 @@ It should be noted that if reflective boundary conditions are specified, then th
 If the user tries to assign reflective boundary conditions to a direction with a non-flat boundary, then the converter utility will throw an error and terminate.
 This check is ignored if all boundary conditions are reflective.
 
+If physical volumes are specified by the Gmsh input file, they will be used as the region indices for elements in the THOR mesh. Elements should not have more than one physical volume assigned to them, as only one physical volume will be read in by the mesh converter per element. Assignment of physical volumes *by material* is recommended in order to simplify region mapping in the THOR input file.
+
 ---
 ## Examples
 ---
diff --git a/src/globals.f90 b/src/globals.f90
index d5a82a7..cc64c62 100644
--- a/src/globals.f90
+++ b/src/globals.f90
@@ -37,6 +37,12 @@ MODULE globals
   !element region tags
   INTEGER, ALLOCATABLE :: el_tag(:)
 
+  !physical volumes
+  INTEGER, ALLOCATABLE :: phys_vol_tag(:)
+
+  !element physical volume tag
+  INTEGER, ALLOCATABLE :: el_phys_group(:)
+
   !adjacency list
   INTEGER, ALLOCATABLE :: adj_list(:,:)
 
diff --git a/src/output_thrm.f90 b/src/output_thrm.f90
index 611c304..d29b326 100644
--- a/src/output_thrm.f90
+++ b/src/output_thrm.f90
@@ -28,7 +28,11 @@ SUBROUTINE output_thrm_file()
     ENDDO
     !print out tet regions and source. Assume each region has its own source. User can change this later
     DO i=1,num_tets
-      WRITE(30,'(I0,A,I0,A,I0)')i,' ',el_tag(i),' ',el_tag(i)
+      IF (el_phys_group(i) .EQ. 0) THEN
+        WRITE(30,'(I0,A,I0,A,I0,A,I0)')i,' ',el_tag(i),' ',el_tag(i)
+      ELSE
+        WRITE(30,'(I0,A,I0,A,I0)')i,' ',el_phys_group(i),' ',el_tag(i)
+      END IF
     ENDDO
     !print out tet composition
     DO i=1,num_tets
@@ -101,6 +105,7 @@ SUBROUTINE calcvols()
       IF(tets_in_reg(i) .GT. 0)THEN
         WRITE(*,'(A,I0,A,I0)')'Region ',i,' tets: ',tets_in_reg(i)
         WRITE(*,'(A,I0,A,ES24.16)')'Region ',i,' volume: ',regvol(i)
+        WRITE(*,'(A,I0,A,I0)')'Region ',i,' physical group: ',phys_vol_tag(i)
       ENDIF
     ENDDO
     WRITE(*,'(A,I0)')'Total number of tets: ',SUM(tets_in_reg)
diff --git a/src/read_gmsh.f90 b/src/read_gmsh.f90
index e706b77..6ade2d8 100644
--- a/src/read_gmsh.f90
+++ b/src/read_gmsh.f90
@@ -26,6 +26,18 @@ SUBROUTINE read_gmsh_file()
     temp_str=TRIM(ADJUSTL(temp_str))
     IF(temp_str(1:2) .NE. '4.')STOP 'gmesh format is wrong, should be 4 version'
 
+    !find the entities
+    REWIND(20)
+    DO
+      READ(20,*,iostat=ios)temp_str
+      IF(temp_str .EQ. '$Entities')THEN
+        !and read them
+        CALL read_entities()
+        EXIT
+      ENDIF
+      IF(ios .NE. 0)STOP 'did not find Entities in gmesh file'
+    ENDDO
+
     !find the nodes
     REWIND(20)
     DO
@@ -78,12 +90,12 @@ SUBROUTINE read_nodes()
 
   SUBROUTINE read_elements()
     INTEGER :: num_entities,numel,minel,maxel
-    INTEGER :: temp_int,loc_num_el,i,j,el_indx,el_dim,ent_tag
+    INTEGER :: temp_int,loc_num_el,i,j,el_indx,el_dim,ent_tag,ent_phys
     INTEGER, ALLOCATABLE :: temp_array(:,:)
     !get the elements info
     READ(20,*)num_entities,numel,minel,maxel
     !last value of array is tag (region/material)
-    ALLOCATE(temp_array(numel,5))
+    ALLOCATE(temp_array(numel,6))
     IF(minel .NE. 1 .OR. maxel .NE. numel)STOP 'elements should be indexed 1 to amount'
     !read in all node data
     num_tets=0
@@ -95,11 +107,13 @@ SUBROUTINE read_elements()
       !actually counting the tets
       IF(el_dim .EQ. 3)THEN
         num_tets=num_tets+loc_num_el
+        ent_phys=phys_vol_tag(ent_tag)
         !get data since this is a tet set
         DO j=1,loc_num_el
           el_indx=el_indx+1
           READ(20,*)temp_int,temp_array(el_indx,1:4)
           temp_array(el_indx,5)=ent_tag
+          temp_array(el_indx,6)=ent_phys
           IF(MOD(el_indx,CEILING(numel*1.0/(max_prog-1.0))) .EQ. 0)THEN
             WRITE(*,'(A)',ADVANCE='NO')'*'
             prog=prog+1
@@ -113,12 +127,13 @@ SUBROUTINE read_elements()
       ENDIF
     ENDDO
     !allocate elements and element tags
-    ALLOCATE(element(num_tets,4),el_tag(num_tets))
+    ALLOCATE(element(num_tets,4),el_tag(num_tets),el_phys_group(num_tets))
     element=0
     el_tag=0
     DO i=1,num_tets
       element(i,:)=temp_array(i,1:4)
       el_tag(i)=temp_array(i,5)
+      el_phys_group(i)=temp_array(i,6)
     ENDDO
     DEALLOCATE(temp_array)
     DO i=prog,max_prog
@@ -126,4 +141,29 @@ SUBROUTINE read_elements()
     ENDDO
     WRITE(*,*)
   ENDSUBROUTINE read_elements
+
+  SUBROUTINE read_entities()
+    INTEGER :: num_points, num_curves, num_surfaces, num_volumes
+    INTEGER :: num_skiplines, i, j
+    REAL :: junk(8)
+
+    READ(20,*) num_points, num_curves, num_surfaces, num_volumes
+    num_skiplines = num_points + num_curves + num_surfaces
+    ! skip to volume entities
+    DO i=1,num_skiplines
+      READ(20,*)
+    END DO
+    ! read first physical group of each volume
+    ALLOCATE(phys_vol_tag(num_volumes))
+    phys_vol_tag = 0
+    DO i=1,num_volumes
+      READ(20,*) junk(1:7), j, phys_vol_tag(i)
+      IF (j==0) THEN 
+        phys_vol_tag(i) = 0
+        WRITE(*,'(A,I0,A)')'Region ',i,' has no physical groups associated.'
+      ELSE IF (j .NE. 1) THEN
+        WRITE(*,'(A,I0,A)')'Warning: Region ',i,' has more than one physical group associated!'
+      END IF
+    END DO
+  ENDSUBROUTINE
 END MODULE read_gmsh