Actual source code: rgring.c
slepc-main 2024-11-09
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: 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>
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;
31: PetscFunctionBegin;
32: ctx->center = center;
33: if (radius == (PetscReal)PETSC_DETERMINE) {
34: ctx->radius = 1.0;
35: } else if (radius != (PetscReal)PETSC_CURRENT) {
36: PetscCheck(radius>0.0,PetscObjectComm((PetscObject)rg),PETSC_ERR_ARG_OUTOFRANGE,"The radius argument must be > 0.0");
37: ctx->radius = radius;
38: }
39: if (vscale == (PetscReal)PETSC_DETERMINE) {
40: ctx->vscale = 1.0;
41: } else if (vscale != (PetscReal)PETSC_CURRENT) {
42: PetscCheck(vscale>0.0,PetscObjectComm((PetscObject)rg),PETSC_ERR_ARG_OUTOFRANGE,"The vscale argument must be > 0.0");
43: ctx->vscale = vscale;
44: }
45: if (start_ang == (PetscReal)PETSC_DETERMINE) {
46: ctx->start_ang = 0.0;
47: } else if (start_ang != (PetscReal)PETSC_CURRENT) {
48: PetscCheck(start_ang>=0.0,PetscObjectComm((PetscObject)rg),PETSC_ERR_ARG_OUTOFRANGE,"The right-hand side angle argument must be >= 0.0");
49: PetscCheck(start_ang<=1.0,PetscObjectComm((PetscObject)rg),PETSC_ERR_ARG_OUTOFRANGE,"The right-hand side angle argument must be <= 1.0");
50: ctx->start_ang = start_ang;
51: }
52: if (end_ang == (PetscReal)PETSC_DETERMINE) {
53: ctx->end_ang = 1.0;
54: } else if (end_ang != (PetscReal)PETSC_CURRENT) {
55: PetscCheck(end_ang>=0.0,PetscObjectComm((PetscObject)rg),PETSC_ERR_ARG_OUTOFRANGE,"The left-hand side angle argument must be >= 0.0");
56: PetscCheck(end_ang<=1.0,PetscObjectComm((PetscObject)rg),PETSC_ERR_ARG_OUTOFRANGE,"The left-hand side angle argument must be <= 1.0");
57: ctx->end_ang = end_ang;
58: }
59: #if !defined(PETSC_USE_COMPLEX)
60: PetscCheck(ctx->start_ang+ctx->end_ang==1.0,PetscObjectComm((PetscObject)rg),PETSC_ERR_ARG_WRONG,"In real scalars the region must be symmetric wrt real axis");
61: #endif
62: if (width == (PetscReal)PETSC_DETERMINE) {
63: ctx->width = 0.1;
64: } else if (width != (PetscReal)PETSC_CURRENT) {
65: PetscCheck(width>0.0,PetscObjectComm((PetscObject)rg),PETSC_ERR_ARG_OUTOFRANGE,"The width argument must be > 0.0");
66: ctx->width = width;
67: }
68: PetscFunctionReturn(PETSC_SUCCESS);
69: }
71: /*@
72: RGRingSetParameters - Sets the parameters defining the ring region.
74: Logically Collective
76: Input Parameters:
77: + rg - the region context
78: . center - center of the ellipse
79: . radius - radius of the ellipse
80: . vscale - vertical scale of the ellipse
81: . start_ang - the right-hand side angle
82: . end_ang - the left-hand side angle
83: - width - width of the ring
85: Options Database Keys:
86: + -rg_ring_center - Sets the center
87: . -rg_ring_radius - Sets the radius
88: . -rg_ring_vscale - Sets the vertical scale
89: . -rg_ring_startangle - Sets the right-hand side angle
90: . -rg_ring_endangle - Sets the left-hand side angle
91: - -rg_ring_width - Sets the width of the ring
93: Notes:
94: The values of center, radius and vscale have the same meaning as in the
95: ellipse region. The startangle and endangle define the span of the ring
96: (by default it is the whole ring), while the width is the separation
97: between the two concentric ellipses (above and below the radius by
98: width/2).
100: The start and end angles are expressed as a fraction of the circumference.
101: The allowed range is [0..1], with 0 corresponding to 0 radians, 0.25 to
102: pi/2 radians, and so on. It is allowed to have startangle>endangle, in
103: which case the ring region crosses over the zero angle.
105: In the case of complex scalars, a complex center can be provided in the
106: command line with [+/-][realnumber][+/-]realnumberi with no spaces, e.g.
107: -rg_ring_center 1.0+2.0i
109: When PETSc is built with real scalars, the center is restricted to a real value,
110: and the start and end angles must be such that the region is symmetric with
111: respect to the real axis.
113: For all arguments except center, you can use PETSC_CURRENT to keep the current
114: value, and PETSC_DETERMINE to set them to a default value (1 for radius, vscale,
115: end_ang, 0 for start_and, and 0.1 for width).
117: Level: advanced
119: .seealso: RGRingGetParameters()
120: @*/
121: PetscErrorCode RGRingSetParameters(RG rg,PetscScalar center,PetscReal radius,PetscReal vscale,PetscReal start_ang,PetscReal end_ang,PetscReal width)
122: {
123: PetscFunctionBegin;
131: PetscTryMethod(rg,"RGRingSetParameters_C",(RG,PetscScalar,PetscReal,PetscReal,PetscReal,PetscReal,PetscReal),(rg,center,radius,vscale,start_ang,end_ang,width));
132: PetscFunctionReturn(PETSC_SUCCESS);
133: }
135: static PetscErrorCode RGRingGetParameters_Ring(RG rg,PetscScalar *center,PetscReal *radius,PetscReal *vscale,PetscReal *start_ang,PetscReal *end_ang,PetscReal *width)
136: {
137: RG_RING *ctx = (RG_RING*)rg->data;
139: PetscFunctionBegin;
140: if (center) *center = ctx->center;
141: if (radius) *radius = ctx->radius;
142: if (vscale) *vscale = ctx->vscale;
143: if (start_ang) *start_ang = ctx->start_ang;
144: if (end_ang) *end_ang = ctx->end_ang;
145: if (width) *width = ctx->width;
146: PetscFunctionReturn(PETSC_SUCCESS);
147: }
149: /*@
150: RGRingGetParameters - Gets the parameters that define the ring region.
152: Not Collective
154: Input Parameter:
155: . rg - the region context
157: Output Parameters:
158: + center - center of the region
159: . radius - radius of the region
160: . vscale - vertical scale of the region
161: . start_ang - the right-hand side angle
162: . end_ang - the left-hand side angle
163: - width - width of the ring
165: Level: advanced
167: .seealso: RGRingSetParameters()
168: @*/
169: PetscErrorCode RGRingGetParameters(RG rg,PetscScalar *center,PetscReal *radius,PetscReal *vscale,PetscReal *start_ang,PetscReal *end_ang,PetscReal *width)
170: {
171: PetscFunctionBegin;
173: PetscUseMethod(rg,"RGRingGetParameters_C",(RG,PetscScalar*,PetscReal*,PetscReal*,PetscReal*,PetscReal*,PetscReal*),(rg,center,radius,vscale,start_ang,end_ang,width));
174: PetscFunctionReturn(PETSC_SUCCESS);
175: }
177: static PetscErrorCode RGView_Ring(RG rg,PetscViewer viewer)
178: {
179: RG_RING *ctx = (RG_RING*)rg->data;
180: int winw,winh;
181: PetscBool isdraw,isascii;
182: PetscDraw draw;
183: PetscDrawAxis axis;
184: PetscReal cx,cy,radius,width,ab,cd,lx,ly,w,end_ang,x1,y1,x2,y2,r,theta,scale=1.2;
185: char str[50];
187: PetscFunctionBegin;
188: PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw));
189: PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
190: if (isascii) {
191: PetscCall(SlepcSNPrintfScalar(str,sizeof(str),ctx->center,PETSC_FALSE));
192: PetscCall(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));
193: } else if (isdraw) {
194: PetscCall(PetscViewerDrawGetDraw(viewer,0,&draw));
195: PetscCall(PetscDrawCheckResizedWindow(draw));
196: PetscCall(PetscDrawGetWindowSize(draw,&winw,&winh));
197: winw = PetscMax(winw,1); winh = PetscMax(winh,1);
198: PetscCall(PetscDrawClear(draw));
199: PetscCall(PetscDrawSetTitle(draw,"Ring region"));
200: PetscCall(PetscDrawAxisCreate(draw,&axis));
201: cx = PetscRealPart(ctx->center)*rg->sfactor;
202: cy = PetscImaginaryPart(ctx->center)*rg->sfactor;
203: radius = ctx->radius*rg->sfactor;
204: width = ctx->width*rg->sfactor;
205: lx = 2*(radius+width);
206: ly = 2*(radius+width)*ctx->vscale;
207: ab = cx;
208: cd = cy;
209: w = scale*PetscMax(lx/winw,ly/winh)/2;
210: PetscCall(PetscDrawAxisSetLimits(axis,ab-w*winw,ab+w*winw,cd-w*winh,cd+w*winh));
211: PetscCall(PetscDrawAxisDraw(axis));
212: PetscCall(PetscDrawAxisDestroy(&axis));
213: /* draw outer ellipse */
214: PetscCall(PetscDrawEllipse(draw,cx,cy,2*(radius+width),2*(radius+width)*ctx->vscale,PETSC_DRAW_ORANGE));
215: /* remove inner part */
216: PetscCall(PetscDrawEllipse(draw,cx,cy,2*(radius-width),2*(radius-width)*ctx->vscale,PETSC_DRAW_WHITE));
217: if (ctx->start_ang!=ctx->end_ang) {
218: /* remove section from end_ang to start_ang */
219: end_ang = (ctx->start_ang<ctx->end_ang)? ctx->end_ang-1: ctx->end_ang;
220: theta = end_ang;
221: r = scale*(radius+width);
222: if (ctx->vscale>1) r *= ctx->vscale;
223: x1 = PetscMin(PetscMax(ab+r*PetscCosReal(2.0*PETSC_PI*theta),ab-w*winw),ab+w*winw);
224: y1 = PetscMin(PetscMax(cd+r*PetscSinReal(2.0*PETSC_PI*theta),cd-w*winh),cd+w*winh);
225: do {
226: theta = PetscMin(PetscFloorReal(8*theta+1)/8,ctx->start_ang);
227: x2 = PetscMin(PetscMax(ab+r*PetscCosReal(2.0*PETSC_PI*theta),ab-w*winw),ab+w*winw);
228: y2 = PetscMin(PetscMax(cd+r*PetscSinReal(2.0*PETSC_PI*theta),cd-w*winh),cd+w*winh);
229: PetscCall(PetscDrawTriangle(draw,cx,cy,x1,y1,x2,y2,PETSC_DRAW_WHITE,PETSC_DRAW_WHITE,PETSC_DRAW_WHITE));
230: x1 = x2; y1 = y2;
231: } while (theta<ctx->start_ang);
232: }
233: PetscCall(PetscDrawFlush(draw));
234: PetscCall(PetscDrawSave(draw));
235: PetscCall(PetscDrawPause(draw));
236: }
237: PetscFunctionReturn(PETSC_SUCCESS);
238: }
240: static PetscErrorCode RGIsTrivial_Ring(RG rg,PetscBool *trivial)
241: {
242: RG_RING *ctx = (RG_RING*)rg->data;
244: PetscFunctionBegin;
245: if (rg->complement) *trivial = PetscNot(ctx->radius);
246: else *trivial = PetscNot(ctx->radius<PETSC_MAX_REAL);
247: PetscFunctionReturn(PETSC_SUCCESS);
248: }
250: static PetscErrorCode RGComputeContour_Ring(RG rg,PetscInt n,PetscScalar *cr,PetscScalar *ci)
251: {
252: RG_RING *ctx = (RG_RING*)rg->data;
253: PetscReal theta,start_ang;
254: PetscInt i,n2=n/2;
256: PetscFunctionBegin;
257: start_ang = (ctx->start_ang>ctx->end_ang)? ctx->start_ang-1: ctx->start_ang;
258: for (i=0;i<n;i++) {
259: if (i < n2) {
260: theta = ((ctx->end_ang-start_ang)*i/n2 + start_ang)*2.0*PETSC_PI;
261: #if defined(PETSC_USE_COMPLEX)
262: cr[i] = ctx->center + (ctx->radius+ctx->width/2.0)*PetscCMPLX(PetscCosReal(theta),ctx->vscale*PetscSinReal(theta));
263: #else
264: if (cr) cr[i] = ctx->center + (ctx->radius+ctx->width/2.0)*PetscCosReal(theta);
265: if (ci) ci[i] = (ctx->radius+ctx->width/2.0)*ctx->vscale*PetscSinReal(theta);
266: #endif
267: } else {
268: theta = ((ctx->end_ang-start_ang)*(n-i)/n2 + start_ang)*2.0*PETSC_PI;
269: #if defined(PETSC_USE_COMPLEX)
270: cr[i] = ctx->center + (ctx->radius-ctx->width/2.0)*PetscCMPLX(PetscCosReal(theta),ctx->vscale*PetscSinReal(theta));
271: #else
272: if (cr) cr[i] = ctx->center + (ctx->radius-ctx->width/2.0)*PetscCosReal(theta);
273: if (ci) ci[i] = (ctx->radius-ctx->width/2.0)*ctx->vscale*PetscSinReal(theta);
274: #endif
275: }
276: }
277: PetscFunctionReturn(PETSC_SUCCESS);
278: }
280: static PetscErrorCode RGComputeBoundingBox_Ring(RG rg,PetscReal *a,PetscReal *b,PetscReal *c,PetscReal *d)
281: {
282: RG_RING *ctx = (RG_RING*)rg->data;
284: PetscFunctionBegin;
285: /* current implementation does not return a tight bounding box */
286: if (a) *a = PetscRealPart(ctx->center) - (ctx->radius+ctx->width/2.0);
287: if (b) *b = PetscRealPart(ctx->center) + (ctx->radius+ctx->width/2.0);
288: if (c) *c = PetscImaginaryPart(ctx->center) - (ctx->radius+ctx->width/2.0)*ctx->vscale;
289: if (d) *d = PetscImaginaryPart(ctx->center) + (ctx->radius+ctx->width/2.0)*ctx->vscale;
290: PetscFunctionReturn(PETSC_SUCCESS);
291: }
293: static PetscErrorCode RGComputeQuadrature_Ring(RG rg,RGQuadRule quad,PetscInt n,PetscScalar *z,PetscScalar *zn,PetscScalar *w)
294: {
295: RG_RING *ctx = (RG_RING*)rg->data;
296: PetscReal max_w=0.0;
297: PetscScalar tmp,tmp2;
298: PetscInt i,j;
300: PetscFunctionBegin;
301: if (quad == RG_QUADRULE_CHEBYSHEV) {
302: #if defined(PETSC_USE_COMPLEX)
303: PetscReal theta;
304: for (i=0;i<n;i++) {
305: theta = PETSC_PI*(i+0.5)/n;
306: zn[i] = PetscCosReal(theta);
307: w[i] = PetscCosReal((n-1)*theta)/n;
308: theta = (ctx->start_ang*2.0+(ctx->end_ang-ctx->start_ang)*(PetscRealPart(zn[i])+1.0))*PETSC_PI;
309: z[i] = rg->sfactor*(ctx->center + ctx->radius*PetscCMPLX(PetscCosReal(theta),ctx->vscale*PetscSinReal(theta)));
310: }
311: #else
312: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Integration points on a vertical line require complex arithmetic");
313: #endif
314: } else { /* RG_QUADRULE_TRAPEZOIDAL */
315: for (i=0;i<n;i++) {
316: zn[i] = (z[i]-rg->sfactor*ctx->center)/(rg->sfactor*ctx->radius);
317: tmp = 1.0; tmp2 = 1.0;
318: for (j=0;j<n;j++) {
319: tmp *= z[j];
320: if (i != j) tmp2 *= z[j]-z[i];
321: }
322: w[i] = tmp/tmp2;
323: max_w = PetscMax(PetscAbsScalar(w[i]),max_w);
324: }
325: for (i=0;i<n;i++) w[i] /= (PetscScalar)max_w;
326: }
327: PetscFunctionReturn(PETSC_SUCCESS);
328: }
330: static PetscErrorCode RGCheckInside_Ring(RG rg,PetscReal px,PetscReal py,PetscInt *inside)
331: {
332: RG_RING *ctx = (RG_RING*)rg->data;
333: PetscReal dx,dy,r;
335: PetscFunctionBegin;
336: /* outer ellipse */
337: #if defined(PETSC_USE_COMPLEX)
338: dx = (px-PetscRealPart(ctx->center))/(ctx->radius+ctx->width/2.0);
339: dy = (py-PetscImaginaryPart(ctx->center))/(ctx->radius+ctx->width/2.0);
340: #else
341: dx = (px-ctx->center)/(ctx->radius+ctx->width/2.0);
342: dy = py/(ctx->radius+ctx->width/2.0);
343: #endif
344: r = 1.0-dx*dx-(dy*dy)/(ctx->vscale*ctx->vscale);
345: *inside = PetscSign(r);
346: /* inner ellipse */
347: #if defined(PETSC_USE_COMPLEX)
348: dx = (px-PetscRealPart(ctx->center))/(ctx->radius-ctx->width/2.0);
349: dy = (py-PetscImaginaryPart(ctx->center))/(ctx->radius-ctx->width/2.0);
350: #else
351: dx = (px-ctx->center)/(ctx->radius-ctx->width/2.0);
352: dy = py/(ctx->radius-ctx->width/2.0);
353: #endif
354: r = -1.0+dx*dx+(dy*dy)/(ctx->vscale*ctx->vscale);
355: *inside *= PetscSign(r);
356: if (*inside == 1) { /* check angles */
357: #if defined(PETSC_USE_COMPLEX)
358: dx = (px-PetscRealPart(ctx->center));
359: dy = (py-PetscImaginaryPart(ctx->center));
360: #else
361: dx = px-ctx->center;
362: dy = py;
363: #endif
364: if (dx == 0) {
365: if (dy == 0) r = -1;
366: else if (dy > 0) r = 0.25;
367: else r = 0.75;
368: } else if (dx > 0) {
369: r = PetscAtanReal((dy/ctx->vscale)/dx);
370: if (dy >= 0) r /= 2*PETSC_PI;
371: else r = r/(2*PETSC_PI)+1;
372: } else r = PetscAtanReal((dy/ctx->vscale)/dx)/(2*PETSC_PI)+0.5;
373: if (ctx->start_ang>ctx->end_ang) {
374: if (r>ctx->end_ang && r<ctx->start_ang) *inside = -1;
375: } else {
376: if (r<ctx->start_ang || r>ctx->end_ang) *inside = -1;
377: }
378: }
379: PetscFunctionReturn(PETSC_SUCCESS);
380: }
382: static PetscErrorCode RGIsAxisymmetric_Ring(RG rg,PetscBool vertical,PetscBool *symm)
383: {
384: RG_RING *ctx = (RG_RING*)rg->data;
386: PetscFunctionBegin;
387: if (vertical) *symm = (PetscRealPart(ctx->center) == 0.0 && PetscAbs(ctx->start_ang+ctx->end_ang-PetscRealConstant(1.0)) == 0.5)? PETSC_TRUE: PETSC_FALSE;
388: else *symm = (PetscImaginaryPart(ctx->center) == 0.0 && ctx->start_ang+ctx->end_ang == 1.0)? PETSC_TRUE: PETSC_FALSE;
389: PetscFunctionReturn(PETSC_SUCCESS);
390: }
392: static PetscErrorCode RGSetFromOptions_Ring(RG rg,PetscOptionItems *PetscOptionsObject)
393: {
394: PetscScalar s;
395: PetscReal r1,r2,r3,r4,r5;
396: PetscBool flg1,flg2,flg3,flg4,flg5,flg6;
398: PetscFunctionBegin;
399: PetscOptionsHeadBegin(PetscOptionsObject,"RG Ring Options");
401: PetscCall(RGRingGetParameters(rg,&s,&r1,&r2,&r3,&r4,&r5));
402: PetscCall(PetscOptionsScalar("-rg_ring_center","Center of ellipse","RGRingSetParameters",s,&s,&flg1));
403: PetscCall(PetscOptionsReal("-rg_ring_radius","Radius of ellipse","RGRingSetParameters",r1,&r1,&flg2));
404: PetscCall(PetscOptionsReal("-rg_ring_vscale","Vertical scale of ellipse","RGRingSetParameters",r2,&r2,&flg3));
405: PetscCall(PetscOptionsReal("-rg_ring_startangle","Right-hand side angle","RGRingSetParameters",r3,&r3,&flg4));
406: PetscCall(PetscOptionsReal("-rg_ring_endangle","Left-hand side angle","RGRingSetParameters",r4,&r4,&flg5));
407: PetscCall(PetscOptionsReal("-rg_ring_width","Width of ring","RGRingSetParameters",r5,&r5,&flg6));
408: if (flg1 || flg2 || flg3 || flg4 || flg5 || flg6) PetscCall(RGRingSetParameters(rg,s,r1,r2,r3,r4,r5));
410: PetscOptionsHeadEnd();
411: PetscFunctionReturn(PETSC_SUCCESS);
412: }
414: static PetscErrorCode RGDestroy_Ring(RG rg)
415: {
416: PetscFunctionBegin;
417: PetscCall(PetscFree(rg->data));
418: PetscCall(PetscObjectComposeFunction((PetscObject)rg,"RGRingSetParameters_C",NULL));
419: PetscCall(PetscObjectComposeFunction((PetscObject)rg,"RGRingGetParameters_C",NULL));
420: PetscFunctionReturn(PETSC_SUCCESS);
421: }
423: SLEPC_EXTERN PetscErrorCode RGCreate_Ring(RG rg)
424: {
425: RG_RING *ring;
427: PetscFunctionBegin;
428: PetscCall(PetscNew(&ring));
429: ring->center = 0.0;
430: ring->radius = PETSC_MAX_REAL;
431: ring->vscale = 1.0;
432: ring->start_ang = 0.0;
433: ring->end_ang = 1.0;
434: ring->width = 0.1;
435: rg->data = (void*)ring;
437: rg->ops->istrivial = RGIsTrivial_Ring;
438: rg->ops->computecontour = RGComputeContour_Ring;
439: rg->ops->computebbox = RGComputeBoundingBox_Ring;
440: rg->ops->computequadrature = RGComputeQuadrature_Ring;
441: rg->ops->checkinside = RGCheckInside_Ring;
442: rg->ops->isaxisymmetric = RGIsAxisymmetric_Ring;
443: rg->ops->setfromoptions = RGSetFromOptions_Ring;
444: rg->ops->view = RGView_Ring;
445: rg->ops->destroy = RGDestroy_Ring;
446: PetscCall(PetscObjectComposeFunction((PetscObject)rg,"RGRingSetParameters_C",RGRingSetParameters_Ring));
447: PetscCall(PetscObjectComposeFunction((PetscObject)rg,"RGRingGetParameters_C",RGRingGetParameters_Ring));
448: PetscFunctionReturn(PETSC_SUCCESS);
449: }