Actual source code: rgpolygon.c
slepc-3.18.0 2022-10-01
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: 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;
56: #if !defined(PETSC_USE_COMPLEX)
58: #endif
59: if (ctx->n) {
60: PetscFree(ctx->vr);
61: #if !defined(PETSC_USE_COMPLEX)
62: PetscFree(ctx->vi);
63: #endif
64: }
65: ctx->n = n;
66: PetscMalloc1(n,&ctx->vr);
67: #if !defined(PETSC_USE_COMPLEX)
68: PetscMalloc1(n,&ctx->vi);
69: #endif
70: for (i=0;i<n;i++) {
71: ctx->vr[i] = vr[i];
72: #if !defined(PETSC_USE_COMPLEX)
73: ctx->vi[i] = vi[i];
74: #endif
75: }
76: return 0;
77: }
79: /*@
80: RGPolygonSetVertices - Sets the vertices that define the polygon region.
82: Logically Collective on rg
84: Input Parameters:
85: + rg - the region context
86: . n - number of vertices
87: . vr - array of vertices
88: - vi - array of vertices (imaginary part)
90: Options Database Keys:
91: + -rg_polygon_vertices - Sets the vertices
92: - -rg_polygon_verticesi - Sets the vertices (imaginary part)
94: Notes:
95: In the case of complex scalars, only argument vr is used, containing
96: the complex vertices; the list of vertices can be provided in the
97: command line with a comma-separated list of complex values
98: [+/-][realnumber][+/-]realnumberi with no spaces.
100: When PETSc is built with real scalars, the real and imaginary parts of
101: the vertices must be provided in two separate arrays (or two lists in
102: the command line). In this case, the region must be symmetric with
103: respect to the real axis.
105: Level: advanced
107: .seealso: RGPolygonGetVertices()
108: @*/
109: PetscErrorCode RGPolygonSetVertices(RG rg,PetscInt n,PetscScalar vr[],PetscScalar vi[])
110: {
114: #if !defined(PETSC_USE_COMPLEX)
116: #endif
117: PetscTryMethod(rg,"RGPolygonSetVertices_C",(RG,PetscInt,PetscScalar*,PetscScalar*),(rg,n,vr,vi));
118: return 0;
119: }
121: static PetscErrorCode RGPolygonGetVertices_Polygon(RG rg,PetscInt *n,PetscScalar **vr,PetscScalar **vi)
122: {
123: RG_POLYGON *ctx = (RG_POLYGON*)rg->data;
124: PetscInt i;
126: if (n) *n = ctx->n;
127: if (vr) {
128: if (!ctx->n) *vr = NULL;
129: else {
130: PetscMalloc1(ctx->n,vr);
131: for (i=0;i<ctx->n;i++) (*vr)[i] = ctx->vr[i];
132: }
133: }
134: #if !defined(PETSC_USE_COMPLEX)
135: if (vi) {
136: if (!ctx->n) *vi = NULL;
137: else {
138: PetscMalloc1(ctx->n,vi);
139: for (i=0;i<ctx->n;i++) (*vi)[i] = ctx->vi[i];
140: }
141: }
142: #endif
143: return 0;
144: }
146: /*@C
147: RGPolygonGetVertices - Gets the vertices that define the polygon region.
149: Not Collective
151: Input Parameter:
152: . rg - the region context
154: Output Parameters:
155: + n - number of vertices
156: . vr - array of vertices
157: - vi - array of vertices (imaginary part)
159: Notes:
160: The values passed by user with RGPolygonSetVertices() are returned (or null
161: pointers otherwise).
162: The returned arrays should be freed by the user when no longer needed.
164: Level: advanced
166: .seealso: RGPolygonSetVertices()
167: @*/
168: PetscErrorCode RGPolygonGetVertices(RG rg,PetscInt *n,PetscScalar **vr,PetscScalar **vi)
169: {
171: PetscUseMethod(rg,"RGPolygonGetVertices_C",(RG,PetscInt*,PetscScalar**,PetscScalar**),(rg,n,vr,vi));
172: return 0;
173: }
175: PetscErrorCode RGView_Polygon(RG rg,PetscViewer viewer)
176: {
177: RG_POLYGON *ctx = (RG_POLYGON*)rg->data;
178: PetscBool isdraw,isascii;
179: int winw,winh;
180: PetscDraw draw;
181: PetscDrawAxis axis;
182: PetscReal a,b,c,d,ab,cd,lx,ly,w,x0,y0,x1,y1,scale=1.2;
183: PetscInt i;
184: char str[50];
186: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
187: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
188: if (isascii) {
189: PetscViewerASCIIPrintf(viewer," vertices: ");
190: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
191: for (i=0;i<ctx->n;i++) {
192: #if defined(PETSC_USE_COMPLEX)
193: SlepcSNPrintfScalar(str,sizeof(str),ctx->vr[i],PETSC_FALSE);
194: #else
195: if (ctx->vi[i]!=0.0) PetscSNPrintf(str,sizeof(str),"%g%+gi",(double)ctx->vr[i],(double)ctx->vi[i]);
196: else PetscSNPrintf(str,sizeof(str),"%g",(double)ctx->vr[i]);
197: #endif
198: PetscViewerASCIIPrintf(viewer,"%s%s",str,(i<ctx->n-1)?", ":"");
199: }
200: PetscViewerASCIIPrintf(viewer,"\n");
201: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
202: } else if (isdraw) {
203: PetscViewerDrawGetDraw(viewer,0,&draw);
204: PetscDrawCheckResizedWindow(draw);
205: PetscDrawGetWindowSize(draw,&winw,&winh);
206: winw = PetscMax(winw,1); winh = PetscMax(winh,1);
207: PetscDrawClear(draw);
208: PetscDrawSetTitle(draw,"Polygonal region");
209: PetscDrawAxisCreate(draw,&axis);
210: RGComputeBoundingBox_Polygon(rg,&a,&b,&c,&d);
211: a *= rg->sfactor;
212: b *= rg->sfactor;
213: c *= rg->sfactor;
214: d *= rg->sfactor;
215: lx = b-a;
216: ly = d-c;
217: ab = (a+b)/2;
218: cd = (c+d)/2;
219: w = scale*PetscMax(lx/winw,ly/winh)/2;
220: PetscDrawAxisSetLimits(axis,ab-w*winw,ab+w*winw,cd-w*winh,cd+w*winh);
221: PetscDrawAxisDraw(axis);
222: PetscDrawAxisDestroy(&axis);
223: for (i=0;i<ctx->n;i++) {
224: #if defined(PETSC_USE_COMPLEX)
225: x0 = PetscRealPart(ctx->vr[i]); y0 = PetscImaginaryPart(ctx->vr[i]);
226: if (i<ctx->n-1) {
227: x1 = PetscRealPart(ctx->vr[i+1]); y1 = PetscImaginaryPart(ctx->vr[i+1]);
228: } else {
229: x1 = PetscRealPart(ctx->vr[0]); y1 = PetscImaginaryPart(ctx->vr[0]);
230: }
231: #else
232: x0 = ctx->vr[i]; y0 = ctx->vi[i];
233: if (i<ctx->n-1) {
234: x1 = ctx->vr[i+1]; y1 = ctx->vi[i+1];
235: } else {
236: x1 = ctx->vr[0]; y1 = ctx->vi[0];
237: }
238: #endif
239: PetscDrawLine(draw,x0*rg->sfactor,y0*rg->sfactor,x1*rg->sfactor,y1*rg->sfactor,PETSC_DRAW_MAGENTA);
240: }
241: PetscDrawFlush(draw);
242: PetscDrawSave(draw);
243: PetscDrawPause(draw);
244: }
245: return 0;
246: }
248: PetscErrorCode RGIsTrivial_Polygon(RG rg,PetscBool *trivial)
249: {
250: RG_POLYGON *ctx = (RG_POLYGON*)rg->data;
252: *trivial = PetscNot(ctx->n);
253: return 0;
254: }
256: PetscErrorCode RGComputeContour_Polygon(RG rg,PetscInt n,PetscScalar *ucr,PetscScalar *uci)
257: {
258: RG_POLYGON *ctx = (RG_POLYGON*)rg->data;
259: PetscReal length,h,d,rem=0.0;
260: PetscInt k=1,idx=ctx->n-1,i;
261: PetscBool ini=PETSC_FALSE;
262: PetscScalar incr,*cr=ucr,*ci=uci;
263: #if !defined(PETSC_USE_COMPLEX)
264: PetscScalar inci;
265: #endif
268: length = SlepcAbsEigenvalue(ctx->vr[0]-ctx->vr[ctx->n-1],ctx->vi[0]-ctx->vi[ctx->n-1]);
269: for (i=0;i<ctx->n-1;i++) length += SlepcAbsEigenvalue(ctx->vr[i]-ctx->vr[i+1],ctx->vi[i]-ctx->vi[i+1]);
270: h = length/n;
271: if (!ucr) PetscMalloc1(n,&cr);
272: if (!uci) PetscMalloc1(n,&ci);
273: cr[0] = ctx->vr[0];
274: #if !defined(PETSC_USE_COMPLEX)
275: ci[0] = ctx->vi[0];
276: #endif
277: incr = ctx->vr[ctx->n-1]-ctx->vr[0];
278: #if !defined(PETSC_USE_COMPLEX)
279: inci = ctx->vi[ctx->n-1]-ctx->vi[0];
280: #endif
281: d = SlepcAbsEigenvalue(incr,inci);
282: incr /= d;
283: #if !defined(PETSC_USE_COMPLEX)
284: inci /= d;
285: #endif
286: while (k<n) {
287: if (ini) {
288: incr = ctx->vr[idx]-ctx->vr[idx+1];
289: #if !defined(PETSC_USE_COMPLEX)
290: inci = ctx->vi[idx]-ctx->vi[idx+1];
291: #endif
292: d = SlepcAbsEigenvalue(incr,inci);
293: incr /= d;
294: #if !defined(PETSC_USE_COMPLEX)
295: inci /= d;
296: #endif
297: if (rem+d>h) {
298: cr[k] = ctx->vr[idx+1]+incr*(h-rem);
299: #if !defined(PETSC_USE_COMPLEX)
300: ci[k] = ctx->vi[idx+1]+inci*(h-rem);
301: #endif
302: k++;
303: ini = PETSC_FALSE;
304: } else {rem += d; idx--;}
305: } else {
306: #if !defined(PETSC_USE_COMPLEX)
307: rem = SlepcAbsEigenvalue(ctx->vr[idx]-cr[k-1],ctx->vi[idx]-ci[k-1]);
308: #else
309: rem = PetscAbsScalar(ctx->vr[idx]-cr[k-1]);
310: #endif
311: if (rem>h) {
312: cr[k] = cr[k-1]+incr*h;
313: #if !defined(PETSC_USE_COMPLEX)
314: ci[k] = ci[k-1]+inci*h;
315: #endif
316: k++;
317: } else {ini = PETSC_TRUE; idx--;}
318: }
319: }
320: if (!ucr) PetscFree(cr);
321: if (!uci) PetscFree(ci);
322: return 0;
323: }
325: PetscErrorCode RGComputeBoundingBox_Polygon(RG rg,PetscReal *a,PetscReal *b,PetscReal *c,PetscReal *d)
326: {
327: RG_POLYGON *ctx = (RG_POLYGON*)rg->data;
328: PetscInt i;
330: if (a) *a = PETSC_MAX_REAL;
331: if (b) *b = -PETSC_MAX_REAL;
332: if (c) *c = PETSC_MAX_REAL;
333: if (d) *d = -PETSC_MAX_REAL;
334: for (i=0;i<ctx->n;i++) {
335: #if defined(PETSC_USE_COMPLEX)
336: if (a) *a = PetscMin(*a,PetscRealPart(ctx->vr[i]));
337: if (b) *b = PetscMax(*b,PetscRealPart(ctx->vr[i]));
338: if (c) *c = PetscMin(*c,PetscImaginaryPart(ctx->vr[i]));
339: if (d) *d = PetscMax(*d,PetscImaginaryPart(ctx->vr[i]));
340: #else
341: if (a) *a = PetscMin(*a,ctx->vr[i]);
342: if (b) *b = PetscMax(*b,ctx->vr[i]);
343: if (c) *c = PetscMin(*c,ctx->vi[i]);
344: if (d) *d = PetscMax(*d,ctx->vi[i]);
345: #endif
346: }
347: return 0;
348: }
350: PetscErrorCode RGCheckInside_Polygon(RG rg,PetscReal px,PetscReal py,PetscInt *inout)
351: {
352: RG_POLYGON *ctx = (RG_POLYGON*)rg->data;
353: PetscReal val,x[VERTMAX],y[VERTMAX];
354: PetscBool mx,my,nx,ny;
355: PetscInt i,j;
357: for (i=0;i<ctx->n;i++) {
358: #if defined(PETSC_USE_COMPLEX)
359: x[i] = PetscRealPart(ctx->vr[i])-px;
360: y[i] = PetscImaginaryPart(ctx->vr[i])-py;
361: #else
362: x[i] = ctx->vr[i]-px;
363: y[i] = ctx->vi[i]-py;
364: #endif
365: }
366: *inout = -1;
367: for (i=0;i<ctx->n;i++) {
368: j = (i+1)%ctx->n;
369: mx = PetscNot(x[i]<0.0);
370: nx = PetscNot(x[j]<0.0);
371: my = PetscNot(y[i]<0.0);
372: ny = PetscNot(y[j]<0.0);
373: if (!((my||ny) && (mx||nx)) || (mx&&nx)) continue;
374: if (((my && ny && (mx||nx)) && (!(mx&&nx)))) {
375: *inout = -*inout;
376: continue;
377: }
378: val = (y[i]*x[j]-x[i]*y[j])/(x[j]-x[i]);
379: if (PetscAbs(val)<10*PETSC_MACHINE_EPSILON) {
380: *inout = 0;
381: return 0;
382: } else if (val>0.0) *inout = -*inout;
383: }
384: return 0;
385: }
387: PetscErrorCode RGSetFromOptions_Polygon(RG rg,PetscOptionItems *PetscOptionsObject)
388: {
389: PetscScalar array[VERTMAX];
390: PetscInt i,k;
391: PetscBool flg,flgi=PETSC_FALSE;
392: #if !defined(PETSC_USE_COMPLEX)
393: PetscScalar arrayi[VERTMAX];
394: PetscInt ki;
395: #else
396: PetscScalar *arrayi=NULL;
397: #endif
399: PetscOptionsHeadBegin(PetscOptionsObject,"RG Polygon Options");
401: k = VERTMAX;
402: for (i=0;i<k;i++) array[i] = 0;
403: PetscOptionsScalarArray("-rg_polygon_vertices","Vertices of polygon","RGPolygonSetVertices",array,&k,&flg);
404: #if !defined(PETSC_USE_COMPLEX)
405: ki = VERTMAX;
406: for (i=0;i<ki;i++) arrayi[i] = 0;
407: PetscOptionsScalarArray("-rg_polygon_verticesi","Vertices of polygon (imaginary part)","RGPolygonSetVertices",arrayi,&ki,&flgi);
409: #endif
410: if (flg || flgi) RGPolygonSetVertices(rg,k,array,arrayi);
412: PetscOptionsHeadEnd();
413: return 0;
414: }
416: PetscErrorCode RGDestroy_Polygon(RG rg)
417: {
418: RG_POLYGON *ctx = (RG_POLYGON*)rg->data;
420: if (ctx->n) {
421: PetscFree(ctx->vr);
422: #if !defined(PETSC_USE_COMPLEX)
423: PetscFree(ctx->vi);
424: #endif
425: }
426: PetscFree(rg->data);
427: PetscObjectComposeFunction((PetscObject)rg,"RGPolygonSetVertices_C",NULL);
428: PetscObjectComposeFunction((PetscObject)rg,"RGPolygonGetVertices_C",NULL);
429: return 0;
430: }
432: SLEPC_EXTERN PetscErrorCode RGCreate_Polygon(RG rg)
433: {
434: RG_POLYGON *polygon;
436: PetscNew(&polygon);
437: rg->data = (void*)polygon;
439: rg->ops->istrivial = RGIsTrivial_Polygon;
440: rg->ops->computecontour = RGComputeContour_Polygon;
441: rg->ops->computebbox = RGComputeBoundingBox_Polygon;
442: rg->ops->checkinside = RGCheckInside_Polygon;
443: rg->ops->setfromoptions = RGSetFromOptions_Polygon;
444: rg->ops->view = RGView_Polygon;
445: rg->ops->destroy = RGDestroy_Polygon;
446: PetscObjectComposeFunction((PetscObject)rg,"RGPolygonSetVertices_C",RGPolygonSetVertices_Polygon);
447: PetscObjectComposeFunction((PetscObject)rg,"RGPolygonGetVertices_C",RGPolygonGetVertices_Polygon);
448: return 0;
449: }