Gautam Bisht avatar Gautam Bisht committed d2eb0a9 Merge

merge

Comments (0)

Files changed (19)

Add a comment to this file

regression_tests/mfd/ugrid_45x1x50/ugrid_45x1x50_mesh.h5

Binary file added.

regression_tests/mfd/ugrid_45x1x50/ugrid_45x1x50_mfd.in

+:Description: 3D toy problem for richards equation
+
+:=========================== flow mode ========================================
+MODE RICHARDS
+
+:=========================== chemistry ========================================
+CHEMISTRY
+PRIMARY_SPECIES
+tracer
+/
+/
+
+:=========================== discretization ===================================
+GRID
+TYPE unstructured_mimetic ugrid_45x1x50_mesh.h5
+END
+
+:=========================== fluid properties =================================
+FLUID_PROPERTY
+DIFFUSION_COEFFICIENT 1.d-14
+/
+
+
+:=========================== saturation functions =============================
+SATURATION_FUNCTION sf1
+SATURATION_FUNCTION_TYPE VAN_GENUCHTEN
+RESIDUAL_SATURATION 0.2d0
+LAMBDA 0.25d0
+ALPHA 1.d-4
+MAX_CAPILLARY_PRESSURE 1.d8
+BETAC 0.d0
+POWER 1.d0
+/
+
+:=========================== material properties ==============================
+MATERIAL_PROPERTY soil1
+ID 1
+POROSITY 0.25d0
+TORTUOSITY 0.5d0
+SATURATION_FUNCTION sf1
+PERMEABILITY
+PERM_X 7.75d-12
+PERM_Y 1.d-12
+PERM_Z 3.25d-12
+PERM_XZ 3.8971d-12
+PERM_XY 0
+PERM_YZ 0
+/
+/
+MATERIAL_PROPERTY soil2
+ID 2
+POROSITY 0.25d0
+TORTUOSITY 0.5d0
+SATURATION_FUNCTION sf1
+PERMEABILITY
+PERM_X 7.75d-12
+PERM_Y 1.d-12
+PERM_Z 3.25d-12
+PERM_XZ -3.8971d-12
+PERM_XY 0
+PERM_YZ 0
+/
+/
+
+:=========================== output options ===================================
+OUTPUT
+  PERIODIC TIME  1 y  
+  :FORMAT TECPLOT BLOCK
+  FORMAT HDF5
+  VELOCITIES
+  PERMEABILITY
+/
+
+:=========================== times ============================================
+TIME
+FINAL_TIME 25 y
+:FINAL_TIME 0.0001 y
+INITIAL_TIMESTEP_SIZE 0.01 y
+MAXIMUM_TIMESTEP_SIZE 1 y
+/
+
+REFERENCE_PRESSURE 101325
+
+:=========================== regions ==========================================
+REGION all
+  FILE ugrid_45x1x50_mesh.h5
+  GRID unstructured
+END 
+
+REGION Top_part
+  FILE ugrid_45x1x50_mesh.h5
+  GRID unstructured
+END 
+
+REGION Bottom_part
+  FILE ugrid_45x1x50_mesh.h5
+  GRID unstructured
+END 
+
+REGION Top
+  FILE ugrid_45x1x50_mesh.h5
+  GRID unstructured
+END 
+
+REGION Bottom
+  FILE ugrid_45x1x50_mesh.h5
+  GRID unstructured
+END 
+
+REGION Top_middle
+  FILE ugrid_45x1x50_mesh.h5
+  GRID unstructured
+END 
+
+:=========================== flow conditions ==================================
+FLOW_CONDITION initial
+  TYPE
+    PRESSURE hydrostatic
+   END
+   INTERPOLATION linear
+   DATUM 0.d0 0.d0 0.d0
+   PRESSURE 155800 
+END
+
+FLOW_CONDITION Top
+TYPE
+FLUX neumann
+END
+FLUX 1.58549d-8 !1.58549d-9 = 5 cm/yr darcy
+!FLUX 0.
+END
+
+FLOW_CONDITION Bottom
+TYPE
+PRESSURE hydrostatic
+END
+INTERPOLATION linear
+DATUM 0.d0 0.d0 0.d0
+PRESSURE 155800 
+END
+
+:=========================== transport conditions =============================
+
+TRANSPORT_CONDITION initial_c
+TYPE dirichlet_zero_gradient
+CONSTRAINT_LIST
+0.d0 initial_c
+END
+END
+
+TRANSPORT_CONDITION source_c
+TYPE dirichlet
+CONSTRAINT_LIST
+0.d0 source_c
+END
+END
+
+:=========================== constraints ======================================
+CONSTRAINT initial_c
+CONCENTRATIONS
+tracer    1.d-5   T
+/
+/
+
+CONSTRAINT source_c
+CONCENTRATIONS
+tracer    1.d0   T
+/
+/
+
+
+:=========================== condition couplers ===============================
+: initial condition
+INITIAL_CONDITION
+FLOW_CONDITION initial
+TRANSPORT_CONDITION initial_c
+REGION all
+END
+
+: top boundary condition
+BOUNDARY_CONDITION
+FLOW_CONDITION Top
+TRANSPORT_CONDITION source_c
+REGION Top_middle
+END
+
+: bottom boundary condition
+BOUNDARY_CONDITION
+FLOW_CONDITION Bottom
+TRANSPORT_CONDITION initial_c
+REGION Bottom
+END
+
+
+:=========================== stratigraphy couplers ============================
+STRATA
+REGION Top_part
+MATERIAL soil1
+END
+
+STRATA
+REGION Bottom_part
+MATERIAL soil2
+END

src/pflotran/condition_control.F90

               
         case default
           ! assign initial conditions values to domain
-          if (discretization%itype == STRUCTURED_GRID_MIMETIC) then
+          if (discretization%itype == STRUCTURED_GRID_MIMETIC.or. &
+              discretization%itype == UNSTRUCTURED_GRID_MIMETIC) then
             call GridVecGetArrayF90(grid,field%flow_xx, xx_p, ierr); CHKERRQ(ierr)
             call VecGetArrayF90(field%flow_xx_faces, xx_faces_p, ierr); CHKERRQ(ierr)        
           else
       
             if (.not.associated(initial_condition)) exit
 
-            if (discretization%itype == STRUCTURED_GRID_MIMETIC) then
+            if (discretization%itype == STRUCTURED_GRID_MIMETIC.or. &
+                discretization%itype == UNSTRUCTURED_GRID_MIMETIC) then
 #ifdef DASVYAT
               if (.not.associated(initial_condition%flow_aux_real_var)) then
                 do icell=1,initial_condition%region%num_cells
   call DiscretizationLocalToLocal(discretization,field%iphas_loc,field%iphas_old_loc,ONEDOF)
 
 #ifdef DASVYAT
-  if (discretization%itype == STRUCTURED_GRID_MIMETIC) then
+  if (discretization%itype == STRUCTURED_GRID_MIMETIC.or. &
+      discretization%itype == UNSTRUCTURED_GRID_MIMETIC) then
 
     call VecRestoreArrayF90(field%flow_xx_faces,xx_faces_p, ierr)
     call RealizationSetUpBC4Faces(realization)

src/pflotran/connection.F90

       allocate(connection%face_id(num_connections))
 #ifdef DASVYAT
       allocate(connection%cntr(1:3, num_connections))
+      connection%cntr = 0.d0
 #endif
       connection%id_dn = 0
       connection%dist = 0.d0

src/pflotran/coupler.F90

   coupler => coupler_list%first
   do
     if (.not.associated(coupler)) exit 
-    if (grid%itype == STRUCTURED_GRID_MIMETIC .and. &
+    if ((grid%itype == STRUCTURED_GRID_MIMETIC .or. &
+         grid%discretization_itype == UNSTRUCTURED_GRID_MIMETIC) .and. &
         (coupler%itype == INITIAL_COUPLER_TYPE .or. &
          coupler%itype == BOUNDARY_COUPLER_TYPE)) then  
        call CouplerComputeConnections(grid,option,coupler)
-       call CouplerComputeConnectionsFaces(grid,option,coupler)      
+       call CouplerComputeConnectionsFaces(grid,option,coupler)
        call CouplerAssignBCtoCells(grid,option,coupler)
  
     else
   allocate(coupler%faces_set(coupler%numfaces_set))
   connection_set => ConnectionCreate(ZERO_INTEGER, &
                                    connection_itype)
-  local_faces = 0
-  iconn = 1
-  do icell = 1, region%num_cells
+
+  if (grid%itype==STRUCTURED_GRID_MIMETIC) then
+    !
+    ! For Structured grid
+    !
+
+    local_faces = 0
+    iconn = 1
+    do icell = 1, region%num_cells
     
-    cell_id_local = region%cell_ids(icell)
-    aux_var => mfd_aux%aux_vars(cell_id_local)
+      cell_id_local = region%cell_ids(icell)
+      aux_var => mfd_aux%aux_vars(cell_id_local)
 
-    do iface = 1,aux_var%numfaces
-      face_id_ghosted = aux_var%face_id_gh(iface)
-      face_id_local = grid%fG2L(face_id_ghosted)
+      do iface = 1,aux_var%numfaces
+        face_id_ghosted = aux_var%face_id_gh(iface)
+        face_id_local = grid%fG2L(face_id_ghosted)
 
-      if (coupler%itype == BOUNDARY_COUPLER_TYPE) then
-        conn_set_ptr => grid%faces(face_id_ghosted)%conn_set_ptr
-        conn_id = grid%faces(face_id_ghosted)%id
+        if (coupler%itype == BOUNDARY_COUPLER_TYPE) then
+          conn_set_ptr => grid%faces(face_id_ghosted)%conn_set_ptr
+          conn_id = grid%faces(face_id_ghosted)%id
 
-        if (conn_set_ptr%itype /= BOUNDARY_CONNECTION_TYPE) cycle
-        if (associated(region%faces)) iface_type = region%faces(icell)
-        bnd_face = PETSC_FALSE
-        select case (iface_type)
-          case(WEST_FACE)
-            if ((conn_set_ptr%dist(1,conn_id) == 1).and.&
-                (conn_set_ptr%dist(2,conn_id) == 0).and.&
-                (conn_set_ptr%dist(3,conn_id) == 0)) bnd_face = PETSC_TRUE
-          case(EAST_FACE)
-            if ((conn_set_ptr%dist(1,conn_id) == -1).and.&
-                (conn_set_ptr%dist(2,conn_id) == 0).and.&
-                (conn_set_ptr%dist(3,conn_id) == 0)) bnd_face = PETSC_TRUE
-          case(SOUTH_FACE)
-            if ((conn_set_ptr%dist(1,conn_id) == 0).and.&
-                (conn_set_ptr%dist(2,conn_id) == 1).and.&
-                (conn_set_ptr%dist(3,conn_id) == 0)) bnd_face = PETSC_TRUE
-          case(NORTH_FACE)
-            if ((conn_set_ptr%dist(1,conn_id) == 0).and.&
-                (conn_set_ptr%dist(2,conn_id) == -1).and.&
-                (conn_set_ptr%dist(3,conn_id) == 0)) bnd_face = PETSC_TRUE
-          case(TOP_FACE)
-            if (grid%structured_grid%invert_z_axis) then
+          if (conn_set_ptr%itype /= BOUNDARY_CONNECTION_TYPE) cycle
+          if (associated(region%faces)) iface_type = region%faces(icell)
+          bnd_face = PETSC_FALSE
+          select case (iface_type)
+            case(WEST_FACE)
+              if ((conn_set_ptr%dist(1,conn_id) == 1).and.&
+                  (conn_set_ptr%dist(2,conn_id) == 0).and.&
+                  (conn_set_ptr%dist(3,conn_id) == 0)) bnd_face = PETSC_TRUE
+            case(EAST_FACE)
+              if ((conn_set_ptr%dist(1,conn_id) == -1).and.&
+                  (conn_set_ptr%dist(2,conn_id) == 0).and.&
+                  (conn_set_ptr%dist(3,conn_id) == 0)) bnd_face = PETSC_TRUE
+            case(SOUTH_FACE)
               if ((conn_set_ptr%dist(1,conn_id) == 0).and.&
-                  (conn_set_ptr%dist(2,conn_id) == 0).and.&
-                  (conn_set_ptr%dist(3,conn_id) == 1)) bnd_face = PETSC_TRUE
-            else
+                  (conn_set_ptr%dist(2,conn_id) == 1).and.&
+                  (conn_set_ptr%dist(3,conn_id) == 0)) bnd_face = PETSC_TRUE
+            case(NORTH_FACE)
               if ((conn_set_ptr%dist(1,conn_id) == 0).and.&
-                  (conn_set_ptr%dist(2,conn_id) == 0).and.&
-                  (conn_set_ptr%dist(3,conn_id) == -1)) bnd_face = PETSC_TRUE
+                  (conn_set_ptr%dist(2,conn_id) == -1).and.&
+                  (conn_set_ptr%dist(3,conn_id) == 0)) bnd_face = PETSC_TRUE
+            case(TOP_FACE)
+              if (grid%structured_grid%invert_z_axis) then
+                if ((conn_set_ptr%dist(1,conn_id) == 0).and.&
+                    (conn_set_ptr%dist(2,conn_id) == 0).and.&
+                    (conn_set_ptr%dist(3,conn_id) == 1)) bnd_face = PETSC_TRUE
+              else
+                if ((conn_set_ptr%dist(1,conn_id) == 0).and.&
+                    (conn_set_ptr%dist(2,conn_id) == 0).and.&
+                    (conn_set_ptr%dist(3,conn_id) == -1)) bnd_face = PETSC_TRUE
+              end if
+            case(BOTTOM_FACE)
+              if (grid%structured_grid%invert_z_axis) then
+                if ((conn_set_ptr%dist(1,conn_id) == 0).and.&
+                    (conn_set_ptr%dist(2,conn_id) == 0).and.&
+                    (conn_set_ptr%dist(3,conn_id) == -1)) bnd_face = PETSC_TRUE
+              else
+                if ((conn_set_ptr%dist(1,conn_id) == 0).and.&
+                    (conn_set_ptr%dist(2,conn_id) == 0).and.&
+                    (conn_set_ptr%dist(3,conn_id) == 1)) bnd_face = PETSC_TRUE
+              end if
+          end select
+
+          if (bnd_face)   then
+            if (face_id_local > 0) then
+              coupler%faces_set(iconn) = face_id_ghosted
+              iconn = iconn + 1
             end if
