Actual source code: inpututils.c

      1: #define PETSCDM_DLL

      3: /*
      4:        Utilities for inputing, creating and managing simple two dimensional grids
      5: */

     7:  #include src/dm/ao/aoimpl.h
     8:  #include petscbt.h
     9:  #include petscdraw.h

     11: /*
     12:     cell_n        - number of cells
     13:     max_cell      - maximum space allocated for cell
     14:     cell_vertex   - vertices of each cell
     15:     cell_edge     - edges of the cell
     16:     cell_cell     - neighbors of cell
     17:     vertex_n      - number of vertices
     18:     vertex_max    - maximum space allocated for vertices
     19:     x,y           - vertex coordinates

     21:     xmin,ymin,xmax,ymax - bounding box of grid

     23:     edge_n        - total edges in the grid
     24:     edge_vertex   - vertex of all edges
     25:     edge_max      - maximum space allocated for edge
     26:     edge_cell     - two neighbor cells who share edge

     28:     vertex_boundary - indicates for each vertex if it is a boundary

     30: */

     36: PetscErrorCode PETSCDM_DLLEXPORT AOData2dGridToAOData(AOData2dGrid agrid,AOData *ao)
     37: {
     39:   PetscInt       *keys,nmax,i;
     40:   AOData         aodata;

     43:   /*
     44:       Create the database
     45:   */
     46:   nmax = PetscMax(agrid->cell_n,agrid->vertex_n);
     47:   nmax = PetscMax(nmax,agrid->edge_n);
     48:   PetscMalloc(nmax*sizeof(PetscInt),&keys);
     49:   for (i=0; i<nmax; i++) {
     50:     keys[i] = i;
     51:   }
     52:   AODataCreateBasic(PETSC_COMM_WORLD,&aodata);
     53:     AODataKeyAdd(aodata,"cell",PETSC_DECIDE,agrid->cell_n);
     54:       AODataSegmentAdd(aodata,"cell","cell",4,agrid->cell_n,keys,agrid->cell_cell,PETSC_INT);
     55:       AODataSegmentAdd(aodata,"cell","vertex",4,agrid->cell_n,keys,agrid->cell_vertex,PETSC_INT);
     56:       AODataSegmentAdd(aodata,"cell","edge",4,agrid->cell_n,keys,agrid->cell_edge,PETSC_INT);
     57:     AODataKeyAdd(aodata,"edge",PETSC_DECIDE,agrid->edge_n);
     58:       AODataSegmentAdd(aodata,"edge","vertex",2,agrid->edge_n,keys,agrid->edge_vertex,PETSC_INT);
     59:       AODataSegmentAdd(aodata,"edge","cell",2,agrid->edge_n,keys,agrid->edge_cell,PETSC_INT);
     60:     AODataKeyAdd(aodata,"vertex",PETSC_DECIDE,agrid->vertex_n);
     61:       AODataSegmentAdd(aodata,"vertex","values",2,agrid->vertex_n,keys,agrid->vertex,PETSC_DOUBLE);
     62:       AODataSegmentAdd(aodata,"vertex","boundary",1,agrid->vertex_n,keys,agrid->vertex_boundary,PETSC_LOGICAL);
     63:   PetscFree(keys);
     64:   *ao = aodata;
     65:   return(0);
     66: }

     70: /*
     71:        User input the cell by drawing them one at a time
     72: */
     73: PetscErrorCode PETSCDM_DLLEXPORT AOData2dGridInput(AOData2dGrid agrid,PetscDraw draw)
     74: {
     75:   PetscDraw       popup;                           /* help window */
     76:   PetscDrawButton button;                          /* mouse button pressed */
     78:   PetscInt             cn, *cell;
     79:   PetscReal       *vertex,cx,cy;
     80:   char            title[120];

     83:   agrid->cell_max      = 500;
     84:   agrid->cell_n        = 0;
     85:   agrid->vertex_max    = 500;
     86:   agrid->vertex_n      = 0;
     87:   agrid->xmin          = PETSC_MAX;
     88:   agrid->xmax          = PETSC_MIN;
     89:   agrid->ymin          = PETSC_MAX;
     90:   agrid->ymax          = PETSC_MIN;

     92:   /*
     93:      Allocate large arrays to hold the nodes and cellrilateral lists
     94:   */
     95:   PetscMalloc(2*agrid->vertex_max*sizeof(PetscReal),&agrid->vertex);
     96:   vertex = agrid->vertex;
     97:   PetscMalloc(4*agrid->cell_max*sizeof(PetscInt),&agrid->cell_vertex);
     98:   cell   = agrid->cell_vertex;

    101:   /*
    102:      Open help window and enter helpful messages
    103:   */
    104:   PetscDrawGetPopup(draw,&popup);
    105:   PetscDrawString(popup,.1,.9,PETSC_DRAW_BLUE,"Use left button to\n   enter cell.");
    106:   PetscDrawString(popup,.1,.7,PETSC_DRAW_BLUE,"Use center button to\n   end.");
    107:   PetscDrawFlush(popup);

    109:   PetscDrawGetMouseButton(draw,&button,&cx,&cy,PETSC_NULL,PETSC_NULL);
    110:   AOData2dGridAddNode(agrid,cx,cy,&cn);
    111:   cell[0] = cn;
    112:   sprintf(title,"Input grid: Number vertex %d Number cell %d",(int)agrid->vertex_n,(int)agrid->cell_n);
    113:   PetscDrawSetTitle(draw,title);
    114:   while (button == BUTTON_LEFT) {
    115:     /* wait for second vertex */
    116:     PetscDrawGetMouseButton(draw,&button,&cx,&cy,PETSC_NULL,PETSC_NULL);
    117:     if (button != BUTTON_LEFT) {
    118:       SETERRQ(PETSC_ERR_USER,"Must press left button to complete quadrilateral");
    119:     }
    120:     AOData2dGridAddNode(agrid,cx,cy,&cn);
    121:     cell[4*agrid->cell_n+1] = cn;
    122:     PetscDrawLine(draw,vertex[2*cell[4*agrid->cell_n]],vertex[1+2*cell[4*agrid->cell_n]],
    123:                          vertex[2*cell[4*agrid->cell_n+1]],vertex[1+2*cell[4*agrid->cell_n+1]],
    124:                          PETSC_DRAW_RED);
    125:     sprintf(title,"Input grid: Number vertex %d Number cell %d",(int)agrid->vertex_n,(int)agrid->cell_n);
    126:     PetscDrawSetTitle(draw,title);
    127:     /* wait for third vertex */
    128:     PetscDrawGetMouseButton(draw,&button,&cx,&cy,PETSC_NULL,PETSC_NULL);
    129:     if (button != BUTTON_LEFT) {
    130:       SETERRQ(PETSC_ERR_USER,"Must press left button to complete quadrilateral");
    131:     }
    132:     AOData2dGridAddNode(agrid,cx,cy,&cn);
    133:     cell[4*agrid->cell_n+2] = cn;
    134:     PetscDrawLine(draw,vertex[2*cell[4*agrid->cell_n+1]],vertex[1+2*cell[4*agrid->cell_n+1]],
    135:                          vertex[2*cell[4*agrid->cell_n+2]],vertex[1+2*cell[4*agrid->cell_n+2]],
    136:                          PETSC_DRAW_RED);
    137:     sprintf(title,"Input grid: Number vertex %d Number cell %d",(int)agrid->vertex_n,(int)agrid->cell_n);
    138:     PetscDrawSetTitle(draw,title);
    139:     /* wait for fourth vertex */
    140:     PetscDrawGetMouseButton(draw,&button,&cx,&cy,PETSC_NULL,PETSC_NULL);
    141:     if (button != BUTTON_LEFT) {
    142:       SETERRQ(PETSC_ERR_USER,"Must press left button to complete quadrilateral");
    143:     }
    144:     AOData2dGridAddNode(agrid,cx,cy,&cn);
    145:     cell[4*agrid->cell_n+3] = cn;
    146:     PetscDrawLine(draw,vertex[2*cell[4*agrid->cell_n+2]],vertex[1+2*cell[4*agrid->cell_n+2]],
    147:                          vertex[2*cell[4*agrid->cell_n+3]],vertex[1+2*cell[4*agrid->cell_n+3]],
    148:                          PETSC_DRAW_RED);
    149:     PetscDrawLine(draw,vertex[2*cell[4*agrid->cell_n]],vertex[1+2*cell[4*agrid->cell_n]],
    150:                          vertex[2*cell[4*agrid->cell_n+3]],vertex[1+2*cell[4*agrid->cell_n+3]],
    151:                          PETSC_DRAW_RED);
    152:     agrid->cell_n++;
    153:     sprintf(title,"Input grid: Number vertex %d Number cell %d",(int)agrid->vertex_n,(int)agrid->cell_n);
    154:     PetscDrawSetTitle(draw,title);

    156:     /* Get the first for the next cellralateral, or BUTTON_CENTER to end */
    157:     PetscDrawGetMouseButton(draw,&button,&cx,&cy,PETSC_NULL,PETSC_NULL);
    158:     if (button != BUTTON_LEFT) {break;}
    159:     AOData2dGridAddNode(agrid,cx,cy,&cn);
    160:     cell[4*agrid->cell_n] = cn;

    162:     sprintf(title,"Input grid: Number vertex %d Number cell %d",(int)agrid->vertex_n,(int)agrid->cell_n);
    163:     PetscDrawSetTitle(draw,title);
    164:   }
    165:   return(0);
    166: }

    170: /*
    171:    Changes the node numbering for the cell to make sure they are all in
    172:    clockwise ordering
    173: */
    174: PetscErrorCode PETSCDM_DLLEXPORT AOData2dGridFlipCells(AOData2dGrid agrid)
    175: {
    176:   PetscInt       i,*cell = agrid->cell_vertex, cell_n = agrid->cell_n;
    177:   PetscReal *vertex = agrid->vertex, sign;

    180:   for (i=0; i<cell_n; i++) {
    181:     /*
    182:        compute the quantity

    184:             x0      x1    x2      x3
    185:             y0      y1    y2      y3
    186:      */

    188:      sign = vertex[2*cell[4*i]]*vertex[1+2*cell[4*i+1]]   + vertex[2*cell[4*i+1]]*vertex[1+2*cell[4*i+2]] +
    189:             vertex[2*cell[4*i+2]]*vertex[1+2*cell[4*i+3]] + vertex[2*cell[4*i+3]]*vertex[1+2*cell[4*i]]   -
    190:             vertex[1+2*cell[4*i]]*vertex[2*cell[4*i+1]]   - vertex[1+2*cell[4*i+1]]*vertex[2*cell[4*i+2]] -
    191:             vertex[1+2*cell[4*i+2]]*vertex[2*cell[4*i+3]] - vertex[1+2*cell[4*i+3]]*vertex[2*cell[4*i]];

    193:      if (sign == 0.0) {
    194:        SETERRQ(PETSC_ERR_USER,"Bad cell, zero area");
    195:      } else if (sign > 0) {
    196:        PetscInt q1tmp = cell[4*i+1];
    197:        cell[4*i+1] = cell[4*i+3];
    198:        cell[4*i+3] = q1tmp;
    199:      }
    200:   }
    201:   return(0);
    202: }

    206: /*
    207:      AOData2dGridAddNode - Maintains a list of nodes given so far
    208: */
    209: PetscErrorCode PETSCDM_DLLEXPORT AOData2dGridAddNode(AOData2dGrid agrid, PetscReal cx, PetscReal cy, PetscInt *cn)
    210: {
    211:   PetscInt i;

    214:   for (i=0; i<agrid->vertex_n; i++) {
    215:     if ((PetscAbsReal(agrid->vertex[2*i] - cx) < 1.e-9) && (PetscAbsReal(agrid->vertex[1+2*i] - cy) < 1.e-9)) {
    216:       *cn = i;
    217:       return(0);
    218:     }
    219:   }
    220:   agrid->vertex[2*agrid->vertex_n] = cx;
    221:   agrid->vertex[1+2*agrid->vertex_n] = cy;
    222:   *cn     = (agrid->vertex_n)++;

    224:   if (cx < agrid->xmin)      agrid->xmin = cx;
    225:   else if (cx > agrid->xmax) agrid->xmax = cx;
    226:   if (cy < agrid->ymin)      agrid->ymin = cy;
    227:   else if (cy > agrid->ymax) agrid->ymax = cy;
    228:   return(0);
    229: }

    233: PetscErrorCode PETSCDM_DLLEXPORT AOData2dGridComputeNeighbors(AOData2dGrid agrid)
    234: {
    236:   PetscInt  i,j,*cell_edge,*edge_cell,*edge,*cell,*neighbors,e;

    239:   agrid->edge_max = 2*agrid->vertex_n;
    240:   agrid->edge_n   = 0;
    241:   PetscMalloc(2*agrid->edge_max*sizeof(PetscInt),&agrid->edge_vertex);
    242:   edge      = agrid->edge_vertex;
    243:   PetscMalloc(4*agrid->cell_max*sizeof(PetscInt),agrid->cell_edge);
    244:   cell_edge = agrid->cell_edge;
    245:   PetscMalloc(2*agrid->edge_max*sizeof(PetscInt),&agrid->edge_cell);
    246:   edge_cell = agrid->edge_cell;

    248:   cell      = agrid->cell_vertex;

    250:   /*
    251:        Mark all neighbors (to start) with -1 to indicate missing neighbor
    252:   */
    253:   for (i=0; i<2*agrid->edge_max; i++) {
    254:     edge_cell[i] = -1;
    255:   }

    257:   for (i=0; i<agrid->cell_n; i++) {
    258:     for (j=0; j<agrid->edge_n; j++) {
    259:       if (cell[4*i] == edge[2*j+1] && cell[4*i+1] == edge[2*j]) {
    260:         cell_edge[4*i]   = j;
    261:         edge_cell[2*j+1] = i;
    262:         goto found0;
    263:       }
    264:     }
    265:     /*
    266:        Add a new edge to the list
    267:     */
    268:     edge_cell[2*agrid->edge_n]   = i;
    269:     edge[2*agrid->edge_n]        = cell[4*i];
    270:     edge[2*agrid->edge_n+1]      = cell[4*i+1];
    271:     cell_edge[4*i]                = agrid->edge_n;
    272:     agrid->edge_n++;
    273:     found0:;
    274:     for (j=0; j<agrid->edge_n; j++) {
    275:       if (cell[4*i+1] == edge[2*j+1] && cell[4*i+2] == edge[2*j]) {
    276:         cell_edge[4*i+1] = j;
    277:         edge_cell[2*j+1] = i;
    278:         goto found1;
    279:       }
    280:     }
    281:     /*
    282:        Add a new edge to the list
    283:     */
    284:     edge_cell[2*agrid->edge_n]   = i;
    285:     edge[2*agrid->edge_n]        = cell[4*i+1];
    286:     edge[2*agrid->edge_n+1]      = cell[4*i+2];
    287:     cell_edge[4*i+1]              = agrid->edge_n;
    288:     agrid->edge_n++;
    289:     found1:;
    290:     for (j=0; j<agrid->edge_n; j++) {
    291:       if (cell[4*i+2] == edge[2*j+1] && cell[4*i+3] == edge[2*j]) {
    292:         cell_edge[4*i+2] = j;
    293:         edge_cell[2*j+1] = i;
    294:         goto found2;
    295:       }
    296:     }
    297:     /*
    298:        Add a new edge to the list
    299:     */
    300:     edge_cell[2*agrid->edge_n]   = i;
    301:     edge[2*agrid->edge_n]        = cell[4*i+2];
    302:     edge[2*agrid->edge_n+1]      = cell[4*i+3];
    303:     cell_edge[4*i+2]              = agrid->edge_n;
    304:     agrid->edge_n++;
    305:     found2:;
    306:     for (j=0; j<agrid->edge_n; j++) {
    307:       if (cell[4*i+3] == edge[2*j+1] && cell[4*i] == edge[2*j]) {
    308:         cell_edge[4*i+3] = j;
    309:         edge_cell[2*j+1] = i;
    310:         goto found3;
    311:       }
    312:     }
    313:     /*
    314:        Add a new edge to the list
    315:     */
    316:     edge_cell[2*agrid->edge_n]   = i;
    317:     edge[2*agrid->edge_n]        = cell[4*i+3];
    318:     edge[2*agrid->edge_n+1]      = cell[4*i];
    319:     cell_edge[4*i+3]              = agrid->edge_n;
    320:     agrid->edge_n++;
    321:     found3:;

    323:   }

    325:   PetscMalloc(4*agrid->cell_n*sizeof(PetscInt),&agrid->cell_cell);
    326:   neighbors = agrid->cell_cell;
    327:   for (i=0; i<agrid->cell_n; i++) {
    328:     for (j=0; j<4; j++) {
    329:       e = 2*agrid->cell_edge[4*i+j];

    331:       /* get the edge neighbor that is not the current cell */
    332:       if (i == agrid->edge_cell[e]) e++;
    333:       neighbors[4*i+j] = agrid->edge_cell[e];
    334:     }
    335:   }

    337:   return(0);
    338: }

    342: PetscErrorCode PETSCDM_DLLEXPORT AOData2dGridComputeVertexBoundary(AOData2dGrid agrid)
    343: {
    345:   PetscInt  i,j,*count,*cell_vertex = agrid->cell_vertex;

    348:   /*
    349:       allocate bitarray for boundary info
    350:   */
    351:   PetscBTCreate(agrid->vertex_n,agrid->vertex_boundary);

    353:   /*
    354:       count contains number of cell that contain the given vertex
    355:   */
    356:   PetscMalloc(agrid->vertex_n*sizeof(PetscInt),&count);
    357:   PetscMemzero(count,agrid->vertex_n*sizeof(PetscInt));

    359:   for (i=0; i<agrid->cell_n; i++) {
    360:     for (j=0; j<4; j++) {
    361:       count[cell_vertex[4*i+j]]++;
    362:     }
    363:   }

    366:   for (i=0; i<agrid->vertex_n; i++) {
    367:     /* UGLY! Just for a quick solution: I want Dirichlet b.c. only at left edge! */
    368:     PetscTruth neumann_bc;
    369:     PetscOptionsHasName(PETSC_NULL,"-dirichlet_on_left",&neumann_bc);
    370:     if (neumann_bc) {
    371:       if ((count[i] < 4) && (agrid->vertex[2*i] == agrid->xmin)) {
    372:         PetscBTSet(agrid->vertex_boundary,i);
    373:       }
    374:     } else {
    375:       if (count[i] < 4) {
    376:         PetscBTSet(agrid->vertex_boundary,i);
    377:       }
    378:     }
    379:   }

    381:   PetscFree(count);

    383:   return(0);
    384: }

    388: /*
    389:      Show the numbering of the vertex, cell and edge
    390: */
    391: PetscErrorCode PETSCDM_DLLEXPORT AOData2dGridDraw(AOData2dGrid agrid,PetscDraw draw)
    392: {
    394:   PetscInt       i, *cell = agrid->cell_vertex, *edge = agrid->edge_vertex;
    395:   char      str[5];
    396:   PetscReal *vertex = agrid->vertex,xx,yy,xmin,xmax,ymin,ymax,h,w;

    399:   w = agrid->xmax - agrid->xmin;
    400:   h = agrid->ymax - agrid->ymin;
    401:   xmin = agrid->xmin - .1*w;
    402:   xmax = agrid->xmax + .1*w;
    403:   ymin = agrid->ymin - .1*h;
    404:   ymax = agrid->ymax + .1*h;
    405:   PetscDrawSetCoordinates(draw,xmin,ymin,xmax,ymax);

    407:   /*
    408:      Number the vertex
    409:   */
    410:   for (i=0; i<agrid->vertex_n; i++) {
    411:     sprintf(str,"%d",(int)i);
    412:     PetscDrawString(draw,vertex[2*i],vertex[1+2*i],PETSC_DRAW_BLUE,str);
    413:   }

    415:   /*
    416:      Number the cell
    417:   */
    418:   for (i=0; i<agrid->cell_n; i++) {
    419:     sprintf(str,"%d",(int)i);
    420:     xx = .25*(vertex[2*cell[4*i]] + vertex[2*cell[4*i+1]] + vertex[2*cell[4*i+2]] + vertex[2*cell[4*i+3]]);
    421:     yy = .25*(vertex[1+2*cell[4*i]] + vertex[1+2*cell[4*i+1]] + vertex[1+2*cell[4*i+2]] + vertex[1+2*cell[4*i+3]]);
    422:     PetscDrawString(draw,xx,yy,PETSC_DRAW_GREEN,str);
    423:   }

    425:   /*
    426:      Number the edge
    427:   */
    428:   for (i=0; i<agrid->edge_n; i++) {
    429:     sprintf(str,"%d",(int)i);
    430:     xx = .5*(vertex[2*edge[2*i]] + vertex[2*edge[2*i+1]]);
    431:     yy = .5*(vertex[1+2*edge[2*i]] + vertex[1+2*edge[2*i+1]]);
    432:     PetscDrawLine(draw,vertex[2*edge[2*i]],vertex[1+2*edge[2*i]],vertex[2*edge[2*i+1]],vertex[1+2*edge[2*i+1]],PETSC_DRAW_BLACK);
    433:     PetscDrawString(draw,xx,yy,PETSC_DRAW_VIOLET,str);
    434:   }

    436:   return(0);
    437: }

    441: /*
    442:     Frees all the memory space allocated in AGrid
    443: */
    444: PetscErrorCode PETSCDM_DLLEXPORT AOData2dGridDestroy(AOData2dGrid agrid)
    445: {

    449:    PetscFree(agrid->vertex);
    450:    PetscFree(agrid->cell_vertex);
    451:    PetscFree(agrid->cell_edge);
    452:    PetscFree(agrid->edge_vertex);
    453:    PetscFree(agrid->edge_cell);
    454:    PetscFree(agrid->cell_cell);
    455:    PetscFree(agrid->vertex_boundary);
    456:    PetscFree(agrid);
    457:    return(0);
    458: }

    462: /*
    463:
    464: */
    465: PetscErrorCode PETSCDM_DLLEXPORT AOData2dGridCreate(AOData2dGrid *agrid)
    466: {
    469:   PetscNew(struct _p_AOData2dGrid,agrid);
    470:   return(0);
    471: }