From f563313ce97c61720d0daafe6f01498278e964fe Mon Sep 17 00:00:00 2001
From: Sharan Roongta <s.roongta@mpie.de>
Date: Sat, 5 Dec 2020 23:52:30 +0100
Subject: [PATCH] PETSc provides subroutine to read physical tags

---
 src/mesh/discretization_mesh.f90 | 47 ++++++++++----------------------
 1 file changed, 14 insertions(+), 33 deletions(-)

diff --git a/src/mesh/discretization_mesh.f90 b/src/mesh/discretization_mesh.f90
index d998c84db..69c86798e 100644
--- a/src/mesh/discretization_mesh.f90
+++ b/src/mesh/discretization_mesh.f90
@@ -83,6 +83,7 @@ subroutine discretization_mesh_init(restart)
   class(tNode), pointer :: &
     num_mesh
   integer :: integrationOrder                                                                       !< order of quadrature rule required
+  
 
   print'(/,a)',   ' <<<+-  discretization_mesh init  -+>>>'
 
@@ -96,12 +97,16 @@ subroutine discretization_mesh_init(restart)
   debug_element = config_debug%get_asInt('element',defaultVal=1)
   debug_ip      = config_debug%get_asInt('integrationpoint',defaultVal=1)
 
+ ! vol_tag = 10
   call DMPlexCreateFromFile(PETSC_COMM_WORLD,interface_geomFile,PETSC_TRUE,globalMesh,ierr)
   CHKERRQ(ierr)
   call DMGetDimension(globalMesh,dimPlex,ierr)
   CHKERRQ(ierr)
   call DMGetStratumSize(globalMesh,'depth',dimPlex,mesh_NcpElemsGlobal,ierr)
   CHKERRQ(ierr)
+  call DMView(globalMesh, PETSC_VIEWER_STDOUT_WORLD,ierr)
+  CHKERRQ(ierr)
+
   ! get number of IDs in face sets (for boundary conditions?)
   call DMGetLabelSize(globalMesh,'Face Sets',mesh_Nboundaries,ierr)
   CHKERRQ(ierr)
@@ -109,6 +114,14 @@ subroutine discretization_mesh_init(restart)
   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)
 
+  if (worldrank == 0) then
+    call DMClone(globalMesh,geomMesh,ierr)
+    CHKERRQ(ierr)
+  else
+    call DMPlexDistribute(globalMesh,0,sf,geomMesh,ierr)
+    CHKERRQ(ierr)
+  endif
+
   allocate(mesh_boundaries(mesh_Nboundaries), source = 0)
   call DMGetLabelSize(globalMesh,'Face Sets',nFaceSets,ierr)
   CHKERRQ(ierr)
@@ -123,38 +136,6 @@ subroutine discretization_mesh_init(restart)
   endif
   call MPI_Bcast(mesh_boundaries,mesh_Nboundaries,MPI_INTEGER,0,PETSC_COMM_WORLD,ierr)
 
-  if (worldrank == 0) then
-    fileContent = IO_readlines(interface_geomFile)
-    l = 0
-    do
-      l = l + 1
-      if (IO_isBlank(fileContent(l))) cycle         ! need also to ignore C and C++ style comments?
-      if (trim(fileContent(l)) == '$Elements') then
-        j = 0
-        l = l + 1
-        do
-          l = l + 1
-          if (trim(fileContent(l)) == '$EndElements') exit
-          chunkPos = IO_stringPos(fileContent(l))
-          if(IO_intValue(fileContent(l),chunkPos,1) == 3) then
-            do k = 1, IO_intValue(fileContent(l),chunkPos,4)
-              call DMSetLabelValue(globalMesh,'material',j,IO_intValue(fileContent(l),chunkPos,2),ierr)
-              CHKERRQ(ierr)
-              j = j + 1
-            enddo
-          endif
-          l = l + IO_intValue(fileContent(l),chunkPos,4)
-        enddo
-        exit
-      endif
-    enddo
-    call DMClone(globalMesh,geomMesh,ierr)
-    CHKERRQ(ierr)
-  else
-    call DMPlexDistribute(globalMesh,0,sf,geomMesh,ierr)
-    CHKERRQ(ierr)
-  endif
-
   call DMDestroy(globalMesh,ierr); CHKERRQ(ierr)
 
   call DMGetStratumSize(geomMesh,'depth',dimPlex,mesh_NcpElems,ierr)
@@ -169,7 +150,7 @@ subroutine discretization_mesh_init(restart)
 
   allocate(materialAt(mesh_NcpElems))
   do j = 1, mesh_NcpElems
-    call DMGetLabelValue(geomMesh,'material',j-1,materialAt(j),ierr)
+    call DMGetLabelValue(geomMesh,'Cell Sets',j-1,materialAt(j),ierr)
     CHKERRQ(ierr)
   end do