-          case(BOTTOM_FACE)
-            if (grid%structured_grid%invert_z_axis) then
-              if ((conn_set_ptr%dist(1,conn_id) == 0).and.&
-                  (conn_set_ptr%dist(2,conn_id) == 0).and.&
-                  (conn_set_ptr%dist(3,conn_id) == -1)) bnd_face = PETSC_TRUE
-            else
-              if ((conn_set_ptr%dist(1,conn_id) == 0).and.&
-                  (conn_set_ptr%dist(2,conn_id) == 0).and.&
-                  (conn_set_ptr%dist(3,conn_id) == 1)) bnd_face = PETSC_TRUE
+          end if
+
+        else
+          if (face_id_local > 0) then
+            if (local_faces(face_id_local)==0) then
+              coupler%faces_set(iconn) = face_id_ghosted
+              iconn = iconn + 1
+              local_faces(face_id_local) = 1
             end if
-        end select
-
-        if (bnd_face)   then
-          if (face_id_local > 0) then
-            coupler%faces_set(iconn) = face_id_ghosted
-            iconn = iconn + 1
           end if
         end if
+        end do
+      end do
 
-      else
-        if (face_id_local > 0) then
-          if (local_faces(face_id_local)==0) then
-            coupler%faces_set(iconn) = face_id_ghosted
-            iconn = iconn + 1
-            local_faces(face_id_local) = 1
+    else
+
+    !
+    ! For Untructured grid
+    !
+
+    local_faces = 0
+    iconn = 1
+    do icell = 1, region%num_cells
+    
+      cell_id_local = region%cell_ids(icell)
+      aux_var => mfd_aux%aux_vars(cell_id_local)
+
+      do iface = 1,aux_var%numfaces
+        face_id_ghosted = aux_var%face_id_gh(iface)
+        face_id_local = grid%fG2L(face_id_ghosted)
+
+        if (coupler%itype == BOUNDARY_COUPLER_TYPE) then
+
+          coupler%faces_set(iconn)=grid%fU2M(region%faces(icell),region%cell_ids(icell))
+          iconn=iconn+1
+          exit
+
+        else
+          if (face_id_local > 0) then
+            if (local_faces(face_id_local)==0) then
+              coupler%faces_set(iconn) = face_id_ghosted
+              iconn = iconn + 1
+              local_faces(face_id_local) = 1
+            end if
           end if
         end if
-      end if
       end do
     end do
 
-    select case(coupler%itype)
-      case(INITIAL_COUPLER_TYPE)
-        connection_set%itype = INITIAL_CONNECTION_TYPE
-      case(SRC_SINK_COUPLER_TYPE)
-        connection_set%itype = SRC_SINK_CONNECTION_TYPE
-      case(BOUNDARY_COUPLER_TYPE)
-        connection_set%itype = BOUNDARY_CONNECTION_TYPE
-    end select
+  endif
+
+  select case(coupler%itype)
+    case(INITIAL_COUPLER_TYPE)
+      connection_set%itype = INITIAL_CONNECTION_TYPE
+    case(SRC_SINK_COUPLER_TYPE)
+      connection_set%itype = SRC_SINK_CONNECTION_TYPE
+    case(BOUNDARY_COUPLER_TYPE)
+      connection_set%itype = BOUNDARY_CONNECTION_TYPE
+  end select
   deallocate(local_faces)
 
 #endif

src/pflotran/discretization.F90

             call InputReadNChars(input,option,discretization%filename,MAXSTRINGLENGTH, &
                                  PETSC_TRUE)
             call InputErrorMsg(input,option,'unstructured filename','GRID')
+          case('unstructured_mimetic')
+            discretization%itype = UNSTRUCTURED_GRID_MIMETIC
+            option%mimetic = PETSC_TRUE
+            word = discretization%ctype
+            discretization%ctype = 'unstructured_mimetic'
+            unstructured_grid_itype = IMPLICIT_UNSTRUCTURED_GRID
+            unstructured_grid_ctype = 'implicit unstructured'
+            call InputReadNChars(input,option,discretization%filename,MAXSTRINGLENGTH, &
+                                 PETSC_TRUE)
+            call InputErrorMsg(input,option,'unstructured filename','GRID')
           case('structured_mimetic')
             discretization%itype = STRUCTURED_GRID_MIMETIC
             option%mimetic = PETSC_TRUE
   
   grid => GridCreate()
   select case(discretization%itype)
-    case(UNSTRUCTURED_GRID)
+    case(UNSTRUCTURED_GRID,UNSTRUCTURED_GRID_MIMETIC)
       un_str_grid => UGridCreate()
       select case(unstructured_grid_itype)
         case(IMPLICIT_UNSTRUCTURED_GRID)
       grid%itype = discretization%itype
       grid%ctype = discretization%ctype
   end select
+  grid%discretization_itype=discretization%itype
   discretization%grid => grid
 
 end subroutine DiscretizationReadRequiredCards
       discretization%dm_index_to_ndof(NPHASEDOF) = option%nphase
       discretization%dm_index_to_ndof(NFLOWDOF) = option%nflowdof
       discretization%dm_index_to_ndof(NTRANDOF) = option%ntrandof
-    case(UNSTRUCTURED_GRID)
+    case(UNSTRUCTURED_GRID,UNSTRUCTURED_GRID_MIMETIC)
 
       ! petsc will call parmetis to calculate the graph/dual
 #if !defined(PETSC_HAVE_PARMETIS)
       discretization%grid%ngmax = discretization%grid%structured_grid%ngmax
       discretization%grid%global_offset = &
         discretization%grid%structured_grid%global_offset
-    case(UNSTRUCTURED_GRID)
+    case(UNSTRUCTURED_GRID,UNSTRUCTURED_GRID_MIMETIC)
       discretization%grid%nmax = discretization%grid%unstructured_grid%nmax
       discretization%grid%nlmax = discretization%grid%unstructured_grid%nlmax
       discretization%grid%ngmax = discretization%grid%unstructured_grid%ngmax
     case(STRUCTURED_GRID, STRUCTURED_GRID_MIMETIC)
       call StructGridCreateDM(discretization%grid%structured_grid, &
                                   dm_ptr%dm,ndof,stencil_width,stencil_type,option)
-    case(UNSTRUCTURED_GRID)
+    case(UNSTRUCTURED_GRID,UNSTRUCTURED_GRID_MIMETIC)
       call UGridCreateUGDM(discretization%grid%unstructured_grid, &
                            dm_ptr%ugdm,ndof,option)
       call DMShellCreate(option%mycomm,dm_ptr%dm,ierr)
         case(NATURAL)
           call DMDACreateNaturalVector(dm_ptr%dm,vector,ierr)
       end select
-    case(UNSTRUCTURED_GRID)
+    case(UNSTRUCTURED_GRID,UNSTRUCTURED_GRID_MIMETIC)
       call UGridDMCreateVector(discretization%grid%unstructured_grid, &
                                dm_ptr%ugdm,vector, &
                                vector_type,option)
   
 end function DiscretizationGetDMPtrFromIndex
 
+! ************************************************************************** !
+! ************************************************************************** !
 function DiscretizationGetDMCPtrFromIndex(discretization,dm_index)
 
   implicit none
                                  dm_ptr%ugdm,mat_type,Jacobian,option)
       call MatSetOption(Jacobian,MAT_KEEP_NONZERO_PATTERN,PETSC_FALSE,ierr)
       call MatSetOption(Jacobian,MAT_ROW_ORIENTED,PETSC_FALSE,ierr)
+    case(UNSTRUCTURED_GRID_MIMETIC)
+      select case(dm_index)
+        case(NFLOWDOF)
+          call MFDCreateJacobianLP(discretization%grid, discretization%MFD, mat_type, Jacobian, option)
+          call MatSetOption(Jacobian,MAT_KEEP_NONZERO_PATTERN,PETSC_FALSE,ierr)
+          call MatSetOption(Jacobian,MAT_ROW_ORIENTED,PETSC_FALSE,ierr)
+        case(NTRANDOF)
+          call UGridDMCreateJacobian(discretization%grid%unstructured_grid, &
+                                     dm_ptr%ugdm,mat_type,Jacobian,option)
+          call MatSetOption(Jacobian,MAT_KEEP_NONZERO_PATTERN,PETSC_FALSE,ierr)
+          call MatSetOption(Jacobian,MAT_ROW_ORIENTED,PETSC_FALSE,ierr)
+      end select
     case(STRUCTURED_GRID_MIMETIC)
 #ifdef DASVYAT
       select case(dm_index)
   
     
   select case(discretization%itype)
-    case(STRUCTURED_GRID_MIMETIC)
+    case(STRUCTURED_GRID_MIMETIC,UNSTRUCTURED_GRID_MIMETIC)
       call DiscretizGlobalToLocalFacesBegin(discretization,global_vec,local_vec,dm_index)
       call DiscretizGlobalToLocalFacesEnd(discretization,global_vec,local_vec,dm_index)
   end select
   
     
   select case(discretization%itype)
-    case(STRUCTURED_GRID_MIMETIC)
+    case(STRUCTURED_GRID_MIMETIC,UNSTRUCTURED_GRID_MIMETIC)
       call DiscretizGlobalToLocalLPBegin(discretization,global_vec,local_vec,dm_index)
       call DiscretizGlobalToLocalLPEnd(discretization,global_vec,local_vec,dm_index)
   end select
     case(STRUCTURED_GRID,STRUCTURED_GRID_MIMETIC)
       call DMLocalToGlobalBegin(dm_ptr%dm,local_vec,INSERT_VALUES,global_vec,ierr)
       call DMLocalToGlobalEnd(dm_ptr%dm,local_vec,INSERT_VALUES,global_vec,ierr)
