1:  #include <petsc/private/dmpleximpl.h>
  2:  #include <petsc/private/isimpl.h>
  3:  #include <petsc/private/vecimpl.h>
  4: #include <petscviewerhdf5.h>

  6: #if defined(PETSC_HAVE_HDF5)
  7: PetscErrorCode DMPlexLoad_HDF5_Xdmf_Internal(DM dm, PetscViewer viewer)
  8: {
  9:   Vec             coordinates;
 10:   IS              cells;
 11:   PetscInt        dim, spatialDim, N, numCells, numVertices, numCorners, bs;
 12:   PetscMPIInt     rank;
 13:   MPI_Comm        comm;
 14:   PetscErrorCode  ierr;

 17:   PetscObjectGetComm((PetscObject)dm, &comm);
 18:   MPI_Comm_rank(comm, &rank);

 20:   /* Read topology */
 21:   PetscViewerHDF5ReadAttribute(viewer, "/viz/topology/cells", "cell_dim", PETSC_INT, (void *) &dim);
 22:   PetscViewerHDF5ReadAttribute(viewer, "/viz/topology/cells", "cell_corners", PETSC_INT, (void *) &numCorners);
 23:   PetscViewerHDF5PushGroup(viewer, "/viz/topology");
 24:   ISCreate(comm, &cells);
 25:   PetscObjectSetName((PetscObject) cells, "cells");
 26:   ISLoad(cells, viewer);
 27:   ISGetLocalSize(cells, &numCells);
 28:   ISGetBlockSize(cells, &bs);
 29:   if (bs != numCorners) SETERRQ(comm, PETSC_ERR_FILE_UNEXPECTED, "cell_corners and 2-dimension of cells not equal");
 30:   PetscViewerHDF5PopGroup(viewer);
 31:   numCells /= numCorners;

 33:   /* Read geometry */
 34:   PetscViewerHDF5PushGroup(viewer, "/geometry");
 35:   VecCreate(comm, &coordinates);
 36:   PetscObjectSetName((PetscObject) coordinates, "vertices");
 37:   PetscViewerHDF5ReadSizes(viewer, "vertices", &spatialDim, &N);
 38:   VecSetSizes(coordinates, PETSC_DECIDE, N);
 39:   VecSetBlockSize(coordinates, spatialDim);
 40:   VecLoad(coordinates, viewer);
 41:   VecGetLocalSize(coordinates, &numVertices);
 42:   VecGetBlockSize(coordinates, &spatialDim);
 43:   PetscViewerHDF5PopGroup(viewer);
 44:   numVertices /= spatialDim;

 46: #if defined(PETSC_USE_DEBUG)
 47:   /* Check that maximum index referred in cells is in line with global number of vertices */
 48:   {
 49:     PetscInt max1, max2;
 50:     ISGetMinMax(cells, NULL, &max1);
 51:     VecGetSize(coordinates, &max2);
 52:     max2 /= spatialDim; max2--;
 53:     if (max1 > max2) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "maximum index in cells = %d > %d = total number of vertices - 1", max1, max2);
 54:   }
 55: #endif

 57:   {
 58:     const PetscScalar *coordinates_arr;
 59:     PetscReal *coordinates_arr_real;
 60:     const PetscInt *cells_arr;
 61:     int *cells_arr_int;
 62:     PetscSF sfVert=NULL;
 63: #if defined(PETSC_USE_64BIT_INDICES) || defined(PETSC_USE_COMPLEX)
 64:     PetscInt i;
 65: #endif

 67:     VecGetArrayRead(coordinates, &coordinates_arr);
 68:     ISGetIndices(cells, &cells_arr);

 70: #if defined(PETSC_USE_64BIT_INDICES)
 71:     /* convert to 32-bit integers if PetscInt is 64-bit */
 72:     /*TODO More systematic would be to change all the function arguments to PetscInt */
 73:     PetscMalloc1(numCells*numCorners, &cells_arr_int);
 74:     for (i = 0; i < numCells*numCorners; ++i) {
 75:       PetscMPIIntCast(cells_arr[i], &cells_arr_int[i]);
 76:     }
 77: #else
 78:     cells_arr_int = (int*)cells_arr;
 79: #endif

 81: #if defined(PETSC_USE_COMPLEX)
 82:     /* convert to real numbers if PetscScalar is complex */
 83:     /*TODO More systematic would be to change all the function arguments to PetscScalar */
 84:     PetscMalloc1(numVertices*spatialDim, &coordinates_arr_real);
 85:     for (i = 0; i < numVertices*spatialDim; ++i) {
 86:       coordinates_arr_real[i] = PetscRealPart(coordinates_arr[i]);
 87: #if defined(PETSC_USE_DEBUG)
 88:       if (PetscImaginaryPart(coordinates_arr[i])) {
 89:         SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Vector of coordinates contains complex numbers but only real vectors are currently supported.");
 90:       }
 91: #endif  /* defined(PETSC_USE_DEBUG) */
 92:     }
 93: #else
 94:     coordinates_arr_real = (PetscReal*)coordinates_arr;
 95: #endif

 97:     DMSetDimension(dm, dim);
 98:     DMPlexBuildFromCellList_Parallel_Internal(dm, spatialDim, numCells, numVertices, numCorners, cells_arr_int, PETSC_TRUE, &sfVert);
 99:     DMPlexBuildCoordinates_Parallel_Internal( dm, spatialDim, numCells, numVertices, sfVert, coordinates_arr_real);
100:     VecRestoreArrayRead(coordinates, &coordinates_arr);
101:     ISRestoreIndices(cells, &cells_arr);
102:     PetscSFDestroy(&sfVert);
103: #if defined(PETSC_USE_64BIT_INDICES)
104:     PetscFree(cells_arr_int);
105: #endif
106: #if defined(PETSC_USE_COMPLEX)
107:     PetscFree(coordinates_arr_real);
108: #endif
109:   }
110:   ISDestroy(&cells);
111:   VecDestroy(&coordinates);

113:   /* scale coordinates - unlike in DMPlexLoad_HDF5_Internal, this can only be done after DM is populated */
114:   {
115:     PetscReal lengthScale;

117:     DMPlexGetScale(dm, PETSC_UNIT_LENGTH, &lengthScale);
118:     DMGetCoordinates(dm, &coordinates);
119:     VecScale(coordinates, 1.0/lengthScale);
120:   }

122:   /* Read Labels */
123:   /* TODO: this probably does not work as elements get permuted */
124:   /* DMPlexLoadLabels_HDF5_Internal(dm, viewer); */
125:   return(0);
126: }
127: #endif