Actual source code: rgpolygon.c
slepc-3.16.2 2022-02-01
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: }