-   case(UNSTRUCTURED_GRID)
+   case(UNSTRUCTURED_GRID,UNSTRUCTURED_GRID_MIMETIC)
       call VecScatterBegin(dm_ptr%ugdm%scatter_ltog,local_vec,global_vec, &
                            INSERT_VALUES,SCATTER_FORWARD,ierr)
       call VecScatterEnd(dm_ptr%ugdm%scatter_ltog,local_vec,global_vec, &
     case(STRUCTURED_GRID,STRUCTURED_GRID_MIMETIC)
       call DMDALocalToLocalBegin(dm_ptr%dm,local_vec1,INSERT_VALUES,local_vec2,ierr)
       call DMDALocalToLocalEnd(dm_ptr%dm,local_vec1,INSERT_VALUES,local_vec2,ierr)
-    case(UNSTRUCTURED_GRID)
+    case(UNSTRUCTURED_GRID,UNSTRUCTURED_GRID_MIMETIC)
       call VecScatterBegin(dm_ptr%ugdm%scatter_ltol,local_vec1,local_vec2, &
                            INSERT_VALUES,SCATTER_FORWARD,ierr)
       call VecScatterEnd(dm_ptr%ugdm%scatter_ltol,local_vec1,local_vec2, &
   
   
   select case(discretization%itype)
-    case(STRUCTURED_GRID,STRUCTURED_GRID_MIMETIC)
+    case(STRUCTURED_GRID,STRUCTURED_GRID_MIMETIC,UNSTRUCTURED_GRID_MIMETIC)
       call DiscretizLocalToLocalFacesBegin(discretization,local_vec1,local_vec2,dm_index)
       call DiscretizLocalToLocalFacesEnd(discretization,local_vec1,local_vec2,dm_index)
     case(UNSTRUCTURED_GRID)
   
   
   select case(discretization%itype)
-    case(STRUCTURED_GRID,STRUCTURED_GRID_MIMETIC)
+    case(STRUCTURED_GRID,STRUCTURED_GRID_MIMETIC,UNSTRUCTURED_GRID_MIMETIC)
       call DiscretizLocalToLocalLPBegin(discretization,local_vec1,local_vec2,dm_index)
       call DiscretizLocalToLocalLPEnd(discretization,local_vec1,local_vec2,dm_index)
     case(UNSTRUCTURED_GRID)
     case(STRUCTURED_GRID,STRUCTURED_GRID_MIMETIC)
       call DMDAGlobalToNaturalBegin(dm_ptr%dm,global_vec,INSERT_VALUES,natural_vec,ierr)
       call DMDAGlobalToNaturalEnd(dm_ptr%dm,global_vec,INSERT_VALUES,natural_vec,ierr)
-    case(UNSTRUCTURED_GRID)
+    case(UNSTRUCTURED_GRID,UNSTRUCTURED_GRID_MIMETIC)
       call VecScatterBegin(dm_ptr%ugdm%scatter_gton,global_vec,natural_vec, &
                            INSERT_VALUES,SCATTER_FORWARD,ierr)
       call VecScatterEnd(dm_ptr%ugdm%scatter_gton,global_vec,natural_vec, &
     case(STRUCTURED_GRID,STRUCTURED_GRID_MIMETIC)
       call DMDANaturalToGlobalBegin(dm_ptr%dm,natural_vec,INSERT_VALUES,global_vec,ierr)
       call DMDANaturalToGlobalEnd(dm_ptr%dm,natural_vec,INSERT_VALUES,global_vec,ierr)
-    case(UNSTRUCTURED_GRID)
+    case(UNSTRUCTURED_GRID,UNSTRUCTURED_GRID_MIMETIC)
       call VecScatterBegin(dm_ptr%ugdm%scatter_ntog,natural_vec,global_vec, &
                            INSERT_VALUES,SCATTER_FORWARD,ierr)
       call VecScatterEnd(dm_ptr%ugdm%scatter_ntog,natural_vec,global_vec, &
   select case(discretization%itype)
     case(STRUCTURED_GRID,STRUCTURED_GRID_MIMETIC)
       call DMGlobalToLocalBegin(dm_ptr%dm,global_vec,INSERT_VALUES,local_vec,ierr)
-    case(UNSTRUCTURED_GRID)
+    case(UNSTRUCTURED_GRID,UNSTRUCTURED_GRID_MIMETIC)
       call VecScatterBegin(dm_ptr%ugdm%scatter_gtol,global_vec,local_vec, &
                            INSERT_VALUES,SCATTER_FORWARD,ierr)
   end select
   select case(discretization%itype)
     case(STRUCTURED_GRID,STRUCTURED_GRID_MIMETIC)
       call DMGlobalToLocalEnd(dm_ptr%dm,global_vec,INSERT_VALUES,local_vec,ierr)
-    case(UNSTRUCTURED_GRID)
+    case(UNSTRUCTURED_GRID,UNSTRUCTURED_GRID_MIMETIC)
       call VecScatterEnd(dm_ptr%ugdm%scatter_gtol,global_vec,local_vec, &
                          INSERT_VALUES,SCATTER_FORWARD,ierr)
   end select
   PetscErrorCode :: ierr
   
   select case(discretization%itype)
-    case(STRUCTURED_GRID_MIMETIC)
+    case(STRUCTURED_GRID_MIMETIC,UNSTRUCTURED_GRID_MIMETIC)
       call VecScatterBegin( discretization%MFD%scatter_gtol_LP, global_vec, local_vec, &
                                 INSERT_VALUES,SCATTER_FORWARD, ierr)
   end select
   
   
   select case(discretization%itype)
-    case(STRUCTURED_GRID_MIMETIC)
+    case(STRUCTURED_GRID_MIMETIC,UNSTRUCTURED_GRID_MIMETIC)
       call VecScatterEnd( discretization%MFD%scatter_gtol_LP, global_vec, local_vec, &
                                 INSERT_VALUES,SCATTER_FORWARD, ierr)
   end select
   select case(discretization%itype)
     case(STRUCTURED_GRID,STRUCTURED_GRID_MIMETIC)
       call DMDALocalToLocalBegin(dm_ptr%dm,local_vec1,INSERT_VALUES,local_vec2,ierr)
-    case(UNSTRUCTURED_GRID)
+    case(UNSTRUCTURED_GRID,UNSTRUCTURED_GRID_MIMETIC)
       call VecScatterBegin(dm_ptr%ugdm%scatter_ltol,local_vec1,local_vec2, &
                            INSERT_VALUES,SCATTER_FORWARD,ierr)
   end select
   select case(discretization%itype)
     case(STRUCTURED_GRID,STRUCTURED_GRID_MIMETIC)
       call DMDALocalToLocalEnd(dm_ptr%dm,local_vec1,INSERT_VALUES,local_vec2,ierr)
-    case(UNSTRUCTURED_GRID)
+    case(UNSTRUCTURED_GRID,UNSTRUCTURED_GRID_MIMETIC)
       call VecScatterEnd(dm_ptr%ugdm%scatter_ltol,local_vec1,local_vec2, &
                          INSERT_VALUES,SCATTER_FORWARD,ierr)    
   end select
   select case(discretization%itype)
     case(STRUCTURED_GRID,STRUCTURED_GRID_MIMETIC)
       call DMDAGlobalToNaturalBegin(dm_ptr%dm,global_vec,INSERT_VALUES,natural_vec,ierr)
-    case(UNSTRUCTURED_GRID)
+    case(UNSTRUCTURED_GRID,UNSTRUCTURED_GRID_MIMETIC)
       call VecScatterBegin(dm_ptr%ugdm%scatter_gton,global_vec,natural_vec, &
                            INSERT_VALUES,SCATTER_FORWARD,ierr)
   end select
   select case(discretization%itype)
     case(STRUCTURED_GRID,STRUCTURED_GRID_MIMETIC)
       call DMDAGlobalToNaturalEnd(dm_ptr%dm,global_vec,INSERT_VALUES,natural_vec,ierr)
-    case(UNSTRUCTURED_GRID)
+    case(UNSTRUCTURED_GRID,UNSTRUCTURED_GRID_MIMETIC)
       call VecScatterEnd(dm_ptr%ugdm%scatter_gton,global_vec,natural_vec, &
                          INSERT_VALUES,SCATTER_FORWARD,ierr)       
   end select
   select case(discretization%itype)
     case(STRUCTURED_GRID,STRUCTURED_GRID_MIMETIC)
       call DMDAGetAO(discretization%dm_1dof,ao,ierr)
-    case(UNSTRUCTURED_GRID)
+    case(UNSTRUCTURED_GRID,UNSTRUCTURED_GRID_MIMETIC)
       ao = discretization%grid%unstructured_grid%ao_natural_to_petsc
   end select
   call AOApplicationToPetsc(ao,size(int_array),int_array,ierr)
         deallocate(discretization%dmc_ntrandof)
         nullify(discretization%dmc_ntrandof)
       endif
-    case(UNSTRUCTURED_GRID)
+    case(UNSTRUCTURED_GRID,UNSTRUCTURED_GRID_MIMETIC)
       if (associated(discretization%dm_1dof%ugdm)) &
         call UGridDMDestroy(discretization%dm_1dof%ugdm)
       if (associated(discretization%dm_nflowdof%ugdm)) &

src/pflotran/grid.F90

   PetscInt, parameter, public :: STRUCTURED_GRID = 1
   PetscInt, parameter, public :: UNSTRUCTURED_GRID = 2
   PetscInt, parameter, public :: STRUCTURED_GRID_MIMETIC = 3
+  PetscInt, parameter, public :: UNSTRUCTURED_GRID_MIMETIC = 6
 
   type, public :: grid_type 
   
     Vec :: e2f             ! global vector to establish connection between global face_id and cell_id
     Vec :: e2n, e2n_LP     ! global cell connectivity vector
 
+    ! For mapping faces between unstrutured grid and mimetic discretization
+    PetscInt, pointer :: fU2M(:,:),fM2U(:)
+    PetscInt :: discretization_itype
+
     PetscReal, pointer :: x(:), y(:), z(:) ! coordinates of ghosted grid cells
 
     PetscReal :: x_min_global, y_min_global, z_min_global
   grid%e2n_LP = 0
 #endif
 
+  nullify(grid%fU2M)
+  nullify(grid%fM2U)
   nullify(grid%hash)
   grid%num_hash_bins = 1000
+  grid%discretization_itype = 0
 
   GridCreate => grid
 
   endif
 #endif
 
+#ifdef MFD_UGRID
+  if(grid%itype==IMPLICIT_UNSTRUCTURED_GRID) then
+
+    call GridPopulateFaces(grid, option)
+
+    ! Note: Allocation of memory has to happen after call to GridPopulateFaces()
+    allocate(grid%fL2G(grid%nlmax_faces))
+    allocate(grid%fG2L(grid%ngmax_faces))
+    grid%fL2G = 0
+    grid%fG2L = 0
+  endif
+#endif
+
 end subroutine GridComputeInternalConnect
 
 ! ************************************************************************** !
    type(face_type), pointer :: faces(:)
    character(len=MAXWORDLENGTH) :: filename
 
+#ifdef MFD_UGRID
+  if (grid%itype==IMPLICIT_UNSTRUCTURED_GRID) then
+    call GridPopulateFacesForUGrid(grid,option)
+    return
+  endif
+#endif
+
    total_faces = grid%ngmax_faces
    allocate(faces(total_faces))
    face_id = 0
    cur_connection_set => connection_set_list%first
    do 
      if (.not.associated(cur_connection_set)) exit
-     do iconn = 1, cur_connection_set%num_connections 
+     do iconn = 1, cur_connection_set%num_connections
        face_id = face_id + 1
        faces(face_id)%conn_set_ptr => cur_connection_set
        faces(face_id)%id = iconn
    enddo
 
    grid%faces => faces
+
 #endif
 end subroutine GridPopulateFaces 
 
 
   PetscInt, pointer :: numfaces(:)
 
+#ifdef MFD_UGRID
+  if (grid%itype==IMPLICIT_UNSTRUCTURED_GRID) then
+    call GridComputeCell2FaceForUGrid(grid,MFD_aux,option)
+    return
+  endif
+#endif
+
   MFD_aux => MFDAuxCreate()
   grid%MFD => MFD_aux
  
     if (conn%itype==BOUNDARY_CONNECTION_TYPE) then
       ghosted_id_dn = conn%id_dn(iface)
       local_id_dn = grid%nG2L(ghosted_id_dn) ! Ghost to local mapping
-   
+
       if (local_id_dn>0) then
         aux_var => MFD_aux%aux_vars(local_id_dn)
         call MFDAuxAddFace(aux_var,option, icount)
 
   VecScatter :: VC_global2ghosted
 
+#ifdef MFD_UGRID
+  if (grid%itype==IMPLICIT_UNSTRUCTURED_GRID) then
+    call GridSetGlobalCell2FaceForUGrid(grid,MFD_aux,DOF,option)
+    return
+  endif
+#endif
+
   select case(DOF)
     case(ONEDOF)
       ndof = 1
       ndof = option%nflowdof
     case(NTRANDOF)
       ndof = option%ntrandof
-   end select
-
-   MFD_aux%ndof = ndof
-
-    allocate(grid%fL2P(grid%nlmax_faces))
-    allocate(grid%fG2P(grid%ngmax_faces))
-
-    global_offset = 0
-    do iface = 1,grid%nlmax_faces
-       grid%fL2P(iface)=grid%global_faces_offset + grid%global_cell_offset + iface - 1
-    end do
-
-    global_offset = grid%global_faces_offset + grid%global_cell_offset
-
-    stride = 6                ! Only for hexagons
-
-    call VecCreate(option%mycomm, grid%e2f, ierr)
-    call VecSetSizes(grid%e2f, grid%nlmax*stride, PETSC_DECIDE, ierr)
-    call VecSetFromOptions(grid%e2f, ierr)
-
-    call VecDuplicate(grid%e2f, grid%e2n, ierr)
-
-    call VecGetArrayF90(grid%e2f, e2f_local_values, ierr)
-    call VecGetArrayF90(grid%e2n, e2n_local_values, ierr)
-
-
-    e2f_local_values = 0
-    e2n_local_values = 0
-
-    do icell = 1, grid%nlmax
-      aux_var => MFD_aux%aux_vars(icell) 
-      do icount = 1, aux_var%numfaces
-        ghost_face_id = aux_var%face_id_gh(icount)
-        local_face_id = grid%fG2L(ghost_face_id)
+  end select
+
+  MFD_aux%ndof = ndof
+
+  allocate(grid%fL2P(grid%nlmax_faces))
+  allocate(grid%fG2P(grid%ngmax_faces))
+
+  global_offset = 0
+  do iface = 1,grid%nlmax_faces
+    grid%fL2P(iface)=grid%global_faces_offset + grid%global_cell_offset + iface - 1
+  enddo
+
+  global_offset = grid%global_faces_offset + grid%global_cell_offset
+
+  stride = 6                ! Only for hexagons
+
+  call VecCreate(option%mycomm, grid%e2f, ierr)
+  call VecSetSizes(grid%e2f, grid%nlmax*stride, PETSC_DECIDE, ierr)
+  call VecSetFromOptions(grid%e2f, ierr)
+
+  call VecDuplicate(grid%e2f, grid%e2n, ierr)
+
+  call VecGetArrayF90(grid%e2f, e2f_local_values, ierr)
+  call VecGetArrayF90(grid%e2n, e2n_local_values, ierr)
+
+  e2f_local_values = 0
+  e2n_local_values = 0
+
+  do icell = 1, grid%nlmax
+    aux_var => MFD_aux%aux_vars(icell)
+    do icount = 1, aux_var%numfaces
+      ghost_face_id = aux_var%face_id_gh(icount)
+      local_face_id = grid%fG2L(ghost_face_id)
+      conn => grid%faces(ghost_face_id)%conn_set_ptr
+      iface = grid%faces(ghost_face_id)%id
+
+      if (conn%itype==INTERNAL_CONNECTION_TYPE) then
+
+        if (local_face_id > 0) then
+
+          e2f_local_values((icell-1)*stride + icount) = grid%fL2P(local_face_id) + 1
+          ghosted_id_up = conn%id_up(iface)
+          ghosted_id_dn = conn%id_dn(iface)
+
+          local_id_up = grid%nG2L(ghosted_id_up) ! = zero for ghost nodes
+          local_id_dn = grid%nG2L(ghosted_id_dn) ! Ghost to local mapping
+          if (local_id_up==icell) then
+            e2n_local_values((icell-1)*stride + icount) = grid%nG2P(ghosted_id_dn) + 1
+          else if (local_id_dn==icell) then
+            e2n_local_values((icell-1)*stride + icount) = grid%nG2P(ghosted_id_up) + 1
+          endif
+
+        else if (local_face_id == 0) then
+
+          ghosted_id_up = conn%id_up(iface)
+          ghosted_id_dn = conn%id_dn(iface)
+
+          local_id_up = grid%nG2L(ghosted_id_up) ! = zero for ghost nodes
+          local_id_dn = grid%nG2L(ghosted_id_dn) ! Ghost to local mapping
+          if (local_id_up==icell) then
+            e2n_local_values((icell-1)*stride + icount) = grid%nG2P(ghosted_id_dn) + 1
+          else if (local_id_dn==icell) then
+            e2n_local_values((icell-1)*stride + icount) = grid%nG2P(ghosted_id_up) + 1
+          endif
+        endif
+
+      else if (conn%itype==BOUNDARY_CONNECTION_TYPE) then
+        e2n_local_values((icell-1)*stride + icount) = 0
+        if (local_face_id > 0) e2f_local_values((icell-1)*stride + icount) = grid%fL2P(local_face_id) + 1
+      endif
+    enddo
+  enddo
+
+  call VecRestoreArrayF90(grid%e2f, e2f_local_values, ierr)
+  call VecRestoreArrayF90(grid%e2n, e2n_local_values, ierr)
+
+  allocate(ghosted_ids(grid%ngmax_faces - grid%nlmax_faces))
+  ghosted_ids = 0
+
+  num_ghosted_upd = 0
+  do icount = 1, grid%ngmax_faces
+    conn => grid%faces(icount)%conn_set_ptr
+    iface = grid%faces(icount)%id
+    if (conn%itype==INTERNAL_CONNECTION_TYPE) then
+      if (conn%local(iface) == 0) then
+        num_ghosted_upd = num_ghosted_upd + 1
+        ghosted_id_up = conn%id_up(iface)
+        ghosted_id_dn = conn%id_dn(iface)
+        if (grid%nG2L(ghosted_id_up)==0) then    ! up_cell is ghosted
+          ghosted_ids(num_ghosted_upd) = grid%nG2P(ghosted_id_up) + 1          ! +1 since global indexes strats from 0
+        else if (grid%nG2L(ghosted_id_dn)==0) then    ! down_cell is ghosted
+          ghosted_ids(num_ghosted_upd) = grid%nG2P(ghosted_id_dn) + 1
+        endif
+      endif
+    endif
+  enddo
+
+  call VecCreateSeq(PETSC_COMM_SELF, num_ghosted_upd*stride, ghosted_e2f, ierr)
+  call VecCreateSeq(PETSC_COMM_SELF, num_ghosted_upd*stride, ghosted_e2n, ierr)
+         
+  allocate(strided_indices_local(num_ghosted_upd))
+  allocate(strided_indices_ghosted(num_ghosted_upd))
+
+  do icount = 1, num_ghosted_upd
+    strided_indices_local(icount) = (icount -1)
+    strided_indices_ghosted(icount) = (ghosted_ids(icount)-1)
+  enddo
+
+  call ISCreateBlock(option%mycomm, stride, num_ghosted_upd, strided_indices_local, &
+                    PETSC_COPY_VALUES, is_local_bl, ierr)
+  call ISCreateBlock(option%mycomm, stride, num_ghosted_upd, strided_indices_ghosted,&
+                    PETSC_COPY_VALUES, is_a2g_bl, ierr)
+
+  call VecScatterCreate(grid%e2f, is_a2g_bl, ghosted_e2f, is_local_bl, VC_global2ghosted, ierr)
+
+  call ISDestroy(is_local_bl, ierr)
+  call ISDestroy(is_a2g_bl, ierr)
+
+  call VecScatterBegin(VC_global2ghosted, grid%e2f, ghosted_e2f, &
+                      INSERT_VALUES,SCATTER_FORWARD,ierr)
+  call VecScatterEnd(VC_global2ghosted, grid%e2f, ghosted_e2f, &
+                      INSERT_VALUES,SCATTER_FORWARD,ierr)
+
+  call VecScatterBegin(VC_global2ghosted, grid%e2n, ghosted_e2n, &
+                      INSERT_VALUES,SCATTER_FORWARD,ierr)
+  call VecScatterEnd(VC_global2ghosted, grid%e2n, ghosted_e2n, &
+                    INSERT_VALUES,SCATTER_FORWARD,ierr)
+
+  call VecGetArrayF90(ghosted_e2n, vec_ptr_e2n_gh, ierr)
+  call VecGetArrayF90(ghosted_e2f, vec_ptr_e2f_gh, ierr)
+
+  call VecGetArrayF90(grid%e2n, vec_ptr_e2n, ierr)
+  call VecGetArrayF90(grid%e2f, vec_ptr_e2f, ierr)
+ 
+  do icell = 1, grid%nlmax
+    aux_var => MFD_aux%aux_vars(icell)
+    do icount = 1, aux_var%numfaces
+      ghost_face_id = aux_var%face_id_gh(icount)
+      local_face_id = grid%fG2L(ghost_face_id)
+
+      if (local_face_id == 0) then
         conn => grid%faces(ghost_face_id)%conn_set_ptr
         iface = grid%faces(ghost_face_id)%id
-        if (conn%itype==INTERNAL_CONNECTION_TYPE) then
-
-          if (local_face_id > 0) then
-
-            e2f_local_values((icell-1)*stride + icount) = grid%fL2P(local_face_id) + 1
-            ghosted_id_up = conn%id_up(iface)
-            ghosted_id_dn = conn%id_dn(iface)
-
-            local_id_up = grid%nG2L(ghosted_id_up) ! = zero for ghost nodes
-            local_id_dn = grid%nG2L(ghosted_id_dn) ! Ghost to local mapping
-            if (local_id_up==icell) then
-              e2n_local_values((icell-1)*stride + icount) = grid%nG2P(ghosted_id_dn) + 1
-            else if (local_id_dn==icell) then
-              e2n_local_values((icell-1)*stride + icount) = grid%nG2P(ghosted_id_up) + 1
-            endif
-
-          else if (local_face_id == 0) then
-
-            ghosted_id_up = conn%id_up(iface)
-            ghosted_id_dn = conn%id_dn(iface)
-
-            local_id_up = grid%nG2L(ghosted_id_up) ! = zero for ghost nodes
-            local_id_dn = grid%nG2L(ghosted_id_dn) ! Ghost to local mapping
-            if (local_id_up==icell) then
-              e2n_local_values((icell-1)*stride + icount) = grid%nG2P(ghosted_id_dn) + 1
-            else if (local_id_dn==icell) then
-              e2n_local_values((icell-1)*stride + icount) = grid%nG2P(ghosted_id_up) + 1
-            endif
+
+        ghosted_id_up = conn%id_up(iface)
+        ghosted_id_dn = conn%id_dn(iface)
+
+        local_id_up = grid%nG2L(ghosted_id_up) ! = zero for ghost nodes
+        local_id_dn = grid%nG2L(ghosted_id_dn) ! Ghost to local mapping
+        if (icell==local_id_up) then
+          global_neigh_id = grid%nG2P(ghosted_id_dn) + 1
+        else if (icell==local_id_dn) then
+          global_neigh_id = grid%nG2P(ghosted_id_up) + 1
+        endif
+        do jcount = 1, num_ghosted_upd
+          if (ghosted_ids(jcount)==global_neigh_id) exit
+        enddo
+        do iface=1, stride             ! assumption that cell has 6 faces
+          if (vec_ptr_e2n_gh((jcount-1)*stride + iface) == grid%nG2P(grid%nL2G(icell)) + 1) then
+            vec_ptr_e2f((icell -1)*stride +icount) = vec_ptr_e2f_gh((jcount-1)*stride + iface)
+            grid%fG2P(ghost_face_id) = int(vec_ptr_e2f_gh((jcount-1)*stride + iface)) - 1
+            exit
           endif
-
-        else if (conn%itype==BOUNDARY_CONNECTION_TYPE) then
-          e2n_local_values((icell-1)*stride + icount) = 0
-          if (local_face_id > 0) e2f_local_values((icell-1)*stride + icount) = grid%fL2P(local_face_id) + 1
-        endif
-      enddo
-   enddo
-
-   call VecRestoreArrayF90(grid%e2f, e2f_local_values, ierr)
-   call VecRestoreArrayF90(grid%e2n, e2n_local_values, ierr)
-
-    allocate(ghosted_ids(grid%ngmax_faces - grid%nlmax_faces))
-    ghosted_ids = 0
-
-    num_ghosted_upd = 0
-    do icount = 1, grid%ngmax_faces
-      conn => grid%faces(icount)%conn_set_ptr
-      iface = grid%faces(icount)%id
-      if (conn%itype==INTERNAL_CONNECTION_TYPE) then
-        if (conn%local(iface) == 0) then
-          num_ghosted_upd = num_ghosted_upd + 1
-          ghosted_id_up = conn%id_up(iface)
-          ghosted_id_dn = conn%id_dn(iface)
-          if (grid%nG2L(ghosted_id_up)==0) then    ! up_cell is ghosted
-            ghosted_ids(num_ghosted_upd) = grid%nG2P(ghosted_id_up) + 1          ! +1 since global indexes strats from 0
-          else if (grid%nG2L(ghosted_id_dn)==0) then    ! down_cell is ghosted
-            ghosted_ids(num_ghosted_upd) = grid%nG2P(ghosted_id_dn) + 1
-          endif
-        endif
+        end do
+      else
+        grid%fG2P(ghost_face_id) = grid%fL2P(local_face_id)
       endif
     enddo
-
-    call VecCreateSeq(PETSC_COMM_SELF, num_ghosted_upd*stride, ghosted_e2f, ierr)
-    call VecCreateSeq(PETSC_COMM_SELF, num_ghosted_upd*stride, ghosted_e2n, ierr)
-         
-    allocate(strided_indices_local(num_ghosted_upd))
-    allocate(strided_indices_ghosted(num_ghosted_upd))
-
-    do icount = 1, num_ghosted_upd
-      strided_indices_local(icount) = (icount -1)
-      strided_indices_ghosted(icount) = (ghosted_ids(icount)-1)
-    enddo
-
-    call ISCreateBlock(option%mycomm, stride, num_ghosted_upd, strided_indices_local, &
-                      PETSC_COPY_VALUES, is_local_bl, ierr)
-    call ISCreateBlock(option%mycomm, stride, num_ghosted_upd, strided_indices_ghosted,&
-                      PETSC_COPY_VALUES, is_a2g_bl, ierr)
-
-    call VecScatterCreate(grid%e2f, is_a2g_bl, ghosted_e2f, is_local_bl, VC_global2ghosted, ierr)
-
-    call ISDestroy(is_local_bl, ierr)
-    call ISDestroy(is_a2g_bl, ierr)
-
-    call VecScatterBegin(VC_global2ghosted, grid%e2f, ghosted_e2f, &
-                        INSERT_VALUES,SCATTER_FORWARD,ierr)
-    call VecScatterEnd(VC_global2ghosted, grid%e2f, ghosted_e2f, &
-                        INSERT_VALUES,SCATTER_FORWARD,ierr) 
-
-    call VecScatterBegin(VC_global2ghosted, grid%e2n, ghosted_e2n, &
-                        INSERT_VALUES,SCATTER_FORWARD,ierr)
-    call VecScatterEnd(VC_global2ghosted, grid%e2n, ghosted_e2n, &
-                      INSERT_VALUES,SCATTER_FORWARD,ierr)
-
-    call VecGetArrayF90(ghosted_e2n, vec_ptr_e2n_gh, ierr)
-    call VecGetArrayF90(ghosted_e2f, vec_ptr_e2f_gh, ierr)
-
-    call VecGetArrayF90(grid%e2n, vec_ptr_e2n, ierr)
-    call VecGetArrayF90(grid%e2f, vec_ptr_e2f, ierr)
- 
-    do icell = 1, grid%nlmax
-      aux_var => MFD_aux%aux_vars(icell) 
-      do icount = 1, aux_var%numfaces
-        ghost_face_id = aux_var%face_id_gh(icount)
-        local_face_id = grid%fG2L(ghost_face_id)
- 
-        if (local_face_id == 0) then
-          conn => grid%faces(ghost_face_id)%conn_set_ptr
-          iface = grid%faces(ghost_face_id)%id
-
-          ghosted_id_up = conn%id_up(iface)
-          ghosted_id_dn = conn%id_dn(iface)
-
-          local_id_up = grid%nG2L(ghosted_id_up) ! = zero for ghost nodes
-          local_id_dn = grid%nG2L(ghosted_id_dn) ! Ghost to local mapping
-          if (icell==local_id_up) then
-            global_neigh_id = grid%nG2P(ghosted_id_dn) + 1
-          else if (icell==local_id_dn) then
-            global_neigh_id = grid%nG2P(ghosted_id_up) + 1
-          endif
-          do jcount = 1, num_ghosted_upd
-            if (ghosted_ids(jcount)==global_neigh_id) exit
-          enddo
-          do iface=1, stride             ! assumption that cell has 6 faces
-            if (vec_ptr_e2n_gh((jcount-1)*stride + iface) == grid%nG2P(grid%nL2G(icell)) + 1) then
-              vec_ptr_e2f((icell -1)*stride +icount) = vec_ptr_e2f_gh((jcount-1)*stride + iface)
-              grid%fG2P(ghost_face_id) = vec_ptr_e2f_gh((jcount-1)*stride + iface) - 1
-              exit
-            endif
-          end do
-        else
-          grid%fG2P(ghost_face_id) = grid%fL2P(local_face_id)
-        endif
-      enddo
-    enddo
-
-    call VecRestoreArrayF90(grid%e2n, vec_ptr_e2n, ierr)
-    call VecRestoreArrayF90(grid%e2f, vec_ptr_e2f, ierr)
-
-    call VecRestoreArrayF90(ghosted_e2n, vec_ptr_e2n_gh, ierr)
-    call VecRestoreArrayF90(ghosted_e2f, vec_ptr_e2f_gh, ierr)
-
-    call VecDestroy(ghosted_e2n, ierr)
-    call VecDestroy(ghosted_e2f, ierr)
-    call VecScatterDestroy(VC_global2ghosted, ierr)
-
-    call CreateMFDStruct4LP(grid, MFD_aux, ndof, option)
-
-    deallocate(ghosted_ids)
-    deallocate(strided_indices_local)
-    deallocate(strided_indices_ghosted)
+  enddo
+
+  call VecRestoreArrayF90(grid%e2n, vec_ptr_e2n, ierr)
+  call VecRestoreArrayF90(grid%e2f, vec_ptr_e2f, ierr)
+  call VecRestoreArrayF90(ghosted_e2n, vec_ptr_e2n_gh, ierr)
+  call VecRestoreArrayF90(ghosted_e2f, vec_ptr_e2f_gh, ierr)
+
+  call VecDestroy(ghosted_e2n, ierr)
+  call VecDestroy(ghosted_e2f, ierr)
+  call VecScatterDestroy(VC_global2ghosted, ierr)
+
+  call CreateMFDStruct4LP(grid, MFD_aux, ndof, option)
+
+  deallocate(ghosted_ids)
+  deallocate(strided_indices_local)
+  deallocate(strided_indices_ghosted)
 #endif
 
 end subroutine GridComputeGlobalCell2FaceConnectivity
        int_array, PETSC_COPY_VALUES, MFD_aux%is_ghosted_local_LP, ierr)
   deallocate(int_array)
 
+  call VecScatterCreate(local_vec_LP,MFD_aux%is_local_local_LP,global_vec_LP, &
+                        MFD_aux%is_local_petsc_LP,MFD_aux%scatter_ltog_LP,ierr)
+
+
   allocate(int_array(NG))
   do id = 1, grid%ngmax_faces
     int_array(id) = (grid%fG2P(id))
   call ISLocalToGlobalMappingBlock(MFD_aux%mapping_ltog_LP, ndof, &
                                    MFD_aux%mapping_ltogb_LP, ierr)
 
-  call VecScatterCreate(local_vec_LP,MFD_aux%is_local_local_LP,global_vec_LP, &
-                        MFD_aux%is_local_petsc_LP,MFD_aux%scatter_ltog_LP,ierr)
-
   call VecScatterCreate(global_vec_LP,MFD_aux%is_ghosted_petsc_LP,local_vec_LP, &
                         MFD_aux%is_ghosted_local_LP, MFD_aux%scatter_gtol_LP, ierr)
 
 
 end subroutine GridSaveBoundaryCellInfo
 
+! ************************************************************************** !
+!> This routine populates faces of a unstructured grid for MIMETIC
+!! discretization. It does the following:
+!! 1) Calculates nlmax_faces, ngmax_faces
+!! 2) If boundary faces are present, adds boundary connection in 
+!!    grid%boundary_connection_set_list
+!! 3) Lastly, save information about ghosted faces: faces(ngmax_faces)
+!!
+!! Note: This subroutine performs functions of StructGridComputeInternConnect(),
+!!       StructGridComputeBoundConnect(), and GridPopulateFaces() for a
+!!       STRUCTURED_GRID_MIMETIC.
+!!
+!> @author
+!! Gautam Bisht, LBNL
+!!
+!! date: 03/29/13
+! ************************************************************************** !
+subroutine GridPopulateFacesForUGrid(grid,option)
+
+  use Option_module
+  use Unstructured_Cell_module
+  use Connection_module
+
+  implicit none
+
+  type(grid_type) :: grid
+  type(option_type) :: option
+
+#if defined(MFD_UGRID)
+  type(face_type), pointer :: faces(:)
+  type(unstructured_grid_type),pointer :: ugrid
+  type(connection_set_list_type), pointer :: connection_set_list
+  type(connection_set_type), pointer :: cur_connection_set
+  type(connection_set_type), pointer :: bnd_connections
+  type(connection_set_type), pointer :: conn
+
+  PetscInt :: iconn
+  PetscInt :: nconn
+  PetscInt :: icell
+  PetscInt :: iface
+  PetscInt :: local_id
+  PetscInt :: ghosted_id
+  PetscInt :: ghosted_id_up
+  PetscInt :: ghosted_id_dn
+  PetscInt :: face_id
+  PetscInt :: nfaces
+  PetscInt :: nfaces_intrn_loc
+  PetscInt :: nfaces_intrn_nonloc
+  PetscInt :: nfaces_bnd
+  PetscInt :: offset
+  PetscInt,pointer::int_array(:)
+  PetscReal, pointer :: vec_ptr(:)
+
+  Vec :: proc_id_loc
+  Vec :: proc_id_nat
+  Vec :: proc_id_ghosted
+  VecScatter :: vec_scat
+
+  IS :: is_tmp1, is_tmp2
+
+  PetscErrorCode :: ierr
+  PetscViewer :: viewer
+
+  ugrid => grid%unstructured_grid
+
+  nfaces_bnd = 0
+  nfaces_intrn_loc = 0
+  nfaces_intrn_nonloc = 0
+
+  ! Compute number of boundary faces
+  do local_id = 1, ugrid%nlmax
+    nfaces = UCellGetNFaces(ugrid%cell_type(local_id),option)
+    do iface = 1, nfaces
+      face_id = ugrid%cell_to_face_ghosted(iface,local_id)
+      if (ugrid%face_to_cell_ghosted(2,face_id) < 1) then
+        ! boundary face, since not connected to 2 cells
+        nfaces_bnd = nfaces_bnd + 1
+      endif
+    enddo
+  enddo
+
+  ! For local+ghost cells, determine processor id on which the cell resides
+  call VecCreateMPI(option%mycomm,ugrid%nlmax,PETSC_DETERMINE,proc_id_loc,ierr)
+  call VecCreateMPI(option%mycomm,ugrid%nlmax,PETSC_DETERMINE,proc_id_nat,ierr)
+  call VecCreateMPI(option%mycomm,ugrid%ngmax,PETSC_DETERMINE,proc_id_ghosted,ierr)
+
+  call VecGetArrayF90(proc_id_loc,vec_ptr,ierr)
+  do local_id = 1,ugrid%nlmax
+    vec_ptr(local_id)=option%myrank
+  enddo
+  call VecRestoreArrayF90(proc_id_loc,vec_ptr,ierr)
+
+  allocate(int_array(ugrid%nlmax))
+  do local_id=1,ugrid%nlmax
+    int_array(local_id)=local_id-1+ugrid%global_offset
+  enddo
+  call ISCreateGeneral(option%mycomm,ugrid%nlmax, &
+                       int_array,PETSC_COPY_VALUES,is_tmp1,ierr)
+
+  do local_id=1,ugrid%nlmax
+    int_array(local_id)=grid%nG2A(grid%nL2G(local_id))-1
+  enddo
+  call ISCreateGeneral(option%mycomm,ugrid%nlmax, &
+                       int_array,PETSC_COPY_VALUES,is_tmp2,ierr)
+  deallocate(int_array)
+
+  call VecScatterCreate(proc_id_loc,is_tmp1,proc_id_nat,is_tmp2,vec_scat,ierr)
+  call ISDestroy(is_tmp1,ierr)
+  call ISDestroy(is_tmp2,ierr)
+
+  call VecScatterBegin(vec_scat,proc_id_loc,proc_id_nat, &
+                       INSERT_VALUES,SCATTER_FORWARD,ierr)
+  call VecScatterEnd(vec_scat,proc_id_loc,proc_id_nat, &
+                     INSERT_VALUES,SCATTER_FORWARD,ierr)
+  call VecScatterDestroy(vec_scat,ierr)
+
+  offset=0
+  call MPI_Exscan(ugrid%ngmax,offset, &
+                  ONE_INTEGER_MPI,MPIU_INTEGER,MPI_SUM,option%mycomm,ierr)
+  allocate(int_array(ugrid%ngmax))
+  do ghosted_id=1,ugrid%ngmax
+    int_array(ghosted_id)=ghosted_id-1+offset
+  enddo
+  call ISCreateGeneral(option%mycomm,ugrid%ngmax, &
+                       int_array,PETSC_COPY_VALUES,is_tmp2,ierr)
+
+  do ghosted_id=1,ugrid%ngmax
+    int_array(ghosted_id)=grid%nG2A(ghosted_id)-1
+  enddo
+  call ISCreateGeneral(option%mycomm,ugrid%ngmax, &
+                       int_array,PETSC_COPY_VALUES,is_tmp1,ierr)
+  deallocate(int_array)
+
+  call VecScatterCreate(proc_id_nat,is_tmp1,proc_id_ghosted,is_tmp2,vec_scat,ierr)
+  call ISDestroy(is_tmp1,ierr)
+  call ISDestroy(is_tmp2,ierr)
+
+  call VecScatterBegin(vec_scat,proc_id_nat,proc_id_ghosted, &
+                       INSERT_VALUES,SCATTER_FORWARD,ierr)
+  call VecScatterEnd(vec_scat,proc_id_nat,proc_id_ghosted, &
+                     INSERT_VALUES,SCATTER_FORWARD,ierr)
+  call VecScatterDestroy(vec_scat,ierr)
+
+  ! Compute number of the two types of internal faces:
+  ! - nfaces_intrn_loc: Both cells sharing the face are local cells
+  ! - nfaces_intrn_nonloc: One of the cells sharing the face is a ghost cell
+  connection_set_list => grid%internal_connection_set_list
+  cur_connection_set => connection_set_list%first
+  call VecGetArrayF90(proc_id_ghosted,vec_ptr,ierr)
+  do 
+    if (.not.associated(cur_connection_set)) exit
+    do iconn = 1, cur_connection_set%num_connections
+      ghosted_id_up = cur_connection_set%id_up(iconn)
+      ghosted_id_dn = cur_connection_set%id_dn(iconn)
+      if (ghosted_id_up<=ugrid%nlmax .and. ghosted_id_dn<=ugrid%nlmax) then
+        ! Both cells sharing the face are local cells
+        nfaces_intrn_loc = nfaces_intrn_loc + 1
+        cur_connection_set%local(iconn) = 1
+      else
+        if((ghosted_id_dn>ugrid%nlmax.and.vec_ptr(ghosted_id_dn)>option%myrank).or. &
+           (ghosted_id_up>ugrid%nlmax.and.vec_ptr(ghosted_id_up)>option%myrank)) then
+          ! If the ghost cell resides on a proc whose rank is greater than
+          ! mine, the face is local
+          nfaces_intrn_loc = nfaces_intrn_loc + 1
+          cur_connection_set%local(iconn) = 1
+        else
+          ! Otherwise, the face is non-local
+          nfaces_intrn_nonloc = nfaces_intrn_nonloc + 1
+        endif
+      endif
+    enddo
+    cur_connection_set => cur_connection_set%next
+  enddo
+  call VecRestoreArrayF90(proc_id_ghosted,vec_ptr,ierr)
+
+  ! Save information about number of faces
+  grid%nlmax_faces = nfaces_intrn_loc+nfaces_bnd
+  grid%ngmax_faces = nfaces_intrn_loc+nfaces_bnd+nfaces_intrn_nonloc
+
+  ! If there are any boundary faces, set up Boundary connections for the grid
+  ! Note: Equivalent subroutine for STRUCTURED_GRID_MIMETIC is 
+  ! StructGridComputeBoundConnect()
+  if (nfaces_bnd>0) then
+    ! allocate memory
+    nullify(bnd_connections)
+    bnd_connections => ConnectionCreate(nfaces_bnd,BOUNDARY_CONNECTION_TYPE)
+
+    ! Populate 'bnd_connections'
+    nconn = 0
+    do local_id = 1, ugrid%nlmax
+      nfaces = UCellGetNFaces(ugrid%cell_type(local_id),option)
+      do iface = 1, nfaces
+        face_id = ugrid%cell_to_face_ghosted(iface,local_id)
+        if (ugrid%face_to_cell_ghosted(2,face_id) < 1) then
+          ! boundary face, since not connected to 2 cells
+          nconn=nconn+1
+          bnd_connections%id_dn(nconn) = local_id
+          call GridPopulateConnection(grid,bnd_connections,iface,nconn, &
+                                    local_id,option)
+          bnd_connections%cntr(:,nconn)=bnd_connections%intercp(:,nconn)
+        endif
+      enddo
+    enddo
+
+    ! Add to list
+    allocate(grid%boundary_connection_set_list)
+    call ConnectionInitList(grid%boundary_connection_set_list)
+    call ConnectionAddToList(bnd_connections,grid%boundary_connection_set_list)
+  endif
+
+  ! Allocate memory to save 'faces'
+  allocate(faces(grid%ngmax_faces))
+
+  ! Add all internal-faces
+  face_id = 0
+  connection_set_list => grid%internal_connection_set_list
+  cur_connection_set => connection_set_list%first
+  do
+    if (.not.associated(cur_connection_set)) exit
+    do iconn = 1, cur_connection_set%num_connections
+      face_id = face_id + 1
+      faces(face_id)%conn_set_ptr => cur_connection_set
+      faces(face_id)%id = iconn
+    enddo
+    cur_connection_set => cur_connection_set%next
+  enddo
+
+  ! Add any boundary-faces
+  connection_set_list => grid%boundary_connection_set_list
+  cur_connection_set => connection_set_list%first
+  do
+    if (.not.associated(cur_connection_set)) exit
+    do iconn = 1, cur_connection_set%num_connections
+      face_id = face_id + 1
+      faces(face_id)%conn_set_ptr => cur_connection_set
+      faces(face_id)%id = iconn
+    enddo
+    cur_connection_set => cur_connection_set%next
+  enddo
+
+  ! Save it
+  grid%faces => faces
+
+  ! Free up memory
+  call VecDestroy(proc_id_loc,ierr)
+  call VecDestroy(proc_id_nat,ierr)
+  call VecDestroy(proc_id_ghosted,ierr)
+#endif
+
+end subroutine GridPopulateFacesForUGrid
+
+! ************************************************************************** !
+!> This routine sets up cell to face connection for an unstructure grid for
+!! MIMETIC discretization.
+!!
+!> @author
+!! Gautam Bisht, LBNL
+!!
+!! date: 03/29/13
+! ************************************************************************** !
+subroutine GridComputeCell2FaceForUGrid(grid,MFD,option)
+
+  use MFD_Aux_module
+  use Option_module
+  use Unstructured_Cell_module
+
+  implicit none
+
+  type(grid_type) :: grid
+  type(mfd_type), pointer :: MFD
+  type(option_type) :: option
+
+#ifdef DASVYAT
+
+  type(mfd_auxvar_type), pointer :: aux_var
+  type(connection_set_type), pointer :: conn
+  type(unstructured_grid_type),pointer :: ugrid
+  type(connection_set_list_type), pointer :: connection_set_list
+  type(connection_set_type), pointer :: cur_connection_set
+
+  PetscInt :: iconn
+  PetscInt :: icell
+  PetscInt :: iface
+  PetscInt :: local_id
+  PetscInt :: local_id_up
+  PetscInt :: local_id_dn
+  PetscInt :: ghosted_id
+  PetscInt :: ghosted_id_up
+  PetscInt :: ghosted_id_dn
+  PetscInt :: face_id
+  PetscInt :: nfaces
+  PetscInt :: nfaces_intrn_loc
+  PetscInt :: nfaces_intrn_nonloc
+  PetscInt :: nfaces_bnd
+  PetscInt :: offset
+  PetscInt :: iface2
+  PetscInt :: face_id2
+  PetscInt :: count
+
+  PetscInt,pointer :: bnd_count(:)
+  
+  PetscBool :: found
+  PetscErrorCode :: ierr
+  
+  character(len=MAXWORDLENGTH) :: filename
+
+  MFD => MFDAuxCreate()
+  grid%MFD => MFD
+  ugrid => grid%unstructured_grid
+ 
+  ! Test
+  allocate(grid%fU2M(MAX_FACE_PER_CELL,grid%nlmax))
+  allocate(grid%fM2U(grid%ngmax_faces))
+  allocate(bnd_count(grid%nlmax))
+  bnd_count=0
+  grid%fM2U = -1
+  grid%fU2M = -1
+
+  ! Find offset for boundary faces
+  offset=0
+  nfaces_bnd=0
+  
+  connection_set_list => grid%internal_connection_set_list
+  cur_connection_set => connection_set_list%first
+  do
+    if (.not.associated(cur_connection_set)) exit
+    offset=offset+cur_connection_set%num_connections
+    cur_connection_set => cur_connection_set%next
+  enddo
+
+  call MFDAuxInit(MFD,grid%nlmax,option)
+
+  ! For each local cell, allocate memory for mfd_auxvar_type that depends on
+  ! on number of faces for a given cell.
+  do local_id = 1, grid%nlmax
+    aux_var => MFD%aux_vars(local_id)
+    nfaces = UCellGetNFaces(ugrid%cell_type(local_id),option)
+    call MFDAuxVarInit(aux_var,nfaces,option)
+  enddo
+
+  ! Compute fG2L and fL2G mapping
+  local_id = 1
+  do iface = 1, grid%ngmax_faces
+    conn => grid%faces(iface)%conn_set_ptr
+    face_id = grid%faces(iface)%id
+
+    if (conn%itype==INTERNAL_CONNECTION_TYPE) then
+
+      ghosted_id_up = conn%id_up(face_id)
+      ghosted_id_dn = conn%id_dn(face_id)
+
+      local_id_up = grid%nG2L(ghosted_id_up) ! = zero for ghost nodes
+      local_id_dn = grid%nG2L(ghosted_id_dn) ! Ghost to local mapping
+
+      if (conn%local(face_id) == 1) then
+        grid%fG2L(iface)=local_id
+        grid%fL2G(local_id) = iface
+        local_id = local_id + 1
+      endif
+
+      if (local_id_dn>0) then
+        aux_var => MFD%aux_vars(local_id_dn)
+        call MFDAuxAddFace(aux_var,option,iface)
+      endif
+
+      if (local_id_up>0) then
+        aux_var => MFD%aux_vars(local_id_up)
+        call MFDAuxAddFace(aux_var,option,iface)
+      endif
+
+      ! For 'local_id_up', find the face-id that is shared by cells
+      ! local_id_up ----- local_id_dn
+      found=PETSC_FALSE
+      do iface2=1,MAX_FACE_PER_CELL
+        face_id2=ugrid%cell_to_face_ghosted(iface2,local_id_up)
+        if( (ugrid%face_to_cell_ghosted(1,face_id2)==ghosted_id_up.and. &
+             ugrid%face_to_cell_ghosted(2,face_id2)==ghosted_id_dn).or. &
+            (ugrid%face_to_cell_ghosted(1,face_id2)==ghosted_id_dn.and. &
+             ugrid%face_to_cell_ghosted(2,face_id2)==ghosted_id_up)) then
+          found=PETSC_TRUE
+          grid%fU2M(iface2,local_id_up)=iface
+          grid%fM2U(iface)=face_id2
+          exit
+        endif
+      enddo
+      if(.not.found) then
+        option%io_buffer='1) UGRID face not found to match face in MIMETIC discretization'
+        call printErrMsg(option)
+      endif
+      
+      ! For 'local_id_dn', find the face-id that is shared by cells
+      ! local_id_up ----- local_id_dn
+      if(ghosted_id_dn<=grid%nlmax) then
+        found=PETSC_FALSE
+        do iface2=1,MAX_FACE_PER_CELL
+          face_id2=ugrid%cell_to_face_ghosted(iface2,local_id_dn)
+          if( (ugrid%face_to_cell_ghosted(1,face_id2)==ghosted_id_up.and. &
+               ugrid%face_to_cell_ghosted(2,face_id2)==ghosted_id_dn).or. &
+              (ugrid%face_to_cell_ghosted(1,face_id2)==ghosted_id_dn.and. &
+               ugrid%face_to_cell_ghosted(2,face_id2)==ghosted_id_up)) then
+            found=PETSC_TRUE
+            grid%fU2M(iface2,local_id_dn)=iface
+            grid%fM2U(iface)=face_id2
+            exit
+          endif
+        enddo
+        if(.not.found) then
+          option%io_buffer='2) UGRID face not found to match face in MIMETIC discretization'
+          call printErrMsg(option)
+        endif
+      endif
+
+    else if (conn%itype==BOUNDARY_CONNECTION_TYPE) then
+
+      ghosted_id_dn = conn%id_dn(face_id)
+      local_id_dn = grid%nG2L(ghosted_id_dn) ! Ghost to local mapping
+   
+      if (local_id_dn>0) then
+        aux_var => MFD%aux_vars(local_id_dn)
+        call MFDAuxAddFace(aux_var,option,iface)
+        grid%fG2L(iface)=local_id
+        grid%fL2G(local_id) = iface
+        local_id = local_id + 1
+      endif
+
+      nfaces=UCellGetNFaces(ugrid%cell_type(local_id_dn),option)
+      count=0
+      found=PETSC_FALSE
+      do iface2=1,nfaces
+        face_id2=ugrid%cell_to_face_ghosted(iface2,local_id_dn)
+        if (ugrid%face_to_cell_ghosted(2,face_id2) < 1) then
+          ! boundary face, since not connected to 2 cells
+          count=count+1
+          if(count>bnd_count(local_id_dn)) then
+            bnd_count(local_id_dn)=bnd_count(local_id_dn)+1
+            found=PETSC_TRUE
+            grid%fM2U(iface)=face_id2
+            grid%fU2M(iface2,local_id_dn)=iface
+            exit
+          endif
+        endif
+        if(found) exit
+      enddo
+      if(.not.found) then
+        option%io_buffer='UGRID boundary face not found to match face in MIMETIC discretization'
+        call printErrMsg(option)
+      endif
+
+    endif
+
+  enddo
+
+#endif
+
+end subroutine GridComputeCell2FaceForUGrid
+
+! ************************************************************************** !
+!> This routine sets up cell to face connection for an unstructure grid for
+!! MIMETIC discretization.
+!!
+!> @author
+!! Gautam Bisht, LBNL
+!!
+!! date: 03/29/13
+! ************************************************************************** !
+subroutine GridSetGlobalCell2FaceForUGrid(grid,MFD,DOF,option)
+
+  use MFD_Aux_module
+  use Option_module
+  use Unstructured_Cell_module
+ 
+  implicit none
+#include "finclude/petscvec.h"
+#include "finclude/petscvec.h90"
+#include "finclude/petscmat.h"
+#include "finclude/petscmat.h90"
+#include "finclude/petscdm.h"
+#include "finclude/petscdm.h90"
+#include "finclude/petscis.h"
+#include "finclude/petscis.h90"
+#include "finclude/petscviewer.h"
+#include "definitions.h"
+#include "finclude/petscsnes.h"
+#include "finclude/petscpc.h"
+
+  type(grid_type) :: grid
+  type(mfd_type), pointer :: MFD
+  PetscInt :: DOF
+  type(option_type) :: option
+
+  type(unstructured_grid_type),pointer :: ugrid
+  type(mfd_auxvar_type), pointer :: aux_var
+  type(connection_set_type), pointer :: conn
+
+  PetscInt :: local_id
+  PetscInt :: local_id_up
+  PetscInt :: local_id_dn
+  PetscInt :: ghosted_id
+  PetscInt :: ghosted_id_up
+  PetscInt :: ghosted_id_dn
+  PetscInt :: iface
+  PetscInt :: face_id
+  PetscInt :: ghost_face_id
+  PetscInt :: local_face_id
+  PetscInt :: nfaces
+  PetscInt :: ndof
+  PetscInt :: global_offset
+  PetscInt :: num_ghosted_upd
+  PetscInt :: icount,icell,stride,global_neigh_id,jcount
+  IS :: is_from
+  IS :: is_to
+  PetscViewer :: viewer
+  VecScatter :: scatter
+
+  PetscInt, pointer :: num_faces_cumm(:)
+
+  PetscInt, allocatable :: ghosted_ids(:)
+  PetscInt, allocatable :: indices_to(:)
+  PetscInt, allocatable :: indices_from(:)
+  PetscInt, allocatable :: int_array(:)
+  PetscInt, pointer :: int_ptr(:)
+
+  PetscScalar, pointer :: e2f_local_values(:)
+  PetscScalar, pointer :: e2n_local_values(:)
+  PetscScalar, pointer :: e2f_ghosted_values(:)
+  PetscScalar, pointer :: e2n_ghosted_values(:)
+  PetscScalar, pointer :: vec_ptr_e2f(:)
+  PetscScalar, pointer :: vec_ptr_e2n(:)
+  PetscScalar, pointer :: vec_ptr_e2f_gh(:)
+  PetscScalar, pointer :: vec_ptr_e2n_gh(:)
+
+  Vec :: e2f
+  Vec :: e2n
+  Vec :: e2f_ghosted
+  Vec :: e2n_ghosted
+
+  PetscErrorCode :: ierr
+
+  select case(DOF)
+    case(ONEDOF)
+      ndof = 1
+    case(NFLOWDOF)
+      ndof = option%nflowdof
+    case(NTRANDOF)
+      ndof = option%ntrandof
+  end select
+
+  MFD%ndof = ndof
+  ugrid => grid%unstructured_grid
+
+  ! Allocate memory
+  allocate(grid%fL2P(grid%nlmax_faces))
+  allocate(grid%fG2P(grid%ngmax_faces))
+
+  ! Set fL2P
+  global_offset = 0
+  do iface=1,grid%nlmax_faces
+    grid%fL2P(iface)=grid%global_faces_offset+grid%global_cell_offset+iface-1
+  enddo
+
+  global_offset = grid%global_faces_offset + grid%global_cell_offset
+
+  ! Create e2f and e2n
+  stride=MAX_FACE_PER_CELL
+  call VecCreate(option%mycomm,e2f,ierr)
+  call VecSetSizes(e2f,grid%nlmax*stride,PETSC_DECIDE,ierr)
+  call VecSetFromOptions(e2f,ierr)
+  call VecDuplicate(e2f,e2n,ierr)
+
+  call VecGetArrayF90(e2f,e2f_local_values,ierr)
+  call VecGetArrayF90(e2n,e2n_local_values,ierr)
+
+  e2f_local_values = 0
+  e2n_local_values = 0
+
+  do local_id = 1, grid%nlmax
+    aux_var => MFD%aux_vars(local_id)
+    do iface = 1, aux_var%numfaces
+      ghost_face_id = aux_var%face_id_gh(iface)
+      local_face_id = grid%fG2L(ghost_face_id)
+      conn => grid%faces(ghost_face_id)%conn_set_ptr
+      face_id = grid%faces(ghost_face_id)%id
+
+      if (conn%itype==INTERNAL_CONNECTION_TYPE) then
+
+        if(local_face_id>0) then
+          ! Face is a 'local' internal connection
+          e2f_local_values((local_id-1)*stride+iface)=grid%fL2P(local_face_id)+1
+
+          ghosted_id_up=conn%id_up(face_id)
+          ghosted_id_dn=conn%id_dn(face_id)
+          local_id_up=grid%nG2L(ghosted_id_up)
+          local_id_dn=grid%nG2L(ghosted_id_dn)
+
+          if (local_id_up==local_id) then
+            e2n_local_values((local_id-1)*stride+iface)=grid%nG2P(ghosted_id_dn)+1
+          else if (local_id_dn==local_id) then
+            e2n_local_values((local_id-1)*stride+iface)=grid%nG2P(ghosted_id_up)+1
+          endif
+
+        else if(local_face_id==0) then
+          ! Face is a 'non-local' internal connection
+          ghosted_id_up=conn%id_up(face_id)
+          ghosted_id_dn=conn%id_dn(face_id)
+          local_id_up=grid%nG2L(ghosted_id_up)
+          local_id_dn=grid%nG2L(ghosted_id_dn)
+
+          if (local_id_up==local_id) then
+            e2n_local_values((local_id-1)*stride+iface)=grid%nG2P(ghosted_id_dn)+1
+          else if (local_id_dn==local_id) then
+            e2n_local_values((local_id-1)*stride+iface)=grid%nG2P(ghosted_id_up)+1
+          endif
+        endif
+
+      else if (conn%itype == BOUNDARY_CONNECTION_TYPE) then
+
+        if (local_face_id>0) then
+          e2f_local_values((local_id-1)*stride+iface)=grid%fL2P(local_face_id)+1
+        endif
+
+        e2n_local_values((local_id-1)*stride+iface)=0
+      endif
+    enddo
+  enddo
+
+  call VecRestoreArrayF90(e2f,e2f_local_values,ierr)
+  call VecRestoreArrayF90(e2n,e2n_local_values,ierr)
+
+  ! Find PETSc ID of cells with which a non-local internal connections is shared
+  allocate(ghosted_ids(grid%ngmax_faces-grid%nlmax_faces))
+  ghosted_ids = 0
+
+  num_ghosted_upd = 0
+  do iface = 1, grid%ngmax_faces
+    conn => grid%faces(iface)%conn_set_ptr
+    face_id=grid%faces(iface)%id
+
+    if (conn%itype==INTERNAL_CONNECTION_TYPE) then
+      if (conn%local(face_id) == 0) then
+        ! Non-local internal connection
+        num_ghosted_upd = num_ghosted_upd + 1
+        ghosted_id_up = conn%id_up(face_id)
+        ghosted_id_dn = conn%id_dn(face_id)
+
+        ! Check if upwind or downwind is a ghost cell
+        if (grid%nG2L(ghosted_id_up)==0) then
+          ghosted_ids(num_ghosted_upd)=grid%nG2P(ghosted_id_up)+1 ! +1 since global indexes strats from 0
+        else if (grid%nG2L(ghosted_id_dn)==0) then
+          ghosted_ids(num_ghosted_upd)=grid%nG2P(ghosted_id_dn)+1
+        endif
+      endif
+    endif
+  enddo
+
+  ! Create sequential vectors with a stride
+  call VecCreateSeq(PETSC_COMM_SELF,num_ghosted_upd*stride,e2f_ghosted,ierr)
+  call VecCreateSeq(PETSC_COMM_SELF,num_ghosted_upd*stride,e2n_ghosted,ierr)
+         
+  allocate(indices_to(num_ghosted_upd))
+  allocate(indices_from(num_ghosted_upd))
+
+  do icount=1,num_ghosted_upd
+    indices_to(icount)=icount-1
+    indices_from(icount)=ghosted_ids(icount)-1
+  enddo
+
+  ! Create IS with a block size of 'stride'
+  call ISCreateBlock(option%mycomm,stride,num_ghosted_upd,indices_to, &
+                    PETSC_COPY_VALUES,is_to,ierr)
+  call ISCreateBlock(option%mycomm,stride,num_ghosted_upd,indices_from,&
+                    PETSC_COPY_VALUES,is_from,ierr)
+  deallocate(indices_from)
+  deallocate(indices_to)
+
+  ! Create scatter context
+  call VecScatterCreate(e2f,is_from,e2f_ghosted,is_to,scatter,ierr)
+  call ISDestroy(is_to,ierr)
+  call ISDestroy(is_from,ierr)
+
+  ! Scatter forward the e2f and e2n
+  call VecScatterBegin(scatter,e2f,e2f_ghosted,INSERT_VALUES,SCATTER_FORWARD,ierr)
+  call VecScatterEnd(scatter,e2f,e2f_ghosted,INSERT_VALUES,SCATTER_FORWARD,ierr)
+  call VecScatterBegin(scatter,e2n,e2n_ghosted,INSERT_VALUES,SCATTER_FORWARD,ierr)
+  call VecScatterEnd(scatter,e2n,e2n_ghosted,INSERT_VALUES,SCATTER_FORWARD,ierr)
+  call VecScatterDestroy(scatter,ierr)
+
+  call VecGetArrayF90(e2n_ghosted,vec_ptr_e2n_gh,ierr)
+  call VecGetArrayF90(e2f_ghosted,vec_ptr_e2f_gh,ierr)
+  call VecGetArrayF90(e2n,vec_ptr_e2n,ierr)
+  call VecGetArrayF90(e2f,vec_ptr_e2f,ierr)
+ 
+  jcount=0
+  do local_id=1,grid%nlmax
+    aux_var => MFD%aux_vars(local_id)
+    do iface=1,aux_var%numfaces
+      ghost_face_id=aux_var%face_id_gh(iface)
+      local_face_id=grid%fG2L(ghost_face_id)
+
+      if (local_face_id==0) then
+        conn => grid%faces(ghost_face_id)%conn_set_ptr
+        face_id=grid%faces(ghost_face_id)%id
+
+        ghosted_id_up=conn%id_up(face_id)
+        ghosted_id_dn=conn%id_dn(face_id)
+
+        local_id_up=grid%nG2L(ghosted_id_up) ! = zero for ghost nodes
+        local_id_dn=grid%nG2L(ghosted_id_dn) ! Ghost to local mapping
+
+        if (local_id==local_id_up) then
+          global_neigh_id=grid%nG2P(ghosted_id_dn)+1
+        else if (local_id==local_id_dn) then
+          global_neigh_id=grid%nG2P(ghosted_id_up)+1
+        endif
+
+        ! GB: Can avoid the this loop
+        !do jcount=1,num_ghosted_upd
+        !  if (ghosted_ids(jcount)==global_neigh_id) exit
+        !enddo
+        jcount=jcount+1
+
+        do face_id=1,stride
+          if (vec_ptr_e2n_gh((jcount-1)*stride+face_id)==grid%nG2P(grid%nL2G(local_id))+1) then
+            vec_ptr_e2f((local_id-1)*stride+iface)=vec_ptr_e2f_gh((jcount-1)*stride+face_id)
+            grid%fG2P(ghost_face_id)=int(vec_ptr_e2f_gh((jcount-1)*stride+face_id))-1
+            exit
+          endif
+        end do
+      else
+        grid%fG2P(ghost_face_id)=grid%fL2P(local_face_id)
+      endif
+    enddo
+  enddo
+
+  call VecRestoreArrayF90(e2n_ghosted,vec_ptr_e2n_gh,ierr)
+  call VecRestoreArrayF90(e2f_ghosted,vec_ptr_e2f_gh,ierr)
+  call VecDestroy(e2f_ghosted,ierr)
+  call VecDestroy(e2n_ghosted,ierr)
+
+  ! Count number of faces for all local cells
+  nfaces=0
+  allocate(num_faces_cumm(ugrid%nlmax))
+  num_faces_cumm=0
+  do local_id=1,ugrid%nlmax
+    nfaces=nfaces+UCellGetNFaces(ugrid%cell_type(local_id),option)
+    if(local_id<ugrid%nlmax) then
+      num_faces_cumm(local_id+1)=num_faces_cumm(local_id)+nfaces
+    endif
+  enddo
+
+  ! Allocate memory for grid%e2f and grid%e2n
+  call VecCreate(option%mycomm,grid%e2f,ierr)
+  call VecSetSizes(grid%e2f,nfaces,PETSC_DECIDE,ierr)
+  call VecSetFromOptions(grid%e2f,ierr)
+  call VecDuplicate(grid%e2f,grid%e2n,ierr)
+
+  call VecGetArrayF90(grid%e2f,e2f_local_values,ierr)
+  call VecGetArrayF90(grid%e2n,e2n_local_values,ierr)
+
+  nfaces=0
+  do local_id=1,grid%nlmax
+    do iface=1,UCellGetNFaces(ugrid%cell_type(local_id),option)
+      e2f_local_values(nfaces+iface)=vec_ptr_e2f((local_id-1)*stride+iface)
+      e2n_local_values(nfaces+iface)=vec_ptr_e2n((local_id-1)*stride+iface)
+    enddo
+    nfaces=nfaces+UCellGetNFaces(ugrid%cell_type(local_id),option)
+  enddo
+
+  call VecRestoreArrayF90(grid%e2f,e2f_local_values,ierr)
+  call VecRestoreArrayF90(grid%e2n,e2n_local_values,ierr)
+  call VecRestoreArrayF90(e2n,vec_ptr_e2n,ierr)
+  call VecRestoreArrayF90(e2f,vec_ptr_e2f,ierr)
+
+  call VecDestroy(e2f,ierr)
+  call VecDestroy(e2n,ierr)
+
+  call CreateMFDStruct4LP(grid,MFD,ndof,option)
+
+  deallocate(num_faces_cumm)
+
+
+end subroutine GridSetGlobalCell2FaceForUGrid
+
 end module Grid_module

