Actual source code: rgring.c
slepc-3.11.2 2019-07-30
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2019, 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: Ring region, similar to the ellipse but with a start and end angle,
12: together with the width
13: */
15: #include <slepc/private/rgimpl.h> /*I "slepcrg.h" I*/
16: #include <petscdraw.h>
18: typedef struct {
19: PetscScalar center; /* center of the ellipse */
20: PetscReal radius; /* radius of the ellipse */
21: PetscReal vscale; /* vertical scale of the ellipse */
22: PetscReal start_ang; /* start angle */
23: PetscReal end_ang; /* end angle */
24: PetscReal width; /* ring width */
25: } RG_RING;
27: static PetscErrorCode RGRingSetParameters_Ring(RG rg,PetscScalar center,PetscReal radius,PetscReal vscale,PetscReal start_ang,PetscReal end_ang,PetscReal width)
28: {
29: RG_RING *ctx = (RG_RING*)rg->data;
32: ctx->center = center;
33: if (radius == PETSC_DEFAULT) {
34: ctx->radius = 1.0;
35: } else {
36: if (radius<=0.0) SETERRQ(PetscObjectComm((PetscObject)rg),PETSC_ERR_ARG_OUTOFRANGE,"The radius argument must be > 0.0");
37: ctx->radius = radius;
38: }
39: if (vscale<=0.0) SETERRQ(PetscObjectComm((PetscObject)rg),PETSC_ERR_ARG_OUTOFRANGE,"The vscale argument must be > 0.0");
40: ctx->vscale = vscale;
41: if (start_ang == PETSC_DEFAULT) {
42: ctx->start_ang = 0.0;
43: } else {
44: if (start_ang<0.0) SETERRQ(PetscObjectComm((PetscObject)rg),PETSC_ERR_ARG_OUTOFRANGE,"The right-hand side angle argument must be >= 0.0");
45: if (start_ang>1.0) SETERRQ(PetscObjectComm((PetscObject)rg),PETSC_ERR_ARG_OUTOFRANGE,"The right-hand side angle argument must be <= 1.0");
46: ctx->start_ang = start_ang;
47: }
48: if (end_ang == PETSC_DEFAULT) {
49: ctx->end_ang = 1.0;
50: } else {
51: if (end_ang<0.0) SETERRQ(PetscObjectComm((PetscObject)rg),PETSC_ERR_ARG_OUTOFRANGE,"The left-hand side angle argument must be >= 0.0");
52: if (end_ang>1.0) SETERRQ(PetscObjectComm((PetscObject)rg),PETSC_ERR_ARG_OUTOFRANGE,"The left-hand side angle argument must be <= 1.0");
53: ctx->end_ang = end_ang;
54: }
55: #if !defined(PETSC_USE_COMPLEX)
56: if (ctx->start_ang+ctx->end_ang!=1.0) SETERRQ(PetscObjectComm((PetscObject)rg),PETSC_ERR_ARG_WRONG,"In real scalars the region must be symmetric wrt real axis");
57: #endif
58: if (width == PETSC_DEFAULT) {
59: ctx->width = 0.1;
60: } else {
61: if (width<=0.0) SETERRQ(PetscObjectComm((PetscObject)rg),PETSC_ERR_ARG_OUTOFRANGE,"The width argument must be > 0.0");
62: ctx->width = width;
63: }
64: return(0);
65: }
67: /*@
68: RGRingSetParameters - Sets the parameters defining the ring region.
70: Logically Collective on RG
72: Input Parameters:
73: + rg - the region context
74: . center - center of the ellipse
75: . radius - radius of the ellipse
76: . vscale - vertical scale of the ellipse
77: . start_ang - the right-hand side angle
78: . end_ang - the left-hand side angle
79: - width - width of the ring
81: Options Database Keys:
82: + -rg_ring_center - Sets the center
83: . -rg_ring_radius - Sets the radius
84: . -rg_ring_vscale - Sets the vertical scale
85: . -rg_ring_startangle - Sets the right-hand side angle
86: . -rg_ring_endangle - Sets the left-hand side angle
87: - -rg_ring_width - Sets the width of the ring
89: Notes:
90: The values of center, radius and vscale have the same meaning as in the
91: ellipse region. The startangle and endangle define the span of the ring
92: (by default it is the whole ring), while the width is the separation
93: between the two concentric ellipses (above and below the radius by
94: width/2).
96: The start and end angles are expressed as a fraction of the circumference:
97: the allowed range is [0..1], with 0 corresponding to 0 radians, 0.25 to
98: pi/2 radians, and so on. It is allowed to have startangle>endangle, in
99: which case the ring region crosses over the zero angle.
101: In the case of complex scalars, a complex center can be provided in the
102: command line with [+/-][realnumber][+/-]realnumberi with no spaces, e.g.
103: -rg_ring_center 1.0+2.0i
105: When PETSc is built with real scalars, the center is restricted to a real value,
106: and the start and end angles must be such that the region is symmetric with
107: respect to the real axis.
109: Level: advanced
111: .seealso: RGRingGetParameters()
112: @*/
113: PetscErrorCode RGRingSetParameters(RG rg,PetscScalar center,PetscReal radius,PetscReal vscale,PetscReal start_ang,PetscReal end_ang,PetscReal width)
114: {
125: PetscTryMethod(rg,"RGRingSetParameters_C",(RG,PetscScalar,PetscReal,PetscReal,PetscReal,PetscReal,PetscReal),(rg,center,radius,vscale,start_ang,end_ang,width));
126: return(0);
127: }
129: static PetscErrorCode RGRingGetParameters_Ring(RG rg,PetscScalar *center,PetscReal *radius,PetscReal *vscale,PetscReal *start_ang,PetscReal *end_ang,PetscReal *width)
130: {
131: RG_RING *ctx = (RG_RING*)rg->data;
134: if (center) *center = ctx->center;
135: if (radius) *radius = ctx->radius;
136: if (vscale) *vscale = ctx->vscale;
137: if (start_ang) *start_ang = ctx->start_ang;
138: if (end_ang) *end_ang = ctx->end_ang;
139: if (width) *width = ctx->width;
140: return(0);
141: }
143: /*@
144: RGRingGetParameters - Gets the parameters that define the ring region.
146: Not Collective
148: Input Parameter:
149: . rg - the region context
151: Output Parameters:
152: + center - center of the region
153: . radius - radius of the region
154: . vscale - vertical scale of the region
155: . start_ang - the right-hand side angle
156: . end_ang - the left-hand side angle
157: - width - width of the ring
159: Level: advanced
161: .seealso: RGRingSetParameters()
162: @*/
163: PetscErrorCode RGRingGetParameters(RG rg,PetscScalar *center,PetscReal *radius,PetscReal *vscale,PetscReal *start_ang,PetscReal *end_ang,PetscReal *width)
164: {
169: PetscUseMethod(rg,"RGRingGetParameters_C",(RG,PetscScalar*,PetscReal*,PetscReal*,PetscReal*,PetscReal*,PetscReal*),(rg,center,radius,vscale,start_ang,end_ang,width));
170: return(0);
171: }
173: PetscErrorCode RGView_Ring(RG rg,PetscViewer viewer)
174: {
176: RG_RING *ctx = (RG_RING*)rg->data;
177: int winw,winh;
178: PetscBool isdraw,isascii;
179: PetscDraw draw;
180: PetscDrawAxis axis;
181: PetscReal cx,cy,radius,width,ab,cd,lx,ly,w,end_ang,x1,y1,x2,y2,r,theta,scale=1.2;
182: char str[50];
185: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
186: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
187: if (isascii) {
188: SlepcSNPrintfScalar(str,50,ctx->center,PETSC_FALSE);
189: PetscViewerASCIIPrintf(viewer," center: %s, radius: %g, vscale: %g, start angle: %g, end angle: %g, ring width: %g\n",str,RGShowReal(ctx->radius),RGShowReal(ctx->vscale),(double)ctx->start_ang,(double)ctx->end_ang,(double)ctx->width);
190: } else if (isdraw) {
191: PetscViewerDrawGetDraw(viewer,0,&draw);
192: PetscDrawCheckResizedWindow(draw);
193: PetscDrawGetWindowSize(draw,&winw,&winh);
194: winw = PetscMax(winw,1); winh = PetscMax(winh,1);
195: PetscDrawClear(draw);
196: PetscDrawSetTitle(draw,"Ring region");
197: PetscDrawAxisCreate(draw,&axis);
198: cx = PetscRealPart(ctx->center)*rg->sfactor;
199: cy = PetscImaginaryPart(ctx->center)*rg->sfactor;
200: radius = ctx->radius*rg->sfactor;
201: width = ctx->width*rg->sfactor;
202: lx = 2*(radius+width);
203: ly = 2*(radius+width)*ctx->vscale;
204: ab = cx;
205: cd = cy;
206: w = scale*PetscMax(lx/winw,ly/winh)/2;
207: PetscDrawAxisSetLimits(axis,ab-w*winw,ab+w*winw,cd-w*winh,cd+w*winh);
208: PetscDrawAxisDraw(axis);
209: PetscDrawAxisDestroy(&axis);
210: /* draw outer ellipse */
211: PetscDrawEllipse(draw,cx,cy,2*(radius+width),2*(radius+width)*ctx->vscale,PETSC_DRAW_ORANGE);
212: /* remove inner part */
213: PetscDrawEllipse(draw,cx,cy,2*(radius-width),2*(radius-width)*ctx->vscale,PETSC_DRAW_WHITE);
214: if (ctx->start_ang!=ctx->end_ang) {
215: /* remove section from end_ang to start_ang */
216: end_ang = (ctx->start_ang<ctx->end_ang)? ctx->end_ang-1: ctx->end_ang;
217: theta = end_ang;
218: r = scale*(radius+width);
219: if (ctx->vscale>1) r *= ctx->vscale;
220: x1 = PetscMin(PetscMax(ab+r*PetscCosReal(2.0*PETSC_PI*theta),ab-w*winw),ab+w*winw);
221: y1 = PetscMin(PetscMax(cd+r*PetscSinReal(2.0*PETSC_PI*theta),cd-w*winh),cd+w*winh);
222: do {
223: theta = PetscMin(PetscFloorReal(8*theta+1)/8,ctx->start_ang);
224: x2 = PetscMin(PetscMax(ab+r*PetscCosReal(2.0*PETSC_PI*theta),ab-w*winw),ab+w*winw);
225: y2 = PetscMin(PetscMax(cd+r*PetscSinReal(2.0*PETSC_PI*theta),cd-w*winh),cd+w*winh);
226: PetscDrawTriangle(draw,cx,cy,x1,y1,x2,y2,PETSC_DRAW_WHITE,PETSC_DRAW_WHITE,PETSC_DRAW_WHITE);
227: x1 = x2; y1 = y2;
228: } while (theta<ctx->start_ang);
229: }
230: PetscDrawFlush(draw);
231: PetscDrawSave(draw);
232: PetscDrawPause(draw);
233: }
234: return(0);
235: }
237: PetscErrorCode RGIsTrivial_Ring(RG rg,PetscBool *trivial)
238: {
239: RG_RING *ctx = (RG_RING*)rg->data;
242: if (rg->complement) *trivial = PetscNot(ctx->radius);
243: else *trivial = PetscNot(ctx->radius<PETSC_MAX_REAL);
244: return(0);
245: }
247: PetscErrorCode RGComputeContour_Ring(RG rg,PetscInt n,PetscScalar *cr,PetscScalar *ci)
248: {
249: RG_RING *ctx = (RG_RING*)rg->data;
250: PetscReal theta,start_ang;
251: PetscInt i,n2=n/2;
254: start_ang = (ctx->start_ang>ctx->end_ang)? ctx->start_ang-1: ctx->start_ang;
255: for (i=0;i<n;i++) {
256: if (i < n2) {
257: theta = ((ctx->end_ang-start_ang)*i/n2 + start_ang)*2.0*PETSC_PI;
258: #if defined(PETSC_USE_COMPLEX)
259: cr[i] = ctx->center + (ctx->radius+ctx->width/2.0)*(PetscCosReal(theta)+ctx->vscale*PetscSinReal(theta)*PETSC_i);
260: #else
261: cr[i] = ctx->center + (ctx->radius+ctx->width/2.0)*PetscCosReal(theta);
262: ci[i] = (ctx->radius+ctx->width/2.0)*ctx->vscale*PetscSinReal(theta);
263: #endif
264: } else {
265: theta = ((ctx->end_ang-start_ang)*(n-i)/n2 + start_ang)*2.0*PETSC_PI;
266: #if defined(PETSC_USE_COMPLEX)
267: cr[i] = ctx->center + (ctx->radius-ctx->width/2.0)*(PetscCosReal(theta)+ctx->vscale*PetscSinReal(theta)*PETSC_i);
268: #else
269: cr[i] = ctx->center + (ctx->radius-ctx->width/2.0)*PetscCosReal(theta);
270: ci[i] = (ctx->radius-ctx->width/2.0)*ctx->vscale*PetscSinReal(theta);
271: #endif
272: }
273: }
274: return(0);
275: }
277: PetscErrorCode RGComputeBoundingBox_Ring(RG rg,PetscReal *a,PetscReal *b,PetscReal *c,PetscReal *d)
278: {
279: RG_RING *ctx = (RG_RING*)rg->data;
282: /* current implementation does not return a tight bounding box */
283: if (a) *a = PetscRealPart(ctx->center) - (ctx->radius+ctx->width/2.0);
284: if (b) *b = PetscRealPart(ctx->center) + (ctx->radius+ctx->width/2.0);
285: if (c) *c = PetscImaginaryPart(ctx->center) - (ctx->radius+ctx->width/2.0)*ctx->vscale;
286: if (d) *d = PetscImaginaryPart(ctx->center) + (ctx->radius+ctx->width/2.0)*ctx->vscale;
287: return(0);
288: }
290: PetscErrorCode RGCheckInside_Ring(RG rg,PetscReal px,PetscReal py,PetscInt *inside)
291: {
292: RG_RING *ctx = (RG_RING*)rg->data;
293: PetscReal dx,dy,r;
296: /* outer ellipse */
297: #if defined(PETSC_USE_COMPLEX)
298: dx = (px-PetscRealPart(ctx->center))/(ctx->radius+ctx->width/2.0);
299: dy = (py-PetscImaginaryPart(ctx->center))/(ctx->radius+ctx->width/2.0);
300: #else
301: dx = (px-ctx->center)/(ctx->radius+ctx->width/2.0);
302: dy = py/(ctx->radius+ctx->width/2.0);
303: #endif
304: r = 1.0-dx*dx-(dy*dy)/(ctx->vscale*ctx->vscale);
305: *inside = PetscSign(r);
306: /* inner ellipse */
307: #if defined(PETSC_USE_COMPLEX)
308: dx = (px-PetscRealPart(ctx->center))/(ctx->radius-ctx->width/2.0);
309: dy = (py-PetscImaginaryPart(ctx->center))/(ctx->radius-ctx->width/2.0);
310: #else
311: dx = (px-ctx->center)/(ctx->radius-ctx->width/2.0);
312: dy = py/(ctx->radius-ctx->width/2.0);
313: #endif
314: r = -1.0+dx*dx+(dy*dy)/(ctx->vscale*ctx->vscale);
315: *inside *= PetscSign(r);
316: if (*inside == 1) { /* check angles */
317: #if defined(PETSC_USE_COMPLEX)
318: dx = (px-PetscRealPart(ctx->center));
319: dy = (py-PetscImaginaryPart(ctx->center));
320: #else
321: dx = px-ctx->center;
322: dy = py;
323: #endif
324: if (dx == 0) {
325: if (dy == 0) r = -1;
326: else if (dy > 0) r = 0.25;
327: else r = 0.75;
328: } else if (dx > 0) {
329: r = PetscAtanReal((dy/ctx->vscale)/dx);
330: if (dy >= 0) r /= 2*PETSC_PI;
331: else r = r/(2*PETSC_PI)+1;
332: } else r = PetscAtanReal((dy/ctx->vscale)/dx)/(2*PETSC_PI)+0.5;
333: if (ctx->start_ang>ctx->end_ang) {
334: if (r>ctx->end_ang && r<ctx->start_ang) *inside = -1;
335: } else {
336: if (r<ctx->start_ang || r>ctx->end_ang) *inside = -1;
337: }
338: }
339: return(0);
340: }
342: PetscErrorCode RGSetFromOptions_Ring(PetscOptionItems *PetscOptionsObject,RG rg)
343: {
345: PetscScalar s;
346: PetscReal r1,r2,r3,r4,r5;
347: PetscBool flg1,flg2,flg3,flg4,flg5,flg6;
350: PetscOptionsHead(PetscOptionsObject,"RG Ring Options");
352: RGRingGetParameters(rg,&s,&r1,&r2,&r3,&r4,&r5);
353: PetscOptionsScalar("-rg_ring_center","Center of ellipse","RGRingSetParameters",s,&s,&flg1);
354: PetscOptionsReal("-rg_ring_radius","Radius of ellipse","RGRingSetParameters",r1,&r1,&flg2);
355: PetscOptionsReal("-rg_ring_vscale","Vertical scale of ellipse","RGRingSetParameters",r2,&r2,&flg3);
356: PetscOptionsReal("-rg_ring_startangle","Right-hand side angle","RGRingSetParameters",r3,&r3,&flg4);
357: PetscOptionsReal("-rg_ring_endangle","Left-hand side angle","RGRingSetParameters",r4,&r4,&flg5);
358: PetscOptionsReal("-rg_ring_width","Width of ring","RGRingSetParameters",r5,&r5,&flg6);
359: if (flg1 || flg2 || flg3 || flg4 || flg5 || flg6) { RGRingSetParameters(rg,s,r1,r2,r3,r4,r5); }
361: PetscOptionsTail();
362: return(0);
363: }
365: PetscErrorCode RGDestroy_Ring(RG rg)
366: {
370: PetscFree(rg->data);
371: PetscObjectComposeFunction((PetscObject)rg,"RGRingSetParameters_C",NULL);
372: PetscObjectComposeFunction((PetscObject)rg,"RGRingGetParameters_C",NULL);
373: return(0);
374: }
376: SLEPC_EXTERN PetscErrorCode RGCreate_Ring(RG rg)
377: {
378: RG_RING *ring;
382: PetscNewLog(rg,&ring);
383: ring->center = 0.0;
384: ring->radius = PETSC_MAX_REAL;
385: ring->vscale = 1.0;
386: ring->start_ang = 0.0;
387: ring->end_ang = 1.0;
388: ring->width = 0.1;
389: rg->data = (void*)ring;
391: rg->ops->istrivial = RGIsTrivial_Ring;
392: rg->ops->computecontour = RGComputeContour_Ring;
393: rg->ops->computebbox = RGComputeBoundingBox_Ring;
394: rg->ops->checkinside = RGCheckInside_Ring;
395: rg->ops->setfromoptions = RGSetFromOptions_Ring;
396: rg->ops->view = RGView_Ring;
397: rg->ops->destroy = RGDestroy_Ring;
398: PetscObjectComposeFunction((PetscObject)rg,"RGRingSetParameters_C",RGRingSetParameters_Ring);
399: PetscObjectComposeFunction((PetscObject)rg,"RGRingGetParameters_C",RGRingGetParameters_Ring);
400: return(0);
401: }