Actual source code: triconvert.c

      2: /*
      3:       Converts triangulated grid data file.node and file.ele generated
      4:    by the Triangle code of Shewchuk to a PETSc AOData base.

      6: */

     8:  #include petscao.h
     9:  #include petscbt.h

     11: int main(int argc,char **args)
     12: {
     13:   int        ierr,nvertex,dim,nstuff,nbound,dummy,bound,i,ncell,*cell,*cell_edge,*edge_vertex;
     14:   int        nedge,j,k,*edge_cell,*cell_cell,nc,ncp;
     15:   int        i0,i1;
     16:   int        *shift0,*shift1;
     17:   char       filebase[PETSC_MAX_PATH_LEN],filename[PETSC_MAX_PATH_LEN];
     18:   FILE       *file;
     19:   AOData     ao;
     20:   PetscReal  *vertex;
     21:   PetscBT    vertex_boundary;
     22:   PetscTruth flag;

     24:   PetscInitialize(&argc,&args,0,0);
     25:   PetscMemzero(filename,PETSC_MAX_PATH_LEN*sizeof(char));

     27:   PetscOptionsGetString(0,"-f",filebase,PETSC_MAX_PATH_LEN-1,&flag);
     28:   if (!flag) {
     29:     SETERRQ(1,"Must provide filebase name with -f");
     30:   }

     32:   /*
     33:      Create empty database
     34:   */
     35:   AODataCreateBasic(PETSC_COMM_SELF,&ao);

     37:   /* -----------------------------------------------------------------------------------
     38:        Read in vertex information
     39:   */
     40:   PetscStrcpy(filename,filebase);
     41:   PetscStrcat(filename,".node");
     42:   file = fopen(filename,"r");
     43:   if (!file) {
     44:     SETERRQ1(1,"Unable to open node file: %s",filename);
     45:   }
     46:   fscanf(file,"%d %d %d %d\n",&nvertex,&dim,&nstuff,&nbound);
     47:   if (dim != 2) {
     48:     SETERRQ(1,"Triangulation is not in two dimensions");
     49:   }
     50:   PetscMalloc(2*nvertex*sizeof(PetscReal),&vertex);
     51:   PetscBTCreate(nvertex,vertex_boundary);

     53:   if (nstuff == 1) {
     54:     double v0,v1,ddummy;
     55:     for (i=0; i<nvertex; i++) {
     56:       fscanf(file,"%d %le %le %le %d\n",&dummy,&v0,&v1,&ddummy,&bound);
     57:       vertex[2*i]   = v0;
     58:       vertex[2*i+1] = v1;
     59:       if (bound) {PetscBTSet(vertex_boundary,i);}
     60:     }
     61:   } else  if (nstuff == 0) {
     62:     double v0,v1;
     63:     for (i=0; i<nvertex; i++) {
     64:       fscanf(file,"%d %le %le %d\n",&dummy,&v0,&v1,&bound);
     65:       vertex[2*i]   = v0;
     66:       vertex[2*i+1] = v1;
     67:       if (bound) {PetscBTSet(vertex_boundary,i);}
     68:     }
     69:   } else SETERRQ(1,"No support yet for that number of vertex quantities");
     70:   fclose(file);

     72:   /*
     73:      Put vertex into database
     74:   */
     75:   AODataKeyAdd(ao,"vertex",nvertex,nvertex);
     76:   AODataSegmentAdd(ao,"vertex","values",2,nvertex,0,vertex,PETSC_DOUBLE);
     77:   AODataSegmentAdd(ao,"vertex","boundary",1,nvertex,0,vertex_boundary,PETSC_LOGICAL);
     78:   PetscBTDestroy(vertex_boundary);

     80:   /* -----------------------------------------------------------------------------------
     81:       Read in triangle information
     82:   */
     83:   PetscStrcpy(filename,filebase);
     84:   PetscStrcat(filename,".ele");
     85:   file = fopen(filename,"r");
     86:   if (!file) {
     87:     SETERRQ(1,"Unable to open element file");
     88:   }
     89:   fscanf(file,"%d %d %d\n",&ncell,&nc,&nstuff);ncp = nc;

     91:   PetscMalloc(nc*ncell*sizeof(int),&cell);
     92:   if (nstuff == 0) {
     93:     if (nc == 3) {
     94:       for (i=0; i<ncell; i++) {
     95:         fscanf(file,"%d %d %d %d\n",&dummy,cell+3*i,cell+3*i+1,cell+3*i+2);
     96:       }
     97:     } else if (nc == 6) {
     98:       for (i=0; i<ncell; i++) {
     99:         fscanf(file,"%d %d %d %d %d %d %d\n",&dummy,cell+6*i,cell+6*i+2,cell+6*i+4,
    100:                cell+6*i+3,cell+6*i+5,cell+6*i+1);
    101:       }
    102:     }
    103:   } else {
    104:     SETERRQ(1,"No support yet for that number of element quantities");
    105:   }
    106:   fclose(file);
    107:   for (i=0; i<nc*ncell; i++) {
    108:     cell[i]--;    /* shift because triangle starts at 1, not 0 */
    109:   }

    111:   /*
    112:      Put cell into database
    113:   */
    114:   AODataKeyAdd(ao,"cell",ncell,ncell);
    115:   AODataSegmentAdd(ao,"cell","vertex",nc,ncell,0,cell,PETSC_INT);

    117:   PetscMalloc(nc*ncell*sizeof(int),&cell_edge);
    118:   PetscMalloc(2*nc*ncell*sizeof(int),&edge_cell);
    119:   PetscMalloc(2*nc*ncell*sizeof(int),&edge_vertex);
    120:   PetscMalloc(3*ncell*sizeof(int),&cell_cell);

    122:   /*
    123:       Determine edges
    124:   */
    125:   PetscMalloc(nc*sizeof(int),&shift0);
    126:   PetscMalloc(nc*sizeof(int),&shift1);
    127:   for (i=0; i<nc; i++) {
    128:     shift0[i] = i;
    129:     shift1[i] = (i + 1) % nc;
    130:   }

    133:   nedge = 0;
    134:   for (i=0; i<ncell; i++) {
    135:     for (k=0; k<nc; k++) {
    136:       i0 = cell[nc*i+shift0[k]];
    137:       i1 = cell[nc*i+shift1[k]];
    138:       for (j=0; j<nedge; j++) {
    139:         if ((i0 == edge_vertex[2*j+1] && i1 == edge_vertex[2*j]) ||
    140:             (i1 == edge_vertex[2*j+1] && i0 == edge_vertex[2*j])) {
    141:           cell_edge[nc*i+k]   = j;
    142:           edge_cell[2*j+1]   = i;
    143:           goto found;
    144:         }
    145:       }
    146:       /*
    147:            Add a new edge to the list
    148:       */
    149:       edge_cell[2*nedge]          = i;
    150:       edge_cell[2*nedge+1]        = -1;
    151:       edge_vertex[2*nedge]        = i0;
    152:       edge_vertex[2*nedge+1]      = i1;
    153:       cell_edge[nc*i+k]           = nedge;
    154:       nedge++;
    155:
    156:       found:;
    157:     }
    158:   }

    160:   /*
    161:        Determine cell neighbors
    162:   */
    163:   PetscStrcpy(filename,filebase);
    164:   PetscStrcat(filename,".neigh");
    165:   file = fopen(filename,"r");
    166:   if (file) {
    167:     fscanf(file,"%d %d\n",&ncell,&nc);
    168:     if (nc != 3) SETERRQ(PETSC_ERR_SUP,"Can only handle three neighbors");
    169:     for (i=0; i<ncell; i++) {
    170:       fscanf(file,"%d %d %d %d\n",&dummy,cell_cell+3*i+1,cell_cell+3*i+2,cell_cell+3*i);
    171:     }
    172:     for (i=0; i<3*ncell; i++) {
    173:       cell_cell[i]--;    /* shift because triangle starts at 1, not 0 */
    174:     }
    175:     fclose(file);

    177:   } else { /* no neighbor list file given, generate manually only works for nc == 3 */
    178:     if (nc != 3) SETERRQ(PETSC_ERR_SUP,"No neighbor file given and cannot determine neighbors");

    180:     for (i=0; i<ncell; i++) {
    181:       for (k=0; k<3; k++) {
    182:         i0 = cell_edge[3*i+k];
    183:         if (edge_cell[2*i0] != i) cell_cell[3*i+k] = edge_cell[2*i0];
    184:         else                      cell_cell[3*i+k] = edge_cell[2*i0+1];
    185:       }
    186:     }
    187:   }

    189:   AODataSegmentAdd(ao,"cell","cell",3,ncell,0,cell_cell,PETSC_INT);
    190:   AODataSegmentAdd(ao,"cell","edge",nc,ncell,0,cell_edge,PETSC_INT);

    192:   AODataKeyAdd(ao,"edge",nedge,nedge);
    193:   AODataSegmentAdd(ao,"edge","vertex",2,nedge,0,edge_vertex,PETSC_INT);

    195:   PetscFree(vertex);
    196:   PetscFree(cell);
    197:   PetscFree(shift0);
    198:   PetscFree(shift1);

    200:   /*
    201:       Add information about cell shape and element type to database
    202:   */
    203:   AODataKeyAdd(ao,"info",1,1);
    204:   AODataSegmentAdd(ao,"info","shape",10,1,0,(void*)"triangular",PETSC_CHAR);
    205:   if (ncp == 3) {
    206:     AODataSegmentAdd(ao,"info","element",6,1,0,(void*)"linear",PETSC_CHAR);
    207:   } else if (ncp == 6) {
    208:     AODataSegmentAdd(ao,"info","element",13,1,0,(void*)"quadratic",PETSC_CHAR);
    209:   }

    211:   /*  AODataView(ao,0); */

    213:   { PetscViewer binary;
    214:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,filebase,PETSC_FILE_CREATE,&binary);
    215:   AODataView(ao,binary);
    216:   PetscViewerDestroy(binary);
    217:   }
    218:
    219:   return 0;
    220: }