src/pflotran/hydrostatic.F90

   dy_conn = 0.d0
   dz_conn = 0.d0
 
-  if (grid%itype==STRUCTURED_GRID_MIMETIC) then 
+  if (grid%itype==STRUCTURED_GRID_MIMETIC .or. &
+      grid%discretization_itype==UNSTRUCTURED_GRID_MIMETIC) then
       num_faces = coupler%numfaces_set
   else 
       num_faces = coupler%connection_set%num_connections
 
 
   do iconn=1, num_faces !geh: this should really be num_faces!
-    if (grid%itype==STRUCTURED_GRID_MIMETIC) then
+    if (grid%itype==STRUCTURED_GRID_MIMETIC .or. &
+        grid%discretization_itype==UNSTRUCTURED_GRID_MIMETIC) then
 #ifdef DASVYAT
       face_id_ghosted = coupler%faces_set(iconn)
       conn_set_ptr => grid%faces(face_id_ghosted)%conn_set_ptr
 
     if (associated(pressure_array)) then
       ipressure = idatum+int(dist_z/delta_z)
-      if (grid%itype==STRUCTURED_GRID_MIMETIC) then
+      if (grid%itype==STRUCTURED_GRID_MIMETIC.or. &
+          grid%discretization_itype==UNSTRUCTURED_GRID_MIMETIC) then
         dist_z_for_pressure = conn_set_ptr%cntr(3,conn_id)-(z(ipressure) + z_offset)
       else 
         dist_z_for_pressure = grid%z(ghosted_id)-dz_conn-(z(ipressure) + z_offset)
 
   enddo
 
