Actual source code: rgpolygon.c
slepc-3.21.0 2024-03-30
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: }