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:       PetscViewer    viewer;
  7:       FILE          *f;
  8:       PetscReal     *coords;
  9:       PetscReal     *normals;
 10:       PetscInt      *verts;
 11:       PetscInt       commRank;
 12:       char           buf[2048];

 15:       MPI_Comm_rank(comm, &commRank);
 16:       if (commRank != 0) return;
 17:       PetscViewerCreate(PETSC_COMM_SELF, &viewer);
 18:       PetscViewerSetType(viewer, PETSCVIEWERASCII);
 19:       PetscViewerFileSetMode(viewer, FILE_MODE_READ);
 20:       PetscViewerFileSetName(viewer, filename.c_str());
 21:       PetscViewerASCIIGetPointer(viewer, &f);
 22:       // Read header
 23:       if (fgets(buf, 2048, f) == NULL) throw ALE::Exception("File ended prematurely");
 24:       // Number of vertices
 25:       const char *x = strtok(buf, " ");
 26:       numVertices = atoi(x);
 27:       // Number of elements
 28:       x = strtok(NULL, " ");
 29:       numElements = atoi(x);
 30:       PetscMalloc(numVertices*(dim+1) * sizeof(PetscReal), &coords);CHKERRXX(ierr);
 31:       PetscMalloc(numVertices*(dim+1) * sizeof(PetscReal), &normals);CHKERRXX(ierr);
 32:       for(PetscInt i = 0; i < numVertices; ++i) {
 33:         if (fgets(buf, 2048, f) == NULL) throw ALE::Exception("File ended prematurely");
 34:         x = strtok(buf, " ");
 35:         for(int c = 0; c < dim+1; c++) {
 36:           coords[i*(dim+1)+c] = atof(x);
 37:           x = strtok(NULL, " ");
 38:         }
 39:         for(int c = 0; c < dim+1; c++) {
 40:           normals[i*(dim+1)+c] = atof(x);
 41:           x = strtok(NULL, " ");
 42:         }
 43:         // Ignore ???
 44:         x = strtok(NULL, " ");
 45:         // Ignore ???
 46:         x = strtok(NULL, " ");
 47:         // Ignore ???
 48:         x = strtok(NULL, " ");
 49:       }
 50:       *coordinates = coords;
 51:       *faceNormals = normals;
 52:       PetscMalloc(numElements*numCorners * sizeof(PetscInt), &verts);CHKERRXX(ierr);
 53:       for(PetscInt i = 0; i < numElements; ++i) {
 54:         if (fgets(buf, 2048, f) == NULL) throw ALE::Exception("File ended prematurely");
 55:         x = strtok(buf, " ");
 56:         for(int c = 0; c < numCorners; c++) {
 57:           verts[i*numCorners+c] = atoi(x) - 1;
 58:           x = strtok(NULL, " ");
 59:         }
 60:         // Ignore ???
 61:         x = strtok(NULL, " ");
 62:         // Ignore ???
 63:         x = strtok(NULL, " ");
 64:       }
 65:       *vertices = verts;
 66:       PetscViewerDestroy(&viewer);
 67:     };
 68:     Obj<Builder::Mesh> Builder::readMesh(MPI_Comm comm, const int dim, const std::string& filename, const bool interpolate = false, const int debug = 0) {
 69:       typedef ALE::Mesh<PetscInt,PetscScalar> FlexMesh;
 70:       Obj<Mesh>       mesh  = new Mesh(comm, dim, debug);
 71:       Obj<sieve_type> sieve = new sieve_type(comm, debug);
 72:       Obj<FlexMesh>             m = new FlexMesh(comm, dim, debug);
 73:       Obj<FlexMesh::sieve_type> s = new FlexMesh::sieve_type(comm, debug);
 74:       std::map<Mesh::point_type,Mesh::point_type> renumbering;
 75:       int       *cells;
 76:       PetscReal *coordinates, *normals;
 77:       int        numCells = 0, numVertices = 0, numCorners = dim+1;

 80:       mesh->setSieve(sieve);
 81:       Builder::readInpFile(comm, filename, dim, numCorners, numCells, &cells, numVertices, &coordinates, &normals);
 82:       ALE::SieveBuilder<FlexMesh>::buildTopology(s, dim, numCells, cells, numVertices, interpolate, numCorners, -1, m->getArrowSection("orientation"));
 83:       m->setSieve(s);
 84:       m->stratify();
 85:       ALE::SieveBuilder<FlexMesh>::buildCoordinates(m, dim+1, coordinates);
 86:       PetscFree(cells);CHKERRXX(ierr);
 87:       PetscFree(coordinates);CHKERRXX(ierr);
 88:       PetscFree(normals);CHKERRXX(ierr);
 89:       ALE::ISieveConverter::convertMesh(*m, *mesh, renumbering, false);
 90:       return mesh;
 91:     };
 92:   }
 93: }