-  if ((grid%itype==STRUCTURED_GRID_MIMETIC).and.(coupler%itype == INITIAL_COUPLER_TYPE)) then
+  if ((grid%itype==STRUCTURED_GRID_MIMETIC.or. &
+       grid%discretization_itype==UNSTRUCTURED_GRID_MIMETIC).and. &
+      (coupler%itype == INITIAL_COUPLER_TYPE)) then
      num_regions = coupler%region%num_cells
      do iconn = 1, num_regions
 

src/pflotran/init.F90

                              realization,ierr)
       case(RICHARDS_MODE)
         select case(realization%discretization%itype)
-          case(STRUCTURED_GRID_MIMETIC)
+          case(STRUCTURED_GRID_MIMETIC,UNSTRUCTURED_GRID_MIMETIC)
             call SNESSetFunction(flow_solver%snes,field%flow_r_faces, &
                                  RichardsResidualMFDLP, &
                                  realization,ierr)
                              THMCJacobian,realization,ierr)
       case(RICHARDS_MODE)
         select case(realization%discretization%itype)
-          case(STRUCTURED_GRID_MIMETIC)
+          case(STRUCTURED_GRID_MIMETIC,UNSTRUCTURED_GRID_MIMETIC)
             call SNESSetJacobian(flow_solver%snes,flow_solver%J,flow_solver%Jpre, &
                              RichardsJacobianMFDLP,realization,ierr)
           case default !sp 
     call GridComputeMinv(realization%discretization%grid, &
                          realization%discretization%stencil_width,option)
 
