Actual source code: rgpolygon.c

      1: /*
      2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
      4:    Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain

      6:    This file is part of SLEPc.
      7:    SLEPc is distributed under a 2-clause BSD license (see LICENSE).
      8:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      9: */
     10: /*
     11:    Polygonal region defined by a set of vertices
     12: */

     14: #include <slepc/private/rgimpl.h>
     15: #include <petscdraw.h>

     17: #define VERTMAX 30

     19: typedef struct {
     20:   PetscInt    n;         /* number of vertices */
     21:   PetscScalar *vr,*vi;   /* array of vertices (vi not used in complex scalars) */
     22: } RG_POLYGON;

     24: static PetscErrorCode RGComputeBoundingBox_Polygon(RG,PetscReal*,PetscReal*,PetscReal*,PetscReal*);

     26: #if !defined(PETSC_USE_COMPLEX)
     27: static PetscBool CheckSymmetry(PetscInt n,PetscScalar *vr,PetscScalar *vi)
     28: {
     29:   PetscInt i,j,k;
     30:   /* find change of sign in imaginary part */
     31:   j = vi[0]!=0.0? 0: 1;
     32:   for (k=j+1;k<n;k++) {
     33:     if (vi[k]!=0.0) {
     34:       if (vi[k]*vi[j]<0.0) break;
     35:       j++;
     36:     }
     37:   }
     38:   if (k==n) return (j==1)? PETSC_TRUE: PETSC_FALSE;
     39:   /* check pairing vertices */
     40:   for (i=0;i<n/2;i++) {
     41:     if (vr[k]!=vr[j] || vi[k]!=-vi[j]) return PETSC_FALSE;
     42:     k = (k+1)%n;
     43:     j = (j-1+n)%n;
     44:   }
     45:   return PETSC_TRUE;
     46: }
     47: #endif

     49: static PetscErrorCode RGPolygonSetVertices_Polygon(RG rg,PetscInt n,PetscScalar *vr,PetscScalar *vi)
     50: {
     51:   PetscInt       i;
     52:   RG_POLYGON     *ctx = (RG_POLYGON*)rg->data;

     54:   PetscFunctionBegin;
     55:   PetscCheck(n>2,PetscObjectComm((PetscObject)rg),PETSC_ERR_ARG_OUTOFRANGE,"At least 3 vertices required, you provided %" PetscInt_FMT,n);
     56:   PetscCheck(n<=VERTMAX,PetscObjectComm((PetscObject)rg),PETSC_ERR_ARG_OUTOFRANGE,"Too many points, maximum allowed is %d",VERTMAX);
     57: #if !defined(PETSC_USE_COMPLEX)
     58:   PetscCheck(CheckSymmetry(n,vr,vi),PetscObjectComm((PetscObject)rg),PETSC_ERR_ARG_WRONG,"In real scalars the region must be symmetric wrt real axis");
     59: #endif
     60:   if (ctx->n) {
     61:     PetscCall(PetscFree(ctx->vr));
     62: #if !defined(PETSC_USE_COMPLEX)
     63:     PetscCall(PetscFree(ctx->vi));
     64: #endif
     65:   }
     66:   ctx->n = n;
     67:   PetscCall(PetscMalloc1(n,&ctx->vr));
     68: #if !defined(PETSC_USE_COMPLEX)
     69:   PetscCall(PetscMalloc1(n,&ctx->vi));
     70: #endif
     71:   for (i=0;i<n;i++) {
     72:     ctx->vr[i] = vr[i];
     73: #if !defined(PETSC_USE_COMPLEX)
     74:     ctx->vi[i] = vi[i];
     75: #endif
     76:   }
     77:   PetscFunctionReturn(PETSC_SUCCESS);
     78: }

     80: /*@
     81:    RGPolygonSetVertices - Sets the vertices that define the polygon region.

     83:    Logically Collective

     85:    Input Parameters:
     86: +  rg - the region context
     87: .  n  - number of vertices
     88: .  vr - array of vertices
     89: -  vi - array of vertices (imaginary part)

     91:    Options Database Keys:
     92: +  -rg_polygon_vertices  vr1,vr2,... - comma-separated list of vertices
     93: -  -rg_polygon_verticesi vi1,vi2,... - comma-separated list of vertices (imaginary part)

     95:    Notes:
     96:    In the case of complex scalars, only argument `vr` is used, containing
     97:    the complex vertices; the list of vertices can be provided in the
     98:    command line with a comma-separated list of complex values
     99:    `[+/-][realnumber][+/-]realnumberi` with no spaces.

    101:    When PETSc is built with real scalars, the real and imaginary parts of
    102:    the vertices must be provided in two separate arrays (or two lists in
    103:    the command line). In this case, the region must be symmetric with
    104:    respect to the real axis.

    106:    Level: advanced

    108: .seealso: [](sec:rg), `RGPOLYGON`, `RGPolygonGetVertices()`
    109: @*/
    110: PetscErrorCode RGPolygonSetVertices(RG rg,PetscInt n,PetscScalar vr[],PetscScalar vi[])
    111: {
    112:   PetscFunctionBegin;
    115:   PetscAssertPointer(vr,3);
    116: #if !defined(PETSC_USE_COMPLEX)
    117:   PetscAssertPointer(vi,4);
    118: #endif
    119:   PetscTryMethod(rg,"RGPolygonSetVertices_C",(RG,PetscInt,PetscScalar*,PetscScalar*),(rg,n,vr,vi));
    120:   PetscFunctionReturn(PETSC_SUCCESS);
    121: }

    123: static PetscErrorCode RGPolygonGetVertices_Polygon(RG rg,PetscInt *n,PetscScalar *vr[],PetscScalar *vi[])
    124: {
    125:   RG_POLYGON     *ctx = (RG_POLYGON*)rg->data;
    126:   PetscInt       i;

    128:   PetscFunctionBegin;
    129:   if (n) *n  = ctx->n;
    130:   if (vr) {
    131:     if (!ctx->n) *vr = NULL;
    132:     else {
    133:       PetscCall(PetscMalloc1(ctx->n,vr));
    134:       for (i=0;i<ctx->n;i++) (*vr)[i] = ctx->vr[i];
    135:     }
    136:   }
    137: #if !defined(PETSC_USE_COMPLEX)
    138:   if (vi) {
    139:     if (!ctx->n) *vi = NULL;
    140:     else {
    141:       PetscCall(PetscMalloc1(ctx->n,vi));
    142:       for (i=0;i<ctx->n;i++) (*vi)[i] = ctx->vi[i];
    143:     }
    144:   }
    145: #endif
    146:   PetscFunctionReturn(PETSC_SUCCESS);
    147: }

    149: /*@C
    150:    RGPolygonGetVertices - Gets the vertices that define the polygon region.

    152:    Not Collective

    154:    Input Parameter:
    155: .  rg     - the region context

    157:    Output Parameters:
    158: +  n  - number of vertices
    159: .  vr - array of vertices
    160: -  vi - array of vertices (imaginary part)

    162:    Notes:
    163:    The values passed by user with `RGPolygonSetVertices()` are returned (or `NULL`
    164:    pointers otherwise).
    165:    The returned arrays should be freed by the user when no longer needed.

    167:    Level: advanced

    169: .seealso: [](sec:rg), `RGPOLYGON`, `RGPolygonSetVertices()`
    170: @*/
    171: PetscErrorCode RGPolygonGetVertices(RG rg,PetscInt *n,PetscScalar *vr[],PetscScalar *vi[]) PeNS
    172: {
    173:   PetscFunctionBegin;
    175:   PetscUseMethod(rg,"RGPolygonGetVertices_C",(RG,PetscInt*,PetscScalar**,PetscScalar**),(rg,n,vr,vi));
    176:   PetscFunctionReturn(PETSC_SUCCESS);
    177: }

    179: static PetscErrorCode RGView_Polygon(RG rg,PetscViewer viewer)
    180: {
    181:   RG_POLYGON     *ctx = (RG_POLYGON*)rg->data;
    182:   PetscBool      isdraw,isascii;
    183:   int            winw,winh;
    184:   PetscDraw      draw;
    185:   PetscDrawAxis  axis;
    186:   PetscReal      a,b,c,d,ab,cd,lx,ly,w,x0,y0,x1,y1,scale=1.2;
    187:   PetscInt       i;
    188:   char           str[50];

    190:   PetscFunctionBegin;
    191:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw));
    192:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
    193:   if (isascii) {
    194:     PetscCall(PetscViewerASCIIPrintf(viewer,"  vertices: "));
    195:     PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_FALSE));
    196:     for (i=0;i<ctx->n;i++) {
    197: #if defined(PETSC_USE_COMPLEX)
    198:       PetscCall(SlepcSNPrintfScalar(str,sizeof(str),ctx->vr[i],PETSC_FALSE));
    199: #else
    200:       if (ctx->vi[i]!=0.0) PetscCall(PetscSNPrintf(str,sizeof(str),"%g%+gi",(double)ctx->vr[i],(double)ctx->vi[i]));
    201:       else PetscCall(PetscSNPrintf(str,sizeof(str),"%g",(double)ctx->vr[i]));
    202: #endif
    203:       PetscCall(PetscViewerASCIIPrintf(viewer,"%s%s",str,(i<ctx->n-1)?", ":""));
    204:     }
    205:     PetscCall(PetscViewerASCIIPrintf(viewer,"\n"));
    206:     PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_TRUE));
    207:   } else if (isdraw) {
    208:     PetscCall(PetscViewerDrawGetDraw(viewer,0,&draw));
    209:     PetscCall(PetscDrawCheckResizedWindow(draw));
    210:     PetscCall(PetscDrawGetWindowSize(draw,&winw,&winh));
    211:     winw = PetscMax(winw,1); winh = PetscMax(winh,1);
    212:     PetscCall(PetscDrawClear(draw));
    213:     PetscCall(PetscDrawSetTitle(draw,"Polygonal region"));
    214:     PetscCall(PetscDrawAxisCreate(draw,&axis));
    215:     PetscCall(RGComputeBoundingBox_Polygon(rg,&a,&b,&c,&d));
    216:     a *= rg->sfactor;
    217:     b *= rg->sfactor;
    218:     c *= rg->sfactor;
    219:     d *= rg->sfactor;
    220:     lx = b-a;
    221:     ly = d-c;
    222:     ab = (a+b)/2;
    223:     cd = (c+d)/2;
    224:     w  = scale*PetscMax(lx/winw,ly/winh)/2;
    225:     PetscCall(PetscDrawAxisSetLimits(axis,ab-w*winw,ab+w*winw,cd-w*winh,cd+w*winh));
    226:     PetscCall(PetscDrawAxisDraw(axis));
    227:     PetscCall(PetscDrawAxisDestroy(&axis));
    228:     for (i=0;i<ctx->n;i++) {
    229: #if defined(PETSC_USE_COMPLEX)
    230:       x0 = PetscRealPart(ctx->vr[i]); y0 = PetscImaginaryPart(ctx->vr[i]);
    231:       if (i<ctx->n-1) {
    232:         x1 = PetscRealPart(ctx->vr[i+1]); y1 = PetscImaginaryPart(ctx->vr[i+1]);
    233:       } else {
    234:         x1 = PetscRealPart(ctx->vr[0]); y1 = PetscImaginaryPart(ctx->vr[0]);
    235:       }
    236: #else
    237:       x0 = ctx->vr[i]; y0 = ctx->vi[i];
    238:       if (i<ctx->n-1) {
    239:         x1 = ctx->vr[i+1]; y1 = ctx->vi[i+1];
    240:       } else {
    241:         x1 = ctx->vr[0]; y1 = ctx->vi[0];
    242:       }
    243: #endif
    244:       PetscCall(PetscDrawLine(draw,x0*rg->sfactor,y0*rg->sfactor,x1*rg->sfactor,y1*rg->sfactor,PETSC_DRAW_MAGENTA));
    245:     }
    246:     PetscCall(PetscDrawFlush(draw));
    247:     PetscCall(PetscDrawSave(draw));
    248:     PetscCall(PetscDrawPause(draw));
    249:   }
    250:   PetscFunctionReturn(PETSC_SUCCESS);
    251: }

    253: static PetscErrorCode RGIsTrivial_Polygon(RG rg,PetscBool *trivial)
    254: {
    255:   RG_POLYGON *ctx = (RG_POLYGON*)rg->data;

    257:   PetscFunctionBegin;
    258:   *trivial = PetscNot(ctx->n);
    259:   PetscFunctionReturn(PETSC_SUCCESS);
    260: }

    262: static PetscErrorCode RGComputeContour_Polygon(RG rg,PetscInt n,PetscScalar *ucr,PetscScalar *uci)
    263: {
    264:   RG_POLYGON     *ctx = (RG_POLYGON*)rg->data;
    265:   PetscReal      length,h,d,rem=0.0;
    266:   PetscInt       k=1,idx=ctx->n-1,i;
    267:   PetscBool      ini=PETSC_FALSE;
    268:   PetscScalar    incr,*cr=ucr,*ci=uci;
    269: #if !defined(PETSC_USE_COMPLEX)
    270:   PetscScalar    inci;
    271: #endif

    273:   PetscFunctionBegin;
    274:   PetscCheck(ctx->n,PetscObjectComm((PetscObject)rg),PETSC_ERR_ARG_WRONGSTATE,"No vertices have been set yet");
    275:   length = SlepcAbsEigenvalue(ctx->vr[0]-ctx->vr[ctx->n-1],ctx->vi[0]-ctx->vi[ctx->n-1]);
    276:   for (i=0;i<ctx->n-1;i++) length += SlepcAbsEigenvalue(ctx->vr[i]-ctx->vr[i+1],ctx->vi[i]-ctx->vi[i+1]);
    277:   h = length/n;
    278:   if (!ucr) PetscCall(PetscMalloc1(n,&cr));
    279:   if (!uci) PetscCall(PetscMalloc1(n,&ci));
    280:   cr[0] = ctx->vr[0];
    281: #if !defined(PETSC_USE_COMPLEX)
    282:   ci[0] = ctx->vi[0];
    283: #endif
    284:   incr = ctx->vr[ctx->n-1]-ctx->vr[0];
    285: #if !defined(PETSC_USE_COMPLEX)
    286:   inci = ctx->vi[ctx->n-1]-ctx->vi[0];
    287: #endif
    288:   d = SlepcAbsEigenvalue(incr,inci);
    289:   incr /= d;
    290: #if !defined(PETSC_USE_COMPLEX)
    291:   inci /= d;
    292: #endif
    293:   while (k<n) {
    294:     if (ini) {
    295:       incr = ctx->vr[idx]-ctx->vr[idx+1];
    296: #if !defined(PETSC_USE_COMPLEX)
    297:       inci = ctx->vi[idx]-ctx->vi[idx+1];
    298: #endif
    299:       d = SlepcAbsEigenvalue(incr,inci);
    300:       incr /= d;
    301: #if !defined(PETSC_USE_COMPLEX)
    302:       inci /= d;
    303: #endif
    304:       if (rem+d>h) {
    305:         cr[k] = ctx->vr[idx+1]+incr*(h-rem);
    306: #if !defined(PETSC_USE_COMPLEX)
    307:         ci[k] = ctx->vi[idx+1]+inci*(h-rem);
    308: #endif
    309:         k++;
    310:         ini = PETSC_FALSE;
    311:       } else {rem += d; idx--;}
    312:     } else {
    313: #if !defined(PETSC_USE_COMPLEX)
    314:       rem = SlepcAbsEigenvalue(ctx->vr[idx]-cr[k-1],ctx->vi[idx]-ci[k-1]);
    315: #else
    316:       rem = PetscAbsScalar(ctx->vr[idx]-cr[k-1]);
    317: #endif
    318:       if (rem>h) {
    319:         cr[k] = cr[k-1]+incr*h;
    320: #if !defined(PETSC_USE_COMPLEX)
    321:         ci[k] = ci[k-1]+inci*h;
    322: #endif
    323:         k++;
    324:       } else {ini = PETSC_TRUE; idx--;}
    325:     }
    326:   }
    327:   if (!ucr) PetscCall(PetscFree(cr));
    328:   if (!uci) PetscCall(PetscFree(ci));
    329:   PetscFunctionReturn(PETSC_SUCCESS);
    330: }

    332: static PetscErrorCode RGComputeBoundingBox_Polygon(RG rg,PetscReal *a,PetscReal *b,PetscReal *c,PetscReal *d)
    333: {
    334:   RG_POLYGON *ctx = (RG_POLYGON*)rg->data;
    335:   PetscInt   i;

    337:   PetscFunctionBegin;
    338:   if (a) *a =  PETSC_MAX_REAL;
    339:   if (b) *b = -PETSC_MAX_REAL;
    340:   if (c) *c =  PETSC_MAX_REAL;
    341:   if (d) *d = -PETSC_MAX_REAL;
    342:   for (i=0;i<ctx->n;i++) {
    343: #if defined(PETSC_USE_COMPLEX)
    344:     if (a) *a = PetscMin(*a,PetscRealPart(ctx->vr[i]));
    345:     if (b) *b = PetscMax(*b,PetscRealPart(ctx->vr[i]));
    346:     if (c) *c = PetscMin(*c,PetscImaginaryPart(ctx->vr[i]));
    347:     if (d) *d = PetscMax(*d,PetscImaginaryPart(ctx->vr[i]));
    348: #else
    349:     if (a) *a = PetscMin(*a,ctx->vr[i]);
    350:     if (b) *b = PetscMax(*b,ctx->vr[i]);
    351:     if (c) *c = PetscMin(*c,ctx->vi[i]);
    352:     if (d) *d = PetscMax(*d,ctx->vi[i]);
    353: #endif
    354:   }
    355:   PetscFunctionReturn(PETSC_SUCCESS);
    356: }

    358: static PetscErrorCode RGCheckInside_Polygon(RG rg,PetscReal px,PetscReal py,PetscInt *inout)
    359: {
    360:   RG_POLYGON *ctx = (RG_POLYGON*)rg->data;
    361:   PetscReal  val,x[VERTMAX],y[VERTMAX];
    362:   PetscBool  mx,my,nx,ny;
    363:   PetscInt   i,j;

    365:   PetscFunctionBegin;
    366:   for (i=0;i<ctx->n;i++) {
    367: #if defined(PETSC_USE_COMPLEX)
    368:     x[i] = PetscRealPart(ctx->vr[i])-px;
    369:     y[i] = PetscImaginaryPart(ctx->vr[i])-py;
    370: #else
    371:     x[i] = ctx->vr[i]-px;
    372:     y[i] = ctx->vi[i]-py;
    373: #endif
    374:   }
    375:   *inout = -1;
    376:   for (i=0;i<ctx->n;i++) {
    377:     j = (i+1)%ctx->n;
    378:     mx = PetscNot(x[i]<0.0);
    379:     nx = PetscNot(x[j]<0.0);
    380:     my = PetscNot(y[i]<0.0);
    381:     ny = PetscNot(y[j]<0.0);
    382:     if (!((my||ny) && (mx||nx)) || (mx&&nx)) continue;
    383:     if (((my && ny && (mx||nx)) && (!(mx&&nx)))) {
    384:       *inout = -*inout;
    385:       continue;
    386:     }
    387:     val = (y[i]*x[j]-x[i]*y[j])/(x[j]-x[i]);
    388:     if (PetscAbs(val)<10*PETSC_MACHINE_EPSILON) {
    389:       *inout = 0;
    390:       PetscFunctionReturn(PETSC_SUCCESS);
    391:     } else if (val>0.0) *inout = -*inout;
    392:   }
    393:   PetscFunctionReturn(PETSC_SUCCESS);
    394: }

    396: static PetscErrorCode RGSetFromOptions_Polygon(RG rg,PetscOptionItems PetscOptionsObject)
    397: {
    398:   PetscScalar    array[VERTMAX];
    399:   PetscInt       i,k;
    400:   PetscBool      flg,flgi=PETSC_FALSE;
    401: #if !defined(PETSC_USE_COMPLEX)
    402:   PetscScalar    arrayi[VERTMAX];
    403:   PetscInt       ki;
    404: #else
    405:   PetscScalar    *arrayi=NULL;
    406: #endif

    408:   PetscFunctionBegin;
    409:   PetscOptionsHeadBegin(PetscOptionsObject,"RG Polygon Options");

    411:     k = VERTMAX;
    412:     for (i=0;i<k;i++) array[i] = 0;
    413:     PetscCall(PetscOptionsScalarArray("-rg_polygon_vertices","Vertices of polygon","RGPolygonSetVertices",array,&k,&flg));
    414: #if !defined(PETSC_USE_COMPLEX)
    415:     ki = VERTMAX;
    416:     for (i=0;i<ki;i++) arrayi[i] = 0;
    417:     PetscCall(PetscOptionsScalarArray("-rg_polygon_verticesi","Vertices of polygon (imaginary part)","RGPolygonSetVertices",arrayi,&ki,&flgi));
    418:     PetscCheck(ki==k,PetscObjectComm((PetscObject)rg),PETSC_ERR_ARG_SIZ,"The number of real %" PetscInt_FMT " and imaginary %" PetscInt_FMT " parts do not match",k,ki);
    419: #endif
    420:     if (flg || flgi) PetscCall(RGPolygonSetVertices(rg,k,array,arrayi));

    422:   PetscOptionsHeadEnd();
    423:   PetscFunctionReturn(PETSC_SUCCESS);
    424: }

    426: static PetscErrorCode RGDestroy_Polygon(RG rg)
    427: {
    428:   RG_POLYGON     *ctx = (RG_POLYGON*)rg->data;

    430:   PetscFunctionBegin;
    431:   if (ctx->n) {
    432:     PetscCall(PetscFree(ctx->vr));
    433: #if !defined(PETSC_USE_COMPLEX)
    434:     PetscCall(PetscFree(ctx->vi));
    435: #endif
    436:   }
    437:   PetscCall(PetscFree(rg->data));
    438:   PetscCall(PetscObjectComposeFunction((PetscObject)rg,"RGPolygonSetVertices_C",NULL));
    439:   PetscCall(PetscObjectComposeFunction((PetscObject)rg,"RGPolygonGetVertices_C",NULL));
    440:   PetscFunctionReturn(PETSC_SUCCESS);
    441: }

    443: /*MC
    444:    RGPOLYGON - RGPOLYGON = "polygon" - A polygonal region defined by its
    445:    vertices, given via `RGPolygonSetVertices()`.

    447:    Note:
    448:    The following figure shows an example of a polygonal region.

    450:    ![Polygonal region](../../_static/images/manual/svg/fig-rg-polygon.svg)

    452:    Level: beginner

    454: .seealso: [](sec:rg), `RG`, `RGType`, `RGSetType()`, `RGPolygonSetVertices()`
    455: M*/

    457: SLEPC_EXTERN PetscErrorCode RGCreate_Polygon(RG rg)
    458: {
    459:   RG_POLYGON     *polygon;

    461:   PetscFunctionBegin;
    462:   PetscCall(PetscNew(&polygon));
    463:   rg->data = (void*)polygon;

    465:   rg->ops->istrivial      = RGIsTrivial_Polygon;
    466:   rg->ops->computecontour = RGComputeContour_Polygon;
    467:   rg->ops->computebbox    = RGComputeBoundingBox_Polygon;
    468:   rg->ops->checkinside    = RGCheckInside_Polygon;
    469:   rg->ops->setfromoptions = RGSetFromOptions_Polygon;
    470:   rg->ops->view           = RGView_Polygon;
    471:   rg->ops->destroy        = RGDestroy_Polygon;
    472:   PetscCall(PetscObjectComposeFunction((PetscObject)rg,"RGPolygonSetVertices_C",RGPolygonSetVertices_Polygon));
    473:   PetscCall(PetscObjectComposeFunction((PetscObject)rg,"RGPolygonGetVertices_C",RGPolygonGetVertices_Polygon));
    474:   PetscFunctionReturn(PETSC_SUCCESS);
    475: }