diff --git a/src/meshFEM.f90 b/src/meshFEM.f90
index 1e6a95449..80ab80302 100644
--- a/src/meshFEM.f90
+++ b/src/meshFEM.f90
@@ -147,14 +147,19 @@ subroutine mesh_init(ip,el)
  write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
 #include "compilation_info.f90"
 
+ ! read in file
  call DMPlexCreateFromFile(PETSC_COMM_WORLD,geometryFile,PETSC_TRUE,globalMesh,ierr)
  CHKERRQ(ierr)
+ ! get spatial dimension (2 or 3?)
  call DMGetDimension(globalMesh,dimPlex,ierr)
  CHKERRQ(ierr)
+ write(6,*) 'dimension',dimPlex;flush(6)
  call DMGetStratumSize(globalMesh,'depth',dimPlex,mesh_NcpElemsGlobal,ierr)
  CHKERRQ(ierr)
+ ! get number of IDs in face sets (for boundary conditions?)
  call DMGetLabelSize(globalMesh,'Face Sets',mesh_Nboundaries,ierr)
  CHKERRQ(ierr)
+ write(6,*) 'number of "Face Sets"',mesh_Nboundaries;flush(6)
  call MPI_Bcast(mesh_Nboundaries,1,MPI_INTEGER,0,PETSC_COMM_WORLD,ierr)
  call MPI_Bcast(mesh_NcpElemsGlobal,1,MPI_INTEGER,0,PETSC_COMM_WORLD,ierr)
  call MPI_Bcast(dimPlex,1,MPI_INTEGER,0,PETSC_COMM_WORLD,ierr)
@@ -170,7 +175,9 @@ subroutine mesh_init(ip,el)
  enddo  
  if (nFaceSets > 0) call ISRestoreIndicesF90(faceSetIS,pFaceSets,ierr)
  call MPI_Bcast(mesh_boundaries,mesh_Nboundaries,MPI_INTEGER,0,PETSC_COMM_WORLD,ierr)
- 
+
+ ! this read in function should ignore C and C++ style comments
+ ! it is used for BC only?
  if (worldrank == 0) then
    j = 0
    flag = .false.
@@ -179,8 +186,8 @@ subroutine mesh_init(ip,el)
      read(FILEUNIT,'(a512)') line
      if (trim(line) == IO_EOF) exit                                                                    ! skip empty lines
      if (trim(line) == '$Elements') then
-       read(FILEUNIT,'(a512)') line
-       read(FILEUNIT,'(a512)') line
+       read(FILEUNIT,'(a512)') line ! number of elements (ignore)
+       read(FILEUNIT,'(a512)') line 
        flag = .true.  
      endif
      if (trim(line) == '$EndElements') exit