-
   call RealizationInitAllCouplerAuxVars(realization)
   if (option%ntrandof > 0) then
     call printMsg(option,"  Setting up TRAN Realization ")
   call DiscretizationReadRequiredCards(discretization,input,option)
   
   select case(discretization%itype)
-    case(STRUCTURED_GRID,UNSTRUCTURED_GRID,STRUCTURED_GRID_MIMETIC)
+    case(STRUCTURED_GRID,UNSTRUCTURED_GRID,STRUCTURED_GRID_MIMETIC,UNSTRUCTURED_GRID_MIMETIC)
       patch => PatchCreate()
       patch%grid => discretization%grid
       if (.not.associated(realization%level_list)) then
         call InputDefaultMsg(input,option,'Initial Condition name') 
         call CouplerRead(coupler,input,option)
         call RealizationAddCoupler(realization,coupler)
-        nullify(coupler)        
+        nullify(coupler)
       
 !....................
       case ('SOURCE_SINK')

src/pflotran/makefile

   MYFLAGS += -DSTORE_FLOWRATES
 endif
 
+ifdef mfd_ugrid
+  MYFLAGS += -DMFD_UGRID
+endif
+
 CFLAGS   = 
 FFLAGS   =
 CPPFLAGS = ${MYFLAGS}
 unstructured_cell.o : option.o utility.o 
 unstructured_explicit.o : geometry.o unstructured_grid_aux.o
 unstructured_grid.o : option.o connection.o unstructured_cell.o \
