Actual source code: rgpolygon.c

slepc-3.16.2 2022-02-01
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-2021, 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: 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: {
 52:   PetscInt       i;
 53:   RG_POLYGON     *ctx = (RG_POLYGON*)rg->data;

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

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

 84:    Logically Collective on rg

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

 92:    Options Database Keys:
 93: +  -rg_polygon_vertices - Sets the vertices
 94: -  -rg_polygon_verticesi - Sets the vertices (imaginary part)

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

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

107:    Level: advanced

109: .seealso: RGPolygonGetVertices()
110: @*/
111: PetscErrorCode RGPolygonSetVertices(RG rg,PetscInt n,PetscScalar vr[],PetscScalar vi[])
112: {

119: #if !defined(PETSC_USE_COMPLEX)
121: #endif
122:   PetscTryMethod(rg,"RGPolygonSetVertices_C",(RG,PetscInt,PetscScalar*,PetscScalar*),(rg,n,vr,vi));
123:   return(0);
124: }

126: static PetscErrorCode RGPolygonGetVertices_Polygon(RG rg,PetscInt *n,PetscScalar **vr,PetscScalar **vi)
127: {
129:   RG_POLYGON     *ctx = (RG_POLYGON*)rg->data;
130:   PetscInt       i;

133:   if (n) *n  = ctx->n;
134:   if (vr) {
135:     if (!ctx->n) *vr = NULL;
136:     else {
137:       PetscMalloc1(ctx->n,vr);
138:       for (i=0;i<ctx->n;i++) (*vr)[i] = ctx->vr[i];
139:     }
140:   }
141: #if !defined(PETSC_USE_COMPLEX)
142:   if (vi) {
143:     if (!ctx->n) *vi = NULL;
144:     else {
145:       PetscMalloc1(ctx->n,vi);
146:       for (i=0;i<ctx->n;i++) (*vi)[i] = ctx->vi[i];
147:     }
148:   }
149: #endif
150:   return(0);
151: }

153: /*@C
154:    RGPolygonGetVertices - Gets the vertices that define the polygon region.

156:    Not Collective

158:    Input Parameter:
159: .  rg     - the region context

161:    Output Parameters:
162: +  n  - number of vertices
163: .  vr - array of vertices
164: -  vi - array of vertices (imaginary part)

166:    Notes:
167:    The values passed by user with RGPolygonSetVertices() are returned (or null
168:    pointers otherwise).
169:    The returned arrays should be freed by the user when no longer needed.

171:    Level: advanced

173: .seealso: RGPolygonSetVertices()
174: @*/
175: PetscErrorCode RGPolygonGetVertices(RG rg,PetscInt *n,PetscScalar **vr,PetscScalar **vi)
176: {

181:   PetscUseMethod(rg,"RGPolygonGetVertices_C",(RG,PetscInt*,PetscScalar**,PetscScalar**),(rg,n,vr,vi));
182:   return(0);
183: }

185: PetscErrorCode RGView_Polygon(RG rg,PetscViewer viewer)
186: {
188:   RG_POLYGON     *ctx = (RG_POLYGON*)rg->data;
189:   PetscBool      isdraw,isascii;
190:   int            winw,winh;
191:   PetscDraw      draw;
192:   PetscDrawAxis  axis;
193:   PetscReal      a,b,c,d,ab,cd,lx,ly,w,x0,y0,x1,y1,scale=1.2;
194:   PetscInt       i;
195:   char           str[50];

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

263: PetscErrorCode RGIsTrivial_Polygon(RG rg,PetscBool *trivial)
264: {
265:   RG_POLYGON *ctx = (RG_POLYGON*)rg->data;

268:   *trivial = PetscNot(ctx->n);
269:   return(0);
270: }

272: PetscErrorCode RGComputeContour_Polygon(RG rg,PetscInt n,PetscScalar *ucr,PetscScalar *uci)
273: {
275:   RG_POLYGON     *ctx = (RG_POLYGON*)rg->data;
276:   PetscReal      length,h,d,rem=0.0;
277:   PetscInt       k=1,idx=ctx->n-1,i;
278:   PetscBool      ini=PETSC_FALSE;
279:   PetscScalar    incr,*cr=ucr,*ci=uci;
280: #if !defined(PETSC_USE_COMPLEX)
281:   PetscScalar    inci;
282: #endif

285:   if (!ctx->n) SETERRQ(PetscObjectComm((PetscObject)rg),PETSC_ERR_ARG_WRONGSTATE,"No vertices have been set yet");
286:   length = SlepcAbsEigenvalue(ctx->vr[0]-ctx->vr[ctx->n-1],ctx->vi[0]-ctx->vi[ctx->n-1]);
287:   for (i=0;i<ctx->n-1;i++) length += SlepcAbsEigenvalue(ctx->vr[i]-ctx->vr[i+1],ctx->vi[i]-ctx->vi[i+1]);
288:   h = length/n;
289:   if (!ucr) { PetscMalloc1(n,&cr); }
290:   if (!uci) { PetscMalloc1(n,&ci); }
291:   cr[0] = ctx->vr[0];
292: #if !defined(PETSC_USE_COMPLEX)
293:   ci[0] = ctx->vi[0];
294: #endif
295:   incr = ctx->vr[ctx->n-1]-ctx->vr[0];
296: #if !defined(PETSC_USE_COMPLEX)
297:   inci = ctx->vi[ctx->n-1]-ctx->vi[0];
298: #endif
299:   d = SlepcAbsEigenvalue(incr,inci);
300:   incr /= d;
301: #if !defined(PETSC_USE_COMPLEX)
302:   inci /= d;
303: #endif
304:   while (k<n) {
305:     if (ini) {
306:       incr = ctx->vr[idx]-ctx->vr[idx+1];
307: #if !defined(PETSC_USE_COMPLEX)
308:       inci = ctx->vi[idx]-ctx->vi[idx+1];
309: #endif
310:       d = SlepcAbsEigenvalue(incr,inci);
311:       incr /= d;
312: #if !defined(PETSC_USE_COMPLEX)
313:       inci /= d;
314: #endif
315:       if (rem+d>h) {
316:         cr[k] = ctx->vr[idx+1]+incr*(h-rem);
317: #if !defined(PETSC_USE_COMPLEX)
318:         ci[k] = ctx->vi[idx+1]+inci*(h-rem);
319: #endif
320:         k++;
321:         ini = PETSC_FALSE;
322:       } else {rem += d; idx--;}
323:     } else {
324: #if !defined(PETSC_USE_COMPLEX)
325:       rem = SlepcAbsEigenvalue(ctx->vr[idx]-cr[k-1],ctx->vi[idx]-ci[k-1]);
326: #else
327:       rem = PetscAbsScalar(ctx->vr[idx]-cr[k-1]);
328: #endif
329:       if (rem>h) {
330:         cr[k] = cr[k-1]+incr*h;
331: #if !defined(PETSC_USE_COMPLEX)
332:         ci[k] = ci[k-1]+inci*h;
333: #endif
334:         k++;
335:       } else {ini = PETSC_TRUE; idx--;}
336:     }
337:   }
338:   if (!ucr) { PetscFree(cr); }
339:   if (!uci) { PetscFree(ci); }
340:   return(0);
341: }

343: PetscErrorCode RGComputeBoundingBox_Polygon(RG rg,PetscReal *a,PetscReal *b,PetscReal *c,PetscReal *d)
344: {
345:   RG_POLYGON *ctx = (RG_POLYGON*)rg->data;
346:   PetscInt   i;

349:   if (a) *a =  PETSC_MAX_REAL;
350:   if (b) *b = -PETSC_MAX_REAL;
351:   if (c) *c =  PETSC_MAX_REAL;
352:   if (d) *d = -PETSC_MAX_REAL;
353:   for (i=0;i<ctx->n;i++) {
354: #if defined(PETSC_USE_COMPLEX)
355:     if (a) *a = PetscMin(*a,PetscRealPart(ctx->vr[i]));
356:     if (b) *b = PetscMax(*b,PetscRealPart(ctx->vr[i]));
357:     if (c) *c = PetscMin(*c,PetscImaginaryPart(ctx->vr[i]));
358:     if (d) *d = PetscMax(*d,PetscImaginaryPart(ctx->vr[i]));
359: #else
360:     if (a) *a = PetscMin(*a,ctx->vr[i]);
361:     if (b) *b = PetscMax(*b,ctx->vr[i]);
362:     if (c) *c = PetscMin(*c,ctx->vi[i]);
363:     if (d) *d = PetscMax(*d,ctx->vi[i]);
364: #endif
365:   }
366:   return(0);
367: }

369: PetscErrorCode RGCheckInside_Polygon(RG rg,PetscReal px,PetscReal py,PetscInt *inout)
370: {
371:   RG_POLYGON *ctx = (RG_POLYGON*)rg->data;
372:   PetscReal  val,x[VERTMAX],y[VERTMAX];
373:   PetscBool  mx,my,nx,ny;
374:   PetscInt   i,j;

377:   for (i=0;i<ctx->n;i++) {
378: #if defined(PETSC_USE_COMPLEX)
379:     x[i] = PetscRealPart(ctx->vr[i])-px;
380:     y[i] = PetscImaginaryPart(ctx->vr[i])-py;
381: #else
382:     x[i] = ctx->vr[i]-px;
383:     y[i] = ctx->vi[i]-py;
384: #endif
385:   }
386:   *inout = -1;
387:   for (i=0;i<ctx->n;i++) {
388:     j = (i+1)%ctx->n;
389:     mx = PetscNot(x[i]<0.0);
390:     nx = PetscNot(x[j]<0.0);
391:     my = PetscNot(y[i]<0.0);
392:     ny = PetscNot(y[j]<0.0);
393:     if (!((my||ny) && (mx||nx)) || (mx&&nx)) continue;
394:     if (((my && ny && (mx||nx)) && (!(mx&&nx)))) {
395:       *inout = -*inout;
396:       continue;
397:     }
398:     val = (y[i]*x[j]-x[i]*y[j])/(x[j]-x[i]);
399:     if (PetscAbs(val)<10*PETSC_MACHINE_EPSILON) {
400:       *inout = 0;
401:       return(0);
402:     } else if (val>0.0) *inout = -*inout;
403:   }
404:   return(0);
405: }

407: PetscErrorCode RGSetFromOptions_Polygon(PetscOptionItems *PetscOptionsObject,RG rg)
408: {
410:   PetscScalar    array[VERTMAX];
411:   PetscInt       i,k;
412:   PetscBool      flg,flgi=PETSC_FALSE;
413: #if !defined(PETSC_USE_COMPLEX)
414:   PetscScalar    arrayi[VERTMAX];
415:   PetscInt       ki;
416: #else
417:   PetscScalar    *arrayi=NULL;
418: #endif

421:   PetscOptionsHead(PetscOptionsObject,"RG Polygon Options");

423:     k = VERTMAX;
424:     for (i=0;i<k;i++) array[i] = 0;
425:     PetscOptionsScalarArray("-rg_polygon_vertices","Vertices of polygon","RGPolygonSetVertices",array,&k,&flg);
426: #if !defined(PETSC_USE_COMPLEX)
427:     ki = VERTMAX;
428:     for (i=0;i<ki;i++) arrayi[i] = 0;
429:     PetscOptionsScalarArray("-rg_polygon_verticesi","Vertices of polygon (imaginary part)","RGPolygonSetVertices",arrayi,&ki,&flgi);
430:     if (ki!=k) SETERRQ2(PetscObjectComm((PetscObject)rg),PETSC_ERR_ARG_SIZ,"The number of real %D and imaginary %D parts do not match",k,ki);
431: #endif
432:     if (flg || flgi) { RGPolygonSetVertices(rg,k,array,arrayi); }

434:   PetscOptionsTail();
435:   return(0);
436: }

438: PetscErrorCode RGDestroy_Polygon(RG rg)
439: {
441:   RG_POLYGON     *ctx = (RG_POLYGON*)rg->data;

444:   if (ctx->n) {
445:     PetscFree(ctx->vr);
446: #if !defined(PETSC_USE_COMPLEX)
447:     PetscFree(ctx->vi);
448: #endif
449:   }
450:   PetscFree(rg->data);
451:   PetscObjectComposeFunction((PetscObject)rg,"RGPolygonSetVertices_C",NULL);
452:   PetscObjectComposeFunction((PetscObject)rg,"RGPolygonGetVertices_C",NULL);
453:   return(0);
454: }

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

462:   PetscNewLog(rg,&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:   PetscObjectComposeFunction((PetscObject)rg,"RGPolygonSetVertices_C",RGPolygonSetVertices_Polygon);
473:   PetscObjectComposeFunction((PetscObject)rg,"RGPolygonGetVertices_C",RGPolygonGetVertices_Polygon);
474:   return(0);
475: }