1: #include <petscdmmesh_formats.hh>   /*I      "petscmesh.h"   I*/

  3: namespace ALE {
  4: namespace Bardhan {
  5: void Builder::readInpFile(MPI_Comm comm, const std::string& filename, const int dim, const int numCorners, int& numElements, int *vertices[], int& numVertices, PetscReal *coordinates[], PetscReal *faceNormals[])
  6: {
  7:   PetscViewer    viewer;
  8:   FILE           *f;
  9:   PetscReal      *coords;
 10:   PetscReal      *normals;
 11:   PetscInt       *verts;
 12:   PetscInt       commRank;
 13:   char           buf[2048];

 16:   MPI_Comm_rank(comm, &commRank);
 17:   if (commRank != 0) return;
 18:   PetscViewerCreate(PETSC_COMM_SELF, &viewer);
 19:   PetscViewerSetType(viewer, PETSCVIEWERASCII);
 20:   PetscViewerFileSetMode(viewer, FILE_MODE_READ);
 21:   PetscViewerFileSetName(viewer, filename.c_str());
 22:   PetscViewerASCIIGetPointer(viewer, &f);
 23:   // Read header
 24:   if (!fgets(buf, 2048, f)) throw ALE::Exception("File ended prematurely");
 25:   // Number of vertices
 26:   const char *x = strtok(buf, " ");
 27:   numVertices = atoi(x);
 28:   // Number of elements
 29:   x           = strtok(NULL, " ");
 30:   numElements = atoi(x);

 32:   PetscMalloc(numVertices*(dim+1) * sizeof(PetscReal), &coords);CHKERRXX(ierr);
 33:   PetscMalloc(numVertices*(dim+1) * sizeof(PetscReal), &normals);CHKERRXX(ierr);

 35:   for (PetscInt i = 0; i < numVertices; ++i) {
 36:     if (!fgets(buf, 2048, f)) throw ALE::Exception("File ended prematurely");
 37:     x = strtok(buf, " ");
 38:     for (int c = 0; c < dim+1; c++) {
 39:       coords[i*(dim+1)+c] = atof(x);

 41:       x = strtok(NULL, " ");
 42:     }
 43:     for (int c = 0; c < dim+1; c++) {
 44:       normals[i*(dim+1)+c] = atof(x);

 46:       x = strtok(NULL, " ");
 47:     }
 48:     // Ignore ???
 49:     x = strtok(NULL, " ");
 50:     // Ignore ???
 51:     x = strtok(NULL, " ");
 52:     // Ignore ???
 53:     x = strtok(NULL, " ");
 54:   }
 55:   *coordinates = coords;
 56:   *faceNormals = normals;
 57:   PetscMalloc(numElements*numCorners * sizeof(PetscInt), &verts);CHKERRXX(ierr);
 58:   for (PetscInt i = 0; i < numElements; ++i) {
 59:     if (!fgets(buf, 2048, f)) throw ALE::Exception("File ended prematurely");
 60:     x = strtok(buf, " ");
 61:     for (int c = 0; c < numCorners; c++) {
 62:       verts[i*numCorners+c] = atoi(x) - 1;

 64:       x = strtok(NULL, " ");
 65:     }
 66:     // Ignore ???
 67:     x = strtok(NULL, " ");
 68:     // Ignore ???
 69:     x = strtok(NULL, " ");
 70:   }
 71:   *vertices = verts;
 72:   PetscViewerDestroy(&viewer);
 73: };
 74: Obj<Builder::Mesh> Builder::readMesh(MPI_Comm comm, const int dim, const std::string& filename, const bool interpolate = false, const int debug = 0)
 75: {
 76:   typedef ALE::Mesh<PetscInt,PetscScalar> FlexMesh;
 77:   Obj<Mesh>                                   mesh  = new Mesh(comm, dim, debug);
 78:   Obj<sieve_type>                             sieve = new sieve_type(comm, debug);
 79:   Obj<FlexMesh>                               m     = new FlexMesh(comm, dim, debug);
 80:   Obj<FlexMesh::sieve_type>                   s     = new FlexMesh::sieve_type(comm, debug);
 81:   std::map<Mesh::point_type,Mesh::point_type> renumbering;
 82:   int                                         *cells;
 83:   PetscReal                                   *coordinates, *normals;
 84:   int                                         numCells = 0, numVertices = 0, numCorners = dim+1;
 85:   PetscErrorCode                              ierr;

 87:   mesh->setSieve(sieve);
 88:   Builder::readInpFile(comm, filename, dim, numCorners, numCells, &cells, numVertices, &coordinates, &normals);
 89:   ALE::SieveBuilder<FlexMesh>::buildTopology(s, dim, numCells, cells, numVertices, interpolate, numCorners, -1, m->getArrowSection("orientation"));
 90:   m->setSieve(s);
 91:   m->stratify();
 92:   ALE::SieveBuilder<FlexMesh>::buildCoordinates(m, dim+1, coordinates);
 93:   PetscFree(cells);CHKERRXX(ierr);
 94:   PetscFree(coordinates);CHKERRXX(ierr);
 95:   PetscFree(normals);CHKERRXX(ierr);
 96:   ALE::ISieveConverter::convertMesh(*m, *mesh, renumbering, false);
 97:   return mesh;
 98: };
 99: }
100: }