-                      unstructured_grid_aux.o geometry.o reaction_aux.o hdf5_aux.o
+                      unstructured_grid_aux.o geometry.o reaction_aux.o hdf5_aux.o \
+                      mfd_aux.o
 unstructured_grid_aux.o : unstructured_cell.o geometry.o option.o utility.o
 utility.o : string.o input.o option.o
 variables.o :

src/pflotran/mfd.F90

 ! ************************************************************************** !
 subroutine MFDCreateJacobianLP(grid, mfd_aux, mat_type, J, option)
 
- use Option_module
- use Grid_module
- use MFD_Aux_module
- use Connection_module
+  use Option_module
+  use Grid_module
+  use MFD_Aux_module
+  use Connection_module
+  use Unstructured_Grid_Aux_module
+  use Unstructured_Cell_module
 
   implicit none
 
   PetscInt :: icount, jcount, jface, loc_id_up, loc_id_dn
   PetscInt :: ndof_local
   PetscErrorCode :: ierr
+  PetscInt :: local_id
+  PetscInt :: ghosted_id
+  type(unstructured_grid_type),pointer :: ugrid
 
   allocate(d_nnz(grid%nlmax_faces + grid%nlmax))
   allocate(o_nnz(grid%nlmax_faces + grid%nlmax))
 
+  ugrid => grid%unstructured_grid
+
   d_nnz = 1 ! start 1 since diagonal connection to self
   o_nnz = 0
 
-  do icell = 1, grid%nlmax
-    aux_var => MFD_aux%aux_vars(icell)
-    do icount = 1, aux_var%numfaces
-      ghost_face_id = aux_var%face_id_gh(icount)
+  if (grid%itype == STRUCTURED_GRID_MIMETIC) then
+