Program TestMeshf

#include "finclude/petscdef.h"
Use petsc
Implicit NONE


#undef __FUNCT__
#define __FUNCT__ "main"
   PetscErrorCode       :: ierr
   Character(len=256)   :: filename,prefix,buffer
   Type(DM)             :: smesh,dmesh
   PetscBool            :: flg
   Type(PetscViewer)    :: viewer
   PetscInt             :: mode=0
   PetscMPIInt          :: rank,numproc
   PetscInt             :: labelsize
   PetscInt, Dimension(:), Pointer :: label
   PetscReal, Dimension(:,:), Pointer           :: Coord
   PetscInt, Dimension(:,:), Pointer            :: connect
   
   Call PetscInitialize(PETSC_NULL_CHARACTER,ierr);CHKERRQ(ierr)
   Call MPI_Comm_size(PETSC_COMM_WORLD,numproc,ierr)
   Call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr)

   Call PetscOptionsGetString(PETSC_NULL_CHARACTER,"-p",prefix,flg,ierr);CHKERRQ(ierr);
   Call PetscOptionsGetInt(PETSC_NULL_CHARACTER,"-mode",mode,flg,ierr);CHKERRQ(ierr);
   
   If (mode == 0) Then
#ifdef PETSC_HAVE_EXODUSII
      filename = Trim(prefix)//'.gen'
      buffer = "Reading " // Trim(filename) // "\n"
      Call PetscPrintf(PETSC_COMM_WORLD,buffer,ierr);CHKERRQ(ierr)
      Call DMMeshCreateExodus(PETSC_COMM_WORLD,filename,smesh,ierr);CHKERRQ(ierr)
      Call DMView(smesh,PETSC_VIEWER_STDOUT_WORLD,ierr);CHKERRQ(ierr)
      filename = Trim(prefix)//'_s.dat'
      Call PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_WRITE,viewer,ierr);CHKERRQ(ierr)
      Call DMView(smesh,viewer,ierr);CHKERRQ(ierr)
      Call PetscViewerDestroy(viewer,ierr);CHKERRQ(ierr)

      Call DMMeshDistribute(smesh,PETSC_NULL,dmesh,ierr);CHKERRQ(ierr)
      Call DMView(dmesh,PETSC_VIEWER_STDOUT_WORLD,ierr);CHKERRQ(ierr)
      filename = Trim(prefix)//'_d.dat'
      Call PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_WRITE,viewer,ierr);CHKERRQ(ierr)
      Call DMView(dmesh,viewer,ierr);CHKERRQ(ierr)
      Call PetscViewerDestroy(viewer,ierr);CHKERRQ(ierr)
      
      !!! Testing MeshGetLabelSize and MeshGetLabelIds
      buffer = "CellBlocks"
      Call DMMeshGetLabelSize(smesh,Buffer,labelsize,ierr);CHKERRQ(ierr);
      Allocate(label(labelsize))
      Call DMMeshGetLabelIds(smesh,buffer,label,ierr);CHKERRQ(ierr);
      Write(buffer, *) "Label size (smesh): ",labelsize, "\n"
      Call PetscSynchronizedPrintf(PETSC_COMM_WORLD,buffer,ierr);CHKERRQ(ierr);
      Call PetscSynchronizedFlush(PETSC_COMM_WORLD);CHKERRQ(ierr);
      Write(buffer, *) "Labels (smesh): ",label, "\n"
      Call PetscSynchronizedPrintf(PETSC_COMM_WORLD,buffer,ierr);CHKERRQ(ierr);
      Call PetscSynchronizedFlush(PETSC_COMM_WORLD);CHKERRQ(ierr);
      DeAllocate(label)
      
      
      !!! Testing MeshGetLabelSize and MeshGetLabelIds
      buffer = "CellBlocks"
      Call DMMeshGetLabelSize(dmesh,Buffer,labelsize,ierr);CHKERRQ(ierr);
      Allocate(label(labelsize))
      Call DMMeshGetLabelIds(dmesh,buffer,label,ierr);CHKERRQ(ierr);
      Write(buffer, *) "Label size (dmesh): ",labelsize, "\n"
      Call PetscSynchronizedPrintf(PETSC_COMM_WORLD,buffer,ierr);CHKERRQ(ierr);
      Call PetscSynchronizedFlush(PETSC_COMM_WORLD);CHKERRQ(ierr);
      Write(buffer, *) "Labels (dmesh): ",label, "\n"
      Call PetscSynchronizedPrintf(PETSC_COMM_WORLD,buffer,ierr);CHKERRQ(ierr);
      Call PetscSynchronizedFlush(PETSC_COMM_WORLD);CHKERRQ(ierr);
      DeAllocate(label)

      Call DMDestroy(smesh,ierr);CHKERRQ(ierr)
      Call DMDestroy(dmesh,ierr);CHKERRQ(ierr)
#else
	  SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP, 'Reading EXO mesh requires EXODUS II support in PETSc',ierr)
#endif 
   Else
      Call DMMeshCreate(PETSC_COMM_WORLD,smesh,ierr);CHKERRQ(ierr)
      Call DMMeshCreate(PETSC_COMM_WORLD,dmesh,ierr);CHKERRQ(ierr)
      If (numproc == 1) Then
         filename = Trim(prefix)//'_s.dat'
      Else
         filename = Trim(prefix)//'_d.dat'
      End If
      buffer = "Reading " // Trim(filename) // "\n"
      Call PetscPrintf(PETSC_COMM_WORLD,buffer,ierr);CHKERRQ(ierr)
      Call PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,viewer,ierr);CHKERRQ(ierr)
      Call DMMeshLoad(viewer,smesh,ierr);CHKERRQ(ierr)
      Call PetscViewerDestroy(viewer,ierr);CHKERRQ(ierr)
      Call DMView(smesh,PETSC_VIEWER_STDOUT_WORLD,ierr);CHKERRQ(ierr)
      
      Call DMMeshGetCoordinatesF90(smesh,Coord,ierr);CHKERRQ(ierr)
      Write(*,*) 'Coordinates: '
      Write(*,*) Coord
      Call DMMeshRestoreCoordinatesF90(smesh,Coord,ierr);CHKERRQ(ierr)
      
      Call DMMeshGetElementsF90(smesh,connect,ierr);CHKERRQ(ierr)
      write(*,*) 'Connect: '
      write(*,*) connect
      Call DMMeshRestoreElementsF90(smesh,connect,ierr);CHKERRQ(ierr)
      Call DMDestroy(smesh,ierr);CHKERRQ(ierr)
      Call PetscFinalize(ierr)
      
      Call DMDestroy(smesh,ierr);CHKERRQ(ierr)
  End If  
  Call PetscFinalize(ierr)
End Program TestMeshf