Actual source code: rgpolygon.c

slepc-main 2024-11-09
Report Typos and Errors
  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 - Sets the vertices
 93: -  -rg_polygon_verticesi - Sets the 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: 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: RGPolygonSetVertices()
170: @*/
171: PetscErrorCode RGPolygonGetVertices(RG rg,PetscInt *n,PetscScalar **vr,PetscScalar **vi)
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: SLEPC_EXTERN PetscErrorCode RGCreate_Polygon(RG rg)
444: {
445:   RG_POLYGON     *polygon;

447:   PetscFunctionBegin;
448:   PetscCall(PetscNew(&polygon));
449:   rg->data = (void*)polygon;

451:   rg->ops->istrivial      = RGIsTrivial_Polygon;
452:   rg->ops->computecontour = RGComputeContour_Polygon;
453:   rg->ops->computebbox    = RGComputeBoundingBox_Polygon;
454:   rg->ops->checkinside    = RGCheckInside_Polygon;
455:   rg->ops->setfromoptions = RGSetFromOptions_Polygon;
456:   rg->ops->view           = RGView_Polygon;
457:   rg->ops->destroy        = RGDestroy_Polygon;
458:   PetscCall(PetscObjectComposeFunction((PetscObject)rg,"RGPolygonSetVertices_C",RGPolygonSetVertices_Polygon));
459:   PetscCall(PetscObjectComposeFunction((PetscObject)rg,"RGPolygonGetVertices_C",RGPolygonGetVertices_Polygon));
460:   PetscFunctionReturn(PETSC_SUCCESS);
461: }